correlated_fields.py 29.2 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 ..domain_tuple import DomainTuple
Philipp Arras's avatar
Philipp Arras committed
25
26
from ..domains.power_space import PowerSpace
from ..domains.unstructured_domain import UnstructuredDomain
Philipp Arras's avatar
Philipp Arras committed
27
from ..field import Field
28
from ..logger import logger
Philipp Arras's avatar
Philipp Arras committed
29
from ..multi_field import MultiField
Philipp Arras's avatar
Philipp Arras committed
30
from ..operators.adder import Adder
31
from ..operators.contraction_operator import ContractionOperator
Philipp Arras's avatar
Philipp Arras committed
32
from ..operators.diagonal_operator import DiagonalOperator
Philipp Arras's avatar
Philipp Arras committed
33
from ..operators.distributors import PowerDistributor
Philipp Arras's avatar
Philipp Arras committed
34
from ..operators.endomorphic_operator import EndomorphicOperator
Martin Reinecke's avatar
Martin Reinecke committed
35
from ..operators.harmonic_operators import HarmonicTransformOperator
Philipp Arras's avatar
Philipp Arras committed
36
from ..operators.linear_operator import LinearOperator
Philipp Arras's avatar
Philipp Arras committed
37
from ..operators.operator import Operator
Philipp Arras's avatar
Philipp Arras committed
38
from ..operators.simple_linear_operators import ducktape
39
from ..probing import StatCalculator
Philipp Arras's avatar
Philipp Arras committed
40
from ..sugar import full, makeDomain, makeField, makeOp
41

42

Philipp Haim's avatar
Philipp Haim committed
43
44
45
def _reshaper(x, N):
    x = np.asfarray(x)
    if x.shape in [(), (1,)]:
Philipp Haim's avatar
Philipp Haim committed
46
        return np.full(N, x) if N != 0 else x.reshape(())
Philipp Haim's avatar
Philipp Haim committed
47
48
    elif x.shape == (N,):
        return x
49
50
51
    else:
        raise TypeError("Shape of parameters cannot be interpreted")

52

Martin Reinecke's avatar
Martin Reinecke committed
53
def _lognormal_moments(mean, sig, N=0):
Philipp Haim's avatar
Philipp Haim committed
54
55
56
57
    if N == 0:
        mean, sig = np.asfarray(mean), np.asfarray(sig)
    else:
        mean, sig = (_reshaper(param, N) for param in (mean, sig))
58
    if not np.all(mean > 0):
59
        raise ValueError("mean must be greater 0; got {!r}".format(mean))
60
    if not np.all(sig > 0):
61
        raise ValueError("sig must be greater 0; got {!r}".format(sig))
62

63
    logsig = np.sqrt(np.log1p((sig/mean)**2))
Philipp Arras's avatar
Philipp Arras committed
64
    logmean = np.log(mean) - logsig**2/2
65
    return logmean, logsig
Philipp Arras's avatar
Philipp Arras committed
66
67


Martin Reinecke's avatar
Martin Reinecke committed
68
def _normal(mean, sig, key, N=0):
Philipp Haim's avatar
Philipp Haim committed
69
    if N == 0:
Philipp Haim's avatar
Philipp Haim committed
70
        domain = DomainTuple.scalar_domain()
Philipp Haim's avatar
Philipp Haim committed
71
        mean, sig = np.asfarray(mean), np.asfarray(sig)
72
73
74
75
76
        return Adder(makeField(domain, mean)) @ (
            sig * ducktape(domain, None, key))

    domain = UnstructuredDomain(N)
    mean, sig = (_reshaper(param, N) for param in (mean, sig))
Martin Reinecke's avatar
Martin Reinecke committed
77
78
    return Adder(makeField(domain, mean)) @ (DiagonalOperator(
        makeField(domain, sig)) @ ducktape(domain, None, key))
Philipp Arras's avatar
Philipp Arras committed
79
80


Philipp Arras's avatar
Philipp Arras committed
81
def _log_k_lengths(pspace):
Philipp Arras's avatar
Philipp Arras committed
82
    """Log(k_lengths) without zeromode"""
Philipp Arras's avatar
Philipp Arras committed
83
84
85
    return np.log(pspace.k_lengths[1:])


Philipp Arras's avatar
Philipp Arras committed
86
def _relative_log_k_lengths(power_space):
Philipp Arras's avatar
Philipp Arras committed
87
88
    """Log-distance to first bin
    logkl.shape==power_space.shape, logkl[0]=logkl[1]=0"""
Philipp Arras's avatar
Philipp Arras committed
89
90
91
92
93
94
    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
95
    return np.insert(logkl, 0, 0)
Philipp Arras's avatar
Philipp Arras committed
96
97


Philipp Arras's avatar
Philipp Arras committed
98
def _log_vol(power_space):
99
    power_space = makeDomain(power_space)
Philipp Arras's avatar
Philipp Arras committed
100
101
102
103
104
    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
105
106
107
108
109
110
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
111
def _total_fluctuation_realized(samples):
Philipp Haim's avatar
Philipp Haim committed
112
113
114
    spaces = _structured_spaces(samples[0].domain)
    co = ContractionOperator(samples[0].domain, spaces)
    size = co.domain.size/co.target.size
115
116
    res = 0.
    for s in samples:
Philipp Haim's avatar
Philipp Haim committed
117
118
        res = res + (s - co.adjoint(co(s)/size))**2
    res = res.mean(spaces)/len(samples)
Philipp Haim's avatar
Philipp Haim committed
119
    return np.sqrt(res if np.isscalar(res) else res.val)
120
121


Philipp Arras's avatar
Philipp Arras committed
122
class _LognormalMomentMatching(Operator):
Philipp Haim's avatar
Philipp Haim committed
123
    def __init__(self, mean, sig, key, N_copies):
Philipp Arras's avatar
Philipp Arras committed
124
        key = str(key)
Philipp Haim's avatar
Philipp Haim committed
125
        logmean, logsig = _lognormal_moments(mean, sig, N_copies)
Philipp Arras's avatar
Philipp Arras committed
126
127
        self._mean = mean
        self._sig = sig
Martin Reinecke's avatar
Martin Reinecke committed
128
        op = _normal(logmean, logsig, key, N_copies).ptw("exp")
Philipp Arras's avatar
Philipp Arras committed
129
130
131
132
133
134
135
136
137
138
        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
139
140


Philipp Frank's avatar
Philipp Frank committed
141
class _SlopeRemover(EndomorphicOperator):
Martin Reinecke's avatar
Martin Reinecke committed
142
    def __init__(self, domain, space=0):
Philipp Frank's avatar
Philipp Frank committed
143
        self._domain = makeDomain(domain)
144
145
        assert isinstance(self._domain[space], PowerSpace)
        logkl = _relative_log_k_lengths(self._domain[space])
146
        self._sc = logkl/float(logkl[-1])
Philipp Arras's avatar
Philipp Arras committed
147

148
        self._space = space
Philipp Haim's avatar
Philipp Haim committed
149
150
151
        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
152
        self._capability = self.TIMES | self.ADJOINT_TIMES
Philipp Arras's avatar
Philipp Arras committed
153

154
155
    def apply(self, x, mode):
        self._check_input(x, mode)
Martin Reinecke's avatar
stage 3    
Martin Reinecke committed
156
        x = x.val
Philipp Frank's avatar
Philipp Frank committed
157
        if mode == self.TIMES:
Philipp Haim's avatar
Philipp Haim committed
158
            res = x - x[self._last]*self._sc[self._extender]
Philipp Frank's avatar
Philipp Frank committed
159
        else:
160
161
            res = x.copy()
            res[self._last] -= (x*self._sc[self._extender]).sum(
Philipp Arras's avatar
Philipp Arras committed
162
                axis=self._space, keepdims=True)
Martin Reinecke's avatar
Martin Reinecke committed
163
        return makeField(self._tgt(mode), res)
Philipp Frank's avatar
Philipp Frank committed
164

Philipp Arras's avatar
Philipp Arras committed
165
166

class _TwoLogIntegrations(LinearOperator):
Martin Reinecke's avatar
Martin Reinecke committed
167
    def __init__(self, target, space=0):
Philipp Arras's avatar
Philipp Arras committed
168
        self._target = makeDomain(target)
169
170
171
172
173
        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
174
        self._log_vol = _log_vol(self._target[space])
Philipp Arras's avatar
Philipp Arras committed
175
176
177
178
        self._capability = self.TIMES | self.ADJOINT_TIMES

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

Martin Reinecke's avatar
Martin Reinecke committed
180
        # Maybe make class properties
181
182
        axis = self._target.axes[self._space][0]
        sl = (slice(None),)*axis
Philipp Haim's avatar
Fixes    
Philipp Haim committed
183
        extender_sl = (None,)*axis + (slice(None),) + (None,)*(self._target.axes[-1][-1] - axis)
184
185
        first = sl + (0,)
        second = sl + (1,)
Martin Reinecke's avatar
Martin Reinecke committed
186
187
188
        from_third = sl + (slice(2, None),)
        no_border = sl + (slice(1, -1),)
        reverse = sl + (slice(None, None, -1),)
189

Philipp Arras's avatar
Philipp Arras committed
190
        if mode == self.TIMES:
Martin Reinecke's avatar
stage 3    
Martin Reinecke committed
191
            x = x.val
Philipp Arras's avatar
Philipp Arras committed
192
            res = np.empty(self._target.shape)
193
            res[first] = res[second] = 0
Martin Reinecke's avatar
Martin Reinecke committed
194
            res[from_third] = np.cumsum(x[second], axis=axis)
195
            res[from_third] = (res[from_third] + res[no_border])/2*self._log_vol[extender_sl] + x[first]
Martin Reinecke's avatar
Martin Reinecke committed
196
            res[from_third] = np.cumsum(res[from_third], axis=axis)
Philipp Arras's avatar
Philipp Arras committed
197
        else:
Martin Reinecke's avatar
Martin Reinecke committed
198
            x = x.val_rw()
Philipp Arras's avatar
Philipp Arras committed
199
            res = np.zeros(self._domain.shape)
Martin Reinecke's avatar
Martin Reinecke committed
200
            x[from_third] = np.cumsum(x[from_third][reverse], axis=axis)[reverse]
201
            res[first] += x[from_third]
202
            x[from_third] *= (self._log_vol/2.)[extender_sl]
203
            x[no_border] += x[from_third]
Martin Reinecke's avatar
Martin Reinecke committed
204
            res[second] += np.cumsum(x[from_third][reverse], axis=axis)[reverse]
Martin Reinecke's avatar
Martin Reinecke committed
205
        return makeField(self._tgt(mode), res)
Philipp Arras's avatar
Philipp Arras committed
206
207
208


class _Normalization(Operator):
Martin Reinecke's avatar
Martin Reinecke committed
209
    def __init__(self, domain, space=0):
Philipp Arras's avatar
Philipp Arras committed
210
        self._domain = self._target = DomainTuple.make(domain)
211
        assert isinstance(self._domain[space], PowerSpace)
212
213
214
        hspace = list(self._domain)
        hspace[space] = hspace[space].harmonic_partner
        hspace = makeDomain(hspace)
Philipp Arras's avatar
Philipp Arras committed
215
216
217
        pd = PowerDistributor(hspace,
                              power_space=self._domain[space],
                              space=space)
Martin Reinecke's avatar
Martin Reinecke committed
218
        mode_multiplicity = pd.adjoint(full(pd.target, 1.)).val_rw()
219
        zero_mode = (slice(None),)*self._domain.axes[space][0] + (0,)
Philipp Haim's avatar
Philipp Haim committed
220
        mode_multiplicity[zero_mode] = 0
Philipp Arras's avatar
Philipp Arras committed
221
        self._mode_multiplicity = makeOp(makeField(self._domain, mode_multiplicity))
222
        self._specsum = _SpecialSum(self._domain, space)
Philipp Arras's avatar
Philipp Arras committed
223
224
225

    def apply(self, x):
        self._check_input(x)
Martin Reinecke's avatar
Martin Reinecke committed
226
        amp = x.ptw("exp")
227
        spec = amp**2
Philipp Arras's avatar
Philipp Arras committed
228
229
        # FIXME This normalizes also the zeromode which is supposed to be left
        # untouched by this operator
Philipp Arras's avatar
Philipp Arras committed
230
        return self._specsum(self._mode_multiplicity(spec))**(-0.5)*amp
Philipp Arras's avatar
Philipp Arras committed
231
232
233


class _SpecialSum(EndomorphicOperator):
Martin Reinecke's avatar
Martin Reinecke committed
234
    def __init__(self, domain, space=0):
Philipp Arras's avatar
Philipp Arras committed
235
236
        self._domain = makeDomain(domain)
        self._capability = self.TIMES | self.ADJOINT_TIMES
237
        self._contractor = ContractionOperator(domain, space)
Philipp Arras's avatar
Philipp Arras committed
238
239
240

    def apply(self, x, mode):
        self._check_input(x, mode)
241
        return self._contractor.adjoint(self._contractor(x))
Philipp Arras's avatar
Philipp Arras committed
242
243


Philipp Haim's avatar
Philipp Haim committed
244
class _Distributor(LinearOperator):
Lukas Platz's avatar
Lukas Platz committed
245
    def __init__(self, dofdex, domain, target):
246
247
248
        self._dofdex = np.array(dofdex)
        self._target = DomainTuple.make(target)
        self._domain = DomainTuple.make(domain)
Philipp Haim's avatar
Philipp Haim committed
249
250
251
252
        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
253
        x = x.val
Philipp Haim's avatar
Philipp Haim committed
254
255
256
        if mode == self.TIMES:
            res = x[self._dofdex]
        else:
257
            res = np.zeros(self._tgt(mode).shape, dtype=x.dtype)
Philipp Haim's avatar
Philipp Haim committed
258
            res[self._dofdex] = x
Martin Reinecke's avatar
Martin Reinecke committed
259
        return makeField(self._tgt(mode), res)
Martin Reinecke's avatar
Martin Reinecke committed
260

261

262
263
class _Amplitude(Operator):
    def __init__(self, target, fluctuations, flexibility, asperity,
264
                 loglogavgslope, azm, totvol, key, dofdex):
