unit ucluster_dem_datos;

interface

uses
  xMatDefs,
  Classes, matreal, matent, SysUtils, uAuxiliares, fddp, ufechas, Math,
  Uniform, useriestemporales, MatBool, uopencalc;

type

  { TDiaDemRec }

  TDiaDemRec = class
    // DATOS
    p: TVectR; // Potencia en p.u. de la media del día
    // Propiedades clasificadoras del dato
    PMedia: NReal; // MW-medios del día

    kAnio: word; // entero representando el año
    kMes: word; // 1..12
    kDiaDelMes: word; // 1..31
    kDiaDelAnio: integer; // 1 .. 366 ordinal del día del año
    kDiaSem: integer;
    tipoDeDia: integer; // 0: hábil, 1: Semi-Feriado, 2: Feriado;

    TempMin, TempMax, TempMed: NReal;

    // dt Fecha del dato
    dt: TDateTime;

    constructor Create(P_: TVectR; dt: TDateTime; TempHoraria: TVectR);
    constructor Create_Vacio(N: integer);
    destructor Destroy; override;


    // retorna la distancia a otra observación en el espacio de datos
    function d_Data(const a: TDiaDemRec): NReal;

    // retorna la distancia a otra observación en el espacio de Propiedades
    function d_Prop(const a: TDiaDemRec): NReal;

    procedure CopyFrom(aRec: TDiaDemRec);
    procedure WriteToODS(f: TLibroOpenCalc); virtual;
  end;

  TDAOfDiaDemRec = array of TDiaDemRec;

  { TClusterRec }

  TClusterRec = class
    centroide: TDiaDemRec;
    cnt_integrantes: integer; // cantidad de integrantes
    d_max, d_prom: NReal;
    constructor Create(NMaxIntegrantes, DimVectP: integer);
    destructor Destroy; override;
    function distancia(aDato: TDiaDemRec): NReal;
    procedure WriteToODS(f: TLibroOpenCalc);
  end;

  TClusters = array of TClusterRec;

  { TClasificador_K_Means }

  TClasificador_K_Means = class
    clusters: TClusters;
    datos: TDAofDiaDemRec;
    iCluster: TVectE;
    flg_convergio: boolean;
    cnt_iteraciones: integer;
    K_Clusters: integer;
    NMaxIteraciones: integer;

    constructor Create(datos: TDAOfDiaDemRec;
      k_clusters, NMaxIteraciones: integer);

    procedure WritToArchiODS(archi: string);

    destructor Destroy; override;

  private
    procedure InicializarCentroidesRND;
    function ReclasificarDatos: integer;
    procedure RecalcularCentroides;

    (**
    function EstimarCluster( dt: TDateTime; Temp: NReal ): integer;
    **)
  end;


function ExtraerCurvasDem(series: TSeriesDeDatos;
  kSerieDem, kSerieTemp: integer;
    kHoraIni, kHoraFin: integer // ej 0 a 23 son todas las horas
  ): TDAOfDiaDemRec;



implementation

function ExtraerCurvasDem(series: TSeriesDeDatos;
  kSerieDem, kSerieTemp: integer;
  kHoraIni, kHoraFin: integer // ej 0 a 23 son todas las horas
  ): TDAOfDiaDemRec;
var
  NCurvas: integer;
  dt_InicioDia: TDateTime;
  kHorasDiff: integer;
  kIni, kFin: integer;
  dt_MuestraSiguienteALaUltima: TDateTime;

  res: TDAOfDiaDemRec;
  sDem, sTemp: TVectR;
  Temp, PDem: TVectR;
  kDia, kDia_, kHora: integer;
  kBase: integer;
  dt_Dia: TDateTime;

