correlated_fields.py 4.03 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
Martin Reinecke's avatar
Martin Reinecke committed
25
from ..operators.simple_linear_operators import ducktape
Martin Reinecke's avatar
Martin Reinecke committed
26

Martin Reinecke's avatar
Martin Reinecke committed
27

28
def CorrelatedField(target, amplitude_operator, name='xi'):
29
    """Constructs an operator which turns a white Gaussian excitation field
Martin Reinecke's avatar
Martin Reinecke committed
30
    into a correlated field.
31 32 33 34 35

    This function returns an operator which implements:

        ht @ (vol * A * xi),

36 37
    where `ht` is a harmonic transform operator, `A` is the square root of the
    prior covariance and `xi` is the excitation field.
38 39 40

    Parameters
    ----------
41
    target : Domain, DomainTuple or tuple of Domain
Martin Reinecke's avatar
Martin Reinecke committed
42
        Target of the operator. Must contain exactly one space.
Philipp Arras's avatar
Philipp Arras committed
43
    amplitude_operator: Operator
44
    name : string
45
        :class:`MultiField` key for the xi-field.
46 47 48 49

    Returns
    -------
    Correlated field : Operator
50
    """
51 52 53 54 55
    tgt = DomainTuple.make(target)
    if len(tgt) > 1:
        raise ValueError
    h_space = tgt[0].get_default_codomain()
    ht = HarmonicTransformOperator(h_space, tgt[0])
Philipp Arras's avatar
Philipp Arras committed
56
    p_space = amplitude_operator.target[0]
Martin Reinecke's avatar
Martin Reinecke committed
57
    power_distributor = PowerDistributor(h_space, p_space)
Philipp Arras's avatar
Philipp Arras committed
58
    A = power_distributor(amplitude_operator)
59
    vol = h_space.scalar_dvol**-0.5
Philipp Arras's avatar
Comment  
Philipp Arras committed
60 61 62 63
    # 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.
64
    return ht(vol*A*ducktape(h_space, None, name))
65 66


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

Martin Reinecke's avatar
Martin Reinecke committed
72
    This operator may be used as a model for multi-frequency reconstructions
73 74 75 76 77
    with a correlation structure in both spatial and energy direction.

    Parameters
    ----------
    target : Domain, DomainTuple or tuple of Domain
Martin Reinecke's avatar
Martin Reinecke committed
78
        Target of the operator. Must contain exactly one space.
79 80 81 82 83 84 85 86
    amplitudes: iterable of Operator
        List of two amplitude operators.
    name : string
        :class:`MultiField` key for xi-field.

    Returns
    -------
    Correlated field : Operator
87
    """
88 89 90 91 92
    tgt = DomainTuple.make(target)
    if len(tgt) != 2:
        raise ValueError
    if len(amplitudes) != 2:
        raise ValueError
93

94 95 96 97
    hsp = DomainTuple.make([tt.get_default_codomain() for tt in tgt])
    ht1 = HarmonicTransformOperator(hsp, target=tgt[0], space=0)
    ht2 = HarmonicTransformOperator(ht1.target, space=1)
    ht = ht2 @ ht1
98

99 100 101 102
    psp = [aa.target[0] for aa in amplitudes]
    pd0 = PowerDistributor(hsp, psp[0], 0)
    pd1 = PowerDistributor(pd0.domain, psp[1], 1)
    pd = pd0 @ pd1
103

104 105 106
    dd0 = ContractionOperator(pd.domain, 1).adjoint
    dd1 = ContractionOperator(pd.domain, 0).adjoint
    d = [dd0, dd1]
107

108
    a = [dd @ amplitudes[ii] for ii, dd in enumerate(d)]
Philipp Arras's avatar
Philipp Arras committed
109
    a = reduce(mul, a)
110
    A = pd @ a
Philipp Arras's avatar
Comment  
Philipp Arras committed
111
    # For `vol` see comment in `CorrelatedField`
Philipp Arras's avatar
Philipp Arras committed
112
    vol = reduce(mul, [sp.scalar_dvol**-0.5 for sp in hsp])
113
    return ht(vol*A*ducktape(hsp, None, name))