From d0f6f719b225e6bd8482fee29e981be8743bc7a6 Mon Sep 17 00:00:00 2001 From: nguyed99 Date: Mon, 6 May 2024 13:33:35 +0200 Subject: [PATCH] outer and not tensor product --- src/hpc/contact_process_stat.py | 19 +++-- src/hpc/contact_process_stat_exc.py | 55 ++++++++------- src/hpc/jobs_script.sh | 6 +- src/simulations/contact_process_stat.py | 17 +++-- src/utilities/utils.py | 94 ++++++++++++------------- tests/test_utilities.py | 22 +++--- 6 files changed, 114 insertions(+), 99 deletions(-) diff --git a/src/hpc/contact_process_stat.py b/src/hpc/contact_process_stat.py index 06ac696..ce033fd 100644 --- 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 index 33c6d55..0573c15 100644 --- 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 index 59ae9a0..52e29f8 100644 --- 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 \ No newline at end of file diff --git a/src/simulations/contact_process_stat.py b/src/simulations/contact_process_stat.py index 9b61ed0..40810d1 100644 --- 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 index aa0606d..bb8fb8f 100644 --- 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 - 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 # compute - 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 - 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 - for i in range(mps.order): - if i != 0 and i != r: - contraction = np.einsum('ikjl->il', mps.cores[i]) + # compute + 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 ² - 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 diff --git a/tests/test_utilities.py b/tests/test_utilities.py index a2c5e17..431b4a3 100644 --- a/tests/test_utilities.py +++ b/tests/test_utilities.py @@ -10,7 +10,6 @@ import scikit_tt.tensor_train as tt import src.utilities.utils as utils -import src.models.contact_process_model as model logger = logging.getLogger(__name__) logging.basicConfig(level=logging.INFO) @@ -51,18 +50,21 @@ def test_mixed_canonical_form(self): def test_num_op(self): mps_tests = [ - utils.construct_basis_mps(L, basis=[np.kron(basis_1, basis_1)] * L), - utils.construct_basis_mps(L, basis=[np.kron(basis_1, basis_0)] * L), - utils.construct_basis_mps(L, basis=[np.kron(basis_0, basis_1)] * L), - utils.construct_basis_mps(L, basis=[np.kron(basis_0, basis_0)] * L) + utils.construct_basis_mps(L, basis=[np.outer(basis_1, basis_1)] * L), + utils.construct_basis_mps(L, basis=[np.outer(basis_1, basis_0)] * L), + utils.construct_basis_mps(L, basis=[np.outer(basis_0, basis_1)] * L), + utils.construct_basis_mps(L, basis=[np.outer(basis_0, basis_0)] * L) ] - num_op_chain = model.construct_num_op(L) + number_op = np.array([[0, 0], [0, 1]]) + number_op.reshape((1, 2, 2, 1)) + number_mpo = [None] + number_mpo = tt.TT(number_op) # expected arrays of single site expectation values - site_vals = [2 * np.ones(L), np.ones(L), np.ones(L), np.zeros(L)] - for i, mps in enumerate(mps_tests): - site_expVal = utils.compute_site_expVal(mps, num_op_chain) - assert np.array_equal(site_expVal, site_vals[i]) + site_vals = [np.ones(L), np.zeros(L), np.zeros(L), np.zeros(L)] + for k, mps in enumerate(mps_tests): + site_expVal = utils.compute_site_expVal(mps, number_mpo) + assert np.array_equal(site_expVal, site_vals[k]) def test_purity(self): mps_1 = utils.construct_basis_mps(L, basis=[np.kron(basis_1, basis_1)] * L)