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