{
	Procedimientos de minimizaci¢n basados en b£squeada en la
	direcci¢n del gradiente.
}

unit miner;
interface
uses
	xMatDefs, MatReal;

type
	TFunc_fx= function(kf: integer; var X: TVectR ): NReal;
	TGrad_fx= procedure(kf: integer; var grad, X: TVectR);

const
{ Control de Paso en X }
	CPXMult: NReal= 3;
	CPXDiv: NReal= 4.2;

{ Control de Paso en S (abcisa curvilinea) }
	CPSMult: NReal= 1.1;
	CPSDiv: NReal= 5;

{ control de paso por Filtro de Pasado Ponderado }
	FPPMult: NReal= 1.1;
	FPPDiv: NReal= 5;

procedure MinerCPX(
					fcosto: TFunc_fx;
					grad: TGrad_fx;
					var X: TVectR;
					paso: NReal;
					errx: NReal;
					MaxNIters: longint;
					var NIters: longint;
					var ValCosto: NReal;
					var Convergio: boolean );

procedure MinerCPS(
					fcosto: TFunc_fx;
					grad: TGrad_fx;
					var X: TVectR;
					paso: NReal;
					errs: NReal;
					MaxNIters: longint;
					var NIters: longint;
					var ValCosto: NReal;
					var Convergio: boolean );


procedure MinerFPP(
					fcosto: TFunc_fx;
					grad: TGrad_fx;
					var X: TVectR;
					paso: NReal;
					errx: NReal;
					MaxNIters: longint;
					var NIters: longint;
					var ValCosto: NReal;
					var Convergio: boolean;
					FiltCoefB: NReal ); { Coeficiente de memoria }

procedure vswap( var x, y: TVectR );

implementation


procedure vswap( var x, y: TVectR );
var
	tmp: ^TDAofNReal;
begin
	tmp:= @x.pv;
	x.pv:= y.pv;
	y.pv:= tmp^;
end;

procedure Versor( var X: TVectR);
var
	nor: NReal;
begin
	nor:= X.NormEuclid;
	if  EsCero(nor) then
	begin
		X.Ceros;
		X.pon_e(1,1);
	end
	else
		X.PorReal(1/nor);
end;

{Minero Con Control del Paso en X
Realiza la busqueda de un minimo en una funcion de Rn dando
los pasos en la direccion del gradiente y con sentido opuesto
(para minimizar).
Inicialmente fija un valor de pasos utiles para actualizar el
paso, arbitrariamente.
Lleva un conteo de los pasos utiles seguidos (disminuyo el
valor de fcosto) que dio, si este valor supera al de los
pasos de actualizacion aumenta el tamaño del paso por CPXMult
y disminuye la cantidad de pasos requeridos para una
actualizacion. De esta forma si vamos bien el proceso se acelera.
A su vez cuenta la cantidad de fallas seguidas (no disminuyo
el valor de fcosto) que dio, si superan a los pasos de
actualizacion reduce el tamaño del paso sobre CPXDiv y aumenta
la cantidad de pasos para una actualizacion. De esta forma
si empeoramos se hara una busqueda con pasos mas cortos para no
"mejorar" con el valor anterior que teniamos.
}
procedure MinerCPX(
					fcosto: TFunc_fx;
					grad: TGrad_fx;
					var X: TVectR;
					paso: NReal;
					errx: NReal;
					MaxNIters: longint;
					var NIters: longint;
					var ValCosto: NReal;
					var Convergio: boolean );
var
	PasoOK: boolean;
	fo, fs: NReal;
	v, Xs: TVectR;
	OksSeguidos, FallasSeguidas, OksAct: longint;
begin
	v:= TVectR.Create_Init(X.n);
	Xs:= TVectR.Create_Init(X.n);
	Convergio:= false;
	PasoOK:= false;
	fo:= fcosto(0, X);
	grad(0, v, X);

	versor(v);
	OksSeguidos:= 0;
	FallasSeguidas:= 0;
	OksAct:= 1;
	NIters:= 0;
	repeat
		Xs.copy(X);
		Xs.SumRPV( -paso, v);
		fs:= fcosto(0, Xs);
			Inc( NIters );
		if fs <= fo then
		begin
			Inc( OksSeguidos );
			FallasSeguidas:= 0;
			vswap(X, Xs);
			grad(0, v, X);
			versor(v);
			fo:= fs;
			if OksAct <= OksSeguidos then
			begin
				paso:= paso*CPXMult;
				dec( OksAct );
				if OksAct= 0 then OksAct:= 1;
			end;
			pasoOk:= true;
		end
		else
		begin
			inc( OksAct );
			inc( FallasSeguidas );
      OksSeguidos:=0 ;
			paso:= paso/CPXDiv;
			if paso < errx then
				Convergio:= true;
		end;
{		writeln(Niters: 6,'  ', fo: 12: 4,' ', paso, ' >',OksAct); }
	until Convergio or(NIters >= MaxNIters );
	ValCosto:= fo;
	v.Free;
	Xs.Free;
end;

{Minero con Control de Paso sobre S
Hace lo mismo que el Minero con Control de Paso sobre X pero
el tamaño del paso es sobre la abscisa curvilinea.
De esta forma se amortigua el efecto sobre el paso en X de un
gradiente muy empinado ya que parte del paso se consume en
delta fcosto
}
procedure MinerCPS(
					fcosto: TFunc_fx;
					grad: TGrad_fx;
					var X: TVectR;
					paso: NReal;
					errs: NReal;
					MaxNIters: longint;
					var NIters: longint;
					var ValCosto: NReal;
					var Convergio: boolean );
var
	PasoOK: boolean;
	fo, fs: NReal;
	v, Xs: TVectR;
	OksSeguidos, FallasSeguidas, OksAct: longint;
	ngf2: NReal;
begin
	v:= TVectR.Create_Init(X.n);
	Xs:= TVectR.Create_Init(X.n);
	Convergio:= false;
	PasoOK:= false;
	fo:= fcosto(0, X);
	grad(0, v, X);
	ngf2:= v.ne2;
	versor(v);
	OksSeguidos:= 0;
	FallasSeguidas:= 0;
	OksAct:= 1;
	NIters:= 0;
	repeat
		Xs.copy(X);
		Xs.SumRPV( -paso/sqrt(1+ngf2), v);
		fs:= fcosto(0, Xs);
    Inc( NIters );
		if fs <= fo then
		begin
			Inc( OksSeguidos );
			FallasSeguidas:= 0;
			vswap(X, Xs);
			grad(0, v, X);
			ngf2:= v.ne2;
			versor(v);
			fo:= fs;
			if OksAct <= OksSeguidos then
			begin
				paso:= paso*CPSMult;
				dec( OksAct );
				if OksAct= 0 then OksAct:= 1;
      end;
			pasoOk:= true;
		end
		else
		begin
			inc( OksAct );
			inc( FallasSeguidas );
			OksSeguidos:= 0 ;
			paso:= paso/CPSDiv;
			if paso < errs then
				Convergio:= true;
		end;
//		writeln(Niters: 6,'  ', fo: 12: 4,' ', paso, ' >',OksAct);
	until Convergio or(NIters >= MaxNIters );
  ValCosto:= fo;
  v.Free;
  Xs.Free;
end;



procedure MinerFPP(
					fcosto: TFunc_fx;
					grad: TGrad_fx;
					var X: TVectR;
					paso: NReal;
					errx: NReal;
					MaxNIters: longint;
					var NIters: longint;
					var ValCosto: NReal;
					var Convergio: boolean;
					FiltCoefB: NReal ); { Coeficiente de memoria }

var
	fo, fs: NReal;
	v, vfilt, Xs: TVectR;
	OksSeguidos, FallasSeguidas, OksAct: longint;
	ngf2: NReal;
	FiltCoefA: NReal;

begin
	FiltCoefA:= 1-FiltCoefB;
	v:= TVectR.Create_Init(X.n);
	vfilt:= TVectR.Create_Init(X.n);
	vfilt.ceros;
	Xs:= TVectR.Create_Init(X.n);
	Convergio:= false;
	fo:= fcosto(0, X);
	OksSeguidos:= 0;
	FallasSeguidas:= 0;
	OksAct:= 1;
	NIters:= 0;
	repeat
		grad(0,  v, X);
		vfilt.PorReal( FiltCoefB );
		vfilt.SumRPV( FiltCoefA, v );
		v.Copy( vfilt );
		ngf2:= v.ne2;
		versor(v);
		Xs.copy(X);
		Xs.SumRPV( -paso/sqrt(1+ngf2), v);
		fs:= fcosto(0, Xs);
		Inc( NIters );
		if fs <= fo then
		begin
			Inc( OksSeguidos );
			FallasSeguidas:= 0;
			if OksAct<= OksSeguidos then
			begin
				paso:= paso*FPPMult;
				dec( OksAct );
				if OksAct= 0 then OksAct:= 1;
			end;
		end
		else
		begin
			inc( OksAct );
			inc( FallasSeguidas );
				 OksSeguidos:=0 ;
			paso:= paso/FPPDiv;
			if paso < errx then
				Convergio:= true;
		end;
		vswap(X, Xs);
		fo:= fs;
		{		writeln(Niters: 6,'  ', fo: 12: 4,' ', paso, ' >',OksAct); }
	until Convergio or(NIters >= MaxNIters );
	ValCosto:= fo;
	vfilt.Free;
	v.Free;
	Xs.Free;
end;



end.

