From b63b59fad6baa7d66284b0d2ffdd6d73cb6e5ca6 Mon Sep 17 00:00:00 2001
From: Orion Martin <40585662+95-martin-orion@users.noreply.github.com>
Date: Tue, 3 May 2022 11:27:28 -0700
Subject: [PATCH 1/3] Set up for no-copy state_vector.

---
 .../cirq/circuits/circuit_operation_test.py   | 34 +++++++++----------
 .../contrib/quantum_volume/quantum_volume.py  |  4 +--
 .../grid_parallel_two_qubit_xeb.py            |  3 +-
 cirq-core/cirq/experiments/xeb_simulation.py  |  3 +-
 .../cirq/ops/boolean_hamiltonian_test.py      |  2 +-
 cirq-core/cirq/ops/common_gates_test.py       | 32 ++++++++++++-----
 cirq-core/cirq/sim/mux.py                     |  2 +-
 cirq-core/cirq/sim/sparse_simulator.py        | 10 ++++--
 cirq-core/cirq/sim/sparse_simulator_test.py   | 28 +++++++--------
 cirq-core/cirq/sim/state_vector.py            | 26 +++++++++-----
 cirq-core/cirq/sim/state_vector_simulator.py  | 12 +++++--
 .../cirq/sim/state_vector_simulator_test.py   |  2 +-
 cirq-core/cirq/sim/state_vector_test.py       |  5 +--
 .../calibration/engine_simulator.py           |  2 +-
 docs/simulation.ipynb                         |  2 +-
 docs/tutorials/educators/intro.ipynb          |  2 +-
 docs/tutorials/fourier_checking.ipynb         |  2 +-
 examples/stabilizer_code.ipynb                |  8 ++---
 examples/stabilizer_code_test.py              |  2 +-
 19 files changed, 109 insertions(+), 72 deletions(-)

diff --git a/cirq-core/cirq/circuits/circuit_operation_test.py b/cirq-core/cirq/circuits/circuit_operation_test.py
index 2820bd1a654..2c7673c3358 100644
--- a/cirq-core/cirq/circuits/circuit_operation_test.py
+++ b/cirq-core/cirq/circuits/circuit_operation_test.py
@@ -265,7 +265,7 @@ def test_recursive_params():
 
     # First example should behave like an X when simulated
     result = cirq.Simulator().simulate(cirq.Circuit(circuitop), param_resolver=outer_params)
-    assert np.allclose(result.state_vector(), [0, 1])
+    assert np.allclose(result.state_vector(copy=False), [0, 1])
 
 
 @pytest.mark.parametrize('add_measurements', [True, False])
@@ -341,9 +341,9 @@ def test_repeat_zero_times(add_measurements, use_repetition_ids, initial_reps):
         subcircuit.freeze(), repetitions=initial_reps, use_repetition_ids=use_repetition_ids
     )
     result = cirq.Simulator().simulate(cirq.Circuit(op))
-    assert np.allclose(result.state_vector(), [0, 1] if initial_reps % 2 else [1, 0])
+    assert np.allclose(result.state_vector(copy=False), [0, 1] if initial_reps % 2 else [1, 0])
     result = cirq.Simulator().simulate(cirq.Circuit(op**0))
-    assert np.allclose(result.state_vector(), [1, 0])
+    assert np.allclose(result.state_vector(copy=False), [1, 0])
 
 
 def test_parameterized_repeat():
@@ -352,13 +352,13 @@ def test_parameterized_repeat():
     assert cirq.parameter_names(op) == {'a'}
     assert not cirq.has_unitary(op)
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 0})
-    assert np.allclose(result.state_vector(), [1, 0])
+    assert np.allclose(result.state_vector(copy=False), [1, 0])
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 1})
-    assert np.allclose(result.state_vector(), [0, 1])
+    assert np.allclose(result.state_vector(copy=False), [0, 1])
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 2})
-    assert np.allclose(result.state_vector(), [1, 0])
+    assert np.allclose(result.state_vector(copy=False), [1, 0])
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': -1})
-    assert np.allclose(result.state_vector(), [0, 1])
+    assert np.allclose(result.state_vector(copy=False), [0, 1])
     with pytest.raises(TypeError, match='Only integer or sympy repetitions are allowed'):
         cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 1.5})
     with pytest.raises(ValueError, match='Circuit contains ops whose symbols were not specified'):
@@ -367,13 +367,13 @@ def test_parameterized_repeat():
     assert cirq.parameter_names(op) == {'a'}
     assert not cirq.has_unitary(op)
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 0})
-    assert np.allclose(result.state_vector(), [1, 0])
+    assert np.allclose(result.state_vector(copy=False), [1, 0])
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 1})
-    assert np.allclose(result.state_vector(), [0, 1])
+    assert np.allclose(result.state_vector(copy=False), [0, 1])
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 2})
-    assert np.allclose(result.state_vector(), [1, 0])
+    assert np.allclose(result.state_vector(copy=False), [1, 0])
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': -1})
-    assert np.allclose(result.state_vector(), [0, 1])
+    assert np.allclose(result.state_vector(copy=False), [0, 1])
     with pytest.raises(TypeError, match='Only integer or sympy repetitions are allowed'):
         cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 1.5})
     with pytest.raises(ValueError, match='Circuit contains ops whose symbols were not specified'):
