correlated_fields.py 3.89 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 19
from functools import reduce

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

Martin Reinecke's avatar
Martin Reinecke committed
26

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

    This function returns an operator which implements:

        ht @ (vol * A * xi),

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

    Parameters
    ----------
40
    target : Domain, DomainTuple or tuple of Domain
Martin Reinecke's avatar
Martin Reinecke committed
41
        Target of the operator. Must contain exactly one space.
Philipp Arras's avatar
Philipp Arras committed
42
    amplitude_operator: Operator
43
    name : string
44
        :class:`MultiField` key for the xi-field.
45 46
    codomain : Domain
        The codomain for target[0]. If not supplied, it is inferred.
47 48 49

    Returns
    -------
Philipp Arras's avatar
Docs  
Philipp Arras committed
50 51
    Operator
        Correlated field
52
    """
53 54 55
    tgt = DomainTuple.make(target)
    if len(tgt) > 1:
        raise ValueError
56 57 58
    if codomain is None:
        codomain = tgt[0].get_default_codomain()
    h_space = codomain
59
    ht = HarmonicTransformOperator(h_space, target=tgt[0])
Philipp Arras's avatar
Philipp Arras committed
60
    p_space = amplitude_operator.target[0]
Martin Reinecke's avatar
Martin Reinecke committed
61
    power_distributor = PowerDistributor(h_space, p_space)
Philipp Arras's avatar
Philipp Arras committed
62
    A = power_distributor(amplitude_operator)
63 64
    vol = h_space.scalar_dvol**-0.5
    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
Philipp Arras's avatar
Docs  
Philipp Arras committed
78
        Target of the operator. Must contain exactly two spaces.
79 80 81 82 83 84 85
    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
86 87
    Operator
        Correlated field
88
    """
89 90 91 92 93
    tgt = DomainTuple.make(target)
    if len(tgt) != 2:
        raise ValueError
    if len(amplitudes) != 2:
        raise ValueError
94

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

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

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

109 110 111 112 113
    a = [dd @ amplitudes[ii] for ii, dd in enumerate(d)]
    a = reduce(lambda x, y: x*y, a)
    A = pd @ a
    vol = reduce(lambda x, y: x*y, [sp.scalar_dvol**-0.5 for sp in hsp])
    return ht(vol*A*ducktape(hsp, None, name))