From 046b8840058978ee4c061145a279f2873e687175 Mon Sep 17 00:00:00 2001
From: Orion Martin <40585662+95-martin-orion@users.noreply.github.com>
Date: Tue, 7 Dec 2021 11:34:20 -0800
Subject: [PATCH 01/17] Remove files for subsequent PRs

---
 cirq-core/cirq/__init__.py                    |   1 +
 cirq-core/cirq/devices/__init__.py            |  12 +
 cirq-core/cirq/devices/noise_utils.py         | 232 ++++++++++++++++++
 cirq-core/cirq/devices/noise_utils_test.py    | 140 +++++++++++
 cirq-core/cirq/json_resolver_cache.py         |   1 +
 .../json_test_data/OpIdentifier.json          |  10 +
 .../json_test_data/OpIdentifier.repr          |   4 +
 7 files changed, 400 insertions(+)
 create mode 100644 cirq-core/cirq/devices/noise_utils.py
 create mode 100644 cirq-core/cirq/devices/noise_utils_test.py
 create mode 100644 cirq-core/cirq/protocols/json_test_data/OpIdentifier.json
 create mode 100644 cirq-core/cirq/protocols/json_test_data/OpIdentifier.repr

diff --git a/cirq-core/cirq/__init__.py b/cirq-core/cirq/__init__.py
index 35b84ca21da..b9e3b5121b1 100644
--- a/cirq-core/cirq/__init__.py
+++ b/cirq-core/cirq/__init__.py
@@ -88,6 +88,7 @@
     NO_NOISE,
     NOISE_MODEL_LIKE,
     NoiseModel,
+    OpIdentifier,
     SymmetricalQidPair,
     UNCONSTRAINED_DEVICE,
     NamedTopology,
diff --git a/cirq-core/cirq/devices/__init__.py b/cirq-core/cirq/devices/__init__.py
index 1060f7aea0e..72f793a7e46 100644
--- a/cirq-core/cirq/devices/__init__.py
+++ b/cirq-core/cirq/devices/__init__.py
@@ -47,3 +47,15 @@
     get_placements,
     draw_placements,
 )
+
+from cirq.devices.noise_utils import (
+    OpIdentifier,
+    decay_constant_to_xeb_fidelity,
+    decay_constant_to_pauli_error,
+    pauli_error_to_decay_constant,
+    xeb_fidelity_to_decay_constant,
+    pauli_error_from_t1,
+    pauli_error_from_depolarization,
+    average_error,
+    decoherence_pauli_error,
+)
diff --git a/cirq-core/cirq/devices/noise_utils.py b/cirq-core/cirq/devices/noise_utils.py
new file mode 100644
index 00000000000..e02ed8baf45
--- /dev/null
+++ b/cirq-core/cirq/devices/noise_utils.py
@@ -0,0 +1,232 @@
+# Copyright 2021 The Cirq Developers
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     https://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from typing import TYPE_CHECKING, Any, Dict, Tuple, Type, Union
+import warnings
+import numpy as np
+
+from cirq import ops, protocols, value
+
+if TYPE_CHECKING:
+    import cirq
+
+
+# Tag for gates to which noise must be applied.
+PHYSICAL_GATE_TAG = 'physical_gate'
+
+
+@value.value_equality(distinct_child_types=True)
+class OpIdentifier:
+    """Identifies an operation by gate and (optionally) target qubits."""
+
+    def __init__(self, gate_type: Type['cirq.Gate'], *qubits: 'cirq.Qid'):
+        self._gate_type = gate_type
+        self._gate_family = ops.GateFamily(gate_type)
+        self._qubits: Tuple['cirq.Qid', ...] = tuple(qubits)
+
+    @property
+    def gate_type(self) -> Type['cirq.Gate']:
+        # set to a type during initialization, never modified
+        return self._gate_type
+
+    @property
+    def qubits(self) -> Tuple['cirq.Qid', ...]:
+        return self._qubits
+
+    def _predicate(self, *args, **kwargs):
+        return self._gate_family._predicate(*args, **kwargs)
+
+    def swapped(self):
+        return OpIdentifier(self.gate_type, *self.qubits[::-1])
+
+    def is_proper_subtype_of(self, op_id: 'OpIdentifier'):
+        """Returns true if this is contained within op_id, but not equal to it.
+
+        If this returns true, (x in self) implies (x in op_id), but the reverse
+        implication does not hold. op_id must be more general than self (either
+        by accepting any qubits or having a more general gate type) for this
+        to return true.
+        """
+        more_specific_qubits = self.qubits and not op_id.qubits
+        more_specific_gate = self.gate_type != op_id.gate_type and issubclass(
+            self.gate_type, op_id.gate_type
+        )
+        return (
+            (more_specific_qubits or more_specific_gate)
+            and (more_specific_qubits or self.qubits == op_id.qubits)
+            and (more_specific_gate or self.gate_type == op_id.gate_type)
+        )
+
+    def __contains__(self, item: Union[ops.Gate, ops.Operation]) -> bool:
+        if isinstance(item, ops.Gate):
+            return (not self._qubits) and self._predicate(item)
+        return (
+            (not self.qubits or (item.qubits == self._qubits))
+            and item.gate is not None
+            and self._predicate(item.gate)
+        )
+
+    def __str__(self):
+        return f'{self.gate_type}{self.qubits}'
+
+    def __repr__(self) -> str:
+        fullname = f'{self.gate_type.__module__}.{self.gate_type.__qualname__}'
+        qubits = ', '.join(map(repr, self.qubits))
+        return f'cirq.devices.noise_utils.OpIdentifier({fullname}, {qubits})'
+
+    def _value_equality_values_(self) -> Any:
+        return (self.gate_type, self.qubits)
+
+    def _json_dict_(self) -> Dict[str, Any]:
+        gate_json = protocols.json_cirq_type(self._gate_type)
+        return {
+            'gate_type': gate_json,
+            'qubits': self._qubits,
+        }
+
+    @classmethod
+    def _from_json_dict_(cls, gate_type, qubits, **kwargs) -> 'OpIdentifier':
+        gate_type = protocols.cirq_type_from_json(gate_type)
+        return cls(gate_type, *qubits)
+
+
+# TODO: expose all from top-level cirq?
+def decay_constant_to_xeb_fidelity(decay_constant: float, num_qubits: int = 2) -> float:
+    """Calculates the XEB fidelity from the depolarization decay constant.
+
+    Args:
+        decay_constant: Depolarization decay constant.
+        num_qubits: Number of qubits.
+
+    Returns:
+        Calculated XEB fidelity.
+    """
+    N = 2 ** num_qubits
+    return 1 - ((1 - decay_constant) * (1 - 1 / N))
+
+
+def decay_constant_to_pauli_error(decay_constant: float, num_qubits: int = 1) -> float:
+    """Calculates pauli error from the depolarization decay constant.
+
+    Args:
+        decay_constant: Depolarization decay constant.
+        num_qubits: Number of qubits.
+
+    Returns:
+        Calculated Pauli error.
+    """
+    N = 2 ** num_qubits
+    return (1 - decay_constant) * (1 - 1 / N / N)
+
+
+def pauli_error_to_decay_constant(pauli_error: float, num_qubits: int = 1) -> float:
+    """Calculates depolarization decay constant from pauli error.
+
+    Args:
+        pauli_error: The pauli error.
+        num_qubits: Number of qubits.
+
+    Returns:
+        Calculated depolarization decay constant.
+    """
+    N = 2 ** num_qubits
+    return 1 - (pauli_error / (1 - 1 / N / N))
+
+
+def xeb_fidelity_to_decay_constant(xeb_fidelity: float, num_qubits: int = 2) -> float:
+    """Calculates the depolarization decay constant from XEB fidelity.
+
+    Args:
+        xeb_fidelity: The XEB fidelity.
+        num_qubits: Number of qubits.
+
+    Returns:
+        Calculated depolarization decay constant.
+    """
+    N = 2 ** num_qubits
+    return 1 - (1 - xeb_fidelity) / (1 - 1 / N)
+
+
+def pauli_error_from_t1(t_ns: float, t1_ns: float) -> float:
+    """Calculates the pauli error from T1 decay constant.
+
+    This computes error for a specific duration, `t`.
+
+    Args:
+        t_ns: The duration of the gate in ns.
+        t1_ns: The T1 decay constant in ns.
+
+    Returns:
+        Calculated Pauli error resulting from T1 decay.
+    """
+    t2 = 2 * t1_ns
+    return (1 - np.exp(-t_ns / t2)) / 2 + (1 - np.exp(-t_ns / t1_ns)) / 4
+
+
+def pauli_error_from_depolarization(t_ns: float, t1_ns: float, pauli_error: float = 0) -> float:
+    """Calculates the amount of pauli error from depolarization.
+
+    This computes non-T1 error for a specific duration, `t`. If pauli error
+    from T1 decay is more than total pauli error, this returns zero; otherwise,
+    it returns the portion of pauli error not attributable to T1 error.
+
+    Args:
+        t_ns: The duration of the gate in ns.
+        t1_ns: The T1 decay constant in ns.
+        pauli_error: The total pauli error.
+
+    Returns:
+        Calculated Pauli error resulting from depolarization.
+    """
+    t1_pauli_error = pauli_error_from_t1(t_ns, t1_ns)
+    if pauli_error >= t1_pauli_error:
+        return pauli_error - t1_pauli_error
+
+    warnings.warn("Pauli error from T1 decay is greater than total Pauli error", RuntimeWarning)
+    return 0
+
+
+def average_error(decay_constant: float, num_qubits: int = 1) -> float:
+    """Calculates the average error from the depolarization decay constant.
+
+    Args:
+        decay_constant: Depolarization decay constant.
+        num_qubits: Number of qubits.
+
+    Returns:
+        Calculated average error.
+    """
+    N = 2 ** num_qubits
+    return (1 - decay_constant) * (1 - 1 / N)
+
+
+def decoherence_pauli_error(t1_ns: float, tphi_ns: float, gate_time_ns: float) -> float:
+    """The component of Pauli error caused by decoherence.
+
+    Args:
+        t1_ns: T1 time in nanoseconds.
+        tphi_ns: Tphi time in nanoseconds.
+        gate_time_ns: Duration in nanoseconds of the gate affected by this error.
+
+    Returns:
+        Calculated Pauli error resulting from decoherence.
+    """
+    Gamma2 = (1 / (2 * t1_ns)) + 1 / tphi_ns
+
+    exp1 = np.exp(-gate_time_ns / t1_ns)
+    exp2 = np.exp(-gate_time_ns * Gamma2)
+    px = 0.25 * (1 - exp1)
+    py = px
+    pz = 0.5 * (1 - exp2) - px
+    return px + py + pz
diff --git a/cirq-core/cirq/devices/noise_utils_test.py b/cirq-core/cirq/devices/noise_utils_test.py
new file mode 100644
index 00000000000..61b5852ced4
--- /dev/null
+++ b/cirq-core/cirq/devices/noise_utils_test.py
@@ -0,0 +1,140 @@
+# Copyright 2021 The Cirq Developers
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     https://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import numpy as np
+import pytest
+
+import cirq
+from cirq.devices.noise_utils import (
+    OpIdentifier,
+    decay_constant_to_xeb_fidelity,
+    decay_constant_to_pauli_error,
+    pauli_error_to_decay_constant,
+    xeb_fidelity_to_decay_constant,
+    pauli_error_from_t1,
+    pauli_error_from_depolarization,
+    average_error,
+    decoherence_pauli_error,
+)
+
+
+def test_op_id():
+    op_id = OpIdentifier(cirq.XPowGate)
+    assert cirq.X(cirq.LineQubit(1)) in op_id
+    assert cirq.Rx(rads=1) in op_id
+
+
+@pytest.mark.parametrize(
+    'decay_constant,num_qubits,expected_output',
+    [
+        (0.01, 1, 1 - (0.99 * 1 / 2)),
+        (0.05, 2, 1 - (0.95 * 3 / 4)),
+    ],
+)
+def test_decay_constant_to_xeb_fidelity(decay_constant, num_qubits, expected_output):
+    val = decay_constant_to_xeb_fidelity(decay_constant, num_qubits)
+    assert val == expected_output
+
+
+@pytest.mark.parametrize(
+    'decay_constant,num_qubits,expected_output',
+    [
+        (0.01, 1, 0.99 * 3 / 4),
+        (0.05, 2, 0.95 * 15 / 16),
+    ],
+)
+def test_decay_constant_to_pauli_error(decay_constant, num_qubits, expected_output):
+    val = decay_constant_to_pauli_error(decay_constant, num_qubits)
+    assert val == expected_output
+
+
+@pytest.mark.parametrize(
+    'pauli_error,num_qubits,expected_output',
+    [
+        (0.01, 1, 1 - (0.01 / (3 / 4))),
+        (0.05, 2, 1 - (0.05 / (15 / 16))),
+    ],
+)
+def test_pauli_error_to_decay_constant(pauli_error, num_qubits, expected_output):
+    val = pauli_error_to_decay_constant(pauli_error, num_qubits)
+    assert val == expected_output
+
+
+@pytest.mark.parametrize(
+    'xeb_fidelity,num_qubits,expected_output',
+    [
+        (0.01, 1, 1 - 0.99 / (1 / 2)),
+        (0.05, 2, 1 - 0.95 / (3 / 4)),
+    ],
+)
+def test_xeb_fidelity_to_decay_constant(xeb_fidelity, num_qubits, expected_output):
+    val = xeb_fidelity_to_decay_constant(xeb_fidelity, num_qubits)
+    assert val == expected_output
+
+
+@pytest.mark.parametrize(
+    't,t1_ns,expected_output',
+    [
+        (20, 1e5, (1 - np.exp(-20 / 2e5)) / 2 + (1 - np.exp(-20 / 1e5)) / 4),
+        (4000, 1e4, (1 - np.exp(-4000 / 2e4)) / 2 + (1 - np.exp(-4000 / 1e4)) / 4),
+    ],
+)
+def test_pauli_error_from_t1(t, t1_ns, expected_output):
+    val = pauli_error_from_t1(t, t1_ns)
+    assert val == expected_output
+
+
+@pytest.mark.parametrize(
+    't,t1_ns,pauli_error,expected_output',
+    [
+        (20, 1e5, 0.01, 0.01 - ((1 - np.exp(-20 / 2e5)) / 2 + (1 - np.exp(-20 / 1e5)) / 4)),
+        # In this case, the formula produces a negative result.
+        (4000, 1e4, 0.01, 0),
+    ],
+)
+def test_pauli_error_from_depolarization(t, t1_ns, pauli_error, expected_output):
+    val = pauli_error_from_depolarization(t, t1_ns, pauli_error)
+    assert val == expected_output
+
+
+@pytest.mark.parametrize(
+    'decay_constant,num_qubits,expected_output',
+    [
+        (0.01, 1, 0.99 * 1 / 2),
+        (0.05, 2, 0.95 * 3 / 4),
+    ],
+)
+def test_average_error(decay_constant, num_qubits, expected_output):
+    val = average_error(decay_constant, num_qubits)
+    assert val == expected_output
+
+
+@pytest.mark.parametrize(
+    'T1_ns,Tphi_ns,gate_time_ns',
+    [
+        (1e4, 2e4, 25),
+        (1e5, 2e3, 25),
+        (1e4, 2e4, 4000),
+    ],
+)
+def test_decoherence_pauli_error(T1_ns, Tphi_ns, gate_time_ns):
+    val = decoherence_pauli_error(T1_ns, Tphi_ns, gate_time_ns)
+    # Expected value is of the form:
+    #
+    #   (1/4) * [1 - e^(-t/T1)] + (1/2) * [1 - e^(-t/(2*T1) - t/Tphi]
+    #
+    expected_output = 0.25 * (1 - np.exp(-gate_time_ns / T1_ns)) + 0.5 * (
+        1 - np.exp(-gate_time_ns * ((1 / (2 * T1_ns)) + 1 / Tphi_ns))
+    )
+    assert val == expected_output
diff --git a/cirq-core/cirq/json_resolver_cache.py b/cirq-core/cirq/json_resolver_cache.py
index f18fbc90c91..3692c3b191d 100644
--- a/cirq-core/cirq/json_resolver_cache.py
+++ b/cirq-core/cirq/json_resolver_cache.py
@@ -114,6 +114,7 @@ def _parallel_gate_op(gate, qubits):
         'NamedQubit': cirq.NamedQubit,
         'NamedQid': cirq.NamedQid,
         'NoIdentifierQubit': cirq.testing.NoIdentifierQubit,
+        'OpIdentifier': cirq.OpIdentifier,
         '_PauliX': cirq.ops.pauli_gates._PauliX,
         '_PauliY': cirq.ops.pauli_gates._PauliY,
         '_PauliZ': cirq.ops.pauli_gates._PauliZ,
diff --git a/cirq-core/cirq/protocols/json_test_data/OpIdentifier.json b/cirq-core/cirq/protocols/json_test_data/OpIdentifier.json
new file mode 100644
index 00000000000..d33b909367d
--- /dev/null
+++ b/cirq-core/cirq/protocols/json_test_data/OpIdentifier.json
@@ -0,0 +1,10 @@
+{
+  "cirq_type": "OpIdentifier",
+  "gate_type": "XPowGate",
+  "qubits": [
+    {
+      "cirq_type": "LineQubit",
+      "x": 1
+    }
+  ]
+}
\ No newline at end of file
diff --git a/cirq-core/cirq/protocols/json_test_data/OpIdentifier.repr b/cirq-core/cirq/protocols/json_test_data/OpIdentifier.repr
new file mode 100644
index 00000000000..6b991bb0b2c
--- /dev/null
+++ b/cirq-core/cirq/protocols/json_test_data/OpIdentifier.repr
@@ -0,0 +1,4 @@
+cirq.devices.noise_utils.OpIdentifier(
+    cirq.ops.common_gates.XPowGate,
+    cirq.LineQubit(1)
+)
\ No newline at end of file

From 55069a968eaa8ef2c55b40b509c4b670bda4036e Mon Sep 17 00:00:00 2001
From: Orion Martin <40585662+95-martin-orion@users.noreply.github.com>
Date: Tue, 7 Dec 2021 11:49:19 -0800
Subject: [PATCH 02/17] Add coverage for OpId

---
 cirq-core/cirq/devices/noise_utils_test.py | 32 +++++++++++++++++++++-
 1 file changed, 31 insertions(+), 1 deletion(-)

diff --git a/cirq-core/cirq/devices/noise_utils_test.py b/cirq-core/cirq/devices/noise_utils_test.py
index 61b5852ced4..d49d9e6508a 100644
--- a/cirq-core/cirq/devices/noise_utils_test.py
+++ b/cirq-core/cirq/devices/noise_utils_test.py
@@ -29,12 +29,42 @@
 )
 
 
-def test_op_id():
+def test_op_identifier():
     op_id = OpIdentifier(cirq.XPowGate)
     assert cirq.X(cirq.LineQubit(1)) in op_id
     assert cirq.Rx(rads=1) in op_id
 
 
+def test_op_identifier_subtypes():
+    gate_id = OpIdentifier(cirq.Gate)
+    xpow_id = OpIdentifier(cirq.XPowGate)
+    x_on_q0_id = OpIdentifier(cirq.XPowGate, cirq.LineQubit(0))
+    assert xpow_id.is_proper_subtype_of(gate_id)
+    assert x_on_q0_id.is_proper_subtype_of(xpow_id)
+    assert x_on_q0_id.is_proper_subtype_of(gate_id)
+    assert not xpow_id.is_proper_subtype_of(xpow_id)
+
+
+def test_op_id_str():
+    op_id = OpIdentifier(cirq.XPowGate, cirq.LineQubit(0))
+    print(op_id)
+    print(repr(op_id))
+    assert str(op_id) == "<class 'cirq.ops.common_gates.XPowGate'>(cirq.LineQubit(0),)"
+    assert repr(op_id) == (
+        "cirq.devices.noise_utils.OpIdentifier(cirq.ops.common_gates.XPowGate, cirq.LineQubit(0))"
+    )
+
+
+def test_op_id_swap():
+    q0, q1 = cirq.LineQubit.range(2)
+    base_id = OpIdentifier(cirq.CZPowGate, q0, q1)
+    swap_id = base_id.swapped()
+    assert cirq.CZ(q0, q1) in base_id
+    assert cirq.CZ(q0, q1) not in swap_id
+    assert cirq.CZ(q1, q0) not in base_id
+    assert cirq.CZ(q1, q0) in swap_id
+
+
 @pytest.mark.parametrize(
     'decay_constant,num_qubits,expected_output',
     [

From 75dc3e84ff76d5ddea6b47c626006ea243beaf06 Mon Sep 17 00:00:00 2001
From: Orion Martin <40585662+95-martin-orion@users.noreply.github.com>
Date: Tue, 7 Dec 2021 14:29:13 -0800
Subject: [PATCH 03/17] Clarify if-case

---
 cirq-core/cirq/devices/noise_utils.py | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/cirq-core/cirq/devices/noise_utils.py b/cirq-core/cirq/devices/noise_utils.py
index e02ed8baf45..35c84ce2a25 100644
--- a/cirq-core/cirq/devices/noise_utils.py
+++ b/cirq-core/cirq/devices/noise_utils.py
@@ -62,11 +62,12 @@ def is_proper_subtype_of(self, op_id: 'OpIdentifier'):
         more_specific_gate = self.gate_type != op_id.gate_type and issubclass(
             self.gate_type, op_id.gate_type
         )
-        return (
-            (more_specific_qubits or more_specific_gate)
-            and (more_specific_qubits or self.qubits == op_id.qubits)
-            and (more_specific_gate or self.gate_type == op_id.gate_type)
-        )
+        if more_specific_qubits:
+            return more_specific_gate or self.gate_type == op_id.gate_type
+        elif more_specific_gate:
+            return more_specific_qubits or self.qubits == op_id.qubits
+        else:
+            return False
 
     def __contains__(self, item: Union[ops.Gate, ops.Operation]) -> bool:
         if isinstance(item, ops.Gate):

From 4f9fdc0c8a156eee1a9cd665b32e6e19a65b5a2e Mon Sep 17 00:00:00 2001
From: Orion Martin <40585662+95-martin-orion@users.noreply.github.com>
Date: Mon, 13 Dec 2021 11:38:45 -0800
Subject: [PATCH 04/17] snake_case_gamma

---
 cirq-core/cirq/devices/noise_utils.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/cirq-core/cirq/devices/noise_utils.py b/cirq-core/cirq/devices/noise_utils.py
index 35c84ce2a25..fe59d94a9c1 100644
--- a/cirq-core/cirq/devices/noise_utils.py
+++ b/cirq-core/cirq/devices/noise_utils.py
@@ -223,10 +223,10 @@ def decoherence_pauli_error(t1_ns: float, tphi_ns: float, gate_time_ns: float) -
     Returns:
         Calculated Pauli error resulting from decoherence.
     """
-    Gamma2 = (1 / (2 * t1_ns)) + 1 / tphi_ns
+    gamma_2 = (1 / (2 * t1_ns)) + 1 / tphi_ns
 
     exp1 = np.exp(-gate_time_ns / t1_ns)
-    exp2 = np.exp(-gate_time_ns * Gamma2)
+    exp2 = np.exp(-gate_time_ns * gamma_2)
     px = 0.25 * (1 - exp1)
     py = px
     pz = 0.5 * (1 - exp2) - px

From 2baaaa7d41fda2bb1619b725e97befbf0b121fa2 Mon Sep 17 00:00:00 2001
From: Orion Martin <40585662+95-martin-orion@users.noreply.github.com>
Date: Fri, 12 Nov 2021 09:47:15 -0800
Subject: [PATCH 05/17] Add insertion noise model.

---
 cirq-core/cirq/devices/__init__.py            |   4 +
 .../cirq/devices/insertion_noise_model.py     |  75 +++++++++++++
 .../devices/insertion_noise_model_test.py     | 104 ++++++++++++++++++
 3 files changed, 183 insertions(+)
 create mode 100644 cirq-core/cirq/devices/insertion_noise_model.py
 create mode 100644 cirq-core/cirq/devices/insertion_noise_model_test.py

diff --git a/cirq-core/cirq/devices/__init__.py b/cirq-core/cirq/devices/__init__.py
index 72f793a7e46..89f5d1e22b5 100644
--- a/cirq-core/cirq/devices/__init__.py
+++ b/cirq-core/cirq/devices/__init__.py
@@ -48,6 +48,10 @@
     draw_placements,
 )
 
+from cirq.devices.insertion_noise_model import (
+    InsertionNoiseModel,
+)
+
 from cirq.devices.noise_utils import (
     OpIdentifier,
     decay_constant_to_xeb_fidelity,
diff --git a/cirq-core/cirq/devices/insertion_noise_model.py b/cirq-core/cirq/devices/insertion_noise_model.py
new file mode 100644
index 00000000000..6fc3d44b1cc
--- /dev/null
+++ b/cirq-core/cirq/devices/insertion_noise_model.py
@@ -0,0 +1,75 @@
+# Copyright 2021 The Cirq Developers
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     https://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from dataclasses import dataclass, field
+from typing import TYPE_CHECKING, Dict, List, Sequence
+
+from cirq import devices, ops
+from cirq.devices.noise_utils import (
+    OpIdentifier,
+    PHYSICAL_GATE_TAG,
+)
+
+if TYPE_CHECKING:
+    import cirq
+
+
+@dataclass
+class InsertionNoiseModel(devices.NoiseModel):
+    """Simple base noise model for inserting operations.
+
+    Operations generated by this model for a given moment are all added into a
+    single "noise moment", which is added before or after the original moment
+    based on `prepend`.
+
+    Args:
+        ops_added: a map of gate types (and optionally, qubits they act on) to
+            operations that should be added.
+        prepend: whether to add the new moment before the current one.
+        require_physical_tag: whether to only apply noise to operations tagged
+            with PHYSICAL_GATE_TAG.
+    """
+
+    ops_added: Dict[OpIdentifier, 'cirq.Operation'] = field(default_factory=dict)
+    prepend: bool = False
+    require_physical_tag: bool = True
+
+    def noisy_moment(
+        self, moment: 'cirq.Moment', system_qubits: Sequence['cirq.Qid']
+    ) -> 'cirq.OP_TREE':
+        noise_ops: List['cirq.Operation'] = []
+        for op in moment:
+            if self.require_physical_tag and PHYSICAL_GATE_TAG not in op.tags:
+                # Only non-virtual gates get noise applied.
+                continue
+            op_id = OpIdentifier(type(op.gate), *op.qubits)
+            if op_id in self.ops_added:
+                noise_ops.append(self.ops_added[op_id])
+                continue
+            # Find the closest match, if one exists.
+            parent_id = OpIdentifier(object, *op.qubits)
+            for added_id in self.ops_added:
+                if added_id.qubits != parent_id.qubits:
+                    continue
+                if not issubclass(op_id.gate_type, added_id.gate_type):
+                    continue
+                if issubclass(added_id.gate_type, parent_id.gate_type):
+                    parent_id = added_id
+            if parent_id.gate_type != object:
+                noise_ops.append(self.ops_added[parent_id])
+        if not noise_ops:
+            return [moment]
+        if self.prepend:
+            return [ops.Moment(noise_ops), moment]
+        return [moment, ops.Moment(noise_ops)]
diff --git a/cirq-core/cirq/devices/insertion_noise_model_test.py b/cirq-core/cirq/devices/insertion_noise_model_test.py
new file mode 100644
index 00000000000..cd11de19b18
--- /dev/null
+++ b/cirq-core/cirq/devices/insertion_noise_model_test.py
@@ -0,0 +1,104 @@
+# Copyright 2021 The Cirq Developers
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     https://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import cirq
+from cirq.devices.insertion_noise_model import InsertionNoiseModel
+from cirq.devices.noise_utils import (
+    PHYSICAL_GATE_TAG,
+    OpIdentifier,
+)
+
+
+def test_insertion_noise():
+    q0, q1 = cirq.LineQubit.range(2)
+    op_id0 = OpIdentifier(cirq.XPowGate, q0)
+    op_id1 = OpIdentifier(cirq.PhasedXZGate, q1)
+    model = InsertionNoiseModel(
+        {op_id0: cirq.T(q0), op_id1: cirq.H(q1)}, require_physical_tag=False
+    )
+    assert not model.prepend
+
+    phased_xz = cirq.PhasedXZGate(x_exponent=1.0, z_exponent=0.5, axis_phase_exponent=0.25)
+    moment_0 = cirq.Moment(cirq.X(q0), cirq.X(q1))
+    assert model.noisy_moment(moment_0, system_qubits=[q0, q1]) == [
+        moment_0,
+        cirq.Moment(cirq.T(q0)),
+    ]
+
+    moment_1 = cirq.Moment(phased_xz.on(q0), phased_xz.on(q1))
+    assert model.noisy_moment(moment_1, system_qubits=[q0, q1]) == [
+        moment_1,
+        cirq.Moment(cirq.H(q1)),
+    ]
+
+    moment_2 = cirq.Moment(cirq.X(q0), phased_xz.on(q1))
+    assert model.noisy_moment(moment_2, system_qubits=[q0, q1]) == [
+        moment_2,
+        cirq.Moment(cirq.T(q0), cirq.H(q1)),
+    ]
+
+    moment_3 = cirq.Moment(phased_xz.on(q0), cirq.X(q1))
+    assert model.noisy_moment(moment_3, system_qubits=[q0, q1]) == [moment_3]
+
+
+def test_prepend():
+    q0, q1 = cirq.LineQubit.range(2)
+    op_id0 = OpIdentifier(cirq.XPowGate, q0)
+    op_id1 = OpIdentifier(cirq.ZPowGate, q1)
+    model = InsertionNoiseModel(
+        {op_id0: cirq.T(q0), op_id1: cirq.H(q1)}, prepend=True, require_physical_tag=False
+    )
+
+    moment_0 = cirq.Moment(cirq.X(q0), cirq.Z(q1))
+    assert model.noisy_moment(moment_0, system_qubits=[q0, q1]) == [
+        cirq.Moment(cirq.T(q0), cirq.H(q1)),
+        moment_0,
+    ]
+
+
+def test_require_physical_tag():
+    q0, q1 = cirq.LineQubit.range(2)
+    op_id0 = OpIdentifier(cirq.XPowGate, q0)
+    op_id1 = OpIdentifier(cirq.ZPowGate, q1)
+    model = InsertionNoiseModel({op_id0: cirq.T(q0), op_id1: cirq.H(q1)})
+    assert model.require_physical_tag
+
+    moment_0 = cirq.Moment(cirq.X(q0).with_tags(PHYSICAL_GATE_TAG), cirq.Z(q1))
+    assert model.noisy_moment(moment_0, system_qubits=[q0, q1]) == [
+        moment_0,
+        cirq.Moment(cirq.T(q0)),
+    ]
+
+
+def test_supertype_matching():
+    # Demonstrate that the model applies the closest matching type
+    # if multiple types match a given gate.
+    q0 = cirq.LineQubit(0)
+    op_id0 = OpIdentifier(cirq.Gate, q0)
+    op_id1 = OpIdentifier(cirq.XPowGate, q0)
+    model = InsertionNoiseModel(
+        {op_id0: cirq.T(q0), op_id1: cirq.S(q0)}, require_physical_tag=False
+    )
+
+    moment_0 = cirq.Moment(cirq.Rx(rads=1).on(q0))
+    assert model.noisy_moment(moment_0, system_qubits=[q0]) == [
+        moment_0,
+        cirq.Moment(cirq.S(q0)),
+    ]
+
+    moment_1 = cirq.Moment(cirq.Y(q0))
+    assert model.noisy_moment(moment_1, system_qubits=[q0]) == [
+        moment_1,
+        cirq.Moment(cirq.T(q0)),
+    ]

From 224ee36c8febac28a33a1b8b9ad1f70706adf75c Mon Sep 17 00:00:00 2001
From: Orion Martin <40585662+95-martin-orion@users.noreply.github.com>
Date: Mon, 13 Dec 2021 14:09:50 -0800
Subject: [PATCH 06/17] Update to use new match methods

---
 .../cirq/devices/insertion_noise_model.py     | 24 ++++++++-----------
 .../devices/insertion_noise_model_test.py     |  9 ++++---
 2 files changed, 14 insertions(+), 19 deletions(-)

diff --git a/cirq-core/cirq/devices/insertion_noise_model.py b/cirq-core/cirq/devices/insertion_noise_model.py
index 6fc3d44b1cc..ed8409c3ca6 100644
--- a/cirq-core/cirq/devices/insertion_noise_model.py
+++ b/cirq-core/cirq/devices/insertion_noise_model.py
@@ -13,7 +13,7 @@
 # limitations under the License.
 
 from dataclasses import dataclass, field
-from typing import TYPE_CHECKING, Dict, List, Sequence
+from typing import TYPE_CHECKING, Dict, List, Optional, Sequence
 
 from cirq import devices, ops
 from cirq.devices.noise_utils import (
@@ -53,21 +53,17 @@ def noisy_moment(
             if self.require_physical_tag and PHYSICAL_GATE_TAG not in op.tags:
                 # Only non-virtual gates get noise applied.
                 continue
-            op_id = OpIdentifier(type(op.gate), *op.qubits)
-            if op_id in self.ops_added:
-                noise_ops.append(self.ops_added[op_id])
-                continue
-            # Find the closest match, if one exists.
-            parent_id = OpIdentifier(object, *op.qubits)
-            for added_id in self.ops_added:
-                if added_id.qubits != parent_id.qubits:
+            match_id: Optional[OpIdentifier] = None
+            for op_id in self.ops_added:
+                if op not in op_id:
                     continue
-                if not issubclass(op_id.gate_type, added_id.gate_type):
+                elif match_id is None:
+                    match_id = op_id
                     continue
-                if issubclass(added_id.gate_type, parent_id.gate_type):
-                    parent_id = added_id
-            if parent_id.gate_type != object:
-                noise_ops.append(self.ops_added[parent_id])
+                elif op_id.is_proper_subtype_of(match_id):
+                    match_id = op_id
+            if match_id is not None:
+                noise_ops.append(self.ops_added[match_id])
         if not noise_ops:
             return [moment]
         if self.prepend:
diff --git a/cirq-core/cirq/devices/insertion_noise_model_test.py b/cirq-core/cirq/devices/insertion_noise_model_test.py
index cd11de19b18..15216def58c 100644
--- a/cirq-core/cirq/devices/insertion_noise_model_test.py
+++ b/cirq-core/cirq/devices/insertion_noise_model_test.py
@@ -23,32 +23,31 @@
 def test_insertion_noise():
     q0, q1 = cirq.LineQubit.range(2)
     op_id0 = OpIdentifier(cirq.XPowGate, q0)
-    op_id1 = OpIdentifier(cirq.PhasedXZGate, q1)
+    op_id1 = OpIdentifier(cirq.ZPowGate, q1)
     model = InsertionNoiseModel(
         {op_id0: cirq.T(q0), op_id1: cirq.H(q1)}, require_physical_tag=False
     )
     assert not model.prepend
 
-    phased_xz = cirq.PhasedXZGate(x_exponent=1.0, z_exponent=0.5, axis_phase_exponent=0.25)
     moment_0 = cirq.Moment(cirq.X(q0), cirq.X(q1))
     assert model.noisy_moment(moment_0, system_qubits=[q0, q1]) == [
         moment_0,
         cirq.Moment(cirq.T(q0)),
     ]
 
-    moment_1 = cirq.Moment(phased_xz.on(q0), phased_xz.on(q1))
+    moment_1 = cirq.Moment(cirq.Z(q0), cirq.Z(q1))
     assert model.noisy_moment(moment_1, system_qubits=[q0, q1]) == [
         moment_1,
         cirq.Moment(cirq.H(q1)),
     ]
 
-    moment_2 = cirq.Moment(cirq.X(q0), phased_xz.on(q1))
+    moment_2 = cirq.Moment(cirq.X(q0), cirq.Z(q1))
     assert model.noisy_moment(moment_2, system_qubits=[q0, q1]) == [
         moment_2,
         cirq.Moment(cirq.T(q0), cirq.H(q1)),
     ]
 
-    moment_3 = cirq.Moment(phased_xz.on(q0), cirq.X(q1))
+    moment_3 = cirq.Moment(cirq.Z(q0), cirq.X(q1))
     assert model.noisy_moment(moment_3, system_qubits=[q0, q1]) == [moment_3]
 
 

From 4d460c20381e443d5a25f7868a501e9f30002f38 Mon Sep 17 00:00:00 2001
From: Orion Martin <40585662+95-martin-orion@users.noreply.github.com>
Date: Tue, 14 Dec 2021 13:24:30 -0800
Subject: [PATCH 07/17] Apply review cleanup

---
 .../cirq/devices/insertion_noise_model.py     | 35 +++++++++----------
 1 file changed, 16 insertions(+), 19 deletions(-)

diff --git a/cirq-core/cirq/devices/insertion_noise_model.py b/cirq-core/cirq/devices/insertion_noise_model.py
index ed8409c3ca6..6d6cf66a068 100644
--- a/cirq-core/cirq/devices/insertion_noise_model.py
+++ b/cirq-core/cirq/devices/insertion_noise_model.py
@@ -12,20 +12,17 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-from dataclasses import dataclass, field
+import dataclasses
 from typing import TYPE_CHECKING, Dict, List, Optional, Sequence
 
 from cirq import devices, ops
-from cirq.devices.noise_utils import (
-    OpIdentifier,
-    PHYSICAL_GATE_TAG,
-)
+from cirq.devices import noise_utils
 
 if TYPE_CHECKING:
     import cirq
 
 
-@dataclass
+@dataclasses.dataclass
 class InsertionNoiseModel(devices.NoiseModel):
     """Simple base noise model for inserting operations.
 
@@ -41,7 +38,9 @@ class InsertionNoiseModel(devices.NoiseModel):
             with PHYSICAL_GATE_TAG.
     """
 
-    ops_added: Dict[OpIdentifier, 'cirq.Operation'] = field(default_factory=dict)
+    ops_added: Dict[noise_utils.OpIdentifier, 'cirq.Operation'] = dataclasses.field(
+        default_factory=dict
+    )
     prepend: bool = False
     require_physical_tag: bool = True
 
@@ -49,18 +48,16 @@ def noisy_moment(
         self, moment: 'cirq.Moment', system_qubits: Sequence['cirq.Qid']
     ) -> 'cirq.OP_TREE':
         noise_ops: List['cirq.Operation'] = []
-        for op in moment:
-            if self.require_physical_tag and PHYSICAL_GATE_TAG not in op.tags:
-                # Only non-virtual gates get noise applied.
-                continue
-            match_id: Optional[OpIdentifier] = None
-            for op_id in self.ops_added:
-                if op not in op_id:
-                    continue
-                elif match_id is None:
-                    match_id = op_id
-                    continue
-                elif op_id.is_proper_subtype_of(match_id):
+        candidate_ops = [
+            op
+            for op in moment
+            if (not self.require_physical_tag) or noise_utils.PHYSICAL_GATE_TAG in op.tags
+        ]
+        for op in candidate_ops:
+            match_id: Optional[noise_utils.OpIdentifier] = None
+            candidate_ids = [op_id for op_id in self.ops_added if op in op_id]
+            for op_id in candidate_ids:
+                if match_id is None or op_id.is_proper_subtype_of(match_id):
                     match_id = op_id
             if match_id is not None:
                 noise_ops.append(self.ops_added[match_id])

From ee546c5c56b1ea27181e10b283c226ad6811ab60 Mon Sep 17 00:00:00 2001
From: Orion Martin <40585662+95-martin-orion@users.noreply.github.com>
Date: Thu, 11 Nov 2021 08:08:39 -0800
Subject: [PATCH 08/17] Introduce Thermal noise model.

---
 cirq-core/cirq/devices/__init__.py            |   9 +
 cirq-core/cirq/devices/thermal_noise_model.py | 237 ++++++++++++++++++
 .../cirq/devices/thermal_noise_model_test.py  | 237 ++++++++++++++++++
 3 files changed, 483 insertions(+)
 create mode 100644 cirq-core/cirq/devices/thermal_noise_model.py
 create mode 100644 cirq-core/cirq/devices/thermal_noise_model_test.py

diff --git a/cirq-core/cirq/devices/__init__.py b/cirq-core/cirq/devices/__init__.py
index 89f5d1e22b5..7c303725331 100644
--- a/cirq-core/cirq/devices/__init__.py
+++ b/cirq-core/cirq/devices/__init__.py
@@ -39,6 +39,11 @@
     ConstantQubitNoiseModel,
 )
 
+from cirq.devices.noise_properties import (
+    NoiseProperties,
+    NoiseModelFromNoiseProperties,
+)
+
 from cirq.devices.named_topologies import (
     NamedTopology,
     draw_gridlike,
@@ -52,6 +57,10 @@
     InsertionNoiseModel,
 )
 
+from cirq.devices.thermal_noise_model import (
+    ThermalNoiseModel,
+)
+
 from cirq.devices.noise_utils import (
     OpIdentifier,
     decay_constant_to_xeb_fidelity,
diff --git a/cirq-core/cirq/devices/thermal_noise_model.py b/cirq-core/cirq/devices/thermal_noise_model.py
new file mode 100644
index 00000000000..016d4286f8a
--- /dev/null
+++ b/cirq-core/cirq/devices/thermal_noise_model.py
@@ -0,0 +1,237 @@
+from dataclasses import dataclass
+from typing import TYPE_CHECKING, Dict, List, Optional, Sequence, Union
+from scipy.linalg import expm
+import numpy as np
+
+from cirq import devices, ops, protocols, qis
+from cirq.devices.noise_utils import (
+    PHYSICAL_GATE_TAG,
+)
+
+if TYPE_CHECKING:
+    import cirq
+
+
+def _left_mul(mat: np.ndarray) -> np.ndarray:
+    """Superoperator associated with left multiplication by a square matrix."""
+    mat = np.asarray(mat)
+    if mat.shape[-1] != mat.shape[-2]:
+        raise ValueError(
+            f'_left_mul only accepts square matrices, but input matrix has shape {mat.shape}.'
+        )
+    dim = mat.shape[-1]
+
+    return np.kron(mat, np.eye(dim))
+
+
+def _right_mul(mat: np.ndarray) -> np.ndarray:
+    """Superoperator associated with right multiplication by a square matrix."""
+    mat = np.asarray(mat)
+    if mat.shape[-1] != mat.shape[-2]:
+        raise ValueError(
+            f'_right_mul only accepts square matrices, but input matrix has shape {mat.shape}.'
+        )
+    dim = mat.shape[-1]
+
+    return np.kron(np.eye(dim), np.swapaxes(mat, -2, -1))
+
+
+def _lindbladian(left_op: np.ndarray) -> np.ndarray:
+    r"""Superoperator representing a Lindbladian.
+
+    The Lindbladian generated by a single operator A is the superoperator
+
+        $$
+        L(\rho) = A \rho A^\dagger - 0.5 (A^\dagger A \rho + \rho A^\dagger A)
+        $$
+
+    Args:
+        left_op: The operator acting on the left in the Lindbladian (A above).
+
+    Returns:
+        Superoperator corresponding to the Lindbladian.
+    """
+    left_op = np.asarray(left_op)
+    right_op = left_op.conj().T
+    square = right_op @ left_op
+    out = _left_mul(left_op) @ _right_mul(right_op)
+    out -= 0.5 * (_left_mul(square) + _right_mul(square))
+    return out
+
+
+def _decoherence_matrix(
+    cool_rate: float, dephase_rate: float, heat_rate: float = 0.0, dim: int = 2
+) -> np.ndarray:
+    """Construct a rate matrix associated with decay and dephasing.
+
+    The units of the matrix match the units of the rates specified.
+    This matrix can be used to construct a ThermalChannel after rescaling
+    by an idling time (to make it dimensionless).
+
+    Args:
+        cool_rate: Decay rate of the system, usually 1 / T_1
+        dephase_rate: Static dephasing rate of the system, usually 1 / T_phi
+        heat_rate: Heating rate of the system (default 0).
+        dim: Number of energy levels to include (default 2).
+
+    Returns:
+        np.ndarray rate matrix for decay and dephasing.
+    """
+    # Heating (related to a^dag)
+    rate_matrix = np.diag(np.arange(1, dim) * heat_rate, 1).T.astype(float)
+    # Cooling (related to a)
+    rate_matrix += np.diag(np.arange(1, dim) * cool_rate, 1)
+    # Dephasing (related to n=a^dag * a)
+    # We specify i^2 since we take square root to get the Lindblad op later.
+    rate_matrix += np.diag(dephase_rate * np.arange(dim) ** 2)
+    return rate_matrix
+
+
+def _validate_rates(qubit_dims: Dict['cirq.Qid', int], rates: Dict['cirq.Qid', np.ndarray]) -> None:
+    """Check all rate matrices are square and of appropriate dimension.
+
+    We check rates are positive in the class validator.
+    """
+    if set(qubit_dims) != set(rates):
+        raise ValueError('qubits for rates inconsistent with those through qubit_dims')
+    for q in rates:
+        if rates[q].shape != (qubit_dims[q], qubit_dims[q]):
+            raise ValueError(
+                f'Invalid shape for rate matrix: should be ({qubit_dims[q]}, {qubit_dims[q]}), '
+                f'but got {rates[q].shape}'
+            )
+
+
+@dataclass
+class ThermalNoiseModel(devices.NoiseModel):
+    """NoiseModel representing simulated thermalization of a qubit.
+
+    This model is designed for qubits which use energy levels as their states.
+    "Heating" and "cooling" here are used to refer to environmental noise which
+    transitions a qubit to higher or lower energy levels, respectively.
+    """
+
+    def __init__(
+        self,
+        qubit_dims: Dict['cirq.Qid', int],
+        gate_durations_ns: Dict[type, float],
+        heat_rate_GHz: Union[float, Dict['cirq.Qid', float], None] = None,
+        cool_rate_GHz: Union[float, Dict['cirq.Qid', float], None] = None,
+        dephase_rate_GHz: Union[float, Dict['cirq.Qid', float], None] = None,
+        require_physical_tag: bool = True,
+        skip_measurements: bool = True,
+    ):
+        """Construct a ThermalNoiseModel data object.
+
+        Required Args:
+            qubit_dims: Dimension for all qubits in the system.
+                        Currently only supports dimension=2 (qubits, not qudits)
+        Optional Args:
+            heat_rate_GHz: single number (units GHz) specifying heating rate,
+                        either per qubit, or global value for all.
+                        Given a rate gh, the Lindblad op will be sqrt(gh)*a^dag
+                        (where a is annihilation),
+                        so that the heating Lindbldian is
+                        gh(a^dag • a - 0.5{a*a^dag, •}).
+            cool_rate_GHz: single number (units GHz) specifying cooling rate,
+                        either per qubit, or global value for all.
+                        Given a rate gc, the Lindblad op will be sqrt(gc)*a
+                        so that the cooling Lindbldian is
+                        gc(a • a^dag - 0.5{n, •})
+                        This number is equivalent to 1/T1.
+            dephase_rate_GHz: single number (units GHz) specifying dephasing
+                        rate, either per qubit, or global value for all.
+                        Given a rate gd, Lindblad op will be sqrt(2*gd)*n where
+                        n = a^dag * a, so that the dephasing Lindbldian is
+                        2 * gd * (n • n - 0.5{n^2, •}).
+                        This number is equivalent to 1/Tphi.
+            require_physical_tag: whether to only apply noise to operations
+                        tagged with PHYSICAL_GATE_TAG.
+            skip_measurements: whether to skip applying noise to measurements.
+
+        Returns:
+            The ThermalNoiseModel with specified parameters.
+        """
+        qubits = set(qubit_dims)
+        rate_dict = {}
+
+        def _as_rate_dict(
+            rate_or_dict: Optional[Union[float, Dict['cirq.Qid', float]]]
+        ) -> Dict['cirq.Qid', float]:
+            # Convert float or None input into dictionary form. Make sure no
+            # qubits are missing from dictionary input.
+            if rate_or_dict is None:
+                return {qb: 0.0 for qb in qubits}
+            elif isinstance(rate_or_dict, dict):
+                out = rate_or_dict.copy()
+                for qb in qubits:
+                    if qb not in rate_or_dict:
+                        out[qb] = 0.0
+                return out
+            else:
+                return {qb: rate_or_dict for qb in qubits}
+
+        heat_rate_GHz = _as_rate_dict(heat_rate_GHz)
+        cool_rate_GHz = _as_rate_dict(cool_rate_GHz)
+        dephase_rate_GHz = _as_rate_dict(dephase_rate_GHz)
+
+        for q, dim in qubit_dims.items():
+            gamma_h = heat_rate_GHz[q]
+            gamma_c = cool_rate_GHz[q]
+            gamma_phi = dephase_rate_GHz[q]
+
+            rate_dict[q] = _decoherence_matrix(gamma_c, gamma_phi, gamma_h, dim)
+
+        _validate_rates(qubit_dims, rate_dict)
+        self.gate_durations_ns: Dict[type, float] = gate_durations_ns
+        self.rate_matrix_GHz: Dict['cirq.Qid', np.ndarray] = rate_dict
+        self.require_physical_tag: bool = require_physical_tag
+        self.skip_measurements: bool = skip_measurements
+
+    def noisy_moment(
+        self, moment: 'cirq.Moment', system_qubits: Sequence['cirq.Qid']
+    ) -> 'cirq.OP_TREE':
+        noise_ops: List['cirq.Operation'] = []
+        moment_ns: float = 0
+        for op in moment:
+            op_duration: Optional[float] = None
+            for key, duration in self.gate_durations_ns.items():
+                if not issubclass(type(op.gate), key):
+                    continue  # gate type doesn't match
+                # TODO: remove assumption of same time across qubits
+                # if len(key) > 1 and op_data[:1] != key[:1]:
+                #     continue  # qubits don't match
+                op_duration = duration
+                break
+            if op_duration is None:
+                if not isinstance(op.gate, ops.WaitGate):
+                    continue
+                # special case for wait gates if not predefined
+                op_duration = op.gate.duration.total_nanos()
+            moment_ns = max(moment_ns, op_duration)
+
+        if moment_ns == 0:
+            return [moment]
+
+        for qubit in system_qubits:
+            qubit_op = moment.operation_at(qubit)
+            if qubit_op is None:
+                continue
+            if self.skip_measurements and protocols.is_measurement(qubit_op):
+                continue
+            if self.require_physical_tag and PHYSICAL_GATE_TAG not in qubit_op.tags:
+                # Only non-virtual gates get noise applied.
+                continue
+            rates = self.rate_matrix_GHz[qubit] * moment_ns
+            num_op = np.diag(np.sqrt(np.diag(rates)))
+            annihilation = np.sqrt(np.triu(rates, 1))
+            creation = np.sqrt(np.triu(rates.T, 1)).T
+            # Lindbladian with three Lindblad ops for the three processes
+            # Note: 'time' parameter already specified implicitly through rates
+            L = _lindbladian(annihilation) + _lindbladian(creation) + 2 * _lindbladian(num_op)
+            superop = expm(L.real)
+            kraus_ops = qis.superoperator_to_kraus(superop)
+            noise_ops.append(ops.KrausChannel(kraus_ops).on(qubit))
+        if not noise_ops:
+            return [moment]
+        return [moment, ops.Moment(noise_ops)]
diff --git a/cirq-core/cirq/devices/thermal_noise_model_test.py b/cirq-core/cirq/devices/thermal_noise_model_test.py
new file mode 100644
index 00000000000..bf1a9b1d352
--- /dev/null
+++ b/cirq-core/cirq/devices/thermal_noise_model_test.py
@@ -0,0 +1,237 @@
+import numpy as np
+import pytest
+
+import cirq
+from cirq.devices.noise_utils import (
+    PHYSICAL_GATE_TAG,
+)
+from cirq.devices.thermal_noise_model import (
+    _left_mul,
+    _right_mul,
+    _validate_rates,
+    ThermalNoiseModel,
+)
+
+
+def test_helper_method_errors():
+    with pytest.raises(ValueError, match='_left_mul only accepts square matrices'):
+        _ = _left_mul(np.array([[1, 2, 3], [4, 5, 6]]))
+
+    with pytest.raises(ValueError, match='_right_mul only accepts square matrices'):
+        _ = _right_mul(np.array([[1, 2, 3], [4, 5, 6]]))
+
+    q0, q1 = cirq.LineQubit.range(2)
+    with pytest.raises(ValueError, match='qubits for rates inconsistent'):
+        _validate_rates({q0: 2, q1: 2}, {q0: np.array([[0.01, 0.1], [0.02, 0.2]])})
+
+    with pytest.raises(ValueError, match='qubits for rates inconsistent'):
+        _validate_rates(
+            {q0: 2},
+            {q0: np.array([[0.01, 0.1], [0.02, 0.2]]), q1: np.array([[0.03, 0.3], [0.04, 0.4]])},
+        )
+
+    with pytest.raises(ValueError, match='Invalid shape for rate matrix'):
+        _validate_rates({q0: 2}, {q0: np.array([[0.001, 0.01, 0.1], [0.002, 0.02, 0.2]])})
+
+
+def test_create_thermal_noise_per_qubit():
+    q0, q1 = cirq.LineQubit.range(2)
+    gate_durations = {cirq.PhasedXZGate: 25.0}
+    heat_rate_GHz = {q0: 1e-5, q1: 2e-5}
+    cool_rate_GHz = {q0: 1e-4, q1: 2e-4}
+    dephase_rate_GHz = {q0: 3e-4, q1: 4e-4}
+    model = ThermalNoiseModel(
+        qubit_dims={q0: 2, q1: 2},
+        gate_durations_ns=gate_durations,
+        heat_rate_GHz=heat_rate_GHz,
+        cool_rate_GHz=cool_rate_GHz,
+        dephase_rate_GHz=dephase_rate_GHz,
+    )
+    assert model.gate_durations_ns == gate_durations
+    assert model.require_physical_tag == True
+    assert model.skip_measurements == True
+    assert np.allclose(model.rate_matrix_GHz[q0], np.array([[0, 1e-4], [1e-5, 3e-4]]))
+    assert np.allclose(model.rate_matrix_GHz[q1], np.array([[0, 2e-4], [2e-5, 4e-4]]))
+
+
+def test_create_thermal_noise_mixed_type():
+    q0, q1 = cirq.LineQubit.range(2)
+    gate_durations = {cirq.PhasedXZGate: 25.0}
+    heat_rate_GHz = None
+    cool_rate_GHz = {q0: 1e-4, q1: 2e-4}
+    dephase_rate_GHz = 3e-4
+    model = ThermalNoiseModel(
+        qubit_dims={q0: 2, q1: 2},
+        gate_durations_ns=gate_durations,
+        heat_rate_GHz=heat_rate_GHz,
+        cool_rate_GHz=cool_rate_GHz,
+        dephase_rate_GHz=dephase_rate_GHz,
+    )
+    assert model.gate_durations_ns == gate_durations
+    assert model.require_physical_tag == True
+    assert model.skip_measurements == True
+    assert np.allclose(model.rate_matrix_GHz[q0], np.array([[0, 1e-4], [0, 3e-4]]))
+    assert np.allclose(model.rate_matrix_GHz[q1], np.array([[0, 2e-4], [0, 3e-4]]))
+
+
+def test_incomplete_rates():
+    q0, q1 = cirq.LineQubit.range(2)
+    gate_durations = {cirq.PhasedXZGate: 25.0}
+    heat_rate_GHz = {q1: 1e-5}
+    cool_rate_GHz = {q0: 1e-4}
+    model = ThermalNoiseModel(
+        qubit_dims={q0: 2, q1: 2},
+        gate_durations_ns=gate_durations,
+        heat_rate_GHz=heat_rate_GHz,
+        cool_rate_GHz=cool_rate_GHz,
+        dephase_rate_GHz=None,
+    )
+    assert model.gate_durations_ns == gate_durations
+    assert model.require_physical_tag == True
+    assert model.skip_measurements == True
+    assert np.allclose(model.rate_matrix_GHz[q0], np.array([[0, 1e-4], [0, 0]]))
+    assert np.allclose(model.rate_matrix_GHz[q1], np.array([[0, 0], [1e-5, 0]]))
+
+
+def test_noise_from_zero_duration():
+    # Verify that a moment with no duration has no noise.
+    q0, q1 = cirq.LineQubit.range(2)
+    gate_durations = {}
+    heat_rate_GHz = {q1: 1e-5}
+    cool_rate_GHz = {q0: 1e-4}
+    model = ThermalNoiseModel(
+        qubit_dims={q0: 2, q1: 2},
+        gate_durations_ns=gate_durations,
+        heat_rate_GHz=heat_rate_GHz,
+        cool_rate_GHz=cool_rate_GHz,
+        dephase_rate_GHz=None,
+        require_physical_tag=False,
+        skip_measurements=False,
+    )
+    moment = cirq.Moment(cirq.Z(q0), cirq.Z(q1))
+    assert model.noisy_moment(moment, system_qubits=[q0, q1]) == [moment]
+
+
+def test_noise_from_virtual_gates():
+    # Verify that a moment with only virtual gates has no noise if
+    # require_physical_tag is True.
+    q0, q1 = cirq.LineQubit.range(2)
+    gate_durations = {cirq.ZPowGate: 25.0}
+    heat_rate_GHz = {q1: 1e-5}
+    cool_rate_GHz = {q0: 1e-4}
+    model = ThermalNoiseModel(
+        qubit_dims={q0: 2, q1: 2},
+        gate_durations_ns=gate_durations,
+        heat_rate_GHz=heat_rate_GHz,
+        cool_rate_GHz=cool_rate_GHz,
+        dephase_rate_GHz=None,
+        require_physical_tag=True,
+        skip_measurements=False,
+    )
+    moment = cirq.Moment(cirq.Z(q0), cirq.Z(q1))
+    assert model.noisy_moment(moment, system_qubits=[q0, q1]) == [moment]
+
+    part_virtual_moment = cirq.Moment(cirq.Z(q0), cirq.Z(q1).with_tags(PHYSICAL_GATE_TAG))
+    assert len(model.noisy_moment(part_virtual_moment, system_qubits=[q0, q1])) == 2
+
+    model.require_physical_tag = False
+    assert len(model.noisy_moment(moment, system_qubits=[q0, q1])) == 2
+
+
+def test_noise_from_measurement():
+    # Verify that a moment with only measurement gates has no noise if
+    # skip_measurements is True.
+    q0, q1 = cirq.LineQubit.range(2)
+    gate_durations = {
+        cirq.ZPowGate: 25.0,
+        cirq.MeasurementGate: 4000.0,
+    }
+    heat_rate_GHz = {q1: 1e-5}
+    cool_rate_GHz = {q0: 1e-4}
+    model = ThermalNoiseModel(
+        qubit_dims={q0: 2, q1: 2},
+        gate_durations_ns=gate_durations,
+        heat_rate_GHz=heat_rate_GHz,
+        cool_rate_GHz=cool_rate_GHz,
+        dephase_rate_GHz=None,
+        require_physical_tag=False,
+        skip_measurements=True,
+    )
+    moment = cirq.Moment(cirq.measure(q0, q1, key='m'))
+    assert model.noisy_moment(moment, system_qubits=[q0, q1]) == [moment]
+
+    part_measure_moment = cirq.Moment(cirq.measure(q0, key='m'), cirq.Z(q1))
+    assert len(model.noisy_moment(part_measure_moment, system_qubits=[q0, q1])) == 2
+
+    model.skip_measurements = False
+    assert len(model.noisy_moment(moment, system_qubits=[q0, q1])) == 2
+
+
+def test_noisy_moment_one_qubit():
+    q0, q1 = cirq.LineQubit.range(2)
+    model = ThermalNoiseModel(
+        qubit_dims={q0: 2, q1: 2},
+        gate_durations_ns={
+            cirq.PhasedXZGate: 25.0,
+            cirq.CZPowGate: 25.0,
+        },
+        heat_rate_GHz={q0: 1e-5, q1: 2e-5},
+        cool_rate_GHz={q0: 1e-4, q1: 2e-4},
+        dephase_rate_GHz={q0: 3e-4, q1: 4e-4},
+        require_physical_tag=False,
+    )
+    gate = cirq.PhasedXZGate(x_exponent=1, z_exponent=0.5, axis_phase_exponent=0.25)
+    moment = cirq.Moment(gate.on(q0))
+    noisy_moment = model.noisy_moment(moment, system_qubits=[q0, q1])
+    assert noisy_moment[0] == moment
+    assert len(noisy_moment[1]) == 1
+    noisy_choi = cirq.kraus_to_choi(cirq.kraus(noisy_moment[1].operations[0]))
+    assert np.allclose(
+        noisy_choi,
+        [
+            [9.99750343e-01, 0, 0, 9.91164267e-01],
+            [0, 2.49656565e-03, 0, 0],
+            [0, 0, 2.49656565e-04, 0],
+            [9.91164267e-01, 0, 0, 9.97503434e-01],
+        ],
+    )
+
+
+def test_noisy_moment_two_qubit():
+    q0, q1 = cirq.LineQubit.range(2)
+    model = ThermalNoiseModel(
+        qubit_dims={q0: 2, q1: 2},
+        gate_durations_ns={
+            cirq.PhasedXZGate: 25.0,
+            cirq.CZPowGate: 25.0,
+        },
+        heat_rate_GHz={q0: 1e-5, q1: 2e-5},
+        cool_rate_GHz={q0: 1e-4, q1: 2e-4},
+        dephase_rate_GHz={q0: 3e-4, q1: 4e-4},
+        require_physical_tag=False,
+    )
+    gate = cirq.CZ ** 0.5
+    moment = cirq.Moment(gate.on(q0, q1))
+    noisy_moment = model.noisy_moment(moment, system_qubits=[q0, q1])
+    assert noisy_moment[0] == moment
+    assert len(noisy_moment[1]) == 2
+    noisy_choi_0 = cirq.kraus_to_choi(cirq.kraus(noisy_moment[1].operations[0]))
+    assert np.allclose(
+        noisy_choi_0,
+        [
+            [9.99750343e-01, 0, 0, 9.91164267e-01],
+            [0, 2.49656565e-03, 0, 0],
+            [0, 0, 2.49656565e-04, 0],
+            [9.91164267e-01, 0, 0, 9.97503434e-01],
+        ],
+    )
+    noisy_choi_1 = cirq.kraus_to_choi(cirq.kraus(noisy_moment[1].operations[1]))
+    assert np.allclose(
+        noisy_choi_1,
+        [
+            [9.99501372e-01, 0, 0, 9.87330937e-01],
+            [0, 4.98627517e-03, 0, 0],
+            [0, 0, 4.98627517e-04, 0],
+            [9.87330937e-01, 0, 0, 9.95013725e-01],
+        ],
+    )

From 2e1cad8fd68bcf38b788d81b487ff179f19f0dae Mon Sep 17 00:00:00 2001
From: Orion Martin <40585662+95-martin-orion@users.noreply.github.com>
Date: Thu, 11 Nov 2021 09:15:13 -0800
Subject: [PATCH 09/17] Align with checks

---
 cirq-core/cirq/devices/__init__.py            |  5 --
 cirq-core/cirq/devices/thermal_noise_model.py | 14 +++++
 .../cirq/devices/thermal_noise_model_test.py  | 58 +++++++++++++++++--
 3 files changed, 66 insertions(+), 11 deletions(-)

diff --git a/cirq-core/cirq/devices/__init__.py b/cirq-core/cirq/devices/__init__.py
index 7c303725331..858526365bf 100644
--- a/cirq-core/cirq/devices/__init__.py
+++ b/cirq-core/cirq/devices/__init__.py
@@ -39,11 +39,6 @@
     ConstantQubitNoiseModel,
 )
 
-from cirq.devices.noise_properties import (
-    NoiseProperties,
-    NoiseModelFromNoiseProperties,
-)
-
 from cirq.devices.named_topologies import (
     NamedTopology,
     draw_gridlike,
diff --git a/cirq-core/cirq/devices/thermal_noise_model.py b/cirq-core/cirq/devices/thermal_noise_model.py
index 016d4286f8a..0bdc135d703 100644
--- a/cirq-core/cirq/devices/thermal_noise_model.py
+++ b/cirq-core/cirq/devices/thermal_noise_model.py
@@ -1,3 +1,17 @@
+# Copyright 2021 The Cirq Developers
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     https://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
 from dataclasses import dataclass
 from typing import TYPE_CHECKING, Dict, List, Optional, Sequence, Union
 from scipy.linalg import expm
diff --git a/cirq-core/cirq/devices/thermal_noise_model_test.py b/cirq-core/cirq/devices/thermal_noise_model_test.py
index bf1a9b1d352..c7c738ff6b9 100644
--- a/cirq-core/cirq/devices/thermal_noise_model_test.py
+++ b/cirq-core/cirq/devices/thermal_noise_model_test.py
@@ -1,3 +1,17 @@
+# Copyright 2021 The Cirq Developers
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     https://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
 import numpy as np
 import pytest
 
@@ -48,8 +62,8 @@ def test_create_thermal_noise_per_qubit():
         dephase_rate_GHz=dephase_rate_GHz,
     )
     assert model.gate_durations_ns == gate_durations
-    assert model.require_physical_tag == True
-    assert model.skip_measurements == True
+    assert model.require_physical_tag
+    assert model.skip_measurements
     assert np.allclose(model.rate_matrix_GHz[q0], np.array([[0, 1e-4], [1e-5, 3e-4]]))
     assert np.allclose(model.rate_matrix_GHz[q1], np.array([[0, 2e-4], [2e-5, 4e-4]]))
 
@@ -68,8 +82,8 @@ def test_create_thermal_noise_mixed_type():
         dephase_rate_GHz=dephase_rate_GHz,
     )
     assert model.gate_durations_ns == gate_durations
-    assert model.require_physical_tag == True
-    assert model.skip_measurements == True
+    assert model.require_physical_tag
+    assert model.skip_measurements
     assert np.allclose(model.rate_matrix_GHz[q0], np.array([[0, 1e-4], [0, 3e-4]]))
     assert np.allclose(model.rate_matrix_GHz[q1], np.array([[0, 2e-4], [0, 3e-4]]))
 
@@ -87,8 +101,8 @@ def test_incomplete_rates():
         dephase_rate_GHz=None,
     )
     assert model.gate_durations_ns == gate_durations
-    assert model.require_physical_tag == True
-    assert model.skip_measurements == True
+    assert model.require_physical_tag
+    assert model.skip_measurements
     assert np.allclose(model.rate_matrix_GHz[q0], np.array([[0, 1e-4], [0, 0]]))
     assert np.allclose(model.rate_matrix_GHz[q1], np.array([[0, 0], [1e-5, 0]]))
 
@@ -197,6 +211,38 @@ def test_noisy_moment_one_qubit():
     )
 
 
+def test_noise_from_wait():
+    # Verify that wait-gate noise is duration-dependent.
+    q0 = cirq.LineQubit(0)
+    gate_durations = {cirq.ZPowGate: 25.0}
+    heat_rate_GHz = {q0: 1e-5}
+    cool_rate_GHz = {q0: 1e-4}
+    model = ThermalNoiseModel(
+        qubit_dims={q0: 2},
+        gate_durations_ns=gate_durations,
+        heat_rate_GHz=heat_rate_GHz,
+        cool_rate_GHz=cool_rate_GHz,
+        dephase_rate_GHz=None,
+        require_physical_tag=False,
+        skip_measurements=True,
+    )
+    moment = cirq.Moment(cirq.wait(q0, nanos=100))
+    noisy_moment = model.noisy_moment(moment, system_qubits=[q0])
+    assert noisy_moment[0] == moment
+    assert len(noisy_moment[1]) == 1
+    noisy_choi = cirq.kraus_to_choi(cirq.kraus(noisy_moment[1].operations[0]))
+    print(noisy_choi)
+    assert np.allclose(
+        noisy_choi,
+        [
+            [9.99005480e-01, 0, 0, 9.94515097e-01],
+            [0, 9.94520111e-03, 0, 0],
+            [0, 0, 9.94520111e-04, 0],
+            [9.94515097e-01, 0, 0, 9.90054799e-01],
+        ],
+    )
+
+
 def test_noisy_moment_two_qubit():
     q0, q1 = cirq.LineQubit.range(2)
     model = ThermalNoiseModel(

From 35939290f0dd21f52af0dc8f0c9a76936d2a06af Mon Sep 17 00:00:00 2001
From: Orion Martin <40585662+95-martin-orion@users.noreply.github.com>
Date: Fri, 12 Nov 2021 10:17:21 -0800
Subject: [PATCH 10/17] remove ThermalChannel ref

---
 cirq-core/cirq/devices/thermal_noise_model.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/cirq-core/cirq/devices/thermal_noise_model.py b/cirq-core/cirq/devices/thermal_noise_model.py
index 0bdc135d703..0406f2d3d4d 100644
--- a/cirq-core/cirq/devices/thermal_noise_model.py
+++ b/cirq-core/cirq/devices/thermal_noise_model.py
@@ -79,7 +79,7 @@ def _decoherence_matrix(
     """Construct a rate matrix associated with decay and dephasing.
 
     The units of the matrix match the units of the rates specified.
-    This matrix can be used to construct a ThermalChannel after rescaling
+    This matrix can be used to construct a noise channel after rescaling
     by an idling time (to make it dimensionless).
 
     Args:

From 5b6f2515a3ca8590ea274461f2b014600364235b Mon Sep 17 00:00:00 2001
From: Orion Martin <40585662+95-martin-orion@users.noreply.github.com>
Date: Wed, 15 Dec 2021 10:29:37 -0800
Subject: [PATCH 11/17] Answer review comments

---
 cirq-core/cirq/devices/__init__.py            |  1 -
 .../cirq/devices/insertion_noise_model.py     | 10 +++++--
 .../devices/insertion_noise_model_test.py     | 14 +++++++++
 cirq-core/cirq/devices/noise_utils.py         | 29 +------------------
 cirq-core/cirq/devices/noise_utils_test.py    | 16 +---------
 5 files changed, 23 insertions(+), 47 deletions(-)

diff --git a/cirq-core/cirq/devices/__init__.py b/cirq-core/cirq/devices/__init__.py
index 858526365bf..b1480a9513d 100644
--- a/cirq-core/cirq/devices/__init__.py
+++ b/cirq-core/cirq/devices/__init__.py
@@ -63,7 +63,6 @@
     pauli_error_to_decay_constant,
     xeb_fidelity_to_decay_constant,
     pauli_error_from_t1,
-    pauli_error_from_depolarization,
     average_error,
     decoherence_pauli_error,
 )
diff --git a/cirq-core/cirq/devices/insertion_noise_model.py b/cirq-core/cirq/devices/insertion_noise_model.py
index 6d6cf66a068..678e422424d 100644
--- a/cirq-core/cirq/devices/insertion_noise_model.py
+++ b/cirq-core/cirq/devices/insertion_noise_model.py
@@ -15,7 +15,7 @@
 import dataclasses
 from typing import TYPE_CHECKING, Dict, List, Optional, Sequence
 
-from cirq import devices, ops
+from cirq import devices
 from cirq.devices import noise_utils
 
 if TYPE_CHECKING:
@@ -63,6 +63,10 @@ def noisy_moment(
                 noise_ops.append(self.ops_added[match_id])
         if not noise_ops:
             return [moment]
+
+        from cirq import circuits
+
+        noise_steps = circuits.Circuit(noise_ops)
         if self.prepend:
-            return [ops.Moment(noise_ops), moment]
-        return [moment, ops.Moment(noise_ops)]
+            return [*noise_steps.moments, moment]
+        return [moment, *noise_steps.moments]
diff --git a/cirq-core/cirq/devices/insertion_noise_model_test.py b/cirq-core/cirq/devices/insertion_noise_model_test.py
index 15216def58c..fea42f315a5 100644
--- a/cirq-core/cirq/devices/insertion_noise_model_test.py
+++ b/cirq-core/cirq/devices/insertion_noise_model_test.py
@@ -51,6 +51,20 @@ def test_insertion_noise():
     assert model.noisy_moment(moment_3, system_qubits=[q0, q1]) == [moment_3]
 
 
+def test_colliding_noise_qubits():
+    # Check that noise affecting other qubits doesn't cause issues.
+    q0, q1, q2, q3 = cirq.LineQubit.range(4)
+    op_id0 = OpIdentifier(cirq.CZPowGate)
+    model = InsertionNoiseModel({op_id0: cirq.CNOT(q1, q2)}, require_physical_tag=False)
+
+    moment_0 = cirq.Moment(cirq.CZ(q0, q1), cirq.CZ(q2, q3))
+    assert model.noisy_moment(moment_0, system_qubits=[q0, q1, q2, q3]) == [
+        moment_0,
+        cirq.Moment(cirq.CNOT(q1, q2)),
+        cirq.Moment(cirq.CNOT(q1, q2)),
+    ]
+
+
 def test_prepend():
     q0, q1 = cirq.LineQubit.range(2)
     op_id0 = OpIdentifier(cirq.XPowGate, q0)
diff --git a/cirq-core/cirq/devices/noise_utils.py b/cirq-core/cirq/devices/noise_utils.py
index fe59d94a9c1..22e21aab89d 100644
--- a/cirq-core/cirq/devices/noise_utils.py
+++ b/cirq-core/cirq/devices/noise_utils.py
@@ -13,7 +13,6 @@
 # limitations under the License.
 
 from typing import TYPE_CHECKING, Any, Dict, Tuple, Type, Union
-import warnings
 import numpy as np
 
 from cirq import ops, protocols, value
@@ -47,9 +46,6 @@ def qubits(self) -> Tuple['cirq.Qid', ...]:
     def _predicate(self, *args, **kwargs):
         return self._gate_family._predicate(*args, **kwargs)
 
-    def swapped(self):
-        return OpIdentifier(self.gate_type, *self.qubits[::-1])
-
     def is_proper_subtype_of(self, op_id: 'OpIdentifier'):
         """Returns true if this is contained within op_id, but not equal to it.
 
@@ -175,29 +171,6 @@ def pauli_error_from_t1(t_ns: float, t1_ns: float) -> float:
     return (1 - np.exp(-t_ns / t2)) / 2 + (1 - np.exp(-t_ns / t1_ns)) / 4
 
 
-def pauli_error_from_depolarization(t_ns: float, t1_ns: float, pauli_error: float = 0) -> float:
-    """Calculates the amount of pauli error from depolarization.
-
-    This computes non-T1 error for a specific duration, `t`. If pauli error
-    from T1 decay is more than total pauli error, this returns zero; otherwise,
-    it returns the portion of pauli error not attributable to T1 error.
-
-    Args:
-        t_ns: The duration of the gate in ns.
-        t1_ns: The T1 decay constant in ns.
-        pauli_error: The total pauli error.
-
-    Returns:
-        Calculated Pauli error resulting from depolarization.
-    """
-    t1_pauli_error = pauli_error_from_t1(t_ns, t1_ns)
-    if pauli_error >= t1_pauli_error:
-        return pauli_error - t1_pauli_error
-
-    warnings.warn("Pauli error from T1 decay is greater than total Pauli error", RuntimeWarning)
-    return 0
-
-
 def average_error(decay_constant: float, num_qubits: int = 1) -> float:
     """Calculates the average error from the depolarization decay constant.
 
@@ -213,7 +186,7 @@ def average_error(decay_constant: float, num_qubits: int = 1) -> float:
 
 
 def decoherence_pauli_error(t1_ns: float, tphi_ns: float, gate_time_ns: float) -> float:
-    """The component of Pauli error caused by decoherence.
+    """The component of Pauli error caused by decoherence on a single qubit.
 
     Args:
         t1_ns: T1 time in nanoseconds.
diff --git a/cirq-core/cirq/devices/noise_utils_test.py b/cirq-core/cirq/devices/noise_utils_test.py
index d49d9e6508a..ac541018bed 100644
--- a/cirq-core/cirq/devices/noise_utils_test.py
+++ b/cirq-core/cirq/devices/noise_utils_test.py
@@ -23,7 +23,6 @@
     pauli_error_to_decay_constant,
     xeb_fidelity_to_decay_constant,
     pauli_error_from_t1,
-    pauli_error_from_depolarization,
     average_error,
     decoherence_pauli_error,
 )
@@ -58,7 +57,7 @@ def test_op_id_str():
 def test_op_id_swap():
     q0, q1 = cirq.LineQubit.range(2)
     base_id = OpIdentifier(cirq.CZPowGate, q0, q1)
-    swap_id = base_id.swapped()
+    swap_id = OpIdentifier(base_id.gate_type, *base_id.qubits[::-1])
     assert cirq.CZ(q0, q1) in base_id
     assert cirq.CZ(q0, q1) not in swap_id
     assert cirq.CZ(q1, q0) not in base_id
@@ -125,19 +124,6 @@ def test_pauli_error_from_t1(t, t1_ns, expected_output):
     assert val == expected_output
 
 
-@pytest.mark.parametrize(
-    't,t1_ns,pauli_error,expected_output',
-    [
-        (20, 1e5, 0.01, 0.01 - ((1 - np.exp(-20 / 2e5)) / 2 + (1 - np.exp(-20 / 1e5)) / 4)),
-        # In this case, the formula produces a negative result.
-        (4000, 1e4, 0.01, 0),
-    ],
-)
-def test_pauli_error_from_depolarization(t, t1_ns, pauli_error, expected_output):
-    val = pauli_error_from_depolarization(t, t1_ns, pauli_error)
-    assert val == expected_output
-
-
 @pytest.mark.parametrize(
     'decay_constant,num_qubits,expected_output',
     [

From e8e6db6418d216e9674596ede0ac13edd776d1e2 Mon Sep 17 00:00:00 2001
From: Orion Martin <40585662+95-martin-orion@users.noreply.github.com>
Date: Mon, 20 Dec 2021 10:45:53 -0800
Subject: [PATCH 12/17] Document mutli-type behavior.

---
 cirq-core/cirq/devices/insertion_noise_model.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/cirq-core/cirq/devices/insertion_noise_model.py b/cirq-core/cirq/devices/insertion_noise_model.py
index 678e422424d..33553216de5 100644
--- a/cirq-core/cirq/devices/insertion_noise_model.py
+++ b/cirq-core/cirq/devices/insertion_noise_model.py
@@ -32,7 +32,11 @@ class InsertionNoiseModel(devices.NoiseModel):
 
     Args:
         ops_added: a map of gate types (and optionally, qubits they act on) to
-            operations that should be added.
+            operations that should be added. If two gate types provided apply
+            to a target gate, the most specific type will match; if neither
+            type is more specific (e.g. A is a subtype of B, but B defines
+            qubits and A does not) then the first one appering in this dict
+            will match.
         prepend: whether to add the new moment before the current one.
         require_physical_tag: whether to only apply noise to operations tagged
             with PHYSICAL_GATE_TAG.

From dc88cb1e99f3867c212cba58a3a6e0171813788b Mon Sep 17 00:00:00 2001
From: Orion Martin <40585662+95-martin-orion@users.noreply.github.com>
Date: Mon, 20 Dec 2021 11:34:54 -0800
Subject: [PATCH 13/17] qubit_dims -> qubits

---
 cirq-core/cirq/devices/thermal_noise_model.py | 19 +++++++--------
 .../cirq/devices/thermal_noise_model_test.py  | 24 +++++++++----------
 2 files changed, 21 insertions(+), 22 deletions(-)

diff --git a/cirq-core/cirq/devices/thermal_noise_model.py b/cirq-core/cirq/devices/thermal_noise_model.py
index 0406f2d3d4d..7409a658cf8 100644
--- a/cirq-core/cirq/devices/thermal_noise_model.py
+++ b/cirq-core/cirq/devices/thermal_noise_model.py
@@ -13,7 +13,7 @@
 # limitations under the License.
 
 from dataclasses import dataclass
-from typing import TYPE_CHECKING, Dict, List, Optional, Sequence, Union
+from typing import TYPE_CHECKING, Dict, List, Optional, Sequence, Set, Union
 from scipy.linalg import expm
 import numpy as np
 
@@ -101,17 +101,17 @@ def _decoherence_matrix(
     return rate_matrix
 
 
-def _validate_rates(qubit_dims: Dict['cirq.Qid', int], rates: Dict['cirq.Qid', np.ndarray]) -> None:
+def _validate_rates(qubits: Set['cirq.Qid'], rates: Dict['cirq.Qid', np.ndarray]) -> None:
     """Check all rate matrices are square and of appropriate dimension.
 
     We check rates are positive in the class validator.
     """
-    if set(qubit_dims) != set(rates):
+    if qubits != set(rates):
         raise ValueError('qubits for rates inconsistent with those through qubit_dims')
     for q in rates:
-        if rates[q].shape != (qubit_dims[q], qubit_dims[q]):
+        if rates[q].shape != (q.dimension, q.dimension):
             raise ValueError(
-                f'Invalid shape for rate matrix: should be ({qubit_dims[q]}, {qubit_dims[q]}), '
+                f'Invalid shape for rate matrix: should be ({q.dimension}, {q.dimension}), '
                 f'but got {rates[q].shape}'
             )
 
@@ -127,7 +127,7 @@ class ThermalNoiseModel(devices.NoiseModel):
 
     def __init__(
         self,
-        qubit_dims: Dict['cirq.Qid', int],
+        qubits: Set['cirq.Qid'],
         gate_durations_ns: Dict[type, float],
         heat_rate_GHz: Union[float, Dict['cirq.Qid', float], None] = None,
         cool_rate_GHz: Union[float, Dict['cirq.Qid', float], None] = None,
@@ -166,7 +166,6 @@ def __init__(
         Returns:
             The ThermalNoiseModel with specified parameters.
         """
-        qubits = set(qubit_dims)
         rate_dict = {}
 
         def _as_rate_dict(
@@ -189,14 +188,14 @@ def _as_rate_dict(
         cool_rate_GHz = _as_rate_dict(cool_rate_GHz)
         dephase_rate_GHz = _as_rate_dict(dephase_rate_GHz)
 
-        for q, dim in qubit_dims.items():
+        for q in qubits:
             gamma_h = heat_rate_GHz[q]
             gamma_c = cool_rate_GHz[q]
             gamma_phi = dephase_rate_GHz[q]
 
-            rate_dict[q] = _decoherence_matrix(gamma_c, gamma_phi, gamma_h, dim)
+            rate_dict[q] = _decoherence_matrix(gamma_c, gamma_phi, gamma_h, q.dimension)
 
-        _validate_rates(qubit_dims, rate_dict)
+        _validate_rates(qubits, rate_dict)
         self.gate_durations_ns: Dict[type, float] = gate_durations_ns
         self.rate_matrix_GHz: Dict['cirq.Qid', np.ndarray] = rate_dict
         self.require_physical_tag: bool = require_physical_tag
diff --git a/cirq-core/cirq/devices/thermal_noise_model_test.py b/cirq-core/cirq/devices/thermal_noise_model_test.py
index c7c738ff6b9..e1c0705a45e 100644
--- a/cirq-core/cirq/devices/thermal_noise_model_test.py
+++ b/cirq-core/cirq/devices/thermal_noise_model_test.py
@@ -36,16 +36,16 @@ def test_helper_method_errors():
 
     q0, q1 = cirq.LineQubit.range(2)
     with pytest.raises(ValueError, match='qubits for rates inconsistent'):
-        _validate_rates({q0: 2, q1: 2}, {q0: np.array([[0.01, 0.1], [0.02, 0.2]])})
+        _validate_rates({q0, q1}, {q0: np.array([[0.01, 0.1], [0.02, 0.2]])})
 
     with pytest.raises(ValueError, match='qubits for rates inconsistent'):
         _validate_rates(
-            {q0: 2},
+            {q0},
             {q0: np.array([[0.01, 0.1], [0.02, 0.2]]), q1: np.array([[0.03, 0.3], [0.04, 0.4]])},
         )
 
     with pytest.raises(ValueError, match='Invalid shape for rate matrix'):
-        _validate_rates({q0: 2}, {q0: np.array([[0.001, 0.01, 0.1], [0.002, 0.02, 0.2]])})
+        _validate_rates({q0}, {q0: np.array([[0.001, 0.01, 0.1], [0.002, 0.02, 0.2]])})
 
 
 def test_create_thermal_noise_per_qubit():
@@ -55,7 +55,7 @@ def test_create_thermal_noise_per_qubit():
     cool_rate_GHz = {q0: 1e-4, q1: 2e-4}
     dephase_rate_GHz = {q0: 3e-4, q1: 4e-4}
     model = ThermalNoiseModel(
-        qubit_dims={q0: 2, q1: 2},
+        qubits={q0, q1},
         gate_durations_ns=gate_durations,
         heat_rate_GHz=heat_rate_GHz,
         cool_rate_GHz=cool_rate_GHz,
@@ -75,7 +75,7 @@ def test_create_thermal_noise_mixed_type():
     cool_rate_GHz = {q0: 1e-4, q1: 2e-4}
     dephase_rate_GHz = 3e-4
     model = ThermalNoiseModel(
-        qubit_dims={q0: 2, q1: 2},
+        qubits={q0, q1},
         gate_durations_ns=gate_durations,
         heat_rate_GHz=heat_rate_GHz,
         cool_rate_GHz=cool_rate_GHz,
@@ -94,7 +94,7 @@ def test_incomplete_rates():
     heat_rate_GHz = {q1: 1e-5}
     cool_rate_GHz = {q0: 1e-4}
     model = ThermalNoiseModel(
-        qubit_dims={q0: 2, q1: 2},
+        qubits={q0, q1},
         gate_durations_ns=gate_durations,
         heat_rate_GHz=heat_rate_GHz,
         cool_rate_GHz=cool_rate_GHz,
@@ -114,7 +114,7 @@ def test_noise_from_zero_duration():
     heat_rate_GHz = {q1: 1e-5}
     cool_rate_GHz = {q0: 1e-4}
     model = ThermalNoiseModel(
-        qubit_dims={q0: 2, q1: 2},
+        qubits={q0, q1},
         gate_durations_ns=gate_durations,
         heat_rate_GHz=heat_rate_GHz,
         cool_rate_GHz=cool_rate_GHz,
@@ -134,7 +134,7 @@ def test_noise_from_virtual_gates():
     heat_rate_GHz = {q1: 1e-5}
     cool_rate_GHz = {q0: 1e-4}
     model = ThermalNoiseModel(
-        qubit_dims={q0: 2, q1: 2},
+        qubits={q0, q1},
         gate_durations_ns=gate_durations,
         heat_rate_GHz=heat_rate_GHz,
         cool_rate_GHz=cool_rate_GHz,
@@ -163,7 +163,7 @@ def test_noise_from_measurement():
     heat_rate_GHz = {q1: 1e-5}
     cool_rate_GHz = {q0: 1e-4}
     model = ThermalNoiseModel(
-        qubit_dims={q0: 2, q1: 2},
+        qubits={q0, q1},
         gate_durations_ns=gate_durations,
         heat_rate_GHz=heat_rate_GHz,
         cool_rate_GHz=cool_rate_GHz,
@@ -184,7 +184,7 @@ def test_noise_from_measurement():
 def test_noisy_moment_one_qubit():
     q0, q1 = cirq.LineQubit.range(2)
     model = ThermalNoiseModel(
-        qubit_dims={q0: 2, q1: 2},
+        qubits={q0, q1},
         gate_durations_ns={
             cirq.PhasedXZGate: 25.0,
             cirq.CZPowGate: 25.0,
@@ -218,7 +218,7 @@ def test_noise_from_wait():
     heat_rate_GHz = {q0: 1e-5}
     cool_rate_GHz = {q0: 1e-4}
     model = ThermalNoiseModel(
-        qubit_dims={q0: 2},
+        qubits={q0},
         gate_durations_ns=gate_durations,
         heat_rate_GHz=heat_rate_GHz,
         cool_rate_GHz=cool_rate_GHz,
@@ -246,7 +246,7 @@ def test_noise_from_wait():
 def test_noisy_moment_two_qubit():
     q0, q1 = cirq.LineQubit.range(2)
     model = ThermalNoiseModel(
-        qubit_dims={q0: 2, q1: 2},
+        qubits={q0, q1},
         gate_durations_ns={
             cirq.PhasedXZGate: 25.0,
             cirq.CZPowGate: 25.0,

From af4ce30fcd11e992d78bdeec01ebaeed125de114 Mon Sep 17 00:00:00 2001
From: Orion Martin <40585662+95-martin-orion@users.noreply.github.com>
Date: Mon, 20 Dec 2021 11:52:24 -0800
Subject: [PATCH 14/17] Extract methods, clean docstrings

---
 cirq-core/cirq/devices/thermal_noise_model.py | 86 +++++++++----------
 1 file changed, 43 insertions(+), 43 deletions(-)

diff --git a/cirq-core/cirq/devices/thermal_noise_model.py b/cirq-core/cirq/devices/thermal_noise_model.py
index 7409a658cf8..0c279d28ff3 100644
--- a/cirq-core/cirq/devices/thermal_noise_model.py
+++ b/cirq-core/cirq/devices/thermal_noise_model.py
@@ -73,6 +73,17 @@ def _lindbladian(left_op: np.ndarray) -> np.ndarray:
     return out
 
 
+def _kraus_ops_from_rates(rates: np.ndarray) -> Sequence[np.ndarray]:
+    num_op = np.diag(np.sqrt(np.diag(rates)))
+    annihilation = np.sqrt(np.triu(rates, 1))
+    creation = np.sqrt(np.triu(rates.T, 1)).T
+    # Lindbladian with three Lindblad ops for the three processes
+    # Note: 'time' parameter already specified implicitly through rates
+    L = _lindbladian(annihilation) + _lindbladian(creation) + 2 * _lindbladian(num_op)
+    superop = expm(L.real)
+    return qis.superoperator_to_kraus(superop)
+
+
 def _decoherence_matrix(
     cool_rate: float, dephase_rate: float, heat_rate: float = 0.0, dim: int = 2
 ) -> np.ndarray:
@@ -101,6 +112,20 @@ def _decoherence_matrix(
     return rate_matrix
 
 
+def _as_rate_dict(
+    rate_or_dict: Optional[Union[float, Dict['cirq.Qid', float]]],
+    qubits: Set['cirq.Qid'],
+) -> Dict['cirq.Qid', float]:
+    # Convert float or None input into dictionary form. Make sure no
+    # qubits are missing from dictionary input.
+    if rate_or_dict is None:
+        return {q: 0.0 for q in qubits}
+    elif isinstance(rate_or_dict, dict):
+        return {**{q: 0.0 for q in qubits}, **rate_or_dict}
+    else:
+        return {q: rate_or_dict for q in qubits}
+
+
 def _validate_rates(qubits: Set['cirq.Qid'], rates: Dict['cirq.Qid', np.ndarray]) -> None:
     """Check all rate matrices are square and of appropriate dimension.
 
@@ -142,25 +167,23 @@ def __init__(
                         Currently only supports dimension=2 (qubits, not qudits)
         Optional Args:
             heat_rate_GHz: single number (units GHz) specifying heating rate,
-                        either per qubit, or global value for all.
-                        Given a rate gh, the Lindblad op will be sqrt(gh)*a^dag
-                        (where a is annihilation),
-                        so that the heating Lindbldian is
-                        gh(a^dag • a - 0.5{a*a^dag, •}).
+                either per qubit, or global value for all.
+                Given a rate gh, the Lindblad op will be sqrt(gh)*a^dag
+                (where a is annihilation), so that the heating Lindbldian is
+                gh(a^dag • a - 0.5{a*a^dag, •}).
             cool_rate_GHz: single number (units GHz) specifying cooling rate,
-                        either per qubit, or global value for all.
-                        Given a rate gc, the Lindblad op will be sqrt(gc)*a
-                        so that the cooling Lindbldian is
-                        gc(a • a^dag - 0.5{n, •})
-                        This number is equivalent to 1/T1.
+                either per qubit, or global value for all.
+                Given a rate gc, the Lindblad op will be sqrt(gc)*a
+                so that the cooling Lindbldian is gc(a • a^dag - 0.5{n, •})
+                This number is equivalent to 1/T1.
             dephase_rate_GHz: single number (units GHz) specifying dephasing
-                        rate, either per qubit, or global value for all.
-                        Given a rate gd, Lindblad op will be sqrt(2*gd)*n where
-                        n = a^dag * a, so that the dephasing Lindbldian is
-                        2 * gd * (n • n - 0.5{n^2, •}).
-                        This number is equivalent to 1/Tphi.
+                rate, either per qubit, or global value for all.
+                Given a rate gd, Lindblad op will be sqrt(2*gd)*n where
+                n = a^dag * a, so that the dephasing Lindbldian is
+                2 * gd * (n • n - 0.5{n^2, •}).
+                This number is equivalent to 1/Tphi.
             require_physical_tag: whether to only apply noise to operations
-                        tagged with PHYSICAL_GATE_TAG.
+                tagged with PHYSICAL_GATE_TAG.
             skip_measurements: whether to skip applying noise to measurements.
 
         Returns:
@@ -168,25 +191,9 @@ def __init__(
         """
         rate_dict = {}
 
-        def _as_rate_dict(
-            rate_or_dict: Optional[Union[float, Dict['cirq.Qid', float]]]
-        ) -> Dict['cirq.Qid', float]:
-            # Convert float or None input into dictionary form. Make sure no
-            # qubits are missing from dictionary input.
-            if rate_or_dict is None:
-                return {qb: 0.0 for qb in qubits}
-            elif isinstance(rate_or_dict, dict):
-                out = rate_or_dict.copy()
-                for qb in qubits:
-                    if qb not in rate_or_dict:
-                        out[qb] = 0.0
-                return out
-            else:
-                return {qb: rate_or_dict for qb in qubits}
-
-        heat_rate_GHz = _as_rate_dict(heat_rate_GHz)
-        cool_rate_GHz = _as_rate_dict(cool_rate_GHz)
-        dephase_rate_GHz = _as_rate_dict(dephase_rate_GHz)
+        heat_rate_GHz = _as_rate_dict(heat_rate_GHz, qubits)
+        cool_rate_GHz = _as_rate_dict(cool_rate_GHz, qubits)
+        dephase_rate_GHz = _as_rate_dict(dephase_rate_GHz, qubits)
 
         for q in qubits:
             gamma_h = heat_rate_GHz[q]
@@ -236,14 +243,7 @@ def noisy_moment(
                 # Only non-virtual gates get noise applied.
                 continue
             rates = self.rate_matrix_GHz[qubit] * moment_ns
-            num_op = np.diag(np.sqrt(np.diag(rates)))
-            annihilation = np.sqrt(np.triu(rates, 1))
-            creation = np.sqrt(np.triu(rates.T, 1)).T
-            # Lindbladian with three Lindblad ops for the three processes
-            # Note: 'time' parameter already specified implicitly through rates
-            L = _lindbladian(annihilation) + _lindbladian(creation) + 2 * _lindbladian(num_op)
-            superop = expm(L.real)
-            kraus_ops = qis.superoperator_to_kraus(superop)
+            kraus_ops = _kraus_ops_from_rates(rates)
             noise_ops.append(ops.KrausChannel(kraus_ops).on(qubit))
         if not noise_ops:
             return [moment]

From 76c2df1d5e07473f7b53d3ee7e6a1095682fd348 Mon Sep 17 00:00:00 2001
From: Orion Martin <40585662+95-martin-orion@users.noreply.github.com>
Date: Mon, 20 Dec 2021 12:21:07 -0800
Subject: [PATCH 15/17] Cache kraus_ops_from_rates

---
 cirq-core/cirq/devices/thermal_noise_model.py | 27 ++++++++++++++-----
 1 file changed, 20 insertions(+), 7 deletions(-)

diff --git a/cirq-core/cirq/devices/thermal_noise_model.py b/cirq-core/cirq/devices/thermal_noise_model.py
index 0c279d28ff3..ef7adb92909 100644
--- a/cirq-core/cirq/devices/thermal_noise_model.py
+++ b/cirq-core/cirq/devices/thermal_noise_model.py
@@ -12,10 +12,11 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-from dataclasses import dataclass
-from typing import TYPE_CHECKING, Dict, List, Optional, Sequence, Set, Union
-from scipy.linalg import expm
+import dataclasses
+import functools
+from typing import TYPE_CHECKING, Dict, List, Optional, Sequence, Set, Tuple, Union
 import numpy as np
+import scipy.linalg
 
 from cirq import devices, ops, protocols, qis
 from cirq.devices.noise_utils import (
@@ -73,14 +74,26 @@ def _lindbladian(left_op: np.ndarray) -> np.ndarray:
     return out
 
 
-def _kraus_ops_from_rates(rates: np.ndarray) -> Sequence[np.ndarray]:
+@functools.lru_cache
+def _kraus_ops_from_rates(
+    flat_rates: Tuple[float, ...], shape: Tuple[int, int]
+) -> Sequence[np.ndarray]:
+    """Generate kraus operators from an array of rates.
+
+    Args:
+        flat_rates: A tuple of rates, flattened from a numpy array with:
+            flat_rates = tuple(rates.reshape(-1))
+            This format is necessary to support caching of inputs.
+        shape: The shape of flat_rates prior to flattening.
+    """
+    rates = np.array(flat_rates).reshape(shape)
     num_op = np.diag(np.sqrt(np.diag(rates)))
     annihilation = np.sqrt(np.triu(rates, 1))
     creation = np.sqrt(np.triu(rates.T, 1)).T
     # Lindbladian with three Lindblad ops for the three processes
     # Note: 'time' parameter already specified implicitly through rates
     L = _lindbladian(annihilation) + _lindbladian(creation) + 2 * _lindbladian(num_op)
-    superop = expm(L.real)
+    superop = scipy.linalg.expm(L.real)
     return qis.superoperator_to_kraus(superop)
 
 
@@ -141,7 +154,7 @@ def _validate_rates(qubits: Set['cirq.Qid'], rates: Dict['cirq.Qid', np.ndarray]
             )
 
 
-@dataclass
+@dataclasses.dataclass
 class ThermalNoiseModel(devices.NoiseModel):
     """NoiseModel representing simulated thermalization of a qubit.
 
@@ -243,7 +256,7 @@ def noisy_moment(
                 # Only non-virtual gates get noise applied.
                 continue
             rates = self.rate_matrix_GHz[qubit] * moment_ns
-            kraus_ops = _kraus_ops_from_rates(rates)
+            kraus_ops = _kraus_ops_from_rates(tuple(rates.reshape(-1)), rates.shape)
             noise_ops.append(ops.KrausChannel(kraus_ops).on(qubit))
         if not noise_ops:
             return [moment]

From bcdbb4a67b9c4e0c13b7f0b83d6172217e16fc66 Mon Sep 17 00:00:00 2001
From: Orion Martin <40585662+95-martin-orion@users.noreply.github.com>
Date: Mon, 20 Dec 2021 13:26:14 -0800
Subject: [PATCH 16/17] Document edge behaviors.

---
 cirq-core/cirq/devices/thermal_noise_model.py | 23 +++++++++++--------
 1 file changed, 13 insertions(+), 10 deletions(-)

diff --git a/cirq-core/cirq/devices/thermal_noise_model.py b/cirq-core/cirq/devices/thermal_noise_model.py
index ef7adb92909..a3190e9b27c 100644
--- a/cirq-core/cirq/devices/thermal_noise_model.py
+++ b/cirq-core/cirq/devices/thermal_noise_model.py
@@ -74,7 +74,7 @@ def _lindbladian(left_op: np.ndarray) -> np.ndarray:
     return out
 
 
-@functools.lru_cache
+@functools.lru_cache(maxsize=256)
 def _kraus_ops_from_rates(
     flat_rates: Tuple[float, ...], shape: Tuple[int, int]
 ) -> Sequence[np.ndarray]:
@@ -129,8 +129,10 @@ def _as_rate_dict(
     rate_or_dict: Optional[Union[float, Dict['cirq.Qid', float]]],
     qubits: Set['cirq.Qid'],
 ) -> Dict['cirq.Qid', float]:
-    # Convert float or None input into dictionary form. Make sure no
-    # qubits are missing from dictionary input.
+    """Convert float or None input into dictionary form.
+
+    This method also ensures that no qubits are missing from dictionary keys.
+    """
     if rate_or_dict is None:
         return {q: 0.0 for q in qubits}
     elif isinstance(rate_or_dict, dict):
@@ -176,8 +178,10 @@ def __init__(
         """Construct a ThermalNoiseModel data object.
 
         Required Args:
-            qubit_dims: Dimension for all qubits in the system.
-                        Currently only supports dimension=2 (qubits, not qudits)
+            qubits: Set of all qubits in the system.
+            gate_durations_ns: Map of gate types to their duration in
+                nanoseconds. These values will override default values for
+                gate duration, if any (e.g. WaitGate).
         Optional Args:
             heat_rate_GHz: single number (units GHz) specifying heating rate,
                 either per qubit, or global value for all.
@@ -225,6 +229,9 @@ def noisy_moment(
         self, moment: 'cirq.Moment', system_qubits: Sequence['cirq.Qid']
     ) -> 'cirq.OP_TREE':
         noise_ops: List['cirq.Operation'] = []
+        # Some devices (including Google hardware) require that all gates have
+        # the same duration, but this does not. Instead, each moment is assumed
+        # to be as long as the longest gate it contains.
         moment_ns: float = 0
         for op in moment:
             op_duration: Optional[float] = None
@@ -232,13 +239,9 @@ def noisy_moment(
                 if not issubclass(type(op.gate), key):
                     continue  # gate type doesn't match
                 # TODO: remove assumption of same time across qubits
-                # if len(key) > 1 and op_data[:1] != key[:1]:
-                #     continue  # qubits don't match
                 op_duration = duration
                 break
-            if op_duration is None:
-                if not isinstance(op.gate, ops.WaitGate):
-                    continue
+            if op_duration is None and isinstance(op.gate, ops.WaitGate):
                 # special case for wait gates if not predefined
                 op_duration = op.gate.duration.total_nanos()
             moment_ns = max(moment_ns, op_duration)

From bd49be5eb0b6e1eebe2d89d780a3c333141b36a9 Mon Sep 17 00:00:00 2001
From: Orion Martin <40585662+95-martin-orion@users.noreply.github.com>
Date: Mon, 20 Dec 2021 13:53:43 -0800
Subject: [PATCH 17/17] Handle other qubits nicely

---
 cirq-core/cirq/devices/thermal_noise_model.py | 22 +++++++++++-----
 .../cirq/devices/thermal_noise_model_test.py  | 25 +++++++++++++++++--
 2 files changed, 39 insertions(+), 8 deletions(-)

diff --git a/cirq-core/cirq/devices/thermal_noise_model.py b/cirq-core/cirq/devices/thermal_noise_model.py
index a3190e9b27c..cb6a0ad9d1c 100644
--- a/cirq-core/cirq/devices/thermal_noise_model.py
+++ b/cirq-core/cirq/devices/thermal_noise_model.py
@@ -228,6 +228,20 @@ def __init__(
     def noisy_moment(
         self, moment: 'cirq.Moment', system_qubits: Sequence['cirq.Qid']
     ) -> 'cirq.OP_TREE':
+        if not moment.operations:
+            return [moment]
+        if self.require_physical_tag:
+            physical_ops = [PHYSICAL_GATE_TAG in op.tags for op in moment]
+            if any(physical_ops):
+                if not all(physical_ops):
+                    raise ValueError(
+                        "Moments are expected to be all physical or all virtual ops, "
+                        f"but found {moment.operations}"
+                    )
+            else:
+                # Only moments with physical operations should have noise.
+                return [moment]
+
         noise_ops: List['cirq.Operation'] = []
         # Some devices (including Google hardware) require that all gates have
         # the same duration, but this does not. Instead, each moment is assumed
@@ -244,20 +258,16 @@ def noisy_moment(
             if op_duration is None and isinstance(op.gate, ops.WaitGate):
                 # special case for wait gates if not predefined
                 op_duration = op.gate.duration.total_nanos()
-            moment_ns = max(moment_ns, op_duration)
+            if op_duration is not None:
+                moment_ns = max(moment_ns, op_duration)
 
         if moment_ns == 0:
             return [moment]
 
         for qubit in system_qubits:
             qubit_op = moment.operation_at(qubit)
-            if qubit_op is None:
-                continue
             if self.skip_measurements and protocols.is_measurement(qubit_op):
                 continue
-            if self.require_physical_tag and PHYSICAL_GATE_TAG not in qubit_op.tags:
-                # Only non-virtual gates get noise applied.
-                continue
             rates = self.rate_matrix_GHz[qubit] * moment_ns
             kraus_ops = _kraus_ops_from_rates(tuple(rates.reshape(-1)), rates.shape)
             noise_ops.append(ops.KrausChannel(kraus_ops).on(qubit))
diff --git a/cirq-core/cirq/devices/thermal_noise_model_test.py b/cirq-core/cirq/devices/thermal_noise_model_test.py
index e1c0705a45e..182fba026fe 100644
--- a/cirq-core/cirq/devices/thermal_noise_model_test.py
+++ b/cirq-core/cirq/devices/thermal_noise_model_test.py
@@ -107,6 +107,25 @@ def test_incomplete_rates():
     assert np.allclose(model.rate_matrix_GHz[q1], np.array([[0, 0], [1e-5, 0]]))
 
 
+def test_noise_from_empty_moment():
+    # Verify that a moment with no duration has no noise.
+    q0, q1 = cirq.LineQubit.range(2)
+    gate_durations = {}
+    heat_rate_GHz = {q1: 1e-5}
+    cool_rate_GHz = {q0: 1e-4}
+    model = ThermalNoiseModel(
+        qubits={q0, q1},
+        gate_durations_ns=gate_durations,
+        heat_rate_GHz=heat_rate_GHz,
+        cool_rate_GHz=cool_rate_GHz,
+        dephase_rate_GHz=None,
+        require_physical_tag=False,
+        skip_measurements=False,
+    )
+    moment = cirq.Moment()
+    assert model.noisy_moment(moment, system_qubits=[q0, q1]) == [moment]
+
+
 def test_noise_from_zero_duration():
     # Verify that a moment with no duration has no noise.
     q0, q1 = cirq.LineQubit.range(2)
@@ -146,7 +165,8 @@ def test_noise_from_virtual_gates():
     assert model.noisy_moment(moment, system_qubits=[q0, q1]) == [moment]
 
     part_virtual_moment = cirq.Moment(cirq.Z(q0), cirq.Z(q1).with_tags(PHYSICAL_GATE_TAG))
-    assert len(model.noisy_moment(part_virtual_moment, system_qubits=[q0, q1])) == 2
+    with pytest.raises(ValueError, match="all physical or all virtual"):
+        _ = model.noisy_moment(part_virtual_moment, system_qubits=[q0, q1])
 
     model.require_physical_tag = False
     assert len(model.noisy_moment(moment, system_qubits=[q0, q1])) == 2
@@ -198,7 +218,8 @@ def test_noisy_moment_one_qubit():
     moment = cirq.Moment(gate.on(q0))
     noisy_moment = model.noisy_moment(moment, system_qubits=[q0, q1])
     assert noisy_moment[0] == moment
-    assert len(noisy_moment[1]) == 1
+    # Noise applies to both qubits, even if only one is acted upon.
+    assert len(noisy_moment[1]) == 2
     noisy_choi = cirq.kraus_to_choi(cirq.kraus(noisy_moment[1].operations[0]))
     assert np.allclose(
         noisy_choi,