unit ulearningcfmultiplano;
{$mode delphi}

{$DEFINE LAMBDA_AMORTIGUADOR}

interface

uses
  Classes, SysUtils,
  math,
  xmatdefs, utrayectorias_xcdpcf,uclasificador,
  uAdministradorDeEstadosMultiPlano, matreal,
  uhiperplanos;


function AjusteCFMultiplanos_1PLano_(
   CF: TAdminEstadosMultiPlano;
   trayectorias: TTrayectorias; learningRate: NReal ): NReal;

function AjusteCFMultiplanos_2PLano_(
   CF: TAdminEstadosMultiPlano;
   trayectorias: TTrayectorias; learningRate: NReal ): NReal;

{
devuelve:
-2 cuando no son suficientes ecuaciones para resolver el sistema
-1 cuando no hay puntos
 1 pudo resolver sin problemas
 2 no pudo resolver el sistema, pero devuelve un plano horizontal en CO
}
function CalcNew_a(
  var ahp: THiperplano;
  trayectorias: TTrayectorias;
  kPasoEnTrayectorias: integer  // 0..NPasos
):integer;

//FFUNCION QUE CLASIFICA LOS PUNTOS
//PONE DATOS NO VALIDOS EN LA TRAYECTORIA EN AQUELLOS LUGARES DONDE EL PUNTO NO PERTENECE A EL CONJUNTO
procedure Clasificador2conjuntos(var trayectorias1,trayectorias2:
  TTrayectorias; kPasoEnTrayectorias:Integer);


implementation

function AjusteCFMultiplanos_2PLano_(CF: TAdminEstadosMultiPlano;
  trayectorias: TTrayectorias; learningRate: NReal): NReal;
var
  kPaso: integer;
  ahp1, ahx1, ahp2, ahx2: THiperplano;
  nPuntosT_actualizables: integer;
  eMaxAbs, aMaxAbs: NReal;
  trayectorias1, trayectorias2: TTrayectorias;
begin
  nPuntosT_Actualizables:=  CF.nPuntosT - 1;
  if nPuntosT_actualizables > length( trayectorias.Trayectorias[0].puntos )  then
    nPuntosT_actualizables:= length( trayectorias.Trayectorias[0].puntos ) ;

  eMaxAbs:= 0;
  trayectorias1:=trayectorias.CreateClone;
  trayectorias2:=trayectorias.CreateClone;

  for kPaso:= 1 to nPuntosT_Actualizables do
  begin
    ahp1:= CF.conjuntosDeHiperPlanos[kPaso].hiperPlanos[0];
    ahx1:= ahp1.Clonar;
    ahp2:= CF.conjuntosDeHiperPlanos[kPaso].hiperPlanos[1];
    ahx2:= ahp2.Clonar;



    {
      Clasificamos los puntos (elementos de trayectorias) en dos conjuntos,
      trayectorias1 y trayectorias2, tienen el mismo largo que trayectorias y
      los puntos que no pertenecen están invalidados (el elemento
      trayectorias.Trayectorias[kTrayectoria].puntos[kPaso].valido:=False)
    }


    Clasificador2conjuntos(trayectorias1,trayectorias2, kPaso - 1 );

    {
      El CalcNew_A cuando el punto es inválido no lo considera
    }
    if CalcNew_A( ahx1, trayectorias1, kPaso-1 ) < 0 then
    begin
     ahx1.Free;
     ahx1:=nil;
    end;
    if CalcNew_A( ahx2, trayectorias2, kPaso-1 ) < 0 then
    begin
     ahx2.Free;
     ahx2:=nil;
    end;

    aMaxAbs:= DistanciaMaxRel( ahp1, ahx1 );
    if aMaxAbs > eMaxAbs then
      eMaxAbs:= aMaxAbs;
    aMaxAbs:= DistanciaMaxRel( ahp2, ahx2 );
    if aMaxAbs > eMaxAbs then
      eMaxAbs:= aMaxAbs;

    if ahx1 =nil then
    begin
      ahp1.Free;
      ahp1:=nil;
      ahp1:=ahp2.Clonar;
      ahp1.combinar( learningrate, ahx2 );
    end
    else
    begin
      ahp1.combinar( learningrate, ahx1 );
    end;

    if ahx2 =nil then
    begin
      ahp2.Free;
      ahp2:=nil;
      ahp2:=ahp1.Clonar;
      ahp2.combinar( learningrate, ahx1 );
    end
    else
    begin
      ahp2.combinar( learningrate, ahx2 );
    end;

  ahx1.free;
  ahx1:=nil;
  ahx2.free;
  ahx2:=nil;
  end;
  result:= eMaxAbs;
end;

