unit Fitteru;
{
  Curve fitting using:
    - point-to-point linear aproximation,
    - linear, quadratic and cubic regression

  Usage:
     1. Get the wanted Fitter using Fitter(ACurveFitType) function
     2. Assign PointCount property (number of exp. points to be added)
        Eventually set OriginUse property
     3. Call PointAdd for each point (i.e. PointCount times)
     4. If PointAdd called PointCount times then
          CalcParam is called automatically when the last point added
        else
          Call CalcParams explicitely
        to calculate parameters for wanted FitType function
     5. The fitting function Func is now ready to be used

   Example:
     f := Fitter(cfLinear);
     f.PointCount = 200;
     for i := 0 to f.PointCount - 1 do
       f.PointAdd(x[i], y[i])

     // now e.g. draw graph of the function using
     // y = f.Func(x) function call
}

interface
uses
  Classes,
  fmath, matrices, regress,
  UlanType ;

type

  TFitter = class(TObject)
  private
    FFitResult: integer; { Result of callint regress.XXXFit function }
    FV: TMatrix;
    FX: TVector;{ array of x values to be used for fitting/regression analysis }
    FY: TVector;{ array of y values .. }
    FW: TVector; { weights; (used for originuse CurvePassThrough) }
    FDimPointCount: integer; { requested size of X, Y arrays, assigned in set PointCount
      to APointCount + 1 (for eventual 0,0 point) }
    FPointCount: integer; { number of experimental points requested/added,
      assigned in set PointCount to APointCount , eventually lowered in
      CalcParams method if less points added using PointAdd function, or increased
      by 1 if zero point added }
    FIndex: integer; { Index of the position where next point should be added
      in PointAdd method. Reset to 0 in set PointCount, increased by
      each PointAdd call. }
    FB: TVector; { The resulting parameters of the regr. function calculated
      in CalcParams method. }
    FBCount: integer; { number of B parameters }
    FFunc: TFunc;
    FFitType: TCurveFitType;
    FDeg: integer; { Degree of polynomial regression }
    FOriginUse: TOriginUse; { how to deal with x,y=0,0: ouIgnore, ouComputeWith,
     ouCurvePassThrough }
    FOriginIndex: integer; { -1 if no X[i]=0 point found or added } 
{ucprtype}
    procedure SetPointCount(APointCount: integer);
    procedure SetFitType(AFitType: TCurveFitType);
  protected
    property FitType: TCurveFitType read FFitType write SetFitType;
    property Deg: integer read FDeg write FDeg;
    function DoFunc(AX: Float): Float; virtual; abstract;
      { called from non-object function assigned to Func property }
  public
    constructor Create; reintroduce; virtual;
    destructor Destroy; override;
    procedure PointAdd(AX: Float; AY: Float);
    { Calculate the function parameters - do the regression. Call inherited
      in overloaded methods as the first statement. }
    procedure CalcParams; virtual;

    { If this property changed after CalcParams, call the CalcParams explicitely
      to take the new OriginUse value into consideration }
    property OriginUse: TOriginUse read FOriginUse write FOriginUse;
    { Number of experimental points in X and Y arrays. Should be assigned
      first, then given number of PointAdd should follow (the last call will
      invoke CalcParams method). Then the Func can be used. }
    property PointCount: integer read FPointCount write SetPointCount;
    { Function(X: Float): Float (corresponding to FitType) that uses
      the coeficients calculated in CalcParams and stored in B array }
    property X: TVector read FX;
    { Array of Y coordinates of experimental points }
    property Y: TVector read FY;
    { Array of X coordinates of experimental points }
    property BCount: integer read FBCount;
    { Array of Func function coeficients. }
    property B: TVector read FB;
    { Number of Func function coeficients. }
    property Func: TFunc read FFunc;
  end;

  TFitterClass = class of TFitter;


function Fitter(cf: TCurveFitType): TFitter;

implementation
var
  FFitters: array[TCurveFitType] of TFitter;

type
  TNoneFitter = class(TFitter)
    constructor Create; override;
    procedure CalcParams; override;
    function DoFunc(AX: Float): Float; override;
  end;

  TPTPFitter = class(TFitter)
  public
    constructor Create; override;
    procedure CalcParams; override;
    function DoFunc(AX: Float): Float; override;
  end;

  TLinearFitter = class(TFitter)
  public
    constructor Create; override;
    procedure CalcParams; override;
    function DoFunc(AX: Float): Float; override;
  end;

  TQuadraticFitter = class(TFitter)
  public
    constructor Create; override;
    procedure CalcParams; override;
    function DoFunc(AX: Float): Float; override;
  end;

  TCubicFitter = class(TFitter)
  public
    constructor Create; override;
    procedure CalcParams; override;
    function DoFunc(AX: Float): Float; override;
  end;

function NoneFunc(X: Float): Float;
begin
  Result := 0;
end;

function PointToPoint(X: Float): Float;
begin
  Result := Fitter(cfPointToPoint).DoFunc(X)
end;

function Linear(X: Float): Float;
{var f: TFitter;}
begin
  Result := Fitter(cfLinear).DoFunc(X);
end;

function Quadratic(X: Float): Float;
{var f: TFitter;}
begin
  Result := Fitter(cfQuadratic).DoFunc(X);
{  f := Fitter(cfQuadratic);
  Result := f.B[0] + f.B[1] * X + f.B[2] * X * X;}
end;

function Cubic(X: Float): Float;
{var
  f: TFitter;}
begin
  Result := Fitter(cfCubic).DoFunc(X);
{
  f := Fitter(cfCubic);
  Result := f.B[0] + f.B[1] * X + f.B[2] * X * X +
    f.B[3] * X * X * X;}
end;

{TFitter.}
constructor TFitter.Create;
begin
  inherited;
end;

procedure TFitter.SetFitType(AFitType: TCurveFitType);
begin
  FFitType := AFitType;
  case FFitType of
    cfPointToPoint: begin
      FBCount := 0;
      FFunc := PointToPoint;
      FDeg := 1;
    end;
    cfLinear: begin
      FBCount := 2;
      FFunc := Linear;
      FDeg := 1;
    end;
    cfQuadratic: begin
      FBCount := 3;
      FFunc := Quadratic;
      FDeg := 2;
    end;
    cfCubic: begin
      FBCount := 4;
      FFunc := Cubic;
      FDeg := 3;
    end;
  else
    FFunc := nil;
    FBCount := 0;
    FDeg := 0;
  end;
  if FBCount > 0 then
    DimVector(FB, FBCount - 1);
end;

procedure TFitter.SetPointCount(APointCount: integer);
var
  aDimPointCount: integer;
begin
  aDimPointCount := APointCount + 1;
  if aDimPointCount > FDimPointCount then begin
    DimVector(FX, aDimPointCount); { dim on more slot for eventual inclusion of 0,0 point }
    DimVector(FY, aDimPointCount);
    DimMatrix(FV, Max(aDimPointCount, Deg), Max(aDimPointCount, Deg));
    FDimPointCount := aDimPointCount;
  end;
  FPointCount := APointCount;
  FIndex := 0;
  FOriginIndex := -1;
end;

procedure TFitter.PointAdd(AX: Float; AY: Float);
begin
  if FIndex = FPointCount then
    exit; { Ignore adding of more points than before assigned PointCount }
  FX[FIndex] := AX;
  FY[FIndex] := AY;
  inc(FIndex);
  if FIndex = FPointCount then
    CalcParams;
end;

{ Insert AValue at AIndex into AVector, that has ALength element. Returns true
  if OK and increases ALength by one. I.e. ALength can be less then Length(AVector). }
function VectorInsert(AVector: TVector; var ALength: integer; AValue: Float; AIndex: integer): boolean;
var i: integer;
begin
  Result := false;
  if AIndex < 0 then
    exit;
  if Length(AVector) < ALength + 1 then
    SetLength(AVector, ALength + 1);
  for i := ALength - 1 downto AIndex do begin
    AVector[i + 1] := AVector[i];
  end;
  AVector[AIndex] := AValue;
  inc(ALength);
  Result := true;
end;

procedure TFitter.CalcParams;
var
  i, j: integer;
begin
  FPointCount := FIndex;
  Regress.SetFirstPoint(0);// set first data point index to 0 for fitting routines
  if (OriginUse <> ouIgnore) and (FOriginIndex = -1) then begin
    if FPointCount > 0 then begin
      for i := 0 to FPointCount - 1 do begin
        if FX[i] = 0 then begin
          FOriginIndex := i;
          break;
        end else begin
          if FX[i] > 0 then begin
            j := FPointCount;
            FOriginIndex := i;
            VectorInsert(FX, j, 0, FOriginIndex);
            VectorInsert(FY, FPointCount, 0, FOriginIndex);
            break;
          end;
        end;
      end;
      if FOriginIndex = -1 then begin
        FOriginIndex := FPointCount;
        inc(FPointCount);
        FX[FOriginIndex] := 0;
        FY[ForiginIndex] := 0;
      end;
    end;
    if OriginUse = ouCurvePassThrough then begin
      if (FW = nil) or (Length(FW) < FPointCount) then begin
        DimVector(FW, FPointCount);
      end;
      for i := 0 to FPointCount - 1 do begin
        FW[i] := 1;
      end;
      FW[FOriginIndex] := 100;
    end;
  end;
end;

