Implemented EigenAccelerator to solve teleport example

examples/quantum/gate/teleport_scaffold.cpp

@@ -31,59 +31,33 @@
 #include "XACC.hpp"
 #include "EigenAccelerator.hpp"
 
-/**
- * FIXME For now, create a fake accelerator
- * This will later come from QVM
- */
-class IBM5Qubit: public xacc::quantum::QPUGate<5> {
-public:
-	virtual xacc::AcceleratorType getType() {
-		return xacc::AcceleratorType::qpu_gate;
-	}
-	virtual std::vector<xacc::IRTransformation> getIRTransformations() {
-		std::vector<xacc::IRTransformation> v;
-		return v;
-	}
-	virtual void execute(const std::shared_ptr<xacc::IR> ir) {
-	}
-	virtual ~IBM5Qubit() {
-	}
-
-protected:
-	bool canAllocate(const int N) {
-		return true;
-	}
-};
-
 // Quantum Kernel executing teleportation of
 // qubit state to another.
 const std::string src("__qpu__ teleport () {\n"
-						"   // Initialize qubit to 1\n"
+						"   cbit creg[2];\n"
 						"   X(qreg[0]);\n"
 						"   H(qreg[1]);\n"
 						"   CNOT(qreg[1],qreg[2]);\n"
 						"   CNOT(qreg[0],qreg[1]);\n"
 						"   H(qreg[0]);\n"
-						"   MeasZ(qreg[0]);\n"
-						"   MeasZ(qreg[1]);\n"
-						"   // cZ\n"
-						"   H(qreg[2]);\n"
-						"   CNOT(qreg[2], qreg[1]);\n"
-						"   H(qreg[2]);\n"
-						"   // cX = CNOT\n"
-						"   CNOT(qreg[2], qreg[0]);\n"
+						"   creg[0] = MeasZ(qreg[0]);\n"
+						"   creg[1] = MeasZ(qreg[1]);\n"
+						"   if(creg[0] == 1) Z(qreg[2]);\n"
+						"   if (creg[1] == 1) X(qreg[2]);\n"
 						"}\n");
 
 int main (int argc, char** argv) {
 
-	// Create a reference to the IBM5Qubit Accelerator
-	auto ibm_qpu = std::make_shared<xacc::quantum::EigenAccelerator<3>>();
+	using Simple3QubitAcc = xacc::quantum::EigenAccelerator<3>;
+
+	// Create a reference to the 3 qubit simulation Accelerator
+	auto qpu = std::make_shared<Simple3QubitAcc>();
 
 	// Allocate some qubits, give them a unique identifier...
-	auto qreg = ibm_qpu->allocate("qreg");
+	auto qreg = qpu->allocate("qreg");
 
 	// Construct a new Program
-	xacc::Program quantumProgram(ibm_qpu, src);
+	xacc::Program quantumProgram(qpu, src);
 
 	// Build the program
 	quantumProgram.build("--compiler scaffold --writeIR teleport.xir");
@@ -96,7 +70,6 @@ int main (int argc, char** argv) {
 
 	// Get the execution result
 	qreg->printState(std::cout);
-	auto bits = qreg->toBits();
 
 	return 0;
 }

quantum/gate/accelerators/Qubits.hpp

@@ -23,14 +23,27 @@ protected:
 public:
 	Qubits() :
 			state((int) std::pow(2, NumberOfQubits)) {
+		// Initialize to |000...000> state
+		state.setZero();
+		state(0) = 1.0;
 	}
 
 	void applyUnitary(Eigen::MatrixXcd& U) {
 		state = U * state;
 	}
 
+	QubitState& getState() {
+		return state;
+	}
+
+	void setState(QubitState& st) {
+		state = st;
+	}
+
 	void printState(std::ostream& stream) {
-		stream << state << "\n";
+		for (int i = 0; i < state.rows(); i++) {
+			stream << std::bitset<NumberOfQubits>(i).to_string() << " -> " << state(i) << "\n";
+		}
 	}
 };
 }

quantum/gate/accelerators/eigenaccelerator/EigenAccelerator.hpp

