From 6f4233bbec4c0077731f043e282618a15492ca85 Mon Sep 17 00:00:00 2001 From: Synge Todo Date: Wed, 8 Nov 2023 21:49:53 +0900 Subject: [PATCH] contraction graph --- .gitignore | 3 + src/qailo/__init__.py | 3 +- src/qailo/mps/apply.py | 2 +- src/qailo/mps/mps_c.py | 6 +- src/qailo/mps/norm.py | 4 +- src/qailo/mps/product_state.py | 4 +- src/qailo/mps/state_vector.py | 4 +- src/qailo/mps/type.py | 9 +- src/qailo/mps_p/mps_p.py | 24 ++-- src/qailo/mps_p/product_state.py | 4 +- src/qailo/mps_t/__init__.py | 9 ++ src/qailo/mps_t/mps_t.py | 193 +++++++++++++++++++++++++++++++ src/qailo/mps_t/product_state.py | 19 +++ test/mps/test_apply.py | 66 +++++++---- test/mps/test_canonical.py | 79 ++++++------- test/mps/test_move.py | 93 ++++++--------- test/mps/test_mps.py | 6 + test/mps_p/test_mps_p.py | 15 --- test/mps_t/test_mps_t.py | 36 ++++++ 19 files changed, 410 insertions(+), 169 deletions(-) create mode 100644 src/qailo/mps_t/__init__.py create mode 100644 src/qailo/mps_t/mps_t.py create mode 100644 src/qailo/mps_t/product_state.py delete mode 100644 test/mps_p/test_mps_p.py create mode 100644 test/mps_t/test_mps_t.py diff --git a/.gitignore b/.gitignore index 12e9e05..e0ae5a4 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,6 @@ build **/*.py[cod] **/__pycache__ **/*.egg-info +**/test_mps_t-generator.pkl +**/test_mps_t-graph.json +**/test_mps_t-tensor.pkl diff --git a/src/qailo/__init__.py b/src/qailo/__init__.py index 6475731..4273535 100644 --- a/src/qailo/__init__.py +++ b/src/qailo/__init__.py @@ -1,4 +1,4 @@ -from . import alg, mps, mps_p, operator, state_vector, util +from . import alg, mps, mps_p, mps_t, operator, state_vector, util from . import operator as op from . import state_vector as sv from ._version import version @@ -8,6 +8,7 @@ alg, mps, mps_p, + mps_t, operator, state_vector, util, diff --git a/src/qailo/mps/apply.py b/src/qailo/mps/apply.py index 57a8f7c..06bd104 100644 --- a/src/qailo/mps/apply.py +++ b/src/qailo/mps/apply.py @@ -9,7 +9,7 @@ def _swap_tensors(m, s, maxdim=None): """ swap neighboring two tensors at s and s+1 """ - assert s in range(0, len(m.tensors) - 1) + assert s in range(0, mps.num_qubits(m) - 1) m._apply_two(swap(), s, maxdim=maxdim) p0, p1 = m.t2q[s], m.t2q[s + 1] m.q2t[p0], m.q2t[p1] = s + 1, s diff --git a/src/qailo/mps/mps_c.py b/src/qailo/mps/mps_c.py index ea490a4..03fb7a1 100644 --- a/src/qailo/mps/mps_c.py +++ b/src/qailo/mps/mps_c.py @@ -4,9 +4,10 @@ from ..operator import type as op from .svd import tensor_svd +from .type import MPS -class MPS_C: +class MPS_C(MPS): """ MPS representation of quantum pure state @@ -29,6 +30,9 @@ def __init__(self, tensors): self.t2q = list(range(n)) self.cp = [0, n - 1] + def _tensor(self, t): + return self.tensors[t] + def _canonicalize(self, p0, p1=None): p1 = p0 if p1 is None else p1 n = len(self.tensors) diff --git a/src/qailo/mps/norm.py b/src/qailo/mps/norm.py index ff3bdc2..2da07e4 100644 --- a/src/qailo/mps/norm.py +++ b/src/qailo/mps/norm.py @@ -6,6 +6,6 @@ def norm(m): A = np.identity(1) for t in range(num_qubits(m)): - A = np.einsum("ij,jkl->ikl", A, m.tensors[t]) - A = np.einsum("ijk,ijl->kl", A, m.tensors[t].conj()) + A = np.einsum("ij,jkl->ikl", A, m._tensor(t)) + A = np.einsum("ijk,ijl->kl", A, m._tensor(t).conj()) return np.sqrt(np.trace(A)) diff --git a/src/qailo/mps/product_state.py b/src/qailo/mps/product_state.py index 7aa2b39..a80556b 100644 --- a/src/qailo/mps/product_state.py +++ b/src/qailo/mps/product_state.py @@ -1,6 +1,6 @@ +from ..dispatch import num_qubits from ..state_vector.state_vector import one as sv_one from ..state_vector.state_vector import zero as sv_zero -from ..state_vector.type import num_qubits from ..state_vector.vector import vector from .mps_c import MPS_C from .svd import tensor_svd @@ -9,7 +9,7 @@ def tensor_decomposition(v, nkeep=None, tol=1e-12): if is_mps(v): - return v.tensors + return [v._tensor(s) for s in range(num_qubits(v))] else: n = num_qubits(v) w = vector(v).reshape((1, 2**n)) diff --git a/src/qailo/mps/state_vector.py b/src/qailo/mps/state_vector.py index ac02734..dd0b5ff 100644 --- a/src/qailo/mps/state_vector.py +++ b/src/qailo/mps/state_vector.py @@ -6,10 +6,10 @@ def state_vector(m): assert is_mps(m) n = num_qubits(m) - v = m.tensors[0] + v = m._tensor(0) for t in range(1, n): ss0 = list(range(t + 1)) + [t + 3] ss1 = [t + 3, t + 1, t + 2] - v = np.einsum(v, ss0, m.tensors[t], ss1) + v = np.einsum(v, ss0, m._tensor(t), ss1) v = v.reshape((2,) * n) return np.einsum(v, m.t2q).reshape((2,) * n + (1,)) diff --git a/src/qailo/mps/type.py b/src/qailo/mps/type.py index e21c188..b96217f 100644 --- a/src/qailo/mps/type.py +++ b/src/qailo/mps/type.py @@ -1,10 +1,15 @@ +class MPS(object): + def __init__(self): + True + + def is_canonical(m): return m._is_canonical() def is_mps(m): - return hasattr(m, "tensors") + return isinstance(m, MPS) def num_qubits(m): - return len(m.tensors) + return len(m.q2t) diff --git a/src/qailo/mps_p/mps_p.py b/src/qailo/mps_p/mps_p.py index 28d598a..70b211c 100644 --- a/src/qailo/mps_p/mps_p.py +++ b/src/qailo/mps_p/mps_p.py @@ -3,11 +3,12 @@ import numpy as np from ..mps.svd import tensor_svd +from ..mps.type import MPS from ..operator import type as op from .projector import compact_projector -class MPS_P: +class MPS_P(MPS): """ MPS representation of quantum pure state @@ -31,7 +32,10 @@ def __init__(self, tensors): self.cp = [0, n - 1] # canonicalization matrices # put sentinels (1x1 identities) at t = 0 and t = n - self.cmat = [np.identity(1)] + [None] * (n - 1) + [np.identity(1)] + self.env = [np.identity(1)] + [None] * (n - 1) + [np.identity(1)] + + def _tensor(self, t): + return self.tensors[t] def _canonicalize(self, p0, p1=None): p1 = p0 if p1 is None else p1 @@ -39,14 +43,14 @@ def _canonicalize(self, p0, p1=None): assert 0 <= p0 and p0 <= p1 and p1 < n if self.cp[0] < p0: for t in range(self.cp[0], p0): - A = np.einsum(self.cmat[t], [0, 3], self.tensors[t], [3, 1, 2]) - _, self.cmat[t + 1] = tensor_svd(A, [[0, 1], [2]], "left") + A = np.einsum(self.env[t], [0, 3], self.tensors[t], [3, 1, 2]) + _, self.env[t + 1] = tensor_svd(A, [[0, 1], [2]], "left") self.cp[0] = p0 self.cp[1] = max(p0, self.cp[1]) if self.cp[1] > p1: for t in range(self.cp[1], p1, -1): - A = np.einsum(self.tensors[t], [0, 1, 3], self.cmat[t + 1], [3, 2]) - self.cmat[t], _ = tensor_svd(A, [[0], [1, 2]], "right") + A = np.einsum(self.tensors[t], [0, 1, 3], self.env[t + 1], [3, 2]) + self.env[t], _ = tensor_svd(A, [[0], [1, 2]], "right") self.cp[1] = p1 def _is_canonical(self): @@ -76,14 +80,14 @@ def _is_canonical(self): for t in range(0, self.cp[0]): A = np.einsum(A, [0, 3], self.tensors[t], [3, 1, 2]) A = np.einsum(A, [2, 3, 1], self.tensors[t].conj(), [2, 3, 0]) - B = np.einsum(self.cmat[t + 1], [2, 1], self.cmat[t + 1].conj(), [2, 0]) + B = np.einsum(self.env[t + 1], [2, 1], self.env[t + 1].conj(), [2, 0]) assert A.shape == B.shape assert np.allclose(A, B) A = np.identity(1) for t in range(n - 1, self.cp[1], -1): A = np.einsum(self.tensors[t], [0, 1, 3], A, [3, 2]) A = np.einsum(self.tensors[t].conj(), [1, 2, 3], A, [0, 2, 3]) - B = np.einsum(self.cmat[t], [0, 2], self.cmat[t].conj(), [1, 2]) + B = np.einsum(self.env[t], [0, 2], self.env[t].conj(), [1, 2]) assert np.allclose(A, B) return True @@ -107,8 +111,8 @@ def _apply_two(self, p, s, maxdim=None, reverse=False): else: t0 = np.einsum(t0, [0, 4, 3], p1, [2, 1, 4]) t1 = np.einsum(t1, [0, 4, 3], p0, [2, 4, 1]) - tt0 = np.einsum(self.cmat[s], [0, 4], t0, [4, 1, 2, 3]) - tt1 = np.einsum(t1, [0, 1, 2, 4], self.cmat[s + 2], [4, 3]) + tt0 = np.einsum(self.env[s], [0, 4], t0, [4, 1, 2, 3]) + tt1 = np.einsum(t1, [0, 1, 2, 4], self.env[s + 2], [4, 3]) _, WLh, WR = compact_projector(tt0, [0, 1, 4, 5], tt1, [5, 4, 2, 3], maxdim) self.tensors[s] = np.einsum(t0, [0, 1, 3, 4], WR, [3, 4, 2]) self.tensors[s + 1] = np.einsum(WLh, [3, 4, 0], t1, [4, 3, 1, 2]) diff --git a/src/qailo/mps_p/product_state.py b/src/qailo/mps_p/product_state.py index 148f540..f53606c 100644 --- a/src/qailo/mps_p/product_state.py +++ b/src/qailo/mps_p/product_state.py @@ -11,9 +11,9 @@ def product_state(states): return MPS_P(tensors) -def zero(n=1, logger=None): +def zero(n=1): return product_state([sv_zero()] * n) -def one(n=1, logger=None): +def one(n=1): return product_state([sv_one()] * n) diff --git a/src/qailo/mps_t/__init__.py b/src/qailo/mps_t/__init__.py new file mode 100644 index 0000000..17bff05 --- /dev/null +++ b/src/qailo/mps_t/__init__.py @@ -0,0 +1,9 @@ +from .mps_t import MPS_T +from .product_state import one, product_state, zero + +__all__ = [ + MPS_T, + one, + product_state, + zero, +] diff --git a/src/qailo/mps_t/mps_t.py b/src/qailo/mps_t/mps_t.py new file mode 100644 index 0000000..f2db236 --- /dev/null +++ b/src/qailo/mps_t/mps_t.py @@ -0,0 +1,193 @@ +from copy import deepcopy + +import numpy as np + +from ..mps.svd import tensor_svd +from ..mps.type import MPS +from ..mps_p.projector import _full_projector +from ..operator import type as op + + +class MPS_T(MPS): + """ + MPS representation of quantum pure state + + shape of tensors: [du, dp, dl] + du: dimension of upper leg (1 for top tensor) + dp: dimension of physical leg (typically 2) + dl: dimension of lower leg (1 for bottom tensor) + + canonical position: cp in range(n) + 0 <= cp(0) <= cp(1) < n + tensors [0...cp(0)-1]: top canonical + tensors [cp(1)+1...n-1]: bottom canonical + """ + + def __init__(self, tensors): + assert isinstance(tensors, list) + n = len(tensors) + self.tpool = [] + self.tcurrent = [] + for t in tensors: + self.tcurrent.append(self._register(deepcopy(t))) + self.gpool = [] + self.q2t = list(range(n)) + self.t2q = list(range(n)) + self.cp = [0, n - 1] + # canonicalization matrices + # put sentinels (1x1 identities) at t = 0 and t = n + self.env = [np.identity(1)] + [None] * (n - 1) + [np.identity(1)] + + def _num_qubits(self): + return len(self.tcurrent) + + def _tensor(self, t): + return self.tpool[self.tcurrent[t]][0] + + def _canonicalize(self, p0, p1=None): + p1 = p0 if p1 is None else p1 + assert 0 <= p0 and p0 <= p1 and p1 < self._num_qubits() + if self.cp[0] < p0: + for t in range(self.cp[0], p0): + A = np.einsum(self.env[t], [0, 3], self._tensor(t), [3, 1, 2]) + _, self.env[t + 1] = tensor_svd(A, [[0, 1], [2]], "left") + self.cp[0] = p0 + self.cp[1] = max(p0, self.cp[1]) + if self.cp[1] > p1: + for t in range(self.cp[1], p1, -1): + A = np.einsum(self._tensor(t), [0, 1, 3], self.env[t + 1], [3, 2]) + self.env[t], _ = tensor_svd(A, [[0], [1, 2]], "right") + self.cp[1] = p1 + + def _is_canonical(self): + # tensor shape + n = len(self.tcurrent) + dims = [] + assert self._tensor(0).shape[0] == 1 + dims.append(self._tensor(0).shape[0]) + for t in range(1, n - 1): + dims.append(self._tensor(t).shape[0]) + assert self._tensor(t).shape[0] == self._tensor(t - 1).shape[2] + assert self._tensor(t).shape[2] == self._tensor(t + 1).shape[0] + assert self._tensor(n - 1).shape[2] == 1 + dims.append(self._tensor(n - 1).shape[0]) + dims.append(self._tensor(n - 1).shape[2]) + + # qubit <-> tensor mapping + for q in range(n): + assert self.t2q[self.q2t[q]] == q + for t in range(n): + assert self.q2t[self.t2q[t]] == t + + # canonicality + assert self.cp[0] in range(n) + assert self.cp[1] in range(n) + A = np.identity(1) + for t in range(0, self.cp[0]): + A = np.einsum(A, [0, 3], self._tensor(t), [3, 1, 2]) + A = np.einsum(A, [2, 3, 1], self._tensor(t).conj(), [2, 3, 0]) + B = np.einsum(self.env[t + 1], [2, 1], self.env[t + 1].conj(), [2, 0]) + assert A.shape == B.shape + assert np.allclose(A, B) + A = np.identity(1) + for t in range(n - 1, self.cp[1], -1): + A = np.einsum(self._tensor(t), [0, 1, 3], A, [3, 2]) + A = np.einsum(self._tensor(t).conj(), [1, 2, 3], A, [0, 2, 3]) + B = np.einsum(self.env[t], [0, 2], self.env[t].conj(), [1, 2]) + assert np.allclose(A, B) + return True + + def _apply_one(self, p, s): + assert op.num_qubits(p) == 1 + pid = self._register(p) + self.tcurrent[s] = self._contract(self.tcurrent[s], [0, 3, 2], pid, [1, 3]) + self.cp[0] = min(self.cp[0], s) + self.cp[1] = max(self.cp[1], s) + + def _apply_two(self, p, s, maxdim=None, reverse=False): + """ + apply 2-qubit operator on neighboring tensors, s and s+1 + """ + self._canonicalize(s, s + 1) + tid0 = self.tcurrent[s] + tid1 = self.tcurrent[s + 1] + p0, p1 = tensor_svd(p, [[0, 2], [1, 3]]) + pid0 = self._register(p0) + pid1 = self._register(p1) + if not reverse: + tid0 = self._contract(tid0, [0, 4, 3], pid0, [1, 4, 2]) + tid1 = self._contract(tid1, [0, 4, 3], pid1, [1, 2, 4]) + else: + tid0 = self._contract(tid0, [0, 4, 3], pid1, [2, 1, 4]) + tid1 = self._contract(tid1, [0, 4, 3], pid0, [2, 4, 1]) + tt0 = np.einsum(self.env[s], [0, 4], self.tpool[tid0][0], [4, 1, 2, 3]) + tt1 = np.einsum(self.tpool[tid1][0], [0, 1, 2, 4], self.env[s + 2], [4, 3]) + _, WLhid, WRid = self._projector(tt0, [0, 1, 4, 5], tt1, [5, 4, 2, 3], maxdim) + self.tcurrent[s] = self._contract(tid0, [0, 1, 3, 4], WRid, [3, 4, 2]) + self.tcurrent[s + 1] = self._contract(WLhid, [3, 4, 0], tid1, [4, 3, 1, 2]) + + def _register(self, tensor): + id = len(self.tpool) + self.tpool.append([tensor, "initial", None]) + return id + + def _contract(self, tid0, ss0, tid1, ss1, ss2=None): + id = len(self.tpool) + if ss2 is None: + t = np.einsum(self.tpool[tid0][0], ss0, self.tpool[tid1][0], ss1) + else: + np.einsum(self.tpool[tid0][0], ss0, self.tpool[tid1][0], ss1, ss2) + self.tpool.append([t, "product", [tid0, tid1]]) + return id + + def _projector(self, t0, ss0, t1, ss1, maxdim=None): + S, d, WLh, WR = _full_projector(t0, ss0, t1, ss1) + d = d if maxdim is None else min(d, maxdim) + gid = len(self.gpool) + self.gpool.append([S, d, WLh, WR]) + lid = len(self.tpool) + shape = WLh.shape + WLh = WLh.reshape((np.prod(shape[:-1]), shape[-1])) + WLh = WLh[:, :d].reshape(shape[:-1] + (d,)) + self.tpool.append([WLh, "squeezer", [gid]]) + rid = len(self.tpool) + shape = WR.shape + WR = WR.reshape((np.prod(shape[:-1]), shape[-1])) + WR = WR[:, :d].reshape(shape[:-1] + (d,)) + self.tpool.append([WR, "squeezer", [gid]]) + return gid, lid, rid + + def _dump(self, prefix): + import json + import pickle + + dic = {"prefix": prefix} + tlist = [] + with open(f"{prefix}-tensor.pkl", "wb") as f: + for id, tp in enumerate(self.tpool): + m = {} + m["id"] = id + m["shape"] = tp[0].shape + m["type"] = tp[1] + m["from"] = tp[2] + tlist.append(m) + if tp[1] == "initial": + pickle.dump(tp[0], f) + dic["tensor"] = tlist + glist = [] + with open(f"{prefix}-generator.pkl", "wb") as f: + for id, gp in enumerate(self.gpool): + m = {} + m["id"] = id + m["d"] = gp[1] + m["shape S"] = gp[0].shape + m["shape L"] = gp[2].shape + m["shape R"] = gp[3].shape + glist.append(m) + pickle.dump(gp[0], f) + pickle.dump(gp[2], f) + pickle.dump(gp[3], f) + dic["generator"] = glist + + with open(prefix + "-graph.json", mode="w") as f: + json.dump(dic, f, indent=2) diff --git a/src/qailo/mps_t/product_state.py b/src/qailo/mps_t/product_state.py new file mode 100644 index 0000000..adfa10d --- /dev/null +++ b/src/qailo/mps_t/product_state.py @@ -0,0 +1,19 @@ +from ..mps.product_state import tensor_decomposition +from ..state_vector.state_vector import one as sv_one +from ..state_vector.state_vector import zero as sv_zero +from .mps_t import MPS_T + + +def product_state(states): + tensors = [] + for s in states: + tensors = tensors + tensor_decomposition(s) + return MPS_T(tensors) + + +def zero(n=1): + return product_state([sv_zero()] * n) + + +def one(n=1): + return product_state([sv_one()] * n) diff --git a/test/mps/test_apply.py b/test/mps/test_apply.py index 84b4814..057d9f1 100644 --- a/test/mps/test_apply.py +++ b/test/mps/test_apply.py @@ -2,18 +2,18 @@ import qailo as q from pytest import approx -from qailo.mps.mps_c import MPS_C -from qailo.mps_p.mps_p import MPS_P -def apply(op, m0, m1, m2, m3, v, seq, pos, maxdim=None): +def apply(op, m0, m1, m2, m3, m4, m5, v, seq, pos, maxdim=None): m0 = q.mps.apply(m0, op, pos) m1 = q.mps.apply(m1, op, pos) - m2 = q.mps.apply(m2, op, pos, maxdim) + m2 = q.mps.apply(m2, op, pos) m3 = q.mps.apply(m3, op, pos, maxdim) + m4 = q.mps.apply(m4, op, pos, maxdim) + m5 = q.mps.apply(m5, op, pos, maxdim) v = q.sv.apply(v, op, pos) seq.append([op, pos]) - return m0, m1, m2, m3, v, seq + return m0, m1, m2, m3, m4, m5, v, seq def test_apply(): @@ -21,17 +21,21 @@ def test_apply(): p = 64 maxdim = 2 - m0 = q.mps.zero(n, MPS_C) - m1 = q.mps.zero(n, MPS_P) - m2 = q.mps.zero(n, MPS_C) - m3 = q.mps.zero(n, MPS_P) + m0 = q.mps.zero(n, q.mps.MPS_C) + m1 = q.mps.zero(n, q.mps_p.MPS_P) + m2 = q.mps.zero(n, q.mps_t.MPS_T) + m3 = q.mps.zero(n, q.mps.MPS_C) + m4 = q.mps.zero(n, q.mps_p.MPS_P) + m5 = q.mps.zero(n, q.mps_t.MPS_T) v = q.sv.zero(n) seq = [] i = 4 j = 0 print("apply cz on {} and {}".format(i, j)) - m0, m1, m2, m3, v, seq = apply(q.op.cz(), m0, m1, m2, m3, v, seq, [i, j], maxdim) + m0, m1, m2, m3, m4, m5, v, seq = apply( + q.op.cz(), m0, m1, m2, m3, m4, m5, v, seq, [i, j], maxdim + ) for _ in range(p): i = random.randrange(n) @@ -40,55 +44,67 @@ def test_apply(): t = random.randrange(3) if t == 0: print("apply h on {}".format(i)) - m0, m1, m2, m3, v, seq = apply( - q.op.h(), m0, m1, m2, m3, v, seq, [i], maxdim + m0, m1, m2, m3, m4, m5, v, seq = apply( + q.op.h(), m0, m1, m2, m3, m4, m5, v, seq, [i], maxdim ) elif t == 1: print("apply x on {}".format(i)) - m0, m1, m2, m3, v, seq = apply( - q.op.s(), m0, m1, m2, m3, v, seq, [i], maxdim + m0, m1, m2, m3, m4, m5, v, seq = apply( + q.op.s(), m0, m1, m2, m3, m4, m5, v, seq, [i], maxdim ) elif t == 2: print("apply s on {}".format(i)) - m0, m1, m2, m3, v, seq = apply( - q.op.t(), m0, m1, m2, m3, v, seq, [i], maxdim + m0, m1, m2, m3, m4, m5, v, seq = apply( + q.op.t(), m0, m1, m2, m3, m4, m5, v, seq, [i], maxdim ) else: t = random.randrange(2) if t == 0: print("apply cx on {} and {}".format(i, j)) - m0, m1, m2, m3, v, seq = apply( - q.op.cx(), m0, m1, m2, m3, v, seq, [i, j], maxdim + m0, m1, m2, m3, m4, m5, v, seq = apply( + q.op.cx(), m0, m1, m2, m3, m4, m5, v, seq, [i, j], maxdim ) elif t == 1: print("apply cz on {} and {}".format(i, j)) - m0, m1, m2, m3, v, seq = apply( - q.op.cz(), m0, m1, m2, m3, v, seq, [i, j], maxdim + m0, m1, m2, m3, m4, m5, v, seq = apply( + q.op.cz(), m0, m1, m2, m3, m4, m5, v, seq, [i, j], maxdim ) m0._is_canonical() m1._is_canonical() m2._is_canonical() m3._is_canonical() + m4._is_canonical() + m5._is_canonical() f0 = q.sv.fidelity(q.mps.state_vector(m0), v) f1 = q.sv.fidelity(q.mps.state_vector(m1), v) f2 = q.sv.fidelity(q.mps.state_vector(m2), v) f3 = q.sv.fidelity(q.mps.state_vector(m3), v) - print("fidelity = {} {} {} {}".format(f0, f1, f2, f3)) + f4 = q.sv.fidelity(q.mps.state_vector(m4), v) + f5 = q.sv.fidelity(q.mps.state_vector(m5), v) + print("fidelity = {} {} {} {} {} {}".format(f0, f1, f2, f3, f4, f5)) assert f0 == approx(1) assert f1 == approx(1) + assert f2 == approx(1) f0 = q.sv.fidelity(q.mps.state_vector(m0), v) f1 = q.sv.fidelity(q.mps.state_vector(m1), v) f2 = q.sv.fidelity(q.mps.state_vector(m2), v) f3 = q.sv.fidelity(q.mps.state_vector(m3), v) - print("final fidelity = {} {} {} {}".format(f0, f1, f2, f3)) + f4 = q.sv.fidelity(q.mps.state_vector(m4), v) + f5 = q.sv.fidelity(q.mps.state_vector(m5), v) + print("final fidelity = {} {} {} {} {} {}".format(f0, f1, f2, f3, f4, f5)) assert f0 == approx(1) assert f1 == approx(1) + assert f2 == approx(1) # assert f2 == approx(f3) - m4 = q.mps.zero(n, MPS_P) - f4 = q.sv.fidelity(q.mps.state_vector(q.mps.apply_seq(m4, seq)), v) - assert f4 == approx(1) + m6 = q.mps.zero(n, q.mps_p.MPS_P) + f6 = q.sv.fidelity(q.mps.state_vector(q.mps.apply_seq(m6, seq)), v) + assert f6 == approx(1) + + m7 = q.mps.zero(n, q.mps_t.MPS_T) + f7 = q.sv.fidelity(q.mps.state_vector(q.mps.apply_seq(m7, seq)), v) + assert f7 == approx(1) if __name__ == "__main__": diff --git a/test/mps/test_canonical.py b/test/mps/test_canonical.py index 0176660..d8053d1 100644 --- a/test/mps/test_canonical.py +++ b/test/mps/test_canonical.py @@ -14,56 +14,43 @@ def test_canonical(): tensors.append(np.random.random((d, 2, dn))) d = dn tensors.append(np.random.random((d, 2, 1))) - m0 = q.mps.MPS_C(tensors) - m1 = q.mps_p.MPS_P(tensors) - norm = q.mps.norm(m0) - assert q.mps.norm(m0) == approx(norm) - assert q.mps.norm(m1) == approx(norm) - assert q.mps.is_canonical(m0) - assert q.mps.is_canonical(m1) - - for _ in range(n): - p = np.random.randint(n) - m0._canonicalize(p) - m1._canonicalize(p) - assert q.mps.norm(m0) == approx(norm) - assert q.mps.norm(m1) == approx(norm) - assert q.mps.is_canonical(m0) - assert q.mps.is_canonical(m1) - - for _ in range(n): - p = np.random.randint(n - 1) - m0._canonicalize(p, p + 1) - m1._canonicalize(p, p + 1) - assert q.mps.norm(m0) == approx(norm) - assert q.mps.norm(m1) == approx(norm) - assert q.mps.is_canonical(m0) - assert q.mps.is_canonical(m1) + for mps in [q.mps.MPS_C, q.mps_p.MPS_P, q.mps_t.MPS_T]: + m = mps(tensors) + norm = q.mps.norm(m) + assert q.mps.norm(m) == approx(norm) + assert q.mps.is_canonical(m) + + for _ in range(n): + p = np.random.randint(n) + m._canonicalize(p) + assert q.mps.norm(m) == approx(norm) + assert q.mps.is_canonical(m) + assert q.mps.is_canonical(m) + + for _ in range(n): + p = np.random.randint(n - 1) + m._canonicalize(p, p + 1) + assert q.mps.norm(m) == approx(norm) + assert q.mps.is_canonical(m) v = np.random.random(2**n).reshape((2,) * n + (1,)) v /= np.linalg.norm(v) tensors = q.mps.tensor_decomposition(v, maxdim) - m0 = q.mps.MPS_C(tensors) - m1 = q.mps_p.MPS_P(tensors) - norm = q.mps.norm(m0) - - for _ in range(n): - p = np.random.randint(n) - m0._canonicalize(p) - m1._canonicalize(p) - assert q.mps.norm(m0) == approx(norm) - assert q.mps.norm(m1) == approx(norm) - assert q.mps.is_canonical(m0) - assert q.mps.is_canonical(m1) - - for _ in range(n): - p = np.random.randint(n - 1) - m0._canonicalize(p, p + 1) - m1._canonicalize(p, p + 1) - assert q.mps.norm(m0) == approx(norm) - assert q.mps.norm(m1) == approx(norm) - assert q.mps.is_canonical(m0) - assert q.mps.is_canonical(m1) + for mps in [q.mps.MPS_C, q.mps_p.MPS_P, q.mps_t.MPS_T]: + m = mps(tensors) + norm = q.mps.norm(m) + + for _ in range(n): + p = np.random.randint(n) + m._canonicalize(p) + assert q.mps.norm(m) == approx(norm) + assert q.mps.is_canonical(m) + + for _ in range(n): + p = np.random.randint(n - 1) + m._canonicalize(p, p + 1) + assert q.mps.norm(m) == approx(norm) + assert q.mps.is_canonical(m) if __name__ == "__main__": diff --git a/test/mps/test_move.py b/test/mps/test_move.py index 1cc1c36..9d8a3bd 100644 --- a/test/mps/test_move.py +++ b/test/mps/test_move.py @@ -16,36 +16,23 @@ def test_swap(): tensors.append(np.random.random((d, 2, dn))) d = dn tensors.append(np.random.random((d, 2, 1))) - m0 = q.mps.MPS_C(tensors) - m1 = q.mps_p.MPS_P(tensors) - q.mps.is_canonical(m0) - q.mps.is_canonical(m1) - norm = q.mps.norm(m0) - v = q.sv.vector(q.mps.state_vector(m0)) + for mps in [q.mps.MPS_C, q.mps_p.MPS_P, q.mps_t.MPS_T]: + m = mps(tensors) + q.mps.is_canonical(m) + norm = q.mps.norm(m) + v = q.sv.vector(q.mps.state_vector(m)) + for _ in range(64): + s = np.random.randint(n - 1) + print(f"swap tensors at {s} and {s+1}") + m._canonicalize(s) + _swap_tensors(m, s) + print(q.sv.vector(q.mps.state_vector(m))) + q.mps.is_canonical(m) + assert q.mps.norm(m) == approx(norm) - for _ in range(64): - s = np.random.randint(n - 1) - print(f"swap tensors at {s} and {s+1}") - m0._canonicalize(s) - m1._canonicalize(s) - _swap_tensors(m0, s) - _swap_tensors(m1, s) - print(q.sv.vector(q.mps.state_vector(m0))) - print(q.sv.vector(q.mps.state_vector(m1))) - q.mps.is_canonical(m0) - q.mps.is_canonical(m1) - assert q.mps.norm(m0) == approx(norm) - assert q.mps.norm(m1) == approx(norm) - - v0 = q.sv.vector(q.mps.state_vector(m0)) - v1 = q.sv.vector(q.mps.state_vector(m1)) - assert len(v) == len(v0) - assert len(v) == len(v1) - print(v) - print(v0) - print(v1) - assert q.sv.is_close(v, v0) - assert q.sv.is_close(v, v1) + vn = q.sv.vector(q.mps.state_vector(m)) + assert len(v) == len(vn) + assert q.sv.is_close(v, vn) def test_move(): @@ -61,39 +48,25 @@ def test_move(): tensors.append(np.random.random((d, 2, dn))) d = dn tensors.append(np.random.random((d, 2, 1))) - m0 = q.mps.MPS_C(tensors) - m1 = q.mps_p.MPS_P(tensors) - q.mps.is_canonical(m0) - q.mps.is_canonical(m1) - norm = q.mps.norm(m0) - v = q.sv.vector(q.mps.state_vector(m0)) - print(q.sv.vector(q.mps.state_vector(m0))) - print(q.sv.vector(q.mps.state_vector(m1))) + for mps in [q.mps.MPS_C, q.mps_p.MPS_P, q.mps_t.MPS_T]: + m = mps(tensors) + q.mps.is_canonical(m) + norm = q.mps.norm(m) + v = q.sv.vector(q.mps.state_vector(m)) - for _ in range(16): - p = np.random.randint(n) - s = np.random.randint(n) - print(f"move qubit at {p} to {s}") - q.mps.is_canonical(m0) - q.mps.is_canonical(m1) - _move_qubit(m0, p, s) - _move_qubit(m1, p, s) - print(q.sv.vector(q.mps.state_vector(m0))) - print(q.sv.vector(q.mps.state_vector(m1))) - q.mps.is_canonical(m0) - q.mps.is_canonical(m1) - assert q.mps.norm(m0) == approx(norm) - assert q.mps.norm(m1) == approx(norm) + for _ in range(16): + p = np.random.randint(n) + s = np.random.randint(n) + print(f"move qubit at {p} to {s}") + q.mps.is_canonical(m) + _move_qubit(m, p, s) + print(q.sv.vector(q.mps.state_vector(m))) + q.mps.is_canonical(m) + assert q.mps.norm(m) == approx(norm) - v0 = q.sv.vector(q.mps.state_vector(m0)) - v1 = q.sv.vector(q.mps.state_vector(m1)) - assert len(v) == len(v0) - assert len(v) == len(v1) - print(v) - print(v0) - print(v1) - assert q.sv.is_close(v, v0) - assert q.sv.is_close(v, v1) + vn = q.sv.vector(q.mps.state_vector(m)) + assert len(v) == len(vn) + assert q.sv.is_close(v, vn) if __name__ == "__main__": diff --git a/test/mps/test_mps.py b/test/mps/test_mps.py index 0a4a53e..70888f6 100644 --- a/test/mps/test_mps.py +++ b/test/mps/test_mps.py @@ -6,8 +6,14 @@ def test_mps(): states = [q.sv.one(), q.sv.one(), q.sv.zero(), q.sv.zero()] v = q.sv.product_state(states) m0 = q.mps.product_state(states) + m1 = q.mps_p.product_state(states) + m2 = q.mps_t.product_state(states) v0 = q.mps.state_vector(m0) + v1 = q.mps.state_vector(m1) + v2 = q.mps.state_vector(m2) assert np.allclose(v0, v) + assert np.allclose(v1, v) + assert np.allclose(v2, v) if __name__ == "__main__": diff --git a/test/mps_p/test_mps_p.py b/test/mps_p/test_mps_p.py deleted file mode 100644 index 707cd2b..0000000 --- a/test/mps_p/test_mps_p.py +++ /dev/null @@ -1,15 +0,0 @@ -import numpy as np -import qailo as q -from qailo.mps_p.mps_p import MPS_P - - -def test_mps(): - states = [q.sv.one(), q.sv.one(), q.sv.zero(), q.sv.zero()] - v = q.sv.product_state(states) - m0 = q.mps.product_state(states, MPS_P) - v0 = q.mps.state_vector(m0) - assert np.allclose(v0, v) - - -if __name__ == "__main__": - test_mps() diff --git a/test/mps_t/test_mps_t.py b/test/mps_t/test_mps_t.py new file mode 100644 index 0000000..66a411d --- /dev/null +++ b/test/mps_t/test_mps_t.py @@ -0,0 +1,36 @@ +import qailo as q + + +def test_mps_t(): + v = q.mps_t.zero(3) + + print("input:") + print("state vector:", q.vector(v)) + print("probabitily:", q.probability(v)) + + v = q.apply(v, q.op.h(), [0]) + v = q.apply(v, q.op.h(), [2]) + v = q.apply(v, q.op.cx(), [0, 1]) + v = q.apply(v, q.op.cz(), [1, 2]) + v = q.apply(v, q.op.h(), [2]) + + print("output:") + print("state vector:", q.vector(v)) + print("probabitily:", q.probability(v)) + + print("# tensor pool") + for id, tp in enumerate(v.tpool): + print(f"{id} {tp[0].shape} {tp[1]} {tp[2]}") + assert len(v.tpool) == 25 + + print("# generator pool") + for id, gp in enumerate(v.gpool): + print(f"{id} {gp[0].shape} {gp[1]} {gp[2].shape} {gp[3].shape}") + assert len(v.gpool) == 2 + + prefix = "test_mps_t" + v._dump(prefix) + + +if __name__ == "__main__": + test_mps_t()