destructor TFitter.Destroy;
begin
  FB := nil;
  FX := nil;
  FY := nil;
  FV := nil;
  FW := nil;
  inherited;
end;
{/TFitter.}

{TNoneFitter.}

constructor TNoneFitter.Create;
begin
  inherited;
  FitType := cfNone;
  FFunc := NoneFunc;
end;

procedure TNoneFitter.CalcParams;
begin
  inherited;
end;

function TNoneFitter.DoFunc(AX: Float): Float;
begin
  Result := 0;
end;
{/TNoneFitter.}

{TPTPFitter.}
constructor TPTPFitter.Create;
begin
  inherited;
  FitType := cfPointToPoint;
  //FFunc := PointToPoint;
end;

procedure TPTPFitter.CalcParams;
begin
  inherited;
end;

function TPTPFitter.DoFunc(AX: Float): Float;
var
  i, i1, i2: integer;
begin
  i1 := -1;
  for i := 0 to PointCount - 1 do begin
    { find the interval to which the X belong }
    if AX <= X[i] then begin
      if i = 0 then begin
        { if below first point, extrapolate using first two points }
        i1 := i;
        i2 := i + 1;
      end else begin
        i1 := i - 1;
        i2 := i;
      end;
      break;
    end;
  end;
  if i1 = -1 then begin
    { X is more then last point's X, or there are no points }
    if PointCount >= 2 then begin
      { extrapolate using last two points }
      i1 := PointCount - 2;
      i2 := PointCount - 1;
    end else begin
      { if there is less then two points return 0 }
      Result := 0;
      exit;
    end;
  end;
  if X[i2] = X[i1] then begin
    { division by zero would result, return 0 }
    Result := 0;
    exit;
  end;
  Result :=
    Y[i1] + { y offset }
    (Y[i2] - Y[i1])/(X[i2] - X[i1])
    * (AX - X[i1]);
end;
{/TPTPFitter.}

{TLinearFitter.}
constructor TLinearFitter.Create;
begin
  inherited;
  FitType := cfLinear;
end;

function TLinearFitter.DoFunc(AX: Float): Float;
begin
  Result := B[0] + B[1] * AX;
end;

procedure TLinearFitter.CalcParams;
begin
  inherited;
  if OriginUse = ouCurvePassThrough then begin
    FFitResult := WLinFit(FX, FY, FW, FPointCount, FB, FV);
  end else begin
    FFitResult := LinFit(FX, FY, FPointCount - 1, FB, FV);
  end;
end;
{/TLinearFitter.}

{TQuadraticFitter.}
constructor TQuadraticFitter.Create;
begin
  inherited;
  FitType := cfQuadratic;
  FFunc := Quadratic;
end;

procedure TQuadraticFitter.CalcParams;
begin
  inherited;
  if OriginUse = ouCurvePassThrough then begin
    FFitResult := WPolFit(FX, FY, FW, FPointCount - 1, 2, FB, FV);
  end else begin
    FFitResult := PolFit(FX, FY, FPointCount - 1, 2, FB, FV);
  end;
end;

function TQuadraticFitter.DoFunc(AX: Float): Float;
begin
  Result := B[0] + B[1] * AX + B[2] * AX * AX;
end;
{/TQuadraticFitter.}

{TCubicFitter.}
constructor TCubicFitter.Create;
begin
  inherited;
  FitType := cfCubic;
  FFunc := Linear;
end;

procedure TCubicFitter.CalcParams;
begin
  inherited;
  if OriginUse = ouCurvePassThrough then begin
    FFitResult := WPolFit(FX, FY, FW, FPointCount - 1, 3, FB, FV);
  end else begin
    FFitResult := PolFit(FX, FY, FPointCount - 1, 3, FB, FV);
  end;
end;

function TCubicFitter.DoFunc(AX: Float): Float;
begin
  Result := B[0] + B[1] * AX + B[2] * AX * AX + B[3] * AX * AX * AX;
end;
{/TCubicFitter.}

const
  FFitterClasses: array[TCurveFitType] of TFitterClass =
  (TNoneFitter, TPTPFitter, TLinearFitter, TQuadraticFitter, TCubicFitter);

function Fitter(cf: TCurveFitType): TFitter;
var f: TFitter;
begin
  if FFitters[cf] = nil then begin
    FFitters[cf] := FFitterClasses[cf].Create;
  end;
  Result := FFitters[cf];
end;

procedure FittersClear;
begin
  FillChar(FFitters, sizeof(FFitters), 0);
end;

procedure FittersFree;
var
  i: TCurveFitType;
begin
  for i := low(TCurveFitType) to high(TCurveFitType) do begin
    FFitters[i].Free;
    FFitters[i] := nil;
  end;
end;

initialization
  FittersClear;

finalization
  FittersFree;

end.
