% inverse parameter problem
% ODE modelling

clear global
global xgrid Bet Betmin u nn dt Y YM e dz w

xgrid=[-2; 0.; 0.1; 0.2; 0.3; 0.4; 0.5;1.0];
p=8;Ip=eye(p);
Betexact=[0.25;0.25;0.25;0.4;1.3;0.4;0.25;0.25];
Betmin=0.15;Bet=Betexact;

nn=101;dt=0.02;tc=1.;
u0=-(1+0.25+(sqrt(2*pi))/15)/tc;
t=0:dt:2;Y=zeros(1,nn);w=zeros(1,p);
% compute the input
for n=1:nn
   if t(n)<tc
      u(n)=u0;
   else
      u(n)=0;
   end   
end

% compute the output
ls0=phinl(0); Y=YM;

% conjugate gradient algorithm

Bet(2:p)=0.4;

itmx=20;

for iter=1:itmx
   
   % the ls-criterion
   lsit=phinl(0);
   ls(iter)=lsit;
   
   %  the adjoint 
   
   pw=[0; 0];g=zeros(p,1);psi(nn)=0;
   for n=nn:-1:2
      T0=max(YM(n-1),-2);
      aw=interp1(xgrid,Bet,T0);
      AT=[1+dt -dt/aw;-dt 1+dt/aw];
      pw(2)=pw(2)-dt*e(n-1);pnew=AT\pw;
      psi(n-1)=-pnew(2)*dz(n)*dt;
   end       
   pw=pnew;
 
    % the gradient
    for n=1:nn
       YM(n)=max(YM(n),-2);
    end
        
    sigmi=interp1(xgrid,Ip,YM);
    g=psi*sigmi; ng=norm(g)^2;

% the descent direction

   if mod(iter,p)==1
      w=-g;
   else
      w=-g+ng/ng1*w1;
    end
   
% the descent length (line search)
   
   ro=fminbnd(@phinl,0,100);
   fprintf('iter %2d  S=%g  ro=%f\n', iter, lsit,ro)

   Bet=Bet+ro*w';
   for i=2:p
      Bet(i)=max(Bet(i),Betmin);
   end
      g1=g;w1=w;ng1=ng;

end
figure(1); clf
plot(xgrid,Betexact,'k-',xgrid,Bet,'ro','linewidth', 1.0)
xlabel('x'); ylabel('a(x)');
legend('exact','estimate','location','northwest')
figure(2); clf;
plot(t,Y,'k-',t,YM,'ro','linewidth', 1.0)
xlabel('x'); ylabel('Temperature');
legend('exact','estimate')
figure(3); clf;
semilogy([1:itmx],ls,'k-','linewidth', 1.0)
xlabel('iteration'); ylabel('S');
grid on

% if( Betmin == 0.25 )
%     fig=figure(2);
%     print(fig,fullfile(pwd, 'Figure6_30.tif'),'-dtiff','-r600')
%     fig=figure(1);
%     print(fig,fullfile(pwd, 'Figure6_31.tif'),'-dtiff','-r600')
%     fig=figure(3);
%     axis( [1 itmx 1e-4 1] )
%     print(fig,fullfile(pwd, 'Figure6_32.tif'),'-dtiff','-r600')
% else
%     fig=figure(1);
%     print(fig,fullfile(pwd, 'Figure6_33.tif'),'-dtiff','-r600')
%     fig=figure(3);
%     axis( [1 itmx 1e-3 1] )
%     print(fig,fullfile(pwd, 'Figure6_34.tif'),'-dtiff','-r600')
% end
    


function ls=phinl(r)

global xgrid Bet Betmin u nn dt Y YM e dz w

p=size(Bet);
Betr=Bet+r*w';
 for i=1:p
      Betr(i)=max(Betr(i),Betmin);
 end
zw=[1; 1];YM(1)=1.;
   for n=2:nn
      T0=max(zw(2),-2);
      aw=interp1(xgrid,Betr,T0);
      A=[1+dt -dt;-dt/aw 1+dt/aw]; 
      zw(1)=zw(1)+dt*u(n);znew=A\zw; 
      YM(n)=znew(2);zw=znew;
      dz(n)=(zw(1)-zw(2))/aw^2;
   end
e=Y-YM ; ls=norm(e)^2;
end