diff --git a/ChangeLog b/ChangeLog
index 2716926b4a3745150aff6cfd36d63281d5f9ee3f..4cc565de5306bb40a3264dd6eb219fa4a42a2f17 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,20 @@
Changes since NIFTy 5:
+MPI parallelisation over samples in MetricGaussianKL
+====================================================
+
+The classes `MetricGaussianKL` and `MetricGaussianKL_MPI` have been unified
+into one `MetricGaussianKL` class which has MPI support built in.
+
+New approach for random number generation
+=========================================
+
+The code now uses `numpy`'s new `SeedSequence` and `Generator` classes for the
+production of random numbers (introduced in numpy 1.17. This greatly simplifies
+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 non-linear Operators
=========================================
@@ -17,7 +32,6 @@ behaviour since both `Operator._check_input()` and
`extra.check_jacobian_consistency()` tests for the new conditions to be
fulfilled.
-
Special functions for complete Field reduction operations
=========================================================
@@ -66,12 +80,3 @@ User-visible changes:
replaced by a single function called `makeField`
- the property `local_shape` has been removed from `Domain` (and subclasses)
and `DomainTuple`.
-
-Transfer of MPI parallelization into operators:
-===============================================
-
-As was already the case with the `MetricGaussianKL_MPI` in NIFTy5, MPI
-parallelization in NIFTy6 is handled by specialized MPI-enabled operators.
-
-They are accessible via the `nifty6.mpi` namespace, from which they can be
-imported directly: `from nifty6.mpi import MPIenabledOperator`.
diff --git a/demos/Wiener_Filter.ipynb b/demos/Wiener_Filter.ipynb
index b0ef1f8b38a3e372d0c120a45bceaabe287e35c0..ba46b6958af7cb8283eb71f5b1d401f166fc9447 100644
--- a/demos/Wiener_Filter.ipynb
+++ b/demos/Wiener_Filter.ipynb
@@ -140,7 +140,6 @@
"outputs": [],
"source": [
"import numpy as np\n",
- "np.random.seed(40)\n",
"import nifty6 as ift\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline"
diff --git a/demos/bench_gridder.py b/demos/bench_gridder.py
index 9bf476737238bc400566ce8e118aee9db3272977..d07b823c4d81938536ba3d4085f113f9c6141de1 100644
--- a/demos/bench_gridder.py
+++ b/demos/bench_gridder.py
@@ -5,8 +5,6 @@ import numpy as np
import nifty6 as ift
-np.random.seed(40)
-
N0s, a0s, b0s, c0s = [], [], [], []
for ii in range(10, 26):
@@ -15,15 +13,16 @@ for ii in range(10, 26):
N = int(2**ii)
print('N = {}'.format(N))
- uv = np.random.rand(N, 2) - 0.5
- vis = np.random.randn(N) + 1j*np.random.randn(N)
+ rng = ift.random.current_rng()
+ uv = rng.uniform(-.5, .5, (N,2))
+ vis = rng.normal(0., 1., N) + 1j*rng.normal(0., 1., N)
uvspace = ift.RGSpace((nu, nv))
visspace = ift.TMP_UnstructuredSpace(N)
- img = np.random.randn(nu*nv)
- img = img.reshape((nu, nv))
+ img = rng.standard_normal((nu, nv))
+ img = ift.makeField(uvspace, img)
img = ift.makeTMP_fld(uvspace, img)
t0 = time()
diff --git a/demos/bernoulli_demo.py b/demos/bernoulli_demo.py
index fc761fb18e3cacd209fc93daf4424bab8318d612..de728a4695096b03cd48951af08150873b7387ff 100644
--- a/demos/bernoulli_demo.py
+++ b/demos/bernoulli_demo.py
@@ -27,8 +27,6 @@ import numpy as np
import nifty6 as ift
if __name__ == '__main__':
- np.random.seed(41)
-
# Set up the position space of the signal
mode = 2
if mode == 0:
@@ -62,7 +60,7 @@ if __name__ == '__main__':
p = R(sky)
mock_position = ift.from_random('normal', harmonic_space)
tmp = p(mock_position).val.astype(np.float64)
- data = np.random.binomial(1, tmp)
+ data = ift.random.current_rng().binomial(1, tmp)
data = ift.TMP_fld.from_raw(R.target, data)
# Compute likelihood and Hamiltonian
diff --git a/demos/getting_started_1.py b/demos/getting_started_1.py
index 6fb71b5cf52cfb1973ac7e73c5da7921bbf88842..7a02bbd4542a17146480a36784e76781471947ae 100644
--- a/demos/getting_started_1.py
+++ b/demos/getting_started_1.py
@@ -46,8 +46,6 @@ def make_random_mask():
if __name__ == '__main__':
- np.random.seed(42)
-
# Choose space on which the signal field is defined
if len(sys.argv) == 2:
mode = int(sys.argv[1])
diff --git a/demos/getting_started_2.py b/demos/getting_started_2.py
index 44d3e507401fe4d0bed394a4eca87b19c0460e0d..e1da73ef28fa55052fe39cf0599eb841f2a1d87e 100644
--- a/demos/getting_started_2.py
+++ b/demos/getting_started_2.py
@@ -44,8 +44,6 @@ def exposure_2d():
if __name__ == '__main__':
- np.random.seed(42)
-
# Choose space on which the signal field is defined
if len(sys.argv) == 2:
mode = int(sys.argv[1])
@@ -94,7 +92,7 @@ if __name__ == '__main__':
lamb = R(sky)
mock_position = ift.from_random('normal', domain)
data = lamb(mock_position)
- data = np.random.poisson(data.val.astype(np.float64))
+ data = ift.random.current_rng().poisson(data.val.astype(np.float64))
data = ift.TMP_fld.from_raw(d_space, data)
likelihood = ift.PoissonianEnergy(data) @ lamb
diff --git a/demos/getting_started_3.py b/demos/getting_started_3.py
index ff2241a2c98127c211d810182e8d2e9885f00ce5..ac252905bea3c798c7a40315f00da4bd3a0484f2 100644
--- a/demos/getting_started_3.py
+++ b/demos/getting_started_3.py
@@ -33,20 +33,18 @@ import nifty6 as ift
def random_los(n_los):
- starts = list(np.random.uniform(0, 1, (n_los, 2)).T)
- ends = list(np.random.uniform(0, 1, (n_los, 2)).T)
+ starts = list(ift.random.current_rng().random((n_los, 2)).T)
+ ends = list(ift.random.current_rng().random((n_los, 2)).T)
return starts, ends
def radial_los(n_los):
- starts = list(np.random.uniform(0, 1, (n_los, 2)).T)
- ends = list(0.5 + 0*np.random.uniform(0, 1, (n_los, 2)).T)
+ starts = list(ift.random.current_rng().random((n_los, 2)).T)
+ ends = list(0.5 + 0*ift.random.current_rng().random((n_los, 2)).T)
return starts, ends
if __name__ == '__main__':
- np.random.seed(420)
-
# Choose between random line-of-sight response (mode=0) and radial lines
# of sight (mode=1)
if len(sys.argv) == 2:
diff --git a/demos/getting_started_mf.py b/demos/getting_started_mf.py
index b7ad94432981ce77205d36cab6123fb46534c5fb..97cc035ed71d258d6155a6097185e6797721967e 100644
--- a/demos/getting_started_mf.py
+++ b/demos/getting_started_mf.py
@@ -44,20 +44,18 @@ class SingleTMP_Space(ift.LinearOperator):
def random_los(n_los):
- starts = list(np.random.uniform(0, 1, (n_los, 2)).T)
- ends = list(np.random.uniform(0, 1, (n_los, 2)).T)
+ starts = list(ift.random.current_rng().random((n_los, 2)).T)
+ ends = list(ift.random.current_rng().random((n_los, 2)).T)
return starts, ends
def radial_los(n_los):
- starts = list(np.random.uniform(0, 1, (n_los, 2)).T)
- ends = list(0.5 + 0*np.random.uniform(0, 1, (n_los, 2)).T)
+ starts = list(ift.random.current_rng().random((n_los, 2)).T)
+ ends = list(0.5 + 0*ift.random.current_rng().random((n_los, 2)).T)
return starts, ends
if __name__ == '__main__':
- np.random.seed(43)
-
# Choose between random line-of-sight response (mode=0) and radial lines
# of sight (mode=1)
if len(sys.argv) == 2:
diff --git a/demos/polynomial_fit.py b/demos/polynomial_fit.py
index e49a3fb22b47fa38b20975281e78595fffafd1e3..ceb5e7d458eabbec5a84ff9cf34ff2607ec9c8fa 100644
--- a/demos/polynomial_fit.py
+++ b/demos/polynomial_fit.py
@@ -19,7 +19,6 @@ import matplotlib.pyplot as plt
import numpy as np
import nifty6 as ift
-np.random.seed(12)
def polynomial(coefficients, sampling_points):
@@ -86,7 +85,7 @@ if __name__ == '__main__':
N_params = 10
N_samples = 100
size = (12,)
- x = np.random.random(size) * 10
+ x = ift.random.current_rng().random(size) * 10
y = np.sin(x**2) * x**3
var = np.full_like(y, y.var() / 10)
var[-2] *= 4
diff --git a/nifty6/__init__.py b/nifty6/__init__.py
index ece070d7cf01385b180adcc0b3762b5fe97b8a24..a68839a60ef8aaacb06cb7bdf37faaaf576071b8 100644
--- a/nifty6/__init__.py
+++ b/nifty6/__init__.py
@@ -1,5 +1,7 @@
from .version import __version__
+from . import random
+
from .domains.domain import TMP_Space
from .domains.structured_domain import TMP_StructuredSpace
from .domains.unstructured_domain import TMP_UnstructuredSpace
diff --git a/nifty6/extra.py b/nifty6/extra.py
index 60d411b79a9686c2e436c6c29a55eb03c173a54c..7523547206fb2ee2d05d5d025a5c6de4d9bef60d 100644
--- a/nifty6/extra.py
+++ b/nifty6/extra.py
@@ -278,17 +278,21 @@ def check_jacobian_consistency(op, loc, tol=1e-8, ntries=100, perf_check=True):
dir = loc2-loc
locnext = loc2
dirnorm = dir.norm()
+ hist = []
for i in range(50):
locmid = loc + 0.5*dir
linmid = op(Linearization.make_var(locmid))
dirder = linmid.jac(dir)
numgrad = (lin2.val-lin.val)
xtol = tol * dirder.norm() / np.sqrt(dirder.size)
+ hist.append((numgrad-dirder).norm())
+# print(len(hist),hist[-1])
if (abs(numgrad-dirder) <= xtol).s_all():
break
dir = dir*0.5
dirnorm *= 0.5
loc2, lin2 = locmid, linmid
else:
+ print(hist)
raise ValueError("gradient and value seem inconsistent")
loc = locnext
diff --git a/nifty6/library/special_distributions.py b/nifty6/library/special_distributions.py
index 4f115c0e70ec0b58070ef244acd6e1f468a5da09..95d9775080d51c8598818038b13583502ccac4a3 100644
--- a/nifty6/library/special_distributions.py
+++ b/nifty6/library/special_distributions.py
@@ -25,6 +25,7 @@ from ..field import TMP_fld
from ..linearization import Linearization
from ..operators.operator import Operator
from ..sugar import makeOp
+from .. import random
def _f_on_np(f, arr):
@@ -67,7 +68,8 @@ class _InterpolationOperator(Operator):
if table_func is not None:
if inv_table_func is None:
raise ValueError
- a = func(np.random.randn(10))
+# MR FIXME: not sure whether we should have this in production code
+ a = func(random.current_rng().random(10))
a1 = _f_on_np(lambda x: inv_table_func(table_func(x)), a)
np.testing.assert_allclose(a, a1)
self._table = _f_on_np(table_func, self._table)
diff --git a/nifty6/minimization/metric_gaussian_kl.py b/nifty6/minimization/metric_gaussian_kl.py
index 85b1696eb78c54357f46eb328bcb3b59ce692222..7ca0c06e00d511c93cd14710d803e6bce2fb038e 100644
--- a/nifty6/minimization/metric_gaussian_kl.py
+++ b/nifty6/minimization/metric_gaussian_kl.py
@@ -11,19 +11,66 @@
# 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 .. import utilities
+from .. import random, utilities
+from ..field import TMP_fld
from ..linearization import Linearization
+from ..multi_field import TMP_Field
+from ..operators.endomorphic_operator import EndomorphicOperator
from ..operators.energy_operators import StandardHamiltonian
from ..probing import approximation2endo
-from ..sugar import makeOp
+from ..sugar import full, makeOp
from .energy import Energy
+def _shareRange(nwork, nshares, myshare):
+ nbase = nwork//nshares
+ additional = nwork % nshares
+ lo = myshare*nbase + min(myshare, additional)
+ hi = lo + nbase + int(myshare < additional)
+ return lo, hi
+
+
+def _np_allreduce_sum(comm, arr):
+ if comm is None:
+ return arr
+ from mpi4py import MPI
+ arr = np.array(arr)
+ res = np.empty_like(arr)
+ comm.Allreduce(arr, res, MPI.SUM)
+ return res
+
+
+def _allreduce_sum_field(comm, fld):
+ if comm is None:
+ return fld
+ if isinstance(fld, TMP_fld):
+ return TMP_fld(fld.domain, _np_allreduce_sum(fld.val))
+ res = tuple(
+ TMP_fld(f.domain, _np_allreduce_sum(comm, f.val))
+ for f in fld.values())
+ return TMP_Field(fld.domain, res)
+
+
+class _KLMetric(EndomorphicOperator):
+ def __init__(self, KL):
+ self._KL = KL
+ self._capability = self.TIMES | self.ADJOINT_TIMES
+ self._domain = KL.position.domain
+
+ def apply(self, x, mode):
+ self._check_input(x, mode)
+ return self._KL.apply_metric(x)
+
+ def draw_sample(self, from_inverse=False, dtype=np.float64):
+ return self._KL._metric_sample(from_inverse, dtype)
+
+
class MetricGaussianKL(Energy):
"""Provides the sampled Kullback-Leibler divergence between a distribution
and a Metric Gaussian.
@@ -38,6 +85,7 @@ class MetricGaussianKL(Energy):
typically nonlinear structure of the true distribution these samples have
to be updated eventually by intantiating `MetricGaussianKL` again. For the
true probability distribution the standard parametrization is assumed.
+ The samples of this class can be distributed among MPI tasks.
Parameters
----------
@@ -62,6 +110,17 @@ class MetricGaussianKL(Energy):
napprox : int
Number of samples for computing preconditioner for sampling. No
preconditioning is done by default.
+ comm : MPI communicator or None
+ 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.
_samples : None
Only a parameter for internal uses. Typically not to be set by users.
@@ -79,7 +138,8 @@ class MetricGaussianKL(Energy):
def __init__(self, mean, hamiltonian, n_samples, constants=[],
point_estimates=[], mirror_samples=False,
- napprox=0, _samples=None, lh_sampling_dtype=np.float64):
+ napprox=0, comm=None, _samples=None,
+ lh_sampling_dtype=np.float64):
super(MetricGaussianKL, self).__init__(mean)
if not isinstance(hamiltonian, StandardHamiltonian):
@@ -88,47 +148,72 @@ class MetricGaussianKL(Energy):
raise ValueError
if not isinstance(n_samples, int):
raise TypeError
- self._constants = list(constants)
- self._point_estimates = list(point_estimates)
+ self._constants = tuple(constants)
+ self._point_estimates = tuple(point_estimates)
if not isinstance(mirror_samples, bool):
raise TypeError
self._hamiltonian = hamiltonian
+ self._n_samples = int(n_samples)
+ if comm is not None:
+ self._comm = comm
+ ntask = self._comm.Get_size()
+ rank = self._comm.Get_rank()
+ self._lo, self._hi = _shareRange(self._n_samples, ntask, rank)
+ else:
+ self._comm = None
+ self._lo, self._hi = 0, self._n_samples
+
+ self._mirror_samples = bool(mirror_samples)
+ self._n_eff_samples = self._n_samples
+ if self._mirror_samples:
+ self._n_eff_samples *= 2
+
if _samples is None:
met = hamiltonian(Linearization.make_partial_var(
- mean, point_estimates, True)).metric
- if napprox > 1:
+ mean, self._point_estimates, True)).metric
+ if napprox >= 1:
met._approximation = makeOp(approximation2endo(met, napprox))
- _samples = tuple(met.draw_sample(from_inverse=True,
- dtype=lh_sampling_dtype)
- for _ in range(n_samples))
- if mirror_samples:
- _samples += tuple(-s for s in _samples)
+ _samples = []
+ sseq = random.spawn_sseq(self._n_samples)
+ for i in range(self._lo, self._hi):
+ random.push_sseq(sseq[i])
+ _samples.append(met.draw_sample(from_inverse=True,
+ dtype=lh_sampling_dtype))
+ random.pop_sseq()
+ _samples = tuple(_samples)
+ else:
+ if len(_samples) != self._hi-self._lo:
+ raise ValueError("# of samples mismatch")
self._samples = _samples
-
- # FIXME Use simplify for constant input instead
- self._lin = Linearization.make_partial_var(mean, constants)
+ self._lin = Linearization.make_partial_var(mean, self._constants)
v, g = None, None
- for s in self._samples:
- tmp = self._hamiltonian(self._lin+s)
- if v is None:
- v = tmp.val.val[()]
- g = tmp.gradient
- else:
- v += tmp.val.val[()]
- g = g + tmp.gradient
- self._val = v / len(self._samples)
- self._grad = g * (1./len(self._samples))
+ if len(self._samples) == 0: # hack if there are too many MPI tasks
+ tmp = self._hamiltonian(self._lin)
+ v = 0. * tmp.val.val
+ g = 0. * tmp.gradient
+ else:
+ for s in self._samples:
+ tmp = self._hamiltonian(self._lin+s)
+ if self._mirror_samples:
+ tmp = tmp + self._hamiltonian(self._lin-s)
+ if v is None:
+ v = tmp.val.val_rw()
+ g = tmp.gradient
+ else:
+ v += tmp.val.val
+ g = g + tmp.gradient
+ 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._napprox = napprox
self._sampdt = lh_sampling_dtype
def at(self, position):
- return MetricGaussianKL(position, self._hamiltonian, 0,
- self._constants, self._point_estimates,
- napprox=self._napprox, _samples=self._samples,
- lh_sampling_dtype=self._sampdt)
+ return MetricGaussianKL(
+ position, self._hamiltonian, self._n_samples, self._constants,
+ self._point_estimates, self._mirror_samples, comm=self._comm,
+ _samples=self._samples, lh_sampling_dtype=self._sampdt)
@property
def value(self):
@@ -139,30 +224,48 @@ class MetricGaussianKL(Energy):
return self._grad
def _get_metric(self):
+ lin = self._lin.with_want_metric()
if self._metric is None:
- lin = self._lin.with_want_metric()
- mymap = map(lambda v: self._hamiltonian(lin+v).metric,
- self._samples)
- self._unscaled_metric = utilities.my_sum(mymap)
- self._metric = self._unscaled_metric.scale(1./len(self._samples))
-
- def unscaled_metric(self):
- self._get_metric()
- return self._unscaled_metric, 1/len(self._samples)
+ if len(self._samples) == 0: # hack if there are too many MPI tasks
+ self._metric = self._hamiltonian(lin).metric.scale(0.)
+ else:
+ mymap = map(lambda v: self._hamiltonian(lin+v).metric,
+ self._samples)
+ self.unscaled_metric = utilities.my_sum(mymap)
+ self._metric = self.unscaled_metric.scale(1./self._n_eff_samples)
def apply_metric(self, x):
self._get_metric()
- return self._metric(x)
+ return _allreduce_sum_field(self._comm, self._metric(x))
@property
def metric(self):
- self._get_metric()
- return self._metric
+ return _KLMetric(self)
@property
def samples(self):
- return self._samples
+ if self._comm is not None:
+ res = self._comm.allgather(self._samples)
+ res = [item for sublist in res for item in sublist]
+ else:
+ res = self._samples
+ if self._mirror_samples:
+ res = res + tuple(-item for item in res)
+ return res
+
+ def _unscaled_metric_sample(self, from_inverse=False, dtype=np.float64):
+ if from_inverse:
+ raise NotImplementedError()
+ lin = self._lin.with_want_metric()
+ samp = full(self._hamiltonian.domain, 0.)
+ sseq = random.spawn_sseq(self._n_samples)
+ for i, v in enumerate(self._samples):
+ random.push_sseq(sseq[self._lo+i])
+ samp = samp + self._hamiltonian(lin+v).metric.draw_sample(from_inverse=False, dtype=dtype)
+ if self._mirror_samples:
+ samp = samp + self._hamiltonian(lin-v).metric.draw_sample(from_inverse=False, dtype=dtype)
+ random.pop_sseq()
+ return _allreduce_sum_field(self._comm, samp)
- def __repr__(self):
- return 'KL ({} samples):\n'.format(len(
- self._samples)) + utilities.indent(self._hamiltonian.__repr__())
+ def _metric_sample(self, from_inverse=False, dtype=np.float64):
+ return self._unscaled_metric_sample(from_inverse, dtype)/self._n_eff_samples
diff --git a/nifty6/minimization/metric_gaussian_kl_mpi.py b/nifty6/minimization/metric_gaussian_kl_mpi.py
deleted file mode 100644
index 615a5cbc7bac0e2254692ec364de772ab3a4f04a..0000000000000000000000000000000000000000
--- a/nifty6/minimization/metric_gaussian_kl_mpi.py
+++ /dev/null
@@ -1,256 +0,0 @@
-# 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-2019 Max-Planck-Society
-#
-# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
-
-from .. import utilities
-from ..linearization import Linearization
-from ..operators.energy_operators import StandardHamiltonian
-from ..operators.endomorphic_operator import EndomorphicOperator
-from .energy import Energy
-from mpi4py import MPI
-import numpy as np
-from ..probing import approximation2endo
-from ..sugar import makeOp, full
-from ..field import TMP_fld
-from ..multi_field import TMP_Field
-
-
-_comm = MPI.COMM_WORLD
-ntask = _comm.Get_size()
-rank = _comm.Get_rank()
-master = (rank == 0)
-
-
-def _shareRange(nwork, nshares, myshare):
- nbase = nwork//nshares
- additional = nwork % nshares
- lo = myshare*nbase + min(myshare, additional)
- hi = lo + nbase + int(myshare < additional)
- return lo, hi
-
-
-def np_allreduce_sum(arr):
- arr = np.array(arr)
- res = np.empty_like(arr)
- _comm.Allreduce(arr, res, MPI.SUM)
- return res
-
-
-def allreduce_sum_field(fld):
- if isinstance(fld, TMP_fld):
- return TMP_fld(fld.domain, np_allreduce_sum(fld.val))
- res = tuple(
- TMP_fld(f.domain, np_allreduce_sum(f.val))
- for f in fld.values())
- return TMP_Field(fld.domain, res)
-
-
-class KLMetric(EndomorphicOperator):
- def __init__(self, KL):
- self._KL = KL
- self._capability = self.TIMES | self.ADJOINT_TIMES
- self._domain = KL.position.domain
-
- def apply(self, x, mode):
- self._check_input(x, mode)
- return self._KL.apply_metric(x)
-
- def draw_sample(self, from_inverse=False, dtype=np.float64):
- return self._KL.metric_sample(from_inverse, dtype)
-
-
-class MetricGaussianKL_MPI(Energy):
- """Provides the sampled Kullback-Leibler divergence between a distribution
- and a Metric Gaussian.
-
- A Metric Gaussian is used to approximate another probability distribution.
- It is a Gaussian distribution that uses the Fisher information metric of
- the other distribution at the location of its mean to approximate the
- variance. In order to infer the mean, a stochastic estimate of the
- Kullback-Leibler divergence is minimized. This estimate is obtained by
- sampling the Metric Gaussian at the current mean. During minimization
- these samples are kept constant; only the mean is updated. Due to the
- typically nonlinear structure of the true distribution these samples have
- to be updated eventually by intantiating `MetricGaussianKL` again. For the
- true probability distribution the standard parametrization is assumed.
- The samples of this class are distributed among MPI tasks.
-
- Parameters
- ----------
- mean : TMP_fld
- Mean of the Gaussian probability distribution.
- hamiltonian : StandardHamiltonian
- Hamiltonian of the approximated probability distribution.
- n_samples : integer
- Number of samples used to stochastically estimate the KL.
- constants : list
- List of parameter keys that are kept constant during optimization.
- Default is no constants.
- point_estimates : list
- List of parameter keys for which no samples are drawn, but that are
- (possibly) optimized for, corresponding to point estimates of these.
- Default is to draw samples for the complete domain.
- mirror_samples : boolean
- Whether the negative of the drawn samples are also used,
- as they are equally legitimate samples. If true, the number of used
- samples doubles. Mirroring samples stabilizes the KL estimate as
- extreme sample variation is counterbalanced. Default is False.
- napprox : int
- Number of samples for computing preconditioner for sampling. No
- preconditioning is done by default.
- _samples : None
- Only a parameter for internal uses. Typically not to be set by users.
- seed_offset : int
- A parameter with which one can controll from which seed the samples
- are drawn. Per default, the seed is different for MPI tasks, but the
- same every time this class is initialized.
-
- Note
- ----
- The two lists `constants` and `point_estimates` are independent from each
- other. It is possible to sample along domains which are kept constant
- during minimization and vice versa.
-
- See also
- --------
- `Metric Gaussian Variational Inference`, Jakob Knollmüller,
- Torsten A. Enßlin, `<https://arxiv.org/abs/1901.11033>`_
- """
-
- def __init__(self, mean, hamiltonian, n_samples, constants=[],
- point_estimates=[], mirror_samples=False,
- napprox=0, _samples=None, seed_offset=0,
- lh_sampling_dtype=np.float64):
- super(MetricGaussianKL_MPI, self).__init__(mean)
-
- if not isinstance(hamiltonian, StandardHamiltonian):
- raise TypeError
- if hamiltonian.domain is not mean.domain:
- raise ValueError
- if not isinstance(n_samples, int):
- raise TypeError
- self._constants = list(constants)
- self._point_estimates = list(point_estimates)
- if not isinstance(mirror_samples, bool):
- raise TypeError
-
- self._hamiltonian = hamiltonian
-
- if _samples is None:
- if mirror_samples:
- lo, hi = _shareRange(n_samples*2, ntask, rank)
- else:
- lo, hi = _shareRange(n_samples, ntask, rank)
- met = hamiltonian(Linearization.make_partial_var(
- mean, point_estimates, True)).metric
- if napprox > 1:
- met._approximation = makeOp(approximation2endo(met, napprox))
- _samples = []
- rand_state = np.random.get_state()
- for i in range(lo, hi):
- if mirror_samples:
- np.random.seed(i//2+seed_offset)
- if (i % 2) and (i-1 >= lo):
- _samples.append(-_samples[-1])
-
- else:
- _samples.append(((i % 2)*2-1) *
- met.draw_sample(from_inverse=True,
- dtype=lh_sampling_dtype))
- else:
- np.random.seed(i+seed_offset)
- _samples.append(met.draw_sample(from_inverse=True,
- dtype=lh_sampling_dtype))
- np.random.set_state(rand_state)
- _samples = tuple(_samples)
- if mirror_samples:
- n_samples *= 2
- self._samples = _samples
- self._seed_offset = seed_offset
- self._n_samples = n_samples
- self._lin = Linearization.make_partial_var(mean, constants)
- v, g = None, None
- if len(self._samples) == 0: # hack if there are too many MPI tasks
- tmp = self._hamiltonian(self._lin)
- v = 0. * tmp.val.val
- g = 0. * tmp.gradient
- else:
- for s in self._samples:
- tmp = self._hamiltonian(self._lin+s)
- if v is None:
- v = tmp.val.val_rw()
- g = tmp.gradient
- else:
- v += tmp.val.val
- g = g + tmp.gradient
- self._val = np_allreduce_sum(v)[()] / self._n_samples
- self._grad = allreduce_sum_field(g) / self._n_samples
- self._metric = None
- self._sampdt = lh_sampling_dtype
-
- def at(self, position):
- return MetricGaussianKL_MPI(
- position, self._hamiltonian, self._n_samples, self._constants,
- self._point_estimates, _samples=self._samples,
- seed_offset=self._seed_offset, lh_sampling_dtype=self._sampdt)
-
- @property
- def value(self):
- return self._val
-
- @property
- def gradient(self):
- return self._grad
-
- def _get_metric(self):
- lin = self._lin.with_want_metric()
- if self._metric is None:
- if len(self._samples) == 0: # hack if there are too many MPI tasks
- self._metric = self._hamiltonian(lin).metric.scale(0.)
- else:
- mymap = map(lambda v: self._hamiltonian(lin+v).metric,
- self._samples)
- self.unscaled_metric = utilities.my_sum(mymap)
- self._metric = self.unscaled_metric.scale(1./self._n_samples)
-
- def apply_metric(self, x):
- self._get_metric()
- return allreduce_sum_field(self._metric(x))
-
- @property
- def metric(self):
- return KLMetric(self)
-
- @property
- def samples(self):
- res = _comm.allgather(self._samples)
- res = [item for sublist in res for item in sublist]
- return res
-
- def unscaled_metric_sample(self, from_inverse=False, dtype=np.float64):
- if from_inverse:
- raise NotImplementedError()
- lin = self._lin.with_want_metric()
- samp = full(self._hamiltonian.domain, 0.)
- rand_state = np.random.get_state()
- np.random.seed(rank+np.random.randint(99999))
- for v in self._samples:
- samp = samp + self._hamiltonian(lin+v).metric.draw_sample(from_inverse=False, dtype=dtype)
- np.random.set_state(rand_state)
- return allreduce_sum_field(samp)
-
- def metric_sample(self, from_inverse=False, dtype=np.float64):
- return self.unscaled_metric_sample(from_inverse, dtype)/self._n_samples
diff --git a/nifty6/mpi.py b/nifty6/mpi.py
deleted file mode 100644
index b5a3d58ff23faaabf691f31d6ae3ff2164f033c7..0000000000000000000000000000000000000000
--- a/nifty6/mpi.py
+++ /dev/null
@@ -1 +0,0 @@
-from .minimization.metric_gaussian_kl_mpi import MetricGaussianKL_MPI
diff --git a/nifty6/operators/operator.py b/nifty6/operators/operator.py
index ccc6433f4bea4496e5ad107fc117ce26a468576b..b92f02b589380a1e600ba5397d0088433c7f3075 100644
--- a/nifty6/operators/operator.py
+++ b/nifty6/operators/operator.py
@@ -27,10 +27,6 @@ class Operator(metaclass=NiftyMeta):
domain, and can also provide the Jacobian.
"""
- VALUE_ONLY = 0
- WITH_JAC = 1
- WITH_METRIC = 2
-
@property
def domain(self):
"""The domain on which the Operator's input TMP_fld is defined.
@@ -202,7 +198,7 @@ class Operator(metaclass=NiftyMeta):
return self.apply(x.trivial_jac()).prepend_jac(x.jac)
elif isinstance(x, (TMP_fld, TMP_Field)):
return self.apply(x)
- raise TypeError('Operator can only consume TMP_fld, TMP_Fields and Linearizations')
+ return self @ x
def ducktape(self, name):
from .simple_linear_operators import ducktape
diff --git a/nifty6/random.py b/nifty6/random.py
index 673bf6596444d826f149520a26759f8dd423a603..37b0d414343c8e46ac78c9de81fa12e76c846c24 100644
--- a/nifty6/random.py
+++ b/nifty6/random.py
@@ -11,21 +11,150 @@
# 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.
+"""
+Some remarks on NIFTy's treatment of random numbers
+
+NIFTy makes use of the `Generator` and `SeedSequence` classes introduced to
+`numpy.random` in numpy 1.17.
+
+On first load of the `nifty6.random` module, it creates a stack of
+`SeedSequence` objects which contains a single `SeedSequence` with a fixed seed,
+and also a stack of `Generator` objects, which contains a single generator
+derived from the above seed sequence. Without user intervention, this generator
+will be used for all random number generation tasks within NIFTy. This means
+
+- that random numbers drawn by NIFTy will be reproducible across multiple runs
+ (assuming there are no complications like MPI-enabled runs with a varying
+ number of tasks), and
+- that trying to change random seeds via `numpy.random.seed` will have no
+ effect on the random numbers drawn by NIFTy.
+
+Users who want to change the random seed for a given run can achieve this
+by calling `push_sseq_from_seed()` with a seed of their choice. This will push
+a new seed sequence generated from that seed onto the seed sequence stack, and a
+generator derived from this seed sequence onto the generator stack. Since all
+NIFTy RNG-related calls will use the generator on the top of the stack, all
+calls from this point on will use the new generator.
+If the user already has a `SeedSequence` object at hand, they can pass this to
+NIFTy via `push_sseq`. A new generator derived from this sequence will then also
+be pushed onto the generator stack.
+These operations can be reverted (and should be, as soon as the new generator is
+no longer needed) by a call to `pop_sseq()`.
+When users need direct access to the RNG currently in use, they can access it
+via the `current_rng` function.
+
+
+Example for using multiple seed sequences:
+
+Assume that N samples are needed to compute a KL, which are distributed over
+a variable number of MPI tasks. In this situation, whenever random numbers
+need to be drawn for these samples:
+- each MPI task should spawn as many seed sequences as there are samples
+ _in total_, using `sseq = spawn_sseq(N)`
+- each task loops over the local samples
+ - first pushing the seed sequence for the global(!) index of the
+ current sample via `push_sseq(sseq[iglob])`
+ - drawing the required random numbers
+ - then popping the seed sequence again via `pop_sseq()`
+
+That way, random numbers should be reproducible and independent of the number
+of MPI tasks.
+
+WARNING: do not push/pop the same `SeedSequence` object more than once - this
+will lead to repeated random sequences! Whenever you have to push `SeedSequence`
+objects, generate new ones via `spawn_sseq()`.
+"""
+
import numpy as np
+# Stack of SeedSequence objects. Will always start out with a well-defined
+# default. Users can change the "random seed" used by a calculation by pushing
+# a different SeedSequence before invoking any other nifty6.random calls
+_sseq = [np.random.SeedSequence(42)]
+# Stack of random number generators associated with _sseq.
+_rng = [np.random.default_rng(_sseq[-1])]
+
+
+def spawn_sseq(n, parent=None):
+ """Returns a list of `n` SeedSequence objects which are children of `parent`
+
+ Parameters
+ ----------
+ n : int
+ number of requested SeedSequence objects
+ parent : SeedSequence
+ the object from which the returned objects will be derived
+ If `None`, the top of the current SeedSequence stack will be used
+
+ Returns
+ -------
+ list(SeedSequence)
+ the requested SeedSequence objects
+ """
+ if parent is None:
+ global _sseq
+ parent = _sseq[-1]
+ return parent.spawn(n)
+
+
+def current_rng():
+ """Returns the RNG object currently in use by NIFTy
+
+ Returns
+ -------
+ Generator
+ the current Generator object (top of the generatir stack)
+ """
+ return _rng[-1]
+
+
+def push_sseq(sseq):
+ """Pushes a new SeedSequence object onto the SeedSequence stack.
+ This also pushes a new Generator object built from the new SeedSequence
+ to the generator stack.
+
+ Parameters
+ ----------
+ sseq: SeedSequence
+ the SeedSequence object to be used from this point
+ """
+ _sseq.append(sseq)
+ _rng.append(np.random.default_rng(_sseq[-1]))
+
+
+def push_sseq_from_seed(seed):
+ """Pushes a new SeedSequence object derived from an integer seed onto the
+ SeedSequence stack.
+ This also pushes a new Generator object built from the new SeedSequence
+ to the generator stack.
+
+ Parameters
+ ----------
+ seed: int
+ the seed from which the new SeedSequence will be built
+ """
+ _sseq.append(np.random.SeedSequence(seed))
+ _rng.append(np.random.default_rng(_sseq[-1]))
+
+
+def pop_sseq():
+ """Pops the top of the SeedSequence and generator stacks."""
+ _sseq.pop()
+ _rng.pop()
+
class Random(object):
@staticmethod
def pm1(dtype, shape):
if np.issubdtype(dtype, np.complexfloating):
x = np.array([1+0j, 0+1j, -1+0j, 0-1j], dtype=dtype)
- x = x[np.random.randint(4, size=shape)]
+ x = x[_rng[-1].integers(0, 4, size=shape)]
else:
- x = 2*np.random.randint(2, size=shape) - 1
+ x = 2*_rng[-1].integers(0, 2, size=shape)-1
return x.astype(dtype, copy=False)
@staticmethod
@@ -42,10 +171,10 @@ 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 = np.random.normal(mean.real, std*np.sqrt(0.5), shape)
- x.imag = np.random.normal(mean.imag, std*np.sqrt(0.5), shape)
+ 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)
else:
- x = np.random.normal(mean, std, shape).astype(dtype, copy=False)
+ x = _rng[-1].normal(mean, std, shape).astype(dtype, copy=False)
return x
@staticmethod
@@ -57,13 +186,13 @@ class Random(object):
raise TypeError("low and high must not be complex")
if np.issubdtype(dtype, np.complexfloating):
x = np.empty(shape, dtype=dtype)
- x.real = np.random.uniform(low, high, shape)
- x.imag = np.random.uniform(low, high, shape)
+ x.real = _rng[-1].uniform(low, high, shape)
+ x.imag = _rng[-1].uniform(low, high, shape)
elif np.issubdtype(dtype, np.integer):
if not (np.issubdtype(type(low), np.integer) and
np.issubdtype(type(high), np.integer)):
raise TypeError("low and high must be integer")
- x = np.random.randint(low, high+1, shape)
+ x = _rng[-1].integers(low, high+1, shape)
else:
- x = np.random.uniform(low, high, shape)
+ x = _rng[-1].uniform(low, high, shape)
return x.astype(dtype, copy=False)
diff --git a/setup.py b/setup.py
index f909590dcc64d07e03d03ced3a89608e3af76b48..f90a59031a1039c69408275c7f707a95f62eeb32 100644
--- a/setup.py
+++ b/setup.py
@@ -39,8 +39,8 @@ setup(name="nifty6",
packages=find_packages(include=["nifty6", "nifty6.*"]),
zip_safe=True,
license="GPLv3",
- setup_requires=['scipy>=1.4.1'],
- install_requires=['scipy>=1.4.1'],
+ setup_requires=['scipy>=1.4.1', 'numpy>=1.17'],
+ install_requires=['scipy>=1.4.1', 'numpy>=1.17'],
python_requires='>=3.5',
classifiers=[
"Development Status :: 4 - Beta",
diff --git a/test/common.py b/test/common.py
index 9b95150d7ad0ebf1ea0922348fae3fba3a791fd1..ba5224fe61a222cd9ba87f9ac5d6c6f46d78116d 100644
--- a/test/common.py
+++ b/test/common.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.
@@ -24,3 +24,13 @@ def list2fixture(lst):
return request.param
return myfixture
+
+
+def setup_function():
+ import nifty6 as ift
+ ift.random.push_sseq_from_seed(42)
+
+
+def teardown_function():
+ import nifty6 as ift
+ ift.random.pop_sseq()
diff --git a/test/test_energy_gradients.py b/test/test_energy_gradients.py
index 0cd8490fc9225794fea24ec9f4d7b69a74af3e7a..7e78bc1b5452822e1929509ecd599fa325f3b3c6 100644
--- a/test/test_energy_gradients.py
+++ b/test/test_energy_gradients.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.
@@ -20,6 +20,7 @@ import pytest
import nifty6 as ift
from itertools import product
+from .common import setup_function, teardown_function
# Currently it is not possible to parametrize fixtures. But this will
# hopefully be fixed in the future.
@@ -37,9 +38,11 @@ pmp = pytest.mark.parametrize
@pytest.fixture(params=PARAMS)
def field(request):
- np.random.seed(request.param[0])
+ ift.random.push_sseq_from_seed(request.param[0])
S = ift.ScalingOperator(request.param[1], 1.)
- return S.draw_sample()
+ res = S.draw_sample()
+ ift.random.pop_sseq()
+ return res
def test_gaussian(field):
@@ -102,7 +105,7 @@ def test_inverse_gamma(field):
return
field = field.exp()
space = field.domain
- d = np.random.normal(10, size=space.shape)**2
+ d = ift.random.current_rng().normal(10, size=space.shape)**2
d = ift.TMP_fld(space, d)
energy = ift.InverseGammaLikelihood(d)
ift.extra.check_jacobian_consistency(energy, field, tol=1e-5)
@@ -113,10 +116,10 @@ def testPoissonian(field):
return
field = field.exp()
space = field.domain
- d = np.random.poisson(120, size=space.shape)
+ d = ift.random.current_rng().poisson(120, size=space.shape)
d = ift.TMP_fld(space, d)
energy = ift.PoissonianEnergy(d)
- ift.extra.check_jacobian_consistency(energy, field, tol=1e-7)
+ ift.extra.check_jacobian_consistency(energy, field, tol=1e-6)
def test_bernoulli(field):
@@ -124,7 +127,7 @@ def test_bernoulli(field):
return
field = field.sigmoid()
space = field.domain
- d = np.random.binomial(1, 0.1, size=space.shape)
+ d = ift.random.current_rng().binomial(1, 0.1, size=space.shape)
d = ift.TMP_fld(space, d)
energy = ift.BernoulliEnergy(d)
ift.extra.check_jacobian_consistency(energy, field, tol=1e-5)
diff --git a/test/test_field.py b/test/test_field.py
index 5a33dd0dd87efcf5e7c02c2ce36062ead045d907..c8ac577ba9ce015c251944161a848f8f4d5b88ba 100644
--- a/test/test_field.py
+++ b/test/test_field.py
@@ -20,6 +20,7 @@ import pytest
from numpy.testing import assert_allclose, assert_equal, assert_raises
import nifty6 as ift
+from .common import setup_function, teardown_function
pmp = pytest.mark.parametrize
SPACES = [ift.RGSpace((4,)), ift.RGSpace((5))]
@@ -51,8 +52,6 @@ def _spec2(k):
])
@pmp('space2', [ift.RGSpace((8,), harmonic=True), ift.LMSpace(12)])
def test_power_synthesize_analyze(space1, space2):
- np.random.seed(11)
-
p1 = ift.PowerSpace(space1)
fp1 = ift.PS_field(p1, _spec1)
p2 = ift.PowerSpace(space2)
@@ -82,8 +81,6 @@ def test_power_synthesize_analyze(space1, space2):
])
@pmp('space2', [ift.RGSpace((8,), harmonic=True), ift.LMSpace(12)])
def test_DiagonalOperator_power_analyze2(space1, space2):
- np.random.seed(11)
-
fp1 = ift.PS_field(ift.PowerSpace(space1), _spec1)
fp2 = ift.PS_field(ift.PowerSpace(space2), _spec2)
diff --git a/test/test_gaussian_energy.py b/test/test_gaussian_energy.py
index d216f7a4e78b469d1429a74dd15f2258f57901ff..e2cd2535497613508d3d1c89c79e3104ef8a9a78 100644
--- a/test/test_gaussian_energy.py
+++ b/test/test_gaussian_energy.py
@@ -19,6 +19,7 @@ import numpy as np
import pytest
import nifty6 as ift
+from .common import setup_function, teardown_function
def _flat_PS(k):
@@ -37,7 +38,7 @@ pmp = pytest.mark.parametrize
@pmp('noise', [1, 1e-2, 1e2])
@pmp('seed', [4, 78, 23])
def test_gaussian_energy(space, nonlinearity, noise, seed):
- np.random.seed(seed)
+ ift.random.push_sseq_from_seed(seed)
dim = len(space.shape)
hspace = space.get_default_codomain()
ht = ift.HarmonicTransformOperator(hspace, target=space)
@@ -70,4 +71,5 @@ def test_gaussian_energy(space, nonlinearity, noise, seed):
energy = ift.GaussianEnergy(d, N) @ d_model()
ift.extra.check_jacobian_consistency(
- energy, xi0, ntries=10, tol=5e-8)
+ energy, xi0, ntries=10, tol=1e-6)
+ ift.random.pop_sseq()
diff --git a/test/test_kl.py b/test/test_kl.py
index 31ff187558ca7562aeb119626849f2c310a6b257..425da1ca61efda1179cc8ae86002312940dfbd2d 100644
--- a/test/test_kl.py
+++ b/test/test_kl.py
@@ -20,6 +20,7 @@ import numpy as np
import nifty6 as ift
from numpy.testing import assert_, assert_allclose
import pytest
+from .common import setup_function, teardown_function
pmp = pytest.mark.parametrize
@@ -28,7 +29,6 @@ pmp = pytest.mark.parametrize
@pmp('point_estimates', ([], ['a'], ['b'], ['a', 'b']))
@pmp('mirror_samples', (True, False))
def test_kl(constants, point_estimates, mirror_samples):
- np.random.seed(42)
dom = ift.RGSpace((12,), (2.12))
op0 = ift.HarmonicSmoothingOperator(dom, 3)
op = ift.ducktape(dom, None, 'a')*(op0.ducktape('b'))
@@ -45,12 +45,13 @@ def test_kl(constants, point_estimates, mirror_samples):
point_estimates=point_estimates,
mirror_samples=mirror_samples,
napprox=0)
+ samp_full = kl.samples
klpure = ift.MetricGaussianKL(mean0,
h,
- nsamps,
- mirror_samples=mirror_samples,
+ len(samp_full),
+ mirror_samples=False,
napprox=0,
- _samples=kl.samples)
+ _samples=samp_full)
# Test value
assert_allclose(kl.value, klpure.value)
diff --git a/test/test_linearization.py b/test/test_linearization.py
index 3e2605e765533a443c17be52abf87683c4f45347..b44e5b09765db9e01a9d68eb8865307c21f555c2 100644
--- a/test/test_linearization.py
+++ b/test/test_linearization.py
@@ -20,6 +20,7 @@ import pytest
from numpy.testing import assert_, assert_allclose
import nifty6 as ift
+from .common import setup_function, teardown_function
pmp = pytest.mark.parametrize
diff --git a/test/test_minimizers.py b/test/test_minimizers.py
index 9ae0491f53b2013da95f46ec814ce104ba99916b..01bfcbefe3c64002dba1c813b7c596d6acc667c8 100644
--- a/test/test_minimizers.py
+++ b/test/test_minimizers.py
@@ -22,6 +22,7 @@ import pytest
from numpy.testing import assert_allclose, assert_equal
import nifty6 as ift
+from .common import setup_function, teardown_function
pmp = pytest.mark.parametrize
IC = ift.GradientNormController(tol_abs_gradnorm=1e-5, iteration_limit=1000)
@@ -51,9 +52,14 @@ slow_minimizers = ['ift.SteepestDescent(IC)']
slow_minimizers)
@pmp('space', spaces)
def test_quadratic_minimization(minimizer, space):
+<<<<<<< HEAD
np.random.seed(42)
starting_point = ift.TMP_fld.from_random('normal', domain=space)*10
covariance_diagonal = ift.TMP_fld.from_random('uniform', domain=space) + 0.5
+=======
+ starting_point = ift.Field.from_random('normal', domain=space)*10
+ covariance_diagonal = ift.Field.from_random('uniform', domain=space) + 0.5
+>>>>>>> NIFTy_6
covariance = ift.DiagonalOperator(covariance_diagonal)
required_result = ift.full(space, 1.)
@@ -76,7 +82,6 @@ def test_quadratic_minimization(minimizer, space):
@pmp('space', spaces)
def test_WF_curvature(space):
- np.random.seed(42)
required_result = ift.full(space, 1.)
s = ift.TMP_fld.from_random('uniform', domain=space) + 0.5
@@ -121,10 +126,8 @@ def test_rosenbrock(minimizer):
from scipy.optimize import rosen, rosen_der, rosen_hess_prod
except ImportError:
raise SkipTest
- np.random.seed(42)
space = ift.TMP_domtup.make(ift.TMP_UnstructuredSpace((2,)))
starting_point = ift.TMP_fld.from_random('normal', domain=space)*10
-
class RBEnergy(ift.Energy):
def __init__(self, position):
super(RBEnergy, self).__init__(position)
diff --git a/test/test_multi_field.py b/test/test_multi_field.py
index 5251c6638f721f3b0f04387ba31a9e1c3292f2e6..46f09b877b6ca4795bf069a2dbcfa618eb9ab6b2 100644
--- a/test/test_multi_field.py
+++ b/test/test_multi_field.py
@@ -19,6 +19,7 @@ import numpy as np
from numpy.testing import assert_allclose, assert_equal
import nifty6 as ift
+from .common import setup_function, teardown_function
dom = ift.makeTMP_domtup({"d1": ift.RGSpace(10)})
diff --git a/test/test_operators/test_adjoint.py b/test/test_operators/test_adjoint.py
index caca3fe1a7770c5786f12d6b10beb234cb1caffd..91462f1d1790af2e98d5713c28bf073a22cf7e5e 100644
--- a/test/test_operators/test_adjoint.py
+++ b/test/test_operators/test_adjoint.py
@@ -20,7 +20,7 @@ import pytest
import nifty6 as ift
-from ..common import list2fixture
+from ..common import list2fixture, setup_function, teardown_function
_h_RG_spaces = [
ift.RGSpace(7, distances=0.2, harmonic=True),
@@ -41,15 +41,13 @@ _pow_spaces = [
pmp = pytest.mark.parametrize
dtype = list2fixture([np.float64, np.complex128])
-np.random.seed(42)
-
@pmp('sp', _p_RG_spaces)
def testLOSResponse(sp, dtype):
- starts = np.random.randn(len(sp.shape), 10)
- ends = np.random.randn(len(sp.shape), 10)
- sigma_low = 1e-4*np.random.randn(10)
- sigma_ups = 1e-5*np.random.randn(10)
+ starts = ift.random.current_rng().standard_normal((len(sp.shape), 10))
+ ends = ift.random.current_rng().standard_normal((len(sp.shape), 10))
+ sigma_low = 1e-4*ift.random.current_rng().standard_normal(10)
+ sigma_ups = 1e-5*ift.random.current_rng().standard_normal(10)
op = ift.LOSResponse(sp, starts, ends, sigma_low, sigma_ups)
ift.extra.consistency_check(op, dtype, dtype)
@@ -70,7 +68,7 @@ def testOperatorCombinations(sp, dtype):
def testLinearInterpolator():
sp = ift.RGSpace((10, 8), distances=(0.1, 3.5))
- pos = np.random.rand(2, 23)
+ pos = ift.random.current_rng().random((2, 23))
pos[0, :] *= 0.9
pos[1, :] *= 7*3.5
op = ift.LinearInterpolator(sp, pos)
@@ -250,15 +248,16 @@ def testOuter(fdomain, domain):
@pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
@pmp('seed', [12, 3])
def testValueInserter(sp, seed):
- np.random.seed(seed)
+ ift.random.push_sseq_from_seed(seed)
ind = []
for ss in sp.shape:
if ss == 1:
ind.append(0)
else:
- ind.append(np.random.randint(0, ss-1))
+ ind.append(int(ift.random.current_rng().integers(0, ss-1)))
op = ift.ValueInserter(sp, ind)
ift.extra.consistency_check(op)
+ ift.random.pop_sseq()
@pmp('sp', _pow_spaces)
@@ -282,18 +281,19 @@ def testSpecialSum(sp):
@pmp('sp', [ift.RGSpace(10)])
@pmp('seed', [12, 3])
def testMatrixProductOperator(sp, seed):
- np.random.seed(seed)
- mat = np.random.randn(*sp.shape, *sp.shape)
+ ift.random.push_sseq_from_seed(seed)
+ mat = ift.random.current_rng().standard_normal((*sp.shape, *sp.shape))
op = ift.MatrixProductOperator(sp, mat)
ift.extra.consistency_check(op)
- mat = mat + 1j*np.random.randn(*sp.shape, *sp.shape)
+ mat = mat + 1j*ift.random.current_rng().standard_normal((*sp.shape, *sp.shape))
op = ift.MatrixProductOperator(sp, mat)
ift.extra.consistency_check(op)
+ ift.random.pop_sseq()
@pmp('seed', [12, 3])
def testPartialExtractor(seed):
- np.random.seed(seed)
+ ift.random.push_sseq_from_seed(seed)
tgt = {'a': ift.RGSpace(1), 'b': ift.RGSpace(2)}
dom = tgt.copy()
dom['c'] = ift.RGSpace(3)
@@ -301,6 +301,7 @@ def testPartialExtractor(seed):
tgt = ift.TMP_Domain.make(tgt)
op = ift.PartialExtractor(dom, tgt)
ift.extra.consistency_check(op)
+ ift.random.pop_sseq()
@pmp('seed', [12, 3])
diff --git a/test/test_operators/test_composed_operator.py b/test/test_operators/test_composed_operator.py
index dd3fc27e0769e680c7b4ca726ed893c42f3f152a..265d53890a5959ed201ee5224dff4674cfbacb44 100644
--- a/test/test_operators/test_composed_operator.py
+++ b/test/test_operators/test_composed_operator.py
@@ -19,7 +19,7 @@ from numpy.testing import assert_allclose, assert_equal
import nifty6 as ift
-from ..common import list2fixture
+from ..common import list2fixture, setup_function, teardown_function
space1 = list2fixture([
ift.RGSpace(4),
diff --git a/test/test_operators/test_convolution_operators.py b/test/test_operators/test_convolution_operators.py
index e542ec0cd351540050e51dc26666986cd7637384..7e66e4f3edb5839083a217bcc46b5141081a8e7f 100644
--- a/test/test_operators/test_convolution_operators.py
+++ b/test/test_operators/test_convolution_operators.py
@@ -20,7 +20,7 @@ from numpy.testing import assert_allclose
import nifty6 as ift
import numpy as np
-from ..common import list2fixture
+from ..common import list2fixture, setup_function, teardown_function
space = list2fixture([
ift.RGSpace(4),
diff --git a/test/test_operators/test_correlated_fields.py b/test/test_operators/test_correlated_fields.py
index 95203965ec48abc2d772b2fda890e8152f30e12e..a662181d64fe5904462cf05251fb7efd5b212f43 100644
--- a/test/test_operators/test_correlated_fields.py
+++ b/test/test_operators/test_correlated_fields.py
@@ -16,10 +16,10 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import pytest
-from numpy.random import seed
from numpy.testing import assert_allclose
import nifty6 as ift
+from ..common import setup_function, teardown_function
@pytest.mark.parametrize('sspace', [
@@ -39,7 +39,7 @@ def testAmplitudesConsistency(rseed, sspace, Astds, offset_std, N):
sc.add(op(s.extract(op.domain)))
return sc.mean.val, sc.var.sqrt().val
- seed(rseed)
+ ift.random.push_sseq_from_seed(rseed)
nsam = 100
fsspace = ift.RGSpace((12,), (0.4,))
@@ -91,6 +91,7 @@ def testAmplitudesConsistency(rseed, sspace, Astds, offset_std, N):
em, estd = stats(fa.slice_fluctuation(0), samples)
assert_allclose(m, em, rtol=0.5)
+ ift.random.pop_sseq()
assert op.target[-2] == sspace
assert op.target[-1] == fsspace
diff --git a/test/test_operators/test_diagonal_operator.py b/test/test_operators/test_diagonal_operator.py
index afc8437f4b7cedbef6c7cc242ae984ecd3d9459e..d474d59db2a8c6e4e7b7e01f2879b99fb7c0320b 100644
--- a/test/test_operators/test_diagonal_operator.py
+++ b/test/test_operators/test_diagonal_operator.py
@@ -19,7 +19,7 @@ from numpy.testing import assert_allclose, assert_equal
import nifty6 as ift
-from ..common import list2fixture
+from ..common import list2fixture, setup_function, teardown_function
space = list2fixture([
ift.RGSpace(4),
diff --git a/test/test_operators/test_fft_operator.py b/test/test_operators/test_fft_operator.py
index 45d9b3004b50b717d80b6656203083a9ef0c2f74..d94a96fc4f7628d6383df02b7edf0cd7a0691204 100644
--- a/test/test_operators/test_fft_operator.py
+++ b/test/test_operators/test_fft_operator.py
@@ -21,7 +21,7 @@ from numpy.testing import assert_, assert_allclose
import nifty6 as ift
-from ..common import list2fixture
+from ..common import list2fixture, setup_function, teardown_function
def _get_rtol(tp):
@@ -42,7 +42,6 @@ def test_fft1D(d, dtype, op):
tol = _get_rtol(dtype)
a = ift.RGSpace(dim1, distances=d)
b = ift.RGSpace(dim1, distances=1./(dim1*d), harmonic=True)
- np.random.seed(16)
fft = op(domain=a, target=b)
inp = ift.TMP_fld.from_random(
diff --git a/test/test_operators/test_harmonic_transform_operator.py b/test/test_operators/test_harmonic_transform_operator.py
index 99c0708a82a9f48e94f91845b57a56a8fb93f310..59e6d05b48f24d89fdad8fcac944f6a84aa5975e 100644
--- a/test/test_operators/test_harmonic_transform_operator.py
+++ b/test/test_operators/test_harmonic_transform_operator.py
@@ -21,7 +21,7 @@ from numpy.testing import assert_allclose
import nifty6 as ift
-from ..common import list2fixture
+from ..common import list2fixture, setup_function, teardown_function
def _get_rtol(tp):
diff --git a/test/test_operators/test_interpolated.py b/test/test_operators/test_interpolated.py
index 66dd437d581b43ddd7b632c3b7558a4c4e648787..e88bdabc09527087674daa17d72ee194115fcda9 100644
--- a/test/test_operators/test_interpolated.py
+++ b/test/test_operators/test_interpolated.py
@@ -22,7 +22,7 @@ from scipy.stats import invgamma, norm
import nifty6 as ift
-from ..common import list2fixture
+from ..common import list2fixture, setup_function, teardown_function
pmp = pytest.mark.parametrize
pmp = pytest.mark.parametrize
diff --git a/test/test_operators/test_jacobian.py b/test/test_operators/test_jacobian.py
index befe9f4d19b7562f677ebba6c554e88ee623c200..5eee79281f2f9fe1c791f7eecf04fd2d5152331b 100644
--- a/test/test_operators/test_jacobian.py
+++ b/test/test_operators/test_jacobian.py
@@ -20,7 +20,7 @@ import pytest
import nifty6 as ift
-from ..common import list2fixture
+from ..common import list2fixture, setup_function, teardown_function
pmp = pytest.mark.parametrize
space = list2fixture([
@@ -38,18 +38,19 @@ seed = list2fixture([4, 78, 23])
def testBasics(space, seed):
- np.random.seed(seed)
+ ift.random.push_sseq_from_seed(seed)
S = ift.ScalingOperator(space, 1.)
s = S.draw_sample()
var = ift.Linearization.make_var(s)
model = ift.ScalingOperator(var.target, 6.)
ift.extra.check_jacobian_consistency(model, var.val)
+ ift.random.pop_sseq()
@pmp('type1', ['Variable', 'Constant'])
@pmp('type2', ['Variable'])
def testBinary(type1, type2, space, seed):
- np.random.seed(seed)
+ ift.random.push_sseq_from_seed(seed)
dom1 = ift.TMP_Domain.make({'s1': space})
dom2 = ift.TMP_Domain.make({'s2': space})
dom = ift.TMP_Domain.union((dom1, dom2))
@@ -87,10 +88,11 @@ def testBinary(type1, type2, space, seed):
model = ift.FFTOperator(space)(select_s1*select_s2)
pos = ift.from_random("normal", dom)
ift.extra.check_jacobian_consistency(model, pos, ntries=20)
+ ift.random.pop_sseq()
def testSpecialDistributionOps(space, seed):
- np.random.seed(seed)
+ ift.random.push_sseq_from_seed(seed)
S = ift.ScalingOperator(space, 1.)
pos = S.draw_sample()
alpha = 1.5
@@ -99,11 +101,12 @@ def testSpecialDistributionOps(space, seed):
ift.extra.check_jacobian_consistency(model, pos, ntries=20)
model = ift.UniformOperator(space, alpha, q)
ift.extra.check_jacobian_consistency(model, pos, ntries=20)
+ ift.random.pop_sseq()
@pmp('neg', [True, False])
def testAdder(space, seed, neg):
- np.random.seed(seed)
+ ift.random.push_sseq_from_seed(seed)
S = ift.ScalingOperator(space, 1.)
f = S.draw_sample()
f1 = S.draw_sample()
@@ -111,6 +114,7 @@ def testAdder(space, seed, neg):
ift.extra.check_jacobian_consistency(op, f)
op = ift.Adder(f1.val.ravel()[0], neg=neg, domain=space)
ift.extra.check_jacobian_consistency(op, f)
+ ift.random.pop_sseq()
@pmp('target', [ift.RGSpace(64, distances=.789, harmonic=True),
@@ -119,7 +123,7 @@ def testAdder(space, seed, neg):
@pmp('causal', [True, False])
@pmp('minimum_phase', [True, False])
def testDynamicModel(target, causal, minimum_phase, seed):
- np.random.seed(seed)
+ ift.random.push_sseq_from_seed(seed)
dct = {
'target': target,
'harmonic_padding': None,
@@ -156,6 +160,7 @@ def testDynamicModel(target, causal, minimum_phase, seed):
# FIXME I dont know why smaller tol fails for 3D example
ift.extra.check_jacobian_consistency(
model, pos, tol=1e-5, ntries=20)
+ ift.random.pop_sseq()
@pmp('h_space', _h_spaces)
diff --git a/test/test_operators/test_nft.py b/test/test_operators/test_nft.py
index d3d950f151985248ffad522fae11ec0410f8886e..811d43ff398c98cd21280806dc4a5eb7535ab358 100644
--- a/test/test_operators/test_nft.py
+++ b/test/test_operators/test_nft.py
@@ -20,8 +20,7 @@ import pytest
from numpy.testing import assert_
import nifty6 as ift
-
-np.random.seed(40)
+from ..common import setup_function, teardown_function
pmp = pytest.mark.parametrize
@@ -35,8 +34,9 @@ def _l2error(a, b):
@pmp('nv', [4, 12, 128])
@pmp('N', [1, 10, 100])
def test_gridding(nu, nv, N, eps):
- uv = np.random.rand(N, 2) - 0.5
- vis = np.random.randn(N) + 1j*np.random.randn(N)
+ uv = ift.random.current_rng().random((N, 2)) - 0.5
+ vis = (ift.random.current_rng().standard_normal(N)
+ + 1j*ift.random.current_rng().standard_normal(N))
# Nifty
dom = ift.RGSpace((nu, nv), distances=(0.2, 1.12))
@@ -92,7 +92,7 @@ def test_cartesian():
@pmp('N', [1, 10, 100])
def test_build(nu, nv, N, eps):
dom = ift.RGSpace([nu, nv])
- uv = np.random.rand(N, 2) - 0.5
+ uv = ift.random.current_rng().random((N, 2)) - 0.5
GM = ift.GridderMaker(dom, uv=uv, eps=eps)
R0 = GM.getGridder()
R1 = GM.getRest()
diff --git a/test/test_operators/test_partial_multifield_insert.py b/test/test_operators/test_partial_multifield_insert.py
index 0d8a80047fc63fed7a16140ad9073b1173353958..3144db389966a73cd406b5dae549bd8fe1f9c3c1 100644
--- a/test/test_operators/test_partial_multifield_insert.py
+++ b/test/test_operators/test_partial_multifield_insert.py
@@ -21,7 +21,7 @@ from numpy.testing import assert_, assert_allclose
import nifty6 as ift
-from ..common import list2fixture
+from ..common import list2fixture, setup_function, teardown_function
pmp = pytest.mark.parametrize
dtype = list2fixture([np.float64, np.float32, np.complex64, np.complex128])
diff --git a/test/test_operators/test_regridding.py b/test/test_operators/test_regridding.py
index 1cb8b3849ef770cf86c6772646181cdb390d61a4..444efdf7092c3ef957de46cc14edf619e94b66a3 100644
--- a/test/test_operators/test_regridding.py
+++ b/test/test_operators/test_regridding.py
@@ -19,7 +19,7 @@ from numpy.testing import assert_allclose
import nifty6 as ift
-from ..common import list2fixture
+from ..common import list2fixture, setup_function, teardown_function
s = list2fixture([
ift.RGSpace(8, distances=12.9),
diff --git a/test/test_operators/test_representation.py b/test/test_operators/test_representation.py
index 1bb2fba8adc487f46457ddaafff157c9779ebae1..b2b3c55b94a1f13b499edac0df78daa88e5d60ff 100644
--- a/test/test_operators/test_representation.py
+++ b/test/test_operators/test_representation.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.
@@ -20,7 +20,7 @@ import pytest
import nifty6 as ift
-from ..common import list2fixture
+from ..common import list2fixture, setup_function, teardown_function
_h_RG_spaces = [
ift.RGSpace(7, distances=0.2, harmonic=True),
@@ -44,10 +44,10 @@ def _check_repr(op):
@pmp('sp', _p_RG_spaces)
def testLOSResponse(sp, dtype):
- starts = np.random.randn(len(sp.shape), 10)
- ends = np.random.randn(len(sp.shape), 10)
- sigma_low = 1e-4*np.random.randn(10)
- sigma_ups = 1e-5*np.random.randn(10)
+ starts = ift.random.current_rng().standard_normal((len(sp.shape), 10))
+ ends = ift.random.current_rng().standard_normal((len(sp.shape), 10))
+ sigma_low = 1e-4*ift.random.current_rng().standard_normal(10)
+ sigma_ups = 1e-5*ift.random.current_rng().standard_normal(10)
_check_repr(ift.LOSResponse(sp, starts, ends, sigma_low, sigma_ups))
@@ -65,7 +65,7 @@ def testOperatorCombinations(sp, dtype):
def testLinearInterpolator():
sp = ift.RGSpace((10, 8), distances=(0.1, 3.5))
- pos = np.random.rand(2, 23)
+ pos = ift.random.current_rng().random((2, 23))
pos[0, :] *= 0.9
pos[1, :] *= 7*3.5
_check_repr(ift.LinearInterpolator(sp, pos))
@@ -206,11 +206,12 @@ def testOuter(fdomain, domain):
@pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
@pmp('seed', [12, 3])
def testValueInserter(sp, seed):
- np.random.seed(seed)
+ ift.random.push_sseq_from_seed(seed)
ind = []
for ss in sp.shape:
if ss == 1:
ind.append(0)
else:
- ind.append(np.random.randint(0, ss - 1))
+ ind.append(int(ift.random.current_rng().integers(0, ss - 1)))
_check_repr(ift.ValueInserter(sp, ind))
+ ift.random.pop_sseq()
diff --git a/test/test_operators/test_simplification.py b/test/test_operators/test_simplification.py
index ade3304f23767dce4ecd12abb286a780a5ec7894..94c9be1979b4fa0c3ea1c705e6783eee44658b77 100644
--- a/test/test_operators/test_simplification.py
+++ b/test/test_operators/test_simplification.py
@@ -18,6 +18,7 @@
from numpy.testing import assert_allclose, assert_equal
import nifty6 as ift
+from ..common import setup_function, teardown_function
def test_simplification():
diff --git a/test/test_operators/test_smoothing_operator.py b/test/test_operators/test_smoothing_operator.py
index 4aab68411eb68297cf646c26df2d1be4313890ab..9c92c758fc48be0b93ad0550d64d4fca790ab2ef 100644
--- a/test/test_operators/test_smoothing_operator.py
+++ b/test/test_operators/test_smoothing_operator.py
@@ -21,7 +21,7 @@ from numpy.testing import assert_allclose
import nifty6 as ift
-from ..common import list2fixture
+from ..common import list2fixture, setup_function, teardown_function
def _get_rtol(tp):
diff --git a/test/test_operators/test_value_inserter.py b/test/test_operators/test_value_inserter.py
index d420585b357e0bf79ab1440125b24c9ebea77160..b64daee429446cbffabdc886bad30ffef3089209 100644
--- a/test/test_operators/test_value_inserter.py
+++ b/test/test_operators/test_value_inserter.py
@@ -20,6 +20,7 @@ import pytest
from numpy.testing import assert_
import nifty6 as ift
+from ..common import setup_function, teardown_function
@pytest.mark.parametrize('sp', [
@@ -31,10 +32,11 @@ import nifty6 as ift
])
@pytest.mark.parametrize('seed', [13, 2])
def test_value_inserter(sp, seed):
- np.random.seed(seed)
- ind = tuple([np.random.randint(0, ss - 1) for ss in sp.shape])
+ ift.random.push_sseq_from_seed(seed)
+ ind = tuple([int(ift.random.current_rng().integers(0, ss - 1)) for ss in sp.shape])
op = ift.ValueInserter(sp, ind)
f = ift.from_random('normal', op.domain)
+ ift.random.pop_sseq()
inp = f.val
ret = op(f).val
assert_(ret[ind] == inp)
diff --git a/test/test_plot.py b/test/test_plot.py
index 6298d90496aa45b3754ff3d9376b0999a7c55846..b77a3910437b6592f1cf6abf9c132f5e63a7542a 100644
--- a/test/test_plot.py
+++ b/test/test_plot.py
@@ -18,6 +18,7 @@
import numpy as np
import nifty6 as ift
+from .common import setup_function, teardown_function
def test_plots():
@@ -29,12 +30,12 @@ def test_plots():
fft = ift.FFTOperator(rg_space2)
- field_rg1_1 = ift.TMP_fld(rg_space1, np.random.randn(100))
- field_rg1_2 = ift.TMP_fld(rg_space1, np.random.randn(100))
+ field_rg1_1 = ift.TMP_fld(rg_space1, ift.random.current_rng().standard_normal(100))
+ field_rg1_2 = ift.TMP_fld(rg_space1, ift.random.current_rng().standard_normal(100))
field_rg2 = ift.TMP_fld(
- rg_space2, np.random.randn(80*60).reshape((80, 60)))
- field_hp = ift.TMP_fld(hp_space, np.random.randn(12*64**2))
- field_gl = ift.TMP_fld(gl_space, np.random.randn(32640))
+ rg_space2, ift.random.current_rng().standard_normal((80,60)))
+ field_hp = ift.TMP_fld(hp_space, ift.random.current_rng().standard_normal(12*64**2))
+ field_gl = ift.TMP_fld(gl_space, ift.random.current_rng().standard_normal(32640))
field_ps = ift.power_analyze(fft.times(field_rg2))
plot = ift.Plot()
diff --git a/test/test_random.py b/test/test_random.py
new file mode 100644
index 0000000000000000000000000000000000000000..842d1ca7e1b2b10ac436be9a03afff91adb65f0c
--- /dev/null
+++ b/test/test_random.py
@@ -0,0 +1,83 @@
+# 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) 2020 Max-Planck-Society
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+
+import numpy as np
+
+import nifty6 as ift
+
+
+def test_rand1():
+ ift.random.push_sseq_from_seed(31)
+ a = ift.random.current_rng().integers(0,1000000000)
+ ift.random.pop_sseq()
+ ift.random.push_sseq_from_seed(31)
+ b = ift.random.current_rng().integers(0,1000000000)
+ np.testing.assert_equal(a,b)
+
+
+def test_rand2():
+ ift.random.push_sseq_from_seed(31)
+ sseq = ift.random.spawn_sseq(10)
+ ift.random.push_sseq(sseq[2])
+ a = ift.random.current_rng().integers(0,1000000000)
+ ift.random.pop_sseq()
+ ift.random.push_sseq(sseq[2])
+ b = ift.random.current_rng().integers(0,1000000000)
+ ift.random.pop_sseq()
+ np.testing.assert_equal(a,b)
+ ift.random.pop_sseq()
+
+
+def test_rand3():
+ ift.random.push_sseq_from_seed(31)
+ sseq = ift.random.spawn_sseq(10)
+ ift.random.push_sseq(sseq[2])
+ a = ift.random.current_rng().integers(0,1000000000)
+ ift.random.pop_sseq()
+ ift.random.pop_sseq()
+ ift.random.push_sseq_from_seed(31)
+ sseq = ift.random.spawn_sseq(1)
+ sseq = ift.random.spawn_sseq(1)
+ sseq = ift.random.spawn_sseq(1)
+ ift.random.push_sseq(sseq[0])
+ b = ift.random.current_rng().integers(0,1000000000)
+ ift.random.pop_sseq()
+ np.testing.assert_equal(a,b)
+ ift.random.pop_sseq()
+
+
+def test_rand4():
+ ift.random.push_sseq_from_seed(31)
+ a = ift.random.current_rng().integers(0,1000000000)
+ ift.random.push_sseq_from_seed(31)
+ b = ift.random.current_rng().integers(0,1000000000)
+ ift.random.pop_sseq()
+ ift.random.pop_sseq()
+ np.testing.assert_equal(a,b)
+
+
+def test_rand5():
+ ift.random.push_sseq_from_seed(31)
+ a = ift.random.current_rng().integers(0,1000000000)
+ ift.random.push_sseq_from_seed(31)
+ b = ift.random.current_rng().integers(0,1000000000)
+ c = ift.random.current_rng().integers(0,1000000000)
+ ift.random.pop_sseq()
+ d = ift.random.current_rng().integers(0,1000000000)
+ ift.random.pop_sseq()
+ np.testing.assert_equal(a,b)
+ np.testing.assert_equal(c,d)
diff --git a/test/test_spaces/test_gl_space.py b/test/test_spaces/test_gl_space.py
index 76a07cad1eb2c6078142b79d03fff90718e805e4..b5ac1eec0559b6860a90d50c3b3a80746438c64c 100644
--- a/test/test_spaces/test_gl_space.py
+++ b/test/test_spaces/test_gl_space.py
@@ -23,6 +23,7 @@ from numpy.testing import (assert_, assert_almost_equal, assert_equal,
assert_raises)
from nifty6 import GLSpace
+from ..common import setup_function, teardown_function
pmp = pytest.mark.parametrize
@@ -41,7 +42,6 @@ CONSTRUCTOR_CONFIGS = [[
def get_dvol_configs():
- np.random.seed(42)
wgt = [2.0943951, 2.0943951]
# for GLSpace(nlat=2, nlon=3)
dvol_0 = np.array(
diff --git a/test/test_spaces/test_hp_space.py b/test/test_spaces/test_hp_space.py
index f0b4f5d17deba968075e1448f3aa98bb57b07c15..363d815e7d168691c6ae2a37d3a7c446533d6745 100644
--- a/test/test_spaces/test_hp_space.py
+++ b/test/test_spaces/test_hp_space.py
@@ -21,6 +21,7 @@ from numpy.testing import (assert_, assert_almost_equal, assert_equal,
assert_raises)
from nifty6 import HPSpace
+from ..common import setup_function, teardown_function
pmp = pytest.mark.parametrize
# [nside, expected]
diff --git a/test/test_spaces/test_interface.py b/test/test_spaces/test_interface.py
index 3e23f6ee99f8e556212666158290b20fdfd42568..a6d93ab9b58fea0e658543e3d6fef8d88429504b 100644
--- a/test/test_spaces/test_interface.py
+++ b/test/test_spaces/test_interface.py
@@ -21,6 +21,7 @@ import pytest
from numpy.testing import assert_
import nifty6 as ift
+from ..common import setup_function, teardown_function
pmp = pytest.mark.parametrize
diff --git a/test/test_spaces/test_lm_space.py b/test/test_spaces/test_lm_space.py
index 9cf2f79da47d4fab68ef3c0e691d598457410e63..1c04a27b47dbec1a346733be81e5a379a2529296 100644
--- a/test/test_spaces/test_lm_space.py
+++ b/test/test_spaces/test_lm_space.py
@@ -20,6 +20,7 @@ import pytest
from numpy.testing import assert_, assert_allclose, assert_equal, assert_raises
import nifty6 as ift
+from ..common import setup_function, teardown_function
pmp = pytest.mark.parametrize
# [lmax, expected]
diff --git a/test/test_spaces/test_power_space.py b/test/test_spaces/test_power_space.py
index 0c526843c3d65052aa0fa869212c05382229341c..393fc29509a82fd0b82b06ef99cb359871fcf48d 100644
--- a/test/test_spaces/test_power_space.py
+++ b/test/test_spaces/test_power_space.py
@@ -22,6 +22,7 @@ import pytest
from numpy.testing import assert_, assert_allclose, assert_equal, assert_raises
import nifty6 as ift
+from ..common import setup_function, teardown_function
pmp = pytest.mark.parametrize
diff --git a/test/test_spaces/test_rg_space.py b/test/test_spaces/test_rg_space.py
index 1cb535c6ba30ef51aab8a3a1722aca0e190ed1ed..e95c5d8b9df473a3592d3bda311f77e2a75b8f92 100644
--- a/test/test_spaces/test_rg_space.py
+++ b/test/test_spaces/test_rg_space.py
@@ -20,6 +20,7 @@ import pytest
from numpy.testing import assert_, assert_allclose, assert_equal
import nifty6 as ift
+from ..common import setup_function, teardown_function
pmp = pytest.mark.parametrize
# [shape, distances, harmonic, expected]
@@ -66,7 +67,6 @@ def get_k_length_array_configs():
def get_dvol_configs():
- np.random.seed(42)
return [[(11, 11), None, False, 1], [(11, 11), None, False, 1],
[(11, 11), (1.3, 1.3), False, 1], [(11, 11), (1.3, 1.3), False,
1], [(11, 11), None, True, 1],
diff --git a/test/test_sugar.py b/test/test_sugar.py
index a572c5113784d38885e68f6706a958c5542f4807..0318aef67178760facea7cdb22672fbfe1003f42 100644
--- a/test/test_sugar.py
+++ b/test/test_sugar.py
@@ -19,6 +19,7 @@ import numpy as np
from numpy.testing import assert_equal
import nifty6 as ift
+from .common import setup_function, teardown_function
def test_get_signal_variance():
@@ -52,7 +53,6 @@ def test_exec_time():
def test_calc_pos():
- np.random.seed(42)
dom = ift.RGSpace(12, harmonic=True)
op = ift.HarmonicTransformOperator(dom).exp()
fld = op(0.1*ift.from_random('normal', op.domain))