inv_solve_core

PURPOSE ^

INV_SOLVE_CORE Solver using a generic iterative algorithm

SYNOPSIS ^

function img= inv_solve_core( inv_model, data0, data1);

DESCRIPTION ^

INV_SOLVE_CORE Solver using a generic iterative algorithm
 img        => output image (or vector of images)
 inv_model  => inverse model struct
 data0      => EIT data
 data0, data1 => difference EIT data

 This function is parametrized and uses function pointers where possible to
 allow its use as a general iterative solver framework. There are a large
 number of parameters and functions contained here. Sensible defaults are
 used throughout. You do not need to set every parameter.

 The solver operates as an absolute Gauss-Newton iterative solver by default.
 Wrapper functions are available to call this function in its various forms.
 Look forward to the "See also" section at the end of this help.

 Argument matrices to the internal functions (measurement inverse covariance,
 for example) are only calculated if required. Functions that are supplied to
 this INV_SOLVE_CORE must be able to survive being probed: they will have
 each parameter set to either 0 or 1 to determine if the function is sensitive
 to that argument. This works cleanly for most matrix multiplication based
 functions but for more abstract code, some handling of this behaviour may
 need to be implemented.

 In the following parameters, r_k is the current residual, r_{k-1} is the
 previous iteration's residual. k is the iteration count.

 Parameters denoted with a ** to the right of their default values are
 deprecated legacy parameters, some of which formerly existed under
 'inv_model.parameters.*'.

 Parameters (inv_model.inv_solve_core.*):
   verbose (show progress)                (default 2)
      0: quiet
    >=1: print iteration count
    >=2: print details as the algorithm progresses
    >=3: plot residuals versus iteration count
    >=4: plot result at each iteration, see show_fem
    >=5: plot line search per iteration
   plot_residuals                         (default 0)
    plot residuals without verbose output
   fig_prefix                       (default: <none>)
    figure file prefix; figures not saved if <none>
   fwd_solutions                          (default 0)
    0: ignore
    1: count fwd_solve(), generally the most
       computationally expensive component of
       the iterations
   residual_func =             (default @GN_residual)
    NOTE: @meas_residual exists to maintain
    compatibility with some older code
   max_iterations                        (default 10)  **
   ntol (estimate of machine precision) (default eps)
   tol (stop iter if r_k < tol)           (default 0)
   dtol                              (default -0.01%)
    stop iter if (r_k - r_{k-1})/r_1 < dtol AND
                 k >= dtol_iter
   dtol_iter                              (default 0)
    apply dtol stopping criteria if k >= dtol_iter
   min_value                           (default -inf)  **
   max_value                           (default +inf)  **
   line_optimize_func                (default <none>)  ** TODO
     [next,fmin,res]=f(org,dx,data0,opt);
     opt=line_optimize.* + objective_func
     Deprecated, use line_search_func instead.
   line_optimize.perturb
                   (default line_search_args.perturb)  ** TODO
     Deprecated, use line_search_args.perturb instead.
   update_func                         (default TODO)  ** TODO
     [img,opt]=f(org,next,dx,fmin,res,opt)
     Deprecated, use <TODO> instead.
   update_method                         (default lu)
     Method to use for solving dx: 'pcg' or 'lu'
     If 'lu', will fall back to 'pcg' if out-of-memory.
   do_starting_estimate                   (default 1)  ** TODO
     Deprecated, use <TODO> instead.
   line_search_func       (default @line_search_onm2)
   line_search_dv_func      (default @update_dv_core)
   line_search_de_func      (default @update_de_core)
   line_search_args.perturb
                     (default [0 1/16 1/8 1/4 1/2 1])
    line search for alpha by these steps along sx
   line_search_args.plot                  (default 0)
   c2f_background                         (default 0)
    if > 0, this is additional elem_data
    if a c2f map exists, the default is to decide
    based on an estimate of c2f overlap whether a
    background value is required. If a background is
    required, it is added as the last element of that
    type.
   c2f_background_fixed                   (default 1)
    hold the background estimate fixed or allow it
    to vary as any other elem_data
   elem_fixed                            (default [])
    meas_select already handles selecting from the
    valid measurements. we want the same for the
    elem_data, so we only work on modifying the
    legal values.
    Note that c2f_background's elements are added to
    this list if c2f_background_fixed == 1.
   prior_data             (default to jacobian_bkgnd)
    Sets the priors of type elem_prior. May be
    scalar, per elem_prior, or match the working
    length of each elem_data type. Note that for priors
    using the c2f a background element may be added
    to the end of that range when required; see
    c2f_background.
   elem_len                (default to all elem_data)
    A cell array list of how many of each
    elem_working there are in elem_data.
      prior_data = { 32.1, 10*ones(10,1) };
      elem_prior = {'conductivity', 'movement'};
      elem_len = { 20001, 10 };
   elem_prior                (default 'conductivity')
    Input 'prior_data' type; immediately converted to
    'elem_working' type before first iteration.
   elem_working              (default 'conductivity')
   elem_output               (default 'conductivity')
    The working and output units for 'elem_data'.
    Valid types are 'conductivity' and 'resistivity'
    as plain units or with the prefix 'log_' or
    'log10_'. Conversions are handled internally.
    Scaling factors are applied to the Jacobian
    (calculated in units of 'conductivity') as
    appropriate.
    If elem_working == elem_output, then no
    conversions take place.
    For multiple types, use cell array.
    ex: elem_output = {'log_resistivity', 'movement'}
   meas_input                     (default 'voltage')
   meas_working                   (default 'voltage')
    Similarly to elem_working/output, conversion
    between 'voltage' and 'apparent_resistivity' and
    their log/log10 variants are handled internally.
    If meas_input == meas_working no conversions take
    place. The normalization factor 'N' is calculated
    if 'apparent_resistivity' is used.
   update_img_func             (default: pass-through)
    Called prior to calc_jacobian and update_dv.
    Elements are converted to their "base types"
    before this function is called. For example,
    'log_resistivity' becomes 'conductivity'.
    It is a hook to allow additional updates to the
    model before the Jacobian, or a new set of
    measurements are calculated via fwd_solve.
   return_working_variables               (default: 0)
    If 1, return the last working variables to the user
     img.inv_solve_{type}.J   Jacobian
     img.inv_solve_{type}.dx  descent direction
     img.inv_solve_{type}.sx  search direction
     img.inv_solve_{type}.alpha  line search result
     img.inv_solve_{type}.beta   conjugation parameter
     img.inv_solve_{type}.r   as:
       [ residual, measurement misfit, element misfit ]
       with one row per iteration
   show_fem                       (default: @show_fem)
    Function with which to plot each iteration's
    current parameters.

   Signature for residual_func
    [r,m,e] = f(dv, de, W, hp2RtR)
   where
    r   - the residual = m + e
    m   - measurement error
    e   - prior misfit
    dv  - change in voltage
    de  - change in image elements
    W   - measurement inverse covariance matrix
    hp2 - hyperparameter squared, see CALC_HYPERPARAMETER
    RtR - regularization matrix squared --> hp2RtR = hp2*RtR

   Signature for line_optimize_func
    [alpha, img, dv, opt] = f(img, sx, data0, img0, N, W, hp2RtR, dv, opt)
   where:
    alpha - line search result
    img   - the current image
            (optional, recalculated if not available)
    sx    - the search direction to which alpha should be applied
    data0 - the true measurements     (dv = N*data - N*data0)
    img0  - the image background (de = img - img0)
    N     - a measurement normalization factor, N*dv
    W     - measurement inverse covariance matrix
    hp2   - hyperparameter squared, see CALC_HYPERPARAMETER
    RtR   - regularization matrix squared --> hp2RtR = hp2*RtR
    dv    - change in voltage
            (optional, recalculated if not available)
    opt   - additional arguments, updated at each call

   Signature for line_search_dv_func
    [dv, opt] = update_dv_core(img, data0, N, opt)
   where:
    dv    - change in voltage
    opt   - additional arguments, updated at each call
    data  - the estimated measurements
    img   - the current image
    data0 - the true measurements
    N     - a measurement normalization factor, N*dv

   Signature for line_search_de_func
    de = f(img, img0, opt)
   where:
    de    - change in image elements
    img   - the current image
    img0  - the image background (de = img - img0)
    opt   - additional arguments

   Signature for calc_jacobian_scaling_func
    S = f(x)
   where:
    S - to be used to scale the Jacobian
    x - current img.elem_data in units of 'conductivity'

   Signature for  update_img_func
    img2 = f(img1, opt)
   where
    img1 - an input image, the current working image
    img2 - a (potentially) modified version to be used
    for the fwd_solve/Jacobian calculations

 NOTE that the default line search is very crude. For
 my test problems it seems to amount to an expensive grid
 search. Much more efficient line search algorithms exist
 and some fragments already are coded elsewhere in the
 EIDORS code-base.

 See also: INV_SOLVE_ABS_GN, INV_SOLVE_ABS_GN_LOGC,
           INV_SOLVE_ABS_CG, INV_SOLVE_ABS_CG_LOGC,
           LINE_SEARCH_O2, LINE_SEARCH_ONM2

 (C) 2010-2016 Alistair Boyle, Nolwenn Lesparre, Andy Adler, Bartłomiej Grychtol.
 License: GPL version 2 or version 3

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function img= inv_solve_core( inv_model, data0, data1);
0002 %INV_SOLVE_CORE Solver using a generic iterative algorithm
0003 % img        => output image (or vector of images)
0004 % inv_model  => inverse model struct
0005 % data0      => EIT data
0006 % data0, data1 => difference EIT data
0007 %
0008 % This function is parametrized and uses function pointers where possible to
0009 % allow its use as a general iterative solver framework. There are a large
0010 % number of parameters and functions contained here. Sensible defaults are
0011 % used throughout. You do not need to set every parameter.
0012 %
0013 % The solver operates as an absolute Gauss-Newton iterative solver by default.
0014 % Wrapper functions are available to call this function in its various forms.
0015 % Look forward to the "See also" section at the end of this help.
0016 %
0017 % Argument matrices to the internal functions (measurement inverse covariance,
0018 % for example) are only calculated if required. Functions that are supplied to
0019 % this INV_SOLVE_CORE must be able to survive being probed: they will have
0020 % each parameter set to either 0 or 1 to determine if the function is sensitive
0021 % to that argument. This works cleanly for most matrix multiplication based
0022 % functions but for more abstract code, some handling of this behaviour may
0023 % need to be implemented.
0024 %
0025 % In the following parameters, r_k is the current residual, r_{k-1} is the
0026 % previous iteration's residual. k is the iteration count.
0027 %
0028 % Parameters denoted with a ** to the right of their default values are
0029 % deprecated legacy parameters, some of which formerly existed under
0030 % 'inv_model.parameters.*'.
0031 %
0032 % Parameters (inv_model.inv_solve_core.*):
0033 %   verbose (show progress)                (default 2)
0034 %      0: quiet
0035 %    >=1: print iteration count
0036 %    >=2: print details as the algorithm progresses
0037 %    >=3: plot residuals versus iteration count
0038 %    >=4: plot result at each iteration, see show_fem
0039 %    >=5: plot line search per iteration
0040 %   plot_residuals                         (default 0)
0041 %    plot residuals without verbose output
0042 %   fig_prefix                       (default: <none>)
0043 %    figure file prefix; figures not saved if <none>
0044 %   fwd_solutions                          (default 0)
0045 %    0: ignore
0046 %    1: count fwd_solve(), generally the most
0047 %       computationally expensive component of
0048 %       the iterations
0049 %   residual_func =             (default @GN_residual)
0050 %    NOTE: @meas_residual exists to maintain
0051 %    compatibility with some older code
0052 %   max_iterations                        (default 10)  **
0053 %   ntol (estimate of machine precision) (default eps)
0054 %   tol (stop iter if r_k < tol)           (default 0)
0055 %   dtol                              (default -0.01%)
0056 %    stop iter if (r_k - r_{k-1})/r_1 < dtol AND
0057 %                 k >= dtol_iter
0058 %   dtol_iter                              (default 0)
0059 %    apply dtol stopping criteria if k >= dtol_iter
0060 %   min_value                           (default -inf)  **
0061 %   max_value                           (default +inf)  **
0062 %   line_optimize_func                (default <none>)  ** TODO
0063 %     [next,fmin,res]=f(org,dx,data0,opt);
0064 %     opt=line_optimize.* + objective_func
0065 %     Deprecated, use line_search_func instead.
0066 %   line_optimize.perturb
0067 %                   (default line_search_args.perturb)  ** TODO
0068 %     Deprecated, use line_search_args.perturb instead.
0069 %   update_func                         (default TODO)  ** TODO
0070 %     [img,opt]=f(org,next,dx,fmin,res,opt)
0071 %     Deprecated, use <TODO> instead.
0072 %   update_method                         (default lu)
0073 %     Method to use for solving dx: 'pcg' or 'lu'
0074 %     If 'lu', will fall back to 'pcg' if out-of-memory.
0075 %   do_starting_estimate                   (default 1)  ** TODO
0076 %     Deprecated, use <TODO> instead.
0077 %   line_search_func       (default @line_search_onm2)
0078 %   line_search_dv_func      (default @update_dv_core)
0079 %   line_search_de_func      (default @update_de_core)
0080 %   line_search_args.perturb
0081 %                     (default [0 1/16 1/8 1/4 1/2 1])
0082 %    line search for alpha by these steps along sx
0083 %   line_search_args.plot                  (default 0)
0084 %   c2f_background                         (default 0)
0085 %    if > 0, this is additional elem_data
0086 %    if a c2f map exists, the default is to decide
0087 %    based on an estimate of c2f overlap whether a
0088 %    background value is required. If a background is
0089 %    required, it is added as the last element of that
0090 %    type.
0091 %   c2f_background_fixed                   (default 1)
0092 %    hold the background estimate fixed or allow it
0093 %    to vary as any other elem_data
0094 %   elem_fixed                            (default [])
0095 %    meas_select already handles selecting from the
0096 %    valid measurements. we want the same for the
0097 %    elem_data, so we only work on modifying the
0098 %    legal values.
0099 %    Note that c2f_background's elements are added to
0100 %    this list if c2f_background_fixed == 1.
0101 %   prior_data             (default to jacobian_bkgnd)
0102 %    Sets the priors of type elem_prior. May be
0103 %    scalar, per elem_prior, or match the working
0104 %    length of each elem_data type. Note that for priors
0105 %    using the c2f a background element may be added
0106 %    to the end of that range when required; see
0107 %    c2f_background.
0108 %   elem_len                (default to all elem_data)
0109 %    A cell array list of how many of each
0110 %    elem_working there are in elem_data.
0111 %      prior_data = { 32.1, 10*ones(10,1) };
0112 %      elem_prior = {'conductivity', 'movement'};
0113 %      elem_len = { 20001, 10 };
0114 %   elem_prior                (default 'conductivity')
0115 %    Input 'prior_data' type; immediately converted to
0116 %    'elem_working' type before first iteration.
0117 %   elem_working              (default 'conductivity')
0118 %   elem_output               (default 'conductivity')
0119 %    The working and output units for 'elem_data'.
0120 %    Valid types are 'conductivity' and 'resistivity'
0121 %    as plain units or with the prefix 'log_' or
0122 %    'log10_'. Conversions are handled internally.
0123 %    Scaling factors are applied to the Jacobian
0124 %    (calculated in units of 'conductivity') as
0125 %    appropriate.
0126 %    If elem_working == elem_output, then no
0127 %    conversions take place.
0128 %    For multiple types, use cell array.
0129 %    ex: elem_output = {'log_resistivity', 'movement'}
0130 %   meas_input                     (default 'voltage')
0131 %   meas_working                   (default 'voltage')
0132 %    Similarly to elem_working/output, conversion
0133 %    between 'voltage' and 'apparent_resistivity' and
0134 %    their log/log10 variants are handled internally.
0135 %    If meas_input == meas_working no conversions take
0136 %    place. The normalization factor 'N' is calculated
0137 %    if 'apparent_resistivity' is used.
0138 %   update_img_func             (default: pass-through)
0139 %    Called prior to calc_jacobian and update_dv.
0140 %    Elements are converted to their "base types"
0141 %    before this function is called. For example,
0142 %    'log_resistivity' becomes 'conductivity'.
0143 %    It is a hook to allow additional updates to the
0144 %    model before the Jacobian, or a new set of
0145 %    measurements are calculated via fwd_solve.
0146 %   return_working_variables               (default: 0)
0147 %    If 1, return the last working variables to the user
0148 %     img.inv_solve_{type}.J   Jacobian
0149 %     img.inv_solve_{type}.dx  descent direction
0150 %     img.inv_solve_{type}.sx  search direction
0151 %     img.inv_solve_{type}.alpha  line search result
0152 %     img.inv_solve_{type}.beta   conjugation parameter
0153 %     img.inv_solve_{type}.r   as:
0154 %       [ residual, measurement misfit, element misfit ]
0155 %       with one row per iteration
0156 %   show_fem                       (default: @show_fem)
0157 %    Function with which to plot each iteration's
0158 %    current parameters.
0159 %
0160 %   Signature for residual_func
0161 %    [r,m,e] = f(dv, de, W, hp2RtR)
0162 %   where
0163 %    r   - the residual = m + e
0164 %    m   - measurement error
0165 %    e   - prior misfit
0166 %    dv  - change in voltage
0167 %    de  - change in image elements
0168 %    W   - measurement inverse covariance matrix
0169 %    hp2 - hyperparameter squared, see CALC_HYPERPARAMETER
0170 %    RtR - regularization matrix squared --> hp2RtR = hp2*RtR
0171 %
0172 %   Signature for line_optimize_func
0173 %    [alpha, img, dv, opt] = f(img, sx, data0, img0, N, W, hp2RtR, dv, opt)
0174 %   where:
0175 %    alpha - line search result
0176 %    img   - the current image
0177 %            (optional, recalculated if not available)
0178 %    sx    - the search direction to which alpha should be applied
0179 %    data0 - the true measurements     (dv = N*data - N*data0)
0180 %    img0  - the image background (de = img - img0)
0181 %    N     - a measurement normalization factor, N*dv
0182 %    W     - measurement inverse covariance matrix
0183 %    hp2   - hyperparameter squared, see CALC_HYPERPARAMETER
0184 %    RtR   - regularization matrix squared --> hp2RtR = hp2*RtR
0185 %    dv    - change in voltage
0186 %            (optional, recalculated if not available)
0187 %    opt   - additional arguments, updated at each call
0188 %
0189 %   Signature for line_search_dv_func
0190 %    [dv, opt] = update_dv_core(img, data0, N, opt)
0191 %   where:
0192 %    dv    - change in voltage
0193 %    opt   - additional arguments, updated at each call
0194 %    data  - the estimated measurements
0195 %    img   - the current image
0196 %    data0 - the true measurements
0197 %    N     - a measurement normalization factor, N*dv
0198 %
0199 %   Signature for line_search_de_func
0200 %    de = f(img, img0, opt)
0201 %   where:
0202 %    de    - change in image elements
0203 %    img   - the current image
0204 %    img0  - the image background (de = img - img0)
0205 %    opt   - additional arguments
0206 %
0207 %   Signature for calc_jacobian_scaling_func
0208 %    S = f(x)
0209 %   where:
0210 %    S - to be used to scale the Jacobian
0211 %    x - current img.elem_data in units of 'conductivity'
0212 %
0213 %   Signature for  update_img_func
0214 %    img2 = f(img1, opt)
0215 %   where
0216 %    img1 - an input image, the current working image
0217 %    img2 - a (potentially) modified version to be used
0218 %    for the fwd_solve/Jacobian calculations
0219 %
0220 % NOTE that the default line search is very crude. For
0221 % my test problems it seems to amount to an expensive grid
0222 % search. Much more efficient line search algorithms exist
0223 % and some fragments already are coded elsewhere in the
0224 % EIDORS code-base.
0225 %
0226 % See also: INV_SOLVE_ABS_GN, INV_SOLVE_ABS_GN_LOGC,
0227 %           INV_SOLVE_ABS_CG, INV_SOLVE_ABS_CG_LOGC,
0228 %           LINE_SEARCH_O2, LINE_SEARCH_ONM2
0229 %
0230 % (C) 2010-2016 Alistair Boyle, Nolwenn Lesparre, Andy Adler, Bartłomiej Grychtol.
0231 % License: GPL version 2 or version 3
0232 
0233 % $Id: inv_solve_core.m 5563 2017-06-19 02:33:27Z alistair_boyle $
0234 
0235 %--------------------------
0236 % UNIT_TEST?
0237 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST') && (nargin == 1); img=do_unit_test; return; end
0238 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST') && (nargin == 2); img=do_unit_test(data0); return; end
0239 
0240 %--------------------------
0241 opt = parse_options(inv_model);
0242 if opt.verbose > 1
0243    fprintf('  verbose = %d\n', opt.verbose);
0244 end
0245 %if opt.do_starting_estimate
0246 %    img = initial_estimate( inv_model, data0 ); % TODO
0247 %%%    AB->NL this is Nolwenn's homogeneous estimate...
0248 %%%    calc_background_resistivity is my version of this code
0249 %%%    that is working for my data set
0250 %else
0251 [inv_model, opt] = append_c2f_background(inv_model, opt);
0252 % calc_jacobian_bkgnd, used by mk_image does not understand
0253 % the course-to-fine mapping and explodes when it is fed a
0254 % prior based on the coarse model. Here we give that
0255 % function something it can swallow, then create then plug
0256 % in the correct prior afterwards.
0257 if isfield(inv_model, 'jacobian_bkgnd')
0258   inv_model = rmfield(inv_model,'jacobian_bkgnd');
0259 end
0260 inv_model.jacobian_bkgnd.value = 1;
0261 img = mk_image( inv_model );
0262 img.inv_model = inv_model; % stash the inverse model
0263 img = data_mapper(img); % move data from whatever 'params' to img.elem_data
0264 
0265 % insert the prior data
0266 img = init_elem_data(img, opt);
0267 
0268 % map data and measurements to working types
0269 %  convert elem_data
0270 img = map_img(img, opt.elem_working);
0271 
0272 % solve for difference data?
0273 if nargin == 3
0274    assert(all(~strcmp(opt.meas_working, {'apparent_resistivity','log_apparent_resistivity', 'log10_apparent_resitivity'})), ...
0275           ['meas_working = ''' opt.meas_working ''' not yet supported for difference solutions']);
0276    for i = 1:length(opt.elem_output)
0277       assert(any(strcmp(opt.elem_output{i}, {'conductivity','resistivity','movement'})), ...
0278              ['elem_output = {' strjoin(opt.elem_output,', ') '} but difference solver log normal outputs are not supported']);
0279    end
0280    % dv = (meas1 - meas0) + meas@backgnd
0281    nil = struct;
0282    if isstruct(data0)
0283       nil.meas = data0.meas*0;
0284    else
0285       nil.meas = data0*0;
0286    end
0287    nil.type = 'data';
0288    nil.current_params = opt.meas_input;
0289    [dv, opt, err] = update_dv_core(img, nil, 1, opt);
0290    data0 = calc_difference_data( data0, data1, inv_model.fwd_model) + dv;
0291    assert(strcmp(inv_model.reconst_type, 'difference'), ...
0292           ['expected inv_model.reconst_type = ''difference'' not ' inv_model.reconst_type]);
0293 else
0294    assert(strcmp(inv_model.reconst_type, 'absolute'), ...
0295           ['expected inv_model.reconst_type = ''absolute'' not ' inv_model.reconst_type]);
0296 end
0297 
0298 %  convert measurement data
0299 if ~isstruct(data0)
0300    d = data0;
0301    data0 = struct;
0302    data0.meas = d;
0303    data0.type = 'data';
0304 end
0305 data0.current_params = opt.meas_input;
0306 
0307 % precalculate some of our matrices if required
0308 W  = init_meas_icov(inv_model, opt);
0309 [N, dN] = init_normalization(inv_model.fwd_model, data0, opt);
0310 
0311 % now get on with
0312 img0 = img;
0313 hp2RtR = 0; alpha = 0; k = 0; sx = 0; r = 0; stop = 0; % general init
0314 residuals = zeros(opt.max_iterations,3); % for residuals plots
0315 dxp = 0; % previous step's slope was... nothing
0316 [dv, opt] = update_dv([], img, data0, N, opt);
0317 de = update_de([], img, img0, opt);
0318 if opt.verbose >= 5 % we only save the measurements at each iteration if we are being verbose
0319   dvall = ones(size(data0.meas,1),opt.max_iterations+1)*NaN;
0320 end
0321 while 1
0322   if opt.verbose >= 1
0323      if k == 0
0324         fprintf('  inv_solve_core: start up\n');
0325      else
0326         fprintf('  inv_solve_core: iteration %d (residual = %g)\n', k, r(k,1));
0327      end
0328   end
0329 
0330   % calculate the Jacobian
0331   %  - Jacobian before RtR because it is needed for Noser prior
0332   [J, opt] = update_jacobian(img, dN, k, opt);
0333 
0334   % update RtR, if required (depends on prior)
0335   hp2RtR = update_hp2RtR(inv_model, J, k, img, opt);
0336 
0337   % determine the next search direction sx
0338   %  dx is specific to the algorithm, generally "downhill"
0339   dx = update_dx(J, W, hp2RtR, dv, de, opt);
0340   % choose beta, beta=0 unless doing Conjugate Gradient
0341   beta = update_beta(dx, dxp, sx, opt);
0342   beta_all(k+1)=beta; % save for debug
0343   % sx_k = dx_k + beta * sx_{k-1}
0344   sx = update_sx(dx, beta, sx, opt);
0345   if k ~= 0
0346      dxp = dx; % saved for next iteration if using beta
0347   end
0348 
0349   % plot dx and SVD of Jacobian (before and after regularization)
0350   plot_dx_and_svd_elem(J, W, hp2RtR, k, sx, dx, img, opt);
0351 
0352   % line search for alpha, leaving the final selection as img
0353   % x_n = img.elem_data
0354   % x_{n+1} = x_n + \alpha sx
0355   % img.elem_data = x_{n+1}
0356   [alpha, img, dv, opt] = update_alpha(img, sx, data0, img0, N, W, hp2RtR, k, dv, opt);
0357   alpha_all(k+1) = alpha;
0358   % fix max/min values for x, clears dx if limits are hit, where
0359   % a cleared dv will trigger a recalculation of dv at the next update_dv()
0360   [img, dv] = update_img_using_limits(img, img0, data0, N, dv, opt);
0361 
0362   % update change in element data from the prior de and
0363   % the measurement error dv
0364   [dv, opt] = update_dv(dv, img, data0, N, opt);
0365   de = update_de(de, img, img0, opt);
0366   if opt.verbose >= 5
0367     dvall(:,k+1) = dv;
0368     show_meas_err(dvall, data0, k, N, W, opt);
0369   end
0370   show_fem_iter(k, img, inv_model, stop, opt);
0371 
0372   % now find the residual, quit if we're done
0373   [stop, k, r, img] = update_residual(dv, img, de, W, hp2RtR, k, r, alpha, sx, opt);
0374   if stop
0375      if stop == -1
0376         alpha_all(k) = 0;
0377      end
0378      break;
0379   end
0380 end
0381 [img, opt] = strip_c2f_background(img, opt);
0382 % check we're returning the right size of data
0383 if isfield(inv_model, 'rec_model')
0384   img.fwd_model = inv_model.rec_model;
0385 end
0386 if opt.verbose >= 1
0387    if k==1; itrs=''; else itrs='s'; end
0388    fprintf('  %d fwd_solves required for this solution in %d iteration%s\n', ...
0389            opt.fwd_solutions, k, itrs);
0390 end
0391 % convert data for output
0392 img = map_img(img, opt.elem_output);
0393 if strcmp(inv_model.reconst_type, 'difference')
0394    img0 = map_img(img0, opt.elem_output);
0395    img.elem_data = img.elem_data - img0.elem_data;
0396 end
0397 img.meas_err = dv;
0398 if opt.return_working_variables
0399   img.inv_solve_core.J = J;
0400   img.inv_solve_core.dx = dx;
0401   img.inv_solve_core.sx = sx;
0402   img.inv_solve_core.alpha = alpha_all;
0403   img.inv_solve_core.beta = beta_all;
0404   img.inv_solve_core.k = k;
0405   img.inv_solve_core.r = r(1:(k+1),:); % trim r to n-iterations' rows
0406   img.inv_solve_core.N = N;
0407   img.inv_solve_core.W = W;
0408   img.inv_solve_core.hp2RtR = hp2RtR;
0409   img.inv_solve_core.dv = dv;
0410   img.inv_solve_core.de = de;
0411   if opt.verbose >= 5
0412     img.inv_solve_core.dvall = dvall;
0413   end
0414 end
0415 %img = data_mapper(img, 1); % move data from img.elem_data to whatever 'params'
0416 
0417 function show_meas_err(dvall, data0, k, N, W, opt)
0418    clf;
0419    subplot(211); bar(dvall(:,k+1)); ylabel(sprintf('dv_k [%s]',opt.meas_working)); xlabel('meas #'); title(sprintf('measurement error @ iter=%d',k));
0420    subplot(212); bar(map_meas(dvall(:,k+1),N,opt.meas_working, 'voltage')); ylabel('dv_k [V]'); xlabel('meas #'); title('');
0421    drawnow;
0422    if isfield(opt,'fig_prefix')
0423       print('-dpdf',sprintf('%s-meas_err%d',opt.fig_prefix,k));
0424       print('-dpng',sprintf('%s-meas_err%d',opt.fig_prefix,k));
0425       saveas(gcf,sprintf('%s-meas_err%d.fig',opt.fig_prefix,k));
0426    end
0427    drawnow;
0428 
0429 function img = init_elem_data(img, opt)
0430   if opt.verbose > 1
0431     fprintf('  setting prior elem_data\n');
0432   end
0433   ne2 = 0; % init
0434   img.elem_data = zeros(sum([opt.elem_len{:}]),1); % preallocate
0435   for i=1:length(opt.elem_prior)
0436     ne1 = ne2+1; % next start idx ne1
0437     ne2 = ne1+opt.elem_len{i}-1; % this set ends at idx ne2
0438     if opt.verbose > 1
0439       if length(opt.prior_data{i}) == 1
0440         fprintf('    %d x %s: %0.1f\n',opt.elem_len{i},opt.elem_prior{i}, opt.prior_data{i});
0441       else
0442         fprintf('    %d x %s: ...\n',opt.elem_len{i},opt.elem_prior{i});
0443         if length(opt.prior_data{i}) ~= opt.elem_len{i}
0444            error(sprintf('expected %d elem, got %d elem in elem_prior', ...
0445                          opt.elem_len{i}, length(opt.prior_data{i})));
0446         end
0447       end
0448     end
0449     img.params_sel(i) = {ne1:ne2};
0450     img.elem_data(img.params_sel{i}) = opt.prior_data{i};
0451   end
0452   img.current_params = opt.elem_prior;
0453 
0454 function W = init_meas_icov(inv_model, opt)
0455    W = 1;
0456    if opt.calc_meas_icov
0457       if opt.verbose > 1
0458          disp('  calc measurement inverse covariance W');
0459       end
0460       W   = calc_meas_icov( inv_model );
0461    end
0462    err_if_inf_or_nan(W, 'init_meas_icov');
0463 
0464 function [N, dN] = init_normalization(fmdl, data0, opt)
0465    % precalculate the normalization of the data if required (apparent resistivity)
0466    N = 1;
0467    dN = 1;
0468    vh1.meas = 1;
0469    if ~ischar(opt.meas_input) || ~ischar(opt.meas_working)
0470       error('expected strings for meas_input and meas_working');
0471    end
0472    go =       any(strcmp({opt.meas_input, opt.meas_working},'apparent_resistivity'));
0473    go = go || any(strcmp({opt.meas_input, opt.meas_working},'log_apparent_resistivity'));
0474    go = go || any(strcmp({opt.meas_input, opt.meas_working},'log10_apparent_resistivity'));
0475    if go
0476       if opt.verbose > 1
0477          disp(['  calc measurement normalization matrix N (voltage -> ' opt.meas_working ')']);
0478       end
0479       % calculate geometric factor for apparent_resitivity conversions
0480       img1 = mk_image(fmdl,1);
0481       vh1  = fwd_solve(img1);
0482       N    = spdiag(1./vh1.meas);
0483       err_if_inf_or_nan(N,  'init_normalization: N');
0484    end
0485    if go && (opt.verbose > 1)
0486       disp(['  calc Jacobian normalization matrix   dN (voltage -> ' opt.meas_working ')']);
0487    end
0488    % calculate the normalization factor for the Jacobian
0489    data0 = map_meas_struct(data0, N, 'voltage'); % to voltage
0490    switch opt.meas_working
0491       case 'apparent_resistivity'
0492          dN = da_dv(data0.meas, vh1.meas);
0493       case 'log_apparent_resistivity'
0494          dN = dloga_dv(data0.meas, vh1.meas);
0495       case 'log10_apparent_resistivity'
0496          dN = dlog10a_dv(data0.meas, vh1.meas);
0497       case 'voltage'
0498          dN = dv_dv(data0.meas, vh1.meas);
0499       case 'log_voltage'
0500          dN = dlogv_dv(data0.meas, vh1.meas);
0501       case 'log10_voltage'
0502          dN = dlog10v_dv(data0.meas, vh1.meas);
0503       otherwise
0504          error('hmm');
0505    end
0506    err_if_inf_or_nan(dN, 'init_normalization: dN');
0507 
0508 % r_km1: previous residual, if its the first iteration r_km1 = inf
0509 % r_k: new residual
0510 function [stop, k, r, img] = update_residual(dv, img, de, W, hp2RtR, k, r, alpha, sx, opt)
0511   stop = 0;
0512 
0513   % update residual estimate
0514   if k == 0
0515      r = ones(opt.max_iterations, 3)*NaN;
0516      r_km1 = inf;
0517      r_1   = inf;
0518   else
0519      r_km1 = r(k, 1);
0520      r_1   = r(1, 1);
0521   end
0522   [r_k m_k e_k] = feval(opt.residual_func, dv, de, W, hp2RtR);
0523   % save residual for next iteration
0524   r(k+1,1:3) = [r_k m_k e_k];
0525 
0526   % now do something with that information
0527   if opt.verbose > 1
0528      fprintf('    calc residual\n');
0529      if k == 0
0530         fprintf('    stop @ max iter = %d, tol = %0.3g (%0.3g%%), dtol = %0.3g%% (after %d iter)\n', ...
0531                 opt.max_iterations, opt.tol, opt.tol/r_k*100, opt.dtol*100, opt.dtol_iter);
0532         fprintf('      r =%0.3g = %0.3g meas + %0.3g elem\n', r_k, m_k, e_k);
0533      else
0534         fprintf('      r =%0.3g (%0.03g%%) = %0.3g meas (%0.03g%%) + %0.3g elem (%0.3g%%)\n', ...
0535                 r_k, r_k/r(1,1)*100, m_k, m_k/r(1,2)*100, e_k, e_k/r(1,3)*100);
0536         dr = (r_k - r_km1);
0537         fprintf('      dr=%0.3g (%0.3g%%)\n', dr, dr/r_1*100);
0538      end
0539   end
0540   if opt.plot_residuals
0541      %         optimization_criteria, data misfit, roughness
0542      r(k+1,2:3) = [(dv'*dv)/2 (de'*de)/2];
0543      if k > 0
0544         clf;
0545         x = 1:(k+1);
0546         y = r(x, :);
0547         y = y ./ repmat(max(y,[],1),size(y,1),1) * 100;
0548         plot(x-1, y, 'o-', 'linewidth', 2, 'MarkerSize', 10);
0549         title('residuals');
0550         axis tight;
0551         ylabel('residual (% of max)');
0552         xlabel('iteration');
0553         set(gca, 'xtick', x);
0554         set(gca, 'xlim', [0 max(x)-1]);
0555         legend('residual','meas. misfit','prior misfit');
0556         legend('Location', 'EastOutside');
0557         drawnow;
0558         if isfield(opt,'fig_prefix')
0559            print('-dpdf',sprintf('%s-r',opt.fig_prefix));
0560            print('-dpng',sprintf('%s-r',opt.fig_prefix));
0561            saveas(gcf,sprintf('%s-r.fig',opt.fig_prefix));
0562         end
0563      end
0564   end
0565 
0566   % evaluate stopping criteria
0567   if r_k > r_km1 % bad step
0568      if opt.verbose > 1
0569         fprintf('  terminated at iteration %d (bad step, returning previous iteration''s result)\n',k);
0570      end
0571      img.elem_data = img.elem_data - alpha * sx; % undo the last step
0572      stop = -1;
0573   elseif k >= opt.max_iterations
0574      if opt.verbose > 1
0575         fprintf('  terminated at iteration %d (max iterations)\n',k);
0576      end
0577      stop = 1;
0578   elseif r_k < opt.tol + opt.ntol
0579      if opt.verbose > 1
0580         fprintf('  terminated at iteration %d\n',k);
0581         fprintf('    residual tolerance (%0.3g) achieved\n', opt.tol + opt.ntol);
0582      end
0583      stop = 1;
0584   elseif (k >= opt.dtol_iter) && ((r_k - r_km1)/r_1 > opt.dtol - 2*opt.ntol)
0585      if opt.verbose > 1
0586         fprintf('  terminated at iteration %d (iterations not improving)\n', k);
0587         fprintf('    residual slope tolerance (%0.3g%%) exceeded\n', (opt.dtol - 2*opt.ntol)*100);
0588      end
0589      stop = 1;
0590   end
0591   if ~stop
0592      % update iteration count
0593      k = k+1;
0594   end
0595 
0596 % for Conjugate Gradient, else beta = 0
0597 %  dx_k, dx_{k-1}, sx_{k-1}
0598 function beta = update_beta(dx_k, dx_km1, sx_km1, opt);
0599    if isfield(opt, 'beta_func')
0600       if opt.verbose > 1
0601          try beta_str = func2str(opt.beta_func);
0602          catch
0603             try 
0604                 beta_str = opt.beta_func;
0605             catch
0606                 beta_str = 'unknown';
0607             end
0608          end
0609       end
0610       beta= feval(opt.beta_func, dx_k, dx_km1, sx_km1);
0611    else
0612      beta_str = '<none>';
0613      beta = 0;
0614    end
0615    if opt.verbose > 1
0616       str = sprintf('    calc beta (%s)=%0.3f\n', beta_str, beta);
0617    end
0618 
0619 % update the search direction
0620 % for Gauss-Newton
0621 %   sx_k = dx_k
0622 % for Conjugate-Gradient
0623 %   sx_k = dx_k + beta * sx_{k-1}
0624 function sx = update_sx(dx, beta, sx_km1, opt);
0625    sx = dx + beta * sx_km1;
0626    if(opt.verbose > 1)
0627       nsx = norm(sx);
0628       nsxk = norm(sx_km1);
0629       fprintf( '    update step dx, beta=%0.3g, ||dx||=%0.3g\n', beta, nsx);
0630       if nsxk ~= 0
0631          fprintf( '      acceleration     d||dx||=%0.3g\n', nsx-nsxk);
0632          % ||ddx|| = chord_len = 2 sin(theta/2)
0633          ddx = norm(sx/nsx-sx_km1/nsxk);
0634          fprintf( '      direction change ||ddx||=%0.3g (%0.3g°)\n', ddx, 2*asind(ddx/2)); 
0635       end
0636    end
0637 
0638 % this function constructs the blockwise RtR matrix
0639 function hp2RtR = update_hp2RtR(inv_model, J, k, img, opt)
0640    if k==0 % first the start up iteration use the initial hyperparameter
0641       k=1;
0642    end
0643    % TODO sometimes (with Noser?) this requires the Jacobian, could this be done more efficiently?
0644    % add a test function to determine if img.elem_data affects RtR, skip this if independant
0645    % TODO we could detect in the opt_parsing whether the calc_RtR_prior depends on 'x' and skip this if no
0646    if ~opt.calc_RtR_prior
0647       error('no RtR calculation mechanism, set imdl.inv_solve_core.RtR_prior or imdl.RtR_prior');
0648    end
0649    if opt.verbose > 1
0650       disp('    calc hp^2 R^t R');
0651    end
0652    hp2 = calc_hyperparameter( inv_model )^2; % = \lambda^2
0653    net = sum([opt.elem_len{:}]); % Number of Elements, Total
0654    RtR = sparse(net,net); % init RtR = sparse(zeros(net,net));
0655    esi = 0; eei = 0; % element start, element end
0656    for i = 1:size(opt.RtR_prior,1) % row i
0657       esi = eei + 1;
0658       eei = eei + opt.elem_len{i};
0659       esj = 0; eej = 0; % element start, element end
0660       for j = 1:size(opt.RtR_prior,2) % column j
0661          esj = eej + 1;
0662          eej = eej + opt.elem_len{j};
0663          if isempty(opt.RtR_prior{i,j}) % null entries
0664             continue; % no need to explicitly create zero block matrices
0665          end
0666 
0667          % select a hyperparameter, potentially, per iteration
0668          % if we're at the end of the list, select the last entry
0669          hp=opt.hyperparameter{i,j};
0670          if length(hp) > k
0671             hp=hp(k);
0672          else
0673             hp=hp(end);
0674          end
0675 
0676          try RtR_str = func2str(opt.RtR_prior{i,j});
0677          catch
0678             try 
0679                 RtR_str = opt.RtR_prior{i,j};
0680             catch
0681                 RtR_str = 'unknown';
0682             end
0683          end
0684          if opt.verbose > 1
0685             fprintf('      {%d,%d} regularization RtR (%s), ne=%dx%d, hp=%0.4g\n', i,j,RtR_str,eei-esi+1,eej-esj+1,hp*sqrt(hp2));
0686          end
0687          imgt = map_img(img, opt.elem_working{i});
0688          inv_modelt = inv_model;
0689          inv_modelt.RtR_prior = opt.RtR_prior{i,j};
0690          if strcmp(RtR_str, 'prior_noser') % intercept prior_noser: we already have the Jacobian calculated
0691             if i ~= j
0692                error('noser for off diagonal prior (RtR) is undefined')
0693             end
0694             exponent= 0.5; try exponent = inv_model.prior_noser.exponent; end; try exponent = inv_model.prior_noser.exponent{j}; end
0695             ss = sum(J(:,esj:eej).^2,1)';
0696             Reg = spdiags( ss.^exponent, 0, opt.elem_len{j}, opt.elem_len{j}); % = diag(J'*J)^0.5
0697             RtR(esi:eei, esj:eej) = hp.^2 * Reg;
0698          else
0699             RtR(esi:eei, esj:eej) = hp.^2 * calc_RtR_prior_wrapper(inv_modelt, imgt, opt);
0700          end
0701       end
0702    end
0703    hp2RtR = hp2*RtR;
0704 
0705 function plot_dx_and_svd_elem(J, W, hp2RtR, k, sx, dx, img, opt)
0706    if(opt.verbose >= 5)
0707       % and try a show_fem with the pixel search direction
0708       clf;
0709       imgb=img;
0710       imgb.elem_data = dx;
0711       imgb.current_params = opt.elem_working;
0712       imgb.is_dx_plot = 1; % hint to show_fem that this is a dx plot (for handling log scale if anything clever is going on)
0713       if isfield(imgb.inv_model,'rec_model')
0714          imgb.fwd_model = imgb.inv_model.rec_model;
0715       end
0716       feval(opt.show_fem,imgb,1);
0717       title(sprintf('dx @ iter=%d',k));
0718       drawnow;
0719       if isfield(opt,'fig_prefix')
0720          print('-dpng',sprintf('%s-dx%d',opt.fig_prefix,k));
0721          print('-dpdf',sprintf('%s-dx%d',opt.fig_prefix,k));
0722          saveas(gcf,sprintf('%s-dx%d.fig',opt.fig_prefix,k));
0723       end
0724    end
0725    if opt.verbose < 8
0726       return; % do nothing if not verbose
0727    end
0728    % canonicalize the structures so we don't have to deal with a bunch of scenarios below
0729    if ~isfield(img, 'params_sel')
0730       img.params_sel = {1:length(img.elem_data)};
0731    end
0732    if ~isfield(img, 'current_params')
0733       img.current_params = 'conductivity';
0734    end
0735    if ~iscell(img.current_params)
0736       img.current_params = {img.current_params};
0737    end
0738    % go
0739    cols=length(opt.elem_working);
0740    if norm(sx - dx) < range(dx)/max(dx)*0.01 % sx and dx are within 1%
0741       rows=2;
0742    else
0743       rows=3;
0744    end
0745    clf; % individual SVD plots
0746    for i=1:cols
0747       if 1 % if Tikhonov
0748          hp=opt.hyperparameter{i};
0749          if k ~= 0 && k < length(hp)
0750             hp = hp(k);
0751          else
0752             hp = hp(end);
0753          end
0754       else
0755          hp = [];
0756       end
0757       sel=img.params_sel{i};
0758       str=strrep(img.current_params{i},'_',' ');
0759       plot_svd(J(:,sel), W, hp2RtR(sel,sel), k, hp); xlabel(str);
0760       drawnow;
0761       if isfield(opt,'fig_prefix')
0762          print('-dpdf',sprintf('%s-svd%d-%s',opt.fig_prefix,k,img.current_params{i}));
0763          print('-dpng',sprintf('%s-svd%d-%s',opt.fig_prefix,k,img.current_params{i}));
0764          saveas(gcf,sprintf('%s-svd%d-%s.fig',opt.fig_prefix,k,img.current_params{i}));
0765       end
0766    end
0767    clf; % combo plot
0768    for i=1:cols
0769       if 1 % if Tikhonov
0770          hp=opt.hyperparameter{i};
0771          if k ~= 0 && k < length(hp)
0772             hp = hp(k);
0773          else
0774             hp = hp(end);
0775          end
0776       else
0777          hp = [];
0778       end
0779       subplot(rows,cols,i);
0780       sel=img.params_sel{i};
0781       str=strrep(img.current_params{i},'_',' ');
0782       plot_svd(J(:,sel), W, hp2RtR(sel,sel), k, hp); xlabel(str);
0783       subplot(rows,cols,cols+i);
0784       bar(dx(sel)); ylabel(['dx: ' str]);
0785       if rows > 2
0786          subplot(rows,cols,2*cols+i);
0787          bar(sx(sel)); ylabel(['sx: ' str]);
0788       end
0789    end
0790    drawnow;
0791    if isfield(opt,'fig_prefix')
0792       print('-dpdf',sprintf('%s-svd%d',opt.fig_prefix,k));
0793       print('-dpng',sprintf('%s-svd%d',opt.fig_prefix,k));
0794       saveas(gcf,sprintf('%s-svd%d.fig',opt.fig_prefix,k));
0795    end
0796 
0797 function plot_svd(J, W, hp2RtR, k, hp)
0798    if nargin < 5
0799       hp = [];
0800    end
0801    % calculate the singular values before and after regularization
0802    [~,s1,~]=svd(J'*W*J); s1=sqrt(diag(s1));
0803    [~,s2,~]=svd(J'*W*J + hp2RtR); s2=sqrt(diag(s2));
0804    h=semilogy(s1,'bx'); axis tight; set(h,'LineWidth',2);
0805    hold on; h=semilogy(s2,'go'); axis tight; set(h,'LineWidth',2); hold off;
0806    xlabel('k'); ylabel('value \sigma');
0807    title(sprintf('singular values of J at iteration %d',k));
0808    legend('J^T J', 'J^T J + \lambda^2 R^T R'); legend location best;
0809    % line for \lambda
0810 %     if regularization == 2 % Noser
0811 %        hp_scaled = hp*sqrt(norm(full(RtR)));
0812 %        h=line([1 length(s1)],[hp_scaled hp_scaled]);
0813 %        text(length(s1)/2,hp_scaled*0.9,sprintf('\\lambda ||R^T R||^{%0.1f}= %0.4g; \\lambda = %0.4g', noser_p, hp_scaled, hp));
0814 %        fprintf('  affecting %d of %d singular values\k', length(find(s1<hp_scaled)), length(s1));
0815 %     else % Tikhonov
0816    if length(hp)==1
0817         h=line([1 length(s1)],[hp hp]);
0818         ly=10^(log10(hp)-0.05*range(log10([s1;s2])));
0819         text(length(s1)/2,ly,sprintf('\\lambda = %0.4g', hp));
0820 %       fprintf('  affecting %d of %d singular values\k', length(find(s1<hp)), length(s1));
0821      end
0822      set(h,'LineStyle','-.'); set(h,'LineWidth',2);
0823    set(gca,'YMinorTick','on', 'YMinorGrid', 'on', 'YGrid', 'on');
0824 
0825 
0826 % TODO this function is one giant HACK around broken RtR generation with c2f matrices
0827 function RtR = calc_RtR_prior_wrapper(inv_model, img, opt)
0828    RtR = calc_RtR_prior( inv_model );
0829    if size(RtR,1) < length(img.elem_data)
0830      ne = length(img.elem_data) - size(RtR,1);
0831      % we are correcting for the added background element
0832      for i=1:ne
0833        RtR(end+1:end+1, end+1:end+1) = RtR(1,1);
0834      end
0835      if opt.verbose > 1
0836         fprintf('      c2f: adjusting RtR by appending %d rows/cols\n', ne);
0837         disp(   '      TODO move this fix, or something like it to calc_RtR_prior -- this fix is a quick HACK to get things to run...');
0838      end
0839    end
0840 
0841 % opt is only updated for the fwd_solve count
0842 function [J, opt] = update_jacobian(img, dN, k, opt)
0843    base_types = map_img_base_types(img);
0844    imgb = map_img(img, base_types);
0845    imgb = feval(opt.update_img_func, imgb, opt);
0846    % if the electrodes/geometry moved, we need to recalculate dN if it depends on vh
0847    % note that only apparent_resisitivity needs vh; all others depend on data0 measurements
0848    if any(strcmp(map_img_base_types(img), 'movement')) && any(strcmp(opt.meas_working, 'apparent_resistivity'))
0849       imgh = map_img(imgb, 'conductivity'); % drop everything but conductivity
0850       imgh.elem_data = imgh.elem_data*0 +1; % conductivity = 1
0851       vh = fwd_solve(imgh); vh = vh.meas;
0852       dN = da_dv(1,vh); % = diag(1/vh)
0853       opt.fwd_solutions = opt.fwd_solutions +1;
0854    end
0855    ee = 0; % element select, init
0856    pp = fwd_model_parameters(imgb.fwd_model);
0857    J = zeros(pp.n_meas,sum([opt.elem_len{:}]));
0858    for i=1:length(opt.jacobian)
0859       if(opt.verbose > 1)
0860          if isnumeric(opt.jacobian{i})
0861             J_str = '[fixed]';
0862          elseif isa(opt.jacobian{i}, 'function_handle')
0863             J_str = func2str(opt.jacobian{i});
0864          else
0865             J_str = opt.jacobian{i};
0866          end
0867          if i == 1 fprintf('    calc Jacobian J(x) = ');
0868          else      fprintf('                       + '); end
0869          fprintf('(%s,', J_str);
0870       end
0871       % start and end of these Jacobian columns
0872       es = ee+1;
0873       ee = es+opt.elem_len{i}-1;
0874       % scaling if we are working in something other than direct conductivity
0875       S = feval(opt.calc_jacobian_scaling_func{i}, imgb.elem_data(es:ee)); % chain rule
0876       % finalize the jacobian
0877       % Note that if a normalization (i.e. apparent_resistivity) has been applied
0878       % to the measurements, it needs to be applied to the Jacobian as well!
0879       imgt = imgb;
0880       if  strcmp(base_types{i}, 'conductivity') % make legacy jacobian calculators happy... only conductivity on imgt.elem_data
0881          imgt = map_img(img, 'conductivity');
0882       end
0883       imgt.fwd_model.jacobian = opt.jacobian{i};
0884       Jn = calc_jacobian( imgt ); % unscaled natural units (i.e. conductivity)
0885       J(:,es:ee) = dN * Jn * S; % scaled and normalized
0886       if opt.verbose > 1
0887          tmp = zeros(1,size(J,2));
0888          tmp(es:ee) = 1;
0889          tmp(opt.elem_fixed) = 0;
0890          fprintf(' %d DoF, %d meas, %s)\n', sum(tmp), size(J,1), func2str(opt.calc_jacobian_scaling_func{i}));
0891       end
0892       if opt.verbose >= 5
0893          clf;
0894          t=axes('Position',[0 0 1 1],'Visible','off'); % something to put our title on after we're done
0895          text(0.03,0.1,sprintf('update\\_jacobian (%s), iter=%d', strrep(J_str,'_','\_'), k),'FontSize',20,'Rotation',90);
0896          for y=0:1
0897             if y == 0; D = Jn; else D = J(:,es:ee); end
0898             axes('units', 'normalized', 'position', [ 0.13 0.62-y/2 0.8 0.3 ]);
0899             imagesc(D);
0900             if y == 0; ylabel('meas (1)'); xlabel(['elem (' strrep(base_types{i},'_','\_') ')']);
0901             else       ylabel('meas (dN)'); xlabel(['elem (' strrep(opt.elem_working{i},'_','\_') ')']);
0902             end
0903             os = get(gca, 'Position'); c=colorbar('southoutside'); % colorbar start...
0904             set(gca, 'Position', os); % fix STUPID colorbar resizing
0905             % reduce height, this has to be done after the axes fix or STUPID matlab messes things up real good
0906             cP = get(c,'Position'); set(c,'Position', [0.13    0.54-y/2    0.8    0.010]);
0907             axes('units', 'normalized', 'position', [ 0.93 0.62-y/2 0.05 0.3 ]);
0908             barh(sqrt(sum(D.^2,2))); axis tight; axis ij; set(gca, 'ytick', [], 'yticklabel', []);
0909             axes('units', 'normalized', 'position', [ 0.13 0.92-y/2 0.8 0.05 ]);
0910             bar(sqrt(sum(D.^2,1))); axis tight; set(gca, 'xtick', [], 'xticklabel', []);
0911          end
0912          drawnow;
0913          if isfield(opt,'fig_prefix')
0914             print('-dpng',sprintf('%s-J%d-%s',opt.fig_prefix,k,strrep(J_str,'_','')));
0915             print('-dpdf',sprintf('%s-J%d-%s',opt.fig_prefix,k,strrep(J_str,'_','')));
0916             saveas(gcf,sprintf('%s-J%d-%s.fig',opt.fig_prefix,k,strrep(J_str,'_','')));
0917          end
0918       end
0919    end
0920    if opt.verbose >= 4
0921       % and try a show_fem with the pixel sensitivity
0922       try
0923          clf;
0924          imgb.elem_data = log(sqrt(sum(J.^2,1)));
0925          for i = 1:length(imgb.current_params)
0926             imgb.current_params{i} = [ 'log_sensitivity_' imgb.current_params{i} ];
0927          end
0928          if isfield(imgb.inv_model,'rec_model')
0929             imgb.fwd_model = imgb.inv_model.rec_model;
0930          end
0931          feval(opt.show_fem,imgb,1);
0932          title(sprintf('sensitivity @ iter=%d',k));
0933          drawnow;
0934          if isfield(opt,'fig_prefix')
0935             print('-dpng',sprintf('%s-Js%d-%s',opt.fig_prefix,k,strrep(J_str,'_','')));
0936             print('-dpdf',sprintf('%s-Js%d-%s',opt.fig_prefix,k,strrep(J_str,'_','')));
0937             saveas(gcf,sprintf('%s-Js%d-%s.fig',opt.fig_prefix,k,strrep(J_str,'_','')));
0938          end
0939       catch % no worries if it didn't work... carry on
0940       end
0941    end
0942 
0943 % -------------------------------------------------
0944 % Chain Rule Products for Jacobian Translations
0945 % x is conductivity, we want the chain rule to translate the
0946 % Jacobian of conductivity to conductivity on resistivity or
0947 % logs of either.
0948 % This chain rule works out to a constant.
0949 %
0950 % d log_b(x)     1          d x
0951 % ---------- = ------- , ---------- = x ln(b)
0952 %     d x      x ln(b)   d log_b(x)
0953 function S = dx_dlogx(x);
0954    S = spdiag(x);
0955 function S = dx_dlog10x(x);
0956    S = spdiag(x * log(10));
0957 % resistivity 'y'
0958 % d x     d x   -1
0959 % ----- = --- = ---, y = 1/x --> -(x^2)
0960 % d 1/x   d y   y^2
0961 function S = dx_dy(x);
0962    S = spdiag(-(x.^2));
0963 % then build the log versions of conductivity by combining chain rule products
0964 function S = dx_dlogy(x);
0965 %   S = dx_dy(x) * dy_dlogy(x);
0966 %     = -(x^2) * 1/x = -x
0967    S = spdiag(-x);
0968 function S = dx_dlog10y(x);
0969 %   S = dx_dy(x) * dy_dlog10y(x);
0970 %     = -(x^2) * 1/(ln(10) x) = -x / ln(10)
0971    S = spdiag(-x/log(10));
0972 % ... some renaming to make things understandable above: x = 1/y
0973 %function S = dy_dlogy(x);
0974 %   S = dx_dlogx(1./x);
0975 %function S = dy_dlog10y(x);
0976 %   S = dx_dlog10x(1./x);
0977 % -------------------------------------------------
0978 % apparent_resistivity 'a' versus voltage 'x'
0979 % d a    1  d v    1         v
0980 % --- = --- --- = --- ; a = ---
0981 % d v    vh d v    vh        vh
0982 % log_apparent_resistivity
0983 % d loga   d loga d a    1   1     vh  1     1
0984 % ------ = ------ --- = --- --- = --- --- = ---
0985 % d v       d a   d v    a   vh    v   vh    v
0986 function dN = da_dv(v,vh)
0987    dN = spdiag(1./vh); % N == dN for apparent_resistivity
0988 function dN = dloga_dv(v,vh)
0989    dN = spdiag(1./v);
0990 function dN = dlog10a_dv(v,vh)
0991    dN = spdiag( 1./(v * log(10)) );
0992 function dN = dv_dv(v,vh)
0993    dN = 1;
0994 function dN = dlogv_dv(v,vh) % same as dloga_dv
0995    dN = dloga_dv(v,vh);
0996 function dN = dlog10v_dv(v,vh) % same as dlog10a_dv
0997    dN = dlog10a_dv(v, vh);
0998 % -------------------------------------------------
0999 
1000 
1001 function [alpha, img, dv, opt] = update_alpha(img, sx, data0, img0, N, W, hp2RtR, k, dv, opt)
1002   if k == 0 % first iteration, just setting up, no line search happens
1003      alpha = 0;
1004      return;
1005   end
1006 
1007   if(opt.verbose > 1)
1008      try 
1009          ls_str = func2str(opt.line_search_func);
1010      catch
1011          ls_str = opt.line_search_func;
1012      end
1013      fprintf('    line search, alpha = %s\n', ls_str);
1014   end
1015 
1016   % some sanity checks before we feed this information to the line search
1017   err_if_inf_or_nan(sx, 'sx (pre-line search)');
1018   err_if_inf_or_nan(img.elem_data, 'img.elem_data (pre-line search)');
1019 
1020   if any(size(img.elem_data) ~= size(sx))
1021      error(sprintf('mismatch on elem_data[%d,%d] vs. sx[%d,%d] vector sizes, check c2f_background_fixed',size(img.elem_data), size(sx)));
1022   end
1023 
1024   optt = opt;
1025   if isfield(optt,'fig_prefix') % set fig_prefix for the iteration#
1026     optt.fig_prefix = [opt.fig_prefix '-k' num2str(k)];
1027   end
1028   [alpha, imgo, dv, opto] = feval(opt.line_search_func, img, sx, data0, img0, N, W, hp2RtR, dv, optt);
1029   if isfield(opt,'fig_prefix') % set fig_prefix back to it's old value
1030     opto.fig_prefix = opt.fig_prefix;
1031   end
1032   if ~isempty(imgo)
1033      img = imgo;
1034   else
1035      img.elem_data = img.elem_data + alpha*sx;
1036   end
1037   if ~isempty(opto)
1038      opt = opto;
1039   end
1040 
1041   if(opt.verbose > 1)
1042      fprintf('      selected alpha=%0.3g\n', alpha);
1043   end
1044 
1045   if (alpha == 0) && (k == 1)
1046     error('first iteration failed to advance solution');
1047   end
1048 
1049 function err_if_inf_or_nan(x, str);
1050   if any(any(isnan(x) | isinf(x)))
1051       error(sprintf('bad %s (%d NaN, %d Inf of %d)', ...
1052                     str, ...
1053                     length(find(isnan(x))), ...
1054                     length(find(isinf(x))), ...
1055                     length(x(:))));
1056   end
1057 
1058 
1059 function [img, dv] = update_img_using_limits(img, img0, data0, N, dv, opt)
1060   % fix max/min values for x
1061   if opt.max_value ~= +inf
1062      lih = find(img.elem_data > opt.max_value);
1063      img.elem_data(lih) = opt.max_value;
1064      if opt.verbose > 1
1065         fprintf('    limit max(x)=%g for %d elements\n', opt.max_value, length(lih));
1066      end
1067      dv = []; % dv is now invalid since we changed the conductivity
1068   end
1069   if opt.min_value ~= -inf
1070      lil = find(img.elem_data < opt.min_value);
1071      img.elem_data(lil) = opt.min_value;
1072      if opt.verbose > 1
1073         fprintf('    limit min(x)=%g for %d elements\n', opt.min_value, length(lil));
1074      end
1075      dv = []; % dv is now invalid since we changed the conductivity
1076   end
1077   % update voltage change estimate if the limit operation changed the img data
1078   [dv, opt] = update_dv(dv, img, data0, N, opt, '(dv out-of-date)');
1079 
1080 function  de = update_de(de, img, img0, opt)
1081    img0 = map_img(img0, opt.elem_working);
1082    img  = map_img(img,  opt.elem_working);
1083    err_if_inf_or_nan(img0.elem_data, 'de img0');
1084    err_if_inf_or_nan(img.elem_data,  'de img');
1085    % probably not the most robust check for whether this is the first update
1086    % but this ensures that we get exactly zero for the first iteration and not
1087    % a set of values that has numeric floating point errors that are nearly zero
1088    if isempty(de) % first iteration
1089       % data hasn't changed yet!
1090       de = zeros(size(img0.elem_data));
1091    else
1092       de = img0.elem_data - img.elem_data;
1093    end
1094    de(opt.elem_fixed) = 0; % TODO is this redundant... delete me?
1095    err_if_inf_or_nan(de, 'de out');
1096 
1097 function [dv, opt] = update_dv(dv, img, data0, N, opt, reason)
1098    % estimate current error as a residual
1099    if ~isempty(dv) % need to calculate dv...
1100       return;
1101    end
1102    if nargin < 7
1103       reason = '';
1104    end
1105    if opt.verbose > 1
1106       disp(['    fwd_solve b=Ax ', reason]);
1107    end
1108    [dv, opt, err] = update_dv_core(img, data0, N, opt);
1109 % TODO AB inject the img.error here, so it doesn't need to be recalculated when calc_solution_error=1
1110 %   img.error = err;
1111 
1112 function data = map_meas_struct(data, N, out)
1113    try   
1114        current_meas_params = data.current_params;
1115    catch
1116        current_meas_params = 'voltage';
1117    end
1118    data.meas = map_meas(data.meas, N, current_meas_params, out);
1119    data.current_params = out;
1120    err_if_inf_or_nan(data.meas, 'dv meas');
1121 
1122 % also used by the line search as opt.line_search_dv_func
1123 function [dv, opt, err] = update_dv_core(img, data0, N, opt)
1124    data0 = map_meas_struct(data0, N, 'voltage');
1125    img = map_img(img, map_img_base_types(img));
1126    img = feval(opt.update_img_func, img, opt);
1127    img = map_img(img, 'conductivity'); % drop everything but conductivity
1128    % if the electrodes/geometry moved, we need to recalculate N if it's being used
1129    if any(any(N ~= 1)) && any(strcmp(map_img_base_types(img), 'movement'))
1130       % note: data0 is mapped back to 'voltage' before N is modified
1131       imgh=img; imgh.elem_data = imgh.elem_data*0 +1; % conductivity = 1
1132       vh = fwd_solve(imgh); vh = vh.meas;
1133       N = spdiag(1./vh);
1134       opt.fwd_solutions = opt.fwd_solutions +1;
1135    end
1136    data = fwd_solve(img);
1137    opt.fwd_solutions = opt.fwd_solutions +1;
1138    dv = calc_difference_data(data0, data, img.fwd_model);
1139 %   clf;subplot(211); h=show_fem(img,1); set(h,'EdgeColor','none'); subplot(212); xx=1:length(dv); plot(xx,[data0.meas,data.meas,dv]); legend('d0','d1','dv'); drawnow; pause(1);
1140    if nargout >= 3
1141       err = norm(dv)/norm(data0.meas);
1142    else
1143       err = NaN;
1144    end
1145    dv = map_meas(dv, N, 'voltage', opt.meas_working);
1146    err_if_inf_or_nan(dv, 'dv out');
1147 
1148 function show_fem_iter(k, img, inv_model, stop, opt)
1149   if (opt.verbose < 4) || (stop == -1)
1150      return; % if verbosity is low OR we're dropping the last iteration because it was bad, do nothing
1151   end
1152   if opt.verbose > 1
1153      str=opt.show_fem;
1154      if isa(str,'function_handle')
1155         str=func2str(str);
1156      end
1157      disp(['    ' str '()']);
1158   end
1159   if isequal(opt.show_fem,@show_fem) % opt.show_fem == @show_fem... so we need to try to be smart enough to make show_fem not explode
1160      img = map_img(img, 'resistivity'); % TODO big fat hack to make this work at the expense of an actual function...
1161      [img, opt] = strip_c2f_background(img, opt, '    ');
1162      % check we're returning the right size of data
1163      if isfield(inv_model, 'rec_model')
1164        img.fwd_model = inv_model.rec_model;
1165      end
1166    %  bg = 1;
1167    %  img.calc_colours.ref_level = bg;
1168    %  img.calc_colours.clim = bg;
1169      img.calc_colours.cb_shrink_move = [0.3,0.6,0.02]; % move color bars
1170      if size(img.elem_data,1) ~= size(img.fwd_model.elems,1)
1171         warning(sprintf('img.elem_data has %d elements, img.fwd_model.elems has %d elems\n', ...
1172                         size(img.elem_data,1), ...
1173                         size(img.fwd_model.elems,1)));
1174      end
1175   else % do the "final clean up", same as when we quit
1176      img = map_img(img, opt.elem_output);
1177      [img, opt] = strip_c2f_background(img, opt, '    ');
1178      if isfield(inv_model, 'rec_model')
1179        img.fwd_model = inv_model.rec_model;
1180      end
1181   end
1182   clf; feval(opt.show_fem, img, 1);
1183   title(sprintf('x @ iter=%d',k));
1184   drawnow;
1185   if isfield(opt,'fig_prefix')
1186      print('-dpdf',sprintf('%s-x%d',opt.fig_prefix,k));
1187      print('-dpng',sprintf('%s-x%d',opt.fig_prefix,k));
1188      saveas(gcf,sprintf('%s-x%d.fig',opt.fig_prefix,k));
1189   end
1190 
1191 % TODO confirm that GN line_search_onm2 is using this residual calculation (preferably, directly)
1192 function [ residual meas elem ] = GN_residual(dv, de, W, hp2RtR)
1193 %   [size(dv); size(W); size(de); size(hp2RtR)]
1194    % we operate on whatever the iterations operate on (log data, resistance, etc) + perturb(i)*dx
1195    meas = 0.5*( dv'*W*dv);
1196    elem = 0.5*(de'*hp2RtR*de);
1197    residual = meas + elem;
1198 
1199 function residual = meas_residual(dv, de, W, hp2RtR)
1200    residual = norm(dv);
1201 
1202 %function img = initial_estimate( imdl, data )
1203 %   img = calc_jacobian_bkgnd( imdl );
1204 %   vs = fwd_solve(img);
1205 %
1206 %   if isstruct(data)
1207 %      data = data.meas;
1208 %   else
1209 %     meas_select = [];
1210 %     try
1211 %        meas_select = imdl.fwd_model.meas_select;
1212 %     end
1213 %     if length(data) == length(meas_select)
1214 %        data = data(meas_select);
1215 %     end
1216 %   end
1217 %
1218 %   pf = polyfit(data,vs.meas,1);
1219 %
1220 %   % create elem_data
1221 %   img = data_mapper(img);
1222 %
1223 %   if isfield(img.fwd_model,'coarse2fine');
1224 %      % TODO: the whole coarse2fine needs work here.
1225 %      %   what happens if c2f doesn't cover the whole region
1226 %
1227 %      % TODO: the two cases are very different. c2f case should match other
1228 %      nc = size(img.fwd_model.coarse2fine,2);
1229 %      img.elem_data = mean(img.elem_data)*ones(nc,1)*pf(1);
1230 %   else
1231 %      img.elem_data = img.elem_data*pf(1);
1232 %   end
1233 %
1234 %   % remove elem_data
1235 %%   img = data_mapper(img,1);
1236 %
1237 %function [img opt] = update_step(org, next, dx, fmin,res, opt)
1238 %   if isfield(opt, 'update_func')
1239 %      [img opt] = feval(opt.update_func,org,next,dx,fmin,res,opt);
1240 %   else
1241 %      img = next;
1242 %   end
1243 %
1244 % function bg = calc_background_resistivity(fmdl, va)
1245 %   % have a look at what we've created
1246 %   % compare data to homgeneous (how good is the model?)
1247 %   % NOTE background conductivity is set by matching amplitude of
1248 %   % homogeneous data against the measurements to get rough matching
1249 %   if(opt.verbose>1)
1250 %     fprintf('est. background resistivity\n');
1251 %   end
1252 %   cache_obj = { fmdl, va };
1253 %   BACKGROUND_R = eidors_obj('get-cache', cache_obj, 'calc_background_resistivity');
1254 %   if isempty(BACKGROUND_R);
1255 %     imgh = mk_image(fmdl, 1); % conductivity = 1 S/m
1256 %     vh = fwd_solve(imgh);
1257 %     % take the best fit of the data
1258 %     BACKGROUND_R = vh.meas \ va; % 32 Ohm.m ... agrees w/ Wilkinson's papers
1259 %     % update cache
1260 %     eidors_obj('set-cache', cache_obj, 'calc_background_resistivity', BACKGROUND_R);
1261 %   else
1262 %     if(opt.verbose > 1)
1263 %       fprintf('  ... cache hit\n');
1264 %     end
1265 %   end
1266 %   if(opt.verbose > 1)
1267 %     fprintf('estimated background resistivity: %0.1f Ohm.m\n', BACKGROUND_R);
1268 %   end
1269 
1270 function opt = parse_options(imdl)
1271    % merge legacy options locations
1272 %   imdl = deprecate_imdl_opt(imdl, 'parameters');
1273 %   imdl = deprecate_imdl_opt(imdl, 'inv_solve');
1274 
1275    % for any general options
1276    if isfield(imdl, 'inv_solve_core')
1277       opt = imdl.inv_solve_core;
1278    else
1279       opt = struct;
1280    end
1281 
1282    % verbosity, debug output
1283    % 0: quiet
1284    % 1: print iteration count
1285    % 2: print details as the algorithm progresses
1286    if ~isfield(opt,'verbose')
1287       opt.verbose = 1;
1288       fprintf('  selecting inv_model.inv_solve_core.verbosity=1\n');
1289    end
1290    if opt.verbose > 1
1291       fprintf('  setting default parameters\n');
1292    end
1293    % we track how many fwd_solves we do since they are the most expensive part of the iterations
1294    opt.fwd_solutions = 0;
1295 
1296    if ~isfield(opt, 'show_fem')
1297       opt.show_fem = @show_fem;
1298    end
1299 
1300    if ~isfield(opt, 'residual_func') % the objective function
1301       opt.residual_func = @GN_residual; % r = f(dv, de, W, hp2RtR)
1302       % NOTE: the meas_residual function exists to maintain
1303       % compatibility with Nolwenn's code, the GN_residual
1304       % is a better choice
1305       %opt.residual_func = @meas_residual; % [r,m,e] = f(dv, de, W, hp2RtR)
1306    end
1307 
1308    % calculation of update components
1309    if ~isfield(opt, 'update_func')
1310       opt.update_func = @GN_update; % dx = f(J, W, hp2RtR, dv, de, opt)
1311    end
1312    if ~isfield(opt, 'update_method')
1313       opt.update_method = 'lu';
1314    end
1315 
1316    % figure out if things need to be calculated
1317    if ~isfield(opt, 'calc_meas_icov') % derivative of the objective function
1318       opt.calc_meas_icov = 0; % W
1319    end
1320    if ~isfield(opt, 'calc_RtR_prior') % derivative of the objective function
1321       opt.calc_RtR_prior = 0; % RtR
1322    end
1323    if ~isfield(opt, 'calc_hyperparameter')
1324       opt.calc_hyperparameter = 0; % hp2
1325    end
1326 
1327 %   try
1328       if opt.verbose > 1
1329          fprintf('    examining function %s(...) for required arguments\n', func2str(opt.update_func));
1330       end
1331       % ensure that necessary components are calculated
1332       % opt.update_func: dx = f(J, W, hp2RtR, dv, de, opt)
1333 %TODO BROKEN      args = function_depends_upon(opt.update_func, 6);
1334       args = ones(4,1); % TODO BROKEN
1335       if args(2) == 1
1336          opt.calc_meas_icov = 1;
1337       end
1338       if args(3) == 1
1339          opt.calc_hyperparameter = 1;
1340       end
1341       if args(4) == 1
1342          opt.calc_RtR_prior = 1;
1343       end
1344 %   catch
1345 %      error('exploration of function %s via function_depends_upon() failed', func2str(opt.update_func));
1346 %   end
1347 
1348    % stopping criteria, solution limits
1349    if ~isfield(opt, 'max_iterations')
1350       opt.max_iterations = 10;
1351    end
1352    if ~isfield(opt, 'ntol')
1353       opt.ntol = eps; % attempt to quantify numeric machine precision
1354    end
1355    if ~isfield(opt, 'tol')
1356       opt.tol = 0; % terminate iterations if residual is less than tol
1357    end
1358    if ~isfield(opt, 'dtol')
1359       % terminate iterations if residual slope is greater than dtol
1360       % generally, we would want dtol to be -0.01 (1% decrease) or something similar
1361       %  ... as progress levels out, stop working
1362       opt.dtol = -1e-4; % --> -0.01% slope (really slow)
1363    end
1364    if opt.dtol > 0
1365       error('dtol must be less than 0 (residual decreases at every iteration)');
1366       % otherwise we won't converge...
1367    end
1368    if ~isfield(opt, 'dtol_iter')
1369       %opt.dtol_iter = inf; % ignore dtol for dtol_iter iterations
1370       opt.dtol_iter = 0; % use dtol from the beginning
1371    end
1372    if ~isfield(opt, 'min_value')
1373       opt.min_value = -inf; % min elem_data value
1374    end
1375    if ~isfield(opt, 'max_value')
1376       opt.max_value = +inf; % max elem_data value
1377    end
1378    % provide a graphical display of the line search values & fit
1379    if ~isfield(opt, 'plot_residuals')
1380       if opt.verbose > 2
1381          opt.plot_residuals = 1;
1382       else
1383          opt.plot_residuals = 0;
1384       end
1385    end
1386    if opt.plot_residuals ~= 0
1387       disp('  residual plot (updated per iteration) are enabled, to disable them set');
1388       disp('    inv_model.inv_solve_core.plot_residuals=0');
1389    end
1390 
1391    % line search
1392    if ~isfield(opt,'line_search_func')
1393       % [alpha, img, dv, opt] = f(img, sx, data0, img0, N, W, hp2RtR, dv, opt);
1394       opt.line_search_func = @line_search_onm2;
1395    end
1396    if ~isfield(opt,'line_search_dv_func')
1397       opt.line_search_dv_func = @update_dv_core;
1398       % [dv, opt] = update_dv_core(img, data0, N, opt)
1399    end
1400    if ~isfield(opt,'line_search_de_func')
1401       % we create an anonymous function to skip the first input argument since
1402       % we always want to calculate de in the line search
1403       opt.line_search_de_func = @(img, img0, opt) update_de(1, img, img0, opt);
1404       % de = f(img, img0, opt)
1405    end
1406    % an initial guess for the line search step sizes, may be modified by line search
1407    % TODO this 'sensible default' should be moved to the line_search code since it is not generic to any other line searches
1408    if ~isfield(opt,'line_search_args') || ...
1409       ~isfield(opt.line_search_args, 'perturb')
1410       fmin = 1/4; % arbitrary starting guess
1411       opt.line_search_args.perturb = [0 fmin/4 fmin/2 fmin fmin*2 fmin*4];
1412       %opt.line_search_args.perturb = [0 fmin/4 fmin fmin*4];
1413       %opt.line_search_args.perturb = [0 0.1 0.35 0.7 1.0];
1414       %opt.line_search_args.perturb = [0 0.1 0.5 0.7 1.0];
1415       %pt.line_search_args.perturb = [0 0.1 0.7 0.9 1.0];
1416       %opt.line_search_args.perturb = [0 0.1 0.9 1.0];
1417    end
1418    % provide a graphical display of the line search values & fit
1419    if ~isfield(opt,'line_search_args') || ...
1420       ~isfield(opt.line_search_args, 'plot')
1421       if opt.verbose >= 5
1422          opt.line_search_args.plot = 1;
1423       else
1424          opt.line_search_args.plot = 0;
1425       end
1426    end
1427    % pass fig_prefix to the line search as well, unless they are supposed to go somewhere elese
1428    if isfield(opt,'fig_prefix') && ...
1429       isfield(opt,'line_search_args') && ...
1430       ~isfield(opt.line_search_args, 'fig_prefix')
1431       opt.line_search_args.fig_prefix = opt.fig_prefix;
1432    end
1433    % some help on how to turn off the line search plots if we don't want to see them
1434    if opt.line_search_args.plot ~= 0
1435       disp('  line search plots (per iteration) are enabled, to disable them set');
1436       disp('    inv_model.inv_solve_core.line_search_args.plot=0');
1437    end
1438 
1439    % background
1440    % if > 0, this is the elem_data that holds the background
1441    % this is stripped off just before the iterations complete
1442    if ~isfield(opt, 'c2f_background')
1443      if isfield(imdl, 'fwd_model') && isfield(imdl.fwd_model, 'coarse2fine')
1444         opt.c2f_background = -1; % possible: check if its required later
1445      else
1446         opt.c2f_background = 0;
1447      end
1448    end
1449    if ~isfield(opt, 'c2f_background_fixed')
1450       opt.c2f_background_fixed = 1; % generally, don't touch the background
1451    end
1452 
1453 
1454    % DATA CONVERSION settings
1455    % elem type for the initial estimate is based on calc_jacobian_bkgnd which returns an img
1456    if ~isfield(opt, 'elem_working')
1457       opt.elem_working = {'conductivity'};
1458    end
1459    if ~isfield(opt, 'elem_prior')
1460       opt.elem_prior = {'conductivity'};
1461    end
1462    if ~isfield(opt, 'elem_output')
1463       opt.elem_output = {'conductivity'};
1464    end
1465    if ~isfield(opt, 'meas_input')
1466       opt.meas_input = 'voltage';
1467    end
1468    if ~isfield(opt, 'meas_working')
1469       opt.meas_working = 'voltage';
1470    end
1471    % if the user didn't put these into cell arrays, do
1472    % so here so there is less error checking later in
1473    % the code
1474    for i = {'elem_working', 'elem_prior', 'elem_output'}; %, 'meas_input', 'meas_working'}
1475      % MATLAB voodoo: deincapsulate a cell containing a
1476      % string, then use that to access a struct eleemnt
1477      x = opt.(i{1});
1478      if ~iscell(x)
1479         opt.(i{1}) = {x};
1480      end
1481    end
1482 
1483    if ~isfield(opt, 'prior_data')
1484       if isfield(imdl, 'jacobian_bkgnd') && ...
1485          isfield(imdl.jacobian_bkgnd, 'value') && ...
1486          length(opt.elem_prior) == 1
1487          opt.prior_data = {imdl.jacobian_bkgnd.value};
1488       else
1489          error('requires inv_model.inv_solve_core.prior_data');
1490       end
1491    end
1492 
1493    if ~isfield(opt, 'elem_len')
1494       if length(opt.elem_working) == 1
1495          if isfield(imdl.fwd_model, 'coarse2fine')
1496             c2f = imdl.fwd_model.coarse2fine; % coarse-to-fine mesh mapping
1497             opt.elem_len = { size(c2f,2) };
1498          else
1499             opt.elem_len = { size(imdl.fwd_model.elems,1) };
1500          end
1501       else
1502         error('requires inv_model.inv_solve_core.elem_len');
1503       end
1504    end
1505 
1506    % meas_select already handles selecting from the valid measurements
1507    % we want the same for the elem_data, so we only work on modifying the legal values
1508    % Note that c2f_background's elements are added to this list if opt.c2f_background_fixed == 1
1509    if ~isfield(opt, 'elem_fixed') % give a safe default, if none has been provided
1510       opt.elem_fixed = [];
1511    elseif iscell(opt.elem_fixed) % if its a cell-array, we convert it to absolute
1512      % numbers in elem_data
1513      %  -- requires: opt.elem_len to already be constructed if it was missing
1514       offset=0;
1515       ef=[];
1516       for i=1:length(opt.elem_fixed)
1517          ef = [ef, opt.elem_fixed{i} + offset];
1518          offset = offset + opt.elem_len{i};
1519       end
1520       opt.elem_fixed = ef;
1521    end
1522 
1523    % allow a cell array of jacobians
1524    if ~isfield(opt, 'jacobian')
1525       if ~iscell(imdl.fwd_model.jacobian)
1526          opt.jacobian = {imdl.fwd_model.jacobian};
1527       else
1528          opt.jacobian = imdl.fwd_model.jacobian;
1529       end
1530    elseif isfield(imdl.fwd_model, 'jacobian')
1531       imdl.fwd_model
1532       imdl
1533       error('inv_model.fwd_model.jacobian and inv_model.inv_solve_core.jacobian should not both exist: it''s ambiguous');
1534    end
1535    % defaul hyperparameter is 1
1536    if ~isfield(opt, 'hyperparameter')
1537       opt.hyperparameter = {[]};
1538       for i=1:length(opt.elem_working)
1539          opt.hyperparameter{i} = 1;
1540       end
1541    end
1542    % if the user didn't put these into cell arrays, do
1543    % so here so there is less error checking later in
1544    % the code
1545    for i = {'elem_len', 'prior_data', 'jacobian', 'hyperparameter'}
1546      % MATLAB voodoo: deincapsulate a cell containing a
1547      % string, then use that to access a struct eleemnt
1548      x = opt.(i{1});
1549      if ~iscell(x)
1550         opt.(i{1}) = {x};
1551      end
1552    end
1553    % show what the hyperparameters are configured to when logging
1554    if opt.verbose > 1
1555       fprintf('  hyperparameters\n');
1556       try 
1557           hp_global = imdl.hyperparameter.value;
1558           hp_global_str = sprintf(' x %0.4g',hp_global);
1559       catch
1560           hp_global = 1;
1561           hp_global_str = '';
1562       end
1563       for i=1:length(opt.elem_working)
1564          if isnumeric(opt.hyperparameter{i}) && length(opt.hyperparameter{i}) == 1
1565             fprintf('    %s: %0.4g\n',opt.elem_working{i}, opt.hyperparameter{i}*hp_global);
1566          elseif isa(opt.hyperparameter{i}, 'function_handle')
1567             fprintf('    %s: @%s%s\n',opt.elem_working{i}, func2str(opt.hyperparameter{i}), hp_global_str);
1568          elseif ischar(opt.hyperparameter{i})
1569             fprintf('    %s: @%s%s\n',opt.elem_working{i}, opt.hyperparameter{i}, hp_global_str);
1570          else
1571             fprintf('    %s: ...\n',opt.elem_working{i});
1572          end
1573       end
1574    end
1575 
1576    % REGULARIZATION RtR
1577    % for constructing the blockwise RtR matrix
1578    % can be: explicit matrix, blockwise matrix diagonal, or full blockwise matrix
1579    % blockwise matrices can be function ptrs or explicit
1580    if ~isfield(opt, 'RtR_prior')
1581       if isfield(imdl, 'RtR_prior')
1582          opt.RtR_prior = {imdl.RtR_prior};
1583       else
1584          opt.RtR_prior = {[]}; % null matrix (all zeros)
1585          warning('missing imdl.inv_solve_core.RtR_prior or imdl.RtR_prior: assuming NO regularization RtR=0');
1586       end
1587    end
1588    % bit of a make work project but if its actually a full numeric matrix we
1589    % canoncialize it by breaking it up into the blockwise components
1590    if isnumeric(opt.RtR_prior)
1591       if size(opt.RtR_prior, 1) ~= size(opt.RtR_prior, 2)
1592          error('expecting square matrix for imdl.RtR_prior or imdl.inv_solve_core.RtR_prior');
1593       end
1594       if length(opt.RtR_prior) == 1
1595          opt.RtR_prior = {opt.RtR_prior}; % encapsulate directly into a cell array
1596       else
1597          RtR = opt.RtR_prior;
1598          opt.RtR_prior = {[]};
1599          esi = 0; eei = 0;
1600          for i=1:length(opt.elem_len)
1601             esi = eei +1;
1602             eei = eei +opt.elem_len{i};
1603             esj = 0; eej = 0;
1604             for j=1:length(opt.elem_len)
1605                esj = eej +1;
1606                eej = eej +opt.elem_len{j};
1607                opt.RtR_prior(i,j) = RtR(esi:eei, esj:eej);
1608             end
1609          end
1610       end
1611    elseif ~iscell(opt.RtR_prior) % not a cell array? encapsulate it
1612       opt.RtR_prior = {opt.RtR_prior};
1613    end
1614    % if not square then expand the block matrix
1615    % single row/column: this is our diagonal --> expand to full blockwise matrix
1616    if any(size(opt.RtR_prior) ~= ([1 1]*length(opt.elem_len)))
1617       if (size(opt.RtR_prior, 1) ~= 1) && ...
1618          (size(opt.RtR_prior, 2) ~= 1)
1619          error('odd imdl.RtR_prior or imdl.inv_solve_core.RtR_prior, cannot figure out how to expand it blockwise');
1620       end
1621       if (size(opt.RtR_prior, 1) ~= length(opt.elem_len)) && ...
1622          (size(opt.RtR_prior, 2) ~= length(opt.elem_len))
1623          error('odd imdl.RtR_prior or imdl.inv_solve_core.RtR_prior, not enough blockwise components vs. elem_working types');
1624       end
1625       RtR_diag = opt.RtR_prior;
1626       opt.RtR_prior = {[]}; % delete and start again
1627       for i=1:length(opt.elem_len)
1628          opt.RtR_prior(i,i) = RtR_diag(i);
1629       end
1630    end
1631 
1632    % now sort out the hyperparameter for the "R^T R" (RtR) matrix
1633    hp=opt.hyperparameter;
1634    if size(hp,1) == 1 % one row
1635       hp = hp'; % ... now one col
1636    end
1637    if iscell(hp)
1638       % if it's a cell array that matches size of the RtR, then we're done
1639       if all(size(hp) == size(opt.RtR_prior))
1640          opt.hyperparameter = hp;
1641       % if the columns match, then we can expand on the diagonal, everything else gets '1'
1642       elseif length(hp) == length(opt.RtR_prior)
1643          opt.hyperparameter = opt.RtR_prior;
1644          [opt.hyperparameter{:}] = deal(1); % hp = 1 everywhere
1645          opt.hyperparameter(logical(eye(size(opt.RtR_prior)))) = hp; % assign to diagonal
1646       else
1647          error('hmm, don''t understand this opt.hyperparameter cellarray');
1648       end
1649    % if it's a single hyperparameter, that's the value everywhere
1650    elseif isnumeric(hp)
1651       opt.hyperparameter = opt.RtR_prior;
1652       [opt.hyperparameter{:}] = deal({hp});
1653    else
1654       error('don''t understand this opt.hyperparameter');
1655    end
1656 
1657    % JACOBIAN CHAIN RULE conductivity -> whatever
1658    % where x = conductivity at this iteration
1659    %       S = a scaling matrix, generally a diagonal matrix of size matching Jacobian columns
1660    % Jn = J * S;
1661    % if not provided, determine based on 'elem_working' type
1662    if ~isfield(opt, 'calc_jacobian_scaling_func')
1663       pinv = strfind(opt.elem_working, 'resistivity');
1664       plog = strfind(opt.elem_working, 'log_');
1665       plog10 = strfind(opt.elem_working, 'log10_');
1666       for i = 1:length(opt.elem_working)
1667         prefix = '';
1668         if plog{i}
1669            prefix = 'log';
1670         elseif plog10{i}
1671            prefix = 'log10';
1672         else
1673            prefix = '';
1674         end
1675         if pinv{i}
1676            prefix = [prefix '_inv'];
1677         end
1678         switch(prefix)
1679           case ''
1680              opt.calc_jacobian_scaling_func{i} = @ret1_func;  % S = f(x)
1681           case 'log'
1682              opt.calc_jacobian_scaling_func{i} = @dx_dlogx;   % S = f(x)
1683           case 'log10'
1684              opt.calc_jacobian_scaling_func{i} = @dx_dlog10x; % S = f(x)
1685           case '_inv'
1686              opt.calc_jacobian_scaling_func{i} = @dx_dy;      % S = f(x)
1687           case 'log_inv'
1688              opt.calc_jacobian_scaling_func{i} = @dx_dlogy;   % S = f(x)
1689           case 'log10_inv'
1690              opt.calc_jacobian_scaling_func{i} = @dx_dlog10y; % S = f(x)
1691           otherwise
1692              error('oops');
1693        end
1694      end
1695    end
1696 
1697    if ~isfield(opt, 'update_img_func')
1698       opt.update_img_func = @null_func; % img = f(img, opt)
1699    end
1700 
1701    if ~isfield(opt, 'return_working_variables')
1702       opt.return_working_variables = 0;
1703    end
1704 
1705 function check_matrix_sizes(J, W, hp2RtR, dv, de, opt)
1706    % assuming our equation looks something like
1707    % dx = (J'*W*J + hp2RtR)\(J'*dv + hp2RtR*de);
1708    % check that all the matrix sizes are correct
1709    ne = size(de,1);
1710    nv = size(dv,1);
1711    if size(de,2) ~= 1
1712       error('de cols (%d) not equal 1', size(de,2));
1713    end
1714    if size(dv,2) ~= 1
1715       error('dv cols (%d) not equal 1', size(dv,2));
1716    end
1717    if opt.calc_meas_icov && ...
1718       any(size(W) ~= [nv nv])
1719       error('W size (%d rows, %d cols) is incorrect (%d rows, %d cols)', size(W), nv, nv);
1720    end
1721    if any(size(J) ~= [nv ne])
1722       error('J size (%d rows, %d cols) is incorrect (%d rows, %d cols)', size(J), nv, ne);
1723    end
1724    if opt.calc_RtR_prior && ...
1725       any(size(hp2RtR) ~= [ne ne])
1726       error('hp2RtR size (%d rows, %d cols) is incorrect (%d rows, %d cols)', size(hp2RtR), ne, ne);
1727    end
1728 
1729 function dx = update_dx(J, W, hp2RtR, dv, de, opt)
1730    if(opt.verbose > 1)
1731       fprintf( '    calc step size dx\n');
1732    end
1733    % don't penalize for fixed elements
1734    de(opt.elem_fixed) = 0;
1735 
1736    % TODO move this outside the inner loop of the iterations, it only needs to be done once
1737    check_matrix_sizes(J, W, hp2RtR, dv, de, opt)
1738 
1739    % zero out the appropriate things so that we can get a dx=0 for the elem_fixed
1740    J(:,opt.elem_fixed) = 0;
1741    de(opt.elem_fixed,:) = 0;
1742    hp2RtR(opt.elem_fixed,:) = 0;
1743    V=opt.elem_fixed;
1744    N=size(hp2RtR,1)+1;
1745    hp2RtR(N*(V-1)+1) = 1; % set diagonals to 1 to avoid divide by zero
1746    % do the update step direction calculation
1747    dx = feval(opt.update_func, J, W, hp2RtR, dv, de, opt);
1748 
1749    % check that our elem_fixed stayed fixed
1750    if any(dx(opt.elem_fixed) ~= 0)
1751       error('elem_fixed did''t give dx=0 at update_dx')
1752    end
1753 
1754    if(opt.verbose > 1)
1755       fprintf('      ||dx||=%0.3g\n', norm(dx));
1756       es = 0; ee = 0;
1757       for i=1:length(opt.elem_working)
1758           es = ee +1; ee = ee + opt.elem_len{i};
1759           nd = norm(dx(es:ee));
1760           fprintf( '      ||dx_%d||=%0.3g (%s)\n',i, nd, opt.elem_working{i});
1761       end
1762    end
1763 
1764 function dx = GN_update(J, W, hp2RtR, dv, de, opt)
1765    if ~any(strcmp(opt.update_method, {'lu','pcg'}))
1766       error(['unsupported update_method: ',opt.update_method]);
1767    end
1768    if strcmp(opt.update_method, 'lu')
1769       try
1770          % the actual update
1771          dx = -(J'*W*J + hp2RtR)\(J'*W*dv + hp2RtR*de); % LU decomp
1772       catch ME % boom
1773          fprintf('      LU decomp failed: ');
1774          disp(ME.message)
1775          fprintf('      try opt.update_method = ''pcg'' on future runs\n');
1776          opt.update_method = 'pcg';
1777       end
1778    end
1779    if strcmp(opt.update_method, 'pcg')
1780       tol = 1e-6; % default 1e-6
1781       maxit = []; % default [] --> min(n,20)
1782       M = []; % default [] --> no preconditioner
1783       x0 = []; % default [] --> zeros(n,1)
1784 
1785       % try Preconditioned Conjugate Gradient: A x = b, solve for x
1786       % avoids J'*J for n x m matrix with large number of m cols --> J'*J becomes an m x m dense matrix
1787       LHS = @(x) (J'*(W*(J*x)) + hp2RtR*x);
1788       RHS = -(J'*W*dv + hp2RtR*de);
1789 
1790       tol=100*eps*size(J,2)^2; % rough estimate based on multiply-accumulates
1791 %      maxit = 10;
1792 
1793       [dx, flag, relres, iter, resvec] = pcg(LHS, RHS, tol, maxit, M', M', x0);
1794       % TODO if verbose...
1795       switch flag
1796          case 0
1797             if opt.verbose > 1
1798                fprintf('      PCG: relres=%g < tol=%g @ iter#%d\n',relres,tol,iter);
1799             end
1800          case 1
1801             if opt.verbose > 1
1802                fprintf('      PCG: relres=%g > tol=%g @ iter#%d (max#%d) [max iter]\n',relres,tol,iter,maxit);
1803             end
1804          case 2
1805             error('error: PCG ill-conditioned preconditioner M');
1806          case 3
1807             if opt.verbose > 1
1808                fprintf('      PCG: relres=%g > tol=%g @ iter#%d (max#%d) [stagnated]\n',relres,tol,iter,maxit);
1809             end
1810          case 4
1811             error('error: PCG a scalar quantity became too large or small to continue');
1812          otherwise
1813             error(sprintf('error: PCG unrecognized flag=%d',flag));
1814       end
1815 % NOTE PCG is still a work in progress and generally problem specific
1816 %      % plot convergence for pcg()
1817 %      clf;
1818 %         xlabel('iteration #');
1819 %         ylabel('relative residual');
1820 %         xx = 0:length(resvec)-1;
1821 %         semilogy(xx,resvec/norm(RHS),'b.');
1822 %         hold on;
1823 %         legend('no preconditioner');
1824 %hold on
1825 %% sparsity strategies: www.cerfacs.fr/algor/reports/Dissertations/TH_PA_02_48.pdf
1826 %sJ = J; sJ(J/max(J(:)) > 0.005) = 0; sJ=sparse(sJ); % sparsify J
1827 %M = ichol(sJ'*W*sJ + hp2RtR + speye(length(J)));
1828 %      [dx, flag, relres, iter, resvec] = pcg(LHS, RHS, tol, maxit, M);
1829 %         semilogy(xx,resvec/norm(RHS),'r.');
1830 %         legend('no P', 'IC(sp(J'')*W*sp(J) + hp2RtR)');
1831 %
1832 %
1833 %clf; Jj=J(:,1:800); imagesc(Jj'*Jj);
1834 
1835       % compare with gmres, bicgstab, lsqr
1836       % try preconditioners ilu, ichol (incomplete LU or Cholesky decomp.)
1837 
1838       %rethrow(ME); % we assume this is an 'excessive memory requested' failure
1839    end
1840 
1841 % for each argument, returns 1 if the function depends on it, 0 otherwise
1842 % 'zero' arguments do not need to be calculated since they don't get used
1843 function args = function_depends_upon(func, argn)
1844    % build function call
1845    str = sprintf('%s(',func2str(func));
1846    args = zeros(argn,1);
1847    for i = 1:argn-1
1848       str = [str sprintf('a(%d),',i)];
1849    end
1850    str = [str sprintf('a(%d))',argn)];
1851    % baseline
1852    a = ones(argn,1)*2;
1853    x = eval(str);
1854    % now check for a difference at each argument
1855    for i = 1:argn
1856       a = ones(argn,1)*2;
1857       a(i) = 0;
1858       y = eval(str);
1859       if any(x ~= y)
1860          args(i) = 1;
1861       end
1862    end
1863 
1864 % this function just passes data from its input to its output
1865 function out = null_func(in, varargin);
1866    out = in;
1867 
1868 % this function always returns one
1869 function [out, x, y, z] = ret1_func(varargin);
1870    out = 1;
1871    x = [];
1872    y = [];
1873    z = [];
1874 
1875 % if required, expand the coarse-to-fine matrix to cover the background of the image
1876 % this is removed at the end of the iterations
1877 function [inv_model, opt] = append_c2f_background(inv_model, opt)
1878     % either there is already a background
1879     % or none is required --> -1 means we go to work building one
1880     if opt.c2f_background >= 0
1881       return
1882     end
1883     % check that all elements get assigned a conductivity
1884     % through the c2f conversion
1885     c2f = inv_model.fwd_model.coarse2fine; % coarse-to-fine mesh mapping
1886     nf = size(inv_model.fwd_model.elems,1); % number of fine elements
1887     nc = size(c2f,2); % number of coarse elements
1888     % ... each element of fel aught to sum to '1' since the
1889     % elem_data is being assigned from a continuous unit
1890     % surface value
1891     % now, find any fine elements that are not fully
1892     % mapped between the two meshes (<1) w/in a tolerance
1893     % related to the number of additions in the summation
1894     fel = sum(c2f,2); % collapse mapping onto the fwd_model elements
1895     n = find(fel < 1 - (1e3+nc)*eps);
1896     % ... 1e3 is a fudge factor since we don't care too much
1897     %     about small area mapping errors
1898     % if we do have some unassigned elements,
1899     % expand c2f and add a background element to the 'elem_data'
1900     if length(n) ~= 0
1901       if(opt.verbose > 1)
1902         fprintf('  c2f: adding background conductivity to %d\n    fwd_model elements not covered by rec_model\n', length(n));
1903       end
1904       c2f(n,nc+1) = 1 - fel(n);
1905       inv_model.fwd_model.coarse2fine = c2f;
1906       opt.c2f_background = nc+1;
1907       if opt.c2f_background_fixed
1908          % TODO assumes conductivity/resistivity is the *first* parameterization
1909          opt.elem_fixed(end+1) = nc+1;
1910       end
1911     end
1912     % TODO assumes conductivity/resistivity is the *first* parameterization
1913     opt.elem_len(1) = {size(c2f,2)}; % elem_len +1
1914 
1915 function [img, opt] = strip_c2f_background(img, opt, indent)
1916     if nargin < 3
1917        indent = '';
1918     end
1919     % nothing to do?
1920     if opt.c2f_background <= 0
1921       return;
1922     end
1923 
1924     % if there are multiple 'params' (parametrizations), we assume its the first
1925     % TODO -- this isn't a great assumption but it'll work for now,
1926     %         we should add a better (more general) mechanism
1927     in = img.current_params;
1928     out = opt.elem_output;
1929     if iscell(in)
1930        in = in{1};
1931     end
1932     if iscell(out)
1933        out = out{1};
1934     end
1935 
1936     % go about cleaning up the background
1937     e = opt.c2f_background;
1938     % take backgtround elements and convert to output
1939     % 'params' (resistivity, etc)
1940     bg = map_data(img.elem_data(e), in, out);
1941     img.elem_data_background = bg;
1942     % remove elements from elem_data & c2f
1943     img.elem_data(e) = [];
1944     img.fwd_model.coarse2fine(:,e) = [];
1945     % remove our element from the lists
1946     opt.c2f_background = 0;
1947     ri = find(opt.elem_fixed == e);
1948     opt.elem_fixed(ri) = [];
1949     if isfield(img, 'params_sel')
1950        for i = 1:length(img.params_sel)
1951           t = img.params_sel{i};
1952           ti = find(t == e);
1953           t(ti) = []; % rm 'e' from the list of params_sel
1954           ti = find(t > e);
1955           t(ti) = t(ti)-1; % down-count element indices greater than our deleted one
1956           img.params_sel{i} = t;
1957        end
1958     end
1959 
1960     % show what we got for a background value
1961     if(opt.verbose > 1)
1962        bg = img.elem_data_background;
1963        bg = map_data(bg, in, 'resistivity');
1964        fprintf('%s  background conductivity: %0.1f Ohm.m\n', indent, bg);
1965     end
1966 
1967 function b = has_params(s)
1968 b = false;
1969 if isstruct(s)
1970    b = any(ismember(fieldnames(s),supported_params));
1971 end
1972 
1973 % wrapper function for to_base_types
1974 function out = map_img_base_types(img)
1975   out = to_base_types(img.current_params);
1976 
1977 % convert from know types to their base types
1978 % A helper function for getting to a basic paramterization
1979 % prior to any required scaling, etc.
1980 function type = to_base_types(type)
1981   if ~iscell(type)
1982      type = {type};
1983   end
1984   for i = 1:length(type);
1985      type(i) = {strrep(type{i}, 'log_', '')};
1986      type(i) = {strrep(type{i}, 'log10_', '')};
1987      type(i) = {strrep(type{i}, 'resistivity', 'conductivity')};
1988      type(i) = {strrep(type{i}, 'apparent_resistivity', 'voltage')};
1989   end
1990 
1991 function img = map_img(img, out);
1992    err_if_inf_or_nan(img.elem_data, 'img-pre');
1993    try 
1994        in = img.current_params;
1995    catch
1996        in = {'conductivity'};
1997    end
1998    % make cell array of strings
1999    if ischar(in)
2000       in = {in};
2001       img.current_params = in;
2002    end
2003    if ischar(out)
2004       out = {out};
2005    end
2006 
2007    % if we have mixed data, check that we have a selector to differentiate between them
2008    if ~isfield(img, 'params_sel')
2009       if length(in(:)) == 1
2010          img.params_sel = {1:size(img.elem_data,1)};
2011       else
2012          error('found multiple parametrizations (params) but no params_sel cell array in img');
2013       end
2014    end
2015 
2016    % create data?! we don't know how
2017    if length(out(:)) > length(in(:))
2018       error('missing data (more out types than in types)');
2019    elseif length(out(:)) < length(in(:))
2020       % delete data: we can do that
2021       % TODO we could genearlize this into a reorganizing tool BUT we're just
2022       % interested in something that works, so if we have more than one out(:),
2023       % we don't know what to do currently and error out
2024       if length(out(:)) ~= 1
2025          error('map_img can convert ALL parameters or select a SINGLE output type from multiple input types');
2026       end
2027       inm  = to_base_types(in);
2028       outm = to_base_types(out);
2029       del = sort(find(~strcmp(outm(1), inm(:))), 'descend'); % we do this loop backwards in the hopes of avoiding shuffling data that is about to be deleted
2030       if length(del)+1 ~= length(in)
2031          error('Confused about what to remove from the img. You''ll need to sort the parametrizations out yourself when removing data.');
2032       end
2033       for i = del(:)' % delete each of the extra indices
2034          ndel = length(img.params_sel{i}); % number of deleted elements
2035          for j = i+1:length(img.params_sel)
2036             img.params_sel{j} = img.params_sel{j} - ndel;
2037          end
2038          img.elem_data(img.params_sel{i}) = []; % rm elem_data
2039          img.params_sel(i) = []; % rm params_sel
2040          img.current_params(i) = []; % rm current_params
2041       end
2042       in = img.current_params;
2043    end
2044 
2045    % the sizes now match, we can do the mapping
2046    for i = 1:length(out(:))
2047       % map the data
2048       x = img.elem_data(img.params_sel{i});
2049       x = map_data(x, in{i}, out{i});
2050       img.elem_data(img.params_sel{i}) = x;
2051       img.current_params{i} = out{i};
2052    end
2053    err_if_inf_or_nan(img.elem_data, 'img-post');
2054 
2055    % clean up params_sel/current_params if we only have one parametrization
2056    if length(img.current_params(:)) == 1
2057       img.current_params = img.current_params{1};
2058       img = rmfield(img, 'params_sel'); % unnecessary since we know its all elem_data
2059    end
2060 
2061 function x = map_data(x, in, out)
2062    % check that in and out are single strings, not lists of strings
2063    if ~ischar(in)
2064       if iscell(in) && (length(in(:)) == 1)
2065          in = in{1};
2066       else
2067          error('expecting single string for map_data() "in" type');
2068       end
2069    end
2070    if ~ischar(out)
2071       if iscell(out) && (length(out(:)) == 1)
2072          out = out{1};
2073       else
2074          error('expecting single string for map_data() "out" type');
2075       end
2076    end
2077 
2078    % quit early if there is nothing to do
2079    if strcmp(in, out) % in == out
2080       return; % do nothing
2081    end
2082 
2083    % resistivity to conductivity conversion
2084    % we can't get here if in == out
2085    % we've already checked for log convserions on input or output
2086    if any(strcmp(in,  {'resistivity', 'conductivity'})) && ...
2087       any(strcmp(out, {'resistivity', 'conductivity'}))
2088       x = 1./x; % conductivity <-> resistivity
2089    % log conversion
2090    elseif any(strcmp({in(1:3), out(1:3)}, 'log'))
2091       % log_10 x -> x
2092       if strcmp(in(1:6), 'log10_')
2093          if any(x >= log10(realmax)-eps) warning('loss of precision -> inf'); end
2094          x = map_data(10.^x, in(7:end), out);
2095       % ln x -> x
2096       elseif strcmp(in(1:4), 'log_')
2097          if any(x >= log(realmax)-eps) warning('loss of precision -> inf'); end
2098          x = map_data(exp(x), in(5:end), out);
2099       % x -> log_10 x
2100       elseif strcmp(out(1:6), 'log10_')
2101          if any(x <= 0 + eps) warning('loss of precision -> -inf'); end
2102          x = log10(map_data(x, in, out(7:end)));
2103       % x -> ln x
2104       elseif strcmp(out(1:4), 'log_')
2105          if any(x <= 0 + eps) warning('loss of precision -> -inf'); end
2106          x = log(map_data(x, in, out(5:end)));
2107       else
2108          error(sprintf('unknown conversion (log conversion?) %s - > %s', in, out));
2109       end
2110    else
2111       error('unknown conversion %s -> %s', in, out);
2112    end
2113    x(x == +inf) = +realmax;
2114    x(x == -inf) = -realmax;
2115    err_if_inf_or_nan(x, 'map_data-post');
2116 
2117 function b = map_meas(b, N, in, out)
2118    err_if_inf_or_nan(b, 'map_meas-pre');
2119    if strcmp(in, out) % in == out
2120       return; % do nothing
2121    end
2122 
2123    % resistivity to conductivity conversion
2124    % we can't get here if in == out
2125    if     strcmp(in, 'voltage') && strcmp(out, 'apparent_resistivity')
2126       if N == 1
2127          error('missing apparent resistivity conversion factor N');
2128       end
2129       b = N * b; % voltage -> apparent resistivity
2130    elseif strcmp(in, 'apparent_resistivity') && strcmp(out, 'voltage')
2131       if N == 1
2132          error('missing apparent resistivity conversion factor N');
2133       end
2134       b = N \ b; % apparent resistivity -> voltage
2135    % log conversion
2136    elseif any(strcmp({in(1:3), out(1:3)}, 'log'))
2137       % log_10 b -> b
2138       if strcmp(in(1:6), 'log10_')
2139          if any(b > log10(realmax)-eps) warning('loss of precision -> inf'); end
2140          b = map_meas(10.^b, N, in(7:end), out);
2141       % ln b -> b
2142       elseif strcmp(in(1:4), 'log_')
2143          if any(b > log(realmax)-eps) warning('loss of precision -> inf'); end
2144          b = map_meas(exp(b), N, in(5:end), out);
2145       % b -> log_10 b
2146       elseif strcmp(out(1:6), 'log10_')
2147          if any(b <= 0+eps) warning('loss of precision -> -inf'); end
2148          b = log10(map_meas(b, N, in, out(7:end)));
2149       % b -> ln b
2150       elseif strcmp(out(1:4), 'log_')
2151          if any(b <= 0+eps) warning('loss of precision -> -inf'); end
2152          b = log(map_meas(b, N, in, out(5:end)));
2153       else
2154          error(sprintf('unknown conversion (log conversion?) %s - > %s', in, out));
2155       end
2156    else
2157       error('unknown conversion %s -> %s', in, out);
2158    end
2159    err_if_inf_or_nan(b, 'map_meas-post');
2160 
2161 function x=range(y)
2162 x=max(y)-min(y);
2163 
2164 function pass=do_unit_test(solver)
2165    pass=1;
2166    if nargin == 0
2167       solver = 'inv_solve_core';
2168       test = 'all'
2169    elseif ((length(solver) < 9) || ~strcmp(solver(1:9),'inv_solve'))
2170       test = solver;
2171       solver = 'inv_solve_core';
2172    else
2173       test = 'all'
2174    end
2175    switch(test)
2176       case 'all'
2177          do_unit_test_sub;
2178          do_unit_test_diff;
2179          do_unit_test_rec1(solver);
2180          do_unit_test_rec2(solver);
2181          do_unit_test_rec_mv(solver);
2182       case 'sub'
2183          do_unit_test_sub;
2184       case 'diff'
2185          do_unit_test_diff;
2186       case 'rec1'
2187          do_unit_test_rec1(solver);
2188       case 'rec2'
2189          do_unit_test_rec2(solver);
2190       case 'rec_mv'
2191          do_unit_test_rec_mv(solver);
2192       otherwise
2193          error(['unrecognized solver or tests: ' test]);
2194    end
2195 %pass = pass & do_unit_test_rec2(solver);
2196 % TODO the ..._rec2 unit test is very, very slow... what can we do to speed it up... looks like the perturbations get kinda borked when using the line_search_onm2
2197 
2198 function [imdl, vh, imgi, vi] = unit_test_imdl()
2199    imdl = mk_geophysics_model('h22c',32);
2200    imdl.solve = 'inv_solve_core';
2201    imdl.inv_solve_core.max_iterations = 1; % see that we get the same as the GN 1-step difference soln
2202    imdl.inv_solve_core.verbose = 0;
2203    imdl.reconst_type = 'difference';
2204    imgh = mk_image(imdl,1);
2205    vh = fwd_solve(imgh);
2206 
2207    if nargout > 2
2208       imgi = imgh;
2209       ctrs = interp_mesh(imdl.fwd_model);
2210       x = ctrs(:,1);
2211       y = ctrs(:,2);
2212       r1 = sqrt((x+15).^2 + (y+25).^2);
2213       r2 = sqrt((x-85).^2 + (y+65).^2);
2214       imgi.elem_data(r1<25)= 1/2;
2215       imgi.elem_data(r2<30)= 2;
2216       clf; show_fem(imgi,1);
2217       vi = fwd_solve(imgi);
2218    end
2219 
2220 function do_unit_test_diff()
2221    [imdl, vh, imgi, vi] = unit_test_imdl();
2222 
2223    imdl.reconst_type = 'absolute';
2224    img_abs = inv_solve(imdl,vi);
2225    imdl.reconst_type = 'difference';
2226    img_itr = inv_solve(imdl,vh,vi);
2227    imdl.inv_solve_core.max_iterations = 1; % see that we get the same as the GN 1-step difference soln
2228    imdl.inv_solve_core.line_search_func = @ret1_func;
2229    img_it1 = inv_solve(imdl,vh,vi);
2230    imdl.solve = 'inv_solve_diff_GN_one_step';
2231    img_gn1 = inv_solve(imdl,vh,vi);
2232    clf; subplot(223); show_fem(img_it1,1); title('GN 1-iter');
2233         subplot(224); show_fem(img_gn1,1); title('GN 1-step');
2234         subplot(221); h=show_fem(imgi,1);  title('fwd model'); set(h,'EdgeColor','none');
2235         subplot(222); show_fem(img_itr,1); title('GN 10-iter');
2236         subplot(222); show_fem(img_abs,1); title('GN abs 10-iter');
2237    unit_test_cmp('core (1-step) vs. diff_GN_one_step', img_it1.elem_data, img_gn1.elem_data, eps*1e3);
2238    unit_test_cmp('core (1-step) vs. core (N-step)   ', img_it1.elem_data, img_itr.elem_data, eps);
2239    unit_test_cmp('core (N-step) vs. abs  (N-step)   ', img_it1.elem_data, img_abs.elem_data-1, eps);
2240 
2241 % test sub-functions
2242 % map_meas, map_data
2243 % jacobian scalings
2244 function do_unit_test_sub
2245 d = 1;
2246 while d ~= 1 & d ~= 0
2247   d = rand(1);
2248 end
2249 disp('TEST: map_data()');
2250 elem_types = {'conductivity', 'log_conductivity', 'log10_conductivity', ...
2251               'resistivity',  'log_resistivity',  'log10_resistivity'};
2252 expected = [d         log(d)         log10(d)      1./d      log(1./d)      log10(1./d); ...
2253             exp(d)    d              log10(exp(d)) 1./exp(d) log(1./exp(d)) log10(1./exp(d)); ...
2254             10.^d     log(10.^d )    d             1./10.^d  log(1./10.^d ) log10(1./10.^d ); ...
2255             1./d      log(1./d  )    log10(1./d)   d         log(d)         log10(d); ...
2256             1./exp(d) log(1./exp(d)) log10(1./exp(d)) exp(d) d              log10(exp(d)); ...
2257             1./10.^d  log(1./10.^d)  log10(1./10.^d)  10.^d  log(10.^d)     d ];
2258 for i = 1:length(elem_types)
2259   for j = 1:length(elem_types)
2260     test_map_data(d, elem_types{i}, elem_types{j}, expected(i,j));
2261   end
2262 end
2263 
2264 disp('TEST: map_meas()');
2265 N = 1/15;
2266 Ninv = 1/N;
2267 % function b = map_meas(b, N, in, out)
2268 elem_types = {'voltage', 'log_voltage', 'log10_voltage', ...
2269               'apparent_resistivity',  'log_apparent_resistivity',  'log10_apparent_resistivity'};
2270 expected = [d         log(d)         log10(d)      N*d      log(N*d)      log10(N*d); ...
2271             exp(d)    d              log10(exp(d)) N*exp(d) log(N*exp(d)) log10(N*exp(d)); ...
2272             10.^d     log(10.^d )    d             N*10.^d  log(N*10.^d ) log10(N*10.^d ); ...
2273             Ninv*d      log(Ninv*d  )    log10(Ninv*d)   d         log(d)         log10(d); ...
2274             Ninv*exp(d) log(Ninv*exp(d)) log10(Ninv*exp(d)) exp(d) d              log10(exp(d)); ...
2275             Ninv*10.^d  log(Ninv*10.^d)  log10(Ninv*10.^d)  10.^d  log(10.^d)     d ];
2276 for i = 1:length(elem_types)
2277   for j = 1:length(elem_types)
2278     test_map_meas(d, N, elem_types{i}, elem_types{j}, expected(i,j));
2279   end
2280 end
2281 
2282 disp('TEST: Jacobian scaling');
2283 d = [d d]';
2284 unit_test_cmp( ...
2285    sprintf('Jacobian scaling (%s)', 'conductivity'), ...
2286    ret1_func(d), 1);
2287 
2288 unit_test_cmp( ...
2289    sprintf('Jacobian scaling (%s)', 'log_conductivity'), ...
2290    dx_dlogx(d), diag(d));
2291 
2292 unit_test_cmp( ...
2293    sprintf('Jacobian scaling (%s)', 'log10_conductivity'), ...
2294    dx_dlog10x(d), diag(d)*log(10));
2295 
2296 unit_test_cmp( ...
2297    sprintf('Jacobian scaling (%s)', 'resistivity'), ...
2298    dx_dy(d), diag(-d.^2));
2299 
2300 unit_test_cmp( ...
2301    sprintf('Jacobian scaling (%s)', 'log_resistivity'), ...
2302    dx_dlogy(d), diag(-d));
2303 
2304 unit_test_cmp( ...
2305    sprintf('Jacobian scaling (%s)', 'log10_resistivity'), ...
2306    dx_dlog10y(d), diag(-d)/log(10));
2307 
2308 
2309 function test_map_data(data, in, out, expected)
2310 %fprintf('TEST: map_data(%s -> %s)\n', in, out);
2311    calc_val = map_data(data, in, out);
2312    str = sprintf('map_data(%s -> %s)', in, out);
2313    unit_test_cmp(str, calc_val, expected)
2314 
2315 function test_map_meas(data, N, in, out, expected)
2316 %fprintf('TEST: map_meas(%s -> %s)\n', in, out);
2317    calc_val = map_meas(data, N, in, out);
2318    str = sprintf('map_data(%s -> %s)', in, out);
2319    unit_test_cmp(str, calc_val, expected)
2320 
2321 
2322 % a couple easy reconstructions
2323 % check c2f, apparent_resistivity, log_conductivity, verbosity don't error out
2324 function do_unit_test_rec1(solver)
2325 % -------------
2326 % ADAPTED FROM
2327 % Create simulation data $Id: inv_solve_core.m 5563 2017-06-19 02:33:27Z alistair_boyle $
2328 %  http://eidors3d.sourceforge.net/tutorial/adv_image_reconst/basic_iterative.shtml
2329 % 3D Model
2330 [imdl, vh, imgi, vi] = unit_test_imdl();
2331 imdl.solve = solver;
2332 imdl.reconst_type = 'absolute';
2333 imdl.inv_solve_core.prior_data = 1;
2334 imdl.inv_solve_core.elem_prior = 'conductivity';
2335 imdl.inv_solve_core.elem_working = 'log_conductivity';
2336 imdl.inv_solve_core.meas_working = 'apparent_resistivity';
2337 imdl.inv_solve_core.calc_solution_error = 0;
2338 imdl.inv_solve_core.verbose = 0;
2339 
2340 imgp = map_img(imgi, 'log10_conductivity');
2341 % add noise
2342 %Add 30dB SNR noise to data
2343 %noise_level= std(vi.meas - vh.meas)/10^(30/20);
2344 %vi.meas = vi.meas + noise_level*randn(size(vi.meas));
2345 % Reconstruct Images
2346 img1= inv_solve(imdl, vi);
2347 % -------------
2348 disp('TEST: previous solved at default verbosity');
2349 disp('TEST: now solve same at verbosity=0 --> should be silent');
2350 imdl.inv_solve_core.verbose = 0;
2351 %imdl.inv_solve_core.meas_working = 'apparent_resistivity';
2352 img2= inv_solve(imdl, vi);
2353 % -------------
2354 imdl.inv_solve_core.elem_output = 'log10_resistivity'; % resistivity output works
2355 img3= inv_solve(imdl, vi);
2356 
2357 disp('TEST: try coarse2fine mapping with background');
2358 ctrs= interp_mesh(imdl.rec_model);
2359 x= ctrs(:,1); y= ctrs(:,2);
2360 r=sqrt((x+5).^2 + (y+5).^2);
2361 imdl.rec_model.elems(x-y > 300, :) = [];
2362 c2f = mk_approx_c2f(imdl.fwd_model,imdl.rec_model);
2363 imdl.fwd_model.coarse2fine = c2f;
2364 % solve
2365 %imdl.inv_solve_core.verbose = 10;
2366 imdl.inv_solve_core.elem_output = 'conductivity';
2367 img4= inv_solve(imdl, vi);
2368 
2369 clf; subplot(221); h=show_fem(imgp,1); set(h,'EdgeColor','none'); axis tight; title('synthetic data, logC');
2370    subplot(222); show_fem(img1,1); axis tight; title('#1 verbosity=default');
2371    subplot(223); show_fem(img2,1); axis tight; title('#2 verbosity=0');
2372    subplot(223); show_fem(img3,1); axis tight; title('#3 c2f + log10 resistivity out');
2373    subplot(224); show_fem(img4,1); axis tight; title('#4 c2f cut');
2374    drawnow;
2375 d1 = img1.elem_data; d2 = img2.elem_data; d3 = img3.elem_data; d4 = img4.elem_data;
2376 unit_test_cmp('img1 == img2', d1, d2, eps);
2377 unit_test_cmp('img1 == img3', d1, 1./(10.^d3), eps);
2378 unit_test_cmp('img1 == img4', (d1(x-y <=300)-d4)./d4, 0, 0.05); %  +-5% error
2379 
2380 %clf; subplot(211); plot((img3.elem_data - img1.elem_data)./img1.elem_data);
2381 
2382 % a couple easy reconstructions with movement or similar
2383 function do_unit_test_rec_mv(solver)
2384 disp('TEST: conductivity and movement --> baseline conductivity only');
2385 % -------------
2386 % ADAPTED FROM
2387 % Create simulation data $Id: inv_solve_core.m 5563 2017-06-19 02:33:27Z alistair_boyle $
2388 %  http://eidors3d.sourceforge.net/tutorial/adv_image_reconst/basic_iterative.shtml
2389 % 3D Model
2390 imdl= mk_common_model('c2t4',16); % 576 elements
2391 ne = length(imdl.fwd_model.electrode);
2392 nt = length(imdl.fwd_model.elems);
2393 imdl.solve = solver;
2394 imdl.reconst_type = 'absolute';
2395 % specify the units to work in
2396 imdl.inv_solve_core.meas_input   = 'voltage';
2397 imdl.inv_solve_core.meas_working = 'apparent_resistivity';
2398 imdl.inv_solve_core.elem_prior   = {   'conductivity'   };
2399 imdl.inv_solve_core.prior_data   = {        1           };
2400 imdl.inv_solve_core.elem_working = {'log_conductivity'};
2401 imdl.inv_solve_core.elem_output  = {'log10_conductivity'};
2402 imdl.inv_solve_core.calc_solution_error = 0;
2403 imdl.inv_solve_core.verbose = 0;
2404 imdl.hyperparameter.value = 0.1;
2405 
2406 % set homogeneous conductivity and simulate
2407 imgsrc= mk_image( imdl.fwd_model, 1);
2408 vh=fwd_solve(imgsrc);
2409 % set inhomogeneous conductivity
2410 ctrs= interp_mesh(imdl.fwd_model);
2411 x= ctrs(:,1); y= ctrs(:,2);
2412 r1=sqrt((x+5).^2 + (y+5).^2); r2 = sqrt((x-45).^2 + (y-35).^2);
2413 imgsrc.elem_data(r1<50)= 0.05;
2414 imgsrc.elem_data(r2<30)= 100;
2415 
2416 % inhomogeneous data
2417 vi=fwd_solve( imgsrc );
2418 % add noise
2419 %Add 30dB SNR noise to data
2420 noise_level= std(vi.meas - vh.meas)/10^(30/20);
2421 vi.meas = vi.meas + noise_level*randn(size(vi.meas));
2422 
2423 % show model
2424 hh=clf; subplot(221); imgp = map_img(imgsrc, 'log10_conductivity'); show_fem(imgp,1); axis tight; title('synth baseline, logC');
2425 
2426 % Reconstruct Images
2427 img0= inv_solve(imdl, vi);
2428 figure(hh); subplot(222);
2429  img0 = map_img(img0, 'log10_conductivity');
2430  show_fem(img0, 1); axis tight;
2431 
2432 disp('TEST: conductivity + movement');
2433 imdl.fwd_model = rmfield(imdl.fwd_model, 'jacobian');
2434 % specify the units to work in
2435 imdl.inv_solve_core.elem_prior   = {   'conductivity'   , 'movement'};
2436 imdl.inv_solve_core.prior_data   = {        1           ,     0     };
2437 imdl.inv_solve_core.RtR_prior    = {     @eidors_default, @prior_movement_only};
2438 imdl.inv_solve_core.elem_len     = {       nt           ,   ne*2    };
2439 imdl.inv_solve_core.elem_working = {  'log_conductivity', 'movement'};
2440 imdl.inv_solve_core.elem_output  = {'log10_conductivity', 'movement'};
2441 imdl.inv_solve_core.jacobian     = { @jacobian_adjoint  , @jacobian_movement_only};
2442 imdl.inv_solve_core.hyperparameter = {   [1 1.1 0.9]    ,  sqrt(2e-3)     }; % multiplied by imdl.hyperparameter.value
2443 imdl.inv_solve_core.verbose = 0;
2444 
2445 % electrode positions before
2446 nn = [imgsrc.fwd_model.electrode(:).nodes];
2447 elec_orig = imgsrc.fwd_model.nodes(nn,:);
2448 % set 2% radial movement
2449 nn = imgsrc.fwd_model.nodes;
2450 imgsrc.fwd_model.nodes = nn * [1-0.02 0; 0 1+0.02]; % 1% compress X, 1% expansion Y, NOT conformal
2451 % electrode positions after
2452 nn = [imgsrc.fwd_model.electrode(:).nodes];
2453 elec_mv = imgsrc.fwd_model.nodes(nn,:);
2454 
2455 % inhomogeneous data
2456 vi=fwd_solve( imgsrc );
2457 % add noise
2458 %Add 30dB SNR noise to data
2459 noise_level= std(vi.meas - vh.meas)/10^(30/20);
2460 %vi.meas = vi.meas + noise_level*randn(size(vi.meas));
2461 
2462 % show model
2463 nn = [imgsrc.fwd_model.electrode(1:4).nodes];
2464 figure(hh); subplot(223); imgp = map_img(imgsrc, 'log10_conductivity'); show_fem_move(imgp,elec_mv-elec_orig,10,1); axis tight; title('synth mvmt, logC');
2465 
2466 % Reconstruct Images
2467 img1= inv_solve(imdl, vi);
2468 figure(hh); subplot(224);
2469  imgm = map_img(img1, 'movement');
2470  img1 = map_img(img1, 'log10_conductivity');
2471  show_fem_move(img1,reshape(imgm.elem_data,16,2), 10, 1); axis tight;
2472 
2473 d0 = img0.elem_data; d1 = img1.elem_data;
2474 unit_test_cmp('img0 == img1 + mvmt', d0-d1, 0, 0.40);
2475 
2476 % helper function: calculate jacobian movement by itself
2477 function Jm = jacobian_movement_only (fwd_model, img);
2478   pp = fwd_model_parameters(img.fwd_model);
2479   szJm = pp.n_elec * pp.n_dims; % number of electrodes * dimensions
2480   img = map_img(img, 'conductivity'); % expect conductivity only
2481   Jcm = jacobian_movement(fwd_model, img);
2482   Jm = Jcm(:,(end-szJm+1):end);
2483 %% this plot shows we are grabing the right section of the Jacobian
2484 %  figure();
2485 %  subplot(311); imagesc(Jcm); axis ij equal tight; xlabel(sprintf('||Jcm||=%g',norm(Jcm))); colorbar;
2486 %  Jc = jacobian_adjoint(fwd_model, img);
2487 %  subplot(312); imagesc([Jc Jm]); axis ij equal tight; xlabel(sprintf('||[Jc Jm]||=%g',norm([Jc Jm]))); colorbar;
2488 %  dd = abs([Jc Jm]-Jcm); % difference
2489 %  subplot(313); imagesc(dd); axis ij equal tight; xlabel(sprintf('|| |[Jc Jm]-Jcm| ||=%g',norm(dd))); colorbar;
2490 
2491 function RtR = prior_movement_only(imdl);
2492   imdl.image_prior.parameters(1) = 1; % weighting of movement vs. conductivity ... but we're dropping conductivity here
2493   pp = fwd_model_parameters(imdl.fwd_model);
2494   szPm = pp.n_elec * pp.n_dims; % number of electrodes * dimensions
2495   RtR = prior_movement(imdl);
2496   RtR = RtR((end-szPm+1):end,(end-szPm+1):end);
2497 
2498 function do_unit_test_rec2(solver)
2499 disp('TEST: reconstruct a discontinuity');
2500 [imdl, vh] = unit_test_imdl();
2501 imgi = mk_image(imdl, 1);
2502 ctrs = interp_mesh(imdl.fwd_model);
2503 x = ctrs(:,1); y = ctrs(:,2);
2504 imgi.elem_data(x-y>100)= 1/10;
2505 vi = fwd_solve(imgi); vi = vi.meas;
2506 clf; show_fem(imgi,1);
2507 
2508 imdl.solve = solver;
2509 imdl.hyperparameter.value = 1e2; % was 0.1
2510 
2511 imdl.reconst_type = 'absolute';
2512 imdl.inv_solve_core.elem_working = 'log_conductivity';
2513 imdl.inv_solve_core.elem_output = 'log10_resistivity';
2514 imdl.inv_solve_core.meas_working = 'apparent_resistivity';
2515 imdl.inv_solve_core.dtol_iter = 4; % default 1 -> start checking on the first iter
2516 imdl.inv_solve_core.max_iterations = 20; % default 10
2517 imdl.inv_solve_core.calc_solution_error = 0;
2518 imdl.inv_solve_core.verbose = 10;
2519 imgr= inv_solve_core(imdl, vi);
2520 
2521 clf; show_fem(imgr,1);
2522 img = mk_image( imdl );
2523 img.elem_data= 1./(10.^imgr.elem_data);
2524 
2525 %I = 1; % TODO FIXME -> I is diag(1./vh) the conversion to apparent resistivity
2526 %% TODO these plots are useful, get them built into the solver!
2527 %vCG= fwd_solve(img); vCG = vCG.meas; fmdl = imdl.fwd_model;
2528 %clf; plot(I*(vi-vCG)); title('data misfit');
2529 %clf; hist(abs(I*(vi-vCG)),50); title('|data misfit|, histogram'); xlabel('|misfit|'); ylabel('count');
2530 %clf; show_pseudosection( fmdl, I*vi); title('measurement data');
2531 %clf; show_pseudosection( fmdl, I*vCG); title('reconstruction data');
2532 %clf; show_pseudosection( fmdl, (vCG-vi)/vi*100); title('data misfit');

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