hippylib.modeling package

Submodules

hippylib.modeling.PDEProblem module

class hippylib.modeling.PDEProblem.PDEProblem[source]

Bases: object

Consider the PDE problem: Given \(m\), find \(u\) such that

\[F(u, m, p) = ( f(u, m), p) = 0, \quad \forall p.\]

Here \(F\) is linear in \(p\), but it may be non linear in \(u\) and \(m\).

apply_ij(i, j, dir, out)[source]

Given \(u, m, p\); compute \(\delta_{ij} F(u, m, p; \hat{i}, \tilde{j})\) in the direction \(\tilde{j} =\) dir, \(\forall \hat{i}.\)

apply_ijk(i, j, k, x, jdir, kdir, out)[source]

Given x = [u,a,p]; compute \(\delta_{ijk} F(u,a,p; \hat{i}, \tilde{j}, \tilde{k})\) in the direction \((\tilde{j},\tilde{k}) = (\)jdir,kdir), \(\forall \hat{i}.\)

evalGradientParameter(x, out)[source]

Given \(u, m, p\); evaluate \(\delta_m F(u, m, p; \hat{m}),\, \forall \hat{m}.\)

generate_parameter()[source]

Return a vector in the shape of the parameter.

generate_state()[source]

Return a vector in the shape of the state.

init_parameter(m)[source]

Initialize the parameter.

setLinearizationPoint(x, gauss_newton_approx)[source]

Set the values of the state and parameter for the incremental forward and adjoint solvers. Set whether Gauss Newton approximation of the Hessian should be used.

solveAdj(adj, x, adj_rhs)[source]

Solve the linear adjoint problem: Given \(m\), \(u\); find \(p\) such that

\[\delta_u F(u, m, p;\hat{u}) = 0, \quad \forall \hat{u}.\]
solveFwd(state, x)[source]

Solve the possibly nonlinear forward problem: Given \(m\), find \(u\) such that

\[\delta_p F(u, m, p;\hat{p}) = 0, \quad \forall \hat{p}.\]
solveIncremental(out, rhs, is_adj)[source]

If is_adj = False:

Solve the forward incremental system: Given \(u, m\), find \(\tilde{u}\) such that

\[\delta_{pu} F(u, m, p; \hat{p}, \tilde{u}) = \mbox{rhs}, \quad \forall \hat{p}.\]

If is_adj = True:

Solve the adjoint incremental system: Given \(u, m\), find \(\tilde{p}\) such that

\[\delta_{up} F(u, m, p; \hat{u}, \tilde{p}) = \mbox{rhs}, \quad \forall \hat{u}.\]
class hippylib.modeling.PDEProblem.PDEVariationalProblem(Vh, varf_handler, bc, bc0, is_fwd_linear=False)[source]

Bases: hippylib.modeling.PDEProblem.PDEProblem

_createLUSolver()[source]
apply_ij(i, j, dir, out)[source]

Given \(u, m, p\); compute \(\delta_{ij} F(u, m, p; \hat{i}, \tilde{j})\) in the direction \(\tilde{j} =\) dir, \(\forall \hat{i}\).

apply_ijk(i, j, k, x, jdir, kdir, out)[source]

Given x = [u,a,p]; compute \(\delta_{ijk} F(u,a,p; \hat{i}, \tilde{j}, \tilde{k})\) in the direction \((\tilde{j},\tilde{k}) = (\)jdir,kdir), \(\forall \hat{i}.\)

evalGradientParameter(x, out)[source]

Given \(u, m, p\); evaluate \(\delta_m F(u, m, p; \hat{m}),\, \forall \hat{m}.\)

generate_parameter()[source]

Return a vector in the shape of the parameter.

generate_state()[source]

Return a vector in the shape of the state.

init_parameter(m)[source]

Initialize the parameter.

setLinearizationPoint(x, gauss_newton_approx)[source]

Set the values of the state and parameter for the incremental forward and adjoint solvers.

solveAdj(adj, x, adj_rhs)[source]

Solve the linear adjoint problem: Given \(m, u\); find \(p\) such that

\[\delta_u F(u, m, p;\hat{u}) = 0, \quad \forall \hat{u}.\]
solveFwd(state, x)[source]

Solve the possibly nonlinear forward problem: Given \(m\), find \(u\) such that

\[\delta_p F(u, m, p;\hat{p}) = 0,\quad \forall \hat{p}.\]
solveIncremental(out, rhs, is_adj)[source]

