{ implementation du fit en Differential Evolution                               }
{ 1ere version le 15.06.06                                                      }
{ NOTES                                                                         }
{ * pondération : Q^4 * (Rmeas-Rcalc) * ( log (Rmeas) - log (Rcalc) ) pour les multicouches  }
{ ********* A FAIRE *********************************************************** }
{ * DE avec des multicouches, est-ce que ca marche                              }
{ * suivi du E0 pendant le fit DE                                               }

unit DifferentialEvolution;

interface

uses math387,MtxVec,DataReflec;

const
  PopMultiplier = 10; { population number = PopMultiplier * Nbr Fitted Parameters }
  km = 0.7;     { constante de mutation }
  kr = 0.5;     { recombination constant }

  
procedure Minimisation_DifferentialEvolution(var MyEP:TExperimentalParameters;
                                             var MyML:TMultiLayer;
                                             var MyDataExp,MyDataCalc:TDataReflectivity);

implementation

uses forms,calculReflec,SimulReflecForm,sysutils,minimisation, EditData;


{ minimisation utilisant un algorithme de differential evolution }
procedure Minimisation_DifferentialEvolution(var MyEP:TExperimentalParameters;
                                             var MyML:TMultiLayer;
                                             var MyDataExp,MyDataCalc:TDataReflectivity);

var i,j,i0,j0,jidx,n,m,nbLayers,niter,nbrIter,indEF:integer;
    EFmin,EF0:double;
    a1,a2:integer;
    P:array of TMultilayer;
    b,b1,t:TMultilayer; { il faut créer ces vecteurs quelque part probablement }

function CalculDifference(var ML:TMultiLayer):double;
var diff,diffSmooth:double;
begin
   CalculateReflectivity(MyEP,ML,MyDataCalc);
   { application du facteur de normalisation }
   MyDataCalc.Normalise(MyEP.NormalisationFactor);
   { Calcul de la différence aux moindres carrés }
{   diff:=difference2(MyDataExp.R,MyDataExp.Rs,MyDataCalc.R,MyEP.PonderationType); }
   diff:=Difference(MyDataExp,MyDataCalc,MyEP);
   result:=diff;
end;

