FireTensorAccelerator.hpp 10.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
/***********************************************************************************
 * Copyright (c) 2016, UT-Battelle
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *   * Redistributions of source code must retain the above copyright
 *     notice, this list of conditions and the following disclaimer.
 *   * Redistributions in binary form must reproduce the above copyright
 *     notice, this list of conditions and the following disclaimer in the
 *     documentation and/or other materials provided with the distribution.
 *   * Neither the name of the xacc nor the
 *     names of its contributors may be used to endorse or promote products
 *     derived from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 * Contributors:
 *   Initial API and implementation - Alex McCaskey
 *
 **********************************************************************************/
31
32
#ifndef QUANTUM_GATE_ACCELERATORS_EIGENACCELERATOR_HPP_
#define QUANTUM_GATE_ACCELERATORS_EIGENACCELERATOR_HPP_
33

34
#include "QPUGate.hpp"
35
#include "GraphIR.hpp"
36
#include "QuantumCircuit.hpp"
37
#include "SimulatedQubits.hpp"
38
#include <random>
39
40
41
42

namespace xacc {
namespace quantum {

43
double sqrt2 = std::sqrt(2.0);
44
using QuantumGraphIR = xacc::GraphIR<QuantumCircuit>;
45

46
/**
47
48
49
50
51
 * 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.
52
 */
53
template<const int NQubits>
54
class FireTensorAccelerator : virtual public QPUGate<SimulatedQubits<NQubits>> {
55
56
public:

57
58
59
	/**
	 * The constructor, create tensor gates
	 */
60
	FireTensorAccelerator() {
61
		fire::Tensor<2, fire::EigenProvider, std::complex<double>> h(2,2), cnot(4,4), I(2,2), x(2,2), p0(2,2), p1(2,2), z(2,2);
62
63
64
65
66
67
68
69
70
71
72
73
74
75
		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));
76
77
	}

78
	/**
79
80
	 * Execute the simulation. Requires both a valid SimulatedQubits buffer and
	 * Graph IR instance modeling the quantum circuit.
81
82
83
	 *
	 * @param ir
	 */
84
	virtual void execute(const std::string& bufferId, const std::shared_ptr<xacc::IR> ir) {
85

86
		// Get the requested qubit buffer
87
88
		auto qubits = this->allocatedBuffers[bufferId];
		if (!qubits) {
89
			XACCError("Invalid buffer id. Could not get qubit buffer.");
90
91
		}

92
		// Set the size
93
		int nQubits = qubits->size();
94
95

		// Cast to a GraphIR, if we can...
96
		auto graphir = std::dynamic_pointer_cast<QuantumGraphIR>(ir);
97
		if (!graphir) {
98
			XACCError("Invalid IR - this Accelerator only accepts GraphIR<Graph<CircuitNode>>.");
99
100
		}

101
		// Get the Graph and related info
102
		std::vector<CircuitNode> gateOperations;
103
		std::map<int, int> qubitIdToMeasuredResult;
104
105
106
107
		auto graph = graphir->getGraph();
		int nNodes = graph.order(), layer = 1;
		int finalLayer = graph.getVertexProperty<1>(nNodes - 1);

108
		// Get a vector of all gates
109
110
111
112
113
114
		for (int i = 0; i < nNodes; i++) {
			CircuitNode n;
			n.properties = graph.getVertexProperties(i);
			gateOperations.emplace_back(n);
		}

115
		// Loop over all gates in the circuit
116
		for (auto gate : gateOperations) {
117

118
119
120
121
			// Skip disabled gates...
			if (!std::get<4>(gate.properties)) {
				continue;
			}
122

123
			// Create a list of nQubits Identity gates
124
			std::vector<fire::Tensor<2, fire::EigenProvider, std::complex<double>>> productList;
125
			for (int i = 0; i < nQubits; i++) {
126
				productList.push_back(gates.at("I"));
127
128
			}

129
			// Create a local U gate, initialized to identity
130
			fire::Tensor<2, fire::EigenProvider, std::complex<double>> localU = gates.at("I");
131

132
133
			// Get the current gate anme
			auto gateName = std::get<0>(gate.properties);
134

135
136
			// Get the qubits we're acting on...
			auto actingQubits = std::get<3>(gate.properties);
137

138
			if (gateName == "conditional") {
139

140
141
142
143
144
145
146
				// 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()) {
147
					XACCError("Invalid conditional node - this qubit has not been measured.");
148
149
150
151
152
153
154
155
156
157
				}

				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
Mccaskey, Alex's avatar
Mccaskey, Alex committed
158
						if (std::get<0>(gateOperations[i].properties) == "FinalState") {
159
160
							break;
						}
Mccaskey, Alex's avatar
Mccaskey, Alex committed
161
						std::get<4>(gateOperations[i].properties) = true;
162
163
164
					}
				}

165
166
			} else if (gateName == "measure") {

167
168
				// get rho - outer product of state with itself
				auto rho = qubits->getState() * qubits->getState();
169

170
				// Create a total unitary for this layer of the circuit
171
172
				productList.at(actingQubits[0]) = gates.at("p0");
				auto Pi0 = productList.at(0);
173
				for (int i = 1; i < productList.size(); i++) {
174
					Pi0 = Pi0.kronProd(productList.at(i));
175
176
				}

Mccaskey, Alex's avatar
Mccaskey, Alex committed
177
				// Get probability qubit is a 0
178
179
180
181
				double probZero = 0.0;
				std::array<IndexPair, 1> contractionIndices;
				contractionIndices[0] = std::make_pair(1, 0);
				auto Prob0 = Pi0.contract(rho, contractionIndices);
182
				for (int i = 0; i < Prob0.dimension(0); i++) probZero += std::real(Prob0(i,i));
Mccaskey, Alex's avatar
Mccaskey, Alex committed
183
184

				// Make the measurement random...
185
186
187
188
189
190
191
				std::random_device rd;
				std::mt19937 mt(rd());
				std::uniform_real_distribution<double> dist(0, 1.0);
				int result;
				auto val = dist(mt);
				if (val < std::real(probZero)) {
					result = 0;
192
193
					qubits->applyUnitary(Pi0);
					qubits->normalize();
194
195
				} else {
					result = 1;
196
					productList.at(actingQubits[0]) = gates.at("p1");
197
					// Create a total unitary for this layer of the circuit
198
					auto Pi1 = productList.at(0);
199
					for (int i = 1; i < productList.size(); i++) {
200
						Pi1 = Pi1.kronProd(productList.at(i));
201
					}
202
203
					qubits->applyUnitary(Pi1);
					qubits->normalize();
204
				}
205

206
				// Record the measurement result
207
				qubitIdToMeasuredResult.insert(std::make_pair(actingQubits[0], result));
208

209
210
211
			} else {
				if (gateName != "FinalState" && gateName != "InitialState") {
					// Regular Gate operations...
212

213
214
215
216
217
					if (isParameterized(gate)) {
						auto g = getParameterizedGate(gate);
						gates.insert(std::make_pair(gateName, g));
					}

218
219
220
221
					if (actingQubits.size() == 1) {

						// If this is a one qubit gate, just replace
						// the currect I in the list with the gate
222
						productList.at(actingQubits[0]) = gates.at(gateName);
223
224

						// Create a total unitary for this layer of the circuit
225
						localU = productList.at(0);
226
						for (int i = 1; i < productList.size(); i++) {
227
							localU = localU.kronProd(productList.at(i));
228
229
230
231
232
						}

					} 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)
233
234
						productList.at(actingQubits[0]) = gates.at("p0");
						localU = productList.at(0);
235
						for (int i = 1; i < productList.size(); i++) {
236
							localU = localU.kronProd(productList.at(i));
237
238
239
						}

						// Now create the second term in the sum
240
241
242
						productList.at(actingQubits[0]) = gates.at("p1");
						productList.at(actingQubits[1]) = gates.at(gateName == "cnot" ? "x" : "I");
						auto temp = productList.at(0);
243
						for (int i = 1; i < productList.size(); i++) {
244
							temp = temp.kronProd(productList.at(i));
245
246
247
						}

						// Sum them up
248
						localU = localU + temp;
249
					} else {
250
						XACCError("Can only simulate one and two qubit gates.");
251
252
253
					}

					// Make sure that localU is the correct size
