//
// methode de gradient ou de Newton
// our minimiser la fonction de Rosenbrock (ou banane)
//
function y=J(x)
  y=100*(x(2)-x(1)^2)^2+(x(1)-1)^2;
endfunction  
//
function y=nablaJ(x)
  y=[-400*x(1)*(x(2)-x(1)^2)+2*(x(1)-1);200*(x(2)-x(1)^2)];
endfunction
//
function y=hessienJ(x)
  y=[-400*x(2)+1200*x(1)^2+2,-400*x(1);-400*x(1),200];
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  
//
betaA=0.0001;
alphainit=1;
tau=0.7;
epsilon=1E-3;
//
n=2;
//
xini=ones(n,1)-2*rand(n,1); // starting point
x=xini;
disp('initial value of x');disp(x);
gx=nablaJ(x);
Hx=hessienJ(x);
Jx=J(x);
count=1;
Xtot=[x];
B=eye(n,n);
//
while (norm(gx)>epsilon)&(count<2000)
  ///
  /// A CHOISIR PAR L'UTILISATEUR
  //
  p=-gx;    // methode de gradient
//    p=linsolve(Hx,gx) // methode de Newton
  //
 alpha=backtrackA(x,J,Jx,p,gx,betaA,tau,alphainit); // linesearch 
 xold=x;
 gxold=gx;
 x=x+alpha*(p);  
 Jx=J(x);
 gx=nablaJ(x);
 Hx=hessienJ(x);
 count=count+1;
 Xtot=[Xtot,x];
end
disp('final value of x:');disp(x)
disp('iteration number');disp(count)
//
function f = rosenbrockC ( x1 , x2 )
  x = [x1 x2];
  f = 100.0 *(x(2)-x(1)^2)^2 + (1-x(1))^2;
endfunction
clf()
xdata = linspace(-2,2,100);
ydata = linspace(-2,2,100);
contour ( xdata , ydata , rosenbrockC , [1 10 100 500 1000])
plot2d(Xtot(1,:),Xtot(2,:))
//
[ fopt , xopt ] = optim ( rosenbrock , xini)
disp('solution Scilab (optim):');disp(xopt);




