INV_SOLVE_DIFF_GN_ONE_STEP inverse solver using approach of Adler&Guardo 1996
img= inv_solve_diff_GN_one_step( 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
By default, the correct scaling of the solution that best fits the data
is not calculated, possibly resulting in high solution errors reported by
inv_solve.
To calculate the correct scaling, specify
inv_model.inv_solve_diff_GN_one_step.calc_step_size = 1;
To provide a pre-calculated scaling, specify
inv_model.inv_solve_diff_GN_one_step.calc_step_size = 0;
inv_model.inv_solve_diff_GN_one_step.step_size = 0.8;
The search for correct step_size is performed using FMINBND. The default
search interval is [1e-5 1e1]. You can modify it by specifying:
inv_model.inv_solve_diff_GN_one_step.bounds = [10 200];
Additional options for FMINBD can be passed as:
inv_model.inv_solve_diff_GN_one_step.fminbnd.MaxIter = 10;
The optimal step_size is returned in img.info.step_size.
See also INV_SOLVE, CALC_SOLUTION_ERROR, FMINBND0001 function img= inv_solve_diff_GN_one_step( inv_model, data1, data2) 0002 % INV_SOLVE_DIFF_GN_ONE_STEP inverse solver using approach of Adler&Guardo 1996 0003 % img= inv_solve_diff_GN_one_step( 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 % By default, the correct scaling of the solution that best fits the data 0015 % is not calculated, possibly resulting in high solution errors reported by 0016 % inv_solve. 0017 % To calculate the correct scaling, specify 0018 % inv_model.inv_solve_diff_GN_one_step.calc_step_size = 1; 0019 % To provide a pre-calculated scaling, specify 0020 % inv_model.inv_solve_diff_GN_one_step.calc_step_size = 0; 0021 % inv_model.inv_solve_diff_GN_one_step.step_size = 0.8; 0022 % The search for correct step_size is performed using FMINBND. The default 0023 % search interval is [1e-5 1e1]. You can modify it by specifying: 0024 % inv_model.inv_solve_diff_GN_one_step.bounds = [10 200]; 0025 % Additional options for FMINBD can be passed as: 0026 % inv_model.inv_solve_diff_GN_one_step.fminbnd.MaxIter = 10; 0027 % 0028 % The optimal step_size is returned in img.info.step_size. 0029 % 0030 % See also INV_SOLVE, CALC_SOLUTION_ERROR, FMINBND 0031 0032 % (C) 2005-2013 Andy Adler and Bartlomiej Grychtol. 0033 % License: GPL version 2 or version 3 0034 % $Id: inv_solve_diff_GN_one_step.m 5611 2017-06-27 12:02:50Z htregidgo $ 0035 0036 % TODO: 0037 % Test whether Wiener filter form or Tikhonov form are faster 0038 % Tikhonov: RM= (J'*W*J + hp^2*RtR)\J'*W; 0039 % Wiener: P= inv(RtR); V = inv(W); RM = P*J'/(J*P*J' + hp^2*V) 0040 0041 dv = calc_difference_data( data1, data2, inv_model.fwd_model); 0042 sol = eidors_cache(@get_RM, inv_model,'inv_solve_diff_GN_one_step' ) * dv; 0043 0044 0045 0046 img = data_mapper(calc_jacobian_bkgnd( inv_model )); 0047 img.name= 'solved by inv_solve_diff_GN_one_step'; 0048 img.elem_data = sol; 0049 img.fwd_model= inv_model.fwd_model; 0050 img = data_mapper(img,1); 0051 0052 img = scale_to_fit_data(img, inv_model, data1, data2); 0053 0054 0055 function RM = get_RM( inv_model ) 0056 img_bkgnd= calc_jacobian_bkgnd( inv_model ); 0057 J = calc_jacobian( img_bkgnd); 0058 0059 RtR = calc_RtR_prior( inv_model ); 0060 W = calc_meas_icov( inv_model ); 0061 hp = calc_hyperparameter( inv_model ); 0062 0063 % left_divide now has handling to force symmetric methods when matrices 0064 % are symmetric up to floating point error. 0065 RM = left_divide((J'*W*J + hp^2*RtR),J'*W); 0066 0067 0068 0069 0070 0071 function [img, step_size] = scale_to_fit_data(img, inv_model, data1, data2) 0072 % find the step size to multiply sol by to best fit data 0073 step_size = 1; 0074 do_step = false; 0075 % If calc_step_size, ignore specified step_size 0076 try do_step = inv_model.inv_solve_diff_GN_one_step.calc_step_size; end 0077 0078 if do_step 0079 eidors_msg('inv_solve_diff_GN_one_step: Calculating optimal step size to fit data',2); 0080 % options for fminbnd 0081 try 0082 opt = inv_model.inv_solve_diff_GN_one_step.fminbnd; 0083 catch 0084 opt.Display = 'iter'; 0085 end 0086 % range for fminbnd 0087 try 0088 range = inv_model.inv_solve_diff_GN_one_step.bounds; 0089 catch 0090 range = [1e-5 1e1]; 0091 end 0092 step_size = fminbnd(@(x) to_optimize(img,inv_model,data1,data2, x), ... 0093 range(1), range(2), opt); 0094 else 0095 % if not calculating, check if step_size provided 0096 try 0097 step_size = inv_model.inv_solve_diff_GN_one_step.step_size; 0098 end 0099 end 0100 img.elem_data = img.elem_data * step_size; 0101 img.info.step_size = step_size; 0102 0103 function out = to_optimize(img, inv_model, data1, data2, x) 0104 img.elem_data = img.elem_data*x; 0105 out = calc_solution_error(img, inv_model, data1, data2); 0106