correlated_fields.py 28.4 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-2020 Max-Planck-Society
Philipp Arras's avatar
Philipp Arras committed
15
# Authors: Philipp Frank, Philipp Arras, Philipp Haim
Martin Reinecke's avatar
Martin Reinecke committed
16
#
17
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
Philipp Arras's avatar
Philipp Arras committed
18

19
20
21
from functools import reduce
from operator import mul

Philipp Arras's avatar
Philipp Arras committed
22
import numpy as np
23

Philipp Arras's avatar
Philipp Arras committed
24
from .. import utilities
Philipp Arras's avatar
Philipp Arras committed
25
from ..domain_tuple import DomainTuple
Philipp Arras's avatar
Philipp Arras committed
26
27
from ..domains.power_space import PowerSpace
from ..domains.unstructured_domain import UnstructuredDomain
Philipp Arras's avatar
Philipp Arras committed
28
from ..field import Field
29
from ..logger import logger
Philipp Arras's avatar
Philipp Arras committed
30
from ..multi_field import MultiField
Philipp Arras's avatar
Philipp Arras committed
31
from ..operators.adder import Adder
32
from ..operators.contraction_operator import ContractionOperator
Philipp Arras's avatar
Philipp Arras committed
33
from ..operators.diagonal_operator import DiagonalOperator
Philipp Arras's avatar
Philipp Arras committed
34
from ..operators.distributors import PowerDistributor
Philipp Arras's avatar
Philipp Arras committed
35
from ..operators.endomorphic_operator import EndomorphicOperator
Martin Reinecke's avatar
Martin Reinecke committed
36
from ..operators.harmonic_operators import HarmonicTransformOperator
Philipp Arras's avatar
Philipp Arras committed
37
from ..operators.linear_operator import LinearOperator
Philipp Arras's avatar
Philipp Arras committed
38
from ..operators.operator import Operator
Philipp Arras's avatar
Philipp Arras committed
39
from ..operators.simple_linear_operators import ducktape
40
from ..operators.normal_operators import NormalTransform, LognormalTransform
41
from ..probing import StatCalculator
Philipp Arras's avatar
Philipp Arras committed
42
from ..sugar import full, makeDomain, makeField, makeOp
43

44

Philipp Arras's avatar
Philipp Arras committed
45
def _log_k_lengths(pspace):
Philipp Arras's avatar
Philipp Arras committed
46
    """Log(k_lengths) without zeromode"""
Philipp Arras's avatar
Philipp Arras committed
47
48
49
    return np.log(pspace.k_lengths[1:])


Philipp Arras's avatar
Philipp Arras committed
50
def _relative_log_k_lengths(power_space):
Philipp Arras's avatar
Philipp Arras committed
51
52
    """Log-distance to first bin
    logkl.shape==power_space.shape, logkl[0]=logkl[1]=0"""
Philipp Arras's avatar
Philipp Arras committed
53
54
55
56
57
58
    power_space = DomainTuple.make(power_space)
    assert isinstance(power_space[0], PowerSpace)
    assert len(power_space) == 1
    logkl = _log_k_lengths(power_space[0])
    assert logkl.shape[0] == power_space[0].shape[0] - 1
    logkl -= logkl[0]
Philipp Arras's avatar
Philipp Arras committed
59
    return np.insert(logkl, 0, 0)
Philipp Arras's avatar
Philipp Arras committed
60
61


Philipp Arras's avatar
Philipp Arras committed
62
def _log_vol(power_space):
63
    power_space = makeDomain(power_space)
Philipp Arras's avatar
Philipp Arras committed
64
65
66
67
68
    assert isinstance(power_space[0], PowerSpace)
    logk_lengths = _log_k_lengths(power_space[0])
    return logk_lengths[1:] - logk_lengths[:-1]


Philipp Haim's avatar
Philipp Haim committed
69
70
71
72
73
74
def _structured_spaces(domain):
    if isinstance(domain[0], UnstructuredDomain):
        return np.arange(1, len(domain))
    return np.arange(len(domain))


Philipp Haim's avatar
Philipp Haim committed
75
def _total_fluctuation_realized(samples):
Philipp Haim's avatar
Philipp Haim committed
76
77
78
    spaces = _structured_spaces(samples[0].domain)
    co = ContractionOperator(samples[0].domain, spaces)
    size = co.domain.size/co.target.size
79
80
    res = 0.
    for s in samples:
Philipp Haim's avatar
Philipp Haim committed
81
82
        res = res + (s - co.adjoint(co(s)/size))**2
    res = res.mean(spaces)/len(samples)
Philipp Haim's avatar
Philipp Haim committed
83
    return np.sqrt(res if np.isscalar(res) else res.val)
84
85


Philipp Frank's avatar
Philipp Frank committed
86
class _SlopeRemover(EndomorphicOperator):
Martin Reinecke's avatar
Martin Reinecke committed
87
    def __init__(self, domain, space=0):
Philipp Frank's avatar
Philipp Frank committed
88
        self._domain = makeDomain(domain)
89
90
        assert isinstance(self._domain[space], PowerSpace)
        logkl = _relative_log_k_lengths(self._domain[space])
91
        self._sc = logkl/float(logkl[-1])
Philipp Arras's avatar
Philipp Arras committed
92

93
        self._space = space
Philipp Haim's avatar
Philipp Haim committed
94
95
96
        axis = self._domain.axes[space][0]
        self._last = (slice(None),)*axis + (-1,) + (None,)
        self._extender = (None,)*(axis) + (slice(None),) + (None,)*(self._domain.axes[-1][-1]-axis)
Philipp Frank's avatar
Philipp Frank committed
97
        self._capability = self.TIMES | self.ADJOINT_TIMES
Philipp Arras's avatar
Philipp Arras committed
98

99
100
    def apply(self, x, mode):
        self._check_input(x, mode)
Martin Reinecke's avatar
stage 3  
Martin Reinecke committed
101
        x = x.val
Philipp Frank's avatar
Philipp Frank committed
102
        if mode == self.TIMES:
Philipp Haim's avatar
Philipp Haim committed
103
            res = x - x[self._last]*self._sc[self._extender]
Philipp Frank's avatar
Philipp Frank committed
104
        else:
105
106
            res = x.copy()
            res[self._last] -= (x*self._sc[self._extender]).sum(
Philipp Arras's avatar
Philipp Arras committed
107
                axis=self._space, keepdims=True)
Martin Reinecke's avatar
Martin Reinecke committed
108
        return makeField(self._tgt(mode), res)
Philipp Frank's avatar
Philipp Frank committed
109

Philipp Arras's avatar
Philipp Arras committed
110
111

class _TwoLogIntegrations(LinearOperator):
Martin Reinecke's avatar
Martin Reinecke committed
112
    def __init__(self, target, space=0):
Philipp Arras's avatar
Philipp Arras committed
113
        self._target = makeDomain(target)
114
115
116
117
118
        assert isinstance(self.target[space], PowerSpace)
        dom = list(self._target)
        dom[space] = UnstructuredDomain((2, self.target[space].shape[0]-2))
        self._domain = makeDomain(dom)
        self._space = space
119
        self._log_vol = _log_vol(self._target[space])
Philipp Arras's avatar
Philipp Arras committed
120
121
122
123
        self._capability = self.TIMES | self.ADJOINT_TIMES

    def apply(self, x, mode):
        self._check_input(x, mode)
124

Martin Reinecke's avatar
Martin Reinecke committed
125
        # Maybe make class properties
126
127
        axis = self._target.axes[self._space][0]
        sl = (slice(None),)*axis
Philipp Haim's avatar
Fixes  
Philipp Haim committed
128
        extender_sl = (None,)*axis + (slice(None),) + (None,)*(self._target.axes[-1][-1] - axis)
129
130
        first = sl + (0,)
        second = sl + (1,)
Martin Reinecke's avatar
Martin Reinecke committed
131
132
133
        from_third = sl + (slice(2, None),)
        no_border = sl + (slice(1, -1),)
        reverse = sl + (slice(None, None, -1),)
134

Philipp Arras's avatar
Philipp Arras committed
135
        if mode == self.TIMES:
Martin Reinecke's avatar
stage 3  
Martin Reinecke committed
136
            x = x.val
Philipp Arras's avatar
Philipp Arras committed
137
            res = np.empty(self._target.shape)
138
            res[first] = res[second] = 0
Martin Reinecke's avatar
Martin Reinecke committed
139
            res[from_third] = np.cumsum(x[second], axis=axis)
140
            res[from_third] = (res[from_third] + res[no_border])/2*self._log_vol[extender_sl] + x[first]
Martin Reinecke's avatar
Martin Reinecke committed
141
            res[from_third] = np.cumsum(res[from_third], axis=axis)
Philipp Arras's avatar
Philipp Arras committed
142
        else:
Martin Reinecke's avatar
Martin Reinecke committed
143
            x = x.val_rw()
Philipp Arras's avatar
Philipp Arras committed
144
            res = np.zeros(self._domain.shape)
Martin Reinecke's avatar
Martin Reinecke committed
145
            x[from_third] = np.cumsum(x[from_third][reverse], axis=axis)[reverse]
146
            res[first] += x[from_third]
147
            x[from_third] *= (self._log_vol/2.)[extender_sl]
148
            x[no_border] += x[from_third]
Martin Reinecke's avatar
Martin Reinecke committed
149
            res[second] += np.cumsum(x[from_third][reverse], axis=axis)[reverse]
Martin Reinecke's avatar
Martin Reinecke committed
150
        return makeField(self._tgt(mode), res)
Philipp Arras's avatar
Philipp Arras committed
151
152
153


class _Normalization(Operator):
Martin Reinecke's avatar
Martin Reinecke committed
154
    def __init__(self, domain, space=0):
Philipp Arras's avatar
Philipp Arras committed
155
        self._domain = self._target = DomainTuple.make(domain)
156
        assert isinstance(self._domain[space], PowerSpace)
157
158
159
        hspace = list(self._domain)
        hspace[space] = hspace[space].harmonic_partner
        hspace = makeDomain(hspace)
Philipp Arras's avatar
Philipp Arras committed
160
161
162
        pd = PowerDistributor(hspace,
                              power_space=self._domain[space],
                              space=space)
Martin Reinecke's avatar
Martin Reinecke committed
163
        mode_multiplicity = pd.adjoint(full(pd.target, 1.)).val_rw()
164
        zero_mode = (slice(None),)*self._domain.axes[space][0] + (0,)
Philipp Haim's avatar
Philipp Haim committed
165
        mode_multiplicity[zero_mode] = 0
Philipp Arras's avatar
Philipp Arras committed
166
        self._mode_multiplicity = makeOp(makeField(self._domain, mode_multiplicity))
167
        self._specsum = _SpecialSum(self._domain, space)
Philipp Arras's avatar
Philipp Arras committed
168
169
170

    def apply(self, x):
        self._check_input(x)
Martin Reinecke's avatar
Martin Reinecke committed
171
        amp = x.ptw("exp")
172
        spec = amp**2
Philipp Arras's avatar
Philipp Arras committed
173
174
        # FIXME This normalizes also the zeromode which is supposed to be left
        # untouched by this operator
Philipp Arras's avatar
Philipp Arras committed
175
        return self._specsum(self._mode_multiplicity(spec))**(-0.5)*amp
Philipp Arras's avatar
Philipp Arras committed
176
177
178


class _SpecialSum(EndomorphicOperator):
Martin Reinecke's avatar
Martin Reinecke committed
179
    def __init__(self, domain, space=0):
Philipp Arras's avatar
Philipp Arras committed
180
181
        self._domain = makeDomain(domain)
        self._capability = self.TIMES | self.ADJOINT_TIMES
182
        self._contractor = ContractionOperator(domain, space)
Philipp Arras's avatar
Philipp Arras committed
183
184
185

    def apply(self, x, mode):
        self._check_input(x, mode)
186
        return self._contractor.adjoint(self._contractor(x))
Philipp Arras's avatar
Philipp Arras committed
187
188


Philipp Haim's avatar
Philipp Haim committed
189
class _Distributor(LinearOperator):
Lukas Platz's avatar
Lukas Platz committed
190
    def __init__(self, dofdex, domain, target):
191
192
193
        self._dofdex = np.array(dofdex)
        self._target = DomainTuple.make(target)
        self._domain = DomainTuple.make(domain)
Philipp Haim's avatar
Philipp Haim committed
194
195
196
197
        self._capability = self.TIMES | self.ADJOINT_TIMES

    def apply(self, x, mode):
        self._check_input(x, mode)
Martin Reinecke's avatar
stage 3  
Martin Reinecke committed
198
        x = x.val
Philipp Haim's avatar
Philipp Haim committed
199
200
201
        if mode == self.TIMES:
            res = x[self._dofdex]
        else:
202
            res = np.zeros(self._tgt(mode).shape, dtype=x.dtype)
203
            res = utilities.special_add_at(res, 0, self._dofdex, x)
Martin Reinecke's avatar
Martin Reinecke committed
204
        return makeField(self._tgt(mode), res)
Martin Reinecke's avatar
Martin Reinecke committed
205

206

207
208
class _Amplitude(Operator):
    def __init__(self, target, fluctuations, flexibility, asperity,
209
                 loglogavgslope, azm, totvol, key, dofdex):
Philipp Arras's avatar
Philipp Arras committed
210
211
212
213
214
215
216
217
218
219
220
        """
        fluctuations > 0
        flexibility > 0
        asperity > 0
        loglogavgslope probably negative
        """
        assert isinstance(fluctuations, Operator)
        assert isinstance(flexibility, Operator)
        assert isinstance(asperity, Operator)
        assert isinstance(loglogavgslope, Operator)

Philipp Haim's avatar
Philipp Haim committed
221
222
        if len(dofdex) > 0:
            N_copies = max(dofdex) + 1
Philipp Haim's avatar
Philipp Haim committed
223
            space = 1
Philipp Frank's avatar
cleanup  
Philipp Frank committed
224
225
            distributed_tgt = makeDomain((UnstructuredDomain(len(dofdex)),
                                          target))
Philipp Haim's avatar
Philipp Haim committed
226
            target = makeDomain((UnstructuredDomain(N_copies), target))
Lukas Platz's avatar
Lukas Platz committed
227
            Distributor = _Distributor(dofdex, target, distributed_tgt)
Philipp Haim's avatar
Philipp Haim committed
228
        else:
Philipp Haim's avatar
Philipp Haim committed
229
            N_copies = 0
Philipp Haim's avatar
Philipp Haim committed
230
            space = 0
Philipp Haim's avatar
Philipp Haim committed
231
            distributed_tgt = target = makeDomain(target)
Martin Reinecke's avatar
Martin Reinecke committed
232
        azm_expander = ContractionOperator(distributed_tgt, spaces=space).adjoint
Philipp Haim's avatar
Philipp Haim committed
233
        assert isinstance(target[space], PowerSpace)
Martin Reinecke's avatar
Martin Reinecke committed
234

235
        twolog = _TwoLogIntegrations(target, space)
Philipp Arras's avatar
Philipp Arras committed
236
        dom = twolog.domain
237

238
        shp = dom[space].shape
Martin Reinecke's avatar
Martin Reinecke committed
239
240
        expander = ContractionOperator(dom, spaces=space).adjoint
        ps_expander = ContractionOperator(twolog.target, spaces=space).adjoint
Philipp Arras's avatar
Philipp Arras committed
241
242
243

        # Prepare constant fields
        foo = np.zeros(shp)
244
        foo[0] = foo[1] = np.sqrt(_log_vol(target[space]))
Martin Reinecke's avatar
Martin Reinecke committed
245
        vflex = DiagonalOperator(makeField(dom[space], foo), dom, space)
Philipp Arras's avatar
Philipp Arras committed
246
247
248

        foo = np.zeros(shp, dtype=np.float64)
        foo[0] += 1
Martin Reinecke's avatar
Martin Reinecke committed
249
        vasp = DiagonalOperator(makeField(dom[space], foo), dom, space)
Philipp Arras's avatar
Philipp Arras committed
250
251

        foo = np.ones(shp)
252
        foo[0] = _log_vol(target[space])**2/12.
Martin Reinecke's avatar
Martin Reinecke committed
253
        shift = DiagonalOperator(makeField(dom[space], foo), dom, space)
Martin Reinecke's avatar
Martin Reinecke committed
254

255
        vslope = DiagonalOperator(
Philipp Arras's avatar
Philipp Arras committed
256
            makeField(target[space], _relative_log_k_lengths(target[space])),
Martin Reinecke's avatar
Martin Reinecke committed
257
            target, space)
258
259

        foo, bar = [np.zeros(target[space].shape) for _ in range(2)]
Philipp Arras's avatar
Philipp Arras committed
260
        bar[1:] = foo[0] = totvol
Philipp Arras's avatar
Philipp Arras committed
261
262
263
264
        vol0, vol1 = [
            DiagonalOperator(makeField(target[space], aa), target, space)
            for aa in (foo, bar)
        ]
265

Martin Reinecke's avatar
Martin Reinecke committed
266
        # Prepare fields for Adder
267
        shift, vol0 = [op(full(op.domain, 1)) for op in (shift, vol0)]
Philipp Arras's avatar
Philipp Arras committed
268
269
        # End prepare constant fields

270
271
272
273
        slope = vslope @ ps_expander @ loglogavgslope
        sig_flex = vflex @ expander @ flexibility
        sig_asp = vasp @ expander @ asperity
        sig_fluc = vol1 @ ps_expander @ fluctuations
Philipp Haim's avatar
Philipp Haim committed
274
        sig_fluc = vol1 @ ps_expander @ fluctuations
Philipp Arras's avatar
Philipp Arras committed
275
276

        xi = ducktape(dom, None, key)
Martin Reinecke's avatar
Martin Reinecke committed
277
        sigma = sig_flex*(Adder(shift) @ sig_asp).ptw("sqrt")
278
279
        smooth = _SlopeRemover(target, space) @ twolog @ (sigma*xi)
        op = _Normalization(target, space) @ (slope + smooth)
Philipp Haim's avatar
Philipp Haim committed
280
        if N_copies > 0:
Philipp Haim's avatar
Philipp Haim committed
281
282
            op = Distributor @ op
            sig_fluc = Distributor @ sig_fluc
Martin Reinecke's avatar
Martin Reinecke committed
283
            op = Adder(Distributor(vol0)) @ (sig_fluc*(azm_expander @ azm.ptw("reciprocal"))*op)
Philipp Arras's avatar
Philipp Arras committed
284
285
            self._fluc = (_Distributor(dofdex, fluctuations.target,
                                       distributed_tgt[0]) @ fluctuations)
Philipp Haim's avatar
Philipp Haim committed
286
        else:
Martin Reinecke's avatar
Martin Reinecke committed
287
            op = Adder(vol0) @ (sig_fluc*(azm_expander @ azm.ptw("reciprocal"))*op)
Philipp Frank's avatar
fixup  
Philipp Frank committed
288
            self._fluc = fluctuations
Philipp Arras's avatar
Philipp Arras committed
289

Philipp Arras's avatar
Philipp Arras committed
290
291
        self.apply = op.apply
        self._domain, self._target = op.domain, op.target
292
        self._space = space
293
        self._repr_str = "_Amplitude: " + op.__repr__()
Philipp Arras's avatar
Philipp Arras committed
294

Philipp Arras's avatar
Philipp Arras committed
295
296
297
298
    @property
    def fluctuation_amplitude(self):
        return self._fluc

299
300
301
    def __repr__(self):
        return self._repr_str

302
303

class CorrelatedFieldMaker:
Lukas Platz's avatar
Lukas Platz committed
304
305
306
    """Constrution helper for hirarchical correlated field models.

    The correlated field models are parametrized by creating
307
308
    power spectrum operators ("amplitudes") via calls to
    :func:`add_fluctuations` that act on the targeted field subdomains.
Lukas Platz's avatar
Lukas Platz committed
309
    During creation of the :class:`CorrelatedFieldMaker` via
310
311
312
    :func:`make`, a global offset from zero of the field model
    can be defined and an operator applying fluctuations
    around this offset is parametrized.
Lukas Platz's avatar
Lukas Platz committed
313
314

    The resulting correlated field model operator has a
Martin Reinecke's avatar
Martin Reinecke committed
315
    :class:`~nifty7.multi_domain.MultiDomain` as its domain and
Lukas Platz's avatar
Lukas Platz committed
316
317
318
    expects its input values to be univariately gaussian.

    The target of the constructed operator will be a
Martin Reinecke's avatar
merge  
Martin Reinecke committed
319
    :class:`~nifty7.domain_tuple.DomainTuple` containing the
320
321
    `target_subdomains` of the added fluctuations in the order of
    the `add_fluctuations` calls.
Lukas Platz's avatar
Lukas Platz committed
322

323
    Creation of the model operator is completed by calling the method
Lukas Platz's avatar
Lukas Platz committed
324
325
    :func:`finalize`, which returns the configured operator.

326
327
328
329
330
331
332
333
334
335
    An operator representing an array of correlated field models
    can be constructed by setting the `total_N` parameter of
    :func:`make`. It will have an
    :class:`~nifty.domains.unstructucture_domain.UnstructureDomain`
    of shape `(total_N,)` prepended to its target domain and represent
    `total_N` correlated fields simulataneously.
    The degree of information sharing between the correlated field
    models can be configured via the `dofdex` parameters
    of :func:`make` and :func:`add_fluctuations`.

Lukas Platz's avatar
Lukas Platz committed
336
    See the methods :func:`make`, :func:`add_fluctuations`
337
    and :func:`finalize` for further usage information."""
338
339
340
    def __init__(self, offset_mean, offset_fluctuations_op, prefix, total_N):
        if not isinstance(offset_fluctuations_op, Operator):
            raise TypeError("offset_fluctuations_op needs to be an operator")
341
        self._a = []
342
        self._target_subdomains = []
Philipp Arras's avatar
Formats  
Philipp Arras committed
343

344
345
        self._offset_mean = offset_mean
        self._azm = offset_fluctuations_op
346
        self._prefix = prefix
Philipp Haim's avatar
Philipp Haim committed
347
        self._total_N = total_N
Philipp Arras's avatar
Formats  
Philipp Arras committed
348

349
    @staticmethod
Philipp Arras's avatar
Philipp Arras committed
350
    def make(offset_mean, offset_std_mean, offset_std_std, prefix, total_N=0,
Martin Reinecke's avatar
Martin Reinecke committed
351
             dofdex=None):
Lukas Platz's avatar
Lukas Platz committed
352
353
354
355
356
357
358
        """Returns a CorrelatedFieldMaker object.

        Parameters
        ----------
        offset_mean : float
            Mean offset from zero of the correlated field to be made.
        offset_std_mean : float
Lukas Platz's avatar
Lukas Platz committed
359
            Mean standard deviation of the offset.
Lukas Platz's avatar
Lukas Platz committed
360
        offset_std_std : float
Lukas Platz's avatar
Lukas Platz committed
361
            Standard deviation of the stddev of the offset.
Lukas Platz's avatar
Lukas Platz committed
362
363
        prefix : string
            Prefix to the names of the domains of the cf operator to be made.
Lukas Platz's avatar
Lukas Platz committed
364
            This determines the names of the operator domain.
365
366
        total_N : integer, optional
            Number of field models to create.
Lukas Platz's avatar
Lukas Platz committed
367
368
369
            If not 0, the first entry of the operators target will be an
            :class:`~nifty.domains.unstructured_domain.UnstructuredDomain`
            with length `total_N`.
370
        dofdex : np.array of integers, optional
Philipp Arras's avatar
Philipp Arras committed
371
372
373
            An integer array specifying the zero mode models used if
            total_N > 1. It needs to have length of total_N. If total_N=3 and
            dofdex=[0,0,1], that means that two models for the zero mode are
374
375
376
377
            instanciated; the first one is used for the first and second
            field model and the second is used for the third field model.
            *If not specified*, use the same zero mode model for all
            constructed field models.
Lukas Platz's avatar
Lukas Platz committed
378
        """
Philipp Frank's avatar
Philipp Frank committed
379
380
        if dofdex is None:
            dofdex = np.full(total_N, 0)
381
382
        elif len(dofdex) != total_N:
            raise ValueError("length of dofdex needs to match total_N")
Philipp Frank's avatar
Philipp Frank committed
383
        N = max(dofdex) + 1 if total_N > 0 else 0
384
385
        zm = LognormalTransform(offset_std_mean, offset_std_std,
                                prefix + 'zeromode', N)
Philipp Frank's avatar
fixup  
Philipp Frank committed
386
        if total_N > 0:
Martin Reinecke's avatar
Martin Reinecke committed
387
            zm = _Distributor(dofdex, zm.target, UnstructuredDomain(total_N)) @ zm
388
        return CorrelatedFieldMaker(offset_mean, zm, prefix, total_N)
389
390

    def add_fluctuations(self,
391
                         target_subdomain,
392
393
394
395
396
397
398
399
                         fluctuations_mean,
                         fluctuations_stddev,
                         flexibility_mean,
                         flexibility_stddev,
                         asperity_mean,
                         asperity_stddev,
                         loglogavgslope_mean,
                         loglogavgslope_stddev,
Martin Reinecke's avatar
Martin Reinecke committed
400
401
402
403
                         prefix='',
                         index=None,
                         dofdex=None,
                         harmonic_partner=None):
Lukas Platz's avatar
Lukas Platz committed
404
405
406
407
408
409
        """Function to add correlation structures to the field to be made.

        Correlations are described by their power spectrum and the subdomain
        on which they apply.

        The parameters `fluctuations`, `flexibility`, `asperity` and
410
411
        `loglogavgslope` configure the power spectrum model ("amplitude")
        used on the target field subdomain `target_subdomain`.
Lukas Platz's avatar
Lukas Platz committed
412
413
414
415
416
417
418
419
420
421
422
        It is assembled as the sum of a power law component
        (linear slope in log-log power-frequency-space),
        a smooth varying component (integrated wiener process) and
        a ragged componenent (unintegrated wiener process).

        Multiple calls to `add_fluctuations` are possible, in which case
        the constructed field will have the outer product of the individual
        power spectra as its global power spectrum.

        Parameters
        ----------
Martin Reinecke's avatar
Martin Reinecke committed
423
424
        target_subdomain : :class:`~nifty7.domain.Domain`, \
                           :class:`~nifty7.domain_tuple.DomainTuple`
Lukas Platz's avatar
Lukas Platz committed
425
426
427
428
429
            Target subdomain on which the correlation structure defined
            in this call should hold.
        fluctuations_{mean,stddev} : float
            Total spectral energy -> Amplitude of the fluctuations
        flexibility_{mean,stddev} : float
430
            Amplitude of the non-power-law power spectrum component
Lukas Platz's avatar
Lukas Platz committed
431
        asperity_{mean,stddev} : float
432
            Roughness of the non-power-law power spectrum component
Lukas Platz's avatar
Lukas Platz committed
433
434
435
436
437
            Used to accomodate single frequency peaks
        loglogavgslope_{mean,stddev} : float
            Power law component exponent
        prefix : string
            prefix of the power spectrum parameter domain names
Philipp Arras's avatar
Philipp Arras committed
438
439
440
        index : int
            Position target_subdomain in the final total domain of the
            correlated field operator.
441
442
        dofdex : np.array, optional
            An integer array specifying the power spectrum models used if
Philipp Arras's avatar
Philipp Arras committed
443
            total_N > 1. It needs to have length of total_N. If total_N=3 and
444
445
446
447
448
            dofdex=[0,0,1], that means that two power spectrum models are
            instanciated; the first one is used for the first and second
            field model and the second one is used for the third field model.
            *If not given*, use the same power spectrum model for all
            constructed field models.
Martin Reinecke's avatar
Martin Reinecke committed
449
450
        harmonic_partner : :class:`~nifty7.domain.Domain`, \
                           :class:`~nifty7.domain_tuple.DomainTuple`
Lukas Platz's avatar
Lukas Platz committed
451
452
            In which harmonic space to define the power spectrum
        """
