/////////////////////////////////////////////////
//     LSMA 651, 2010-2011
//        Mr Dumas
//
//gradient or Newton method with linesearch
// by backtracking and Armijo condition
//////////////////////////////////////////////////
function y=J(x)
  y=100*(x(2)-x(1)^2)^2+(x(1)-1)^2; //  fonction 1 : banana shape
 // y=0.1*(x(1))^2+(x(2))^2;          // fonction 2: quadratique
endfunction  
//
function y=nablaJ(x)   // gradient
  y=[-400*x(1)*(x(2)-x(1)^2)+2*(x(1)-1);200*(x(2)-x(1)^2)];    //  fonction 1 
 // y=[2*0.1*(x(1));2*(x(2))];                                 // fonction 2
endfunction
//
function H=hessianJ(x)  // hessian
  H=[-400*(x(2)-x(1)^2)+800*x(1)^2+2,-400*x(1);-400*x(1),200]; //  fonction 1 
 // H=[0.2,0;0,2];                                             // fonction 2
endfunction
//
function z=my_surface(x,y)
  z=100*(y-x^2)^2+(x-1)^2;      // fonction 1
 // z=0.1*(x)^2+(y)^2;          // fonction 2
endfunction
////////////////////////////////////////////////////////////////////
function alpha=backtrackA(x,J,Jx,p,gx,betaA,tau,alphainit) 
alpha=alphainit;
count=1;
while(J(x+alpha*(p))>Jx+betaA*alpha*(p'*gx))&(count<20)
  count=count+1;
  alpha=tau*alpha;
end  
endfunction  
// 
function p=newton(B,gx)
p=linsolve(B,gx); // solve Bp=-gx
endfunction
//////////////////////////
betaA=0.0001;
alphainit=1;
tau=0.7;
epsilon=1E-8;
countmax=2000;
///////////////////////////
//
n=2;
//
x=ones(n,1)-2*rand(n,1); // starting point
disp('initial value of x');disp(x);
gx=nablaJ(x);
Hx=hessianJ(x);
Jx=J(x);
count=1;
Xtot=[x];
B=eye(n,n);
//
while (norm(gx)>epsilon)&(count<countmax)
  p=newton(Hx,gx);                                  //newton
//  p=-gx;                                            // gradient
 alpha=backtrackA(x,J,Jx,p,gx,betaA,tau,alphainit); // linesearch 
 //alpha=alphainit;                                 // no line search
 xold=x;//
 gxold=gx;
 x=x+alpha*(p);  
 Jx=J(x);
 gx=nablaJ(x);
 Hx=hessianJ(x);
 count=count+1;
 Xtot=[Xtot,x];
end
disp('final value of x:');disp(x)
disp('iteration number');disp(count)
//
t=linspace(-2,2,30);
clf()
contour(t,t,my_surface,0:2:20)
plot2d(Xtot(1,:),Xtot(2,:))

