diff --git a/lib/BUILD b/lib/BUILD
index 5daea28e..3bc63597 100644
--- a/lib/BUILD
+++ b/lib/BUILD
@@ -8,6 +8,9 @@ cc_library(
     hdrs = [
         "bits.h",
         "bitstring.h",
+        "channel.h",
+        "channels_cirq.h",
+        "circuit_noisy.h",
         "circuit_qsim_parser.h",
         "circuit.h",
         "expect.h",
@@ -395,12 +398,42 @@ cc_library(
     ],
 )
 
+### Channel and noisy circuit libraries ###
+
+cc_library(
+    name = "channel",
+    hdrs = ["channel.h"],
+    deps = [":gate"],
+)
+
+cc_library(
+    name = "circuit_noisy",
+    hdrs = ["circuit_noisy.h"],
+    deps = [
+        ":channel",
+        ":circuit",
+    ],
+)
+
+cc_library(
+    name = "channels_cirq",
+    hdrs = ["channels_cirq.h"],
+    deps = [
+        ":channel",
+        ":gates_cirq",
+    ],
+)
+
 ### Quantum trajectory simulator ###
 
 cc_library(
     name = "qtrajectory",
     hdrs = ["qtrajectory.h"],
-    deps = [":gate"],
+    deps = [
+        ":circuit_noisy",
+        ":gate",
+        ":gate_appl",
+    ],
 )
 
 ### UnitarySpace libraries ###
diff --git a/lib/channel.h b/lib/channel.h
new file mode 100644
index 00000000..5f2a187b
--- /dev/null
+++ b/lib/channel.h
@@ -0,0 +1,81 @@
+// Copyright 2019 Google LLC. All Rights Reserved.
+//
+// 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.
+
+#ifndef CHANNEL_H_
+#define CHANNEL_H_
+
+#include "gate.h"
+
+namespace qsim {
+
+/**
+ * Kraus operator.
+ */
+template <typename Gate>
+struct KrausOperator {
+  enum Kind {
+    kNormal = 0,
+    kMeasurement = gate::kMeasurement,
+  };
+
+  /**
+   * Kraus operator type;
+   */
+  Kind kind;
+
+  /**
+   * If true, the Kraus operator is a unitary operator times a constant.
+   */
+  bool unitary;
+
+  /**
+   * Lower bound on Kraus operator probability.
+   */
+  double prob;
+
+  /**
+   * Sequence of operations that represent the Kraus operator. This can be just
+   * one operation.
+   */
+  std::vector<Gate> ops;
+};
+
+/**
+ * Quantum channel.
+ */
+template <typename Gate>
+using Channel = std::vector<KrausOperator<Gate>>;
+
+/**
+ * Makes a channel from the gate.
+ * @param time The time to place the channel at.
+ * @param gate The input gate.
+ * @return The output channel.
+ */
+template <typename Gate>
+Channel<Gate> MakeChannelFromGate(unsigned time, const Gate& gate) {
+  auto normal = KrausOperator<Gate>::kNormal;
+  auto measurement = KrausOperator<Gate>::kMeasurement;
+
+  auto kind = gate.kind == gate::kMeasurement ? measurement : normal;
+
+  Channel<Gate> channel = {{kind, true, 1, {gate}}};
+  channel[0].ops[0].time = time;
+
+  return channel;
+}
+
+}  // namespace qsim
+
+#endif  // CHANNEL_H_
diff --git a/lib/channels_cirq.h b/lib/channels_cirq.h
new file mode 100644
index 00000000..a8fd87ce
--- /dev/null
+++ b/lib/channels_cirq.h
@@ -0,0 +1,442 @@
+// Copyright 2019 Google LLC. All Rights Reserved.
+//
+// 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.
+
+#ifndef CHANNELS_CIRQ_H_
+#define CHANNELS_CIRQ_H_
+
+#include <cmath>
+#include <cstdint>
+#include <vector>
+
+#include "channel.h"
+#include "gates_cirq.h"
+
+namespace qsim {
+
+namespace Cirq {
+
+template <typename fp_type>
+using Channel = qsim::Channel<GateCirq<fp_type>>;
+
+/**
+ * Asymmetric depolarizing channel factory.
+ */
+template <typename fp_type>
+struct AsymmetricDepolarizingChannel {
+  static constexpr char name[] = "asymmetric_depolarize";
+
+  AsymmetricDepolarizingChannel(double p_x, double p_y, double p_z)
+      : p_x(p_x), p_y(p_y), p_z(p_z) {}
+
+  static Channel<fp_type> Create(unsigned time, unsigned q,
+                                 double p_x, double p_y, double p_z) {
+    double p1 = 1 - p_x - p_y - p_z;
+
+    auto normal = KrausOperator<GateCirq<fp_type>>::kNormal;
+
+    return {{normal, 1, p1, {}},
+            {normal, 1, p_x, {X<fp_type>::Create(time, q)}},
+            {normal, 1, p_y, {Y<fp_type>::Create(time, q)}},
+            {normal, 1, p_z, {Z<fp_type>::Create(time, q)}}};
+  }
+
+  static Channel<fp_type> Create(unsigned time,
+                                 const std::vector<unsigned>& qubits,
+                                 double p_x, double p_y, double p_z) {
+    double p1 = 1 - p_x - p_y - p_z;
+
+    auto normal = KrausOperator<GateCirq<fp_type>>::kNormal;
+
+    uint64_t size = uint64_t{1} << (2 * qubits.size());
+
+    Channel<fp_type> channel;
+    channel.reserve(size);
+
+    for (uint64_t i = 0; i < size; ++i) {
+      channel.push_back({normal, 1, 0, {}});
+      auto& kop = channel.back();
+
+      kop.ops.reserve(qubits.size());
+
+      double prob = 1;
+
+      for (unsigned q = 0; q < qubits.size(); ++q) {
+        unsigned pauli_index = (i >> (2 * q)) & 3;
+
+        switch (pauli_index) {
+        case 0:
+          prob *= p1;
+          break;
+        case 1:
+          prob *= p_x;
+          kop.ops.push_back(X<fp_type>::Create(time, q));
+          break;
+        case 2:
+          prob *= p_y;
+          kop.ops.push_back(Y<fp_type>::Create(time, q));
+          break;
+        case 3:
+          prob *= p_z;
+          kop.ops.push_back(Z<fp_type>::Create(time, q));
+          break;
+        }
+      }
+
+      kop.prob = prob;
+    }
+
+    return channel;
+  }
+
+  Channel<fp_type> Create(unsigned time, unsigned q) const {
+    return Create(time, q, p_x, p_y, p_z);
+  }
+
+  Channel<fp_type> Create(
+      unsigned time, const std::vector<unsigned>& qubits) const {
+    return Create(time, qubits, p_x, p_y, p_z);
+  }
+
+  double p_x = 0;
+  double p_y = 0;
+  double p_z = 0;
+};
+
+/**
+ * Returns an asymmetric depolarizing channel factory object.
+ */
+template <typename fp_type>
+inline AsymmetricDepolarizingChannel<fp_type> asymmetric_depolarize(
+    double p_x, double p_y, double p_z) {
+  return AsymmetricDepolarizingChannel<fp_type>(p_x, p_y, p_z);
+}
+
+/**
+ * Depolarizing channel factory.
+ */
+template <typename fp_type>
+struct DepolarizingChannel {
+  static constexpr char name[] = "depolarize";
+
+  DepolarizingChannel(double p) : p(p) {}
+
+  static Channel<fp_type> Create(unsigned time, unsigned q, double p) {
+    double p1 = 1 - p;
+    double p2 = p / 3;
+
+    auto normal = KrausOperator<GateCirq<fp_type>>::kNormal;
+
+    return {{normal, 1, p1, {}},
+            {normal, 1, p2, {X<fp_type>::Create(time, q)}},
+            {normal, 1, p2, {Y<fp_type>::Create(time, q)}},
+            {normal, 1, p2, {Z<fp_type>::Create(time, q)}}};
+  }
+
+  static Channel<fp_type> Create(
+      unsigned time, const std::vector<unsigned>& qubits, double p) {
+    double p1 = 1 - p;
+    double p2 = p / 3;
+
+    auto normal = KrausOperator<GateCirq<fp_type>>::kNormal;
+
+    uint64_t size = uint64_t{1} << (2 * qubits.size());
+
+    Channel<fp_type> channel;
+    channel.reserve(size);
+
+    for (uint64_t i = 0; i < size; ++i) {
+      channel.push_back({normal, 1, 0, {}});
+      auto& kop = channel.back();
+
+      kop.ops.reserve(qubits.size());
+
+      double prob = 1;
+
+      for (unsigned q = 0; q < qubits.size(); ++q) {
+        unsigned pauli_index = (i >> (2 * q)) & 3;
+
+        switch (pauli_index) {
+        case 0:
+          prob *= p1;
+          break;
+        case 1:
+          prob *= p2;
+          kop.ops.push_back(X<fp_type>::Create(time, q));
+          break;
+        case 2:
+          prob *= p2;
+          kop.ops.push_back(Y<fp_type>::Create(time, q));
+          break;
+        case 3:
+          prob *= p2;
+          kop.ops.push_back(Z<fp_type>::Create(time, q));
+          break;
+        }
+      }
+
+      kop.prob = prob;
+    }
+
+    return channel;
+  }
+
+  Channel<fp_type> Create(unsigned time, unsigned q) const {
+    return Create(time, q, p);
+  }
+
+  Channel<fp_type> Create(
+      unsigned time, const std::vector<unsigned>& qubits) const {
+    return Create(time, qubits, p);
+  }
+
+  double p = 0;
+};
+
+/**
+ * Returns a depolarizing channel factory object.
+ */
+template <typename fp_type>
+inline DepolarizingChannel<fp_type> depolarize(double p) {
+  return DepolarizingChannel<fp_type>(p);
+}
+
+/**
+ * Generalized amplitude damping channel factory.
+ */
+template <typename fp_type>
+struct GeneralizedAmplitudeDampingChannel {
+  static constexpr char name[] = "generalized_amplitude_damp";
+
+  GeneralizedAmplitudeDampingChannel(double p, double gamma)
+      : p(p), gamma(gamma) {}
+
+  static Channel<fp_type> Create(
+      unsigned time, unsigned q, double p, double gamma) {
+    double p1 = p * (1 - gamma);
+    double p2 = (1 - p) * (1 - gamma);
+    double p3 = 0;
+
+    fp_type t1 = std::sqrt(p);
+    fp_type r1 = std::sqrt(p * (1 - gamma));
+    fp_type s1 = std::sqrt(p * gamma);
+    fp_type t2 = std::sqrt(1 - p);
+    fp_type r2 = std::sqrt((1 - p) * (1 - gamma));
+    fp_type s2 = std::sqrt((1 - p) * gamma);
+
+    using M = Cirq::MatrixGate1<fp_type>;
+    auto normal = KrausOperator<GateCirq<fp_type>>::kNormal;
+
+
+    return {{normal, 0, p1, {M::Create(time, q, {t1, 0,  0, 0,  0, 0, r1, 0})}},
+            {normal, 0, p2, {M::Create(time, q, {r2, 0,  0, 0,  0, 0, t2, 0})}},
+            {normal, 0, p3, {M::Create(time, q, { 0, 0, s1, 0,  0,  0, 0, 0})}},
+            {normal, 0, p3, {M::Create(time, q, { 0, 0,  0, 0, s2, 0,  0, 0})}},
+           };
+  }
+
+  Channel<fp_type> Create(unsigned time, unsigned q) const {
+    return Create(time, q, p, gamma);
+  }
+
+  double p = 1;
+  double gamma = 0;
+};
+
+/**
+ * Returns a generalized amplitude damping channel factory object.
+ */
+template <typename fp_type>
+inline GeneralizedAmplitudeDampingChannel<fp_type> generalized_amplitude_damp(
+    double p, double gamma) {
+  return GeneralizedAmplitudeDampingChannel<fp_type>(p, gamma);
+}
+
+/**
+ * Amplitude damping channel factory.
+ */
+template <typename fp_type>
+struct AmplitudeDampingChannel {
+  static constexpr char name[] = "amplitude_damp";
+
+  AmplitudeDampingChannel(double gamma) : gamma(gamma) {}
+
+  static Channel<fp_type> Create(unsigned time, unsigned q, double gamma) {
+    double p1 = 1 - gamma;
+    double p2 = 0;
+
+    fp_type r = std::sqrt(p1);
+    fp_type s = std::sqrt(gamma);
+
+    using M = Cirq::MatrixGate1<fp_type>;
+    auto normal = KrausOperator<GateCirq<fp_type>>::kNormal;
+
+    return {{normal, 0, p1, {M::Create(time, q, {1, 0, 0, 0, 0, 0, r, 0})}},
+            {normal, 0, p2, {M::Create(time, q, {0, 0, s, 0, 0, 0, 0, 0})}},
+           };
+  }
+
+  Channel<fp_type> Create(unsigned time, unsigned q) const {
+    return Create(time, q, gamma);
+  }
+
+  double gamma = 0;
+};
+
+/**
+ * Returns an amplitude damping channel factory object.
+ */
+template <typename fp_type>
+inline AmplitudeDampingChannel<fp_type> amplitude_damp(double gamma) {
+  return AmplitudeDampingChannel<fp_type>(gamma);
+}
+
+/**
+ *  Phase damping channel factory.
+ */
+template <typename fp_type>
+struct PhaseDampingChannel {
+  static constexpr char name[] = "phase_dump";
+
+  PhaseDampingChannel(double gamma) : gamma(gamma) {}
+
+  static Channel<fp_type> Create(unsigned time, unsigned q, double gamma) {
+    double p1 = 1 - gamma;
+    double p2 = 0;
+
+    fp_type r = std::sqrt(p1);
+    fp_type s = std::sqrt(gamma);
+
+    using M = Cirq::MatrixGate1<fp_type>;
+    auto normal = KrausOperator<GateCirq<fp_type>>::kNormal;
+
+    return {{normal, 0, p1, {M::Create(time, q, {1, 0, 0, 0, 0, 0, r, 0})}},
+            {normal, 0, p2, {M::Create(time, q, {0, 0, 0, 0, 0, 0, s, 0})}},
+           };
+  }
+
+  Channel<fp_type> Create(unsigned time, unsigned q) const {
+    return Create(time, q, gamma);
+  }
+
+  double gamma = 0;
+};
+
+/**
+ * Returns a phase damping channel factory object.
+ */
+template <typename fp_type>
+inline PhaseDampingChannel<fp_type> phase_damp(double gamma) {
+  return PhaseDampingChannel<fp_type>(gamma);
+}
+
+/**
+ *  Reset channel factory.
+ */
+template <typename fp_type>
+struct ResetChannel {
+  static constexpr char name[] = "reset";
+
+  static Channel<fp_type> Create(unsigned time, unsigned q) {
+    using M = Cirq::MatrixGate1<fp_type>;
+    auto normal = KrausOperator<GateCirq<fp_type>>::kNormal;
+
+    return {{normal, 0, 0, {M::Create(time, q, {1, 0, 0, 0, 0, 0, 0, 0})}},
+            {normal, 0, 0, {M::Create(time, q, {0, 0, 1, 0, 0, 0, 0, 0})}},
+           };
+  }
+};
+
+/**
+ * Returns a reset channel factory object.
+ */
+template <typename fp_type>
+inline ResetChannel<fp_type> reset() {
+  return ResetChannel<fp_type>();
+}
+
+/**
+ *  Phase flip channel factory.
+ */
+template <typename fp_type>
+struct PhaseFlipChannel {
+  static constexpr char name[] = "phase_flip";
+
+  PhaseFlipChannel(double p) : p(p) {}
+
+  static Channel<fp_type> Create(unsigned time, unsigned q, double p) {
+    double p1 = 1 - p;
+    double p2 = p;
+
+    auto normal = KrausOperator<GateCirq<fp_type>>::kNormal;
+
+    return {{normal, 1, p1, {}},
+            {normal, 1, p2, {Z<fp_type>::Create(time, q)}}
+           };
+  }
+
+  Channel<fp_type> Create(unsigned time, unsigned q) const {
+    return Create(time, q, p);
+  }
+
+  double p = 0;
+};
+
+/**
+ * Returns a phase flip channel factory object.
+ */
+template <typename fp_type>
+inline PhaseFlipChannel<fp_type> phase_flip(double p) {
+  return PhaseFlipChannel<fp_type>(p);
+}
+
+/**
+ *  Bit flip channel factory.
+ */
+template <typename fp_type>
+struct BitFlipChannel {
+  static constexpr char name[] = "bit_flip";
+
+  BitFlipChannel(double p) : p(p) {}
+
+  static Channel<fp_type> Create(unsigned time, unsigned q, double p) {
+    double p1 = 1 - p;
+    double p2 = p;
+
+    auto normal = KrausOperator<GateCirq<fp_type>>::kNormal;
+
+    return {{normal, 1, p1, {}},
+            {normal, 1, p2, {X<fp_type>::Create(time, q)}}
+           };
+  }
+
+  Channel<fp_type> Create(unsigned time, unsigned q) const {
+    return Create(time, q, p);
+  }
+
+  double p = 0;
+};
+
+/**
+ * Returns a bit flip channel factory object.
+ */
+template <typename fp_type>
+inline BitFlipChannel<fp_type> bit_flip(double p) {
+  return BitFlipChannel<fp_type>(p);
+}
+
+}  // namesapce Cirq
+
+}  // namespace qsim
+
+#endif  // CHANNELS_CIRQ_H_
diff --git a/lib/circuit_noisy.h b/lib/circuit_noisy.h
new file mode 100644
index 00000000..646f5be5
--- /dev/null
+++ b/lib/circuit_noisy.h
@@ -0,0 +1,79 @@
+// Copyright 2019 Google LLC. All Rights Reserved.
+//
+// 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.
+
+#ifndef CIRCUIT_NOISY_H_
+#define CIRCUIT_NOISY_H__
+
+#include <vector>
+
+#include "circuit.h"
+#include "channel.h"
+
+namespace qsim {
+
+/**
+ * Noisy circuit.
+ */
+template <typename Gate>
+using NoisyCircuit = std::vector<Channel<Gate>>;
+
+/**
+ * Makes a noisy circuit from the clean circuit.
+ * Channels are added after each qubit of each gate of the clean cicuit.
+ * Roughly equivalent to cirq.Circuit.with_noise.
+ * @param num_qubits The number of circuit qubits.
+ * @param gates The circuit gates.
+ * @param A channel factory to construct channels.
+ * @return The output noisy circuit.
+ */
+template <typename Gate, typename ChannelFactory>
+inline NoisyCircuit<Gate> MakeNoisy(unsigned num_qubits,
+                                    const std::vector<Gate>& gates,
+                                    const ChannelFactory& channel_factory) {
+  NoisyCircuit<Gate> ncircuit;
+  ncircuit.reserve(4 * gates.size());
+
+  for (const auto& gate : gates) {
+    ncircuit.push_back(MakeChannelFromGate(2 * gate.time, gate));
+
+    for (auto q : gate.qubits) {
+      ncircuit.push_back(channel_factory.Create(2 * gate.time + 1, q));
+    }
+
+    for (auto q : gate.controlled_by) {
+      ncircuit.push_back(channel_factory.Create(2 * gate.time + 1, q));
+    }
+  }
+
+  return ncircuit;
+}
+
+/**
+ * Makes a noisy circuit from the clean circuit.
+ * Channels are added after each qubit of each gate of the clean cicuit.
+ * Roughly equivalent to cirq.Circuit.with_noise.
+ * @param num_qubits The number of circuit qubits.
+ * @param circuit The input cicuit.
+ * @param A channel factory to construct channels.
+ * @return The output noisy circuit.
+ */
+template <typename Gate, typename ChannelFactory>
+inline NoisyCircuit<Gate> MakeNoisy(const Circuit<Gate>& circuit,
+                                    const ChannelFactory& channel_factory) {
+  return MakeNoisy(circuit.num_qubits, circuit.gates, channel_factory);
+}
+
+}  // namespace qsim
+
+#endif  // CIRCUIT_NOISY_H_
diff --git a/lib/qtrajectory.h b/lib/qtrajectory.h
index eed5dc8d..2dc09daf 100644
--- a/lib/qtrajectory.h
+++ b/lib/qtrajectory.h
@@ -20,54 +20,12 @@
 #include <random>
 #include <vector>
 
+#include "circuit_noisy.h"
 #include "gate.h"
+#include "gate_appl.h"
 
 namespace qsim {
 
-/**
- * Kraus operator.
- */
-template <typename Gate>
-struct KrausOperator {
-  enum Kind {
-    kNormal = 0,
-    kMeasurement = gate::kMeasurement,
-  };
-
-  /**
-   * Kraus operator type;
-   */
-  Kind kind;
-
-  /**
-   * If true, the Kraus operator is a unitary operator times a constant.
-   */
-  bool unitary;
-
-  /**
-   * Lower bound on Kraus operator probability.
-   */
-  double prob;
-
-  /**
-   * Sequence of operations that represent the Kraus operator. This can be just
-   * one operation.
-   */
-  std::vector<Gate> ops;
-};
-
-/**
- * Quantum channel.
- */
-template <typename Gate>
-using Channel = std::vector<KrausOperator<Gate>>;
-
-/**
- * Noisy circuit.
- */
-template <typename Gate>
-using NoisyCircuit = std::vector<Channel<Gate>>;
-
 /**
  * Quantum trajectory simulator.
  */
diff --git a/tests/BUILD b/tests/BUILD
index 479d298e..88f7c0ab 100644
--- a/tests/BUILD
+++ b/tests/BUILD
@@ -25,6 +25,26 @@ cc_test(
     ],
 )
 
+cc_test(
+    name = "channels_cirq_test",
+    srcs = ["channels_cirq_test.cc"],
+    deps = [
+        "@com_google_googletest//:gtest_main",
+        "//lib:channels_cirq",
+        "//lib:circuit",
+        "//lib:formux",
+        "//lib:fuser_mqubit",
+        "//lib:gates_cirq",
+        "//lib:io",
+        "//lib:qtrajectory",
+        "//lib:simulator",
+    ],
+    copts = select({
+        ":windows": windows_copts,
+        "//conditions:default": [],
+    }),
+)
+
 cc_test(
     name = "circuit_qsim_parser_test",
     srcs = ["circuit_qsim_parser_test.cc"],
diff --git a/tests/channels_cirq_test.cc b/tests/channels_cirq_test.cc
new file mode 100644
index 00000000..3dc01b80
--- /dev/null
+++ b/tests/channels_cirq_test.cc
@@ -0,0 +1,482 @@
+// Copyright 2019 Google LLC. All Rights Reserved.
+//
+// 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.
+
+#include <cmath>
+#include <cstdint>
+#include <vector>
+
+#include "gtest/gtest.h"
+
+#include "../lib/channels_cirq.h"
+#include "../lib/circuit.h"
+#include "../lib/formux.h"
+#include "../lib/fuser_mqubit.h"
+#include "../lib/gates_cirq.h"
+#include "../lib/io.h"
+#include "../lib/qtrajectory.h"
+#include "../lib/simmux.h"
+
+namespace qsim {
+
+namespace {
+
+using StateSpace = Simulator<For>::StateSpace;
+using State = StateSpace::State;
+using fp_type = StateSpace::fp_type;
+using Gate = Cirq::GateCirq<fp_type>;
+using QTSimulator = QuantumTrajectorySimulator<IO, Gate, MultiQubitGateFuser,
+                                               Simulator<For>>;
+
+Circuit<Gate> CleanCircuit() {
+  using Hd = Cirq::H<fp_type>;
+  using IS = Cirq::ISWAP<fp_type>;
+  using Rx = Cirq::rx<fp_type>;
+  using Ry = Cirq::ry<fp_type>;
+
+  return {
+    2,
+    {
+      Hd::Create(0, 0),
+      Hd::Create(0, 1),
+      IS::Create(1, 0, 1),
+      Rx::Create(2, 0, 0.7),
+      Ry::Create(2, 1, 0.1),
+      IS::Create(3, 0, 1),
+      Ry::Create(4, 0, 0.4),
+      Rx::Create(4, 1, 0.7),
+      IS::Create(5, 0, 1),
+      gate::Measurement<Gate>::Create(6, {0, 1}),
+    },
+  };
+}
+
+void RunBatch(const NoisyCircuit<Gate>& ncircuit,
+              const std::vector<double>& expected_results,
+              unsigned num_reps = 25000) {
+  unsigned num_qubits = 2;
+
+  auto measure = [](uint64_t r, const State& state,
+                    const std::vector<uint64_t>& stat,
+                    std::vector<unsigned>& histogram) {
+    ASSERT_EQ(stat.size(), 1);
+    ++histogram[stat[0]];
+  };
+
+  std::vector<unsigned> histogram(1 << num_qubits, 0);
+
+  QTSimulator::Parameter param;
+  param.collect_mea_stat = true;
+
+  EXPECT_TRUE(QTSimulator::Run(param, num_qubits, ncircuit,
+                                      0, num_reps, measure, histogram));
+
+  for (std::size_t i = 0; i < histogram.size(); ++i) {
+    EXPECT_NEAR(double(histogram[i]) / num_reps, expected_results[i], 0.005);
+  }
+}
+
+template <typename ChannelFactory>
+inline NoisyCircuit<Gate> MakeNoisy2(
+    const std::vector<Gate>& gates, const ChannelFactory& channel_factory) {
+  NoisyCircuit<Gate> ncircuit;
+  ncircuit.reserve(2 * gates.size());
+
+  unsigned prev_time = 0;
+
+  for (const auto& gate : gates) {
+    if (gate.time > prev_time) {
+      ncircuit.push_back(channel_factory.Create(2 * prev_time + 1, {0, 1}));
+      prev_time = gate.time;
+    }
+
+    ncircuit.push_back(MakeChannelFromGate(2 * gate.time, gate));
+  }
+
+  return ncircuit;
+}
+
+}  // namespace
+
+TEST(ChannelsCirqTest, AsymmetricDepolarizingChannel) {
+/* The expected results are obtained with the following Cirq code.
+
+import cirq
+
+qs = cirq.LineQubit.range(2)
+
+channel = cirq.asymmetric_depolarize(0.01, 0.02, 0.05)
+
+ncircuit = cirq.Circuit(
+  cirq.H(qs[0]),
+  cirq.H(qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.ISWAP(qs[0], qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.rx(0.7)(qs[0]),
+  cirq.ry(0.1)(qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.ISWAP(qs[0], qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.ry(0.4)(qs[0]),
+  cirq.rx(0.7)(qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.ISWAP(qs[0], qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.measure(*[qs[0], qs[1]], key='m'),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+)
+
+reps = 10000000
+
+sim = cirq.Simulator()
+res = sim.run(ncircuit, repetitions=reps)
+
+for key, val in sorted(res.histogram(key='m').items()):
+  print(f'{key} {float(val) / reps}')
+*/
+
+  std::vector<double> expected_results = {
+    0.3249745, 0.2416734, 0.1637337, 0.2696184,
+  };
+
+  auto channel = Cirq::asymmetric_depolarize<fp_type>(0.01, 0.02, 0.05);
+  auto circuit = CleanCircuit();
+
+  auto ncircuit = MakeNoisy(circuit, channel);
+  RunBatch(ncircuit, expected_results);
+
+  auto ncircuit2 = MakeNoisy2(circuit.gates, channel);
+  RunBatch(ncircuit2, expected_results);
+}
+
+TEST(ChannelsCirqTest, DepolarizingChannel) {
+/* The expected results are obtained with the following Cirq code.
+
+import cirq
+
+qs = cirq.LineQubit.range(2)
+
+circuit = cirq.Circuit(
+  cirq.H(qs[0]),
+  cirq.H(qs[1]),
+  cirq.ISWAP(qs[0], qs[1]),
+  cirq.rx(0.7)(qs[0]),
+  cirq.ry(0.1)(qs[1]),
+  cirq.ISWAP(qs[0], qs[1]),
+  cirq.ry(0.4)(qs[0]),
+  cirq.rx(0.7)(qs[1]),
+  cirq.ISWAP(qs[0], qs[1]),
+  cirq.measure(*[qs[0], qs[1]], key='m'),
+)
+
+ncircuit = circuit.with_noise(cirq.depolarize(0.01))
+
+reps = 10000000
+
+sim = cirq.Simulator()
+res = sim.run(ncircuit, repetitions=reps)
+
+for key, val in sorted(res.histogram(key='m').items()):
+  print(f'{key} {float(val) / reps}')
+
+*/
+
+  std::vector<double> expected_results = {
+    0.3775343, 0.2423451, 0.0972634, 0.2828572,
+  };
+
+  auto channel = Cirq::depolarize<fp_type>(0.02);
+  auto circuit = CleanCircuit();
+
+  auto ncircuit = MakeNoisy(circuit, channel);
+  RunBatch(ncircuit, expected_results);
+
+  auto ncircuit2 = MakeNoisy2(circuit.gates, channel);
+  RunBatch(ncircuit2, expected_results);
+}
+
+TEST(ChannelsCirqTest, GeneralizedAmplitudeDampingChannel) {
+/* The expected results are obtained with the following Cirq code.
+
+import cirq
+
+qs = cirq.LineQubit.range(2)
+
+channel = cirq.generalized_amplitude_damp(0.5, 0.1)
+
+ncircuit = cirq.Circuit(
+  cirq.H(qs[0]),
+  cirq.H(qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.ISWAP(qs[0], qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.rx(0.7)(qs[0]),
+  cirq.ry(0.1)(qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.ISWAP(qs[0], qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.ry(0.4)(qs[0]),
+  cirq.rx(0.7)(qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.ISWAP(qs[0], qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.measure(*[qs[0], qs[1]], key='m'),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+)
+
+reps = 10000000
+
+sim = cirq.Simulator()
+res = sim.run(ncircuit, repetitions=reps)
+
+for key, val in sorted(res.histogram(key='m').items()):
+  print(f'{key} {float(val) / reps}')
+*/
+
+  std::vector<double> expected_results = {
+    0.318501, 0.260538, 0.164616, 0.256345,
+  };
+
+  auto circuit = CleanCircuit();
+  auto ncircuit = MakeNoisy(
+      circuit, Cirq::generalized_amplitude_damp<fp_type>(0.5, 0.1));
+
+  RunBatch(ncircuit, expected_results);
+}
+
+TEST(ChannelsCirqTest, AmplitudeDampingChannel) {
+/* The expected results are obtained with the following Cirq code.
+
+import cirq
+
+qs = cirq.LineQubit.range(2)
+
+channel = cirq.amplitude_damp(0.05)
+
+ncircuit = cirq.Circuit(
+  cirq.H(qs[0]),
+  cirq.H(qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.ISWAP(qs[0], qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.rx(0.7)(qs[0]),
+  cirq.ry(0.1)(qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.ISWAP(qs[0], qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.ry(0.4)(qs[0]),
+  cirq.rx(0.7)(qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.ISWAP(qs[0], qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.measure(*[qs[0], qs[1]], key='m'),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+)
+
+reps = 10000000
+
+sim = cirq.Simulator()
+res = sim.run(ncircuit, repetitions=reps)
+
+for key, val in sorted(res.histogram(key='m').items()):
+  print(f'{key} {float(val) / reps}')
+*/
+
+  std::vector<double> expected_results = {
+    0.500494, 0.235273, 0.090879, 0.173354,
+  };
+
+  auto circuit = CleanCircuit();
+  auto ncircuit = MakeNoisy(circuit, Cirq::amplitude_damp<fp_type>(0.05));
+
+  RunBatch(ncircuit, expected_results);
+}
+
+TEST(ChannelsCirqTest, PhaseDampingChannel) {
+/* The expected results are obtained with the following Cirq code.
+
+import cirq
+
+qs = cirq.LineQubit.range(2)
+
+channel = cirq.phase_damp(0.02)
+
+ncircuit = cirq.Circuit(
+  cirq.H(qs[0]),
+  cirq.H(qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.ISWAP(qs[0], qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.rx(0.7)(qs[0]),
+  cirq.ry(0.1)(qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.ISWAP(qs[0], qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.ry(0.4)(qs[0]),
+  cirq.rx(0.7)(qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.ISWAP(qs[0], qs[1]),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+  cirq.measure(*[qs[0], qs[1]], key='m'),
+  channel.on(qs[0]),
+  channel.on(qs[1]),
+)
+
+reps = 10000000
+
+sim = cirq.Simulator()
+res = sim.run(ncircuit, repetitions=reps)
+
+for key, val in sorted(res.histogram(key='m').items()):
+  print(f'{key} {float(val) / reps}')
+*/
+
+  std::vector<double> expected_results = {
+    0.412300, 0.230500, 0.057219, 0.299982,
+  };
+
+  auto circuit = CleanCircuit();
+  auto ncircuit = MakeNoisy(circuit, Cirq::phase_damp<fp_type>(0.02));
+
+  RunBatch(ncircuit, expected_results);
+}
+
+TEST(ChannelsCirqTest, ResetChannel) {
+  std::vector<double> expected_results = {
+    1, 0, 0, 0,
+  };
+
+  auto circuit = CleanCircuit();
+  auto ncircuit = MakeNoisy(circuit, Cirq::reset<fp_type>());
+
+  RunBatch(ncircuit, expected_results, 1000);
+}
+
+TEST(ChannelsCirqTest, PhaseFlipChannel) {
+/* The expected results are obtained with the following Cirq code.
+
+import cirq
+
+qs = cirq.LineQubit.range(2)
+
+circuit = cirq.Circuit(
+  cirq.H(qs[0]),
+  cirq.H(qs[1]),
+  cirq.ISWAP(qs[0], qs[1]),
+  cirq.rx(0.7)(qs[0]),
+  cirq.ry(0.1)(qs[1]),
+  cirq.ISWAP(qs[0], qs[1]),
+  cirq.ry(0.4)(qs[0]),
+  cirq.rx(0.7)(qs[1]),
+  cirq.ISWAP(qs[0], qs[1]),
+  cirq.measure(*[qs[0], qs[1]], key='m'),
+)
+
+ncircuit = circuit.with_noise(cirq.phase_flip(0.01))
+
+reps = 10000000
+
+sim = cirq.Simulator()
+res = sim.run(ncircuit, repetitions=reps)
+
+for key, val in sorted(res.histogram(key='m').items()):
+  print(f'{key} {float(val) / reps}')
+
+*/
+
+  std::vector<double> expected_results = {
+    0.3790965, 0.2183726, 0.1037091, 0.2988218,
+  };
+
+  auto circuit = CleanCircuit();
+  auto ncircuit = MakeNoisy(circuit, Cirq::phase_flip<fp_type>(0.05));
+
+  RunBatch(ncircuit, expected_results);
+}
+
+TEST(ChannelsCirqTest, BitFlipChannel) {
+/* The expected results are obtained with the following Cirq code.
+
+import cirq
+
+qs = cirq.LineQubit.range(2)
+
+circuit = cirq.Circuit(
+  cirq.H(qs[0]),
+  cirq.H(qs[1]),
+  cirq.ISWAP(qs[0], qs[1]),
+  cirq.rx(0.7)(qs[0]),
+  cirq.ry(0.1)(qs[1]),
+  cirq.ISWAP(qs[0], qs[1]),
+  cirq.ry(0.4)(qs[0]),
+  cirq.rx(0.7)(qs[1]),
+  cirq.ISWAP(qs[0], qs[1]),
+  cirq.measure(*[qs[0], qs[1]], key='m'),
+)
+
+ncircuit = circuit.with_noise(cirq.bit_flip(0.01))
+
+reps = 10000000
+
+sim = cirq.Simulator()
+res = sim.run(ncircuit, repetitions=reps)
+
+for key, val in sorted(res.histogram(key='m').items()):
+  print(f'{key} {float(val) / reps}')
+
+*/
+
+  std::vector<double> expected_results = {
+    0.389352, 0.242790, 0.081009, 0.286850,
+  };
+
+  auto circuit = CleanCircuit();
+  auto ncircuit = MakeNoisy(circuit, Cirq::bit_flip<fp_type>(0.01));
+
+  RunBatch(ncircuit, expected_results);
+}
+
+}  // namespace qsim
+
+int main(int argc, char** argv) {
+  testing::InitGoogleTest(&argc, argv);
+  return RUN_ALL_TESTS();
+}
diff --git a/tests/make.sh b/tests/make.sh
index 7a992f9d..63d4fcfd 100755
--- a/tests/make.sh
+++ b/tests/make.sh
@@ -21,6 +21,7 @@ path_to_include=googletest/googletest/include
 path_to_lib=googletest/lib
 
 g++ -O3 -I$path_to_include -L$path_to_lib -o bitstring_test.x bitstring_test.cc -lgtest -lpthread
+g++ -O3 -I$path_to_include -L$path_to_lib -mavx2 -mfma -fopenmp -o channels_cirq_test.x channels_cirq_test.cc -lgtest -lpthread
 g++ -O3 -I$path_to_include -L$path_to_lib -o circuit_qsim_parser_test.x circuit_qsim_parser_test.cc -lgtest -lpthread
 g++ -O3 -I$path_to_include -L$path_to_lib -mavx2 -mfma -fopenmp -o expect_test.x expect_test.cc -lgtest -lpthread
 g++ -O3 -I$path_to_include -L$path_to_lib -o fuser_basic_test.x fuser_basic_test.cc -lgtest -lpthread
diff --git a/tests/qtrajectory_test.cc b/tests/qtrajectory_test.cc
index 99c8126b..2b095e2d 100644
--- a/tests/qtrajectory_test.cc
+++ b/tests/qtrajectory_test.cc
@@ -30,38 +30,30 @@ namespace qsim {
 
 namespace {
 
-namespace types {
-
-using StateSpace = Simulator<For>::StateSpace;
-using State = StateSpace::State;
-using fp_type = StateSpace::fp_type;
-using Gate = Cirq::GateCirq<fp_type>;
-using QTSimulator = QuantumTrajectorySimulator<IO, Gate, MultiQubitGateFuser,
+using State = Simulator<For>::State;
+using fp_type = Simulator<For>::fp_type;
+using GateCirq = Cirq::GateCirq<fp_type>;
+using QTSimulator = QuantumTrajectorySimulator<IO, GateCirq,
+                                               MultiQubitGateFuser,
                                                Simulator<For>>;
-using NoisyCircuit = NoisyCircuit<Gate>;
-
-}  // namespace types
 
 void AddBitFlipNoise1(
-    unsigned time, unsigned q, double p, types::NoisyCircuit& ncircuit) {
-  using fp_type = types::Gate::fp_type;
-
+    unsigned time, unsigned q, double p, NoisyCircuit<GateCirq>& ncircuit) {
   double p1 = 1 - p;
   double p2 = p;
 
-  auto normal = KrausOperator<types::Gate>::kNormal;
+  auto normal = KrausOperator<GateCirq>::kNormal;
 
   ncircuit.push_back({{normal, 1, p1, {Cirq::I1<fp_type>::Create(time, q)}},
                       {normal, 1, p2, {Cirq::X<fp_type>::Create(time, q)}}});
 }
 
-void AddBitFlipNoise2(unsigned time, double p, types::NoisyCircuit& ncircuit) {
-  using fp_type = types::Gate::fp_type;
-
+void AddBitFlipNoise2(unsigned time, double p,
+                      NoisyCircuit<GateCirq>& ncircuit) {
   double p1 = 1 - p;
   double p2 = p;
 
-  auto normal = KrausOperator<types::Gate>::kNormal;
+  auto normal = KrausOperator<GateCirq>::kNormal;
 
   ncircuit.push_back({
     {normal, 1, p1 * p1, {Cirq::I1<fp_type>::Create(time, 0),
@@ -82,91 +74,8 @@ void AddBitFlipNoise2(unsigned time, double p, types::NoisyCircuit& ncircuit) {
 //                      {normal, 1, p2, {Cirq::X<fp_type>::Create(time, 1)}}});
 }
 
-void AddPhaseDumpNoise1(
-    unsigned time, unsigned q, double g, types::NoisyCircuit& ncircuit) {
-  using fp_type = types::Gate::fp_type;
-
-  double p1 = 1 - g;
-  double p2 = 0;
-
-  fp_type r = std::sqrt(p1);
-  fp_type s = std::sqrt(g);
-
-  auto normal = KrausOperator<types::Gate>::kNormal;
-
-  using M = Cirq::MatrixGate1<fp_type>;
-
-  ncircuit.push_back(
-      {{normal, 0, p1, {M::Create(time, q, {1, 0, 0, 0, 0, 0, r, 0})}},
-       {normal, 0, p2, {M::Create(time, q, {0, 0, 0, 0, 0, 0, s, 0})}}});
-}
-
-void AddPhaseDumpNoise2(
-    unsigned time, double g, types::NoisyCircuit& ncircuit) {
-  using fp_type = types::Gate::fp_type;
-
-  double p1 = 1 - g;
-  double p2 = 0;
-
-  fp_type r = std::sqrt(p1);
-  fp_type s = std::sqrt(g);
-
-  auto normal = KrausOperator<types::Gate>::kNormal;
-
-  using M = Cirq::MatrixGate1<fp_type>;
-
-  ncircuit.push_back(
-      {{normal, 0, p1, {M::Create(time, 0, {1, 0, 0, 0, 0, 0, r, 0})}},
-       {normal, 0, p2, {M::Create(time, 0, {0, 0, 0, 0, 0, 0, s, 0})}}});
-  ncircuit.push_back(
-      {{normal, 0, p1, {M::Create(time, 1, {1, 0, 0, 0, 0, 0, r, 0})}},
-       {normal, 0, p2, {M::Create(time, 1, {0, 0, 0, 0, 0, 0, s, 0})}}});
-}
-
-void AddAmplDumpNoise1(
-    unsigned time, unsigned q, double g, types::NoisyCircuit& ncircuit) {
-  using fp_type = types::Gate::fp_type;
-
-  double p1 = 1 - g;
-  double p2 = 0;
-
-  fp_type r = std::sqrt(p1);
-  fp_type s = std::sqrt(g);
-
-  auto normal = KrausOperator<types::Gate>::kNormal;
-
-  using M = Cirq::MatrixGate1<fp_type>;
-
-  ncircuit.push_back(
-      {{normal, 0, p1, {M::Create(time, q, {1, 0, 0, 0, 0, 0, r, 0})}},
-       {normal, 0, p2, {M::Create(time, q, {0, 0, s, 0, 0, 0, 0, 0})}}});
-}
-
-void AddAmplDumpNoise2(unsigned time, double g, types::NoisyCircuit& ncircuit) {
-  using fp_type = types::Gate::fp_type;
-
-  double p1 = 1 - g;
-  double p2 = 0;
-
-  fp_type r = std::sqrt(p1);
-  fp_type s = std::sqrt(g);
-
-  auto normal = KrausOperator<types::Gate>::kNormal;
-
-  using M = Cirq::MatrixGate1<fp_type>;
-
-  ncircuit.push_back(
-      {{normal, 0, p1, {M::Create(time, 0, {1, 0, 0, 0, 0, 0, r, 0})}},
-       {normal, 0, p2, {M::Create(time, 0, {0, 0, s, 0, 0, 0, 0, 0})}}});
-  ncircuit.push_back(
-      {{normal, 0, p1, {M::Create(time, 1, {1, 0, 0, 0, 0, 0, r, 0})}},
-       {normal, 0, p2, {M::Create(time, 1, {0, 0, s, 0, 0, 0, 0, 0})}}});
-}
-
 void AddGenAmplDumpNoise1(
-    unsigned time, unsigned q, double g, types::NoisyCircuit& ncircuit) {
-  using fp_type = types::Gate::fp_type;
-
+    unsigned time, unsigned q, double g, NoisyCircuit<GateCirq>& ncircuit) {
   // Probability of exchanging energy with the environment.
   double p = 0.5;
 
@@ -181,7 +90,7 @@ void AddGenAmplDumpNoise1(
   fp_type r2 = std::sqrt((1 - p) * (1 - g));
   fp_type s2 = std::sqrt((1 - p) * g);
 
-  auto normal = KrausOperator<types::Gate>::kNormal;
+  auto normal = KrausOperator<GateCirq>::kNormal;
 
   using M = Cirq::MatrixGate1<fp_type>;
 
@@ -193,9 +102,7 @@ void AddGenAmplDumpNoise1(
 }
 
 void AddGenAmplDumpNoise2(
-    unsigned time, double g, types::NoisyCircuit& ncircuit) {
-  using fp_type = types::Gate::fp_type;
-
+    unsigned time, double g, NoisyCircuit<GateCirq>& ncircuit) {
   // Probability of exchanging energy with the environment.
   double p = 0.5;
 
@@ -210,7 +117,7 @@ void AddGenAmplDumpNoise2(
   fp_type r2 = std::sqrt((1 - p) * (1 - g));
   fp_type s2 = std::sqrt((1 - p) * g);
 
-  auto normal = KrausOperator<types::Gate>::kNormal;
+  auto normal = KrausOperator<GateCirq>::kNormal;
 
   using M = Cirq::MatrixGate1<fp_type>;
 
@@ -227,11 +134,9 @@ void AddGenAmplDumpNoise2(
 }
 
 template <typename AddNoise1, typename AddNoise2>
-types::NoisyCircuit GenerateNoisyCircuit(
+NoisyCircuit<GateCirq> GenerateNoisyCircuit(
     double p, AddNoise1&& add_noise1, AddNoise2&& add_noise2) {
-  using fp_type = types::Gate::fp_type;
-
-  types::NoisyCircuit ncircuit;
+  NoisyCircuit<GateCirq> ncircuit;
   ncircuit.reserve(24);
 
   using Hd = Cirq::H<fp_type>;
@@ -239,7 +144,7 @@ types::NoisyCircuit GenerateNoisyCircuit(
   using Rx = Cirq::rx<fp_type>;
   using Ry = Cirq::ry<fp_type>;
 
-  auto normal = KrausOperator<types::Gate>::kNormal;
+  auto normal = KrausOperator<GateCirq>::kNormal;
 
   ncircuit.push_back({{normal, 1, 1.0, {Hd::Create(0, 0)}}});
   add_noise1(1, 0, p, ncircuit);
@@ -259,19 +164,19 @@ types::NoisyCircuit GenerateNoisyCircuit(
   add_noise1(9, 1, p, ncircuit);
   ncircuit.push_back({{normal, 1, 1.0, {IS::Create(10, 0, 1)}}});
   add_noise2(11, p, ncircuit);
-  ncircuit.push_back({{KrausOperator<types::Gate>::kMeasurement, 1, 1.0,
-                       {gate::Measurement<types::Gate>::Create(12, {0, 1})}}});
+  ncircuit.push_back({{KrausOperator<GateCirq>::kMeasurement, 1, 1.0,
+                       {gate::Measurement<GateCirq>::Create(12, {0, 1})}}});
   add_noise2(13, p, ncircuit);
 
   return ncircuit;
 }
 
-void RunBatch(const types::NoisyCircuit& ncircuit,
-          const std::vector<double>& expected_results) {
+void RunBatch(const NoisyCircuit<GateCirq>& ncircuit,
+              const std::vector<double>& expected_results) {
   unsigned num_qubits = 2;
   unsigned num_reps = 25000;
 
-  auto measure = [](uint64_t r, const types::State& state,
+  auto measure = [](uint64_t r, const State& state,
                     const std::vector<uint64_t>& stat,
                     std::vector<unsigned>& histogram) {
     ASSERT_EQ(stat.size(), 1);
@@ -280,26 +185,26 @@ void RunBatch(const types::NoisyCircuit& ncircuit,
 
   std::vector<unsigned> histogram(1 << num_qubits, 0);
 
-  types::QTSimulator::Parameter param;
+  QTSimulator::Parameter param;
   param.collect_mea_stat = true;
 
-  EXPECT_TRUE(types::QTSimulator::Run(param, num_qubits, ncircuit,
-                                      0, num_reps, measure, histogram));
+  EXPECT_TRUE(QTSimulator::Run(
+      param, num_qubits, ncircuit, 0, num_reps, measure, histogram));
 
   for (std::size_t i = 0; i < histogram.size(); ++i) {
     EXPECT_NEAR(double(histogram[i]) / num_reps, expected_results[i], 0.005);
   }
 }
 
-void RunOnceRepeatedly(const types::NoisyCircuit& ncircuit,
-          const std::vector<double>& expected_results) {
+void RunOnceRepeatedly(const NoisyCircuit<GateCirq>& ncircuit,
+                       const std::vector<double>& expected_results) {
   unsigned num_qubits = 2;
   unsigned num_reps = 25000;
 
-  types::StateSpace state_space(1);
+  Simulator<For>::StateSpace state_space(1);
 
-  types::State scratch = state_space.Null();
-  types::State state = state_space.Create(num_qubits);
+  State scratch = state_space.Null();
+  State state = state_space.Create(num_qubits);
   EXPECT_FALSE(state_space.IsNull(state));
 
   auto state_pointer = state.get();
@@ -308,14 +213,14 @@ void RunOnceRepeatedly(const types::NoisyCircuit& ncircuit,
 
   std::vector<unsigned> histogram(1 << num_qubits, 0);
 
-  types::QTSimulator::Parameter param;
+  QTSimulator::Parameter param;
   param.collect_mea_stat = true;
 
   for (unsigned i = 0; i < num_reps; ++i) {
     state_space.SetStateZero(state);
 
-    EXPECT_TRUE(types::QTSimulator::Run(param, num_qubits, ncircuit, i,
-                                        scratch, state, stat));
+    EXPECT_TRUE(QTSimulator::Run(
+        param, num_qubits, ncircuit, i, scratch, state, stat));
 
     EXPECT_EQ(state_pointer, state.get());
 
@@ -371,116 +276,6 @@ for key, val in sorted(res.histogram(key='m').items()):
   RunBatch(ncircuit1, expected_results);
 }
 
-TEST(QTrajectoryTest, PhaseDump) {
-/* The expected results are obtained with the following Cirq code.
-
-import cirq
-
-qs = cirq.LineQubit.range(2)
-
-channel = cirq.phase_damp(0.02)
-
-ncircuit = cirq.Circuit(
-  cirq.H(qs[0]),
-  cirq.H(qs[1]),
-  channel.on(qs[0]),
-  channel.on(qs[1]),
-  cirq.ISWAP(qs[0], qs[1]),
-  channel.on(qs[0]),
-  channel.on(qs[1]),
-  cirq.rx(0.7)(qs[0]),
-  cirq.ry(0.1)(qs[1]),
-  channel.on(qs[0]),
-  channel.on(qs[1]),
-  cirq.ISWAP(qs[0], qs[1]),
-  channel.on(qs[0]),
-  channel.on(qs[1]),
-  cirq.ry(0.4)(qs[0]),
-  cirq.rx(0.7)(qs[1]),
-  channel.on(qs[0]),
-  channel.on(qs[1]),
-  cirq.ISWAP(qs[0], qs[1]),
-  channel.on(qs[0]),
-  channel.on(qs[1]),
-  cirq.measure(*[qs[0], qs[1]], key='m'),
-  channel.on(qs[0]),
-  channel.on(qs[1]),
-)
-
-reps = 10000000
-
-sim = cirq.Simulator()
-res = sim.run(ncircuit, repetitions=reps)
-
-for key, val in sorted(res.histogram(key='m').items()):
-  print(f'{key} {float(val) / reps}')
-
-*/
-
-  std::vector<double> expected_results = {
-    0.412300, 0.230500, 0.057219, 0.299982,
-  };
-
-  auto ncircuit = GenerateNoisyCircuit(0.02, AddPhaseDumpNoise1,
-                                       AddPhaseDumpNoise2);
-  RunOnceRepeatedly(ncircuit, expected_results);
-}
-
-TEST(QTrajectoryTest, AmplDump) {
-/* The expected results are obtained with the following Cirq code.
-
-import cirq
-
-qs = cirq.LineQubit.range(2)
-
-channel = cirq.amplitude_damp(0.05)
-
-ncircuit = cirq.Circuit(
-  cirq.H(qs[0]),
-  cirq.H(qs[1]),
-  channel.on(qs[0]),
-  channel.on(qs[1]),
-  cirq.ISWAP(qs[0], qs[1]),
-  channel.on(qs[0]),
-  channel.on(qs[1]),
-  cirq.rx(0.7)(qs[0]),
-  cirq.ry(0.1)(qs[1]),
-  channel.on(qs[0]),
-  channel.on(qs[1]),
-  cirq.ISWAP(qs[0], qs[1]),
-  channel.on(qs[0]),
-  channel.on(qs[1]),
-  cirq.ry(0.4)(qs[0]),
-  cirq.rx(0.7)(qs[1]),
-  channel.on(qs[0]),
-  channel.on(qs[1]),
-  cirq.ISWAP(qs[0], qs[1]),
-  channel.on(qs[0]),
-  channel.on(qs[1]),
-  cirq.measure(*[qs[0], qs[1]], key='m'),
-  channel.on(qs[0]),
-  channel.on(qs[1]),
-)
-
-reps = 10000000
-
-sim = cirq.Simulator()
-res = sim.run(ncircuit, repetitions=reps)
-
-for key, val in sorted(res.histogram(key='m').items()):
-  print(f'{key} {float(val) / reps}')
-
-*/
-
-  std::vector<double> expected_results = {
-    0.500494, 0.235273, 0.090879, 0.173354,
-  };
-
-  auto ncircuit = GenerateNoisyCircuit(0.05, AddAmplDumpNoise1,
-                                       AddAmplDumpNoise2);
-  RunBatch(ncircuit, expected_results);
-}
-
 TEST(QTrajectoryTest, GenDump) {
 /* The expected results are obtained with the following Cirq code.
 
@@ -544,14 +339,13 @@ TEST(QTrajectoryTest, CollectKopStat) {
   double p1 = 1 - p;
   double p2 = p;
 
-  using fp_type = types::Gate::fp_type;
   using Hd = Cirq::H<fp_type>;
   using I = Cirq::I1<fp_type>;
   using X = Cirq::X<fp_type>;
 
-  auto normal = KrausOperator<types::Gate>::kNormal;
+  auto normal = KrausOperator<GateCirq>::kNormal;
 
-  types::NoisyCircuit ncircuit;
+  NoisyCircuit<GateCirq> ncircuit;
   ncircuit.reserve(8);
 
   ncircuit.push_back({{normal, 1, 1.0, {Hd::Create(0, 0)}}});
@@ -569,7 +363,7 @@ TEST(QTrajectoryTest, CollectKopStat) {
   ncircuit.push_back({{normal, 1, p1, {I::Create(1, 3)}},
                       {normal, 1, p2, {X::Create(1, 3)}}});
 
-  auto measure = [](uint64_t r, const types::State& state,
+  auto measure = [](uint64_t r, const State& state,
                     const std::vector<uint64_t>& stat,
                     std::vector<std::vector<unsigned>>& histogram) {
     ASSERT_EQ(stat.size(), histogram.size());
@@ -580,11 +374,11 @@ TEST(QTrajectoryTest, CollectKopStat) {
 
   std::vector<std::vector<unsigned>> histogram(8, std::vector<unsigned>(2, 0));
 
-  types::QTSimulator::Parameter param;
+  QTSimulator::Parameter param;
   param.collect_kop_stat = true;
 
-  EXPECT_TRUE(types::QTSimulator::Run(param, num_qubits, ncircuit,
-                                      0, num_reps, measure, histogram));
+  EXPECT_TRUE(QTSimulator::Run(
+      param, num_qubits, ncircuit, 0, num_reps, measure, histogram));
 
   for (std::size_t i = 0; i < 4; ++i) {
     EXPECT_EQ(histogram[i][0], num_reps);
@@ -601,11 +395,9 @@ TEST(QTrajectoryTest, CleanCircuit) {
   unsigned num_qubits = 4;
   auto size = uint64_t{1} << num_qubits;
 
-  std::vector<types::Gate> circuit;
+  std::vector<GateCirq> circuit;
   circuit.reserve(16);
 
-  using fp_type = types::Gate::fp_type;
-
   circuit.push_back(Cirq::H<fp_type>::Create(0, 0));
   circuit.push_back(Cirq::H<fp_type>::Create(0, 1));
   circuit.push_back(Cirq::H<fp_type>::Create(0, 2));
@@ -632,19 +424,19 @@ TEST(QTrajectoryTest, CleanCircuit) {
   circuit.push_back(Cirq::YPowGate<fp_type>::Create(5, 2, 0.9, 0.4));
   circuit.push_back(Cirq::ZPowGate<fp_type>::Create(5, 3, 1.0, 0.5));
 
-  types::NoisyCircuit ncircuit;
+  NoisyCircuit<GateCirq> ncircuit;
   ncircuit.reserve(16);
 
-  auto normal = KrausOperator<types::Gate>::kNormal;
+  auto normal = KrausOperator<GateCirq>::kNormal;
 
   for (std::size_t i = 0; i < circuit.size(); ++i) {
     ncircuit.push_back({{normal, 1, 1.0, {circuit[i]}}});
   }
 
   Simulator<For> simulator(1);
-  types::StateSpace state_space(1);
+  Simulator<For>::StateSpace state_space(1);
 
-  types::State state = state_space.Create(num_qubits);
+  State state = state_space.Create(num_qubits);
   EXPECT_FALSE(state_space.IsNull(state));
 
   state_space.SetStateZero(state);
@@ -654,19 +446,19 @@ TEST(QTrajectoryTest, CleanCircuit) {
     ApplyGate(simulator, gate, state);
   }
 
-  types::State scratch = state_space.Null();
-  types::State nstate = state_space.Create(num_qubits);
+  State scratch = state_space.Null();
+  State nstate = state_space.Create(num_qubits);
   EXPECT_FALSE(state_space.IsNull(nstate));
 
   std::vector<uint64_t> stat;
 
-  types::QTSimulator::Parameter param;
+  QTSimulator::Parameter param;
 
   state_space.SetStateZero(nstate);
 
   // Run quantum trajectory simulator.
-  EXPECT_TRUE(types::QTSimulator::Run(param, num_qubits, ncircuit, 0,
-                                      scratch, nstate, stat));
+  EXPECT_TRUE(QTSimulator::Run(
+      param, num_qubits, ncircuit, 0, scratch, nstate, stat));
 
   EXPECT_EQ(stat.size(), 0);
 
@@ -682,31 +474,30 @@ TEST(QTrajectoryTest, CleanCircuit) {
 TEST(QTrajectoryTest, InitialState) {
   unsigned num_qubits = 3;
 
-  types::NoisyCircuit ncircuit;
+  NoisyCircuit<GateCirq> ncircuit;
   ncircuit.reserve(2);
 
-  using fp_type = types::Gate::fp_type;
-  auto normal = KrausOperator<types::Gate>::kNormal;
+  auto normal = KrausOperator<GateCirq>::kNormal;
 
   ncircuit.push_back({{normal, 1, 1.0, {Cirq::X<fp_type>::Create(0, 0)}}});
   ncircuit.push_back({{normal, 1, 1.0, {Cirq::X<fp_type>::Create(0, 1)}}});
   ncircuit.push_back({{normal, 1, 1.0, {Cirq::X<fp_type>::Create(0, 2)}}});
 
-  types::StateSpace state_space(1);
+  Simulator<For>::StateSpace state_space(1);
 
-  types::State scratch = state_space.Null();
-  types::State state = state_space.Create(num_qubits);
+  State scratch = state_space.Null();
+  State state = state_space.Create(num_qubits);
   EXPECT_FALSE(state_space.IsNull(state));
 
-  types::QTSimulator::Parameter param;
+  QTSimulator::Parameter param;
   std::vector<uint64_t> stat;
 
   for (unsigned i = 0; i < 8; ++i) {
     state_space.SetAmpl(state, i, 1 + i, 0);
   }
 
-  EXPECT_TRUE(types::QTSimulator::Run(param, num_qubits, ncircuit, 0,
-                                      scratch, state, stat));
+  EXPECT_TRUE(QTSimulator::Run(
+      param, num_qubits, ncircuit, 0, scratch, state, stat));
 
   // Expect reversed order of amplitudes.
   for (unsigned i = 0; i < 8; ++i) {