Philipp Frank's avatar
Philipp Frank committed
453
        if harmonic_partner is None:
454
            harmonic_partner = target_subdomain.get_default_codomain()
Philipp Frank's avatar
Fixup  
Philipp Frank committed
455
        else:
456
457
            target_subdomain.check_codomain(harmonic_partner)
            harmonic_partner.check_codomain(target_subdomain)
458

Philipp Haim's avatar
Philipp Haim committed
459
460
        if dofdex is None:
            dofdex = np.full(self._total_N, 0)
461
462
        elif len(dofdex) != self._total_N:
            raise ValueError("length of dofdex needs to match total_N")
Philipp Haim's avatar
Philipp Haim committed
463

Philipp Haim's avatar
Philipp Haim committed
464
465
        if self._total_N > 0:
            N = max(dofdex) + 1
466
            target_subdomain = makeDomain((UnstructuredDomain(N), target_subdomain))
Philipp Haim's avatar
Philipp Haim committed
467
        else:
Philipp Haim's avatar
Philipp Haim committed
468
            N = 0
469
            target_subdomain = makeDomain(target_subdomain)
Philipp Arras's avatar
Philipp Arras committed
470
        prefix = str(prefix)
471
        # assert isinstance(target_subdomain[space], (RGSpace, HPSpace, GLSpace)
Philipp Arras's avatar
Philipp Arras committed
472

473
474
475
476
477
478
479
480
481
        fluct = LognormalTransform(fluctuations_mean, fluctuations_stddev,
                                   self._prefix + prefix + 'fluctuations', N)
        flex = LognormalTransform(flexibility_mean, flexibility_stddev,
                                  self._prefix + prefix + 'flexibility', N)
        asp = LognormalTransform(asperity_mean, asperity_stddev,
                                 self._prefix + prefix + 'asperity', N)
        avgsl = NormalTransform(loglogavgslope_mean, loglogavgslope_stddev,
                                self._prefix + prefix + 'loglogavgslope', N)

Philipp Arras's avatar
Philipp Arras committed
482
        amp = _Amplitude(PowerSpace(harmonic_partner), fluct, flex, asp, avgsl,
483
                         self._azm, target_subdomain[-1].total_volume,
484
                         self._prefix + prefix + 'spectrum', dofdex)
Philipp Haim's avatar
Philipp Haim committed
485

486
487
        if index is not None:
            self._a.insert(index, amp)
488
            self._target_subdomains.insert(index, target_subdomain)
489
490
        else:
            self._a.append(amp)
491
            self._target_subdomains.append(target_subdomain)
492

Philipp Arras's avatar
Philipp Arras committed
493
494
495
496
497
498
499
500
501
502
    def finalize(self, prior_info=100):
        """Finishes model construction process and returns the constructed
        operator.

        Parameters
        ----------
        prior_info : integer
            How many prior samples to draw for property verification statistics
            If zero, skips calculating and displaying statistics.
        """
Philipp Haim's avatar
Philipp Haim committed
503
        n_amplitudes = len(self._a)
Philipp Haim's avatar
Philipp Haim committed
504
        if self._total_N > 0:
Philipp Arras's avatar
Philipp Arras committed
505
506
507
            hspace = makeDomain(
                [UnstructuredDomain(self._total_N)] +
                [dd.target[-1].harmonic_partner for dd in self._a])
Philipp Haim's avatar
Philipp Haim committed
508
509
            spaces = tuple(range(1, n_amplitudes + 1))
            amp_space = 1
Philipp Haim's avatar
Philipp Haim committed
510
511
        else:
            hspace = makeDomain(
Philipp Arras's avatar
Philipp Arras committed
512
                [dd.target[0].harmonic_partner for dd in self._a])
Philipp Haim's avatar
Philipp Haim committed
513
            spaces = tuple(range(n_amplitudes))
Philipp Haim's avatar
Philipp Haim committed
514
            amp_space = 0
515

Martin Reinecke's avatar
Martin Reinecke committed
516
        expander = ContractionOperator(hspace, spaces=spaces).adjoint
Philipp Frank's avatar
fixup  
Philipp Frank committed
517
        azm = expander @ self._azm