If is_adj == False:

Solve the forward incremental system: Given \(u, m\), find \(\tilde{u}\) such that

\[\delta_{pu} F(u, m, p; \hat{p}, \tilde{u}) = \mbox{rhs},\quad \forall \hat{p}.\]

If is_adj == True:

Solve the adjoint incremental system: Given \(u, m\), find \(\tilde{p}\) such that

\[\delta_{up} F(u, m, p; \hat{u}, \tilde{p}) = \mbox{rhs},\quad \forall \hat{u}.\]

hippylib.modeling.expression module

hippylib.modeling.misfit module

class hippylib.modeling.misfit.ContinuousStateObservation(Vh, dX, bcs, form=None)[source]

Bases: hippylib.modeling.misfit.Misfit

This class implements continuous state observations in a subdomain \(X \subset \Omega\) or \(X \subset \partial \Omega\).

Constructor:

Vh: the finite element space for the state variable.

dX: the integrator on subdomain X where observation are presents. E.g. dX = ufl.dx means observation on all \(\Omega\) and dX = ufl.ds means observations on all \(\partial \Omega\).

bcs: If the forward problem imposes Dirichlet boundary conditions \(u = u_D \mbox{ on } \Gamma_D\); bcs is a list of dolfin.DirichletBC object that prescribes homogeneuos Dirichlet conditions \(u = 0 \mbox{ on } \Gamma_D\).

form: if form = None we compute the \(L^2(X)\) misfit: \(\int_X (u - u_d)^2 dX,\) otherwise the integrand specified in the given form will be used.

apply_ij(i, j, dir, out)[source]

Apply the second variation \(\delta_{ij}\) (i,j = STATE,PARAMETER) of the cost in direction dir.

cost(x)[source]

Given x evaluate the cost functional. Only the state u and (possibly) the parameter m are accessed.

grad(i, x, out)[source]

Given the state and the paramter in x, compute the partial gradient of the misfit functional in with respect to the state (i == STATE) or with respect to the parameter (i == PARAMETER).

setLinearizationPoint(x, gauss_newton_approx=False)[source]

Set the point for linearization.

Inputs:

x=[u, m, p] - linearization point

gauss_newton_approx (bool) - whether to use Gauss Newton approximation

class hippylib.modeling.misfit.DiscreteStateObservation(B, data=None, noise_variance=None)[source]

Bases: hippylib.modeling.misfit.Misfit

This class define a misfit function for a discrete linear observation operator B

Constructor:
B is the observation operator data is the data noise_variance is the variance of the noise
apply_ij(i, j, dir, out)[source]

Apply the second variation \(\delta_{ij}\) (i,j = STATE,PARAMETER) of the cost in direction dir.

cost(x)[source]

Given x evaluate the cost functional. Only the state u and (possibly) the parameter m are accessed.

grad(i, x, out)[source]

Given the state and the paramter in x, compute the partial gradient of the misfit functional in with respect to the state (i == STATE) or with respect to the parameter (i == PARAMETER).

setLinearizationPoint(x, gauss_newton_approx=False)[source]

Set the point for linearization.

Inputs:

x=[u, m, p] - linearization point

gauss_newton_approx (bool) - whether to use Gauss Newton approximation

class hippylib.modeling.misfit.Misfit[source]

Bases: object

Abstract class to model the misfit component of the cost functional. In the following x will denote the variable [u, m, p], denoting respectively the state u, the parameter m, and the adjoint variable p.

The methods in the class misfit will usually access the state u and possibly the parameter m. The adjoint variables will never be accessed.

apply_ij(i, j, dir, out)[source]

Apply the second variation \(\delta_{ij}\) (i,j = STATE,PARAMETER) of the cost in direction dir.

cost(x)[source]

Given x evaluate the cost functional. Only the state u and (possibly) the parameter m are accessed.

grad(i, x, out)[source]

Given the state and the paramter in x, compute the partial gradient of the misfit functional in with respect to the state (i == STATE) or with respect to the parameter (i == PARAMETER).

setLinearizationPoint(x, gauss_newton_approx=False)[source]

Set the point for linearization.

Inputs:

x=[u, m, p] - linearization point

gauss_newton_approx (bool) - whether to use Gauss Newton approximation

class hippylib.modeling.misfit.MultDiscreteStateObservation(B, data=None, Mpar=1.0)[source]

Bases: hippylib.modeling.misfit.Misfit

