%
%  This program computes the regularization value for 
%  CRESO, maximum curvarture,
%  and GCV and the `optimal' value, and then plots them
%  on an L-curve
%
%  SNR is the signal to noise level
%  N is the number of points to search over between zero and one
%    to find the optimal regularization parameter (b)
%  bmax is the maximum value of b allowed
%  which is 1, 2, or 3, depending on the case to be run
%
function lcurves(N,bmin_o,bmax_o,SNR,which);
%
%  for now, keep the random numbers the same from iteration to iteration
%
  randn('state',0);
%
%  convert to logarithms
%
  bmin = log10(bmin_o);
  bmax = log10(bmax_o);
%
%  Open a file to write the results to  
%
fp = fopen('lcurves.dat','wt');
figure;
hold off
hold on
if (which == 1) 
%
%================
%  Test 1
%================
%
  Q1 = [-0.2042 0.0559 0.3715 -0.2077;
         0.1347 0.2208 0.5294 0.1336;
         0.4803 0.3889 0.6862 0.4670;
         0.6078 0.6088 1.2159 0.6081];
  x1 = [1; 1; 1; 1];
  Y  = Q1*x1;
  RMS = norm(Y)/sqrt(length(Y));
  sigma = RMS/(10^(0.05*SNR));
  noise = sigma*randn(length(Y),1);
  d = Y + noise;
  Q = Q1;
%
  fprintf(fp,'=================================== \n\n');
  fprintf(fp,'\n Test Case I SNR = %f \n',SNR);  
%
%  plot the lcurve
%
  t = logspace(bmin,bmax,N);
  x = zeros(N-1,1);
  y = zeros(N-1,1);
%
  for i=2:N
     est = inv(Q'*Q+t(i)*eye(size(Q'*Q)))*Q'*d;
     x(i-1) = norm(est)^2;
     y(i-1) = norm(Q*est - d)^2;
  end;
%
  plot(log(x),log(y)); grid; title(['Test I: SNR = ' num2str(SNR) ' dB, bmin = ' num2str(bmin_o) ' bmax = ' num2str(bmax_o)]);
  xlabel('log(||g||^2)');
  ylabel('log(||Qg-d||^2');
%
%  try the gcv method
%
  b = gcv_tik(Q,d,N,bmin,bmax,fp);
  est = inv(Q'*Q+b*eye(size(Q'*Q)))*Q'*d;
  xb = norm(est)^2;
  yb = norm(Q*est - d)^2;
  plot(log(xb),log(yb),'s');
%
%  CRESO with Zero Order Tikhnov
% 
  b = creso_tik(Q,d,N,bmin,bmax,fp);
  est = inv(Q'*Q+b*eye(size(Q'*Q)))*Q'*d;
  xb = norm(est)^2;
  yb = norm(Q*est - d)^2;
  plot(log(xb),log(yb),'*');
%
%  Maximum Curvature with Zero Order Tikhnov
% 
  b = max_curve_tik(Q,d,N,bmin,bmax,fp);
  est = inv(Q'*Q+b*eye(size(Q'*Q)))*Q'*d;
  xb = norm(est)^2;
  yb = norm(Q*est - d)^2;
  plot(log(xb),log(yb),'d');
%
%  Find the `optimal' value of b
%
  b = optimal(Q,d,N,bmin,bmax,x1,fp);
  est = inv(Q'*Q+b*eye(size(Q'*Q)))*Q'*d;
  xb = norm(est)^2;
  yb = norm(Q*est - d)^2;
  plot(log(xb),log(yb),'o');
  legend('L-curve','CRESO','GCV','Curvature','Optimal');
end
if (which == 2)
%
%================
%  Test 2
%================
%
  Q2 = [-0.0505 0.0683 0.1175 0.0538 0.1709;
         0.2215 0.1729 0.0005 0.2794 0.2963;
         0.4663 0.4675 -0.0623 0.4237 0.3403;
        -0.1871 0.0701 0.2650 0.0460 0.3136;
         0.2109 0.0203 0.5130 -0.0603 0.7171];
  x2 = [1; 1; 1; 1; 1];
  Y = Q2*x2;
  RMS = norm(Y)/sqrt(length(Y));
  sigma = RMS/(10^(0.05*SNR));
  noise = sigma*randn(length(Y),1);
  d = Y + noise;
  Q = Q2;
%
  fprintf(fp,'=================================== \n\n');  
  fprintf(fp,'\n Test Case II SNR = %f \n', SNR); 
%
%  plot the lcurve
%
  t = logspace(bmin,bmax,N);
  x = zeros(N-1,1);
  y = zeros(N-1,1);
%
  for i=2:N
     est = inv(Q'*Q+t(i)*eye(size(Q'*Q)))*Q'*d;
     x(i-1) = norm(est)^2;
     y(i-1) = norm(Q*est - d)^2;
  end;
%
  plot(log(x),log(y)); grid; title(['Test II:  SNR = ' num2str(SNR) ' dB, bmin = ' num2str(bmin_o) ' bmax = ' num2str(bmax_o)]);
  xlabel('log(||g||^2)');
  ylabel('log(||Qg-d||^2');
%
%  try the gcv method
%
  b = gcv_tik(Q,d,N,bmin,bmax,fp);
  est = inv(Q'*Q+b*eye(size(Q'*Q)))*Q'*d;
  xb = norm(est)^2;
  yb = norm(Q*est - d)^2;
  plot(log(xb),log(yb),'s');
%
%  CRESO with Zero Order Tikhnov
% 
  b = creso_tik(Q,d,N,bmin,bmax,fp);
  est = inv(Q'*Q+b*eye(size(Q'*Q)))*Q'*d;
  xb = norm(est)^2;
  yb = norm(Q*est - d)^2;
  plot(log(xb),log(yb),'*');
%
%  Maximum Curvature with Zero Order Tikhnov
% 
  b = max_curve_tik(Q,d,N,bmin,bmax,fp);
  est = inv(Q'*Q+b*eye(size(Q'*Q)))*Q'*d;
  xb = norm(est)^2;
  yb = norm(Q*est - d)^2;
  plot(log(xb),log(yb),'d');
%
%  Find the `optimal' value of b
%
  b = optimal(Q,d,N,bmin,bmax,x2,fp);
  est = inv(Q'*Q+b*eye(size(Q'*Q)))*Q'*d;
  xb = norm(est)^2;
  yb = norm(Q*est - d)^2;
  plot(log(xb),log(yb),'o');
  legend('L-curve','CRESO','GCV','Curvature','Optimal');
end
if (which == 3)
%
%==================
%  Test 3
%=================
%
  Q3 = [1.001  1.967  3.055  3.985  4.991  6.001;
        7.000  8.026  8.957 10.010 11.007 12.001;
       13.000 14.010 14.982 16.007 17.004 17.998;
       19.000 20.011 20.982 22.002 29.004 30.000;
       24.998 26.005 26.990 28.008 29.004 30.000];
  x3 = [1; 1; 1; 1; 1; 1];
  Y = Q3*x3;
  RMS = norm(Y)/sqrt(length(Y));
  sigma = RMS/(10^(0.05*SNR));
  noise = sigma*randn(length(Y),1);
  d = Y + noise;
  Q = Q3;
%
  fprintf(fp,'=================================== \n\n');
  fprintf(fp,'\n  Test Case III SNR = %f \n',SNR); 
%
%  plot the lcurve
%
  t = logspace(bmin,bmax,N);
  x = zeros(N-1,1);
  y = zeros(N-1,1);
%
  for i=2:N
     est = inv(Q'*Q+t(i)*eye(size(Q'*Q)))*Q'*d;
     x(i-1) = norm(est)^2;
     y(i-1) = norm(Q*est - d)^2;
  end;
%
  plot(log(x),log(y)); grid; title(['Test III: SNR = ' num2str(SNR) ' dB, bmin = ' num2str(bmin_o) ' bmax = ' num2str(bmax_o)]);
  xlabel('log(||g||^2)');
  ylabel('log(||Qg-d||^2');
%
%  try the gcv method
%
  b = gcv_tik(Q,d,N,bmin,bmax,fp);
  est = inv(Q'*Q+b*eye(size(Q'*Q)))*Q'*d;
  xb = norm(est)^2;
  yb = norm(Q*est - d)^2;
  plot(log(xb),log(yb),'s');
%
%  CRESO with Zero Order Tikhnov
% 
  b = creso_tik(Q,d,N,bmin,bmax,fp);
  est = inv(Q'*Q+b*eye(size(Q'*Q)))*Q'*d;
  xb = norm(est)^2;
  yb = norm(Q*est - d)^2;
  plot(log(xb),log(yb),'*');
%
%  Maximum Curvature with Zero Order Tikhnov
% 
  b = max_curve_tik(Q,d,N,bmin,bmax,fp);
  est = inv(Q'*Q+b*eye(size(Q'*Q)))*Q'*d;
  xb = norm(est)^2;
  yb = norm(Q*est - d)^2;
  plot(log(xb),log(yb),'d');
%
%  Find the `optimal' value of b
%
  b = optimal(Q,d,N,bmin,bmax,x3,fp);
  est = inv(Q'*Q+b*eye(size(Q'*Q)))*Q'*d;
  xb = norm(est)^2;
  yb = norm(Q*est - d)^2;
  plot(log(xb),log(yb),'o');
  legend('L-curve','CRESO','GCV','Curvature','Optimal');
