unit uexplorador_eg;

{$mode delphi}

interface

uses
  Classes, SysUtils,
  {$IFDEF USE_DB_CON}
  lib_dbmysql,
  {$ELSE}
  uDataSetGenerico,
  urosx,
  {$ENDIF}
  {$IFDEF data_NO_}
  udbconpg,
  {$ENDIF}
  math,
  uoddface,
  utipos_ga,
  fddp_conmatr,
  fddp,
  xmatdefs, matreal, matent, uexplorador;


type

  { TExplorador_EG }

  TExplorador_EG = class( TExplorador )
    grad: TVectR;
    paso: NReal;
    f0: NReal;
    cnt_OkSeguidos, cnt_FallosSeguidos: integer;

    constructor Create(dbconx: TDBrosxCon; Problema: TProblema );
    function ProponerNuevoIndividuo: TIndividuo; override;
    destructor destroy; override;
  end;

(********************

Xs = X0 - paso * G0  (ec.1)

Fs = evaluar( Xs )

Ahora bien, seguramente Fs <> F0 + G0 ( Xs-X0 )
Para estimar un nuevo paso, supondremos que el gadiente G0 evoluciona
durante el paso variando de forma tal que el producto G0( Xs - X0 )
cambia del valor L0 = G0( Xs -X0 ) al valor Ls = Gs (Xs - X0).

G ( X - X0 ) = L0 * ( 1 - alfa ) + Ls * alfa  (ec.1)
donde alfa = 0 implica X = X0 y alfa = 1 implica X = Xs.

suponiendo esa variación de la poyección del gradiente sobre el recorrido de X
el valor Fs es calculable cómo:
Fs - F0 = integral( G dX ) (ec.2)

X - X0 = alfa ( Xs-X0 )  (ec.3)
con lo cual dX = (Xs - X0) dAlfa  y haciendo el cambio de variables, la integral
se puede escribir como:
Fs - F0 = integral( ( L0 (1-alfa) + Ls alfa ) (Xs-X0 ) dAlfa )  (ec.4)

Fs - F0 = ( L0 ( alfa - 0.5 alfa² )+ 0.5 Ls alfa² ) ) (Xs-X0) | [1, 0]
(Fs - F0)/(Xs-X0) = L0 0.5 + 0.5 Ls de donde podemos despejar Ls como:
Ls = 2 (Fs-F0)/(Xs-X0) - L0     (ec.5)

Con esa estimación de Ls, podemos calcular el alfa que anule el integrando de
la ec.4. (con lo cuál se determina un mínimo o un máximo de Ls)

L0 ( 1- alfa ) + Ls alfa = 0
.: L0 = (L0 - Ls ) alfa
.: alfa = L0 / ( L0 - Ls)  (ec.5)

Observaciones.
1) Nos interesa buscar el mínimo de F.
1.1) si L0 = Ls no hay variación detectable del gradiente y simplemente se puede
    continuar con el paso actual o incrementarlo levemente.

2) Por construcción L0 < 0 (pues se buscó un movimiento de X con -Gradiente)
3) Si Ls < 0 quiere decir, que al final del intervalo continúa bajando
3.1) abs( Ls ) < abs( L0)  => alfa  > 1  (de la ec.5)
     Se podría incrementar el paso sin supear el valor correspondiente a Alfa.
3.2) abs( Ls ) > abs( L0 )
     Significa que la función decrece más rápido en la direción del movimiento
     con lo cual se podría aumentar el paso
4) Si Ls > 0 , se llegó a un punto dónde parecería que nos pasamos del mínimo
     y se cumplira 0 < alfa < 1 con lo cual tenemos una propuesta para el nuevo
     paso.
     Estimando el nuevo punto a partir de Xs se tendría:
      Xss = Xs - NuevoPaso Gs
      siendo NuevoPaso = (1-alfa) * Paso
     Estimando el nuevo punto a partir de X se tendría:
      Xss = X - alfa * Paso * G0

     Como Gs y G0, son efectivamente calculados, seguramente no resulten co-lineales
     y por tanto ambas formas de aproximación no resultan necesariamente en el mismo
     resultado; por lo cual, podría considerarse el promedio de ambas.




******************)
implementation


function TExplorador_EG.ProponerNuevoIndividuo: TIndividuo;
var
  sql: string;
  ds: TResultadoQuery;
  r: TDataRecord;
  a: TIndividuo;
  xxt: TMatR;
  xf: TMatR;
  x0, x: TVectR;
  k, j, DimXR, DimXE: integer;
  x_k: NReal;

  {$IFDEF EXPLORAR_ZONAS_OSCURAS}
  res_det: NReal;
  res_invertible: boolean;
  res_exp10: integer;
  dim_nucleo: integer;
  fr: NReal;
  maux: TMatR;
  nucleo_xxt: TMatR;
  {$ENDIF}


  DimX, res_exp10: integer;
  beta, res_det: NReal;
  res_invertible: boolean;
  nPoblacion: integer;
begin

  nPoblacion := Problema.CantidadDeIndividuos;
  if nPoblacion < 200 then
  begin
    a:= rnd_Individuo;
    problema.decodificar_adn( a );
    if problema.viable(a) then
      Result := a
    else
    begin
      a.Free;
      Result := nil;
    end;
    exit;
  end;

  DimXR:= length( problema.DescriptoresR );
  DimXE:= length( problema.DescriptoresE );
  DimX:= DimXR + DimXE;
  // matrices para el problema lineal de mínimos cuadrados
  xxt:= TMatR.Create_Init( DimX+1, DimX + 1);
  xf:= TMatR.Create_Init( DimX +1, 1 );

  // vector para cargar muestras
  x:= TVectR.Create_Init( DimX +1 );
  x0:= TVectR.Create_Init( DimX );

writeln( 'cnt:', cnt_evals,', oks:', cnt_OkSeguidos, ', fas:', cnt_FallosSeguidos,', f0:', f0,', fl:', f_VE_Last, ', Paso: ', paso);
  if self.cnt_evals >= 1 then
  begin
    if self.f_VE_Last > f0 then
    begin
      cnt_OkSeguidos:= 0;
      inc( cnt_FallosSeguidos );
      if cnt_FallosSeguidos > 1 then
      begin
        paso:= paso * InvGoldenRatio;
        cnt_FallosSeguidos:= 0;
      end;
    end
    else if self.f_VE_Last < f0 then
    begin
      cnt_FallosSeguidos:= 0;
      inc( cnt_OkSeguidos );
      if cnt_OkSeguidos > 2 then
      begin
        paso:= paso * GoldenRatio;
        cnt_OkSeguidos:= 0;
      end;
    end;
  end;


    sql:= 'SELECT * FROM '+problema.db_table+
        ' WHERE cnt_evaluaciones > 0 ORDER BY f_objetivo LIMIT 1000 ';
    ds:= dbcon.sql_query( sql );

    ds.First;
    r:=ds.getFilaActiva;

    a := TIndividuo.CreateFromRec(dbcon, problema, r);
    problema.decodificar_adn( a );


  // cargamos la muestra en X0
  for k:= 1 to DimXR do
    x0.pon_e( k, a.XR.e( k ) );
  for k:= 1 to a.XE.n do
    x0.pon_e( k + DimXR, a.XE.e( k ) );

  f0:= a.f_objetivo;
  a.Free;

  while not ds.EOF do
  begin
    r:=ds.getFilaActiva;
    a := TIndividuo.CreateFromRec(dbcon, problema, r);
    problema.decodificar_adn( a );

    // cargamos la muestra en X
    for k:= 1 to DimXR do
      x.pon_e( k, a.XR.e( k ) - x0.e(k) );
    for k:= 1 to a.XE.n do
      x.pon_e( k + DimXR, a.XE.e( k ) - x0.e( DimXR + k ) );
    x.pon_e( x.n, 1 );



    // agregamos la muestra X
   for k:= 1 to x.n do
   begin
    x_k:= x.e(k);
    xxt.acum_e( k, k, sqr( x_k ));
    for j:= k+1 to x.n do
      xxt.acum_e(k, j, x_k * x.e( j ) );
    xf.acum_e( k, 1, x_k * (a.f_objetivo - f0) );
   end;
   a.Free;
   ds.Next;
  end;

  ds.Free;

  // Ahora completamos el tríangulo inferior de XXt.
  for k:= 2 to xxt.nf do
    for j:= 1 to k-1 do
       xxt.pon_e( k, j, xxt.e( j, k ) );


