unit Fittingu;
{
  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)
     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
    FX: TVector;{ array of x values to be used for fitting/regression analysis }
    FY: TVector;{ array of y values .. }
    FDimPointCount: integer; { requested size of X, Y arrays, assigned in set PointCount }
    FPointCount: integer; { number of experimental points requested/added,
      assigned in set PointCount, eventually lowered in CalcParams method
      if less points added using PointAdd function }
    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;

    procedure SetPointCount(APointCount: integer);
    procedure SetFitType(AFitType: TCurveFitType);
  protected
    property FitType: TCurveFitType read FFitType write SetFitType;
    function DoFunc(X: Float): Float; { called from non-object function
      assigned to Func property }
    property B: TVector read FB;
    property BCount: integer read FBCount;
    property X: TVector read FX;
    property Y: TVector read FY;
  public
    destructor Destroy; override;
    procedure PointAdd(AX: Float; AY: Float);
    { Calculate the function parameters - do the regression. }
    procedure CalcParams; virtual; abstract;
    property PointCount: integer read FPointCount write SetPointCount;
    property Func: TFunc read FFunc;
  end;

function Fitter(cf: TCurveFitType): TFitter;

implementation
var
  FFitters: array[TCurveFitType] of TFitter;

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

function Linear(X: Float): Float;
var f: TFitter;
begin
  f := Fitter(cfLinear);
  Result := f.B[0] + f.B[1] * X;
end;

function Quadratic(X: Float): Float;
var f: TFitter;
begin
  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
  f := Fitter(cfCubic);
  Result := f.B[0] + f.B[1] * X + f.B[2] * X * X +
    f.B[3] * X * X * X;
end;

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

procedure TFitter.SetPointCount(APointCount: integer);
begin
  FPointCount := APointCount;
  FDimPointCount := FPointCount;
  DimVector(FX, PointCount - 1);
  DimVector(FY, PointCount - 1);
  FIndex := 0;
end;

procedure TFitter.PointAdd(AX: Float; AY: Float);
begin
  FX[FIndex] := AX;
  FY[FIndex] := AY;
  inc(FIndex);
  if FIndex = FDimPointCount then
    CalcParams;
end;

procedure TFitter.CalcParams;
begin
  FPointCount := FIndex;
  {...}
end;

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

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

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

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

initialization
  FittersClear;

finalization
  FittersFree;

end.
