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;
    function schrittZuGrosz(const params,schritt: 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;
    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);
    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;
    function schrittZuGrosz(const params,schritt: tExtendedArray): boolean; override;
  end;

// TODOs:
// - 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;

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
  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];
    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;

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
  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!');

  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);
      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
        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,schrittLaenge);
        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 (schrittLaenge<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);
      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 **************************************************

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:=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];
      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;

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.