diff --git a/examples/quantum/gate/CMakeLists.txt b/examples/quantum/gate/CMakeLists.txt index f0fb15b997030b0acac4447ff4b0749cb7152254..e014e7d8f07bfc4d3d36f1d01f7cdeaecd515b69 100644 --- a/examples/quantum/gate/CMakeLists.txt +++ b/examples/quantum/gate/CMakeLists.txt @@ -29,7 +29,7 @@ # #**********************************************************************************/ include_directories(${CMAKE_SOURCE_DIR}/quantum/gate/accelerators) -include_directories(${CMAKE_SOURCE_DIR}/quantum/gate/accelerators/eigenaccelerator) +include_directories(${CMAKE_SOURCE_DIR}/quantum/gate/accelerators/firetensoraccelerator) include_directories(${CMAKE_SOURCE_DIR}/tpls/common/tpls/fire/tpls/eigen) include_directories(${CMAKE_SOURCE_DIR}/tpls/common/tpls/fire/tensors) include_directories(${CMAKE_SOURCE_DIR}/tpls/common/tpls/fire/tensors/impl) diff --git a/examples/quantum/gate/teleport_scaffold.cpp b/examples/quantum/gate/teleport_scaffold.cpp index 59bfb6bdbb765c8d31c3846edba1a27c9991b1f0..1b89f0381f6ffc07712fdf8750238d840bd5e4f2 100644 --- a/examples/quantum/gate/teleport_scaffold.cpp +++ b/examples/quantum/gate/teleport_scaffold.cpp @@ -29,7 +29,7 @@ * **********************************************************************************/ #include "XACC.hpp" -#include "EigenAccelerator.hpp" +#include "FireTensorAccelerator.hpp" // Quantum Kernel executing teleportation of // qubit state to another. @@ -51,7 +51,7 @@ const std::string src("__qpu__ teleport (qbit qreg) {\n" int main (int argc, char** argv) { // Create a convenient alias... - using Simple6QubitAcc = xacc::quantum::EigenAccelerator<6>; + using Simple6QubitAcc = xacc::quantum::FireTensorAccelerator<6>; // Create a reference to the 3 qubit simulation Accelerator auto qpu = std::make_shared<Simple6QubitAcc>(); diff --git a/quantum/gate/accelerators/SimulatedQubits.hpp b/quantum/gate/accelerators/SimulatedQubits.hpp index b93e9967fdd5b5abbb133c6f8e086974cd180fab..d5e8b73759afe1f101b3a348fa2a5cdc28e151b7 100644 --- a/quantum/gate/accelerators/SimulatedQubits.hpp +++ b/quantum/gate/accelerators/SimulatedQubits.hpp @@ -3,35 +3,55 @@ #include "Accelerator.hpp" #include <complex> -#include <Eigen/Dense> +#include "Tensor.hpp" namespace xacc { namespace quantum { -using QubitState = Eigen::VectorXcd; +using QubitState = fire::Tensor<1>; +using IndexPair = std::pair<int,int>; /** + * SimulatedQubits is an AcceleratorBuffer that + * models simulated qubits. As such, it keeps track of the + * state of the qubit buffer using a Rank 1 Fire Tensor. * + * It provides an interface for applying unitary operations on the + * qubit buffer state. */ template<const int TotalNumberOfQubits> class SimulatedQubits: public AcceleratorBuffer { protected: + /** + * The qubit buffer state. + */ QubitState bufferState; public: + /** + * The Constructor, creates a state of size + * 2^TotalNumberOfQubits. + * + * @param str + */ SimulatedQubits(const std::string& str) : - AcceleratorBuffer(str), bufferState((int) std::pow(2, TotalNumberOfQubits)) { + AcceleratorBuffer(str), bufferState( + (int) std::pow(2, TotalNumberOfQubits)) { // Initialize to |000...000> state - bufferState.setZero(); bufferState(0) = 1.0; } + /** + * The Constructor, creates a state with given size + * N. + * @param str + * @param N + */ SimulatedQubits(const std::string& str, const int N) : AcceleratorBuffer(str, N), bufferState((int) std::pow(2, N)) { - bufferState.setZero(); bufferState(0) = 1.0; } @@ -41,30 +61,62 @@ public: (int) std::pow(2, sizeof...(indices) + 1)) { } - void applyUnitary(Eigen::MatrixXcd& U) { - bufferState = U * bufferState; + /** + * Apply the given unitary matrix to this qubit + * buffer state. + * + * @param U + */ + void applyUnitary(fire::Tensor<2>& U) { + assert( + U.dimension(0) == bufferState.dimension(0) + && U.dimension(1) == bufferState.dimension(0)); + std::array<IndexPair, 1> contractionIndices; + contractionIndices[0] = std::make_pair(1, 0); + bufferState = U.contract(bufferState, contractionIndices); } + /** + * Normalize the state. + */ + void normalize() { + double sum = 0.0; + for (int i = 0; i < bufferState.dimension(0); i++) + sum += bufferState(i) * bufferState(i); + for (int i = 0; i < bufferState.dimension(0); i++) + bufferState(i) = bufferState(i) / std::sqrt(sum); + } + + /** + * Return the current state + * + * @return + */ QubitState& getState() { return bufferState; } + /** + * Set the state. + * @param st + */ void setState(QubitState& st) { bufferState = st; } + /** + * Print the state to the provided output stream. + * + * @param stream + */ void printBufferState(std::ostream& stream) { - for (int i = 0; i < bufferState.rows(); i++) { + for (int i = 0; i < bufferState.dimension(0); i++) { stream << std::bitset<TotalNumberOfQubits>(i).to_string().substr( size(), TotalNumberOfQubits) << " -> " << bufferState(i) << "\n"; } } - - void mapBufferStateToSystemState() { - - } }; } } diff --git a/quantum/gate/accelerators/eigenaccelerator/EigenAccelerator.hpp b/quantum/gate/accelerators/firetensoraccelerator/FireTensorAccelerator.hpp similarity index 65% rename from quantum/gate/accelerators/eigenaccelerator/EigenAccelerator.hpp rename to quantum/gate/accelerators/firetensoraccelerator/FireTensorAccelerator.hpp index 488391dad576d8b18808849441edd9e8b6ca968e..2f03bd9d0ff7464a8e5421b8665cdade8cedd3f8 100644 --- a/quantum/gate/accelerators/eigenaccelerator/EigenAccelerator.hpp +++ b/quantum/gate/accelerators/firetensoraccelerator/FireTensorAccelerator.hpp @@ -35,7 +35,6 @@ #include "QasmToGraph.hpp" #include "GraphIR.hpp" #include "SimulatedQubits.hpp" -#include <unsupported/Eigen/KroneckerProduct> #include <random> namespace xacc { @@ -46,44 +45,52 @@ using GraphType = qci::common::Graph<CircuitNode>; using QuantumGraphIR = xacc::GraphIR<GraphType>; /** - * + * The FireTensorAccelerator is an XACC Accelerator that simulates + * gate based quantum computing circuits. It models the QPUGate Accelerator + * with SimulatedQubit AcceleratorBuffer. It relies on the Fire Scientific Computing + * Framework's tensor module to model a specific set of quantum gates. It uses these + * tensors to build up the unitary matrix described by the circuit. */ template<const int NQubits> -class EigenAccelerator : virtual public QPUGate<SimulatedQubits<NQubits>> { +class FireTensorAccelerator : virtual public QPUGate<SimulatedQubits<NQubits>> { public: /** * The constructor, create tensor gates */ - EigenAccelerator() { - Eigen::MatrixXcd h(2,2), cnot(4,4), I(2,2), x(2,2), p0(2,2), p1(2,2), z(2,2); - h << 1.0/sqrt2,1.0/sqrt2, 1.0/sqrt2,-1.0/sqrt2; - cnot << 1, 0, 0, 0,0, 1, 0, 0,0, 0, 0, 1, 0, 0, 1, 0; - x << 0, 1, 1, 0; - I << 1,0,0,1; - p0 << 1, 0, 0, 0; - p1 << 0, 0, 0, 1; - z << 1, 0, 0, -1; - gates.insert(std::map<std::string, Eigen::MatrixXcd>::value_type("h",h)); - gates.insert(std::map<std::string, Eigen::MatrixXcd>::value_type("cnot",cnot)); - gates.insert(std::map<std::string, Eigen::MatrixXcd>::value_type("I",I)); - gates.insert(std::map<std::string, Eigen::MatrixXcd>::value_type("x",x)); - gates.insert(std::map<std::string, Eigen::MatrixXcd>::value_type("p0",p0)); - gates.insert(std::map<std::string, Eigen::MatrixXcd>::value_type("p1",p1)); - gates.insert(std::map<std::string, Eigen::MatrixXcd>::value_type("z",z)); + FireTensorAccelerator() { + fire::Tensor<2> h(2,2), cnot(4,4), I(2,2), x(2,2), p0(2,2), p1(2,2), z(2,2); + h.setValues({{1.0/sqrt2, 1.0/sqrt2},{1.0/sqrt2,-1.0/sqrt2}}); + cnot.setValues({{1,0,0,0},{0,0,1,0},{0,0,0,1},{0,0,1,0}}); + x.setValues({{0, 1},{1, 0}}); + I.setValues({{1,0},{0,1}}); + p0.setValues({{1,0},{0,0}}); + p1.setValues({{0,0},{0,1}}); + z.setValues({{1,0},{0,-1}}); + gates.insert(std::make_pair("h",h)); + gates.insert(std::make_pair("cont",cnot)); + gates.insert(std::make_pair("x",x)); + gates.insert(std::make_pair("I",I)); + gates.insert(std::make_pair("p0",p0)); + gates.insert(std::make_pair("p1",p1)); + gates.insert(std::make_pair("z",z)); } /** + * Execute the simulation. Requires both a valid SimulatedQubits buffer and + * Graph IR instance modeling the quantum circuit. * * @param ir */ virtual void execute(const std::string& bufferId, const std::shared_ptr<xacc::IR> ir) { + // Get the requested qubit buffer auto qubits = this->allocatedBuffers[bufferId]; if (!qubits) { QCIError("Invalid buffer id. Could not get qubit buffer."); } + // Set the size int nQubits = qubits->size(); // Cast to a GraphIR, if we can... @@ -92,23 +99,21 @@ public: QCIError("Invalid IR - this Accelerator only accepts GraphIR<Graph<CircuitNode>>."); } - // Get the Graph + // Get the Graph and related info std::vector<CircuitNode> gateOperations; std::map<int, int> qubitIdToMeasuredResult; auto graph = graphir->getGraph(); int nNodes = graph.order(), layer = 1; int finalLayer = graph.getVertexProperty<1>(nNodes - 1); + // Get a vector of all gates for (int i = 0; i < nNodes; i++) { CircuitNode n; n.properties = graph.getVertexProperties(i); gateOperations.emplace_back(n); } -// std::cout << "Initial State:\n"; -// qubits->printState(std::cout); -// std::cout << "\n"; - + // Loop over all gates in the circuit for (auto gate : gateOperations) { // Skip disabled gates... @@ -117,16 +122,17 @@ public: } // Create a list of nQubits Identity gates - std::vector<Eigen::MatrixXcd> productList(nQubits); + std::vector<fire::Tensor<2>> productList; for (int i = 0; i < nQubits; i++) { - productList[i] = gates["I"]; + productList.push_back(gates.at("I")); } - // Create a local U gate - Eigen::MatrixXcd localU; + // Create a local U gate, initialized to identity + fire::Tensor<2> localU = gates.at("I"); // Get the current gate anme auto gateName = std::get<0>(gate.properties); + // Get the qubits we're acting on... auto actingQubits = std::get<3>(gate.properties); @@ -153,26 +159,28 @@ public: if (std::get<0>(gateOperations[i].properties) == "FinalState") { break; } -// std::cout << "Enabling " << graph.getVertexProperty<0>(i) << "\n"; std::get<4>(gateOperations[i].properties) = true; } } } else if (gateName == "measure") { - // get rho - auto rho = qubits->getState() * qubits->getState().transpose(); + // get rho - outer product of state with itself + auto rho = qubits->getState() * qubits->getState(); - productList[actingQubits[0]] = gates["p0"]; // Create a total unitary for this layer of the circuit - auto temp = productList[0]; + productList.at(actingQubits[0]) = gates.at("p0"); + auto Pi0 = productList.at(0); for (int i = 1; i < productList.size(); i++) { - temp = kroneckerProduct(temp, productList[i]).eval(); + Pi0 = Pi0.kronProd(productList.at(i)); } // Get probability qubit is a 0 - auto temp2 = temp * rho; - auto probZero = temp2.trace(); + double probZero = 0.0; + std::array<IndexPair, 1> contractionIndices; + contractionIndices[0] = std::make_pair(1, 0); + auto Prob0 = Pi0.contract(rho, contractionIndices); + for (int i = 0; i < Prob0.dimension(0); i++) probZero += Prob0(i,i); // Make the measurement random... std::random_device rd; @@ -182,23 +190,21 @@ public: auto val = dist(mt); if (val < std::real(probZero)) { result = 0; - Eigen::VectorXcd newState = (temp * qubits->getState()); - newState.normalize(); - qubits->setState(newState); + qubits->applyUnitary(Pi0); + qubits->normalize(); } else { result = 1; - productList[actingQubits[0]] = gates["p1"]; + productList.at(actingQubits[0]) = gates.at("p1"); // Create a total unitary for this layer of the circuit - auto temp = productList[0]; + auto Pi1 = productList.at(0); for (int i = 1; i < productList.size(); i++) { - temp = kroneckerProduct(temp, productList[i]).eval(); + Pi1 = Pi1.kronProd(productList.at(i)); } - Eigen::VectorXcd newState = (temp * qubits->getState()); - newState.normalize(); - qubits->setState(newState); + qubits->applyUnitary(Pi1); + qubits->normalize(); } -// std::cout << "Measured qubit " << actingQubits[0] << " to be a " << result << ": prob was " << probZero << "\n"; + // Record the measurement result qubitIdToMeasuredResult.insert(std::make_pair(actingQubits[0], result)); } else { @@ -209,65 +215,62 @@ public: // If this is a one qubit gate, just replace // the currect I in the list with the gate - productList[actingQubits[0]] = gates[gateName]; + productList.at(actingQubits[0]) = gates.at(gateName); // Create a total unitary for this layer of the circuit - localU = productList[0]; + localU = productList.at(0); for (int i = 1; i < productList.size(); i++) { - localU = - kroneckerProduct(localU, productList[i]).eval(); + localU = localU.kronProd(productList.at(i)); } } else if (actingQubits.size() == 2) { // If this is a 2 qubit gate, then we need t // to construct Kron(P0, I, ..., I) + Kron(P1, I, ..., Gate, ..., I) - productList[actingQubits[0]] = gates["p0"]; - localU = productList[0]; + productList.at(actingQubits[0]) = gates.at("p0"); + localU = productList.at(0); for (int i = 1; i < productList.size(); i++) { - localU = - kroneckerProduct(localU, productList[i]).eval(); + localU = localU.kronProd(productList.at(i)); } // Now create the second term in the sum - productList[actingQubits[0]] = gates["p1"]; - productList[actingQubits[1]] = gates[gateName == "cnot" ? "x" : "I"]; - auto temp = productList[0]; + productList.at(actingQubits[0]) = gates.at("p1"); + productList.at(actingQubits[1]) = gates.at(gateName == "cnot" ? "x" : "I"); + auto temp = productList.at(0); for (int i = 1; i < productList.size(); i++) { - temp = - kroneckerProduct(temp, productList[i]).eval(); + temp = temp.kronProd(productList.at(i)); } // Sum them up - localU += temp; + localU = localU + temp; } else { QCIError("Can only simulate one and two qubit gates."); } // Make sure that localU is the correct size - assert(localU.rows() == std::pow(2, nQubits) && localU.cols() == std::pow(2,nQubits)); + assert( + localU.dimension(0) == std::pow(2, nQubits) + && localU.dimension(1) + == std::pow(2, nQubits)); + // Appy the unitary and update th state qubits->applyUnitary(localU); -// std::cout << "Current State after " << gateName << ":\n"; -// qubits->printState(std::cout); -// std::cout << "\n" << localU << "\n"; } } } - - // Map buffer state to accelerator system state - } /** * The destructor */ - virtual ~EigenAccelerator() { - } + virtual ~FireTensorAccelerator() {} protected: - std::map<std::string, Eigen::MatrixXcd> gates; + /** + * Mapping of gate names to actual gate matrices. + */ + std::map<std::string, fire::Tensor<2>> gates; }; } } diff --git a/quantum/gate/accelerators/tests/EigenAcceleratorTester.cpp b/quantum/gate/accelerators/tests/FireTensorAcceleratorTester.hpp similarity index 98% rename from quantum/gate/accelerators/tests/EigenAcceleratorTester.cpp rename to quantum/gate/accelerators/tests/FireTensorAcceleratorTester.hpp index 0c8dbdfc53548ca7c2c8ab1a6a7d02424866cc4e..37211f99054cae91504bd1f7515a766098a4fa46 100644 --- a/quantum/gate/accelerators/tests/EigenAcceleratorTester.cpp +++ b/quantum/gate/accelerators/tests/FireTensorAcceleratorTester.hpp @@ -33,14 +33,14 @@ #include <memory> #include <boost/test/included/unit_test.hpp> -#include "EigenAccelerator.hpp" +#include "FireTensorAccelerator.hpp" #include "GraphIR.hpp" using namespace xacc::quantum; BOOST_AUTO_TEST_CASE(checkConstruction) { - EigenAccelerator<3> acc; + FireTensorAccelerator<3> acc; // auto qreg = acc.allocate("qreg"); diff --git a/xacc/accelerator/Accelerator.hpp b/xacc/accelerator/Accelerator.hpp index 73a1e34300d3985e0721b6193344b4534dfee54e..562100c52880a44f850b7495252ae5c4a1b359f2 100644 --- a/xacc/accelerator/Accelerator.hpp +++ b/xacc/accelerator/Accelerator.hpp @@ -55,14 +55,6 @@ struct is_valid_bitstype<T, decltype(std::declval<T>().N, void())> : std::true_t enum AcceleratorType { qpu_gate, qpu_aqc, npu }; /** - * The AcceleratorBits class provides a common - * base class for allocating accelerator-specific - * bit resources (for example, qubits). It takes an - * integer template parameter at construction that indicates - * the number of bits this AcceleratorBits models. - * - * Derived Accelerators should define a subclass of this - * class that models the hardware. * * @author Alex McCaskey */ @@ -91,39 +83,6 @@ protected: int bufferSize = 0; std::string bufferId; }; -//template<const int Number> -//class AcceleratorBits { -//public: -// /** -// * Reference to the number of bits -// */ -// static constexpr int N = Number; -// -// /** -// * Return the current state of the bits -// * @return -// */ -// virtual std::bitset<(size_t) Number> measure() { -// return bits; -// } -// -// template<typename... SubBits> -// auto allocateSubset( -// SubBits ... subset) -> -// decltype(std::shared_ptr<AcceleratorBits<sizeof...(SubBits)>>()) { -// } -// virtual ~AcceleratorBits() {} -// -//protected: -// -// /** -// * The bits themselves -// */ -// std::bitset<(size_t)Number> bits; -// -// std::vector<int> activeBits; -// -//}; class IAccelerator : public qci::common::QCIObject { public: