{$DEFINE PERMITIR_DEFORMADORES_NULOS}
unit uanalisisserial;

interface

uses
  SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
  xmatdefs,
  fddp, fddp_weibull,
  ufechas,
  matreal, matent,
  uhermitaprox,
  uregresion,
  StdCtrls, udisnormcan,
  uAuxiliares,
  uverdoc,
  useriestemporales,
  umodelosintcegh,
  wrffti01, wrfftf01,
  Math, uestimadoryrellenado,
  uProbabilidadTransicion,
  uradiacionsolar,
  ucurvasaeros,
  uopencalc,
  ExtCtrls, ComCtrls, uhurst, autoco, Types;

type

  { TForm1 }

  TForm1 = class(TForm)
    BAnalizar: TButton;
    BBuscar: TButton;
    BCalcularPromediosMensuales: TButton;
    BLeerDatosBasicos: TButton;
    BPredictorNino34: TButton;
    btAutoco: TButton;
    btCargarCEGH: TButton;
    btClustering: TButton;
    btDestendenciar: TButton;
    btCompletarSeries: TButton;
    BTestVarianza: TButton;
    btGenerarSeries: TButton;
    btH2f: TButton;
    btHurst: TButton;
    btSelSeries: TButton;
    btWeibulizar: TButton;
    btAyuda_principlal: TButton;
    BtPotVelAeroGen: TButton;
    btCalcularErrores: TButton;
    Button1: TButton;
    btCrearSAS: TButton;
    Button2: TButton;
    btModeladoAportes: TButton;
    cbAniosAnalogos: TCheckBox;
    cbCopiarUltimomCicloComoEstadoinicial: TCheckBox;
    cbFiltradoKTPorIh0: TCheckBox;
    cbFiltrarMenores: TCheckBox;
    cbFiltroCiclico: TCheckBox;
    cb_FiltrarAlGaussianizar: TCheckBox;
    cbCurvaAero: TComboBox;
    cb_Filtrarmenores_CalcErrores: TCheckBox;
    cbTendenciaExponencial: TCheckBox;
    cb_InfoParcial: TCheckBox;
    cbCompletarHuecosConRuidoBlanco: TCheckBox;
    eArchiSeriesDeDatos: TEdit;
    eMinicicloCaracteristico: TLabeledEdit;
    eMinicicloCaracteristicoTendencia: TLabeledEdit;
    eNClusters: TLabeledEdit;
    eGradoTendencia: TLabeledEdit;
    eNCronicas: TLabeledEdit;
    eNDecimalesDesformadores: TEdit;
    eNDigitosDesformadores: TEdit;
    eNPuntosPorMiniCiclo: TEdit;
    eNP_Deformadores: TLabeledEdit;
    eOrdenDelFiltro: TEdit;
    eOrdenESHY: TEdit;
    eOverlapping: TEdit;
    eSemanaInicioVentanaAnaloga: TEdit;
    ESHY: TButton;
    eTraslapping: TEdit;
    eUmbralAlGaussianizar: TEdit;
    eUmbralFiltro: TEdit;
    eVarianzaDeLoNoExplicado: TEdit;
    gbClustering: TGroupBox;
    gbDestendenciador: TGroupBox;
    gb_GapFilling: TGroupBox;
    GenerarRuidoRectificado: TButton;
    GroupBox1: TGroupBox;
    GroupBox2: TGroupBox;
    gbNeuroGaussian: TGroupBox;
    GroupBox3: TGroupBox;
    GroupBox4: TGroupBox;
    Label1: TLabel;
    Label10: TLabel;
    eUmral_calcerrores: TLabeledEdit;
    eNRetardos: TLabeledEdit;
    eNNeuronasPorCanal: TLabeledEdit;
    eNCronicasRuido: TLabeledEdit;
    ePrecisionMedida: TLabeledEdit;
    lbl_ModeloAerogen: TLabel;
    Label2: TLabel;
    Label3: TLabel;
    Label4: TLabel;
    Label5: TLabel;
    Label6: TLabel;
    Label7: TLabel;
    Label8: TLabel;
    Label9: TLabel;
    LArchiSeries: TLabel;
    LOrdenFiltro: TLabel;
    lstAniosAnalogos: TMemo;
    memo_model: TMemo;
    OpenDialog1: TOpenDialog;
    PageControl1: TPageControl;
    P_Transicion: TButton;
    rgTipoModelo: TRadioGroup;
    rbg_calculoB: TRadioGroup;
    Splitter1: TSplitter;
    TestDistribucionGaussiana: TButton;
    tsAniosAnalogos: TTabSheet;
    tsCEGH: TTabSheet;
    tsExperimental: TTabSheet;
    tsModelador: TTabSheet;
    tsRellenado: TTabSheet;


    procedure BAnalizarClick(Sender: TObject);
    procedure btCalcularErroresClick(Sender: TObject);
    procedure btCargarCEGHClick(Sender: TObject);
    procedure btCompletarSeriesClick(Sender: TObject);
    procedure btCrearSASClick(Sender: TObject);
    procedure btGenerarSeriesClick(Sender: TObject);
    procedure btModeladoAportesClick(Sender: TObject);
    procedure btSelSeriesClick(Sender: TObject);
    procedure btWeibulizarClick(Sender: TObject);
    procedure BtPotVelAeroGenClick(Sender: TObject);
    procedure btClusteringClick(Sender: TObject);
    procedure Button1Click(Sender: TObject);
    procedure Button2Click(Sender: TObject);
    procedure FormCreate(Sender: TObject);
    procedure memo_modelChange(Sender: TObject);
    procedure P_TransicionClick(Sender: TObject);
    procedure ESHYClick(Sender: TObject);
    procedure rgTipoModeloClick(Sender: TObject);
    procedure TestDistribucionGaussianaClick(Sender: TObject);
    procedure GenerarRuidoRectificadoClick(Sender: TObject);
    procedure BBuscarClick(Sender: TObject);
    procedure BLeerDatosBasicosClick(Sender: TObject);
    procedure BPredictorNino34Click(Sender: TObject);
    procedure BTestVarianzaClick(Sender: TObject);
    procedure BCalcularPromediosMensualesClick(Sender: TObject);
    procedure btAyuda_principlalClick(Sender: TObject);
    procedure btHurstClick(Sender: TObject);
    procedure btAutocoClick(Sender: TObject);
    procedure btH2fClick(Sender: TObject);

  private
    seriesDeDatos: TSeriesDeDatos;
    //seriesDeDatos1:  TSeriesDeDatos;
    aniosanalogos: array of integer;
    semanaDeInicio: integer;
    kBuscarDesde: integer;

  public
    { Public declarations }
    cegh: TModeloCEGH;

    curvasAeros_lst: TList;

    procedure GetAniosAnalogos;
    procedure MascararAniosNoAnalogos(datos: TSeriesDeDatos; mascara: double);
    function Mascarado(anio, semana: integer): boolean;
    function buscarAnio(anio: integer): integer;

  end;

var
  Form1: TForm1;


implementation

{$R *.lfm}


// recorre las series y quita los tramos en que halla actuado el filtro
// acortando las series.
procedure AcortarTramosFiltrados(var series: TDAOfVectR; umbralFiltro: NReal);

var
  j, jval, NSems: integer;
  flg_Filtrado: boolean;
  iserie: integer;

begin
  NSems := Series[0].n;
  j := 1;
  while (j <= nSems) do
  begin
    //detecto si actuó el filtro en alguna de las series.
    flg_Filtrado := False;
    for iserie := 0 to high(series) do
      if series[iserie].e(j) <= umbralFiltro then
      begin
        flg_Filtrado := True;
        break;
      end;

    // si actuó el filtro acortamos todas las series.
    if flg_Filtrado then
    begin
      for iserie := 0 to high(series) do
      begin
        for jval := j to nSems - 1 do
          series[iserie].pon_e(jval, series[iserie].e(jval + 1));
      end;
      Dec(nSems);
    end
    else
      Inc(j);
  end;

  // si corresponde acortamos las series.
  if nSems < series[0].n then
  begin
    for iserie := 0 to high(series) do
    begin
      setlength(series[iserie].pv, nSems + 1);
      series[iserie].n := nSems;
    end;
  end;
end;


// retorna la matriz B triangular Inferior de la factorización de
// Cholesky que tal que B*transpuesta(B) es la matriz de covarianzas
// de las series.
function descuadrador_Cholesky(var series: TDAOfVectR; filtrar: boolean;
  umbralFiltro: NReal): TMatR;

var
  A: TMatR; // matriz de covarianzas
  B: TMatR;
  k, j: integer;
  m: NReal;
  NDatos: integer;
  DimRes: Integer;

begin

  // SI la orden es filtrar, acortamos las series
  // quitando los tramos en los que actuó el filtro.
  if Filtrar then
    AcortarTramosFiltrados(series, umbralFiltro);


  NDatos := series[0].n;

  A := TMatR.Create_Init(length(series), length(series));
  // construimos la matriz de covarianzas
  // como es simétrica lleno solo el triangulo superior
  for k := 1 to A.nf do
    for j := k to A.nc do
    begin
      // se supone que las series tienen valor esperado NULO
      m := series[k - 1].PEV(series[j - 1]) / NDatos;
      A.pon_e(k, j, m);
    end;

  B := A.Raiz_Cholesky;
  if B = nil then
    B := A.RaizPorPotenciaIterada( DimRes, true );
  A.Free;
  Result := B;
end;

function descuadrador(var series: TDAOfVectR; filtrar: boolean;
  umbralFiltro: NReal): TMatR;
const
  MaxNIters = 1000;

var
  NSems, NSeries, NDescuadres: integer;

var
  x_series, x_autovects: TVectR;
  xs_series, xs_autovects: TVectR;
  autovects: TDAOfVectR;
  autovals: array of double;
  convergio: array of boolean;
  kvect: integer;
  NAutovects: integer;
  resB: TMatR;

var
  sal{, f}: textfile;
  //  r: string;
  //  ksem: integer;
  //  ap_b,ap_p,ap_s: double;

  k, j: integer;
  me2: double;

  buscando: boolean;

  descuadrando: boolean;
  niter: integer;
  factorCambioNorma: NReal;

  iserie: integer;

  varianzas: TVectR;

  sb2: NReal;

