correlated_fields.py 21.8 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
Philipp Arras's avatar
Philipp Arras committed
15
# Authors: Philipp Frank, Philipp Arras
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

Philipp Arras's avatar
Philipp Arras committed
19
import numpy as np
20

Philipp Arras's avatar
Philipp Arras committed
21
from ..domain_tuple import DomainTuple
Philipp Arras's avatar
Philipp Arras committed
22
23
from ..domains.power_space import PowerSpace
from ..domains.unstructured_domain import UnstructuredDomain
24
from ..field import Field
Philipp Arras's avatar
Philipp Arras committed
25
from ..operators.adder import Adder
26
from ..operators.contraction_operator import ContractionOperator
Philipp Arras's avatar
Philipp Arras committed
27
from ..operators.distributors import PowerDistributor
Philipp Arras's avatar
Philipp Arras committed
28
from ..operators.endomorphic_operator import EndomorphicOperator
Martin Reinecke's avatar
Martin Reinecke committed
29
from ..operators.harmonic_operators import HarmonicTransformOperator
Philipp Arras's avatar
Philipp Arras committed
30
from ..operators.linear_operator import LinearOperator
31
from ..operators.diagonal_operator import DiagonalOperator
Philipp Arras's avatar
Philipp Arras committed
32
33
from ..operators.operator import Operator
from ..operators.simple_linear_operators import VdotOperator, ducktape
Philipp Arras's avatar
Philipp Arras committed
34
from ..operators.value_inserter import ValueInserter
35
from ..probing import StatCalculator
36
from ..sugar import from_global_data, full, makeDomain, get_default_codomain
37

38

39
def _reshaper(x, shape):
40
41
42
43
44
45
46
47
    x = np.array(x)
    if x.shape == shape:
        return np.asfarray(x)
    elif x.shape in [(), (1,)]:
        return np.full(shape, x, dtype=np.float)
    else:
        raise TypeError("Shape of parameters cannot be interpreted")

48
49
50
51

def _lognormal_moments(mean, sig, shape = ()):
    mean, sig = (_reshaper(param, shape) for param in (mean, sig))
    assert np.all(mean > 0 )
52
    assert np.all(sig > 0)
Philipp Arras's avatar
Philipp Arras committed
53
54
    logsig = np.sqrt(np.log((sig/mean)**2 + 1))
    logmean = np.log(mean) - logsig**2/2
55
    return logmean, logsig
Philipp Arras's avatar
Philipp Arras committed
56
57


58
def _normal(mean, sig, key, domain = DomainTuple.scalar_domain()):
59
    domain = makeDomain(domain)
60
    mean, sig = (_reshaper(param, domain.shape) for param in (mean, sig))
61
    return Adder(from_global_data(domain, mean)) @ (
62
63
        DiagonalOperator(from_global_data(domain,sig))
        @ ducktape(domain, None, key))
Philipp Arras's avatar
Philipp Arras committed
64
65


Philipp Arras's avatar
Philipp Arras committed
66
def _log_k_lengths(pspace):
Philipp Arras's avatar
Philipp Arras committed
67
    """Log(k_lengths) without zeromode"""
Philipp Arras's avatar
Philipp Arras committed
68
69
70
    return np.log(pspace.k_lengths[1:])


Philipp Arras's avatar
Philipp Arras committed
71
def _relative_log_k_lengths(power_space):
Philipp Arras's avatar
Philipp Arras committed
72
73
    """Log-distance to first bin
    logkl.shape==power_space.shape, logkl[0]=logkl[1]=0"""
Philipp Arras's avatar
Philipp Arras committed
74
75
76
77
78
79
    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
80
    return np.insert(logkl, 0, 0)
Philipp Arras's avatar
Philipp Arras committed
81
82


Philipp Arras's avatar
Philipp Arras committed
83
def _log_vol(power_space):
84
    power_space = makeDomain(power_space)
Philipp Arras's avatar
Philipp Arras committed
85
86
87
88
89
    assert isinstance(power_space[0], PowerSpace)
    logk_lengths = _log_k_lengths(power_space[0])
    return logk_lengths[1:] - logk_lengths[:-1]


