Repository: dask/dask

Issue:
Unable to replicate random samples between dask and numpy
<!-- Please include a self-contained copy-pastable example that generates the issue if possible.

Please be concise with code posted. See guidelines below on how to provide a good bug report:

- Craft Minimal Bug Reports http://matthewrocklin.com/blog/work/2018/02/28/minimal-bug-reports
- Minimal Complete Verifiable Examples https://stackoverflow.com/help/mcve

Bug reports that follow these guidelines are easier to diagnose, and so are often handled much more quickly.
-->

**What happened**:
When the same random seed is used for both dask and numpy, the output of random.normal() do not match each other. 

**What you expected to happen**:
I would have expected the output to match verbatim. When setting a dask random state, and then calling `<random_state>._numpy_state.normal()` the results will match those from numpy. However, this returns a np.ndarray() versus a dask array. I am not sure if this is the intended behavior. 

**Minimal Complete Verifiable Example**:
Using depreciated random.seed():
```python
import dask.array as da
da.random.seed(42)
da.random.normal(0, 1, size=10).compute()

array([-1.31918189,  0.11348963,  1.2427626 , -1.26858024,  2.27840845,
        0.57524507, -1.14656226,  0.90500015, -1.17130718, -0.89184311])

import numpy as np
np.random.seed(42)
np.random.normal(0, 1 ,size=10)

array([ 0.49671415, -0.1382643 ,  0.64768854,  1.52302986, -0.23415337,
       -0.23413696,  1.57921282,  0.76743473, -0.46947439,  0.54256004])
```
Using RandomState():
```python
import dask.array as da
rnd = da.random.RandomState(42)
rnd.normal(0, 1, size=10).compute()

array([-1.31918189,  0.11348963,  1.2427626 , -1.26858024,  2.27840845,
        0.57524507, -1.14656226,  0.90500015, -1.17130718, -0.89184311])

import numpy as np
nprnd = np.random.RandomState(42)
nprnd.normal(0, 1, size=10)

array([ 0.49671415, -0.1382643 ,  0.64768854,  1.52302986, -0.23415337,
       -0.23413696,  1.57921282,  0.76743473, -0.46947439,  0.54256004])
```

Numpy array returned when `_numpy_state` is used (this may be intentional) :
```python
import dask.array as da
rnd = da.random.RandomState(42)
foo = rnd._numpy_state.normal(0, 1, size=10)
type(foo)

numpy.ndarray
```

**Anything else we need to know?**:

**Environment**:

- Dask version: 2021.12.0
- Python version: 3.9
- Operating System: Ubuntu 20.04.3 LTS (Focal Fossa)
- Install method (conda, pip, source): conda

<!-- If you are reporting an issue such as scale stability, cluster deadlock.
Please provide a cluster dump state with this issue, by running client.dump_cluster_state()

https://distributed.dask.org/en/stable/api.html?highlight=dump_cluster_state#distributed.Client.dump_cluster_state

-->

<details>
<summary>Cluster Dump State:</summary>

</details>

Discussion:
erayaslan: good comments. thank you @rjzamora  Done!

Patch:
diff --git a/dask/array/random.py b/dask/array/random.py
index a98959177..a96e85969 100644
--- a/dask/array/random.py
+++ b/dask/array/random.py
@@ -28,15 +28,16 @@ class Generator:
     Container for the BitGenerators.
 
     ``Generator`` exposes a number of methods for generating random
-    numbers drawn from a variety of probability distributions. In addition to
-    the distribution-specific arguments, each method takes a keyword argument
-    `size` that defaults to ``None``. If `size` is ``None``, then a single
-    value is generated and returned. If `size` is an integer, then a 1-D
-    array filled with generated values is returned. If `size` is a tuple,
-    then an array with that shape is filled and returned.
+    numbers drawn from a variety of probability distributions and serves
+    as a replacement for ``RandomState``. The main difference between the
+    two is that ``Generator`` relies on an additional ``BitGenerator`` to
+    manage state and generate the random bits, which are then transformed
+    into random values from useful distributions. The default ``BitGenerator``
+    used by ``Generator`` is ``PCG64``. The ``BitGenerator`` can be changed
+    by passing an instantiated ``BitGenerator`` to ``Generator``.
 
-    The function :func:`dask.array.random.default_rng` will instantiate
-    a `Generator` with NumPy's default `BitGenerator`.
+    The function :func:`dask.array.random.default_rng` is the recommended way
+    to instantiate a ``Generator``.
 
     .. warning::
 
@@ -45,9 +46,6 @@ class Generator:
        ``Generator`` does not provide a version compatibility guarantee. In
        particular, as better algorithms evolve the bit stream may change.
 
-       This object is identical to ``np.random.Generator`` except that all
-       functions also take a ``chunks=`` keyword argument.
-
     Parameters
     ----------
     bit_generator : BitGenerator
@@ -55,6 +53,13 @@ class Generator:
 
     Notes
     -----
+    In addition to the distribution-specific arguments, each ``Generator``
+    method takes a keyword argument `size` that defaults to ``None``. If
+    `size` is ``None``, then a single value is generated and returned. If
+    `size` is an integer, then a 1-D array filled with generated values is
+    returned. If `size` is a tuple, then an array with that shape is
+    filled and returned.
+
     The Python stdlib module `random` contains pseudo-random number generator
     with a number of methods that are similar to the ones available in
     ``Generator``. It uses Mersenne Twister, and this bit generator can
@@ -62,6 +67,13 @@ class Generator:
     Dask-aware, has the advantage that it provides a much larger number
     of probability distributions to choose from.
 
+    All ``Generator`` methods are identical to ``np.random.Generator`` except
+    that they also take a `chunks=` keyword argument.
+
+    ``Generator`` does not guarantee parity in the generated numbers
+    with any third party library. In particular, numbers generated by
+    `Dask` and `NumPy` will differ even if they use the same seed.
+
     Examples
     --------
     >>> from numpy.random import PCG64
@@ -387,8 +399,8 @@ def default_rng(seed=None):
         A seed to initialize the `BitGenerator`. If None, then fresh,
         unpredictable entropy will be pulled from the OS. If an ``int`` or
         ``array_like[ints]`` is passed, then it will be passed to
-        `SeedSequence` to derive the initial `BitGenerator` state. One may also
-        pass in a `SeedSequence` instance.
+        `SeedSequence` to derive the initial `BitGenerator` state. One may
+        also pass in a `SeedSequence` instance.
         Additionally, when passed a `BitGenerator`, it will be wrapped by
         `Generator`. If passed a `Generator`, it will be returned unaltered.
 
@@ -399,13 +411,14 @@ def default_rng(seed=None):
 
     Notes
     -----
-    If ``seed`` is not a `BitGenerator` or a `Generator`, a new `BitGenerator`
-    is instantiated. This function does not manage a default global instance.
+    If ``seed`` is not a `BitGenerator` or a `Generator`, a new
+    `BitGenerator` is instantiated. This function does not manage a default
+    global instance.
 
     Examples
     --------
-    ``default_rng`` is the recommended constructor for the random number class
-    ``Generator``. Here are several ways we can construct a random
+    ``default_rng`` is the recommended constructor for the random number
+    class ``Generator``. Here are several ways we can construct a random
     number generator using ``default_rng`` and the ``Generator`` class.
 
     Here we use ``default_rng`` to generate a random float:
diff --git a/dask/array/ufunc.py b/dask/array/ufunc.py
index ffb59fb9c..8554edf3e 100644
--- a/dask/array/ufunc.py
+++ b/dask/array/ufunc.py
@@ -327,7 +327,7 @@ def modf(x):
         for key in core.flatten(tmp.__dask_keys__())
     }
 
-    a = np.empty_like(getattr(x, "_meta", x), shape=(1,) * x.ndim, dtype=x.dtype)
+    a = np.ones_like(getattr(x, "_meta", x), shape=(1,) * x.ndim, dtype=x.dtype)
     l, r = np.modf(a)
 
     graph = HighLevelGraph.from_collections(left, ldsk, dependencies=[tmp])
diff --git a/dask/config.py b/dask/config.py
index c3a1836a7..6aee66f17 100644
--- a/dask/config.py
+++ b/dask/config.py
@@ -79,7 +79,12 @@ def canonical_name(k: str, config: dict) -> str:
     return k
 
 
