//
// Rozenbrock facile
//
function y = rozen1(x)
   y = (x(2)-x(1)^2)^2 + (1-x(1))^2
endfunction

function df = grd_rozen1(x)
   cte = 2*(x(2)-x(1)^2)
   df = [ -2*x(1)*cte-2*(1-x(1)) ;  cte ];
endfunction
 
//
// Rozenbrock plus difficile
//
function y = rozen2(x)
   y = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2
endfunction

function df = grd_rozen2(x)
   cte = 200*(x(2)-x(1)^2)
   df = [ -2*x(1)*cte-2*(1-x(1)) ;  cte ];
endfunction

// fonction helical
function y = helical(x)
   // xopt = [1;0;0]; x0 = [-1;0;0]
   if x(1) > 0 then, theta = 0, else, theta = 5, end
   v = [10*(x(3) - (5/%pi)*atan(x(2)/x(1)) - theta);...
	10*(sqrt(x(1)^2 + x(2)^2) - 1);...
	x(3)];
   y = v'*v
endfunction

function df = grd_helical(x)
   // 
   if x(1) > 0 then, theta = 0, else, theta = 5, end
   r = x(2)/x(1)
   d = sqrt(x(1)^2 + x(2)^2)
   v = [10*(x(3) - (5/%pi)*atan(r) - theta);...
	10*(d - 1);...
	x(3)];

   coef1 = (50/%pi)/((1+r^2)*x(1));
   coef2 = 10/d;
   Jac = [  coef1*r   ,-coef1      , 10;...
	    coef2*x(1), coef2*x(2) , 0 ;...
	    0         , 0          , 1];
   df = 2*Jac'*v
endfunction

// fonction expsquares (dimension variable)
function f = expsquares(x)
   n = length(x);
   f = exp(-sum(x)) + 0.5*sum(((1:n)'.*x).^2);
endfunction

function g = grd_expsquares(x)
   n = length(x);
   vect = (1:n)';
   c1 = exp(-sum(x));
   g = -c1 + vect.^2 .*x;
endfunction

function [fopt,xopt] = sol_expsquares(n)
// xopt(j) = exp(-y))/j^2
// y sol de (1 + ... + 1/n^2)exp(-y) = y
   C = sum(1./(1:n).^2);
   y = 1 + 0.5*(log(C)-1);
   while %t
      temp = exp(-y)
      dy = (C*temp-y)/(C*temp+1)
      y = y + dy
      if ( abs(dy/y) <= 2*%eps ) then, break, end
   end
   xopt = exp(-y)./((1:n)').^2;
   fopt = expsquares(xopt);
endfunction