AmplitudeOperator: Add docu, do not create PowerSpace, remove zero_mode

Furthermore, the zero_mode keyword is gone for AmplitudeOperator. It is not used
in the demos and is not really appropriate from a theoretical point of view
anyways. We will add in a better option later.

demos/getting_started_3.py

@@ -46,32 +46,32 @@ if __name__ == '__main__':
     mode = 1
 
     position_space = ift.RGSpace([128, 128])
+    harmonic_space = position_space.get_default_codomain()
+    ht = ift.HarmonicTransformOperator(harmonic_space, position_space)
+    power_space = ift.PowerSpace(harmonic_space)
 
     # Set up an amplitude operator for the field
-    # The parameters mean:
-    # 64 spectral bins
-    #
-    # Spectral smoothness (affects Gaussian process part)
-    # 3 = relatively high variance of spectral curbvature
-    # 0.4 = quefrency mode below which cepstrum flattens
-    #
-    # Power-law part of spectrum:
-    # -5 = preferred power-law slope
-    # 0.5 = low variance of power-law slope
-    # 0.4 = y-intercept mean
-    # 0.3 = relatively high y-intercept variance
-    A = ift.AmplitudeOperator(position_space, 64, 3, 0.4, -5., 0.5, 0.4, 0.3)
+    dct = {
+        'target': power_space,
+        'n_pix': 64,  # 64 spectral bins
+
+        # Spectral smoothness (affects Gaussian process part)
+        'a': 3,  # relatively high variance of spectral curbvature
+        'k0': .4,  # quefrency mode below which cepstrum flattens
+
+        # Power-law part of spectrum:
+        'sm': -5,  # preferred power-law slope
+        'sv': .5,  # low variance of power-law slope
+        'im': .4,  # y-intercept mean
+        'iv': .3  # relatively high y-intercept variance
+    }
+    A = ift.AmplitudeOperator(**dct)
 
     # 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]
     power_distributor = ift.PowerDistributor(harmonic_space, power_space)
-
-    vol = ift.ScalingOperator(harmonic_space.scalar_dvol**(-0.5),
-                              harmonic_space)
-    correlated_field = ht(
-        vol(power_distributor(A))*ift.ducktape(harmonic_space, None, 'xi'))
+    vol = harmonic_space.scalar_dvol**-0.5
+    xi = ift.ducktape(harmonic_space, None, 'xi')
+    correlated_field = ht(vol*power_distributor(A)*xi)
     # Alternatively, one can use:
     # correlated_field = ift.CorrelatedField(position_space, A)
 

nifty5/library/amplitude_operator.py

@@ -19,76 +19,51 @@ import numpy as np
 
 from ..domains.power_space import PowerSpace
 from ..field import Field
-from ..sugar import makeOp, sqrt
+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):
-    """Creates a ...
-    Writes the sum of all modes into the zero-mode.
-
-    Parameters
-    ----------
-    domain: ???
-        ???
-    cepstrum: Callable
-        ???
-    """
-
+def _create_cepstrum_amplitude_field(domain, cepstrum):
     dim = len(domain.shape)
     shape = domain.shape
-
     q_array = domain.get_k_array()
 
-    # Fill cepstrum field (all non-zero modes)
-    no_zero_modes = (slice(1, None),)*dim
-    ks = q_array[(slice(None),) + no_zero_modes]
+    # 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 cepstrum field (zero-mode subspaces)
+    # Fill zero-mode subspaces
     for i in range(dim):
