unit myMath;

interface

uses math, sysutils;

type TExtPoint = record
                  case byte of
                   0: (x,y: extended;);     // coordinates or
                   1: (val,err: extended;); // value and relative error
                 end;

function myfloattostr(f: double): string;

function multiply(a,b: TExtPoint): TExtPoint;
function divide(a,b: TExtPoint): TExtPoint;
function valuetostr(a: TExtPoint; separator: String): String;
function valuetostr(a: TExtPoint): String;
function valuetostr(a: TExtPoint; separator: String; mitFehler: boolean): String;
function valuetostr(a: TExtPoint; mitFehler: boolean): String;

procedure GaussFitDownhillSimplex(W: array of TExtPoint; var y0,dt,delta: extended);
procedure GaussFitGradientenverfahren(W: array of TExtPoint; var y0,dt,delta: extended);

var Komma: char = ',';

implementation

function myfloattostr(f: double): string;
var s,t: string;
begin
 s:=floattostr(f);
 t:='';
 while pos('.',s)>0 do
  begin
   t:=t+copy(s,1,pos('.',s)-1)+komma;
   delete(s,1,pos('.',s));
  end;
 myfloattostr:=t+s;
end;

//*****************************************************************************

function multiply(a,b: TExtPoint): TExtPoint;
var erg: TExtPoint;
begin
 erg.val:=a.val*b.val;
 erg.err:=a.err+b.err;
 multiply:=erg;
end;

function divide(a,b: TExtPoint): TExtPoint;
var erg: TExtPoint;
begin
 if b.val=0 then writeln('Division durch Null!');
 erg.val:=a.val/b.val;
 erg.err:=a.err+b.err;
 divide:=erg;
end;

function valuetostr(a: TExtPoint; separator: String): String;
begin
 valuetostr:=valuetostr(a,separator,true);
end;

function valuetostr(a: TExtPoint): String;
begin
 valuetostr:=valuetostr(a,#9);
end;

function valuetostr(a: TExtPoint; separator: String; mitFehler: boolean): String;
begin
 if mitFehler then valuetostr:=myfloattostr(a.val)+separator+myfloattostr(a.val*a.err)+separator+myfloattostr(a.err)
              else valuetostr:=myfloattostr(a.val);
end;

function valuetostr(a: TExtPoint; mitFehler: boolean): String;
begin
 valuetostr:=valuetostr(a,#9,mitFehler);
end;

//*****************************************************************************

function qerr(W: array of TExtPoint; Param: TExtPoint): extended;
var I:   Integer;
    erg: extended;
begin
 erg:=0;
 For I:=0 to length(W)-1 do
  erg:=erg +
       sqr(W[I].y - Param.y*exp(-sqr(W[I].x/Param.x)));
 qerr:=erg;
end;

procedure GaussFitDownhillSimplex(W: array of TExtPoint; var y0,dt,delta: extended);
var Simplex:    array[0..2] of TExtPoint;
    Errs:       array[0..2] of Extended;
    mitt,np:    TExtPoint;
    ma,mi,nerr: extended;
    I,J:        Integer;
const gutFaktor = 2.5;
      schlechtFaktor = 0.5;
      mittelFaktor = 1.0;
begin
 Simplex[0].x:=dt;
 Simplex[0].y:=y0;
 Simplex[1].x:=dt+delta*1000;
 Simplex[1].y:=y0;
 Simplex[2].x:=dt;
 Simplex[2].y:=y0+delta*1000;
 For I:=0 to 2 do
  Errs[I]:=qerr(W,Simplex[I]);
 repeat
  J:=0;
  For I:=1 to 2 do
   if Errs[I]>Errs[J] then
    J:=I;
  mitt.x:=0;
  mitt.y:=0;
  Ma:=Errs[Byte(J=0)];
  Mi:=Errs[Byte(J=0)];
  For I:=0 to 2 do
   if J<>I then
    begin
     Ma:=Max(Ma,Errs[I]);
     Mi:=Min(Mi,Errs[I]);
     mitt.x:=mitt.x+Simplex[I].x;
     mitt.y:=mitt.y+Simplex[I].y;
    end;
  mitt.x:=mitt.x/2;
  mitt.y:=mitt.y/2;
  np.x:=2*mitt.x-Simplex[J].x;
  np.y:=2*mitt.y-Simplex[J].y;
  nerr:=qerr(W,np);
  np.x:=mitt.x-Simplex[J].x;
  np.y:=mitt.y-Simplex[J].y;
  if nerr<Mi then // Minimum !
   begin
    Simplex[J].x:=mitt.x+gutFaktor*np.x;
    Simplex[J].y:=mitt.y+gutFaktor*np.y;
   end
    else
   if nerr>Ma then // Maximum !
    begin
     Simplex[J].x:=mitt.x+schlechtFaktor*np.x;
     Simplex[J].y:=mitt.y+schlechtFaktor*np.y;
    end
     else
    begin
     Simplex[J].x:=mitt.x+mittelFaktor*np.x;
     Simplex[J].y:=mitt.y+mittelFaktor*np.y;
    end;
  Errs[J]:=qerr(W,Simplex[J]);  
 until (abs((Simplex[0].x-Simplex[1].x)*(Simplex[0].y-Simplex[2].y) -
            (Simplex[0].x-Simplex[2].x)*(Simplex[0].y-Simplex[1].y)) <= delta);
 J:=0;
 For I:=1 to 2 do
  if Errs[I]<Errs[J] then
   J:=I;
 delta:=Errs[J];
 dt:=Simplex[J].x;
 y0:=Simplex[J].y;
end;

procedure GaussFitGradientenverfahren(W: array of TExtPoint; var y0,dt,delta: extended);
var ex,tmp,err: extended;
    i,vz,vzz:   integer;
    dy0,
    ddt,fac:    extended;
    weiter:     boolean;
begin
 dt:=1/dt;
 err:=0;
 fac:=2;
 vz:=0;
 vzz:=0;
 for i:=0 to length(W)-1 do
  err:=err+sqr(W[i].y-y0*exp(-sqr(W[i].x*dt)));
 repeat
  dy0:=0;
  ddt:=0;
  for i:=0 to length(W)-1 do
   begin
    ex:=exp(-sqr(W[i].x*dt));
    tmp:=W[i].y*ex-y0*sqr(ex);
    dy0:=dy0 - tmp;
    ddt:=ddt + sqr(W[i].x)*tmp;
   end;
  ddt:=ddt*dt*y0*2;
  weiter:=sqr(ddt)+sqr(dy0)>1;
  if weiter then tmp:=1/sqrt(sqr(ddt)+sqr(dy0))
            else tmp:=1;
  if (vz<>2*byte(ddt>0)+byte(dy0>0)) or
     (ddt=0) or (dy0=0) then
   begin
    fac:=fac/2;
    vzz:=0;
   end
    else
   begin
    inc(vzz);
    for i:=5 to vzz do
     fac:=min(1,fac*1.2);
   end;
  vz:=2*byte(ddt>0)+byte(dy0>0);
  tmp:=tmp*fac;
  dt:=dt-ddt*tmp;
  y0:=y0-dy0*tmp;
  tmp:=err;
  err:=0;
  for i:=0 to length(W)-1 do
   err:=err+sqr(W[i].y-y0*exp(-sqr(W[i].x*dt)));
 until (tmp-err<delta*err) and (fac<delta);
 dt:=1/dt;
 delta:=err;
end;

end.