correlated_fields.py 30.1 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 ..operators.normal_operators import NormalTransform, LognormalTransform
41
from ..probing import StatCalculator
Philipp Arras's avatar
Philipp Arras committed
42
from ..sugar import full, makeDomain, makeField, makeOp
43

44

Philipp Arras's avatar
Philipp Arras committed
45
def _log_k_lengths(pspace):
Philipp Arras's avatar
Philipp Arras committed
46
    """Log(k_lengths) without zeromode"""
Philipp Arras's avatar
Philipp Arras committed
47 48 49
    return np.log(pspace.k_lengths[1:])


Philipp Arras's avatar
Philipp Arras committed
50
def _relative_log_k_lengths(power_space):
Philipp Arras's avatar
Philipp Arras committed
51 52
    """Log-distance to first bin
    logkl.shape==power_space.shape, logkl[0]=logkl[1]=0"""
Philipp Arras's avatar
Philipp Arras committed
53 54 55 56 57 58
    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
59
    return np.insert(logkl, 0, 0)
Philipp Arras's avatar
Philipp Arras committed
60 61


Philipp Arras's avatar
Philipp Arras committed
62
def _log_vol(power_space):
63
    power_space = makeDomain(power_space)
Philipp Arras's avatar
Philipp Arras committed
64 65 66 67 68
    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
69 70 71 72 73 74
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
75
def _total_fluctuation_realized(samples):
Philipp Haim's avatar
Philipp Haim committed
76 77 78
    spaces = _structured_spaces(samples[0].domain)
    co = ContractionOperator(samples[0].domain, spaces)
    size = co.domain.size/co.target.size
79 80
    res = 0.
    for s in samples:
Philipp Haim's avatar
Philipp Haim committed
81 82
        res = res + (s - co.adjoint(co(s)/size))**2
    res = res.mean(spaces)/len(samples)
Philipp Haim's avatar
Philipp Haim committed
83
    return np.sqrt(res if np.isscalar(res) else res.val)
84 85


Philipp Frank's avatar
Philipp Frank committed
86
class _SlopeRemover(EndomorphicOperator):
Martin Reinecke's avatar
Martin Reinecke committed
87
    def __init__(self, domain, space=0):
Philipp Frank's avatar
Philipp Frank committed
88
        self._domain = makeDomain(domain)
89 90
        assert isinstance(self._domain[space], PowerSpace)
        logkl = _relative_log_k_lengths(self._domain[space])
91
        self._sc = logkl/float(logkl[-1])
Philipp Arras's avatar
Philipp Arras committed
92

93
        self._space = space
Philipp Haim's avatar
Philipp Haim committed
94 95 96
        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
97
        self._capability = self.TIMES | self.ADJOINT_TIMES
Philipp Arras's avatar
Philipp Arras committed
98

99 100
    def apply(self, x, mode):
        self._check_input(x, mode)
Martin Reinecke's avatar
stage 3  
Martin Reinecke committed
101
        x = x.val
Philipp Frank's avatar
Philipp Frank committed
102
        if mode == self.TIMES:
Philipp Haim's avatar
Philipp Haim committed
103
            res = x - x[self._last]*self._sc[self._extender]
Philipp Frank's avatar
Philipp Frank committed
104
        else:
105 106
            res = x.copy()
            res[self._last] -= (x*self._sc[self._extender]).sum(
Philipp Arras's avatar
Philipp Arras committed
107
                axis=self._space, keepdims=True)
Martin Reinecke's avatar
Martin Reinecke committed
108
        return makeField(self._tgt(mode), res)
Philipp Frank's avatar
Philipp Frank committed
109

Philipp Arras's avatar
Philipp Arras committed
110 111

class _TwoLogIntegrations(LinearOperator):
Martin Reinecke's avatar
Martin Reinecke committed
112
    def __init__(self, target, space=0):
Philipp Arras's avatar
Philipp Arras committed
113
        self._target = makeDomain(target)
114 115 116 117 118
        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
119
        self._log_vol = _log_vol(self._target[space])
Philipp Arras's avatar
Philipp Arras committed
120 121 122 123
        self._capability = self.TIMES | self.ADJOINT_TIMES

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

Martin Reinecke's avatar
Martin Reinecke committed
125
        # Maybe make class properties
126 127
        axis = self._target.axes[self._space][0]
        sl = (slice(None),)*axis
Philipp Haim's avatar
Fixes  
Philipp Haim committed
128
        extender_sl = (None,)*axis + (slice(None),) + (None,)*(self._target.axes[-1][-1] - axis)
129 130
        first = sl + (0,)
        second = sl + (1,)
Martin Reinecke's avatar
Martin Reinecke committed
131 132 133
        from_third = sl + (slice(2, None),)
        no_border = sl + (slice(1, -1),)
        reverse = sl + (slice(None, None, -1),)
134

Philipp Arras's avatar
Philipp Arras committed
135
        if mode == self.TIMES:
Martin Reinecke's avatar
stage 3  
Martin Reinecke committed
136
            x = x.val
Philipp Arras's avatar
Philipp Arras committed
137
            res = np.empty(self._target.shape)
138
            res[first] = res[second] = 0
Martin Reinecke's avatar
Martin Reinecke committed
139
            res[from_third] = np.cumsum(x[second], axis=axis)
140
            res[from_third] = (res[from_third] + res[no_border])/2*self._log_vol[extender_sl] + x[first]
Martin Reinecke's avatar
Martin Reinecke committed
141
            res[from_third] = np.cumsum(res[from_third], axis=axis)
Philipp Arras's avatar
Philipp Arras committed
142
        else:
Martin Reinecke's avatar
Martin Reinecke committed
143
            x = x.val_rw()
Philipp Arras's avatar
Philipp Arras committed
144
            res = np.zeros(self._domain.shape)
Martin Reinecke's avatar
Martin Reinecke committed
145
            x[from_third] = np.cumsum(x[from_third][reverse], axis=axis)[reverse]
146
            res[first] += x[from_third]
147
            x[from_third] *= (self._log_vol/2.)[extender_sl]
148
            x[no_border] += x[from_third]
Martin Reinecke's avatar
Martin Reinecke committed
149
            res[second] += np.cumsum(x[from_third][reverse], axis=axis)[reverse]
Martin Reinecke's avatar
Martin Reinecke committed
150
        return makeField(self._tgt(mode), res)
Philipp Arras's avatar
Philipp Arras committed
151 152 153


class _Normalization(Operator):
Martin Reinecke's avatar
Martin Reinecke committed
154
    def __init__(self, domain, space=0):
Philipp Arras's avatar
Philipp Arras committed
155
        self._domain = self._target = DomainTuple.make(domain)
156
        assert isinstance(self._domain[space], PowerSpace)
157 158 159
        hspace = list(self._domain)
        hspace[space] = hspace[space].harmonic_partner
        hspace = makeDomain(hspace)
Philipp Arras's avatar
Philipp Arras committed
160 161 162
        pd = PowerDistributor(hspace,
                              power_space=self._domain[space],
                              space=space)
Martin Reinecke's avatar
Martin Reinecke committed
163
        mode_multiplicity = pd.adjoint(full(pd.target, 1.)).val_rw()
164
        zero_mode = (slice(None),)*self._domain.axes[space][0] + (0,)
Philipp Haim's avatar
Philipp Haim committed
165
        mode_multiplicity[zero_mode] = 0
Philipp Arras's avatar
Philipp Arras committed
166
        self._mode_multiplicity = makeOp(makeField(self._domain, mode_multiplicity))
167
        self._specsum = _SpecialSum(self._domain, space)
Philipp Arras's avatar
Philipp Arras committed
168 169 170

    def apply(self, x):
        self._check_input(x)
Martin Reinecke's avatar
Martin Reinecke committed
171
        amp = x.ptw("exp")
172
        spec = amp**2
Philipp Arras's avatar
Philipp Arras committed
173 174
        # FIXME This normalizes also the zeromode which is supposed to be left
        # untouched by this operator
Philipp Arras's avatar
Philipp Arras committed
175
        return self._specsum(self._mode_multiplicity(spec))**(-0.5)*amp
Philipp Arras's avatar
Philipp Arras committed
176 177 178


class _SpecialSum(EndomorphicOperator):
Martin Reinecke's avatar
Martin Reinecke committed
179
    def __init__(self, domain, space=0):
Philipp Arras's avatar
Philipp Arras committed
180 181
        self._domain = makeDomain(domain)
        self._capability = self.TIMES | self.ADJOINT_TIMES
182
        self._contractor = ContractionOperator(domain, space)
Philipp Arras's avatar
Philipp Arras committed
183 184 185

    def apply(self, x, mode):
        self._check_input(x, mode)
186
        return self._contractor.adjoint(self._contractor(x))
Philipp Arras's avatar
Philipp Arras committed
187 188


Philipp Haim's avatar
Philipp Haim committed
189
class _Distributor(LinearOperator):
Lukas Platz's avatar
Lukas Platz committed
190
    def __init__(self, dofdex, domain, target):
191 192 193
        self._dofdex = np.array(dofdex)
        self._target = DomainTuple.make(target)
        self._domain = DomainTuple.make(domain)
Philipp Haim's avatar
Philipp Haim committed
194 195 196 197
        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
198
        x = x.val
Philipp Haim's avatar
Philipp Haim committed
199 200 201
        if mode == self.TIMES:
            res = x[self._dofdex]
        else:
202
            res = np.zeros(self._tgt(mode).shape, dtype=x.dtype)
203
            res = utilities.special_add_at(res, 0, self._dofdex, x)
Martin Reinecke's avatar
Martin Reinecke committed
204
        return makeField(self._tgt(mode), res)
Martin Reinecke's avatar
Martin Reinecke committed
205

206

207 208
class _Amplitude(Operator):
    def __init__(self, target, fluctuations, flexibility, asperity,
209
                 loglogavgslope, azm, totvol, key, dofdex):
Philipp Arras's avatar
Philipp Arras committed
210 211
        """
        fluctuations > 0
212 213
        flexibility > 0 or None
        asperity > 0 or None
Philipp Arras's avatar
Philipp Arras committed
214 215 216
        loglogavgslope probably negative
        """
        assert isinstance(fluctuations, Operator)
217 218
        assert isinstance(flexibility, Operator) or flexibility is None
        assert isinstance(asperity, Operator) or asperity is None
Philipp Arras's avatar
Philipp Arras committed
219 220
        assert isinstance(loglogavgslope, Operator)

Philipp Haim's avatar
Philipp Haim committed
221 222
        if len(dofdex) > 0:
            N_copies = max(dofdex) + 1
Philipp Haim's avatar
Philipp Haim committed
223
            space = 1
Philipp Frank's avatar
cleanup  
Philipp Frank committed
224 225
            distributed_tgt = makeDomain((UnstructuredDomain(len(dofdex)),
                                          target))
Philipp Haim's avatar
Philipp Haim committed
226
            target = makeDomain((UnstructuredDomain(N_copies), target))
Lukas Platz's avatar
Lukas Platz committed
227
            Distributor = _Distributor(dofdex, target, distributed_tgt)
Philipp Haim's avatar
Philipp Haim committed
228
        else:
Philipp Haim's avatar
Philipp Haim committed
229
            N_copies = 0
Philipp Haim's avatar
Philipp Haim committed
230
            space = 0
Philipp Haim's avatar
Philipp Haim committed
231
            distributed_tgt = target = makeDomain(target)
Martin Reinecke's avatar
Martin Reinecke committed
232
        azm_expander = ContractionOperator(distributed_tgt, spaces=space).adjoint
Philipp Haim's avatar
Philipp Haim committed
233
        assert isinstance(target[space], PowerSpace)
Martin Reinecke's avatar
Martin Reinecke committed
234

235
        twolog = _TwoLogIntegrations(target, space)
Philipp Arras's avatar
Philipp Arras committed
236
        dom = twolog.domain
237

238
        shp = dom[space].shape
Martin Reinecke's avatar
Martin Reinecke committed
239 240
        expander = ContractionOperator(dom, spaces=space).adjoint
        ps_expander = ContractionOperator(twolog.target, spaces=space).adjoint
Philipp Arras's avatar
Philipp Arras committed
241 242 243

        # Prepare constant fields
        foo = np.zeros(shp)
244
        foo[0] = foo[1] = np.sqrt(_log_vol(target[space]))
Martin Reinecke's avatar
Martin Reinecke committed
245
        vflex = DiagonalOperator(makeField(dom[space], foo), dom, space)
Philipp Arras's avatar
Philipp Arras committed
246 247 248

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

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

255
        vslope = DiagonalOperator(
Philipp Arras's avatar
Philipp Arras committed
256
            makeField(target[space], _relative_log_k_lengths(target[space])),
Martin Reinecke's avatar
Martin Reinecke committed
257
            target, space)
258 259

        foo, bar = [np.zeros(target[space].shape) for _ in range(2)]
Philipp Arras's avatar
Philipp Arras committed
260
        bar[1:] = foo[0] = totvol
Philipp Arras's avatar
Philipp Arras committed
261 262 263 264
        vol0, vol1 = [
            DiagonalOperator(makeField(target[space], aa), target, space)
            for aa in (foo, bar)
        ]
265

Martin Reinecke's avatar
Martin Reinecke committed
266
        # Prepare fields for Adder
267
        shift, vol0 = [op(full(op.domain, 1)) for op in (shift, vol0)]
Philipp Arras's avatar
Philipp Arras committed
268 269
        # End prepare constant fields

270
        slope = vslope @ ps_expander @ loglogavgslope
271 272
        sig_flex = vflex @ expander @ flexibility if flexibility is not None else None
        sig_asp = vasp @ expander @ asperity if asperity is not None else None
273
        sig_fluc = vol1 @ ps_expander @ fluctuations
Philipp Haim's avatar
Philipp Haim committed
274
        sig_fluc = vol1 @ ps_expander @ fluctuations
Philipp Arras's avatar
Philipp Arras committed
275

276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296
        if sig_asp is None and sig_flex is None:
            op = _Normalization(target, space) @ slope
        elif sig_flex is None:
            xi = ducktape(dom, None, key)
            # TODO: this is wrong; but I know how to correct it as
            # `sigma = sig_flex * ([...] @ sig_asp)` should always be zero
            # for a `sig_flex` of zero
            sigma = (Adder(shift) @ sig_asp).ptw("sqrt")
            smooth = _SlopeRemover(target, space) @ twolog @ (sigma * xi)
            op = _Normalization(target, space) @ (slope + smooth)
        elif sig_asp is None:
            xi = ducktape(dom, None, key)
            sigma = DiagonalOperator(shift.ptw("sqrt"), dom, space) @ sig_flex
            smooth = _SlopeRemover(target, space) @ twolog @ (sigma * xi)
            op = _Normalization(target, space) @ (slope + smooth)
        else:
            xi = ducktape(dom, None, key)
            sigma = sig_flex * (Adder(shift) @ sig_asp).ptw("sqrt")
            smooth = _SlopeRemover(target, space) @ twolog @ (sigma * xi)
            op = _Normalization(target, space) @ (slope + smooth)

Philipp Haim's avatar
Philipp Haim committed
297
        if N_copies > 0:
Philipp Haim's avatar
Philipp Haim committed
298 299
            op = Distributor @ op
            sig_fluc = Distributor @ sig_fluc
Martin Reinecke's avatar
Martin Reinecke committed
300
            op = Adder(Distributor(vol0)) @ (sig_fluc*(azm_expander @ azm.ptw("reciprocal"))*op)
Philipp Arras's avatar
Philipp Arras committed
301 302
            self._fluc = (_Distributor(dofdex, fluctuations.target,
                                       distributed_tgt[0]) @ fluctuations)
Philipp Haim's avatar
Philipp Haim committed
303
        else:
Martin Reinecke's avatar
Martin Reinecke committed
304
            op = Adder(vol0) @ (sig_fluc*(azm_expander @ azm.ptw("reciprocal"))*op)
Philipp Frank's avatar
fixup  
Philipp Frank committed
305
            self._fluc = fluctuations
Philipp Arras's avatar
Philipp Arras committed
306

Philipp Arras's avatar
Philipp Arras committed
307 308
        self.apply = op.apply
        self._domain, self._target = op.domain, op.target
309
        self._space = space
310
        self._repr_str = "_Amplitude: " + op.__repr__()
Philipp Arras's avatar
Philipp Arras committed
311

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

316 317 318
    def __repr__(self):
        return self._repr_str

319 320

class CorrelatedFieldMaker:
321
    """Construction helper for hierarchical correlated field models.
Lukas Platz's avatar
Lukas Platz committed
322 323

    The correlated field models are parametrized by creating
324 325
    power spectrum operators ("amplitudes") via calls to
    :func:`add_fluctuations` that act on the targeted field subdomains.
Lukas Platz's avatar
Lukas Platz committed
326
    During creation of the :class:`CorrelatedFieldMaker` via
327 328 329
    :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
330 331

    The resulting correlated field model operator has a
Martin Reinecke's avatar
Martin Reinecke committed
332
    :class:`~nifty7.multi_domain.MultiDomain` as its domain and
Lukas Platz's avatar
Lukas Platz committed
333 334 335
    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
336
    :class:`~nifty7.domain_tuple.DomainTuple` containing the
337 338
    `target_subdomains` of the added fluctuations in the order of
    the `add_fluctuations` calls.
Lukas Platz's avatar
Lukas Platz committed
339

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

343 344 345 346 347 348 349 350 351 352
    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
353
    See the methods :func:`make`, :func:`add_fluctuations`
354
    and :func:`finalize` for further usage information."""
355 356 357
    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")
358
        self._a = []
359
        self._target_subdomains = []
Philipp Arras's avatar
Formats  
Philipp Arras committed
360

361 362
        self._offset_mean = offset_mean
        self._azm = offset_fluctuations_op
363
        self._prefix = prefix
Philipp Haim's avatar
Philipp Haim committed
364
        self._total_N = total_N
Philipp Arras's avatar
Formats  
Philipp Arras committed
365

366
    @staticmethod
Philipp Arras's avatar
Philipp Arras committed
367
    def make(offset_mean, offset_std_mean, offset_std_std, prefix, total_N=0,
Martin Reinecke's avatar
Martin Reinecke committed
368
             dofdex=None):
Lukas Platz's avatar
Lukas Platz committed
369 370 371 372 373 374 375
        """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
376
            Mean standard deviation of the offset.
Lukas Platz's avatar
Lukas Platz committed
377
        offset_std_std : float
Lukas Platz's avatar
Lukas Platz committed
378
            Standard deviation of the stddev of the offset.
Lukas Platz's avatar
Lukas Platz committed
379 380
        prefix : string
            Prefix to the names of the domains of the cf operator to be made.
Lukas Platz's avatar
Lukas Platz committed
381
            This determines the names of the operator domain.
382 383
        total_N : integer, optional
            Number of field models to create.
Lukas Platz's avatar
Lukas Platz committed
384 385 386
            If not 0, the first entry of the operators target will be an
            :class:`~nifty.domains.unstructured_domain.UnstructuredDomain`
            with length `total_N`.
387
        dofdex : np.array of integers, optional
Philipp Arras's avatar
Philipp Arras committed
388 389 390
            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
391
            instantiated; the first one is used for the first and second
392 393 394
            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
395
        """
Philipp Frank's avatar
Philipp Frank committed
396 397
        if dofdex is None:
            dofdex = np.full(total_N, 0)
398 399
        elif len(dofdex) != total_N:
            raise ValueError("length of dofdex needs to match total_N")
Philipp Frank's avatar
Philipp Frank committed
400
        N = max(dofdex) + 1 if total_N > 0 else 0
401 402
        zm = LognormalTransform(offset_std_mean, offset_std_std,
                                prefix + 'zeromode', N)
Philipp Frank's avatar
fixup  
Philipp Frank committed
403
        if total_N > 0:
Martin Reinecke's avatar
Martin Reinecke committed
404
            zm = _Distributor(dofdex, zm.target, UnstructuredDomain(total_N)) @ zm
405
        return CorrelatedFieldMaker(offset_mean, zm, prefix, total_N)
406 407

    def add_fluctuations(self,
408
                         target_subdomain,
409 410 411 412 413 414 415 416
                         fluctuations_mean,
                         fluctuations_stddev,
                         flexibility_mean,
                         flexibility_stddev,
                         asperity_mean,
                         asperity_stddev,
                         loglogavgslope_mean,
                         loglogavgslope_stddev,
Martin Reinecke's avatar
Martin Reinecke committed
417 418 419 420
                         prefix='',
                         index=None,
                         dofdex=None,
                         harmonic_partner=None):
Lukas Platz's avatar
Lukas Platz committed
421 422 423 424 425 426
        """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
427 428
        `loglogavgslope` configure the power spectrum model ("amplitude")
        used on the target field subdomain `target_subdomain`.
Lukas Platz's avatar
Lukas Platz committed
429 430 431
        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
432
        a ragged component (un-integrated wiener process).
Lukas Platz's avatar
Lukas Platz committed
433 434 435 436 437 438 439

        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
440 441
        target_subdomain : :class:`~nifty7.domain.Domain`, \
                           :class:`~nifty7.domain_tuple.DomainTuple`
Lukas Platz's avatar
Lukas Platz committed
442 443 444 445
            Target subdomain on which the correlation structure defined
            in this call should hold.
        fluctuations_{mean,stddev} : float
            Total spectral energy -> Amplitude of the fluctuations
446
        flexibility_{mean,stddev} : float or None
447
            Amplitude of the non-power-law power spectrum component
448
        asperity_{mean,stddev} : float or None
449
            Roughness of the non-power-law power spectrum component
450
            Used to accommodate single frequency peaks
Lukas Platz's avatar
Lukas Platz committed
451 452 453 454
        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
455 456 457
        index : int
            Position target_subdomain in the final total domain of the
            correlated field operator.
458 459
        dofdex : np.array, optional
            An integer array specifying the power spectrum models used if
Philipp Arras's avatar
Philipp Arras committed
460
            total_N > 1. It needs to have length of total_N. If total_N=3 and
461
            dofdex=[0,0,1], that means that two power spectrum models are
462
            instantiated; the first one is used for the first and second
463 464 465
            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
466 467
        harmonic_partner : :class:`~nifty7.domain.Domain`, \
                           :class:`~nifty7.domain_tuple.DomainTuple`
Lukas Platz's avatar
Lukas Platz committed
468 469
            In which harmonic space to define the power spectrum
        """
Philipp Frank's avatar
Philipp Frank committed
470
        if harmonic_partner is None:
471
            harmonic_partner = target_subdomain.get_default_codomain()
Philipp Frank's avatar
Fixup  
Philipp Frank committed
472
        else:
473 474
            target_subdomain.check_codomain(harmonic_partner)
            harmonic_partner.check_codomain(target_subdomain)
475

Philipp Haim's avatar
Philipp Haim committed
476 477
        if dofdex is None:
            dofdex = np.full(self._total_N, 0)
478 479
        elif len(dofdex) != self._total_N:
            raise ValueError("length of dofdex needs to match total_N")
Philipp Haim's avatar
Philipp Haim committed
480

Philipp Haim's avatar
Philipp Haim committed
481 482
        if self._total_N > 0:
            N = max(dofdex) + 1
483
            target_subdomain = makeDomain((UnstructuredDomain(N), target_subdomain))
Philipp Haim's avatar
Philipp Haim committed
484
        else:
Philipp Haim's avatar
Philipp Haim committed
485
            N = 0
486
            target_subdomain = makeDomain(target_subdomain)
Philipp Arras's avatar
Philipp Arras committed
487
        prefix = str(prefix)
488
        # assert isinstance(target_subdomain[space], (RGSpace, HPSpace, GLSpace)
Philipp Arras's avatar
Philipp Arras committed
489

490
        ve = "{0}_mean and {0}_stddev must be strictly positive (or both None)"
491 492 493 494 495
        if fluctuations_mean > 0. and fluctuations_stddev > 0.:
            fluct = LognormalTransform(fluctuations_mean, fluctuations_stddev,
                                       self._prefix + prefix + 'fluctuations', N)
        else:
            raise ValueError(ve.format("fluctuations"))
496 497 498
        if flexibility_mean is None and flexibility_stddev is None:
            flex = None
        elif flexibility_mean > 0. and flexibility_stddev > 0.:
499 500 501 502
            flex = LognormalTransform(flexibility_mean, flexibility_stddev,
                                      self._prefix + prefix + 'flexibility', N)
        else:
            raise ValueError(ve.format("flexibility"))
503 504 505
        if asperity_mean is None and asperity_stddev is None:
            asp = None
        elif asperity_mean > 0. and asperity_stddev > 0.:
506 507 508 509
            asp = LognormalTransform(asperity_mean, asperity_stddev,
                                     self._prefix + prefix + 'asperity', N)
        else:
            raise ValueError(ve.format("asperity"))
510 511 512
        avgsl = NormalTransform(loglogavgslope_mean, loglogavgslope_stddev,
                                self._prefix + prefix + 'loglogavgslope', N)

Philipp Arras's avatar
Philipp Arras committed
513
        amp = _Amplitude(PowerSpace(harmonic_partner), fluct, flex, asp, avgsl,
514
                         self._azm, target_subdomain[-1].total_volume,
515
                         self._prefix + prefix + 'spectrum', dofdex)
Philipp Haim's avatar
Philipp Haim committed
516

517 518
        if index is not None:
            self._a.insert(index, amp)
519
            self._target_subdomains.insert(index, target_subdomain)
520 521
        else:
            self._a.append(amp)
522
            self._target_subdomains.append(target_subdomain)
523

Philipp Arras's avatar
Philipp Arras committed
524 525 526 527 528 529 530 531 532 533
    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
534
        n_amplitudes = len(self._a)
Philipp Haim's avatar
Philipp Haim committed
535
        if self._total_N > 0:
Philipp Arras's avatar
Philipp Arras committed
536 537 538
            hspace = makeDomain(
                [UnstructuredDomain(self._total_N)] +
                [dd.target[-1].harmonic_partner for dd in self._a])
Philipp Haim's avatar
Philipp Haim committed
539 540
            spaces = tuple(range(1, n_amplitudes + 1))
            amp_space = 1
Philipp Haim's avatar
Philipp Haim committed
541 542
        else:
            hspace = makeDomain(
Philipp Arras's avatar
Philipp Arras committed
543
                [dd.target[0].harmonic_partner for dd in self._a])
Philipp Haim's avatar
Philipp Haim committed
544
            spaces = tuple(range(n_amplitudes))
Philipp Haim's avatar
Philipp Haim committed
545
            amp_space = 0
546

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

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

566 567
        if self._offset_mean is not None:
            offset = self._offset_mean
568 569 570 571 572 573 574
            # 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
575
        self.statistics_summary(prior_info)
576 577
        return op

578 579 580 581 582 583
    def statistics_summary(self, prior_info):
        from ..sugar import from_random

        if prior_info == 0:
            return

584 585
        lst = [('Offset amplitude', self.amplitude_total_offset),
               ('Total fluctuation amplitude', self.total_fluctuation)]
586
        namps = len(self._a)
587 588 589 590 591 592 593 594
        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:
595 596
            sc = StatCalculator()
            for _ in range(prior_info):