{$IFDEF EXPLORAR_ZONAS_OSCURAS}
// Esto lo implementé, pero después me arrepentí. Como en la consulta
// se trae on conjunto de Los Mejores, puede ocurrir que en dicho conjunto
// determinadas dimensiones del espacio de búsqueda estén FIJAS pues ya se
// determinó el valor óptimo y por tanto no ayuda que este explorador inyecte
// ruido en las mismas. Este mecanismo de inyección de ruido, podría tener
// sentido con otra implementación que seleccione un conjunto de individuos en base
// a otro criterio que no sea "Los Mejorcitos". O simplemente, si se utilizan TODOS
// los individuos.
//
// Calculamos una base del NUCLEO para poder "inyectar ruido" para explorar
// la zona no explorada.
  maux:= TMatR.Create_Clone( xxt );
  dim_nucleo:= maux.CalcBSEN( nucleo_xxt );

  if dim_nucleo < xxt.nf then
  begin
    // bien, ahora sumamos algo a la diagonal para regularizar.
    for k:= 1 to xxt.nf do
      for j:= 1 to xxt.nc do
        xxt.pon_e(  k, j , max( 1e-6, xxt.e( k,j) * 1.0001 ) );

    res_det:= xxt.Escaler( xf, res_invertible, res_exp10 );
    // en las primeras componentes de xf tenemos el gradiente, en la última
    // el término constante.
    for k:= 1 to xf.nf-1 do
      grad.pon_e( k, xf.e(k, 1 ) );
  end
  else
  begin
   // No hay información del gradiente, será todo Random.
    grad.Ceros;
  end;

  // ahora inyectamos el rudio
  // las filas de nucleo_xxt son vectores de la base del núlceo de xxt
  for k:= 1 to nucleo_xxt.nf do
  begin
   fr:= 2* fuente.rnd -1 ; // va de -1 a +1
   for j:= 1 to DimXR do
     grad.acum_e(j, nucleo_xxt.e( k, j ) * fr * ( problema.DescriptoresR[j-1].x_max - problema.DescriptoresR[j-1].x_min )/4 );
   for j:= 1 to DimXE do
     grad.acum_e(DimXR+j,  nucleo_xxt.e( k, DimXR+j ) * fr * ( problema.DescriptoresE[j-1].K_max - problema.DescriptoresE[j-1].k_min )/4 );
  end;

  maux.Free;
  nucleo_xxt.Free;


{$ELSE}


  // Bien, aquí, tendríamos que sumar alfa * I a xxt para "regularizar" y que la solución
  // exista imponiendo un CERO en las dimensiones que no es "fuerte" el gadiente.
  // En definitiva, que en las direcciones que no tiene información NO intente moverse.

  // buscamos el máximo elemento diagonal
  beta:= xxt.e(1,1);
  for k:= 2 to xxt.nf do
    if xxt.e(k,k) > beta then beta := xxt.e(k,k);

  // sumamos en la diagonal 0.1/1000 del máximo o 1E-10, lo que resulte
  // mayor.
  beta:= max( 1.0e-10, beta * 1e-4 );
  for k:= 1 to xxt.nf do
     xxt.acum_e(k,k, beta );

   res_det:= xxt.Escaler( xf, res_invertible, res_exp10 );
   // en las primeras componentes de xf tenemos el gradiente, en la última
   // el término constante.
   for k:= 1 to xf.nf-1 do
      grad.pon_e( k, xf.e(k, 1 ) );
{$ENDIF}

  a:= TIndividuo.CreateNew( dbcon, problema );
  for k:= 1 to DimXR do
   a.XR.pon_e( k, X0.e( k ) - paso * grad.e( k ) );
  for k:= 1 to DimXE do
   a.XE.pon_e( k, round( X0.e( DimXR+k ) - paso * grad.e( DimXR+k ) ) );

  problema.encajonar( a );
  problema.codificar_adn( a );
  problema.decodificar_adn( a );


  // Aquí tenemos un sistema xxt a = xf
  // el sistema permite aproxima f = a' x  donde x es el conjunto de parámetros y un
  // 1 al final
  //
  {
  intuitivamente, como xxt puede no expandir el espacio completo, el sistema
  puede no ser invertible. En ese caso, sea W la matriz cuyas columnas son
  autovectores (de norma 1) de xxt asociados a los autovalores l1 .. ln

  xxt W = diag( l1, ... ln )

  aqui la tengo que pensar un poco, pero hay que calcular una base ortonormal
  de los autovectores de xxt. Los que estan asociados al núlceo para mi
  expanden "La Zona Oscura" del espacio ...mientras que los asociados
  a utovalores <> 0 (>0 por construcción) explanden la parte explorada
  "La Zona Iluminada" del espacio.
  El gradiente de la función estimado solo aporte información sobre la Zona Iluminada

  Para dal el paso, la estrategia sería
   x_s = ( Wt xo ) - ( Wt dFdx )* pasoGrad + sum( j , wj * rnd( xmin, xmax)  ) * AtenuadorRuido

   y vuelta a calcular con una estrategia de variación de pasoGrad.

  }

  X0.Free;
  x.Free;

  if problema.viable(a) then
    Result := a
  else
  begin
    a.Free;
    Result := nil;
  end;

end;

destructor TExplorador_EG.destroy;
begin
  grad.Free;
  inherited destroy;
end;



{ TExplorador_EG }

constructor TExplorador_EG.Create(dbconx: TDBrosxCon; Problema: TProblema);
begin
  inherited Create(dbconx, Problema);
  cnt_OkSeguidos:= 0;
  cnt_FallosSeguidos:= 0;
  paso:= 1;
  grad:= TVectR.Create_init( length( Problema.DescriptoresE) + length( Problema.DescriptoresR ) );
end;



end.
