outer and not tensor product

4psireal2 · May 6, 2024 · d0f6f71 · d0f6f71
1 parent 5e91ea1
commit d0f6f71
Show file tree

Hide file tree

Showing 6 changed files with 114 additions and 99 deletions.
diff --git a/src/hpc/contact_process_stat.py b/src/hpc/contact_process_stat.py
@@ -6,7 +6,7 @@
 import scikit_tt.tensor_train as tt
 from scikit_tt.solvers.evp import als
 
-from src.models.contact_process_model import (construct_lindblad, construct_num_op)
+from src.models.contact_process_model import (construct_lindblad)
 from src.utilities.utils import (canonicalize_mps, compute_correlation, compute_dens_dens_corr, compute_purity,
                                  compute_site_expVal, compute_expVal, compute_entanglement_spectrum,
                                  construct_basis_mps, compute_overlap)
@@ -39,6 +39,12 @@
 correlations = np.zeros((len(OMEGAS), L - 1))
 dens_dens_corr = np.zeros((len(OMEGAS), L - 1))
 
+### observable operator
+number_op = np.array([[0, 0], [0, 1]])
+number_op.reshape((1, 2, 2, 1))
+number_mpo = [None]
+number_mpo = tt.TT(number_op)
+
 for i, OMEGA in enumerate(OMEGAS):
     for j, bond_dim in enumerate(bond_dims):
         print(f"Run ALS for {L=}, {OMEGA=} and {bond_dim=}")
@@ -75,7 +81,7 @@
 
         # compute observables
         print("Compute particle numbers")
-        particle_nums = compute_site_expVal(hermit_mps, construct_num_op(L))
+        particle_nums = compute_site_expVal(hermit_mps, number_mpo)
         print(f"Particle number/site: {particle_nums}")
         n_s[j, i] = np.mean(particle_nums)
         print(f"Mean Particle number: {n_s[j, i]}")
@@ -89,20 +95,19 @@
             purities[i] = compute_purity(gs_mps)
             print(f"Purity: {purities[-1]}")
 
-            print("Compute two-point correlation for largest bond dimension")
-            an_op = construct_num_op(1)
+            logger.info("Compute two-point correlation for largest bond dimension")
             for k in range(L - 1):
-                correlations[i, k] = compute_correlation(gs_mps, an_op, r0=0, r1=k + 1)
+                correlations[i, k] = compute_correlation(gs_mps, number_mpo, r0=0, r1=k + 1)
 
             print("Compute density-density correlation for largest bond dimension")
             for k in range(L - 1):
-                dens_dens_corr[i, k] = compute_dens_dens_corr(gs_mps, an_op, r=k + 1)
+                dens_dens_corr[i, k] = compute_dens_dens_corr(gs_mps, number_mpo, r=k + 1)
 
             print("Compute half-chain entanglement spectrum for largest bond dimension")
             entanglement_spectrum[i, :] = compute_entanglement_spectrum(gs_mps)
 
             basis_0 = np.array([1, 0])
-            dark_state = construct_basis_mps(L, basis=[np.kron(basis_0, basis_0)] * L)
+            dark_state = construct_basis_mps(L, basis=[np.outer(basis_0, basis_0)] * L)
             print(f"Overlap with dark state: {compute_overlap(dark_state, gs_mps)}")  #NOTE: negative?
 
 time3 = "{:%Y_%m_%d_%H_%M_%S}".format(datetime.now())

diff --git a/src/hpc/contact_process_stat_exc.py b/src/hpc/contact_process_stat_exc.py
@@ -10,7 +10,7 @@
 import scikit_tt.tensor_train as tt
 from scikit_tt.solvers.evp import als
 
-from src.models.contact_process_model import (construct_lindblad, construct_num_op)
+from src.models.contact_process_model import (construct_lindblad)
 from src.utilities.utils import (canonicalize_mps, compute_correlation, compute_dens_dens_corr, compute_purity,
                                  compute_site_expVal, compute_expVal, compute_entanglement_spectrum,
                                  construct_basis_mps, compute_overlap)
@@ -32,7 +32,7 @@
 conv_eps = 1e-6
 
 ### Stationary simulation
-print("Stationary simulation")
+logger.info("Stationary simulation")
 eval_0 = np.zeros((bond_dims.shape[0], OMEGAS.shape[0]))
 eval_1 = np.zeros((bond_dims.shape[0], OMEGAS.shape[0]))
 evp_residual = np.zeros((bond_dims.shape[0], OMEGAS.shape[0]))
@@ -43,9 +43,15 @@
 correlations = np.zeros((len(OMEGAS), L - 1))
 dens_dens_corr = np.zeros((len(OMEGAS), L - 1))
 
+### observable operator
+number_op = np.array([[0, 0], [0, 1]])
+number_op.reshape((1, 2, 2, 1))
+number_mpo = [None]
+number_mpo = tt.TT(number_op)
+
 for i, OMEGA in enumerate(OMEGAS):
     for j, bond_dim in enumerate(bond_dims):
-        print(f"Run ALS for {L=}, {OMEGA=} and {bond_dim=}")
+        logger.info(f"Run ALS for {L=}, {OMEGA=} and {bond_dim=}")
         lindblad = construct_lindblad(gamma=1.0, omega=OMEGA, L=L)
         lindblad_hermitian = lindblad.transpose(conjugate=True) @ lindblad
 
@@ -55,59 +61,58 @@
         time1 = time.time()
         eigenvalues, eigentensors, _ = als(lindblad_hermitian, mps, number_ev=2, repeats=10, conv_eps=conv_eps, sigma=0)
         time2 = time.time()
-        print(f"Elapsed time: {time2 - time1} seconds")
+        logger.info(f"Elapsed time: {time2 - time1} seconds")
 
         evp_residual[j, i] = (lindblad_hermitian @ eigentensors[1] - eigenvalues[1] * eigentensors[1]).norm()**2
         eval_0[j, i] = eigenvalues[0]
         eval_1[j, i] = eigenvalues[1]
-        print(f"Eigensolver error for first excited state: {evp_residual[j, i]}")
+        logger.info(f"Eigensolver error for first excited state: {evp_residual[j, i]}")
 
-        print(f"Ground state energy per site E: {eigenvalues[1]/L}")
-        print(f"Norm of first excited state: {eigentensors[1].norm()**2}")
+        logger.info(f"First excited state energy per site: {eigenvalues[1]/L}")
+        logger.info(f"Norm of first excited state: {eigentensors[1].norm()**2}")
 
         mps = canonicalize_mps(eigentensors[1])
         mps_dag = mps.transpose(conjugate=True)
 
         # compute non-Hermitian part of mps
         non_hermit_mps = (1 / 2) * (mps - mps_dag)
-        print(f"The norm of the non-Hermitian part: {non_hermit_mps.norm()**2}")
+        logger.info(f"The norm of the non-Hermitian part: {non_hermit_mps.norm()**2}")
 
-        print(f"expVal_1, eigenvalue_1: {compute_expVal(mps, lindblad_hermitian)}, {eigenvalues[1]}")
+        logger.info(f"expVal_0, eigenvalue_0: {compute_expVal(mps, lindblad_hermitian)}, {eigenvalues[1]}")
 
         # compute Hermitian part of mps
         hermit_mps = (1 / 2) * (mps + mps_dag)
 
         # compute observables
-        print("Compute particle numbers")
-        particle_nums = compute_site_expVal(hermit_mps, construct_num_op(L))
-        print(f"Particle number/site: {particle_nums}")
+        logger.info("Compute particle numbers")
+        particle_nums = compute_site_expVal(hermit_mps, number_mpo)
+        logger.info(f"Particle number/site: {particle_nums}")
         n_s[j, i] = np.mean(particle_nums)
-        print(f"Mean Particle number: {n_s[j, i]}")
+        logger.info(f"Mean Particle number: {n_s[j, i]}")
 
         if bond_dim == bond_dims[-1]:
-            print(f"{bond_dim=}")
-            print("Compute spectral gap of L†L for largest bond dimension")
+            logger.info(f"{bond_dim=}")
+            logger.info("Compute spectral gap of L†L for largest bond dimension")
             spectral_gaps[i] = abs(eigenvalues[1] - eigenvalues[0])
 
-            print("Compute purity of state for largest bond dimension")
+            logger.info("Compute purity of state for largest bond dimension")
             purities[i] = compute_purity(mps)
-            print(f"Purity: {purities[-1]}")
+            logger.info(f"Purity: {purities[-1]}")
 
-            print("Compute two-point correlation for largest bond dimension")
-            an_op = construct_num_op(1)
+            logger.info("Compute two-point correlation for largest bond dimension")
             for k in range(L - 1):