597
                sc.add(op(from_random(op.domain, 'normal')))
Martin Reinecke's avatar
merge  
Martin Reinecke committed
598
            mean = sc.mean.val
Martin Reinecke's avatar
Martin Reinecke committed
599
            stddev = sc.var.ptw("sqrt").val
600
            for m, s in zip(mean.flatten(), stddev.flatten()):
601
                logger.info('{}: {:.02E} ± {:.02E}'.format(kk, m, s))
602 603 604

    def moment_slice_to_average(self, fluctuations_slice_mean, nsamples=1000):
        fluctuations_slice_mean = float(fluctuations_slice_mean)
605 606 607
        if not fluctuations_slice_mean > 0:
            msg = "fluctuations_slice_mean must be greater zero; got {!r}"
            raise ValueError(msg.format(fluctuations_slice_mean))
608
        from ..sugar import from_random
609 610
        scm = 1.
        for a in self._a:
Martin Reinecke's avatar
Martin Reinecke committed
611
            op = a.fluctuation_amplitude*self._azm.ptw("reciprocal")
612
            res = np.array([op(from_random(op.domain, 'normal')).val
613 614
                            for _ in range(nsamples)])
            scm *= res**2 + 1.
615
        return fluctuations_slice_mean/np.mean(np.sqrt(scm))
616

Philipp Arras's avatar
Philipp Arras committed
617
    @property
Philipp Haim's avatar
Philipp Haim committed
618
    def normalized_amplitudes(self):
619
        """Returns the power spectrum operators used in the model"""
620
        return self._a
Philipp Arras's avatar
Philipp Arras committed
621

Philipp Haim's avatar
Philipp Haim committed
622 623 624 625 626 627 628
    @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
629 630
        dom = self._a[0].target
        expand = ContractionOperator(dom, len(dom)-1).adjoint
Philipp Haim's avatar
Philipp Haim committed
631 632
        return self._a[0]*(expand @ self.amplitude_total_offset)

633 634 635
    @property
    def amplitude_total_offset(self):
        return self._azm
Philipp Arras's avatar
Philipp Arras committed
636 637

    @property
638
    def total_fluctuation(self):
639
        """Returns operator which acts on prior or posterior samples"""
640
        if len(self._a) == 0:
Philipp Arras's avatar
Philipp Arras committed
641
            raise NotImplementedError
642
        if len(self._a) == 1:
643
            return self.average_fluctuation(0)
644 645
        q = 1.
        for a in self._a:
Martin Reinecke's avatar
Martin Reinecke committed
646
            fl = a.fluctuation_amplitude*self._azm.ptw("reciprocal")
Philipp Arras's avatar
Philipp Arras committed
647
            q = q*(Adder(full(fl.target, 1.)) @ fl**2)
Martin Reinecke's avatar
Martin Reinecke committed
648
        return (Adder(full(q.target, -1.)) @ q).ptw("sqrt")*self._azm
649

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

    def average_fluctuation(self, space):
668
        """Returns operator which acts on prior or posterior samples"""
669
        if len(self._a) == 0:
Philipp Arras's avatar
Philipp Arras committed
670
            raise NotImplementedError
671
        if space >= len(self._a):
672
            raise ValueError("invalid space specified; got {!r}".format(space))
673
        if len(self._a) == 1:
Philipp Haim's avatar
Philipp Haim committed
674 675
            return self._a[0].fluctuation_amplitude
        return self._a[space].fluctuation_amplitude
676

677 678
    @staticmethod
    def offset_amplitude_realized(samples):
Philipp Haim's avatar
Philipp Haim committed
679
        spaces = _structured_spaces(samples[0].domain)
680 681
        res = 0.
        for s in samples:
Philipp Haim's avatar
Philipp Haim committed
682
            res = res + s.mean(spaces)**2
Philipp Haim's avatar
Philipp Haim committed
683 684
        res = res/len(samples)
        return np.sqrt(res if np.isscalar(res) else res.val)
Philipp Arras's avatar
Philipp Arras committed
685

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

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