Fix seed in StatevectorSampler (#98)

* Fix seed in `StatevectorSampler`

* update doctests

* add release note

* renew release note

* remove release note

* add release note

* remove release note

povm_toolbox/post_processor/median_of_means.py

@@ -58,17 +58,18 @@ class MedianOfMeans(POVMPostProcessor):
     >>> from qiskit.circuit import QuantumCircuit
     >>> from qiskit.primitives import StatevectorSampler
     >>> from qiskit.quantum_info import SparsePauliOp
+    >>> from numpy.random import default_rng
     >>> circ = QuantumCircuit(2)
     >>> _ = circ.h(0)
     >>> _ = circ.cx(0, 1)
     >>> povm = ClassicalShadows(2, seed=42)
-    >>> sampler = StatevectorSampler(seed=42)
+    >>> sampler = StatevectorSampler(seed=default_rng(42))
     >>> povm_sampler = POVMSampler(sampler)
     >>> job = povm_sampler.run([circ], povm=povm, shots=16)
     >>> result = job.result()
     >>> post_processor = MedianOfMeans(result[0], num_batches=4, seed=42)
     >>> post_processor.get_expectation_value(SparsePauliOp("ZI"))  # doctest: +FLOAT_CMP
-    (-0.75, 2.9154759474226504)
+    (-1.6653345369377348e-16, 2.9154759474226504)
     """
 
     def __init__(

povm_toolbox/post_processor/povm_post_processor.py

@@ -37,17 +37,18 @@ class POVMPostProcessor:
     >>> from qiskit.circuit import QuantumCircuit
     >>> from qiskit.primitives import StatevectorSampler
     >>> from qiskit.quantum_info import SparsePauliOp
+    >>> from numpy.random import default_rng
     >>> circ = QuantumCircuit(2)
     >>> _ = circ.h(0)
     >>> _ = circ.cx(0, 1)
     >>> povm = ClassicalShadows(2, seed=42)
-    >>> sampler = StatevectorSampler(seed=42)
+    >>> sampler = StatevectorSampler(seed=default_rng(42))
     >>> povm_sampler = POVMSampler(sampler)
     >>> job = povm_sampler.run([circ], povm=povm, shots=16)
     >>> result = job.result()
     >>> post_processor = POVMPostProcessor(result[0])
     >>> post_processor.get_expectation_value(SparsePauliOp("ZI"))  # doctest: +FLOAT_CMP
-    (-0.75, 0.33541019662496846)
+    (-2.7755575615628914e-17, 0.3872983346207416)
 
     Additionally, this post-processor also supports the customization of the Dual frame in which the
     decomposition weights of the provided observable are obtained. Check out

povm_toolbox/sampler/povm_sampler.py

@@ -42,11 +42,12 @@ class POVMSampler:
     >>> from povm_toolbox.sampler import POVMSampler
     >>> from qiskit.circuit import QuantumCircuit
     >>> from qiskit.primitives import StatevectorSampler
+    >>> from numpy.random import default_rng
     >>> circ = QuantumCircuit(2)
     >>> _ = circ.h(0)
     >>> _ = circ.cx(0, 1)
     >>> povm = ClassicalShadows(2, seed=42)
-    >>> sampler = StatevectorSampler(seed=42)
+    >>> sampler = StatevectorSampler(seed=default_rng(42))
     >>> povm_sampler = POVMSampler(sampler)
     >>> job = povm_sampler.run([circ], povm=povm, shots=16)
     >>> result = job.result()

test/library/test_classical_shadows.py

@@ -13,6 +13,7 @@
 from unittest import TestCase
 
 import numpy as np
+from numpy.random import default_rng
 from povm_toolbox.library import ClassicalShadows
 from povm_toolbox.post_processor import POVMPostProcessor
 from povm_toolbox.quantum_info.single_qubit_povm import SingleQubitPOVM
@@ -54,7 +55,7 @@ def test_init(self):
             num_qubits,
             seed=self.SEED,
         )
-        sampler = StatevectorSampler(seed=self.SEED)
+        sampler = StatevectorSampler(seed=default_rng(self.SEED))
         povm_sampler = POVMSampler(sampler=sampler)
 
         job = povm_sampler.run([qc], shots=32, povm=measurement)

test/library/test_locally_biased_classical_shadows.py

@@ -13,6 +13,7 @@
 from unittest import TestCase
 
 import numpy as np
+from numpy.random import default_rng
 from povm_toolbox.library import LocallyBiasedClassicalShadows
 from povm_toolbox.post_processor import POVMPostProcessor
 from povm_toolbox.quantum_info.single_qubit_povm import SingleQubitPOVM
@@ -56,7 +57,7 @@ def test_init(self):
                 bias=np.array([0.2, 0.3, 0.5]),
                 seed=self.SEED,
             )
-            sampler = StatevectorSampler(seed=self.SEED)
+            sampler = StatevectorSampler(seed=default_rng(self.SEED))
             povm_sampler = POVMSampler(sampler=sampler)
 
             job = povm_sampler.run([qc], shots=32, povm=measurement)
@@ -79,7 +80,7 @@ def test_init(self):
                 bias=np.array([[0.5, 0.1, 0.4], [0.3, 0.4, 0.3]]),
                 seed=self.SEED,
             )
-            sampler = StatevectorSampler(seed=self.SEED)
+            sampler = StatevectorSampler(seed=default_rng(self.SEED))
             povm_sampler = POVMSampler(sampler=sampler)
 
             job = povm_sampler.run([qc], shots=32, povm=measurement)
@@ -137,7 +138,7 @@ def test_zero_bias(self):
         with self.subTest("Test that the POVM is not IC."):
             self.assertFalse(measurement.definition().informationally_complete)
 
-        sampler = StatevectorSampler(seed=self.SEED)
+        sampler = StatevectorSampler(seed=default_rng(self.SEED))
         povm_sampler = POVMSampler(sampler=sampler)
 
         job = povm_sampler.run([qc], shots=32, povm=measurement)
@@ -151,7 +152,7 @@ def test_zero_bias(self):
         with self.subTest("Test with compatible observable."):
             observable = SparsePauliOp(["ZZ"], coeffs=[1.0])
             exp_value, std = post_processor.get_expectation_value(observable)
-            self.assertAlmostEqual(exp_value, -0.6249999999999992)
+            self.assertAlmostEqual(exp_value, 0.6249999999999993)
             self.assertAlmostEqual(std, 0.4347552147751572)
 
         with self.subTest("Test with incompatible observable."):

test/library/test_mutually_unbiased_bases_measurements.py

@@ -40,7 +40,7 @@ def test_init(self):
                 angles=np.array([0.75, -np.pi / 3, 0.2]),
                 seed=self.SEED,
             )
-            sampler = StatevectorSampler(seed=self.SEED)
+            sampler = StatevectorSampler(seed=default_rng(self.SEED))
             povm_sampler = POVMSampler(sampler=sampler)
 
             job = povm_sampler.run([qc], shots=128, povm=measurement)
@@ -50,12 +50,12 @@ def test_init(self):
 
             observable = SparsePauliOp(["ZI"], coeffs=[1.0])
             exp_value, std = post_processor.get_expectation_value(observable)
-            self.assertAlmostEqual(exp_value, 0.8203852056374071)
-            self.assertAlmostEqual(std, 0.110717917472513)
+            self.assertAlmostEqual(exp_value, 0.6518233926221875)
+            self.assertAlmostEqual(std, 0.11921601584589436)
             observable = SparsePauliOp(["ZY"], coeffs=[1.0])
             exp_value, std = post_processor.get_expectation_value(observable)
-            self.assertAlmostEqual(exp_value, 0.5861179063820388)
-            self.assertAlmostEqual(std, 0.22424900800106443)
+            self.assertAlmostEqual(exp_value, -1.0553339936080581)
+            self.assertAlmostEqual(std, 0.2102928553415571)
 
         with self.subTest("Test specific angles for each qubit."):
             measurement = MutuallyUnbiasedBasesMeasurements(
@@ -64,7 +64,7 @@ def test_init(self):
                 angles=np.array([[1.2, 0.0, 0.4], [3.5, -0.4, 0.8]]),
                 seed=self.SEED,
             )
-            sampler = StatevectorSampler(seed=self.SEED)
+            sampler = StatevectorSampler(seed=default_rng(self.SEED))
             povm_sampler = POVMSampler(sampler=sampler)
 
             job = povm_sampler.run([qc], shots=128, povm=measurement)
@@ -74,12 +74,12 @@ def test_init(self):
 
             observable = SparsePauliOp(["ZI"], coeffs=[1.0])
             exp_value, std = post_processor.get_expectation_value(observable)
-            self.assertAlmostEqual(exp_value, 0.7697473374029422)
-            self.assertAlmostEqual(std, 0.10914516505579042)
+            self.assertAlmostEqual(exp_value, 0.7504820624371005)
+            self.assertAlmostEqual(std, 0.11019654428864213)
             observable = SparsePauliOp(["ZY"], coeffs=[1.0])
             exp_value, std = post_processor.get_expectation_value(observable)
-            self.assertAlmostEqual(exp_value, -1.5095011005732455)
-            self.assertAlmostEqual(std, 0.2080399186206604)
+            self.assertAlmostEqual(exp_value, -0.8423974419138216)
+            self.assertAlmostEqual(std, 0.23586993121676594)
 
     def test_init_errors(self):
         """Test that the ``__init__`` method raises errors correctly."""

test/library/test_povm_implementation.py

@@ -12,6 +12,7 @@
 
 from unittest import TestCase
 
+from numpy.random import default_rng
 from povm_toolbox.library import ClassicalShadows
 from qiskit.circuit import ClassicalRegister, QuantumCircuit
 from qiskit.circuit.exceptions import CircuitError
@@ -87,7 +88,7 @@ def test_errors_raised(self):
     def test_composed_circuits(self):
         """Test the composition of the input circuit with the measurement circuit."""
 
-        sampler = StatevectorSampler(seed=self.SEED)
+        sampler = StatevectorSampler(seed=default_rng(self.SEED))
 
         with self.subTest("Composed circuit."):
             pvm = ClassicalShadows(3, seed=self.SEED)
@@ -104,7 +105,7 @@ def test_composed_circuits(self):
             result = pvm._get_bitarray(job.result()[0].data)
 
             # validate outcome
-            expected = {"101": 28, "110": 66, "111": 17, "100": 17}
+            expected = {"000": 25, "001": 15, "010": 7, "100": 34, "101": 17, "110": 15, "111": 15}
             self.assertEqual(result.get_counts(), expected)
 
         with self.subTest("Composed after transpilation of input circuit."):
@@ -125,7 +126,16 @@ def test_composed_circuits(self):
             result = pvm._get_bitarray(job.result()[0].data)
 
             # validate outcome
-            expected = {"101": 68, "110": 26, "111": 17, "100": 17}
+            expected = {
+                "000": 27,
+                "001": 8,
+                "010": 2,
+                "011": 9,
+                "100": 30,
+                "101": 17,
+                "110": 21,
+                "111": 14,
+            }
             self.assertEqual(result.get_counts(), expected)
 
         with self.subTest("With a TranspileLayout present"):
@@ -152,7 +162,7 @@ def test_composed_circuits(self):
             result = pvm._get_bitarray(job.result()[0].data)
 
             # validate outcome
-            expected = {"101": 68, "110": 26, "111": 17, "100": 17}
+            expected = {"001": 15, "010": 4, "011": 5, "100": 62, "101": 17, "110": 17, "111": 8}
             self.assertEqual(result.get_counts(), expected)
 
         with self.subTest("Using measurement_layout"):
@@ -176,7 +186,7 @@ def test_composed_circuits(self):
             result = pvm._get_bitarray(job.result()[0].data)
 
             # validate outcome
-            expected = {"10": 43, "11": 85}
+            expected = {"10": 67, "11": 27, "00": 20, "01": 14}
             self.assertEqual(result.get_counts(), expected)
 
         with self.subTest("Test the insert_barriers option"):

test/post_processor/test_dual_from_empirical_frequencies.py

@@ -14,6 +14,7 @@
 from unittest import TestCase
 
 import numpy as np
+from numpy.random import default_rng
 from povm_toolbox.library import (
     ClassicalShadows,
     RandomizedProjectiveMeasurements,
@@ -59,7 +60,7 @@ def setUp(self) -> None:
         measurement = RandomizedProjectiveMeasurements(
             num_qubits, bias=bias, angles=angles, seed=self.SEED
         )
-        povm_sampler = POVMSampler(sampler=StatevectorSampler(seed=self.SEED))
+        povm_sampler = POVMSampler(sampler=StatevectorSampler(seed=default_rng(self.SEED)))
         job = povm_sampler.run([qc], shots=127, povm=measurement)
         pub_result = job.result()[0]
         self.post_processor = POVMPostProcessor(pub_result)
@@ -70,16 +71,16 @@ def test_empirical_dual(self):
 
         with self.subTest("Test canonical dual."):
             exp_value, std = self.post_processor.get_expectation_value(observable)
-            self.assertAlmostEqual(exp_value, -1.920301957227794)
-            self.assertAlmostEqual(std, 0.468920055605158)
+            self.assertAlmostEqual(exp_value, 0.7168743782790395)
+            self.assertAlmostEqual(std, 0.4724232310929573)
 
         with self.subTest("Test empirical dual with default arguments."):
             self.post_processor.dual = dual_from_empirical_frequencies(
                 povm_post_processor=self.post_processor
             )
             exp_value, std = self.post_processor.get_expectation_value(observable)
-            self.assertAlmostEqual(exp_value, -4.156136105226559)
-            self.assertAlmostEqual(std, 0.10312825405831737)
+            self.assertAlmostEqual(exp_value, 0.10431839053361033)
+            self.assertAlmostEqual(std, 0.3117383078877236)
 
         with self.subTest("Test empirical dual with lists arguments."):
             self.post_processor.dual = dual_from_empirical_frequencies(
@@ -92,8 +93,8 @@ def test_empirical_dual(self):
                 ],
             )
             exp_value, std = self.post_processor.get_expectation_value(observable)
-            self.assertAlmostEqual(exp_value, -4.16623022578294)
-            self.assertAlmostEqual(std, 0.1035227891358374)
+            self.assertAlmostEqual(exp_value, 0.1106660978854775)
+            self.assertAlmostEqual(std, 0.3118348520393772)
 
         with self.subTest(
             "Test empirical dual with `bias` and `ansatz` arguments repeated for all qubits."
@@ -105,8 +106,8 @@ def test_empirical_dual(self):
                 ansatz=SparsePauliOp(["I"], coeffs=np.array([0.5])),
             )
             exp_value, std = self.post_processor.get_expectation_value(observable)
-            self.assertAlmostEqual(exp_value, -4.0202588245449125)
-            self.assertAlmostEqual(std, 0.11724423989411743)
+            self.assertAlmostEqual(exp_value, 0.14024741778087743)
+            self.assertAlmostEqual(std, 0.3148832990774694)
 
     def test_errors_raised(self):
         """Test that the method raises the appropriate errors when suitable."""
@@ -134,7 +135,7 @@ def test_errors_raised(self):
         ):
             qc = QuantumCircuit(2)
             qc.ry(theta=Parameter("theta"), qubit=0)
-            povm_sampler = POVMSampler(sampler=StatevectorSampler(seed=self.SEED))
+            povm_sampler = POVMSampler(sampler=StatevectorSampler(seed=default_rng(self.SEED)))
             measurement = ClassicalShadows(num_qubits=2, seed=self.SEED)
             job = povm_sampler.run([(qc, np.array([0, np.pi]))], shots=16, povm=measurement)
             pub_result = job.result()[0]

test/post_processor/test_dual_from_marginal_probabilities.py

@@ -13,6 +13,7 @@
 from unittest import TestCase
 
 import numpy as np
+from numpy.random import default_rng
 from povm_toolbox.library import (
     RandomizedProjectiveMeasurements,
 )
@@ -84,7 +85,7 @@ def test_marginal_dual(self):
         measurement = RandomizedProjectiveMeasurements(
             num_qubits, bias=bias, angles=angles, seed=self.SEED
         )
-        sampler = StatevectorSampler(seed=self.SEED)
+        sampler = StatevectorSampler(seed=default_rng(self.SEED))
         povm_sampler = POVMSampler(sampler=sampler)
         job = povm_sampler.run([qc], shots=127, povm=measurement)
         pub_result = job.result()[0]
@@ -97,20 +98,20 @@ def test_marginal_dual(self):
 
         with self.subTest("Test canonical dual."):
             exp_value, std = post_processor.get_expectation_value(observable)
-            self.assertAlmostEqual(exp_value, -2.0140231870395082)
-            self.assertAlmostEqual(std, 0.664625983884081)
+            self.assertAlmostEqual(exp_value, -1.125222785367572)
+            self.assertAlmostEqual(std, 0.5305722525114711)
 
         with self.subTest("Test marginal dual."):
             post_processor.dual = dual_from_marginal_probabilities(
                 povm=post_processor.povm, state=Statevector(qc)
             )
             exp_value, std = post_processor.get_expectation_value(observable)
-            self.assertAlmostEqual(exp_value, -2.431562926033147)
-            self.assertAlmostEqual(std, 0.6951590669925843)
+            self.assertAlmostEqual(exp_value, -1.4247893013196669)
+            self.assertAlmostEqual(std, 0.4097543480314448)
 
         with self.subTest("Test threshold on marginal dual."):
             post_processor.dual = dual_from_marginal_probabilities(
                 povm=post_processor.povm, state=Statevector(qc), threshold=1.0
             )
-            self.assertAlmostEqual(exp_value, -2.4315629260331466)
-            self.assertAlmostEqual(std, 0.6951590669925845)
+            self.assertAlmostEqual(exp_value, -1.4247893013196669)
+            self.assertAlmostEqual(std, 0.4097543480314448)