-                correlations[i, k] = compute_correlation(mps, an_op, r0=0, r1=k + 1)
+                correlations[i, k] = compute_correlation(mps, number_mpo, r0=0, r1=k + 1)
 
-            print("Compute density-density correlation for largest bond dimension")
+            logger.info("Compute density-density correlation for largest bond dimension")
             for k in range(L - 1):
-                dens_dens_corr[i, k] = compute_dens_dens_corr(mps, an_op, r=k + 1)
+                dens_dens_corr[i, k] = compute_dens_dens_corr(mps, number_mpo, r=k + 1)
 
-            print("Compute half-chain entanglement spectrum for largest bond dimension")
+            logger.info("Compute half-chain entanglement spectrum for largest bond dimension")
             entanglement_spectrum[i, :] = compute_entanglement_spectrum(mps)
 
             basis_0 = np.array([1, 0])
-            dark_state = construct_basis_mps(L, basis=[np.kron(basis_0, basis_0)] * L)
-            print(f"Overlap with dark state: {compute_overlap(dark_state, mps)}")  #NOTE: negative?
+            dark_state = construct_basis_mps(L, basis=[np.outer(basis_0, basis_0)] * L)
+            logger.info(f"Overlap with dark state: {compute_overlap(dark_state, mps)}")  #NOTE: negative?
 
 time3 = "{:%Y_%m_%d_%H_%M_%S}".format(datetime.now())
 

