unit ucalibrarconoensamble;

{$mode delphi}

interface

uses
  Classes,
  FileUtil,
  fpsTypes,
  Math,
  matreal,
  uauxiliares,
  ucalibrarconopronosticos,
  ucosa,
//  ucosaConNombre,
//  ucosanubeseable,
//  uDataSetGenerico,
  uensembles,
  ufechas,
  uFuenteSintetizador,
  ufxgx,
  {$IFDEF GRAFICAR_ENSEMBLES}
  ugraficador,
  {$ENDIF}
  uInicioYFinal,
  umodelosintcegh,
  uopencalc,
//  uparseutils,
  uproblema,
  upronostico,
//  uresfxgx,
  usalasdejuego,
  useriestemporales,
  {$IFDEF PRONO_HIDRO_IMFIA}
  udbconpg,
  uoddface_roscon_cfg,
  uestadosinuy,
  {$ENDIF}
  SysUtils,
  xmatdefs;


CONST
  CARPETA_RESULTADOS_HIDRO = 'rchaer:chinchulin@192.168.0.35:/home/rchaer/prono_hidro_imfia/ensambles';

  NCAPA_ = 0;
  NESCENARIO_ = 0;
  NPRONOSTICOS_ = 0;
  NPRONOSTICO_ = 0;

  NCANALES_ = 6;
  DUR_PASO_ = 1/24;
  SUMA_0 = 1e-70;
  BETA_0 = 1e+20;
  BETA_MIN = 0.01;
  BETA_MAX = 1e+50;
  BETA_INI = 100;
  INVERSO_DE_CERO = 1e+60;



{GFlieller 29-11-2018
Ajustar el cono de un CEGH a ensamble de pronosticos
1- Buscar las matrices A y B en el CEGH. Calcular la matriz C=(B B^T)^-1
2- Transformar el ensamble a espacio gaussiano
3- Calcular los sesgos
4- Calcular los desvios
5- Calcular la matriz de covarianzas
6- Calcular la matriz de ajuste con el paso anterior "Z"
7- Calcular la matriz M(covar,Z)
8- Resolver MX = "1/X"
9- Volver al paso 6
}

type

    { TcalibrarConoEnsamble }

    TcalibrarConoEnsamble = class
      //  Filtros A y B.
      A : TMatR;
      B : TMatR;
      //  C = (B B^T)-1 .
      C : TMatR;

      //  Fecha inicio del cono.
      fechaIni : TFecha;

      //  CEGH que se va a calibrar.
      fcegh : TFuenteSintetizadorCEGH;

      //  Ensembles de pronosticos en Espacio Real y en Espacio Gausssiano.
      Ensambles : array of TMatR;
      EnsamblesGauss : array of TMatR;

      //  Xg es una representacion del los ensembles en array (por miembro) de
      //  array (por paso de tiempo) de matrices de Ncanales x 1 para poder
      //  multiplicarlos como multiplicacion de matrices.
      Xg : array of array of TMatR;
      //  Xesp es el valor esperado de Xg representado como array (por paso de
      //  tiempo), con el promedio de todos los miembros de cada canal, como
      //  matriz de Ncanales x 1.
      Xesp : array of TMatR;

      //  Numero de canales de los ensembles, miembros de cada ensemble y pasos
      //  de tiempo de cada ensemble. Se asume que todos los ensembles tienen el
      //  el mismo numero de miembros y que cada miembro tiene el mismo numero
      // pasos de tiempo.
      Ncanales : integer; // Series [Bonete,Palmar,Salto]
      Nmiembros : integer;
      Npasos : integer;

      //  Numero de pasos de optimizacion
      pasosTotales : integer;

      //  Sesgos calculados. Array (por paso de tiempo) de matrices Ncanales x 1.
      S : array of TMatR;

      //  Desvios por miembro y por paso de tiempo respresntados en matrices
      //  de Ncanales x 1, para crear las matrices de covarianzas de
      //  Ncanales x Ncanales
      des : array of array of TMatR;
      covar : array of TMatR;

      //  Matriz auxiliar para crear M = (d*d^T - Z). C
      //  Z = A *( Z0 + L * B * B^T * L^T )* A^T
      Z : array of TMatR;

      //  Matriz M para resolver el problema (M x ('1/Lambda') ) = Lambda
      M : array of TMatR;
      Mreducida : TMatR;

      //  Matriz de Betas calculados para el paso de tiempo 'actual' como matriz
      //  de Ncalaes x 1. Lambdas del paso en el mismo formato.
      Beta : TMatR;
      Lambda : TMatR;

      //  Registro de los betas y lambdas de todos los pasos de tiempo en el
      //  formato de 'Beta' y 'Lambda'. Para cargar al CEGH se pasa Lambdas de
      //  formato array de matrices a array de TDA of NReal
      betas : array of TVectR;
      Lambdas : array of TMatR;

      //  Factores y sesgos que se cargan al TPronostico del TPronosticos del
      //  escenario de pronosticos del CEGH que se quiere modificar
      factores : array of TDAofNReal;
      sesgos: array of TDAofNReal;

      //  El escenario de pronosticos del CEGH se arma con TPronosticos,
      //  compuesto por varios TPronostico. Se cargan los datos de factores y
      //  sesgos a aPronostico que se carga a aPronosticos que reemplaza el
      //  primer elemento del escenario de pronosticos del cegh.
      aPronosticos: TPronosticos;
      aPronostico :array of TPronostico;


      constructor Create( acegh: TFuenteSintetizadorCEGH;
                          aEnsambles: array of TMatR; aFecha : TFecha);
      //  Crea el problema en base al CEGH que se va a calibrar y el ensembles
      //  que se va a usar. Aca se leen los filtros, se lee el cegh, se
      //  gaussianiza el ensemble, se calculan las guias, se calculan los
      //  sesgos, se calcula las dseviaciones y las covarianzas.

      function calibrarConoEnsamble() : boolean;
      //  Loop de resolucion de los factores por paso.

      procedure cargarDimensiones(aEnsambles: array of TMatR);
      // carga las dimensiones: Ncanales, Nmiembros, Nensambles

      procedure cargarEnsamblesGausiano(aEnsambles: array of TMatR);
      // carga todos los ensambles a la clase

      procedure calcularGuias;
      // Calcula las guias en mundo gaussiano y real

      procedure calcularSesgos;
      // Calcula los sesgos S_k

      procedure calcularCovar;
      // Con los desvios calculados y cargados es des[][] calcula las matrices de covarianza

      function calcularZ(Z0 : TMatR; kpaso : integer) : TMatR;
      // Calcula la matriz que en la doc se denota como Z. Corresponde a un
      // ajuste de la varianza entre pasos

      procedure cargarFiltros;
      // Carga los filtros A y B del modelo. Ademas calcula la matriz C = ( B x B^{T} )^{-1}

      procedure armarProblema(Z : TMatR; kpaso : integer);
      // Arma la matriz M del problema M * X = 1/X

      procedure ImprimirLambdas;
      // Guarda las matrices de lambdas

      procedure guardarSala(sala: TSalaDeJuego);

      function TransformarEnsamble(aModelo:TModeloCEGH; aEnsambles:array of TMatR; hcanal: integer) : TMatR;
      // Trasforma el ensamble a espacio gaussiano

      private procedure ArmarResultados();
      // Carga al modelo CEGH los valores de guia, guia gaussina, sesgos,
      // factores, NPAC, NPSA, NPCC

      procedure calcularDesvios;
      // En Tpronosticos, que es la clase que utiliza el Sintetizador, la guia,
      // los sesgos y los factores son TDAofNReal, pero para manipular estos
      // valores en la unidad se usan arrays of TMatR. Acá lo que se hace es
      // pasar los valores que se tienen al formato que requiere el Sintetizador.

      procedure free;

      procedure publicarEnsemble;

      function ReducirProblema(MatrizOrig: TMatR;
                  var indices :TDAofNInt; var L : TMatR;
                  var MatrizNula :Boolean): TMatR;
{ Toma la matriz M del problema, calcula los factores que son tribiales,
indeterminados o solo se pueden resolver analiticamente. Elimina las variables
calculadas de la matriz. Devuelve la matriz reducida, un indice con el lugar de
las variables calculadas, un vector con las variables calculadas y un flag que
indica si la matriz reducida es de tamano nulo.
}

end;

type

  { T_fobjetivoCEGH }

  T_fobjetivoCEGH = class(Tfx)


  Error: NReal;
  M : TMatR;
  FX : TVectR;
  Ncanales : integer;

  constructor Create(aM: TMatR);
  // Crea la funcion objetivo usando la matriz M del problema

  function f(const X: TVectR): NReal; override;
  // Calculo de costo.

  destructor destroy; override;

  procedure acum_g(var grad: TVectR; const X: TVectR); override;
  // calculo de gradiente. El gradiente depende de FX, X y M
  end;

type

    { TProblemaResolverFactores }
    TProblemaResolverFactores = class(TProblema_m01)

    M: TMatR;                 // Matriz del problema
    Betas: TVectR;
    Lambdas: TMatR;
    Ncanales:integer;

    constructor Create(aM:TMatR);
    // Crea el problema M * X = 1/X

    procedure free;

    function resolverFactores(var factores:TMatR): boolean;
    // Resuelve el problema M * X = 1/X

    end;


//  Funciones sueltas
procedure WriteToArchiEnsembles(EnsembleProno,EnsembleCEGH:TEnsembles;
            nombreArchi:string; formato:TsSpreadsheetFormat=sfOpenDocument);

procedure WriteToArchiODSResultdosSim(DatosEnsamble,DatosSimCEGH:array of TMatR;
          dtIni, paso : TDateTime; archi: string;
          formato:TsSpreadsheetFormat=sfOpenDocument);

procedure WriteToArchiXLTResultdosSim(DatosEnsamble,DatosSimCEGH:array of TMatR;
          archi : string);



{$IFDEF PRONO_HIDRO_IMFIA}
procedure cargarEnsamblesLocales(var SeriesDeEnsambles : array of TMatR;archidir: string;Nseries: integer);
          // Carga los ramilletes de pronosticos del IMIFIA a partir de archivo .xlt
          // guardados localmente. Se usa para debuguear el programa.

procedure cargarEnsamblesDesdeEstado(var SeriesDeEnsambles : array of TMatR;estado:TEstado_SIN_UY;Nseries: integer);
          // Carga los ramilletes de pronosticos del IMFIA a partir del estado. En el
          // estado debe estar precargados los ramilletes a traves de consulta sql y
          // corrida del modelo del IMFIA.


procedure graficarResultados;


function getEnsemblesAportesVates(dirSalida: string;
                                   dbpqPronos: TDBPQCon; estado :TEstado_SIN_UY) : TEnsembles;
{$ENDIF}




implementation

procedure WriteToArchiXLTResultdosSim(DatosEnsamble,
  DatosSimCEGH: array of TMatR; archi: string);
var
  SeriesPE_CEGH, SeriesPE_IMFIA : array of TMatR;               // Series de probabilidad de excedencia
  PE5_CEGH,PE10_CEGH,PE20_CEGH,PE50_CEGH,PE80_CEGH,PE90_CEGH,PE95_CEGH, PE5_IMFIA,PE10_IMFIA,PE20_IMFIA,PE50_IMFIA,PE80_IMFIA,PE90_IMFIA,PE95_IMFIA : array of TVectR;   // Cortes de probabilidad
  kSerie: Integer;
begin
  // 5.1  Paso los ensambles a .xlt
  DatosEnsamble[0].WriteArchiXLT('/home/'+GetEnvironmentVariable('USER')+'/prono_hidro_imfia/Resultados/Bonete/Ensamble/'+archi+'.xlt');
  if length(DatosEnsamble)>1 then
    DatosEnsamble[1].WriteArchiXLT('/home/'+GetEnvironmentVariable('USER')+'/prono_hidro_imfia/Resultados/Palmar/Ensamble/'+archi+'.xlt');
  if length(DatosEnsamble)>2 then
    DatosEnsamble[2].WriteArchiXLT('/home/'+GetEnvironmentVariable('USER')+'/prono_hidro_imfia/Resultados/SG/Ensamble/'+archi+'.xlt');

  // 5.2  Paso los ramilletes de sinteticas a .xlt
  DatosSimCEGH[0].WriteArchiXLT('/home/'+GetEnvironmentVariable('USER')+'/prono_hidro_imfia/Resultados/Bonete/Sinteticas/'+archi+'.xlt');
  if length(DatosSimCEGH)>1 then
    DatosSimCEGH[1].WriteArchiXLT('/home/'+GetEnvironmentVariable('USER')+'/prono_hidro_imfia/Resultados/Palmar/Sinteticas/'+archi+'.xlt');
  if length(DatosSimCEGH)>2 then
    DatosSimCEGH[2].WriteArchiXLT('/home/'+GetEnvironmentVariable('USER')+'/prono_hidro_imfia/Resultados/SG/Sinteticas/'+archi+'.xlt');

  for kSerie := 0 to (length(DatosSimCEGH)-1) do
  begin
  // 5.3  Ordeno las series para sacar cada corte de probabilidad
  DatosSimCEGH[kSerie].Sort();
  // 5.4  Cargo cada corte de porbabilidad a la serie correspondiente
  PE5_CEGH[kSerie]:= DatosSimCEGH[kSerie].Columna(round(0.05*DatosSimCEGH[kSerie].nc));
  PE10_CEGH[kSerie]:=DatosSimCEGH[kSerie].Columna(round(0.10*DatosSimCEGH[kSerie].nc));
  PE20_CEGH[kSerie]:=DatosSimCEGH[kSerie].Columna(round(0.20*DatosSimCEGH[kSerie].nc));
  PE50_CEGH[kSerie]:=DatosSimCEGH[kSerie].Columna(round(0.50*DatosSimCEGH[kSerie].nc));
  PE80_CEGH[kSerie]:=DatosSimCEGH[kSerie].Columna(round(0.80*DatosSimCEGH[kSerie].nc));
  PE90_CEGH[kSerie]:=DatosSimCEGH[kSerie].Columna(round(0.90*DatosSimCEGH[kSerie].nc));
  PE95_CEGH[kSerie]:=DatosSimCEGH[kSerie].Columna(round(0.95*DatosSimCEGH[kSerie].nc));
  // 5.5  Creo matrices con todos los cortes de probabilidad para cada canal
  SeriesPE_CEGH[kSerie]:=TMatR.Create_Init(336,7);
  SeriesPE_CEGH[kSerie].pon_columna(1,1,PE5_CEGH[kSerie].toTDAOfNReal());
  SeriesPE_CEGH[kSerie].pon_columna(1,2,PE10_CEGH[kSerie].toTDAOfNReal());
  SeriesPE_CEGH[kSerie].pon_columna(1,3,PE20_CEGH[kSerie].toTDAOfNReal());
  SeriesPE_CEGH[kSerie].pon_columna(1,4,PE50_CEGH[kSerie].toTDAOfNReal());
  SeriesPE_CEGH[kSerie].pon_columna(1,5,PE80_CEGH[kSerie].toTDAOfNReal());
  SeriesPE_CEGH[kSerie].pon_columna(1,6,PE90_CEGH[kSerie].toTDAOfNReal());
  SeriesPE_CEGH[kSerie].pon_columna(1,7,PE95_CEGH[kSerie].toTDAOfNReal());
  end;

  SeriesPE_CEGH[0].WriteArchiXLT('/home/'+GetEnvironmentVariable('USER')+'/prono_hidro_imfia/Resultados/Bonete/Sinteticas_PE/'+archi+'.xlt');
  if length(DatosEnsamble)>1 then
    SeriesPE_CEGH[1].WriteArchiXLT('/home/'+GetEnvironmentVariable('USER')+'/prono_hidro_imfia/Resultados/Palmar/Sinteticas_PE/'+archi+'.xlt');
  if length(DatosEnsamble)>2 then
    SeriesPE_CEGH[2].WriteArchiXLT('/home/'+GetEnvironmentVariable('USER')+'/prono_hidro_imfia/Resultados/SG/Sinteticas_PE/'+archi+'.xlt');


  for kSerie := 0 to (length(DatosEnsamble)-1) do
    begin
    // 5.5  Ordeno las series para sacar cada corte de probabilidad
    DatosEnsamble[kSerie].Sort();
    // 5.6  Cargo cada corte de porbabilidad a la serie correspondiente
    PE5_IMFIA[kSerie] := DatosEnsamble[kSerie].Columna(round(0.05*DatosEnsamble[kSerie].nc));
    PE10_IMFIA[kSerie]:=DatosEnsamble[kSerie].Columna(round(0.10*DatosEnsamble[kSerie].nc));
    PE20_IMFIA[kSerie]:=DatosEnsamble[kSerie].Columna(round(0.20*DatosEnsamble[kSerie].nc));
    PE50_IMFIA[kSerie]:=DatosEnsamble[kSerie].Columna(round(0.50*DatosEnsamble[kSerie].nc));
    PE80_IMFIA[kSerie]:=DatosEnsamble[kSerie].Columna(round(0.80*DatosEnsamble[kSerie].nc));
    PE90_IMFIA[kSerie]:=DatosEnsamble[kSerie].Columna(round(0.90*DatosEnsamble[kSerie].nc));
    PE95_IMFIA[kSerie]:=DatosEnsamble[kSerie].Columna(round(0.95*DatosEnsamble[kSerie].nc));
    // 5.7  Creo matrices con todos los cortes de probabilidad para cada canal
    SeriesPE_IMFIA[kSerie]:=TMatR.Create_Init(336,7);
    SeriesPE_IMFIA[kSerie].pon_columna(1,1,PE5_IMFIA[kSerie].toTDAOfNReal());
    SeriesPE_IMFIA[kSerie].pon_columna(1,2,PE10_IMFIA[kSerie].toTDAOfNReal());
    SeriesPE_IMFIA[kSerie].pon_columna(1,3,PE20_IMFIA[kSerie].toTDAOfNReal());
    SeriesPE_IMFIA[kSerie].pon_columna(1,4,PE50_IMFIA[kSerie].toTDAOfNReal());
    SeriesPE_IMFIA[kSerie].pon_columna(1,5,PE80_IMFIA[kSerie].toTDAOfNReal());
    SeriesPE_IMFIA[kSerie].pon_columna(1,6,PE90_IMFIA[kSerie].toTDAOfNReal());
    SeriesPE_IMFIA[kSerie].pon_columna(1,7,PE95_IMFIA[kSerie].toTDAOfNReal());
    end;


  SeriesPE_IMFIA[0].WriteArchiXLT('/home/'+GetEnvironmentVariable('USER')+'/prono_hidro_imfia/Resultados/Bonete/Ensambles_PE/'+archi+'.xlt');
  if length(SeriesPE_IMFIA)>1 then
    SeriesPE_IMFIA[1].WriteArchiXLT('/home/'+GetEnvironmentVariable('USER')+'/prono_hidro_imfia/Resultados/Palmar/Ensambles_PE/'+archi+'.xlt');
  if length(SeriesPE_IMFIA)>2 then
    SeriesPE_IMFIA[2].WriteArchiXLT('/home/'+GetEnvironmentVariable('USER')+'/prono_hidro_imfia/Resultados/SG/Ensambles_PE/'+archi+'.xlt');

  // 5.8  Libero memoria
  for kSerie := 0 to (length(DatosSimCEGH)-1) do
  begin
    PE5_CEGH[kserie].free;
    PE10_CEGH[kserie].free;
    PE20_CEGH[kserie].free;
    PE50_CEGH[kserie].free;
    PE80_CEGH[kserie].free;
    PE90_CEGH[kserie].free;
    PE95_CEGH[kserie].free;
    SeriesPE_CEGH[kSerie].free;
  end;
  for kSerie := 0 to (length(DatosEnsamble)-1) do
  begin
    PE5_IMFIA[kserie].free;
    PE10_IMFIA[kserie].free;
    PE20_IMFIA[kserie].free;
    PE50_IMFIA[kserie].free;
    PE80_IMFIA[kserie].free;
    PE90_IMFIA[kserie].free;
    PE95_IMFIA[kserie].free;
    SeriesPE_IMFIA[kSerie].free;
  end;

  SetLength(PE5_CEGH ,0);
  SetLength(PE10_CEGH,0);
  SetLength(PE20_CEGH,0);
  SetLength(PE50_CEGH,0);
  SetLength(PE80_CEGH,0);
  SetLength(PE90_CEGH,0);
  SetLength(PE95_CEGH,0);
  SetLength(SeriesPE_CEGH,0);
  SetLength(PE5_IMFIA,0);
  SetLength(PE10_IMFIA,0);
  SetLength(PE20_IMFIA,0);
  SetLength(PE50_IMFIA,0);
  SetLength(PE80_IMFIA,0);
  SetLength(PE90_IMFIA,0);
  SetLength(PE95_IMFIA,0);
  SetLength(SeriesPE_IMFIA,0);



end;


procedure WriteToArchiEnsembles(EnsembleProno,EnsembleCEGH:TEnsembles; nombreArchi:string; formato:TsSpreadsheetFormat=sfOpenDocument);
// Este procedimiento guarda los cortes de probabilidad de los resultyado de la
// simulacion del cegh y los ensembles de pronosticos en el archivo
// nombreArchi. nombreArchi debe incluir la extension xlsx o ods.
{
1.  Creo el archivo ods o xlsx.
2.  Creo cada una de las hojas para los pronosticos. Una hoja para cada canal
    ensemble completo con todas las cronicas
3.  Creo cada una de las hojas para los ensembles de cegh. Una hoja para cada
    canal, solo incluyo los cortes de probabilidad PE5, PE10, PE20, PE50, PE80,
    PE90 y PE95.
}


var
  ods: TLibroOpenCalc;
  nombresPE : array[0..6] of string =
                             ('PE5','PE10','PE20','PE50','PE80','PE90','PE95');
  PE : array[0..6] of NReal = ( 0.05, 0.1, 0.2, 0.5, 0.8, 0.9, 0.95 );
  jmiembro, hcanal, kpaso, colPE: Integer;
  NombreHoja: String;
  matAux: TMatR;
  columna: TDAofNReal;


begin

  //  1.
  if (formato = TsSpreadsheetFormat.sfOpenDocument) or (formato = TsSpreadsheetFormat.sfOOXML) then
     ods := TLibroOpenCalc.Create(TRUE,nombreArchi,False,formato)
  else
     writeln('El formato de los resultados debe ser .xlsx o .ods');

  //  2.
  for hcanal := 0 to high(EnsembleProno.nombresEnsembles) do
  begin
      //  2.1  Crear hoja
      NombreHoja := EnsembleProno.nombresEnsembles[hcanal];
      ods.agregoHoja(NombreHoja);
      //  2.2  Escribir el encabezado
      ods.go(NombreHoja,1,1);
      ods.write('Paso');
      for jmiembro := 1 to EnsembleProno.Nmiembros do
        ods.write('member_'+IntToStr(jmiembro));
      //  2.3  Escribir el paso de tiempo
      for kpaso := 1 to EnsembleProno.Npasos do
        begin
        ods.go(NombreHoja,1+kpaso,1);
        ods.Write( EnsembleProno.dtPrimeraMuestra + (kpaso-1)*EnsembleProno.dtPaso );
        end;
      //  2.4  Escribir las matrices del ensemble
      ods.Write_MatR(EnsembleProno.EnsemblesMat[hcanal],NombreHoja,2,2);
  end;

  //  3.
  for hcanal := 0 to high(EnsembleCEGH.nombresEnsembles) do
  begin
      //  3.1  Crear la hoja
      NombreHoja := EnsembleCEGH.nombresEnsembles[hcanal]+'_CEGH';
      ods.agregoHoja(NombreHoja);
      //  3.2  Escribir el encabezado
      ods.go(NombreHoja,1,1);
      ods.write('Paso');
      for jmiembro := 0 to high(nombresPE) do
        ods.write(nombresPE[jmiembro]);
      //  3.3  Escribir el paso de tiempo
      for kpaso := 1 to EnsembleCEGH.Npasos do
        begin
        ods.go(NombreHoja,1+kpaso,1);
        ods.Write( EnsembleCEGH.dtPrimeraMuestra + (kpaso-1)*EnsembleCEGH.dtPaso );
        end;

      //  3.4  Ordenar el ensemble
      EnsembleCEGH.EnsemblesMat[hcanal].Sort();
      //  3.5  Escribir las columnas correspondientes a los PE
      for jmiembro := 0 to high(nombresPE) do
      begin
        matAux := TMatR.Create_Init(EnsembleCEGH.EnsemblesMat[hcanal].nf,1);
        colPE := round( PE[jmiembro]*EnsembleCEGH.EnsemblesMat[hcanal].nc );
        columna := EnsembleCEGH.EnsemblesMat[hcanal].Columna(colPE).toTDAOfNReal();
        matAux.pon_columna( 1, 1, columna);

        ods.Write_MatR(matAux,NombreHoja,2,jmiembro+2);
        matAux.free;

      end;

  end;

  ods.guardar;
  ods.free;


end;



procedure WriteToArchiODSResultdosSim(DatosEnsamble,DatosSimCEGH: array of TMatR; dtIni, paso : TDateTime; archi: string; formato:TsSpreadsheetFormat=sfOpenDocument);
// 5. Guardo en un ods los conos del ensamble y los conos creados con el CEGH
//    como cortes de probabilidad
var
ods: TLibroOpenCalc;
Nombres : array of string;
NombreHoja : string;
SeriesPE_CEGH, SeriesPE_IMFIA : array of TMatR;               // Series de probabilidad de excedencia
PE5_CEGH,PE10_CEGH,PE20_CEGH,PE50_CEGH,PE80_CEGH,PE90_CEGH,PE95_CEGH, PE5_IMFIA,PE10_IMFIA,PE20_IMFIA,PE50_IMFIA,PE80_IMFIA,PE90_IMFIA,PE95_IMFIA : array of TVectR;   // Cortes de probabilidad
kSerie, NcanalesCEGH, hcanal, kpaso, jmiembro: Integer;

// GFlieller@2019-05-23:
// A partir de ahora no guardo mas todas las cronicas del CEGH, sino solo los
// cortes de probabilidad.

begin

// 5.0  Cargo los nombres
//NcanalesIMFIA := length(DatosEnsamble);

NcanalesCEGH := trunc(length(DatosSimCEGH));
SetLength(Nombres,NcanalesCEGH);
Nombres[0] := 'Bonete';
  if NcanalesCEGH > 1 then
     Nombres[1] := 'Palmar';
  if NcanalesCEGH > 2 then
     Nombres[2] := 'SG';
  if NcanalesCEGH > 3 then
     Nombres[3] := 'in34';
  if NcanalesCEGH > 4 then
     Nombres[4] := 'cmo1';
  if NcanalesCEGH > 5 then
     Nombres[5] := 'cmo3';

// 5.1  Creo el libro ods que va contener todas las hojas
if (formato = TsSpreadsheetFormat.sfOpenDocument) or (formato = TsSpreadsheetFormat.sfOOXML) then
   ods := TLibroOpenCalc.Create(TRUE,archi,False,formato)
else
   writeln('El formato de los resultados debe ser .xlsx o .ods');

// 5.2  Agrego cada hoja de acuerdo a la cantidad de datos que tenga
for hcanal := 0 to NcanalesCEGH-1 do
    begin

    //  5.2.1  Escribo el encabezado
    NombreHoja := Nombres[hcanal] + '_ensemble';
    ods.agregoHoja(NombreHoja);
    ods.go(NombreHoja,1,1);
    ods.write('Paso');
    for jmiembro := 1 to DatosEnsamble[hcanal].nc do
        ods.write('member_'+IntToStr(jmiembro));

    //  5.2.2  Escribir el paso de tiempo
    for kpaso := 1 to DatosEnsamble[hcanal].nf do
        begin
        ods.go(NombreHoja,1+kpaso,1);
        ods.Write( dtIni + (kpaso-1)*paso   );
        end;

    //  5.2.3  Escribir las cronicas del ensemble
    ods.Write_MatR(DatosEnsamble[hcanal],NombreHoja,2,2);

    end;


{ GFlieller@2019-05-23
// 5.3  Agrego cada hoja de acuerdo a la cantidad de datos que tenga
for hcanal := 0 to NcanalesCEGH-1 do
    begin

    //  5.3.1  Escribo el encabezado
    NombreHoja := Nombres[hcanal] + '_CEGH';
    ods.agregoHoja(NombreHoja);
    ods.go(NombreHoja,1,2);
    for jmiembro := 1 to DatosSimCEGH[hcanal].nc do
        ods.write('cron_'+IntToStr(jmiembro));

    //  5.3.2  Escribir el paso de tiempo
    for kpaso := 1 to DatosSimCEGH[hcanal].nf do
        begin
        ods.go(NombreHoja,1+kpaso,1);
        ods.Write( DateTimeToIsoStr(dtIni + (kpaso-1)*dtPaso  ) );
        end;

    //  5.3.3  Escribir las cronicas del ensemble
    ods.Write_MatR(DatosSimCEGH[hcanal],NombreHoja,2,2);

    end;

}

/////////////////////////////////////////////////////////////////

// 5.4  Fijo la cantidad de canales de cada corte de probabilidad
setlength(PE5_CEGH,length(DatosSimCEGH));
setlength(PE10_CEGH,length(DatosSimCEGH));
setlength(PE20_CEGH,length(DatosSimCEGH));
setlength(PE50_CEGH,length(DatosSimCEGH));
setlength(PE80_CEGH,length(DatosSimCEGH));
setlength(PE90_CEGH,length(DatosSimCEGH));
setlength(PE95_CEGH,length(DatosSimCEGH));
setlength(SeriesPE_CEGH,length(DatosSimCEGH));

for kSerie := 0 to (length(DatosSimCEGH)-1) do
  begin
  // 5.5  Ordeno las series para sacar cada corte de probabilidad
  DatosSimCEGH[kSerie].Sort();
  // 5.6  Cargo cada corte de porbabilidad a la serie correspondiente
  PE5_CEGH[kSerie]:= DatosSimCEGH[kSerie].Columna(round(0.05*DatosSimCEGH[kSerie].nc));
  PE10_CEGH[kSerie]:=DatosSimCEGH[kSerie].Columna(round(0.10*DatosSimCEGH[kSerie].nc));
  PE20_CEGH[kSerie]:=DatosSimCEGH[kSerie].Columna(round(0.20*DatosSimCEGH[kSerie].nc));
  PE50_CEGH[kSerie]:=DatosSimCEGH[kSerie].Columna(round(0.50*DatosSimCEGH[kSerie].nc));
  PE80_CEGH[kSerie]:=DatosSimCEGH[kSerie].Columna(round(0.80*DatosSimCEGH[kSerie].nc));
  PE90_CEGH[kSerie]:=DatosSimCEGH[kSerie].Columna(round(0.90*DatosSimCEGH[kSerie].nc));
  PE95_CEGH[kSerie]:=DatosSimCEGH[kSerie].Columna(round(0.95*DatosSimCEGH[kSerie].nc));
  // 5.7  Creo matrices con todos los cortes de probabilidad para cada canal
  SeriesPE_CEGH[kSerie]:=TMatR.Create_Init(336,7);
  SeriesPE_CEGH[kSerie].pon_columna(1,1,PE5_CEGH[kSerie].toTDAOfNReal());
  SeriesPE_CEGH[kSerie].pon_columna(1,2,PE10_CEGH[kSerie].toTDAOfNReal());
  SeriesPE_CEGH[kSerie].pon_columna(1,3,PE20_CEGH[kSerie].toTDAOfNReal());
  SeriesPE_CEGH[kSerie].pon_columna(1,4,PE50_CEGH[kSerie].toTDAOfNReal());
  SeriesPE_CEGH[kSerie].pon_columna(1,5,PE80_CEGH[kSerie].toTDAOfNReal());
  SeriesPE_CEGH[kSerie].pon_columna(1,6,PE90_CEGH[kSerie].toTDAOfNReal());
  SeriesPE_CEGH[kSerie].pon_columna(1,7,PE95_CEGH[kSerie].toTDAOfNReal());
  end;


