Add simple correlated field maker

src/__init__.py

@@ -88,6 +88,7 @@ from .library.adjust_variances import (make_adjust_variances_hamiltonian,
                                        do_adjust_variances)
 from .library.gridder import Gridder
 from .library.correlated_fields import CorrelatedFieldMaker
+from .library.correlated_fields_simple import SimpleCorrelatedFieldMaker
 
 from . import extra
 

src/library/correlated_fields_simple.py

+# 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
+# Authors: Philipp Arras
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+
+from functools import reduce
+from operator import mul
+
+import numpy as np
+
+from .. import utilities
+from ..domain_tuple import DomainTuple
+from ..domains.power_space import PowerSpace
+from ..domains.unstructured_domain import UnstructuredDomain
+from ..field import Field
+from ..logger import logger
+from ..multi_field import MultiField
+from ..operators.adder import Adder
+from ..operators.contraction_operator import ContractionOperator
+from ..operators.diagonal_operator import DiagonalOperator
+from ..operators.distributors import PowerDistributor
+from ..operators.endomorphic_operator import EndomorphicOperator
+from ..operators.harmonic_operators import HarmonicTransformOperator
+from ..operators.linear_operator import LinearOperator
+from ..operators.operator import Operator
+from ..operators.simple_linear_operators import ducktape
+from ..operators.normal_operators import NormalTransform, LognormalTransform
+from ..probing import StatCalculator
+from ..sugar import full, makeDomain, makeField, makeOp
+
+
+class SimpleCorrelatedFieldMaker:
+    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):
+        from .correlated_fields import _Amplitude
+        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)
+        amp = _Amplitude(PowerSpace(harmonic_partner), fluct, flex, asp, avgsl,
+                         zm, target.total_volume, prefix + 'spectrum', [])
+        ht = HarmonicTransformOperator(harmonic_partner, target)
+        pd = PowerDistributor(harmonic_partner, amp.target[0])
+        expander = ContractionOperator(harmonic_partner, spaces=0).adjoint
+        self._op = ht(
+            expander(zm)*pd(amp)*
+            ducktape(harmonic_partner, None, prefix + 'xi'))
+        if offset_mean is not None:
+            self._op = Adder(full(self._op.target,
+                                  float(offset_mean))) @ self._op
+
+    def finalize(self):
+        return self._op

test/test_operators/test_correlated_fields.py

@@ -15,6 +15,7 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
+import numpy as np
 import pytest
 from numpy.testing import assert_, assert_allclose
 
@@ -23,6 +24,12 @@ import nifty7 as ift
 from ..common import setup_function, teardown_function
 
 pmp = pytest.mark.parametrize
+spaces = [
+    ift.RGSpace(4),
+    ift.RGSpace((4, 4), (0.123, 0.4)),
+    ift.HPSpace(8),
+    ift.GLSpace(4)
+]
 
 
 def _stats(op, samples):
@@ -32,6 +39,14 @@ def _stats(op, samples):
     return sc.mean.val, sc.var.ptw("sqrt").val
 
 
+def _rand():
+    return ift.random.current_rng().normal()
+
+
+def _posrand():
+    return np.exp(_rand())
+
+
 @pmp('dofdex', [[0, 0], [0, 1]])
 @pmp('seed', [12, 3])
 def testDistributor(dofdex, seed):
@@ -46,12 +61,7 @@ def testDistributor(dofdex, seed):
         ift.extra.check_linear_operator(op)
 
 
-@pmp('sspace', [
-    ift.RGSpace(4),
-    ift.RGSpace((4, 4), (0.123, 0.4)),
-    ift.HPSpace(8),
-    ift.GLSpace(4)
-])
+@pmp('sspace', spaces)
 @pmp('N', [0, 2])
 def testAmplitudesInvariants(sspace, N):
     fsspace = ift.RGSpace((12,), (0.4,))
@@ -110,3 +120,54 @@ def testAmplitudesInvariants(sspace, N):
     for ampl in fa.normalized_amplitudes:
         ift.extra.check_operator(ampl, 0.1*ift.from_random(ampl.domain), ntries=10)
     ift.extra.check_operator(op, 0.1*ift.from_random(op.domain), ntries=10)
+
+
+@pmp('seed', [42, 31])
+@pmp('domain', spaces)
+def test_complicated_vs_simple(seed, domain):
+    with ift.random.Context(seed):
+        offset_mean = _rand()
+        offset_std_mean = _posrand()
+        offset_std_std = _posrand()
+        fluctuations_mean = _posrand()
+        fluctuations_stddev = _posrand()
+        flexibility_mean = _posrand()
+        flexibility_stddev = _posrand()
+        asperity_mean = _posrand()
+        asperity_stddev = _posrand()
+        loglogavgslope_mean = _posrand()
+        loglogavgslope_stddev = _posrand()
+        prefix = 'foobar'
+        hspace = domain.get_default_codomain()
+        cfm_simple = ift.SimpleCorrelatedFieldMaker(domain,
+                                                    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=prefix,
+                                                    harmonic_partner=hspace)
+        cfm = ift.CorrelatedFieldMaker.make(offset_mean, offset_std_mean,
+                                            offset_std_std, prefix)
+        cfm.add_fluctuations(domain,
+                             fluctuations_mean,
+                             fluctuations_stddev,
+                             flexibility_mean,
+                             flexibility_stddev,
+                             asperity_mean,
+                             asperity_stddev,
+                             loglogavgslope_mean,
+                             loglogavgslope_stddev,
+                             prefix='',
+                             harmonic_partner=hspace)
+        op1 = cfm.finalize()
+        op0 = cfm_simple.finalize()
+        assert_(op0.domain is op1.domain)
+        inp = ift.from_random(op0.domain)
+        ift.extra.assert_allclose(op0(inp), op1(inp))