@@ -382,11 +382,11 @@ def test_parameterized_repeat():
     assert cirq.parameter_names(op) == {'a', 'b'}
     assert not cirq.has_unitary(op)
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 1, 'b': 1})
-    assert np.allclose(result.state_vector(), [0, 1])
+    assert np.allclose(result.state_vector(copy=False), [0, 1])
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 2, 'b': 1})
-    assert np.allclose(result.state_vector(), [1, 0])
+    assert np.allclose(result.state_vector(copy=False), [1, 0])
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 1, 'b': 2})
-    assert np.allclose(result.state_vector(), [1, 0])
+    assert np.allclose(result.state_vector(copy=False), [1, 0])
     with pytest.raises(TypeError, match='Only integer or sympy repetitions are allowed'):
         cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 1.5, 'b': 1})
     with pytest.raises(ValueError, match='Circuit contains ops whose symbols were not specified'):
@@ -395,11 +395,11 @@ def test_parameterized_repeat():
     assert cirq.parameter_names(op) == {'a', 'b'}
     assert not cirq.has_unitary(op)
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 1, 'b': 1})
-    assert np.allclose(result.state_vector(), [1, 0])
+    assert np.allclose(result.state_vector(copy=False), [1, 0])
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 1.5, 'b': 1})
-    assert np.allclose(result.state_vector(), [0, 1])
+    assert np.allclose(result.state_vector(copy=False), [0, 1])
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 1, 'b': 1.5})
-    assert np.allclose(result.state_vector(), [0, 1])
+    assert np.allclose(result.state_vector(copy=False), [0, 1])
     with pytest.raises(TypeError, match='Only integer or sympy repetitions are allowed'):
         cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 1.5, 'b': 1.5})
     with pytest.raises(ValueError, match='Circuit contains ops whose symbols were not specified'):
diff --git a/cirq-core/cirq/contrib/quantum_volume/quantum_volume.py b/cirq-core/cirq/contrib/quantum_volume/quantum_volume.py
index 042a895d4d4..0c5ca5a1fa8 100644
--- a/cirq-core/cirq/contrib/quantum_volume/quantum_volume.py
+++ b/cirq-core/cirq/contrib/quantum_volume/quantum_volume.py
@@ -81,12 +81,12 @@ def compute_heavy_set(circuit: cirq.Circuit) -> List[int]:
     # output is defined in terms of probabilities, where our wave function is in
     # terms of amplitudes. We convert it by using the Born rule: squaring each
     # amplitude and taking their absolute value
-    median = np.median(np.abs(results.state_vector() ** 2))
+    median = np.median(np.abs(results.state_vector(copy=False) ** 2))
 
     # The output wave function is a vector from the result value (big-endian) to
     # the probability of that bit-string. Return all of the bit-string
     # values that have a probability greater than the median.
-    return [idx for idx, amp in enumerate(results.state_vector()) if np.abs(amp**2) > median]
+    return [idx for idx, amp in enumerate(results.state_vector(copy=False)) if np.abs(amp**2) > median]
 
 
 @dataclass