//  5.8  Agrego cada hoja de acuerdo a la cantidad de datos que tenga
for hcanal := 0 to NcanalesCEGH-1 do
  begin

  //  5.8.1  Agrego la hoja
  NombreHoja := Nombres[hcanal] + '_PE_CEGH';
  ods.agregoHoja( NombreHoja );

  //  5.8.2  Escribo el encabezado
  ods.go(NombreHoja,1,1);
  ods.write('Paso');
  ods.write('PE5');
  ods.write('PE10');
  ods.write('PE20');
  ods.write('PE50');
  ods.write('PE80');
  ods.write('PE90');
  ods.write('PE95');

  //  5.8.3  Escribo los pasos de tiempo
  for kpaso := 1 to DatosEnsamble[hcanal].nf do
      begin
      ods.go(NombreHoja,1+kpaso,1);
      ods.Write( dtIni + (kpaso-1)*paso   );
      end;

  //  5.8.4  Escribo los cortes de probabilidad
  ods.Write_MatR(SeriesPE_CEGH[hcanal],NombreHoja,2,2);

  end;


// 5.9  Fijo la cantidad de canales de cada corte de probabilidad
setlength(PE5_IMFIA,length(DatosEnsamble));
setlength(PE10_IMFIA,length(DatosEnsamble));
setlength(PE20_IMFIA,length(DatosEnsamble));
setlength(PE50_IMFIA,length(DatosEnsamble));
setlength(PE80_IMFIA,length(DatosEnsamble));
setlength(PE90_IMFIA,length(DatosEnsamble));
setlength(PE95_IMFIA,length(DatosEnsamble));
setlength(SeriesPE_IMFIA,length(DatosEnsamble));

for kSerie := 0 to (length(DatosEnsamble)-1) do
  begin
  // 5.9  Ordeno las series para sacar cada corte de probabilidad
  DatosEnsamble[kSerie].Sort();
  // 5.10  Cargo cada corte de porbabilidad a la serie correspondiente
  PE5_IMFIA[kSerie] := DatosEnsamble[kSerie].Columna(round(0.05*DatosEnsamble[kSerie].nc));
  PE10_IMFIA[kSerie]:=DatosEnsamble[kSerie].Columna(round(0.10*DatosEnsamble[kSerie].nc));
  PE20_IMFIA[kSerie]:=DatosEnsamble[kSerie].Columna(round(0.20*DatosEnsamble[kSerie].nc));
  PE50_IMFIA[kSerie]:=DatosEnsamble[kSerie].Columna(round(0.50*DatosEnsamble[kSerie].nc));
  PE80_IMFIA[kSerie]:=DatosEnsamble[kSerie].Columna(round(0.80*DatosEnsamble[kSerie].nc));
  PE90_IMFIA[kSerie]:=DatosEnsamble[kSerie].Columna(round(0.90*DatosEnsamble[kSerie].nc));
  PE95_IMFIA[kSerie]:=DatosEnsamble[kSerie].Columna(round(0.95*DatosEnsamble[kSerie].nc));
  // 5.11  Creo matrices con todos los cortes de probabilidad para cada canal
  SeriesPE_IMFIA[kSerie]:=TMatR.Create_Init(PE50_IMFIA[kSerie].n,7);
  SeriesPE_IMFIA[kSerie].pon_columna(1,1,PE5_IMFIA[kSerie].toTDAOfNReal());
  SeriesPE_IMFIA[kSerie].pon_columna(1,2,PE10_IMFIA[kSerie].toTDAOfNReal());
  SeriesPE_IMFIA[kSerie].pon_columna(1,3,PE20_IMFIA[kSerie].toTDAOfNReal());
  SeriesPE_IMFIA[kSerie].pon_columna(1,4,PE50_IMFIA[kSerie].toTDAOfNReal());
  SeriesPE_IMFIA[kSerie].pon_columna(1,5,PE80_IMFIA[kSerie].toTDAOfNReal());
  SeriesPE_IMFIA[kSerie].pon_columna(1,6,PE90_IMFIA[kSerie].toTDAOfNReal());
  SeriesPE_IMFIA[kSerie].pon_columna(1,7,PE95_IMFIA[kSerie].toTDAOfNReal());
  end;


// 5.11  Creo las hojas de los cortes de probabilidad
for hcanal := 0 to NcanalesCEGH-1 do
  begin

  //  5.11.1  Agrego la hoja
  NombreHoja := Nombres[hcanal] + '_PE_IMFIA';
  ods.agregoHoja( NombreHoja );

  //  5.11.2  Escribo el encabezado
  ods.go(NombreHoja,1,1);
  ods.write('Paso');
  ods.write('PE5');
  ods.write('PE10');
  ods.write('PE20');
  ods.write('PE50');
  ods.write('PE80');
  ods.write('PE90');
  ods.write('PE95');

  //  5.11.3  Escribo los pasos de tiempo
  for kpaso := 1 to DatosEnsamble[hcanal].nf do
      begin
      ods.go(NombreHoja,1+kpaso,1);
      ods.Write( dtIni + (kpaso-1)*paso   );
      end;

  //  5.11.4  Escribo los cortes de probabilidad
  ods.Write_MatR(SeriesPE_IMFIA[hcanal],NombreHoja,2,2);

  end;


///////////////////////////////////////////

// 5.12  Guardo el ods
ods.Guardar;
// 5.13  Libero la memoria usada
ods.free;
for kSerie := 0 to (length(DatosSimCEGH)-1) do
begin
  PE5_CEGH[kserie].free;
  PE10_CEGH[kserie].free;
  PE20_CEGH[kserie].free;
  PE50_CEGH[kserie].free;
  PE80_CEGH[kserie].free;
  PE90_CEGH[kserie].free;
  PE95_CEGH[kserie].free;
  SeriesPE_CEGH[kSerie].free;
end;
for kSerie := 0 to (length(DatosEnsamble)-1) do
begin
  PE5_IMFIA[kserie].free;
  PE10_IMFIA[kserie].free;
  PE20_IMFIA[kserie].free;
  PE50_IMFIA[kserie].free;
  PE80_IMFIA[kserie].free;
  PE90_IMFIA[kserie].free;
  PE95_IMFIA[kserie].free;
  SeriesPE_IMFIA[kSerie].free;
end;

SetLength(PE5_CEGH ,0);
SetLength(PE10_CEGH,0);
SetLength(PE20_CEGH,0);
SetLength(PE50_CEGH,0);
SetLength(PE80_CEGH,0);
SetLength(PE90_CEGH,0);
SetLength(PE95_CEGH,0);
SetLength(SeriesPE_CEGH,0);
SetLength(PE5_IMFIA,0);
SetLength(PE10_IMFIA,0);
SetLength(PE20_IMFIA,0);
SetLength(PE50_IMFIA,0);
SetLength(PE80_IMFIA,0);
SetLength(PE90_IMFIA,0);
SetLength(PE95_IMFIA,0);
SetLength(SeriesPE_IMFIA,0);


end;

procedure cargarEnsamblesLocales(var SeriesDeEnsambles: array of TMatR;
  archidir: string; Nseries: integer);
  // Carga los ramilletes de pronosticos del IMIFIA a partir de archivo .xlt
  // guardados localmente. Se usa para debuguear el programa.
begin
  SeriesDeEnsambles[0] := TMatR.CreateFromArchiXLT(archidir+'BonseriesIMFIA.xlt');
  if NSeries>1 then
     SeriesDeEnsambles[1] := TMatR.CreateFromArchiXLT(archidir+'PalseriesIMFIA.xlt');
  if NSeries>2 then
     SeriesDeEnsambles[2] := TMatR.CreateFromArchiXLT(archidir+'SGseriesIMFIA.xlt');

end;


{$IFDEF PRONO_HIDRO_IMFIA}
procedure cargarEnsamblesDesdeEstado(var SeriesDeEnsambles: array of TMatR;
  estado: TEstado_SIN_UY; Nseries: integer);
  // Carga los ramilletes de pronosticos del IMFIA a partir del estado. En el
  // estado debe estar precargados los ramilletes a traves de consulta sql y
  // corrida del modelo del IMFIA.
begin
  SeriesDeEnsambles[0] := TMatR.Create_Clone(estado.QBonIMFIAmatHora);
  if NSeries>1 then
     SeriesDeEnsambles[1] := TMatR.Create_Clone(estado.QPalIMFIAmatHora);
  if NSeries>2 then
     SeriesDeEnsambles[2] := TMatR.Create_Clone(estado.QSGOficialMatHr);
end;
{$ENDIF}

procedure graficarResultados;
var
parametros : array of string;
begin
      // correr el script de octave
      {$IFDEF WINDOWS}
        RunChildAndWait('octave-cli', ['graficar_resultados.m']);
      {$ELSE}
        setlength(parametros, 2);
        parametros[0] := '--no-window-system';
        parametros[1] := '--no-gui';
        parametros[1] := '/home'+ DirectorySeparator +
                      GetEnvironmentVariable('USER') + DirectorySeparator +
                      'prono_hidro_imfia/graficar_resultados.m';

        RunChildAndWait('/usr/bin/octave', parametros);
      {$ENDIF}
end;


{$IFDEF PRONO_HIDRO_IMFIA}
function getEnsemblesAportesVates(dirSalida: string; dbpqPronos: TDBPQCon;
  estado: TEstado_SIN_UY): TEnsembles;
{
Este procedimiento es especifico para obtenr los ensembles de aportes completos
para ser usados en la sala Vates.
Al final se tendra un archivo .xlt, que indica en el encabezado el nombre de la
fuente, el numero de canales, el numero de miembros que tiene el ensemble,
el numero de paso de tiempo, el DT de la primera muestra y el paso de tiempo.
Cada columna correspondera a un miembro y cada fila a un paso de tiempo. El
archivo salida se guarda en el directorio dirSalida.
Este archivo puede ser leido usando la funcion TEnsembles.CreateFromArchiXLT. El
robot de Vates debería correr la funcion getEnsemblesAportesVates, la funcion
TEnsembles.CreateFromArchiXLT y luego la funcion calibrarCEGHconEnsemble.

En el caso del archivo de salida de este procedimiento:
Nombre: Aportes
Ncanales: 6
Nmiembros: 21 (cntl + 20 miembros perturbados)
Npasos: 14*24
DtPrimeraMuestra: ()
Paso: 1/24

Bonete

Palmar

SG CTM

in34

cmo1

cmo3
}

var
//  dtPronoIni : TDateTime;
  diaInicioProno : integer;
  ensembles:TEnsembles;
  NombreFuente : string;

begin

  NombreFuente := 'Aportes';

  //  Crear Estado

//  dtPronoIni := estado.opciones.dtProno;
  diaInicioProno := trunc(Now);
//  OJO ACA ESTOY HACIENDO PRUEBAS GFLIELLER@2019-09-17
//  diaInicioProno := 43714;

  writeln('Obteniendo ensambles de pronosticos de aportes.');

  //  1-  Bajar el ensemble de Salto Grande CTM al estado
  estado.getEnsembleAportesSaltoOficial;

  estado.getiN34(dbcon_rosx);

  estado.getCMO(dbpqPronos);

  //  2-  Bajar Pronosticos, Calcular arpotes y cargarlos al estado
  estado.getPronoHidroIMFIA(dbpqPronos,diaInicioProno, true);   //

