Description of inv_solve_abs

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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
