optimierung neu

2 files changed, 427 insertions, 0 deletions

diff --git a/optimierung.inc b/optimierung.inc
new file mode 100644
index 0000000..80a5cb2
--- /dev/null
+++ b/optimierung.inc
@@ -0,0 +1,22 @@
+{$IFDEF tGausz2dFitParameterBuffer.initSamples}
+//procedure tGausz2dFitParameterBuffer.initSamples(qu: tExtendedArray; xSteps,tSiz,zoom: longint);
+var
+  x,y: longint;
+begin
+  sXSteps:=xSteps div zoom;
+  sTSiz:=tSiz div zoom;
+  setLength(samples,sXSteps*sTSiz);
+  for x:=0 to sXSteps-1 do
+    for y:=0 to sTSiz-1 do
+      samples[x+y*sXSteps]:=0;
+  for x:=0 to xSteps-1 do
+    for y:=0 to tSiz-1 do
+      samples[round(x/(xSteps-1)*(sXSteps-1)) + round(y/(tSiz-1)*(sTSiz-1))*sTSiz]:=
+        samples[round(x/(xSteps-1)*(sXSteps-1)) + round(y/(tSiz-1)*(sTSiz-1))*sTSiz] +
+        qu[x + y*xSteps];
+  for x:=0 to sXSteps-1 do
+    for y:=0 to sTSiz-1 do
+      samples[x+y*sXSteps]:=
+        samples[x+y*sXSteps] / xSteps * sXSteps / tSiz * sTSiz;
+end;
+{$ENDIF}
diff --git a/optimierung.pas b/optimierung.pas
new file mode 100644
index 0000000..5c8e91e
--- /dev/null
+++ b/optimierung.pas
@@ -0,0 +1,405 @@
+unit optimierung;
+
+{$mode objfpc}{$H+}
+
+interface
+
+uses
+  Classes, SysUtils, lowlevelunit;
+
+type
+  tParameterBuffer = class(tObject)
+    function abstandsQuadrat(const params: tExtendedArray): extended; virtual; abstract;
+    function aktiviereWeitereParameter(const params: tExtendedArray; out neu,neuD: tExtendedArray): boolean; virtual; abstract;
+  end;
+
+  tOptimierer = class(tObject)
+  private
+    buf: tParameterBuffer;
+    dim: longint;
+    function rPar(i: longint): extended; virtual; abstract;
+    function rAbst: extended; virtual; abstract;
+  public
+    property parameter[i: longint]: extended
+      read rPar;
+    property abstand: extended
+      read rAbst;
+    constructor create(buffer: tParameterBuffer; dimension: longint);
+    function optimiere(minVerbesserung, minSchritt: extended; maxSchritte: int64): int64; virtual; abstract;
+  end;
+
+  tDownHillSimplexOptimierer = class(tOptimierer)
+  private
+    pars:                 array of tExtendedArray;
+    absts:                tExtendedArray;
+    bester,schlechtester: longint;
+    function rPar(i: longint): extended; override;
+    function rAbst: extended; override;
+    function simplexSchritt(wer: longint; wieWeit: extended): extended; inline;
+    procedure ermittleRaenge(schlechtesten, besten: boolean); inline;
+  public
+    constructor create(buffer: tParameterBuffer; startParams,dStartParams: tExtendedArray);
+    destructor destroy; override;
+    function optimiere(minVerbesserung, minSchritt: extended; maxSchritte: int64): int64; override;
+  end;
+
+  tGausz2dFitParameterBuffer = class(tParameterBuffer)
+  private
+    anz,iteration: longint;
+    xxs,yys,xys:   tExtendedArray;
+    sXSteps,sTSiz: longint;
+    samples:       tExtendedArray;
+    procedure setzeBufferParameter(const params: tExtendedArray); inline;
+  public
+    // offset, { xx, yy, x0, y0, a, e }
+    // exp ( -(xx^2 * cos(a)^2 + yy^2 * sin(a)^2) * (X-x0)^2 - (yy^2 * cos(a)^2 + xx^2 * sin(a)^2) * (Y-y0)^2 + 2(xx^2 - yy^2) * cos(a)sin(a) * (X-x0) * (Y-y0) + e)
+    geschaetzteFitParams,
+      geschaetzteDFitParams: tExtendedArray;
+    constructor create(anzahl: longint);
+    destructor destroy; override;
+    procedure initSamples(qu: array of single; xSteps,tSiz,zoom: longint); overload;
+    procedure initSamples(qu: array of double; xSteps,tSiz,zoom: longint); overload;
+    procedure initSamples(qu: array of extended; xSteps,tSiz,zoom: longint); overload;
+    procedure initParamsSchaetzung;
+    procedure initOffset(offset: extended);
+    function aktiviereWeitereParameter(const params: tExtendedArray; out neu,neuD: tExtendedArray): boolean; override;
+    function abstandsQuadrat(const params: tExtendedArray): extended; override;
+  end;
+
+
+// TODOs:
+// - Grenzen für Parameter
+// - einstellbare Schrittlängenfaktoren
+
+implementation
+
+uses
+  math;
+
+// tOptimierer *****************************************************************
+
+constructor tOptimierer.create(buffer: tParameterBuffer; dimension: longint);
+begin
+  inherited create;
+  buf:=buffer;
+  dim:=dimension;
+end;
+
+// tDownHillSimplexOptimierer **************************************************
+
+function tDownHillSimplexOptimierer.rPar(i: longint): extended;
+begin
+  result:=pars[bester][i];
+end;
+
+function tDownHillSimplexOptimierer.rAbst: extended;
+begin
+  result:=absts[bester];
+end;
+
+function tDownHillSimplexOptimierer.simplexSchritt(wer: longint; wieWeit: extended): extended;
+var
+  i,j: longint;
+  tmp: extended;
+begin           // neu - Schwerpunkt = wieWeit * (Schwerpunkt - alt)   =>  wieWeit = 1 <=> Spiegelung
+  result:=0;
+  for i:=0 to dim-1 do begin
+    tmp:=0;
+    for j:=0 to dim do
+      if j<>wer then
+        tmp:=tmp+pars[j,i];
+    result:=result + sqr((1 + wieWeit) * (tmp/dim - pars[wer,i]));
+    pars[wer,i]:=tmp * (1 + wieWeit)/dim - wieWeit * pars[wer,i];
+  end;
+end;
+
+procedure tDownHillSimplexOptimierer.ermittleRaenge(schlechtesten, besten: boolean);
+var
+  i: longint;
+begin
+  if besten then
+    for i:=0 to dim do
+      if absts[i] < absts[bester] then
+        bester:=i;
+  if schlechtesten then
+    for i:=0 to dim do
+      if absts[i] > absts[schlechtester] then
+        schlechtester:=i;
+end;
+
+constructor tDownHillSimplexOptimierer.create(buffer: tParameterBuffer; startParams,dStartParams: tExtendedArray);
+var
+  i,j: longint;
+begin
+  inherited create(buffer, length(startParams));
+  if length(dStartParams)<>dim then
+    fehler('tDownHillSimplexOptimierer.create: die Längen von startParams und dStartParams müssen übereinstimmen!');
+  setLength(pars,dim+1);
+  for i:=0 to dim do begin
+    setLength(pars[i],dim);
+    for j:=0 to dim-1 do
+      pars[i,j]:=startParams[j] + dStartParams[j] * byte(j=i);
+  end;
+  setLength(absts,dim+1);
+  bester:=0;
+  schlechtester:=0;
+  for i:=0 to length(absts)-1 do
+    absts[i]:=buf.abstandsQuadrat(pars[i]);
+  ermittleRaenge(true,true);
+end;
+
+destructor tDownHillSimplexOptimierer.destroy;
+var
+  i: longint;
+begin
+  for i:=0 to length(pars)-1 do
+    setLength(pars[i],0);
+  setLength(pars,0);
+  setLength(absts,0);
+  inherited destroy;
+end;
+
+function tDownHillSimplexOptimierer.optimiere(minVerbesserung, minSchritt: extended; maxSchritte: int64): int64;
+var
+  schritt,verbesserung,tmp: extended;
+  i,j,altDim:               longint;
+  neu,neuD:                 tExtendedArray;
+begin
+  if (minVerbesserung<=0) and (minSchritt<=0) and (maxSchritte<=0) then
+    fehler('tDownHillSimplexOptimierer.optimiere benötigt ein Abbruchkriterium!');
+
+  if minSchritt>0 then
+    minSchritt:=sqr(minSchritt);
+
+  result:=0;
+  verbesserung:=0;
+  repeat
+    repeat
+      inc(result);
+      schritt:=simplexSchritt(schlechtester,1);
+
+      tmp:=absts[schlechtester];
+      absts[schlechtester]:=buf.abstandsQuadrat(pars[schlechtester]);
+      verbesserung:=tmp - absts[schlechtester];
+      if absts[schlechtester] < absts[bester] then begin // ist jetzt der beste
+        // => wir gehen noch weiter in dieser Richtung
+        simplexSchritt(schlechtester,-3);
+        schritt:=schritt * sqr(2); // statt von -1 zu 1 sind wir von -1 zu 3 gelaufen => Faktor 2
+        absts[schlechtester]:=buf.abstandsQuadrat(pars[schlechtester]);
+        verbesserung:=tmp - absts[schlechtester];
+        ermittleRaenge(true,true);
+        continue;
+      end;
+      if verbesserung < 0 then begin // verschlechtert
+        // => wir gehen einen kleineren Schritt
+        simplexSchritt(schlechtester,-0.5);
+        schritt:=schritt * sqr(3/4); // statt von -1 zu 1 sind wir von -1 zu 0.5 gelaufen => Faktor 3/4
+        absts[schlechtester]:=buf.abstandsQuadrat(pars[schlechtester]);
+        verbesserung:=tmp - absts[schlechtester];
+        ermittleRaenge(true,false);
+        continue;
+      end;
+      ermittleRaenge(true,false);
+    until
+      ((maxSchritte>0) and (result>=maxSchritte)) or
+      ((minVerbesserung>0) and (verbesserung<minVerbesserung)) or
+      ((minSchritt>0) and (schritt<minSchritt));
+    if buf.aktiviereWeitereParameter(pars[bester],neu,neuD) then begin
+      altDim:=dim;
+      dim:=dim+length(neu);
+      for i:=0 to altDim do begin // alte Eckpunkte um neue Dimensionen erweitern
+        setLength(pars[i],dim);
+        for j:=altDim to dim-1 do
+          pars[i,j]:=neu[j-altDim];
+      end;
+      setLength(pars,dim+1);
+      for i:=altDim+1 to dim do begin // neue Eckpunkte an besten Punkt anfügen
+        setLength(pars[i],dim);
+        for j:=0 to dim-1 do
+          pars[i,j]:=pars[bester,j];
+        for j:=altDim to dim-1 do
+          pars[i,j]:=pars[i,j] + neuD[j-altDim] * byte(j+1=i);
+      end;
+      setLength(absts,dim+1);
+      for i:=0 to dim do // auch alte Eckpunkte müssen aktualisiert werden, da sich die Fitfunktion geändert hat!
+        absts[i]:=buf.abstandsQuadrat(pars[i]);
+      ermittleRaenge(true,true);
+    end
+    else
+      break;
+  until false;
+end;
+
+// tGausz2dFitParameterBuffer **************************************************
+
+const
+  xx_i  = 1;
+  yy_i  = 2;
+  x0_i  = 3;
+  y0_i  = 4;
+  a_i   = 5;
+  e_i   = 6;
+  p_l   = 6;
+
+procedure tGausz2dFitParameterBuffer.setzeBufferParameter(const params: tExtendedArray);
+var
+  i:     longint;
+  ca,sa: extended;
+begin
+  for i:=0 to iteration-1 do begin
+    ca:=cos(params[p_l*i + a_i]);
+    sa:=sin(params[p_l*i + a_i]);
+    xxs[i]:=-sqr(params[p_l*i + xx_i]*ca)-sqr(params[p_l*i + yy_i]*sa);
+    yys[i]:=-sqr(params[p_l*i + yy_i]*ca)-sqr(params[p_l*i + xx_i]*sa);
+    xys[i]:=(sqr(params[p_l*i + xx_i])-sqr(params[p_l*i + yy_i]))*2*ca*sa;
+  end;
+end;
+
+constructor tGausz2dFitParameterBuffer.create(anzahl: longint);
+begin
+  inherited create;
+  setLength(geschaetzteFitParams,0);
+  setLength(geschaetzteDFitParams,0);
+  setLength(xxs,0);
+  setLength(yys,0);
+  setLength(xys,0);
+  setLength(samples,0);
+  anz:=anzahl;
+end;
+
+destructor tGausz2dFitParameterBuffer.destroy;
+begin
+  setLength(geschaetzteFitParams,0);
+  setLength(geschaetzteDFitParams,0);
+  setLength(xxs,0);
+  setLength(yys,0);
+  setLength(xys,0);
+  setLength(samples,0);
+  inherited destroy;
+end;
+
+{$DEFINE tGausz2dFitParameterBuffer.initSamples}
+procedure tGausz2dFitParameterBuffer.initSamples(qu: array of single; xSteps,tSiz,zoom: longint);
+{$I optimierung.inc}
+procedure tGausz2dFitParameterBuffer.initSamples(qu: array of double; xSteps,tSiz,zoom: longint);
+{$I optimierung.inc}
+procedure tGausz2dFitParameterBuffer.initSamples(qu: array of extended; xSteps,tSiz,zoom: longint);
+{$I optimierung.inc}
+{$UNDEF tGausz2dFitParameterBuffer.initSamples}
+
+procedure tGausz2dFitParameterBuffer.initParamsSchaetzung;
+begin
+  setLength(geschaetzteFitParams,0);
+  setLength(geschaetzteDFitParams,0);
+  iteration:=0;
+end;
+
+procedure tGausz2dFitParameterBuffer.initOffset(offset: extended);
+begin
+  setLength(geschaetzteFitParams,1);
+  setLength(geschaetzteDFitParams,1);
+  geschaetzteFitParams[0]:=offset;
+  geschaetzteDFitParams[0]:=1;
+end;
+
+function tGausz2dFitParameterBuffer.aktiviereWeitereParameter(const params: tExtendedArray; out neu,neuD: tExtendedArray): boolean;
+var
+  i,x,y:     longint;
+  dx,dy,tmp: extended;
+  samples2:  tExtendedArray;
+begin
+  result:= iteration < anz;
+  if not result then
+    exit;
+
+  setzeBufferParameter(params);
+
+  inc(iteration);
+  setLength(xxs,length(xxs)+1);
+  setLength(yys,length(yys)+1);
+  setLength(xys,length(xys)+1);
+  setLength(neu,6);
+  setLength(neuD,6);
+  for i:=0 to 3 do // xx,yy,x0,y0
+    neuD[i]:=1;
+  neuD[a_i-1]:=0.1; // a
+  neuD[e_i-1]:=0.1; // e
+
+  setLength(samples2,length(samples));
+
+  for y:=0 to sTSiz-1 do
+    for x:=0 to sXSteps-1 do begin
+      samples2[x+y*sXSteps]:=samples[x+y*sXSteps] - params[0];
+      for i:=0 to iteration-2 do begin
+        dx:= x - params[p_l*i + x0_i];
+        dy:= y - params[p_l*i + y0_i];
+        samples2[x+y*sXSteps]:= samples2[x+y*sXSteps] -
+          exp(
+            params[p_l*i + e_i] +
+            xxs[i] * sqr(dx) +
+            yys[i] * sqr(dy) +
+            xys[i] * dx * dy
+          );
+      end;
+    end;
+  neu[e_i-1]:=0;
+  neu[xx_i-1]:=(10 + 0.5 * random) / sXSteps;
+  neu[yy_i-1]:=(10 + 0.5 * random) / sTSiz;
+  neu[x0_i-1]:=0;
+  neu[y0_i-1]:=0;
+  neu[a_i-1]:=0.01 * random - 0.005;
+  tmp:=-1;
+  for y:=0 to sTSiz-1 do
+    for x:=0 to sXSteps-1 do
+      if samples2[x+y*sXSteps]>tmp then begin
+        tmp:= samples2[x+y*sXSteps];
+        neu[x0_i-1]:=x;
+        neu[y0_i-1]:=y;
+      end;
+  if tmp>0 then
+    neu[e_i-1]:=ln(tmp)
+  else
+    neu[e_i-1]:=0;
+  x:=round(neu[x0_i-1]);
+  y:=round(neu[y0_i-1]);
+  for i:=1 to min(x,sXSteps-1-x) do
+    if samples2[x+i+y*sXSteps] + samples2[x-i+y*sXSteps] < 2 * exp(-1) * samples2[x+y*sXSteps] then begin
+      neu[xx_i-1]:=1/i;
+      break;
+    end;
+  for i:=1 to min(y,sTSiz-1-y) do
+    if samples2[x+(y+i)*sXSteps] + samples2[x+(y-i)*sXSteps] - 2 * params[0] < 2 * exp(-1) * samples2[x+y*sXSteps] then begin
+      neu[yy_i-1]:=1/i;
+      break;
+    end;
+
+  setLength(samples2,0);
+end;
+
+function tGausz2dFitParameterBuffer.abstandsQuadrat(const params: tExtendedArray): extended;
+var
+  i,x,y:     longint;
+  tmp,dx,dy: extended;
+begin
+  setzeBufferParameter(params);
+
+  result:=0;
+  for x:=0 to sXSteps-1 do
+    for y:=0 to sTSiz-1 do begin
+      tmp:= params[0] - samples[x + y*sXSteps];
+      for i:=0 to iteration-1 do begin
+        dx:= x - params[p_l*i + x0_i];
+        dy:= y - params[p_l*i + y0_i];
+        tmp:= tmp +
+          exp(
+            params[p_l*i + e_i] +
+            xxs[i] * sqr(dx) +
+            yys[i] * sqr(dy) +
+            xys[i] * dx * dy
+          );
+      end;
+      result:= result + sqr(tmp);
+    end;
+end;
+
+end.