{

En el intento de tener un buen minimizador de problemas
con restricciones volvemos a la idea de sumar las restricciones
solamente fuera del dominio y con multiplicadores. Los multiplicadores
calculados de manera tal que resulte que el gradiente compuesto
es paralelo a la frontera del dominio, sobre el dominio y con una
componente en el sentido de llevar el punto a la frontera cuando nos
alejamos del dominio.

Problema:
 min r0(x)
 rk(x) < 0; con k = 1..nr

Problema relajado:
 min L = r0 + suma(ak* rk(x), k=1..nr)

-----------------------------------------
El gradiente de L es:
 g=  dL/dx = g0+ suma( ak* gk(x), k= 1..nr)

 donde ak > 0 si rk(x) >= 0 y 0 si rk(x) < 0,
 y gk= drk/dx con k= 0..nr

Si la soluci¢n del problema est  dentro del dominio, todos los
ak son CERO.
Si la solución est  en la frontera, los ak correspondientes a las
restricciones activas serán > 0 y el resto = 0.

1) Supongamos que en determinado punto de la iteración está activa
solamente la restricción (j) entonces,

g= g0 + aj * gj(x)

para que se cumpla lo requerido respecto del gradiente compuesto (g)
imponemos que: (g, gj(x)) = rj(x)  ( >0 al estar activa la restricción)

=> (g0, gj) + sum( ai * ( gi, gj )) = rj

=> aj =[(gi, gj)]^(-1). (rj - (g0, gj))

}
unit miner10;
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
	PFichaFunc= ^TFichaFun;
	TFichaFunc= record
		activa: boolean;
		rel: NReal;
		uv: NReal;
		ug: TVectR;
		f: TFunc_fx;
		g: TGrad_fx;
		nextActiva: integer; { índice de la siguiente activa }
	end;

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

var
	ProblemaPtr: PProblema;
	NRestriccionesActivas: integer;
	NVs, NRs: integer;

{ Computa las funciones, actualizando el campo (activa) para
las restricciones y guardando en (uv) el valor de la funci¢n.
Retorna la cantidad de restricciones activas }
function ComputarFunciones( var X: TVectR ): NReal; far;
var
	NActivas: integer;
	k: integer;
	w: NReal;
	pr: PFichaFunc;
	AntAct: integer;

begin
	NActivas:= 0;
	AntAct:= 0;

{ Primero recorremos las restricciones }
	for k:= 1 to NRs do
	begin
		pr:= @ProblemaPtr^[k];
		w:=pr^.f(k, x);
		pr^.uv:= w;
		if w > 0 then
		begin
			ProblemaPtr^[AntAct].NextActiva:= k; { encadenamos la proxima activa }
			AntAct:= k;
			inc(NActivas);
			pr^.Activa:= true;
		end
		else
			pr^.Activa:= false;
	end;

{ Terminamos con (-1) la lista de activas }
	ProblemaPtr^[AntAct].NextActiva:= -1;

{ Llamamos la funci¢n de costos del problema no relajado }
	pr:= @ProblemaPtr^[0];
	pr^.uv:=	pr^.f(0, x);
	NRestriccionesActivas:= NActivas;
	ComputarFunciones:= NActivas;
end;



{ Calculo los gradientes para la funci¢n de costo y para
las restricciones activas }
procedure ComputarGradientes( X: TVectR ); far;
var
	k: integer;
	norm_r, norm_f: NReal;
	pr: PFichaFunc;

begin
	{ barrido de las restricciones activas calculando sus
	gradientes }
	na:= ProblemaPtr^[0].NextActiva;
	while na > 0 do
	begin
		pr:= @ProblemaPtr^[na];
		pr^.g(k,  pr^.ug, x);
		na:= pr^.NextActiva;
	end;

	{ c lculo del gradiente de la funci¢n de costo original }
	ProblemaPtr^[0].g(0, ug, X);
end;



procedure CalcularMultiplicadores;

{=> ai =[(gi, gj)]^(-1). (rj - (g0, gj))}
var
	M: TMatR;
	a: TVectR;
	k, j: integer;

begin
	N:= NRestriccionesActivas;
	M.Init(N, N);
	a.Init(N);
	for k:= 1 to N do
	begin
		a.pon_e( r( k ) - gkgj( 0, k )
		m.pon_e(k, k, gkgj(k,k));
		for j:= k+1 to N do
		begin
			m.pon_e(k, j, gkgj(k,j));
			m.pon_e(j, k, m.e(k,j));
		end;
	end;
	res:= M.Escaler( a );
	if res .....
	for j:= 1 to N do
		ponMultiplicador(j, a.e(j));
	M.done;
	a.done;
end;


{ C lculo de la funci¢n del problema relajado }
function frel(kf: integer; var X: TVectR ): NReal; far;
var
	k: integer;
	m, w: NReal;
	pr: PFichaFunc;

begin
	m:=0;

{ Primero recorremos las restricciones }
	for k:= 1 to NRs do
	begin
		pr:= @ProblemaPtr^[k];
		if pr^.Activa then

		w:=pr^.f(k, x);
		pr^.uv:= w;
		if w > 0 then
		begin
			pr^.Activa:= true;
			m:= m+pr^.rel*w;
		end
		else
			pr^.Activa:= false;
	end;

{ Llamamos la funci¢n de costos del problema no relajado }
	m:=	ProblemaPtr^[0].f(0, x)+ m;
	frel:= m;
end;


{ computa el gradiente de la funci¢n relajada }
procedure grel(kf: integer; var grad, X: TVectR ); far;
var
	k: integer;
	mg: TVectR;
	TempRel: NReal;
	norm_r, norm_f: NReal;
	pr: PFichaFunc;

begin
	mg.init( X.N );
	mg.Ceros;
	for k:= 1 to NRs do
	begin
		pr:= ProblemaPtr^[k];
		if pr^.activa> 0 then
		begin
			mg.sumRPV( pr^.rel,  ug );
		end;
	end;
	pr:= ProblemaPtr^[0];
	grad.igual(pr^.ug);
	grad.Sum( 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.

