{ unité de fit des spectres 1D, 2D et 3D }
{ 100*(B-A*A)*(B-A*A)+(1-A)*(1-A) }
{ 100*exp(-sqr(X-64)/10000)*exp(-sqr(Y-64)/10000); }
{ ============================================================================= }
{ fenetre de fit de spectres 1D, 2D et 3D                                       }
{ 07.05.02 : FO : * possibilité de fitter des spectres 1D                       }
{                 * mise a jour automatique de la liste des spectres dès que le focus passe sur la fenetre de fit }
{                 * refresh automatique                                         }
{ ============================================================================= }
{ A FAIRE                                                                       }
{ * verifier la validité de la somme circulaire !!!                             }
{   la différence ne fonctionne pas avec le spectre "circulaire"                }
{ * ce n'est pas raisonable de fitter systématiquement avec un vecteur pars fixe de 5 paramètres !!?? }
{   il faut que cette taille soit variable                                      }
unit FitForm;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls, Parser10, Math387, optimization, EditReal, Menus,
  Spin, ExtCtrls, editData, MtxVec;

type
  TFittingForm = class(TForm)
    MyFunction: TEdit;
    Label1: TLabel;
    ParserFitFunction: TParser;
    MainMenu1: TMainMenu;
    Fitting1: TMenuItem;
    Simplex1: TMenuItem;
    GroupBox1: TGroupBox;
    Abox: TEditReal;
    Label5: TLabel;
    Label6: TLabel;
    Bbox: TEditReal;
    Label8: TLabel;
    Cbox: TEditReal;
    Label9: TLabel;
    Dbox: TEditReal;
    Ebox: TEditReal;
    Label10: TLabel;
    ChoixDim: TRadioGroup;
    GroupBox3: TGroupBox;
    ChoixSpectre2D: TComboBox;
    ChoixSpectre3D: TComboBox;
    ChoixSpectre1D: TComboBox;
    Label3: TLabel;
    Label7: TLabel;
    Label11: TLabel;
    GroupBox4: TGroupBox;
    Label2: TLabel;
    toleranceBox: TEditReal;
    Label4: TLabel;
    IterMaxBox: TSpinEdit;
    Label12: TLabel;
    EcartBox: TEditReal;
    CalculateDifference1: TMenuItem;
    ests1: TMenuItem;
    OutputSpectrum: TEdit;
    Label13: TLabel;
    FittingStatus: TLabel;
    BFGS1: TMenuItem;
    Marquardt1: TMenuItem;
    Label14: TLabel;
    MyFunctionSpectrum: TEdit;
    procedure Simplex1Click(Sender: TObject);
    procedure CalculateDifference1Click(Sender: TObject);
    procedure UpdateCombos1Click(Sender: TObject);
    procedure Updateavailablespectra1Click(Sender: TObject);
    procedure Plotdiff2D1Click(Sender: TObject);
    procedure FormCreate(Sender: TObject);
    procedure BFGS1Click(Sender: TObject);
    procedure Marquardt1Click(Sender: TObject);
    procedure FormActivate(Sender: TObject);
    procedure ChoixDimClick(Sender: TObject);
  private
    { Private declarations }
    procedure Plot1D;
    procedure PlotDiff1D;
    procedure Plot2D;
    procedure PlotDiff2D;
    procedure PlotDiff;
    procedure MiseAJour;
    function  CalculDifference1D(ind:integer):Tsample;
    function  CalculDifference2D(ind:integer):Tsample;
    function  CalculDifference3D(ind:integer):Tsample;
    procedure MiseAjourChoixSpectres;
    procedure fitterSimplex;
    procedure fitterBFGS;
    procedure fitterMarquardt;    
  public
    { Public declarations }
  end;

var
  FittingForm: TFittingForm;
  DimSpectrum,IndexSpectrum:integer;
  Pars : Array [0..4] of TSample;

implementation

uses MainForm;

{$R *.dfm}

procedure TFittingForm.MiseAJour;
begin
    ParserFitFunction.expression:=MyFunction.Text;
    { initial estimates for x1 and x2}
    Pars[0] := Abox.value;
    Pars[1] := Bbox.value;
    Pars[2] := Cbox.value;
    Pars[3] := Dbox.value;
    Pars[4] := Ebox.value;
    DimSpectrum:=ChoixDim.itemIndex+1;
    if DimSpectrum=1 then IndexSpectrum:=und.indexof(ChoixSpectre1D.Text);
    if DimSpectrum=2 then IndexSpectrum:=deuxd.indexof(ChoixSpectre2D.Text);
    if DimSpectrum=3 then IndexSpectrum:=troisd.indexof(ChoixSpectre3D.Text);
end;

procedure TFittingForm.Plot1D;
var ind,i:integer;
begin
   ind:=und.addspectre(MyFunctionSpectrum.text,und.spect[indexSpectrum].X,1,1);
   und.spect[ind].typeX:=und.spect[indexSpectrum].typeX;
   for i:=0 to und.spect[indexSpectrum].X-1 do
      begin
         ParserFitFunction.X:=und.spect[indexSpectrum].data[i,0,0].X;
         und.spect[ind].data[i,0,0].X:=und.spect[indexSpectrum].data[i,0,0].X;
         und.spect[ind].data[i,0,0].Z:=ParserFitFunction.value;
      end;
end;

procedure TFittingForm.PlotDiff1D;
var ind,i:integer;
begin
   ind:=und.addspectre(outputSpectrum.text,und.spect[indexSpectrum].X,1,1);
   und.spect[ind].typeX:=und.spect[indexSpectrum].typeX;
   for i:=0 to und.spect[indexSpectrum].X-1 do
      begin
         ParserFitFunction.X:=und.spect[indexSpectrum].data[i,0,0].X;
         und.spect[ind].data[i,0,0].X:=und.spect[indexSpectrum].data[i,0,0].X;
         und.spect[ind].data[i,0,0].Z:=(und.spect[indexSpectrum].data[i,0,0].Z-ParserFitFunction.value);
      end;
end;

procedure TFittingForm.PlotDiff2D;
var ind,i,j:integer;
    aux,ecart:double;
begin
   ind:=deuxd.addspectre(outputSpectrum.text,deuxd.spect[indexSpectrum].X,deuxd.spect[indexSpectrum].Y,1);
   deuxd.spect[ind].typeX:=deuxd.spect[indexSpectrum].typeX;
   deuxd.spect[ind].typeY:=deuxd.spect[indexSpectrum].typeY;
   for i:=0 to deuxd.spect[indexSpectrum].X-1 do
      for j:=0 to deuxd.spect[indexSpectrum].Y-1 do
      begin
         ParserFitFunction.X:=deuxd.spect[indexSpectrum].data[i,j,0].X;
         deuxd.spect[ind].data[i,j,0].X:=deuxd.spect[indexSpectrum].data[i,j,0].X;
         ParserFitFunction.Y:=deuxd.spect[indexSpectrum].data[i,j,0].Y;
         deuxd.spect[ind].data[i,j,0].Y:=deuxd.spect[indexSpectrum].data[i,j,0].Y;
         deuxd.spect[ind].data[i,j,0].Z:=(deuxd.spect[indexSpectrum].data[i,j,0].Z-ParserFitFunction.value);
      end;
end;

procedure TFittingForm.Plot2D;
var ind,i,j:integer;
    aux,ecart:double;
begin
   ind:=deuxd.addspectre(MyFunctionSpectrum.text,deuxd.spect[indexSpectrum].X,deuxd.spect[indexSpectrum].Y,1);
   deuxd.spect[ind].typeX:=deuxd.spect[indexSpectrum].typeX;
   deuxd.spect[ind].typeY:=deuxd.spect[indexSpectrum].typeY;
   for i:=0 to deuxd.spect[indexSpectrum].X-1 do
      for j:=0 to deuxd.spect[indexSpectrum].Y-1 do
      begin
         ParserFitFunction.X:=deuxd.spect[indexSpectrum].data[i,j,0].X;
         deuxd.spect[ind].data[i,j,0].X:=deuxd.spect[indexSpectrum].data[i,j,0].X;
         ParserFitFunction.Y:=deuxd.spect[indexSpectrum].data[i,j,0].Y;
         deuxd.spect[ind].data[i,j,0].Y:=deuxd.spect[indexSpectrum].data[i,j,0].Y;
         deuxd.spect[ind].data[i,j,0].Z:=ParserFitFunction.value;
      end;
end;

procedure TFittingForm.PlotDiff;
begin
   case dimSpectrum of
     1:begin plot1D;plotDiff1D; end;
     2:begin plot2D;plotDiff2D; end;
   end;
end;

{ INPUT : ind, index du spectre 1D a fitter }
function TFittingForm.CalculDifference1D(ind:integer):Tsample;
var i,j:integer;
    aux,ecart:double;
begin
   ecart:=0;
   for i:=0 to und.spect[ind].X-1 do
      begin
         ParserFitFunction.X:=und.spect[ind].data[i,0,0].X;
         aux:=(und.spect[ind].data[i,0,0].Z-ParserFitFunction.value);
         ecart:=ecart+aux*aux;
      end;
   result:=ecart;
end;

{ INPUT : ind, index du spectre 2D a fitter }
function TFittingForm.CalculDifference2D(ind:integer):Tsample;
var i,j:integer;
    aux,ecart:double;
begin
   ecart:=0;
   for i:=0 to deuxd.spect[ind].X-1 do
      for j:=0 to deuxd.spect[ind].Y-1 do
      begin
         ParserFitFunction.X:=deuxd.spect[ind].data[i,j,0].X;
         ParserFitFunction.Y:=deuxd.spect[ind].data[i,j,0].Y;
         aux:=(deuxd.spect[ind].data[i,j,0].Z-ParserFitFunction.value);
         ecart:=ecart+aux*aux;
      end;
   result:=ecart;
end;

{ INPUT : ind, index du spectre 2D a fitter }
function TFittingForm.CalculDifference3D(ind:integer):Tsample;
var i,j,k:integer;
    aux,ecart:double;
begin
   ecart:=0;
   for i:=0 to troisd.spect[ind].X-1 do
      for j:=0 to troisd.spect[ind].Y-1 do
         for k:=0 to troisd.spect[ind].LT-1 do
         begin
         ParserFitFunction.X:=troisd.spect[ind].data[i,j,k].X;
         ParserFitFunction.Y:=troisd.spect[ind].data[i,j,k].Y;
         ParserFitFunction.Z:=troisd.spect[ind].data[i,j,k].LT;
         aux:=(troisd.spect[ind].data[i,j,k].Z-ParserFitFunction.value);
         ecart:=ecart+aux*aux;
         end;
   result:=ecart;
end;

function Banana(Pars,Consts : Array of TSample;PConsts: Array of Pointer): TSample;
begin
   with FittingForm do
   begin
    ParserFitFunction.A:=pars[0];
    ParserFitFunction.B:=pars[1];
    ParserFitFunction.C:=pars[2];
    ParserFitFunction.D:=pars[3];
    ParserFitFunction.E:=pars[4];
    case dimSpectrum of
      1:result:= CalculDifference1D(indexSpectrum);
      2:result:= CalculDifference2D(indexSpectrum);
      3:result:= CalculDifference3D(indexSpectrum);
    end;
{    result:= 100*Sqr(ParserFitFunction.B-Sqr(ParserFitFunction.A))+Sqr(1-ParserFitFunction.A); }
   end;
end;

{ mise a jour de la combo ChoixSpectre à partir de la liste des spectres 2D }
{ INPUT: la liste de spectres "deuxd" en général                            }
{        a passer en var parce que la taille risque de devenir gigantesque  }
procedure TFittingForm.MiseAjourChoixSpectres;
var i:integer;
begin
{reinitialisation des items a chaque ouverture du combo}
   choixSpectre1D.items.clear;choixSpectre2D.items.clear;choixSpectre3D.items.clear;
   for i:=1 to totalspectre do
   begin
      if und.spect[i]<>nil    then choixSpectre1D.items.add(und.spect[i].nom);
      if deuxd.spect[i]<>nil  then choixSpectre2D.items.add(deuxd.spect[i].nom);
      if troisd.spect[i]<>nil then choixSpectre3D.items.add(troisd.spect[i].nom);
   end;
   case ChoixDim.itemIndex of
      0:begin choixSpectre1D.enabled:=true;choixSpectre2D.enabled:=false;choixSpectre3D.enabled:=false; end;
      1:begin choixSpectre1D.enabled:=false;choixSpectre2D.enabled:=true;choixSpectre3D.enabled:=false; end;
      2:begin choixSpectre1D.enabled:=false;choixSpectre2D.enabled:=false;choixSpectre3D.enabled:=true; end;
   end;
end;

procedure TFittingForm.fitterSimplex;
var Iters : integer;
    Fmin:double;
begin
   FittingStatus.caption:='Fitting in progress';FittingStatus.update;
   MiseAJour;
{   Iters := Simplex(Banana,Pars,[],[],FMin,IterMaxBox.value,toleranceBox.value); }
   (* stop if Iters >1000 or Tolerance < 1e-8
   Returns Pars = [1,1] and FMin = 0, meaning x1=1, x2=1 and minimum value is 0 *)
   Abox.value:=pars[0];Bbox.value:=pars[1];Cbox.value:=pars[2];Dbox.value:=pars[3];Ebox.value:=pars[4];
   EcartBox.value:=Fmin;
   FittingStatus.caption:='Idle';FittingStatus.update;
   PlotDiff;
end;




procedure TFittingForm.fitterBFGS;
var Iters : integer;
    Fmin:double;
    IHess:TMtx;
begin
   FittingStatus.caption:='Fitting in progress';FittingStatus.update;
   MiseAJour;
{   Iters := BFGS(Banana,nil,Pars,[],[],FMin,IHess,IterMaxBox.value,toleranceBox.value); { use numeric grad procedure}
   (* stop if Iters >1000 or Tolerance < 1e-8
   Returns Pars = [1,1] and FMin = 0, meaning x1=1, x2=1 and minimum value is 0 *)
   Abox.value:=pars[0];Bbox.value:=pars[1];Cbox.value:=pars[2];Dbox.value:=pars[3];Ebox.value:=pars[4];
   EcartBox.value:=Fmin;
   FittingStatus.caption:='Idle';FittingStatus.update;
   PlotDiff;
end;

procedure TFittingForm.fitterMarquardt;
var Iters : integer;
    Fmin:double;
    IHess:TMtx;
begin
   createIt(IHess);
   FittingStatus.caption:='Fitting in progress';FittingStatus.update;
   MiseAJour;
{   Iters := Marquardt(Banana,nil,Pars,[],[],FMin,IHess,IterMaxBox.value,toleranceBox.value); }
   (* stop if Iters >1000 or Tolerance < 1e-8
   Returns Pars = [1,1] and FMin = 0, meaning x1=1, x2=1 and minimum value is 0 *)
   Abox.value:=pars[0];Bbox.value:=pars[1];Cbox.value:=pars[2];Dbox.value:=pars[3];Ebox.value:=pars[4];
   EcartBox.value:=Fmin;
   FittingStatus.caption:='Idle';FittingStatus.update;
   PlotDiff;
end;

procedure TFittingForm.Simplex1Click(Sender: TObject);
begin
   fitterSimplex;
end;

procedure TFittingForm.BFGS1Click(Sender: TObject);
begin
   fitterBFGS;
end;

procedure TFittingForm.Marquardt1Click(Sender: TObject);
begin
   fitterMarquardt;
end;

procedure TFittingForm.CalculateDifference1Click(Sender: TObject);
var Iters : integer;
    Fmin:double;
begin
  MiseAJour;
  EcartBox.Value:=banana(pars,[],[]);
  PlotDiff;
  { update des différents graphs }
  sendmessage(principalform.handle,wm_redessine,0,0);
end;

procedure TFittingForm.UpdateCombos1Click(Sender: TObject);
begin
   MiseAjourChoixSpectres
end;

procedure TFittingForm.Updateavailablespectra1Click(Sender: TObject);
begin
   MiseAjourChoixSpectres
end;

procedure TFittingForm.Plotdiff2D1Click(Sender: TObject);
begin
   PlotDiff2D
end;

procedure TFittingForm.FormCreate(Sender: TObject);
begin
left:=0;top:=0;width:=523;height:=250;
WindowState:=wsMinimized;
end;

{ on met les combos a jour a chaque fois que la fenetre de fit reprend le focus }
procedure TFittingForm.FormActivate(Sender: TObject);
var i:integer;
begin
{reinitialisation des items a chaque ouverture du combo}
MiseAjourChoixSpectres
{   choixSpectre1D.items.clear;
   for i:=1 to totalspectre do
   begin
      if und<>nil then
      if und.spect[i]<>nil then choixSpectre1D.items.add(und.spect[i].nom);
   end;
   choixSpectre2D.items.clear;
   for i:=1 to totalspectre do
   begin
      if deuxd<>nil then
      if deuxd.spect[i]<>nil then choixSpectre2D.items.add(deuxd.spect[i].nom);
   end;
   choixSpectre3D.items.clear;
   for i:=1 to totalspectre do
   begin
      if troisd<>nil then
      if troisd.spect[i]<>nil then choixSpectre3D.items.add(troisd.spect[i].nom);
   end; }
end;

procedure TFittingForm.ChoixDimClick(Sender: TObject);
begin
   MiseAjourChoixSpectres
end;

end.
