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.