This class implements discrete state observations defined by the linear operator B. A multiplicative Gamma(M,M) noise model is assumed

Constructor:

B is the observation operator

Mpar Gamma distribution parameter

apply_ij(i, j, dir, out)[source]

Apply the second variation \(\delta_{ij}\) (i,j = STATE,PARAMETER) of the cost in direction dir.

cost(x)[source]

Given x evaluate the cost functional. Only the state u and (possibly) the parameter m are accessed.

grad(i, x, out)[source]

Given the state and the paramter in x, compute the partial gradient of the misfit functional in with respect to the state (i == STATE) or with respect to the parameter (i == PARAMETER).

setLinearizationPoint(x, gauss_newton_approx=False)[source]

Set the point for linearization.

Inputs:

x=[u, m, p] - linearization point

gauss_newton_approx (bool) - whether to use Gauss Newton approximation

hippylib.modeling.misfit.MultPointwiseStateObservation(Vh, obs_points, Mpar)[source]

This function returns an instance of DiscreteStateObservation for pointwise state observations at given locations. A multiplicative Gamma(M,M) noise model is assumed Inputs:

Vh is the finite element space for the state variable data is the data obs_points is a 2D array number of points by geometric dimensions that stores the location of the observations.

Mpar Gamma distribution parameter

class hippylib.modeling.misfit.MultiStateMisfit(misfits)[source]

Bases: hippylib.modeling.misfit.Misfit

append(misfit)[source]
apply_ij(i, j, dir, out)[source]

Apply the second variation \(\delta_{ij}\) (i,j = STATE,PARAMETER) of the cost in direction dir.

cost(x)[source]

Given x evaluate the cost functional. Only the state u and (possibly) the parameter m are accessed.

grad(i, x, out)[source]

Given the state and the paramter in x, compute the partial gradient of the misfit functional in with respect to the state (i == STATE) or with respect to the parameter (i == PARAMETER).

setLinearizationPoint(x, gauss_newton_approx=False)[source]

Set the point for linearization.

Inputs:

x=[u, m, p] - linearization point

gauss_newton_approx (bool) - whether to use Gauss Newton approximation

hippylib.modeling.misfit.PointwiseStateObservation(Vh, obs_points)[source]

This function returns an instance of DiscreteStateObservation for pointwise state observations at given locations. Inputs:

Vh is the finite element space for the state variable

obs_points is a 2D array number of points by geometric dimensions that stores the location of the observations.

hippylib.modeling.model module

class hippylib.modeling.model.Model(problem, prior, misfit)[source]

Bases: object

This class contains the full description of the inverse problem. As inputs it takes a PDEProblem object, a Prior object, and a Misfit object.

In the following we will denote with

  • u the state variable
  • m the (model) parameter variable
  • p the adjoint variable

Create a model given:

  • problem: the description of the forward/adjoint problem and all the sensitivities
  • prior: the prior component of the cost functional
  • misfit: the misfit componenent of the cost functional
Rsolver()[source]

Return an object Rsovler that is a suitable solver for the regularization operator \(R\).

The solver object should implement the method Rsolver.solve(z,r) such that \(Rz pprox r\).

applyC(dm, out)[source]

Apply the \(C\) block of the Hessian to a (incremental) parameter variable, i.e. out = \(C dm\)

Parameters:

  • dm the (incremental) parameter variable
  • out the action of the \(C\) block on dm

Note

This routine assumes that out has the correct shape.

applyCt(dp, out)[source]

Apply the transpose of the \(C\) block of the Hessian to a (incremental) adjoint variable. out = \(C^t dp\)

Parameters:

  • dp the (incremental) adjoint variable
  • out the action of the \(C^T\) block on dp

..note:: This routine assumes that out has the correct shape.

applyR(dm, out)[source]

Apply the regularization \(R\) to a (incremental) parameter variable. out = \(R dm\)

Parameters:

  • dm the (incremental) parameter variable
  • out the action of \(R\) on dm

Note

This routine assumes that out has the correct shape.

applyWmm(dm, out)[source]

Apply the \(W_{mm}\) block of the Hessian to a (incremental) parameter variable. out = \(W_{mm} dm\)

Parameters:

  • dm the (incremental) parameter variable
  • out the action of the \(W_{mm}\) on block dm

Note

This routine assumes that out has the correct shape.

applyWmu(du, out)[source]

Apply the \(W_{mu}\) block of the Hessian to a (incremental) state variable. out = \(W_{mu} du\)

