unit uestimadoryrellenado;
{$mode delphi}
interface

uses
  Classes, SysUtils, xmatdefs,
  umodelosintcegh, useriestemporales,
  // usparsematreal,
  matent, matreal, matbool;

type
  // Lee de series la posición k y la guarda en X
  // En H marca con TRUE si hay un hueco y en cntHuecos queda la
  // cuenta de la cantidad de huecos.

  { TXHRec }

  TXHRec = class
    X: TVectR;
    H: TVectBool;
    cntHuecos: integer;
    constructor Create(NSeries: integer);
    function ReadPos(series: TDAOfVectR; const k: integer; const umbral_Hueco: NReal): integer;
    procedure WritePos(series: TDAOfVectR; const k: integer);
    procedure Free; virtual;
  end;



  { TRelleandorDeHuecos_CEGH }

  TRelleandorDeHuecos_CEGH = class
    umbral_hueco: NReal;
    cegh: TModeloCEGH;
    A, B: TMatR; // matrices del CEGH X_k+1 = sum( Ah, X_k-h ) + B R_k

    constructor Create( cegh: TModeloCEGH; umbral_hueco: NReal);
    procedure Free;
    // revisa las series que supone son las que corresponden al vecto X
    // y completa los huecos
    function Rellenar(series: TDAofVectR; MaxNPasosT: integer;
      series_ptr: TSeriesDeDatos): integer;

  private
    NRetardos: integer; // A.nc div A.nf se calcula en el Crete
    series: TDAOfVectR; // puntero a las series pasadas en Rellenar.
    NSeries: integer; // cantidad de series
    NDatos: integer; // largo de las series
    kIni, kFin: integer; // índices del explorador de series.
    NPasosT: integer; // Transiciones a cubrir.
    //MA, MB: TSparseMatR; // para armar el problema

    MA, MB: TMatR; // para armar el problema

    top_h, left_h: TVectBool; // marca los huecos
    pemp_A: TMatR;
    pemp_B: TVectR;

    // busca un tramo con huecos entre dos pasos completos
    // retorna True si encontró, y False si ya no quedan
    function buscarNuevoTramo: boolean;

    // Armado Problema Entrada Más Probable
    // Usa Tramo para armar MA y MB
    procedure ArmardoProblemaEntradaMasProbable;

    // Sustituye en la variable explicada por la fila kFil
    procedure SustituirFilaVar(kFil: integer);

    // Intercambia filas de MA y MB y los casilleros correspondientes
    // de left y left_h
    procedure SwapFilas(k1, k2: integer);


    function xAt_jCol(jCol: integer): NReal;
    function xAt_kFil(kFil: integer): NReal;

    procedure Calc_XS(Xs: TXHRec; X_ant: array of TXHRec; R: TVectR; iPaso: Integer
      );
    procedure CompletarSeries(R: TVectR);

  end;

(** rch@201602071749
  Encuentra el R de máxima verosimilitud para la expresión A R = B
    R = A' inv( A A' ) B
  se supone que R es un vector de gaussianas independientes con pdf = N(0,1)
   **)
function EntradaMasProbable(A: TMatR; B: TVectR): TVectR;




implementation

function EntradaMasProbable(A: TMatR; B: TVectR): TVectR;
var
  M, N, M_adj, M_inv: TMatR;
  va: NReal;
  flg_invertible: boolean;
  det: NReal;
  exp10: integer;
  R: TVectR;
  k, j: integer;

begin

//  A.WriteM;
//  A.WriteArchiXLT( 'c:\basura\pemp_A.xlt' );
//  B.WriteXLTSimple_archi('c:\basura\pemp_B.xlt' );
  // M = A A'
  M := TMatr.Create_Init(A.nf, A.nf);
  for k := 1 to M.nf do
  begin
    va := A.Fila(k).PEV(A.Fila(k)); // Diagonal
    M.pon_e(k, k, va);
    for j := k + 1 to M.nc do  // Laterales
    begin
      va := A.Fila(k).PEV(A.Fila(j));
      M.Pon_e(k, j, va);
      M.Pon_e(j, k, va);
    end;
  end;

  // N = inv( M ) B
  N := TMatR.Create_init(B.n, 1);
  N.CopyVectCol(B, 1);
  det := M.Escaler(N, flg_invertible, exp10);
  if not flg_invertible then
  begin
    M.Free;
    N.Free;
    raise Exception.Create('EntradaMasProbable ... AA'' no invertible ');
  end;

  //  R = A' N
  R := TVectR.Create_init(A.nc);
  for k := 1 to R.n do
  begin
    va := 0;
    for j := 1 to A.nf do
      va := va + A.e(j, k) * N.e(j, 1);
    R.pon_e(k, va);
  end;
  M.Free;
  N.Free;
  Result := R;
end;


(*** Métodos de TXHRec ***********)


constructor TXHRec.Create(NSeries: integer);
begin
  inherited Create;
  X := TVectR.Create_init(NSeries);
  H := TVectBool.Create_init(NSeries);
end;

function TXHRec.ReadPos(series: TDAOfVectR; const k: integer;
  const umbral_Hueco: NReal): integer;
var
  j: integer;
  a: NReal;
begin
  cntHuecos := 0;
  for j := 1 to X.n do
  begin
    a := series[j - 1].e(k);
    if a <= umbral_hueco then
    begin
      H.pon_e(j, True);
      Inc(cntHuecos);
    end
    else
      H.pon_e(j, False);
    X.pon_e(j, a);
  end;
  Result := cntHuecos;
end;



procedure TXHRec.WritePos(series: TDAOfVectR; const k: integer);
var
  j: integer;
begin
  for j := 1 to X.n do
    if  H.e(j) then
    series[j - 1].pon_e(k, X.e(j));

end;



procedure TXHRec.Free;
begin
  X.Free;
  H.Free;
  inherited Free;
end;


(*** Métodos de TRellenadorDeHuecos_CEGH *****)
constructor TRelleandorDeHuecos_CEGH.Create(cegh: TModeloCEGH;
  umbral_hueco: NReal);
begin
  inherited Create;
  self.umbral_hueco := umbral_hueco;
  self.cegh:= cegh;
  self.A := cegh.A_cte;
  self.B := cegh.B_cte;
  MA := nil;
  MB := nil;
  top_h:= nil;
  left_h:= nil;
  pemp_A:= nil;
  pemp_B:= nil;


  kIni := 0;
  kFin := 0;


  NRetardos := A.nc div A.nf;
end;

procedure TRelleandorDeHuecos_CEGH.Free;
begin
  if MA <> nil then
    MA.Free;
  if MB <> nil then
    MB.Free;
  if Top_h <> nil then
    Top_h.Free;
  if Left_h <> nil then
    Left_h.Free;
  if pemp_A <> nil then
    pemp_A.Free;
  if pemp_B <> nil then
    pemp_B.Free;
  inherited Free;
end;

procedure TRelleandorDeHuecos_CEGH.Calc_XS(Xs: TXHRec; X_ant: array of TXHRec;
  R: TVectR; iPaso: Integer);
var
  k, j, jRetardo, jBase, jBase_ruido: integer;
begin
  jBase_ruido:= R.n-iPaso*NSeries;
  for k := 1 to NSeries do
  begin
    Xs.X.pon_e(k, 0);
    for jRetardo := 0 to high(X_ant) do
    begin
      jBase := jRetardo * NSeries;
      for j := 1 to NSeries do
        Xs.X.acum_e(k, A.e(k, jBase + j) * X_ant[jRetardo].X.e(j));
    end;
    for j := 1 to B.nc do
      Xs.X.acum_e(k, B.e(k, j) * R.e(jBase_ruido + j));
  end;
end;

procedure TRelleandorDeHuecos_CEGH.CompletarSeries(R: TVectR);
var
  X_ant: array of TXHRec;
  Xs: TXHRec;
  k, j: integer;
  kPasoT: integer;
  rt: TXHRec;
begin
  setlength(X_ant, NRetardos);
  for k := 0 to NRetardos - 1 do
  begin
    x_ant[k] := TXHRec.Create(NSeries);
    x_ant[k].ReadPos(series, kIni + k, umbral_hueco);
  end;

  xs := TXHRec.Create(NSeries);

  for kPasoT := kIni + NRetardos to kFin do
  begin
    xs.ReadPos(series, kPasoT, umbral_hueco);
    Calc_Xs(Xs, X_ant, R, kPasoT-(kIni + NRetardos)+1);
    Xs.WritePos(series, kPasoT);
    rt := X_ant[0];
    for j := 0 to high(X_Ant) - 1 do
      X_ant[j] := X_ant[j + 1];
    X_ant[high(X_ant)] := xs;
    xs := rt;
  end;
  for j := 0 to high(X_Ant) do
    x_ant[j].Free;
  setlength(X_ant,0);
  xs.Free;
end;

// revisa las series que supone son las que corresponden al vecto X
// y completa los huecos
function TRelleandorDeHuecos_CEGH.Rellenar(series: TDAofVectR;
  MaxNPasosT: integer; series_ptr: TSeriesDeDatos): integer;
var
  pemp_R: TVectR;
  cnt: integer;
begin
  self.Series := series;
  NSeries := length(series);
  NDatos := series[0].n;

  cnt := 0;
  while buscarNuevoTramo do
  begin
    if ( NPasosT <= MaxNPasosT ) then
    begin
      Inc(cnt);
      if cnt = 14 then
         writeln('hola');
      writeln( 'Rellenando tramo: ', cnt, ', kIni: ', kIni, ', kFin: ', kFin, ', NPasosT: ', NPasosT );
    if cnt = 14 then
    series_ptr.WriteToArchi('c:\basura\series_gaussianas_rellenadas_paso_'+IntToStr( cnt )+'_1.txt');
      ArmardoProblemaEntradaMasProbable;
    if cnt = 14 then
    series_ptr.WriteToArchi('c:\basura\series_gaussianas_rellenadas_paso_'+IntToStr( cnt )+'_2.txt');
      {$IFDEF VERBOSO}
      writeln('Armando .... pemp');
      {$ENDIF}

      pemp_R := EntradaMasProbable(pemp_A, pemp_B);
    if cnt = 14 then
    series_ptr.WriteToArchi('c:\basura\series_gaussianas_rellenadas_paso_'+IntToStr( cnt )+'_3.txt');

      {$IFDEF VERBOSO}
      writeln('Completando series' );
      {$ENDIF}
      CompletarSeries(pemp_R);
    if cnt = 14 then
    series_ptr.WriteToArchi('c:\basura\series_gaussianas_rellenadas_paso_'+IntToStr( cnt )+'_4.txt');
      pemp_R.Free;

      series_ptr.WriteToArchi('c:\basura\series_gaussianas_rellenadas_paso_'+IntToStr( cnt )+'.txt');
    end
    else
      writeln( '.... salteando tramo NPasos. '+IntToStr( NPasosT ) );
  end;
  Result := cnt;
end;

// busca un tramo con huecos entre dos pasos completos
// retorna True si encontró, y False si ya no quedan
function TRelleandorDeHuecos_CEGH.buscarNuevoTramo: boolean;
var
  buscando: boolean;
  r: TXHRec;
  k: integer;
  cntFirmes: integer;

begin
  NPasosT := 0;
  r := TXHRec.Create(length(series));
  if kFin = 0 then
  begin
    // Si es la primera vez busco "base X firme"
    k := kFin;
    writeln(kFin);
    buscando := True;
    cntFirmes := 0;
    while buscando and (k < NDatos) do
    begin
      Inc(k);
      r.ReadPos(series, k, umbral_hueco);
      if r.cntHuecos = 0 then
      begin
        Inc(cntFirmes);
        buscando := cntFirmes < NRetardos;
      end
      else
        cntFirmes := 0;
    end;
    if buscando then
    begin
      r.Free;
      Result := False;
      exit;
    end
    else
    begin
      kIni := k;
    end;
  end
  else
  begin
    kIni := kFin;
  end;


  // Llegados aquí, kIni es el último de un conjunto de NRetardos X firmes.
  // ahora bauscamos el primer pozo.
  buscando := True;
  while buscando and (kIni < NDatos) do
  begin
    r.ReadPos(series, kIni, umbral_hueco);
    if kIni=kFin then
      Assert(r.cntHuecos=0, 'TRelleandorDeHuecos_CEGH.buscarNuevoTramo: ...');
    if r.cntHuecos > 0 then
      buscando := False
    else
      Inc(kIni);
  end;

  if buscando then
  begin
    r.Free;
    Result := False; // no hay más huecos
    exit;
  end;


  kFin := kIni;
  cntFirmes := 0;
  buscando := True; // buscamos un final firme
  while buscando and (kFin < NDatos) do
  begin
    Inc(kFin);
    r.ReadPos(series, kFin, umbral_hueco);
    if r.cntHuecos = 0 then
    begin
      //si hay un bloque de firmes menor a NRetardos, lo mete adentro del tramo de huecos
      Inc(cntFirmes);
      buscando := cntFirmes < NRetardos;
    end
    else
      cntFirmes := 0;
  end;

  if buscando then
    Result := False // no se encontró base firme al final
  else
  begin
    // Si encontramos base firme ajustamos kIni para comenzar al inicio
    // del primer paso firme y copiamos todos los registros

    NPasosT := kFin - kIni + 1;
    kIni := kIni - NRetardos;

    //kIni := kIni - NRetardos;
    //NPasosT := kFin - kIni + 1 - NRetardos;


    Result := True;
  end;
  r.Free;

end;


// Sustituye en la variable explicada por la fila kFil
procedure TRelleandorDeHuecos_CEGH.SustituirFilaVar(kFil: integer);
var
  jCol: integer;
  k, j: integer;
  m: NReal;
  iPasoT: integer;
  jSerie, jCol1: integer;
begin
  iPasoT:= (kFil -1 ) div NSeries;

  jSerie:= ((kFil -1 ) mod NSeries) + 1;

  jCol := ( NPasosT - iPasoT -2 ) * NSeries + jSerie;

  jCol1:= MA.nc - (NRetardos+1 +iPasoT)*NSeries + jSerie ;

//  writeln('iPasoT ',iPasoT,'jSerie ',jSerie,'kFil ',kFil,'jCol ',jCol,'jCol1 ',jCol1);

  for k := 1 to MA.nf do
  begin
    if k <> kFil then
    begin
      m := MA.e(k, jCol);
      if abs(m) > 1e-10 then
      begin
        for j := 1 to MA.nc do
          if j <> jCol then
            MA.acum_e(k, j, MA.e(kFil, j) * m);
        for j := 1 to MB.nc do
          MB.acum_e(k, j, MB.e(kFil, j) * m);
      end;
      MA.pon_e(k, jCol, 0);
    end;
  end;
end;

// Intercambia filas de MA y MB y los casilleros correspondientes
// de left y left_h
procedure TRelleandorDeHuecos_CEGH.SwapFilas(k1, k2: integer);
begin
  MA.IntercambieFilas(k1, k2);
  MB.IntercambieFilas(k1, k2);
end;


function TRelleandorDeHuecos_CEGH.xAt_jCol(jCol: integer): NReal;
var
  kRec, iSerie: integer;
begin
  iSerie := (jCol - 1) mod NSeries;
  kRec := kIni + ((MA.nc - jCol) div NSeries);
  Result := series[iSerie].e(kRec);
end;


function TRelleandorDeHuecos_CEGH.xAt_kFil(kFil: integer): NReal;
var
  kRec, iSerie: integer;
begin
  iSerie := (kFil - 1) mod NSeries;
  kRec := kIni + NRetardos + (kFil - 1) div NSeries;
  Result := series[iSerie].e(kRec);
end;

procedure TRelleandorDeHuecos_CEGH.ArmardoProblemaEntradaMasProbable;
var
  NX, MX, NR, MR: integer;
  kRec, k, j: integer;
  r: TXHRec;
  jColBase, kFilBase, jCol, kFil: integer;
  jColBaseMA, jColBaseMB: integer;
  cntRestricciones: integer;
  kPasoT: integer;

