From f3b457a12c4dad3c03cd6e0a6493555a0dde6678 Mon Sep 17 00:00:00 2001
From: Rene Zander <rene.zander@fokus.fraunhofer.de>
Date: Fri, 15 Nov 2024 16:08:16 +0100
Subject: [PATCH] Implementation of change_of_basis for commuting groups of
 QubitOperators

---
 src/qrisp/operators/qubit/__init__.py         |   1 +
 .../operators/qubit/commutativity_tools.py    | 256 ++++++++++++++++++
 src/qrisp/operators/qubit/qubit_operator.py   | 136 ++++++++--
 src/qrisp/operators/qubit/qubit_term.py       |  34 ++-
 4 files changed, 395 insertions(+), 32 deletions(-)
 create mode 100644 src/qrisp/operators/qubit/commutativity_tools.py

diff --git a/src/qrisp/operators/qubit/__init__.py b/src/qrisp/operators/qubit/__init__.py
index df4f7e94..b84b0e41 100644
--- a/src/qrisp/operators/qubit/__init__.py
+++ b/src/qrisp/operators/qubit/__init__.py
@@ -22,3 +22,4 @@
 from qrisp.operators.qubit.bound_qubit_operator import *
 #from qrisp.operators.qubit.pauli_measurement import *
 from qrisp.operators.qubit.operator_factors import *