-def update(old: dict, new: Mapping, priority: Literal["old", "new"] = "new") -> dict:
+def update(
+    old: dict,
+    new: Mapping,
+    priority: Literal["old", "new", "new-defaults"] = "new",
+    defaults: Mapping | None = None,
+) -> dict:
     """Update a nested dictionary with values from another
 
     This is like dict.update except that it smoothly merges nested values
@@ -88,9 +93,12 @@ def update(old: dict, new: Mapping, priority: Literal["old", "new"] = "new") ->
 
     Parameters
     ----------
-    priority: string {'old', 'new'}
+    priority: string {'old', 'new', 'new-defaults'}
         If new (default) then the new dictionary has preference.
         Otherwise the old dictionary does.
+        If 'new-defaults', a mapping should be given of the current defaults.
+        Only if a value in ``old`` matches the current default, it will be
+        updated with ``new``.
 
     Examples
     --------
@@ -104,6 +112,12 @@ def update(old: dict, new: Mapping, priority: Literal["old", "new"] = "new") ->
     >>> update(a, b, priority='old')  # doctest: +SKIP
     {'x': 1, 'y': {'a': 2, 'b': 3}}
 
+    >>> d = {'x': 0, 'y': {'a': 2}}
+    >>> a = {'x': 1, 'y': {'a': 2}}
+    >>> b = {'x': 2, 'y': {'a': 3, 'b': 3}}
+    >>> update(a, b, priority='new-defaults', defaults=d)  # doctest: +SKIP
+    {'x': 1, 'y': {'a': 3, 'b': 3}}
+
     See Also
     --------
     dask.config.merge
@@ -114,9 +128,23 @@ def update(old: dict, new: Mapping, priority: Literal["old", "new"] = "new") ->
         if isinstance(v, Mapping):
             if k not in old or old[k] is None:
                 old[k] = {}
-            update(old[k], v, priority=priority)
+            update(
+                old[k],
+                v,
+                priority=priority,
+                defaults=defaults.get(k) if defaults else None,
+            )
         else:
-            if priority == "new" or k not in old:
+            if (
+                priority == "new"
+                or k not in old
+                or (
+                    priority == "new-defaults"
+                    and defaults
+                    and k in defaults
+                    and defaults[k] == old[k]
+                )
+            ):
                 old[k] = v
 
     return old
@@ -568,11 +596,13 @@ def update_defaults(
     It does two things:
 
     1.  Add the defaults to a global collection to be used by refresh later
-    2.  Updates the global config with the new configuration
-        prioritizing older values over newer ones
+    2.  Updates the global config with the new configuration.
+        Old values are prioritized over new ones, unless the current value
+        is the old default, in which case it's updated to the new default.
     """
+    current_defaults = merge(*defaults)
     defaults.append(new)
-    update(config, new, priority="old")
+    update(config, new, priority="new-defaults", defaults=current_defaults)
 
 
 def expand_environment_variables(config: Any) -> Any:
diff --git a/docs/source/array-chunks.rst b/docs/source/array-chunks.rst
index 0d98d94dd..f490d6295 100644
--- a/docs/source/array-chunks.rst
+++ b/docs/source/array-chunks.rst
@@ -65,7 +65,8 @@ following:
 
 .. code-block:: python
 
-   >>> x = da.from_array(np.random.randn(100), chunks=20)
+   >>> rng = np.random.default_rng()
+   >>> x = da.from_array(rng.standard_normal(100), chunks=20)
    >>> x += 0.1
    >>> y = x[x > 0]  # don't know how many values are greater than 0 ahead of time
 
diff --git a/docs/source/array-creation.rst b/docs/source/array-creation.rst
index 504d90fb3..34a2aca9a 100644
--- a/docs/source/array-creation.rst
+++ b/docs/source/array-creation.rst
@@ -63,7 +63,8 @@ list see the :doc:`Array API <array-api>`:
 .. code-block:: python
 
    >>> import dask.array as da
-   >>> x = da.random.random((10000, 10000), chunks=(1000, 1000))
+   >>> rng = da.random.default_rng()
+   >>> x = rng.random((10000, 10000), chunks=(1000, 1000))
 
 
 Concatenation and Stacking
