Merge remote-tracking branch 'origin/NIFTy_7' into constant_support_pa

ChangeLog.md

@@ -76,6 +76,19 @@ print(met)
 print(met.draw_sample())
 ```
 
+New approach for sampling complex numbers
+=========================================
+
+When calling draw_sample_with_dtype with a complex dtype,
+the variance is now used for the imaginary part and real part separately.
+This is done in order to be consistent with the Hamiltonian.
+Note that by this, 
+```
+np.std(ift.from_random(domain, 'normal', dtype=np.complex128).val)
+```` 
+does not give 1, but sqrt(2) as a result.
+
+
 MPI parallelisation over samples in MetricGaussianKL
 ----------------------------------------------------
 
@@ -91,6 +104,7 @@ the generation of reproducible random numbers in the presence of MPI parallelism
 and leads to cleaner code overall. Please see the documentation of
 `nifty7.random` for details.
 
+
 Interface Change for from_random and OuterProduct
 -------------------------------------------------
 

src/field.py

@@ -136,6 +136,8 @@ class Field(Operator):
             The domain of the output random Field.
         dtype : type
             The datatype of the output random Field.
+            If the datatype is complex, each real and imaginary part
+            have variance 1
 
         Returns
         -------

src/multi_field.py

@@ -112,6 +112,8 @@ class MultiField(Operator):
             The domain of the output random Field.
         dtype : type
             The datatype of the output random Field.
+            If the datatype is complex, each real an imaginary part have
+            variance 1.
 
         Returns
         -------

src/operators/energy_operators.py

@@ -236,6 +236,13 @@ class GaussianEnergy(EnergyOperator):
     domain : Domain, DomainTuple, tuple of Domain or MultiDomain
         Operator domain. By default it is inferred from `mean` or
         `covariance` if specified
+    sampling_dtype : type
+        Here one can specify whether the distribution is a complex Gaussian or
+        not. Note that for a complex Gaussian the inverse_covariance is
+        .. math ::
+        (<ff^dagger>)^{-1}_P(f)/2,
+        where the additional factor of 2 is necessary because the 
+        domain of s has double as many dimensions as in the real case.
 
     Note
     ----

src/random.py

@@ -238,8 +238,8 @@ class Random(object):
             raise TypeError("mean must not be complex for a real result field")
         if np.issubdtype(dtype, np.complexfloating):
             x = np.empty(shape, dtype=dtype)
-            x.real = _rng[-1].normal(mean.real, std*np.sqrt(0.5), shape)
-            x.imag = _rng[-1].normal(mean.imag, std*np.sqrt(0.5), shape)
+            x.real = _rng[-1].normal(mean.real, std, shape)
+            x.imag = _rng[-1].normal(mean.imag, std, shape)
         else:
             x = _rng[-1].normal(mean, std, shape).astype(dtype, copy=False)
         return x

src/sugar.py

@@ -266,6 +266,8 @@ def from_random(domain, random_type='normal', dtype=np.float64, **kwargs):
             The random distribution to use.
     dtype : type
         data type of the output field (e.g. numpy.float64)
+        If the datatype is complex, each real an imaginary part have
+        variance 1.
     **kwargs : additional parameters for the random distribution
         ('mean' and 'std' for 'normal', 'low' and 'high' for 'uniform')
 

test/test_operators/test_complex_consistency.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
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+
+import numpy as np
+
+import nifty7 as ift
+
+Nsamp = 20000
+np.random.seed(42)
+
+
+def _to_array(d):
+    if isinstance(d, np.ndarray):
+        return d
+    assert isinstance(d, dict)
+    return np.concatenate(list(d.values()))
+
+
+def test_GaussianEnergy():
+    sp = ift.UnstructuredDomain(Nsamp)
+    S = ift.ScalingOperator(sp, 1.)
+    samp = S.draw_sample_with_dtype(dtype=np.complex128)
+    real_std = np.std(samp.val.real)
+    imag_std = np.std(samp.val.imag)
+    np.testing.assert_allclose(real_std, imag_std,
+                               atol=5./np.sqrt(Nsamp))
+    sp1 = ift.UnstructuredDomain(1)
+    mean = ift.full(sp1, 0.)
+    real_icov = ift.ScalingOperator(sp1, real_std**(-2))
+    imag_icov = ift.ScalingOperator(sp1, imag_std**(-2))
+    real_energy = ift.GaussianEnergy(mean, inverse_covariance=real_icov)
+    imag_energy = ift.GaussianEnergy(mean, inverse_covariance=imag_icov)
+    icov = ift.ScalingOperator(sp1, 1.)
+    complex_energy = ift.GaussianEnergy(mean+0.j, inverse_covariance=icov)
+    for i in range(min(10, Nsamp)):
+        fld = ift.full(sp1, samp.val[i])
+        val1 = (real_energy(fld.real) + imag_energy(fld.imag)).val
+        val2 = complex_energy(fld).val
+        np.testing.assert_allclose(val1, val2, atol=10./np.sqrt(Nsamp))

