diff --git a/qiskit/providers/aer/backends/qasm_simulator.py b/qiskit/providers/aer/backends/qasm_simulator.py index c080785464..57805d581e 100644 --- a/qiskit/providers/aer/backends/qasm_simulator.py +++ b/qiskit/providers/aer/backends/qasm_simulator.py @@ -16,7 +16,7 @@ import logging from math import log2 from qiskit.util import local_hardware_info -from qiskit.providers.models import BackendConfiguration +from qiskit.providers.models import QasmBackendConfiguration from .aerbackend import AerBackend # pylint: disable=import-error from .qasm_controller_wrapper import qasm_controller_execute @@ -173,7 +173,7 @@ class QasmSimulator(AerBackend): def __init__(self, configuration=None, provider=None): super().__init__( qasm_controller_execute, - BackendConfiguration.from_dict(self.DEFAULT_CONFIGURATION), + QasmBackendConfiguration.from_dict(self.DEFAULT_CONFIGURATION), provider=provider) def _validate(self, qobj, backend_options, noise_model): diff --git a/qiskit/providers/aer/backends/statevector_simulator.py b/qiskit/providers/aer/backends/statevector_simulator.py index 7afbd85f59..98bb4e216b 100644 --- a/qiskit/providers/aer/backends/statevector_simulator.py +++ b/qiskit/providers/aer/backends/statevector_simulator.py @@ -19,7 +19,7 @@ import logging from math import log2 from qiskit.util import local_hardware_info -from qiskit.providers.models import BackendConfiguration +from qiskit.providers.models import QasmBackendConfiguration from .aerbackend import AerBackend # pylint: disable=import-error from .statevector_controller_wrapper import statevector_controller_execute @@ -99,7 +99,7 @@ class StatevectorSimulator(AerBackend): def __init__(self, configuration=None, provider=None): super().__init__(statevector_controller_execute, - BackendConfiguration.from_dict(self.DEFAULT_CONFIGURATION), + QasmBackendConfiguration.from_dict(self.DEFAULT_CONFIGURATION), provider=provider) def _validate(self, qobj, backend_options, noise_model): diff --git a/qiskit/providers/aer/backends/unitary_simulator.py b/qiskit/providers/aer/backends/unitary_simulator.py index 7abcf113a4..b9286368ba 100644 --- a/qiskit/providers/aer/backends/unitary_simulator.py +++ b/qiskit/providers/aer/backends/unitary_simulator.py @@ -19,7 +19,7 @@ import logging from math import log2, sqrt from qiskit.util import local_hardware_info -from qiskit.providers.models import BackendConfiguration +from qiskit.providers.models import QasmBackendConfiguration from .aerbackend import AerBackend from ..aererror import AerError @@ -105,7 +105,7 @@ class UnitarySimulator(AerBackend): def __init__(self, configuration=None, provider=None): super().__init__(unitary_controller_execute, - BackendConfiguration.from_dict(self.DEFAULT_CONFIGURATION), + QasmBackendConfiguration.from_dict(self.DEFAULT_CONFIGURATION), provider=provider) def _validate(self, qobj, backend_options, noise_model): diff --git a/src/noise/abstract_error.hpp b/src/noise/abstract_error.hpp deleted file mode 100644 index 5e88848f29..0000000000 --- a/src/noise/abstract_error.hpp +++ /dev/null @@ -1,76 +0,0 @@ -/** - * This code is part of Qiskit. - * - * (C) Copyright IBM 2018, 2019. - * - * 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. - */ - -#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: - - // Constructors - AbstractError() = default; - virtual ~AbstractError() = default; - - // Alias for return type - using NoiseOps = std::vector; - - // 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) const = 0; - - // Load from a JSON file - virtual void load_from_json(const json_t &js) = 0; - - // Set number of qubits or memory bits for error - inline void set_num_qubits(uint_t num_qubits) {num_qubits_ = num_qubits;} - - // Get number of qubits or memory bits for error - inline uint_t get_num_qubits() const {return num_qubits_;} - - // 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; - - // Number of qubits / memory bits the error applies to - uint_t num_qubits_ = 0; -}; - - -//------------------------------------------------------------------------- -} // end namespace Noise -//------------------------------------------------------------------------- -} // end namespace AER -//------------------------------------------------------------------------- -#endif \ No newline at end of file diff --git a/src/noise/noise_model.hpp b/src/noise/noise_model.hpp index 0ad76c3c65..adaf1acb34 100644 --- a/src/noise/noise_model.hpp +++ b/src/noise/noise_model.hpp @@ -22,9 +22,6 @@ #include "framework/types.hpp" #include "framework/rng.hpp" #include "framework/circuit.hpp" -#include "noise/abstract_error.hpp" - -// For JSON parsing of specific error types #include "noise/quantum_error.hpp" #include "noise/readout_error.hpp" @@ -44,6 +41,7 @@ namespace Noise { class NoiseModel { public: + using Method = QuantumError::Method; using NoiseOps = std::vector; NoiseModel() = default; @@ -51,8 +49,18 @@ class NoiseModel { // Sample a noisy implementation of a full circuit // An RngEngine is passed in as a reference so that sampling - // can be done in a thread-safe manner - Circuit sample_noise(const Circuit &circ, RngEngine &rng) const; + // can be done in a thread-safe manner. + // Sample methods are: + // standard: each noisy op will be returned along with additional noise ops + // superop: each noisy gate or reset will be returned as a single superop + Circuit sample_noise(const Circuit &circ, + RngEngine &rng) const; + + // Set sample mode to superoperator + // This will cause all QuantumErrors stored in the noise model + // to calculate their superoperator representations and raise + // an exception if they cannot be converted. + void activate_superop_method(); //----------------------------------------------------------------------- // Checking if errors types are in noise model @@ -99,7 +107,7 @@ class NoiseModel { // Add a ReadoutError to the noise model void add_readout_error(const ReadoutError &error, const std::vector &op_qubits = {}); - + // Set which single qubit gates should use the X90 waltz error model inline void set_x90_gates(const stringset_t &x90_gates) { x90_gates_ = x90_gates; @@ -130,8 +138,9 @@ class NoiseModel { private: - // Sample noise for the current operation - NoiseOps sample_noise(const Operations::Op &op, RngEngine &rng) const; + // Sample noise for the current operation. + NoiseOps sample_noise(const Operations::Op &op, + RngEngine &rng) const; // Sample noise for the current operation void sample_readout_noise(const Operations::Op &op, @@ -144,9 +153,9 @@ class NoiseModel { RngEngine &rng) const; void sample_nonlocal_quantum_noise(const Operations::Op &op, - NoiseOps &noise_before, + NoiseOps &noise_ops, NoiseOps &noise_after, - RngEngine &rng) const; + RngEngine &rng) const; // Sample noise for the current operation NoiseOps sample_noise_helper(const Operations::Op &op, @@ -204,6 +213,14 @@ class NoiseModel { std::string remap_string(const std::string key, const std::unordered_map &mapping) const; + // Helper function to try and convert an instruciton to superop matrix + // If conversion isn't possible this returns an empty matrix + cmatrix_t op2superop(const Operations::Op &op) const; + + // Try and convert an instruciton to unitary matrix + // If conversion isn't possible this returns an empty matrix + cmatrix_t op2unitary(const Operations::Op &op) const; + // Table of single-qubit gates to use a Waltz X90 based error model stringset_t x90_gates_; @@ -216,11 +233,14 @@ class NoiseModel { // Joint OpSet of all errors Operations::OpSet opset_; + + // Sampling method + Method method_ = Method::standard; }; //========================================================================= -// Noise Model class +// Noise sampling //========================================================================= NoiseModel::NoiseOps NoiseModel::sample_noise(const Operations::Op &op, @@ -236,9 +256,13 @@ NoiseModel::NoiseOps NoiseModel::sample_noise(const Operations::Op &op, if (gate != waltz_gate_table_.end()) { switch (gate->second) { case WaltzGate::u3: - return sample_noise_x90_u3(op.qubits[0], op.params[0], op.params[1], op.params[2], rng); + return sample_noise_x90_u3(op.qubits[0], + op.params[0], op.params[1], op.params[2], + rng); case WaltzGate::u2: - return sample_noise_x90_u2(op.qubits[0], op.params[0], op.params[1], rng); + return sample_noise_x90_u2(op.qubits[0], + op.params[0], op.params[1], + rng); case WaltzGate::x: return sample_noise_x90_u3(op.qubits[0], M_PI, 0., M_PI, rng); case WaltzGate::y: @@ -256,7 +280,8 @@ NoiseModel::NoiseOps NoiseModel::sample_noise(const Operations::Op &op, } -Circuit NoiseModel::sample_noise(const Circuit &circ, RngEngine &rng) const { +Circuit NoiseModel::sample_noise(const Circuit &circ, + RngEngine &rng) const { bool noise_active = true; // set noise active to on-state Circuit noisy_circ = circ; // copy input circuit noisy_circ.measure_sampling_flag = false; // disable measurement opt flag @@ -275,6 +300,9 @@ Circuit NoiseModel::sample_noise(const Circuit &circ, RngEngine &rng) const { case Operations::OpType::kraus: noisy_circ.ops.push_back(op); break; + case Operations::OpType::superop: + noisy_circ.ops.push_back(op); + break; case Operations::OpType::roerror: noisy_circ.ops.push_back(op); break; @@ -297,8 +325,18 @@ Circuit NoiseModel::sample_noise(const Circuit &circ, RngEngine &rng) const { } +void NoiseModel::activate_superop_method() { + // Set internal sampling method + method_ = Method::superop; + // Compute superoperators + for (auto& qerror : quantum_errors_) { + qerror.compute_superoperator(); + } +} + + void NoiseModel::add_readout_error(const ReadoutError &error, - const std::vector &op_qubits) { + const std::vector &op_qubits) { // Add roerror to noise model ops opset_.optypes.insert(Operations::OpType::roerror); // Add error term as unique pointer @@ -402,10 +440,51 @@ NoiseModel::NoiseOps NoiseModel::sample_noise_helper(const Operations::Op &op, sample_readout_noise(op, noise_after, rng); } - // combine the original op with the noise ops before and after + // Combine errors noise_before.reserve(noise_before.size() + noise_after.size() + 1); noise_before.push_back(op); noise_before.insert(noise_before.end(), noise_after.begin(), noise_after.end()); + if (op.type != Operations::OpType::measure && + noise_before.size() == 2 && + noise_before[0].qubits == noise_before[1].qubits) { + // Try and fuse operations + // If either are superoperators combine superoperators + // else if either are unitaries combine unitaries + // otherwise return the full list + auto& first_op = noise_before[0]; + auto& second_op = noise_before[1]; + + if (second_op.type == Operations::OpType::superop) { + auto& current = second_op; + const auto mat = op2superop(first_op); + if (!mat.empty()) { + current.mats[0] = current.mats[0] * mat; + return NoiseOps({current}); + } + } else if (first_op.type == Operations::OpType::superop) { + auto& current = first_op; + const auto mat = op2superop(second_op); + if (!mat.empty()) { + current.mats[0] = mat * current.mats[0]; + return NoiseOps({current}); + } + } else if (second_op.type == Operations::OpType::matrix) { + auto& current = noise_before[1]; + const auto mat = op2unitary(first_op); + if (!mat.empty()) { + current.mats[0] = current.mats[0] * mat; + return NoiseOps({current}); + } + } else if (first_op.type == Operations::OpType::matrix) { + auto& current = first_op; + const auto mat = op2unitary(second_op); + if (!mat.empty()) { + current.mats[0] = mat * current.mats[0]; + return NoiseOps({current}); + } + } + } + // Otherwise return the list of ops return noise_before; } @@ -530,7 +609,8 @@ void NoiseModel::sample_local_quantum_noise(const Operations::Op &op, ? iter_qubits->second : iter_default->second; for (auto &pos : error_positions) { - auto noise_ops = quantum_errors_[pos].sample_noise(string2reg(qubit_keys[qs]), rng); + auto noise_ops = quantum_errors_[pos].sample_noise(string2reg(qubit_keys[qs]), rng, + method_); // Duplicate same sampled error operations if (quantum_errors_[pos].errors_after()) noise_after.insert(noise_after.end(), noise_ops.begin(), noise_ops.end()); @@ -587,7 +667,8 @@ void NoiseModel::sample_nonlocal_quantum_noise(const Operations::Op &op, auto &target_qubits = target_pair.first; auto &error_positions = target_pair.second; for (auto &pos : error_positions) { - auto ops = quantum_errors_[pos].sample_noise(string2reg(target_qubits), rng); + auto ops = quantum_errors_[pos].sample_noise(string2reg(target_qubits), rng, + method_); if (quantum_errors_[pos].errors_after()) noise_after.insert(noise_after.end(), ops.begin(), ops.end()); else @@ -610,43 +691,167 @@ NoiseModel::waltz_gate_table_ = { NoiseModel::NoiseOps NoiseModel::sample_noise_x90_u3(uint_t qubit, - complex_t theta, - complex_t phi, - complex_t lambda, - RngEngine &rng) const { - NoiseOps ret; + complex_t theta, + complex_t phi, + complex_t lambda, + RngEngine &rng) const { + // sample noise for single X90 const auto x90 = Operations::make_unitary({qubit}, Utils::Matrix::X90, "x90"); - if (std::abs(lambda) > u1_threshold_ - && std::abs(lambda - 2 * M_PI) > u1_threshold_ - && std::abs(lambda + 2 * M_PI) > u1_threshold_) - ret.push_back(Operations::make_u1(qubit, lambda)); // add 1st U1 - auto sample = sample_noise_helper(x90, rng); // sample noise for 1st X90 - ret.insert(ret.end(), sample.begin(), sample.end()); // add 1st noisy X90 - if (std::abs(theta + M_PI) > u1_threshold_ - && std::abs(theta - M_PI) > u1_threshold_) - ret.push_back(Operations::make_u1(qubit, theta + M_PI)); // add 2nd U1 - sample = sample_noise_helper(x90, rng); // sample noise for 2nd X90 - ret.insert(ret.end(), sample.begin(), sample.end()); // add 2nd noisy X90 - if (std::abs(phi + M_PI) > u1_threshold_ - && std::abs(phi - M_PI) > u1_threshold_) - ret.push_back(Operations::make_u1(qubit, phi + M_PI)); // add 3rd U1 - return ret; + switch (method_) { + case Method::superop: { + // The first element of the sample should be the superoperator to combine + auto sample = sample_noise_helper(x90, rng); + // The first element of the sample should be the superoperator to combine + if (sample[0].type != Operations::OpType::superop) { + throw std::runtime_error("Sampling superoperator noise failed."); + } + cmatrix_t& current = sample[0].mats[0]; + // Combine with middle u1 gate with two noisy x90 superops + auto mat = Utils::Matrix::u1(theta + M_PI); + auto super = Utils::tensor_product(AER::Utils::conjugate(mat), mat); + current = current * super * current; + + // Prepend with first u1 matrix with superop + mat = Utils::Matrix::u1(lambda); + super = Utils::tensor_product(AER::Utils::conjugate(mat), + mat); + current = current * super; + // Append third u1 matrix with superop + mat = Utils::Matrix::u1(phi + M_PI); + super = Utils::tensor_product(AER::Utils::conjugate(mat), mat); + current = super * current; + return sample; + } + default: { + NoiseOps ret; + if (std::abs(lambda) > u1_threshold_ + && std::abs(lambda - 2 * M_PI) > u1_threshold_ + && std::abs(lambda + 2 * M_PI) > u1_threshold_) { + ret.push_back(Operations::make_u1(qubit, lambda)); // add 1st U1 + } + auto sample = sample_noise_helper(x90, rng); // sample noise for 1st X90 + ret.insert(ret.end(), sample.begin(), sample.end()); // add 1st noisy X90 + if (std::abs(theta + M_PI) > u1_threshold_ + && std::abs(theta - M_PI) > u1_threshold_) { + ret.push_back(Operations::make_u1(qubit, theta + M_PI)); // add 2nd U1 + } + sample = sample_noise_helper(x90, rng); // sample noise for 2nd X90 + ret.insert(ret.end(), sample.begin(), sample.end()); // add 2nd noisy X90 + if (std::abs(phi + M_PI) > u1_threshold_ + && std::abs(phi - M_PI) > u1_threshold_) { + ret.push_back(Operations::make_u1(qubit, phi + M_PI)); // add 3rd U1 + } + return ret; + } + } } NoiseModel::NoiseOps NoiseModel::sample_noise_x90_u2(uint_t qubit, - complex_t phi, - complex_t lambda, - RngEngine &rng) const { - NoiseOps ret; + complex_t phi, + complex_t lambda, + RngEngine &rng) const { + // sample noise for single X90 const auto x90 = Operations::make_unitary({qubit}, Utils::Matrix::X90, "x90"); - if (std::abs(lambda - 0.5 * M_PI) > u1_threshold_) - ret.push_back(Operations::make_u1(qubit, lambda - 0.5 * M_PI)); // add 1st U1 - auto sample = sample_noise_helper(x90, rng); // sample noise for 1st X90 - ret.insert(ret.end(), sample.begin(), sample.end()); // add 1st noisy X90 - if (std::abs(phi + 0.5 * M_PI) > u1_threshold_) - ret.push_back(Operations::make_u1(qubit, phi + 0.5 * M_PI)); // add 2nd U1 - return ret; + auto sample = sample_noise_helper(x90, rng); + switch (method_) { + case Method::superop: { + // The first element of the sample should be the superoperator to combine + if (sample[0].type != Operations::OpType::superop) { + throw std::runtime_error("Sampling superoperator noise failed."); + } + cmatrix_t ¤t = sample[0].mats[0]; + // Combine first u1 matrix with superop + auto mat = Utils::Matrix::u1(lambda - 0.5 * M_PI); + auto super = Utils::tensor_product(AER::Utils::conjugate(mat), + mat); + current = current * super; + // Combine second u1 matrix with superop + mat = Utils::Matrix::u1(phi + 0.5 * M_PI); + super = Utils::tensor_product(AER::Utils::conjugate(mat), mat); + current = super * current; + return sample; + } + default: { + NoiseOps ret; + // Standard method doesn't combine any ops + if (std::abs(lambda - 0.5 * M_PI) > u1_threshold_) { + // add 1st u1 + ret.push_back(Operations::make_u1(qubit, lambda - 0.5 * M_PI)); + } + // add 1st noisy x90 + ret.insert(ret.end(), sample.begin(), sample.end()); + if (std::abs(phi + 0.5 * M_PI) > u1_threshold_) { + // add 2nd u1 + ret.push_back(Operations::make_u1(qubit, phi + 0.5 * M_PI)); + } + return ret; + } + } +} + + +cmatrix_t NoiseModel::op2superop(const Operations::Op &op) const { + switch (op.type) { + case Operations::OpType::superop: + return op.mats[0]; + case Operations::OpType::kraus: { + const auto dim = op.mats[0].GetRows(); + cmatrix_t mat(dim * dim, dim * dim); + for (const auto& kraus : op.mats) { + mat += Utils::unitary_superop(kraus); + } + return mat; + } + case Operations::OpType::reset: + return Utils::SMatrix::reset(1ULL << op.qubits.size()); + case Operations::OpType::matrix: + return Utils::unitary_superop(op.mats[0]); + case Operations::OpType::gate: { + // Check if a parameterized gate + if (op.name == "u1") { + return Utils::SMatrix::u1(op.params[0]); + } + if (op.name == "u2") { + return Utils::SMatrix::u2(op.params[0], op.params[1]); + } + if (op.name == "u3") { + return Utils::SMatrix::u3(op.params[0], op.params[1], op.params[2]); + } + if (Utils::SMatrix::allowed_name(op.name)) { + // Check if we can convert this gate to a standard superoperator matrix + return Utils::SMatrix::from_name(op.name); + } + } + default: + return cmatrix_t(); + } +} + + +cmatrix_t NoiseModel::op2unitary(const Operations::Op &op) const { + switch (op.type) { + case Operations::OpType::matrix: + return op.mats[0]; + case Operations::OpType::gate: { + // Check if a parameterized gate + if (op.name == "u1") { + return Utils::SMatrix::u1(op.params[0]); + } + if (op.name == "u2") { + return Utils::SMatrix::u2(op.params[0], op.params[1]); + } + if (op.name == "u3") { + return Utils::SMatrix::u3(op.params[0], op.params[1], op.params[2]); + } + if (Utils::SMatrix::allowed_name(op.name)) { + // Check if we can convert this gate to a standard superoperator matrix + return Utils::SMatrix::from_name(op.name); + } + } + default: + return cmatrix_t(); + } } diff --git a/src/noise/quantum_error.hpp b/src/noise/quantum_error.hpp index a06be42870..f466ccfd10 100644 --- a/src/noise/quantum_error.hpp +++ b/src/noise/quantum_error.hpp @@ -15,7 +15,7 @@ #ifndef _aer_noise_quantum_error_hpp_ #define _aer_noise_quantum_error_hpp_ -#include "noise/abstract_error.hpp" +#include "simulators/superoperator/superoperator_state.hpp" namespace AER { namespace Noise { @@ -27,27 +27,39 @@ namespace Noise { // Quantum error class that can model any error that is expressed as a // qobj instruction acting on qubits. -class QuantumError : public AbstractError { +class QuantumError { public: + // Methods for sampling + enum class Method {standard, superop}; + // Alias for return type using NoiseOps = std::vector; //----------------------------------------------------------------------- - // Error base required methods + // Sampling method //----------------------------------------------------------------------- // Sample a noisy implementation of op NoiseOps sample_noise(const reg_t &qubits, - RngEngine &rng) const override; + RngEngine &rng, + Method method = Method::standard) const; - // Load a QuantumError object from a JSON Error object - void load_from_json(const json_t &js) override; + // Return the opset for the quantum error + const Operations::OpSet& opset() const {return opset_;} + + // Return the superoperator matrix representation of the error. + // If the error cannot be converted to a superoperator and error + // will be raised. + const cmatrix_t& superoperator() const; //----------------------------------------------------------------------- - // Additional class methods + // Initialization //----------------------------------------------------------------------- + // Load a QuantumError object from a JSON Error object + void load_from_json(const json_t &js); + // Sets the sub-circuits and probabilities to be sampled from. // The length of the circuits vector and probability vector must be equal. void set_circuits(const std::vector &circuits, @@ -58,15 +70,37 @@ class QuantumError : public AbstractError { // non-kraus subcircuits. void set_from_kraus(const std::vector &mats); + // Compute the superoperator representation of the quantum error + void compute_superoperator(); + + //----------------------------------------------------------------------- + // Utility + //----------------------------------------------------------------------- + + // Set number of qubits or memory bits for error + inline void set_num_qubits(uint_t num_qubits) {num_qubits_ = num_qubits;} + + // Get number of qubits or memory bits for error + inline uint_t get_num_qubits() const {return num_qubits_;} + + // 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_;} + // Set threshold for checking probabilities and matrices void set_threshold(double); - const Operations::OpSet& opset() const {return opset_;} - protected: + // Number of qubits sthe error applies to + uint_t num_qubits_ = 0; + // Probabilities, first entry is no-error (identity) rvector_t probabilities_; - //std::discrete_distribution probabilities_; // List of unitary error matrices std::vector circuits_; @@ -76,37 +110,56 @@ class QuantumError : public AbstractError { // threshold for validating if matrices are unitary double threshold_ = 1e-10; + + // Superoperator matrix representation of the error + cmatrix_t superoperator_; + + // flag for where errors should be applied relative to the sampled op + bool errors_after_op_ = true; }; //------------------------------------------------------------------------- // Implementation: Mixed unitary error subclass //------------------------------------------------------------------------- + QuantumError::NoiseOps QuantumError::sample_noise(const reg_t &qubits, - RngEngine &rng) const { + RngEngine &rng, + Method method) const { if (qubits.size() < get_num_qubits()) { std::stringstream msg; msg << "QuantumError: qubits size (" << qubits.size() << ")"; msg << " < error qubits (" << get_num_qubits() << ")."; throw std::invalid_argument(msg.str()); } - auto r = rng.rand_int(probabilities_); - // Check for invalid arguments - if (r + 1 > circuits_.size()) { - std::stringstream msg; - msg << "QuantumError: probability outcome (" << r << ")"; - msg << " is greater than number of circuits (" << circuits_.size() << ")."; - throw std::invalid_argument(msg.str()); - } - NoiseOps noise_ops = circuits_[r]; - // Add qubits to noise op commands; - for (auto &op : noise_ops) { - // Update qubits based on position in qubits list - for (size_t q=0; q < op.qubits.size(); q++) { - op.qubits[q] = qubits[op.qubits[q]]; + switch (method) { + case Method::superop: { + // Truncate qubits to size of the actual error + reg_t op_qubits = qubits; + op_qubits.resize(get_num_qubits()); + auto op = Operations::make_superop(op_qubits, superoperator()); + return NoiseOps({op}); + } + default: { + auto r = rng.rand_int(probabilities_); + // Check for invalid arguments + if (r + 1 > circuits_.size()) { + throw std::invalid_argument( + "QuantumError: probability outcome (" + std::to_string(r) + ")" + ") is greater than number of circuits (" + std::to_string(circuits_.size()) + ")." + ); + } + NoiseOps noise_ops = circuits_[r]; + // Add qubits to noise op commands; + for (auto &op : noise_ops) { + // Update qubits based on position in qubits list + for (size_t q=0; q < op.qubits.size(); q++) { + op.qubits[q] = qubits[op.qubits[q]]; + } + } + return noise_ops; } } - return noise_ops; } void QuantumError::set_threshold(double threshold) { @@ -116,11 +169,11 @@ void QuantumError::set_threshold(double threshold) { void QuantumError::set_circuits(const std::vector &circuits, const rvector_t &probs) { if (probs.size() != circuits.size()) { - std::stringstream msg; - msg << "QuantumError: invalid input, number of circuits ("; - msg << circuits.size() << ")" << "and number of probabilities ("; - msg << probs.size() << ") are not equal."; - throw std::invalid_argument(msg.str()); + throw std::invalid_argument( + "QuantumError: invalid input, number of circuits (" + + std::to_string(circuits.size()) + ") and number of probabilities (" + + std::to_string(probs.size()) + ") are not equal." + ); } // Check probability vector double total = 0.; @@ -259,6 +312,35 @@ void QuantumError::set_from_kraus(const std::vector &mats) { } +const cmatrix_t& QuantumError::superoperator() const { + // Check the superoperator is actually computed + // If not raise an exception + if (superoperator_.empty()) { + throw std::runtime_error("QuantumError: superoperator is empty."); + } + return superoperator_; +} + + +void QuantumError::compute_superoperator() { + // Initialize superoperator matrix to correct size + size_t dim = 1 << (2 * get_num_qubits()); + superoperator_.initialize(dim, dim); + // We use the superoperator simulator state to do this + QubitSuperoperator::State<> superop; + for (size_t j=0; j; + //----------------------------------------------------------------------- - // Error base required methods + // Error sampling //----------------------------------------------------------------------- // Sample a noisy implementation of op NoiseOps sample_noise(const reg_t &memory, - RngEngine &rng) const override; - - // Load a ReadoutError object from a JSON Error object - void load_from_json(const json_t &js) override; + RngEngine &rng) const; //----------------------------------------------------------------------- - // Additional class methods + // Initialization //----------------------------------------------------------------------- + // Load a ReadoutError object from a JSON Error object + void load_from_json(const json_t &js); + // Set the default assignment probabilities for measurement of each qubit // Each vector is a list of probabilities {P(0|q), P(1|q), ... P(n-1|q)} // For the no-error case this list viewed as a matrix should be an // identity matrix void set_probabilities(const std::vector &probs); + //----------------------------------------------------------------------- + // Utility + //----------------------------------------------------------------------- + + // Set number of qubits or memory bits for error + inline void set_num_qubits(uint_t num_qubits) {num_qubits_ = num_qubits;} + + // Get number of qubits or memory bits for error + inline uint_t get_num_qubits() const {return num_qubits_;} + protected: + + // Number of qubits/memory bits the error applies to + uint_t num_qubits_ = 0; + + // Vector of assignment probability vectors std::vector assignment_probabilities_; // threshold for checking probabilities diff --git a/src/simulators/densitymatrix/densitymatrix.hpp b/src/simulators/densitymatrix/densitymatrix.hpp index 402a921b4b..76216a34a1 100755 --- a/src/simulators/densitymatrix/densitymatrix.hpp +++ b/src/simulators/densitymatrix/densitymatrix.hpp @@ -61,6 +61,9 @@ class DensityMatrix : public UnitaryMatrix { // number of qubits an exception is raised. void initialize_from_vector(const cvector_t &data); + // Returns the number of qubits for the superoperator + virtual uint_t num_qubits() const override {return BaseMatrix::num_qubits_;} + //----------------------------------------------------------------------- // Apply Matrices //----------------------------------------------------------------------- @@ -118,7 +121,7 @@ class DensityMatrix : public UnitaryMatrix { // Convert qubit indicies to vectorized-density matrix qubitvector indices // For the QubitVector apply matrix function - reg_t superop_qubits(const reg_t &qubits) const; + virtual reg_t superop_qubits(const reg_t &qubits) const; // Construct a vectorized superoperator from a vectorized matrix // This is equivalent to vec(tensor(conj(A), A)) @@ -181,9 +184,9 @@ template reg_t DensityMatrix::superop_qubits(const reg_t &qubits) const { reg_t superop_qubits = qubits; // Number of qubits - const auto num_qubits = BaseMatrix::num_qubits(); + const auto nq = num_qubits(); for (const auto q: qubits) { - superop_qubits.push_back(q + num_qubits); + superop_qubits.push_back(q + nq); } return superop_qubits; } @@ -219,10 +222,10 @@ void DensityMatrix::apply_unitary_matrix(const reg_t &qubits, // Check if we apply as two N-qubit matrix multiplications or a single 2N-qubit matrix mult. if (qubits.size() > apply_unitary_threshold_) { // Apply as two N-qubit matrix mults - auto num_qubits = BaseMatrix::num_qubits(); + auto nq = num_qubits(); reg_t conj_qubits; for (const auto q: qubits) { - conj_qubits.push_back(q + num_qubits); + conj_qubits.push_back(q + nq); } // Apply id \otimes U BaseVector::apply_matrix(qubits, mat); @@ -250,7 +253,7 @@ void DensityMatrix::apply_cnot(const uint_t qctrl, const uint_t qtrgt) { std::vector> pairs = { {{1, 3}, {4, 12}, {5, 15}, {6, 14}, {7, 13}, {9, 11}} }; - const size_t nq = BaseMatrix::num_qubits(); + const size_t nq = num_qubits(); const reg_t qubits = {{qctrl, qtrgt, qctrl + nq, qtrgt + nq}}; BaseVector::apply_permutation_matrix(qubits, pairs); } @@ -266,7 +269,7 @@ void DensityMatrix::apply_cz(const uint_t q0, const uint_t q1) { BaseVector::data_[inds[13]] *= -1.; BaseVector::data_[inds[14]] *= -1.; }; - const auto nq = BaseMatrix::num_qubits(); + const auto nq = num_qubits(); const areg_t<4> qubits = {{q0, q1, q0 + nq, q1 + nq}}; BaseVector::apply_lambda(lambda, qubits); } @@ -276,7 +279,7 @@ void DensityMatrix::apply_swap(const uint_t q0, const uint_t q1) { std::vector> pairs = { {{1, 2}, {4, 8}, {5, 10}, {6, 9}, {7, 11}, {13, 14}} }; - const size_t nq = BaseMatrix::num_qubits(); + const size_t nq = num_qubits(); const reg_t qubits = {{q0, q1, q0 + nq, q1 + nq}}; BaseVector::apply_permutation_matrix(qubits, pairs); } @@ -289,7 +292,7 @@ void DensityMatrix::apply_x(const uint_t qubit) { std::swap(BaseVector::data_[inds[1]], BaseVector::data_[inds[2]]); }; // Use the lambda function - const areg_t<2> qubits = {{qubit, qubit + BaseMatrix::num_qubits()}}; + const areg_t<2> qubits = {{qubit, qubit + num_qubits()}}; BaseVector::apply_lambda(lambda, qubits); } @@ -303,7 +306,7 @@ void DensityMatrix::apply_y(const uint_t qubit) { BaseVector::data_[inds[2]] = cache; }; // Use the lambda function - const areg_t<2> qubits = {{qubit, qubit + BaseMatrix::num_qubits()}}; + const areg_t<2> qubits = {{qubit, qubit + num_qubits()}}; BaseVector::apply_lambda(lambda, qubits); } @@ -315,7 +318,7 @@ void DensityMatrix::apply_z(const uint_t qubit) { BaseVector::data_[inds[2]] *= -1; }; // Use the lambda function - const areg_t<2> qubits = {{qubit, qubit + BaseMatrix::num_qubits()}}; + const areg_t<2> qubits = {{qubit, qubit + num_qubits()}}; BaseVector::apply_lambda(lambda, qubits); } @@ -327,7 +330,7 @@ void DensityMatrix::apply_toffoli(const uint_t qctrl0, {{3, 7}, {11, 15}, {19, 23}, {24, 56}, {25, 57}, {26, 58}, {27, 63}, {28, 60}, {29, 61}, {30, 62}, {31, 59}, {35, 39}, {43,47}, {51, 55}} }; - const size_t nq = BaseMatrix::num_qubits(); + const size_t nq = num_qubits(); const reg_t qubits = {{qctrl0, qctrl1, qtrgt, qctrl0 + nq, qctrl1 + nq, qtrgt + nq}}; BaseVector::apply_permutation_matrix(qubits, pairs); diff --git a/src/simulators/qasm/qasm_controller.hpp b/src/simulators/qasm/qasm_controller.hpp index 4afa8ddd12..747071745c 100755 --- a/src/simulators/qasm/qasm_controller.hpp +++ b/src/simulators/qasm/qasm_controller.hpp @@ -23,6 +23,7 @@ #include "simulators/stabilizer/stabilizer_state.hpp" #include "simulators/tensor_network/tensor_network_state.hpp" #include "simulators/densitymatrix/densitymatrix_state.hpp" +#include "simulators/superoperator/superoperator_state.hpp" namespace AER { @@ -225,7 +226,6 @@ class QasmController : public Base::Controller { uint_t shots, State_t &state, const Initstate_t &initial_state, - const Method method, OutputData &data, RngEngine &rng) const; @@ -623,14 +623,27 @@ OutputData QasmController::run_circuit_helper(const Circuit &circ, // Note: this will set to `true` if sampling is enabled for the circuit data.add_additional_data("metadata", json_t::object({{"measure_sampling", false}})); - // Check if there is noise for the implementation + + // Choose execution method based on noise and method if (noise.is_ideal()) { run_circuit_without_noise(circ, shots, state, initial_state, method, data, rng); - } else if (!noise.has_quantum_errors()) { - run_circuit_without_noise(noise.sample_noise(circ, rng), - shots, state, initial_state, method, data, rng); + } + else if (method == Method::density_matrix && noise.has_quantum_errors()) { + // We can sample the noise model using superoperator method + // and then execute the resulting circuit containing superoperators + Noise::NoiseModel noise_cpy = noise; + noise_cpy.activate_superop_method(); + Circuit noise_circ = noise_cpy.sample_noise(circ, rng); + run_circuit_without_noise(noise_circ, shots, state, initial_state, method, data, rng); + } + else if (noise.has_quantum_errors() == false) { + // We can insert the readout errors from the noise model and then + // execute the resulting circuit + Circuit noise_circ = noise.sample_noise(circ, rng); + run_circuit_without_noise(noise_circ, shots, state, initial_state, method, data, rng); } else { - run_circuit_with_noise(circ, noise, shots, state, initial_state, method, data, rng); + // Run sampling a noisy instance of the circuit for each shot + run_circuit_with_noise(circ, noise, shots, state, initial_state, data, rng); } return data; } @@ -654,10 +667,8 @@ void QasmController::run_circuit_with_noise(const Circuit &circ, uint_t shots, State_t &state, const Initstate_t &initial_state, - const Method method, OutputData &data, RngEngine &rng) const { - // TODO: noise sampling for density matrix method // Sample a new noise circuit and optimize for each shot while(shots-- > 0) { Circuit noise_circ = noise.sample_noise(circ, rng); @@ -666,7 +677,7 @@ void QasmController::run_circuit_with_noise(const Circuit &circ, optimize_circuit(noise_circ, dummy, state, data); } run_single_shot(noise_circ, state, initial_state, data, rng); - } + } } diff --git a/src/simulators/superoperator/superoperator.hpp b/src/simulators/superoperator/superoperator.hpp new file mode 100755 index 0000000000..e742238f34 --- /dev/null +++ b/src/simulators/superoperator/superoperator.hpp @@ -0,0 +1,152 @@ +/** + * This code is part of Qiskit. + * + * (C) Copyright IBM 2018, 2019. + * + * 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. + */ + +#ifndef _qv_superoperator_hpp_ +#define _qv_superoperator_hpp_ + + +#include "framework/utils.hpp" +#include "simulators/densitymatrix/densitymatrix.hpp" + +namespace QV { + +//============================================================================ +// Superoperator class +//============================================================================ + +// This class is derived from the DensityMatrix class and stores an N-qubit +// superoperator as a 2 * N-qubit vector. +// The vector is formed using column-stacking vectorization of the +// superoperator (itself with respect to column-stacking vectorization). + +template +class Superoperator : public DensityMatrix { + +public: + // Parent class aliases + using BaseVector = QubitVector; + using BaseDensity = DensityMatrix; + using BaseUnitary = UnitaryMatrix; + + //----------------------------------------------------------------------- + // Constructors and Destructor + //----------------------------------------------------------------------- + + Superoperator() : Superoperator(0) {}; + explicit Superoperator(size_t num_qubits); + Superoperator(const Superoperator& obj) = delete; + Superoperator &operator=(const Superoperator& obj) = delete; + + //----------------------------------------------------------------------- + // Utility functions + //----------------------------------------------------------------------- + + // Set the size of the vector in terms of qubit number + virtual void set_num_qubits(size_t num_qubits) override; + + // Returns the number of qubits for the superoperator + virtual uint_t num_qubits() const override {return num_qubits_;} + + // Initialize to the identity superoperator + void initialize(); + + // Initializes the vector to a custom initial state. + // The matrix can either be superoperator matrix or unitary matrix. + // The type is inferred by the dimensions of the input matrix. + void initialize_from_matrix(const AER::cmatrix_t &data); + +protected: + // Number of qubits for the superoperator + size_t num_qubits_; +}; + +/******************************************************************************* + * + * Implementations + * + ******************************************************************************/ + + +//------------------------------------------------------------------------------ +// Constructors & Destructor +//------------------------------------------------------------------------------ + +template +Superoperator::Superoperator(size_t num_qubits) { + set_num_qubits(num_qubits); +} + +//------------------------------------------------------------------------------ +// Utility +//------------------------------------------------------------------------------ + +template +void Superoperator::set_num_qubits(size_t num_qubits) { + num_qubits_ = num_qubits; + // Superoperator is same size matrix as a unitary matrix + // of twice as many qubits + BaseDensity::set_num_qubits(2 * num_qubits); +} + + +template +void Superoperator::initialize() { + // Set underlying unitary matrix to identity + BaseUnitary::initialize(); +} + + +template +void Superoperator::initialize_from_matrix(const AER::cmatrix_t &mat) { + if (AER::Utils::is_square(mat)) { + const size_t nrows = mat.GetRows(); + if (nrows == BaseUnitary::rows_) { + // The matrix is the same size as the superoperator matrix so we + // initialze as the matrix. + BaseUnitary::initialize_from_matrix(mat); + return; + } else if (nrows * nrows == BaseUnitary::rows_) { + // If the input matrix has half the number of rows we assume it is + // A unitary matrix input so we convert to a superoperator + BaseUnitary::initialize_from_matrix( + AER::Utils::tensor_product(AER::Utils::conjugate(mat), mat) + ); + return; + } + } + // Throw an exception if the input matrix is the wrong size for + // unitary or superoperator input + throw std::runtime_error( + "Superoperator::initial matrix is wrong size (" + + std::to_string(BaseUnitary::rows_) + "," + + std::to_string(BaseUnitary::rows_) + ")!=(" + + std::to_string(mat.GetRows()) + "," + std::to_string(mat.GetColumns()) + ")." + ); +}; + + +//------------------------------------------------------------------------------ +} // end namespace QV +//------------------------------------------------------------------------------ + +// ostream overload for templated qubitvector +template +inline std::ostream &operator<<(std::ostream &out, const QV::Superoperator&m) { + out << m.matrix(); + return out; +} + +//------------------------------------------------------------------------------ +#endif // end module + diff --git a/src/simulators/superoperator/superoperator_state.hpp b/src/simulators/superoperator/superoperator_state.hpp new file mode 100755 index 0000000000..b8b192f1eb --- /dev/null +++ b/src/simulators/superoperator/superoperator_state.hpp @@ -0,0 +1,477 @@ +/** + * This code is part of Qiskit. + * + * (C) Copyright IBM 2018, 2019. + * + * 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. + */ + +#ifndef _superoperator_state_hpp +#define _superoperator_state_hpp + +#include +#define _USE_MATH_DEFINES +#include + +#include "framework/utils.hpp" +#include "framework/json.hpp" +#include "base/state.hpp" +#include "superoperator.hpp" + + +namespace AER { +namespace QubitSuperoperator { + +// Allowed gates enum class +enum class Gates { + u1, u2, u3, id, x, y, z, h, s, sdg, t, tdg, // single qubit + cx, cz, swap, // two qubit + ccx // three qubit +}; + +// Allowed snapshots enum class +enum class Snapshots {superop}; + +//========================================================================= +// QubitUnitary State subclass +//========================================================================= + +template > +class State : public Base::State { +public: + using BaseState = Base::State; + + State() = default; + virtual ~State() = default; + + //----------------------------------------------------------------------- + // Base class overrides + //----------------------------------------------------------------------- + + // Return the string name of the State class + virtual std::string name() const override {return "superoperator";} + + // Return the set of qobj instruction types supported by the State + virtual Operations::OpSet::optypeset_t allowed_ops() const override { + return Operations::OpSet::optypeset_t({ + Operations::OpType::gate, + Operations::OpType::reset, + Operations::OpType::snapshot, + Operations::OpType::barrier, + Operations::OpType::matrix, + Operations::OpType::kraus, + Operations::OpType::superop + }); + } + + // Return the set of qobj gate instruction names supported by the State + virtual stringset_t allowed_gates() const override { + return {"U", "CX", "u1", "u2", "u3", "cx", "cz", "swap", + "id", "x", "y", "z", "h", "s", "sdg", "t", "tdg", "ccx"}; + } + + // Return the set of qobj snapshot types supported by the State + virtual stringset_t allowed_snapshots() const override { + return {"superoperator"}; + } + + // Apply a sequence of operations by looping over list + // If the input is not in allowed_ops an exeption will be raised. + virtual void apply_ops(const std::vector &ops, + OutputData &data, + RngEngine &rng) override; + + // Initializes an n-qubit unitary to the identity matrix + virtual void initialize_qreg(uint_t num_qubits) override; + + // Initializes to a specific n-qubit unitary superop + virtual void initialize_qreg(uint_t num_qubits, + const data_t &unitary) override; + + // Returns the required memory for storing an n-qubit state in megabytes. + // For this state the memory is indepdentent of the number of ops + // and is approximately 16 * 1 << 4 * num_qubits bytes + virtual size_t required_memory_mb(uint_t num_qubits, + const std::vector &ops) override; + + // Load the threshold for applying OpenMP parallelization + // if the controller/engine allows threads for it + // Config: {"omp_qubit_threshold": 3} + virtual void set_config(const json_t &config) override; + + //----------------------------------------------------------------------- + // Additional methods + //----------------------------------------------------------------------- + + // Initializes to a specific n-qubit unitary given as a complex matrix + virtual void initialize_qreg(uint_t num_qubits, const cmatrix_t &unitary); + + // Initialize OpenMP settings for the underlying QubitVector class + void initialize_omp(); + +protected: + + //----------------------------------------------------------------------- + // Apply Instructions + //----------------------------------------------------------------------- + + // Applies a Gate operation to the state class. + // This should support all and only the operations defined in + // allowed_operations. + void apply_gate(const Operations::Op &op); + + // Apply a supported snapshot instruction + // If the input is not in allowed_snapshots an exeption will be raised. + virtual void apply_snapshot(const Operations::Op &op, OutputData &data); + + // Apply a matrix to given qubits (identity on all other qubits) + void apply_matrix(const reg_t &qubits, const cmatrix_t & mat); + + // Apply a matrix to given qubits (identity on all other qubits) + void apply_matrix(const reg_t &qubits, const cvector_t & vmat); + + // Reset the specified qubits to the |0> state by simulating + // a measurement, applying a conditional x-gate if the outcome is 1, and + // then discarding the outcome. + void apply_reset(const reg_t &qubits); + + // Apply a Kraus error operation + void apply_kraus(const reg_t &qubits, + const std::vector &krausops); + + //----------------------------------------------------------------------- + // 1-Qubit Gates + //----------------------------------------------------------------------- + + // Optimize phase gate with diagonal [1, phase] + void apply_gate_phase(const uint_t qubit, const complex_t phase); + + //----------------------------------------------------------------------- + // Multi-controlled u3 + //----------------------------------------------------------------------- + + // Apply N-qubit multi-controlled single qubit waltz gate specified by + // parameters u3(theta, phi, lambda) + // NOTE: if N=1 this is just a regular u3 gate. + void apply_gate_u3(const uint_t qubit, + const double theta, + const double phi, + const double lambda); + + //----------------------------------------------------------------------- + // Config Settings + //----------------------------------------------------------------------- + + // OpenMP qubit threshold + int omp_qubit_threshold_ = 3; + + // Threshold for chopping small values to zero in JSON + double json_chop_threshold_ = 1e-10; + + // Table of allowed gate names to gate enum class members + const static stringmap_t gateset_; +}; + + +//============================================================================ +// Implementation: Allowed ops and gateset +//============================================================================ + +template +const stringmap_t State::gateset_({ + // Single qubit gates + {"id", Gates::id}, // Pauli-Identity gate + {"x", Gates::x}, // Pauli-X gate + {"y", Gates::y}, // Pauli-Y gate + {"z", Gates::z}, // Pauli-Z gate + {"s", Gates::s}, // Phase gate (aka sqrt(Z) gate) + {"sdg", Gates::sdg}, // Conjugate-transpose of Phase gate + {"h", Gates::h}, // Hadamard gate (X + Z / sqrt(2)) + {"t", Gates::t}, // T-gate (sqrt(S)) + {"tdg", Gates::tdg}, // Conjguate-transpose of T gate + // Waltz Gates + {"u1", Gates::u1}, // zero-X90 pulse waltz gate + {"u2", Gates::u2}, // single-X90 pulse waltz gate + {"u3", Gates::u3}, // two X90 pulse waltz gate + {"U", Gates::u3}, // two X90 pulse waltz gate + // Two-qubit gates + {"CX", Gates::cx}, // Controlled-X gate (CNOT) + {"cx", Gates::cx}, // Controlled-X gate (CNOT) + {"cz", Gates::cz}, // Controlled-Z gate + {"swap", Gates::swap}, // SWAP gate + // Three-qubit gates + {"ccx", Gates::ccx} // Controlled-CX gate (Toffoli) +}); + +//============================================================================ +// Implementation: Base class method overrides +//============================================================================ + +template +void State::apply_ops(const std::vector &ops, + OutputData &data, + RngEngine &rng) { + // Simple loop over vector of input operations + for (const auto op: ops) { + switch (op.type) { + case Operations::OpType::barrier: + break; + case Operations::OpType::gate: + // Note conditionals will always fail since no classical registers + if (BaseState::creg_.check_conditional(op)) + apply_gate(op); + break; + case Operations::OpType::reset: + apply_reset(op.qubits); + break; + case Operations::OpType::matrix: + apply_matrix(op.qubits, op.mats[0]); + break; + case Operations::OpType::kraus: + apply_kraus(op.qubits, op.mats); + break; + case Operations::OpType::superop: + BaseState::qreg_.apply_superop_matrix(op.qubits, Utils::vectorize_matrix(op.mats[0])); + break; + case Operations::OpType::snapshot: + apply_snapshot(op, data); + break; + default: + throw std::invalid_argument("QubitSuperoperator::State::invalid instruction \'" + + op.name + "\'."); + } + } +} + +template +size_t State::required_memory_mb(uint_t num_qubits, + const std::vector &ops) { + // An n-qubit unitary as 2^4n complex doubles + // where each complex double is 16 bytes + (void)ops; // avoid unused variable compiler warning + size_t shift_mb = std::max(0, num_qubits + 4 - 20); + size_t mem_mb = 1ULL << (4 * shift_mb); + return mem_mb; +} + + +template +void State::set_config(const json_t &config) { + // Set OMP threshold for state update functions + JSON::get_value(omp_qubit_threshold_, "superoperator_parallel_threshold", config); + + // Set threshold for truncating snapshots + JSON::get_value(json_chop_threshold_, "zero_threshold", config); + BaseState::qreg_.set_json_chop_threshold(json_chop_threshold_); +} + + +template +void State::initialize_qreg(uint_t num_qubits) { + initialize_omp(); + BaseState::qreg_.set_num_qubits(num_qubits); + BaseState::qreg_.initialize(); +} + + +template +void State::initialize_qreg(uint_t num_qubits, + const data_t &supermat) { + // Check dimension of state + if (supermat.num_qubits() != num_qubits) { + throw std::invalid_argument("QubitSuperoperator::State::initialize: initial state does not match qubit number"); + } + initialize_omp(); + BaseState::qreg_.set_num_qubits(num_qubits); + const size_t sz = 1ULL << BaseState::qreg_.size(); + BaseState::qreg_.initialize_from_data(supermat.data(), sz); +} + + +template +void State::initialize_qreg(uint_t num_qubits, + const cmatrix_t &mat) { + + // Check dimension of unitary + const auto sz_uni = 1ULL << (2 * num_qubits); + const auto sz_super = 1ULL << (4 * num_qubits); + if (mat.size() != sz_uni && mat.size() != sz_super) { + throw std::invalid_argument( + "QubitSuperoperator::State::initialize: initial state does not match qubit number"); + } + initialize_omp(); + BaseState::qreg_.set_num_qubits(num_qubits); + BaseState::qreg_.initialize_from_matrix(mat); +} + + +template +void State::initialize_omp() { + BaseState::qreg_.set_omp_threshold(omp_qubit_threshold_); + if (BaseState::threads_ > 0) + BaseState::qreg_.set_omp_threads(BaseState::threads_); // set allowed OMP threads in qubitvector +} + +//========================================================================= +// Implementation: Reset +//========================================================================= + +template +void State::apply_reset(const reg_t &qubits) { + // TODO: This can be more efficient by adding reset + // to base class rather than doing a matrix multiplication + // where all but 1 row is zeros. + const auto reset_op = Utils::SMatrix::reset(1ULL << qubits.size()); + BaseState::qreg_.apply_superop_matrix(qubits, Utils::vectorize_matrix(reset_op)); +} + +//========================================================================= +// Implementation: Kraus Noise +//========================================================================= + +template +void State::apply_kraus(const reg_t &qubits, + const std::vector &kmats) { + // Convert to Superoperator + const auto nrows = kmats[0].GetRows(); + cmatrix_t superop(nrows * nrows, nrows * nrows); + for (const auto kraus : kmats) { + superop += Utils::tensor_product(Utils::conjugate(kraus), kraus); + } + BaseState::qreg_.apply_superop_matrix(qubits, Utils::vectorize_matrix(superop)); +} + +//========================================================================= +// Implementation: Gates +//========================================================================= + +template +void State::apply_gate(const Operations::Op &op) { + // Look for gate name in gateset + auto it = gateset_.find(op.name); + if (it == gateset_.end()) + throw std::invalid_argument("Unitary::State::invalid gate instruction \'" + + op.name + "\'."); + switch (it -> second) { + case Gates::u3: + apply_gate_u3(op.qubits[0], + std::real(op.params[0]), + std::real(op.params[1]), + std::real(op.params[2])); + break; + case Gates::u2: + apply_gate_u3(op.qubits[0], + M_PI / 2., + std::real(op.params[0]), + std::real(op.params[1])); + break; + case Gates::u1: + apply_gate_phase(op.qubits[0], std::exp(complex_t(0., 1.) * op.params[0])); + break; + case Gates::cx: + BaseState::qreg_.apply_cnot(op.qubits[0], op.qubits[1]); + break; + case Gates::cz: + BaseState::qreg_.apply_cz(op.qubits[0], op.qubits[1]); + break; + case Gates::id: + break; + case Gates::x: + BaseState::qreg_.apply_x(op.qubits[0]); + break; + case Gates::y: + BaseState::qreg_.apply_y(op.qubits[0]); + break; + case Gates::z: + BaseState::qreg_.apply_z(op.qubits[0]); + break; + case Gates::h: + apply_gate_u3(op.qubits[0], M_PI / 2., 0., M_PI); + break; + case Gates::s: + apply_gate_phase(op.qubits[0], complex_t(0., 1.)); + break; + case Gates::sdg: + apply_gate_phase(op.qubits[0], complex_t(0., -1.)); + break; + case Gates::t: { + const double isqrt2{1. / std::sqrt(2)}; + apply_gate_phase(op.qubits[0], complex_t(isqrt2, isqrt2)); + } break; + case Gates::tdg: { + const double isqrt2{1. / std::sqrt(2)}; + apply_gate_phase(op.qubits[0], complex_t(isqrt2, -isqrt2)); + } break; + case Gates::swap: { + BaseState::qreg_.apply_swap(op.qubits[0], op.qubits[1]); + } break; + case Gates::ccx: + BaseState::qreg_.apply_toffoli(op.qubits[0], op.qubits[1], op.qubits[2]); + break; + default: + // We shouldn't reach here unless there is a bug in gateset + throw std::invalid_argument("Superoperator::State::invalid gate instruction \'" + + op.name + "\'."); + } +} + +template +void State::apply_matrix(const reg_t &qubits, const cmatrix_t &mat) { + if (qubits.empty() == false && mat.size() > 0) { + BaseState::qreg_.apply_unitary_matrix(qubits, Utils::vectorize_matrix(mat)); + } +} + +template +void State::apply_matrix(const reg_t &qubits, const cvector_t &vmat) { + // Check if diagonal matrix + if (vmat.size() == 1ULL << qubits.size()) { + BaseState::qreg_.apply_diagonal_unitary_matrix(qubits, vmat); + } else { + BaseState::qreg_.apply_unitary_matrix(qubits, vmat); + } +} + + +template +void State::apply_gate_phase(uint_t qubit, complex_t phase) { + cvector_t diag(2); + diag[0] = 1.0; + diag[1] = phase; + BaseState::qreg_.apply_diagonal_unitary_matrix(reg_t({qubit}), diag); +} + + +template +void State::apply_gate_u3(const uint_t qubit, + double theta, + double phi, + double lambda) { + const auto u3 = Utils::VMatrix::u3(theta, phi, lambda); + BaseState::qreg_.apply_unitary_matrix(reg_t({qubit}), u3); +} + + +template +void State::apply_snapshot(const Operations::Op &op, + OutputData &data) { + // Look for snapshot type in snapshotset + if (op.name == "superopertor" || op.name == "state") { + BaseState::snapshot_state(op, data, "superoperator"); + } else { + throw std::invalid_argument("QubitSuperoperator::State::invalid snapshot instruction \'" + + op.name + "\'."); + } +} + +//------------------------------------------------------------------------------ +} // end namespace QubitSuperoperator +} // end namespace AER +//------------------------------------------------------------------------------ +#endif diff --git a/test/terra/backends/qasm_simulator/qasm_measure.py b/test/terra/backends/qasm_simulator/qasm_measure.py index 1a34946ebc..1933bb8912 100644 --- a/test/terra/backends/qasm_simulator/qasm_measure.py +++ b/test/terra/backends/qasm_simulator/qasm_measure.py @@ -19,6 +19,7 @@ from qiskit.providers.aer.noise import NoiseModel from qiskit.providers.aer.noise.errors import ReadoutError, depolarizing_error + class QasmMeasureTests: """QasmSimulator measure tests.""" @@ -92,6 +93,60 @@ def test_measure_nondeterministic_without_sampling(self): self.assertIn("measure_sampling", res.metadata) self.assertEqual(res.metadata["measure_sampling"], False) + def test_measure_sampling_with_readouterror(self): + """Test QasmSimulator measure with deterministic counts with sampling and readout-error""" + readout_error = [0.01, 0.1] + noise_model = NoiseModel() + readout = [[1.0 - readout_error[0], readout_error[0]], + [readout_error[1], 1.0 - readout_error[1]]] + noise_model.add_all_qubit_readout_error(ReadoutError(readout)) + + shots = 1000 + circuits = ref_measure.measure_circuits_deterministic( + allow_sampling=True) + targets = ref_measure.measure_counts_deterministic(shots) + qobj = assemble(circuits, self.SIMULATOR, shots=shots) + result = self.SIMULATOR.run( + qobj, noise_model=noise_model, + backend_options=self.BACKEND_OPTS).result() + self.is_completed(result) + + # Test sampling was enable + for res in result.results: + self.assertIn("measure_sampling", res.metadata) + self.assertEqual(res.metadata["measure_sampling"], True) + + def test_measure_sampling_with_quantum_noise(self): + """Test QasmSimulator measure with deterministic counts with sampling and readout-error""" + readout_error = [0.01, 0.1] + noise_model = NoiseModel() + depolarizing = {'u3': (1, 0.001), 'cx': (2, 0.02)} + readout = [[1.0 - readout_error[0], readout_error[0]], + [readout_error[1], 1.0 - readout_error[1]]] + noise_model.add_all_qubit_readout_error(ReadoutError(readout)) + for gate, (num_qubits, gate_error) in depolarizing.items(): + noise_model.add_all_qubit_quantum_error( + depolarizing_error(gate_error, num_qubits), gate) + + shots = 1000 + circuits = ref_measure.measure_circuits_deterministic( + allow_sampling=True) + targets = ref_measure.measure_counts_deterministic(shots) + qobj = assemble(circuits, self.SIMULATOR, shots=shots) + result = self.SIMULATOR.run( + qobj, noise_model=noise_model, + backend_options=self.BACKEND_OPTS).result() + self.is_completed(result) + + method = self.BACKEND_OPTS.get("method") + for res in result.results: + self.assertIn("measure_sampling", res.metadata) + if method == "density_matrix": + self.assertEqual(res.metadata["measure_sampling"], True) + else: + self.assertEqual(res.metadata["measure_sampling"], False) + + class QasmMultiQubitMeasureTests: """QasmSimulator measure tests.""" @@ -160,64 +215,7 @@ def test_measure_nondeterministic_multi_qubit_without_sampling(self): qobj, backend_options=self.BACKEND_OPTS).result() self.is_completed(result) self.compare_counts(result, circuits, targets, delta=0.05 * shots) - - # Test sampling was disabled - for res in result.results: - self.assertIn("measure_sampling", res.metadata) - self.assertEqual(res.metadata["measure_sampling"], False) - def test_measure_sampling_with_readouterror(self): - """Test QasmSimulator measure with deterministic counts with sampling and readout-error""" - readout_error = [0.01, 0.1] - noise_model = NoiseModel() - readout = [[1.0 - readout_error[0], readout_error[0]], - [readout_error[1], 1.0 - readout_error[1]]] - noise_model.add_all_qubit_readout_error(ReadoutError(readout)) - - shots = 1000 - circuits = ref_measure.measure_circuits_deterministic( - allow_sampling=True) - targets = ref_measure.measure_counts_deterministic(shots) - qobj = assemble(circuits, self.SIMULATOR, shots=shots) - result = self.SIMULATOR.run( - qobj, - noise_model=noise_model, - backend_options=self.BACKEND_OPTS).result() - self.is_completed(result) - - # Test sampling was disabled - for res in result.results: - self.assertIn("measure_sampling", res.metadata) - self.assertEqual(res.metadata["measure_sampling"], True) - - def test_measure_sampling_with_quantum_noise(self): - """Test QasmSimulator measure with deterministic counts with sampling and readout-error""" - readout_error = [0.01, 0.1] - noise_model = NoiseModel() - depolarizing = {'u3': (1, 0.001), 'cx': (2, 0.02)} - readout = [[1.0 - readout_error[0], readout_error[0]], - [readout_error[1], 1.0 - readout_error[1]]] - noise_model.add_all_qubit_readout_error(ReadoutError(readout)) - for gate, (num_qubits, gate_error) in depolarizing.items(): - noise_model.add_all_qubit_quantum_error( - depolarizing_error(gate_error, num_qubits), gate) - - shots = 1000 - circuits = ref_measure.measure_circuits_deterministic( - allow_sampling=True) - targets = ref_measure.measure_counts_deterministic(shots) - qobj = assemble(circuits, self.SIMULATOR, shots=shots) - result = self.SIMULATOR.run( - qobj, - noise_model=noise_model, - backend_options=self.BACKEND_OPTS).result() - self.is_completed(result) - - for res in result.results: - self.assertIn("measure_sampling", res.metadata) - self.assertEqual(res.metadata["measure_sampling"], False) - - # Test sampling was disabled for res in result.results: self.assertIn("measure_sampling", res.metadata) self.assertEqual(res.metadata["measure_sampling"], False)