smooth_linear_amplitude.py 5.26 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1 2 3 4 5 6 7 8 9 10 11 12 13
# 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/>.
#
14
# Copyright(C) 2013-2019 Max-Planck-Society
Martin Reinecke's avatar
Martin Reinecke committed
15
#
16
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
Philipp Arras's avatar
Philipp Arras committed
17

18
import numpy as np
Philipp Arras's avatar
Philipp Arras committed
19

Martin Reinecke's avatar
Martin Reinecke committed
20
from ..domains.power_space import PowerSpace
21
from ..field import Field
22
from ..sugar import makeOp
23 24 25


def _ceps_kernel(dof_space, k, a, k0):
26
    return a**2/(1 + (k/k0)**2)**2
27 28


29
def _create_cepstrum_amplitude_field(domain, cepstrum):
30 31
    dim = len(domain.shape)
    shape = domain.shape
32
    q_array = domain.get_k_array()
33

34
    # Fill all non-zero modes
Philipp Arras's avatar
Philipp Arras committed
35 36
    no_zero_modes = (slice(1, None),)*dim
    ks = q_array[(slice(None),) + no_zero_modes]
37 38 39
    cepstrum_field = np.zeros(shape)
    cepstrum_field[no_zero_modes] = cepstrum(ks)

40
    # Fill zero-mode subspaces
41
    for i in range(dim):
Philipp Arras's avatar
Philipp Arras committed
42 43 44
        fst_dims = (slice(None),)*i
        sl = fst_dims + (slice(1, None),)
        sl2 = fst_dims + (0,)
45
        cepstrum_field[sl2] = np.sum(cepstrum_field[sl], axis=i)
46
    return Field.from_global_data(domain, cepstrum_field)
Martin Reinecke's avatar
Martin Reinecke committed
47

Martin Reinecke's avatar
Martin Reinecke committed
48

49
def CepstrumOperator(domain, a, k0):
50
    '''
51 52
    .. math::
        C(k) = \\left(\\frac{a}{1+(k/k0)^2}\\right)^2
53
    '''
Philipp Arras's avatar
Philipp Arras committed
54 55
    from ..operators.qht_operator import QHTOperator
    from ..operators.symmetrizing_operator import SymmetrizingOperator
56

Philipp Arras's avatar
Philipp Arras committed
57 58 59
    if a <= 0 or k0 <= 0:
        raise ValueError

60
    qht = QHTOperator(target=domain)
Philipp Arras's avatar
Philipp Arras committed
61
    dof_space = qht.domain[0]
62 63 64 65 66 67
    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())


68 69 70
def SLAmplitude(target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
    '''Operator for parametrizing smooth amplitudes (square roots of power
    spectra).
71 72 73

    The general guideline for setting up generative models in IFT is to
    transform the problem into the eigenbase of the prior and formulate the
74 75
    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
76 77 78
    double-logarithmic scale).

    This function assembles an :class:`Operator` which maps two a-priori white
79
    Gaussian random fields to a smooth amplitude which is composed out of
80 81 82 83 84 85
    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)
Philipp Arras's avatar
Philipp Arras committed
86

87
    This is then exponentiated and exponentially binned (in this order).
88 89 90 91 92 93 94

    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`).
Martin Reinecke's avatar
Martin Reinecke committed
95 96 97

    Parameters
    ----------
98 99 100 101 102 103 104 105
    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`).
Philipp Arras's avatar
Philipp Arras committed
106
    sm : float
107
        Expected exponent of power law. FIXME
Philipp Arras's avatar
Philipp Arras committed
108
    sv : float
109
        Prior variance of exponent of power law.
Philipp Arras's avatar
Philipp Arras committed
110
    im : float
111
        Expected y-intercept of power law. FIXME
Philipp Arras's avatar
Philipp Arras committed
112
    iv : float
113
        Prior variance of y-intercept of power law.
114 115 116 117 118 119 120

    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.
Martin Reinecke's avatar
Martin Reinecke committed
121
    '''
122
    from ..operators.exp_transform import ExpTransform
123 124
    from ..operators.slope_operator import SlopeOperator
    from ..operators.offset_operator import OffsetOperator
Martin Reinecke's avatar
Martin Reinecke committed
125

126 127 128 129 130
    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)
131 132
    if sv <= 0 or iv <= 0:
        raise ValueError
133 134 135 136

    et = ExpTransform(target, n_pix)
    dom = et.domain[0]

137
    # Smooth component
138 139
    dct = {'a': a, 'k0': k0}
    smooth = CepstrumOperator(dom, **dct).ducktape(keys[0])
Martin Reinecke's avatar
Martin Reinecke committed
140

141 142 143 144 145 146 147 148 149 150
    # Linear component
    phi_mean = np.array([sm, im + sm*dom.t_0[0]])
    phi_sig = np.array([sv, iv])
    slope = SlopeOperator(dom)
    phi_mean = Field.from_global_data(slope.domain, phi_mean)
    phi_sig = Field.from_global_data(slope.domain, phi_sig)
    linear = slope(OffsetOperator(phi_mean)(makeOp(phi_sig))).ducktape(keys[1])

    # Combine linear and smooth component
    loglog_ampl = 0.5*(smooth + linear)
Philipp Arras's avatar
Changes  
Philipp Arras committed
151

152 153
    # Go from loglog-space to linear-linear-space
    return et(loglog_ampl.exp())