unit uMapMakingMiner;

interface

{
Se trata de hacer un minero que no conozca el gradiente de la funcion
a optimizar, ni de sus restricciones pero ademas evaluar la funcion es
costoso. El minero guarda los valores por los que va pasando y calcula
un "gradiente promedio" como la suma de la aproximación del gradiente
por (f(X) - f(Xi))/dist(X, Xi) ponderado por 

}

uses
  MatReal, xMatDefs, Miner;

//xmaxDistanciaMapa es la maxima distancia entre X y Xi para la cual se
//considerara el valor guardado de f(Xi)
procedure Init(NVariables, NRestricciones: integer; xmaxDistanciaMapa: NReal);
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;

implementation

uses Math;

type
  TNodoValorFuncional = class
    public
      X: TVectR;
      FdeX: NReal;
      sig: TNodoValorFuncional;
      Constructor Create(X: TVectR; FdeX: NReal; siguiente: TNodoValorFuncional);
      procedure Free;
  end;

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

var
  mapa: TNodoValorFuncional;
  rel, maxDistanciaMapa: NReal;
  NVs, NRs, NPuntosMapa: Integer;
  problema: array of TFichaFunc;

Constructor TNodoValorFuncional.Create(X: TVectR; FdeX: NReal; siguiente: TNodoValorFuncional);
begin
  inherited Create;
  self.X:= TVectR.Create_Clone(X);
  self.FdeX:= FdeX;
  Self.sig:= siguiente;
end;

procedure TNodoValorFuncional.Free;
begin
  if sig <> NIL then
    sig.Free;
  X.Free;
  inherited Free;
end;

procedure Init(NVariables, NRestricciones: integer; xmaxDistanciaMapa: NReal);
begin
  mapa:= NIL;
	NVs:= NVariables;
	NRs:= NRestricciones;
  NPuntosMapa:= 0;
  SetLength(problema, NRs + 1);
  maxDistanciaMapa:= xmaxDistanciaMapa;
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;
  //se crea el nodo nuevo y se agrega el mapa viejo como cola
  //de la lista
  mapa:= TNodoValorFuncional.Create(X, m, mapa);
  Inc(NPuntosMapa);
end;

{Para cada variable Xi sortea un valor entre -1 y 1 y se lo asigna
}
procedure grel(kf: integer; var grad, X: TVectR );
var
  i: Integer;
  sumAlfas, dist, alfa: NReal;
  iter: TNodoValorFuncional;
  FdeX: NReal;
begin
  grad.Ceros;
  if mapa <> NIL then
  begin
    iter:= mapa.sig; //salteo el primero porque fue el valor que calcule recien
    sumAlfas:= 0;
    FdeX:= problema[0].uv;
    while iter <> NIL do
    begin
      dist:= X.distancia(iter.X);
      if dist < maxDistanciaMapa then
      begin
        alfa:= (FdeX - iter.FdeX)/dist;
        sumAlfas:= sumAlfas + alfa;
        for i:= 1 to NVs do
          grad.acum_e(i, (X.e(i) - mapa.X.e(i)) * alfa);
      end;
      iter:= iter.sig;
    end;
    grad.PorReal(1/sumAlfas);
  end;
  //sorteos para explorar
  for i:= 1 to NVs do
    grad.acum_e(i, (Random*2 - 1) * 1/sqrt(NPuntosMapa) );
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;
    grel(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
    grel(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);
  if mapa <> NIL then
    mapa.Free;
end;


end.
