Add fidelity and isclose methods to Statevec

CHANGELOG.md

@@ -9,6 +9,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+- #386, #433: Added `Statevec.fidelity` and `Statevec.isclose` methods for pure-state fidelity computation and equality check up to global phase.
+
 ### Fixed
 
 - #429

graphix/sim/statevec.py

@@ -388,6 +388,45 @@ def expectation_value(self, op: Matrix, qargs: Sequence[int]) -> complex:
         st1.evolve(op, qargs)
         return complex(np.dot(st2.psi.flatten().conjugate(), st1.psi.flatten()))
 
+    def fidelity(self, other: Statevec) -> float:
+        r"""Calculate the fidelity against another statevector.
+
+        The fidelity is defined as :math:`|\langle\psi_1|\psi_2\rangle|^2`.
+
+        Parameters
+        ----------
+        other : :class:`Statevec`
+            statevector to compare with
+
+        Returns
+        -------
+        float
+            Fidelity between the two statevectors.
+        """
+        inner = np.dot(self.flatten().conjugate(), other.flatten())
+        return float(np.abs(inner) ** 2)
+
+    def isclose(self, other: Statevec, *, rtol: float = 1e-09, atol: float = 0.0) -> bool:
+        """Check if two quantum states are equal up to global phase.
+
+        Two states are considered close if their fidelity is close to 1.
+
+        Parameters
+        ----------
+        other : :class:`Statevec`
+            statevector to compare with
+        rtol : float
+            relative tolerance for :func:`math.isclose`
+        atol : float
+            absolute tolerance for :func:`math.isclose`
+
+        Returns
+        -------
+        bool
+            ``True`` if the states are equal up to global phase.
+        """
+        return math.isclose(self.fidelity(other), 1, rel_tol=rtol, abs_tol=atol)
+
     def subs(self, variable: Parameter, substitute: ExpressionOrSupportsFloat) -> Statevec:
         """Return a copy of the state vector where all occurrences of the given variable in measurement angles are substituted by the given value."""
         result = Statevec()

tests/test_flow_core.py

@@ -4,7 +4,6 @@
 from typing import TYPE_CHECKING, NamedTuple
 
 import networkx as nx
-import numpy as np
 import pytest
 
 from graphix.command import E, M, N, X, Z
@@ -412,8 +411,7 @@ def test_corrections_to_pattern(self, test_case: XZCorrectionsTestCase, fx_rng:
 
                 for _ in range(n_shots):
                     state = pattern.simulate_pattern(input_state=PlanarState(plane, alpha))
-                    result = np.abs(np.dot(state.flatten().conjugate(), state_ref.flatten()))
-                    assert result == pytest.approx(1)
+                    assert state.isclose(state_ref)
 
 
 class TestFlow:

tests/test_graphsim.py

@@ -5,7 +5,6 @@
 import networkx as nx
 import numpy as np
 import numpy.typing as npt
-import pytest
 
 from graphix.clifford import Clifford
 from graphix.fundamentals import ANGLE_PI, Plane, angle_to_rad
@@ -87,21 +86,21 @@ def test_fig2(self) -> None:
         gstate.normalize()
         gstate.remove_qubit(0)
         gstate2 = graph_state_to_statevec(g)
-        assert np.abs(np.dot(gstate.flatten().conjugate(), gstate2.flatten())) == pytest.approx(1)
+        assert gstate.isclose(gstate2)
 
         g.measure_y(1, choice=0)
         gstate.evolve_single(meas_op(0.5 * ANGLE_PI), 0)  # y meas
         gstate.normalize()
         gstate.remove_qubit(0)
         gstate2 = graph_state_to_statevec(g)
-        assert np.abs(np.dot(gstate.flatten().conjugate(), gstate2.flatten())) == pytest.approx(1)
+        assert gstate.isclose(gstate2)
 
         g.measure_z(3)
         gstate.evolve_single(meas_op(0.5 * ANGLE_PI, plane=Plane.YZ), 1)  # z meas
         gstate.normalize()
         gstate.remove_qubit(1)
         gstate2 = graph_state_to_statevec(g)
-        assert np.abs(np.dot(gstate.flatten().conjugate(), gstate2.flatten())) == pytest.approx(1)
+        assert gstate.isclose(gstate2)
 
     def test_e2(self) -> None:
         nqubit = 6
@@ -112,23 +111,23 @@ def test_e2(self) -> None:
 
         g.equivalent_graph_e2(3, 4)
         gstate2 = graph_state_to_statevec(g)
-        assert np.abs(np.dot(gstate.flatten().conjugate(), gstate2.flatten())) == pytest.approx(1)
+        assert gstate.isclose(gstate2)
 
         g.equivalent_graph_e2(4, 0)
         gstate3 = graph_state_to_statevec(g)
-        assert np.abs(np.dot(gstate.flatten().conjugate(), gstate3.flatten())) == pytest.approx(1)
+        assert gstate.isclose(gstate3)
 
         g.equivalent_graph_e2(4, 5)
         gstate4 = graph_state_to_statevec(g)
-        assert np.abs(np.dot(gstate.flatten().conjugate(), gstate4.flatten())) == pytest.approx(1)
+        assert gstate.isclose(gstate4)
 
         g.equivalent_graph_e2(0, 3)
         gstate5 = graph_state_to_statevec(g)
-        assert np.abs(np.dot(gstate.flatten().conjugate(), gstate5.flatten())) == pytest.approx(1)
+        assert gstate.isclose(gstate5)
 
         g.equivalent_graph_e2(0, 3)
         gstate6 = graph_state_to_statevec(g)
-        assert np.abs(np.dot(gstate.flatten().conjugate(), gstate6.flatten())) == pytest.approx(1)
+        assert gstate.isclose(gstate6)
 
     def test_e1(self) -> None:
         nqubit = 6
@@ -139,15 +138,15 @@ def test_e1(self) -> None:
         g.equivalent_graph_e1(3)
 
         gstate2 = graph_state_to_statevec(g)
-        assert np.abs(np.dot(gstate.flatten().conjugate(), gstate2.flatten())) == pytest.approx(1)
+        assert gstate.isclose(gstate2)
         g.z(4)
         gstate = graph_state_to_statevec(g)
         g.equivalent_graph_e1(4)
         gstate2 = graph_state_to_statevec(g)
-        assert np.abs(np.dot(gstate.flatten().conjugate(), gstate2.flatten())) == pytest.approx(1)
+        assert gstate.isclose(gstate2)
         g.equivalent_graph_e1(4)
         gstate3 = graph_state_to_statevec(g)
-        assert np.abs(np.dot(gstate.flatten().conjugate(), gstate3.flatten())) == pytest.approx(1)
+        assert gstate.isclose(gstate3)
 
     def test_local_complement(self) -> None:
         nqubit = 6

tests/test_opengraph.py

@@ -10,7 +10,6 @@
 from typing import TYPE_CHECKING, NamedTuple
 
 import networkx as nx
-import numpy as np
 import pytest
 
 from graphix.command import E
@@ -800,11 +799,8 @@ def check_determinism(pattern: Pattern, fx_rng: Generator, n_shots: int = 3) ->
 
         for _ in range(n_shots):
             state = pattern.simulate_pattern(input_state=PlanarState(plane, alpha))
-            result = np.abs(np.dot(state.flatten().conjugate(), state_ref.flatten()))
-
-            if result == pytest.approx(1):
-                continue
-            return False
+            if not state.isclose(state_ref):
+                return False
 
     return True
 
@@ -872,7 +868,7 @@ def test_double_entanglement(self) -> None:
         pattern2 = pattern.extract_opengraph().to_pattern()
         state = pattern.simulate_pattern()
         state2 = pattern2.simulate_pattern()
-        assert np.abs(np.dot(state.flatten().conjugate(), state2.flatten())) == pytest.approx(1)
+        assert state.isclose(state2)
 
     def test_from_to_pattern(self, fx_rng: Generator) -> None:
         n_qubits = 2
@@ -885,7 +881,7 @@ def test_from_to_pattern(self, fx_rng: Generator) -> None:
             alpha = 2 * ANGLE_PI * fx_rng.random()
             state_ref = pattern_ref.simulate_pattern(input_state=PlanarState(plane, alpha))
             state = pattern.simulate_pattern(input_state=PlanarState(plane, alpha))
-            assert np.abs(np.dot(state.flatten().conjugate(), state_ref.flatten())) == pytest.approx(1)
+            assert state.isclose(state_ref)
 
     def test_isclose_measurement(self) -> None:
         og_1 = OpenGraph(

tests/test_optimization.py

@@ -1,6 +1,5 @@
 from __future__ import annotations
 
-import numpy as np
 import pytest
 from numpy.random import PCG64, Generator
 
@@ -50,7 +49,7 @@ def test_standardize_clifford_entanglement(fx_rng: Generator) -> None:
 
                 state_ref = p_ref.simulate_pattern(input_state=PlanarState(Plane.XY, alpha))
                 state_p = p.simulate_pattern(input_state=PlanarState(Plane.XY, alpha))
-                assert np.abs(np.dot(state_p.flatten().conjugate(), state_ref.flatten())) == pytest.approx(1)
+                assert state_p.isclose(state_ref)
 
 
 @pytest.mark.parametrize("jumps", range(1, 11))
@@ -67,7 +66,7 @@ def test_incorporate_pauli_results(fx_bg: PCG64, jumps: int) -> None:
     pattern2 = incorporate_pauli_results(pattern)
     state = pattern.simulate_pattern(rng=rng)
     state2 = pattern2.simulate_pattern(rng=rng)
-    assert np.abs(np.dot(state.flatten().conjugate(), state2.flatten())) == pytest.approx(1)
+    assert state.isclose(state2)
 
 
 @pytest.mark.parametrize("jumps", range(1, 11))
@@ -101,7 +100,7 @@ def test_remove_useless_domains(fx_bg: PCG64, jumps: int) -> None:
     pattern2 = remove_useless_domains(pattern)
     state = pattern.simulate_pattern(rng=rng)
     state2 = pattern2.simulate_pattern(rng=rng)
-    assert np.abs(np.dot(state.flatten().conjugate(), state2.flatten())) == pytest.approx(1)
+    assert state.isclose(state2)
 
 
 def test_to_space_optimal_pattern() -> None:
@@ -123,4 +122,4 @@ def test_to_space_optimal_pattern() -> None:
     pattern2 = StandardizedPattern.from_pattern(pattern).to_space_optimal_pattern()
     state = pattern.simulate_pattern()
     state2 = pattern2.simulate_pattern()
-    assert np.abs(np.dot(state.flatten().conjugate(), state2.flatten())) == pytest.approx(1)
+    assert state.isclose(state2)

tests/test_parameter.py

@@ -180,7 +180,7 @@ def test_random_circuit_with_parameters(fx_bg: PCG64, jumps: int, use_xreplace:
     else:
         state = circuit.subs(alpha, assignment[alpha]).subs(beta, assignment[beta]).simulate_statevector().statevec
         state_mbqc = pattern.subs(alpha, assignment[alpha]).subs(beta, assignment[beta]).simulate_pattern()
-    assert np.abs(np.dot(state_mbqc.flatten().conjugate(), state.flatten())) == pytest.approx(1)
+    assert state_mbqc.isclose(state)
 
 
 def test_visualization() -> None: