amplitude_operator.py 5.5 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
35
36
    # Fill all non-zero modes
    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):
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
57
58

    # FIXME a>0 k0>0
    qht = QHTOperator(target=domain)
Philipp Arras's avatar
Philipp Arras committed
59
    dof_space = qht.domain[0]
60
61
62
63
64
65
66
    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):
67
68
69
70
71
72
73
    '''
    Parameters
    ----------

    sm, sv : slope_mean = expected exponent of power law (e.g. -4),
                slope_variance (default=1)

74
    im, iv : y-intercept_mean, y-intercept_std  of power_slope
75

76
    '''
77
    from ..operators.slope_operator import SlopeOperator
Philipp Arras's avatar
Cleanup    
Philipp Arras committed
78
    from ..operators.offset_operator import OffsetOperator
79
80
81
82

    # sv, iv>0

    phi_mean = np.array([sm, im + sm*domain.t_0[0]])
83
    phi_sig = np.array([sv, iv])
84
    slope = SlopeOperator(domain)
85
86
    phi_mean = Field.from_global_data(slope.domain, phi_mean)
    phi_sig = Field.from_global_data(slope.domain, phi_sig)
Philipp Arras's avatar
Philipp Arras committed
87
    return slope(OffsetOperator(phi_mean)(makeOp(phi_sig)))
88
89


90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
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)
Philipp Arras's avatar
Philipp Arras committed
108

109
110
111
112
113
114
115
116
    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`).
Martin Reinecke's avatar
Martin Reinecke committed
117
118
119

    Parameters
    ----------
120
121
122
123
124
125
126
127
    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
128
    sm : float
129
        Expected exponent of power law (see :class:`SlopeOperator`).
Philipp Arras's avatar
Philipp Arras committed
130
    sv : float
131
        Prior variance of exponent of power law (see :class:`SlopeOperator`).
Philipp Arras's avatar
Philipp Arras committed
132
    im : float
133
        Expected y-intercept of power law (see :class:`SlopeOperator`).
Philipp Arras's avatar
Philipp Arras committed
134
    iv : float
135
136
137
138
139
140
141
142
        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.
Martin Reinecke's avatar
Martin Reinecke committed
143
    '''
144
    from ..operators.exp_transform import ExpTransform
Martin Reinecke's avatar
Martin Reinecke committed
145

146
147
148
149
150
151
152
153
154
155
156
    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])
Martin Reinecke's avatar
Martin Reinecke committed
157

158
159
    dct = {'sm': sm, 'sv': sv, 'im': im, 'iv': iv}
    linear = SlopeOperator(dom, **dct).ducktape(keys[1])
Philipp Arras's avatar
Changes    
Philipp Arras committed
160

161
    return et((0.5*(smooth + linear)).exp())