unit ufraccionadorcegh;

{$AppType CONSOLE}

interface

uses
  Classes, SysUtils,
  umodelosintcegh,
  matreal, xMatDefs, fddp_conmatr, fddp;


// Anterior desarrollo del Fraccionador que toma un archivo CEGH y un factor
// de Fraccionacion y calcula las Matrices A y B del filtro asi como las
// funciones deformantes del nuevo sintetizador (ver si se borra el botón)
procedure Viejo_Fraccionar(paso_discr: integer; archOrigen, archDestino: string);

// Nuevo Programa principal del Fraccionador que toma un archivo CEGH y un factor
// de Fraccionacion y calcula las Matrices A y B del filtro asi como las
// funciones deformantes del nuevo sintetizador
procedure Fraccionar(paso_discr: integer; arch_org, arch_dest: string);

// Funcion para interpolar las funciones deformantes cuando se Fracciona
function Interpolacion_funciones_deformantes(alf, Inferior, Superior: NReal): NReal;

// Procedimiento auxiliar para mostrar los datos del sintetizador origen
procedure muestra(arch: string; var a, b, c: integer);

//Cálculo de H = 1/N^2 sum( i = 1; i = N-1 ; (N-i) A^i )
//ec.14)
function CalcMatriz_H_(MA: TMatR; paso_discr: integer): TMatR;

// Si X_{k+1} = A X_k + B R_k es un proceso y consideramos el Promedio Movil
// de N Muestras, Y = PromedioMovil( X, N ) la matriz de covarianzas de Y
// es calculada como vec(Syy) = M vec(Sxx)
// La función retonra la matriz M
function MatTranfSigmaPromedioMovil(MA: TMatR; paso_discr: integer): TMatR;

// Devuelve la matriz diagnoal 1/sqrt(SigmaXX(j,j))
function CalcLamda(SigmaXX: TMatR): TMatR;

implementation



(*
ec.14)
Cálculo de H = 1/N^2 sum( i = 1; i = N-1 ; (N-i) A^i )
*)
function CalcMatriz_H_(MA: TMatR; paso_discr: integer): TMatR;
var
  H, aux: TMatR;
  i: integer;
  aux2: TMatR;

begin
  H := TMatR.Create_Init(MA.nf, MA.nc);
  H.Ceros;

  aux := TMatR.Create_Clone(MA);
  aux2 := TMatR.Create_CLone(MA);
  for i := 1 to paso_discr - 1 do
  begin
    aux2.Igual(aux);
    aux2.PorReal(paso_discr - i);
    H.Suma(H, aux2);
    if i < (paso_discr - 1) then
      aux.Mult(aux, MA);
  end;
  aux.Free;
  aux2.Free;

  H.PorReal(1 / (paso_discr * paso_discr));
  Result := H;
end;

// Matriz corchetes rectos ec.19
function MatTranfSigmaPromedioMovil(MA: TMatR; paso_discr: integer): TMatR;
var
  H, Mcron, aux1, aux2: TMatR;
  nn: integer;
  I: TMatR;
  k, j: integer;

begin
  H := CalcMatriz_H_(MA, paso_discr);

  I := TMatR.Create_Init(H.nf, H.nc);
  I.identidad;

  nn := H.nf * H.nc;
  Mcron := TMatR.Create_Init(nn, nn);
  Mcron.Ceros;

  aux1 := TMatR.Create_Kron(I, H);
  aux2 := TMatR.Create_Kron(H, I);

  for k := 1 to nn do
  begin
    for j := 1 to nn do
      if j = k then
        Mcron.pon_e(k, j, 1 / paso_discr + (aux1.e(k, j) + aux2.e(k, j)))
      else
        Mcron.pon_e(k, j, aux1.e(k, j) + aux2.e(k, j));
  end;

  aux1.Free;
  aux2.Free;
  I.Free;
  Result := Mcron;

end;

function Interpolacion_funciones_deformantes(alf, Inferior, Superior: NReal): NReal;
begin
  Result := alf * Inferior + (1 - alf) * Superior;
end;

function TransformarB(MA, MB: TMatR; n: integer): TMatR;
var
  M, Aj, BBT: TMatR;
  VecBBT: TvectR;
  mVecBBT: TMatR;

  i: integer;
  aux: TMatR;
  invertible: boolean;
  exp10: integer;
  res: TMatR;

begin
  Aj := TMatR.Create_Identidad(MA.nf);
  M := TMatR.Create_kron(Aj, Aj);
  for i := 1 to n - 1 do
  begin
    Aj.Mult(Aj, MA);
    aux := TMatR.Create_Kron(Aj, Aj);
    M.suma(M, aux);
    aux.Free;
  end;
  Aj.Free;

  BBT := MB.Crear_Transpuesta;
  BBT.Mult(MB, aux);
  VecBBT := BBT.vec;
  mVecBBT := VecBBT.reshape(VecBBT.n, 1);

  M.Escaler(mVecBBT, invertible, exp10);
  if not invertible then
    raise Exception.Create('ERror, no pude invertir sum( kron( aj^j, aj^j; j= 0..n-1 )'
      );
  mVecBBT.CopyColVect(VecBBT, 1);
  mVecBBT.Free;
  M.Free;

  BBT.Free;
  BBT := VecBBT.reshape(MB.nf, MB.nc);
  res := BBT.raiz_Cholesky;
  BBT.Free;
  VecBBT.Free;
  Result := res;

end;

// Teniendo en cuenta la relación Vec(DDT)= I x Vec(SigmaYY) - Kron(C,C)Vec(SigmaYY)
// En esta función se calcula Vec(SigmaYY) = Vec(DDT) x inv [I - Kron(C,C)]
// El resultado es el vector Vec (SigmaYY)

procedure Calcular_SigmaXXyB(MA, MC, MD: TMatR; paso_discr: integer;
  var SigmaXX, MatB: TMatR);
var
  Cj, AuxMC, DDT, BBT, MAT: TMatR;
  VecDDT: TVectR;

  invertible: boolean;
  exp10: integer;
  vaux: TVectR;

  SigmaYY: TMatR;
  M: TMatR;
  dimres: integer;

begin

  // Cálculo de SigmaYY ec.6
  // Vec(DDT)= I x Vect(SigmaYY) - Kron(C,C)Vect(SigmaYY)
  AuxMC := TMatR.Create_kron(MC, MC);
  Cj := TMatR.Create_Identidad(AuxMC.nf);
  AuxMC.PorReal(-1);
  AuxMC.Suma(Cj, AuxMC);
  Cj.Free;

  DDT := MD.Crear_Transpuesta;
  DDT.Mult(MD, DDT);
  VecDDT := DDT.vec;
  SigmaYY := VecDDT.reshape(VecDDT.n, 1);
  VecDDT.Free;
  DDT.Free;

  // Hallo Inv(I - kron(C,C)) x vec(DDT), lo que es igual a vec(SigmaYY)

  AuxMC.Escaler(SigmaYY, invertible, exp10);
  if not invertible then
    raise Exception.Create('Error, no pude invertir la función (I - kron(C,C))');
  // fin cálculo ec.6
  AuxMC.Free;


  // Ahora que tenemos SigmaYY, calculamos SigmaXX usando ec.19
  M := MatTranfSigmaPromedioMovil(MA, paso_discr);
  M.Escaler(SigmaYY, invertible, exp10);
  M.Free;

  vaux := SigmaYY.vec;
  SigmaXX := vaux.reshape(MA.nf, MA.nc);
  vaux.Free;
  SigmaYY.Free;

  BBT := TMatR.Create_Init(MA.nf, MA.nc);

  MAT := MA.Crear_Transpuesta;
  BBT.Mult(MA, SigmaXX);
  BBT.Mult(BBT, MAT);
  BBT.PorReal(-1);
  BBT.Suma(SigmaXX, BBT);
  MAT.Free;

  MatB := BBT.raiz_Cholesky;

  if MatB = nil then
  begin
    MatB := BBt.RaizPorPotenciaIterada(dimRes, False);
    if dimRes < 0 then
    begin
      BBt.WriteArchiXLT('bbt_autovalneg_.xlt');
      writeln('Ojo ... imposiblre raíz de bbt ');
    end;
  end;
  BBt.Free;

end;


function CalcLamda(SigmaXX: TMatR): TMatR;
var
  MatLamda: TMatR;
  aux: NReal;
  i: integer;

begin
  MatLamda := TMatR.Create_Init(SigmaXX.nf, SigmaXX.nc);
  MatLamda.Ceros;
  for i := 1 to SigmaXX.nf do
  begin
    aux := SigmaXX.e(i, i);
    aux := 1 / (sqrt(aux));
    MatLamda.pon_e(i, i, aux);
  end;
  Result := MatLamda;
end;

procedure Viejo_Fraccionar(paso_discr: integer; archOrigen, archDestino: string
  );
var
  sintetizador_origen, sintetizador_destino: TModeloCEGH;
  MatA, MatB: TMatR;
  NSS_origen, NFD_origen, NPFD_origen, NFD_destino, kserie, kpaso,
  k, j, DurPasoSorteo_origen, DurPasoSorteo_destino: integer;
  F_inf, F_sup, alfa: NReal;
  aux_Vect: TVectR;
  ir: NReal;
  funcionesDeformantesI, funcionesDeformantesOrg: TDAOf_ddp_VectDeMuestras;
  res: TMatR;

begin
  sintetizador_origen := TModeloCEGH.CreateFromArchi(archOrigen);
  sintetizador_destino := TModeloCEGH.CreateFromArchi(archOrigen);

  MatA := sintetizador_origen.A_cte.power_r(1.0 / paso_Discr);

  MatB := TMatR.Create_Init(sintetizador_origen.B_cte.nf, sintetizador_origen.B_cte.nc);
  MatB.Igual(sintetizador_origen.B_cte);

  sintetizador_destino.A_cte.Igual(MatA);

  res := TransformarB(MatA, MatB, paso_discr);
  if res <> nil then
  begin
    sintetizador_destino.B_cte.Igual(res);
    res.Free;
  end
  else
    raise Exception.Create('Fraccionar ... FALLO CHOLESKY al calcular B ');

  // Interpolacion de los datos entre muestras de mayor frecuencia

  NSS_origen := sintetizador_origen.nBornesSalida;

  // Cantidad de dformadores por cada serie (lo miro sobre la primer serie)
  NFD_origen := length(sintetizador_origen.funcionesDeformantes[0]);

  // Cantidad de puntos de discretización de un deformador
  NPFD_origen := sintetizador_origen.funcionesDeformantes[0][0].a.n;

  DurPasoSorteo_origen := sintetizador_origen.durPasoDeSorteoEnHoras;

  NFD_destino := NFD_origen * paso_discr;
  DurPasoSorteo_destino := trunc(DurPasoSorteo_origen / paso_discr);
  sintetizador_destino.durPasoDeSorteoEnHoras := DurPasoSorteo_destino;

  for kserie := 0 to NSS_origen - 1 do
  begin

    SetLength(sintetizador_destino.funcionesDeformantes[kserie], NFD_destino);

    funcionesDeformantesI := sintetizador_destino.funcionesDeformantes[kserie];

    funcionesDeformantesOrg := sintetizador_Origen.funcionesDeformantes[kserie];
    for kpaso := 0 to NFD_Destino - 1 do
    begin
      ir := kPaso / paso_Discr;
      k := trunc(ir);
      alfa := 1 - frac(ir);

      aux_Vect := TVectR.Create_Init(NPFD_origen);
      for j := 0 to NPFD_origen do
      begin
        F_inf := funcionesDeformantesOrg[k].a.e(j);

        if (k+1) > high( funcionesDeformantesOrg ) then
          F_sup := funcionesDeformantesOrg[0].a.e(j)
        else
          F_sup := funcionesDeformantesOrg[k + 1].a.e(j);

        aux_Vect.pon_e(j, F_inf * alfa + F_sup * (1 - alfa));

      end;

      funcionesDeformantesI[kpaso] :=
        Tf_ddp_VectDeMuestras.Create_SinClonarMuestras(aux_Vect, nil, 0);
    end;
  end;

  sintetizador_destino.WriteToArchi(archDestino, 12, 4);
  sintetizador_origen.NombresDeBornes_Publicados.Clear;
  sintetizador_destino.NombresDeBornes_Publicados.Clear;
  sintetizador_origen.Free;
end;

procedure Calc_A_B_DiagInvLambda(var A_s, B_s: TMatR; var DiagInvLambda: TVectR;
  A, B: TMatR; paso_Discr: integer);

var
  MatA, MatB, MatC, MatD, InvLambda: TMatR;
  SigmaXX, Lambda: TMatR;
  i: integer;
  inv_aux: NReal;

begin
  MatC := A;
  MatA := MatC.power_r(1.0 / paso_Discr);  //ec.15
  MatD := TMatR.Create_Clone(B);

  {$IFDEF dbg}
  MatA.WriteArchiXLT('matA.xlt');
  MatC.WriteArchiXLT('matC.xlt');
  MatD.WriteArchiXLT('matD.xlt');
  MatB.WriteArchiXLT('matB.xlt');
  SigmaXX.WriteArchiXLT('sigmaXX.xlt');
  {$ENDIF}

  Calcular_SigmaXXyB(MatA, MatC, MatD, paso_discr, SigmaXX, MatB); // ec.16

  {$IFDEF dbg}
  MatA.WriteArchiXLT('matA.xlt');
  MatC.WriteArchiXLT('matC.xlt');
  MatD.WriteArchiXLT('matD.xlt');
  MatB.WriteArchiXLT('matB.xlt');
  SigmaXX.WriteArchiXLT('sigmaXX.xlt');
  {$ENDIF}

  // Hallo la matriz A* y B* tal que A* = Lamda x A x InvLamda y B* = Lamda x B
  // La matrices del filtro CEGH van a ser A* y B*

  Lambda := CalcLamda(SigmaXX);
  InvLambda := TMatR.Create_identidad(Lambda.nc);
  DiagInvLambda := TVectR.Create_init(lambda.nc);
  for i := 1 to Lambda.nf do
  begin
    inv_aux := 1 / Lambda.e(i, i);
    DiagInvLambda.pon_e(i, inv_aux);
    InvLambda.pon_e(i, i, inv_aux);
  end;

  MatA.Mult(Lambda, MatA);
  MatA.Mult(MatA, InvLambda);
  A_s.Igual(MatA);
  MatA.Free;

  if MatB <> nil then
  begin
    MatB.Mult(Lambda, MatB);
    B_s.Igual(MatB);
    MatB.Free;
  end
  else
    raise Exception.Create('Fraccionar ... FALLO CHOLESKY al calcular B ');

  Lambda.Free;
  InvLambda.Free;
end;

procedure Calc_A_B_C_DiagInvLambda_(var A_s, B_s, C_s: TMatR; var DiagInvLambda: TVectR;
  A, B, C: TMatR; paso_Discr: integer);

var
  MatA, MatB, MatC, MatD, InvLambda: TMatR;
  SigmaXX, Lambda: TMatR;
  i: integer;
  inv_aux: NReal;
  AC: TMatR;
  nfilasAC: integer;
begin
  //para el modelo del CEGH con entradas aparte (C != 0) uno la matriz A con la C
  //para calcular lambda, despues las separo
  nfilasAC := A.nf + C.nc;
  if C <> Nil then
  begin
    AC := TMatR.Create_Init(nfilasAC, A.nc);
    for i := 1 to A.nf do
    begin
      AC.pon_fila(1, 1, A.Fila(i).toTDAOfNReal());
    end;
    for i := 1 to C.nc do
    begin
      AC.pon_fila(1, 1, C.Columna(i).toTDAOfNReal());
    end;
  end
  else
    AC := A;

  MatC := AC;
  MatA := MatC.power_r(1.0 / paso_Discr);  //ec.15
  MatD := TMatR.Create_Clone(B);


  Calcular_SigmaXXyB(MatA, MatC, MatD, paso_discr, SigmaXX, MatB); // ec.16

  // Hallo la matriz A* y B* tal que A* = Lamda x A x InvLamda y B* = Lamda x B
  // La matrices del filtro CEGH van a ser A* y B*

  Lambda := CalcLamda(SigmaXX);
  InvLambda := TMatR.Create_identidad(Lambda.nc);
  DiagInvLambda := TVectR.Create_init(lambda.nc);
  for i := 1 to Lambda.nf do
  begin
    inv_aux := 1 / Lambda.e(i, i);
    DiagInvLambda.pon_e(i, inv_aux);
    InvLambda.pon_e(i, i, inv_aux);
  end;

  MatA.Mult(Lambda, MatA);
  MatA.Mult(MatA, InvLambda);
  A_s.Igual(MatA);
  MatA.Free;

  if MatB <> nil then
  begin
    MatB.Mult(Lambda, MatB);
    B_s.Igual(MatB);
    MatB.Free;
  end
  else
    raise Exception.Create('Fraccionar ... FALLO CHOLESKY al calcular B ');

  Lambda.Free;
  InvLambda.Free;
end;


procedure Fraccionar(paso_discr: integer; arch_org, arch_dest: string);
var
  sintetizador_origen, sintetizador_destino: TModeloCEGH;
  NSS_origen, NFD_origen, NPFD_origen, NFD_destino, kserie, kpaso,
  k_inf, k_sup, j, resultado, DurPasoSorteo_origen, DurPasoSorteo_destino: integer;
  F_inf, F_sup, alfa: NReal;
  aux_Vect, nuevosPesos: TVectR;
  ir, my, stdy, ratio: NReal;
  funcionesDeformantesDes, funcionesDeformantesOrg: TDAOf_ddp_VectDeMuestras;
  cnt_inters: integer;
  nuevasmuestras: TVectR;
  auxInvLambda: NReal;

  aDeformador: Tf_ddp_VectDeMuestras;

  DiagInvLambda: TVectR;
  DiagsInvLambda: array of TVectR;
  flg_MonoFiltro: boolean;

  mca_Destino, mcb_Destino: TDAOfMatR;
begin
  arch_dest := arch_org + '111.txt';


  sintetizador_origen := TModeloCEGH.CreateFromArchi(arch_org);
  sintetizador_destino := TModeloCEGH.CreateFromArchi(arch_org);


  // Interpolacion de los datos entre muestras de mayor frecuencia
  NSS_origen := sintetizador_origen.nBornesSalida;

  // Cantidad de deformadores por cada serie (lo miro sobre la primer serie)
  NFD_origen := length(sintetizador_origen.funcionesDeformantes[0]);

  // Cantidad de puntos de discretización de un deformador
  NPFD_origen := sintetizador_origen.funcionesDeformantes[0][0].a.n;

  {$IFDEF dbg}
  sintetizador_origen.WriteToArchi('prueba.xlt', NFD_origen, NPFD_origen);
  {$ENDIF}

  DurPasoSorteo_origen := sintetizador_origen.durPasoDeSorteoEnHoras;

  DurPasoSorteo_destino := trunc(DurPasoSorteo_origen / paso_discr + 0.01);
  sintetizador_destino.durPasoDeSorteoEnHoras := DurPasoSorteo_destino;

  NFD_destino := NFD_origen * paso_discr;

  flg_MonoFiltro := sintetizador_destino.A_cte <> nil;
  DiagInvLambda := nil;

  if flg_MonoFiltro then
  begin
    {$IFDEF dbg_VC}
    Calc_A_B_C_DiagInvLambda_(
      sintetizador_destino.A_cte, sintetizador_destino.B_cte, sintetizador_destino.C_cte,
      DiagInvLambda,
      sintetizador_origen.A_cte, sintetizador_origen.B_cte, sintetizador_origen.C_cte, paso_discr);
    {$ELSE}
    Calc_A_B_DiagInvLambda(
      sintetizador_destino.A_cte, sintetizador_destino.B_cte,
      DiagInvLambda,
      sintetizador_origen.A_cte, sintetizador_origen.B_cte, paso_discr);
    {$ENDIF}

  end
  else
  begin
    setlength(DiagsInvLambda, NFD_Destino);
    for kPaso := 0 to NFD_Origen - 1 do
    begin
      writeln('kPaso: ', kPaso);
      DiagsInvLambda[kPaso] := nil;
      Calc_A_B_DiagInvLambda(
        sintetizador_destino.mcA[kPaso], sintetizador_destino.mcB[kPaso],
        DiagsInvLambda[kPaso],
        sintetizador_origen.mcA[kPaso], sintetizador_origen.mcB[kPaso], paso_discr);
    end;

    setlength(mcA_Destino, NFD_Destino);
    setlength(mcB_Destino, NFD_Destino);
    for kpaso := 0 to NFD_Destino - 1 do
    begin
      ir := kPaso / paso_Discr;
      k_inf := trunc(ir) mod NFD_Origen;
      k_sup := (k_inf + 1) mod NFD_Origen;
      alfa := 1 - frac(ir);
      mcA_Destino[kPaso] := TMatR.Create_Combinar(
        alfa, sintetizador_destino.mcA[k_inf], (1 - alfa),
        sintetizador_destino.mcA[k_sup]);
      mcB_Destino[kPaso] := TMatR.Create_Combinar(
        alfa, sintetizador_destino.mcB[k_inf], (1 - alfa),
        sintetizador_destino.mcB[k_sup]);
    end;
    for kPaso := 0 to NFD_Origen - 1 do
    begin
      sintetizador_destino.mcA[kPaso].Free;
      sintetizador_destino.mcB[kPaso].Free;
    end;
    setlength(sintetizador_destino.mcA, 0);
    setlength(sintetizador_destino.mcB, 0);
    sintetizador_destino.mcA := mcA_Destino;
    sintetizador_destino.mcB := mcB_Destino;
  end;


  for kserie := 0 to NSS_origen - 1 do
  begin
    SetLength(sintetizador_destino.funcionesDeformantes[kserie], NFD_destino);
    funcionesDeformantesDes := sintetizador_destino.funcionesDeformantes[kserie];
    funcionesDeformantesOrg := sintetizador_Origen.funcionesDeformantes[kserie];

    for kpaso := 0 to NFD_Destino - 1 do
    begin
      ir := kPaso / paso_Discr;
      k_inf := trunc(ir) mod NFD_Origen;
      k_sup := (k_inf + 1) mod NFD_Origen;
      alfa := 1 - frac(ir);
      aux_Vect := TVectR.Create_Init(NPFD_origen);
      for j := 1 to NPFD_origen do
      begin
        F_inf := funcionesDeformantesOrg[k_inf].a.e(j);
        F_sup := funcionesDeformantesOrg[k_sup].a.e(j);
        aux_Vect.pon_e(j, F_inf * alfa + F_sup * (1 - alfa));
      end;
      funcionesDeformantesDes[kpaso] :=
        Tf_ddp_VectDeMuestras.Create_SinClonarMuestras(aux_Vect, nil, 0);
    end;

    for kPaso := 0 to NFD_Destino - 1 do
    begin
      aDeformador := funcionesDeformantesDes[kpaso];
      aDeformador.a.promedioDesvEst(my, stdy, False);
      if flg_MonoFiltro then
        auxInvLambda := DiagInvLambda.e(kSerie + 1)
      else
      begin
        ir := kPaso / paso_Discr;
        k_inf := trunc(ir) mod NFD_Origen;
        k_sup := (k_inf + 1) mod NFD_Origen;
        alfa := 1 - frac(ir);
        auxInvLambda := alfa * DiagsInvLambda[k_inf].e(kSerie + 1) +
          (1 - alfa) * DiagsInvLambda[k_sup].e(kSerie + 1);
      end;

      aDeformador.a.PorReal(auxInvLambda);
      stdy := stdy * auxInvLambda;

      resultado := cambiarHistograma(aDeformador.a, nil, my, stdy,
        nuevosPesos, ratio, cnt_inters);

      if resultado <> 0 then
        raise Exception.Create('No converge CambiarHistograma');
      nuevasMuestras := muestrasEquiprobables(aDeformador.a,
        nuevosPesos, aDeformador.a.n);
      nuevosPesos.Free;
      aDeformador.a.Free;
      aDeformador.a := nuevasMuestras;

    end;
  end;

  if flg_MonoFiltro then
    DiagInvLambda.Free
  else
  begin
    for kPaso := 0 to NFD_Destino - 1 do
      DiagsInvLambda[kPaso].Free;
    setlength(DiagsInvLambda, 0);
  end;

  sintetizador_destino.WriteToArchi(arch_dest, 0, 3);
  sintetizador_origen.NombresDeBornes_Publicados.Clear;
  sintetizador_destino.NombresDeBornes_Publicados.Clear;
  sintetizador_origen.Free;
end;


procedure muestra(arch: string; var a, b, c: integer);
var
  sintetizador: TModeloCEGH;
begin
  sintetizador := TModeloCEGH.CreateFromArchi(arch);
  a := sintetizador.nBornesSalida;
  b := sintetizador.nPuntosPorPeriodo;
  c := sintetizador.funcionesDeformantes[0][0].a.n;
  sintetizador.Free;
end;

end.
