function [xmin,fmin]=Newton(dim,eps)

fun=3;             % selection the function we are going to optimize fun=3 is the one used in the exercise
FLAG=1; 
xcur=randn(dim,1); % initial point - drawn from Gaussian distribution

counter=0;
display(['Iteration: ',num2str(counter),' - current value: ',num2str(f(xcur,fun))]);
fiter =[f(xcur,fun)];
xiter =[xcur(1:2)];

while(FLAG)
 % get gradient and Hessian at current point
 gradfval = gradf(xcur,fun);
 Hessfval = Hessf(xcur,fun);
 
 % compute the Newton descent direction
 descent = -Hessfval\gradfval;

 % get the stepsize with the Armijo rule
 alpha = getStepSize(xcur,gradfval,descent,fun);
   
 % update the current point
 xcur = xcur + alpha*descent;
 counter=counter+1;
 
 % display function value during iteration
 if(counter<100 | rem(counter,100)==0)
  display(['Iteration: ',num2str(counter),' - current value: ',num2str(f(xcur,fun)),'; current stepsize: ',num2str(alpha) ]);
 end
 
 % saves function values (note that this way of growing an array is quite
 % inefficient in Matlab since it COPIES the whole array -> super slow if
 % the array is quite large)
 fiter=[fiter f(xcur,fun)];
 xiter=[xiter xcur(1:2)];
 
 % stopping condition with the Newton decrement 
 %lambda=gradfval * Hessfval^{-1} *gradfval 
 %      =-gradfval*desent
 if(abs(gradfval'*descent)<eps) FLAG=0; end
end

% first we plot the difference to the optimum (exact) during the iteration
figure 
hold on
switch(fun)
    case 1, plot(log10(fiter-f([0 0]',fun)),'k.','MarkerSize',10);
            plot(log10(fiter-f([0 0]',fun)),'k-');
    case 2, plot(log10(fiter-f([0 0]',fun)),'k.','MarkerSize',10);
            plot(log10(fiter-f([0 0]',fun)),'k-');
    case 3, plot(log10(fiter-f([-log(2)/2 0]',fun)),'k.','MarkerSize',10);
            plot(log10(fiter-f([-log(2)/2 0]',fun)),'k-');
end
ylabel('$log(f(x^k)-p^*)$','FontSize',10,'Interpreter','latex')
xlabel('k','FontSize',10)
hold off

% plots the path of the sequence
figure, 
hold on
[X,Y]=meshgrid(-2:0.05:2,-1:0.05:1);

switch(fun)
    case 1, Z = exp(1/dim*(X.^2 +Y.^2));
    case 2, u1=[1 -0.2]'; u1=u1/norm(u1); u2=[0.2 1]'; u2=u2/norm(u2); r1=0.1; r2=20;
            Z= r1 * (u1(1)*X+u1(2)*Y).^2 + r2*(u2(1)*X + u2(2)*Y).^2;
    case 3, Z=exp(X+3*Y-0.1) + exp(X-3*Y-0.1)+exp(-X-0.1);
end
colormap summer
%surf(X,Y,Z); 
contour(X,Y,Z,100);
samples = [xiter', fiter'];
%scatter3(samples(:,1),samples(:,2),samples(:,3)+1E-1,10,'k','filled');
plot(samples(:,1),samples(:,2),'-k');
plot(samples(:,1),samples(:,2),'.k','MarkerSize',10);
title('Path of the sequence','F);
axis equal
hold off







function fval = f(point,fun)
  switch(fun)
    case 1, dim=length(point); fval = exp(norm(point).^2/dim);
    case 2, u1=[1 -0.2]'; u1=u1/norm(u1);
            u2=[0.2 1]'; u2=u2/norm(u2);
            r1=0.1; r2=20;
            fval = r1*sum(u1.*point)^2 + r2*sum(u2.*point)^2;
    case 3, fval = exp(point(1)+3*point(2)-0.1) + exp(point(1)-3*point(2)-0.1)+exp(-point(1)-0.1);
  end

function gradfval = gradf(point,fun)
  switch(fun)
    case 1, dim=length(point); gradfval = 2*point*1/dim*exp(1/dim*norm(point).^2);
    case 2, u1=[1 -0.2]'; u1=u1/norm(u1);
            u2=[0.2 1]'; u2=u2/norm(u2);
            r1=0.1; r2=20;
            gradfval = 2*r1*sum(u1.*point)*u1 + 2*r2*sum(u2.*point)*u2;
    case 3, gradfval(1,1)=exp(point(1)+3*point(2)-0.1) + exp(point(1)-3*point(2)-0.1)-exp(-point(1)-0.1); 
            gradfval(2,1)=3*exp(point(1)+3*point(2)-0.1) -3*exp(point(1)-3*point(2)-0.1);
  end
 
function Hessfval = Hessf(point,fun)
  switch(fun)
      case 1, dim=length(point); Hessfval =(4*point*point'/dim^2 + 2*eye(dim,dim)/dim)*exp(norm(point).^2/dim);
      case 2, u1=[1 -0.2]'; u1=u1/norm(u1);
              u2=[0.2 1]'; u2=u2/norm(u2);
              r1=0.1; r2=20;
              Hessfval = r1*u1*u1' + r2*u2*u2';
      case 3, Hessfval(1,1)=exp(point(1)+3*point(2)-0.1) + exp(point(1)-3*point(2)-0.1)+exp(-point(1)-0.1);
              Hessfval(1,2)=3*exp(point(1)+3*point(2)-0.1) -3*exp(point(1)-3*point(2)-0.1);
              Hessfval(2,1)=Hessfval(1,2);
              Hessfval(2,2)=9*exp(point(1)+3*point(2)-0.1) +9*exp(point(1)-3*point(2)-0.1);
  end
  

function stepsize=getStepSize(xcur,gradfval,descent,fun)
  alpha = 0.1; beta=0.5;
  t=1; FLAG=1;

  while(FLAG)
    xnew = xcur + t*descent;
    if( f(xnew,fun) > f(xcur,fun) + alpha*t*gradfval'*descent)
     t=beta*t;
    else
     FLAG=0;
    end
  end
  stepsize=t;