{$A+,B-,D+,E+,F+,G-,I+,L+,N+,O-,P-,Q-,R-,S+,T-,V+,X+}

unit MinRKF;
{
  Minimizador de Fun‡äes usando RKF
}

INTERFACE

uses RKF;

type

  TMinimizando = procedure(x: PVector; var f:real; gradf:PVector);
  {
    Dado um estado x, retorna o valor da fun‡Æo sendo minimizada
    em f e seu gradiente em gradf
  }

  P = ^T;
  T = object
      constructor init(An: integer); {An ‚ o numero de variaveis do sistema}
      destructor finish;
      procedure minimize(
            f:TMinimizando;
            Ax0,Ax:PVector;
            dtmin,dtmax,tol,limit:real
            );
      {
       Minimiza a funcao "f(x)" tomando inicialmento o estado "Ax0",
       retornando em "Ax" o estado para o qual
         "gradiente(f)<limit"
       tomando como parametros "dtmin" e "dtmax" que ‚ o intervalo
       em que necessariamente estar  o passo de tempo "dt" fazendo
       o poss¡vel para que o erro nÆo ultrapasse "tol" durante a
       integra‡Æo.
      }
      procedure SetDebugProcedure(AUseState:TUseState);
      PRIVATE
        G:RKF.T;
        UseState: TUseState;
      end;

function NormOfVector(n:integer; v:Pvector):real;


IMPLEMENTATION


function NormOfVector(n:integer; v:Pvector):real;
var i:integer;
    sumsqr:real;
begin
  sumsqr:=0.0;
  for i:=0 to n-1 do sumsqr:=sumsqr+sqr(v^[i]);
  NormOfVector:=sqrt(sumsqr);
end;

VAR FCT:TMinimizando;
    ZLimit:real;
    ZN:integer;

PROCEDURE Derivative(t: REAL; n: integer; s, v: PVector); FAR;
   {
     Dado um tempo "t" e o estado correspondente "s[0..n-1]",
     calcula a derivada do estado "v[0..n-1]".
   }
 var i:integer;
     k:real;
 begin
   FCT(s,k,v); {v contem gradiente de s}
   k:=NormOfVector(n,v);
   if k < ZLimit then begin
     k:=0.0;
   end
   else
     k:=-1.0/k;
   for i:=0 to n-1 do
     v^[i]:=v^[i]*k;
 end;

FUNCTION DetectEvent(t0, t1: REAL; n: integer; s0, s1: PVector): REAL; FAR;
   {
     Dados dois tempos "t0" e "t1", e os estados correspondentes
     "s0[0..n-1]" e "s1[1..n-1]", determina o instante "t"
     no intervalo "[t0 __ t1]" onde ocorre o proximo evento. }
   var x:PVector;
       l:real;
   begin
     DetectEvent:=0.0;
     x:=NEWVECTOR(n);
     FCT(s1,l,x);
     if NormOfVector(n,x)<Zlimit then DetectEvent:=0.5*(t1+t0);
     GETRIDOFVECTOR(n,X);
   end;

PROCEDURE DefaultUseState(t: REAL; n: integer; s: PVector); FAR;
    {
      Usa o estado "s[0..n-1]" correspondente ao tempo "t". }
   begin
   end;

constructor T.init;
begin
  G.init(An);
  Zn:=An;
  UseState:=DefaultUseState;
end;

destructor  T.finish;
begin
  G.finish;
end;

procedure T.SetDebugProcedure;
begin
  UseState:=AUseState;
end;

procedure T.minimize(
            f:TMinimizando;
            Ax0,Ax:PVector;
            dtmin,dtmax,tol,limit:real
            );
var FinishingTime:real;
begin
  FCT:=f;
  ZLimit:=limit;
  CopyState(Zn, Ax0, Ax);
  FinishingTime:= G.Integrate(0.0, 1.0, dtmin, dtmax, tol,
                     Ax, Derivative, DetectEvent, UseState);
  if FinishingTime>=1.0 then writeln('MAIOR QUE 1'#7);
end;

begin
end.