Updated noise model class for local and nonlocal noise, added reset error error type

legacy-build/aer_simulator.cpp

@@ -22,14 +22,11 @@
 #include "base/engine.hpp"
 #include "simulators/qubitvector/qubitvector.hpp"
 #include "simulators/qubitvector/qv_state.hpp"
+#include "framework/interface.hpp"
 
 // Noise
-#include "base/noise.hpp"
+#include "base/noise_model.hpp"
-#include "noise/simple_model.hpp"
-#include "noise/unitary_error.hpp"
-#include "noise/gate_error.hpp"
 
-#include "framework/interface.hpp"
 /*******************************************************************************
  *
  * Main
@@ -75,33 +72,21 @@ int main(int argc, char **argv) {
     using namespace AER;
     using State = QubitVector::State;       // State class
     using Engine = Base::Engine<QV::QubitVector>; // Optimized Engine class
-    using NoiseModel = Noise::SimpleModel;
+    using Base::NoiseModel;
-    
+
     // Initialize simulator
     Base::Controller<Engine, State> sim;
   
     // Check for noise_params
     if (JSON::check_key("config", qobj) &&
         JSON::check_key("noise_params", qobj["config"])) {
-      NoiseModel noise(qobj["config"]["noise_params"]);
+      NoiseModel noise1(qobj["config"]["noise_params"]);
-      out << sim.execute(qobj, &noise).dump(4) << std::endl;
+      out << sim.execute(qobj, &noise1).dump(4) << std::endl;
     } else {
       // execute without noise
       out << sim.execute(qobj).dump(4) << std::endl;
     }
 
-    // Amplitude damping channel
-    /*
-    NoiseModel kraus_noise;
-    std::vector<cmatrix_t> amp_damp(2);
-    double gamma = 0.4;
-    amp_damp[0] = Utils::make_matrix<complex_t>({{{1, 0}, {0, 0}},
-                                                 {{0, 0}, {std::sqrt(gamma), 0}}});
-    amp_damp[1] = Utils::make_matrix<complex_t>({{{0, 0}, {std::sqrt(1-gamma), 0}},
-                                                 {{0, 0}, {0, 0}}});
-    kraus_noise.add_error(Noise::GateError(amp_damp), {"x", "y", "z", "s", "sdg", "h", "t", "tdg", "u1", "u2", "u3"});
-    out << sim.execute(qobj, &kraus_noise).dump(4) << std::endl;
-    */
     return 0;
   } catch (std::exception &e) {
     std::stringstream msg;

src/base/controller.hpp

@@ -29,7 +29,7 @@
 
 // Base Controller
 #include "framework/qobj.hpp"
-#include "base/noise.hpp"
+#include "base/noise_model.hpp"
 
 namespace AER {
 namespace Base {
@@ -72,12 +72,12 @@ public:
   // Load and execute a qobj
   inline json_t execute(const json_t &qobj) {return execute(qobj, nullptr);}
 
-  json_t execute(const json_t &qobj, Noise::Model *noise_ptr);
+  json_t execute(const json_t &qobj, Base::NoiseModel *noise_ptr);
 
   // Execute a single circuit
   json_t execute_circuit(Circuit &circ, 
                          int max_shot_threads,
-                         Noise::Model *noise_ptr);
+                         Base::NoiseModel *noise_ptr);
 
   // Check if measurement sampling can be performed for a circuit
   // and set circuit flag if it is compatible                               
@@ -134,7 +134,7 @@ void Controller<Engine_t, State_t>::load_config(const json_t &config) {
 
 template < class Engine_t, class State_t>
 json_t Controller<Engine_t, State_t>::execute(const json_t &qobj_js,
-                                              Noise::Model *noise_ptr) {
+                                              Base::NoiseModel *noise_ptr) {
   
   auto timer_start = myclock_t::now(); // start timer
   Qobj qobj; // Load QOBJ from json
@@ -209,7 +209,7 @@ json_t Controller<Engine_t, State_t>::execute(const json_t &qobj_js,
 template <class Engine_t, class State_t>
 json_t Controller<Engine_t, State_t>::execute_circuit(Circuit &circ,
                                                       int max_shot_threads,
-                                                      Noise::Model *noise_ptr) {
+                                                      Base::NoiseModel *noise_ptr) {
   
   // Initialize Return
   auto timer_start = myclock_t::now(); // state circuit timer

src/base/engine.hpp

@@ -20,7 +20,7 @@
 #include "framework/snapshot.hpp"
 #include "framework/utils.hpp"
 #include "base/state.hpp"
-#include "base/noise.hpp"
+#include "base/noise_model.hpp"
 
 namespace AER {
 namespace Base {
@@ -56,7 +56,7 @@ public:
   void execute(const Circuit &circ,
                uint_t shots,
                State<state_t> *state_ptr,
-               Noise::Model *noise_ptr = nullptr);
+               Base::NoiseModel *noise_ptr = nullptr);
 
   // This method performs the same function as 'execute', except
   // that it only simulates a single shot and then generates samples
@@ -104,7 +104,7 @@ protected:
   // Apply an operation to the state
   void apply_op(const Operations::Op &op,
                 State<state_t> *state_ptr,
-                Noise::Model *noise_ptr);
+                Base::NoiseModel *noise_ptr);
 
   // Initialize an engine and circuit
   void initialize(State<state_t> *state, const Circuit &circ);
@@ -212,7 +212,7 @@ template <class state_t>
 void Engine<state_t>::execute(const Circuit &circ,
                               uint_t shots,
                               State<state_t> *state_ptr,
-                              Noise::Model *noise_ptr) {
+                              Base::NoiseModel *noise_ptr) {
   // Check if ideal simulation check if sampling is possible
   if (noise_ptr == nullptr && circ.measure_sampling_flag) {
     execute_with_measure_sampling(circ, shots, state_ptr);
@@ -233,7 +233,7 @@ void Engine<state_t>::execute(const Circuit &circ,
 template <class state_t>
 void Engine<state_t>::apply_op(const Operations::Op &op,
                                State<state_t> *state_ptr,
-                               Noise::Model *noise_ptr) {
+                               Base::NoiseModel *noise_ptr) {
   auto it = engine_ops_.find(op.name);
   if (it == engine_ops_.end()) {
     // check if op passes conditional

src/base/noise.hpp

@@ -1,86 +0,0 @@
-/**
- * Copyright 2018, IBM.
- *
- * This source code is licensed under the Apache License, Version 2.0 found in
- * the LICENSE.txt file in the root directory of this source tree.
- */
-
-/**
- * @file    noise.hpp
- * @brief   Noise model base class for qiskit-aer simulator engines
- * @author  Christopher J. Wood <cjwood@us.ibm.com>
- */
-
-#ifndef _aer_base_noise_hpp_
-#define _aer_base_noise_hpp_
-
-#include "framework/operations.hpp"
-#include "framework/types.hpp"
-#include "framework/rng.hpp"
-
-// Work in progress
-
-/* Schema idea?
-
-{
-  "type": "gate_error",   // string
-  "operations": ["x", "y", "z"], // list string
-  "qubits": [[0], [1], ...], // if null or missing gate applies to all qubits
-  "targets": [[qs,..], ...], // if null or missing gate targets are gate qubits
-  "params": [mat1, mat2, mat3,...] // kraus CPTP ops
-}
-
-*/
-
-namespace AER {
-namespace Noise {
-
-using NoiseOps = std::vector<Operations::Op>;
-using NoisePair = std::pair<double, NoiseOps>;
-
-//=========================================================================
-// Error abstract base class
-//=========================================================================
-
-class Error {
-public:
-  // Sample an realization of the error from the error model using the passed
-  // in RNG engine.
-  virtual NoiseOps sample_noise(const Operations::Op &op,
-                                const reg_t &qubits,
-                                RngEngine &rng) = 0;
-
-  // Return the number of error terms
-  inline size_t num_qubits() const {return num_qubits_;}
-
-  // Return the number of error terms
-  inline void set_num_qubits(size_t nq) {num_qubits_ = nq;}
-
-protected:
-  size_t num_qubits_;
-};
-
-//=========================================================================
-// Noise Model abstract base class
-//=========================================================================
-
-class Model {
-public:
-  // Sample noise for the current operation
-  // Base class returns the input operation with no noise (ideal)
-  // derived classes should implement this method
-  virtual NoiseOps sample_noise(const Operations::Op &op) {return {op};}
-
-  // Set the RngEngine seed to a fixed value
-  inline void set_rng_seed(uint_t seed) { rng_ = RngEngine(seed);}
-
-protected:
-  RngEngine rng_; // initialized with random seed
-};
-
-//-------------------------------------------------------------------------
-} // end namespace Noise
-//-------------------------------------------------------------------------
-} // end namespace AER
-//-------------------------------------------------------------------------
-#endif

src/base/noise_model.hpp

@@ -0,0 +1,297 @@
+/**
+ * Copyright 2018, IBM.
+ *
+ * This source code is licensed under the Apache License, Version 2.0 found in
+ * the LICENSE.txt file in the root directory of this source tree.
+ */
+
+/**
+ * @file    noise_model.hpp
+ * @brief   Noise model base class for qiskit-aer simulator engines
+ * @author  Christopher J. Wood <cjwood@us.ibm.com>
+ */
+
+#ifndef _aer_base_noise_model_hpp_
+#define _aer_base_noise_model_hpp_
+
+#include "framework/operations.hpp"
+#include "framework/types.hpp"
+#include "framework/rng.hpp"
+#include "framework/circuit.hpp"
+#include "noise/abstract_error.hpp"
+
+// move JSON parsing
+#include "noise/gate_error.hpp"
+
+
+namespace AER {
+namespace Base {
+
+//=========================================================================
+// Noise Model class
+//=========================================================================
+
+// This allows specification of default and individual noise parameters for
+// each circuit operation that may either be local (applied to the qubits
+// in the operation) or nonlocal (applied to different qubits). The errors
+// may each also be applied either before or after the operation as per
+// the specification of the AbstractError subclass.
+
+class NoiseModel {
+public:
+
+  using NoiseOps = std::vector<Operations::Op>;
+
+  NoiseModel() = default;
+  NoiseModel(const json_t &js) {load_from_json(js);}
+
+  // Sample noise for the current operation
+  virtual NoiseOps sample_noise(const Operations::Op &op);
+
+  // Sample a noisy implementation of a full circuit
+  virtual Circuit sample_noise(const Circuit &circ);
+
+  // Load a noise model from JSON
+  void load_from_json(const json_t &js);
+
+  // Set the RngEngine seed to a fixed value
+  inline void set_rng_seed(uint_t seed) { rng_ = RngEngine(seed);}
+
+  // Add a local error to the noise model for specific qubits
+  template <class DerivedError>
+  void add_local_error(const DerivedError &error, 
+                 const std::vector<std::string> &op_labels,
+                 const std::vector<reg_t> &op_qubit_sets);
+
+
+  // Add a local Error type to the model default for all qubits
+  // without a specific error present
+  template <class DerivedError>
+  void add_local_error(const DerivedError &error, 
+                 const std::vector<std::string> &op_labels) {
+    add_local_error(error, op_labels, {reg_t()});
+  }
+
+  // Add a non-local Error type to the model for specific qubits
+  template <class DerivedError>
+  void add_nonlocal_error(const DerivedError &error, 
+                          const std::vector<std::string> &op_labels,
+                          const std::vector<reg_t> &op_qubit_sets,
+                          const std::vector<reg_t> &noise_qubit_sets);
+
+
+  // Set which single qubit gates should use the X90 waltz error model
+  inline void set_waltz_gates(const std::set<std::string> &waltz_gates) {
+    waltz_gates_ = waltz_gates;
+  }
+
+  // Set threshold for applying u1 rotation angles.
+  // an Op for u1(theta) will only be added if |theta| > 0 and |theta - 2*pi| > 0
+  inline void set_u1_threshold(double threshold) {
+    u1_threshold_ = threshold;
+  }
+
+private:
+  // Flags which say whether the local or nonlocal error tables are used
+  bool local_errors_ = false;
+  bool nonlocal_errors_ = false;
+
+  // Table of errors
+  std::vector<std::unique_ptr<Noise::AbstractError>> error_ptrs_;
+
+  // Table indexes a name with a vector of the position of noise operations
+  // Cant make inner map and unordered map due to lack of hashing for vector<uint_t>
+  using qubit_map_t = std::map<reg_t, std::vector<size_t>>;
+  std::unordered_map<std::string, qubit_map_t> local_error_table_;
+
+  // Nonlocal noise lookup table. Things get messy here...
+  // the outer map is a table from gate strings to gate qubit maps
+  // the gate qubit map is a map from qubits to another map of target qubits
+  // which is then the final qubit map from target qubits to error_ptr positions
+  using nonlocal_qubit_map_t = std::map<reg_t, qubit_map_t>;
+  std::unordered_map<std::string, nonlocal_qubit_map_t> nonlocal_error_table_;
+
+  // Table of single-qubit gates to use a Waltz X90 based error model
+  std::set<std::string> waltz_gates_;
+
+  // Lookup table for gate strings to enum
+  enum class Gate {id, x, y, z, h, s, sdg, t, tdg, u0, u1, u2, u3};
+  const static std::unordered_map<std::string, Gate> waltz_gate_table_;
+
+  // waltz threshold for applying u1 rotations if |theta - 2n*pi | > threshold
+  double u1_threshold_ = 1e-10;
+
+  // Rng engine
+  RngEngine rng_; // initialized with random seed
+
+  // Sample a noisy implementation of a two-X90 pulse u3 gate
+  NoiseOps sample_u3(uint_t qubit, complex_t theta,
+                     complex_t phi, complex_t lam); // TODO
+  
+  // Sample a noisy implementation of a single-X90 pulse u2 gate
+  NoiseOps sample_u2(uint_t qubit, complex_t phi, complex_t lam); // TODO
+};
+
+
+//-------------------------------------------------------------------------
+// Implementation
+//-------------------------------------------------------------------------
+
+NoiseModel::NoiseOps NoiseModel::sample_noise(const Operations::Op &op) {
+
+  // Return operator set
+  NoiseOps noise_before;
+  NoiseOps noise_after;
+  // Apply local errors first
+  if (local_errors_) {
+    // Get the qubit error map for  gate name
+    auto iter = local_error_table_.find(op.name);
+    if (iter != local_error_table_.end()) {
+      // Check if the qubits are listed in the inner model
+      const auto qubit_map = iter->second;
+      auto iter_qubits = qubit_map.find(op.qubits);
+      auto iter_default = qubit_map.find({});
+      if (iter_qubits != qubit_map.end() ||
+          iter_default != qubit_map.end()) {
+        auto &error_positions = (iter_qubits != qubit_map.end())
+          ? iter_qubits->second
+          : iter_default->second;
+        for (auto &pos : error_positions) {
+          auto ops = error_ptrs_[pos]->sample_noise(op.qubits, rng_);
+          if (error_ptrs_[pos]->errors_after())
+            noise_after.insert(noise_after.end(), ops.begin(), ops.end());
+          else
+            noise_before.insert(noise_before.end(), ops.begin(), ops.end());
+        }
+      }
+    }
+  }
+  // Apply nonlocal errors second
+  if (nonlocal_errors_) {
+    // Get the inner error map for  gate name
+    auto iter = nonlocal_error_table_.find(op.name);
+    if (iter != nonlocal_error_table_.end()) {
+      // Check if the qubits are listed in the inner model
+      const auto qubit_map = iter->second;
+      auto iter_qubits = qubit_map.find(op.qubits);
+      if (iter_qubits != qubit_map.end()) {
+        for (auto &target_pair : iter_qubits->second) {
+          auto &target_qubits = target_pair.first;
+          auto &error_positions = target_pair.second;
+          for (auto &pos : error_positions) {
+            auto ops = error_ptrs_[pos]->sample_noise(target_qubits, rng_);
+            if (error_ptrs_[pos]->errors_after())
+              noise_after.insert(noise_after.end(), ops.begin(), ops.end());
+            else
+              noise_before.insert(noise_before.end(), ops.begin(), ops.end());
+          }
+        }
+      }
+    }
+  }
+  // combine the original op with the noise ops before and after
+  noise_before.reserve(noise_before.size() + noise_after.size() + 1);
+  noise_before.push_back(op);
+  noise_before.insert(noise_before.begin(), noise_after.begin(), noise_after.end());
+  return noise_before;
+}
+
+
+Circuit NoiseModel::sample_noise(const Circuit &circ) {
+    Circuit noisy_circ = circ; // copy input circuit
+    noisy_circ.measure_sampling_flag = false; // disable measurement opt flag 
+    noisy_circ.ops.clear(); // delete ops
+    // Sample a noisy realization of the circuit
+    for (const auto &op: circ.ops) {
+      NoiseOps noisy_op = sample_noise(op);
+      // insert noisy op sequence into the circuit
+      noisy_circ.ops.insert(noisy_circ.ops.begin(), noisy_op.begin(), noisy_op.end());
+    }
+    return noisy_circ;
+}
+
+
+template <class DerivedError>
+void NoiseModel::add_local_error(const DerivedError &error,
+                           const std::vector<std::string> &op_labels,
+                           const std::vector<reg_t> &qubit_sets) {  
+  // Turn on local error flag
+  if (!op_labels.empty()) {
+    local_errors_ = true;
+  }
+  // Add error term as unique pointer                                          
+  error_ptrs_.push_back(std::unique_ptr<Noise::AbstractError>(std::make_unique<DerivedError>(error)));
+  // position is error vector is length - 1
+  const auto error_pos = error_ptrs_.size() - 1;
+  // Add error index to the error table
+  for (const auto &gate: op_labels)
+    for (const auto &qubits : qubit_sets)
+      local_error_table_[gate][qubits].push_back(error_pos); 
+}
+
+
+template <class DerivedError>
+void NoiseModel::add_nonlocal_error(const DerivedError &error,
+                                    const std::vector<std::string> &op_labels,
+                                    const std::vector<reg_t> &qubit_sets,
+                                    const std::vector<reg_t> &noise_qubit_sets) {  
+  // Turn on nonlocal error flag
+  if (!op_labels.empty() && !qubit_sets.empty() && !noise_qubit_sets.empty()) {
+    nonlocal_errors_ = true;
+  }
+  // Add error term as unique pointer                                          
+  error_ptrs_.push_back(std::unique_ptr<Noise::AbstractError>(std::make_unique<DerivedError>(error)));
+  // position is error vector is length - 1
+  const auto error_pos = error_ptrs_.size() - 1;
+  // Add error index to the error table
+  for (const auto &gate: op_labels)
+    for (const auto &qubits_gate : qubit_sets)
+      for (const auto &qubits_noise : noise_qubit_sets)
+        nonlocal_error_table_[gate][qubits_gate][qubits_noise].push_back(error_pos); 
+}
+
+
+// Currently we only support a Kraus and unitary gate errors
+// the noise config should be an array of gate_error objects:
+//  [gate_error1, gate_error2, ...]
+// where each gate_error object is of the form
+//  {
+//    "type": "gate_error",   // string
+//    "operations": ["x", "y", "z"], // list string
+//    "params": [mat1, mat2, mat3]
+//  }
+// TODO add reset errors and X90 based errors
+
+void NoiseModel::load_from_json(const json_t &js) {
+  // Check json is an array
+  if (!js.is_array()) {
+    throw std::invalid_argument("Noise params JSON is not an array");
+  }
+
+  for (const auto &gate_js : js) {
+    std::string type;
+    JSON::get_value(type, "type", gate_js);
+    std::vector<std::string> ops;
+    JSON::get_value(ops, "operations", gate_js);
+    std::vector<cmatrix_t> mats;
+    JSON::get_value(mats, "params", gate_js);
+
+    if (type == "gate_error") {
+      Noise::GateError error(mats);
+      add_local_error(error, ops);
+    } else {
+      throw std::invalid_argument("Invalid noise type (" + type + ")");
+    }
+  }
+}
+
+void from_json(const json_t &js, NoiseModel &model) {
+  model = NoiseModel(js);
+}
+
+//-------------------------------------------------------------------------
+} // end namespace Noise
+//-------------------------------------------------------------------------
+} // end namespace AER
+//-------------------------------------------------------------------------
+#endif

src/framework/operations.hpp

@@ -63,10 +63,10 @@ struct Op {
   using qubit_set_t = std::set<uint_t, std::greater<uint_t>>;
   using pauli_component_t = std::tuple<complex_t,    // component coefficient
                                        qubit_set_t,  // component qubits
-                                       std::string>; // component pauli string
+                                       std::string>; // component Pauli string
   using matrix_component_t = std::tuple<complex_t,          // component coefficient
-                                        std::vector<reg_t>, // qubit subregisters Ex: [[2], [1, 0]]
+                                        std::vector<reg_t>, // qubit sub-registers Ex: [[2], [1, 0]]
-                                        std::vector<cmatrix_t>>; // submatrices Ex: [M2, M10]
+                                        std::vector<cmatrix_t>>; // sub-matrices Ex: [M2, M10]
   std::vector<pauli_component_t> params_pauli_obs;
   std::vector<matrix_component_t> params_mat_obs; // not that diagonal matrices are stored as
                                                   // 1 x M row-matrices
@@ -97,6 +97,55 @@ inline void check_qubits(const reg_t &qubits) {
   }
 }
 
+//------------------------------------------------------------------------------
+// Generator functions
+//------------------------------------------------------------------------------
+
+Op make_mat(const reg_t &qubits, const cmatrix_t &mat, std::string label = "") {
+  Op op;
+  op.name = "mat";
+  op.qubits = qubits;
+  op.mats = {mat};
+  if (label != "")
+    op.string_params = {label};
+  return op;
+};
+
+template <typename T> // real or complex numeric type
+Op make_u1(uint_t qubit, T lam) {
+  Op op;
+  op.name = "u1";
+  op.qubits = {qubit};
+  op.params = {lam};
+  return op;
+};
+
+template <typename T> // real or complex numeric type
+Op make_u2(uint_t qubit, T phi, T lam) {
+  Op op;
+  op.name = "u2";
+  op.qubits = {qubit};
+  op.params = {phi, lam};
+  return op;
+};
+
+template <typename T> // real or complex numeric type
+Op make_u3(uint_t qubit, T theta, T phi, T lam) {
+  Op op;
+  op.name = "u3";
+  op.qubits = {qubit};
+  op.params = {theta, phi, lam};
+  return op;
+};
+
+Op make_reset(uint_t qubit, uint_t state = 0) {
+  Op op;
+  op.name = "reset";
+  op.qubits = {qubit};
+  op.params = {state};
+  return op;
+};
+
 //------------------------------------------------------------------------------
 // JSON conversion
 //------------------------------------------------------------------------------
@@ -208,10 +257,9 @@ Op json_to_op_reset(const json_t &js) {
   Op op;
   op.name = "reset";
   JSON::get_value(op.qubits, "qubits", js);
-  JSON::get_value(op.params, "params", js);
+  op.params = {0}; // default reset to 0 state
-  // If params is missing default reset to 0
+  if (JSON::check_key("params", js)) {
-  if (op.params.empty()) {
+    op.params[0] = js["params"].get<rvector_t>()[0];
-    op.params = cvector_t(op.qubits.size(), 0.);
   }
   // Validation
   check_qubits(op.qubits);
@@ -277,12 +325,18 @@ Op json_to_op_mat(const json_t &js) {
   cmatrix_t mat;
   JSON::get_value(mat, "params", js);
   op.mats.push_back(mat);
-
+  // Check for a label
+  std::string label;
+  JSON::get_value(label, "label", js);
+  if (!label.empty()) {
+    op.string_params.push_back(label);
+  }
   // Validation
   check_qubits(op.qubits);
   return op;
 }
 
+// TODO: Remove and treat as special case of mat.
 Op json_to_op_dmat(const json_t &js) {
   Op op;
   op.name = "dmat";

src/framework/utils.hpp

@@ -28,21 +28,23 @@ namespace Utils {
 class Matrix {
 public:
   // Single-qubit gates
-  const static cmatrix_t I;
+  const static cmatrix_t I;     // name: "id"
-  const static cmatrix_t X;
+  const static cmatrix_t X;     // name: "x"
-  const static cmatrix_t Y;
+  const static cmatrix_t Y;     // name: "y"
-  const static cmatrix_t Z;
+  const static cmatrix_t Z;     // name: "z"
-  const static cmatrix_t H;
+  const static cmatrix_t H;     // name: "h"
-  const static cmatrix_t S;
+  const static cmatrix_t S;     // name: "s"
-  const static cmatrix_t T;
+  const static cmatrix_t Sdg;   // name: "sdg"
-  const static cmatrix_t X90;
+  const static cmatrix_t T;     // name: "t"
+  const static cmatrix_t Tdg;   // name: "tdg"
+  const static cmatrix_t X90;   // name: "x90"
 
   // Two-qubit gates
-  const static cmatrix_t CX;
+  const static cmatrix_t CX;    // name: "cx"
-  const static cmatrix_t CZ;
+  const static cmatrix_t CZ;    // name: "cz"
-  const static cmatrix_t SWAP;
+  const static cmatrix_t SWAP;  // name: "swap"
-  const static cmatrix_t CR; // TODO
+  const static cmatrix_t CR;    // TODO
-  const static cmatrix_t CR90; // TODO
+  const static cmatrix_t CR90;  // TODO
   
   // Identity Matrix
   static cmatrix_t Identity(size_t dim);
@@ -51,6 +53,33 @@ public:
   static cmatrix_t U1(double lam);
   static cmatrix_t U2(double phi, double lam);
   static cmatrix_t U3(double theta, double phi, double lam);
+
+  // Complex arguments are implemented by taking std::real
+  // of the input
+  inline static cmatrix_t U1(complex_t lam) {return U1(std::real(lam));}
+  inline static cmatrix_t U2(complex_t phi, complex_t lam) {
+    return U2(std::real(phi), std::real(lam));
+  }
+  inline static cmatrix_t U3(complex_t theta, complex_t phi, complex_t lam) {
+    return U3(std::real(theta), std::real(phi), std::real(lam));
+  };
+
+  // Return the matrix for a named matrix string
+  // Allowed names correspond to all the const static single-qubit
+  // and two-qubit gate members
+  inline static const cmatrix_t from_name(const std::string &name) {
+    return *label_map_.at(name);
+  }
+
+  // Check if the input name string is allowed
+  inline static bool allowed_name(const std::string &name) {
+    return (label_map_.find(name) != label_map_.end());
+  }
+
+
+private:
+  // Lookup table that returns a pointer to the static data member
+  const static std::unordered_map<std::string, const cmatrix_t*> label_map_;
 };
 
 //------------------------------------------------------------------------------
@@ -202,9 +231,15 @@ const cmatrix_t Matrix::Z = make_matrix<complex_t>({{{1, 0}, {0, 0}},
 const cmatrix_t Matrix::S = make_matrix<complex_t>({{{1, 0}, {0, 0}},
                                                     {{0, 0}, {0, 1}}});
 
+const cmatrix_t Matrix::Sdg = make_matrix<complex_t>({{{1, 0}, {0, 0}},
+                                                     {{0, 0}, {0, -1}}});
+
 const cmatrix_t Matrix::T = make_matrix<complex_t>({{{1, 0}, {0, 0}},
                                                     {{0, 0}, {1 / std::sqrt(2), 1 / std::sqrt(2)}}});
 
+const cmatrix_t Matrix::Tdg = make_matrix<complex_t>({{{1, 0}, {0, 0}},
+                                                      {{0, 0}, {1 / std::sqrt(2), -1 / std::sqrt(2)}}});
+
 const cmatrix_t Matrix::H = make_matrix<complex_t>({{{1 / std::sqrt(2.), 0}, {1 / std::sqrt(2.), 0}},
                                                     {{1 / std::sqrt(2.), 0}, {-1 / std::sqrt(2.), 0}}});
 
@@ -212,23 +247,31 @@ const cmatrix_t Matrix::X90 = make_matrix<complex_t>({{{1. / std::sqrt(2.), 0},
                                                       {{0, -1. / std::sqrt(2.)}, {1. / std::sqrt(2.), 0}}});
 
 const cmatrix_t Matrix::CX = make_matrix<complex_t>({{{1, 0}, {0, 0}, {0, 0}, {0, 0}},
-                                                      {{0, 0}, {0, 0}, {0, 0}, {1, 0}},
+                                                     {{0, 0}, {0, 0}, {0, 0}, {1, 0}},
-                                                      {{0, 0}, {0, 0}, {1, 0}, {0, 0}},
+                                                     {{0, 0}, {0, 0}, {1, 0}, {0, 0}},
-                                                      {{0, 0}, {1, 0}, {0, 0}, {0, 0}}});
+                                                     {{0, 0}, {1, 0}, {0, 0}, {0, 0}}});
 
 const cmatrix_t Matrix::CZ = make_matrix<complex_t>({{{1, 0}, {0, 0}, {0, 0}, {0, 0}},
-                                                      {{0, 0}, {1, 0}, {0, 0}, {0, 0}},
+                                                     {{0, 0}, {1, 0}, {0, 0}, {0, 0}},
-                                                      {{0, 0}, {0, 0}, {1, 0}, {0, 0}},
+                                                     {{0, 0}, {0, 0}, {1, 0}, {0, 0}},
-                                                      {{0, 0}, {0, 0}, {0, 0}, {-1, 0}}});
+                                                     {{0, 0}, {0, 0}, {0, 0}, {-1, 0}}});
 
 const cmatrix_t Matrix::SWAP = make_matrix<complex_t>({{{1, 0}, {0, 0}, {0, 0}, {0, 0}},
-                                                        {{0, 0}, {0, 0}, {1, 0}, {0, 0}},
+                                                       {{0, 0}, {0, 0}, {1, 0}, {0, 0}},
-                                                        {{0, 0}, {1, 0}, {0, 0}, {0, 0}},
+                                                       {{0, 0}, {1, 0}, {0, 0}, {0, 0}},
-                                                        {{0, 0}, {0, 0}, {0, 0}, {1, 0}}});
+                                                       {{0, 0}, {0, 0}, {0, 0}, {1, 0}}});
 
 // TODO const cmatrix_t Matrix::CR = ...
 // TODO const cmatrix_t Matrix::CR90 = ...
 
+// Lookup table
+const std::unordered_map<std::string, const cmatrix_t*> Matrix::label_map_ = {
+  {"id", &Matrix::I}, {"x", &Matrix::X}, {"y", &Matrix::Y}, {"z", &Matrix::Z},
+  {"h", &Matrix::H}, {"s", &Matrix::S}, {"sdg", &Matrix::Sdg},
+  {"t", &Matrix::T}, {"tdg", &Matrix::Tdg}, {"x90", &Matrix::X90},
+  {"cx", &Matrix::CX}, {"cz", &Matrix::CZ}, {"swap", &Matrix::SWAP}
+};
+
 cmatrix_t Matrix::Identity(size_t dim) {
   cmatrix_t mat(dim, dim);
   for (size_t j=0; j<dim; j++)

src/noise/abstract_error.hpp

@@ -0,0 +1,59 @@
+/**
+ * Copyright 2018, IBM.
+ *
+ * This source code is licensed under the Apache License, Version 2.0 found in
+ * the LICENSE.txt file in the root directory of this source tree.
+ */
+
+/**
+ * @file    abstract_error.hpp
+ * @brief   Abstract Error base class for Noise model errors
+ * @author  Christopher J. Wood <cjwood@us.ibm.com>
+ */
+
+#ifndef _aer_noise_abstract_error_hpp_
+#define _aer_noise_abstract_error_hpp_
+
+#include "framework/operations.hpp"
+#include "framework/types.hpp"
+#include "framework/rng.hpp"
+
+namespace AER {
+namespace Noise {
+
+//=========================================================================
+// Error abstract base class
+//=========================================================================
+
+class AbstractError {
+public:
+
+  // Alias for return type
+  using NoiseOps = std::vector<Operations::Op>;
+
+  // Sample an realization of the error from the error model using the passed
+  // in RNG engine.
+  virtual NoiseOps sample_noise(const reg_t &qubits,
+                                RngEngine &rng) = 0;
+  
+  // Set the sampled errors to be applied after the original operation
+  inline void set_errors_after() {errors_after_op_ = true;}
+
+  // Set the sampled errors to be applied before the original operation
+  inline void set_errors_before() {errors_after_op_ = false;}
+
+  // Returns true if the errors are to be applied after the operation
+  inline bool errors_after() const {return errors_after_op_;}
+
+private:
+  // flag for where errors should be applied relative to the sampled op
+  bool errors_after_op_ = true;
+};
+
+
+//-------------------------------------------------------------------------
+} // end namespace Noise
+//-------------------------------------------------------------------------
+} // end namespace AER
+//-------------------------------------------------------------------------
+#endif

src/noise/gate_error.hpp

@@ -29,7 +29,7 @@ namespace Noise {
 // This combines unitary and Kraus errors into one error class to prevent
 // the inefficient use of unitary matrices in a Kraus decomposition.
 
-class GateError : public Error {
+class GateError : public AbstractError {
 public:
 
   GateError() = default;
@@ -40,8 +40,7 @@ public:
   //-----------------------------------------------------------------------
 
   // Sample a noisy implementation of op
-  NoiseOps sample_noise(const Operations::Op &op,
+  NoiseOps sample_noise(const reg_t &qubits,
-                        const reg_t &qubits,
                         RngEngine &rng) override;
 
   //-----------------------------------------------------------------------
@@ -81,18 +80,17 @@ protected:
 // Implementation
 //-----------------------------------------------------------------------
 
-NoiseOps GateError::sample_noise(const Operations::Op &op,
+GateError::NoiseOps GateError::sample_noise(const reg_t &qubits,
-                                 const reg_t &qubits,
+                                            RngEngine &rng) {
-                                 RngEngine &rng) {
 
   auto noise_type = rng.rand_int(probabilities_);
   switch (noise_type) {
     case 0:
-      return {op};
+      return {};
     case 1:
-      return unitary_error_.sample_noise(op, qubits, rng);
+      return unitary_error_.sample_noise(qubits, rng);
     case 2:
-      return kraus_error_.sample_noise(op, qubits, rng);
+      return kraus_error_.sample_noise(qubits, rng);
     default:
       // We shouldn't get here, but just in case...
       throw std::invalid_argument("GateError type is out of range.");

src/noise/kraus_error.hpp

@@ -14,7 +14,7 @@
 #ifndef _aer_noise_kraus_error_hpp_
 #define _aer_noise_kraus_error_hpp_
 
-#include "base/noise.hpp"
+#include "noise/abstract_error.hpp"
 
 // TODO: optimization for gate fusion with original operator matrix
 
@@ -32,7 +32,7 @@ namespace Noise {
 // of the Kraus error applying. If this is set the error map is
 // Error(rho) = (1-p) rho + p E(rho), with E defined as in (1).
 
-class KrausError : public Error {
+class KrausError : public AbstractError {
 public:
 
   //-----------------------------------------------------------------------
@@ -40,24 +40,48 @@ public:
   //-----------------------------------------------------------------------
 
   // Return the operator followed by a Kraus operation
-  NoiseOps sample_noise(const Operations::Op &op,
+  NoiseOps sample_noise(const reg_t &qubits,
-                        const reg_t &qubits,
                         RngEngine &rng) override;
 
   //-----------------------------------------------------------------------
   // Additional class methods
   //-----------------------------------------------------------------------
 
+  // Set the Kraus matrices for the error
   void set_kraus(const std::vector<cmatrix_t> &ops);
+
+  // Set the probability of the Kraus error channel being applied
+  // the probability of no error (identity) will be 1 - p_Kraus
   void set_probability(double p_kraus);
 
+  // Set the sampled errors to be applied after the original operation
+  inline void set_errors_after() {errors_after_op_ = true;}
+
+  // Set the sampled errors to be applied before the original operation
+  inline void set_errors_before() {errors_after_op_ = false;}
+
+  // Set whether the input operator should be combined into the error
+  // term (if compatible). The default is True
+  inline void combine_error(bool val) {combine_error_ = val;}
+
 protected:
   // Probabilities of Kraus error being applied
   // When sampled the outcomes correspond to:
   // 0 -> Kraus error map applied
   // 1 -> no error
   std::discrete_distribution<uint_t> probabilities_;
-  Operations::Op kraus_; // CPTP Kraus map op
+  
+  // The Kraus operation for the error channel
+  // This still needs to have the qubits it is applied to added to the Op
+  // which is done during the sample_noise method
+  Operations::Op kraus_;
+
+  // Flag for applying errors before or after the operation
+  bool errors_after_op_ = true;
+
+  // Combine error with input matrix to return a single operation
+  // if it acts on the same qubits
+  bool combine_error_ = true;
 };
 
 
@@ -66,9 +90,8 @@ protected:
 //=========================================================================
 
 // TODO add gate fusion optimization
-NoiseOps KrausError::sample_noise(const Operations::Op &op,
+KrausError::NoiseOps KrausError::sample_noise(const reg_t &qubits,
-                                  const reg_t &qubits,
+                                              RngEngine &rng) {
-                                  RngEngine &rng) {
 
   // check we have the correct number of qubits in the op for the error
   //if (qubits.size() != num_qubits_) {
@@ -76,13 +99,13 @@ NoiseOps KrausError::sample_noise(const Operations::Op &op,
   //}
   auto r = rng.rand_int(probabilities_);
   if (r == 1) {
-    // no error
+    return {}; // no error
-    return {op};
   }
-  Operations::Op err = kraus_;
+  Operations::Op error = kraus_;
-  err.qubits = qubits;
+  error.qubits = qubits;
-  return {op, err};
+  return {error};
-};
+}
+
 
 void KrausError::set_kraus(const std::vector<cmatrix_t> &ops) {
   kraus_ = Operations::Op();
@@ -90,6 +113,7 @@ void KrausError::set_kraus(const std::vector<cmatrix_t> &ops) {
   kraus_.mats = ops;
 }
 
+
 void KrausError::set_probability(double p_kraus) {
   if (p_kraus < 0 || p_kraus > 1) {
     throw std::invalid_argument("Invalid KrausError probability.");

src/noise/reset_error.hpp

@@ -0,0 +1,127 @@
+/**
+ * Copyright 2018, IBM.
+ *
+ * This source code is licensed under the Apache License, Version 2.0 found in
+ * the LICENSE.txt file in the root directory of this source tree.
+ */
+
+/**
+ * @file    reset_error.hpp
+ * @brief   Reset Error class for Qiskit-Aer simulator
+ * @author  Christopher J. Wood <cjwood@us.ibm.com>
+ */
+
+#ifndef _aer_noise_reset_error_hpp_
+#define _aer_noise_reset_error_hpp_
+
+#include "noise/unitary_error.hpp"
+
+
+namespace AER {
+namespace Noise {
+
+//=========================================================================
+// Reset Error class
+//=========================================================================
+
+// Single qubit reset error:
+// E(rho) = (1 - P) rho + \sum_j p_j |j><j|
+// where: P = sum_j p_j
+//        p_j are the probabilities of reset to each Z-basis state
+// If this error is applied to multiple qubits (eg for a multi-qubit gate)
+// The probabilities can either be independent of which qubit is being 
+// reset, or be specified for each qubit being reset.
+
+class ResetError : public AbstractError {
+public:
+
+  //-----------------------------------------------------------------------
+  // Error base required methods
+  //-----------------------------------------------------------------------
+
+  // Sample a noisy implementation of op
+  NoiseOps sample_noise(const reg_t &qubits,
+                        RngEngine &rng) override;
+
+  //-----------------------------------------------------------------------
+  // Additional class methods
+  //-----------------------------------------------------------------------
+
+  // Set the default probabilities for reset error for a qubit.
+  // The input vector should be length-2 with 
+  // probs[0] being the probability of reseting to |0> state
+  // probs[1] the probability of resetting to the |1> state
+  // The probability of no-reset error is given by 1 - p[0] - p[1]
+  inline void set_default_probabilities(const rvector_t &probs) {
+    default_probabilities_ = format_probabilities(probs);
+  }
+
+  // Set the reset error probabilities for a specific qubit.
+  // This will override the default probabilities for that qubit.
+  inline void set_qubit_probabilities(uint_t qubit, const rvector_t &probs) {
+    multi_probabilities_[qubit] = format_probabilities(probs);
+  }
+
+protected:
+  using probs_t = std::discrete_distribution<uint_t>;
+
+  // Flag for applying errors before or after the operation
+  bool errors_after_op_ = true;
+  
+  // Reset error probabilities for each qubit
+  // 0 -> no reset
+  // 1 -> reset to |0>
+  // 2 -> reset to |1>
+  probs_t default_probabilities_; 
+
+  // Map for getting reset parameters for each qubit
+  std::unordered_map<uint_t, probs_t> multi_probabilities_; 
+
+  // Helper function to convert input to internal probability format
+  probs_t format_probabilities(const rvector_t &probs);
+};
+
+//-------------------------------------------------------------------------
+// Implementation: Mixed unitary error subclass
+//-------------------------------------------------------------------------
+
+ResetError::NoiseOps ResetError::sample_noise(const reg_t &qubits,
+                                              RngEngine &rng) {
+  // Initialize return ops,
+  NoiseOps ret;
+  // Sample reset error for each qubit
+  for (const auto &q: qubits) {
+    auto it = multi_probabilities_.find(q);
+    auto r = (it == multi_probabilities_.end())
+      ? rng.rand_int(default_probabilities_)
+      : rng.rand_int(it->second);
+    if (r > 0)
+      ret.push_back(Operations::make_reset(q, r));
+  }
+  return ret;
+}
+
+
+ResetError::probs_t ResetError::format_probabilities(const rvector_t &probs) {
+  // error probability vector with [0] = p_identity
+  rvector_t error_probs = {1.}; 
+  for (const auto &p : probs) {
+    if (p < 0 || p > 1) {
+      throw std::invalid_argument("ResetError probability is not valid (p=" + std::to_string(p) +").");
+    }
+    error_probs[0] -= p;
+    error_probs.push_back(p);
+  }
+  // Check the p_identity value is valid
+  if (error_probs[0] < 0) {
+    throw std::invalid_argument("ResetError probabilities > 1.");
+  }
+  return probs_t(error_probs.begin(), error_probs.end());
+}
+
+//-------------------------------------------------------------------------
+} // end namespace Noise
+//-------------------------------------------------------------------------
+} // end namespace AER
+//-------------------------------------------------------------------------
+#endif

src/noise/unitary_error.hpp

@@ -14,7 +14,7 @@
 #ifndef _aer_noise_unitary_error_hpp_
 #define _aer_noise_unitary_error_hpp_
 
-#include "base/noise.hpp"
+#include "noise/abstract_error.hpp"
 
 // TODO optimization for gate fusion with original operator matrix
 
@@ -28,16 +28,14 @@ namespace Noise {
 // Unitary error class that can model mixed unitary error channels, and
 // coherent unitary errors.
 
-class UnitaryError : public Error {
+class UnitaryError : public AbstractError {
 public:
-
   //-----------------------------------------------------------------------
   // Error base required methods
   //-----------------------------------------------------------------------
 
   // Sample a noisy implementation of op
-  NoiseOps sample_noise(const Operations::Op &op,
+  NoiseOps sample_noise(const reg_t &qubits,
-                        const reg_t &qubits,
                         RngEngine &rng) override;
 
   //-----------------------------------------------------------------------
@@ -53,27 +51,42 @@ public:
   // Sets the unitary error matrices for the error map
   void set_unitaries(const std::vector<cmatrix_t> &mats);
 
+  // Set the sampled errors to be applied after the original operation
+  inline void set_errors_after() {errors_after_op_ = true;}
+
+  // Set the sampled errors to be applied before the original operation
+  inline void set_errors_before() {errors_after_op_ = false;}
+
+  // Set whether the input operator should be combined into the error
+  // term (if compatible). The default is True
+  inline void combine_error(bool val) {combine_error_ = val;}
+
 protected:
   // Probabilities, first entry is no-error (identity)
   std::discrete_distribution<uint_t> probabilities_;
 
   // List of unitary error matrices
   std::vector<cmatrix_t> unitaries_;
+
+  // Flag for applying errors before or after the operation
+  bool errors_after_op_ = true; 
+
+  // Combine error with input matrix to return a single operation
+  // if it acts on the same qubits
+  bool combine_error_ = true;
 };
 
 //-------------------------------------------------------------------------
 // Implementation: Mixed unitary error subclass
 //-------------------------------------------------------------------------
 
-// TODO combine terms into single matrix
+UnitaryError::NoiseOps UnitaryError::sample_noise(const reg_t &qubits,
-NoiseOps UnitaryError::sample_noise(const Operations::Op &op,
+                                                  RngEngine &rng) {
-                                    const reg_t &qubits,
-                                    RngEngine &rng) {
   auto r = rng.rand_int(probabilities_);
   // Check for no error case
   if (r == 0) {
     // no error
-    return {op};
+    return {};
   }
   // Check for invalid arguments
   if (r > unitaries_.size()) {
@@ -82,17 +95,15 @@ NoiseOps UnitaryError::sample_noise(const Operations::Op &op,
   if (unitaries_.empty()) {
     throw std::invalid_argument("Unitary error matrices are not set.");
   }
-  //if (qubits.size() != num_qubits_) {
+  // Get the error operation for the unitary
-  //  throw std::invalid_argument("Incorrect number of qubits for the error.");
-  //}
-  // TODO: combine op and error into single matrix operation
   Operations::Op error;
   error.name = "mat";
   error.mats.push_back(unitaries_[r - 1]);
   error.qubits = qubits;
-  return {op, error};
+  return {error};
 }
 
+
 void UnitaryError::set_probabilities(const rvector_t &probs) {
   rvector_t probs_with_identity({1.});
   bool probs_valid = true;

src/simulators/qubitvector/qv_state.hpp

@@ -424,7 +424,7 @@ void State::apply_op(const Operations::Op &op) {
     data_.apply_cz(op.qubits[0], op.qubits[1]);
     break;
   case Gates::reset:
-    apply_reset(op.qubits, 0);
+    apply_reset(op.qubits, uint_t(std::real(op.params[0])));
     break;
   case Gates::barrier:
     break;