Philipp Arras's avatar
Philipp Arras committed
265
266
267
268
269
270
271
272
273
274
275
        """
        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
276
277
        if len(dofdex) > 0:
            N_copies = max(dofdex) + 1
Philipp Haim's avatar
Philipp Haim committed
278
            space = 1
Philipp Frank's avatar
cleanup    
Philipp Frank committed
279
280
            distributed_tgt = makeDomain((UnstructuredDomain(len(dofdex)),
                                          target))
Philipp Haim's avatar
Philipp Haim committed
281
            target = makeDomain((UnstructuredDomain(N_copies), target))
Lukas Platz's avatar
Lukas Platz committed
282
            Distributor = _Distributor(dofdex, target, distributed_tgt)
Philipp Haim's avatar
Philipp Haim committed
283
        else:
Philipp Haim's avatar
Philipp Haim committed
284
            N_copies = 0
Philipp Haim's avatar
Philipp Haim committed
285
            space = 0
Philipp Haim's avatar
Philipp Haim committed
286
            distributed_tgt = target = makeDomain(target)
Martin Reinecke's avatar
Martin Reinecke committed
287
        azm_expander = ContractionOperator(distributed_tgt, spaces=space).adjoint
Philipp Haim's avatar
Philipp Haim committed
288
        assert isinstance(target[space], PowerSpace)
Martin Reinecke's avatar
Martin Reinecke committed
289

290
        twolog = _TwoLogIntegrations(target, space)
Philipp Arras's avatar
Philipp Arras committed
291
        dom = twolog.domain
292

293
        shp = dom[space].shape
Martin Reinecke's avatar
Martin Reinecke committed
294
295
        expander = ContractionOperator(dom, spaces=space).adjoint
        ps_expander = ContractionOperator(twolog.target, spaces=space).adjoint
Philipp Arras's avatar
Philipp Arras committed
296
297
298

        # Prepare constant fields
        foo = np.zeros(shp)
299
        foo[0] = foo[1] = np.sqrt(_log_vol(target[space]))
Martin Reinecke's avatar
Martin Reinecke committed
300
        vflex = DiagonalOperator(makeField(dom[space], foo), dom, space)
Philipp Arras's avatar
Philipp Arras committed
301
302
303

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

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

310
        vslope = DiagonalOperator(
Philipp Arras's avatar
Philipp Arras committed
311
            makeField(target[space], _relative_log_k_lengths(target[space])),
Martin Reinecke's avatar
Martin Reinecke committed
312
            target, space)
313
314

        foo, bar = [np.zeros(target[space].shape) for _ in range(2)]
Philipp Arras's avatar
Philipp Arras committed
315
        bar[1:] = foo[0] = totvol
Philipp Arras's avatar
Philipp Arras committed
316
317
318
319
        vol0, vol1 = [
            DiagonalOperator(makeField(target[space], aa), target, space)
            for aa in (foo, bar)
        ]
320

Martin Reinecke's avatar
Martin Reinecke committed
321
        # Prepare fields for Adder
322
        shift, vol0 = [op(full(op.domain, 1)) for op in (shift, vol0)]
Philipp Arras's avatar
Philipp Arras committed
323
324
        # End prepare constant fields

325
326
327
328
        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
329
        sig_fluc = vol1 @ ps_expander @ fluctuations
Philipp Arras's avatar
Philipp Arras committed
330
331

        xi = ducktape(dom, None, key)
Martin Reinecke's avatar
Martin Reinecke committed
332
        sigma = sig_flex*(Adder(shift) @ sig_asp).ptw("sqrt")
333
334
        smooth = _SlopeRemover(target, space) @ twolog @ (sigma*xi)
        op = _Normalization(target, space) @ (slope + smooth)
Philipp Haim's avatar
Philipp Haim committed
335
        if N_copies > 0:
Philipp Haim's avatar
Philipp Haim committed
336
337
            op = Distributor @ op
            sig_fluc = Distributor @ sig_fluc
Martin Reinecke's avatar
Martin Reinecke committed
338
            op = Adder(Distributor(vol0)) @ (sig_fluc*(azm_expander @ azm.ptw("reciprocal"))*op)
Philipp Arras's avatar
Philipp Arras committed
339
340
            self._fluc = (_Distributor(dofdex, fluctuations.target,
                                       distributed_tgt[0]) @ fluctuations)
Philipp Haim's avatar
Philipp Haim committed
341
        else:
Martin Reinecke's avatar
Martin Reinecke committed
342
            op = Adder(vol0) @ (sig_fluc*(azm_expander @ azm.ptw("reciprocal"))*op)
Philipp Frank's avatar
fixup    
Philipp Frank committed
343
            self._fluc = fluctuations
Philipp Arras's avatar
Philipp Arras committed
344

Philipp Arras's avatar
Philipp Arras committed
345
346
        self.apply = op.apply
        self._domain, self._target = op.domain, op.target
347
        self._space = space
Philipp Arras's avatar
Philipp Arras committed
348

Philipp Arras's avatar
Philipp Arras committed
349
350
351
352
    @property
    def fluctuation_amplitude(self):
        return self._fluc

353
354

class CorrelatedFieldMaker:
Lukas Platz's avatar
Lukas Platz committed
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
    """Constrution helper for hirarchical correlated field models.

    The correlated field models are parametrized by creating
    powerspectrum operators acting on their target subdomains
    via calls to :func:`add_fluctuations`.
    During creation of the :class:`CorrelatedFieldMaker` via
    :func:`make`, a global offset from zero can be added to the
    field to be created and an operator applying gaussian fluctuations
    around this offset needs to be parametrized.

    The resulting correlated field model operator has a
    :class:`~nifty6.multi_domain.MultiDomain` as its domain and
    expects its input values to be univariately gaussian.

    The target of the constructed operator will be a
    :class:`~nifty6.domain_tuple.DomainTuple`
371
    containing the `target_subdomains` of the added fluctuations in the
Lukas Platz's avatar
Lukas Platz committed
372
373
374
375
376
377
378
    order of the `add_fluctuations` calls.

    Creation of the model operator is finished by calling the method
    :func:`finalize`, which returns the configured operator.

    See the methods :func:`make`, :func:`add_fluctuations`
    and :func:`finalize` for usage instructions."""
379
380
381
    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")
382
        self._a = []
383
        self._target_subdomains = []
Philipp Arras's avatar
Formats    
Philipp Arras committed
384

385
386
        self._offset_mean = offset_mean
        self._azm = offset_fluctuations_op
387
        self._prefix = prefix
Philipp Haim's avatar
Philipp Haim committed
388
        self._total_N = total_N
Philipp Arras's avatar
Formats    
Philipp Arras committed
389

390
    @staticmethod
Philipp Arras's avatar
Philipp Arras committed
391
    def make(offset_mean, offset_std_mean, offset_std_std, prefix, total_N=0,
Martin Reinecke's avatar
Martin Reinecke committed
392
             dofdex=None):
Lukas Platz's avatar
Lukas Platz committed
393
394
395
396
397
398
399
        """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
400
            Mean standard deviation of the offset.
Lukas Platz's avatar
Lukas Platz committed
401
        offset_std_std : float
Lukas Platz's avatar
Lukas Platz committed
402
            Standard deviation of the stddev of the offset.
Lukas Platz's avatar
Lukas Platz committed
403
404
        prefix : string
            Prefix to the names of the domains of the cf operator to be made.
Lukas Platz's avatar
Lukas Platz committed
405
            This determines the names of the operator domain.
Lukas Platz's avatar
Lukas Platz committed
406
        total_N : integer
Lukas Platz's avatar
Lukas Platz committed
407
408
409
410
            Number of copies with of the field to return.
            If not 0, the first entry of the operators target will be an
            :class:`~nifty.domains.unstructured_domain.UnstructuredDomain`
            with length `total_N`.
Philipp Arras's avatar
Philipp Arras committed
411
412
413
414
415
416
        dofdex : np.array of integers
            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
            instanciated; the first one is used for the first and second copy
            of the model and the second is used for the third.
Lukas Platz's avatar
Lukas Platz committed
417
        """
Philipp Frank's avatar
Philipp Frank committed
418
419
        if dofdex is None:
            dofdex = np.full(total_N, 0)
420
421
        elif len(dofdex) != total_N:
            raise ValueError("length of dofdex needs to match total_N")
Philipp Frank's avatar
Philipp Frank committed
422
        N = max(dofdex) + 1 if total_N > 0 else 0
Lukas Platz's avatar
Lukas Platz committed
423
424
        zm = _LognormalMomentMatching(offset_std_mean,
                                      offset_std_std,
Philipp Haim's avatar
Philipp Haim committed
425
                                      prefix + 'zeromode',
Philipp Frank's avatar
Philipp Frank committed
426
                                      N)
Philipp Frank's avatar
fixup    
Philipp Frank committed
427
        if total_N > 0:
Martin Reinecke's avatar
Martin Reinecke committed
428
            zm = _Distributor(dofdex, zm.target, UnstructuredDomain(total_N)) @ zm
429
        return CorrelatedFieldMaker(offset_mean, zm, prefix, total_N)
430
431

    def add_fluctuations(self,
432
                         target_subdomain,
433
434
435
436
437
438
439
440
                         fluctuations_mean,
                         fluctuations_stddev,
                         flexibility_mean,
                         flexibility_stddev,
                         asperity_mean,
                         asperity_stddev,
                         loglogavgslope_mean,
                         loglogavgslope_stddev,
Martin Reinecke's avatar
Martin Reinecke committed
441
442
443
444
                         prefix='',
                         index=None,
                         dofdex=None,
                         harmonic_partner=None):
Lukas Platz's avatar
Lukas Platz committed
445
446
447
448
449
450
451
        """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
        `loglogavgslope` configure the power spectrum model used on the
