correlated_fields_simple.py 4.97 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14
# 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/>.
#
# Copyright(C) 2013-2020 Max-Planck-Society
Philipp Arras's avatar
Philipp Arras committed
15
# Author: Philipp Arras
16 17 18 19 20 21 22 23 24 25
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.

import numpy as np

from ..domains.power_space import PowerSpace
from ..operators.adder import Adder
from ..operators.contraction_operator import ContractionOperator
from ..operators.distributors import PowerDistributor
from ..operators.harmonic_operators import HarmonicTransformOperator
Philipp Arras's avatar
Philipp Arras committed
26
from ..operators.normal_operators import LognormalTransform, NormalTransform
Philipp Arras's avatar
Philipp Arras committed
27
from .correlated_fields import _TwoLogIntegrations
28 29
from ..operators.operator import Operator
from ..operators.simple_linear_operators import ducktape
Philipp Arras's avatar
Philipp Arras committed
30 31
from ..operators.value_inserter import ValueInserter
from ..sugar import full, makeField, makeOp
Philipp Arras's avatar
Philipp Arras committed
32
from .correlated_fields import (_log_vol, _Normalization,
Philipp Arras's avatar
Philipp Arras committed
33 34 35
                                _relative_log_k_lengths, _SlopeRemover)


Philipp Arras's avatar
Philipp Arras committed
36
class SimpleCorrelatedField(Operator):
Philipp Arras's avatar
Philipp Arras committed
37 38 39 40 41 42
    """Simplified version of :class:`~nifty7.library.correlated_fields.CorrelatedFieldMaker`.

    Assumes `total_N = 0`, `dofdex = None` and the presence of only one power
    spectrum, i.e. only one call of
    :func:`~nifty7.library.correlated_fields.CorrelatedFieldMaker.add_fluctuations`.
    """
43 44
    def __init__(self, target, offset_mean, offset_std, fluctuations,
                 flexibility, asperity, loglogavgslope, prefix='',
45 46 47 48 49 50
                 harmonic_partner=None):
        if harmonic_partner is None:
            harmonic_partner = target.get_default_codomain()
        else:
            target.check_codomain(harmonic_partner)
            harmonic_partner.check_codomain(target)
51 52 53 54
        for kk in [offset_std, fluctuations, loglogavgslope]:
            if len(kk) != 2:
                raise TypeError
        for kk in [flexibility, asperity]:
55
            if not (kk is None or len(kk) == 2):
56 57 58 59 60 61
                raise TypeError
        if flexibility is None and asperity is not None:
            raise ValueError
        fluct = LognormalTransform(*fluctuations, prefix + 'fluctuations', 0)
        avgsl = NormalTransform(*loglogavgslope, prefix + 'loglogavgslope', 0)
        zm = LognormalTransform(*offset_std, prefix + 'zeromode', 0)
Philipp Arras's avatar
Philipp Arras committed
62

Philipp Arras's avatar
Philipp Arras committed
63 64
        pspace = PowerSpace(harmonic_partner)
        twolog = _TwoLogIntegrations(pspace)
Philipp Arras's avatar
Philipp Arras committed
65
        expander = ContractionOperator(twolog.domain, 0).adjoint
Philipp Arras's avatar
Philipp Arras committed
66
        ps_expander = ContractionOperator(pspace, 0).adjoint
67
        vslope = makeOp(makeField(pspace, _relative_log_k_lengths(pspace)))
Philipp Arras's avatar
Philipp Arras committed
68
        slope = vslope @ ps_expander @ avgsl
69 70 71 72 73 74 75 76 77 78 79 80 81 82
        a = slope

        if flexibility is not None:
            flex = LognormalTransform(*flexibility, prefix + 'flexibility', 0)
            dom = twolog.domain[0]
            vflex = np.zeros(dom.shape)
            vflex[0] = vflex[1] = np.sqrt(_log_vol(pspace))
            vflex = makeOp(makeField(dom, vflex))
            sig_flex = vflex @ expander @ flex
            xi = ducktape(dom, None, prefix + 'spectrum')

            shift = np.ones(dom.shape)
            shift[0] = _log_vol(pspace)**2/12.
            shift = makeField(dom, shift)
Philipp Frank's avatar
Philipp Frank committed
83
            if asperity is None:
Philipp Arras's avatar
Philipp Arras committed
84
                asp = makeOp(shift.ptw("sqrt")) @ (xi*sig_flex)
Philipp Frank's avatar
Philipp Frank committed
85 86 87 88 89 90
            else:
                asp = LognormalTransform(*asperity, prefix + 'asperity', 0)
                vasp = np.zeros(dom.shape)
                vasp[0] = 1
                vasp = makeOp(makeField(dom, vasp))
                sig_asp = vasp @ expander @ asp
Philipp Arras's avatar
Philipp Arras committed
91 92
                asp = xi*sig_flex*(Adder(shift) @ sig_asp).ptw("sqrt")
            a = a + _SlopeRemover(pspace, 0) @ twolog @ asp
93
        a = _Normalization(pspace, 0) @ a
Philipp Arras's avatar
Philipp Arras committed
94
        maskzm = np.ones(pspace.shape)
Philipp Arras's avatar
Philipp Arras committed
95
        maskzm[0] = 0
Philipp Arras's avatar
Philipp Arras committed
96 97
        maskzm = makeOp(makeField(pspace, maskzm))
        insert = ValueInserter(pspace, (0,))
Philipp Arras's avatar
Philipp Arras committed
98
        a = (maskzm @ ((ps_expander @ fluct)*a)) + insert(zm)
Philipp Arras's avatar
Philipp Arras committed
99
        self._a = a.scale(target.total_volume)
Philipp Arras's avatar
Philipp Arras committed
100

101
        ht = HarmonicTransformOperator(harmonic_partner, target)
Philipp Arras's avatar
Philipp Arras committed
102
        pd = PowerDistributor(harmonic_partner, pspace)
Philipp Arras's avatar
Philipp Arras committed
103
        xi = ducktape(harmonic_partner, None, prefix + 'xi')
Philipp Arras's avatar
Philipp Arras committed
104
        op = ht(pd(self._a)*xi)
105
        if offset_mean is not None:
Philipp Arras's avatar
Philipp Arras committed
106 107 108 109
            op = Adder(full(op.target, float(offset_mean))) @ op
        self.apply = op.apply
        self._domain = op.domain
        self._target = op.target
110 111 112

    @property
    def amplitude(self):
Philipp Arras's avatar
Philipp Arras committed
113
        """Analoguous to :func:`~nifty7.library.correlated_fields.CorrelatedFieldMaker.amplitude`."""
114
        return self._a