From 64311f45318d93ba596c5acac0b318bdf6e5d1de Mon Sep 17 00:00:00 2001 From: Hiroshi Horii Date: Fri, 3 Mar 2023 00:46:53 +0900 Subject: [PATCH] apply clang-format only src/simulators/*stabilizer --- .../extended_stabilizer/ch_runner.hpp | 684 ++++---- .../chlib/chstabilizer.hpp | 1449 ++++++++--------- .../extended_stabilizer/chlib/core.hpp | 490 +++--- .../extended_stabilizer_state.hpp | 1017 ++++++------ src/simulators/extended_stabilizer/gates.hpp | 265 ++- src/simulators/stabilizer/binary_vector.hpp | 84 +- src/simulators/stabilizer/clifford.hpp | 546 ++++--- src/simulators/stabilizer/pauli.hpp | 62 +- .../stabilizer/stabilizer_state.hpp | 494 +++--- 9 files changed, 2284 insertions(+), 2807 deletions(-) diff --git a/src/simulators/extended_stabilizer/ch_runner.hpp b/src/simulators/extended_stabilizer/ch_runner.hpp index b9303e727e..0029fadfdb 100644 --- a/src/simulators/extended_stabilizer/ch_runner.hpp +++ b/src/simulators/extended_stabilizer/ch_runner.hpp @@ -15,8 +15,8 @@ #ifndef chsimulator_runner_hpp #define chsimulator_runner_hpp -#include "chlib/core.hpp" #include "chlib/chstabilizer.hpp" +#include "chlib/core.hpp" #include "gates.hpp" #include "framework/json.hpp" @@ -24,14 +24,13 @@ #include "framework/rng.hpp" #include "framework/types.hpp" - #define _USE_MATH_DEFINES #include #include +#include #include #include -#include #include #ifdef _OPENMP @@ -42,8 +41,8 @@ namespace CHSimulator { using chstabilizer_t = StabilizerState; -const double T_ANGLE = M_PI/4.; -const double TDG_ANGLE = -1.*M_PI/4.; +const double T_ANGLE = M_PI / 4.; +const double TDG_ANGLE = -1. * M_PI / 4.; const U1Sample t_sample(T_ANGLE); const U1Sample tdg_sample(TDG_ANGLE); @@ -52,8 +51,7 @@ const uint_t ZERO = 0ULL; thread_local std::unordered_map Z_ROTATIONS; -class Runner -{ +class Runner { private: uint_t n_qubits_ = 0; uint_t num_states_ = 0; @@ -74,34 +72,31 @@ class Runner json_t serialize_state(uint_t rank) const; public: - Runner(): n_qubits_(0), num_states_(0) {}; + Runner() : n_qubits_(0), num_states_(0){}; Runner(uint_t n_qubits); virtual ~Runner() = default; void initialize(uint_t n_qubits); void initialize_omp(uint_t n_threads, uint_t threshold_rank); - bool empty() const - { - return (n_qubits_ == 0 || num_states_ == 0); - } + bool empty() const { return (n_qubits_ == 0 || num_states_ == 0); } uint_t get_num_states() const; uint_t get_n_qubits() const; bool check_omp_threshold(); - //Convert each state to a json object and return it. + // Convert each state to a json object and return it. std::vector serialize_decomposition() const; - //Creates num_states_ copies of the 'base state' in the runner - //This will be either the |0>^n state, or a stabilizer state - //produced by applying the first m Clifford gates of the - //circuit. + // Creates num_states_ copies of the 'base state' in the runner + // This will be either the |0>^n state, or a stabilizer state + // produced by applying the first m Clifford gates of the + // circuit. void initialize_decomposition(uint_t n_states, double delta); - //Check if the coefficient omega is 0 + // Check if the coefficient omega is 0 bool check_eps(uint_t rank); - //Methods for applying gates + // Methods for applying gates void apply_cx(uint_t control, uint_t target, uint_t rank); void apply_cz(uint_t control, uint_t target, uint_t rank); void apply_swap(uint_t qubit_1, uint_t qubit_2, uint_t rank); @@ -113,422 +108,363 @@ class Runner void apply_x(uint_t qubit, uint_t rank); void apply_y(uint_t qubit, uint_t rank); void apply_z(uint_t qubit, uint_t rank); - //Methods for non-clifford gates, including parameters for the Sampling + // Methods for non-clifford gates, including parameters for the Sampling void apply_t(uint_t qubit, double r, int rank); void apply_tdag(uint_t qubit, double r, int rank); void apply_u1(uint_t qubit, complex_t lambda, double r, int rank); - void apply_ccx(uint_t control_1, uint_t control_2, uint_t target, uint_t branch, int rank); - void apply_ccz(uint_t control_1, uint_t control_2, uint_t target, uint_t branch, int rank); - //Measure a Pauli projector on each term in the decomposition and update their coefficients - // omega. + void apply_ccx(uint_t control_1, uint_t control_2, uint_t target, + uint_t branch, int rank); + void apply_ccz(uint_t control_1, uint_t control_2, uint_t target, + uint_t branch, int rank); + // Measure a Pauli projector on each term in the decomposition and update + // their coefficients + // omega. void apply_pauli(pauli_t &P); void apply_pauli_projector(const std::vector &generators); - void apply_pauli_projector(const std::vector &generators, uint_t rank); - //Routine for Norm Estimation, thin wrapper for the CHSimulator method that uses AER::RngEngine - //to set up the estimation routine. - double norm_estimation(uint_t n_samples, uint_t n_repetitions, AER::RngEngine &rng); - double norm_estimation(uint_t n_samples, uint_t n_repetitions, std::vector generators, AER::RngEngine &rng); - - //Metropolis Estimation for sampling from the output distribution + void apply_pauli_projector(const std::vector &generators, + uint_t rank); + // Routine for Norm Estimation, thin wrapper for the CHSimulator method that + // uses AER::RngEngine to set up the estimation routine. + double norm_estimation(uint_t n_samples, uint_t n_repetitions, + AER::RngEngine &rng); + double norm_estimation(uint_t n_samples, uint_t n_repetitions, + std::vector generators, AER::RngEngine &rng); + + // Metropolis Estimation for sampling from the output distribution uint_t metropolis_estimation(uint_t n_steps, AER::RngEngine &rng); - std::vector metropolis_estimation(uint_t n_steps, uint_t n_shots, AER::RngEngine &rng); + std::vector metropolis_estimation(uint_t n_steps, uint_t n_shots, + AER::RngEngine &rng); // NE Sampling Routines - uint_t ne_single_sample(uint_t default_samples, uint_t repetitions, bool preserve_states, - const AER::reg_t &qubits, AER::RngEngine &rng); - - std::vector ne_probabilities(uint_t default_samples, uint_t repetitions, - const AER::reg_t &qubits, AER::RngEngine &rng); - //Efficient Sampler for the output distribution of a stabilizer state + uint_t ne_single_sample(uint_t default_samples, uint_t repetitions, + bool preserve_states, const AER::reg_t &qubits, + AER::RngEngine &rng); + + std::vector ne_probabilities(uint_t default_samples, + uint_t repetitions, + const AER::reg_t &qubits, + AER::RngEngine &rng); + // Efficient Sampler for the output distribution of a stabilizer state uint_t stabilizer_sampler(AER::RngEngine &rng); std::vector stabilizer_sampler(uint_t n_shots, AER::RngEngine &rng); - //Utilities for the state-vector snapshot. + // Utilities for the state-vector snapshot. complex_t amplitude(uint_t x_measure); AER::Vector statevector(); - }; //========================================================================= // Implementation //========================================================================= -Runner::Runner(uint_t num_qubits) -{ - initialize(num_qubits); -} +Runner::Runner(uint_t num_qubits) { initialize(num_qubits); } -void Runner::initialize(uint_t num_qubits) -{ +void Runner::initialize(uint_t num_qubits) { states_.clear(); coefficients_.clear(); n_qubits_ = num_qubits; num_states_ = 1; num_threads_ = 1; states_ = std::vector(1, chstabilizer_t(num_qubits)); - coefficients_.push_back(complex_t(1.,0.)); + coefficients_.push_back(complex_t(1., 0.)); } -void Runner::initialize_decomposition(uint_t n_states, double delta) -{ +void Runner::initialize_decomposition(uint_t n_states, double delta) { num_states_ = n_states; states_.reserve(num_states_); coefficients_.reserve(num_states_); - if(states_.size() > 1 || coefficients_.size() > 1) - { - throw std::runtime_error(std::string("CHSimulator::Runner was initialized without") + - std::string("being properly cleared since the last ") + - std::string("experiment.")); + if (states_.size() > 1 || coefficients_.size() > 1) { + throw std::runtime_error( + std::string("CHSimulator::Runner was initialized without") + + std::string("being properly cleared since the last ") + + std::string("experiment.")); } - coefficients_[0] = complex_t(1./ delta,0.); + coefficients_[0] = complex_t(1. / delta, 0.); chstabilizer_t base_sate(states_[0]); complex_t coeff(coefficients_[0]); - for(uint_t i=1; i omp_threshold_; -} +bool Runner::check_omp_threshold() { return num_states_ > omp_threshold_; } //------------------------------------------------------------------------- // Operations on the decomposition //------------------------------------------------------------------------- -void Runner::apply_pauli(pauli_t &P) -{ +void Runner::apply_pauli(pauli_t &P) { const int_t END = num_states_; - #pragma omp parallel for if (num_states_ > omp_threshold_ && num_threads_ > 1) num_threads(num_threads_) - for(int_t i=0; i omp_threshold_ && num_threads_ > 1) \ + num_threads(num_threads_) + for (int_t i = 0; i < END; i++) states_[i].MeasurePauli(P); - } } -void Runner::apply_pauli_projector(const std::vector &generators) -{ +void Runner::apply_pauli_projector(const std::vector &generators) { const int_t END = num_states_; - #pragma omp parallel for if(num_states_ > omp_threshold_ && num_threads_ > 1) num_threads(num_threads_) - for(int_t i=0; i omp_threshold_ && num_threads_ > 1) \ + num_threads(num_threads_) + for (int_t i = 0; i < END; i++) apply_pauli_projector(generators, i); - } } -void Runner::apply_pauli_projector(const std::vector &generators, uint_t rank) -{ +void Runner::apply_pauli_projector(const std::vector &generators, + uint_t rank) { states_[rank].MeasurePauliProjector(generators); } -bool Runner::check_eps(uint_t rank) -{ - return (states_[rank].Omega().eps == 1); -} +bool Runner::check_eps(uint_t rank) { return (states_[rank].Omega().eps == 1); } -void Runner::apply_cx(uint_t control, uint_t target, uint_t rank) -{ +void Runner::apply_cx(uint_t control, uint_t target, uint_t rank) { states_[rank].CX(control, target); } -void Runner::apply_cz(uint_t control, uint_t target, uint_t rank) -{ +void Runner::apply_cz(uint_t control, uint_t target, uint_t rank) { states_[rank].CZ(control, target); } -void Runner::apply_swap(uint_t qubit_1, uint_t qubit_2, uint_t rank) -{ +void Runner::apply_swap(uint_t qubit_1, uint_t qubit_2, uint_t rank) { states_[rank].CX(qubit_1, qubit_2); states_[rank].CX(qubit_2, qubit_1); states_[rank].CX(qubit_1, qubit_2); } +void Runner::apply_h(uint_t qubit, uint_t rank) { states_[rank].H(qubit); } -void Runner::apply_h(uint_t qubit, uint_t rank) -{ - states_[rank].H(qubit); -} +void Runner::apply_s(uint_t qubit, uint_t rank) { states_[rank].S(qubit); } -void Runner::apply_s(uint_t qubit, uint_t rank) -{ - states_[rank].S(qubit); -} - -void Runner::apply_sdag(uint_t qubit, uint_t rank) -{ +void Runner::apply_sdag(uint_t qubit, uint_t rank) { states_[rank].Sdag(qubit); } -void Runner::apply_sx(uint_t qubit, uint_t rank) -{ +void Runner::apply_sx(uint_t qubit, uint_t rank) { states_[rank].Sdag(qubit); states_[rank].H(qubit); states_[rank].Sdag(qubit); } -void Runner::apply_sxdg(uint_t qubit, uint_t rank) -{ +void Runner::apply_sxdg(uint_t qubit, uint_t rank) { states_[rank].S(qubit); states_[rank].H(qubit); states_[rank].S(qubit); } -void Runner::apply_x(uint_t qubit, uint_t rank) -{ - states_[rank].X(qubit); -} +void Runner::apply_x(uint_t qubit, uint_t rank) { states_[rank].X(qubit); } -void Runner::apply_y(uint_t qubit, uint_t rank) -{ - states_[rank].Y(qubit); -} +void Runner::apply_y(uint_t qubit, uint_t rank) { states_[rank].Y(qubit); } -void Runner::apply_z(uint_t qubit, uint_t rank) -{ - states_[rank].Z(qubit); -} +void Runner::apply_z(uint_t qubit, uint_t rank) { states_[rank].Z(qubit); } -void Runner::apply_t(uint_t qubit, double r, int rank) -{ +void Runner::apply_t(uint_t qubit, double r, int rank) { sample_branch_t branch = t_sample.sample(r); coefficients_[rank] *= branch.first; if (branch.second == Gates::s) - { states_[rank].S(qubit); - } } -void Runner::apply_tdag(uint_t qubit, double r, int rank) -{ +void Runner::apply_tdag(uint_t qubit, double r, int rank) { sample_branch_t branch = tdg_sample.sample(r); coefficients_[rank] *= branch.first; if (branch.second == Gates::sdg) - { states_[rank].Sdag(qubit); - } } -void Runner::apply_u1(uint_t qubit, complex_t param, double r, int rank) -{ +void Runner::apply_u1(uint_t qubit, complex_t param, double r, int rank) { double lambda = std::real(param); - auto it = Z_ROTATIONS.find(lambda); //Look for cached z_rotations + auto it = Z_ROTATIONS.find(lambda); // Look for cached z_rotations sample_branch_t branch; - if (it == Z_ROTATIONS.end()) - { + if (it == Z_ROTATIONS.end()) { U1Sample rotation(lambda); Z_ROTATIONS.insert({lambda, rotation}); branch = rotation.sample(r); - } - else - { + } else { branch = it->second.sample(r); } coefficients_[rank] *= branch.first; - switch(branch.second) - { - case Gates::s: - states_[rank].S(qubit); - break; - case Gates::sdg: - states_[rank].Sdag(qubit); - break; - case Gates::z: - states_[rank].Z(qubit); - break; - default: - break; - } -} - -void Runner::apply_ccx(uint_t control_1, uint_t control_2, uint_t target, uint_t branch, int rank) -{ - switch(branch) //Decomposition of the CCX gate into Cliffords - { - case 1: - states_[rank].CZ(control_1, control_2); - break; - case 2: - states_[rank].CX(control_1, target); - break; - case 3: - states_[rank].CX(control_2, target); - break; - case 4: - states_[rank].CZ(control_1, control_2); - states_[rank].CX(control_1, target); - states_[rank].Z(control_1); - break; - case 5: - states_[rank].CZ(control_1, control_2); - states_[rank].CX(control_2, target); - states_[rank].Z(control_2); - break; - case 6: - states_[rank].CX(control_1, target); - states_[rank].CX(control_2, target); - states_[rank].X(target); - break; - case 7: - states_[rank].CZ(control_1, control_2); - states_[rank].CX(control_1, target); - states_[rank].CX(control_2, target); - states_[rank].Z(control_1); - states_[rank].Z(control_2); - states_[rank].X(target); - coefficients_[rank] *= -1; //Additional phase - break; - default: //Identity - break; - } -} - -void Runner::apply_ccz(uint_t control_1, uint_t control_2, uint_t target, uint_t branch, int rank) -{ - switch(branch) //Decomposition of the CCZ gate into Cliffords - { - case 1: - states_[rank].CZ(control_1, control_2); - break; - case 2: - states_[rank].CZ(control_1, target); - break; - case 3: - states_[rank].CZ(control_2, target); - break; - case 4: - states_[rank].CZ(control_1, control_2); - states_[rank].CZ(control_1, target); - states_[rank].Z(control_1); - break; - case 5: - states_[rank].CZ(control_1, control_2); - states_[rank].CZ(control_2, target); - states_[rank].Z(control_2); - break; - case 6: - states_[rank].CZ(control_1, target); - states_[rank].CZ(control_2, target); - states_[rank].Z(target); - break; - case 7: - states_[rank].CZ(control_1, control_2); - states_[rank].CZ(control_1, target); - states_[rank].CZ(control_2, target); - states_[rank].Z(control_1); - states_[rank].Z(control_2); - states_[rank].Z(target); - coefficients_[rank] *= -1; //Additional phase - break; - default: //Identity - break; + switch (branch.second) { + case Gates::s: + states_[rank].S(qubit); + break; + case Gates::sdg: + states_[rank].Sdag(qubit); + break; + case Gates::z: + states_[rank].Z(qubit); + break; + default: + break; + } +} + +void Runner::apply_ccx(uint_t control_1, uint_t control_2, uint_t target, + uint_t branch, int rank) { + switch (branch) // Decomposition of the CCX gate into Cliffords + { + case 1: + states_[rank].CZ(control_1, control_2); + break; + case 2: + states_[rank].CX(control_1, target); + break; + case 3: + states_[rank].CX(control_2, target); + break; + case 4: + states_[rank].CZ(control_1, control_2); + states_[rank].CX(control_1, target); + states_[rank].Z(control_1); + break; + case 5: + states_[rank].CZ(control_1, control_2); + states_[rank].CX(control_2, target); + states_[rank].Z(control_2); + break; + case 6: + states_[rank].CX(control_1, target); + states_[rank].CX(control_2, target); + states_[rank].X(target); + break; + case 7: + states_[rank].CZ(control_1, control_2); + states_[rank].CX(control_1, target); + states_[rank].CX(control_2, target); + states_[rank].Z(control_1); + states_[rank].Z(control_2); + states_[rank].X(target); + coefficients_[rank] *= -1; // Additional phase + break; + default: // Identity + break; + } +} + +void Runner::apply_ccz(uint_t control_1, uint_t control_2, uint_t target, + uint_t branch, int rank) { + switch (branch) // Decomposition of the CCZ gate into Cliffords + { + case 1: + states_[rank].CZ(control_1, control_2); + break; + case 2: + states_[rank].CZ(control_1, target); + break; + case 3: + states_[rank].CZ(control_2, target); + break; + case 4: + states_[rank].CZ(control_1, control_2); + states_[rank].CZ(control_1, target); + states_[rank].Z(control_1); + break; + case 5: + states_[rank].CZ(control_1, control_2); + states_[rank].CZ(control_2, target); + states_[rank].Z(control_2); + break; + case 6: + states_[rank].CZ(control_1, target); + states_[rank].CZ(control_2, target); + states_[rank].Z(target); + break; + case 7: + states_[rank].CZ(control_1, control_2); + states_[rank].CZ(control_1, target); + states_[rank].CZ(control_2, target); + states_[rank].Z(control_1); + states_[rank].Z(control_2); + states_[rank].Z(target); + coefficients_[rank] *= -1; // Additional phase + break; + default: // Identity + break; } } //------------------------------------------------------------------------- -//Measurement +// Measurement //------------------------------------------------------------------------- -double Runner::norm_estimation(uint_t n_samples, uint_t repetitions, AER::RngEngine &rng) -{ +double Runner::norm_estimation(uint_t n_samples, uint_t repetitions, + AER::RngEngine &rng) { const int_t NSAMPLES = n_samples; const int_t NQUBITS = n_qubits_; std::vector xi_samples(repetitions, 0.); - for(uint_t m=0; m < repetitions; m++) - { + for (uint_t m = 0; m < repetitions; m++) { std::vector adiag_1(n_samples, 0ULL); std::vector adiag_2(n_samples, 0ULL); - std::vector< std::vector > a(n_samples, std::vector(n_qubits_, 0ULL)); - #pragma omp parallel if (num_threads_ > 1) num_threads(num_threads_) + std::vector> a(n_samples, + std::vector(n_qubits_, 0ULL)); +#pragma omp parallel if (num_threads_ > 1) num_threads(num_threads_) { - #ifdef _WIN32 - #pragma omp for - #else - #pragma omp for collapse(2) - #endif - for (int_t l=0; l generators, AER::RngEngine &rng) -{ +double Runner::norm_estimation(uint_t n_samples, uint_t repetitions, + std::vector generators, + AER::RngEngine &rng) { apply_pauli_projector(generators); return norm_estimation(n_samples, repetitions, rng); } -uint_t Runner::ne_single_sample(uint_t default_samples, - uint_t repetitions, - bool preserve_states, - const AER::reg_t &qubits, - AER::RngEngine &rng) -{ +uint_t Runner::ne_single_sample(uint_t default_samples, uint_t repetitions, + bool preserve_states, const AER::reg_t &qubits, + AER::RngEngine &rng) { uint_t n_samples = std::llrint(4 * std::pow(qubits.size(), 2)); if (default_samples > n_samples) - { n_samples = default_samples; - } double denominator = norm_estimation(n_samples, repetitions, rng); std::vector generators; std::vector states_cache(states_); uint_t out_string = ZERO; - for (uint_t i=0; i Runner::ne_probabilities(uint_t default_samples, uint_t repetitions, - const AER::reg_t &qubits, AER::RngEngine &rng) -{ +std::vector Runner::ne_probabilities(uint_t default_samples, + uint_t repetitions, + const AER::reg_t &qubits, + AER::RngEngine &rng) { uint_t n_probs = 1ULL << qubits.size(); std::vector probs(n_probs, 0.); std::vector states_cache(states_); std::vector generators; - for (uint_t i=0; i> j) & 1ULL) { + for (uint_t j = 0; j < qubits.size(); j++) + if ((i >> j) & 1ULL) generators[j].e = 2; - } else - { generators[j].e = 0; - } - } - double p_estimate = norm_estimation(default_samples, repetitions, generators, rng); + double p_estimate = + norm_estimation(default_samples, repetitions, generators, rng); probs[i] = p_estimate / norm; states_ = states_cache; } return probs; } -uint_t Runner::metropolis_estimation(uint_t n_steps, AER::RngEngine &rng) -{ +uint_t Runner::metropolis_estimation(uint_t n_steps, AER::RngEngine &rng) { init_metropolis(rng); - for (uint_t i=0; i Runner::metropolis_estimation(uint_t n_steps, uint_t n_shots, AER::RngEngine &rng) -{ +std::vector Runner::metropolis_estimation(uint_t n_steps, + uint_t n_shots, + AER::RngEngine &rng) { std::vector shots(n_shots, zer); shots[0] = metropolis_estimation(n_steps, rng); - for (uint_t i=1; i omp_threshold_ && num_threads_ > 1) num_threads(num_threads_) reduction(+:local_real) reduction(+:local_imag) - for (int_t i=0; i omp_threshold_ && num_threads_ > 1) num_threads(num_threads_) reduction(+:local_real) reduction(+:local_imag) + for (int_t i = 0; i < END; i++) { scalar_t amp = states_[i].Amplitude(x_string_); - if(amp.eps == 1) - { + if (amp.eps == 1) { complex_t local = (amp.to_complex() * coefficients_[i]); local_real += local.real(); local_imag += local.imag(); @@ -623,39 +547,29 @@ void Runner::init_metropolis(AER::RngEngine &rng) old_ampsum_ = complex_t(local_real, local_imag); } -void Runner::metropolis_step(AER::RngEngine &rng) -{ +void Runner::metropolis_step(AER::RngEngine &rng) { uint_t proposal = rng.rand(0ULL, n_qubits_); - if(accept_) - { + if (accept_) x_string_ ^= (one << last_proposal_); - } - double real_part = 0.,imag_part =0.; - if (accept_ == 0) - { + double real_part = 0., imag_part = 0.; + if (accept_ == 0) { const int_t END = num_states_; - #pragma omp parallel for if(num_states_ > omp_threshold_ && num_threads_ > 1) num_threads(num_threads_) reduction(+:real_part) reduction(+:imag_part) - for (int_t i=0; i omp_threshold_ && num_threads_ > 1) num_threads(num_threads_) reduction(+:real_part) reduction(+:imag_part) + for (int_t i = 0; i < END; i++) { scalar_t amp = states_[i].ProposeFlip(proposal); - if(amp.eps == 1) - { + if (amp.eps == 1) { complex_t local = (amp.to_complex() * coefficients_[i]); real_part += local.real(); imag_part += local.imag(); } } - } - else - { + } else { const int_t END = num_states_; - #pragma omp parallel for if(num_states_ > omp_threshold_ && num_threads_ > 1) num_threads(num_threads_) reduction(+:real_part) reduction(+:imag_part) - for (int_t i=0; i omp_threshold_ && num_threads_ > 1) num_threads(num_threads_) reduction(+:real_part) reduction(+:imag_part) + for (int_t i = 0; i < END; i++) { states_[i].AcceptFlip(); scalar_t amp = states_[i].ProposeFlip(proposal); - if(amp.eps == 1) - { + if (amp.eps == 1) { complex_t local = (amp.to_complex() * coefficients_[i]); real_part += local.real(); imag_part += local.imag(); @@ -663,58 +577,49 @@ void Runner::metropolis_step(AER::RngEngine &rng) } } complex_t ampsum(real_part, imag_part); - double p_threshold = std::norm(ampsum)/std::norm(old_ampsum_); - if(std::isinf(p_threshold) || std::isnan(p_threshold)) - { + double p_threshold = std::norm(ampsum) / std::norm(old_ampsum_); + if (std::isinf(p_threshold) || std::isnan(p_threshold)) { accept_ = true; old_ampsum_ = ampsum; - last_proposal_ = proposal; //We try to move away from node with 0 probability. - } - else - { + last_proposal_ = + proposal; // We try to move away from node with 0 probability. + } else { double rand = rng.rand(); - if (rand < p_threshold) - { + if (rand < p_threshold) { accept_ = true; old_ampsum_ = ampsum; last_proposal_ = proposal; - } - else - { + } else { accept_ = false; } } } -uint_t Runner::stabilizer_sampler(AER::RngEngine &rng) -{ - uint_t max = (1ULL << n_qubits_) -1; +uint_t Runner::stabilizer_sampler(AER::RngEngine &rng) { + uint_t max = (1ULL << n_qubits_) - 1; return states_[0].Sample(rng.rand_int(ZERO, max)); } -std::vector Runner::stabilizer_sampler(uint_t n_shots, AER::RngEngine &rng) -{ - if(num_states_ > 1) - { - throw std::invalid_argument("CH::Runner::stabilizer_sampler: This method can only be used for a single Stabilizer state.\n"); +std::vector Runner::stabilizer_sampler(uint_t n_shots, + AER::RngEngine &rng) { + if (num_states_ > 1) { + throw std::invalid_argument( + "CH::Runner::stabilizer_sampler: This method can only be used for a " + "single Stabilizer state.\n"); } std::vector shots; shots.reserve(n_shots); - for(uint_t i=0; i omp_threshold_ && num_threads_ > 1) num_threads(num_threads_) reduction(+:real_part) reduction(+:imag_part) - for(int_t i=0; i omp_threshold_ && num_threads_ > 1) num_threads(num_threads_) reduction(+:real_part) reduction(+:imag_part) + for (int_t i = 0; i < END; i++) { complex_t amplitude = states_[i].Amplitude(x_measure).to_complex(); amplitude *= coefficients_[i]; real_part += amplitude.real(); @@ -723,24 +628,21 @@ complex_t Runner::amplitude(uint_t x_measure) return {real_part, imag_part}; } -AER::Vector Runner::statevector() -{ +AER::Vector Runner::statevector() { uint_t ceil = 1ULL << n_qubits_; AER::Vector svector(ceil, false); double norm_squared = 0; - for(uint_t i=0; i Runner::statevector() // Utility //========================================================================= -inline void to_json(json_t &js, const Runner &rn) -{ +inline void to_json(json_t &js, const Runner &rn) { js["num_qubits"] = rn.get_n_qubits(); js["num_states"] = rn.get_num_states(); js["decomposition"] = rn.serialize_decomposition(); } -std::vector Runner::serialize_decomposition() const -{ +std::vector Runner::serialize_decomposition() const { std::vector serialized_states(num_states_); const int_t END = num_states_; - #pragma omp parallel for if(num_threads_ > 1 && num_states_ > omp_threshold_) num_threads(num_threads_) - for(int_t i=0; i 1 && num_states_ > omp_threshold_) \ + num_threads(num_threads_) + for (int_t i = 0; i < END; i++) serialized_states[i] = serialize_state(i).dump(); - } return serialized_states; } -json_t Runner::serialize_state(uint_t rank) const -{ +json_t Runner::serialize_state(uint_t rank) const { json_t js = json_t::object(); std::vector gamma; std::vector M; @@ -780,10 +678,8 @@ json_t Runner::serialize_state(uint_t rank) const G = states_[rank].GMatrix(); uint_t gamma1 = states_[rank].Gamma1(); uint_t gamma2 = states_[rank].Gamma2(); - for(uint_t i=0; i> i) & 1ULL) + 2*((gamma2 >> i) & 1ULL)); - } + for (uint_t i = 0; i < n_qubits_; i++) + gamma.push_back(((gamma1 >> i) & 1ULL) + 2 * ((gamma2 >> i) & 1ULL)); js["gamma"] = gamma; js["M"] = M; js["F"] = F; @@ -793,5 +689,5 @@ json_t Runner::serialize_state(uint_t rank) const return js; } -} // Close namespace CHSimulator +} // namespace CHSimulator #endif diff --git a/src/simulators/extended_stabilizer/chlib/chstabilizer.hpp b/src/simulators/extended_stabilizer/chlib/chstabilizer.hpp index e55c3988da..72f52d70e8 100644 --- a/src/simulators/extended_stabilizer/chlib/chstabilizer.hpp +++ b/src/simulators/extended_stabilizer/chlib/chstabilizer.hpp @@ -26,110 +26,85 @@ #include "core.hpp" #include "framework/utils.hpp" -namespace CHSimulator -{ -// Clifford simulator based on the CH-form -class StabilizerState -{ +namespace CHSimulator { +// Clifford simulator based on the CH-form +class StabilizerState { public: - // constructor creates basis state |phi>=|00...0> + // constructor creates basis state |phi>=|00...0> StabilizerState(const unsigned n_qubits); - //Copy constructor - StabilizerState(const StabilizerState& rhs); + // Copy constructor + StabilizerState(const StabilizerState &rhs); // n = number of qubits // qubits are numbered as q=0,1,...,n-1 - - uint_fast64_t NQubits() const - { - return n; - } - scalar_t Omega() const - { - return omega; - } - uint_fast64_t Gamma1() const - { - return gamma1; - } - uint_fast64_t Gamma2() const - { - return gamma2; - } - std::vector GMatrix() const - { - return G; - } - std::vector FMatrix() const - { - return F; - } - std::vector MMatrix() const - { - return M; - } - + uint_fast64_t NQubits() const { return n; } + scalar_t Omega() const { return omega; } + uint_fast64_t Gamma1() const { return gamma1; } + uint_fast64_t Gamma2() const { return gamma2; } + std::vector GMatrix() const { return G; } + std::vector FMatrix() const { return F; } + std::vector MMatrix() const { return M; } // Clifford gates void CX(unsigned q, unsigned r); // q=control, r=target void CZ(unsigned q, unsigned r); - void H(unsigned q); - void S(unsigned q); - void Z(unsigned q); + void H(unsigned q); + void S(unsigned q); + void Z(unsigned q); void X(unsigned q); void Y(unsigned q); void Sdag(unsigned q); // S inverse - // state preps - void CompBasisVector(uint_fast64_t x); // prepares |x> + void CompBasisVector(uint_fast64_t x); // prepares |x> void HadamardBasisVector(uint_fast64_t x); // prepares H(x)|00...0> - // measurements scalar_t Amplitude(uint_fast64_t x); // computes the amplitude uint_fast64_t Sample(); // returns a sample from the distribution ||^2 uint_fast64_t Sample(uint_fast64_t v_mask); - void MeasurePauli(const pauli_t P); // applies a gate (I+P)/2 - // where P is an arbitrary Pauli operator - void MeasurePauliProjector(const std::vector& generators); + void MeasurePauli(const pauli_t P); // applies a gate (I+P)/2 + // where P is an arbitrary Pauli operator + void MeasurePauliProjector(const std::vector &generators); + + inline scalar_t ScalarPart() { return omega; } - inline scalar_t ScalarPart() {return omega;} + // InnerProduct & Norm Estimation + scalar_t InnerProduct(const uint_fast64_t &A_diag1, + const uint_fast64_t &A_diag2, + const std::vector &A); - - //InnerProduct & Norm Estimation - scalar_t InnerProduct(const uint_fast64_t& A_diag1, const uint_fast64_t& A_diag2, const std::vector& A); - // Metropolis updates: // To initialize Metropolis by a state x call Amplitude(x) - scalar_t ProposeFlip(unsigned flip_pos); // returns the amplitude + scalar_t ProposeFlip(unsigned flip_pos); // returns the amplitude // where x'=bitflip(x,q) // x = current Metropolis state - inline void AcceptFlip() {P=Q;} // accept the proposed bit flip - - friend double NormEstimate(std::vector& states, - const std::vector< std::complex >& phases, - const std::vector& Samples_d1, - const std::vector &Samples_d2, - const std::vector< std::vector >& Samples); - #ifdef _OPENMP - friend double ParallelNormEstimate(std::vector& states, - const std::vector< std::complex >& phases, - const std::vector& Samples_d1, - const std::vector &Samples_d2, - const std::vector< std::vector >& Samples, - int n_threads); - #endif + inline void AcceptFlip() { P = Q; } // accept the proposed bit flip + + friend double + NormEstimate(std::vector &states, + const std::vector> &phases, + const std::vector &Samples_d1, + const std::vector &Samples_d2, + const std::vector> &Samples); +#ifdef _OPENMP + friend double + ParallelNormEstimate(std::vector &states, + const std::vector> &phases, + const std::vector &Samples_d1, + const std::vector &Samples_d2, + const std::vector> &Samples, + int n_threads); +#endif private: - - unsigned n; + unsigned n; // define CH-data (F,G,M,gamma,v,s,omega) see Section IIIA for notations // stabilizer tableaux of the C-layer - uint_fast64_t gamma1; //phase vector gamma = gamma1 + 2*gamma2 (mod 4) - uint_fast64_t gamma2; + uint_fast64_t gamma1; // phase vector gamma = gamma1 + 2*gamma2 (mod 4) + uint_fast64_t gamma2; // each column of F,G,M is represented by uint_fast64_t integer std::vector F; // F[i] = i-th column of F std::vector G; // G[i] = i-th column of G @@ -137,7 +112,7 @@ class StabilizerState uint_fast64_t v; // H-layer U_H=H(v) uint_fast64_t s; // initial state |s> - scalar_t omega; // scalar factor such that |phi>=omega*U_C*U_H|s> + scalar_t omega; // scalar factor such that |phi>=omega*U_C*U_H|s> // multiplies the C-layer on the right by a C-type gate void RightCX(unsigned q, unsigned r); // q=control, r=target @@ -145,16 +120,15 @@ class StabilizerState void RightS(unsigned q); void RightZ(unsigned q); void RightSdag(unsigned q); - - // computes a Pauli operator U_C^{-1}X(x)U_C + + // computes a Pauli operator U_C^{-1}X(x)U_C pauli_t GetPauliX(uint_fast64_t x); // replace the initial state |s> in the CH-form by a superposition // (|t> + i^b |u>)*sqrt(1/2) as described in Proposition 3 - // update the CH-form + // update the CH-form void UpdateSvector(uint_fast64_t t, uint_fast64_t u, unsigned b); - // F-transposed and M-transposed. Need this to compute amplitudes std::vector FT; // FT[i] = i-th row of F std::vector MT; // MT[i] = i-th row of M @@ -164,14 +138,12 @@ class StabilizerState bool isReadyFT; // true if F-transposed is available bool isReadyMT; // true if M-transposed is available - // auxiliary Pauli operators pauli_t P; pauli_t Q; - }; -using cdouble = std::complex; +using cdouble = std::complex; //-------------------------------// // Implementation // @@ -179,328 +151,263 @@ using cdouble = std::complex; // Clifford simulator based on the CH-form for for n<=64 qubits -StabilizerState::StabilizerState(const unsigned n_qubits): -n(n_qubits), // number of qubits -gamma1(zer), -gamma2(zer), -F(n,zer), -G(n,zer), -M(n,zer), -v(zer), -s(zer), -omega(), -FT(n,zer), -MT(n,zer), -isReadyFT(false), -isReadyMT(false) -{ +StabilizerState::StabilizerState(const unsigned n_qubits) + : n(n_qubits), // number of qubits + gamma1(zer), gamma2(zer), F(n, zer), G(n, zer), M(n, zer), v(zer), s(zer), + omega(), FT(n, zer), MT(n, zer), isReadyFT(false), isReadyMT(false) { // initialize by the basis vector 00...0 - if(n>63) - { - throw std::invalid_argument("The CH simulator only supports up to 63 qubits.\n"); + if (n > 63) { + throw std::invalid_argument( + "The CH simulator only supports up to 63 qubits.\n"); } // G and F are identity matrices - for (unsigned q=0; q>q) & one)*C; - } +void StabilizerState::S(unsigned q) { + isReadyMT = false; // we are going to change M + uint_fast64_t C = (one << q); + for (unsigned p = 0; p < n; p++) + M[p] ^= ((G[p] >> q) & one) * C; // update phase vector: gamma[q] gets gamma[q] - 1 - gamma1^=C; - gamma2^=((gamma1 >> q) & one )*C; + gamma1 ^= C; + gamma2 ^= ((gamma1 >> q) & one) * C; } -void StabilizerState::Sdag(unsigned q) -{ - isReadyMT=false;// we are going to change M - uint_fast64_t C=(one<>q) & one)*C; - } +void StabilizerState::Sdag(unsigned q) { + isReadyMT = false; // we are going to change M + uint_fast64_t C = (one << q); + for (unsigned p = 0; p < n; p++) + M[p] ^= ((G[p] >> q) & one) * C; // update phase vector: gamma[q] gets gamma[q] + 1 - gamma2^=((gamma1 >> q) & one )*C; - gamma1^=C; + gamma2 ^= ((gamma1 >> q) & one) * C; + gamma1 ^= C; } -void StabilizerState::Z(unsigned q) -{ +void StabilizerState::Z(unsigned q) { // update phase vector: gamma[q] gets gamma[q] + 2 - gamma2^=(one<> q)&one) + 4*((gamma2>>q) & one); - //Commute the z_string through the hadamard layer - // Each z that hits a hadamard becomes a Pauli X + // Initial phase correction + int phase = 2 * ((gamma1 >> q) & one) + 4 * ((gamma2 >> q) & one); + // Commute the z_string through the hadamard layer + // Each z that hits a hadamard becomes a Pauli X s ^= (z_string & v); - //Remaining Z gates add a global phase from their action on the s string - phase += 4*(AER::Utils::hamming_parity((z_string & ~(v)) & s)); - //Commute the x_string through the hadamard layer - // Any remaining X gates update s + // Remaining Z gates add a global phase from their action on the s string + phase += 4 * (AER::Utils::hamming_parity((z_string & ~(v)) & s)); + // Commute the x_string through the hadamard layer + // Any remaining X gates update s s ^= (x_string & ~(v)); - //New z gates add a global phase from their action on the s string - phase += 4*(AER::Utils::hamming_parity((x_string & v) & s)); - //Update the global phase - omega.e = (omega.e + phase)%8; + // New z gates add a global phase from their action on the s string + phase += 4 * (AER::Utils::hamming_parity((x_string & v) & s)); + // Update the global phase + omega.e = (omega.e + phase) % 8; } -void StabilizerState::Y(unsigned q) -{ +void StabilizerState::Y(unsigned q) { Z(q); X(q); - //Add a global phase of -i - omega.e = (omega.e + 2)%8; + // Add a global phase of -i + omega.e = (omega.e + 2) % 8; } -void StabilizerState::RightCX(unsigned q, unsigned r) -{ - G[q]^=G[r]; - F[r]^=F[q]; - M[q]^=M[r]; +void StabilizerState::RightCX(unsigned q, unsigned r) { + G[q] ^= G[r]; + F[r] ^= F[q]; + M[q] ^= M[r]; } -void StabilizerState::CX(unsigned q, unsigned r) -{ - isReadyMT=false;// we are going to change M and F - isReadyFT=false; - uint_fast64_t C=(one<>q) & one)*T; - F[p]^=((F[p]>>r) & one)*C; - M[p]^=((M[p]>>r) & one)*C; +void StabilizerState::CX(unsigned q, unsigned r) { + isReadyMT = false; // we are going to change M and F + isReadyFT = false; + uint_fast64_t C = (one << q); + uint_fast64_t T = (one << r); + bool b = false; + for (unsigned p = 0; p < n; p++) { + b = (b != ((M[p] & C) && (F[p] & T))); + G[p] ^= ((G[p] >> q) & one) * T; + F[p] ^= ((F[p] >> r) & one) * C; + M[p] ^= ((M[p] >> r) & one) * C; } // update phase vector as // gamma[q] gets gamma[q] + gamma[r] + 2*b (mod 4) - if (b) gamma2^=C; - b= ( (gamma1 & C) && (gamma1 & T)); - gamma1^=((gamma1 >> r) & one)*C; - gamma2^=((gamma2 >> r) & one)*C; - if (b) gamma2^=C; + if (b) + gamma2 ^= C; + b = ((gamma1 & C) && (gamma1 & T)); + gamma1 ^= ((gamma1 >> r) & one) * C; + gamma2 ^= ((gamma2 >> r) & one) * C; + if (b) + gamma2 ^= C; } -void StabilizerState::RightCZ(unsigned q, unsigned r) -{ - isReadyMT=false;// we are going to change M - M[q]^=F[r]; - M[r]^=F[q]; - gamma2^=(F[q] & F[r]); +void StabilizerState::RightCZ(unsigned q, unsigned r) { + isReadyMT = false; // we are going to change M + M[q] ^= F[r]; + M[r] ^= F[q]; + gamma2 ^= (F[q] & F[r]); } -void StabilizerState::CZ(unsigned q, unsigned r) -{ - isReadyMT=false;// we are going to change M - uint_fast64_t C=(one<>r) & one)*C; - M[p]^=((G[p]>>q) & one)*T; +void StabilizerState::CZ(unsigned q, unsigned r) { + isReadyMT = false; // we are going to change M + uint_fast64_t C = (one << q); + uint_fast64_t T = (one << r); + for (unsigned p = 0; p < n; p++) { + M[p] ^= ((G[p] >> r) & one) * C; + M[p] ^= ((G[p] >> q) & one) * T; } } -pauli_t StabilizerState::GetPauliX(uint_fast64_t x) -{ +pauli_t StabilizerState::GetPauliX(uint_fast64_t x) { // make sure that M-transposed and F-transposed have been already computed if (!isReadyMT) - { TransposeM(); - } if (!isReadyFT) - { TransposeF(); - } pauli_t R; - for (unsigned pos=0; pos>pos) & one); - P1.e+=2*((gamma2>>pos) & one); - P1.X=FT[pos]; - P1.Z=MT[pos]; - R*=P1; + P1.e = 1 * ((gamma1 >> pos) & one); + P1.e += 2 * ((gamma2 >> pos) & one); + P1.X = FT[pos]; + P1.Z = MT[pos]; + R *= P1; } } - return R; + return R; } -scalar_t StabilizerState::Amplitude(uint_fast64_t x) -{ +scalar_t StabilizerState::Amplitude(uint_fast64_t x) { // compute transposed matrices if needed if (!isReadyMT) - { TransposeM(); - } if (!isReadyFT) - { TransposeF(); - } // compute Pauli U_C^{-1} X(x) U_C - P=GetPauliX(x); + P = GetPauliX(x); - if (!omega.eps) return omega; // the state is zero + if (!omega.eps) + return omega; // the state is zero // now the amplitude = complex conjugate of // Z-part of P is absorbed into 0^n scalar_t amp; - amp.e=2*P.e; - int p = (int) AER::Utils::popcount(v); - amp.p= -1 * p;// each Hadamard gate contributes 1/sqrt(2) - bool isNonZero=true; - - - for (unsigned q=0; q + amp.e = 2 * P.e; + int p = (int)AER::Utils::popcount(v); + amp.p = -1 * p; // each Hadamard gate contributes 1/sqrt(2) + bool isNonZero = true; + + for (unsigned q = 0; q < n; q++) { + uint_fast64_t pos = (one << q); + if (v & pos) { + amp.e += + 4 * ((s & pos) && (P.X & pos)); // minus sign that comes from <1|H|1> + } else { + isNonZero = (((P.X & pos) == (s & pos))); } - else - { - isNonZero=( ((P.X & pos)==(s & pos)) ); - } - if (!isNonZero) break; + if (!isNonZero) + break; } - - amp.e%=8; - if (isNonZero) - { + amp.e %= 8; + if (isNonZero) { amp.conjugate(); - } - else - { - amp.eps=0; + } else { + amp.eps = 0; return amp; } // multiply amp by omega - amp.p+=omega.p; - amp.e=(amp.e + omega.e) % 8; + amp.p += omega.p; + amp.e = (amp.e + omega.e) % 8; return amp; } -scalar_t StabilizerState::ProposeFlip(unsigned flip_pos) -{ +scalar_t StabilizerState::ProposeFlip(unsigned flip_pos) { // Q gets Pauli operator U_C^{-1} X_{flip_pos} U_C - Q.e=1*((gamma1>>flip_pos) & one); - Q.e+=2*((gamma2>>flip_pos) & one); - Q.X=FT[flip_pos]; - Q.Z=MT[flip_pos]; - Q*=P; + Q.e = 1 * ((gamma1 >> flip_pos) & one); + Q.e += 2 * ((gamma2 >> flip_pos) & one); + Q.X = FT[flip_pos]; + Q.Z = MT[flip_pos]; + Q *= P; - if (!omega.eps) return omega; // the state is zero + if (!omega.eps) + return omega; // the state is zero // the rest is the same as Amplitude() except that P becomes Q @@ -510,561 +417,503 @@ scalar_t StabilizerState::ProposeFlip(unsigned flip_pos) // the rest is the same as Amplitude() except that P is replaced by Q scalar_t amp; - amp.e=2*Q.e; - amp.p=-1*(AER::Utils::popcount(v));// each Hadamard gate contributes 1/sqrt(2) - bool isNonZero=true; - - for (unsigned q=0; q + amp.e = 2 * Q.e; + amp.p = -1 * + (AER::Utils::popcount(v)); // each Hadamard gate contributes 1/sqrt(2) + bool isNonZero = true; + + for (unsigned q = 0; q < n; q++) { + uint_fast64_t pos = (one << q); + if (v & pos) { + amp.e += + 4 * ((s & pos) && (Q.X & pos)); // minus sign that comes from <1|H|1> + } else { + isNonZero = (((Q.X & pos) == (s & pos))); } - else - { - isNonZero=( ((Q.X & pos)==(s & pos)) ); - } - if (!isNonZero) break; + if (!isNonZero) + break; } - amp.e%=8; - if (isNonZero) - { + amp.e %= 8; + if (isNonZero) { amp.conjugate(); - } - else - { - amp.eps=0; + } else { + amp.eps = 0; return amp; } // multiply amp by omega - amp.p+=omega.p; - amp.e=(amp.e + omega.e) % 8; + amp.p += omega.p; + amp.e = (amp.e + omega.e) % 8; return amp; } -uint_fast64_t StabilizerState::Sample() -{ - uint_fast64_t x=zer; - for (unsigned q=0; q1 then multiply U_C on the right by the first half of the circuit VC - nu0^=qpos; // set q-th bit to zero - if (nu0) - for (unsigned q1=q+1; q10 then apply the second half of the circuit VC - if (nu1) - for (unsigned q1=0; q11 then apply the circuit VC - nu1^=qpos; - if (nu1) - for (unsigned q1=q+1; q1 = eta^{e1} S^{e2} H^{e3} |e4> - // here eta=exp( i (pi/4) ) - // a=0,1 - // b=0,1,2,3 - // - // H^0 S^0 |+> = eta^0 S^0 H^1 |0> - // H^0 S^1 |+> = eta^0 S^1 H^1 |0> - // H^0 S^2 |+> = eta^0 S^0 H^1 |1> - // H^0 S^3 |+> = eta^0 S^1 H^1 |1> - // - // H^1 S^0 |+> = eta^0 S^0 H^0 |0> - // H^1 S^1 |+> = eta^1 S^1 H^1 |1> - // H^1 S^2 |+> = eta^0 S^0 H^0 |1> - // H^1 S^3 |+> = eta^{7} S^1 H^1 |0> - // - // "analytic" formula: - // e1 = a * (b mod 2) * ( 3*b -2 ) - // e2 = b mod 2 - // e3 = not(a) + a * (b mod 2) - // e4 = not(a)*(b>=2) + a*( (b==1) || (b==2) ) - - bool a=((v & qpos)>0); - unsigned e1=a*(b % 2)*( 3*b -2); - unsigned e2 = b % 2; - bool e3 = ( (!a) != (a && ((b % 2)>0) ) ); - bool e4 = ( ( (!a) && (b>=2) ) != (a && ((b==1) || (b==2))) ); - - // update CH-form - // set q-th bit of s to e4 - s&=~qpos; - s^=e4*qpos; - - // set q-th bit of v to e3 - v&=~qpos; - v^=e3*qpos; - - // update the scalar factor omega - omega.e=(omega.e + e1) % 8; - - // multiply the C-layer on the right by S^{e2} on the q-th qubit - if (e2) RightS(q); -} - - -void StabilizerState::H(unsigned q) -{ - isReadyMT=false;// we are going to change M and F - isReadyFT=false; - // extract the q-th row of F,G,M - uint_fast64_t rowF=zer; - uint_fast64_t rowG=zer; - uint_fast64_t rowM=zer; - for (unsigned j=0; j>q) & one ); - rowG^=pos*( (G[j]>>q) & one ); - rowM^=pos*( (M[j]>>q) & one ); + // now t and u are distinct + isReadyFT = false; // below we are going to change F and M + isReadyMT = false; + // naming of variables roughly follows Section IIIA + uint_fast64_t ut = u ^ t; + uint_fast64_t nu0 = (~v) & ut; + uint_fast64_t nu1 = v & ut; + // + b %= 4; + unsigned q = 0; + uint_fast64_t qpos = zer; + if (nu0) { + + // the subset nu0 is non-empty + // find the first element of nu0 + q = 0; + while (!(nu0 & (one << q))) + q++; + + qpos = (one << q); + if (!(nu0 & qpos)) { + throw std::logic_error( + "Failed to find first bit of nu despite being non-empty."); } - // after commuting H through the C and H laters it maps |s> to a state - // sqrt(0.5)*[ (-1)^alpha |t> + i^{gamma[p]} (-1)^beta |u> ] - // - // compute t,s,alpha,beta - uint_fast64_t t = s ^ (rowG & v); - uint_fast64_t u = s ^ (rowF & (~v)) ^ (rowM & v); - - unsigned alpha = AER::Utils::popcount( rowG & (~v) & s ); - unsigned beta = AER::Utils::popcount( (rowM & (~v) & s) ^ (rowF & v & (rowM ^ s)) ); - - if (alpha % 2) omega.e=(omega.e+4) % 8; - // get the phase gamma[q] - unsigned phase = ((gamma1>>q) & one) + 2*((gamma2>>q) & one); - unsigned b=(phase + 2*alpha + 2*beta) % 4; - - - // now the initial state is sqrt(0.5)*(|t> + i^b |u>) + // if nu0 has size >1 then multiply U_C on the right by the first half of + // the circuit VC + nu0 ^= qpos; // set q-th bit to zero + if (nu0) + for (unsigned q1 = q + 1; q1 < n; q1++) + if (nu0 & (one << q1)) + RightCX(q, q1); + + // if nu1 has size >0 then apply the second half of the circuit VC + if (nu1) + for (unsigned q1 = 0; q1 < n; q1++) + if (nu1 & (one << q1)) + RightCZ(q, q1); + + } // if (nu0) + else { + // if we got here when nu0 is empty + // find the first element of nu1 + q = 0; + while (!(nu1 & (one << q))) + q++; + + qpos = (one << q); + if (!(nu1 & qpos)) + throw std::logic_error("Expect at least one non-zero element in nu1.\n"); + + // if nu1 has size >1 then apply the circuit VC + nu1 ^= qpos; + if (nu1) + for (unsigned q1 = q + 1; q1 < n; q1++) + if (nu1 & (one << q1)) + RightCX(q1, q); + + } // if (nu0) else + + // update the initial state + // if t_q=1 then switch t_q and u_q + // uint_fast64_t y,z; + if (t & qpos) { + s = u; + omega.e = (omega.e + 2 * b) % 8; + b = (4 - b) % 4; + if (!!(u & qpos)) { + } + } else + s = t; + + // change the order of H and S gates + // H^{a} S^{b} |+> = eta^{e1} S^{e2} H^{e3} |e4> + // here eta=exp( i (pi/4) ) + // a=0,1 + // b=0,1,2,3 + // + // H^0 S^0 |+> = eta^0 S^0 H^1 |0> + // H^0 S^1 |+> = eta^0 S^1 H^1 |0> + // H^0 S^2 |+> = eta^0 S^0 H^1 |1> + // H^0 S^3 |+> = eta^0 S^1 H^1 |1> + // + // H^1 S^0 |+> = eta^0 S^0 H^0 |0> + // H^1 S^1 |+> = eta^1 S^1 H^1 |1> + // H^1 S^2 |+> = eta^0 S^0 H^0 |1> + // H^1 S^3 |+> = eta^{7} S^1 H^1 |0> + // + // "analytic" formula: + // e1 = a * (b mod 2) * ( 3*b -2 ) + // e2 = b mod 2 + // e3 = not(a) + a * (b mod 2) + // e4 = not(a)*(b>=2) + a*( (b==1) || (b==2) ) + + bool a = ((v & qpos) > 0); + unsigned e1 = a * (b % 2) * (3 * b - 2); + unsigned e2 = b % 2; + bool e3 = ((!a) != (a && ((b % 2) > 0))); + bool e4 = (((!a) && (b >= 2)) != (a && ((b == 1) || (b == 2)))); + + // update CH-form + // set q-th bit of s to e4 + s &= ~qpos; + s ^= e4 * qpos; + + // set q-th bit of v to e3 + v &= ~qpos; + v ^= e3 * qpos; + + // update the scalar factor omega + omega.e = (omega.e + e1) % 8; + + // multiply the C-layer on the right by S^{e2} on the q-th qubit + if (e2) + RightS(q); +} - // take care of the trivial case - if (t==u) +void StabilizerState::H(unsigned q) { + isReadyMT = false; // we are going to change M and F + isReadyFT = false; + // extract the q-th row of F,G,M + uint_fast64_t rowF = zer; + uint_fast64_t rowG = zer; + uint_fast64_t rowM = zer; + for (unsigned j = 0; j < n; j++) { + uint_fast64_t pos = (one << j); + rowF ^= pos * ((F[j] >> q) & one); + rowG ^= pos * ((G[j] >> q) & one); + rowM ^= pos * ((M[j] >> q) & one); + } + + // after commuting H through the C and H laters it maps |s> to a state + // sqrt(0.5)*[ (-1)^alpha |t> + i^{gamma[p]} (-1)^beta |u> ] + // + // compute t,s,alpha,beta + uint_fast64_t t = s ^ (rowG & v); + uint_fast64_t u = s ^ (rowF & (~v)) ^ (rowM & v); + + unsigned alpha = AER::Utils::popcount(rowG & (~v) & s); + unsigned beta = + AER::Utils::popcount((rowM & (~v) & s) ^ (rowF & v & (rowM ^ s))); + + if (alpha % 2) + omega.e = (omega.e + 4) % 8; + // get the phase gamma[q] + unsigned phase = ((gamma1 >> q) & one) + 2 * ((gamma2 >> q) & one); + unsigned b = (phase + 2 * alpha + 2 * beta) % 4; + + // now the initial state is sqrt(0.5)*(|t> + i^b |u>) + + // take care of the trivial case + if (t == u) { + s = t; + if (!((b == 1) || (b == 3))) // otherwise the state is not normalized { - s=t; - if(!((b==1) || (b==3))) // otherwise the state is not normalized - { - throw std::logic_error("State is not properly normalised, b should be 1 or 3.\n"); - } - if (b==1) - omega.e=(omega.e + 1) % 8; - else - omega.e=(omega.e + 7) % 8; + throw std::logic_error( + "State is not properly normalised, b should be 1 or 3.\n"); } + if (b == 1) + omega.e = (omega.e + 1) % 8; else - UpdateSvector(t,u,b); - + omega.e = (omega.e + 7) % 8; + } else + UpdateSvector(t, u, b); } -void StabilizerState::MeasurePauli(pauli_t PP) -{ - // compute Pauli R = U_C^{-1} P U_C - pauli_t R; - R.e = PP.e; - - for (unsigned j=0; j>j) & one) - { - // multiply R by U_C^{-1} X_j U_C - // extract the j-th rows of F and M - uint_fast64_t rowF=zer; - uint_fast64_t rowM=zer; - for (unsigned i=0; i>j) & one); - rowM^=(one<>j) & one); - } - R.e+= 2*AER::Utils::popcount(R.Z & rowF); // extra sign from Pauli commutation - R.Z^=rowM; - R.X^=rowF; - R.e+= ((gamma1>>j) & one) + 2*((gamma2>>j) & one); - } - for (unsigned q=0; q becomes 0.5*(|s> + R |s>) = 0.5*(|s> + i^b |s ^ R.X>) - unsigned b = (R.e + 2*AER::Utils::popcount(R.Z & s) ) % 4; - UpdateSvector(s, s ^ R.X, b); - // account for the extra factor sqrt(1/2) - omega.p-=1; - - isReadyMT=false;// we have changed change M and F - isReadyFT=false; +void StabilizerState::MeasurePauli(pauli_t PP) { + // compute Pauli R = U_C^{-1} P U_C + pauli_t R; + R.e = PP.e; + + for (unsigned j = 0; j < n; j++) + if ((PP.X >> j) & one) { + // multiply R by U_C^{-1} X_j U_C + // extract the j-th rows of F and M + uint_fast64_t rowF = zer; + uint_fast64_t rowM = zer; + for (unsigned i = 0; i < n; i++) { + rowF ^= (one << i) * ((F[i] >> j) & one); + rowM ^= (one << i) * ((M[i] >> j) & one); + } + R.e += 2 * AER::Utils::popcount( + R.Z & rowF); // extra sign from Pauli commutation + R.Z ^= rowM; + R.X ^= rowF; + R.e += ((gamma1 >> j) & one) + 2 * ((gamma2 >> j) & one); + } + for (unsigned q = 0; q < n; q++) + R.Z ^= (one << q) * (AER::Utils::popcount(PP.Z & G[q]) % 2); + + // now R=U_C^{-1} PP U_C + // next conjugate R by U_H + uint_fast64_t tempX = ((~v) & R.X) ^ (v & R.Z); + uint_fast64_t tempZ = ((~v) & R.Z) ^ (v & R.X); + // the sign flips each time a Hadamard hits Y on some qubit + R.e = (R.e + 2 * AER::Utils::popcount(v & R.X & R.Z)) % 4; + R.X = tempX; + R.Z = tempZ; + + // now the initial state |s> becomes 0.5*(|s> + R |s>) = 0.5*(|s> + i^b |s ^ + // R.X>) + unsigned b = (R.e + 2 * AER::Utils::popcount(R.Z & s)) % 4; + UpdateSvector(s, s ^ R.X, b); + // account for the extra factor sqrt(1/2) + omega.p -= 1; + + isReadyMT = false; // we have changed change M and F + isReadyFT = false; } -void StabilizerState::MeasurePauliProjector(const std::vector& generators) -//Measure generators of a projector. +void StabilizerState::MeasurePauliProjector( + const std::vector &generators) +// Measure generators of a projector. { - for (uint_fast64_t i=0; iMeasurePauli(generators[i]); - if (omega.eps == 0) - { - break; - } - } + for (uint_fast64_t i = 0; i < generators.size(); i++) { + this->MeasurePauli(generators[i]); + if (omega.eps == 0) + break; + } } -scalar_t StabilizerState::InnerProduct(const uint_fast64_t& A_diag1, - const uint_fast64_t& A_diag2, - const std::vector& A) -{ - uint_fast64_t K_diag1 = zer, K_diag2=zer, J_diag1=gamma1, J_diag2=gamma2; - std::vector J(n, zer); - std::vector K(n, zer); - std::vector placeholder(n, zer); - if(!isReadyMT) - { - TransposeM(); - } - if (!isReadyFT) - { - TransposeF(); - } - //Setup the J matrix - for (size_t i=0; i &A) { + uint_fast64_t K_diag1 = zer, K_diag2 = zer, J_diag1 = gamma1, + J_diag2 = gamma2; + std::vector J(n, zer); + std::vector K(n, zer); + std::vector placeholder(n, zer); + if (!isReadyMT) + TransposeM(); + if (!isReadyFT) + TransposeF(); + // Setup the J matrix + for (size_t i = 0; i < n; i++) { + for (size_t j = i; j < n; j++) { + if (AER::Utils::hamming_parity(MT[i] & FT[j])) { + J[i] |= (one << j); + J[j] |= (one << i); + } } - //placeholder = J*G) - for (size_t i=0; i> r) & one) - { - uint_fast64_t one_bit = (J_diag1 >> r) & one; - if ((K_diag1 >> i) & one_bit) - { - K_diag2 ^= shift; - } - K_diag1 ^= one_bit*shift; - K_diag2 ^= ((J_diag2 >> r)&one) * shift; - } - for (size_t k=r+1; k> r) & (col_i >> r) & (col_i >> k) & one) - { - K_diag2 ^= shift; - } - } - } + for (size_t r = 0; r < n; r++) { + if ((col_i >> r) & one) { + uint_fast64_t one_bit = (J_diag1 >> r) & one; + if ((K_diag1 >> i) & one_bit) + K_diag2 ^= shift; + K_diag1 ^= one_bit * shift; + K_diag2 ^= ((J_diag2 >> r) & one) * shift; + } + for (size_t k = r + 1; k < n; k++) + if ((J[k] >> r) & (col_i >> r) & (col_i >> k) & one) + K_diag2 ^= shift; } - unsigned col=0; - QuadraticForm q(AER::Utils::popcount(v)); - //We need to setup a quadratic form to evaluate the Exponential Sum - for (size_t i=0; i>i) & one) - { - uint_fast64_t shift = (one << col); - //D = Diag(K(1,1)) + 2*[s + s*K](1) - // J = K(1,1); - q.D1 ^= ((K_diag1 >> i) & one) * shift; - q.D2 ^= (( ((K_diag2 >> i) ^ (s >> i)) & one) - ^ AER::Utils::hamming_parity(K[i] & s)) * shift; - // q.D2 ^= ((s >> i) & one) * shift; - // q.D2 ^= AER::Utils::hamming_parity(K[i] & s) * shift; - unsigned row=0; - for (size_t j=0; j>j) & one) - { - q.J[col] |= (((K[i] >> j) & one) << row); - row++; - } - } - col++; + } + unsigned col = 0; + QuadraticForm q(AER::Utils::popcount(v)); + // We need to setup a quadratic form to evaluate the Exponential Sum + for (size_t i = 0; i < n; i++) { + if ((v >> i) & one) { + uint_fast64_t shift = (one << col); + // D = Diag(K(1,1)) + 2*[s + s*K](1) + // J = K(1,1); + q.D1 ^= ((K_diag1 >> i) & one) * shift; + q.D2 ^= ((((K_diag2 >> i) ^ (s >> i)) & one) ^ + AER::Utils::hamming_parity(K[i] & s)) * + shift; + // q.D2 ^= ((s >> i) & one) * shift; + // q.D2 ^= AER::Utils::hamming_parity(K[i] & s) * shift; + unsigned row = 0; + for (size_t j = 0; j < n; j++) { + if ((v >> j) & one) { + q.J[col] |= (((K[i] >> j) & one) << row); + row++; } + } + col++; } - // Q = 4* (s.v) + sKs - q.Q = AER::Utils::hamming_parity(s&v)*4; - for (size_t i=0; i>i) & one) - { - q.Q = (q.Q + 4*((K_diag2 >> i) & one) + 2*((K_diag1 >> i) & one))%8;// + 2*((K_diag1 >> i) & one))%8; - for (size_t j=i+1; j>j) & (K[j] >> i) & one) - { - q.Q ^= 4; - } - } - } + } + // Q = 4* (s.v) + sKs + q.Q = AER::Utils::hamming_parity(s & v) * 4; + for (size_t i = 0; i < n; i++) { + if ((s >> i) & one) { + q.Q = (q.Q + 4 * ((K_diag2 >> i) & one) + 2 * ((K_diag1 >> i) & one)) % + 8; // + 2*((K_diag1 >> i) & one))%8; + for (size_t j = i + 1; j < n; j++) + if ((s >> j) & (K[j] >> i) & one) + q.Q ^= 4; } - scalar_t amp = q.ExponentialSum(); - // Reweight by 2^{-(n+|v|)}/2 - amp.p -= (n+AER::Utils::popcount(v)); - // We need to further multiply by omega* - scalar_t psi_amp(omega); - psi_amp.conjugate(); - amp *= psi_amp; - return amp; + } + scalar_t amp = q.ExponentialSum(); + // Reweight by 2^{-(n+|v|)}/2 + amp.p -= (n + AER::Utils::popcount(v)); + // We need to further multiply by omega* + scalar_t psi_amp(omega); + psi_amp.conjugate(); + amp *= psi_amp; + return amp; } -double NormEstimate(std::vector& states, const std::vector< std::complex >& phases, - const std::vector& Samples_d1, const std::vector &Samples_d2, - const std::vector< std::vector >& Samples) -{ - // Norm estimate for a state |psi> = \sum_{i} c_{i}|phi_{i}> - double xi=0; - unsigned L = Samples_d1.size(); - // std::vector data = (L,0.); - for (size_t i=0; i &states, + const std::vector> &phases, + const std::vector &Samples_d1, + const std::vector &Samples_d2, + const std::vector> &Samples) { + // Norm estimate for a state |psi> = \sum_{i} c_{i}|phi_{i}> + double xi = 0; + unsigned L = Samples_d1.size(); + // std::vector data = (L,0.); + for (size_t i = 0; i < L; i++) { + double re_eta = 0., im_eta = 0.; + const int_t END = states.size(); +#pragma omp parallel for reduction(+ : re_eta) reduction(+ : im_eta) + for (int_t j = 0; j < END; j++) { + if (states[j].ScalarPart().eps != 0) { + scalar_t amp = + states[j].InnerProduct(Samples_d1[i], Samples_d2[i], Samples[i]); + if (amp.eps != 0) { + if (amp.e % 2) + amp.p--; + double mag = pow(2, amp.p / (double)2); + cdouble phase(RE_PHASE[amp.e], IM_PHASE[amp.e]); + phase *= conj(phases[j]); + re_eta += (mag * real(phase)); + im_eta += (mag * imag(phase)); + } } - xi += (pow(re_eta,2) + pow(im_eta, 2)); } - return std::pow(2, states[0].NQubits()) * (xi/L); + xi += (pow(re_eta, 2) + pow(im_eta, 2)); + } + return std::pow(2, states[0].NQubits()) * (xi / L); } -double ParallelNormEstimate(std::vector& states, const std::vector< std::complex >& phases, - const std::vector& Samples_d1, const std::vector& Samples_d2, - const std::vector< std::vector >& Samples, int n_threads) -{ - double xi=0; - unsigned L = Samples_d1.size(); - unsigned chi = states.size(); - for (uint_fast64_t i=0; i &states, + const std::vector> &phases, + const std::vector &Samples_d1, + const std::vector &Samples_d2, + const std::vector> &Samples, + int n_threads) { + double xi = 0; + unsigned L = Samples_d1.size(); + unsigned chi = states.size(); + for (uint_fast64_t i = 0; i < L; i++) { + double re_eta = 0., im_eta = 0.; +#pragma omp parallel for reduction(+:re_eta) reduction(+:im_eta) num_threads(n_threads) + for (int_t j = 0; j < chi; j++) { + if (states[j].ScalarPart().eps != 0) { + scalar_t amp = + states[j].InnerProduct(Samples_d1[i], Samples_d2[i], Samples[i]); + if (amp.eps != 0) { + if (amp.e % 2) + amp.p--; + double mag = pow(2, amp.p / (double)2); + cdouble phase(RE_PHASE[amp.e], IM_PHASE[amp.e]); + phase *= conj(phases[j]); + re_eta += (mag * real(phase)); + im_eta += (mag * imag(phase)); } - xi += (pow(re_eta,2) + pow(im_eta, 2)); + } } - return pow(2., states[0].NQubits())*(xi/L); + xi += (pow(re_eta, 2) + pow(im_eta, 2)); + } + return pow(2., states[0].NQubits()) * (xi / L); } -} +} // namespace CHSimulator #endif diff --git a/src/simulators/extended_stabilizer/chlib/core.hpp b/src/simulators/extended_stabilizer/chlib/core.hpp index 1fac013e18..729e0b8930 100644 --- a/src/simulators/extended_stabilizer/chlib/core.hpp +++ b/src/simulators/extended_stabilizer/chlib/core.hpp @@ -17,15 +17,14 @@ #include #include -#include #include +#include #include #include #include "framework/utils.hpp" -namespace CHSimulator -{ +namespace CHSimulator { static const std::array RE_PHASE = {1, 1, 0, -1, -1, -1, 0, 1}; static const std::array IM_PHASE = {0, 1, 1, 1, 0, -1, -1, -1}; @@ -46,91 +45,79 @@ struct scalar_t { int e = 0; // constructor makes number 1 scalar_t() = default; - scalar_t(const scalar_t& rhs) = default; - scalar_t(const std::complex coeff) - { + scalar_t(const scalar_t &rhs) = default; + scalar_t(const std::complex coeff) { double abs_val = std::abs(coeff); - if(std::abs(abs_val-0.)<1e-8) - { + if (std::abs(abs_val - 0.) < 1e-8) { eps = 0; - } - else - { - p = (int) std::log2(std::pow(abs_val, 2)); - bool is_real_zero = (std::abs(coeff.real()-0)<1e-8); - bool is_imag_zero = (std::abs(coeff.imag()-0)<1e-8); + } else { + p = (int)std::log2(std::pow(abs_val, 2)); + bool is_real_zero = (std::abs(coeff.real() - 0) < 1e-8); + bool is_imag_zero = (std::abs(coeff.imag() - 0) < 1e-8); bool is_real_posi = (coeff.real() > 0) && !(is_real_zero); bool is_imag_posi = (coeff.imag() > 0) && !(is_imag_zero); - char switch_var = (is_real_zero * 1) + (is_imag_zero*2) + (is_real_posi * 4) + (is_imag_posi * 8); - switch(switch_var) - { - case 0: - e = 5; - break; - case 1: - e = 6; - break; - case 2: - e = 4; - break; - case 4: - e=7; - break; - case 6: - e=0; - break; - case 8: - e=3; - break; - case 9: - e=2; - break; - case 12: - e=1; - break; - default: - throw std::runtime_error("Unsure what to do here chief."); - break; + char switch_var = (is_real_zero * 1) + (is_imag_zero * 2) + + (is_real_posi * 4) + (is_imag_posi * 8); + switch (switch_var) { + case 0: + e = 5; + break; + case 1: + e = 6; + break; + case 2: + e = 4; + break; + case 4: + e = 7; + break; + case 6: + e = 0; + break; + case 8: + e = 3; + break; + case 9: + e = 2; + break; + case 12: + e = 1; + break; + default: + throw std::runtime_error("Unsure what to do here chief."); + break; } } } // complex conjugation - inline void conjugate() - { - e%=8; - e=(8-e) % 8; + inline void conjugate() { + e %= 8; + e = (8 - e) % 8; } - inline void makeOne() - { - eps=1; - p=0; - e=0; + inline void makeOne() { + eps = 1; + p = 0; + e = 0; } - scalar_t& operator*=(const scalar_t& rhs); - scalar_t operator*(const scalar_t& rhs) const; + scalar_t &operator*=(const scalar_t &rhs); + scalar_t operator*(const scalar_t &rhs) const; - std::complex to_complex() const - { - if (eps==0) - { + std::complex to_complex() const { + if (eps == 0) return {0., 0.}; - } - std::complex mag(std::pow(2, p/(double)2), 0.); + std::complex mag(std::pow(2, p / (double)2), 0.); std::complex phase(RE_PHASE[e], IM_PHASE[e]); - if(e % 2){ + if (e % 2) phase /= std::sqrt(2); - } - return mag*phase; + return mag * phase; } - void Print() - { - std::cout<<"eps="<63) - { +QuadraticForm::QuadraticForm(unsigned int n_qubits) + : n(n_qubits), Q(0), D1(zer), D2(zer), J(n_qubits, zer) { + if (n > 63) throw qubex; - } } -QuadraticForm::QuadraticForm(const QuadraticForm &rhs) : J(rhs.n, zer) -{ - n=rhs.n; - Q=rhs.Q; - D1=rhs.D1; - D2=rhs.D2; - for (size_t i=0; i> i) & one) + ((q.D1 >> i) & one))); - } + for (size_t i = 0; i < q.n; i++) + os << " " << (2 * (2 * ((q.D2 >> i) & one) + ((q.D1 >> i) & one))); os << std::endl; - os <<"J:\n"; - for (size_t i=0; i> j) & one) << " "; - } os << "\n"; } return os; @@ -311,96 +268,80 @@ std::ostream& operator<<(std::ostream& os, const QuadraticForm& q) // Output: integers p>=0, m=0,1,2,..,7, and eps=0,1 such that // Z = sum_x exp(j*(pi/4)*q(x)) = eps*2^(p/2)*exp(j*pi*m/4) // if eps=0 then p and m can be ignored -scalar_t QuadraticForm::ExponentialSum() -{ - +scalar_t QuadraticForm::ExponentialSum() { + // Variables for Z2-valued exponential sums // Z=(2^pow2)*sigma // where pow2 is integer // sigma=0,1 // if Z=0 then pow2,sigma contain junk and isZero=1 - int pow2_real=0; - bool sigma_real=false; - bool isZero_real=false; + int pow2_real = 0; + bool sigma_real = false; + bool isZero_real = false; - int pow2_imag=0; - bool sigma_imag=false; - bool isZero_imag=false; + int pow2_imag = 0; + bool sigma_imag = false; + bool isZero_imag = false; scalar_t amp; amp.makeOne(); - int m = n+1; + int m = n + 1; // Lreal, Limag define the linear part of the Z2-valued forms // M defines the quadratic part of the Z2-valued forms // each column of M is represented by long integer - uint_fast64_t Lreal, c=zer; - std::vector M(n+1,zer); + uint_fast64_t Lreal, c = zer; + std::vector M(n + 1, zer); Lreal = D2; - for (size_t i=0; i> i) & one ) - { + for (size_t i = 0; i < n; i++) { + if ((D1 >> i) & one) { c |= (one << i); M[n] |= (one << i); } } uint_fast64_t Limag = Lreal; Limag ^= (one << n); - for (size_t i=0; i> j) & one); - bool c1 = !!((c>>i) & one); - bool c2 = !!((c>>j) & one); - if (x ^ (c1 & c2)) - { - M[j] ^= (one << i); - } + for (size_t i = 0; i < n; i++) { + for (size_t j = (i + 1); j < n; j++) { + bool x = !!((J[i] >> j) & one); + bool c1 = !!((c >> i) & one); + bool c2 = !!((c >> j) & one); + if (x ^ (c1 & c2)) + M[j] ^= (one << i); } } uint_fast64_t active = zer; active = ~(active); - int nActive=m; - int firstActive=0; - while (nActive>=1) - { + int nActive = m; + int firstActive = 0; + while (nActive >= 1) { // find the first active variable int i1; - for (i1=firstActive; i1> i1) & one) - { break; - } - } - firstActive=i1; + firstActive = i1; // find i2 such that M(i1,i2)!=M(i2,i1) int i2; - bool isFound=false; - for (i2=0; i2>i2) & one) != ((M[i2]>>i1) & one) ); - if (isFound) - { - break; - } + bool isFound = false; + for (i2 = 0; i2 < m; i2++) { + isFound = (((M[i1] >> i2) & one) != ((M[i2] >> i1) & one)); + if (isFound) + break; } - bool L1real = !!(((Lreal >> i1) & one) ^ ((M[i1]>>i1) & one)); - bool L1imag = !!(((Limag >> i1) & one) ^ ((M[i1]>>i1) & one)); + bool L1real = !!(((Lreal >> i1) & one) ^ ((M[i1] >> i1) & one)); + bool L1imag = !!(((Limag >> i1) & one) ^ ((M[i1] >> i1) & one)); // take care of the trivial cases - if (!isFound) - { + if (!isFound) { // the form is linear in the variable i1 if (L1real) - isZero_real=true; - + isZero_real = true; + pow2_real++; - + if (L1imag) - isZero_imag=true; - + isZero_imag = true; + pow2_imag++; nActive--; @@ -410,33 +351,29 @@ scalar_t QuadraticForm::ExponentialSum() // set column i1 to zero M[i1] = zer; // set row i1 to zero - for (int j=0; j>i2) & one) ^ ((M[i2]>>i2) & one)); - bool L2imag = !!(((Limag>>i2) & one) ^ ((M[i2]>>i2) & one)); - Lreal&=~(one<> i2) & one) ^ ((M[i2] >> i2) & one)); + bool L2imag = !!(((Limag >> i2) & one) ^ ((M[i2] >> i2) & one)); + Lreal &= ~(one << i1); + Lreal &= ~(one << i2); + Limag &= ~(one << i1); + Limag &= ~(one << i2); // Extract rows i1 and i2 of M - uint_fast64_t m1=zer; - uint_fast64_t m2=zer; - for (int j=0; j>i1) & one)<>i2) & one)<> i1) & one) << j; + m2 ^= ((M[j] >> i2) & one) << j; } m1 ^= M[i1]; m2 ^= M[i2]; @@ -448,111 +385,80 @@ scalar_t QuadraticForm::ExponentialSum() M[i1] = zer; M[i2] = zer; // set rows i1,i2 to zero - for (int j=0; j> j) & one) - { - M[j]^=m1; - } - } + for (int j = 0; j < m; j++) + if ((m2 >> j) & one) + M[j] ^= m1; active &= ~(one << i1); active &= ~(one << i2); - nActive-=2; + nActive -= 2; } // int J0 = Q; // int J0 = inJ0[0]; - Q%=8; + Q %= 8; // Combine the real and the imaginary parts - if (isZero_imag) - { - amp.p = 2*pow2_real - 2; - amp.e = (4*sigma_real + Q) % 8; - } else if (isZero_real) - { - amp.p = 2*pow2_imag - 2; - amp.e = (2 + 4*sigma_imag + Q) % 8; + if (isZero_imag) { + amp.p = 2 * pow2_real - 2; + amp.e = (4 * sigma_real + Q) % 8; + } else if (isZero_real) { + amp.p = 2 * pow2_imag - 2; + amp.e = (2 + 4 * sigma_imag + Q) % 8; } else { // now both parts are nonzero - amp.p = 2*pow2_real - 1; + amp.p = 2 * pow2_real - 1; amp.eps = 1; if (sigma_real == 0) - { - if (sigma_imag == 0) - { - amp.e=(1 + Q) % 8; - } else { - amp.e=(7 + Q) % 8; - } - } else { if (sigma_imag == 0) - { - amp.e=(3 + Q) % 8; - } else { - amp.e=(5 + Q) % 8; - } - } + amp.e = (1 + Q) % 8; + else + amp.e = (7 + Q) % 8; + else if (sigma_imag == 0) + amp.e = (3 + Q) % 8; + else + amp.e = (5 + Q) % 8; } return amp; } // prints a binary form of x -void Print(uint_fast64_t x, unsigned n) -{ - for (unsigned i=0; i0); - } - std::cout< 0); + std::cout << std::endl; } -// prints a binary form of each element of A +// prints a binary form of each element of A // each element of A represents column of a binary matrix -void Print(std::vector A, unsigned n) -{ - for (unsigned row=0; row0); - } - std::cout< A, unsigned n) { + for (unsigned row = 0; row < n; row++) { + for (unsigned col = 0; col < n; col++) + std::cout << ((A[col] & (one << row)) > 0); + std::cout << std::endl; } } -} +} // namespace CHSimulator #endif diff --git a/src/simulators/extended_stabilizer/extended_stabilizer_state.hpp b/src/simulators/extended_stabilizer/extended_stabilizer_state.hpp index c8e648d3c8..f5a61ceed2 100644 --- a/src/simulators/extended_stabilizer/extended_stabilizer_state.hpp +++ b/src/simulators/extended_stabilizer/extended_stabilizer_state.hpp @@ -18,30 +18,35 @@ #include #include -#include "simulators/state.hpp" #include "framework/json.hpp" #include "framework/types.hpp" +#include "simulators/state.hpp" -#include "chlib/core.hpp" -#include "chlib/chstabilizer.hpp" #include "ch_runner.hpp" +#include "chlib/chstabilizer.hpp" +#include "chlib/core.hpp" #include "gates.hpp" -namespace AER{ +namespace AER { namespace ExtendedStabilizer { // OpSet of supported instructions const Operations::OpSet StateOpSet( - // Op types - {Operations::OpType::gate, Operations::OpType::measure, - Operations::OpType::reset, Operations::OpType::barrier, - Operations::OpType::roerror, Operations::OpType::bfunc, - Operations::OpType::qerror_loc, - Operations::OpType::save_statevec, - }, //Operations::OpType::save_expval, Operations::OpType::save_expval_var}, - // Gates - {"CX", "u0", "u1", "p", "cx", "cz", "swap", "id", "x", "y", "z", "h", - "s", "sdg", "sx", "sxdg", "t", "tdg", "ccx", "ccz", "delay", "pauli"}); + // Op types + { + Operations::OpType::gate, + Operations::OpType::measure, + Operations::OpType::reset, + Operations::OpType::barrier, + Operations::OpType::roerror, + Operations::OpType::bfunc, + Operations::OpType::qerror_loc, + Operations::OpType::save_statevec, + }, // Operations::OpType::save_expval, Operations::OpType::save_expval_var}, + // Gates + {"CX", "u0", "u1", "p", "cx", "cz", "swap", "id", + "x", "y", "z", "h", "s", "sdg", "sx", "sxdg", + "t", "tdg", "ccx", "ccz", "delay", "pauli"}); using chpauli_t = CHSimulator::pauli_t; using chstate_t = CHSimulator::Runner; @@ -50,85 +55,72 @@ using Gates = CHSimulator::Gates; uint_t zero = 0ULL; uint_t toff_branch_max = 7ULL; -enum class SamplingMethod { - metropolis, - resampled_metropolis, - norm_estimation -}; +enum class SamplingMethod { metropolis, resampled_metropolis, norm_estimation }; -class State: public QuantumState::State -{ +class State : public QuantumState::State { public: using BaseState = QuantumState::State; - + State() : BaseState(StateOpSet) {} virtual ~State() = default; - std::string name() const override {return "extended_stabilizer";} + std::string name() const override { return "extended_stabilizer"; } - //Apply a sequence of operations to the cicuit. For each operation, - //we loop over the terms in the decomposition in parallel + // Apply a sequence of operations to the cicuit. For each operation, + // we loop over the terms in the decomposition in parallel template void apply_ops(InputIterator first, InputIterator last, - ExperimentResult &result, - RngEngine &rng, - bool final_ops = false); + ExperimentResult &result, RngEngine &rng, + bool final_ops = false); // Apply a single operation // If the op is not in allowed_ops an exeption will be raised. - virtual void apply_op(const Operations::Op &op, - ExperimentResult &result, - RngEngine &rng, - bool final_op = false) override; + virtual void apply_op(const Operations::Op &op, ExperimentResult &result, + RngEngine &rng, bool final_op = false) override; void initialize_qreg(uint_t num_qubits) override; - size_t required_memory_mb(uint_t num_qubits, - const std::vector &ops) - const override; + size_t + required_memory_mb(uint_t num_qubits, + const std::vector &ops) const override; void set_config(const Config &config) override; - std::vector sample_measure(const reg_t& qubits, - uint_t shots, + std::vector sample_measure(const reg_t &qubits, uint_t shots, RngEngine &rng) override; protected: - - //Alongside the sample measure optimisaiton, we can parallelise - //circuit applicaiton over the states. This reduces the threading overhead - //as we only have to fork once per circuit. + // Alongside the sample measure optimisaiton, we can parallelise + // circuit applicaiton over the states. This reduces the threading overhead + // as we only have to fork once per circuit. template void apply_ops_parallel(InputIterator first, InputIterator last, - ExperimentResult &result, - RngEngine &rng); + ExperimentResult &result, RngEngine &rng); - //Small routine that eschews any parallelisation/decomposition and applies a stabilizer - //circuit to a single state. This is used to optimize a circuit with a large - //initial clifford fraction, or for running stabilizer circuits. + // Small routine that eschews any parallelisation/decomposition and applies a + // stabilizer circuit to a single state. This is used to optimize a circuit + // with a large initial clifford fraction, or for running stabilizer circuits. template void apply_stabilizer_circuit(InputIterator first, InputIterator last, - ExperimentResult &result, - RngEngine &rng); + ExperimentResult &result, RngEngine &rng); // Applies a sypported Gate operation to the state class. // If the input is not in allowed_gates an exeption will be raised. - // TODO: Investigate OMP synchronisation over stattes to remove these different versions - // One option would be tasks, but the memory overhead isn't clear + // TODO: Investigate OMP synchronisation over stattes to remove these + // different versions One option would be tasks, but the memory overhead isn't + // clear void apply_gate(const Operations::Op &op, RngEngine &rng); void apply_gate(const Operations::Op &op, RngEngine &rng, uint_t rank); // Apply a multi-qubit Pauli gate - void apply_pauli(const reg_t &qubits, const std::string& pauli, uint_t rank); + void apply_pauli(const reg_t &qubits, const std::string &pauli, uint_t rank); // Measure qubits and return a list of outcomes [q0, q1, ...] - // If a state subclass supports this function then "measure" + // If a state subclass supports this function then "measure" // should be contained in the set returned by the 'allowed_ops' // method. - void apply_measure(const reg_t &qubits, - const reg_t &cmemory, - const reg_t &cregister, - RngEngine &rng); + void apply_measure(const reg_t &qubits, const reg_t &cmemory, + const reg_t &cregister, RngEngine &rng); // Reset the specified qubits to the |0> state by measuring the // projectors Id+Z_{i} for each qubit i @@ -145,26 +137,24 @@ class State: public QuantumState::State ExperimentResult &result); // Compute and save the expval for the current simulator state - void apply_save_expval(const Operations::Op &op, - ExperimentResult &result, - RngEngine &rng); + void apply_save_expval(const Operations::Op &op, ExperimentResult &result, + RngEngine &rng); - // Helper function for computing expectation value - double expval_pauli(const reg_t &qubits, - const std::string& pauli, - RngEngine &rng); + // Helper function for computing expectation value + double expval_pauli(const reg_t &qubits, const std::string &pauli, + RngEngine &rng); // Helper function for computing expectation value virtual double expval_pauli(const reg_t &qubits, - const std::string& pauli) override; + const std::string &pauli) override; //----------------------------------------------------------------------- - //Parameters and methods specific to the Stabilizer Rank Decomposition + // Parameters and methods specific to the Stabilizer Rank Decomposition //----------------------------------------------------------------------- - //Allowed error in the stabilizer rank decomposition. - //The required number of states scales as \delta^{-2} - //for allowed error \delta + // Allowed error in the stabilizer rank decomposition. + // The required number of states scales as \delta^{-2} + // for allowed error \delta double approximation_error_ = 0.05; uint_t norm_estimation_samples_ = 100; @@ -176,7 +166,7 @@ class State: public QuantumState::State // output distribution uint_t metropolis_mixing_steps_ = 5000; - //Minimum number of states before we try to parallelise + // Minimum number of states before we try to parallelise uint_t omp_threshold_rank_ = 100; double snapshot_chop_threshold_ = 1e-10; @@ -192,14 +182,15 @@ class State: public QuantumState::State void compute_extent(const Operations::Op &op, double &xi) const; template - std::pair decomposition_parameters(InputIterator first, InputIterator last); - + std::pair decomposition_parameters(InputIterator first, + InputIterator last); + template - std::pair check_stabilizer_opt(InputIterator first, InputIterator last) const; + std::pair check_stabilizer_opt(InputIterator first, + InputIterator last) const; template bool check_measurement_opt(InputIterator first, InputIterator last) const; - }; //========================================================================= @@ -207,149 +198,135 @@ class State: public QuantumState::State //========================================================================= const stringmap_t State::gateset_({ - // Single qubit gates - {"delay", Gates::id}, // Delay gate - {"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)) - {"sx", Gates::sx}, // sqrt(X) gate - {"sxdg", Gates::sxdg}, // Inverse sqrt(X) gate - {"t", Gates::t}, // T-gate (sqrt(S)) - {"tdg", Gates::tdg}, // Conjguate-transpose of T gate - // Waltz Gates - {"u0", Gates::u0}, // idle gate in multiples of X90 - {"u1", Gates::u1}, // zero-X90 pulse waltz gate - {"p", Gates::u1}, // zero-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) - {"ccz", Gates::ccz}, // Constrolled-CZ gate (H3 Toff H3) - // Multi-qubit Pauli - {"pauli", Gates::pauli} // Multi-qubit Pauli gate + // Single qubit gates + {"delay", Gates::id}, // Delay gate + {"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)) + {"sx", Gates::sx}, // sqrt(X) gate + {"sxdg", Gates::sxdg}, // Inverse sqrt(X) gate + {"t", Gates::t}, // T-gate (sqrt(S)) + {"tdg", Gates::tdg}, // Conjguate-transpose of T gate + // Waltz Gates + {"u0", Gates::u0}, // idle gate in multiples of X90 + {"u1", Gates::u1}, // zero-X90 pulse waltz gate + {"p", Gates::u1}, // zero-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) + {"ccz", Gates::ccz}, // Constrolled-CZ gate (H3 Toff H3) + // Multi-qubit Pauli + {"pauli", Gates::pauli} // Multi-qubit Pauli gate }); //------------------------------------------------------------------------- // Implementation: Initialisation and Config //------------------------------------------------------------------------- -void State::initialize_qreg(uint_t num_qubits) -{ +void State::initialize_qreg(uint_t num_qubits) { BaseState::qreg_.initialize(num_qubits); BaseState::qreg_.initialize_omp(BaseState::threads_, omp_threshold_rank_); } -void State::set_config(const Config &config) -{ +void State::set_config(const Config &config) { // Set the error upper bound in the stabilizer rank approximation approximation_error_ = config.extended_stabilizer_approximation_error; // Set the number of samples used in the norm estimation routine if (config.extended_stabilizer_norm_estimation_default_samples.has_value()) - norm_estimation_samples_ = config.extended_stabilizer_norm_estimation_default_samples.value(); - // Set the desired number of repetitions of the norm estimation step. If not explicitly set, we - // compute a default basd on the approximation error - norm_estimation_repetitions_ = std::llrint(std::log2(1. / approximation_error_)); - norm_estimation_repetitions_ = config.extended_stabilizer_norm_estimation_repetitions; + norm_estimation_samples_ = + config.extended_stabilizer_norm_estimation_default_samples.value(); + // Set the desired number of repetitions of the norm estimation step. If not + // explicitly set, we compute a default basd on the approximation error + norm_estimation_repetitions_ = + std::llrint(std::log2(1. / approximation_error_)); + norm_estimation_repetitions_ = + config.extended_stabilizer_norm_estimation_repetitions; // Set the number of steps used in the metropolis sampler before we // consider the distribution as approximating the output metropolis_mixing_steps_ = config.extended_stabilizer_metropolis_mixing_time; - //Set the threshold of the decomposition before we use omp + // Set the threshold of the decomposition before we use omp omp_threshold_rank_ = config.extended_stabilizer_parallel_threshold; - //Set the truncation threshold for the probabilities snapshot. + // Set the truncation threshold for the probabilities snapshot. snapshot_chop_threshold_ = config.zero_threshold; - //Set the number of samples for the probabilities snapshot - probabilities_snapshot_samples_ = config.extended_stabilizer_probabilities_snapshot_samples; - //Set the measurement strategy + // Set the number of samples for the probabilities snapshot + probabilities_snapshot_samples_ = + config.extended_stabilizer_probabilities_snapshot_samples; + // Set the measurement strategy std::string sampling_method_str = "resampled_metropolis"; sampling_method_str = config.extended_stabilizer_sampling_method; if (sampling_method_str == "metropolis") { sampling_method_ = SamplingMethod::metropolis; - } - else if (sampling_method_str == "resampled_metropolis") - { + } else if (sampling_method_str == "resampled_metropolis") { sampling_method_ = SamplingMethod::resampled_metropolis; - } - else if (sampling_method_str == "norm_estimation") { + } else if (sampling_method_str == "norm_estimation") { sampling_method_ = SamplingMethod::norm_estimation; - } - else { + } else { throw std::runtime_error( - std::string("Unrecognised sampling method ") + sampling_method_str + - std::string("for the extended stabilizer simulator.") - ); + std::string("Unrecognised sampling method ") + sampling_method_str + + std::string("for the extended stabilizer simulator.")); } } template -std::pair State::decomposition_parameters(InputIterator first, InputIterator last) -{ - double xi=1.; +std::pair State::decomposition_parameters(InputIterator first, + InputIterator last) { + double xi = 1.; unsigned three_qubit_gate_count = 0; - for (auto op = first; op != last; op++) - { - if (op->type == Operations::OpType::gate) - { + for (auto op = first; op != last; op++) { + if (op->type == Operations::OpType::gate) { compute_extent(op, xi); auto it = CHSimulator::gate_types_.find(op->name); - if (it->second == CHSimulator::Gatetypes::non_clifford && op->qubits.size() == 3) - { //We count the number of 3 qubit gates for normalisation purposes + if (it->second == CHSimulator::Gatetypes::non_clifford && + op->qubits.size() == 3) { // We count the number of 3 qubit gates for + // normalisation purposes three_qubit_gate_count++; } } } - uint_t chi=1; - if (xi >1) - { + uint_t chi = 1; + if (xi > 1) { double err_scaling = std::pow(approximation_error_, -2); - chi = std::llrint(std::ceil(xi*err_scaling)); + chi = std::llrint(std::ceil(xi * err_scaling)); } return std::pair({chi, three_qubit_gate_count}); } template -std::pair State::check_stabilizer_opt(InputIterator first, InputIterator last) const -{ - for(auto op = first; op != last; op++) - { +std::pair State::check_stabilizer_opt(InputIterator first, + InputIterator last) const { + for (auto op = first; op != last; op++) { if (op->type != Operations::OpType::gate) - { continue; - } auto it = CHSimulator::gate_types_.find(op->name); - if(it == CHSimulator::gate_types_.end()) - { - throw std::invalid_argument("CHState::check_measurement_opt doesn't recognise a the operation \'"+ - op->name + "\'."); + if (it == CHSimulator::gate_types_.end()) { + throw std::invalid_argument("CHState::check_measurement_opt doesn't " + "recognise a the operation \'" + + op->name + "\'."); } - if(it->second == CHSimulator::Gatetypes::non_clifford) - { + if (it->second == CHSimulator::Gatetypes::non_clifford) return std::pair({false, op - first}); - } } return std::pair({true, 0}); } template -bool State::check_measurement_opt(InputIterator first, InputIterator last) const -{ - for (auto op = first; op != last; op++) - { +bool State::check_measurement_opt(InputIterator first, + InputIterator last) const { + for (auto op = first; op != last; op++) { if (op->conditional) - { return false; - } if (op->type == Operations::OpType::measure || op->type == Operations::OpType::bfunc || op->type == Operations::OpType::save_statevec || - op->type == Operations::OpType::save_expval) - { + op->type == Operations::OpType::save_expval) { return false; } } @@ -361,79 +338,73 @@ bool State::check_measurement_opt(InputIterator first, InputIterator last) const //------------------------------------------------------------------------- void State::apply_op(const Operations::Op &op, ExperimentResult &result, RngEngine &rng, bool final_op) { - apply_ops(&op, &op+1, result, rng, final_op); + apply_ops(&op, &op + 1, result, rng, final_op); } template -void State::apply_ops(InputIterator first, InputIterator last, ExperimentResult &result, - RngEngine &rng, bool final_ops) -{ +void State::apply_ops(InputIterator first, InputIterator last, + ExperimentResult &result, RngEngine &rng, + bool final_ops) { std::pair stabilizer_opts = check_stabilizer_opt(first, last); bool is_stabilizer = stabilizer_opts.first; - if(is_stabilizer) - { + if (is_stabilizer) { apply_stabilizer_circuit(first, last, result, rng); - } - else - { - //Split the circuit into stabilizer and non-stabilizer fractions + } else { + // Split the circuit into stabilizer and non-stabilizer fractions size_t first_non_clifford = stabilizer_opts.second; - if (first_non_clifford > 0) - { - //Apply the stabilizer circuit first. This optimisaiton avoids duplicating the application - //of the initial stabilizer circuit chi times. - apply_stabilizer_circuit(first, first+first_non_clifford, result, rng); + if (first_non_clifford > 0) { + // Apply the stabilizer circuit first. This optimisaiton avoids + // duplicating the application of the initial stabilizer circuit chi + // times. + apply_stabilizer_circuit(first, first + first_non_clifford, result, rng); } - auto it_nonstab_begin = first+first_non_clifford; + auto it_nonstab_begin = first + first_non_clifford; uint_t chi = compute_chi(it_nonstab_begin, last); double delta = std::pow(approximation_error_, -2); BaseState::qreg_.initialize_decomposition(chi, delta); - //Check for measurement optimisaitons + // Check for measurement optimisaitons bool measurement_opt = check_measurement_opt(first, last); - if(measurement_opt) - { + if (measurement_opt) { apply_ops_parallel(it_nonstab_begin, last, result, rng); - } - else - { - for (auto it = it_nonstab_begin; it != last; it++) - { + } else { + for (auto it = it_nonstab_begin; it != last; it++) { const auto op = *it; - if(BaseState::creg().check_conditional(op)) { + if (BaseState::creg().check_conditional(op)) { switch (op.type) { - case Operations::OpType::gate: - apply_gate(op, rng); - break; - case Operations::OpType::reset: - apply_reset(op.qubits, rng); - break; - case Operations::OpType::barrier: - case Operations::OpType::qerror_loc: - break; - case Operations::OpType::measure: - apply_measure(op.qubits, op.memory, op.registers, rng); - break; - case Operations::OpType::roerror: - BaseState::creg().apply_roerror(op, rng); - break; - case Operations::OpType::bfunc: - BaseState::creg().apply_bfunc(op); - break; - case Operations::OpType::save_statevec: - apply_save_statevector(op, result); - break; - // Disabled until can fix bug in expval - // case Operations::OpType::save_expval: - // case Operations::OpType::save_expval_var: - // apply_save_expval(op, result, rng); - // break; - default: - throw std::invalid_argument("CH::State::apply_ops does not support operations of the type \'" + - op.name + "\'."); - break; + case Operations::OpType::gate: + apply_gate(op, rng); + break; + case Operations::OpType::reset: + apply_reset(op.qubits, rng); + break; + case Operations::OpType::barrier: + case Operations::OpType::qerror_loc: + break; + case Operations::OpType::measure: + apply_measure(op.qubits, op.memory, op.registers, rng); + break; + case Operations::OpType::roerror: + BaseState::creg().apply_roerror(op, rng); + break; + case Operations::OpType::bfunc: + BaseState::creg().apply_bfunc(op); + break; + case Operations::OpType::save_statevec: + apply_save_statevector(op, result); + break; + // Disabled until can fix bug in expval + // case Operations::OpType::save_expval: + // case Operations::OpType::save_expval_var: + // apply_save_expval(op, result, rng); + // break; + default: + throw std::invalid_argument("CH::State::apply_ops does not support " + "operations of the type \'" + + op.name + "\'."); + break; } } } @@ -441,90 +412,69 @@ void State::apply_ops(InputIterator first, InputIterator last, ExperimentResult } } -std::vector State::sample_measure(const reg_t& qubits, - uint_t shots, - RngEngine &rng) -{ +std::vector State::sample_measure(const reg_t &qubits, uint_t shots, + RngEngine &rng) { std::vector output_samples; - if(BaseState::qreg_.get_num_states() == 1) - { + if (BaseState::qreg_.get_num_states() == 1) { output_samples = BaseState::qreg_.stabilizer_sampler(shots, rng); - } - else - { - if (sampling_method_ == SamplingMethod::metropolis) - { - output_samples = BaseState::qreg_.metropolis_estimation(metropolis_mixing_steps_, shots, rng); + } else if (sampling_method_ == SamplingMethod::metropolis) { + output_samples = BaseState::qreg_.metropolis_estimation( + metropolis_mixing_steps_, shots, rng); + } else if (sampling_method_ == SamplingMethod::resampled_metropolis) { + output_samples.reserve(shots); + for (uint_t i = 0; i < shots; i++) { + output_samples.push_back(BaseState::qreg_.metropolis_estimation( + metropolis_mixing_steps_, rng)); } - else if (sampling_method_ == SamplingMethod::resampled_metropolis) - { - output_samples.reserve(shots); - for (uint_t i=0; i all_samples; all_samples.reserve(shots); - for(uint_t sample: output_samples) - { + for (uint_t sample : output_samples) { reg_t sample_bits(qubits.size(), 0ULL); - for(size_t i=0; i> qubits[i]) & 1ULL) - { + for (size_t i = 0; i < qubits.size(); i++) + if ((sample >> qubits[i]) & 1ULL) sample_bits[i] = 1ULL; - } - } all_samples.push_back(sample_bits); } return all_samples; } - //------------------------------------------------------------------------- // Implemenation: Protected Methods //------------------------------------------------------------------------- -//Method with slighty optimized parallelisation for the case of a sample_measure circuit +// Method with slighty optimized parallelisation for the case of a +// sample_measure circuit template -void State::apply_ops_parallel(InputIterator first, InputIterator last, ExperimentResult &result, RngEngine &rng) -{ +void State::apply_ops_parallel(InputIterator first, InputIterator last, + ExperimentResult &result, RngEngine &rng) { const int_t NUM_STATES = BaseState::qreg_.get_num_states(); - #pragma omp parallel for if(BaseState::qreg_.check_omp_threshold() && BaseState::threads_>1) num_threads(BaseState::threads_) - for(int_t i=0; i < NUM_STATES; i++) - { - if(!BaseState::qreg_.check_eps(i)) - { +#pragma omp parallel for if (BaseState::qreg_.check_omp_threshold() && \ + BaseState::threads_ > 1) \ + num_threads(BaseState::threads_) + for (int_t i = 0; i < NUM_STATES; i++) { + if (!BaseState::qreg_.check_eps(i)) continue; - } - for(auto it = first; it != last; it++) - { - switch (it->type) - { - case Operations::OpType::gate: - apply_gate(*it, rng, i); - break; - case Operations::OpType::barrier: - case Operations::OpType::qerror_loc: - break; - default: - throw std::invalid_argument("CH::State::apply_ops_parallel does not support operations of the type \'" + - it->name + "\'."); - break; + for (auto it = first; it != last; it++) { + switch (it->type) { + case Operations::OpType::gate: + apply_gate(*it, rng, i); + break; + case Operations::OpType::barrier: + case Operations::OpType::qerror_loc: + break; + default: + throw std::invalid_argument("CH::State::apply_ops_parallel does not " + "support operations of the type \'" + + it->name + "\'."); + break; } } } @@ -532,100 +482,88 @@ void State::apply_ops_parallel(InputIterator first, InputIterator last, Experime template void State::apply_stabilizer_circuit(InputIterator first, InputIterator last, - ExperimentResult &result, RngEngine &rng) -{ - for (auto it = first; it != last; ++it) - { + ExperimentResult &result, RngEngine &rng) { + for (auto it = first; it != last; ++it) { const Operations::Op op = *it; - if(BaseState::creg().check_conditional(op)) { - switch (op.type) - { - case Operations::OpType::gate: - apply_gate(op, rng, 0); - break; - case Operations::OpType::reset: - apply_reset(op.qubits, rng); - break; - case Operations::OpType::barrier: - case Operations::OpType::qerror_loc: - break; - case Operations::OpType::measure: - apply_measure(op.qubits, op.memory, op.registers, rng); - break; - case Operations::OpType::roerror: - BaseState::creg().apply_roerror(op, rng); - break; - case Operations::OpType::bfunc: - BaseState::creg().apply_bfunc(op); - break; - case Operations::OpType::save_statevec: - apply_save_statevector(op, result); - break; - case Operations::OpType::save_expval: - case Operations::OpType::save_expval_var: - apply_save_expval(op, result, rng); - break; - default: - throw std::invalid_argument("CH::State::apply_stabilizer_circuit does not support operations of the type \'" + - op.name + "\'."); - break; + if (BaseState::creg().check_conditional(op)) { + switch (op.type) { + case Operations::OpType::gate: + apply_gate(op, rng, 0); + break; + case Operations::OpType::reset: + apply_reset(op.qubits, rng); + break; + case Operations::OpType::barrier: + case Operations::OpType::qerror_loc: + break; + case Operations::OpType::measure: + apply_measure(op.qubits, op.memory, op.registers, rng); + break; + case Operations::OpType::roerror: + BaseState::creg().apply_roerror(op, rng); + break; + case Operations::OpType::bfunc: + BaseState::creg().apply_bfunc(op); + break; + case Operations::OpType::save_statevec: + apply_save_statevector(op, result); + break; + case Operations::OpType::save_expval: + case Operations::OpType::save_expval_var: + apply_save_expval(op, result, rng); + break; + default: + throw std::invalid_argument("CH::State::apply_stabilizer_circuit does " + "not support operations of the type \'" + + op.name + "\'."); + break; } } } } -void State::apply_measure(const reg_t &qubits, const reg_t &cmemory, const reg_t &cregister, RngEngine &rng) -{ +void State::apply_measure(const reg_t &qubits, const reg_t &cmemory, + const reg_t &cregister, RngEngine &rng) { uint_t out_string; // Flag if the Pauli projector is applied already as part of the sampling bool do_projector_correction = true; // Prepare an output register for the qubits we are measurig reg_t outcome(qubits.size(), 0ULL); - if(BaseState::qreg_.get_num_states() == 1) - { - //For a single state, we use the efficient sampler defined in Sec IV.A ofarxiv:1808.00128 + if (BaseState::qreg_.get_num_states() == 1) { + // For a single state, we use the efficient sampler defined in Sec IV.A + // ofarxiv:1808.00128 out_string = BaseState::qreg_.stabilizer_sampler(rng); - } - else - { - if (sampling_method_ == SamplingMethod::norm_estimation) - { - do_projector_correction = false; - //Run the norm estimation routine - out_string = BaseState::qreg_.ne_single_sample( - norm_estimation_samples_, norm_estimation_repetitions_, false, qubits, rng - ); - } - else - { - // We use the metropolis algorithm to sample an output string non-destructively - // This is a single measure step so we do the same for metropolis or resampled_metropolis - out_string = BaseState::qreg_.metropolis_estimation(metropolis_mixing_steps_, rng); - } - } - if (do_projector_correction) - { - //We prepare the Pauli projector corresponding to the measurement result - std::vectorpaulis(qubits.size(), chpauli_t()); - for (uint_t i=0; i paulis(qubits.size(), chpauli_t()); + for (uint_t i = 0; i < qubits.size(); i++) { + // Create a Pauli projector onto 1+-Z_{i} on qubit i paulis[i].Z = (1ULL << qubits[i]); - if ((out_string >> qubits[i]) & 1ULL) - { - //Additionally, store the output bit for this qubit + if ((out_string >> qubits[i]) & 1ULL) { + // Additionally, store the output bit for this qubit paulis[i].e = 2; } } - //Project the decomposition onto the measurement outcome + // Project the decomposition onto the measurement outcome BaseState::qreg_.apply_pauli_projector(paulis); } - for (uint_t i=0; i> qubits[i]) & 1ULL) - { - //Additionally, store the output bit for this qubit + for (uint_t i = 0; i < qubits.size(); i++) { + // Create a Pauli projector onto 1+-Z_{i} on qubit i + if ((out_string >> qubits[i]) & 1ULL) { + // Additionally, store the output bit for this qubit outcome[i] = 1ULL; } } @@ -633,143 +571,132 @@ void State::apply_measure(const reg_t &qubits, const reg_t &cmemory, const reg_t BaseState::creg().store_measure(outcome, cmemory, cregister); } -void State::apply_reset(const reg_t &qubits, AER::RngEngine &rng) -{ +void State::apply_reset(const reg_t &qubits, AER::RngEngine &rng) { uint_t measure_string; const int_t NUM_STATES = BaseState::qreg_.get_num_states(); - if(BaseState::qreg_.get_num_states() == 1) - { + if (BaseState::qreg_.get_num_states() == 1) { measure_string = BaseState::qreg_.stabilizer_sampler(rng); - } - else - { - measure_string = BaseState::qreg_.metropolis_estimation(metropolis_mixing_steps_, rng); + } else { + measure_string = + BaseState::qreg_.metropolis_estimation(metropolis_mixing_steps_, rng); } std::vector paulis(qubits.size(), chpauli_t()); - for(size_t i=0; i 1 && BaseState::qreg_.check_omp_threshold()) num_threads(BaseState::threads_) - for(int_t i=0; i< NUM_STATES; i++) - { - for (auto qubit: qubits) - { - if ((measure_string>>qubit) & 1ULL) - { +#pragma omp parallel for if (BaseState::threads_ > 1 && \ + BaseState::qreg_.check_omp_threshold()) \ + num_threads(BaseState::threads_) + for (int_t i = 0; i < NUM_STATES; i++) { + for (auto qubit : qubits) + if ((measure_string >> qubit) & 1ULL) BaseState::qreg_.apply_x(qubit, i); - } - } } } -void State::apply_gate(const Operations::Op &op, RngEngine &rng) -{ +void State::apply_gate(const Operations::Op &op, RngEngine &rng) { const int_t NUM_STATES = BaseState::qreg_.get_num_states(); - #pragma omp parallel for if (BaseState::threads_ > 1 && BaseState::qreg_.check_omp_threshold()) num_threads(BaseState::threads_) - for(int_t i=0; i < NUM_STATES; i++) - { - if(BaseState::qreg_.check_eps(i)) - { +#pragma omp parallel for if (BaseState::threads_ > 1 && \ + BaseState::qreg_.check_omp_threshold()) \ + num_threads(BaseState::threads_) + for (int_t i = 0; i < NUM_STATES; i++) + if (BaseState::qreg_.check_eps(i)) apply_gate(op, rng, i); - } - } } -void State::apply_gate(const Operations::Op &op, RngEngine &rng, uint_t rank) -{ +void State::apply_gate(const Operations::Op &op, RngEngine &rng, uint_t rank) { auto it = gateset_.find(op.name); - if (it == gateset_.end()) - { - throw std::invalid_argument("CH::State: Invalid gate operation \'" - +op.name + "\'."); - } - switch(it->second) - { - case Gates::x: - BaseState::qreg_.apply_x(op.qubits[0], rank); - break; - case Gates::y: - BaseState::qreg_.apply_y(op.qubits[0], rank); - break; - case Gates::z: - BaseState::qreg_.apply_z(op.qubits[0], rank); - break; - case Gates::s: - BaseState::qreg_.apply_s(op.qubits[0], rank); - break; - case Gates::sdg: - BaseState::qreg_.apply_sdag(op.qubits[0], rank); - break; - case Gates::h: - BaseState::qreg_.apply_h(op.qubits[0], rank); - break; - case Gates::sx: - BaseState::add_global_phase(M_PI / 4.); - BaseState::qreg_.apply_sx(op.qubits[0], rank); - break; - case Gates::sxdg: - BaseState::add_global_phase(-M_PI / 4.); - BaseState::qreg_.apply_sxdg(op.qubits[0], rank); - break; - case Gates::cx: - BaseState::qreg_.apply_cx(op.qubits[0], op.qubits[1], rank); - break; - case Gates::cz: - BaseState::qreg_.apply_cz(op.qubits[0], op.qubits[1], rank); - break; - case Gates::swap: - BaseState::qreg_.apply_swap(op.qubits[0], op.qubits[1], rank); - break; - case Gates::t: - BaseState::qreg_.apply_t(op.qubits[0], rng.rand(), rank); - break; - case Gates::tdg: - BaseState::qreg_.apply_tdag(op.qubits[0], rng.rand(), rank); - break; - case Gates::ccx: - BaseState::qreg_.apply_ccx(op.qubits[0], op.qubits[1], op.qubits[2], rng.rand_int(zero, toff_branch_max), rank); - break; - case Gates::ccz: - BaseState::qreg_.apply_ccz(op.qubits[0], op.qubits[1], op.qubits[2], rng.rand_int(zero, toff_branch_max), rank); - break; - case Gates::u1: - BaseState::qreg_.apply_u1(op.qubits[0], op.params[0], rng.rand(), rank); - break; - case Gates::pauli: - apply_pauli(op.qubits, op.string_params[0], rank); - break; - default: //u0 or Identity - break; + if (it == gateset_.end()) { + throw std::invalid_argument("CH::State: Invalid gate operation \'" + + op.name + "\'."); + } + switch (it->second) { + case Gates::x: + BaseState::qreg_.apply_x(op.qubits[0], rank); + break; + case Gates::y: + BaseState::qreg_.apply_y(op.qubits[0], rank); + break; + case Gates::z: + BaseState::qreg_.apply_z(op.qubits[0], rank); + break; + case Gates::s: + BaseState::qreg_.apply_s(op.qubits[0], rank); + break; + case Gates::sdg: + BaseState::qreg_.apply_sdag(op.qubits[0], rank); + break; + case Gates::h: + BaseState::qreg_.apply_h(op.qubits[0], rank); + break; + case Gates::sx: + BaseState::add_global_phase(M_PI / 4.); + BaseState::qreg_.apply_sx(op.qubits[0], rank); + break; + case Gates::sxdg: + BaseState::add_global_phase(-M_PI / 4.); + BaseState::qreg_.apply_sxdg(op.qubits[0], rank); + break; + case Gates::cx: + BaseState::qreg_.apply_cx(op.qubits[0], op.qubits[1], rank); + break; + case Gates::cz: + BaseState::qreg_.apply_cz(op.qubits[0], op.qubits[1], rank); + break; + case Gates::swap: + BaseState::qreg_.apply_swap(op.qubits[0], op.qubits[1], rank); + break; + case Gates::t: + BaseState::qreg_.apply_t(op.qubits[0], rng.rand(), rank); + break; + case Gates::tdg: + BaseState::qreg_.apply_tdag(op.qubits[0], rng.rand(), rank); + break; + case Gates::ccx: + BaseState::qreg_.apply_ccx(op.qubits[0], op.qubits[1], op.qubits[2], + rng.rand_int(zero, toff_branch_max), rank); + break; + case Gates::ccz: + BaseState::qreg_.apply_ccz(op.qubits[0], op.qubits[1], op.qubits[2], + rng.rand_int(zero, toff_branch_max), rank); + break; + case Gates::u1: + BaseState::qreg_.apply_u1(op.qubits[0], op.params[0], rng.rand(), rank); + break; + case Gates::pauli: + apply_pauli(op.qubits, op.string_params[0], rank); + break; + default: // u0 or Identity + break; } } -void State::apply_pauli(const reg_t &qubits, const std::string& pauli, uint_t rank) { +void State::apply_pauli(const reg_t &qubits, const std::string &pauli, + uint_t rank) { const auto size = qubits.size(); for (size_t i = 0; i < qubits.size(); ++i) { const auto qubit = qubits[size - 1 - i]; switch (pauli[i]) { - case 'I': - break; - case 'X': - BaseState::qreg_.apply_x(qubit, rank); - break; - case 'Y': - BaseState::qreg_.apply_y(qubit, rank); - break; - case 'Z': - BaseState::qreg_.apply_z(qubit, rank); - break; - default: - throw std::invalid_argument("invalid Pauli \'" + std::to_string(pauli[i]) + "\'."); + case 'I': + break; + case 'X': + BaseState::qreg_.apply_x(qubit, rank); + break; + case 'Y': + BaseState::qreg_.apply_y(qubit, rank); + break; + case 'Z': + BaseState::qreg_.apply_z(qubit, rank); + break; + default: + throw std::invalid_argument("invalid Pauli \'" + + std::to_string(pauli[i]) + "\'."); } } } @@ -782,18 +709,14 @@ void State::apply_save_statevector(const Operations::Op &op, " Only the full statevector can be saved."); } auto statevec = BaseState::qreg_.statevector(); - if (BaseState::has_global_phase_) { + if (BaseState::has_global_phase_) statevec *= BaseState::global_phase_; - } - result.save_data_pershot( - creg(), op.string_params[0], - std::move(statevec), - op.type, op.save_type); + result.save_data_pershot(creg(), op.string_params[0], std::move(statevec), + op.type, op.save_type); } void State::apply_save_expval(const Operations::Op &op, - ExperimentResult &result, - RngEngine& rng) { + ExperimentResult &result, RngEngine &rng) { // Check empty edge case if (op.expval_params.empty()) { throw std::invalid_argument( @@ -809,135 +732,123 @@ void State::apply_save_expval(const Operations::Op &op, // param is tuple (pauli, coeff, sq_coeff) const auto val = expval_pauli(op.qubits, std::get<0>(param), rng); expval += std::get<1>(param) * val; - if (variance) { + if (variance) sq_expval += std::get<2>(param) * val; - } } if (variance) { std::vector expval_var(2); - expval_var[0] = expval; // mean - expval_var[1] = sq_expval - expval * expval; // variance - result.save_data_average(creg(), op.string_params[0], expval_var, op.type, op.save_type); + expval_var[0] = expval; // mean + expval_var[1] = sq_expval - expval * expval; // variance + result.save_data_average(creg(), op.string_params[0], expval_var, op.type, + op.save_type); } else { - result.save_data_average(creg(), op.string_params[0], expval, op.type, op.save_type); + result.save_data_average(creg(), op.string_params[0], expval, op.type, + op.save_type); } } -double State::expval_pauli(const reg_t &qubits, - const std::string& pauli, +double State::expval_pauli(const reg_t &qubits, const std::string &pauli, RngEngine &rng) { - // Compute expval components - auto state_cpy = BaseState::qreg_; - auto phi_norm = state_cpy.norm_estimation(norm_estimation_samples_, norm_estimation_repetitions_, rng); - std::vectorpaulis(1, chpauli_t()); - for (uint_t pos = 0; pos < qubits.size(); ++pos) { - switch (pauli[pauli.size() - 1 - pos]) { - case 'I': - break; - case 'X': - paulis[0].X += (1ULL << qubits[pos]); - break; - case 'Y': - paulis[0].X += (1ULL << qubits[pos]); - paulis[0].Z += (1ULL << qubits[pos]); - break; - case 'Z': - paulis[0].Z += (1ULL << qubits[pos]); - break; - default: { - std::stringstream msg; - msg << "QubitVectorState::invalid Pauli string \'" << pauli[pos] - << "\'."; - throw std::invalid_argument(msg.str()); - } - } + // Compute expval components + auto state_cpy = BaseState::qreg_; + auto phi_norm = state_cpy.norm_estimation(norm_estimation_samples_, + norm_estimation_repetitions_, rng); + std::vector paulis(1, chpauli_t()); + for (uint_t pos = 0; pos < qubits.size(); ++pos) { + switch (pauli[pauli.size() - 1 - pos]) { + case 'I': + break; + case 'X': + paulis[0].X += (1ULL << qubits[pos]); + break; + case 'Y': + paulis[0].X += (1ULL << qubits[pos]); + paulis[0].Z += (1ULL << qubits[pos]); + break; + case 'Z': + paulis[0].Z += (1ULL << qubits[pos]); + break; + default: { + std::stringstream msg; + msg << "QubitVectorState::invalid Pauli string \'" << pauli[pos] << "\'."; + throw std::invalid_argument(msg.str()); } - auto g_norm = state_cpy.norm_estimation(norm_estimation_samples_, norm_estimation_repetitions_, paulis, rng); - return (2*g_norm - phi_norm); + } + } + auto g_norm = state_cpy.norm_estimation( + norm_estimation_samples_, norm_estimation_repetitions_, paulis, rng); + return (2 * g_norm - phi_norm); } -double State::expval_pauli(const reg_t &qubits, - const std::string& pauli) { - // empty implementation of base class virtual method - // since in the extended stabilizer, expval relies on RNG - return 0; +double State::expval_pauli(const reg_t &qubits, const std::string &pauli) { + // empty implementation of base class virtual method + // since in the extended stabilizer, expval relies on RNG + return 0; } //------------------------------------------------------------------------- // Implementation: Utility //------------------------------------------------------------------------- -inline void to_json(json_t &js, cvector_t vec) -{ +inline void to_json(json_t &js, cvector_t vec) { js = json_t(); - for (uint_t j=0; j < vec.size(); j++) { + for (uint_t j = 0; j < vec.size(); j++) js.push_back(vec[j]); - } } template -uint_t State::compute_chi(InputIterator first, InputIterator last) const -{ +uint_t State::compute_chi(InputIterator first, InputIterator last) const { double xi = 1; for (auto op = first; op != last; op++) - { compute_extent(*op, xi); - } double err_scaling = std::pow(approximation_error_, -2); - return std::llrint(std::ceil(xi*err_scaling)); + return std::llrint(std::ceil(xi * err_scaling)); } -void State::compute_extent(const Operations::Op &op, double &xi) const -{ - if(op.type == Operations::OpType::gate) - { +void State::compute_extent(const Operations::Op &op, double &xi) const { + if (op.type == Operations::OpType::gate) { auto it = gateset_.find(op.name); - if (it == gateset_.end()) - { - throw std::invalid_argument("CH::State: Invalid gate operation \'" - +op.name + "\'."); + if (it == gateset_.end()) { + throw std::invalid_argument("CH::State: Invalid gate operation \'" + + op.name + "\'."); } - switch (it->second) - { - case Gates::t: - xi *= CHSimulator::t_extent; - break; - case Gates::tdg: - xi *= CHSimulator::t_extent; - break; - case Gates::ccx: - xi *= CHSimulator::ccx_extent; - break; - case Gates::ccz: - xi *= CHSimulator::ccx_extent; - break; - case Gates::u1: - xi *= CHSimulator::u1_extent(std::real(op.params[0])); - break; - default: - break; + switch (it->second) { + case Gates::t: + xi *= CHSimulator::t_extent; + break; + case Gates::tdg: + xi *= CHSimulator::t_extent; + break; + case Gates::ccx: + xi *= CHSimulator::ccx_extent; + break; + case Gates::ccz: + xi *= CHSimulator::ccx_extent; + break; + case Gates::u1: + xi *= CHSimulator::u1_extent(std::real(op.params[0])); + break; + default: + break; } } } size_t State::required_memory_mb(uint_t num_qubits, - const std::vector &ops) const -{ + const std::vector &ops) const { size_t required_chi = compute_chi(ops.cbegin(), ops.cend()); // 5 vectors of num_qubits*8byte words // Plus 2*CHSimulator::scalar_t which has 3 4 byte words // Plus 2*CHSimulator::pauli_t which has 2 8 byte words and one 4 byte word; - double mb_per_state = 5e-5*num_qubits;// - size_t required_mb = std::llrint(std::ceil(mb_per_state*required_chi)); + double mb_per_state = 5e-5 * num_qubits; // + size_t required_mb = std::llrint(std::ceil(mb_per_state * required_chi)); if (sampling_method_ == SamplingMethod::norm_estimation) - { required_mb *= 2; - } return required_mb; } -} -} +} // namespace ExtendedStabilizer +} // namespace AER #endif diff --git a/src/simulators/extended_stabilizer/gates.hpp b/src/simulators/extended_stabilizer/gates.hpp index 29004f2621..98be8717c0 100644 --- a/src/simulators/extended_stabilizer/gates.hpp +++ b/src/simulators/extended_stabilizer/gates.hpp @@ -24,24 +24,35 @@ #include "framework/operations.hpp" #include "framework/types.hpp" -namespace CHSimulator -{ - using uint_t = uint_fast64_t; - using complex_t = std::complex; - - enum class Gates { - u0, u1, id, x, y, z, h, s, sdg, sx, sxdg, t, tdg, - cx, cz, swap, - ccx, ccz, pauli - }; - - enum class Gatetypes { - pauli, - clifford, - non_clifford - }; - - const AER::stringmap_t gate_types_ = { +namespace CHSimulator { +using uint_t = uint_fast64_t; +using complex_t = std::complex; + +enum class Gates { + u0, + u1, + id, + x, + y, + z, + h, + s, + sdg, + sx, + sxdg, + t, + tdg, + cx, + cz, + swap, + ccx, + ccz, + pauli +}; + +enum class Gatetypes { pauli, clifford, non_clifford }; + +const AER::stringmap_t gate_types_ = { // One-qubit gates {"u0", Gatetypes::pauli}, // Pauli-Identity gate {"delay", Gatetypes::pauli}, // Pauli-Identity gate @@ -59,51 +70,50 @@ namespace CHSimulator {"u1", Gatetypes::non_clifford}, // zero-X90 pulse waltz gate {"p", Gatetypes::non_clifford}, // zero-X90 pulse waltz gate // Two-qubit gates - {"CX", Gatetypes::clifford}, // Controlled-X gate (CNOT) - {"cx", Gatetypes::clifford}, // Controlled-X gate (CNOT) - {"cz", Gatetypes::clifford}, // Controlled-Z gate - {"swap", Gatetypes::clifford}, // SWAP gate + {"CX", Gatetypes::clifford}, // Controlled-X gate (CNOT) + {"cx", Gatetypes::clifford}, // Controlled-X gate (CNOT) + {"cz", Gatetypes::clifford}, // Controlled-Z gate + {"swap", Gatetypes::clifford}, // SWAP gate // Three-qubit gates {"ccx", Gatetypes::non_clifford}, // Controlled-CX gate (Toffoli) {"ccz", Gatetypes::non_clifford}, // Controlled-CZ gate (H3 Toff H3) // N-qubit gates - {"pauli", Gatetypes::pauli}, // N-qubit Pauli gate - }; - - using sample_branch_t = std::pair; - - const double root2 = std::sqrt(2); - const double root1_2 = 1./root2; - const complex_t pi_over_8_phase(0., M_PI/8); - const complex_t omega(root1_2, root1_2); - const complex_t omega_star(root1_2, -1*root1_2); - const complex_t root_omega = std::exp(pi_over_8_phase); - const complex_t root_omega_star = std::conj(root_omega); - - const double tan_pi_over_8 = std::tan(M_PI/8.); - - //Base class for sampling over non-Clifford gates in the Sum Over Cliffords routine. - struct Sample { - public: - Sample() = default; - virtual ~Sample() = default; - Sample(const Sample& other) : branches(other.branches) {}; - std::vector branches; - - virtual sample_branch_t sample(double r) const = 0; - }; - - //Functor class that defines how to sample branches over a U1 operation - //Used for caching each U1 gate angle we encounter in the circuit -struct U1Sample : public Sample -{ + {"pauli", Gatetypes::pauli}, // N-qubit Pauli gate +}; + +using sample_branch_t = std::pair; + +const double root2 = std::sqrt(2); +const double root1_2 = 1. / root2; +const complex_t pi_over_8_phase(0., M_PI / 8); +const complex_t omega(root1_2, root1_2); +const complex_t omega_star(root1_2, -1 * root1_2); +const complex_t root_omega = std::exp(pi_over_8_phase); +const complex_t root_omega_star = std::conj(root_omega); + +const double tan_pi_over_8 = std::tan(M_PI / 8.); + +// Base class for sampling over non-Clifford gates in the Sum Over Cliffords +// routine. +struct Sample { +public: + Sample() = default; + virtual ~Sample() = default; + Sample(const Sample &other) : branches(other.branches){}; + std::vector branches; + + virtual sample_branch_t sample(double r) const = 0; +}; + +// Functor class that defines how to sample branches over a U1 operation +// Used for caching each U1 gate angle we encounter in the circuit +struct U1Sample : public Sample { double p_threshold; U1Sample() = default; U1Sample(double lambda); - U1Sample(const U1Sample &other) : Sample(other) - { + U1Sample(const U1Sample &other) : Sample(other) { p_threshold = other.p_threshold; } @@ -112,135 +122,94 @@ struct U1Sample : public Sample sample_branch_t sample(double r) const override; }; -U1Sample::U1Sample(double lambda) -{ +U1Sample::U1Sample(double lambda) { // Shift parameter into +- 2 Pi - uint_t shift_factor = std::floor(std::abs(lambda)/(2*M_PI)); - if (shift_factor != 0) - { - if(lambda < 0) - { - lambda += shift_factor*(2*M_PI); - } + uint_t shift_factor = std::floor(std::abs(lambda) / (2 * M_PI)); + if (shift_factor != 0) { + if (lambda < 0) + lambda += shift_factor * (2 * M_PI); else - { - lambda -= shift_factor*(2*M_PI); - } + lambda -= shift_factor * (2 * M_PI); } // Shift parameter into +- Pi if (lambda > M_PI) - { - lambda -= 2*M_PI; - } - else if (lambda < -1*M_PI) - { - lambda += 2*M_PI; - } + lambda -= 2 * M_PI; + else if (lambda < -1 * M_PI) + lambda += 2 * M_PI; // Compute the coefficients double angle = std::abs(lambda); - bool s_z_quadrant = (angle > M_PI/2); - if(s_z_quadrant) - { - angle = angle - M_PI/2; - } + bool s_z_quadrant = (angle > M_PI / 2); + if (s_z_quadrant) + angle = angle - M_PI / 2; angle /= 2; - complex_t coeff_0 = std::cos(angle)-std::sin(angle); - complex_t coeff_1 = root2*std::sin(angle); + complex_t coeff_0 = std::cos(angle) - std::sin(angle); + complex_t coeff_1 = root2 * std::sin(angle); complex_t phase_0, phase_1; std::array gates; - if(lambda < 0) - { + if (lambda < 0) { coeff_0 *= root_omega_star; coeff_1 = coeff_1 * root_omega; - if(s_z_quadrant) - { + if (s_z_quadrant) { gates[0] = Gates::sdg; gates[1] = Gates::z; - } - else - { + } else { gates[0] = Gates::id; gates[1] = Gates::sdg; } - } - else - { + } else { coeff_0 *= root_omega; coeff_1 = coeff_1 * root_omega_star; - if(s_z_quadrant) - { + if (s_z_quadrant) { gates[0] = Gates::s; gates[1] = Gates::z; - } - else - { + } else { gates[0] = Gates::id; gates[1] = Gates::s; } } phase_0 = std::polar(1.0, std::arg(coeff_0)); phase_1 = std::polar(1.0, std::arg(coeff_1)); - branches = - { - sample_branch_t(phase_0, gates[0]), - sample_branch_t(phase_1, gates[1]) - }; - p_threshold = std::abs(coeff_0) / (std::abs(coeff_0)+std::abs(coeff_1)); + branches = {sample_branch_t(phase_0, gates[0]), + sample_branch_t(phase_1, gates[1])}; + p_threshold = std::abs(coeff_0) / (std::abs(coeff_0) + std::abs(coeff_1)); } -sample_branch_t U1Sample::sample(double r) const -{ - if (r M_PI) - { - lambda -= 2*M_PI; - } - else if (lambda < -1*M_PI) - { - lambda += 2*M_PI; - } - double angle = std::abs(lambda); - bool s_z_quadrant = (angle > M_PI/2); - if(s_z_quadrant) - { - angle = angle - M_PI/2; - } - angle /= 2; - return std::pow(std::cos(angle) + tan_pi_over_8*std::sin(angle), 2.); +// Stabilizer extent for different non-Clifford unitaries +// Special case of Eq. 28 in arXiv:1808.00128 +const double t_extent = std::pow(1. / (std::cos(M_PI / 8.)), 2); +// Equation 32 in arXiv: 1808.00128 +const double ccx_extent = 16. / 9.; +const double ccx_coeff = 1. / 6.; +// General result for z rotations, Eq. 28 in arXiv 1809.00128 +inline double u1_extent(double lambda) { + // Shift parameter into +- 2 Pi + uint_t shift_factor = std::floor(std::abs(lambda) / (2 * M_PI)); + if (shift_factor != 0) { + if (lambda < 0) + lambda += shift_factor * (2 * M_PI); + else + lambda -= shift_factor * (2 * M_PI); } - + // Shift parameter into +- Pi + if (lambda > M_PI) + lambda -= 2 * M_PI; + else if (lambda < -1 * M_PI) + lambda += 2 * M_PI; + double angle = std::abs(lambda); + bool s_z_quadrant = (angle > M_PI / 2); + if (s_z_quadrant) + angle = angle - M_PI / 2; + angle /= 2; + return std::pow(std::cos(angle) + tan_pi_over_8 * std::sin(angle), 2.); } +} // namespace CHSimulator + #endif diff --git a/src/simulators/stabilizer/binary_vector.hpp b/src/simulators/stabilizer/binary_vector.hpp index 8520084e7c..8e587da19e 100644 --- a/src/simulators/stabilizer/binary_vector.hpp +++ b/src/simulators/stabilizer/binary_vector.hpp @@ -12,15 +12,14 @@ * that they have been altered from the originals. */ - #ifndef _binary_vector_hpp_ #define _binary_vector_hpp_ #include +#include #include #include #include -#include namespace AER { namespace BV { @@ -37,10 +36,11 @@ class BinaryVector { const static size_t BLOCK_BITS = 6; const static size_t BLOCK_MASK = (1ull << BLOCK_BITS) - 1; - BinaryVector() : m_length(0) {}; + BinaryVector() : m_length(0){}; explicit BinaryVector(uint64_t length) - : m_length(length), m_data((length + BLOCK_SIZE - 1) >> BLOCK_BITS, ZERO_){}; + : m_length(length), + m_data((length + BLOCK_SIZE - 1) >> BLOCK_BITS, ZERO_){}; BinaryVector(std::vector mdata) : m_length(mdata.size() << BLOCK_BITS), m_data(mdata){}; @@ -65,7 +65,9 @@ class BinaryVector { uint64_t getLength() const { return m_length; }; - void makeZero() { m_data.assign((m_length + BLOCK_SIZE - 1) >> BLOCK_BITS, ZERO_); } + void makeZero() { + m_data.assign((m_length + BLOCK_SIZE - 1) >> BLOCK_BITS, ZERO_); + } bool isZero() const; @@ -73,27 +75,12 @@ class BinaryVector { bool isSame(const BinaryVector &rhs, bool pad) const; std::vector nonzeroIndices() const; - const std::vector& getData() const - { - return m_data; - } + const std::vector &getData() const { return m_data; } - size_t blockSize(void) - { - return BLOCK_SIZE; - } - size_t blockLength(void) const - { - return m_data.size(); - } - uint64_t& operator()(const uint64_t pos) - { - return m_data[pos]; - } - uint64_t operator()(const uint64_t pos) const - { - return m_data[pos]; - } + size_t blockSize(void) { return BLOCK_SIZE; } + size_t blockLength(void) const { return m_data.size(); } + uint64_t &operator()(const uint64_t pos) { return m_data[pos]; } + uint64_t operator()(const uint64_t pos) const { return m_data[pos]; } protected: uint64_t m_length; @@ -102,7 +89,6 @@ class BinaryVector { static const uint64_t ONE_; }; - /******************************************************************************* * * Related Functions @@ -113,7 +99,6 @@ inline bool operator==(const BinaryVector &lhs, const BinaryVector &rhs) { return lhs.isSame(rhs, true); } - inline int64_t gauss_eliminate(std::vector &M, const int64_t start_col = 0) // returns the rank of M. @@ -136,17 +121,14 @@ inline int64_t gauss_eliminate(std::vector &M, M[r] += M[i]; } } - if (i >= rank) { + if (i >= rank) M[i].swap(M[rank - 1]); - } } return rank; } - -inline std::vector string_to_bignum(std::string val, - uint64_t blockSize, - uint64_t base) { +inline std::vector +string_to_bignum(std::string val, uint64_t blockSize, uint64_t base) { std::vector ret; if (blockSize * log2(base) > 64) { throw std::runtime_error( @@ -163,7 +145,6 @@ inline std::vector string_to_bignum(std::string val, return ret; } - inline std::vector string_to_bignum(std::string val) { std::string type = val.substr(0, 2); if (type == "0b" || type == "0B") @@ -179,7 +160,6 @@ inline std::vector string_to_bignum(std::string val) { } } - /******************************************************************************* * * BinaryVector Class Methods @@ -189,22 +169,17 @@ inline std::vector string_to_bignum(std::string val) { const uint64_t BinaryVector::ZERO_ = 0ULL; const uint64_t BinaryVector::ONE_ = 1ULL; -BinaryVector::BinaryVector(std::string val) -{ +BinaryVector::BinaryVector(std::string val) { m_data = string_to_bignum(val); m_length = m_data.size() << BLOCK_BITS; } - -void BinaryVector::setLength(uint64_t length) -{ +void BinaryVector::setLength(uint64_t length) { m_length = length; m_data.assign((length + BLOCK_SIZE - 1) >> BLOCK_BITS, ZERO_); } - -void BinaryVector::setValue(bool value, uint64_t pos) -{ +void BinaryVector::setValue(bool value, uint64_t pos) { auto q = pos >> BLOCK_BITS; auto r = pos & BLOCK_MASK; if (value) @@ -213,15 +188,12 @@ void BinaryVector::setValue(bool value, uint64_t pos) m_data[q] &= ~(ONE_ << r); } - -void BinaryVector::flipAt(const uint64_t pos) -{ +void BinaryVector::flipAt(const uint64_t pos) { auto q = pos >> BLOCK_BITS; auto r = pos & BLOCK_MASK; m_data[q] ^= (ONE_ << r); } - BinaryVector &BinaryVector::operator+=(const BinaryVector &rhs) { const auto size = m_data.size(); for (size_t i = 0; i < size; i++) @@ -229,15 +201,12 @@ BinaryVector &BinaryVector::operator+=(const BinaryVector &rhs) { return (*this); } - -bool BinaryVector::operator[](const uint64_t pos) const -{ +bool BinaryVector::operator[](const uint64_t pos) const { auto q = pos >> BLOCK_BITS; auto r = pos & BLOCK_MASK; return ((m_data[q] & (ONE_ << r)) != 0); } - void BinaryVector::swap(BinaryVector &rhs) { uint64_t tmp; tmp = rhs.m_length; @@ -247,7 +216,6 @@ void BinaryVector::swap(BinaryVector &rhs) { m_data.swap(rhs.m_data); } - bool BinaryVector::isZero() const { const size_t size = m_data.size(); for (size_t i = 0; i < size; i++) @@ -256,19 +224,16 @@ bool BinaryVector::isZero() const { return true; } - bool BinaryVector::isSame(const BinaryVector &rhs) const { if (m_length != rhs.m_length) return false; const size_t size = m_data.size(); - for (size_t q = 0; q < size; q++) { + for (size_t q = 0; q < size; q++) if (m_data[q] != rhs.m_data[q]) return false; - } return true; } - bool BinaryVector::isSame(const BinaryVector &rhs, bool pad) const { if (!pad) return isSame(rhs); @@ -292,7 +257,6 @@ bool BinaryVector::isSame(const BinaryVector &rhs, bool pad) const { return true; } - std::vector BinaryVector::nonzeroIndices() const { std::vector result; size_t i = 0; @@ -305,9 +269,8 @@ std::vector BinaryVector::nonzeroIndices() const { auto m = m_data[i]; size_t r = 0; while (r < BLOCK_SIZE) { - while ((m & (ONE_ << r)) == 0) { + while ((m & (ONE_ << r)) == 0) r++; - } if (r >= BLOCK_SIZE) break; result.push_back(static_cast((i * BLOCK_SIZE) + r)); @@ -318,9 +281,8 @@ std::vector BinaryVector::nonzeroIndices() const { return result; } - //------------------------------------------------------------------------------ } // end namespace BV -} // AER +} // namespace AER //------------------------------------------------------------------------------ #endif \ No newline at end of file diff --git a/src/simulators/stabilizer/clifford.hpp b/src/simulators/stabilizer/clifford.hpp index b2e30ccedd..00bbaf5b96 100644 --- a/src/simulators/stabilizer/clifford.hpp +++ b/src/simulators/stabilizer/clifford.hpp @@ -18,8 +18,8 @@ #include "framework/types.hpp" #include "framework/utils.hpp" -#include "pauli.hpp" #include "framework/json_parser.hpp" +#include "pauli.hpp" #include @@ -48,27 +48,28 @@ class Clifford { //----------------------------------------------------------------------- // Utility functions //----------------------------------------------------------------------- - //initialize + // initialize void initialize(const uint64_t nqubit); // Get number of qubits of the Clifford table - uint64_t num_qubits() const {return num_qubits_;} + uint64_t num_qubits() const { return num_qubits_; } // Return true if the number of qubits is 0 - bool empty() const {return (num_qubits_ == 0);} + bool empty() const { return (num_qubits_ == 0); } // Return JSON serialization of QubitVector; json_t json() const; // Access stabilizer table -// Pauli::Pauli &operator[](uint64_t j) {return table_[j];} -// const Pauli::Pauli& operator[](uint64_t j) const {return table_[j];} + // Pauli::Pauli &operator[](uint64_t j) {return table_[j];} + // const Pauli::Pauli& operator[](uint64_t j) const {return + // table_[j];} - //set stabilizer - void set_destabilizer(const int i, const Pauli::Pauli& P); - void set_stabilizer(const int i, const Pauli::Pauli& P); + // set stabilizer + void set_destabilizer(const int i, const Pauli::Pauli &P); + void set_stabilizer(const int i, const Pauli::Pauli &P); - //set phase + // set phase void set_destabilizer_phases(const int i, const bool p); void set_stabilizer_phases(const int i, const bool p); @@ -101,17 +102,16 @@ class Clifford { // Measurement //----------------------------------------------------------------------- - // If we perform a single qubit Z measurement, + // If we perform a single qubit Z measurement, // will the outcome be random or deterministic. - bool is_deterministic_outcome(const uint64_t& qubit) const; + bool is_deterministic_outcome(const uint64_t &qubit) const; // Return the outcome (0 or 1) of a single qubit Z measurement, and // update the stabilizer to the conditional (post measurement) state if // the outcome was random. bool measure_and_update(const uint64_t qubit, const uint64_t randint); - double expval_pauli(const reg_t &qubits, - const std::string& pauli); + double expval_pauli(const reg_t &qubits, const std::string &pauli); //----------------------------------------------------------------------- // Configuration settings @@ -121,23 +121,22 @@ class Clifford { void set_json_chop_threshold(double threshold); // Set the threshold for chopping values to 0 in JSON - double get_json_chop_threshold() {return json_chop_threshold_;} + double get_json_chop_threshold() { return json_chop_threshold_; } // Set the maximum number of OpenMP thread for operations. void set_omp_threads(int n); // Get the maximum number of OpenMP thread for operations. - uint64_t get_omp_threads() {return omp_threads_;} + uint64_t get_omp_threads() { return omp_threads_; } // Set the qubit threshold for activating OpenMP. // If self.qubits() > threshold OpenMP will be activated. void set_omp_threshold(int n); // Get the qubit threshold for activating OpenMP. - uint64_t get_omp_threshold() {return omp_threshold_;} + uint64_t get_omp_threshold() { return omp_threshold_; } protected: - //----------------------------------------------------------------------- // Protected data members //----------------------------------------------------------------------- @@ -151,10 +150,11 @@ class Clifford { // Config settings //----------------------------------------------------------------------- - uint64_t omp_threads_ = 1; // Disable multithreading by default - uint64_t omp_threshold_ = 1000; // Qubit threshold for multithreading when enabled - double json_chop_threshold_ = 0; // Threshold for chopping small values - // in JSON serialization + uint64_t omp_threads_ = 1; // Disable multithreading by default + uint64_t omp_threshold_ = + 1000; // Qubit threshold for multithreading when enabled + double json_chop_threshold_ = 0; // Threshold for chopping small values + // in JSON serialization //----------------------------------------------------------------------- // Helper functions @@ -197,30 +197,27 @@ void Clifford::set_omp_threshold(int n) { // Constructors & Destructor //------------------------------------------------------------------------------ -Clifford::Clifford(uint64_t nq) : num_qubits_(nq) -{ - initialize(nq); -} +Clifford::Clifford(uint64_t nq) : num_qubits_(nq) { initialize(nq); } -void Clifford::initialize(uint64_t nq) -{ +void Clifford::initialize(uint64_t nq) { num_qubits_ = nq; destabilizer_table_.resize(nq); stabilizer_table_.resize(nq); int nid = omp_get_num_threads(); - auto init_func = [this, nq](AER::int_t i) - { + auto init_func = [this, nq](AER::int_t i) { destabilizer_table_[i].X.setLength(nq); destabilizer_table_[i].Z.setLength(nq); - destabilizer_table_[i].X.setValue(1,i); + destabilizer_table_[i].X.setValue(1, i); stabilizer_table_[i].X.setLength(nq); stabilizer_table_[i].Z.setLength(nq); - stabilizer_table_[i].Z.setValue(1,i); + stabilizer_table_[i].Z.setValue(1, i); }; - AER::Utils::apply_omp_parallel_for((num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0, nq, init_func, omp_threads_); + AER::Utils::apply_omp_parallel_for( + (num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0, nq, + init_func, omp_threads_); // Add phases destabilizer_phases_.setLength(nq); @@ -231,33 +228,42 @@ void Clifford::initialize(uint64_t nq) // Apply Clifford gates //------------------------------------------------------------------------------ -void Clifford::append_cx(const uint64_t qcon, const uint64_t qtar) -{ +void Clifford::append_cx(const uint64_t qcon, const uint64_t qtar) { const uint64_t mask = (~0ull); int nid = omp_get_num_threads(); - auto cx_func = [this, qtar, qcon, mask](AER::int_t i) - { - destabilizer_phases_(i) = destabilizer_phases_(i) ^ (destabilizer_table_[qcon].X(i) & destabilizer_table_[qtar].Z(i) & - (destabilizer_table_[qtar].X(i) ^ destabilizer_table_[qcon].Z(i) ^ mask)); - stabilizer_phases_(i) = stabilizer_phases_(i) ^ (stabilizer_table_[qcon].X(i) & stabilizer_table_[qtar].Z(i) & - (stabilizer_table_[qtar].X(i) ^ stabilizer_table_[qcon].Z(i) ^ mask)); - - destabilizer_table_[qtar].X(i) = destabilizer_table_[qtar].X(i) ^ destabilizer_table_[qcon].X(i); - destabilizer_table_[qcon].Z(i) = destabilizer_table_[qtar].Z(i) ^ destabilizer_table_[qcon].Z(i); - stabilizer_table_[qtar].X(i) = stabilizer_table_[qtar].X(i) ^ stabilizer_table_[qcon].X(i); - stabilizer_table_[qcon].Z(i) = stabilizer_table_[qtar].Z(i) ^ stabilizer_table_[qcon].Z(i); + auto cx_func = [this, qtar, qcon, mask](AER::int_t i) { + destabilizer_phases_(i) = + destabilizer_phases_(i) ^ + (destabilizer_table_[qcon].X(i) & destabilizer_table_[qtar].Z(i) & + (destabilizer_table_[qtar].X(i) ^ destabilizer_table_[qcon].Z(i) ^ + mask)); + stabilizer_phases_(i) = + stabilizer_phases_(i) ^ + (stabilizer_table_[qcon].X(i) & stabilizer_table_[qtar].Z(i) & + (stabilizer_table_[qtar].X(i) ^ stabilizer_table_[qcon].Z(i) ^ mask)); + + destabilizer_table_[qtar].X(i) = + destabilizer_table_[qtar].X(i) ^ destabilizer_table_[qcon].X(i); + destabilizer_table_[qcon].Z(i) = + destabilizer_table_[qtar].Z(i) ^ destabilizer_table_[qcon].Z(i); + stabilizer_table_[qtar].X(i) = + stabilizer_table_[qtar].X(i) ^ stabilizer_table_[qcon].X(i); + stabilizer_table_[qcon].Z(i) = + stabilizer_table_[qtar].Z(i) ^ stabilizer_table_[qcon].Z(i); }; - AER::Utils::apply_omp_parallel_for((num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0, destabilizer_phases_.blockLength(), cx_func, omp_threads_); + AER::Utils::apply_omp_parallel_for( + (num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0, + destabilizer_phases_.blockLength(), cx_func, omp_threads_); } -void Clifford::append_h(const uint64_t qubit) -{ +void Clifford::append_h(const uint64_t qubit) { int nid = omp_get_num_threads(); - auto h_func = [this, qubit](AER::int_t i) - { - destabilizer_phases_(i) ^= (destabilizer_table_[qubit].X(i) & destabilizer_table_[qubit].Z(i)); - stabilizer_phases_(i) ^= (stabilizer_table_[qubit].X(i) & stabilizer_table_[qubit].Z(i)); + auto h_func = [this, qubit](AER::int_t i) { + destabilizer_phases_(i) ^= + (destabilizer_table_[qubit].X(i) & destabilizer_table_[qubit].Z(i)); + stabilizer_phases_(i) ^= + (stabilizer_table_[qubit].X(i) & stabilizer_table_[qubit].Z(i)); // exchange X and Z uint64_t t = destabilizer_table_[qubit].X(i); destabilizer_table_[qubit].X(i) = destabilizer_table_[qubit].Z(i); @@ -266,115 +272,115 @@ void Clifford::append_h(const uint64_t qubit) stabilizer_table_[qubit].X(i) = stabilizer_table_[qubit].Z(i); stabilizer_table_[qubit].Z(i) = t; }; - AER::Utils::apply_omp_parallel_for((num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0, destabilizer_phases_.blockLength(), h_func, omp_threads_); + AER::Utils::apply_omp_parallel_for( + (num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0, + destabilizer_phases_.blockLength(), h_func, omp_threads_); } -void Clifford::append_s(const uint64_t qubit) -{ +void Clifford::append_s(const uint64_t qubit) { int nid = omp_get_num_threads(); - auto s_func = [this, qubit](AER::int_t i) - { - destabilizer_phases_(i) ^= (destabilizer_table_[qubit].X(i) & destabilizer_table_[qubit].Z(i)); + auto s_func = [this, qubit](AER::int_t i) { + destabilizer_phases_(i) ^= + (destabilizer_table_[qubit].X(i) & destabilizer_table_[qubit].Z(i)); destabilizer_table_[qubit].Z(i) ^= destabilizer_table_[qubit].X(i); - stabilizer_phases_(i) ^= (stabilizer_table_[qubit].X(i) & stabilizer_table_[qubit].Z(i)); + stabilizer_phases_(i) ^= + (stabilizer_table_[qubit].X(i) & stabilizer_table_[qubit].Z(i)); stabilizer_table_[qubit].Z(i) ^= stabilizer_table_[qubit].X(i); }; - AER::Utils::apply_omp_parallel_for((num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0, destabilizer_phases_.blockLength(), s_func, omp_threads_); + AER::Utils::apply_omp_parallel_for( + (num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0, + destabilizer_phases_.blockLength(), s_func, omp_threads_); } -void Clifford::append_x(const uint64_t qubit) -{ +void Clifford::append_x(const uint64_t qubit) { int nid = omp_get_num_threads(); - auto x_func = [this, qubit](AER::int_t i) - { + auto x_func = [this, qubit](AER::int_t i) { destabilizer_phases_(i) ^= destabilizer_table_[qubit].Z(i); stabilizer_phases_(i) ^= stabilizer_table_[qubit].Z(i); }; - AER::Utils::apply_omp_parallel_for((num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0, destabilizer_phases_.blockLength(), x_func, omp_threads_); + AER::Utils::apply_omp_parallel_for( + (num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0, + destabilizer_phases_.blockLength(), x_func, omp_threads_); } -void Clifford::append_y(const uint64_t qubit) -{ +void Clifford::append_y(const uint64_t qubit) { int nid = omp_get_num_threads(); - auto y_func = [this, qubit](AER::int_t i) - { - destabilizer_phases_(i) ^= (destabilizer_table_[qubit].Z(i) ^ destabilizer_table_[qubit].X(i)); - stabilizer_phases_(i) ^= (stabilizer_table_[qubit].Z(i) ^ stabilizer_table_[qubit].X(i)); + auto y_func = [this, qubit](AER::int_t i) { + destabilizer_phases_(i) ^= + (destabilizer_table_[qubit].Z(i) ^ destabilizer_table_[qubit].X(i)); + stabilizer_phases_(i) ^= + (stabilizer_table_[qubit].Z(i) ^ stabilizer_table_[qubit].X(i)); }; - AER::Utils::apply_omp_parallel_for((num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0, destabilizer_phases_.blockLength(), y_func, omp_threads_); + AER::Utils::apply_omp_parallel_for( + (num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0, + destabilizer_phases_.blockLength(), y_func, omp_threads_); } -void Clifford::append_z(const uint64_t qubit) -{ +void Clifford::append_z(const uint64_t qubit) { int nid = omp_get_num_threads(); - auto z_func = [this, qubit](AER::int_t i) - { + auto z_func = [this, qubit](AER::int_t i) { destabilizer_phases_(i) ^= destabilizer_table_[qubit].X(i); stabilizer_phases_(i) ^= stabilizer_table_[qubit].X(i); }; - AER::Utils::apply_omp_parallel_for((num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0, destabilizer_phases_.blockLength(), z_func, omp_threads_); + AER::Utils::apply_omp_parallel_for( + (num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0, + destabilizer_phases_.blockLength(), z_func, omp_threads_); } //------------------------------------------------------------------------------ // Utility //------------------------------------------------------------------------------ -std::pair Clifford::z_anticommuting(const uint64_t qubit) const -{ +std::pair +Clifford::z_anticommuting(const uint64_t qubit) const { for (uint_t i = 0; i < stabilizer_table_[qubit].X.blockLength(); i++) { - if (stabilizer_table_[qubit].X(i) != 0){ + if (stabilizer_table_[qubit].X(i) != 0) { uint_t p = i << stabilizer_table_[qubit].X.BLOCK_BITS; - for (uint_t j = 0; j < stabilizer_table_[qubit].X.BLOCK_SIZE; j++) { - if (stabilizer_table_[qubit].X[p+j]) - return std::make_pair(true, p+j); - } + for (uint_t j = 0; j < stabilizer_table_[qubit].X.BLOCK_SIZE; j++) + if (stabilizer_table_[qubit].X[p + j]) + return std::make_pair(true, p + j); } } return std::make_pair(false, 0); } - -std::pair Clifford::x_anticommuting(const uint64_t qubit) const -{ +std::pair +Clifford::x_anticommuting(const uint64_t qubit) const { for (uint_t i = 0; i < stabilizer_table_[qubit].Z.blockLength(); i++) { - if (stabilizer_table_[qubit].Z(i) != 0){ + if (stabilizer_table_[qubit].Z(i) != 0) { uint_t p = i << stabilizer_table_[qubit].Z.BLOCK_BITS; - for (uint_t j = 0; j < stabilizer_table_[qubit].Z.BLOCK_SIZE; j++) { - if (stabilizer_table_[qubit].Z[p+j]) - return std::make_pair(true, p+j); - } + for (uint_t j = 0; j < stabilizer_table_[qubit].Z.BLOCK_SIZE; j++) + if (stabilizer_table_[qubit].Z[p + j]) + return std::make_pair(true, p + j); } } return std::make_pair(false, 0); } -void Clifford::set_destabilizer(const int idx, const Pauli::Pauli& P) -{ - for (int64_t i = 0; i < static_cast( num_qubits_); i++){ +void Clifford::set_destabilizer(const int idx, + const Pauli::Pauli &P) { + for (int64_t i = 0; i < static_cast(num_qubits_); i++) { destabilizer_table_[i].X.setValue(P.X[i], idx); destabilizer_table_[i].Z.setValue(P.Z[i], idx); } } -void Clifford::set_stabilizer(const int idx, const Pauli::Pauli& P) -{ - for (int64_t i = 0; i < static_cast( num_qubits_); i++){ +void Clifford::set_stabilizer(const int idx, + const Pauli::Pauli &P) { + for (int64_t i = 0; i < static_cast(num_qubits_); i++) { stabilizer_table_[i].X.setValue(P.X[i], idx); stabilizer_table_[i].Z.setValue(P.Z[i], idx); } } -void Clifford::set_destabilizer_phases(const int i, const bool p) -{ +void Clifford::set_destabilizer_phases(const int i, const bool p) { destabilizer_phases_.setValue(p, i); } -void Clifford::set_stabilizer_phases(const int i, const bool p) -{ +void Clifford::set_stabilizer_phases(const int i, const bool p) { stabilizer_phases_.setValue(p, i); } -void Clifford::apply_set_stabilizer(const Clifford &clifford) -{ +void Clifford::apply_set_stabilizer(const Clifford &clifford) { destabilizer_table_ = clifford.destabilizer_table_; stabilizer_table_ = clifford.stabilizer_table_; destabilizer_phases_ = clifford.destabilizer_phases_; @@ -385,15 +391,15 @@ void Clifford::apply_set_stabilizer(const Clifford &clifford) // Measurement //------------------------------------------------------------------------------ -bool Clifford::is_deterministic_outcome(const uint64_t& qubit) const { +bool Clifford::is_deterministic_outcome(const uint64_t &qubit) const { // Clifford state measurements only have three probabilities: // (p0, p1) = (0.5, 0.5), (1, 0), or (0, 1) // The random case happens if there is a row anti-commuting with Z[qubit] return !z_anticommuting(qubit).first; } -bool Clifford::measure_and_update(const uint64_t qubit, const uint64_t randint) -{ +bool Clifford::measure_and_update(const uint64_t qubit, + const uint64_t randint) { // Clifford state measurements only have three probabilities: // (p0, p1) = (0.5, 0.5), (1, 0), or (0, 1) // The random case happens if there is a row anti-commuting with Z[qubit] @@ -407,25 +413,25 @@ bool Clifford::measure_and_update(const uint64_t qubit, const uint64_t randint) uint64_t rS = 0ull - (uint64_t)stabilizer_phases_[row]; - auto measure_non_determinisitic_func = [this, rS, row, qubit](AER::int_t i) - { + auto measure_non_determinisitic_func = [this, rS, row, + qubit](AER::int_t i) { uint64_t row_mask = ~0ull; - if((row >> destabilizer_phases_.BLOCK_BITS) == i) + if ((row >> destabilizer_phases_.BLOCK_BITS) == i) row_mask ^= (1ull << (row & destabilizer_phases_.BLOCK_MASK)); uint64_t d_mask = row_mask & destabilizer_table_[qubit].X(i); uint64_t s_mask = row_mask & stabilizer_table_[qubit].X(i); - if(d_mask != 0 || s_mask != 0){ - //calculating exponents by 2-bits integer * 64-qubits at once + if (d_mask != 0 || s_mask != 0) { + // calculating exponents by 2-bits integer * 64-qubits at once uint64_t d0 = 0, d1 = 0; uint64_t s0 = 0, s1 = 0; - for (size_t q = 0; q < num_qubits_; q++){ - uint64_t t0,t1; + for (size_t q = 0; q < num_qubits_; q++) { + uint64_t t0, t1; uint64_t rX = 0ull - (uint64_t)stabilizer_table_[q].X[row]; uint64_t rZ = 0ull - (uint64_t)stabilizer_table_[q].Z[row]; - //destabilizer + // destabilizer t0 = destabilizer_table_[q].X(i) & rZ; t1 = destabilizer_table_[q].Z(i) ^ rX; @@ -444,7 +450,7 @@ bool Clifford::measure_and_update(const uint64_t qubit, const uint64_t randint) destabilizer_table_[q].X(i) ^= (d_mask & rX); destabilizer_table_[q].Z(i) ^= (d_mask & rZ); - //stabilizer + // stabilizer t0 = stabilizer_table_[q].X(i) & rZ; t1 = stabilizer_table_[q].Z(i) ^ rX; @@ -464,40 +470,45 @@ bool Clifford::measure_and_update(const uint64_t qubit, const uint64_t randint) stabilizer_table_[q].Z(i) ^= (s_mask & rZ); } d1 ^= (rS ^ destabilizer_phases_(i)); - destabilizer_phases_(i) = (destabilizer_phases_(i) & (~d_mask)) | (d1 & d_mask); + destabilizer_phases_(i) = + (destabilizer_phases_(i) & (~d_mask)) | (d1 & d_mask); s1 ^= (rS ^ stabilizer_phases_(i)); - stabilizer_phases_(i) = (stabilizer_phases_(i) & (~s_mask)) | (s1 & s_mask); + stabilizer_phases_(i) = + (stabilizer_phases_(i) & (~s_mask)) | (s1 & s_mask); } }; - AER::Utils::apply_omp_parallel_for((num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0, destabilizer_phases_.blockLength(), measure_non_determinisitic_func, omp_threads_); + AER::Utils::apply_omp_parallel_for( + (num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0, + destabilizer_phases_.blockLength(), measure_non_determinisitic_func, + omp_threads_); // Update state - auto measure_update_func = [this, row](AER::int_t q) - { + auto measure_update_func = [this, row](AER::int_t q) { destabilizer_table_[q].X.setValue(stabilizer_table_[q].X[row], row); destabilizer_table_[q].Z.setValue(stabilizer_table_[q].Z[row], row); stabilizer_table_[q].X.setValue(0, row); stabilizer_table_[q].Z.setValue(0, row); }; - AER::Utils::apply_omp_parallel_for((num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0, num_qubits_, measure_update_func, omp_threads_); + AER::Utils::apply_omp_parallel_for( + (num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0, + num_qubits_, measure_update_func, omp_threads_); destabilizer_phases_.setValue(stabilizer_phases_[row], row); stabilizer_table_[qubit].Z.setValue(1, row); stabilizer_phases_.setValue(outcome, row); return outcome; - } - else { + } else { // Deterministic outcome bool outcome = false; Pauli::Pauli accum(num_qubits_); uint64_t blocks = destabilizer_phases_.blockLength(); - if(blocks < 2){ + if (blocks < 2) { for (uint64_t i = 0; i < num_qubits_; i++) { if (destabilizer_table_[qubit].X[i]) { bool b0 = false, b1 = false; - for (size_t q = 0; q < num_qubits_; q++){ - bool t0,t1,add; + for (size_t q = 0; q < num_qubits_; q++) { + bool t0, t1, add; bool accumX = accum.X[q]; bool accumZ = accum.Z[q]; @@ -521,84 +532,88 @@ bool Clifford::measure_and_update(const uint64_t qubit, const uint64_t randint) } b1 ^= (stabilizer_phases_[i] ^ outcome); - if(b0){ + if (b0) throw std::runtime_error("Clifford: rowsum error"); - } outcome = b1; } } - } - else{ + } else { uint64_t blockSize = destabilizer_phases_.blockSize(); - //loop for cache blocking + // loop for cache blocking for (uint64_t ii = 0; ii < blocks; ii++) { uint64_t destabilizer_mask = destabilizer_table_[qubit].X(ii); - if(destabilizer_mask == 0) + if (destabilizer_mask == 0) continue; uint64_t exponent_l = 0; uint64_t exponent_lc = 0; uint64_t exponent_h = 0; - auto measure_determinisitic_func = [this, &accum, &exponent_l, &exponent_lc, &exponent_h, blocks, blockSize, destabilizer_mask, ii](AER::int_t qq) - { - uint64_t qs = qq*blockSize; - uint64_t qe = qs + blockSize; - if(qe > num_qubits_) - qe = num_qubits_; - - uint64_t local_exponent_l = 0; - uint64_t local_exponent_h = 0; - - for (uint64_t q = qs; q < qe; q++){ - uint64_t sX = stabilizer_table_[q].X(ii); - uint64_t sZ = stabilizer_table_[q].Z(ii); - - //set accum for this block - uint64_t accumX = destabilizer_mask & sX; - uint64_t accumZ = destabilizer_mask & sZ; - for(int b=1;b> (blockSize - 1)), q); - accum.Z.setValue((accumZ >> (blockSize - 1)), q); - - accumX ^= sX; - accumZ ^= sZ; - - //exponents for this block - uint64_t t0,t1; - - t0 = accumX & sZ; - t1 = accumZ ^ sX; - - local_exponent_h ^= (t0 & local_exponent_l); - local_exponent_l ^= t0; - local_exponent_h ^= (t0 & t1); - - t0 = sX & accumZ; - t1 = sZ ^ accumX; - t1 ^= t0; - - local_exponent_h ^= (t0 & local_exponent_l); - local_exponent_l ^= t0; - local_exponent_h ^= (t0 & t1); - } + auto measure_determinisitic_func = + [this, &accum, &exponent_l, &exponent_lc, &exponent_h, blocks, + blockSize, destabilizer_mask, ii](AER::int_t qq) { + uint64_t qs = qq * blockSize; + uint64_t qe = qs + blockSize; + if (qe > num_qubits_) + qe = num_qubits_; + + uint64_t local_exponent_l = 0; + uint64_t local_exponent_h = 0; + + for (uint64_t q = qs; q < qe; q++) { + uint64_t sX = stabilizer_table_[q].X(ii); + uint64_t sZ = stabilizer_table_[q].Z(ii); + + // set accum for this block + uint64_t accumX = destabilizer_mask & sX; + uint64_t accumZ = destabilizer_mask & sZ; + for (int b = 1; b < blockSize; b *= 2) { + accumX ^= (accumX << b); + accumZ ^= (accumZ << b); + } + accumX ^= (0ull - (uint64_t)accum.X[q]); + accumZ ^= (0ull - (uint64_t)accum.Z[q]); + accum.X.setValue((accumX >> (blockSize - 1)), q); + accum.Z.setValue((accumZ >> (blockSize - 1)), q); + + accumX ^= sX; + accumZ ^= sZ; + + // exponents for this block + uint64_t t0, t1; + + t0 = accumX & sZ; + t1 = accumZ ^ sX; + + local_exponent_h ^= (t0 & local_exponent_l); + local_exponent_l ^= t0; + local_exponent_h ^= (t0 & t1); + + t0 = sX & accumZ; + t1 = sZ ^ accumX; + t1 ^= t0; + + local_exponent_h ^= (t0 & local_exponent_l); + local_exponent_l ^= t0; + local_exponent_h ^= (t0 & t1); + } #pragma omp atomic - exponent_lc |= local_exponent_l; + exponent_lc |= local_exponent_l; #pragma omp atomic - exponent_l ^= local_exponent_l; + exponent_l ^= local_exponent_l; #pragma omp atomic - exponent_h ^= local_exponent_h; - }; - AER::Utils::apply_omp_parallel_for((num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0, blocks, measure_determinisitic_func, omp_threads_); + exponent_h ^= local_exponent_h; + }; + AER::Utils::apply_omp_parallel_for( + (num_qubits_ > omp_threshold_ && omp_threads_ > 1 && nid == 1), 0, + blocks, measure_determinisitic_func, omp_threads_); - exponent_h ^= (exponent_lc ^ exponent_l); //if exponent_l is 0 and any of local_exponent_l is 1, then flip exponent_h + exponent_h ^= + (exponent_lc ^ + exponent_l); // if exponent_l is 0 and any of local_exponent_l is + // 1, then flip exponent_h exponent_h ^= (stabilizer_phases_(ii) & destabilizer_mask); outcome ^= ((AER::Utils::popcount(exponent_h) & 1) != 0); @@ -608,62 +623,56 @@ bool Clifford::measure_and_update(const uint64_t qubit, const uint64_t randint) } } -double Clifford::expval_pauli(const reg_t &qubits, - const std::string& pauli) -{ +double Clifford::expval_pauli(const reg_t &qubits, const std::string &pauli) { // Construct Pauli on N-qubits Pauli::Pauli P(num_qubits_); uint_t phase = 0; - for (size_t i = 0; i < qubits.size(); ++i) { + for (size_t i = 0; i < qubits.size(); ++i) switch (pauli[pauli.size() - 1 - i]) { - case 'X': - P.X.set1(qubits[i]); - break; - case 'Y': - P.X.set1(qubits[i]); - P.Z.set1(qubits[i]); - phase += 1; - break; - case 'Z': - P.Z.set1(qubits[i]); - break; - default: - break; + case 'X': + P.X.set1(qubits[i]); + break; + case 'Y': + P.X.set1(qubits[i]); + P.Z.set1(qubits[i]); + phase += 1; + break; + case 'Z': + P.Z.set1(qubits[i]); + break; + default: + break; }; - } - // Check if there is a stabilizer that anti-commutes with an odd number of qubits - // If so expectation value is 0 + // Check if there is a stabilizer that anti-commutes with an odd number of + // qubits If so expectation value is 0 for (size_t i = 0; i < num_qubits_; i++) { size_t num_anti = 0; - for (const auto& qubit : qubits) { - if (P.Z[qubit] & stabilizer_table_[qubit].X[i]) { - num_anti++; - } - if (P.X[qubit] & stabilizer_table_[qubit].Z[i]) { - num_anti++; - } + for (const auto &qubit : qubits) { + if (P.Z[qubit] & stabilizer_table_[qubit].X[i]) + num_anti++; + if (P.X[qubit] & stabilizer_table_[qubit].Z[i]) + num_anti++; } - if(num_anti % 2 == 1) + if (num_anti % 2 == 1) return 0.0; } // Otherwise P is (-1)^a prod_j S_j^b_j for Clifford stabilizers // If P anti-commutes with D_j then b_j = 1. // Multiply P by stabilizers with anti-commuting destabilizers - auto PZ = P.Z; // Make a copy of P.Z + auto PZ = P.Z; // Make a copy of P.Z for (size_t i = 0; i < num_qubits_; i++) { // Check if destabilizer anti-commutes size_t num_anti = 0; - for (const auto& qubit : qubits) { - if (P.Z[qubit] & destabilizer_table_[qubit].X[i]) { - num_anti++; - } - if (P.X[qubit] & destabilizer_table_[qubit].Z[i]) { - num_anti++; - } + for (const auto &qubit : qubits) { + if (P.Z[qubit] & destabilizer_table_[qubit].X[i]) + num_anti++; + if (P.X[qubit] & destabilizer_table_[qubit].Z[i]) + num_anti++; } - if (num_anti % 2 == 0) continue; + if (num_anti % 2 == 0) + continue; // If anti-commutes multiply Pauli by stabilizer phase += 2 * (uint_t)stabilizer_phases_[i]; @@ -676,13 +685,11 @@ double Clifford::expval_pauli(const reg_t &qubits, return (phase % 4) ? -1.0 : 1.0; } - //------------------------------------------------------------------------------ // JSON Serialization //------------------------------------------------------------------------------ -json_t Clifford::json() const -{ +json_t Clifford::json() const { json_t js = json_t::object(); // Add destabilizers json_t stab; @@ -710,14 +717,12 @@ json_t Clifford::json() const return js; } -inline void to_json(json_t &js, const Clifford &clif) { - js = clif.json(); -} +inline void to_json(json_t &js, const Clifford &clif) { js = clif.json(); } template -void build_from(const inputdata_t& input, Clifford& clif) -{ - bool has_keys = AER::Parser::check_keys({"stabilizer", "destabilizer"}, input); +void build_from(const inputdata_t &input, Clifford &clif) { + bool has_keys = AER::Parser::check_keys( + {"stabilizer", "destabilizer"}, input); if (!has_keys) throw std::invalid_argument("Invalid Clifford JSON."); @@ -727,7 +732,8 @@ void build_from(const inputdata_t& input, Clifford& clif) const auto nq = stab.size(); if (nq != destab.size()) { - throw std::invalid_argument("Invalid Clifford JSON: stabilizer and destabilizer lengths do not match."); + throw std::invalid_argument("Invalid Clifford JSON: stabilizer and " + "destabilizer lengths do not match."); } clif.initialize(nq); @@ -736,57 +742,61 @@ void build_from(const inputdata_t& input, Clifford& clif) // Get destabilizer label = destab[i]; switch (label[0]) { - case '-': - clif.set_destabilizer_phases(i,1); - clif.set_destabilizer(i, Pauli::Pauli(label.substr(1, nq))); - break; - case '+': - clif.set_destabilizer(i, Pauli::Pauli(label.substr(1, nq))); - break; - case 'I': - case 'X': - case 'Y': - case 'Z': - clif.set_destabilizer(i, Pauli::Pauli(label)); - break; - default: - throw std::invalid_argument("Invalid Stabilizer JSON string."); + case '-': + clif.set_destabilizer_phases(i, 1); + clif.set_destabilizer( + i, Pauli::Pauli(label.substr(1, nq))); + break; + case '+': + clif.set_destabilizer( + i, Pauli::Pauli(label.substr(1, nq))); + break; + case 'I': + case 'X': + case 'Y': + case 'Z': + clif.set_destabilizer(i, Pauli::Pauli(label)); + break; + default: + throw std::invalid_argument("Invalid Stabilizer JSON string."); } // Get stabilizer label = stab[i]; switch (label[0]) { - case '-': - clif.set_stabilizer_phases(i, 1); - clif.set_stabilizer(i, Pauli::Pauli(label.substr(1, nq))); - break; - case '+': - clif.set_stabilizer(i, Pauli::Pauli(label.substr(1, nq))); - break; - case 'I': - case 'X': - case 'Y': - case 'Z': - clif.set_stabilizer(i, Pauli::Pauli(label)); - break; - default: - throw std::invalid_argument("Invalid Stabilizer JSON string."); + case '-': + clif.set_stabilizer_phases(i, 1); + clif.set_stabilizer(i, + Pauli::Pauli(label.substr(1, nq))); + break; + case '+': + clif.set_stabilizer(i, + Pauli::Pauli(label.substr(1, nq))); + break; + case 'I': + case 'X': + case 'Y': + case 'Z': + clif.set_stabilizer(i, Pauli::Pauli(label)); + break; + default: + throw std::invalid_argument("Invalid Stabilizer JSON string."); } } } -inline void from_json(const json_t &js, Clifford &clif) -{ +inline void from_json(const json_t &js, Clifford &clif) { build_from(js, clif); } //------------------------------------------------------------------------------ } // end namespace Clifford -} // AER +} // namespace AER //------------------------------------------------------------------------------ // ostream overload for templated qubitvector template -std::ostream &operator<<(std::ostream &out, const AER::Clifford::Clifford &clif) { +std::ostream &operator<<(std::ostream &out, + const AER::Clifford::Clifford &clif) { out << clif.json().dump(); return out; } diff --git a/src/simulators/stabilizer/pauli.hpp b/src/simulators/stabilizer/pauli.hpp index 83ac23d787..39b3aae5fd 100644 --- a/src/simulators/stabilizer/pauli.hpp +++ b/src/simulators/stabilizer/pauli.hpp @@ -15,8 +15,8 @@ #ifndef _pauli_hpp_ #define _pauli_hpp_ -#include #include "binary_vector.hpp" +#include namespace AER { namespace Pauli { @@ -35,10 +35,10 @@ class Pauli { bv_type Z; // Construct an empty pauli - Pauli() {}; + Pauli(){}; // Construct an n-qubit identity Pauli - explicit Pauli(uint64_t len) : X(len), Z(len) {}; + explicit Pauli(uint64_t len) : X(len), Z(len){}; // Construct an n-qubit Pauli from a string label; explicit Pauli(const std::string &label); @@ -47,10 +47,10 @@ class Pauli { std::string str() const; // exponent g of i such that P(x1,z1) P(x2,z2) = i^g P(x1+x2,z1+z2) - static int8_t phase_exponent(const Pauli& pauli1, const Pauli& pauli2); + static int8_t phase_exponent(const Pauli &pauli1, + const Pauli &pauli2); }; - /******************************************************************************* * * Implementations @@ -64,23 +64,23 @@ Pauli::Pauli(const std::string &label) { Z = bv_type(num_qubits); // The label corresponds to tensor product order // So the first element of label is the last qubit and vice versa - for (size_t i =0; i < num_qubits; i++) { + for (size_t i = 0; i < num_qubits; i++) { const auto qubit_i = num_qubits - 1 - i; switch (label[i]) { - case 'I': - break; - case 'X': - X.set1(qubit_i); - break; - case 'Y': - X.set1(qubit_i); - Z.set1(qubit_i); - break; - case 'Z': - Z.set1(qubit_i); - break; - default: - throw std::invalid_argument("Invalid Pauli label"); + case 'I': + break; + case 'X': + X.set1(qubit_i); + break; + case 'Y': + X.set1(qubit_i); + Z.set1(qubit_i); + break; + case 'Z': + Z.set1(qubit_i); + break; + default: + throw std::invalid_argument("Invalid Pauli label"); } } } @@ -90,9 +90,10 @@ std::string Pauli::str() const { // Check X and Z are same length const auto num_qubits = X.getLength(); if (Z.getLength() != num_qubits) - throw std::runtime_error("Pauli::str X and Z vectors are different length."); + throw std::runtime_error( + "Pauli::str X and Z vectors are different length."); std::string label; - for (size_t i =0; i < num_qubits; i++) { + for (size_t i = 0; i < num_qubits; i++) { const auto qubit_i = num_qubits - 1 - i; if (!X[qubit_i] && !Z[qubit_i]) label.push_back('I'); @@ -107,30 +108,33 @@ std::string Pauli::str() const { } template -int8_t Pauli::phase_exponent(const Pauli& pauli1, const Pauli& pauli2) { +int8_t Pauli::phase_exponent(const Pauli &pauli1, + const Pauli &pauli2) { int8_t exponent = 0; for (size_t q = 0; q < pauli1.X.getLength(); q++) { - exponent += pauli2.X[q] * pauli1.Z[q] * (1 + 2 * pauli2.Z[q] + 2 * pauli1.X[q]); - exponent -= pauli1.X[q] * pauli2.Z[q] * (1 + 2 * pauli1.Z[q] + 2 * pauli2.X[q]); + exponent += + pauli2.X[q] * pauli1.Z[q] * (1 + 2 * pauli2.Z[q] + 2 * pauli1.X[q]); + exponent -= + pauli1.X[q] * pauli2.Z[q] * (1 + 2 * pauli1.Z[q] + 2 * pauli2.X[q]); exponent %= 4; } if (exponent < 0) - exponent += 4; + exponent += 4; return exponent; } //------------------------------------------------------------------------------ } // end namespace Pauli -} // AER +} // namespace AER //------------------------------------------------------------------------------ // ostream overload for templated qubitvector template -inline std::ostream &operator<<(std::ostream &out, const AER::Pauli::Pauli &pauli) { +inline std::ostream &operator<<(std::ostream &out, + const AER::Pauli::Pauli &pauli) { out << pauli.str(); return out; } //------------------------------------------------------------------------------ #endif - diff --git a/src/simulators/stabilizer/stabilizer_state.hpp b/src/simulators/stabilizer/stabilizer_state.hpp index b236d0b9cc..8f8fabbcee 100644 --- a/src/simulators/stabilizer/stabilizer_state.hpp +++ b/src/simulators/stabilizer/stabilizer_state.hpp @@ -15,11 +15,11 @@ #ifndef _aer_stabilizer_state_hpp #define _aer_stabilizer_state_hpp -#include "framework/utils.hpp" -#include "framework/json.hpp" +#include "clifford.hpp" #include "framework/config.hpp" +#include "framework/json.hpp" +#include "framework/utils.hpp" #include "simulators/state.hpp" -#include "clifford.hpp" namespace AER { namespace Stabilizer { @@ -31,22 +31,17 @@ using OpType = Operations::OpType; // OpSet of supported instructions const Operations::OpSet StateOpSet( - // Op types - {OpType::gate, OpType::measure, - OpType::reset, - OpType::barrier, OpType::bfunc, OpType::qerror_loc, - OpType::roerror, OpType::save_expval, - OpType::save_expval_var, OpType::save_probs, - OpType::save_probs_ket, OpType::save_amps_sq, - OpType::save_stabilizer, OpType::save_clifford, - OpType::save_state, OpType::set_stabilizer, - OpType::jump, OpType::mark - }, - // Gates - {"CX", "cx", "cy", "cz", "swap", "id", "x", "y", "z", "h", "s", "sdg", - "sx", "sxdg", "delay", "pauli"}); - -enum class Gates {id, x, y, z, h, s, sdg, sx, sxdg, cx, cy, cz, swap, pauli}; + // Op types + {OpType::gate, OpType::measure, OpType::reset, OpType::barrier, + OpType::bfunc, OpType::qerror_loc, OpType::roerror, OpType::save_expval, + OpType::save_expval_var, OpType::save_probs, OpType::save_probs_ket, + OpType::save_amps_sq, OpType::save_stabilizer, OpType::save_clifford, + OpType::save_state, OpType::set_stabilizer, OpType::jump, OpType::mark}, + // Gates + {"CX", "cx", "cy", "cz", "swap", "id", "x", "y", "z", "h", "s", "sdg", "sx", + "sxdg", "delay", "pauli"}); + +enum class Gates { id, x, y, z, h, s, sdg, sx, sxdg, cx, cy, cz, swap, pauli }; //============================================================================ // Stabilizer Table state class @@ -66,35 +61,31 @@ class State : public QuantumState::State { //----------------------------------------------------------------------- // Return the string name of the State class - virtual std::string name() const override {return "stabilizer";} + virtual std::string name() const override { return "stabilizer"; } // Apply an operation // If the op is not in allowed_ops an exeption will be raised. - virtual void apply_op(const Operations::Op &op, - ExperimentResult &result, - RngEngine &rng, - bool final_op = false) override; + virtual void apply_op(const Operations::Op &op, ExperimentResult &result, + RngEngine &rng, bool final_op = false) override; // Initializes an n-qubit state to the all |0> state virtual void initialize_qreg(uint_t num_qubits) override; // TODO: currently returns 0 // Returns the required memory for storing an n-qubit state in megabytes. - virtual size_t required_memory_mb(uint_t num_qubits, - const std::vector &ops) - const override; + virtual size_t + required_memory_mb(uint_t num_qubits, + const std::vector &ops) const override; // Load any settings for the State class from a config JSON virtual void set_config(const Config &config) override; // Sample n-measurement outcomes without applying the measure operation // to the system state - virtual std::vector sample_measure(const reg_t& qubits, - uint_t shots, + virtual std::vector sample_measure(const reg_t &qubits, uint_t shots, RngEngine &rng) override; protected: - //----------------------------------------------------------------------- // Apply instructions //----------------------------------------------------------------------- @@ -105,16 +96,14 @@ class State : public QuantumState::State { // Applies a sypported Gate operation to the state class. // If the input is not in allowed_gates an exeption will be raised. - void apply_pauli(const reg_t &qubits, const std::string& pauli); + void apply_pauli(const reg_t &qubits, const std::string &pauli); // Measure qubits and return a list of outcomes [q0, q1, ...] // If a state subclass supports this function then "measure" // should be contained in the set returned by the 'allowed_ops' // method. - virtual void apply_measure(const reg_t &qubits, - const reg_t &cmemory, - const reg_t &cregister, - RngEngine &rng); + virtual void apply_measure(const reg_t &qubits, const reg_t &cmemory, + const reg_t &cregister, RngEngine &rng); // Reset the specified qubits to the |0> state by simulating // a measurement, applying a conditional x-gate if the outcome is 1, and @@ -129,7 +118,8 @@ class State : public QuantumState::State { //----------------------------------------------------------------------- // Save Clifford state of simulator - void apply_save_stabilizer(const Operations::Op &op, ExperimentResult &result); + void apply_save_stabilizer(const Operations::Op &op, + ExperimentResult &result); // Save probabilities void apply_save_probs(const Operations::Op &op, ExperimentResult &result); @@ -140,22 +130,18 @@ class State : public QuantumState::State { // Helper function for computing expectation value virtual double expval_pauli(const reg_t &qubits, - const std::string& pauli) override; + const std::string &pauli) override; // Return the probability of an outcome bitstring. double get_probability(const reg_t &qubits, const std::string &outcome); template - void get_probabilities_auxiliary(const reg_t& qubits, - std::string outcome, - double outcome_prob, - T& probs); - - void get_probability_helper(const reg_t& qubits, - const std::string &outcome, - std::string &outcome_carry, - double &prob_carry); - + void get_probabilities_auxiliary(const reg_t &qubits, std::string outcome, + double outcome_prob, T &probs); + + void get_probability_helper(const reg_t &qubits, const std::string &outcome, + std::string &outcome_carry, double &prob_carry); + //----------------------------------------------------------------------- // Measurement Helpers //----------------------------------------------------------------------- @@ -176,33 +162,31 @@ class State : public QuantumState::State { // Table of allowed gate names to gate enum class members const static stringmap_t gateset_; - }; - //============================================================================ // Implementation: Allowed ops and gateset //============================================================================ const stringmap_t State::gateset_({ - // Single qubit gates - {"delay", Gates::id},// Delay gate - {"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)) - {"sx", Gates::sx}, // Sqrt X gate. - {"sxdg", Gates::sxdg}, // Inverse Sqrt X gate. - // Two-qubit gates - {"CX", Gates::cx}, // Controlled-X gate (CNOT) - {"cx", Gates::cx}, // Controlled-X gate (CNOT), - {"cy", Gates::cy}, // Controlled-Y gate - {"cz", Gates::cz}, // Controlled-Z gate - {"swap", Gates::swap}, // SWAP gate - {"pauli", Gates::pauli} // Pauli gate + // Single qubit gates + {"delay", Gates::id}, // Delay gate + {"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)) + {"sx", Gates::sx}, // Sqrt X gate. + {"sxdg", Gates::sxdg}, // Inverse Sqrt X gate. + // Two-qubit gates + {"CX", Gates::cx}, // Controlled-X gate (CNOT) + {"cx", Gates::cx}, // Controlled-X gate (CNOT), + {"cy", Gates::cy}, // Controlled-Y gate + {"cz", Gates::cz}, // Controlled-Z gate + {"swap", Gates::swap}, // SWAP gate + {"pauli", Gates::pauli} // Pauli gate }); //============================================================================ @@ -222,8 +206,7 @@ void State::initialize_qreg(uint_t num_qubits) { //------------------------------------------------------------------------- size_t State::required_memory_mb(uint_t num_qubits, - const std::vector &ops) - const { + const std::vector &ops) const { (void)ops; // avoid unused variable compiler warning // The Clifford object requires very little memory. // A Pauli vector consists of 2 binary vectors each with @@ -232,7 +215,7 @@ size_t State::required_memory_mb(uint_t num_qubits, size_t mem = 16 * (4 + num_qubits); // Pauli bytes // Clifford = 2n * Pauli + 2n phase ints mem = 2 * num_qubits * (mem + 16); // Clifford bytes - mem = mem >> 20; // Clifford mb + mem = mem >> 20; // Clifford mb return mem; } @@ -249,142 +232,142 @@ void State::set_config(const Config &config) { // Implementation: apply operations //========================================================================= -void State::apply_op(const Operations::Op &op, - ExperimentResult &result, +void State::apply_op(const Operations::Op &op, ExperimentResult &result, RngEngine &rng, bool final_op) { if (BaseState::creg().check_conditional(op)) { switch (op.type) { - case OpType::barrier: - case OpType::qerror_loc: - break; - case OpType::reset: - apply_reset(op.qubits, rng); - break; - case OpType::measure: - apply_measure(op.qubits, op.memory, op.registers, rng); - break; - case OpType::bfunc: - BaseState::creg().apply_bfunc(op); - break; - case OpType::roerror: - BaseState::creg().apply_roerror(op, rng); - break; - case OpType::gate: - apply_gate(op); - break; - case OpType::set_stabilizer: - apply_set_stabilizer(op.clifford); - break; - case OpType::save_expval: - case OpType::save_expval_var: - apply_save_expval(op, result); - break; - case OpType::save_probs: - case OpType::save_probs_ket: - apply_save_probs(op, result); - break; - case OpType::save_amps_sq: - apply_save_amplitudes_sq(op, result); - break; - case OpType::save_state: - case OpType::save_stabilizer: - case OpType::save_clifford: - apply_save_stabilizer(op, result); - break; - default: - throw std::invalid_argument("Stabilizer::State::invalid instruction \'" + - op.name + "\'."); - } - } -} - -void State::apply_gate(const Operations::Op &op) { - // Check Op is supported by State - auto it = gateset_.find(op.name); - if (it == gateset_.end()) - throw std::invalid_argument("Stabilizer::State::invalid gate instruction \'" + - op.name + "\'."); - switch (it -> second) { - case Gates::id: - break; - case Gates::x: - BaseState::qreg_.append_x(op.qubits[0]); + case OpType::barrier: + case OpType::qerror_loc: break; - case Gates::y: - BaseState::qreg_.append_y(op.qubits[0]); + case OpType::reset: + apply_reset(op.qubits, rng); break; - case Gates::z: - BaseState::qreg_.append_z(op.qubits[0]); + case OpType::measure: + apply_measure(op.qubits, op.memory, op.registers, rng); break; - case Gates::h: - BaseState::qreg_.append_h(op.qubits[0]); + case OpType::bfunc: + BaseState::creg().apply_bfunc(op); break; - case Gates::s: - BaseState::qreg_.append_s(op.qubits[0]); + case OpType::roerror: + BaseState::creg().apply_roerror(op, rng); break; - case Gates::sdg: - BaseState::qreg_.append_z(op.qubits[0]); - BaseState::qreg_.append_s(op.qubits[0]); + case OpType::gate: + apply_gate(op); break; - case Gates::sx: - BaseState::qreg_.append_z(op.qubits[0]); - BaseState::qreg_.append_s(op.qubits[0]); - BaseState::qreg_.append_h(op.qubits[0]); - BaseState::qreg_.append_z(op.qubits[0]); - BaseState::qreg_.append_s(op.qubits[0]); + case OpType::set_stabilizer: + apply_set_stabilizer(op.clifford); break; - case Gates::sxdg: - BaseState::qreg_.append_s(op.qubits[0]); - BaseState::qreg_.append_h(op.qubits[0]); - BaseState::qreg_.append_s(op.qubits[0]); + case OpType::save_expval: + case OpType::save_expval_var: + apply_save_expval(op, result); break; - case Gates::cx: - BaseState::qreg_.append_cx(op.qubits[0], op.qubits[1]); + case OpType::save_probs: + case OpType::save_probs_ket: + apply_save_probs(op, result); break; - case Gates::cz: - BaseState::qreg_.append_h(op.qubits[1]); - BaseState::qreg_.append_cx(op.qubits[0], op.qubits[1]); - BaseState::qreg_.append_h(op.qubits[1]); + case OpType::save_amps_sq: + apply_save_amplitudes_sq(op, result); break; - case Gates::cy: - BaseState::qreg_.append_z(op.qubits[1]); - BaseState::qreg_.append_s(op.qubits[1]); - BaseState::qreg_.append_cx(op.qubits[0], op.qubits[1]); - BaseState::qreg_.append_s(op.qubits[1]); - break; - case Gates::swap: - BaseState::qreg_.append_cx(op.qubits[0], op.qubits[1]); - BaseState::qreg_.append_cx(op.qubits[1], op.qubits[0]); - BaseState::qreg_.append_cx(op.qubits[0], op.qubits[1]); - break; - case Gates::pauli: - apply_pauli(op.qubits, op.string_params[0]); + case OpType::save_state: + case OpType::save_stabilizer: + case OpType::save_clifford: + apply_save_stabilizer(op, result); break; default: - // We shouldn't reach here unless there is a bug in gateset - throw std::invalid_argument("Stabilizer::State::invalid gate instruction \'" + + throw std::invalid_argument("Stabilizer::State::invalid instruction \'" + op.name + "\'."); + } } } -void State::apply_pauli(const reg_t &qubits, const std::string& pauli) { +void State::apply_gate(const Operations::Op &op) { + // Check Op is supported by State + auto it = gateset_.find(op.name); + if (it == gateset_.end()) + throw std::invalid_argument( + "Stabilizer::State::invalid gate instruction \'" + op.name + "\'."); + switch (it->second) { + case Gates::id: + break; + case Gates::x: + BaseState::qreg_.append_x(op.qubits[0]); + break; + case Gates::y: + BaseState::qreg_.append_y(op.qubits[0]); + break; + case Gates::z: + BaseState::qreg_.append_z(op.qubits[0]); + break; + case Gates::h: + BaseState::qreg_.append_h(op.qubits[0]); + break; + case Gates::s: + BaseState::qreg_.append_s(op.qubits[0]); + break; + case Gates::sdg: + BaseState::qreg_.append_z(op.qubits[0]); + BaseState::qreg_.append_s(op.qubits[0]); + break; + case Gates::sx: + BaseState::qreg_.append_z(op.qubits[0]); + BaseState::qreg_.append_s(op.qubits[0]); + BaseState::qreg_.append_h(op.qubits[0]); + BaseState::qreg_.append_z(op.qubits[0]); + BaseState::qreg_.append_s(op.qubits[0]); + break; + case Gates::sxdg: + BaseState::qreg_.append_s(op.qubits[0]); + BaseState::qreg_.append_h(op.qubits[0]); + BaseState::qreg_.append_s(op.qubits[0]); + break; + case Gates::cx: + BaseState::qreg_.append_cx(op.qubits[0], op.qubits[1]); + break; + case Gates::cz: + BaseState::qreg_.append_h(op.qubits[1]); + BaseState::qreg_.append_cx(op.qubits[0], op.qubits[1]); + BaseState::qreg_.append_h(op.qubits[1]); + break; + case Gates::cy: + BaseState::qreg_.append_z(op.qubits[1]); + BaseState::qreg_.append_s(op.qubits[1]); + BaseState::qreg_.append_cx(op.qubits[0], op.qubits[1]); + BaseState::qreg_.append_s(op.qubits[1]); + break; + case Gates::swap: + BaseState::qreg_.append_cx(op.qubits[0], op.qubits[1]); + BaseState::qreg_.append_cx(op.qubits[1], op.qubits[0]); + BaseState::qreg_.append_cx(op.qubits[0], op.qubits[1]); + break; + case Gates::pauli: + apply_pauli(op.qubits, op.string_params[0]); + break; + default: + // We shouldn't reach here unless there is a bug in gateset + throw std::invalid_argument( + "Stabilizer::State::invalid gate instruction \'" + op.name + "\'."); + } +} + +void State::apply_pauli(const reg_t &qubits, const std::string &pauli) { const auto size = qubits.size(); for (size_t i = 0; i < qubits.size(); ++i) { const auto qubit = qubits[size - 1 - i]; switch (pauli[i]) { - case 'I': - break; - case 'X': - BaseState::qreg_.append_x(qubit); - break; - case 'Y': - BaseState::qreg_.append_y(qubit); - break; - case 'Z': - BaseState::qreg_.append_z(qubit); - break; - default: - throw std::invalid_argument("invalid Pauli \'" + std::to_string(pauli[i]) + "\'."); + case 'I': + break; + case 'X': + BaseState::qreg_.append_x(qubit); + break; + case 'Y': + BaseState::qreg_.append_y(qubit); + break; + case 'Z': + BaseState::qreg_.append_z(qubit); + break; + default: + throw std::invalid_argument("invalid Pauli \'" + + std::to_string(pauli[i]) + "\'."); } } } @@ -393,34 +376,26 @@ void State::apply_pauli(const reg_t &qubits, const std::string& pauli) { // Implementation: Reset and Measurement Sampling //========================================================================= - -void State::apply_measure(const reg_t &qubits, - const reg_t &cmemory, - const reg_t &cregister, - RngEngine &rng) { +void State::apply_measure(const reg_t &qubits, const reg_t &cmemory, + const reg_t &cregister, RngEngine &rng) { // Apply measurement and get classical outcome reg_t outcome = apply_measure_and_update(qubits, rng); // Add measurement outcome to creg BaseState::creg().store_measure(outcome, cmemory, cregister); } - void State::apply_reset(const reg_t &qubits, RngEngine &rng) { // Apply measurement and get classical outcome reg_t outcome = apply_measure_and_update(qubits, rng); // Use the outcome to apply X gate to any qubits left in the // |1> state after measure, then discard outcome. - for (size_t j=0; j < qubits.size(); j++) { - if (outcome[j] == 1) { + for (size_t j = 0; j < qubits.size(); j++) + if (outcome[j] == 1) qreg_.append_x(qubits[j]); - } - } } - -reg_t State::apply_measure_and_update(const reg_t &qubits, - RngEngine &rng) { +reg_t State::apply_measure_and_update(const reg_t &qubits, RngEngine &rng) { // Measurement outcome probabilities in the clifford // table are either deterministic or random. // We generate the distribution for the random case @@ -436,10 +411,10 @@ reg_t State::apply_measure_and_update(const reg_t &qubits, return outcome; } -std::vector State::sample_measure(const reg_t &qubits, - uint_t shots, +std::vector State::sample_measure(const reg_t &qubits, uint_t shots, RngEngine &rng) { - // TODO: see if we can improve efficiency by directly sampling from Clifford table + // TODO: see if we can improve efficiency by directly sampling from Clifford + // table auto qreg_cache = BaseState::qreg_; std::vector samples; samples.reserve(shots); @@ -453,9 +428,9 @@ std::vector State::sample_measure(const reg_t &qubits, void State::apply_set_stabilizer(const Clifford::Clifford &clifford) { if (clifford.num_qubits() != BaseState::qreg_.num_qubits()) { throw std::invalid_argument( - "set stabilizer must be defined on full width of qubits (" + - std::to_string(clifford.num_qubits()) + " != " + - std::to_string(BaseState::qreg_.num_qubits()) + ")."); + "set stabilizer must be defined on full width of qubits (" + + std::to_string(clifford.num_qubits()) + + " != " + std::to_string(BaseState::qreg_.num_qubits()) + ")."); } BaseState::qreg_.apply_set_stabilizer(clifford); } @@ -465,30 +440,30 @@ void State::apply_set_stabilizer(const Clifford::Clifford &clifford) { //========================================================================= void State::apply_save_stabilizer(const Operations::Op &op, - ExperimentResult &result) { + ExperimentResult &result) { std::string key = op.string_params[0]; OpType op_type = op.type; switch (op_type) { - case OpType::save_clifford: { - if (key == "_method_") { - key = "clifford"; - } - break; - } - case OpType::save_stabilizer: - case OpType::save_state: { - if (key == "_method_") { - key = "stabilizer"; - } - op_type = OpType::save_stabilizer; - break; - } - default: - // We shouldn't ever reach here... - throw std::invalid_argument("Invalid save state instruction for stabilizer"); + case OpType::save_clifford: { + if (key == "_method_") + key = "clifford"; + break; + } + case OpType::save_stabilizer: + case OpType::save_state: { + if (key == "_method_") + key = "stabilizer"; + op_type = OpType::save_stabilizer; + break; + } + default: + // We shouldn't ever reach here... + throw std::invalid_argument( + "Invalid save state instruction for stabilizer"); } json_t clifford = BaseState::qreg_; - result.save_data_pershot(creg(), key, std::move(clifford), op_type, op.save_type); + result.save_data_pershot(creg(), key, std::move(clifford), op_type, + op.save_type); } void State::apply_save_probs(const Operations::Op &op, @@ -500,68 +475,67 @@ void State::apply_save_probs(const Operations::Op &op, // to store. const size_t num_qubits = op.qubits.size(); if (num_qubits > max_qubits_snapshot_probs_) { - std::string msg = - "Stabilizer::State::snapshot_probabilities: " - "cannot return measure probabilities for " + - std::to_string(num_qubits) + "-qubit measurement. Maximum is set to " + - std::to_string(max_qubits_snapshot_probs_); + std::string msg = "Stabilizer::State::snapshot_probabilities: " + "cannot return measure probabilities for " + + std::to_string(num_qubits) + + "-qubit measurement. Maximum is set to " + + std::to_string(max_qubits_snapshot_probs_); throw std::runtime_error(msg); } if (op.type == OpType::save_probs_ket) { std::map probs; - get_probabilities_auxiliary( - op.qubits, std::string(op.qubits.size(), 'X'), 1, probs); - result.save_data_average(creg(), op.string_params[0], - std::move(probs), op.type, op.save_type); + get_probabilities_auxiliary(op.qubits, std::string(op.qubits.size(), 'X'), + 1, probs); + result.save_data_average(creg(), op.string_params[0], std::move(probs), + op.type, op.save_type); } else { std::vector probs(1ULL << op.qubits.size(), 0.); - get_probabilities_auxiliary( - op.qubits, std::string(op.qubits.size(), 'X'), 1, probs); - result.save_data_average(creg(), op.string_params[0], - std::move(probs), op.type, op.save_type); + get_probabilities_auxiliary(op.qubits, std::string(op.qubits.size(), 'X'), + 1, probs); + result.save_data_average(creg(), op.string_params[0], std::move(probs), + op.type, op.save_type); } } void State::apply_save_amplitudes_sq(const Operations::Op &op, ExperimentResult &result) { if (op.int_params.empty()) { - throw std::invalid_argument("Invalid save_amplitudes_sq instructions (empty params)."); + throw std::invalid_argument( + "Invalid save_amplitudes_sq instructions (empty params)."); } uint_t num_qubits = op.qubits.size(); if (num_qubits != BaseState::qreg_.num_qubits()) { - throw std::invalid_argument("Save amplitude square must be defined on full width of qubits."); + throw std::invalid_argument( + "Save amplitude square must be defined on full width of qubits."); } - rvector_t amps_sq(op.int_params.size(), 1.0); // Must be initialized in 1 for helper func + rvector_t amps_sq(op.int_params.size(), + 1.0); // Must be initialized in 1 for helper func for (size_t i = 0; i < op.int_params.size(); i++) { - amps_sq[i] = get_probability(op.qubits, Utils::int2bin(op.int_params[i], num_qubits)); + amps_sq[i] = get_probability(op.qubits, + Utils::int2bin(op.int_params[i], num_qubits)); } - result.save_data_average(creg(), op.string_params[0], - std::move(amps_sq), op.type, op.save_type); + result.save_data_average(creg(), op.string_params[0], std::move(amps_sq), + op.type, op.save_type); } -double State::expval_pauli(const reg_t &qubits, - const std::string& pauli) -{ +double State::expval_pauli(const reg_t &qubits, const std::string &pauli) { return BaseState::qreg_.expval_pauli(qubits, pauli); } -static void set_value_helper(std::map& probs, - const std::string &outcome, - double value) { +static void set_value_helper(std::map &probs, + const std::string &outcome, double value) { probs[Utils::bin2hex(outcome)] = value; } -static void set_value_helper(std::vector& probs, - const std::string &outcome, - double value) { +static void set_value_helper(std::vector &probs, + const std::string &outcome, double value) { probs[std::stoull(outcome, 0, 2)] = value; } template void State::get_probabilities_auxiliary(const reg_t &qubits, std::string outcome, - double outcome_prob, - T &probs) { + double outcome_prob, T &probs) { uint_t qubit_for_branching = -1; for (uint_t i = 0; i < qubits.size(); ++i) { uint_t qubit = qubits[qubits.size() - i - 1]; @@ -569,11 +543,10 @@ void State::get_probabilities_auxiliary(const reg_t &qubits, if (BaseState::qreg_.is_deterministic_outcome(qubit)) { bool single_qubit_outcome = BaseState::qreg_.measure_and_update(qubit, 0); - if (single_qubit_outcome) { + if (single_qubit_outcome) outcome[i] = '1'; - } else { + else outcome[i] = '0'; - } } else { qubit_for_branching = i; } @@ -588,11 +561,10 @@ void State::get_probabilities_auxiliary(const reg_t &qubits, for (uint_t single_qubit_outcome = 0; single_qubit_outcome < 2; ++single_qubit_outcome) { std::string new_outcome = outcome; - if (single_qubit_outcome) { + if (single_qubit_outcome) new_outcome[qubit_for_branching] = '1'; - } else { + else new_outcome[qubit_for_branching] = '0'; - } auto copy_of_qreg = BaseState::qreg_; BaseState::qreg_.measure_and_update( @@ -620,11 +592,10 @@ void State::get_probability_helper(const reg_t &qubits, if (BaseState::qreg_.is_deterministic_outcome(qubit)) { bool single_qubit_outcome = BaseState::qreg_.measure_and_update(qubit, 0); - if (single_qubit_outcome) { + if (single_qubit_outcome) outcome_carry[i] = '1'; - } else { + else outcome_carry[i] = '0'; - } if (outcome[i] != outcome_carry[i]) { prob_carry = 0.0; return; @@ -634,9 +605,8 @@ void State::get_probability_helper(const reg_t &qubits, } } } - if (qubit_for_branching == -1) { + if (qubit_for_branching == -1) return; - } outcome_carry[qubit_for_branching] = outcome[qubit_for_branching]; uint_t single_qubit_outcome = (outcome[qubit_for_branching] == '1') ? 1 : 0; auto cached_qreg = BaseState::qreg_;