@@ -35,6 +35,7 @@
 #include "QasmToGraph.hpp"
 #include "GraphIR.hpp"
 #include <unsupported/Eigen/KroneckerProduct>
+#include <random>
 
 namespace xacc {
 namespace quantum {
@@ -54,15 +55,19 @@ public:
 	 * The constructor, create tensor gates
 	 */
 	EigenAccelerator() {
-		Eigen::MatrixXcd h(2,2), cnot(4,4), I(2,2), x(2,2);
+		Eigen::MatrixXcd h(2,2), cnot(4,4), I(2,2), x(2,2), p0(2,2), p1(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;
 		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));
 	}
 
 	/**
@@ -71,7 +76,7 @@ public:
 	 */
 	virtual void execute(const std::shared_ptr<xacc::IR> ir) {
 
-		auto qubitsType = std::dynamic_pointer_cast<Qubits<NQubits>>(this->bits);
+		auto qubits = std::dynamic_pointer_cast<Qubits<NQubits>>(this->bits);
 
 		// Cast to a GraphIR, if we can...
 		auto graphir = std::dynamic_pointer_cast<QuantumGraphIR>(ir);
@@ -81,6 +86,7 @@ public:
 
 		// Get the Graph
 		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);
@@ -91,59 +97,208 @@ public:
 			gateOperations.emplace_back(n);
 		}
 
+		std::cout << "Initial State:\n";
+		qubits->printState(std::cout);
+		std::cout << "\n";
 
