unit maximize;

(*
   PROJECT          : StatPascal library
   FILE-SPEC.       : maximize.sp
   AUTHOR           : R.-D. Reiss, M. ThomasM
   CREATION DATE    : Apr-5-2000
   LAST UPDATE      : Dec-6-2000
   VERSION          : 1.0
   RELATED FILES    : See tmle.sp for usage of routine.
   FUNCTIONAL DESCR.: Unit implementing maximization using the
                      Newton-Raphson iteration. Derivates are
                      calculated numerically.
   SPECIAL FEATURES : requires Xtremes Version 3.0
*)

interface

type realvecfunc = function (realvector): real;
type boolvecfunc = function (realvector): boolean;

function searchmax (f: realvecfunc; start: realvector; 
                    constraints: boolvecfunc; dimension: integer): realvector;

implementation

function searchmax;
    const eps = 1e-3;
    var xnew, xold, hv, fs, d: realvector;
        fnew, fold, s : real;
	steps, i, j: integer;
	H: realmatrix;
	valid: boolean;
    begin
	xnew := start;
	fnew := f (xnew);
	steps := 0;	
	repeat
	    fold := fnew;
	    xold := xnew;
	    if size (hv) <> 0 then hv := hv [combine (false)];	% delete vectors hv and fs 
	    if size (fs) <> 0 then fs := fs [combine (false)];
	    for i := 1 to dimension do begin
	 	fs := combine (fs, (f (xold + eps * unitvector (i, dimension)) - fold) / eps);
		for j := 1 to dimension do 
		    hv := combine (hv, (f (xold + eps * (unitvector (i, dimension) + unitvector (j, dimension)))
                                       -f (xold + eps * unitvector (i, dimension))
                 		       -f (xold + eps * unitvector (j, dimension)) + fold) / sqr (eps))
	    end;
	    H := makematrix (hv, dimension, dimension);
	    d := invert (H) * fs;
	    s := 1.0;
	    repeat
		xnew := xold - s * d;
		s := s / 2.0;
		valid := constraints (xnew);
		if valid then fnew := f (xnew)
	    until (valid and (fnew > fold)) or (s < 1e-10);
	    steps := steps + 1
	until (steps = 200) or (abs (fnew - fold) < 1e-10);
	return xnew
    end;

end.