function CalcNew_a(var ahp: THiperplano; trayectorias: TTrayectorias;
  kPasoEnTrayectorias: integer): integer;

var
  XtX, XtY: TMatR;
  kTrayectoria: integer;
  aTrayectoria: TTrayectoria;
  aPuntoTrayectoria: TPuntoTrayectoria; //   X: TEstadoX; CF: NReal;  cdp: NReal;
  X: TVectR;
  av: NReal;
  jFil, jCol: integer;
  //res: boolean;
  e10: longint;
  dimX, k, CantPuntosValidos: integer;
  CFMedio, ve_xi: NReal;
  desv_xi: ValReal;
  xi, yi: NReal;
  lambda_amortiguador: NReal;
  resEcaler: boolean;

begin
  dimX:=  length( ahp.A ) + 1;
  XtX:= TMatR.Create_init( dimX, dimX );
  XtX.Ceros;
  XtY:= TMatR.Create_init( dimX, 1);
  XtY.Ceros;

  Result:=0;
  CantPuntosValidos:=0;

  X:= TVectR.Create_init( dimX );
  X.pon_e( dimX, 1 );

  // agregamos los productos xtx y xty
  // de xtx solo el triángulo superior.
  for kTrayectoria:= 0 to trayectorias.nTrayectorias-1 do
  begin
    aTrayectoria:= Trayectorias.Trayectorias[kTrayectoria];
    aPuntoTrayectoria:= aTrayectoria.puntos[ kPasoEnTrayectorias ];
    //***ACA PREGUNTAR SI EL DATO ES VALIDO
    if aPuntoTrayectoria.valido then
    begin
      aPuntoTrayectoria.X.CopyToVectR( X ); // copia dimx-1 elementos y deja el 1 al final
      //***CONTADOR CON CANTIDAD DE DATOS VALIDOS
      inc(CantPuntosValidos);
      for jFil:= 1 to dimX do
      begin
        av:= X.e( jFil );
        XtX.acum_e(jFil, jFil, av*av );
        for jCol:= jFil + 1 to dimX do
          XtX.acum_e( jFil, jCol, av * X.e(jcol ) );
        XtY.acum_e( jFil, 1, av * aPuntoTrayectoria.CF );
      end;
    end;
  end;

  if CantPuntosValidos < dimX then
     Result:=-2 // devuelve -2 cuando no son suficientes ecuaciones para resolver el sistema
  else
    if CantPuntosValidos=0 then
       Result:=-1; //devuelve -1 cuando no hay puntos

  if Result < 0 then
    Exit;
  //***DIVIDIR ENTRE CANTIDAD DE DATOS VALIDOS
  CFMedio:= XtY.e( dimx, 1 )/CantPuntosValidos;

  //Promediamos y le damos el peso correspondiente al LearningRate
  XtX.PorRealTriangSup( 1/CantPuntosValidos );
  XtY.PorReal( 1/CantPuntosValidos );


  //XtX.WriteArchiXLT('XtX_'+ intToStr(kPasoEnTrayectorias) +'.xlt');
  //XtY.WriteArchiXLT('XtY_'+ intToStr(kPasoEnTrayectorias) +'.xlt');

  {$IFDEF LAMBDA_AMORTIGUADOR}
  lambda_amortiguador:= 1.0E-12;
  for k:= 1 to dimx do
    XtX.acum_e(k, k, lambda_amortiguador );
  {$ENDIF}

  XtX.SimetrizarTriangSup;

  XtX.Escaler( XtY, resEcaler, e10 );

  if resEcaler then
  begin
   for k:= 1 to dimx-1 do
      ahp.A[k-1]:= XtY.e( k, 1 );
   ahp.C0:= XtY.e( dimx, 1 );
   Result:=1;// pudo resolver sin problemas
  end
  else
  begin
    vclear( ahp.A );
    ahp.C0:= CFMedio;
    Result:=2;// no pudo resolver el sistema, pero devuelve un plano horizontal en CO
  end;


  XtX.Free;
  XtY.Free;
  X.Free;
end;

procedure Clasificador2conjuntos( var trayectorias1,
  trayectorias2: TTrayectorias; kPasoEnTrayectorias:Integer);
var
  kTrayectoria, kEstado, nEstados, nTrayectorias, nEstadosDiscretos,
    nEstadosContinuos, nEstadosMasCF: Integer;
    aTrayectoria: TTrayectoria;
  aPuntoTrayectoria: TPuntoTrayectoria;
  puntosTrayectoria: TVectR;
  puntosTrayectorias_kPaso:TMatR;
  clasificador: TClasificadorKMeansPorPlanos;
  conjuntos: TDAofMatR;