//  estado.getPronoBaygorriaIMFIA;

  //  3-  Escribir los ensembles de aportes en formato xlt
  //  3.1- Cargar los ensambles desde el estado. Los miembrso son Nmiembros + el miembro de control.
  ensembles:=TEnsembles.Create( NombreFuente,  NCANALES_,NMIEMBROS+1,
                                (NDIASPRONOSHIDRO-1)*24, diaInicioProno,
                                DUR_PASO_);
  ensembles.CargarDesdeEstado(estado);

  ensembles.WriteToArchiXLT(dirSalida);

  Result := ensembles;

end;
{$ENDIF}


{ TProblemaResolverFactores }

constructor TProblemaResolverFactores.Create(aM: TMatR);
var
  hcanal: Integer;
begin
// 3.4.1  Creo el problema como hijo de la clase TProblema_m01
inherited Create_init(1, aM.nc + 1, 0, nil, nil);

// 3.4.2  Cargo la matriz M y Ncanales para tenerlos disponibles
M := aM;
Ncanales := M.nc;

// 3.4.3  Creo la función objetivo en base a la matriz M
f := T_fobjetivoCEGH.Create(M);

// 3.4.4  Cargo el valor inicial de B
Betas:= TVectR.Create_InitVal(Ncanales,BETA_INI);  // Lambda base 0.01

// 3.4.5  Fijo las restricciones de caja
for hcanal := 1 to Ncanales do
    begin
    cota_inf_set(hcanal, BETA_MIN);   //  Lambda Maximo 100
    cota_sup_set(hcanal, BETA_MAX);  //  Lambda Minimo 1e-40
    end;

end;

procedure TProblemaResolverFactores.free;
begin
//  Betas.free;
  inherited free;
end;

function TProblemaResolverFactores.resolverFactores(var factores:TMatR): boolean;
//  Loop de resolucion de factores por paso
var
X : TVectR;
l : TVectR;
Convergio: boolean;
hcanal: integer;
cont_Iters: longint;
fCosto, lcosto, SumaFila: NReal;
begin

// 3.5.1
// Creo un vector de restricciones lambda (l) vacio
l := TVectR.Create_init(1);
l.pon_e(1, 0);

// 3.5.2
// Creo el vector X de resolucion del problema. Lo inicio en el valor de Beta.
// En el primer paso tiene un valor precargado, para el resto utiliza como punto
// de partida el resultado en el paso anterior
X := TVectR.Create_Init(Ncanales);
for hcanal:= 1 to Ncanales do
begin
{
     if M.Fila(hcanal).normSumAbs < 1e-20 then
        if M.Fila(hcanal).normSumAbs = 0 then
           X.pon_e(hcanal, 1e+60)
        else
        X.pon_e(hcanal, 1e+9)
     else
     X.pon_e(hcanal,Betas.e(hcanal));
}
     SumaFila := M.Fila(hcanal).normSumAbs;
     if SumaFila < SUMA_0 then
        X.pon_e(hcanal, BETA_MAX)
     else
     X.pon_e( hcanal,sqrt( 1/SumaFila ) )  ;

end;

cont_Iters:=0;
// 3.5.3
{GFlieller@2019-06-14 20:29 Cambio a MinInBox: Mejora considerablemente}
// Como no tiene restricciones es mejor utilizar el MinInBox
MinInBox_CPX(X, 1, 1e-60, 500000, cont_Iters, lcosto, fcosto, convergio, false);

// 3.5.4
// Cargo en Betas el resultado de la optimizacion
Betas := TVectR.Create_Clone(X);

// 3.5.5
// Cargo los factores: 1/Beta
for hcanal:=1 to Ncanales do
begin
     if Betas.e(hcanal) > BETA_0 then
        factores.pon_e(hcanal,1,0)
     else
     factores.pon_e(hcanal,1,1/Betas.e(hcanal));
end;

// 3.5.6
// Libero la memoria usada
l.free;
X.free;
result:= convergio;

end;


{ T_fobjetivoCEGH }

constructor T_fobjetivoCEGH.Create(aM: TMatR);
// 3.4.3
begin
// 3.4.3.0 Creo el T_fobjetivoCEGH en base al padre
     inherited Create;

// 3.4.3.1  Cargo la matriz del problema
self.M := aM;

self.Ncanales := M.nc;

// 3.4.3.2  Construyo FX que es el vector FX = M * X -"1/X", que buscamos minimizar.
self.FX  := TVectR.Create_Init(Ncanales);


end;

function T_fobjetivoCEGH.f(const X: TVectR): NReal;
// Calculo del error
var
X_1: TVectR;
hcanal: Integer;
error : NReal;
begin

// Esto se usa adentro del MaxInBox que es el 3.5.3
// 3.5.3.1  Creo X_1 que es lo que tambien se llamo "1/X"
X_1 := TVectR.Create_Init(Ncanales);
for hcanal:=1 to Ncanales do
    begin
    X_1.pon_e(hcanal,-1/X.e(hcanal));
    end;

// 3.5.3.2  Calculo M * X - "1/X" que es la funcion de error
FX.Mult(M,X);                        // FX = M * X
FX.Suma(FX,X_1);                     // FX = M * X - "1/X"

// 3.5.3.3  Se da como resultado la norma euclidea por un factor amplificador
//          para poder manejar el error evitando errores numericos.
error:= FX.normEuclid;         // error = 1e+08*RAIZ(FX * FX^T)

X_1.free;
result := error;

end;

destructor T_fobjetivoCEGH.destroy;
begin
  FX.free;
end;

procedure T_fobjetivoCEGH.acum_g(var grad: TVectR; const X: TVectR);
// Calculo del gradiente del error
var
  kderivada, hcanal: Integer;
  derivada: NReal;
begin

// 3.5.3.4
// error = 1e+08*RAIZ(FX * FX^T) = 1e+08*RAIZ( suma(f_j(x)^2) )
// derror/dx_i = ( 1e+08/(2*RAIZ(error)) ) * suma( 2*f_j(x)*df_j/dx_i))
// derror/dx_i = ( 1e+08/(2*RAIZ(error)) ) * suma( 2*f_j(x)*(
for kderivada:=1 to Ncanales do
    begin
    derivada := 0;
    if FX.normEuclid <> 0 then
    begin
        for hcanal:=1 to Ncanales do
             derivada := derivada + 2*FX.e(hcanal)*M.e(hcanal,kderivada);

        derivada := derivada + 2*FX.e(kderivada)*(1/power(X.e(kderivada),2));

        derivada := (derivada/(2*FX.normEuclid));
    end;
    grad.acum_e(kderivada,derivada);

    end;

end;


{ TcalibrarConoEnsamble }

constructor TcalibrarConoEnsamble.Create(acegh: TFuenteSintetizadorCEGH;
            aEnsambles: array of TMatR; aFecha :TFecha);

var
  hcanal: Integer;
begin
inherited Create;


self.fechaIni:= aFecha;

self.fcegh := acegh;

cargarDimensiones(aEnsambles);

SetLength(Ensambles,Ncanales);
for hcanal:=0 to Ncanales-1 do
    Ensambles[hcanal] := TMatR.Create_Clone(aEnsambles[hcanal]);

// 2.1 C es la misma matriz que se menciona en la documentación
cargarFiltros;

// 2.2 Pasa ensamble a Gaussiano
cargarEnsamblesGausiano(aEnsambles);

publicarEnsemble;

// 2.3 Calcula y carga las guias
calcularGuias;

// 2.4 Calcula y carga los desvios respecto a la guia
calcularDesvios;

// 2.5 Calculo la matriz de covarianzas
// Para cada miembro del ensamble y para cada paso se calcula des_desT := des * des^T. desT = des^T
calcularCovar;

// 2.6 Calcula los sesgos
calcularSesgos;

end;

function TcalibrarConoEnsamble.calibrarConoEnsamble(): boolean;

var
Z0, Z, LambdaAux, LambdaResolvedor: TMatR;
kpaso, k, elemento : integer;
convergio, matrizRedNula: boolean;
Problema: TProblemaResolverFactores;
idx : TDAofNInt;

begin
// 3.1 Inicializo
Setlength(betas, Npasos);
setlength(Lambdas,Npasos);
setlength(M,Npasos);
Z0 := TMatR.Create_Init(Ncanales,Ncanales);
for kPaso := 0 to Npasos-2 do
  begin
  // 3.2 Calculo la matriz de ajuste del paso anterior Z
  Z := calcularZ(Z0,kPaso);
  // 3.3 Calculo M en base a la covarianza y Z
  armarProblema(Z, kPaso);

  //  Si los desvios fueran todos 0 entonces no se puede resolver por el metodo numerico
  // 3.4 Creo el problema calculo de Beta en base a la matriz M
  LambdaAux := TMatR.Create_InitVal(Ncanales,1,1e-20);
  {$IFDEF REDUCIR_SISTEMA}
  // 3.4.1 Si arma el problema
  // 3.4.2 Reduce la el problama. Filas de 0, tribialidades o incompatibilidades.

  if kPaso = 0 then
    convergio := true
  else
    begin

    setlength(idx,Ncanales);
    matrizRedNula:= false;
    Mreducida := ReducirProblema(M[kPaso], idx, LambdaAux, matrizRedNula);



    if matrizRedNula = false then
      begin
      Problema := TProblemaResolverfactores.Create(MReducida);

      LambdaResolvedor := TMatR.Create_InitVal(MReducida.nf,1,1e-20);
      // 3.5 Resolver los valores de Lambda en base a resolver los Beta
      if kPaso = 0 then
        convergio := true
      else
        convergio := Problema.resolverFactores(LambdaResolvedor);

{      elemento := 0; }

      elemento := 1;
      for k:=1 to Ncanales do
        begin
        if (kInArray(idx,k)<>-1) then
            begin
            LambdaAux.pon_e( k, 1, LambdaResolvedor.e( elemento, 1) );
            elemento := elemento + 1;
            end;

{          if k= idx[elemento] then
            begin
            elemento := elemento+1;
            LambdaAux.pon_e(k,1, LambdaResolvedor.e(j,1));
            end;  }
        end;

      Problema.free;
      end;
  end;
  {$ELSE}
  Problema := TProblemaResolverfactores.Create(M[kpaso]);
  // 3.5 Resolver los valores de Lambda en base a resolver los Beta
  convergio := Problema.resolverFactores(LambdaAux);
  Problema.free;
  {$ENDIF}
  writeln( FloatToStr( 100*kPaso/Npasos ) + ' %' );
  Lambdas[kpaso] := TMatR.Create_Clone(LambdaAux);
  LambdaAux.free;

  // 3.6 Preservo Z para el paso siguiente como informacion del paso previo
  Z0 := Z;
  end;
