correlated_fields.py 30.4 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
2
3
4
5
6
7
8
9
10
11
12
13
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
14
# Copyright(C) 2013-2020 Max-Planck-Society
Philipp Arras's avatar
Philipp Arras committed
15
# Authors: Philipp Frank, Philipp Arras, Philipp Haim
Martin Reinecke's avatar
Martin Reinecke committed
16
#
17
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
Philipp Arras's avatar
Philipp Arras committed
18

19
20
21
from functools import reduce
from operator import mul

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

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

43

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

53

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

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


Martin Reinecke's avatar
Martin Reinecke committed
69
def _normal(mean, sig, key, N=0):
Philipp Haim's avatar
Philipp Haim committed
70
    if N == 0:
Philipp Haim's avatar
Philipp Haim committed
71
        domain = DomainTuple.scalar_domain()
Philipp Haim's avatar
Philipp Haim committed
72
        mean, sig = np.asfarray(mean), np.asfarray(sig)
73
74
75
76
77
        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
78
79
    return Adder(makeField(domain, mean)) @ (DiagonalOperator(
        makeField(domain, sig)) @ ducktape(domain, None, key))
Philipp Arras's avatar
Philipp Arras committed
80
81


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


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


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


Philipp Arras's avatar
Philipp Arras committed
123
class _LognormalMomentMatching(Operator):
Philipp Haim's avatar
Philipp Haim committed
124
    def __init__(self, mean, sig, key, N_copies):
Philipp Arras's avatar
Philipp Arras committed
125
        key = str(key)
Philipp Haim's avatar
Philipp Haim committed
126
        logmean, logsig = _lognormal_moments(mean, sig, N_copies)
Philipp Arras's avatar
Philipp Arras committed
127
128
        self._mean = mean
        self._sig = sig
Martin Reinecke's avatar
Martin Reinecke committed
129
        op = _normal(logmean, logsig, key, N_copies).ptw("exp")
Philipp Arras's avatar
Philipp Arras committed
130
131
        self._domain, self._target = op.domain, op.target
        self.apply = op.apply
132
        self._repr_str = f"_LognormalMomentMatching: " + op.__repr__()
Philipp Arras's avatar
Philipp Arras committed
133
134
135
136
137
138
139
140

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

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

142
143
144
    def __repr__(self):
        return self._repr_str

Philipp Arras's avatar
Philipp Arras committed
145

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

153
        self._space = space
Philipp Haim's avatar
Philipp Haim committed
154
155
156
        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
157
        self._capability = self.TIMES | self.ADJOINT_TIMES
Philipp Arras's avatar
Philipp Arras committed
158

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

Philipp Arras's avatar
Philipp Arras committed
170
171

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

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

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

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


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

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


class _SpecialSum(EndomorphicOperator):
Martin Reinecke's avatar
Martin Reinecke committed
239
    def __init__(self, domain, space=0):
Philipp Arras's avatar
Philipp Arras committed
240
241
        self._domain = makeDomain(domain)
        self._capability = self.TIMES | self.ADJOINT_TIMES
242
        self._contractor = ContractionOperator(domain, space)
Philipp Arras's avatar
Philipp Arras committed
243
244
245

    def apply(self, x, mode):
        self._check_input(x, mode)
246
        return self._contractor.adjoint(self._contractor(x))
Philipp Arras's avatar
Philipp Arras committed
247
248


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

266

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

295
        twolog = _TwoLogIntegrations(target, space)
Philipp Arras's avatar
Philipp Arras committed
296
        dom = twolog.domain
297

298
        shp = dom[space].shape
Martin Reinecke's avatar
Martin Reinecke committed
299
300
        expander = ContractionOperator(dom, spaces=space).adjoint
        ps_expander = ContractionOperator(twolog.target, spaces=space).adjoint
Philipp Arras's avatar
Philipp Arras committed
301
302
303

        # Prepare constant fields
        foo = np.zeros(shp)
304
        foo[0] = foo[1] = np.sqrt(_log_vol(target[space]))
Martin Reinecke's avatar
Martin Reinecke committed
305
        vflex = DiagonalOperator(makeField(dom[space], foo), dom, space)
Philipp Arras's avatar
Philipp Arras committed
306
307
308

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

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

315
        vslope = DiagonalOperator(
Philipp Arras's avatar
Philipp Arras committed
316
            makeField(target[space], _relative_log_k_lengths(target[space])),
Martin Reinecke's avatar
Martin Reinecke committed
317
            target, space)
