diff --git a/examples/tenpy_sz_convention.py b/examples/tenpy_sz_convention.py new file mode 100644 index 00000000..bca7b84a --- /dev/null +++ b/examples/tenpy_sz_convention.py @@ -0,0 +1,163 @@ +""" +Demonstrates the different internal basis conventions between +TeNPy's TFIChain and XXZChain models and showcases a robust +method for handling these inconsistencies when converting to TensorCircuit. + +1. TFIChain: Shows a direct conversion works perfectly. +2. XXZChain: + a. Runs DMRG to obtain a non-trivial ground state. + b. Shows that direct conversion leads to incorrect expectation values for correlation functions. + c. Demonstrates that applying a layer of X-gates in TensorCircuit +3. Tensor Dissection: Provides definitive proof of the differing internal basis conventions between the two models. +""" + +import numpy as np +from tenpy.networks.mps import MPS +from tenpy.models.tf_ising import TFIChain +from tenpy.models.xxz_chain import XXZChain +from tenpy.algorithms import dmrg + +import tensorcircuit as tc + +print("Scenario 1: Antiferromagnetic State (TFIChain)") +L = 10 +afm_model_params = {"L": L, "bc_MPS": "finite", "J": 1.0, "g": 0.0, "conserve": None} +afm_M = TFIChain(afm_model_params) +afm_state = ["up", "down"] * (L // 2) +print(f"Testing with a simple product state: {afm_state}") + +psi_afm = MPS.from_product_state(afm_M.lat.mps_sites(), afm_state, bc=afm_M.lat.bc_MPS) +tc_afm_state = tc.quantum.tenpy2qop(psi_afm) +circuit_afm = tc.MPSCircuit(L, wavefunction=tc_afm_state) + +mag_z_afm_tenpy = psi_afm.expectation_value("Sz") +mag_z_afm_tc = np.array( + [ + tc.backend.numpy(tc.backend.real(circuit_afm.expectation((tc.gates.z(), [i])))) + / 2.0 + for i in range(L) + ] +) +print("\nAntiferromagnetic state site-by-site magnetization comparison:") +print("TeNPy:", np.round(mag_z_afm_tenpy, 8)) +print("TC: ", np.round(mag_z_afm_tc, 8)) +np.testing.assert_allclose(mag_z_afm_tenpy, mag_z_afm_tc, atol=1e-5) +print( + "\n[SUCCESS] TFI-based Antiferromagnetic state matches perfectly with the pure converter." +) + + +print("Scenario 2: XXZChain Model") +xxz_model_params = {"L": L, "bc_MPS": "finite", "Jxx": 1.0, "Jz": 0.5, "hz": 0.1} +xxz_M = XXZChain(xxz_model_params) +example_state = ["up", "down", "up", "up", "down", "down", "up", "down", "down", "up"] +print(f"Testing with a random product state: {example_state}") +psi_rand_xxz = MPS.from_product_state( + xxz_M.lat.mps_sites(), example_state, bc=xxz_M.lat.bc_MPS +) +tc_rand_xxz_state = tc.quantum.tenpy2qop(psi_rand_xxz) +circuit_rand_xxz = tc.MPSCircuit(L, wavefunction=tc_rand_xxz_state) +mag_z_rand_xxz_tenpy = psi_rand_xxz.expectation_value("Sz") +mag_z_rand_xxz_tc = np.array( + [ + tc.backend.numpy( + tc.backend.real(circuit_rand_xxz.expectation((tc.gates.z(), [i]))) + ) + / 2.0 + for i in range(L) + ] +) +print("\nXXZ-based random state site-by-site magnetization comparison:") +print("TeNPy:", np.round(mag_z_rand_xxz_tenpy, 8)) +print("TC: ", np.round(mag_z_rand_xxz_tc, 8)) +try: + np.testing.assert_allclose(mag_z_rand_xxz_tenpy, mag_z_rand_xxz_tc, atol=1e-5) +except AssertionError as e: + print("\n[SUCCESS] As expected, the direct comparison fails for XXZChain.") + print( + "This is because the pure converter does not handle its inverted basis convention." + ) + print("\nVerifying that the values match after correcting the sign:") + np.testing.assert_allclose(mag_z_rand_xxz_tenpy, -mag_z_rand_xxz_tc, atol=1e-5) + print( + "[SUCCESS] Test passes after applying the sign correction for the XXZChain model." + ) + + +print("Scenario 3: Tensor Dissection for Both Models") +simple_L = 2 +simple_labels = ["up", "down"] +print("\nDissecting TFIChain-based Tensors") +sites_tfi = afm_M.lat.mps_sites()[:simple_L] +psi_simple_tfi = MPS.from_product_state(sites_tfi, simple_labels, bc="finite") +for i, label in enumerate(simple_labels): + B_tensor = psi_simple_tfi.get_B(i).to_ndarray() + print( + f"For '{label}', TFIChain internal tensor has non-zero at physical index {np.where(B_tensor[0,:,0] != 0)[0][0]}" + ) +print("\nDissecting XXZChain-based Tensors") +sites_xxz = xxz_M.lat.mps_sites()[:simple_L] +psi_simple_xxz = MPS.from_product_state(sites_xxz, simple_labels, bc="finite") +for i, label in enumerate(simple_labels): + B_tensor = psi_simple_xxz.get_B(i).to_ndarray() + print( + f"For '{label}', XXZChain internal tensor has non-zero at physical index {np.where(B_tensor[0,:,0] != 0)[0][0]}" + ) +print("\n Conclusion") +print("The dissection above shows the root cause of the different behaviors:") +print( + " - TFIChain's 'up' maps to index 0, 'down' to index 1. This matches TC's standard." +) +print( + " - XXZChain's 'up' maps to index 1, 'down' to index 0. This is INVERTED from TC's standard." +) +print("\nTherefore, a single, universal converter is not feasible without context.") +print( + "The robust solution is to use a pure converter and apply corrections on a case-by-case basis," +) +print("or to create model-specific converters.") + + +print("--- Scenario 3: Correcting XXZChain DMRG state with X-gates ---") + +L = 10 +xxz_model_params = {"L": L, "bc_MPS": "finite", "Jxx": 1.0, "Jz": 1.0, "conserve": None} +xxz_M = XXZChain(xxz_model_params) +psi0_xxz = MPS.from_product_state( + xxz_M.lat.mps_sites(), ["up", "down"] * (L // 2), bc=xxz_M.lat.bc_MPS +) +dmrg_params = {"max_sweeps": 10, "trunc_params": {"chi_max": 64}} +eng = dmrg.TwoSiteDMRGEngine(psi0_xxz, xxz_M, dmrg_params) +E, psi_gs_xxz = eng.run() +print(f"XXZ DMRG finished. Ground state energy: {E:.10f}") + +state_raw_quvector = tc.quantum.tenpy2qop(psi_gs_xxz) + +i, j = L // 2 - 1, L // 2 +corr_tenpy = psi_gs_xxz.correlation_function("Sz", "Sz", sites1=[i], sites2=[j])[0, 0] +print("\nApplying X-gate to each qubit to correct the basis convention...") +circuit_to_be_corrected = tc.MPSCircuit(L, wavefunction=state_raw_quvector) + +for k in range(L): + circuit_to_be_corrected.x(k) + +corr_tc_corrected = ( + tc.backend.real( + circuit_to_be_corrected.expectation((tc.gates.z(), [i]), (tc.gates.z(), [j])) + ) + / 4.0 +) + +print(f"\nComparing correlation function for the DMRG ground state:") +print(f"TeNPy (Ground Truth): {corr_tenpy:.8f}") +print(f"TC (after X-gate correction): {corr_tc_corrected:.8f}") +np.testing.assert_allclose(corr_tenpy, corr_tc_corrected, atol=1e-5) +print( + "\n[SUCCESS] The correlation functions match perfectly after applying the X-gate correction." +) +print( + "This demonstrates the recommended physical approach to handle the XXZChain's inverted basis convention." +) + + +print("\n\nWorkflow demonstration and analysis complete!") diff --git a/requirements/requirements-extra.txt b/requirements/requirements-extra.txt index 51370e2d..22416202 100644 --- a/requirements/requirements-extra.txt +++ b/requirements/requirements-extra.txt @@ -5,6 +5,8 @@ qiskit-nature mitiq cirq torch==2.2.2 +physics-tenpy==1.0.6 +quimb==1.11.1 # jupyter mthree==1.1.0 openfermion diff --git a/tensorcircuit/quantum.py b/tensorcircuit/quantum.py index 1266ef48..6ec08aa1 100644 --- a/tensorcircuit/quantum.py +++ b/tensorcircuit/quantum.py @@ -23,6 +23,7 @@ ) import numpy as np +import tensornetwork as tn from tensornetwork.network_components import AbstractNode, CopyNode, Edge, Node, connect from tensornetwork.network_operations import ( copy, @@ -1151,33 +1152,281 @@ def generate_local_hamiltonian( return hop -def tn2qop(tn_mpo: Any) -> QuOperator: +# TODO(@Charlespkuer): Add more conversion functions for other packages +def extract_tensors_from_qop(qop: QuOperator) -> Tuple[List[Node], bool, int]: """ - Convert MPO in TensorNetwork package to QuOperator. + Extract and sort tensors from QuOperator for conversion to other tensor network formats. - :param tn_mpo: MPO in the form of TensorNetwork package - :type tn_mpo: ``tn.matrixproductstates.mpo.*`` - :return: MPO in the form of QuOperator + :param qop: Input QuOperator to extract tensors from + :type qop: QuOperator + :return: Tuple containing (sorted_nodes, is_mps, nwires) where: + - sorted_nodes: List of Node objects sorted in linear chain order + - is_mps: Boolean flag indicating if the structure is MPS (True) or MPO (False) + - nwires: Integer number of physical edges/qubits in the system + :rtype: Tuple[List[Node], bool, int] + """ + is_mps = len(qop.in_edges) == 0 + nwires = len(qop.out_edges) + if not is_mps and len(qop.in_edges) != nwires: + raise ValueError( + "MPO must have the same number of input and output edges. " + f"Got {len(qop.in_edges)} and {nwires}." + ) + + # Collect all nodes from edges + nodes_for_sorting = qop.nodes + if len(nodes_for_sorting) != nwires: + raise ValueError(f"Number of nodes does not match number of wires.") + + # Find endpoint nodes + endpoint_nodes = set() + physical_edges = set(qop.out_edges) if is_mps else set(qop.in_edges + qop.out_edges) + + if is_mps: + rank_2_nodes = {node for node in nodes_for_sorting if len(node.edges) == 2} + if len(rank_2_nodes) == 2: + endpoint_nodes = rank_2_nodes + + if not endpoint_nodes: + endpoint_nodes = {edge.node1 for edge in qop.ignore_edges if edge.node1} + + if not endpoint_nodes and len(nodes_for_sorting) > 1: + virtual_bond_counts = {} + virtual_bond_dim_sums = {} + + for node in nodes_for_sorting: + virtual_bonds = 0 + virtual_dim_sum = 0 + + for edge in node.edges: + if edge not in physical_edges and not edge.is_dangling(): + virtual_bonds += 1 + virtual_dim_sum += edge.dimension + + virtual_bond_counts[node] = virtual_bonds + virtual_bond_dim_sums[node] = virtual_dim_sum + + min_dim_sum = min(virtual_bond_dim_sums.values()) + min_dim_nodes = { + node + for node, dim_sum in virtual_bond_dim_sums.items() + if dim_sum == min_dim_sum + } + + if len(min_dim_nodes) == 2: + endpoint_nodes = min_dim_nodes + + if not endpoint_nodes: + if len(nodes_for_sorting) == 1: + raise ValueError("Cannot determine chain structure: only one node found.") + elif len(nodes_for_sorting) >= 2: + raise ValueError(f"Cannot identify endpoint nodes for your nodes.") + + # Sort nodes along the chain + sorted_nodes: list[Node] = [] + if endpoint_nodes and len(endpoint_nodes) >= 1: + current = next(iter(endpoint_nodes)) + while current and len(sorted_nodes) < nwires: + sorted_nodes.append(current) + current = next( + ( + e.node2 if e.node1 is current else e.node1 + for e in current.edges + if not e.is_dangling() + and e not in physical_edges + and (e.node2 if e.node1 is current else e.node1) not in sorted_nodes + ), + None, + ) + + if not sorted_nodes: + raise ValueError("No valid chain structure found in the QuOperator. ") + if len(sorted_nodes) > 0 and len(qop.ignore_edges) > 0: + if sorted_nodes[0] is not qop.ignore_edges[0].node1: + sorted_nodes = sorted_nodes[::-1] + + return sorted_nodes, is_mps, nwires + + +def tenpy2qop(tenpy_obj: Any) -> QuOperator: + """ + Converts a TeNPy MPO or MPS to a TensorCircuit QuOperator. + This definitive version correctly handles axis ordering and boundary + conditions to be compatible with `eval_matrix`. + + :param tenpy_obj: A MPO or MPS object from the TeNPy package. + :type tenpy_obj: Union[tenpy.networks.mpo.MPO, tenpy.networks.mps.MPS] + :return: The corresponding state or operator as a QuOperator. :rtype: QuOperator """ - tn_mpo = tn_mpo.tensors - nwires = len(tn_mpo) - mpo = [] - for i in range(nwires): - mpo.append(Node(tn_mpo[i], name=f"mpo_{i}")) + # MPO objects have _W attribute containing tensor list (documented in tenpy.networks.mpo.MPO) + # MPS objects have _B attribute containing tensor list (documented in tenpy.networks.mps.MPS) + # These are internal attributes that store the actual tensor data for each site + # Reference: https://tenpy.readthedocs.io/en/latest/reference/tenpy.networks.mpo.html + # https://tenpy.readthedocs.io/en/latest/reference/tenpy.networks.mps.html + is_mpo = hasattr(tenpy_obj, "_W") + tenpy_tensors = tenpy_obj._W if is_mpo else tenpy_obj._B + nwires = len(tenpy_tensors) + if nwires == 0: + return quantum_constructor([], [], []) + + nodes = [] + if is_mpo: + original_tensors_obj = tenpy_tensors + + for i, W_obj in enumerate(original_tensors_obj): + arr = W_obj.to_ndarray() + labels = W_obj.get_leg_labels() + wL_idx = labels.index("wL") + p_idx = labels.index("p") + p_star_idx = labels.index("p*") + wR_idx = labels.index("wR") + + arr_reordered = arr.transpose((wL_idx, p_idx, p_star_idx, wR_idx)) + if nwires == 1: + arr_reordered = arr_reordered[[0], :, :, :] + arr_reordered = arr_reordered[:, :, :, [-1]] + else: + if i == 0: + arr_reordered = arr_reordered[[0], :, :, :] + elif i == nwires - 1: + arr_reordered = arr_reordered[:, :, :, [-1]] + + node = Node( + arr_reordered, name=f"mpo_{i}", axis_names=["wL", "p", "p*", "wR"] + ) + nodes.append(node) + + if nwires > 1: + for i in range(nwires - 1): + nodes[i][3] ^ nodes[i + 1][0] + + out_edges = [n[2] for n in nodes] + in_edges = [n[1] for n in nodes] + ignore_edges = [nodes[0][0], nodes[-1][3]] + else: # MPS + for i in range(nwires): + B_obj = tenpy_obj.get_B(i) + arr = B_obj.to_ndarray() + labels = B_obj.get_leg_labels() + vL_idx = labels.index("vL") + p_idx = labels.index("p") + vR_idx = labels.index("vR") + arr_reordered = arr.transpose((vL_idx, p_idx, vR_idx)) + node = Node(arr_reordered, name=f"mps_{i}", axis_names=["vL", "p", "vR"]) + nodes.append(node) + + if nwires > 1: + for i in range(nwires - 1): + nodes[i][2] ^ nodes[i + 1][0] + + out_edges = [n[1] for n in nodes] + in_edges = [] + ignore_edges = [nodes[0][0], nodes[-1][2]] + + qop = quantum_constructor(out_edges, in_edges, [], ignore_edges) - for i in range(nwires - 1): - connect(mpo[i][1], mpo[i + 1][0]) - # TODO(@refraction-ray): whether in and out edge is in the correct order require further check - qop = quantum_constructor( - [mpo[i][-1] for i in range(nwires)], # out_edges - [mpo[i][-2] for i in range(nwires)], # in_edges - [], - [mpo[0][0], mpo[-1][1]], # ignore_edges - ) return qop +def qop2tenpy(qop: QuOperator) -> Any: + """ + Convert TensorCircuit QuOperator to MPO or MPS from TeNPy. + + Requirements: QuOperator must represent valid MPS/MPO structure: + - Linear chain topology with open boundaries only + - MPS: no input edges, consistent virtual bonds, rank-3 or 4(with empty input edges) tensors + - MPO: equal input/output edges, rank-4 tensors + - Cyclic boundary conditions NOT supported + + :param qop: The corresponding state/operator as a QuOperator. + :return: MPO or MPS object from the TeNPy package. + :rtype: Union[tenpy.networks.mpo.MPO, tenpy.networks.mps.MPS] + """ + try: + from tenpy.networks import MPO, MPS, Site + from tenpy.linalg import np_conserved as npc + from tenpy.linalg import LegCharge + except ImportError: + raise ImportError("Please install TeNPy package to use this function.") + + sorted_nodes, is_mps, nwires = extract_tensors_from_qop(qop) + + physical_dim = qop.out_edges[0].dimension if is_mps else qop.in_edges[0].dimension + sites = [Site(LegCharge.from_trivial(physical_dim), "q") for _ in range(nwires)] + + # MPS Conversion + if is_mps: + tensors = [] + for i, node in enumerate(sorted_nodes): + tensor = np.asarray(node.tensor) + if tensor.ndim == 3: + if i == 0: + if tensor.shape[0] > 1: + tensor = tensor[0:1, :, :] + elif i == len(sorted_nodes) - 1: + if tensor.shape[2] > 1: + tensor = tensor[:, :, 0:1] + tensors.append( + npc.Array.from_ndarray( + tensor, + legcharges=[LegCharge.from_trivial(s) for s in tensor.shape], + labels=["vL", "p", "vR"], + ) + ) + + SVs = ( + [np.ones([1])] + + [np.ones(tensors[i].get_leg("vR").ind_len) for i in range(nwires - 1)] + + [np.ones([1])] + ) + return MPS(sites, tensors, SVs, bc="finite") + + # MPO Conversion + raw_tensors = [np.asarray(node.tensor) for node in sorted_nodes] + + if nwires == 1: + chi = 1 + IdL = IdR = 0 + reconstructed_tensors = raw_tensors + else: + chi = max( + raw_tensors[0].shape[3] if raw_tensors[0].ndim > 3 else 1, + raw_tensors[-1].shape[0] if raw_tensors[-1].ndim > 3 else 1, + ) + IdL = 0 + IdR = chi - 1 if chi > 1 else 0 + + reconstructed_tensors = [] + for i, tensor in enumerate(raw_tensors): + if i == 0 and tensor.shape[0] < chi: + new_shape = (chi,) + tensor.shape[1:] + padded_tensor = np.zeros(new_shape, dtype=tensor.dtype) + padded_tensor[IdL, ...] = tensor[0, ...] + reconstructed_tensors.append(padded_tensor) + elif i == nwires - 1 and len(tensor.shape) > 3 and tensor.shape[3] < chi: + new_shape = tensor.shape[:3] + (chi,) + padded_tensor = np.zeros(new_shape, dtype=tensor.dtype) + padded_tensor[..., IdR] = tensor[..., 0] + reconstructed_tensors.append(padded_tensor) + else: + reconstructed_tensors.append(tensor) + + tenpy_Ws = [] + for tensor in reconstructed_tensors: + labels = ["wL", "wR", "p", "p*"] + tensor = np.transpose(tensor, (0, 3, 1, 2)) + tenpy_Ws.append( + npc.Array.from_ndarray( + tensor, + legcharges=[LegCharge.from_trivial(s) for s in tensor.shape], + labels=labels, + ) + ) + + return MPO(sites, tenpy_Ws, bc="finite", IdL=IdL, IdR=IdR) + + def quimb2qop(qb_mpo: Any) -> QuOperator: """ Convert MPO in Quimb package to QuOperator. @@ -1219,6 +1468,151 @@ def quimb2qop(qb_mpo: Any) -> QuOperator: return qop +def qop2quimb(qop: QuOperator) -> Any: + """ + Convert QuOperator to MPO or MPS in Quimb package. + + Requirements: QuOperator must represent valid MPS/MPO structure: + - Linear chain topology with open boundaries only + - MPS: no input edges, consistent virtual bonds between adjacent tensors + - MPO: equal input/output edges, rank-4 tensors + - Edge connectivity: each internal node connected to exactly 2 neighbors + - Cyclic boundary conditions NOT supported + + :param qop: MPO in the form of QuOperator + :return: MPO in the form of Quimb package + :rtype: quimb.tensor.tensor_gen.MatrixProductOperator + """ + try: + import quimb.tensor as qtn + except ImportError: + raise ImportError("Please install Quimb package to use this function.") + + sorted_nodes, is_mps, _ = extract_tensors_from_qop(qop) + + quimb_tensors = [] + node_map = {node: i for i, node in enumerate(sorted_nodes)} + + for i, node in enumerate(sorted_nodes): + tensor_data = node.tensor + inds: List[str] = [] + + for axis, edge in enumerate(node.edges): + if edge in qop.out_edges: + site_index = qop.out_edges.index(edge) + inds.append(f"k{site_index}") + elif edge in qop.in_edges and not is_mps: + site_index = qop.in_edges.index(edge) + inds.append(f"b{site_index}") + elif edge in qop.ignore_edges: + if i == 0: + inds.append("_left_dangling") + elif i == len(sorted_nodes) - 1: + inds.append("_right_dangling") + else: + inds.append(f"_ignore_{i}_{axis}") + else: + neighbor = edge.node1 if edge.node2 is node else edge.node2 + if neighbor in node_map: + j = node_map[neighbor] + left, right = min(i, j), max(i, j) + inds.append(f"v{left}_{right}") + else: + inds.append(f"_unconnected_{i}_{axis}") + + quimb_tensors.append(qtn.Tensor(tensor_data, inds=inds, tags=f"I{i}")) + + tn = qtn.TensorNetwork(quimb_tensors) + + if is_mps: + return tn.as_network(qtn.MatrixProductState) + else: + return tn.as_network(qtn.MatrixProductOperator) + + +def tn2qop(tn_obj: Any) -> QuOperator: + """ + Convert MPO in TensorNetwork package to QuOperator. + + :param tn_mpo: MPO in the form of TensorNetwork package + :type tn_mpo: ``tn.matrixproductstates.mpo.*`` + :return: MPO in the form of QuOperator + :rtype: QuOperator + """ + tn_tensors = tn_obj.tensors + nwires = len(tn_tensors) + + if nwires == 0: + return quantum_constructor([], [], []) + + is_mps = all(len(t.shape) <= 3 for t in tn_tensors) + + nodes = [] + for i in range(nwires): + nodes.append(Node(tn_tensors[i], name=f"tensor_{i}")) + + if is_mps: + for i in range(nwires - 1): + connect(nodes[i][-1], nodes[i + 1][0]) + + out_edges = [] + for i, node in enumerate(nodes): + if len(node.edges) == 2: + physical_edge = next(e for e in node.edges if e.is_dangling()) + out_edges.append(physical_edge) + else: + out_edges.append(node[1]) + + in_edges = [] + + ignore_edges = [] + left_dangling = next( + (e for e in nodes[0].edges if e.is_dangling() and e not in out_edges), None + ) + right_dangling = next( + (e for e in nodes[-1].edges if e.is_dangling() and e not in out_edges), None + ) + + if left_dangling: + ignore_edges.append(left_dangling) + if right_dangling: + ignore_edges.append(right_dangling) + + else: + for i in range(nwires - 1): + connect(nodes[i][1], nodes[i + 1][0]) + + out_edges = [nodes[i][-1] for i in range(nwires)] + in_edges = [nodes[i][-2] for i in range(nwires)] + ignore_edges = [nodes[0][0], nodes[-1][1]] + + qop = quantum_constructor( + out_edges, + in_edges, + [], + ignore_edges, + ) + return qop + + +def qop2tn(qop: QuOperator) -> Any: + """ + Convert QuOperator back to MPO or MPS in TensorNetwork package. + + :param qop: MPO or MPS in the form of QuOperator, param in docstring + :return: MPO or MPS in the form of TensorNetwork + :rtype: Union[tn.matrixproductstates.MPO, tn.matrixproductstates.MPS] + """ + sorted_nodes, is_mps, _ = extract_tensors_from_qop(qop) + + tensors = [node.tensor for node in sorted_nodes] + + if is_mps: + return tn.FiniteMPS(tensors, canonicalize=False) + else: + return tn.matrixproductstates.mpo.FiniteMPO(tensors) + + # TODO(@refraction-ray): Z2 analogy or more general analogies for the following u1 functions diff --git a/tests/test_quantum.py b/tests/test_quantum.py index 633b3e2c..816f6dba 100644 --- a/tests/test_quantum.py +++ b/tests/test_quantum.py @@ -398,30 +398,245 @@ def test_negativity(backend, highp): @pytest.mark.parametrize("backend", [lf("npb"), lf("tfb"), lf("jaxb")]) -def test_tn2qop(backend): - nwires = 6 - dtype = np.complex64 - # only obc is supported, even if you supply nwires Jx terms - Jx = np.array([1.0 for _ in range(nwires - 1)]) # strength of xx interaction (OBC) - Bz = np.array([-1.0 for _ in range(nwires)]) # strength of transverse field - tn_mpo = tn.matrixproductstates.mpo.FiniteTFI(Jx, Bz, dtype=dtype) - qu_mpo = tc.quantum.tn2qop(tn_mpo) - h1 = qu_mpo.eval_matrix() +def test_tenpy2qop(backend): + """ + Tests the conversion from TeNPy MPO/MPS to TensorCircuit QuOperator. + This test verifies the numerical correctness against a precisely matched + Hamiltonian and state vector. + """ + try: + from tenpy.models.tf_ising import TFIChain + from tenpy.networks.mps import MPS + except ImportError: + pytest.skip("TeNPy is not installed") + + nwires = 4 + Jx = 1.0 + Bz = -1.0 + model_params = {"L": nwires, "J": Jx, "g": Bz, "bc_MPS": "finite"} + model = TFIChain(model_params) + + # MPO conversion + tenpy_mpo = model.H_MPO + qu_mpo = tc.quantum.tenpy2qop(tenpy_mpo) + h_actual = qu_mpo.eval_matrix() g = tc.templates.graphs.Line1D(nwires, pbc=False) - h2 = tc.quantum.heisenberg_hamiltonian( - g, hzz=0, hxx=1, hyy=0, hz=1, hx=0, hy=0, sparse=False, numpy=True + h_expected = tc.quantum.heisenberg_hamiltonian( + g, hzz=0, hxx=-Jx, hyy=0, hz=Bz, hx=0, hy=0, sparse=False, numpy=True ) - np.testing.assert_allclose(h1, h2, atol=1e-5) + np.testing.assert_allclose(h_actual, h_expected, atol=1e-5) + + # MPS conversion + psi_tenpy = MPS.from_product_state( + model.lat.mps_sites(), [0, 1] * (nwires // 2), bc=model.lat.bc_MPS + ) + qu_mps = tc.quantum.tenpy2qop(psi_tenpy) + psi_actual = np.reshape(qu_mps.eval_matrix(), -1) + psi_expected = np.zeros(2**nwires) + psi_expected[0b1010] = 1.0 + np.testing.assert_allclose(psi_actual, psi_expected, atol=1e-5) + + +@pytest.mark.parametrize("backend", [lf("npb"), lf("tfb"), lf("jaxb")]) +def test_qop2tenpy(backend): + try: + from tenpy.networks.mps import MPS + except ImportError: + pytest.skip("TeNPy is not installed") + # MPS conversion + nwires_mps = 4 + chi_mps = 2 + phys_dim = 2 + mps_tensors = ( + [np.random.rand(1, phys_dim, chi_mps)] + + [np.random.rand(chi_mps, phys_dim, chi_mps) for _ in range(nwires_mps - 2)] + + [np.random.rand(chi_mps, phys_dim, 1)] + ) + mps_nodes = [tn.Node(t) for t in mps_tensors] + for i in range(nwires_mps - 1): + mps_nodes[i][2] ^ mps_nodes[i + 1][0] + qop_mps = tc.QuOperator( + out_edges=[n[1] for n in mps_nodes], + in_edges=[], + ref_nodes=[], + ignore_edges=[mps_nodes[0][0], mps_nodes[-1][2]], + ) + tenpy_mps = tc.quantum.qop2tenpy(qop_mps) + mat = qop_mps.eval_matrix() + vec_from_qop = np.ravel(mat) + full_wavefunction_tensor = tenpy_mps.get_theta(0, tenpy_mps.L) + vec_from_tenpy = np.ravel(full_wavefunction_tensor.to_ndarray()) + np.testing.assert_allclose(vec_from_qop, vec_from_tenpy, atol=1e-5) + + # MPO conversion + nwires_mpo = 4 + chi_mpo = 3 + t_left = np.random.rand(1, 2, 2, chi_mpo) + t_bulk = [np.random.rand(chi_mpo, 2, 2, chi_mpo) for _ in range(nwires_mpo - 2)] + t_right = np.random.rand(chi_mpo, 2, 2, 1) + mpo_tensors = [t_left] + t_bulk + [t_right] + mpo_nodes = [tn.Node(t) for t in mpo_tensors] + for i in range(nwires_mpo - 1): + mpo_nodes[i][3] ^ mpo_nodes[i + 1][0] + qop_mpo = tc.QuOperator( + [n[2] for n in mpo_nodes], + [n[1] for n in mpo_nodes], + mpo_nodes, + [mpo_nodes[0][0], mpo_nodes[-1][3]], + ) + tenpy_mpo = tc.quantum.qop2tenpy(qop_mpo) + mat_from_qop = qop_mpo.eval_matrix() + + L = tenpy_mpo.L + sites = tenpy_mpo.sites + num_tests = 3 + + for _ in range(num_tests): + random_config = [np.random.randint(0, 2) for _ in range(L)] + psi = MPS.from_product_state(sites, random_config, bc="finite") + + exp_tenpy = tenpy_mpo.expectation_value(psi) + + psi_vector = psi.get_theta(0, L).to_ndarray().flatten() + exp_qop = np.real(np.dot(np.conj(psi_vector), np.dot(mat_from_qop, psi_vector))) + + np.testing.assert_allclose(exp_tenpy, exp_qop, atol=1e-5) @pytest.mark.parametrize("backend", [lf("npb"), lf("tfb"), lf("jaxb")]) -def test_qb2qop(backend): +def test_tenpy_roundtrip(backend): try: - import quimb + from tenpy.networks.mps import MPS + from tenpy.networks.mpo import MPO + from tenpy.networks.site import Site + from tenpy.linalg.charges import LegCharge + from tenpy.linalg import np_conserved as npc + except ImportError: + pytest.skip("TeNPy is not installed") + + tc.set_backend(backend) + + # MPO roundtrip test + nwires_mpo = 3 + chi_mpo = 4 + phys_dim = 2 + sites = [Site(LegCharge.from_trivial(phys_dim), "q") for _ in range(nwires_mpo)] + + t_left = np.random.rand(1, chi_mpo, phys_dim, phys_dim) + Ws = [ + npc.Array.from_ndarray( + t_left, + labels=["wL", "wR", "p", "p*"], + legcharges=[LegCharge.from_trivial(s) for s in t_left.shape], + ) + ] + for _ in range(nwires_mpo - 2): + t_bulk = np.random.rand(chi_mpo, chi_mpo, phys_dim, phys_dim) + Ws.append( + npc.Array.from_ndarray( + t_bulk, + labels=["wL", "wR", "p", "p*"], + legcharges=[LegCharge.from_trivial(s) for s in t_bulk.shape], + ) + ) + t_right = np.random.rand(chi_mpo, 1, phys_dim, phys_dim) + Ws.append( + npc.Array.from_ndarray( + t_right, + labels=["wL", "wR", "p", "p*"], + legcharges=[LegCharge.from_trivial(s) for s in t_right.shape], + ) + ) + mpo_original = MPO(sites, Ws, IdL=0, IdR=chi_mpo - 1) + + qop_mpo = tc.quantum.tenpy2qop(mpo_original) + mpo_roundtrip = tc.quantum.qop2tenpy(qop_mpo) + + assert mpo_original.L == mpo_roundtrip.L + + IdL_rt = mpo_roundtrip.IdL + IdR_rt = mpo_roundtrip.IdR + if isinstance(IdL_rt, (list, np.ndarray)): + IdL_rt = IdL_rt[0] if len(IdL_rt) > 0 else 0 + if isinstance(IdR_rt, (list, np.ndarray)): + IdR_rt = IdR_rt[0] if len(IdR_rt) > 0 else 0 + + for i in range(mpo_original.L): + tensor_orig = mpo_original.get_W(i).to_ndarray() + tensor_rt = mpo_roundtrip.get_W(i).to_ndarray() + + if i == 0: + if tensor_orig.shape[0] < tensor_rt.shape[0]: + tensor_rt_effective = tensor_rt[IdL_rt : IdL_rt + 1, ...] + else: + tensor_rt_effective = tensor_rt + np.testing.assert_allclose(tensor_orig, tensor_rt_effective, atol=1e-5) + elif i == mpo_original.L - 1: + if tensor_orig.shape[1] < tensor_rt.shape[1]: + tensor_rt_effective = tensor_rt[..., IdR_rt : IdR_rt + 1, :, :] + else: + tensor_rt_effective = tensor_rt + np.testing.assert_allclose(tensor_orig, tensor_rt_effective, atol=1e-5) + else: + np.testing.assert_allclose(tensor_orig, tensor_rt, atol=1e-5) + + # MPS roundtrip test + nwires_mps = 4 + chi_mps = 5 + sites_mps = [Site(LegCharge.from_trivial(phys_dim), "q") for _ in range(nwires_mps)] + + b_left = np.random.rand(1, phys_dim, chi_mps) + Bs = [ + npc.Array.from_ndarray( + b_left, + labels=["vL", "p", "vR"], + legcharges=[LegCharge.from_trivial(s) for s in b_left.shape], + ) + ] + for _ in range(nwires_mps - 2): + b_bulk = np.random.rand(chi_mps, phys_dim, chi_mps) + Bs.append( + npc.Array.from_ndarray( + b_bulk, + labels=["vL", "p", "vR"], + legcharges=[LegCharge.from_trivial(s) for s in b_bulk.shape], + ) + ) + b_right = np.random.rand(chi_mps, phys_dim, 1) + Bs.append( + npc.Array.from_ndarray( + b_right, + labels=["vL", "p", "vR"], + legcharges=[LegCharge.from_trivial(s) for s in b_right.shape], + ) + ) + + SVs = [np.ones([1])] + for B in Bs[:-1]: + sv_dim = B.get_leg("vR").ind_len + SVs.append(np.ones(sv_dim)) + SVs.append(np.ones([1])) + + mps_original = MPS(sites_mps, Bs, SVs) + + qop_mps = tc.quantum.tenpy2qop(mps_original) + mps_roundtrip = tc.quantum.qop2tenpy(qop_mps) + + assert mps_original.L == mps_roundtrip.L + for i in range(mps_original.L): + tensor_orig = mps_original.get_B(i, form="B").to_ndarray() + tensor_rt = mps_roundtrip.get_B(i, form="B").to_ndarray() + np.testing.assert_allclose(tensor_orig, tensor_rt, atol=1e-5) + + +@pytest.mark.parametrize("backend", [lf("npb"), lf("tfb"), lf("jaxb")]) +def test_quimb2qop(backend): + try: + import quimb.tensor.tensor_builder as qtb except ImportError: pytest.skip("quimb is not installed") nwires = 6 - qb_mpo = quimb.tensor.tensor_builder.MPO_ham_ising(nwires, 4, 2, cyclic=True) + qb_mpo = qtb.MPO_ham_ising(nwires, 4, 2, cyclic=True) qu_mpo = tc.quantum.quimb2qop(qb_mpo) h1 = qu_mpo.eval_matrix() g = tc.templates.graphs.Line1D(nwires, pbc=True) @@ -431,8 +646,7 @@ def test_qb2qop(backend): np.testing.assert_allclose(h1, h2, atol=1e-5) # in out edge order test - builder = quimb.tensor.tensor_builder.SpinHam1D() - # new version quimb breaking API change: SpinHam1D -> SpinHam + builder = qtb.SpinHam1D() builder += 1, "Y" builder += 1, "X" H = builder.build_mpo(3) @@ -444,15 +658,288 @@ def test_qb2qop(backend): ) np.testing.assert_allclose(m1, m2, atol=1e-5) - # test mps case - - s1 = quimb.tensor.tensor_builder.MPS_rand_state(3, 4) + s1 = qtb.MPS_rand_state(3, 4) s2 = tc.quantum.quimb2qop(s1) m1 = s1.to_dense() m2 = s2.eval_matrix() np.testing.assert_allclose(m1, m2, atol=1e-5) +@pytest.mark.parametrize("backend", [lf("npb"), lf("tfb"), lf("jaxb")]) +def test_qop2quimb(backend): + try: + import quimb.tensor as qtn # pylint: disable=unused-import + except ImportError: + pytest.skip("quimb is not installed") + + # MPO Conversion + nwires_mpo = 4 + chi_mpo = 3 + phys_dim = 2 + + t_left = np.random.rand(1, phys_dim, phys_dim, chi_mpo) + t_bulk = [ + np.random.rand(chi_mpo, phys_dim, phys_dim, chi_mpo) + for _ in range(nwires_mpo - 2) + ] + t_right = np.random.rand(chi_mpo, phys_dim, phys_dim, 1) + mpo_tensors = [t_left] + t_bulk + [t_right] + mpo_nodes = [tn.Node(t) for t in mpo_tensors] + + for i in range(nwires_mpo - 1): + mpo_nodes[i][3] ^ mpo_nodes[i + 1][0] + + qop_mpo = tc.QuOperator( + out_edges=[n[2] for n in mpo_nodes], + in_edges=[n[1] for n in mpo_nodes], + ref_nodes=mpo_nodes, + ignore_edges=[mpo_nodes[0][0], mpo_nodes[-1][3]], + ) + + quimb_mpo = tc.quantum.qop2quimb(qop_mpo) + + mat_from_qop = qop_mpo.eval_matrix() + + ket_inds_mpo = [f"k{i}" for i in range(nwires_mpo)] + bra_inds_mpo = [f"b{i}" for i in range(nwires_mpo)] + mat_from_quimb = quimb_mpo.to_dense(ket_inds_mpo, bra_inds_mpo) + + np.testing.assert_allclose(mat_from_qop, mat_from_quimb, atol=1e-5) + + # MPS Conversion + nwires_mps = 5 + chi_mps = 8 + + mps_tensors = ( + [np.random.rand(1, phys_dim, chi_mps)] + + [np.random.rand(chi_mps, phys_dim, chi_mps) for _ in range(nwires_mps - 2)] + + [np.random.rand(chi_mps, phys_dim, 1)] + ) + mps_nodes = [tn.Node(t) for t in mps_tensors] + + for i in range(nwires_mps - 1): + mps_nodes[i][2] ^ mps_nodes[i + 1][0] + + qop_mps = tc.QuOperator( + out_edges=[n[1] for n in mps_nodes], + in_edges=[], + ref_nodes=mps_nodes, + ignore_edges=[mps_nodes[0][0], mps_nodes[-1][2]], + ) + + quimb_mps = tc.quantum.qop2quimb(qop_mps) + + mat_from_qop = qop_mps.eval_matrix() + vec_from_qop = np.ravel(mat_from_qop) + + ket_inds_mps = [f"k{i}" for i in range(nwires_mps)] + vec_from_quimb = np.ravel(quimb_mps.to_dense(ket_inds_mps)) + + np.testing.assert_allclose(vec_from_qop, vec_from_quimb, atol=1e-5) + + +@pytest.mark.parametrize("backend", [lf("npb"), lf("tfb"), lf("jaxb")]) +def test_quimb_roundtrip(backend): + try: + import quimb.tensor as qtn + except ImportError: + pytest.skip("quimb is not installed") + # MPO roundtrip test + nwires_mpo = 4 + mpo_original = qtn.MPO_ham_ising(nwires_mpo) + qop_mpo = tc.quantum.quimb2qop(mpo_original) + mpo_roundtrip = tc.quantum.qop2quimb(qop_mpo) + ket_inds = [f"k{i}" for i in range(nwires_mpo)] + bra_inds = [f"b{i}" for i in range(nwires_mpo)] + mat_original = mpo_original.to_dense(ket_inds, bra_inds) + mat_roundtrip = mpo_roundtrip.to_dense(ket_inds, bra_inds) + np.testing.assert_allclose(mat_original, mat_roundtrip, atol=1e-5) + + # MPS roundtrip test + nwires_mps = 5 + bond_dim = 8 + mps_original = qtn.MPS_rand_state(nwires_mps, bond_dim) + mps_original.normalize() + qop_mps = tc.quantum.quimb2qop(mps_original) + mps_roundtrip = tc.quantum.qop2quimb(qop_mps) + ket_inds_mps = [f"k{i}" for i in range(nwires_mps)] + vec_original = np.ravel(mps_original.to_dense(ket_inds_mps)) + vec_roundtrip = np.ravel(mps_roundtrip.to_dense(ket_inds_mps)) + np.testing.assert_allclose( + abs(np.vdot(vec_original, vec_roundtrip)), 1.0, atol=1e-5 + ) + + +@pytest.mark.parametrize("backend", [lf("npb"), lf("tfb"), lf("jaxb")]) +def test_tn2qop(backend): + nwires = 6 + dtype = np.complex64 + + # MPO test + Jx = np.array([1.0 for _ in range(nwires - 1)]) + Bz = np.array([-1.0 for _ in range(nwires)]) + tn_mpo = tn.matrixproductstates.mpo.FiniteTFI(Jx, Bz, dtype=dtype) + qu_mpo = tc.quantum.tn2qop(tn_mpo) + h1 = qu_mpo.eval_matrix() + g = tc.templates.graphs.Line1D(nwires, pbc=False) + h2 = tc.quantum.heisenberg_hamiltonian( + g, hzz=0, hxx=1, hyy=0, hz=1, hx=0, hy=0, sparse=False, numpy=True + ) + np.testing.assert_allclose(h1, h2, atol=1e-5) + + # MPS test + mps_tensors = [] + bond_dim = 1 + phys_dim = 2 + + first_tensor = np.zeros((phys_dim, bond_dim), dtype=dtype) + first_tensor[0, 0] = 1.0 + mps_tensors.append(first_tensor) + + for _ in range(nwires - 2): + middle_tensor = np.zeros((bond_dim, phys_dim, bond_dim), dtype=dtype) + middle_tensor[0, 0, 0] = 1.0 + mps_tensors.append(middle_tensor) + + if nwires > 1: + last_tensor = np.zeros((bond_dim, phys_dim), dtype=dtype) + last_tensor[0, 0] = 1.0 + mps_tensors.append(last_tensor) + + mps_tensors = [np.array(t, dtype=dtype) for t in mps_tensors] + + tn_mps = tn.FiniteMPS(mps_tensors, canonicalize=False) + + # turn to QuOperator + qu_mps = tc.quantum.tn2qop(tn_mps) + assert qu_mps.is_vector() + assert len(qu_mps.out_edges) == nwires + assert len(qu_mps.in_edges) == 0 + state_tensor = qu_mps.eval() + state_vector = np.ravel(state_tensor) + norm = tc.backend.norm(state_vector) + np.testing.assert_allclose(norm, 1.0, atol=1e-5) + expected_state = np.zeros(2**nwires, dtype=dtype) + expected_state[0] = 1.0 + + np.testing.assert_allclose(state_vector, expected_state, atol=1e-5) + + +@pytest.mark.parametrize("backend", [lf("npb"), lf("tfb"), lf("jaxb")]) +def test_qop2tn(backend): + nwires = 4 + chi = 5 + phys_dim = 2 + dtype = np.complex64 + + # MPO conversion test + tensors_mpo_tc = [np.random.rand(1, chi, phys_dim, phys_dim).astype(dtype)] + for _ in range(nwires - 2): + tensors_mpo_tc.append( + np.random.rand(chi, chi, phys_dim, phys_dim).astype(dtype) + ) + tensors_mpo_tc.append(np.random.rand(chi, 1, phys_dim, phys_dim).astype(dtype)) + + nodes_mpo_tc = [tn.Node(t) for t in tensors_mpo_tc] + for i in range(nwires - 1): + nodes_mpo_tc[i][1] ^ nodes_mpo_tc[i + 1][0] + + qop_mpo = tc.QuOperator( + out_edges=[n[3] for n in nodes_mpo_tc], + in_edges=[n[2] for n in nodes_mpo_tc], + ref_nodes=nodes_mpo_tc, + ignore_edges=[nodes_mpo_tc[0][0], nodes_mpo_tc[-1][1]], + ) + + tn_mpo_rt = tc.quantum.qop2tn(qop_mpo) + qop_mpo_rt = tc.quantum.tn2qop(tn_mpo_rt) + + mat_original = qop_mpo.eval_matrix() + mat_roundtrip = qop_mpo_rt.eval_matrix() + np.testing.assert_allclose(mat_original, mat_roundtrip, atol=1e-5) + + # MPS conversion test + tensors_mps_tc = [np.random.rand(1, phys_dim, chi).astype(dtype)] + for _ in range(nwires - 2): + tensors_mps_tc.append(np.random.rand(chi, phys_dim, chi).astype(dtype)) + tensors_mps_tc.append(np.random.rand(chi, phys_dim, 1).astype(dtype)) + + nodes_mps_tc = [tn.Node(t) for t in tensors_mps_tc] + for i in range(nwires - 1): + nodes_mps_tc[i][2] ^ nodes_mps_tc[i + 1][0] + + qop_mps = tc.QuOperator( + out_edges=[n[1] for n in nodes_mps_tc], + in_edges=[], + ref_nodes=nodes_mps_tc, + ignore_edges=[nodes_mps_tc[0][0], nodes_mps_tc[-1][2]], + ) + + tn_mps_rt = tc.quantum.qop2tn(qop_mps) + qop_mps_rt = tc.quantum.tn2qop(tn_mps_rt) + + vec_original = np.ravel(qop_mps.eval_matrix()) + vec_roundtrip = np.ravel(qop_mps_rt.eval_matrix()) + + vec_original = vec_original / np.linalg.norm(vec_original) + vec_roundtrip = vec_roundtrip / np.linalg.norm(vec_roundtrip) + + overlap = abs(np.vdot(vec_original, vec_roundtrip)) + np.testing.assert_allclose(overlap, 1.0, atol=1e-5) + + +@pytest.mark.parametrize("backend", [lf("npb"), lf("tfb"), lf("jaxb")]) +def test_tn_roundtrip(backend): + # MPO roundtrip test + nwires_mpo = 4 + dtype = np.complex64 + + Jx = np.array([1.0 for _ in range(nwires_mpo - 1)]) + Bz = np.array([-1.0 for _ in range(nwires_mpo)]) + tn_mpo_original = tn.matrixproductstates.mpo.FiniteTFI(Jx, Bz, dtype=dtype) + + qop_mpo = tc.quantum.tn2qop(tn_mpo_original) + tn_mpo_roundtrip = tc.quantum.qop2tn(qop_mpo) + + qop_mpo_original = tc.quantum.tn2qop(tn_mpo_original) + qop_mpo_roundtrip = tc.quantum.tn2qop(tn_mpo_roundtrip) + + mat_original = qop_mpo_original.eval_matrix() + mat_roundtrip = qop_mpo_roundtrip.eval_matrix() + + np.testing.assert_allclose(mat_original, mat_roundtrip, atol=1e-5) + + # MPS roundtrip test + nwires_mps = 4 + chi_mps = 2 + phys_dim = 2 + + mps_tensors = ( + [np.random.rand(1, phys_dim, chi_mps).astype(dtype)] + + [ + np.random.rand(chi_mps, phys_dim, chi_mps).astype(dtype) + for _ in range(nwires_mps - 2) + ] + + [np.random.rand(chi_mps, phys_dim, 1).astype(dtype)] + ) + + tn_mps_original = tn.FiniteMPS(mps_tensors, canonicalize=False) + qop_mps = tc.quantum.tn2qop(tn_mps_original) + tn_mps_roundtrip = tc.quantum.qop2tn(qop_mps) + + qop_mps_original = tc.quantum.tn2qop(tn_mps_original) + qop_mps_roundtrip = tc.quantum.tn2qop(tn_mps_roundtrip) + + vec_original = np.ravel(qop_mps_original.eval_matrix()) + vec_roundtrip = np.ravel(qop_mps_roundtrip.eval_matrix()) + + vec_original = vec_original / np.linalg.norm(vec_original) + vec_roundtrip = vec_roundtrip / np.linalg.norm(vec_roundtrip) + + overlap = abs(np.vdot(vec_original, vec_roundtrip)) + np.testing.assert_allclose(overlap, 1.0, atol=1e-5) + + @pytest.mark.parametrize("backend", [lf("npb"), lf("tfb"), lf("jaxb")]) def test_counts_2(backend): z0 = tc.backend.convert_to_tensor(np.array([0.1, 0, -0.3, 0]))