{ Обобщенный поиск минимума функции многих переменных.}

{---------------------------------------------------------------------------
The control units for mass-spectrometer MI1201-AGM
(c) Copyright	 Aleksandrov O.E., 1998

Модуль управления масс-спектрометром МИ1201-АГМ
(c) Собственность     Александрова О.Е., 1998

   Molecular Physics department	 	620002, Екатеринбург, К-2
   USTU, Ekaterinsburg, K-2, 620002     УГТУ,
   RUSSIA     	   			Кафедра молекулярной физики
   phone 75-48-39                       тел. 75-48-39

   E-mail: aleks@dpt.ustu.ru
----------------------------------------------------------------------------}

{Основная поцедура: GenFitX или GenFit. }

unit uGenFit;

INTERFACE

const
  cFPMin:Double=1.0e-300; { минимальное значение для арифметики
			  плавающей точки (1.0e-30; - for Real) }
  cItMax:word=3000;     { максимальное количество итераций }

  cCheckCnt:word=100; {}


const
{$IfDef Seg16}
  cMaxDataSize=($FFFF-4*SizeOf(Word)) div SizeOf(Double);
{$Else}
  cMaxDataSize=(($FFFFFFFF div 2)-4*SizeOf(Word)) div SizeOf(Double);
{$EndIf Def Seg16}

  cMaxParameterSize=32;

type
  tParameterNumber=1..cMaxParameterSize;
  tParametersArray=array[tParameterNumber] of Double;
  tParametersSet=set of tParameterNumber;

  tDataNumber=1..cMaxDataSize;
  tDataArray=array[tDataNumber] of Double;

  tParameters=record
    Size:tParameterNumber;
    Par:tParametersArray;
  end;
  tPtrParameters=^tParameters;

  tData=record
    Size:tDataNumber;
    Data:tDataArray;
  end;
  tPtrData=^tData;

  tErrorCode=(ecOK,
	      ecMaxIterationNumberExeeded,
	      ecNotCovering);

{ Процедура вычисления значения аппроксимирующей функции должна
иметь такой заголовок. P - значения параметров и X - значение аргумента,
при которых вычисляется функция. }
  tFunc  = function(const P:tParameters; X:Double):Double;

{ Процедура вычисления значений всех частных производных (градиента)
аппроксимирующей функции должна иметь такой заголовок.
Значения производных возвращаются в переменной GradF.
P - значения параметров и X - значение аргумента, при которых вычисляется
градиент. }
  tGradF = procedure(const P:tParameters; X:Double; var GradF:tParameters);

{ аргументы см. GeneralFitXX }
procedure GeneralFit(aF:tFunc; aGradF:tGradF;
		     var P:tParameters; const X,Y:tData;
			 Eps:Double;    var ErrCode:tErrorCode);

{ аргументы см. GeneralFitXX }
procedure GeneralFitX(aF:tFunc; aGradF:tGradF;
		      var P:tParameters; const X,Y:tData;
			  Eps:Double; var ReachedEps:Double;
		      var ErrCode:tErrorCode );

{ Поиск значений параметров произвольной функции, при которых наилучшим
образом (минимум суммы квадратов отклонений) аппроксимируется набор заданных
точек Xi,Yi}
procedure GeneralFitXX(aF:tFunc; { вычисляет значение аппроксимирующей
				(см. тип tFunc)}
		      aGradF:tGradF; { вычисляет значения всех частных
				     производных аппроксимирующей функции
				     (см. тип tGradF)}
		      var   P:tParameters; { ВХОД: значения начального
					   приближения для параметров;
					     ВЫХОД: значения найденного
					   наилучшего приближения для параметров.
					   (см. тип tParameters)}
		      const X,Y:tData;     { массивы данных для Xi и Yi
					   (см. тип tData)}
			    Eps:Double; { желаемая погрешность }
			    MaxIters:word; { максимальное число итераций,
					    если =0 используется умолчание }
		      var   ReachedEps:Double; { реально достигнутая погрешность }
		      var   ErrCode:tErrorCode { код ошибки });

{ Поиск значений параметров произвольной функции, при которых наилучшим
образом (минимум суммы квадратов отклонений) аппроксимируется набор заданных
точек Xi,Yi
! С возможностью исключения параметров из процесса изменения
см. аргумент SkippedP.
Другие аргументы см. GeneralFitXX }
procedure GeneralFitXXX(aF:tFunc; aGradF:tGradF;
		      var P:tParameters;
		      const SkippedP:tParametersSet; { Параметры,
				не включаемые в процесс оптимизации }
		      const X,Y:tData;
		      Eps:Double; MaxIters:word;
		      var   ReachedEps:Double;  var ErrCode:tErrorCode );

procedure AllocateP(var pP:tPtrParameters; Size:tParameterNumber);
procedure AllocateD(var pD:tPtrData;       Size:tDataNumber);

procedure FreeP(var pP:tPtrParameters);
procedure FreeD(var pD:tPtrData);

procedure MoveP(const p1:tParameters; var p2:tParameters);
procedure MoveD(const d1:tData;       var d2:tData);

{ Вычисление суммы квадратов отклонений }
function SumOfSqrDF(aF:tFunc; const P:tParameters; const X,Y:tData):Double;

IMPLEMENTATION

type
  tExtParsNumber=1..cMaxParameterSize;
  tExtParsArray=array[tExtParsNumber] of Extended;

  tExtPars=record
    Size:tExtParsNumber;
    Par:tExtParsArray;
  end;
  tPtrExtPars=^tExtPars;

function FullSizeP(Size:tParameterNumber):word;
  var x:tPtrParameters;
begin
  FullSizeP:=Size*SizeOf(x^.Par[1])+SizeOf(x^.Size);
end;

function FullSizeD(Size:tDataNumber):word;
  var x:tPtrData;
begin
  FullSizeD:=Size*SizeOf(x^.Data[1])+SizeOf(x^.Size);
end;

function FullSizeEP(Size:tExtParsNumber):word;
  var x:tPtrExtPars;
begin
  FullSizeEP:=Size*SizeOf(x^.Par[1])+SizeOf(x^.Size);
end;

procedure AllocateP;
begin
  GetMem(pP, FullSizeP(Size));
  if pP<>NIL then pP^.Size:=Size;
end;

procedure AllocateD;
begin
  GetMem(pD, FullSizeD(Size));
  if pD<>NIL then pD^.Size:=Size;
end;

procedure AllocateEP(var pEP:tPtrExtPars;  Size:tExtParsNumber);
begin
  GetMem(pEP, FullSizeEP(Size));
  if pEP<>NIL then pEP^.Size:=Size;
end;

procedure FreeP;
begin
  if pP<>NIL then FreeMem(pP, FullSizeP(pP^.Size));
  pP:=NIL;
end;

procedure FreeD;
begin
  if pD<>NIL then FreeMem(pD, FullSizeD(pD^.Size));
  pD:=NIL;
end;

procedure FreeEP(var pEP:tPtrExtPars);
begin
  if pEP<>NIL then FreeMem(pEP, FullSizeEP(pEP^.Size));
  pEP:=NIL;
end;

procedure MoveP;
begin
{$IfOpt R+}
    If p1.Size<>p2.Size then RunError(201);
{$EndIf}
  Move(p1.Par, p2.Par, SizeOf(p1.Par[1])*p1.Size);
end;

procedure MoveD;
begin
{$IfOpt R+}
    If d1.Size<>d2.Size then RunError(201);
{$EndIf}
  Move(d1.Data,d2.Data,SizeOf(d1.Data[1])*d1.Size);
end;

function VectorLengthSqr(const v: tParameters):Double;
{ VectorLengthSqr=(v*v)}
  var vl:Extended; i:word;
begin
  vl:=0;
  for i:=1 to v.Size do begin
    vl:=vl+Sqr(v.Par[i]);
  end;
  VectorLengthSqr:=vl;
end;

function VectorLength(const v: tParameters):Double;
{ VectorLength=Sqrt(v*v)}
begin
  VectorLength:=Sqrt(VectorLengthSqr(v));
end;

function FindMaxElementNum(const x:tParameters):word;
  var i:word; Max:Double;
begin
  FindMaxElementNum:=1;
  with X do begin
    Max:=Abs(Par[1]);
    for i:=2 to Size do if Abs(Par[i])>Max then begin
      FindMaxElementNum:=i; Max:=Abs(Par[i]);
    end;
  end;
end;

function SumOfSqrDF;
  var i:Word; z:Extended;
begin
{$IfOpt R+}
   If P.Size<1 then RunError(201);
   If P.Size>cMaxParameterSize then RunError(201);
   If (X.Size<>Y.Size) then RunError(201);
   If (X.Size<=P.Size-1) then RunError(201);
   If (X.Size>cMaxDataSize) then RunError(201);
{$EndIf}
  z:=0;
  for i:=1 to X.Size do begin
    z:=z+Sqr( Y.Data[i] - aF(P, X.Data[i]) );
  end;
  SumOfSqrDF:=z;
end;

procedure GeneralFit;
begin
  GeneralFitXXX(aF, aGradF, P,[], X,Y, Eps, cItMax, Eps, ErrCode);
end;

procedure GeneralFitX;
begin
  GeneralFitXXX(aF, aGradF, P,[], X,Y, Eps, cItMax, ReachedEps, ErrCode);
end;

procedure GeneralFitXX;
begin
  GeneralFitXXX(aF, aGradF, P,[], X,Y, Eps, MaxIters, ReachedEps, ErrCode);
end;

procedure GeneralFitXXX;

  { Вычисление функционала Ф(P)}
  function F(const P:tParameters):Double;
    var i:Word; z:Extended;
  begin
    z:=0;
    for i:=1 to X.Size do begin
      z:=z+Sqr( Y.Data[i] - aF(P, X.Data[i]) );
    end;
    F:=z/X.Size;
  end;

  var
    jMax:word;
    P1, Grad, S:tPtrParameters;
    ss:double;

  procedure DoPartialStep;
  begin
    P1^.Par[jMax]:=P.Par[jMax]-ss*Grad^.Par[jMax];
  end;

  procedure DoFullStep;
    var i:word;
  begin
    for i:=1 to P.Size do if not (i in SkippedP) then
      P1^.Par[i]:=P.Par[i]-ss*Grad^.Par[i];
  end;

  var pGr:tPtrExtPars;

  { Вычисление параметрического градиента функционала Ф(P)}
  procedure GradF(const P:tParameters; var G:tParameters);
    var i,j:word; F_minus_Y:Extended; tmp:double;
  begin
    for j:=1 to p.Size do pGr^.Par[j]:=0;
    for i:=1 to X.Size do begin
      aGradF(P, X.Data[i], G);
      F_minus_Y:=(aF(P, X.Data[i])-Y.Data[i]);
      for j:=1 to p.Size do if not (i in SkippedP) then
	pGr^.Par[j]:=pGr^.Par[j]+G.Par[j]*F_minus_Y;
    end;
    tmp:=2/X.Size;
    for j:=1 to p.Size do G.Par[j]:=pGr^.Par[j]*tmp;
  end;

  var
    ICounter:integer;
    LCounter,LCounterMax, CheckCnt:word;
    i:word;
    ExitRepeat:boolean;
    Func,Func1:Double;
    DP,DF:Extended;
    PrevEps, GLen: Double;

  procedure ProcessExitRepeat;
  begin
    ReachedEps:=GLen*ss;
    Dec(ICounter);
    Dec(CheckCnt);
    if CheckCnt=0 then begin
      ExitRepeat:=ExitRepeat or (PrevEps<ReachedEps);
      PrevEps:=ReachedEps;
      CheckCnt:=cCheckCnt;
    end;
    ExitRepeat:=ExitRepeat or (ReachedEps<=Eps) or (ICounter<=0);
  end;