begin

  // Inicializo punteros para definir variables
  aTrayectoria:= trayectorias1.Trayectorias[0];
  aPuntoTrayectoria:= aTrayectoria.puntos[ kPasoEnTrayectorias ];
  // dichas variables
  nEstadosContinuos:= aPuntoTrayectoria.X.Xr.n;
  nEstadosDiscretos:= aPuntoTrayectoria.X.Xd.n;;
  nEstados:= nEstadosContinuos + nEstadosDiscretos;
  nTrayectorias:=trayectorias1.nTrayectorias;

  puntosTrayectoria:=TVectR.Create_Init(nTrayectorias);
  nEstadosMasCF:= nEstados + 1 ;// Ncontinuas + Ndiscretas + 1 (por el CF)

  {$IFDEF DEBUG}
  trayectorias1.WriteToXLTArchi('aux2HPfinal\trayectoria_'+  IntToStr(kPasoEnTrayectorias) + '.xlt');
  {$ENDIF}

  puntosTrayectorias_kPaso:=TMatR.Create_Init(nTrayectorias,nEstadosMasCF);

  for kTrayectoria:=1 to nTrayectorias  do
  begin
    aTrayectoria:= trayectorias1.Trayectorias[ kTrayectoria - 1 ];
    aPuntoTrayectoria:= aTrayectoria.puntos[ kPasoEnTrayectorias ];

    //Guardamos en puntosTrayectoria_kPaso =[Xr_1,....., Xr_nr, Xd_1,.....,Xd_nd, CF]
    puntosTrayectorias_kPaso.pm[kTrayectoria]:= TVectR.Create_Init(nEstadosMasCF);
    for kEstado:=1 to aPuntoTrayectoria.X.Xr.n do
      puntosTrayectorias_kPaso.pm[kTrayectoria].pon_e(kEstado,aPuntoTrayectoria.X.Xr.e(kEstado));
    for kEstado:=aPuntoTrayectoria.X.Xr.n + 1 to nEstadosMasCF-1 do
        puntosTrayectorias_kPaso.pm[kTrayectoria].pon_e(kEstado,aPuntoTrayectoria.X.Xd.e(kEstado - aPuntoTrayectoria.X.Xr.n));

    puntosTrayectorias_kPaso.pm[kTrayectoria].pon_e(nEstadosMasCF,aPuntoTrayectoria.CF);

  end;

   clasificador:=TClasificadorKMeansPorPlanos.Create(puntosTrayectorias_kPaso,2,kPasoEnTrayectorias);
   conjuntos:=clasificador.clasificar;

  {$IFDEF DEBUG}
   conjuntos[0].WriteArchiXLT('aux2HPfinal\Conjunto1_final'+IntToStr(kPasoEnTrayectorias)+'.xlt');
   conjuntos[1].WriteArchiXLT('aux2HPfinal\Conjunto2_final'+IntToStr(kPasoEnTrayectorias)+'.xlt');
  {$ENDIF}

  for kTrayectoria:=0 to nTrayectorias - 1 do
  begin
    if conjuntos[0].e(kTrayectoria + 1,1) = NO_VALIDO  then
      trayectorias1.Trayectorias[kTrayectoria].puntos[kPasoEnTrayectorias].anularPunto;
    if conjuntos[1].e(kTrayectoria + 1,1) = NO_VALIDO  then
      trayectorias2.Trayectorias[kTrayectoria].puntos[kPasoEnTrayectorias].anularPunto;
  end;

end;

function AjusteCFMultiplanos_1PLano_(CF: TAdminEstadosMultiPlano;
  trayectorias: TTrayectorias; learningRate: NReal): NReal;
var
  kPaso: integer;
  ahp, ahx: THiperplano;
  nPuntosT_actualizables: integer;
  eMaxAbs, aMaxAbs: NReal;
begin
  nPuntosT_Actualizables:=  CF.nPuntosT - 1;
  if nPuntosT_actualizables > length( trayectorias.Trayectorias[0].puntos )  then
    nPuntosT_actualizables:= length( trayectorias.Trayectorias[0].puntos ) ;

  eMaxAbs:= 0;

  for kPaso:= 1 to nPuntosT_Actualizables do
  begin
    ahp:= CF.conjuntosDeHiperPlanos[kPaso].hiperPlanos[0];
    ahx:= ahp.Clonar;
    CalcNew_A( ahx, trayectorias, kPaso-1 );
    aMaxAbs:= DistanciaMaxRel( ahp, ahx );
    if aMaxAbs > eMaxAbs then
      eMaxAbs:= aMaxAbs;
    ahp.combinar( learningrate, ahx );
    ahx.free;
  end;
  result:= eMaxAbs;
end;

end.

