Repository: blaze/dask (python)

Issue:
random states may have same seed
The current implementation of seeding the random number generator is:

```
seeds = [self._numpy_state.randint(np.iinfo(np.int32).max) for i in sizes]
```

This can result in repetition of seeds => random sequences are the same, which is probably not what users intend. The reasoning is the same as the [birthday problem](https://en.wikipedia.org/wiki/Birthday_problem) in probability.

```
from __future__ import division

import numpy as np
from scipy.misc import comb, factorial


def prob_same(k, n=np.iinfo(np.int32).max):
    '''
    Compute the probability that at least two seeds are the same
    '''
    different = comb(n, k, exact=True) * factorial(k, exact=True)
    total = n ** k
    prob_all_different = different / total
    return 1 - prob_all_different
```

Here's a few figures:
```
In [246]: prob_same(1000)
Out[246]: 0.00023257080042193845

In [247]: prob_same(10000)
Out[247]: 0.023011865852016755

In [248]: prob_same(50000)
Out[248]: 0.4412591367423576
```

So if there are 50000 chunks, there's a 44% chance that a seed will be repeated. Probably no big deal for most use cases, but it would be better to not have to make this assumption.

It's possible to eliminate this by just sampling with replacement. Ie. change the implementation to something like:
```
seeds = self._numpy_state.choice(np.iinfo(np.int32).max, size=len(sizes), replace=True)
```

I don't know if this alternative implementation would be totally valid, since it seems like more may be required. See this [paper by L'Ecuyer](http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wsc15rng.pdf). But it would effectively guarantee different seeds.

Pull request description:
Ensure random seeds are different.
 Close #719.



Patch:
diff --git a/dask/array/random.py b/dask/array/random.py
index 6a8660a9d..76fd87cd4 100644
--- a/dask/array/random.py
+++ b/dask/array/random.py
@@ -4,6 +4,7 @@ import numpy as np
 from itertools import product
 from .core import normalize_chunks, Array
 from ..base import tokenize
+from ..utils import different_seeds
 
 def doc_wraps(func):
     """ Copy docstring from one function to another """
@@ -56,7 +57,7 @@ class RandomState(object):
 
         # Build graph
         sizes = list(product(*chunks))
-        seeds = [self._numpy_state.randint(np.iinfo(np.int32).max) for i in sizes]
+        seeds = different_seeds(len(sizes), self._numpy_state)
         token = tokenize(seeds, size, chunks, args, kwargs)
         name = 'da.random.{0}-{1}'.format(func.__name__, token)
         keys = product([name], *[range(len(bd)) for bd in chunks])
diff --git a/dask/dataframe/core.py b/dask/dataframe/core.py
index a5baae2be..bb4b5c96d 100644
--- a/dask/dataframe/core.py
+++ b/dask/dataframe/core.py
@@ -25,7 +25,7 @@ from ..array.core import partial_by_order
 from .. import threaded
 from ..compatibility import unicode, apply, operator_div, bind_method
 from ..utils import (repr_long_list, IndexCallable,
-                     pseudorandom, derived_from)
+                     pseudorandom, derived_from, different_seeds)
 from ..base import Base, compute, tokenize, normalize_token
 
 no_default = '__no_default__'
@@ -282,9 +282,17 @@ class _Frame(Base):
             columns = self.column_info
         return map_partitions(func, columns, self, *args, **kwargs)
 
-    def random_split(self, p, seed=None):
+    def random_split(self, p, random_state=None):
         """ Pseudorandomly split dataframe into different pieces row-wise
 
+        Parameters
+        ----------
+        frac : float, optional
+            Fraction of axis items to return.
+        random_state: int or np.random.RandomState
+            If int create a new RandomState with this as the seed
+        Otherwise draw from the passed RandomState
+
         Examples
         --------
 
@@ -292,12 +300,16 @@ class _Frame(Base):
 
         >>> a, b = df.random_split([0.5, 0.5])  # doctest: +SKIP
 
-        80/10/10 split, consistent seed
+        80/10/10 split, consistent random_state
+
+        >>> a, b, c = df.random_split([0.8, 0.1, 0.1], random_state=123)  # doctest: +SKIP
+
+        See Also:
+        ---------
 
-        >>> a, b, c = df.random_split([0.8, 0.1, 0.1], seed=123)  # doctest: +SKIP
+            dask.DataFrame.sample
         """
-        seeds = np.random.RandomState(seed).randint(0, np.iinfo(np.int32).max,
-                                                    self.npartitions)
+        seeds = different_seeds(self.npartitions, random_state)
         dsk_full = dict(((self._name + '-split-full', i),
                          (pd_split, (self._name, i), p, seed))
                        for i, seed in enumerate(seeds))
@@ -457,21 +469,27 @@ class _Frame(Base):
     def sample(self, frac, random_state=None):
         """ Random sample of items
 
-        This only implements the ``frac`` option from pandas.
+        Parameters
+        ----------
+        frac : float, optional
+            Fraction of axis items to return.
+        random_state: int or np.random.RandomState
+            If int create a new RandomState with this as the seed
+        Otherwise draw from the passed RandomState
 
         See Also:
-            pd.DataFrame.sample
+        ---------
+
+            dask.DataFrame.random_split, pd.DataFrame.sample
         """
+
         if random_state is None:
             random_state = np.random.randint(np.iinfo(np.int32).max)
 
         name = 'sample-' + tokenize(self, frac, random_state)
         func = getattr(self._partition_type, 'sample')
 
-        if not isinstance(random_state, np.random.RandomState):
-            random_state = np.random.RandomState(random_state)
-        seeds = random_state.randint(np.iinfo(np.int32).max,
-                                     size=self.npartitions)
+        seeds = different_seeds(self.npartitions, random_state)
 
         dsk = dict(((name, i),
                    (apply, func, (tuple, [(self._name, i)]),
@@ -2121,13 +2139,13 @@ def quantile(df, q):
     return return_type(dsk, name3, df.name, new_divisions)
 
 
-def pd_split(df, p, seed=0):
+def pd_split(df, p, random_state=None):
     """ Split DataFrame into multiple pieces pseudorandomly
 
     >>> df = pd.DataFrame({'a': [1, 2, 3, 4, 5, 6],
     ...                    'b': [2, 3, 4, 5, 6, 7]})
 
-    >>> a, b = pd_split(df, [0.5, 0.5], seed=123)  # roughly 50/50 split
+    >>> a, b = pd_split(df, [0.5, 0.5], random_state=123)  # roughly 50/50 split
     >>> a
        a  b
     1  2  3
@@ -2141,7 +2159,7 @@ def pd_split(df, p, seed=0):
     4  5  6
     """
     p = list(p)
-    index = pseudorandom(len(df), p, seed)
+    index = pseudorandom(len(df), p, random_state)
     return [df.iloc[index == i] for i in range(len(p))]
 
 
diff --git a/dask/dataframe/demo.py b/dask/dataframe/demo.py
index 2c2e94793..01df9d8bc 100644
--- a/dask/dataframe/demo.py
+++ b/dask/dataframe/demo.py
@@ -2,6 +2,7 @@ import pandas as pd
 import numpy as np
 
 from .core import tokenize, DataFrame
+from ..utils import different_seeds
 
 __all__ = ['make_timeseries']
 
@@ -69,16 +70,16 @@ def make_timeseries(start, end, dtypes, freq, partition_freq, seed=None):
     ...                              freq='2H', partition_freq='1D', seed=1)
     >>> df.head()
                            id    name     value
-    2000-01-01 00:00:00   971   Edith -0.053963
-    2000-01-01 02:00:00   937   Kevin -0.625998
-    2000-01-01 04:00:00   965  Xavier -0.922559
-    2000-01-01 06:00:00  1011   Sarah  0.866421
-    2000-01-01 08:00:00  1005   Zelda -0.319609
+    2000-01-01 00:00:00   960     Dan  0.824008
+    2000-01-01 02:00:00  1033  Xavier  0.575390
+    2000-01-01 04:00:00   986  George  0.693842
+    2000-01-01 06:00:00  1073   Sarah  0.900580
+    2000-01-01 08:00:00   976  Hannah -0.373847
     """
     divisions = list(pd.DatetimeIndex(start=start, end=end,
                                       freq=partition_freq))
     state = np.random.RandomState(seed)
-    seeds = state.randint(0, 2**31 - 1, size=len(divisions))
+    seeds = different_seeds(len(divisions), state)
     name = 'make-timeseries-' + tokenize(start, end, dtypes, freq, partition_freq)
     dsk = dict(((name, i), (make_timeseries_part, divisions[i], divisions[i + 1],
                                                  dtypes, freq, seeds[i]))
diff --git a/dask/utils.py b/dask/utils.py
index a9683e07e..0f84e589e 100644
--- a/dask/utils.py
+++ b/dask/utils.py
@@ -256,13 +256,13 @@ def skip(func):
     pass
 
 
-def pseudorandom(n, p, key):
+def pseudorandom(n, p, random_state=None):
     """ Pseudorandom array of integer indexes
 
-    >>> pseudorandom(5, [0.5, 0.5], key=123)
+    >>> pseudorandom(5, [0.5, 0.5], random_state=123)
     array([1, 0, 0, 1, 1], dtype=int8)
 
-    >>> pseudorandom(10, [0.5, 0.2, 0.2, 0.1], key=5)
+    >>> pseudorandom(10, [0.5, 0.2, 0.2, 0.1], random_state=5)
     array([0, 2, 0, 3, 0, 1, 2, 1, 0, 0], dtype=int8)
     """
     import numpy as np
@@ -271,7 +271,10 @@ def pseudorandom(n, p, key):
     assert np.allclose(1, cp[-1])
     assert len(p) < 256
 
-    x = np.random.RandomState(key).random_sample(n)
+    if not isinstance(random_state, np.random.RandomState):
+        random_state = np.random.RandomState(random_state)
+
+    x = random_state.random_sample(n)
     out = np.empty(n, dtype='i1')
 
     for i, (low, high) in enumerate(zip(cp[:-1], cp[1:])):
@@ -279,6 +282,32 @@ def pseudorandom(n, p, key):
     return out
 
 
+def different_seeds(n, random_state=None):
+    """ A list of different 32 bit integer seeds
+
+    Parameters
+    ----------
+    n: int
+        Number of distinct seeds to return
+    random_state: int or np.random.RandomState
+        If int create a new RandomState with this as the seed
+    Otherwise draw from the passed RandomState
+    """
+    import numpy as np
+
+    if not isinstance(random_state, np.random.RandomState):
+        random_state = np.random.RandomState(random_state)
+
+    big_n = np.iinfo(np.int32).max
+
+    seeds = set(random_state.randint(big_n, size=n))
+    while len(seeds) < n:
+        seeds.add(random_state.randint(big_n))
+
+    # Sorting makes it easier to know what seeds are for what chunk
+    return sorted(seeds)
+
+
 def is_integer(i):
     """
     >>> is_integer(6)


Test patch:
diff --git a/dask/array/tests/test_linalg.py b/dask/array/tests/test_linalg.py
index ebce5cd07..0c9e86e7f 100644
--- a/dask/array/tests/test_linalg.py
+++ b/dask/array/tests/test_linalg.py
@@ -103,7 +103,7 @@ def test_linalg_consistent_names():
 
 
 def test_svd_compressed():
-    m, n = 300, 250
+    m, n = 500, 250
     r = 10
     np.random.seed(4321)
     mat1 = np.random.randn(m, r)
@@ -111,23 +111,16 @@ def test_svd_compressed():
     mat = mat1.dot(mat2)
     data = da.from_array(mat, chunks=(50, 50))
 
-    n_iter = 6
-    for i in range(n_iter):
-        u, s, vt = svd_compressed(data, r, seed=4321)
-        u, s, vt = da.compute(u, s, vt)
-        if i == 0:
-            usvt = np.dot(u, np.dot(np.diag(s), vt))
-        else:
-            usvt += np.dot(u, np.dot(np.diag(s), vt))
-    usvt /= n_iter
+    u, s, vt = svd_compressed(data, r, seed=4321)
+    u, s, vt = da.compute(u, s, vt)
 
-    tol = 2e-1
+    usvt = np.dot(u, np.dot(np.diag(s), vt))
+
+    tol = 0.2
     assert np.allclose(np.linalg.norm(mat - usvt),
                        np.linalg.norm(mat),
                        rtol=tol, atol=tol)  # average accuracy check
 
-    u, s, vt = svd_compressed(data, r, seed=4321)
-    u, s, vt = da.compute(u, s, vt)
     u = u[:, :r]
     s = s[:r]
     vt = vt[:r, :]
diff --git a/dask/tests/test_utils.py b/dask/tests/test_utils.py
index 1b7701117..df9e0f69d 100644
--- a/dask/tests/test_utils.py
+++ b/dask/tests/test_utils.py
@@ -1,9 +1,10 @@
 import os
 from itertools import product
 
+import numpy as np
 
 from dask.utils import (textblock, filetext, takes_multiple_arguments,
-                        Dispatch, tmpfile, next_linesep)
+                        Dispatch, tmpfile, next_linesep, different_seeds)
 
 
 def test_textblock():
@@ -113,3 +114,21 @@ def test_gh606():
             stop = len(bin_euro) * 10 + len(bin_linesep)
             res = textblock(f, 0, stop, encoding=encoding)
             assert res == ((euro * 10) + linesep + (yen * 10) + linesep).encode(encoding)
+
+
+def test_different_seeds():
+    seed = 37
+    state = np.random.RandomState(seed)
+    n = 100000
+
+    # Use an integer
+    seeds = set(different_seeds(n, seed))
+    assert len(seeds) == n
+
+    # Use RandomState object
+    seeds2 = set(different_seeds(n, state))
+    assert seeds == seeds2
+
+    # Should be sorted
+    smallseeds = different_seeds(10, 1234)
+    assert smallseeds == sorted(smallseeds)