Merge pull request #124 from UCL-CCS/utils_tests

updates to symmer utils and tests added. Not poetry version needed to be fixed to 1.4.0 in github workflow.

.github/workflows/pull_request.yaml

@@ -31,7 +31,8 @@ jobs:
           pip install .
       - name: Poetry install package
         run: |
-          pip install poetry
+          # poetry 1.4.1 has problems with debugpy for some reason.
+          pip install poetry==1.4.0
           poetry install
       - name: Linting
         run: |

symmer/operators/independent_op.py

@@ -88,7 +88,9 @@ def symmetry_generators(cls,
         S_symp = cref_matrix[PwordOp.n_terms:,np.all(~cref_matrix[:PwordOp.n_terms], axis=0)].T
         S = cls(S_symp, np.ones(S_symp.shape[0]))
         if S.n_terms==0:
-            raise RuntimeError('The input PauliwordOp has no Z2 symmetries.')
+            warnings.warn('The input PauliwordOp has no Z2 symmetries.')
+            return S
+            # raise RuntimeError('The input PauliwordOp has no Z2 symmetries.')
         if np.all(S.adjacency_matrix) or commuting_override:
             return S
         else:

symmer/utils.py

@@ -1,7 +1,7 @@
 from symmer.operators import PauliwordOp, QuantumState
 import numpy as np
 import scipy as sp
-from typing import Union, List, Tuple
+from typing import List, Tuple, Union
 from functools import reduce
 import py3Dmol
 
@@ -48,51 +48,80 @@ def exact_gs_energy(
                     )
                 )
                 expval_n_particle += Z_coeff * np.sum(sign * np.square(abs(psi.state_op.coeff_vec)))
-            if round(expval_n_particle) == n_particles:
+            if np.round(expval_n_particle) == n_particles:
                 return evl, QuantumState.from_array(evc.reshape([-1,1]))
         # if a solution is not found within the first n_eig eigenvalues then error
         raise RuntimeError('No eigenvector of the correct particle number was identified - try increasing n_eigs.')
 
 
-def random_anitcomm_2n_1_PauliwordOp(n_qubits, complex_coeff=True, apply_unitary=True):
+def random_anitcomm_2n_1_PauliwordOp(n_qubits, complex_coeff=False, apply_clifford=True):
     """ Generate a anticommuting PauliOperator of size 2n+1 on n qubits (max possible size)
         with normally distributed coefficients. Generates in structured way then uses Clifford rotation (default)
         to try and make more random (can stop this to allow FAST build, but inherenet structure
          will be present as operator is formed in specific way!)
     """
-    base = 'X' * n_qubits
-    I_term = 'I' * n_qubits
+    # base = 'X' * n_qubits
+    # I_term = 'I' * n_qubits
+    # P_list = [base]
+    # for i in range(n_qubits):
+    #     # Z_term
+    #     P_list.append(base[:i] + 'Z' + I_term[i + 1:])
+    #     # Y_term
+    #     P_list.append(base[:i] + 'Y' + I_term[i + 1:])
+    # coeff_vec = np.random.randn(len(P_list)).astype(complex)
+    # if complex_coeff:
+    #     coeff_vec += 1j * np.random.randn((len(P_list)))
+    # P_anticomm = PauliwordOp.from_dictionary((dict(zip(P_list, coeff_vec))))
+
+    Y_base = np.hstack((np.eye(n_qubits), np.tril(np.ones(n_qubits))))
+    X_base = Y_base.copy()
+    X_base[:, n_qubits:] = np.tril(np.ones(n_qubits), -1)
+
+    ac_symp = np.vstack((Y_base, X_base))
+
+    Z_symp = np.zeros(2 * n_qubits)
+    Z_symp[n_qubits:] = np.ones(n_qubits)
+
+    ac_symp = np.vstack((ac_symp, Z_symp)).astype(bool)
+
+    coeff_vec = np.random.randn(ac_symp.shape[0]).astype(complex)
+    if complex_coeff:
+        coeff_vec += 1j * np.random.randn(2 * n_qubits + 1).astype(complex)
 
