% lininva.m
%
% Solve M604D HW07 Question A
%
% Philip D. Loewen
% 2 May 2001

saveplots = 0;      % 0-1 variable: Make *.eps files (y/n)?
format compact; format long;

% Initialize constants specified on problem sheet
a = 1/4; b = 1/2;   % Interval endpoints
c = 1;              % Exponential decay rate
m = 9;              % Problem dimension

% PART A(b):  Generate Data Vector d
disp(' ');
disp('---------------------------------------------------------');
disp(' Part A(b):  Calculate data vector.');
disp('---------------------------------------------------------');
disp(' ');

k = (1:m)';
d = ( gamma(k) ./ (c*ones(size(k))).^k ) ...
     .* (gammainc(c*b,k)-gammainc(c*a,k))
  
% PART A(c):  Solve the Minimum-Norm Problem assuming data are exact.

disp(' ');
disp('---------------------------------------------------------');
disp(' Part A(c):  Fit data exactly with minimum-norm function.');
disp('---------------------------------------------------------');
disp(' ');

% Build matrix Gamma = G * G' (generalized interpretation of G')
ss = 1:(2*m-1);
ggtvals = gamma(ss) .* gammainc(2*c,ss) ./ (2*c*ones(size(ss))).^ss;
J = ones(m,1)*(1:m) + (1:m)'*ones(1,m) - 1;
GGT = ggtvals(J);

