correlated_fields.py 7.09 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
from functools import reduce
Philipp Arras's avatar
Philipp Arras committed
19
from operator import mul
20

Philipp Arras's avatar
Philipp Arras committed
21
from ..domain_tuple import DomainTuple
22
from ..operators.contraction_operator import ContractionOperator
Philipp Arras's avatar
Philipp Arras committed
23
from ..operators.distributors import PowerDistributor
Martin Reinecke's avatar
Martin Reinecke committed
24
from ..operators.harmonic_operators import HarmonicTransformOperator
25
from ..operators.simple_linear_operators import ducktape
26

27 28 29 30 31 32
# import numpy as np
# from ..operators.value_inserter import ValueInserter
# from ..operators.simple_linear_operators import FieldAdapter
# from ..operators.diagonal_operator import DiagonalOperator
# from ..sugar import from_global_data
# from ..library.inverse_gamma_operator import InverseGammaOperator
Martin Reinecke's avatar
Martin Reinecke committed
33

34 35

def CorrelatedField(target, amplitude_operator, name='xi', codomain=None):
36
    """Constructs an operator which turns a white Gaussian excitation field
Martin Reinecke's avatar
Martin Reinecke committed
37
    into a correlated field.
38 39 40 41 42

    This function returns an operator which implements:

        ht @ (vol * A * xi),

43 44
    where `ht` is a harmonic transform operator, `A` is the square root of the
    prior covariance and `xi` is the excitation field.
45 46 47

    Parameters
    ----------
48
    target : Domain, DomainTuple or tuple of Domain
Martin Reinecke's avatar
Martin Reinecke committed
49
        Target of the operator. Must contain exactly one space.
Philipp Arras's avatar
Philipp Arras committed
50
    amplitude_operator: Operator
51
    name : string
52
        :class:`MultiField` key for the xi-field.
53 54
    codomain : Domain
        The codomain for target[0]. If not supplied, it is inferred.
55 56 57

    Returns
    -------
Philipp Arras's avatar
Philipp Arras committed
58 59
    Operator
        Correlated field
Philipp Arras's avatar
Philipp Arras committed
60

Martin Reinecke's avatar
Martin Reinecke committed
61 62 63 64 65 66
    Notes
    -----
    In NIFTy, non-harmonic RGSpaces are by definition periodic. Therefore
    the operator constructed by this method will output a correlated field
    with *periodic* boundary conditions. If a non-periodic field is needed,
    one needs to combine this operator with a :class:`FieldZeroPadder`.
67
    """
68 69 70
    tgt = DomainTuple.make(target)
    if len(tgt) > 1:
        raise ValueError
71 72 73
    if codomain is None:
        codomain = tgt[0].get_default_codomain()
    h_space = codomain
74
    ht = HarmonicTransformOperator(h_space, target=tgt[0])
Philipp Arras's avatar
Philipp Arras committed
75
    p_space = amplitude_operator.target[0]
Martin Reinecke's avatar
Martin Reinecke committed
76
    power_distributor = PowerDistributor(h_space, p_space)
Philipp Arras's avatar
Philipp Arras committed
77
    A = power_distributor(amplitude_operator)
78
    vol = h_space.scalar_dvol**-0.5
79

Philipp Arras's avatar
Philipp Arras committed
80 81 82 83
    # When doubling the resolution of `tgt` the value of the highest k-mode
    # will scale with a square root. `vol` cancels this effect such that the
    # same power spectrum can be used for the spaces with the same volume,
    # different resolutions and the same object in them.
84
    return ht(vol*(A)*ducktape(h_space, None, name))
85 86


87
def MfCorrelatedField(target, amplitudes, name='xi'):
88
    """Constructs an operator which turns white Gaussian excitation fields
Martin Reinecke's avatar
Martin Reinecke committed
89 90
    into a correlated field defined on a DomainTuple with two entries and two
    separate correlation structures.
91

Martin Reinecke's avatar
Martin Reinecke committed
92
    This operator may be used as a model for multi-frequency reconstructions
93 94 95 96 97
    with a correlation structure in both spatial and energy direction.

    Parameters
    ----------
    target : Domain, DomainTuple or tuple of Domain
Philipp Arras's avatar
Philipp Arras committed
98
        Target of the operator. Must contain exactly two spaces.
99 100 101 102 103 104 105
    amplitudes: iterable of Operator
        List of two amplitude operators.
    name : string
        :class:`MultiField` key for xi-field.

    Returns
    -------
Philipp Arras's avatar
Philipp Arras committed
106 107
    Operator
        Correlated field
Philipp Arras's avatar
Philipp Arras committed
108

Martin Reinecke's avatar
Martin Reinecke committed
109 110 111 112 113 114 115
    Notes
    -----
    In NIFTy, non-harmonic RGSpaces are by definition periodic. Therefore
    the operator constructed by this method will output a correlated field
    with *periodic* boundary conditions. If a non-periodic field is needed,
    one needs to combine this operator with a :class:`FieldZeroPadder` or even
    two (one for the energy and one for the spatial subdomain)
116
    """
117 118 119 120 121
    tgt = DomainTuple.make(target)
    if len(tgt) != 2:
        raise ValueError
    if len(amplitudes) != 2:
        raise ValueError
122

123 124
    hsp = DomainTuple.make([tt.get_default_codomain() for tt in tgt])
    ht1 = HarmonicTransformOperator(hsp, target=tgt[0], space=0)
125
    ht2 = HarmonicTransformOperator(ht1.target, target=tgt[1], space=1)
126
    ht = ht2 @ ht1
127

128 129 130 131
    psp = [aa.target[0] for aa in amplitudes]
    pd0 = PowerDistributor(hsp, psp[0], 0)
    pd1 = PowerDistributor(pd0.domain, psp[1], 1)
    pd = pd0 @ pd1
132

133 134 135
    dd0 = ContractionOperator(pd.domain, 1).adjoint
    dd1 = ContractionOperator(pd.domain, 0).adjoint
    d = [dd0, dd1]
136

137
    a = [dd @ amplitudes[ii] for ii, dd in enumerate(d)]
Philipp Arras's avatar
Philipp Arras committed
138
    a = reduce(mul, a)
139
    A = pd @ a
Philipp Arras's avatar
Philipp Arras committed
140
    # For `vol` see comment in `CorrelatedField`
Philipp Arras's avatar
Philipp Arras committed
141
    vol = reduce(mul, [sp.scalar_dvol**-0.5 for sp in hsp])
142

143
    return ht(vol*A*ducktape(hsp, None, name))
144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194


def MfPartiallyCorrelatedField(target, energy_amplitude, name='xi_u'):
    """Constructs an operator which turns white Gaussian excitation fields
    into a correlated field defined on a DomainTuple with two entries.
    On the first domain, a white correlation structure is assumed. On
    the second domain, the correlation structures given by energy_amplitude
    is used.

    This operator may be used as a model for multi-frequency reconstructions
    with correlation structure only in the energy direction.

    Parameters
    ----------
    target : Domain, DomainTuple or tuple of Domain
        Target of the operator. Must contain exactly two spaces.
        It is assumed that the second space is the energy domain.
    energy_amplitude: Operator
        amplitude operator for the energy correlation structure
    name : string
        :class:`MultiField` key for xi-field.

    Returns
    -------
    Operator
        Correlated field

    Notes
    -----
    In NIFTy, non-harmonic RGSpaces are by definition periodic. Therefore
    the operator constructed by this method will output a correlated field
    with *periodic* boundary conditions. If a non-periodic field is needed,
    one needs to combine this operator with a :class:`FieldZeroPadder` or even
    two (one for the energy and one for the spatial subdomain)
    """
    tgt = DomainTuple.make(target)
    if len(tgt) != 2:
        raise ValueError

    h_space = DomainTuple.make([dom.get_default_codomain() for dom in tgt])
    ht1 = HarmonicTransformOperator(h_space, target=tgt[0], space=0)
    ht2 = HarmonicTransformOperator(ht1.target, target=tgt[1], space=1)
    ht = ht2 @ ht1

    p_space = energy_amplitude.target[0]
    power_distributor = PowerDistributor(h_space[-1], p_space)
    A = power_distributor(energy_amplitude)

    dd = ContractionOperator(h_space, 0).adjoint

    return ht((dd @ A)*ducktape(h_space, None, name))