-    P_list = [base]
-    for i in range(n_qubits):
-        # Z_term
-        P_list.append(base[:i] + 'Z' + I_term[i + 1:])
-        # Y_term
-        P_list.append(base[:i] + 'Y' + I_term[i + 1:])
+    P_anticomm = PauliwordOp(ac_symp, coeff_vec)
 
-    coeff_vec = np.random.randn(len(P_list)).astype(complex)
-    if complex_coeff:
-        coeff_vec += 1j * np.random.randn((len(P_list)))
+    if apply_clifford:
+        # apply clifford rotations to get rid of structure
+        U_cliff_rotations = []
+        for _ in range(n_qubits + 5):
+            P_rand = PauliwordOp.random(n_qubits, n_terms=1)
+            P_rand.coeff_vec = [1]
+            U_cliff_rotations.append((P_rand, None))
 
-    P_anticomm = PauliwordOp.from_dictionary((dict(zip(P_list, coeff_vec))))
+        P_anticomm = P_anticomm.perform_rotations(U_cliff_rotations)
 
-    # random rotations to get rid of structure
-    if apply_unitary:
-        U = PauliwordOp.haar_random(n_qubits=n_qubits)
-        P_anticomm = U * P_anticomm * U.dagger
+    assert P_anticomm.n_terms == 2 * n_qubits + 1
 
-    anti_comm_check = P_anticomm.adjacency_matrix.astype(int) - np.eye(P_anticomm.adjacency_matrix.shape[0])
-    assert(np.sum(anti_comm_check) == 0), 'operator needs to be made of anti-commuting Pauli operators'
+    ## expensive check
+    # anti_comm_check = P_anticomm.adjacency_matrix.astype(int) - np.eye(P_anticomm.adjacency_matrix.shape[0])
+    # assert(np.sum(anti_comm_check) == 0), 'operator needs to be made of anti-commuting Pauli operators'
 
     return P_anticomm
 
+
 def tensor_list(factor_list:List[PauliwordOp]) -> PauliwordOp:
     """ Given a list of PauliwordOps, recursively tensor from the right
     """
     return reduce(lambda x,y:x.tensor(y), factor_list)
 
 
-def gram_schmidt_from_quantum_state(state) ->np.array:
+def product_list(product_list:List[PauliwordOp]) -> PauliwordOp:
+    """ Given a list of PauliwordOps, recursively take product from the right
+    """
+    return reduce(lambda x,y:x*y, product_list)
+
+
+def gram_schmidt_from_quantum_state(state:Union[np.array, list, QuantumState]) ->np.array:
     """
     build a unitary to build a quantum state from the zero state (aka state defines first column of unitary)
     uses gram schmidt to find other (orthogonal) columns of matrix
@@ -102,14 +131,17 @@ def gram_schmidt_from_quantum_state(state) ->np.array:
     Returns:
         M (np.array): unitary matrix preparing input state from zero state
     """
-    state = np.asarray(state).reshape([-1])
-
-    N_qubits = round(np.log2(state.shape[0]))
 
-    missing_amps = 2**N_qubits - state.shape[0]
-    state = np.hstack((state, np.zeros(missing_amps, dtype=complex)))
+    if isinstance(state, QuantumState):
+        N_qubits = state.n_qubits
+        state = state.to_sparse_matrix.toarray().reshape([-1])
+    else:
+        state = np.asarray(state).reshape([-1])
+        N_qubits = round(np.log2(state.shape[0]))
+        missing_amps = 2**N_qubits - state.shape[0]
+        state = np.hstack((state, np.zeros(missing_amps, dtype=complex)))
 
-    assert len(state) == 2**N_qubits, 'state is not defined on power of two'
+    assert state.shape[0] == 2**N_qubits, 'state is not defined on power of two'
     assert np.isclose(np.linalg.norm(state), 1), 'state is not normalized'
 
     M = np.eye(2**N_qubits, dtype=complex)