318
319

        foo, bar = [np.zeros(target[space].shape) for _ in range(2)]
Philipp Arras's avatar
Philipp Arras committed
320
        bar[1:] = foo[0] = totvol
Philipp Arras's avatar
Philipp Arras committed
321
322
323
324
        vol0, vol1 = [
            DiagonalOperator(makeField(target[space], aa), target, space)
            for aa in (foo, bar)
        ]
325

Martin Reinecke's avatar
Martin Reinecke committed
326
        # Prepare fields for Adder
327
        shift, vol0 = [op(full(op.domain, 1)) for op in (shift, vol0)]
Philipp Arras's avatar
Philipp Arras committed
328
329
        # End prepare constant fields

330
331
332
333
        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
334
        sig_fluc = vol1 @ ps_expander @ fluctuations
Philipp Arras's avatar
Philipp Arras committed
335
336

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

Philipp Arras's avatar
Philipp Arras committed
350
351
        self.apply = op.apply
        self._domain, self._target = op.domain, op.target
352
        self._space = space
353
        self._repr_str = "_Amplitude: " + op.__repr__()
Philipp Arras's avatar
Philipp Arras committed
354

Philipp Arras's avatar
Philipp Arras committed
355
356
357
358
    @property
    def fluctuation_amplitude(self):
        return self._fluc

359
360
361
    def __repr__(self):
        return self._repr_str

362
363

class CorrelatedFieldMaker:
Lukas Platz's avatar
Lukas Platz committed
364
365
366
    """Constrution helper for hirarchical correlated field models.

    The correlated field models are parametrized by creating
367
368
    power spectrum operators ("amplitudes") via calls to
    :func:`add_fluctuations` that act on the targeted field subdomains.
Lukas Platz's avatar
Lukas Platz committed
369
    During creation of the :class:`CorrelatedFieldMaker` via
370
371
372
    :func:`make`, a global offset from zero of the field model
    can be defined and an operator applying fluctuations
    around this offset is parametrized.
Lukas Platz's avatar
Lukas Platz committed
373
374

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

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

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

386
387
388
389
390
391
392
393
394
395
    An operator representing an array of correlated field models
    can be constructed by setting the `total_N` parameter of
    :func:`make`. It will have an
    :class:`~nifty.domains.unstructucture_domain.UnstructureDomain`
    of shape `(total_N,)` prepended to its target domain and represent
    `total_N` correlated fields simulataneously.
    The degree of information sharing between the correlated field
    models can be configured via the `dofdex` parameters
    of :func:`make` and :func:`add_fluctuations`.

Lukas Platz's avatar
Lukas Platz committed
396
    See the methods :func:`make`, :func:`add_fluctuations`
397
    and :func:`finalize` for further usage information."""
398
399
400
    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")
401
        self._a = []
402
        self._target_subdomains = []
Philipp Arras's avatar
Formats    
Philipp Arras committed
403

404
405
        self._offset_mean = offset_mean
        self._azm = offset_fluctuations_op
406
        self._prefix = prefix
Philipp Haim's avatar
Philipp Haim committed
407
        self._total_N = total_N
Philipp Arras's avatar
Formats    
Philipp Arras committed
408

409
    @staticmethod
Philipp Arras's avatar
Philipp Arras committed
410
    def make(offset_mean, offset_std_mean, offset_std_std, prefix, total_N=0,
Martin Reinecke's avatar
Martin Reinecke committed
411
             dofdex=None):
Lukas Platz's avatar
Lukas Platz committed
412
413
414
415
416
417
418
        """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
419
            Mean standard deviation of the offset.
Lukas Platz's avatar
Lukas Platz committed
420
        offset_std_std : float
Lukas Platz's avatar
Lukas Platz committed
421
            Standard deviation of the stddev of the offset.
Lukas Platz's avatar
Lukas Platz committed
422
423
        prefix : string
            Prefix to the names of the domains of the cf operator to be made.
Lukas Platz's avatar
Lukas Platz committed
424
            This determines the names of the operator domain.
425
426
        total_N : integer, optional
            Number of field models to create.
Lukas Platz's avatar
Lukas Platz committed
427
428
429
            If not 0, the first entry of the operators target will be an
            :class:`~nifty.domains.unstructured_domain.UnstructuredDomain`
            with length `total_N`.
430
        dofdex : np.array of integers, optional
Philipp Arras's avatar
Philipp Arras committed
431
432
433
            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
434
435
436
437
            instanciated; the first one is used for the first and second
            field model and the second is used for the third field model.
            *If not specified*, use the same zero mode model for all
            constructed field models.
Lukas Platz's avatar
Lukas Platz committed
438
        """
