Commit aebaafd3 authored by Mccaskey, Alex's avatar Mccaskey, Alex
Browse files

Finished 2 site assignment error preprocessor


Signed-off-by: Mccaskey, Alex's avatarAlex McCaskey <mccaskeyaj@ornl.gov>
parent 5156d6a9
......@@ -24,43 +24,237 @@ std::vector<std::shared_ptr<AcceleratorBuffer>> ReadoutErrorAcceleratorBufferPos
std::string zeroStr = "";
for (int i = 0; i < nQubits; i++) zeroStr += "0";
std::map<std::string, double> errorRates;
for (auto& kv : measurements) {
auto kernelIdx = kv.second;
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 < buffers.size(); i++) {
auto localBitStr = zeroStr;
std::stringstream s;
s << kv.first[1] << kv.first[2];
if (boost::contains(s.str(), "01")) {
auto kernel = ir.getKernels()[kernelIdx];
if (first) {
// we have a p01 buffer
auto kernel = ir.getKernels()[i];
auto bit = kernel->getInstruction(0)->bits()[0];
localBitStr[nQubits - bit - 1] = '1';
first = false;
} else {
// we have a p10 buffer
first = true;
}
std::cout << "HI: " << kv.first << ", " << kv.second << ", " << localBitStr << "\n";
// get bit string from
auto prob = buffers[kv.second]->computeMeasurementProbability(localBitStr);
probs.push_back(buffers[i]->computeMeasurementProbability(localBitStr));
counter++;
errorRates.insert({kv.first, std::isnan(prob) ? 0.0 : prob});
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();
}
}
for (auto& kv : errorRates) {
std::stringstream s, s2;
s2 << kv.first[1] << kv.first[2];
s << kv.first[0]
<< (s2.str() == "01" ?
" probability expected 0 but got 1 error rate = " :
" probability expected 1 but got 0 error rate = ");
xacc::info("Qubit " + s.str() + std::to_string(kv.second));
std::stringstream s, s2, s3;
s << "Qubit " << kv.first << " p01 = " << kv.second.first << ", p10 = " << kv.second.second;
xacc::info(s.str());
}
// Return new AcceleratorBuffers subtype, StaticExpValAcceleratorBuffer that has static
std::map<std::string, double> oldExpects;
for (int i = 0; i < allTerms.size(); i++) {
std::cout << "OLDEXPECTS: " << allTerms[i] << ", " << buffers[i]->getExpectationValueZ() << "\n";
oldExpects.insert({allTerms[i], buffers[i]->getExpectationValueZ()});
}
auto fixed = fix_assignments(oldExpects, sites, errorRates);
// constant fixed expectation value from the calculation
return std::vector<std::shared_ptr<AcceleratorBuffer>>{};
std::vector<std::shared_ptr<AcceleratorBuffer>> fixedBuffers;
for (int i = 0; i < allTerms.size(); i++) {
auto staticBuffer = std::make_shared<StaticExpectationValueBuffer>(buffers[i]->name(), buffers[i]->size(), fixed[allTerms[i]]);
fixedBuffers.push_back(staticBuffer);
}
return fixedBuffers;
}
std::map<std::string, double> ReadoutErrorAcceleratorBufferPostprocessor::fix_assignments(
std::map<std::string, double> oldExpects,
std::map<std::string, std::vector<int>> sites,
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;
newExpects.insert({kv.first, exptZZ(oldExpects[kv.first], oldExpects[k0], oldExpects[k1], ap01, ap10, bp01, bp10)});
} 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;
}
}
}
......
......@@ -39,15 +39,30 @@ class ReadoutErrorAcceleratorBufferPostprocessor : public AcceleratorBufferPostp
protected:
std::map<std::string, int> measurements;
std::map<std::string, int> extraKernels;
std::map<std::string, std::vector<int>> sites;
std::vector<std::string> allTerms;
IR& ir;
double exptZ(double E_Z, double p01, double p10) {
auto p_p = p01 + p10;
auto p_m = p01 - p10;
return (E_Z-p_m)/(1.0 - p_p);
}
double exptZZ(double E_ZZ, double E_ZI, double E_IZ, double ap01, double ap10,
double bp01, double bp10, bool averaged = false);
std::map<std::string, double> fix_assignments(
std::map<std::string, double> oldExpects,
std::map<std::string, std::vector<int>> sites,
std::map<int, std::pair<double, double>> errorRates);
public:
ReadoutErrorAcceleratorBufferPostprocessor(IR& i,
std::map<std::string, int> extraKernelIndexMap,
std::map<std::string, int> measurementIndexMap) : ir(i),
extraKernels(extraKernelIndexMap), measurements(measurementIndexMap) {
std::map<std::string, std::vector<int>> sitesMap,
std::vector<std::string> orderedTerms) :
ir(i), sites(sitesMap), allTerms(
orderedTerms) {
}
virtual std::vector<std::shared_ptr<AcceleratorBuffer>> process(std::vector<std::shared_ptr<AcceleratorBuffer>> buffers);
......
......@@ -28,28 +28,32 @@ std::shared_ptr<AcceleratorBufferPostprocessor> ReadoutErrorIRPreprocessor::proc
auto gateRegistry = GateInstructionRegistry::instance();
// Search IR Functions and construct Pauli Term strings, then add any Pauli that is not there
std::vector<std::string> pauliTerms;
std::vector<std::string> orderedPauliTerms;
std::vector<std::map<int, std::string>> pauliTerms;
for (auto kernel : ir.getKernels()) {
std::string pauliStr = "";
CountGatesOfTypeVisitor<Measure> v(kernel);
bool allZTerm = false;
if (kernel->nInstructions() == v.countGates()) {
allZTerm = true;
}
std::string pauliTerm = "";
std::map<int, std::string> pauliTerm;
for (auto inst : kernel->getInstructions()) {
auto bit = inst->bits()[0];
if (allZTerm) {
pauliTerm = "Z" + std::to_string(bit) + " " + pauliTerm;
pauliTerm[bit] = "Z";
pauliStr = "Z" + std::to_string(bit) + pauliStr;
continue;
}
if (inst->getName() == "H") {
pauliTerm = "X" + std::to_string(bit) + " " + pauliTerm;
pauliStr = "X" + std::to_string(bit) + pauliStr;
pauliTerm[bit] = "X";
} else if (inst->getName() == "Rx") {
pauliTerm = "Y" + std::to_string(bit) + " " + pauliTerm;
pauliStr = "X" + std::to_string(bit) + pauliStr;
pauliTerm[bit] = "Y";
} else if (inst->getName() == "Measure") {
// do nothing
continue;
......@@ -60,18 +64,41 @@ std::shared_ptr<AcceleratorBufferPostprocessor> ReadoutErrorIRPreprocessor::proc
}
}
boost::trim(pauliTerm);
orderedPauliTerms.push_back(pauliStr);
pauliTerms.push_back(pauliTerm);
}
std::cout << "TERMSMaps:\n";
std::map<std::string, std::vector<int>> sites;
for (auto a : pauliTerms) {
std::stringstream s;
std::vector<int> tmp;
for (auto& kv : a) {
s << kv.second << kv.first;
tmp.push_back(kv.first);
std::cout << kv.second << kv.first << " ";
}
sites.insert({s.str(), tmp});
std::cout << "\n";
}
std::cout << "SITES:\n";
for (auto & kv : sites) {
std::cout << kv.first << " -> (";
for (auto i : kv.second) {
std::cout << i << " ";
}
std::cout << ")\n";
}
std::vector<std::string> extraKernelsNeeded;
for (auto t : pauliTerms) {
std::vector<std::string> ops;
boost::split(ops, t, boost::is_any_of(" "));
if (ops.size() > 1) {
for (auto o : ops) {
if (std::find(pauliTerms.begin(), pauliTerms.end(), o) == pauliTerms.end()) {
extraKernelsNeeded.push_back(o);
if (t.size() > 1) {
for (auto& kv : t) {
auto termStr = kv.second + std::to_string(kv.first);
if (!sites.count(termStr)) {
extraKernelsNeeded.push_back(termStr);
}
}
}
......@@ -80,9 +107,8 @@ std::shared_ptr<AcceleratorBufferPostprocessor> ReadoutErrorIRPreprocessor::proc
std::sort( extraKernelsNeeded.begin(), extraKernelsNeeded.end() );
extraKernelsNeeded.erase( std::unique( extraKernelsNeeded.begin(), extraKernelsNeeded.end() ), extraKernelsNeeded.end() );
std::map<std::string, int> extraKernelsMap, measureMap;
for (auto o : extraKernelsNeeded) {
std::cout << "EXTRA: " << o << "\n";
int nKernels = ir.getKernels().size();
auto extraKernel = std::make_shared<GateFunction>(o, "readout-error-extra");
......@@ -109,7 +135,7 @@ std::shared_ptr<AcceleratorBufferPostprocessor> ReadoutErrorIRPreprocessor::proc
ir.addKernel(extraKernel);
extraKernelsMap.insert({o, nKernels});
orderedPauliTerms.push_back(o);
}
int nKernels = ir.getKernels().size();
......@@ -136,22 +162,11 @@ std::shared_ptr<AcceleratorBufferPostprocessor> ReadoutErrorIRPreprocessor::proc
ir.addKernel(f10);
ir.addKernel(f01);
measureMap.insert({std::to_string(qbit) + "01", nKernels});
measureMap.insert({std::to_string(qbit) + "10", nKernels+1});
qbit++;
}
for (auto& kv : extraKernelsMap) {
std::cout << kv.first << ", " << kv.second << "\n";
}
for (auto& kv : measureMap) {
std::cout << kv.first << ", " << kv.second << "\n";
}
// Construct a ReadoutErrorABPostprocessor
return std::make_shared<ReadoutErrorAcceleratorBufferPostprocessor>(ir, extraKernelsMap, measureMap);
return std::make_shared<ReadoutErrorAcceleratorBufferPostprocessor>(ir, sites, orderedPauliTerms);
}
}
......
......@@ -88,15 +88,21 @@ std::shared_ptr<IR> createXACCIR(std::unordered_map<std::string, Term> terms) {
return newIr;
}
class FakeAB : public AcceleratorBuffer {
public:
FakeAB(const std::string& str, const int N) : AcceleratorBuffer(str, N) {}
virtual const double getExpectationValueZ() {
return 1.0;
}
};
BOOST_AUTO_TEST_CASE(checkSimple) {
// (-2.143303525+0j)*X0*X1 + (-3.91311896+0j)*X1*X2 +
// (-2.143303525+0j)*Y0*Y1 + (-3.91311896+0j)*Y1*Y2 + (0.218290555+0j)*Z0 + (-6.125+0j)*Z1 + (-9.625+0j)*Z2
// needs x0, x1, x2, y0, y1, y2
xacc::Initialize();
xacc::setOption("n-qubits", "3");
std::unordered_map<std::string, Term> test {{"X0X1", {{0,"X"}, {1,"X"}}},
{"X1X2", {{1,"X"}, {2,"X"}}},
{"Y0Y1", {{0,"Y"}, {1,"Y"}}},
......@@ -118,8 +124,15 @@ BOOST_AUTO_TEST_CASE(checkSimple) {
std::vector<std::shared_ptr<AcceleratorBuffer>> buffers;
for (int i = 0; i < ir->getKernels().size(); i++) {
buffers.push_back(std::make_shared<AcceleratorBuffer>("b"+std::to_string(i), 3));
buffers.push_back(std::make_shared<FakeAB>("b"+std::to_string(i), 3));
}
auto fixedBuffers = bufferProcessor->process(buffers);
BOOST_VERIFY(test.size() + nExtraKernels == fixedBuffers.size());
for (int i = 0; i < test.size(); i++) {
BOOST_VERIFY(1.0 == fixedBuffers[i]->getExpectationValueZ());
}
bufferProcessor->process(buffers);
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment