unit uRandMiner2;

interface

{
La idea es poder optimizar una función sin necesitar conocer su gradiente.
Se trata de no asumir nada sobre la funcion y evitar tener que calcular
f repetidas veces por lo que el programa sortea una direccion al azar para
usar como gradiente.
Se larga desde X NExploraciones veces y en cada iteracion conserva la mejor
solucion que haya encontrado, luego se larga desde la mejor solución con
pasos cada vez mas pequeños intentando buscar en la localidad de la mejor
solucion algun X que mejore.
Puede ser util usar este metodo para encontrar algun valor cercano al
optimo y luego usar uBatMiner para afinar la solucion.
}

uses
  miner, xMatDefs, MatReal, Math;

procedure Init(NVariables, NRestricciones: integer);
procedure DefFunc(irestric: integer;
            	    frestric: TFunc_fx    { dominio= frestric(X) < 0 });
procedure Minimizar(var X: TVectR;
                  	pasoini: NReal;
                	  errs: NReal;
                  	MaxNItersRel: longint;
                  	var NItersRel: longint;
                  	MaxNIters: longint;
                	  var NIters: longint;
                  	var ValCosto: NReal;
                  	var Convergio: boolean;
                    NExploraciones: Integer );
procedure Done;
procedure printXFX(X: TVectR; fx: NReal);

implementation

type
	TFichaFunc= record
		uv: NReal;
		f: TFunc_fx;
	end;

var
  rel: NReal;
  NVs, NRs: Integer;
  problema: array of TFichaFunc;

procedure Init(NVariables, NRestricciones: integer);
begin
	NVs:= NVariables;
	NRs:= NRestricciones;
  SetLength(problema, NRs + 1);
end;

procedure DefFunc(irestric: integer;
             	    frestric: TFunc_fx{ dominio= frestric(X) < 0 });
begin
  problema[irestric].f:= frestric;
  problema[irestric].uv:= 0;
end;

//Acumula el valor de las restricciones violadas y se lo suma
//ponderado por rel al valor de la funcion objetivo
function frel(kf: integer; var X: TVectR ): NReal;
var
	k: integer;
	m, w: NReal;
begin
	m:= 0;
	for k:= 1 to NRs do
	begin
		w:= problema[k].f(k, x);
		problema[k].uv:= w;
		if w > 0 then m:= m + w;
	end;
	m:=	problema[0].f(0, x) + m*rel;
//  printXFX(X, m);
  result:= m;
end;

{Para cada variable Xi sortea un valor entre -1 y 1 y se lo asigna
}
procedure grelExploracion(kf: integer; var grad, X: TVectR );
var
  i: Integer;
begin
  for i:= 1 to NVs do
    grad.pon_e(i, RandomRange(1, -1));
end;

procedure grelAfinacion(kf: integer; var grad, X: TVectR );
var
  i: Integer;
begin
  for i:= 1 to NVs do
    grad.pon_e(i, Random*2 - 1 );
end;

procedure Minimizar(var X: TVectR;
                    pasoini: NReal;
                    errs: NReal;
                    MaxNItersRel: longint;
                    var NItersRel: longint;
                    MaxNIters: longint;
                    var NIters: longint;
                    var ValCosto: NReal;
                    var Convergio: boolean;
                    NExploraciones: Integer );
var
  i: Integer;
  xClone, g: TVectR;
  mejorf0: NReal;
  paso: NReal;
begin
  rel:= 10;
  g:= TVectR.Create_Init(X.n);
  xClone:= TVectR.Create_Clone(X);
  mejorf0:= MaxNReal;
  NIters:= 0;
  for i:= 1 to NExploraciones do
  begin
    ValCosto:= frel(0, xClone);
    if ValCosto < mejorf0 then
    begin
      X.Copy(xClone);
      mejorf0:= ValCosto;
    end;
    grelExploracion(0, g, X);
    xClone.sumRPV(pasoini, g);
  end;

  xClone.Copy(X);
  ValCosto:= mejorf0;
  paso:= pasoini/2;
  NItersRel:= 1;
  while (paso > errs) and (NItersRel <= MaxNItersRel) do
  begin
    grelAfinacion(0, g, X);
    xClone.sumRPV(pasoini, g);
    ValCosto:= frel(0, xClone);
    if ValCosto < mejorf0 then
    begin
      X.Copy(xClone);
      mejorf0:= ValCosto;
      paso:= paso / 2;
      NItersRel:= 1;
    end
    else
    begin
      Inc(NItersRel);
      //vuelvo a como estaba antes para mantenerme en la localidad
      //de la mejor solucion encontrada
      xClone.Copy(X);
    end
  end;

  ValCosto:= mejorf0;
  Convergio:= paso <= errs;
	{
	write('Rx: ', relax:12,'  NI: ', NIters:4);
	if convergio then
		writeln( '  x1: ',X.e(1):8:-4, '  x2: ', X.e(2):8:-4,' f(x): ', ValCosto:8:-4 )
	else
		writeln(' NO CONVERGIO ');	 }
	xClone.Free;
  g.Free;
end;

procedure Done;
begin
  SetLength(problema, 0);
end;

procedure printXFX(X: TVectR; fx: NReal);
begin
  X.Print;
  writeln('f(X)= ', fx:15:3);
end;

end.