test/test_operators/test_fisher_metric.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
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+
+import numpy as np
+import pytest
+
+import nifty7 as ift
+
+from ..common import list2fixture, setup_function, teardown_function
+
+spaces = [ift.GLSpace(5),
+          ift.MultiDomain.make({'': ift.RGSpace(5, distances=.789)}),
+          (ift.RGSpace(3, distances=.789), ift.UnstructuredDomain(2))]
+pmp = pytest.mark.parametrize
+field = list2fixture([ift.from_random(sp, 'normal') for sp in spaces] +
+                     [ift.from_random(sp, 'normal', dtype=np.complex128) for sp in spaces])
+
+dtype = list2fixture([np.float64,
+                     np.complex128])
+Nsamp = 2000
+np.random.seed(42)
+
+
+def _to_array(d):
+    if isinstance(d, np.ndarray):
+        return d
+    assert isinstance(d, dict)
+    return np.concatenate(list(d.values()))
+
+
+def _complex2real(sp):
+    tup = tuple([d for d in sp])
+    rsp = ift.DomainTuple.make((ift.UnstructuredDomain(2),) + tup)
+    rl = ift.DomainTupleFieldInserter(rsp, 0, (0,))
+    im = ift.DomainTupleFieldInserter(rsp, 0, (1,))
+    x = ift.ScalingOperator(sp, 1)
+    return rl(x.real)+im(x.imag)
+
+
+def test_complex2real():
+    sp = ift.UnstructuredDomain(3)
+    op = _complex2real(ift.makeDomain(sp))
+    f = ift.from_random(op.domain, 'normal', dtype=np.complex128)
+    assert np.all((f == op.adjoint_times(op(f))).val)
+    assert op(f).dtype == np.float64
+    f = ift.from_random(op.target, 'normal')
+    assert np.all((f == op(op.adjoint_times(f))).val)
+
+
+def energy_tester(pos, get_noisy_data, energy_initializer):
+    if isinstance(pos, ift.Field):
+        if np.issubdtype(pos.dtype, np.complexfloating):
+            op = _complex2real(pos.domain)
+        else:
+            op = ift.ScalingOperator(pos.domain, 1.)
+    else:
+        ops = []
+        for k,dom in pos.domain.items():
+            if np.issubdtype(pos[k].dtype, np.complexfloating):
+                ops.append(_complex2real(dom).ducktape(k).ducktape_left(k))
+            else:
+                FA = ift.FieldAdapter(dom, k)
+                ops.append(FA.adjoint @ FA)
+        realizer = ift.utilities.my_sum(ops)
+        from nifty7.operator_spectrum import _DomRemover
+        flattener = _DomRemover(realizer.target)
+        op = flattener @ realizer
+
+    npos = op(pos)
+    nget_noisy_data = lambda mean: get_noisy_data(op.adjoint_times(mean))
+    nenergy_initializer = lambda mean: energy_initializer(mean) @ op.adjoint
+    _actual_energy_tester(npos, nget_noisy_data, nenergy_initializer)
+
+
+def _actual_energy_tester(pos, get_noisy_data, energy_initializer):
+    domain = pos.domain
+    test_vec = ift.from_random(domain, 'normal')
+    results = []
+    lin = ift.Linearization.make_var(pos)
+    for i in range(Nsamp):
+        data = get_noisy_data(pos)
+        energy = energy_initializer(data)
+        grad = energy(lin).jac.adjoint(ift.full(energy.target, 1.))
+        results.append(_to_array((grad*grad.s_vdot(test_vec)).val))
+    print(energy)
+    print(grad)
+    res = np.mean(np.array(results), axis=0)
+    std = np.std(np.array(results), axis=0)/np.sqrt(Nsamp)
+    energy = energy_initializer(data)
+    lin = ift.Linearization.make_var(pos, want_metric=True)
+    res2 = _to_array(energy(lin).metric(test_vec).val)
+    np.testing.assert_allclose(res/std, res2/std, atol=6)
+
+
+def test_GaussianEnergy(field):
+    dtype = field.dtype
+    icov = ift.from_random(field.domain, 'normal')**2
+    icov = ift.makeOp(icov)
+    get_noisy_data = lambda mean: mean + icov.draw_sample_with_dtype(
+        from_inverse=True, dtype=dtype)
+    E_init = lambda data: ift.GaussianEnergy(mean=data, inverse_covariance=icov)
+    energy_tester(field, get_noisy_data, E_init)
+
+
+def test_PoissonEnergy(field):
+    if not isinstance(field, ift.Field):
+        pytest.skip("MultiField Poisson energy  not supported")
+    if np.iscomplexobj(field.val):
+        pytest.skip("Poisson energy not defined for complex flux")
+    get_noisy_data = lambda mean: ift.makeField(mean.domain, np.random.poisson(mean.val))
+    # Make rate positive and high enough to avoid bad statistic
+    lam = 10*(field**2).clip(0.1, None)
+    E_init = lambda data: ift.PoissonianEnergy(data)
+    energy_tester(lam, get_noisy_data, E_init)
+
+def test_VariableCovarianceGaussianEnergy(dtype):
+    dom = ift.UnstructuredDomain(3)
+    res = ift.from_random(dom, 'normal', dtype=dtype)
+    ivar = ift.from_random(dom, 'normal')**2+4.
+    mf = ift.MultiField.from_dict({'res':res, 'ivar':ivar})
+    energy = ift.VariableCovarianceGaussianEnergy(dom, 'res', 'ivar', dtype)
+    def get_noisy_data(mean):
+        samp = ift.from_random(dom, 'normal', dtype)
+        samp = samp/mean['ivar'].sqrt()
+        return samp + mean['res']
+    def E_init(data):
+        adder = ift.Adder(ift.MultiField.from_dict({'res':data}), neg=True)
+        return energy.partial_insert(adder)
+    energy_tester(mf, get_noisy_data, E_init)