diff --git a/docs/source/array-gufunc.rst b/docs/source/array-gufunc.rst
index 9b00e2646..1dbbe507a 100644
--- a/docs/source/array-gufunc.rst
+++ b/docs/source/array-gufunc.rst
@@ -25,7 +25,7 @@ NumPy Generalized UFuncs
     import dask.array as da
     import numpy as np
 
-    x = da.random.normal(size=(3, 10, 10), chunks=(2, 10, 10))
+    x = da.random.default_rng().normal(size=(3, 10, 10), chunks=(2, 10, 10))
 
     w, v = np.linalg._umath_linalg.eig(x, output_dtypes=(float, float))
 
@@ -47,7 +47,7 @@ Wrap your own Python function
 
 .. code-block:: python
 
-    x = da.random.normal(size=(10, 5), chunks=(2, 5))
+    x = da.random.default_rng().normal(size=(10, 5), chunks=(2, 5))
 
     def foo(x):
         return np.mean(x, axis=-1)
diff --git a/docs/source/array-random.rst b/docs/source/array-random.rst
new file mode 100644
index 000000000..a3ae71d73
--- /dev/null
+++ b/docs/source/array-random.rst
@@ -0,0 +1,76 @@
+Random Number Generation
+========================
+
+Dask's random number routines produce pseudo random numbers using combinations
+of a ``BitGenerator`` to create sequences and a ``Generator`` to use those
+sequences to sample from different statistical distributions:
+
+    **BitGenerators**: Objects that generate random sequences. These are
+    provided by a backend library such as NumPy or CuPy and are typically
+    unsigned integer words filled with sequences of either 32 or 64 random
+    bits.
+
+    ``Generators:`` Objects that transform sequences of random bits from a
+    ``BitGenerator`` into sequences of numbers that follow a specific probability
+    distribution (such as uniform, Normal or Binomial) within a specified
+    interval.
+
+Since Dask version 2023.2.1 the ``Generator`` can be initialized with a number
+of different ``BitGenerator`` classes. It exposes many different probability
+distributions. The legacy ``RandomState`` random number routines are still
+available, but are considered frozen and will not be getting any updates.
+
+Differences with NumPy
+----------------------
+
+Dask follows the NumPy interface for random number generation with some
+differences:
+
+- Methods under ``dask.array.random`` take a ``chunks`` keyword.
+- Dask tries to be backend agnostic. In other words, you can mostly use CuPy
+  and NumPy interchangeably as a backend for random number generation. Any
+  library providing a similar interface should also work with some effort.
+
+Notes
+-----
+
+- Dask does not guarantee that the same number generator is used across versions.
+  This means that numbers generated by ``dask.array.random`` by a new version may
+  not be the same as the previous one, even when the same seed and distribution
+  are used. As better algorithms evolve, the bit stream may change.
+- Dask does not guarantee parity in the generated numbers with any third party
+  library. In particular, numbers generated by Dask and NumPy or CuPy will differ
+  even when given the same seed and ``BitGenerator``. Dask tends to spawn ``SeedSequence``
+  children to produce independent random number streams in parallel.
+- ``Generator.integers`` is now the canonical way to generate integer random numbers
+  from a discrete uniform distribution. The `endpoint` keyword can be used to
+  specify open or closed intervals. This replaces both `randint` and `random_integers`.
+- ``Generator.random`` is now the canonical way to generate floating-point random
+  numbers, which replaces `random_sample`. The ``dask.array.random.random``
+  method still uses ``RandomState`` for backwards compatibility and should be
+  avoided for new code. Please use ``Generator.random`` instead.
+
+Quick Start
+-----------
+
+Call ``default_rng`` to get a new instance of a Generator, then call its methods to
+obtain samples from different distributions. By default, ``Generator`` uses bits
+provided by PCG64 which has better statistical properties than the legacy MT19937
+used in ``RandomState``.
+
+.. code-block:: python
+
+    # Do this (new version)
+    import dask.array as da
+    rng = da.random.default_rng()
+    vals = rng.standard_normal(10)
+    more_vals = rng.standard_normal(10)
+
+    # instead of this (legacy version)
+    import dask.array as da
+    vals = da.random.standard_normal(10)
+    more_vals = da.random.standard_normal(10)
+
+
+
+For further info, please see `NumPy docs <https://numpy.org/devdocs/reference/random/index.html>`_ 
diff --git a/docs/source/array-sparse.rst b/docs/source/array-sparse.rst
index f92984e38..6eb3f2dc3 100644
--- a/docs/source/array-sparse.rst
+++ b/docs/source/array-sparse.rst
@@ -19,7 +19,8 @@ Say we have a Dask array with mostly zeros:
 
 .. code-block:: python
 