Parameters:

  • du the (incremental) state variable
  • out the action of the \(W_{mu}\) block on du

Note

This routine assumes that out has the correct shape.

applyWum(dm, out)[source]

Apply the \(W_{um}\) block of the Hessian to a (incremental) parameter variable. out = \(W_{um} dm\)

Parameters:

  • dm the (incremental) parameter variable
  • out the action of the \(W_{um}\) block on du

Note

This routine assumes that out has the correct shape.

applyWuu(du, out)[source]

Apply the \(W_{uu}\) block of the Hessian to a (incremental) state variable. out = \(W_{uu} du\)

Parameters:

  • du the (incremental) state variable
  • out the action of the \(W_{uu}\) block on du

Note

This routine assumes that out has the correct shape.

apply_ij(i, j, d, out)[source]
cost(x)[source]

Given the list x = [u,m,p] which describes the state, parameter, and adjoint variable compute the cost functional as the sum of the misfit functional and the regularization functional.

Return the list [cost functional, regularization functional, misfit functional]

Note

p is not needed to compute the cost functional

evalGradientParameter(x, mg, misfit_only=False)[source]

Evaluate the gradient for the variational parameter equation at the point x=[u,m,p].

Parameters:

  • x = [u,m,p] the point at which to evaluate the gradient.
  • mg the variational gradient \((g, mtest)\), mtest being a test function in the parameter space (Output parameter)

Returns the norm of the gradient in the correct inner product \(g_norm = sqrt(g,g)\)

generate_vector(component='ALL')[source]

By default, return the list [u,m,p] where:

  • u is any object that describes the state variable
  • m is a dolfin.Vector object that describes the parameter variable. (Needs to support linear algebra operations)
  • p is any object that describes the adjoint variable

If component = STATE return only u

If component = PARAMETER return only m

If component = ADJOINT return only p

init_parameter(m)[source]

Reshape m so that it is compatible with the parameter variable

setPointForHessianEvaluations(x, gauss_newton_approx=False)[source]

Specify the point x = [u,m,p] at which the Hessian operator (or the Gauss-Newton approximation) needs to be evaluated.

Parameters:

  • x = [u,m,p]: the point at which the Hessian or its Gauss-Newton approximation needs to be evaluated.
  • gauss_newton_approx (bool): whether to use Gauss-Newton approximation (default: use Newton)

Note

This routine should either:

  • simply store a copy of x and evaluate action of blocks of the Hessian on the fly
  • or partially precompute the block of the hessian (if feasible)
solveAdj(out, x)[source]

Solve the linear adjoint problem.

Parameters:

  • out: is the solution of the adjoint problem (i.e. the adjoint p) (Output parameter)

  • x = [u, m, p] provides

    1. the parameter variable m for assembling the adjoint operator
    2. the state variable u for assembling the adjoint right hand side

    Note

    p is not accessed

solveAdjIncremental(sol, rhs)[source]

Solve the incremental adjoint problem for a given right-hand side

Parameters:

  • sol the solution of the incremental adjoint problem (Output)
  • rhs the right hand side of the linear system
solveFwd(out, x)[source]

Solve the (possibly non-linear) forward problem.

Parameters:

  • out: is the solution of the forward problem (i.e. the state) (Output parameters)

  • x = [u,m,p] provides

    1. the parameter variable m for the solution of the forward problem
    2. the initial guess u if the forward problem is non-linear

    Note

    p is not accessed.

solveFwdIncremental(sol, rhs)[source]

Solve the linearized (incremental) forward problem for a given right-hand side

Parameters:

  • sol the solution of the linearized forward problem (Output)
  • rhs the right hand side of the linear system

hippylib.modeling.modelVerify module

hippylib.modeling.modelVerify.modelVerify(model, m0, is_quadratic=False, misfit_only=False, verbose=True, eps=None)[source]

Verify the reduced Gradient and the Hessian of a model. It will produce two loglog plots of the finite difference checks for the gradient and for the Hessian. It will also check for symmetry of the Hessian.

hippylib.modeling.modelVerify.modelVerifyPlotErrors(is_quadratic, eps, err_grad, err_H)[source]

hippylib.modeling.pointwiseObservation module

hippylib.modeling.pointwiseObservation.assemblePointwiseObservation(Vh, targets, prune_and_sort=False)[source]

Assemble the pointwise observation matrix:

Inputs

  • Vh: FEniCS finite element space.
  • targets: observation points (numpy array).
hippylib.modeling.pointwiseObservation.exportPointwiseObservation(Vh, B, data, fname, varname='observation')[source]

This function writes a VTK PolyData file to visualize pointwise data.

Inputs:

  • B: observation operator.
  • mesh: mesh.
  • data: dolfin.Vector containing the data.
  • fname: filename for the file to export (without extension).
  • varname: name of the variable for the .vtk file.
hippylib.modeling.pointwiseObservation.write_vtk(points, data, fname, varname='observation')[source]

This function writes a VTK PolyData file to visualize pointwise data.

Inputs:

  • points: locations of the points (numpy array of size equal to number of points times space dimension).
  • data: pointwise values (numpy array of size equal to number of points).
  • fname: filename for the .vtk file to export.
  • varname: name of the variable for the .vtk file.

hippylib.modeling.posterior module

class hippylib.modeling.posterior.GaussianLRPosterior(prior, d, U, mean=None)[source]

Bases: object

Class for the low rank Gaussian Approximation of the Posterior. This class provides functionality for approximate Hessian apply, solve, and Gaussian sampling based on the low rank factorization of the Hessian.

In particular if \(d\) and \(U\) are the dominant eigenpairs of \(H_{\mbox{misfit}} U[:,i] = d[i] R U[:,i]\) then we have:

  • low rank Hessian apply: \(y = ( R + RU D U^{T}) x.\)
  • low rank Hessian solve: \(y = (R^-1 - U (I + D^{-1})^{-1} U^T) x.\)
  • low rank Hessian Gaussian sampling: \(y = ( I - U S U^{T}) x\), where \(S = I - (I + D)^{-1/2}\) and \(x \sim \mathcal{N}(0, R^{-1}).\)

Construct the Gaussian approximation of the posterior. Input: - prior: the prior mode. - d: the dominant generalized eigenvalues of the Hessian misfit. - U: the dominant generalized eigenvector of the Hessian misfit \(U^T R U = I.\) - mean: the MAP point.

_sample_given_prior(s_prior, s_post)[source]
_sample_given_white_noise(noise, s_prior, s_post)[source]
cost(m)[source]
init_vector(x, dim)[source]

Inizialize a vector x to be compatible with the range/domain of \(H\). If dim == "noise" inizialize x to be compatible with the size of white noise used for sampling.

klDistanceFromPrior(sub_comp=False)[source]
pointwise_variance(**kwargs)[source]

Compute/estimate the pointwise variance of the posterior, prior distribution and the pointwise variance reduction informed by the data.

See _Prior.pointwise_variance for more details.

sample(*args, **kwargs)[source]

possible calls:

  1. sample(s_prior, s_post, add_mean=True)

    Given a prior sample s_prior compute a sample s_post from the posterior.

    • s_prior is a sample from the prior centered at 0 (input).
    • s_post is a sample from the posterior (output).
    • if add_mean=True (default) then the samples will be centered at the map point.
  2. sample(noise, s_prior, s_post, add_mean=True)

    Given noise \(\sim \mathcal{N}(0, I)\) compute a sample s_prior from the prior and s_post from the posterior.

    • noise is a realization of white noise (input).
    • s_prior is a sample from the prior (output).
    • s_post is a sample from the posterior.
    • if add_mean=True (default) then the prior and posterior samples will be centered at the respective means.
trace(**kwargs)[source]

Compute/estimate the trace of the posterior, prior distribution and the trace of the data informed correction.

See _Prior.trace for more details.

trace_update()[source]
class hippylib.modeling.posterior.LowRankHessian(prior, d, U)[source]

Bases: object

Operator that represents the action of the low rank approximation of the Hessian and of its inverse.

init_vector(x, dim)[source]
inner(x, y)[source]
mult(x, y)[source]
solve(sol, rhs)[source]
class hippylib.modeling.posterior.LowRankPosteriorSampler(prior, d, U)[source]

Bases: object

Object to sample from the low rank approximation of the posterior.

\[y = ( I - U S U^{T}) x,\]

where

\(S = I - (I + D)^{-1/2}, x \sim \mathcal{N}(0, R^{-1}).\)

init_vector(x, dim)[source]
sample(noise, s)[source]

hippylib.modeling.prior module

hippylib.modeling.prior.BiLaplacianComputeCoefficients(sigma2, rho, ndim)[source]

