correlated_fields.py 21.3 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
62
    assert np.all(sig > 0)
    return Adder(from_global_data(domain, mean)) @ (
63
64
        DiagonalOperator(from_global_data(domain,sig))
        @ ducktape(domain, None, key))
Philipp Arras's avatar
Philipp Arras committed
65
66


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


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


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


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


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
105
class _LognormalMomentMatching(Operator):
106
107
    def __init__(self, mean, sig, key,
            domain = DomainTuple.scalar_domain()):
Philipp Arras's avatar
Philipp Arras committed
108
        key = str(key)
109
        logmean, logsig = _lognormal_moments(mean, sig, domain.shape)
Philipp Arras's avatar
Philipp Arras committed
110
111
        self._mean = mean
        self._sig = sig
112
        op = _normal(logmean, logsig, key, domain).exp()
Philipp Arras's avatar
Philipp Arras committed
113
114
115
116
117
118
119
120
121
122
        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
123
124


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

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

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

Philipp Arras's avatar
Philipp Arras committed
149
150

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

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

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

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


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

    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
213
        return self._specsum(self._mode_multiplicity*spec)**(-0.5)*amp
Philipp Arras's avatar
Philipp Arras committed
214
215
216


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

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


227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
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)
    

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

263
        target = makeDomain(target)
264
        twolog = _TwoLogIntegrations(target, space)
Philipp Arras's avatar
Philipp Arras committed
265
        dom = twolog.domain
266
267
268
269
        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
270
271
272

        # Prepare constant fields
        foo = np.zeros(shp)
273
274
        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
275
276
277

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

        foo = np.ones(shp)
281
282
        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
283
        
284
285
286
287
288
        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
289
        bar[1:] = foo[0] = totvol
290
291
292
293
294
295
296
        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
297
298
        # End prepare constant fields

299
300
301
302
        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
303
304

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

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

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

319
320
321
322

class CorrelatedFieldMaker:
    def __init__(self):
        self._a = []
323
        self._spaces = []
324
        self._azm = None
325
        self._position_spaces = []
326
327

    def add_fluctuations(self,
328
                         position_space,
329
330
331
332
333
334
335
336
                         fluctuations_mean,
                         fluctuations_stddev,
                         flexibility_mean,
                         flexibility_stddev,
                         asperity_mean,
                         asperity_stddev,
                         loglogavgslope_mean,
                         loglogavgslope_stddev,
337
                         prefix='',
338
339
340
341
342
343
                         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
344
        prefix = str(prefix)
345
346
347
        #assert isinstance(position_space[space], (RGSpace, HPSpace, GLSpace)
        auxdom = makeDomain(tuple(dom for i, dom in enumerate(position_space)
                                    if i != space))
Philipp Arras's avatar
Philipp Arras committed
348

Philipp Arras's avatar
Philipp Arras committed
349
350
        fluct = _LognormalMomentMatching(fluctuations_mean,
                                         fluctuations_stddev,
351
352
                                         prefix + 'fluctuations',
                                         auxdom)
Philipp Arras's avatar
Philipp Arras committed
353
        flex = _LognormalMomentMatching(flexibility_mean, flexibility_stddev,
354
355
                                        prefix + 'flexibility',
                                        auxdom)
Philipp Arras's avatar
Philipp Arras committed
356
        asp = _LognormalMomentMatching(asperity_mean, asperity_stddev,
357
358
                                       prefix + 'asperity', 
                                       auxdom)
359
        avgsl = _normal(loglogavgslope_mean, loglogavgslope_stddev,
360
361
362
                        prefix + 'loglogavgslope', auxdom)
        amp = _Amplitude(power_space,
                         fluct, flex, asp, avgsl, prefix + 'spectrum', space)
363
364
        if index is not None:
            self._a.insert(index, amp)
365
            self._position_spaces.insert(index, position_space)
366
            self._spaces.insert(index, space)
367
368
        else:
            self._a.append(amp)
369
            self._position_spaces.append(position_space)
370
            self._spaces.append(spaces)
371

372
    def finalize_from_op(self, zeromode, prefix='', space = 0):
373
        assert isinstance(zeromode, Operator)
374
        self._azm = zeromode
375
376
377
378
379
380
381
382
383
        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):
               zeroind += (slice(None),)*self._spaces[i]
               zeroind += (0,)*len(dd[self._spaces[i]].shape) 
        #tuple(zeroind = zeroind + (slice(None),)*space +  len(dd.shape)*(0,)
384
385
        foo = np.ones(hspace.shape)
        foo[zeroind] = 0
386

387
388
389
        ZeroModeInserter = _slice_extractor(hspace, 
                     self._azm.target, zeroind).adjoint
        azm = Adder(from_global_data(hspace, foo)) @ ZeroModeInserter @ zeromode
390
391

        n_amplitudes = len(self._a)
392
393
394
        ht = HarmonicTransformOperator(hspace,
                                   self._position_spaces[0][self._spaces[0]],
                                   space=self._spaces[0])
395
        for i in range(1, n_amplitudes):
396
            ht = (HarmonicTransformOperator(ht.target,
397
398
                                    self._position_spaces[i][self._spaces[i]],
                                    space=self._spaces[i]) @ ht)
399

400
        pd = PowerDistributor(hspace, self._a[0].target[self._spaces[0]], self._spaces[0])
401
        for i in range(1, n_amplitudes):
402
403
404
            pd = (PowerDistributor(pd.domain,
                                   self._a[i].target[self._spaces[i]],
                                   space=self._spaces[i]) @ pd)
Philipp Arras's avatar
Philipp Arras committed
405

406
        spaces = np.cumsum(self._spaces) + np.arange(len(self._spaces))
407
408
409
410
        a = ContractionOperator(pd.domain, spaces[1:]).adjoint @ self._a[0]
        for i in range(1, n_amplitudes):
            co = ContractionOperator(pd.domain, spaces[:i] + spaces[(i + 1):])
            a = a*(co.adjoint @ self._a[i])
Philipp Arras's avatar
Philipp Arras committed
411

412
        return ht(azm*(pd @ a)*ducktape(hspace, None, prefix + 'xi'))
Philipp Arras's avatar
Philipp Arras committed
413
414
415
416

    def finalize(self,
                 offset_amplitude_mean,
                 offset_amplitude_stddev,
417
                 prefix='',
418
419
                 offset=None,
                 prior_info=100):
Philipp Arras's avatar
Philipp Arras committed
420
421
422
        """
        offset vs zeromode: volume factor
        """
423
        prefix = str(prefix)
Philipp Arras's avatar
Philipp Arras committed
424
        if offset is not None:
425
            raise NotImplementedError
Philipp Arras's avatar
Philipp Arras committed
426
            offset = float(offset)
427
428
429
430

        auxdom = makeDomain((dom for i, domT in enumerate(self._position_spaces)
                                    for j, dom in enumerate(domT)
                                        if j != self._spaces[i]))
Philipp Arras's avatar
Philipp Arras committed
431
432
        azm = _LognormalMomentMatching(offset_amplitude_mean,
                                       offset_amplitude_stddev,
433
434
                                       prefix + 'zeromode',
                                       auxdom)
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
        op = self.finalize_from_op(azm, prefix)
        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)
458
459
            for m, s in zip(mean.flatten(), stddev.flatten()):
                print('{}: {:.02E} ± {:.02E}'.format(kk, m, s))
460
461
462
463
464
465
466
467
468
469
470

    def moment_slice_to_average(self, fluctuations_slice_mean, nsamples=1000):
        fluctuations_slice_mean = float(fluctuations_slice_mean)
        assert fluctuations_slice_mean > 0
        scm = 1.
        for a in self._a:
            m, std = a.fluctuation_amplitude.mean, a.fluctuation_amplitude.std
            mu, sig = _lognormal_moments(m, std)
            flm = np.exp(mu + sig*np.random.normal(size=nsamples))
            scm *= flm**2 + 1.
        return fluctuations_slice_mean/np.mean(np.sqrt(scm))
471

Philipp Arras's avatar
Philipp Arras committed
472
473
    @property
    def amplitudes(self):
474
        return self._a
Philipp Arras's avatar
Philipp Arras committed
475

476
477
478
    @property
    def amplitude_total_offset(self):
        return self._azm
Philipp Arras's avatar
Philipp Arras committed
479
480

    @property
481
    def total_fluctuation(self):
482
        """Returns operator which acts on prior or posterior samples"""
483
        if len(self._a) == 0:
Philipp Arras's avatar
Philipp Arras committed
484
            raise NotImplementedError
485
486
487
488
489
        if len(self._a) == 1:
            return self._a[0].fluctuation_amplitude
        q = 1.
        for a in self._a:
            fl = a.fluctuation_amplitude
Philipp Arras's avatar
Philipp Arras committed
490
491
            q = q*(Adder(full(fl.target, 1.)) @ fl**2)
        return (Adder(full(q.target, -1.)) @ q).sqrt()
492

Philipp Arras's avatar
Philipp Arras committed
493
    def slice_fluctuation(self, space):
494
        """Returns operator which acts on prior or posterior samples"""
495
        if len(self._a) == 0:
Philipp Arras's avatar
Philipp Arras committed
496
            raise NotImplementedError
497
498
499
500
501
502
503
        assert space < len(self._a)
        if len(self._a) == 1:
            return self._a[0].fluctuation_amplitude
        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
504
                q = q*fl**2
505
            else:
Philipp Arras's avatar
Philipp Arras committed
506
                q = q*(Adder(full(fl.target, 1.)) @ fl**2)
507
        return q.sqrt()
Philipp Arras's avatar
Philipp Arras committed
508
509

    def average_fluctuation(self, space):
510
        """Returns operator which acts on prior or posterior samples"""
511
        if len(self._a) == 0:
Philipp Arras's avatar
Philipp Arras committed
512
            raise NotImplementedError
513
514
515
516
517
        assert space < len(self._a)
        if len(self._a) == 1:
            return self._a[0].fluctuation_amplitude
        return self._a[space].fluctuation_amplitude

518
519
    @staticmethod
    def offset_amplitude_realized(samples):
520
521
        res = 0.
        for s in samples:
Philipp Frank's avatar
fixes    
Philipp Frank committed
522
            res = res + s.mean()**2
523
        return np.sqrt(res/len(samples))
Philipp Arras's avatar
Philipp Arras committed
524

525
526
527
528
529
530
531
532
    @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."""
533
534
535
        ldom = len(samples[0].domain)
        assert space < ldom
        if ldom == 1:
536
            return _total_fluctuation_realized(samples)
Philipp Arras's avatar
Philipp Arras committed
537
        res1, res2 = 0., 0.
538
        for s in samples:
Philipp Frank's avatar
fixes    
Philipp Frank committed
539
540
541
542
543
544
545
            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)

546

Philipp Arras's avatar
Philipp Arras committed
547
    @staticmethod
548
549
550
551
552
553
554
555
556
557
558
    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
559
560
        res = 0.
        for s in samples:
561
562
563
564
            r = s.mean(spaces)
            res = res + (r - r.mean())**2
        res = res/len(samples)
        return np.sqrt(res.mean())