smooth_linear_amplitude.py 6.42 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 24 25 26 27
from ..operators.exp_transform import ExpTransform
from ..operators.offset_operator import OffsetOperator
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(dof_space, k, a, k0):
32
    return a**2/(1 + (k/k0)**2)**2
33 34


Philipp Arras's avatar
Docs  
Philipp Arras committed
35 36
def CepstrumOperator(target, a, k0):
    '''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 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85
    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::
        \\mathrm{sqrt\_ceps}(k) = \\frac{a}{1+(k/k0)^2}

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

    FIXME The prior on the zero mode is ...

    Parameters
    ----------
    target : LogRGSpace
        Target domain of the operator, needs to be non-harmonic and
        one-dimensional.
    a : float
        Strength of smoothness prior (positive only). FIXME
    k0 : float
        Cutoff of smothness prior in quefrency space (positive only). FIXME
    '''
    a, k0 = float(a), float(k0)
    target = DomainTuple.make(target)
    if a <= 0 or k0 <= 0:
        raise ValueError
    if len(target) > 1 or target[0].harmonic or len(target[0].shape) > 1:
        raise TypeError

    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()
86
    # Fill all non-zero modes
Philipp Arras's avatar
Philipp Arras committed
87 88
    no_zero_modes = (slice(1, None),)*dim
    ks = q_array[(slice(None),) + no_zero_modes]
89
    cepstrum_field = np.zeros(shape)
Philipp Arras's avatar
Docs  
Philipp Arras committed
90
    cepstrum_field[no_zero_modes] = _ceps_kernel(dom, ks, a, k0)
91
    # Fill zero-mode subspaces
92
    for i in range(dim):
Philipp Arras's avatar
Philipp Arras committed
93 94 95
        fst_dims = (slice(None),)*i
        sl = fst_dims + (slice(1, None),)
        sl2 = fst_dims + (0,)
96
        cepstrum_field[sl2] = np.sum(cepstrum_field[sl], axis=i)
Philipp Arras's avatar
Docs  
Philipp Arras committed
97
    cepstrum = Field.from_global_data(dom, cepstrum_field)
Philipp Arras's avatar
Philipp Arras committed
98

99 100 101
    return sym @ qht @ makeOp(cepstrum.sqrt())


102 103 104
def SLAmplitude(target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
    '''Operator for parametrizing smooth amplitudes (square roots of power
    spectra).
105 106 107

    The general guideline for setting up generative models in IFT is to
    transform the problem into the eigenbase of the prior and formulate the
108 109
    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
110 111 112
    double-logarithmic scale).

    This function assembles an :class:`Operator` which maps two a-priori white
113
    Gaussian random fields to a smooth amplitude which is composed out of
114 115 116 117 118 119
    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
120

121
    This is then exponentiated and exponentially binned (in this order).
122 123 124 125 126 127 128

    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
129 130 131

    Parameters
    ----------
132 133 134 135 136 137 138 139
    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
140
    sm : float
141
        Expected exponent of power law. FIXME
Philipp Arras's avatar
Philipp Arras committed
142
    sv : float
143
        Prior standard deviation of exponent of power law.
Philipp Arras's avatar
Philipp Arras committed
144
    im : float
145
        Expected y-intercept of power law. FIXME
Philipp Arras's avatar
Philipp Arras committed
146
    iv : float
147
        Prior standard deviation of y-intercept of power law.
148 149 150 151 152 153 154

    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
155
    '''
156 157 158 159 160
    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)
161 162
    if sv <= 0 or iv <= 0:
        raise ValueError
163 164 165 166

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

167
    # Smooth component
168 169
    dct = {'a': a, 'k0': k0}
    smooth = CepstrumOperator(dom, **dct).ducktape(keys[0])
Martin Reinecke's avatar
Martin Reinecke committed
170

171
    # Linear component
172 173 174 175 176 177
    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)
    linear = (sl @ OffsetOperator(mean) @ makeOp(sig)).ducktape(keys[1])
178 179 180

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

182
    # Go from loglog-space to linear-linear-space
Philipp Arras's avatar
Philipp Arras committed
183
    return et @ loglog_ampl.exp()