Merge pull request #152 from UCL-CCS/bug_fixes

fixed error in noncontextual op (noncommuting Z2 symmetries)

symmer/operators/base.py

@@ -1,3 +1,7 @@
+import os
+import quimb
+from symmer.operators.utils import *
+import ray
 import numpy as np
 import pandas as pd
 import networkx as nx
@@ -9,7 +13,7 @@
 import multiprocessing as mp
 from cached_property import cached_property
 from scipy.sparse import csr_matrix, csc_matrix, coo_matrix, dok_matrix
-from symmer.operators.utils import *
+
 from openfermion import QubitOperator, count_qubits
 import matplotlib.pyplot as plt
 from qiskit.opflow import PauliSumOp as ibm_PauliSumOp
@@ -21,8 +25,6 @@
 warnings.simplefilter('ignore', category=NumbaDeprecationWarning)
 warnings.simplefilter('ignore', category=NumbaPendingDeprecationWarning)
 
-import ray
-
 class PauliwordOp:
  """ 
  A class thats represents an operator defined over the Pauli group in the symplectic representation.
@@ -557,7 +559,6 @@ def jordan_generator_reconstruction(self, generators: "PauliwordOp"):
  """
  assert check_jordan_independent(generators), 'The non-symmetry elements do not pairwise anticommute.'
 
-
  # first, separate symmetry elements from anticommuting ones
  symmetry_mask = np.all(generators.commutes_termwise(generators), axis=1)
 
@@ -569,23 +570,22 @@ def jordan_generator_reconstruction(self, generators: "PauliwordOp"):
  op_reconstruction = np.zeros([self.n_terms, generators.n_terms])
  successfully_reconstructed = np.zeros(self.n_terms, dtype=bool)
 
- Z2_terms = generators[symmetry_mask]
  ac_terms = generators[~symmetry_mask]
 
  # loop over anticommuting elements to enforce Jordan condition (no two anticommuting elements multiplied)
  for _, clq in ac_terms.clique_cover(edge_relation='C').items():
- clq_indices = [np.where(np.all(generators.symp_matrix == t, axis=1))[0][0] for t in clq.symp_matrix] 
+ clq_indices = [np.where(np.all(generators.symp_matrix == t, axis=1))[0][0] for t in clq.symp_matrix]
  mask_symmetries_with_P = symmetry_mask.copy()
- mask_symmetries_with_P[np.array(clq_indices)] = True 
+ mask_symmetries_with_P[np.array(clq_indices)] = True
  # reconstruct this PauliwordOp in the augemented symmetry + single anticommuting term generating set
  augmented_symmetries = generators[mask_symmetries_with_P]
  recon_mat_P, successful_P = self.generator_reconstruction(augmented_symmetries)
  # np.ix_ needed to correctly slice op_reconstruction as mask method does not work
- row, col = np.ix_(successful_P,mask_symmetries_with_P)
+ row, col = np.ix_(successful_P, mask_symmetries_with_P)
  op_reconstruction[row, col] = recon_mat_P[successful_P]
  # will have duplicate succesful reconstruction of symmetries, so only sets True once in logical OR
  successfully_reconstructed = np.logical_or(successfully_reconstructed, successful_P)
- 
+
  return op_reconstruction.astype(int), successfully_reconstructed
 
  @cached_property
@@ -1089,6 +1089,8 @@ def is_noncontextual(self) -> bool:
  bool: True if the operator is noncontextual, False if contextual.
  """
  return check_adjmat_noncontextual(self.generators.adjacency_matrix)
+ # from symmer.utils import get_generators_including_xz_products
+ # return check_adjmat_noncontextual(get_generators_including_xz_products(self).adjacency_matrix)
 
  def _rotate_by_single_Pword(self, 
  Pword: "PauliwordOp", 
@@ -1415,7 +1417,7 @@ def generators(self) -> "PauliwordOp":
  row_red = _rref_binary(self.symp_matrix)
  non_zero_rows = row_red[np.sum(row_red, axis=1).astype(bool)]
  generators = PauliwordOp(non_zero_rows,
- np.ones(non_zero_rows.shape[0]))
+ np.ones(non_zero_rows.shape[0], dtype=complex))
 
  assert check_independent(generators), 'generators are not independent'
  assert generators.n_terms <= 2*self.n_qubits, 'cannot have an independent generating set of size greaterthan 2 time num qubits'
@@ -2350,7 +2352,7 @@ def get_ij_operator(i:int, j:int, n_qubits:int,
  else:
  return ij_symp_matrix, coeffs
 
-@ray.remote
+@ray.remote(num_cpus=os.cpu_count())
 def single_term_expval(P_op: PauliwordOp, psi: QuantumState) -> float:
  """ 
  Expectation value calculation for a single Pauli operator given a QuantumState psi

symmer/operators/independent_op.py

@@ -29,7 +29,7 @@ def __init__(self,
  target_sqp (str): Target SQP (Single-Qubit Pauli). By default it is to 'Z'.
  """
  if coeff_vec is None:
- coeff_vec = np.ones(symp_matrix.shape[0])
+ coeff_vec = np.ones(symp_matrix.shape[0], dtype=complex)
  super().__init__(symp_matrix, coeff_vec)
  self._check_stab()
  self.coeff_vec = self.coeff_vec.real.astype(int)
@@ -141,7 +141,7 @@ def symmetry_generators(cls,
  S_commuting = S.clique_cover(edge_relation='C')[0]
  warnings.warn('Greedy method may identify non-optimal commuting symmetry terms; might be able to taper again.')
 
- return cls(S_commuting.symp_matrix, np.ones(S_commuting.n_terms))
+ return cls(S_commuting.symp_matrix, np.ones(S_commuting.n_terms, dtype=complex))
 
  def _check_stab(self) -> None:
  """ 
@@ -169,7 +169,7 @@ def __str__(self) -> str:
  out_string = ''
  for pauli_vec, coeff in zip(self.symp_matrix, self.coeff_vec):
  p_string = symplectic_to_string(pauli_vec)
- out_string += (f'{coeff: d} {p_string} \n')
+ out_string += (f'{coeff} {p_string} \n')
  return out_string[:-2]
 
  def __repr__(self) -> str:

symmer/operators/noncontextual_op.py

@@ -1,4 +1,5 @@
 import warnings
+import os
 import itertools
 import numpy as np
 import networkx as nx
@@ -243,12 +244,12 @@ def _from_generators_noncontextual_op(cls,
  NoncontextualOp: A NoncontextualOp instance constructed from the given Hamiltonian and generators.
  """
  assert generators is not None, 'Must specify a noncontextual generating set.'
+ assert generators.is_noncontextual, 'Generating set is contextual.'
  if use_jordan_product:
  _, noncontextual_terms_mask = H.jordan_generator_reconstruction(generators)
  else:
- assert generators.is_noncontextual, 'Generating set is contextual.'
  _, noncontextual_terms_mask = H.generator_reconstruction(generators, override_independence_check=True)
- 
+
  return cls.from_PauliwordOp(H[noncontextual_terms_mask])
 
  @classmethod
@@ -266,7 +267,12 @@ def _from_stabilizers_noncontextual_op(cls,
  NoncontextualOp: A NoncontextualOp instance constructed from the given PauliwordOp and stabilizers.
  """
  symmetries = IndependentOp.symmetry_generators(stabilizers, commuting_override=True)
- generators = NoncontextualOp.from_hamiltonian(symmetries, strategy='DFS_magnitude')
+ noncon = NoncontextualOp.from_hamiltonian(symmetries, strategy='DFS_magnitude')
+ generators = noncon.symmetry_generators
+ if noncon.clique_operator.n_terms>0:
+ generators+=noncon.clique_operator
+ use_jordan_product=True
+
  return cls._from_generators_noncontextual_op(H=H, generators=generators, use_jordan_product=use_jordan_product)
 
  def draw_graph_structure(self, 
@@ -291,7 +297,6 @@ def draw_graph_structure(self,
  include_symmetries (bool, optional): Determines whether to include symmetry edges in the visualization. Default is True.
  """
  adjmat = self.adjacency_matrix.copy()
- adjmat = self.adjacency_matrix.copy()
  index_symmetries = np.where(np.all(adjmat, axis=1))[0]
  np.fill_diagonal(adjmat, False)
 
@@ -314,17 +319,31 @@ def noncontextual_generators(self) -> None:
  """ 
  Find an independent generating set for the noncontextual operator.
  """
+ Z2_symmerties = IndependentOp.symmetry_generators(self, commuting_override=True)
 
- ## rather than searching over self (operator iteself) use generators instead ==> much faster
- # get Z2 symmetry of generators (may anticommute among themselves, but will commute with all generators)
- Z2_symmerties = IndependentOp.symmetry_generators(self.generators, commuting_override=True)
+ if not np.all(Z2_symmerties.commutes_termwise(Z2_symmerties)):
+ # need to account for Z2_symmerties not commuting with themselves
+ sym_gens = self.generators
+ z2_mask = np.sum(sym_gens.adjacency_matrix, axis=1) == sym_gens.n_terms
 
- # find terms not generated these Z2 symmerties
- # Note: CANNOT just use commuting subset from generators as this misses parts of problem off
- _, z2_mask = self.generator_reconstruction(Z2_symmerties)
+ Z2_incomplete = sym_gens[z2_mask]
+ _, missing_mask = sym_gens.generator_reconstruction(Z2_incomplete)
+ Z2_missing = sym_gens[~missing_mask]
+
+ cover = Z2_missing.clique_cover('C')
+ clique_rep_list = [C.sort()[0] for C in cover.values()]
+
+ sym_from_cliques = sum((cover[n] - C_rep) * C_rep for n, C_rep in enumerate(clique_rep_list) if
+ cover[n].n_terms > 1)
+
+ Z2_symmerties = (sym_from_cliques + Z2_incomplete).generators
+ _, z2_mask = self.generator_reconstruction(Z2_symmerties)
+ else:
+ _, z2_mask = self.generator_reconstruction(Z2_symmerties)
+
+ remaining = self[~z2_mask]
+ self.decomposed = remaining.clique_cover('C')
 
- ## noncon structure means remaining terms should be disjoint commuting cliques (graph colouring is trivial here!)
- self.decomposed = self[~z2_mask].clique_cover('C')
  self.n_cliques = len(self.decomposed)
  if self.n_cliques > 0:
  # choose clique representatives with the greatest coefficient
@@ -333,24 +352,20 @@ def noncontextual_generators(self) -> None:
  self.clique_operator = AntiCommutingOp.from_PauliwordOp(
  sum(clique_rep_list)
  )
+ self.clique_operator.coeff_vec = np.ones_like(self.clique_operator.coeff_vec)
 
  ## cliques can form new Z2 syms
- # sym_from_cliques = sum((self.decomposed[n]-C_rep) * C_rep for n, C_rep in enumerate(clique_rep_list) if
- # self.decomposed[n].n_terms > 1)
- # if sym_from_cliques:
- # _, mask = sym_from_cliques.generator_reconstruction(Z2_symmerties)
- # sym_from_cliques.coeff_vec = np.ones_like(sym_from_cliques.coeff_vec)
- # Z2_symmerties += sym_from_cliques[~mask]
-
+ sym_from_cliques = sum((self.decomposed[n] - C_rep) * C_rep for n, C_rep in enumerate(clique_rep_list) if
+ self.decomposed[n].n_terms > 1)
+ if sym_from_cliques:
+ Z2_symmerties = (sym_from_cliques + Z2_symmerties).generators
  else:
  self.clique_operator = PauliwordOp.empty(self.n_qubits).cleanup()
 
  self.symmetry_generators = IndependentOp.from_PauliwordOp(Z2_symmerties)
  _, Z2_mask = self.generator_reconstruction(Z2_symmerties)
  self.decomposed['symmetry'] = self[Z2_mask]
 
-
-
  def noncontextual_reconstruction(self) -> None:
  """ 
  Reconstruct the noncontextual operator in each independent basis GuCi - one for every clique.
@@ -360,6 +375,7 @@ def noncontextual_reconstruction(self) -> None:
  np.vstack([self.symmetry_generators.symp_matrix, self.clique_operator.symp_matrix]),
  np.ones(self.symmetry_generators.n_terms + self.n_cliques)
  )
+
  # Cannot simultaneously know eigenvalues of cliques so we peform a generator reconstruction
  # that respects the jordan product A*B = {A, B}/2, i.e. anticommuting elements are zeroed out
  jordan_recon_matrix, successful = self.jordan_generator_reconstruction(noncon_generators)#, override_independence_check=True)
@@ -749,12 +765,12 @@ def energy_xUSO(self) -> Tuple[float, np.array, np.array]:
  nu_vec = np.ones(self.NC_op.symmetry_generators.n_terms, dtype=int)
  nu_vec[self.fixed_ev_mask] = self.fixed_eigvals
  # must ensure the binary variables are correctly ordered in the solution:
- nu_vec[~self.fixed_ev_mask] = np.array([solution[f'x{i}'] for i in range(COST.num_binary_variables)])
+ nu_vec[~self.fixed_ev_mask] = np.array([solution[x_i] for x_i in sorted(COST.variables)])
 
  return self.NC_op.get_energy(nu_vec), nu_vec
 
 
-@ray.remote
+@ray.remote(num_cpus=os.cpu_count())
 def get_noncon_energy(noncon_H:NoncontextualOp, nu: np.array) -> float:
  """
  The classical objective function that encodes the noncontextual energies.

symmer/operators/utils.py

@@ -431,17 +431,22 @@ def check_jordan_independent(operators):
  # check_independent(Symmetries) & 
  # np.all(Anticommuting.adjacency_matrix == np.eye(Anticommuting.n_terms))
  # )
+ from symmer.operators import IndependentOp
 
- symmetry_mask = np.all(operators.commutes_termwise(operators), axis=1)
- Z2_terms = operators[symmetry_mask]
- ac_terms = operators[~symmetry_mask]
+ Z2_gen_mask = np.sum(operators.commutes_termwise(operators), axis=1) == operators.n_terms
+ Z2_terms = operators[Z2_gen_mask]
 
+ # find terms not generated these Z2 symmerties
+ _, z2_mask = operators.generator_reconstruction(Z2_terms)
+
+ ac_terms = operators[~z2_mask]
  if ac_terms.n_terms > 0:
  decomposed = ac_terms.clique_cover(edge_relation='C')
  clique_rep_list = [C.sort()[0] for C in decomposed.values()]
- ac_op = sum(clique_rep_list)
- assert (np.sum(ac_op.adjacency_matrix.astype(int)
- - np.eye(ac_op.adjacency_matrix.shape[0])) == 0), 'ac_op is not anticommuting'
+ ac_op = sum(clique_rep_list).cleanup()
+ if ac_op.n_terms>0:
+ assert (np.sum(ac_op.adjacency_matrix.astype(int)
+ - np.eye(ac_op.adjacency_matrix.shape[0])) == 0), f'ac_op is not anticommuting: {ac_op}'
  del ac_op
 
  Z2_from_cliques = sum((decomposed[n] - C_rep) * C_rep for n, C_rep in enumerate(clique_rep_list) if
@@ -517,6 +522,10 @@ def perform_noncontextual_sweep(operator) -> "PauliwordOp":
  np.ones(non_zero_rows.shape[0]))
  if not check_adjmat_noncontextual(generators.adjacency_matrix):
  mask[ind] = ~mask[ind]
+
+ # noncon_H = operator[mask]
+ # assert noncon_H.is_noncontextual, 'DFS output is not noncontextual'
+ # return noncon_H
  return operator[mask]
 
 def binary_array_to_int(bin_arr):

symmer/projection/contextual_subspace.py

@@ -1,7 +1,8 @@
 import numpy as np
 from symmer.operators import PauliwordOp, IndependentOp, NoncontextualOp, QuantumState
 from symmer.projection.utils import (
- update_eigenvalues, StabilizerIdentification, ObservableBiasing, stabilizer_walk
+ update_eigenvalues, StabilizerIdentification, ObservableBiasing, stabilizer_walk,
+ # get_noncon_generators_from_commuting_stabilizers
 )
 from symmer.projection import S3_projection
 from symmer.evolution import trotter

symmer/projection/utils.py

@@ -2,6 +2,8 @@
 from copy import deepcopy
 from scipy.optimize import differential_evolution
 from symmer.operators import PauliwordOp, IndependentOp
+from typing import Union, Optional
+from symmer.utils import random_anitcomm_2n_1_PauliwordOp, product_list
 
 def norm(vector: np.array) -> float:
  """
@@ -268,4 +270,70 @@ def objective(x):
  if print_info:
  print(f'Optimal score w(S)={stab_score} for HOMO/LUMO bias {bias_param}')
 
- return S
+ return S
+
+def get_noncon_generators_from_commuting_stabilizers(stabilizers: Union[PauliwordOp, IndependentOp],
+ weighting_operator: PauliwordOp,
+ return_clique_only: Optional[bool] = False) -> IndependentOp:
+ """
+ Given a set of commuting stabilizers and weighting operator find best noncontextual generating set
+ (ie works out best anticommuting addition to generators that reconstructs most of the weighting_operator)
+
+ Args:
+ stabilizers (PauliwordOp): operator containing commuting symmetries
+ weighting_operator (PauliwordOp): operator to inform ac choice
+
+ Returns:
+ new_stabilizers (IndependentOp): noncontextual generators that contain ac component (generates noncon op under Jordan product)
+ """
+ if not np.all(stabilizers.commutes_termwise(stabilizers)):
+ # stabilizers already contain ac component
+ return stabilizers #, PauliwordOp.empty(stabilizers.n_qubits).cleanup()
+ else:
+ # below generates the generators of inout stabilizers
+ generators = stabilizers.generators
+
+ best_l1_norm = -1
+ # find qubits uniquely defined by generators
+ unique_q_inds = ~(np.sum(np.logical_xor(generators.Z_block, generators.X_block), axis=0)-1).astype(bool)
+ for stab in generators:
+ # find unique non identity positions
+ act_positions = np.logical_and(np.logical_xor(stab.Z_block, stab.X_block)[0], unique_q_inds)
+
+ # work out number of qubits on these positions
+ n_act_qubits = np.sum(act_positions)
+
+ # find AC clique of size 2n containing given stabilizer
+ ac_basis = random_anitcomm_2n_1_PauliwordOp(n_act_qubits, apply_clifford=False)[1:]
+ new_basis = PauliwordOp(np.zeros((n_act_qubits * 2, stab.n_qubits * 2), dtype=bool),
+ np.ones(n_act_qubits * 2))
+
+ new_basis.symp_matrix[:, [*act_positions, *act_positions]] = ac_basis.symp_matrix
+
+ # ensure stab is in new_basis
+ gen, mask = stab.generator_reconstruction(new_basis)
+ required_products = gen[0].nonzero()[0][1:]
+ if len(required_products) > 0:
+ prod = product_list(new_basis[required_products])
+ new_basis = (new_basis * prod).cleanup()
+ new_basis.coeff_vec = np.ones_like(new_basis.coeff_vec)
+
+ # find best reconstruction
+ _, mask = weighting_operator.generator_reconstruction(new_basis)
+ success = weighting_operator[mask]
+ l1_norm = np.linalg.norm(success.coeff_vec, ord=1)
+
+ if l1_norm > best_l1_norm:
+ new_stabilizers = generators - stab + new_basis
+ best_l1_norm = l1_norm
+ stab_used = stab.copy()
+
+ assert new_stabilizers.is_noncontextual, 'new stabilizers are not noncontextual'
+
+ # commuting_stabs = IndependentOp.from_PauliwordOp(stabilizers)
+ # anticommuting_stabs = IndependentOp.from_PauliwordOp(new_stabilizers) - commuting_stabs
+ # return commuting_stabs, anticommuting_stabs
+ if return_clique_only:
+ return IndependentOp.from_PauliwordOp(new_stabilizers) - generators, stab_used
+ else:
+ return IndependentOp.from_PauliwordOp(new_stabilizers)