From bb7c7d7598224889cb8e7e129a69e12042a63ee1 Mon Sep 17 00:00:00 2001 From: Sora Shiratani Date: Sat, 23 Dec 2023 05:01:41 +0900 Subject: [PATCH 01/40] :label: Annotate bitops.py --- src/qailo/util/bitops.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/qailo/util/bitops.py b/src/qailo/util/bitops.py index 7a69aa8..7e58e46 100644 --- a/src/qailo/util/bitops.py +++ b/src/qailo/util/bitops.py @@ -1,15 +1,18 @@ -def bit(n, s, k): +from __future__ import annotations + + +def bit(n: int, s: int, k: int) -> int: """extract k-th bit from integer s""" return (s >> (n - k - 1)) & 1 -def binary2str(n, i): +def binary2str(n: int, i: int) -> str: return bin(i)[2:].zfill(n)[::-1] -def str2binary(s): +def str2binary(s: str) -> int: n = len(s) - c = 0 + c: int = 0 for i in range(n): assert s[i] == "0" or s[i] == "1" if s[i] == "1": From a6d6f28456f62b5f5ff0f5e7efca1ecfdcf1142b Mon Sep 17 00:00:00 2001 From: Sora Shiratani Date: Sat, 23 Dec 2023 04:51:18 +0900 Subject: [PATCH 02/40] :art: Remove len(arr.shape) --- src/qailo/mps/svd.py | 6 +++--- src/qailo/operator/type.py | 8 ++++---- src/qailo/state_vector/type.py | 2 +- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/qailo/mps/svd.py b/src/qailo/mps/svd.py index 934c7e2..deb701d 100644 --- a/src/qailo/mps/svd.py +++ b/src/qailo/mps/svd.py @@ -2,7 +2,7 @@ def compact_svd(A, nkeep=None, tol=1e-12): - assert len(A.shape) == 2 + assert A.ndim == 2 U, S, Vh = np.linalg.svd(A, full_matrices=False) V = Vh.conj().T dimS = sum([1 if x > tol * S[0] else 0 for x in S]) @@ -13,8 +13,8 @@ def compact_svd(A, nkeep=None, tol=1e-12): def tensor_svd(T, partition, canonical="center", nkeep=None, tol=1e-12): legsL = len(partition[0]) legsR = len(partition[1]) - assert len(T.shape) == legsL + legsR - assert sorted(partition[0] + partition[1]) == list(range(len(T.shape))) + assert T.ndim == legsL + legsR + assert sorted(partition[0] + partition[1]) == list(range(T.ndim)) dimsL = [T.shape[i] for i in partition[0]] dimsR = [T.shape[i] for i in partition[1]] m = np.einsum(T, partition[0] + partition[1]).reshape( diff --git a/src/qailo/operator/type.py b/src/qailo/operator/type.py index 4c35c3f..261ccff 100644 --- a/src/qailo/operator/type.py +++ b/src/qailo/operator/type.py @@ -7,20 +7,20 @@ def is_density_matrix(v): v.shape[-1] == 1 and v.shape[-2] == 1 and v.shape[-3] > 1 - and len(v.shape) % 2 == 0 + and v.ndim % 2 == 0 ) return False def is_operator(v): if isinstance(v, np.ndarray): - return v.shape[-1] > 1 and len(v.shape) % 2 == 0 + return v.shape[-1] > 1 and v.ndim % 2 == 0 return False def num_qubits(v): if is_density_matrix(v): - return (len(v.shape) - 2) // 2 + return (v.ndim - 2) // 2 elif is_operator(v): - return len(v.shape) // 2 + return v.ndim // 2 raise ValueError diff --git a/src/qailo/state_vector/type.py b/src/qailo/state_vector/type.py index 49fe5fe..bca64e3 100644 --- a/src/qailo/state_vector/type.py +++ b/src/qailo/state_vector/type.py @@ -9,4 +9,4 @@ def is_state_vector(v): def num_qubits(v): - return len(v.shape) - 1 + return v.ndim - 1 From 2298e62ef139fd753565ffd8b14eb5092ac0b163 Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 02:29:07 +0900 Subject: [PATCH 03/40] :technologist: Add typeutil --- src/qailo/typeutil/__init__.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 src/qailo/typeutil/__init__.py diff --git a/src/qailo/typeutil/__init__.py b/src/qailo/typeutil/__init__.py new file mode 100644 index 0000000..e69de29 From bea91a1c15bb48a10f7ba08a3e2d6f4e01f2c272 Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 02:31:56 +0900 Subject: [PATCH 04/40] :construction: Annotate state_vector/type.py --- src/qailo/state_vector/type.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/qailo/state_vector/type.py b/src/qailo/state_vector/type.py index bca64e3..19aaf8b 100644 --- a/src/qailo/state_vector/type.py +++ b/src/qailo/state_vector/type.py @@ -1,12 +1,18 @@ +from __future__ import annotations + +from typing import Any + import numpy as np +import numpy.typing as npt +from typing_extensions import TypeGuard -def is_state_vector(v): +def is_state_vector(v: Any) -> TypeGuard[np.ndarray]: if isinstance(v, np.ndarray): return v.shape[-1] == 1 and v.shape[-2] > 1 else: return False -def num_qubits(v): +def num_qubits(v: npt.NDArray) -> int: return v.ndim - 1 From aba083887ee84f41dc46e188998d7d16de5dd8e0 Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 02:34:23 +0900 Subject: [PATCH 05/40] :construction: Annotate operator/type.py --- src/qailo/operator/type.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/src/qailo/operator/type.py b/src/qailo/operator/type.py index 261ccff..1984bda 100644 --- a/src/qailo/operator/type.py +++ b/src/qailo/operator/type.py @@ -1,7 +1,13 @@ +from __future__ import annotations + +from typing import Any + import numpy as np +import numpy.typing as npt +from typing_extensions import TypeGuard -def is_density_matrix(v): +def is_density_matrix(v: Any) -> TypeGuard[np.ndarray]: if isinstance(v, np.ndarray): return ( v.shape[-1] == 1 @@ -12,13 +18,13 @@ def is_density_matrix(v): return False -def is_operator(v): +def is_operator(v: Any) -> TypeGuard[np.ndarray]: if isinstance(v, np.ndarray): return v.shape[-1] > 1 and v.ndim % 2 == 0 return False -def num_qubits(v): +def num_qubits(v: npt.NDArray) -> int: if is_density_matrix(v): return (v.ndim - 2) // 2 elif is_operator(v): From bd679cdf52d8ff25624aed081a68513a14ee505d Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 03:11:39 +0900 Subject: [PATCH 06/40] :safety_vest: Add type validation --- src/qailo/state_vector/fidelity.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/qailo/state_vector/fidelity.py b/src/qailo/state_vector/fidelity.py index 7fe2bcb..fb47d26 100644 --- a/src/qailo/state_vector/fidelity.py +++ b/src/qailo/state_vector/fidelity.py @@ -4,4 +4,6 @@ def fidelity(v0, v1): v0 = v0 / np.linalg.norm(v0) v1 = v1 / np.linalg.norm(v1) - return np.abs(np.vdot(v0, v1)) ** 2 + ret = np.abs(np.vdot(v0, v1)) ** 2 + assert isinstance(ret, float) + return ret From a24b604e23c70a7379bb4b3dc0baebe16e36fb0b Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 04:20:27 +0900 Subject: [PATCH 07/40] :safety_vest: Add eincheck.py --- src/qailo/typeutil/eincheck.py | 70 ++++++++++++++++++++++++++++++++++ 1 file changed, 70 insertions(+) create mode 100644 src/qailo/typeutil/eincheck.py diff --git a/src/qailo/typeutil/eincheck.py b/src/qailo/typeutil/eincheck.py new file mode 100644 index 0000000..b172926 --- /dev/null +++ b/src/qailo/typeutil/eincheck.py @@ -0,0 +1,70 @@ +"""Wrap np.einsum for better type checking.""" + +from __future__ import annotations + +from typing import overload + +import numpy as np +import numpy.typing as npt + + +@overload +def einsum(op1: np.ndarray, ss1: list[int]) -> npt.ArrayLike: ... + + +@overload +def einsum(op1: np.ndarray, ss1: list[int], ss_out: list[int]) -> npt.ArrayLike: ... + + +@overload +def einsum( + op1: np.ndarray, + ss1: list[int], + op2: np.ndarray, + ss2: list[int], +) -> npt.ArrayLike: ... + + +@overload +def einsum( + op1: np.ndarray, + ss1: list[int], + op2: np.ndarray, + ss2: list[int], + ss_out: list[int], +) -> npt.ArrayLike: ... + + +def einsum(*args, **kwargs): + return np.einsum(*args, **kwargs) + + +@overload +def einsum_cast(op1: np.ndarray, ss1: list[int]) -> npt.NDArray: ... + + +@overload +def einsum_cast(op1: np.ndarray, ss1: list[int], ss_out: list[int]) -> npt.NDArray: ... + + +@overload +def einsum_cast( + op1: np.ndarray, + ss1: list[int], + op2: np.ndarray, + ss2: list[int], +) -> npt.NDArray: ... + + +@overload +def einsum_cast( + op1: np.ndarray, + ss1: list[int], + op2: np.ndarray, + ss2: list[int], + ss_out: list[int], +) -> npt.NDArray: ... + + +def einsum_cast(*args, **kwargs): + return np.asarray(np.einsum(*args, **kwargs)) From 20e1d704e279c5f50a98e66d40befafacb44d85f Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 03:48:30 +0900 Subject: [PATCH 08/40] :construction: Annotate trace.py As a dependency of state_vector --- src/qailo/operator/trace.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/qailo/operator/trace.py b/src/qailo/operator/trace.py index 39dffd4..1731811 100644 --- a/src/qailo/operator/trace.py +++ b/src/qailo/operator/trace.py @@ -1,9 +1,14 @@ -import numpy as np +from __future__ import annotations +from typing import Collection + +import numpy.typing as npt + +from ..typeutil import eincheck as ec from .type import is_density_matrix, is_operator, num_qubits -def trace(q, pos=None): +def trace(q: npt.NDArray, pos: Collection[int] | None = None) -> npt.ArrayLike: assert is_density_matrix(q) or is_operator(q) n = num_qubits(q) if pos is None: @@ -16,4 +21,4 @@ def trace(q, pos=None): ss = list(range(2 * n)) for i in pos: ss[n + i] = ss[i] - return np.einsum(q.reshape((2,) * (2 * n)), ss) + return ec.einsum(q.reshape((2,) * (2 * n)), ss) From 4f3d619fa8259d45bd4fe3b6ff575b20ce1da256 Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 03:42:05 +0900 Subject: [PATCH 09/40] :art: Use np.real to use implicit cast --- src/qailo/state_vector/probability.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/qailo/state_vector/probability.py b/src/qailo/state_vector/probability.py index df6852a..23568df 100644 --- a/src/qailo/state_vector/probability.py +++ b/src/qailo/state_vector/probability.py @@ -17,4 +17,5 @@ def probability(v, pos=None): for k in range(num_qubits(w)): if k not in pos: tpos.append(k) - return np.diag(matrix(trace(density_matrix(w), tpos).real)) + traced = np.real(trace(density_matrix(w), tpos)) + return np.diag(matrix(traced)) From 8a2680f81272b7d47374f2a51978415fe368fee1 Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 03:44:59 +0900 Subject: [PATCH 10/40] :rotating_light: Cast to fix type --- src/qailo/state_vector/probability.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/qailo/state_vector/probability.py b/src/qailo/state_vector/probability.py index 23568df..4b64a44 100644 --- a/src/qailo/state_vector/probability.py +++ b/src/qailo/state_vector/probability.py @@ -9,7 +9,7 @@ def probability(v, pos=None): assert is_state_vector(v) - w = v / np.linalg.norm(v) + w = v / float(np.linalg.norm(v)) if pos is None: return abs(vector(w)) ** 2 else: From d71f1517e2187577a8726e83951dd51eb9225ef8 Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 03:46:18 +0900 Subject: [PATCH 11/40] :construction: Annotate state_vector module --- src/qailo/state_vector/apply.py | 20 ++++++++++++++++---- src/qailo/state_vector/density_matrix.py | 8 ++++++-- src/qailo/state_vector/fidelity.py | 5 ++++- src/qailo/state_vector/is_close.py | 5 ++++- src/qailo/state_vector/probability.py | 7 ++++++- src/qailo/state_vector/state_vector.py | 14 ++++++++++---- src/qailo/state_vector/vector.py | 6 +++++- 7 files changed, 51 insertions(+), 14 deletions(-) diff --git a/src/qailo/state_vector/apply.py b/src/qailo/state_vector/apply.py index 8899978..e3d4d6a 100644 --- a/src/qailo/state_vector/apply.py +++ b/src/qailo/state_vector/apply.py @@ -1,10 +1,19 @@ -import numpy as np +from __future__ import annotations + +from typing import Iterable, Sequence + +import numpy.typing as npt from ..operator import type as op +from ..typeutil import eincheck as ec from . import type as sv -def apply(v, p, pos=None): +def apply( + v: npt.NDArray, + p: npt.NDArray, + pos: Sequence[int] | None = None, +) -> npt.NDArray: assert sv.is_state_vector(v) and op.is_operator(p) n = sv.num_qubits(v) m = op.num_qubits(p) @@ -19,10 +28,13 @@ def apply(v, p, pos=None): for i in range(m): ss_v[pos[i]] = ss_op[m + i] ss_to[pos[i]] = ss_op[i] - return np.einsum(v, ss_v, p, ss_op, ss_to) + return ec.einsum_cast(v, ss_v, p, ss_op, ss_to) -def apply_seq(v, seq): +def apply_seq( + v: npt.NDArray, + seq: Iterable[tuple[npt.NDArray, Sequence[int]]], +) -> npt.NDArray: for p, qubit in seq: v = apply(v, p, qubit) return v diff --git a/src/qailo/state_vector/density_matrix.py b/src/qailo/state_vector/density_matrix.py index c0a6753..34176a3 100644 --- a/src/qailo/state_vector/density_matrix.py +++ b/src/qailo/state_vector/density_matrix.py @@ -1,9 +1,13 @@ +from __future__ import annotations + import numpy as np +import numpy.typing as npt +from ..typeutil import eincheck as ec from .type import is_state_vector, num_qubits -def density_matrix(v): +def density_matrix(v: npt.NDArray) -> npt.NDArray: assert is_state_vector(v) n = num_qubits(v) w = v.copy() @@ -11,4 +15,4 @@ def density_matrix(v): ss0 = list(range(n)) ss1 = list(range(n, 2 * n)) shape = (2,) * (2 * n) + (1, 1) - return np.einsum(w, ss0, w.conj(), ss1).reshape(shape) + return ec.einsum_cast(w, ss0, w.conj(), ss1).reshape(shape) diff --git a/src/qailo/state_vector/fidelity.py b/src/qailo/state_vector/fidelity.py index fb47d26..e3f2664 100644 --- a/src/qailo/state_vector/fidelity.py +++ b/src/qailo/state_vector/fidelity.py @@ -1,7 +1,10 @@ +from __future__ import annotations + import numpy as np +import numpy.typing as npt -def fidelity(v0, v1): +def fidelity(v0: npt.NDArray, v1: npt.NDArray) -> float: v0 = v0 / np.linalg.norm(v0) v1 = v1 / np.linalg.norm(v1) ret = np.abs(np.vdot(v0, v1)) ** 2 diff --git a/src/qailo/state_vector/is_close.py b/src/qailo/state_vector/is_close.py index 2c6c3a7..ef6f5cb 100644 --- a/src/qailo/state_vector/is_close.py +++ b/src/qailo/state_vector/is_close.py @@ -1,5 +1,8 @@ +from __future__ import annotations + import numpy as np +import numpy.typing as npt -def is_close(p0, p1): +def is_close(p0: npt.NDArray, p1: npt.NDArray) -> bool: return p0.shape == p1.shape and np.allclose(p0, p1) diff --git a/src/qailo/state_vector/probability.py b/src/qailo/state_vector/probability.py index 4b64a44..da63b5c 100644 --- a/src/qailo/state_vector/probability.py +++ b/src/qailo/state_vector/probability.py @@ -1,4 +1,9 @@ +from __future__ import annotations + +from typing import Container + import numpy as np +import numpy.typing as npt from ..operator.matrix import matrix from ..operator.trace import trace @@ -7,7 +12,7 @@ from .vector import vector -def probability(v, pos=None): +def probability(v: npt.NDArray, pos: Container[int] | None = None) -> npt.NDArray: assert is_state_vector(v) w = v / float(np.linalg.norm(v)) if pos is None: diff --git a/src/qailo/state_vector/state_vector.py b/src/qailo/state_vector/state_vector.py index 4b9f2e2..f976a05 100644 --- a/src/qailo/state_vector/state_vector.py +++ b/src/qailo/state_vector/state_vector.py @@ -1,9 +1,15 @@ +from __future__ import annotations + +from typing import Sequence + import numpy as np +import numpy.typing as npt +from ..typeutil import eincheck as ec from .type import is_state_vector, num_qubits -def product_state(states): +def product_state(states: Sequence[npt.NDArray]) -> npt.NDArray: m = len(states) assert m > 0 v = states[0] @@ -11,7 +17,7 @@ def product_state(states): for i in range(1, m): n0 = num_qubits(v) n1 = num_qubits(states[i]) - v = np.einsum( + v = ec.einsum_cast( v.reshape((2,) * n0), list(range(n0)), states[i], @@ -20,7 +26,7 @@ def product_state(states): return v -def zero(n=1): +def zero(n: int = 1) -> npt.NDArray: assert n > 0 if n == 1: return np.array((1.0, 0.0)).reshape((2, 1)) @@ -28,7 +34,7 @@ def zero(n=1): return product_state([zero(1)] * n) -def one(n=1): +def one(n: int = 1) -> npt.NDArray: assert n > 0 if n == 1: return np.array((0.0, 1.0)).reshape((2, 1)) diff --git a/src/qailo/state_vector/vector.py b/src/qailo/state_vector/vector.py index 09af217..1d6af1e 100644 --- a/src/qailo/state_vector/vector.py +++ b/src/qailo/state_vector/vector.py @@ -1,7 +1,11 @@ +from __future__ import annotations + +import numpy.typing as npt + from . import type as sv -def vector(v, c=None): +def vector(v: npt.NDArray, c: list | None = None) -> npt.NDArray: assert sv.is_state_vector(v) if c is None: n = sv.num_qubits(v) From 119172f2f73ad286badfa5addb5ddb7add83ed4e Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 04:13:56 +0900 Subject: [PATCH 12/40] :art: Add OPAutomaton --- src/qailo/operator/type.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/qailo/operator/type.py b/src/qailo/operator/type.py index 1984bda..36aee5a 100644 --- a/src/qailo/operator/type.py +++ b/src/qailo/operator/type.py @@ -1,6 +1,6 @@ from __future__ import annotations -from typing import Any +from typing import Any, NamedTuple import numpy as np import numpy.typing as npt @@ -30,3 +30,8 @@ def num_qubits(v: npt.NDArray) -> int: elif is_operator(v): return v.ndim // 2 raise ValueError + + +class OPAutomaton(NamedTuple): + op: npt.NDArray + pos: list[int] From 65ee3c4abba4abd25ebf5a54b96aac0597ec21d1 Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 04:17:54 +0900 Subject: [PATCH 13/40] :art: Remove inverse slice --- src/qailo/operator/inverse.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/qailo/operator/inverse.py b/src/qailo/operator/inverse.py index 2001619..34cedf1 100644 --- a/src/qailo/operator/inverse.py +++ b/src/qailo/operator/inverse.py @@ -4,7 +4,7 @@ def inverse_seq(seq): res = [] - for p, pos in seq[::-1]: + for p, pos in reversed(seq): assert is_operator(p) res.append([hconj(p), pos]) return res From ebc547319ceda4586a4f5fe33181fea88854b2db Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 04:27:51 +0900 Subject: [PATCH 14/40] :boom: Use OPAutomaton instead of list Co-authored-by: Synge Todo --- src/qailo/operator/controlled.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/qailo/operator/controlled.py b/src/qailo/operator/controlled.py index 1bf9b88..8cf5815 100644 --- a/src/qailo/operator/controlled.py +++ b/src/qailo/operator/controlled.py @@ -73,10 +73,10 @@ def controlled_seq(u, pos): if n <= 1: raise ValueError elif n == 2: - seq.append([controlled(u), pos]) + seq.append(OPAutomaton(controlled(u), pos)) else: - seq.append([control_begin(), [pos[0], pos[1]]]) + seq.append(OPAutomaton(control_begin(), [pos[0], pos[1]])) for i in range(1, n - 2): - seq.append([control_propagate(), [pos[i], pos[i + 1]]]) - seq.append([control_end(u), [pos[-2], pos[-1]]]) + seq.append(OPAutomaton(control_propagate(), [pos[i], pos[i + 1]])) + seq.append(OPAutomaton(control_end(u), [pos[-2], pos[-1]])) return seq From e15802bb275eca95067d648c66785caa0c54ddae Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 04:29:36 +0900 Subject: [PATCH 15/40] :construction: Annotate operator module --- src/qailo/operator/controlled.py | 23 +++++++++++++---------- src/qailo/operator/hconj.py | 9 ++++++--- src/qailo/operator/identity.py | 5 ++++- src/qailo/operator/inverse.py | 12 ++++++++---- src/qailo/operator/matrix.py | 13 ++++++++++++- src/qailo/operator/multiply.py | 13 ++++++++++--- src/qailo/operator/one_qubit.py | 23 +++++++++++++---------- src/qailo/operator/property.py | 13 ++++++++----- src/qailo/operator/swap.py | 5 ++++- 9 files changed, 78 insertions(+), 38 deletions(-) diff --git a/src/qailo/operator/controlled.py b/src/qailo/operator/controlled.py index 8cf5815..23d8e15 100644 --- a/src/qailo/operator/controlled.py +++ b/src/qailo/operator/controlled.py @@ -1,11 +1,14 @@ +from __future__ import annotations + import numpy as np +import numpy.typing as npt from .matrix import matrix from .one_qubit import p, x, z -from .type import is_operator, num_qubits +from .type import OPAutomaton, is_operator, num_qubits -def controlled(u): +def controlled(u: npt.NDArray) -> npt.NDArray: assert is_operator(u) m = num_qubits(u) n = m + 1 @@ -14,18 +17,18 @@ def controlled(u): return op.reshape((2,) * (2 * n)) -def cp(phi): +def cp(phi: float) -> npt.NDArray: return controlled(p(phi)) -def cx(n=2): +def cx(n: int = 2) -> npt.NDArray: op = x() for _ in range(n - 1): op = controlled(op) return op -def cz(n=2): +def cz(n: int = 2) -> npt.NDArray: op = z() for _ in range(n - 1): op = controlled(op) @@ -33,7 +36,7 @@ def cz(n=2): # mpo representation of controlled gates -def control_begin(): +def control_begin() -> npt.NDArray: op = np.zeros([2, 2, 2, 2, 2]) op[0, 0, 0, 0, 0] = 1 op[0, 1, 0, 0, 1] = 1 @@ -42,7 +45,7 @@ def control_begin(): return op.reshape([2, 4, 2, 2]) -def control_propagate(): +def control_propagate() -> npt.NDArray: op = np.zeros([2, 2, 2, 2, 2, 2]) op[0, 0, 0, 0, 0, 0] = 1 op[1, 0, 0, 1, 0, 0] = 1 @@ -55,7 +58,7 @@ def control_propagate(): return op.reshape([2, 4, 4, 2]) -def control_end(u): +def control_end(u: npt.NDArray) -> npt.NDArray: assert is_operator(u) m = num_qubits(u) n = m + 1 @@ -67,9 +70,9 @@ def control_end(u): return op.reshape((2,) * n + (4,) + (2,) * m) -def controlled_seq(u, pos): +def controlled_seq(u: npt.NDArray, pos: list[int]) -> list[OPAutomaton]: n = len(pos) - seq = [] + seq: list[OPAutomaton] = [] if n <= 1: raise ValueError elif n == 2: diff --git a/src/qailo/operator/hconj.py b/src/qailo/operator/hconj.py index 5b5f7b3..c5a554c 100644 --- a/src/qailo/operator/hconj.py +++ b/src/qailo/operator/hconj.py @@ -1,11 +1,14 @@ -import numpy as np +from __future__ import annotations +import numpy.typing as npt + +from ..typeutil import eincheck as ec from .type import is_operator, num_qubits -def hconj(op): +def hconj(op: npt.NDArray) -> npt.NDArray: assert is_operator(op) n = num_qubits(op) ss_from = list(range(2 * n)) ss_to = list(range(n, 2 * n)) + list(range(n)) - return np.einsum(op, ss_from, ss_to).conjugate() + return ec.einsum_cast(op, ss_from, ss_to).conjugate() diff --git a/src/qailo/operator/identity.py b/src/qailo/operator/identity.py index 1d7c13c..a5ba1d8 100644 --- a/src/qailo/operator/identity.py +++ b/src/qailo/operator/identity.py @@ -1,5 +1,8 @@ +from __future__ import annotations + import numpy as np +import numpy.typing as npt -def identity(n): +def identity(n: int) -> npt.NDArray: return np.identity(2**n).reshape((2,) * (2 * n)) diff --git a/src/qailo/operator/inverse.py b/src/qailo/operator/inverse.py index 34cedf1..a325efe 100644 --- a/src/qailo/operator/inverse.py +++ b/src/qailo/operator/inverse.py @@ -1,10 +1,14 @@ +from __future__ import annotations + +from typing import Reversible + from .hconj import hconj -from .type import is_operator +from .type import OPAutomaton, is_operator -def inverse_seq(seq): - res = [] +def inverse_seq(seq: Reversible[OPAutomaton]) -> list[OPAutomaton]: + res: list[OPAutomaton] = [] for p, pos in reversed(seq): assert is_operator(p) - res.append([hconj(p), pos]) + res.append(OPAutomaton(hconj(p), pos)) return res diff --git a/src/qailo/operator/matrix.py b/src/qailo/operator/matrix.py index a22d737..58fc40e 100644 --- a/src/qailo/operator/matrix.py +++ b/src/qailo/operator/matrix.py @@ -1,7 +1,18 @@ +from __future__ import annotations + +import typing +from typing import TypeVar + +import numpy as np +import numpy.typing as npt + from .type import is_density_matrix, is_operator, num_qubits +T_co = TypeVar("T_co", covariant=True, bound=np.generic) + -def matrix(op): +def matrix(op: npt.NDArray[T_co]) -> npt.NDArray[T_co]: assert is_density_matrix(op) or is_operator(op) + op = typing.cast(npt.NDArray[T_co], op) n = num_qubits(op) return op.reshape([2**n, 2**n]) diff --git a/src/qailo/operator/multiply.py b/src/qailo/operator/multiply.py index 915369b..a3feccb 100644 --- a/src/qailo/operator/multiply.py +++ b/src/qailo/operator/multiply.py @@ -1,9 +1,16 @@ -import numpy as np +from __future__ import annotations +from typing import Sequence + +import numpy.typing as npt + +from ..typeutil import eincheck as ec from .type import is_operator, num_qubits -def multiply(pin, p, pos=None): +def multiply( + pin: npt.NDArray, p: npt.NDArray, pos: Sequence[int] | None = None +) -> npt.NDArray: assert is_operator(pin) assert is_operator(p) n = num_qubits(pin) @@ -21,4 +28,4 @@ def multiply(pin, p, pos=None): for i in range(m): ss_pin[pos[i]] = ss_p[m + i] ss_to[pos[i]] = ss_p[i] - return np.einsum(pin, ss_pin, p, ss_p, ss_to) + return ec.einsum_cast(pin, ss_pin, p, ss_p, ss_to) diff --git a/src/qailo/operator/one_qubit.py b/src/qailo/operator/one_qubit.py index e09bfee..25149d4 100644 --- a/src/qailo/operator/one_qubit.py +++ b/src/qailo/operator/one_qubit.py @@ -1,45 +1,48 @@ +from __future__ import annotations + import numpy as np +import numpy.typing as npt -def h(): +def h() -> npt.NDArray: return np.array([[1.0, 1.0], [1.0, -1.0]]) / np.sqrt(2) -def p(phi): +def p(phi: float) -> npt.NDArray: return np.array([[1.0, 0.0], [0.0, np.exp(1.0j * phi)]]) -def rx(phi): +def rx(phi: float) -> npt.NDArray: c = np.cos(phi / 2) s = np.sin(phi / 2) return np.array([[c, -1.0j * s], [-1.0j * s, c]]) -def ry(phi): +def ry(phi: float) -> npt.NDArray: c = np.cos(phi / 2) s = np.sin(phi / 2) return np.array([[c, -s], [s, c]]) -def rz(phi): +def rz(phi: float) -> npt.NDArray: return np.array([[np.exp(-1.0j * phi / 2), 0], [0.0, np.exp(1.0j * phi / 2)]]) -def s(): +def s() -> npt.NDArray: return np.array([[1.0, 0.0], [0.0, 1.0j]]) -def t(): +def t() -> npt.NDArray: return np.array([[1.0, 0.0], [0.0, np.exp(1.0j * np.pi / 4)]]) -def x(): +def x() -> npt.NDArray: return np.array([[0.0, 1.0], [1.0, 0.0]]) -def y(): +def y() -> npt.NDArray: return np.array([[0.0, -1.0j], [1.0j, 0.0]]) -def z(): +def z() -> npt.NDArray: return np.array([[1.0, 0.0], [0.0, -1.0]]) diff --git a/src/qailo/operator/property.py b/src/qailo/operator/property.py index 2b71100..47891d7 100644 --- a/src/qailo/operator/property.py +++ b/src/qailo/operator/property.py @@ -1,4 +1,7 @@ +from __future__ import annotations + import numpy as np +import numpy.typing as npt from .hconj import hconj from .identity import identity @@ -7,21 +10,21 @@ from .type import is_operator, num_qubits -def is_close(p0, p1): +def is_close(p0: npt.NDArray, p1: npt.NDArray) -> bool: return p0.shape == p1.shape and np.allclose(p0, p1) -def is_hermitian(op): +def is_hermitian(op: npt.NDArray) -> bool: return np.allclose(op, hconj(op)) -def is_identity(op): +def is_identity(op: npt.NDArray) -> bool: assert is_operator(op) n = num_qubits(op) return np.allclose(op, identity(n)) -def is_semi_positive(op): +def is_semi_positive(op: npt.NDArray) -> bool: if not is_hermitian(op): return False evs = np.linalg.eigvalsh(matrix(op)) @@ -31,7 +34,7 @@ def is_semi_positive(op): return True -def is_unitary(op): +def is_unitary(op: npt.NDArray) -> bool: assert is_operator(op) n = num_qubits(op) return is_identity(multiply(hconj(op), op, range(n))) and is_identity( diff --git a/src/qailo/operator/swap.py b/src/qailo/operator/swap.py index 6be4eb2..37bad4b 100644 --- a/src/qailo/operator/swap.py +++ b/src/qailo/operator/swap.py @@ -1,6 +1,9 @@ +from __future__ import annotations + import numpy as np +import numpy.typing as npt -def swap(): +def swap() -> npt.NDArray: op = np.array([[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]]) return op.reshape([2, 2, 2, 2]) From 619f24a8b1cded7736c6c653971d94cda89a69c1 Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 04:37:46 +0900 Subject: [PATCH 16/40] :safety_vest: Add string overload --- src/qailo/typeutil/eincheck.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/qailo/typeutil/eincheck.py b/src/qailo/typeutil/eincheck.py index b172926..2c326b1 100644 --- a/src/qailo/typeutil/eincheck.py +++ b/src/qailo/typeutil/eincheck.py @@ -8,6 +8,10 @@ import numpy.typing as npt +@overload +def einsum(ss: str, *op: npt.NDArray) -> npt.ArrayLike: ... + + @overload def einsum(op1: np.ndarray, ss1: list[int]) -> npt.ArrayLike: ... @@ -39,6 +43,10 @@ def einsum(*args, **kwargs): return np.einsum(*args, **kwargs) +@overload +def einsum_cast(ss: str, *op: npt.NDArray) -> npt.NDArray: ... + + @overload def einsum_cast(op1: np.ndarray, ss1: list[int]) -> npt.NDArray: ... From 94d0e28956871c5f54a9060046b5f96268ff6427 Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 04:54:54 +0900 Subject: [PATCH 17/40] :safety_vest: Update mps as ABC --- src/qailo/mps/type.py | 38 +++++++++++++++++++++++++++++++++++--- 1 file changed, 35 insertions(+), 3 deletions(-) diff --git a/src/qailo/mps/type.py b/src/qailo/mps/type.py index dd433a7..a182452 100644 --- a/src/qailo/mps/type.py +++ b/src/qailo/mps/type.py @@ -1,6 +1,38 @@ -class mps(object): - def __init__(self): - True +from __future__ import annotations + +from abc import ABC, abstractmethod + +import numpy.typing as npt +from typing_extensions import SupportsIndex + + +class mps(ABC): + @abstractmethod + def __init__(self, tensors: list[npt.NDArray], nkeep: int | None) -> None: ... + + @abstractmethod + def _num_qubits(self) -> int: ... + + @abstractmethod + def _norm(self) -> float: ... + + @abstractmethod + def _state_vector(self) -> npt.NDArray: ... + + @abstractmethod + def _tensor(self, t: SupportsIndex | slice) -> npt.NDArray: ... + + @abstractmethod + def _canonicalize(self, p0: int, p1: int | None) -> None: ... + + @abstractmethod + def _is_canonical(self) -> bool: ... + + @abstractmethod + def _apply_one(self, p: npt.NDArray, s: int) -> None: ... + + @abstractmethod + def _apply_two(self, p: npt.NDArray, s: int, reverse: bool | None) -> None: ... def is_canonical(m): From 29c6fa0ccb66fb453840f8dfa09b4197fec5cba8 Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 05:09:53 +0900 Subject: [PATCH 18/40] :bug: Change index type --- src/qailo/mps/type.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/qailo/mps/type.py b/src/qailo/mps/type.py index a182452..e26cb85 100644 --- a/src/qailo/mps/type.py +++ b/src/qailo/mps/type.py @@ -3,7 +3,7 @@ from abc import ABC, abstractmethod import numpy.typing as npt -from typing_extensions import SupportsIndex +from typing_extensions import TypeGuard class mps(ABC): @@ -20,7 +20,7 @@ def _norm(self) -> float: ... def _state_vector(self) -> npt.NDArray: ... @abstractmethod - def _tensor(self, t: SupportsIndex | slice) -> npt.NDArray: ... + def _tensor(self, t: int) -> npt.NDArray: ... @abstractmethod def _canonicalize(self, p0: int, p1: int | None) -> None: ... From e188a1c6d433f6ad474ca99c02cd6e497ae45b5e Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 05:15:09 +0900 Subject: [PATCH 19/40] :safety_vest: Add typecheck --- src/qailo/mps/mps_c.py | 4 +++- src/qailo/mps/mps_p.py | 6 +++--- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/qailo/mps/mps_c.py b/src/qailo/mps/mps_c.py index 6ddf013..27abc7a 100644 --- a/src/qailo/mps/mps_c.py +++ b/src/qailo/mps/mps_c.py @@ -39,7 +39,9 @@ def _norm(self): for t in range(self._num_qubits()): A = np.einsum("ij,jkl->ikl", A, self._tensor(t)) A = np.einsum("ijk,ijl->kl", A, self._tensor(t).conj()) - return np.sqrt(np.trace(A)) + ret = np.sqrt(np.trace(A)) + assert isinstance(ret, float) + return ret def _state_vector(self): n = self._num_qubits() diff --git a/src/qailo/mps/mps_p.py b/src/qailo/mps/mps_p.py index bb09661..40d6c7b 100644 --- a/src/qailo/mps/mps_p.py +++ b/src/qailo/mps/mps_p.py @@ -41,9 +41,9 @@ def _num_qubits(self): def _norm(self): A = np.identity(1) for t in range(self._num_qubits()): - A = np.einsum("ij,jkl->ikl", A, self._tensor(t)) - A = np.einsum("ijk,ijl->kl", A, self._tensor(t).conj()) - return np.sqrt(np.trace(A)) + ret = np.sqrt(np.trace(A)) + assert isinstance(ret, float) + return ret def _state_vector(self): n = self._num_qubits() From ef48b6ff42ccf3495aaa9805d952d8df198c6e27 Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 05:32:00 +0900 Subject: [PATCH 20/40] :bug: Add default value --- src/qailo/mps/type.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/qailo/mps/type.py b/src/qailo/mps/type.py index e26cb85..1834a54 100644 --- a/src/qailo/mps/type.py +++ b/src/qailo/mps/type.py @@ -32,7 +32,7 @@ def _is_canonical(self) -> bool: ... def _apply_one(self, p: npt.NDArray, s: int) -> None: ... @abstractmethod - def _apply_two(self, p: npt.NDArray, s: int, reverse: bool | None) -> None: ... + def _apply_two(self, p: npt.NDArray, s: int, reverse: bool = False) -> None: ... def is_canonical(m): From 1c3462b95e67d263b53bff857f7b7b3f6f82086b Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 05:41:09 +0900 Subject: [PATCH 21/40] :bug: Expose members --- src/qailo/mps/type.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/qailo/mps/type.py b/src/qailo/mps/type.py index 1834a54..6a71d3b 100644 --- a/src/qailo/mps/type.py +++ b/src/qailo/mps/type.py @@ -7,6 +7,9 @@ class mps(ABC): + q2t: list[int] + t2q: list[int] + @abstractmethod def __init__(self, tensors: list[npt.NDArray], nkeep: int | None) -> None: ... From 8e7d8b1eb6d2572697c453085aacf44ba576ee60 Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 05:47:50 +0900 Subject: [PATCH 22/40] :construction: Annotate mps module --- src/qailo/mps/apply.py | 21 ++++++--- src/qailo/mps/mps_c.py | 58 ++++++++++++++--------- src/qailo/mps/mps_p.py | 85 +++++++++++++++++++++------------- src/qailo/mps/product_state.py | 32 +++++++++---- src/qailo/mps/projector.py | 41 ++++++++++------ src/qailo/mps/state_vector.py | 8 +++- src/qailo/mps/svd.py | 33 ++++++++++--- src/qailo/mps/type.py | 7 +-- 8 files changed, 191 insertions(+), 94 deletions(-) diff --git a/src/qailo/mps/apply.py b/src/qailo/mps/apply.py index 52bdbf5..ceef6cd 100644 --- a/src/qailo/mps/apply.py +++ b/src/qailo/mps/apply.py @@ -1,11 +1,16 @@ +from __future__ import annotations + from copy import deepcopy +from typing import Iterable, Sequence + +import numpy.typing as npt from ..operator import type as op from ..operator.swap import swap from . import type as mps -def _swap_tensors(m, s): +def _swap_tensors(m: mps.mps, s: int) -> None: """ swap neighboring two tensors at s and s+1 """ @@ -16,7 +21,7 @@ def _swap_tensors(m, s): m.t2q[s], m.t2q[s + 1] = p1, p0 -def _move_qubit(m, p, s): +def _move_qubit(m: mps.mps, p: int, s: int) -> None: if m.q2t[p] != s: # print(f"moving qubit {p} at {m.q2t[p]} to {s}") for u in range(m.q2t[p], s): @@ -27,7 +32,7 @@ def _move_qubit(m, p, s): _swap_tensors(m, u - 1) -def _apply(m, p, pos=None): +def _apply(m: mps.mps, p: npt.NDArray, pos: Sequence[int] | None = None) -> mps.mps: assert op.is_operator(p) n = mps.num_qubits(m) if pos is None: @@ -50,15 +55,19 @@ def _apply(m, p, pos=None): return m -def apply(m, p, pos=None): +def apply(m: mps.mps, p: npt.NDArray, pos: Sequence[int] | None = None) -> mps.mps: return _apply(deepcopy(m), p, pos) -def _apply_seq(m, seq): +def _apply_seq( + m: mps.mps, seq: Iterable[tuple[npt.NDArray, Sequence[int] | None]] +) -> mps.mps: for p, qubit in seq: _apply(m, p, qubit) return m -def apply_seq(m, seq): +def apply_seq( + m: mps.mps, seq: Iterable[tuple[npt.NDArray, Sequence[int] | None]] +) -> mps.mps: return _apply_seq(deepcopy(m), seq) diff --git a/src/qailo/mps/mps_c.py b/src/qailo/mps/mps_c.py index 27abc7a..029c5a9 100644 --- a/src/qailo/mps/mps_c.py +++ b/src/qailo/mps/mps_c.py @@ -1,9 +1,13 @@ +from __future__ import annotations + from copy import deepcopy import numpy as np +import numpy.typing as npt from ..operator import type as op -from .svd import tensor_svd +from ..typeutil import eincheck as ec +from .svd import LegPartition, tensor_svd from .type import mps @@ -22,7 +26,11 @@ class canonical_mps(mps): tensors [cp(1)+1...n-1]: right canonical """ - def __init__(self, tensors, nkeep=None): + tensors: list[npt.NDArray] + nkeep: int | None + cp: list[int] + + def __init__(self, tensors: list[npt.NDArray], nkeep: int | None = None) -> None: assert isinstance(tensors, list) self.tensors = deepcopy(tensors) self.nkeep = nkeep @@ -31,38 +39,38 @@ def __init__(self, tensors, nkeep=None): self.t2q = list(range(n)) self.cp = [0, n - 1] - def _num_qubits(self): + def _num_qubits(self) -> int: return len(self.tensors) - def _norm(self): + def _norm(self) -> float: A = np.identity(1) for t in range(self._num_qubits()): - A = np.einsum("ij,jkl->ikl", A, self._tensor(t)) - A = np.einsum("ijk,ijl->kl", A, self._tensor(t).conj()) + A = ec.einsum_cast("ij,jkl->ikl", A, self._tensor(t)) + A = ec.einsum_cast("ijk,ijl->kl", A, self._tensor(t).conj()) ret = np.sqrt(np.trace(A)) assert isinstance(ret, float) return ret - def _state_vector(self): + def _state_vector(self) -> npt.NDArray: n = self._num_qubits() v = self._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, self._tensor(t), ss1) + v = ec.einsum_cast(v, ss0, self._tensor(t), ss1) v = v.reshape((2,) * n) - return np.einsum(v, self.t2q).reshape((2,) * n + (1,)) + return ec.einsum_cast(v, self.t2q).reshape((2,) * n + (1,)) - def _tensor(self, t): + def _tensor(self, t: int) -> npt.NDArray: return self.tensors[t] - def _canonicalize(self, p0, p1=None): + def _canonicalize(self, p0: int, p1: int | None) -> None: p1 = p0 if p1 is None else p1 n = len(self.tensors) assert 0 <= p0 and p0 <= p1 and p1 < n if self.cp[0] < p0: for t in range(self.cp[0], p0): - L, R = tensor_svd(self.tensors[t], [[0, 1], [2]], "left") + L, R = tensor_svd(self.tensors[t], LegPartition([0, 1], [2]), "left") self.tensors[t] = L self.tensors[t + 1] = np.einsum( R, [0, 3], self.tensors[t + 1], [3, 1, 2] @@ -71,14 +79,14 @@ def _canonicalize(self, p0, p1=None): self.cp[1] = max(p0, self.cp[1]) if self.cp[1] > p1: for t in range(self.cp[1], p1, -1): - L, R = tensor_svd(self.tensors[t], [[0], [1, 2]], "right") + L, R = tensor_svd(self.tensors[t], LegPartition([0], [1, 2]), "right") self.tensors[t - 1] = np.einsum( self.tensors[t - 1], [0, 1, 3], L, [3, 2] ) self.tensors[t] = R self.cp[1] = p1 - def _is_canonical(self): + def _is_canonical(self) -> bool: # tensor shape n = len(self.tensors) dims = [] @@ -102,23 +110,27 @@ def _is_canonical(self): assert self.cp[0] in range(n) assert self.cp[1] in range(n) for t in range(0, self.cp[0]): - A = np.einsum(self.tensors[t], [2, 3, 1], self.tensors[t].conj(), [2, 3, 0]) + A = ec.einsum_cast( + self.tensors[t], [2, 3, 1], self.tensors[t].conj(), [2, 3, 0] + ) assert np.allclose(A, np.identity(A.shape[0])) for t in range(self.cp[1] + 1, n): - A = np.einsum(self.tensors[t], [1, 3, 2], self.tensors[t].conj(), [0, 3, 2]) + A = ec.einsum_cast( + self.tensors[t], [1, 3, 2], self.tensors[t].conj(), [0, 3, 2] + ) assert np.allclose(A, np.identity(A.shape[0])) return True - def _apply_one(self, p, s): + def _apply_one(self, p: npt.NDArray, s: int) -> None: """ apply 1-qubit operator on tensor at s """ assert op.num_qubits(p) == 1 - self.tensors[s] = np.einsum(self.tensors[s], [0, 3, 2], p, [1, 3]) + self.tensors[s] = ec.einsum_cast(self.tensors[s], [0, 3, 2], p, [1, 3]) self.cp[0] = min(self.cp[0], s) self.cp[1] = max(self.cp[1], s) - def _apply_two(self, p, s, reverse=False): + def _apply_two(self, p: npt.NDArray, s: int, reverse: bool = False) -> None: """ apply 2-qubit operator on neighboring tensors at s and s+1 """ @@ -126,11 +138,11 @@ def _apply_two(self, p, s, reverse=False): self._canonicalize(s, s + 1) t0 = self.tensors[s] t1 = self.tensors[s + 1] - t = np.einsum(t0, [0, 1, 4], t1, [4, 2, 3]) + t = ec.einsum_cast(t0, [0, 1, 4], t1, [4, 2, 3]) if not reverse: - t = np.einsum(t, [0, 4, 5, 3], p, [1, 2, 4, 5]) + t = ec.einsum_cast(t, [0, 4, 5, 3], p, [1, 2, 4, 5]) else: - t = np.einsum(t, [0, 4, 5, 3], p, [2, 1, 5, 4]) - L, R = tensor_svd(t, [[0, 1], [2, 3]], nkeep=self.nkeep) + t = ec.einsum_cast(t, [0, 4, 5, 3], p, [2, 1, 5, 4]) + L, R = tensor_svd(t, LegPartition([0, 1], [2, 3]), nkeep=self.nkeep) self.tensors[s] = L self.tensors[s + 1] = R diff --git a/src/qailo/mps/mps_p.py b/src/qailo/mps/mps_p.py index 40d6c7b..958217b 100644 --- a/src/qailo/mps/mps_p.py +++ b/src/qailo/mps/mps_p.py @@ -1,10 +1,14 @@ +from __future__ import annotations + from copy import deepcopy import numpy as np +import numpy.typing as npt -from ..mps.svd import tensor_svd +from ..mps.svd import LegPartition, tensor_svd from ..mps.type import mps from ..operator import type as op +from ..typeutil import eincheck as ec from .projector import projector @@ -23,7 +27,12 @@ class projector_mps(mps): tensors [cp(1)+1...n-1]: bottom canonical """ - def __init__(self, tensors, nkeep=None): + tensors: list[npt.NDArray] + nkeep: int | None + cp: list[int] + env: list[npt.NDArray | None] + + def __init__(self, tensors: list[npt.NDArray], nkeep: int | None = None) -> None: assert isinstance(tensors, list) self.tensors = deepcopy(tensors) self.nkeep = nkeep @@ -35,17 +44,19 @@ def __init__(self, tensors, nkeep=None): # 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): + def _num_qubits(self) -> int: return len(self.tensors) - def _norm(self): + def _norm(self) -> float: A = np.identity(1) for t in range(self._num_qubits()): + A = ec.einsum_cast("ij,jkl->ikl", A, self._tensor(t)) + A = ec.einsum_cast("ijk,ijl->kl", A, self._tensor(t).conj()) ret = np.sqrt(np.trace(A)) assert isinstance(ret, float) return ret - def _state_vector(self): + def _state_vector(self) -> npt.NDArray: n = self._num_qubits() v = self._tensor(0) for t in range(1, n): @@ -53,31 +64,35 @@ def _state_vector(self): ss1 = [t + 3, t + 1, t + 2] v = np.einsum(v, ss0, self._tensor(t), ss1) v = v.reshape((2,) * n) - return np.einsum(v, self.t2q).reshape((2,) * n + (1,)) + return ec.einsum_cast(v, self.t2q).reshape((2,) * n + (1,)) - def _tensor(self, t): + def _tensor(self, t: int) -> npt.NDArray: return self.tensors[t] - def _canonicalize(self, p0, p1=None): + def _canonicalize(self, p0: int, p1: int | None) -> None: p1 = p0 if p1 is None else p1 n = len(self.tensors) 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.env[t], [0, 3], self.tensors[t], [3, 1, 2]) - _, self.env[t + 1] = tensor_svd(A, [[0, 1], [2]], "left") + env = self.env[t] + assert env is not None + A = ec.einsum_cast(env, [0, 3], self.tensors[t], [3, 1, 2]) + _, self.env[t + 1] = tensor_svd(A, LegPartition([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.env[t + 1], [3, 2]) - self.env[t], _ = tensor_svd(A, [[0], [1, 2]], "right") + env = self.env[t + 1] + assert env is not None + A = ec.einsum_cast(self.tensors[t], [0, 1, 3], env, [3, 2]) + self.env[t], _ = tensor_svd(A, LegPartition([0], [1, 2]), "right") self.cp[1] = p1 - def _is_canonical(self): + def _is_canonical(self) -> bool: # tensor shape n = len(self.tensors) - dims = [] + dims: list[int] = [] assert self.tensors[0].shape[0] == 1 dims.append(self.tensors[0].shape[0]) for t in range(1, n - 1): @@ -99,41 +114,49 @@ def _is_canonical(self): 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.tensors[t], [3, 1, 2]) - A = np.einsum(A, [2, 3, 1], self.tensors[t].conj(), [2, 3, 0]) - B = np.einsum(self.env[t + 1], [2, 1], self.env[t + 1].conj(), [2, 0]) + A = ec.einsum_cast(A, [0, 3], self.tensors[t], [3, 1, 2]) + A = ec.einsum_cast(A, [2, 3, 1], self.tensors[t].conj(), [2, 3, 0]) + env = self.env[t + 1] + assert env is not None + B = ec.einsum_cast(env, [2, 1], env.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.env[t], [0, 2], self.env[t].conj(), [1, 2]) + A = ec.einsum_cast(self.tensors[t], [0, 1, 3], A, [3, 2]) + A = ec.einsum_cast(self.tensors[t].conj(), [1, 2, 3], A, [0, 2, 3]) + env = self.env[t] + assert env is not None + B = ec.einsum_cast(env, [0, 2], env.conj(), [1, 2]) assert np.allclose(A, B) return True - def _apply_one(self, p, s): + def _apply_one(self, p: npt.NDArray, s: int) -> None: assert op.num_qubits(p) == 1 self.tensors[s] = np.einsum(self.tensors[s], [0, 3, 2], p, [1, 3]) self.cp[0] = min(self.cp[0], s) self.cp[1] = max(self.cp[1], s) - def _apply_two(self, p, s, reverse=False): + def _apply_two(self, p: npt.NDArray, s: int, reverse: bool = False) -> None: """ apply 2-qubit operator on neighboring tensors, s and s+1 """ self._canonicalize(s, s + 1) t0 = self.tensors[s] t1 = self.tensors[s + 1] - p0, p1 = tensor_svd(p, [[0, 2], [1, 3]]) + p0, p1 = tensor_svd(p, LegPartition([0, 2], [1, 3])) if not reverse: - t0 = np.einsum(t0, [0, 4, 3], p0, [1, 4, 2]) - t1 = np.einsum(t1, [0, 4, 3], p1, [1, 2, 4]) + t0 = ec.einsum_cast(t0, [0, 4, 3], p0, [1, 4, 2]) + t1 = ec.einsum_cast(t1, [0, 4, 3], p1, [1, 2, 4]) 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.env[s], [0, 4], t0, [4, 1, 2, 3]) - tt1 = np.einsum(t1, [0, 1, 2, 4], self.env[s + 2], [4, 3]) + t0 = ec.einsum_cast(t0, [0, 4, 3], p1, [2, 1, 4]) + t1 = ec.einsum_cast(t1, [0, 4, 3], p0, [2, 4, 1]) + env = self.env[s] + assert env is not None + tt0 = ec.einsum_cast(env, [0, 4], t0, [4, 1, 2, 3]) + env = self.env[s + 2] + assert env is not None + tt1 = ec.einsum_cast(t1, [0, 1, 2, 4], env, [4, 3]) _, PLh, PR = projector(tt0, [0, 1, 4, 5], tt1, [5, 4, 2, 3], nkeep=self.nkeep) - self.tensors[s] = np.einsum(t0, [0, 1, 3, 4], PR, [3, 4, 2]) - self.tensors[s + 1] = np.einsum(PLh, [3, 4, 0], t1, [4, 3, 1, 2]) + self.tensors[s] = ec.einsum_cast(t0, [0, 1, 3, 4], PR, [3, 4, 2]) + self.tensors[s + 1] = ec.einsum_cast(PLh, [3, 4, 0], t1, [4, 3, 1, 2]) diff --git a/src/qailo/mps/product_state.py b/src/qailo/mps/product_state.py index 230135e..16ca942 100644 --- a/src/qailo/mps/product_state.py +++ b/src/qailo/mps/product_state.py @@ -1,38 +1,52 @@ +from __future__ import annotations + +from typing import Iterable + +import numpy as np +import numpy.typing as npt + 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.vector import vector from .mps_c import canonical_mps -from .svd import tensor_svd -from .type import is_mps +from .svd import LegPartition, tensor_svd +from .type import is_mps, mps -def tensor_decomposition(v, nkeep=None, tol=1e-12): +def tensor_decomposition( + v: mps | npt.NDArray, nkeep: int | None = None, tol: float = 1e-12 +) -> list[npt.NDArray]: if is_mps(v): return [v._tensor(s) for s in range(num_qubits(v))] else: + assert isinstance(v, np.ndarray) n = num_qubits(v) w = vector(v).reshape((1, 2**n)) - tensors = [] + tensors: list[npt.NDArray] = [] for t in range(n - 1): dims = w.shape w = w.reshape(dims[0], 2, dims[1] // 2) - t, w = tensor_svd(w, [[0, 1], [2]], "left", nkeep, tol) + t, w = tensor_svd(w, LegPartition([0, 1], [2]), "left", nkeep, tol) tensors.append(t) tensors.append(w.reshape(w.shape + (1,))) return tensors -def product_state(states, nkeep=None, mps=canonical_mps): - tensors = [] +def product_state( + states: Iterable[mps | npt.NDArray], + nkeep: int | None = None, + mps: type[mps] = canonical_mps, +) -> mps: + tensors: list[npt.NDArray] = [] for s in states: tensors = tensors + tensor_decomposition(s, nkeep) return mps(tensors, nkeep) -def zero(n=1, nkeep=None, mps=canonical_mps): +def zero(n: int = 1, nkeep: int | None = None, mps: type[mps] = canonical_mps) -> mps: return product_state([sv_zero()] * n, nkeep, mps) -def one(n=1, nkeep=None, mps=canonical_mps): +def one(n: int = 1, nkeep: int | None = None, mps: type[mps] = canonical_mps) -> mps: return product_state([sv_one()] * n, nkeep, mps) diff --git a/src/qailo/mps/projector.py b/src/qailo/mps/projector.py index b8b14d8..36d909f 100644 --- a/src/qailo/mps/projector.py +++ b/src/qailo/mps/projector.py @@ -1,9 +1,13 @@ +from __future__ import annotations + import numpy as np +import numpy.typing as npt from ..mps.svd import compact_svd +from ..typeutil import eincheck as ec -def normalize_ss(ss0, ss1): +def normalize_ss(ss0: list[int], ss1: list[int]) -> tuple[list[int], list[int]]: """ Order the subscripts that are not contracted from 0. Subscripts that will be contracted are assigned a larger value. @@ -28,17 +32,19 @@ def normalize_ss(ss0, ss1): return ss0_out, ss1_out -def collect_legs(ss0, ss1): +def collect_legs( + ss0: list[int], ss1: list[int] +) -> tuple[list[int], list[int], list[int], list[int]]: """ legs0L: list of legs in T0 that are NOT contracted legs0R: list of legs in T0 that will be contracted legs1L: list of legs in T1 that will be contracted legs1R: list of legs in T1 that are NOT contracted """ - legs0L = [] - legs0R = [] - legs1L = [] - legs1R = [] + legs0L: list[int] = [] + legs0R: list[int] = [] + legs1L: list[int] = [] + legs1R: list[int] = [] for k in range(len(ss0)): if ss0[k] not in ss1: legs0L.append(k) @@ -52,7 +58,14 @@ def collect_legs(ss0, ss1): return legs0L, legs0R, legs1L, legs1R -def projector(T0, ss0_in, T1, ss1_in, nkeep=None, tol=1e-12): +def projector( + T0: npt.NDArray, + ss0_in: list[int], + T1: npt.NDArray, + ss1_in: list[int], + nkeep: int | None = None, + tol: float = 1e-12, +) -> tuple[npt.NDArray, npt.NDArray, npt.NDArray]: ss0, ss1 = normalize_ss(ss0_in, ss1_in) legs0L, legs0R, legs1L, legs1R = collect_legs(ss0, ss1) dim0L = np.prod([T0.shape[i] for i in legs0L]) @@ -60,13 +73,13 @@ def projector(T0, ss0_in, T1, ss1_in, nkeep=None, tol=1e-12): dims1L = [T1.shape[i] for i in legs1L] dim1R = np.prod([T1.shape[i] for i in legs1R]) assert len(dims0R) == len(dims1L) - TT0 = np.einsum(T0, ss0).reshape([dim0L] + dims0R) - TT1 = np.einsum(T1, ss1).reshape([dim1R] + dims0R) + TT0 = ec.einsum_cast(T0, ss0).reshape([dim0L] + dims0R) + TT1 = ec.einsum_cast(T1, ss1).reshape([dim1R] + dims0R) ss_sum = list(range(2, len(dims0R) + 2)) - A = np.einsum(TT0, [0] + ss_sum, TT1, [1] + ss_sum) + A = ec.einsum_cast(TT0, [0] + ss_sum, TT1, [1] + ss_sum) S, U, V = compact_svd(A, nkeep=nkeep, tol=tol) - U = np.einsum(U, [0, 1], np.sqrt(1 / S), [1], [0, 1]) - V = np.einsum(V, [0, 1], np.sqrt(1 / S), [1], [0, 1]) - PL = np.einsum(TT0.conj(), [0] + ss_sum, U, [0, max(ss_sum) + 1]) - PR = np.einsum(TT1, [0] + ss_sum, V, [0, max(ss_sum) + 1]) + U = ec.einsum_cast(U, [0, 1], np.sqrt(1 / S), [1], [0, 1]) + V = ec.einsum_cast(V, [0, 1], np.sqrt(1 / S), [1], [0, 1]) + PL = ec.einsum_cast(TT0.conj(), [0] + ss_sum, U, [0, max(ss_sum) + 1]) + PR = ec.einsum_cast(TT1, [0] + ss_sum, V, [0, max(ss_sum) + 1]) return S, PL.conj(), PR diff --git a/src/qailo/mps/state_vector.py b/src/qailo/mps/state_vector.py index ad0d653..71b7822 100644 --- a/src/qailo/mps/state_vector.py +++ b/src/qailo/mps/state_vector.py @@ -1,6 +1,10 @@ -from .type import is_mps +from __future__ import annotations +import numpy.typing as npt -def state_vector(m): +from .type import is_mps, mps + + +def state_vector(m: mps) -> npt.NDArray: assert is_mps(m) return m._state_vector() diff --git a/src/qailo/mps/svd.py b/src/qailo/mps/svd.py index deb701d..22de2b8 100644 --- a/src/qailo/mps/svd.py +++ b/src/qailo/mps/svd.py @@ -1,7 +1,22 @@ +from __future__ import annotations + +from typing import NamedTuple + import numpy as np +import numpy.typing as npt +from typing_extensions import Literal + +from ..typeutil import eincheck as ec + + +class LegPartition(NamedTuple): + leg0: list[int] + leg1: list[int] -def compact_svd(A, nkeep=None, tol=1e-12): +def compact_svd( + A: npt.NDArray, nkeep: int | None = None, tol: float = 1e-12 +) -> tuple[npt.NDArray, npt.NDArray, npt.NDArray]: assert A.ndim == 2 U, S, Vh = np.linalg.svd(A, full_matrices=False) V = Vh.conj().T @@ -10,7 +25,13 @@ def compact_svd(A, nkeep=None, tol=1e-12): return S[:dimS], U[:, :dimS], V[:, :dimS] -def tensor_svd(T, partition, canonical="center", nkeep=None, tol=1e-12): +def tensor_svd( + T: npt.NDArray, + partition: LegPartition, + canonical: Literal["center", "left", "right"] = "center", + nkeep: int | None = None, + tol: float = 1e-12, +) -> tuple[npt.NDArray, npt.NDArray]: legsL = len(partition[0]) legsR = len(partition[1]) assert T.ndim == legsL + legsR @@ -24,12 +45,12 @@ def tensor_svd(T, partition, canonical="center", nkeep=None, tol=1e-12): L = U R = V.conj().T if canonical == "center": - L = np.einsum("ij,j->ij", L, np.sqrt(S)) - R = np.einsum("i,ij->ij", np.sqrt(S), R) + L = ec.einsum_cast("ij,j->ij", L, np.sqrt(S)) + R = ec.einsum_cast("i,ij->ij", np.sqrt(S), R) elif canonical == "left": - R = np.einsum("i,ij->ij", S, R) + R = ec.einsum_cast("i,ij->ij", S, R) elif canonical == "right": - L = np.einsum("ij,j->ij", L, S) + L = ec.einsum_cast("ij,j->ij", L, S) else: raise ValueError L = L.reshape(dimsL + [S.shape[0]]) diff --git a/src/qailo/mps/type.py b/src/qailo/mps/type.py index 6a71d3b..c41d9e7 100644 --- a/src/qailo/mps/type.py +++ b/src/qailo/mps/type.py @@ -1,6 +1,7 @@ from __future__ import annotations from abc import ABC, abstractmethod +from typing import Any import numpy.typing as npt from typing_extensions import TypeGuard @@ -38,13 +39,13 @@ def _apply_one(self, p: npt.NDArray, s: int) -> None: ... def _apply_two(self, p: npt.NDArray, s: int, reverse: bool = False) -> None: ... -def is_canonical(m): +def is_canonical(m: mps) -> bool: return m._is_canonical() -def is_mps(m): +def is_mps(m: Any) -> TypeGuard[mps]: return isinstance(m, mps) -def num_qubits(m): +def num_qubits(m: mps) -> int: return len(m.q2t) From 9d342146f0132117b8504a4175a591a0d8fd7538 Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 06:01:16 +0900 Subject: [PATCH 23/40] :art: Improve flow --- src/qailo/dispatch.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/qailo/dispatch.py b/src/qailo/dispatch.py index f7bd8e5..c2f9f5d 100644 --- a/src/qailo/dispatch.py +++ b/src/qailo/dispatch.py @@ -27,9 +27,11 @@ def apply_seq(v, seq): def norm(v): if sv.is_state_vector(v): - return np.linalg.norm(v) + return float(np.linalg.norm(v)) elif mps.is_mps(v): return v._norm() + else: + assert False def num_qubits(v): From d0033fe5b08930a75ffc3ef8095c45223bbadc54 Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 06:06:33 +0900 Subject: [PATCH 24/40] :wheelchair: Add OPSeqElement --- src/qailo/mps/apply.py | 9 +++------ src/qailo/state_vector/apply.py | 3 ++- src/qailo/util/helpertype.py | 10 ++++++++++ 3 files changed, 15 insertions(+), 7 deletions(-) create mode 100644 src/qailo/util/helpertype.py diff --git a/src/qailo/mps/apply.py b/src/qailo/mps/apply.py index ceef6cd..b0906df 100644 --- a/src/qailo/mps/apply.py +++ b/src/qailo/mps/apply.py @@ -7,6 +7,7 @@ from ..operator import type as op from ..operator.swap import swap +from ..util.helpertype import OPSeqElement from . import type as mps @@ -59,15 +60,11 @@ def apply(m: mps.mps, p: npt.NDArray, pos: Sequence[int] | None = None) -> mps.m return _apply(deepcopy(m), p, pos) -def _apply_seq( - m: mps.mps, seq: Iterable[tuple[npt.NDArray, Sequence[int] | None]] -) -> mps.mps: +def _apply_seq(m: mps.mps, seq: Iterable[OPSeqElement]) -> mps.mps: for p, qubit in seq: _apply(m, p, qubit) return m -def apply_seq( - m: mps.mps, seq: Iterable[tuple[npt.NDArray, Sequence[int] | None]] -) -> mps.mps: +def apply_seq(m: mps.mps, seq: Iterable[OPSeqElement]) -> mps.mps: return _apply_seq(deepcopy(m), seq) diff --git a/src/qailo/state_vector/apply.py b/src/qailo/state_vector/apply.py index e3d4d6a..5596e73 100644 --- a/src/qailo/state_vector/apply.py +++ b/src/qailo/state_vector/apply.py @@ -6,6 +6,7 @@ from ..operator import type as op from ..typeutil import eincheck as ec +from ..util.helpertype import OPSeqElement from . import type as sv @@ -33,7 +34,7 @@ def apply( def apply_seq( v: npt.NDArray, - seq: Iterable[tuple[npt.NDArray, Sequence[int]]], + seq: Iterable[OPSeqElement], ) -> npt.NDArray: for p, qubit in seq: v = apply(v, p, qubit) diff --git a/src/qailo/util/helpertype.py b/src/qailo/util/helpertype.py new file mode 100644 index 0000000..5e9bb8b --- /dev/null +++ b/src/qailo/util/helpertype.py @@ -0,0 +1,10 @@ +from __future__ import annotations + +from typing import NamedTuple + +import numpy.typing as npt + + +class OPSeqElement(NamedTuple): + p: npt.NDArray + qubit: list[int] From 213bdfecb5a75b79070f24628fa49d34330a679a Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 06:07:01 +0900 Subject: [PATCH 25/40] :construction: Annotate dispatch.py --- src/qailo/dispatch.py | 45 +++++++++++++++++++++++++++++++++++++------ 1 file changed, 39 insertions(+), 6 deletions(-) diff --git a/src/qailo/dispatch.py b/src/qailo/dispatch.py index c2f9f5d..2f02c06 100644 --- a/src/qailo/dispatch.py +++ b/src/qailo/dispatch.py @@ -1,11 +1,32 @@ +from __future__ import annotations + +from typing import Container, Iterable, overload + import numpy as np +import numpy.typing as npt from . import mps from . import operator as op from . import state_vector as sv +from .mps import type as mpstype +from .util.helpertype import OPSeqElement + + +@overload +def apply( + v: npt.NDArray, p: npt.NDArray, pos: list[int] | None = None +) -> npt.NDArray: ... + +@overload +def apply( + v: mpstype.mps, p: npt.NDArray, pos: list[int] | None = None +) -> mpstype.mps: ... -def apply(v, p, pos=None): + +def apply( + v: npt.NDArray | mpstype.mps, p: npt.NDArray, pos: list[int] | None = None +) -> npt.NDArray | mpstype.mps: if sv.is_state_vector(v): v = sv.apply(v, p, pos) elif mps.is_mps(v): @@ -15,7 +36,17 @@ def apply(v, p, pos=None): return v -def apply_seq(v, seq): +@overload +def apply_seq(v: npt.NDArray, seq: Iterable[OPSeqElement]) -> npt.NDArray: ... + + +@overload +def apply_seq(v: mpstype.mps, seq: Iterable[OPSeqElement]) -> mpstype.mps: ... + + +def apply_seq( + v: npt.NDArray | mpstype.mps, seq: Iterable[OPSeqElement] +) -> npt.NDArray | mpstype.mps: if sv.is_state_vector(v): v = sv.apply_seq(v, seq) elif mps.is_mps(v): @@ -25,7 +56,7 @@ def apply_seq(v, seq): return v -def norm(v): +def norm(v: npt.NDArray | mpstype.mps) -> float: if sv.is_state_vector(v): return float(np.linalg.norm(v)) elif mps.is_mps(v): @@ -34,7 +65,7 @@ def norm(v): assert False -def num_qubits(v): +def num_qubits(v: npt.NDArray | mpstype.mps) -> int: if sv.is_state_vector(v): return sv.num_qubits(v) elif op.is_operator(v): @@ -44,7 +75,9 @@ def num_qubits(v): assert False -def probability(v, pos=None): +def probability( + v: npt.NDArray | mpstype.mps, pos: Container[int] | None = None +) -> npt.NDArray: if sv.is_state_vector(v): return sv.probability(v, pos) elif mps.is_mps(v): @@ -52,7 +85,7 @@ def probability(v, pos=None): assert False -def vector(v, c=None): +def vector(v: npt.NDArray | mpstype.mps, c: list | None = None) -> npt.NDArray: if sv.is_state_vector(v): return sv.vector(v, c) elif mps.is_mps(v): From c7c3d9dd5d9f02ad4ab06a49161bb80748c026cb Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 06:39:44 +0900 Subject: [PATCH 26/40] :construction: Annotate alg module --- src/qailo/alg/qft.py | 21 ++++++++++++--------- src/qailo/alg/qpe.py | 5 ++++- 2 files changed, 16 insertions(+), 10 deletions(-) diff --git a/src/qailo/alg/qft.py b/src/qailo/alg/qft.py index adbd967..9f75174 100644 --- a/src/qailo/alg/qft.py +++ b/src/qailo/alg/qft.py @@ -1,36 +1,39 @@ # ref: https://learn.qiskit.org/course/ch-algorithms/quantum-fourier-transform +from __future__ import annotations + import numpy as np import qailo as q +from qailo.util.helpertype import OPSeqElement -def qft_rotations_seq(n): - seq = [] +def qft_rotations_seq(n: int) -> list[OPSeqElement]: + seq: list[OPSeqElement] = [] if n == 0: return seq n -= 1 # print(f"H on [{n}]") - seq.append([q.op.h(), [n]]) + seq.append(OPSeqElement(q.op.h(), [n])) for p in range(n): # print(f"CP(pi/{2**(n-p)} on [{p}, {n}]") - seq.append([q.op.cp(np.pi / 2 ** (n - p)), [p, n]]) + seq.append(OPSeqElement(q.op.cp(np.pi / 2 ** (n - p)), [p, n])) seq += qft_rotations_seq(n) return seq -def swap_registers_seq(n): - seq = [] +def swap_registers_seq(n: int) -> list[OPSeqElement]: + seq: list[OPSeqElement] = [] for p in range(n // 2): # print(f"swap on [{p}, {n-p-1}]") - seq.append([q.op.swap(), [p, n - p - 1]]) + seq.append(OPSeqElement(q.op.swap(), [p, n - p - 1])) return seq -def qft_seq(n): +def qft_seq(n: int) -> list[OPSeqElement]: """QFT on the first n qubits in circuit""" return qft_rotations_seq(n) + swap_registers_seq(n) -def inverse_qft_seq(n): +def inverse_qft_seq(n: int) -> list[OPSeqElement]: return q.op.inverse_seq(qft_seq(n)) diff --git a/src/qailo/alg/qpe.py b/src/qailo/alg/qpe.py index dc704c7..6f880fa 100644 --- a/src/qailo/alg/qpe.py +++ b/src/qailo/alg/qpe.py @@ -1,11 +1,14 @@ +from __future__ import annotations + from copy import deepcopy import numpy as np +import numpy.typing as npt import qailo as q -def qpe(n, u, v): +def qpe(n: int, u: npt.NDArray, v: npt.NDArray) -> npt.NDArray: m = q.num_qubits(u) w = deepcopy(v) assert q.num_qubits(w) == m + n From 1a9f91ef467e6a829bdc5aaee58aa8d301f81464 Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 06:15:34 +0900 Subject: [PATCH 27/40] :heavy_plus_sign: Add typing-extensions --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index b052bb9..5a175b8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -26,7 +26,7 @@ write_to = "src/qailo/_version.py" package-dir = { "" = "src" } [project.optional-dependencies] -dev = ["pytest", "black", "ruff"] +dev = ["pytest", "black", "ruff", "typing-extensions"] [tool.black] target-version = ["py38"] From d4c74979440ea5a3f3711f7c6e6c155ca81f283b Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 06:19:01 +0900 Subject: [PATCH 28/40] :fire: Remove OPAutomaton --- src/qailo/operator/controlled.py | 15 ++++++++------- src/qailo/operator/inverse.py | 9 +++++---- src/qailo/operator/type.py | 7 +------ 3 files changed, 14 insertions(+), 17 deletions(-) diff --git a/src/qailo/operator/controlled.py b/src/qailo/operator/controlled.py index 23d8e15..caedb46 100644 --- a/src/qailo/operator/controlled.py +++ b/src/qailo/operator/controlled.py @@ -3,9 +3,10 @@ import numpy as np import numpy.typing as npt +from ..util.helpertype import OPSeqElement from .matrix import matrix from .one_qubit import p, x, z -from .type import OPAutomaton, is_operator, num_qubits +from .type import is_operator, num_qubits def controlled(u: npt.NDArray) -> npt.NDArray: @@ -70,16 +71,16 @@ def control_end(u: npt.NDArray) -> npt.NDArray: return op.reshape((2,) * n + (4,) + (2,) * m) -def controlled_seq(u: npt.NDArray, pos: list[int]) -> list[OPAutomaton]: +def controlled_seq(u: npt.NDArray, pos: list[int]) -> list[OPSeqElement]: n = len(pos) - seq: list[OPAutomaton] = [] + seq: list[OPSeqElement] = [] if n <= 1: raise ValueError elif n == 2: - seq.append(OPAutomaton(controlled(u), pos)) + seq.append(OPSeqElement(controlled(u), pos)) else: - seq.append(OPAutomaton(control_begin(), [pos[0], pos[1]])) + seq.append(OPSeqElement(control_begin(), [pos[0], pos[1]])) for i in range(1, n - 2): - seq.append(OPAutomaton(control_propagate(), [pos[i], pos[i + 1]])) - seq.append(OPAutomaton(control_end(u), [pos[-2], pos[-1]])) + seq.append(OPSeqElement(control_propagate(), [pos[i], pos[i + 1]])) + seq.append(OPSeqElement(control_end(u), [pos[-2], pos[-1]])) return seq diff --git a/src/qailo/operator/inverse.py b/src/qailo/operator/inverse.py index a325efe..a8907dc 100644 --- a/src/qailo/operator/inverse.py +++ b/src/qailo/operator/inverse.py @@ -2,13 +2,14 @@ from typing import Reversible +from ..util.helpertype import OPSeqElement from .hconj import hconj -from .type import OPAutomaton, is_operator +from .type import is_operator -def inverse_seq(seq: Reversible[OPAutomaton]) -> list[OPAutomaton]: - res: list[OPAutomaton] = [] +def inverse_seq(seq: Reversible[OPSeqElement]) -> list[OPSeqElement]: + res: list[OPSeqElement] = [] for p, pos in reversed(seq): assert is_operator(p) - res.append(OPAutomaton(hconj(p), pos)) + res.append(OPSeqElement(hconj(p), pos)) return res diff --git a/src/qailo/operator/type.py b/src/qailo/operator/type.py index 36aee5a..1984bda 100644 --- a/src/qailo/operator/type.py +++ b/src/qailo/operator/type.py @@ -1,6 +1,6 @@ from __future__ import annotations -from typing import Any, NamedTuple +from typing import Any import numpy as np import numpy.typing as npt @@ -30,8 +30,3 @@ def num_qubits(v: npt.NDArray) -> int: elif is_operator(v): return v.ndim // 2 raise ValueError - - -class OPAutomaton(NamedTuple): - op: npt.NDArray - pos: list[int] From 4754b4170ea9a2b6ebd7cc65289989ff844685da Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 06:24:07 +0900 Subject: [PATCH 29/40] :bug: Add default arg --- src/qailo/mps/mps_c.py | 2 +- src/qailo/mps/mps_p.py | 2 +- src/qailo/mps/type.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/qailo/mps/mps_c.py b/src/qailo/mps/mps_c.py index 029c5a9..17434ea 100644 --- a/src/qailo/mps/mps_c.py +++ b/src/qailo/mps/mps_c.py @@ -64,7 +64,7 @@ def _state_vector(self) -> npt.NDArray: def _tensor(self, t: int) -> npt.NDArray: return self.tensors[t] - def _canonicalize(self, p0: int, p1: int | None) -> None: + def _canonicalize(self, p0: int, p1: int | None = None) -> None: p1 = p0 if p1 is None else p1 n = len(self.tensors) assert 0 <= p0 and p0 <= p1 and p1 < n diff --git a/src/qailo/mps/mps_p.py b/src/qailo/mps/mps_p.py index 958217b..cb79bdd 100644 --- a/src/qailo/mps/mps_p.py +++ b/src/qailo/mps/mps_p.py @@ -69,7 +69,7 @@ def _state_vector(self) -> npt.NDArray: def _tensor(self, t: int) -> npt.NDArray: return self.tensors[t] - def _canonicalize(self, p0: int, p1: int | None) -> None: + def _canonicalize(self, p0: int, p1: int | None = None) -> None: p1 = p0 if p1 is None else p1 n = len(self.tensors) assert 0 <= p0 and p0 <= p1 and p1 < n diff --git a/src/qailo/mps/type.py b/src/qailo/mps/type.py index c41d9e7..dd6338d 100644 --- a/src/qailo/mps/type.py +++ b/src/qailo/mps/type.py @@ -27,7 +27,7 @@ def _state_vector(self) -> npt.NDArray: ... def _tensor(self, t: int) -> npt.NDArray: ... @abstractmethod - def _canonicalize(self, p0: int, p1: int | None) -> None: ... + def _canonicalize(self, p0: int, p1: int | None = None) -> None: ... @abstractmethod def _is_canonical(self) -> bool: ... From 355d2bcc53cc73d0a848dce56624ea50a965767a Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 06:30:43 +0900 Subject: [PATCH 30/40] :art: Remove np.einsum --- src/qailo/mps/mps_c.py | 4 ++-- src/qailo/mps/mps_p.py | 4 ++-- src/qailo/mps/svd.py | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/qailo/mps/mps_c.py b/src/qailo/mps/mps_c.py index 17434ea..ec1534f 100644 --- a/src/qailo/mps/mps_c.py +++ b/src/qailo/mps/mps_c.py @@ -72,7 +72,7 @@ def _canonicalize(self, p0: int, p1: int | None = None) -> None: for t in range(self.cp[0], p0): L, R = tensor_svd(self.tensors[t], LegPartition([0, 1], [2]), "left") self.tensors[t] = L - self.tensors[t + 1] = np.einsum( + self.tensors[t + 1] = ec.einsum_cast( R, [0, 3], self.tensors[t + 1], [3, 1, 2] ) self.cp[0] = p0 @@ -80,7 +80,7 @@ def _canonicalize(self, p0: int, p1: int | None = None) -> None: if self.cp[1] > p1: for t in range(self.cp[1], p1, -1): L, R = tensor_svd(self.tensors[t], LegPartition([0], [1, 2]), "right") - self.tensors[t - 1] = np.einsum( + self.tensors[t - 1] = ec.einsum_cast( self.tensors[t - 1], [0, 1, 3], L, [3, 2] ) self.tensors[t] = R diff --git a/src/qailo/mps/mps_p.py b/src/qailo/mps/mps_p.py index cb79bdd..c251a7f 100644 --- a/src/qailo/mps/mps_p.py +++ b/src/qailo/mps/mps_p.py @@ -62,7 +62,7 @@ def _state_vector(self) -> npt.NDArray: 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, self._tensor(t), ss1) + v = ec.einsum_cast(v, ss0, self._tensor(t), ss1) v = v.reshape((2,) * n) return ec.einsum_cast(v, self.t2q).reshape((2,) * n + (1,)) @@ -133,7 +133,7 @@ def _is_canonical(self) -> bool: def _apply_one(self, p: npt.NDArray, s: int) -> None: assert op.num_qubits(p) == 1 - self.tensors[s] = np.einsum(self.tensors[s], [0, 3, 2], p, [1, 3]) + self.tensors[s] = ec.einsum_cast(self.tensors[s], [0, 3, 2], p, [1, 3]) self.cp[0] = min(self.cp[0], s) self.cp[1] = max(self.cp[1], s) diff --git a/src/qailo/mps/svd.py b/src/qailo/mps/svd.py index 22de2b8..c85cfc3 100644 --- a/src/qailo/mps/svd.py +++ b/src/qailo/mps/svd.py @@ -38,7 +38,7 @@ def tensor_svd( assert sorted(partition[0] + partition[1]) == list(range(T.ndim)) dimsL = [T.shape[i] for i in partition[0]] dimsR = [T.shape[i] for i in partition[1]] - m = np.einsum(T, partition[0] + partition[1]).reshape( + m = ec.einsum_cast(T, partition[0] + partition[1]).reshape( np.prod(dimsL), np.prod(dimsR) ) S, U, V = compact_svd(m, nkeep=nkeep, tol=tol) From 225d13b01d38cdfe720e6b3e596e0d891a5336bc Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 06:44:23 +0900 Subject: [PATCH 31/40] :art: Remove untyped list --- src/qailo/mps/mps_c.py | 2 +- src/qailo/state_vector/probability.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/qailo/mps/mps_c.py b/src/qailo/mps/mps_c.py index ec1534f..5abec66 100644 --- a/src/qailo/mps/mps_c.py +++ b/src/qailo/mps/mps_c.py @@ -89,7 +89,7 @@ def _canonicalize(self, p0: int, p1: int | None = None) -> None: def _is_canonical(self) -> bool: # tensor shape n = len(self.tensors) - dims = [] + dims: list[int] = [] assert self.tensors[0].shape[0] == 1 dims.append(self.tensors[0].shape[0]) for t in range(1, n - 1): diff --git a/src/qailo/state_vector/probability.py b/src/qailo/state_vector/probability.py index da63b5c..7dbe078 100644 --- a/src/qailo/state_vector/probability.py +++ b/src/qailo/state_vector/probability.py @@ -18,7 +18,7 @@ def probability(v: npt.NDArray, pos: Container[int] | None = None) -> npt.NDArra if pos is None: return abs(vector(w)) ** 2 else: - tpos = [] + tpos: list[int] = [] for k in range(num_qubits(w)): if k not in pos: tpos.append(k) From 1578b7af651da913a754748249c15623fc841f64 Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 06:47:44 +0900 Subject: [PATCH 32/40] :safety_vest: Annotate fallfack --- src/qailo/typeutil/eincheck.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/qailo/typeutil/eincheck.py b/src/qailo/typeutil/eincheck.py index 2c326b1..e352309 100644 --- a/src/qailo/typeutil/eincheck.py +++ b/src/qailo/typeutil/eincheck.py @@ -74,5 +74,5 @@ def einsum_cast( ) -> npt.NDArray: ... -def einsum_cast(*args, **kwargs): +def einsum_cast(*args, **kwargs) -> npt.NDArray: return np.asarray(np.einsum(*args, **kwargs)) From 817e52bb3976e17e7d48d0a9cb33a72141bf664a Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 07:01:47 +0900 Subject: [PATCH 33/40] :bug: Use relative import --- test/mps/test_projector.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/mps/test_projector.py b/test/mps/test_projector.py index 5ae40f5..632f621 100644 --- a/test/mps/test_projector.py +++ b/test/mps/test_projector.py @@ -42,7 +42,7 @@ def test_projector(): w2 = np.random.random((2, 2)) + 1.0j * np.random.random((2, 2)) p = np.random.random((2, 2, 2, 2)) + 1.0j * np.random.random((2, 2, 2, 2)) B = np.einsum(t0, [0, 4, 6], t1, [6, 5, 3], p, [1, 2, 4, 5]) - p0, p1 = q.mps.svd.tensor_svd(p, [[0, 2], [1, 3]]) + p0, p1 = tensor_svd(p, [[0, 2], [1, 3]]) assert np.allclose(np.einsum(p0, [0, 2, 4], p1, [4, 1, 3]), p) t0 = np.einsum(t0, [0, 4, 3], p0, [1, 4, 2]) t1 = np.einsum(t1, [0, 4, 3], p1, [1, 2, 4]) @@ -66,7 +66,7 @@ def test_projector(): w2 = np.random.random((2, 2)) + 1.0j * np.random.random((2, 2)) p = np.random.random((2, 2, 2, 2)) + 1.0j * np.random.random((2, 2, 2, 2)) B = np.einsum(t0, [0, 4, 6], t1, [6, 5, 3], p, [1, 2, 4, 5]) - p0, p1 = q.mps.svd.tensor_svd(p, [[0, 2], [1, 3]]) + p0, p1 = tensor_svd(p, [[0, 2], [1, 3]]) assert np.allclose(np.einsum(p0, [0, 2, 4], p1, [4, 1, 3]), p) t0 = np.einsum(t0, [0, 4, 3], p0, [1, 4, 2]) t1 = np.einsum(t1, [1, 4, 3], p1, [0, 2, 4]) From 91e390a196f8dcbe9882803638b67a87ae96bd5b Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 07:05:35 +0900 Subject: [PATCH 34/40] :white_check_mark: Update tests --- test/mps/test_apply.py | 20 +++++++++++++++++--- test/mps/test_canonical.py | 8 +++++++- test/mps/test_move.py | 10 ++++++++-- test/mps/test_projector.py | 10 ++++++---- test/mps/test_svd.py | 26 +++++++++++++------------- test/state_vector/test_apply_seq.py | 5 +++-- 6 files changed, 54 insertions(+), 25 deletions(-) diff --git a/test/mps/test_apply.py b/test/mps/test_apply.py index c6296a2..c9d7f83 100644 --- a/test/mps/test_apply.py +++ b/test/mps/test_apply.py @@ -1,16 +1,30 @@ +from __future__ import annotations + import random +import numpy.typing as npt import qailo as q from pytest import approx +from qailo.mps.type import mps +from qailo.util.helpertype import OPSeqElement -def apply(op, m0, m1, m2, m3, v, seq, pos): +def apply( + op: npt.NDArray, + m0: mps, + m1: mps, + m2: mps, + m3: mps, + v: npt.NDArray, + seq: list[OPSeqElement], + pos: list[int], +) -> tuple[mps, mps, mps, mps, npt.NDArray, list[OPSeqElement]]: m0 = q.mps.apply(m0, op, pos) m1 = q.mps.apply(m1, op, pos) m2 = q.mps.apply(m2, op, pos) m3 = q.mps.apply(m3, op, pos) v = q.sv.apply(v, op, pos) - seq.append([op, pos]) + seq.append(OPSeqElement(op, pos)) return m0, m1, m2, m3, v, seq @@ -24,7 +38,7 @@ def test_apply(): m2 = q.mps.zero(n, nkeep=nkeep, mps=q.mps.canonical_mps) m3 = q.mps.zero(n, nkeep=nkeep, mps=q.mps.projector_mps) v = q.sv.zero(n) - seq = [] + seq: list[OPSeqElement] = [] i = 4 j = 0 diff --git a/test/mps/test_canonical.py b/test/mps/test_canonical.py index 98e0e40..fed2611 100644 --- a/test/mps/test_canonical.py +++ b/test/mps/test_canonical.py @@ -1,12 +1,16 @@ +from __future__ import annotations + import numpy as np +import numpy.typing as npt import qailo as q from pytest import approx +from qailo.mps.type import mps as mpstype def test_canonical(): n = 16 nkeep = 4 - tensors = [] + tensors: list[npt.NDArray] = [] d = np.random.randint(2, nkeep) tensors.append(np.random.random((1, 2, d))) for _ in range(n - 2): @@ -16,6 +20,7 @@ def test_canonical(): tensors.append(np.random.random((d, 2, 1))) for mps in [q.mps.canonical_mps, q.mps.projector_mps]: m = mps(tensors) + assert isinstance(m, mpstype) norm = q.norm(m) assert q.norm(m) == approx(norm) assert q.mps.is_canonical(m) @@ -38,6 +43,7 @@ def test_canonical(): tensors = q.mps.tensor_decomposition(v, nkeep) for mps in [q.mps.canonical_mps, q.mps.projector_mps]: m = mps(tensors) + assert isinstance(m, mpstype) norm = q.norm(m) for _ in range(n): diff --git a/test/mps/test_move.py b/test/mps/test_move.py index b420416..893c667 100644 --- a/test/mps/test_move.py +++ b/test/mps/test_move.py @@ -1,14 +1,18 @@ +from __future__ import annotations + import numpy as np +import numpy.typing as npt import qailo as q from pytest import approx from qailo.mps.apply import _move_qubit, _swap_tensors +from qailo.mps.type import mps as mpstype def test_swap(): np.random.seed(1234) n = 12 nkeep = 4 - tensors = [] + tensors: list[npt.NDArray] = [] d = np.random.randint(2, nkeep) tensors.append(np.random.random((1, 2, d))) for _ in range(n - 2): @@ -18,6 +22,7 @@ def test_swap(): tensors.append(np.random.random((d, 2, 1))) for mps in [q.mps.canonical_mps, q.mps.projector_mps]: m = mps(tensors) + assert isinstance(m, mpstype) q.mps.is_canonical(m) norm = q.norm(m) v = q.sv.vector(q.mps.state_vector(m)) @@ -40,7 +45,7 @@ def test_move(): # n = 12 n = 4 nkeep = 4 - tensors = [] + tensors: list[npt.NDArray] = [] d = np.random.randint(2, nkeep) tensors.append(np.random.random((1, 2, d))) for _ in range(n - 2): @@ -50,6 +55,7 @@ def test_move(): tensors.append(np.random.random((d, 2, 1))) for mps in [q.mps.canonical_mps, q.mps.projector_mps]: m = mps(tensors) + assert isinstance(m, mpstype) q.mps.is_canonical(m) norm = q.norm(m) v = q.sv.vector(q.mps.state_vector(m)) diff --git a/test/mps/test_projector.py b/test/mps/test_projector.py index 632f621..ad8b487 100644 --- a/test/mps/test_projector.py +++ b/test/mps/test_projector.py @@ -1,6 +1,6 @@ import numpy as np import qailo as q -from qailo.mps.svd import compact_svd, tensor_svd +from qailo.mps.svd import LegPartition, compact_svd, tensor_svd def test_projector(): @@ -30,7 +30,9 @@ def test_projector(): A = np.einsum(T0, [0, 1, 4, 5], PR, [4, 5, 6], PLh, [7, 8, 6], T1, [8, 7, 2, 3]) L, R = tensor_svd( - np.einsum(T0, [0, 1, 4, 5], T1, [5, 4, 2, 3]), [[0, 1], [2, 3]], nkeep=d + np.einsum(T0, [0, 1, 4, 5], T1, [5, 4, 2, 3]), + LegPartition([0, 1], [2, 3]), + nkeep=d, ) B = np.einsum(L, [0, 1, 4], R, [4, 2, 3]) assert np.allclose(A, B) @@ -42,7 +44,7 @@ def test_projector(): w2 = np.random.random((2, 2)) + 1.0j * np.random.random((2, 2)) p = np.random.random((2, 2, 2, 2)) + 1.0j * np.random.random((2, 2, 2, 2)) B = np.einsum(t0, [0, 4, 6], t1, [6, 5, 3], p, [1, 2, 4, 5]) - p0, p1 = tensor_svd(p, [[0, 2], [1, 3]]) + p0, p1 = tensor_svd(p, LegPartition([0, 2], [1, 3])) assert np.allclose(np.einsum(p0, [0, 2, 4], p1, [4, 1, 3]), p) t0 = np.einsum(t0, [0, 4, 3], p0, [1, 4, 2]) t1 = np.einsum(t1, [0, 4, 3], p1, [1, 2, 4]) @@ -66,7 +68,7 @@ def test_projector(): w2 = np.random.random((2, 2)) + 1.0j * np.random.random((2, 2)) p = np.random.random((2, 2, 2, 2)) + 1.0j * np.random.random((2, 2, 2, 2)) B = np.einsum(t0, [0, 4, 6], t1, [6, 5, 3], p, [1, 2, 4, 5]) - p0, p1 = tensor_svd(p, [[0, 2], [1, 3]]) + p0, p1 = tensor_svd(p, LegPartition([0, 2], [1, 3])) assert np.allclose(np.einsum(p0, [0, 2, 4], p1, [4, 1, 3]), p) t0 = np.einsum(t0, [0, 4, 3], p0, [1, 4, 2]) t1 = np.einsum(t1, [1, 4, 3], p1, [0, 2, 4]) diff --git a/test/mps/test_svd.py b/test/mps/test_svd.py index b4a6137..de99d35 100644 --- a/test/mps/test_svd.py +++ b/test/mps/test_svd.py @@ -1,5 +1,5 @@ import numpy as np -from qailo.mps.svd import compact_svd, tensor_svd +from qailo.mps.svd import LegPartition, compact_svd, tensor_svd def test_compact_svd(): @@ -42,7 +42,7 @@ def test_svd_left(): for _ in range(nt): m, n, p, d = np.random.randint(2, maxn, size=(4,)) T = np.random.random((m, n, p)) - L, R = tensor_svd(T, [[0, 1], [2]], "left") + L, R = tensor_svd(T, LegPartition([0, 1], [2]), "left") assert len(L.shape) == 3 assert len(R.shape) == 2 assert L.shape[0] == m @@ -50,7 +50,7 @@ def test_svd_left(): assert L.shape[2] == R.shape[0] assert R.shape[1] == p assert np.allclose(T, np.einsum("ijl,lk->ijk", L, R)) - L, R = tensor_svd(T, [[0, 1], [2]], "left", nkeep=d) + L, R = tensor_svd(T, LegPartition([0, 1], [2]), "left", nkeep=d) assert len(L.shape) == 3 assert len(R.shape) == 2 assert L.shape[0] == m @@ -65,7 +65,7 @@ def test_svd_left(): for _ in range(nt): m, n, p, d = np.random.randint(2, maxn, size=(4,)) T = np.random.random((m, n, p)) + 1j * np.random.random((m, n, p)) - L, R = tensor_svd(T, [[0, 1], [2]], "left") + L, R = tensor_svd(T, LegPartition([0, 1], [2]), "left") assert len(L.shape) == 3 assert len(R.shape) == 2 assert L.shape[0] == m @@ -73,7 +73,7 @@ def test_svd_left(): assert L.shape[2] == R.shape[0] assert R.shape[1] == p assert np.allclose(T, np.einsum("ijl,lk->ijk", L, R)) - L, R = tensor_svd(T, [[0, 1], [2]], "left", nkeep=d) + L, R = tensor_svd(T, LegPartition([0, 1], [2]), "left", nkeep=d) assert len(L.shape) == 3 assert len(R.shape) == 2 assert L.shape[0] == m @@ -92,7 +92,7 @@ def test_svd_right(): for _ in range(nt): m, n, p, d = np.random.randint(2, maxn, size=(4,)) T = np.random.random((m, n, p)) - L, R = tensor_svd(T, [[0], [1, 2]], "right") + L, R = tensor_svd(T, LegPartition([0], [1, 2]), "right") assert len(L.shape) == 2 assert len(R.shape) == 3 assert L.shape[0] == m @@ -100,7 +100,7 @@ def test_svd_right(): assert R.shape[1] == n assert R.shape[2] == p assert np.allclose(T, np.einsum("il,ljk->ijk", L, R)) - L, R = tensor_svd(T, [[0], [1, 2]], "right", nkeep=d) + L, R = tensor_svd(T, LegPartition([0], [1, 2]), "right", nkeep=d) assert len(L.shape) == 2 assert len(R.shape) == 3 assert L.shape[0] == m @@ -115,7 +115,7 @@ def test_svd_right(): for _ in range(nt): m, n, p, d = np.random.randint(2, maxn, size=(4,)) T = np.random.random((m, n, p)) + 1j * np.random.random((m, n, p)) - L, R = tensor_svd(T, [[0], [1, 2]], "right") + L, R = tensor_svd(T, LegPartition([0], [1, 2]), "right") assert len(L.shape) == 2 assert len(R.shape) == 3 assert L.shape[0] == m @@ -123,7 +123,7 @@ def test_svd_right(): assert R.shape[1] == n assert R.shape[2] == p assert np.allclose(T, np.einsum("il,ljk->ijk", L, R)) - L, R = tensor_svd(T, [[0], [1, 2]], "right", nkeep=d) + L, R = tensor_svd(T, LegPartition([0], [1, 2]), "right", nkeep=d) assert len(L.shape) == 2 assert len(R.shape) == 3 assert L.shape[0] == m @@ -142,7 +142,7 @@ def test_svd_two(): for _ in range(nt): m, n, p, r, d = np.random.randint(2, maxn, size=(5,)) T = np.random.random((m, n, p, r)) - L, R = tensor_svd(T, [[0, 1], [2, 3]]) + L, R = tensor_svd(T, LegPartition([0, 1], [2, 3])) assert len(L.shape) == 3 assert len(R.shape) == 3 assert L.shape[0] == m @@ -151,7 +151,7 @@ def test_svd_two(): assert R.shape[1] == p assert R.shape[2] == r assert np.allclose(T, np.einsum("ijm,mkl->ijkl", L, R)) - L, R = tensor_svd(T, [[0, 1], [2, 3]], nkeep=d) + L, R = tensor_svd(T, LegPartition([0, 1], [2, 3]), nkeep=d) assert len(L.shape) == 3 assert len(R.shape) == 3 assert L.shape[0] == m @@ -167,7 +167,7 @@ def test_svd_two(): for _ in range(nt): m, n, p, r, d = np.random.randint(2, maxn, size=(5,)) T = np.random.random((m, n, p, r)) + 1j * np.random.random((m, n, p, r)) - L, R = tensor_svd(T, [[0, 1], [2, 3]]) + L, R = tensor_svd(T, LegPartition([0, 1], [2, 3])) assert len(L.shape) == 3 assert len(R.shape) == 3 assert L.shape[0] == m @@ -176,7 +176,7 @@ def test_svd_two(): assert R.shape[1] == p assert R.shape[2] == r assert np.allclose(T, np.einsum("ijm,mkl->ijkl", L, R)) - L, R = tensor_svd(T, [[0, 1], [2, 3]], nkeep=d) + L, R = tensor_svd(T, LegPartition([0, 1], [2, 3]), nkeep=d) assert len(L.shape) == 3 assert len(R.shape) == 3 assert L.shape[0] == m diff --git a/test/state_vector/test_apply_seq.py b/test/state_vector/test_apply_seq.py index bd65ae6..cd0aa33 100644 --- a/test/state_vector/test_apply_seq.py +++ b/test/state_vector/test_apply_seq.py @@ -1,15 +1,16 @@ import qailo as q from pytest import approx +from qailo.util.helpertype import OPSeqElement def test_apply_seq(): for n in range(1, 8): v0 = q.sv.zero(n) v1 = q.sv.zero(n) - seq = [] + seq: list[OPSeqElement] = [] for i in range(n): v1 = q.apply(v1, q.op.h(), [i]) - seq.append([q.op.h(), [i]]) + seq.append(OPSeqElement(q.op.h(), [i])) v2 = q.sv.apply_seq(v0, seq) assert q.sv.fidelity(v0, v0) == approx(1) assert q.sv.fidelity(v1, v1) == approx(1) From 6556fa57ba2e50aa91b826fa939cbda4f9e2fcc0 Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 07:11:34 +0900 Subject: [PATCH 35/40] :technologist: Remove np.einsum --- test/mps/test_projector.py | 71 +++++++++++++++++++++++--------------- 1 file changed, 44 insertions(+), 27 deletions(-) diff --git a/test/mps/test_projector.py b/test/mps/test_projector.py index ad8b487..eaf917f 100644 --- a/test/mps/test_projector.py +++ b/test/mps/test_projector.py @@ -1,8 +1,22 @@ +from __future__ import annotations + import numpy as np +import numpy.typing as npt import qailo as q from qailo.mps.svd import LegPartition, compact_svd, tensor_svd +def checked_einsum(*args: npt.NDArray | list[int]) -> npt.NDArray: + for i, arg in enumerate(args): + if i % 2 == 0: + assert isinstance(arg, np.ndarray) + else: + assert isinstance(arg, list) + ret = np.einsum(*args) + assert isinstance(ret, np.ndarray) + return ret + + def test_projector(): maxn = 4 nt = 16 @@ -13,10 +27,11 @@ def test_projector(): _, PLh, PR = q.mps.projector(T0, [0, 2], T1, [2, 1], nkeep=d) assert PLh.shape[1] <= d and PR.shape[1] <= d assert PLh.shape[0] == n1 and PR.shape[0] == n1 - A = np.einsum(T0, [0, 2], PR, [2, 3], PLh.T, [3, 4], T1, [4, 1]) + A = checked_einsum(T0, [0, 2], PR, [2, 3], PLh.T, [3, 4], T1, [4, 1]) - S, U, V = compact_svd(np.einsum(T0, [0, 2], T1, [2, 1]), nkeep=d) + S, U, V = compact_svd(checked_einsum(T0, [0, 2], T1, [2, 1]), nkeep=d) B = np.einsum("k,ik,jk->ij", S, U, V) + assert isinstance(B, np.ndarray) assert np.allclose(A, B) for _ in range(nt): @@ -27,14 +42,16 @@ def test_projector(): assert PLh.shape[2] <= d and PR.shape[2] <= d assert PLh.shape[0] == n2 and PLh.shape[1] == n3 assert PR.shape[0] == n2 and PR.shape[1] == n3 - A = np.einsum(T0, [0, 1, 4, 5], PR, [4, 5, 6], PLh, [7, 8, 6], T1, [8, 7, 2, 3]) + A = checked_einsum( + T0, [0, 1, 4, 5], PR, [4, 5, 6], PLh, [7, 8, 6], T1, [8, 7, 2, 3] + ) L, R = tensor_svd( - np.einsum(T0, [0, 1, 4, 5], T1, [5, 4, 2, 3]), + checked_einsum(T0, [0, 1, 4, 5], T1, [5, 4, 2, 3]), LegPartition([0, 1], [2, 3]), nkeep=d, ) - B = np.einsum(L, [0, 1, 4], R, [4, 2, 3]) + B = checked_einsum(L, [0, 1, 4], R, [4, 2, 3]) assert np.allclose(A, B) for _ in range(nt): @@ -43,21 +60,21 @@ def test_projector(): t1 = np.random.random((3, 2, 2)) + 1.0j * np.random.random((3, 2, 2)) w2 = np.random.random((2, 2)) + 1.0j * np.random.random((2, 2)) p = np.random.random((2, 2, 2, 2)) + 1.0j * np.random.random((2, 2, 2, 2)) - B = np.einsum(t0, [0, 4, 6], t1, [6, 5, 3], p, [1, 2, 4, 5]) + B = checked_einsum(t0, [0, 4, 6], t1, [6, 5, 3], p, [1, 2, 4, 5]) p0, p1 = tensor_svd(p, LegPartition([0, 2], [1, 3])) - assert np.allclose(np.einsum(p0, [0, 2, 4], p1, [4, 1, 3]), p) - t0 = np.einsum(t0, [0, 4, 3], p0, [1, 4, 2]) - t1 = np.einsum(t1, [0, 4, 3], p1, [1, 2, 4]) - assert np.allclose(np.einsum(t0, [0, 1, 4, 5], t1, [5, 4, 2, 3]), B) - tt0 = np.einsum(w0, [0, 4], t0, [4, 1, 2, 3]) - tt1 = np.einsum(t1, [0, 1, 2, 4], w2, [4, 3]) + assert np.allclose(checked_einsum(p0, [0, 2, 4], p1, [4, 1, 3]), p) + t0 = checked_einsum(t0, [0, 4, 3], p0, [1, 4, 2]) + t1 = checked_einsum(t1, [0, 4, 3], p1, [1, 2, 4]) + assert np.allclose(checked_einsum(t0, [0, 1, 4, 5], t1, [5, 4, 2, 3]), B) + tt0 = checked_einsum(w0, [0, 4], t0, [4, 1, 2, 3]) + tt1 = checked_einsum(t1, [0, 1, 2, 4], w2, [4, 3]) _, PLh, PR = q.mps.projector(tt0, [0, 1, 4, 5], tt1, [5, 4, 2, 3]) assert np.allclose( - np.einsum(PLh, [2, 3, 0], PR, [2, 3, 1]), np.identity(PLh.shape[2]) + checked_einsum(PLh, [2, 3, 0], PR, [2, 3, 1]), np.identity(PLh.shape[2]) ) - tt0 = np.einsum(t0, [0, 1, 3, 4], PR, [3, 4, 2]) - tt1 = np.einsum(PLh, [3, 4, 0], t1, [4, 3, 1, 2]) - A = np.einsum(tt0, [0, 1, 4], tt1, [4, 2, 3]) + tt0 = checked_einsum(t0, [0, 1, 3, 4], PR, [3, 4, 2]) + tt1 = checked_einsum(PLh, [3, 4, 0], t1, [4, 3, 1, 2]) + A = checked_einsum(tt0, [0, 1, 4], tt1, [4, 2, 3]) print(np.linalg.norm(A - B)) assert np.allclose(A, B) @@ -67,21 +84,21 @@ def test_projector(): t1 = np.random.random((3, 2, 2)) + 1.0j * np.random.random((3, 2, 2)) w2 = np.random.random((2, 2)) + 1.0j * np.random.random((2, 2)) p = np.random.random((2, 2, 2, 2)) + 1.0j * np.random.random((2, 2, 2, 2)) - B = np.einsum(t0, [0, 4, 6], t1, [6, 5, 3], p, [1, 2, 4, 5]) + B = checked_einsum(t0, [0, 4, 6], t1, [6, 5, 3], p, [1, 2, 4, 5]) p0, p1 = tensor_svd(p, LegPartition([0, 2], [1, 3])) - assert np.allclose(np.einsum(p0, [0, 2, 4], p1, [4, 1, 3]), p) - t0 = np.einsum(t0, [0, 4, 3], p0, [1, 4, 2]) - t1 = np.einsum(t1, [1, 4, 3], p1, [0, 2, 4]) - assert np.allclose(np.einsum(t0, [0, 1, 4, 5], t1, [4, 5, 2, 3]), B) - tt0 = np.einsum(w0, [0, 4], t0, [4, 1, 2, 3]) - tt1 = np.einsum(t1, [0, 1, 2, 4], w2, [4, 3]) + assert np.allclose(checked_einsum(p0, [0, 2, 4], p1, [4, 1, 3]), p) + t0 = checked_einsum(t0, [0, 4, 3], p0, [1, 4, 2]) + t1 = checked_einsum(t1, [1, 4, 3], p1, [0, 2, 4]) + assert np.allclose(checked_einsum(t0, [0, 1, 4, 5], t1, [4, 5, 2, 3]), B) + tt0 = checked_einsum(w0, [0, 4], t0, [4, 1, 2, 3]) + tt1 = checked_einsum(t1, [0, 1, 2, 4], w2, [4, 3]) _, PLh, PR = q.mps.projector(tt0, [0, 1, 4, 5], tt1, [4, 5, 2, 3]) assert np.allclose( - np.einsum(PLh, [2, 3, 0], PR, [2, 3, 1]), np.identity(PLh.shape[2]) + checked_einsum(PLh, [2, 3, 0], PR, [2, 3, 1]), np.identity(PLh.shape[2]) ) - tt0 = np.einsum(t0, [0, 1, 3, 4], PR, [3, 4, 2]) - tt1 = np.einsum(PLh, [3, 4, 0], t1, [3, 4, 1, 2]) - A = np.einsum(tt0, [0, 1, 4], tt1, [4, 2, 3]) + tt0 = checked_einsum(t0, [0, 1, 3, 4], PR, [3, 4, 2]) + tt1 = checked_einsum(PLh, [3, 4, 0], t1, [3, 4, 1, 2]) + A = checked_einsum(tt0, [0, 1, 4], tt1, [4, 2, 3]) print(np.linalg.norm(A - B)) assert np.allclose(A, B) From c24215256b01a99cecd15160d42160d89e4d1216 Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 07:32:54 +0900 Subject: [PATCH 36/40] :label: Add py.typed --- pyproject.toml | 11 +++++------ src/qailo/py.typed | 0 2 files changed, 5 insertions(+), 6 deletions(-) create mode 100644 src/qailo/py.typed diff --git a/pyproject.toml b/pyproject.toml index 5a175b8..f739bdd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -6,13 +6,9 @@ build-backend = "setuptools.build_meta" name = "qailo" description = "Simplest Quantum Circuit Simulator" readme = "README.md" -authors = [ - { name = "Synge Todo", email = "wistaria@phys.s.u-tokyo.ac.jp" }, -] +authors = [{ name = "Synge Todo", email = "wistaria@phys.s.u-tokyo.ac.jp" }] license = { file = "LICENSE" } -classifiers = [ - "Programming Language :: Python :: 3", -] +classifiers = ["Programming Language :: Python :: 3"] dynamic = ["version"] dependencies = ["matplotlib", "numpy"] @@ -25,6 +21,9 @@ write_to = "src/qailo/_version.py" [tool.setuptools] package-dir = { "" = "src" } +[tool.setuptools.package-data] +"qailo" = ["py.typed"] + [project.optional-dependencies] dev = ["pytest", "black", "ruff", "typing-extensions"] diff --git a/src/qailo/py.typed b/src/qailo/py.typed new file mode 100644 index 0000000..e69de29 From b8eacf9e650fe86e228af3111f5229590235d57c Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 07:42:51 +0900 Subject: [PATCH 37/40] :bug: Add __future__ --- test/state_vector/test_apply_seq.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/state_vector/test_apply_seq.py b/test/state_vector/test_apply_seq.py index cd0aa33..9cb0d7e 100644 --- a/test/state_vector/test_apply_seq.py +++ b/test/state_vector/test_apply_seq.py @@ -1,3 +1,5 @@ +from __future__ import annotations + import qailo as q from pytest import approx from qailo.util.helpertype import OPSeqElement From 5debf738e787fc2175b08f87ebd156a5a0aeb6d4 Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 08:00:27 +0900 Subject: [PATCH 38/40] :heavy_minus_sign: Import from typing --- src/qailo/mps/svd.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/qailo/mps/svd.py b/src/qailo/mps/svd.py index c85cfc3..a8a3c80 100644 --- a/src/qailo/mps/svd.py +++ b/src/qailo/mps/svd.py @@ -1,10 +1,9 @@ from __future__ import annotations -from typing import NamedTuple +from typing import Literal, NamedTuple import numpy as np import numpy.typing as npt -from typing_extensions import Literal from ..typeutil import eincheck as ec From 66214ec3e40db7bd843346dc750bc80e7285e5b5 Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 08:00:54 +0900 Subject: [PATCH 39/40] :heavy_plus_sign: Move typing-extensions to deps --- pyproject.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index f739bdd..cd666a4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -10,7 +10,7 @@ authors = [{ name = "Synge Todo", email = "wistaria@phys.s.u-tokyo.ac.jp" }] license = { file = "LICENSE" } classifiers = ["Programming Language :: Python :: 3"] dynamic = ["version"] -dependencies = ["matplotlib", "numpy"] +dependencies = ["matplotlib", "numpy", "typing-extensions"] [tool.setuptools.dynamic] version = { attr = "qailo._version.version" } @@ -25,7 +25,7 @@ package-dir = { "" = "src" } "qailo" = ["py.typed"] [project.optional-dependencies] -dev = ["pytest", "black", "ruff", "typing-extensions"] +dev = ["pytest", "black", "ruff"] [tool.black] target-version = ["py38"] From f7538ae6acff95213cbcd6d8f7a2662e5e7132df Mon Sep 17 00:00:00 2001 From: SS Date: Mon, 4 Mar 2024 08:15:39 +0900 Subject: [PATCH 40/40] :bug: Add mps compatibility --- src/qailo/alg/qpe.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/qailo/alg/qpe.py b/src/qailo/alg/qpe.py index 6f880fa..1147ad2 100644 --- a/src/qailo/alg/qpe.py +++ b/src/qailo/alg/qpe.py @@ -1,14 +1,24 @@ from __future__ import annotations from copy import deepcopy +from typing import overload import numpy as np import numpy.typing as npt import qailo as q +from qailo.mps.type import mps -def qpe(n: int, u: npt.NDArray, v: npt.NDArray) -> npt.NDArray: +@overload +def qpe(n: int, u: npt.NDArray, v: npt.NDArray) -> npt.NDArray: ... + + +@overload +def qpe(n: int, u: npt.NDArray, v: mps) -> mps: ... + + +def qpe(n: int, u: npt.NDArray, v: npt.NDArray | mps) -> npt.NDArray | mps: m = q.num_qubits(u) w = deepcopy(v) assert q.num_qubits(w) == m + n