Z0.free;
// 3.7. Cargo los factores lambda a un ods
ImprimirLambdas;

// 3.8 Convierte los resultados al formato que usa el CEGH
ArmarResultados();

result := convergio;

end;


procedure TcalibrarConoEnsamble.cargarDimensiones(aEnsambles: array of TMatR);
begin
Ncanales := length(aEnsambles);
Nmiembros := aEnsambles[0].nc;
Npasos := aEnsambles[0].nf;
pasosTotales := round( (fcegh.globs.fechaFinOpt.dt - fcegh.globs.fechaIniOpt.dt)
                                                / fcegh.globs.dt_DelPaso);
end;

procedure TcalibrarConoEnsamble.cargarEnsamblesGausiano(aEnsambles: array of TMatR);
var
hcanal,jmiembro,kpaso , pasoInicial: integer;

begin

// 2.2.1 Paso cada ensamble a mundo gaussiano, en forma de array de matrices.

setlength(EnsamblesGauss, length(aEnsambles));
pasoInicial := fechaIni.horasDesdeElInicioDelAnio;
writeln( 'Paso inicial: '+ IntToSTr(pasoInicial) + ' ('+
         DateToStr(fechaIni.dt)+').');

for hcanal := 1 to Ncanales do
    begin
        EnsamblesGauss[hcanal-1] := TransformarEnsamble(fcegh.datosModelo_Sim,
                                    aEnsambles, hcanal);
    end;

// 2.2.2 Convierto los arrays de matrices en arrays de array de matrices columna
//       para poder realizar multiplicaciones mas facilmente
// Tengo m (M = Ncanales) matrices de r*p (r=Nmiembros,p=Npasos). Los vectores
// X_k deben ser de m*1 y ser p vectores, 1 por paso, para los r miembros.
// X[miembro][paso] = (X_jkc1, X_jkc2,..., X_jkcm)
setlength(Xg, Nmiembros);
for jmiembro := 1 to Nmiembros do
begin
    setlength(Xg[jmiembro-1],Npasos);
    for kpaso := 1 to Npasos do
    begin
        Xg[jmiembro-1][kpaso-1] := TMatR.Create_InitVal(Ncanales,1,0);
        for hcanal := 1 to Ncanales do
        begin
            Xg[jmiembro-1][kpaso-1].pon_e( hcanal,1, EnsamblesGauss[hcanal-1].e(kpaso,jmiembro) );
        end;
    end;
end;

end;

procedure TcalibrarConoEnsamble.calcularGuias;

var
guias: array of TDAofNReal;
guiasgaussianas : array of TDAofNReal;
hcanal, kpaso , kPasoHorarioDelAnio, horasIniCono: integer;

begin

// 2.3.1
// Creo el pronostico en base a los datos del modelo de simulacion del CEGH
aPronosticos:= Tpronosticos.Create_Default(fcegh.datosModelo_Sim);
setlength(aPronostico,Ncanales);
setlength(guias,Ncanales);
setlength(guiasGaussianas,Ncanales);

// 2.3.2 creo y cargo las guias gaussianas y reales y las cargo al pronostico
for hcanal:=0 to (Ncanales-1) do
begin

  setlength(guiasGaussianas[hcanal],pasosTotales);
  setlength(guias[hcanal],pasosTotales);
  // Armo la guia en el mundo gaussiano y paso su transformada al mundo real
  horasIniCono := fechaIni.horasDesdeElInicioDelAnio;
  for kpaso:=0 to Npasos-1 do
  begin
    kPasoHorarioDelAnio := horasIniCono + kpaso;
    kPasoHorarioDelAnio := kPasoCiclico(kPasoHorarioDelAnio,24*364,0);
    guiasGaussianas[hcanal][kpaso] := EnsamblesGauss[hcanal].Fila(kpaso+1).promedio;
    guias[hcanal][kpaso] := fcegh.datosModelo_Sim.gTox(guiasGaussianas[hcanal][kpaso],
                              hcanal+1,kPasoHorarioDelAnio);
  end;

end;

// 2.3.3 Calculo Xesp = <X_k>, p vectores de m*1 entradas, promedios por canal de miembros
setlength(Xesp,pasosTotales);
for kpaso := 0 to Npasos-1 do
begin
    Xesp[kpaso] := TMatR.Create_Init(Ncanales,1);
    for hcanal := 0 to Ncanales-1 do
        Xesp[kpaso].pon_e(hcanal+1,1,guiasgaussianas[hcanal][kpaso]);
end;

// 2.3.4 Completo el resto de los pasos para que decaiga con A
for kpaso := Npasos to pasosTotales-1 do
begin
    Xesp[kpaso] := TMatR.Create_Clone(Xesp[kpaso - 1]);
    Xesp[kpaso].Mult(A,Xesp[kpaso-1]);
end;

// 2.3.5 Completo las guias
for hcanal:=0 to (Ncanales-1) do
begin
  for kpaso:= Npasos to pasosTotales - 1 do
  begin
      kPasoHorarioDelAnio := horasIniCono + kpaso;
      kPasoHorarioDelAnio := kPasoCiclico(kPasoHorarioDelAnio,24*364,0);
      guiasGaussianas[hcanal][kpaso] := Xesp[kpaso].e(hcanal+1,1);
      guias[hcanal][kpaso] := fcegh.datosModelo_Sim.gTox(guiasGaussianas[hcanal][kpaso],
                                hcanal+1,kPasoHorarioDelAnio);
  end;

  // 2.3.6 cargo las guias
  // Cargo la guia en el aPronosticos de la clase TcalibrarConoEnsamble
  aPronostico[hcanal] := TPronostico.Create_Default(fcegh.datosModelo_Sim);
  aPronostico[hcanal].guia_pe:=0.5;
  aPronostico[hcanal].guia_eg:=guiasgaussianas[hcanal];
  aPronostico[hcanal].guia:=guias[hcanal];

end;

end;


procedure TcalibrarConoEnsamble.calcularSesgos;
// Calculo los sesgos
var
kpaso: integer;
Saux: TMatR;

begin
// 2.6
setlength(S,pasosTotales-1);
for kpaso := 1 to pasosTotales-1 do
  begin
    Saux:= TMatR.Create_Clone(Xesp[kPaso-1]);
    Saux.Mult(A,Saux);                                          // A*X0
    S[kpaso-1]:= TMatR.Create_Combinar(1,Xesp[kpaso],-1,Saux);  // S = X1 - A*X0
    Saux.free;
  end;

end;


procedure TcalibrarConoEnsamble.calcularCovar;
// 2.5.
// Calculo las matrices covarianzas
// p matrices de covarianza de dimesion m x m, promedio de las matrices de des * des^T
// En la documentación se le llamo Xi
var
jmiembro,kpaso : integer;
des_desT, desT: array of array of TMatR;

begin
// 2.5.1 Calculo d*d^T para todos los pasos
setlength(desT, Nmiembros);
setlength(des_desT, Nmiembros);
for jmiembro := 0 to (Nmiembros-1) do
begin
  setlength(desT[jmiembro],Npasos);
  setlength(des_desT[jmiembro],Npasos);
  for kpaso := 0 to (Npasos-1) do
  begin
    desT[jmiembro][kpaso] := TMatR.Create_Clone(des[jmiembro][kpaso]);
    desT[jmiembro][kpaso].Transponer;
    des_desT[jmiembro][kpaso] := TMatR.Create_InitVal(Ncanales,Ncanales,0);
    des_desT[jmiembro][kpaso].Mult(des[jmiembro][kpaso],desT[jmiembro][kpaso]);
  end;
end;

// 2.5.2 Calculo las matrices de covarianza para todos los pasos
setlength(covar, Npasos);
for kpaso := 0 to (Npasos-1) do
begin
  covar[kpaso] := TMatR.Create_InitVal(Ncanales,Ncanales,0);
  for jmiembro := 0 to Nmiembros-1 do
    covar[kpaso].Suma(covar[kpaso],des_desT[jmiembro][kpaso]);

  covar[kpaso].PorReal(1/(Nmiembros-1));

end;  // Debug: Esta igual que en el metodo B -- GF@2019-02-12 11:58


// 2.5.3 Libero la memoria que utilice
for jmiembro:=0 to (Nmiembros-1) do
begin
    for kpaso:=0 to (Npasos-1) do
    begin
        desT[jmiembro][kpaso].Free;
        des_desT[jmiembro][kpaso].Free;
    end;
    SetLength(desT[jmiembro],0);
    SetLength(des_desT[jmiembro],0);
end;
SetLength(desT,0);
SetLength(des_desT,0);

end;


procedure TcalibrarConoEnsamble.calcularDesvios;
//  Calculo de los desvios por paso de tiempo
var
X_normal_0 : array of array of TMatR;
X_normal_0_promed : array of TMatR;
p : NReal;
jmiembro, kpaso, hcanal: Integer;


// 2.4.1 Normalizo los valores en el mundo gaussiano con respecto al promedio
begin
setlength(X_normal_0, Ncanales);
for hcanal:=0 to Ncanales-1 do
begin
SetLength(X_normal_0[hcanal],Npasos);
    for kpaso:=0 to Npasos-1 do
    begin
        X_normal_0[hcanal][kpaso] := TMatR.Create_Init(Nmiembros,1);
        for jmiembro := 0 to Nmiembros-1 do
          X_normal_0[hcanal][kpaso].pon_e( jmiembro+1,1,Xg[jmiembro][kpaso].e(
                                           hcanal+1,1)-Xesp[kpaso].e(hcanal+1,1) );
        X_normal_0[hcanal][kpaso].Sort();
    end;
end;    // Esto esta igual que en el metodo b. (X-guia) GF@2019-02-11-11:46


// 2.4.2 Calculo el promedio por paso de la X normalizada
SetLength(X_normal_0_promed,Npasos);
for kpaso := 0 to Npasos-1 do
begin
X_normal_0_promed[kpaso] := TMatR.Create_Init(Ncanales,1);
    for hcanal := 0 to Ncanales-1 do
      begin
      p:=0;
           for jmiembro:=0 to Nmiembros-1 do
               p:= p + X_normal_0[hcanal][kpaso].e(jmiembro+1,1);
      p:=p/Nmiembros;
      X_normal_0_promed[kpaso].pon_e(hcanal+1,1, p )  ;
      end;
end;

// 2.4.3 Cargo las desviaciones
SetLength(des,Nmiembros);
for jmiembro := 0 to Nmiembros-1 do
  begin
    SetLength(des[jmiembro],Npasos);
    For kpaso := 0 to Npasos-1 do
    begin
      des[jmiembro][kpaso]:= TMatR.Create_Init(Ncanales,1);
      for hcanal:=1 to Ncanales do
        des[jmiembro][kpaso].pon_e( hcanal,1,
                                    X_normal_0[hcanal-1][kpaso].e(jmiembro+1,1)
                                    - X_normal_0_promed[kpaso].e(hcanal,1) );
    end;
  end;

