Merge remote-tracking branch 'origin/NIFTy_5' into spectral_plotting

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)

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

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,

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

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`.

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`.

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`.

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

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.

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
     ----------

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.

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")

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))

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
     -----

nifty5/library/amplitude_operator.py → nifty5/library/smooth_linear_amplitude.py

@@ -17,88 +17,113 @@
 
 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(dof_space, k, a, k0):
-    return a**2/(1 + (k/k0)**2)**2
+def _ceps_kernel(k, a, k0):
+    return (a/(1+np.sum((k.T/k0)**2, axis=-1).T))**2
 
 
-def _create_cepstrum_amplitude_field(domain, cepstrum):
-    dim = len(domain.shape)
-    shape = domain.shape
-    q_array = domain.get_k_array()
+def CepstrumOperator(target, a, k0):
+    '''Turns a white Gaussian random field into a smooth field on a LogRGSpace.
 
-    # 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)
+    Composed out of three operators:
 
-    # 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)
+        sym @ qht @ diag(sqrt_ceps),
 
+    where sym is a :class:`SymmetrizingOperator`, qht is a :class:`QHTOperator`
+    and ceps is the so-called cepstrum:
 
-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())
+        \\mathrm{sqrt\_ceps}(k) = \\frac{a}{1+(k/k0)^2}
 
+    These operators are combined in this fashion in order to generate:
 
-def SlopeOperator(domain, sm, sv, im, iv):
-    '''
-    Parameters
-    ----------
+    - A field which is smooth, i.e. second derivatives are punished (note
+      that the sqrt-cepstrum is essentially proportional to 1/k**2).
 
-    sm, sv : slope_mean = expected exponent of power law (e.g. -4),
-                slope_variance (default=1)
+    - 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`.
 
-    im, iv : y-intercept_mean, y-intercept_std  of power_slope
+    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.
     '''
-    from ..operators.slope_operator import SlopeOperator
-    from ..operators.offset_operator import OffsetOperator
-
-    # sv, iv>0
+    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)
 
-    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)))
+    return sym @ qht @ makeOp(cepstrum.sqrt())
 
 
-def AmplitudeOperator(
-        target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
-    '''Operator for parametrizing smooth power spectra.
+def SLAmplitude(target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):