///////////////////////////////////////
//        Data Master 2003           //
//   Copyright (c) 1993-2003 RRR     //
///////////////////////////////////////

unit CoFitter;

{$WARN SYMBOL_PLATFORM OFF}
{$B-}

interface

uses
  ComObj, ActiveX, {AxCtrls}Events, Classes, Fitter_TLB, StdVcl;

type
  TDMFitter = class(TAutoObject, IConnectionPointContainer, IDMFitter)
  private
    { Private declarations }
    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; // string expression buffer
    function SetVariantArray(const Value: variant): variant;
  public
    { Public declarations }
    procedure Initialize; override;
  protected
    { Protected declarations }
    property ConnectionPoints: TConnectionPoints read FConnectionPoints
      implements IConnectionPointContainer;
    function Get_Expression: OleVariant; safecall;
    function LinearFit: WordBool; safecall;
    function LMFit: WordBool; safecall;
    procedure Set_Expression(Value: OleVariant); safecall;
    function Get_Deviation: Double; 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;
    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

// Notice: a lot of code is just copy & paste from CoMain.pas

uses ComServ, Windows, Variants, SysUtils, Parser, LMFit, Wrapper;

threadvar // Warning: global variables may be NOT thread-safe!
  Fitter: variant; // reference to EXTERNAL object used for LM fitting 
  DMFitter: TDMFitter; // reference to the current instance of DMFitter
  Pars, Sig, XVector: TRealArray;
  UseXVector: boolean; // used in LMFitFunc3; true if parsing CXn

function GetX(X: TReal): variant; // returns X - single value or array
var                               // Warning: uses global DMFitter!
  I: integer;
begin                             
  if DMFitter.FXDimension>0 then  // X is integer index
  begin
    Result:=VarArrayCreate([0, DMFitter.FXDimension], varVariant);
    Assert((Round(X)>=0) and (Round(X)<=(VarArrayHighBound(DMFitter.FX, 1))));
    for I:=0 to DMFitter.FXDimension do
    Result[I]:=DMFitter.FX[Round(X), I]; 
  end else 
    Result:=X; // "normal" real X
end;
  
// Following 2 callbacks support EXTERNAL object used for LM fitting

function LMFitFunc(X: TReal; A: FitArray; NParam: Integer): TReal; stdcall;
var
  vP: variant;
  I: integer;
begin
  vP:=VarArrayCreate([0, NParam-1], varVariant);
  for I:=0 to NParam-1 do
  vP[I]:=A[I+1];
  Result:=Fitter.CalculateFunction(GetX(X), vP);
end;

procedure LMFitDeriv(X: TReal; A, SigA: FitArray; NParam: Integer;
  var Deriv: FitArray); stdcall;
var
  vPar, vSig, vRes: variant;
  I: integer;
begin
  vPar:=VarArrayCreate([0, NParam-1], varVariant);
  vSig:=VarArrayCreate([0, NParam-1], varVariant);
  for I:=0 to NParam-1 do
  begin
    vPar[I]:=A[I+1];
    vSig[I]:=SigA[I+1];
  end;
  vRes:=Fitter.CalculateDerivative(GetX(X), vPar, vSig);
  for I:=0 to NParam-1 do
  Deriv[I+1]:=vRes[I]; // no any checks!!!
end;

// LM fitting progress and linear fitting with user-defined basis

procedure LMFitProgress; stdcall;
var
  I: integer;
begin
  if Assigned(DMFitter) then
  begin
    for I:=0 to VarArrayHighBound(DMFitter.FParameters,1) do
    begin // allows to create ANIMATED fitter
      DMFitter.FParameters[I]:=Pars[I+1];
      DMFitter.FSigmas[I]:=Sig[I+1];
    end;
    if Assigned(DMFitter.FConnectionPoint) then
    for I:=0 to DMFitter.FConnectionPoint.SinkList.Count-1 do
      (IUnknown(DMFitter.FConnectionPoint.SinkList[I]) as 
      IDMFitterEvents).OnProgress;
  end;
end;

function LinearFitBasis(X: double; Term: integer): double; stdcall;
var
  I: integer;
  R: double; 
begin
  Result:=0;
  if Assigned(DMFitter.FConnectionPoint) then
  for I:=0 to DMFitter.FConnectionPoint.SinkList.Count-1 do
  begin
    R:=(IUnknown(DMFitter.FConnectionPoint.SinkList[I]) as 
    IDMFitterEvents).OnGetLinearBasis(X, Term);
    if R<>0 then
    Result:=R;
  end;
end;

// This callback supports multidimensional X in expressions

function LMFitFunc3(X: TReal; A: FitArray; NParam: Integer): TReal; stdcall;
var
  I: integer;
begin
  Assert(UseXVector);
  Assert(DMFitter.FXDimension>0);
  Assert((Round(X)>=0) and (Round(X)<=(VarArrayHighBound(DMFitter.FX, 1))));
  for I:=0 to DMFitter.FXDimension do // initialize CXn array
  XVector[I+1]:=DMFitter.FX[Round(X), I]; // warning: X - integer INDEX!
  Result:=NLSFParse(DMFitter.FStrExpr, XVector[1], TRealArray(A));
end;

// Next 3 callbacks support event-handled LM fitting (using global DMFitter var)

function LMFitFunc2(X: TReal; A: FitArray; NParam: Integer): TReal; stdcall;
var
  VA: variant;
  I: integer;
  R: double; 
  {We use result buffer because under .Net event handler may be called several 
  times, and only one call is useful (others return zero)}
begin
  VA:=VarArrayCreate([0, NParam-1], varVariant);
  for I:=0 to NParam-1 do
  VA[I]:=A[I+1];
  Result:=0;
  if Assigned(DMFitter.FConnectionPoint) then
  for I:=0 to DMFitter.FConnectionPoint.SinkList.Count-1 do
  begin
    R:=(IUnknown(DMFitter.FConnectionPoint.SinkList[I]) as 
    IDMFitterEvents).OnGetLMFunction(GetX(X), VA);
    if R<>0 then // zero result may mean "unhandled" event
    Result:=R;
  end;
end;

procedure LMFitDeriv2(X: TReal; A, SigA: FitArray; NParam: Integer;
  var Deriv: FitArray); stdcall;
var
  VA, VSigA, Res, R: variant;
  I: integer;
begin
  VA:=VarArrayCreate([0, NParam-1], varVariant);
  VSigA:=VarArrayCreate([0, NParam-1], varVariant);
  for I:=0 to NParam-1 do
  begin
    VA[I]:=A[I+1];
    VSigA[I]:=SigA[I+1];
  end;
  R:=Unassigned;
  Res:=Unassigned;
  if Assigned(DMFitter.FConnectionPoint) then
  for I:=0 to DMFitter.FConnectionPoint.SinkList.Count-1 do
  begin
    R:=(IUnknown(DMFitter.FConnectionPoint.SinkList[I]) as 
    IDMFitterEvents).OnGetLMDerivative(GetX(X), VA, VSigA);
    if not VarIsEmpty(R) then // empty result may mean unhandled event
    Res:=DMFitter.SetVariantArray(R);
  end;
  Assert(VarIsArray(Res)); // some checks
  Assert(VarArrayHighBound(Res,1)=NParam-1);
  for I:=0 to VarArrayHighBound(Res,1) do
  Deriv[I+1]:=Res[I];
end;

procedure LMFitDeriv3(X: TReal; A, SigA: FitArray; NParam: Integer;
  var Deriv: FitArray); stdcall; // as in Wrapper.LMDerivative()
var
  I: integer;
  A1: FitArray;
  Y1, Y2, dP, NLSFDelta: TReal;
begin
  // calculate derivative Delta
  if VarIsArray(DMFitter.FOptions) and (VarArrayHighBound(DMFitter.FOptions,1)=2)
  then NLSFDelta:=DMFitter.FOptions[2]
  else NLSFDelta:=1e-10;
  if NLSFDelta<=0 then NLSFDelta:=1e-10; // value of FOptions not checked!
  // calculate derivatives
  for I:=1 to Nparam do
  if SigA[I]<0 then Deriv[I]:=0 else
  begin
    if Abs(A[I])<NLSFDelta
    then dP:=NLSFDelta
    else dP:=Abs(A[I])*NLSFDelta;
    A1:=A;
    A1[I]:=A1[I]-dP;
    if UseXVector 
    then Y1:=LMFitFunc3(X, A1, Nparam) 
    else Y1:=LMFitFunc2(X, A1, Nparam);
    A1:=A;
    A1[I]:=A1[I]+dP;
    if UseXVector 
    then Y2:=LMFitFunc3(X, A1, Nparam) 
    else Y2:=LMFitFunc2(X, A1, Nparam);
    Deriv[I]:=(Y2-Y1)/2/dP;
  end;
end;

{ TDMFitter }

procedure TDMFitter.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 TDMFitter.LinearFit: WordBool;
var
  DX, DY, Res: array of double;
  Va: double;
  I, NumPoints, FitType: integer;
