From 9eb93cf0a23beb18235927a1cca5201bdb87f5d0 Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Wed, 12 Jul 2023 15:09:47 +0800
Subject: [PATCH 01/31] Initial commit

---
 src/qibotn/MpsHelper.py | 270 ++++++++++++++++++++++++++++++++++++++++
 src/qibotn/cutn.py      |  68 +++++++++-
 2 files changed, 336 insertions(+), 2 deletions(-)
 create mode 100644 src/qibotn/MpsHelper.py

diff --git a/src/qibotn/MpsHelper.py b/src/qibotn/MpsHelper.py
new file mode 100644
index 00000000..526cd3c3
--- /dev/null
+++ b/src/qibotn/MpsHelper.py
@@ -0,0 +1,270 @@
+# Copyright (c) 2021-2023, NVIDIA CORPORATION & AFFILIATES
+#
+# SPDX-License-Identifier: BSD-3-Clause
+
+import itertools
+
+import cupy as cp
+import numpy as np
+
+import cuquantum
+from cuquantum import cutensornet as cutn
+
+class MPSHelper:
+    
+    """
+    MPSHelper(num_sites, phys_extent, max_virtual_extent, initial_state, data_type, compute_type)
+
+    Create an MPSHelper object for gate splitting algorithm.
+           i              j       
+        -------A-------B-------                      i        j        k
+              p|       |q            ------->     -------A`-------B`-------
+               GGGGGGGGG                                r|        |s
+              r|       |s
+    
+    Args:
+        num_sites: The number of sites in the MPS.
+        phys_extents: The extent for the physical mode where the gate tensors are acted on.
+        max_virtual_extent: The maximal extent allowed for the virtual mode shared between adjacent MPS tensors.
+        initial_state: A sequence of :class:`cupy.ndarray` representing the initial state of the MPS.
+        data_type (cuquantum.cudaDataType): The data type for all tensors and gates. 
+        compute_type (cuquantum.ComputeType): The compute type for all gate splitting.
+    
+    """
+
+    def __init__(self, num_sites, phys_extent, max_virtual_extent, initial_state, data_type, compute_type):
+        self.num_sites = num_sites
+        self.phys_extent = phys_extent
+        self.data_type = data_type
+        self.compute_type = compute_type
+        
+        self.phys_modes = []
+        self.virtual_modes = []
+        self.new_mode = itertools.count(start=0, step=1)
+        
+        for i in range(num_sites+1):
+            self.virtual_modes.append(next(self.new_mode))
+            if i != num_sites:
+                self.phys_modes.append(next(self.new_mode))
+        
+        untruncated_max_extent = phys_extent ** (num_sites // 2)
+        if max_virtual_extent == 0:
+            self.max_virtual_extent = untruncated_max_extent
+        else:
+            self.max_virtual_extent = min(max_virtual_extent, untruncated_max_extent)
+        
+        self.handle = cutn.create()
+        self.work_desc = cutn.create_workspace_descriptor(self.handle)
+        self.svd_config = cutn.create_tensor_svd_config(self.handle)
+        self.svd_info = cutn.create_tensor_svd_info(self.handle)
+        self.gate_algo = cutn.GateSplitAlgo.DIRECT
+
+        self.desc_tensors = []
+        self.state_tensors = []
+
+        # create tensor descriptors
+        for i in range(self.num_sites):
+            temp = initial_state[i]
+            self.state_tensors.append(initial_state[i].astype(temp.dtype, order="F"))
+            extent = self.get_tensor_extent(i)
+            modes = self.get_tensor_modes(i)
+            desc_tensor = cutn.create_tensor_descriptor(self.handle, 3, extent, 0, modes, self.data_type)
+            self.desc_tensors.append(desc_tensor)
+    
+    def get_tensor(self, site):
+        """Get the tensor operands for a specific site."""
+        return self.state_tensors[site]
+    
+    def get_tensor_extent(self, site):
+        """Get the extent of the MPS tensor at a specific site."""
+        return self.state_tensors[site].shape
+    
+    def get_tensor_modes(self, site):
+        """Get the current modes of the MPS tensor at a specific site."""
+        return (self.virtual_modes[site], self.phys_modes[site], self.virtual_modes[site+1])
+
+    def set_svd_config(self, abs_cutoff, rel_cutoff, renorm, partition):
+        """Update the SVD truncation setting.
+        
+        Args:
+            abs_cutoff: The cutoff value for absolute singular value truncation.
+            rel_cutoff: The cutoff value for relative singular value truncation.
+            renorm (cuquantum.cutensornet.TensorSVDNormalization): The option for renormalization of the truncated singular values.
+            partition (cuquantum.cutensornet.TensorSVDPartition): The option for partitioning of the singular values.
+        """        
+        
+        if partition != cutn.TensorSVDPartition.UV_EQUAL:
+            raise NotImplementedError("this basic example expects partition to be cutensornet.TensorSVDPartition.UV_EQUAL")
+
+        svd_config_attributes = [cutn.TensorSVDConfigAttribute.ABS_CUTOFF, 
+                                 cutn.TensorSVDConfigAttribute.REL_CUTOFF, 
+                                 cutn.TensorSVDConfigAttribute.S_NORMALIZATION,
+                                 cutn.TensorSVDConfigAttribute.S_PARTITION]
+            
+        for (attr, value) in zip(svd_config_attributes, [abs_cutoff, rel_cutoff, renorm, partition]):
+            dtype = cutn.tensor_svd_config_get_attribute_dtype(attr)
+            value = np.array([value], dtype=dtype)
+            cutn.tensor_svd_config_set_attribute(self.handle, 
+                self.svd_config, attr, value.ctypes.data, value.dtype.itemsize)
+    
+    def set_gate_algorithm(self, gate_algo):
+        """Set the algorithm to use for all gate split operations.
+        
+        Args:
+            gate_algo (cuquantum.cutensornet.GateSplitAlgo): The gate splitting algorithm to use.
+        """
+
+        self.gate_algo = gate_algo
+
+    def compute_max_workspace_sizes(self):
+        """Compute the maximal workspace needed for MPS gating algorithm."""
+        modes_in_A = [ord(c) for c in ('i', 'p', 'j')]
+        modes_in_B = [ord(c) for c in ('j', 'q', 'k')]
+        modes_in_G = [ord(c) for c in ('p', 'q', 'r', 's')]
+        modes_out_A = [ord(c) for c in ('i', 'r', 'j')]
+        modes_out_B = [ord(c) for c in ('j', 's', 'k')]
+        
+        max_extents_AB = (self.max_virtual_extent, self.phys_extent, self.max_virtual_extent)
+        extents_in_G = (self.phys_extent, self.phys_extent, self.phys_extent, self.phys_extent)
+        
+        desc_tensor_in_A = cutn.create_tensor_descriptor(self.handle, 3, max_extents_AB, 0, modes_in_A, self.data_type)
+        desc_tensor_in_B = cutn.create_tensor_descriptor(self.handle, 3, max_extents_AB, 0, modes_in_B, self.data_type)
+        desc_tensor_in_G = cutn.create_tensor_descriptor(self.handle, 4, extents_in_G, 0, modes_in_G, self.data_type)
+        desc_tensor_out_A = cutn.create_tensor_descriptor(self.handle, 3, max_extents_AB, 0, modes_out_A, self.data_type)
+        desc_tensor_out_B = cutn.create_tensor_descriptor(self.handle, 3, max_extents_AB, 0, modes_out_B, self.data_type)
+        
+        cutn.workspace_compute_gate_split_sizes(self.handle, 
+            desc_tensor_in_A, desc_tensor_in_B, desc_tensor_in_G, 
+            desc_tensor_out_A, desc_tensor_out_B, 
+            self.gate_algo, self.svd_config, self.compute_type, self.work_desc)
+        
+        workspace_size = cutn.workspace_get_memory_size(self.handle, self.work_desc, cutn.WorksizePref.MIN, cutn.Memspace.DEVICE, cutn.WorkspaceKind.SCRATCH)
+
+        # free resources
+        cutn.destroy_tensor_descriptor(desc_tensor_in_A)
+        cutn.destroy_tensor_descriptor(desc_tensor_in_B)
+        cutn.destroy_tensor_descriptor(desc_tensor_in_G)
+        cutn.destroy_tensor_descriptor(desc_tensor_out_A)
+        cutn.destroy_tensor_descriptor(desc_tensor_out_B)
+        return workspace_size
+
+    def set_workspace(self, work, workspace_size):
+        """Compute the maximal workspace needed for MPS gating algorithm.
+        
+        Args:
+            work: Pointer to the allocated workspace.
+            workspace_size: The required workspace size on the device.
+        """
+        cutn.workspace_set_memory(self.handle, self.work_desc, cutn.Memspace.DEVICE, cutn.WorkspaceKind.SCRATCH, work.ptr, workspace_size)
+
+    def apply_gate(self, site_A, site_B, gate, verbose, stream):
+        """Inplace execution of the apply gate algoritm on site A and site B.
+        
+        Args:
+            site_A: The first site on which the gate is applied to.
+            site_B: The second site on which the gate is applied to.
+            gate (cupy.ndarray): The input data for the gate tensor.
+            verbose: Whether to print out the runtime information during truncation.
+            stream (cupy.cuda.Stream): The CUDA stream on which the computation is performed.
+        """
+        if site_B - site_A != 1:
+            raise ValueError("Site B must be the right site of site A")
+        if site_B >= self.num_sites:
+            raise ValueError("Site index cannot exceed maximum number of sites")
+        
+        desc_tensor_in_A = self.desc_tensors[site_A]
+        desc_tensor_in_B = self.desc_tensors[site_B]
+        
+        phys_mode_in_A = self.phys_modes[site_A]
+        phys_mode_in_B = self.phys_modes[site_B]
+        phys_mode_out_A = next(self.new_mode)
+        phys_mode_out_B = next(self.new_mode)
+        modes_G = (phys_mode_in_A, phys_mode_in_B, phys_mode_out_A, phys_mode_out_B)
+        extent_G = (self.phys_extent, self.phys_extent, self.phys_extent, self.phys_extent)
+        desc_tensor_in_G = cutn.create_tensor_descriptor(self.handle, 4, extent_G, 0, modes_G, self.data_type)
+
+        # construct and initialize the expected output A and B
+        tensor_in_A = self.state_tensors[site_A]
+        tensor_in_B = self.state_tensors[site_B]
+        left_extent_A = tensor_in_A.shape[0]
+        extent_AB_in = tensor_in_A.shape[2]
+        right_extent_B = tensor_in_B.shape[2]
+        combined_extent_left = min(left_extent_A, extent_AB_in * self.phys_extent) * self.phys_extent
+        combined_extent_right = min(right_extent_B, extent_AB_in * self.phys_extent) * self.phys_extent
+        extent_Aout_B = min(combined_extent_left, combined_extent_right, self.max_virtual_extent)
+        
+        extent_out_A = (left_extent_A, self.phys_extent, extent_Aout_B)
+        extent_out_B = (extent_Aout_B, self.phys_extent, right_extent_B)
+
+        tensor_out_A = cp.zeros(extent_out_A, dtype=tensor_in_A.dtype, order="F")
+        tensor_out_B = cp.zeros(extent_out_B, dtype=tensor_in_B.dtype, order="F")
+
+        # create tensor descriptors for output A and B
+        modes_out_A = (self.virtual_modes[site_A], phys_mode_out_A, self.virtual_modes[site_A+1])
+        modes_out_B = (self.virtual_modes[site_B], phys_mode_out_B, self.virtual_modes[site_B+1])
+
+        desc_tensor_out_A = cutn.create_tensor_descriptor(self.handle, 3, extent_out_A, 0, modes_out_A, self.data_type)
+        desc_tensor_out_B = cutn.create_tensor_descriptor(self.handle, 3, extent_out_B, 0, modes_out_B, self.data_type)
+
+        cutn.gate_split(self.handle, 
+                        desc_tensor_in_A, tensor_in_A.data.ptr, 
+                        desc_tensor_in_B, tensor_in_B.data.ptr, 
+                        desc_tensor_in_G, gate.data.ptr, 
+                        desc_tensor_out_A, tensor_out_A.data.ptr, 
+                        0, # we factorize singular values equally onto output A and B.
+                        desc_tensor_out_B, tensor_out_B.data.ptr, 
+                        self.gate_algo, self.svd_config, self.compute_type, 
+                        self.svd_info, self.work_desc, stream.ptr)
+        
+        if verbose:
+            full_extent = np.array([0], dtype=cutn.tensor_svd_info_get_attribute_dtype(cutn.TensorSVDInfoAttribute.FULL_EXTENT))
+            reduced_extent = np.array([0], dtype=cutn.tensor_svd_info_get_attribute_dtype(cutn.TensorSVDInfoAttribute.REDUCED_EXTENT))
+            discarded_weight = np.array([0], dtype=cutn.tensor_svd_info_get_attribute_dtype(cutn.TensorSVDInfoAttribute.DISCARDED_WEIGHT))
+            
+            cutn.tensor_svd_info_get_attribute(
+                self.handle, self.svd_info, cutn.TensorSVDInfoAttribute.FULL_EXTENT, 
+                full_extent.ctypes.data, full_extent.dtype.itemsize)
+            cutn.tensor_svd_info_get_attribute(
+                self.handle, self.svd_info, cutn.TensorSVDInfoAttribute.REDUCED_EXTENT, 
+                reduced_extent.ctypes.data, reduced_extent.dtype.itemsize)
+            cutn.tensor_svd_info_get_attribute(
+                self.handle, self.svd_info, cutn.TensorSVDInfoAttribute.DISCARDED_WEIGHT, 
+                discarded_weight.ctypes.data, discarded_weight.dtype.itemsize)
+            
+            print("Virtual bond truncated from {0} to {1} with a discarded weight of {2:.6f}".format(full_extent[0], reduced_extent[0], discarded_weight[0]))
+
+        self.phys_modes[site_A] = phys_mode_out_A
+        self.phys_modes[site_B] = phys_mode_out_B
+        self.desc_tensors[site_A] = desc_tensor_out_A
+        self.desc_tensors[site_B] = desc_tensor_out_B
+        
+        extent_out_A = np.zeros((3,), dtype=np.int64)
+        extent_out_B = np.zeros((3,), dtype=np.int64)
+        extent_out_A, strides_out_A = cutn.get_tensor_details(self.handle, desc_tensor_out_A)[2:]
+        extent_out_B, strides_out_B = cutn.get_tensor_details(self.handle, desc_tensor_out_B)[2:]
+
+        # Recall that `cutensornet.gate_split` can potentially find reduced extent during SVD truncation when value-based truncation is used.
+        # Therefore we here update the container for output tensor A and B.
+        if extent_out_A[2] != extent_Aout_B:
+            # note strides in cutensornet are in the unit of count and strides in cupy/numpy are in the unit of nbytes
+            strides_out_A = [i * tensor_out_A.itemsize for i in strides_out_A]
+            strides_out_B = [i * tensor_out_B.itemsize for i in strides_out_B]
+            tensor_out_A = cp.ndarray(extent_out_A, dtype=tensor_out_A.dtype, memptr=tensor_out_A.data, strides=strides_out_A)
+            tensor_out_B = cp.ndarray(extent_out_B, dtype=tensor_out_B.dtype, memptr=tensor_out_B.data, strides=strides_out_B)
+
+        self.state_tensors[site_A] = tensor_out_A
+        self.state_tensors[site_B] = tensor_out_B
+
+        cutn.destroy_tensor_descriptor(desc_tensor_in_A)
+        cutn.destroy_tensor_descriptor(desc_tensor_in_B)
+        cutn.destroy_tensor_descriptor(desc_tensor_in_G)
+    
+    def __del__(self):
+        """Free all resources owned by the object."""
+        for desc_tensor in self.desc_tensors:
+            cutn.destroy_tensor_descriptor(desc_tensor)
+        cutn.destroy(self.handle)
+        cutn.destroy_workspace_descriptor(self.work_desc)
+        cutn.destroy_tensor_svd_config(self.svd_config)
+        cutn.destroy_tensor_svd_info(self.svd_info)
+
diff --git a/src/qibotn/cutn.py b/src/qibotn/cutn.py
index e6f3e8c3..6e4a98bd 100644
--- a/src/qibotn/cutn.py
+++ b/src/qibotn/cutn.py
@@ -1,8 +1,72 @@
-# from qibotn import quimb as qiboquimb
 from qibotn.QiboCircuitConvertor import QiboCircuitToEinsum
 from cuquantum import contract
-
+from MpsHelper import MPSHelper
+import cuquantum 
+from cuquantum import cutensornet as cutn
+import cupy as cp
+import numpy as np
 
 def eval(qibo_circ, datatype):
     myconvertor = QiboCircuitToEinsum(qibo_circ, dtype=datatype)
     return contract(*myconvertor.state_vector_operands())
+
+def eval_mps(qibo_circ, datatype):
+    #Create MPS
+    cutensornet.create()
+    return contract()
+
+
+if __name__ == "__main__":
+    print("cuTensorNet-vers:", cutn.get_version())
+    dev = cp.cuda.Device()  # get current device
+    props = cp.cuda.runtime.getDeviceProperties(dev.id)
+    print("===== device info ======")
+    print("GPU-name:", props["name"].decode())
+    print("GPU-clock:", props["clockRate"])
+    print("GPU-memoryClock:", props["memoryClockRate"])
+    print("GPU-nSM:", props["multiProcessorCount"])
+    print("GPU-major:", props["major"])
+    print("GPU-minor:", props["minor"])
+    print("========================")
+    
+    data_type = cuquantum.cudaDataType.CUDA_C_64F
+    compute_type = cuquantum.ComputeType.COMPUTE_64F
+    num_sites = 16
+    phys_extent = 2
+    max_virtual_extent = 12
+    
+    ## we initialize the MPS state as a product state |000...000>
+    initial_state = []
+    for i in range(num_sites):
+        # we create dummpy indices for MPS tensors on the boundary for easier bookkeeping
+        # we'll use Fortran layout throughout this example
+        tensor = cp.zeros((1,2,1), dtype=np.complex128, order="F")
+        tensor[0,0,0] = 1.0
+        initial_state.append(tensor)
+    
+    mps_helper = MPSHelper(num_sites, phys_extent, max_virtual_extent, initial_state, data_type, compute_type)
+
+    ##################################
+    # Setup options for gate operation
+    ##################################
+    
+    abs_cutoff = 1e-2
+    rel_cutoff = 1e-2
+    renorm = cutn.TensorSVDNormalization.L2
+    partition = cutn.TensorSVDPartition.UV_EQUAL
+    mps_helper.set_svd_config(abs_cutoff, rel_cutoff, renorm, partition)
+
+    gate_algo = cutn.GateSplitAlgo.REDUCED
+    mps_helper.set_gate_algorithm(gate_algo)
+    
+    #####################################
+    # Workspace estimation and allocation
+    #####################################
+
+    free_mem, total_mem = dev.mem_info
+    worksize = free_mem *.7
+    required_workspace_size = mps_helper.compute_max_workspace_sizes()
+    work = cp.cuda.alloc(worksize)
+    print(f"Maximal workspace size requried: {required_workspace_size / 1024 ** 3:.3f} GB")
+    mps_helper.set_workspace(work, required_workspace_size)
+    

From 76f61bc9fe0754d00c0d43cb33f605146a91d9c0 Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Fri, 14 Jul 2023 09:50:43 +0800
Subject: [PATCH 02/31] Updated cuquantum requirement

---
 setup.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/setup.py b/setup.py
index 2b1184ce..43761fc5 100644
--- a/setup.py
+++ b/setup.py
@@ -43,7 +43,7 @@ def version():
         "qibojit>=0.0.7",
         "quimb[tensor]>=1.4.0",
         "cupy>=11.6.0",
-        "cuquantum-python-cu11",
+        "cuquantum-python-cu11>=23.3.0",
     ],
     extras_require={
         "docs": [],

From 3cb0fec99c2b54a2d48c7401de224736d04091f9 Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Fri, 14 Jul 2023 09:51:06 +0800
Subject: [PATCH 03/31] Added MPS codes

---
 src/qibotn/MPSContractionHelper.py | 155 +++++++++++++++++++++++++++++
 src/qibotn/MPSUtils.py             |  78 +++++++++++++++
 src/qibotn/QiboCircuitToMPS.py     |  24 +++++
 src/qibotn/cutn.py                 |  78 ++++-----------
 4 files changed, 278 insertions(+), 57 deletions(-)
 create mode 100644 src/qibotn/MPSContractionHelper.py
 create mode 100644 src/qibotn/MPSUtils.py
 create mode 100644 src/qibotn/QiboCircuitToMPS.py

diff --git a/src/qibotn/MPSContractionHelper.py b/src/qibotn/MPSContractionHelper.py
new file mode 100644
index 00000000..cd0905e7
--- /dev/null
+++ b/src/qibotn/MPSContractionHelper.py
@@ -0,0 +1,155 @@
+from cuquantum import contract, contract_path, CircuitToEinsum, tensor
+
+class MPSContractionHelper:
+    """
+    A helper class to compute various quantities for a given MPS.
+    
+    Interleaved format is used to construct the input args for `cuquantum.contract`. 
+    A concrete example on how the modes are populated for a 7-site MPS is provided below:
+    
+          0     2     4     6     8    10     12    14        
+    bra -----A-----B-----C-----D-----E-----F-----G-----
+             |     |     |     |     |     |     |     
+            1|    3|    5|    7|    9|   11|   13|     
+             |     |     |     |     |     |     |     
+    ket -----a-----b-----c-----d-----e-----f-----g-----
+          15    16    17    18    19    20    21    22
+    
+    
+    The follwing compute quantities are supported:
+    
+        - the norm of the MPS.
+        - the equivalent state vector from the MPS.
+        - the expectation value for a given operator.
+        - the equivalent state vector after multiplying an MPO to an MPS.
+    
+    Note that for the nth MPS tensor (rank-3), the modes of the tensor are expected to be `(i,p,j)` 
+    where i denotes the bonding mode with the (n-1)th tensor, p denotes the physical mode for the qubit and 
+    j denotes the bonding mode with the (n+1)th tensor.
+    
+    Args:
+        num_qubits: The number of qubits for the MPS.
+    """
+    
+    def __init__(self, num_qubits):
+        self.num_qubits = num_qubits
+        self.path_cache = dict()
+        self.bra_modes = [(2*i, 2*i+1, 2*i+2) for i in range(num_qubits)]
+        offset = 2*num_qubits+1
+        self.ket_modes = [(i+offset, 2*i+1, i+1+offset) for i in range(num_qubits)]
+    
+    def contract_norm(self, mps_tensors, options=None):
+        """
+        Contract the corresponding tensor network to form the norm of the MPS.
+
+        Args:
+            mps_tensors: A list of rank-3 ndarray-like tensor objects. 
+                The indices of the ith tensor are expected to be bonding index to the i-1 tensor, 
+                the physical mode, and then the bonding index to the i+1th tensor.
+            options: Specify the contract and decompose options. 
+
+        Returns:
+            The norm of the MPS.
+        """
+        interleaved_inputs = []
+        for i, o in enumerate(mps_tensors):
+            interleaved_inputs.extend([o, self.bra_modes[i], o.conj(), self.ket_modes[i]])
+        interleaved_inputs.append([]) # output
+        return self._contract('norm', interleaved_inputs, options=options).real
+
+    def contract_state_vector(self, mps_tensors, options=None):
+        """
+        Contract the corresponding tensor network to form the state vector representation of the MPS.
+
+        Args:
+            mps_tensors: A list of rank-3 ndarray-like tensor objects. 
+                The indices of the ith tensor are expected to be bonding index to the i-1 tensor, 
+                the physical mode, and then the bonding index to the i+1th tensor.
+            options: Specify the contract and decompose options. 
+
+        Returns:
+            An ndarray-like object as the state vector.
+        """
+        interleaved_inputs = []
+        for i, o in enumerate(mps_tensors):
+            interleaved_inputs.extend([o, self.bra_modes[i]])
+        output_modes = tuple([bra_modes[1] for bra_modes in self.bra_modes])
+        interleaved_inputs.append(output_modes) # output
+        return self._contract('sv', interleaved_inputs, options=options)
+    
+    def contract_expectation(self, mps_tensors, operator, qubits, options=None, normalize=False):
+        """
+        Contract the corresponding tensor network to form the state vector representation of the MPS.
+
+        Args:
+            mps_tensors: A list of rank-3 ndarray-like tensor objects. 
+                The indices of the ith tensor are expected to be bonding index to the i-1 tensor, 
+                the physical mode, and then the bonding index to the i+1th tensor.
+            operator: A ndarray-like tensor object. 
+                The modes of the operator are expected to be output qubits followed by input qubits, e.g, 
+                ``A, B, a, b`` where `a, b` denotes the inputs and `A, B'` denotes the outputs. 
+            qubits: A sequence of integers specifying the qubits that the operator is acting on. 
+            options: Specify the contract and decompose options. 
+            normalize: Whether to scale the expectation value by the normalization factor.
+
+        Returns:
+            An ndarray-like object as the state vector.
+        """
+        
+        interleaved_inputs = []
+        extra_mode = 3 * self.num_qubits + 2
+        operator_modes = [None] * len(qubits) + [self.bra_modes[q][1] for q in qubits]
+        qubits = list(qubits)
+        for i, o in enumerate(mps_tensors):
+            interleaved_inputs.extend([o, self.bra_modes[i]])
+            k_modes = self.ket_modes[i]
+            if i in qubits:
+                k_modes = (k_modes[0], extra_mode, k_modes[2])
+                q = qubits.index(i)
+                operator_modes[q] = extra_mode # output modes
+                extra_mode += 1
+            interleaved_inputs.extend([o.conj(), k_modes])
+        interleaved_inputs.extend([operator, tuple(operator_modes)])
+        interleaved_inputs.append([]) # output
+        if normalize:
+            norm = self.contract_norm(mps_tensors, options=options)
+        else:
+            norm = 1
+        return self._contract(f'exp{qubits}', interleaved_inputs, options=options) / norm
+    
+    def contract_mps_mpo_to_state_vector(self, mps_tensors, mpo_tensors, options=None):
+        """
+        Contract the corresponding tensor network to form the output state vector from applying the MPO to the MPS.
+
+        Args:
+            mps_tensors: A list of rank-3 ndarray-like tensor objects. 
+                The indices of the ith tensor are expected to be the bonding index to the i-1 tensor, 
+                the physical mode, and then the bonding index to the i+1th tensor.
+            mpo_tensors: A list of rank-4 ndarray-like tensor objects.
+                The indics of the ith tensor are expected to be the bonding index to the i-1 tensor, 
+                the output physical mode, the bonding index to the i+1th tensor and then the inputput physical mode.
+            options: Specify the contract and decompose options. 
+
+        Returns:
+            An ndarray-like object as the output state vector.
+        """
+        interleaved_inputs = []
+        for i, o in enumerate(mps_tensors):
+            interleaved_inputs.extend([o, self.bra_modes[i]])
+        output_modes = []
+        offset = 2 * self.num_qubits + 1
+        for i, o in enumerate(mpo_tensors):
+            mpo_modes = (2*i+offset, 2*i+offset+1, 2*i+offset+2, 2*i+1)
+            output_modes.append(2*i+offset+1)
+            interleaved_inputs.extend([o, mpo_modes])
+        interleaved_inputs.append(output_modes)
+        return self._contract('mps_mpo', interleaved_inputs, options=options)
+    
+    def _contract(self, key, interleaved_inputs, options=None):
+        """
+        Perform the contraction task given interleaved inputs. Path will be cached.
+        """
+        if key not in self.path_cache:
+            self.path_cache[key] = contract_path(*interleaved_inputs, options=options)[0]
+        path = self.path_cache[key]
+        return contract(*interleaved_inputs, options=options, optimize={'path':path})
diff --git a/src/qibotn/MPSUtils.py b/src/qibotn/MPSUtils.py
new file mode 100644
index 00000000..cb77e4b9
--- /dev/null
+++ b/src/qibotn/MPSUtils.py
@@ -0,0 +1,78 @@
+import cupy as cp
+from cuquantum.cutensornet.experimental import contract_decompose
+from cuquantum import contract
+
+def get_initial_mps(num_qubits, dtype='complex128'):
+    """
+    Generate the MPS with an initial state of |00...00> 
+    """
+    state_tensor = cp.asarray([1, 0], dtype=dtype).reshape(1,2,1)
+    mps_tensors = [state_tensor] * num_qubits
+    return mps_tensors
+
+def mps_site_right_swap(
+    mps_tensors, 
+    i, 
+    algorithm=None, 
+    options=None
+):
+    """
+    Perform the swap operation between the ith and i+1th MPS tensors.
+    """
+    # contraction followed by QR decomposition
+    a, _, b = contract_decompose('ipj,jqk->iqj,jpk', *mps_tensors[i:i+2], algorithm=algorithm, options=options)
+    mps_tensors[i:i+2] = (a, b)
+    return mps_tensors
+
+def apply_gate(
+    mps_tensors, 
+    gate, 
+    qubits, 
+    algorithm=None, 
+    options=None
+):
+    """
+    Apply the gate operand to the MPS tensors in-place.
+    
+    Args:
+        mps_tensors: A list of rank-3 ndarray-like tensor objects. 
+            The indices of the ith tensor are expected to be the bonding index to the i-1 tensor, 
+            the physical mode, and then the bonding index to the i+1th tensor.
+        gate: A ndarray-like tensor object representing the gate operand. 
+            The modes of the gate is expected to be output qubits followed by input qubits, e.g, 
+            ``A, B, a, b`` where ``a, b`` denotes the inputs and ``A, B`` denotes the outputs. 
+        qubits: A sequence of integers denoting the qubits that the gate is applied onto.
+        algorithm: The contract and decompose algorithm to use for gate application. 
+            Can be either a `dict` or a `ContractDecomposeAlgorithm`.
+        options: Specify the contract and decompose options. 
+    
+    Returns:
+        The updated MPS tensors.
+    """
+    
+    n_qubits = len(qubits)
+    if n_qubits == 1:
+        # single-qubit gate
+        i = qubits[0]
+        mps_tensors[i] = contract('ipj,qp->iqj', mps_tensors[i], gate, options=options) # in-place update
+    elif n_qubits == 2:
+        # two-qubit gate
+        i, j = qubits
+        if i > j:
+            # swap qubits order
+            return apply_gate(mps_tensors, gate.transpose(1,0,3,2), (j, i), algorithm=algorithm, options=options)
+        elif i+1 == j:
+            # two adjacent qubits
+            a, _, b = contract_decompose('ipj,jqk,rspq->irj,jsk', *mps_tensors[i:i+2], gate, algorithm=algorithm, options=options)
+            mps_tensors[i:i+2] = (a, b) # in-place update
+        else:
+            # non-adjacent two-qubit gate
+            # step 1: swap i with i+1
+            mps_site_right_swap(mps_tensors, i, algorithm=algorithm, options=options)
+            # step 2: apply gate to (i+1, j) pair. This amounts to a recursive swap until the two qubits are adjacent
+            apply_gate(mps_tensors, gate, (i+1, j), algorithm=algorithm, options=options) 
+            # step 3: swap back i and i+1
+            mps_site_right_swap(mps_tensors, i, algorithm=algorithm, options=options)
+    else:
+        raise NotImplementedError("Only one- and two-qubit gates supported")
+    return mps_tensors
\ No newline at end of file
diff --git a/src/qibotn/QiboCircuitToMPS.py b/src/qibotn/QiboCircuitToMPS.py
new file mode 100644
index 00000000..f66804de
--- /dev/null
+++ b/src/qibotn/QiboCircuitToMPS.py
@@ -0,0 +1,24 @@
+import cupy as cp
+import numpy as np
+
+from cuquantum import cutensornet as cutn
+from QiboCircuitConvertor import QiboCircuitToEinsum
+from MPSUtils import get_initial_mps, apply_gate
+
+class QiboCircuitToMPS:
+    def __init__(self,circ_qibo, gate_algo, dtype = 'complex128',rand_seed=0,):
+        np.random.seed(rand_seed)
+        cp.random.seed(rand_seed)
+        
+        self.num_qubits = circ_qibo.nqubits
+        self.handle = cutn.create()
+        self.options = {'handle': self.handle}
+        self.dtype = dtype
+        self.mps_tensors = get_initial_mps(self.num_qubits, dtype=dtype)
+        circuitconvertor = QiboCircuitToEinsum(circ_qibo)
+
+        for (gate, qubits) in circuitconvertor.gate_tensors:
+            # mapping from qubits to qubit indices
+            # apply the gate in-place
+            apply_gate(self.mps_tensors, gate, qubits, algorithm=gate_algo, options=self.options)
+
diff --git a/src/qibotn/cutn.py b/src/qibotn/cutn.py
index 6e4a98bd..f29a62ed 100644
--- a/src/qibotn/cutn.py
+++ b/src/qibotn/cutn.py
@@ -5,68 +5,32 @@
 from cuquantum import cutensornet as cutn
 import cupy as cp
 import numpy as np
+from qibo.models import QFT
+from QiboCircuitToMPS import QiboCircuitToMPS
+from MPSContractionHelper import MPSContractionHelper
+
 
 def eval(qibo_circ, datatype):
     myconvertor = QiboCircuitToEinsum(qibo_circ, dtype=datatype)
     return contract(*myconvertor.state_vector_operands())
 
-def eval_mps(qibo_circ, datatype):
-    #Create MPS
-    cutensornet.create()
-    return contract()
-
+def eval_mps(qibo_circ, gate_algo, datatype):
+    myconvertor = QiboCircuitToMPS(qibo_circ, gate_algo, datatype)  
+    mps_helper = MPSContractionHelper(myconvertor.num_qubits)
+    sv_mps = mps_helper.contract_state_vector(myconvertor.mps_tensors,myconvertor.options)
+    return sv_mps
 
 if __name__ == "__main__":
-    print("cuTensorNet-vers:", cutn.get_version())
-    dev = cp.cuda.Device()  # get current device
-    props = cp.cuda.runtime.getDeviceProperties(dev.id)
-    print("===== device info ======")
-    print("GPU-name:", props["name"].decode())
-    print("GPU-clock:", props["clockRate"])
-    print("GPU-memoryClock:", props["memoryClockRate"])
-    print("GPU-nSM:", props["multiProcessorCount"])
-    print("GPU-major:", props["major"])
-    print("GPU-minor:", props["minor"])
-    print("========================")
-    
-    data_type = cuquantum.cudaDataType.CUDA_C_64F
-    compute_type = cuquantum.ComputeType.COMPUTE_64F
-    num_sites = 16
-    phys_extent = 2
-    max_virtual_extent = 12
-    
-    ## we initialize the MPS state as a product state |000...000>
-    initial_state = []
-    for i in range(num_sites):
-        # we create dummpy indices for MPS tensors on the boundary for easier bookkeeping
-        # we'll use Fortran layout throughout this example
-        tensor = cp.zeros((1,2,1), dtype=np.complex128, order="F")
-        tensor[0,0,0] = 1.0
-        initial_state.append(tensor)
-    
-    mps_helper = MPSHelper(num_sites, phys_extent, max_virtual_extent, initial_state, data_type, compute_type)
-
-    ##################################
-    # Setup options for gate operation
-    ##################################
-    
-    abs_cutoff = 1e-2
-    rel_cutoff = 1e-2
-    renorm = cutn.TensorSVDNormalization.L2
-    partition = cutn.TensorSVDPartition.UV_EQUAL
-    mps_helper.set_svd_config(abs_cutoff, rel_cutoff, renorm, partition)
-
-    gate_algo = cutn.GateSplitAlgo.REDUCED
-    mps_helper.set_gate_algorithm(gate_algo)
-    
-    #####################################
-    # Workspace estimation and allocation
-    #####################################
-
-    free_mem, total_mem = dev.mem_info
-    worksize = free_mem *.7
-    required_workspace_size = mps_helper.compute_max_workspace_sizes()
-    work = cp.cuda.alloc(worksize)
-    print(f"Maximal workspace size requried: {required_workspace_size / 1024 ** 3:.3f} GB")
-    mps_helper.set_workspace(work, required_workspace_size)
+    num_qubits  = 25
+    swaps = True
+    circ_qibo = QFT(num_qubits, swaps)
     
+    exact_gate_algorithm = {'qr_method': False, 
+                            'svd_method':{'partition': 'UV', 'abs_cutoff':1e-12}}
+    dtype="complex128"
+    sv_mps = eval_mps(circ_qibo, exact_gate_algorithm, dtype)
+    sv_reference = eval(circ_qibo, dtype)
+    state_vec = np.array(circ_qibo())
+    print(f"State vector difference: {abs(sv_mps-sv_reference).max():0.3e}")
+    assert cp.allclose(sv_mps, sv_reference)
+    assert cp.allclose(sv_mps.flatten(), state_vec)

From a17d8e6b7828439233d502e4b1b6d3a5f251af03 Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Mon, 17 Jul 2023 14:33:38 +0800
Subject: [PATCH 04/31] free handle

---
 src/qibotn/QiboCircuitToMPS.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/qibotn/QiboCircuitToMPS.py b/src/qibotn/QiboCircuitToMPS.py
index f66804de..49aa125c 100644
--- a/src/qibotn/QiboCircuitToMPS.py
+++ b/src/qibotn/QiboCircuitToMPS.py
@@ -22,3 +22,5 @@ def __init__(self,circ_qibo, gate_algo, dtype = 'complex128',rand_seed=0,):
             # apply the gate in-place
             apply_gate(self.mps_tensors, gate, qubits, algorithm=gate_algo, options=self.options)
 
+    def __del__(self):
+        cutn.destroy(self.handle)

From b043e6a0173dc508a3cb9cbef10ca863be478064 Mon Sep 17 00:00:00 2001
From: Liwei Yang <yang0345@e.ntu.edu.sg>
Date: Mon, 24 Jul 2023 17:29:29 +0800
Subject: [PATCH 05/31] Add the pytest function for MPS in cuquantum.

---
 src/qibotn/QiboCircuitToMPS.py           | 12 ++++----
 src/qibotn/cutn.py                       | 23 +++++++-------
 tests/test_cuquantum_cutensor_backend.py | 38 ++++++++++++++++++++++++
 3 files changed, 57 insertions(+), 16 deletions(-)

diff --git a/src/qibotn/QiboCircuitToMPS.py b/src/qibotn/QiboCircuitToMPS.py
index 49aa125c..4c18aac0 100644
--- a/src/qibotn/QiboCircuitToMPS.py
+++ b/src/qibotn/QiboCircuitToMPS.py
@@ -2,14 +2,15 @@
 import numpy as np
 
 from cuquantum import cutensornet as cutn
-from QiboCircuitConvertor import QiboCircuitToEinsum
-from MPSUtils import get_initial_mps, apply_gate
+from qibotn.QiboCircuitConvertor import QiboCircuitToEinsum
+from qibotn.MPSUtils import get_initial_mps, apply_gate
+
 
 class QiboCircuitToMPS:
-    def __init__(self,circ_qibo, gate_algo, dtype = 'complex128',rand_seed=0,):
+    def __init__(self, circ_qibo, gate_algo, dtype='complex128', rand_seed=0,):
         np.random.seed(rand_seed)
         cp.random.seed(rand_seed)
-        
+
         self.num_qubits = circ_qibo.nqubits
         self.handle = cutn.create()
         self.options = {'handle': self.handle}
@@ -20,7 +21,8 @@ def __init__(self,circ_qibo, gate_algo, dtype = 'complex128',rand_seed=0,):
         for (gate, qubits) in circuitconvertor.gate_tensors:
             # mapping from qubits to qubit indices
             # apply the gate in-place
-            apply_gate(self.mps_tensors, gate, qubits, algorithm=gate_algo, options=self.options)
+            apply_gate(self.mps_tensors, gate, qubits,
+                       algorithm=gate_algo, options=self.options)
 
     def __del__(self):
         cutn.destroy(self.handle)
diff --git a/src/qibotn/cutn.py b/src/qibotn/cutn.py
index f29a62ed..21c316fc 100644
--- a/src/qibotn/cutn.py
+++ b/src/qibotn/cutn.py
@@ -1,33 +1,34 @@
 from qibotn.QiboCircuitConvertor import QiboCircuitToEinsum
 from cuquantum import contract
-from MpsHelper import MPSHelper
-import cuquantum 
 from cuquantum import cutensornet as cutn
 import cupy as cp
 import numpy as np
 from qibo.models import QFT
-from QiboCircuitToMPS import QiboCircuitToMPS
-from MPSContractionHelper import MPSContractionHelper
+from qibotn.QiboCircuitToMPS import QiboCircuitToMPS
+from qibotn.MPSContractionHelper import MPSContractionHelper
 
 
 def eval(qibo_circ, datatype):
     myconvertor = QiboCircuitToEinsum(qibo_circ, dtype=datatype)
     return contract(*myconvertor.state_vector_operands())
 
+
 def eval_mps(qibo_circ, gate_algo, datatype):
-    myconvertor = QiboCircuitToMPS(qibo_circ, gate_algo, datatype)  
+    myconvertor = QiboCircuitToMPS(qibo_circ, gate_algo, dtype=datatype)
     mps_helper = MPSContractionHelper(myconvertor.num_qubits)
-    sv_mps = mps_helper.contract_state_vector(myconvertor.mps_tensors,myconvertor.options)
+    sv_mps = mps_helper.contract_state_vector(
+        myconvertor.mps_tensors, myconvertor.options)
     return sv_mps
 
+
 if __name__ == "__main__":
-    num_qubits  = 25
+    num_qubits = 25
     swaps = True
     circ_qibo = QFT(num_qubits, swaps)
-    
-    exact_gate_algorithm = {'qr_method': False, 
-                            'svd_method':{'partition': 'UV', 'abs_cutoff':1e-12}}
-    dtype="complex128"
+
+    exact_gate_algorithm = {'qr_method': False,
+                            'svd_method': {'partition': 'UV', 'abs_cutoff': 1e-12}}
+    dtype = "complex128"
     sv_mps = eval_mps(circ_qibo, exact_gate_algorithm, dtype)
     sv_reference = eval(circ_qibo, dtype)
     state_vec = np.array(circ_qibo())
diff --git a/tests/test_cuquantum_cutensor_backend.py b/tests/test_cuquantum_cutensor_backend.py
index e7f28043..1df3d0f9 100644
--- a/tests/test_cuquantum_cutensor_backend.py
+++ b/tests/test_cuquantum_cutensor_backend.py
@@ -2,6 +2,7 @@
 
 import config
 import numpy as np
+import cupy as cp
 import pytest
 import qibo
 from qibo.models import QFT
@@ -46,3 +47,40 @@ def test_eval(nqubits: int, dtype="complex128"):
     assert 1e-2 * qibo_time < cutn_time < 1e2 * qibo_time
     assert np.allclose(
         result_sv, result_tn), "Resulting dense vectors do not match"
+
+
+@pytest.mark.gpu
+@pytest.mark.parametrize("nqubits", [2, 5, 10])
+def test_mps(nqubits: int, dtype="complex128"):
+    """Evaluate MPS with cuQuantum.
+
+    Args:
+        nqubits (int): Total number of qubits in the system.
+        dtype (str): The data type for precision, 'complex64' for single,
+            'complex128' for double.
+    """
+    import qibotn.cutn
+
+    # Test qibo
+    qibo.set_backend(backend=config.qibo.backend,
+                     platform=config.qibo.platform)
+
+    qibo_time, (circ_qibo, result_sv) = time(
+        lambda: qibo_qft(nqubits, swaps=True))
+
+    result_sv_cp = cp.asarray(result_sv)
+
+    # Test of MPS
+    gate_algo = {'qr_method': False,
+                 'svd_method': {
+                     'partition': 'UV',
+                     'abs_cutoff': 1e-12,
+                 }}
+
+    cutn_time, result_tn = time(
+        lambda: qibotn.cutn.eval_mps(circ_qibo, gate_algo, dtype).flatten())
+
+    print(
+        f"State vector difference: {abs(result_tn - result_sv_cp).max():0.3e}")
+
+    assert cp.allclose(result_tn, result_sv_cp)

From 4da70db97cdd8bb13e0428022af605fae2d55c9a Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Mon, 24 Jul 2023 17:55:39 +0800
Subject: [PATCH 06/31] remove MPO operation

---
 src/qibotn/MPSContractionHelper.py | 28 ----------------------------
 1 file changed, 28 deletions(-)

diff --git a/src/qibotn/MPSContractionHelper.py b/src/qibotn/MPSContractionHelper.py
index cd0905e7..60192989 100644
--- a/src/qibotn/MPSContractionHelper.py
+++ b/src/qibotn/MPSContractionHelper.py
@@ -117,34 +117,6 @@ def contract_expectation(self, mps_tensors, operator, qubits, options=None, norm
             norm = 1
         return self._contract(f'exp{qubits}', interleaved_inputs, options=options) / norm
     
-    def contract_mps_mpo_to_state_vector(self, mps_tensors, mpo_tensors, options=None):
-        """
-        Contract the corresponding tensor network to form the output state vector from applying the MPO to the MPS.
-
-        Args:
-            mps_tensors: A list of rank-3 ndarray-like tensor objects. 
-                The indices of the ith tensor are expected to be the bonding index to the i-1 tensor, 
-                the physical mode, and then the bonding index to the i+1th tensor.
-            mpo_tensors: A list of rank-4 ndarray-like tensor objects.
-                The indics of the ith tensor are expected to be the bonding index to the i-1 tensor, 
-                the output physical mode, the bonding index to the i+1th tensor and then the inputput physical mode.
-            options: Specify the contract and decompose options. 
-
-        Returns:
-            An ndarray-like object as the output state vector.
-        """
-        interleaved_inputs = []
-        for i, o in enumerate(mps_tensors):
-            interleaved_inputs.extend([o, self.bra_modes[i]])
-        output_modes = []
-        offset = 2 * self.num_qubits + 1
-        for i, o in enumerate(mpo_tensors):
-            mpo_modes = (2*i+offset, 2*i+offset+1, 2*i+offset+2, 2*i+1)
-            output_modes.append(2*i+offset+1)
-            interleaved_inputs.extend([o, mpo_modes])
-        interleaved_inputs.append(output_modes)
-        return self._contract('mps_mpo', interleaved_inputs, options=options)
-    
     def _contract(self, key, interleaved_inputs, options=None):
         """
         Perform the contraction task given interleaved inputs. Path will be cached.

From ac712d241e80b0758c30a46b7f2fadbdaa057607 Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Mon, 24 Jul 2023 18:04:59 +0800
Subject: [PATCH 07/31] removed additional helper file

---
 src/qibotn/MpsHelper.py | 270 ----------------------------------------
 1 file changed, 270 deletions(-)
 delete mode 100644 src/qibotn/MpsHelper.py

diff --git a/src/qibotn/MpsHelper.py b/src/qibotn/MpsHelper.py
deleted file mode 100644
index 526cd3c3..00000000
--- a/src/qibotn/MpsHelper.py
+++ /dev/null
@@ -1,270 +0,0 @@
-# Copyright (c) 2021-2023, NVIDIA CORPORATION & AFFILIATES
-#
-# SPDX-License-Identifier: BSD-3-Clause
-
-import itertools
-
-import cupy as cp
-import numpy as np
-
-import cuquantum
-from cuquantum import cutensornet as cutn
-
-class MPSHelper:
-    
-    """
-    MPSHelper(num_sites, phys_extent, max_virtual_extent, initial_state, data_type, compute_type)
-
-    Create an MPSHelper object for gate splitting algorithm.
-           i              j       
-        -------A-------B-------                      i        j        k
-              p|       |q            ------->     -------A`-------B`-------
-               GGGGGGGGG                                r|        |s
-              r|       |s
-    
-    Args:
-        num_sites: The number of sites in the MPS.
-        phys_extents: The extent for the physical mode where the gate tensors are acted on.
-        max_virtual_extent: The maximal extent allowed for the virtual mode shared between adjacent MPS tensors.
-        initial_state: A sequence of :class:`cupy.ndarray` representing the initial state of the MPS.
-        data_type (cuquantum.cudaDataType): The data type for all tensors and gates. 
-        compute_type (cuquantum.ComputeType): The compute type for all gate splitting.
-    
-    """
-
-    def __init__(self, num_sites, phys_extent, max_virtual_extent, initial_state, data_type, compute_type):
-        self.num_sites = num_sites
-        self.phys_extent = phys_extent
-        self.data_type = data_type
-        self.compute_type = compute_type
-        
-        self.phys_modes = []
-        self.virtual_modes = []
-        self.new_mode = itertools.count(start=0, step=1)
-        
-        for i in range(num_sites+1):
-            self.virtual_modes.append(next(self.new_mode))
-            if i != num_sites:
-                self.phys_modes.append(next(self.new_mode))
-        
-        untruncated_max_extent = phys_extent ** (num_sites // 2)
-        if max_virtual_extent == 0:
-            self.max_virtual_extent = untruncated_max_extent
-        else:
-            self.max_virtual_extent = min(max_virtual_extent, untruncated_max_extent)
-        
-        self.handle = cutn.create()
-        self.work_desc = cutn.create_workspace_descriptor(self.handle)
-        self.svd_config = cutn.create_tensor_svd_config(self.handle)
-        self.svd_info = cutn.create_tensor_svd_info(self.handle)
-        self.gate_algo = cutn.GateSplitAlgo.DIRECT
-
-        self.desc_tensors = []
-        self.state_tensors = []
-
-        # create tensor descriptors
-        for i in range(self.num_sites):
-            temp = initial_state[i]
-            self.state_tensors.append(initial_state[i].astype(temp.dtype, order="F"))
-            extent = self.get_tensor_extent(i)
-            modes = self.get_tensor_modes(i)
-            desc_tensor = cutn.create_tensor_descriptor(self.handle, 3, extent, 0, modes, self.data_type)
-            self.desc_tensors.append(desc_tensor)
-    
-    def get_tensor(self, site):
-        """Get the tensor operands for a specific site."""
-        return self.state_tensors[site]
-    
-    def get_tensor_extent(self, site):
-        """Get the extent of the MPS tensor at a specific site."""
-        return self.state_tensors[site].shape
-    
-    def get_tensor_modes(self, site):
-        """Get the current modes of the MPS tensor at a specific site."""
-        return (self.virtual_modes[site], self.phys_modes[site], self.virtual_modes[site+1])
-
-    def set_svd_config(self, abs_cutoff, rel_cutoff, renorm, partition):
-        """Update the SVD truncation setting.
-        
-        Args:
-            abs_cutoff: The cutoff value for absolute singular value truncation.
-            rel_cutoff: The cutoff value for relative singular value truncation.
-            renorm (cuquantum.cutensornet.TensorSVDNormalization): The option for renormalization of the truncated singular values.
-            partition (cuquantum.cutensornet.TensorSVDPartition): The option for partitioning of the singular values.
-        """        
-        
-        if partition != cutn.TensorSVDPartition.UV_EQUAL:
-            raise NotImplementedError("this basic example expects partition to be cutensornet.TensorSVDPartition.UV_EQUAL")
-
-        svd_config_attributes = [cutn.TensorSVDConfigAttribute.ABS_CUTOFF, 
-                                 cutn.TensorSVDConfigAttribute.REL_CUTOFF, 
-                                 cutn.TensorSVDConfigAttribute.S_NORMALIZATION,
-                                 cutn.TensorSVDConfigAttribute.S_PARTITION]
-            
-        for (attr, value) in zip(svd_config_attributes, [abs_cutoff, rel_cutoff, renorm, partition]):
-            dtype = cutn.tensor_svd_config_get_attribute_dtype(attr)
-            value = np.array([value], dtype=dtype)
-            cutn.tensor_svd_config_set_attribute(self.handle, 
-                self.svd_config, attr, value.ctypes.data, value.dtype.itemsize)
-    
-    def set_gate_algorithm(self, gate_algo):
-        """Set the algorithm to use for all gate split operations.
-        
-        Args:
-            gate_algo (cuquantum.cutensornet.GateSplitAlgo): The gate splitting algorithm to use.
-        """
-
-        self.gate_algo = gate_algo
-
-    def compute_max_workspace_sizes(self):
-        """Compute the maximal workspace needed for MPS gating algorithm."""
-        modes_in_A = [ord(c) for c in ('i', 'p', 'j')]
-        modes_in_B = [ord(c) for c in ('j', 'q', 'k')]
-        modes_in_G = [ord(c) for c in ('p', 'q', 'r', 's')]
-        modes_out_A = [ord(c) for c in ('i', 'r', 'j')]
-        modes_out_B = [ord(c) for c in ('j', 's', 'k')]
-        
-        max_extents_AB = (self.max_virtual_extent, self.phys_extent, self.max_virtual_extent)
-        extents_in_G = (self.phys_extent, self.phys_extent, self.phys_extent, self.phys_extent)
-        
-        desc_tensor_in_A = cutn.create_tensor_descriptor(self.handle, 3, max_extents_AB, 0, modes_in_A, self.data_type)
-        desc_tensor_in_B = cutn.create_tensor_descriptor(self.handle, 3, max_extents_AB, 0, modes_in_B, self.data_type)
-        desc_tensor_in_G = cutn.create_tensor_descriptor(self.handle, 4, extents_in_G, 0, modes_in_G, self.data_type)
-        desc_tensor_out_A = cutn.create_tensor_descriptor(self.handle, 3, max_extents_AB, 0, modes_out_A, self.data_type)
-        desc_tensor_out_B = cutn.create_tensor_descriptor(self.handle, 3, max_extents_AB, 0, modes_out_B, self.data_type)
-        
-        cutn.workspace_compute_gate_split_sizes(self.handle, 
-            desc_tensor_in_A, desc_tensor_in_B, desc_tensor_in_G, 
-            desc_tensor_out_A, desc_tensor_out_B, 
-            self.gate_algo, self.svd_config, self.compute_type, self.work_desc)
-        
-        workspace_size = cutn.workspace_get_memory_size(self.handle, self.work_desc, cutn.WorksizePref.MIN, cutn.Memspace.DEVICE, cutn.WorkspaceKind.SCRATCH)
-
-        # free resources
-        cutn.destroy_tensor_descriptor(desc_tensor_in_A)
-        cutn.destroy_tensor_descriptor(desc_tensor_in_B)
-        cutn.destroy_tensor_descriptor(desc_tensor_in_G)
-        cutn.destroy_tensor_descriptor(desc_tensor_out_A)
-        cutn.destroy_tensor_descriptor(desc_tensor_out_B)
-        return workspace_size
-
-    def set_workspace(self, work, workspace_size):
-        """Compute the maximal workspace needed for MPS gating algorithm.
-        
-        Args:
-            work: Pointer to the allocated workspace.
-            workspace_size: The required workspace size on the device.
-        """
-        cutn.workspace_set_memory(self.handle, self.work_desc, cutn.Memspace.DEVICE, cutn.WorkspaceKind.SCRATCH, work.ptr, workspace_size)
-
-    def apply_gate(self, site_A, site_B, gate, verbose, stream):
-        """Inplace execution of the apply gate algoritm on site A and site B.
-        
-        Args:
-            site_A: The first site on which the gate is applied to.
-            site_B: The second site on which the gate is applied to.
-            gate (cupy.ndarray): The input data for the gate tensor.
-            verbose: Whether to print out the runtime information during truncation.
-            stream (cupy.cuda.Stream): The CUDA stream on which the computation is performed.
-        """
-        if site_B - site_A != 1:
-            raise ValueError("Site B must be the right site of site A")
-        if site_B >= self.num_sites:
-            raise ValueError("Site index cannot exceed maximum number of sites")
-        
-        desc_tensor_in_A = self.desc_tensors[site_A]
-        desc_tensor_in_B = self.desc_tensors[site_B]
-        
-        phys_mode_in_A = self.phys_modes[site_A]
-        phys_mode_in_B = self.phys_modes[site_B]
-        phys_mode_out_A = next(self.new_mode)
-        phys_mode_out_B = next(self.new_mode)
-        modes_G = (phys_mode_in_A, phys_mode_in_B, phys_mode_out_A, phys_mode_out_B)
-        extent_G = (self.phys_extent, self.phys_extent, self.phys_extent, self.phys_extent)
-        desc_tensor_in_G = cutn.create_tensor_descriptor(self.handle, 4, extent_G, 0, modes_G, self.data_type)
-
-        # construct and initialize the expected output A and B
-        tensor_in_A = self.state_tensors[site_A]
-        tensor_in_B = self.state_tensors[site_B]
-        left_extent_A = tensor_in_A.shape[0]
-        extent_AB_in = tensor_in_A.shape[2]
-        right_extent_B = tensor_in_B.shape[2]
-        combined_extent_left = min(left_extent_A, extent_AB_in * self.phys_extent) * self.phys_extent
-        combined_extent_right = min(right_extent_B, extent_AB_in * self.phys_extent) * self.phys_extent
-        extent_Aout_B = min(combined_extent_left, combined_extent_right, self.max_virtual_extent)
-        
-        extent_out_A = (left_extent_A, self.phys_extent, extent_Aout_B)
-        extent_out_B = (extent_Aout_B, self.phys_extent, right_extent_B)
-
-        tensor_out_A = cp.zeros(extent_out_A, dtype=tensor_in_A.dtype, order="F")
-        tensor_out_B = cp.zeros(extent_out_B, dtype=tensor_in_B.dtype, order="F")
-
-        # create tensor descriptors for output A and B
-        modes_out_A = (self.virtual_modes[site_A], phys_mode_out_A, self.virtual_modes[site_A+1])
-        modes_out_B = (self.virtual_modes[site_B], phys_mode_out_B, self.virtual_modes[site_B+1])
-
-        desc_tensor_out_A = cutn.create_tensor_descriptor(self.handle, 3, extent_out_A, 0, modes_out_A, self.data_type)
-        desc_tensor_out_B = cutn.create_tensor_descriptor(self.handle, 3, extent_out_B, 0, modes_out_B, self.data_type)
-
-        cutn.gate_split(self.handle, 
-                        desc_tensor_in_A, tensor_in_A.data.ptr, 
-                        desc_tensor_in_B, tensor_in_B.data.ptr, 
-                        desc_tensor_in_G, gate.data.ptr, 
-                        desc_tensor_out_A, tensor_out_A.data.ptr, 
-                        0, # we factorize singular values equally onto output A and B.
-                        desc_tensor_out_B, tensor_out_B.data.ptr, 
-                        self.gate_algo, self.svd_config, self.compute_type, 
-                        self.svd_info, self.work_desc, stream.ptr)
-        
-        if verbose:
-            full_extent = np.array([0], dtype=cutn.tensor_svd_info_get_attribute_dtype(cutn.TensorSVDInfoAttribute.FULL_EXTENT))
-            reduced_extent = np.array([0], dtype=cutn.tensor_svd_info_get_attribute_dtype(cutn.TensorSVDInfoAttribute.REDUCED_EXTENT))
-            discarded_weight = np.array([0], dtype=cutn.tensor_svd_info_get_attribute_dtype(cutn.TensorSVDInfoAttribute.DISCARDED_WEIGHT))
-            
-            cutn.tensor_svd_info_get_attribute(
-                self.handle, self.svd_info, cutn.TensorSVDInfoAttribute.FULL_EXTENT, 
-                full_extent.ctypes.data, full_extent.dtype.itemsize)
-            cutn.tensor_svd_info_get_attribute(
-                self.handle, self.svd_info, cutn.TensorSVDInfoAttribute.REDUCED_EXTENT, 
-                reduced_extent.ctypes.data, reduced_extent.dtype.itemsize)
-            cutn.tensor_svd_info_get_attribute(
-                self.handle, self.svd_info, cutn.TensorSVDInfoAttribute.DISCARDED_WEIGHT, 
-                discarded_weight.ctypes.data, discarded_weight.dtype.itemsize)
-            
-            print("Virtual bond truncated from {0} to {1} with a discarded weight of {2:.6f}".format(full_extent[0], reduced_extent[0], discarded_weight[0]))
-
-        self.phys_modes[site_A] = phys_mode_out_A
-        self.phys_modes[site_B] = phys_mode_out_B
-        self.desc_tensors[site_A] = desc_tensor_out_A
-        self.desc_tensors[site_B] = desc_tensor_out_B
-        
-        extent_out_A = np.zeros((3,), dtype=np.int64)
-        extent_out_B = np.zeros((3,), dtype=np.int64)
-        extent_out_A, strides_out_A = cutn.get_tensor_details(self.handle, desc_tensor_out_A)[2:]
-        extent_out_B, strides_out_B = cutn.get_tensor_details(self.handle, desc_tensor_out_B)[2:]
-
-        # Recall that `cutensornet.gate_split` can potentially find reduced extent during SVD truncation when value-based truncation is used.
-        # Therefore we here update the container for output tensor A and B.
-        if extent_out_A[2] != extent_Aout_B:
-            # note strides in cutensornet are in the unit of count and strides in cupy/numpy are in the unit of nbytes
-            strides_out_A = [i * tensor_out_A.itemsize for i in strides_out_A]
-            strides_out_B = [i * tensor_out_B.itemsize for i in strides_out_B]
-            tensor_out_A = cp.ndarray(extent_out_A, dtype=tensor_out_A.dtype, memptr=tensor_out_A.data, strides=strides_out_A)
-            tensor_out_B = cp.ndarray(extent_out_B, dtype=tensor_out_B.dtype, memptr=tensor_out_B.data, strides=strides_out_B)
-
-        self.state_tensors[site_A] = tensor_out_A
-        self.state_tensors[site_B] = tensor_out_B
-
-        cutn.destroy_tensor_descriptor(desc_tensor_in_A)
-        cutn.destroy_tensor_descriptor(desc_tensor_in_B)
-        cutn.destroy_tensor_descriptor(desc_tensor_in_G)
-    
-    def __del__(self):
-        """Free all resources owned by the object."""
-        for desc_tensor in self.desc_tensors:
-            cutn.destroy_tensor_descriptor(desc_tensor)
-        cutn.destroy(self.handle)
-        cutn.destroy_workspace_descriptor(self.work_desc)
-        cutn.destroy_tensor_svd_config(self.svd_config)
-        cutn.destroy_tensor_svd_info(self.svd_info)
-

From 6a22db07ad96cea5a95c049d7a499c48631b8837 Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Mon, 24 Jul 2023 18:06:19 +0800
Subject: [PATCH 08/31] remove __main__

---
 src/qibotn/cutn.py | 16 ----------------
 1 file changed, 16 deletions(-)

diff --git a/src/qibotn/cutn.py b/src/qibotn/cutn.py
index 21c316fc..222ebc9e 100644
--- a/src/qibotn/cutn.py
+++ b/src/qibotn/cutn.py
@@ -19,19 +19,3 @@ def eval_mps(qibo_circ, gate_algo, datatype):
     sv_mps = mps_helper.contract_state_vector(
         myconvertor.mps_tensors, myconvertor.options)
     return sv_mps
-
-
-if __name__ == "__main__":
-    num_qubits = 25
-    swaps = True
-    circ_qibo = QFT(num_qubits, swaps)
-
-    exact_gate_algorithm = {'qr_method': False,
-                            'svd_method': {'partition': 'UV', 'abs_cutoff': 1e-12}}
-    dtype = "complex128"
-    sv_mps = eval_mps(circ_qibo, exact_gate_algorithm, dtype)
-    sv_reference = eval(circ_qibo, dtype)
-    state_vec = np.array(circ_qibo())
-    print(f"State vector difference: {abs(sv_mps-sv_reference).max():0.3e}")
-    assert cp.allclose(sv_mps, sv_reference)
-    assert cp.allclose(sv_mps.flatten(), state_vec)

From 73d65a6af1cf73aacd612e1bc053a4435bd84467 Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Thu, 3 Aug 2023 10:23:56 +0800
Subject: [PATCH 09/31] Removed unused imports

---
 src/qibotn/cutn.py | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/src/qibotn/cutn.py b/src/qibotn/cutn.py
index 222ebc9e..9b9ac7a1 100644
--- a/src/qibotn/cutn.py
+++ b/src/qibotn/cutn.py
@@ -1,9 +1,6 @@
 from qibotn.QiboCircuitConvertor import QiboCircuitToEinsum
 from cuquantum import contract
-from cuquantum import cutensornet as cutn
-import cupy as cp
-import numpy as np
-from qibo.models import QFT
+
 from qibotn.QiboCircuitToMPS import QiboCircuitToMPS
 from qibotn.MPSContractionHelper import MPSContractionHelper
 

From fc3e0c24bdd1f51769d1cb5b792ad4505a5a8bf6 Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Thu, 3 Aug 2023 10:27:19 +0800
Subject: [PATCH 10/31] Return function output directly

---
 src/qibotn/cutn.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/qibotn/cutn.py b/src/qibotn/cutn.py
index 9b9ac7a1..fc603bac 100644
--- a/src/qibotn/cutn.py
+++ b/src/qibotn/cutn.py
@@ -13,6 +13,7 @@ def eval(qibo_circ, datatype):
 def eval_mps(qibo_circ, gate_algo, datatype):
     myconvertor = QiboCircuitToMPS(qibo_circ, gate_algo, dtype=datatype)
     mps_helper = MPSContractionHelper(myconvertor.num_qubits)
-    sv_mps = mps_helper.contract_state_vector(
-        myconvertor.mps_tensors, myconvertor.options)
-    return sv_mps
+
+    return mps_helper.contract_state_vector(
+        myconvertor.mps_tensors, myconvertor.options
+    )

From ca6d5796dcba7dd4cc4b5e4f608678062b4b2fd3 Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Thu, 3 Aug 2023 11:26:20 +0800
Subject: [PATCH 11/31] Removed duplicated self.options

---
 src/qibotn/QiboCircuitToMPS.py | 20 +++++++++++++++-----
 1 file changed, 15 insertions(+), 5 deletions(-)

diff --git a/src/qibotn/QiboCircuitToMPS.py b/src/qibotn/QiboCircuitToMPS.py
index 4c18aac0..f504d5ad 100644
--- a/src/qibotn/QiboCircuitToMPS.py
+++ b/src/qibotn/QiboCircuitToMPS.py
@@ -7,22 +7,32 @@
 
 
 class QiboCircuitToMPS:
-    def __init__(self, circ_qibo, gate_algo, dtype='complex128', rand_seed=0,):
+    def __init__(
+        self,
+        circ_qibo,
+        gate_algo,
+        dtype="complex128",
+        rand_seed=0,
+    ):
         np.random.seed(rand_seed)
         cp.random.seed(rand_seed)
 
         self.num_qubits = circ_qibo.nqubits
         self.handle = cutn.create()
-        self.options = {'handle': self.handle}
         self.dtype = dtype
         self.mps_tensors = get_initial_mps(self.num_qubits, dtype=dtype)
         circuitconvertor = QiboCircuitToEinsum(circ_qibo)
 
-        for (gate, qubits) in circuitconvertor.gate_tensors:
+        for gate, qubits in circuitconvertor.gate_tensors:
             # mapping from qubits to qubit indices
             # apply the gate in-place
-            apply_gate(self.mps_tensors, gate, qubits,
-                       algorithm=gate_algo, options=self.options)
+            apply_gate(
+                self.mps_tensors,
+                gate,
+                qubits,
+                algorithm=gate_algo,
+                options={"handle": self.handle},
+            )
 
     def __del__(self):
         cutn.destroy(self.handle)

From ce018abb16b8053f155381cce6330f3a31970add Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Thu, 3 Aug 2023 13:32:20 +0800
Subject: [PATCH 12/31] Removed default dtype

---
 src/qibotn/MPSUtils.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/qibotn/MPSUtils.py b/src/qibotn/MPSUtils.py
index cb77e4b9..b5f472ab 100644
--- a/src/qibotn/MPSUtils.py
+++ b/src/qibotn/MPSUtils.py
@@ -2,7 +2,7 @@
 from cuquantum.cutensornet.experimental import contract_decompose
 from cuquantum import contract
 
-def get_initial_mps(num_qubits, dtype='complex128'):
+def initial(num_qubits, dtype):
     """
     Generate the MPS with an initial state of |00...00> 
     """
@@ -75,4 +75,4 @@ def apply_gate(
             mps_site_right_swap(mps_tensors, i, algorithm=algorithm, options=options)
     else:
         raise NotImplementedError("Only one- and two-qubit gates supported")
-    return mps_tensors
\ No newline at end of file
+    return mps_tensors

From 32611cdd74d9064df2add1eea84e76bf1fe3aa97 Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Thu, 3 Aug 2023 13:34:51 +0800
Subject: [PATCH 13/31] Update the mps initialization function name

---
 src/qibotn/QiboCircuitToMPS.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/qibotn/QiboCircuitToMPS.py b/src/qibotn/QiboCircuitToMPS.py
index f504d5ad..d51093f5 100644
--- a/src/qibotn/QiboCircuitToMPS.py
+++ b/src/qibotn/QiboCircuitToMPS.py
@@ -3,7 +3,7 @@
 
 from cuquantum import cutensornet as cutn
 from qibotn.QiboCircuitConvertor import QiboCircuitToEinsum
-from qibotn.MPSUtils import get_initial_mps, apply_gate
+from qibotn.MPSUtils import initial, apply_gate
 
 
 class QiboCircuitToMPS:
@@ -20,7 +20,7 @@ def __init__(
         self.num_qubits = circ_qibo.nqubits
         self.handle = cutn.create()
         self.dtype = dtype
-        self.mps_tensors = get_initial_mps(self.num_qubits, dtype=dtype)
+        self.mps_tensors = initial(self.num_qubits, dtype=dtype)
         circuitconvertor = QiboCircuitToEinsum(circ_qibo)
 
         for gate, qubits in circuitconvertor.gate_tensors:

From c15e4ca9677b0f88f7454e845236303487e5aa5e Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Thu, 3 Aug 2023 15:32:51 +0800
Subject: [PATCH 14/31] Update function input as options attribute removed

---
 src/qibotn/cutn.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/qibotn/cutn.py b/src/qibotn/cutn.py
index fc603bac..ab0e9c12 100644
--- a/src/qibotn/cutn.py
+++ b/src/qibotn/cutn.py
@@ -15,5 +15,5 @@ def eval_mps(qibo_circ, gate_algo, datatype):
     mps_helper = MPSContractionHelper(myconvertor.num_qubits)
 
     return mps_helper.contract_state_vector(
-        myconvertor.mps_tensors, myconvertor.options
+        myconvertor.mps_tensors, {"handle": myconvertor.handle}
     )

From 32bee13044da17305708d6f2f0d7b625de940907 Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Thu, 3 Aug 2023 15:33:34 +0800
Subject: [PATCH 15/31] Removed return of tensor as it is modified in place

---
 src/qibotn/MPSUtils.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/qibotn/MPSUtils.py b/src/qibotn/MPSUtils.py
index b5f472ab..1a9dad5d 100644
--- a/src/qibotn/MPSUtils.py
+++ b/src/qibotn/MPSUtils.py
@@ -75,4 +75,4 @@ def apply_gate(
             mps_site_right_swap(mps_tensors, i, algorithm=algorithm, options=options)
     else:
         raise NotImplementedError("Only one- and two-qubit gates supported")
-    return mps_tensors
+    #return mps_tensors

From 139748dcefbaae163fddcd620ae1564baf40b962 Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Tue, 15 Aug 2023 14:29:09 +0800
Subject: [PATCH 16/31] change to literal for dict

---
 src/qibotn/MPSContractionHelper.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/qibotn/MPSContractionHelper.py b/src/qibotn/MPSContractionHelper.py
index 60192989..fe2e98e3 100644
--- a/src/qibotn/MPSContractionHelper.py
+++ b/src/qibotn/MPSContractionHelper.py
@@ -33,7 +33,7 @@ class MPSContractionHelper:
     
     def __init__(self, num_qubits):
         self.num_qubits = num_qubits
-        self.path_cache = dict()
+        self.path_cache = {}
         self.bra_modes = [(2*i, 2*i+1, 2*i+2) for i in range(num_qubits)]
         offset = 2*num_qubits+1
         self.ket_modes = [(i+offset, 2*i+1, i+1+offset) for i in range(num_qubits)]

From 696f0524cc14fda8feedb08a24b0b1b18e745219 Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Tue, 15 Aug 2023 16:03:09 +0800
Subject: [PATCH 17/31] change list conversion method

---
 src/qibotn/MPSContractionHelper.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/qibotn/MPSContractionHelper.py b/src/qibotn/MPSContractionHelper.py
index fe2e98e3..3bf11062 100644
--- a/src/qibotn/MPSContractionHelper.py
+++ b/src/qibotn/MPSContractionHelper.py
@@ -99,7 +99,8 @@ def contract_expectation(self, mps_tensors, operator, qubits, options=None, norm
         interleaved_inputs = []
         extra_mode = 3 * self.num_qubits + 2
         operator_modes = [None] * len(qubits) + [self.bra_modes[q][1] for q in qubits]
-        qubits = list(qubits)
+        #qubits = list(qubits)
+        qubits = [q for q in qubits]
         for i, o in enumerate(mps_tensors):
             interleaved_inputs.extend([o, self.bra_modes[i]])
             k_modes = self.ket_modes[i]

From 7abb82c1fb91109a8b5d5f1e28aa639276e0b875 Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Tue, 15 Aug 2023 16:03:39 +0800
Subject: [PATCH 18/31] Change list conversion method

---
 src/qibotn/MPSContractionHelper.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/qibotn/MPSContractionHelper.py b/src/qibotn/MPSContractionHelper.py
index 3bf11062..b4ac3660 100644
--- a/src/qibotn/MPSContractionHelper.py
+++ b/src/qibotn/MPSContractionHelper.py
@@ -99,7 +99,6 @@ def contract_expectation(self, mps_tensors, operator, qubits, options=None, norm
         interleaved_inputs = []
         extra_mode = 3 * self.num_qubits + 2
         operator_modes = [None] * len(qubits) + [self.bra_modes[q][1] for q in qubits]
-        #qubits = list(qubits)
         qubits = [q for q in qubits]
         for i, o in enumerate(mps_tensors):
             interleaved_inputs.extend([o, self.bra_modes[i]])

From a42ba15e88f47ab89e90e110a48942bd3db60d4f Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Wed, 16 Aug 2023 18:18:27 +0800
Subject: [PATCH 19/31] use **kwargs

---
 src/qibotn/MPSUtils.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/qibotn/MPSUtils.py b/src/qibotn/MPSUtils.py
index 1a9dad5d..8b8cdb3f 100644
--- a/src/qibotn/MPSUtils.py
+++ b/src/qibotn/MPSUtils.py
@@ -13,14 +13,13 @@ def initial(num_qubits, dtype):
 def mps_site_right_swap(
     mps_tensors, 
     i, 
-    algorithm=None, 
-    options=None
+    **kwargs
 ):
     """
     Perform the swap operation between the ith and i+1th MPS tensors.
     """
     # contraction followed by QR decomposition
-    a, _, b = contract_decompose('ipj,jqk->iqj,jpk', *mps_tensors[i:i+2], algorithm=algorithm, options=options)
+    a, _, b = contract_decompose('ipj,jqk->iqj,jpk', *mps_tensors[i:i+2], algorithm=list(kwargs.items())[0][1], options=list(kwargs.items())[1][1])
     mps_tensors[i:i+2] = (a, b)
     return mps_tensors
 

From c3809d36e66e84e4ff4da800d97212323cc3c983 Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Wed, 16 Aug 2023 18:19:09 +0800
Subject: [PATCH 20/31] remove return in apply_gate()

---
 src/qibotn/MPSUtils.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/qibotn/MPSUtils.py b/src/qibotn/MPSUtils.py
index 8b8cdb3f..767ece6f 100644
--- a/src/qibotn/MPSUtils.py
+++ b/src/qibotn/MPSUtils.py
@@ -74,4 +74,3 @@ def apply_gate(
             mps_site_right_swap(mps_tensors, i, algorithm=algorithm, options=options)
     else:
         raise NotImplementedError("Only one- and two-qubit gates supported")
-    #return mps_tensors

From cb21d1d1c0b080de2483c720eeba304aa0b21199 Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Thu, 17 Aug 2023 10:41:32 +0800
Subject: [PATCH 21/31] Update fucntions to use **kwargs

---
 src/qibotn/MPSUtils.py | 17 ++++++++---------
 1 file changed, 8 insertions(+), 9 deletions(-)

diff --git a/src/qibotn/MPSUtils.py b/src/qibotn/MPSUtils.py
index 767ece6f..0ad53e16 100644
--- a/src/qibotn/MPSUtils.py
+++ b/src/qibotn/MPSUtils.py
@@ -19,7 +19,7 @@ def mps_site_right_swap(
     Perform the swap operation between the ith and i+1th MPS tensors.
     """
     # contraction followed by QR decomposition
-    a, _, b = contract_decompose('ipj,jqk->iqj,jpk', *mps_tensors[i:i+2], algorithm=list(kwargs.items())[0][1], options=list(kwargs.items())[1][1])
+    a, _, b = contract_decompose('ipj,jqk->iqj,jpk', *mps_tensors[i:i+2], algorithm=kwargs.get('algorithm',None), options=kwargs.get('options',None))
     mps_tensors[i:i+2] = (a, b)
     return mps_tensors
 
@@ -27,8 +27,7 @@ def apply_gate(
     mps_tensors, 
     gate, 
     qubits, 
-    algorithm=None, 
-    options=None
+    **kwargs
 ):
     """
     Apply the gate operand to the MPS tensors in-place.
@@ -53,24 +52,24 @@ def apply_gate(
     if n_qubits == 1:
         # single-qubit gate
         i = qubits[0]
-        mps_tensors[i] = contract('ipj,qp->iqj', mps_tensors[i], gate, options=options) # in-place update
+        mps_tensors[i] = contract('ipj,qp->iqj', mps_tensors[i], gate, options=kwargs.get('options',None)) # in-place update
     elif n_qubits == 2:
         # two-qubit gate
         i, j = qubits
         if i > j:
             # swap qubits order
-            return apply_gate(mps_tensors, gate.transpose(1,0,3,2), (j, i), algorithm=algorithm, options=options)
+            return apply_gate(mps_tensors, gate.transpose(1,0,3,2), (j, i), **kwargs)
         elif i+1 == j:
             # two adjacent qubits
-            a, _, b = contract_decompose('ipj,jqk,rspq->irj,jsk', *mps_tensors[i:i+2], gate, algorithm=algorithm, options=options)
+            a, _, b = contract_decompose('ipj,jqk,rspq->irj,jsk', *mps_tensors[i:i+2], gate, algorithm=kwargs.get('algorithm',None), options=kwargs.get('options',None))
             mps_tensors[i:i+2] = (a, b) # in-place update
         else:
             # non-adjacent two-qubit gate
             # step 1: swap i with i+1
-            mps_site_right_swap(mps_tensors, i, algorithm=algorithm, options=options)
+            mps_site_right_swap(mps_tensors, i, **kwargs)
             # step 2: apply gate to (i+1, j) pair. This amounts to a recursive swap until the two qubits are adjacent
-            apply_gate(mps_tensors, gate, (i+1, j), algorithm=algorithm, options=options) 
+            apply_gate(mps_tensors, gate, (i+1, j), **kwargs) 
             # step 3: swap back i and i+1
-            mps_site_right_swap(mps_tensors, i, algorithm=algorithm, options=options)
+            mps_site_right_swap(mps_tensors, i, **kwargs)
     else:
         raise NotImplementedError("Only one- and two-qubit gates supported")

From 3fafe2b3ff120f5ba6e3155f0f2ff8575da74ee0 Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Thu, 17 Aug 2023 10:53:56 +0800
Subject: [PATCH 22/31] Remove path cache

---
 src/qibotn/MPSContractionHelper.py | 18 +++++++-----------
 1 file changed, 7 insertions(+), 11 deletions(-)

diff --git a/src/qibotn/MPSContractionHelper.py b/src/qibotn/MPSContractionHelper.py
index b4ac3660..0dd6007f 100644
--- a/src/qibotn/MPSContractionHelper.py
+++ b/src/qibotn/MPSContractionHelper.py
@@ -33,7 +33,6 @@ class MPSContractionHelper:
     
     def __init__(self, num_qubits):
         self.num_qubits = num_qubits
-        self.path_cache = {}
         self.bra_modes = [(2*i, 2*i+1, 2*i+2) for i in range(num_qubits)]
         offset = 2*num_qubits+1
         self.ket_modes = [(i+offset, 2*i+1, i+1+offset) for i in range(num_qubits)]
@@ -55,7 +54,7 @@ def contract_norm(self, mps_tensors, options=None):
         for i, o in enumerate(mps_tensors):
             interleaved_inputs.extend([o, self.bra_modes[i], o.conj(), self.ket_modes[i]])
         interleaved_inputs.append([]) # output
-        return self._contract('norm', interleaved_inputs, options=options).real
+        return self._contract(interleaved_inputs, options=options).real
 
     def contract_state_vector(self, mps_tensors, options=None):
         """
@@ -75,7 +74,7 @@ def contract_state_vector(self, mps_tensors, options=None):
             interleaved_inputs.extend([o, self.bra_modes[i]])
         output_modes = tuple([bra_modes[1] for bra_modes in self.bra_modes])
         interleaved_inputs.append(output_modes) # output
-        return self._contract('sv', interleaved_inputs, options=options)
+        return self._contract(interleaved_inputs, options=options)
     
     def contract_expectation(self, mps_tensors, operator, qubits, options=None, normalize=False):
         """
@@ -115,13 +114,10 @@ def contract_expectation(self, mps_tensors, operator, qubits, options=None, norm
             norm = self.contract_norm(mps_tensors, options=options)
         else:
             norm = 1
-        return self._contract(f'exp{qubits}', interleaved_inputs, options=options) / norm
+        return self._contract(interleaved_inputs, options=options) / norm
     
-    def _contract(self, key, interleaved_inputs, options=None):
-        """
-        Perform the contraction task given interleaved inputs. Path will be cached.
-        """
-        if key not in self.path_cache:
-            self.path_cache[key] = contract_path(*interleaved_inputs, options=options)[0]
-        path = self.path_cache[key]
+    def _contract(self, interleaved_inputs, options=None):
+
+        path = contract_path(*interleaved_inputs, options=options)[0]
+        
         return contract(*interleaved_inputs, options=options, optimize={'path':path})

From 89bdbfbe68a98003338312d92c713d9b818f7ec5 Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Thu, 17 Aug 2023 13:23:29 +0800
Subject: [PATCH 23/31] Format with Black

---
 src/qibotn/MPSContractionHelper.py | 94 ++++++++++++++++--------------
 src/qibotn/MPSUtils.py             | 69 ++++++++++++----------
 src/qibotn/QiboCircuitConvertor.py |  3 +-
 3 files changed, 89 insertions(+), 77 deletions(-)

diff --git a/src/qibotn/MPSContractionHelper.py b/src/qibotn/MPSContractionHelper.py
index 0dd6007f..9c71200a 100644
--- a/src/qibotn/MPSContractionHelper.py
+++ b/src/qibotn/MPSContractionHelper.py
@@ -1,59 +1,64 @@
 from cuquantum import contract, contract_path, CircuitToEinsum, tensor
 
+
 class MPSContractionHelper:
     """
     A helper class to compute various quantities for a given MPS.
-    
-    Interleaved format is used to construct the input args for `cuquantum.contract`. 
+
+    Interleaved format is used to construct the input args for `cuquantum.contract`.
     A concrete example on how the modes are populated for a 7-site MPS is provided below:
-    
-          0     2     4     6     8    10     12    14        
+
+          0     2     4     6     8    10     12    14
     bra -----A-----B-----C-----D-----E-----F-----G-----
-             |     |     |     |     |     |     |     
-            1|    3|    5|    7|    9|   11|   13|     
-             |     |     |     |     |     |     |     
+             |     |     |     |     |     |     |
+            1|    3|    5|    7|    9|   11|   13|
+             |     |     |     |     |     |     |
     ket -----a-----b-----c-----d-----e-----f-----g-----
           15    16    17    18    19    20    21    22
-    
-    
+
+
     The follwing compute quantities are supported:
-    
+
         - the norm of the MPS.
         - the equivalent state vector from the MPS.
         - the expectation value for a given operator.
         - the equivalent state vector after multiplying an MPO to an MPS.
-    
-    Note that for the nth MPS tensor (rank-3), the modes of the tensor are expected to be `(i,p,j)` 
-    where i denotes the bonding mode with the (n-1)th tensor, p denotes the physical mode for the qubit and 
+
+    Note that for the nth MPS tensor (rank-3), the modes of the tensor are expected to be `(i,p,j)`
+    where i denotes the bonding mode with the (n-1)th tensor, p denotes the physical mode for the qubit and
     j denotes the bonding mode with the (n+1)th tensor.
-    
+
     Args:
         num_qubits: The number of qubits for the MPS.
     """
-    
+
     def __init__(self, num_qubits):
         self.num_qubits = num_qubits
-        self.bra_modes = [(2*i, 2*i+1, 2*i+2) for i in range(num_qubits)]
-        offset = 2*num_qubits+1
-        self.ket_modes = [(i+offset, 2*i+1, i+1+offset) for i in range(num_qubits)]
-    
+        self.bra_modes = [(2 * i, 2 * i + 1, 2 * i + 2) for i in range(num_qubits)]
+        offset = 2 * num_qubits + 1
+        self.ket_modes = [
+            (i + offset, 2 * i + 1, i + 1 + offset) for i in range(num_qubits)
+        ]
+
     def contract_norm(self, mps_tensors, options=None):
         """
         Contract the corresponding tensor network to form the norm of the MPS.
 
         Args:
-            mps_tensors: A list of rank-3 ndarray-like tensor objects. 
-                The indices of the ith tensor are expected to be bonding index to the i-1 tensor, 
+            mps_tensors: A list of rank-3 ndarray-like tensor objects.
+                The indices of the ith tensor are expected to be bonding index to the i-1 tensor,
                 the physical mode, and then the bonding index to the i+1th tensor.
-            options: Specify the contract and decompose options. 
+            options: Specify the contract and decompose options.
 
         Returns:
             The norm of the MPS.
         """
         interleaved_inputs = []
         for i, o in enumerate(mps_tensors):
-            interleaved_inputs.extend([o, self.bra_modes[i], o.conj(), self.ket_modes[i]])
-        interleaved_inputs.append([]) # output
+            interleaved_inputs.extend(
+                [o, self.bra_modes[i], o.conj(), self.ket_modes[i]]
+            )
+        interleaved_inputs.append([])  # output
         return self._contract(interleaved_inputs, options=options).real
 
     def contract_state_vector(self, mps_tensors, options=None):
@@ -61,10 +66,10 @@ def contract_state_vector(self, mps_tensors, options=None):
         Contract the corresponding tensor network to form the state vector representation of the MPS.
 
         Args:
-            mps_tensors: A list of rank-3 ndarray-like tensor objects. 
-                The indices of the ith tensor are expected to be bonding index to the i-1 tensor, 
+            mps_tensors: A list of rank-3 ndarray-like tensor objects.
+                The indices of the ith tensor are expected to be bonding index to the i-1 tensor,
                 the physical mode, and then the bonding index to the i+1th tensor.
-            options: Specify the contract and decompose options. 
+            options: Specify the contract and decompose options.
 
         Returns:
             An ndarray-like object as the state vector.
@@ -73,28 +78,30 @@ def contract_state_vector(self, mps_tensors, options=None):
         for i, o in enumerate(mps_tensors):
             interleaved_inputs.extend([o, self.bra_modes[i]])
         output_modes = tuple([bra_modes[1] for bra_modes in self.bra_modes])
-        interleaved_inputs.append(output_modes) # output
+        interleaved_inputs.append(output_modes)  # output
         return self._contract(interleaved_inputs, options=options)
-    
-    def contract_expectation(self, mps_tensors, operator, qubits, options=None, normalize=False):
+
+    def contract_expectation(
+        self, mps_tensors, operator, qubits, options=None, normalize=False
+    ):
         """
         Contract the corresponding tensor network to form the state vector representation of the MPS.
 
         Args:
-            mps_tensors: A list of rank-3 ndarray-like tensor objects. 
-                The indices of the ith tensor are expected to be bonding index to the i-1 tensor, 
+            mps_tensors: A list of rank-3 ndarray-like tensor objects.
+                The indices of the ith tensor are expected to be bonding index to the i-1 tensor,
                 the physical mode, and then the bonding index to the i+1th tensor.
-            operator: A ndarray-like tensor object. 
-                The modes of the operator are expected to be output qubits followed by input qubits, e.g, 
-                ``A, B, a, b`` where `a, b` denotes the inputs and `A, B'` denotes the outputs. 
-            qubits: A sequence of integers specifying the qubits that the operator is acting on. 
-            options: Specify the contract and decompose options. 
+            operator: A ndarray-like tensor object.
+                The modes of the operator are expected to be output qubits followed by input qubits, e.g,
+                ``A, B, a, b`` where `a, b` denotes the inputs and `A, B'` denotes the outputs.
+            qubits: A sequence of integers specifying the qubits that the operator is acting on.
+            options: Specify the contract and decompose options.
             normalize: Whether to scale the expectation value by the normalization factor.
 
         Returns:
             An ndarray-like object as the state vector.
         """
-        
+
         interleaved_inputs = []
         extra_mode = 3 * self.num_qubits + 2
         operator_modes = [None] * len(qubits) + [self.bra_modes[q][1] for q in qubits]
@@ -105,19 +112,18 @@ def contract_expectation(self, mps_tensors, operator, qubits, options=None, norm
             if i in qubits:
                 k_modes = (k_modes[0], extra_mode, k_modes[2])
                 q = qubits.index(i)
-                operator_modes[q] = extra_mode # output modes
+                operator_modes[q] = extra_mode  # output modes
                 extra_mode += 1
             interleaved_inputs.extend([o.conj(), k_modes])
         interleaved_inputs.extend([operator, tuple(operator_modes)])
-        interleaved_inputs.append([]) # output
+        interleaved_inputs.append([])  # output
         if normalize:
             norm = self.contract_norm(mps_tensors, options=options)
         else:
             norm = 1
         return self._contract(interleaved_inputs, options=options) / norm
-    
-    def _contract(self, interleaved_inputs, options=None):
 
+    def _contract(self, interleaved_inputs, options=None):
         path = contract_path(*interleaved_inputs, options=options)[0]
-        
-        return contract(*interleaved_inputs, options=options, optimize={'path':path})
+
+        return contract(*interleaved_inputs, options=options, optimize={"path": path})
diff --git a/src/qibotn/MPSUtils.py b/src/qibotn/MPSUtils.py
index 0ad53e16..fd1b4c73 100644
--- a/src/qibotn/MPSUtils.py
+++ b/src/qibotn/MPSUtils.py
@@ -2,73 +2,80 @@
 from cuquantum.cutensornet.experimental import contract_decompose
 from cuquantum import contract
 
+
 def initial(num_qubits, dtype):
     """
-    Generate the MPS with an initial state of |00...00> 
+    Generate the MPS with an initial state of |00...00>
     """
-    state_tensor = cp.asarray([1, 0], dtype=dtype).reshape(1,2,1)
+    state_tensor = cp.asarray([1, 0], dtype=dtype).reshape(1, 2, 1)
     mps_tensors = [state_tensor] * num_qubits
     return mps_tensors
 
-def mps_site_right_swap(
-    mps_tensors, 
-    i, 
-    **kwargs
-):
+
+def mps_site_right_swap(mps_tensors, i, **kwargs):
     """
     Perform the swap operation between the ith and i+1th MPS tensors.
     """
     # contraction followed by QR decomposition
-    a, _, b = contract_decompose('ipj,jqk->iqj,jpk', *mps_tensors[i:i+2], algorithm=kwargs.get('algorithm',None), options=kwargs.get('options',None))
-    mps_tensors[i:i+2] = (a, b)
+    a, _, b = contract_decompose(
+        "ipj,jqk->iqj,jpk",
+        *mps_tensors[i : i + 2],
+        algorithm=kwargs.get("algorithm", None),
+        options=kwargs.get("options", None)
+    )
+    mps_tensors[i : i + 2] = (a, b)
     return mps_tensors
 
-def apply_gate(
-    mps_tensors, 
-    gate, 
-    qubits, 
-    **kwargs
-):
+
+def apply_gate(mps_tensors, gate, qubits, **kwargs):
     """
     Apply the gate operand to the MPS tensors in-place.
-    
+
     Args:
-        mps_tensors: A list of rank-3 ndarray-like tensor objects. 
-            The indices of the ith tensor are expected to be the bonding index to the i-1 tensor, 
+        mps_tensors: A list of rank-3 ndarray-like tensor objects.
+            The indices of the ith tensor are expected to be the bonding index to the i-1 tensor,
             the physical mode, and then the bonding index to the i+1th tensor.
-        gate: A ndarray-like tensor object representing the gate operand. 
-            The modes of the gate is expected to be output qubits followed by input qubits, e.g, 
-            ``A, B, a, b`` where ``a, b`` denotes the inputs and ``A, B`` denotes the outputs. 
+        gate: A ndarray-like tensor object representing the gate operand.
+            The modes of the gate is expected to be output qubits followed by input qubits, e.g,
+            ``A, B, a, b`` where ``a, b`` denotes the inputs and ``A, B`` denotes the outputs.
         qubits: A sequence of integers denoting the qubits that the gate is applied onto.
-        algorithm: The contract and decompose algorithm to use for gate application. 
+        algorithm: The contract and decompose algorithm to use for gate application.
             Can be either a `dict` or a `ContractDecomposeAlgorithm`.
-        options: Specify the contract and decompose options. 
-    
+        options: Specify the contract and decompose options.
+
     Returns:
         The updated MPS tensors.
     """
-    
+
     n_qubits = len(qubits)
     if n_qubits == 1:
         # single-qubit gate
         i = qubits[0]
-        mps_tensors[i] = contract('ipj,qp->iqj', mps_tensors[i], gate, options=kwargs.get('options',None)) # in-place update
+        mps_tensors[i] = contract(
+            "ipj,qp->iqj", mps_tensors[i], gate, options=kwargs.get("options", None)
+        )  # in-place update
     elif n_qubits == 2:
         # two-qubit gate
         i, j = qubits
         if i > j:
             # swap qubits order
-            return apply_gate(mps_tensors, gate.transpose(1,0,3,2), (j, i), **kwargs)
-        elif i+1 == j:
+            return apply_gate(mps_tensors, gate.transpose(1, 0, 3, 2), (j, i), **kwargs)
+        elif i + 1 == j:
             # two adjacent qubits
-            a, _, b = contract_decompose('ipj,jqk,rspq->irj,jsk', *mps_tensors[i:i+2], gate, algorithm=kwargs.get('algorithm',None), options=kwargs.get('options',None))
-            mps_tensors[i:i+2] = (a, b) # in-place update
+            a, _, b = contract_decompose(
+                "ipj,jqk,rspq->irj,jsk",
+                *mps_tensors[i : i + 2],
+                gate,
+                algorithm=kwargs.get("algorithm", None),
+                options=kwargs.get("options", None)
+            )
+            mps_tensors[i : i + 2] = (a, b)  # in-place update
         else:
             # non-adjacent two-qubit gate
             # step 1: swap i with i+1
             mps_site_right_swap(mps_tensors, i, **kwargs)
             # step 2: apply gate to (i+1, j) pair. This amounts to a recursive swap until the two qubits are adjacent
-            apply_gate(mps_tensors, gate, (i+1, j), **kwargs) 
+            apply_gate(mps_tensors, gate, (i + 1, j), **kwargs)
             # step 3: swap back i and i+1
             mps_site_right_swap(mps_tensors, i, **kwargs)
     else:
diff --git a/src/qibotn/QiboCircuitConvertor.py b/src/qibotn/QiboCircuitConvertor.py
index c30cfb64..1a4d688d 100644
--- a/src/qibotn/QiboCircuitConvertor.py
+++ b/src/qibotn/QiboCircuitConvertor.py
@@ -44,8 +44,7 @@ def state_vector_operands(self):
         for key in qubits_frontier:
             out_list.append(qubits_frontier[key])
 
-        operand_exp_interleave = [x for y in zip(
-            operands, mode_labels) for x in y]
+        operand_exp_interleave = [x for y in zip(operands, mode_labels) for x in y]
         operand_exp_interleave.append(out_list)
         return operand_exp_interleave
 

From 78246cf9a0e8cd43ae76e227b5a456eb51126856 Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Thu, 17 Aug 2023 13:48:40 +0800
Subject: [PATCH 24/31] Change filename to that of naming convention

---
 src/qibotn/{MPSContractionHelper.py => mps_contraction_helper.py} | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename src/qibotn/{MPSContractionHelper.py => mps_contraction_helper.py} (100%)

diff --git a/src/qibotn/MPSContractionHelper.py b/src/qibotn/mps_contraction_helper.py
similarity index 100%
rename from src/qibotn/MPSContractionHelper.py
rename to src/qibotn/mps_contraction_helper.py

From 1dba4c3e94daa226ea09a2d56a0b68e0dd33fddf Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Thu, 17 Aug 2023 14:00:11 +0800
Subject: [PATCH 25/31] Update lib name

---
 src/qibotn/cutn.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/qibotn/cutn.py b/src/qibotn/cutn.py
index ab0e9c12..ae36497e 100644
--- a/src/qibotn/cutn.py
+++ b/src/qibotn/cutn.py
@@ -2,7 +2,7 @@
 from cuquantum import contract
 
 from qibotn.QiboCircuitToMPS import QiboCircuitToMPS
-from qibotn.MPSContractionHelper import MPSContractionHelper
+from qibotn.mps_contraction_helper import MPSContractionHelper
 
 
 def eval(qibo_circ, datatype):

From d81cce5e00fb2d565631c5ea26b2f0882f5f33b9 Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Thu, 17 Aug 2023 15:40:35 +0800
Subject: [PATCH 26/31] Change list conversion  using constructor

---
 src/qibotn/mps_contraction_helper.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/qibotn/mps_contraction_helper.py b/src/qibotn/mps_contraction_helper.py
index 9c71200a..ee8e4e42 100644
--- a/src/qibotn/mps_contraction_helper.py
+++ b/src/qibotn/mps_contraction_helper.py
@@ -105,7 +105,7 @@ def contract_expectation(
         interleaved_inputs = []
         extra_mode = 3 * self.num_qubits + 2
         operator_modes = [None] * len(qubits) + [self.bra_modes[q][1] for q in qubits]
-        qubits = [q for q in qubits]
+        qubits = list(qubits)
         for i, o in enumerate(mps_tensors):
             interleaved_inputs.extend([o, self.bra_modes[i]])
             k_modes = self.ket_modes[i]

From d558609bf6477a8adb270d342882a1c8874e5a52 Mon Sep 17 00:00:00 2001
From: Liwei Yang <yang0345@e.ntu.edu.sg>
Date: Tue, 29 Aug 2023 17:29:50 +0800
Subject: [PATCH 27/31] Minor fix of typo and format.

---
 src/qibotn/QiboCircuitConvertor.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/qibotn/QiboCircuitConvertor.py b/src/qibotn/QiboCircuitConvertor.py
index 1a4d688d..2753b9df 100644
--- a/src/qibotn/QiboCircuitConvertor.py
+++ b/src/qibotn/QiboCircuitConvertor.py
@@ -6,7 +6,7 @@ class QiboCircuitToEinsum:
     """Convert a circuit to a Tensor Network (TN) representation.
     The circuit is first processed to an intermediate form by grouping each gate
     matrix with its corresponding qubit it is acting on to a list. It is then
-    converted it to an equivalent TN expression through the class function
+    converted to an equivalent TN expression through the class function
     state_vector_operands() following the Einstein summation convention in the
     interleave format.
 
@@ -44,7 +44,8 @@ def state_vector_operands(self):
         for key in qubits_frontier:
             out_list.append(qubits_frontier[key])
 
-        operand_exp_interleave = [x for y in zip(operands, mode_labels) for x in y]
+        operand_exp_interleave = [x for y in zip(
+            operands, mode_labels) for x in y]
         operand_exp_interleave.append(out_list)
         return operand_exp_interleave
 

From ba442c98307ce0e124d8f2d8f7b5c50f0e498dc4 Mon Sep 17 00:00:00 2001
From: tankya2 <tankya2@ihpc.a-star.edu.sg>
Date: Wed, 30 Aug 2023 10:46:05 +0800
Subject: [PATCH 28/31] Format with black

---
 src/qibotn/QiboCircuitConvertor.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/qibotn/QiboCircuitConvertor.py b/src/qibotn/QiboCircuitConvertor.py
index 2753b9df..9212f70a 100644
--- a/src/qibotn/QiboCircuitConvertor.py
+++ b/src/qibotn/QiboCircuitConvertor.py
@@ -44,8 +44,7 @@ def state_vector_operands(self):
         for key in qubits_frontier:
             out_list.append(qubits_frontier[key])
 
-        operand_exp_interleave = [x for y in zip(
-            operands, mode_labels) for x in y]
+        operand_exp_interleave = [x for y in zip(operands, mode_labels) for x in y]
         operand_exp_interleave.append(out_list)
         return operand_exp_interleave
 

From 49021959025698e4b8fa7a94cda5487454a8ae2b Mon Sep 17 00:00:00 2001
From: Liwei Yang <yang0345@e.ntu.edu.sg>
Date: Wed, 30 Aug 2023 15:58:32 +0800
Subject: [PATCH 29/31] Move cuquantum-python-cu11 package to extras to avoid
 the issue of missing CUDA

---
 setup.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/setup.py b/setup.py
index 43761fc5..38c647c7 100644
--- a/setup.py
+++ b/setup.py
@@ -43,7 +43,6 @@ def version():
         "qibojit>=0.0.7",
         "quimb[tensor]>=1.4.0",
         "cupy>=11.6.0",
-        "cuquantum-python-cu11>=23.3.0",
     ],
     extras_require={
         "docs": [],
@@ -55,6 +54,9 @@ def version():
         "analysis": [
             "pylint>=2.16.0",
         ],
+        "cuquantum": [
+            "cuquantum-python-cu11>=23.3.0",
+        ],
     },
     python_requires=">=3.8.0",
     long_description=(HERE / "README.md").read_text(encoding="utf-8"),

From b84ccadb208143d558d0cfd2483eedb763f3e581 Mon Sep 17 00:00:00 2001
From: Liwei Yang <yang0345@e.ntu.edu.sg>
Date: Wed, 30 Aug 2023 16:41:17 +0800
Subject: [PATCH 30/31] Move cupy package to extras to avoid the issue of
 missing CUDA

---
 setup.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/setup.py b/setup.py
index 38c647c7..e7819456 100644
--- a/setup.py
+++ b/setup.py
@@ -42,7 +42,6 @@ def version():
         "qibo>=0.1.10",
         "qibojit>=0.0.7",
         "quimb[tensor]>=1.4.0",
-        "cupy>=11.6.0",
     ],
     extras_require={
         "docs": [],
@@ -54,7 +53,8 @@ def version():
         "analysis": [
             "pylint>=2.16.0",
         ],
-        "cuquantum": [
+        "cuda": [
+            "cupy>=11.6.0",
             "cuquantum-python-cu11>=23.3.0",
         ],
     },

From e91045f5a3bc978a64dbaa7215e7f52dd2843fc2 Mon Sep 17 00:00:00 2001
From: Liwei Yang <yang0345@e.ntu.edu.sg>
Date: Tue, 12 Sep 2023 12:13:23 +0800
Subject: [PATCH 31/31] Check CUDA_PATH to determine if CUDA is installed and
 if the build job needs to run

---
 .github/workflows/rules.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/rules.yml b/.github/workflows/rules.yml
index c0f29179..d1c06323 100644
--- a/.github/workflows/rules.yml
+++ b/.github/workflows/rules.yml
@@ -8,7 +8,7 @@ on:
 
 jobs:
   build:
-    if: contains(github.event.pull_request.labels.*.name, 'run-workflow') || github.event_name == 'push'
+    if: contains(github.event.pull_request.labels.*.name, 'run-workflow') || github.event_name == 'push' && {{ $CUDA_PATH != '' }}
     strategy:
       matrix:
         os: [ubuntu-latest]