Unverified Commit 7e3f7cab authored by Mccaskey, Alex's avatar Mccaskey, Alex Committed by GitHub
Browse files

Merge pull request #460 from tnguyen-ornl/tnguyen/obs-grouping

Fixed https://github.com/eclipse/xacc/issues/399
parents 58e84c95 b2ca0b11
......@@ -270,6 +270,182 @@ PauliOperator::observe(std::shared_ptr<CompositeInstruction> function) {
return observed;
}
std::vector<std::shared_ptr<CompositeInstruction>>
PauliOperator::observe(std::shared_ptr<CompositeInstruction> function, const HeterogeneousMap &grouping_options) {
if (grouping_options.size() == 0) {
return observe(function);
}
// For this, we require the Accelerator info:
if (!grouping_options.pointerLikeExists<Accelerator>("accelerator")) {
xacc::error("'accelerator' is required for Observable grouping");
}
auto qpu = grouping_options.getPointerLike<Accelerator>("accelerator");
const bool shots_enabled = [&qpu](){
// Remote QPU's (IBM, etc.) -> shots
if (qpu->isRemote()) {
return true;
}
// This is not fully-compatible yet, but we need Accelerators to expose
// the 'shots' info to their getProperties() for this to work.
auto qpu_props = qpu->getProperties();
if (qpu_props.keyExists<int>("shots")) {
return qpu_props.get<int>("shots") > 0;
}
return false;
}();
if (!shots_enabled) {
// Log that we cannot do grouping (since shots is not set)
xacc::info(qpu->name() + " accelerator is not running 'shots' execution. "
"Observable grouping will be ignored.");
return observe(function);
}
// For this grouping, we only support *single* grouping,
// i.e. all sub-terms commute.
const int nbQubits = std::max<int>(function->nPhysicalBits(), nQubits());
const bool all_terms_commute = [this, nbQubits]() {
// Check that each qubit location has a **unique** basis
// across all terms:
// This is much faster than checking the commutes()
std::unordered_map<int, int> qubit_to_basis;
for (auto &inst : terms) {
Term spinInst = inst.second;
std::vector<int> meas_bits;
if (!spinInst.isIdentity()) {
auto [v, w] = spinInst.toBinaryVector(nbQubits);
assert(v.size() == w.size());
for (int i = 0; i < v.size(); ++i) {
if (v[i] != 0 || w[i] != 0) {
// 1, 2, 3 => X, Z, Y
int op_id = v[i] + 2 * w[i];
if (qubit_to_basis.find(i) == qubit_to_basis.end()) {
// Not seen this:
qubit_to_basis[i] = op_id;
} else {
if (qubit_to_basis[i] != op_id) {
// Different basis at this location.
return false;
}
}
}
}
}
}
// unique basis at each qubit location.
return true;
}();
if (!all_terms_commute) {
xacc::info("Not all terms commute with each other. Grouping cannot be done.");
return observe(function);
}
// Grouping can be done:
// all terms in this Observable will generate a **single** observed circuit,
// i.e. if there is a Pauli-X at qubit 3 (add Hadamard gate),
// there is no chance that another term with Pauli-Z or Pauli-Y at that location.
// Create a new GateQIR to hold the spin based terms
auto gateRegistry = xacc::getService<IRProvider>("quantum");
std::vector<std::shared_ptr<CompositeInstruction>> observed;
int counter = 0;
auto pi = xacc::constants::pi;
// If the incoming function has instructions that have
// their buffer_names set, then we need to set the
// new measurement instructions buffer names to be the same.
// Here we assume that all instructions have the same buffer name
std::string buf_name = "";
if (function->nInstructions() > 0 &&
!function->getInstruction(0)->getBufferNames().empty()) {
buf_name = function->getInstruction(0)->getBufferNames()[0];
}
// Specify the bit-ordering of the accelerator
// to decode later.
const std::string compositeName =
std::string("GroupObserve_") +
(qpu->getBitOrder() == Accelerator::BitOrder::MSB ? "MSB" : "LSB");
// Single observed circuit.
auto gateFunction =
gateRegistry->createComposite(compositeName, function->getVariables());
if (function->hasChildren()) {
gateFunction->addInstruction(function->clone());
}
for (auto arg : function->getArguments()) {
gateFunction->addArgument(arg, 0);
}
// Track the basis change per qubit
std::unordered_map<size_t, std::string> qubits_to_basis_change_inst;
// Populate GateQIR now...
for (auto &inst : terms) {
Term spinInst = inst.second;
auto termsMap = spinInst.ops();
std::vector<std::pair<int, std::string>> terms;
for (auto &kv : termsMap) {
if (kv.second != "I" && !kv.second.empty()) {
terms.push_back({kv.first, kv.second});
}
}
for (int i = terms.size() - 1; i >= 0; i--) {
std::size_t qbit = terms[i].first;
auto gateName = terms[i].second;
if (qubits_to_basis_change_inst.find(qbit) !=
qubits_to_basis_change_inst.end()) {
// Have seen this qubit before
const std::string previous_gate = qubits_to_basis_change_inst[qbit];
if (gateName != previous_gate) {
// Something wrong with the commute check:
std::stringstream error_ss;
error_ss << "Internal error: Qubit " << qbit << " was observed with "
<< previous_gate << "; but requested to change to "
<< gateName << " during grouping.";
xacc::error(error_ss.str());
}
// Nothing to do: same basis.
continue;
}
// First time see this qubit:
qubits_to_basis_change_inst[qbit] = gateName;
// Adding the gate:
if (gateName == "X") {
auto hadamard =
gateRegistry->createInstruction("H", {qbit});
if (!buf_name.empty()) {
hadamard->setBufferNames({buf_name});
}
gateFunction->addInstruction(hadamard);
} else if (gateName == "Y") {
auto rx =
gateRegistry->createInstruction("Rx", {qbit});
if (!buf_name.empty()) {
rx->setBufferNames({buf_name});
}
InstructionParameter p(pi / 2.0);
rx->setParameter(0, p);
gateFunction->addInstruction(rx);
}
}
}
// Add measure (all qubits)
for (size_t i = 0; i < nbQubits; ++i) {
auto meas = gateRegistry->createInstruction("Measure", {i});
if (!buf_name.empty())
meas->setBufferNames({buf_name});
xacc::InstructionParameter classicalIdx((int)i);
meas->setParameter(0, classicalIdx);
gateFunction->addInstruction(meas);
}
// std::cout << "Group observed circuit:\n" << gateFunction->toString() << "\n";
return {gateFunction};
}
std::pair<std::vector<int>, std::vector<int>>
Term::toBinaryVector(const int nQubits) {
// return v,w
......@@ -811,14 +987,60 @@ void PauliOperator::normalize() {
return;
}
double PauliOperator::calcExpValFromGroupedExecution(
std::shared_ptr<AcceleratorBuffer> buffer) {
assert(buffer->nChildren() == 1);
auto resultBuffer = buffer->getChildren()[0];
std::complex<double> energy =
getIdentitySubTerm() ? getIdentitySubTerm()->coefficient() : 0.0;
const auto bit_order = resultBuffer->name().find("MSB") != std::string::npos
? AcceleratorBuffer::BitOrder::MSB
: AcceleratorBuffer::BitOrder::LSB;
auto temp_buffer = xacc::qalloc(resultBuffer->size());
for (auto &inst : terms) {
Term spinInst = inst.second;
std::vector<int> meas_bits;
if (!spinInst.isIdentity()) {
// std::cout << "Term: " << inst.first << "\n";
auto [v, w] = spinInst.toBinaryVector(resultBuffer->size());
assert(v.size() == w.size());
for (int i = 0; i < v.size(); ++i) {
// std::cout << "v = " << v[i] << "; w = " << w[i] << "\n";
if (v[i] != 0 || w[i] != 0) {
// Has an operator here:
meas_bits.emplace_back(i);
}
}
temp_buffer->setMeasurements(resultBuffer->getMarginalCounts(meas_bits, bit_order));
const auto coeff = spinInst.coeff();
const double term_exp_val = temp_buffer->getExpectationValueZ();
// temp_buffer->print();
// std::cout << "Exp = " << term_exp_val << "\n";
// std::cout << "Coeff = " << coeff << "\n";
energy += (term_exp_val * coeff);
}
}
return energy.real();
}
double PauliOperator::postProcess(std::shared_ptr<AcceleratorBuffer> buffer,
const std::string &postProcessTask,
const HeterogeneousMap &extra_data) {
if (buffer->nChildren() == 1 &&
buffer->getChildren()[0]->name().find("GroupObserve") !=
std::string::npos &&
!buffer->getChildren()[0]->getMeasurementCounts().empty()) {
// std::cout << "Grouping post processing!\n";
return calcExpValFromGroupedExecution(buffer);
}
if (buffer->nChildren() < getNonIdentitySubTerms().size()) {
xacc::error(
"The buffer doesn't contain enough sub-buffers as expected. Expect: " +
std::to_string(getNonIdentitySubTerms().size()) +
"; Received: " + std::to_string(buffer->nChildren()));
xacc::error("The buffer doesn't contain enough sub-buffers as expected. "
"Expect: " +
std::to_string(getNonIdentitySubTerms().size()) +
"; Received: " + std::to_string(buffer->nChildren()));
return 0.0;
}
......
......@@ -247,6 +247,9 @@ public:
std::vector<std::shared_ptr<CompositeInstruction>>
observe(std::shared_ptr<CompositeInstruction> function) override;
std::vector<std::shared_ptr<CompositeInstruction>>
observe(std::shared_ptr<CompositeInstruction> function,
const HeterogeneousMap &grouping_options) override;
std::vector<std::shared_ptr<Observable>> getSubTerms() override {
std::vector<std::shared_ptr<Observable>> ret;
......@@ -349,6 +352,9 @@ public:
virtual double postProcess(std::shared_ptr<AcceleratorBuffer> buffer,
const std::string &postProcessTask,
const HeterogeneousMap &extra_data) override;
private:
double
calcExpValFromGroupedExecution(std::shared_ptr<AcceleratorBuffer> buffer);
};
} // namespace quantum
......
......@@ -490,6 +490,108 @@ TEST(PauliOperatorTester, checkNormalize) {
}
TEST(PauliOperatorTester, checkGroupingCommuteCheck) {
PauliOperator op0;
op0.fromString("(0, -1) Z0 Z1 + (0, -1) Z1 Z2 + (0, -1) Z2 Z0");
std::cout << op0.toString() << "\n";
PauliOperator op1;
op1.fromString("(0, -1) Z0 + (0, -1) Z1 + (0, -1) X2");
std::cout << op1.toString() << "\n";
PauliOperator op2;
op2.fromString("(0, -1) Z0 Z1 + (0, -1) Z1 Z2 + (0, -1) X3 Y4");
std::cout << op2.toString() << "\n";
PauliOperator op3;
op3.fromString("(0, -1) Z0 X1 + (0, -1) X1 Y2 + (0, -1) Y2 Z3");
std::cout << op3.toString() << "\n";
// No optimization allowed here:
PauliOperator op4_no_opt_1;
op4_no_opt_1.fromString(
"(0, -1) Z0 X1 + (0, -1) X1 Y2 + (0, -1) Y2 Z3 + (0, -1) X0");
std::cout << op4_no_opt_1.toString() << "\n";
PauliOperator op4_no_opt_2;
op4_no_opt_2.fromString("(0, -1) Z0 Z1 + (0, -1) Z1 Z2 + (0, -1) Z2 X0");
std::cout << op4_no_opt_2.toString() << "\n";
std::vector<PauliOperator> opt_cases{op0, op1, op2, op3};
for (auto &test_case : opt_cases) {
auto qpu = xacc::getAccelerator("qpp", {{"shots", 1024}});
auto gateRegistry = xacc::getService<xacc::IRProvider>("quantum");
auto f = gateRegistry->createComposite("f");
auto h0 = gateRegistry->createInstruction("H", 0);
auto h1 = gateRegistry->createInstruction("H", 1);
auto h2 = gateRegistry->createInstruction("H", 2);
f->addInstructions({h0, h1, h2});
auto observed = test_case.observe(f, {{"accelerator", qpu}});
EXPECT_EQ(observed.size(), 1);
}
// Check non-commute as well:
std::vector<PauliOperator> no_opt_cases{op4_no_opt_1, op4_no_opt_2};
for (auto &test_case : no_opt_cases) {
auto qpu = xacc::getAccelerator("qpp", {{"shots", 1024}});
auto gateRegistry = xacc::getService<xacc::IRProvider>("quantum");
auto f = gateRegistry->createComposite("f");
auto h0 = gateRegistry->createInstruction("H", 0);
auto h1 = gateRegistry->createInstruction("H", 1);
auto h2 = gateRegistry->createInstruction("H", 2);
f->addInstructions({h0, h1, h2});
auto observed = test_case.observe(f, {{"accelerator", qpu}});
EXPECT_EQ(observed.size(), test_case.getNonIdentitySubTerms().size());
}
}
TEST(PauliOperatorTester, checkGroupingQaoaPostProcessLSB) {
PauliOperator op;
op.fromString("(1.5,0) + (-0.5,0) Z0 Z1 + (-0.5,0) Z0 Z2 + (-0.5,0) Z1 Z2");
std::cout << op.toString() << "\n";
xacc::Observable *obs = &op;
auto qpu = xacc::getAccelerator("qpp", {{"shots", 1024}});
auto qaoa_ansatz =
xacc::createComposite("qaoa", {{"nbQubits", 3},
{"nbSteps", 1},
{"cost-ham", obs},
{"parameter-scheme", "Standard"}});
const std::vector<double> opt_params{0.308, 0.614205};
auto f = (*qaoa_ansatz)(opt_params);
auto observed = op.observe(f, {{"accelerator", qpu}});
EXPECT_EQ(observed.size(), 1);
auto buffer = xacc::qalloc(3);
qpu->execute(buffer, observed);
buffer->print();
auto exp_val = op.postProcess(
buffer, xacc::Observable::PostProcessingTask::EXP_VAL_CALC, {});
EXPECT_NEAR(exp_val, 2.0, 0.1);
}
TEST(PauliOperatorTester, checkGroupingQaoaPostProcessMSB) {
PauliOperator op;
op.fromString("(1.5,0) + (-0.5,0) Z0 Z1 + (-0.5,0) Z0 Z2 + (-0.5,0) Z1 Z2");
std::cout << op.toString() << "\n";
xacc::Observable *obs = &op;
// Aer => MSB
auto qpu = xacc::getAccelerator("aer", {{"shots", 1024}});
auto qaoa_ansatz =
xacc::createComposite("qaoa", {{"nbQubits", 3},
{"nbSteps", 1},
{"cost-ham", obs},
{"parameter-scheme", "Standard"}});
const std::vector<double> opt_params{0.308, 0.614205};
auto f = (*qaoa_ansatz)(opt_params);
auto observed = op.observe(f, {{"accelerator", qpu}});
EXPECT_EQ(observed.size(), 1);
auto buffer = xacc::qalloc(3);
qpu->execute(buffer, observed);
buffer->print();
auto exp_val = op.postProcess(
buffer, xacc::Observable::PostProcessingTask::EXP_VAL_CALC, {});
EXPECT_NEAR(exp_val, 2.0, 0.1);
}
int main(int argc, char **argv) {
xacc::Initialize(argc, argv);
::testing::InitGoogleTest(&argc, argv);
......
......@@ -152,10 +152,67 @@ void QAOA::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const {
kernel->expand(m);
}
// Handle Max-cut optimization on shots-based backends (including physical
// backends). We only want to execute a single circuit for observable with all
// commuting terms such as the maxcut Hamiltonian.
// Limitation: this grouping cannot handle gradient strategy at the moment.
// Observe the cost Hamiltonian:
auto kernels = m_costHamObs->observe(kernel);
auto kernels = [&] {
if (dynamic_cast<xacc::quantum::PauliOperator *>(m_costHamObs)) {
return m_costHamObs->observe(kernel, {{"accelerator", m_qpu}});
} else {
return m_costHamObs->observe(kernel);
}
}();
// Grouping is possible (no gradient strategy)
// TODO: Gradient strategy to handle grouping as well.
int iterCount = 0;
if (m_costHamObs->getNonIdentitySubTerms().size() > 1 &&
kernels.size() == 1 && !gradientStrategy) {
OptFunction f(
[&, this](const std::vector<double> &x, std::vector<double> &dx) {
auto tmpBuffer = xacc::qalloc(buffer->size());
std::vector<std::shared_ptr<CompositeInstruction>> fsToExec{
kernels[0]->operator()(x)};
m_qpu->execute(tmpBuffer, fsToExec);
double energy = m_costHamObs->postProcess(tmpBuffer);
// We will only have one child buffer for each parameter set.
assert(tmpBuffer->getChildren().size() == 1);
auto result_buf = tmpBuffer->getChildren()[0];
result_buf->addExtraInfo("parameters", x);
result_buf->addExtraInfo("energy", energy);
buffer->appendChild("Iter" + std::to_string(iterCount), result_buf);
std::stringstream ss;
ss << "Iter " << iterCount << ": E("
<< (!x.empty() ? std::to_string(x[0]) : "");
for (int i = 1; i < x.size(); i++) {
ss << "," << std::setprecision(3) << x[i];
if (i > 4) {
// Don't print too many params
ss << ", ...";
break;
}
}
ss << ") = " << std::setprecision(12) << energy;
xacc::info(ss.str());
iterCount++;
if (m_maximize)
energy *= -1.0;
return energy;
}, kernel->nVariables());
auto result = m_optimizer->optimize(f);
// Reports the final cost:
double finalCost = result.first;
if (m_maximize)
finalCost *= -1.0;
buffer->addExtraInfo("opt-val", ExtraInfo(finalCost));
buffer->addExtraInfo("opt-params", ExtraInfo(result.second));
return;
}
// Construct the optimizer/minimizer:
OptFunction f(
[&, this](const std::vector<double> &x, std::vector<double> &dx) {
......@@ -315,8 +372,28 @@ QAOA::execute(const std::shared_ptr<AcceleratorBuffer> buffer,
m_single_exec_kernel = kernel;
}
// Observe the cost Hamiltonian:
auto kernels = m_costHamObs->observe(kernel);
// Observe the cost Hamiltonian, with the input Accelerator:
// i.e. perform grouping (e.g. max-cut QAOA, Pauli) if possible:
auto kernels = [&] {
if (dynamic_cast<xacc::quantum::PauliOperator *>(m_costHamObs)) {
return m_costHamObs->observe(kernel, {{"accelerator", m_qpu}});
} else {
return m_costHamObs->observe(kernel);
}
}();
if (m_costHamObs->getNonIdentitySubTerms().size() > 1 &&
kernels.size() == 1) {
// Grouping was done:
// just execute the single observed kernel:
std::vector<std::shared_ptr<CompositeInstruction>> fsToExec{
kernels[0]->operator()(x)};
m_qpu->execute(buffer, fsToExec);
const double finalCost = m_costHamObs->postProcess(buffer);
// std::cout << "Compute energy from grouping: " << finalCost << "\n";
return { finalCost };
}
std::vector<double> coefficients;
std::vector<std::string> kernelNames;
std::vector<std::shared_ptr<CompositeInstruction>> fsToExec;
......
......@@ -124,7 +124,7 @@ TEST(QAOATester, checkInitialStateConstruction) {
auto acc = xacc::getAccelerator("qpp");
auto buffer = xacc::qalloc(2);
auto optimizer = xacc::getOptimizer("nlopt", {{"initial-parameters", random_vector(-2., 2., 8)}});
auto optimizer = xacc::getOptimizer("nlopt", {{"initial-parameters", random_vector(-2., 2., 12)}});
auto qaoa = xacc::getService<Algorithm>("QAOA");
// Create deuteron Hamiltonian
auto H_N_2 = xacc::quantum::getObservable(
......@@ -142,7 +142,7 @@ TEST(QAOATester, checkInitialStateConstruction) {
std::make_pair("optimizer", optimizer),
std::make_pair("observable", H_N_2),
// number of time steps (p) param
std::make_pair("steps", 4),
std::make_pair("steps", 6),
std::make_pair("initial-state", initial_program),
std::make_pair("parameter-scheme", "Standard")}));
qaoa->execute(buffer);
......@@ -243,6 +243,75 @@ TEST(QAOATester, checkMaxCut) {
// EXPECT_NEAR((*buffer)["opt-val"].as<double>(), 2.0, 1e-3);
}
TEST(QAOATester, checkMaxCutGrouping) {
auto acc = xacc::getAccelerator("qpp", {{"shots", 8192}});
auto buffer = xacc::qalloc(3);
xacc::set_verbose(true);
auto optimizer = xacc::getOptimizer(
"nlopt", {{"ftol", 0.001},
{"maximize", true},
{"initial-parameters", random_vector(-2., 2., 2)}});
auto qaoa = xacc::getService<Algorithm>("maxcut-qaoa");
auto graph = xacc::getService<xacc::Graph>("boost-digraph");
// Triangle graph
for (int i = 0; i < 3; i++) {
graph->addVertex();
}
graph->addEdge(0, 1);
graph->addEdge(0, 2);
graph->addEdge(1, 2);
const bool initOk = qaoa->initialize(
{std::make_pair("accelerator", acc),
std::make_pair("optimizer", optimizer),
std::make_pair("graph", graph),
// number of time steps (p) param
std::make_pair("steps", 1),
// "Standard" or "Extended"
std::make_pair("parameter-scheme", "Standard")});
qaoa->execute(buffer);
buffer->print();
std::cout << "Opt-val: " << (*buffer)["opt-val"].as<double>() << "\n";
// There seems to be a local minima at 1.5 as well...
EXPECT_NEAR((*buffer)["opt-val"].as<double>(), 2.0, 0.1);
}
TEST(QAOATester, checkP1TriangleGraphGroupingExpVal) {
auto acc = xacc::getAccelerator("aer", {{"shots", 8192}});
auto optimizer = xacc::getOptimizer(
"nlopt",
{{"maximize", true}, {"initial-parameters", random_vector(-2., 2., 2)}});
auto H = xacc::quantum::getObservable(
"pauli", std::string("1.5 I - 0.5 Z0 Z1 - 0.5 Z0 Z2 - 0.5 Z1 Z2"));
auto qaoa = xacc::getAlgorithm("QAOA", {{"accelerator", acc},
{"optimizer", optimizer},
{"observable", H},
{"steps", 1},
{"parameter-scheme", "Standard"}});
auto all_betas =
xacc::linspace(-xacc::constants::pi / 4., xacc::constants::pi / 4., 20);
auto all_gammas =
xacc::linspace(-xacc::constants::pi, xacc::constants::pi, 20);
for (auto gamma : all_gammas) {
for (auto beta : all_betas) {
auto buffer = xacc::qalloc(3);
auto cost = qaoa->execute(buffer, std::vector<double>{gamma, beta})[0];
auto d = 1;
auto e = 1;
auto f = 1;
auto theory = 3 * (.5 +
.25 * std::sin(4 * beta) * std::sin(gamma) *
(std::cos(gamma) + std::cos(gamma)) -
.25 * std::sin(2 * beta) * std::sin(2 * beta) *
(std::pow(std::cos(gamma), d + e - 2 * f)) *
(1 - std::cos(2 * gamma)));
// std::cout << "Cost = " << cost << "; expected = " << theory << "\n";
EXPECT_NEAR(cost, theory, 0.1);
}
}
}
int main(int argc, char **argv) {
xacc::Initialize(argc, argv);
::testing::InitGoogleTest(&argc, argv);
......
......@@ -77,6 +77,9 @@ public:
virtual std::vector<std::pair<int, int>> getConnectivity() override;
virtual BitOrder getBitOrder() override { return BitOrder::MSB; }
virtual HeterogeneousMap getProperties() override {
return {{"shots", m_shots}};