90
def _total_fluctuation_realized(samples, space = 0):
91
92
    res = 0.
    for s in samples:
93
94
        res = res + (s - s.mean(space, keepdims = True))**2
    return np.sqrt((res/len(samples)).mean(space))
95
96
97
98
99
100
101
102
103


def _stats(op, samples):
    sc = StatCalculator()
    for s in samples:
        sc.add(op(s.extract(op.domain)))
    return sc.mean.to_global_data(), sc.var.sqrt().to_global_data()


Philipp Arras's avatar
Philipp Arras committed
104
class _LognormalMomentMatching(Operator):
105
106
    def __init__(self, mean, sig, key,
            domain = DomainTuple.scalar_domain()):
Philipp Arras's avatar
Philipp Arras committed
107
        key = str(key)
108
        logmean, logsig = _lognormal_moments(mean, sig, domain.shape)
Philipp Arras's avatar
Philipp Arras committed
109
110
        self._mean = mean
        self._sig = sig
111
        op = _normal(logmean, logsig, key, domain).exp()
Philipp Arras's avatar
Philipp Arras committed
112
113
114
115
116
117
118
119
120
121
        self._domain, self._target = op.domain, op.target
        self.apply = op.apply

    @property
    def mean(self):
        return self._mean

    @property
    def std(self):
        return self._sig
Philipp Arras's avatar
Philipp Arras committed
122
123


Philipp Frank's avatar
Philipp Frank committed
124
class _SlopeRemover(EndomorphicOperator):
125
    def __init__(self, domain, space = 0):
Philipp Frank's avatar
Philipp Frank committed
126
        self._domain = makeDomain(domain)
127
128
        assert isinstance(self._domain[space], PowerSpace)
        logkl = _relative_log_k_lengths(self._domain[space])
129
        self._sc = logkl/float(logkl[-1])
Philipp Arras's avatar
Philipp Arras committed
130

131
        self._space = space
Philipp Haim's avatar
Philipp Haim committed
132
133
134
        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
135
        self._capability = self.TIMES | self.ADJOINT_TIMES
Philipp Arras's avatar
Philipp Arras committed
136

137
138
    def apply(self, x, mode):
        self._check_input(x, mode)
Philipp Frank's avatar
Philipp Frank committed
139
140
        x = x.to_global_data()
        if mode == self.TIMES:
Philipp Haim's avatar
Philipp Haim committed
141
            res = x - x[self._last]*self._sc[self._extender]
Philipp Frank's avatar
Philipp Frank committed
142
        else:
143
144
145
            res = x.copy()
            res[self._last] -= (x*self._sc[self._extender]).sum(
                    axis = self._space, keepdims = True)
146
        return from_global_data(self._tgt(mode), res)
Philipp Frank's avatar
Philipp Frank committed
147

Philipp Arras's avatar
Philipp Arras committed
148
149

class _TwoLogIntegrations(LinearOperator):
150
    def __init__(self, target, space = 0):
Philipp Arras's avatar
Philipp Arras committed
151
        self._target = makeDomain(target)
152
153
154
155
156
157
        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
        logk_lengths = _log_k_lengths(self._target[space])
158
        self._log_vol = _log_vol(self._target[space])
Philipp Arras's avatar
Philipp Arras committed
159
160
161
162
        self._capability = self.TIMES | self.ADJOINT_TIMES

    def apply(self, x, mode):
        self._check_input(x, mode)
163
164
165
166

        #Maybe make class properties
        axis = self._target.axes[self._space][0]
        sl = (slice(None),)*axis
Philipp Haim's avatar
Fixes    
Philipp Haim committed
167
        extender_sl = (None,)*axis + (slice(None),) + (None,)*(self._target.axes[-1][-1] - axis)
168
169
170
171
172
        first = sl + (0,)
        second = sl + (1,)
        from_third = sl + (slice(2,None),)
        no_border = sl + (slice(1,-1),)
        reverse = sl + (slice(None,None,-1),)
173
174

        x = x.to_global_data_rw()
Philipp Arras's avatar
Philipp Arras committed
175
176
        if mode == self.TIMES:
            res = np.empty(self._target.shape)
177
            res[first] = res[second] = 0
178
            res[from_third] = np.cumsum(x[second], axis = axis)
179
            res[from_third] = (res[from_third] + res[no_border])/2*self._log_vol[extender_sl] + x[first]
180
            res[from_third] = np.cumsum(res[from_third], axis = axis)
Philipp Arras's avatar
Philipp Arras committed
181
182
        else:
            res = np.zeros(self._domain.shape)
183
184
            x[from_third] = np.cumsum(x[from_third][reverse], axis = axis)[reverse]
            res[first] += x[from_third]
185
            x[from_third] *= (self._log_vol/2.)[extender_sl]
186
187
188
            x[no_border] += x[from_third]
            res[second] += np.cumsum(x[from_third][reverse], axis = axis)[reverse]
        return from_global_data(self._tgt(mode), res)
Philipp Arras's avatar
Philipp Arras committed
189
190
191


class _Normalization(Operator):
192
    def __init__(self, domain, space = 0):
Philipp Arras's avatar
Philipp Arras committed
193
        self._domain = self._target = makeDomain(domain)
194
        assert isinstance(self._domain[space], PowerSpace)
195
196
197
198
        hspace = list(self._domain)
        hspace[space] = hspace[space].harmonic_partner
        hspace = makeDomain(hspace)
        pd = PowerDistributor(hspace, power_space=self._domain[space], space = space)
Philipp Arras's avatar
Philipp Arras committed
199
        # TODO Does not work on sphere yet
200
        mode_multiplicity = pd.adjoint(full(pd.target, 1.)).to_global_data_rw()
201
        zero_mode = (slice(None),)*self._domain.axes[space][0] + (0,)
Philipp Haim's avatar
Philipp Haim committed
202
        mode_multiplicity[zero_mode] = 0
203
204
        self._mode_multiplicity = from_global_data(self._domain, mode_multiplicity)
        self._specsum = _SpecialSum(self._domain, space)
Philipp Arras's avatar
Philipp Arras committed
205
206
207
208
209
210
211

    def apply(self, x):
        self._check_input(x)
        amp = x.exp()
        spec = (2*x).exp()
        # FIXME This normalizes also the zeromode which is supposed to be left
        # untouched by this operator
212
        return self._specsum(self._mode_multiplicity*spec)**(-0.5)*amp
Philipp Arras's avatar
Philipp Arras committed
213
214
215


class _SpecialSum(EndomorphicOperator):
216
    def __init__(self, domain, space = 0):
Philipp Arras's avatar
Philipp Arras committed
217
218
        self._domain = makeDomain(domain)
        self._capability = self.TIMES | self.ADJOINT_TIMES
219
        self._contractor = ContractionOperator(domain, space)
Philipp Arras's avatar
Philipp Arras committed
220
221
222

    def apply(self, x, mode):
        self._check_input(x, mode)
223
        return self._contractor.adjoint(self._contractor(x))
Philipp Arras's avatar
Philipp Arras committed
224
225


226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
class _slice_extractor(LinearOperator):
    #FIXME it should be tested if the the domain and target are consistent with the slice
    def __init__(self, domain, target, sl):
        self._domain = makeDomain(domain)
        self._target = makeDomain(target)
        self._sl = sl
        self._capability = self.TIMES | self.ADJOINT_TIMES

    def apply(self, x, mode):
        self._check_input(x, mode)
        x = x.to_global_data()
        if mode == self.TIMES:
            res = x[self._sl]
            res = res.reshape(self._target.shape)
        else:
            res = np.zeros(self._domain.shape)
            res[self._sl] = x
        return from_global_data(self._tgt(mode), res)
    

246
247
class _Amplitude(Operator):
    def __init__(self, target, fluctuations, flexibility, asperity,
248
                 loglogavgslope, key, space = 0):
