new file mode 100644

@@ -0,0 +1,111 @@

+from functools import partial

+import string

+

+import hypothesis.strategies as st

+from hypothesis.extra import numpy as hnp

+import numpy as np

+

+# Numbers

+

+

+def is_unit(c: complex) -> bool:

+    return np.isclose(abs(c), 1)

+

+

+MAX_QUBITS = 6

+

+phases = st.floats(

+    min_value=0, max_value=2 * np.pi, allow_nan=False, allow_infinity=False

+)

+valid_complex = st.complex_numbers(

+    max_magnitude=1e10, allow_infinity=False, allow_nan=False

+)

+unit_complex = phases.map(lambda p: np.exp(1j * p))

+nonunit_complex = valid_complex.filter(lambda c: not is_unit(c))

+

+n_qubits = st.shared(st.integers(min_value=1, max_value=MAX_QUBITS))

+

+rng = st.integers(0, 2**32 - 1).map(lambda n: np.random.default_rng(n))

+

+

+# Gates

+

+

+# Choose which qubits from 'n_qubits' to operate on with a gate that

+# operates on 'gate_size' qubits

+def select_n_qubits(gate_size: int):

+    """Return a function that, given 'n_qubits', select 'gate_size' qubits from 0-'n_qubits'."""

+

+    def _strat(n_qubits):

+        assert n_qubits >= gate_size

+        possible_qubits = st.integers(0, n_qubits - 1)

+        return st.lists(

+            possible_qubits, min_size=gate_size, max_size=gate_size, unique=True

+        ).map(tuple)

+

+    return _strat

+

+

+def unitary(n_qubits):

+    """Return a strategy for generating unitary matrices of size 2**n_qubits."""

+    size = 1 << n_qubits

+    return (

+        hnp.arrays(complex, (size, size), elements=valid_complex)

+        .map(lambda a: np.linalg.qr(a)[0])

+        .filter(lambda u: np.all(np.isfinite(u)))

+    )

+

+

+single_qubit_gates = unitary(1)

+two_qubit_gates = unitary(2)

+n_qubit_gates = n_qubits.flatmap(unitary)

+

+

+# States

+

+

+def ket(n_qubits):

+    return normalized_array(1 << n_qubits)

+

+

+def is_normalizable_to_full_precision(v):

+    # If the squared norm is > 0 but sub-normal

+    # (https://en.wikipedia.org/wiki/Subnormal_number)

+    # then normalizing the vector with v/|v| will not yield a vector

+    # that is normalized to machine precision.

+    return np.linalg.norm(v) ** 2 > np.finfo(float).smallest_normal

+

+

+def normalized_array(size):

+    return (

+        hnp.arrays(complex, (size,), elements=valid_complex)

+        .filter(is_normalizable_to_full_precision)  # vectors must be normalizable

+        .map(lambda v: v / np.linalg.norm(v))

+    )

+

+

+state_dimensions = n_qubits.map(lambda n: 2**n)

+state_shapes = state_dimensions.map(lambda x: (x,))

+

+valid_states = n_qubits.flatmap(ket)

+zero_states = state_shapes.map(lambda s: np.zeros(s, complex))

+

+invalid_state_dimensions = st.integers(0, 2**MAX_QUBITS).filter(

+    lambda n: not np.log2(n).is_integer()

+)

+invalid_shape_states = invalid_state_dimensions.flatmap(normalized_array)

+invalid_norm_states = valid_states.flatmap(

+    lambda x: nonunit_complex.map(lambda c: c * x)

+)

+

+invalid_states = st.one_of(invalid_shape_states, invalid_norm_states, zero_states)

+

+

+# Bitstrings

+

+classical_bitstrings = n_qubits.flatmap(

+    lambda n: st.integers(0, 2**n).map(partial("{:0{n}b}".format, n=n))

+)

+invalid_bitstrings = n_qubits.flatmap(

+    lambda n: st.text(alphabet=string.digits, min_size=n, max_size=n)

+).filter(lambda s: any(b not in "01" for b in s))

@@ -2,74 +2,30 @@ from functools import reduce

 from hypothesis import given

 import hypothesis.strategies as st

-import hypothesis.extra.numpy as hnp

 import numpy as np

 import pytest

 import qsim.gate

-

-# -- Strategies for generating values --

-

-

-n_qubits = st.shared(st.integers(min_value=1, max_value=6))

-

-

-# Choose which qubits from 'n_qubits' to operate on with a gate that

-# operates on 'gate_size' qubits

-def select_n_qubits(gate_size):

-    def _strat(n_qubits):

-        assert n_qubits >= gate_size

-        possible_qubits = st.integers(0, n_qubits - 1)

-        return st.lists(

-            possible_qubits, min_size=gate_size, max_size=gate_size, unique=True

-        ).map(tuple)

-

-    return _strat

-

-

-valid_complex = st.complex_numbers(

-    max_magnitude=1e10, allow_infinity=False, allow_nan=False

-)

-phases = st.floats(

-    min_value=0, max_value=2 * np.pi, allow_nan=False, allow_infinity=False

+from .common import (

+    ket,

+    single_qubit_gates,

+    n_qubit_gates,

+    valid_states,

+    n_qubits,

+    phases,

+    select_n_qubits,

 )

-def unitary(n_qubits):

-    size = 1 << n_qubits

-    return (

-        hnp.arrays(complex, (size, size), elements=valid_complex)

-        .map(lambda a: np.linalg.qr(a)[0])

-        .filter(lambda u: np.all(np.isfinite(u)))

-    )

-

-

-def ket(n_qubits):

-    size = 1 << n_qubits

-    return (

-        hnp.arrays(complex, (size,), elements=valid_complex)

-        .filter(lambda v: np.linalg.norm(v) > 0)  # vectors must be normalizable

-        .map(lambda v: v / np.linalg.norm(v))

-    )

-

-

-single_qubit_gates = unitary(1)

-two_qubit_gates = unitary(2)

-n_qubit_gates = n_qubits.flatmap(unitary)

-

-# Projectors on the single qubit computational basis

-project_zero = np.array([[1, 0], [0, 0]])

-project_one = np.array([[0, 0], [0, 1]])

-

-

 def product_gate(single_qubit_gates):

     # We reverse so that 'single_qubit_gates' can be indexed by the qubit

     # identifier; e.g. qubit #0 is actually the least-significant qubit

     return reduce(np.kron, reversed(single_qubit_gates))

-# -- Tests --

+project_zero = np.array([[1, 0], [0, 0]])

+project_one = np.array([[0, 0], [0, 1]])

 @given(n_qubits, n_qubit_gates)

@@ -145,7 +101,7 @@ def test_swap():

     assert np.all(qsim.gate.swap @ qsim.gate.swap == np.identity(4))

-@given(single_qubit_gates, n_qubits.flatmap(ket), n_qubits.flatmap(select_n_qubits(1)))

+@given(single_qubit_gates, valid_states, n_qubits.flatmap(select_n_qubits(1)))

 def test_applying_single_gates(gate, state, selected):

     (qubit,) = selected

     n_qubits = state.shape[0].bit_length() - 1

@@ -167,11 +123,13 @@ def test_applying_single_gates(gate, state, selected):

 def test_applying_controlled_single_qubit_gates(gate, state, selected):

     control, qubit = selected

     n_qubits = state.shape[0].bit_length() - 1

-    # When control qubit is |0⟩ the controlled gate acts like the identity on the other qubit

+    # When control qubit is |0⟩ the controlled gate acts

+    # like the identity on the other qubit

     parts_zero = [np.identity(2)] * n_qubits

     parts_zero[control] = project_zero

     parts_zero[qubit] = np.identity(2)

-    # When control qubit is |1⟩ the controlled gate acts like the original gate on the other qubit

+    # When control qubit is |1⟩ the controlled gate acts

+    # like the original gate on the other qubit

     parts_one = [np.identity(2)] * n_qubits

     parts_one[control] = project_one

     parts_one[qubit] = gate

@@ -1,67 +1,21 @@

-import string

-from functools import partial

-

 from hypothesis import given

-from hypothesis.extra import numpy as hnp

-import hypothesis.strategies as st

 import numpy as np

 import pytest

 import qsim.state

-

-MAX_QUBITS = 5

-

-

-def is_unit(c: complex):

-    return np.isclose(abs(c), 1)

-

-

-n_qubits = st.integers(1, MAX_QUBITS)

-state_dimensions = n_qubits.map(lambda n: 2 ** n)

-state_shapes = state_dimensions.map(lambda x: (x,))

-classical_bitstrings = n_qubits.flatmap(

-    lambda n: st.integers(0, 2 ** n).map(partial("{:0{n}b}".format, n=n))

+from .common import (

+    classical_bitstrings,

+    invalid_bitstrings,

+    invalid_states,

+    invalid_norm_states,

+    n_qubits,

+    valid_states,

 )

-invalid_bitstrings = n_qubits.flatmap(

-    lambda n: st.text(alphabet=string.digits, min_size=n, max_size=n)

-).filter(lambda s: any(b not in "01" for b in s))

-

-invalid_state_dimensions = st.integers(0, 2 ** MAX_QUBITS).filter(

-    lambda n: not np.log2(n).is_integer()

-)

-

-

-valid_complex = st.complex_numbers(

-    max_magnitude=1e10, allow_infinity=False, allow_nan=False

-)

-nonunit_complex = valid_complex.filter(lambda c: not is_unit(c))

-

-

-def ket(n_qubits):

-    return normalized_array(1 << n_qubits)

-

-

-def normalized_array(size):

-    return (

-        hnp.arrays(complex, (size,), elements=valid_complex)

-        .filter(lambda v: np.linalg.norm(v) > 0)  # vectors must be normalizable

-        .map(lambda v: v / np.linalg.norm(v))

-    )

-

-

-valid_states = n_qubits.flatmap(ket)

-invalid_shape_states = invalid_state_dimensions.flatmap(normalized_array)

-invalid_norm_states = valid_states.flatmap(

-    lambda x: nonunit_complex.map(lambda c: c * x)

-)

-zero_states = state_shapes.map(lambda s: np.zeros(s, complex))

-invalid_states = st.one_of(invalid_shape_states, invalid_norm_states, zero_states)

 @given(classical_bitstrings)

 def test_from_classical(bitstring):

-

     state = qsim.state.from_classical(bitstring)

     i = int(bitstring, base=2)

@@ -74,7 +28,6 @@ def test_from_classical(bitstring):

 @given(classical_bitstrings)

 def test_from_classical_works_on_integer_lists(bitstring):

-

     int_list = [int(b) for b in bitstring]

     assert np.all(

@@ -88,6 +41,11 @@ def test_from_classical_raises_on_bad_input(bitstring):

         qsim.state.from_classical(bitstring)

+@given(n_qubits)

+def test_zero(n):

+    assert np.array_equal(qsim.state.zero(n), qsim.state.from_classical("0" * n))

+

+

 @given(classical_bitstrings)

 def test_num_qubits(s):

     state = qsim.state.from_classical(s)