// 2.4.4  Libero memoria
for kpaso := 0 to Npasos-1 do
begin
   X_normal_0_promed[kpaso].free;
end;
SetLength(X_normal_0_promed,0);

for hcanal:=0 to Ncanales-1 do
begin
    for kpaso:=0 to Npasos-1 do
    begin
         X_normal_0[hcanal][kpaso].Free;
    end;
    SetLength(X_normal_0[hcanal],0);
end;
SetLength(X_normal_0,0);

end;



function TcalibrarConoEnsamble.calcularZ(Z0: TMatR; kpaso: integer) : TMatR;
var
Z, AT, BT, L, Zaux:TMatR;
hcanal: Integer;
begin
  if kpaso = 0 then
     // 3.2.1 El Z0 se calcula como Z0 = A * d * d^T * A^T
     begin
     Z := TMatR.Create_Clone(A);
     Z.Transponer;
     Z.Mult(covar[kpaso],Z);
     Z.Mult(A,Z);
     end
  else
     begin
     // 3.2.2 Z = A *( Z0 + L * B * B^T * L^T )* A^T
     // 3.2.2.1 Calculo las traspuestas
     AT := TMatR.Create_Clone(A);
     AT.Transponer;
     BT := TMatR.Create_Clone(B);
     BT.Transponer;

     // 3.2.2.2 Calculo la matriz L, la cual es diagonal. L^T = L
     L := TMatR.Create_InitVal(Ncanales,Ncanales,0);
     for hcanal := 1 to Ncanales do
         L.pon_e(hcanal,hcanal,Lambdas[kpaso-1].e(hcanal,1));

     // 3.2.2.3 Calculo Zaux = L * B * B^T * L^T
     Zaux:= TMatR.Create_Clone(Z0);
     Zaux.Mult(BT,L);
     Zaux.Mult(B,Zaux);
     Zaux.Mult(L,Zaux);

     // 3.2.2.4 Calculo Zaux = Z0 + L * B * B^T * L^T
     Zaux.Suma(Z0,Zaux);

     // 3.2.2.4 Calculo Z = A *( Z0 + L * B * B^T * L^T )*A^T
     Zaux.Mult(Zaux,AT);
     Z := TMatR.Create_InitVal(Ncanales,Ncanales,0);
     Z.Mult(A,Zaux);

     // 3.2.2.5 Libero la memoria que use
     AT.Free;
     BT.Free;
     Zaux.Free;
     L.free;
     end;

result := Z;

end;


procedure TcalibrarConoEnsamble.cargarFiltros;
// Cargo las matrices A y B. Calculo C = (B * B^T)^-1.
begin
  A := fcegh.datosModelo_Sim.A_cte;
  B := fcegh.datosModelo_Sim.B_cte;

  C := TMatR.create_clone(B);
  C.Transponer;
  C.Mult(B,C);
  C.Inv;

end;


procedure TcalibrarConoEnsamble.armarProblema(Z : TMatR; kpaso : integer);
// 3.3 Calculo la matriz M del problema

//var
//  Cinv, M_: TMatR;
begin
if kpaso=0 then
   // 3.3.1 M0 = d*d^T . C
   begin
   M[kpaso] := TMatR.Create_Clone(covar[kpaso]);
   M[kpaso].MaximizarElementoAElemento(0);
   M[kpaso].porMatrizElementoAElemento(C);

   end
else
   // 3.3.2 M = (d*d^T - Z). C
   begin
   M[kpaso] := TMatR.Create_Combinar(1,covar[kpaso],-1,Z);
   M[kpaso].MaximizarElementoAElemento(0);
   M[kpaso].porMatrizElementoAElemento(C);
   end;

end;


procedure TcalibrarConoEnsamble.ImprimirLambdas;
// Guarda en un archivo xlt los valores de los factores lambda

var
Lambda : TMatR;
hcanal, kpaso : integer;
fechaAhora: TFecha;
begin

// 3.7.1  Creo una matriz para conter en cada columna los lambda de cada canal
//       de cada paso
Lambda:=TMatR.Create_Init(Ncanales,Npasos);
for kpaso:=1 to Npasos-1 do
    for hcanal:=1 to Ncanales do
    begin
        Lambda.pon_e(hcanal,kpaso, Lambdas[kpaso-1].e(hcanal,1));
    end;

// 3.7.2  Guardo la matriz en xlt
fechaAhora:= TFecha.Create_Dt(Now);
Lambda.WriteArchiXLT('/home/' + GetEnvironmentVariable('USER') + '/prono_hidro_imfia/Resultados/Lambdas/Lambdas_'+fechaAhora.AsAAAAMMDDhhmmtr+'.xlt');

//  3.7.3  Libero la memoria usada
Lambda.free;
fechaAhora.free;
end;


function TcalibrarConoEnsamble.TransformarEnsamble(aModelo:TModeloCEGH ;aEnsambles: array of TMatR; hcanal: integer): TMatR;
var
jmiembro,kpaso,kPasoHorarioDelAnio, horasIniCono:integer;
MatAux:TMatR;
begin
MatAux:=TMatR.Create_Init(Npasos,Nmiembros);

// writeln('Paso inicial: '+ IntToSTr(fcegh.escenariosDePronosticos[0].fechaIniCono.horasDesdeElInicioDelAnio));
  for jmiembro := 1 to Nmiembros do
  begin
    for kpaso := 0 to Npasos-1 do
    begin
    horasIniCono := fechaIni.horasDesdeElInicioDelAnio;
    kPasoHorarioDelAnio := horasIniCono + kpaso;
    kPasoHorarioDelAnio := kPasoCiclico(kPasoHorarioDelAnio,24*364,0);
    MatAux.pon_e(kpaso+1,jmiembro, aModelo.xTog(
                                 aEnsambles[hcanal-1].e(kpaso+1,jmiembro),
                                 hcanal, kPasoHorarioDelAnio) );
    end;
  end;

result := TMatR.Create_Clone(MatAux);
MatAux.free;
end;


procedure TcalibrarConoEnsamble.ArmarResultados();
//  Pasa los resultados de la calibracion a
//  los formatos en que los utiliza el CEGH
var
hcanal, kpaso , pasosTotales, pasosCono: integer;
l: NReal;
// 3.8   Los sesgos, guias y factores tienen que ir cargados en Pronostico como
//       TDAofNReal. Aca nos ocupamos de cargar unicamente los valores de los
//       sesgos y factores, las guias se cargan aparte.
begin
// 3.8.1  Para cada canal se crea un TDAofNReal.
setlength(sesgos, Ncanales);
setlength(factores,Ncanales);

// 3.8.2  Pasos de optimizacion del CEGH. No importa si los del pronostico son
//        mas, no se van a usar, y si fueran menos los tengo que rellenar
pasosTotales := round( (fcegh.globs.fechaFinOpt.dt - fcegh.globs.fechaIniOpt.dt)
                                                / fcegh.globs.dt_DelPaso);

// por si los pasos del pronostico fueran mas que los del cegh
pasosCono := min(Npasos,pasosTotales);

for hcanal:=0 to (Ncanales-1) do
begin
// 3.8.2  Cada TDofNReal tiene como largo la cantidad de pasos de opt del CEGH
  setlength(sesgos[hcanal],pasosTotales);
  setlength(factores[hcanal],pasosTotales);
  //  3.8.2.1  Escribo los pasos que conozco
  for kpaso:=0 to pasosCono-2 do
  begin
//      sesgos[hcanal][kpaso] := S[kpaso].e(hcanal+1,1);
      factores[hcanal][kpaso] := Lambdas[kpaso].e(hcanal+1,1);
  end;
  //  3.8.2.2  Relleno los pasos que no conozco
  for kpaso:=pasosCono-1 to pasosTotales-2 do
  begin
     // sesgos[hcanal][kpaso] := 0;
     // factores[hcanal][kpaso] := Lambdas[Npasos-2].e(hcanal+1,1) +
     //                            ( 1 - Lambdas[Npasos-2].e(hcanal+1,1) )*
     //                            (( kPaso - Npasos )/( pasosTotales - Npasos ));
     l := Lambdas[Npasos-2].e(hcanal+1,1)*(kPaso - pasosCono)/(24*7);
     factores[hcanal][kpaso] := min(1, l );
  end;

  for kpaso:=0 to pasosTotales-2 do
  begin
      sesgos[hcanal][kpaso] := S[kpaso].e(hcanal+1,1);
  end;
end;

// 3.8.3  Cargo al Pronostico los parametros del cono. NPSA = -1 indica que
//        los parametros no se deben considerar. Lo cargo en NPSA porque los
//        otros parametros pueden usarse para otras cosas
//        //  GF@2019-07-12 Agrego flg_modo_cono. flg= 1 sustituye a NPSA= -1
for hcanal := 0 to Ncanales-1 do
begin
 aPronostico[hcanal].NPCC := pasosTotales;
 aPronostico[hcanal].NPSA := -1;
 aPronostico[hcanal].NPAC := pasosTotales;
 aPronostico[hcanal].flg_modo_cono := 1;
// aPronostico[hcanal].NPCC := 181;
// aPronostico[hcanal].NPLC := 0;
// aPronostico[hcanal].NPSA := 72;
// aPronostico[hcanal].NPAC := 200;
// aPronostico[hcanal].flg_modo_cono := 1;

// 3.8.4  Cargo los sesgos y factores a cada Pronostico
setlength(aPronostico[hcanal].factor, pasosTotales);
vclear(aPronostico[hcanal].factor);
aPronostico[hcanal].factor:= factores[hcanal];
vclear(aPronostico[hcanal].sesgo);
aPronostico[hcanal].sesgo := sesgos[hcanal];

// 3.8.5  Cargo cada Pronostico al Pronosticos
aPronosticos.setItem(hcanal,aPronostico[hcanal]);
end;
aPronosticos.Probabilidad:=1;
aPronosticos.fechaIniCono.dt := fcegh.escenariosDePronosticos[0].fechaIniCono.dt;
aPronosticos.prepararse(1);

// 3.8.6  Libero memoria innecesaria
for kpaso:=0 to Npasos-2 do
begin
     S[kpaso].free;
     Lambdas[kpaso].free;
end;
setlength(S,0);
setlength(Lambdas,0);

end;