begin
  if MA <> nil then
    MA.Free;
  if MB <> nil then
    MB.Free;
  if Top_h <> nil then
    Top_h.Free;
  if Left_h <> nil then
    Left_h.Free;
  if pemp_A <> nil then
    pemp_A.Free;
  if pemp_B <> nil then
    pemp_B.Free;

  MX := NPasosT * A.nf;
  NX := (NPasosT - 1) * A.nf + A.nc;

  MR := NPasosT * B.nf;
  NR := NPasosT * B.nc;

  assert(MR = MX, 'Atención A.nf <> B.nf no es posible!!');
  assert(A.nf = length(series), 'Atención, A.nf <> length( series )!!');
  assert(A.nf = B.nf, 'Atención, A.nf <> B.nf!!');
  //MA := TSparseMatR.Create_Init(MX, NX);
  //MB := TSparseMatR.Create_Init(MR, NR);

  MA := TMatR.Create_Init(MX, NX);
  MB := TMatR.Create_Init(MR, NR);

  top_h := TVectBool.Create_init(MA.nc);
  left_h := TVectBool.Create_init(MA.nf + 1);

  {$IFDEF VERBOSO}
  writeln( 'NPasosT: ', NPasosT );
  {$ENDIF}
  // ahora leemos todos las Fichas involucradas para determinar cuales
  // tienen huecos y los anotamos en top_h y left_h
  r := TXHRec.Create(NSeries);
  // Primer Tramo solo top_h
  for kRec := kIni to (kIni + NRetardos - 1) do
  begin
    r.ReadPos(series, kRec, umbral_hueco);
    jColBase := MA.nc - (kRec - kIni + 1) * NSeries;
    for jCol := 1 to NSeries do
      top_h.pon_e(jColBase + jCol, r.H.e(jCol));
  end;

  // Tramo Central top_h y left_h
  for kRec := kIni + NRetardos to (kFin - 1) do
  begin
    r.ReadPos(series, kRec, umbral_hueco);
    jColBase := MA.nc - (kRec - kIni + 1) * NSeries;
    kFilbase := (kRec - (kIni + NRetardos)) * NSeries;
    for j := 1 to NSeries do
    begin
      top_h.pon_e(jColBase + j, r.H.e(j));
      left_h.pon_e(kFilBase + j, r.H.e(j));
    end;
  end;

  // Ultimo Tramo sólo left_h
  for kRec := kFin to kFin do
  begin
    r.ReadPos(series, kRec, umbral_hueco);
    kFilbase := (kRec - (kIni + NRetardos)) * NSeries;
    for j := 1 to NSeries do
    begin
      left_h.pon_e(kFilBase + j, r.H.E(j));
    end;
  end;
  r.Free;

  {$IFDEF VERBOSO}
  writeln( 'Estampillando A y B' );
  {$ENDIF}
  // Ahora estampillamos las matrices A y B en la diagonl inversa de MA y MB
  for kPasoT := 1 to NPasosT do
  begin
    kFilBase := (kPasoT - 1) * A.nf;
    jColBaseMA := (NPasosT - kPasoT) * NSeries;
    jColBaseMB := (NPasosT - kPasoT) * B.nc;
    for kFil := 1 to A.nf do
    begin
      for jCol := 1 to A.nc do
        MA.pon_e(kFilBase + kFil, jColBaseMA + jCol, A.e(kFil, jCol));
      for jCol := 1 to B.nc do
        MB.pon_e(kFilBase + kFil, jColBaseMB + jCol, B.e(kFil, jCol));
    end;
  end;


  // [X_0+NRetardos, X_1+NRetardos, ... Xn+1] = MA [Xn,  ...X2, X1, X0] + MB [R0, R1, ... Rn]
  // Ahora recorremos lso Xs y donde hay un hueco, sustituimos
  // esa variable en los X
  // Si el hueco está en la fila k, la variable a sustituir está en la
  // columna k+ A.nc
  cntRestricciones := 0;
  for kFil := 1 to MA.nf do
    if left_h.e(kFil) then
    begin
      {$IFDEF VERBOSO}
      writeln( 'Sustituyenco kFil: '+IntToStr( kFil ) );
      {$ENDIF}
      SustituirFilaVar(kFil);
    end
    else
      Inc(cntRestricciones);

  // Ahora, para terminar de armar el problema de Entrada Más probable tenemos
  // que la matriz A del problema son las filas de MB (las que no explican huecos)
  // y la matriz B del problema es calculable
  // como I - MA X sobre las filas de MA y(que no explican huecos) teniendo
  // en cuenta que las columnas correspondientes a las X desconocidas tienen
  // los coeficients nulos. Se MI el término.
  pemp_A := TMatR.Create_init(cntRestricciones, MB.nc);
  pemp_B := TVectR.Create_init(cntRestricciones);
  k := 1;
  for kFil := 1 to pemp_A.nf do
  begin
    {$IFDEF VERBOSO}
    writeln( 'completando kFil: '+IntToSTr( kFil ) );
    {$ENDIF}
    while left_h.e(k) do
    begin
      if k < left_h.n then
       Inc(k)
      else
       raise Exception.Create( 'NO puede ser ... left_h['+IntToStr(k)+']= TRUE' );
    end;

    for jCol := 1 to pemp_A.nc do
      pemp_A.pon_e(kFil, jCol, MB.e(k, jCol));

    pemp_B.pon_e(kFil, XAt_kFil(k));
    for jCol := 1 to MA.nc do
    begin
      if not top_h.e(jCol) then
        pemp_B.acum_e(kFil, -MA.e(k, jCol) * XAt_jCol(jCol));
    end;
    Inc(k);
  end;
end;

end.
