diff --git a/demos/getting_started_0.ipynb b/demos/getting_started_0.ipynb
index ff2cddb7037914f12f93000bac39acb9999fef02..cce66e1f22b5151cb76e5fdcc6654962b45ea45c 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 22d7d5c80ca073fbedd641b3bd6cca26c6b95719..c00bbe7436f5ad2ddcd14041db142536e5ec4e14 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 7fb3f716c248f189fbba364319980642601dca7f..1d5ec727fce130eaa4e13bfbf5bbfba0926425c1 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 68e9b78c83cac32d5c7987bac57ac0af7a4dbd5b..c66467cbd8bce9eb490fb533ce9226a228bb2106 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 173e1b0ce786a13c3f976b187c577796fea33bbc..3d0e6eabc7d88f88e5b3c68420381a4e56023c9f 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 83156aa2323218b67304b3a06b61594edb2af668..abd63c2bd8e0914af3406600694feafe1c346d39 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 6c30ac0870f47b1cf8a344043b1e47a758f48125..7872b462cbea7969467e481a2454ce2af6a99764 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 341234904debbe83082ee2c2da965a47507e68a5..d6a5aadbfdb8f4f792f326a11ed2a0bb222e2c98 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 4b4e94a1c7c48ae870b56c6e45d20b892b38d39c..d2fb275e1caf273c37d1063163ad544040c727f0 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 69ece97a683b34b2e2e2a52ffb8a4615b10f99c1..21055550d2180e9f6d8b46380068019b6b44b691 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 2c423c6df415771369ec246bf37e7f5769417ad0..0f68e96b1bc5152cf669d64564b1268cbcd66d56 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 52eff6b3d2589c502d27d18a8651fea73b45b179..def9e7d0e822276e4496d9215d5b1366f709d172 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 a82c3b45ae796ff71d7814efba6fa46ec5135b19..0a787ce8eb4c9efcfe54b251bcc880ad6ee91a59 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 f9f24b50487bbfefb88c6883b7077267c60abeee..3e0bec38a9931c13fd4ab349e7506e39fc56aa7d 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 2887588148e0172f02e032e8ffdeabd98e6b524e..9cad550aabe223aa462f43a312d1f15cd3e05cd4 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 a3ad95df95c6759fe99391046fc752ad1261a064..fe9d5f203f0aae0ced3e4f87c129270ded1aafed 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 5a9ecb2042ee379f49960ca126ed47cd5a67f6b2..26e59bee9ee2b5c335452812bf2ab892d5efb9b8 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 11f1175af0cdf23be1dbef60cd7a5b07c2cb9950..29573ace69b435c80a60e825039fcc6519838381 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 47da4003779f5870a7439f9f2c020b452e85a432..e758d57d6b8a9f0be11401e56141fd3de0813af6 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 cf08e9d13691343d128969056304cb58321f7203..7aa43402ca2f6647cc4da5d53c1a3cb3a583c131 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 975c3332be94c0f0fa648e532dc2009f23fade7b..b9fe53d12287815bfae2b3cc51ab4b24b2191db4 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 10bbfd70bd2f228a4dfaaa44ccfac06cf53c5a38..5c9c0da77c5615bf3bf1cb9c413f1e724493af49 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 94150f24f7671f785cede5699b4d73b31799e8af..69a5e9efe1accf4c2174660fa77551378ac3d4fc 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 4cf89662134b1a3741bd6350ab2b1c8ae184e952..d0c656a8aeeacd25690e6056b43f4510d91ffab9 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 369dffd3e8697cf1bc94ba5d950dbedd004c071d..df537c723651427cf4f708d926c37457d84a1d40 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 7cd35e1306a983a06a86ac931230bbd17bf8b8c2..ccd713641424425a06b232169984c42eece89a01 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 08c6abfd81fe9b570890a888268674451daad2ba..985f69c8ef15ea2310fb860a7a827b285dbbcc0f 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 cd9259430877a96b61d754e3d3b37ca7ed120371..6b0aff862ff58760f446a5157ffc3fb75446daff 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 90e0f98824a6ac61f9fbe1d704913724a67a9d44..174a3f0d07f3cf32f1cedcf5fc168ff2c8b99b12 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)