Repository: sdepy/sdepy

Issue:
Use a default_rng() instance instead of random.normal
Feature request:
wiener source is using np.random.normal and np.random.multivariate_normal.

Both need np.random.seed(SEED) to be set for tests or repeatable experiments.
random.seed is known to not work as expected in multithreaded environments (some test frameworks, REST-Implementations).
Numpy strongly recommends to use an instance of a random number generator instead:
https://numpy.org/doc/stable/reference/random/generated/numpy.random.normal.html?highlight=random%20normal#numpy.random.normal

I'd like to have the opportunity to pass a random generator as optional parameter to the class which is then used as random source.

Discussion:

Patch:
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index f80437d..175260a 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -43,8 +43,8 @@ jobs:
       fail-fast: true
       matrix:
         include:
-          - name: Fast python 3.5 tests with minimal requirements
-            py: 3.5
+          - name: Fast python 3.6 tests with minimal requirements
+            py: 3.6
             dependencies: numpy==1.15.2 scipy==0.19.1 pytest==3.8.1
             fast: true
             full: false
diff --git a/doc/quickguide.rst b/doc/quickguide.rst
index 392b2e3..dd5651a 100644
--- a/doc/quickguide.rst
+++ b/doc/quickguide.rst
@@ -65,7 +65,6 @@ It is best explained by examples, involving ``my_process``, ``myp`` and
 1. **Scalar process** in 100000 paths, with default parameters, computed
    at 5 time points (``coarse_timeline``), using 100 steps in between::
 
-    >>> np.random.seed(1)  # make doctests predictable
     >>> x = my_process(x0=1, paths=100*1000,
     ...                steps=100)(coarse_timeline)
     >>> x.shape
@@ -168,7 +167,7 @@ It is best explained by examples, involving ``my_process``, ``myp`` and
     >>> for i in range(10):
     ...     gr = plt.plot(k[0, :, 0], x[-1, i, :, :].mean(axis=-1))
     >>> lb = plt.xlabel('k'), plt.ylabel('x(t=2).mean()')
-    >>> plt.show() # doctest: +SKIP
+    >>> plt.show()  # doctest: +SKIP
 
    In the example above, set ``steps>=100`` to go from inaccurate and fast,
    to meaningful and slow.
@@ -292,7 +291,7 @@ It is best explained by examples, involving ``my_process``, ``myp`` and
     >>> gr = plt.plot(t, x1(t)[:, i], t, x2(t)[:, i], t, x3(t)[:, i])
     >>> plt.show() # doctest: +SKIP
     >>> gr = plt.plot(t, y1[:, i], t, y2[:, i], t, y3[:, i])
-    >>> plt.show() # doctest: +SKIP
+    >>> plt.show()  # doctest: +SKIP
 
 
 ------------------------------------
@@ -504,7 +503,6 @@ an exact result irrespective of the number and size
 of the integration steps (this happens since, by implementation,
 it integrates the linear SDE for ``log(x)``)::
 
-    >>> np.random.seed(1)
     >>> args = dict(dw=sdepy.true_wiener_source(paths=100),
     ...             paths=100, x0=10)
     >>> timeline = (0, 1)
@@ -524,7 +522,8 @@ it integrates the linear SDE for ``log(x)``)::
     >>> ax = plt.axes(label=0); ax.set_xscale('log'); ax.set_yscale('log')
     >>> gr = ax.plot(steps, errors)
     >>> plt.show()  # doctest: +SKIP
-    >>> print('euler error: {:.2e}\n   rk error: {:.2e}'.format(errors[-1,0], errors[-1,1]))
+    >>> print(f'euler error: {errors[-1,0]:.2e}\n'
+    ...       f'   rk error: {errors[-1,1]:.2e}')  # doctest: +SKIP
     euler error: 1.70e-03
        rk error: 8.80e-06
 
@@ -547,7 +546,6 @@ Setup::
     >>> from scipy.special import erf
     >>> from scipy.integrate import quad
 
-    >>> np.random.seed(1)
     >>> k = .5
     >>> x0, x1 = 0, 10;
     >>> t0, t1 = 0, 1
@@ -608,9 +606,9 @@ Comparison against exact values::
     >>> gr = plt.plot(xgrid, u2, y, u2_exact(y, t1), ':')
     >>> plt.show()  # doctest: +SKIP
 
-    >>> print('u1 error: {:.2e}\nu2 error: {:.2e}'.format(
-    ...     np.abs(u1 - u1_exact(xgrid, t1)).mean(),
-    ...     np.abs(u2 - u2_exact(xgrid, t1)).mean()))
+    >>> print(f'u1 error: {np.abs(u1 - u1_exact(xgrid, t1)).mean():.2e}\n'
+    ...       f'u2 error: {np.abs(u2 - u2_exact(xgrid, t1)).mean():.2e}'
+    ... )  # doctest: +SKIP
     u1 error: 2.49e-03
     u2 error: 5.51e-03
 
@@ -653,7 +651,6 @@ simulated across 50000 independent paths and their antithetics
     ...     mu=(riskfree - dividend_yield),
     ...     sigma=sigma,
     ...     vshape=4, paths=100*1000, steps=maturity*250)
-    >>> np.random.seed(1)
     >>> x = p(timeline)
     >>> x.shape
     (9, 4, 100000)
diff --git a/sdepy/infrastructure.py b/sdepy/infrastructure.py
index 2d69924..d1b2d55 100644
--- a/sdepy/infrastructure.py
+++ b/sdepy/infrastructure.py
@@ -18,10 +18,36 @@ import inspect
 import warnings
 
 
+########################################
+#  Default random number generator
+########################################
+
+# To use as sdepy default the legacy numpy random generator,
+# seeded with SEED, set
+# >>> sdepy.infrastructure.default_rng = np.random.RandomState(SEED)
+# To use the legacy global state, seeded via `np.random.seed(SEED)`, set
+# >>> sdepy.infrastructure.default_rng = np.random
+try:
+    default_rng = np.random.default_rng()
+    _rng_types = (np.random.Generator, np.random.RandomState,
+                  type(None))
+
+except AttributeError:
+    # ensure compatibility with legacy numpy and scipy versions
+    default_rng = np.random
+    _rng_types = (type(None))
+
+
 ########################################
 #  Private functions for recurring tasks
 ########################################
 
+def _get_default_rng():
+    """Get current default random number generator, as stored in
+    sdepy.infrastructure.default_rng."""
+    return default_rng
+
+
 def _shape_setup(shape):
     """Array shape preprocessing, return (shape,) if shape is an integer."""
     return (shape,) if isinstance(shape, int) else shape
@@ -1282,6 +1308,10 @@ class source:
         Shape of source values.
     dtype : data-type
         Data type of source values. Defaults to ``None``.
+    rng : numpy.random.Generator, or numpy.random.RandomState, or None
+        Random numbers generator used. If ``None``, defaults to
+        ``sdepy.infrastructure.default_rng``, a global variabile
+        initialized on import to ``numpy.random.default_rng()``.
 
     Returns
     -------
@@ -1311,13 +1341,20 @@ class source:
 
     Attributes
     ----------
+    rng
     size
     t
     """
 
-    def __init__(self, *, paths=1, vshape=(), dtype=None):
-        self.paths, self.vshape, self.dtype = paths, vshape, dtype
-        self.vshape = _shape_setup(self.vshape)
+    def __init__(self, *, paths=1, vshape=(), dtype=None, rng=None):
+        self.paths, self.dtype = paths, dtype
+        self.vshape = _shape_setup(vshape)
+        if not isinstance(rng, _rng_types):
+            raise TypeError(
+                f'rng should be an instance of `numpy.random.Generator` '
+                f'or `numpy.random.RandomState`, not of {type(rng)}. '
+                f'For legacy numpy versions, only `rng=None` is allowed.')
+        self._rng = default_rng if (rng is None) else rng
 
     def __call__(self, t, dt=None):
         """Realization of stochasticity source values or increments."""
@@ -1325,6 +1362,16 @@ class source:
         t, dt = np.asarray(t), np.asarray(dt)
         return t + dt + np.nan
 
+    @property
+    def rng(self):
+        """Read-only access to the random number generator used
+        by the stochasticity source.
+        """
+        # prevent modifications of the `rng` attribute (for `true_source`
+        # subclasses, such changes would silently fail to propagate
+        # to sources stored in `_dw`, `_dn` or `_dj` attributes).
+        return self._rng
+
     @property
     def size(self):
         """
@@ -1358,6 +1405,10 @@ class wiener_source(source):
         Shape of source values.
     dtype : data-type
         Data type of source values. Defaults to ``None``.
+    rng : numpy.random.Generator, or numpy.random.RandomState, or None
+        Random numbers generator used. If ``None``, defaults to
+        ``sdepy.infrastructure.default_rng``, a global variabile
+        initialized on import to ``numpy.random.default_rng()``.
     corr : array-like, or callable, or None
         Correlation matrix of the standard Wiener process increments,
         possibly time-dependent, or ``None`` for no correlations,
@@ -1435,9 +1486,9 @@ class wiener_source(source):
     ``dt`` intervals.
     """
 
-    def __init__(self, *, paths=1, vshape=(), dtype=None,
+    def __init__(self, *, paths=1, vshape=(), dtype=None, rng=None,
                  corr=None, rho=None):
-        super().__init__(paths=paths, vshape=vshape, dtype=dtype)
+        super().__init__(paths=paths, vshape=vshape, dtype=dtype, rng=rng)
 
         # get the correlation matrix from 'corr' and 'rho'
         self.corr = corr = _get_corr_matrix(corr, rho)
@@ -1463,15 +1514,15 @@ class wiener_source(source):
         corr = self.corr
         t, dt = np.broadcast_arrays(t, dt)
         tshape = dt.shape
-        # using np.random.normal and np.random.multivariate_normal,
+        # using numpy normal and multivariate_normal,
         # instead of scipy.stats.norm and scipy.stats.multivariate_normal
         # to imporve speed (aviod overhead of scipy.stats random variable
         # instantiation at each call)
         if corr is None:
             # --- handle uncorrelated samples (vshape may be any shape)
-            dz = np.random.normal(
+            dz = self.rng.normal(
                 0., 1., size=tshape + vshape + (paths,)
-               ).astype(dtype, copy=False)
+                ).astype(dtype, copy=False)
         else:
             # --- handle correlated samples
             M = vshape[-1]
@@ -1486,7 +1537,7 @@ class wiener_source(source):
                             'invalid correlation matrix shape {}'
                             .format(cov.shape))
                     cov = cov[..., 0]  # remove paths axis if present
-                dz = np.random.multivariate_normal(
+                dz = self.rng.multivariate_normal(
                     mean=mean, cov=cov, size=tshape + vshape[:-1] + (paths,)
                     ).astype(dtype, copy=False)
             else:
@@ -1502,7 +1553,7 @@ class wiener_source(source):
                                 'invalid correlation matrix shape {}'
                                 .format(cov.shape))
                         cov = cov[..., 0]  # remove paths axis if present
-                    dz[i] = np.random.multivariate_normal(
+                    dz[i] = self.rng.multivariate_normal(
                         mean=mean, cov=cov, size=vshape[:-1] + (paths,)
                         )
             # reshape dz to ``tshape + vshape + (paths,)``
@@ -1531,6 +1582,10 @@ class poisson_source(source):
         Shape of source values.
     dtype : data-type
         Data type of source values. Defaults to ``int``.
+    rng : numpy.random.Generator, or numpy.random.RandomState, or None
+        Random numbers generator used. If ``None``, defaults to
+        ``sdepy.infrastructure.default_rng``, a global variabile
+        initialized on import to ``numpy.random.default_rng()``.
     lam : array-like, or callable
         Intensity of the Poisson process, possibly time-dependent.
         Should be an array of non-negative values, broadcastable to shape
@@ -1552,8 +1607,8 @@ class poisson_source(source):
     source
     """
 
-    def __init__(self, *, paths=1, vshape=(), dtype=int, lam=1.):
-        super().__init__(paths=paths, vshape=vshape, dtype=dtype)
+    def __init__(self, *, paths=1, vshape=(), dtype=int, rng=None, lam=1.):
+        super().__init__(paths=paths, vshape=vshape, dtype=dtype, rng=rng)
         self.lam = lam = _variable_param_setup(lam)
         lam_shape = _get_param_shape(lam)
         if lam_shape is not None:
@@ -1576,13 +1631,13 @@ class poisson_source(source):
         lam = self.lam
         dn = np.empty(tshape + vshape + (paths,),
                       dtype=dtype)
-        # using np.random.poisson instead of scipy.stats.poisson
+        # using numpy poisson instead of scipy.stats.poisson
         # to imporve speed (aviod overhead of scipy.stats random variable
         # instantiation at each call)
         for i in np.ndindex(tshape):
             L = (lam(t[i] + dt[i]/2) if callable(lam) else lam)
-            dn[i] = sign_dt[i]*np.random.poisson(abs_dt[i]*L,
-                                                 vshape + (paths,))
+            dn[i] = sign_dt[i]*self.rng.poisson(abs_dt[i]*L,
+                                                vshape + (paths,))
         return dn
 
 
@@ -1690,8 +1745,8 @@ class _exp_rv:
         self._s = np.sign(a)
         self._rv = scipy.stats.expon(scale=np.abs(a))
 
-    def rvs(self, size):
-        return self._s*self._rv.rvs(size)
+    def rvs(self, size, random_state=None):
+        return self._s*self._rv.rvs(size=size, random_state=random_state)
 
     def mean(self):
         return self._s*self._rv.mean()
@@ -1721,11 +1776,11 @@ class _double_exp_rv:
         self._rvxb = scipy.stats.expon(scale=b)
         self._rvu = scipy.stats.uniform(scale=1.)
 
-    def rvs(self, size):
+    def rvs(self, size, random_state=None):
         pa = self._pa
-        rvs_plus = self._rvxa.rvs(size)
-        rvs_minus = self._rvxb.rvs(size)
-        uniform = self._rvu.rvs(size)
+        rvs_plus = self._rvxa.rvs(size=size, random_state=random_state)
+        rvs_minus = self._rvxb.rvs(size=size, random_state=random_state)
+        uniform = self._rvu.rvs(size=size, random_state=random_state)
         return np.where(uniform <= pa, rvs_plus, -rvs_minus) + 0
 
     def mean(self):
@@ -1796,9 +1851,10 @@ def rvmap(f, y):
             yt = y(t) if callable(y) else y
 
             class new_yt_class:
-                def rvs(self, size):
-                    return (f(t, yt.rvs(size)) if time_dependent_f
-                            else f(yt.rvs(size)))
+                def rvs(self, size, random_state=None):
+                    yt_rvs = yt.rvs(size=size, random_state=random_state)
+                    return (f(t, yt_rvs) if time_dependent_f
+                            else f(yt_rvs))
 
             new_yt = new_yt_class()
             return new_yt
@@ -1806,8 +1862,8 @@ def rvmap(f, y):
     else:
 
         class new_y_class:
-            def rvs(self, size):
-                return f(y.rvs(size))
+            def rvs(self, size, random_state=None):
+                return f(y.rvs(size=size, random_state=random_state))
 
         new_y = new_y_class()
         return new_y
@@ -1842,9 +1898,17 @@ class cpoisson_source(source):
         Shape of source values.
     dtype : data-type
         Data type of source values. Defaults to ``None``.
-    dn : source or source class, optional
-        If given, ``dn`` is used as the underlying source of Poisson process
-        increments, overriding the ``ptype`` and ``lam`` parameters.
+    rng : numpy.random.Generator, or numpy.random.RandomState, or None
+        Random numbers generator used. If ``None``, defaults to
+        ``sdepy.infrastructure.default_rng``, a global variabile
+        initialized on import to ``numpy.random.default_rng()``.
+        Used to generate Poisson process increments (unless ``dn`` is
+        explicilty given with its own ``dn.rng``) and ``y`` random variates.
+    dn : source or source class, or None
+        Underlying source of Poisson process increments. If a class, it
+        is used to instantiate the source; if a source, it is used as it is,
+        overriding the given ``ptype`` and ``lam`` parameters.
+        If ``None``, it is instantiated as a ``sdepy.poisson_source``.
     ptype : data-type
         Data type of Poisson process increments. Defaults to ``int``.
     lam : array-like, or callable
@@ -1853,12 +1917,13 @@ class cpoisson_source(source):
     y : distribution, or callable, or None
         Distribution of random variates to be compounded with the
         Poisson process increments, possibly time-dependent.
-        May be any ``scipy.stats`` distribution instance,
-        or any object exposing an ``rvs(shape)`` method
-        to generate independent random variates of the given shape,
-        or a callable with ``y(t)`` evaluating to such object.
-        The following preset distributions may be specified, possibly
-        with time-varying parameters:
+        May be any ``scipy.stats`` distribution instance, or any object
+        exposing an ``rvs`` method to be invoked with signature
+        ``rvs(size=output_shape, random_state=rng)`` to generate independent
+        random variates with given shape and generator, or a callable
+        with ``y(t)`` evaluating to such object. The following
+        preset distributions may be specified, possibly with time-varying
+        parameters:
           -  ``y=norm_rv(a, b)`` - normal distribution with mean ``a``
              and standard deviation ``b``.
           -  ``y=uniform_rv(a, b)`` - uniform distribution
@@ -1926,15 +1991,15 @@ class cpoisson_source(source):
         ``y`` random variates.
     """
 
-    def __init__(self, *, paths=1, vshape=(), dtype=None,
+    def __init__(self, *, paths=1, vshape=(), dtype=None, rng=None,
                  dn=None, ptype=int, lam=1.,
                  y=None):
-        super().__init__(paths=paths, vshape=vshape, dtype=dtype)
+        super().__init__(paths=paths, vshape=vshape, dtype=dtype, rng=rng)
 
         # setup of poisson source
         self.dn = _source_setup(dn, poisson_source,
-                                paths=self.paths, vshape=self.vshape,
-                                dtype=ptype, lam=lam)
+                                paths=paths, vshape=vshape,
+                                dtype=ptype, rng=rng, lam=lam)
 
         # mind not breaking the source protocol
         self.ptype = self.dn.dtype if hasattr(dn, 'dtype') else ptype
@@ -1943,6 +2008,20 @@ class cpoisson_source(source):
         # setup of random variable sampling source
         self.y = uniform_rv(a=0, b=1) if y is None else y
 
+    @staticmethod
+    def _y_rvs(rvs, size, rng):
+        """Wrap calls to `rvs` method of `y` for backward compatibility"""
+        try:
+            return rvs(size=size, random_state=rng)
+        except TypeError:
+            rv_sample = rvs(size=size)
+            warnings.warn(
+                'The use of cpoission_source with distributions '
+                'not accepting a `random_state` keyword argument '
+                'is deprecated, and will not be supported in future releases.',
+                DeprecationWarning)
+            return rv_sample
+
     def __call__(self, t, dt):
         """See cpoisson_source class documentation."""
 
@@ -1961,7 +2040,7 @@ class cpoisson_source(source):
             for j in range(1, nmax+1):
                 index = (dn_positive_i == j)
                 if index.any():
-                    y_sample = rv.rvs(size=(index.sum(), j))
+                    y_sample = self._y_rvs(rv.rvs, (index.sum(), j), self.rng)
                     dz[i][index] = sign_dt[i]*y_sample.sum(axis=-1)
                     y_value.append(y_sample)
         self.dn_value = dn
@@ -1983,14 +2062,14 @@ def _antithetics(source_class, transform):
     """
 
     class antithetics_class(source):
-        def __init__(self, *, paths=2, vshape=(), **args):
-            self.paths = paths
-            self.vshape = _shape_setup(vshape)
+        def __init__(self, *, paths=2, vshape=(),
+                     dtype=None, rng=None, **args):
             if paths % 2:
                 raise ValueError(
                     'the number of paths for sources with antithetics '
                     'should be even, not {}'.format(paths))
-            self._dz = source_class(paths=paths//2, vshape=vshape, **args)
+            self._dz = source_class(paths=paths//2, vshape=vshape,
+                                    dtype=dtype, rng=rng, **args)
 
         __init__.__doc__ = ("See {} class documentation"
                             .format(source_class.__name__))
@@ -1999,6 +2078,22 @@ def _antithetics(source_class, transform):
             dz = self._dz(t, dt)
             return np.concatenate((dz, transform(dz)), axis=-1)
 
+        @property
+        def paths(self):
+            return 2*self._dz.paths
+
+        @property
+        def vshape(self):
+            return self._dz.vshape
+
+        @property
+        def dtype(self):
+            return self._dz.dtype
+
+        @property
+        def rng(self):
+            return self._dz.rng
+
     return antithetics_class
 
 
@@ -2098,7 +2193,7 @@ class true_source(source):
 
     Parameters
     ----------
-    paths, vshape, dtype
+    paths, vshape, dtype, rng
         See ``source`` class documentation.
     rtol : float, or 'max'
         relative tolerance used in assessing the coincidence
@@ -2131,10 +2226,10 @@ class true_source(source):
     new_outside
     """
 
-    def __init__(self, *, paths=1, vshape=(), dtype=None,
+    def __init__(self, *, paths=1, vshape=(), dtype=None, rng=None,
                  rtol='max', t0=0., z0=0.):
-        self.paths, self.vshape, self.dtype = paths, vshape, dtype
-        self.vshape = _shape_setup(self.vshape)
+        super().__init__(paths=paths, vshape=vshape, dtype=dtype, rng=rng)
+
         # time points are handled one by one,
         # their type is set to float
         self.rtol = np.finfo(float).resolution \
@@ -2359,16 +2454,16 @@ dw, source of standard Wiener process (brownian motion) increments with memory.
     See source and true_source methods.
     """
 
-    def __init__(self, *, paths=1, vshape=(), dtype=None,
+    def __init__(self, *, paths=1, vshape=(), dtype=None, rng=None,
                  corr=None, rho=None,
                  rtol='max', t0=0., z0=0.):
 
+        super().__init__(paths=paths, vshape=vshape, dtype=dtype, rng=rng,
+                         rtol=rtol, t0=t0, z0=z0)
         self._dw = wiener_source(paths=paths,
-                                 vshape=vshape, dtype=dtype,
+                                 vshape=vshape, dtype=dtype, rng=rng,
                                  corr=corr, rho=rho)
         self.corr = self._dw.corr
-        super().__init__(paths=paths, vshape=vshape, dtype=dtype,
-                         rtol=rtol, t0=t0, z0=z0)
 
     def new_outside(self, w, t, s):
         # approximate in case of time depentent correlations
@@ -2448,15 +2543,13 @@ class true_poisson_source(true_source):
     -------
     See ``source`` and ``true_source`` methods.
     """
-    def __init__(self, *, paths=1, vshape=(), dtype=int, lam=1.,
+    def __init__(self, *, paths=1, vshape=(), dtype=int, rng=None, lam=1.,
                  rtol='max', t0=0., z0=0):
 
-        super().__init__(paths=paths, vshape=vshape, dtype=dtype,
+        super().__init__(paths=paths, vshape=vshape, dtype=dtype, rng=rng,
                          rtol=rtol, t0=t0, z0=z0)
-        self._dn = poisson_source(paths=self.paths,
-                                  vshape=self.vshape,
-                                  dtype=self.dtype,
-                                  lam=lam)
+        self._dn = poisson_source(paths=paths, vshape=vshape, dtype=dtype,
+                                  rng=rng, lam=lam)
         self.lam = self._dn.lam
 
     def new_outside(self, n, t, s):
@@ -2467,7 +2560,7 @@ class true_poisson_source(true_source):
             return None
         p = (s - t1)/(t2 - t1)
         n = n2 - n1
-        return n1 + np.random.binomial(n, p, n.shape)
+        return n1 + self.rng.binomial(n, p, n.shape)
 
 
 class true_cpoisson_source(true_source):
@@ -2480,6 +2573,11 @@ class true_cpoisson_source(true_source):
         See ``cpoisson_source`` class documentation.
     rtol, t0, z0
         See ``true_source`` class documentation.
+    dn : true_poission_source, or None
+        If provided, it is used as the underlying source of Poisson process
+        increments with memory, overriding the given ``ptype`` and ``lam``.
+        If ``None`` (default), it is instantiated as
+        a ``sdepy.true_poission_source``.
 
     Returns
     -------
@@ -2507,18 +2605,18 @@ class true_cpoisson_source(true_source):
     -------
     See ``source`` and ``true_source`` methods.
     """
-    def __init__(self, *, paths=1, vshape=(), dtype=None,
+    def __init__(self, *, paths=1, vshape=(), dtype=None, rng=None,
                  rtol='max', t0=0., z0=0.,
                  dn=None, ptype=int, lam=1.,
                  y=None):
 
-        super().__init__(paths=paths, vshape=vshape, dtype=dtype,
+        super().__init__(paths=paths, vshape=vshape, dtype=dtype, rng=rng,
                          rtol=rtol, t0=t0, z0=z0)
         if dn is None:
             dn = true_poisson_source(paths=paths, vshape=vshape,
-                                     dtype=ptype, lam=lam)
-        self._dj = cpoisson_source(paths=self.paths,
-                                   vshape=self.vshape, dtype=self.dtype,
+                                     dtype=ptype, rng=rng, lam=lam)
+        self._dj = cpoisson_source(paths=paths, vshape=vshape,
+                                   dtype=dtype, rng=rng,
                                    dn=dn, y=y)  # ptype, lam are set by dn
         self.ptype, self.lam = self._dj.ptype, self._dj.lam
         self.dn = self._dj.dn
diff --git a/sdepy/integration.py b/sdepy/integration.py
index 8bcae83..ffe0a07 100644
--- a/sdepy/integration.py
+++ b/sdepy/integration.py
@@ -23,6 +23,7 @@ from .infrastructure import (
 ########################################
 
 from .infrastructure import (
+    _get_default_rng,
     _shape_setup,
     _const_param_setup,
     _variable_param_setup,
@@ -832,6 +833,11 @@ class SDE:
         Shape of the values of the process.
     dtype : data-type, optional
         Data-type of the process. Defaults to the numpy default.
+    rng : numpy.random.Generator, or numpy.random.RandomState, or None
+        Random numbers generator used to instantiate sources.
+        If ``None``, defaults to
+        ``sdepy.infrastructure.default_rng``, a global variabile
+        initialized on import to ``numpy.random.default_rng()``.
     steps : iterable, or int, or None
         Specification of the time points to be touched during integration
         (as accepted by a cooperating ``integrator`` class).
@@ -865,7 +871,8 @@ class SDE:
         of ``sde``, ``init`` and ``more`` methods, and stochasticity sources
         parameters, as implied by the signature of ``source_{id}`` methods.
         Each keyword should be used once (e.g. ``corr``, a ``source_dw``
-        parameter, cannot be used as the name of a SDE parameter).
+        parameter, should not be used as the name of a SDE parameter) or,
+        if repeated, must have consistent default values.
 
     Returns
     -------
@@ -1071,9 +1078,8 @@ class SDE:
     # initialization
     # --------------
 
-    def __init__(self, *, paths=1,
-                 vshape=(),
-                 dtype=None, steps=None, i0=0,
+    def __init__(self, *, paths=1, vshape=(), dtype=None, rng=None,
+                 steps=None, i0=0,
                  info=None, getinfo=True,
                  method='euler',
                  **args):
@@ -1115,6 +1121,14 @@ class SDE:
         self.paths = paths
         self.dtype = dtype
 
+        # set rng passed upon source instantiation
+        # (needed for compatibility with legacy numpy versions,
+        # for which sources accept only rng=None)
+        self._rng_asis = rng
+        # set rng accessed via the rng attribute (in case rng is None,
+        # should be the default_rng at the time of instantiation)
+        self._rng = _get_default_rng() if (rng is None) else rng
+
         # setup stochasticity sources
         # (this uses paths, shapes and dtype attributes)
         self.sources = self._sources_setup()
@@ -1222,6 +1236,15 @@ class SDE:
     # public interface vs subclasses
     # ------------------------------
 
+    @property
+    def rng(self):
+        """Read-only access to the random number generator used
+        to instantiate stochasticity sources.
+        """
+        # prevent modifications of the `rng` attribute (changes would silently
+        # fail to propagate to sources).
+        return self._rng
+
     @property
     def args(self):
         """
@@ -1294,7 +1317,8 @@ class SDE:
         ----------
         dw : source, or source subclass, or None
             If an object complying with the ``source`` protocol,
-            it is returned (``corr`` and ``rho`` are ignored).
+            it is returned (``corr`` and ``rho``, as well as
+            ``self.dtype`` and ``self.rng``, are ignored).
             If a source subclass, it is instantiated with the
             given parameters, and returned. If None, a
             new instance of ``wiener_source`` is returned, with
@@ -1305,7 +1329,10 @@ class SDE:
         -------
         An object complying with the ``source`` protocol,
         instantiating the requested stochasticity source.
-        The shape of source values is set to ``wshape``.
+        The shape and paths of source values are set to
+        ``self.paths``, ``self.wshape`` respectively.
+        If instantiated, ``dtype=self.dtype`` and ``rng=self.rng``
+        keywords are used upon instantiation.
 
         See Also
         --------
@@ -1313,7 +1340,7 @@ class SDE:
         """
         return _source_setup(dw, wiener_source,
                              paths=self.paths, vshape=self.wshape,
-                             dtype=self.dtype,
+                             dtype=self.dtype, rng=self._rng_asis,
                              corr=corr, rho=rho)
 
     def source_dn(self,
@@ -1326,7 +1353,8 @@ class SDE:
         ----------
         dn : source, or source subclass, or None
             If an object complying with the ``source`` protocol,
-            it is returned (``ptype`` and ``lam`` are ignored).
+            it is returned (``ptype`` and ``lam``, as well as
+            ``self.rng``, are ignored).
             If a source subclass, it is instantiated with the
             given parameters, and returned. If None, a
             new instance of ``poisson_source`` is returned, with
@@ -1337,7 +1365,10 @@ class SDE:
         -------
         An object complying with the ``source`` protocol,
         instantiating the requested stochasticity source.
-        The shape of source values is set to ``wshape``.
+        The shape and paths of source values are set to
+        ``self.paths``, ``self.wshape`` respectively.
+        If instantiated, ``rng=self.rng`` keyword is used
+        upon instantiation.
 
         See Also
         --------
@@ -1346,7 +1377,7 @@ class SDE:
         return _source_setup(dn, poisson_source,
                              paths=self.paths, vshape=self.wshape,
                              # dtype of the poisson source is ptype
-                             dtype=ptype,
+                             dtype=ptype, rng=self._rng_asis,
                              lam=lam)
 
     def source_dj(self,
@@ -1361,7 +1392,8 @@ class SDE:
         ----------
         dj : source, or source subclass, or None
             If an object complying with the ``source`` protocol,
-            it is returned (``ptype``, ``lam`` and ``y`` are ignored).
+            it is returned (``ptype``, ``lam`` and ``y``, as well as
+            ``self.dtype`` and ``self.rng``, are ignored).
             If a source subclass, it is instantiated with the
             given parameters, and returned. If None, a
             new instance of ``cpoisson_source`` is returned, with
@@ -1372,7 +1404,10 @@ class SDE:
         -------
         An object complying with the ``source`` protocol,
         instantiating the requested stochasticity source.
-        The shape of source values is set to ``wshape``.
+        The shape and paths of source values are set to
+        ``self.paths``, ``self.wshape`` respectively.
+        If instantiated, ``dtype=self.dtype`` and ``rng=self.rng``
+        keywords are used upon instantiation.
 
         See Also
         --------
@@ -1381,6 +1416,7 @@ class SDE:
         return _source_setup(dj, cpoisson_source,
                              paths=self.paths, vshape=self.wshape,
                              dtype=self.dtype,  # dtype of the source
+                             rng=self._rng_asis,
                              dn=dn,
                              ptype=ptype,  # dtype of underlying poisson source
                              lam=lam,
@@ -2036,8 +2072,8 @@ class wiener_SDE(SDE):
 
 class wiener_process(wiener_SDE, integrator):
     """
-    wiener_process(paths=1, vshape=(), dtype=None, steps=None, i0=0,
-    info=None, getinfo=True, method='euler',
+    wiener_process(paths=1, vshape=(), dtype=None, rng=None,
+    steps=None, i0=0, info=None, getinfo=True, method='euler',
     x0=0., mu=0., sigma=1., dw=None, corr=None, rho=None)
 
     Wiener process (Brownian motion) with drift.
@@ -2053,7 +2089,7 @@ class wiener_process(wiener_SDE, integrator):
 
     Parameters
     ----------
-    paths, vshape, dtype, steps, i0, info, getinfo, method
+    paths, vshape, dtype, rng, steps, i0, info, getinfo, method
         See ``SDE`` class documentation.
     x0 : array-like
         Initial condition.
@@ -2096,8 +2132,8 @@ class lognorm_SDE(SDE):
 
 class lognorm_process(lognorm_SDE, integrator):
     """
-    lognorm_process(paths=1, vshape=(), dtype=None, steps=None, i0=0,
-    info=None, getinfo=True, method='euler',
+    lognorm_process(paths=1, vshape=(), dtype=None, rng=None,
+    steps=None, i0=0, info=None, getinfo=True, method='euler',
     x0=1., mu=0., sigma=1., dw=None, corr=None, rho=None)
 
     Lognormal process.
@@ -2113,7 +2149,7 @@ class lognorm_process(lognorm_SDE, integrator):
 
     Parameters
     ----------
-    paths, vshape, dtype, steps, i0, info, getinfo, method
+    paths, vshape, dtype, rng, steps, i0, info, getinfo, method
         See ``SDE`` class documentation.
     x0 : array-like
         Initial condition.
@@ -2165,8 +2201,8 @@ class ornstein_uhlenbeck_SDE(SDE):
 
 class ornstein_uhlenbeck_process(ornstein_uhlenbeck_SDE, integrator):
     """
-    ornstein_uhlenbeck_process(paths=1, vshape=(), dtype=None, steps=None,
-    i0=0, info=None, getinfo=True, method='euler',
+    ornstein_uhlenbeck_process(paths=1, vshape=(), dtype=None, rng=None,
+    steps=None, i0=0, info=None, getinfo=True, method='euler',
     x0=0., theta=0., k=1., sigma=1., dw=None, corr=None, rho=None)
 
     Ornstein-Uhlenbeck process (mean-reverting Brownian motion).
@@ -2182,7 +2218,7 @@ class ornstein_uhlenbeck_process(ornstein_uhlenbeck_SDE, integrator):
 
     Parameters
     ----------
-    paths, vshape, dtype, steps, i0, info, getinfo, method
+    paths, vshape, dtype, rng, steps, i0, info, getinfo, method
         See ``SDE`` class documentation.
     x0 : array-like
         Initial condition.
@@ -2238,8 +2274,8 @@ class hull_white_SDE(SDE):
 
 class hull_white_process(hull_white_SDE, integrator):
     """
-    hull_white_process(paths=1, vshape=(), dtype=None, steps=None, i0=0,
-    info=None, getinfo=True, method='euler',
+    hull_white_process(paths=1, vshape=(), dtype=None, rng=None,
+    steps=None, i0=0, info=None, getinfo=True, method='euler',
     factors=1, x0=0., theta=0., k=1., sigma=1., dw=None, corr=None, rho=None)
 
     F-factors Hull-White process (sum of F correlated mean-reverting Brownian
@@ -2258,7 +2294,7 @@ class hull_white_process(hull_white_SDE, integrator):
 
     Parameters
     ----------
-    paths, vshape, dtype, steps, i0, info, getinfo, method
+    paths, vshape, dtype, rng, steps, i0, info, getinfo, method
         See ``SDE`` class documentation.
     x0 : array-like
         Initial condition.
@@ -2281,8 +2317,8 @@ class hull_white_process(hull_white_SDE, integrator):
 
 class hull_white_1factor_process(ornstein_uhlenbeck_process):
     """
-hull_white_1factor_process(paths=1, vshape=(), dtype=None, steps=None, i0=0,
-info=None, getinfo=True, method='euler',
+hull_white_1factor_process(paths=1, vshape=(), dtype=None, rng=None,
+steps=None, i0=0, info=None, getinfo=True, method='euler',
 x0=0., theta=0., k=1., sigma=1., dw=None, corr=None, rho=None)
 
 1-factor Hull-White process (F=1 Hull-White process with F-index
@@ -2320,8 +2356,8 @@ class cox_ingersoll_ross_SDE(SDE):
 
 class cox_ingersoll_ross_process(cox_ingersoll_ross_SDE, integrator):
     """
-    cox_ingersoll_ross_process(paths=1, vshape=(), dtype=None, steps=None,
-    i0=0, info=None, getinfo=True, method='euler',
+    cox_ingersoll_ross_process(paths=1, vshape=(), dtype=None, rng=None,
+    steps=None, i0=0, info=None, getinfo=True, method='euler',
     x0=1., theta=1., k=1., xi=1., dw=None, corr=None, rho=None)
 
     Cox-Ingersoll-Ross mean reverting process.
@@ -2337,7 +2373,7 @@ class cox_ingersoll_ross_process(cox_ingersoll_ross_SDE, integrator):
 
     Parameters
     ----------
-    paths, vshape, dtype, steps, i0, info, getinfo, method
+    paths, vshape, dtype, rng, steps, i0, info, getinfo, method
         See ``SDE`` class documentation.
     x0 : array-like
         Initial condition.
@@ -2409,8 +2445,8 @@ class full_heston_SDE(SDEs):
 
 
 class full_heston_process(full_heston_SDE, integrator):
-    """full_heston_process(paths=1, vshape=(), dtype=None, steps=None,
-    i0=0, info=None, getinfo=True, method='euler',
+    """full_heston_process(paths=1, vshape=(), dtype=None, rng=None,
+    steps=None, i0=0, info=None, getinfo=True, method='euler',
     x0=1., mu=0., sigma=1., y0=1., theta=1., k=1., xi=1.,
     dw=None, corr=None, rho=None)
 
@@ -2435,7 +2471,7 @@ class full_heston_process(full_heston_SDE, integrator):
 
     Parameters
     ----------
-    paths, vshape, dtype, steps, i0, info, getinfo, method
+    paths, vshape, dtype, rng, steps, i0, info, getinfo, method
         See ``SDE`` class documentation.
     x0, y0 : array-like
         Initial conditions for ``x(t)`` and ``y(t)`` processes respectively.
@@ -2506,8 +2542,8 @@ class heston_SDE(full_heston_SDE):
 
 
 class heston_process(heston_SDE, integrator):
-    """heston_process(paths=1, vshape=(), dtype=None, steps=None,
-    i0=0, info=None, getinfo=True, method='euler',
+    """heston_process(paths=1, vshape=(), dtype=None, rng=None,
+    steps=None, i0=0, info=None, getinfo=True, method='euler',
     x0=1., mu=0., sigma=1., y0=1., theta=1., k=1., xi=1.,
     dw=None, corr=None, rho=None)
 
@@ -2518,7 +2554,7 @@ class heston_process(heston_SDE, integrator):
 
     Parameters
     ----------
-    paths, vshape, dtype, steps, i0, info, getinfo, method
+    paths, vshape, dtype, rng, steps, i0, info, getinfo, method
         See ``SDE`` class documentation.
     x0, mu, sigma, y0, theta, k, xi, dw, corr, rho
         See ``full_heston_process`` class documentation.
@@ -2588,8 +2624,8 @@ class jumpdiff_SDE(SDE):
 
 
 class jumpdiff_process(jumpdiff_SDE, integrator):
-    """jumpdiff_process(paths=1, vshape=(), dtype=None, steps=None,
-    i0=0, info=None, getinfo=True, method='euler',
+    """jumpdiff_process(paths=1, vshape=(), dtype=None, rng=None,
+    steps=None, i0=0, info=None, getinfo=True, method='euler',
     x0=1., mu=0., sigma=1., dw=None, corr=None, rho=None, dj=None, dn=None,
     ptype=int, lam=1., y=None)
 
@@ -2609,7 +2645,7 @@ class jumpdiff_process(jumpdiff_SDE, integrator):
 
     Parameters
     ----------
-    paths, vshape, dtype, steps, i0, info, getinfo, method
+    paths, vshape, dtype, rng, steps, i0, info, getinfo, method
         See ``SDE`` class documentation.
     x0 : array-like
         Initial condition.
@@ -2675,8 +2711,8 @@ class merton_jumpdiff_SDE(jumpdiff_SDE):
 
 
 class merton_jumpdiff_process(merton_jumpdiff_SDE, integrator):
-    """merton_jumpdiff_process(paths=1, vshape=(), dtype=None, steps=None,
-    i0=0, info=None, getinfo=True, method='euler',
+    """merton_jumpdiff_process(paths=1, vshape=(), dtype=None, rng=None,
+    steps=None, i0=0, info=None, getinfo=True, method='euler',
     x0=1., mu=0., sigma=1., dw=None, corr=None, rho=None, dj=None, dn=None,
     ptype=int, lam=1., a=0., b=1.)
 
@@ -2712,8 +2748,8 @@ class kou_jumpdiff_SDE(jumpdiff_SDE):
 
 
 class kou_jumpdiff_process(kou_jumpdiff_SDE, integrator):
-    """kou_jumpdiff_process(paths=1, vshape=(), dtype=None, steps=None,
-    i0=0, info=None, getinfo=True, method='euler',
+    """kou_jumpdiff_process(paths=1, vshape=(), dtype=None, rng=None,
+    steps=None, i0=0, info=None, getinfo=True, method='euler',
     x0=1., mu=0., sigma=1., dw=None, corr=None, rho=None, dj=None, dn=None,
     ptype=int, lam=1., a=0.5, b=0.5, pa=0.5)
 


Test patch:
diff --git a/sdepy/tests/shared.py b/sdepy/tests/shared.py
index aee744d..7328213 100644
--- a/sdepy/tests/shared.py
+++ b/sdepy/tests/shared.py
@@ -84,6 +84,7 @@ class _pytest_tester:
 
 # in all test modules, 'sp' is the package
 # to be tested
+import sdepy
 import sdepy as sp
 from sdepy import _config
 
@@ -109,7 +110,8 @@ KFUNC = _config.KFUNC
 # (see eps function below)
 EPS_FACTOR = 16
 
-# each testing routine should call np.random.seed(SEED)
+# each testing routine should call legacy_seed(SEED) and
+# access random numbers via rng().random, rng().normal etc.
 SEED = 1234
 
 # DIR is the directory used, in quantitative tests,
@@ -155,6 +157,24 @@ else:
     quant = slow = lambda f: f
 
 
+def legacy_seed(seed):
+    """Use instead of np.random.seed(seed), to use the legacy
+    numpy random numbers, with the given seed, as sdepy default.
+    """
+    if sdepy.infrastructure.default_rng == np.random:
+        np.random.seed(seed)
+    else:
+        sdepy.infrastructure.default_rng = np.random.RandomState(seed)
+    # uncomment this to test using PCG64:
+    # sdepy.infrastructure.default_rng = np.random.default_rng(seed)
+
+
+def rng():
+    """Use `rng()` instead of `np.random` to access the current
+    sdepy default random number generator."""
+    return sdepy.infrastructure.default_rng
+
+
 def do(testing_function, *case_iterators, **args):
     """call a test function f across all test cases exposed
     in case_iterators, passing keyword arguments in args"""
@@ -175,10 +195,6 @@ def eps(dtype):
         return 0
 
 
-def choice(list_):
-    return list_[np.random.randint(len(list_))]
-
-
 class const_errors:
 
     def __init__(self, error):
@@ -247,7 +263,7 @@ def plot_histogram(hist, **kwargs):
 message = """
 A quantitative test has failed, running in {} definition with {} paths.
 This test relies on the reproducibility of expected errors, once
-random numbers have been seeded with np.random.seed({}), and its failure
+random numbers have been seeded with np.random.RandomState({}), and its failure
 may not necessarily indicate that the package is broken.
 Consider changing the test configuration in
 {} by setting
diff --git a/sdepy/tests/test_analytical.py b/sdepy/tests/test_analytical.py
index d0ce299..a0bbcfd 100644
--- a/sdepy/tests/test_analytical.py
+++ b/sdepy/tests/test_analytical.py
@@ -135,7 +135,7 @@ def test_warnings():
 
 
 def test_analytical():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     # all functions excluding hw2f stats
     # ----------------------------------
@@ -217,7 +217,7 @@ def analytical(f, tshape, xushape, params_shape, nvar, hw2f):
 
     # generate variables and parameters values
     def make(shape):
-        return 0.1*np.random.random(shape)
+        return 0.1*rng().random(shape)
     t = 1 + make(tshape)
     xu = 1 + make(xushape)
     if hw2f:
diff --git a/sdepy/tests/test_integrator.py b/sdepy/tests/test_integrator.py
index eae6769..d6237a5 100644
--- a/sdepy/tests/test_integrator.py
+++ b/sdepy/tests/test_integrator.py
@@ -61,7 +61,7 @@ class generator_cls(paths_generator):
 
 # tests with output timeline == integration timeline
 def test_generator_t():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     integr = generator_cls(paths=1, vshape=(), xw0=1)
     t = np.linspace(0, 10, 11)
@@ -113,7 +113,7 @@ def test_generator_t():
 
 # tests with output timeline != integration timeline
 def test_generator_steps():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     integr = generator_cls(paths=1, vshape=(), xw0=1,
                            info={})
@@ -168,7 +168,7 @@ class deep_cls(paths_generator):
 
 def test_generator_depth():
 
-    np.random.seed(SEED)
+    legacy_seed(SEED)
     integr = deep_cls(dtype=int)
     x = integr((1, 2, 3))
     assert_((x == np.array((10, 20, 30)).reshape(-1, 1)).all())
@@ -183,7 +183,7 @@ def test_generator_depth():
 # ------------------------------
 
 def test_SDE():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     paths = 3
     vshapes = ((), 2, (2, 5))
@@ -370,7 +370,7 @@ def test_SDE():
     # -------------------------------------
     def f(t, x=0, y=0, z=0, w=0):
         return ({'dw': x}, {'dw': y}, {'dw': z}, {'dw': w})
-    rm4 = np.random.random((4, 4))
+    rm4 = rng().random((4, 4))
     corr4 = np.eye(4) + 0.1 * (rm4 + rm4.T)
     xs = integrate(f)(x0=(1,)*4, corr=corr4, paths=11, steps=30)(t)
     for u in xs:
diff --git a/sdepy/tests/test_kfunc.py b/sdepy/tests/test_kfunc.py
index b252350..21cddf5 100644
--- a/sdepy/tests/test_kfunc.py
+++ b/sdepy/tests/test_kfunc.py
@@ -273,32 +273,32 @@ def test_kfunc_params():
 
     for X in (sp.wiener_source, sp.dw):
         assert_(kf(X)(paths=2).params ==
-                {'paths': 2, 'vshape': (), 'dtype': None,
+                {'paths': 2, 'vshape': (), 'dtype': None, 'rng': None,
                  'corr': None, 'rho': None})
 
     for X in (sp.odd_wiener_source, sp.odd_dw):
-        # defaults beyond paths and vshape are missed
+        # defaults beyond paths, vshape, dtype, rng are missed
         # due to _antithetics wrapping
         assert_(kf(X)(paths=2).params ==
-                {'paths': 2, 'vshape': ()})
+                {'paths': 2, 'vshape': (), 'dtype': None, 'rng': None})
 
     for X in (sp.true_wiener_source, sp.true_dw):
         assert_(kf(X)(paths=2).params ==
-                {'paths': 2, 'vshape': (), 'dtype': None,
+                {'paths': 2, 'vshape': (), 'dtype': None, 'rng': None,
                  'rtol': 'max', 't0': 0., 'z0': 0.,
                 'corr': None, 'rho': None})
 
     for X in (sp.poisson_source, sp.dn):
         assert_(kf(X)(lam=2, dtype=np.float32).params ==
                 {'paths': 1, 'vshape': (),
-                 'dtype': np.float32, 'lam': 2})
+                 'dtype': np.float32, 'rng': None, 'lam': 2})
 
     # processes
     # ---------
 
     for X in (sp.wiener_process, sp.wiener):
         assert_(kf(X)(x0=2).params ==
-                {'paths': 1, 'vshape': (), 'dtype': None,
+                {'paths': 1, 'vshape': (), 'dtype': None, 'rng': None,
                  'steps': None, 'i0': 0, 'info': None,
                  'getinfo': True, 'method': 'euler',
                  'dw': None, 'corr': None, 'rho': None,
@@ -306,7 +306,7 @@ def test_kfunc_params():
 
     for X in (sp.lognorm_process, sp.lognorm):
         assert_(kf(X)(x0=2).params ==
-                {'paths': 1, 'vshape': (), 'dtype': None,
+                {'paths': 1, 'vshape': (), 'dtype': None, 'rng': None,
                  'steps': None, 'i0': 0, 'info': None,
                  'getinfo': True, 'method': 'euler',
                  'dw': None, 'corr': None, 'rho': None,
@@ -314,7 +314,7 @@ def test_kfunc_params():
 
     for X in (sp.ornstein_uhlenbeck_process, sp.oruh):
         assert_(kf(X)(theta=2, steps=20).params ==
-                {'paths': 1, 'vshape': (), 'dtype': None,
+                {'paths': 1, 'vshape': (), 'dtype': None, 'rng': None,
                  'steps': 20, 'i0': 0, 'info': None,
                  'getinfo': True, 'method': 'euler', 'theta': 2,
                  'dw': None, 'corr': None, 'rho': None,
@@ -322,7 +322,7 @@ def test_kfunc_params():
 
     for X in (sp.hull_white_process, sp.hwff):
         assert_(kf(X)(vshape=3, factors=2).params ==
-                {'paths': 1, 'vshape': 3, 'dtype': None,
+                {'paths': 1, 'vshape': 3, 'dtype': None, 'rng': None,
                  'steps': None, 'i0': 0, 'info': None,
                  'getinfo': True, 'method': 'euler', 'theta': 0,
                  'dw': None, 'corr': None, 'rho': None,
@@ -331,7 +331,7 @@ def test_kfunc_params():
 
     for X in (sp.hull_white_1factor_process, sp.hw1f):
         assert_(kf(X)(vshape=3).params ==
-                {'paths': 1, 'vshape': 3, 'dtype': None,
+                {'paths': 1, 'vshape': 3, 'dtype': None, 'rng': None,
                  'steps': None, 'i0': 0, 'info': None,
                  'getinfo': True, 'method': 'euler', 'theta': 0.,
                  'dw': None, 'corr': None, 'rho': None,
@@ -339,7 +339,7 @@ def test_kfunc_params():
 
     for X in (sp.cox_ingersoll_ross_process, sp.cir):
         assert_(kf(X)(vshape=3).params ==
-                {'paths': 1, 'vshape': 3, 'dtype': None,
+                {'paths': 1, 'vshape': 3, 'dtype': None, 'rng': None,
                  'steps': None, 'i0': 0, 'info': None,
                  'getinfo': True, 'method': 'euler',
                  'dw': None, 'corr': None, 'rho': None,
@@ -348,7 +348,7 @@ def test_kfunc_params():
     for X in (sp.heston_process, sp.heston,
               sp.full_heston_process, sp.heston_xy):
         assert_(kf(X)(vshape=2).params ==
-                {'paths': 1, 'vshape': 2, 'dtype': None,
+                {'paths': 1, 'vshape': 2, 'dtype': None, 'rng': None,
                  'steps': None, 'i0': 0, 'info': None,
                  'getinfo': True, 'method': 'euler',
                  'dw': None, 'corr': None, 'rho': None,
@@ -357,7 +357,7 @@ def test_kfunc_params():
 
     for X in (sp.jumpdiff_process, sp.jumpdiff):
         assert_(kf(X)(vshape=2).params ==
-                {'paths': 1, 'vshape': 2, 'dtype': None,
+                {'paths': 1, 'vshape': 2, 'dtype': None, 'rng': None,
                  'steps': None, 'i0': 0, 'info': None,
                  'getinfo': True, 'method': 'euler', 'ptype': int,
                  'dw': None, 'corr': None, 'rho': None,
@@ -366,7 +366,7 @@ def test_kfunc_params():
 
     for X in (sp.merton_jumpdiff_process, sp.mjd):
         assert_(kf(X)(vshape=2).params ==
-                {'paths': 1, 'vshape': 2, 'dtype': None,
+                {'paths': 1, 'vshape': 2, 'dtype': None, 'rng': None,
                  'steps': None, 'i0': 0, 'info': None,
                  'getinfo': True, 'method': 'euler', 'ptype': int,
                  'dw': None, 'corr': None, 'rho': None,
@@ -375,7 +375,7 @@ def test_kfunc_params():
 
     for X in (sp.kou_jumpdiff_process, sp.kou):
         assert_(kf(X)(vshape=2).params ==
-                {'paths': 1, 'vshape': 2, 'dtype': None,
+                {'paths': 1, 'vshape': 2, 'dtype': None, 'rng': None,
                  'steps': None, 'i0': 0, 'info': None,
                  'getinfo': True, 'method': 'euler', 'ptype': int,
                  'dw': None, 'corr': None, 'rho': None,
diff --git a/sdepy/tests/test_montecarlo.py b/sdepy/tests/test_montecarlo.py
index 95473e5..aaceaaa 100644
--- a/sdepy/tests/test_montecarlo.py
+++ b/sdepy/tests/test_montecarlo.py
@@ -33,7 +33,7 @@ test_aaa_config.config = True
 
 # main test
 def test_montecarlo_workflow():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     # do cases
     shape = [(), (2,), (3, 2)]
@@ -54,8 +54,8 @@ def montecarlo_workflow(shape, paths, dtype):
     a = montecarlo(bins='auto')
     assert_raises(ValueError, a.histogram)
     for i in range(10):
-        sample = np.random.normal(size=shape + (PATHS,)).astype(dtype)
-        sample *= (1+np.arange(size)).reshape(shape + (1,))
+        sample = (100*rng().normal(size=shape + (PATHS,))).astype(dtype)
+        sample *= (1 + np.arange(size)).reshape(shape + (1,))
         a.update(sample)
     assert_(a.paths == PATHS*10)
 
@@ -63,7 +63,7 @@ def montecarlo_workflow(shape, paths, dtype):
     str(a)
     x = np.linspace(-3*size, 3*size, 100)
     for k in range(5):
-        i = choice(indexes)
+        i = (indexes)[-1]
         a[i].pdf(x)  # default
         a[i].pdf(x, bandwidth=.7)  # gaussian kde with explicit bandwidth
         a[i].pdf(x, kind='linear', method='interp')
@@ -115,7 +115,7 @@ def montecarlo_workflow(shape, paths, dtype):
     # use paths along different axes
     e = None
     for k in range(len(shape)):
-        axis_sample = np.random.normal(
+        axis_sample = rng().normal(
             size=shape[:k] + (PATHS,) + shape[k:]).astype(dtype)
         # first run
         e = montecarlo(axis_sample, axis=k, bins=20, range=(-3, 3))
@@ -133,7 +133,7 @@ def montecarlo_workflow(shape, paths, dtype):
                     rtol=eps(f.mean().dtype))
 
     # antithetic sampling
-    sample = np.random.normal(size=(2*PATHS,) + shape)
+    sample = rng().normal(size=(2*PATHS,) + shape)
     even = montecarlo(sample, axis=0, use='even')
     odd = montecarlo(sample, axis=0, use='odd')
     for a in (even, odd):
@@ -141,7 +141,7 @@ def montecarlo_workflow(shape, paths, dtype):
         assert_(a.paths == PATHS)
     assert_allclose(even.mean(), montecarlo(sample, 0).mean())
     assert_allclose(odd.mean(), (sample[:PATHS]-sample[PATHS:]).mean(axis=0)/2)
-    sample2 = np.random.normal(size=(4*PATHS,) + shape)
+    sample2 = rng().normal(size=(4*PATHS,) + shape)
     even.update(sample2, axis=0)
     odd.update(sample2, axis=0)
     for a in (even, odd):
@@ -168,7 +168,7 @@ def montecarlo_workflow(shape, paths, dtype):
 
 # main test
 def test_montecarlo_cumulate():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     shape = [(), (2,), (3, 2)]
     paths = [1, 10]
@@ -180,7 +180,7 @@ def test_montecarlo_cumulate():
 # cases
 def montecarlo_cumulate(shape, paths, dtype):
     PATHS = paths
-    sample = np.random.normal(size=shape + (4*PATHS,)).astype(dtype)
+    sample = rng().normal(size=shape + (4*PATHS,)).astype(dtype)
 
     # full sample as first run
     a = montecarlo(sample)
@@ -220,7 +220,7 @@ def montecarlo_cumulate(shape, paths, dtype):
 @slow
 @quant
 def test_montecarlo_histogram():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     context = 'histogram'
     err_expected = load_errors(context)
@@ -277,7 +277,7 @@ def montecarlo_histogram(context, test_id, err_expected, err_realized,
     dist, support, bins, _ = dist_info[test_id]
     test_key = context + '_' + test_id
 
-    sample = dist.rvs(size=shape + (paths,))
+    sample = dist.rvs(size=shape + (paths,), random_state=rng())
     sample_bins = np.empty(shape + bins.shape)
     for i in np.ndindex(shape):
         sample_bins[i] = bins
@@ -335,7 +335,7 @@ def montecarlo_histogram(context, test_id, err_expected, err_realized,
 @slow
 @quant
 def test_montecarlo_dist():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     context = 'distr'
     err_expected = load_errors(context)
@@ -394,7 +394,7 @@ def montecarlo_dist(context, test_id, err_expected, err_realized,
     test_key = context + '_' + test_id
 
     # prepare sample
-    sample = dist.rvs(size=shape + (paths,))
+    sample = dist.rvs(size=shape + (paths,), random_state=rng())
     sample_bins = np.empty(shape + bins.shape)
     for i in np.ndindex(shape):
         sample_bins[i] = bins
diff --git a/sdepy/tests/test_process.py b/sdepy/tests/test_process.py
index e40153f..eefd155 100644
--- a/sdepy/tests/test_process.py
+++ b/sdepy/tests/test_process.py
@@ -20,7 +20,7 @@ def test_import():
 
 # enumerate test cases and launch tests
 def test_constructor():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     # do cases
     t = [2., (5.,), (0, 1.), np.linspace(0., 4., 12)]
@@ -45,7 +45,7 @@ def process_constructor(t, paths, vshape, dtype):
         tt = tt.reshape(1)
 
     # x constructor
-    x = np.random.random(tt.shape + vshape + (paths,))
+    x = rng().random(tt.shape + vshape + (paths,))
     p = process(t, x=x, dtype=dtype)
     assert_array_equal(p.t, tt)
     assert_(isinstance(p.t, np.ndarray))
@@ -54,7 +54,7 @@ def process_constructor(t, paths, vshape, dtype):
     assert_(p.dtype == np.dtype(dtype))
 
     # v constructor
-    v = np.random.random(tt.shape + vshape)
+    v = rng().random(tt.shape + vshape)
     p = process(t, v=v, dtype=dtype)
     assert_array_equal(p.t, tt)
     assert_(isinstance(p.t, np.ndarray))
@@ -63,7 +63,7 @@ def process_constructor(t, paths, vshape, dtype):
     assert_(p.dtype == np.dtype(dtype))
 
     # c constructor
-    c = np.random.random(vshape)
+    c = rng().random(vshape)
     p = process(t, c=c, dtype=dtype)
     assert_array_equal(p.t, tt)
     assert_(isinstance(p.t, np.ndarray))
@@ -98,7 +98,7 @@ def process_constructor(t, paths, vshape, dtype):
 
 # enumerate test cases and launch tests
 def test_broadcasting():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     # do cases
     t = [np.linspace(0., 4., 12)]
@@ -110,9 +110,9 @@ def test_broadcasting():
 
 # case testing
 def process_broadcasting(t, paths, vshape):
-    p = process(t, x=np.random.random(t.shape + vshape + (paths,)))
+    p = process(t, x=rng().random(t.shape + vshape + (paths,)))
     q = process(t.copy(),
-                x=np.random.random(t.shape + vshape + (paths,)))
+                x=rng().random(t.shape + vshape + (paths,)))
 
     # same processes
     a1 = p + q
@@ -195,7 +195,7 @@ def process_broadcasting(t, paths, vshape):
 
 # enumerate test cases and launch tests
 def test_interp():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     # do cases on t, s, paths and vshape
     t = [2+np.zeros(1), np.linspace(1., 4., 12)]
@@ -261,7 +261,7 @@ def process_interp(t, s, paths, vshape, dtype, kind):
 
 # simple stand alone test on interpolated values
 def test_interp_values():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     # check interpolated values
     p = process(t=(1, 2, 3.), v=(1, 4, 9.))
@@ -277,7 +277,7 @@ def test_interp_values():
 
 # enumerate test cases and launch tests
 def test_t_getitem():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     # do cases
     t = [np.linspace(1., 4., 12)]
@@ -288,7 +288,7 @@ def test_t_getitem():
 
 # case testing
 def process_t_getitem(t, paths, vshape):
-    p = process(t, x=np.random.random(t.shape + vshape + (paths,)))
+    p = process(t, x=rng().random(t.shape + vshape + (paths,)))
 
     def verify(q, newt, newx):
         assert_(isinstance(q, process))
@@ -321,7 +321,7 @@ def process_t_getitem(t, paths, vshape):
 
 # enumerate test cases and launch tests
 def test_p_getitem():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     # do cases
     t = [np.zeros(1) + 2, np.linspace(1., 4., 12)]
@@ -332,7 +332,7 @@ def test_p_getitem():
 
 # case testing
 def process_p_getitem(t, paths, vshape):
-    p = process(t, x=np.random.random(t.shape + vshape + (paths,)))
+    p = process(t, x=rng().random(t.shape + vshape + (paths,)))
 
     def verify(q, newt, newx):
         assert_(isinstance(q, process))
@@ -368,7 +368,7 @@ def process_p_getitem(t, paths, vshape):
 
 # enumerate test cases and launch tests
 def test_v_getitem():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     # do cases
     t = [np.zeros(1) + 2, np.linspace(1., 4., 11)]
@@ -389,7 +389,7 @@ def process_v_getitem(t, paths):
             verify(p[('v',) + i], p.t, p[(slice(None),) + i + (slice(None),)])
 
     # vshape = ()
-    p = process(t, x=np.random.random(t.shape + () + (paths,)))
+    p = process(t, x=rng().random(t.shape + () + (paths,)))
     q = p['v', ()]
     verify(q, p.t, p)
     iall(p)
@@ -399,7 +399,7 @@ def process_v_getitem(t, paths):
     assert_raises(IndexError, ierr)
 
     # vshape = (5,)
-    p = process(t, x=np.random.random(t.shape + (5,) + (paths,)))
+    p = process(t, x=rng().random(t.shape + (5,) + (paths,)))
     iall(p)
     q = p['v', 0:2]
     verify(q, p.t, p[:, 0:2, :])
@@ -420,7 +420,7 @@ def process_v_getitem(t, paths):
     verify(q, p.t, p[:, i, :])
 
     # vshape = (5, 7)
-    p = process(t, x=np.random.random(t.shape + (5, 7) + (paths,)))
+    p = process(t, x=rng().random(t.shape + (5, 7) + (paths,)))
     iall(p)
     q = p['v', 0:2]
     verify(q, p.t, p[:, 0:2, :])
@@ -463,8 +463,8 @@ def process_v_getitem(t, paths):
 
 # enumerate test cases and launch tests
 def test_properties():
-    np.random.seed(SEED)
-    p = process(t=(1., 2, 4), x=np.random.random((3, 5, 7, 11)))
+    legacy_seed(SEED)
+    p = process(t=(1., 2, 4), x=rng().random((3, 5, 7, 11)))
 
     def sharemem(a, b):
         return a.__array_interface__['data'][0] == \
@@ -498,7 +498,7 @@ def process_properties(t, paths, vshape, dtype):
     if tt.ndim == 0:
         tt = tt.reshape(1)
 
-    p = process(t, x=np.random.random(tt.shape + vshape + (paths,)),
+    p = process(t, x=rng().random(tt.shape + vshape + (paths,)),
                 dtype=dtype)
     assert_(p.shape == p.t.shape + p.vshape + (p.paths,))
     assert_array_equal(p.t, p.tx.flatten())
@@ -517,10 +517,10 @@ def process_properties(t, paths, vshape, dtype):
 # ----------------
 
 def test_shapeas():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     def mkp(vshape):
-        return process(t=(0, 1), x=np.random.random((2,) + vshape + (3,)))
+        return process(t=(0, 1), x=rng().random((2,) + vshape + (3,)))
 
     for vshape in ((), (5,), (5, 7)):
         p = mkp(vshape)
@@ -555,9 +555,9 @@ def test_shapeas():
 # ----------------
 
 def test_copy():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
-    p = process(t=(1, 2, 3), x=np.random.random((3, 5, 7, 11)))
+    p = process(t=(1, 2, 3), x=rng().random((3, 5, 7, 11)))
     q = p.pcopy()
     qt = p.tcopy()
     qx = p.xcopy()
@@ -581,9 +581,9 @@ def test_copy():
 # -----------------------
 
 def test_summary():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
-    p = process(t=(1, 2, 3), x=1 + np.random.random((3, 5, 7, 11)),
+    p = process(t=(1, 2, 3), x=1 + rng().random((3, 5, 7, 11)),
                 dtype=np.float64)
 
     funcs = ('min', 'max', 'sum', 'mean', 'var', 'std')
@@ -686,11 +686,11 @@ def test_summary():
 # ----------------------------------------
 
 def test_increments():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     t = np.linspace(1, 2, 100)
     t *= t
-    p = process(t, x=1 + np.random.random(t.shape + (2, 3) + (5,)),
+    p = process(t, x=1 + rng().random(t.shape + (2, 3) + (5,)),
                 dtype=np.float64)
 
     dp = p.tdiff()
@@ -727,7 +727,7 @@ def test_increments():
 
 # enumerate test cases and launch tests
 def test_chf_cdf():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     t = [2., (5.,), (0, 1.), np.linspace(0., 4., 17)]
     paths = [1, 19]
@@ -744,7 +744,7 @@ def process_chf_cdf(t, s, u, paths, vshape):
         tt = tt.reshape(1)
     ss, uu = np.asarray(s), np.asarray(u)
 
-    p = process(t, x=np.random.random(tt.shape + vshape + (paths,)))
+    p = process(t, x=rng().random(tt.shape + vshape + (paths,)))
 
     f = p.chf(u)
     assert_(f.shape == tt.shape + uu.shape + vshape)
@@ -775,14 +775,14 @@ def process_chf_cdf(t, s, u, paths, vshape):
 
 # enumerate test cases and launch tests
 def test_piecewise():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     dtype = [np.float64, np.float32, np.float16]
     paths = [None, 5]
     vshape = [(3,), (2, 3)]
     mode = [None, 'mid', 'forward', 'backward']
     shift = [0, -2, -4]
-    do(tst_piecewise, [np.float], paths, vshape, mode, shift)
+    do(tst_piecewise, [float], paths, vshape, mode, shift)
     do(tst_piecewise,      dtype, paths,   [()], mode,   [0])
 
     with assert_raises(ValueError):
@@ -797,11 +797,11 @@ def tst_piecewise(dtype, paths, vshape, mode, shift):
     t = np.array((1, 2, 3), dtype=dtype) + shift
     if paths is None:
         paths = 1
-        v = vv = 1 + np.random.random((3,) + vshape).astype(dtype)
+        v = vv = 1 + rng().random((3,) + vshape).astype(dtype)
         x = v.reshape(v.shape + (1,))
         xx = None
     else:
-        x = xx = 1 + np.random.random((3,) + vshape + (paths,)).astype(dtype)
+        x = xx = 1 + rng().random((3,) + vshape + (paths,)).astype(dtype)
         vv = None
 
     if mode is None:
diff --git a/sdepy/tests/test_processes.py b/sdepy/tests/test_processes.py
index 9c5bc38..86745c3 100644
--- a/sdepy/tests/test_processes.py
+++ b/sdepy/tests/test_processes.py
@@ -79,13 +79,13 @@ all_shortcuts = (wiener, lognorm,
 # enumerate test cases with constant parameters and launch tests
 # --------------------------------------------------------------
 def test_processes():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     # setup some parameter test values
     # --------------------------------
     CM2 = ((1, .2), (.2, 1))
     CM3 = ((1, .2, -.3), (.2, 1, .1), (-.3, .1, 1))
-    CM6 = np.eye(6) + 0.1*np.random.random((6, 6))
+    CM6 = np.eye(6) + 0.1*rng().random((6, 6))
     CM6 = (CM6 + CM6.T)/2
 
     # processes to act as sources
@@ -97,11 +97,14 @@ def test_processes():
                          )(np.linspace(0, 1, 100))
 
     # a 'bare' object with source protocol
-    def S2(s, ds):
-        return np.random.normal(size=(2, 11))*sqrt(np.abs(ds))
-
-    S2.vshape = (2,)
-    S2.paths = 11
+    def make_S2(rng_):
+        if rng_ is None:
+            rng_ = rng()
+        def S2(s, ds):
+            return rng_.normal(size=(2, 11))*sqrt(np.abs(ds))
+        S2.vshape = (2,)
+        S2.paths = 11
+        return S2
 
     # a generic scipy random variable
     rv = scipy.stats.trapz(c=.1, d=.2)
@@ -119,9 +122,10 @@ def test_processes():
     do(processes, cls, params, paths, dtype)
 
     paths, dtype = [11], [None]
-    params = (dict(vshape=2, dw=Z2),
-              dict(vshape=2, dw=S2),
-              )
+    params = (
+        dict(vshape=2, dw=Z2),
+        lambda rng=None: dict(vshape=2, dw=make_S2(rng)),
+    )
     cls = list(cls)
     for p in (heston_process, heston,
               full_heston_process, heston_xy,
@@ -142,8 +146,10 @@ def test_processes():
         dict(vshape=(2,), corr=CM2),
         dict(vshape=(2, 3)),
         dict(vshape=(2, 3), corr=CM3),
-        dict(vshape=(2, 3), dw=wiener_source(paths=11, vshape=(2, 3))),
-        dict(vshape=(), dw=true_wiener_source(paths=11, vshape=())),
+        lambda rng=None: dict(
+            vshape=(2, 3), dw=wiener_source(paths=11, vshape=(2, 3), rng=rng)),
+        lambda rng=None: dict(
+            vshape=(), dw=true_wiener_source(paths=11, vshape=(), rng=rng)),
         dict(vshape=(2, 3), dw=Z23),
         )
     do(processes, cls, params, paths, dtype)
@@ -153,9 +159,11 @@ def test_processes():
         dict(vshape=2, x0=.1, theta=.2, k=.4, sigma=.5),
         dict(vshape=(2,), corr=CM2),
         dict(vshape=(2,), dw=Z2),
-        dict(vshape=(2, 3), dw=wiener_source(paths=11, vshape=(2, 3))),
+        lambda rng=None: dict(
+            vshape=(2, 3), dw=wiener_source(paths=11, vshape=(2, 3), rng=rng)),
         dict(vshape=(2, 3), dw=Z23),
-        dict(vshape=(), dw=true_wiener_source(paths=11, vshape=()))
+        lambda rng=None: dict(
+            vshape=(), dw=true_wiener_source(paths=11, vshape=(), rng=rng))
     )
     do(processes, cls, params, paths, dtype)
 
@@ -166,25 +174,28 @@ def test_processes():
         dict(vshape=(), factors=2,
              x0=((.1,), (.2,)), theta=((.4,), (.5,)), k=((.6,), (.7,)),
              sigma=((.8,), (.9,)), rho=0.25),
-        dict(vshape=(), factors=2,
-             dw=wiener_source(paths=11, vshape=(2,),
-                              corr=((1, .25,), (.25, 1)))),
+        lambda rng=None: dict(
+            vshape=(), factors=2,
+            dw=wiener_source(paths=11, vshape=(2,), rng=rng,
+                             corr=((1, .25,), (.25, 1)))),
         dict(vshape=(3,), factors=2, dw=Z32),
-        dict(vshape=(), factors=3,
-             dw=wiener_source(paths=11, vshape=(3,), corr=CM3)),
+        lambda rng=None: dict(vshape=(), factors=3,
+             dw=wiener_source(paths=11, vshape=(3,), corr=CM3, rng=rng)),
         dict(vshape=(2, 3), factors=5,
              sigma=np.arange(2*3*5).reshape(2, 3, 5, 1)/(300)),
         dict(vshape=2, factors=3, dw=Z23),
         dict(vshape=(), factors=2, dw=Z2),
-        dict(vshape=(), factors=2, dw=S2),
+        lambda rng=None: dict(vshape=(), factors=2, dw=make_S2(rng)),
     )
     do(processes, cls, params, paths, dtype)
 
     cls = (hull_white_1factor_process, hw1f)
     params = (
         dict(vshape=2, x0=.1, theta=.2, k=.4, sigma=.5),
-        dict(vshape=(2, 3), dw=wiener_source(paths=11, vshape=(2, 3))),
-        dict(vshape=(), dw=true_wiener_source(paths=11, vshape=())),
+        lambda rng=None: dict(
+            vshape=(2, 3), dw=wiener_source(paths=11, vshape=(2, 3), rng=rng)),
+        lambda rng=None: dict(
+            vshape=(), dw=true_wiener_source(paths=11, vshape=(), rng=rng)),
     )
     do(processes, cls, params, paths, dtype)
 
@@ -193,8 +204,10 @@ def test_processes():
         dict(vshape=2, x0=.1, theta=.2, k=.3, xi=.4),
         dict(vshape=(2,), corr=CM2),
         dict(vshape=(2,), rho=-.1),
-        dict(vshape=(2, 3), dw=wiener_source(paths=11, vshape=(2, 3))),
-        dict(vshape=(), dw=true_wiener_source(paths=11, vshape=())),
+        lambda rng=None: dict(
+            vshape=(2, 3), dw=wiener_source(paths=11, vshape=(2, 3), rng=rng)),
+        lambda rng=None: dict(
+            vshape=(), dw=true_wiener_source(paths=11, vshape=(), rng=rng)),
         dict(vshape=(2, 3), dw=Z23)
     )
     do(processes, cls, params, paths, dtype)
@@ -208,10 +221,12 @@ def test_processes():
              y0=.4, theta=.5, k=.6, xi=.7),
         dict(vshape=(3,), corr=CM6),
         dict(vshape=(3,), rho=(.1, .2, .3)),
-        dict(vshape=(5,), dw=wiener_source(paths=11, vshape=(10,))),
-        dict(vshape=10, dw=true_wiener_source(paths=11, vshape=(20,))),
+        lambda rng=None: dict(
+            vshape=(5,), dw=wiener_source(paths=11, vshape=(10,), rng=rng)),
+        lambda rng=None: dict(
+            vshape=10, dw=true_wiener_source(paths=11, vshape=(20,), rng=rng)),
         dict(vshape=(), dw=Z2),
-        dict(vshape=(), dw=S2),
+        lambda rng=None: dict(vshape=(), dw=make_S2(rng)),
     )
     do(processes, cls, params, paths, dtype)
 
@@ -229,11 +244,12 @@ def test_processes():
         dict(vshape=(), corr=CM2),
         dict(vshape=(1,), corr=CM2),
         dict(vshape=(2, 3, 1), rho=.25),
-        dict(vshape=(3), dw=wiener_source(paths=11, vshape=(6,))),
-        dict(vshape=(2, 3, 1),
-             dw=true_wiener_source(paths=11, vshape=(2, 3, 2))),
+        lambda rng=None: dict(
+            vshape=(3), dw=wiener_source(paths=11, vshape=(6,), rng=rng)),
+        lambda rng=None: dict(vshape=(2, 3, 1),
+            dw=true_wiener_source(paths=11, vshape=(2, 3, 2), rng=rng)),
         dict(vshape=(), dw=Z2),
-        dict(vshape=(), dw=S2),
+        lambda rng=None: dict(vshape=(), dw=make_S2(rng)),
     )
     do(processes, cls, params, paths, dtype)
 
@@ -245,17 +261,24 @@ def test_processes():
         dict(vshape=2, x0=.1, mu=.2, sigma=.3,
              lam=.4),
         dict(vshape=(2,), corr=CM2),
-        dict(vshape=(2, 3), dw=wiener_source(paths=11, vshape=(2, 3))),
-        dict(vshape=(3,), dw=true_wiener_source(paths=11, vshape=(3,)),
-             dj=cpoisson_source(paths=11, vshape=(3,))),
-        dict(vshape=(3,), dw=true_wiener_source(paths=11, vshape=(3,)),
-             dj=true_cpoisson_source(paths=11, vshape=(3,))),
-        dict(vshape=2, dj=poisson_source(paths=11, vshape=2)),
-        dict(vshape=2, dj=true_poisson_source(paths=11, vshape=2)),
-        dict(vshape=2, dj=cpoisson_source(paths=11, vshape=2, y=rv)),
-        dict(vshape=2, dj=true_cpoisson_source(paths=11, vshape=2, y=rv)),
+        lambda rng=None: dict(
+            vshape=(2, 3), dw=wiener_source(paths=11, vshape=(2, 3), rng=rng)),
+        lambda rng=None: dict(
+            vshape=(3,), dw=true_wiener_source(paths=11, vshape=(3,), rng=rng),
+            dj=cpoisson_source(paths=11, vshape=(3,), rng=rng)),
+        lambda rng=None: dict(
+            vshape=(3,), dw=true_wiener_source(paths=11, vshape=(3,), rng=rng),
+            dj=true_cpoisson_source(paths=11, vshape=(3,), rng=rng)),
+        lambda rng=None: dict(
+            vshape=2, dj=poisson_source(paths=11, vshape=2, rng=rng)),
+        lambda rng=None: dict(
+            vshape=2, dj=true_poisson_source(paths=11, vshape=2, rng=rng)),
+        lambda rng=None: dict(
+            vshape=2, dj=cpoisson_source(paths=11, vshape=2, y=rv, rng=rng)),
+        lambda rng=None: dict(
+            vshape=2, dj=true_cpoisson_source(paths=11, vshape=2, y=rv, rng=rng)),
         dict(vshape=2, dw=Z2, dj=Z2),
-        dict(vshape=2, dw=S2, dj=S2),
+        lambda rng=None: dict(vshape=2, dw=make_S2(rng), dj=make_S2(rng)),
     )
     do(processes, cls, params, paths, dtype)
 
@@ -287,7 +310,7 @@ def test_processes():
 # enumerate test cases with time-dependent parameters and launch tests
 # --------------------------------------------------------------------
 def test_processes_local():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     # setup some parameter test values
     # --------------------------------
@@ -304,9 +327,14 @@ def test_processes_local():
 
     def CM2(t): return ((1, .2*t/10), (.2*t/10, 1))
 
-    def CM6(t):
-        C = np.eye(6) + t*0.1*np.random.random((6, 6))
-        return (C + C.T)/2
+    rng_ = rng()
+    def make_CM6(rng_):
+        if rng_ is None:
+            rng_ = rng()
+        def CM6(t):
+            C = np.eye(6) + t*0.1*rng_.random((6, 6))
+            return (C + C.T)/2
+        return CM6
 
     def R(t): return 0.5 - 0.1*t
 
@@ -325,19 +353,21 @@ def test_processes_local():
         dict(vshape=(2,), mu=A, sigma=C, corr=CM2P),
         dict(vshape=(2,), mu=A, sigma=C, rho=R),
         dict(vshape=(2,), mu=A, sigma=C, rho=RP),
-        dict(vshape=(6,), mu=A, sigma=C, corr=CM6),
+        lambda rng=None: dict(vshape=(6,), mu=A, sigma=C, corr=make_CM6(rng)),
         dict(vshape=(6,), mu=A, sigma=C, rho=R3P),
-        dict(vshape=(6,), mu=A, sigma=C,
-             dw=true_wiener_source(paths=11, vshape=6, rho=R3P)),
-        )
+        lambda rng=None: dict(
+            vshape=(6,), mu=A, sigma=C,
+            dw=true_wiener_source(paths=11, vshape=6, rho=R3P, rng=rng)),
+    )
     do(processes, cls, params, paths, dtype)
 
     cls = (ornstein_uhlenbeck_process, )
     params = (
         dict(vshape=(), x0=.1, theta=A, k=B, sigma=C),
         dict(vshape=(2,), theta=A, k=B, sigma=C, corr=CM2),
-        dict(vshape=(2,), theta=A, k=B, sigma=C,
-             dw=true_wiener_source(paths=11, vshape=2, corr=CM2)),
+        lambda rng=None: dict(
+            vshape=(2,), theta=A, k=B, sigma=C,
+            dw=true_wiener_source(paths=11, vshape=2, corr=CM2, rng=rng)),
     )
     do(processes, cls, params, paths, dtype)
 
@@ -360,7 +390,8 @@ def test_processes_local():
     params = (
         dict(vshape=(), x0=.1, mu=A, sigma=B,
              y0=.2, theta=C, k=B, xi=A, rho=B),
-        dict(vshape=(3,), theta=A, k=B, xi=C, corr=CM6),
+        lambda rng=None: dict(
+            vshape=(3,), theta=A, k=B, xi=C, corr=make_CM6(rng)),
         dict(vshape=(3,), theta=A, k=B, xi=C, rho=R3P)
     )
     do(processes, cls, params, paths, dtype)
@@ -369,7 +400,8 @@ def test_processes_local():
     params = (
         dict(vshape=(), x0=.1, mu=A, sigma=B,
              y0=.2, theta=A, k=B, xi=C, rho=B),
-        dict(vshape=(3,), theta=A, k=B, xi=C, corr=CM6)
+        lambda rng=None: dict(
+            vshape=(3,), theta=A, k=B, xi=C, corr=make_CM6(rng))
     )
     do(processes, cls, params, paths, dtype)
 
@@ -402,6 +434,12 @@ def test_processes_local():
 # case testing
 def processes(cls, params, paths, dtype):
 
+    # store params maker to be used for rng tests
+    make_params = params
+    if callable(make_params):
+        # recover params with default rng for other tests
+        params = make_params()
+
     # print(cls.__name__, sep='', end='')
     P = cls(**params, paths=paths, dtype=dtype)
 
@@ -437,6 +475,43 @@ def processes(cls, params, paths, dtype):
                     # initial values are the same for all paths
                     assert_((p[i0] == p[i0, ..., 0][..., np.newaxis]).all())
 
+                # extensive test of rng parameter (too slow to be included
+                # in the test suite)
+                if False:
+                    try:
+                        make_rngs = (
+                            np.random.default_rng,
+                            np.random.RandomState,
+                            lambda z: np.random.Generator(
+                                np.random.PCG64(z)),
+                        )
+                    except AttributeError:
+                        continue
+                    for make_rng in make_rngs:
+                        rng1 = make_rng(SEED)
+                        rng2 = make_rng(SEED)
+                        assert rng1 is not rng2
+                        if callable(make_params):
+                            params1 = make_params(rng1)
+                            params2 = make_params(rng2)
+                        else:
+                            params1 = params2 = params
+                        P1 = cls(**params1, paths=paths,
+                                 dtype=dtype, rng=rng1,
+                                 steps=steps, i0=i0)
+                        P2 = cls(**params2, paths=paths,
+                                 dtype=dtype, rng=rng2,
+                                 steps=steps, i0=i0)
+                        assert P1.rng is rng1
+                        assert P2.rng is rng2
+                        ps1, ps2 = P1(t), P2(t)
+                        ps1, ps2 = [
+                            ps if isinstance(ps, (tuple, list)) else (ps,)
+                            for ps in (ps1, ps2)]
+                        assert len(ps1) == len(ps2)
+                        for z, w in zip(ps1, ps2):
+                            assert_allclose(z, w)
+
     # check initial conditions on the last computed p
     if cls in {wiener_process, lognorm_process,
                ornstein_uhlenbeck_process,
@@ -465,6 +540,51 @@ def processes(cls, params, paths, dtype):
         assert_allclose(p0, x0, rtol=eps(p.dtype))
         assert_allclose(q0, y0, rtol=eps(p.dtype))
 
+    # test rng parameter
+    try:
+        make_rngs = (
+            np.random.default_rng,
+            np.random.RandomState,
+            lambda z: np.random.Generator(np.random.PCG64(z)),
+        )
+    except AttributeError:
+        return
+    for make_rng in make_rngs:
+        rng1 = make_rng(SEED)
+        rng2 = make_rng(SEED)
+        assert rng1 is not rng2
+        if callable(make_params):
+            params1 = make_params(rng1)
+            params2 = make_params(rng2)
+        else:
+            params1 = params2 = params
+        t = tlist[-1]
+        P1 = cls(**params1, paths=paths, dtype=dtype, rng=rng1)
+        P2 = cls(**params2, paths=paths, dtype=dtype, rng=rng2)
+        assert P1.rng is rng1
+        assert P2.rng is rng2
+        ps1, ps2 = P1(t), P2(t)
+        ps1, ps2 = [ps if isinstance(ps, (tuple, list)) else (ps,)
+                    for ps in (ps1, ps2)]
+        assert len(ps1) == len(ps2)
+        for z, w in zip(ps1, ps2):
+            assert_allclose(z, w)
+        # test that global sdepy.infrastructure.default_rng
+        # propagates correctly as a default
+        tmp = sdepy.infrastructure.default_rng
+        sdepy.infrastructure.default_rng = make_rng(SEED)
+        if callable(make_params):
+            params3 = make_params()
+        else:
+            params3 = params
+        P3 = cls(**params3, paths=paths, dtype=dtype, rng=None)
+        assert P3.rng is sdepy.infrastructure.default_rng
+        ps3 = P3(t)
+        ps3 = ps3 if isinstance(ps3, (tuple, list)) else (ps3,)
+        assert len(ps1) == len(ps3)
+        assert_allclose(ps1, ps3)
+        sdepy.infrastructure.default_rng = tmp
+
 
 # stand alone test
 def test_processes_exceptions():
@@ -476,7 +596,7 @@ def test_processes_exceptions():
 
 
 def test_jumps():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     # jump diffusion process integrated as is
     # (jumpdiff_process integrates the logarithm)
@@ -550,7 +670,7 @@ def test_processes_misc():
     # test exactness of wiener_process and lognorm_process
     # with constant parameters
 
-    np.random.seed(SEED)
+    legacy_seed(SEED)
     paths = 31
     x0 = 10
     mu, sigma = .2, .7
diff --git a/sdepy/tests/test_quant.py b/sdepy/tests/test_quant.py
index 71ea125..975b103 100644
--- a/sdepy/tests/test_quant.py
+++ b/sdepy/tests/test_quant.py
@@ -296,7 +296,7 @@ def quant_case(case, context, err_expected, err_realized, PATHS):
     # plot paths for visual inspection (no testing)
     if PLOT:
         print('plotting...')
-        np.random.seed(SEED)
+        legacy_seed(SEED)
         fig = bfig()
         plt.title('{}: 30 sample paths'.format(Xid))
         plt.xlabel('t')
@@ -306,7 +306,7 @@ def quant_case(case, context, err_expected, err_realized, PATHS):
 
     # check mean, std and variance
     if Xmean is not None:
-        np.random.seed(SEED)
+        legacy_seed(SEED)
         x = X(t, paths=PATHS)
         if 'heston' in Xid:
             log(x, out=x)
@@ -323,7 +323,7 @@ def quant_case(case, context, err_expected, err_realized, PATHS):
         del x
 
     # check pdf and cdf
-    np.random.seed(SEED)
+    legacy_seed(SEED)
     a = sp.montecarlo(bins=200)
     for i in range(10):
         x = X(s, paths=PATHS)
@@ -533,7 +533,7 @@ def params_case(case, context, err_expected, err_realized, PATHS):
         pvalues = args[param]
 
     # process to be tested
-    np.random.seed(SEED)
+    legacy_seed(SEED)
     if iskfunc(Xclass):
         # test kfunc call with parameters
         x = Xclass(**args)(s, paths=PATHS)
@@ -587,7 +587,7 @@ def params_case(case, context, err_expected, err_realized, PATHS):
 
 def test_bs():
     """a bare-bone test on Black-Scholes call and put valuation"""
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     Kc, Kp = 1.2, 0.6
     T = 2.
diff --git a/sdepy/tests/test_source.py b/sdepy/tests/test_source.py
index ae16769..4a0e71b 100644
--- a/sdepy/tests/test_source.py
+++ b/sdepy/tests/test_source.py
@@ -45,7 +45,7 @@ true_wiener_source = sp.true_wiener_source
 
 # main test
 def test_source_general():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     # do cases
     paths = [10]
@@ -194,6 +194,36 @@ def source_general(paths, dtype, source_and_params):
     else:
         assert_(False)
 
+    # test rng parameter
+    try:
+        make_rngs = (
+            np.random.default_rng,
+            np.random.RandomState,
+            lambda z: np.random.Generator(np.random.PCG64(z)),
+        )
+    except AttributeError:
+        return
+    for make_rng in make_rngs:
+        rng1 = make_rng(SEED)
+        rng2 = make_rng(SEED)
+        assert rng1 is not rng2
+        src1 = cls(paths=paths, dtype=dtype, rng=rng1, **params)
+        src2 = cls(paths=paths, dtype=dtype, rng=rng2, **params)
+        assert src1.rng is rng1
+        assert src2.rng is rng2
+        s1, s2 = src1(0., 1.), src2(0., 1.)
+        assert_allclose(s1, s2)
+
+        # test that global sdepy.infrastructure.default_rng
+        # propagates correctly as a default
+        tmp = sdepy.infrastructure.default_rng
+        sdepy.infrastructure.default_rng = make_rng(SEED)
+        src3 = cls(paths=paths, dtype=dtype, rng=None, **params)
+        assert src3.rng is sdepy.infrastructure.default_rng
+        s3 = src3(0., 1.)
+        assert_allclose(s1, s3)
+        sdepy.infrastructure.default_rng = tmp
+
 
 # ---------------------
 # source specific tests
@@ -201,7 +231,7 @@ def source_general(paths, dtype, source_and_params):
 
 # main test
 def test_source_specific():
-    np.random.seed(SEED)
+    legacy_seed(SEED)
 
     # wiener tests
     src = wiener_source(vshape=3, paths=5)
@@ -226,7 +256,7 @@ def test_source_specific():
     val = 2.
 
     class const_rv:
-        def rvs(self, size):
+        def rvs(self, size, random_state):
             return np.full(size, fill_value=val)
 
     src = cpoisson_source(lam=1., paths=100, ptype=np.int16,
@@ -236,6 +266,21 @@ def test_source_specific():
     s, n = src(0, 100), src.dn_value
     assert_allclose(n*val, s, rtol=eps(s.dtype))
 
+    # compound poisson tests with legacy rv signature
+    val = 2.
+
+    class const_rv_legacy:
+        def rvs(self, size):
+            return np.full(size, fill_value=val)
+
+    src = cpoisson_source(lam=1., paths=100, ptype=np.int16,
+                              y=const_rv_legacy())
+    with assert_warns(DeprecationWarning):
+        s, n = src(0, 1), src.dn_value
+        assert_allclose(n*val, s, rtol=eps(s.dtype))
+        s, n = src(0, 100), src.dn_value
+        assert_allclose(n*val, s, rtol=eps(s.dtype))
+
     # true sources tests
     for cls in (true_wiener_source, true_poisson_source, true_cpoisson_source,
                 true_dw, true_dn, true_dj):
@@ -301,7 +346,7 @@ def test_source_true_wiener():
     rho2, irho2 = (lambda t: 0.1 + t/6, lambda t, t0: 0.1 + (t + t0)/12)
     rho3, irho3 = (None, lambda t, t0: 0)
 
-    np.random.seed(SEED)
+    legacy_seed(SEED)
     for t0 in (0, 1):
         for rho, irho in ((rho0, irho0), (rho1, irho1),
                           (rho2, irho2), (rho3, irho3)):