Merge branch 'MPS-to-vector-feature' into Noise_Characterization

Author: MaxFroehlich1410

examples/1_Noisy_Ising_Model.py

@@ -11,11 +11,10 @@
 
 # Define the system Hamiltonian
 L = 10
-d = 2
 J = 1
 g = 0.5
 H_0 = MPO()
-H_0.init_Ising(L, d, J, g)
+H_0.init_Ising(L, J, g)
 
 # Define the initial state
 state = MPS(L, state='zeros')

examples/2_Equivalence_Checking.py

@@ -4,7 +4,7 @@
 from qiskit.circuit.library.n_local import TwoLocal
 import qiskit.compiler
 
-from yaqs.circuits.equivalence_checking import equivalence_checker
+from yaqs.circuits import equivalence_checker
 
 # Define the initial circuit
 num_qubits = 5

examples/5_Benchmarker.py

@@ -5,6 +5,6 @@
 num_qubits = 10
 model = {'name': 'Ising', 'L': num_qubits, 'J': 1, 'g': 0.5}
 demo_circuit = create_Ising_circuit(model, dt=0.1, timesteps=10)
-
-# Run the benchmark on the demo circuit.
-benchmarker.run(demo_circuit, style='dots')
+if __name__ == "__main__":
+    # Run the benchmark on the demo circuit.
+    benchmarker.run(demo_circuit, style='dots')

src/yaqs/circuits/CircuitTJM.py

@@ -10,7 +10,7 @@
 from yaqs.core.methods.dissipation import apply_dissipation
 from yaqs.core.methods.stochastic_process import stochastic_process
 from yaqs.core.methods.operations import measure
-from yaqs.circuits.dag.dag_utils import convert_dag_to_tensor_algorithm
+from yaqs.circuits.utils.dag_utils import convert_dag_to_tensor_algorithm
 
 
 from typing import TYPE_CHECKING, Union
@@ -78,7 +78,6 @@ def construct_generator_MPO(gate, length: int) -> MPO | int | int:
             W = np.zeros((1, 1, 2, 2), dtype=complex)
             W[0, 0] = gate.generator[second_gen]
             tensors.append(W)
-            # break
         else:
             W = np.zeros((1, 1, 2, 2), dtype=complex)
             W[0, 0] = np.eye(2)

…ivalence_checking/equivalence_checker.py src/yaqs/circuits/equivalence_checker.pysrc/yaqs/circuits/equivalence_checking/equivalence_checker.py renamed to src/yaqs/circuits/equivalence_checker.py src/yaqs/circuits/equivalence_checking/equivalence_checker.py renamed to src/yaqs/circuits/equivalence_checker.py

@@ -3,7 +3,7 @@
 import time
 
 from yaqs.core.data_structures.networks import MPO
-from yaqs.circuits.equivalence_checking.mpo_utils import iterate
+from yaqs.circuits.utils.mpo_utils import iterate
 
 from typing import TYPE_CHECKING
 if TYPE_CHECKING:

src/yaqs/circuits/dag/dag_utils.py src/yaqs/circuits/utils/dag_utils.pysrc/yaqs/circuits/dag/dag_utils.py renamed to src/yaqs/circuits/utils/dag_utils.py

@@ -3,7 +3,7 @@
 import opt_einsum as oe
 from qiskit.converters import dag_to_circuit
 
-from yaqs.circuits.dag.dag_utils import check_longest_gate, convert_dag_to_tensor_algorithm, get_temporal_zone, select_starting_point
+from yaqs.circuits.utils.dag_utils import check_longest_gate, convert_dag_to_tensor_algorithm, get_temporal_zone, select_starting_point
 
 from typing import TYPE_CHECKING
 if TYPE_CHECKING:

…rcuits/equivalence_checking/mpo_utils.py src/yaqs/circuits/utils/mpo_utils.pysrc/yaqs/circuits/equivalence_checking/mpo_utils.py renamed to src/yaqs/circuits/utils/mpo_utils.py src/yaqs/circuits/equivalence_checking/mpo_utils.py renamed to src/yaqs/circuits/utils/mpo_utils.py

@@ -13,6 +13,7 @@
 # Convention (sigma, chi_l-1, chi_l)
 class MPS:
     def __init__(self, length: int, tensors: list = None, physical_dimensions: list = None, state: str = 'zeros'):
+        self.flipped = False
         if tensors is not None:
             assert len(tensors) == length
             self.tensors = tensors
