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 4833 2015-03-29 21:32:08Z bgrychtol-ipa $
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    RM= (J'*W*J +  hp^2*RtR)\J'*W;
0064 
0065    
0066    
0067 function [img, step_size] = scale_to_fit_data(img, inv_model, data1, data2)
0068    % find the step size to multiply sol by to best fit data
0069    step_size = 1;
0070    do_step   = false;
0071    % If calc_step_size, ignore specified step_size
0072    try do_step = inv_model.inv_solve_diff_GN_one_step.calc_step_size; end
0073    
0074    if do_step
0075       eidors_msg('inv_solve_diff_GN_one_step: Calculating optimal step size to fit data',2);
0076       % options for fminbnd
0077       try 
0078          opt = inv_model.inv_solve_diff_GN_one_step.fminbnd;
0079       catch
0080          opt.Display = 'iter';
0081       end
0082       % range for fminbnd
0083       try
0084          range = inv_model.inv_solve_diff_GN_one_step.bounds;
0085       catch
0086          range = [1e-5 1e1];
0087       end
0088       step_size = fminbnd(@(x) to_optimize(img,inv_model,data1,data2, x), ...
0089                            range(1), range(2), opt);
0090    else
0091       % if not calculating, check if step_size provided
0092       try
0093          step_size = inv_model.inv_solve_diff_GN_one_step.step_size;
0094       end
0095    end
0096    img.elem_data = img.elem_data * step_size;
0097    img.info.step_size = step_size;
0098 
0099 function out = to_optimize(img, inv_model, data1, data2, x)
0100    img.elem_data = img.elem_data*x;
0101    out = calc_solution_error(img, inv_model, data1, data2);
0102

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005