smooth_linear_amplitude.py 8.23 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

Philipp Arras's avatar
Docs  
Philipp Arras committed
20
from ..domain_tuple import DomainTuple
Martin Reinecke's avatar
Martin Reinecke committed
21
from ..domains.power_space import PowerSpace
22
from ..field import Field
Philipp Arras's avatar
Philipp Arras committed
23
from ..operators.adder import Adder
Philipp Arras's avatar
Philipp Arras committed
24
25
26
27
from ..operators.exp_transform import ExpTransform
from ..operators.qht_operator import QHTOperator
from ..operators.slope_operator import SlopeOperator
from ..operators.symmetrizing_operator import SymmetrizingOperator
28
from ..sugar import makeOp
29
30
31
from ..library.inverse_gamma_operator import InverseGammaOperator
from ..operators.value_inserter import ValueInserter

32

33
def _ceps_kernel(k, a, k0):
Philipp Arras's avatar
Philipp Arras committed
34
    return (a/(1 + np.sum((k.T/k0)**2, axis=-1).T))**2
35
36


Philipp Arras's avatar
Docs  
Philipp Arras committed
37
def CepstrumOperator(target, a, k0):
38
    """Turns a white Gaussian random field into a smooth field on a LogRGSpace.
39

Philipp Arras's avatar
Docs  
Philipp Arras committed
40
41
42
43
44
45
46
47
    Composed out of three operators:

        sym @ qht @ diag(sqrt_ceps),

    where sym is a :class:`SymmetrizingOperator`, qht is a :class:`QHTOperator`
    and ceps is the so-called cepstrum:

    .. math::
Martin Reinecke's avatar
Martin Reinecke committed
48
        \\mathrm{sqrt\\_ceps}(k) = \\frac{a}{1+(k/k0)^2}
Philipp Arras's avatar
Docs  
Philipp Arras committed
49
50
51
52
53
54
55

    These operators are combined in this fashion in order to generate:

    - A field which is smooth, i.e. second derivatives are punished (note
      that the sqrt-cepstrum is essentially proportional to 1/k**2).

    - A field which is symmetric around the pixel in the middle of the space.
Philipp Arras's avatar
Docs  
Philipp Arras committed
56
57
      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
Philipp Arras's avatar
Docs  
Philipp Arras committed
58
59
60
      amplitude whenever :class:`CepstrumOperator` is used as in
      :class:`SLAmplitude`.

Philipp Arras's avatar
Docs  
Philipp Arras committed
61
62
    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
63
    axis.
Philipp Arras's avatar
Docs  
Philipp Arras committed
64
65
66
67

    Parameters
    ----------
    target : LogRGSpace
68
        Target domain of the operator, needs to be non-harmonic.
Philipp Arras's avatar
Docs  
Philipp Arras committed
69
    a : float
70
71
72
73
        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
74
        Strength of smoothness prior in quefrency space (positive only) along
75
76
        each axis. If float then the strength is the same along each axis.
        Larger values result in a weaker constraining prior.
77
    """
78
    a = float(a)
Philipp Arras's avatar
Docs  
Philipp Arras committed
79
    target = DomainTuple.make(target)
80
    if a <= 0:
Philipp Arras's avatar
Docs  
Philipp Arras committed
81
        raise ValueError
82
    if len(target) > 1 or target[0].harmonic:
Philipp Arras's avatar
Docs  
Philipp Arras committed
83
        raise TypeError
Philipp Arras's avatar
Docs  
Philipp Arras committed
84
85
86
87
88
    if isinstance(k0, (float, int)):
        k0 = np.array([k0]*len(target.shape))
    else:
        k0 = np.array(k0)
    if len(k0) != len(target.shape):
89
90
91
        raise ValueError
    if np.any(np.array(k0) <= 0):
        raise ValueError
Philipp Arras's avatar
Docs  
Philipp Arras committed
92
93
94
95
96
97
98
99
100

    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()
101
    # Fill all non-zero modes
Philipp Arras's avatar
Philipp Arras committed
102
103
    no_zero_modes = (slice(1, None),)*dim
    ks = q_array[(slice(None),) + no_zero_modes]
104
    cepstrum_field = np.zeros(shape)
Philipp Arras's avatar
Docs  
Philipp Arras committed
105
    cepstrum_field[no_zero_modes] = _ceps_kernel(ks, a, k0)
106
    # Fill zero-mode subspaces
107
    for i in range(dim):
Philipp Arras's avatar
Philipp Arras committed
108
109
110
        fst_dims = (slice(None),)*i
        sl = fst_dims + (slice(1, None),)
        sl2 = fst_dims + (0,)