Philipp Frank's avatar
Philipp Frank committed
439
440
        if dofdex is None:
            dofdex = np.full(total_N, 0)
441
442
        elif len(dofdex) != total_N:
            raise ValueError("length of dofdex needs to match total_N")
Philipp Frank's avatar
Philipp Frank committed
443
        N = max(dofdex) + 1 if total_N > 0 else 0
Lukas Platz's avatar
Lukas Platz committed
444
445
        zm = _LognormalMomentMatching(offset_std_mean,
                                      offset_std_std,
Philipp Haim's avatar
Philipp Haim committed
446
                                      prefix + 'zeromode',
Philipp Frank's avatar
Philipp Frank committed
447
                                      N)
Philipp Frank's avatar
fixup    
Philipp Frank committed
448
        if total_N > 0:
Martin Reinecke's avatar
Martin Reinecke committed
449
            zm = _Distributor(dofdex, zm.target, UnstructuredDomain(total_N)) @ zm
450
        return CorrelatedFieldMaker(offset_mean, zm, prefix, total_N)
451
452

    def add_fluctuations(self,
453
                         target_subdomain,
454
455
456
457
458
459
460
461
                         fluctuations_mean,
                         fluctuations_stddev,
                         flexibility_mean,
                         flexibility_stddev,
                         asperity_mean,
                         asperity_stddev,
                         loglogavgslope_mean,
                         loglogavgslope_stddev,
Martin Reinecke's avatar
Martin Reinecke committed
462
463
464
465
                         prefix='',
                         index=None,
                         dofdex=None,
                         harmonic_partner=None):
Lukas Platz's avatar
Lukas Platz committed
466
467
468
469
470
471
        """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
472
473
        `loglogavgslope` configure the power spectrum model ("amplitude")
        used on the target field subdomain `target_subdomain`.
Lukas Platz's avatar
Lukas Platz committed
474
475
476
477
478
479
480
481
482
483
484
        It is assembled as the sum of a power law component
        (linear slope in log-log power-frequency-space),
        a smooth varying component (integrated wiener process) and
        a ragged componenent (unintegrated wiener process).

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

        Parameters
        ----------
Martin Reinecke's avatar
Martin Reinecke committed
485
486
        target_subdomain : :class:`~nifty7.domain.Domain`, \
                           :class:`~nifty7.domain_tuple.DomainTuple`
Lukas Platz's avatar
Lukas Platz committed
487
488
489
490
491
            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
492
            Amplitude of the non-power-law power spectrum component
Lukas Platz's avatar
Lukas Platz committed
493
        asperity_{mean,stddev} : float
494
            Roughness of the non-power-law power spectrum component
Lukas Platz's avatar
Lukas Platz committed
495
496
497
498
499
            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
500
501
502
        index : int
            Position target_subdomain in the final total domain of the
            correlated field operator.
503
504
        dofdex : np.array, optional
            An integer array specifying the power spectrum models used if
Philipp Arras's avatar
Philipp Arras committed
505
            total_N > 1. It needs to have length of total_N. If total_N=3 and
506
507
508
509
510
            dofdex=[0,0,1], that means that two power spectrum models are
            instanciated; the first one is used for the first and second
            field model and the second one is used for the third field model.
            *If not given*, use the same power spectrum model for all
            constructed field models.
Martin Reinecke's avatar
Martin Reinecke committed
511
512
        harmonic_partner : :class:`~nifty7.domain.Domain`, \
                           :class:`~nifty7.domain_tuple.DomainTuple`
Lukas Platz's avatar
Lukas Platz committed
513
514
            In which harmonic space to define the power spectrum
        """
Philipp Frank's avatar
Philipp Frank committed
515
        if harmonic_partner is None:
516
            harmonic_partner = target_subdomain.get_default_codomain()
Philipp Frank's avatar
Fixup    
Philipp Frank committed
517
        else:
518
519
            target_subdomain.check_codomain(harmonic_partner)
            harmonic_partner.check_codomain(target_subdomain)
520

Philipp Haim's avatar
Philipp Haim committed
521
522
        if dofdex is None:
            dofdex = np.full(self._total_N, 0)
523
524
        elif len(dofdex) != self._total_N:
            raise ValueError("length of dofdex needs to match total_N")
Philipp Haim's avatar
Philipp Haim committed
525

Philipp Haim's avatar
Philipp Haim committed
526
527
        if self._total_N > 0:
            N = max(dofdex) + 1
528
            target_subdomain = makeDomain((UnstructuredDomain(N), target_subdomain))
Philipp Haim's avatar
Philipp Haim committed
529
        else:
Philipp Haim's avatar
Philipp Haim committed
530
            N = 0
531
            target_subdomain = makeDomain(target_subdomain)
Philipp Arras's avatar
Philipp Arras committed
532
        prefix = str(prefix)
533
        # assert isinstance(target_subdomain[space], (RGSpace, HPSpace, GLSpace)
Philipp Arras's avatar
Philipp Arras committed
534

Philipp Arras's avatar
Philipp Arras committed
535
536
        fluct = _LognormalMomentMatching(fluctuations_mean,
                                         fluctuations_stddev,
537
                                         self._prefix + prefix + 'fluctuations',
Philipp Haim's avatar
Philipp Haim committed
538
                                         N)
Philipp Arras's avatar
Philipp Arras committed
539
        flex = _LognormalMomentMatching(flexibility_mean, flexibility_stddev,
540
                                        self._prefix + prefix + 'flexibility',
Philipp Haim's avatar
Philipp Haim committed
541
                                        N)
Philipp Arras's avatar
Philipp Arras committed
542
        asp = _LognormalMomentMatching(asperity_mean, asperity_stddev,
Philipp Arras's avatar
Philipp Arras committed
543
                                       self._prefix + prefix + 'asperity', N)
544
        avgsl = _normal(loglogavgslope_mean, loglogavgslope_stddev,
545
                        self._prefix + prefix + 'loglogavgslope', N)
Philipp Arras's avatar
Philipp Arras committed
546
        amp = _Amplitude(PowerSpace(harmonic_partner), fluct, flex, asp, avgsl,
547
                         self._azm, target_subdomain[-1].total_volume,
548
                         self._prefix + prefix + 'spectrum', dofdex)
Philipp Haim's avatar
Philipp Haim committed
549

550
551
        if index is not None:
            self._a.insert(index, amp)
552
            self._target_subdomains.insert(index, target_subdomain)
553
554
        else:
            self._a.append(amp)
555
            self._target_subdomains.append(target_subdomain)
556

Philipp Arras's avatar
Philipp Arras committed
557
558
559
560
561
562
563
564
565
566
    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
567
        n_amplitudes = len(self._a)
Philipp Haim's avatar
Philipp Haim committed
568
        if self._total_N > 0:
Philipp Arras's avatar
Philipp Arras committed
569
570
571
            hspace = makeDomain(
                [UnstructuredDomain(self._total_N)] +
                [dd.target[-1].harmonic_partner for dd in self._a])
Philipp Haim's avatar
Philipp Haim committed
572
573
            spaces = tuple(range(1, n_amplitudes + 1))
            amp_space = 1
Philipp Haim's avatar
Philipp Haim committed
574
575
        else:
            hspace = makeDomain(
Philipp Arras's avatar
Philipp Arras committed
576
                [dd.target[0].harmonic_partner for dd in self._a])
Philipp Haim's avatar
Philipp Haim committed
577
            spaces = tuple(range(n_amplitudes))
Philipp Haim's avatar
Philipp Haim committed
578
            amp_space = 0
579

Martin Reinecke's avatar
Martin Reinecke committed
580
        expander = ContractionOperator(hspace, spaces=spaces).adjoint
Philipp Frank's avatar
fixup    
Philipp Frank committed
581
        azm = expander @ self._azm
582

583
        ht = HarmonicTransformOperator(hspace,
584
                                       self._target_subdomains[0][amp_space],
Martin Reinecke's avatar
Martin Reinecke committed
585
                                       space=spaces[0])
586
        for i in range(1, n_amplitudes):
587
            ht = HarmonicTransformOperator(ht.target,
588
                                           self._target_subdomains[i][amp_space],
589
590
591
592
593
                                           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
594
            pd = PowerDistributor(co.target, pp, amp_space)
595
596
            a.append(co.adjoint @ pd @ self._a[ii])
        corr = reduce(mul, a)
Philipp Arras's avatar
Philipp Arras committed
597
        op = ht(azm*corr*ducktape(hspace, None, self._prefix + 'xi'))
Philipp Arras's avatar
Philipp Arras committed
598

599
600
        if self._offset_mean is not None:
            offset = self._offset_mean
601
602
603
604
605
606
607
            # 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
608
        self.statistics_summary(prior_info)
609
610
        return op

611
612
613
614
615
616
    def statistics_summary(self, prior_info):
        from ..sugar import from_random

        if prior_info == 0:
            return

617
618
        lst = [('Offset amplitude', self.amplitude_total_offset),
               ('Total fluctuation amplitude', self.total_fluctuation)]
619
        namps = len(self._a)
620
621
622
623
624
625
626
627
        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:
628
629
            sc = StatCalculator()
            for _ in range(prior_info):
630
                sc.add(op(from_random(op.domain, 'normal')))
Martin Reinecke's avatar
merge    
Martin Reinecke committed
631
            mean = sc.mean.val
Martin Reinecke's avatar
Martin Reinecke committed
632
            stddev = sc.var.ptw("sqrt").val
633
            for m, s in zip(mean.flatten(), stddev.flatten()):
634
                logger.info('{}: {:.02E} ± {:.02E}'.format(kk, m, s))
635
636
637

    def moment_slice_to_average(self, fluctuations_slice_mean, nsamples=1000):
        fluctuations_slice_mean = float(fluctuations_slice_mean)
638
639
640
        if not fluctuations_slice_mean > 0:
            msg = "fluctuations_slice_mean must be greater zero; got {!r}"
            raise ValueError(msg.format(fluctuations_slice_mean))
641
        from ..sugar import from_random
642
643
        scm = 1.
        for a in self._a:
Martin Reinecke's avatar
Martin Reinecke committed
644
            op = a.fluctuation_amplitude*self._azm.ptw("reciprocal")
645
            res = np.array([op(from_random(op.domain, 'normal')).val
646
647
                            for _ in range(nsamples)])
            scm *= res**2 + 1.
648
        return fluctuations_slice_mean/np.mean(np.sqrt(scm))
649

Philipp Arras's avatar
Philipp Arras committed
650
    @property
Philipp Haim's avatar
Philipp Haim committed
651
    def normalized_amplitudes(self):
652
        """Returns the power spectrum operators used in the model"""
653
        return self._a
Philipp Arras's avatar
Philipp Arras committed
654

Philipp Haim's avatar
Philipp Haim committed
655
656
657
658
659
660
661
    @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
662
663
        dom = self._a[0].target
        expand = ContractionOperator(dom, len(dom)-1).adjoint
Philipp Haim's avatar
Philipp Haim committed
664
665
        return self._a[0]*(expand @ self.amplitude_total_offset)

666
667
668
    @property
    def amplitude_total_offset(self):
        return self._azm
Philipp Arras's avatar
Philipp Arras committed
669
670

    @property
671
    def total_fluctuation(self):
672
        """Returns operator which acts on prior or posterior samples"""
673
        if len(self._a) == 0:
Philipp Arras's avatar
Philipp Arras committed
674
            raise NotImplementedError
675
        if len(self._a) == 1:
676
            return self.average_fluctuation(0)
677
678
        q = 1.
        for a in self._a:
Martin Reinecke's avatar
Martin Reinecke committed
679
            fl = a.fluctuation_amplitude*self._azm.ptw("reciprocal")
Philipp Arras's avatar
Philipp Arras committed
680
            q = q*(Adder(full(fl.target, 1.)) @ fl**2)
Martin Reinecke's avatar
Martin Reinecke committed
681
        return (Adder(full(q.target, -1.)) @ q).ptw("sqrt")*self._azm
682

Philipp Arras's avatar
Philipp Arras committed
683
    def slice_fluctuation(self, space):
684
        """Returns operator which acts on prior or posterior samples"""
685
        if len(self._a) == 0:
Philipp Arras's avatar
Philipp Arras committed
686
            raise NotImplementedError
687
        if space >= len(self._a):
688
            raise ValueError("invalid space specified; got {!r}".format(space))
689
        if len(self._a) == 1:
690
            return self.average_fluctuation(0)
691
692
        q = 1.
        for j in range(len(self._a)):
Martin Reinecke's avatar
Martin Reinecke committed
693
            fl = self._a[j].fluctuation_amplitude*self._azm.ptw("reciprocal")
694
            if j == space:
Philipp Arras's avatar
Philipp Arras committed
695
                q = q*fl**2
696
            else:
Philipp Arras's avatar
Philipp Arras committed
697
                q = q*(Adder(full(fl.target, 1.)) @ fl**2)
Martin Reinecke's avatar
Martin Reinecke committed
698
        return q.ptw("sqrt")*self._azm
Philipp Arras's avatar
Philipp Arras committed
699
700

    def average_fluctuation(self, space):
701
        """Returns operator which acts on prior or posterior samples"""
702
        if len(self._a) == 0:
Philipp Arras's avatar
Philipp Arras committed
703
            raise NotImplementedError
704
        if space >= len(self._a):
705
            raise ValueError("invalid space specified; got {!r}".format(space))
706
        if len(self._a) == 1:
Philipp Haim's avatar
Philipp Haim committed
707
708
            return self._a[0].fluctuation_amplitude
        return self._a[space].fluctuation_amplitude
709

710
711
    @staticmethod
    def offset_amplitude_realized(samples):
Philipp Haim's avatar
Philipp Haim committed
712
        spaces = _structured_spaces(samples[0].domain)
713
714
        res = 0.
        for s in samples:
Philipp Haim's avatar
Philipp Haim committed
715
            res = res + s.mean(spaces)**2
Philipp Haim's avatar
Philipp Haim committed
716
717
        res = res/len(samples)
        return np.sqrt(res if np.isscalar(res) else res.val)
Philipp Arras's avatar
Philipp Arras committed
718

719
720
721
722
723
724
725
726
    @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
727
728
        spaces = _structured_spaces(samples[0].domain)
        if space >= len(spaces):
729
            raise ValueError("invalid space specified; got {!r}".format(space))
Philipp Haim's avatar
Philipp Haim committed
730
        if len(spaces) == 1:
731
            return _total_fluctuation_realized(samples)
Philipp Haim's avatar
Philipp Haim committed
732
        space = space + spaces[0]
Philipp Arras's avatar
Philipp Arras committed
733
        res1, res2 = 0., 0.
734
        for s in samples:
Philipp Frank's avatar
fixes    
Philipp Frank committed
735
736
737
738
            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
739
        res = res1.mean(spaces) - res2.mean(spaces[:-1])
Philipp Haim's avatar
Philipp Haim committed
740
        return np.sqrt(res if np.isscalar(res) else res.val)
Philipp Frank's avatar
fixes    
Philipp Frank committed
741

Philipp Arras's avatar
Philipp Arras committed
742
    @staticmethod
743
744
745
    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
746
747
        spaces = _structured_spaces(samples[0].domain)
        if space >= len(spaces):
748
            raise ValueError("invalid space specified; got {!r}".format(space))
Philipp Haim's avatar
Philipp Haim committed
749
        if len(spaces) == 1:
750
            return _total_fluctuation_realized(samples)
Philipp Haim's avatar
Philipp Haim committed
751
752
753
        space = space + spaces[0]
        sub_spaces = set(spaces)
        sub_spaces.remove(space)
Philipp Arras's avatar
Philipp Arras committed
754
        # Domain containing domain[space] and domain[0] iff total_N>0
Philipp Haim's avatar
Philipp Haim committed
755
        sub_dom = makeDomain([samples[0].domain[ind]
Philipp Arras's avatar
Philipp Arras committed
756
                              for ind in (set([0])-set(spaces)) | set([space])])
Philipp Haim's avatar
Philipp Haim committed
757
        co = ContractionOperator(sub_dom, len(sub_dom)-1)
758
        size = co.domain.size/co.target.size
Philipp Arras's avatar
Philipp Arras committed
759
760
        res = 0.
        for s in samples:
Philipp Haim's avatar
Philipp Haim committed
761
            r = s.mean(sub_spaces)
762
            res = res + (r - co.adjoint(co(r)/size))**2
Philipp Haim's avatar
Philipp Haim committed
763
        res = res.mean(spaces[0])/len(samples)
Philipp Haim's avatar
Philipp Haim committed
764
        return np.sqrt(res if np.isscalar(res) else res.val)