% try several minimization methods

function minimize
    global outfile
    format compact
%    outfile = fopen('min.dat','w');
%	fprintf(outfile,'Hello\n');

    xgrid = -3:0.1:2
    ygrid = -2:0.1:3
    Z = gridf(xgrid,ygrid);
	meshc(xgrid,ygrid,Z);
    stop
	X0 = [1;4/5];
%	r  = newton(@myfuncGrad,@myfuncH,X0,1.e-8,7); % white line

	r = steepest_descent(@myfunc,@myfuncGrad,X0,2,1.e-8,30); % blue line

end

function Z=gridf(X,Y)
	fprintf('Starting gridf\n');
	nx=numel(X)
	ny=numel(Y)
	Z = zeros(ny,nx);
    k=0;
    for i=1:nx
       for j=1:ny
          k=k+1;
          Z(j,i) = myfunc([X(i),Y(j)]);
	   end
	end
	fprintf('Done with gridf\n');
end

function z = myfunc(X)
	x = X(1);
    y = X(2);
    z = x^3 + y^3 + 3*x^2 - 3*y^2 - 8;
end

function Df = myfuncGrad(X)
	x = X(1);
    y = X(2);
    Dfx = 3*x^2 + 6*x;
	Dfy = 3*y^2 - 6*y;
    Df = [Dfx,Dfy]';
end

function D2f = myfuncH(X)
	x = X(1);
    y = X(2);
    Dfxx = 6*x + 6;
	Dfxy = 0;
	Dfyy = 6*y - 6;
    D2f = [Dfxx,Dfxy; Dfxy,Dfyy];
end


function r = newton(f,df,x0,abstol,maxits)
	r = x0
	for i = 1:maxits
		dx = -df(r)\f(r);
		r = r + dx
		if norm(dx) < abstol
			fprintf('Converged to desired accuracy in %i steps\n',i)
			return;
		end		
	end
	error('Failed to converge in %i steps',maxits)
end

function r = steepest_descent(f,df,x0,dist,abstol,maxits)
	global outfile
    global theX0 theU theObjective
	theObjective = f;
	r = x0
	range = dist;
	for i = 1:maxits
		% do a line minimization along the gradient of f at r0
		theX0 = r;
		gradf = df(r);
		theU = -gradf/norm(gradf);  % unit downhill vector
%    	myfplot(@func_on_line,[0,range]);
%       hold on;
        t = gss(@func_on_line,0,range,abstol/10,20);
		r = r + t*theU
		fprintf(outfile,'%f24.15 %f24.15 \n',r(1),r(2));
		if abs(t) < abstol
			fprintf('Converged to desired accuracy in %i steps\n',i)
			return;
		end		
	end
%   hold off;
end

function z = func_on_line(t)
    global theX0 theU theObjective
	z = theObjective(theX0 + t*theU);
%    xandf = [x,z]
end


function myfplot(f,interval)
	n = 200;
	x = linspace(interval(1),interval(2),n);
	y = zeros(1,n);
    for i=1:n
	   y(i) = f(x(i));
    end
	plot(x,y,'b');
end