jInv.InverseSolve
The InverseSolve submodule contains misfit functions, regularizers, optimization and other tools for solving PDE parameter estimation problems.
#
jInv.InverseSolve.GlobalToLocal — Type.
type jInv.InverseSolve.GlobalToLocal
Maps global model to local model
sigLocal = PForInv*sigGlobal + sigmaBackground
Fields: PForInv::SparseMatrixCSC - linear operator between inverse mesh and forward mesh, e.g., linear interpolation matrix and projection on active set. sigmaBackground::Vector - background model
Constructors: getGlobalToLocal(P::SparseMatrixCSC) getGlobalToLocal(P::SparseMatrixCSC,sigmaBack::Vector)
Example: Mesh2Mesh = getInterpolationMatrix(Minv,Mfwd) # Minv and Mfwd have different resolutions sigmaBack = 1.2*ones(Minv.nc) # put background conductivity gloc = getGlobalToLocal(Mesh2Mesh,sigmaBack)
sigLocal = gloc.PForInv' * sigGlobal + gloc.sigmaBack
# this call is equivalent, but needed in case the Mesh2Mesh matrix is compressed
sigLocalFast = interpGlobalToLocal(sigGlobal,gloc.PForInv,gloc.sigmaBack)
#
jInv.InverseSolve.InverseParam — Type.
type jInv.InverseSolve.InverseParam
Type storing parameters for Inversion.
Fields:
Minv::AbstractMesh
modelfun::Function - model function (evaluated by main worker), see models.jl
regularizer::Function - regularizer, see regularizer.jl
alpha::Real - regularization parameter
mref::Array - reference model
boundsLow::Vector - lower bounds for model
boundsHigh::Vector - upper bounds for model
maxStep::Real - maximum step in optimization
pcgMaxIter::Int - maximum number of PCG iterations
pcgTol::Real - tolerance for PCG
minUpdate::Real - stopping criteria
maxIter::Int - maximum number of iterations
HesPrec - A preconditioner for the Hessian.
Constructor:
getInverseParam
Example:
Minv = getRegularMesh(domain,n)
modelfun = expMod
regularizer(m,mref,Minv) = wdiffusionReg(m,mref,Minv)
alpha = 1e-3
mref = zeros(Minv.nc)
pInv = getInverseParam(Minv,modelfun,regularizer,alpha,mref)
#
jInv.InverseSolve.MisfitParam — Type.
type jInv.InverseSolve.MisfitParam
Type storing information about one term in the misfit
F(m) = sum_i^n phi_i(pFor(model(m)),dobs,Wd)
Fields:
pFor::ForwardProbType - forward problem
Wd - inverse standard deviation
dobs - observed data
misfit::Function - misfit function
modelfun::Function - model function (evaluated locally)
gloc - mapping from inverse to forward mesh (including active cells projection).
Constructors:
getMisfitParam(pFor,Wd,dobs,misfit,model,gloc=getGlobalToLocal(1.0))
getMisfitParam(pForRFs::Array{RemoteChannel}, Wd::Array, dobs::Array, misfit::Function, Iact,sigmaBack::Vector, Mesh2MeshRFs::Union{Array{RemoteChannel},Array{Float64}}=ones(length(pForRFs)), modelfun::Function=identityMod,fname="")
#
jInv.InverseSolve.HuberFun — Function.
mis,dmis,d2mis = HuberFun(dc,dobs,Wd,C)
Computes misfit via
misfit(dc,dobs) = sqrt(abs(Wd*res).^2 + eps)
Input:
dc::Array - simulated data
dobs::Array - measured data
Wd::Array - diagional weighting
eps - conditioning parameter (default=1e-3)
Output:
mis::Real - misfit
dmis - gradient
d2mis - diagonal of Hessian
#
jInv.InverseSolve.SSDFun — Method.
For complex data misfit is computed as 0.5*|real(dc)-(dobs)|_Wd^2 + 0.5*|complex(dc)-complex(dobs)|_W^2
#
jInv.InverseSolve.SSDFun — Method.
mis,dmis,d2mis = SSDFun(dc,dobs,Wd)
Input:
dc::Array - simulated data
dobs::Array - measured data
Wd::Array - diagonal weighting
Output:
mis::Real - misfit, 0.5*|dc-dobs|_Wd^2
dmis - gradient
d2mis - diagonal of Hessian
#
jInv.InverseSolve.boundMod — Method.
sigma,dsigma = boundMod(m;boundLow=0.0,boundHigh=1.0)
maps model parameter to conductivity via
sigma = 0.5*(boundHigh-boundLow) * (tanh(m)+1.0 ) + boundLow
#
jInv.InverseSolve.computeMisfit — Function.
Dc,F,dF,d2F = computeMisfit(...)
Computes misfit for PDE parameter estimation problem.
computeMisfit has several options.
#
jInv.InverseSolve.diffusionReg — Method.
Rc,dR,d2R = diffusionReg(m,mref,M,Iact=1.0)
Compute diffusion regularizer
0.5*||GRAD*(m-mref)||_V^2
Input:
m - model
mref - reference model
M - Mesh
Iact - projector on active cells
Output
Rc - value of regularizer
dR - gradient w.r.t. m
d2R - Hessian
#
jInv.InverseSolve.expMod — Method.
sigma,dsigma = expMod(model)
maps model parameter to conductivity via
sigma(m) = exp(m)
#
jInv.InverseSolve.fMod — Method.
sigma,dsigma = fMod(model;f::Function=identity,df::Function=m->speye(length(m)))
maps model parameter to conductivity via
sigma(m) = f(m) and dsigma(m) = sdiag(df(m))
#
jInv.InverseSolve.getInverseParam — Method.
function jInv.InverseSolve.getInverseParam(...)
Constructs an InverseParam
Required Input:
Minv::AbstractMesh - mesh of model
modFun::Function - model
regularizer::Function - regularizer, see regularizer.jl
alpha::Real - regularization parameter
mref - reference model
boundsLow::Vector - lower bounds for model
boundsHigh::Vector - upper bounds for model
Optional Inputs:
maxStep::Real=1.0 - maximum step in optimization
pcgMaxIter::Int=10 - maximum number of PCG iterations
pcgTol::Real - tolerance for PCG
minUpdate::Real=1e-4 - stopping criteria
maxIter::Int=10 - maximum number of iterations
#
jInv.InverseSolve.getMisfitParam — Function.
function jInv.InverseSolve.getMisfitParam
Required Input:
Input for a call: getMisfitParam(pFor,Wd,dobs,misfit,model,gloc=getGlobalToLocal(1.0))
pFor::ForwardProbType - forward problem
Wd - inverse standard deviation of data
dobs - observed data
misfit::Function - Misfit function
modelfun::Function - model function. If all misfits have same model function,
put modelfun in inverseParam and identityMod here for efficiency.
gloc - mapping from inverse to forward mesh (default: identity). See globalToLocal.jl
Input for a call: getMisfitParam(pForRFs::Array{RemoteChannel}, Wd::Array, dobs::Array, misfit::Function, Iact,sigmaBack::Vector, Mesh2MeshRFs::Union{Array{RemoteChannel},Array{Float64}}=ones(length(pForRFs)), modelfun::Function=identityMod,fname="")
pForRFs::Array{RemoteChannel} - Array of remote references for forward problem parameters on each worker.
Wd::Array - Array of inverse standard deviation for the data on each worker.
dobs::Array - Array of observed data for each worker.
misfit::Function - Misfit function
Iact - Projector to active cells.
sigmaBack::Vector - Background model ("frozen" cells).
Mesh2MeshRFs::Union{Array{RemoteChannel},Array{Float64}} - mapping from inverse to forward mesh (default: identity)
modelfun::Function=identityMod - model function. If all misfits have same model function,
put modelfun in inverseParam and identityMod here for efficiency.
fname="" - (optional)
#
jInv.InverseSolve.iteratedTikhonov — Method.
mc,DC = iteratedTikhonov(mc,pInv::InverseParam,pMis,nAlpha,alphaFac,
targetMisfit;indCredit=[],dumpResults::Function=dummy)
Perform (Projected) Gauss-NewtonCG using iterated Tikhonov procedure to decrease
regularization parameter and update reference model after fixed number of GN iterations
set in pInv.
Input:
mc::Vector - Initial guess for model
pInv::InverseParam - parameter for inversion
pMis - misfit terms
nAlpha - maximum number of allowed regularization parameter (alpha) values
alphaFac - alpha decrease factor. alpha_(i+1) = alpha_i/alphaFac
targetMisfit - Termination criterion. iteratedTikhonov will exit
-when data misfit < targetMisfit
indCredit - indices of forward problems to work on
dumpResults - A function pointer for saving the results throughout the iterations.
- We assume that dumpResults is dumpResults(mc,Dc,iter,pInv,pMis),
- where mc is the recovered model, Dc is the predicted data.
- If dumpResults is not given, nothing is done (dummy() is called).
Output:
mc - final model
Dc - final computed data
tikhonovFlag - Data misfit convergence flag
hist - Iteration history. This is a vector. Each entry is
- a structure containing the projGNCG history for each
- alpha value.
#
jInv.InverseSolve.logBarrier — Function.
Rc,dR,d2R = logBarrier(m::Vector, z::Vector, M::AbstractMesh,low::Vector,high::Vector, epsilon)
Computes logBarrier regularizer
R = -log(1 - ((m-high)/epsilon).^2) if high-epsilon < m < high
0 if low+epsilon <= m <= high-epsilon
-log(1 - ((m-low)/epsilon).^2) if low < m < low+epsilon
Input:
m - model
z - not being used. Here for compatibility.
M - Mesh. not being used. Here for compatibility.
low - low bound for each coordinate.
high - high bound for each coordinate.
epsilon - layer width of the barier.
Output
g - value of regularizer
dg - gradient w.r.t. m
d2g - Hessian (diagonal matrix). Second derivative is not continous.
#
jInv.InverseSolve.logBarrierSquared — Function.
Rc,dR,d2R = logBarrierSquared(m::Vector, z::Vector, M::AbstractMesh,low::Vector,high::Vector, epsilon)
Computes logBarrier regularizer
R = (log(1 - ((m-high)/epsilon).^2))^2 if high-epsilon < m < high
0 if low+epsilon <= m <= high-epsilon
(log(1 - ((m-low)/epsilon).^2))^2 if low < m < low+epsilon
Input:
m - model
z - not being used. Here for compatibility.
M - Mesh. not being used. Here for compatibility.
low - low bound for each coordinate.
high - high bound for each coordinate.
epsilon - layer width of the barier.
Output
g - value of regularizer
dg - gradient w.r.t. m
d2g - Gauss Newton Hessian approximation (diagonal matrix). Second derivative approx is continous.
#
jInv.InverseSolve.projGNCG — Method.
mc,Dc,outerFlag = projGNCG(mc,pInv::InverseParam,pMis, indFor = [], dumpResults::Function = dummy)
(Projected) Gauss-NewtonCG
Input:
mc::Vector - intial guess for model
pInv::InverseParam - parameter for inversion
pMis - misfit terms
indCredit - indices of forward problems to work on
dumpResults - A function pointer for saving the results throughout the iterations.
- We assume that dumpResults is dumpResults(mc,Dc,iter,pInv,pMis),
- where mc is the recovered model, Dc is the predicted data.
- If dumpResults is not given, nothing is done (dummy() is called).
out::Int - flag for output (-1: no output, 1: final status, 2: residual norm at each iteration)
Output:
mc - final model
Dc - data
outerFlag - flag for convergence
His - iteration history
#
jInv.InverseSolve.projPCG — Method.
dm = projPCG(H,g,Active,Precond,cgTol,maxIter)
Projected Preconditioned Conjugate Gradient method for solving
H*dm = g subject to dm[!Active] == 0
Input:
H::Function - computes action of Hessian
g::Vector - right hand side
Active - describes active cells
Precond::Function - preconditioner
cgTol - tolerance
maxIter - maximum number of iterations
#
jInv.InverseSolve.projSD — Method.
mc,Dc,outerFlag = projSD(mc,pInv::InverseParam,pMis, indFor = [], dumpResults::Function = dummy)
(Projected) Steepest Descent method for solving
min_x misfit(x) + regularizer(x) subject to x in C
where C is a convex set and a projection operator proj(x) needs to be provided.
Input:
mc::Vector - intial guess for model
pInv::InverseParam - parameter for inversion
pMis - misfit terms
indCredit - indices of forward problems to work on
dumpResults - A function pointer for saving the results throughout the iterations.
- We assume that dumpResults is dumpResults(mc,Dc,iter,pInv,pMis),
- where mc is the recovered model, Dc is the predicted data.
- If dumpResults is not given, nothing is done (dummy() is called).
out::Int - flag for output (-1: no output, 1: final status, 2: residual norm at each iteration)
Output:
mc - final model
Dc - data
outerFlag - flag for convergence
His - iteration history
#
jInv.InverseSolve.smallnessReg — Method.
Rc,dR,d2R = smallnessReg(m,mref,M,Iact=1.0)
Compute smallness regularizer (L2 difference to reference model)
R(m) = 0.5*||m-mref||_V^2
Input:
m - model
mref - reference model
M - Mesh
Output
Rc - value of regularizer
dR - gradient w.r.t. m
d2R - Hessian
#
jInv.InverseSolve.wTVReg — Method.
Rc,dR,d2R = wTVReg(m,mref,M,Iact,C=[])
Compute weighted total variation regularizer
Input:
m - model
mref - reference model
M - Mesh
Iact - projector on active cells
C - anisotropy parameters (default: [1 1 1])
eps - conditioning parameter for TV norm (default: 1e-3)
Output
Rc - value of regularizer
dR - gradient w.r.t. m
d2R - Hessian
#
jInv.InverseSolve.wdiffusionRegNodal — Method.
Rc,dR,d2R = wdiffusionRegNodal(m::Vector, mref::Vector, M::AbstractMesh; Iact=1.0, C=[])
Computes weighted diffusion regularizer for nodal model
Input:
m - model
mref - reference model
M - Mesh
Iact - projector on active cells
C - optional parameters
Output
Rc - value of regularizer
dR - gradient w.r.t. m
d2R - Hessian
#
jInv.InverseSolve.projGrad — Method.
function projGrad
Projects gradient, i.e.,
| gc[i], xl[i] < x[i] < xh[i]
gc[i] = | max(gc[i],0) xc[i] == xh[i]
| min(gc[i],0) xc[i] == xl[i]
Input:
gc - gradient vector
mc - model
boundsLow - lower bounds
boundsHigh - upper bounds