inv_solve_diff_kalman

 INV_SOLVE_DIFF_KALMAN inverse solver for difference EIT
 img= inv_solve_diff_kalman( inv_model, data1, data2)

 img        => output image
 inv_model  => inverse model struct
 data1      => differential data at earlier time
 data2      => differential data at later time

 if inv_model.fwd_model.stimulation(:).delta_time
   exists and is non_zero, then the kalman filter will
   be applied to each data measurement separately

 Note that the classic Kalman filter assumes that the
   time step between each measurement is constant
   (ie as part of the state update eqn). inv_solve_diff_kalman
   cannot work with non-constant time steps

 if inv_model.inv_solve_diff_kalman.keep_K_k1 = 1
  then img outputs img.inv_solve_diff_kalman.K_k1 = K_k1
  this can be used to estimate noise properties

0001 function img= inv_solve_diff_kalman( inv_model, data1, data2)
0002 % INV_SOLVE_DIFF_KALMAN inverse solver for difference EIT
0003 % img= inv_solve_diff_kalman( inv_model, data1, data2)
0004 %
0005 % img        => output image
0006 % inv_model  => inverse model struct
0007 % data1      => differential data at earlier time
0008 % data2      => differential data at later time
0009 %
0010 % if inv_model.fwd_model.stimulation(:).delta_time
0011 %   exists and is non_zero, then the kalman filter will
0012 %   be applied to each data measurement separately
0013 %
0014 % Note that the classic Kalman filter assumes that the
0015 %   time step between each measurement is constant
0016 %   (ie as part of the state update eqn). inv_solve_diff_kalman
0017 %   cannot work with non-constant time steps
0018 %
0019 % if inv_model.inv_solve_diff_kalman.keep_K_k1 = 1
0020 %  then img outputs img.inv_solve_diff_kalman.K_k1 = K_k1
0021 %  this can be used to estimate noise properties
0022  
0023 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0024 % $Id: inv_solve_diff_kalman.m 3661 2012-11-20 21:18:01Z bgrychtol $
0025 
0026 fwd_model= inv_model.fwd_model;
0027 pp= fwd_model_parameters( fwd_model );
0028 
0029 img_bkgnd= calc_jacobian_bkgnd( inv_model );
0030 J = calc_jacobian(img_bkgnd);
0031 
0032 RtR = calc_RtR_prior( inv_model );
0033 Q   = calc_meas_icov( inv_model );
0034 hp  = calc_hyperparameter( inv_model );
0035 
0036 if isfield(fwd_model.stimulation(1),'delta_time')
0037    delta_time= [fwd_model.stimulation(:).delta_time];
0038    if diff(delta_time) ~= 0;
0039       error('All time steps must be same for kalman filter');
0040    end
0041 else
0042    delta_time=0;
0043 end
0044 
0045 % sequence is a vector location of each stimulation in the frame
0046 if delta_time == 0
0047    sequence = size(J,1);
0048 else
0049    for i=1:length(fwd_model.stimulation)
0050       sequence(i) = size(fwd_model.stimulation(i).meas_pattern,1);
0051    end
0052    sequence= cumsum( sequence );
0053 end
0054 
0055 
0056 if pp.normalize
0057    dva= data2 ./ data1 - 1;
0058 else   
0059    dva= data2 - data1;
0060 end
0061 
0062 [sol, K_k1] = kalman_inv( J, hp^2*Q, RtR, dva, sequence );
0063 
0064 % create a data structure to return
0065 img.name= 'solved by inv_solve_diff_kalman';
0066 img.elem_data = sol;
0067 img.inv_model= inv_model;
0068 img.fwd_model= fwd_model;
0069 
0070 try % keep parameter if requested
0071    if inv_model.inv_solve_diff_kalman.keep_K_k1
0072       img.inv_solve_diff_kalman.K_k1 = K_k1;
0073    end
0074 end
0075 
0076 % Kalman filter - estimates x where z_k = H_k*x_k + noise
0077 % J - Jacobian NxM
0078 % RegM - Regularization on the Measurements
0079 % RegI - Regularization on the Image
0080 % y is vector of measurements
0081 % seq is sequence vector
0082 %
0083 % K_k1 is the linearized reconstruction matrix for
0084 %   the final step. It can be used to estimate noise
0085 %   properties of the algorithm
0086 %
0087 function [x, K_k1]= kalman_inv( J, RegM, RegI, y, seq);
0088 %x = (J'*RegM*J + RegI )\J'*RegM*y; return;
0089 %Notation x_k1_k is x_{k+1|k}
0090 
0091 % n is nmeas, m is ndata
0092 [m,n]=  size(J);
0093 
0094 % F is the state transition matrix (I for random walk)
0095 F_k= speye(n);
0096 % Q is state noise covariance (model with I)
0097 Q_k= RegI;
0098 
0099 % Initial error covariance estimate.
0100 C_k1_k1= speye(n);
0101 
0102 % mean x_priori image - assume 0
0103 x0= zeros(n,1);
0104 x_k1_k1= x0;
0105 
0106 ll= size(y,2);
0107 x= zeros(n,ll*length(seq));
0108 
0109 seq= [0;seq(:)];
0110 iter=0;
0111 for i=1:ll
0112    for ss= 2:length(seq);
0113       eidors_msg('inv_solve_diff_kalman: iteration %d.%d',i,ss-1,2);
0114 
0115       seq_i= (seq(ss-1)+1) : seq(ss);
0116 
0117 % The Kalman filter doesn't need the regularization at all
0118       H_k1= J(seq_i,:);
0119       yi= y(seq_i,i);
0120       G = RegM(seq_i,seq_i);
0121       [x_k1_k1, C_k1_k1, K_k1] = ...
0122              kalman_step( x_k1_k1, C_k1_k1, ...
0123                           H_k1, yi, F_k, Q_k, G );
0124       iter=iter+1;
0125       x(:,iter) = x_k1_k1;
0126    end
0127 end
0128 
0129 function [x_k1_k1, C_k1_k1, K_k1] = ...
0130                   kalman_step( x_k_k, C_k_k, ...
0131                                H_k1, yi, F_k, Q_k, G )
0132    n= size(H_k1,2);
0133 
0134    % Prediction
0135    x_k1_k = F_k * x_k_k;
0136    C_k1_k = F_k * C_k_k * F_k' + Q_k;
0137    % Correction
0138    HCHt   = H_k1 * C_k1_k * H_k1';
0139    K_k1   = C_k1_k * H_k1' / (HCHt + G);
0140    yerr   = yi - H_k1 * x_k1_k;
0141    x_k1_k1= x_k1_k + K_k1 * yerr; 
0142    C_k1_k1= (speye(n) - K_k1 * H_k1) * C_k1_k;
0143    
0144 
0145 function x= kalman_inv_cgls( J, RegM, RegI, y);
0146    [m,n]=  size(J);
0147    Rx0= zeros(n,1);
0148    H= [chol(W)*J;RegI];
0149    b= [y;Rx0];
0150    x= H\b;
0151 
0152 % Adapted from code in Hansen's regularization tools
0153 function x= cg_ls_inv0( J, R, y, Rx0, maxiter, tol )
0154 %  A = [J;R];
0155    b=[y;Rx0];
0156    [m,n]= size(J);
0157    m_idx= 1:m; n_idx = m+(1:n);
0158    Jt = J.';
0159    Rt = R.';
0160    x = zeros(n,1); 
0161 %  d = A'*b;
0162    d = Jt*b(m_idx) + Rt*b(n_idx);
0163    r = b; 
0164    normr2 = d'*d; 
0165     
0166    k=0; % Iterate.
0167    x_delta_filt= 1; x_filt= .1;
0168    while 1 
0169      % Update x and r vectors.
0170 %    Ad = A*d;
0171      Ad = [J*d;R*d];
0172      Alpha = normr2/(Ad'*Ad); 
0173      xpre= x;
0174      x  = x + Alpha*d; 
0175 
0176      k= k+1; if k==maxiter; break ; end
0177 
0178      x_delta= norm(xpre-x)/norm(x);
0179      x_delta_filt= x_delta_filt*(1-x_filt) + x_filt*x_delta;
0180 
0181      if x_delta_filt<tol; break ; end
0182 
0183      r  = r - Alpha*Ad; 
0184 %    s  = A'*r;
0185      s  = Jt*r(m_idx) + Rt*r(n_idx);
0186     
0187      % Update d vector.
0188      normr2_new = s'*s; 
0189      Beta = normr2_new/normr2; 
0190      normr2 = normr2_new; 
0191      d = s + Beta*d; 
0192       
0193    end 
0194 %     disp([k, x_delta, x_delta_filt]);
0195 
0196