452
        target field subdomain `target_subdomain`.
Lukas Platz's avatar
Lukas Platz committed
453
454
455
456
457
458
459
460
461
462
463
        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
        ----------
464
465
        target_subdomain : :class:`~nifty6.domain.Domain`, \
                           :class:`~nifty6.domain_tuple.DomainTuple`
Lukas Platz's avatar
Lukas Platz committed
466
467
468
469
470
471
472
473
474
475
476
477
478
            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
            Smooth variation speed of the power spectrum
        asperity_{mean,stddev} : float
            Strength of unsmoothed power spectrum variations
            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
479
480
481
482
483
484
485
486
487
        index : int
            Position target_subdomain in the final total domain of the
            correlated field operator.
        dofdex : np.array
            An integer array specifying the amplitude 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
            instanciated; the first one is used for the first and second copy
            of the model and the second is used for the third.
Lukas Platz's avatar
Lukas Platz committed
488
489
490
491
        harmonic_partner : :class:`~nifty6.domain.Domain`, \
                           :class:`~nifty6.domain_tuple.DomainTuple`
            In which harmonic space to define the power spectrum
        """
Philipp Frank's avatar
Philipp Frank committed
492
        if harmonic_partner is None:
493
            harmonic_partner = target_subdomain.get_default_codomain()
Philipp Frank's avatar
Fixup    
Philipp Frank committed
494
        else:
495
496
            target_subdomain.check_codomain(harmonic_partner)
            harmonic_partner.check_codomain(target_subdomain)
497

Philipp Haim's avatar
Philipp Haim committed
498
499
        if dofdex is None:
            dofdex = np.full(self._total_N, 0)
500
501
        elif len(dofdex) != self._total_N:
            raise ValueError("length of dofdex needs to match total_N")
Philipp Haim's avatar
Philipp Haim committed
502

Philipp Haim's avatar
Philipp Haim committed
503
504
        if self._total_N > 0:
            N = max(dofdex) + 1
505
            target_subdomain = makeDomain((UnstructuredDomain(N), target_subdomain))
Philipp Haim's avatar
Philipp Haim committed
506
        else:
Philipp Haim's avatar
Philipp Haim committed
507
            N = 0
508
            target_subdomain = makeDomain(target_subdomain)
Philipp Arras's avatar
Philipp Arras committed
509
        prefix = str(prefix)
510
        # assert isinstance(target_subdomain[space], (RGSpace, HPSpace, GLSpace)
Philipp Arras's avatar
Philipp Arras committed
511

Philipp Arras's avatar
Philipp Arras committed
512
513
        fluct = _LognormalMomentMatching(fluctuations_mean,
                                         fluctuations_stddev,
514
                                         self._prefix + prefix + 'fluctuations',
Philipp Haim's avatar
Philipp Haim committed
515
                                         N)
Philipp Arras's avatar
Philipp Arras committed
516
        flex = _LognormalMomentMatching(flexibility_mean, flexibility_stddev,
517
                                        self._prefix + prefix + 'flexibility',
Philipp Haim's avatar
Philipp Haim committed
518
                                        N)
Philipp Arras's avatar
Philipp Arras committed
519
        asp = _LognormalMomentMatching(asperity_mean, asperity_stddev,
Philipp Arras's avatar
Philipp Arras committed
520
                                       self._prefix + prefix + 'asperity', N)
521
        avgsl = _normal(loglogavgslope_mean, loglogavgslope_stddev,
522
                        self._prefix + prefix + 'loglogavgslope', N)
Philipp Arras's avatar
Philipp Arras committed
523
        amp = _Amplitude(PowerSpace(harmonic_partner), fluct, flex, asp, avgsl,
524
                         self._azm, target_subdomain[-1].total_volume,
525
                         self._prefix + prefix + 'spectrum', dofdex)
Philipp Haim's avatar
Philipp Haim committed
526

527
528
        if index is not None:
            self._a.insert(index, amp)
529
            self._target_subdomains.insert(index, target_subdomain)
530
531
        else:
            self._a.append(amp)
532
            self._target_subdomains.append(target_subdomain)
533

Philipp Arras's avatar
Philipp Arras committed
534
535
536
537
538
539
540
541
542
543
    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
544
        n_amplitudes = len(self._a)
Philipp Haim's avatar
Philipp Haim committed
545
        if self._total_N > 0:
Philipp Arras's avatar
Philipp Arras committed
546
547
548
            hspace = makeDomain(
                [UnstructuredDomain(self._total_N)] +
                [dd.target[-1].harmonic_partner for dd in self._a])
Philipp Haim's avatar
Philipp Haim committed
549
550
            spaces = tuple(range(1, n_amplitudes + 1))
            amp_space = 1
Philipp Haim's avatar
Philipp Haim committed
551
552
        else:
            hspace = makeDomain(
Philipp Arras's avatar
Philipp Arras committed
553
                [dd.target[0].harmonic_partner for dd in self._a])
Philipp Haim's avatar
Philipp Haim committed
554
            spaces = tuple(range(n_amplitudes))
Philipp Haim's avatar
Philipp Haim committed
555
            amp_space = 0
556

Martin Reinecke's avatar
Martin Reinecke committed
557
        expander = ContractionOperator(hspace, spaces=spaces).adjoint
Philipp Frank's avatar
fixup    
Philipp Frank committed
558
        azm = expander @ self._azm
559

560
        ht = HarmonicTransformOperator(hspace,
561
                                       self._target_subdomains[0][amp_space],
Martin Reinecke's avatar
Martin Reinecke committed
562
                                       space=spaces[0])
563
        for i in range(1, n_amplitudes):
564
            ht = HarmonicTransformOperator(ht.target,
565
                                           self._target_subdomains[i][amp_space],
566
567
568
569
570
                                           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
571
            pd = PowerDistributor(co.target, pp, amp_space)
572
573
            a.append(co.adjoint @ pd @ self._a[ii])
        corr = reduce(mul, a)
Philipp Arras's avatar
Philipp Arras committed
574
        op = ht(azm*corr*ducktape(hspace, None, self._prefix + 'xi'))
Philipp Arras's avatar
Philipp Arras committed
575

576
577
        if self._offset_mean is not None:
            offset = self._offset_mean
578
579
580
581
582
583
584
            # 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
585
        self.statistics_summary(prior_info)
586
587
        return op

588
589
590
591
592
593
    def statistics_summary(self, prior_info):
        from ..sugar import from_random

        if prior_info == 0:
            return

594
595
        lst = [('Offset amplitude', self.amplitude_total_offset),
               ('Total fluctuation amplitude', self.total_fluctuation)]
596
        namps = len(self._a)
597
598
599
600
601
602
603
604
        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:
605
606
            sc = StatCalculator()
            for _ in range(prior_info):
607
                sc.add(op(from_random(op.domain, 'normal')))
Martin Reinecke's avatar
merge    
Martin Reinecke committed
608
            mean = sc.mean.val
Martin Reinecke's avatar
Martin Reinecke committed
609
            stddev = sc.var.ptw("sqrt").val
610
            for m, s in zip(mean.flatten(), stddev.flatten()):
611
                logger.info('{}: {:.02E} ± {:.02E}'.format(kk, m, s))
612
613
614

    def moment_slice_to_average(self, fluctuations_slice_mean, nsamples=1000):
        fluctuations_slice_mean = float(fluctuations_slice_mean)
615
616
617
        if not fluctuations_slice_mean > 0:
            msg = "fluctuations_slice_mean must be greater zero; got {!r}"
            raise ValueError(msg.format(fluctuations_slice_mean))
618
        from ..sugar import from_random
619
620
        scm = 1.
        for a in self._a:
Martin Reinecke's avatar
Martin Reinecke committed
621
            op = a.fluctuation_amplitude*self._azm.ptw("reciprocal")
622
            res = np.array([op(from_random(op.domain, 'normal')).val
623
624
                            for _ in range(nsamples)])
            scm *= res**2 + 1.
625
        return fluctuations_slice_mean/np.mean(np.sqrt(scm))
626

Philipp Arras's avatar
Philipp Arras committed
627
    @property
Philipp Haim's avatar
Philipp Haim committed
628
    def normalized_amplitudes(self):
629
        return self._a
Philipp Arras's avatar
Philipp Arras committed
630

Philipp Haim's avatar
Philipp Haim committed
631
632
633
634
635
636
637
    @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
638
639
        dom = self._a[0].target
        expand = ContractionOperator(dom, len(dom)-1).adjoint
Philipp Haim's avatar
Philipp Haim committed
640
641
        return self._a[0]*(expand @ self.amplitude_total_offset)

642
643
644
    @property
    def amplitude_total_offset(self):
        return self._azm
Philipp Arras's avatar
Philipp Arras committed
645
646

    @property
647
    def total_fluctuation(self):
648
        """Returns operator which acts on prior or posterior samples"""
649
        if len(self._a) == 0:
Philipp Arras's avatar
Philipp Arras committed
650
            raise NotImplementedError
651
        if len(self._a) == 1:
652
            return self.average_fluctuation(0)
653
654
        q = 1.
        for a in self._a:
Martin Reinecke's avatar
Martin Reinecke committed
655
            fl = a.fluctuation_amplitude*self._azm.ptw("reciprocal")
Philipp Arras's avatar
Philipp Arras committed
656
            q = q*(Adder(full(fl.target, 1.)) @ fl**2)
Martin Reinecke's avatar
Martin Reinecke committed
657
        return (Adder(full(q.target, -1.)) @ q).ptw("sqrt")*self._azm
658

Philipp Arras's avatar
Philipp Arras committed
659
    def slice_fluctuation(self, space):
660
        """Returns operator which acts on prior or posterior samples"""
661
        if len(self._a) == 0:
Philipp Arras's avatar
Philipp Arras committed
662
            raise NotImplementedError
663
        if space >= len(self._a):
664
            raise ValueError("invalid space specified; got {!r}".format(space))
665
        if len(self._a) == 1:
666
            return self.average_fluctuation(0)
667
668
        q = 1.
        for j in range(len(self._a)):
Martin Reinecke's avatar
Martin Reinecke committed
669
            fl = self._a[j].fluctuation_amplitude*self._azm.ptw("reciprocal")
670
            if j == space:
Philipp Arras's avatar
Philipp Arras committed
671
                q = q*fl**2
672
            else:
Philipp Arras's avatar
Philipp Arras committed
673
                q = q*(Adder(full(fl.target, 1.)) @ fl**2)
Martin Reinecke's avatar
Martin Reinecke committed
674
        return q.ptw("sqrt")*self._azm
Philipp Arras's avatar
Philipp Arras committed
675
676

    def average_fluctuation(self, space):
677
        """Returns operator which acts on prior or posterior samples"""
678
        if len(self._a) == 0:
Philipp Arras's avatar
Philipp Arras committed
679
            raise NotImplementedError
680
        if space >= len(self._a):
681
            raise ValueError("invalid space specified; got {!r}".format(space))
682
        if len(self._a) == 1:
Philipp Haim's avatar
Philipp Haim committed
683
684
            return self._a[0].fluctuation_amplitude
        return self._a[space].fluctuation_amplitude
685

686
687
    @staticmethod
    def offset_amplitude_realized(samples):
Philipp Haim's avatar
Philipp Haim committed
688
        spaces = _structured_spaces(samples[0].domain)
689
690
        res = 0.
        for s in samples:
Philipp Haim's avatar
Philipp Haim committed
691
            res = res + s.mean(spaces)**2
Philipp Haim's avatar
Philipp Haim committed
692
693
        res = res/len(samples)
        return np.sqrt(res if np.isscalar(res) else res.val)
Philipp Arras's avatar
Philipp Arras committed
694

695
696
697
698
699
700
701
702
    @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
703
704
        spaces = _structured_spaces(samples[0].domain)
        if space >= len(spaces):
705
            raise ValueError("invalid space specified; got {!r}".format(space))
Philipp Haim's avatar
Philipp Haim committed
706
        if len(spaces) == 1:
707
            return _total_fluctuation_realized(samples)
Philipp Haim's avatar
Philipp Haim committed
708
        space = space + spaces[0]
Philipp Arras's avatar
Philipp Arras committed
709
        res1, res2 = 0., 0.
710
        for s in samples:
Philipp Frank's avatar
fixes    
Philipp Frank committed
711
712
713
714
            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
715
        res = res1.mean(spaces) - res2.mean(spaces[:-1])
Philipp Haim's avatar
Philipp Haim committed
716
        return np.sqrt(res if np.isscalar(res) else res.val)
Philipp Frank's avatar
fixes    
Philipp Frank committed
717

Philipp Arras's avatar
Philipp Arras committed
718
    @staticmethod
719
720
721
    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
722
723
        spaces = _structured_spaces(samples[0].domain)
        if space >= len(spaces):
724
            raise ValueError("invalid space specified; got {!r}".format(space))
Philipp Haim's avatar
Philipp Haim committed
725
        if len(spaces) == 1:
726
            return _total_fluctuation_realized(samples)
Philipp Haim's avatar
Philipp Haim committed
727
728
729
        space = space + spaces[0]
        sub_spaces = set(spaces)
        sub_spaces.remove(space)
Philipp Arras's avatar
Philipp Arras committed
730
        # Domain containing domain[space] and domain[0] iff total_N>0
Philipp Haim's avatar
Philipp Haim committed
731
        sub_dom = makeDomain([samples[0].domain[ind]
Philipp Arras's avatar
Philipp Arras committed
732
                              for ind in (set([0])-set(spaces)) | set([space])])
Philipp Haim's avatar
Philipp Haim committed
733
        co = ContractionOperator(sub_dom, len(sub_dom)-1)
734
        size = co.domain.size/co.target.size
Philipp Arras's avatar
Philipp Arras committed
735
736
        res = 0.
        for s in samples:
Philipp Haim's avatar
Philipp Haim committed
737
            r = s.mean(sub_spaces)
738
            res = res + (r - co.adjoint(co(r)/size))**2
Philipp Haim's avatar
Philipp Haim committed
739
        res = res.mean(spaces[0])/len(samples)
Philipp Haim's avatar
Philipp Haim committed
740
        return np.sqrt(res if np.isscalar(res) else res.val)