Commit ec2f007a authored by Philipp Arras's avatar Philipp Arras

Operators -> Models

parent 4c0bd584
......@@ -53,7 +53,7 @@ if __name__ == '__main__':
A = ift.create_power_operator(harmonic_space, sqrtpspec)
# Set up a sky model and instrumental response
# Set up a sky operator and instrumental response
sky = ift.positive_tanh(HT(A))
GR = ift.GeometryRemover(position_space)
R = GR
......
......@@ -76,7 +76,7 @@ if __name__ == '__main__':
a = ift.PS_field(p_space, sqrtpspec)
A = pd(a)
# Define sky model
# Define sky operator
sky = ift.exp(HT(ift.makeOp(A)))
M = ift.DiagonalOperator(exposure)
......
......@@ -47,7 +47,7 @@ if __name__ == '__main__':
position_space = ift.RGSpace([128, 128])
# Set up an amplitude model for the field
# Set up an amplitude operator for the field
# The parameters mean:
# 64 spectral bins
#
......@@ -60,9 +60,9 @@ if __name__ == '__main__':
# 0.5 = low variance of power-law slope
# 0.4 = y-intercept mean
# 0.3 = relatively high y-intercept variance
A = ift.AmplitudeModel(position_space, 64, 3, 0.4, -5., 0.5, 0.4, 0.3)
A = ift.AmplitudeOperator(position_space, 64, 3, 0.4, -5., 0.5, 0.4, 0.3)
# Build the model for a correlated signal
# Build the operator for a correlated signal
harmonic_space = position_space.get_default_codomain()
ht = ift.HarmonicTransformOperator(harmonic_space, position_space)
power_space = A.target[0]
......@@ -98,7 +98,7 @@ if __name__ == '__main__':
name='Newton', tol=1e-7, iteration_limit=35)
minimizer = ift.NewtonCG(ic_newton)
# Set up model likelihood and information Hamiltonian
# Set up likelihood and information Hamiltonian
likelihood = ift.GaussianEnergy(mean=data, covariance=N)(signal_response)
H = ift.Hamiltonian(likelihood, ic_sampling)
......
......@@ -15,7 +15,7 @@ From such a perspective,
- IFT problems largely consist of the combination of several high dimensional
*minimization* problems.
- Within NIFTy, *models* are used to define the characteristic equations and
- 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.
......@@ -99,13 +99,13 @@ 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 models also require a set of several such fields.
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 :class:`HPSpace` (for location) with an :class:`RGSpace`
(for energy).
- a polarised field, which could be modeled as a product of any structured
- a polarized field, which could be modeled as a product of any structured
domain (representing location) with a four-element
:class:`UnstructuredDomain` holding Stokes I, Q, U and V components.
- a model for the sky emission, which holds both the current realization
......@@ -269,59 +269,59 @@ The properties :attr:`~LinearOperator.adjoint` and
were the original operator's adjoint or inverse, respectively.
Models
======
Operators
=========
Model classes (represented by NIFTy5's abstract :class:`Model` class) are used to construct
Operator classes (represented by NIFTy5's abstract :class:`Operator` class) are used to construct
the equations of a specific inference problem.
Most models are defined via a position, which is a :class:`MultiField` object,
Most operators are defined via a position, which is a :class:`MultiField` object,
their value at this position, which is again a :class:`MultiField` object and a Jacobian derivative,
which is a :class:`LinearOperator` and is needed for the minimization procedure.
Using the existing basic model classes one can construct more complicated models, as
Using the existing basic operator classes one can construct more complicated operators, as
NIFTy allows for easy and self-consinstent combination via point-wise multiplication,
addition and subtraction. The model resulting from these operations then automatically
addition and subtraction. The operator resulting from these operations then automatically
contains the correct Jacobians, positions and values.
Notably, :class:`Constant` and :class:`Variable` allow for an easy way to turn
inference of specific quantities on and off.
The basic model classes also allow for more complex operations on models such as
The basic operator classes also allow for more complex operations on operators such as
the application of :class:`LinearOperators` or local non-linearities.
As an example one may consider the following combination of ``x``, which is a model of type
:class:`Variable` and ``y``, which is a model of type :class:`Constant`::
As an example one may consider the following combination of ``x``, which is an operator of type
:class:`Variable` and ``y``, which is an operator of type :class:`Constant`::
z = x*x + y
``z`` will then be a model with the following properties::
``z`` will then be an operator with the following properties::
z.value = x.value*x.value + y.value
z.position = Union(x.position, y.position)
z.jacobian = 2*makeOp(x.value)
Basic models
Basic operators
------------
# FIXME All this is outdated!
Basic model classes provided by NIFTy are
Basic operator classes provided by NIFTy are
- :class:`Constant` contains a constant value and has a zero valued Jacobian.
Like other models, it has a position, but its value does not depend on it.
Like other operators, it has a position, but its value does not depend on it.
- :class:`Variable` returns the position as its value, its derivative is one.
- :class:`LinearModel` applies a :class:`LinearOperator` on the model.
- :class:`LocalModel` applies a non-linearity locally on the model.
- :class:`MultiModel` combines various models into one. In this case the position,
value and Jacobian are combined into corresponding :class:`MultiFields` and operators.
Advanced models
---------------
Advanced operators
------------------
NIFTy also provides a library of more sophisticated models which are used for more
NIFTy also provides a library of more sophisticated operators which are used for more
specific inference problems. Currently these are:
- :class:`AmplitudeModel`, which returns a smooth power spectrum.
- :class:`PointModel`, which models point sources which follow a inverse gamma distribution.
- :class:`SmoothSkyModel`, which models a diffuse lognormal field. It takes an amplitude model
- :class:`AmplitudeOperator`, which returns a smooth power spectrum.
- :class:`InverseGammaOperator`, which models point sources which follow a inverse gamma distribution.
- :class:`CorrelatedField`, which models a diffuse log-normal field. It takes an amplitude operator
to specify the correlation structure of the field.
......
......@@ -174,7 +174,7 @@ NIFTy takes advantage of this formulation in several ways:
1) All prior degrees of freedom have unit covariance which improves the condition number of operators which need to be inverted.
2) The amplitude operator can be regarded as part of the response, :math:`{R'=R\,A}`. In general, more sophisticated responses can be constructed out of the composition of simpler operators.
3) The response can be non-linear, e.g. :math:`{R'(s)=R \exp(A\,\xi)}`, see demos/getting_started_2.py.
4) The amplitude operator can be made dependent on unknowns as well, e.g. :math:`A=A(\tau)= F\, \widehat{e^\tau}` represents an amplitude model with a positive definite, unknown spectrum defined in Fourier domain. The amplitude field :math:`{\tau}` would get its own amplitude model, with a cepstrum (spectrum of a log spectrum) defined in quefrency space (harmonic space of a logarithmically binned harmonic space) to regularize its degrees of freedom by imposing some (user-defined level of) spectral smoothness.
4) The amplitude operator can be made dependent on unknowns as well, e.g. :math:`A=A(\tau)= F\, \widehat{e^\tau}` represents an amplitude operator with a positive definite, unknown spectrum defined in Fourier domain. The amplitude field :math:`{\tau}` would get its own amplitude operator, with a cepstrum (spectrum of a log spectrum) defined in quefrency space (harmonic space of a logarithmically binned harmonic space) to regularize its degrees of freedom by imposing some (user-defined level of) spectral smoothness.
5) NIFTy can calculate the gradient of the information Hamiltonian and the Fischer information metric with respect to all unknown parameters, here :math:`{\xi}` and :math:`{\tau}`, by automatic differentiation. The gradients are used for MAP and HMCF estimates, and the Fischer matrix is required in addition to the gradient by Metric Gaussian Variational Inference (MGVI), which is available in NIFTy as well. MGVI is an implicit operator extension of Automatic Differentiation Variational Inference (ADVI).
The reconstruction of a non-Gaussian signal with unknown covarinance from a non-trivial (tomographic) response is demonstrated in demos/getting_started_3.py. Here, the uncertainty of the field and the power spectrum of its generating process are probed via posterior samples provided by the MGVI algorithm.
......
......@@ -73,8 +73,8 @@ from .minimization.kl_energy import KL_Energy
from .sugar import *
from .plot import Plot
from .library.amplitude_model import AmplitudeModel
from .library.inverse_gamma_model import InverseGammaModel
from .library.amplitude_operator import AmplitudeOperator
from .library.inverse_gamma_operator import InverseGammaOperator
from .library.los_response import LOSResponse
from .library.wiener_filter_curvature import WienerFilterCurvature
......
......@@ -75,14 +75,14 @@ def make_adjust_variances(a,
def do_adjust_variances(position,
amplitude_model,
amplitude_operator,
minimizer,
xi_key='xi',
samples=[]):
h_space = position[xi_key].domain[0]
pd = PowerDistributor(h_space, amplitude_model.target[0])
a = pd(amplitude_model)
pd = PowerDistributor(h_space, amplitude_operator.target[0])
a = pd(amplitude_operator)
xi = ducktape(None, position.domain, xi_key)
ham = make_adjust_variances(a, xi, position, samples=samples)
......
......@@ -62,7 +62,7 @@ def create_cepstrum_amplitude_field(domain, cepstrum):
return Field.from_global_data(domain, cepstrum_field)
def CepstrumOperator(logk_space, ceps_a, ceps_k, zero_mode=True):
def _CepstrumOperator(logk_space, ceps_a, ceps_k, zero_mode=True):
'''
Parameters
----------
......@@ -88,7 +88,7 @@ def CepstrumOperator(logk_space, ceps_a, ceps_k, zero_mode=True):
return res
def SlopeModel(logk_space, sm, sv, im, iv):
def _SlopePowerSpectrum(logk_space, sm, sv, im, iv):
'''
Parameters
----------
......@@ -109,8 +109,8 @@ def SlopeModel(logk_space, sm, sv, im, iv):
return slope(OffsetOperator(phi_mean)(makeOp(phi_sig)))
def AmplitudeModel(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
keys=['tau', 'phi'], zero_mode=True):
def AmplitudeOperator(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
keys=['tau', 'phi'], zero_mode=True):
'''
Computes a smooth power spectrum.
Output is defined on a PowerSpace.
......@@ -137,9 +137,9 @@ def AmplitudeModel(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
et = ExpTransform(PowerSpace(h_space), Npixdof)
logk_space = et.domain[0]
smooth = CepstrumOperator(logk_space, ceps_a, ceps_k, zero_mode)
smooth = _CepstrumOperator(logk_space, ceps_a, ceps_k, zero_mode)
smooth = smooth.ducktape(keys[0])
linear = SlopeModel(logk_space, sm, sv, im, iv)
linear = _SlopePowerSpectrum(logk_space, sm, sv, im, iv)
linear = linear.ducktape(keys[1])
fac = ScalingOperator(0.5, smooth.target)
......
......@@ -16,7 +16,6 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from ..domain_tuple import DomainTuple
from ..multi_domain import MultiDomain
from ..operators.contraction_operator import ContractionOperator
from ..operators.distributors import PowerDistributor
from ..operators.harmonic_operators import HarmonicTransformOperator
......@@ -24,7 +23,7 @@ from ..operators.simple_linear_operators import ducktape
from ..operators.scaling_operator import ScalingOperator
def CorrelatedField(s_space, amplitude_model, name='xi'):
def CorrelatedField(s_space, amplitude_operator, name='xi'):
'''
Function for construction of correlated fields
......@@ -32,23 +31,26 @@ def CorrelatedField(s_space, amplitude_model, name='xi'):
----------
s_space : Domain
Field domain
amplitude_model: Operator
model for correlation structure
amplitude_operator: Operator
operator for correlation structure
name : string
MultiField component name
'''
h_space = s_space.get_default_codomain()
ht = HarmonicTransformOperator(h_space, s_space)
p_space = amplitude_model.target[0]
p_space = amplitude_operator.target[0]
power_distributor = PowerDistributor(h_space, p_space)
A = power_distributor(amplitude_model)
A = power_distributor(amplitude_operator)
vol = h_space.scalar_dvol
vol = ScalingOperator(vol**(-0.5), h_space)
return ht(vol(A)*ducktape(h_space, None, name))
def MfCorrelatedField(s_space_spatial, s_space_energy, amplitude_model_spatial,
amplitude_model_energy, name="xi"):
def MfCorrelatedField(s_space_spatial,
s_space_energy,
amplitude_operator_spatial,
amplitude_operator_energy,
name="xi"):
'''
Method for construction of correlated multi-frequency fields
'''
......@@ -59,8 +61,8 @@ def MfCorrelatedField(s_space_spatial, s_space_energy, amplitude_model_spatial,
ht2 = HarmonicTransformOperator(ht1.target, space=1)
ht = ht2(ht1)
p_space_spatial = amplitude_model_spatial.target[0]
p_space_energy = amplitude_model_energy.target[0]
p_space_spatial = amplitude_operator_spatial.target[0]
p_space_energy = amplitude_operator_energy.target[0]
pd_spatial = PowerDistributor(h_space, p_space_spatial, 0)
pd_energy = PowerDistributor(pd_spatial.domain, p_space_energy, 1)
......@@ -69,8 +71,8 @@ def MfCorrelatedField(s_space_spatial, s_space_energy, amplitude_model_spatial,
dom_distr_spatial = ContractionOperator(pd.domain, 1).adjoint
dom_distr_energy = ContractionOperator(pd.domain, 0).adjoint
a_spatial = dom_distr_spatial(amplitude_model_spatial)
a_energy = dom_distr_energy(amplitude_model_energy)
a_spatial = dom_distr_spatial(amplitude_operator_spatial)
a_energy = dom_distr_energy(amplitude_operator_energy)
a = a_spatial*a_energy
A = pd(a)
return ht(A*ducktape(h_space, None, name))
......@@ -25,9 +25,9 @@ from ..operators.operator import Operator
from ..sugar import makeOp
class InverseGammaModel(Operator):
class InverseGammaOperator(Operator):
def __init__(self, domain, alpha, q, delta=0.001):
"""Model which transforms a Gaussian into an inverse gamma distribution.
"""Operator which transforms a Gaussian into an inverse gamma distribution.
The pdf of the inverse gamma distribution is defined as follows:
......
......@@ -24,7 +24,7 @@ import numpy as np
class Energy_Tests(unittest.TestCase):
def make_model(self, **kwargs):
def make_operator(self, **kwargs):
np.random.seed(kwargs['seed'])
S = ift.ScalingOperator(1., kwargs['space'])
s = S.draw_sample()
......@@ -37,10 +37,10 @@ class Energy_Tests(unittest.TestCase):
[4, 78, 23]
))
def testGaussian(self, space, seed):
model = self.make_model(
op = self.make_operator(
space_key='s1', space=space, seed=seed)['s1']
energy = ift.GaussianEnergy(domain=space)
ift.extra.check_value_gradient_consistency(energy, model)
ift.extra.check_value_gradient_consistency(energy, op)
# @expand(product(
# [ift.GLSpace(15),
......@@ -62,12 +62,12 @@ class Energy_Tests(unittest.TestCase):
ift.RGSpace([32, 32], distances=.789)
], [4, 78, 23]))
def testInverseGammaLikelihood(self, space, seed):
model = self.make_model(space_key='s1', space=space, seed=seed)['s1']
model = model.exp()
op = self.make_operator(space_key='s1', space=space, seed=seed)['s1']
op = op.exp()
d = np.random.normal(10, size=space.shape)**2
d = ift.Field.from_global_data(space, d)
energy = ift.InverseGammaLikelihood(d)
ift.extra.check_value_gradient_consistency(energy, model, tol=1e-7)
ift.extra.check_value_gradient_consistency(energy, op, tol=1e-7)
@expand(product(
[ift.GLSpace(15),
......@@ -76,13 +76,13 @@ class Energy_Tests(unittest.TestCase):
[4, 78, 23]
))
def testPoissonian(self, space, seed):
model = self.make_model(
op = self.make_operator(
space_key='s1', space=space, seed=seed)['s1']
model = model.exp()
op = op.exp()
d = np.random.poisson(120, size=space.shape)
d = ift.Field.from_global_data(space, d)
energy = ift.PoissonianEnergy(d)
ift.extra.check_value_gradient_consistency(energy, model, tol=1e-7)
ift.extra.check_value_gradient_consistency(energy, op, tol=1e-7)
@expand(product(
[ift.GLSpace(15),
......@@ -91,16 +91,16 @@ class Energy_Tests(unittest.TestCase):
[4, 78, 23]
))
def testHamiltonian_and_KL(self, space, seed):
model = self.make_model(
op = self.make_operator(
space_key='s1', space=space, seed=seed)['s1']
model = model.exp()
op = op.exp()
lh = ift.GaussianEnergy(domain=space)
hamiltonian = ift.Hamiltonian(lh)
ift.extra.check_value_gradient_consistency(hamiltonian, model)
ift.extra.check_value_gradient_consistency(hamiltonian, op)
S = ift.ScalingOperator(1., space)
samps = [S.draw_sample() for i in range(3)]
kl = ift.SampledKullbachLeiblerDivergence(hamiltonian, samps)
ift.extra.check_value_gradient_consistency(kl, model)
ift.extra.check_value_gradient_consistency(kl, op)
@expand(product(
[ift.GLSpace(15),
......@@ -109,10 +109,10 @@ class Energy_Tests(unittest.TestCase):
[4, 78, 23]
))
def testBernoulli(self, space, seed):
model = self.make_model(
op = self.make_operator(
space_key='s1', space=space, seed=seed)['s1']
model = model.positive_tanh()
op = op.positive_tanh()
d = np.random.binomial(1, 0.1, size=space.shape)
d = ift.Field.from_global_data(space, d)
energy = ift.BernoulliEnergy(d)
ift.extra.check_value_gradient_consistency(energy, model, tol=1e-6)
ift.extra.check_value_gradient_consistency(energy, op, tol=1e-6)
......@@ -23,7 +23,7 @@ import nifty5 as ift
import numpy as np
class Model_Tests(unittest.TestCase):
class OperatorTests(unittest.TestCase):
@staticmethod
def make_linearization(type, space, seed):
np.random.seed(seed)
......@@ -43,8 +43,8 @@ class Model_Tests(unittest.TestCase):
))
def testBasics(self, space, seed):
var = self.make_linearization("Variable", space, seed)
model = ift.ScalingOperator(6., var.target)
ift.extra.check_value_gradient_consistency(model, var.val)
op = ift.ScalingOperator(6., var.target)
ift.extra.check_value_gradient_consistency(op, var.val)
@expand(product(
['Variable', 'Constant'],
......@@ -63,29 +63,29 @@ class Model_Tests(unittest.TestCase):
dom = ift.MultiDomain.union((dom1, dom2))
select_s1 = ift.ducktape(None, dom, "s1")
select_s2 = ift.ducktape(None, dom, "s2")
model = select_s1*select_s2
op = select_s1*select_s2
pos = ift.from_random("normal", dom)
ift.extra.check_value_gradient_consistency(model, pos, ntries=20)
model = select_s1+select_s2
ift.extra.check_value_gradient_consistency(op, pos, ntries=20)
op = select_s1+select_s2
pos = ift.from_random("normal", dom)
ift.extra.check_value_gradient_consistency(model, pos, ntries=20)
model = select_s1.scale(3.)
ift.extra.check_value_gradient_consistency(op, pos, ntries=20)
op = select_s1.scale(3.)
pos = ift.from_random("normal", dom1)
ift.extra.check_value_gradient_consistency(model, pos, ntries=20)
model = ift.ScalingOperator(2.456, space)(select_s1*select_s2)
ift.extra.check_value_gradient_consistency(op, pos, ntries=20)
op = ift.ScalingOperator(2.456, space)(select_s1*select_s2)
pos = ift.from_random("normal", dom)
ift.extra.check_value_gradient_consistency(model, pos, ntries=20)
model = ift.positive_tanh(ift.ScalingOperator(2.456, space)(
ift.extra.check_value_gradient_consistency(op, pos, ntries=20)
op = ift.positive_tanh(ift.ScalingOperator(2.456, space)(
select_s1*select_s2))
pos = ift.from_random("normal", dom)
ift.extra.check_value_gradient_consistency(model, pos, ntries=20)
ift.extra.check_value_gradient_consistency(op, pos, ntries=20)
pos = ift.from_random("normal", dom)
model = ift.OuterProduct(pos['s1'], ift.makeDomain(space))
ift.extra.check_value_gradient_consistency(model, pos['s2'], ntries=20)
op = ift.OuterProduct(pos['s1'], ift.makeDomain(space))
ift.extra.check_value_gradient_consistency(op, pos['s2'], ntries=20)
if isinstance(space, ift.RGSpace):
model = ift.FFTOperator(space)(select_s1*select_s2)
op = ift.FFTOperator(space)(select_s1*select_s2)
pos = ift.from_random("normal", dom)
ift.extra.check_value_gradient_consistency(model, pos, ntries=20)
ift.extra.check_value_gradient_consistency(op, pos, ntries=20)
@expand(product(
[ift.GLSpace(15),
......@@ -100,45 +100,32 @@ class Model_Tests(unittest.TestCase):
[1.3],
[4, 78, 23],
))
def testModelLibrary(self, space, Npixdof, ceps_a,
ceps_k, sm, sv, im, iv, seed):
# tests amplitude model and coorelated field model
def testOperatorLibrary(self, space, Npixdof, ceps_a,
ceps_k, sm, sv, im, iv, seed):
# tests amplitude operator and coorelated field operator
np.random.seed(seed)
model = ift.AmplitudeModel(space, Npixdof, ceps_a, ceps_k, sm,
op = ift.AmplitudeOperator(space, Npixdof, ceps_a, ceps_k, sm,
sv, im, iv)
S = ift.ScalingOperator(1., model.domain)
S = ift.ScalingOperator(1., op.domain)
pos = S.draw_sample()
ift.extra.check_value_gradient_consistency(model, pos, ntries=20)
ift.extra.check_value_gradient_consistency(op, pos, ntries=20)
model2 = ift.CorrelatedField(space, model)
S = ift.ScalingOperator(1., model2.domain)
op2 = ift.CorrelatedField(space, op)
S = ift.ScalingOperator(1., op2.domain)
pos = S.draw_sample()
ift.extra.check_value_gradient_consistency(model2, pos, ntries=20)
ift.extra.check_value_gradient_consistency(op2, pos, ntries=20)
@expand(product(
[ift.GLSpace(15),
ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)],
[4, 78, 23]))
def testPointModel(self, space, seed):
def testInvGammaOperator(self, space, seed):
S = ift.ScalingOperator(1., space)
pos = S.draw_sample()
alpha = 1.5
q = 0.73
model = ift.InverseGammaModel(space, alpha, q)
op = ift.InverseGammaOperator(space, alpha, q)
# FIXME All those cdfs and ppfs are not very accurate
ift.extra.check_value_gradient_consistency(model, pos, tol=1e-2,
ift.extra.check_value_gradient_consistency(op, pos, tol=1e-2,
ntries=20)
# @expand(product(
# ['Variable', 'Constant'],
# [ift.GLSpace(15),
# ift.RGSpace(64, distances=.789),
# ift.RGSpace([32, 32], distances=.789)],
# [4, 78, 23]
# ))
# def testMultiModel(self, type, space, seed):
# model = self.make_model(
# type, space_key='s', space=space, seed=seed)['s']
# mmodel = ift.MultiModel(model, 'g')
# ift.extra.check_value_gradient_consistency(mmodel)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment