correlated_fields.py 29.3 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1 2 3 4 5 6 7 8 9 10 11 12 13
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
14
# Copyright(C) 2013-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 Frank's avatar
fixup  
Philipp Frank committed
534
    def _finalize_from_op(self):
Philipp Haim's avatar
Philipp Haim committed
535
        n_amplitudes = len(self._a)
Philipp Haim's avatar
Philipp Haim committed
536
        if self._total_N > 0:
Philipp Arras's avatar
Philipp Arras committed
537 538 539
            hspace = makeDomain(
                [UnstructuredDomain(self._total_N)] +
                [dd.target[-1].harmonic_partner for dd in self._a])
Philipp Haim's avatar
Philipp Haim committed
540 541
            spaces = tuple(range(1, n_amplitudes + 1))
            amp_space = 1
Philipp Haim's avatar
Philipp Haim committed
542 543
        else:
            hspace = makeDomain(
Philipp Arras's avatar
Philipp Arras committed
544
                [dd.target[0].harmonic_partner for dd in self._a])
Philipp Haim's avatar
Philipp Haim committed
545
            spaces = tuple(range(n_amplitudes))
Philipp Haim's avatar
Philipp Haim committed
546
            amp_space = 0
547

Martin Reinecke's avatar
Martin Reinecke committed
548
        expander = ContractionOperator(hspace, spaces=spaces).adjoint
Philipp Frank's avatar
fixup  
Philipp Frank committed
549
        azm = expander @ self._azm
550

551
        ht = HarmonicTransformOperator(hspace,
552
                                       self._target_subdomains[0][amp_space],
Martin Reinecke's avatar
Martin Reinecke committed
553
                                       space=spaces[0])
554
        for i in range(1, n_amplitudes):
555
            ht = HarmonicTransformOperator(ht.target,
556
                                           self._target_subdomains[i][amp_space],
557 558 559 560 561
                                           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
562
            pd = PowerDistributor(co.target, pp, amp_space)
563 564 565
            a.append(co.adjoint @ pd @ self._a[ii])
        corr = reduce(mul, a)
        return ht(azm*corr*ducktape(hspace, None, self._prefix + 'xi'))
Philipp Arras's avatar
Philipp Arras committed
566

567
    def finalize(self, prior_info=100):
Lukas Platz's avatar
Lukas Platz committed
568 569 570 571 572 573 574 575
        """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 Arras's avatar
Philipp Arras committed
576
        """
Philipp Frank's avatar
fixup  
Philipp Frank committed
577
        op = self._finalize_from_op()
578 579
        if self._offset_mean is not None:
            offset = self._offset_mean
580 581 582 583 584 585 586
            # 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
587
        self.statistics_summary(prior_info)
588 589
        return op

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

        if prior_info == 0:
            return

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

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

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

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

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

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

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

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

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

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

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