unit umixturedensitynetworks;

{$mode delphi}

interface

uses
  Classes, SysUtils, xmatdefs, uneuronas, Math, math01, utipos_xcf, matreal,
  usigmoide,
  fddp_weibull, udisnormcan;

type

  { TKernel_PDF }

  TKernel_PDF = class
    constructor Create(nParametros: integer);
    function nparams: integer; // cuenta la cantidad de parámetros necesarios.
    procedure set_input(x: NReal); virtual; abstract; // guarda x y calcula y
    function get_pdf: NReal; // retorna y
    function get_dpdk(kParam: integer): NReal; virtual;
      abstract;// retorna la derivada dp(x)/dparam(k)
    destructor Destroy; override;

    function GetParam(k: integer): NReal;
    procedure SetParam(k: integer; val: NReal);

    // auxiliares para cálculo
    //
    // Distribución
    function F( x: NReal ): NReal; virtual; abstract;
    // invesa de la distribución
    function invF( u: NReal ): NReal; virtual; abstract;


  private
    x: NReal; // último valro fijado de la entrada
    y: NReal; // último valor calculado de pdf por eficiencia del cálculo de las derivadas.
    params: TDAOfNReal;
    procedure CalcAuxs; virtual;

  end;

  TDAOfKPDF = array of TKernel_PDF;

  { TKernel_Gaussian }
  TKernel_Gaussian = class(TKernel_PDF)
    constructor Create(mu, sigma: NReal);
    procedure set_input(x: NReal); override;
    function get_dpdk(kParam: integer): NReal; override;

    // Distribución
    function F( x: NReal ): NReal; override;
    // invesa de la distribución
    function invF( u: NReal ): NReal; override;

  private
    u, s, s2, s3, inv_s, C_: NReal;
    procedure CalcAuxs; override;
  end;

  { TKernel_Weibull }
  TKernel_Weibull = class(TKernel_PDF)
    constructor Create(lambda, k: NReal); // escala y forma
    procedure set_input(x: NReal); override;
    function get_dpdk(kParam: integer): NReal; override;

    // Distribución
    function F( x: NReal ): NReal; override;
    // invesa de la distribución
    function invF( u: NReal ): NReal; override;

  private
    k, l, inv_k: NReal;
    z: NReal;
    procedure CalcAuxs; override;
    function dw_dk: NReal;
    function dw_dl: NReal;
  end;


  { TKernel_DevSig } // Derivado de Sigmoide

  TKernel_DevSig = class(TKernel_PDF)
    constructor Create( x0, alfa: NReal);
    procedure set_input( x: NReal ); override;
    function get_dpdk( kParam: integer ): NReal; override;

    // Distribución
    function F( x: NReal ): NReal; override;
    // invesa de la distribución
    function invF( u: NReal ): NReal; override;

  private
    x0, alfa, inv_alfa, F_, dFdz: NReal;
    procedure CalcAuxs; override;
  end;




  { TMixtureDensityNetwork }

  TMixtureDensityNetwork = class
    Cerebro: TCerebro;
    kernels: array of TDAOfKPDF;
    pesos: array of TDAOfNReal; // pesos para combinar las funciones del kernel
    aData: TDAOfRecXCF;
    DimKernel: integer;

    // auxiliares para indexar las salidas de la red
    NEntradas, NSalidasRN, NSalidas: integer;
    Offsetparam1, Offsetparam2: integer;

    constructor Create( aData: TDAOfRecXCF; nCapasOcultas, DimKernel: integer);


    procedure SetEntradas( X: TVectR );
    procedure SetEntradas_prnd( X: TVectR; prnd: NReal );
    procedure evaluar;

    // Calcula el valor de la varialbe explicada y del vector explicado Y
    // que es excedido con un Probabilidad de Excedencia PE, para el último
    // valor del vector explicativo X fijado con SetEntradas y evaluado con evaluar
    // jSalida = 0 .. NSalidas-1
    function yOfPE( PE: NReal; jSalida: integer ): NReal;

    // Calcula la probabilidad de excedencia del valor y del vector explicado Y
    // para el valor del vector explicativo X fijado con SetEntradas y evaluado con evaluar.
    function PEOfy( y: NReal; jSalida: integer ): NReal;


    function EntrenarCerebro_EBP_minibatch( var grad: TVectR;
      paso, momento, alfa_Regularizacion: NReal; NMuestras: integer;
    // si NMuestras <= 0 ejecuta todas las muestras una vez.
      semilla: integer; prnd: NReal =
      0.0  // porcion de ruido para afectar las entradas
      ): NReal;


    destructor Destroy; override;

    // función auxiliar para calcular yOfPE
    function PEMenosPEObjetivo( y: NReal ): NReal;

    private
      y,  rn_y: TVectR;
      uObjetivo: NReal;
      jSalidaObjetivo: integer;


  end;