end;
%
  hold off
  fclose(fp);  
%
%  all done
%
  return;
%
%==========================
%
%  This function computes the gcv the elegant way
%
function  b = gcv_tik(Q,d,N,bmin,bmax,fp)
%
%  find the svd of Q
%
  [U,S,V]=svd(Q);
%
%  set the range of regularization parameters (t)
%
  t = logspace(bmin,bmax,N);

  
  
%
%  Now find the minimum index
%
  [dummy,i] = min(V);
  b = t(i);
%
  x = inv(Q'*Q+b*eye(size(Q'*Q)))*Q'*d;  
  fprintf(fp,'\n');
  fprintf(fp,' Elegant Method: b = %f \n',b);
  fprintf(fp,' The estimated x = %f \n',x);
%
%  all done
%
  return;
%
%================================  
%
%  This function used creso to compute b
%
function  b = creso_tik(Q,d,N,bmin,bmax,fp)
%

%
  t = logspace(bmin,bmax,N);
%
  
%
%  Now find the maximum index
%
  [dummy,i] = max(creso);
  b = t(i);
%
  x = inv(Q'*Q+b*eye(size(Q'*Q)))*Q'*d;  
  fprintf(fp,'\n');
  fprintf(fp,' CRESO Method: b = %f \n',b);
  fprintf(fp,' The estimated x = %f \n',x);
%
%  all done
%
  return;
%
%
%================================  
%
%  This function used the maximum curvature to compute b
%
function  b = max_curve_tik(Q,d,N,bmin,bmax,fp)
%

%  set the range of regularization parameters (t)
%
  t = logspace(bmin,bmax,N);

%
%  Now find the maximum index
%
   [dummy,i] = max(curve);
   b = t(i);
%
  x = inv(Q'*Q+b*eye(size(Q'*Q)))*Q'*d;  
  fprintf(fp,'\n');
  fprintf(fp,' Maximum Curvature Method: b = %f \n',b);
  fprintf(fp,' The estimated x = %f \n',x);
%
%  all done
%
  return;
%
%========================================
%
%  This routine computes the optimal value of b
%  For now, the optimal value is the one which produces
%  the minimin error in any component
%
function b = optimal(Q,d,N,bmin,bmax,x,fp);
%
  t = logspace(bmin,bmax,N);
  V = zeros(N,1);
  for i=2:N
     est = inv(Q'*Q+t(i)*eye(size(Q'*Q)))*Q'*d;
     V(i) = max(abs(est-x));
  end;
  V(1) = max(V);
  [dummy,i]=min(V);
  b = t(i);
  est = inv(Q'*Q+b*eye(size(Q'*Q)))*Q'*d;
  fprintf(fp,'\n');
  fprintf(fp,' Optimal: b = %f \n',b);
  fprintf(fp,' The estimated x = %f \n',est);
%
%  all done
%
  return;
