{ optimización con restricciones }
unit miner5;
interface
uses
	xMatDefs, miner, MatReal;

procedure Init( NVariables, NRestricciones: integer );
procedure DefFunc(
	irestric: integer;
	frestric: TFunc_fx;    { dominio= frestric(X) < 0 }
	grestric: TGrad_fx );

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

const
	MaxNRestric= 1000;

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

	PProblema= ^TProblema;
	TProblema= array[0..MaxNRestric] of TFichaFunc;

var
	ProblemaPtr: PProblema;
	rel: NReal;
	NVs, NRs: integer;

//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; far;
var
	k: integer;
	m, w: NReal;
begin
	m:= 0;
	for k:= 1 to NRs do
	begin
		w:= ProblemaPtr^[k].f(k, x);
		ProblemaPtr^[k].uv:= w;
		if w > 0 then m:= m + w;
	end;
	m:=	ProblemaPtr^[0].f(0, x) + rel*m;
	frel:= m;
end;

{Acumula los gradientes de las restricciones violadas en mg, y
se lo suma al gradiente de la funcion objetivo ponderado por
max(rel, (|df||mg| - df.mg)/(|mg|^2))
}
procedure grel(kf: integer; var grad, X: TVectR ); far;
var
	k: integer;
	mg: TVectR;
	TempRel: NReal;
	norm_r, norm_f: NReal;
begin
	mg.init( X.N );
	mg.Ceros;
	for k:= 1 to NRs do
	begin
		ProblemaPtr^[k].g(k, grad, x);
		if ProblemaPtr^[k].uv > 0 then
			mg.Sum( grad );
	end;
	ProblemaPtr^[0].g(0, grad, X);
	norm_r:= mg.NormEuclid;
	norm_f:= grad.NormEuclid;
	if norm_r > AsumaCero then
	begin
	  TempRel:= (norm_f*norm_r - grad.PEV(mg))/ (norm_r*norm_r);
	  if TempRel > rel then rel:= TempRel;
	end;
	grad.SumRPV( rel, mg );
	mg.done;
end;

procedure Init( NVariables, NRestricciones: integer );
begin
	NVs:= NVariables;
	NRs:= NRestricciones;
	GetMem( ProblemaPtr, SizeOF(TFichaFunc)*(NRs+1) );
end;

procedure DefFunc(
	irestric: integer;
	frestric: TFunc_fx;
	grestric: TGrad_fx );
begin
	ProblemaPtr^[irestric].f:= frestric;
	ProblemaPtr^[irestric].g:= grestric;
	ProblemaPtr^[irestric].uv:= 0;
end;

procedure Done;
begin
	FreeMem( ProblemaPtr, SizeOF(TFichaFunc)*(NRs+1) );
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;
	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
	relax: NReal;
begin
	relax:= 1;

		MinRelax(
				relax,
				X,
				pasoini,
				errs,
				MaxNIters,
				NIters,
				ValCosto,
				Convergio);

	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 ');
end;

end.