implementation

{ TKernel_DevSig }

constructor TKernel_DevSig.Create( x0, alfa: NReal);
begin
  inherited Create(2); // x0 = params[0], alfa = params[1]
  params[0] := x0;
  params[1] := alfa;
  CalcAuxs;
end;

procedure TKernel_DevSig.set_input(x: NReal);
begin
  self.x:= x;
  F_:= sigmoide( alfa * ( x - x0 ) );
  y:= F_*( 1.0 - F_ );
  dFdz :=  y * ( 1.0 - 2.0 * F_);
end;

function TKernel_DevSig.get_dpdk(kParam: integer): NReal;
begin
  case kParam of
    0: // dp/dxo
      result:= -alfa * dFdz;
    1: // dp/dalfa
      result:= dFdz * ( x - x0 );
  end;
end;

function TKernel_DevSig.F(x: NReal): NReal;
begin
  result:= F_;
end;

function TKernel_DevSig.invF( u: NReal ): NReal;
begin
  result:= inv_sigmoide( u ) * inv_alfa + x0;
end;

procedure TKernel_DevSig.CalcAuxs;
begin
  x0:= params[0];
  alfa:= params[1];
  inv_alfa:= 1.0 / alfa;
end;

{ TMixtureDensityNetwork }

constructor TMixtureDensityNetwork.Create(aData: TDAOfRecXCF;
  nCapasOcultas, DimKernel: integer);

var
  k, j: integer;
  CapasDef: TDAofNInt;
  TiposSat: TDaOfTipoSaturacion;
  aRec: TRecXCF;

begin
  self.aData := aData;
  self.DimKernel := DimKernel;

  aRec := aData[0];
  NEntradas := aRec.X.n;
  NSalidas := aRec.Y.n;
  NSalidasRN := NSalidas * (DimKernel * 3); // peso y dos parámetros para la PDF

  Offsetparam1 := DimKernel * NSalidas;
  Offsetparam2 := Offsetparam1 * 2;

  setlength(CapasDef, nCapasOcultas + 1);
  setlength(TiposSat, nCapasOcultas + 1);
  for k := 0 to high(CapasDef) do
    CapasDef[k] := NSalidasRN;

  for k := 0 to high(TiposSat) - 1 do
    TiposSat[k] := CTS_Tanh;
  TiposSat[high(TiposSat)] := CTS_Lineal;

  Cerebro := TCerebro.Create(nEntradas, CapasDef, TiposSat, False);
  randomize;
  Cerebro.RandomInit(0.1);
  Cerebro.EscalarEntradas(aData);

  // supongamos que en las salidas de la RN están primero
  // las N salidas en base a las cuales por SoftMax endremos Los Pesos a
  // considerar de las f del kernel.
  // luego las N correspondientes al parámetro CERO de las f del kernel
  // y luego las N correspondientes al parámetro UNO del las f del kernel.

  // para el caso del kernel gaussiano, el param[0] es el valor esperado y
  // el param[1] el desvío estándar (que es siempre > 0).
  // para tener en cuenta este aspecto, la salida de la RN correspondiente al desvío
  // pa pasaremos por exp( salida ) para tener así el desvío.
  // En el caso de que el kernel sea con funciones de weibull, ambos parámtros son > 0 por lo
  // ambas salidas de la RN la pasaremos por (exp( salida ) ) para calcular el parámetro.

  setlength(Kernels, NSalidas);
  setlength(pesos, NSalidas);
  for k := 0 to high(kernels) do
  begin
    setlength(kernels[k], DimKernel);
    setlength(pesos[k], DimKernel);
    for j := 0 to high(kernels[k]) do
