   unit uDatosCronicosMensuales;
interface
uses
  Classes, SysUtils, Math;

type

TTira = array of double;

TDatosCronicos = class
  nMesesInicialesSinDato: integer;
  priCron, ultCron, priAEco, ultAEco: integer;
  cronicaInicioPrimerTira: integer;

  NAEcos, NMeses: integer;
  NCronicasS: integer; // cantidad de crónicas simuladas <= NCronicasN
  NCronicasH: integer; // cantidad de crónicas históricas de aportes.

  cronicas: array of TTira;
  promedio: TTira;
  tendencia: TTira;

  // almacena una serie larga desde priCron al ultCron
  // de valores mensuales sin tendencia.
  serie_mensual: TTira;
  serie_semanal: TTira;

  constructor Create( priAEco_, ultAEco_, priCron_, ultCron_: integer;
        cronicaInicioPrimerTira_: integer );
  procedure Free;

//  procedure Calc_Diff( cmo_pat_anterior: TDatosCronicos );
  procedure Calc_Promedio;
  procedure Determinar_MesesInicialesSinDatos;
  procedure Calc_Tendencia( ovlp_neg, ovlp_pos: integer; desplazarCaja: boolean );
  procedure QuitarTendencia;
  procedure Calc_SerieMensual;
  procedure Calc_SerieSemanal;
end;


function Resampling( v1: TTira; n2: integer ): TTira;

procedure GenerarSeriesYTendencia( archi: string );

implementation




function nextpal( var s: string; sep: string ): string; overload;
var
  r: string;
  i: integer;
begin
  i:= pos( sep, s );
  if i = 0 then
  begin
    r:= s;
    s:= '';
  end
  else
  begin
    r:= copy( s, 1, i-1 );
    delete( s, 1, i + length( sep ) - 1 );
  end;
  result:= r;
end;

function nextpal( var s: string ): string;  overload;
begin
  result:= nextpal( s, #9 );
end;

function nextint( var s: string ): integer;
begin
  result:= strToInt( trim( nextpal( s )) );
end;

function Resampling( v1: TTira; n2: integer ): TTira;
var
  v2: TTira;
  n1: integer;
  jr: double;
  j0, j1: integer;
  alfa: double;
  k: integer;

begin
  n1:= length( v1 );
  setlength( v2, n2 );
  (**
  for k:= 0 to n2 - 1 do
  begin
    jr:= k / n2 * n1;
    j0:= trunc( jr );
    j1:= min( j0 +1, ( n1 - 1 ) );
    if ( j1 > j0 ) then
    begin
      alfa:= jr - j0;
      v2[k]:= v1[ j0] * ( 1- alfa ) + v1[j1] * alfa;
    end
    else
      v2[k]:= v1[ j0 ];
  end;
  **)
  for k:= 0 to n2 - 1 do
  begin
    jr:= k / n2 * n1;
    j0:= trunc( jr );
    v2[k]:= v1[ j0 ];
  end;

  result:= v2;
end;


(*
procedure TDatosCronicos.Calc_Diff( cmo_pat_anterior: TDatosCronicos );
var
  icron, imes: integer;
begin
  for imes:= 0 to NMeses - 1 do
    for icron:= 0 to NCronicasS - 1 do
      cronicas[icron][imes]:= cmo_pat_anterior.cronicas[icron][imes] - cronicas[icron][imes];
end;
  *)
procedure TDatosCronicos.Calc_Promedio;
var
  icron, imes: integer;
  m: double;
begin
  for imes:= 0 to NMeses - 1 do
  begin
    m:= 0;
    for icron:= 0 to NCronicasS - 1 do
      m:= m + cronicas[icron][imes];
    m:= m / NCronicasS;
    promedio[imes]:= m;
  end;
end;


procedure TDatosCronicos.Determinar_MesesInicialesSinDatos;
var
  imes:  integer;
  buscando: boolean;
begin
  buscando:= true;
  imes:= 0;
  while buscando and (imes < NMeses) do
  begin
    if abs( promedio[imes] ) < 1.0E-18 then
      inc( imes )
    else
      buscando:= false;
  end;
  nMesesInicialesSinDato:= imes;
end;



procedure TDatosCronicos.Calc_Tendencia( ovlp_neg, ovlp_pos : integer; desplazarCaja: boolean );
var
  icron, imes: integer;
  m: double;
  cnt: integer;
  iov: integer;
  jmes: integer;
  desp: integer;
begin
// primero rellenamos el inicio sin datos con -11111 solo para que quede claro que ahí no vale
  for imes:= 0 to nMesesInicialesSinDato - 1 do
    tendencia[imes]:= -11111;

// ahora calculamos la tendencia como el promedio en el conjunto de las crónicas
// para el mismo casillero de mes y considerando que cada muestra desparrama hacia
// sua adyasentes de acuerdo al overlapping pasado como parámetro.
  for imes:= nMesesInicialesSinDato to NMeses -1 do
  begin
    cnt:= 0;
    m:= 0;

    desp:= 0; // caja de overlapping centrada
    if desplazarCaja  then
    begin
      if (imes - ovlp_neg ) < nMesesInicialesSinDato then
        desp:=  nMesesInicialesSinDato - ( imes - ovlp_neg )
      else
        if (imes + ovlp_pos ) > (NMeses - 1 ) then
          desp:= (NMeses - 1 )  - (imes + ovlp_pos );
    end;

    for iov:= -ovlp_neg to ovlp_pos do
    begin
      jmes:= imes + iov + desp;
      if ( jmes >= nMesesInicialesSinDato ) and ( jmes <= (NMeses - 1 ) ) then
      begin
        inc( cnt );
        m:= m + promedio[ jmes ];
      end;
    end;
    m:= m / cnt;
    tendencia[imes]:= m;
  end;

end;


procedure TDatosCronicos.QuitarTendencia;
var
  icron, imes: integer;

begin
  for icron:= 0 to NCronicasS -1 do
  begin
    for imes:= 0 to nMesesInicialesSinDato - 1 do
      cronicas[icron][imes]:= -11111;
   for imes:= nMesesInicialesSinDato to NMeses - 1 do
   begin
     cronicas[icron][imes]:= cronicas[icron][imes] / tendencia[imes];
   end;
  end;
end;



procedure TDatosCronicos.Calc_serieMensual;
var
  cnts: array of integer;
  icron: integer;
  jcron: integer;
  imes, jmes: integer;
begin
  setlength( cnts, length( serie_mensual ) );
  for imes:= 0 to high( cnts ) do cnts[imes]:= 0;

  for icron:= 0 to nCronicasS - 1 do
  begin
   for imes:= nMesesInicialesSinDato to nMeses - 1 do
   begin
      jcron:= cronicaInicioPrimerTira + icron + imes div 12;
      while jcron < priCron do jcron := jcron + NCronicasH;
      while jcron > ultCron do jcron := jcron - NCronicasH;
      jmes:= ( jcron -  priCron ) * 12 + imes mod 12;
      inc( cnts[ jmes ] );
      serie_mensual[jmes]:= serie_mensual[jmes] + cronicas[icron][imes];
   end;
  end;

  for imes:= 0 to high( cnts ) do
  begin
   serie_mensual[imes]:= serie_mensual[imes] / cnts[imes];
  end;
  setlength( cnts, 0 );
end;

procedure TDatosCronicos.Calc_serieSemanal;
begin
  setlength( serie_semanal, 0 );
  serie_semanal:= Resampling( serie_mensual,  (length( serie_mensual ) div 12 ) *52 );
end;



constructor TDatosCronicos.Create(
    priAEco_, ultAEco_, priCron_, ultCron_: integer;
    cronicaInicioPrimerTira_: integer );
var
  icron: integer;
begin
  priAEco:= priAEco_;
  ultAEco:= ultAEco_;
  priCron:= priCron_;
  ultCron:= ultCron_;
  cronicaInicioPrimerTira:= cronicaInicioPrimerTira_;

  NAEcos:= ultAEco - priAEco + 1;
  NCronicasH:= ultCron - priCron + 1;
  NCronicasS:= NCronicasH - ( cronicaInicioPrimerTira - priCron );
  NMeses:= nAEcos * 12;
  setlength( cronicas, nCronicasS );
  for icron := 0 to nCronicasS - 1 do
    setlength( cronicas[icron], nAEcos *12 );
  setlength( promedio, nAEcos * 12  );
  setlength( tendencia, nAEcos * 12 );
  setlength( serie_mensual, nCronicasH *12 );
  setlength( serie_semanal, nCronicasH * 52 );
end;

procedure TDatosCronicos.Free;
var
  icron: integer;
begin
  for icron := 0 to high( cronicas ) do
    setlength( cronicas[icron], 0 );
  setlength( cronicas, 0 );
  setlength( promedio, 0 );
  setlength( tendencia, 0 );
  setlength( serie_mensual, 0 );
  setlength( serie_semanal, 0 );
end;

procedure GenerarSeriesYTendencia( archi: string );
var
  fent: textfile;
  r: string;
  priAEco, ultAEco: integer;
  priCron, ultCron: integer;
  kpat: integer;
  cronicaInicioTira: integer;
  cronicaInicioPrimerTira: integer;
  
  NPats, NAecos, NCrons: integer;
  cmo: array of TDatosCronicos;
  pal: string;
  icron, iaeco, ipat, imes, kmes: integer;
  m: double;

  fsal: textfile;
  isem: integer;

begin
  priCron:= 1931;
  ultCron:= 2009;
  cronicaInicioPrimerTira:= -1;

  NCrons:= ultCron - priCron +1;
  NPats:= 3;

    assignfile( fent, archi );
    reset( fent );
    readln( fent, r );
    pal:= nextpal( r );
    if pal <> 'PEN' then
    begin
      raise Exception.Create( 'Error, pal <> PEN en linea 1' );
    end;

    pal:= nextpal( r );
    delete( pal, 1,1 );
    priAEco:= StrToInt( pal );
    pal:= trim( nextpal( r, '-' ) );
    ultAEco:= 2000 + StrToInt( pal );
    NAEcos:= ultAEco - priAEco + 1;

    readln( fent, r ); //	CUSTO	MARGINAL DE	DEMANDA		SISTEMA	: SUL

    for iAeco := 0 to nAEcos -1 do
    begin
      readln( fent, r ); // --  línea en blanco --
      readln( fent, r ); //  A	NO: 2	11
      writeln( r );
      readln( fent, r ); //	PAT	1	2	3	4	5	6	7	8	9	10	11	12	MEDIA

      iCron:= 0;
      while
        ( cronicaInicioPrimerTira < 0 )
        or ( iCron <  (NCrons - (CronicaInicioPrimerTira - priCron ))) do
      begin
        for iPat:= 0 to NPats - 1 do
        begin
            readln( fent, r );
            pal:= trim(nextpal( r ));
            if pal <> '' then
              cronicaInicioTira:= StrToInt( pal ); // sino dejo el valor anterior.
            if cronicaInicioPrimerTira < 0 then
            begin
              cronicaInicioPrimerTira:= cronicaInicioTira;
              setlength( cmo, nPats );
              for kpat:= 0 to npats - 1 do
                cmo[kpat]:= TDatosCronicos.Create( priAEco, ultAEco, priCron, ultCron, cronicaInicioPrimerTira );
            end;

            kPat:= nextInt( r );
            if (kPat-1) <> iPat then
              raise Exception.Create( 'No coincide kPat con iPat ');

            for kmes:= 0 to 11 do
            begin
                pal:= nextpal( r );
                m:= StrToFloat( pal );
                imes:= kmes + iAeco * 12;
                cmo[ipat].cronicas[iCron][ imes ]:= m;
            end;
        end;
        inc( iCron );
      end;
      readln( fent, r ); // MEDI	A	0	0	0	0	16.56	19.42	27.94	27.72	38.04	48.93	40.14	46.18	33.12
    end;
    closefile( fent );

    // bien ahora hacemos el cuenterío
(*
    for ipat:= NPats -1 downto 1 do
      cmo[ipat].Calc_Diff( cmo[ipat-1] );
*)

    for ipat:= 0 to NPats -1 do
      cmo[ipat].Calc_Promedio;

    cmo[0].Determinar_MesesInicialesSinDatos;
    for ipat:= 1 to NPats -1 do
      cmo[ipat].nMesesInicialesSinDato:= cmo[0].nMesesInicialesSinDato;

    for ipat:= 0 to NPats -1 do
    begin
      cmo[ipat].Calc_Tendencia( 6, 5, false );
      cmo[ipat].QuitarTendencia;
      cmo[ipat].Calc_SerieMensual;
      cmo[ipat].Calc_SerieSemanal;
    end;

    // ahora guardamos la series y sus tendencias.
    assignfile( fsal , 'cmo_series_semanales.xlt' );
    rewrite( fsal );
    writeln( fsal, 'priCron:',#9, pricron, #9, 'ultCron:', ultCron );
    write( fsal, 'Año', #9, 'Semana' );

(*
    for ipat := 1 to NPats do
      if ipat = 1 then
        write( fsal, #9, 'cmo', ipat )
      else
        write( fsal, #9, 'cmo', ipat-1,'-',ipat );
        *)
    for ipat := 1 to NPats do
        write( fsal, #9, 'cmo', ipat );

    writeln( fsal );

    for isem:= 0 to high( cmo[0].serie_semanal ) do
    begin
      write( fsal, priCron + isem div 52, #9, (isem mod 52) +1 );
      for ipat:= 0 to high( cmo ) do
        write( fsal, #9, cmo[ipat].serie_semanal[isem]:12:6 );
      writeln( fsal );
    end;
    closefile( fsal );

    assignfile( fsal , 'cmo_tendencias.xlt' );
    rewrite( fsal );

    write( fsal, 'año',#9,'mes' );

(*
    for ipat := 1 to NPats do
      if ipat = 1 then
        write( fsal, #9, 'cmo', ipat )
      else
        write( fsal, #9, 'cmo', ipat-1,'-',ipat );
*)
    for ipat := 1 to NPats do
        write( fsal, #9, 'cmo', ipat );
    writeln( fsal );

    for imes:= 0 to high( cmo[0].tendencia ) do
    begin
      write(fsal, priAeco, #9, (imes mod 12)+1 );
      for ipat:= 0 to high( cmo ) do
        write( fsal, #9, cmo[ipat].tendencia[imes]:12:6 );
      writeln( fsal );
    end;
    closefile( fsal );
end;

end.

