diff --git a/ChangeLog b/ChangeLog
index f9be1d69d1228bfb28b7615db495a24f83def34b..06ae02cfecdfe5f2451619755a7c5e5be79e7412 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -58,6 +58,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
====================================================
@@ -73,6 +86,7 @@ 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
=================================================
diff --git a/nifty6/field.py b/nifty6/field.py
index 4c94d58d2eba2e1e4ea288bc92fd66249b8c496b..1cd9253841ba80126bd49e2c0e956316e1e6d123 100644
--- a/nifty6/field.py
+++ b/nifty6/field.py
@@ -135,6 +135,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
-------
diff --git a/nifty6/multi_field.py b/nifty6/multi_field.py
index 1051aaec7cd49a96e46e4f00c58c4257c81dd7df..e880fa9f84945befe66743fb6124e106e711c008 100644
--- a/nifty6/multi_field.py
+++ b/nifty6/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
-------
diff --git a/nifty6/operators/energy_operators.py b/nifty6/operators/energy_operators.py
index ff05301c0072d8729945d5bbc977d48461f6fd11..21ece3648ad781d20b154722d6906baf92e0068f 100644
--- a/nifty6/operators/energy_operators.py
+++ b/nifty6/operators/energy_operators.py
@@ -27,7 +27,7 @@ from .linear_operator import LinearOperator
from .operator import Operator
from .sampling_enabler import SamplingDtypeSetter, SamplingEnabler
from .scaling_operator import ScalingOperator
-from .simple_linear_operators import VdotOperator
+from .simple_linear_operators import VdotOperator, FieldAdapter
def _check_sampling_dtype(domain, dtypes):
@@ -187,6 +187,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
----
diff --git a/nifty6/random.py b/nifty6/random.py
index 85278324a02e7eb1472da68a30ba20aecf9236d5..82b9d718f42c774be73c37da0d53eed3a19ab470 100644
--- a/nifty6/random.py
+++ b/nifty6/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
diff --git a/nifty6/sugar.py b/nifty6/sugar.py
index 683b970af363b757bbceadd8ec4c5eb140705b2a..d384af8598828dd7fc27132a9f9e8bb2d4c418bc 100644
--- a/nifty6/sugar.py
+++ b/nifty6/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')
diff --git a/test/test_operators/test_complex_consistency.py b/test/test_operators/test_complex_consistency.py
new file mode 100644
index 0000000000000000000000000000000000000000..dbb34fb303a86ed10f493d93022a99011eb71078
--- /dev/null
+++ b/test/test_operators/test_complex_consistency.py
@@ -0,0 +1,77 @@
+# 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 nifty6 as ift
+
+from ..common import list2fixture, setup_function, teardown_function
+from scipy.stats import norm
+
+
+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))
+
+
+def test_WienerFilter():
+ sp = ift.UnstructuredDomain(30)
+ S = ift.ScalingOperator(sp, 1.)
+ R = ift.ContractionOperator(sp, None)
+ N = ift.ScalingOperator(R.target, 1/5**2)
+ IC = ift.GradientNormController(iteration_limit=10, tol_abs_gradnorm=1e-7)
+ D = ift.WienerFilterCurvature(R, N, S, IC, IC).inverse
+ s = S.draw_sample_with_dtype(dtype=np.complex128)
+ d = R(s)+N.draw_sample_with_dtype(dtype=np.complex128)
+ j = R.adjoint_times(N.inverse_times(d))
+ m = D(j)
+ sr = s.real
+ jr = j.real
+ mr = D(jr)
+ si = s.imag
+ ji = j.imag
+ mi = D(ji)
+ np.testing.assert_allclose((mr+mi*1.j).val, m.val, rtol=1e-7)
diff --git a/test/test_operators/test_fisher_metric.py b/test/test_operators/test_fisher_metric.py
new file mode 100644
index 0000000000000000000000000000000000000000..bacb89244f0418a9888929e04e52aaaf400fd435
--- /dev/null
+++ b/test/test_operators/test_fisher_metric.py
@@ -0,0 +1,107 @@
+# 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 nifty6 as ift
+
+from ..common import list2fixture, setup_function, teardown_function
+
+spaces = [ift.GLSpace(5),
+ (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])
+
+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_complex(pos, get_noisy_data, energy_initializer):
+ op = _complex2real(pos.domain)
+ npos = op(pos)
+ nget_noisy_data = lambda mean : get_noisy_data(op.adjoint_times(mean))
+ nenergy_initializer = lambda mean : energy_initializer(mean) @ op.adjoint
+ energy_tester(npos, nget_noisy_data, nenergy_initializer)
+
+def energy_tester(pos, get_noisy_data, energy_initializer):
+ if np.issubdtype(pos.dtype, np.complexfloating):
+ energy_tester_complex(pos, get_noisy_data, energy_initializer)
+ return
+ 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 mean : ift.GaussianEnergy(mean=mean,
+ inverse_covariance=icov)
+ energy_tester(field, get_noisy_data, E_init)
+
+def test_PoissonEnergy(field):
+ if not isinstance(field, ift.Field):
+ return
+ if np.iscomplexobj(field.val):
+ return
+ def get_noisy_data(mean):
+ return ift.makeField(mean.domain, np.random.poisson(mean.val))
+ lam = 10*(field**2).clip(0.1,None) # make rate positive and high enough to avoid bad statistic
+ E_init = lambda mean : ift.PoissonianEnergy(mean)
+ energy_tester(lam, get_noisy_data, E_init)
+