diff --git a/cirq-core/cirq/experiments/grid_parallel_two_qubit_xeb.py b/cirq-core/cirq/experiments/grid_parallel_two_qubit_xeb.py
index c39f35bbef1..54eaee5b59b 100644
--- a/cirq-core/cirq/experiments/grid_parallel_two_qubit_xeb.py
+++ b/cirq-core/cirq/experiments/grid_parallel_two_qubit_xeb.py
@@ -496,7 +496,8 @@ def _get_xeb_result(
             while moment_index < 2 * depth:
                 step_result = next(step_results)
                 moment_index += 1
-            amplitudes = step_result.state_vector()
+            # copy=False is safe because state_vector_to_probabilities will copy anyways
+            amplitudes = step_result.state_vector(copy=False)
             probabilities = value.state_vector_to_probabilities(amplitudes)
             _, counts = np.unique(measurements, return_counts=True)
             empirical_probs = counts / len(measurements)
diff --git a/cirq-core/cirq/experiments/xeb_simulation.py b/cirq-core/cirq/experiments/xeb_simulation.py
index 72bbf252217..7fd19c530a7 100644
--- a/cirq-core/cirq/experiments/xeb_simulation.py
+++ b/cirq-core/cirq/experiments/xeb_simulation.py
@@ -65,7 +65,8 @@ def __call__(self, task: _Simulate2qXEBTask) -> List[Dict[str, Any]]:
             if cycle_depth not in cycle_depths:
                 continue
 
-            psi = step_result.state_vector()
+            # copy=False is safe because state_vector_to_probabilities will copy anyways
+            psi = step_result.state_vector(copy=False)
             pure_probs = value.state_vector_to_probabilities(psi)
 
             records += [
diff --git a/cirq-core/cirq/ops/boolean_hamiltonian_test.py b/cirq-core/cirq/ops/boolean_hamiltonian_test.py
index c4a3f781720..c48fa431f0f 100644
--- a/cirq-core/cirq/ops/boolean_hamiltonian_test.py
+++ b/cirq-core/cirq/ops/boolean_hamiltonian_test.py
@@ -80,7 +80,7 @@ def test_circuit(boolean_str):
 
     circuit.append(hamiltonian_gate.on(*qubits))
 
-    phi = cirq.Simulator().simulate(circuit, qubit_order=qubits, initial_state=0).state_vector()
+    phi = cirq.Simulator().simulate(circuit, qubit_order=qubits, initial_state=0).state_vector(copy=False)
     actual = np.arctan2(phi.real, phi.imag) - math.pi / 2.0 > 0.0
 
     # Compare the two:
diff --git a/cirq-core/cirq/ops/common_gates_test.py b/cirq-core/cirq/ops/common_gates_test.py
index bead8592348..6cdf69c53c7 100644
--- a/cirq-core/cirq/ops/common_gates_test.py
+++ b/cirq-core/cirq/ops/common_gates_test.py
@@ -1131,7 +1131,7 @@ def test_xpow_dim_3():
 
     sim = cirq.Simulator()
     circuit = cirq.Circuit([x(cirq.LineQid(0, 3)) ** 0.5] * 6)
-    svs = [step.state_vector() for step in sim.simulate_moment_steps(circuit)]
+    svs = [step.state_vector(copy=True) for step in sim.simulate_moment_steps(circuit)]
     # fmt: off
     expected = [
         [0.67, 0.67, 0.33],
@@ -1159,7 +1159,7 @@ def test_xpow_dim_4():
 
     sim = cirq.Simulator()
     circuit = cirq.Circuit([x(cirq.LineQid(0, 4)) ** 0.5] * 8)
-    svs = [step.state_vector() for step in sim.simulate_moment_steps(circuit)]
+    svs = [step.state_vector(copy=True) for step in sim.simulate_moment_steps(circuit)]
     # fmt: off
     expected = [
         [0.65, 0.65, 0.27, 0.27],
@@ -1190,11 +1190,15 @@ def test_zpow_dim_3():
 
     sim = cirq.Simulator()
     circuit = cirq.Circuit([z(cirq.LineQid(0, 3)) ** 0.5] * 6)
-    svs = [step.state_vector() for step in sim.simulate_moment_steps(circuit, initial_state=0)]
+    svs = [
+        step.state_vector(copy=True) for step in sim.simulate_moment_steps(circuit, initial_state=0)
+    ]
     expected = [[1, 0, 0]] * 6
     assert np.allclose((svs), expected)
 
-    svs = [step.state_vector() for step in sim.simulate_moment_steps(circuit, initial_state=1)]
+    svs = [
+        step.state_vector(copy=True) for step in sim.simulate_moment_steps(circuit, initial_state=1)
+    ]
     # fmt: off
     expected = [
         [0, L**0.5, 0],
@@ -1207,7 +1211,9 @@ def test_zpow_dim_3():
     # fmt: on
     assert np.allclose((svs), expected)
 
-    svs = [step.state_vector() for step in sim.simulate_moment_steps(circuit, initial_state=2)]
+    svs = [
+        step.state_vector(copy=True) for step in sim.simulate_moment_steps(circuit, initial_state=2)
+    ]
     # fmt: off
     expected = [
         [0, 0, L],
@@ -1235,11 +1241,15 @@ def test_zpow_dim_4():
 
     sim = cirq.Simulator()
     circuit = cirq.Circuit([z(cirq.LineQid(0, 4)) ** 0.5] * 8)
-    svs = [step.state_vector() for step in sim.simulate_moment_steps(circuit, initial_state=0)]
+    svs = [
+        step.state_vector(copy=True) for step in sim.simulate_moment_steps(circuit, initial_state=0)
+    ]
     expected = [[1, 0, 0, 0]] * 8
     assert np.allclose((svs), expected)
 
-    svs = [step.state_vector() for step in sim.simulate_moment_steps(circuit, initial_state=1)]
+    svs = [
+        step.state_vector(copy=True) for step in sim.simulate_moment_steps(circuit, initial_state=1)
+    ]
     # fmt: off
     expected = [
         [0, 1j**0.5, 0, 0],
@@ -1254,7 +1264,9 @@ def test_zpow_dim_4():
     # fmt: on
     assert np.allclose(svs, expected)
 
-    svs = [step.state_vector() for step in sim.simulate_moment_steps(circuit, initial_state=2)]
+    svs = [
+        step.state_vector(copy=True) for step in sim.simulate_moment_steps(circuit, initial_state=2)
+    ]
     # fmt: off
     expected = [
         [0, 0, 1j, 0],
@@ -1269,7 +1281,9 @@ def test_zpow_dim_4():
     # fmt: on
     assert np.allclose(svs, expected)
 
-    svs = [step.state_vector() for step in sim.simulate_moment_steps(circuit, initial_state=3)]
+    svs = [
+        step.state_vector(copy=True) for step in sim.simulate_moment_steps(circuit, initial_state=3)
+    ]
     # fmt: off
     expected = [
         [0, 0, 0, 1j**1.5],
diff --git a/cirq-core/cirq/sim/mux.py b/cirq-core/cirq/sim/mux.py
index 8de78f285f1..bc609b8d6b0 100644
--- a/cirq-core/cirq/sim/mux.py
+++ b/cirq-core/cirq/sim/mux.py
@@ -159,7 +159,7 @@ def final_state_vector(
         param_resolver=param_resolver,
     )
 
-    return result.state_vector()
+    return result.state_vector(copy=False)
 
 
 def sample_sweep(
diff --git a/cirq-core/cirq/sim/sparse_simulator.py b/cirq-core/cirq/sim/sparse_simulator.py
index f15a9f44703..00361bafac1 100644
--- a/cirq-core/cirq/sim/sparse_simulator.py
+++ b/cirq-core/cirq/sim/sparse_simulator.py
@@ -18,7 +18,7 @@
 
 import numpy as np
 
-from cirq import ops
+from cirq import _compat, ops
 from cirq._compat import deprecated_parameter
 from cirq.sim import simulator, state_vector, state_vector_simulator, state_vector_simulation_state
 
@@ -247,7 +247,7 @@ def __init__(
         self._dtype = dtype
         self._state_vector: Optional[np.ndarray] = None
 
-    def state_vector(self, copy: bool = True):
+    def state_vector(self, copy: Optional[bool] = None):
         """Return the state vector at this point in the computation.
 
         The state is returned in the computational basis with these basis
@@ -280,6 +280,12 @@ def state_vector(self, copy: bool = True):
                 parameters from the state vector and store then using False
                 can speed up simulation by eliminating a memory copy.
         """
+        if copy is None:
+            _compat._warn_or_error(
+                "Starting in v0.16, state_vector will not copy the state by default. "
+                "Explicitly set copy=True to copy the state."
+            )
+            copy = True
         if self._state_vector is None:
             self._state_vector = np.array([1])
             state = self._merged_sim_state
diff --git a/cirq-core/cirq/sim/sparse_simulator_test.py b/cirq-core/cirq/sim/sparse_simulator_test.py
index b63e6ce0b57..8b7c0131dc8 100644
--- a/cirq-core/cirq/sim/sparse_simulator_test.py
+++ b/cirq-core/cirq/sim/sparse_simulator_test.py
@@ -546,9 +546,9 @@ def test_simulate_moment_steps(dtype: Type[np.number], split: bool):
     simulator = cirq.Simulator(dtype=dtype, split_untangled_states=split)
     for i, step in enumerate(simulator.simulate_moment_steps(circuit)):
         if i == 0:
-            np.testing.assert_almost_equal(step.state_vector(), np.array([0.5] * 4))
+            np.testing.assert_almost_equal(step.state_vector(copy=True), np.array([0.5] * 4))
         else:
-            np.testing.assert_almost_equal(step.state_vector(), np.array([1, 0, 0, 0]))
+            np.testing.assert_almost_equal(step.state_vector(copy=True), np.array([1, 0, 0, 0]))
 
 
 @pytest.mark.parametrize('dtype', [np.complex64, np.complex128])
@@ -559,7 +559,7 @@ def test_simulate_moment_steps_empty_circuit(dtype: Type[np.number], split: bool
     step = None
     for step in simulator.simulate_moment_steps(circuit):
         pass
-    assert np.allclose(step.state_vector(), np.array([1]))
+    assert np.allclose(step.state_vector(copy=True), np.array([1]))
     assert not step.qubit_map
 
 
@@ -595,10 +595,10 @@ def test_simulate_moment_steps_intermediate_measurement(dtype: Type[np.number],
             result = int(step.measurements['0'][0])
             expected = np.zeros(2)
             expected[result] = 1
-            np.testing.assert_almost_equal(step.state_vector(), expected)
+            np.testing.assert_almost_equal(step.state_vector(copy=True), expected)
         if i == 2:
             expected = np.array([np.sqrt(0.5), np.sqrt(0.5) * (-1) ** result])
-            np.testing.assert_almost_equal(step.state_vector(), expected)
+            np.testing.assert_almost_equal(step.state_vector(copy=True), expected)
 
 
 @pytest.mark.parametrize('dtype', [np.complex64, np.complex128])
@@ -706,8 +706,8 @@ def _apply_unitary_(self, args: cirq.ApplyUnitaryArgs):
 
     initial_state = np.array([np.sqrt(0.5), np.sqrt(0.5)], dtype=np.complex64)
     result = simulator.simulate(circuit, initial_state=initial_state)
-    np.testing.assert_array_almost_equal(result.state_vector(), initial_state)
-    assert not initial_state is result.state_vector()
+    np.testing.assert_array_almost_equal(result.state_vector(copy=False), initial_state)
+    assert not initial_state is result.state_vector(copy=False)
 
 
 def test_does_not_modify_initial_state():
@@ -731,7 +731,7 @@ def _apply_unitary_(self, args: cirq.ApplyUnitaryArgs):
     result = simulator.simulate(circuit, initial_state=initial_state)
     np.testing.assert_array_almost_equal(np.array([1, 0], dtype=np.complex64), initial_state)
     np.testing.assert_array_almost_equal(
-        result.state_vector(), np.array([0, 1], dtype=np.complex64)
+        result.state_vector(copy=False), np.array([0, 1], dtype=np.complex64)
     )
 
 
@@ -762,7 +762,7 @@ def test_simulates_composite():
     c = cirq.Circuit(MultiHTestGate().on(*cirq.LineQubit.range(2)))
     expected = np.array([0.5] * 4)
     np.testing.assert_allclose(c.final_state_vector(), expected)
-    np.testing.assert_allclose(cirq.Simulator().simulate(c).state_vector(), expected)
+    np.testing.assert_allclose(cirq.Simulator().simulate(c).state_vector(copy=False), expected)
 
 
 def test_simulate_measurement_inversions():
@@ -779,7 +779,7 @@ def test_works_on_pauli_string_phasor():
     a, b = cirq.LineQubit.range(2)
     c = cirq.Circuit(np.exp(0.5j * np.pi * cirq.X(a) * cirq.X(b)))
     sim = cirq.Simulator()
-    result = sim.simulate(c).state_vector()
+    result = sim.simulate(c).state_vector(copy=False)
     np.testing.assert_allclose(result.reshape(4), np.array([0, 0, 0, 1j]), atol=1e-8)
 
 
@@ -787,7 +787,7 @@ def test_works_on_pauli_string():
     a, b = cirq.LineQubit.range(2)
     c = cirq.Circuit(cirq.X(a) * cirq.X(b))
     sim = cirq.Simulator()
-    result = sim.simulate(c).state_vector()
+    result = sim.simulate(c).state_vector(copy=False)
     np.testing.assert_allclose(result.reshape(4), np.array([0, 0, 0, 1]), atol=1e-8)
 
 
@@ -1297,9 +1297,9 @@ def test_final_state_vector_is_not_last_object():
     initial_state = np.array([1, 0], dtype=np.complex64)
     circuit = cirq.Circuit(cirq.wait(q))
     result = sim.simulate(circuit, initial_state=initial_state)
-    assert result.state_vector() is not initial_state
-    assert not np.shares_memory(result.state_vector(), initial_state)
-    np.testing.assert_equal(result.state_vector(), initial_state)
+    assert result.state_vector(copy=False) is not initial_state
+    assert not np.shares_memory(result.state_vector(copy=False), initial_state)
+    np.testing.assert_equal(result.state_vector(copy=False), initial_state)
 
 
 def test_deterministic_gate_noise():
diff --git a/cirq-core/cirq/sim/state_vector.py b/cirq-core/cirq/sim/state_vector.py
index ebccd000ca0..dc1dec16a38 100644
--- a/cirq-core/cirq/sim/state_vector.py
+++ b/cirq-core/cirq/sim/state_vector.py
@@ -57,7 +57,7 @@ def _qid_shape_(self) -> Tuple[int, ...]:
         return self._qid_shape
 
     @abc.abstractmethod
-    def state_vector(self) -> np.ndarray:
+    def state_vector(self, copy: Optional[bool] = False) -> np.ndarray:
         """Return the state vector (wave function).
 
         The vector is returned in the computational basis with these basis
@@ -83,6 +83,12 @@ def state_vector(self) -> np.ndarray:
                 |  6  |   1    |   1    |   0    |
                 |  7  |   1    |   1    |   1    |
 
+        Args:
+            copy: If True, the returned state vector will be a copy of that
+            stored by the object. This is potentially expensive for large
+            state vectors, but prevents mutation of the object state, e.g. for
+            operating on intermediate states of a circuit.
+            Defaults to False.
         """
         raise NotImplementedError()
 
@@ -95,19 +101,21 @@ def dirac_notation(self, decimals: int = 2) -> str:
         Returns:
             A pretty string consisting of a sum of computational basis kets
             and non-zero floats of the specified accuracy."""
-        return qis.dirac_notation(self.state_vector(), decimals, qid_shape=self._qid_shape)
+        return qis.dirac_notation(
+            self.state_vector(copy=False), decimals, qid_shape=self._qid_shape
+        )
 
     def density_matrix_of(self, qubits: List['cirq.Qid'] = None) -> np.ndarray:
         r"""Returns the density matrix of the state.
 
         Calculate the density matrix for the system on the list, qubits.
-        Any qubits not in the list that are present in self.state_vector() will
+        Any qubits not in the list that are present in self.state_vector(copy=False) will
         be traced out. If qubits is None the full density matrix for
-        self.state_vector() is returned, given self.state_vector() follows
+        self.state_vector(copy=False) is returned, given self.state_vector(copy=False) follows
         standard Kronecker convention of numpy.kron.
 
         For example:
-        self.state_vector() = np.array([1/np.sqrt(2), 1/np.sqrt(2)],
+        self.state_vector(copy=False) = np.array([1/np.sqrt(2), 1/np.sqrt(2)],
             dtype=np.complex64)
         qubits = None
         gives us
@@ -132,7 +140,7 @@ def density_matrix_of(self, qubits: List['cirq.Qid'] = None) -> np.ndarray:
                 corresponding to the state.
         """
         return qis.density_matrix_from_state_vector(
-            self.state_vector(),
+            self.state_vector(copy=False),
             [self.qubit_map[q] for q in qubits] if qubits is not None else None,
             qid_shape=self._qid_shape,
         )
@@ -141,8 +149,8 @@ def bloch_vector_of(self, qubit: 'cirq.Qid') -> np.ndarray:
         """Returns the bloch vector of a qubit in the state.
 
         Calculates the bloch vector of the given qubit
-        in the state given by self.state_vector(), given that
-        self.state_vector() follows the standard Kronecker convention of
+        in the state given by self.state_vector(copy=False), given that
+        self.state_vector(copy=False) follows the standard Kronecker convention of
         numpy.kron.
 
         Args:
@@ -157,7 +165,7 @@ def bloch_vector_of(self, qubit: 'cirq.Qid') -> np.ndarray:
                 corresponding to the state.
         """
         return qis.bloch_vector_from_state_vector(
-            self.state_vector(), self.qubit_map[qubit], qid_shape=self._qid_shape
+            self.state_vector(copy=False), self.qubit_map[qubit], qid_shape=self._qid_shape
         )
 
 
diff --git a/cirq-core/cirq/sim/state_vector_simulator.py b/cirq-core/cirq/sim/state_vector_simulator.py
index cbf965dc99d..9fd1379a776 100644
--- a/cirq-core/cirq/sim/state_vector_simulator.py
+++ b/cirq-core/cirq/sim/state_vector_simulator.py
@@ -29,7 +29,7 @@
 
 import numpy as np
 
-from cirq import ops, value, qis
+from cirq import _compat, ops, value, qis
 from cirq._compat import deprecated_class, proper_repr
 from cirq.sim import simulator, state_vector, simulator_base
 
@@ -163,7 +163,7 @@ def final_state_vector(self) -> np.ndarray:
             self._final_state_vector = tensor
         return self._final_state_vector
 
-    def state_vector(self) -> np.ndarray:
+    def state_vector(self, copy: bool = None) -> np.ndarray:
         """Return the state vector at the end of the computation.
 
         The state is returned in the computational basis with these basis
@@ -189,7 +189,13 @@ def state_vector(self) -> np.ndarray:
                 |  6  |   1    |   1    |   0    |
                 |  7  |   1    |   1    |   1    |
         """
-        return self.final_state_vector.copy()
+        if copy is None:
+            _compat._warn_or_error(
+                "Starting in v0.16, state_vector will not copy the state by default. "
+                "Explicitly set copy=True to copy the state."
+            )
+            copy = True
+        return self.final_state_vector.copy() if copy else self.final_state_vector
 
     def _value_equality_values_(self):
         measurements = {k: v.tolist() for k, v in sorted(self.measurements.items())}
diff --git a/cirq-core/cirq/sim/state_vector_simulator_test.py b/cirq-core/cirq/sim/state_vector_simulator_test.py
index 6bff10df026..774817f38bc 100644
--- a/cirq-core/cirq/sim/state_vector_simulator_test.py
+++ b/cirq-core/cirq/sim/state_vector_simulator_test.py
@@ -140,7 +140,7 @@ def test_state_vector_trial_state_vector_is_copy():
     trial_result = cirq.StateVectorTrialResult(
         params=cirq.ParamResolver({}), measurements={}, final_simulator_state=final_simulator_state
     )
-    assert trial_result.state_vector() is not final_simulator_state.target_tensor
+    assert trial_result.state_vector(copy=True) is not final_simulator_state.target_tensor
 
 
 def test_str_big():
diff --git a/cirq-core/cirq/sim/state_vector_test.py b/cirq-core/cirq/sim/state_vector_test.py
index d8b4df663de..4f5f2ea342c 100644
--- a/cirq-core/cirq/sim/state_vector_test.py
+++ b/cirq-core/cirq/sim/state_vector_test.py
@@ -14,6 +14,7 @@
 """Tests for state_vector.py"""
 
 import itertools
+from typing import Optional
 import pytest
 
 import numpy as np
@@ -24,7 +25,7 @@
 
 def test_state_mixin():
     class TestClass(cirq.StateVectorMixin):
-        def state_vector(self) -> np.ndarray:
+        def state_vector(self, copy: Optional[bool] = None) -> np.ndarray:
             return np.array([0, 0, 1, 0])
 
     qubits = cirq.LineQubit.range(2)
@@ -330,7 +331,7 @@ def test_measure_state_empty_state():
 
 
 class BasicStateVector(cirq.StateVectorMixin):
-    def state_vector(self) -> np.ndarray:
+    def state_vector(self, copy: Optional[bool] = None) -> np.ndarray:
         return np.array([0, 1, 0, 0])
 
 
diff --git a/cirq-google/cirq_google/calibration/engine_simulator.py b/cirq-google/cirq_google/calibration/engine_simulator.py
index 58986a6eb24..00ccb6e36fd 100644
--- a/cirq-google/cirq_google/calibration/engine_simulator.py
+++ b/cirq-google/cirq_google/calibration/engine_simulator.py
@@ -386,7 +386,7 @@ def create_from_characterizations_sqrt_iswap(
 
     def final_state_vector(self, program: cirq.Circuit) -> np.ndarray:
         result = self.simulate(program)
-        return result.state_vector()
+        return result.state_vector(copy=False)
 
     def get_calibrations(
         self, requests: Sequence[PhasedFSimCalibrationRequest]
diff --git a/docs/simulation.ipynb b/docs/simulation.ipynb
index b1349748e82..6c62a8ab261 100644
--- a/docs/simulation.ipynb
+++ b/docs/simulation.ipynb
@@ -518,7 +518,7 @@
     "circuit = cirq.Circuit()\n",
     "circuit.append(basic_circuit())    \n",
     "for i, step in enumerate(simulator.simulate_moment_steps(circuit)):\n",
-    "    print('state at step %d: %s' % (i, np.around(step.state_vector(), 3)))"
+    "    print('state at step %d: %s' % (i, np.around(step.state_vector(copy=True), 3)))"
    ]
   },
   {
diff --git a/docs/tutorials/educators/intro.ipynb b/docs/tutorials/educators/intro.ipynb
index ca371aa4f7e..12634b3a2d0 100644
--- a/docs/tutorials/educators/intro.ipynb
+++ b/docs/tutorials/educators/intro.ipynb
@@ -1811,7 +1811,7 @@
     "\n",
     "# Step through the simulation results.\n",
     "for step in simulator.simulate_moment_steps(circuit):\n",
-    "    prob = np.abs(step.state_vector()) ** 2\n",
+    "    prob = np.abs(step.state_vector(copy=True)) ** 2\n",
     "    probs.append(prob[0])\n",
     "\n",
     "# Plot the probability of the ground state at each simulation step.\n",
diff --git a/docs/tutorials/fourier_checking.ipynb b/docs/tutorials/fourier_checking.ipynb
index 35e30098377..3e3d13938b7 100644
--- a/docs/tutorials/fourier_checking.ipynb
+++ b/docs/tutorials/fourier_checking.ipynb
@@ -743,7 +743,7 @@
     "for step in s.simulate_moment_steps(circuit):\n",
     "    print(step.dirac_notation())\n",
     "    print(\"|0> state probability to observe: \",\n",
-    "          np.abs(step.state_vector()[0])**2)"
+    "          np.abs(step.state_vector(copy=True)[0])**2)"
    ]
   },
   {
diff --git a/examples/stabilizer_code.ipynb b/examples/stabilizer_code.ipynb
index 63edf1492f6..d2e99448681 100644
--- a/examples/stabilizer_code.ipynb
+++ b/examples/stabilizer_code.ipynb
@@ -101,7 +101,7 @@
    "source": [
     "for initial_state in [0, 1]:\n",
     "    results = cirq.Simulator().simulate(circuit, qubit_order=qubits, initial_state=initial_state)\n",
-    "    print('%d becomes %s' % (initial_state, cirq.dirac_notation(results.state_vector())))\n",
+    "    print('%d becomes %s' % (initial_state, cirq.dirac_notation(results.state_vector(copy=False))))\n",
     "    \n",
     "    "
    ]
@@ -160,7 +160,7 @@
    ],
    "source": [
     "results = cirq.Simulator().simulate(circuit, qubit_order=(qubits+ancillas), initial_state=0)\n",
-    "print('%s' % (cirq.dirac_notation(results.state_vector())))"
+    "print('%s' % (cirq.dirac_notation(results.state_vector(copy=False))))"
    ]
   },
   {
@@ -216,7 +216,7 @@
     }
    ],
    "source": [
-    "decoded = code.decode(qubits, ancillas, results.state_vector())\n",
+    "decoded = code.decode(qubits, ancillas, results.state_vector(copy=False))\n",
     "print('decoded=%s' % (decoded))"
    ]
   },
@@ -256,7 +256,7 @@
     "\n",
     "for initial_state in [0, 1]:\n",
     "    results = cirq.Simulator().simulate(circuit, qubit_order=qubits, initial_state=initial_state)\n",
-    "    print('%d becomes:\\n %s\\n' % (initial_state, cirq.dirac_notation(4.0 * results.state_vector())))"
+    "    print('%d becomes:\\n %s\\n' % (initial_state, cirq.dirac_notation(4.0 * results.state_vector(copy=False))))"
    ]
   }
  ],
diff --git a/examples/stabilizer_code_test.py b/examples/stabilizer_code_test.py
index ddd28e2b71a..8a8615c73a8 100644
--- a/examples/stabilizer_code_test.py
+++ b/examples/stabilizer_code_test.py
@@ -37,7 +37,7 @@ def encode_corrupt_correct(
         .simulate(
             circuit, qubit_order=(qubits + ancillas), initial_state=(input_val * 2 ** len(ancillas))
         )
-        .state_vector()
+        .state_vector(copy=False)
     )
 
     decoded = code.decode(qubits, ancillas, state_vector)

From 6efcdb99d251ceae79721aa698cd94d81d158916 Mon Sep 17 00:00:00 2001
From: Orion Martin <40585662+95-martin-orion@users.noreply.github.com>
Date: Tue, 3 May 2022 11:36:07 -0700
Subject: [PATCH 2/3] Format and docstring fixes.

---
 .../cirq/contrib/quantum_volume/quantum_volume.py      |  4 +++-
 cirq-core/cirq/ops/boolean_hamiltonian_test.py         |  6 +++++-
 cirq-core/cirq/sim/state_vector.py                     | 10 +++++-----
 3 files changed, 13 insertions(+), 7 deletions(-)

diff --git a/cirq-core/cirq/contrib/quantum_volume/quantum_volume.py b/cirq-core/cirq/contrib/quantum_volume/quantum_volume.py
index 0c5ca5a1fa8..9dd947727e5 100644
--- a/cirq-core/cirq/contrib/quantum_volume/quantum_volume.py
+++ b/cirq-core/cirq/contrib/quantum_volume/quantum_volume.py
@@ -86,7 +86,9 @@ def compute_heavy_set(circuit: cirq.Circuit) -> List[int]:
     # The output wave function is a vector from the result value (big-endian) to
     # the probability of that bit-string. Return all of the bit-string
     # values that have a probability greater than the median.
-    return [idx for idx, amp in enumerate(results.state_vector(copy=False)) if np.abs(amp**2) > median]
+    return [
+        idx for idx, amp in enumerate(results.state_vector(copy=False)) if np.abs(amp**2) > median
+    ]
 
 
 @dataclass
diff --git a/cirq-core/cirq/ops/boolean_hamiltonian_test.py b/cirq-core/cirq/ops/boolean_hamiltonian_test.py
index c48fa431f0f..426f7a6b9eb 100644
--- a/cirq-core/cirq/ops/boolean_hamiltonian_test.py
+++ b/cirq-core/cirq/ops/boolean_hamiltonian_test.py
@@ -80,7 +80,11 @@ def test_circuit(boolean_str):
 
     circuit.append(hamiltonian_gate.on(*qubits))
 
-    phi = cirq.Simulator().simulate(circuit, qubit_order=qubits, initial_state=0).state_vector(copy=False)
+    phi = (
+        cirq.Simulator()
+        .simulate(circuit, qubit_order=qubits, initial_state=0)
+        .state_vector(copy=False)
+    )
     actual = np.arctan2(phi.real, phi.imag) - math.pi / 2.0 > 0.0
 
     # Compare the two:
diff --git a/cirq-core/cirq/sim/state_vector.py b/cirq-core/cirq/sim/state_vector.py
index dc1dec16a38..3b9ab25effe 100644
--- a/cirq-core/cirq/sim/state_vector.py
+++ b/cirq-core/cirq/sim/state_vector.py
@@ -109,13 +109,13 @@ def density_matrix_of(self, qubits: List['cirq.Qid'] = None) -> np.ndarray:
         r"""Returns the density matrix of the state.
 
         Calculate the density matrix for the system on the list, qubits.
-        Any qubits not in the list that are present in self.state_vector(copy=False) will
+        Any qubits not in the list that are present in self.state_vector() will
         be traced out. If qubits is None the full density matrix for
-        self.state_vector(copy=False) is returned, given self.state_vector(copy=False) follows
+        self.state_vector() is returned, given self.state_vector() follows
         standard Kronecker convention of numpy.kron.
 
         For example:
-        self.state_vector(copy=False) = np.array([1/np.sqrt(2), 1/np.sqrt(2)],
+        self.state_vector() = np.array([1/np.sqrt(2), 1/np.sqrt(2)],
             dtype=np.complex64)
         qubits = None
         gives us
@@ -149,8 +149,8 @@ def bloch_vector_of(self, qubit: 'cirq.Qid') -> np.ndarray:
         """Returns the bloch vector of a qubit in the state.
 
         Calculates the bloch vector of the given qubit
-        in the state given by self.state_vector(copy=False), given that
-        self.state_vector(copy=False) follows the standard Kronecker convention of
+        in the state given by self.state_vector(), given that
+        self.state_vector() follows the standard Kronecker convention of
         numpy.kron.
 
         Args:

From 5219eb92bed384bfbff093adfef322d7715b586b Mon Sep 17 00:00:00 2001
From: Orion Martin <40585662+95-martin-orion@users.noreply.github.com>
Date: Tue, 3 May 2022 12:58:48 -0700
Subject: [PATCH 3/3] Catch missed tests.

---
 cirq-core/cirq/sim/simulator_test.py             |  4 +++-
 cirq-core/cirq/sim/sparse_simulator_test.py      | 11 +++++++++++
 .../cirq/sim/state_vector_simulator_test.py      | 16 ++++++++++++++++
 3 files changed, 30 insertions(+), 1 deletion(-)

diff --git a/cirq-core/cirq/sim/simulator_test.py b/cirq-core/cirq/sim/simulator_test.py
index 59b4d99e5d9..c71e35138fb 100644
--- a/cirq-core/cirq/sim/simulator_test.py
+++ b/cirq-core/cirq/sim/simulator_test.py
@@ -416,7 +416,9 @@ def _kraus_(self):
             cirq.Circuit(Reset11To00().on(*cirq.LineQubit.range(2))), initial_state=k
         )
         np.testing.assert_allclose(
-            out.state_vector(), cirq.one_hot(index=k % 3, shape=4, dtype=np.complex64), atol=1e-8
+            out.state_vector(copy=False),
+            cirq.one_hot(index=k % 3, shape=4, dtype=np.complex64),
+            atol=1e-8,
         )
 
 
diff --git a/cirq-core/cirq/sim/sparse_simulator_test.py b/cirq-core/cirq/sim/sparse_simulator_test.py
index 8b7c0131dc8..905e3238f2c 100644
--- a/cirq-core/cirq/sim/sparse_simulator_test.py
+++ b/cirq-core/cirq/sim/sparse_simulator_test.py
@@ -551,6 +551,17 @@ def test_simulate_moment_steps(dtype: Type[np.number], split: bool):
             np.testing.assert_almost_equal(step.state_vector(copy=True), np.array([1, 0, 0, 0]))
 
 
+def test_simulate_moment_steps_implicit_copy_deprecated():
+    q0 = cirq.LineQubit(0)
+    simulator = cirq.Simulator()
+    steps = list(simulator.simulate_moment_steps(cirq.Circuit(cirq.X(q0))))
+
+    with cirq.testing.assert_deprecated(
+        "state_vector will not copy the state by default", deadline="v0.16"
+    ):
+        _ = steps[0].state_vector()
+
+
 @pytest.mark.parametrize('dtype', [np.complex64, np.complex128])
 @pytest.mark.parametrize('split', [True, False])
 def test_simulate_moment_steps_empty_circuit(dtype: Type[np.number], split: bool):
diff --git a/cirq-core/cirq/sim/state_vector_simulator_test.py b/cirq-core/cirq/sim/state_vector_simulator_test.py
index 774817f38bc..b98353a544a 100644
--- a/cirq-core/cirq/sim/state_vector_simulator_test.py
+++ b/cirq-core/cirq/sim/state_vector_simulator_test.py
@@ -143,6 +143,22 @@ def test_state_vector_trial_state_vector_is_copy():
     assert trial_result.state_vector(copy=True) is not final_simulator_state.target_tensor
 
 
+def test_implicit_copy_deprecated():
+    final_state_vector = np.array([0, 1], dtype=np.complex64)
+    qubit_map = {cirq.NamedQubit('a'): 0}
+    final_simulator_state = cirq.StateVectorSimulationState(
+        qubits=list(qubit_map), initial_state=final_state_vector
+    )
+    trial_result = cirq.StateVectorTrialResult(
+        params=cirq.ParamResolver({}), measurements={}, final_simulator_state=final_simulator_state
+    )
+
+    with cirq.testing.assert_deprecated(
+        "state_vector will not copy the state by default", deadline="v0.16"
+    ):
+        _ = trial_result.state_vector()
+
+
 def test_str_big():
     qs = cirq.LineQubit.range(10)
     final_simulator_state = cirq.StateVectorSimulationState(