% inverse input problem
% PDE model

clear
plotopts={'linewidth', 1.0};
sigma_noise = 0.005;
rng(129);  % seed random number generator - always use the same one


dx=0.05;nx=21;x=0:dx:1;
tf=2;nt=51;dt=tf/(nt-1);
t=0:dt:tf;

% two sensor locations
C=zeros(2,nx);
C(1,2)=1;
C(2,11)=1; 
b=10;
fo=dt/(b*dx*dx);A=(1+2*fo)*eye(nx);
for i=2:nx-1
   A(i,i+1)=-fo;A(i,i-1)=-fo;
end
A(1,2)=-2*fo;A(nx,nx-1)=-2*fo;

% the input heat flux
% to be determined
for k=1:nt
   if k<6
      u(k)=k-1;
   elseif k<16
      u(k)=11-k;
   elseif k<21
      u(k)=-21+k;
   else
      u(k)=0;
   end     
end
uex=u;

% the output data  
% with an additive noise
Tw=zeros(nx,1);Y(:,1)=C*Tw;
for k=2:nt
   Tw(1)=Tw(1)+2*fo*dx*u(k);
   Tnew=A\Tw;Y(:,k)=C*Tnew;
   Tw=Tnew;
end
noise=sigma_noise*randn(size(Y));
Y=Y+noise; 
eps=norm(noise)^2

% conjugate gradient algorithm

u=zeros(nt,1);lsi=1e6;it=1;

while lsi>eps
   
   % the output model
   % and the LS-criterion
  Tw=zeros(nx,1);Y(:,1)=C*Tw;
  for k=2:nt
     Tw(1)=Tw(1)+2*fo*dx*u(k);
      Tnew=A\Tw;YM(:,k)=C*Tnew;
      Tw=Tnew;
   end
   e=Y-YM;
   lsi=norm(e)^2;
   ls(it)=lsi;
   
   % the adjoint variable
   % and the gradient
   Pw=zeros(nx,1);g=Pw;
   g(nt)=0;
   for k=nt:-1:2
      Pw=Pw+2*dt*C'*e(:,k-1);
      Pnew=A\Pw; Pw=Pnew;
      g(k-1)=Pw(1);
   end
  ng=norm(g)^2;
   
   % the descent direction
   if it==1
      w=-g;
   else
      w=-g+ng/ng1*w1;
   end
   
   % the descent step length
   % quadratic case
   dTw=zeros(nx,1);dY(:,1)=C*dTw;
   for k=2:nt
      dTw(1)=dTw(1)+2*fo*dx*w(k);
      dTnew=A\dTw;dY(:,k)=C*dTnew;
      dTw=dTnew;
   end
   
   ro=g'*w/norm(dY)^2;
   fprintf('iter %2d  S=%g  ro=%f\n', it, lsi,ro)

   
   % the new iterate
   
   u=u+ro*w;
   g1=g;ng1=ng;w1=w;
   it=it+1;
end

figure(1); clf;
semilogy(ls, plotopts{:})
grid on
xlabel('iteration'); ylabel('S');

figure(2); clf;
plot(t,uex,'k-',t,u,'ro', plotopts{:})
legend('Exact', 'estimate')
xlabel('time t'); ylabel('Heat flux u(t)')

figure(3); clf;
plot(t,Y,'k-',t,YM,'ro', plotopts{:})
legend('Exact', 'Exact', 'estimate', 'estimate')
xlabel('time t'); ylabel('Temperature output Y(t)');
text(0.5,0.35,'x_1=0.05');
text(0.75, 0.1, 'x_2=0.5');
 

% if( sigma_noise == 0.005 )
%     fig=figure(3);
%     print(fig,fullfile(pwd, 'Figure6_44.tif'),'-dtiff','-r600')
%     fig=figure(2);
%     print(fig,fullfile(pwd, 'Figure6_45.tif'),'-dtiff','-r600')
%     fig=figure(1);
%     print(fig,fullfile(pwd, 'Figure6_48.tif'),'-dtiff','-r600')
% else
%     fig=figure(3);
%     print(fig,fullfile(pwd, 'Figure6_46.tif'),'-dtiff','-r600')
%     fig=figure(2);
%     print(fig,fullfile(pwd, 'Figure6_47.tif'),'-dtiff','-r600')
%     fig=figure(1);
%     print(fig,fullfile(pwd, 'Figure6_49.tif'),'-dtiff','-r600')
% end