begin
  if abs(series.dtEntreMuestras - 1 / 24.0) > dt_1ms / 10.0 then
    raise Exception.Create('ExtraerCurvaDem ... solo es aplicable a muestras horarias ');

  kIni := 1;
  kFin := series.NPuntos;
  dt_InicioDia := series.dtPrimeraMuestra_InicioDelDia;

  // solo sacaremos días enteros por lo cual verificamos si es necesario
  // cambiar kIni y/o kFin;
  kHorasDiff := round((series.dtPrimeraMuestra - dt_InicioDia) * 24);
  if kHorasDiff > 0 then
    kIni := (24 - kHOrasDiff) + 1
  else
    kIni := 1;

  dt_MuestraSiguienteALaUltima := series.dtPrimeraMuestra + series.NPuntos / 24.0;
  kHorasDiff := round(frac(dt_MuestraSiguienteALaUltima) * 24.0);
  kFin := series.NPuntos - kHorasDiff;

  NCurvas := Round((kFin - kIni + 1) / 24.0);
  if NCurvas <= 0 then
    raise Exception.Create('No hay datos ni para una curva');

  setlength(res, NCurvas);

  Temp := TVectR.Create_Init( 24 );

  sDem := series.series[kSerieDem];
  sTemp := series.series[kSerieTemp];

  kDia_:= 0;
  for kDia := 0 to NCurvas - 1 do
  begin
    PDem := TVectR.Create_Init(kHoraFin - kHoraIni + 1);
    kBase := kIni + kDia * 24;
    dt_Dia := series.dtPrimeraMuestra + (kBase - 1) / 24.0;
    for kHora := 0 to 23 do
    begin
      if (kHora >= kHoraIni) and (kHora <= kHoraFin ) then
        begin
         pDem.pon_e( kHora + 1 - kHoraIni, sDem.e(kBase + kHora));
        end;
      Temp.pon_e( kHora + 1, sTemp.e(kBase + kHora));
    end;
    if ( pDem.minVal > -1100 ) and ( Temp.minVal > -1100 ) then
    begin
       res[kDia_] := TDiaDemRec.Create(PDem, dt_Dia, Temp);
       inc( kDia_ );
    end;
  end;
  setlength( res, kDia_ );
  Temp.Free;
  Result := res;
end;

function TClasificador_K_Means.ReclasificarDatos: integer;
var
  kDato: integer;
  kCluster: integer;
  kMejor: integer;
  dMejor: NReal;
  aCluster: TClusterRec;
  aDato: TDiaDemRec;
  cnt_cambios: integer;
  dTest: NReal;
begin
  for kCluster := 0 to high(Clusters) do
  begin
    Clusters[kCluster].cnt_integrantes := 0;
    Clusters[kCluster].d_prom:= 0;
    Clusters[kCluster].d_max:= 0;
  end;

  cnt_cambios := 0;

  // Reclasifico y cuento cambios de clase
  for kDato := 0 to high(datos) do
  begin
    aDato := datos[kDato];
    aCluster := clusters[0];
    kMejor := 0;
    dMejor := aDato.d_Data(aCluster.centroide);
    for kCluster := 1 to high(Clusters) do
    begin
      aCluster := clusters[kCluster];
      dTest := aDato.d_Data(aCluster.centroide);
      if dTest < dMejor then
      begin
        kMejor := kCluster;
        dMejor := dTest;
      end;
    end;

    if kMejor <> iCluster.e(kDato + 1) then
    begin
      iCluster.pon_e(kDato + 1, kMejor);
      Inc(cnt_cambios);
    end;
    clusters[kMejor].d_prom:= clusters[kMejor].d_prom + dMejor;
    if dMejor > clusters[kMejor].d_max then
      clusters[kMejor].d_max:= dMejor;
    Inc(clusters[kMejor].cnt_integrantes);
  end;

writeln( 'cnt_iters: ', cnt_iteraciones, ' : cambios: ', cnt_cambios );
  if cnt_cambios > 0 then
    RecalcularCentroides;
  Result := cnt_cambios;
end;

procedure TClasificador_K_Means.RecalcularCentroides;
var
  kCluster, kDato: integer;
  aCluster: TClusterRec;
begin
  // limpiamos los acumuladores
  for kCluster := 0 to high(Clusters) do
    Clusters[kCluster].centroide.p.Ceros;

  // recorremos los datos y los sumamos al que corresponda
  for kDato := 0 to high(datos) do
  begin
    kCluster := iCluster.e(kDato + 1);
    clusters[kCluster].centroide.p.sum(datos[kDato].p);
  end;

  // ahora promediamos
  for kCluster := 0 to high(Clusters) do
  begin
    aCluster := Clusters[kCluster];
    aCluster.centroide.p.PorReal(1 / aCluster.cnt_integrantes);
    aCluster.d_prom:= aCluster.d_prom / aCluster.cnt_integrantes;
  end;

end;