% Solve linear system to get minimum-norm solution in absence of noise
lambda = GGT\d(1:m);
disp(['Minimum-norm solution assuming data exact is  x=G''*lambda,']);
disp(['where lambda is given below:',sprintf('\n  %13.5e',lambda')]);

disp(' ');
disp(['This x has N(x) = ',sprintf('%11.3e',0.5*lambda'*GGT*lambda),',']);
mv = GGT*lambda - d;                  % Misfit vector Gx-d
disp(['           E(x) = ',sprintf('%11.3e',0.5*mv'*mv),'.']);

% Plot minimum-norm solution
numt = 201;                           % Number of t-values to plot with
tt = linspace(0,1,numt);              % Vector of t-values
ttmtx = ones(m,1)*tt;                 % Big matrix of t-values, all rows equal
kkmtx = (1:m)'*ones(1,numt);          % Big matrix of k-values, all cols equal
GT = ttmtx.^(kkmtx-1) .* exp(-c*ttmtx); % Row j gives values of fcn gj(t)
xx = lambda'*GT;                      % Values of function G'*lambda
figure('Name',['Exact data fitting (mu=infinity)']);
xc = plot(tt,xx);
hold on;
xe = plot([0,a,a,b,b,1],[0,0,1,1,0,0],':');
title('Exact data fitting (mu=infinity)');
xlabel('t');
ylabel('x');
legend([xc,xe],'Reconstruction','Exact');
set(gca,'YLim',[-0.4,1.2]);

if saveplots,
  set(gcf,'PaperPosition',[1,1,4,3]);
  disp(['Figure number ',num2str(gcf),' saved as minnormx.eps.'])
  print -deps minnormx.eps
end;

% PART A(d):  Modify approach to account for possible data errors.
disp(' ');
disp('---------------------------------------------------------');
disp(' Part A(d):  Tradeoff between E and N using various mu.');
disp('---------------------------------------------------------');
disp(' ');


% Build covariance matrix S:  important to *square* the std dev vals.
S    = diag([abs(d)/10].^2);
Sinv = diag(1./[abs(d)/10].^2);

nummu = 30;
muENc = zeros(4,nummu+1);

disp(['    k   mu          cond(LHS)   E           N     ',...
     '      |(E,N)|']);
k=1;
for mu=sort([logspace(-5,3,nummu),0.1873357]),
  mtx = S/mu + GGT;
  cn  = cond(mtx);
  lambda = (S/mu + GGT)\d(1:m);
  N   = 0.5*lambda'*GGT*lambda;
  mv  = GGT*lambda - d;                  % Misfit vector Gx-d
  E   = 0.5*mv'*Sinv*mv;
  disp([sprintf('   %2.0f',k),...
        sprintf(' %11.3e',[mu,cn,E,N,sqrt(E*E+N*N)])]);
  muENc(:,k) = [mu;E;N;cn];
  k = k+1;
end
disp(' ');

figure('Name','Tradeoff Curve');
plot(muENc(2,:),muENc(3,:),'r','LineWidth',2);
xlabel('E(x)');
ylabel('N(x)');

if saveplots,
  set(gcf,'PaperPosition',[1,1,4,3]);
  disp(['Figure number ',num2str(gcf),' saved as tradeoff.eps.'])
  print -deps tradeoff.eps
end;

figure('Name','Condition Number vs Regularization');
plot(muENc(1,:),muENc(4,:),'r','LineWidth',2);
xlabel('\mu');
ylabel('cond(LHS)');

if saveplots,
  set(gcf,'PaperPosition',[1,1,4,3]);
  disp(['Figure number ',num2str(gcf),' saved as regcond.eps.'])
  print -deps regcond.eps
end;

% PART A(e):  Choose mu to hit target misfits of 8 and 10.
disp(' ');
disp('---------------------------------------------------------');
disp(' Part A(e):  Choose mu to hit misfit targets.');
disp('---------------------------------------------------------');

figure('Name',['Misfit reconstructions with 2E = 8 and 10']);
crvhndls = plot([0,a,a,b,b,1],[0,0,1,1,0,0],':');
hold on;
title(['Misfit reconstructions with 2E = 8 and 10']);
xlabel('t');
ylabel('x');

for mftar = [8,10],  % do all this stuff for each target.
  mulo   = 8.e-4;
  lambda = (S/mulo + GGT)\d(1:m);
  mv     = GGT*lambda - d;         % Data error vector Gx-d
  mflo   = mv'*Sinv*mv;            % Weighted norm of data error

  muhi   = 3.e-3;
  lambda = (S/muhi + GGT)\d(1:m);
  mv     = GGT*lambda - d;         % Data error vector Gx-d
  mfhi   = mv'*Sinv*mv;            % Weighted norm of data error

  it = 1;
  while muhi-mulo>0.0001*mulo & it<20,
    mu = (mulo+muhi)/2;
    lambda = (S/mu + GGT)\d(1:m);
    mv  = GGT*lambda - d;         % Data error vector Gx-d
    mf  = mv'*Sinv*mv;            % Weighted norm of data error
%    disp(['it=',sprintf('%2.0f',it),...
%        ': [mulo,muhi]=',sprintf('[%11.3e,%11.3e]',[mulo,muhi]),...
%        ', [mflo,mfhi]=',sprintf('[%8.4f,%8.4f]',[mflo,mfhi]),'.']);
    if mf>mftar,
      mulo = mu; mflo = mf;
    else
      muhi = mu; mfhi = mf;
    end;
    it = it+1;
  end;
  mu = (mulo+muhi)/2;
  disp(' ');
  disp(['For target misfit of ',num2str(mftar),...
      sprintf(' use mu = %12.4e.',mu)]);
  lambda = (S/mu + GGT)\d(1:m);
  mv  = GGT*lambda - d;         % Data error vector Gx-d
  mf  = mv'*Sinv*mv;            % Weighted norm of data error
  disp(['Actual misfit then equals',...
      sprintf(' %12.4e, with lambda shown below:',mf),...
      sprintf('\n   %12.4e',lambda')]);
  
  xx = lambda'*GT;                      % Values of function G'*lambda
  xc = plot(tt,xx);
  crvhndls = [crvhndls,xc];

end;

set(crvhndls(2),'LineStyle','--');
set(gca,'YLim',[-0.4,1.2]);
legend(crvhndls,'Exact','Misfit 8','Misfit 10');
if saveplots,
  set(gcf,'PaperPosition',[1,1,4,3]);
  disp(['Figure number ',num2str(gcf),' saved as misfits.eps.'])
  print -deps misfits.eps
end;

% PART A(f):  Choose mu to hit corner of tradeoff curve.
disp(' ');
disp('---------------------------------------------------------');
disp(' Part A(f):  Choose mu to hit corner of tradeoff curve.');
disp('---------------------------------------------------------');

mulo   = 8.e-4;
lambda = (S/mulo + GGT)\d(1:m);
mv     = GGT*lambda - d;         % Data error vector Gx-d
Elo    = 0.5*mv'*Sinv*mv;        % Weighted norm of data error
Nlo    = 0.5*lambda'*GGT*lambda;
Flo    = mulo*Nlo-Elo;

muhi   = 3.e-1;
lambda = (S/muhi + GGT)\d(1:m);
mv     = GGT*lambda - d;         % Data error vector Gx-d
Ehi    = 0.5*mv'*Sinv*mv;            % Weighted norm of data error
Nhi    = 0.5*lambda'*GGT*lambda;
Fhi    = muhi*Nhi-Ehi;

it = 1;
while Fhi-Flo>0.0001*Flo & it<20,
  mu = (mulo+muhi)/2;
  lambda = (S/mu + GGT)\d(1:m);
  mv  = GGT*lambda - d;         % Data error vector Gx-d
  E   = 0.5*mv'*Sinv*mv;        % Weighted norm of data error
  N   = 0.5*lambda'*GGT*lambda;
  F   = mu*N-E;
  disp(['it=',sprintf('%2.0f',it),...
      ': [mulo,muhi]=',sprintf('[%11.3e,%11.3e]',[mulo,muhi]),...
      ', [Flo,Fhi]=',sprintf('[%8.4f,%8.4f]',[Flo,Fhi]),'.']);
  if F<0,
    mulo = mu; Flo = F;
  else
    muhi = mu; Fhi = F;
  end;
  it = it+1;
end;
mu = (mulo+muhi)/2;
disp(' ');
disp(['For corner of tradeoff curve',...
    sprintf(' use mu = %12.4e.',mu)]);
lambda = (S/mu + GGT)\d(1:m);
mv  = GGT*lambda - d;         % Data error vector Gx-d
E   = 0.5*mv'*Sinv*mv;            % Weighted norm of data error
N   = 0.5*lambda'*GGT*lambda;
disp(['Then |(E,N)| equals',...
    sprintf(' %12.4e, with lambda shown below:',sqrt(E*E+N*N)),...
    sprintf('\n   %12.4e',lambda')]);

figure('Name','Reconstruction from corner of tradeoff curve');
hold on;
title('Reconstruction from corner of tradeoff curve');
plot([0,a,a,b,b,1],[0,0,1,1,0,0],':');
set(gca,'YLim',[-0.4,1.2]);
xlabel('t');
ylabel('x');

xx = lambda'*GT;                      % Values of function G'*lambda
xc = plot(tt,xx);

if saveplots,
  set(gcf,'PaperPosition',[1,1,4,3]);
  disp(['Figure number ',num2str(gcf),' saved as cornfit.eps.'])
  print -deps cornfit.eps
end;

figure('Name','Corner Point'); hold on;
title('Tradeoff Point Nearest the Origin'); 
xlabel('E');
ylabel('N');
plot(muENc(2,:),muENc(3,:),'r');
plot(muENc(2,:),muENc(3,:),'k+');
axis equal
set(gca,'XLim',0.1*[0,1],'YLim',0.1*[0,1]);
th = linspace(0,pi/2,90);
rr = 8.8877e-002;
plot(rr*cos(th),rr*sin(th),'g');
zoom on;

if saveplots,
  set(gcf,'PaperPosition',[1,1,4,3]);
  disp(['Figure number ',num2str(gcf),' saved as cornerpt.eps.'])
  print -deps cornerpt.eps
end;