procedure TcalibrarConoEnsamble.free;
var
k, j:integer;
begin
  if A<>nil then
     A.free;
  if B<>nil then
     B.free;
  if C<>nil then
     C.free;
  if aPronosticos<> nil then
     aPronosticos.Free;
  if fcegh<> nil then
     fcegh.Free;
  if Beta<>nil then
     Beta.free;
  if Lambda<>nil then
       Lambda.free;

  if length(Z)<>0 then
     begin
     for k:=0 to High(Z) do
         if Z[k] <> nil then
            Z[k].Free;
     SetLength(Z,0);
     end;

  if length(M)<>0 then
     begin
     for k:=0 to High(M) do
         if M[k] <> nil then
            M[k].Free;
     SetLength(M,0);
     end;

  if length(S) <> 0 then
    begin
    for k:=0 to High(S) do
         if S[k] <> nil then
            S[k].Free;
    SetLength(S,0);
    end;

  if length(EnsamblesGauss) <> 0 then
      begin
      for k:=0 to High(EnsamblesGauss) do
           if EnsamblesGauss[k] <> nil then
              EnsamblesGauss[k].Free;
      SetLength(EnsamblesGauss,0);
      end;

  if length(Xesp) <> 0 then
      begin
      for k:=0 to High(Xesp) do
           if Xesp[k] <> nil then
              Xesp[k].Free;
      SetLength(Xesp,0);
      end;

  if length(betas) <> 0 then
     begin
     for k:=0 to High(betas) do
         if betas[k] <> nil then
            betas[k].free;
     SetLength(betas,0);
     end;

  if length(Lambdas) <> 0 then
     begin
     for k:=0 to High(Lambdas) do
         if Lambdas[k] <> nil then
            Lambdas[k].free;
     SetLength(Lambdas,0);
     end;

  if length(covar) <> 0 then
     begin
     for k:=0 to High(covar) do
         if covar[k] <> nil then
            covar[k].free;
     SetLength(covar,0);
     end;

  if length(Xg) <> 0 then
      begin
      for k:=0 to High(Xg) do
         if length(Xg[k]) <> 0 then
         begin
              for j:=0 to High(Xg[k]) do
              if Xg[k][j] <> nil then
                 Xg[k][j].Free;
              setlength(Xg[k],0);
         end;
      SetLength(Xg,0);
      end;

  if length(des) <> 0 then
      begin
      for k:=0 to High(des) do
         if length(des[k]) <> 0 then
         begin
              for j:=0 to High(des[k]) do
              if des[k][j] <> nil then
                 des[k][j].Free;
              setlength(des[k],0);
         end;
      SetLength(Xg,0);
      end;

end;

procedure TcalibrarConoEnsamble.publicarEnsemble;

// Se guardan los ensambles en la carpeta local de pronosticos y se guarda el
// ultimo pronostico accesible desde la pagina. Se almacenan 30 pronosticos y se
// eliminan los viejos
// 1- Guardar los pronosticos con la fecha del día
// 2- Se copia el ultimo pronostico a /html_vates/hidro
// 3- Se elimina los pronosticos dejasndo solo los 30 mas nuevos

var
   Nombres : array of string;
   jmiembro, hcanal, kpaso, karchi: Integer;
   archi, archi_publico,carpeta: string;
   texto: textFile;
   archivosXLT: TStringList;


begin

// 1- Guardo el ultimo pronostico

   setlength(Nombres,6);
   Nombres[0]:='Bonete';
   Nombres[1]:='Palmar';
   Nombres[2]:='SaltoGrande';
   Nombres[3]:='IN34';
   Nombres[4]:='CMO1';
   Nombres[5]:='CMO3';

   carpeta := '/home/' + GetEnvironmentVariable('USER') +
              '/prono_hidro_imfia/ensembles';

   for hcanal:=0 to high(Ensambles) do
       begin
       archi := carpeta + DirectorySeparator + Nombres[hcanal] + '_' +
                DateToStr(fechaIni.dt) + '.xlt';
       AssignFile(texto, archi);
       Rewrite(texto);

       writeln(texto, 'Nombre: ' + Nombres[hcanal] );
       writeln(texto, 'dt_Inicio: ' + DateTimeToIsoStr(fechaIni.dt) );
       writeln(texto,'Ncronicas: ' + IntToStr(Nmiembros));
       writeln(texto,'NPasos: ' + IntToStr(Npasos));
       write(texto,'dt');

       for jmiembro:=1 to Nmiembros do
           write(texto,#9 + 'member_' + IntToStr(jmiembro) );

       writeln(texto,'');

       for kpaso := 1 to Npasos do
           begin
           write(texto, FloatToStr(fechaIni.dt+kpaso/24) );
           for jmiembro:=1 to Nmiembros do
               write(texto,#9 + FloatToStr(Ensambles[hcanal].e(kpaso,jmiembro) ) );
           writeln(texto,'');
           end;
       closefile(texto);

       // 2- Guardo el ultimo en /html_vates/hidro
       archi_publico := '/home/' + GetEnvironmentVariable('USER') +
                          '/html_vates/hidro/'+ Nombres[hcanal] +'.xlt';
       cp(archi,archi_publico);

       end;


// 3- Busco los archivos si fileAge > 30 dias lo elimino
   for hcanal:=0 to high(Ensambles) do
   begin
       archivosXLT := FindAllFiles(carpeta, Nombres[hcanal] + '*' );
       archivosXLT.Sort;
       for karchi:=0 to archivosXLT.Count - 30 do
           DeleteFile(archivosXLT.Strings[karchi]);
   end;


end;

function TcalibrarConoEnsamble.ReducirProblema(MatrizOrig: TMatR;
            var indices :TDAofNInt; var L : TMatR;
            var MatrizNula :Boolean): TMatR;
{ Toma la matriz M del problema, calcula los factores que son tribiales,
indeterminados o solo se pueden resolver analiticamente. Elimina las variables
calculadas de la matriz. Devuelve la matriz reducida, un indice con el lugar de
las variables calculadas, un vector con las variables calculadas y un flag que
indica si la matriz reducida es de tamano nulo.
}
var
MatReducida, Xmat, Ymat:TMatR;
FlagM, X : array of NReal;
kfila, elementosFila, jcol, NMatReducida, indice, fila1, fila2, k:integer;
idx : array of integer;
m11, m12, m22, sumafila2: NReal;

begin
setlength(X,MatrizOrig.nf);

setlength(FlagM,MatrizOrig.nf);


for k:=0 to high(FlagM) do
    begin
    X[k] := INVERSO_DE_CERO;
    FlagM[k] := 0;
    end;

for kfila:=1 to MatrizOrig.nf do
begin


//  Reviso si la fila es nula
//------------------------------------------------------------------------------

  if MatrizOrig.Fila(kfila).normMaxAbs = 0 then
     begin
     X[kfila-1] := INVERSO_DE_CERO;
     FlagM[kfila-1] := 1;
     end
//  Reviso si la fila es diagonal
//------------------------------------------------------------------------------
  else
  if MatrizOrig.Fila(kfila).normSumAbs = abs( MatrizOrig.e(kfila,kfila) )  then
     begin
     X[kfila-1] := sqrt( 1 / MatrizOrig.e(kfila,kfila)  );
     FlagM[kfila-1] := 1;
     end

//------------------------------------------------------------------------------

//  Reviso si la fila forma un sistema de 2x2 con otra fila hacia adelante
//------------------------------------------------------------------------------
  else
    begin
    elementosFila:=0;
    fila1 := kfila;
    //  Chequeo las columnas para adelante
    for jcol:=1 to MatrizOrig.nf do
      begin
      if ( (MatrizOrig.e(kfila,jcol)<>0) and (kfila<>jcol) ) then
        begin
          elementosFila:=elementosFila+1;
          fila2 := jcol;
        end;
      end;

    if ( (elementosFila = 1) and (fila2 > fila1) ) then
    //  Entonces hay un sistema de 2x2 y todavia no se habia registrado
      begin

      sumafila2 := MatrizOrig.Fila(fila2).normSumAbs - abs(MatrizOrig.e(fila2,fila1))
                                             - abs(MatrizOrig.e(fila2,fila2)) ;

      if (sumafila2 = 0) then
        begin
        //  Calculo el lambda de kfila y columna
        m11 := MatrizOrig.e(fila1,fila1);
        m12 := MatrizOrig.e(fila1,fila2);
        m22 := MatrizOrig.e(fila2,fila2);

        if ((m22 = 0 ) or (m11 = 0)) then
          begin
          X[fila1-1] := sqrt(  1 / abs( m11 + m12));
          X[fila2-1] := sqrt(  1 / abs( m22 + m12));
          end
        else
          begin
          X[fila1-1] := sqrt( abs( 1 / ( m11 + m12 * sqrt(abs( m11 / m22) ) ) ) );
          X[fila2-1] := sqrt( abs( 1 / ( m22 + m12 * sqrt( abs ( m22 / m11) ) ) ) );
        end;


        // Marco las filas fila1 y fila2 para sacarlas
        FlagM[fila1-1] := 1;
        FlagM[fila2-1] := 1;
        end;

      end;

    end;

end;
//------------------------------------------------------------------------------
//  A este punto ya se cargaron los valores de los X tribiales y se marcaron las
//  filas/columnas a eliminar



//------------------------------------------------------------------------------
//  Elimiar las filas y columnas marcadas
NMatReducida := 0;
for k:=0 to high(FlagM) do
  NMatReducida := NMatReducida + round(FlagM[k]) ;

NMatReducida := MatrizOrig.nf - NMatReducida;

if NMatReducida > 0 then
begin
  MatReducida := TMatR.Create_InitVal(NMatReducida,NMatReducida,0);
  MatrizNula := false;

  setlength(idx,NMatReducida);
  indice:=0;
  for k:=1 to MatrizOrig.nf do
    begin
      if FlagM[k-1] = 0 then
        begin
        idx[indice] := k;
        indice := indice + 1;
        end;
    end;

  for kfila := 1 to NMatReducida do
    begin
    for jcol := 1 to NMatReducida do
      begin
      MatReducida.pon_e(kfila, jcol, MatrizOrig.e( idx[kfila-1], idx[jcol-1] ) );
      end;
    end;

  end
  else
    begin
    MatrizNula := true;
    setlength(idx,0);
    end;
//------------  Armo los vectores Beta y Lambda
  Xmat := TMatR.Create_Init(Ncanales,1);
  Xmat.pon_columna(1,1,X) ;
  Ymat := TMatR.Create_InitVal(Ncanales,1,0.1);
  for k:=1 to Ncanales do
      Ymat.pon_e(k,1,1/Xmat.e(k,1));

//------------ Cargo la salida
  indices := idx;
  L:=Ymat;

  result := MatReducida;

end;

procedure TcalibrarConoEnsamble.guardarSala(sala: TSalaDeJuego);
// Guarda la sala con los cambios que se hicieron durante la iteracion actual

var
  f: TArchiTexto;
  dirTemporales: array [0 .. MAX_PATH] of char;
  archi: String;

begin

  dirTemporales := GetTempDir;

  f := TArchiTexto.CreateForWrite(dirTemporales + 'guardar_temp.ese', True, 1);
  try
    begin
      f.wr('sala', TSalaDeJuego(sala));
      f.Free;

      archi := sala.archisala_;
      f := TArchiTexto.CreateForWrite(archi, True, 1);
      f.wr('sala', TSalaDeJuego(sala));

      sala.setDirCorrida(archi);
    end
  finally
    f.Free;
  end;
end;





end.

