program fitdemo;
(*  compile this demo program from the command-line as
    DCC32 -u"C:\OPTIVEC\LIB4D" FITDEMO.DPR
    or open as a project in the Delphi IDE
    and add C:\OPTIVEC\LIB4D to the search path in
    Project/Options/Directories and Conditionals
    For multi-core support, write "LIB7M" instead of "LIB4D"

    Copyright 1996-2008 OptiCode - Dr. Martin Sander Software Dev.
*)

{$N+} {$D+}
uses SysUtils, WinTypes, WinProcs, Messages, Classes, Graphics, Controls,
     Forms, Dialogs, ExtCtrls, StdCtrls,
     VecLib, VDstd, VDmath, VIstd, MDstd, VDnlfit, VDmnlfit, Vgraph;

type
      TForm1 = class(TForm)
        PaintBox1: TPaintBox;
        BNextView: TButton;
        procedure ShowView(Sender: TObject);
        procedure NextView(Sender: TObject);
    end;


var
   Form1:                      TForm1;
   XExp, YExp, YFit,
   YExp2, YFit2, YExp3, YFit3: dVector;
   FitPars:                    array[0..6] of double;
       (* the highest number of parameters we will have in our examples *)
   ParStatus:                  array[0..6] of Integer;
   Opt:                        VD_NONLINFITOPTIONS;
   ExpList:                    array[0..2] of VD_EXPERIMENT;
   DataText:                   string[20];
   vview:                      Integer;
const
   ProgName    = 'OptiVec Data-Fitting Demo';
   polydeg     = 5;     (* refers to the polynomial fitting example *)
   sizex: UInt = 200;
   {$R *.DFM}

    (* the following function is used with VD_linfit: *)
procedure PolyModel( BasFuncs:dVector; x:double; nfuncs:UInt );
var i:UInt;
begin
   (* This function fills the vector BasFuncs with powers of x.
      VD_linfit will then determine the coefficient for each of these
      powers.
      Note that the coefficients do not occur in the model function!
      The basis functions with known coefficients (whose fitting has
      been disabled) are still part of the model and must be calculated.
      You will see below that the values of the known (disabled)
      coefficients  must be set prior to calling VD_linfit.          *)
     BasFuncs^ := 1.0;
     for i:=1 to nfuncs-1 do
          VD_Pelement( BasFuncs,i )^ := VD_element( BasFuncs, i-1 ) * x;
end;

   (* the following function is used with VD_nonlinfit: *)
procedure VPolyModel( Y, X:dVector; size:UInt );
    (* Here, the model function has to fill a whole result vector,
       using your first guess of FitPars. In contrast to the
       linear case, now the coefficients are explicitly used in the
       model function. You must initialize FitPars with something,
       even if you have no idea about the result.
       FitPars must be global, so that the model function can access the
       parameters. With VD_nonlinfit, you can use just any functions
       in your model.
       For better comparison, we use the same polynomial approximation
       as before, but now we code it as if we didn't know that a
       polynomial actually is linear in its coefficients (and as if
       we didn't know either how to efficiently code a polynomial at all). *)
var i:  UInt;
    xi: double;
begin
    for i:=0 to size-1 do
    begin
        xi  := VD_element( X, i );
        VD_Pelement( Y, i )^ := FitPars[0]
            + FitPars[1] * xi
            + FitPars[2] * xi * xi
            + FitPars[3] * xi * xi * xi
            + FitPars[4] * xi * xi * xi * xi
            + FitPars[5] * xi * xi * xi * xi * xi;
                  (* FitPars is a static array and can be accessed
                     with the [ ] operator  *)
   end;
end;

   (* the following function is used with VD_multiNonlinfit: *)
procedure VSineModel( Y, X:dVector; size, iExperiment:UInt );
    (* According to the parameter iExperiment, the model function must
       choose the correct parameters for the calculation of the model
       function. The model function itself is the same for all experiments. *)
var omega, phase, amp: double;
begin
    case iExperiment of
        0: begin
               phase := FitPars[1];
               amp   := FitPars[4];
           end;
        1: begin
               phase := FitPars[2];
               amp   := FitPars[5];
           end;
        2: begin
               phase := FitPars[3];
               amp   := FitPars[6];
           end;
    end;
    omega := FitPars[0];
              (* we assume this parameter to be the same for all *)
    VDx_sin( Y, X, size, omega, phase, amp );
end;

procedure TForm1.NextView(Sender: TObject);
begin
       inc(vview);
       PaintBox1.Invalidate;
end;

procedure TForm1.ShowView(Sender: TObject);
var   StrDum:       array[0..80] of char;
      vDC:          HDC;
begin
   V_initPlot( PaintBox1 );
   with PaintBox1 do
   begin
       vDC := Canvas.Handle;
       V_setPlotRegion( Left + Width div 4, Top, Left+Width, Top+Height-Height div 16 );
   end;
   case vview  of
       0: begin   TextOut( vDC, 10, 10,
                  'This is a series of graphs illustrating OptiVec data-fitting functions.', 71 );
                  TextOut( vDC, 10, 40,
                  'Press [Alt] [F4] to end the demonstration.', 41 );
          end;
       1: begin   TextOut( vDC, 2, 10, 'Suppose these are', 17 );
                  TextOut( vDC, 2, 30, 'your experimental', 17 );
                  TextOut( vDC, 2, 50, 'data points.', 12 );
                  TextOut( vDC, 2, 70, '(Actually, they consist', 23 );
                  TextOut( vDC, 2, 90, 'of a simple cubic with', 22 );
                  TextOut( vDC, 2,110, '1% added noise.)', 16 );

                  VD_ramp( XExp, sizex, 0, 1.0/(sizex-1) ); (* 'experimental' x-axis from 0 to 1 *)
                  VD_cubic( YExp, XExp, sizex );        (* fake 'measured' y-data as y = x^3 *)
                  VD_noise( YFit, sizex, 1, 0.005 );
                  VD_addV( YExp, YExp, YFit, sizex ); (* add 1% peak-to-peak 'experimental noise' *)
                  VD_xyAutoPlot( XExp, YExp, sizex, PS_NULL+SY_CROSS, GREEN );
          end;
       2: begin  // fit your data to one of the simplest models, a polynomial
                  TextOut( vDC, 2, 10, 'The red curve is a', 18 );
                  TextOut( vDC, 2, 30, 'fifth-order polynomial', 22 );
                  TextOut( vDC, 2, 50, 'fitted to your data.', 20 );
                  TextOut( vDC, 2, 70, 'Without noise, the', 18 );
                  TextOut( vDC, 2, 90, 'coefficients should', 19 );
                  TextOut( vDC, 2,110, 'have been:', 10 );
                  TextOut( vDC, 2,130, '{0, 0, 0, 1.0, 0, 0}', 20 );

                  V_setBoundaryCheckWarn( FALSE );
                  (* Using static array FitPars; don't want to see warnings about that *)
                  VD_polyfit( @FitPars, polydeg, XExp, YExp, sizex );
                  VD_poly( YFit, XExp, sizex, @FitPars, polydeg ); // calculate fit curve
                  TextOut( vDC, 2,160, 'Actually, we got:', 17 );
                  Str( FitPars[0]:8:5, DataText );
                  TextOut( vDC, 2,178, StrPCopy(StrDum, 'a0 = ' + DataText), 13 );
                  Str( FitPars[1]:8:5, DataText );
                  TextOut( vDC, 2,196, StrPCopy(StrDum, 'a1 = ' + DataText), 13 );
                  Str( FitPars[2]:8:5, DataText );
                  TextOut( vDC, 2,214, StrPCopy(StrDum, 'a2 = ' + DataText), 13 );
                  Str( FitPars[3]:8:5, DataText );
                  TextOut( vDC, 2,232, StrPCopy(StrDum, 'a3 = ' + DataText), 13 );
                  Str( FitPars[4]:8:5, DataText );
                  TextOut( vDC, 2,250, StrPCopy(StrDum, 'a4 = ' + DataText), 13 );
                  Str( FitPars[5]:8:5, DataText );
                  TextOut( vDC, 2,268, StrPCopy(StrDum, 'a5 = ' + DataText), 13 );
                  TextOut( vDC, 2,290, 'Note how even moderate', 22 );
                  TextOut( vDC, 2,308, 'noise leads to rather', 21 );
                  TextOut( vDC, 2,326, 'large errors in the', 19 );
                  TextOut( vDC, 2,344, 'fit parameters, if', 18 );
                  TextOut( vDC, 2,362, 'there are too many', 18 );
                  TextOut( vDC, 2,380, '"free" parameters.', 18 );
                  VD_xy2AutoPlot( XExp, YExp, sizex, PS_NULL or SY_CROSS, GREEN,
                                  XExp, YFit, sizex, PS_SOLID, LIGHTRED );
           end;
       3: begin // refine your fit by switching to a general linear model,
                        // giving you the chance to consider only the uneven terms
                  TextOut( vDC, 2, 10, 'Suppose you know that', 21 );
                  TextOut( vDC, 2, 28, 'the coefficients of', 19 );
                  TextOut( vDC, 2, 46, 'all even terms are 0.', 21 );
                  TextOut( vDC, 2, 64, 'Then you fit to your', 20 );
                  TextOut( vDC, 2, 82, 'own linear model,', 17 );
                  TextOut( vDC, 2,100, 'consisting only of', 18 );
                  TextOut( vDC, 2,118, 'uneven terms.', 13 );
                  TextOut( vDC, 2,140, 'Now we get:', 11 );

                  ParStatus[0] := 0;
                  ParStatus[2] := 0;
                  ParStatus[4] := 0;  (* disable fitting of even terms *)
                  FitPars[0] := 0.0;  (* the disabled fitting parameters must be *)
                  FitPars[2] := 0.0;  (*  initialized before calling VD_linfit ! *)
                  FitPars[4] := 0.0;  (* set them to the known value, 0.0 *)
                  ParStatus[1] := 1;
                  ParStatus[3] := 1;
                  ParStatus[5] := 1;  (* enable fitting of uneven terms *)
                  VD_linfit( @FitPars, @ParStatus, polydeg+1,
                             XExp, YExp, sizex,
                             @PolyModel );
                  VD_poly( YFit, XExp, sizex, @FitPars, polydeg );
                                 (* calculate new fit curve  *)
                  TextOut( vDC, 2,164, 'a0 =  0 (fix)', 13 );
                  Str( FitPars[1]:8:5, DataText );
                  TextOut( vDC, 2,182, StrPCopy(StrDum, 'a1 = ' + DataText), 13 );
                  TextOut( vDC, 2,200, 'a2 =  0 (fix)', 13 );
                  Str( FitPars[3]:8:5, DataText );
                  TextOut( vDC, 2,218, StrPCopy(StrDum, 'a3 = ' + DataText), 13 );
                  TextOut( vDC, 2,236, 'a4 =  0 (fix)', 13 );
                  Str( FitPars[5]:8:5, DataText );
                  TextOut( vDC, 2,254, StrPCopy(StrDum, 'a5 = ' + DataText), 13 );
                  TextOut( vDC, 2,280, 'This is about as close', 22 );
                  TextOut( vDC, 2,298, 'as we can get in the', 20 );
                  TextOut( vDC, 2,316, 'presence of noise.', 18 );
                  VD_xy2AutoPlot( XExp, YExp, sizex, PS_NULL or SY_CROSS, GREEN,
                                  XExp, YFit, sizex, PS_SOLID, LIGHTRED );
           end;
       4: begin  // here, we mis-use a non-linear fitting algorithm
                        // for our simple problem.
                  TextOut( vDC, 2, 10, 'Let us fire with the', 19 );
                  TextOut( vDC, 2, 30, '"cannon" of a non-', 18 );
                  TextOut( vDC, 2, 50, 'linear fit on our', 17 );
                  TextOut( vDC, 2, 70, 'simple "sparrow"', 16 );
                  TextOut( vDC, 2, 90, 'problem.', 8 );
                  TextOut( vDC, 2,110, 'It takes much longer', 20 );
                  TextOut( vDC, 2,130, 'to find the result...', 21 );

                  ParStatus[0] := 0;
                  ParStatus[2] := 0;
                  ParStatus[4] := 0;  (* disable fitting of even terms *)
                  ParStatus[1] := 1;
                  ParStatus[3] := 1;
                  ParStatus[5] := 1;  (* enable fitting of uneven terms *)
                  FitPars[0] := 0.0;
                  FitPars[2] := 0.0;
                  FitPars[4] := 0.0;  (* set known coefficients to the value, 0.0, as before *)
                  FitPars[1] := 1.5;
                  FitPars[3] := 1.5;
                  FitPars[5] := 1.5;  (* you must provide some guess here!
                       all fitting parameters must be initialized before calling VD_nonlinfit ! *)
                  VD_getNonlinfitOptions( Opt );
                  Opt.FigureOfMerit := 0;  (* choose least-square fitting *)
                  Opt.AbsTolChi     := 1.e-6;
                  Opt.FracTolChi    := 1.e-3; (* makes the fit fast, but not very accurate *)
                  Opt.LevelOfMethod := 3;  (* if you fear you might jump into a
                        local rather than the true global parameter optimum, try
                        LevelOfMethod := 7 - but only if you have time to wait for the result *)
                  VD_setNonlinfitOptions( Opt );

                  VD_nonlinfit( @FitPars, @ParStatus, polydeg+1,
                                XExp, YExp, sizex,
                                @VPolyModel, nil );
                  (* If you know the derivatives with respect to each parameter, put
                     your knowledge into a DerivModel function and replace the 'nil'
                     parameter with its address.
                     (Actually, here we do know; but let's pretend we don't,
                     and have VD_nonlinfit call the numeric differentiation
                     procedure. *)
                  VPolyModel( YFit, XExp, sizex );  (* get fit curve from model *)

                  TextOut( vDC, 2,160, 'But finally we get:', 19 );
                  TextOut( vDC, 2,180, 'a0 =  0 (fix)', 13 );
                  Str( FitPars[1]:8:5, DataText );
                  TextOut( vDC, 2,200, StrPCopy(StrDum, 'a1 = ' + DataText), 13 );
                  TextOut( vDC, 2,220, 'a2 =  0 (fix)', 13 );
                  Str( FitPars[3]:8:5, DataText );
                  TextOut( vDC, 2,240, StrPCopy(StrDum, 'a3 = ' + DataText), 13 );
                  TextOut( vDC, 2,260, 'a4 =  0 (fix)', 13 );
                  Str( FitPars[5]:8:5, DataText );
                  TextOut( vDC, 2,280, StrPCopy(StrDum, 'a5 = ' + DataText), 13 );
                  TextOut( vDC, 2,310, 'That is virtually the', 21 );
                  TextOut( vDC, 2,330, 'same as before.', 15 );
                  VD_xy2AutoPlot( XExp, YExp, sizex, PS_NULL or SY_CROSS, GREEN,
                                  XExp, YFit, sizex, PS_SOLID, LIGHTRED );
           end;
       5: begin  (* finally, let's suppose you have several experimental
                    curves, measuring the same physical process under
                    slightly different conditions.
                    Say, we have a vibration, and each measurement
                    begins with a different phase and has a somewhat
                    different amplitude, but the same frequency.  *)

                  SetViewportOrgEx( vDC, 0, 0, nil );
                  TextOut( vDC, 2, 10, 'Finally, let us fit', 19 );
                  TextOut( vDC, 2, 28, 'several sets of experi-', 23 );
                  TextOut( vDC, 2, 46, 'mental data at once', 19 );
                  TextOut( vDC, 2, 64, '(sine waves with the', 20 );
                  TextOut( vDC, 2, 82, 'same frequency, but', 19 );
                  TextOut( vDC, 2,100, 'different phases and', 20 );
                  TextOut( vDC, 2,118, 'amplitudes):', 12 );
                  TextOut( vDC, 2,136, 'First you see the', 17 );
                  TextOut( vDC, 2,154, '"experimental" data', 19 );
                  TextOut( vDC, 2,172, 'Please wait...', 14 );

                  VD_ramp( XExp, sizex, 0, 1.0/(sizex-1) ); (* x-axis again from 0 to 1 *)
                  VDx_sin( YExp,  XExp, sizex, 15.0,  0.0, 1.2 );  (* first 'measurement' *)
                  VDx_sin( YExp2, XExp, sizex, 15.0,  0.5, 1.0 );  (* second 'measurement'*)
                  VDx_sin( YExp3, XExp, sizex, 15.0, -1.8, 0.75 ); (* third 'measurement' *)

                  VD_xy2AutoPlot( XExp, YExp,  sizex, PS_NULL + SY_CROSS, GREEN,
                                  XExp, YExp2, sizex, PS_NULL + SY_CROSS, BLUE );
                  VD_xyDataPlot(  XExp, YExp3, sizex, PS_NULL + SY_CROSS, MAGENTA );

                    (* cram your experiments into the array of VD_EXPERIMENT structs: *)
                  ExpList[0].X := XExp;   ExpList[0].Y := YExp;   ExpList[0].size := sizex;
                  ExpList[1].X := XExp;   ExpList[1].Y := YExp2;  ExpList[1].size := sizex;
                  ExpList[2].X := XExp;   ExpList[2].Y := YExp3;  ExpList[2].size := sizex;
                         (* we are not using the InvVar and WeightOfExperiment fields
                            of ExpList, as we are not weighting the data. *)

                  VI_equC( @ParStatus, 7, 1 ); (* we have 1 frequency, 3 phases, and 3 amplitudes,
                                 and all these 7 parameters are unknown.
                                 We must provide a first guess for each of them: *)
                  FitPars[0] := 10.0;  (* the frequency term *)
                  FitPars[1] :=  0.0;
                  FitPars[2] :=  0.0;
                  FitPars[3] :=  0.0;  (* the three phases  *)
                  FitPars[4] :=  1.5;
                  FitPars[5] :=  1.5;
                  FitPars[6] :=  1.5;  (* the three amplitudes *)
                  VD_getNonlinfitOptions( Opt );
                  Opt.AbsTolChi  := 1.e-8;
                  Opt.FracTolChi := 1.e-6; (* force higher accuracy to avoid premature break-off *)
                    (*  Unlike almost every other fitting routine available, you
                        can get a result even for input parameters much farther off
                        from the true value than the 'guesses' chosen above.
                        But then you must run VD_multiNonlinfit at 'full power' and
                        enable the following line:                  *)
                    (* Opt.LevelOfMethod = 7; *)
                  VD_setNonlinfitOptions( Opt );

                  VD_multiNonlinfit( @FitPars, @ParStatus, 7,
                                     @ExpList, 3,
                                     @VSineModel,
                                     nil ); (* Again, we pretend we don't know the derivatives
                               bring the phases into the range -PI < phase < + PI *)
                  VD_modC( @(FitPars[1]), @(FitPars[1]), 3, 2.0*PI );
                  if FitPars[1] > PI then FitPars[1] := FitPars[1] - 2.0*PI
                  else if FitPars[1] < -PI then FitPars[1] := FitPars[1] + 2.0*PI;
                  if FitPars[2] > PI then FitPars[2] := FitPars[2] - 2.0*PI
                  else if FitPars[2] < -PI then FitPars[2] := FitPars[2] + 2.0*PI;
                  if FitPars[3] > PI then FitPars[3] := FitPars[3] - 2.0*PI
                  else if FitPars[3] < -PI then FitPars[3] := FitPars[3] + 2.0*PI;
                  VSineModel( YFit,  XExp, sizex, 0 );  (* get fit curves from your model *)
                  VSineModel( YFit2, XExp, sizex, 1 );
                  VSineModel( YFit3, XExp, sizex, 2 );

                  SetViewportOrgEx( vDC, 0, 200, nil );
                            (* go to lower part of text printing region *)
                  TextOut( vDC, 2,0, 'Here are the results', 20 );
                  TextOut( vDC, 2,18, '(in brackets: "true")', 21 );
                  Str( FitPars[0]:9:5, DataText );
                  TextOut( vDC, 2,36, StrPCopy(StrDum, 'freq =' + DataText + ' (15.0)'), 22 );
                  Str( FitPars[1]:9:5, DataText );
                  TextOut( vDC, 2,54, StrPCopy(StrDum, 'ph1  =' + DataText + ' (0.0)'), 21 );
                  Str( FitPars[2]:9:5, DataText );
                  TextOut( vDC, 2,72, StrPCopy(StrDum, 'ph2  =' + DataText + ' (0.5)'), 21 );
                  Str( FitPars[3]:9:5, DataText );
                  TextOut( vDC, 2,90, StrPCopy(StrDum, 'ph3  =' + DataText + ' (-1.8)'), 22 );
                  Str( FitPars[4]:9:5, DataText );
                  TextOut( vDC, 2,108, StrPCopy(StrDum, 'amp1 =' + DataText + ' (1.2)'), 21 );
                  Str( FitPars[5]:9:5, DataText );
                  TextOut( vDC, 2,126, StrPCopy(StrDum, 'amp2 =' + DataText + ' (1.0)'), 21 );
                  Str( FitPars[6]:9:5, DataText );
                  TextOut( vDC, 2,144, StrPCopy(StrDum, 'amp3 =' + DataText + ' (0.75)'), 22 );

                  V_continuePlot();   (* go back to plot viewport *)
                  VD_xyDataPlot(  XExp, YFit, sizex, PS_SOLID, LIGHTRED );
                  VD_xyDataPlot(  XExp, YFit2, sizex, PS_SOLID, LIGHTRED );
                  VD_xyDataPlot(  XExp, YFit3, sizex, PS_SOLID, LIGHTRED );

                  vview := -1;
           end;
   end;
end;

{$R *.RES}

begin
    V_initMT( 2 );
       (* multi-core initialization for 2 processors; no effect for single-thread *)
    XExp  := VD_vector( sizex );
    YExp  := VD_vector( sizex );
    YFit  := VD_vector( sizex );
    YExp2 := VD_vector( sizex );
    YExp3 := VD_vector( sizex );
    YFit2 := VD_vector( sizex );
    YFit3 := VD_vector( sizex );
    Application.Initialize;
    Application.CreateForm(TForm1, Form1);
    Application.Run;
    V_freeAll;
end.