begin
  StopFittingFlag:=false;
  ReplotFlag:=false;
  with SimulReflecWindow do
  begin FittingStatus.caption:='Fitting...';FittingStatus.update;
        nbrIter:=spinEditNbrIter.value end;
  { création d'un spectre des Energies de minimisation }
  indEF:=und.addspectre('EF',nbrIter,1,1);

  { number of ajustable parameters }
  n:=MyML.nbrParam;
  nbLayers:=MyML.nbrLayers;
  { creation du vecteur de population }
  { on utilise le type MyMultiLayer }
  m:=PopMultiplier*n;
  setlength(P,m);
  for i:=0 to m-1 do P[i]:=TMultiLayer.create;
  b:=TMultiLayer.create;b1:=TMultiLayer.create;t:=TMultiLayer.create;
  { on clone le point le départ }
  for i:=0 to m-1 do P[i].clone(MyML);
  { on tire les paramètres à fitter au hasard dans les vecteurs P1 à Pm }
  for i:=1 to m-1 do
  begin
    { on utilise le step comme paramètre qui définit le range +/- step }
    { scan sur les différentes couches }
    with P[i] do
    begin
      for j:=0 to nbrLayers-2 do
      begin { on se limite au cas des RX pour le moment }
        { ici le range utilisé est % au vecteur P[j], ce n'est pas forcement juste }
        if zYN[j]=0 then z[j]:=z[j]+(random-0.5)*zStep[j];
        if dYN[j]=0 then d[j]:=d[j]+(random-0.5)*dStep[j];
        if fYN[j]=0 then f[j]:=f[j]+(random-0.5)*fStep[j];
        if fiYN[j]=0 then fi[j]:=fi[j]+(random-0.5)*fiStep[j];
        if bYN[j]=0 then bn[j]:=bn[j]+(random-0.5)*bStep[j];
        if biYN[j]=0 then bni[j]:=bni[j]+(random-0.5)*biStep[j];
        if MYN[j]=0 then Aim[j]:=Aim[j]+(random-0.5)*AimStep[j];        
        if sigYN[j]=0 then sigma[j]:=sigma[j]+(random-0.5)*sigmaStep[j];
      end;
    P[i].UpdatePotential(MyEP);
   end; { fin du with }
  end;

{ boucle principale de génération }
for niter:=1 to nbrIter do
begin
  { évaluation de la function d'erreur pour chaque population }
  EFmin:=1e10;i0:=0;
  for i:=0 to m-1 do
  begin EF0:=CalculDifference(P[i]);
        Application.processMessages;
        if StopFittingFlag=true then begin MyML.clone(b);MyML.UpdatePotential(MyEP);exit end;
        if ReplotFlag=true then begin CalculateReflectivity(MyEP,b,MyDataCalc);SimulReflecWindow.PlotReflectivity(MyEP,MyDataCalc);replotFlag:=false end;
        if EF0<EFmin then begin EFmin:=EF0;i0:=i; end;
        end;
  b.clone(P[i0]); { on utilise le meilleur vecteur disponible }

{ boucle sur les populations pour effectuer les mutations }
for i:=0 to m-1 do
begin
  { mutation }
  { choix de 2 populations au hasard }
  a1:=trunc(random(m));
  a2:=a1;
  while a2=a1 do a2:=trunc(random(m));
  { mutation du vecteur b -> b1 }
  b1.clone(b);
  for j:=0 to nbLayers-2 do
    begin { on se limite au cas des RX pour le moment }
    if b1.zYN[j]=0 then b1.z[j]:=b.z[j]+km*(P[a1].z[j]-P[a2].z[j]);
    if b1.dYN[j]=0 then b1.d[j]:=b.d[j]+km*(P[a1].d[j]-P[a2].d[j]);
    if b1.fYN[j]=0 then b1.f[j]:=b.f[j]+km*(P[a1].f[j]-P[a2].f[j]);
    if b1.fiYN[j]=0 then b1.fi[j]:=b.fi[j]+km*(P[a1].fi[j]-P[a2].fi[j]);
    if b1.bYN[j]=0 then b1.bn[j]:=b.bn[j]+km*(P[a1].bn[j]-P[a2].bn[j]);
    if b1.biYN[j]=0 then b1.bni[j]:=b.bni[j]+km*(P[a1].bni[j]-P[a2].bni[j]);
    if b1.MYN[j]=0 then b1.Aim[j]:=b.Aim[j]+km*(P[a1].Aim[j]-P[a2].Aim[j]);
    if b1.sigYN[j]=0 then b1.sigma[j]:=b.sigma[j]+km*(P[a1].sigma[j]-P[a2].sigma[j]);
    end;
  { creation of the trial vector t }
  { choice of a random parameter j0 }
  t.clone(P[i]);
  j0:=trunc(random(n));
  jidx:=0; { index du paramètre fittable en cours }
  for j:=0 to nbLayers-2 do
    begin { on se limite au cas des RX pour le moment }
    if t.zYN[j]=0 then begin
                     if jidx<>j0 then if (random<0.5) then t.z[j]:=b1.z[j] else t.z[j]:=P[i].z[j]
                                 else t.z[j]:=b1.z[j];
                     inc(jidx);end;
    if t.dYN[j]=0 then begin
                     if jidx<>j0 then if (random<0.5) then t.d[j]:=b1.d[j] else t.d[j]:=P[i].d[j]
                                 else t.d[j]:=b1.d[j];
                     inc(jidx);end;
    if t.fYN[j]=0 then begin
                     if jidx<>j0 then if (random<0.5) then t.f[j]:=b1.f[j] else t.f[j]:=P[i].f[j]
                                 else t.f[j]:=b1.f[j];
                     inc(jidx);end;
    if t.fiYN[j]=0 then begin
                     if jidx<>j0 then if (random<0.5) then t.fi[j]:=b1.fi[j] else t.fi[j]:=P[i].fi[j]
                                 else t.fi[j]:=b1.fi[j];
                     inc(jidx);end;
    if t.bYN[j]=0 then begin
                     if jidx<>j0 then if (random<0.5) then t.bn[j]:=b1.bn[j] else t.bn[j]:=P[i].bn[j]
                                 else t.bn[j]:=b1.bn[j];
                     inc(jidx);end;
    if t.biYN[j]=0 then begin
                     if jidx<>j0 then if (random<0.5) then t.bni[j]:=b1.bni[j] else t.bni[j]:=P[i].bni[j]
                                 else t.bni[j]:=b1.bni[j];
                     inc(jidx);end;
    if t.MYN[j]=0 then begin
                     if jidx<>j0 then if (random<0.5) then t.Aim[j]:=b1.Aim[j] else t.Aim[j]:=P[i].Aim[j]
                                 else t.Aim[j]:=b1.Aim[j];
                     inc(jidx);end;
    if t.sigYN[j]=0 then begin
                     if jidx<>j0 then if (random<0.5) then t.z[j]:=b1.sigma[j] else t.sigma[j]:=P[i].sigma[j]
                                 else t.z[j]:=b1.sigma[j];
                     inc(jidx);end;
    end;
  { application des contraintes au cas ou certains paramètres sont allés trop loin }
  for j:=0 to nbLayers-2 do
  begin { on se limite au cas des RX pour le moment }
    { ici le range utilisé est % au vecteur P[j], ce n'est pas forcement juste }
    if t.zYN[j]=0 then if (t.z[j]>(MyML.z[j]+MyML.zStep[j])) or (t.z[j]<(MyML.z[j]-MyML.zStep[j]))
                       then t.z[j]:=MyML.z[j]+(random-0.5)*MyML.zStep[j];
    if t.dYN[j]=0 then if (t.d[j]>(MyML.d[j]+MyML.dStep[j])) or (t.d[j]<(MyML.d[j]-MyML.dStep[j]))
                       then t.d[j]:=MyML.d[j]+(random-0.5)*MyML.dStep[j];
    if t.fYN[j]=0 then if (t.f[j]>(MyML.f[j]+MyML.fStep[j])) or (t.f[j]<(MyML.f[j]-MyML.fStep[j]))
                       then t.f[j]:=MyML.f[j]+(random-0.5)*MyML.fStep[j];
    if t.fiYN[j]=0 then if (t.fi[j]>(MyML.fi[j]+MyML.fiStep[j])) or (t.f[j]<(MyML.fi[j]-MyML.fiStep[j]))
                       then t.fi[j]:=MyML.fi[j]+(random-0.5)*MyML.fiStep[j];
    if t.bYN[j]=0 then if (t.bn[j]>(MyML.bn[j]+MyML.bStep[j])) or (t.bn[j]<(MyML.bn[j]-MyML.bStep[j]))
                       then t.bn[j]:=MyML.bn[j]+(random-0.5)*MyML.bStep[j];
    if t.biYN[j]=0 then if (t.bni[j]>(MyML.bni[j]+MyML.biStep[j])) or (t.bni[j]<(MyML.bni[j]-MyML.biStep[j]))
                       then t.bni[j]:=MyML.bni[j]+(random-0.5)*MyML.biStep[j];
    if t.MYN[j]=0 then if (t.Aim[j]>(MyML.Aim[j]+MyML.AimStep[j])) or (t.Aim[j]<(MyML.Aim[j]-MyML.AimStep[j]))
                       then t.Aim[j]:=MyML.Aim[j]+(random-0.5)*MyML.AimStep[j];
    if t.sigYN[j]=0 then if (t.sigma[j]>(MyML.sigma[j]+MyML.sigmaStep[j])) or (t.sigma[j]<(MyML.sigma[j]-MyML.sigmaStep[j]))
                       then t.sigma[j]:=MyML.sigma[j]+(random-0.5)*MyML.sigmaStep[j];
  end;
  { calcul de l'erreur avec le vecteur t }
  t.UpdatePotential(MyEP);
  EF0:=CalculDifference(t);
  if EF0<EFmin then P[i].clone(t);

  Application.processMessages;
  if StopFittingFlag=true then begin MyML.clone(b);MyML.UpdatePotential(MyEP);exit end;
  if ReplotFlag=true then begin CalculateReflectivity(MyEP,b,MyDataCalc);SimulReflecWindow.PlotReflectivity(MyEP,MyDataCalc);replotFlag:=false end;
end; { fin de la boucle sur les populations }

  { mise à jour de l'affichage }
  SimulReflecWindow.updateDisplayWhileFitting(SimulReflecWindow.MultiLayerGrid,b,MyEP,EFmin);
  with SimulReflecWindow do IterCrt.caption:=IntToStr(nIter);
  und.spect[indEF].data[niter-1,0,0].x:=niter;und.spect[indEF].data[niter-1,0,0].z:=EFmin;

end; { fin de la boucle sur le nombre d'itérations }

{ on sort le meilleur vecteur }
MyML.clone(b);
MyML.UpdatePotential(MyEP);

{ purge des tableaux }
  for i:=0 to m-1 do P[i].free;
  b.free;b1.free;t.free;

end;

end.
