unit ubibfun_onl;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils, xMatDefs, MatReal, Math;

type

  //** TIPOS ESPECIFICOS **
  //unidades básicas para formar bloques

  //todas heredan de una función genérica : TFunc
  TFunc = class
    nvars: NInt;
    //cantidad de variables involucradas en la función  = length(idVars)
    idVars: TDAOfNInt;
    //variables del vector "X" que pertenecen al dominio de la función
    //i.e. aquí X es el vector de la descripción global del sistema (es en general de dimensión >= nvars)
    procedure setIdVars(idV: TDAOfNInt);
    function Evaluar(var X: TVectR): NReal; virtual; abstract;
    function Gradiente(var X: TVectR): TVectR; virtual; abstract;
  end;

  TDAOfTFunc = array of TFunc;

  //POLINOMIO : sum_{i=0}^{i=grado} coef(i) * x^i
  TFunc_Polinomio = class(TFunc)
    grado: NInt;
    coef: TVectR;

    constructor Create_Init(gr: NInt; vcoef: TVectR);
    procedure Free;

    function Evaluar(var X: TVectR): NReal; override;
    function Gradiente(var X: TVectR): TVectR; override;
  end;

  TDAOfTFunc_Polinomio = array of TFunc_Polinomio;

  //RECTA : a*x + b
  TFunc_Recta = class(TFunc)
    a, b: NReal;

    constructor Create_Init(a_, b_: NReal);
    procedure Free;

    function Evaluar(var X: TVectR): NReal; override;
    function Gradiente(var X: TVectR): TVectR; override;
  end;

  TDAOfTFunc_Recta = array of TFunc_Recta;

  //COMBINACION LINEAL : ( sum_{i=1}^{i=nvar} coef(i) * x_i ) + k
  //la constante k es un término independiente opcional
  TFunc_CombinacionLineal = class(TFunc)
    coef: TVectR;
    k: NReal;

    constructor Create_Init(n: NInt; vcoef: TVectR); overload;
    constructor Create_Init(n: NInt; vcoef: TVectR; ko: NReal); overload;
    procedure Free;

    function Evaluar(var X: TVectR): NReal; override;
    function Gradiente(var X: TVectR): TVectR; override;
  end;

  TDAOfTFunc_CombinacionLineal = array of TFunc_CombinacionLineal;


  //PRODUCTORIA : k * prod_{i=1}^{i=nvars} x(i)^{exponente(i)}
  TFunc_Productoria = class(TFunc)
    k: NReal;
    exponentes: TVectR;

    constructor Create_Init(n: NInt; ko: NReal; expon: TVectR); overload;
    constructor Create_Init(n: NInt; expon: TVectR); overload;
    procedure Free;

    function Evaluar(var X: TVectR): NReal; override;
    function Gradiente(var X: TVectR): TVectR; override;
  end;

  TDAOfTFunc_Productoria = array of TFunc_Productoria;

  //TIPO UD (User Defined): función de Rn en R ingresada por el usuario
  TFunc_iRNenR = function(x: TVectR): NReal;
  TFunc_ikRNenR = function(x: TVectR; k: NInt): NReal;

  TFunc_UD = class(TFunc)
    f: TFunc_iRNenR;
    g: array of TFunc_ikRNenR;

    constructor Create_Init(nvar: NInt);
    procedure Free;

    procedure InscribirEcuacion(fi: TFunc_iRNenR);
    procedure InscribirDerivada(gij: TFunc_ikRNenR; j: NInt);

    function Evaluar(var X: TVectR): NReal; override;
    function Gradiente(var X: TVectR): TVectR; override;
  end;

  //** TIPO BLOQUE **
  //los tipos específicos se pueden ensamblar en un bloque: combinación lineal de funciones
  //por defecto: Agregar_f(fun) asume el coeficiente es 1, de lo contrario se debe especificar el vector en coef
  TFunc_Bloque = class(TFunc)
    f: array of TFunc;
    nbloques: NInt;
    coef: array of NReal;

    constructor Create;
    procedure Free;

    procedure Agregar_f(fun: TFunc); overload;
    procedure Agregar_f(fun: TFunc; c: NReal); overload;
    function Evaluar(var X: TVectR): NReal; override;
    function Gradiente(var X: TVectR): TVectR; override;
  end;

  TDAOfTFunc_Bloque = array of TFunc_Bloque;

  //** TIPO BLOQUE PRODUCTORIA **
  //una función que es un producto de funciones
  // c0 * prod_{i=1}^{i=nbloques} fi
  TFunc_Bloque_Productoria = class(TFunc)
    f: array of TFunc;
    nbloques: NInt;
    c: NReal;

    constructor Create_Init(c0: NReal);
    constructor Create();
    procedure Free;

    procedure Agregar_f(fun: TFunc);
    function Evaluar(var X: TVectR): NReal; override;
    function Gradiente(var X: TVectR): TVectR; override;
  end;

  TDAOfTFunc_Bloque_Productoria = array of TFunc_Bloque_Productoria;

  //= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
  //esta clase se utiliza en los algoritmos de optimización para especificar
  //un conjunto arbitrario de parámetros, que pueden ser vectores (Vec)
  //o números (Num)
  TParams_onl = class
    Vec: array of TVectR;
    Num: array of NReal;

    constructor Create_Init(nvec: NInt; nnum: NInt);
    procedure Free;
  end;

//______________________________________________________________________________

implementation

procedure TFunc.setIdVars(idV: TDAOfNInt);
var
  i: NInt;
  n: NInt;
begin
  n := length(idV);
  setlength(self.idVars, n);
  self.nvars := n;
  for i := 1 to n do
    self.idVars[i - 1] := idV[i - 1];
end;

constructor TFunc_Polinomio.Create_Init(gr: NInt; vcoef: TVectR);
begin
  inherited Create;
  Self.grado := gr;
  self.coef := TvectR.Create_Clone(vcoef);
  self.nvars := 1;
end;

procedure TFunc_Polinomio.Free;
begin
  setlength(self.idVars, 0);
  Self.coef.Free;
  inherited Free;
end;

function TFunc_Polinomio.Evaluar(var x: TVectR): NReal;
var
  a: TVectR;
  b: TVectR;
  i: NInt;
  xe: NReal;
  gr: NInt;
begin
  xe := x.e(self.idVars[0]);
  gr := Self.grado;
  a := TvectR.Create_Clone(self.coef);
  b := TvectR.Create_Init(gr);

  // método de Horner
  b.pon_e(1, a.e(gr + 1) * xe + a.e(gr));
  for i := 1 to gr - 1 do
  begin
    b.pon_e(i + 1, b.e(i) * xe + a.e(gr - i));
  end;
  Result := b.e(gr);
end;

function TFunc_Polinomio.Gradiente(var X: TVectR): TVectR;
var
  xa: TVectR;
  xe: NReal;
  a, b: TVectR;
  gr, i: NInt;
begin
  xe := x.e(self.idVars[0]);
  xa := TVectR.Create_Init(1);
  gr := Self.grado;

  //construyo un polinomio y luego lo evalúo con el método de Horner
  a := TVectR.Create_Init(self.grado); //tiene un coef menos que el original
  for i := 1 to gr do
  begin
    a.pon_e(i, self.coef.e(i + 1) * i);
  end;

  gr := gr - 1;
  if (gr > 0) then
  begin
    b := TvectR.Create_Init(gr);
    // método de Horner
    b.pon_e(1, a.e(gr + 1) * xe + a.e(gr));
    for i := 1 to gr - 1 do
    begin
      b.pon_e(i + 1, b.e(i) * xe + a.e(gr - i));
    end;
    xa.pon_e(1, b.e(gr));
    Result := xa;
  end
  else
  begin
    xa.pon_e(1, a.e(1));
    Result := xa;
  end;
end;

constructor TFunc_Recta.Create_Init(a_, b_: NReal);
begin
  inherited Create;
  self.nvars := 1;
  self.a := a_;
  self.b := b_;
end;

procedure TFunc_Recta.Free;
begin
  setlength(idVars, 0);
  inherited Free;
end;

function TFunc_Recta.Evaluar(var x: TVectR): NReal;
var
  i: NInt;
  xa: NReal;
begin
  xa := x.e(self.idVars[0]);
  Result := self.a * xa + self.b;
end;

function TFunc_Recta.Gradiente(var X: TVectR): TVectR;
var
  xa: TVectR;
begin
  xa := TVectR.Create_Init(1);
  xa.pon_e(1, self.a);
  Result := xa;
end;

constructor TFunc_CombinacionLineal.Create_Init(n: NInt; vcoef: TVectR);
begin
  inherited Create;
  self.nvars := n;
  self.coef := TvectR.Create_Clone(vcoef);
  self.k := 0;
end;

constructor TFunc_CombinacionLineal.Create_Init(n: NInt; vcoef: TVectR; ko: NReal);
begin
  inherited Create;
  self.nvars := n;
  self.coef := TvectR.Create_Clone(vcoef);
  self.k := ko;
end;

procedure TFunc_CombinacionLineal.Free;
begin
  Self.coef.Free;
  setlength(idVars, 0);
  inherited Free;
end;

function TFunc_CombinacionLineal.Evaluar(var x: TVectR): NReal;
var
  i: NInt;
  xa: NReal;
begin
  xa := 0;
  for i := 1 to self.nvars do
  begin
    xa := xa + self.coef.e(i) * x.e(self.idVars[i - 1]);
  end;
  Result := xa + self.k;
end;

function TFunc_CombinacionLineal.Gradiente(var X: TVectR): TVectR;
var
  xa: TVectR;
begin
  xa := TVectR.Create_Clone(self.coef);
  Result := xa;
end;

constructor TFunc_Productoria.Create_Init(n: NInt; ko: NReal; expon: TVectR);
begin
  inherited Create;
  Self.k := ko;
  self.nvars := n;
  self.exponentes := TVectR.Create_Clone(expon);
end;

constructor TFunc_Productoria.Create_Init(n: NInt; expon: TVectR);
begin
  inherited Create;
  Self.k := 1;
  self.nvars := n;
  self.exponentes := TVectR.Create_Clone(expon);
end;

procedure TFunc_Productoria.Free;
begin
  Self.exponentes.Free;
  setlength(idVars, 0);
  inherited Free;
end;

function TFunc_Productoria.Evaluar(var X: TVectR): NReal;
var
  ar: NReal;
  i: NInt;
  xe: NReal;
begin
  ar := 1;
  for i := 1 to self.nvars do
  begin
    xe := X.e(idVars[i - 1]);
    if not (xe = 0) then
      ar := ar * power(xe, Self.exponentes.e(i))
    else
    begin
      ar := 0;
      break;
    end;
  end;
  ar := ar * Self.k;
  Result := ar;
end;

function TFunc_Productoria.Gradiente(var X: TVectR): TVectR;
var
  ax: TVectR;
  i, j: NInt;
  ar: NReal;
begin
  ax := TVectR.Create_Init(length(idVars));
  for i := 1 to self.nvars do
  begin
    ar := 1;
    for j := 1 to self.nvars do
    begin
      if not (i = j) then
      begin
        ar := ar * power(X.e(self.idVars[j - 1]), Self.exponentes.e(j));
      end
      else
        ar := ar * self.exponentes.e(j) * power(x.e(self.idVars[j - 1]),
          self.exponentes.e(j) - 1);
    end;
    ar := ar * self.k;
    ax.pon_e(i, ar);
  end;
  Result := ax;
end;

function fnula_i(x: TVectR): NReal;
begin
  Result := 0;
end;

function fnula_ij(x: TVectR; k: NInt): NReal;
begin
  Result := 0;
end;

constructor TFunc_UD.Create_Init(nvar: NInt);
var
  k: NInt;
begin
  inherited Create;
  self.nvars := nvar;
  setlength(Self.g, nvar);
  self.f := @fnula_i;
  for k := low(self.g) to high(self.g) do
    self.g[k] := @fnula_ij;

end;

procedure TFunc_UD.Free;
begin
  setlength(Self.g, 0);
  setlength(idVars, 0);
  inherited Free;
end;

procedure TFunc_UD.InscribirEcuacion(fi: TFunc_iRNenR);
begin
  Self.f := fi;
end;

procedure TFunc_UD.InscribirDerivada(gij: TFunc_ikRNenR; j: NInt);
begin
  Self.g[j - 1] := gij;
end;

function TFunc_UD.Evaluar(var X: TVectR): NReal;
begin
  Result := Self.f(X);
end;

function TFunc_UD.Gradiente(var X: TVectR): TVectR;
var
  ax: TVectR;
  k: NInt;
begin
  ax := TVectR.Create_Init(self.nvars);

  for k := 1 to self.nvars do
    ax.pon_e(k, Self.g[k - 1](X, self.idVars[k - 1]));

  Result := ax;
end;

constructor TFunc_Bloque.Create;
begin
  inherited Create;
  self.nvars := 0;
  Self.nbloques := 0;
  setlength(Self.f, 0);
  setlength(self.coef, 0);
end;

procedure TFunc_Bloque.Free;
begin
  setlength(Self.f, 0);
  setlength(self.coef, 0);
  setlength(idVars, 0);
  inherited Free;
end;

procedure TFunc_Bloque.Agregar_f(fun: TFunc);
var
  i, j: NInt;
  ES_NUEVA: boolean;
begin
  //¿fun introduce variables nuevas al bloque?
  for i := 1 to fun.nvars do
  begin
    ES_NUEVA := True;
    for j := 1 to self.nvars do
    begin
      if (self.idVars[j - 1] = fun.idVars[i - 1]) then
      begin
        ES_NUEVA := False;
        break;
      end;
    end;

    if ES_NUEVA then
    begin
      setlength(self.idVars, self.nvars + 1);
      self.idVars[high(idVars)] := fun.idVars[i - 1];
      self.nvars := self.nvars + 1;
    end;
  end;
  //agregar función al bloque
  setlength(Self.f, Self.nbloques + 1);
  Self.f[High(Self.f)] := fun;
  Self.nbloques := Self.nbloques + 1;
  setlength(self.coef, length(self.coef) + 1);
  self.coef[high(self.coef)] := 1;
end;

procedure TFunc_Bloque.Agregar_f(fun: TFunc; c: NReal);
var
  i, j: NInt;
  ES_NUEVA: boolean;
begin
  //¿fun introduce variables nuevas al bloque?
  for i := 1 to fun.nvars do
  begin
    ES_NUEVA := True;
    for j := 1 to self.nvars do
    begin
      if (self.idVars[j - 1] = fun.idVars[i - 1]) then
      begin
        ES_NUEVA := False;
        break;
      end;
    end;

    if ES_NUEVA then
    begin
      setlength(self.idVars, self.nvars + 1);
      self.idVars[high(idVars)] := fun.idVars[i - 1];
      self.nvars := self.nvars + 1;
    end;
  end;
  //agregar función al bloque
  setlength(Self.f, Self.nbloques + 1);
  Self.f[High(Self.f)] := fun;
  Self.nbloques := Self.nbloques + 1;
  setlength(self.coef, length(self.coef) + 1);
  self.coef[high(self.coef)] := c;
end;

function TFunc_Bloque.Evaluar(var X: TVectR): NReal;
var
  k: NInt;
  ax: NReal;
begin
  ax := 0;
  for k := 1 to Self.nbloques do
  begin
    ax := ax + self.coef[k - 1] * self.f[k - 1].Evaluar(X);
  end;
  Result := ax;
end;

function TFunc_Bloque.Gradiente(var X: TVectR): TVectR;
var
  i, j, k: NInt;
  xa, xb: TVectR;
  pos: NInt;
  xx1, xx2: NReal;
  xint1: NInt;
begin
  xa := TVectR.Create_Init(self.nvars);

  //recorro todos los bloques: lo importante es ubicar correctamente a la derivada
  //en el vector xa
  for i := 1 to self.nbloques do
  begin

    xb := TVectR.Create_Init(self.f[i - 1].nvars);
    xb := self.f[i - 1].Gradiente(X);

    xx1 := xb.e(1); //DEBUGGG


    xb.PorReal(self.coef[i - 1]);


    xint1 := xb.n; //DEBUGUGUGUG
    xx1 := xb.e(1); //DEBUGGG

    //recorrer las variables de la función que derivé
    for j := 1 to self.f[i - 1].nvars do
    begin
      //encontrar dónde las tengo que ubicar en el vector de idvars de TFunc
      for k := 1 to self.nvars do
      begin
        if (self.f[i - 1].idVars[j - 1] = self.idVars[k - 1]) then
        begin
          //lo encontré
          pos := k;
          break;
        end;
      end;
      xa.acum_e(pos, xb.e(j));
    end;

  end;
  Result := xa;
end;

constructor TFunc_Bloque_Productoria.Create();
begin
  inherited Create;
  self.nvars := 0;
  Self.nbloques := 0;
  setlength(Self.f, 0);
  self.c := 1;
end;

constructor TFunc_Bloque_Productoria.Create_Init(c0: NReal);
begin
  inherited Create;
  self.nvars := 0;
  setlength(self.idVars, 0);
  Self.nbloques := 0;
  setlength(Self.f, 0);
  self.c := c0;
end;

procedure TFunc_Bloque_Productoria.Free;
begin
  setlength(Self.f, 0);
  setlength(self.idVars, 0);
  inherited Free;
end;

procedure TFunc_Bloque_Productoria.Agregar_f(fun: TFunc);
var
  i, j: NInt;
  ES_NUEVA: boolean;
begin
  //¿fun introduce variables nuevas al bloque?
  for i := 1 to fun.nvars do
  begin
    ES_NUEVA := True;
    for j := 1 to self.nvars do
    begin
      if (self.idVars[j - 1] = fun.idVars[i - 1]) then
      begin
        ES_NUEVA := False;
        break;
      end;
    end;

    if ES_NUEVA then
    begin
      setlength(self.idVars, self.nvars + 1);
      self.idVars[high(idVars)] := fun.idVars[i - 1];
      self.nvars := self.nvars + 1;
    end;
  end;
  //agregar función al bloque
  setlength(Self.f, Self.nbloques + 1);
  Self.f[High(Self.f)] := fun;
  Self.nbloques := Self.nbloques + 1;

end;

function TFunc_Bloque_Productoria.Evaluar(var X: TVectR): NReal;
var
  k: NInt;
  ax, bx: NReal;
begin

  ax := 1;
  for k := 1 to Self.nbloques do
  begin
    bx := self.f[k - 1].Evaluar(x);
    if not (bx = 0) then
      ax := ax * bx
    else
    begin
      ax := 0;
      break;
    end;
  end;
  Result := ax * self.c;
end;

function TFunc_Bloque_Productoria.Gradiente(var X: TVectR): TVectR;
var
  i, j, k: NInt;
  xa, xb: TVectR;
  ax, bx: NReal;
begin
  xa := TVectR.Create_Init(self.nvars);

  for i := 1 to self.nbloques do
  begin

    xb := TVectR.Create_Init(self.f[i - 1].nvars);
    xb := self.f[i - 1].Gradiente(X);

    //a cada elemento de xb lo tengo que multiplicar por el resto de las funciones
    //evaluadas en X
    ax := 1;
    for k := 1 to Self.nbloques do
    begin
      if not (i = k) then
      begin
        bx := self.f[k - 1].Evaluar(X);
        if not (bx = 0) then
          ax := ax * bx
        else
        begin
          ax := 0;
          break;
        end;
      end;
    end;
    xb.PorReal(ax);
    //calculo el aporte de la función self.f[i-1] a la componente j del
    //gradiente
    for j := 1 to self.f[i - 1].nvars do
    begin
      //xa.acum_e(self.f[i-1].idVars[j-1], xb.e(j));
      xa.acum_e(j, xb.e(j));
    end;
  end;
  xa.PorReal(self.c);
  Result := xa;
end;

constructor TParams_onl.Create_Init(nvec: NInt; nnum: NInt);
begin
  inherited Create;
  setlength(Self.Vec, nvec);
  setlength(Self.Num, nnum);
end;

procedure TParams_onl.Free;
begin
  setlength(Self.Vec, 0);
  setlength(Self.Num, 0);
  inherited Free;
end;

end.
