unit CoMLFitter;
{$WARN SYMBOL_PLATFORM OFF}
{$B-,R-}
interface
uses
ComObj, ActiveX, Events, Classes, Fitter_TLB, StdVcl;
type
TLMFitter = class(TAutoObject, IConnectionPointContainer, IDMFitter)
private
FConnectionPoints: TConnectionPoints;
FConnectionPoint: TConnectionPoint;
FX,
FY,
FWeight,
FExpression,
FParameters,
FSigmas,
FOptions: OleVariant;
FXDimension,
FWeightType,
FParamCount,
FIterations,
FResultCode: integer;
FDeviation: double;
FStrExpr: array[0..4096] of char; FSig: array of boolean; procedure UnFixParameters(const DP: array of double);
procedure doProgress;
public
procedure Initialize; override;
protected
property ConnectionPoints: TConnectionPoints read FConnectionPoints
implements IConnectionPointContainer;
function Get_Deviation: Double; safecall;
function Get_Expression: OleVariant; safecall;
function Get_Iterations: Integer; safecall;
function Get_Options: OleVariant; safecall;
function Get_ParamCount: Integer; safecall;
function Get_Parameters: OleVariant; safecall;
function Get_ResultCode: Integer; safecall;
function Get_Sigmas: OleVariant; safecall;
function Get_Weight: OleVariant; safecall;
function Get_WeightType: Integer; safecall;
function Get_X: OleVariant; safecall;
function Get_Y: OleVariant; safecall;
function LinearFit: WordBool; safecall;
function LMFit: WordBool; safecall;
procedure Set_Expression(Value: OleVariant); safecall;
procedure Set_Iterations(Value: Integer); safecall;
procedure Set_Options(Value: OleVariant); safecall;
procedure Set_ParamCount(Value: Integer); safecall;
procedure Set_Parameters(Value: OleVariant); safecall;
procedure Set_Sigmas(Value: OleVariant); safecall;
procedure Set_Weight(Value: OleVariant); safecall;
procedure Set_WeightType(Value: Integer); safecall;
procedure Set_X(Value: OleVariant); safecall;
procedure Set_Y(Value: OleVariant); safecall;
end;
implementation
uses
Windows, ComServ, Variants, SysUtils, CoFitter, Parser, Wrapper{$endif}; const LM_DIFF_DELTA=1e-6; LM_STOP_THRESH=1E-17; LM_INIT_MU=1E-03; MLFITTER_MAXITER=9000;
MLFITTER_MAXPARAM=100;
type
DoubleArray=array[0..1] of double;
PDoubleArray=^DoubleArray;
TLMOptions=record
iterations: integer;
lambda, threshold1, threshold2, threshold3, delta, chisqr: double;
resultcode, n_func, n_jac, n_lin: integer;
end;
TLMNC=function(FitFunc, Reserved: pointer; DataVector, Parameters,
Sigmas: PDoubleArray; DataSize, ParamSize: integer;
var Options: TLMOptions): boolean; stdcall;
TLMJNC=function(FitFunc, DerFunc, Reserved: pointer; DataVector, Parameters,
Sigmas: PDoubleArray; DataSize, ParamSize: integer;
var Options: TLMOptions): boolean; stdcall;
TLMBC=function(FitFunc, Reserved: pointer; DataVector, Parameters,
Sigmas, LBound, UBound: PDoubleArray; DataSize, ParamSize: integer;
var Options: TLMOptions): boolean; stdcall;
TLMJBC=function(FitFunc, DerFunc, Reserved: pointer; DataVector, Parameters,
Sigmas, LBound, UBound: PDoubleArray; DataSize, ParamSize: integer;
var Options: TLMOptions): boolean; stdcall;
var
lh: THandle; LMNC: TLMNC; LMJNC: TLMJNC;
LMBC: TLMBC;
LMJBC: TLMJBC;
procedure ApplyWeight(vect: PDoubleArray; vsize: integer; FitObj: TLMFitter);
var
i: integer;
begin
if FitObj.FWeightType<>0 then
for i:=0 to vsize-1 do
vect^[i]:=vect^[i]*FitObj.FWeight[i];
end;
procedure ObjFitFunc(par, vect: PDoubleArray; psize, vsize: integer;
reserved: pointer); cdecl;
var
i: integer;
vY: variant;
begin
with TLMFitter(reserved) do
begin
UnFixParameters(par^);
FIterations:=FIterations+1;
vY:=FExpression.CalculateFunction(FX, FParameters);
end;
if VarIsArray(vY) and (VarArrayHighBound(vY,1)+1=vsize) then for i:=0 to vsize-1 do
vect^[i]:=vY[i];
ApplyWeight(vect, vsize, TLMFitter(reserved));
end;
procedure ObjDerFunc(par, vect: PDoubleArray; psize, vsize: integer;
reserved: pointer); cdecl;
var
i: integer;
vY: variant;
begin
with TLMFitter(reserved) do
begin
UnFixParameters(par^);
vY:=FExpression.CalculateDerivative(FX, FParameters, FSigmas);
end;
if VarIsArray(vY) and (VarArrayHighBound(vY,1)+1=vsize*psize) then for i:=0 to vsize*psize-1 do
vect^[i]:=vY[i];
ApplyWeight(vect, vsize*psize, TLMFitter(reserved));
end;
procedure StrFitFunc(par, vect: PDoubleArray; psize, vsize: integer;
reserved: pointer); cdecl;
var
i,j: integer;
PA: TRealArray;
X: TReal;
DA: array [0..MaxCols-1] of double;
begin
TLMFitter(reserved).UnFixParameters(par^);
TLMFitter(reserved).doProgress;
for i:=0 to TLMFitter(reserved).FParamCount-1 do
if i<MaxCols
then PA[i+1]:=TLMFitter(reserved).FParameters[i];
if vsize=VarArrayHighBound(TLMFitter(reserved).FX,1)+1 then for i:=0 to vsize-1 do begin
if TLMFitter(reserved).FXDimension=0
then X:=TLMFitter(reserved).FX[i] else
begin
for j:=0 to TLMFitter(reserved).FXDimension do if j<MaxCols
then DA[j]:=TLMFitter(reserved).FX[i,j];
SetXVector(DA);
X:=DA[0]; end;
vect^[i]:=NLSFParse(TLMFitter(reserved).FStrExpr, X, PA);
end;
ApplyWeight(vect, vsize, TLMFitter(reserved));
end;
procedure EvtFitFunc(par, vect: PDoubleArray; psize, vsize: integer;
reserved: pointer); cdecl;
var
i,j: integer;
X: variant;
R: double;
begin
with TLMFitter(reserved) do begin
UnFixParameters(par^);
doProgress;
end;
for i:=0 to vsize-1 do
begin
vect^[i]:=0; if TLMFitter(reserved).FXDimension=0
then X:=TLMFitter(reserved).FX[i] else
begin
X:=VarArrayCreate([0, TLMFitter(reserved).FXDimension], varVariant);
for j:=0 to TLMFitter(reserved).FXDimension do
X[j]:=TLMFitter(reserved).FX[i, j];
end;
with TLMFitter(reserved) do
if Assigned(FConnectionPoint) then
for j:=0 to FConnectionPoint.SinkList.Count-1 do
begin
R:=(IUnknown(FConnectionPoint.SinkList[j]) as
IDMFitterEvents).OnGetLMFunction(X, FParameters);
if R<>0 then vect^[i]:=R;
end;
end;
ApplyWeight(vect, vsize, TLMFitter(reserved));
end;
procedure EvtDerFunc(par, vect: PDoubleArray; psize, vsize: integer;
reserved: pointer); cdecl;
var
i,j,v_idx: integer;
X,Res,R: variant;
begin
with TLMFitter(reserved) do UnFixParameters(par^);
v_idx:=0; for i:=0 to vsize-1 do
begin
if TLMFitter(reserved).FXDimension=0
then X:=TLMFitter(reserved).FX[i] else
begin
X:=VarArrayCreate([0, TLMFitter(reserved).FXDimension], varVariant);
for j:=0 to TLMFitter(reserved).FXDimension do
X[j]:=TLMFitter(reserved).FX[i, j];
end;
R:=Unassigned;
Res:=Unassigned;
with TLMFitter(reserved) do
if Assigned(FConnectionPoint) then
for j:=0 to FConnectionPoint.SinkList.Count-1 do
begin
R:=(IUnknown(FConnectionPoint.SinkList[j]) as
IDMFitterEvents).OnGetLMDerivative(X, FParameters, FSigmas);
if not VarIsEmpty(R) then Res:=SetVariantArray(R); end;
Assert(VarIsArray(Res)); Assert(VarArrayHighBound(Res,1)=TLMFitter(reserved).FParamCount-1);
with TLMFitter(reserved) do
for j:=0 to FParamCount-1 do
if FSig[j] then
begin
Assert(v_idx<vsize*psize);
vect^[v_idx]:=Res[j];
Inc(v_idx);
end;
end; ApplyWeight(vect, vsize*psize, TLMFitter(reserved));
end;
procedure TLMFitter.UnFixParameters(const DP: array of double);
var
i, N: integer;
begin
N:=0;
Assert(Length(FSig)=FParamCount);
Assert(VarArrayHighBound(FParameters,1)=FParamCount-1);
for i:=0 to FParamCount-1 do
if FSig[i] then
begin
FParameters[i]:=DP[N];
Inc(N);
end;
end;
procedure TLMFitter.doProgress;
var
i: integer;
begin
FIterations:=FIterations+1;
if Assigned(FConnectionPoint) then
for i:=0 to FConnectionPoint.SinkList.Count-1 do
(IUnknown(FConnectionPoint.SinkList[I]) as
IDMFitterEvents).OnProgress;
end;
procedure TLMFitter.Initialize;
begin
inherited Initialize;
FConnectionPoints := TConnectionPoints.Create(Self);
if AutoFactory.EventTypeInfo <> nil then
FConnectionPoint := FConnectionPoints.CreateConnectionPoint(
AutoFactory.EventIID, ckMulti, EventConnect)
else FConnectionPoint := nil;
end;
function TLMFitter.Get_Deviation: Double;
begin
Result:=FDeviation;
end;
function TLMFitter.Get_Expression: OleVariant;
begin
Result:=FExpression;
end;
function TLMFitter.Get_Iterations: Integer;
begin
Result:=FIterations;
end;
function TLMFitter.Get_Options: OleVariant;
begin
Result:=FOptions;
end;
function TLMFitter.Get_ParamCount: Integer;
begin
Result:=FParamCount;
end;
function TLMFitter.Get_Parameters: OleVariant;
begin
Result:=FParameters;
end;
function TLMFitter.Get_ResultCode: Integer;
begin
Result:=FResultCode;
end;
function TLMFitter.Get_Sigmas: OleVariant;
begin
Result:=FSigmas;
end;
function TLMFitter.Get_Weight: OleVariant;
begin
Result:=FWeight;
end;
function TLMFitter.Get_WeightType: Integer;
begin
Result:=FWeightType;
end;
function TLMFitter.Get_X: OleVariant;
begin
Result:=FX;
end;
function TLMFitter.Get_Y: OleVariant;
begin
Result:=FY;
end;
function TLMFitter.LinearFit: WordBool;
begin
Result:=false; end;
function TLMFitter.LMFit: WordBool;
var
i, Fl, NumPoints: integer;
ffu: (ffBad, ffClass, ffClassDer, ffString, ffEvent, ffEventDer);
Dat, Par, Sig, Ubo, Lbo: array of double;
Opt: TLMOptions;
bc, bc1: boolean; Wt: double;
begin
Result:=false; FResultCode:=7;
if (lh=0) or (not Assigned(LMNC)) or (not Assigned(LMJNC))
or (not Assigned(LMBC)) or (not Assigned(LMJBC))
then Exit; FResultCode:=3;
if (FIterations<1) or (FIterations>MLFITTER_MAXITER)
then Exit;
if (FParamCount<1) or (FParamCount>MLFITTER_MAXPARAM)
then Exit; if not VarIsArray(FY)
then Exit;
NumPoints:=VarArrayHighBound(FY,1)+1;
if not (VarIsArray(FParameters) and VarIsArray(FSigmas))
then Exit;
if VarArrayHighBound(FParameters,1)<>VarArrayHighBound(FSigmas,1)
then Exit;
if VarArrayHighBound(FParameters,1)<>FParamCount-1
then Exit;
if FWeightType<>0 then
begin
if not VarIsArray(FWeight)
then Exit;
if VarArrayHighBound(FWeight,1)+1<>NumPoints
then Exit;
end;
ffu:=ffBad;
if VarIsType(FExpression, varDispatch) then ffu:=ffClass;
if ffu=ffClass then if FExpression.CanCalculateDerivative
then ffu:=ffClassDer; if VarIsStr(FExpression) then
begin
if FParamCount>MaxCols
then Exit; if (not VarIsArray(FX)) or (VarArrayHighBound(FX,1)+1<>NumPoints)
then Exit; StrPLCopy(FStrExpr, FExpression, SizeOf(FStrExpr)-1);
if (Pos('CX2', UpperCase(FStrExpr))>0) and (FXDimension=0)
then Exit; ffu:=ffString;
end;
if VarIsOrdinal(FExpression) then
begin
if (not VarIsArray(FX)) or (VarArrayHighBound(FX,1)+1<>NumPoints)
then Exit; if FExpression=-1
then ffu:=ffEventDer else ffu:=ffEvent; end;
if ffu=ffBad
then Exit; SetLength(Dat, NumPoints);
SetLength(Par, FParamCount);
SetLength(Sig, FParamCount);
SetLength(FSig, FParamCount);
for i:=0 to Length(Dat)-1 do
Dat[i]:=FY[i];
ApplyWeight(@Dat[0], NumPoints, Self);
bc:=(VarArrayDimCount(FSigmas)=2) and (VarArrayHighBound(FSigmas,2)=2);
bc1:=(VarArrayDimCount(FSigmas)=1) and VarIsStr(FSigmas[0]); if bc or bc1 then begin
SetLength(Ubo, FParamCount);
SetLength(Lbo, FParamCount);
end;
for i:=0 to FParamCount-1 do
begin
Par[i]:=FParameters[i];
if not (bc or bc1)
then Sig[i]:=FSigmas[i];
if bc then begin
Sig[i]:=FSigmas[i,0];
Lbo[i]:=FSigmas[i,1];
Ubo[i]:=FSigmas[i,2];
end;
if bc1 then with TStringList.Create do
try
Delimiter:=' '; if VarIsStr(FSigmas[i])
then DelimitedText:=FSigmas[i]
else Exit; if Count<>3
then Exit;
try
Lbo[i]:=StrToFloat(Strings[1]); except
Val(Strings[1], Lbo[i], Fl); end;
if Fl<>0 then Exit; try
Ubo[i]:=StrToFloat(Strings[2]);
except
Val(Strings[2], Ubo[i], Fl);
end;
if Fl<>0 then Exit;
try
Sig[i]:=StrToFloat(Strings[0]);
except
Val(Strings[0], Sig[i], Fl);
end;
if Fl<>0 then Exit;
finally
Free;
end;
FSig[i]:=Sig[i]=0;
end;
Fl:=0; for i:=0 to FParamCount-1 do
if FSig[i] then
begin
Par[Fl]:=Par[i];
if bc or bc1 then begin
Ubo[Fl]:=Ubo[i];
Lbo[Fl]:=Lbo[i];
end;
Inc(Fl);
end;
Opt.iterations:=FIterations;
FIterations:=0; if VarIsArray(FOptions) and (VarArrayHighBound(FOptions,1)=4) then
begin
Opt.lambda:=FOptions[0];
if Opt.lambda=0
then Opt.lambda:=LM_INIT_MU; Opt.threshold1:=FOptions[1];
Opt.threshold2:=FOptions[2];
Opt.threshold3:=FOptions[3];
Opt.delta:=FOptions[4];
if Opt.delta=0
then Opt.delta:=LM_DIFF_DELTA; end else
begin Opt.lambda:=LM_INIT_MU;
Opt.threshold1:=LM_STOP_THRESH;
Opt.threshold2:=LM_STOP_THRESH;
Opt.threshold3:=LM_STOP_THRESH;
Opt.delta:=LM_DIFF_DELTA;
end;
FResultCode:=1;
if bc or bc1 then
case ffu of
ffString: Result:=LMBC(@StrFitFunc, Self, @Dat[0], @Par[0], @Sig[0],
@Lbo[0], @Ubo[0], Length(Dat), Fl, Opt);
ffClass: Result:=LMBC(@ObjFitFunc, Self, @Dat[0], @Par[0], @Sig[0],
@Lbo[0], @Ubo[0], Length(Dat), Fl, Opt);
ffClassDer: Result:=LMJBC(@ObjFitFunc, @ObjDerFunc, Self, @Dat[0],
@Par[0], @Sig[0], @Lbo[0], @Ubo[0], Length(Dat), Fl, Opt);
ffEvent: Result:=LMBC(@EvtFitFunc, Self, @Dat[0], @Par[0], @Sig[0],
@Lbo[0], @Ubo[0], Length(Dat), Fl, Opt);
ffEventDer: Result:=LMJBC(@EvtFitFunc, @EvtDerFunc, Self, @Dat[0],
@Par[0], @Sig[0], @Lbo[0], @Ubo[0], Length(Dat), Fl, Opt);
else
Exit; end else
case ffu of
ffString: Result:=LMNC(@StrFitFunc, Self, @Dat[0], @Par[0], @Sig[0],
Length(Dat), Fl, Opt);
ffClass: Result:=LMNC(@ObjFitFunc, Self, @Dat[0], @Par[0], @Sig[0],
Length(Dat), Fl, Opt);
ffClassDer: Result:=LMJNC(@ObjFitFunc, @ObjDerFunc, Self,
@Dat[0], @Par[0], @Sig[0], Length(Dat), Fl, Opt);
ffEvent: Result:=LMNC(@EvtFitFunc, Self, @Dat[0], @Par[0], @Sig[0],
Length(Dat), Fl, Opt);
ffEventDer: Result:=LMJNC(@EvtFitFunc, @EvtDerFunc, Self,
@Dat[0], @Par[0], @Sig[0], Length(Dat), Fl, Opt);
else
Exit;
end;
FIterations:=Opt.iterations;
Fl:=0;
for i:=0 to FParamCount-1 do
if FSig[i] then
begin
if bc
then FSigmas[i,0]:=Sig[Fl]
else FSigmas[i]:=Sig[Fl];
Inc(Fl);
end;
case Opt.resultcode of 1: FResultCode:=4; 2: FResultCode:=4; 3: FResultCode:=6; 4: FResultCode:=2; 5: FResultCode:=0; 6: FResultCode:=5; 7: FResultCode:=2; end;
if FWeightType<>0 then
begin
Wt:=0;
for i:=0 to NumPoints-1 do
Wt:=Wt+FWeight[i];
Assert(Wt<>0);
FDeviation:=Opt.chisqr*Abs(NumPoints/Wt);
end else
FDeviation:=Opt.chisqr;
end;
procedure TLMFitter.Set_Expression(Value: OleVariant);
begin
FExpression:=Value;
end;
procedure TLMFitter.Set_Iterations(Value: Integer);
begin
FIterations:=Value;
end;
procedure TLMFitter.Set_Options(Value: OleVariant);
begin
FOptions:=SetVariantArray(Value);
end;
procedure TLMFitter.Set_ParamCount(Value: Integer);
begin
FParamCount:=Value;
end;
procedure TLMFitter.Set_Parameters(Value: OleVariant);
begin
FParameters:=SetVariantArray(Value);
end;
procedure TLMFitter.Set_Sigmas(Value: OleVariant);
begin
if VarArrayDimCount(Value)=2
then FSigmas:=Value
else FSigmas:=SetVariantArray(Value);
end;
procedure TLMFitter.Set_Weight(Value: OleVariant);
begin
FWeight:=SetVariantArray(Value);
end;
procedure TLMFitter.Set_WeightType(Value: Integer);
begin
FWeightType:=Value;
end;
procedure TLMFitter.Set_X(Value: OleVariant);
begin if VarArrayDimCount(Value)=2 then
begin
FX:=Value; FXDimension:=VarArrayHighBound(Value, 2);
end else
begin
FX:=SetVariantArray(Value);
FXDimension:=0;
end;
end;
procedure TLMFitter.Set_Y(Value: OleVariant);
begin
FY:=SetVariantArray(Value);
end;
var
FNBuf: array[0..MAX_PATH] of char;
initialization
TAutoObjectFactory.Create(ComServer, TLMFitter, Class_LMFitter,
ciMultiInstance, tmApartment);
GetModuleFilename(hInstance, FNBuf, SizeOf(FNBuf)-1);
StrPLCopy(FNBuf, ExtractFilePath(StrPas(FNBuf))+'lmhelper.dll',
SizeOf(FNBuf)-1);
lh:=LoadLibrary(FNBuf);
if lh<>0 then
begin
LMNC:=GetProcAddress(lh, 'LMNC'); LMJNC:=GetProcAddress(lh, 'LMJNC');
LMBC:=GetProcAddress(lh, 'LMBC');
LMJBC:=GetProcAddress(lh, 'LMJBC');
end;
finalization
if lh<>0
then FreeLibrary(lh);
end.