@@ -65,14 +66,19 @@ def __init__(self, length: int, tensors: list = None, physical_dimensions: list
                         vector[0] = 1
                     else:
                         vector[1] = 1
+                elif state == 'random':
+                    vector[0] = np.random.rand()
+                    vector[1] = 1 - vector[0]
                 else:
                     raise ValueError("Invalid state string")
 
                 tensor = np.expand_dims(vector, axis=(0, 1))
 
                 tensor = np.transpose(tensor, (2, 0, 1))
                 self.tensors.append(tensor)
-        self.flipped = False
+
+            if state == 'random':
+                self.normalize()
 
     def write_max_bond_dim(self) -> int:
         global_max = 0
@@ -279,10 +285,44 @@ def convert_to_vector(self) -> np.ndarray:
         vec = vec.flatten()
         return vec
 
+    def convert_to_vector(self) -> np.ndarray:
+        """
+        Converts the MPS to a full state vector representation.
+
+        Returns:
+                A one-dimensional NumPy array of length \(\prod_{\ell=1}^L d_\ell\) 
+                representing the state vector.
+        """
+        # Start with the first tensor.
+        # Assume each tensor has shape (d, chi_left, chi_right) with chi_left=1 for the first tensor.
+        vec = self.tensors[0]  # shape: (d_1, 1, chi_1)
+
+        # Contract sequentially with the remaining tensors.
+        for i in range(1, self.length):
+            # Contract the last bond of vec with the middle index (left bond) of the next tensor.
+            vec = np.tensordot(vec, self.tensors[i], axes=([-1], [1]))
+            # After tensordot, if vec had shape (..., chi_i) and the new tensor has shape (d_{i+1}, chi_i, chi_{i+1}),
+            # then vec now has shape (..., d_{i+1}, chi_{i+1}).
+            # Reshape to merge all physical indices into one index.
+            new_shape = (-1, vec.shape[-1])
+            vec = np.reshape(vec, new_shape)
+
+        # At the end, the final bond dimension should be 1.
+        vec = np.squeeze(vec, axis=-1)
+        # Flatten the resulting multi-index into a one-dimensional state vector.
+        vec = vec.flatten()
+        return vec
+
+
+
+
+    
+
 
 # Convention (sigma, sigma', chi_l,  chi_l+1)
 class MPO:
-    def init_Ising(self, length: int, physical_dimension: int, J: float, g: float):
+    def init_Ising(self, length: int, J: float, g: float):
+        physical_dimension = 2
         zero = np.zeros((physical_dimension, physical_dimension), dtype=complex)
         identity = np.eye(physical_dimension, dtype=complex)
         X = getattr(GateLibrary, "x")().matrix
@@ -323,7 +363,8 @@ def init_Ising(self, length: int, physical_dimension: int, J: float, g: float):
         self.length = length
         self.physical_dimension = physical_dimension
 
-    def init_Heisenberg(self, length: int, physical_dimension: int, Jx: float, Jy: float, Jz: float, h: float):
+    def init_Heisenberg(self, length: int, Jx: float, Jy: float, Jz: float, h: float):
+        physical_dimension = 2
         zero = np.zeros((physical_dimension, physical_dimension), dtype=complex)
         identity = np.eye(physical_dimension, dtype=complex)
         X = getattr(GateLibrary, "x")().matrix

src/yaqs/core/data_structures/networks.py

@@ -11,6 +11,7 @@ def __init__(self, name: str, site: int):
         self.trajectories = None
 
     def initialize(self, sim_params):
+        self.results = None
         if type(sim_params) == PhysicsSimParams:
             if sim_params.sample_timesteps:
                 self.trajectories = np.empty((sim_params.N, len(sim_params.times)), dtype=float)

src/yaqs/core/data_structures/simulation_parameters.py

@@ -17,19 +17,22 @@ def create_Ising_circuit(model, dt, timesteps):
             circ.rx(theta=alpha, qubit=site)
 
         for site in range(model['L'] // 2):
-            circ.cx(control_qubit=2*site, target_qubit=2*site+1)
-            circ.rz(phi=beta, qubit=2*site+1)
-            circ.cx(control_qubit=2*site, target_qubit=2*site+1)
+            circ.rzz(beta, qubit1=2*site, qubit2=2*site+1)
+            # circ.cx(control_qubit=2*site, target_qubit=2*site+1)
+            # circ.rz(phi=beta, qubit=2*site+1)
+            # circ.cx(control_qubit=2*site, target_qubit=2*site+1)
 
         for site in range(1, model['L'] // 2):
-            circ.cx(control_qubit=2*site-1, target_qubit=2*site)
-            circ.rz(phi=beta, qubit=2*site)
-            circ.cx(control_qubit=2*site-1, target_qubit=2*site)
+            circ.rzz(beta, qubit1=2*site-1, qubit2=2*site)
+            # circ.cx(control_qubit=2*site-1, target_qubit=2*site)
+            # circ.rz(phi=beta, qubit=2*site)
+            # circ.cx(control_qubit=2*site-1, target_qubit=2*site)
 
         if model['L'] % 2 != 0 and model['L'] != 1:
-            circ.cx(control_qubit=model['L']-2, target_qubit=model['L']-1)
-            circ.rz(phi=beta, qubit=model['L']-1)
-            circ.cx(control_qubit=model['L']-2, target_qubit=model['L']-1)
+            circ.rzz(beta, qubit1=model['L']-2, qubit2=model['L']-1)
+            # circ.cx(control_qubit=model['L']-2, target_qubit=model['L']-1)
+            # circ.rz(phi=beta, qubit=model['L']-1)
+            # circ.cx(control_qubit=model['L']-2, target_qubit=model['L']-1)
 
     return circ