-        # Prepare indices
-        fst_dims = (slice(None),)*i
-        sl = fst_dims + (slice(1, None),)
-        sl2 = fst_dims + (0,)
-
-        # Do summation
+        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(logk_space, ceps_a, ceps_k, zero_mode=True):
+def CepstrumOperator(domain, a, k0):
     '''
-    Parameters
-    ----------
-    ceps_a, ceps_k0 : Smoothness parameters in ceps_kernel
-                        eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2
-                        a = ceps_a,  k0 = ceps_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
-    qht = QHTOperator(target=logk_space)
+
+    # FIXME a>0 k0>0
+    qht = QHTOperator(target=domain)
     dof_space = qht.domain[0]
-    sym = SymmetrizingOperator(logk_space)
-    kern = lambda k: _ceps_kernel(dof_space, k, ceps_a, ceps_k)
-    cepstrum = create_cepstrum_amplitude_field(dof_space, kern)
-    res = sym(qht(makeOp(sqrt(cepstrum))))
-    if not zero_mode:
-        shp = res.target.shape
-        mask = np.ones(shp)
-        mask[(0,)*len(shp)] = 0.
-        mask = makeOp(Field.from_global_data(res.target, mask))
-        res = mask(res)
-    return res
-
-
-def _SlopePowerSpectrum(logk_space, sm, sv, im, iv):
+    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
     ----------
@@ -97,54 +72,90 @@ def _SlopePowerSpectrum(logk_space, sm, sv, im, iv):
                 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
-    phi_mean = np.array([sm, im + sm*logk_space.t_0[0]])
+
+    # sv, iv>0
+
+    phi_mean = np.array([sm, im + sm*domain.t_0[0]])
     phi_sig = np.array([sv, iv])
-    slope = SlopeOperator(logk_space)
+    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(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
-                      keys=['tau', 'phi'], zero_mode=True):
-    ''' Operator for parametrizing smooth power spectra.
+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)
 
-    Computes a smooth power spectrum.
-    Output is defined on a PowerSpace.
+    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
     ----------
-    Npixdof : int
-        #pix in dof_space
-    ceps_a : float
-        Smoothness parameters in ceps_kernel eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2 a = ceps_a,  k0 = ceps_k0
-    ceps_k0 : float
-        Smoothness parameters in ceps_kernel eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2 a = ceps_a,  k0 = ceps_k0
+    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
-        slope_mean = expected exponent of power law (e.g. -4)
+        Expected exponent of power law (see :class:`SlopeOperator`).
     sv : float
-        slope_variance (default=1)
+        Prior variance of exponent of power law (see :class:`SlopeOperator`).
     im : float
-        y-intercept_mean
+        Expected y-intercept of power law (see :class:`SlopeOperator`).
     iv : float
-        y-intercept_variance  of power_slope
+        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
-    from ..operators.scaling_operator import ScalingOperator
 
-    h_space = s_space.get_default_codomain()
-    et = ExpTransform(PowerSpace(h_space), Npixdof)
-    logk_space = et.domain[0]
+    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])
 
-    smooth = _CepstrumOperator(logk_space, ceps_a, ceps_k, zero_mode)
-    smooth = smooth.ducktape(keys[0])
-    linear = _SlopePowerSpectrum(logk_space, sm, sv, im, iv)
-    linear = linear.ducktape(keys[1])
+    dct = {'sm': sm, 'sv': sv, 'im': im, 'iv': iv}
+    linear = SlopeOperator(dom, **dct).ducktape(keys[1])
 
-    fac = ScalingOperator(0.5, smooth.target)
-    return et((fac(smooth + linear)).exp())
+    return et((0.5*(smooth + linear)).exp())

test/test_model_gradients.py

@@ -17,6 +17,7 @@
 
 import numpy as np
 import pytest
+from numpy.testing import assert_
 
 import nifty5 as ift
 
@@ -53,11 +54,11 @@ def testBasics(space, seed):
 @pmp('type2', ['Variable'])
 def testBinary(type1, type2, space, seed):
     dom1 = ift.MultiDomain.make({'s1': space})
-    # FIXME Remove?
-    lin1 = _make_linearization(type1, dom1, seed)
     dom2 = ift.MultiDomain.make({'s2': space})
