Merge branch 'simple_cfm' into 'NIFTy_7'

Simple cfm

See merge request !552

.gitlab-ci.yml

@@ -53,6 +53,7 @@ test_mpi:
 pages:
   stage: release
   script:
+    - rm -rf docs/build docs/source/mod
     - sh docs/generate.sh
     - mv docs/build/ public/
   artifacts:

demos/getting_started_3.py

@@ -52,20 +52,17 @@ def main():
     else:
         mode = 0
     filename = "getting_started_3_mode_{}_".format(mode) + "{}.png"
-
     position_space = ift.RGSpace([128, 128])
 
     #  For a detailed showcase of the effects the parameters
     #  of the CorrelatedField model have on the generated fields,
     #  see 'getting_started_4_CorrelatedFields.ipynb'.
 
-    cfmaker = ift.CorrelatedFieldMaker.make(
-        offset_mean=      0.0,
-        offset_std_mean= 1e-3,
-        offset_std_std=  1e-6,
-        prefix='')
+    args = {
+        'offset_mean': 0,
+        'offset_std_mean': 1e-3,
+        'offset_std_std': 1e-6,
 
-    fluctuations_dict = {
         # Amplitude of field fluctuations
         'fluctuations_mean':   2.0,  # 1.0
         'fluctuations_stddev': 1.0,  # 1e-2
@@ -82,10 +79,9 @@ def main():
         'asperity_mean':   0.5,  # 0.1
         'asperity_stddev': 0.5  # 0.5
     }
-    cfmaker.add_fluctuations(position_space, **fluctuations_dict)
 
-    correlated_field = cfmaker.finalize()
-    A = cfmaker.amplitude
+    correlated_field = ift.SimpleCorrelatedField(position_space, **args)
+    A = correlated_field.amplitude
 
     # Apply a nonlinearity
     signal = ift.sigmoid(correlated_field)
@@ -143,8 +139,6 @@ def main():
         plot.output(ny=1, ysize=6, xsize=16,
                     name=filename.format("loop_{:02d}".format(i)))
 
-    # Draw posterior samples
-    KL = ift.MetricGaussianKL.make(mean, H, N_samples)
     sc = ift.StatCalculator()
     for sample in KL.samples:
         sc.add(signal(sample + KL.position))
@@ -156,11 +150,15 @@ def main():
     plot.add(ift.sqrt(sc.var), title="Posterior Standard Deviation")
 
     powers = [A.force(s + KL.position) for s in KL.samples]
+    sc = ift.StatCalculator()
+    for pp in powers:
+        sc.add(pp)
     plot.add(
         powers + [A.force(mock_position),
-                  A.force(KL.position)],
+                  A.force(KL.position), sc.mean],
         title="Sampled Posterior Power Spectrum",
-        linewidth=[1.]*len(powers) + [3., 3.])
+        linewidth=[1.]*len(powers) + [3., 3., 3.],
+        label=[None]*len(powers) + ['Ground truth', 'Posterior latent mean', 'Posterior mean'])
     plot.output(ny=1, nx=3, xsize=24, ysize=6, name=filename_res)
     print("Saved results as '{}'.".format(filename_res))
 

docs/generate.sh

-rm -rf docs/build docs/source/mod
 EXCLUDE="nifty7/logger.py nifty7/git_version.py"
 sphinx-apidoc -e -o docs/source/mod nifty7 ${EXCLUDE}
 sphinx-build -b html docs/source/ docs/build/

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 SimpleCorrelatedField
 
 from . import extra
 

src/extra.py

@@ -20,7 +20,6 @@ from itertools import combinations
 import numpy as np
 from numpy.testing import assert_
 
-from . import random
 from .domain_tuple import DomainTuple
 from .field import Field
 from .linearization import Linearization
@@ -124,7 +123,7 @@ def check_operator(op, loc, tol=1e-12, ntries=100, perf_check=True,
                                metric_sampling)
 
 
-def assert_allclose(f1, f2, atol, rtol):
+def assert_allclose(f1, f2, atol=0, rtol=1e-7):
     if isinstance(f1, Field):
         return np.testing.assert_allclose(f1.val, f2.val, atol=atol, rtol=rtol)
     for key, val in f1.items():

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
+# Author: 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
+from ..operators.normal_operators import LognormalTransform, NormalTransform
+from .correlated_fields import _TwoLogIntegrations
+from ..operators.operator import Operator
+from ..operators.simple_linear_operators import ducktape
+from ..operators.value_inserter import ValueInserter
+from ..sugar import full, makeField, makeOp
+from .correlated_fields import (_log_vol, _Normalization,
+                                _relative_log_k_lengths, _SlopeRemover)
+
+
+class SimpleCorrelatedField(Operator):
+    """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`.
+    """
+    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)
+
+        pspace = PowerSpace(harmonic_partner)
+        twolog = _TwoLogIntegrations(pspace)
+        dom = twolog.domain[0]
+        vflex = np.zeros(dom.shape)
+        vasp = np.zeros(dom.shape)
+        shift = np.ones(dom.shape)
+        vflex[0] = vflex[1] = np.sqrt(_log_vol(pspace))
+        vasp[0] = 1
+        shift[0] = _log_vol(pspace)**2/12.
+        vflex = makeOp(makeField(dom, vflex))
+        vasp = makeOp(makeField(dom, vasp))
+        shift = makeField(dom, shift)
+        vslope = makeOp(makeField(pspace, _relative_log_k_lengths(pspace)))
+
+        expander = ContractionOperator(twolog.domain, 0).adjoint
+        ps_expander = ContractionOperator(pspace, 0).adjoint
+        slope = vslope @ ps_expander @ avgsl
+        sig_flex = vflex @ expander @ flex
+        sig_asp = vasp @ expander @ asp
+        xi = ducktape(dom, None, prefix + 'spectrum')
+        smooth = xi*sig_flex*(Adder(shift) @ sig_asp).ptw("sqrt")
+        smooth = _SlopeRemover(pspace, 0) @ twolog @ smooth
+        a = _Normalization(pspace, 0) @ (slope + smooth)
+
+        maskzm = np.ones(pspace.shape)
+        maskzm[0] = 0
+        maskzm = makeOp(makeField(pspace, maskzm))
+        insert = ValueInserter(pspace, (0,))
+        a = (maskzm @ ((ps_expander @ fluct)*a)) + insert(zm)
+        self._a = a.scale(target.total_volume)
+
+        ht = HarmonicTransformOperator(harmonic_partner, target)
+        pd = PowerDistributor(harmonic_partner, pspace)
+        xi = ducktape(harmonic_partner, None, prefix + 'xi')
+        op = ht(pd(self._a)*xi)
+        if offset_mean is not None:
+            op = Adder(full(op.target, float(offset_mean))) @ op
+        self.apply = op.apply
+        self._domain = op.domain
+        self._target = op.target
+
+    @property
+    def amplitude(self):
+        """Analoguous to :func:`~nifty7.library.correlated_fields.CorrelatedFieldMaker.amplitude`."""
+        return self._a

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,59 @@ 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()
+        scf = ift.SimpleCorrelatedField(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)
+        inp = ift.from_random(scf.domain)
+        op1 = cfm.finalize()
+        assert_(scf.domain is op1.domain)
+        ift.extra.assert_allclose(scf(inp), op1(inp))
+        ift.extra.check_operator(scf, inp, ntries=10)
+
+        op1 = cfm.amplitude
+        op0 = scf.amplitude
+        assert_(op0.domain is op1.domain)
+        ift.extra.assert_allclose(op0.force(inp), op1.force(inp))