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