-   x = da.random.random((100000, 100000), chunks=(1000, 1000))
+   rng = da.random.default_rng()
+   x = rng.random((100000, 100000), chunks=(1000, 1000))
    x[x < 0.95] = 0
 
 We can convert each of these chunks of NumPy arrays into a **sparse.COO** array:
diff --git a/docs/source/array-stats.rst b/docs/source/array-stats.rst
index 6b7f6400e..ad630a7b7 100644
--- a/docs/source/array-stats.rst
+++ b/docs/source/array-stats.rst
@@ -11,7 +11,8 @@ You can calculate various measures of an array including skewness, kurtosis, and
 .. code-block:: python
 
    >>> from dask.array import stats
-   >>> x = da.random.beta(1, 1, size=(1000,), chunks=10)
+   >>> rng = da.random.default_rng()
+   >>> x = rng.beta(1, 1, size=(1000,), chunks=10)
    >>> k, s, m = [stats.kurtosis(x), stats.skew(x), stats.moment(x, 5)]
    >>> dask.compute(k, s, m)
    (1.7612340817172787, -0.064073498030693302, -0.00054523780628304799)
@@ -27,8 +28,9 @@ results.
 
 .. code-block:: python
 
-   >>> a = da.random.uniform(size=(50,), chunks=(25,))
-   >>> b = a + da.random.uniform(low=-0.15, high=0.15, size=(50,), chunks=(25,))
+   >>> rng = da.random.default_rng()
+   >>> a = rng.uniform(size=(50,), chunks=(25,))
+   >>> b = a + rng.uniform(low=-0.15, high=0.15, size=(50,), chunks=(25,))
    >>> result = stats.ttest_rel(a, b)
    >>> result.compute()
    Ttest_relResult(statistic=-1.5102104380013242, pvalue=0.13741197274874514)
diff --git a/docs/source/array.rst b/docs/source/array.rst
index c7b903c36..e301285d1 100644
--- a/docs/source/array.rst
+++ b/docs/source/array.rst
@@ -16,6 +16,7 @@ Array
    array-assignment.rst
    array-stack.rst
    array-gufunc.rst
+   array-random.rst
    array-api.rst
 
 Dask Array implements a subset of the NumPy ndarray interface using blocked
diff --git a/docs/source/develop.rst b/docs/source/develop.rst
index f7c6138f2..5161e9c6b 100644
--- a/docs/source/develop.rst
+++ b/docs/source/develop.rst
@@ -198,7 +198,8 @@ typically tested directly against the NumPy or pandas libraries using the
    from dask.array.utils import assert_eq
 
    def test_aggregations():
-       nx = np.random.random(100)
+       rng = np.random.default_rng()
+       nx = rng.random(100)
        dx = da.from_array(nx, chunks=(10,))
 
        assert_eq(nx.sum(), dx.sum())
diff --git a/docs/source/diagnostics-local.rst b/docs/source/diagnostics-local.rst
index 2ebceded2..2a12b426b 100644
--- a/docs/source/diagnostics-local.rst
+++ b/docs/source/diagnostics-local.rst
@@ -35,7 +35,7 @@ computation:
 
     >>> import dask.array as da
     >>> from dask.diagnostics import ProgressBar
-    >>> a = da.random.normal(size=(10000, 10000), chunks=(1000, 1000))
+    >>> a = da.random.default_rng().normal(size=(10000, 10000), chunks=(1000, 1000))
     >>> res = a.dot(a.T).mean(axis=0)
 
     >>> with ProgressBar():
