correlated_fields.py 7.79 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
Philipp Arras's avatar
Philipp Arras committed
22
from ..operators.adder import Adder
23
from ..operators.contraction_operator import ContractionOperator
Philipp Arras's avatar
Philipp Arras committed
24
from ..operators.distributors import PowerDistributor
Martin Reinecke's avatar
Martin Reinecke committed
25
from ..operators.harmonic_operators import HarmonicTransformOperator
Philipp Arras's avatar
Philipp Arras committed
26 27 28
from ..operators.operator import Operator
from ..operators.simple_linear_operators import VdotOperator, ducktape
from ..sugar import full
Martin Reinecke's avatar
Martin Reinecke committed
29

Martin Reinecke's avatar
Martin Reinecke committed
30

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

    This function returns an operator which implements:

        ht @ (vol * A * xi),

39 40
    where `ht` is a harmonic transform operator, `A` is the square root of the
    prior covariance and `xi` is the excitation field.
41 42 43

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

    Returns
    -------
Philipp Arras's avatar
Docs  
Philipp Arras committed
54 55
    Operator
        Correlated field
Philipp Arras's avatar
Docs  
Philipp Arras committed
56

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


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

Martin Reinecke's avatar
Martin Reinecke committed
87
    This operator may be used as a model for multi-frequency reconstructions
88 89 90 91 92
    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
93
        Target of the operator. Must contain exactly two spaces.
94 95 96 97 98 99 100
    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
101 102
    Operator
        Correlated field
Philipp Arras's avatar
Docs  
Philipp Arras committed
103

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

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

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

128 129 130
    dd0 = ContractionOperator(pd.domain, 1).adjoint
    dd1 = ContractionOperator(pd.domain, 0).adjoint
    d = [dd0, dd1]
131

132
    a = [dd @ amplitudes[ii] for ii, dd in enumerate(d)]
Philipp Arras's avatar
Philipp Arras committed
133
    a = reduce(mul, a)
134
    A = pd @ a
Philipp Arras's avatar
Comment  
Philipp Arras committed
135
    # For `vol` see comment in `CorrelatedField`
Philipp Arras's avatar
Philipp Arras committed
136
    vol = reduce(mul, [sp.scalar_dvol**-0.5 for sp in hsp])
137
    return ht(vol*A*ducktape(hsp, None, name))
Philipp Arras's avatar
Philipp Arras committed
138 139 140 141 142 143 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 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215


def CorrelatedFieldNormAmplitude(target,
                                 amplitudes,
                                 stdmean,
                                 stdstd,
                                 names=['xi', 'std']):
    """Constructs an operator which turns white Gaussian excitation fields
    into a correlated field defined on a DomainTuple with n entries and n
    separate correlation structures.

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

    Parameters
    ----------
    target : Domain, DomainTuple or tuple of Domain
        Target of the operator. Must contain exactly n spaces.
    amplitudes: Opertor, iterable of Operator
        Amplitude operator if n = 1 or list of n amplitude operators.
    stdmean : float
        Prior mean of the overall standart deviation.
    stdstd : float
        Prior standart deviation of the overall standart deviation.
    names : iterable of string
        :class:`MultiField` keys for xi-field and std-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)
    """

    amps = [
        amplitudes,
    ] if isinstance(amplitudes, Operator) else amplitudes

    tgt = DomainTuple.make(target)
    if len(tgt) != len(amps):
        raise ValueError
    stdmean, stdstd = float(stdmean), float(stdstd)
    if stdstd <= 0:
        raise ValueError

    psp = [aa.target[0] for aa in amps]
    hsp = DomainTuple.make([tt.get_default_codomain() for tt in tgt])

    ht = HarmonicTransformOperator(hsp, target=tgt[0], space=0)
    pd = PowerDistributor(hsp, psp[0], 0)

    for i in range(1, len(amps)):
        ht = HarmonicTransformOperator(ht.target, target=tgt[i], space=i) @ ht
        pd = pd @ PowerDistributor(pd.domain, psp[i], space=i)

    spaces = tuple(range(len(amps)))

    a = ContractionOperator(pd.domain, spaces[1:]).adjoint @ amps[0]
    for i in range(1, len(amps)):
        a = a*(ContractionOperator(pd.domain, spaces[:i] + spaces[
            (i + 1):]).adjoint @ amps[i])

    expander = VdotOperator(full(a.target, 1.)).adjoint

    Std = stdstd*ducktape(expander.domain, None, names[1])
    Std = expander @ (Adder(full(expander.domain, stdmean)) @ Std).exp()

    A = pd @ (Std*a)
    # For `vol` see comment in `CorrelatedField`
    vol = reduce(mul, [sp.scalar_dvol**-0.5 for sp in hsp])
    return ht(vol*A*ducktape(hsp, None, names[0]))