@@ -130,6 +162,7 @@ def gram_schmidt_from_quantum_state(state) ->np.array:
 
     return M
 
+
 def Draw_molecule(
         xyz_string: str, width: int = 400, height: int = 400, style: str = "sphere"
     ) -> py3Dmol.view:

tests/test_operators/test_independent_op.py

@@ -42,10 +42,16 @@ def test_from_dictionary():
     )
     assert op1 == op2
 
-def test_no_symmetry_generator_error():
+# def test_no_symmetry_generator_error():
+#     op = PauliwordOp.from_list(['X', 'Y','Z'])
+#     with pytest.raises(RuntimeError):
+#         IndependentOp.symmetry_generators(op)
+
+def test_no_symmetry_generators():
     op = PauliwordOp.from_list(['X', 'Y','Z'])
-    with pytest.raises(RuntimeError):
-        IndependentOp.symmetry_generators(op)
+    ind_op = IndependentOp.symmetry_generators(op)
+    assert ind_op.n_terms == 0
+
 
 def test_commuting_overide_symmetry_generators():
     op = PauliwordOp.from_list(

tests/test_operators/test_noncontextual_op.py

@@ -190,6 +190,17 @@ def test_from_hamiltonian_basis():
     assert H_noncon_basis2.n_terms >= H_noncon_basis1.n_terms
 
 
+def test_noncon_no_symmertry_generators():
+    Pwords = PauliwordOp.from_list(['X', 'Y', 'Z'])
+    E_ground = -1.7320508075688772
+
+    with pytest.warns():
+        # capture warning of no Z2 symmetries
+        H_noncon = NoncontextualOp.from_PauliwordOp(Pwords)
+    assert H_noncon.n_terms == Pwords.n_terms
+    H_noncon.solve()
+    assert np.isclose(H_noncon.energy, E_ground)
+
 
 ####################################
 # Testing noncontextual optimizers #
@@ -348,10 +359,12 @@ def test_solve_brute_force_discrete_partial_ref():
 
     partial_reference_state = noncon_problem['partial_reference_state']
 
-    H_noncon_op.solve(strategy='brute_force',
-                   ref_state=partial_reference_state,
-                   discrete_optimization_order=None,
-                   num_anneals=None)
+    with pytest.warns():
+        # capture warning when Z stabilizers measured give zero expec value
+        H_noncon_op.solve(strategy='brute_force',
+                       ref_state=partial_reference_state,
+                       discrete_optimization_order=None,
+                       num_anneals=None)
     assert np.isclose(H_noncon_op.energy, noncon_problem['E'])
 
 
@@ -374,10 +387,12 @@ def test_solve_brute_force_PUSO_discrete_partial_ref():
 
     partial_reference_state = noncon_problem['partial_reference_state']
 
-    H_noncon_op.solve(strategy='brute_force_PUSO',
-                   ref_state=partial_reference_state,
-                   discrete_optimization_order='first', # <- HERE
-                   num_anneals=None)
+    with pytest.warns():
+        # capture warning when Z stabilizers measured give zero expec value
+        H_noncon_op.solve(strategy='brute_force_PUSO',
+                       ref_state=partial_reference_state,
+                       discrete_optimization_order='first', # <- HERE
+                       num_anneals=None)
     assert np.isclose(H_noncon_op.energy, noncon_problem['E'])
 
 
@@ -387,10 +402,12 @@ def test_solve_brute_force_QUSO_discrete_partial_ref():
 
     partial_reference_state = noncon_problem['partial_reference_state']
 
-    H_noncon_op.solve(strategy='brute_force_QUSO',
-                   ref_state=partial_reference_state,
-                   discrete_optimization_order='first', # <- HERE
-                   num_anneals=None)
+    with pytest.warns():
+        # capture warning when Z stabilizers measured give zero expec value
+        H_noncon_op.solve(strategy='brute_force_QUSO',
+                       ref_state=partial_reference_state,
+                       discrete_optimization_order='first', # <- HERE
+                       num_anneals=None)
     assert np.isclose(H_noncon_op.energy, noncon_problem['E'])
 
 
@@ -400,10 +417,12 @@ def test_solve_annealing_PUSO_discrete_partial_ref():
 
     partial_reference_state = noncon_problem['partial_reference_state']
 
-    H_noncon_op.solve(strategy='annealing_PUSO',
-                   ref_state=partial_reference_state,
-                   discrete_optimization_order='first', # <- HERE
-                   num_anneals=1_000)
+    with pytest.warns():
+        # capture warning when Z stabilizers measured give zero expec value
+        H_noncon_op.solve(strategy='annealing_PUSO',
+                       ref_state=partial_reference_state,
+                       discrete_optimization_order='first', # <- HERE
+                       num_anneals=1_000)
     assert np.isclose(H_noncon_op.energy, noncon_problem['E'])
 
 
@@ -413,10 +432,12 @@ def test_solve_annealing_QUSO_discrete_partial_ref():
 
     partial_reference_state = noncon_problem['partial_reference_state']
 
-    H_noncon_op.solve(strategy='annealing_QUSO',
-                   ref_state=partial_reference_state,
-                   discrete_optimization_order='first', # <- HERE
-                   num_anneals=1_000)
+    with pytest.warns():
+        # capture warning when Z stabilizers measured give zero expec value
+        H_noncon_op.solve(strategy='annealing_QUSO',
+                       ref_state=partial_reference_state,
+                       discrete_optimization_order='first', # <- HERE
+                       num_anneals=1_000)
     assert np.isclose(H_noncon_op.energy, noncon_problem['E'])
 
 
@@ -486,12 +507,17 @@ def test_get_qaoa_qubo_with_partial_reference():
     H_noncon_op = NoncontextualOp.from_PauliwordOp(H_noncon)
     reference = noncon_problem['partial_reference_state']
 
-    H_noncon_op.symmetry_generators.update_sector(reference)
+    with pytest.warns():
+        # capture warning when Z stabilizers measured give zero expec value
+        H_noncon_op.symmetry_generators.update_sector(reference)
+
     ev_assignment = H_noncon_op.symmetry_generators.coeff_vec
     fixed_ev_mask = ev_assignment != 0
     fixed_eigvals = (ev_assignment[fixed_ev_mask]).astype(int)
 
-    QAOA_dict = H_noncon_op.get_qaoa(ref_state=reference, type='qubo')
+    with pytest.warns():
+        # capture warning when Z stabilizers measured give zero expec value
+        QAOA_dict = H_noncon_op.get_qaoa(ref_state=reference, type='qubo')
 
     for key in QAOA_dict.keys():
         exp = QAOA_dict[key]
@@ -575,12 +601,17 @@ def test_get_qaoa_pubo_with_partial_reference():
     H_noncon_op = NoncontextualOp.from_PauliwordOp(H_noncon)
     reference = noncon_problem['partial_reference_state']
 
-    H_noncon_op.symmetry_generators.update_sector(reference)
+    with pytest.warns():
+        # capture warning when Z stabilizers measured give zero expec value
+        H_noncon_op.symmetry_generators.update_sector(reference)
+
     ev_assignment = H_noncon_op.symmetry_generators.coeff_vec
     fixed_ev_mask = ev_assignment != 0
     fixed_eigvals = (ev_assignment[fixed_ev_mask]).astype(int)
 
-    QAOA_dict = H_noncon_op.get_qaoa(ref_state=reference, type='pubo')
+    with pytest.warns():
+        # capture warning when Z stabilizers measured give zero expec value
+        QAOA_dict = H_noncon_op.get_qaoa(ref_state=reference, type='pubo')
 
     for key in QAOA_dict.keys():
         exp = QAOA_dict[key]