(**


function TClasificador_K_Means.EstimarCluster(dt: TDateTime; Temp: NReal
  ): integer;

var
  f: TFecha;
  kDia, kDiaSem, tipoDeDia: integer;
  Candidatos: TList;
  aDato: TDiaDemRec;

function DistanciaCandidato( aDato: TDiaDemRec ): NReal;
var
  dDias: NReal;
  dTemp: NReal;

begin
  if aDato.tipoDeDia <> tipoDeDia then
    result:= -1
  else


end;

begin
  f := TFecha.Create_Dt(dt);
  kDia := f.getDiaDelAnio;
  kDiaSem:= diasem( dt );
  tipoDeDia := Ord(TipoDeDiaUruguay(dt));
  f.Free;

  Candidatos:= TList.Create;

  // filtrado de los candidatos
  for kDato:= 0 to high( datos ) do
  begin
    aDato:= datos[kDato];
    if aDato.tipoDeDia = tipoDeDia then
      Candidatos.add( aDato );
  end;
  Candidatos.Sort( Compara_kDia );


end;

   **)
procedure TClasificador_K_Means.InicializarCentroidesRND;
var
  flg_Seleccionado: TVectBool;
  kCluster, kDato: integer;
  NDatos: integer;
  unif: uniform.T_UNIFORM_RND_GENERATOR;
  buscando: boolean;

begin
  NDatos := length(datos);
  flg_Seleccionado := TVectBool.Create_Init(NDatos);
  flg_Seleccionado.Ceros;
  unif := T_UNIFORM_RND_GENERATOR.Create( 31 );
  // Para inicializar cargamos al azar un centroide en cada cluster
  // cuidando de que no se repitan
  for kCluster := 0 to high(clusters) do
  begin
    kDato := trunc(unif.Call_UNI * NDatos);
    buscando := True;
    while buscando do
    begin
      if kDato >= NDatos then
        kDato := 0;
      if flg_Seleccionado.e(kDato + 1) then
        Inc(kDato)
      else
      begin
        buscando := False;
        Clusters[kCluster].centroide.CopyFrom(datos[kDato]);
        flg_Seleccionado.pon_e(kDato + 1, True);
      end;
    end;
  end;
  flg_Seleccionado.Free;
  unif.Free;
end;



{ TClasificador_K_Means }

constructor TClasificador_K_Means.Create(datos: TDAOfDiaDemRec;
  k_clusters, NMaxIteraciones: integer);
var
  // Conjunto de vectores para almacenar los centroides y los vectores asociados
  NDatos: integer;
  kDato, kCluster: integer;

begin
  inherited Create;
  self.datos := datos;
  self.K_Clusters := k_CLusters;
  self.NMaxIteraciones := NMaxIteraciones;
  flg_convergio := False;

  NDatos := length(datos);
  if k_Clusters > NDatos then
    raise Exception.Create('No tiene sentdio k_Clusters > NDatos ');

  setlength(Clusters, k_clusters);
  for kCluster := 0 to k_Clusters - 1 do
    Clusters[kCluster] := TClusterRec.Create(NDatos, datos[0].p.n);

  // Iniciamos con centroides elegidos al azar
  InicializarCentroidesRND;

  // Inicializamos el vector de asignación en -1
  iCluster := TVectE.Create_Init(NDatos);
  for kDato := 1 to length(datos) do
    iCluster.pon_e(kDato, -1); // no clasificado

  cnt_iteraciones := 0;
  flg_convergio := False;
  while (not flg_convergio) and (cnt_iteraciones < NMaxIteraciones) do
  begin
    if ReclasificarDatos > 0 then
      Inc(cnt_iteraciones)
    else
      flg_convergio := True;
  end;
end;

procedure TClasificador_K_Means.WritToArchiODS(archi: string);
var
  f: TLibroOpenCalc;
  kCluster, kDato: integer;
  aDato: TDiaDemRec;
  aCluster: TClusterRec;
begin
  f := TLibroOpenCalc.Create(True, archi);
  f.write('NIteraciones:' );
  f.write( cnt_iteraciones );
  f.write( 'NMaxIteraciones:' );
  f.Write( NMaxIteraciones );
  f.writeln;

  f.writeln('Clusters: ');
  for kCluster := 0 to high(clusters) do
  begin
    aCluster := clusters[kCluster];
    f.Write(kCluster);
    aCluster.writeToOds(f);
    f.writeln;
  end;

  f.writeln;
  f.writeln('Datos:');
  for kDato := 0 to high(datos) do
  begin
    aDato := datos[kDato];
    f.Write(kDato);
    f.Write(iCluster.e(kDato + 1));
    aDato.WriteToODS(F);
    f.writeln;
  end;
  f.Free;

end;

destructor TClasificador_K_Means.Destroy;
var
  kDato, kCluster: integer;
