//módulo que tiene la clase que define el problema de optimización no lineal.
//para resolver el problema se deben ejecutar algoritmos externos.

unit udescripcionproblema_onl;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils, xMatDefs, matreal,
  ubibfun_onl, uresolvedor_onl, usimplex;

type
  TProblema_onl = class
public
    spx : TSimplex;   //problema simplex hijo

    mensajeDeError: string;

    nvars, nres, nres_desigualdad, nres_igualdad : NInt;

    x_inf, x_sup: TVectR; // restricciones de caja

    (*
       en esta implementación se minimiza la función de costo.
       Ojo! el simplex está programado para maximizar la función objetivo.
       No obstante, para los actores que cargan en simplex es
       transpartente porque el cambio de signo se realiza en este módulo, cuando
       se pasa la información del simplex a func bloques (en resolver()).

       Por tanto los actores que quieran cargar bloques de funciones no
       lineales deben cargar directamente su función de costo sin cambiarle
       el signo.
    *)
    f_costo               : TFunc_Bloque;        //función de costo
    g_resDesigualdad      : TDAOfTFunc_Bloque;   // >= 0     restricciones
    h_resIgualdad         : TDAOfTFunc_Bloque;   // = 0     restricciones

    XOpt     : TVectR;

    cnt_res_desigualdad : NInt;
    cnt_res_igualdad    : NInt;

	flg_x: array of shortint; // es 0 si no hay cota superior ni inferior fijadas
	  						  // es 1 si hay cotas superior e inferior fijadas
	 				    	  // es 2 si hay cota superior fijada y no cota inferior
                              // es 3 si no hay cota superior fijada pero sí hay cota inferior
                              // es 4 si esta variable fue fijada.

    constructor Create_init( ires, ivar: NInt);
    procedure Free;

    (* Fijamos que la restricción kfila es de igualdad *)
    procedure FijarRestriccionIgualdad( kfila: NInt );

    (* Fija el valor de una variable. Esto permite escribir las ecuaciones
	considerando la variable pero luego imponerle un valor.
    Debe ser llamado con el indice que tenía la variable cuando se cargo en el
    Simplex
    *)
    procedure FijarVariable( ivar: NInt; valor: NReal );

    // método para menejo de las restricciones de caja
    (*
    	Fija el valor de la cota inferior
    *)
    procedure cota_inf_set( ivar: NInt; vxinf: NReal );

    (*
    	Fija el valor de la cota superior
    *)
    procedure cota_sup_set( ivar: NInt; vxsup: NReal );

    (*
        Funciones que se incluyen por compatibilidad con Actores que usan el Simplex
    *)
    function e(k, j: integer): NReal;
    procedure pon_e(k, j: integer; x: NReal);
    procedure acum_e(k, j: integer; x: NReal);

    (*
        Funciones auxiliares para leer los resultados
    *)
    function xval( ix: NInt ): NReal;    //retorna el resutaldo de la variable x
    function yval( iy: NInt ): NReal;    //retorna el valor que toma la restricicón
    function xmult( ix: NInt ): NReal;   //retorna el multiplicador de Lagrange asociado a la variable
    function ymult( iy: NInt ): NReal;   //retorna el multiplicador asociado a la restricción irest
    function fval: NReal;                //retorna el valor de la función objetivo

    //procedure DumpSistemaToXLT( var f: textfile ); overload; virtual;
    //procedure DumpSistemaToXLT( archi: string; InfoAdicional: string ); overload; virtual;

    (*
    	Busca un punto de arranque factible y si lo encuentra maximiza la función fval
    	El resultado es 0 si logró encontrar un punto factible y realizar la maximización
    	si no se puede resolver el resultado es <> 0 y se guarda en la variable
    	"MensajeDeError" del objeto la causa encontrada.
    *)
    function resolver: NInt; virtual;

    procedure AgregarFObjetivo(f: TFunc);

    procedure AgregarRestriccionDeDesigualdad(g : TFunc);

    procedure AgregarRestriccionDeIgualdad(h : TFunc);

 end;

implementation

constructor TProblema_onl.Create_init(ires, ivar : NInt);
var
   i : NInt;
begin
      {$IFDEF SIMPLEX_HIJO}
  	 	Self.spx := usimplex.TSimplex.Create_init(ires,ivar);
     {$ENDIF}

     inherited Create;
     Self.nvars:=ivar-1;
     Self.nres:=ires-1;  //involucra: restricciones de igualdad, restricciones de desigualdad
                         //no involucra: restricciones de caja

     cnt_res_desigualdad:=0;
     cnt_res_igualdad:=0;

     x_inf := TVectR.Create_Init( ivar - 1 );
	 x_sup := TVectR.Create_Init( ivar - 1 );

     setlength(Self.flg_x, nvars);
     for i:=0 to nvars-1 do
       flg_x[i] := 0;

     f_costo := TFunc_Bloque.Create;

     //las restricciones de caja se van agregando a medida que se llaman a las funciones
     //apropiadas

     //todavía no sé cuáles de las restricciones son de igualdad y cuáles de desigualadad,
     //el setlength de los array con restricciones los voy actualizando a medida que
     //van ingresando restricciones.

     setlength(g_resDesigualdad, 0);

     setlength(h_resIgualdad, 0);

     XOpt := TVectR.Create_Init(ivar-1);

end;

procedure TProblema_onl.Free;
begin
   {$IFDEF SIMPLEX_HIJO}
       Self.spx.Free;
   {$ENDIF}

   x_inf.Free;
   x_sup.Free;

   f_costo.Free;
   setlength(g_resDesigualdad,0);
   setlength(h_resIgualdad,0);

   XOpt.Free;

   inherited Free;
end;

procedure TProblema_onl.FijarRestriccionIgualdad( kfila: NInt );
begin
     {$IFDEF SIMPLEX_HIJO}
         Self.spx.FijarRestriccionIgualdad(kfila);
     {$ENDIF}

     //paso una restricción que asumí era de desigualdad a una restricción de igualdad
     //esto es opcional no se si vale la pena programarlo para el ONL


end;

procedure TProblema_onl.FijarVariable( ivar: NInt; valor: NReal );
begin
     {$IFDEF SIMPLEX_HIJO}
         Self.spx.FijarVariable(ivar, valor);
     {$ENDIF}

     cota_inf_set( ivar, valor );
     cota_sup_set( ivar, valor );
end;

procedure TProblema_onl.cota_inf_set( ivar: NInt; vxinf: NReal );
begin
     {$IFDEF SIMPLEX_HIJO}
         Self.spx.cota_inf_set(ivar, vxinf);
     {$ENDIF}

     //rellenar algun vector de flags que despues se le pasa al optimizador
     //*sin terminar*
end;

procedure TProblema_onl.cota_sup_set( ivar: NInt; vxsup: NReal );
begin
    // {$IFDEF SIMPLEX_HIJO}
     //    Self.spx.cota_sup_set(ivar, vxinf);
     //{$ENDIF}

     //rellenar algun vector de flags que despues se le pasa al optimizador
     //*sin terminar*
end;

function TProblema_onl.e(k, j: integer): NReal;
begin
     result := Self.spx.e(k, j);
end;

procedure TProblema_onl.pon_e(k, j: integer; x: NReal);
begin
     Self.spx.pon_e(k, k, x);
end;

procedure TProblema_onl.acum_e(k, j: integer; x: NReal);
begin
     Self.spx.acum_e(k, j, x);
end;

function TProblema_onl.xval( ix: NInt ): NReal;
begin
    // {$IFDEF SIMPLEX_HIJO}
  //       Result := Self.spx.xval(ix);
  //       exit();
  //   {$ENDIF}

     Result := XOpt.e(ix);
end;

function TProblema_onl.yval( iy: NInt ): NReal;
begin
  //   {$IFDEF SIMPLEX_HIJO}
   //      Result := Self.spx.yval(iy);
  //       exit();
    // {$ENDIF}

     // ¿la restriccion es de desigualdad o de igualdad?
     if (iy <= nres_desigualdad) then begin
        Result := g_resDesigualdad[iy-1].Evaluar(XOpt);
     end else begin
        Result := h_resIgualdad[iy-1].Evaluar(XOpt);
     end;

end;

function TProblema_onl.xmult(ix: NInt) : NReal;
begin
  //   {$IFDEF SIMPLEX_HIJO}
   //      Result := Self.spx.xmult(ix);
   //      exit();
 //    {$ENDIF}
     //*sin terminar* para el onl
end;

function TProblema_onl.ymult(iy: NInt) : NReal;
begin
 //    {$IFDEF SIMPLEX_HIJO}
   //      Result := Self.spx.ymult(iy);
//         exit();
//     {$ENDIF}
     //*sin terminar* para el onl
end;

function TProblema_onl.fval : NReal;
begin
   //  {$IFDEF SIMPLEX_HIJO}
   //      Result := Self.spx.fval;
  //       exit();
  //   {$ENDIF}

     Result := Self.f_costo.Evaluar(XOpt);
end;

function TProblema_onl.Resolver: integer;
var
  i, j : NInt;
  vparams : TParams_onl;
  TolX,
  TolFun,
  alfa0,
  MaxIter : NReal;
  fcosto_spx: TFunc_CombinacionLineal;
  fc_coefs : TVectR;
  fci : NReal;
  cc : TDAOfNReal;  //coeficientes combinación lineal
  vidv : TDAOFNInt;
  idV : TVectR;
  rij : NReal;
  coef_res_ : TVectR;
  res_ : TFunc_CombinacionLineal;
  res_cajaTFR : TFunc_Recta;
  res_caja : TFunc_Polinomio;
  vc : TVectR;
  a : NReal;
  xxaux,xxaux1, xxaux2, xxaux3, xxaux4, xxaux5, xxaux6, xxaux7, xxaux8 : TVectR;
  xxaux1g1,xxaux1g2,
  xxaux2g1,xxaux2g2,
  xxaux3g1,xxaux3g2,
  xxaux4g1,xxaux4g2,
  xxaux5g1,xxaux5g2,
  xxaux6g1,xxaux6g2 :NReal;
  lxxaux,lxxaux1, lxxaux2, lxxaux3, lxxaux4, lxxaux5, lxxaux6, lxxaux7, lxxaux8 : NInt;
  xxr1: NReal;
  bi, acum : NReal;
  XX1, XX2, XX3, XX4, XX5, XX6, XX7, XX8 : NReal;

begin

   {$IFDEF RESOLVER_CON_SIMPLEX}
   if not (self.spx = nil) then begin
   //intento resolver
  if spx.resolver = 0 then
  begin
    //ok, encontró solución
    Writeln('Solucion optima encontrada:');
    //spx.fval obtiene el valor de z
    Writeln('z= ', FloatToStrF(-spx.fval, ffGeneral, 8, 4));
    Writeln;
    for i:= 1 to spx.nc-1 do
	//spx.xval(i) obtiene el valor de la variable i
      Writeln(#9, spx.fGetNombreVar(i), '= ', FloatToStrF(spx.xval(i), ffGeneral, 8, 3));
    Writeln;
    for i:= 1 to spx.nf-1 do
//spx.yval(i) obtiene el valor de la restriccion i
      Writeln(#9, spx.fGetNombreRes(i), '= ', FloatToStrF(spx.yval(i), ffGeneral, 8, 3));
  end
  else
    //Error, lanzamos la excepción
    raise Exception.Create('Error resolviendo simplex: ' + spx.mensajeDeError);

       exit();
       end;
   {$ENDIF}

   //Cargar la información del simplex hijo al problema de optimización no lineal
   //OJO: considero que las columnas del simplex se refieren a posiciones absolutas de la
   //descripción global del problema. . . de lo contrario los actores me deben dar esa
   //información
if not (self.spx = nil) then begin
   //asumo por ahora que el simplex es la única parte del problema de optimización..

   //1. Inscribo la función de costo
   setlength(cc,0);
   setlength(vidv, 0);
   lxxaux:=self.spx.nc-1; //cantidad de variables
   lxxaux2 := self.spx.nf-1; //cantidad de restricciones

   for j:= 1 to self.spx.nc-1 do begin
           fci := self.spx.e(self.spx.nf,j);
           if (abs(fci) > AsumaCero) then begin
              setlength(cc, length(cc)+1);
              cc[high(cc)]:=-fci;
              setlength(vidv, length(vidv)+1);
              vidv[high(vidv)]:=j;
           end;
   end;

   lxxaux:=length(vidv);
   lxxaux2:=length(cc);

   fc_coefs:=TVectR.Create_Init(length(cc));
   fc_coefs.pon_ev(1,cc);
   fcosto_spx := TFunc_CombinacionLineal.Create_Init(fc_coefs.n, fc_coefs);
   fcosto_spx.setIdVars(vidv);
   self.AgregarFObjetivo(fcosto_spx);

   //2. Recorro todas las restricciones: igualdad, desigualdad
   for i:=1 to self.spx.nf-1 do begin
       setlength(cc,0);
       setlength(vidv, 0);

       //generar una TFunc con la restricción
       for j:= 1 to self.spx.nc-1 do begin
           rij := self.spx.e(i,j);
           if (abs(rij) > AsumaCero) then begin
              setlength(cc, length(cc)+1);
              cc[high(cc)]:=rij;
              setlength(vidv, length(vidv)+1);
              vidv[high(vidv)]:=j;
           end;
       end;
       //la columna de los términos indep se trata aparte porque el usimplex
       //ya modificó dichos elementos en caso que se hayan seteado cotas inferiores

       bi := self.spx.e(i, spx.nc);
       acum:=0;
       for j:=1 to self.spx.nc-1 do begin
           acum:=acum+self.spx.e(i,j)*self.spx.x_inf.e(j);
       end;
       bi:=bi-acum;
       //solo agrego el término indep si es distinto de cero, de lo contrario CombinacionLineal ya se está
       //asumiendo con térm indep cero
       if (abs(bi) > AsumaCero) then begin
       coef_res_:=TVectR.Create_Init(length(cc));
       coef_res_.pon_ev(1, cc);
       res_:=TFunc_CombinacionLineal.Create_Init(coef_res_.n, coef_res_, bi);
       res_.setIdVars(vidv);
       end else begin
       coef_res_:=TVectR.Create_Init(length(cc));
       coef_res_.pon_ev(1, cc);
       res_:=TFunc_CombinacionLineal.Create_Init(coef_res_.n, coef_res_);
       res_.setIdVars(vidv);

       end;


       //de desigualdad: >= 0
       if (self.spx.flg_y[i] = 0) then begin
          self.AgregarRestriccionDeDesigualdad(res_)
       //de igualdad: = 0
       end else if (self.spx.flg_y[i] = 2) then begin
           self.AgregarRestriccionDeIgualdad(res_);
       end;


   end;



   //restricciones de caja: las agrego como restricciones de desigualdad
   for j:=1 to spx.nc-1 do begin

           //tipo : x <= a
           if (self.spx.flg_x[j] <> 0) then begin //hay cota superior
           a := self.spx.x_sup.e(j)+self.spx.x_inf.e(j); //el segundo término lo agrego porque
           //el simplex modifica la cota superior inicial restándole la cota inferior (no debería hacerlo antes de ponerse a calcular!)

           res_cajaTFR:=TFunc_Recta.Create_Init(-1,a);
           setlength(vidv,1);
           vidv[0]:=j;
           res_cajaTFR.setIdVars(vidv);
           self.AgregarRestriccionDeDesigualdad(res_cajaTFR);



           end;

           //tipo : a <= x
           a := self.spx.x_inf.e(j);

           res_cajaTFR:=TFunc_Recta.Create_Init(1,-a);
           setlength(vidv,1);
           vidv[0]:=j;
           res_cajaTFR.setIdVars(vidv);
           self.AgregarRestriccionDeDesigualdad(res_cajaTFR);


     end;


end;
    (*
   //PARTE PARA PROBAR SI LAS FUNCIONES ESTÁN BIEN INICIALIZADAS
   xxaux:=TVectR.Create_Init(2);
   xxaux.pon_e(1, 30); xxaux.pon_e(2, -10);
   XX1 := self.g_resDesigualdad[1-1].Evaluar(xxaux);

   lxxaux2:=length(self.h_resIgualdad);
   XX2 := self.h_resIgualdad[1-1].Evaluar(xxaux);

   XX3 := self.g_resDesigualdad[2-1].Evaluar(xxaux);
   XX4 := self.g_resDesigualdad[3-1].Evaluar(xxaux);
   XX5 := self.g_resDesigualdad[4-1].Evaluar(xxaux);
   XX6 := self.g_resDesigualdad[5-1].Evaluar(xxaux);

   //chequeo los gradientes
   xxaux1 := TVectR.Create_Clone(self.g_resDesigualdad[1-1].Gradiente(xxaux));
   lxxaux1:=xxaux1.n; xxaux1g1 := xxaux1.e(1); xxaux1g2 := xxaux1.e(2);

   xxaux2 := TVectR.Create_Clone(self.h_resIgualdad[1-1].Gradiente(xxaux));
   lxxaux2:=xxaux2.n;  xxaux2g1 := xxaux2.e(1); xxaux2g2 := xxaux2.e(2);

   xxaux3 := TVectR.Create_Clone(self.g_resDesigualdad[2-1].Gradiente(xxaux));
   lxxaux3:=xxaux3.n;  xxaux3g1 := xxaux3.e(1);

   xxaux4 := TVectR.Create_Clone(self.g_resDesigualdad[3-1].Gradiente(xxaux));
   lxxaux4:=xxaux4.n;  xxaux4g1 := xxaux4.e(1);

   xxaux5 := TVectR.Create_Clone(self.g_resDesigualdad[4-1].Gradiente(xxaux));
   lxxaux5:=xxaux5.n;  xxaux5g1 := xxaux5.e(1);

   xxaux6 := TVectR.Create_Clone(self.g_resDesigualdad[5-1].Gradiente(xxaux));
   lxxaux6:=xxaux6.n;  xxaux6g1 := xxaux6.e(1);
  *)
   //Escoger el método que utilizo para resolver el problema
   //despues se puede poner un campo seteado por el usuario que elige qué algoritmo
   //se utiliza y con qué parámetros...

     //el resolvedor_onl se encarga de tomar los parámetros y llamar a la función correspondiente


     //HACER EL DUMP DEL SIMPLEX : : : DESPUES BORRAR ESTO
     //self.spx.DumpSistemaToXLT('uex_mixed_spx_c2_DUMP.xlt','');


     vparams := TParams_onl.Create_Init(1,4);
     //punto inicial: el origen del SC
     vparams.Vec[0] := TVectR.Create_Init(nvars);
     vparams.Vec[0].Ceros;

     //ojo, esto no puede ser pasado asi nomas; o bien se le pide a cada actor que mande un punto inicial
     //dentro de su rango de validez de la variable, se ejecuta un algoritmo para arrancar con un punto factible

     //set TolX
     TolX := 1e-6;
     vparams.Num[0] := TolX;

     //set TolFun
     TolFun := 1e-6;
     vparams.Num[1] := TolFun;

     //set alfa0
     alfa0 := 1;
     vparams.Num[2] := alfa0;

     //set MaxIter
     MaxIter := 1e4;
     vparams.Num[3] := MaxIter;

     writeln('en udescripcionproblemaonl , high(g_resDesigualdad)=', high(g_resDesigualdad));
     writeln('en udescripcionproblemaonl , high(h_resIgualdad)=', high(h_resIgualdad));

     uresolvedor_onl.optimizar(f_costo,
                               g_resDesigualdad,
                               h_resIgualdad,
                               nvars,
                               nres_desigualdad,
                               nres_igualdad,
                               XOpt,
                               vparams);
   end;

procedure TProblema_onl.AgregarFObjetivo(f: TFunc);
begin
   Self.f_costo.Agregar_f(f);
end;

procedure TProblema_onl.AgregarRestriccionDeDesigualdad(g: TFunc);
begin
   //redimensionar
   setlength(self.g_resDesigualdad, cnt_res_desigualdad+1);
   self.g_resDesigualdad[cnt_res_desigualdad] := TFunc_Bloque.Create;
   self.g_resDesigualdad[cnt_res_desigualdad].Agregar_f(g);
   Inc(cnt_res_desigualdad);
end;

procedure TProblema_onl.AgregarRestriccionDeIgualdad(h: TFunc);
begin
   //redimensionar
   setlength(self.h_resIgualdad, cnt_res_igualdad+1);
   self.h_resIgualdad[cnt_res_igualdad] := TFunc_Bloque.Create;
   self.h_resIgualdad[cnt_res_igualdad].Agregar_f(h);
   Inc(cnt_res_igualdad);

end;


end.
