unit uBatMiner2;

{
La idea es poder optimizar una función sin necesitar conocer su gradiente.
Para estimar el valor del gradiente el programa avanza o retrocede
cada variable un paso y calcula el valor de f, la dirección del gradiente
se fija hacia donde la función mejore.
Como la optimizacion se hace variable a variable se precisa que la funcion
relajada sea continua en X.
Si f es costosa de calcular este metodo no resulta eficiente puesto que
calcula el valor de f NVariables*3 veces por calculo de gradiente. 
}

interface

uses
  miner, xMatDefs, MatReal;

procedure Init(NVariables, NRestricciones: integer);
procedure DefFunc(irestric: integer;
            	    frestric: TFunc_fx    { dominio= frestric(X) < 0 });
procedure MinRelax( relax: NReal;
                  	var X: TVectR;
                  	paso: NReal;
                  	errs: NReal;
                	  MaxNIters: longint;
                  	var NIters: longint;
                  	var ValCosto: NReal;
                  	var Convergio: boolean );
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 );
procedure Done;

implementation

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

var
  globPaso, 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;
	frel:= m;
end;

{Para cada variable Xi del problema prueba los valores de frel que se
se obtienen dejando Xi igual, aumentandola paso y disminuyendola paso.
Como grad[i] devuelve 0 si el mejor valor de frel es sin cambiar Xi,
paso si el mejor valor de frel es con Xi = Xi + paso y -paso si el
mejor valor de frel es con Xi = Xi - paso.
NOTA: para calcular el valor de frel para Xi+1 asume que X1..Xi seran
las que obtuvieron el mejor valor de frel hasta ahora
}
procedure grel(kf: integer; var grad, X: TVectR );
var
  xClone: TVectR;
  i: Integer;
  mejorf0, f0, mejorXi: NReal;
begin
  xClone:= TVectR.Create_Clone(X);
	mejorf0:= MaxNReal;
  mejorXi:= 0;
  for i:= 1 to NVs do
  begin
    f0:= frel(0, xClone);
    //Dejar Xi como esta
    if f0 < mejorf0 then
    begin
      mejorf0:= f0;
      grad.pon_e(i, 0);
      mejorXi:= xClone.e(i);
    end;
    //Xi = Xi + paso
    xClone.acum_e(i, globPaso);
    f0:= frel(0, xClone);
    if f0 < mejorf0 then
    begin
      mejorf0:= f0;
      grad.pon_e(i, -globPaso); //se pone - porque se toma la direccion opuesta a grad
      mejorXi:= xClone.e(i);
    end;
    //Xi = Xi - paso
    xClone.acum_e(i, -2*globPaso);
    f0:= frel(0, xClone);
    if f0 < mejorf0 then
    begin
      mejorf0:= f0;
      grad.pon_e(i, globPaso);  //se pone + porque se toma la direccion opuesta a grad
      mejorXi:= xClone.e(i);
    end;
    xClone.pon_e(i, mejorXi);//el menos es porque se toma la direccion
                              //opuesta del gradiente
  end;
end;

procedure MinRelax(relax: NReal;
                   var X: TVectR;
                   paso: NReal;
                   errs: NReal;
                   MaxNIters: longint;
                   var NIters: longint;
                   var ValCosto: NReal;
                   var Convergio: boolean );
begin
	rel:= relax;
  globPaso:= paso;
	MinerCPS(	frel, grel, X, paso, errs, MaxNITers, NIters, ValCosto, Convergio );
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 );
var
	krel: integer;
	relax: NReal;
	xa: TVectR;
	errsrel, pasorel: NReal;
	dx: NReal;
begin
	xa:= TVectR.Create_Clone(X);
	nitersrel:= 0;
	krel:= 0;
	repeat
		inc(krel);
		inc(nitersrel);
		relax:= krel*sqrt(krel);
		pasorel:= {dx;} pasoini/sqrt(relax);
		if pasorel < errs then pasorel:= 1.1*errs;
		errsrel:= pasorel/10;
		if errsrel < errs then errsrel:= errs;
		MinRelax(
				relax,
				X,
				pasorel,
				errsrel,
				MaxNIters,
				NIters,
				ValCosto,
				Convergio);
		Xa.SumRPV( -1, X);
		dx:= Xa.NormEuclid;
		Xa.copy( X );
	{
	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 ');	 }
	until (dx < errs ) or ( nitersrel > MaxNItersRel ) or not Convergio;
	xa.Free;
end;

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

end.
 