Code Overview#
Executive summary#
The fundamental building blocks required for IFT computations are best recognized from a large distance, ignoring all technical details.
From such a perspective,
IFT problems largely consist of the combination of several high dimensional minimization problems.
Within NIFTy, operators are used to define the characteristic equations and properties of the problems.
The equations are built mostly from the application of linear operators, but there may also be nonlinear functions involved.
The unknowns in the equations represent either continuous physical fields, or they are simply individual measured data points.
Discretized fields have geometrical information (like locations and volume elements) associated with every entry; this information is called the field’s domain.
In the following sections, the concepts briefly presented here will be discussed in more detail; this is done in reversed order of their introduction, to avoid forward references.
Domains#
Abstract base class#
One of the fundamental building blocks of the NIFTy8 framework is the domain.
Its required capabilities are expressed by the abstract Domain
class.
A domain must be able to answer the following queries:
Unstructured domains#
Domains can be either structured (i.e. there is geometrical information associated with them, like position in space and volume factors), or unstructured (meaning that the data points have no associated manifold).
Unstructured domains can be described by instances of NIFTy’s
UnstructuredDomain
class.
Structured domains#
In contrast to unstructured domains, these domains have an assigned geometry.
NIFTy requires them to provide the volume elements of their grid cells.
The additional methods are specified in the abstract class
StructuredDomain
:
The properties
scalar_dvol
,dvol
, andtotal_volume
provide information about the domain’s pixel volume(s) and its total volume.The property
harmonic
specifies whether a domain is harmonic (i.e. describes a frequency space) or notIf (and only if) the domain is harmonic, the methods
get_k_length_array()
,get_unique_k_lengths()
, andget_fft_smoothing_kernel_function()
provide absolute distances of the individual grid cells from the origin and assist with Gaussian convolution.
NIFTy comes with several concrete subclasses of StructuredDomain
:
RGSpace
represents a regular Cartesian grid with an arbitrary number of dimensions, which is supposed to be periodic in each dimension.LogRGSpace
implements a Cartesian grid with logarithmically spaced bins and an arbitrary number of dimensions.HPSpace
andGLSpace
describe pixelisations of the 2-sphere; their counterpart in harmonic space isLMSpace
, which contains spherical harmonic coefficients.PowerSpace
is used to describe one-dimensional power spectra.
Among these, RGSpace
and LogRGSpace
can
be harmonic or not (depending on constructor arguments),
GLSpace
, HPSpace
,
and PowerSpace
are pure position domains (i.e.
nonharmonic), and LMSpace
is always harmonic.
Combinations of domains#
The fundamental classes described above are often sufficient to specify the domain of a field. In some cases, however, it will be necessary to define the field on a product of elementary domains instead of a single one. More sophisticated operators also require a set of several such fields. Some examples are:
sky emission depending on location and energy. This could be represented by a product of an
HPSpace
(for location) with anRGSpace
(for energy).a polarized field, which could be modelled as a product of any structured domain (representing location) with a four-element
UnstructuredDomain
holding Stokes I, Q, U and V components.a model for the sky emission, which holds both the current realisation (on a harmonic domain) and a few inferred model parameters (e.g. on an unstructured grid).
Consequently, NIFTy defines a class called DomainTuple
holding a sequence of Domain
objects. The full domain is
specified as the product of all elementary domains. Thus, an instance of
DomainTuple
would be suitable to describe the first two
examples above. In principle, a DomainTuple
can even be empty, which implies that the field living on it is a scalar.
A DomainTuple
supports iteration and indexing, and also
provides the properties shape
and
size
in analogy to the elementary
Domain
.
An aggregation of several DomainTuple
s, each member
identified by a name, is described by the MultiDomain
class. In contrast to a DomainTuple
a
MultiDomain
is a collection and does not define the
product space of its elements. It would be the adequate space to use in the
last of the above examples.
Fields#
Fields on a single DomainTuple#
A Field
object consists of the following components:
a domain in form of a
DomainTuple
objecta data type (e.g. numpy.float64)
an array containing the actual values
Usually, the array is stored in the form of a numpy.ndarray
, but for very
resource-intensive tasks NIFTy also provides an alternative storage method to
be used with distributed memory processing.
Fields support a wide range of arithmetic operations, either involving two fields of equal domains or a field and a scalar. Arithmetic operations are performed point-wise, and the returned field has the same domain as the input field(s). Available operators are addition (“+”), subtraction (“-“), multiplication (“*”), division (“/”), floor division (“//”) and exponentiation (“**”). Inplace operators (“+=”, etc.) are not supported. Further, boolean operators, performing a point-wise comparison of a field with either another field of equal domain or a scalar, are available as well. These include equals (“==”), not equals (“!=”), less (“<”), less or equal (“<=”), greater (“>”) and greater or equal (“>=). The domain of the field returned equals that of the input field(s), while the stored data is of boolean type.
Contractions (like summation, integration, minimum/maximum, computation of
statistical moments) can be carried out either over an entire field (producing
a scalar result) or over sub-domains (resulting in a field defined on a smaller
domain). Scalar products of two fields can also be computed easily as well.
See the documentation of Field
for details.
There is also a set of convenience functions to generate fields with constant
values or fields filled with random numbers according to a user-specified
distribution: full
, from_random
.
Like almost all NIFTy objects, fields are immutable: their value or any other attribute cannot be modified after construction. To manipulate a field in ways that are not covered by the provided standard operations, its data content must be extracted first, then changed, and a new field has to be created from the result.
Fields defined on a MultiDomain#
The MultiField
class can be seen as a dictionary of
individual Field
s, each identified by a name, which is defined
on a MultiDomain
.
Operators#
All transformations between different NIFTy fields are expressed in the form of
Operator
objects. The interface of this class is
rather minimalistic: it has a property called
domain
which returns a
DomainTuple
or MultiDomain
object
specifying the structure of the Field
or
MultiField
it expects as input, another property
target
describing its output, and finally
an overloaded apply
method, which can take:
a
Field
/MultiField
object, in which case it returns the transformedField
/MultiField
.a
Linearization
object, in which case it returns the transformedLinearization
.
This is the interface that all objects derived from
Operator
must implement. In addition,
Operator
objects can be added/subtracted,
multiplied, chained (via the __call__
method or the @ operator) and
support point-wise application of functions like exp()
, log()
,
sqrt()
, conjugate()
.
Advanced operators#
NIFTy provides a library of commonly employed operators which can be used for specific inference problems. Currently these are:
SLAmplitude
, which returns a smooth power spectrum.InverseGammaOperator
, which models point sources which are distributed according to a inverse-gamma distribution.CorrelatedField
, which models a diffuse field whose correlation structure is described by an amplitude operator.
Linear Operators#
A linear operator (represented by NIFTy8’s abstract
LinearOperator
class) is derived from
Operator
and can be interpreted as an (implicitly defined)
matrix. Since its operation is linear, it can provide some additional
functionality which is not available for the more generic
Operator
class.
Linear Operator basics#
There are four basic ways of applying an operator to a field :
direct application:
adjoint application:
inverse application:
adjoint inverse application:
Note: The inverse of the adjoint of a linear map and the adjoint of the inverse of a linear map (if all those exist) are the same.
These different actions of a linear operator Op
on a field f
can be
invoked in various ways:
direct multiplication:
Op(f)
orOp.times(f)
orOp.apply(f, Op.TIMES)
adjoint multiplication:
Op.adjoint_times(f)
orOp.apply(f, Op.ADJOINT_TIMES)
inverse multiplication:
Op.inverse_times(f)
orOp.apply(f, Op.INVERSE_TIMES)
adjoint inverse multiplication:
Op.adjoint_inverse_times(f)
orOp.apply(f, Op.ADJOINT_INVERSE_TIMES)
Operator classes defined in NIFTy may implement an arbitrary subset of these
four operations. This subset can be queried using the
capability
property.
If needed, the set of supported operations can be enhanced by iterative
inversion methods; for example, an operator defining direct and adjoint
multiplication could be enhanced by this approach to support the complete set.
This functionality is provided by NIFTy’s
InversionEnabler
class, which is itself a linear
operator.
Direct multiplication and adjoint inverse multiplication transform a field
defined on the operator’s domain
to one defined on the
operator’s target
, whereas adjoint multiplication and inverse
multiplication transform from target
to
domain
.
Operators with identical domain and target can be derived from
EndomorphicOperator
. Typical examples for this
category are the ScalingOperator
, which simply
multiplies its input by a scalar value, and
DiagonalOperator
, which multiplies every value of
its input field with potentially different values.
Further operator classes provided by NIFTy are
HarmonicTransformOperator
for transforms from a harmonic domain to its counterpart in position space, and their adjointPowerDistributor
for transforms from aPowerSpace
to an associated harmonic domain, and their adjoint.GeometryRemover
, which transforms from structured domains to unstructured ones. This is typically needed when building instrument response operators.
Syntactic sugar#
NIFTy allows simple and intuitive construction of altered and combined
operators.
As an example, if A
, B
and C
are of type LinearOperator
and f1
and f2
are of type Field
, writing:
X = A(B.inverse(A.adjoint)) + C
f2 = X(f1)
will perform the operation suggested intuitively by the notation, checking
domain compatibility while building the composed operator.
The properties adjoint
and
inverse
return a new operator which behaves as if it
were the original operator’s adjoint or inverse, respectively.
The combined operator infers its domain and target from its constituents,
as well as the set of operations it can support.
Instantiating operator adjoints or inverses by adjoint
and similar methods is to be distinguished from the instant application of
operators performed by adjoint_times
and similar
methods.
Minimization#
Most problems in IFT are solved by (possibly nested) minimizations of high-dimensional functions, which are often nonlinear.
Energy functionals#
In NIFTy8 such functions are represented by objects of type
Energy
. These hold the prescription how to calculate the
function’s value
, gradient
and
(optionally) metric
at any given
position
in parameter space. Function values are
floating-point scalars, gradients have the form of fields defined on the energy’s
position domain, and metrics are represented by linear operator objects.
Energies are classes that typically have to be provided by the user when
tackling new IFT problems. An example of concrete energy classes delivered with
NIFTy8 is QuadraticEnergy
(with
position-independent metric, mainly used with conjugate gradient minimization).
For MGVI and GeoVI, NIFTy provides
SampledKLEnergy()
that instantiate objects
containing the sampled estimate of the KL divergence, its gradient and the
Fisher metric. These constructors require an instance
of StandardHamiltonian
, an operator to
compute the negative log-likelihood of the problem in standardized coordinates
at a given position in parameter space.
Finally, the StandardHamiltonian
can be constructed from the likelihood, represented by a
LikelihoodEnergyOperator
instance.
Several commonly used forms of the likelihoods are already provided in
NIFTy, such as GaussianEnergy
,
PoissonianEnergy
,
InverseGammaEnergy
or
BernoulliEnergy
, but the user
is free to implement any likelihood customized to the problem at hand.
The demo code demos/getting_started_3.py illustrates how to set up an energy
functional for MGVI and minimize it.
Iteration control#
Iterative minimization of an energy requires some means of checking the quality of the current solution estimate and stopping once it is sufficiently accurate. In case of numerical problems, the iteration needs to be terminated as well, returning a suitable error description.
In NIFTy8, this functionality is encapsulated in the abstract
IterationController
class, which is provided with the initial energy
object before starting the minimization, and is updated with the improved
energy after every iteration. Based on this information, it can either continue
the minimization or return the current estimate indicating convergence or
failure.
Sensible stopping criteria can vary significantly with the problem being
solved; NIFTy provides a concrete sub-class of IterationController
called GradientNormController
, which should be appropriate in many
circumstances. A full list of the available IterationController
s
in NIFTy can be found below, but users have complete freedom to implement custom
IterationController
sub-classes for their specific applications.
Minimization algorithms#
All minimization algorithms in NIFTy inherit from the abstract
Minimizer
class, which presents a minimalistic interface
consisting only of a __call__()
method taking an
Energy
object and optionally a preconditioning operator, and
returning the energy at the discovered minimum and a status code.
For energies with a quadratic form (i.e. which can be expressed by means of a
QuadraticEnergy
object), an obvious choice of
algorithm is the ConjugateGradient
minimizer.
A similar algorithm suited for nonlinear problems is provided by
NonlinearCG
.
Many minimizers for nonlinear problems can be characterized as
First deciding on a direction for the next step.
Then finding a suitable step length along this direction, resulting in the next energy estimate.
This family of algorithms is encapsulated in NIFTy’s
DescentMinimizer
class, which currently has three
generally usable concrete implementations:
NewtonCG
, L_BFGS
and
VL_BFGS
. Of these algorithms, only
NewtonCG
requires the energy object to provide
a metric
property, the others only need energy values and
gradients. Further available descent minimizers are
RelaxedNewton
and SteepestDescent
.
The flexibility of NIFTy’s design allows using externally provided minimizers.
With only small effort, adaptors for two SciPy minimizers were written; they are
available under the names ScipyCG
and
L_BFGS_B
.
Application to operator inversion#
The machinery presented here cannot only be used for minimizing functionals derived from IFT, but also for the numerical inversion of linear operators, if the desired application mode is not directly available. A classical example is the information propagator whose inverse is defined as:
.
It needs to be applied in forward direction in order to calculate the Wiener
filter solution, but only its inverse application is straightforward.
To use it in forward direction, we make use of NIFTy’s
InversionEnabler
class, which internally
applies the (approximate) inverse of the given operator by
solving the equation for .
This is accomplished by minimizing a suitable
QuadraticEnergy
with the ConjugateGradient
algorithm. An example is provided in
WienerFilterCurvature()
.
Posterior analysis and visualization#
After the minimization of an energy functional has converged, samples can be drawn
from the posterior distribution at the current position to investigate the result.
The probing module offers class called StatCalculator
which allows to evaluate the mean
and the unbiased
variance var
of these samples.
Fields can be visualized using the Plot
class, which invokes
matplotlib for plotting.