unit uProbabilidadTransicion;

{$mode delphi}

interface

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

type

  // Lee la posición k y k+1 de X
  // Calcula en Z la transición Xk+1 - Suma(Ah*Xk-h).
  // En H marca con TRUE si hay un hueco y en cntHuecos queda la
  // cuenta de la cantidad de huecos.

  { TXH }

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


  { TTransicion_CEGH }

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

    InvSigma: TMatR;
    inv_divisor_pdf: NReal;

    constructor Create(cegh: TModeloCEGH );
    procedure Free;

  public


    kIni, kFin: integer; // índices del explorador de series.
    NPasosT: integer; // Transiciones a cubrir.

    function Calc_PDF_Serie( series: TSeriesDeDatos ): TVectR;
    procedure Calc_Zs( var Zs: TVectR;  Xs: TXH; X_ant: array of TXH);
    function Calc_DensProb( Zs: TVectR ): NReal;


  end;

implementation

{ TXH }

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

function TXH.ReadPos(series: TDAOfVectR; k: integer; 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 TXH.WritePos(series: TDAOfVectR; k: integer);
var
  j: integer;
begin
  for j := 1 to X.n do
    series[j - 1].pon_e(k, X.e(j));
end;

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

{ TTransicion_CEGH }

constructor TTransicion_CEGH.Create(cegh: TModeloCEGH );
var
  Ma: TMatR;
  detSigma: NReal;

begin
  inherited Create;
  self.cegh := cegh;
  self.A := cegh.A_cte;
  self.B := cegh.B_cte;

  kIni := 0;
  kFin := 0;

  InvSigma := TMatR.Create_init(cegh.B_cte.nf, cegh.B_cte.nc);
  Ma := TMatR.Create_Clone(cegh.B_cte);
  Ma.Transponer;
  InvSigma.Mult(cegh.B_cte, Ma);
  if not InvSigma.Inv(detSigma) then
    raise Exception.Create('Sigma no invertible!!!');
  Ma.Free;
  inv_divisor_pdf := 1/( power( 2*pi, invSigma.nf / 2.0 ) * sqrt( detSigma ) );
end;

procedure TTransicion_CEGH.Free;
begin
  InvSigma.Free;
  inherited Free;
end;

function TTransicion_CEGH.Calc_PDF_Serie( series: TSeriesDeDatos ): TVectR;
 var
  buscando: boolean;
  X_ant: array of TXH;
  Xs: TXH;
  k: integer;
  cntFirmes: integer;
  res: TVectR;
  cnt_huecos: integer;
  Zs: TVectR;
  kPaso: integer;

procedure ShiftEstado;
var
 t: TXH;
 j: integer;
begin
  t:= 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:= t;
end;

begin
  res:= TVectR.Create_init( series.NPuntos );

  setlength( X_ant, cegh.nRetardos );
  for k:= 0 to high( X_ant ) do
   X_ant[k] := TXH.Create( series.NSeries );
  Xs:=TXH.Create( series.NSeries );
  Xs.ReadPos( series.get_series_cronActiva, 1, -11000);

  cnt_huecos:= cegh.nRetardos* series.NSeries + Xs.cntHuecos;

  Zs:= TVectR.Create_Init( series.NSeries );

  for kPaso:= 1 to res.n do
  begin
    if cnt_Huecos = 0 then
    begin
      Calc_Zs( Zs, Xs, X_ant );
      res.pon_e(kPaso, Calc_DensProb( Zs ));
    end
    else
      res.pon_e(kPaso, -11000);

    shiftEstado;
    cnt_huecos:= cnt_huecos - Xs.cntHuecos;
    if kPaso < res.n then
      Xs.ReadPos( series.get_series_cronActiva, kPaso+1, -11000);
    cnt_huecos:= cnt_huecos + Xs.cntHuecos;

  end;

  for k:= 0 to high( X_ant ) do
   X_ant[k].free;
  setlength( X_ant, 0 );
  Xs.Free;
  Zs.Free;

  result:= res;
end;

procedure TTransicion_CEGH.Calc_Zs( var Zs: TVectR;  Xs: TXH; X_ant: array of TXH);
var
    k, j, jRetardo, jBase: integer;
    m: NREal;

begin
   for k:=1 to Zs.n do
   begin
     m:= 0;
     For jRetardo:=0 to high (X_ant) do
     begin
       jBase := jRetardo * Zs.n;
       for j := 1 to Zs.n do
         m:= m + A.e(k, jBase + j) * X_ant[jRetardo].X.e(j);
     end;
     m:= Xs.X.e(k) - m;
     Zs.pon_e( k, m );
    end;
end;


function TTransicion_CEGH.Calc_DensProb(Zs: TVectR): NReal;
var
  cons: NReal;
  aux: NReal;
  v: TVectR;

begin
 v:= TVectR.Create_init( InvSigma.nf );
 InvSigma.Transformar( v, Zs );
 aux:= exp( - 0.5 * Zs.PEV(v) );
 v.Free;
 Result:= aux* inv_divisor_pdf;
end;

end.