+from qrisp.operators.qubit.commutativity_tools import *
diff --git a/src/qrisp/operators/qubit/commutativity_tools.py b/src/qrisp/operators/qubit/commutativity_tools.py
new file mode 100644
index 00000000..7e8cade3
--- /dev/null
+++ b/src/qrisp/operators/qubit/commutativity_tools.py
@@ -0,0 +1,256 @@
+"""
+\********************************************************************************
+* Copyright (c) 2024 the Qrisp authors
+*
+* This program and the accompanying materials are made available under the
+* terms of the Eclipse Public License 2.0 which is available at
+* http://www.eclipse.org/legal/epl-2.0.
+*
+* This Source Code may also be made available under the following Secondary
+* Licenses when the conditions for such availability set forth in the Eclipse
+* Public License, v. 2.0 are satisfied: GNU General Public License, version 2
+* with the GNU Classpath Exception which is
+* available at https://www.gnu.org/software/classpath/license.html.
+*
+* SPDX-License-Identifier: EPL-2.0 OR GPL-2.0 WITH Classpath-exception-2.0
+********************************************************************************/
+"""
+
+import numpy as np
+
+
+def inverse_mod2(matrix):
+    """
+    Calculates the inverse of a binary matrix over the field of integers modulo 2.
+
+    Parameters
+    ----------
+    matrix : numpy.array
+        An invertible binary matrix.
+    result : numpy.array
+        The inverse of the given matrix over the field of integers modulo 2.
+
+    """
+    rows, columns = matrix.shape
+    if rows != columns:
+        raise Exception("The matrix must be a square matrix")
+    
+    matrix = matrix.copy()
+    result = np.eye(rows, dtype=int)
+
+    for i in range(rows):
+        # Find the pivot row
+        max_row = i + np.argmax(matrix[i:,i])
+
+         # Swap the current row with the pivot row
+        matrix[[i, max_row]] = matrix[[max_row, i]]
+        result[[i, max_row]] = result[[max_row, i]]
+        
+        # Eliminate all rows below the pivot
+        for j in range(i + 1, rows):
+            if matrix[j, i] == 1:
+                matrix[j] = (matrix[j] + matrix[i]) % 2
+                result[j] = (result[j] + result[i]) % 2
+
+    # Backward elimination
+    for i in range(rows - 1, -1, -1):
+        for j in range(i - 1, -1, -1):
+            if matrix[j, i] == 1:
+                matrix[j] = (matrix[j] + matrix[i]) % 2
+                result[j] = (result[j] + result[i]) % 2
+
+    return result
+
+
+def gaussian_elimination_mod2(matrix, type='row', show_pivots=False):
+    r"""
+    Performs Gaussian elimination in the field $F_2$.
+    
+    Parameters
+    ----------
+    matrix : numpy.array
+        A binary matrix.
+    type : str, optional
+        Available aore ``row`` for row echelon form, and ``column`` for column echelon form. 
+        The default is ``row``.
+    show_pivots : Boolean, optional
+        If ``True``, the pivot columns (rows) are returned.
+
+    Returns
+    -------
+    matrix : numpy.array
+        The row (column) echelon form of the given matrix.
+    pivots : list[int], optional
+        A list of indices for the pivot columns (rows).
+    
+    """
+    
+    matrix = matrix.copy()
+    
+    rows, columns = matrix.shape
+    column_index = 0
+    row_index = 0
+    pivots = [] # the pivot columns (type='row') or rows (type='column')
+
+    if type=='row':
+        while row_index<rows and column_index<columns:
+            # Find the pivot row
+            max_row = row_index + np.argmax(matrix[row_index:,column_index])
+
+            if matrix[max_row,column_index]==0:
+                column_index+=1
+            else:
+                # Pivot
+                pivots.append(column_index)
+
+                # Swap the current row with the pivot row
+                matrix[[row_index, max_row]] = matrix[[max_row, row_index]]
+        
+                # Eliminate all rows below the pivot
+                for j in range(row_index + 1, rows):
+                    if matrix[j, column_index] == 1:
+                        matrix[j] = (matrix[j] + matrix[row_index]) % 2
+
+                row_index+=1
+                column_index+=1
+
+    elif type=='column':
+        while row_index<rows and column_index<columns:
+            # Find the pivot column
+            max_column = column_index + np.argmax(matrix[row_index,column_index:])
+
+            if matrix[row_index,max_column]==0:
+                row_index+=1
+            else:
+                # Pivot
+                pivots.append(row_index)
+
+                # Swap the current column with the pivot column
+                matrix[:,[column_index, max_column]] = matrix[:,[max_column, column_index]]
+        
+                # Eliminate all rows below the pivot
+                for j in range(column_index + 1, columns):
+                    if matrix[row_index, j] == 1:
+                        matrix[:,j] = (matrix[:,j] + matrix[:,column_index]) % 2
+
+                row_index+=1
+                column_index+=1     
+
+    if show_pivots:
+        return matrix, pivots #(pivot_rows,pivot_cols)
+    else:
+        return matrix
+
+
+def construct_change_of_basis(S):
+    """
+    Implements the CZ construction outlined in https://quantum-journal.org/papers/q-2021-01-20-385/.
+
+    Parameters
+    ----------
+    S : numpy.array
+        A matrix representing a list of commuting Pauli operators: Each column is a Pauli operator in binary representation.
+
+    Returns
+    -------
+    A : numpy.array
+        Adjacency matrix for the graph state.
+    R_inv : numpy.array
+        A matrix representing a list of new Pauli operators.
+    h_list : numpy.array
+        A list indicating the qubits on which a Hadamard gate is applied.
+    s_list : numpy.array
+        A list indicating the qubits on which an S gate is applied.
+    perm_vec : numpy.array
+        A vector repesenting a permutation.
+
+    """
+
+    n = int(S.shape[0]/2)
+
+    ####################
+    # Step 0: Calculate S_0: Independent columns (i.e., Pauli terms) of S
+    ####################
+
+    S_reduced, independent_cols = gaussian_elimination_mod2(S, show_pivots=True)
+    k = len(independent_cols)
+
+    #S0 = np.vstack((z_matrix[:,independent_cols], x_matrix[:,independent_cols]))
+    S0 = S[:,independent_cols]
+
+    R0_inv = S_reduced[:k, :]
+
+    ####################
+    # Step 1: Calculate S_1: The first k rows of the X component have full rank k
+    ####################
+
+    # Find independent rows in X component of S0
+    S0X_reduced, independent_rows = gaussian_elimination_mod2(S0[-n:, :], type='column', show_pivots=True)
+
+    # Construnct S1 by applying a Hadamard (i.e., a swap) to the rows of S0 not in independent_rows
+    h_list = [i for i in range(n) if i not in independent_rows]
+    S1 = S0.copy()
+    for i in h_list:
+        S1[[i, n+i]] = S1[[n+i, i]]
+
+    # Find independent rows in X component of S1
+    S1X_reduced, independent_rows = gaussian_elimination_mod2(S1[-n:, :], type="column", show_pivots=True)
+
+    # Construct permutation achieving that the first k rows in X component of S1 are independent
+    perm = np.arange(0,n)
+    for index1, index2 in enumerate(independent_rows):
+        perm[index1] = index2
+        perm[index2] = index1
+
+    # Apply permutation to rows of S1 for Z and X component
+    S1 = np.vstack((S1[:n,:][perm],S1[-n:,:][perm]))
+
+    #################### 
+    # Step 2: Calculate S2: The first k rows of the X component are the identity matrix
+    ####################
+
+    R1_inv = S1[n:n+k,:]
+    R1 = inverse_mod2(R1_inv)
+    S2 = S1 @ R1 % 2
+
+    ####################
+    # Step 3: Calculate S3: Basis extension
+    ####################
+
+    C = S2[:k, :]
+    D = S2[k:n, :]
+    F = S2[-(n-k):, :]
+
+    S3 = np.block([[C, np.transpose(D)],
+                    [D, np.zeros((n-k,n-k), dtype=int)],
+                    [np.eye(k, dtype=int), np.zeros((k,n-k), dtype=int)],
+                    [F, np.eye(n-k, dtype=int)]])
+    R2_inv = np.block([[np.eye(k, dtype=int)],
+                       [np.zeros((n-k,k), dtype=int)]])   
+
+    ####################
+    # Step 4: Calculate S4 
+    ####################
+
+    R3_inv = S3[-n:,:]
+    R3 = inverse_mod2(R3_inv)
+
+    S4 = S3 @ R3 % 2
+
+    # Remove diagonal entries in upper left block 
+    s_list = []
+    for i in range(n):
+        if S4[:n, :][i,i]==1:
+            s_list.append(i)
+
+    Q2 = np.block([[np.eye(n, dtype=int),S4[:n, :]*np.eye(n, dtype=int)],
+                   [np.zeros((n,n), dtype=int),np.eye(n, dtype=int)]])
+    
+    S4 = Q2 @ S4 % 2
+
+    R_inv = R3_inv @ R2_inv @ R1_inv @ R0_inv % 2
+
+    # Adjacency matrix for the graph 
+    A = S4[:n, :]
+
+    return A, R_inv, h_list, s_list, perm
\ No newline at end of file
diff --git a/src/qrisp/operators/qubit/qubit_operator.py b/src/qrisp/operators/qubit/qubit_operator.py
index cc2f04a8..5e3c512c 100644
--- a/src/qrisp/operators/qubit/qubit_operator.py
+++ b/src/qrisp/operators/qubit/qubit_operator.py
@@ -19,7 +19,8 @@
 from qrisp.operators.hamiltonian import Hamiltonian
 from qrisp.operators.qubit.qubit_term import QubitTerm
 from qrisp.operators.qubit.measurement import get_measurement
