inv_solve_core

PURPOSE

INV_SOLVE_CORE Solver using a generic iterative algorithm

SYNOPSIS

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

DESCRIPTION

INV_SOLVE_CORE Solver using a generic iterative algorithm
img => output image (or vector of images)
inv_model => inverse model struct
data0 => EIT data
data0, data1 => difference EIT data
This function is parametrized and uses function pointers where possible to
allow its use as a general iterative solver framework. There are a large
number of parameters and functions contained here. Sensible defaults are
used throughout. You do not need to set every parameter.
The solver operates as an absolute Gauss-Newton iterative solver by default.
Wrapper functions are available to call this function in its various forms.
Look forward to the "See also" section at the end of this help.
Argument matrices to the internal functions (measurement inverse covariance,
for example) are only calculated if required. Functions that are supplied to
this INV_SOLVE_CORE must be able to survive being probed: they will have
each parameter set to either 0 or 1 to determine if the function is sensitive
to that argument. This works cleanly for most matrix multiplication based
functions but for more abstract code, some handling of this behaviour may
need to be implemented.
In the following parameters, r_k is the current residual, r_{k-1} is the
previous iteration's residual. k is the iteration count.
Parameters denoted with a ** to the right of their default values are
deprecated legacy parameters, some of which formerly existed under
'inv_model.parameters.*'.
Parameters (inv_model.inv_solve_core.*):
verbose (show progress) (default 2)
0: quiet
>=1: print iteration count
>=2: print details as the algorithm progresses
>=3: plot residuals versus iteration count
>=4: plot result at each iteration, see show_fem
>=5: plot line search per iteration
plot_residuals (default 0)
plot residuals without verbose output
fig_prefix (default: <none>)
figure file prefix; figures not saved if <none>
fwd_solutions (default 0)
0: ignore
1: count fwd_solve(), generally the most
computationally expensive component of
the iterations
residual_func = (default @GN_residual)
NOTE: @meas_residual exists to maintain
compatibility with some older code
max_iterations (default 10) **
ntol (estimate of machine precision) (default eps)
tol (stop iter if r_k < tol) (default 0)
dtol (default -0.01%)
stop iter if (r_k - r_{k-1})/r_1 < dtol AND
k >= dtol_iter
dtol_iter (default 0)
apply dtol stopping criteria if k >= dtol_iter
min_value (default -inf) **
max_value (default +inf) **
line_optimize_func (default <none>) ** TODO
[next,fmin,res]=f(org,dx,data0,opt);
opt=line_optimize.* + objective_func
Deprecated, use line_search_func instead.
line_optimize.perturb
(default line_search_args.perturb) ** TODO
Deprecated, use line_search_args.perturb instead.
update_func (default TODO) ** TODO
[img,opt]=f(org,next,dx,fmin,res,opt)
Deprecated, use <TODO> instead.
update_method (default lu)
Method to use for solving dx: 'pcg' or 'lu'
If 'lu', will fall back to 'pcg' if out-of-memory.
do_starting_estimate (default 1) ** TODO
Deprecated, use <TODO> instead.
line_search_func (default @line_search_onm2)
line_search_dv_func (default @update_dv_core)
line_search_de_func (default @update_de_core)
line_search_args.perturb
(default [0 1/16 1/8 1/4 1/2 1])
line search for alpha by these steps along sx
line_search_args.plot (default 0)
c2f_background (default 0)
if > 0, this is additional elem_data
if a c2f map exists, the default is to decide
based on an estimate of c2f overlap whether a
background value is required. If a background is
required, it is added as the last element of that
type.
c2f_background_fixed (default 1)
hold the background estimate fixed or allow it
to vary as any other elem_data
elem_fixed (default [])
meas_select already handles selecting from the
valid measurements. we want the same for the
elem_data, so we only work on modifying the
legal values.
Note that c2f_background's elements are added to
this list if c2f_background_fixed == 1.
prior_data (default to jacobian_bkgnd)
Sets the priors of type elem_prior. May be
scalar, per elem_prior, or match the working
length of each elem_data type. Note that for priors
using the c2f a background element may be added
to the end of that range when required; see
c2f_background.
elem_len (default to all elem_data)
A cell array list of how many of each
elem_working there are in elem_data.
prior_data = { 32.1, 10*ones(10,1) };
elem_prior = {'conductivity', 'movement'};
elem_len = { 20001, 10 };
elem_prior (default 'conductivity')
Input 'prior_data' type; immediately converted to
'elem_working' type before first iteration.
elem_working (default 'conductivity')
elem_output (default 'conductivity')
The working and output units for 'elem_data'.
Valid types are 'conductivity' and 'resistivity'
as plain units or with the prefix 'log_' or
'log10_'. Conversions are handled internally.
Scaling factors are applied to the Jacobian
(calculated in units of 'conductivity') as
appropriate.
If elem_working == elem_output, then no
conversions take place.
For multiple types, use cell array.
ex: elem_output = {'log_resistivity', 'movement'}
meas_input (default 'voltage')
meas_working (default 'voltage')
Similarly to elem_working/output, conversion
between 'voltage' and 'apparent_resistivity' and
their log/log10 variants are handled internally.
If meas_input == meas_working no conversions take
place. The normalization factor 'N' is calculated
if 'apparent_resistivity' is used.
update_img_func (default: pass-through)
Called prior to calc_jacobian and update_dv.
Elements are converted to their "base types"
before this function is called. For example,
'log_resistivity' becomes 'conductivity'.
It is a hook to allow additional updates to the
model before the Jacobian, or a new set of
measurements are calculated via fwd_solve.
return_working_variables (default: 0)
If 1, return the last working variables to the user
img.inv_solve_{type}.J Jacobian
img.inv_solve_{type}.dx descent direction
img.inv_solve_{type}.sx search direction
img.inv_solve_{type}.alpha line search result
img.inv_solve_{type}.beta conjugation parameter
img.inv_solve_{type}.r as:
[ residual, measurement misfit, element misfit ]
with one row per iteration
show_fem (default: @show_fem)
Function with which to plot each iteration's
current parameters.
Signature for residual_func
[r,m,e] = f(dv, de, W, hp2RtR)
where
r - the residual = m + e
m - measurement error
e - prior misfit
dv - change in voltage
de - change in image elements
W - measurement inverse covariance matrix
hp2 - hyperparameter squared, see CALC_HYPERPARAMETER
RtR - regularization matrix squared --> hp2RtR = hp2*RtR
Signature for line_optimize_func
[alpha, img, dv, opt] = f(img, sx, data0, img0, N, W, hp2RtR, dv, opt)
where:
alpha - line search result
img - the current image
(optional, recalculated if not available)
sx - the search direction to which alpha should be applied
data0 - the true measurements (dv = N*data - N*data0)
img0 - the image background (de = img - img0)
N - a measurement normalization factor, N*dv
W - measurement inverse covariance matrix
hp2 - hyperparameter squared, see CALC_HYPERPARAMETER
RtR - regularization matrix squared --> hp2RtR = hp2*RtR
dv - change in voltage
(optional, recalculated if not available)
opt - additional arguments, updated at each call
Signature for line_search_dv_func
[dv, opt] = update_dv_core(img, data0, N, opt)
where:
dv - change in voltage
opt - additional arguments, updated at each call
data - the estimated measurements
img - the current image
data0 - the true measurements
N - a measurement normalization factor, N*dv
Signature for line_search_de_func
de = f(img, img0, opt)
where:
de - change in image elements
img - the current image
img0 - the image background (de = img - img0)
opt - additional arguments
Signature for calc_jacobian_scaling_func
S = f(x)
where:
S - to be used to scale the Jacobian
x - current img.elem_data in units of 'conductivity'
Signature for update_img_func
img2 = f(img1, opt)
where
img1 - an input image, the current working image
img2 - a (potentially) modified version to be used
for the fwd_solve/Jacobian calculations
NOTE that the default line search is very crude. For
my test problems it seems to amount to an expensive grid
search. Much more efficient line search algorithms exist
and some fragments already are coded elsewhere in the
EIDORS code-base.
See also: INV_SOLVE_ABS_GN, INV_SOLVE_ABS_GN_LOGC,
INV_SOLVE_ABS_CG, INV_SOLVE_ABS_CG_LOGC,
LINE_SEARCH_O2, LINE_SEARCH_ONM2
(C) 2010-2016 Alistair Boyle, Nolwenn Lesparre, Andy Adler, Bartłomiej Grychtol.
License: GPL version 2 or version 3