smooth_linear_amplitude.py 7.48 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
def _ceps_kernel(k, a, k0):
Philipp Arras's avatar
Philipp Arras committed
32
    return (a/(1 + np.sum((k.T/k0)**2, axis=-1).T))**2
33 34


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

Philipp Arras's avatar
Docs  
Philipp Arras committed
38 39 40 41 42 43 44 45
    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
46
        \\mathrm{sqrt\\_ceps}(k) = \\frac{a}{1+(k/k0)^2}
Philipp Arras's avatar
Docs  
Philipp Arras committed
47 48 49 50 51 52 53

    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
54 55
      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
56 57 58
      amplitude whenever :class:`CepstrumOperator` is used as in
      :class:`SLAmplitude`.

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

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

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

112 113
    return sym @ qht @ makeOp(cepstrum.sqrt())

Martin Reinecke's avatar
Martin Reinecke committed
114

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

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

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

134
    This is then exponentiated and exponentially binned (in this order).
135 136 137 138 139 140

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

    Parameters
    ----------
146 147 148 149 150 151 152
    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
153 154
        Cutoff of smothness prior in quefrency space (see
        :class:`CepstrumOperator`).
Philipp Arras's avatar
Philipp Arras committed
155
    sm : float
Philipp Arras's avatar
Docs  
Philipp Arras committed
156
        Expected exponent of power law.
Philipp Arras's avatar
Philipp Arras committed
157
    sv : float
158
        Prior standard deviation of exponent of power law.
Philipp Arras's avatar
Philipp Arras committed
159
    im : float
Philipp Arras's avatar
Docs  
Philipp Arras committed
160 161
        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
162
    iv : float
163
        Prior standard deviation of y-intercept of power law.
164 165 166 167 168 169 170

    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
171
    '''
Martin Reinecke's avatar
Martin Reinecke committed
172 173
    return LinearSLAmplitude(target=target, n_pix=n_pix, a=a, k0=k0, sm=sm,
                             sv=sv, im=im, iv=iv, keys=keys).exp()
174 175


Martin Reinecke's avatar
Martin Reinecke committed
176 177 178 179 180
def LinearSLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv,
                      keys=['tau', 'phi']):
    '''LinearOperator for parametrizing smooth log-amplitudes (square roots of
    power spectra).

181 182 183
    Logarithm of SLAmplitude
    See documentation of SLAmplitude for more details
    '''
184 185 186 187 188
    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)
189 190
    if sv <= 0 or iv <= 0:
        raise ValueError
191 192 193 194

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

195
    # Smooth component
196 197
    dct = {'a': a, 'k0': k0}
    smooth = CepstrumOperator(dom, **dct).ducktape(keys[0])
Martin Reinecke's avatar
Martin Reinecke committed
198

199
    # Linear component
200 201 202 203 204
    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
205
    linear = sl @ Adder(mean) @ makeOp(sig).ducktape(keys[1])
206 207 208

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

210
    # Go from loglog-space to linear-linear-space
211
    return et @ loglog_ampl