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

Martin Reinecke's avatar
Martin Reinecke committed
31

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

    This function returns an operator which implements:

        ht @ (vol * A * xi),

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

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

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

Martin Reinecke's avatar
Martin Reinecke committed
58 59 60 61 62 63
    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`.
64
    """
65 66 67
    tgt = DomainTuple.make(target)
    if len(tgt) > 1:
        raise ValueError
68 69 70
    if codomain is None:
        codomain = tgt[0].get_default_codomain()
    h_space = codomain
71
    ht = HarmonicTransformOperator(h_space, target=tgt[0])
72 73 74 75 76 77
    if isinstance(amplitude_operator, Operator):
        p_space = amplitude_operator.target[0]
    elif isinstance(amplitude_operator, Field):
        p_space = amplitude_operator.domain[0]
    else:
        raise TypeError
Martin Reinecke's avatar
Martin Reinecke committed
78
    power_distributor = PowerDistributor(h_space, p_space)
Philipp Arras's avatar
Philipp Arras committed
79
    A = power_distributor(amplitude_operator)
80
    vol = h_space.scalar_dvol**-0.5
81
    xi = ducktape(h_space, None, name)
Philipp Arras's avatar
Comment  
Philipp Arras committed
82 83 84 85
    # 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.
86
    if isinstance(amplitude_operator, Field):
Philipp Arras's avatar
Philipp Arras committed
87
        return ht(makeOp(A) @ xi).scale(vol)
88
    return ht(A*xi).scale(vol)
89 90


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

Martin Reinecke's avatar
Martin Reinecke committed
96
    This operator may be used as a model for multi-frequency reconstructions
97 98 99 100 101
    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
102
        Target of the operator. Must contain exactly two spaces.
103 104 105 106 107 108 109
    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
110 111
    Operator
        Correlated field
Philipp Arras's avatar
Docs  
Philipp Arras committed
112

Martin Reinecke's avatar
Martin Reinecke committed
113 114 115 116 117 118 119
    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)
120
    """
121 122 123 124 125
    tgt = DomainTuple.make(target)
    if len(tgt) != 2:
        raise ValueError
    if len(amplitudes) != 2:
        raise ValueError
126

127 128
    hsp = DomainTuple.make([tt.get_default_codomain() for tt in tgt])
    ht1 = HarmonicTransformOperator(hsp, target=tgt[0], space=0)
129
    ht2 = HarmonicTransformOperator(ht1.target, target=tgt[1], space=1)
130
    ht = ht2 @ ht1
131

132 133 134 135
    psp = [aa.target[0] for aa in amplitudes]
    pd0 = PowerDistributor(hsp, psp[0], 0)
    pd1 = PowerDistributor(pd0.domain, psp[1], 1)
    pd = pd0 @ pd1
136

137 138 139
    dd0 = ContractionOperator(pd.domain, 1).adjoint
    dd1 = ContractionOperator(pd.domain, 0).adjoint
    d = [dd0, dd1]
140

141
    a = [dd @ amplitudes[ii] for ii, dd in enumerate(d)]
Philipp Arras's avatar
Philipp Arras committed
142
    a = reduce(mul, a)
143
    A = pd @ a
Philipp Arras's avatar
Comment  
Philipp Arras committed
144
    # For `vol` see comment in `CorrelatedField`
Philipp Arras's avatar
Philipp Arras committed
145
    vol = reduce(mul, [sp.scalar_dvol**-0.5 for sp in hsp])
146
    return ht(vol*A*ducktape(hsp, None, name))
Philipp Arras's avatar
Philipp Arras committed
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


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,
Philipp Arras's avatar
Philipp Arras committed
190 191 192 193
    ] if isinstance(amplitudes, (Operator, Field)) else amplitudes

    cls = Operator if isinstance(amps[0], Operator) else Field
    assert all([isinstance(aa, cls) for aa in amps])
Philipp Arras's avatar
Philipp Arras committed
194 195 196 197 198 199 200 201

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

Philipp Arras's avatar
Philipp Arras committed
202 203 204 205
    if cls == Field:
        psp = [aa.domain[0] for aa in amps]
    else:
        psp = [aa.target[0] for aa in amps]
Philipp Arras's avatar
Philipp Arras committed
206 207 208 209 210 211 212 213 214 215 216
    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)))

Philipp Arras's avatar
Philipp Arras committed
217
    a = ContractionOperator(pd.domain, spaces[1:]).adjoint(amps[0])
Philipp Arras's avatar
Philipp Arras committed
218
    for i in range(1, len(amps)):
Philipp Arras's avatar
Philipp Arras committed
219 220
        a = a*(ContractionOperator(
            pd.domain, spaces[:i] + spaces[(i + 1):]).adjoint(amps[i]))
Philipp Arras's avatar
Philipp Arras committed
221

Philipp Arras's avatar
Philipp Arras committed
222
    expander = VdotOperator(full(pd.domain, 1.)).adjoint
Philipp Arras's avatar
Philipp Arras committed
223 224 225 226

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

Philipp Arras's avatar
Philipp Arras committed
227 228 229 230
    if cls == Field:
        A = pd @ (makeOp(a) @ Std)
    else:
        A = pd @ (Std*a)
Philipp Arras's avatar
Philipp Arras committed
231 232 233
    # 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]))