correlated_fields_simple.py 4.55 KB
 Philipp Arras committed Jun 27, 2020 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 ``````# 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 . # # Copyright(C) 2013-2020 Max-Planck-Society # Authors: Philipp Arras # # 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 committed Jun 27, 2020 26 ``````from ..operators.normal_operators import LognormalTransform, NormalTransform `````` Philipp Arras committed Jun 27, 2020 27 ``````from .correlated_fields import _TwoLogIntegrations `````` Philipp Arras committed Jun 27, 2020 28 29 ``````from ..operators.operator import Operator from ..operators.simple_linear_operators import ducktape `````` Philipp Arras committed Jun 27, 2020 30 31 ``````from ..operators.value_inserter import ValueInserter from ..sugar import full, makeField, makeOp `````` Philipp Arras committed Jun 27, 2020 32 ``````from .correlated_fields import (_log_vol, _Normalization, `````` Philipp Arras committed Jun 27, 2020 33 34 35 `````` _relative_log_k_lengths, _SlopeRemover) `````` Philipp Arras committed Jun 27, 2020 36 ``````class SimpleCorrelatedField(Operator): `````` Philipp Arras committed Jun 27, 2020 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 `````` def __init__(self, target, offset_mean, offset_std_mean, offset_std_std, fluctuations_mean, fluctuations_stddev, flexibility_mean, flexibility_stddev, asperity_mean, asperity_stddev, loglogavgslope_mean, loglogavgslope_stddev, prefix='', 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) fluct = LognormalTransform(fluctuations_mean, fluctuations_stddev, prefix + 'fluctuations', 0) flex = LognormalTransform(flexibility_mean, flexibility_stddev, prefix + 'flexibility', 0) asp = LognormalTransform(asperity_mean, asperity_stddev, prefix + 'asperity', 0) avgsl = NormalTransform(loglogavgslope_mean, loglogavgslope_stddev, prefix + 'loglogavgslope', 0) zm = LognormalTransform(offset_std_mean, offset_std_std, prefix + 'zeromode', 0) `````` Philipp Arras committed Jun 27, 2020 67 68 `````` tgt = PowerSpace(harmonic_partner) `````` Philipp Arras committed Jun 27, 2020 69 `````` twolog = _TwoLogIntegrations(tgt) `````` Philipp Arras committed Jun 27, 2020 70 71 72 73 74 75 76 77 78 `````` dom = twolog.domain[0] vflex = np.zeros(dom.shape) vasp = np.zeros(dom.shape, dtype=np.float64) shift = np.ones(dom.shape, dtype=np.float64) vflex[0] = vflex[1] = np.sqrt(_log_vol(tgt)) vasp[0] = 1 shift[0] = _log_vol(tgt)**2/12. vflex = makeOp(makeField(dom, vflex)) vasp = makeOp(makeField(dom, vasp)) `````` Philipp Arras committed Jun 27, 2020 79 `````` shift = makeField(dom, shift) `````` Philipp Arras committed Jun 27, 2020 80 `````` vslope = makeOp(makeField(tgt, _relative_log_k_lengths(tgt))) `````` Philipp Arras committed Jun 27, 2020 81 82 `````` expander = ContractionOperator(twolog.domain, 0).adjoint `````` Philipp Arras committed Jun 27, 2020 83 `````` ps_expander = ContractionOperator(tgt, 0).adjoint `````` Philipp Arras committed Jun 27, 2020 84 85 86 87 `````` slope = vslope @ ps_expander @ avgsl sig_flex = vflex @ expander @ flex sig_asp = vasp @ expander @ asp xi = ducktape(dom, None, prefix + 'spectrum') `````` Philipp Arras committed Jun 27, 2020 88 89 `````` smooth = xi*sig_flex*(Adder(shift) @ sig_asp).ptw("sqrt") smooth = _SlopeRemover(tgt, 0) @ twolog @ smooth `````` Philipp Arras committed Jun 27, 2020 90 `````` op = _Normalization(tgt, 0) @ (slope + smooth) `````` Philipp Arras committed Jun 27, 2020 91 `````` `````` Philipp Arras committed Jun 27, 2020 92 93 94 95 96 97 `````` maskzm = np.ones(tgt.shape) maskzm[0] = 0 maskzm = makeOp(makeField(tgt, maskzm)) insert = ValueInserter(tgt, (0,)) a = (maskzm @ ((ps_expander @ fluct)*op)) + insert(zm) self._a = a.scale(target.total_volume) `````` Philipp Arras committed Jun 27, 2020 98 `````` `````` Philipp Arras committed Jun 27, 2020 99 `````` ht = HarmonicTransformOperator(harmonic_partner, target) `````` Philipp Arras committed Jun 27, 2020 100 `````` pd = PowerDistributor(harmonic_partner, tgt) `````` Philipp Arras committed Jun 27, 2020 101 `````` xi = ducktape(harmonic_partner, None, prefix + 'xi') `````` Philipp Arras committed Jun 27, 2020 102 `````` op = ht(pd(self._a)*xi) `````` Philipp Arras committed Jun 27, 2020 103 `````` if offset_mean is not None: `````` Philipp Arras committed Jun 27, 2020 104 105 106 107 `````` op = Adder(full(op.target, float(offset_mean))) @ op self.apply = op.apply self._domain = op.domain self._target = op.target``````