esys.downunder.splitregularizations Package¶
Classes¶
- class esys.downunder.splitregularizations.CostFunction¶
A function f(x) that can be minimized (base class).
Example of usage:
cf=DerivedCostFunction() # ... calculate x ... args=cf.getArguments(x) # this could be potentially expensive! f=cf.getValue(x, *args) # ... it could be required to update x without using the gradient... # ... but then ... gf=cf.getGradient(x, *args)
The class distinguishes between the representation of the solution x (x-type) and the gradients (r-type).
- Note
The provides_inverse_Hessian_approximation class member should be set to
True
in subclasses that provide a valid implementation ofgetInverseHessianApproximation()
- __init__()¶
Constructor. Initializes logger.
- getArguments(x)¶
returns precalculated values that are shared in the calculation of f(x) and grad f(x) and the Hessian operator. The default implementation returns an empty tuple.
Note
The tuple returned by this call will be passed back to this
CostFunction
in other calls(eg:getGradient
). Its contents are not specified at this level because no code, other than theCostFunction
which created it, will be interacting with it. That is, the implementor can put whatever information they find useful in it.- Parameters
x (x-type) – location of derivative
- Return type
tuple
- getDualProduct(x, r)¶
returns the dual product of
x
andr
- Return type
float
- getGradient(x, *args)¶
returns the gradient of f at x using the precalculated values for x.
- Parameters
x (x-type) – location of derivative
args – pre-calculated values for
x
fromgetArguments()
- Return type
r-type
- getInverseHessianApproximation(x, r, *args)¶
returns an approximative evaluation p of the inverse of the Hessian operator of the cost function for a given gradient r at a given location x: H(x) p = r
- Parameters
x (x-type) – location of Hessian operator to be evaluated
r (r-type) – a given gradient
args – pre-calculated values for
x
fromgetArguments()
- Return type
x-type
- Note
In general it is assumed that the Hessian H(x) needs to be calculated in each call for a new location x. However, the solver may suggest that this is not required, typically when the iteration is close to completeness.
- Note
Subclasses that implement this method should set the class variable
provides_inverse_Hessian_approximation
toTrue
to enable the solver to call this method.
- getNorm(x)¶
returns the norm of
x
- Return type
float
- getValue(x, *args)¶
returns the value f(x) using the precalculated values for x.
- Parameters
x (x-type) – a solution approximation
- Return type
float
- provides_inverse_Hessian_approximation = False¶
- updateHessian()¶
notifies the class that the Hessian operator needs to be updated. This method is called by the solver class.
- class esys.downunder.splitregularizations.SplitRegularization(domain, numLevelSets=1, w0=None, w1=None, wc=None, location_of_set_m=<esys.escriptcore.escriptcpp.Data object>, useDiagonalHessianApproximation=False, tol=1e-08, coordinates=None, scale=None, scale_c=None)¶
The regularization term for the level set function
m
within the cost function J for an inversion:J(m)=1/2 * sum_k integrate( mu[k] * ( w0[k] * m_k**2 * w1[k,i] * m_{k,i}**2) + sum_l<k mu_c[l,k] wc[l,k] * | curl(m_k) x curl(m_l) |^2
where w0[k], w1[k,i] and wc[k,l] are non-negative weighting factors and mu[k] and mu_c[l,k] are trade-off factors which may be altered during the inversion. The weighting factors are normalized such that their integrals over the domain are constant:
integrate(w0[k] + inner(w1[k,:],1/L[:]**2))=scale[k] volume(domain)* integrate(wc[l,k]*1/L**4)=scale_c[k] volume(domain) *
- __init__(domain, numLevelSets=1, w0=None, w1=None, wc=None, location_of_set_m=<esys.escriptcore.escriptcpp.Data object>, useDiagonalHessianApproximation=False, tol=1e-08, coordinates=None, scale=None, scale_c=None)¶
initialization.
- Parameters
domain (
Domain
) – domainnumLevelSets (
int
) – number of level setsw0 (
Scalar
ifnumLevelSets
== 1 orData
object of shape (numLevelSets
,) ifnumLevelSets
> 1) – weighting factor for the m**2 term. If not set zero is assumed.w1 (
Vector
ifnumLevelSets
== 1 orData
object of shape (numLevelSets
, DIM) ifnumLevelSets
> 1) – weighting factor for the grad(m_i) terms. If not set zero is assumedwc (
Data
object of shape (numLevelSets
,numLevelSets
)) – weighting factor for the cross gradient terms. If not set zero is assumed. Used for the case ifnumLevelSets
> 1 only. Only valueswc[l,k]
in the lower triangle (l<k) are used.location_of_set_m (
Scalar
ifnumLevelSets
== 1 orData
object of shape (numLevelSets
,) ifnumLevelSets
> 1) – marks location of zero values of the level set functionm
by a positive entry.useDiagonalHessianApproximation (
bool
) – if True cross gradient terms between level set components are ignored when calculating approximations of the inverse of the Hessian Operator. This can speed-up the calculation of the inverse but may lead to an increase of the number of iteration steps in the inversion.tol (positive
float
) – tolerance when solving the PDE for the inverse of the Hessian Operatorcoordinates (ReferenceSystem` or
SpatialCoordinateTransformation
) – defines coordinate system to be usedscale (
Scalar
ifnumLevelSets
== 1 orData
object of shape (numLevelSets
,) ifnumLevelSets
> 1) – weighting factor for level set function variation terms. If not set one is used.scale_c (
Data
object of shape (numLevelSets
,``numLevelSets``)) – scale for the cross gradient terms. If not set one is assumed. Used for the case ifnumLevelSets
> 1 only. Only valuesscale_c[l,k]
in the lower triangle (l<k) are used.
- getArguments(m)¶
- getCoordinateTransformation()¶
returns the coordinate transformation being used
- Return type
CoordinateTransformation
- getDomain()¶
returns the domain of the regularization term
- Return type
Domain
- getDualProduct(m, r)¶
returns the dual product of a gradient represented by X=r[1] and Y=r[0] with a level set function m:
Y_i*m_i + X_ij*m_{i,j}
- Return type
float
- getGradient()¶
returns the gradient of f at x using the precalculated values for x.
- Parameters
x (x-type) – location of derivative
args – pre-calculated values for
x
fromgetArguments()
- Return type
r-type
- getGradientAtPoint()¶
returns the gradient of the cost function J with respect to m.
- Note
This implementation returns Y_k=dPsi/dm_k and X_kj=dPsi/dm_kj
- getInverseHessianApproximationAtPoint(r, solve=True)¶
- getNorm(m)¶
returns the norm of
m
.- Parameters
m (
Data
) – level set function- Return type
float
- getNumLevelSets()¶
returns the number of level set functions
- Return type
int
- getNumTradeOffFactors()¶
returns the number of trade-off factors being used.
- Return type
int
- getPDE()¶
returns the linear PDE to be solved for the Hessian Operator inverse
- Return type
linearPDEs.LinearPDE
- getValue(m, grad_m)¶
returns the value of the cost function J with respect to m. This equation is specified in the inversion cookbook.
- Return type
float
- getValueAtPoint()¶
returns the value of the cost function J with respect to m. This equation is specified in the inversion cookbook.
- Return type
float
- setPoint(m)¶
sets the point which this function will work with
- Parameters
m (
Data
) – level set function
- setTradeOffFactors(mu=None)¶
sets the trade-off factors for the level-set variation and the cross-gradient.
- Parameters
mu (
list
offloat
or`numpy.array`
) – new values for the trade-off factors where values mu[:numLevelSets] are the trade-off factors for the level-set variation and the remaining values for the cross-gradient part with mu_c[l,k]=mu[numLevelSets+l+((k-1)*k)/2] (l<k). If no values for mu are given ones are used. Values must be positive.
- setTradeOffFactorsForCrossGradient(mu_c=None)¶
sets the trade-off factors for the cross-gradient terms.
- Parameters
mu_c (
float
,list
offloat
ornumpy.array
) – new values for the trade-off factors for the cross-gradient terms. Values must be positive. If no value is given ones are used. Only value mu_c[l,k] for l<k are used.
- setTradeOffFactorsForVariation(mu=None)¶
sets the trade-off factors for the level-set variation part.
- Parameters
mu (
float
,list
offloat
or`numpy.array`
) – new values for the trade-off factors. Values must be positive.
- updateHessian()¶
notifies the class to recalculate the Hessian operator.
Functions¶
- esys.downunder.splitregularizations.makeTransformation(domain, coordinates=None)¶
returns a
SpatialCoordinateTransformation
for the given domain- Parameters
domain (
esys.escript.AbstractDomain
) – domain in the domain of the coordinate transformationcoordinates (
ReferenceSystem
orSpatialCoordinateTransformation
) – the reference system or spatial coordinate system.
- Returns
the spatial coordinate system for the given domain of the specified reference system
coordinates
. Ifcoordinates
is already spatial coordinate system based on the riven domaincoordinates
is returned. Otherwise an appropriate spatial coordinate system is created.- Return type
SpatialCoordinateTransformation