diff --git a/src/hpc/jobs_script.sh b/src/hpc/jobs_script.sh
@@ -75,11 +75,11 @@ log_time() {
 
 # activate virtualenv
 source /scratch/nguyed99/tensor/bin/activate
- 
-# launch Python script
+
+ # launch Python script
 log_memory_cpu_usage & 
 LOG_PID=$!
 
-log_time python3 contact_process_stat.py > "$LOG_PATH/contact_process_stat_L_${L}_${timestamp}.out"
+log_time python3 contact_process_stat_exc.py > "$LOG_PATH/contact_process_stat_exc_L_${L}_${timestamp}.out"
 
 kill $LOG_PID
diff --git a/src/simulations/contact_process_stat.py b/src/simulations/contact_process_stat.py
@@ -5,7 +5,7 @@
 import scikit_tt.tensor_train as tt
 from scikit_tt.solvers.evp import als
 
-from src.models.contact_process_model import (construct_lindblad, construct_num_op)
+from src.models.contact_process_model import (construct_lindblad)
 from src.utilities.utils import (canonicalize_mps, compute_correlation, compute_dens_dens_corr, compute_purity,
                                  compute_site_expVal, compute_expVal, compute_entanglement_spectrum,
                                  construct_basis_mps, compute_overlap)
@@ -36,6 +36,12 @@
 correlations = np.zeros((len(OMEGAS), L - 1))
 dens_dens_corr = np.zeros((len(OMEGAS), L - 1))
 
+### observable operator
+number_op = np.array([[0, 0], [0, 1]])
+number_op.reshape((1, 2, 2, 1))
+number_mpo = [None]
+number_mpo = tt.TT(number_op)
+
 for i, OMEGA in enumerate(OMEGAS):
     for j, bond_dim in enumerate(bond_dims):
         print(f"Run ALS for {L=}, {OMEGA=} and {bond_dim=}")
@@ -72,7 +78,7 @@
 
         # compute observables
         print("Compute particle numbers")
-        particle_nums = compute_site_expVal(hermit_mps, construct_num_op(L))
+        particle_nums = compute_site_expVal(hermit_mps, number_mpo)
         print(f"Particle number/site: {particle_nums}")
         n_s[j, i] = np.mean(particle_nums)
         print(f"Mean Particle number: {n_s[j, i]}")
@@ -87,19 +93,18 @@
             print(f"Purity: {purities[-1]}")
 
             print("Compute two-point correlation for largest bond dimension")
-            an_op = construct_num_op(1)
             for k in range(L - 1):
-                correlations[i, k] = compute_correlation(gs_mps, an_op, r0=0, r1=k + 1)
+                correlations[i, k] = compute_correlation(gs_mps, number_mpo, r0=0, r1=k + 1)
 
             print("Compute density-density correlation for largest bond dimension")
             for k in range(L - 1):
-                dens_dens_corr[i, k] = compute_dens_dens_corr(gs_mps, an_op, r=k + 1)
+                dens_dens_corr[i, k] = compute_dens_dens_corr(gs_mps, number_mpo, r=k + 1)
 
             print("Compute half-chain entanglement spectrum for largest bond dimension")
             entanglement_spectrum[i, :] = compute_entanglement_spectrum(gs_mps)
 
             basis_0 = np.array([1, 0])
-            dark_state = construct_basis_mps(L, basis=[np.kron(basis_0, basis_0)] * L)
+            dark_state = construct_basis_mps(L, basis=[np.outer(basis_0, basis_0)] * L)
             print(f"Overlap with dark state: {compute_overlap(dark_state, gs_mps)}")  #NOTE: negative?
 
 time3 = "{:%Y_%m_%d_%H_%M_%S}".format(datetime.now())

diff --git a/src/utilities/utils.py b/src/utilities/utils.py
@@ -77,20 +77,20 @@ def compute_site_expVal(mps: tt.TT, mpo: tt.TT) -> np.ndarray:
     site_vals = np.zeros(mps.order, dtype=float)
 
     for i in range(mps.order):
-        left_boundary = np.ones(1)
-        right_boundary = np.ones(1)
+        left_boundary = np.ones((1, 1))
+        right_boundary = np.ones((1, 1))
 
         for j in range(mps.order):
+
             if j == i:
-                contraction = np.tensordot(mps.cores[j], mpo.cores[0], axes=([1, 2], [0, 1]))
-                contraction = np.einsum('ijlk->ij', contraction)  # tracing over the physical indices
+                contraction = np.tensordot(mps.cores[j], mpo.cores[0], axes=([1, 2], [2, 1]))
             else:
-                contraction = np.einsum('ikjl->il', mps.cores[j])  # tracing over the physical indices
+                contraction = np.tensordot(mps.cores[j], np.eye(2).reshape(1, 2, 2, 1), axes=([1, 2], [2, 1]))
 
-            left_boundary = np.tensordot(left_boundary, contraction, axes=([0], [0]))
+            left_boundary = np.tensordot(left_boundary, contraction, axes=([0, 1], [0, 2]))
 
-        if (left_boundary @ right_boundary).imag < 1e-12:
-            site_vals[i] = (left_boundary @ right_boundary).real
+        if (left_boundary @ right_boundary).item().imag < 1e-12:
+            site_vals[i] = (left_boundary @ right_boundary).item().real
         else:
             raise ValueError("Complex expectation value is found.")
 
@@ -129,50 +129,49 @@ def compute_correlation(mps: tt.TT, mpo: tt.TT, r0: int, r1: int) -> float:
     - r0, r1: first, second index of sites on the TT 
     """
 
-    left_boundary = np.ones(1)
-    right_boundary = np.ones(1)
+    left_boundary = np.ones((1, 1))
+    right_boundary = np.ones((1, 1))
 
     # compute <O_{r0} . O_{r1}>
-    for i in range(mps.order):
-        if i != r0 and i != r1:
-            contraction = np.einsum('ikjl->il', mps.cores[i])
+    for j in range(mps.order):
+        if j != r0 and j != r1:
+            contraction = np.tensordot(mps.cores[j], np.eye(2).reshape(1, 2, 2, 1), axes=([1, 2], [2, 1]))
         else:
-            contraction = np.tensordot(mps.cores[i], mpo.cores[0], axes=([1, 2], [0, 1]))
-            contraction = np.einsum('ijlk->ij', contraction)  # tracing over the physical indices
-        left_boundary = np.tensordot(left_boundary, contraction, axes=([0], [0]))
+            contraction = np.tensordot(mps.cores[j], mpo.cores[0], axes=([1, 2], [2, 1]))
+
+        left_boundary = np.tensordot(left_boundary, contraction, axes=([0, 1], [0, 2]))
 
-    mean_product = left_boundary @ right_boundary
+    mean_product = (left_boundary @ right_boundary).item()
 
     # compute <O_{r0}><O_{r1}>
     # compute <O_{r0}>
-    left_boundary = np.ones(1)
-    right_boundary = np.ones(1)
+
+    left_boundary = np.ones((1, 1))
+    right_boundary = np.ones((1, 1))
 
     for j in range(mps.order):
         if j == r0:
-            contraction = np.tensordot(mps.cores[j], mpo.cores[0], axes=([1, 2], [0, 1]))
-            contraction = np.einsum('ijlk->ij', contraction)  # tracing over the physical indices
+            contraction = np.tensordot(mps.cores[j], mpo.cores[0], axes=([1, 2], [2, 1]))
         else:
-            contraction = np.einsum('ikjl->il', mps.cores[j])  # tracing over the physical indices
+            contraction = np.tensordot(mps.cores[j], np.eye(2).reshape(1, 2, 2, 1), axes=([1, 2], [2, 1]))
 
-        left_boundary = np.tensordot(left_boundary, contraction, axes=([0], [0]))
+        left_boundary = np.tensordot(left_boundary, contraction, axes=([0, 1], [0, 2]))
 
-    product_mean = left_boundary @ right_boundary
+    product_mean = (left_boundary @ right_boundary).item()
 
     # compute <O_{r1}>
-    left_boundary = np.ones(1)
-    right_boundary = np.ones(1)
+    left_boundary = np.ones((1, 1))
+    right_boundary = np.ones((1, 1))
 
     for j in range(mps.order):
         if j == r1:
-            contraction = np.tensordot(mps.cores[j], mpo.cores[0], axes=([1, 2], [0, 1]))
-            contraction = np.einsum('ijlk->ij', contraction)  # tracing over the physical indices
+            contraction = np.tensordot(mps.cores[j], mpo.cores[0], axes=([1, 2], [2, 1]))
         else:
-            contraction = np.einsum('ikjl->il', mps.cores[j])  # tracing over the physical indices
+            contraction = np.tensordot(mps.cores[j], np.eye(2).reshape(1, 2, 2, 1), axes=([1, 2], [2, 1]))
 
-        left_boundary = np.tensordot(left_boundary, contraction, axes=([0], [0]))
+        left_boundary = np.tensordot(left_boundary, contraction, axes=([0, 1], [0, 2]))
 
-    product_mean *= left_boundary @ right_boundary
+    product_mean *= (left_boundary @ right_boundary).item()
 
     return mean_product - product_mean
 
@@ -186,34 +185,33 @@ def compute_dens_dens_corr(mps: tt.TT, mpo: tt.TT, r: int) -> float:
     - mpo: 1 core with shape (2,2,2,2)
     """
 
-    left_boundary = np.ones(1)
-    right_boundary = np.ones(1)
+    left_boundary = np.ones((1, 1))
+    right_boundary = np.ones((1, 1))
 
-    # compute <O_{r0} . O_{r1}>
-    for i in range(mps.order):
-        if i != 0 and i != r:
-            contraction = np.einsum('ikjl->il', mps.cores[i])
+    # compute <O_{0} . O_{r}>
+    for j in range(mps.order):
+        if j != 0 and j != r:
+            contraction = np.tensordot(mps.cores[j], np.eye(2).reshape(1, 2, 2, 1), axes=([1, 2], [2, 1]))
         else:
-            contraction = np.tensordot(mps.cores[i], mpo.cores[0], axes=([1, 2], [0, 1]))
-            contraction = np.einsum('ijlk->ij', contraction)  # tracing over the physical indices
-        left_boundary = np.tensordot(left_boundary, contraction, axes=([0], [0]))
+            contraction = np.tensordot(mps.cores[j], mpo.cores[0], axes=([1, 2], [2, 1]))
+
+        left_boundary = np.tensordot(left_boundary, contraction, axes=([0, 1], [0, 2]))
 
-    mean_product = left_boundary @ right_boundary
+    mean_product = (left_boundary @ right_boundary).item()
 
     # compute <O_{0}>²
-    left_boundary = np.ones(1)
-    right_boundary = np.ones(1)
+    left_boundary = np.ones((1, 1))
+    right_boundary = np.ones((1, 1))
 
     for j in range(mps.order):
         if j == 0:
-            contraction = np.tensordot(mps.cores[j], mpo.cores[0], axes=([1, 2], [0, 1]))
-            contraction = np.einsum('ijlk->ij', contraction)  # tracing over the physical indices
+            contraction = np.tensordot(mps.cores[j], mpo.cores[0], axes=([1, 2], [2, 1]))
         else:
-            contraction = np.einsum('ikjl->il', mps.cores[j])  # tracing over the physical indices
+            contraction = np.tensordot(mps.cores[j], np.eye(2).reshape(1, 2, 2, 1), axes=([1, 2], [2, 1]))
 
-        left_boundary = np.tensordot(left_boundary, contraction, axes=([0], [0]))
+        left_boundary = np.tensordot(left_boundary, contraction, axes=([0, 1], [0, 2]))
 
-    product_mean = (left_boundary @ right_boundary)**2
+    product_mean = ((left_boundary @ right_boundary).item())**2
     return mean_product - product_mean