%
%  This program implements first order Tikhonov regularization
%  with a constant regularization parameter
%  for solving the inverse problem of electrocardiography.
%  This version is for the endocardial data
%
%  Outputs: 
%      est = the estimated endocardial electrograms
%  Inputs:
%      mu = the constant regularization parameter
%
function est = probe_tikk(mu);
%  
%  read in the necessary data
%
%  the transfer matrix
%
  load ZHB.ASCII
  Z = ZHB;
%
%  TMAT selects those points on the endocardium that correspond to 
%  where we have sensors
%
  load TMAT
%
%  load the gradient matrix for first order regularization
%
  load ASCIIKGRAD
  R = sqrtm(ASCIIKGRAD);
  R = real(R);
%
%  load in the endocardial data measured at the sensors
%
  load endo_data
  endo = endo_data';
%  
%  load in the probe data
%
  load probe_data
  probe = probe_data';
%
%  determine some of the parameters we need
%
  N_electrodes = 66;
  N = 201;
%
%  compute the gsvd of Z and R
%
  [U,V,X,C,S] = gsvd(Z,R);
  [m,p] = size(Z);
  temp_s = diag(S);
  ts = temp_s(p-m+1:p);
  tc = diag(C,p-m);
  s = ts(1:m);
  c = tc(1:m);
  s2 = s.*s;
  c2 = c.*c;
%
  Xin = inv(X');
  Xin_trunc = Xin(:,p-m+1:p);
%
  est = [];
%
%  compute the endocardial estimate at each time
%
  for time = 1:N
    up = U'*probe(:,time);
    top = c.*up;
    bot = c2 + mu*s2;
    th = top./bot;
    full_est = Xin_trunc*th;
    est(:,time) = TMAT*full_est;
  end;
%
%  Now make the epicardial estimate plot
% 
  time = [1:N];
  ave_cc = plot_15(est,endo,time,mu);
  fprintf(' The average correlation coefficient for endoardial electrograms is %f \n', ave_cc);
%  
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Plot the estimated endocardial electrograms
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
function ave_cc = plot_15(ya,ytrue,time,mu);
%
%  
  scale = 1.2;
%
%  plot the data
%
  figure;
  fs = 14;
%   whitebg(1,[0 0 0]);
%   set(1,'Color',[0 0 0]);
%   set(1,'InvertHardCopy','off');
%
  orient landscape;
%
  total_cc = [];
  for k=1:4
      p = k;
      for j=1:4
         ii = (k-1)*4+j;
         if (ii ~= 16)
           ymin = scale*min(ytrue(ii,:));
           ymax = scale*max(ytrue(ii,:));
           subplot(4,4,p);
           axis([time(1) time(end) ymin ymax]);
           hold on
           plot(time,ya(ii,:),'b-','Linewidth',1.5);
           plot(time,ytrue(ii,:),'r-','Linewidth',1.5);
           axis off;
           cc = corrcoef(ya(ii,:),ytrue(ii,:));
           total_cc = [total_cc cc(1,2)];
           title([num2str(cc(1,2),2)]);
           hold off
         end;
         p = p+4;
      end;
  end;  
  ave_cc = mean(total_cc);
  text(10,-12,['\mu  = ' num2str(mu)]);
  text(10,-15,[' ave cc = ' num2str(ave_cc)]);
  print -deps2 endo.eps
  return;