-    # FIXME Remove?
-    lin2 = _make_linearization(type2, dom2, seed)
+
+    # FIXME Remove this?
+    _make_linearization(type1, dom1, seed)
+    _make_linearization(type2, dom2, seed)
 
     dom = ift.MultiDomain.union((dom1, dom2))
     select_s1 = ift.ducktape(None, dom, "s1")
@@ -74,8 +75,7 @@ def testBinary(type1, type2, space, seed):
     model = 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.sigmoid(
-        ift.ScalingOperator(2.456, space)(select_s1*select_s2))
+    model = ift.sigmoid(2.456*(select_s1*select_s2))
     pos = ift.from_random("normal", dom)
     ift.extra.check_value_gradient_consistency(model, pos, ntries=20)
     pos = ift.from_random("normal", dom)
@@ -91,8 +91,10 @@ def testModelLibrary(space, seed):
     # Tests amplitude model and coorelated field model
     Npixdof, ceps_a, ceps_k, sm, sv, im, iv = 4, 0.5, 2., 3., 1.5, 1.75, 1.3
     np.random.seed(seed)
-    model = ift.AmplitudeOperator(space, Npixdof, ceps_a, ceps_k, sm, sv, im,
+    domain = ift.PowerSpace(space.get_default_codomain())
+    model = ift.AmplitudeOperator(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()
     ift.extra.check_value_gradient_consistency(model, pos, ntries=20)
@@ -112,28 +114,38 @@ def testPointModel(space, seed):
     # FIXME All those cdfs and ppfs are not very accurate
     ift.extra.check_value_gradient_consistency(model, pos, tol=1e-2, ntries=20)
 
-@pmp('domain', [ift.RGSpace(64, distances=.789),
-                ift.RGSpace([32, 32], distances=.789),
-                ift.RGSpace([32, 32, 8], distances=.789)])
+
+@pmp('domain', [
+    ift.RGSpace(64, distances=.789),
+    ift.RGSpace([32, 32], distances=.789),
+    ift.RGSpace([32, 32, 8], distances=.789)
+])
 @pmp('causal', [True, False])
 @pmp('minimum_phase', [True, False])
 @pmp('seed', [4, 78, 23])
 def testDynamicModel(domain, causal, minimum_phase, seed):
-    model, _ = ift.dynamic_operator(domain,None,1.,1.,'f',
-                                    causal = causal,
-                                    minimum_phase = minimum_phase)
+    model, _ = ift.dynamic_operator(
+        domain, None, 1., 1., 'f', causal=causal, minimum_phase=minimum_phase)
     S = ift.ScalingOperator(1., model.domain)
     pos = S.draw_sample()
     # FIXME I dont know why smaller tol fails for 3D example
-    ift.extra.check_value_gradient_consistency(model, pos, tol=1e-5,
-                                               ntries=20)
+    ift.extra.check_value_gradient_consistency(model, pos, tol=1e-5, ntries=20)
     if len(domain.shape) > 1:
-        model, _ = ift.dynamic_lightcone_operator(domain,None,3.,1.,
-                                                  'f','c',1.,5,
-                                                  causal = causal,
-                                                  minimum_phase = minimum_phase)
+        dct = {
+            'domain': domain,
+            'harmonic_padding': None,
+            'sm_s0': 3.,
+            'sm_x0': 1.,
+            'key': 'f',
+            'lightcone_key': 'c',
+            'sigc': 1.,
+            'quant': 5,
+            'causal': causal,
+            'minimum_phase': minimum_phase
+        }
+        model, _ = ift.dynamic_lightcone_operator(**dct)
         S = ift.ScalingOperator(1., model.domain)
         pos = S.draw_sample()
         # FIXME I dont know why smaller tol fails for 3D example
-        ift.extra.check_value_gradient_consistency(model, pos, tol=1e-5,
-                                                   ntries=20)
+        ift.extra.check_value_gradient_consistency(
+            model, pos, tol=1e-5, ntries=20)