correlated_fields.py 8.1 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
87
88
    if isinstance(amplitude_operator, Field):
        return ht(makeOp(A)@xi).scale(vol)
    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
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
216
217
218
219
220
221
222
223
224


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]))