From 3879673330f9b6031ca0cce174f4080d9c7fdcfe Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Tue, 12 May 2020 13:09:52 +0200
Subject: [PATCH] Introduce SamplingDtypeEnabler
---
demos/getting_started_0.ipynb | 21 ++--
demos/getting_started_1.py | 4 +-
demos/getting_started_3.py | 2 +-
demos/getting_started_mf.py | 2 +-
nifty6/library/wiener_filter_curvature.py | 23 +++-
nifty6/linearization.py | 2 -
nifty6/minimization/metric_gaussian_kl.py | 26 ++---
nifty6/multi_field.py | 19 ++--
nifty6/operators/block_diagonal_operator.py | 8 +-
nifty6/operators/diagonal_operator.py | 4 +-
nifty6/operators/endomorphic_operator.py | 9 +-
nifty6/operators/energy_operators.py | 119 ++++++++++++++++++--
nifty6/operators/inversion_enabler.py | 11 +-
nifty6/operators/linear_operator.py | 4 +-
nifty6/operators/operator_adapter.py | 15 ++-
nifty6/operators/sampling_enabler.py | 13 ++-
nifty6/operators/sandwich_operator.py | 6 +-
nifty6/operators/scaling_operator.py | 4 +-
nifty6/operators/sum_operator.py | 4 +-
nifty6/probing.py | 4 +-
test/test_energy_gradients.py | 16 +--
test/test_field.py | 4 +-
test/test_gaussian_energy.py | 2 +-
test/test_kl.py | 3 +-
test/test_minimizers.py | 16 ++-
test/test_mpi/test_kl.py | 5 +-
test/test_operators/test_interpolated.py | 3 +-
test/test_operators/test_jacobian.py | 17 +--
test/test_sugar.py | 2 +-
29 files changed, 237 insertions(+), 131 deletions(-)
diff --git a/demos/getting_started_0.ipynb b/demos/getting_started_0.ipynb
index ff2cddb70..cce66e1f2 100644
--- a/demos/getting_started_0.ipynb
+++ b/demos/getting_started_0.ipynb
@@ -171,7 +171,7 @@
" tol_abs_gradnorm=0.1)\n",
" # WienerFilterCurvature is (R.adjoint*N.inverse*R + Sh.inverse) plus some handy\n",
" # helper methods.\n",
- " return ift.WienerFilterCurvature(R,N,Sh,iteration_controller=IC,iteration_controller_sampling=IC)"
+ " return ift.WienerFilterCurvature(R,N,Sh,iteration_controller=IC,iteration_controller_sampling=IC, data_sampling_dtype=np.float64, prior_sampling_dtype=np.float64)"
]
},
{
@@ -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(dtype=np.float64)\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(dtype=np.float64)\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)"
@@ -571,7 +571,7 @@
"N = ift.ScalingOperator(s_space, sigma2)\n",
"\n",
"# Fields and data\n",
- "sh = Sh.draw_sample(dtype=np.float64)\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",
@@ -652,7 +652,7 @@
"ma = np.max(s_data)\n",
"\n",
"fig, axes = plt.subplots(3, 2, figsize=(15, 22.5))\n",
- "sample = HT(curv.draw_sample(dtype=np.float64, from_inverse=True)+m).val\n",
+ "sample = HT(curv.draw_sample(from_inverse=True)+m).val\n",
"post_mean = (m_mean + HT(m)).val\n",
"\n",
"data = [s_data, m_data, post_mean, sample, s_data - m_data, uncertainty]\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,
diff --git a/demos/getting_started_1.py b/demos/getting_started_1.py
index 22d7d5c80..c00bbe743 100644
--- a/demos/getting_started_1.py
+++ b/demos/getting_started_1.py
@@ -114,8 +114,8 @@ if __name__ == '__main__':
N = ift.ScalingOperator(data_space, noise)
# Create mock data
- MOCK_SIGNAL = S.draw_sample(dtype=np.float64)
- MOCK_NOISE = N.draw_sample(dtype=np.float64)
+ 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
diff --git a/demos/getting_started_3.py b/demos/getting_started_3.py
index 7fb3f716c..1d5ec727f 100644
--- a/demos/getting_started_3.py
+++ b/demos/getting_started_3.py
@@ -99,7 +99,7 @@ if __name__ == '__main__':
# Generate mock signal and data
mock_position = ift.from_random('normal', signal_response.domain)
- data = signal_response(mock_position) + N.draw_sample(dtype=np.float64)
+ data = signal_response(mock_position) + N.draw_sample_with_dtype(dtype=np.float64)
# Minimization parameters
ic_sampling = ift.AbsDeltaEnergyController(
diff --git a/demos/getting_started_mf.py b/demos/getting_started_mf.py
index 68e9b78c8..c66467cbd 100644
--- a/demos/getting_started_mf.py
+++ b/demos/getting_started_mf.py
@@ -98,7 +98,7 @@ if __name__ == '__main__':
# Generate mock signal and data
mock_position = ift.from_random('normal', signal_response.domain)
- data = signal_response(mock_position) + N.draw_sample(dtype=np.float64)
+ data = signal_response(mock_position) + N.draw_sample_with_dtype(dtype=np.float64)
plot = ift.Plot()
plot.add(signal(mock_position), title='Ground Truth')
diff --git a/nifty6/library/wiener_filter_curvature.py b/nifty6/library/wiener_filter_curvature.py
index 173e1b0ce..3d0e6eabc 100644
--- a/nifty6/library/wiener_filter_curvature.py
+++ b/nifty6/library/wiener_filter_curvature.py
@@ -21,7 +21,9 @@ from ..operators.sandwich_operator import SandwichOperator
def WienerFilterCurvature(R, N, S, iteration_controller=None,
- iteration_controller_sampling=None):
+ iteration_controller_sampling=None,
+ data_sampling_dtype=None,
+ prior_sampling_dtype=None):
"""The curvature of the WienerFilterEnergy.
This operator implements the second derivative of the
@@ -43,10 +45,19 @@ def WienerFilterCurvature(R, N, S, iteration_controller=None,
iteration_controller_sampling : IterationController
The iteration controller to use for sampling.
"""
- M = SandwichOperator.make(R, N.inverse)
+ Ninv = N.inverse
+ Sinv = S.inverse
+ if data_sampling_dtype is not None:
+ from ..operators.energy_operators import SamplingDtypeEnabler
+ Ninv = SamplingDtypeEnabler(Ninv, data_sampling_dtype)
+ if prior_sampling_dtype is not None:
+ from ..operators.energy_operators import SamplingDtypeEnabler
+ Sinv = SamplingDtypeEnabler(Sinv, data_sampling_dtype)
+ M = SandwichOperator.make(R, Ninv)
if iteration_controller_sampling is not None:
- op = SamplingEnabler(M, S.inverse, iteration_controller_sampling,
- S.inverse)
+ op = SamplingEnabler(M, Sinv, iteration_controller_sampling,
+ Sinv)
else:
- op = M + S.inverse
- return InversionEnabler(op, iteration_controller, S.inverse)
+ op = M + Sinv
+ op = InversionEnabler(op, iteration_controller, Sinv)
+ return op
diff --git a/nifty6/linearization.py b/nifty6/linearization.py
index 83156aa23..abd63c2bd 100644
--- a/nifty6/linearization.py
+++ b/nifty6/linearization.py
@@ -18,7 +18,6 @@
import numpy as np
from .sugar import makeOp
-from . import utilities
from .operators.operator import Operator
@@ -277,7 +276,6 @@ class Linearization(Operator):
ContractionOperator(self._jac.target, spaces, 1)(self._jac))
def ptw(self, op, *args, **kwargs):
- from .pointwise import ptw_dict
t1, t2 = self._val.ptw_with_deriv(op, *args, **kwargs)
return self.new(t1, makeOp(t2)(self._jac))
diff --git a/nifty6/minimization/metric_gaussian_kl.py b/nifty6/minimization/metric_gaussian_kl.py
index 6c30ac087..7872b462c 100644
--- a/nifty6/minimization/metric_gaussian_kl.py
+++ b/nifty6/minimization/metric_gaussian_kl.py
@@ -67,8 +67,8 @@ class _KLMetric(EndomorphicOperator):
self._check_input(x, mode)
return self._KL.apply_metric(x)
- def draw_sample(self, dtype, from_inverse=False):
- return self._KL._metric_sample(dtype, from_inverse)
+ def draw_sample(self, from_inverse=False):
+ return self._KL._metric_sample(from_inverse)
class MetricGaussianKL(Energy):
@@ -114,13 +114,6 @@ class MetricGaussianKL(Energy):
If not None, samples will be distributed as evenly as possible
across this communicator. If `mirror_samples` is set, then a sample and
its mirror image will always reside on the same task.
- lh_sampling_dtype : type
- Determines which dtype in data space shall be used for drawing samples
- from the metric. If the inference is based on complex data,
- lh_sampling_dtype shall be set to complex accordingly. The reason for
- the presence of this parameter is that metric of the likelihood energy
- is just an `Operator` which does not know anything about the dtype of
- the fields on which it acts. Default is float64.
_local_samples : None
Only a parameter for internal uses. Typically not to be set by users.
@@ -138,8 +131,7 @@ class MetricGaussianKL(Energy):
def __init__(self, mean, hamiltonian, n_samples, constants=[],
point_estimates=[], mirror_samples=False,
- napprox=0, comm=None, _local_samples=None,
- lh_sampling_dtype=np.float64):
+ napprox=0, comm=None, _local_samples=None):
super(MetricGaussianKL, self).__init__(mean)
if not isinstance(hamiltonian, StandardHamiltonian):
@@ -179,8 +171,7 @@ class MetricGaussianKL(Energy):
sseq = random.spawn_sseq(self._n_samples)
for i in range(self._lo, self._hi):
with random.Context(sseq[i]):
- _local_samples.append(met.draw_sample(
- dtype=lh_sampling_dtype, from_inverse=True))
+ _local_samples.append(met.draw_sample(from_inverse=True))
_local_samples = tuple(_local_samples)
else:
if len(_local_samples) != self._hi-self._lo:
@@ -206,13 +197,12 @@ class MetricGaussianKL(Energy):
self._val = _np_allreduce_sum(self._comm, v)[()] / self._n_eff_samples
self._grad = _allreduce_sum_field(self._comm, g) / self._n_eff_samples
self._metric = None
- self._sampdt = lh_sampling_dtype
def at(self, position):
return MetricGaussianKL(
position, self._hamiltonian, self._n_samples, self._constants,
self._point_estimates, self._mirror_samples, comm=self._comm,
- _local_samples=self._local_samples, lh_sampling_dtype=self._sampdt)
+ _local_samples=self._local_samples)
@property
def value(self):
@@ -264,7 +254,7 @@ class MetricGaussianKL(Energy):
if self._mirror_samples:
yield -s
- def _metric_sample(self, dtype, from_inverse=False):
+ def _metric_sample(self, from_inverse=False):
if from_inverse:
raise NotImplementedError()
lin = self._lin.with_want_metric()
@@ -272,7 +262,7 @@ class MetricGaussianKL(Energy):
sseq = random.spawn_sseq(self._n_samples)
for i, v in enumerate(self._local_samples):
with random.Context(sseq[self._lo+i]):
- samp = samp + self._hamiltonian(lin+v).metric.draw_sample(dtype=dtype, from_inverse=False)
+ samp = samp + self._hamiltonian(lin+v).metric.draw_sample(from_inverse=False)
if self._mirror_samples:
- samp = samp + self._hamiltonian(lin-v).metric.draw_sample(dtype=dtype, from_inverse=False)
+ samp = samp + self._hamiltonian(lin-v).metric.draw_sample(from_inverse=False)
return _allreduce_sum_field(self._comm, samp)/self._n_eff_samples
diff --git a/nifty6/multi_field.py b/nifty6/multi_field.py
index 341234904..d6a5aadbf 100644
--- a/nifty6/multi_field.py
+++ b/nifty6/multi_field.py
@@ -11,7 +11,7 @@
# 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-2019 Max-Planck-Society
+# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
@@ -49,15 +49,15 @@ class MultiField(Operator):
self._val = val
@staticmethod
- def from_dict(dict, domain=None):
+ def from_dict(dct, domain=None):
if domain is None:
- for dd in dict.values():
+ for dd in dct.values():
if not isinstance(dd.domain, DomainTuple):
raise TypeError('Values of dictionary need to be Fields '
'defined on DomainTuples.')
domain = MultiDomain.make({key: v._domain
- for key, v in dict.items()})
- res = tuple(dict[key] if key in dict else Field(dom, 0.)
+ for key, v in dct.items()})
+ res = tuple(dct[key] if key in dct else Field(dom, 0.)
for key, dom in zip(domain.keys(), domain.domains()))
return MultiField(domain, res)
@@ -103,10 +103,11 @@ class MultiField(Operator):
@staticmethod
def from_random(random_type, domain, dtype=np.float64, **kwargs):
domain = MultiDomain.make(domain)
-# dtype = MultiField.build_dtype(dtype, domain)
- return MultiField(
- domain, tuple(Field.from_random(random_type, dom, dtype, **kwargs)
- for dom in domain._domains))
+ if dtype in [np.float64, np.complex128]:
+ dtype = {kk: dtype for kk in domain.keys()}
+ dct = {kk: Field.from_random(random_type, domain[kk], dtype[kk], **kwargs)
+ for kk in domain.keys()}
+ return MultiField.from_dict(dct)
def _check_domain(self, other):
if other._domain != self._domain:
diff --git a/nifty6/operators/block_diagonal_operator.py b/nifty6/operators/block_diagonal_operator.py
index 4b4e94a1c..d2fb275e1 100644
--- a/nifty6/operators/block_diagonal_operator.py
+++ b/nifty6/operators/block_diagonal_operator.py
@@ -11,12 +11,10 @@
# 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-2019 Max-Planck-Society
+# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
-import numpy as np
-
from ..multi_domain import MultiDomain
from ..multi_field import MultiField
from .endomorphic_operator import EndomorphicOperator
@@ -48,10 +46,10 @@ class BlockDiagonalOperator(EndomorphicOperator):
for op, v in zip(self._ops, x.values()))
return MultiField(self._domain, val)
- def draw_sample(self, dtype, from_inverse=False):
+ def draw_sample_with_dtype(self, dtype, from_inverse=False):
from ..sugar import from_random
val = tuple(
- op.draw_sample(dtype, from_inverse)
+ op.draw_sample_with_dtype(dtype[key], from_inverse)
if op is not None
else from_random('normal', self._domain[key], dtype=dtype)
for op, key in zip(self._ops, self._domain.keys()))
diff --git a/nifty6/operators/diagonal_operator.py b/nifty6/operators/diagonal_operator.py
index 69ece97a6..21055550d 100644
--- a/nifty6/operators/diagonal_operator.py
+++ b/nifty6/operators/diagonal_operator.py
@@ -11,7 +11,7 @@
# 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-2019 Max-Planck-Society
+# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
@@ -160,7 +160,7 @@ class DiagonalOperator(EndomorphicOperator):
res = samp.val*np.sqrt(self._ldiag)
return Field(self._domain, res)
- def draw_sample(self, dtype, from_inverse=False):
+ def draw_sample_with_dtype(self, dtype, from_inverse=False):
res = Field.from_random(random_type="normal", domain=self._domain,
dtype=dtype)
return self.process_sample(res, from_inverse)
diff --git a/nifty6/operators/endomorphic_operator.py b/nifty6/operators/endomorphic_operator.py
index 2c423c6df..0f68e96b1 100644
--- a/nifty6/operators/endomorphic_operator.py
+++ b/nifty6/operators/endomorphic_operator.py
@@ -11,12 +11,10 @@
# 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-2019 Max-Planck-Society
+# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
-import numpy as np
-
from .linear_operator import LinearOperator
@@ -32,15 +30,16 @@ class EndomorphicOperator(LinearOperator):
for endomorphic operators."""
return self._domain
- def draw_sample(self, dtype, from_inverse=False):
+ def draw_sample_with_dtype(self, dtype, from_inverse=False):
"""Generate a zero-mean sample
+ FIXME
Generates a sample from a Gaussian distribution with zero mean and
covariance given by the operator.
Parameters
----------
- dtype : numpy datatype
+ dtype : numpy datatype FIXME
the data type to be used for the sample
from_inverse : bool (default : False)
if True, the sample is drawn from the inverse of the operator
diff --git a/nifty6/operators/energy_operators.py b/nifty6/operators/energy_operators.py
index 52eff6b3d..def9e7d0e 100644
--- a/nifty6/operators/energy_operators.py
+++ b/nifty6/operators/energy_operators.py
@@ -11,7 +11,7 @@
# 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-2019 Max-Planck-Society
+# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
@@ -28,6 +28,73 @@ from .operator import Operator
from .sampling_enabler import SamplingEnabler
from .scaling_operator import ScalingOperator
from .simple_linear_operators import VdotOperator
+from .endomorphic_operator import EndomorphicOperator
+
+
+def _check_sampling_dtype(domain, dtypes):
+ if dtypes is None:
+ return
+ if isinstance(domain, DomainTuple):
+ dtypes = {'': dtypes}
+ elif isinstance(domain, MultiDomain):
+ if dtypes in [np.float64, np.complex128]:
+ return
+ dtypes = dtypes.values()
+ if set(domain.keys()) != set(dtypes.keys()):
+ raise ValueError
+ else:
+ raise TypeError
+ for dt in dtypes.values():
+ if dt not in [np.float64, np.complex128]:
+ raise ValueError
+
+
+def _field_to_dtype(field):
+ if isinstance(field, Field):
+ dt = field.dtype
+ elif isinstance(field, MultiField):
+ dt = {kk: ff.dtype for kk, ff in field.items()}
+ else:
+ raise TypeError
+ _check_sampling_dtype(field.domain, dt)
+ return dt
+
+
+class SamplingDtypeEnabler(EndomorphicOperator):
+ def __init__(self, endomorphic_operator, dtype):
+ if not isinstance(endomorphic_operator, EndomorphicOperator):
+ raise TypeError
+ if not hasattr(endomorphic_operator, 'draw_sample_with_dtype'):
+ raise TypeError
+ dom = endomorphic_operator.domain
+ if isinstance(dom, MultiDomain):
+ if dtype in [np.float64, np.complex128]:
+ dtype = {kk: dtype for kk in dom.keys()}
+ if set(dtype.keys()) != set(dom.keys()):
+ raise TypeError
+ self._dtype = dtype
+ self._domain = dom
+ self._capability = endomorphic_operator._capability
+ self.apply = endomorphic_operator.apply
+ self._op = endomorphic_operator
+
+ def draw_sample(self, from_inverse=False):
+ """Generate a zero-mean sample
+
+ Generates a sample from a Gaussian distribution with zero mean and
+ covariance given by the operator.
+
+ Parameters
+ ----------
+ from_inverse : bool (default : False)
+ if True, the sample is drawn from the inverse of the operator
+
+ Returns
+ -------
+ Field
+ A sample from the Gaussian of given covariance.
+ """
+ return self._op.draw_sample_with_dtype(self._dtype, from_inverse=from_inverse)
class EnergyOperator(Operator):
@@ -117,11 +184,13 @@ class VariableCovarianceGaussianEnergy(EnergyOperator):
Inverse covariance diagonal key of the Gaussian.
"""
- def __init__(self, domain, residual_key, inverse_covariance_key):
+ def __init__(self, domain, residual_key, inverse_covariance_key, sampling_dtype):
self._r = str(residual_key)
self._icov = str(inverse_covariance_key)
dom = DomainTuple.make(domain)
self._domain = MultiDomain.make({self._r: dom, self._icov: dom})
+ self._sampling_dtype = sampling_dtype
+ _check_sampling_dtype(self._domain, sampling_dtype)
def apply(self, x):
self._check_input(x)
@@ -129,7 +198,8 @@ class VariableCovarianceGaussianEnergy(EnergyOperator):
if not x.want_metric:
return res
mf = {self._r: x.val[self._icov], self._icov: .5*x.val[self._icov]**(-2)}
- return res.add_metric(makeOp(MultiField.from_dict(mf)))
+ met = makeOp(MultiField.from_dict(mf))
+ return res.add_metric(SamplingDtypeEnabler(met, self._sampling_dtype))
class GaussianEnergy(EnergyOperator):
@@ -158,7 +228,7 @@ class GaussianEnergy(EnergyOperator):
At least one of the arguments has to be provided.
"""
- def __init__(self, mean=None, inverse_covariance=None, domain=None):
+ def __init__(self, mean=None, inverse_covariance=None, domain=None, sampling_dtype=None):
if mean is not None and not isinstance(mean, (Field, MultiField)):
raise TypeError
if inverse_covariance is not None and not isinstance(inverse_covariance, LinearOperator):
@@ -174,12 +244,25 @@ class GaussianEnergy(EnergyOperator):
if self._domain is None:
raise ValueError("no domain given")
self._mean = mean
+
+ # Infer sampling dtype
+ if self._mean is None:
+ _check_sampling_dtype(self._domain, sampling_dtype)
+ else:
+ if sampling_dtype is None:
+ sampling_dtype = _field_to_dtype(self._mean)
+ else:
+ if sampling_dtype != _field_to_dtype(self._mean):
+ raise ValueError("Sampling dtype and mean not compatible")
+
if inverse_covariance is None:
self._op = Squared2NormOperator(self._domain).scale(0.5)
self._met = ScalingOperator(self._domain, 1)
else:
self._op = QuadraticFormOperator(inverse_covariance)
self._met = inverse_covariance
+ if sampling_dtype is not None:
+ self._met = SamplingDtypeEnabler(self._met, sampling_dtype)
def _checkEquivalence(self, newdom):
newdom = makeDomain(newdom)
@@ -230,7 +313,7 @@ class PoissonianEnergy(EnergyOperator):
res = x.sum() - x.ptw("log").vdot(self._d)
if not x.want_metric:
return res
- return res.add_metric(makeOp(1./x.val))
+ return res.add_metric(SamplingDtypeEnabler(makeOp(1./x.val), np.float64))
class InverseGammaLikelihood(EnergyOperator):
@@ -264,13 +347,20 @@ class InverseGammaLikelihood(EnergyOperator):
elif not isinstance(alpha, Field):
raise TypeError
self._alphap1 = alpha+1
+ if not self._beta.dtype == np.float64:
+ # FIXME Add proper complex support for this energy
+ raise TypeError
+ self._sampling_dtype = _field_to_dtype(self._beta)
def apply(self, x):
self._check_input(x)
res = x.ptw("log").vdot(self._alphap1) + x.ptw("reciprocal").vdot(self._beta)
if not x.want_metric:
return res
- return res.add_metric(makeOp(self._alphap1/(x.val**2)))
+ met = makeOp(self._alphap1/(x.val**2))
+ if self._sampling_dtype is not None:
+ met = SamplingDtypeEnabler(met, self._sampling_dtype)
+ return res.add_metric(met)
class StudentTEnergy(EnergyOperator):
@@ -290,9 +380,12 @@ class StudentTEnergy(EnergyOperator):
Degree of freedom parameter for the student t distribution
"""
- def __init__(self, domain, theta):
+ def __init__(self, domain, theta, sampling_dtype=np.float64):
self._domain = DomainTuple.make(domain)
self._theta = theta
+ self._sampling_dtype = sampling_dtype
+ if sampling_dtype == np.complex128:
+ raise NotImplementedError('Complex data not supported yet')
def apply(self, x):
self._check_input(x)
@@ -300,6 +393,8 @@ class StudentTEnergy(EnergyOperator):
if not x.want_metric:
return res
met = makeOp((self._theta+1) / (self._theta+3), self.domain)
+ if self._sampling_dtype is not None:
+ met = SamplingDtypeEnabler(met, self._sampling_dtype)
return res.add_metric(met)
@@ -323,7 +418,7 @@ class BernoulliEnergy(EnergyOperator):
def __init__(self, d):
if not isinstance(d, Field) or not np.issubdtype(d.dtype, np.integer):
raise TypeError
- if not np.all(np.logical_or(d.val == 0, d.val == 1)):
+ if np.any(np.logical_and(d.val != 0, d.val != 1)):
raise ValueError
self._d = d
self._domain = DomainTuple.make(d.domain)
@@ -333,7 +428,8 @@ class BernoulliEnergy(EnergyOperator):
res = -x.ptw("log").vdot(self._d) + (1.-x).ptw("log").vdot(self._d-1.)
if not x.want_metric:
return res
- return res.add_metric(makeOp(1./(x.val*(1. - x.val))))
+ met = makeOp(1./(x.val*(1. - x.val)))
+ return res.add_metric(SamplingDtypeEnabler(met, np.float64))
class StandardHamiltonian(EnergyOperator):
@@ -362,6 +458,7 @@ class StandardHamiltonian(EnergyOperator):
ic_samp : IterationController
Tells an internal :class:`SamplingEnabler` which convergence criterion
to use to draw Gaussian samples.
+ sampling_dtype : FIXME
See also
--------
@@ -370,9 +467,9 @@ class StandardHamiltonian(EnergyOperator):
`<https://arxiv.org/abs/1812.04403>`_
"""
- def __init__(self, lh, ic_samp=None, _c_inp=None):
+ def __init__(self, lh, ic_samp=None, _c_inp=None, sampling_dtype=np.float64):
self._lh = lh
- self._prior = GaussianEnergy(domain=lh.domain)
+ self._prior = GaussianEnergy(domain=lh.domain, sampling_dtype=sampling_dtype)
if _c_inp is not None:
_, self._prior = self._prior.simplify_for_constant_input(_c_inp)
self._ic_samp = ic_samp
diff --git a/nifty6/operators/inversion_enabler.py b/nifty6/operators/inversion_enabler.py
index a82c3b45a..0a787ce8e 100644
--- a/nifty6/operators/inversion_enabler.py
+++ b/nifty6/operators/inversion_enabler.py
@@ -11,12 +11,10 @@
# 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-2019 Max-Planck-Society
+# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
-import numpy as np
-
from ..logger import logger
from ..minimization.conjugate_gradient import ConjugateGradient
from ..minimization.iteration_controllers import IterationController
@@ -78,5 +76,8 @@ class InversionEnabler(EndomorphicOperator):
logger.warning("Error detected during operator inversion")
return r.position
- def draw_sample(self, dtype, from_inverse=False):
- return self._op.draw_sample(dtype, from_inverse)
+ def draw_sample(self, from_inverse=False):
+ return self._op.draw_sample(from_inverse)
+
+ def draw_sample_with_dtype(self, dtype, from_inverse=False):
+ return self._op.draw_sample_with_dtype(dtype, from_inverse)
diff --git a/nifty6/operators/linear_operator.py b/nifty6/operators/linear_operator.py
index f9f24b504..3e0bec38a 100644
--- a/nifty6/operators/linear_operator.py
+++ b/nifty6/operators/linear_operator.py
@@ -11,12 +11,10 @@
# 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-2019 Max-Planck-Society
+# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
-from ..field import Field
-from ..multi_field import MultiField
from .operator import Operator
diff --git a/nifty6/operators/operator_adapter.py b/nifty6/operators/operator_adapter.py
index 288758814..9cad550aa 100644
--- a/nifty6/operators/operator_adapter.py
+++ b/nifty6/operators/operator_adapter.py
@@ -11,12 +11,10 @@
# 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-2019 Max-Planck-Society
+# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
-import numpy as np
-
from .linear_operator import LinearOperator
@@ -58,10 +56,15 @@ class OperatorAdapter(LinearOperator):
return self._op.apply(x,
self._modeTable[self._trafo][self._ilog[mode]])
- def draw_sample(self, dtype, from_inverse=False):
+ def draw_sample(self, from_inverse=False):
+ if self._trafo & self.INVERSE_BIT:
+ return self._op.draw_sample(not from_inverse)
+ return self._op.draw_sample(from_inverse)
+
+ def draw_sample_with_dtype(self, dtype, from_inverse=False):
if self._trafo & self.INVERSE_BIT:
- return self._op.draw_sample(dtype, not from_inverse)
- return self._op.draw_sample(dtype, from_inverse)
+ return self._op.draw_sample_with_dtype(dtype, not from_inverse)
+ return self._op.draw_sample_with_dtype(dtype, from_inverse)
def __repr__(self):
from ..utilities import indent
diff --git a/nifty6/operators/sampling_enabler.py b/nifty6/operators/sampling_enabler.py
index a3ad95df9..fe9d5f203 100644
--- a/nifty6/operators/sampling_enabler.py
+++ b/nifty6/operators/sampling_enabler.py
@@ -62,19 +62,19 @@ class SamplingEnabler(EndomorphicOperator):
self._domain = self._op.domain
self._capability = self._op.capability
- def draw_sample(self, dtype, from_inverse=False):
+ def draw_sample(self, from_inverse=False):
try:
- return self._op.draw_sample(dtype, from_inverse)
+ return self._op.draw_sample(from_inverse)
except NotImplementedError:
if not from_inverse:
raise ValueError("from_inverse must be True here")
if self._start_from_zero:
- b = self._op.draw_sample(dtype)
+ b = self._op.draw_sample()
energy = QuadraticEnergy(0*b, self._op, b)
else:
- s = self._prior.draw_sample(dtype, from_inverse=True)
+ s = self._prior.draw_sample(from_inverse=True)
sp = self._prior(s)
- nj = self._likelihood.draw_sample(dtype)
+ nj = self._likelihood.draw_sample()
energy = QuadraticEnergy(s, self._op, sp + nj,
_grad=self._likelihood(s) - nj)
inverter = ConjugateGradient(self._ic)
@@ -85,6 +85,9 @@ class SamplingEnabler(EndomorphicOperator):
energy, convergence = inverter(energy)
return energy.position
+ def draw_sample_with_dtype(self, dtype, from_inverse=False):
+ return self._op.draw_sample_with_dtype(dtype, from_inverse)
+
def apply(self, x, mode):
return self._op.apply(x, mode)
diff --git a/nifty6/operators/sandwich_operator.py b/nifty6/operators/sandwich_operator.py
index 5a9ecb204..26e59bee9 100644
--- a/nifty6/operators/sandwich_operator.py
+++ b/nifty6/operators/sandwich_operator.py
@@ -74,12 +74,12 @@ class SandwichOperator(EndomorphicOperator):
def apply(self, x, mode):
return self._op.apply(x, mode)
- def draw_sample(self, dtype, from_inverse=False):
+ def draw_sample(self, from_inverse=False):
# Inverse samples from general sandwiches are not possible
if from_inverse:
if self._bun.capability & self._bun.INVERSE_TIMES:
try:
- s = self._cheese.draw_sample(dtype, from_inverse)
+ s = self._cheese.draw_sample(from_inverse)
return self._bun.inverse_times(s)
except NotImplementedError:
pass
@@ -88,7 +88,7 @@ class SandwichOperator(EndomorphicOperator):
# Samples from general sandwiches
return self._bun.adjoint_times(
- self._cheese.draw_sample(dtype, from_inverse))
+ self._cheese.draw_sample(from_inverse))
def __repr__(self):
from ..utilities import indent
diff --git a/nifty6/operators/scaling_operator.py b/nifty6/operators/scaling_operator.py
index 11f1175af..29573ace6 100644
--- a/nifty6/operators/scaling_operator.py
+++ b/nifty6/operators/scaling_operator.py
@@ -11,7 +11,7 @@
# 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-2019 Max-Planck-Society
+# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
@@ -91,7 +91,7 @@ class ScalingOperator(EndomorphicOperator):
raise ValueError("operator not positive definite")
return 1./np.sqrt(fct) if from_inverse else np.sqrt(fct)
- def draw_sample(self, dtype, from_inverse=False):
+ def draw_sample_with_dtype(self, dtype, from_inverse=False):
from ..sugar import from_random
return from_random(random_type="normal", domain=self._domain,
std=self._get_fct(from_inverse), dtype=dtype)
diff --git a/nifty6/operators/sum_operator.py b/nifty6/operators/sum_operator.py
index 47da40037..e758d57d6 100644
--- a/nifty6/operators/sum_operator.py
+++ b/nifty6/operators/sum_operator.py
@@ -190,7 +190,7 @@ class SumOperator(LinearOperator):
res = res.flexible_addsub(tmp, neg)
return res
- def draw_sample(self, dtype, from_inverse=False):
+ def draw_sample(self, from_inverse=False):
if from_inverse:
raise NotImplementedError(
"cannot draw from inverse of this operator")
@@ -199,7 +199,7 @@ class SumOperator(LinearOperator):
from .simple_linear_operators import NullOperator
if isinstance(op, NullOperator):
continue
- tmp = op.draw_sample(dtype, from_inverse)
+ tmp = op.draw_sample(from_inverse)
res = tmp if res is None else res.unite(tmp)
return res
diff --git a/nifty6/probing.py b/nifty6/probing.py
index cf08e9d13..7aa43402c 100644
--- a/nifty6/probing.py
+++ b/nifty6/probing.py
@@ -100,9 +100,9 @@ def probe_with_posterior_samples(op, post_op, nprobes, dtype):
sc = StatCalculator()
for i in range(nprobes):
if post_op is None:
- sc.add(op.draw_sample(dtype=dtype, from_inverse=True))
+ sc.add(op.draw_sample(from_inverse=True))
else:
- sc.add(post_op(op.draw_sample(dtype=dtype, from_inverse=True)))
+ sc.add(post_op(op.draw_sample(from_inverse=True)))
if nprobes == 1:
return sc.mean, None
diff --git a/test/test_energy_gradients.py b/test/test_energy_gradients.py
index 975c3332b..b9fe53d12 100644
--- a/test/test_energy_gradients.py
+++ b/test/test_energy_gradients.py
@@ -40,7 +40,7 @@ pmp = pytest.mark.parametrize
def field(request):
with ift.random.Context(request.param[0]):
S = ift.ScalingOperator(request.param[1], 1.)
- return S.draw_sample(dtype=np.float64)
+ return S.draw_sample_with_dtype(dtype=np.float64)
def test_gaussian(field):
@@ -50,7 +50,7 @@ def test_gaussian(field):
def test_ScaledEnergy(field):
icov = ift.ScalingOperator(field.domain, 1.2)
- energy = ift.GaussianEnergy(inverse_covariance=icov)
+ energy = ift.GaussianEnergy(inverse_covariance=icov, sampling_dtype=np.float64)
ift.extra.check_jacobian_consistency(energy.scale(0.3), field)
lin = ift.Linearization.make_var(field, want_metric=True)
@@ -59,12 +59,13 @@ def test_ScaledEnergy(field):
res1 = met1(field)
res2 = met2(field)/0.3
ift.extra.assert_allclose(res1, res2, 0, 1e-12)
- met2.draw_sample(dtype=np.float64)
+ met1.draw_sample()
+ met2.draw_sample()
def test_QuadraticFormOperator(field):
op = ift.ScalingOperator(field.domain, 1.2)
- endo = ift.makeOp(op.draw_sample(dtype=np.float64))
+ endo = ift.makeOp(op.draw_sample_with_dtype(dtype=np.float64))
energy = ift.QuadraticFormOperator(endo)
ift.extra.check_jacobian_consistency(energy, field)
@@ -85,8 +86,7 @@ def test_hamiltonian_and_KL(field):
lh = ift.GaussianEnergy(domain=space)
hamiltonian = ift.StandardHamiltonian(lh)
ift.extra.check_jacobian_consistency(hamiltonian, field)
- S = ift.ScalingOperator(space, 1.)
- samps = [S.draw_sample(dtype=np.float64) for i in range(3)]
+ samps = [ift.from_random('normal', space) for i in range(3)]
kl = ift.AveragedEnergy(hamiltonian, samps)
ift.extra.check_jacobian_consistency(kl, field)
@@ -96,9 +96,9 @@ def test_variablecovariancegaussian(field):
return
dc = {'a': field, 'b': field.ptw("exp")}
mf = ift.MultiField.from_dict(dc)
- energy = ift.VariableCovarianceGaussianEnergy(field.domain, 'a', 'b')
+ energy = ift.VariableCovarianceGaussianEnergy(field.domain, 'a', 'b', np.float64)
ift.extra.check_jacobian_consistency(energy, mf, tol=1e-6)
- energy(ift.Linearization.make_var(mf, want_metric=True)).metric.draw_sample(dtype=np.float64)
+ energy(ift.Linearization.make_var(mf, want_metric=True)).metric.draw_sample()
def test_inverse_gamma(field):
diff --git a/test/test_field.py b/test/test_field.py
index 10bbfd70b..5c9c0da77 100644
--- a/test/test_field.py
+++ b/test/test_field.py
@@ -65,7 +65,7 @@ def test_power_synthesize_analyze(space1, space2):
sc1 = ift.StatCalculator()
sc2 = ift.StatCalculator()
for ii in range(samples):
- sk = opfull.draw_sample(dtype=np.float64)
+ sk = opfull.draw_sample_with_dtype(dtype=np.float64)
sp = ift.power_analyze(sk, spaces=(0, 1), keep_phase_information=False)
sc1.add(sp.sum(spaces=1)/fp2.s_sum())
@@ -93,7 +93,7 @@ def test_DiagonalOperator_power_analyze2(space1, space2):
sc2 = ift.StatCalculator()
for ii in range(samples):
- sk = S_full.draw_sample(dtype=np.float64)
+ sk = S_full.draw_sample_with_dtype(dtype=np.float64)
sp = ift.power_analyze(sk, spaces=(0, 1), keep_phase_information=False)
sc1.add(sp.sum(spaces=1)/fp2.s_sum())
sc2.add(sp.sum(spaces=0)/fp1.s_sum())
diff --git a/test/test_gaussian_energy.py b/test/test_gaussian_energy.py
index 94150f24f..69a5e9efe 100644
--- a/test/test_gaussian_energy.py
+++ b/test/test_gaussian_energy.py
@@ -53,7 +53,7 @@ def test_gaussian_energy(space, nonlinearity, noise, seed):
pspec = ift.PS_field(pspace, pspec)
A = Dist(pspec.ptw("sqrt"))
N = ift.ScalingOperator(space, noise)
- n = N.draw_sample(dtype=np.float64)
+ n = N.draw_sample_with_dtype(dtype=np.float64)
R = ift.ScalingOperator(space, 10.)
def d_model():
diff --git a/test/test_kl.py b/test/test_kl.py
index 4cf896621..d0c656a8a 100644
--- a/test/test_kl.py
+++ b/test/test_kl.py
@@ -34,7 +34,8 @@ def test_kl(constants, point_estimates, mirror_samples, mf):
op = ift.HarmonicSmoothingOperator(dom, 3)
if mf:
op = ift.ducktape(dom, None, 'a')*(op.ducktape('b'))
- lh = ift.GaussianEnergy(domain=op.target) @ op
+ import numpy as np
+ lh = ift.GaussianEnergy(domain=op.target, sampling_dtype=np.float64) @ op
ic = ift.GradientNormController(iteration_limit=5)
h = ift.StandardHamiltonian(lh, ic_samp=ic)
mean0 = ift.from_random('normal', h.domain)
diff --git a/test/test_minimizers.py b/test/test_minimizers.py
index 369dffd3e..df537c723 100644
--- a/test/test_minimizers.py
+++ b/test/test_minimizers.py
@@ -86,15 +86,17 @@ def test_WF_curvature(space):
N = ift.DiagonalOperator(n)
all_diag = 1./s + r**2/n
curv = ift.WienerFilterCurvature(R, N, S, iteration_controller=IC,
- iteration_controller_sampling=IC)
+ iteration_controller_sampling=IC,
+ data_sampling_dtype=np.float64,
+ prior_sampling_dtype=np.float64)
m = curv.inverse(required_result)
assert_allclose(
m.val,
1./all_diag.val,
rtol=1e-3,
atol=1e-3)
- curv.draw_sample(dtype=np.float64)
- curv.draw_sample(dtype=np.float64, from_inverse=True)
+ curv.draw_sample()
+ curv.draw_sample(from_inverse=True)
if len(space.shape) == 1:
R = ift.ValueInserter(space, [0])
@@ -103,15 +105,17 @@ def test_WF_curvature(space):
all_diag = 1./s + R(1/n)
curv = ift.WienerFilterCurvature(R.adjoint, N, S,
iteration_controller=IC,
- iteration_controller_sampling=IC)
+ iteration_controller_sampling=IC,
+ data_sampling_dtype=np.float64,
+ prior_sampling_dtype=np.float64)
m = curv.inverse(required_result)
assert_allclose(
m.val,
1./all_diag.val,
rtol=1e-3,
atol=1e-3)
- curv.draw_sample(dtype=np.float64)
- curv.draw_sample(dtype=np.float64, from_inverse=True)
+ curv.draw_sample()
+ curv.draw_sample(from_inverse=True)
@pmp('minimizer', minimizers + newton_minimizers)
diff --git a/test/test_mpi/test_kl.py b/test/test_mpi/test_kl.py
index 7cd35e130..ccd713641 100644
--- a/test/test_mpi/test_kl.py
+++ b/test/test_mpi/test_kl.py
@@ -11,10 +11,11 @@
# 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-2019 Max-Planck-Society
+# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+import numpy as np
import pytest
from mpi4py import MPI
from numpy.testing import assert_, assert_allclose
@@ -46,7 +47,7 @@ def test_kl(constants, point_estimates, mirror_samples, mode, mf):
op = ift.HarmonicSmoothingOperator(dom, 3)
if mf:
op = ift.ducktape(dom, None, 'a')*(op.ducktape('b'))
- lh = ift.GaussianEnergy(domain=op.target) @ op
+ lh = ift.GaussianEnergy(domain=op.target, sampling_dtype=np.float64) @ op
ic = ift.GradientNormController(iteration_limit=5)
h = ift.StandardHamiltonian(lh, ic_samp=ic)
mean0 = ift.from_random('normal', h.domain)
diff --git a/test/test_operators/test_interpolated.py b/test/test_operators/test_interpolated.py
index 08c6abfd8..985f69c8e 100644
--- a/test/test_operators/test_interpolated.py
+++ b/test/test_operators/test_interpolated.py
@@ -33,8 +33,7 @@ seed = list2fixture([4, 78, 23])
def testInterpolationAccuracy(space, seed):
- S = ift.ScalingOperator(space, 1.)
- pos = S.draw_sample(dtype=np.float64)
+ pos = ift.from_random('normal', space)
alpha = 1.5
qs = [0.73, pos.ptw("exp").val]
for q in qs:
diff --git a/test/test_operators/test_jacobian.py b/test/test_operators/test_jacobian.py
index cd9259430..6b0aff862 100644
--- a/test/test_operators/test_jacobian.py
+++ b/test/test_operators/test_jacobian.py
@@ -39,8 +39,7 @@ seed = list2fixture([4, 78, 23])
def testBasics(space, seed):
with ift.random.Context(seed):
- S = ift.ScalingOperator(space, 1.)
- s = S.draw_sample(dtype=np.float64)
+ s = ift.from_random('normal', space)
var = ift.Linearization.make_var(s)
model = ift.ScalingOperator(var.target, 6.)
ift.extra.check_jacobian_consistency(model, var.val)
@@ -91,8 +90,7 @@ def testBinary(type1, type2, space, seed):
def testSpecialDistributionOps(space, seed):
with ift.random.Context(seed):
- S = ift.ScalingOperator(space, 1.)
- pos = S.draw_sample(dtype=np.float64)
+ pos = ift.from_random('normal', space)
alpha = 1.5
q = 0.73
model = ift.InverseGammaOperator(space, alpha, q)
@@ -104,9 +102,8 @@ def testSpecialDistributionOps(space, seed):
@pmp('neg', [True, False])
def testAdder(space, seed, neg):
with ift.random.Context(seed):
- S = ift.ScalingOperator(space, 1.)
- f = S.draw_sample(dtype=np.float64)
- f1 = S.draw_sample(dtype=np.float64)
+ f = ift.from_random('normal', space)
+ f1 = ift.from_random('normal', space)
op = ift.Adder(f1, neg)
ift.extra.check_jacobian_consistency(op, f)
op = ift.Adder(f1.val.ravel()[0], neg=neg, domain=space)
@@ -130,8 +127,7 @@ def testDynamicModel(target, causal, minimum_phase, seed):
'minimum_phase': minimum_phase
}
model, _ = ift.dynamic_operator(**dct)
- S = ift.ScalingOperator(model.domain, 1.)
- pos = S.draw_sample(dtype=np.float64)
+ pos = ift.from_random('normal', model.domain)
# FIXME I dont know why smaller tol fails for 3D example
ift.extra.check_jacobian_consistency(model, pos, tol=1e-5, ntries=20)
if len(target.shape) > 1:
@@ -151,8 +147,7 @@ def testDynamicModel(target, causal, minimum_phase, seed):
dct['sigc'] = 1.
dct['quant'] = 5
model, _ = ift.dynamic_lightcone_operator(**dct)
- S = ift.ScalingOperator(model.domain, 1.)
- pos = S.draw_sample(dtype=np.float64)
+ pos = ift.from_random('normal', model.domain)
# FIXME I dont know why smaller tol fails for 3D example
ift.extra.check_jacobian_consistency(
model, pos, tol=1e-5, ntries=20)
diff --git a/test/test_sugar.py b/test/test_sugar.py
index 90e0f9882..174a3f0d0 100644
--- a/test/test_sugar.py
+++ b/test/test_sugar.py
@@ -42,7 +42,7 @@ def test_exec_time():
dom = ift.RGSpace(12, harmonic=True)
op = ift.HarmonicTransformOperator(dom)
op1 = op.ptw("exp")
- lh = ift.GaussianEnergy(domain=op.target) @ op1
+ lh = ift.GaussianEnergy(domain=op.target, sampling_dtype=np.float64) @ op1
ic = ift.GradientNormController(iteration_limit=2)
ham = ift.StandardHamiltonian(lh, ic_samp=ic)
kl = ift.MetricGaussianKL(ift.full(ham.domain, 0.), ham, 1)
--
GitLab