begin
  (*** OJO NO USE ESTE METODO ***)


  // SI la orden es filtrar, acortamos las series
  // quitando los tramos en los que actuó el filtro.
  if Filtrar then
    AcortarTramosFiltrados(series, umbralFiltro);

  // más vale hacer esto después de acortar las series.
  NSeries := length(Series);
  NSems := Series[0].n;

  // calculamos las varianzas
  varianzas := TVectR.Create_Init(NSeries);
  for iserie := 1 to NSeries do
    varianzas.pv[iserie] := series[iserie - 1].varianza;

  NDescuadres := NSeries;

  setlength(autovects, NDescuadres);
  setlength(autovals, NDescuadres);
  setlength(convergio, NDescuadres);


  for k := 1 to NDescuadres do
    autovects[k - 1] := TVectR.Create_Init(NSems);

  x_series := TVectR.Create_Init(NSeries);
  x_autovects := TVectR.Create_Init(NDescuadres);
  xs_series := TVectR.Create_Init(NSeries);
  xs_autovects := TVectR.Create_Init(NDescuadres);

  kvect := 1;
  descuadrando := True;

  NAutovects := NDescuadres;

  while (kvect <= NDescuadres) and (descuadrando) do
  begin
    // cargamos un versor cualquiera como candidato
    // a ser el autovalor principal del conjunto de series.
    // elegimos la serie como inicio de iteración, pero no es imperativo.
    autovects[kvect - 1].Copy(Series[kvect - 1]);
    autovects[kvect - 1].hacerunitario;

    // inicializamos variables auxiliares de la iteración
    x_series.Ceros;
    x_autovects.Ceros;
    buscando := True;
    niter := 1;


    while buscando and (niter <= MaxNIters) do
    begin
      // proyecciones sobre las series
      for k := 1 to NSeries do
        xs_series.pon_e(k, series[k - 1].PEV(autovects[kvect - 1]));
      // proyecciones sobre autovectores ya estimados
      for k := 1 to kvect - 1 do
        xs_autovects.pon_e(k, autovects[k - 1].PEV(autovects[kvect - 1]));

      // calculo de variación para dar fin a la iteración
      me2 := (xs_series.distancia2(x_series) + xs_autovects.distancia2(
        x_autovects)) / NSems;
      if (me2 < 1e-14) then
        buscando := False;
      vswap(xs_series, x_series);
      vswap(xs_autovects, x_autovects);
      // formación del vector
      autovects[kvect - 1].ceros;
      for k := 1 to NSeries do
        autovects[kvect - 1].sumRPV(x_series.e(k), series[k - 1]);

      for k := 1 to kvect - 1 do
        autovects[kvect - 1].sumRPV(-autovals[k - 1] * x_autovects.e(k),
          autovects[k - 1]);

      // usamos la raiz de la varianza como distancia
      autovals[kvect - 1] := autovects[kvect - 1].normEuclid;
      // lo hacemos unitario
      autovects[kvect - 1].PorReal(1 / autovals[kvect - 1]);
      Inc(niter);
    end;
    if buscando then
    begin
      descuadrando := False;
      NAutovects := kvect - 1;
    end
    else
    begin
      writeln('me2: ', me2, ' landa: ', autovals[kvect - 1]);
      writeln('niter: ', niter);
      Inc(kvect);
    end;
  end;


(* Escalamos los autovectores para que sean de varianza 1
en lugar de Norma Euclidea 1 *)
  factorCambioNorma := sqrt(NSems);
  for j := 1 to NAutovects do
    autovects[j - 1].PorReal(factorCambioNorma);

(* generalmente NAutovects = NSeries, sería demasiada casualidad
que las muestras dieran con un autovalor NULO *)
  resB := TMatR.Create_Init(NSeries, NAutovects);
  for k := 1 to NSeries do
    for j := 1 to NAutovects do
      resB.pon_e(k, j, series[k - 1].PEV(autovects[j - 1]) / NSems);
