diff --git a/demos/getting_started_3.py b/demos/getting_started_3.py
index e088f427ea8fdd6e2d34d815055442311a56f944..9ca2994970f3c825f232f3457c9ce2216a9f0a11 100644
--- a/demos/getting_started_3.py
+++ b/demos/getting_started_3.py
@@ -65,7 +65,7 @@ if __name__ == '__main__':
'im': .4, # y-intercept mean
'iv': .3 # relatively high y-intercept variance
}
- A = ift.AmplitudeOperator(**dct)
+ A = ift.SLAmplitude(**dct)
# Build the operator for a correlated signal
power_distributor = ift.PowerDistributor(harmonic_space, power_space)
diff --git a/docs/source/code.rst b/docs/source/code.rst
index cf462ced86ad79a8bc38c8c2ed11849855c39cb6..06bc84a559f6e0f5a8297a2f2ab5e924a7718435 100644
--- a/docs/source/code.rst
+++ b/docs/source/code.rst
@@ -212,7 +212,7 @@ specific inference problems. Currently these are:
.. currentmodule:: nifty5.library
-- :class:`~amplitude_operator.AmplitudeOperator`, which returns a smooth power spectrum.
+- :class:`~smooth_linear_amplitude.SLAmplitude`, which returns a smooth power spectrum.
- :class:`~inverse_gamma_operator.InverseGammaOperator`, which models point sources which are
distributed according to a inverse-gamma distribution.
- :class:`~correlated_fields.CorrelatedField`, which models a diffuse log-normal field. It takes an
diff --git a/nifty5/__init__.py b/nifty5/__init__.py
index 30f7d5b41c647e481537150f198a3f8bd67f4978..026e69f31f656a4a1b3289eccbefd1c7378c1ed1 100644
--- a/nifty5/__init__.py
+++ b/nifty5/__init__.py
@@ -73,7 +73,7 @@ from .minimization.kl_energy import KL_Energy
from .sugar import *
from .plot import Plot
-from .library.amplitude_operator import AmplitudeOperator
+from .library.smooth_linear_amplitude import (SLAmplitude, CepstrumOperator)
from .library.inverse_gamma_operator import InverseGammaOperator
from .library.los_response import LOSResponse
from .library.dynamic_operator import (dynamic_operator,
diff --git a/nifty5/data_objects/numpy_do.py b/nifty5/data_objects/numpy_do.py
index 8160f3b34163283a146c977cdd14d4aa36f0d370..e582e36023bb914a0195d5e6deb9379a62a68548 100644
--- a/nifty5/data_objects/numpy_do.py
+++ b/nifty5/data_objects/numpy_do.py
@@ -15,7 +15,7 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
-# Data object module for NIFTy that uses simple numpy ndarrays.
+# Data object module that uses simple numpy ndarrays.
import numpy as np
from numpy import absolute, clip, cos, cosh, empty, empty_like, exp, full, log
diff --git a/nifty5/domains/gl_space.py b/nifty5/domains/gl_space.py
index 7f6a48d97fae0686e91614a1e562aca53dcd08d3..d6091d5388129f14515d1f5c2e0ce4659d688e25 100644
--- a/nifty5/domains/gl_space.py
+++ b/nifty5/domains/gl_space.py
@@ -21,7 +21,7 @@ from .structured_domain import StructuredDomain
class GLSpace(StructuredDomain):
- """NIFTy subclass for Gauss-Legendre pixelizations of the two-sphere.
+ """Represents a 2-sphere with Gauss-Legendre pixelization.
Its harmonic partner domain is the
:class:`~nifty5.domains.lm_space.LMSpace`.
diff --git a/nifty5/domains/hp_space.py b/nifty5/domains/hp_space.py
index 998fa7ac5548a170cd1f5a252817ab252c2a265e..6058647330d67e56de7aa0a212c3387a4edf0eda 100644
--- a/nifty5/domains/hp_space.py
+++ b/nifty5/domains/hp_space.py
@@ -21,7 +21,7 @@ from .structured_domain import StructuredDomain
class HPSpace(StructuredDomain):
- """NIFTy subclass for HEALPix discretizations of the two-sphere.
+ """Represents 2-sphere with HEALPix discretization.
Its harmonic partner domain is the
:class:`~nifty5.domains.lm_space.LMSpace`.
diff --git a/nifty5/domains/lm_space.py b/nifty5/domains/lm_space.py
index 80223a63053d1def496082a8387c0431c90ab9a5..4c224a6dec96764f18bc27019a1f88b2d33630b8 100644
--- a/nifty5/domains/lm_space.py
+++ b/nifty5/domains/lm_space.py
@@ -22,7 +22,7 @@ from .structured_domain import StructuredDomain
class LMSpace(StructuredDomain):
- """NIFTy subclass for sets of spherical harmonic coefficients.
+ """Represents a set of spherical harmonic coefficients.
Its harmonic partner spaces are :class:`~nifty5.domains.hp_space.HPSpace`
and :class:`~nifty5.domains.gl_space.GLSpace`.
diff --git a/nifty5/domains/log_rg_space.py b/nifty5/domains/log_rg_space.py
index 096e18ed213043db894638c612e36ccb940e1483..65cfb46a50aa3a9e16b26a32b6d6d5990ae798d1 100644
--- a/nifty5/domains/log_rg_space.py
+++ b/nifty5/domains/log_rg_space.py
@@ -16,30 +16,29 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from functools import reduce
+
import numpy as np
-from .. import dobj
from ..field import Field
-from ..sugar import exp
from .structured_domain import StructuredDomain
class LogRGSpace(StructuredDomain):
- """NIFTy subclass for logarithmic Cartesian grids.
+ '''Represents a logarithmic Cartesian grid.
Parameters
----------
shape : int or tuple of int
Number of grid points or numbers of gridpoints along each axis.
bindistances : float or tuple of float
- Distance between two grid points along each axis. These are
- measured on logarithmic scale and are constant therfore.
+ Logarithmic distance between two grid points along each axis.
+ Equidistant spacing of bins on logarithmic scale is assumed.
t_0 : float or tuple of float
- FIXME
+ Coordinate of pixel ndim*(1,).
harmonic : bool, optional
Whether the space represents a grid in position or harmonic space.
Default: False.
- """
+ '''
_needed_for_hash = ['_shape', '_bindistances', '_t_0', '_harmonic']
def __init__(self, shape, bindistances, t_0, harmonic=False):
@@ -77,6 +76,7 @@ class LogRGSpace(StructuredDomain):
@property
def t_0(self):
+ """np.ndarray : array of coordinates of pixel ndim*(1,)."""
return np.array(self._t_0)
def __repr__(self):
@@ -84,16 +84,49 @@ class LogRGSpace(StructuredDomain):
self.shape, self.harmonic))
def get_default_codomain(self):
+ """Returns a :class:`LogRGSpace` object representing the (position or
+ harmonic) partner domain of `self`, depending on `self.harmonic`. The
+ `bindistances` are transformed and `t_0` stays the same.
+
+ Returns
+ -------
+ LogRGSpace
+ The parter domain
+ """
codomain_bindistances = 1./(self.bindistances*self.shape)
return LogRGSpace(self.shape, codomain_bindistances, self._t_0, True)
def get_k_length_array(self):
+ """Generates array of distances to origin of the space.
+
+ Returns
+ -------
+ numpy.ndarray
+ Distances to origin of the space. If any index of the array is
+ zero then the distance is np.nan if self.harmonic True.
+ The dtype is float64, the shape is `self.shape`.
+
+ Raises
+ ------
+ NotImplementedError
+ If `self.harmonic` is False.
+ """
if not self.harmonic:
raise NotImplementedError
ks = self.get_k_array()
return Field.from_global_data(self, np.linalg.norm(ks, axis=0))
def get_k_array(self):
+ """Generates coordinates of the space.
+
+ Returns
+ -------
+ numpy.ndarray
+ Coordinates of the space. If one index of the array is zero the
+ corresponding coordinate is -np.inf (np.nan) if self.harmonic is
+ False (True).
+ The dtype is float64 and shape: `(len(self.shape),) + self.shape`.
+ """
ndim = len(self.shape)
k_array = np.zeros((ndim,) + self.shape)
dist = self.bindistances
diff --git a/nifty5/domains/power_space.py b/nifty5/domains/power_space.py
index 22828e73a2a184b0e014be6ebf42be211057d001..f6c357fd34dfcd7dfd6f28d68d89f25f81fd2645 100644
--- a/nifty5/domains/power_space.py
+++ b/nifty5/domains/power_space.py
@@ -22,7 +22,7 @@ from .structured_domain import StructuredDomain
class PowerSpace(StructuredDomain):
- """NIFTy class for spaces of power spectra.
+ """Represents non-equidistantly binned spaces for power spectra.
A power space is the result of a projection of a harmonic domain where
k-modes of equal length get mapped to one power index.
diff --git a/nifty5/domains/rg_space.py b/nifty5/domains/rg_space.py
index af0d44e7e6ecb5ffeb0b964b6a262f1631555515..2e1d2c0296c7a80bbdfe2d53bd9c849069323bee 100644
--- a/nifty5/domains/rg_space.py
+++ b/nifty5/domains/rg_space.py
@@ -24,7 +24,7 @@ from .structured_domain import StructuredDomain
class RGSpace(StructuredDomain):
- """NIFTy subclass for regular Cartesian grids.
+ """Represents a regular Cartesian grid.
Parameters
----------
diff --git a/nifty5/domains/structured_domain.py b/nifty5/domains/structured_domain.py
index c2b858ee939f4aae49f9b582caab6ee475d2caf2..be5deb75b23ca6f7bbac1c3160da250557aacd51 100644
--- a/nifty5/domains/structured_domain.py
+++ b/nifty5/domains/structured_domain.py
@@ -21,7 +21,7 @@ from .domain import Domain
class StructuredDomain(Domain):
- """The abstract base class for all structured NIFTy domains.
+ """The abstract base class for all structured domains.
An instance of a space contains information about the manifold's
geometry and enhances the functionality of Domain by methods that
@@ -50,7 +50,7 @@ class StructuredDomain(Domain):
@property
def total_volume(self):
- """float : Total domain volume
+ """float : Total domain volume.
Returns the sum over all the domain's pixel volumes.
"""
@@ -63,7 +63,7 @@ class StructuredDomain(Domain):
raise NotImplementedError
def get_k_length_array(self):
- """k vector lengths, if applicable,
+ """k vector lengths, if applicable.
Returns the length of the k vector for every pixel.
This method is only implemented for harmonic domains.
diff --git a/nifty5/field.py b/nifty5/field.py
index 8876fb03d0e0122adbaadf506e11976f2c2890df..722c2a9c0fd34787377009d0fff771f0d58b0fa3 100644
--- a/nifty5/field.py
+++ b/nifty5/field.py
@@ -23,18 +23,15 @@ from .domain_tuple import DomainTuple
class Field(object):
- _scalar_dom = DomainTuple.scalar_domain()
-
"""The discrete representation of a continuous field over multiple spaces.
- In NIFTy, Fields are used to store data arrays and carry all the needed
- metainformation (i.e. the domain) for operators to be able to work on them.
+ Stores data arrays and carries all the needed metainformation (i.e. the
+ domain) for operators to be able to operate on them.
Parameters
----------
domain : DomainTuple
the domain of the new Field
-
val : data_object
This object's global shape must match the domain shape
After construction, the object will no longer be writeable!
@@ -42,9 +39,11 @@ class Field(object):
Notes
-----
If possible, do not invoke the constructor directly, but use one of the
- many convenience functions for Field construction!
+ many convenience functions for instantiation!
"""
+ _scalar_dom = DomainTuple.scalar_domain()
+
def __init__(self, domain, val):
if not isinstance(domain, DomainTuple):
raise TypeError("domain must be of type DomainTuple")
@@ -337,7 +336,7 @@ class Field(object):
"""
if not isinstance(x, Field):
raise TypeError("The multiplier must be an instance of " +
- "the NIFTy field class")
+ "the Field class")
from .operators.outer_product_operator import OuterProduct
return OuterProduct(self, x.domain)(x)
@@ -360,7 +359,7 @@ class Field(object):
"""
if not isinstance(x, Field):
raise TypeError("The dot-partner must be an instance of " +
- "the NIFTy field class")
+ "the Field class")
if x._domain != self._domain:
raise ValueError("Domain mismatch")
diff --git a/nifty5/library/amplitude_operator.py b/nifty5/library/amplitude_operator.py
deleted file mode 100644
index c9e6172447709b4d1e5cb6d4816e719ea1c55d76..0000000000000000000000000000000000000000
--- a/nifty5/library/amplitude_operator.py
+++ /dev/null
@@ -1,161 +0,0 @@
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
-#
-# Copyright(C) 2013-2019 Max-Planck-Society
-#
-# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
-
-import numpy as np
-
-from ..domains.power_space import PowerSpace
-from ..field import Field
-from ..sugar import makeOp
-
-
-def _ceps_kernel(dof_space, k, a, k0):
- return a**2/(1 + (k/k0)**2)**2
-
-
-def _create_cepstrum_amplitude_field(domain, cepstrum):
- dim = len(domain.shape)
- shape = domain.shape
- q_array = domain.get_k_array()
-
- # Fill all non-zero modes
- no_zero_modes = (slice(1, None), )*dim
- ks = q_array[(slice(None), ) + no_zero_modes]
- cepstrum_field = np.zeros(shape)
- cepstrum_field[no_zero_modes] = cepstrum(ks)
-
- # Fill zero-mode subspaces
- for i in range(dim):
- fst_dims = (slice(None), )*i
- sl = fst_dims + (slice(1, None), )
- sl2 = fst_dims + (0, )
- cepstrum_field[sl2] = np.sum(cepstrum_field[sl], axis=i)
- return Field.from_global_data(domain, cepstrum_field)
-
-
-def CepstrumOperator(domain, a, k0):
- '''
- .. math::
- C(k) = \\left(\\frac{a}{1+(k/k0)^2}\\right)^2
- '''
- from ..operators.qht_operator import QHTOperator
- from ..operators.symmetrizing_operator import SymmetrizingOperator
-
- # FIXME a>0 k0>0
- qht = QHTOperator(target=domain)
- dof_space = qht.domain[0]
- sym = SymmetrizingOperator(domain)
- kern = lambda k: _ceps_kernel(dof_space, k, a, k0)
- cepstrum = _create_cepstrum_amplitude_field(dof_space, kern)
- return sym @ qht @ makeOp(cepstrum.sqrt())
-
-
-def SlopeOperator(domain, sm, sv, im, iv):
- '''
- Parameters
- ----------
-
- sm, sv : slope_mean = expected exponent of power law (e.g. -4),
- slope_variance (default=1)
-
- im, iv : y-intercept_mean, y-intercept_std of power_slope
-
- '''
- from ..operators.slope_operator import SlopeOperator
- from ..operators.offset_operator import OffsetOperator
-
- # sv, iv>0
-
- phi_mean = np.array([sm, im + sm*domain.t_0[0]])
- phi_sig = np.array([sv, iv])
- slope = SlopeOperator(domain)
- phi_mean = Field.from_global_data(slope.domain, phi_mean)
- phi_sig = Field.from_global_data(slope.domain, phi_sig)
- return slope(OffsetOperator(phi_mean)(makeOp(phi_sig)))
-
-
-def AmplitudeOperator(
- target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
- '''Operator for parametrizing smooth power spectra.
-
- The general guideline for setting up generative models in IFT is to
- transform the problem into the eigenbase of the prior and formulate the
- generative model in this base. This is done here for the case of a power
- spectrum which is smooth and has a linear component (both on
- double-logarithmic scale).
-
- This function assembles an :class:`Operator` which maps two a-priori white
- Gaussian random fields to a smooth power spectrum which is composed out of
- a linear and a smooth component.
-
- On double-logarithmic scale, i.e. both x and y-axis on logarithmic scale,
- the output of the generated operator is:
-
- AmplitudeOperator = 0.5*(smooth_component + linear_component)
-
- This is then exponentiated and exponentially binned (via a :class:`ExpTransform`).
-
- The prior on the linear component is parametrized by four real numbers,
- being expected value and prior variance on the slope and the y-intercept
- of the linear function.
-
- The prior on the smooth component is parametrized by two real numbers: the
- strength and the cutoff of the smoothness prior (see :class:`CepstrumOperator`).
-
- Parameters
- ----------
- n_pix : int
- Number of pixels of the space in which the .
- target : PowerSpace
- Target of the Operator.
- a : float
- Strength of smoothness prior (see :class:`CepstrumOperator`).
- k0 : float
- Cutoff of smothness prior in quefrency space (see :class:`CepstrumOperator`).
- sm : float
- Expected exponent of power law (see :class:`SlopeOperator`).
- sv : float
- Prior variance of exponent of power law (see :class:`SlopeOperator`).
- im : float
- Expected y-intercept of power law (see :class:`SlopeOperator`).
- iv : float
- Prior variance of y-intercept of power law (see :class:`SlopeOperator`).
-
- Returns
- -------
- Operator
- Operator which is defined on the space of white excitations fields and
- which returns on its target a power spectrum which consists out of a
- smooth and a linear part.
- '''
- from ..operators.exp_transform import ExpTransform
-
- if not (isinstance(n_pix, int) and isinstance(target, PowerSpace)):
- raise TypeError
-
- a, k0 = float(a), float(k0)
- sm, sv, im, iv = float(sm), float(sv), float(im), float(iv)
-
- et = ExpTransform(target, n_pix)
- dom = et.domain[0]
-
- dct = {'a': a, 'k0': k0}
- smooth = CepstrumOperator(dom, **dct).ducktape(keys[0])
-
- dct = {'sm': sm, 'sv': sv, 'im': im, 'iv': iv}
- linear = SlopeOperator(dom, **dct).ducktape(keys[1])
-
- return et((0.5*(smooth + linear)).exp())
diff --git a/nifty5/library/correlated_fields.py b/nifty5/library/correlated_fields.py
index a0e85ed0439a186f3ea82221d482ab030eb7928f..71567cda325c982fd514d7001763770a5691738a 100644
--- a/nifty5/library/correlated_fields.py
+++ b/nifty5/library/correlated_fields.py
@@ -15,64 +15,93 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+from functools import reduce
+
from ..domain_tuple import DomainTuple
from ..operators.contraction_operator import ContractionOperator
from ..operators.distributors import PowerDistributor
from ..operators.harmonic_operators import HarmonicTransformOperator
from ..operators.simple_linear_operators import ducktape
-from ..operators.scaling_operator import ScalingOperator
-def CorrelatedField(s_space, amplitude_operator, name='xi'):
- '''
- Function for construction of correlated fields
+def CorrelatedField(target, amplitude_operator, name='xi'):
+ '''Constructs an operator which turns a white Gaussian excitation field
+ into a correlated field.
+
+ This function returns an operator which implements:
+
+ ht @ (vol * A * xi),
+
+ where `ht` is a harmonic transform operator, `A` is the sqare root of the
+ prior covariance an `xi` is the excitation field.
Parameters
----------
- s_space : Domain
- Field domain
+ target : Domain, DomainTuple or tuple of Domain
+ Target of the operator. Must contain exactly one space.
amplitude_operator: Operator
- operator for correlation structure
name : string
- MultiField component name
+ :class:`MultiField` key for xi-field.
+
+ Returns
+ -------
+ Correlated field : Operator
'''
- h_space = s_space.get_default_codomain()
- ht = HarmonicTransformOperator(h_space, s_space)
+ tgt = DomainTuple.make(target)
+ if len(tgt) > 1:
+ raise ValueError
+ h_space = tgt[0].get_default_codomain()
+ ht = HarmonicTransformOperator(h_space, tgt[0])
p_space = amplitude_operator.target[0]
power_distributor = PowerDistributor(h_space, p_space)
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))
+ vol = h_space.scalar_dvol**-0.5
+ return ht(vol*A*ducktape(h_space, None, name))
-def MfCorrelatedField(s_space_spatial,
- s_space_energy,
- amplitude_operator_spatial,
- amplitude_operator_energy,
- name="xi"):
- '''
- Method for construction of correlated multi-frequency fields
+def MfCorrelatedField(target, amplitudes, name='xi'):
+ '''Constructs an operator which turns white Gaussian excitation fields
+ into a correlated field defined on a DomainTuple with two entries and two
+ separate correlation structures.
+
+ This operator may be used as a model for multi-frequency reconstructions
+ with a correlation structure in both spatial and energy direction.
+
+ Parameters
+ ----------
+ target : Domain, DomainTuple or tuple of Domain
+ Target of the operator. Must contain exactly one space.
+ amplitudes: iterable of Operator
+ List of two amplitude operators.
+ name : string
+ :class:`MultiField` key for xi-field.
+
+ Returns
+ -------
+ Correlated field : Operator
'''
- h_space_spatial = s_space_spatial.get_default_codomain()
- h_space_energy = s_space_energy.get_default_codomain()
- h_space = DomainTuple.make((h_space_spatial, h_space_energy))
- ht1 = HarmonicTransformOperator(h_space, target=s_space_spatial, space=0)
- ht2 = HarmonicTransformOperator(ht1.target, space=1)
- ht = ht2(ht1)
+ tgt = DomainTuple.make(target)
+ if len(tgt) != 2:
+ raise ValueError
+ if len(amplitudes) != 2:
+ raise ValueError
- p_space_spatial = amplitude_operator_spatial.target[0]
- p_space_energy = amplitude_operator_energy.target[0]
+ hsp = DomainTuple.make([tt.get_default_codomain() for tt in tgt])
+ ht1 = HarmonicTransformOperator(hsp, target=tgt[0], space=0)
+ ht2 = HarmonicTransformOperator(ht1.target, space=1)
+ ht = ht2 @ ht1
- pd_spatial = PowerDistributor(h_space, p_space_spatial, 0)
- pd_energy = PowerDistributor(pd_spatial.domain, p_space_energy, 1)
- pd = pd_spatial(pd_energy)
+ psp = [aa.target[0] for aa in amplitudes]
+ pd0 = PowerDistributor(hsp, psp[0], 0)
+ pd1 = PowerDistributor(pd0.domain, psp[1], 1)
+ pd = pd0 @ pd1
- dom_distr_spatial = ContractionOperator(pd.domain, 1).adjoint
- dom_distr_energy = ContractionOperator(pd.domain, 0).adjoint
+ dd0 = ContractionOperator(pd.domain, 1).adjoint
+ dd1 = ContractionOperator(pd.domain, 0).adjoint
+ d = [dd0, dd1]
- 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))
+ a = [dd @ amplitudes[ii] for ii, dd in enumerate(d)]
+ a = reduce(lambda x, y: x*y, a)
+ A = pd @ a
+ vol = reduce(lambda x, y: x*y, [sp.scalar_dvol**-0.5 for sp in hsp])
+ return ht(vol*A*ducktape(hsp, None, name))
diff --git a/nifty5/library/dynamic_operator.py b/nifty5/library/dynamic_operator.py
index ea385106e4eaf4d88b248f6099e4e508f67f669c..aad8fb2afcf5c8c37ac46cee0e349a6b0f94df74 100644
--- a/nifty5/library/dynamic_operator.py
+++ b/nifty5/library/dynamic_operator.py
@@ -147,7 +147,7 @@ def dynamic_operator(domain,
Parameters
----------
domain : RGSpace
- The space under consideration.
+ The position space in which the Green's function shall be constructed.
harmonic_padding : None, int, list of int
Amount of central padding in harmonic space in pixels. If None the
field is not padded at all.
@@ -158,9 +158,9 @@ def dynamic_operator(domain,
key : String
key for dynamics encoding parameter.
causal : boolean
- Whether or not the reconstructed dynamics should be causal in time.
+ Whether or not the Green's function shall be causal in time.
minimum_phase: boolean
- Whether or not the reconstructed dynamics should be minimum phase.
+ Whether or not the Green's function shall be a minimum phase filter.
Returns
-------
@@ -197,14 +197,14 @@ def dynamic_lightcone_operator(domain,
quant,
causal=True,
minimum_phase=False):
- '''Constructs an operator encoding the Green's function of a linear
- homogeneous dynamic system. The Green's function is constrained to be
- within a light cone.
+ '''Extends the functionality of :function: dynamic_operator to a Green's
+ function which is constrained to be within a light cone.
Parameters
----------
domain : RGSpace
- The space under consideration. Must have dim > 1.
+ The position space in which the Green's function shall be constructed.
+ It needs to have at least two dimensions.
harmonic_padding : None, int, list of int
Amount of central padding in harmonic space in pixels. If None the
field is not padded at all.
@@ -221,18 +221,17 @@ def dynamic_lightcone_operator(domain,
quant : float
Quantization of the light cone in pixels.
causal : boolean
- Whether or not the reconstructed dynamics should be causal in time.
+ Whether or not the Green's function shall be causal in time.
minimum_phase: boolean
- Whether or not the reconstructed dynamics should be minimum phase.
+ Whether or not the Green's function shall be a minimum phase filter.
Returns
-------
Operator
- The Operator encoding the dynamic Green's function in harmonic space
- when evaluated.
- dict
- A collection of sub-chains of Operator which can be used
- for plotting and evaluation.
+ The Operator encoding the dynamic Green's function in harmonic space.
+ Dictionary of Operator
+ A collection of sub-chains of Operator which can be used for plotting
+ and evaluation.
Notes
-----
diff --git a/nifty5/library/smooth_linear_amplitude.py b/nifty5/library/smooth_linear_amplitude.py
new file mode 100644
index 0000000000000000000000000000000000000000..331b270c87db7e99a80dc9422eedb74398c412ad
--- /dev/null
+++ b/nifty5/library/smooth_linear_amplitude.py
@@ -0,0 +1,198 @@
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+#
+# Copyright(C) 2013-2019 Max-Planck-Society
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+
+import numpy as np
+
+from ..domain_tuple import DomainTuple
+from ..domains.power_space import PowerSpace
+from ..field import Field
+from ..operators.exp_transform import ExpTransform
+from ..operators.offset_operator import OffsetOperator
+from ..operators.qht_operator import QHTOperator
+from ..operators.slope_operator import SlopeOperator
+from ..operators.symmetrizing_operator import SymmetrizingOperator
+from ..sugar import makeOp
+
+
+def _ceps_kernel(k, a, k0):
+ return (a/(1+np.sum((k.T/k0)**2, axis=-1).T))**2
+
+
+def CepstrumOperator(target, a, k0):
+ '''Turns a white Gaussian random field into a smooth field on a LogRGSpace.
+
+ Composed out of three operators:
+
+ sym @ qht @ diag(sqrt_ceps),
+
+ where sym is a :class:`SymmetrizingOperator`, qht is a :class:`QHTOperator`
+ and ceps is the so-called cepstrum:
+
+ .. math::
+ \\mathrm{sqrt\_ceps}(k) = \\frac{a}{1+(k/k0)^2}
+
+ These operators are combined in this fashion in order to generate:
+
+ - A field which is smooth, i.e. second derivatives are punished (note
+ that the sqrt-cepstrum is essentially proportional to 1/k**2).
+
+ - A field which is symmetric around the pixel in the middle of the space.
+ This is result of the :class:`SymmetrizingOperator` and needed in order
+ to decouple the degrees of freedom at the beginning and the end of the
+ amplitude whenever :class:`CepstrumOperator` is used as in
+ :class:`SLAmplitude`.
+
+ The prior on the zero mode (or zero subspaces for more than one dimensions)
+ is the integral of the prior over all other modes along the corresponding
+ axis.
+
+ Parameters
+ ----------
+ target : LogRGSpace
+ Target domain of the operator, needs to be non-harmonic.
+ a : float
+ Cutoff of smoothness prior (positive only). Controls the
+ regularization of the inverse laplace operator to be finite at zero.
+ Larger values for the cutoff results in a weaker constraining prior.
+ k0 : float, list of float
+ Strength of smothness prior in quefrency space (positive only) along
+ each axis. If float then the strength is the same along each axis.
+ Larger values result in a weaker constraining prior.
+ '''
+ a = float(a)
+ target = DomainTuple.make(target)
+ if a <= 0:
+ raise ValueError
+ if len(target) > 1 or target[0].harmonic:
+ raise TypeError
+ if isinstance(k0, (float, int)):
+ k0 = np.array([k0]*len(target.shape))
+ else:
+ k0 = np.array(k0)
+ if len(k0) != len(target.shape):
+ raise ValueError
+ if np.any(np.array(k0) <= 0):
+ raise ValueError
+
+ qht = QHTOperator(target)
+ dom = qht.domain[0]
+ sym = SymmetrizingOperator(target)
+
+ # Compute cepstrum field
+ dim = len(dom.shape)
+ shape = dom.shape
+ q_array = dom.get_k_array()
+ # Fill all non-zero modes
+ no_zero_modes = (slice(1, None),)*dim
+ ks = q_array[(slice(None),) + no_zero_modes]
+ cepstrum_field = np.zeros(shape)
+ cepstrum_field[no_zero_modes] = _ceps_kernel(ks, a, k0)
+ # Fill zero-mode subspaces
+ for i in range(dim):
+ fst_dims = (slice(None),)*i
+ sl = fst_dims + (slice(1, None),)
+ sl2 = fst_dims + (0,)
+ cepstrum_field[sl2] = np.sum(cepstrum_field[sl], axis=i)
+ cepstrum = Field.from_global_data(dom, cepstrum_field)
+
+ return sym @ qht @ makeOp(cepstrum.sqrt())
+
+
+def SLAmplitude(target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
+ '''Operator for parametrizing smooth amplitudes (square roots of power
+ spectra).
+
+ The general guideline for setting up generative models in IFT is to
+ transform the problem into the eigenbase of the prior and formulate the
+ generative model in this base. This is done here for the case of an
+ amplitude which is smooth and has a linear component (both on
+ double-logarithmic scale).
+
+ This function assembles an :class:`Operator` which maps two a-priori white
+ Gaussian random fields to a smooth amplitude which is composed out of
+ a linear and a smooth component.
+
+ On double-logarithmic scale, i.e. both x and y-axis on logarithmic scale,
+ the output of the generated operator is:
+
+ AmplitudeOperator = 0.5*(smooth_component + linear_component)
+
+ This is then exponentiated and exponentially binned (in this order).
+
+ The prior on the linear component is parametrized by four real numbers,
+ being expected value and prior variance on the slope and the y-intercept
+ of the linear function.
+
+ The prior on the smooth component is parametrized by two real numbers: the
+ strength and the cutoff of the smoothness prior (see :class:`CepstrumOperator`).
+
+ Parameters
+ ----------
+ n_pix : int
+ Number of pixels of the space in which the .
+ target : PowerSpace
+ Target of the Operator.
+ a : float
+ Strength of smoothness prior (see :class:`CepstrumOperator`).
+ k0 : float
+ Cutoff of smothness prior in quefrency space (see
+ :class:`CepstrumOperator`).
+ sm : float
+ Expected exponent of power law.
+ sv : float
+ Prior standard deviation of exponent of power law.
+ im : float
+ Expected y-intercept of power law. This is the value at t_0 of the
+ LogRGSpace (see :class:`ExpTransform`).
+ iv : float
+ Prior standard deviation of y-intercept of power law.
+
+ Returns
+ -------
+ Operator
+ Operator which is defined on the space of white excitations fields and
+ which returns on its target a power spectrum which consists out of a
+ smooth and a linear part.
+ '''
+ if not (isinstance(n_pix, int) and isinstance(target, PowerSpace)):
+ raise TypeError
+
+ a, k0 = float(a), float(k0)
+ sm, sv, im, iv = float(sm), float(sv), float(im), float(iv)
+ if sv <= 0 or iv <= 0:
+ raise ValueError
+
+ et = ExpTransform(target, n_pix)
+ dom = et.domain[0]
+
+ # Smooth component
+ dct = {'a': a, 'k0': k0}
+ smooth = CepstrumOperator(dom, **dct).ducktape(keys[0])
+
+ # Linear component
+ sl = SlopeOperator(dom)
+ mean = np.array([sm, im + sm*dom.t_0[0]])
+ sig = np.array([sv, iv])
+ mean = Field.from_global_data(sl.domain, mean)
+ sig = Field.from_global_data(sl.domain, sig)
+ linear = (sl @ OffsetOperator(mean) @ makeOp(sig)).ducktape(keys[1])
+
+ # Combine linear and smooth component
+ loglog_ampl = 0.5*(smooth + linear)
+
+ # Go from loglog-space to linear-linear-space
+ return et @ loglog_ampl.exp()
diff --git a/nifty5/operators/diagonal_operator.py b/nifty5/operators/diagonal_operator.py
index ccb2d4196c3d1a0b725d9b069b53f50ea056ff9f..3b2812182e3ebdf4dbf22cabae289527e3579c7b 100644
--- a/nifty5/operators/diagonal_operator.py
+++ b/nifty5/operators/diagonal_operator.py
@@ -24,7 +24,7 @@ from .endomorphic_operator import EndomorphicOperator
class DiagonalOperator(EndomorphicOperator):
- """NIFTy class for diagonal operators.
+ """Represents a linear operator which is diagonal.
The NIFTy DiagonalOperator class is a subclass derived from the
EndomorphicOperator. It multiplies an input field pixel-wise with its
diff --git a/nifty5/operators/endomorphic_operator.py b/nifty5/operators/endomorphic_operator.py
index e97e98e7e7181293a09f75f90a6a16a65c8f67d1..9b1d797f4d2dc0bdaaac15486506f74b2f3b03e0 100644
--- a/nifty5/operators/endomorphic_operator.py
+++ b/nifty5/operators/endomorphic_operator.py
@@ -21,12 +21,9 @@ from .linear_operator import LinearOperator
class EndomorphicOperator(LinearOperator):
- """NIFTy class for endomorphic operators.
-
- The NIFTy EndomorphicOperator describes linear operators with identical
- domain and target.
+ """Represents a linear operator which is endomorphic, i.e. one which
+ has identical domain and target.
"""
-
@property
def target(self):
"""DomainTuple : returns :attr:`domain`
diff --git a/nifty5/operators/exp_transform.py b/nifty5/operators/exp_transform.py
index 16c26ef67f4fb72fd007eb70fc7c0ad8d56d439d..78673d8160f3d422dd95163e25f7493de8eff2c0 100644
--- a/nifty5/operators/exp_transform.py
+++ b/nifty5/operators/exp_transform.py
@@ -34,6 +34,8 @@ class ExpTransform(LinearOperator):
Then transforms between this log-space and its target, which is defined in
normal units.
+ FIXME Write something on t_0 of domain space
+
E.g: A field in log-log-space can be transformed into log-norm-space,
that is the y-axis stays logarithmic, but the x-axis is transfromed.
diff --git a/nifty5/operators/offset_operator.py b/nifty5/operators/offset_operator.py
index ca4ffea26aaf7ccaea04bd49218f82ac505e3e21..a8ebba807008a2e7a4f340614b101d54b0c4ce0f 100644
--- a/nifty5/operators/offset_operator.py
+++ b/nifty5/operators/offset_operator.py
@@ -19,6 +19,12 @@ from .operator import Operator
class OffsetOperator(Operator):
+ '''Shifts the input by a fixed field.
+
+ Parameters
+ ----------
+ field : Field
+ The field by which the input is shifted.'''
def __init__(self, field):
self._field = field
self._domain = self._target = field.domain
diff --git a/nifty5/operators/operator.py b/nifty5/operators/operator.py
index 9893c35fbf04305b789e499f08ab4a0de34fccb4..3b2c8a102bb8992576bf7134d9056dd35b23fa7f 100644
--- a/nifty5/operators/operator.py
+++ b/nifty5/operators/operator.py
@@ -26,16 +26,23 @@ class Operator(NiftyMetaBase()):
@property
def domain(self):
- """DomainTuple or MultiDomain : the operator's input domain
+ """The domain on which the Operator's input Field is defined.
- The domain on which the Operator's input Field is defined."""
+ Returns
+ -------
+ domain : DomainTuple or MultiDomain
+ """
return self._domain
@property
def target(self):
- """DomainTuple or MultiDomain : the operator's output domain
+ """The domain on which the Operator's output Field is defined.
+
+ Returns
+ -------
+ target : DomainTuple or MultiDomain
+ """
- The domain on which the Operator's output Field is defined."""
return self._target
@staticmethod
@@ -103,10 +110,19 @@ class Operator(NiftyMetaBase()):
return _OpChain.make((_Clipper(self.target, min, max), self))
def apply(self, x):
+ '''Applies the operator to a Field or MultiField.
+
+ Parameters
+ ----------
+ x : Field or MultiField
+ Input on which the operator shall act. Needs to be defined on
+ :attr:`domain`.
+ '''
raise NotImplementedError
def force(self, x):
- """Extract correct subset of domain of x and apply operator."""
+ """Extract subset of domain of x according to `self.domain` and apply
+ operator."""
return self.apply(x.extract(self.domain))
def _check_input(self, x):
diff --git a/nifty5/operators/qht_operator.py b/nifty5/operators/qht_operator.py
index e327d6758a92e3b81eefbc172cf25a22e151fc5b..ff942521614f07a5944b0dd2547bc4e6c8d17758 100644
--- a/nifty5/operators/qht_operator.py
+++ b/nifty5/operators/qht_operator.py
@@ -26,9 +26,9 @@ from .linear_operator import LinearOperator
class QHTOperator(LinearOperator):
"""Does a Hartley transform on LogRGSpace
- This operator takes a field on a LogRGSpace and transforms it
- according to the Hartley transform. The zero modes are not transformed
- because they are infinitely far away.
+ This operator takes a field on a LogRGSpace and performs a Hartley
+ transform. The zero modes are not transformed because they are infinitely
+ far away in LogRGSpaces.
Parameters
----------
diff --git a/nifty5/operators/regridding_operator.py b/nifty5/operators/regridding_operator.py
index b6efcfac840ab720e71812c9d61fa51bcf645764..aad93b5cb18310c657c4ea34fe1e0cce83a3447c 100644
--- a/nifty5/operators/regridding_operator.py
+++ b/nifty5/operators/regridding_operator.py
@@ -26,6 +26,19 @@ from .linear_operator import LinearOperator
class RegriddingOperator(LinearOperator):
+ '''Linearly interpolates a RGSpace to an RGSpace with coarser resolution.
+
+ Parameters
+ ----------
+ domain : Domain, DomainTuple or tuple of Domain
+ domain[space] needs to be an :class:`RGSpace`.
+ new_shape : tuple of int
+ Shape of the space which domain[space] is replaced by. Each entry must
+ be smaller or equal to the respective entry in `domain[space].shape`.
+ space : int
+ Index of space in `domain` on which the operator shall act.
+ Default is 0.
+ '''
def __init__(self, domain, new_shape, space=0):
self._domain = DomainTuple.make(domain)
self._space = infer_space(self._domain, space)
@@ -38,6 +51,8 @@ class RegriddingOperator(LinearOperator):
raise ValueError("Shape mismatch")
if any([a > b for a, b in zip(new_shape, dom.shape)]):
raise ValueError("New shape must not be larger than old shape")
+ if any([ii <= 0 for ii in new_shape]):
+ raise ValueError('New shape must not be zero or negative.')
newdist = tuple(dom.distances[i]*dom.shape[i]/new_shape[i]
for i in range(len(dom.shape)))
diff --git a/nifty5/operators/scaling_operator.py b/nifty5/operators/scaling_operator.py
index d0748d08806806495b8119460831197e81d1ce32..708301a1011a488a070662383690977f0ebd51fc 100644
--- a/nifty5/operators/scaling_operator.py
+++ b/nifty5/operators/scaling_operator.py
@@ -24,9 +24,6 @@ from .endomorphic_operator import EndomorphicOperator
class ScalingOperator(EndomorphicOperator):
"""Operator which multiplies a Field with a scalar.
- The NIFTy ScalingOperator class is a subclass derived from the
- EndomorphicOperator. It multiplies an input field with a given factor.
-
Parameters
----------
factor : scalar
diff --git a/nifty5/operators/slope_operator.py b/nifty5/operators/slope_operator.py
index ac4aa8e18e39ed59769642123f2b2d4e62aa4fd0..1c84f6d5469684ccfccf0fb664c79ffc0787f8d0 100644
--- a/nifty5/operators/slope_operator.py
+++ b/nifty5/operators/slope_operator.py
@@ -25,30 +25,30 @@ from .linear_operator import LinearOperator
class SlopeOperator(LinearOperator):
- """Creates a slope on target.
+ """Evaluates a line on a LogRGSpace given slope and y-intercept
- This operator creates a field on a LogRGSpace, which is created according
- to a slope of given entries, (mean, y-intercept). The slope mean is the
- power law of the field in normal-space.
+ Slope and y-intercept of this line are the two parameters which are
+ defined on an UnstructeredDomain (in this order) which is the domain of
+ the operator. Being a LogRGSpace each pixel has a well-defined coordinate
+ value.
+
+ y-intercept is defined to be the value at t_0 of the target.
Parameters
----------
- domain : domain or DomainTuple, shape=(2,)
- It has to be an UnstructuredDomain.
- The domain of the slope mean and the y-intercept mean.
- target : domain or DomainTuple
- The output domain has to a LogRGSpace
- sigmas : np.array, shape=(2,)
- The slope variance and the y-intercept variance.
+ target : LogRGSpace
+ The target of the operator which needs to be one-dimensional.
"""
def __init__(self, target):
- if not isinstance(target, LogRGSpace):
+ self._target = DomainTuple.make(target)
+ if len(self._target) > 1:
+ raise TypeError
+ if len(self._target[0].shape) > 1:
+ raise TypeError
+ if not isinstance(self._target[0], LogRGSpace):
raise TypeError
- if len(target.shape) != 1:
- raise ValueError("Slope Operator only works for ndim == 1")
self._domain = DomainTuple.make(UnstructuredDomain((2,)))
- self._target = DomainTuple.make(target)
self._capability = self.TIMES | self.ADJOINT_TIMES
pos = self.target[0].get_k_array() - self.target[0].t_0[0]
self._pos = pos[0, 1:]
diff --git a/nifty5/operators/symmetrizing_operator.py b/nifty5/operators/symmetrizing_operator.py
index 6b149c20c165e1f0338457fdf6e23f8eed1c99dc..4f38603a4a2ad773737fafbdc8b0a989f81b128d 100644
--- a/nifty5/operators/symmetrizing_operator.py
+++ b/nifty5/operators/symmetrizing_operator.py
@@ -23,6 +23,17 @@ from .endomorphic_operator import EndomorphicOperator
class SymmetrizingOperator(EndomorphicOperator):
+ '''Subtracts the field axes-wise in reverse order from itself. The slice of
+ all elements with at least one index being zero is not touched.
+
+ Parameters
+ ----------
+ domain : Domain, DomainTuple or tuple of Domain
+ Domain of the operator. `domain[space]` needs to be a non-harmonic
+ :class:`LogRGSpace`.
+ space : int
+ Index of space in domain on which the operator shall act. Default is 0.
+ '''
def __init__(self, domain, space=0):
self._domain = DomainTuple.make(domain)
self._capability = self.TIMES | self.ADJOINT_TIMES
diff --git a/test/test_model_gradients.py b/test/test_model_gradients.py
index 90fbdd4e5fe2dda2d341f69caaaf8a704cd92eb6..e30e0ce5497683551c5c301c20301ca9ec41c667 100644
--- a/test/test_model_gradients.py
+++ b/test/test_model_gradients.py
@@ -92,8 +92,7 @@ def testModelLibrary(space, seed):
Npixdof, ceps_a, ceps_k, sm, sv, im, iv = 4, 0.5, 2., 3., 1.5, 1.75, 1.3
np.random.seed(seed)
domain = ift.PowerSpace(space.get_default_codomain())
- model = ift.AmplitudeOperator(domain, Npixdof, ceps_a, ceps_k, sm, sv, im,
- iv)
+ model = ift.SLAmplitude(domain, Npixdof, ceps_a, ceps_k, sm, sv, im, iv)
assert_(isinstance(model, ift.Operator))
S = ift.ScalingOperator(1., model.domain)
pos = S.draw_sample()
@@ -104,6 +103,12 @@ def testModelLibrary(space, seed):
pos = S.draw_sample()
ift.extra.check_value_gradient_consistency(model2, pos, ntries=20)
+ domtup = ift.DomainTuple.make((space, space))
+ model3 = ift.MfCorrelatedField(domtup, [model, model])
+ S = ift.ScalingOperator(1., model3.domain)
+ pos = S.draw_sample()
+ ift.extra.check_value_gradient_consistency(model3, pos, ntries=20)
+
def testPointModel(space, seed):
S = ift.ScalingOperator(1., space)