correlated_fields.py 4.91 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', codomain=None):
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
    codomain : Domain
        The codomain for target[0]. If not supplied, it is inferred.
48 49 50

    Returns
    -------
Philipp Arras's avatar
Docs  
Philipp Arras committed
51 52
    Operator
        Correlated field
Philipp Arras's avatar
Docs  
Philipp Arras committed
53

Martin Reinecke's avatar
Martin Reinecke committed
54 55 56 57 58 59
    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`.
60
    """
61 62 63
    tgt = DomainTuple.make(target)
    if len(tgt) > 1:
        raise ValueError
64 65 66
    if codomain is None:
        codomain = tgt[0].get_default_codomain()
    h_space = codomain
67
    ht = HarmonicTransformOperator(h_space, target=tgt[0])
Philipp Arras's avatar
Philipp Arras committed
68
    p_space = amplitude_operator.target[0]
Martin Reinecke's avatar
Martin Reinecke committed
69
    power_distributor = PowerDistributor(h_space, p_space)
Philipp Arras's avatar
Philipp Arras committed
70
    A = power_distributor(amplitude_operator)
71
    vol = h_space.scalar_dvol**-0.5
Philipp Arras's avatar
Comment  
Philipp Arras committed
72 73 74 75
    # 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.
76
    return ht(vol*A*ducktape(h_space, None, name))
77 78


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

Martin Reinecke's avatar
Martin Reinecke committed
84
    This operator may be used as a model for multi-frequency reconstructions
85 86 87 88 89
    with a correlation structure in both spatial and energy direction.

    Parameters
    ----------
    target : Domain, DomainTuple or tuple of Domain
Philipp Arras's avatar
Docs  
Philipp Arras committed
90
        Target of the operator. Must contain exactly two spaces.
91 92 93 94 95 96 97
    amplitudes: iterable of Operator
        List of two amplitude operators.
    name : string
        :class:`MultiField` key for xi-field.

    Returns
    -------
Philipp Arras's avatar
Docs  
Philipp Arras committed
98 99
    Operator
        Correlated field
Philipp Arras's avatar
Docs  
Philipp Arras committed
100

Martin Reinecke's avatar
Martin Reinecke committed
101 102 103 104 105 106 107
    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)
108
    """
109 110 111 112 113
    tgt = DomainTuple.make(target)
    if len(tgt) != 2:
        raise ValueError
    if len(amplitudes) != 2:
        raise ValueError
114

115 116
    hsp = DomainTuple.make([tt.get_default_codomain() for tt in tgt])
    ht1 = HarmonicTransformOperator(hsp, target=tgt[0], space=0)
117
    ht2 = HarmonicTransformOperator(ht1.target, target=tgt[1], space=1)
118
    ht = ht2 @ ht1
119

120 121 122 123
    psp = [aa.target[0] for aa in amplitudes]
    pd0 = PowerDistributor(hsp, psp[0], 0)
    pd1 = PowerDistributor(pd0.domain, psp[1], 1)
    pd = pd0 @ pd1
124

125 126 127
    dd0 = ContractionOperator(pd.domain, 1).adjoint
    dd1 = ContractionOperator(pd.domain, 0).adjoint
    d = [dd0, dd1]
128

129
    a = [dd @ amplitudes[ii] for ii, dd in enumerate(d)]
Philipp Arras's avatar
Philipp Arras committed
130
    a = reduce(mul, a)
131
    A = pd @ a
Philipp Arras's avatar
Comment  
Philipp Arras committed
132
    # For `vol` see comment in `CorrelatedField`
Philipp Arras's avatar
Philipp Arras committed
133
    vol = reduce(mul, [sp.scalar_dvol**-0.5 for sp in hsp])
134
    return ht(vol*A*ducktape(hsp, None, name))