Added state fidelity

CHANGELOG.md

@@ -10,6 +10,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added
 
 ### Fixed
+- #435
+  - Added `graphix.sim.statevec.Statevec.fidelity` and `graphix.sim.statevec.Statevec.isclose` methods to calculate and check fidelity between two states.
+  - Updated old code with uses of the above methods.
+  - Solves [#386](https://github.com/TeamGraphix/graphix/issues/386).
 
 - #429
   - Modify `graphix.noise_models.noise_model.ApplyNoise` to handle conditionality based on a `domain` attribute (like `command.X` and `command.Z`).

graphix/sim/statevec.py

@@ -400,6 +400,47 @@ def xreplace(self, assignment: Mapping[Parameter, ExpressionOrSupportsFloat]) ->
         result.psi = np.vectorize(lambda value: parameter.xreplace(value, assignment))(self.psi)
         return result
 
+    def fidelity(self, other: Statevec) -> float:
+        r"""Return the fidelity with another statevector.
+
+        This is calculated as
+
+        .. math::
+            F(\psi_{1}, \psi_{2}) = |\langle \psi_{1} | \psi_{2} \rangle|^2
+
+        Parameters
+        ----------
+        other : :class:`graphix.sim.statevec.Statevec`
+            other statevector
+
+        Returns
+        -------
+        float : fidelity between self and other
+        """
+        st1 = copy.copy(self)
+        st1.normalize()
+        st2 = copy.copy(other)
+        st2.normalize()
+        return float(np.abs(np.vdot(st1.psi.flatten(), st2.psi.flatten()))) ** 2  # type: ignore[arg-type]
+
+    def isclose(self, other: Statevec, atol: float = 1e-15, rtol: float = 1e-10) -> bool:
+        r"""Return whether the state is equal to another state up to global phase.
+
+        Parameters
+        ----------
+        other : :class:`graphix.sim.statevec.Statevec`
+            other statevector
+        atol : float, optional
+            absolute tolerance for numerical comparison, defaults to 1e-15
+        rtol : float, optional
+            relative tolerance for numerical comparison, defaults to 1e-10
+
+        Returns
+        -------
+        bool : whether the states are equal up to global phase
+        """
+        return math.isclose(self.fidelity(other), 1, abs_tol=atol, rel_tol=rtol)
+
 
 @dataclass(frozen=True)
 class StatevectorBackend(DenseStateBackend[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
@@ -799,9 +798,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):
+            if state.isclose(state_ref):
                 continue
             return False
 
@@ -871,7 +869,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
@@ -884,7 +882,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
 
@@ -49,7 +48,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_ref.isclose(state_p)
 
 
 @pytest.mark.parametrize("jumps", range(1, 11))
@@ -66,7 +65,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))
@@ -100,7 +99,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:
@@ -122,4 +121,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

@@ -168,7 +168,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.isclose(state_mbqc)
 
 
 def test_visualization() -> None: