unit uMIPSimplexIteradorNoLineal;

interface

uses
  xmatdefs, usimplex, sysutils, TVectors,
{$IFDEF SPXMEJORCAMINO}
  umipsimplex_mejorcamino,
{$ELSE}
  umipsimplex
{$ENDIF};

type
  TMIPSimplexIteradorNoLineal = class(TMIPSimplex)
    private
      //Funcion externas para obtener nombres de restricciones
      //Si hay iteraciones el nombre de las restricciones será reemplazado por
      //extGetNombreRes(ires) + '_NroIter'
      extGetNombreRes: TFuncNombre;

      nMaxIteraciones: Integer;
      iterActual: Integer;
      //matrices[iter] contiene la matriz de la iteración iter
      //Si aun no se llego a la iteración iter matrices[iter] contiene NIL
      //Las filas de matrices[iter] se cargan a demanda, es decir sus valores son NIL
      //hasta que alguien ponga un valor en la fila. En caso de requerir un elemento
      //de una fila NIL se utilizara el elemento de la ultima iteración no NIL
      //Al llamar a resolver en la iteración iter se crea una única matriz con la
      //todas las filas no nulas en todas las matrices hasta iter
      //las matrices se indexan desde 1,1 como en MatR
      matrices: TDAOfMatOfNReal;
      //nFilasNoNulasIter[iter] contiene la cantidad de filas no nulas en matrices[iter]
      nFilasNoNulasIter: TDAofNInt;
      matrizActual: TMatOfNReal; //matrizActual = matrices[iterActual]

      //filas[fila_jres][filas_spx]
      //Filas contiene para la fila fila_jres (por la que preguntarán los actores)
      //cuales son las filas en el simplex que la representan con los sucesivos cortes
      filas: TDAOfTIntVector;

      spxActual: TMIPSimplex;

      function GetNombreResCIters(i: Integer) : String;      
    public
			constructor Create_init( mfilas, ncolumnas, nenteras, nMaxIteraciones: integer;
                               xfGetNombreVar, xfGetNombreRes : TFuncNombre); reintroduce; virtual;
      procedure Free;

      function e(k,j:integer):NReal; override;
      procedure pon_e(k,j: integer; x: NReal); override;
      procedure acum_e(k,j:integer; x: NReal); override;

      function xval( ix: integer ): NReal; override;
      function yval( iy: integer ): NReal; override;
      function xmult( ix: integer ): NReal; override;
      function ymult( iy: integer ): NReal; override;
      function fval: NReal; override;

      //Copia todos los datos que no sean la matriz DE self A spx.
      //Cajas, variables enteras, acoples y cualquier otra cosa que pueda surgir
      //y deba saberse al momento de resolver
      procedure copiarDatosNoMatriz(var spx: TMIPSimplex);

      function resolver: integer; override;
      procedure limpiar; override;

      procedure IniciarIters;
      procedure incIter;
  end;

procedure test();

implementation

constructor TMIPSimplexIteradorNoLineal.Create_init( mfilas, ncolumnas, nenteras, nMaxIteraciones: integer;
                                                     xfGetNombreVar, xfGetNombreRes : TFuncNombre);
var
  i: Integer;
begin
  inherited Create_init(mfilas, ncolumnas, nenteras, xfGetNombreVar, GetNombreResCIters);
  if not Assigned(xfGetNombreRes) then
    extGetNombreRes:= defaultGetNombreRes
  else
    extGetNombreRes:= fGetNombreRes;  

  self.nMaxIteraciones:= nMaxIteraciones;
  SetLength(matrices, nMaxIteraciones + 1);
  SetLength(nFilasNoNulasIter, nMaxIteraciones + 1);
  for i:= 0 to High(nFilasNoNulasIter) do
    nFilasNoNulasIter[i]:= 0;
  SetLength(filas, mfilas + 1);//Las columnas de jresOriginal se inicializan en cada resolver
  for i:= 0 to high(filas) do
    filas[i]:= TIntVector.Create(10);
end;

procedure TMIPSimplexIteradorNoLineal.Free;
var
  i: Integer;
  j: Integer;
begin
  if spxActual <> NIL then
    spxActual.Free;
  for i:= 0 to High(matrices) do
  begin
    for j:= 0 to high(matrices[i]) do
      SetLength(matrices[i][j], 0);
    SetLength(matrices[i], 0);
  end;
  SetLength(matrices, 0);
  SetLength(nFilasNoNulasIter, 0);
  for i:= 0 to high(filas) do
    filas[i].Free;
  SetLength(filas, 0);
  inherited Free;
end;

function TMIPSimplexIteradorNoLineal.e(k, j:integer):NReal;
var
  i_iter: Integer;
begin
  i_iter:= iterActual;
  while matrices[i_iter][k] = NIL do
  begin
    i_iter:= i_iter - 1;
    assert(i_iter >= 0, 'TMIPSimplexIteradorNoLineal.e: pidieron un elemento para una fila que nunca se cargó')        
  end;
  
  result:= matrices[i_iter][k][j];
end;

procedure TMIPSimplexIteradorNoLineal.pon_e(k, j: integer; x: NReal);
var
  i: Integer;
begin
  if matrizActual[k] = NIL then
  begin
    SetLength(matrizActual[k], nc + 1);
    for i:= 0 to high(matrizActual[k]) do
      matrizActual[k][i]:= 0;
    nFilasNoNulasIter[iterActual]:= nFilasNoNulasIter[iterActual] + 1;
  end;

  matrizActual[k][j]:= x;
end;

procedure TMIPSimplexIteradorNoLineal.acum_e(k, j:integer; x: NReal);
var
  i: Integer;
begin
  if matrizActual[k] = NIL then
  begin
    SetLength(matrizActual[k], nc + 1);
    for i:= 0 to high(matrizActual[k]) do
      matrizActual[k][i]:= 0;
    matrizActual[k][j]:= x;
    nFilasNoNulasIter[iterActual]:= nFilasNoNulasIter[iterActual] + 1;
  end
  else
    matrizActual[k][j]:= matrizActual[k][j] + x;
end;

function TMIPSimplexIteradorNoLineal.xval( ix: integer ): NReal;
begin
  result:= spxActual.xval(ix);
end;

function TMIPSimplexIteradorNoLineal.yval( iy: integer ): NReal;
var
  i: Integer;
  val, minYVal: NReal;
begin
  minYVal:= 0;
  for i:= 0 to filas[iy].size do
  begin
    val:= spxActual.yval(filas[iy][i]);
    if val = 0 then
      break
    else if val < minYVal then
      minYVal:= val;
  end;
  result:= minYVal;
end;

function TMIPSimplexIteradorNoLineal.xmult( ix: integer ): NReal;
begin
  result:= spxActual.xmult(ix);
end;

function TMIPSimplexIteradorNoLineal.ymult( iy: integer ): NReal;
var
  i: Integer;
  res: NReal;
begin
  res:= 0;
  for i:= 0 to filas[iy].size do
    res:= res + spxActual.ymult(filas[iy][i]);
  result:= res;
end;

function TMIPSimplexIteradorNoLineal.fval: NReal;
begin
  result:= spxActual.fval;
end;

procedure TMIPSimplexIteradorNoLineal.copiarDatosNoMatriz(var spx: TMIPSimplex);
var
  i: Integer;
begin
  for i:= 1 to nc - 1 do
  begin
    spx.cota_inf_set(i, self.x_inf.pv[i]);
    spx.cota_sup_set(i, self.x_inf.pv[i] + self.x_sup.pv[i]);
  end;
  for i:= 0 to nvents - 1 do
  begin
    spx.lstvents[i]:= lstvents[i];
    spx.lstAcoplesVEnts[i]:= copy(self.lstAcoplesVEnts[i], 0, length(self.lstAcoplesVEnts[i]));
  end;
end;

function TMIPSimplexIteradorNoLineal.resolver: integer;
var
  i_iter, j_fila, jfila_spx: Integer;
  nFilasIter: Integer;
begin
  nFilasIter:= 0;
  for i_iter:= 0 to iterActual do
    nFilasIter:= nFilasIter + nFilasNoNulasIter[i_iter];

  //TODO cambiar las funciones para obtener los nombres
  if spxActual <> NIL then
    spxActual.Free;
  spxActual:= TMIPSimplex.Create_init(nFilasIter, nc, nvents, fGetNombreVar, fGetNombreRes);;

  jfila_spx:= 1;
  for j_fila:= 1 to nf do
  begin
    filas[j_fila].Clear;
    for i_iter:= 0 to iterActual do
    begin
      //TODO en el else se podría hacer un break si se sabe que una vez que
      //no se agrega una restricción en una iteración ya no se va a agregar
      if matrices[i_iter][j_fila] <> NIL then
      begin
        spxActual.pon_fila(jfila_spx, 0, matrices[i_iter][j_fila]);
        filas[j_fila].add(jfila_spx);
        jfila_spx:= jfila_spx + 1;
      end;
    end;
  end;

  copiarDatosNoMatriz(spxActual);

//spxActual.DumpSistemaToXLT('debugIterador_iter' + IntToStr(iterActual) + '.xlt', '');
  result:= spxActual.resolver;
end;

procedure TMIPSimplexIteradorNoLineal.limpiar;
var
  i, j: Integer;
begin
  inherited limpiar;
  for i:= 0 to High(matrices) do
    for j:= 0 to high(matrices[i]) do
      matrices[i][j]:= NIL;
  for i:= 0 to high(nFilasNoNulasIter) do
    nFilasNoNulasIter[i]:= 0;

  matrizActual:= NIL;
  if spxActual <> NIL then
    spxActual.Free;
end;

procedure TMIPSimplexIteradorNoLineal.IniciarIters;
var
  i: Integer;
begin
  iterActual:= 0;
  SetLength(matrices[0], nf + 1);
  matrizActual:= matrices[0];
  for i:= 0 to high(matrizActual) do
    matrizActual[i]:= NIL;
end;

procedure TMIPSimplexIteradorNoLineal.incIter;
var
  i: Integer;
begin
  iterActual:= iterActual + 1;
  SetLength(matrices[iterActual], nf + 1);
  matrizActual:= matrices[iterActual];
  for i:= 0 to high(matrizActual) do
    matrizActual[i]:= NIL;
end;

function TMIPSimplexIteradorNoLineal.GetNombreResCIters(i: Integer) : String;
var
  iFila, iIter, k, j: Integer;
begin
  iFila:= -1;
  iIter:= -1;
  //Busco la fila original, la que indico el actor
  for k:= 1 to high(filas) do
  begin
    for j:= 0 to filas[k].size - 1 do
      if filas[k][j] = i then
      begin
        iFila:= k;
        if filas[k].size <> 1 then  
          iIter:= j;
        break;
      end;
    if (iFila <> -1) then
      break;
  end;

  if iIter <> -1 then
    result:= extGetNombreRes(iFila) + '_' + IntToStr(iIter + 1)
  else
    result:= extGetNombreRes(iFila);
end;

procedure test();
var
  iIteracion: Integer;
  res, xi, xiAnterior, yi: NReal;
  spx: TMIPSimplexIteradorNoLineal;
  zxi, dzxi: NReal;
begin
  //max 8x1 + 19x2 + 7x3
  //s.a.
  //3x1 + 4x2 + x3 <= 25  ==>   -3x1 -4x2 -x3 +25 >= 0
  //x1 + x2 + 3x3 <= 50   ==>   -x1 -x2 -3x3 +50 >= 0
  //x1, x2 >= 0 reales, 16 >= x3 >= 0 entera
  spx:= TMIPSimplexIteradorNoLineal.Create_init(3, 4, 1, 5, NIL, NIL);
  spx.IniciarIters;

  spx.pon_e(1, 1, -3);
  spx.pon_e(1, 2, -4);
  spx.pon_e(1, 3, -1);
  spx.pon_e(1, 4, 25);

  spx.pon_e(2, 1, -1);
  spx.pon_e(2, 2, -1);
  spx.pon_e(2, 3, -3);
  spx.pon_e(2, 4, 50);

  spx.pon_e(spx.nf, 1, 8);
  spx.pon_e(spx.nf, 2, 19);
  spx.pon_e(spx.nf, 3, 7);

  spx.cota_inf_set(1, 0);
  spx.cota_inf_set(2, 0);
  spx.cota_inf_set(3, 0);
  spx.cota_sup_set(1, MaxNReal);
  spx.cota_sup_set(2, MaxNReal);
  spx.cota_sup_set(3, MaxNReal);
  spx.set_entera(1, 2, 16);

  spx.resolver;
//  res:= spx.fval;
  spx.Free;

  //max -x^2 = z
  //s.a.
  // -2 <= x <= 2 real

  //Para resolverlo resolveremos la siguiente serie de problemas
  //      {max y
  //P(i)  {s.a.
  //      {y <= z(x0) + z'(x0)*(x-x0) ==> -y +x * z'(x0) -z'(x0)*xi+z(x0) >= 0
  //      {...
  //      {y <= z(xi) + z'(xi)*(x-xi) ==> -y +x * z'(xi) -z'(xi)*xi+z(xi) >= 0
  //      {-2 <= x <= 2 real

  //Donde x e y son variables de control, y es limitada por las restricciones,
  //z' es la derivada de z y xi es el valor resultado de x en la iteración
  //anterior

  spx:= TMIPSimplexIteradorNoLineal.Create_init(2, 3, 0, 20, NIL, NIL);
  spx.IniciarIters;
  xiAnterior:= -MaxNReal;
  xi:= -2;

  //Cargamos los datos (filas, cotas, variables enteras y acoples) que no participan
  //en la iteración
  spx.pon_e(2, 1, 1);

  spx.cota_inf_set(1, -MaxInt);
  spx.cota_sup_set(1, MaxInt);
  spx.cota_inf_set(2, -2);
  spx.cota_sup_set(2, 2);  

  iIteracion:= 0;
  while (iIteracion < spx.nMaxIteraciones) and (abs(xi - xiAnterior) > 0.001) do
  begin
    zxi:= -sqr(xi);  //z(xi) = -(xi^2)
    dzxi:= -2 * xi;  //z'(xi) = -2 * xi

    //y será la variable 1 y x la variable 2
    spx.pon_e(1, 1, -1);
    spx.pon_e(1, 2, dzxi);
    spx.pon_e(1, 3, -dzxi * xi + zxi);

    spx.resolver;
    res:= spx.fval;
    xiAnterior:= xi;
    yi:= spx.yval(1);
    xi:= spx.xval(2);
    writeln('Iter' + IntToStr(iIteracion) + ' xi=' + FloatToStrF(xi, fffixed, 10, 3) + ' yi=' + FloatToStrF(yi, fffixed, 10, 3) + ' fval= ' + FloatToStrF(res, fffixed, 10, 3));
    spx.incIter;
    iIteracion:= iIteracion + 1;
  end;

//  res:= spx.fval;
  spx.Free;
end;

end.
