merge NIFTy_6

.gitlab-ci.yml

@@ -43,6 +43,13 @@ test_serial:
     - >
       grep TOTAL coverage.txt | awk '{ print "TOTAL: "$4; }'
 
+test_mpi:
+  stage: test
+  variables:
+    OMPI_MCA_btl_vader_single_copy_mechanism: none
+  script:
+    - mpiexec -n 2 --bind-to none pytest-3 -q test/test_mpi
+
 pages:
   stage: release
   script:
@@ -61,7 +68,7 @@ before_script:
 run_ipynb:
   stage: demo_runs
   script:
-    - jupyter nbconvert --execute --ExecutePreprocessor.timeout=None demos/Wiener_Filter.ipynb
+    - jupyter nbconvert --execute --ExecutePreprocessor.timeout=None demos/getting_started_0.ipynb
 
 run_getting_started_1:
   stage: demo_runs

ChangeLog

 Changes since NIFTy 5:
 
+Minimum Python version increased to 3.6
+=======================================
+
+
+New operators
+=============
+
+In addition to the below changes, the following operators were introduced:
+
+* UniformOperator: Transforms a Gaussian into a uniform distribution
+* VariableCovarianceGaussianEnergy: Energy operator for inferring covariances
+* MultiLinearEinsum: Multi-linear version of numpy's einsum with derivates
+* LinearEinsum: Linear version of numpy's einsum with one free field
+* PartialConjugate: Conjugates parts of a multi-field
+* SliceOperator: Geometry preserving mask operator
+* SplitOperator: Splits a single field into a multi-field
+
+FFT convention adjusted
+=======================
+
+When going to harmonic space, NIFTy's FFT operator now uses a minus sign in the
+exponent (and, consequently, a plus sign on the adjoint transform). This
+convention is consistent with almost all other numerical FFT libraries.
+
+Interface change in EndomorphicOperator.draw_sample()
+=====================================================
+
+Both complex-valued and real-valued Gaussian probability distributions have
+Hermitian and positive endomorphisms as covariance. Just by looking at an
+endomorphic operator itself it is not clear whether it is viewed as covariance
+for real or complex Gaussians when a sample of the respective distribution shall
+be drawn. Therefore, we introduce the method `draw_sample_with_dtype()` which
+needs to be given the data type of the probability distribution. This function
+is implemented for all operators which actually draw random numbers
+(`DiagonalOperator` and `ScalingOperator`). The class `SamplingDtypeSetter` acts
+as a wrapper for this kind of operators in order to fix the data type of the
+distribution. Samples from these operators can be drawn with `.draw_sample()`.
+In order to dive into those subtleties I suggest running the following code and
+playing around with the dtypes.
+
+```
+import nifty6 as ift
+import numpy as np
+
+dom = ift.UnstructuredDomain(5)
+dtype = [np.float64, np.complex128][1]
+invcov = ift.ScalingOperator(dom, 3)
+e = ift.GaussianEnergy(mean=ift.from_random(dom, 'normal', dtype=dtype),
+                       inverse_covariance=invcov)
+pos = ift.from_random(dom, 'normal', dtype=np.complex128)
+lin = e(ift.Linearization.make_var(pos, want_metric=True))
+met = lin.metric
+print(met)
+print(met.draw_sample())
+```
+
 MPI parallelisation over samples in MetricGaussianKL
 ====================================================
 
@@ -15,6 +71,13 @@ the generation of reproducible random numbers in the presence of MPI parallelism
 and leads to cleaner code overall. Please see the documentation of
 `nifty6.random` for details.
 
+Interface Change for from_random and OuterProduct
+=================================================
+
+The sugar.from_random, Field.from_random, MultiField.from_random now take domain
+as the first argument and default to 'normal' for the second argument.
+Likewise OuterProduct takes domain as the first argument and a field as the second.
+
 Interface Change for non-linear Operators
 =========================================
 

README.md

@@ -45,7 +45,7 @@ Installation
 
 ### Requirements
 
-- [Python 3](https://www.python.org/) (3.5.x or later)
+- [Python 3](https://www.python.org/) (3.6.x or later)
 - [SciPy](https://www.scipy.org/)
 
 Optional dependencies:

demos/bernoulli_demo.py

@@ -43,7 +43,7 @@ if __name__ == '__main__':
     harmonic_space = position_space.get_default_codomain()
     HT = ift.HarmonicTransformOperator(harmonic_space, position_space)
 
-    position = ift.from_random('normal', harmonic_space)
+    position = ift.from_random(harmonic_space, 'normal')
 
     # Define power spectrum and amplitudes
     def sqrtpspec(k):
@@ -58,13 +58,13 @@ if __name__ == '__main__':
 
     # Generate mock data
     p = R(sky)
-    mock_position = ift.from_random('normal', harmonic_space)
+    mock_position = ift.from_random(harmonic_space, 'normal')
     tmp = p(mock_position).val.astype(np.float64)
     data = ift.random.current_rng().binomial(1, tmp)
     data = ift.Field.from_raw(R.target, data)
 
     # Compute likelihood and Hamiltonian
-    position = ift.from_random('normal', harmonic_space)
+    position = ift.from_random(harmonic_space, 'normal')
     likelihood = ift.BernoulliEnergy(data) @ p
     ic_newton = ift.DeltaEnergyController(
         name='Newton', iteration_limit=100, tol_rel_deltaE=1e-8)

demos/Wiener_Filter.ipynb → demos/getting_started_0.ipynb

@@ -236,7 +236,7 @@
     "R = HT #*ift.create_harmonic_smoothing_operator((h_space,), 0, 0.02)\n",
     "\n",
     "# Fields and data\n",
-    "sh = Sh.draw_sample()\n",
+    "sh = Sh.draw_sample_with_dtype(dtype=np.float64)\n",
     "noiseless_data=R(sh)\n",
     "noise_amplitude = np.sqrt(0.2)\n",
     "N = ift.ScalingOperator(s_space, noise_amplitude**2)\n",
@@ -394,7 +394,7 @@
     "# R is defined below\n",
     "\n",
     "# Fields\n",
-    "sh = Sh.draw_sample()\n",
+    "sh = Sh.draw_sample_with_dtype(dtype=np.float64)\n",
     "s = HT(sh)\n",
     "n = ift.Field.from_random(domain=s_space, random_type='normal',\n",
     "                      std=noise_amplitude, mean=0)"
@@ -471,7 +471,7 @@
    },
    "outputs": [],
    "source": [
-    "m_mean, m_var = ift.probe_with_posterior_samples(curv, HT, 200)"
+    "m_mean, m_var = ift.probe_with_posterior_samples(curv, HT, 200, np.float64)"
    ]
   },
   {
@@ -571,7 +571,7 @@
     "N = ift.ScalingOperator(s_space, sigma2)\n",
     "\n",
     "# Fields and data\n",
-    "sh = Sh.draw_sample()\n",
+    "sh = Sh.draw_sample_with_dtype(dtype=np.float64)\n",
     "n = ift.Field.from_random(domain=s_space, random_type='normal',\n",
     "                      std=np.sqrt(sigma2), mean=0)\n",
     "\n",
@@ -598,7 +598,7 @@
     "m = D(j)\n",
     "\n",
     "# Uncertainty\n",
-    "m_mean, m_var = ift.probe_with_posterior_samples(curv, HT, 20)\n",
+    "m_mean, m_var = ift.probe_with_posterior_samples(curv, HT, 20, np.float64)\n",
     "\n",
     "# Get data\n",
     "s_data = HT(sh).val\n",
@@ -709,8 +709,15 @@
     "\n",
     "https://gitlab.mpcdf.mpg.de/ift/NIFTy\n",
     "\n",
-    "NIFTy v5 **more or less stable!**"
+    "NIFTy v6 **more or less stable!**"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
@@ -730,7 +737,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.5"
+   "version": "3.8.2"
   }
  },
  "nbformat": 4,

demos/getting_started_1.py

@@ -40,7 +40,7 @@ def make_checkerboard_mask(position_space):
 
 def make_random_mask():
     # Random mask for spherical mode
-    mask = ift.from_random('pm1', position_space)
+    mask = ift.from_random(position_space, 'pm1')
     mask = (mask + 1)/2
     return mask.val
 
@@ -114,8 +114,8 @@ if __name__ == '__main__':
     N = ift.ScalingOperator(data_space, noise)
 
     # Create mock data
-    MOCK_SIGNAL = S.draw_sample()
-    MOCK_NOISE = N.draw_sample()
+    MOCK_SIGNAL = S.draw_sample_with_dtype(dtype=np.float64)
+    MOCK_NOISE = N.draw_sample_with_dtype(dtype=np.float64)
     data = R(MOCK_SIGNAL) + MOCK_NOISE
 
     # Build inverse propagator D and information source j

demos/getting_started_2.py

@@ -90,7 +90,7 @@ if __name__ == '__main__':
     # Generate mock data and define likelihood operator
     d_space = R.target[0]
     lamb = R(sky)
-    mock_position = ift.from_random('normal', domain)
+    mock_position = ift.from_random(domain, 'normal')
     data = lamb(mock_position)
     data = ift.random.current_rng().poisson(data.val.astype(np.float64))
     data = ift.Field.from_raw(d_space, data)
@@ -103,7 +103,7 @@ if __name__ == '__main__':
 
     # Compute MAP solution by minimizing the information Hamiltonian
     H = ift.StandardHamiltonian(likelihood)
-    initial_position = ift.from_random('normal', domain)
+    initial_position = ift.from_random(domain, 'normal')
     H = ift.EnergyAdapter(initial_position, H, want_metric=True)
     H, convergence = minimizer(H)
 

demos/getting_started_3.py

@@ -55,9 +55,32 @@ if __name__ == '__main__':
 
     position_space = ift.RGSpace([128, 128])
 
-    cfmaker = ift.CorrelatedFieldMaker.make(1e-3, 1e-6, '')
-    cfmaker.add_fluctuations(position_space,
-                             1., 1e-2, 1, .5, .1, .5, -3, 0.5, '')
+    cfmaker = ift.CorrelatedFieldMaker.make(
+        offset_mean =      0.0,  # 0.
+        offset_std_mean = 1e-3,  # 1e-3
+        offset_std_std =  1e-6,  # 1e-6
+        prefix = '')
+
+    fluctuations_dict = {
+        # Amplitude of the fluctuations
+        'fluctuations_mean':   2.0,  # 1.0
+        'fluctuations_stddev': 1.0,  # 1e-2
+
+        # Smooth variation speed
+        'flexibility_mean':   2.5,  # 1.0
+        'flexibility_stddev': 1.0,  # 0.5
+
+        # How strong the ragged component of the spectrum is
+        # (Ratio of Wiener process and integrated Wiener process ?)
+        'asperity_mean':   0.5,  # 0.1
+        'asperity_stddev': 0.5,  # 0.5
+
+        # Slope of linear spectrum component
+        'loglogavgslope_mean': -2.0,  # -3.0
+        'loglogavgslope_stddev': 0.5  #  0.5
+    }
+    cfmaker.add_fluctuations(position_space, **fluctuations_dict)
+
     correlated_field = cfmaker.finalize()
     A = cfmaker.amplitude
 
@@ -75,8 +98,8 @@ if __name__ == '__main__':
     N = ift.ScalingOperator(data_space, noise)
 
     # Generate mock signal and data
-    mock_position = ift.from_random('normal', signal_response.domain)
-    data = signal_response(mock_position) + N.draw_sample()
+    mock_position = ift.from_random(signal_response.domain, 'normal')
+    data = signal_response(mock_position) + N.draw_sample_with_dtype(dtype=np.float64)
 
     # Minimization parameters
     ic_sampling = ift.AbsDeltaEnergyController(

demos/getting_started_mf.py

@@ -73,7 +73,7 @@ if __name__ == '__main__':
     sp2 = ift.RGSpace(npix2)
 
     # Set up signal model
-    cfmaker = ift.CorrelatedFieldMaker.make(1e-2, 1e-6, '')
+    cfmaker = ift.CorrelatedFieldMaker.make(0., 1e-2, 1e-6, '')
     cfmaker.add_fluctuations(sp1, 0.1, 1e-2, 1, .1, .01, .5, -2, 1., 'amp1')
     cfmaker.add_fluctuations(sp2, 0.1, 1e-2, 1, .1, .01, .5,
                              -1.5, .5, 'amp2')
@@ -97,8 +97,8 @@ if __name__ == '__main__':
     N = ift.ScalingOperator(data_space, noise)
 
     # Generate mock signal and data
-    mock_position = ift.from_random('normal', signal_response.domain)
-    data = signal_response(mock_position) + N.draw_sample()
+    mock_position = ift.from_random(signal_response.domain, 'normal')
+    data = signal_response(mock_position) + N.draw_sample_with_dtype(dtype=np.float64)
 
     plot = ift.Plot()
     plot.add(signal(mock_position), title='Ground Truth')
@@ -114,7 +114,9 @@ if __name__ == '__main__':
     ic_newton = ift.AbsDeltaEnergyController(name='Newton',
                                              deltaE=0.01,
                                              iteration_limit=35)
-    minimizer = ift.NewtonCG(ic_newton)
+    ic_sampling.enable_logging()
+    ic_newton.enable_logging()
+    minimizer = ift.NewtonCG(ic_newton, activate_logging=True)
 
     ## number of samples used to estimate the KL
     N_samples = 20
@@ -143,10 +145,15 @@ if __name__ == '__main__':
         plot.add([A2.force(KL.position),
                   A2.force(mock_position)],
                  title="power2")
-        plot.output(nx=2,
+        plot.add((ic_newton.history, ic_sampling.history,
+                  minimizer.inversion_history),
+                 label=['KL', 'Sampling', 'Newton inversion'],
+                 title='Cumulative energies', s=[None, None, 1],
+                 alpha=[None, 0.2, None])
+        plot.output(nx=3,
                     ny=2,
                     ysize=10,
-                    xsize=10,
+                    xsize=15,
                     name=filename.format("loop_{:02d}".format(i)))
 
     # Done, draw posterior samples

nifty6/__init__.py

@@ -25,7 +25,8 @@ from .operators.adder import Adder
 from .operators.diagonal_operator import DiagonalOperator
 from .operators.distributors import DOFDistributor, PowerDistributor
 from .operators.domain_tuple_field_inserter import DomainTupleFieldInserter
-from .operators.contraction_operator import ContractionOperator
+from .operators.einsum import LinearEinsum, MultiLinearEinsum
+from .operators.contraction_operator import ContractionOperator, IntegrationOperator
 from .operators.linear_interpolation import LinearInterpolator
 from .operators.endomorphic_operator import EndomorphicOperator
 from .operators.harmonic_operators import (
@@ -35,15 +36,16 @@ from .operators.field_zero_padder import FieldZeroPadder
 from .operators.inversion_enabler import InversionEnabler
 from .operators.mask_operator import MaskOperator
 from .operators.regridding_operator import RegriddingOperator
-from .operators.sampling_enabler import SamplingEnabler
+from .operators.sampling_enabler import SamplingEnabler, SamplingDtypeSetter
 from .operators.sandwich_operator import SandwichOperator
 from .operators.scaling_operator import ScalingOperator
+from .operators.selection_operators import SliceOperator, SplitOperator
 from .operators.block_diagonal_operator import BlockDiagonalOperator
 from .operators.outer_product_operator import OuterProduct
 from .operators.simple_linear_operators import (
-    VdotOperator, ConjugationOperator, Realizer,
-    FieldAdapter, ducktape, GeometryRemover, NullOperator,
-    MatrixProductOperator, PartialExtractor)
+    VdotOperator, ConjugationOperator, Realizer, FieldAdapter, ducktape,
+    GeometryRemover, NullOperator, PartialExtractor)
+from .operators.matrix_product_operator import MatrixProductOperator
 from .operators.value_inserter import ValueInserter
 from .operators.energy_operators import (
     EnergyOperator, GaussianEnergy, PoissonianEnergy, InverseGammaLikelihood,

nifty6/domains/rg_space.py

@@ -142,8 +142,7 @@ class RGSpace(StructuredDomain):
 
     @staticmethod
     def _kernel(x, sigma):
-        from ..sugar import exp
-        return exp(x*x * (-2.*np.pi*np.pi*sigma*sigma))
+        return (x*x * (-2.*np.pi*np.pi*sigma*sigma)).ptw("exp")
 
     def get_fft_smoothing_kernel_function(self, sigma):
         if (not self.harmonic):

nifty6/extra.py

@@ -42,8 +42,8 @@ def _adjoint_implementation(op, domain_dtype, target_dtype, atol, rtol,
     needed_cap = op.TIMES | op.ADJOINT_TIMES
     if (op.capability & needed_cap) != needed_cap:
         return
-    f1 = from_random("normal", op.domain, dtype=domain_dtype)
-    f2 = from_random("normal", op.target, dtype=target_dtype)
+    f1 = from_random(op.domain, "normal", dtype=domain_dtype)
+    f2 = from_random(op.target, "normal", dtype=target_dtype)
     res1 = f1.s_vdot(op.adjoint_times(f2))
     res2 = op.times(f1).s_vdot(f2)
     if only_r_linear:
@@ -55,11 +55,11 @@ def _inverse_implementation(op, domain_dtype, target_dtype, atol, rtol):
     needed_cap = op.TIMES | op.INVERSE_TIMES
     if (op.capability & needed_cap) != needed_cap:
         return
-    foo = from_random("normal", op.target, dtype=target_dtype)
+    foo = from_random(op.target, "normal", dtype=target_dtype)
     res = op(op.inverse_times(foo))
     assert_allclose(res, foo, atol=atol, rtol=rtol)
 
-    foo = from_random("normal", op.domain, dtype=domain_dtype)
+    foo = from_random(op.domain, "normal", dtype=domain_dtype)
     res = op.inverse_times(op(foo))
     assert_allclose(res, foo, atol=atol, rtol=rtol)
 
@@ -75,8 +75,8 @@ def _check_linearity(op, domain_dtype, atol, rtol):
     needed_cap = op.TIMES
     if (op.capability & needed_cap) != needed_cap:
         return
-    fld1 = from_random("normal", op.domain, dtype=domain_dtype)
-    fld2 = from_random("normal", op.domain, dtype=domain_dtype)
+    fld1 = from_random(op.domain, "normal", dtype=domain_dtype)
+    fld2 = from_random(op.domain, "normal", dtype=domain_dtype)
     alpha = np.random.random()  # FIXME: this can break badly with MPI!
     val1 = op(alpha*fld1+fld2)
     val2 = alpha*op(fld1)+op(fld2)
@@ -88,7 +88,7 @@ def _actual_domain_check_linear(op, domain_dtype=None, inp=None):
     if (op.capability & needed_cap) != needed_cap:
         return
     if domain_dtype is not None:
-        inp = from_random("normal", op.domain, dtype=domain_dtype)
+        inp = from_random(op.domain, "normal", dtype=domain_dtype)
     elif inp is None:
         raise ValueError('Need to specify either dtype or inp')
     assert_(inp.domain is op.domain)
@@ -219,7 +219,7 @@ def consistency_check(op, domain_dtype=np.float64, target_dtype=np.float64,
 def _get_acceptable_location(op, loc, lin):
     if not np.isfinite(lin.val.s_sum()):
         raise ValueError('Initial value must be finite')
-    dir = from_random("normal", loc.domain)
+    dir = from_random(loc.domain, "normal")
     dirder = lin.jac(dir)
     if dirder.norm() == 0:
         dir = dir * (lin.val.norm()*1e-5)
@@ -296,3 +296,11 @@ def check_jacobian_consistency(op, loc, tol=1e-8, ntries=100, perf_check=True):
             print(hist)
             raise ValueError("gradient and value seem inconsistent")
         loc = locnext
+
+        # FIXME The following code shows that we need prober tests for complex
+        # derivatives
+        ddtype = loc.values()[0].dtype if isinstance(loc, MultiField) else loc.dtype
+        tdtype = dirder.values()[0].dtype if isinstance(dirder, MultiField) else dirder.dtype
+        only_r_linear = ddtype != tdtype
+        consistency_check(linmid.jac, domain_dtype=ddtype, target_dtype=tdtype,
+                          only_r_linear=only_r_linear)

nifty6/field.py

@@ -20,9 +20,10 @@ import numpy as np
 
 from . import utilities
 from .domain_tuple import DomainTuple
+from .operators.operator import Operator
 
 
-class Field(object):
+class Field(Operator):
     """The discrete representation of a continuous field over multiple spaces.
 
     Stores data arrays and carries all the needed meta-information (i.e. the
@@ -49,7 +50,7 @@ class Field(object):
             raise TypeError("domain must be of type DomainTuple")
         if not isinstance(val, np.ndarray):
             if np.isscalar(val):
-                val = np.full(domain.shape, val)
+                val = np.broadcast_to(val, domain.shape)
             else:
                 raise TypeError("val must be of type numpy.ndarray")
         if domain.shape != val.shape:
@@ -123,7 +124,7 @@ class Field(object):
         return Field(DomainTuple.make(new_domain), self._val)
 
     @staticmethod
-    def from_random(random_type, domain, dtype=np.float64, **kwargs):
+    def from_random(domain, random_type='normal', dtype=np.float64, **kwargs):
         """Draws a random field with the given parameters.
 
         Parameters
@@ -282,7 +283,7 @@ class Field(object):
             raise TypeError("The multiplier must be an instance of " +
                             "the Field class")
         from .operators.outer_product_operator import OuterProduct
-        return OuterProduct(self, x.domain)(x)
+        return OuterProduct(x.domain, self)(x)
 
     def vdot(self, x, spaces=None):
         """Computes the dot product of 'self' with x.
@@ -634,10 +635,9 @@ class Field(object):
         Field
             The result of the operation.
         """
-        from .sugar import sqrt
         if self.scalar_weight(spaces) is not None:
             return self._contraction_helper('std', spaces)
-        return sqrt(self.var(spaces))
+        return self.var(spaces).ptw("sqrt")
 
     def s_std(self):
         """Determines the standard deviation of the Field.
@@ -677,17 +677,6 @@ class Field(object):
     def flexible_addsub(self, other, neg):
         return self-other if neg else self+other
 
-    def sigmoid(self):
-        return 0.5*(1.+self.tanh())
-
-    def clip(self, min=None, max=None):
-        min = min.val if isinstance(min, Field) else min
-        max = max.val if isinstance(max, Field) else max
-        return Field(self._domain, np.clip(self._val, min, max))
-
-    def one_over(self):
-        return 1/self
-
     def _binary_op(self, other, op):
         # if other is a field, make sure that the domains match
         f = getattr(self._val, op)
@@ -699,6 +688,26 @@ class Field(object):
             return Field(self._domain, f(other))
         return NotImplemented
 
+    def _prep_args(self, args, kwargs):
+        for arg in args + tuple(kwargs.values()):
+            if not (arg is None or np.isscalar(arg) or arg.jac is None):
+                raise TypeError("bad argument")
+        argstmp = tuple(arg if arg is None or np.isscalar(arg) else arg._val
+                        for arg in args)
+        kwargstmp = {key: val if val is None or np.isscalar(val) else val._val
+                     for key, val in kwargs.items()}
+        return argstmp, kwargstmp
+
+    def ptw(self, op, *args, **kwargs):
+        from .pointwise import ptw_dict
+        argstmp, kwargstmp = self._prep_args(args, kwargs)
+        return Field(self._domain, ptw_dict[op][0](self._val, *argstmp, **kwargstmp))
+
+    def ptw_with_deriv(self, op, *args, **kwargs):
+        from .pointwise import ptw_dict
+        argstmp, kwargstmp = self._prep_args(args, kwargs)
+        tmp = ptw_dict[op][1](self._val, *argstmp, **kwargstmp)
+        return (Field(self._domain, tmp[0]), Field(self._domain, tmp[1]))
 
 for op in ["__add__", "__radd__",
            "__sub__", "__rsub__",
@@ -721,11 +730,3 @@ for op in ["__iadd__", "__isub__", "__imul__", "__idiv__",
                 "In-place operations are deliberately not supported")
         return func2
     setattr(Field, op, func(op))
-
-for f in ["sqrt", "exp", "log", "sin", "cos", "tan", "sinh", "cosh", "tanh",
-          "absolute", "sinc", "sign", "log10", "log1p", "expm1"]:
-    def func(f):
-        def func2(self):
-            return Field(self._domain, getattr(np, f)(self.val))
-        return func2
-    setattr(Field, f, func(f))

nifty6/library/correlated_fields.py

@@ -61,7 +61,7 @@ def _lognormal_moments(mean, sig, N=0):
     if not np.all(sig > 0):
         raise ValueError("sig must be greater 0; got {!r}".format(sig))
 
-    logsig = np.sqrt(np.log((sig/mean)**2 + 1))
+    logsig = np.sqrt(np.log1p((sig/mean)**2))
     logmean = np.log(mean) - logsig**2/2
     return logmean, logsig
 
@@ -126,7 +126,7 @@ class _LognormalMomentMatching(Operator):
         logmean, logsig = _lognormal_moments(mean, sig, N_copies)
         self._mean = mean
         self._sig = sig
-        op = _normal(logmean, logsig, key, N_copies).exp()
+        op = _normal(logmean, logsig, key, N_copies).ptw("exp")
         self._domain, self._target = op.domain, op.target
         self.apply = op.apply
 
@@ -224,8 +224,8 @@ class _Normalization(Operator):
 
     def apply(self, x):
         self._check_input(x)
-        amp = x.exp()
-        spec = (2*x).exp()
+        amp = x.ptw("exp")
+        spec = amp**2
         # FIXME This normalizes also the zeromode which is supposed to be left
         # untouched by this operator
         return self._specsum(self._mode_multiplicity(spec))**(-0.5)*amp
@@ -243,12 +243,11 @@ class _SpecialSum(EndomorphicOperator):
 
 
 class _Distributor(LinearOperator):
-    def __init__(self, dofdex, domain, target, space=0):
+    def __init__(self, dofdex, domain, target):
         self._dofdex = dofdex
 
         self._target = makeDomain(target)
         self._domain = makeDomain(domain)
-        self._sl = (slice(None),)*space
         self._capability = self.TIMES | self.ADJOINT_TIMES
 
     def apply(self, x, mode):
@@ -282,7 +281,7 @@ class _Amplitude(Operator):
             distributed_tgt = makeDomain((UnstructuredDomain(len(dofdex)),
                                           target))
             target = makeDomain((UnstructuredDomain(N_copies), target))
-            Distributor = _Distributor(dofdex, target, distributed_tgt, 0)
+            Distributor = _Distributor(dofdex, target, distributed_tgt)
         else:
             N_copies = 0
             space = 0
@@ -332,17 +331,17 @@ class _Amplitude(Operator):
         sig_fluc = vol1 @ ps_expander @ fluctuations
 
         xi = ducktape(dom, None, key)
-        sigma = sig_flex*(Adder(shift) @ sig_asp).sqrt()
+        sigma = sig_flex*(Adder(shift) @ sig_asp).ptw("sqrt")
         smooth = _SlopeRemover(target, space) @ twolog @ (sigma*xi)
         op = _Normalization(target, space) @ (slope + smooth)
         if N_copies > 0:
             op = Distributor @ op