V1 upgrade

pyproject.toml

@@ -17,7 +17,6 @@ flake8 = "^4.0.1"
 py3Dmol = "^1.8.0"
 ncon = "^1.0.0"
 opt-einsum = "^3.3.0"
-qubovert = "^1.2.5"
 pytest = "^7.2.2"
 numba = "0.56"
 quimb = "^1.5.1"

symmer/operators/anticommuting_op.py

@@ -167,14 +167,31 @@ def unitary_partitioning(self, s_index: int=None, up_method: Optional[str]='seq_
             AC_op (AntiCommutingOp): normalized clique - i.e. self == gamma_l * AC_op
         """
         assert up_method in ['LCU', 'seq_rot'], f'unknown unitary partitioning method: {up_method}'
-        AC_op = self.copy()
+
+        if s_index is None:
+            s_index = self.get_least_dense_term_index()
+
+        if np.isclose(self.coeff_vec[s_index], 0):
+            # need to correct for s_index having zero coeff...
+            s_index = np.argmax(abs(self.coeff_vec))
+            warnings.warn(f's indexed term has zero coeff, s_index set to {s_index} so that nonzero operator is rotated onto')
+
+        s_index = int(s_index)
+        BsPs    = self[s_index]
+
+        # NOTE: term to reduce to is at the top of sym matrix i.e. s_index of ZERO now!
+        no_BsPs = (self - BsPs).cleanup()
+        if (len(no_BsPs.coeff_vec)==1 and no_BsPs.coeff_vec[0]==0):
+            AC_op = BsPs
+        else:
+            AC_op = BsPs.append(no_BsPs)
 
         if AC_op.n_terms == 1:
-            rotations = None
+            rotations = []
             gamma_l = np.linalg.norm(AC_op.coeff_vec)
             AC_op.coeff_vec = AC_op.coeff_vec / gamma_l
-            Ps = PauliwordOp(AC_op.symp_matrix, [1])
-            return Ps, rotations, gamma_l, AC_op
+            Ps = AC_op
+            return Ps, rotations, gamma_l, self.multiply_by_constant(1/gamma_l)
 
         else:
 
@@ -183,23 +200,6 @@ def unitary_partitioning(self, s_index: int=None, up_method: Optional[str]='seq_
             gamma_l = np.linalg.norm(AC_op.coeff_vec)
             AC_op.coeff_vec = AC_op.coeff_vec / gamma_l
 
-            if s_index is None:
-                s_index = self.get_least_dense_term_index()
-
-            if s_index!=0:
-                # re-order so s term is ALWAYS at top of symplectic matrix and thus is index as 0!
-                ### assert s_index <= AC_op.n_terms-1, 's_index out of range'
-                AC_op.coeff_vec[[0, s_index]] = AC_op.coeff_vec[[s_index, 0]]
-                AC_op.symp_matrix[[0, s_index]] = AC_op.symp_matrix[[s_index, 0]]
-                AC_op = AntiCommutingOp(AC_op.symp_matrix, AC_op.coeff_vec) # need to reinit otherwise Z and X blocks wrong
-
-            # assert not np.isclose(AC_op.coeff_vec[0], 0), f's_index cannot have zero coefficent: {AC_op.coeff_vec[0]}'
-            if np.isclose(AC_op[0].coeff_vec, 0):
-                # need to correct for s_index having zero coeff... then need to swap to nonzero index
-                non_zero_index = np.argmax(abs(AC_op.coeff_vec))
-                AC_op.coeff_vec[[0, non_zero_index]] = AC_op.coeff_vec[[non_zero_index, 0]]
-                AC_op.symp_matrix[[0, non_zero_index]] = AC_op.symp_matrix[[non_zero_index, 0]]
-
             if up_method=='seq_rot':
                 if len(self.X_sk_rotations)!=0:
                     self.X_sk_rotations = []
@@ -210,11 +210,31 @@ def unitary_partitioning(self, s_index: int=None, up_method: Optional[str]='seq_
                     self.R_LCU = None
 
                 Ps = self.generate_LCU_operator(AC_op)
-                rotations = self.R_LCU
+                rotations = LCU_as_seq_rot(self.R_LCU)
             else:
                 raise ValueError(f'unknown unitary partitioning method: {up_method}!')
 
-            return Ps, rotations, gamma_l, AC_op
+            return Ps, rotations, gamma_l, self.multiply_by_constant(1/gamma_l)
+
+    def multiply_by_constant(self, constant: float) -> "AntiCommutingOp":
+        """ Return AntiCommutingOp under constant multiplication
+        """
+        AC_op_copy = self.copy()
+        AC_op_copy.coeff_vec *= constant
+        return AC_op_copy
+
+    @classmethod
+    def random(cls, n_qubits: int, n_terms: Union[None, int]=None, apply_clifford=True) -> "AntiCommutingOp":
+        """
+        generate a random real coefficient anticommuting op
+
+        """
+        from symmer.utils import random_anitcomm_2n_1_PauliwordOp
+        if n_terms is None:
+            n_terms = 2*n_qubits+1
+
+        assert n_terms<= 2*n_qubits+1, f'cannot have {n_terms} Pops on {n_qubits} qubits'
+        return cls.from_PauliwordOp( random_anitcomm_2n_1_PauliwordOp(n_qubits, apply_clifford=apply_clifford)[:n_terms])
 
     def generate_LCU_operator(self, AC_op) -> PauliwordOp:
         """
@@ -235,74 +255,45 @@ def generate_LCU_operator(self, AC_op) -> PauliwordOp:
             R_LCU (PauliwordOp): PauliwordOp that is a linear combination of unitaries
             P_s (PauliwordOp): single PauliwordOp that has been reduced too.
         """
-        # need to remove zero coeff terms
-        AC_op_cpy = AC_op.copy()
-        before_cleanup = AC_op_cpy.n_terms
-        AC_op = AC_op_cpy[np.where(abs(AC_op.coeff_vec)>1e-15)[0]]
-        post_cleanup = AC_op.n_terms
-        # AC_op = AC_op.cleanup(zero_threshold=1e-15)  ## cleanup re-orders which is BAD for s_index
+        ## s_index is ensured to be in zero position in unitary_partitioning method! 
+        ## if using function without this method need to ensure term to rotate onto is the zeroth index of AC_op!
+        s_index=0
 
+        # note gamma_l norm applied on init!
+        Ps_LCU = PauliwordOp(AC_op.symp_matrix[s_index], [1])
+        βs = AC_op.coeff_vec[s_index]
 
-        if (before_cleanup>1 and post_cleanup==1):
-            if AC_op.coeff_vec[0]<0:
-                # need to fix neg sign (use Pauli multiplication)
+        #  ∑ β_k 𝑃_k  ... note this doesn't contain 𝛽_s 𝑃_s
+        no_βsPs = AC_op - (Ps_LCU.multiply_by_constant(βs))
 
-                # as s index defaults to 0, take the next term (in CS-VQE this will commute with symmetries)!
-                if np.isclose(AC_op_cpy[0].coeff_vec, 0):
-                    # need to correct for s_index having zero coeff... then need to swap to nonzero index
-                    non_zero_index = np.argmax(abs(AC_op_cpy.coeff_vec))
+        # Ω_𝑙 ∑ 𝛿_k 𝑃_k  ... renormalized!
+        omega_l = np.linalg.norm(no_βsPs.coeff_vec)
+        no_βsPs.coeff_vec = no_βsPs.coeff_vec / omega_l
 
-                    AC_op_cpy.coeff_vec[[0, non_zero_index]] = AC_op_cpy.coeff_vec[[non_zero_index, 0]]
-                    AC_op_cpy.symp_matrix[[0, non_zero_index]] = AC_op_cpy.symp_matrix[[non_zero_index, 0]]
+        phi_n_1 = np.arccos(βs)
+        # require sin(𝜙_{𝑛−1}) to be positive...
+        if (phi_n_1 > np.pi):
+            phi_n_1 = 2 * np.pi - phi_n_1
 
+        alpha = phi_n_1
+        I_term = 'I' * Ps_LCU.n_qubits
+        self.R_LCU = PauliwordOp.from_dictionary({I_term: np.cos(alpha / 2)})
 
-                sign_correction = PauliwordOp(AC_op_cpy.symp_matrix[1],[1])
+        sin_term = -np.sin(alpha / 2)
 
-                self.R_LCU = sign_correction
-                Ps_LCU = PauliwordOp(AC_op.symp_matrix, [1])
-            else:
-                self.R_LCU = PauliwordOp.from_list(['I'*AC_op.n_qubits])
-                Ps_LCU = PauliwordOp(AC_op.symp_matrix, AC_op.coeff_vec)
-        else:
-            s_index=0
-
-            # note gamma_l norm applied on init!
-            Ps_LCU = PauliwordOp(AC_op.symp_matrix[s_index], [1])
-            βs = AC_op.coeff_vec[s_index]
-
-            #  ∑ β_k 𝑃_k  ... note this doesn't contain 𝛽_s 𝑃_s
-            no_βsPs = AC_op - (Ps_LCU.multiply_by_constant(βs))
-
-            # Ω_𝑙 ∑ 𝛿_k 𝑃_k  ... renormalized!
-            omega_l = np.linalg.norm(no_βsPs.coeff_vec)
-            no_βsPs.coeff_vec = no_βsPs.coeff_vec / omega_l
-
-            phi_n_1 = np.arccos(βs)
-            # require sin(𝜙_{𝑛−1}) to be positive...
-            if (phi_n_1 > np.pi):
-                phi_n_1 = 2 * np.pi - phi_n_1
-
-            alpha = phi_n_1
-            I_term = 'I' * Ps_LCU.n_qubits
-            self.R_LCU = PauliwordOp.from_dictionary({I_term: np.cos(alpha / 2)})
-
-            sin_term = -np.sin(alpha / 2)
-
-            for dkPk in no_βsPs:
-                dk_PkPs = dkPk * Ps_LCU
-                self.R_LCU += dk_PkPs.multiply_by_constant(sin_term)
+        for dkPk in no_βsPs:
+            dk_PkPs = dkPk * Ps_LCU
+            self.R_LCU += dk_PkPs.multiply_by_constant(sin_term)
 
         return Ps_LCU
 
-def LCU_as_seq_rot(AC_op: PauliwordOp, include_global_phase_correction=False):
+def LCU_as_seq_rot(R_LCU: PauliwordOp):
     """
     Convert a unitary composed of a
     See equations 18 and 19 of https://arxiv.org/pdf/1907.09040.pdf
-    Note global phase fix is not necessary
 
     Args:
-        AC_op (PauliwordOp): unitary composed as a linear combination of anticommuting Pauli operators (excluding identity)
-        include_global_phase_correction (optional): whether to fix global phase to matrix input operator matrix exactly
+        R_LCU (PauliwordOp): unitary composed as a normalized linear combination of imaginary anticommuting Pauli operators (excluding identity)
     Returns:
         expon_p_terms (list): list of rotations generated by Pauli operators to implement AC_op unitary
 
@@ -326,40 +317,34 @@ def LCU_as_seq_rot(AC_op: PauliwordOp, include_global_phase_correction=False):
     check = reduce(lambda a,b: a*b, [exponentiate_single_Pop(x.multiply_by_constant(1j*y/2)) for x, y in exp_terms])
     print(check == rotations_LCU)
     """
-
-    assert AC_op.n_terms > 1, 'AC_op must have more than 1 term'
-    assert np.isclose(np.linalg.norm(AC_op.coeff_vec), 1), 'AC_op must be l2 normalized'
+    if isinstance(R_LCU, list) and len(R_LCU)==0:
+        # case where there are no rotations
+        return list()
+
+    assert R_LCU.n_terms > 1, 'AC_op must have more than 1 term'
+    assert np.isclose(np.linalg.norm(R_LCU.coeff_vec), 1), 'AC_op must be l2 normalized'
 
     expon_p_terms = []
 
     # IF imaginary components the this makes real (but need phase correction later!)
-    coeff_vec = AC_op.coeff_vec.real + AC_op.coeff_vec.imag
+    coeff_vec = R_LCU.coeff_vec.real + R_LCU.coeff_vec.imag
     for k, c_k in enumerate(coeff_vec):
-        P_k = AC_op[k]
+        P_k = R_LCU[k]
         theta_k = np.arcsin(c_k / np.linalg.norm(coeff_vec[:(k + 1)]))
         P_k.coeff_vec[0] = 1
         expon_p_terms.append(tuple((P_k, theta_k)))
 
     expon_p_terms = [*expon_p_terms, *expon_p_terms[::-1]]
 
-    ### check
-    # from symmer.evolution.exponentiation import exponentiate_single_Pop
-    # terms = [exponentiate_single_Pop(op.multiply_by_constant(1j*angle/2)) for op,angle in expon_p_terms]
-    # final_op = reduce(lambda x,y: x*y, terms) * PauliwordOp.from_dictionary({'I'*AC_op.n_qubits: -1j})
-    # assert AC_op == final_op
-
-    # in circuit this would be done with Z * Y * X gate series:
-    # global phase correction (not necessary)
-    ## phase_correction = PauliwordOp.from_dictionary({'I'*AC_op.n_qubits: -1j})
-
-    if include_global_phase_correction:
-        ## multiply by -1j Identity term!
-        phase_rot = (PauliwordOp.from_dictionary({'I' * AC_op.n_qubits: 1}), -np.pi)
-        expon_p_terms.append(phase_rot)
-
-        # check1 = reduce(lambda a,b: a*b, [exponentiate_single_Pop(x.multiply_by_constant(1j*y/2)) for x, y in expon_p_terms])
-        # assert check1 == AC_op
-
+    ##### manual phase correction with a rotation!
+    # if include_global_phase_correction:
+    #     ## multiply by -1j Identity term!
+    #     phase_rot = (PauliwordOp.from_dictionary({'I' * R_LCU.n_qubits: 1}), -np.pi)
+    #     expon_p_terms.append(phase_rot)
+
+    ## phase correction - change angle by -pi in first rotation!
+    expon_p_terms[0] = (expon_p_terms[0][0], expon_p_terms[0][1]-np.pi)
+
     return expon_p_terms
 
 # from symmer.operators.utils import mul_symplectic

symmer/operators/base.py

@@ -945,7 +945,7 @@ def __getitem__(self,
         elif isinstance(key, (list, np.ndarray)):
             mask = np.asarray(key)
         else:
-            raise ValueError('Unrecognised input, must be an integer, slice, list or np.array')
+            raise ValueError(f'Unrecognised input {type(key)}, must be an integer, slice, list or np.array')
 
         symp_items = self.symp_matrix[mask]
         coeff_items = self.coeff_vec[mask]
@@ -1098,7 +1098,8 @@ def adjacency_matrix_qwc(self) -> np.array:
     def is_noncontextual(self) -> bool:
         """ 
         Returns True if the operator is noncontextual, False if contextual
-        Scales as O(N^2), compared with the O(N^3) algorithm of https://doi.org/10.1103/PhysRevLett.123.200501
+        Scales as O(M^2), compared with the O(M^3) algorithm of https://doi.org/10.1103/PhysRevLett.123.200501
+        where M is the number of terms in the operator.
         Constructing the adjacency matrix is by far the most expensive part - very fast once that has been built.
 
         Returns:
@@ -1111,7 +1112,8 @@ def is_noncontextual(self) -> bool:
 
     def _rotate_by_single_Pword(self,
             Pword: "PauliwordOp", 
-            angle: float = None
+            angle: float = None,
+            threshold: float = 1e-18
         ) -> "PauliwordOp":
         """ 
         Let R(t) = e^{i t/2 Q} = cos(t/2)*I + i*sin(t/2)*Q, then one of the following can occur:
@@ -1127,10 +1129,17 @@ def _rotate_by_single_Pword(self,
         Args:
             Pword (PauliwordOp): The Pauliword to rotate by.
             angle (float): The rotation angle in radians. If None, a Clifford rotation (angle=pi/2) is assumed.
-
+            threshold (float): Angle threshold for Clifford rotation (precision to which the angle is a multiple of pi/2)
         Returns:
             PauliwordOp: The rotated operator.
         """
+        if angle is None: # for legacy compatibility
+            angle = np.pi/2
+
+        if angle.imag != 0:
+            warnings.warn('Complex component in angle: this will be ignored.')
+        angle = angle.real
+
         assert(Pword.n_terms==1), 'Only rotation by single Pauliword allowed here'
         if Pword.coeff_vec[0] != 1:
             # non-1 coefficients will affect the sign and angle in the exponent of R(t)
@@ -1150,15 +1159,25 @@ def _rotate_by_single_Pword(self,
             commute_self = PauliwordOp(self.symp_matrix[commute_vec], self.coeff_vec[commute_vec])
             anticom_self = PauliwordOp(self.symp_matrix[~commute_vec], self.coeff_vec[~commute_vec])
 
-            if angle is None:
-                # assumes pi/2 rotation so Clifford
-                anticom_part = (anticom_self*Pword_copy).multiply_by_constant(-1j)
+            multiple = angle * 2 / np.pi
+            int_part = round(multiple)
+            if abs(int_part - multiple)<=threshold:
+                if int_part % 2 == 0:
+                    # no rotation for angle congruent to 0 or pi disregarding sign, fixed in next line
+                    anticom_part = anticom_self
+                else:
+                    # rotation of -pi/2 disregarding sign, fixed in next line
+                    anticom_part = (anticom_self*Pword_copy).multiply_by_constant(-1j)
+                if int_part in [2,3]:
+                    anticom_part = anticom_part.multiply_by_constant(-1)
                 # if rotation is Clifford cannot produce duplicate terms so cleanup not necessary
                 return PauliwordOp(
                     np.vstack([anticom_part.symp_matrix, commute_self.symp_matrix]), 
                     np.hstack([anticom_part.coeff_vec, commute_self.coeff_vec])
                 )
             else:
+                if abs(angle)>1e6:
+                    warnings.warn('Large angle can lead to precision errors: recommend using high-precision math library such as mpmath or redefine angle in range [-pi, pi]')
                 # if angle is specified, performs non-Clifford rotation
                 anticom_part = (anticom_self.multiply_by_constant(np.cos(angle)) + 
                                 (anticom_self*Pword_copy).multiply_by_constant(-1j*np.sin(angle)))
@@ -1182,9 +1201,12 @@ def perform_rotations(self,
             PauliwordOp: The operator after performing the rotations.
         """
         op_copy = self.copy()
-        for pauli_rotation, angle in rotations:
-            op_copy = op_copy._rotate_by_single_Pword(pauli_rotation, angle).cleanup()
-        return op_copy
+        if rotations == []:
+            return op_copy.cleanup()
+        else:
+            for pauli_rotation, angle in rotations:
+                op_copy = op_copy._rotate_by_single_Pword(pauli_rotation, angle).cleanup()
+            return op_copy
 
     def tensor(self, right_op: "PauliwordOp") -> "PauliwordOp":
         """