unit linesearchCuadratico;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils, MatReal, xMatDefs;

function evaluar(x: NReal): NReal;
function opt_quad(x0, x1, x2: NReal; f0, f1, f2: NReal; TolX, TolFun: NReal; kmax: NInt;var kniters: NInt): NReal;
procedure opt_quad0(var x0,x1,x2: NReal; var f0,f1,f2: NReal; TolX, TolFun: NReal; var k: NInt; var xopt: NReal);

implementation

function evaluar(x: NReal): NReal;
begin
     Result:=(x*x-4)*(x*x-4)/8-1;      //eq. p325   [2] Applied Numerical Methods Using MATLAB
     //Result:=exp(-x+3)*sin(x) + x*x/17;
end;

function opt_quad(x0, x1, x2: NReal; f0, f1, f2: NReal; TolX, TolFun: NReal; kmax: NInt;var kniters: NInt): NReal;
var
   xopt: NReal;
   k: NInt;
begin
     //control de concavidad, otros.

     //optimizar
     k:=kmax;
     opt_quad0(x0,x1,x2,f0,f1,f2,TolX,TolFun,k,xopt);
     kniters:=kmax-k+1;
     Result:=xopt;
end;

procedure opt_quad0(var x0,x1,x2: NReal; var f0,f1,f2: NReal; TolX, TolFun: NReal; var k: NInt; var xopt: NReal);
var
  num, den : NReal;
  x3, f3: NReal;
begin
     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));
     x3:=num/den;
     f3:=evaluar(x3);

     writeln('paso ',k, '  :');
     writeln('            x0 = ',x0:3:3,' x1 = ',x1:3:3,' x2 = ',x2:3:3);
     writeln('            f0 = ',f0:3:3,' f1 = ',f1:3:3,' f2 = ',f2:3:3);

     if (k <= 0) or (abs(x3-x1) < TolX) or (abs(f3-f1) < TolFun) then
     begin
          x0:=x3;
          f0:=f3;
          xopt:=x0;
     end
     else begin
          if (x3 < x1) then
             if (f3 < f1) then
             begin
                  f0:=f0; f1:=f3; f2:=f1;
                  x0:=x0; x1:=x3; x2:=x1;
             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;

          k:=k-1;
          opt_quad0(x0,x1,x2,f0,f1,f2,TolX,TolFun,k,xopt);
     end;
end;

end.