This class is responsible to compute the parameters gamma and delta for the BiLaplacianPrior given the marginal variance sigma2 and correlation length rho. ndim is the dimension of the domain 2D or 3D

hippylib.modeling.prior.BiLaplacianPrior(Vh, gamma, delta, Theta=None, mean=None, rel_tol=1e-12, max_iter=1000, robin_bc=False)[source]

This function construct an instance of :code”SqrtPrecisionPDE_Prior with covariance matrix \(C = (\delta I + \gamma \mbox{div } \Theta \nabla) ^ {-2}\).

The magnitude of \(\delta\gamma\) governs the variance of the samples, while the ratio \(\frac{\gamma}{\delta}\) governs the correlation lenght.

Here \(\Theta\) is a SPD tensor that models anisotropy in the covariance kernel.

Input:

  • Vh: the finite element space for the parameter
  • gamma and delta: the coefficient in the PDE (floats, dl.Constant, dl.Expression, or dl.Function)
  • Theta: the SPD tensor for anisotropic diffusion of the PDE
  • mean: the prior mean
  • rel_tol: relative tolerance for solving linear systems involving covariance matrix
  • max_iter: maximum number of iterations for solving linear systems involving covariance matrix
  • robin_bc: whether to use Robin boundary condition to remove boundary artifacts
class hippylib.modeling.prior.GaussianRealPrior(Vh, covariance, mean=None)[source]

Bases: hippylib.modeling.prior._Prior

This class implements a finite-dimensional Gaussian prior, \(\mathcal{N}(\boldsymbol{m}, \boldsymbol{C})\), where \(\boldsymbol{m}\) is the mean of the Gaussian distribution, and \(\boldsymbol{C}\) is its covariance. The underlying finite element space is assumed to be the “R” space.

Constructor

Inputs: - Vh: Finite element space on which the prior is

defined. Must be the Real space with one global degree of freedom
  • covariance: The covariance of the prior. Must be a
    numpy.ndarray of appropriate size
  • :code:`mean`(optional): Mean of the prior distribution. Must be of
    type dolfin.Vector()
init_vector(x, dim)[source]

Inizialize a vector x to be compatible with the range/domain of \(R\).

If dim == "noise" inizialize x to be compatible with the size of white noise used for sampling.

sample(noise, s, add_mean=True)[source]

Given noise \(\sim \mathcal{N}(0, I)\) compute a sample s from the prior.

If add_mean == True add the prior mean value to s.

class hippylib.modeling.prior.LaplacianPrior(Vh, gamma, delta, mean=None, rel_tol=1e-12, max_iter=100)[source]

Bases: hippylib.modeling.prior._Prior

This class implements a prior model with covariance matrix \(C = (\delta I - \gamma \Delta) ^ {-1}\).

The magnitude of \(\gamma\) governs the variance of the samples, while the ratio \(\frac{\gamma}{\delta}\) governs the correlation length.

Note

\(C\) is a trace class operator only in 1D while it is not a valid prior in 2D and 3D.

Construct the prior model. Input:

  • Vh: the finite element space for the parameter
  • gamma and delta: the coefficient in the PDE
  • Theta: the SPD tensor for anisotropic diffusion of the PDE
  • mean: the prior mean
init_vector(x, dim)[source]

Inizialize a vector x to be compatible with the range/domain of \(R\).

If dim == "noise" inizialize x to be compatible with the size of white noise used for sampling.

sample(noise, s, add_mean=True)[source]

Given noise \(\sim \mathcal{N}(0, I)\) compute a sample s from the prior.

If add_mean == True add the prior mean value to s.

hippylib.modeling.prior.MollifiedBiLaplacianPrior(Vh, gamma, delta, locations, m_true, Theta=None, pen=10.0, order=2, rel_tol=1e-12, max_iter=1000)[source]

This function construct an instance of :code”SqrtPrecisionPDE_Prior with covariance matrix \(C = \left( [\delta + \mbox{pen} \sum_i m(x - x_i) ] I + \gamma \mbox{div } \Theta \nabla\right) ^ {-2}\),

where

  • \(\Theta\) is a SPD tensor that models anisotropy in the covariance kernel.
  • \(x_i (i=1,...,n)\) are points were we assume to know exactly the value of the parameter (i.e., \(m(x_i) = m_{\mbox{true}}( x_i) \mbox{ for } i=1,...,n).\)
  • \(m\) is the mollifier function: \(m(x - x_i) = \exp\left( - \left[\frac{\gamma}{\delta}\| x - x_i \|_{\Theta^{-1}}\right]^{\mbox{order}} \right).\)
  • pen is a penalization parameter.