254
255
256
257
					assert(
							localU.dimension(0) == std::pow(2, nQubits)
									&& localU.dimension(1)
											== std::pow(2, nQubits));
258

259
					// Appy the unitary and update th state
260
261
262
263
264
					qubits->applyUnitary(localU);

				}
			}
		}
265
266
	}

267
268
269
	/**
	 * The destructor
	 */
270
	virtual ~FireTensorAccelerator() {}
271
272
273

protected:

274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
	bool isParameterized(CircuitNode& node) {
		return !std::get<5>(node.properties).empty();
	}

	fire::Tensor<2, fire::EigenProvider, std::complex<double>> getParameterizedGate(
			CircuitNode& node) {
		fire::Tensor<2, fire::EigenProvider, std::complex<double>> g(2, 2);
		if (std::get<0>(node.properties) == "rz") {

			// Fixme... How to avoid the double here???

			auto param = this->template getRuntimeParameter<double>(
					std::get<5>(node.properties)[0]);
			std::complex<double> i(0, 1);
			auto rotation = std::exp(i * param);
			g.setValues( { { 1, 0 }, { 0, rotation } });
		} else {
291
			XACCError("We don't know what this gate is... yet.");
292
293
294
295
296
		}

		return g;

	}
297
298
299
	/**
	 * Mapping of gate names to actual gate matrices.
	 */
300
	std::map<std::string, fire::Tensor<2, fire::EigenProvider, std::complex<double>>> gates;
301
};
302
}
303
304
305
306
}



307

308
#endif