ReadoutErrorAcceleratorBufferPostprocessor.cpp 11.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
/*******************************************************************************
 * Copyright (c) 2017 UT-Battelle, LLC.
 * All rights reserved. This program and the accompanying materials
 * are made available under the terms of the Eclipse Public License v1.0
 * and Eclipse Distribution License v1.0 which accompanies this
 * distribution. The Eclipse Public License is available at
 * http://www.eclipse.org/legal/epl-v10.html and the Eclipse Distribution License
 * is available at https://eclipse.org/org/documents/edl-v10.php
 *
 * Contributors:
 *   Alexander J. McCaskey - initial API and implementation
 *******************************************************************************/
#include "ReadoutErrorAcceleratorBufferPostprocessor.hpp"
#include <boost/algorithm/string.hpp>
#include "XACC.hpp"

17
18
#include <boost/range/adaptor/sliced.hpp>

19
20
namespace xacc {
namespace quantum {
21

22
23
24
25
26
27
28
29
std::vector<std::shared_ptr<AcceleratorBuffer>> ReadoutErrorAcceleratorBufferPostprocessor::process(
		std::vector<std::shared_ptr<AcceleratorBuffer>> buffers) {

	// Goal here is to get all ap01, ap10, ..., and fix all expectation values
	int nQubits = ir.maxBit() + 1;
	std::string zeroStr = "";
	for (int i = 0; i < nQubits; i++) zeroStr += "0";

30
31
32
	int nKernels = ir.getKernels().size();
	int nRepititions = buffers.size() / nKernels;

Mccaskey, Alex's avatar
Mccaskey, Alex committed
33
34
35
36
37
	int nonIdentityKernels = 0;
	for (auto& k : ir.getKernels()) {
		if (k->nInstructions() > 0) {
			nonIdentityKernels++;
		}
38
39
	}

Mccaskey, Alex's avatar
Mccaskey, Alex committed
40
41
42
43
	if (buffers.size() % nonIdentityKernels != 0) {
		xacc::error("ReadoutError Postprocessor: Invalid number of buffers and kernels - " + std::to_string(buffers.size()) + ", " + std::to_string(nonIdentityKernels) );

	}
44
45

	std::vector<std::vector<std::shared_ptr<AcceleratorBuffer>>> bufvec;
Mccaskey, Alex's avatar
Mccaskey, Alex committed
46
47
	for (int i = 0; i < buffers.size(); i+=nonIdentityKernels) {
		std::vector<std::shared_ptr<AcceleratorBuffer>> sub(buffers.begin() + i, buffers.begin() + i + nonIdentityKernels);
48
49
50
51
52
53
		bufvec.push_back(sub);
	}

	std::vector<std::shared_ptr<AcceleratorBuffer>> fixedBuffers;

	for ( auto subList : bufvec) {
Mccaskey, Alex's avatar
Mccaskey, Alex committed
54
55
56
57
58
		std::vector<std::shared_ptr<Function>> nonIdentityKernels;
		for (int i = 0; i < nKernels; i++) {
			if (ir.getKernels()[i]->nInstructions() > 0) {
				nonIdentityKernels.push_back(ir.getKernels()[i]);
			}
59
60
		}

Mccaskey, Alex's avatar
Mccaskey, Alex committed
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
		auto nIk = nonIdentityKernels.size();

		std::map<int, std::pair<double,double>> errorRates;
		bool first = true;
		int counter = 0, qbitCount=0;
		std::vector<double> probs;
		for (int i = allTerms.size(); i < nIk; i++) {
			auto localBitStr = zeroStr;
			auto kernel = nonIdentityKernels[i];
			if (first) {
				// we have a p01 buffer
				auto bit = kernel->getInstruction(0)->bits()[0];
				localBitStr[nQubits - bit - 1] = '1';
				first = false;
			} else {
				// we have a p10 buffer
				first = true;
			}

80
			xacc::info(kernel->name() + " - Computing measurement probability for bit string = " + localBitStr);
Mccaskey, Alex's avatar
Mccaskey, Alex committed
81
82
83
84
85
86
87
88
89
90
91
92

			probs.push_back(subList[i]->computeMeasurementProbability(localBitStr));
			counter++;

			if (counter == 2) {
				errorRates.insert(
						{ qbitCount, { std::isnan(probs[0]) ? 0.0 : probs[0],
								std::isnan(probs[1]) ? 0.0 : probs[1] } });
				counter = 0;
				qbitCount++;
				probs.clear();
			}
93
94
		}

95

Mccaskey, Alex's avatar
Mccaskey, Alex committed
96
97
98
99
		for (auto& kv : errorRates) {
			std::stringstream s, s2, s3;
			s << "Qubit " << kv.first << ": p01 = " << kv.second.first << ", p10 = " << kv.second.second;
			xacc::info(s.str());
100
		}
101

102

103
	// Return new AcceleratorBuffers subtype, StaticExpValAcceleratorBuffer that has static
104

Mccaskey, Alex's avatar
Mccaskey, Alex committed
105
106
107
108
109
		std::map<std::string, double> oldExpects;
		for (int i = 0; i < allTerms.size(); i++) {
			xacc::info("Raw Expectatations: " + allTerms[i] + " = " + std::to_string(subList[i]->getExpectationValueZ()));
			oldExpects.insert({allTerms[i], subList[i]->getExpectationValueZ()});
		}
110

Mccaskey, Alex's avatar
Mccaskey, Alex committed
111
		auto fixed = fix_assignments(oldExpects, sites, errorRates);
112

Mccaskey, Alex's avatar
Mccaskey, Alex committed
113
		// constant fixed expectation value from the calculation
114

Mccaskey, Alex's avatar
Mccaskey, Alex committed
115
116
117
118
		for (int i = 0; i < allTerms.size(); i++) {
			auto staticBuffer = std::make_shared<StaticExpectationValueBuffer>(subList[i]->name(), subList[i]->size(), fixed[allTerms[i]]);
			fixedBuffers.push_back(staticBuffer);
		}
119

120
	}
121
	return fixedBuffers;
122
}
123
124
125

std::map<std::string, double> ReadoutErrorAcceleratorBufferPostprocessor::fix_assignments(
		std::map<std::string, double> oldExpects,
126
		std::map<std::string, std::vector<int>> sites, // GIVES KEY X0X1 -> [0,1]
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
		std::map<int, std::pair<double, double>> errorRates) {

	std::map<std::string, double> newExpects;

	for (auto& kv : sites) {

		if (kv.first != "I") {
			if (kv.second.size() == 1) {
				double p01 = errorRates[kv.second[0]].first;
				double p10 = errorRates[kv.second[0]].second;
				newExpects.insert({kv.first, exptZ(oldExpects[kv.first], p01, p10)});
			} else if (kv.second.size() == 2) {

				std::stringstream s,t;
				s << kv.first[0] << kv.first[1];
				t << kv.first[2] << kv.first[3];
				auto k0 = s.str();
				auto k1 = t.str();

				double ap01 = errorRates[kv.second[0]].first;
				double ap10 = errorRates[kv.second[0]].second;
				double bp01 = errorRates[kv.second[1]].first;
				double bp10 = errorRates[kv.second[1]].second;

151
152
				std::stringstream s2;
				newExpects.insert({kv.first, exptZZ(oldExpects[kv.first], oldExpects[k0], oldExpects[k1], ap10, ap01, bp10, bp01)});
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296

			} else {
				xacc::error("Correction for paulis with support on > 2 sites not implemented.");
			}
		}

	}

	return newExpects;
}

double ReadoutErrorAcceleratorBufferPostprocessor::exptZZ(double E_ZZ, double E_ZI, double E_IZ, double a01, double a10,
		double b01, double b10, bool averaged) {

    double a00 = 1 - a10;
    double a11 = 1 - a01;
    double b00 = 1 - b10;
    double b11 = 1 - b01;
    double E_II = 1;

	double E_ZZ_fixed = 0;
    if(averaged) {
        E_ZZ_fixed = (E_ZZ + (a01 - a10)*(- E_ZI - E_IZ + a01 - a10))/std::pow((-1 + a01 + a10),2);
    } else {
        E_ZZ_fixed = (-(2 - a00 * b01 - a11 * b01 - a00 * b10 - a11 * b10 +
        a01*a01 * (b00 * b01 + b10 * (2 * b01 + b11)) +
        a10*a10 * (b00 * b01 + b10 * (2 * b01 + b11)) +
        a01 * ((a00 + 2 * a10) * b01*b01 - 2 * b10 + a00 * b10*b10 + 2 * a10 * b10*b10 +
           b00 * (-1 + a11 * b01 + a00 * b10 + 2 * a10 * b10) - b11 +
           a11 * b10 * b11 + b01 * (-2 + 2 * a11 * b10 + a00 * b11 + 2 * a10 * b11)) +
         a10 * (a11 * b01*b01 - 2 * b10 + a11 * b10*b10 +
           b00 * (-1 + a00 * b01 + a11 * b10) - b11 + a00 * b10 * b11 +
           b01 * (-2 + 2 * a00 * b10 + a11 * b11))) * E_ZZ -
    a10 * b00 * E_ZI + a00 * b01 * E_ZI +
    a11 * b01 * E_ZI + a00 * a10 * b00 * b01 * E_ZI +
    a10*a10 * b00 * b01 * E_ZI - 2 * a00 * a11 * b01*b01 * E_ZI -
    a10 * a11 * b01*b01 * E_ZI - a00 * b10 * E_ZI -
    a11 * b10 * E_ZI + a10 * a11 * b00 * b10 * E_ZI +
    2 * a00 * a11 * b10*b10 * E_ZI + a10 * a11 * b10*b10 * E_ZI +
    a10 * b11 * E_ZI - a10 * a11 * b01 * b11 * E_ZI -
    a00 * a10 * b10 * b11 * E_ZI - a10*a10 * b10 * b11 * E_ZI -
    a10 * b00 * E_IZ - a00 * b01 * E_IZ +
    a11 * b01 * E_IZ + a00 * a10 * b00 * b01 * E_IZ +
    a10*a10 * b00 * b01 * E_IZ - a10 * a11 * b01*b01 * E_IZ -
    a00 * b10 * E_IZ + a11 * b10 * E_IZ -
    a10 * a11 * b00 * b10 * E_IZ +
    2 * a00 * a10 * b01 * b10 * E_IZ -
    a10 * a11 * b10*b10 * E_IZ - a10 * b11 * E_IZ +
    2 * a10*a10 * b00 * b11 * E_IZ -
    a10 * a11 * b01 * b11 * E_IZ +
    a00 * a10 * b10 * b11 * E_IZ + a10*a10 * b10 * b11 * E_IZ -
    a10 * b00 * E_II +
    a00 * b01 * E_II -
    a11 * b01 * E_II +
    a00 * a10 * b00 * b01 * E_II +
    a10*a10 * b00 * b01 * E_II +
    2 * a10 * a11 * b00 * b01 * E_II +
    a10 * a11 * b01*b01 * E_II -
    2 * a00 * a10 * a11 * b00 * b01*b01 * E_II -
    2 * a10*a10 * a11 * b00 * b01*b01 * E_II -
    a00 * b10 * E_II +
    a11 * b10 * E_II -
    a10 * a11 * b00 * b10 * E_II +
    2 * a00 * a10 * a11 * b00 * b01 * b10 * E_II +
    2 * a10*a10 * a11 * b00 * b01 * b10 * E_II -
    4 * a00 * a10 * a11 * b01*b01 * b10 * E_II -
    2 * a10*a10 * a11 * b01*b01 * b10 * E_II -
    a10 * a11 * b10*b10 * E_II +
    4 * a00 * a10 * a11 * b01 * b10*b10 * E_II +
    2 * a10*a10 * a11 * b01 * b10*b10 * E_II +
    a10 * b11 * E_II +
    a10 * a11 * b01 * b11 * E_II -
    2 * a10*a10 * a11 * b00 * b01 * b11 * E_II -
    a00 * a10 * b10 * b11 * E_II -
    a10*a10 * b10 * b11 * E_II -
    2 * a10 * a11 * b10 * b11 * E_II +
    2 * a10*a10 * a11 * b00 * b10 * b11 * E_II -
    2 * a00 * a10 * a11 * b01 * b10 * b11 * E_II -
    2 * a10*a10 * a11 * b01 * b10 * b11 * E_II +
    2 * a00 * a10 * a11 * b10*b10 * b11 * E_II +
    2 * a10*a10 * a11 * b10*b10 * b11 * E_II +
    a01*a01 * (-b10 * (2 * a00 * b01 * (-b01 + b10) * E_II +
           b11 * (E_ZI + E_IZ + (-1 - 2 * a00 * b01 -
                2 * a10 * b01 + 2 * a00 * b10 +
                2 * a10 * b10) * E_II)) +
       b00 * (2 * (a00 + a10) * b01*b01 * E_II -
          2 * b11 * (E_IZ + (a00 +
                2 * a10) * b10 * E_II) +
          b01 * (E_ZI - E_IZ + (-1 - 2 * a00 * b10 -
                2 * a10 * b10 + 2 * a00 * b11 +
                4 * a10 * b11) * E_II))) +
    a01 * (b11 * E_ZI - 2 * a10 * b01 * b11 * E_ZI -
       a11 * b10 * b11 * E_ZI - 2 * a11 * b01 * b10 * E_IZ +
       b11 * E_IZ - a11 * b10 * b11 * E_IZ +
       2 * a10 * a11 * b01*b01 * b10 * E_II -
       2 * a10 * a11 * b01 * b10*b10 * E_II -
       b11 * E_II +
       a11 * b10 * b11 * E_II -
       2 * a10*a10 * b01 * b10 * b11 * E_II +
       2 * a10 * a11 * b01 * b10 * b11 * E_II +
       2 * a10*a10 * b10*b10 * b11 * E_II -
       2 * a10 * a11 * b10*b10 * b11 * E_II +
       b00 * ((-1 + a11 * b01 + a00 * b10 +
             2 * a10 * b10) * E_ZI + (1 - a11 * b01 +
             a00 * b10) * E_IZ + (1 - 2 * a10*a10 * b01*b01 +
             2 * a10*a10 * b01 * b10 - 4 * a10*a10 * b01 * b11 + 4 * a10*a10 * b10 * b11 +
             a11 * (2 * a10 * b01*b01 - 2 * a10 * b10 * b11 +
                b01 * (-1 - 2 * a10 * b10 + 2 * a10 * b11)) +
             a00 * (-2 * (a10 - a11) * b01*b01 + b10 + 2 * a10 * b10 * b11 -
                2 * b01 * (1 + a11 * b10 +
                   a10 * (-b10 + b11)))) * E_II) +
        a00 * (-b01*b01 * (E_ZI - E_IZ + E_II +
             2 * a10 * b10 * E_II -
             4 * a11 * b10 * E_II) -
          b01 * (-2 * (a10 - 2 * a11) * b10*b10 * E_II +
             b11 * (E_ZI - E_IZ + E_II + 2 * a10 * b10 * E_II -
                2 * a11 * b10 * E_II)) +
          b10 * (2 * b11 * E_II +
             b10 * (E_ZI + E_IZ + E_II + 2 * a10 * b11 * E_II -
                2 * a11 * b11 * E_II)))))/(2 * ((-1 +
         a11 * (b01 + b10)) * (1 - a00 * (b01 + b10) +
         a10*a10 * (b00 + b10) * (b01 + b11) +
         a10 * (-b01 + b00 * (-1 + a00 * b01) - b10 + 2 * a00 * b01 * b10 - b11 +
            a00 * b10 * b11)) +
      a01*a01 * (b10 * (a00 * b01*b01 + (-1 + a00 * b10 + a10 * b10) * b11 +
            b01 * (-1 + a10 * b11 + a00 * (b10 + b11))) +
         b00 * ((a00 + a10) * b01*b01 + (-1 + a00 * b10 + 2 * a10 * b10) * b11 +
            b01 * (-1 + a10 * b10 + 2 * a10 * b11 + a00 * (b10 + b11)))) +
      a01 * (b01*b01 * (a10 * (-1 + a11 * b10) +
            a00 * (-1 + a10 * b10 + 2 * a11 * b10)) + (-1 + a00 * b10 +
            a10 * b10) * (-b11 + b10 * (-1 + a10 * b11 + a11 * b11)) +
         b01 * (1 - 2 * a10 * b11 + a10*a10 * b10 * b11 +
            a11 * b10 * (-2 + a10 * (b10 + b11)) +
            a00 * ((a10 + 2 * a11) * b10*b10 - b11 +
               b10 * (-2 + a10 * b11 + a11 * b11))) +
         b00 * (1 - a11 * b01 + a10*a10 * (b01 + b10) * (b01 + 2 * b11) +
            a00 * (b01 + b10) * (-1 + a11 * b01 + a10 * (b01 + b11)) +
            a10 * (a11 * b01*b01 - 2 * b10 - 2 * b11 + a11 * b10 * b11 +
               b01 * (-2 + a11 * (b10 + b11)))))));
    }
    return E_ZZ_fixed;

}

297
298
299
300
301
}
}