time_prior_solve

 TIME_PRIOR_SOLVE inverse solver to account for time differences
 img= aa_inv_solve( inv_model, data1, data2)
 img        => output image (or vector of images)
 inv_model  => inverse model struct
 data1      => differential data at earlier time
 data2      => differential data at later time

 both data1 and data2 may be matrices (MxT) each of
  M measurements at T times
 if either data1 or data2 is a vector, then it is expanded
  to be the same size matrix

0001 function img= time_prior_solve( inv_model, data1, data2)
0002 % TIME_PRIOR_SOLVE inverse solver to account for time differences
0003 % img= aa_inv_solve( inv_model, data1, data2)
0004 % img        => output image (or vector of images)
0005 % inv_model  => inverse model struct
0006 % data1      => differential data at earlier time
0007 % data2      => differential data at later time
0008 %
0009 % both data1 and data2 may be matrices (MxT) each of
0010 %  M measurements at T times
0011 % if either data1 or data2 is a vector, then it is expanded
0012 %  to be the same size matrix
0013 
0014 % TODO: This function really should be calling the proper
0015 %   prior calculator functions, and not reimplementing
0016 %   them internally
0017 
0018 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0019 % $Id: time_prior_solve.html 2819 2011-09-07 16:43:11Z aadler $
0020 
0021 fwd_model= inv_model.fwd_model;
0022 time_steps = inv_model.time_prior_solve.time_steps;
0023 l_ts  = time_steps*2 + 1;
0024 
0025 % The one_step reconstruction matrix is cached
0026 one_step_inv = eidors_obj('get-cache', inv_model, 'time_prior_solve');
0027 if ~isempty(one_step_inv)
0028     eidors_msg('time_prior_solve: using cached value', 2);
0029 else
0030     img_bkgnd= calc_jacobian_bkgnd( inv_model );
0031     J = calc_jacobian( fwd_model, img_bkgnd);
0032 
0033 %   one_step_inv= standard_form( inv_model, J );
0034     one_step_inv= data_form( inv_model, J );
0035 
0036     eidors_obj('set-cache', inv_model, 'time_prior_solve', one_step_inv);
0037     eidors_msg('time_prior_solve: setting cached value', 2);
0038 end
0039 
0040 dva = calc_difference_data( data1, data2, inv_model.fwd_model);
0041 
0042 l_dva = size( dva, 2);
0043 
0044 idx= [-time_steps:time_steps]'*ones(1,l_dva) + ...
0045      ones(l_ts,1)*(1:l_dva);
0046 % replicate first and last measurements
0047 idx(idx<1) = 1;
0048 idx(idx>l_dva) = l_dva;
0049 
0050 dvat= reshape(dva(:,idx),[],l_dva);
0051  
0052 sol = one_step_inv * dvat;
0053 
0054 % create a data structure to return
0055 img.name= 'solved by time_prior_solve';
0056 img.elem_data = sol;
0057 img.inv_model= inv_model;
0058 img.fwd_model= fwd_model;
0059 
0060 % calculate the one_step_inverse using the standard
0061 % formulation (JtWJ + hp^2*RtR)\JtW
0062 function one_step_inv= standard_form( inv_model, J )
0063     RtR = calc_RtR_prior( inv_model );
0064     W   = calc_meas_icov( inv_model );
0065     hp  = calc_hyperparameter( inv_model );
0066 
0067     time_steps = inv_model.time_prior_solve.time_steps;
0068     l_ts  = time_steps*2 + 1;
0069 
0070     JtWJ = kron( speye(l_ts), J'*W*J);
0071     JtW  = kron( speye(l_ts), J'*W);
0072     one_step_inv= (JtWJ +  hp^2*RtR)\JtW;
0073 
0074     n_el= size(J,2);
0075     one_step_inv= one_step_inv(n_el*time_steps + (1:n_el),:);
0076 
0077 % calculate the one_step_inverse using the data form
0078 % CovX * J' * inv(J*CovX*J' + CovZ)
0079 %   iRtR*Jt/(Ji*RtR*Jt +  hp^2*iW);
0080 function one_step_inv= data_form( inv_model, J );
0081     space_prior= inv_model.time_smooth_prior.space_prior;
0082     time_weight= inv_model.time_smooth_prior.time_weight;
0083     ts         = inv_model.time_prior_solve.time_steps;
0084 
0085     space_Reg= feval(space_prior, inv_model);
0086 
0087     iRtRJt_frac=  (space_Reg\J');
0088     JiRtRJt_frac= J*iRtRJt_frac;
0089 
0090     % JiRtRJt_mult accounts for different parts of the
0091     % frame being taken at different times
0092     delta_vec= calc_delta( inv_model, J);
0093     delta_vec1= delta_vec*ones(1,length(delta_vec));
0094     JiRtRJt_mult = time_weight.^abs(delta_vec1 - delta_vec1');
0095 
0096     [x,y]= meshgrid(-ts:ts,  -ts:ts);
0097     time_w_mat= time_weight.^abs(x-y);
0098 
0099     JiRtRJt= kron( time_w_mat, JiRtRJt_frac .* JiRtRJt_mult );
0100 
0101     %FIXME: do we multiply be JiRtRJt_mult here?
0102     iRtRJt=  kron( time_w_mat(ts+1,:), iRtRJt_frac );
0103 
0104     iW   = kron( speye(1+2*ts), inv( ...
0105                  calc_meas_icov( inv_model ) ));
0106     hp   = calc_hyperparameter( inv_model );
0107 
0108     one_step_inv= iRtRJt/(JiRtRJt +  hp^2*iW);
0109 
0110 % if measurements are taken at different times,
0111 % then calculate a delta of each wrt the centre time
0112 function delta_vec= calc_delta( inv_model, J)
0113    stimulation= inv_model.fwd_model.stimulation;
0114    n_N= size(J,1);
0115 
0116    if isfield(stimulation(1),'delta_time')
0117       delta_time= [stimulation(:).delta_time];
0118       if diff(delta_time) ~= 0;
0119          error('All time steps must be same for kalman filter');
0120       end
0121    else
0122       delta_time=0;
0123    end
0124 
0125    % sequence is a vector location of each stimulation in the frame
0126    if delta_time == 0
0127       seq= size(J,1);
0128    else
0129       for i=1:length(stimulation)
0130          seq(i) = size(stimulation(i).meas_pattern,1);
0131       end
0132       seq= cumsum( seq );
0133    end
0134 
0135    delta_time= cumsum(delta_time);
0136 
0137    delta_vec= zeros(size(J,1),1);
0138    seq= [0;seq(:)];
0139    for i=1:length(seq)-1
0140       delta_vec( (seq(i)+1):seq(i+1) )= delta_time(i);
0141    end
0142 
0143    % normalize so middle time is centre, and max time is 1
0144    delta_vec= (delta_vec - mean(delta_vec)) / ...
0145               (sum(delta_time) + eps );
0146 
0147