diff --git a/qiskit/providers/aer/backends/aer_simulator.py b/qiskit/providers/aer/backends/aer_simulator.py index 4846a9e0f6..8fb4ea05f3 100644 --- a/qiskit/providers/aer/backends/aer_simulator.py +++ b/qiskit/providers/aer/backends/aer_simulator.py @@ -139,6 +139,8 @@ class AerSimulator(AerBackend): +--------------------------+---------------+ | ``extended_stabilizer`` | No | +--------------------------+---------------+ + | ``compute`` | No | + +--------------------------+---------------+ | ``unitary`` | Yes | +--------------------------+---------------+ | ``superop`` | No | @@ -409,6 +411,9 @@ class AerSimulator(AerBackend): 'extended_stabilizer': sorted([ 'quantum_channel', 'qerror_loc', 'roerror', 'snapshot', 'save_statevector' ]), + 'clifford_phase_compute': sorted([ + 'save_specific_probability' + ]), 'unitary': sorted([ 'snapshot', 'save_state', 'save_unitary', 'set_unitary' ]), @@ -449,7 +454,7 @@ class AerSimulator(AerBackend): _SIMULATION_METHODS = [ 'automatic', 'statevector', 'density_matrix', 'stabilizer', 'matrix_product_state', 'extended_stabilizer', - 'unitary', 'superop' + 'unitary', 'superop', 'clifford_phase_compute' ] _AVAILABLE_METHODS = None diff --git a/qiskit/providers/aer/backends/backend_utils.py b/qiskit/providers/aer/backends/backend_utils.py index 110ed4383d..e62ae4af0d 100644 --- a/qiskit/providers/aer/backends/backend_utils.py +++ b/qiskit/providers/aer/backends/backend_utils.py @@ -79,6 +79,10 @@ 'cx', 'cz', 'id', 'x', 'y', 'z', 'h', 's', 'sdg', 'sx', 'sxdg', 'swap', 'u0', 't', 'tdg', 'u1', 'p', 'ccx', 'ccz', 'delay', 'pauli' ]), + 'clifford_phase_compute': sorted([ + 'cx', 'cz', 'id', 'x', 'y', 'z', 'h', 's', 'sdg', + 'swap', 't', 'tdg', 'delay', 'p', 'u1', 'rz' + ]), 'unitary': sorted([ 'u1', 'u2', 'u3', 'u', 'p', 'r', 'rx', 'ry', 'rz', 'id', 'x', 'y', 'z', 'h', 's', 'sdg', 'sx', 'sxdg', 't', 'tdg', 'swap', 'cx', diff --git a/qiskit/providers/aer/library/save_instructions/save_probabilities.py b/qiskit/providers/aer/library/save_instructions/save_probabilities.py index caabbe53c6..34701d10dc 100644 --- a/qiskit/providers/aer/library/save_instructions/save_probabilities.py +++ b/qiskit/providers/aer/library/save_instructions/save_probabilities.py @@ -74,6 +74,38 @@ def __init__(self, conditional=conditional) +class SaveSpecificProbability(SaveAverageData): + """Save a probability for a specific measurement probability.""" + def __init__(self, num_qubits, + states, qubits, + label="specific-probabilities", + pershot=False, + conditional=False): + """Instruction to save specific probabilities. + + Args: + states (list): list of ints indicating the outcome to compute the probability for. + qubits (list): list of ints indicating which qubits the measurement is on. + num_qubits (int): the number of qubits for the snapshot type. + label (str): the key for retrieving saved data from results. + pershot (bool): if True save a list of probabilities for each shot + of the simulation rather than the average over + all shots [Default: False]. + conditional (bool): if True save the probabilities data conditional + on the current classical register values + [Default: False]. + e.g: + if states = [0,1,0], qubits = [0,1,2] + we compute the probability of the outcome 0 on qubit 0, 1 on qubit 1 and 0 on qubit 2 + if states = [0,1], qubits = [5,1] + we compute the probability of the outcome 0 on qubit 5 and 1 on qubit 0 + """ + super().__init__("save_specific_prob", num_qubits, label, + pershot=pershot, + conditional=conditional, + params=[qubits, states]) + + def save_probabilities(self, qubits=None, label="probabilities", @@ -140,5 +172,38 @@ def save_probabilities_dict(self, return self.append(instr, qubits) +def save_specific_probability(self, states, qubits, label="specific_probability", + pershot=False, + conditional=False): + """Instruction to save specific probabilities. + + Args: + states (list): list of ints indicating the outcome to compute the probability for + qubits (list): list of ints indicating which qubits the measurement is on + label (str): the key for retrieving saved data from results. + pershot (bool): if True save a list of probabilities for each shot + of the simulation rather than the average over + all shots [Default: False]. + conditional (bool): if True save the probabilities data conditional + on the current classical register values + [Default: False]. + e.g: + if states = [0,1,0], qubits = [0,1,2] + we compute the probability of 0 on qubit 0, 1 on qubit 1 and 0 on qubit 2 + if states = [0,1], qubits = [5,1] + we compute the probability of 0 on qubit 5 and 1 on qubit 0 + + Returns: + QuantumCircuit: with attached instruction. + """ + if qubits is None: + qubits = default_qubits(self) + instr = SaveSpecificProbability(len(qubits), states, qubits, label=label, + pershot=pershot, + conditional=conditional) + return self.append(instr, qubits) + + QuantumCircuit.save_probabilities = save_probabilities QuantumCircuit.save_probabilities_dict = save_probabilities_dict +QuantumCircuit.save_specific_probability = save_specific_probability diff --git a/src/controllers/aer_controller.hpp b/src/controllers/aer_controller.hpp index 2f1f35639f..39363e09cc 100755 --- a/src/controllers/aer_controller.hpp +++ b/src/controllers/aer_controller.hpp @@ -59,6 +59,7 @@ #include "simulators/statevector/statevector_state.hpp" #include "simulators/superoperator/superoperator_state.hpp" #include "simulators/unitary/unitary_state.hpp" +#include "simulators/clifford_plus_phase/compute.hpp" namespace AER { @@ -111,6 +112,7 @@ class Controller { matrix_product_state, stabilizer, extended_stabilizer, + clifford_phase_compute, unitary, superop }; @@ -127,6 +129,7 @@ class Controller { {Method::matrix_product_state, "matrix_product_state"}, {Method::stabilizer, "stabilizer"}, {Method::extended_stabilizer, "extended_stabilizer"}, + {Method::clifford_phase_compute, "clifford_phase_compute"}, {Method::unitary, "unitary"}, {Method::superop, "superop"} }; @@ -477,6 +480,8 @@ void Controller::set_config(const json_t &config) { method_ = Method::stabilizer; } else if (method == "extended_stabilizer") { method_ = Method::extended_stabilizer; + } else if (method == "clifford_phase_compute") { + method_ = Method::clifford_phase_compute; } else if (method == "matrix_product_state") { method_ = Method::matrix_product_state; } else if (method == "unitary") { @@ -673,6 +678,8 @@ void Controller::set_parallelization_circuit(const Circuit &circ, } case Method::extended_stabilizer: break; + case Method::clifford_phase_compute: + break; default: throw std::invalid_argument("Cannot set parallelization for unresolved method."); } @@ -1136,7 +1143,10 @@ void Controller::run_circuit(const Circuit &circ, const Noise::NoiseModel &noise circ, noise, config, Method::extended_stabilizer, result); case Method::matrix_product_state: return run_circuit_helper( - circ, noise, config, Method::matrix_product_state, result); + circ, noise, config, Method::matrix_product_state, result); + case Method::clifford_phase_compute: + return run_circuit_helper( + circ, noise, config, Method::clifford_phase_compute, result); default: throw std::runtime_error("Controller:Invalid simulation method"); } @@ -1145,7 +1155,6 @@ void Controller::run_circuit(const Circuit &circ, const Noise::NoiseModel &noise //------------------------------------------------------------------------- // Utility methods //------------------------------------------------------------------------- - size_t Controller::required_memory_mb(const Circuit &circ, const Noise::NoiseModel &noise, const Method method) const { @@ -1198,6 +1207,10 @@ size_t Controller::required_memory_mb(const Circuit &circ, MatrixProductState::State state; return state.required_memory_mb(circ.num_qubits, circ.ops); } + case Method::clifford_phase_compute: { + CliffPhaseCompute::State state; + return state.required_memory_mb(circ.num_qubits, circ.ops); + } default: // We shouldn't get here, so throw an exception if we do throw std::runtime_error("Controller: Invalid simulation method"); @@ -1674,6 +1687,10 @@ bool Controller::check_measure_sampling_opt(const Circuit &circ, method == Method::unitary) { return true; } + + if(method == Method::clifford_phase_compute){ + return false; + } // If circuit contains a non-initial initialize that is not a full width // instruction we can't sample @@ -1884,6 +1901,8 @@ bool Controller::validate_method(Method method, return validate_state(Stabilizer::State(), circ, noise_model, throw_except); case Method::extended_stabilizer: return validate_state(ExtendedStabilizer::State(), circ, noise_model, throw_except); + case Method::clifford_phase_compute: + return validate_state(CliffPhaseCompute::State(), circ, noise_model, throw_except); case Method::matrix_product_state: return validate_state(MatrixProductState::State(), circ, noise_model, throw_except); case Method::statevector: diff --git a/src/framework/circuit.hpp b/src/framework/circuit.hpp index af6148c3a5..82de091874 100755 --- a/src/framework/circuit.hpp +++ b/src/framework/circuit.hpp @@ -177,12 +177,15 @@ Circuit::Circuit(const inputdata_t &circ, const json_t &qobj_config, bool trunca // conversion we could call `get_reversed_ops` on the inputdata without first // converting. std::vector converted_ops; + int i = 0; for(auto the_op: input_ops){ + i += 1; converted_ops.emplace_back(Operations::input_to_op(the_op)); } - ops = std::move(converted_ops); + + ops = std::move(converted_ops); set_params(truncation); - + // Check for specified memory slots uint_t memory_slots = 0; Parser::get_value(memory_slots, "memory_slots", config); @@ -205,7 +208,7 @@ Circuit::Circuit(const inputdata_t &circ, const json_t &qobj_config, bool trunca // is explicitly disabled. num_qubits = n_qubits; } - } + } } //------------------------------------------------------------------------- @@ -354,6 +357,7 @@ void Circuit::set_params(bool truncation) { case OpType::save_probs: case OpType::save_probs_ket: case OpType::save_amps: + case OpType::save_specific_prob: case OpType::save_amps_sq: case OpType::save_stabilizer: case OpType::save_clifford: @@ -490,6 +494,7 @@ bool Circuit::check_result_ancestor(const Op& op, std::unordered_set& an case OpType::save_densmat: case OpType::save_probs: case OpType::save_probs_ket: + case OpType::save_specific_prob: case OpType::save_amps: case OpType::save_amps_sq: case OpType::save_stabilizer: diff --git a/src/framework/operations.hpp b/src/framework/operations.hpp index b64e256728..cbbe51b8a8 100755 --- a/src/framework/operations.hpp +++ b/src/framework/operations.hpp @@ -44,7 +44,7 @@ enum class OpType { // Save instructions save_state, save_expval, save_expval_var, save_statevec, save_statevec_dict, save_densmat, save_probs, save_probs_ket, save_amps, save_amps_sq, - save_stabilizer, save_clifford, save_unitary, save_mps, save_superop, + save_stabilizer, save_clifford, save_unitary, save_mps, save_superop, save_specific_prob, // Set instructions set_statevec, set_densmat, set_unitary, set_superop, set_stabilizer, set_mps, @@ -61,8 +61,8 @@ static const std::unordered_set SAVE_TYPES = { OpType::save_statevec, OpType::save_statevec_dict, OpType::save_densmat, OpType::save_probs, OpType::save_probs_ket, OpType::save_amps, OpType::save_amps_sq, OpType::save_stabilizer, - OpType::save_clifford, - OpType::save_unitary, OpType::save_mps, OpType::save_superop + OpType::save_clifford, OpType::save_unitary, OpType::save_mps, + OpType::save_superop, OpType::save_specific_prob }; inline std::ostream& operator<<(std::ostream& stream, const OpType& type) { @@ -126,6 +126,9 @@ inline std::ostream& operator<<(std::ostream& stream, const OpType& type) { case OpType::save_superop: stream << "save_superop"; break; + case OpType::save_specific_prob: + stream << "save_specific_prob"; + break; case OpType::set_statevec: stream << "set_statevector"; break; @@ -576,6 +579,8 @@ template Op input_to_op_save_expval(const inputdata_t& input, bool variance); template Op input_to_op_save_amps(const inputdata_t& input, bool squared); +template +Op input_to_op_save_specific_prob(const inputdata_t& input); // Snapshots template @@ -678,6 +683,9 @@ Op input_to_op(const inputdata_t& input) { return input_to_op_save_amps(input, false); if (name == "save_amplitudes_sq") return input_to_op_save_amps(input, true); + if (name == "save_specific_prob"){ + return input_to_op_save_specific_prob(input); + } // Set if (name == "set_statevector") return input_to_op_set_vector(input, OpType::set_statevec); @@ -1237,6 +1245,17 @@ Op input_to_op_save_amps(const inputdata_t& input, bool squared) { return op; } +template +Op input_to_op_save_specific_prob(const inputdata_t& input) { + Op op = input_to_op_save_default(input, OpType::save_specific_prob); + const inputdata_t& params = Parser::get_value("params", input); + op.qubits = Parser::template get_list_elem>(params, 0); + op.int_params = Parser::template get_list_elem>(params, 1); + + return op; +} + + //------------------------------------------------------------------------------ // Implementation: Snapshot deserialization //------------------------------------------------------------------------------ diff --git a/src/simulators/clifford_plus_phase/ag_state.hpp b/src/simulators/clifford_plus_phase/ag_state.hpp new file mode 100644 index 0000000000..b610ef6a84 --- /dev/null +++ b/src/simulators/clifford_plus_phase/ag_state.hpp @@ -0,0 +1,914 @@ +#ifndef _aer_extended_stabilizer_estimator_ag_hpp +#define _aer_extended_stabilizer_estimator_ag_hpp + +#include +#include +#include +#include +#include "framework/types.hpp" +#include "framework/operations.hpp" +#include "simulators/stabilizer/pauli.hpp" + +namespace AER{ +namespace CliffPhaseCompute{ + +const double T_ANGLE = M_PI/4.; +const double AG_CHOP_THRESHOLD = 1e-10; //if phases are closer to Clifford phases than this we say they are Clifford +const uint_t ONE = 1u; +/* + * We need a variant of the stabilizer tableau described in + * Improved Simulation of Stabilizer Circuits + * Scott Aaronson, Daniel Gottesman (2004) + * + * We only store the stabilizer part of the tableau because we don't need the destabilizer part + * We also add some extra subroutines they don't have which are useful for the Estimate and Compute algorithms + */ + +class AGState{ +public: + //initialize such that the jth stabilizer is the Pauli Z_j + AGState(uint_t qubits, uint_t stabilizers) : num_qubits(qubits), num_stabilizers(stabilizers) {}; + AGState() : num_qubits(0), num_stabilizers(0) {}; + AGState(uint_t qubits) : num_qubits(qubits), num_stabilizers(qubits) {}; + + void initialize(); + void initialize(uint_t qubits); + + size_t num_qubits; + size_t num_stabilizers; //we can represent mixed states so we may have fewer stabilizers than qubits + + std::vector table; + std::vector< unsigned char > phases; + + std::vector magic_phases; //each T or non-Clifford Z-rotation gate we add comes with a phase in [0, pi/4) + virtual ~AGState() = default; //TODO I'm not actually certain the default destructor is correct here + + void Print(); + + // CX controlled on a and targetted on b + void applyCX(size_t a, size_t b); + // CZ controlled on a and targetted on b + void applyCZ(size_t a, size_t b); + void applyH(size_t a); + void applyS(size_t a); + void applyX(size_t a); + void applyY(size_t a); + void applyZ(size_t a); + void applySwap(size_t a, size_t b); + //this adds a new magic qubit initialised in the |0> state and applies the CX gate required for out gadget + void gadgetized_phase_gate(size_t a, double phase); + + // "add" row i onto row h + // this group operation is equivalent to the multiplication of Pauli matrices + void rowsum(size_t h, size_t i); + + /* Does not change the table at all + * updates the "tableau row" given to be row * table_i + * returns the update to the phase of row (0 or 1) + */ + bool rowsum2(Pauli::Pauli row, bool phase, size_t i); + + //return the bool you get at index (stabilizer, column) of the matrix obtained by joining the X and Z matrix as two blocks + //i.e. if column < num_qubits return X[stabilizer][column] otherwise return Z[stabilizer][column-num_qubits] + bool tableau_element(size_t stabilizer, size_t column); + + std::pair first_non_zero_in_col(size_t col, size_t startpoint); + std::pair first_non_zero_in_row(size_t col, size_t startpoint); + + void swap_rows(size_t i, size_t j); + void delete_last_row(); +/* + * create a QCircuit that brings the state represented by this table to the state |0><0|^k \otimes I^(n-k) / (2^(n-k)) + * Explicitly num_stabilizers stabilisers on num_qubits qubits with the jth stabilzier being a +z on the jth qubit + * Note that in addition to returning the circuit that does the simplification this method also brings the state into the simplified form + * if you need the state after calling this method then either use the circuit to rebuild it, or make a copy + */ + std::vector simplifying_unitary(); + + /* + * attempt to find a Pauli X operator acting on qubit q in a stabilizer with an index at least a and at most this->num_stabilizer-2 (we ignore the last stabilizer) + */ + std::pair find_x_on_q(size_t q, size_t a); + /* + * attempt to find a Pauli Y operator acting on qubit q in a stabilizer with an index at least a and at most this->num_stabilizer-2 (we ignore the last stabilizer) + */ + std::pair find_y_on_q(size_t q, size_t a); + /* + * attempt to find a Pauli Z operator acting on qubit q in a stabilizer with an index at least a and at most this->num_stabilizer-2 (we ignore the last stabilizer) + */ + std::pair find_z_on_q(size_t q, size_t a); + + bool independence_test(int q); + + /* + * Apply constraints arising from the fact that we measure the first w qubits and project the last t onto T gates + * In particular we remove stabilisers if qubits [0, w) get killed by taking the expectation value <0| P |0> + * and we remove stabilisers if qubits in [w, this->num_qubits-t) aren't the identity + * we do not remove any qubits from the table + */ + std::pair apply_constraints(size_t w, size_t t); + size_t apply_T_constraints(); + /* + * Go through our stabilizer table and delete every qubit for which every stabilizer is the identity on that qubit + * In other words delete every column from the X and Z matrices if both are 0 for every element in that column + * intended only for use when we have restricted our table to only have magic qubits + * also deletes magic_phases elements to reflect deletion of the magic qubits + */ + void delete_identity_magic_qubits(); + +}; + +// Implementation +void AGState::initialize() +{ + this->table = std::vector(); + + for(size_t i = 0; i < this->num_stabilizers; i++){ + this->table.push_back(Pauli::Pauli(num_qubits)); + this->table[i].Z.set1(i); + this->phases.push_back(0); + } +} + +// Implementation +void AGState::initialize(uint_t qubits) +{ + this->num_qubits = qubits; + this->num_stabilizers = qubits; + this->initialize(); +} + + +void AGState::Print(){ + for(size_t i = 0; i < this->num_stabilizers; i++){ + for(size_t j = 0; j < this->num_qubits; j++){ + std::cout << this->table[i].X[j]; + } + std::cout << "|"; + for(size_t j = 0; j < this->num_qubits; j++){ + std::cout << this->table[i].Z[j]; + } + std::cout << " " << this->phases[i] << std::endl; + } +} + +void AGState::applyCX(size_t a, size_t b){ + for(size_t i = 0; i < this->num_stabilizers; i++){ + this->phases[i] ^= this->table[i].X[a] & this->table[i].Z[b] & (this->table[i].X[b] ^ this->table[i].Z[a] ^ true); + this->table[i].X.xorAt(this->table[i].X[a], b); + this->table[i].Z.xorAt(this->table[i].Z[b], a); + } +} + +void AGState::applyCZ(size_t a, size_t b){ + for(size_t i = 0; i < this->num_stabilizers; i++){ + this->phases[i] ^= this->table[i].X[a] & this->table[i].X[b] & (this->table[i].Z[a] ^ this->table[i].Z[b]); + this->table[i].Z.xorAt(this->table[i].X[b], a); + this->table[i].Z.xorAt(this->table[i].X[a], b); + } +} + +void AGState::applySwap(size_t a, size_t b){ + this->applyCX(a,b); + this->applyCX(b,a); + this->applyCX(a,b); + // for(size_t i = 0; i < this->num_stabilizers; i++){ + // this->table[i].X.xorAt(this->table[i].X[a], b); + // this->table[i].X.xorAt(this->table[i].X[b], a); + // this->table[i].X.xorAt(this->table[i].X[a], b); + + // this->table[i].Z.xorAt(this->table[i].Z[b], a); + // this->table[i].Z.xorAt(this->table[i].Z[a], b); + // this->table[i].Z.xorAt(this->table[i].Z[b], a); + // } +} + +void AGState::applyH(size_t a){ + bool scratch; + for(size_t i = 0; i < this->num_stabilizers; i++){ + this->phases[i] ^= this->table[i].X[a] & this->table[i].Z[a]; + scratch = this->table[i].X[a]; + this->table[i].X.setValue(this->table[i].Z[a], a); + this->table[i].Z.setValue(scratch, a); + } +} + +void AGState::applyS(size_t a){ + for(size_t i = 0; i < this->num_stabilizers; i++){ + this->phases[i] ^= (this->table[i].X[a] & this->table[i].Z[a]); + this->table[i].Z.xorAt(this->table[i].X[a], a); + } +} + +void AGState::applyX(size_t a){ + //TODO write a proper X implementation + this->applyH(a); + this->applyS(a); + this->applyS(a); + this->applyH(a); +} + +void AGState::applyY(size_t a){ + //TODO write a proper Y implementation + this->applyS(a); + this->applyX(a); + this->applyS(a); + this->applyS(a); + this->applyS(a); +} + +void AGState::applyZ(size_t a){ + //TODO write a proper Z implementation + this->applyS(a); + this->applyS(a); +} + +void AGState::gadgetized_phase_gate(size_t a, double phase){ + //phase = fmod(phase , M_PI*2); + while(phase > (M_PI*2)){ + phase -= (M_PI*2); + } + while(phase < 0){ + phase += (M_PI*2); + } + + //now phase is in [0, 2*pi) + while(phase > M_PI/2.){ + phase -= (M_PI/2.); + this->applyS(a); + } + + //now phase is in [0, M_PI/2.] + if(fabs(phase) < AG_CHOP_THRESHOLD){ + //phase on gate is zero so it is an identity + }else if(fabs(phase - M_PI/2.) < AG_CHOP_THRESHOLD){ + //phase on gate is pi/2 so it is an S + this->applyS(a); + + }else{ + //its actually non-Clifford + for(size_t i = 0; i < this->num_stabilizers; i++){ + this->table[i].X.resize(this->num_qubits + 1); + this->table[i].Z.resize(this->num_qubits + 1); + } + this->table.push_back(Pauli::Pauli(this->num_qubits + 1)); + this->table[this->num_stabilizers].Z.set0(this->num_qubits); + this->table[this->num_stabilizers].Z.set1(this->num_qubits); + this->phases.push_back(0); + this->num_stabilizers += 1; + this->num_qubits += 1; + + //we want our phases to be in [0, pi/4] + if(phase > T_ANGLE){ + phase = M_PI/2 - phase; + this->applyS(a); + this->applyX(a); + this->applyCX(a, this->num_qubits-1); + this->applyX(a); + }else{ + this->applyCX(a, this->num_qubits-1); + } + this->magic_phases.push_back(phase); + } +} + +void AGState::rowsum(size_t h, size_t i){ + unsigned char sum = 2*(this->phases[h] + this->phases[i]) + Pauli::Pauli::phase_exponent(this->table[i], this->table[h]); + + sum %= 4u; + + if(sum == 0){ + this->phases[h] = false; + } + if(sum == 2){ + this->phases[h] = true; + } + + this->table[h].X += this->table[i].X; + this->table[h].Z += this->table[i].Z; + +} +bool AGState::rowsum2(Pauli::Pauli row, bool phase, size_t i){ + unsigned char sum = 2*(phase + this->phases[i]) + Pauli::Pauli::phase_exponent(this->table[i], row); + sum %= 4u; + row.X += this->table[i].X; + row.Z += this->table[i].Z; + + if(sum == 0){ + return false; + } + if(sum == 2){ + return true; + } + + //we should never reach here - maybe printing a warning or throwing an exception would be sensible + return false; +} + + +void AGState::delete_identity_magic_qubits(){ + //indended for use when we've already restricted our table to only contain the magic qubits + size_t qubits_deleted = 0; + + for(size_t q = 0; q < this->num_qubits; q++){ + size_t non_identity_paulis = 0; + for(size_t s = 0; (s < this->num_stabilizers) && (non_identity_paulis == 0); s++){ + if(this->table[s].X[q] || this->table[s].Z[q]){ + non_identity_paulis += 1; + } + } + if(non_identity_paulis == 0){ + //every stabiliser is identity on this qubit + //so we can just delete this qubit + qubits_deleted += 1; + }else{ + if(qubits_deleted > 0){ + for(size_t s = 0; s < this->num_stabilizers; s++){ + this->table[s].X.setValue(this->table[s].X[q], q-qubits_deleted); + this->table[s].Z.setValue(this->table[s].Z[q], q-qubits_deleted); + } + magic_phases[q - qubits_deleted] = magic_phases[q]; + } + } + } + this->num_qubits = this->num_qubits - qubits_deleted; + + for(size_t s = 0; s < this->num_stabilizers; s++){ + this->table[s].X.resize(this->num_qubits); + this->table[s].Z.resize(this->num_qubits); + } + magic_phases.resize(this->num_qubits); +} + +bool AGState::tableau_element(size_t stabilizer, size_t column){ + if(column < this->num_qubits){ + return this->table[stabilizer].X[column]; + }else{ + return this->table[stabilizer].Z[column-num_qubits]; + } +} + +std::pair AGState::first_non_zero_in_col(size_t col, size_t startpoint){ + for(size_t i = startpoint; i < this->num_stabilizers; i++){ + if(this->tableau_element(i, col)){ + return std::pair(true, i); + } + } + return std::pair(false, 0); +} +std::pair AGState::first_non_zero_in_row(size_t row, size_t startpoint){ + for(size_t i = startpoint; i < 2*this->num_qubits; i++){ + if(this->tableau_element(row, i) != 0){ + return std::pair(true, i); + } + } + + return std::pair(false, 0); +} + +void AGState::swap_rows(size_t i, size_t j){ + + //this->table[i].X.swap(this->table[j].X); + //this->table[i].Z.swap(this->table[j].Z); + if(i != j){ + std::swap(this->table[i], this->table[j]); + std::swap(this->phases[i], this->phases[j]); + } +} + +void AGState::delete_last_row(){ + this->num_stabilizers -= 1; + this->table.resize(this->num_stabilizers); + this->phases.resize(this->num_stabilizers); +} + +Operations::Op make_H(uint_t a){ + Operations::Op op; + op.type = Operations::OpType::gate; + op.name = "h"; + op.qubits = {a}; + return op; +} + +Operations::Op make_S(uint_t a){ + Operations::Op op; + op.type = Operations::OpType::gate; + op.name = "s"; + op.qubits = {a}; + return op; +} + +Operations::Op make_CX(uint_t a, uint_t b){ + Operations::Op op; + op.type = Operations::OpType::gate; + op.name = "cx"; + op.qubits = {a, b}; + return op; +} + +Operations::Op make_CZ(uint_t a, uint_t b){ + Operations::Op op; + op.type = Operations::OpType::gate; + op.name = "cz"; + op.qubits = {a, b}; + return op; +} + + +std::vector AGState::simplifying_unitary(){ + std::vector circuit; + //first we do "augmented gaussian elimination" on the circuit + //augmented here means that if we don't find a pivot in our column + //we will hadamard that qubit and try again + //in other words bring in a pivot from the z part if necessary + + size_t h = 0; + size_t k = 0; + while(h < this->num_stabilizers && k < this->num_qubits){ + std::pair poss_pivot = this->first_non_zero_in_col(k, h); + if(!poss_pivot.first){ + this->applyH(k); + circuit.push_back(make_H(k)); + poss_pivot = this->first_non_zero_in_col(k, h); + } + if(!poss_pivot.first){ + k += 1; + }else{ + size_t pivot = poss_pivot.second; //now known to exist + if(pivot != h){ + //swap rows h and pivot of the table + this->swap_rows(h,pivot); + } + for(size_t j = 0; j < this->num_stabilizers; j++){ + if((j != h) && this->table[j].X[k]){ + this->rowsum(j,h); + } + } + h += 1; + k += 1; + } + } + + //so now we have a reduced row echelon form with the X part of the table having full rank + + //we swap columns (using CX) to make the X part into a kxk identity followed by a "junk" block + + for(size_t r = 0; r < this->num_stabilizers; r++){ + if(!this->table[r].X[r]){ + size_t col = this->first_non_zero_in_row(r, 0).second; + + this->applyCX(r, col); + circuit.push_back(make_CX(r,col)); + this->applyCX(col, r); + circuit.push_back(make_CX(col,r)); + this->applyCX(r, col); + circuit.push_back(make_CX(r,col)); + } + } + + + //now we use CX to clear out the "junk" block at the end of the X table + for(size_t r = 0; r < this->num_stabilizers; r++){ + for(size_t col = this->num_stabilizers; col < this->num_qubits; col++){ + if(this->table[r].X[col]){ + this->applyCX(r, col); + circuit.push_back(make_CX(r,col)); + } + } + } + + //now we clear the leading diagonal of the z block + for(size_t r = 0; r < this->num_stabilizers; r++){ + if(this->table[r].Z[r]){ + this->applyS(r); + circuit.push_back(make_S(r)); + } + } + + //clear out the last k x (n-k) block of the z matrix + for(size_t col = this->num_stabilizers; col < this->num_qubits; col++){ + for(size_t r = 0; r < this->num_stabilizers; r++){ + if(this->table[r].Z[col]){ + this->applyCZ(r, col); + circuit.push_back(make_CZ(r, col)); + } + } + } + + //clear out the first k x k block of the Z matrix, using that it is symmetric + for(size_t col = 0; col < this->num_stabilizers; col++){ + for(size_t r = 0; r < col; r++){ + if(this->table[r].Z[col]){ + this->applyCZ(r, col); + circuit.push_back(make_CZ(r, col)); + } + } + } + + //fix the phases + for(size_t r = 0; r < this->num_stabilizers; r++){ + if(this->phases[r]){ + this->applyS(r); + circuit.push_back(make_S(r)); + this->applyS(r); + circuit.push_back(make_S(r)); + } + } + + //swap the identity matrix to the z part + for(size_t r = 0; r < this->num_stabilizers; r++){ + this->applyH(r); + circuit.push_back(make_H(r)); + } + + return circuit; +} + +std::pair AGState::find_x_on_q(size_t q, size_t a){ + for(size_t row = a; row < this->num_stabilizers-1; row++){ + if(this->table[row].X[q] && !this->table[row].Z[q]){ + return std::pair(true, row); + } + } + return std::pair(false, 0); +} + +std::pair AGState::find_y_on_q(size_t q, size_t a){ + for(size_t row = a; row < this->num_stabilizers-1; row++){ + if(this->table[row].X[q] && this->table[row].Z[q]){ + return std::pair(true, row); + } + } + return std::pair(false, 0); +} + +std::pair AGState::find_z_on_q(size_t q, size_t a){ + for(size_t row = a; row < this->num_stabilizers-1; row++){ + if(!this->table[row].X[q] && this->table[row].Z[q]){ + return std::pair(true, row); + } + } + return std::pair(false, 0); +} + +/* + * we test that if you ignore the first q+1 qubits whether that last (n-q-1) part of the last stabiliser in the table can be generated (up to phase) + * by the other stabilisers + * for use in the apply_constraints code + */ +bool AGState::independence_test(int q){ + //we basically do gaussian elimination + + size_t a = 0; + size_t b = q+1; + while(a < this->num_stabilizers-1 && b < this->num_qubits){ + std::pair x = this->find_x_on_q(b, a); + std::pair y = this->find_y_on_q(b, a); + std::pair z = this->find_z_on_q(b, a); + + if(y.first && x.first){ + this->rowsum(y.second, x.second); + z = y; + y = std::pair(false,0); + } + if(y.first && z.first){ + this->rowsum(y.second, z.second); + x = y; + y = std::pair(false,0); + } + + if(x.first){ + if(x.second != a){ + this->swap_rows(a, x.second); + } + if(z.first && (z.second == a)){ + z = x; + } + for(size_t j = 0; j < this->num_stabilizers; j++){ + if((j != a) && this->table[j].X[b]){ + this->rowsum(j,a); + } + } + a += 1; + } + if(y.first){ + if(y.second != a){ + this->swap_rows(a,y.second); + } + for(size_t j = 0; j < this->num_stabilizers; j++){ + if((j != a) && this->table[j].X[b] && this->table[j].Z[b]){ + this->rowsum(j,a); + } + } + a += 1; + } + if(z.first){ + if(z.second != a){ + this->swap_rows(a,z.second); + } + for(size_t j = 0; j < this->num_stabilizers; j++){ + if((j != a) && this->table[j].Z[b]){ + this->rowsum(j,a); + } + } + a += 1; + } + b += 1; + } + + for(size_t p = q+1; p < this->num_qubits; p++){ + if((this->table[this->num_stabilizers-1].X[p] == 1) || (this->table[this->num_stabilizers-1].Z[p] == 1)){ + return true; + } + } + return false; +} + + +/* + * Apply constraints arising from the fact that we measure the first w qubits and project the last t onto T gates + * In particular we kill stabilisers if qubits [0, w) get killed by taking the expectation value <0| P |0> + * and we kill stabilisers if qubits in [w, table->n-t) aren't the identity + * we do not remove any qubits from the table + * note that it is possible for the region a constraints to be "inconsistent" - this means that the probability we were calculating is zero + * in that case we return pair(false, 0) + * in other cases we return pair(true, v) where v is as defined in the paper + */ +std::pair AGState::apply_constraints(size_t w, size_t t){ + size_t v = 0; + + //first apply region a constraints (measurement) + for(size_t q=0; q < w; q++){ //iterate over all the measured qubits + std::pair y_stab = std::pair(false, 0); //store the index of the first stab we come to with both x and z = 1 on this qubit + std::pair x_stab = std::pair(false, 0); //store the index of the first stab we come to with x=1, z=0 + std::pair z_stab = std::pair(false, 0); //store the index of the first stab we come to with z=1, x=0 + + for(size_t s=0; s < this->num_stabilizers && (((!y_stab.first) + (!x_stab.first) + (!z_stab.first)) > 1); s++){//iterate over all stabilisers and find interesting stabilisers + + if(this->table[s].X[q] && this->table[s].Z[q]){ + y_stab.first = true; + y_stab.second = s; + } + if(this->table[s].X[q] && !this->table[s].Z[q]){ + x_stab.first = true; + x_stab.second = s; + } + if(!this->table[s].X[q] && this->table[s].Z[q]){ + z_stab.first = true; + z_stab.second = s; + } + } + //there are several cases here + //either a single z, a single x, a single y or we can generate the whole Pauli group on this qubit + + //case 1) we generate the whole group + //put things in standard form (first stab is x then z) + if((y_stab.first + x_stab.first + z_stab.first) >= 2){ //we have at least two of the set + if(!x_stab.first){//we don't have a generator for x alone, but we can make one + this->rowsum(y_stab.second, z_stab.second); + //now we have a z and an x but not a y + x_stab = y_stab; + y_stab = std::pair(false,0); + }else if(!z_stab.first){//we don't have a generator for z alone, but we can make one + this->rowsum(y_stab.second, x_stab.second); + //now we have a z and an x but not a y + z_stab = y_stab; + y_stab = std::pair(false,0); + } + } + + if(y_stab.first && x_stab.first && z_stab.first){ //we have all 3 + //ignore the y one if we have all 3 + y_stab = std::pair(false,0); + } + + //now the only possibilities are that we have an an x, y or z or both an x and a z + //zero everything else on this qubit + for(size_t s = 0; s < this->num_stabilizers; s++){ + if((!y_stab.first || s != y_stab.second) && (!x_stab.first || s != x_stab.second) && (!z_stab.first || s != z_stab.second)){ + if(this->table[s].X[q] && this->table[s].Z[q] && y_stab.first){ + this->rowsum(s, y_stab.second); + } + + if(this->table[s].X[q]){ + this->rowsum(s, x_stab.second); + } + + if(this->table[s].Z[q]){ + this->rowsum(s, z_stab.second); + } + } + } + + //case 1 - there is a generator which does not commute with Z_q + if(y_stab.first || x_stab.first){ + //we can't have both >= 0 + size_t non_commuting_generator = y_stab.first ? y_stab.second : x_stab.second; + //we delete the non-commuting guy + this->swap_rows(non_commuting_generator, this->num_stabilizers-1); + this->delete_last_row(); + }else{ + //case 2 - all generators commute with Z_q + //our generating set contains either Z_q or -Z_q + //we need to work out which one it is + + //swap our Z_q guy to the end + this->swap_rows(z_stab.second, this->num_stabilizers-1); + bool independent = this->independence_test(q); + + if(!independent){ + if(this->phases[this->num_stabilizers-1] == 0){ + // +Z_q + v += 1; + this->delete_last_row(); + }else{ + //our chosen measurement outcome is impossible + return std::pair(false,0); + } + }else{ + //if we get here there has been an error + //TODO decide if we're going to throw an exception or print an error message here + } + } + } + + //time to impose region b constraints + for(size_t q=w; q < this->num_qubits - t ; q++){ //iterate over all the non-measured non-magic qubits + std::pair y_stab = std::pair(false,0); //store the index of the first stab we come to with both x and z = 1 on this qubit + std::pair x_stab = std::pair(false,0); //store the index of the first stab we come to with x=1, z=0 + std::pair z_stab = std::pair(false,0); //store the index of the first stab we come to with z=1, x=0 + + for(size_t s=0; s < this->num_stabilizers && (((!y_stab.first) + (!x_stab.first) + (!z_stab.first)) > 1); s++){//iterate over all stabilisers and find interesting stabilisers + if(this->table[s].X[q] && this->table[s].Z[q]){ + y_stab.first = true; + y_stab.second = s; + } + if(this->table[s].X[q] && !this->table[s].Z[q]){ + x_stab.first = true; + x_stab.second = s; + } + if(!this->table[s].X[q] && this->table[s].Z[q]){ + z_stab.first = true; + z_stab.second = s; + } + } + + //there are several cases here + //either a single z, a single x, a single y or we can generate the whole Pauli group on this qubit + + //case 1) we generate the whole group + //put things in standard form (first stab is x then z) + if((y_stab.first + x_stab.first + z_stab.first) >= 2){ //we have at least two of the set + if(!x_stab.first){//we don't have a generator for x alone, but we can make one + this->rowsum(y_stab.second, z_stab.second); + //now we have a z and an x but not a y + x_stab = y_stab; + y_stab = std::pair(false,0); + }else if(!z_stab.first){//we don't have a generator for z alone, but we can make one + this->rowsum(y_stab.second, x_stab.second); + //now we have a z and an x but not a y + z_stab = y_stab; + y_stab = std::pair(false,0); + } + } + + if(y_stab.first && x_stab.first && z_stab.first){ //we have all 3 + //ignore the y one if we have all 3 + y_stab = std::pair(false,0); + } + + //now the only possibilities are that we have an x_and_z, an x a z or an x and a z + //zero everything else on this qubit + for(size_t s = 0; s < this->num_stabilizers; s++){ + if((!y_stab.first || s != y_stab.second) && (!x_stab.first || s != x_stab.second) && (!z_stab.first || s != z_stab.second)){ + if(this->table[s].X[q] && this->table[s].Z[q] && y_stab.first){ + this->rowsum(s, y_stab.second); + } + + if(this->table[s].X[q]){ + this->rowsum(s, x_stab.second); + } + + if(this->table[s].Z[q]){ + this->rowsum(s, z_stab.second); + } + } + } + + //now we just delete the non-identity guys on this qubit + uint_t num_to_delete = 0; + if(y_stab.first){ + //if we have a Y stab we don't have either of the others + this->swap_rows(y_stab.second, this->num_stabilizers-1); + num_to_delete += 1; + }else{ + if(x_stab.first){ + this->swap_rows(x_stab.second, this->num_stabilizers-1); + if(z_stab.first && (this->num_stabilizers - 1 == z_stab.second)){ + z_stab = x_stab; + } + num_to_delete += 1; + } + if(z_stab.first){ + this->swap_rows(z_stab.second, this->num_stabilizers-1-num_to_delete); + num_to_delete += 1; + } + } + + //delete the last num_to_delete rows + //TODO should we implement an erase method that works like the std::vector one so you can pass a range? + for(size_t deletes = 0; deletes < num_to_delete; deletes++){ + this->delete_last_row(); + } + } + + return std::pair(true, v); +} + + +/* + * Our magic states are equatorial + * so = 0 + * here we delete any stabilisers with a Z in the magic region + * which we assume now all the qubits + * we return the number of qubits which have identities on them in every generator after this deletion + */ +size_t AGState::apply_T_constraints(){ + size_t starting_rows = this->num_stabilizers; + size_t deleted_rows = 0; + for(size_t reps = 0; reps < starting_rows; reps++){ + for(size_t q=0; q < this->num_qubits; q++){ //iterate over all the magic qubits + std::pair y_stab = std::pair(false,0); //store the index of the first stab we come to with both x and z = 1 on this qubit + std::pair x_stab = std::pair(false,0); //store the index of the first stab we come to with x=1, z=0 + std::pair z_stab = std::pair(false,0); //store the index of the first stab we come to with z=1, x=0 + + for(size_t s=0; s < this->num_stabilizers && (((!y_stab.first) + (!x_stab.first) + (!z_stab.first)) > 1); s++){//iterate over all stabilisers and find interesting stabi lisers + if(this->table[s].X[q] && this->table[s].Z[q]){ + y_stab.first = true; + y_stab.second = s; + } + if(this->table[s].X[q] && !this->table[s].Z[q]){ + x_stab.first = true; + x_stab.second = s; + } + if(!this->table[s].X[q] && this->table[s].Z[q]){ + z_stab.first = true; + z_stab.second = s; + } + } + + //there are several cases here + //either a single z, a single x, a single y or we can generate the whole Pauli group on this qubit + + //case 1) we generate the whole group + //put things in standard form (first stab is x then z) + if((y_stab.first + x_stab.first + z_stab.first) >= 2){ //we have at least two of the set + if(!x_stab.first){//we don't have a generator for x alone, but we can make one + this->rowsum(y_stab.second, z_stab.second); + //now we have a z and an x but not a y + x_stab = y_stab; + y_stab = std::pair(false,0); + }else if(!z_stab.first){//we don't have a generator for z alone, but we can make one + this->rowsum(y_stab.second, x_stab.second); + //now we have a z and an x but not a y + z_stab = y_stab; + y_stab = std::pair(false,0); + } + } + + if(y_stab.first && x_stab.first && z_stab.first){ //we have all 3 + //ignore the y one if we have all 3 + y_stab = std::pair(false,0); + } + + if(z_stab.first && !x_stab.first){ + //kill all other z stuff on this qubit + for(size_t s = 0; s < this->num_stabilizers; s++){ + if((s != z_stab.second) && this->table[s].Z[q]){ + this->rowsum(s, z_stab.second); + } + } + //now delete the z guy + + if(z_stab.second != this->num_stabilizers-1){ + this->swap_rows(z_stab.second, this->num_stabilizers-1); + } + this->delete_last_row(); + deleted_rows += 1; + } + } + } + return deleted_rows; +} + +inline void to_json(json_t &js, const AGState &state) +{ + js["num_qubits"] = state.num_qubits; + js["num_stabilizers"] = state.num_stabilizers; + js["phases"] = state.phases; + js["magic_phases"] = state.magic_phases; +} + +} +} +#endif diff --git a/src/simulators/clifford_plus_phase/compute.hpp b/src/simulators/clifford_plus_phase/compute.hpp new file mode 100644 index 0000000000..5389a5e8e3 --- /dev/null +++ b/src/simulators/clifford_plus_phase/compute.hpp @@ -0,0 +1,565 @@ +#ifndef _aer_extended_stabilizer_compute_hpp +#define _aer_extended_stabilizer_compute_hpp + +#define _USE_MATH_DEFINES +#include + +#include +#include +#include + +#include "simulators/state.hpp" +#include "framework/json.hpp" +#include "framework/types.hpp" + +#include "ag_state.hpp" +#include "simulators/stabilizer/pauli.hpp" + +namespace AER{ +namespace CliffPhaseCompute{ + +// OpSet of supported instructions +const Operations::OpSet StateOpSet( + // Op types + {Operations::OpType::gate, Operations::OpType::save_specific_prob}, + // Gates + {"CX", "cx", "cz", "swap", "id", "x", "y", "z", "h", + "s", "sdg", "t","p", "rz", "u1"}, + // Snapshots + {} + ); + + + +enum class Gates { + id, x, y, z, h, s, sdg, sx, t, tdg, cx, cz, swap, p, rz, u1, +}; + + +using agstate_t = CliffPhaseCompute::AGState; + + +class State: public Base::State{ +public: + using BaseState = Base::State; + + State() : BaseState(StateOpSet) {}; + virtual ~State() = default; + + std::string name() const override {return "clifford_phase_compute";} + + //Apply a sequence of operations to the cicuit + template + void apply_ops(InputIterator first, InputIterator last, + ExperimentResult &result, + RngEngine &rng, + bool final_ops = false); + + virtual void apply_op(const Operations::Op &op, + ExperimentResult &result, + RngEngine &rng, + bool final_op = false) override; + + void apply_gate(const Operations::Op &op); + + void initialize_qreg(uint_t num_qubits) override; + void initialize_qreg(uint_t num_qubits, const agstate_t &state) override; + size_t required_memory_mb(uint_t num_qubits, const std::vector &ops) const override; + void apply_save_specific_prob(const Operations::Op &op, ExperimentResult &result); + double expval_pauli(const reg_t &qubits, const std::string& pauli); +private: + const static stringmap_t gateset_; + double compute_algorithm_all_phases_T(AGState &state); + double compute_algorithm_arbitrary_phases(AGState &state); + size_t num_code_qubits; //our AG state has code+magic qubits + double compute_probability(std::vector measured_qubits, std::vector outcomes); + template + uint_t count_magic_gates(InputIterator first, InputIterator last) const; +}; + + +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)) + {"t", Gates::t}, // T-gate (sqrt(S)) + {"rz", Gates::rz}, // Pauli-Z rotation gate + {"p", Gates::rz}, // Parameterized phase gate + {"u1", Gates::rz}, + {"tdg", Gates::tdg}, // Conjguate-transpose of T gate + // Two-qubit gates + {"CX", Gates::cx}, // Controlled-X gate (CNOT) + {"cx", Gates::cx}, // Controlled-X gate (CNOT) + {"CZ", Gates::cz}, // Controlled-Z gate + {"cz", Gates::cz}, // Controlled-Z gate + {"swap", Gates::swap}, // SWAP gate + // Three-qubit gates + }); + +template +void State::apply_ops(InputIterator first, InputIterator last, ExperimentResult &result, RngEngine &rng, bool final_ops){ + for(auto it = first; it != last; ++it){ + apply_op(*it, result, rng, final_ops); + } +} + +void State::apply_op(const Operations::Op &op, ExperimentResult &result, + RngEngine &rng, bool final_op) { + switch(op.type){ + case Operations::OpType::gate: + this->apply_gate(op); + break; + case Operations::OpType::save_specific_prob: + this->apply_save_specific_prob(op, result); + break; + default: + throw std::invalid_argument("Compute::State::invalid instruction \'" + op.name + "\'."); + } +} + +void State::apply_gate(const Operations::Op &op){ + auto it = gateset_.find(op.name); + if (it == gateset_.end()) + { + throw std::invalid_argument("Compute::State: Invalid gate operation \'" + +op.name + "\'."); + } + switch(it->second){ + case Gates::id: + break; + case Gates::x: + this->qreg_.applyX(op.qubits[0]); + break; + case Gates::y: + this->qreg_.applyY(op.qubits[0]); + break; + case Gates::z: + this->qreg_.applyZ(op.qubits[0]); + break; + case Gates::s: + this->qreg_.applyS(op.qubits[0]); + break; + case Gates::sdg: + this->qreg_.applyZ(op.qubits[0]); + this->qreg_.applyS(op.qubits[0]); + break; + case Gates::h: + this->qreg_.applyH(op.qubits[0]); + break; + case Gates::t: + this->qreg_.gadgetized_phase_gate(op.qubits[0], T_ANGLE); + break; + case Gates::tdg: + this->qreg_.gadgetized_phase_gate(op.qubits[0], -T_ANGLE); + break; + case Gates::rz: + this->qreg_.gadgetized_phase_gate(op.qubits[0], op.params[0].real()); + break; + case Gates::cx: + this->qreg_.applyCX(op.qubits[0],op.qubits[1]); + break; + case Gates::cz: + this->qreg_.applyCZ(op.qubits[0],op.qubits[1]); + break; + case Gates::swap: + this->qreg_.applySwap(op.qubits[0],op.qubits[1]); + break; + default: //u0 or Identity + break; + } +} + +void State::apply_save_specific_prob(const Operations::Op &op, ExperimentResult &result){ + double p = this->compute_probability(op.qubits, op.int_params); + save_data_average(result, op.string_params[0], p, op.type, op.save_type); +} + +// This function converts an unsigned binary number to reflected binary Gray code. +uint_t BinaryToGray(uint_t num) +{ + return num ^ (num >> 1); // The operator >> is shift right. The operator ^ is exclusive or. +} + + +double State::compute_algorithm_all_phases_T(AGState &state){ + if(state.num_stabilizers > 63){ + //we use an int_t == int_fast64_t variable to store our loop counter, this means we can deal with at most 63 stabilizers + //in the case of 63 stabilizers we will iterate over all length 63 bitstrings + //this is the most we can store in a signed 64 bit integer + //the integer has to be signed due to openMP's requirements + //realistically a computation with 63 stabilizers = 2^63 - 1 iterations is not going to terminate anyway so this restriction shouldn't matter + std::stringstream msg; + msg << "CliffPhaseCompute::State::compute_algorithm_all_phases_T called with " << state.num_stabilizers << " stabilizers. Maximum possible is 63."; + throw std::runtime_error(msg.str()); + } + + if(state.num_stabilizers == 0){ + return 1; + } + + uint_t num_iterations = (ONE << state.num_stabilizers); + + double acc = 0.; + Pauli::Pauli row(state.num_qubits); + unsigned char phase = 0; + + uint_t num_threads_ = BaseState::threads_; + if(num_threads_ > num_iterations){ + num_threads_ = num_iterations; + } + + uint_t chunk_size = num_iterations / num_threads_; + if(num_iterations % num_threads_ > 0){ + chunk_size += 1; + } +#pragma omp parallel for if(num_threads_ > 1) num_threads(num_threads_) \ + firstprivate(row, phase) shared(state) reduction(+:acc) + for(int_t thread = 0; thread < num_threads_; ++thread) { + uint_t start_value = thread * chunk_size; + uint_t end_value = (thread + 1) * chunk_size; + int_t mask_with_bits_to_flip = BinaryToGray(start_value); // note BinaryToGray(0u) == 0u + //set up the row and phase variables + for(int_t bit = 0; bit < state.num_stabilizers; bit++){ + if((mask_with_bits_to_flip >> bit) & ONE){ + phase += state.phases[bit] + (Pauli::Pauli::phase_exponent(row, state.table[bit]) / 2); + phase %= 2; + row += state.table[bit]; + } + } + + size_t XCount = 0; + size_t YCount = 0; + size_t ZCount = 0; + + for(size_t j = 0; j < state.num_qubits; j++){ + if(row.X[j] && !row.Z[j]){ + XCount += 1; + } + if(!row.X[j] && row.Z[j]){ + ZCount += 1; + break; + } + if(row.X[j] && row.Z[j]){ + YCount += 1; + } + } + + if(ZCount == 0){ + if(((phase + YCount) % 2) == 0){ + acc += powl(1./2., (XCount + YCount)/2.);; + }else{ + acc -= powl(1./2., (XCount + YCount)/2.);; + } + } + + + for(uint_t mask = start_value+1; mask < end_value && mask < num_iterations; ++mask){ + uint_t mask_with_bit_to_flip = BinaryToGray((uint_t)mask) ^ BinaryToGray((uint_t)(mask - 1)); + uint_t bit_to_flip = 0; + for(size_t bit = 0; bit < state.num_stabilizers; bit++){ + if((mask_with_bit_to_flip >> bit) & ONE){ + bit_to_flip = bit; + break; + } + } + phase += state.phases[bit_to_flip] + (Pauli::Pauli::phase_exponent(row, state.table[bit_to_flip]) / 2); + phase %= 2; + row += state.table[bit_to_flip]; + + size_t XCount = 0; + size_t YCount = 0; + size_t ZCount = 0; + + for(size_t j = 0; j < state.num_qubits; j++){ + if(row.X[j] && !row.Z[j]){ + XCount += 1; + } + if(!row.X[j] && row.Z[j]){ + ZCount += 1; + break; + } + if(row.X[j] && row.Z[j]){ + YCount += 1; + } + } + + if(ZCount == 0){ + if(((phase + YCount) % 2) == 0){ + acc += powl(1./2., (XCount + YCount)/2.);; + }else{ + acc -= powl(1./2., (XCount + YCount)/2.);; + } + } + } + } + + return acc; +} + +double State::compute_algorithm_arbitrary_phases(AGState &state){ + if(state.num_stabilizers > 63){ + //we use an int_t == int_fast64_t variable to store our loop counter, this means we can deal with at most 63 stabilizers + //in the case of 63 stabilizers we will iterate over all length 63 bitstrings + //this is the most we can store in a signed 64 bit integer + //the integer has to be signed due to openMP's requirements + //realistically a computation with 63 stabilizers = 2^63 - 1 iterations is not going to terminate anyway so this restriction shouldn't matter + std::stringstream msg; + msg << "CliffPhaseCompute::State::compute_algorithm_all_phases_T called with " << state.num_stabilizers << " stabilizers. Maximum possible is 63."; + throw std::runtime_error(msg.str()); + } + + if(state.num_stabilizers == 0){ + return 1; + } + + + std::vector x_values; + std::vector y_values; + + for(double phase : state.magic_phases){ + x_values.push_back(cos(phase)); + y_values.push_back(-sin(phase)); + } + + uint_t num_iterations = (ONE << state.num_stabilizers); + + double acc = 0.; + + uint_t num_threads_ = BaseState::threads_; + if(num_threads_ > num_iterations){ + num_threads_ = num_iterations; + } + uint_t chunk_size = num_iterations / num_threads_; + if((num_iterations % num_threads_) > 0){ + chunk_size += 1; + } + //std::cout << state.num_stabilizers << ", " << num_iterations << ", " << chunk_size << ", " << num_threads_ << std::endl; + +#pragma omp parallel for if(num_threads_ > 1) num_threads(num_threads_) shared(state) reduction(+:acc) + for(int_t thread = 0; thread < num_threads_; ++thread) { + uint_t start_value = thread * chunk_size; + uint_t end_value = (thread + 1) * chunk_size; + + //we do the first iteration of each thread separately because it has to initialise the row and phase variables + Pauli::Pauli row(state.num_qubits); + unsigned char phase = 0; + //if(start_value >= num_iterations){ + // std::cout << "wtf " << thread << " " << start_value << " " << chunk_size << std::endl; + // + uint_t mask_with_bits_to_flip = BinaryToGray(start_value); // note BinaryToGray(0u) == 0u + for(int_t bit = 0; bit < state.num_stabilizers; bit++){ + if((mask_with_bits_to_flip >> bit) & ONE){ + phase += state.phases[bit] + (Pauli::Pauli::phase_exponent(row, state.table[bit]) / 2); + phase %= 2; + row += state.table[bit]; + } + } + double prod = (phase == 0) ? 1. : -1.; + int ZCount = 0; + for(size_t j = 0; j < state.num_qubits; j++){ + if(row.X[j] && !row.Z[j]){ + prod *= x_values[j]; + } + if(!row.X[j] && row.Z[j]){ + ZCount += 1; + break; + } + if(row.X[j] && row.Z[j]){ + prod *= y_values[j]; + } + } + + if(ZCount == 0){ + acc += prod; + } + + for(uint_t mask = start_value+1; (mask < end_value) && (mask < num_iterations); ++mask){ + uint_t mask_with_bit_to_flip = BinaryToGray(mask) ^ BinaryToGray(mask - 1); + uint_t bit_to_flip = 0; + for(size_t bit = 0; bit < state.num_stabilizers; bit++){ + if((mask_with_bit_to_flip >> bit) & ONE){ + bit_to_flip = bit; + break; + } + } + phase += state.phases[bit_to_flip] + (Pauli::Pauli::phase_exponent(row, state.table[bit_to_flip]) / 2); + phase %= 2; + row += state.table[bit_to_flip]; + + double prod = (phase == 0) ? 1. : -1.; + int ZCount = 0; + for(size_t j = 0; j < state.num_qubits; j++){ + if(row.X[j] && !row.Z[j]){ + prod *= x_values[j]; + } + if(!row.X[j] && row.Z[j]){ + ZCount += 1; + break; + } + if(row.X[j] && row.Z[j]){ + prod *= y_values[j]; + } + } + + if(ZCount == 0){ + acc += prod; + } + } + } + + return acc; +} + + +double State::compute_probability(std::vector measured_qubits, std::vector outcomes){ + AGState copied_ag(this->qreg_); //copy constructor TODO check this + + //first reorder things so the first w qubits are measured + std::vector measured_qubits_sorted(measured_qubits); + + std::vector qubit_indexes; + for(uint_t i = 0; i < this->qreg_.num_qubits; i++){ + qubit_indexes.push_back(i); + } + std::sort(measured_qubits_sorted.begin(), measured_qubits_sorted.end()); + for(uint_t i = 0; i < measured_qubits.size(); i++){ + uint_t w = measured_qubits_sorted[i]; + uint_t idx1 = 0; + uint_t idx2 = 0; + for(uint_t j = 0; j < qubit_indexes.size(); j++){ + if(qubit_indexes[j] == w){ + idx1 = j; + break; + } + } + for(uint_t j = 0; j < measured_qubits.size(); j++){ + if(measured_qubits[j] == w){ + idx2 = j; + break; + } + } + + if(idx1 != idx2){ + copied_ag.applySwap(idx1, idx2); + } + } + + //now all the measured qubits are at the start and the magic qubits are at the end + //from this point on we will assume we're looking for the measurement outcome 0 on all measured qubits + //so apply X gates to measured qubits where we're looking for outcome 1 to correct this + for(uint_t i = 0; i < outcomes.size(); i++){ + if(outcomes[i] == 1){ + copied_ag.applyX(i); + } + } + + uint_t w = measured_qubits.size(); + uint_t t = copied_ag.magic_phases.size(); + + //now all the measured qubits are at the start and the magic qubits are at the end + + std::pair v_pair = copied_ag.apply_constraints(w, t); + + if(!v_pair.first){ + return 0.; + } + uint_t v = v_pair.second; + + //at this point we can delete all the non-magic qubits + for(uint_t q = 0; q < t; q++){ + copied_ag.applySwap(q, q+(copied_ag.num_qubits - t)); + } + + for(uint_t s = 0; s < copied_ag.num_stabilizers; s++){ + copied_ag.table[s].X.resize(t); + copied_ag.table[s].Z.resize(t); + } + + copied_ag.num_qubits = t; + copied_ag.apply_T_constraints(); + copied_ag.delete_identity_magic_qubits(); + + //we can make the compute algorithm much faster if all of our non-Clifford gates are in fact T gates (pi/4 rotations) + bool all_phases_are_T = true; + for(uint_t i = 0; i < copied_ag.num_qubits; i++){ + if(fabs(copied_ag.magic_phases[i] - T_ANGLE) > AG_CHOP_THRESHOLD){ + all_phases_are_T = false; + break; + } + } + + if(copied_ag.num_qubits == 0){ + return pow(2., ((double)v) - ((double)w)); + } + + if(all_phases_are_T){ + return compute_algorithm_all_phases_T(copied_ag) * powl(2., ((double)v) - ((double)w)); + } else { + return compute_algorithm_arbitrary_phases(copied_ag) * powl(2., ((double)v) - ((double)w)); + } +} + +void State::initialize_qreg(uint_t num_qubits) { + this->qreg_.initialize(num_qubits); + this->num_code_qubits = num_qubits; +} +void State::initialize_qreg(uint_t num_qubits, const agstate_t &state) { + if(BaseState::qreg_.num_qubits != num_qubits){ + throw std::invalid_argument("CH::State::initialize: initial state does not match qubit number."); + } + BaseState::qreg_ = state; +} + +template +uint_t State::count_magic_gates(InputIterator first, InputIterator last) const { + uint_t count = 0; + for (auto op = first; op != last; op++) + { + auto it = gateset_.find(op->name); + if (it != gateset_.end()) + { + if(it->second == Gates::t || it->second == Gates::rz || it->second == Gates::tdg){ + count += 1; + } + } + } + return count; +} + + +size_t State::required_memory_mb(uint_t num_qubits, const std::vector &ops) const { + uint_t t = count_magic_gates(ops.cbegin(), ops.cend()); + + uint_t total_qubits = num_qubits + t; + + // we store a stabilizer tableau with n+t stabilizers on n+t qubits + // each stabilizer (each "row" in the table) is a Pauli::Pauli on n+t qubits + // each Pauli:: Pauli stores two BV::BinaryVectors of length n+t and in addition needs one byte of space we use to store the phase + // each BinaryVector stores 8 * ceil((n+t)/64) bytes + //note that this is an upper bound, in reality the compress algorithm will give us a state with fewer qubits and stabilizers + + //in addition each working thread gets a Pauli (two binary vectors and a phase) of space to do its work in + uint_t bv_size_bytes = 8*((num_qubits+t)/64 + ((((num_qubits+t) % 64)) ? 1 : 0)); //add an extra 8 byte int that we partially use if there is a remainder + uint_t pauli_size_bytes = 2*bv_size_bytes + 1; + + // State::compute_probability copies the stabilizer tableau before operating on it so double the memory required for that + size_t mb = (pauli_size_bytes * (2*(num_qubits+t)+BaseState::threads_) + t*sizeof(double))/(1<<20); + + return mb; +} + + +double State::expval_pauli(const reg_t &qubits, const std::string& pauli) { + return 0; //TODO fix this +} + +} //close namespace CliffPhaseCompute +} //close namespace AER + +#endif diff --git a/src/simulators/stabilizer/binary_vector.hpp b/src/simulators/stabilizer/binary_vector.hpp index a982097c17..d8b2b7775d 100644 --- a/src/simulators/stabilizer/binary_vector.hpp +++ b/src/simulators/stabilizer/binary_vector.hpp @@ -46,6 +46,10 @@ class BinaryVector { void setLength(uint64_t length); + //preserves the vector up to the new length + //and pads with zeros if necessary + void resize(uint64_t length); + void setVector(std::string); void setValue(bool value, uint64_t pos); @@ -53,7 +57,8 @@ class BinaryVector { void set1(uint64_t pos) { setValue(ONE_, pos); }; void flipAt(uint64_t pos); - + void xorAt(bool value, uint64_t pos); + BinaryVector &operator+=(const BinaryVector &rhs); bool operator[](const uint64_t pos) const; @@ -180,6 +185,21 @@ void BinaryVector::setLength(uint64_t length) { m_data.assign((length - 1) / BLOCK_SIZE + 1, ZERO_); } +void BinaryVector::resize(uint64_t new_length) { + if(new_length == 0){ + m_data.resize(0); + } else { + m_data.resize((new_length - 1) / BLOCK_SIZE + 1, ZERO_); + } + + //zero the rest of the last block if necessary + if((new_length < m_length) && (new_length % BLOCK_SIZE) > 0){ + for(size_t i = (new_length % BLOCK_SIZE) ; i < BLOCK_SIZE; i++){ + m_data[m_data.size() - 1] &= ~(ONE_ << i); + } + } + m_length = new_length; +} void BinaryVector::setValue(bool value, uint64_t pos) { auto q = pos / BLOCK_SIZE; @@ -197,6 +217,11 @@ void BinaryVector::flipAt(const uint64_t pos) { m_data[q] ^= (ONE_ << r); } +void BinaryVector::xorAt(bool value, uint64_t pos) { + auto q = pos / BLOCK_SIZE; + auto r = pos % BLOCK_SIZE; + m_data[q] ^= ((value & ONE_) << r); +} BinaryVector &BinaryVector::operator+=(const BinaryVector &rhs) { const auto size = m_data.size(); @@ -297,4 +322,4 @@ std::vector BinaryVector::nonzeroIndices() const { //------------------------------------------------------------------------------ } // end namespace BV //------------------------------------------------------------------------------ -#endif \ No newline at end of file +#endif diff --git a/src/simulators/stabilizer/pauli.hpp b/src/simulators/stabilizer/pauli.hpp index ecf0ec5673..0bb10c34f1 100644 --- a/src/simulators/stabilizer/pauli.hpp +++ b/src/simulators/stabilizer/pauli.hpp @@ -46,6 +46,8 @@ class Pauli { // 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); + + Pauli &operator+=(const Pauli &rhs); }; @@ -114,6 +116,12 @@ int8_t Pauli::phase_exponent(const Pauli& pauli1, const Pauli& pauli2) { return exponent; } +Pauli &Pauli::operator+=(const Pauli &rhs) { + this->X += rhs.X; + this->Z += rhs.Z; + return (*this); +} + //------------------------------------------------------------------------------ } // end namespace Pauli //------------------------------------------------------------------------------