(*
 las filas de B son las cordenadas de las series en la base de los
 autovectores.*)
  Result := resB;

  for k := 1 to NSeries do
  begin
    sb2 := resB.Fila(k).ne2;
    resB.Fila(k).PorReal(sqrt(varianzas.e(k) / sb2));
  end;

  (*** IMPRIMIMOS el archivo de control de resultados ****)
  assignfile(sal, 'descuadre.xlt');
  rewrite(sal);
  writeln(sal, 'MATRIZ B ');
  for k := 1 to NSeries do
  begin
    Write(sal, 'serie: ', #9, k);
    for j := 1 to NAutovects do
      Write(sal, #9, resB.e(k, j));
    writelN(sal);
  end;

  writeln(sal);
  Write(sal, 'autovalores:', #9);
  for j := 1 to NAutovects do
    Write(sal, #9, autovals[j - 1]);
  writeln(sal);


  writeln(sal, ' SERIES y AUTOVECTORES ');
  for j := 1 to NSeries do
    Write(sal, #9, 's', j);
  for j := 1 to NAutovects do
    Write(sal, #9, 'v', j);
  writeln(sal);

  for k := 1 to NSems do
  begin
    Write(sal, k);
    for j := 1 to NSeries do
      Write(sal, #9, series[j - 1].e(k));
    for j := 1 to NAutovects do
      Write(sal, #9, autovects[j - 1].e(k));
    writeln(sal);
  end;
  closefile(sal);

end;




function ivertirfiltro(filtro: TVectR; series: array of TVectR;
  kSerie: integer; filtrar: boolean; umbralFiltro: NReal;
  var varianzaDelResultado: NReal; var cntSecuenciasUsadas: integer): TVectR;
var
  k: integer;
  orden: integer;
  res: TVectR;
  a, av: NReal;
  h: integer;
  jSerie: integer;
  NSeries: integer;
  NPuntos: integer;
  flgFiltrado: boolean;

begin
  NSeries := length(series);
  orden := filtro.n div NSeries;
  NPuntos := series[0].n;
  res := TVectR.Create_Init(NPuntos - orden);

  varianzaDelResultado := 0;
  cntSecuenciasUsadas := 0;

  if not Filtrar then
  begin
    for k := 1 + orden to NPuntos do
    begin
      a := series[kSerie - 1].pv[k];
      for jSerie := 1 to NSeries do
        for h := 1 to orden do
          if (k - h) > 0 then
            a := a - filtro.e((h - 1) * NSeries + jSerie) * series[jSerie - 1].pv[k - h];
      res.pv[k - orden] := a;
      varianzaDelResultado := a * a + varianzaDelResultado;
      Inc(cntSecuenciasUsadas);
    end;
  end
  else
  begin
    for k := 1 + orden to NPuntos do
    begin
      a := series[kSerie - 1].pv[k];
      flgFiltrado := a <= umbralFiltro;
      for jSerie := 1 to NSeries do
        for h := 1 to orden do
          if ((k - h) > 0) and (not FlgFiltrado) then
          begin
            av := series[jSerie - 1].pv[k - h];
            if av > umbralFiltro then
              a := a - filtro.e((h - 1) * NSeries + jSerie) * av
            else
              flgFiltrado := True;
          end;

      if flgFiltrado then
        res.pv[k - orden] := umbralFiltro
      else
      begin
        res.pv[k - orden] := a;
        varianzaDelResultado := a * a + varianzaDelResultado;
        Inc(cntSecuenciasUsadas);
      end;
    end;
  end;

  varianzaDelResultado := varianzaDelResultado / (cntSecuenciasUsadas - 1);
  Result := res;
end;


procedure Write_SeriesGaussianas(seriesDeDatos: TSeriesDeDatos);
begin
  seriesDeDatos.WriteToArchi('series_gaussianas.xlt');
end;



procedure Write_FFTSeriesGaussianas(seriesDeDatos: TSeriesDeDatos);
var
  sal: textfile;
  kSerie, kPunto: integer;
  coA, coB: TDAOfNreal;
  A0: NReal;
begin
  setlength(coA, seriesDeDatos.NPuntos div 2);
  setlength(coB, seriesDeDatos.NPuntos div 2);
  WRFFTI01.Init(seriesDeDatos.NPuntos);
  for kSerie := 1 to seriesDeDatos.NSeries do
  begin
    WRFFTF01.FFTF(seriesDeDatos.Series[kSerie - 1].pv[1], A0, coA[0], coB[0]);
    assignfile(sal, 'FFT_seriegaussiana_' + IntToStr(kSerie) + '.txt');
    rewrite(sal);
    writeln(sal, 'A0: ', A0);
    writeln(sal, 'k', #9, 'CA', #9, 'CB');
    for kPunto := 0 to high(coA) do
      writeln(sal, kPunto + 1, #9, coA[kPunto], #9, coB[kPunto]);
    closefile(sal);
  end;
  WRFFTI01.done;
  setlength(coA, 0);
  setlength(coB, 0);
end;




// calcula los deformadores y el filtro óptimo en el espacio gaussiano.
// los lentes son por paso mientras que el filtro es único.
// La versión Analizar2 identifica un filtro para cada paso.
function Analizar(archi: string; // ARchivo de salida. 'SintetizadorCEGH.txt'
  seriesDeDatos: TSeriesDeDatos;
  NCronicasRuido: integer;
  PrecisionMedida_pu: NReal;
  ordenfiltro: integer; overlapping: integer;
  traslapping: integer; nPuntosPorMiniciclo: integer; horasEntreMuestras: integer;
  FiltrarMenores: boolean; umbralFiltro: NReal; FiltrarAlGaussianizar: boolean;
  umbralFiltroAlGaussianizar: NReal; kMetodoCalculoB: integer;
  NDigitosDeformadores, NDecimalesDeformadores: integer;
  // formato escritura
  NP_Deformadores: integer; // Cantidad de puntos de los deformadores.
  flg_InfoParcial: boolean; anular_negativos_hasta: integer =
  0 // ej. si = 4 anula los deformadores donde sean negativos para series 1,..4
  ): TVectR;

var
  //  r: string;
  ruidos: TDAOfVectR;
  kserie: integer;
  k, j, h: integer;
  //  sal: TextFile;
  sintetizadorSimSEE: TModeloCEGH;
  nVEred: integer;
  varNoExplicada: NReal;
  d_varNoExplicada: NReal;
  cntSecuenciasUsadas: integer;
  At, BtB: TMatR;
  res: TVectR;
  filas: integer;
  Ae: TMatR;
  Ap: TMatR;
  f: integer;
  c: integer;
  flg_CorrCorPla: boolean;
  o: integer;
  nseries: integer;
  Vdif: TVectR;
  Ve: TVectR;
  vpc: TVectR;
  Mcegh: TModeloCEGH;
  series: TSeriesDeDatos;
  kpaso: integer;
  filtrar_huecos: boolean;
  valmin: integer;
  idxhuecos: TDAofNInt;
  ek: TMatR;
  Xpk: TMatR;
  archiDatos: string;
  idx: integer;
  prod: extended;
  archiSalidaGauss: string;
  archiSalida: string;
  DimRes: integer;
begin

  nVEred := seriesDeDatos.NSeries * ordenfiltro;


  // creamos el sintetizadorSimSEE
  sintetizadorSimSEE := TModeloCEGH.Create(seriesDeDatos.NSeries,
    seriesDeDatos.nombresSeries, ordenfiltro, trunc(seriesDeDatos.rNPPorCiclo),
    NP_Deformadores, horasEntreMuestras, nVEred);

  // **********************************************
  // PASO 1) Creamos deformadores y gaussianizamos
  // ----------------------------------------------

  // Creamos los deformadores
  sintetizadorSimSEE.CrearDeformadores(seriesDeDatos,
    overlapping, traslapping, NPuntosPorMiniciclo,
    FiltrarMenores, umbralFiltro, NCronicasRuido, PrecisionMedida_pu );

  // Transformamos las series al espacio gaussiano
  sintetizadorSimSEE.GaussianizarSeries(
    seriesDeDatos,
    FiltrarMenores,
    umbralFiltro,
    filtrarAlGaussianizar,
    umbralFiltroAlGaussianizar, 1, 0 );

  // Antes de identificar el filtro, normalizamos las series por las dudas.
  sintetizadorSimSEE.Normalizar(seriesDeDatos, filtrarMenores, umbralFiltro);

  // Escribimos las series gaussianas
  write_SeriesGaussianas(seriesDeDatos);

  // Escribimos archivo con covarianzas(tau)
  seriesDeDatos.WriteCOVARS('covars_k_gaussianas.xlt', umbralFiltro);

  // Escribimos FFT de las series gaussianas
  Write_FFTSeriesGaussianas(seriesDeDatos);

  // ================================================




  (**************************************************************)
  (* Identificación del filtro lineal (en el espacio gaussiano) *)
  (*------------------------------------------------------------*)
  // Primero identificamos la matriz A del filtro por mínimos cuadrados
  At := calcmincuad_mrh(seriesDeDatos.get_series_cronActiva, ordenfiltro, filtrarMenores,
    umbralFiltro, flg_InfoParcial, BtB);
  SintetizadorSimSEE.A_cte := At.Crear_transpuesta;
  At.Free;


  setlength(ruidos, seriesDeDatos.NSeries);
  // inversión de las series por el filtro
  // para tener las gaussianas de entrada
  // El largo de las series de ruidos es NPuntos-Orden pues
  // para los primeros puntos de las series (de 1 a Orden) no hay
  // datos para calcular el ruido.

    (*
    varNoExplicada := 0;
    for kserie := 1 to seriesDeDatos.NSeries do
    begin
      ruidos[kserie - 1] := ivertirfiltro(SintetizadorSimSee.A_cte.fila(kserie),
        seriesDeDatos.series, kserie, filtrarMenores, umbralFiltro,
        d_varNoExplicada, cntSecuenciasUsadas);
      varNoExplicada     := varNoExplicada + d_varNoExplicada / seriesDeDatos.NSeries;
    end;
    *)


  res := TVectR.Create_INit(seriesDeDatos.NSeries);
  for kSerie := 1 to res.n do
    res.pon_e(kSerie, BtB.e(kSerie, kSerie));

  varNoExplicada := BtB.traza / seriesDeDatos.NSeries;

  if kMetodoCalculoB = 0 then
  begin
    BtB.WriteArchiXLT('btb.xlt');
    SintetizadorSimSee.B_cte := BtB.raiz_Cholesky;
    if SintetizadorSimSee.B_cte = nil then
      SintetizadorSimSee.B_cte := BtB.RaizPorPotenciaIterada( DimRes, true );
  end
  else
  begin
    // OJO, como efecto secundario el dscuadrador puede acortar los vectores
    // de ruidos, eliminando los que fueron filtrados.
    if kMetodoCalculoB = 2 then
      SintetizadorSimSee.B_cte :=
        descuadrador_Cholesky(ruidos, FiltrarMenores, umbralFiltro)
    else
      SintetizadorSimSee.B_cte := descuadrador(ruidos, FiltrarMenores, umbralFiltro);
  end;
  BtB.Free;


  // bien ahora vamos a completar la información de reducción de estados.
  for j := 1 to ordenfiltro do
    for k := 1 to seriesDeDatos.NSeries do
    begin

      if ordenfiltro > 1 then
        SintetizadorSimSEE.nombreVarE[(k - 1) + (j - 1) * SeriesDeDatos.NSeries] :=
          seriesDeDatos.nombresSeries[k - 1] + '_' + IntToStr(j)
      else
        SintetizadorSimSEE.nombreVarE[(k - 1) + (j - 1) * SeriesDeDatos.NSeries] :=
          seriesDeDatos.nombresSeries[k - 1];

      //SintetizadorSimSEE.MRed.pon_e(k, (k - 1) * ordenfiltro + j, 1);
      SintetizadorSimSEE.MRed_cte.pon_e(k + (j - 1) * seriesDeDatos.NSeries,
        k + (j - 1) * seriesDeDatos.NSeries, 1);
      SintetizadorSimSEE.nDiscsVsE[(k - 1) * ordenfiltro + j - 1] := 5;
      setlength(SintetizadorSimSEE.ProbsVsE[(k - 1) * ordenfiltro + j - 1], 5);
      for h := 0 to 4 do
        SintetizadorSimSEE.ProbsVsE[(k - 1) * ordenfiltro + j - 1][h] := 0.2;
    end;

  if anular_negativos_hasta > 0 then
    for k := 1 to anular_negativos_hasta do
      SintetizadorSimSEE.Anular_Negativos_DelDeformador(k);

  SintetizadorSimSEE.WriteToArchi(archi, NDigitosDeformadores,
    NDecimalesDeformadores);
  //   Result := VarNoExplicada;
  Result := res;
end;




// calcula los deformadores y el filtro óptimo en el espacio gaussiano.
// los lentes son por paso mientras que el filtro es único.
// Genera modelo del tipo Y[k] = sum( A[i] X[k-i] + B R[k] ; i = 0, ..., NRetardos )
function Analizar_Y_AX_B(archi: string; // ARchivo de salida. 'SintetizadorCEGH.txt'
  seriesDeDatos: TSeriesDeDatos;
  NCronicasRuido: integer;
  PrecisionMedida_pu: NReal;
  ordenfiltro: integer; overlapping: integer;
  traslapping: integer; nPuntosPorMiniciclo: integer; horasEntreMuestras: integer;
  FiltrarMenores: boolean; umbralFiltro: NReal; FiltrarAlGaussianizar: boolean;
  umbralFiltroAlGaussianizar: NReal; kMetodoCalculoB: integer;
  NDigitosDeformadores, NDecimalesDeformadores: integer;
  // formato escritura
  NP_Deformadores: integer; // Cantidad de puntos de los deformadores.
  flg_InforParcial: boolean; anular_negativos_hasta: integer =
  0 // ej. si = 4 anula los deformadores donde sean negativos para series 1,..4
  ): TVectR;

var
  //  r: string;
  ruidos: TDAOfVectR;
  kserie: integer;
  k, j, h: integer;
  //  sal: TextFile;
  {$IFDEF CEGH_ENTRADA_SALIDA}
  sintetizadorSimSEE: TModeloCEGH_EntradaSalida;
  {$ELSE}
  sintetizadorSimSEE: TModeloCEGH;
  {$ENDIF}
  nVEred: integer;
  varNoExplicada: NReal;
  d_varNoExplicada: NReal;
  At: TMatR;
  BtB: TMatR;
  res: TVectR;
  DimRes: integer;

begin

  if not flg_InforParcial then
    raise Exception.Create('este modelo está disponible solo con INFORMACION PARCIAL');

  nVEred := seriesDeDatos.NSeries * ordenfiltro;


  // creamos el sintetizadorSimSEE
  {$IFDEF CEGH_ENTRADA_SALIDA}
  sintetizadorSimSEE := TModeloCEGH_EntradaSalida.Create(
    seriesDeDatos.NSeriesX, seriesDeDatos.NSeriesY, seriesDeDatos.nombresSeries,
    ordenfiltro, trunc(seriesDeDatos.rNPPorCiclo), NP_Deformadores,
    horasEntreMuestras, nVEred);
  {$ELSE}
  sintetizadorSimSEE := TModeloCEGH.Create(seriesDeDatos.NSeries,
    seriesDeDatos.nombresSeries, ordenfiltro, trunc(seriesDeDatos.rNPPorCiclo),
    NP_Deformadores, horasEntreMuestras, nVEred);
  {$ENDIF}


  // **********************************************
  // PASO 1) Creamos deformadores y gaussianizamos
  // ----------------------------------------------

  // Creamos los deformadores
  sintetizadorSimSEE.CrearDeformadores(seriesDeDatos,
    overlapping, traslapping, NPuntosPorMiniciclo,
    FiltrarMenores, umbralFiltro, NCronicasRuido, PrecisionMedida_pu );

  // Transformamos las series al espacio gaussiano
  sintetizadorSimSEE.GaussianizarSeries(
    seriesDeDatos,
    FiltrarMenores,
    umbralFiltro,
    filtrarAlGaussianizar,
    umbralFiltroAlGaussianizar, 1, 0);

  // Antes de identificar el filtro, normalizamos las series por las dudas.
  sintetizadorSimSEE.Normalizar(seriesDeDatos, filtrarMenores, umbralFiltro);

  // Escribimos las series gaussianas
  write_SeriesGaussianas(seriesDeDatos);

  // Escribimos archivo con covarianzas(tau)
  seriesDeDatos.WriteCOVARS('covars_k_gaussianas.xlt', umbralFiltro);

  // Escribimos FFT de las series gaussianas
  //  Write_FFTSeriesGaussianas( seriesDeDatos );

  // ================================================


  (**************************************************************)
  (* Identificación del filtro lineal (en el espacio gaussiano) *)
  (*------------------------------------------------------------*)
  // Primero identificamos la matriz A del filtro por mínimos cuadrados
  At := calcmincuad_mrh_Y_AX_B(seriesDeDatos.get_seriesY_cronActiva, seriesDeDatos.get_seriesX_cronActiva,
    ordenfiltro - 1, filtrarMenores, umbralFiltro, BtB);
  SintetizadorSimSEE.A_cte := At.Crear_Transpuesta;
  At.Free;


  res := TVectR.Create_INit(seriesDeDatos.NSeriesY);
  for kSerie := 1 to res.n do
    res.pon_e(kSerie, BtB.e(kSerie, kSerie));

  varNoExplicada := BtB.traza / seriesDeDatos.NSeriesY;

  if kMetodoCalculoB = 0 then
  begin
    BtB.WriteArchiXLT('btb.xlt');
    SintetizadorSimSee.B_cte := BtB.raiz_Cholesky;
    if SintetizadorSimSee.B_cte = nil then
      SintetizadorSimSee.B_cte := BtB.RaizPorPotenciaIterada( DimRes, true );
  end
  else
  begin
    // OJO, como efecto secundario el dscuadrador puede acortar los vectores
    // de ruidos, eliminando los que fueron filtrados.
    if kMetodoCalculoB = 2 then
      SintetizadorSimSee.B_cte :=
        descuadrador_Cholesky(ruidos, FiltrarMenores, umbralFiltro)
    else
      SintetizadorSimSee.B_cte := descuadrador(ruidos, FiltrarMenores, umbralFiltro);
  end;
  BtB.Free;


  // bien ahora vamos a completar la información de reducción de estados.
  for j := 1 to ordenfiltro do
    for k := 1 to seriesDeDatos.NSeriesX do
    begin

      if ordenfiltro > 1 then
        SintetizadorSimSEE.nombreVarE[(k - 1) + (j - 1) * SeriesDeDatos.NSeriesX] :=
          seriesDeDatos.nombresSeriesX[k - 1] + '_' + IntToStr(j)
      else
        SintetizadorSimSEE.nombreVarE[(k - 1) + (j - 1) * SeriesDeDatos.NSeriesX] :=
          seriesDeDatos.nombresSeriesX[k - 1];

      SintetizadorSimSEE.MRed_cte.pon_e(k, (k - 1) * ordenfiltro + j, 1);
      SintetizadorSimSEE.nDiscsVsE[(k - 1) * ordenfiltro + j - 1] := 5;
      setlength(SintetizadorSimSEE.ProbsVsE[(k - 1) * ordenfiltro + j - 1], 5);
      for h := 0 to 4 do
        SintetizadorSimSEE.ProbsVsE[(k - 1) * ordenfiltro + j - 1][h] := 0.2;
    end;

  if anular_negativos_hasta > 0 then
    for k := 1 to anular_negativos_hasta do
      SintetizadorSimSEE.Anular_Negativos_DelDeformador(k);

  SintetizadorSimSEE.WriteToArchi(archi, NDigitosDeformadores,
    NDecimalesDeformadores);
  //   Result := VarNoExplicada;
  Result := res;

end;



// construye los deformadores (uno por paso) y los filtros óptimos
// para cada paso en el espacio gaussiano.
function Analizar_FiltroVariable(archi: string;
  // Archivo de salida 'SintetizadorCEGH_mc.txt'
  seriesDeDatos: TSeriesDeDatos;
  NCronicasRuido: integer;
  PrecisionMedida_pu: NReal;
  ordenfiltro: integer; overlapping: integer;
  traslapping: integer; NPuntosPorMiniciclo: integer; horasEntreMuestras: integer;
  FiltrarMenores: boolean; umbralFiltro: NReal; FiltrarAlGaussianizar: boolean;
  umbralFiltroAlGaussianizar: NReal;
  NDigitosDeformadores, NDecimalesDeformadores: integer; // formato escritura
  NP_Deformadores: integer;   // cantidad de puntos usados en la discretización de los deformadores
  flg_InfoParcial: boolean;
  flg_CompletarHuecosConRuidoBlanco: boolean ): TVectR;

var
  k, j, h: integer;
  //  sal: TextFile;

  nVEred: integer;

  sintetizadorSimSEE: TModeloCEGH;

  mcA, mcB: TDAOfMatR;

  fsal: textfile;
  kfil, jcol: integer;
  varNoExplicado: NReal;
  res: TVectR;

  nFilas, nRetardo: integer;
  kSerie: integer;

begin

  nVEred := seriesDeDatos.NSeries * ordenfiltro;

  // creamos el sintetizadorSimSEE
  sintetizadorSimSEE := TModeloCEGH.Create_MultiCiclo(
    seriesDeDatos.NSeries, seriesDeDatos.nombresSeries, ordenfiltro,
    trunc(seriesDeDatos.rNPPorCiclo), NP_Deformadores, horasEntreMuestras, nVERed);


  // **********************************************
  // PASO 1) Creamos deformadores y gaussianizamos
  // ----------------------------------------------

  // Creamos los deformadores
  sintetizadorSimSEE.CrearDeformadores(seriesDeDatos,
    overlapping, traslapping, NPuntosPorMiniciclo,
    FiltrarMenores, umbralFiltro, NCronicasRuido, PrecisionMedida_pu );

  // Transformamos las series al espacio gaussiano
  sintetizadorSimSEE.GaussianizarSeries(
    seriesDeDatos,
    FiltrarMenores,
    umbralFiltro,
    filtrarAlGaussianizar,
    umbralFiltroAlGaussianizar, 1, 0);

  // Antes de identificar el filtro, normalizamos las series por las dudas.
  sintetizadorSimSEE.Normalizar(seriesDeDatos, filtrarMenores, umbralFiltro);

  // Escribimos las series gaussianas
  write_SeriesGaussianas(seriesDeDatos);

  // Escribimos archivo con covarianzas(tau)
  seriesDeDatos.WriteCOVARS('covars_k_gaussianas.xlt', umbralFiltro);

  // Escribimos FFT de las series gaussianas
  //  Write_FFTSeriesGaussianas( seriesDeDatos );

  // ================================================


  (**************************************************************)
  (* Identificación del filtro lineal (en el espacio gaussiano) *)
  (*------------------------------------------------------------*)
  // Primero identificamos la matriz A del filtro por mínimos cuadrados
  res := calcmincuad_mrh_mc(seriesDeDatos.get_series_cronActiva, ordenfiltro,
    filtrarMenores, umbralFiltro, seriesDeDatos.rNPPorCiclo,
    seriesDeDatos.rOffsetCiclo, overlapping, traslapping,
    NPuntosPorMiniciclo, mcA{%H-}, mcB{%H-}, True, seriesDeDatos.nombresSeries,
    flg_InfoParcial, flg_CompletarHuecosConRuidoBlanco, VarNoExplicado{%H-});


  SintetizadorSimSEE.mcA := mcA;
  SintetizadorSimSEE.mcB := mcB;

  assignfile(fsal, 'mab.xlt');
  rewrite(fsal);
  k := 0;
  for kfil := 1 to mcA[k].nf do
  begin
    nFilas:= mcA[k].nf;
    for jcol := 1 to mcA[k].nc do
    begin
      nRetardo:= (( jcol -1 ) div nFilas )+1;
      kSerie:= ( jcol - 1)  mod  nFilas;
      Write(fsal, #9, seriesDeDatos.nombresSeries[kFil-1] + '-' +
        seriesDeDatos.nombresSeries[ kSerie ]+'_'+IntToStr( nRetardo ) );
    end;
  end;
  Write(fsal, #9, '|');


  for kfil := 1 to mcB[k].nf do
    for jcol := 1 to mcB[k].nc do
    begin
      Write(fsal, #9, seriesDeDatos.nombresSeries[kfil - 1] + '-R_' +
        IntToStr( jCol ));
    end;
  writeln(fsal);

  for k := 0 to high(mcA) do
  begin

    writeln('k: ', k);
    Write(fsal, k);
    if mcA[k] <> nil then
    begin
      for kfil := 1 to mcA[k].nf do
        for jcol := 1 to mcA[k].nc do
          Write(fsal, #9, mcA[k].e(kfil, jcol));
      Write(fsal, #9, '|');
      if (mcB[k] <> nil) then
      begin
        for kfil := 1 to mcB[k].nf do
          for jcol := 1 to mcB[k].nc do
            Write(fsal, #9, mcB[k].e(kfil, jcol));
      end
      else
        Write(fsal, #9, 'MCB[' + IntToStr(k) + ']=NIL');
    end;
    writeln(fsal);
  end;
  closefile(fsal);

  Assignfile(fsal, 'matricesA.xlt');
  rewrite(fsal);
  for k := 0 to high(mcA) do
  begin
    writeln(fsal, 'paso: ', #9, k);
    if mcA[k] <> nil then
      mcA[k].WriteXLT(fsal)
    else
      writeln(fsal, 'SISTEMA NO INVERTIBLE');
  end;
  closefile(fsal);


  Assignfile(fsal, 'matricesB.xlt');
  rewrite(fsal);
  for k := 0 to high(mcB) do
  begin
    writeln(fsal, 'paso: ', #9, k);
    if mcB[k] <> nil then
      mcB[k].WriteXLT(fsal)
    else
      writeln(fsal, 'Falló CHOLESKY');
  end;
  closefile(fsal);



  // bien ahora vamos a completar la información de reducción de estados.
  for j := 1 to ordenfiltro do
    for k := 1 to seriesDeDatos.NSeries do
    begin

      if ordenfiltro > 1 then
        SintetizadorSimSEE.nombreVarE[(k - 1) + (j - 1) * SeriesDeDatos.NSeries] :=
          seriesDeDatos.nombresSeries[k - 1] + '_' + IntToStr(j)
      else
        SintetizadorSimSEE.nombreVarE[(k - 1) + (j - 1) * SeriesDeDatos.NSeries] :=
          seriesDeDatos.nombresSeries[k - 1];

      SintetizadorSimSEE.MRed_cte.pon_e(k, (k - 1) * ordenfiltro + j, 1);
      SintetizadorSimSEE.nDiscsVsE[(k - 1) * ordenfiltro + j - 1] := 5;
      setlength(SintetizadorSimSEE.ProbsVsE[(k - 1) * ordenfiltro + j - 1], 5);
      for h := 0 to 4 do
        SintetizadorSimSEE.ProbsVsE[(k - 1) * ordenfiltro + j - 1][h] := 0.2;
    end;


  SintetizadorSimSEE.WriteToArchi(archi,
    NDigitosDeformadores, NDecimalesDeformadores);
  Result := res;
end;


procedure TForm1.BBuscarClick(Sender: TObject);
begin
  if OpenDialog1.Execute then
    eArchiSeriesDeDatos.Text := OpenDialog1.FileName;
end;


function TForm1.buscarAnio(anio: integer): integer;
var
  buscando: boolean;
  kres: integer;
begin
  kres := -1;
  buscando := True;
  while buscando and (kBuscarDesde < length(aniosanalogos)) do
    if aniosanalogos[kBuscarDesde] = anio then
    begin
      buscando := False;
      kres := kBuscarDesde;
    end
    else
    if aniosanalogos[kBuscarDesde] < anio then
      Inc(kBuscarDesde)
    else
      buscando := False;

  Result := kres;

end;

function TForm1.Mascarado(anio, semana: integer): boolean;
var
  anioInicio: integer;
begin
  if semana >= self.semanaDeInicio then
    anioInicio := anio
  else
    anioInicio := anio - 1;
  Result := buscarAnio(anioInicio) < 0;
end;


procedure TForm1.MascararAniosNoAnalogos(datos: TSeriesDeDatos; mascara: double);
var
  kserie: integer;
  anio, isem, ksem: integer;
begin
  for kserie := 0 to high(datos.get_series_cronActiva) do
  begin
    anio := datos.PM_Anio;
    kBuscarDesde := 0;
    for isem := 1 to datos.series[0].n do
    begin
      // supongo que empiezan las series en la primer semana del año.
      ksem := ((isem - 1) mod 52) + 1;
      if mascarado(anio, ksem) then
        datos.series[kserie].pon_e(isem, mascara);
      if ksem = 52 then
        Inc(anio);
    end;
  end;
end;


function sino(x: boolean): string;
begin
  if x then
    Result := 'Si'
  else
    Result := 'No';
end;

procedure TForm1.BAnalizarClick(Sender: TObject);
var
  varNoExplicada: NReal;
  varNoExplicada_Vect: TVectR;
  archi: string;
  kSerie: integer;
  j: integer;
  ih0: NReal;
  dtMuestra: TDateTime;
  anular_negativos_hasta: integer;
  v_nCronicasRuido: integer;
  v_PrecisionMedida, v_UmbralFiltro, v_UmbralFiltroAlGaussianizar: NReal;
  v_OrdenDelFiltro, v_Traslapping: integer;
  v_NPuntosPorMiniCiclo, v_HorasEntreMuestras: integer;
  flg_FiltrarMenores, flg_FiltrarAlGaussianizar, flg_InfoParcial: Boolean;
  v_NDigitosDeformadores, v_NDecimalesDeformadores,
  v_NPuntosDeformador: integer;
  s: string;

begin
  varNoExplicada_Vect := nil;

  if FileExists(eArchiSeriesDeDatos.Text) then
  begin
    if seriesDeDatos <> nil then
      seriesDeDatos.Free;

    seriesDeDatos := TSeriesDeDatos.CreateFromArchi(eArchiSeriesDeDatos.Text);

    if seriesDeDatos.rNPPorCiclo < StrToInt(eNPuntosPorMiniciclo.Text) then
    begin
      ShowMessage(
        'NPuntosPorPeriodo < NPuntosPorMiniciclo; Impongo: NPuntosPorPeriodo = NPuntosPorMniciclo='
        + eNPuntosPorMiniciclo.Text);
      seriesDeDatos.rNPPorCiclo := StrToInt(eNPuntosPorMiniciclo.Text);
    end;

    if cbFiltradoKTPorIh0.Checked then
    begin
      (*** REFILTRO LOS KTs****)
      for j := 1 to seriesDeDatos.NPuntos do
      begin
        dtMuestra := seriesDeDatos.dtMuestra(j);
        Ih0 := RadiacionSolarExtraterrestrePlanoHorizontal_Ih0(
          dtMuestra, -31.430612, -57.90961, -3);
        //     seriesDeDatos.series[4].pon_e( j, Ih0 );
        if Ih0 < 0.2 then
          for kSerie := 0 to 3 do
            seriesDeDatos.series[kSerie].pon_e(j, -1234000);
      end;
      archi := eArchiSeriesDeDatos.Text;
      archi := ChangeFileExt(archi, '.fil');
      seriesDeDatos.WriteToArchi(archi);
    end;

    if self.cbAniosAnalogos.Checked then
    begin
      GetAniosAnalogos;
      MascararAniosNoAnalogos(seriesDeDatos, strToFloat(eUmbralFiltro.Text));
    end;

    if not cbFiltroCiclico.Checked then
    begin

      if cbFiltradoKTPorIh0.Checked then
        anular_negativos_hasta := 4
      else
        anular_negativos_hasta := 0;

// Pasamos los valores de los campos del formularios a variables antes de llamar
// a la función para facilitar el debugger de errores de conversión.
      v_nCronicasRuido:=  StrToInt( eNCronicasRuido.text );
      v_PrecisionMedida:= StrToFloat( ePrecisionMedida.Text );
      v_OrdenDelFiltro:= StrToInt(eOrdenDelFiltro.Text);
      v_Traslapping:= StrToInt(eTraslapping.Text);
      v_NPuntosPorMiniCiclo:= StrToInt( eNPuntosPorMiniciclo.Text );
      v_HorasEntreMuestras:= trunc( seriesDeDatos.dtEntreMuestras * 24 + 0.1) ;
      flg_FiltrarMenores:= cbFiltrarMenores.Checked;
      v_UmbralFiltro:= strToFloat(eUmbralFiltro.Text);
      flg_FiltrarAlGaussianizar:= cb_FiltrarAlGaussianizar.Checked;
      v_UmbralFiltroAlGaussianizar:= strToFloat(eUmbralAlGaussianizar.Text);
      v_NDigitosDeformadores:= StrToInt(eNDigitosDesformadores.Text);
      v_NDecimalesDeformadores:= StrToInt(eNDecimalesDesformadores.Text);
      v_NPuntosDeformador:= StrToInt(eNP_Deformadores.Text);
      flg_InfoParcial:= cb_InfoParcial.Checked;

      if rgTipoModelo.ItemIndex = 0 then
        varNoExplicada_Vect :=
          Analizar('SintetizadorCEGH.txt', seriesDeDatos,
          v_nCronicasRuido, v_PrecisionMedida,
          v_OrdenDelFiltro, StrToInt(eOverlapping.Text),
          v_Traslapping, v_NPuntosPorMiniCiclo,
          v_HorasEntreMuestras, flg_FiltrarMenores,
          v_UmbralFiltro, flg_FiltrarAlGaussianizar,
          v_UmbralFiltroAlGaussianizar, self.rbg_calculoB.ItemIndex,
          v_NDigitosDeformadores, v_NDecimalesDeformadores,
          v_NPuntosDeformador, flg_InfoParcial,
          anular_negativos_hasta )
      else
        varNoExplicada_Vect :=
          Analizar_Y_AX_B('SintetizadorCEGH_Y_AX_B.txt', seriesDeDatos,
          StrToInt( eNCronicasRuido.text ), StrToFloat( ePrecisionMedida.Text ),
          StrToInt(eOrdenDelFiltro.Text), StrToInt(eOverlapping.Text),
          StrToInt(eTraslapping.Text), StrToInt(eNPuntosPorMiniciclo.Text),
          v_HorasEntreMuestras, cbFiltrarMenores.Checked,
          strToFloat(eUmbralFiltro.Text), cb_FiltrarAlGaussianizar.Checked,
          strToFloat(eUmbralAlGaussianizar.Text), self.rbg_calculoB.ItemIndex,
          StrToInt(eNDigitosDesformadores.Text), StrToInt(eNDecimalesDesformadores.Text),
          StrToInt(eNP_Deformadores.Text), cb_InfoParcial.Checked,
          anular_negativos_hasta);
    end
    else
    begin
      if rgTipoModelo.ItemIndex = 0 then
        varNoExplicada_vect := Analizar_FiltroVariable('SintetizadorCEGH_mc.txt',
          seriesDeDatos,
          StrToInt( eNCronicasRuido.text ), StrToFloat( ePrecisionMedida.Text ),
          StrToInt(eOrdenDelFiltro.Text),
          StrToInt(eOverlapping.Text), StrToInt(eTraslapping.Text),
          StrToInt(eNPuntosPorMiniciclo.Text), v_HorasEntreMuestras,
          cbFiltrarMenores.Checked, strToFloat(eUmbralFiltro.Text),
          cb_FiltrarAlGaussianizar.Checked, strToFloat(eUmbralAlGaussianizar.Text),
          StrToInt(eNDigitosDesformadores.Text),
          StrToInt(eNDecimalesDesformadores.Text), StrToInt(eNP_Deformadores.Text),
          cb_InfoParcial.Checked, cbCompletarHuecosConRuidoBlanco.Checked )
      else
        ShowMessage('El método Multi Filtro no está implementado para modelado Y_AX_B');
    end;


    if varNoExplicada_Vect = nil then
      eVarianzaDeLoNoExplicado.Text := FloatToStr(varNoExplicada)
    else
    begin
      eVarianzaDeLoNoExplicado.Text := varNoExplicada_Vect.serialize;
      varNoExplicada_Vect.Free;
      varNoExplicada_Vect := nil;
    end;

    if seriesDeDatos <> nil then
    begin
      seriesDeDatos.Free;
      seriesDeDatos := nil;
    end;
  end
  else
    ShowMessage('Nombre de Archivo Inválido: ' + eArchiSeriesDeDatos.Text);
end;

procedure TForm1.btCalcularErroresClick(Sender: TObject);
var
  archi: string;
  SeriesDeErrores: TSeriesDeDatos;

begin
  OpenDialog1.Title := 'Seleccione el CEGH';
  if not OpenDialog1.Execute then
    exit;
  archi := OpenDialog1.FileName;
  if cegh <> nil then
    cegh.Free;
  cegh := TModeloCEGH.CreateFromArchi(archi);

  OpenDialog1.Title := 'Seleccione las series';
  if not OpenDialog1.Execute then
    exit;
  archi := OpenDialog1.FileName;

  if seriesDeDatos <> nil then
    seriesDeDatos.Free;
  seriesDeDatos := TSeriesDeDatos.CreateFromArchi(archi);

  SeriesDeErrores := cegh.Calc_SeriesError(seriesDeDatos,
    cb_Filtrarmenores_CalcErrores.Checked, StrToFloat(eUmral_calcerrores.Text));

  SeriesDeErrores.WriteToArchi('SeriesDeErrores.sas');
end;

procedure TForm1.btCargarCEGHClick(Sender: TObject);
var
  archi: string;
begin
  if not OpenDialog1.Execute then
    exit;
  archi := OpenDialog1.FileName;
  cegh := TModeloCEGH.CreateFromArchi(archi);

end;

procedure TForm1.btCompletarSeriesClick(Sender: TObject);
var
  rellenador: TRelleandorDeHuecos_CEGH;
  series_originales: TSeriesDeDatos;//serie a rellenar
  k, j: integer;
begin
  if FileExists(eArchiSeriesDeDatos.Text) then
    //se cargó el archivo de series de datos
  begin
    if seriesDeDatos <> nil then
      //si ya habia una serie de datos cargada la borro
      seriesDeDatos.Free;

    //cargo la serie de datos desde el archivo
    seriesDeDatos := TSeriesDeDatos.CreateFromArchi(eArchiSeriesDeDatos.Text);
    rellenador := TRelleandorDeHuecos_CEGH.Create(cegh, strToFloat(eUmbralFiltro.Text));

    //dbg: guardo la serie de datos original en este archivo para comparar luego de rellenar
    series_originales := TSeriesDeDatos.CreateFromArchi(eArchiSeriesDeDatos.Text);

    //pasa las series al mundo gaussiano, a menos que el valor sea menor al umbral pasado
    cegh.xTog_Series(seriesDeDatos, True, -11000);
    //dbg: guardo la serie de datos gaussiana
    seriesDeDatos.WriteToArchi('c:\basura\series_gaussianas_originales.txt');
    //Se rellenan los huecos en las series de a tramos
    rellenador.Rellenar(seriesDeDatos.get_series_cronActiva, 1800, seriesDeDatos);
    //dbg: guardo la serie gaussiana ya rellenada
    seriesDeDatos.WriteToArchi('c:\basura\series_gaussianas_rellenadas.txt');
    //vuelvo al espacio real
    cegh.gTox_Series(seriesDeDatos, True, -11000);
    for k := 1 to SeriesDeDatos.NPuntos do
      for j := 1 to SeriesDeDatos.NSeries do
        if series_originales.series[j - 1].e(k) > -10000 then
          seriesDeDatos.series[j - 1].pon_e(k, series_originales.series[j - 1].e(k));
    seriesDeDatos.WriteToArchi('c:\basura\series_completadas.txt');
  end;
end;

procedure TForm1.btCrearSASClick(Sender: TObject);
var
  aSeries: TSeriesDeDatos;

begin
  aSeries := TSeriesDeDatos.CreateVacia(now, 1, 100);
  aSeries.AddSerie('S1');
  aSeries.AddSerie('S2');
  aSeries.WriteToArchi('series_ejemplo.sas');
  aSeries.Free;
end;

procedure TForm1.btGenerarSeriesClick(Sender: TObject);
var
  kRealizacion: integer;
  nCiclos: integer;
  kUltimoCiclo: integer;
  DimX: integer;
  kDato, kPaso: integer;
  w, v: TVectR;
  semilla, kSerie: integer;
  kCronica: integer;
  nCronicas: integer;
  archi: string;
  cronicas: TList;
  aSeries: TSeriesDeDatos;
  estadistico: TSeriesDeDatos;
  mv: TVectR;
  NSeries: integer;

begin
  if FileExists(eArchiSeriesDeDatos.Text) then
  begin
    if seriesDeDatos = nil then
      seriesDeDatos := TSeriesDeDatos.CreateFromArchi(eArchiSeriesDeDatos.Text);

    NSeries := seriesDeDatos.NSeries;
    nCiclos := trunc(seriesDeDatos.NPuntos / seriesDeDatos.rNPPorCiclo);
    kUltimoCiclo := round((nCiclos - 1) * seriesDeDatos.rNPPorCiclo);
    cronicas := TList.Create;

    dimX := cegh.Dim_X;

    NCronicas := StrToInt(eNCronicas.Text);
    for kRealizacion := 1 to NCronicas do
    begin
      aSeries := seriesDeDatos.Clone;
      semilla := 30 + kRealizacion;
      cegh.GenRealizacion(aSeries, semilla);
      archi := 'c:\basura\series_sinteticas_' + IntToStr(semilla) + '.xlt';
      writeln('Escribiendo: ' + archi);
      aSeries.WriteToArchi(archi);
      if cbCopiarUltimomCicloComoEstadoinicial.Checked then
      begin
        // ahora copio el último ciclo como primero para que se altere en elstado inicial
        for kSerie := 1 to seriesDeDatos.NSeries do
        begin
          w := seriesDeDatos.series[kSerie - 1];
          v := aSeries.series[kSerie - 1];
          for kDato := 1 to DimX do
            w.pon_e(kDato, v.e(kUltimoCiclo + kDato));
        end;
      end;
      cronicas.add(aSeries);
    end;

    // Ahora procesamos para cada serie los conos de probabilidad 0.1, 0.9 y valor esperado.
    estadistico := seriesDeDatos.Clone;
    while estadistico.NSeries > 0 do
      estadistico.DelSerie(estadistico.NSeries - 1);

    for kSerie := 1 to aSeries.NSeries do
    begin
      estadistico.AddSerie(aSeries.nombresSeries[kSerie - 1] + '_pe10');
      estadistico.AddSerie(aSeries.nombresSeries[kSerie - 1] + '_pe90');
      estadistico.AddSerie(aSeries.nombresSeries[kSerie - 1] + '_ve');
    end;

    mv := TVectR.Create_init(nCronicas);

    for kPaso := 1 to estadistico.NPuntos do
    begin
      for kSerie := 0 to NSeries - 1 do
      begin
        for kCronica := 0 to NCronicas - 1 do
        begin
          aSeries := cronicas[kCronica];
          mv.pon_e(kCronica + 1, aSeries.series[kSerie].e(kPaso));
        end;
        mv.Sort(True);
        estadistico.series[kSerie * 3 + 0].pon_e(kPaso, mv.pe_VaR(0.1));
        estadistico.series[kSerie * 3 + 1].pon_e(kPaso, mv.pe_VaR(0.9));
        estadistico.series[kSerie * 3 + 2].pon_e(kPaso, mv.promedio);
      end;
    end;

    mv.Free;

    archi := 'c:\basura\estadistico_series.xlt';
    writeln('escribiendo: ', archi);
    estadistico.WriteToArchi(archi);

    for kCronica := 0 to cronicas.Count - 1 do
    begin
      aSeries := cronicas[kCronica];
      aSeries.Free;
    end;
    cronicas.Free;
  end;
end;

procedure TForm1.btModeladoAportesClick(Sender: TObject);
var
  archi: string;
begin
  if OpenDialog1.Execute then
  begin
    archi := OpenDialog1.FileName;
    if SeriesDeDatos <> nil then
      SeriesDeDatos.Free;
    memo_model.Lines.add('Cargando ...' + archi);
    SeriesDeDatos := TSeriesDeDatos.CreateFromArchi(archi);


  end;
end;




procedure TForm1.btSelSeriesClick(Sender: TObject);
var
  archi: string;
begin
  if OpenDialog1.Execute then
  begin
    archi := OpenDialog1.FileName;
    if SeriesDeDatos <> nil then
      SeriesDeDatos.Free;
    memo_model.Lines.add('Cargando ...' + archi);
    SeriesDeDatos := TSeriesDeDatos.CreateFromArchi(archi);
  end;
end;




procedure TForm1.btWeibulizarClick(Sender: TObject);
var
  cegh: TModeloCEGH;
  archi: string;
  vdatos: TVectR;
  weibull: Tf_ddp_Weibull;
  sal: textfile;
  kSerie, kPuntoT: integer;

begin
  if not OpenDialog1.Execute then
    exit;

  archi := OpenDialog1.FileName;

  assignfile(sal, 'weibullizator.xlt');
  rewrite(sal);
  writeln(sal, 'CEGH:', #9, archi);
  cegh := TModeloCEGH.CreateFromArchi(archi);

  for kSerie := 1 to cegh.nBornesSalida do
  begin
    writeln(sal);
    writeln(sal, 'serie: ', cegh.NombresDeBornes_Publicados[kSerie - 1]);
    writeln(sal, 'j', #9, 'k', #9, 'lambda');
    for kPuntoT := 1 to cegh.nPuntosPorPeriodo do
    begin
      vdatos := cegh.funcionesDeformantes[kSerie - 1][kPuntoT - 1].a;
      weibull := Tf_ddp_Weibull.Create_vmuestras(vdatos, nil, 0);
      writeln(sal, kPuntoT, #9, weibull.k, #9, weibull.lambda);
      weibull.Free;
    end;
  end;
  closefile(sal);
  weibull.Free;
end;



procedure TForm1.BtPotVelAeroGenClick(Sender: TObject);
var
  CurvaAero: TCurvaAerogen;
  kSerie, kPaso: integer;
  vserie: TVectR;
  P_MW, v: NReal;

begin
  if cbCurvaAero.ItemIndex < 0 then
  begin
    ShowMessage('Debe selección primero el modelo de aerogenerador.');
    exit;
  end;


  CurvaAero := curvasAeros_lst.Items[cbCurvaAero.ItemIndex];

  if FileExists(eArchiSeriesDeDatos.Text) then
  begin
    if seriesDeDatos <> nil then
      seriesDeDatos.Free;

    seriesDeDatos := TSeriesDeDatos.CreateFromArchi(eArchiSeriesDeDatos.Text);


    for kSerie := 0 to seriesDeDatos.NSeries - 1 do
    begin
      vserie := seriesDeDatos.series[kSerie];
      for kPaso := 1 to seriesDeDatos.NPuntos do
      begin
        P_MW := vserie.e(kPaso) / 1000.0;
        v := CurvaAero.v_of_P_MW(P_MW);
        vserie.pon_e(kPaso, v);
      end;
    end;

    seriesDeDatos.WriteToArchi('c:\basura\series_PotToVel.xlt');

  end;

  CurvaAero.Free;
end;

procedure TForm1.btClusteringClick(Sender: TObject);
var
  aSeries: TSeriesDeDatos;
  aSerie: TVectR;
  kSerie: integer;
  j: integer;
  PuntosMiniCiclo: integer;
  a: NReal;
  umbralFiltro: NReal;
  ma, mb: TDAOfNreal;
  flg_Exponencial: boolean;
  fsal: TextFile;
  d_Tendencia: NReal;

begin

  flg_Exponencial := cbTendenciaExponencial.Checked;
  PuntosMiniCiclo := StrToInt(eMinicicloCaracteristicoTendencia.Text);

  if seriesDeDatos = nil then
    btSelSeriesClick(Sender);

  if seriesDeDatos = nil then
    exit;
  aSeries := seriesDeDatos.Clone;

  // Hacemos los promedios móviles filtrados
  for kSerie := 0 to aSeries.NSeries do
    aSeries.series[kSerie].PromedioMovil_Filtrado(PuntosMiniCiclo, -1111);

  // Ahora si corresponde tomas los LN
  if flg_Exponencial then
  begin
    umbralFiltro := -1E30;
    for kSerie := 0 to aSeries.NSeries do
    begin
      aSerie := aSeries.series[kSerie];
      for j := 1 to aSerie.N do
      begin
        a := aSerie.pv[j];
        if a > AsumaCero then
          aSerie.pv[j] := LN(a)
        else
          aSerie.pv[j] := umbralFiltro - 1212;
      end;
    end;
  end;

  setlength(ma, aSeries.NSeries);
  setlength(mb, aSeries.NSeries);
  assignfile(fsal, 'tedencias.xlt');
  if cbTendenciaExponencial.Checked then
    writeln(fsal, 'Tipo: Exp( ma * k + mb )')
  else
    writeln(fsal, 'Tipo: ma * k + mb');
  writeln(fsal, 'NPromedioMiniciclo: ' + IntToStr(PuntosMiniCiclo));
  writeln(fsal, 'iSerie'#9'ma'#9'mb');
  for kSerie := 0 to aSeries.NSeries do
  begin
    aSerie := aSeries.series[kSerie];
    // Ahora calculamos la tendencia lineal
    aSerie.AproximacionLinealFiltrada(ma[kSerie], mb[kSerie], umbralFiltro);
    writeln(fsal, kSerie + 1, #9, ma[kSerie], #9, mb[kSerie]);
  end;
  closefile(fsal);


  // ahora quitamos las tendencias sobre las series originales
  for kSerie := 0 to aSeries.NSeries do
  begin
    aSerie := seriesDeDatos.series[kSerie];
    for j := 1 to aSerie.n do
    begin
      a := aSerie.pv[j];
      if a > -1111 then
      begin
        if flg_Exponencial then
        begin
          d_tendencia := exp(ma[kSerie] * j + mb[kSerie]);
          aSerie.pv[j] := a / d_tendencia;
        end
        else
        begin
          d_tendencia := ma[kSerie] * j + mb[kSerie];
          aSerie.pv[j] := a - d_tendencia;
        end;
      end;
    end;
  end;

  // Ahora salvamos las series destendenciadas con otro nombre
  seriesDeDatos.WriteToArchi('series_sintendencia.sas');

end;

procedure TForm1.Button1Click(Sender: TObject);
begin

end;


procedure CalcFourier( f: TLibroOpenCalc; M, ucos, usin: TMatR; celda: string );
var
  kPaso: integer;
  kArmonico: integer;
  ve: TVectR;
  ca: TMatR;
  cb: TMatR;
  NArmonicos: integer;
  fil, col: integer;
  af: TVectR;
begin
    NArmonicos:= ucos.nf;

    ve:= TVectR.Create_Init( M.nf  );
    ca:= TMatR.Create_Init( M.nf, NArmonicos );
    cb:= TMatR.Create_Init( M.nf, NArmonicos );

    for kPaso:= 1 to M.nf do
    begin
      af:= M.fila( kPaso );
      ve.pon_e( kPaso, af.promedio );
      for kArmonico := 1 to NArmonicos do
      begin
        ca.pon_e( kPaso, kArmonico, ucos.Fila(kArmonico ).PEV( af ) );
        cb.pon_e( kPaso, kArmonico, usin.Fila(kArmonico ).PEV( af ) );
      end;
    end;


    f.CeldaStrToFilCol( fil, col, celda );
    f.ir( fil, col );
    for kPaso:= 1 to M.nf do
    begin
      f.Write( ve.e( kPaso ) );
      for kArmonico := 1 to NArmonicos do
      begin
        f.Write( ca.e( kPaso, kArmonico ));
        f.Write( cb.e( kPaso, kArmonico ));
      end;
      inc( fil );
      f.ir( fil, col );
    end;

    ve.Free;
    ca.Free;
    cb.Free;
end;


procedure TForm1.Button2Click(Sender: TObject);
var
  f: TLibroOpenCalc;
  M_Feriado, M_SemiFeriado, M_Habil: TMatR;

  ve: TVectR;
  ca, cb: TMatR;

  kPaso: integer;

  ucos, usin: TMatR;
  NPuntos: integer;
  kArmonico: integer;
  NArmonicos: integer;
begin

  f:= TLibroOpenCalc.Create( false, '/home/rchaer/SimSEE/datos_comunes/CurvasDiaras_pu_TDD_Temp.ods' );
  M_Feriado:= f.CreateRead_MatR( 'Hoja1', 'C3:Z54' );
  M_SemiFeriado:= f.CreateRead_MatR( 'Hoja1', 'C60:Z111' );
  M_Habil:= f.CreateRead_MatR( 'Hoja1', 'C119:Z380' );

  NPuntos:= M_Feriado.nc;
  NArmonicos:= 5;
  ucos:= TMatR.Create_Init_pm( NArmonicos, NPuntos );
  usin:= TMatR.Create_init_pm( NArmonicos, NPuntos );

  for kArmonico:= 1 to NArmonicos do
  begin
    ucos.pm[kArmonico]:= create_cos_n( kArmonico, NPuntos );
    usin.pm[kArmonico]:= create_sin_n( kArmonico, NPuntos );
  end;

  CalcFourier( f,  M_Feriado, ucos, usin, 'AB3' );
  CalcFourier( f,  M_SemiFeriado, ucos, usin, 'AB60' );
  CalcFourier( f,  M_Habil, ucos, usin, 'AB119' );

  f.WriteToFile( 'salida.ods' );

  f.Free;

  ucos.Free;
  usin.Free;

  M_Feriado.Free;
  M_SemiFeriado.Free;
  M_Habil.Free;
end;

procedure TForm1.FormCreate(Sender: TObject);
var
  k: integer;
  aCurvaAero: TCurvaAerogen;
begin

  curvasAeros_lst := ucurvasaeros.GetListaDeCurvasConocidas;
  cbCurvaAero.Items.Clear;
  for k := 0 to curvasAeros_lst.Count - 1 do
  begin
    aCurvaAero := curvasAeros_lst.items[k];
    cbCurvaAero.Items.Add(aCurvaAero.modelo);
  end;

  ePrecisionMedida.Text:= FloatToStr( 0.03 );

end;

procedure TForm1.memo_modelChange(Sender: TObject);
begin

end;

procedure TForm1.P_TransicionClick(Sender: TObject);
var
  Aux_CEGH: TTransicion_CEGH;
  xs_pdf: TVectR;
begin
  if FileExists(eArchiSeriesDeDatos.Text) then
  begin
    if seriesDeDatos <> nil then
      seriesDeDatos.Free;
    seriesDeDatos := TSeriesDeDatos.CreateFromArchi(eArchiSeriesDeDatos.Text);
    cegh.xTog_Series(seriesDeDatos, True, -11000);
    Aux_CEGH := TTransicion_CEGH.Create(cegh);
    xs_pdf := Aux_CEGH.Calc_PDF_Serie(seriesDeDatos);
    seriesDeDatos.AddSerie('dens_prob', xs_pdf, SALIDA);
    seriesDeDatos.WriteToArchi('c:\basura\series_completadas_con_pdf.txt');
  end;

end;


procedure TForm1.ESHYClick(Sender: TObject);
const
  ceB = 0.62739306;
  ceP = 0.24865988;
  ceS = 0.12394705;
var
  varNoExplicada: NReal;
  varNoExplicada_Vect: TVectR;
  archi: string;
  vectempH: TVectR;
  i, Totalseries, j, k, kpunto, largo: integer;
begin
  if FileExists(eArchiSeriesDeDatos.Text) then
  begin
    if seriesDeDatos <> nil then
      seriesDeDatos.Free;

    for i := 1 to StrToInt(eOrdenESHY.Text) do
    begin

      seriesDeDatos := TSeriesDeDatos.CreateFromArchi(eArchiSeriesDeDatos.Text);
      Totalseries := seriesDeDatos.NSeries;
      seriesDeDatos.NSeries := Totalseries + 1;

      seriesDeDatos.AddSerie('TotalSeries', nil, ENTRADA );
      vectempH := TVectR.Create_Init(seriesDeDatos.NPuntos);
      for kpunto := 1 to seriesDeDatos.NPuntos do
      begin
        vectempH.pv[kpunto] :=
          ceB * seriesDeDatos.series[0].pv[kpunto] + ceP *
          seriesDeDatos.series[1].pv[kpunto] + ceS *
          seriesDeDatos.series[2].pv[kpunto] * 2;
      end;

      seriesDeDatos.nombresSeries.Add('H_' + IntToStr(i));

      for kpunto := 1 to seriesDeDatos.NPuntos do
      begin
        for j := 0 to i - 1 do
        begin
          if (kpunto - j) >= 0 then
          begin
            seriesDeDatos.series[Totalseries].pv[kpunto] :=
              seriesDeDatos.series[Totalseries].pv[kpunto] +
              vectempH.pv[kpunto - j];
          end;
        end;
        if kpunto < i then
          k := kpunto
        else
          k := i;
        seriesDeDatos.series[Totalseries].pv[kpunto] :=
          seriesDeDatos.series[Totalseries].pv[kpunto] / (k);
      end;
      seriesDeDatos.WriteToArchi(eArchiSeriesDeDatos.Text +
        'serie_LikeESHY_H_' + IntToStr(i) + '.txt');
      //seriesDeDatos.n
      if self.cbAniosAnalogos.Checked then
      begin
        GetAniosAnalogos;
        MascararAniosNoAnalogos(seriesDeDatos, strToFloat(eUmbralFiltro.Text));
      end;
      archi := eArchiSeriesDeDatos.Text;
      largo := Length(archi) - 4;
      //archi:=Delete(archi,largo,4);
      archi := 'CEGH_' + archi + '_LikeESHY_H_' + IntToStr(i) + '.txt';
      if not cbFiltroCiclico.Checked then
      begin
        varNoExplicada_Vect :=
          Analizar(archi, seriesDeDatos,
          StrToInt( eNCronicasRuido.text ), StrToFloat( ePrecisionMedida.Text ),
          StrToInt(eOrdenDelFiltro.Text),
          StrToInt(eOverlapping.Text), StrToInt(eTraslapping.Text),
          StrToInt(eNPuntosPorMiniciclo.Text),
          round( seriesDeDatos.dtEntreMuestras*24 ),
          cbFiltrarMenores.Checked, strToFloat(eUmbralFiltro.Text),
          cb_FiltrarAlGaussianizar.Checked, strToFloat(eUmbralAlGaussianizar.Text),
          self.rbg_calculoB.ItemIndex, StrToInt(eNDigitosDesformadores.Text),
          StrToInt(eNDecimalesDesformadores.Text), StrToInt(eNP_Deformadores.Text),
          cb_InfoParcial.Checked);
        varNoExplicada_Vect.Free;
      end
      else
      begin
        varNoExplicada_Vect :=
          Analizar_FiltroVariable(archi, seriesDeDatos,
          StrToInt( eNCronicasRuido.text ), StrToFloat( ePrecisionMedida.Text ),
          StrToInt(eOrdenDelFiltro.Text), StrToInt(eOverlapping.Text),
          StrToInt(eTraslapping.Text), StrToInt(eNPuntosPorMiniciclo.Text),
          round( seriesDeDatos.dtEntreMuestras*24 ), cbFiltrarMenores.Checked,
          strToFloat(eUmbralFiltro.Text), cb_FiltrarAlGaussianizar.Checked,
          strToFloat(eUmbralAlGaussianizar.Text), StrToInt(
          eNDigitosDesformadores.Text), StrToInt(eNDecimalesDesformadores.Text),
          StrToInt(eNP_Deformadores.Text), cb_InfoParcial.Checked, cbCompletarHuecosConRuidoBlanco.Checked);

      end;
    end;

    if seriesDeDatos <> nil then
    begin
      seriesDeDatos.Free;
      seriesDeDatos := nil;
    end;
  end
  else
    ShowMessage('Nombre de Archivo Inválido: ' + eArchiSeriesDeDatos.Text);

end;

procedure TForm1.rgTipoModeloClick(Sender: TObject);
begin

end;

procedure TForm1.BLeerDatosBasicosClick(Sender: TObject);
begin
  if seriesDeDatos <> nil then
    seriesDeDatos.Free;
  seriesDeDatos := TSeriesDeDatos.CreateFromArchi(eArchiSeriesDeDatos.Text);
end;



procedure TForm1.BPredictorNino34Click(Sender: TObject);
var
  predictor: TModeloCEGH;
  //Es de dimensión NPasos X NCronicas.
  //Guarda para cada paso hacia adelante los resultados de los sorteos de cada crónica
  pronostico: array of TVectR;
  NCronicas: integer;
  NPasos: integer;
  kPaso: integer;
  iCronica: integer;
  serieHistorica: TSeriesDeDatos;
  y, x, r: TVectR;
  //Son de dimensión NPasos X (NDatos + NPasos)
  //Guardan para cada paso hacia adelante los valores medios, valores con
  //probabilidad 20% y valores con probabilidad 80% de los resultados obtenidos
  //en las crónicas de pronóstico
  vm, p20, p80: TMatR;
  f: textFile;
  iDesp: integer;
  NDatos: integer;
  NOrdenFiltro: integer;
  vy: NREal;
  i, j: integer;
begin
  serieHistorica := TSeriesDeDatos.CreateFromArchi(
    '\simsee\src-delphi\ax\analisisserial\serieNino34.txt');

  predictor := TModeloCEGH.CreateFromArchi(
    '\simsee\src-delphi\ax\analisisserial\SintCEGH-mensual-nino34-fo1.txt');

  x := predictor.CrearVector_EstadoX;
  y := predictor.CrearVector_Salida;
  r := predictor.CrearVector_RBG;

  NOrdenFiltro := predictor.calcOrdenDelFiltro;

  NDatos := serieHistorica.NPuntos;

  NCronicas := 1000;

  NPasos := 5;

  vm := TMatR.Create_Init(NPasos, NDatos + NPasos);
  p20 := TMatR.Create_Init(NPasos, NDatos + NPasos);
  p80 := TMatR.Create_Init(NPasos, NDatos + NPasos);

  setlength(pronostico, NPasos);
  for kPaso := 0 to high(pronostico) do
    pronostico[kPaso] := TVectR.Create_Init(NCronicas);

  for idesp := 1 to NDatos do
  begin
    for iCronica := 1 to NCronicas do
    begin
      // cargamos el vector de estados desde la serie de datos
      predictor.InicializarDesdeDatosReales(x, serieHistorica, idesp, 0);
      // ?? gaussianar

      for kPaso := 0 to NPasos - 1 do
      begin
        // sortear salidas de ruido blanco gaussiano
        predictor.SortearValores(r);

        // calcula y= Ax + Br
        predictor.CalcularProximasSalidas(y, x, r, 0);


      {$IFDEF DBG_ESTIMG}
        vy := y.pv[1];
      {$ELSE}
        // desgaussianar
        vy := predictor.gTox(y.pv[1], 1, serieHistorica.kDefomador(idesp - 1 + kPaso));
      {$ENDIF}


        // cambia el estado
        predictor.EvolucionarEstado_(x, y);

        pronostico[kPaso].pon_e(iCronica, vy);
      end;
    end;

    for kPaso := 0 to NPasos - 1 do
    begin
      pronostico[kPaso].Sort(True);
      vm.pon_e(1 + kpaso, idesp + kpaso + 1, pronostico[kPaso].promedio);
      p20.pon_e(1 + kpaso, idesp + kpaso + 1, pronostico[kPaso].e(
        trunc(0.2 * NCronicas + 0.5)));
      p80.pon_e(1 + kpaso, idesp + kpaso + 1, pronostico[kPaso].e(
        trunc(0.8 * NCronicas + 0.5)));
    end;
  end;

  assignFile(f, 'PronosticoNino34-0712.xlt');
  rewrite(f);

  Write(f, #9, 'vr');
  for i := 1 to vm.nf do
    Write(f, #9, 'vm' + IntToStr(i));
  for i := 1 to vm.nf do
    Write(f, #9, 'p20' + IntToStr(i));
  for i := 1 to vm.nf do
    Write(f, #9, 'p80' + IntToStr(i));
  writeln(f);

  for j := 1 to vm.nc do
  begin
    Write(f, j);
    if j < NDatos then
    begin
      {$IFDEF DBG_ESTIMG}
      vy := predictor.xTog(serieHistorica.series[0].e(j), 1,
        ((j - 1) mod serieHistorica.NPPorCiclo) + 1);
      Write(f, #9, vy: 12: 4);
      {$ELSE}
      Write(f, #9, serieHistorica.series[0].e(j): 12: 4);
      {$ENDIF}
    end
    else
      Write(f, #9, '');

    for i := 1 to vm.nf do
      Write(f, #9, vm.e(i, j): 12: 4);

    for i := 1 to vm.nf do
      Write(f, #9, p20.e(i, j): 12: 4);

    for i := 1 to vm.nf do
      Write(f, #9, p80.e(i, j): 12: 4);

    writeln(f);
  end;
  closeFile(f);
end;

procedure TForm1.BTestVarianzaClick(Sender: TObject);
var
  predictor: TModeloCEGH;
  NCronicas: integer;
  x, y, r: TVectR;
  icronica: integer;
  y2: TVectR;
  j: integer;
begin
  if not OpenDialog1.Execute then
    exit;

  predictor := TModeloCEGH.CreateFromArchi(OpenDialog1.FileName);

(*
  predictor:= TModeloCEGH.CreateFromArchi(
    'C:\simsee\datos_comunes\sintetizadores\SintCEGH-mensual-nino34-fo3.txt');
    *)
  x := predictor.CrearVector_EstadoX;
  y := predictor.CrearVector_Salida;
  y2 := TVectR.Create_Clone(y);
  r := predictor.CrearVector_RBG;
  NCronicas := 1000000;


  // le damos un poco de manija al filtro para que quede cargado
  for iCronica := 1 to 1000 do
  begin
    predictor.SortearValores(r);
    predictor.CalcularProximasSalidas(y, x, r, 0);
    predictor.EvolucionarEstado_(x, y);
  end;

  // ahora empezamos a calcular las varianzas de las salidas.
  for iCronica := 1 to NCronicas do
  begin
    predictor.SortearValores(r);
    predictor.CalcularProximasSalidas(y, x, r, 0);
    for j := 1 to y.n do
      y2.pv[j] := y2.pv[j] + sqr(y.pv[j]);
    predictor.EvolucionarEstado_(x, y);
  end;

  for j := 1 to y.n do
    y2.pv[j] := y2.pv[j] / NCronicas;

  writeln('Varianza de las series sintéticas.');
  for j := 1 to y.n do
    writeln('ys[', j, ']: ', y2.pv[j]);
end;



procedure TForm1.btAyuda_principlalClick(Sender: TObject);
begin
  verdoc( 'analisisserial-manualdeusuario', '');
end;

procedure TForm1.btHurstClick(Sender: TObject);
var
  hurst: array of THurst_rescaled;
  kserie: integer;
  sal: textfile;
  j: integer;

begin
  if FileExists(eArchiSeriesDeDatos.Text) then
  begin
    if seriesDeDatos <> nil then
      seriesDeDatos.Free;
    seriesDeDatos := TSeriesDeDatos.CreateFromArchi(eArchiSeriesDeDatos.Text);

    setlength(hurst, seriesDeDatos.NSeries);
    assignfile(sal, 'hurst.xlt');
    rewrite(sal);

    writeln(sal, 'serie', #9, 'slope', #9, 'intercept');
    for kserie := 0 to seriesDeDatos.NSeries - 1 do
    begin
      hurst[kserie] := THurst_Rescaled.Create(
        @(seriesDeDatos.series[kserie].pv[1]), seriesDeDatos.series[kserie].n);
      writeln(sal, kserie + 1, #9, hurst[kserie].slope, #9, hurst[kserie].intercept);
    end;

    writeln(sal);

    Write(sal, 'logN');
    for kserie := 0 to seriesDeDatos.NSeries - 1 do
      Write(sal, #9, 'logRS_' + seriesDeDatos.nombresSeries[kserie]);
    writeln(sal);

    for j := 1 to hurst[0].points_x.n do
    begin
      Write(sal, hurst[0].points_x.pv[j]);
      for kserie := 0 to seriesDeDatos.NSeries - 1 do
        Write(sal, #9, hurst[kserie].points_y.pv[j]);
      writeln(sal);
    end;

    for kserie := 0 to seriesDeDatos.NSeries - 1 do
    begin
      hurst[kserie].Free;
    end;

    closefile(sal);
  end;
end;

procedure TForm1.btAutocoClick(Sender: TObject);
var
  kserie: integer;
  sal: textfile;
  j: integer;

begin
  if FileExists(eArchiSeriesDeDatos.Text) then
  begin
    if seriesDeDatos <> nil then
      seriesDeDatos.Free;
    seriesDeDatos := TSeriesDeDatos.CreateFromArchi(eArchiSeriesDeDatos.Text);

    for kserie := 0 to seriesDeDatos.NSeries - 1 do
    begin
      //      seriesDeDatos.series[kserie].aplicar_hanning;
      AutoCorrelacion(seriesDeDatos.series[kserie].n,
        seriesDeDatos.series[kserie].pv[1]);
    end;
    assignfile(sal, 'autoco.xlt');
    rewrite(sal);

    Write(sal, 'k');
    for kserie := 0 to seriesDeDatos.NSeries - 1 do
      Write(sal, #9, 'r_' + seriesDeDatos.nombresSeries[kserie]);
    writeln(sal);

    for j := 1 to seriesDeDatos.series[0].n do
    begin
      Write(sal, j - 1);
      for kserie := 0 to seriesDeDatos.NSeries - 1 do
        Write(sal, #9, seriesDeDatos.series[kserie].pv[j]);
      writeln(sal);
    end;

    wrffti01.done;
    closefile(sal);
  end;
end;

procedure TForm1.btH2fClick(Sender: TObject);
var
  kserie: integer;
  sal: textfile;
  j: integer;
  ncoefs: integer;
  vm: NReal;
begin
  if FileExists(eArchiSeriesDeDatos.Text) then
  begin
    if seriesDeDatos <> nil then
      seriesDeDatos.Free;
    seriesDeDatos := TSeriesDeDatos.CreateFromArchi(eArchiSeriesDeDatos.Text);

    for kserie := 0 to seriesDeDatos.NSeries - 1 do
    begin
      vm := seriesDeDatos.series[kserie].promedio;
      seriesDeDatos.series[kserie].MasReal(-vm);
      seriesDeDatos.series[kserie].aplicar_hanning;
      h2f(seriesDeDatos.series[kserie].n, seriesDeDatos.series[kserie].pv[1]);
    end;

    assignfile(sal, 'h2f.xlt');
    rewrite(sal);

    Write(sal, 'k');
    for kserie := 0 to seriesDeDatos.NSeries - 1 do
      Write(sal, #9, 'h2_' + seriesDeDatos.nombresSeries[kserie]);
    writeln(sal);

    nCoefs := (seriesDeDatos.series[0].n div 2) + 1;
    for j := 1 to nCoefs do
    begin
      Write(sal, j - 1);
      for kserie := 0 to seriesDeDatos.NSeries - 1 do
        Write(sal, #9, seriesDeDatos.series[kserie].pv[j]);
      writeln(sal);
    end;

    wrffti01.done;
    closefile(sal);
  end;
end;




procedure TForm1.GetAniosAnalogos;
var
  s: string;
  k: integer;
  anio: integer;
begin
  s := trim(self.lstAniosAnalogos.Text);
  if s = '' then
  begin
    setlength(aniosanalogos, 0);
    exit;
  end;
  setlength(aniosanalogos, 10000);

  k := 0;
  while s <> '' do
  begin
    anio := nextint(s);
    aniosAnalogos[k] := anio;
    Inc(k);
  end;
  setlength(aniosAnalogos, k);

  semanaDeInicio := StrToInt(eSemanaInicioVentanaAnaloga.Text);
  kBuscarDesde := 0;
end;

procedure TForm1.BCalcularPromediosMensualesClick(Sender: TObject);
var
  seriesDeDatos: TSeriesDeDatos;
  vmm: array of TVectR;
  cnts: array of TVectE;
  dt, dtbase: TDateTime;

  NDatos: integer;
  ipaso: integer;
  dtpaso: TDateTime;
  anio, mes, dia: word;
  v: NReal;
  filtrar: boolean;
  umbralF: NReal;
  flgPasa: boolean;
  fsal: textfile;
  iserie: integer;

begin
  if cbFiltrarMenores.Checked then
  begin
    filtrar := True;
    umbralF := StrToFloat(eUmbralFiltro.Text);
  end
  else
    filtrar := False;

  seriesDeDatos := TSeriesDeDatos.CreateFromArchi(eArchiSeriesDeDatos.Text);

  setlength(vmm, seriesDeDatos.NSeries);
  setlength(cnts, seriesDeDatos.NSeries);
  for iserie := 0 to seriesDeDatos.NSeries - 1 do
  begin
    vmm[iserie] := TVectR.Create_Init(12);
    cnts[iserie] := TVectE.Create_Init(12);
  end;

  NDatos := seriesDeDatos.series[0].n;
  dt := EncodeDate(seriesDeDatos.PM_Anio, seriesDeDatos.PM_Mes, seriesDeDatos.PM_Dia);
  dtbase := dt + (seriesDeDatos.PM_Hora + (seriesDeDatos.PM_Minuto +
    seriesDeDatos.PM_segundo / 60) / 60) / 24;
  dtpaso := seriesDeDatos.PeriodoDeMuestreo_horas / 24.0;

  for iPaso := 1 to NDatos do
  begin
    for iSerie := 0 to seriesDeDatos.NSeries - 1 do
    begin
      flgPasa := True;
      v := seriesDeDatos.series[iSerie].e(ipaso);
      if filtrar then
        flgPasa := v > umbralF;

      if flgPasa then
      begin
        dt := dtbase + (ipaso - 1) * dtpaso;
        DecodeDate(dt, anio, mes, dia);
        vmm[iSerie].acum_e(mes, v);
        cnts[iSerie].acum_e(mes, 1);
      end;
    end;
  end;

  assignfile(fsal, 'promedios_mensuales.xlt');
  rewrite(fsal);
  writeln(fsal, DateTimeToStr(now()));
  writeln(fsal, 'ANALISIS_SERIAL_PROMEDIOS_MENSUALES');
  writeln(fsal, 'Serie de datos: ' + eArchiSeriesDeDatos.Text);
  writeln(fsal, 'Valores esperados por mes:::::');
  for mes := 1 to 12 do
  begin
    Write(fsal, mes: 8);
    for iserie := 0 to SeriesDeDatos.NSeries - 1 do
    begin
      if cnts[iserie].pv[mes] > 0 then
      begin
        vmm[iserie].pv[mes] := vmm[iserie].pv[mes] / cnts[iserie].pv[mes];
        Write(#9, vmm[iserie].pv[mes]: 12: 3);
      end
      else
        Write(fsal, #9, '-SIN VALOR-');
    end;
    writeln(fsal);
  end;
  closefile(fsal);
end;



procedure TForm1.TestDistribucionGaussianaClick(Sender: TObject);
begin
  Test_DistribucionNormalCanonica;
end;

procedure TForm1.GenerarRuidoRectificadoClick(Sender: TObject);
(*

FuenteRuido( 0, 1, AmpLitudRuido );
R( 1, 2 , 1E3 ohm );
C( 0, 2, 0.5E-6 F );
R( 0, 2, 1E3 ohm );
dt= 1e-3 segs;
salida= tensión(0,2)
*)

const
  NPuntos = 1024 * 4;
var
  k: integer;
  uk: NReal;
  Vk, Vks: NReal;
  dt, r1, r2, C: NReal;
  a1, a2: NReal;
  sal: textfile;
begin
  dt := 1E-3;
  r1 := 1E3;
  r2 := 1E3;
  C := 20E-6;
  a1 := 1 - dt / r2 / C;
  a2 := dt / r1 / C;
  Vk := 0;

  assignfile(sal, 'series_ruidorectificado.xls');
  rewrite(sal);
  writeln(sal, 1, #9, 'NSeries');
  writeln(sal, NPuntos, #9, 'NPuntos');
  writeln(sal, 1, #9, 'Puntos por ciclo');
  writeln(sal, #9, 'Ruido');
  for k := 0 to NPuntos - 1 do
  begin
    Vks := a1 * Vk;
    uk := (random - 0.5) * 1000;
    if uk > Vk then
      Vks := Vks + a2 * (uk - Vk);
    writeln(sal, k, #9, Vks: 12: 1);
    Vk := Vks;
  end;
  closefile(sal);
end;

end.
