amplitude_operator.py 4.75 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
Philipp Arras's avatar
Philipp Arras committed
22
from ..sugar import makeOp, sqrt
23
24
25


def _ceps_kernel(dof_space, k, a, k0):
26
    return a**2/(1 + (k/k0)**2)**2
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43


def create_cepstrum_amplitude_field(domain, cepstrum):
    """Creates a ...
    Writes the sum of all modes into the zero-mode.

    Parameters
    ----------
    domain: ???
        ???
    cepstrum: Callable
        ???
    """

    dim = len(domain.shape)
    shape = domain.shape

44
    q_array = domain.get_k_array()
45
46

    # Fill cepstrum field (all non-zero modes)
47
    no_zero_modes = (slice(1, None),)*dim
48
49
50
51
52
53
54
    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)
    for i in range(dim):
        # Prepare indices
Martin Reinecke's avatar
Martin Reinecke committed
55
56
57
        fst_dims = (slice(None),)*i
        sl = fst_dims + (slice(1, None),)
        sl2 = fst_dims + (0,)
58
59
60
61

        # Do summation
        cepstrum_field[sl2] = np.sum(cepstrum_field[sl], axis=i)

62
    return Field.from_global_data(domain, cepstrum_field)
Martin Reinecke's avatar
Martin Reinecke committed
63

Martin Reinecke's avatar
Martin Reinecke committed
64

Philipp Arras's avatar
Philipp Arras committed
65
def _CepstrumOperator(logk_space, ceps_a, ceps_k, zero_mode=True):
66
67
68
69
70
71
72
73
    '''
    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
    '''

Philipp Arras's avatar
Philipp Arras committed
74
75
76
77
78
79
80
    from ..operators.qht_operator import QHTOperator
    from ..operators.symmetrizing_operator import SymmetrizingOperator
    qht = QHTOperator(target=logk_space)
    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)
81
82
83
84
85
86
87
88
    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
89
90


Philipp Arras's avatar
Philipp Arras committed
91
def _SlopePowerSpectrum(logk_space, sm, sv, im, iv):
92
93
94
95
96
97
98
    '''
    Parameters
    ----------

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

99
    im, iv : y-intercept_mean, y-intercept_std  of power_slope
100
101
    '''

102
    from ..operators.slope_operator import SlopeOperator
Philipp Arras's avatar
Cleanup    
Philipp Arras committed
103
    from ..operators.offset_operator import OffsetOperator
104
105
106
107
108
    phi_mean = np.array([sm, im + sm*logk_space.t_0[0]])
    phi_sig = np.array([sv, iv])
    slope = SlopeOperator(logk_space)
    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
109
    return slope(OffsetOperator(phi_mean)(makeOp(phi_sig)))
110
111


Philipp Arras's avatar
Philipp Arras committed
112
113
def AmplitudeOperator(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
                      keys=['tau', 'phi'], zero_mode=True):
Philipp Arras's avatar
Philipp Arras committed
114
115
    ''' Operator for parametrizing smooth power spectra.

Martin Reinecke's avatar
Martin Reinecke committed
116
    Computes a smooth power spectrum.
Philipp Arras's avatar
Philipp Arras committed
117
    Output is defined on a PowerSpace.
Martin Reinecke's avatar
Martin Reinecke committed
118
119
120

    Parameters
    ----------
Philipp Arras's avatar
Philipp Arras committed
121
122
123
124
125
126
127
128
129
130
131
132
133
134
    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
    sm : float
        slope_mean = expected exponent of power law (e.g. -4)
    sv : float
        slope_variance (default=1)
    im : float
        y-intercept_mean
    iv : float
        y-intercept_variance  of power_slope
Martin Reinecke's avatar
Martin Reinecke committed
135
136
    '''

137
138
    from ..operators.exp_transform import ExpTransform
    from ..operators.scaling_operator import ScalingOperator
Martin Reinecke's avatar
Martin Reinecke committed
139

140
141
142
    h_space = s_space.get_default_codomain()
    et = ExpTransform(PowerSpace(h_space), Npixdof)
    logk_space = et.domain[0]
Martin Reinecke's avatar
Martin Reinecke committed
143

Philipp Arras's avatar
Philipp Arras committed
144
    smooth = _CepstrumOperator(logk_space, ceps_a, ceps_k, zero_mode)
Martin Reinecke's avatar
Martin Reinecke committed
145
    smooth = smooth.ducktape(keys[0])
Philipp Arras's avatar
Philipp Arras committed
146
    linear = _SlopePowerSpectrum(logk_space, sm, sv, im, iv)
Martin Reinecke's avatar
Martin Reinecke committed
147
    linear = linear.ducktape(keys[1])
Philipp Arras's avatar
Changes    
Philipp Arras committed
148

149
    fac = ScalingOperator(0.5, smooth.target)
Martin Reinecke's avatar
Martin Reinecke committed
150
    return et((fac(smooth + linear)).exp())