feat: experimental exact simulators (3/5) (#1855)

docs/source/index.rst

@@ -32,6 +32,7 @@ If you’re new to pyQuil, head to the `getting started <getting_started>`_ guid
    wavefunction_simulator
    compiler
    noise
+   simulation_architecture
    advanced_usage
    troubleshooting
    introducing_v4

pyproject.toml

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "pyquil"
-version = "4.21.0"
+version = "4.22.0"
 description = "A Python library for creating Quantum Instruction Language (Quil) programs."
 authors = ["Rigetti Computing <softapps@rigetti.com>"]
 readme = "README.md"

test/unit/test_noise_model.py

@@ -33,6 +33,16 @@
 from pyquil.quil import Program
 from pyquil.quilatom import FormalArgument, Qubit
 from pyquil.quilbase import DefCircuit, Gate, Measurement, ResetQubit
+from pyquil.simulation._simulator import DensityMatrixSimulator
+
+_EMPTY_PARAMS = jnp.array([], dtype=float)
+
+
+def _dm(program, noise_model=None, qubits=None):
+    """Compute density matrix."""
+    sim = DensityMatrixSimulator(program, qubits=qubits, noise_model=noise_model)
+    return sim.compute(_EMPTY_PARAMS)
+
 
 # ──────────────────────────────────────────────────────────
 # Channel tests
@@ -491,6 +501,59 @@ def test_qubits(self):
         ch = ResetChannel.from_reset_fidelity(inst=inst, fidelity=0.99)
         assert ch.qubits == [3]
 
+    def test_ideal_reset_maps_excited_to_ground(self):
+        """An ideal reset on an excited qubit should produce |0><0|."""
+        inst = RESET(0)
+        ch = ResetChannel.from_reset_fidelity(inst=inst, fidelity=1.0)
+        noise_model = NoiseModel.from_channels([ch])
+        # Prepare |1> then reset
+        program = Program(X(0), RESET(0))
+        rho = _dm(program, noise_model=noise_model)
+        target_rho = qx.zero_state_matrix(1)
+        assert qx.fidelity(rho, target_rho) > 0.9999
+
+    def test_ideal_reset_maps_superposition_to_ground(self):
+        """An ideal reset on a superposition state should produce |0><0|."""
+        inst = RESET(0)
+        ch = ResetChannel.from_reset_fidelity(inst=inst, fidelity=1.0)
+        noise_model = NoiseModel.from_channels([ch])
+        # Prepare |+> then reset
+        program = Program(RX(np.pi / 2, 0), RESET(0))
+        rho = _dm(program, noise_model=noise_model)
+        target_rho = qx.zero_state_matrix(1)
+        assert qx.fidelity(rho, target_rho) > 0.9999
+
+    def test_noisy_reset_reduces_fidelity(self):
+        """A noisy reset should produce a state with fidelity < 1 relative to |0><0|."""
+        inst = RESET(0)
+        ch = ResetChannel.from_reset_fidelity(inst=inst, fidelity=0.90)
+        noise_model = NoiseModel.from_channels([ch])
+        program = Program(X(0), RESET(0))
+        rho = _dm(program, noise_model=noise_model)
+        target_rho = qx.zero_state_matrix(1)
+        fid = float(qx.fidelity(rho, target_rho))
+        # Should be less than perfect but still high
+        assert 0.85 < fid < 1.0
+
+    def test_reset_in_multi_qubit_circuit(self):
+        """Reset on one qubit should not affect the other qubit."""
+        inst = RESET(0)
+        ch = ResetChannel.from_reset_fidelity(inst=inst, fidelity=1.0)
+        noise_model = NoiseModel.from_channels([ch])
+        # Prepare |11> then reset qubit 0
+        program = Program(X(0), X(1), RESET(0))
+        rho = _dm(program, noise_model=noise_model)
+        # Expected state: |0> on qubit 0, |1> on qubit 1 → |01>
+        target_rho = (qx.gates.I | qx.gates.X) @ qx.zero_state_matrix(2)
+        assert qx.fidelity(rho, target_rho) > 0.9999
+
+    def test_global_reset(self):
+        """A global RESET (no qubit specified) resets all qubits to |0>."""
+        program = Program(X(0), X(1), RESET())
+        rho = _dm(program)
+        target_rho = qx.zero_state_matrix(2)
+        assert qx.fidelity(rho, target_rho) > 0.9999
+
     def test_global_reset_channel_rejected(self):
         """ResetChannel is intentionally scoped to targeted resets."""
         with pytest.raises(TypeError, match="targeted"):