//    kernels[k][j] := TKernel_Weibull.Create(random, random);
    kernels[k][j] := TKernel_Gaussian.Create(random, random);
//    kernels[k][j] := TKernel_DevSig.Create(random, random);
  end;

  y := TVectR.Create_Init(NSalidas);
  rn_y := TVectR.Create_Init(NSalidasRN);


end;

procedure TMixtureDensityNetwork.SetEntradas(X: TVectR);
begin
  Cerebro.SetEntradas( X );
end;

procedure TMixtureDensityNetwork.SetEntradas_prnd(X: TVectR; prnd: NReal);
begin
  Cerebro.SetEntradas_prnd( X, prnd );
end;

procedure TMixtureDensityNetwork.evaluar;
var
  jSalida, j: integer;
  aKernelFunc: TKernel_PDF;

begin
    Cerebro.evaluar;
    Cerebro.GetSalidas(rn_Y);
    for jSalida := 1 to NSalidas do
    begin
      for j := 1 to DimKernel do
      begin
        // Leemos las salidas para los pesos y los calculamos
        pesos[jSalida - 1][j - 1] := rn_Y.e( (jSalida - 1) * DimKernel + j );
        aKernelFunc := kernels[jSalida - 1][j - 1];
        if aKernelFunc is TKernel_Weibull then
          aKernelFunc.SetParam( 0, exp(rn_y.e(OffsetParam1 + (jSalida - 1) * DimKernel + j)) )
        else
          aKernelFunc.SetParam( 0, rn_y.e( OffsetParam1 + (jSalida - 1) * DimKernel + j ));

        if aKernelFunc is TKernel_DevSig then
          aKernelFunc.SetParam(1, rn_y.e( OffsetParam2 + (jSalida - 1) * DimKernel + j ))
        else
          aKernelFunc.SetParam(1, exp( rn_y.e( OffsetParam2 + (jSalida - 1) * DimKernel + j )));
      end;
      AutoSoftMax( pesos[jSalida - 1] );

      (*
     write( 'Pesos: ', DAOfNRealToStr_( pesos[jSalida-1], 4,3,',') );
      write( ': ');
      for j:= 1 to DimKernel do
      begin
        aKernelFunc := kernels[jSalida - 1][j - 1];
        write( ': ', DAOfNRealToStr_( aKernelFunc.params, 4, 3, ',' ) );
      end;
      writeln;
        *)
    end;

end;


function TMixtureDensityNetwork.yOfPE(PE: NReal; jSalida: integer): NReal;
var
  yLeft, yRight, y, m: NReal;
  jFunc: integer;
  aKernelFunc: TKernel_PDF;
  dicot: T_ProblemaFx_Cero_externaOfObject;
  res: boolean;
  fAtRoot: NReal;
  NoOfIters: word;

begin
  uObjetivo:= 1- PE;
  jSalidaObjetivo:= jSalida;

  aKernelFunc := kernels[jSalida][0];
  yLeft:= aKernelFunc.invF( uObjetivo );
  yRight:= yLeft;
  for jFunc:= 1 to DimKernel-1 do
  begin
    aKernelFunc := kernels[jSalida][jFunc];
    m:= aKernelFunc.invF( uObjetivo );
    if m < yLeft then
       yLeft:= m
    else if m > yRight then
       yRight := m;
  end;

  m:= yRight - yLeft;
  if m <= AsumaCero then
  begin
    result:= ( yRight + yLeft )/ 2.0;
    exit;
  end;

  // bien ahora sabemos que en la combinación de funciones la solución está
  // entre yLeft e yRight
  dicot:= T_ProblemaFx_Cero_externaOfObject.Create( PEMenosPEObjetivo );
  res:= dicot.Dicot( yLeft, yRight, max( AsumaCero, m*1.0E-5 ), 500, y, fAtRoot, NoOfIters );

  result:= y;
end;


function TMixtureDensityNetwork.PEOfy(y: NReal; jSalida: integer): NReal;
var
  aKernelFunc: TKernel_PDF;
  j: integer;
  cdf: NReal;
begin
  cdf := 0;
  for j := 1 to DimKernel do
  begin
    aKernelFunc := kernels[jSalida][j - 1];
    cdf := cdf + pesos[jSalida][j - 1] * aKernelFunc.F( y );
  end;
  result:= 1-cdf;
end;



function TMixtureDensityNetwork.EntrenarCerebro_EBP_minibatch(
  var grad: TVectR;
  paso, momento, alfa_Regularizacion: NReal; NMuestras: integer;
  semilla: integer; prnd: NReal): NReal;
var
  kMuestra: integer;
  xRec: TRecXCF;
  e: TVectR;
  k, j, jSalida, h: integer;
  flg_todas: boolean;
  flg_InyectarRuido: boolean;
  aKernelFunc: TKernel_PDF;
  m, fdp, b, JCosto, dbdp, f_j: NReal;

  function rnd(k: integer): integer;
  begin
    if flg_todas then
      Result := k
    else
      Result := random(length(aData));
  end;

begin
  RandSeed := semilla;
  e := TVectR.Create_Init(Cerebro.NSalidas);
  flg_InyectarRuido := prnd <> 0;

  if NMuestras <= 0 then
  begin
    NMuestras := length(aData);
    flg_todas := True;
  end
  else
    flg_todas := False;

  Cerebro.BackPropagation_init;
  JCosto := 0;

  for k := 1 to NMuestras do
  begin
    Cerebro.BackPropagation_ClearErr;
    kMuestra := rnd(k - 1);
    xRec := aData[kMuestra];

    if flg_InyectarRuido then
      SetEntradas_prnd(xRec.X, prnd)
    else
      SetEntradas(xRec.X);

    evaluar;

    // Ahora calculamos las probabilidades de la salida observada
    for jSalida := 1 to NSalidas do
    begin
      fdp := 0;
      for j := 1 to DimKernel do
      begin
        aKernelFunc := kernels[jSalida - 1][j - 1];
        aKernelFunc.set_input( xRec.Y.e(jSalida));
        fdp := fdp + pesos[jSalida - 1][j - 1] * aKernelFunc.get_pdf;
      end;
      y.pon_e(jSalida, fdp);
    end;


    // Ahora, hay que hacer el gradiente de la función objetivo utilizando la regla
    // de la cadena hasta la entrada de la RN para aplicar backpropagation en la RN.
    // El objetivo es maximimzar la verosimilitud L = prod p(x_k, y_k ) que es equivalente
    // a minimizar J = - sum ln( p(x_k, y_k) )  en el conjunto de muestras (x_k, y_k)
    // Ahora bien, -sum( ln( p(x_k, y_k) ) ) = -sum( ln( p(y_k/x_k) ) ) -sum( ln( p( x_k) ) )
    // y como el término -sum( ln(px_k) ) es independiente de los parámetros de entrenamiento no
    // interviene en la optimización; que se reduce a: min( -sum( ln( p(y_k/x_k) ) ) )
    // Para el cálculo del gradiente, la sumatoria nos permite calcular el gradiente sobre
    // cada muestra debiendo hacer la derivada de b = -ln( p(y/x) ) para cada muestra (x, y)


    // limpiamos el vector de error
    e.Ceros;
    for jSalida := 1 to NSalidas do
    begin
      b := y.e(jSalida);
      JCosto := JCosto - ln(b);
      // d_b/d_p(y/x) = -1/p(y/x)
      dbdp := -1.0/b;

      // A su vez p= p(y/x) es la suma ponderada de las funciones del kernel
      // p = sum_j( pesos_j * f_j )
      // Llamemos z0_j, z1_j y z2_j a las salidas de la N donde el subíndice j identifica
      // la función del kernel.
      // Los pesos son una combinación softmax de los valores z0 por lo cual
      // d_p/d_z0_h = sum_j( d_pesos_j/d_z0_h * f_j )
      // donde d_pesos_j/d_z0_h =  pesos_j * ( (j=h) - pesos_h )
      for j := 1 to dimKernel do
      begin
        aKernelFunc := kernels[jSalida - 1][j - 1];
        f_j := aKernelFunc.get_pdf;

        for h := 1 to dimKernel do
        begin
          if h = j then
            m := pesos[jSalida - 1][j - 1] * (1 - pesos[jSalida - 1][h - 1])
          else
            m := -pesos[jSalida - 1][j - 1]*pesos[jSalida - 1][h - 1];
          e.acum_e((jSalida - 1) * DimKernel + h, m * f_j * dbdp);
        end;
      end;

      // Ahora las derivadas sobre param1
      // en el caso de la weibull hay un exp() en el medio.
      // d_p/d_z1_h = sum_j( pesos_j * df_j/d_z1_h )
      // df_j/d_z1_h = 0 si h <> j y df_j/d_param1 * d_param1/d_z1_j si h = j
      for j := 1 to dimKernel do
      begin
        aKernelFunc := kernels[jSalida - 1][j - 1];
        m := pesos[jSalida - 1][j - 1] * aKernelFunc.get_dpdk(0);

        if aKernelFunc is TKernel_Weibull then
          m := m * aKernelFunc.GetParam(0);
        // para tener en cuenta la derivada de exp( z1_j )

        e.acum_e(Offsetparam1 + (jSalida - 1) * DimKernel + j, m * dbdp);
      end;

      // ahora las derivadas sobre param2
      // tanto en gaussiana como en weibull va un exp( z2 ) en el medio.
      for j := 1 to dimKernel do
      begin
        aKernelFunc := kernels[jSalida - 1][j - 1];
        m := pesos[jSalida - 1][j - 1] * aKernelFunc.get_dpdk(1);

         if not ( aKernelFunc is TKernel_DevSig ) then
          m := m * aKernelFunc.GetParam(1);
          // para tener en cuenta la derivada de exp( z2_j )

        e.acum_e(Offsetparam2 + (jSalida - 1) * DimKernel + j, m * dbdp);
      end;
    end;

    Cerebro.BackPropagation_addMuestra(e);
  end;

  Cerebro.BackPropagation_DarPaso(paso, momento, 1 - momento, alfa_Regularizacion);

  Result := JCosto / NMuestras;
  // writeln( result) ;
  e.Free;
end;

destructor TMixtureDensityNetwork.Destroy;
var
  k, j: integer;

begin
  for k := 0 to high(kernels) do
  begin
    for j := 0 to high(kernels[k]) do
      TKernel_Gaussian(kernels[k][j]).Free;
    setlength(kernels[k], 0);
    setlength(pesos[k], 0);
  end;
  setlength(Kernels, 0);
  setlength(pesos, 0);

  y.Free;
  rn_y.Free;
  inherited Destroy;
end;

function TMixtureDensityNetwork.PEMenosPEObjetivo(y: NReal): NReal;
begin
  result:= (1-PEOfy( y, jSalidaObjetivo )) - uObjetivo;
end;

{ TKernel_Weibull }

constructor TKernel_Weibull.Create(lambda, k: NReal);
begin
  inherited Create(2);
  params[0] := lambda;
  params[1] := k;
  CalcAuxs;
end;

procedure TKernel_Weibull.set_input(x: NReal);
begin
  self.x := x;
  if x > 0 then
  begin
    z := power(x / l, k);
    y := k / x * z * exp(-z);
  end
  else
  begin
    z := 0;
    y := 0;
  end;
end;

function TKernel_Weibull.dw_dl: NReal;
begin
  //    z:= power(x/l,k);
  //    y:=  k/x* z *exp(-z )
  // dy_dl = dy_dz dz_dl
  //  dy_dz = y* ( 1/z - 1 )
  //  dz_dl = k power( x/l, k-1 ) ( -1/(l*l) )= -k*z/(l*x)
  if x > 0 then
    Result := y * k / (l * x) * (z - 1)
  else
    Result := 0;
end;

function TKernel_Weibull.dw_dk: NReal;
begin
  //    z:= power(x/l,k);
  //    y:=  k/x* z *exp(-z )
  // dy_dk = y * (1/k +  ( 1/z - 1 ) dz_dk )
  // dz_dk = k ln( x/l )
  if x > 0 then
    Result := y * (inv_k + (1.0 / z - 1.0) * k * ln(x / l))
  else
    Result := 0;
end;


function TKernel_Weibull.get_dpdk(kParam: integer): NReal;
begin
  if x < 0 then
    Result := 0
  else
    case kParam of
      0: Result := dw_dl;
      1: Result := dw_dk;
    end;
end;

function TKernel_Weibull.F(x: NReal): NReal;
begin
  result:= DistribucionWeibull( x, params[0], params[1] );
end;

function TKernel_Weibull.invF(u: NReal): NReal;
begin
  result:= Inv_DistribucionWeibull(x, params[0], params[1] );
end;

procedure TKernel_Weibull.CalcAuxs;
begin
  l := params[0];
  k := params[1];
end;

{ TKernel_Gaussian }

constructor TKernel_Gaussian.Create(mu, sigma: NReal);
begin
  inherited Create(2); // mu = params[0], sigma = params[1]
  params[0] := mu;
  params[1] := sigma;
  CalcAuxs;
end;

procedure TKernel_Gaussian.CalcAuxs;
begin
  u := params[0];
  s := params[1];
  s2 := s * s;
  s3 := s2 * s;
  inv_s := 1.0 / s;
  C_ := 1.0 / ( sqrt(2.0 * pi) * s  );
end;

procedure TKernel_Gaussian.set_input(x: NReal);
begin
  self.x := x;
  y := C_ * exp(-0.5 * sqr((x - u) * inv_s));
end;

function TKernel_Gaussian.get_dpdk(kParam: integer): NReal;
begin
  case kParam of
    0: Result := y * (x - u) / s2;
    1: Result := y * (sqr(x - u) / s3 - inv_s);
    else
      raise Exception.Create(ClassName + ', get_dpk kParam out of range: ' + IntToStr(
        kParam));
  end;

end;

function TKernel_Gaussian.F(x: NReal): NReal;
begin
  result:= DistribucionNormal( x, params[0], params[1] );
end;

function TKernel_Gaussian.invF(u: NReal): NReal;
begin
  result:= Inv_DistribucionNormal( u, params[0], params[1] );
end;


{ TKernel_PDF }

constructor TKernel_PDF.Create(nParametros: integer);
begin
  inherited Create;
  setlength(params, nParametros);
end;

function TKernel_PDF.nparams: integer;
begin
  Result := length(params);
end;


function TKernel_PDF.get_pdf: NReal;
begin
  Result := y;
end;


destructor TKernel_PDF.Destroy;
begin
  setlength(params, 0);
  inherited Destroy;
end;

function TKernel_PDF.GetParam(k: integer): NReal;
begin
  Result := params[k];
end;

procedure TKernel_PDF.SetParam(k: integer; val: NReal);
begin
  params[k] := val;
  CalcAuxs;
end;

procedure TKernel_PDF.CalcAuxs;
begin
  // nada;
end;

end.