-from qrisp import cx, h, s, x, sx_dg, IterationEnvironment, conjugate, merge
+from qrisp.operators.qubit.commutativity_tools import gaussian_elimination_mod2, inverse_mod2, construct_change_of_basis
+from qrisp import cx, cz, h, s, x, sx_dg, IterationEnvironment, conjugate, merge
 
 import sympy as sp
 import numpy as np
@@ -790,7 +791,7 @@ def commuting_qw_groups(self, show_bases=False, use_graph_coloring = True):
     # Measurement settings and measurement
     #
     
-    def change_of_basis(self, qarg):
+    def change_of_basis(self, qarg, method="commuting_qw"):
         """
         Performs several operations on a quantum argument such that the hermitian
         part of self is diagonal when conjugated with these operations.
@@ -799,6 +800,10 @@ def change_of_basis(self, qarg):
         ----------
         qarg : QuantumVariable or list[Qubit]
             The quantum argument to apply the change of basis on.
+        method : str, optional
+            The method for calculating the change of basis. 
+            Available are ``commuting`` (all QubitTerms must mutually commute) and ``commuting_qw`` (all QubitTerms must mutually commute qubit-wise).
+            The default is ``commuting_qw``.
 
         Returns
         -------
@@ -876,46 +881,119 @@ def change_of_basis(self, qarg):
         # We track which qubit is in which basis to raise an error if a
         # violation with the requirement of qubit wise commutativity is detected.
         basis_dict = {}
+
+        new_factor_dicts = []
+        prefactors = []
         
         ladder_conversion = {"A" : "P1", "C" : "P0"}
+
+        if method=="commuting_qw":
         
-        # We iterate through the terms and apply the appropriate basis transformation
-        for term, coeff in self.terms_dict.items():
+            # We iterate through the terms and apply the appropriate basis transformation
+            for term, coeff in self.terms_dict.items():
             
-            factor_dict = term.factor_dict
-            # This dictionary will contain the factors of the new term
-            new_factor_dict = {}
+                factor_dict = term.factor_dict
+                # This dictionary will contain the factors of the new term
+                new_factor_dict = {}
             
-            prefactor = 1
+                prefactor = 1
             
-            for j in range(n):
+                for j in range(n):
                 
-                # If there is no entry in the factor dict, this corresponds to
-                # identity => no basis change required.
-                if j not in factor_dict:
-                    continue
+                    # If there is no entry in the factor dict, this corresponds to
+                    # identity => no basis change required.
+                    if j not in factor_dict:
+                        continue
                 
-                # If j is already in the basis dict, we assert that the bases agree
-                # (otherwise there is a violation of qubit-wise commutativity)
-                if j in basis_dict:
-                    assert basis_dict[j] == factor_dict[j]
-                    continue
+                    # If j is already in the basis dict, we assert that the bases agree
+                    # (otherwise there is a violation of qubit-wise commutativity)
+                    if j in basis_dict:
+                        assert basis_dict[j] == factor_dict[j]
+                        continue
                 
-                # We treat ladder operators in the next section
-                if factor_dict[j] not in ["X", "Y", "Z"]:
-                    continue
+                    # We treat ladder operators in the next section
+                    if factor_dict[j] not in ["X", "Y", "Z"]:
+                        continue
                 
-                # Update the basis dict
-                basis_dict[j] = factor_dict[j]
+                    # Update the basis dict
+                    basis_dict[j] = factor_dict[j]
                 
-                # Append the appropriate basis-change gate
-                if factor_dict[j]=="X":
-                    h(qarg[j])
+                    # Append the appropriate basis-change gate
+                    if factor_dict[j]=="X":
+                        h(qarg[j])
                     
-                if factor_dict[j]=="Y":
-                    sx_dg(qarg[j])
+                    if factor_dict[j]=="Y":
+                        sx_dg(qarg[j])
+            
+                    new_factor_dict[j] = "Z"
+                
+                new_factor_dicts.append(new_factor_dict)
+                prefactors.append(prefactor)
+
+        if method=="commuting":
+
+            # Calculate S: Matrix where the colums correspond to the binary representation (Z/X) of the Pauli terms
+            x_vectors = []
+            z_vectors = []
+            coeffs = []
+            for term, coeff in self.terms_dict.items():
+                x_vector, z_vector = term.binary_representation(n)
+                x_vectors.append(x_vector)
+                z_vectors.append(z_vector)
+                coeffs.append(coeff)
+            x_matrix = np.stack(x_vectors, axis=1)
+            z_matrix = np.stack(z_vectors, axis=1)
+
+            S = np.vstack((z_matrix, x_matrix))
             
-                new_factor_dict[j] = "Z"
+            # Construct and apply change of basis
+            A, R_inv, h_list, s_list, perm = construct_change_of_basis(S)
+
+            def inv_graph_state(qarg):
+                for i in range(n):
+                    for j in range(i):
+                        if A[i,j]==1:
+                            cz(qarg[perm[i]],qarg[perm[j]])
+            h(qarg[:n])
+
+            def change_of_basis(qarg):
+                for i in h_list:
+                    h(qarg[perm[i]])
+                for i in s_list:
+                    s(qarg[perm[i]])
+                inv_graph_state(qarg)
+
+            change_of_basis(qarg)
+
+            # Construct new QubitOperator
+            #
+            # Factor (-1) appears if S gate is applied to X, or Hadamard gate H is applied to Y:
+            # S^dagger X S = -Y
+            # S^dagger Y S = X
+            # S^dagger Z S = Z
+            # H X H = Z
+            # H Y H = -Y
+            # H Z H = X
+            # For the original Pauli terms this translates to: Factor (-1) appears if S gate is applied to Y, or Hadamard gate H is applied to Y. (No factor (-1) occurs if S*H is applied.)
+
+            s_vector = np.zeros(n, dtype=int)
+            s_vector[s_list] = 1
+            h_vector = np.zeros(n, dtype=int)
+            h_vector[h_list] = 1
+            sh_vector = s_vector + h_vector % 2 
+            sign_vector = sh_vector @ (x_matrix*z_matrix) % 2
+
+            for index,z_vector in enumerate(R_inv.T):
+                new_factor_dict = {perm[i]:"Z" for i in range(n) if z_vector[i]==1}
+                new_factor_dicts.append(new_factor_dict)
+                prefactor = (-1)**sign_vector[index]
+                prefactors.append(prefactor)
+
+        # Ladder operators 
+        for term, coeff in self.terms_dict.items():    
+
+            prefactor = prefactors.pop(0)    
+            new_factor_dict = new_factor_dicts.pop(0)
             
             # Next we treat the ladder operators
             ladder_operators = [base for base in term.factor_dict.items() if base[1] in ["A", "C"]]
diff --git a/src/qrisp/operators/qubit/qubit_term.py b/src/qrisp/operators/qubit/qubit_term.py
index a7443854..8e467980 100644
--- a/src/qrisp/operators/qubit/qubit_term.py
+++ b/src/qrisp/operators/qubit/qubit_term.py
@@ -20,8 +20,16 @@
 from qrisp.operators.qubit.visualization import X_,Y_,Z_
 
 from sympy import Symbol
+import numpy as np
 
-PAULI_TABLE = {("I","I"):("I",1),("I","X"):("X",1),("I","Y"):("Y",1),("I","Z"):("Z",1),("I","A"):("A",1),("I","C"):("C",1),("I","P0"):("P0",1),("I","P1"):("P1",1),("X","I"):("X",1),("X","X"):("I",1),("X","Y"):("Z",1j),("X","Z"):("Y",(-0-1j)),("X","A"):("P0",1),("X","C"):("P1",1),("X","P0"):("A",1),("X","P1"):("C",1),("Y","I"):("Y",1),("Y","X"):("Z",(-0-1j)),("Y","Y"):("I",1),("Y","Z"):("X",1j),("Y","A"):("P0",(-0-1j)),("Y","C"):("P1",1j),("Y","P0"):("A",1j),("Y","P1"):("C",(-0-1j)),("Z","I"):("Z",1),("Z","X"):("Y",1j),("Z","Y"):("X",(-0-1j)),("Z","Z"):("I",1),("Z","A"):("A",-1),("Z","C"):("C",1),("Z","P0"):("P0",1),("Z","P1"):("P1",-1),("A","I"):("A",1),("A","X"):("P1",1),("A","Y"):("P1",(-0-1j)),("A","Z"):("A",1),("A","A"):("I",0),("A","C"):("P1",1),("A","P0"):("A",1),("A","P1"):("I",0),("C","I"):("C",1),("C","X"):("P0",1),("C","Y"):("P0",1j),("C","Z"):("C",-1),("C","A"):("P0",1),("C","C"):("I",0),("C","P0"):("I",0),("C","P1"):("C",1),("P0","I"):("P0",1),("P0","X"):("C",1),("P0","Y"):("C",(-0-1j)),("P0","Z"):("P0",1),("P0","A"):("I",0),("P0","C"):("C",1),("P0","P0"):("P0",1),("P0","P1"):("I",0),("P1","I"):("P1",1),("P1","X"):("A",1),("P1","Y"):("A",1j),("P1","Z"):("P1",-1),("P1","A"):("A",1),("P1","C"):("I",0),("P1","P0"):("I",0),("P1","P1"):("P1",1)}
+PAULI_TABLE = {("I","I"):("I",1),("I","X"):("X",1),("I","Y"):("Y",1),("I","Z"):("Z",1),("I","A"):("A",1),("I","C"):("C",1),("I","P0"):("P0",1),("I","P1"):("P1",1),
+            ("X","I"):("X",1),("X","X"):("I",1),("X","Y"):("Z",1j),("X","Z"):("Y",(-0-1j)),("X","A"):("P0",1),("X","C"):("P1",1),("X","P0"):("A",1),("X","P1"):("C",1),
+            ("Y","I"):("Y",1),("Y","X"):("Z",(-0-1j)),("Y","Y"):("I",1),("Y","Z"):("X",1j),("Y","A"):("P0",(-0-1j)),("Y","C"):("P1",1j),("Y","P0"):("A",1j),("Y","P1"):("C",(-0-1j)),
+            ("Z","I"):("Z",1),("Z","X"):("Y",1j),("Z","Y"):("X",(-0-1j)),("Z","Z"):("I",1),("Z","A"):("A",-1),("Z","C"):("C",1),("Z","P0"):("P0",1),("Z","P1"):("P1",-1),
+            ("A","I"):("A",1),("A","X"):("P1",1),("A","Y"):("P1",(-0-1j)),("A","Z"):("A",1),("A","A"):("I",0),("A","C"):("P1",1),("A","P0"):("A",1),("A","P1"):("I",0),
+            ("C","I"):("C",1),("C","X"):("P0",1),("C","Y"):("P0",1j),("C","Z"):("C",-1),("C","A"):("P0",1),("C","C"):("I",0),("C","P0"):("I",0),("C","P1"):("C",1),
+            ("P0","I"):("P0",1),("P0","X"):("C",1),("P0","Y"):("C",(-0-1j)),("P0","Z"):("P0",1),("P0","A"):("I",0),("P0","C"):("C",1),("P0","P0"):("P0",1),("P0","P1"):("I",0),
+            ("P1","I"):("P1",1),("P1","X"):("A",1),("P1","Y"):("A",1j),("P1","Z"):("P1",-1),("P1","A"):("A",1),("P1","C"):("I",0),("P1","P0"):("I",0),("P1","P1"):("P1",1)}
 
 #
 # QubitTerm
@@ -53,6 +61,21 @@ def copy(self):
     def is_identity(self):
         return len(self.factor_dict)==0
     
+    def binary_representation(self, n):
+        x_vector = np.zeros(n, dtype=int)
+        z_vector = np.zeros(n, dtype=int)
+        for index in range(n):
+            curr_factor = self.factor_dict.get(index,"I")
+            if curr_factor=="X":
+                x_vector[index] = 1
+            elif curr_factor=="Z":
+                z_vector[index] = 1
+            elif curr_factor=="Y":
+                x_vector[index] = 1
+                z_vector[index] = 1
+
+        return x_vector, z_vector
+    
     def serialize(self):
         # This function serializes the QubitTerm in a way that facilitates the 
         # measurement post-processing. To learn more details about the reasoning
@@ -479,12 +502,17 @@ def commute(self, other):
         keys.update(set(a.keys()))
         keys.update(set(b.keys()))
 
-        # Count non-commuting operators
+        # Count non-commuting X, Y, Z operators
         commute = True
 
         for key in keys:
-            if a.get(key,"I")!="I" and b.get(key,"I")!="I" and a.get(key,"I")!=b.get(key,"I"):
+            #if a.get(key,"I")!="I" and b.get(key,"I")!="I" and a.get(key,"I")!=b.get(key,"I"):
+            #    commute = not commute
+            if not PAULI_TABLE[a.get(key,"I"), b.get(key,"I")] == PAULI_TABLE[b.get(key,"I"), a.get(key,"I")]:
                 commute = not commute
+                if (a.get(key,"I") not in ["X","Y","Z"]) or (b.get(key,"I") not in ["X","Y","Z"]):
+                    return False
+  
         return commute
 
     def commute_qw(self, other):