{$IFDEF DUODINAMICO}
{$I upronostico_DUODINAMICO_.pas}
{$ELSE}
unit upronostico;

interface

uses
  Classes, SysUtils, xMatDefs, matreal,
  ufechas, ucosa, Math, umatriz_ruida,
  umodelosintcegh, uglobs;

const
  C_DT_1MINUTO = 1.0 / 24.0 / 60.0;

type


  TMatricesRuida = class;

  // definición del Cono de Pronósticos para una serie temporal.
  // La "guia" es la especificación del Cono de Pronósticos en el Espacio de la Señal.
  // en base a la guía se calculan el vector "sesgos" y "factor" que contiene
  // que se aplican en el Espacio Gaussiano para guiar la salida del Filtro Lineal.
  // el ruido que se aplica sobre el canal se calcula como:
  // rs = sesgo[j] + factor[j] * Brbg[i]
  // donde "Brbg" es el Ruido Blanco Gaussiano transformado por la matriz B
  // del filtro lineal que corersponde al canal en cuestión
  // El control del cono de pronósticos se realiza mediante una "guia" del cono
  // e indicaciones de cómo abrir el cono de dispersión.
  // El vector "guia" indica valores para la guía. De ese vector, se utilizan los
  // primeros valores para inicializar el estado del filtro lineal.
  // El parámetro NPSA = Número de Pasos Sin Apertura del cono, determina la
  // cantidad de valores de la guía, posteriores a los usados para fijar el estado
  // inicial que serán utilizados como valores determinísticos.
  // Ese valor puede ir entre 0 hasta el total de valores que queden en la guía.



  { TPronostico }
  T_aux_pronostico = class  // auxilar para lectura
    p: NReal;
    nretardos: integer;
    Serie: string;
    _fechaIniSesgo_: TCosa;
  end;

  TPronosticos = class;

  { TPronostico
    Información del Cono de Pronóstico de una Serie de un modelo CEGH.
  }
  TPronostico = class(TCosa)

  private
    rangoFechaSesgo: NReal; // rango de aplicación del sesgo.
    rangoFechaFactor: NReal; // rango de aplicación del ruido

  public
    (**************************************************************************)
    (*              A T R I B U T O S   P E R S I S T E N T E S               *)
    (**************************************************************************)

    Serie_: string; // Nombre de la serie a la que esta asociado el pronostico
    kSerie_: integer; // A resolver en PrepararMemoria
    nRetardos_: integer;

    // Número de valores de la guía que serán usados para inicializar el filtro lineal
    guia: TDAOfNReal; // Guía en el espacio de la señal real
    guia_eg: TDAOfNReal; // Guía en el espacio gasussiano de la señal.
    guia_pe: NReal;      // probabilidad de excedencia de la guía.
    sesgo: TDAOfNReal; // sesgo del ruido (aditivo)
    factor: TDAOfNReal; // multiplicadores del ruido (multimplicativo)

    // parámetros de calibración control del cono de pronóstico a partir de los determinismos.
    NPCC: integer; // Número de Pasos de Control del Cono
    NPLC: integer; // Número de Pasos de Liberación del Cono
    NPSA: integer; // Número de pasos Sin Apertura del Cono (determinista).
    NPAC: integer; // Número de Pasos de Apertura del Cono

    AperturaInicial_pu: NReal; // Apertura Inicial del cono en p.u.


    // información para consulta a servicio web de pronósticos.
    url_get: string; // ej: http://simsee.org/pronosticos
    nombre_get: string; // ej: PEOL_1MW

    escenario: TPronosticos;


    constructor Create(capa: integer; xNPCC, xNPLC, xNPSA, xNPAC: integer;
      xurl_get, xnombre_get: string);
    constructor Create_Default(cegh: TModeloCEGH);

    function Rec: TCosa_RecLnk; override;
    procedure BeforeRead(version, id_hilo: integer); override;
    procedure AfterRead(f:TArchiTexto); override;

    procedure Cambiar_GUIA(xNPCC, xNPLC, xNPSA, xNPAC: integer; xNuevaGuia: TDAOfNReal);
    procedure Free; override;
    function determinismoSoloEstadoInicial: boolean;

    // retorna max( NPCC+NPLC , NPSA+NPAC)
    // es útil para determinar a partir de qué paso el cono ya no tiene influencia.
    function NPasosCono: integer;

    // debe ser llamado en el preparar memoria de la fuente
    procedure calcular_rangos(DuracionPasoSorteoEnHoras: NReal);

    // determina si para la fechaDeInicioDelPaso pasada como parámetro es válido el cono
    // depronósticos o no. Si el resultado es TRUE indica que la fecha está en el rango del cono
    // y entonces en kSesgo y/o kFactor se devuelve los índices para leer los valores correspondientes
    // de los vectores "sesgo" y "factor".
    function fechaEnRango(
      fechaDeInicioDelPaso: TFecha;
      var kSesgo, kFactor: integer): boolean;


    // extiende con CEROS las guias de pronósticos en espacios gaussianos.
    // la hipótesis es que luego de finalizado el cono, la guía es el CERO
    // en el espacio gaussiano.
    function fguia_extendida_(aGuia: TDAOfNReal; kPaso: integer; DefVal: NReal): NReal;

  private
    px: T_aux_pronostico; // auxiliar para lectura
    fechaIniCono_dummy: TFecha;

  end;


  { TPronosticos
  Conjunto de Pronosticos, uno por SERIE asociados a una probabilidad
  de ocurrencia  Probabilidad.
  }

  TPronosticos = class(TListaDeCosas)
  private

    datosModelo: TModeloCEGH;
    globs: TGlobs;
    matricesRuida: TMatricesRuida;

  public

    Probabilidad: NReal;
    fechaIniCono: TFecha; // fecha a partir de la cual aplica el sesgo

    Bamps:TList;

    constructor Create_Default(cegh: TModeloCEGH);
    procedure prepararse( DuracionPasoSorteoEnHoras: NReal );
    // retorna el máximo largo en pasos de sorteo de la fuente de los conos
    // de pronósticos represnetado.
    function max_NPasosCono: integer;
    constructor Create_ReadFromText(f: TArchiTexto); override;
    procedure WriteToText_(f: TArchiTexto); override;
    function getItem(i: integer): TPronostico;
    procedure setItem(i: integer; cosa: TPronostico);
    property items[i: integer]: TPronostico read getItem write setItem; default;
    procedure Free; override;
    function NSeries: integer;


    procedure CalcularMatricesRuida(datosModelo: TModeloCEGH; globs: TGlobs);

    {$IFDEF DUODINAMICO}
    function GetBamp (kPaso: Integer):TMatR;
    {$ENDIF}
    function Guia_eg (kPaso: Integer):TVectR;
    function Guia (kPaso: Integer):TVectR;
    function Sesgos (kPaso: Integer):TVectR;
    function GetLambda (kPaso: Integer): TMatR;

    function getLargoGuia:Integer;

    // El inicio del pronostico no tiene porque coicidir con el inicio del
    // horizonte de opt/sim. Esta funcion devuelve kPasoProno correspondiente
    // al kPaso de opt o sim.
    // Ej: Si fechaIniSim<fechaIniCono -> kPasoPronos=0
    function kPasoPronos(kPaso: Integer):Integer;

  end;


  { TEscenarioDePronosticos
  Conjunto de TPronosticos.
  }

  TEscenarioDePronosticos = class(TListaDeCosas)
  private

    function getItem(i: integer): TPronosticos;

  public

    constructor Create(capa: integer); override;overload;

    function Add(Pronosticos: TPronosticos; Probabilidad: NReal): integer;
    property items[Index: integer]: TPronosticos read getItem; default;

    // 0<= xP <= 1.0
    // Retorna un pronóstico tal que la probabilidad acumulada de los pronósticos
    // de la lista de escenarios hasta el devuelto inclusive supera o iguala xP
    // y sin considerar el devuelto no la supera.
    // Por ej. si el primer escenario tiene probabilidad 0.2 será devuelto siempre
    // que xP <= 0.2
    function GetEscenarioPorP_(xP: NReal): TPronosticos;

    procedure prepararse(DuracionPasoSorteoEnHoras: NReal);

  public

    procedure CalcularMatricesRuida(datosModelo: TModeloCEGH; globs: TGlobs);

  end;

  { TMatricesRuida }

  TMatricesRuida = class
  private

    modelo: TModeloCEGH;
    escenarioPronos: TPronosticos;

    // Filtro:     X_k+1 = A * X_k + B * R_k
    // Reduccion:  Z_k = M_r * X_k
    // Ampliacion: X_k = M_a * Z_k + B_a * W_a

    // Matrices del filtro
    A, B: TMatR;
    // Matrices de reduccion y ampliacion
    Mr, Ma: TMatR;

    MatrizRuidaMultiRetardos: TMatR;

    //*********** AUXILIARES ********
    // A'
    At,
    // M_a * M_r
    MaMr,
    // (M_a * M_r)'
    MaMr_t,
    // B * B'
    BBt,

    //***********

    BampEE: TMatR;

 {$IFDEF DUODINAMICO}
   lst: TList;
 {$ENDIF}

  public

    constructor Create(aModelo: TModeloCEGH; aEscenarioPronos: TPronosticos);

    (*

    Filtro con pronosticos: X_k+1 = A * X_k + lambda * B * W_k + S_k   ec.1

    Reduccion: Y_k = M_r * X_k  ec.2

    Amplicacion: X_k = M_a * Y_k + Bamp_k * Wamp_k   ec.3

    Aplicando reduccion al filtro:

    Y_k+1 = M_r * ( A * X_k + lambda * B * W_k + S_k ) ec.4

    Sustituyendo ec.4 en ec.3

    X_k+1 = M_a * M_r * ( A * X_k + lambda * B * W_k + S_k ) + Bamp_k+1 * Wamp_k+1 ec.5

    Queremos Bamp_k+1 que recompone la matriz de covarianzas del estado X_k+1


    Def. E([x-E(x)][x-E(x)]')

    Calculando la matriz de covarianzas de ec.1

    <X_k+1 * X_k+1'>  = A * <X_k * X_k'> * A' + lambda_k * B * B' * lambda_k' + A * <X_k> * <X_k'> * A' + S_k * <X_k'> * A' + A * <X_k> * S_k' ec.6

    Calculando la matriz de covarianzas de ec.5

    <X_k+1 * X_k+1'>  = M_a * M_r  * [ A * <X_k * X_k'> * A' + lambda_k * B * B' * lambda_k' + A * <X_k> * <X_k'> * A' + S_k * <X_k'> * A' + A * <X_k> * S_k' ] * M_r' * M_a' + Bamp_k * Bamp_k'  ec.7

    Igualando ec.6 y ec.7, despejando Bamp * Bamp' podemos hallar Bamp_k por algún método.

    *)

    procedure calcularse(globs: TGlobs);

    {$IFDEF DUODINAMICO}
    function Get_Bamp(kPaso: integer): TMatR;
    function Count: integer;
    {$ENDIF}
  end;

const
  ERROR = 1.0e-5;


procedure AlInicio;
procedure AlFinal;

implementation

uses
  uversion_architexto;

procedure TPronosticos.CalcularMatricesRuida(datosModelo: TModeloCEGH;
  globs: TGlobs);
var
  aModelo: TModeloCEGH;
begin

  self.globs:=globs;
  self.datosModelo:=datosModelo;

  matricesRuida:=TMatricesRuida.Create(datosModelo, Self);
 {$IFDEF RUIDA_PRONOS}
   matricesRuida.calcularse(globs);
 {$ELSE}
   lst.Add(Matriz_RUIDA_EE(datosModelo.A_cte, datosModelo.B_cte,
     datosModelo.MAmp_cte, datosModelo.MRed_cte));
 {$ENDIF}
end;



{ TEscenarioDePronosticos }

function TEscenarioDePronosticos.getItem(i: integer): TPronosticos;
begin
  Result := lst.items[i];
end;

constructor TEscenarioDePronosticos.Create(capa: integer);
begin
  inherited Create(capa, '');
end;

function TEscenarioDePronosticos.Add(Pronosticos: TPronosticos;
  Probabilidad: NReal): integer;
begin
  Result := inherited Add(Pronosticos);
  Pronosticos.Probabilidad := Probabilidad;
end;

function TEscenarioDePronosticos.GetEscenarioPorP_(xP: NReal): TPronosticos;
var
  i: integer;
  PAcum: NReal;
begin
  if (xP < 0) or (xP > 1) then
    raise Exception.Create('TEscenarioDePronosticos.GetEscenarioPorP: 0<=xP<=1, xP=' +
      FloatToStrF(xP, ffGeneral, 3, 4));
  PAcum := 0;
  for i := 0 to self.Count - 1 do
  begin
    PAcum := PAcum + (items[i] as TPronosticos).Probabilidad;
    if PAcum >= xp then
    begin
      Result := items[i] as TPronosticos;
      break;
    end;
  end;
end;

procedure TEscenarioDePronosticos.prepararse(DuracionPasoSorteoEnHoras: NReal);
var
  aPronosticos: TPronosticos;
  k: integer;
begin
  for k := 0 to Count - 1 do
  begin
    aPronosticos := items[k];
    aPronosticos.prepararse(DuracionPasoSorteoEnHoras);
  end;
end;

procedure TEscenarioDePronosticos.CalcularMatricesRuida(
  datosModelo: TModeloCEGH; globs: TGlobs);
var
  k: Integer;
begin

  // Matrices de reconstruccion
  for k:=0 to self.Count-1 do
    getItem(k).CalcularMatricesRuida(datosModelo, globs);

end;

{$IFDEF DUODINAMICO}
function TPronosticos.GetBamp(kPaso: Integer): TMatR;
begin
  {$IFDEF RUIDA_PRONOS}
  if (kPasoPronos(kPaso)<matricesRuida.Count) then
    Result:=matricesRuida.Get_Bamp(kPasoPronos(kPaso))
  else
    Result:=matricesRuida.lst.Last;
  {$ELSE}
     Result:= matricesRuida.Get_Bamp(0)
  {$ENDIF}
end;
{$ENDIF}

function TPronosticos.Guia_eg(kPaso: Integer): TVectR;
var
  i: Integer;
begin
  Result:=TVectR.Create_Init(Self.Count);
  for i := 1 to Result.n do
    Result.pon_e(i, getItem(i-1).fguia_extendida_(getItem(i-1).guia_eg, kPasoPronos(kPaso), 0));
end;

function TPronosticos.Guia(kPaso: Integer): TVectR;
var
  i: Integer;
begin
  Result:=TVectR.Create_Init(Self.Count);
  for i := 1 to Result.n do
    Result.pon_e(i, getItem(i-1).fguia_extendida_(getItem(i-1).guia, kPasoPronos(kPaso), 0));
end;

function TPronosticos.Sesgos(kPaso: Integer): TVectR;
var
  i: Integer;
begin
  Result:=TVectR.Create_Init(Self.Count);
  for i := 1 to Result.n do
    Result.pon_e(i, getItem(i-1).fguia_extendida_(getItem(i-1).sesgo, kPasoPronos(kPaso), 0));
end;

function TPronosticos.GetLambda(kPaso: Integer): TMatR;
var
  i: Integer;
begin
  Result:=TMatR.Create_Init(Self.Count, self.Count);
  Result.Ceros;
  for i := 1 to self.Count do
    Result.pon_e(i,i, getItem(i-1).fguia_extendida_(getItem(i-1).factor, kPasoPronos(kPaso), 0));
end;

function TPronosticos.getLargoGuia: Integer;
var
  k: Integer;
begin
  Result:=0;
  for k:=0 to self.Count-1 do
    Result:=max(Result, Length(self[k].guia));
end;

function TPronosticos.kPasoPronos(kPaso: Integer): Integer;
var
  fechaIni: TFecha;
  dt: TDateTime;
begin
  if globs.EstadoDeLaSala=CES_OPTIMIZANDO then
    fechaIni:=globs.fechaIniOpt
  else
    fechaIni:=globs.fechaIniSim;

  dt:=fechaIni.dt+globs.DurPaso_minutos*kPaso/60/24;

  if dt<fechaIniCono.dt then
    Result:= 0
  else
    Result:=trunc( (dt-fechaIniCono.dt)*24*60/globs.DurPaso_minutos );

end;

function TPronosticos.getItem(i: integer): TPronostico;
begin
  Result := TPronostico(inherited items[i]);
end;

procedure TPronosticos.setItem(i: integer; cosa: TPronostico);
begin
  inherited items[i] := cosa;
end;

procedure TPronosticos.Free;
begin
  if fechaIniCono <> nil then fechaIniCono.Free;
  if Assigned(matricesRuida) then matricesRuida.Free;
  inherited Free;
end;

function TPronosticos.NSeries: integer;
begin
  result:= self.Count;
end;

constructor TPronosticos.Create_Default(cegh: TModeloCEGH);
var
  i: integer;
begin
  inherited Create(0, '');
  self.datosModelo := cegh;
  for i := 0 to cegh.NombresDeBornes_Publicados.Count - 1 do
    self.Add(TPronostico.Create_Default(cegh));
  Probabilidad := 0.0;
  fechaIniCono := TFecha.Create_Dt(0);
end;

procedure TPronosticos.prepararse(DuracionPasoSorteoEnHoras: NReal);
var
  k: integer;
begin
  for k := 0 to Count - 1 do
  begin
    items[k].calcular_rangos(DuracionPasoSorteoEnHoras);
    items[k].escenario:= self;
  end;
end;

function TPronosticos.max_NPasosCono: integer;
var
  k: integer;
  m, res: integer;
begin
  res := TPronostico(items[0]).NPasosCono;
  for k := 1 to Count - 1 do
  begin
    m := TPronostico(items[k]).NPasosCono;
    if (m > res) then
      res := m;
  end;
  Result := res;
end;

constructor TPronosticos.Create_ReadFromText(f: TArchiTexto);
begin

  if f.Version >= 143 then
  begin
    if f.Version < 164 then
    begin
      f.rd('P', Probabilidad);
      Probabilidad := Probabilidad / 100.0;
    end
    else
      f.rd('Probabilidad', Probabilidad);
  end
  else
    Probabilidad := 1.0;


  if f.Version >= 166 then
    f.rd( 'fechaIniCono', fechaIniCono )
  else
    fechaIniCono:= TFecha.Create_Dt(0);

  inherited Create_ReadFromText(f);
end;

procedure TPronosticos.WriteToText_(f: TArchiTexto);
var
  P100: integer;

begin
  if uversion_architexto.VERSION_ArchiTexto < 164 then
  begin
    P100 := Round(Probabilidad * 100);
    f.wr('P', P100);
  end
  else
    f.wr('Probabilidad', Probabilidad);

  f.wr( 'fechaIniCono', fechaIniCono );

  inherited WriteToText_(f);
end;


constructor TPronostico.Create(capa: integer; xNPCC, xNPLC, xNPSA, xNPAC: integer;
  xurl_get, xnombre_get: string);
begin
  inherited Create(capa);

  NPCC := xNPCC;
  NPLC := xNPLC;
  NPSA := xNPSA;
  NPAC := xNPAC;

  url_get := xurl_get;
  nombre_get := xnombre_get;

  rangoFechaSesgo := 0;
  guia_pe := 0.5;
  self.AperturaInicial_pu := 0.0;
  //  cantValoresDeterministicosUsados := 0;
end;


procedure TPronostico.Cambiar_GUIA(xNPCC, xNPLC, xNPSA, xNPAC: integer;
  xNuevaGuia: TDAOfNReal);
begin
  NPCC := xNPCC;
  NPLC := xNPLC;
  NPSA := xNPSA;
  NPAC := xNPAC;

  guia := copy(xNuevaGuia);

  //fechaIniSesgo_ := TFecha.Create_Dt(0);
  rangoFechaSesgo := 0;
  guia_pe := 0.5;
end;


constructor TPronostico.Create_Default(cegh: TModeloCEGH);
begin
  Create(0, 0, 0, 0, 0, '', '');
  SetLength(guia, cegh.CalcOrdenDelFiltro);
  vclear(guia);

end;

function TPronostico.Rec: TCosa_RecLnk;
begin
  Result := inherited Rec;
  Result.addCampoDef('Serie', Serie_, 133, 143);
  Result.addCampoDef('p', px.p, 133, 143);
  // Si la version es menor el nombre de los bornes para los pronosticos los carga
  // Create_ReadFromText de la fuente
  Result.addCampoDef('nretardos', px.nretardos, 0, 143);
  Result.addCampoDef('guia', guia);
  Result.addCampoDef('guia_eg', guia_eg, 82, 143);
  Result.addCampoDef('guia_pe', guia_pe, 0, 143 );
  Result.addCampoDef('NPCC', NPCC);
  Result.addCampoDef('NPLC', NPLC);
  Result.addCampoDef('NPSA', NPSA);
  Result.addCampoDef('NPAC', NPAC);
  Result.addCampoDef('sesgo', sesgo);
  Result.addCampoDef('factor', factor);

  Result.addCampoDef('fechaIniSesgo', px._fechaIniSesgo_, 0, 143);
  Result.addCampoDef('rangoFechaSesgo', rangoFechaSesgo, 0, 143);

  Result.addCampoDef('fechaIniCono', fechaIniCono_dummy, 164, 166 );
  Result.addCampoDef('AperturaInicial_pu', AperturaInicial_pu, 164);

  Result.addCampoDef('url_get', url_get, 126);
  Result.addCampoDef('nombre_get', nombre_get, 126);
end;

procedure TPronostico.BeforeRead(version, id_hilo: integer);
begin
  inherited BeforeRead(version, id_hilo);
  url_get := '';
  nombre_get := '';
  px := T_aux_pronostico.Create;
end;

procedure TPronostico.AfterRead(f: TArchiTexto);
begin
  inherited AfterRead(f);
  if f.version < 82 then
    setlength(guia_eg, length(guia)); // no los dejo en CERO.
  //  cantValoresDeterministicosUsados := 0;


  px.Free;
end;

function TPronostico.NPasosCono: integer;
begin
  Result := max(NPCC + NPLC, NPSA + NPAC);
end;


procedure TPronostico.calcular_rangos(DuracionPasoSorteoEnHoras: NReal);
begin
  rangoFechaSesgo := (NPCC + NPLC) * DuracionPasoSorteoEnHoras / 24.0;
  rangoFechaFactor := (NPSA + NPAC) * DuracionPasoSorteoEnHoras / 24.0;
end;

function TPronostico.fechaEnRango(
  fechaDeInicioDelPaso: TFecha;
  var kSesgo, kFactor: integer): boolean;
var
  dxdt: NReal;
begin
  dxdt := fechaDeInicioDelPaso.dt - escenario.fechaIniCono.dt + C_DT_1MINUTO;
  if (dxdt > 0) then
  begin
    if (dxdt < rangoFechaSesgo) then
      kSesgo := trunc(dxdt / rangoFechaSesgo * length(sesgo))
    else
      kSesgo := -1; // indicamos fuera de rango.

    if (dxdt < rangoFechaFactor) then
      kFactor := trunc(dxdt / rangoFechaFactor * length(factor))
    else
      kFactor := -1; // indicamos fuera de rango.
  end
  else
  begin
    kSesgo := -1;
    kFactor := -1;
  end;
  Result := (kSesgo >= 0) or (kFactor >= 0);
end;

function TPronostico.fguia_extendida_(aGuia: TDAOfNReal; kPaso: integer;
  DefVal: NReal): NReal;
begin
  if kPaso < length(aGuia) then
    Result := aGuia[kPaso]
  else
    Result := DefVal;
end;

procedure TPronostico.Free;
begin
  setlength(guia, 0);
  setlength(sesgo, 0);
  setlength(factor, 0);
  inherited Free;
end;

function TPronostico.determinismoSoloEstadoInicial: boolean;
begin
  Result := NPSA = 0;
end;

procedure AlInicio;
begin
  registrarClaseDeCosa(TPronostico.ClassName, TPronostico);
  registrarClaseDeCosa(TPronosticos.ClassName, TPronosticos);
  registrarClaseDeCosa(TEscenarioDePronosticos.ClassName, TEscenarioDePronosticos);
end;

procedure AlFinal;
begin
end;

{ TMatricesRuida }

{$IFDEF DUODINAMICO}
function TMatricesRuida.Get_Bamp(kPaso: integer): TMatR;
begin
  if A.nf <> A.nc then
    Result := MatrizRuidaMultiRetardos
  else
  begin
    if lst.Count > kPaso then
      Result := lst.Last
    else
      Result := lst[kPaso];
  end;
end;

function TMatricesRuida.Count: integer;
begin
  Result := lst.Count;
end;
{$ENDIF}


constructor TMatricesRuida.Create(aModelo: TModeloCEGH;
  aEscenarioPronos: TPronosticos);
var
  Bt: TMatR;
begin

  modelo:=aModelo;
  escenarioPronos:=aEscenarioPronos;
  {$IFDEF DUODINAMICO}
  lst := TList.Create;
  {$ENDIF}

  self.A := modelo.A_cte;
  self.B := modelo.B_cte;

  if modelo.nVE=0 then
  begin
    self.Ma:=nil;
    self.Mr:=nil;
  end
  else
  begin
    self.Ma := modelo.MAmp_cte;
    self.Mr := modelo.MRed_cte;
  end;


  At := TMatR.Create_Clone(A);
  At.Transponer;
  if A.nf <> A.nc then
  begin
  {$IFDEF RUIDA_MULTI_RETARDOS}
    MatrizRuidaMultiRetardos := Matriz_RUIDA_MultiRetardos(A, B, Ma, Mr);
  {$ELSE}
    MatrizRuidaMultiRetardos := nil;
  {$ENDIF}
  end;

  BampEE := Matriz_RUIDA_EE(A, B, Ma, Mr);

  if Mr = nil then
  begin
    MaMr := nil;
    MaMr_t := nil;
  end
  else
  begin
    MaMr := TMatR.Create_Init(Ma.nf, Mr.nc);
    MaMr.Mult(Ma, Mr);
    MaMr_t := TMatR.Create_Clone(MaMr);
    MaMr_t.Transponer;
  end;

  Bt := TMatR.Create_Clone(B);
  Bt.Transponer;
  BBt := TMatR.Create_Init(B.nf, B.nf);
  BBt.Mult(B, Bt);
  Bt.Free;

end;

procedure TMatricesRuida.calcularse(globs: TGlobs);
var
  kIter, dimRes: Integer;
  // Estado actual y estado siguiente
  XXt, XsXst: TMatR;

  XXtAt, AXXtAt, Lambda_k, Lambda_k_t, BBtLambda_k_t,
    Lambda_kBBtLambda_k_t,
    AVE_XVE_X_tA_t, S_kVE_X_tA_t, CovX_s, Kr, Id_menos_Kr, Id,
    S_k_tVE_XA, Bamp_kBamp_k_t, Bamp_k, MrH, Fa, H: TMatR;

  VE_Xa, VE_X, AVE_Xa, AVE_X, VE_X_tA_t, S_k, vecCovX_s, vecBampBamp_t: TVectR;

