function [value,px,py]=project2(n)
% first plot out functions f, g, and h.
npts=500;
x=linspace(-0.8,0.5,npts);
y=linspace(-0.8,0.5,npts);

f=zeros(npts,npts);
for i=1:npts
    for j=1:npts
        f(i,j)=fun(x(i),y(j));
    end
end

% this is the minimal value among 250,000 functions values which should
% give a crude estimate about the global minimum of the function
min(min(f))

% function f is highly oscillatory and its graph gives us very 
% little information. 
% mesh(x,y,f)
% % press the enter key to continue
% pause

g=fun1(x);
h=fun2(y);

% By zooming in, we see that min g(x) <-0.5 and min h(y)<-1.6, and hence
% min f(x,y) < -0.5-1.6+0.5 = -1.6

% For |x|, |y|>2, we have f(x,y)>=e^(-1)-1-1-1-0.5+1/4(4+4)>-1.13>-1.6>min
% f(x,y). Thus we may restrict our search to the rectangular domain
% [-2,2]x[-2,2].

subplot(1,2,1);
plot(x,g,'r')
subplot(1,2,2);
plot(y,h,'b')

% find all local minima of g and h on the interval [-2,2] with g<=1 and
% h<=-0.4.
rootx=findx0(-2,2,1,@fun1,@fun1x,@fun1xx);
rooty=findx0(-2,2,-0.4,@fun2,@fun2x,@fun2xx);

nx=length(rootx)
ny=length(rooty)

value=100;

% Use the positions of local minima of g and h as starting points together
% with Newton's method to find all local minima of the original function f.
for i=1:nx,
    for j=1:ny,
        p=zeros(2,1);
        p(1)=rootx(i);
        p(2)=rooty(j);
        p1=newton(p,1e-14,10,@gfun,@jfun);
        x=p1(1);
        y=p1(2);
        v=fun(x,y);
        
        % compare and keep the smaller value
        if v < value
            value = v;
            px=x;
            py=y;
        end
    end
end

function root=findx0(a,b,fbnd,f,fx,fxx)

n=2000;
h=(b-a)/n;
x=a:h:b;

% divide the whole interval [a,b] into very small intervals, use 
% intermediate value theorem to identify the intervals on which the zeros 
% of fx are located. use the midpoints of these small intervals as
% the starting points of Newton's method. 
y=feval(fx,x);
y1=y(1:end-1).*y(2:end);
ind=find(y1<0);

% x0 is a vector containing starting points for local minima of g and h.
x0=(x(ind)+x(ind+1))/2;

j=1;
for i=1:length(ind),
    r=newton(x0(i),1e-14,10,fx,fxx); 
    f1=feval(f,r);
    f2=feval(fxx,r);
    if (f1<fbnd && f2 > 0) 
        root(j)=r;
        j=j+1;
    end
end


function root = newton(x0,error_bd,max_iterate,f,dfdx)
%
% function newton(x0,error_bd,max_iterate,f,dfdx)
%
% This is Newton's method for solving an equation f(x) = 0.
%
% It works for both 1-dimensional and multidimensional problems.
%
% The functions f(x) and dfdx(x) are given as input parameters.
% The parameter error_bd is used in the error test for the 
% accuracy of each iterate.  The parameter max_iterate
% is an upper limit on the number of iterates to be 
% computed.  An initial guess x0 must also be given.
%

error = ones(size(x0));
it_count = 0;
while norm(error) > error_bd && it_count <= max_iterate
    fx = feval(f,x0);
    dfx = feval(dfdx,x0);
    if norm(dfx) == 0
        disp('The Jacobian is zero.  Stop')
        return
    end
    % Newton's method. For multidimensional case, "\" requires solving a
    % linear system.
    x1 = x0 - dfx\fx;
    error = x1 - x0;

    x0 = x1;
    it_count = it_count + 1;
end

root=x1;

% functions g(x), g'(x), g''(x)
function y=fun1(x)
y = exp(sin(40*x))+sin(60*sin(x))+1/4*(x.^2);

function y=fun1x(x)
y = 40*cos(40*x).*exp(sin(40*x))+60*cos(60*sin(x)).*cos(x)+1/2*x;

function y=fun1xx(x)
y = -1600*sin(40*x).*exp(sin(40*x))+1600*cos(40*x).^2.*exp(sin(40*x))-3600*sin(60*sin(x)).*cos(x).^2-60*cos(60*sin(x)).*sin(x)+1/2;

% functions h(y),h'(y), h''(y)
function y=fun2(x)
y=sin(50*exp(x))+sin(sin(70*x))+1/4*(x.^2);

function y=fun2x(x)
y=50*cos(50*exp(x)).*exp(x)+70*cos(sin(70*x)).*cos(70*x)+1/2*x;

function y=fun2xx(x)
y=-2500*sin(50*exp(x)).*exp(x).^2+50*cos(50*exp(x)).*exp(x)-4900*sin(sin(70*x)).*cos(70*x).^2-4900*cos(sin(70*x)).*sin(70*x)+1/2;

% original function f(x,y)
function z=fun(x,y)
z=exp(sin(40*x))+sin(50*exp(y))+sin(60*sin(x))+sin(sin(70*y))-1/2*sin(10*(x+y))+1/4*(x^2+y^2);

% first derivatives
function z=funx(x,y)
z=40*cos(40*x)*exp(sin(40*x))+60*cos(60*sin(x))*cos(x)-5*cos(10*x+10*y)+1/2*x;

function z=funy(x,y)
z=50*cos(50*exp(y))*exp(y)+70*cos(sin(70*y))*cos(70*y)-5*cos(10*x+10*y)+1/2*y;

% gradient f
function dfdx=gfun(p)
x=p(1);
y=p(2);
dfdx=zeros(2,1);
dfdx(1)=funx(x,y);
dfdx(2)=funy(x,y);

% second derivatives 
function zxx=funxx(x,y)
zxx = -1600*sin(40*x)*exp(sin(40*x))+1600*cos(40*x)^2*exp(sin(40*x))-3600*sin(60*sin(x))*cos(x)^2-60*cos(60*sin(x))*sin(x)+50*sin(10*x+10*y)+1/2;

function zyy=funyy(x,y)
zyy=-2500*sin(50*exp(y))*exp(y)^2+50*cos(50*exp(y))*exp(y)-4900*sin(sin(70*y))*cos(70*y)^2-4900*cos(sin(70*y))*sin(70*y)+50*sin(10*x+10*y)+1/2;

function zxy=funxy(x,y)
zxy=50*sin(10*x+10*y);

% Jacobian of the gradient of f
function J=jfun(p)
x=p(1);
y=p(2);

J=zeros(2);
J(1,1)=funxx(x,y);
J(1,2)=funxy(x,y);
J(2,1)=J(1,2);
J(2,2)=funyy(x,y);