function curvefittingSD

% This code exemplifies the steepest descent applied to 
% a curve fitting problem

close all, clf
global xdata ydata

xdata=-5:5;
ydata=3*exp(-(xdata-1).^2/4)-0.1*rand(1,length(xdata));

tol=10^(-2);

rho =1/4;
c=10^(-4);

x=[1;0;0]; % initial guess

[y,gy] =myfunc(x);
p=-gy;
counter=0;

while norm(gy)>tol
    
    alphak=linesearch(@myfunc,p,x,rho,c);
    x=x+alphak*p;
    p=-gy;
    [y,gy]=myfunc(x);
    counter=counter+1

    clf, plot(xdata,ydata,'o')
    axis([-5 5 0 3.6])
    hold on
    xplot=-5:0.01:5;
    yplot= x(1)*exp(-x(2)*(xplot-x(3)).^2);
    plot(xplot,yplot);
    legend('data set', 'approximation')
    text(-4.5,3.4, strcat('iter=',num2str(counter)))
    text(-4,3.2,strcat('a=',num2str(x(1))))
    text(-4,3,strcat('b=',num2str(x(2))))
    text(-4,2.8,strcat('c=',num2str(x(3))))
    
    text(-4,2.4,strcat('Residual=',num2str(y)))
    text(-4,2.2,strcat('|Gradient|=',num2str(norm(gy))))
    title('Curve fitting via steepest descent, with phi(x)=a exp(-b(x-c)^2)')
    pause(0.01)

end

text(-3,1,strcat('Optimization terminated in ',' ',num2str(counter),' ',' iterations, gradient less than ',' ',num2str(tol)))
end

%%%%%%%%%%%%%%%%%%%%%%%%

function [fx,gx] = myfunc(x)

global xdata ydata

a=x(1); 
b=x(2);
c=x(3);

phi = @(y) a*exp(-b*(y-c).^2);
phia= @(y) exp(-b*(y-c).^2);
phib = @(y) a*exp(-b*(y-c).^2).*(-(y-c).^2);
phic = @(y) a*exp(-b*(y-c).^2).*(2*b*(y-c));


fx =  sum((phi(xdata)-ydata).^2);
gx = [sum(2*(phi(xdata)-ydata).* phia(xdata));
      sum(2*(phi(xdata)-ydata).* phib(xdata));
      sum(2*(phi(xdata)-ydata).* phic(xdata))];
    
end

%%%%%%%%%%%%%%%%%%%%%%%

function alphak = linesearch(f,d,x,rho,c)

% Backtracking line search
% See Algorithm 3.1 on page 37 of Nocedal and Wright
% Input Parameters :
% f: MATLAB file that returns function value
% d: The search direction
% x: previous iterate
% rho :- The backtrack step between (0,1) usually 1/2
% c: parameter between 0 and 1 , usually 10^{-4}
% Output :
% alphak: step length calculated by algorithm

alphak = 1;
[fk, gk] = feval(f,x);
xx = x;
x = x + alphak*d;
fk1 = feval(f,x);
while fk1 > fk + c*alphak*(gk'*d)
  alphak = alphak*rho;
  x = xx + alphak*d;
  fk1 = feval(f,x);
end

end