begin

  //XXt inicialmente es una matriz nula porque el estado es deterministico
  XXt := TMatR.Create_InitVal(A.nf, A.nf,0.0);

  VE_Xa:=TVectR.Create_InitVal(A.nc);

  if Assigned(Mr) then
  begin
    Kr := TMatR.Create_Kron(MaMr, MaMr);
    Id:=TMatR.Create_identidad(Kr.nf);

    Id_menos_Kr:=TMatR.Create_Init(Kr.nf, Kr.nc);
    Kr.PorReal(-1);
    Id_menos_Kr.Suma(Id, Kr);
  end;

  kIter:=0;
  repeat

    Lambda_k :=escenarioPronos.GetLambda(kIter) ;
    Lambda_k_t:=Lambda_k.Crear_Transpuesta;

    S_k:=escenarioPronos.Sesgos(kIter);

    VE_X:=escenarioPronos.Guia_eg(kIter);
    if kIter>=escenarioPronos.getLargoGuia then
    begin
      AVE_Xa := TVectR.Create_Init(A.nf);
      AVE_Xa.Mult(A, VE_Xa);
      VE_X.suma(AVE_Xa, escenarioPronos.Sesgos(kIter));
      AVE_Xa.Free;
    end;

    // A * <X_k * X_k'> A'
    XXtAt := TMatR.Create_Init(XXt.nf, At.nc);
    XXtAt.Mult(XXt, At);
    AXXtAt := TMatR.Create_Init(A.nf, XXtAt.nc);
    AXXtAt.Mult(A, XXtAt);
    XXtAt.Free;

    // Lambda_k * B * B' * Lambda_k'
    BBtLambda_k_t := TMatR.Create_Init(BBt.nf, Lambda_k_t.nc);
    BBtLambda_k_t.Mult(BBt, Lambda_k_t);
    Lambda_kBBtLambda_k_t := TMatR.Create_Init(lambda_k.nf, BBtLambda_k_t.nc);
    Lambda_kBBtLambda_k_t.Mult(lambda_k, BBtLambda_k_t);
    BBtLambda_k_t.Free;

    // A * <X>
    AVE_X:=TVectR.Create_Init(A.nf);
    AVE_X.Mult(A, VE_X);

    // <X'> * A'
    VE_X_tA_t:=AVE_X.clonar;

    // A * <X> * <X'> * A'
    AVE_XVE_X_tA_t := TMatR.Create_Init_pm(AVE_X.n, VE_X_tA_t.n) ;
    AVE_XVE_X_tA_t.Mult(AVE_X, VE_X_tA_t);

    // S_k * <X'> * A'
    S_kVE_X_tA_t:=TMatR.Create_Init_pm(S_k.n, VE_X_tA_t.n);
    S_kVE_X_tA_t.Mult(S_k, VE_X_tA_t);
    VE_X_tA_t.Free;

    // S_k_t * <X> * A
    S_k_tVE_XA:=S_kVE_X_tA_t.Crear_Transpuesta;
    AVE_X.Free;

    CovX_s:=TMatR.Create_InitVal(AXXtAt.nf, AXXtAt.nc, 0.0);
    CovX_s.Suma(AXXtAt, Lambda_kBBtLambda_k_t);
    CovX_s.Suma(CovX_s, AVE_XVE_X_tA_t);
    CovX_s.Suma(CovX_s, S_kVE_X_tA_t);
    CovX_s.Suma(CovX_s, S_k_tVE_XA);

    AXXtAt.Free;
    Lambda_kBBtLambda_k_t.Free;
    AVE_XVE_X_tA_t.Free;
    S_kVE_X_tA_t.Free;
    S_k_tVE_XA.Free;


    if Assigned(Mr) then
    begin
      vecCovX_s:=CovX_s.vec;
      vecBampBamp_t:=TVectR.Create_Init(vecCovX_s.n);
      vecBampBamp_t.Mult(Id_menos_Kr, vecCovX_s);

      Bamp_kBamp_k_t:=vecBampBamp_t.reshape(Ma.nf, Ma.nc);

      CovX_s.Free;
    end
    else
      Bamp_kBamp_k_t:=CovX_s;

    H := Bamp_kBamp_k_t.RaizPorPotenciaIterada(dimRes, False);

    Bamp_kBamp_k_t.Free;

    if (H <> nil) and (dimRes > 0) then
    begin

      MrH := TMatR.Create_Init(Mr.nf, H.nc);
      MrH.Mult(Mr, H);

      MrH.CalcBSEN(Fa);

      Fa.OrtoNormal;
      Fa.Transponer;

      Bamp_k := TMatR.Create_Init(H.nf, Fa.nc);
      Bamp_k.Mult(H, Fa);

      H.Free;
      MrH.Free;
      Fa.Free;

    end
    else
    begin
      Bamp_k := nil;
    end;

 {$IFDEF DUODINAMICO}
    lst.Add(Bamp_k);
 {$ENDIF}
    VE_Xa.Free;
    VE_Xa:=VE_X;

    inc(kIter);

  until (max(globs.calcNPasosSim,globs.calcNPasosOpt)=kIter);


end;


end.
{$ENDIF}
