unit uclasificador;

{$mode delphi}

interface

uses
  Classes, SysUtils, matreal, xmatdefs, math;

const
  NO_VALIDO = 77777;

type


  TDAOfDAOfVectR = array of TDAOfVectR ;



  { TClasificador }

  TClasificador = class
     NConjuntos: integer;
     elementos:TMatR;
     NElementos:integer;
     constructor Create(xElementos: TMatR; xNConjuntos: integer);
     procedure Free;

  end;

  { TClasificadorKMeansPorPlanos }

  TClasificadorKMeansPorPlanos = class(TClasificador)

    planos: TDAOfVectR;
    kPaso: integer;
    constructor Create(xElementos: TMatR; xNConjuntos, xkPaso: integer);
    function clasificar:TDAofMatR;

    function distanciaEuclideaVectR(punto1,punto2: TVectR): NReal;

    {P=(X1*,..., Xn-1*,Y*) --> TVectR = P
    Plano ->   A0 + A1X1 + A2X2 + .....+ An-1Xn-1* + BY = 0 --> TVectR = ( A1, A2,....,An-1, A0, B) RESPETAR EL ORDEN, EL TERMINO INDEP VA PENULTIMO
    Distancia =  abs(A0 + A1X1* + A2X2* + .....BY*)/sqrt(A1^2 + .... + B^2))}
    function distanciaEuclideaAPlano(punto, plano: TVectR): NReal;

    function aproxPlano(elementos: TMatR; out invertible: Boolean): TVectR;

  end;


implementation

{ TClasificadorKMeansPorPlanos }

constructor TClasificadorKMeansPorPlanos.Create(xElementos: TMatR; xNConjuntos,
  xkPaso: integer);
begin
  inherited Create(xElementos,xNConjuntos);
  SetLength(planos,NConjuntos);
  kPaso:=xkPaso;
end;

function TClasificadorKMeansPorPlanos.clasificar: TDAofMatR;
var
  res:TDAofMatR;
  minimo, maximo, largoIntervalo, distancia,distancia_k, coef, dist_ant,
     CF_max, CF_min, C_medio, C_mediana:NReal;
  plano_new, aElemento, plano_repetido, MedianasdeElementos:TVectR;
  kPlano, kElemento, kPlanoElemento, kCoeficiente, NCoeficientes,
    kCoordenada, kIter, kPlano_Ant: Integer;
  salir, esNil, planosSimilares, invertible: Boolean;
  lst_elementos: TList;
  aElemento_puntero: Pointer;

  function Comparar(Item1, Item2: pointer): integer;
  var
    res:integer;
    a1,a2: TVectR;
  begin
    a1 := Item1;
    a2 := Item2;
    if a1.pv[a1.n] > a2.pv[a2.n] then
      res := 1
    else if a1.pv[a1.n] = a2.pv[a2.n] then
      res := 0
    else
      res := -1;
    Result:=res;
  end;