Philipp Arras's avatar
Philipp Arras committed
249
250
251
252
253
254
255
256
257
258
259
        """
        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)
        target = makeDomain(target)
260
        assert isinstance(target[space], PowerSpace)
Philipp Arras's avatar
Philipp Arras committed
261

262
        target = makeDomain(target)
263
        twolog = _TwoLogIntegrations(target, space)
Philipp Arras's avatar
Philipp Arras committed
264
        dom = twolog.domain
265
266
267
268
        shp = dom[space].shape
        totvol = target[space].harmonic_partner.get_default_codomain().total_volume
        expander = ContractionOperator(dom, spaces = space).adjoint
        ps_expander = ContractionOperator(twolog.target, spaces = space).adjoint
Philipp Arras's avatar
Philipp Arras committed
269
270
271

        # Prepare constant fields
        foo = np.zeros(shp)
272
273
        foo[0] = foo[1] = np.sqrt(_log_vol(target[space]))
        vflex = DiagonalOperator(from_global_data(dom[space], foo), dom, space)
Philipp Arras's avatar
Philipp Arras committed
274
275
276

        foo = np.zeros(shp, dtype=np.float64)
        foo[0] += 1
277
        vasp = DiagonalOperator(from_global_data(dom[space], foo), dom, space)
Philipp Arras's avatar
Philipp Arras committed
278
279

        foo = np.ones(shp)
280
281
        foo[0] = _log_vol(target[space])**2/12.
        shift = DiagonalOperator(from_global_data(dom[space], foo), dom, space)
Philipp Frank's avatar
Philipp Frank committed
282
        
283
284
285
286
287
        vslope = DiagonalOperator(
                    from_global_data(target[space], _relative_log_k_lengths(target[space])),
                    target, space)

        foo, bar = [np.zeros(target[space].shape) for _ in range(2)]
Philipp Arras's avatar
Philipp Arras committed
288
        bar[1:] = foo[0] = totvol
289
290
291
292
293
294
295
        vol0, vol1 = [DiagonalOperator(from_global_data(target[space], aa), 
                target, space) for aa in (foo, bar)]

        #Prepare fields for Adder
        #NOTE alternative would be adjoint contraction_operator acting
        #on every space except the specified on
        shift, vol0 = [op(full(op.domain, 1)) for op in (shift, vol0)]
Philipp Arras's avatar
Philipp Arras committed
296
297
        # End prepare constant fields

298
299
300
301
        slope = vslope @ ps_expander @ loglogavgslope
        sig_flex = vflex @ expander @ flexibility
        sig_asp = vasp @ expander @ asperity
        sig_fluc = vol1 @ ps_expander @ fluctuations
Philipp Arras's avatar
Philipp Arras committed
302
303

        xi = ducktape(dom, None, key)
Philipp Arras's avatar
Philipp Arras committed
304
        sigma = sig_flex*(Adder(shift) @ sig_asp).sqrt()
305
306
        smooth = _SlopeRemover(target, space) @ twolog @ (sigma*xi)
        op = _Normalization(target, space) @ (slope + smooth)
Philipp Arras's avatar
Philipp Arras committed
307
308
        op = Adder(vol0) @ (sig_fluc*op)

Philipp Arras's avatar
Philipp Arras committed
309
        self.apply = op.apply
Philipp Arras's avatar
Philipp Arras committed
310
        self._fluc = fluctuations
Philipp Arras's avatar
Philipp Arras committed
311
        self._domain, self._target = op.domain, op.target
312
        self._space = space
Philipp Arras's avatar
Philipp Arras committed
313

Philipp Arras's avatar
Philipp Arras committed
314
315
316
317
    @property
    def fluctuation_amplitude(self):
        return self._fluc

318
319

class CorrelatedFieldMaker:
320
    def __init__(self, amplitude_offset, prefix):
321
        self._a = []
322
        self._spaces = []
323
        self._position_spaces = []
324
325
326
327
328
329
330
331
332
333
334
335
336
337
        
        self._azm = amplitude_offset
        self._prefix = prefix
    
    @staticmethod
    def make(offset_amplitude_mean, offset_amplitude_stddev, prefix):
        offset_amplitude_stddev = float(offset_amplitude_stddev)
        offset_amplitude_mean = float(offset_amplitude_mean)
        assert offset_amplitude_stddev > 0
        assert offset_amplitude_mean > 0
        zm = _LognormalMomentMatching(offset_amplitude_mean,
                                      offset_amplitude_stddev,
                                      prefix + 'zeromode')
        return CorrelatedFieldMaker(zm, prefix)
338
339

    def add_fluctuations(self,
340
                         position_space,
341
342
343
344
345
346
347
348
                         fluctuations_mean,
                         fluctuations_stddev,
                         flexibility_mean,
                         flexibility_stddev,
                         asperity_mean,
                         asperity_stddev,
                         loglogavgslope_mean,
                         loglogavgslope_stddev,
349
                         prefix='',
350
351
352
353
354
355
                         index=None,
                         space=0):
        position_space = makeDomain(position_space)
        power_space = list(position_space)
        power_space[space] = PowerSpace(position_space[space].get_default_codomain())
        power_space = makeDomain(power_space)
Philipp Arras's avatar
Philipp Arras committed
356
        prefix = str(prefix)
357
        #assert isinstance(position_space[space], (RGSpace, HPSpace, GLSpace)
358
359
        #NOTE alternative to get auxilliary domain
        #auxdom = ContractionOperator(position_space, space).domain
360
361
        auxdom = makeDomain(tuple(dom for i, dom in enumerate(position_space)
                                    if i != space))
Philipp Arras's avatar
Philipp Arras committed
362

Philipp Arras's avatar
Philipp Arras committed
363
364
        fluct = _LognormalMomentMatching(fluctuations_mean,
                                         fluctuations_stddev,
365
366
                                         prefix + 'fluctuations',
                                         auxdom)
367
368
        #FIXME How should this work on domain tuples?
        #fluct = fluct*self._azm.one_over()
Philipp Arras's avatar
Philipp Arras committed
369
        flex = _LognormalMomentMatching(flexibility_mean, flexibility_stddev,
370
371
                                        prefix + 'flexibility',
                                        auxdom)
Philipp Arras's avatar
Philipp Arras committed
372
        asp = _LognormalMomentMatching(asperity_mean, asperity_stddev,
373
374
                                       prefix + 'asperity', 
                                       auxdom)
375
        avgsl = _normal(loglogavgslope_mean, loglogavgslope_stddev,
376
377
378
                        prefix + 'loglogavgslope', auxdom)
        amp = _Amplitude(power_space,
                         fluct, flex, asp, avgsl, prefix + 'spectrum', space)
379
380
        if index is not None:
            self._a.insert(index, amp)
381
            self._position_spaces.insert(index, position_space)
382
            self._spaces.insert(index, space)
383
384
        else:
            self._a.append(amp)
385
            self._position_spaces.append(position_space)
Philipp Haim's avatar
Typo    
Philipp Haim committed
386
            self._spaces.append(space)
387

388
    def finalize_from_op(self, zeromode, prefix='', space = 0):
389
        assert isinstance(zeromode, Operator)
390
        self._azm = zeromode
391
392
393
394
395
396
        hspace = []
        tuple(hspace.extend(tuple(get_default_codomain(dd, space)))
                for dd, space in zip(self._position_spaces, self._spaces))
        hspace = makeDomain(hspace)
        zeroind = ()
        for i, dd in enumerate(self._position_spaces):
Philipp Haim's avatar
Philipp Haim committed
397
               zeroind += (slice(None),)*(self._spaces[i])
398
               zeroind += (0,)*len(dd[self._spaces[i]].shape) 
Philipp Haim's avatar
Philipp Haim committed
399
400
               zeroind += (slice(None),)*(len(dd)-self._spaces[i]-1)

401
402
        foo = np.ones(hspace.shape)
        foo[zeroind] = 0
403

404
405
406
        ZeroModeInserter = _slice_extractor(hspace, 
                     self._azm.target, zeroind).adjoint
        azm = Adder(from_global_data(hspace, foo)) @ ZeroModeInserter @ zeromode
407
408

        n_amplitudes = len(self._a)
Philipp Haim's avatar
Philipp Haim committed
409
410
411
412
413
414
        spaces = [self._spaces[0],]
        for i in range(1,n_amplitudes):
               spaces.extend(
                       [len(self._position_spaces[i-1]) 
                       - self._spaces[i-1] + self._spaces[i]])
        spaces = list(np.cumsum(spaces))
415
416
        ht = HarmonicTransformOperator(hspace,
                                   self._position_spaces[0][self._spaces[0]],
Philipp Haim's avatar
Philipp Haim committed
417
                                   space=spaces[0])
418
        for i in range(1, n_amplitudes):
419
            ht = (HarmonicTransformOperator(ht.target,
420
                                    self._position_spaces[i][self._spaces[i]],
Philipp Haim's avatar
Philipp Haim committed
421
                                    space=spaces[i]) @ ht)
422

423
        pd = PowerDistributor(hspace, self._a[0].target[self._spaces[0]], self._spaces[0])
424
        for i in range(1, n_amplitudes):
Philipp Haim's avatar
Philipp Haim committed
425
            pd = (pd @ PowerDistributor(pd.domain,
426
                                   self._a[i].target[self._spaces[i]],
Philipp Haim's avatar
Philipp Haim committed
427
                                   space=spaces[i]))
Philipp Arras's avatar
Philipp Arras committed
428

Philipp Haim's avatar
Philipp Haim committed
429
        all_spaces = list(range(len(hspace)))
430
431
        a = ContractionOperator(pd.domain, spaces[1:]).adjoint @ self._a[0]
        for i in range(1, n_amplitudes):
Philipp Haim's avatar
Philipp Haim committed
432
433
            co = ContractionOperator(pd.domain,
                    all_spaces[:spaces[i]] + all_spaces[spaces[i] + 1:])
434
            a = a*(co.adjoint @ self._a[i])
Philipp Arras's avatar
Philipp Arras committed
435

436
        return ht(azm*(pd @ a)*ducktape(hspace, None, prefix + 'xi'))
Philipp Arras's avatar
Philipp Arras committed
437
438

    def finalize(self,
439
440
                 offset=None,
                 prior_info=100):
Philipp Arras's avatar
Philipp Arras committed
441
442
443
444
        """
        offset vs zeromode: volume factor
        """
        if offset is not None:
445
            raise NotImplementedError
Philipp Arras's avatar
Philipp Arras committed
446
            offset = float(offset)
447

448
        op = self.finalize_from_op(self._azm, self._prefix)
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
        if prior_info > 0:
            from ..sugar import from_random
            samps = [
                from_random('normal', op.domain) for _ in range(prior_info)
            ]
            self.statistics_summary(samps)
        return op

    def statistics_summary(self, samples):
        lst = [('Offset amplitude', self.amplitude_total_offset),
               ('Total fluctuation amplitude', self.total_fluctuation)]

        namps = len(self.amplitudes)
        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:
            mean, stddev = _stats(op, samples)
471
472
            for m, s in zip(mean.flatten(), stddev.flatten()):
                print('{}: {:.02E} ± {:.02E}'.format(kk, m, s))
473
474
475
476

    def moment_slice_to_average(self, fluctuations_slice_mean, nsamples=1000):
        fluctuations_slice_mean = float(fluctuations_slice_mean)
        assert fluctuations_slice_mean > 0
477
        from ..sugar import from_random
478
479
        scm = 1.
        for a in self._a:
480
481
482
483
            op = a.fluctuation_amplitude
            res= np.array([op(from_random('normal',op.domain)).to_global_data()
                            for _ in range(nsamples)])
            scm *= res**2 + 1.
484
        return fluctuations_slice_mean/np.mean(np.sqrt(scm))
485

Philipp Arras's avatar
Philipp Arras committed
486
487
    @property
    def amplitudes(self):
488
        return self._a
Philipp Arras's avatar
Philipp Arras committed
489

490
491
492
    @property
    def amplitude_total_offset(self):
        return self._azm
Philipp Arras's avatar
Philipp Arras committed
493
494

    @property
495
    def total_fluctuation(self):
496
        """Returns operator which acts on prior or posterior samples"""
497
        if len(self._a) == 0:
Philipp Arras's avatar
Philipp Arras committed
498
            raise NotImplementedError
499
        if len(self._a) == 1:
500
            return self.average_fluctuation(0)
501
502
503
        q = 1.
        for a in self._a:
            fl = a.fluctuation_amplitude
Philipp Arras's avatar
Philipp Arras committed
504
            q = q*(Adder(full(fl.target, 1.)) @ fl**2)
505
        return (Adder(full(q.target, -1.)) @ q).sqrt() * self._azm
506

Philipp Arras's avatar
Philipp Arras committed
507
    def slice_fluctuation(self, space):
508
        """Returns operator which acts on prior or posterior samples"""
509
        if len(self._a) == 0:
Philipp Arras's avatar
Philipp Arras committed
510
            raise NotImplementedError
511
512
        assert space < len(self._a)
        if len(self._a) == 1:
513
            return self.average_fluctuation(0)
514
515
516
517
        q = 1.
        for j in range(len(self._a)):
            fl = self._a[j].fluctuation_amplitude
            if j == space:
Philipp Arras's avatar
Philipp Arras committed
518
                q = q*fl**2
519
            else:
Philipp Arras's avatar
Philipp Arras committed
520
                q = q*(Adder(full(fl.target, 1.)) @ fl**2)
521
        return q.sqrt() * self._azm
Philipp Arras's avatar
Philipp Arras committed
522
523

    def average_fluctuation(self, space):
524
        """Returns operator which acts on prior or posterior samples"""
525
        if len(self._a) == 0:
Philipp Arras's avatar
Philipp Arras committed
526
            raise NotImplementedError
527
528
        assert space < len(self._a)
        if len(self._a) == 1:
529
530
            return self._a[0].fluctuation_amplitude*self._azm
        return self._a[space].fluctuation_amplitude*self._azm
531

532
533
    @staticmethod
    def offset_amplitude_realized(samples):
534
535
        res = 0.
        for s in samples:
Philipp Frank's avatar
fixes    
Philipp Frank committed
536
            res = res + s.mean()**2
537
        return np.sqrt(res/len(samples))
Philipp Arras's avatar
Philipp Arras committed
538

539
540
541
542
543
544
545
546
    @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."""
547
548
549
        ldom = len(samples[0].domain)
        assert space < ldom
        if ldom == 1:
550
            return _total_fluctuation_realized(samples)
Philipp Arras's avatar
Philipp Arras committed
551
        res1, res2 = 0., 0.
552
        for s in samples:
Philipp Frank's avatar
fixes    
Philipp Frank committed
553
554
555
556
557
558
559
            res1 = res1 + s**2
            res2 = res2 + s.mean(space)**2
        res1 = res1/len(samples)
        res2 = res2/len(samples)
        res = res1.mean() - res2.mean()
        return np.sqrt(res)

560

Philipp Arras's avatar
Philipp Arras committed
561
    @staticmethod
562
563
564
565
566
567
568
569
570
571
572
    def average_fluctuation_realized(samples, space):
        """Computes average fluctuations from collection of field (defined in signal
        space) realizations."""
        ldom = len(samples[0].domain)
        assert space < ldom
        if ldom == 1:
            return _total_fluctuation_realized(samples)
        spaces = ()
        for i in range(ldom):
            if i != space:
                spaces += (i,)
Philipp Arras's avatar
Philipp Arras committed
573
574
        res = 0.
        for s in samples:
575
576
577
578
            r = s.mean(spaces)
            res = res + (r - r.mean())**2
        res = res/len(samples)
        return np.sqrt(res.mean())