inv_solve_time_prior

 INV_SOLVE_TIME_PRIOR inverse solver to account for time differences
 img= inv_solve_time_prior( 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

 Parameters
   inv_model.inv_solve_time_prior.time_steps => time_steps

0001 function img= inv_solve_time_prior( inv_model, data1, data2)
0002 % INV_SOLVE_TIME_PRIOR inverse solver to account for time differences
0003 % img= inv_solve_time_prior( 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 % Parameters
0015 %   inv_model.inv_solve_time_prior.time_steps => time_steps
0016 
0017 % TODO: This function really should be calling the proper
0018 %   prior calculator functions, and not reimplementing
0019 %   them internally
0020 
0021 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0022 % $Id: inv_solve_time_prior.m 5112 2015-06-14 13:00:41Z aadler $
0023 
0024 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'); do_unit_test; return; end
0025 
0026 fwd_model= inv_model.fwd_model;
0027 time_steps = inv_model.inv_solve_time_prior.time_steps;
0028 l_ts  = time_steps*2 + 1;
0029 
0030 % The one_step reconstruction matrix is cached
0031 
0032 one_step_inv= eidors_cache(@data_form, inv_model, 'inv_solve_time_prior' );
0033 
0034 dva = calc_difference_data( data1, data2, inv_model.fwd_model);
0035 
0036 l_dva = size( dva, 2);
0037 
0038 idx= [-time_steps:time_steps]'*ones(1,l_dva) + ...
0039      ones(l_ts,1)*(1:l_dva);
0040 % replicate first and last measurements
0041 idx(idx<1) = 1;
0042 idx(idx>l_dva) = l_dva;
0043 
0044 dvat= reshape(dva(:,idx),[],l_dva);
0045  
0046 sol = one_step_inv * dvat;
0047 
0048 % create a data structure to return
0049 img.name= 'solved by inv_solve_time_prior';
0050 img.elem_data = sol;
0051 img.inv_model= inv_model;
0052 img.fwd_model= fwd_model;
0053 
0054 % calculate the one_step_inverse using the standard
0055 % formulation (JtWJ + hp^2*RtR)\JtW
0056 function one_step_inv= standard_form( inv_model)
0057     RtR = calc_RtR_prior( inv_model );
0058     W   = calc_meas_icov( inv_model );
0059     hp  = calc_hyperparameter( inv_model );
0060 
0061     time_steps = inv_model.inv_solve_time_prior.time_steps;
0062     l_ts  = time_steps*2 + 1;
0063 
0064     JtWJ = kron( speye(l_ts), J'*W*J);
0065     JtW  = kron( speye(l_ts), J'*W);
0066     one_step_inv= (JtWJ +  hp^2*RtR)\JtW;
0067 
0068     n_el= size(J,2);
0069     one_step_inv= one_step_inv(n_el*time_steps + (1:n_el),:);
0070 
0071 % calculate the one_step_inverse using the data form
0072 % CovX * J' * inv(J*CovX*J' + CovZ)
0073 %   iRtR*Jt/(Ji*RtR*Jt +  hp^2*iW);
0074 function one_step_inv= data_form( inv_model)
0075     img_bkgnd= calc_jacobian_bkgnd( inv_model );
0076     J = calc_jacobian(img_bkgnd);
0077     space_prior= inv_model.prior_time_smooth.space_prior;
0078     time_weight= inv_model.prior_time_smooth.time_weight;
0079     ts         = inv_model.inv_solve_time_prior.time_steps;
0080 
0081     space_Reg= feval(space_prior, inv_model);
0082 
0083     iRtRJt_frac=  (space_Reg\J');
0084     JiRtRJt_frac= J*iRtRJt_frac;
0085 
0086     % JiRtRJt_mult accounts for different parts of the
0087     % frame being taken at different times
0088     delta_vec= calc_delta( inv_model, J);
0089     delta_vec1= delta_vec*ones(1,length(delta_vec));
0090     JiRtRJt_mult = time_weight.^abs(delta_vec1 - delta_vec1');
0091 
0092     [x,y]= meshgrid(-ts:ts,  -ts:ts);
0093     time_w_mat= time_weight.^abs(x-y);
0094 
0095     JiRtRJt= kron( time_w_mat, JiRtRJt_frac .* JiRtRJt_mult );
0096 
0097     %FIXME: do we multiply be JiRtRJt_mult here?
0098     iRtRJt=  kron( time_w_mat(ts+1,:), iRtRJt_frac );
0099 
0100     iW   = kron( speye(1+2*ts), inv( ...
0101                  calc_meas_icov( inv_model ) ));
0102     hp   = calc_hyperparameter( inv_model );
0103 
0104     one_step_inv= iRtRJt/(JiRtRJt +  hp^2*iW);
0105 
0106 % if measurements are taken at different times,
0107 % then calculate a delta of each wrt the centre time
0108 function delta_vec= calc_delta( inv_model, J)
0109    stimulation= inv_model.fwd_model.stimulation;
0110    n_N= size(J,1);
0111 
0112    if isfield(stimulation(1),'delta_time')
0113       delta_time= [stimulation(:).delta_time];
0114       if diff(delta_time) ~= 0;
0115          error('All time steps must be same for kalman filter');
0116       end
0117    else
0118       delta_time=0;
0119    end
0120 
0121    % sequence is a vector location of each stimulation in the frame
0122    if delta_time == 0
0123       seq= size(J,1);
0124    else
0125       for i=1:length(stimulation)
0126          seq(i) = size(stimulation(i).meas_pattern,1);
0127       end
0128       seq= cumsum( seq );
0129    end
0130 
0131    delta_time= cumsum(delta_time);
0132 
0133    delta_vec= zeros(size(J,1),1);
0134    seq= [0;seq(:)];
0135    for i=1:length(seq)-1
0136       delta_vec( (seq(i)+1):seq(i+1) )= delta_time(i);
0137    end
0138 
0139    % normalize so middle time is centre, and max time is 1
0140    delta_vec= (delta_vec - mean(delta_vec)) / ...
0141               (sum(delta_time) + eps );
0142 
0143    
0144 function do_unit_test
0145    do_unit_test_2d
0146    do_unit_test_3d
0147 
0148 function do_unit_test_2d
0149    time_steps=  3;
0150    time_weight= .8;
0151 
0152    [vh,vi,xyr_pt]=simulate_2d_movement(7, [], [],[1,0.5,0.4]);
0153 
0154    imdl_TS = mk_common_model( 'c2c2', 16 ); % 576 element
0155    imdl_TS.fwd_model.normalize_measurements= 0;
0156    imdl_TS.hyperparameter.value= 0.10;
0157 
0158    imdl_TS.RtR_prior= @prior_time_smooth;
0159    imdl_TS.prior_time_smooth.space_prior= @prior_noser;
0160    imdl_TS.prior_noser.exponent= .5;
0161    imdl_TS.prior_time_smooth.time_weight= time_weight;
0162    imdl_TS.prior_time_smooth.time_steps=  time_steps;
0163    imdl_TS.solve= @inv_solve_time_prior;
0164    imdl_TS.inv_solve_time_prior.time_steps=   time_steps;
0165 
0166 image_select= floor(length(xyr_pt)/2)+1;  % this image is at 9 O'Clock
0167 time_steps=  3; ts_expand= 1;
0168 time_weight= .8;
0169 ts_vec= -time_steps:time_steps;
0170 
0171     im_sel= image_select+ ts_vec*ts_expand;
0172     vi_sel= vi(:,im_sel);
0173     sel  = 1 + time_steps; % choose the middle
0174 
0175    img= inv_solve( imdl_TS, vh, vi_sel);
0176    img.elem_data= img.elem_data(:,sel);
0177    show_fem(img);
0178 
0179 
0180 function do_unit_test_3d
0181    n_sims= 20;
0182    stim = mk_stim_patterns(16,2,'{ad}','{ad}',{},1);
0183    fmdl = mk_library_model('cylinder_16x2el_vfine');
0184    fmdl.stimulation = stim;
0185    [vh,vi,xyzr_pt]= simulate_3d_movement( n_sims, fmdl);
0186