begin
  for kDato := 0 to high(datos) do
    datos[kDato].Free;
  setlength(datos, 0);

  for kCluster := 0 to high(clusters) do
    clusters[kCluster].Free;
  setlength(clusters, 0);
end;


{ TClusterRec }

constructor TClusterRec.Create(NMaxIntegrantes, DimVectP: integer);
begin
  inherited Create;
  centroide := TDiaDemRec.Create_Vacio( DimVectP );
  cnt_integrantes := 0;
  d_max:= 0;
  d_prom:= 0;
end;

destructor TClusterRec.Destroy;
begin
  centroide.Free;
  inherited Destroy;
end;

function TClusterRec.distancia(aDato: TDiaDemRec): NReal;
begin
  Result := centroide.d_Data(aDato);
end;


procedure TClusterRec.WriteToODS(f: TLibroOpenCalc);
begin
  f.Write(cnt_integrantes);
  f.write( d_max );
  f.write( d_prom );
  centroide.WriteToODS(f);
end;

{ TDiaDemRec }

constructor TDiaDemRec.Create(P_: TVectR; dt: TDateTime; TempHoraria: TVectR);
var
  f: TFecha;
  kHora: integer;
  aTemp: NReal;
begin
  inherited Create;
  p := P_;
  self.dt := dt;

  PMedia := p.promedio;
  p.PorReal(1 / PMedia);

  f := TFecha.Create_Dt(dt);
  DecodeDate( dt, kAnio, kMes, kDiaDelMes );
  kDiaDelAnio := f.getDiaDelAnio;
  kDiaSem:= diasem( dt );
  tipoDeDia := Ord(TipoDeDiaUruguay(dt));

  TempMin := TempHoraria.e(1);
  TempMax := TempMin;
  TempMed := TempMin;
  for kHora := 2 to 24 do
  begin
    aTemp := TempHoraria.e(kHora);
    if aTemp < TempMin then
      TempMin := aTemp
    else if aTemp > TempMax then
      TempMax := aTemp;
    TempMed := TempMed + aTemp;
  end;
  TempMed := TempMed / 24.0;
end;

constructor TDiaDemRec.Create_Vacio( N: integer );
begin
  inherited Create;
  p := TVectR.Create_Init( N );
  PMedia := 0;
  kAnio :=1900;
  kMes:= 1;
  kDiaDelMes:= 1;
  kDiaDelAnio := 1;
  tipoDeDia := 0;
  TempMin := 0;
  TempMax := 0;
  TempMed := 0;
end;


destructor TDiaDemRec.Destroy;
begin
  p.Free;
  inherited Destroy;
end;

function TDiaDemRec.d_Data(const a: TDiaDemRec): NReal;
begin
  Result := p.distancia(a.p);
end;

function TDiaDemRec.d_Prop(const a: TDiaDemRec): NReal;
var
  res: NReal;

begin
  res := abs(PMedia - a.PMedia);
  res := res + DistanciaCircular( kDiaDelAnio - 1, a.kDiaDelAnio - 1, 366);
  res := res + DistanciaCircular( kDiaSem - 1, a.kDiaSem - 1, 7);
  res := res + abs(tipoDeDia - a.tipoDeDia);
  res := res + abs(TempMin - a.TempMin);
  res := res + abs(TempMax - a.TempMax);
  res := res + abs(TempMed - a.TempMed);
  Result := res;
end;

procedure TDiaDemRec.CopyFrom(aRec: TDiaDemRec);
begin
  p.Copy(aRec.p);
  PMedia := aRec.PMedia;
  kAnio:= aRec.kAnio;
  kMes:= aRec.kMes;
  kDiaDelMes:= aRec.kDiaDelMes;
  kDiaDelAnio := aRec.kDiaDelAnio;
  tipoDeDia := aRec.tipoDeDia;
  TempMin := aRec.TempMin;
  TempMax := aRec.TempMax;
  TempMed := aRec.TempMed;
end;

procedure TDiaDemRec.WriteToODS(f: TLibroOpenCalc);
var
  k: integer;

begin
  f.Write(dt);
  f.Write(PMedia);
  f.write( kAnio );
  f.write( kMes );
  f.write( kDiaDelMes );
  f.Write(kDiaDelAnio);
  f.Write(kDiaSem);
  f.Write(tipoDeDia);
  f.Write(TempMin);
  f.Write(TempMax);
  f.Write(TempMed);

  for k := 1 to p.n do
    f.Write(p.e(k));

end;

end.
