0001 function img= time_prior_solve( inv_model, data1, data2)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
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
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
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
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
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
0061
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
0078
0079
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
0091
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
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
0111
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
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
0144 delta_vec= (delta_vec - mean(delta_vec)) / ...
0145 (sum(delta_time) + eps );
0146
0147