% non quadratic case- conjugate gradient
clear
 
% global Y B w
xx=0:0.1:4; yy=xx;w=[0; 0];
Y=[0.66;1.00;1.20;1.33;1.43];



for i=1:41
    for j=1:41
      B=[xx(i);yy(j)]; 
      q(j,i)=phi2(0,B,w,Y);
   end
end
figure(1); clf;
contour(xx,yy,q,40)
xlabel('B_1')
ylabel('B_2')
hold on
clear xx; clear yy   
xtickformat('%0.1f')
ytickformat('%0.1f') 

hold on
clear xx; clear yy ;  

B=[4; 0];
xx(1)=B(1);yy(1)=B(2);
for k=2:20
   
   ls(k-1)=phi2(0,B,w,Y);
   
   g= [0;0];
   for n=1:5
      e(n)=Y(n)-B(1)*n/(B(2)+n);
      g(1)=g(1)-e(n)*n/(B(2)+n);
      g(2)=g(2)+e(n)*B(1)*n/(B(2)+n)^2;
   end
   ng=g'*g;
   
  if mod(k,2)==0
      w=-g;
   else
      w=-g+ng/ng1*w1;
  end
   
   linesearch = @(r) phi2(r,B,w,Y);
   ro=fminbnd(linesearch,0,100);
   
 B=B+ro*w;
 g1=g;w1=w;ng1=ng;
 xx(k)=B(1);yy(k)=B(2);
    
end

plot(xx,yy,'lineWidth',2)

figure(2); clf;
semilogy(ls, 'linewidth', 1.0)
xlabel('iteration');
ylabel('S=(Y-XB)^T(Y-XB)')
ytickformat('%0.2f')
grid on


% fig=figure(2);
% print(fig,fullfile(pwd, 'Figure6_14.tif'),'-dtiff','-r600')
% 
% fig=figure(1)
% print(fig,fullfile(pwd, 'Figure6_15.tif'),'-dtiff','-r600')



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function S=phi2(r,B,w,Y)

% global Y B w

bw=B+r*w;

for n=1:5
   e(n)=Y(n)-bw(1)*n/(bw(2)+n);
end 

S=0.5*e*e';

end