-		Eigen::MatrixXcd U = Eigen::MatrixXcd::Identity(std::pow(2, NQubits),
-				std::pow(2, NQubits));
+		for (auto gate : gateOperations) {
 
-		while (layer < finalLayer) {
-
-			std::vector<CircuitNode> currentLayerGates;
-			std::copy_if(gateOperations.begin(), gateOperations.end(),
-					std::back_inserter(currentLayerGates),
-					[&](const CircuitNode& c) {return std::get<1>(c.properties) == layer;});
+			// Skip disabled gates...
+			if (!std::get<4>(gate.properties)) {
+				continue;
+			}
 
+			// Create a list of nQubits Identity gates
 			std::vector<Eigen::MatrixXcd> productList(NQubits);
 			for (int i = 0; i < NQubits; i++) {
 				productList[i] = gates["I"];
 			}
 
-			// Can parallize this...
-			for (auto n : currentLayerGates) {
+			// Create a local U gate
+			Eigen::MatrixXcd localU;
 
-				auto gateName = std::get<0>(n.properties);
-				auto actingQubits = std::get<3>(n.properties);
+			// 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);
 
-				if (gateName != "measure") {
-					auto gate = gates[gateName];
+			if (gateName == "conditional") {
 
-					if (actingQubits.size() == 1) {
-						productList[actingQubits[0]] = gate;
-					} else if (actingQubits.size() == 2) {
-						productList[actingQubits[0]] = gate;
-						productList.erase(productList.begin() + actingQubits[1]);
-					} else {
-						QCIError("Can only simulate one and two qubit gates.");
+				// If we hit a measured result then we
+				// need to figure out the measured qubit it
+				// depends on, then if its a 1, enable the disabled
+				// gates from this node to the next FinalState node
+				auto qubitWeNeed = std::get<3>(gate.properties)[0];
+				auto qubitFound = qubitIdToMeasuredResult.find(qubitWeNeed);
+				if (qubitFound == qubitIdToMeasuredResult.end()) {
+					QCIError("Invalid conditional node - this qubit has not been measured.");
+				}
+
+				auto result = qubitIdToMeasuredResult[qubitWeNeed];
+
+				auto currentNodeId = std::get<2>(gate.properties);
+				if (result == 1) {
+					// Walk over the next nodes until we hit a FinalState node
+					// set their enabled state to true
+					for (int i = currentNodeId+1; i < nNodes; i++) {
+						// Once we hit the next final state node, then break out
+						if (graph.getVertexProperty<0>(i) == "FinalState") {
+							break;
+						}
+						std::cout << "Enabling " << graph.getVertexProperty<0>(i) << "\n";
+						graph.setVertexProperty<4>(i, true);
 					}
-				} else {
-					// Setup measurement gate
 				}
 
+			} else if (gateName == "measure") {
+
+				// get rho
+				auto rho = qubits->getState() * qubits->getState().transpose();
+
+				productList[actingQubits[0]] = gates["p0"];
 				// Create a total unitary for this layer of the circuit
-				Eigen::MatrixXcd result = productList[0];
+				auto temp = productList[0];
 				for (int i = 1; i < productList.size(); i++) {
-					result = kroneckerProduct(result, productList[i]).eval();
+					temp = kroneckerProduct(temp, productList[i]).eval();
 				}
-				assert(result.rows() == std::pow(2, NQubits) && result.cols() == std::pow(2,NQubits));
 
-				// Update the circuit unitary matrix
-				U = result * U;
+				auto temp2 = temp * rho;
+				auto probZero = temp2.trace();
+				std::random_device rd;
+				std::mt19937 mt(rd());
+				std::uniform_real_distribution<double> dist(0, 1.0);
+				int result;
+				auto val = dist(mt);
+				std::cout << "Val: " << val << "\n";
+				if (val < std::real(probZero)) {
+					result = 0;
+					Eigen::VectorXcd newState = (temp * qubits->getState());
+					newState.normalize();
+					qubits->setState(newState);
+				} else {
+					result = 1;
+					productList[actingQubits[0]] = gates["p1"];
+					// Create a total unitary for this layer of the circuit
+					auto temp = productList[0];
+					for (int i = 1; i < productList.size(); i++) {
+						temp = kroneckerProduct(temp, productList[i]).eval();
+					}
+					Eigen::VectorXcd newState = (temp * qubits->getState());
+					newState.normalize();
+					qubits->setState(newState);
+				}
 
-			}
+				std::cout << "Measured qubit " << actingQubits[0] << " to be a " << result << ": prob was " << probZero << "\n";
+				qubitIdToMeasuredResult.insert(std::make_pair(actingQubits[0], result));
 
-			layer++;
-		}
+			} else {
+				if (gateName != "FinalState" && gateName != "InitialState") {
+					// Regular Gate operations...
 
-		qubitsType->applyUnitary(U);
+					if (actingQubits.size() == 1) {
+
+						// If this is a one qubit gate, just replace
+						// the currect I in the list with the gate
+						productList[actingQubits[0]] = gates[gateName];
+
+						// Create a total unitary for this layer of the circuit
+						localU = productList[0];
+						for (int i = 1; i < productList.size(); i++) {
+							localU =
+									kroneckerProduct(localU, productList[i]).eval();
+						}
+
+					} 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];
+						for (int i = 1; i < productList.size(); i++) {
+							localU =
+									kroneckerProduct(localU, productList[i]).eval();
+						}
+
+						// 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];
+						for (int i = 1; i < productList.size(); i++) {
+							temp =
+									kroneckerProduct(temp, productList[i]).eval();
+						}
+
+						// Sum them up
+						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));
+
+					qubits->applyUnitary(localU);
+
+					std::cout << "Current State after " << gateName << ":\n";
+					qubits->printState(std::cout);
+					std::cout << "\n" << localU << "\n";
+				}
+			}
+		}
+//
+//		while (layer < finalLayer) {
+//
+//			std::vector<CircuitNode> currentLayerGates;
+//			std::copy_if(gateOperations.begin(), gateOperations.end(),
+//					std::back_inserter(currentLayerGates),
+//					[&](const CircuitNode& c) {return std::get<1>(c.properties) == layer;});
+//
+//			std::vector<Eigen::MatrixXcd> productList(NQubits);
+//			for (int i = 0; i < NQubits; i++) {
+//				productList[i] = gates["I"];
+//			}
+//
+//			// Can parallize this...
+//			for (auto n : currentLayerGates) {
+//
+//				auto gateName = std::get<0>(n.properties);
+//				auto actingQubits = std::get<3>(n.properties);
+//
+//				if (gateName != "measure" || gateName != "FinalState" || gateName != "InitialState") {
+//					auto gate = gates[gateName];
+//
+//					if (actingQubits.size() == 1) {
+//						productList[actingQubits[0]] = gate;
+//					} else if (actingQubits.size() == 2) {
+//						productList[actingQubits[0]] = gate;
+//						productList.erase(productList.begin() + actingQubits[1]);
+//					} else {
+//						QCIError("Can only simulate one and two qubit gates.");
+//					}
+//				} else {
+//
+//					if (gateName == "conditional") {
+//
+//					} else if (gateName == "measure") {
+//
+//					}
+//				}
+//
+//				// Create a total unitary for this layer of the circuit
+//				Eigen::MatrixXcd result = productList[0];
+//				for (int i = 1; i < productList.size(); i++) {
+//					result = kroneckerProduct(result, productList[i]).eval();
+//				}
+//				assert(result.rows() == std::pow(2, NQubits) && result.cols() == std::pow(2,NQubits));
+//
+//				// Update the circuit unitary matrix
+//				U = result * U;
+//
+//			}
+//
+//			layer++;
+//		}
+//
+//		qubitsType->applyUnitary(U);
 	}
 
 	/**

quantum/gate/compilers/scaffold/ScaffoldCompiler.cpp

@@ -49,7 +49,11 @@ void ScaffoldCompiler::modifySource() {
 	kernelSource.erase(kernelSource.find("__qpu__"), 7);
 	kernelSource = std::string("module ") + kernelSource;
 
-	std::string qubitAllocationLine = "   qbit qreg[3];\n";
+	std::string qubitAllocationLine;// = "   qbit qreg[3];\n";
+
+	std::regex qbitName("qbit\\s.*");
+	qubitAllocationLine = (*std::sregex_iterator(kernelSource.begin(), kernelSource.end(),
+				qbitName)).str() + "\n";
 
 	// conditional on measurements
 	// FIXME FOR NOW WE ONLY ACCEPT format
@@ -59,7 +63,6 @@ void ScaffoldCompiler::modifySource() {
 	std::regex ifstmts("if\\s?\\(\\w+\\[\\w+\\]\\s?=.*\\s?\\)\\s?");
 	for (auto i = std::sregex_iterator(kernelSource.begin(), kernelSource.end(),
 			ifstmts); i != std::sregex_iterator(); ++i) {
-		std::cout << "HELLO WORLD: " << (*i).str() << "\n";
 		std::vector<std::string> splitVec;
 		std::string ifLine = (*i).str();
 		boost::trim(ifLine);

quantum/gate/utils/QasmToGraph.hpp

@@ -186,18 +186,18 @@ public:
 		int maxLayer = layer+1;
 
 		// Print info...
-		for (auto cn : gateOperations) {
-			std::cout << "Gate Operation: \n";
-			std::cout << "\tName: " << std::get<0>(cn.properties) << "\n";
-			std::cout << "\tLayer: " << std::get<1>(cn.properties) << "\n";
-			std::cout << "\tGate Vertex Id: " << std::get<2>(cn.properties) << "\n";
-			std::cout << "\tActing Qubits: ";
-			std::vector<int> qubits = std::get<3>(cn.properties);
-			for (auto v : qubits) {
-				std::cout << v << ", ";
-			}
-			std::cout << "\n\n";
-		}
+//		for (auto cn : gateOperations) {
+//			std::cout << "Gate Operation: \n";
+//			std::cout << "\tName: " << std::get<0>(cn.properties) << "\n";
+//			std::cout << "\tLayer: " << std::get<1>(cn.properties) << "\n";
+//			std::cout << "\tGate Vertex Id: " << std::get<2>(cn.properties) << "\n";
+//			std::cout << "\tActing Qubits: ";
+//			std::vector<int> qubits = std::get<3>(cn.properties);
+//			for (auto v : qubits) {
+//				std::cout << v << ", ";
+//			}
+//			std::cout << "\n\n";
+//		}
 
 		generateEdgesFromLayer(1, graph, gateOperations, 0);
 
@@ -250,7 +250,7 @@ public:
 		for (auto g : conditionalGraphs) {
 			CircuitNode node;
 			std::get<0>(node.properties) = "conditional";
-			std::get<2>(node.properties) = measurementIds[counter];
+			std::get<2>(node.properties) = id;//measurementIds[counter];
 			std::get<3>(node.properties) = measurementQubits[counter];
 			mainGraph.addVertex(node);