From d43cd9e54801d0dd23a0bca4c4299795db539e80 Mon Sep 17 00:00:00 2001
From: Erich Eckner <git@eckner.net>
Date: Wed, 19 Dec 2018 13:59:42 +0100
Subject: optimierung.pas: Grenzen für Parameter umgesetzt
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 optimierung.pas | 80 ++++++++++++++++++++++++++++++++++++++++++++-------------
 1 file changed, 62 insertions(+), 18 deletions(-)

diff --git a/optimierung.pas b/optimierung.pas
index 5c8e91e..dd3ef2b 100644
--- a/optimierung.pas
+++ b/optimierung.pas
@@ -11,6 +11,7 @@ type
   tParameterBuffer = class(tObject)
     function abstandsQuadrat(const params: tExtendedArray): extended; virtual; abstract;
     function aktiviereWeitereParameter(const params: tExtendedArray; out neu,neuD: tExtendedArray): boolean; virtual; abstract;
+    function schrittZuGrosz(const params,schritt: tExtendedArray): boolean; virtual; abstract;
   end;
 
   tOptimierer = class(tObject)
@@ -35,7 +36,7 @@ type
     bester,schlechtester: longint;
     function rPar(i: longint): extended; override;
     function rAbst: extended; override;
-    function simplexSchritt(wer: longint; wieWeit: extended): extended; inline;
+    procedure simplexSchritt(wer: longint; wieWeit: extended; var schritt: tExtendedArray; var schrittLaenge: extended); inline;
     procedure ermittleRaenge(schlechtesten, besten: boolean); inline;
   public
     constructor create(buffer: tParameterBuffer; startParams,dStartParams: tExtendedArray);
@@ -64,11 +65,10 @@ type
     procedure initOffset(offset: extended);
     function aktiviereWeitereParameter(const params: tExtendedArray; out neu,neuD: tExtendedArray): boolean; override;
     function abstandsQuadrat(const params: tExtendedArray): extended; override;
+    function schrittZuGrosz(const params,schritt: tExtendedArray): boolean; override;
   end;
 
-
 // TODOs:
-// - Grenzen für Parameter
 // - einstellbare Schrittlängenfaktoren
 
 implementation
@@ -97,18 +97,19 @@ begin
   result:=absts[bester];
 end;
 
-function tDownHillSimplexOptimierer.simplexSchritt(wer: longint; wieWeit: extended): extended;
+procedure tDownHillSimplexOptimierer.simplexSchritt(wer: longint; wieWeit: extended; var schritt: tExtendedArray; var schrittLaenge: extended);
 var
   i,j: longint;
   tmp: extended;
 begin           // neu - Schwerpunkt = wieWeit * (Schwerpunkt - alt)   =>  wieWeit = 1 <=> Spiegelung
-  result:=0;
+  schrittLaenge:=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]));
+    schritt[i]:= (1 + wieWeit) * (tmp/dim - pars[wer,i]);
+    schrittLaenge:=schrittLaenge + sqr(schritt[i]);
     pars[wer,i]:=tmp * (1 + wieWeit)/dim - wieWeit * pars[wer,i];
   end;
 end;
@@ -161,9 +162,9 @@ end;
 
 function tDownHillSimplexOptimierer.optimiere(minVerbesserung, minSchritt: extended; maxSchritte: int64): int64;
 var
-  schritt,verbesserung,tmp: extended;
-  i,j,altDim:               longint;
-  neu,neuD:                 tExtendedArray;
+  schrittLaenge,verbesserung,tmp: extended;
+  i,j,altDim:                     longint;
+  neu,neuD,schritt:               tExtendedArray;
 begin
   if (minVerbesserung<=0) and (minSchritt<=0) and (maxSchritte<=0) then
     fehler('tDownHillSimplexOptimierer.optimiere benötigt ein Abbruchkriterium!');
@@ -171,29 +172,33 @@ begin
   if minSchritt>0 then
     minSchritt:=sqr(minSchritt);
 
+  setLength(schritt,dim);
+  for i:=0 to dim-1 do
+    schritt[i]:=0;
+  schrittLaenge:=0;
   result:=0;
   verbesserung:=0;
   repeat
     repeat
       inc(result);
-      schritt:=simplexSchritt(schlechtester,1);
+      simplexSchritt(schlechtester,1,schritt,schrittLaenge);
 
       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];
+        if not buf.schrittZuGrosz(pars[schlechtester],schritt) then begin
+          // => wir gehen noch etwas weiter in dieser Richtung
+          simplexSchritt(schlechtester,-1.1,schritt,schrittLaenge);
+          absts[schlechtester]:=buf.abstandsQuadrat(pars[schlechtester]);
+          verbesserung:=tmp - absts[schlechtester];
+        end;
         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
+        simplexSchritt(schlechtester,-0.5,schritt,schrittLaenge);
         absts[schlechtester]:=buf.abstandsQuadrat(pars[schlechtester]);
         verbesserung:=tmp - absts[schlechtester];
         ermittleRaenge(true,false);
@@ -203,7 +208,7 @@ begin
     until
       ((maxSchritte>0) and (result>=maxSchritte)) or
       ((minVerbesserung>0) and (verbesserung<minVerbesserung)) or
-      ((minSchritt>0) and (schritt<minSchritt));
+      ((minSchritt>0) and (schrittLaenge<minSchritt));
     if buf.aktiviereWeitereParameter(pars[bester],neu,neuD) then begin
       altDim:=dim;
       dim:=dim+length(neu);
@@ -224,10 +229,19 @@ begin
       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);
+      setLength(schritt,dim);
+      for i:=0 to dim-1 do
+        schritt[i]:=0;
     end
     else
       break;
   until false;
+  tmp:=pars[0,3]; // x0
+  schrittLaenge:=tmp;
+  for i:=1 to dim do begin
+    tmp:=max(tmp,pars[i,3]);
+    schrittLaenge:=min(schrittLaenge,pars[i,3]);
+  end;
 end;
 
 // tGausz2dFitParameterBuffer **************************************************
@@ -384,6 +398,20 @@ begin
   setzeBufferParameter(params);
 
   result:=0;
+  for x:=1 to iteration-1 do
+    for y:=0 to x-1 do begin
+      dx:= params[p_l*x + x0_i] - params[p_l*y + x0_i];
+      dy:= params[p_l*x + y0_i] - params[p_l*y + y0_i];
+      result:=result +
+        exp(
+          params[p_l*x + e_i] +
+          params[p_l*y + e_i] +
+          (xxs[x] + xxs[y]) * sqr(dx) +
+          (yys[x] + yys[y]) * sqr(dy) +
+          (xys[x] + xys[y]) * dx * dy
+        );
+    end;
+  result:= result * sXSteps * sTSiz * 1000;
   for x:=0 to sXSteps-1 do
     for y:=0 to sTSiz-1 do begin
       tmp:= params[0] - samples[x + y*sXSteps];
@@ -402,4 +430,20 @@ begin
     end;
 end;
 
+function tGausz2dFitParameterBuffer.schrittZuGrosz(const params,schritt: tExtendedArray): boolean;
+var
+  i: longint;
+begin
+  result:=true;
+  for i:=0 to iteration-1 do
+    if (abs(schritt[p_l*i + xx_i]*2) > abs(params[p_l*i + xx_i])) or
+      (abs(schritt[p_l*i + yy_i]*2) > abs(params[p_l*i + yy_i])) or
+      (abs(schritt[p_l*i + x0_i]*2) > sXSteps) or
+      (abs(schritt[p_l*i + y0_i]*2) > sTSiz) or
+      (abs(schritt[p_l*i + a_i]) > pi/10) or
+      (abs(schritt[p_l*i + e_i]) > 1) then
+      exit;
+  result:=false;
+end;
+
 end.
-- 
cgit v1.2.3-54-g00ecf