unit usteepestdescent_quad;

(*==============================================================================
  ==============================================================================
 Algoritmo de steepest descent. Refs [1], [2]
 Este algoritmo busca el mínimo de una función objetivo N-dimensional
 en la dirección negativa del gradiente:
 -g(x1,...,xN) = -grad_f(x) = -(df/dx1,...,df/dxN)'
 con el tamaño de paso alfa_k (en la iteración k) ajustado de tal forma
 que la función objetivo se minimiza a lo largo de la dirección por un
 método de búsqueda unidimensional (line search).

 -------------------------------------------------------------------------------
 descripción del algoritmo:
 paso 0)  con el número de iteración k=0, encontrar f0=f(x0):
          k      <- 0
          x_k    <- x
          f_xk   <- f_x0
 paso 1)  encontrar el tamaño de paso alfa_k a lo largo de la dirección del
          gradiente negativo por el método de aproximación cuadrática:
          g_k    <- -grad f_xk
          alfa_k <- argmin_alfa f(x_k - alfa*g_k/norm(g_k))
 paso 2)  obtener el próximo punto y actualizar f:
          x_ka   <- x_k
          f_ka   <- f_xk
          x_k    <- x_k - alfa*g_k/norm(g_k);
          f_xk   <- f(x_k)
 paso 3)  si x_ka \approx x_k y f_ka \approx f_xk, entonces declarar
          a x_k el mínimo y terminar. de lo contrario incrementar k y volver
          al paso 1.
 -------------------------------------------------------------------------------
 Refs.:
 [1] Luenberger
 [2] Applied Numerical Methods Using MATLAB

================================================================================
================================================================================
*)

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils, xMatDefs, matreal,
  ubibfun_onl;

procedure optimizar(f: TFunc_Bloque;
                    nvars : NInt;
                    var XOpt : TVectR;
                    x0: TVectR;
					TolX: NReal;
					TolFun: NReal;
					alfa0: NReal;
					MaxIter: NInt);

function linesearch(f: TFunc_Bloque;
                    x_k: TVectR;
                    alfa0: NReal;
                    g_k: TVectR;
                    MaxIter: NInt;
                    TolX: NReal;
                    TolFun: NReal): NReal;
var
 {$IFDEF DBG}
         fileOut : Text;
 {$ENDIF}
implementation

procedure optimizar(f: TFunc_Bloque;
                    nvars : NInt;
                    var XOpt : TVectR;
                    x0: TVectR;
					TolX: NReal;
					TolFun: NReal;
					alfa0: NReal;
					MaxIter: NInt);
var
 f_xk, f_xka : NReal;
 alfa_k: NReal;
 k: NInt;
 g_k: TVectR;
 x_k: TVectR;
 x_ka : TVectR;
 n: NInt;
 i: NInt;
 g_k_normEuclid: NReal;
 CONVERGIO : Boolean;
begin
     //algunas inicializaciones
     n := f.nvars;
     x_k := TVectR.Create_Clone(x0);
     f_xk := f.Evaluar(x0);
     alfa_k:=alfa0;
     k:=0;
{$IFDEF DBG}
     //imprimir en consola
     writeln('--------------------------------------------------------------------------------');
     writeln('paso = ', k,'  x_k(1) = ',x_k.e(1):1:6, '  x_k(2) = ', x_k.e(2):1:6,'  f_xk = ', f_xk:1:6,'  alfa_k = ',alfa0:1:6);
     writeln('--------------------------------------------------------------------------------');
     //tambien imprimir en archivo
     Assign (fileOut, 'a.txt');
     Rewrite (fileOut);
     writeln(fileOut,'paso    ','x_k(1)    ', 'x_k(2)    ','f_xk    ','alfa_k     ');
     writeln(fileOut,k,'    ',x_k.e(1):1:6,'    ',x_k.e(2):1:6,'    ',f_xk:1:6,'    ',alfa0:1:6);
{$ENDIF}
     repeat begin
            g_k := TVectR.Create_Clone(f.Gradiente(x_k));
            g_k_normEuclid := g_k.normEuclid;
            for i:= 1 to n do   //normalizar
                g_k.pon_e(i, g_k.e(i)/g_k_normEuclid);

            //LINE SEARCH : hallar alfa_k = argmin_alfa f(x_k - alfa*g_k)
            alfa_k := linesearch(f, x_k, alfa_k, g_k,
                                 MaxIter, TolX, TolFun);

            x_ka:=TVectR.Create_Clone(x_k);
            f_xka:=f_xk;
            for i:=1 to n do
                x_k.pon_e(i, x_k.e(i) - alfa_k * g_k.e(i));

            f_xk:=f.Evaluar(x_k);

            CONVERGIO := (x_ka.distancia(x_k) < TolX) and (abs(f_xka-f_xk) < TolFun); //cambiando "and" por "or" se requiere convergencia mas 'dèbil'

            k:=k+1;
{$IFDEF DBG}
        writeln('paso = ', k,'  x_k(1) = ',x_k.e(1):1:6, '  x_k(2) = ', x_k.e(2):1:6,'  f_xk = ', f_xk:1:6,'  alfa_k = ',alfa_k:1:6);
        writeln('-----------------------------------------------------------------------------');
        writeln(fileOut,k,'    ',x_k.e(1):1:6,'    ',x_k.e(2):1:6,'    ',f_xk:1:6,'    ',alfa_k:1:6);
{$ENDIF}
     end until ((CONVERGIO) or (k > MaxIter));
     xOpt := TVectR.Create_Clone(x_k);

{$IFDEF DBG}
     Close(fileOut);
{$ENDIF}

end;

function linesearch(f: TFunc_Bloque;
                    x_k: TVectR;
                    alfa0: NReal;
                    g_k: TVectR;
                    MaxIter: NInt;
                    TolX: NReal;
                    TolFun: NReal): NReal;
var
 num, den : NReal;
 alfa_opt,x0,x1,x2,x3:NReal; //representan valores del parámetro real alfa
 f0,f1,f2,f3:NReal;  //evaluaciones de f(x_k - alfa*g_k), es una función de R->R en alfa
 vaux: TVectR;
 k: NInt;
begin
     k:=MaxIter;

     //generar el initial guess
     x0:=alfa0;
     x1:=alfa0*1.5;
     x2:=alfa0*2.0;

     vaux:=TVectR.Create_Clone(g_k); vaux.PorReal(-x0); vaux.Suma(x_k,vaux);
     f0:=f.Evaluar(vaux);

     vaux:=TVectR.Create_Clone(g_k); vaux.PorReal(-x1); vaux.Suma(x_k,vaux);
     f1:=f.Evaluar(vaux);

     vaux:=TVectR.Create_Clone(g_k); vaux.PorReal(-x2); vaux.Suma(x_k,vaux);
     f2:=f.Evaluar(vaux);

     //iteración
     repeat
           //me calculo los nuevos x3, f3
           num:=f0*(x1*x1-x2*x2)+f1*(x2*x2-x0*x0)+f2*(x0*x0-x1*x1);
           den:=2*(f0*(x1-x2)+f1*(x2-x0)+f2*(x0-x1));

           if (den <> 0) then
              x3:=num/den;

           vaux:=TVectR.Create_Clone(g_k); vaux.PorReal(-x3); vaux.Suma(x_k, vaux);
           f3:=f.Evaluar(vaux);

           if (k <= 0) or (abs(x3-x1) < TolX) or (abs(f3-f1) < TolFun) then begin
              break;
           end
           else begin
                if (x3 < x1) then
                       if (f3 < f1) then begin
                          f0:=f0; f2:=f1; f1:=f3;
                          x0:=x0; x2:=x1; x1:=x3;
                       end else begin
                          f0:=f3; f1:=f1; f2:=f2;
                          x0:=x3; x1:=x1; x2:=x2;
                       end
                else //x3 > x1
                       if (f3 < f1) then
                       begin
                          f0:=f1; f1:=f3; f2:=f2;
                          x0:=x1; x1:=x3; x2:=x2;
                       end else begin
                          f0:=f0; f1:=f1; f2:=f3;
                          x0:=x0; x1:=x1; x2:=x3;
                       end;
           end;
           k:=k-1;
     until (false);  //el escape está en el if de arriba

     x0:=x3;
     f0:=f3;
     alfa_opt:=x0;
     Result := alfa_opt;
end;

end.
