Refactor and improve fusion (#255)

* refactor fusion by deleting matrix_sequence and fix bugs in mcu in qubitvector.

* CH supports nop

* specialize u1 and cu1 fusion

* refactor fusion-related codes

* remove nop operation type

src/base/controller.hpp

@@ -194,7 +194,7 @@ class Controller {
 
   // Generate an equivalent circuit with input_circ as output_circ.
   template <class state_t>
-  Circuit optimize_circuit(const Circuit &input_circ,
+  void optimize_circuit(Circuit &input_circ,
                            state_t& state,
                            OutputData &data) const;
 
@@ -477,20 +477,18 @@ bool Controller::validate_memory_requirements(state_t &state,
 // Circuit optimization
 //-------------------------------------------------------------------------
 template <class state_t>
-Circuit Controller::optimize_circuit(const Circuit &input_circ,
+void Controller::optimize_circuit(Circuit &input_circ,
                                      state_t& state,
                                      OutputData &data) const {
 
-  Circuit working_circ = input_circ;
   Operations::OpSet allowed_opset;
   allowed_opset.optypes = state.allowed_ops();
   allowed_opset.gates = state.allowed_gates();
   allowed_opset.snapshots = state.allowed_snapshots();
 
-  for (std::shared_ptr<Transpile::CircuitOptimization> opt: optimizations_)
-    opt->optimize_circuit(working_circ, allowed_opset, data);
-
-  return working_circ;
+  for (std::shared_ptr<Transpile::CircuitOptimization> opt: optimizations_) {
+    opt->optimize_circuit(input_circ, allowed_opset, data);
+  }
 }
 
 //-------------------------------------------------------------------------

src/framework/operations.hpp

@@ -36,7 +36,7 @@ enum class RegComparison {Equal, NotEqual, Less, LessEqual, Greater, GreaterEqua
 // Enum class for operation types
 enum class OpType {
   gate, measure, reset, bfunc, barrier, snapshot,
-  matrix, matrix_sequence, multiplexer, kraus, roerror, noise_switch, initialize
+  matrix, multiplexer, kraus, roerror, noise_switch, initialize
 };
 
 inline std::ostream& operator<<(std::ostream& stream, const OpType& type) {
@@ -62,9 +62,6 @@ inline std::ostream& operator<<(std::ostream& stream, const OpType& type) {
   case OpType::matrix:
     stream << "matrix";
     break;
-  case OpType::matrix_sequence:
-    stream << "matrix_sequence";
-    break;
   case OpType::multiplexer:
     stream << "multiplexer";
     break;
@@ -138,6 +135,20 @@ inline std::ostream& operator<<(std::ostream& s, const Op& op) {
     s << qubit;
     first = false;
   }
+  s << "],[";
+  first = true;
+  for (reg_t reg: op.regs) {
+    if (!first) s << ",";
+    s << "[";
+    bool first0 = true;
+    for (size_t qubit: reg) {
+      if (!first0) s << ",";
+      s << qubit;
+      first0 = false;
+    }
+    s << "]";
+    first = false;
+  }
   s << "]";
   return s;
 }
@@ -390,14 +401,15 @@ inline Op make_unitary(const reg_t &qubits, const cmatrix_t &mat, std::string la
   return op;
 }
 
-inline Op make_matrix_sequence(const std::vector<reg_t> &regs, const std::vector<cmatrix_t> &mats, std::string label = "") {
+inline Op make_fusion(const reg_t &qubits, const cmatrix_t &mat, const std::vector<Op>& fusioned_ops, std::string label = "") {
   Op op;
-  op.type = OpType::matrix_sequence;
-  op.name = "matrix_sequence";
-  op.regs = regs;
-  op.mats = mats;
+  op.type = OpType::matrix;
+  op.name = "fusion";
+  op.qubits = qubits;
+  op.mats = {mat};
   if (label != "")
     op.string_params = {label};
+
   return op;
 }
 

src/framework/utils.hpp

@@ -92,6 +92,7 @@ class Matrix {
 template<class T> matrix<T> make_matrix(const std::vector<std::vector<T>> &mat);
 template<class T> matrix<T> devectorize_matrix(const std::vector<T> &vec);
 template<class T> std::vector<T> vectorize_matrix(const matrix<T> &mat);
+template<class T> std::vector<T> vectorize_diagonal_matrix(const matrix<T>& mat);
 
 // Transformations
 template <class T> matrix<T> transpose(const matrix<T> &A);
@@ -496,7 +497,6 @@ matrix<T> devectorize_matrix(const std::vector<T>& vec) {
   return mat;
 }
 
-
 template<class T>
 std::vector<T> vectorize_matrix(const matrix<T>& mat) {
   std::vector<T> vec;
@@ -510,6 +510,15 @@ std::vector<T> vectorize_matrix(const matrix<T>& mat) {
   return vec;
 }
 
+template<class T>
+std::vector<T> vectorize_diagonal_matrix(const matrix<T>& mat) {
+  std::vector<T> vec;
+  size_t size = std::min(mat.GetRows(), mat.GetColumns());
+  vec.resize(size, 0.);
+  for (size_t i=0; i < size; i++)
+    vec[i] = mat(i, i);
+  return vec;
+}
 
 template <class T>
 matrix<T> make_matrix(const std::vector<std::vector<T>> & mat) {

src/simulators/qasm/qasm_controller.hpp

@@ -251,7 +251,7 @@ class QasmController : public Base::Controller {
 // Constructor
 //-------------------------------------------------------------------------
 QasmController::QasmController() {
-  add_circuit_optimization(Transpile::ReduceNop());
+  add_circuit_optimization(Transpile::ReduceBarrier());
   add_circuit_optimization(Transpile::Fusion());
   add_circuit_optimization(Transpile::TruncateQubits());
 }
@@ -501,7 +501,7 @@ void QasmController::run_circuit_with_noise(const Circuit &circ,
   while(shots-- > 0) {
     Circuit noise_circ = noise_model_.sample_noise(circ, rng);
     if (noise_circ.num_qubits > circuit_opt_noise_threshold_) {
-      noise_circ = optimize_circuit(noise_circ, state, data);
+      optimize_circuit(noise_circ, state, data);
     }
     run_single_shot(noise_circ, state, initial_state, data, rng);
   }                                   
@@ -518,7 +518,7 @@ void QasmController::run_circuit_without_noise(const Circuit &circ,
   // Optimize circuit for state type
   Circuit opt_circ = circ;
   if (circ.num_qubits > circuit_opt_ideal_threshold_) {
-    opt_circ = optimize_circuit(opt_circ, state, data);
+    optimize_circuit(opt_circ, state, data);
   }
 
   // Check if measure sampler and optimization are valid