518

519
        ht = HarmonicTransformOperator(hspace,
520
                                       self._target_subdomains[0][amp_space],
Martin Reinecke's avatar
Martin Reinecke committed
521
                                       space=spaces[0])
522
        for i in range(1, n_amplitudes):
523
            ht = HarmonicTransformOperator(ht.target,
524
                                           self._target_subdomains[i][amp_space],
525
526
527
528
529
                                           space=spaces[i]) @ ht
        a = []
        for ii in range(n_amplitudes):
            co = ContractionOperator(hspace, spaces[:ii] + spaces[ii + 1:])
            pp = self._a[ii].target[amp_space]
Philipp Haim's avatar
Philipp Haim committed
530
            pd = PowerDistributor(co.target, pp, amp_space)
531
532
            a.append(co.adjoint @ pd @ self._a[ii])
        corr = reduce(mul, a)
Philipp Arras's avatar
Philipp Arras committed
533
        op = ht(azm*corr*ducktape(hspace, None, self._prefix + 'xi'))
Philipp Arras's avatar
Philipp Arras committed
534

535
536
        if self._offset_mean is not None:
            offset = self._offset_mean
537
538
539
540
541
542
543
            # Deviations from this offset must not be considered here as they
            # are learned by the zeromode
            if isinstance(offset, (Field, MultiField)):
                op = Adder(offset) @ op
            else:
                offset = float(offset)
                op = Adder(full(op.target, offset)) @ op
544
        self.statistics_summary(prior_info)
545
546
        return op

547
548
549
550
551
552
    def statistics_summary(self, prior_info):
        from ..sugar import from_random

        if prior_info == 0:
            return

553
554
        lst = [('Offset amplitude', self.amplitude_total_offset),
               ('Total fluctuation amplitude', self.total_fluctuation)]
555
        namps = len(self._a)
556
557
558
559
560
561
562
563
        if namps > 1:
            for ii in range(namps):
                lst.append(('Slice fluctuation (space {})'.format(ii),
                            self.slice_fluctuation(ii)))
                lst.append(('Average fluctuation (space {})'.format(ii),
                            self.average_fluctuation(ii)))

        for kk, op in lst:
564
565
            sc = StatCalculator()
            for _ in range(prior_info):
566
                sc.add(op(from_random(op.domain, 'normal')))
Martin Reinecke's avatar
merge  
Martin Reinecke committed
567
            mean = sc.mean.val
Martin Reinecke's avatar
Martin Reinecke committed
568
            stddev = sc.var.ptw("sqrt").val
569
            for m, s in zip(mean.flatten(), stddev.flatten()):
570
                logger.info('{}: {:.02E} ± {:.02E}'.format(kk, m, s))
571
572
573

    def moment_slice_to_average(self, fluctuations_slice_mean, nsamples=1000):
        fluctuations_slice_mean = float(fluctuations_slice_mean)
574
575
576
        if not fluctuations_slice_mean > 0:
            msg = "fluctuations_slice_mean must be greater zero; got {!r}"
            raise ValueError(msg.format(fluctuations_slice_mean))
577
        from ..sugar import from_random
578
579
        scm = 1.
        for a in self._a:
Martin Reinecke's avatar
Martin Reinecke committed
580
            op = a.fluctuation_amplitude*self._azm.ptw("reciprocal")
581
            res = np.array([op(from_random(op.domain, 'normal')).val
582
583
                            for _ in range(nsamples)])
            scm *= res**2 + 1.
584
        return fluctuations_slice_mean/np.mean(np.sqrt(scm))
585

Philipp Arras's avatar
Philipp Arras committed
586
    @property
Philipp Haim's avatar
Philipp Haim committed
587
    def normalized_amplitudes(self):
588
        """Returns the power spectrum operators used in the model"""
589
        return self._a
Philipp Arras's avatar
Philipp Arras committed
590

Philipp Haim's avatar
Philipp Haim committed
591
592
593
594
595
596
597
    @property
    def amplitude(self):
        if len(self._a) > 1:
            s = ('If more than one spectrum is present in the model,',
                 ' no unique set of amplitudes exist because only the',
                 ' relative scale is determined.')
            raise NotImplementedError(s)
Philipp Haim's avatar
Fix  
Philipp Haim committed
598
599
        dom = self._a[0].target
        expand = ContractionOperator(dom, len(dom)-1).adjoint
Philipp Haim's avatar
Philipp Haim committed
600
601
        return self._a[0]*(expand @ self.amplitude_total_offset)

602
603
604
    @property
    def amplitude_total_offset(self):
        return self._azm
Philipp Arras's avatar
Philipp Arras committed
605
606

    @property
607
    def total_fluctuation(self):
608
        """Returns operator which acts on prior or posterior samples"""
609
        if len(self._a) == 0:
Philipp Arras's avatar
Philipp Arras committed
610
            raise NotImplementedError
611
        if len(self._a) == 1:
612
            return self.average_fluctuation(0)
613
614
        q = 1.
        for a in self._a:
Martin Reinecke's avatar
Martin Reinecke committed
615
            fl = a.fluctuation_amplitude*self._azm.ptw("reciprocal")
Philipp Arras's avatar
Philipp Arras committed
616
            q = q*(Adder(full(fl.target, 1.)) @ fl**2)
Martin Reinecke's avatar
Martin Reinecke committed
617
        return (Adder(full(q.target, -1.)) @ q).ptw("sqrt")*self._azm
618

