inv_solve_core

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 1)
      0: quiet
    >=1: print iteration count and residual
    >=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 cholesky)
     Method to use for solving dx: 'pcg' or 'cholesky'
     If 'cholesky', will fall back to 'pcg' if
     out-of-memory. Matlab has trouble detecting the
     out-of-memory condition on many machines, and is
     likely to just grind to a halt swapping. Beware.
   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, hps2RtR, hpt2LLt)
   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
    hps2  - spatial hyperparameter squared, see CALC_HYPERPARAMETER
    RtR   - regularization matrix squared --> hps2RtR = hps2*RtR
    hpt2  - temporal hyperparameter squared
    LLt   - temporal regularization matrix squared --> hpt2LLt = hpt2*LLt

   Signature for line_optimize_func
    [alpha, img, dv, opt] = f(img, sx, data0, img0, N, W, hps2RtR, hpt2LLt, 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
    hps2  - spatial hyperparameter squared, see CALC_HYPERPARAMETER
    RtR   - regularization matrix squared --> hps2RtR = hps2*RtR
    hpt2  - temporal hyperparameter squared
    LLt   - temporal regularization matrix squared --> hpt2LLt = hpt2*LLt
    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

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