Support Store and storage Var (#2028)

* support store and storage var

* add test to store cr with cr

* set qiskit 1.1.0rc1 for test

* fix lint

* add reno

---------

Co-authored-by: Jun Doi <[email protected]>

qiskit_aer/backends/aer_compiler.py

@@ -19,6 +19,7 @@
 from typing import List
 from warnings import warn
 from concurrent.futures import Executor
+import uuid
 import numpy as np
 
 from qiskit.circuit import QuantumCircuit, Clbit, ClassicalRegister, ParameterExpression
@@ -27,6 +28,7 @@
 from qiskit.circuit.library import Initialize
 from qiskit.providers.options import Options
 from qiskit.pulse import Schedule, ScheduleBlock
+from qiskit.circuit import Store
 from qiskit.circuit.controlflow import (
     WhileLoopOp,
     ForLoopOp,
@@ -61,7 +63,7 @@
 )
 
 from .backend_utils import circuit_optypes
-from ..library.control_flow_instructions import AerMark, AerJump
+from ..library.control_flow_instructions import AerMark, AerJump, AerStore
 
 
 class AerCompiler:
@@ -100,7 +102,9 @@ def compile(self, circuits, optypes=None):
                 if self._is_dynamic(circuit, compiled_optypes[idx]):
                     pm = PassManager([Decompose(["mark", "jump"])])
                     compiled_circ = pm.run(self._inline_circuit(circuit, None, None))
-
+                    # compiled_circ._vars_local = inlined_circ._vars_local
+                    # compiled_circ._vars_input = inlined_circ._vars_input
+                    # compiled_circ._vars_capture = inlined_circ._vars_capture
                     compiled_circuits.append(compiled_circ)
                     # Recompute optype for compiled circuit
                     compiled_optypes[idx] = circuit_optypes(compiled_circ)
@@ -210,6 +214,12 @@ def _inline_circuit(self, circ, continue_label, break_label, bit_map=None):
                 ret._append(
                     AerJump(continue_label, ret.num_qubits, ret.num_clbits), ret.qubits, ret.clbits
                 )
+            elif isinstance(instruction.operation, Store):
+                ret._append(
+                    AerStore(ret.num_qubits, ret.num_clbits, instruction.operation),
+                    ret.qubits,
+                    ret.clbits,
+                )
             else:
                 ret._append(instruction)
         return ret
@@ -647,7 +657,7 @@ def assemble_circuit(circuit: QuantumCircuit, basis_gates=None):
 
     num_qubits = circuit.num_qubits
     num_memory = circuit.num_clbits
-    max_conditional_idx = 0
+    extra_creg_idx = 0
 
     qreg_sizes = []
     creg_sizes = []
@@ -688,6 +698,12 @@ def assemble_circuit(circuit: QuantumCircuit, basis_gates=None):
     aer_circ.num_memory = num_memory
     aer_circ.global_phase_angle = global_phase
 
+    var_heap_map = {}
+    for var in _iter_var_recursive(circuit):
+        memory_pos = num_memory + extra_creg_idx
+        var_heap_map[var.name] = (memory_pos, var.type.width)
+        extra_creg_idx += var.type.width
+
     num_of_aer_ops = 0
     index_map = []
     for inst in circuit.data:
@@ -708,10 +724,10 @@ def assemble_circuit(circuit: QuantumCircuit, basis_gates=None):
                     if clbit in ctrl_reg:
                         mask |= 1 << idx
                         val |= ((ctrl_val >> list(ctrl_reg).index(clbit)) & 1) << idx
-            conditional_reg = num_memory + max_conditional_idx
+            conditional_reg = num_memory + extra_creg_idx
             aer_circ.bfunc(f"0x{mask:X}", f"0x{val:X}", "==", conditional_reg)
             num_of_aer_ops += 1
-            max_conditional_idx += 1
+            extra_creg_idx += 1
         elif hasattr(inst.operation, "condition_expr") and inst.operation.condition_expr:
             conditional_expr = inst.operation.condition_expr
 
@@ -740,11 +756,30 @@ def _assemble_type(expr_type):
         raise AerError(f"unknown type: {expr_type.__class__}")
 
 
+def _iter_var_recursive(circuit):
+    yield from circuit.iter_vars()
+    for data in circuit.data:
+        for param in data[0].params:
+            if isinstance(param, QuantumCircuit):
+                yield from _iter_var_recursive(param)
+
+
+def _find_var_clbits(circuit, var_uuid):
+    clbit_index = circuit.num_clbits
+    for var in _iter_var_recursive(circuit):
+        if var.var == var_uuid:
+            return list(range(clbit_index, clbit_index + var.type.width))
+        clbit_index += var.type.width
+    raise AerError(f"Var is not registed in this circuit: uuid={var_uuid}")
+
+
 def _assemble_clbit_indices(circ, c):
     if isinstance(c, (ClassicalRegister, list)):
         return [circ.find_bit(cbit).index for cbit in c]
     elif isinstance(c, Clbit):
         return [circ.find_bit(c).index]
+    elif isinstance(c, uuid.UUID):
+        return _find_var_clbits(circ, c)
     else:
         raise AerError(f"unknown clibt list: {c.__class__}")
 
@@ -981,7 +1016,18 @@ def _assemble_op(
         raise AerError(
             "control-flow instructions must be converted " f"to jump and mark instructions: {name}"
         )
-
+    elif name == "aer_store":
+        if not isinstance(operation.store.lvalue, Var):
+            raise AerError(f"unsupported lvalue : {operation.store.lvalue.__class__}")
+        aer_circ.store(qubits, _assemble_clbit_indices(circ, operation.store.lvalue.var),
+                       operation.store.rvalue.accept(_AssembleExprImpl(circ)))
+        num_of_aer_ops = 1
+    elif name == "store":
+        if not isinstance(operation.lvalue, Var):
+            raise AerError(f"unsupported lvalue : {operation.lvalue.__class__}")
+        aer_circ.store(qubits, _assemble_clbit_indices(circ, operation.lvalue.var),
+                       operation.rvalue.accept(_AssembleExprImpl(circ)))
+        num_of_aer_ops = 1
     else:
         raise AerError(f"unknown instruction: {name}")
 

qiskit_aer/backends/backend_utils.py

@@ -117,6 +117,7 @@
             "pauli",
             "mcx_gray",
             "ecr",
+            "store",
         ]
     ),
     "density_matrix": sorted(
@@ -198,6 +199,7 @@
             "diagonal",
             "initialize",
             "ecr",
+            "store",
         ]
     ),
     "stabilizer": sorted(
@@ -218,6 +220,7 @@
             "pauli",
             "ecr",
             "rz",
+            "store",
         ]
     ),
     "extended_stabilizer": sorted(
@@ -244,6 +247,7 @@
             "pauli",
             "ecr",
             "rz",
+            "store",
         ]
     ),
     "unitary": sorted(
@@ -308,6 +312,7 @@
             "multiplexer",
             "pauli",
             "ecr",
+            "store",
         ]
     ),
     "superop": sorted(
@@ -346,6 +351,7 @@
             "unitary",
             "diagonal",
             "pauli",
+            "store",
         ]
     ),
     "tensor_network": sorted(
@@ -412,6 +418,7 @@
             "pauli",
             "mcx_gray",
             "ecr",
+            "store",
         ]
     ),
 }

qiskit_aer/backends/wrappers/aer_circuit_binding.hpp

@@ -193,6 +193,7 @@ void bind_aer_circuit(MODULE m) {
   aer_circuit.def("measure", &Circuit::measure);
   aer_circuit.def("reset", &Circuit::reset);
   aer_circuit.def("set_qerror_loc", &Circuit::set_qerror_loc);
+  aer_circuit.def("store", &Circuit::store);
 }
 
 #endif

qiskit_aer/library/control_flow_instructions/__init__.py

@@ -13,3 +13,4 @@
 
 from .jump import AerJump
 from .mark import AerMark
+from .store import AerStore

qiskit_aer/library/control_flow_instructions/store.py

@@ -0,0 +1,29 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2018, 2021.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+"""
+Simulator instruction to set a program counter
+"""
+
+from qiskit.circuit import Instruction
+
+
+class AerStore(Instruction):
+    """
+    Store instruction for Aer to work wround transpilation issue
+    of qiskit.circuit.Store
+    """
+
+    _directive = True
+
+    def __init__(self, num_qubits, num_clbits, store):
+        super().__init__("aer_store", num_qubits, num_clbits, [store.lvalue, store.rvalue])
+        self.store = store

releasenotes/notes/add_var_storage-99ec3509828754d4.yaml

@@ -0,0 +1,6 @@
+---
+fixes:
+  - |
+    Enable new control flow introduced in Qiskit 1.1. ``store`` instructions are
+    newly available since Qiskit 1.1 and Aer failed if this instruction is in a circuit.
+    This PR enables ``store`` and its simulation.

src/framework/circuit.hpp

@@ -275,6 +275,11 @@ class Circuit {
     ops.push_back(Operations::make_reset(qubits, cond_regidx));
   }
 
+  void store(const reg_t &qubits, const reg_t &clbits,
+             const std::shared_ptr<Operations::CExpr> expr) {
+    ops.push_back(Operations::make_store(qubits, clbits, expr));
+  }
+
 private:
   Operations::OpSet opset_;       // Set of operation types contained in circuit
   std::set<uint_t> qubitset_;     // Set of qubits used in the circuit
@@ -450,6 +455,12 @@ void Circuit::set_params(bool truncation) {
     const auto &op = ops[rpos];
     if (op.type == OpType::mark && last_ancestor_pos == 0)
       last_ancestor_pos = rpos;
+    if (op.type == OpType::store) {
+      // Conservertively OpType::store does not allow sampling
+      can_sample = false;
+      if (last_ancestor_pos == 0)
+        last_ancestor_pos = rpos;
+    }
     if (!truncation || check_result_ancestor(op, ancestor_qubits)) {
       add_op_metadata(op);
       ancestor[rpos] = true;
@@ -571,7 +582,7 @@ void Circuit::set_params(bool truncation) {
   }
   for (size_t pos = 0; pos < head_end; ++pos) {
     if (ops_to_remove && !ancestor[pos] && ops[pos].type != OpType::mark &&
-        ops[pos].type != OpType::jump) {
+        ops[pos].type != OpType::jump && ops[pos].type != OpType::store) {
       // Skip if not ancestor
       continue;
     }
@@ -675,7 +686,8 @@ bool Circuit::check_result_ancestor(
   case OpType::save_clifford:
   case OpType::save_unitary:
   case OpType::save_mps:
-  case OpType::save_superop: {
+  case OpType::save_superop:
+  case OpType::store: {
     ancestor_qubits.insert(op.qubits.begin(), op.qubits.end());
     return true;
   }

src/framework/creg.hpp

@@ -76,6 +76,9 @@ class ClassicalRegister {
   // Apply readout error instruction to classical registers
   void apply_roerror(const Operations::Op &op, RngEngine &rng);
 
+  // Apply store instruction to classical registers
+  void apply_store(const Operations::Op &op);
+
   // Store a measurement outcome in the specified memory and register bit
   // locations
   void store_measure(const reg_t &outcome, const reg_t &memory,
@@ -258,6 +261,36 @@ void ClassicalRegister::apply_roerror(const Operations::Op &op,
   }
 }
 
+// Apply store instruction to classical registers
+void ClassicalRegister::apply_store(const Operations::Op &op) {
+  const auto &registers = op.registers;
+  const auto &expr = op.expr;
+
+  uint_t outcome = 0ULL;
+  if (expr->type->type == Operations::ValueType::Bool) {
+    outcome = op.expr->eval_bool(creg_memory_) ? 1ULL : 0ULL;
+  } else if (expr->type->type == Operations::ValueType::Uint) {
+    outcome = op.expr->eval_uint(creg_memory_);
+  }
+
+  reg_t memory;
+  reg_t memory_output;
+
+  for (size_t i = 0; i < registers.size(); i++) {
+    uint_t val = (outcome & 1ULL); // 0 or 1
+    char val_char = val ? '1' : '0';
+    outcome >>= 1;
+    if (registers[i] < creg_memory_.size()) {
+      memory.push_back(registers[i]);
+      memory_output.push_back(val);
+    }
+    const size_t pos = creg_register_.size() - registers[i] - 1;
+    creg_register_[pos] = val_char; // int->string->char
+  }
+
+  store_measure(memory_output, memory, reg_t());
+}
+
 //------------------------------------------------------------------------------
 } // end namespace AER
 //------------------------------------------------------------------------------

src/framework/operations.hpp

@@ -410,7 +410,8 @@ enum class OpType {
   jump,
   mark,
   unary_expr,
-  binary_expr
+  binary_expr,
+  store
 };
 
 enum class DataSubType {
@@ -555,6 +556,9 @@ inline std::ostream &operator<<(std::ostream &stream, const OpType &type) {
   case OpType::binary_expr:
     stream << "binary_expr";
     break;
+  case OpType::store:
+    stream << "store";
+    break;
   default:
     stream << "unknown";
   }
@@ -1033,6 +1037,17 @@ inline Op make_reset(const reg_t &qubits, const int_t conditional) {
   return op;
 }
 
+inline Op make_store(const reg_t &qubits, const reg_t &clbits,
+                     const std::shared_ptr<CExpr> expr) {
+  Op op;
+  op.type = OpType::store;
+  op.name = "store";
+  op.qubits = qubits;
+  op.registers = clbits;
+  op.expr = expr;
+  return op;
+}
+
 inline Op make_multiplexer(const reg_t &qubits,
                            const std::vector<cmatrix_t> &mats,
                            const int_t conditional = -1,

src/simulators/circuit_executor.hpp

@@ -1324,7 +1324,8 @@ bool Executor<state_t>::check_measure_sampling_opt(const Circuit &circ) const {
       circ.opset().contains(Operations::OpType::kraus) ||
       circ.opset().contains(Operations::OpType::superop) ||
       circ.opset().contains(Operations::OpType::jump) ||
-      circ.opset().contains(Operations::OpType::mark)) {
+      circ.opset().contains(Operations::OpType::mark) ||
+      circ.opset().contains(Operations::OpType::store)) {
     return false;
   }
   // Otherwise true

src/simulators/density_matrix/densitymatrix_state.hpp

@@ -45,7 +45,7 @@ const Operations::OpSet StateOpSet(
      OpType::set_densmat,  OpType::save_expval, OpType::save_expval_var,
      OpType::save_densmat, OpType::save_probs,  OpType::save_probs_ket,
      OpType::save_amps_sq, OpType::save_state,  OpType::jump,
-     OpType::mark},
+     OpType::mark,         OpType::store},
     // Gates
     {"U",   "CX", "u1",   "u2",  "u3",    "u",     "cx",  "cy",  "cz",  "swap",
      "id",  "x",  "y",    "z",   "h",     "s",     "sdg", "t",   "tdg", "ccx",