Merge pull request #189 from UCL-CCS/updated_to_sparse

Updated to sparse

symmer/__init__.py

@@ -1,4 +1,6 @@
 """Main init for package."""
+import warnings
+warnings.filterwarnings('ignore', module='cotengra')
 from symmer.process_handler import process
 from symmer.operators  import PauliwordOp, QuantumState
 from symmer.projection import QubitTapering, ContextualSubspace, QubitSubspaceManager

symmer/operators/base.py

@@ -21,6 +21,11 @@
 from qiskit.quantum_info import SparsePauliOp
 warnings.simplefilter('always', UserWarning)
 
+from qiskit._accelerate.sparse_pauli_op import (
+    ZXPaulis,
+    to_matrix_sparse,
+)
+
 from numba.core.errors import NumbaDeprecationWarning, NumbaPendingDeprecationWarning
 warnings.simplefilter('ignore', category=NumbaDeprecationWarning)
 warnings.simplefilter('ignore', category=NumbaPendingDeprecationWarning)
@@ -1464,33 +1469,45 @@ def to_sparse_matrix(self) -> csr_matrix:
         if self.n_qubits == 0:
             return csr_matrix(self.coeff_vec)
 
-        if self.n_qubits>15:
-            from symmer.utils import get_sparse_matrix_large_pauliwordop
-            sparse_matrix = get_sparse_matrix_large_pauliwordop(self)
-            return sparse_matrix
-        else:
-            x_int = binary_array_to_int(self.X_block).reshape(-1, 1)
-            z_int = binary_array_to_int(self.Z_block).reshape(-1, 1)
+        # if self.n_qubits>15:
+        #     from symmer.utils import get_sparse_matrix_large_pauliwordop
+        #     sparse_matrix = get_sparse_matrix_large_pauliwordop(self)
+        #     return sparse_matrix
+        # else:
+        #     x_int = binary_array_to_int(self.X_block).reshape(-1, 1)
+        #     z_int = binary_array_to_int(self.Z_block).reshape(-1, 1)
 
-            Y_number = np.sum(np.bitwise_and(self.X_block, self.Z_block), axis=1)
-            global_phase = (-1j) ** Y_number
+        #     Y_number = np.sum(np.bitwise_and(self.X_block, self.Z_block), axis=1)
+        #     global_phase = (-1j) ** Y_number
 
-            dimension = 2 ** self.n_qubits
-            row_ind = np.repeat(np.arange(dimension).reshape(1, -1), self.X_block.shape[0], axis=0)
-            col_ind = np.bitwise_xor(row_ind, x_int)
+        #     dimension = 2 ** self.n_qubits
+        #     row_ind = np.repeat(np.arange(dimension).reshape(1, -1), self.X_block.shape[0], axis=0)
+        #     col_ind = np.bitwise_xor(row_ind, x_int)
 
-            row_inds_and_Zint = np.bitwise_and(row_ind, z_int)
-            vals = global_phase.reshape(-1, 1) * (-1) ** (
-                        count1_in_int_bitstring(row_inds_and_Zint) % 2)  # .astype(complex))
+        #     row_inds_and_Zint = np.bitwise_and(row_ind, z_int)
+        #     vals = global_phase.reshape(-1, 1) * (-1) ** (
+        #                 count1_in_int_bitstring(row_inds_and_Zint) % 2)  # .astype(complex))
 
-            values_and_coeff = np.einsum('ij,i->ij', vals, self.coeff_vec)
+        #     values_and_coeff = np.einsum('ij,i->ij', vals, self.coeff_vec)
 
-            sparse_matrix = csr_matrix(
-                (values_and_coeff.flatten(), (row_ind.flatten(), col_ind.flatten())),
-                shape=(dimension, dimension),
-                dtype=complex
-            )
-            return sparse_matrix
+        #     sparse_matrix = csr_matrix(
+        #         (values_and_coeff.flatten(), (row_ind.flatten(), col_ind.flatten())),
+        #         shape=(dimension, dimension),
+        #         dtype=complex
+        #     )
+        #     return sparse_matrix
+
+        phase = np.zeros(self.n_terms, dtype=np.uint8)
+        zx = ZXPaulis(
+                    self.X_block[:,::-1],
+                    self.Z_block[:,::-1],
+                    phase,
+                    self.coeff_vec,
+                )
+
+        data, indices, indptr = to_matrix_sparse(zx, force_serial=False)
+        side = 1 << self.n_qubits
+        return csr_matrix((data, indices, indptr), shape=(side, side))
 
     def conjugate_op(self, R: 'PauliwordOp') -> 'PauliwordOp':
         """

symmer/process_handler.py

@@ -5,14 +5,18 @@
 from ray import remote, put, get
 from multiprocessing import Process, Queue, set_start_method
 
-if sys.platform.lower() in ['linux', 'darwin']:
+if sys.platform.lower() in ['linux', 'darwin', 'linux2']:
     set_start_method('fork', force = True)
 else:
     set_start_method('spawn', force = True)
 
 class ProcessHandler:
 
-    method  = 'ray'
+    if sys.platform.lower() in ['linux', 'darwin', 'linux2']:
+        method  = 'mp'
+    else:
+        method  = 'ray'
+
     verbose = False
 
     def __init__(self):