The magnitude of \(\delta \gamma\) governs the variance of the samples, while the ratio \(\frac{\gamma}{\delta}\) governs the correlation length.

The prior mean is computed by solving

\[\left( [\delta + \sum_i m(x - x_i) ] I + \gamma \mbox{div } \Theta \nabla \right) m = \sum_i m(x - x_i) m_{\mbox{true}}.\]

Input:

  • Vh: the finite element space for the parameter
  • gamma and delta: the coefficients in the PDE
  • locations: the points \(x_i\) at which we assume to know the true value of the parameter
  • m_true: the true model
  • Theta: the SPD tensor for anisotropic diffusion of the PDE
  • pen: a penalization parameter for the mollifier
class hippylib.modeling.prior.SqrtPrecisionPDE_Prior(Vh, sqrt_precision_varf_handler, mean=None, rel_tol=1e-12, max_iter=1000)[source]

Bases: hippylib.modeling.prior._Prior

This class implement a prior model with covariance matrix \(C = A^{-1} M A^-1\), where A is the finite element matrix arising from discretization of sqrt_precision_varf_handler

Construct the prior model. Input:

  • Vh: the finite element space for the parameter
  • code:sqrt_precision_varf_handler:
     the PDE representation of the sqrt of the covariance operator
  • mean: the prior mean
init_vector(x, dim)[source]

Inizialize a vector x to be compatible with the range/domain of \(R\).

If dim == "noise" inizialize x to be compatible with the size of white noise used for sampling.

sample(noise, s, add_mean=True)[source]

Given noise \(\sim \mathcal{N}(0, I)\) compute a sample s from the prior.

If add_mean == True add the prior mean value to s.

class hippylib.modeling.prior._BilaplacianR(A, Msolver)[source]

Bases: object

Operator that represent the action of the regularization/precision matrix for the Bilaplacian prior.

init_vector(x, dim)[source]
mpi_comm()[source]
mult(x, y)[source]
class hippylib.modeling.prior._BilaplacianRsolver(Asolver, M)[source]

Bases: object

Operator that represent the action of the inverse the regularization/precision matrix for the Bilaplacian prior.

init_vector(x, dim)[source]
solve(x, b)[source]
class hippylib.modeling.prior._Prior[source]

Bases: object

Abstract class to describe the prior model. Concrete instances of a _Prior class should expose the following attributes and methods.

Attributes:

  • R: an operator to apply the regularization/precision operator.
  • Rsolver: an operator to apply the inverse of the regularization/precision operator.
  • M: the mass matrix in the control space.
  • mean: the prior mean.

Methods:

  • init_vector(self,x,dim): Inizialize a vector x to be compatible with the range/domain of R If dim == "noise" inizialize x to be compatible with the size of white noise used for sampling.
  • sample(self, noise, s, add_mean=True): Given noise \(\sim \mathcal{N}(0, I)\) compute a sample s from the prior. If add_mean==True add the prior mean value to s.
cost(m)[source]
getHessianPreconditioner()[source]

Return the preconditioner for Newton-CG

grad(m, out)[source]
init_vector(x, dim)[source]
pointwise_variance(method, k=1000000, r=200)[source]

Compute/estimate the prior pointwise variance.

  • If method=="Exact" we compute the diagonal entries of \(R^{-1}\) entry by entry. This requires to solve \(n\) linear system in \(R\) (not scalable, but ok for illustration purposes).
sample(noise, s, add_mean=True)[source]
trace(method='Exact', tol=0.1, min_iter=20, max_iter=100, r=200)[source]

Compute/estimate the trace of the prior covariance operator.

  • If method=="Exact" we compute the trace exactly by summing the diagonal entries of \(R^{-1}M\). This requires to solve \(n\) linear system in \(R\) (not scalable, but ok for illustration purposes).
  • If method=="Estimator" use the trace estimator algorithms implemeted in the class TraceEstimator. tol is a relative bound on the estimator standard deviation. In particular, we used enough samples in the Estimator such that the standard deviation of the estimator is less then tol\(tr(\mbox{Prior})\). min_iter and max_iter are the lower and upper bound on the number of samples to be used for the estimation of the trace.
class hippylib.modeling.prior._RinvM(Rsolver, M)[source]