Philipp Arras's avatar
Philipp Arras committed
619
    def slice_fluctuation(self, space):
620
        """Returns operator which acts on prior or posterior samples"""
621
        if len(self._a) == 0:
Philipp Arras's avatar
Philipp Arras committed
622
            raise NotImplementedError
623
        if space >= len(self._a):
624
            raise ValueError("invalid space specified; got {!r}".format(space))
625
        if len(self._a) == 1:
626
            return self.average_fluctuation(0)
627
628
        q = 1.
        for j in range(len(self._a)):
Martin Reinecke's avatar
Martin Reinecke committed
629
            fl = self._a[j].fluctuation_amplitude*self._azm.ptw("reciprocal")
630
            if j == space:
Philipp Arras's avatar
Philipp Arras committed
631
                q = q*fl**2
632
            else:
Philipp Arras's avatar
Philipp Arras committed
633
                q = q*(Adder(full(fl.target, 1.)) @ fl**2)
Martin Reinecke's avatar
Martin Reinecke committed
634
        return q.ptw("sqrt")*self._azm
Philipp Arras's avatar
Philipp Arras committed
635
636

    def average_fluctuation(self, space):
637
        """Returns operator which acts on prior or posterior samples"""
638
        if len(self._a) == 0:
Philipp Arras's avatar
Philipp Arras committed
639
            raise NotImplementedError
640
        if space >= len(self._a):
641
            raise ValueError("invalid space specified; got {!r}".format(space))
642
        if len(self._a) == 1:
Philipp Haim's avatar
Philipp Haim committed
643
644
            return self._a[0].fluctuation_amplitude
        return self._a[space].fluctuation_amplitude
645

646
647
    @staticmethod
    def offset_amplitude_realized(samples):
Philipp Haim's avatar
Philipp Haim committed
648
        spaces = _structured_spaces(samples[0].domain)
649
650
        res = 0.
        for s in samples:
Philipp Haim's avatar
Philipp Haim committed
651
            res = res + s.mean(spaces)**2
Philipp Haim's avatar
Philipp Haim committed
652
653
        res = res/len(samples)
        return np.sqrt(res if np.isscalar(res) else res.val)
Philipp Arras's avatar
Philipp Arras committed
654

655
656
657
658
659
660
661
662
    @staticmethod
    def total_fluctuation_realized(samples):
        return _total_fluctuation_realized(samples)

    @staticmethod
    def slice_fluctuation_realized(samples, space):
        """Computes slice fluctuations from collection of field (defined in signal
        space) realizations."""
Philipp Haim's avatar
Philipp Haim committed
663
664
        spaces = _structured_spaces(samples[0].domain)
        if space >= len(spaces):
665
            raise ValueError("invalid space specified; got {!r}".format(space))
Philipp Haim's avatar
Philipp Haim committed
666
        if len(spaces) == 1:
667
            return _total_fluctuation_realized(samples)
Philipp Haim's avatar
Philipp Haim committed
668
        space = space + spaces[0]
Philipp Arras's avatar
Philipp Arras committed
669
        res1, res2 = 0., 0.
670
        for s in samples:
Philipp Frank's avatar
fixes  
Philipp Frank committed
671
672
673
674
            res1 = res1 + s**2
            res2 = res2 + s.mean(space)**2
        res1 = res1/len(samples)
        res2 = res2/len(samples)
Philipp Haim's avatar
Philipp Haim committed
675
        res = res1.mean(spaces) - res2.mean(spaces[:-1])
Philipp Haim's avatar
Philipp Haim committed
676
        return np.sqrt(res if np.isscalar(res) else res.val)
Philipp Frank's avatar
fixes  
Philipp Frank committed
677

Philipp Arras's avatar
Philipp Arras committed
678
    @staticmethod
679
680
681
    def average_fluctuation_realized(samples, space):
        """Computes average fluctuations from collection of field (defined in signal
        space) realizations."""
Philipp Haim's avatar
Philipp Haim committed
682
683
        spaces = _structured_spaces(samples[0].domain)
        if space >= len(spaces):
684
            raise ValueError("invalid space specified; got {!r}".format(space))
Philipp Haim's avatar
Philipp Haim committed
685
        if len(spaces) == 1:
686
            return _total_fluctuation_realized(samples)
Philipp Haim's avatar
Philipp Haim committed
687
688
689
        space = space + spaces[0]
        sub_spaces = set(spaces)
        sub_spaces.remove(space)
Philipp Arras's avatar
Philipp Arras committed
690
        # Domain containing domain[space] and domain[0] iff total_N>0
Philipp Haim's avatar
Philipp Haim committed
691
        sub_dom = makeDomain([samples[0].domain[ind]
Philipp Arras's avatar
Philipp Arras committed
692
                              for ind in (set([0])-set(spaces)) | set([space])])
Philipp Haim's avatar
Philipp Haim committed
693
        co = ContractionOperator(sub_dom, len(sub_dom)-1)
694
        size = co.domain.size/co.target.size
Philipp Arras's avatar
Philipp Arras committed
695
696
        res = 0.
        for s in samples:
Philipp Haim's avatar
Philipp Haim committed
697
            r = s.mean(sub_spaces)
698
            res = res + (r - co.adjoint(co(r)/size))**2
Philipp Haim's avatar
Philipp Haim committed
699
        res = res.mean(spaces[0])/len(samples)
Philipp Haim's avatar
Philipp Haim committed
700
        return np.sqrt(res if np.isscalar(res) else res.val)