Repository: jcmgray/quimb

Issue:
CircuitMPS.sample with seed returns the same bitstring for every shot
### What happened?

- Create a `CircuitMPS` that flips a coin.
- Draw 50 samples from it.
  - Do this with `seed=None`. Behavior is as expected.
  - Do this with `seed` set to some random integers. In each instance, the same bit is measured every single time.

### What did you expect to happen?

The same bit should not be returned every time. It's a coin flip.

### Minimal Complete Verifiable Example

```Python
import quimb.tensor

circ = quimb.tensor.CircuitMPS(1)
circ.apply_gate("H", 0)

shots = 50
for seed in [None, 1613, 3024, 7943]:
    print(f"seed = {seed}")
    samples = "".join(circ.sample(shots, seed=seed))
    print(f"\t{samples}")
```


### Relevant log output

```text
seed = None
	11110111110001001110000100011101110001101110101001
seed = 1613
	11111111111111111111111111111111111111111111111111
seed = 3024
	00000000000000000000000000000000000000000000000000
seed = 7943
	00000000000000000000000000000000000000000000000000
```


### Anything else we need to know?

It looks like the same Generator is being recreated before each shot:
https://github.com/jcmgray/quimb/blob/1960e054af2b15c99dfa143ac7956f5cf3b13f04/quimb/tensor/tensor_1d.py#L3209
Instead, the Generator should be created only once, before sampling, and reused for every shot.

This is probably unrelated but I noticed that here quimb sets the random number generator seed in a way that is discouraged: https://github.com/jcmgray/quimb/blob/1960e054af2b15c99dfa143ac7956f5cf3b13f04/quimb/tensor/circuit.py#L3260
Instead, a Generator should be created with `np.random.default_rng` and used.

### Environment

- Arch Linux
- Python 3.11.9
- quimb 1960e054af2b15c99dfa143ac7956f5cf3b13f04

Discussion:

Patch:
diff --git a/quimb/tensor/tensor_1d.py b/quimb/tensor/tensor_1d.py
index 3e9572d..0d3c73e 100644
--- a/quimb/tensor/tensor_1d.py
+++ b/quimb/tensor/tensor_1d.py
@@ -6,6 +6,7 @@ import operator
 from math import log, log2
 from numbers import Integral
 
+import numpy as np
 import scipy.sparse.linalg as spla
 from autoray import conj, dag, do, get_dtype_name, reshape, size, transpose
 
@@ -3266,8 +3267,9 @@ class MatrixProductState(TensorNetwork1DVector, TensorNetwork1DFlat):
         # do right canonicalization once (supplying info avoids re-performing)
         psi0 = self.canonicalize(0, info=info)
 
+        rng = np.random.default_rng(seed)
         for _ in range(C):
-            yield psi0.sample_configuration(seed=seed, info=info)
+            yield psi0.sample_configuration(seed=rng, info=info)
 
 class MatrixProductOperator(TensorNetwork1DOperator, TensorNetwork1DFlat):
     """Initialise a matrix product operator, with auto labelling and tagging.


Test patch:
diff --git a/tests/test_tensor/test_circuit.py b/tests/test_tensor/test_circuit.py
index a8a0ca5..b740f9b 100644
--- a/tests/test_tensor/test_circuit.py
+++ b/tests/test_tensor/test_circuit.py
@@ -696,6 +696,13 @@ class TestCircuitMPS:
         for x in circ.sample(10):
             assert x in {"000010", "111101"}
 
+    def test_mps_sampling_seed(self):
+        N = 1
+        circ = qtn.CircuitMPS(N)
+        circ.h(0)
+        samples = list(circ.sample(10, seed=1234))
+        assert len(set(samples)) == 2
+
     def test_permmps_sampling(self):
         N = 6
         circ = qtn.CircuitPermMPS(N)
@@ -710,6 +717,13 @@ class TestCircuitMPS:
         for x in circ.sample(10):
             assert x in {"000010", "111101"}
 
+    def test_permmps_sampling_seed(self):
+        N = 1
+        circ = qtn.CircuitPermMPS(N)
+        circ.h(0)
+        samples = list(circ.sample(10, seed=1234))
+        assert len(set(samples)) == 2
+
 
 class TestCircuitGen:
     @pytest.mark.parametrize(
diff --git a/tests/test_tensor/test_tensor_1d.py b/tests/test_tensor/test_tensor_1d.py
index 296d5dc..67b0fb2 100644
--- a/tests/test_tensor/test_tensor_1d.py
+++ b/tests/test_tensor/test_tensor_1d.py
@@ -767,6 +767,14 @@ class TestMatrixProductState:
             ).contract()
         ) ** 2 == pytest.approx(omega)
 
+    def test_sample_seed(self):
+        psi = qtn.MPS_rand_state(10, 7)
+        configs = [
+            "".join(map(str, config))
+            for config, _ in psi.sample(10, seed=1234)
+        ]
+        assert len(set(configs)) > 1
+
 
 class TestMatrixProductOperator:
     @pytest.mark.parametrize("cyclic", [False, True])