Bases: object

Operator that models the action of \(R^{-1}M\). It is used in the randomized trace estimator.

init_vector(x, dim)[source]
mult(x, y)[source]

hippylib.modeling.reducedHessian module

class hippylib.modeling.reducedHessian.FDHessian(model, m0, h, misfit_only=False)[source]

Bases: object

This class implements matrix free application of the reduced Hessian operator. The constructor takes the following parameters:

  • model: the object which contains the description of the problem.
  • m0: the value of the parameter at which the Hessian needs to be evaluated.
  • h: the mesh size for FD.
  • misfit_only: a boolean flag that describes whenever the full Hessian or only the misfit component of the Hessian is used.

Type help(Template) for more information on which methods model should implement.

Construct the reduced Hessian Operator

init_vector(x, dim)[source]

Reshape the Vector x so that it is compatible with the reduced Hessian operator.

Parameters:

  • x: the vector to reshape
  • dim: if 0 then x will be reshaped to be compatible with the range of the reduced Hessian, if 1 then x will be reshaped to be compatible with the domain of the reduced Hessian

Note

Since the reduced Hessian is a self adjoint operator, the range and the domain is the same. Either way, we choosed to add the parameter dim for consistency with the interface of Matrix in dolfin.

inner(x, y)[source]

Perform the inner product between x and y in the norm induced by the reduced Hessian \(H,\, (x, y)_H = x' H y\).

mult(x, y)[source]

Apply the reduced Hessian (or the Gauss-Newton approximation) to the vector x. Return the result in y.

class hippylib.modeling.reducedHessian.ReducedHessian(model, misfit_only=False)[source]

Bases: object

This class implements matrix free application of the reduced Hessian operator. The constructor takes the following parameters:

  • model: the object which contains the description of the problem.
  • misfit_only: a boolean flag that describes whenever the full Hessian or only the misfit component of the Hessian is used.

Type help(modelTemplate) for more information on which methods model should implement.

Construct the reduced Hessian Operator

GNHessian(x, y)[source]

Apply the Gauss-Newton approximation of the reduced Hessian to the vector x. Return the result in y.

TrueHessian(x, y)[source]

Apply the the reduced Hessian to the vector x. Return the result in y.

init_vector(x, dim)[source]

Reshape the Vector x so that it is compatible with the reduced Hessian operator.

Parameters:

  • x: the vector to reshape.
  • dim: if 0 then x will be reshaped to be compatible with the range of the reduced Hessian, if 1 then x will be reshaped to be compatible with the domain of the reduced Hessian.

Note

Since the reduced Hessian is a self adjoint operator, the range and the domain is the same. Either way, we choosed to add the parameter dim for consistency with the interface of Matrix in dolfin.

inner(x, y)[source]

Perform the inner product between x and y in the norm induced by the reduced Hessian \(H,\,(x, y)_H = x' H y\).

mult(x, y)[source]

Apply the reduced Hessian (or the Gauss-Newton approximation) to the vector x. Return the result in y.

hippylib.modeling.timeDependentVector module

class hippylib.modeling.timeDependentVector.TimeDependentVector(times, tol=1e-10, mpi_comm=<sphinx.ext.autodoc.importer._MockObject object>)[source]

Bases: object

A class to store time dependent vectors. Snapshots are stored/retrieved by specifying the time of the snapshot. Times at which the snapshot are taken must be specified in the constructor.

Constructor:

  • times: time frame at which snapshots are stored.
  • tol : tolerance to identify the frame of the snapshot.
axpy(a, other)[source]

Compute \(x = x + \mbox{a*other}\) snapshot per snapshot.

copy()[source]

Return a copy of all the time frames and snapshots

initialize(M, dim)[source]

Initialize all the snapshot to be compatible with the range/domain of an operator M.

inner(other)[source]

Compute the inner products: \(a+= (\mbox{self[i]},\mbox{other[i]})\) for each snapshot.

norm(time_norm, space_norm)[source]

Compute the space-time norm of the snapshot.

retrieve(u, t)[source]

Retrieve snapshot u relative to time t. If t does not belong to the list of time frame an error is raised.

store(u, t)[source]

Store snapshot u relative to time t. If t does not belong to the list of time frame an error is raised.

zero()[source]

Zero out each snapshot.

hippylib.modeling.variables module

Enumerator for the variables of the inverse problem: - the STATE, PARAMETER, and ADJOINT variables.

Module contents