begin
{$IfOpt R+}
   If P.Size<1 then RunError(201);
   If P.Size>cMaxParameterSize then RunError(201);
   If (X.Size<>Y.Size) then RunError(201);
   If (X.Size<=P.Size-1) then RunError(201);
   If (X.Size>cMaxDataSize) then RunError(201);
{$EndIf}
  AllocateP (P1,  P.Size);
  AllocateP (Grad,P.Size);
  AllocateEP(pGr, P.Size);
  AllocateP (S,   P.Size);
  MoveP(P,P1^);
  Func1:=F(P);
  for i:=1 to P.Size do S^.Par[i]:=0.01;
  ErrCode:=ecOK;
  ExitRepeat:=False;
  if MaxIters=0 then ICounter:=cItMax else ICounter:=cItMax;
  LCounterMax:=cItMax div 10;
  jMax:=1;
  PrevEps:=1/cFPMin;
  CheckCnt:=cCheckCnt;
  repeat
    Func:=Func1;
    P.Par[jMax]:=P1^.Par[jMax];
    GradF(P,Grad^);
    GLen:=VectorLength(Grad^);
    jMax:=FindMaxElementNum(Grad^);
{$IfDef DEBUG}
{    write(jMax:2);}
{$EndIf Def DEBUG}
    ss:=(1+1/7)*S^.Par[jMax];
    DoPartialStep;
    Func1:=F(P1^);

    if Func1>=Func then begin
      LCounter:=LCounterMax;
      repeat
	Dec(LCounter);
	ss:=ss/2;
	DoPartialStep;
	Func1:=F(P1^);
      until (Func1<=Func) or (LCounter=0);
      If (LCounter=0) then begin
	ErrCode:=ecNotCovering;
	ExitRepeat:=True;
      end;
      Dec(ICounter);
    end;

    { Параболическое уточненние шага }
    DP:=Sqr(Grad^.Par[jMax])*ss;
    DF:=(Func1+DP-Func);
    if (Abs(DF)>cFPMin) then begin
      ss:=0.5*ss*(DP/DF);
      DoPartialStep;
      Func1:=F(P1^);
    end;

    S^.Par[jMax]:=ss;
    ProcessExitRepeat;
  until ExitRepeat;

  with S^ do begin
    ss:=Par[1];
    for i:=2 to Size do if (not (i in SkippedP)) and (Par[i]<ss) then begin
      ss:=Par[i];
    end;
  end;
  CheckCnt:=cCheckCnt;
  ExitRepeat:=(ReachedEps<=Eps);
  while not ExitRepeat do begin
    Func:=Func1;
    MoveP(P1^,P);
    GradF(P,Grad^);
    GLen:=VectorLength(Grad^);
    DoFullStep;
    Func1:=F(P1^);

    { Параболическое уточненние шага }
    DP:=ss*Sqr(GLen);
    DF:=(Func1+DP-Func);
    if (Abs(DF)>cFPMin) then begin
      ss:=0.5*ss*DP/DF;
      DoFullStep;
      Func1:=F(P1^);
    end;

    ProcessExitRepeat;
  end;

  if ErrCode=ecOK then begin
    If ICounter<=0 then ErrCode:=ecMaxIterationNumberExeeded
    else if (ReachedEps>Eps) then ErrCode:=ecNotCovering;
  end;

{$IfDef DEBUG}
  writeln;
  writeln(' GenFit - Iterations:   ',cItMax-ICounter);
  writeln(' GenFit - SS:          ',ss);
  writeln(' GenFit - Grad length: ',GLen);
{$EndIf Def DEBUG}

  FreeP(P1);
  FreeP(Grad);
  FreeEP(pGr);
  FreeP(S);
end;

end.