inv_solve_diff_GN_one_step

PURPOSE ^

INV_SOLVE_DIFF_GN_ONE_STEP inverse solver using approach of Adler&Guardo 1996

SYNOPSIS ^

function img= inv_solve_diff_GN_one_step( inv_model, data1, data2)

DESCRIPTION ^

 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, FMINBND

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 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

Generated on Fri 30-Dec-2022 19:44:54 by m2html © 2005