111
        cepstrum_field[sl2] = np.sum(cepstrum_field[sl], axis=i)
Philipp Arras's avatar
Docs  
Philipp Arras committed
112
    cepstrum = Field.from_global_data(dom, cepstrum_field)
Philipp Arras's avatar
Philipp Arras committed
113

114
115
116
    return sym @ qht @ makeOp(cepstrum.sqrt())


117
def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, za=None, zq=None,
118
                keys=['tau', 'phi', 'zero_mode']):
119
120
    '''Operator for parametrizing smooth amplitudes (square roots of power
    spectra).
121
122
123

    The general guideline for setting up generative models in IFT is to
    transform the problem into the eigenbase of the prior and formulate the
124
125
    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
126
127
128
    double-logarithmic scale).

    This function assembles an :class:`Operator` which maps two a-priori white
129
    Gaussian random fields to a smooth amplitude which is composed out of
130
131
132
133
134
135
    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
136

137
    This is then exponentiated and exponentially binned (in this order).
138
139
140
141
142
143

    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
Martin Reinecke's avatar
Martin Reinecke committed
144
145
    strength and the cutoff of the smoothness prior
    (see :class:`CepstrumOperator`).
Martin Reinecke's avatar
Martin Reinecke committed
146
147
148

    Parameters
    ----------
149
150
151
152
153
154
155
    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
Philipp Arras's avatar
Docs  
Philipp Arras committed
156
157
        Cutoff of smothness prior in quefrency space (see
        :class:`CepstrumOperator`).
Philipp Arras's avatar
Philipp Arras committed
158
    sm : float
Philipp Arras's avatar
Docs  
Philipp Arras committed
159
        Expected exponent of power law.
Philipp Arras's avatar
Philipp Arras committed
160
    sv : float
161
        Prior standard deviation of exponent of power law.
Philipp Arras's avatar
Philipp Arras committed
162
    im : float
Philipp Arras's avatar
Docs  
Philipp Arras committed
163
164
        Expected y-intercept of power law. This is the value at t_0 of the
        LogRGSpace (see :class:`ExpTransform`).
Philipp Arras's avatar
Philipp Arras committed
165
    iv : float
166
        Prior standard deviation of y-intercept of power law.
167
    za : float, optional
168
169
        Parameter of the optional zero mode prior (inverse-gamma): alpha
        See :class:`InverseGammaOperator` for interpretation.
170
    zq : float, optional
171
172
        Parameter of the optional zero mode prior (inverse-gamma): q
        See :class:`InverseGammaOperator` for interpretation.
173
174
175
176
177
178
179

    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
180
    '''
181
182
183
184
185
    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)
186
187
    if sv <= 0 or iv <= 0:
        raise ValueError
188

Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
189
190
    separate_zero_mode_prior = za is not None and zq is not None
    if separate_zero_mode_prior:
191
192
        za, zq = float(za), float(zq)
    else:
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
193
        if za is not None or zq is not None:
194
195
            raise ValueError("za and zq have to be given together")

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

199
    # Smooth component
200
201
    dct = {'a': a, 'k0': k0}
    smooth = CepstrumOperator(dom, **dct).ducktape(keys[0])
Martin Reinecke's avatar
Martin Reinecke committed
202

203
    # Linear component
204
205
206
207
208
    sl = SlopeOperator(dom)
    mean = np.array([sm, im + sm*dom.t_0[0]])
    sig = np.array([sv, iv])
    mean = Field.from_global_data(sl.domain, mean)
    sig = Field.from_global_data(sl.domain, sig)
Philipp Arras's avatar
Philipp Arras committed
209
    linear = sl @ Adder(mean) @ makeOp(sig).ducktape(keys[1])
210
211
212

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

214
215
216
217
    if not separate_zero_mode_prior:
        # Go from loglog-space to linear-linear-space
        return et @ loglog_ampl.exp()
    else:
Martin Reinecke's avatar
fix  
Martin Reinecke committed
218
        zero_mode = ValueInserter(et.target, (0,)*len(et.target.shape))
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
219
        zero_mode = (
Martin Reinecke's avatar
fix  
Martin Reinecke committed
220
            zero_mode @
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
221
            InverseGammaOperator(zero_mode.domain, za, zq).ducktape(keys[2]))
222
223
224

        mask = np.ones(et.target.shape)
        mask[(0,)*len(et.target.shape)] = 0.
Martin Reinecke's avatar
tweaks  
Martin Reinecke committed
225
        mask = makeOp(Field.from_global_data(et.target, mask))
226

227
        return mask @ et @ loglog_ampl.exp() + zero_mode