// initial signals
for i=1:20 f(i)=20*sin(i-1)+20; end

for i=1:20 p(i)=0.5+3*i; end

// Helper function for the energy function of exercise 1
function g = ex1function(u,f,h,alpha)
  n=length(f)
  g = 0
  for i=1:n
    g=g+0.5*(u(i)-f(i))^2
  end
  for i=1:n-1
    g=g+alpha*0.5*((u(i+1)-u(i))/h)^2
  end
  
endfunction

// Helper function for the gradient of exercise 1
function g=ex1Gradient(u,f,h,alpha)
  N=length(u)
  g(1)=(u(1)-f(1))-alpha/h^2*(u(2)-u(1))
  for i=2:N-1
    g(i)=(u(i)-f(i))-alpha/h^2*(u(i-1)-2*u(i)+u(i+1))
  end
  g(N)=(u(N)-f(N))+alpha/h^2*(u(N)-u(N-1))
endfunction

// Helper function for the Hessian of exercise 1
function H=ex1Hessian(u,f,h,alpha)
  N=length(u)
  H(N,N)=0
  H(1,1)=1+alpha/h^2
  for i=2:N-1
    H(i,i)=1+2*alpha/h^2
    H(i-1,i)=-alpha/h^2
    H(i,i-1)=-alpha/h^2  
  end
  H(N-1,N)=-alpha/h^2
  H(N,N-1)=-alpha/h^2
  H(N,N)=1+alpha/h^2
endfunction

// Gradient Descent Line Search Algorithm
function u=ex1GDLS(f,h,alpha,alphaGDLS,betaGDLS)
  N=length(f)
  eps = 0.0001
  err = 1000
  u=f
  iter = 0
  // Stopping criterion, if |grad(x^k+1)-grad(x^k)| < eps
  while (err > eps)
    d=-ex1Gradient(u,f,h,alpha)
    // Line Search
    gammaGDLS = 1
    while (ex1function(u+gammaGDLS*d,f,h,alpha)> ex1function(u,f,h,alpha) + alphaGDLS*gammaGDLS*d'*d)
      gammaGDLS= betaGDLS*gammaGDLS
     // printf("%f\n",gammaGDLS)
    end
    uold=u
    u=u+gammaGDLS*d
    iter=iter+1
    err = norm(ex1Gradient(u,f,h,alpha)-ex1Gradient(uold,f,h,alpha))
    //printf("err %f\n",err)
  end
  printf("Number of iterations: %i\n",iter)
  endfunction


// Newton Algorithm without stepsize control

function u=ex1Newton(f,h,alpha)
  N=length(f)
  eps = 0.0001
  err = 1000
  u=f
  // Stopping criterion, if |grad(x^k+1)-grad(x^k)| < eps
  while (err > eps)
    d=-inv(ex1Hessian(u,f,h,alpha))*ex1Gradient(u,f,h,alpha)
    uold=u
    u=u+d
    err = norm(ex1Gradient(u,f,h,alpha)-ex1Gradient(uold,f,h,alpha))
    printf("err %f\n",err)
  end
endfunction

//**************************************************//

//
function f = ex2Banana(v)
  x=v(1)
  y=v(2)
  f=(1-x)^2+100*(y-x^2)^2
endfunction

function f = ex2BaSh(v)
  x=v(1)
  y=v(2)
  f=(x-2)^4 + (x-2*y)^2
endfunction


function g = ex2GradientBanana(v)
  x=v(1)
  y=v(2)
  g(1)=-2*(1-x)-400*x*(y-x^2)
  g(2)=200*y-200*x^2
endfunction

function g = ex2GradientBaSh(v)
  x=v(1)
  y=v(2)
  g(1)=4*(x-2)^3+2*x-4*y
  g(2)=-4*x+8*y
endfunction

function H= ex2HessianBanana(v)
  x=v(1)
  y=v(2)
  H(1,1)=2+1200*x^2-400*y
  H(2,1)=-400*x
  H(1,2)=-400*x
  H(2,2)=200
endfunction

function H= ex2HessianBaSh(v)
  x=v(1)
  y=v(2)
  H(1,1)=12*(x-2)^2+2
  H(2,1)=-4
  H(1,2)=-4
  H(2,2)=8
endfunction

function u=ex2BananaGDLS(f,alphaGDLS,betaGDLS)
  N=length(f)
  eps = 0.0001
  err = 100
  u=f
  numb=0
  // Stopping criterion, if |grad(x^k+1)-grad(x^k)| < eps
  while (err > eps)
    d=-ex2GradientBanana(u)
    // Line Search
    gammaGDLS = 1
    while (ex2Banana(u+gammaGDLS*d)> ex2Banana(u) + alphaGDLS*gammaGDLS*d'*d)
      gammaGDLS= betaGDLS*gammaGDLS
     // printf("%f\n",gammaGDLS)
    end
    uold=u
    u=u+gammaGDLS*d
    err = norm(ex2GradientBanana(u)-ex2GradientBanana(uold))
    numb=numb+1
    //    printf("err %f\n",err)
  end
  printf("Number of iterations: %i\n",numb)
endfunction

function u=ex2BaShGDLS(f,alphaGDLS,betaGDLS)
  N=length(f)
  eps = 0.0001
  err = 100
  u=f
  numb=0
  // Stopping criterion, if |grad(x^k+1)-grad(x^k)| < eps
  while (err > eps)
    d=-ex2GradientBaSh(u)
    // Line Search
    gammaGDLS = 1
    while (ex2BaSh(u+gammaGDLS*d)> ex2BaSh(u) + alphaGDLS*gammaGDLS*d'*d)
      gammaGDLS= betaGDLS*gammaGDLS
     // printf("%f\n",gammaGDLS)
    end
    uold=u
    u=u+gammaGDLS*d
    err = norm(ex2GradientBaSh(u)-ex2GradientBaSh(uold))
    numb=numb+1
    //    printf("err %f\n",err)
  end
  printf("Number of iterations: %i\n",numb)
endfunction

function u=ex2BananaNewton(f)
  N=length(f)
  eps = 0.0001
  err = 100
  u=f
  numb=0
  // Stopping criterion, if |grad(x^k+1)-grad(x^k)| < eps
  while (err > eps)
    d=-inv(ex2HessianBanana(u))*ex2GradientBanana(u)
    uold=u
    u=u+d
    err = norm(ex2GradientBanana(u)-ex2GradientBanana(uold))
    numb=numb+1
    //    printf("err %f\n",err)
  end
  printf("Number of iterations: %i\n",numb)
endfunction

function u=ex2BaShNewton(f)
  N=length(f)
  eps = 0.0001
  err = 100
  u=f
  numb=0
  // Stopping criterion, if |grad(x^k+1)-grad(x^k)| < eps
  while (err > eps)
    d=-inv(ex2HessianBaSh(u))*ex2GradientBaSh(u)
    uold=u
    u=u+d
    err = norm(ex2GradientBaSh(u)-ex2GradientBaSh(uold))
    numb=numb+1
    //    printf("err %f\n",err)
  end
  printf("Number of iterations: %i\n",numb)
endfunction

function u=ex2BananaTrust(f)
  N=length(f)
  eps = 0.0001
  err = 100
  u=f
  numb=0
  // Stopping criterion, if |grad(x^k+1)-grad(x^k)| < eps
  while (err > eps)
    H=ex2HessianBanana(u)
    minEV = min(real(spec(H)))
    if (minEV<=0) 
      alpha = abs(minEV)+1
    else 
      alpha = 0
    end
    printf("Alpha: %f\n",alpha)
    l(1:N)=1
    I=diag(l)
    d=-inv(H+alpha*I)*ex2GradientBanana(u)
    uold=u
    u=u+d
    err = norm(ex2GradientBanana(u)-ex2GradientBanana(uold))
    numb=numb+1
    //    printf("err %f\n",err)
  end
  printf("Number of iterations: %i\n",numb)
endfunction

function u=ex2BaShTrust(f)
  N=length(f)
  eps = 0.0001
  err = 100
  u=f
  numb=0
  // Stopping criterion, if |grad(x^k+1)-grad(x^k)| < eps
  while (err > eps)
    H=ex2HessianBaSh(u)
    minEV = min(real(spec(H)))
    if (minEV<=0) 
      alpha = abs(minEV)+1
    else 
      alpha = 0
    end
    printf("Alpha: %f\n",alpha)
    l(1:N)=1
    I=diag(l)
    d=-inv(H+alpha*I)*ex2GradientBaSh(u)
    uold=u
    u=u+d
    err = norm(ex2GradientBaSh(u)-ex2GradientBaSh(uold))
    numb=numb+1
    //    printf("err %f\n",err)
  end
  printf("Number of iterations: %i\n",numb)
endfunction

function u=ex3b(u)
  err=1000
  iter=0
    while (err > 0.00001)
    fp =4*u^3+9*u^2-6*u-7
    fpp=12*u^2+18*u-6
    uold=u
    u=u-fp/fpp
    iter =iter+1
    err=abs(u-uold)
  end
  printf("Number of iterations: %i\n",iter)
endfunction


function u=ex3c(u,a)
  err=1000
  iter=0
    while (err > 0.000001)
    fp =1/u-a
    fpp=-1/u^2
    uold=u
    u=u-fp/fpp
    iter =iter+1
    err=abs(u-uold)
  end
  printf("Number of iterations: %i\n",iter)
endfunction

function u=ex3d(u)
  err=1000
  iter=0
    while (err > 0.000001)
    fp =atan(u)
    fpp=1/(u^2+1)
    uold=u
    u=u-fp/fpp
    iter =iter+1
    err=abs(u-uold)
    printf("Number of iterations: %i\n",iter)
  end
endfunction


function u=ex4Armijo(f,alphaGDLS,beta1,beta2)
  N=length(f)
  eps = 0.001
  err = 1000
  u=f
  lambda =0.5
  iter = 0
  //gammaGDLS=1
  //sigma=-gammaGDLS*ex1Gradient(u,f,h,alpha)
  
  // Stopping criterion, if |grad(x^k+1)-grad(x^k)| < eps
  while (err > eps)
    gammaGDLS=1
    d=ex2GradientBanana(u)
    sigma=-gammaGDLS*ex2GradientBanana(u)'*d/norm(d)^2
     while (ex2Banana(u+sigma*d)> ex2Banana(u) + alphaGDLS*sigma*d'*d)
    sigma=lambda*beta1*sigma+(1-lambda)*beta2*sigma
    end
    uold=u
    u=u+sigma*d
    iter=iter+1
    err = norm(ex2GradientBanana(u))
    printf("err %f\n",err)
  end
  printf("Number of iterations: %i\n",iter)
endfunction