begin

  NCoeficientes:=elementos.nc + 1; // //elementos.nc=cant variables de estado + 1 (corresponde al CF), y se le suma + 1 (corresponde al termino independiente)

  plano_repetido:=nil;

  SetLength(res, NConjuntos);
  for kPlano:=0 to NConjuntos-1 do
  begin
    res[kPlano]:=TMatR.Create_Init(NElementos, NCoeficientes - 1); //sin termino independiente
    for kElemento:=1 to NElementos  do
      res[kPlano].pm[kElemento]:=nil // Dato no valido
  end;
  distancia:=MaxNReal;

  // Primero elegimos NConjuntos planos inicales con algún criterio
  // En principio RAAAAANDOMMMMM

  RandSeed:=33; //semilla inicial del random

  SetLength(planos,NConjuntos);

  aElemento := elementos.pm[10];
  {$IFDEF DEBUG}
  //////////////////////
  //aElemento.WriteXLTSimple_archi('aux2hp\punto.xlt');
  elementos.WriteArchiXLT('aux2hp\puntosTodos'+IntToStr(kPaso)+'.xlt');
  //////////////////////////////
  {$ENDIF}

  //****ELECCION DE PLANOS INICIALES DE FORMA ALEATORIOS
{  for kPlano:=0 to NConjuntos - 1 do
  begin
    coef:=0;
    planos[kPlano]:= TVectR.Create_Init(NCoeficientes);
    for kCoeficiente:=1 to NCoeficientes - 2 do
    begin
      planos[kPlano].pon_e(kCoeficiente,Random(10000000)-500000);
      coef:= coef + planos[kPlano].e(kCoeficiente)*aElemento.e(kCoeficiente);
    end;
    planos[kPlano].pon_e(NCoeficientes , -1); //coef del CF
    coef:=-coef + aElemento.e(NCoeficientes - 1);
    planos[kPlano].pon_e(NCoeficientes - 1,coef); //coef del termino independiente

    ////////////////////////////
    planos[kPlano].WriteXLTSimple_archi('aux2hp\planoInicial'+ IntToStr(kPlano)+'_'+IntToStr(kPaso) + '.xlt');
    ////////////////////////////

  end;  }


  //****ELECCION DE PLANOS INICIALES SEGUN SU CF


   //OJO! FUNCIONA SOLO PARA DOS PLANOS*************!!!!!!!!!!!!!!!

   MedianasdeElementos:= elementos.medianaFila;

   C_mediana:=MedianasdeElementos.e(NCoeficientes-1);

   for kElemento:=1 to NElementos   do // recorro todos los puntos
       if elementos.pm[kElemento].pv[elementos.nc] <= C_mediana then
         res[0].pm[kElemento]:=elementos.pm[kElemento].clonar
       else
         res[1].pm[kElemento]:=elementos.pm[kElemento].clonar;

   for kPlano :=0 to NConjuntos - 1 do
   begin
      esNil:=true;
      for kElemento:=1 to NElementos do
      begin
        esNil:= esNil and (res[kPlano].Fila(kElemento) = nil);
        if res[kPlano].Fila(kElemento) = nil then
          res[kPlano].pm[kElemento]:= TVectR.Create_InitVal(NCoeficientes-1,NO_VALIDO);
      end;
    end;

   {$IFDEF DEBUG}
      res[0].WriteArchiXLT('aux2hp\conjunto_inicial'+IntToStr(0)+'_'+IntToStr(kPaso) + '.xlt');
      res[1].WriteArchiXLT('aux2hp\conjunto_inicial'+IntToStr(1)+'_'+IntToStr(kPaso) + '.xlt');
   {$ENDIF}

   for kPlano :=0 to NConjuntos - 1 do
   begin
     plano_new:= aproxPlano(res[kPlano],invertible);
     planos[kPlano]:=plano_new.clonar;
     planos[kPlano].WriteXLTSimple_archi('aux2hp\planoInicial'+ IntToStr(kPlano)+'_'+IntToStr(kPaso) + '.xlt');
   end;

  salir:=false;

  // Calculamos la distancia de los puntos a los planos y dependiendo qué distancia sea menor
  // pertenecerá a ese conjunto plano HITLERPLANOSSS

  kIter:=0;
  repeat

  for kElemento:=1 to NElementos   do // recorro todos los puntos
  begin
    kPlanoElemento:=-1;
    kPlano_Ant:=0;
    distancia:=MaxInt;
    dist_ant:= MaxInt;
    for kPlano :=0 to NConjuntos - 1 do // recorro todos los planos
    begin
      if res[kPlano].pm[kElemento]<>nil then
      begin
       if res[kPlano].pm[kElemento].pv[1]<>NO_VALIDO then
       begin
          dist_ant:= distanciaEuclideaAPlano(elementos.pm[kElemento],planos[kPlano]);
          kPlano_Ant:=kPlano;
       end;
        res[kPlano].pm[kElemento].Free; //elimino lo que tenga el conjunto
        res[kPlano].pm[kElemento]:=nil;
      end;

      distancia_k:=distanciaEuclideaAPlano(elementos.pm[kElemento],planos[kPlano]);
      if distancia > distancia_k then
      begin
        distancia:= distancia_k;
        kPlanoElemento:=kPlano;
      end;
    end;

    if ((abs(dist_ant-distancia)) >= abs(0.2*distancia)) then   //para que cambie de conjunto, la distancia nueva tiene que ser mayor a 1.1*distancia anterior
      res[kPlanoElemento].pm[kElemento]:=elementos.pm[kElemento].clonar
    else
      begin
      res[kPlano_Ant].pm[kElemento]:=elementos.pm[kElemento].clonar
      end
  end;

  salir:=true;

  //Elegimos nuevo plano, como el mejor plano que aproxime los puntos del conjunto
  for kPlano :=0 to NConjuntos - 1 do
  begin

    planosSimilares:=true;
    esNil:=true;
    for kElemento:=1 to NElementos do
    begin
      esNil:= esNil and (res[kPlano].Fila(kElemento) = nil);
      if res[kPlano].Fila(kElemento) = nil then
        res[kPlano].pm[kElemento]:= TVectR.Create_InitVal(NCoeficientes-1,NO_VALIDO);
    end;

    if esNil then
    begin
      WriteLn('esNil');
      plano_new:=nil;//planos[kPlano]; //****OJO, NO ES LA MEJOR SOLUCION.
      {$IFDEF DEBUG}
      res[kPlano].WriteArchiXLT('aux2hp\conjunto'+IntToStr(kPlano)+'_'+IntToStr(kPaso)+'_'+IntToStr(kIter) + '.xlt');
      {$ENDIF}
    end
    else
    begin
      {$IFDEF DEBUG}
      res[kPlano].WriteArchiXLT('aux2hp\conjunto'+IntToStr(kPlano)+'_'+IntToStr(kPaso)+ '_'+ IntToStr(kIter)+ '.xlt');
      {$ENDIF}
      plano_new:= aproxPlano(res[kPlano],invertible);
     { if plano_new = nil then
       plano_new:=planos[kPlano]; //****OJO, NO ES LA MEJOR SOLUCION.  }
    end;

    if plano_new <> nil then
    begin
      for kCoordenada:=1 to plano_new.n do
          planosSimilares:=planosSimilares and (abs(plano_new.e(kCoordenada) - planos[kPlano].e(kCoordenada))<1E-3);

       if not planosSimilares then
       begin
          planos[kPlano]:= plano_new.clonar;
          if plano_repetido<>nil then
          begin
            plano_repetido.Free;
            plano_repetido:=nil;
          end;
          plano_repetido:= plano_new.clonar;
       end;

      if not invertible then
        raise Exception.Create('no hay solución, No invertible')
      else
        salir:= (salir and planosSimilares);

      {$IFDEF DEBUG}
       planos[kPlano].WriteXLTSimple_archi('aux2hp\planoNew'+ IntToStr(kPlano) +'_'+IntToStr(kPaso)+ '_'+ IntToStr(kIter)+ '.xlt');
      {$ENDIF}
    end
    else
      planos[kPlano]:=nil;
  end;

  //los planos que dieron nil, se ponen igual a otro
  for kPlano :=0 to NConjuntos - 1 do
    if planos[kPlano]=nil then
    begin
      planos[kPlano]:=TVectR.Create_Init(plano_repetido.n);
      for kCoeficiente:=1 to NCoeficientes - 2 do
      if planos[kPlano].pv[kCoeficiente]<>0 then
        planos[kPlano].pon_e(kCoeficiente,0);
      planos[kPlano].pon_e(NCoeficientes-1,aElemento.pv[elementos.nc]);
      planos[kPlano].pon_e(NCoeficientes,-1);
      {$IFDEF DEBUG}
       planos[kPlano].WriteXLTSimple_archi('aux2hp\planoNew'+ IntToStr(kPlano) +'_'+IntToStr(kPaso)+ '_'+ IntToStr(kIter)+ '.xlt');
      {$ENDIF}
    end;

  Inc(kIter);

  if plano_new<>nil then
  begin
    plano_new.Free;
    plano_new:=nil
  end;
  if plano_repetido<>nil then
  begin
    plano_repetido.Free;
    plano_repetido:=nil;
  end;

  until (salir or (kIter>NElementos));

  Result:=res;

end;

function TClasificadorKMeansPorPlanos.distanciaEuclideaVectR(punto1,
  punto2: TVectR): NReal;
var
  kElemento:integer;
begin
  Result:=0;
  for kElemento:=1 to punto1.n do
    Result:=Result + sqr(punto1.e(kElemento) - punto2.e(kElemento));
  Result:=sqrt(Result);
end;

function TClasificadorKMeansPorPlanos.distanciaEuclideaAPlano(punto, plano: TVectR
  ): NReal;
var
  kElemento,NElementos: Integer;
  numerador, denominador:NReal;
begin

  NElementos:=punto.n; // (cantidad de X) + 1 (correspondiente al Y)
  numerador:=plano.pv[NElementos]; // A0
  denominador:=0;

  {
    P=(X1*,..., Xn-1*,Y*) --> TVectR = P
    Plano ->   A0 + A1X1 + A2X2 + .....+ An-1Xn-1* + BY = 0 --> TVectR = ( A1, A2,....,An-1, A0, B) RESPETAR EL ORDEN, EL TERMINO INDEP VA PENULTIMO
    Distancia =  abs(A0 + A1X1* + A2X2* + .....BY*)/sqrt(A1^2 + .... + B^2))
  }

  for kElemento:=1 to NElementos - 1 do
  begin
    numerador:=numerador + plano.e(kElemento) * punto.e(kElemento);
    denominador:=denominador + sqr(plano.e(kElemento));
  end;

  numerador:=numerador +  plano.e(NElementos+1) * punto.e(NElementos); //BY*
  denominador:=denominador + sqr(plano.e(NElementos + 1)); //B^2
  numerador:=abs(numerador);
  denominador:=sqrt(denominador);

  Result:=numerador/denominador;

end;

function TClasificadorKMeansPorPlanos.aproxPlano(elementos: TMatR; out invertible:Boolean): TVectR;

var
  XtX, XtY: TMatR;
  X, aElemento: TVectR;
  av: NReal;
  jFil, jCol: integer;
  res: TVectR;
  resEscaler:Boolean;
  e10: longint;
  dimX, k, CantPuntosValidos, kElemento, kEstado: integer;
  desv_xi: ValReal;
  xi, yi: NReal;
  lambda_amortiguador, CFMedio, promedio, det: NReal;

begin
  invertible:=false;

  //elementos.WriteArchiXLT('elementos.xlt');

  dimX:= elementos.nc + 1 - 1; //(cant de estados + 1) + 1 por el terimo independiente y -1 para sacar el CF que esta contando en elementos[0].n
  XtX:= TMatR.Create_init( dimX, dimX );
  XtX.Ceros;
  XtY:= TMatR.Create_init( dimX, 1);
  XtY.Ceros;

  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 kElemento:= 1 to elementos.nf  do
  begin
    aElemento:= elementos.pm[kElemento];
    promedio:= aElemento.promedio;
    //***ACA PREGUNTAR SI EL DATO ES VALIDO
    if abs(promedio - NO_VALIDO) > 1E-6 then
    begin
      for kEstado:=1 to dimX - 1  do //el ultimo no, porq es el termino independiente
       X.pon_e(kEstado, aElemento.e(kEstado)); // copia dimx-1 elementos y deja el 1 al final
      X.pon_e(dimX, 1); //termino independiente

      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 * aElemento.e(aElemento.n) ); //CF
      end;
    end
    {else
      WriteLn(aElemento.e(1));  }
  end;

  //*********Ver que hacer en ese caso, por ahora son considerados outliers
  if CantPuntosValidos < dimX then
  begin
    WriteLn('CantPuntosValidos < cantidad de estados + 1, no puede resolver el sistema');
    Result:=nil;
    invertible:=true;
    Exit;
  //raise Exception.Create('CantPuntosValidos < cantidad de estados + 1, no puede resolver el sistema');
  end;

  //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');


  lambda_amortiguador:= 0.000000001;


  for k:= 1 to dimx do
    XtX.acum_e(k, k, lambda_amortiguador );


  XtX.SimetrizarTriangSup;

  XtY.WriteArchiXLT('XtY.xlt');
  XtX.WriteArchiXLT('XtX.xlt');

  det:=XtX.Escaler( XtY, resEscaler, e10 );

  res:=TVectR.Create_Init(dimx + 1);  //cantidad de estados + termino indep + CF
                                      //A0 + A1X1 + A2X2 + ..... + ANXN - Y = 0 con XN=#estados

  XtY.WriteArchiXLT('XtY_solucion.xlt');

  if resEscaler then
  begin
   for k:= 1 to dimx-1 do
      res.pon_e(k,XtY.e( k, 1 ));
   res.pon_e(dimx, XtY.e( dimx, 1 )); //termino independiente
   res.pon_e(dimx+1, -1); // -1 corresondiente a la variable Y
  end
  else
  begin
    res.Ceros;
    res.pon_e(dimx, CFMedio);
    res.pon_e(dimx+1, -1);
  end;

  invertible:=(det*power(10,e10)>1E-5) and (resEscaler) ;

  if not invertible then
     WriteLn('No invertible');

  Result:=res;

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




end;

{ TClasificador }

constructor TClasificador.Create(xElementos: TMatR; xNConjuntos: integer);
var
  kElemento:integer;
begin

  NElementos:=xElementos.nf;
  elementos:=TMatR.Create_Clone(xElementos);
  NConjuntos:=xNConjuntos;
end;

procedure TClasificador.Free;
var
  kElemento:integer;
begin
elementos.Free;
end;

end.