diff --git a/docs/source/futures.rst b/docs/source/futures.rst
index 178e111df..5fddeafb7 100644
--- a/docs/source/futures.rst
+++ b/docs/source/futures.rst
@@ -843,7 +843,7 @@ The client will calculate the gradient of the loss function above.
    ps_future = client.submit(ParameterServer, actor=True)
    ps = ps_future.result()
 
-   ps.put('parameters', np.random.random(1000))
+   ps.put('parameters', np.random.default_rng().random(1000))
    for k in range(20):
        params = ps.get('parameters').result()
        new_params = train(params)
diff --git a/docs/source/graph_manipulation.rst b/docs/source/graph_manipulation.rst
index 3953382ed..3d225443f 100644
--- a/docs/source/graph_manipulation.rst
+++ b/docs/source/graph_manipulation.rst
@@ -19,7 +19,7 @@ Consider the following example:
 .. code-block:: python
 
    >>> import dask.array as da
-   >>> x = da.random.normal(size=500_000_000, chunks=100_000)
+   >>> x = da.random.default_rng().normal(size=500_000_000, chunks=100_000)
    >>> x_mean = x.mean()
    >>> y = (x - x_mean).max().compute()
 
diff --git a/docs/source/user-interfaces.rst b/docs/source/user-interfaces.rst
index ba512c8a3..a85beca6d 100644
--- a/docs/source/user-interfaces.rst
+++ b/docs/source/user-interfaces.rst
@@ -38,8 +38,9 @@ large subset of the API from the original project.
 
    # Arrays
    import dask.array as da
-   x = da.random.uniform(low=0, high=10, size=(10000, 10000),  # normal numpy code
-                         chunks=(1000, 1000))  # break into chunks of size 1000x1000
+   rng = da.random.default_rng()
+   x = rng.uniform(low=0, high=10, size=(10000, 10000),  # normal numpy code
+                   chunks=(1000, 1000))  # break into chunks of size 1000x1000
 
    y = x + x.T - x.mean(axis=0)  # Use normal syntax for high level algorithms
 


Test patch:
diff --git a/dask/tests/test_config.py b/dask/tests/test_config.py
index 2d8836b99..762fccf1d 100644
--- a/dask/tests/test_config.py
+++ b/dask/tests/test_config.py
@@ -52,6 +52,42 @@ def test_update():
     assert b == {"x": 2, "y": {"a": 3, "b": 2}, "z": 3}
 
 
+def test_update_new_defaults():
+    d = {"x": 1, "y": 1, "z": {"a": 1, "b": 1}}
+    o = {"x": 1, "y": 2, "z": {"a": 1, "b": 2}, "c": 2, "c2": {"d": 2}}
+    n = {"x": 3, "y": 3, "z": OrderedDict({"a": 3, "b": 3}), "c": 3, "c2": {"d": 3}}
+    assert update(o, n, priority="new-defaults", defaults=d) == {
+        "x": 3,
+        "y": 2,
+        "z": {"a": 3, "b": 2},
+        "c": 2,
+        "c2": {"d": 2},
+    }
+    assert update(o, n, priority="new-defaults", defaults=o) == update(
+        o, n, priority="new"
+    )
+    assert update(o, n, priority="new-defaults", defaults=None) == update(
+        o, n, priority="old"
+    )
+
+
+def test_update_defaults():
+    defaults = [
+        {"a": 1, "b": {"c": 1}},
+        {"a": 2, "b": {"d": 2}},
+    ]
+    current = {"a": 2, "b": {"c": 1, "d": 3}, "extra": 0}
+    new = {"a": 0, "b": {"c": 0, "d": 0}, "new-extra": 0}
+    update_defaults(new, current, defaults=defaults)
+
+    assert defaults == [
+        {"a": 1, "b": {"c": 1}},
+        {"a": 2, "b": {"d": 2}},
+        {"a": 0, "b": {"c": 0, "d": 0}, "new-extra": 0},
+    ]
+    assert current == {"a": 0, "b": {"c": 0, "d": 3}, "extra": 0, "new-extra": 0}
+
+
 def test_merge():
     a = {"x": 1, "y": {"a": 1}}
     b = {"x": 2, "z": 3, "y": {"b": 2}}