begin
  Result:=false; // failure by default
  FResultCode:=3{frBadParams};
  if not (VarIsArray(FX) and VarIsArray(FY))
  then Exit; // must be arrays
  NumPoints:=VarArrayHighBound(FX, 1)+1;
  if (VarArrayHighBound(FY, 1)+1<>NumPoints)
  then Exit; // check array length
  if VarIsStr(FExpression) then
  begin
    StrPLCopy(FStrExpr, FExpression, SizeOf(FStrExpr)-1);
    FitType:=-2; // -2 means string Expression
  end else
  if VarIsOrdinal(FExpression) then
  begin
    FitType:=FExpression;
    if not ((FitType in [0..5]) or (FitType=-1)) // -1 means use EVENT!
    then Exit; // invalid basis type (see Wrapper.pas)
  end else Exit; // Expression must be integer
  FResultCode:=1{frError};
  SetLength(DX, NumPoints);
  SetLength(DY, NumPoints); // init arrays
  SetLength(Res, 100);
  for I:=0 to NumPoints-1 do  // copy data
  begin
    DX[I]:=FX[I];
    DY[I]:=FY[I];
  end;
  if FitType=-1 then
  try
    DMFitter:=self; // assign global reference used in callback
    FResultCode:=Wrapper.LinearFitEx(NumPoints,
      FParamCount or LinearFitExNumTermsMask, @LinearFitBasis,
      @DX[0], @DY[0], @Res[0], FDeviation, Va)
  finally
    DMFitter:=nil;
  end else
    if FitType=-2
    then FResultCode:=Wrapper.LinearFitEx(NumPoints, FParamCount, @FStrExpr,
           @DX[0], @DY[0], @Res[0], FDeviation, Va)
    else FResultCode:=Wrapper.LinearFit(NumPoints, FParamCount, FitType,
           @DX[0], @DY[0], @Res[0], FDeviation, Va);
  FParameters:=VarArrayCreate([0, FParamCount-1], varVariant);
  for I:=0 to FParamCount-1 do
  FParameters[I]:=Res[I];
  Result:=FResultCode={frOk}0; // result=true means OK
end;

function TDMFitter.LMFit: WordBool;
var
  rX, rY, rW: array of TReal;
  Chi, ParDel, ChiDel, Deriv: TReal;
  I, NumPoints, EventMode: integer;
  WP: pointer;
begin
  Result:=false; // failure by default
  FResultCode:=3{frBadParams};
  Fitter:=Unassigned;
  // check input and initialize parameters
  if not (FWeightType in [0..2])
  then Exit;  // invalid weight
  if (FIterations<1) or (FIterations>10000)
  then Exit;
  if (FParamCount<1) or (FParamCount>25)
  then Exit;
  EventMode:=0; // meaningful values: -1 - both, -2 - numeric differentiation
  if VarIsStr(FExpression)
  then StrPLCopy(FStrExpr, FExpression, SizeOf(FStrExpr)-1)
  else if VarIsType(FExpression, varDispatch) // fitter class?
       then Fitter:=VarAsType(FExpression, varDispatch)
       else if VarIsOrdinal(FExpression) // use events?
            then EventMode:=FExpression
            else Exit;
  if VarIsStr(FExpression) then // check parameters for multi-X expression
  begin
    UseXVector:=Pos('CX2', UpperCase(FStrExpr))>0;
    if UseXVector and (FXDimension=0)
    then Exit; 
  end; 
  if not (VarIsArray(FX) and VarIsArray(FY))
  then Exit;
  if FWeightType>0 then if not VarIsArray(FWeight)
  then Exit;
  NumPoints:=VarArrayHighBound(FX,1)+1;
  if (VarArrayHighBound(FY,1)+1<>NumPoints) or
     ((FWeightType>0) and (VarArrayHighBound(FWeight,1)+1<>NumPoints))
  then Exit;
  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;
  for I:=0 to VarArrayHighBound(FParameters,1) do
  begin
    Pars[I+1]:=FParameters[I];
    Sig[I+1]:=FSigmas[I];
  end;
  if VarIsArray(FOptions) and (VarArrayHighBound(FOptions,1)=2) then
  begin
    ParDel:=FOptions[0];
    ChiDel:=FOptions[1];
    Deriv:=FOptions[2];
  end else
  begin
    ParDel:=0;
    ChiDel:=0;
    Deriv:=0;
  end;
  // copy data
  SetLength(rX, NumPoints);
  SetLength(rY, NumPoints);
  if FWeightType>0 then
  begin
    SetLength(rW, NumPoints);
    WP:=@rW[0];
  end else
    WP:=nil;
  FResultCode:={frError}1;
  for I:=0 to NumPoints-1 do
  begin
    if FXDimension>0
    then rX[I]:=I // remember row INDEX in FX[..] instead of X!!!
    else rX[I]:=FX[I];
    rY[I]:=FY[I];
    if FWeightType>0
    then rW[I]:=FWeight[I];
  end;
  DMFitter:=self; // assign global reference used in callbacks
  try
    case EventMode of
       0: if VarIsEmpty(Fitter) then
            if UseXVector
            then FResultCode:=LMNLSFEx(NumPoints, FParamCount, FWeightType,
              @rX[0], @rY[0], WP, Pars, Sig, FIterations,
              LMFitProgress, LMFitFunc3, LMFitDeriv3, Chi, ParDel, ChiDel)
            else FResultCode:=LMNLSF(FStrExpr, NumPoints, FWeightType, @rX[0], @rY[0],
              WP, Pars, Sig, FIterations, LMFitProgress, Chi, ParDel, ChiDel, Deriv)
          else FResultCode:=LMNLSFEx(NumPoints,