unit uparsimoniat;

{$mode delphi}

interface

uses
  Classes, SysUtils, math, xmatdefs, matreal;

(**
 Una Parsimonia, administra un juego de parámetros PR en base a un juego
 de parámetros de parsiomonia PP de menor dimensión.
 Para ello, cada componente de PR se expresa como una constante
 más una serie de potencias del tiempo t (de inicio del paso) y
 una serie de Fourier (cosenos y senos) de una frecuencia fundamentamental y
 sus armónicos.

 Dim(PP) = 1 + nPotencias + 2* nArmonicos

 La idea es que se dispone de un tipo de modelo CF( X, t ) = M( X, PR(t) )
 al que se le quiere dar Parsimonia imponiendo PR(t) = fp( PP, t )

 En la calibramción de M se calcula el gradiente del e2 respecto de PR
 e2 = sum( i, t; ( M(X_i, t, PR(t)) - Y_i  )^2 ) ; suma en i y en t ( t = k*T)

 Para poder optimizar en base a la Parsimonia, hay que aplicar la regla de la
 cadena y derivar para cada paso t y cada componente de PR respecto de las componentes
 de PP

 de2/dPP_h = sum( i, t; 2 (  M(X_i, t, PR(t) )-Y_i ) dM(X_i, t, PR(t) )/dPR ) dPR(t)/dPP_h

 sum( i, t; 2 ( Y_i - M(X_i, t, PR(t) ) dM(X_i, t, PR(t) )/dPR) es el gradiente del error cuadrático
 calculado sobre los parámetros del modelo sin parsimonia.

 dPR(t)/dPP_h

 son las derivadas de los parámetros del modelo sin parsimonia en el paso t
 respecto del parámetro h de la Parsimonia.

 La Parsimonia debe saber calcular PR = fp( t, PP )
 y también dPR(t)/dPP = gPR = gp( t, PP )
 por eficiencia del cálculo gPR es una matriz, donde la fila h es
 el transpuesto de dPR(t)/dPP_h

 así

 de2/dPP_h =  sum( i, t; gPRPP^T  * 2 ( Y_i - M(X_i, t, PR(t) ) dM(X_i, t, PR(t) )/dPR )

 de2/dPP_h =  sum( t; gPRPP^T  *sum( i; 2 ( Y_i - M(X_i, t, PR(t) ) dM(X_i, t, PR(t) )/dPR ) )


**)

type

  { TParsimonia }

  TParsimonia = class
    dtPasoT, dtCicloPrincipal: NReal;
    nPotencias, nArmonicos, nParametrosPorPaso: integer;

    vpot: TVectR;
    vcos: TVectR;
    vsin: TVectR;

    dimParsimonia: integer; // ( 1 + nPotencias + 2 * nArmonicos )

    constructor Create(
      dtPasoT, dtCicloPrincipal: NReal; nPotencias, nArmonicos,
      nParametrosPorPaso: integer);

    function dim_xp: NReal; // ( 1 + nPotencias + 2* nArmonicos ) * nParametros

    // actualiza los vectores pot, cos y sin
    // de acuerdo al tiempo
    // kPaso = 1 ... nPasosT
    procedure set_dtk( kPaso: integer );

    // calcula vector de parametros parsimoniosos según el vector de parámetros
    procedure fp(var PR: TVectR; const x: TVectR);

    // calcula los gadientes de d_pr/ dxp
    procedure gp( var gPR: TMatR );


    destructor destroy; override;

  end;



implementation

{ TParsimonia }

constructor TParsimonia.Create(dtPasoT, dtCicloPrincipal: NReal; nPotencias,
  nArmonicos, nParametrosPorPaso: integer);
begin
  self.dtPasoT:= dtPasoT;
  self.dtCicloPrincipal:= dtCicloPrincipal;
  self.nPotencias:= nPotencias;
  self.nArmonicos:= nArmonicos;
  self.nParametrosPorPaso:= nParametrosPorPaso;

  dimParsimonia:= 1 + nPotencias + 2 * nArmonicos;

  vpot:= TVectR.Create_Init( nPotencias );
  vcos:= TVectR.Create_Init( nArmonicos );
  vsin:= TVectR.Create_init( nArmonicos );
end;

function TParsimonia.dim_xp: NReal;
begin
  result:= ( 1 + nPotencias + 2 * nArmonicos ) * nParametrosPorPaso;
end;

procedure TParsimonia.set_dtk( kPaso: integer);
var
  k: integer;
  p, w0, alfa: NReal;
  t: NReal;
begin
  t:= ( kPaso - 1 ) * dtPasoT;

  p:= 1;
  for k:= 1 to vpot.n do
  begin
    p:= p * t;
    vpot.pon_e( k, p );
  end;

  w0:= 2 * pi * t / dtCicloPrincipal;
  for k:= 1 to vcos.n do
  begin
    alfa:= w0 * k;
    vcos.pon_e( k, cos( alfa ) );
    vsin.pon_e( k, sin( alfa ) );
  end;

end;

procedure TParsimonia.fp( var PR: TVectR; const x: TVectR);
var
  kPar, kOffset, j: integer;
  a: NReal;
begin
  kOffset:= 1;
  for kPar:= 1 to PR.n do
  begin
    a:= X.e( kOffset );
    inc( kOffset );

    for j:= 1 to nPotencias do
    begin
     a:= a + vpot.e( j ) * x.e( kOffset );
     inc( kOffset );
    end;

    for j:= 1 to nArmonicos do
    begin
      a:= a + vcos.e(j) * x.e( kOffset ) + vsin.e(j) * x.e( kOffset+1 );
      inc( kOffset, 2 );
    end;

    pr.pon_e( kPar, a );
  end;
end;



procedure TParsimonia.gp(var gPR: TMatR);
  var
  kOffset, j: integer;
  kPar: integer;
begin
  for kPar:= 1 to gPR.nc do
  begin
    kOffset:= 1;
    gPR.pon_e( kOffset, kPar, 1.0  );
    inc( kOffset );
    for j:= 1 to nPotencias do
    begin
       gPR.pon_e( kOffset, kPar, vpot.e( j ) );
       inc( kOffset );
    end;
    for j:= 1 to nArmonicos do
    begin
      gpR.pon_e( kOffset, kPar, vcos.e(j) );
      gpR.pon_e( kOffset+1, kPar, vsin.e(j) );
      inc( kOffset, 2 );
    end;
  end;


end;


destructor TParsimonia.destroy;
begin
  inherited destroy;
end;

end.

