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

Merge pull request #201 from danclaudino/adapt-vqe-pr

Implementation of ADAPT-VQE and Parameter Shift
parents a30cc096 159ce1f5
Pipeline #106684 failed with stage
in 8 minutes and 48 seconds
import xacc
qpu = xacc.getAccelerator('tnqvm')
optimizer = xacc.getOptimizer('nlopt',{'nlopt-optimizer':'l-bfgs'})
buffer = xacc.qalloc(4)
geom = '''
H 0.000000 0.0 0.0
H 0.0 0.0 .7474
'''
H = xacc.getObservable('pyscf', {'basis': 'sto-3g', 'geometry': geom})
print(H.toString())
adapt = xacc.getAlgorithm('adapt-vqe', {'accelerator': qpu,
'optimizer': optimizer,
'observable': H,
'nElectrons': 2,
'gradient-strategy': 'parameter-shift-gradient',
'print-operators': True,
'threshold': 1.0e-2,
'print-threshold': 1.0e-10,
'maxiter': 2,
'pool': "singlet-adapted-uccsd"})
# execute
adapt.execute(buffer)
......@@ -24,23 +24,26 @@ namespace quantum {
FermionTerm &FermionTerm::operator*=(const FermionTerm &v) noexcept {
coeff() *= std::get<0>(v);
// std::cout << "FermionTerm: " << id() << ", " << FermionTerm::id(std::get<1>(v)) << "\n";
//std::cout << "FermionTerm: " << id() << ", " << FermionTerm::id(std::get<1>(v)) << "\n";
auto otherOps = std::get<1>(v);
for (auto &kv : otherOps) {
auto site = kv.first;
auto c_or_a = kv.second;
// std::cout << "HELLO: " << site << ", " << std::boolalpha << c_or_a << "\n";
//std::cout << "\n\nHELLO: " << site << ", " << std::boolalpha << c_or_a << "\n";
Operators o = ops();
if (!o.empty()) {
auto it = std::find_if(o.begin(), o.end(),
[&](const std::pair<int, bool> &element) {
return element.first == site;
});
// std::cout << it->first << ", " << std::boolalpha << it->second << "\n";
//std::cout << it->first << ", " << std::boolalpha << it->second << "\n";
if (it->first == site) {
if (it->second && c_or_a) {
if (it->second == c_or_a) {
// zero out this FermionTerm
ops().clear();
} else {//this adds the adjoint of operators whose sites are already in the product
ops().push_back({site, c_or_a});
}
} else {
......@@ -104,6 +107,8 @@ const std::string FermionOperator::toString() {
s << c << " " << std::get<2>(kv.second) << " ";
Operators ops = std::get<1>(kv.second);
// I changed this to get the operators printed in the order they're built
/*
std::vector<int> creations, annhilations;
for (auto &t : ops) {
if (t.second) {
......@@ -112,18 +117,25 @@ const std::string FermionOperator::toString() {
annhilations.push_back(t.first);
}
}
for (auto &t : creations) {
s << t << "^" << std::string(" ");
}
for (auto &t : annhilations) {
s << t << std::string(" ");
}*/
for (auto &t : ops) {
if (t.second) {
s << t.first << "^" << std::string(" ");
}
else{
s << t.first << std::string(" ");
}
}
s << "+ ";
}
// std::cout << "tostring " << s.str() << "\n";
auto r = s.str().substr(0, s.str().size() - 2);
xacc::trim(r);
return r;
......@@ -205,20 +217,31 @@ FermionOperator::operator*=(const FermionOperator &v) noexcept {
std::unordered_map<std::string, FermionTerm> newTerms;
for (auto &kv : terms) {
for (auto &vkv : v.terms) {
auto multTerm = kv.second * vkv.second;
FermionTerm multTerm;
if (kv.first == "I"){
multTerm = vkv.second;
} else {
multTerm = kv.second * vkv.second;
}
if (!multTerm.ops().empty()) {
auto id = multTerm.id();
if (!newTerms.insert({id, multTerm}).second) {
if (!newTerms.insert({id, multTerm}).second && kv.first != "I") {
newTerms.at(id).coeff() += multTerm.coeff();
}
if (!newTerms.insert({id, multTerm}).second && kv.first == "I") {
newTerms.at(id).coeff() = std::get<0>(kv.second) * std::get<0>(vkv.second);
}
if (std::abs(newTerms.at(id).coeff()) < 1e-12) {
newTerms.erase(id);
}
}
}
}
terms = newTerms;
return *this;
}
......@@ -258,6 +281,13 @@ FermionOperator::operator*=(const std::complex<double> v) noexcept {
return *this;
}
std::shared_ptr<Observable> FermionOperator::commutator(std::shared_ptr<Observable> op) {
FermionOperator& A = *std::dynamic_pointer_cast<FermionOperator>(op);
std::shared_ptr<FermionOperator> commutatorHA = std::make_shared<FermionOperator>((*this) * A - A * (*this));
return commutatorHA;
}
} // namespace quantum
} // namespace xacc
......
......@@ -204,6 +204,9 @@ void fromOptions(const HeterogeneousMap& options) override {
return;
}
std::shared_ptr<Observable>
commutator(std::shared_ptr<Observable> obs) override;
};
} // namespace quantum
......
......@@ -69,6 +69,30 @@ TEST(FermionOperatorTester, checkFromStr2) {
FermionOperator op(src);
std::cout << op.toString() << "\n";
}
TEST(FermionOperatorTester, checkMult) {
//FermionOperator op1("(-0.5,0) 3^ 2^ 0 1 + (0.5,-0) 1^ 0^ 2 3 + (0.5,-0) 0^ 1^ 3 2 + (0.5,0) 2^ 3^ 0 1 + (0.5,0) 3^ 3^ 1 1 + (-0.5,0) 1^ 1^ 3 3 + (-0.5,0) 1^ 0^ 3 2 + (-0.5,0) 2^ 3^ 1 0"), op2("(0.708024, 0)");
FermionOperator op1("3^ 2^ 3 2"), op2("0^ 3 1^ 2");
for (auto b : op1.getTerms()){std::cout << "First " << b.second.id() << "\n";}
for (auto b : op2.getTerms()){std::cout << "Second " << b.second.id() << "\n";}
FermionTerm mult;
for (auto a : op1.getTerms()){
for (auto b : op2.getTerms()){
mult = a.second * b.second;
}
}
for (auto k : mult.ops()){
std::cout << k.first << " " << k.second << "\n";
}
//auto mult = op1.getTerms().second * op2.getTerms().second;
std::cout << "Print op =" << (op1*op2).toString() << "\n\n";
}
int main(int argc, char** argv) {
xacc::Initialize(argc,argv);
::testing::InitGoogleTest(&argc, argv);
......
......@@ -760,6 +760,14 @@ void PauliOperator::fromXACCIR(std::shared_ptr<IR> ir) {
}
}
std::shared_ptr<Observable> PauliOperator::commutator(std::shared_ptr<Observable> op) {
PauliOperator& A = *std::dynamic_pointer_cast<PauliOperator>(op);
std::shared_ptr<PauliOperator> commutatorHA = std::make_shared<PauliOperator>((*this) * A - A * (*this));
return commutatorHA;
}
} // namespace quantum
} // namespace xacc
......
......@@ -339,6 +339,9 @@ public:
const std::string name() const override { return "pauli"; }
const std::string description() const override { return ""; }
void fromOptions(const HeterogeneousMap &options) override { return; }
std::shared_ptr<Observable>
commutator(std::shared_ptr<Observable> obs) override;
};
} // namespace quantum
......
......@@ -10,6 +10,7 @@
# Contributors:
# Alexander J. McCaskey - initial API and implementation
# *******************************************************************************/
add_subdirectory(adapt_vqe)
add_subdirectory(vqe)
add_subdirectory(rdm)
add_subdirectory(ml)
......@@ -17,6 +18,7 @@ add_subdirectory(rotoselect)
add_subdirectory(qpt)
add_subdirectory(qaoa)
add_subdirectory(qpe)
add_subdirectory(gradient_strategies)
file(GLOB PYDECORATORS ${CMAKE_CURRENT_SOURCE_DIR}/vqe/python/*.py
${CMAKE_CURRENT_SOURCE_DIR}/ml/ddcl/python/*.py)
......
# *******************************************************************************
# Copyright (c) 2019 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 v.10 which accompany 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
# *******************************************************************************/
set(LIBRARY_NAME xacc-algorithm-adapt_vqe)
file(GLOB SRC *.cpp)
usfunctiongetresourcesource(TARGET ${LIBRARY_NAME} OUT SRC)
usfunctiongeneratebundleinit(TARGET ${LIBRARY_NAME} OUT SRC)
add_library(${LIBRARY_NAME} SHARED ${SRC})
target_include_directories(
${LIBRARY_NAME}
PUBLIC .)
target_link_libraries(${LIBRARY_NAME} PUBLIC xacc CppMicroServices PRIVATE xacc-quantum-gate)
set(_bundle_name xacc_algorithm_adapt_vqe)
set_target_properties(${LIBRARY_NAME}
PROPERTIES COMPILE_DEFINITIONS
US_BUNDLE_NAME=${_bundle_name}
US_BUNDLE_NAME
${_bundle_name})
usfunctionembedresources(TARGET
${LIBRARY_NAME}
WORKING_DIRECTORY
${CMAKE_CURRENT_SOURCE_DIR}
FILES
manifest.json)
if(APPLE)
set_target_properties(${LIBRARY_NAME}
PROPERTIES INSTALL_RPATH "@loader_path/../lib")
set_target_properties(${LIBRARY_NAME}
PROPERTIES LINK_FLAGS "-undefined dynamic_lookup")
else()
set_target_properties(${LIBRARY_NAME}
PROPERTIES INSTALL_RPATH "$ORIGIN/../lib")
set_target_properties(${LIBRARY_NAME} PROPERTIES LINK_FLAGS "-shared")
endif()
if(XACC_BUILD_TESTS)
add_subdirectory(tests)
endif()
install(TARGETS ${LIBRARY_NAME} DESTINATION ${CMAKE_INSTALL_PREFIX}/plugins)
\ No newline at end of file
/*******************************************************************************
* Copyright (c) 2020 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:
* Daniel Claudino - initial API and implementation
*******************************************************************************/
#include "adapt_vqe.hpp"
#include "FermionOperator.hpp"
#include "ObservableTransform.hpp"
#include "xacc.hpp"
#include "xacc_service.hpp"
#include "xacc_observable.hpp"
#include "Circuit.hpp"
#include "AlgorithmGradientStrategy.hpp"
#include <memory>
#include <iomanip>
#include <sstream>
#include <string>
using namespace xacc;
using namespace xacc::quantum;
namespace xacc {
namespace algorithm {
bool ADAPT_VQE::initialize(const HeterogeneousMap &parameters) {
if (!parameters.pointerLikeExists<Observable>("observable")) {
xacc::info("Obs was false\n");
return false;
}
if (!parameters.pointerLikeExists<Accelerator>("accelerator")) {
xacc::info("Acc was false\n");
return false;
}
if(!parameters.stringExists("pool")){
return false;
}
optimizer = parameters.get<std::shared_ptr<Optimizer>>("optimizer");
observable = parameters.get<std::shared_ptr<Observable>>("observable");
accelerator = parameters.get<std::shared_ptr<Accelerator>>("accelerator");
nElectrons = parameters.get<int>("nElectrons");
if (parameters.keyExists<int>("maxiter")) {
_maxIter = parameters.get<int>("maxiter");
}
if (parameters.keyExists<double>("threshold")) {
_gradThreshold = parameters.get<double>("threshold");
}
if (parameters.keyExists<double>("print-threshold")) {
_printThreshold = parameters.get<double>("print-threshold");
}
if (parameters.keyExists<bool>("print-operators")) {
_printOps = parameters.get<bool>("print-operators");
}
pool = parameters.getString("pool");
// Check if Observable is Fermion or Pauli and manipulate accordingly
//
// if string has ^, it's FermionOperator
if (observable->toString().find("^") != std::string::npos){
auto jw = xacc::getService<ObservableTransform>("jw");
if (std::dynamic_pointer_cast<FermionOperator>(observable)) {
observable = jw->transform(observable);
} else {
auto fermionObservable = xacc::quantum::getObservable("fermion", observable->toString());
observable = jw->transform(std::dynamic_pointer_cast<Observable>(fermionObservable));
}
// observable is PauliOperator, but does not cast down to it
// Not sure about the likelihood of this happening, but want to cover all bases
} else if (observable->toString().find("X") != std::string::npos
|| observable->toString().find("Y") != std::string::npos
|| observable->toString().find("Z") != std::string::npos
&& !std::dynamic_pointer_cast<PauliOperator>(observable)){
auto pauliObservable = xacc::quantum::getObservable("pauli", observable->toString());
observable = std::dynamic_pointer_cast<Observable>(pauliObservable);
}
if (parameters.keyExists<std::vector<double>>("checkpoint-parameters")) {
checkpointParams= parameters.get<std::vector<double>>("checkpoint-parameters");
}
if (parameters.keyExists<std::vector<int>>("checkpoint-ops")) {
checkpointOps = parameters.get<std::vector<int>>("checkpoint-ops");
}
if (parameters.stringExists("gradient-strategy")) {
gradStrategyName = parameters.getString("gradient-strategy");
}
return true;
}
const std::vector<std::string> ADAPT_VQE::requiredParameters() const {
return {"observable", "optimizer", "accelerator", "nElectrons", "pool"};
}
void ADAPT_VQE::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const {
auto ansatzRegistry = xacc::getIRProvider("quantum");
auto ansatzInstructions = ansatzRegistry->createComposite("ansatzCircuit");
auto operatorPool = xacc::getService<OperatorPool>(pool);
auto operators = operatorPool->generate(buffer->size(), nElectrons);
std::vector<std::shared_ptr<Observable>> pauliOps;
std::vector<int> ansatzOps;
auto jw = xacc::getService<ObservableTransform>("jw");
std::stringstream ss;
// Mean-field state
std::size_t j;
for (int i = 0; i < nElectrons/2; i++) {
j = (std::size_t) i;
auto alphaXGate = ansatzRegistry->createInstruction("X", std::vector<std::size_t>{j});
ansatzInstructions->addInstruction(alphaXGate);
j = (std::size_t) (i + buffer->size()/2);
auto betaXGate = ansatzRegistry->createInstruction("X", std::vector<std::size_t>{j});
ansatzInstructions->addInstruction(betaXGate);
}
// Vector of commutators, need to compute them only once
std::vector<std::shared_ptr<Observable>> commutators;
for (auto op : operators){
std::shared_ptr<Observable> pauliOp;
if(std::dynamic_pointer_cast<PauliOperator>(op)){
pauliOp = op;
} else {
pauliOp = jw->transform(op);
}
pauliOps.push_back(pauliOp);
commutators.push_back(observable->commutator(pauliOp));
}
int initialIter = 0;
double oldEnergy = 0.0;
std::vector<double> x; // these are the variational parameters
// Resume from a previously optimized ansatz
if (!checkpointOps.empty()){
if (!checkpointParams.empty()){
x = checkpointParams;
} else {
x.resize(checkpointOps.size());
}
initialIter = checkpointOps.size();
for (int i = 0; i < checkpointOps.size(); i++){
auto exp_i_theta = std::dynamic_pointer_cast<quantum::Circuit>(
xacc::getService<Instruction>("exp_i_theta"));
exp_i_theta->expand(
{std::make_pair("pauli", pauliOps[checkpointOps[i]]->toString()),
std::make_pair("param_id", std::string("x") + std::to_string(i)),
std::make_pair("no-i", true)});
ansatzInstructions->addVariable(std::string("x") + std::to_string(i));
for (auto& inst : exp_i_theta->getInstructions()){
ansatzInstructions->addInstruction(inst);
}
}
auto newOptimizer = xacc::getOptimizer(optimizer->name(),
{std::make_pair(optimizer->name() + "-optimizer", optimizer->get_algorithm()),
std::make_pair("initial-parameters", x)});
auto init_vqe = xacc::getAlgorithm(
"vqe", {std::make_pair("observable", observable),
std::make_pair("optimizer", optimizer),
std::make_pair("accelerator", accelerator),
std::make_pair("ansatz", ansatzInstructions)});
auto tmp_buffer = xacc::qalloc(buffer->size());
oldEnergy = init_vqe->execute(tmp_buffer, x)[0];
ss << std::setprecision(12) << oldEnergy << "\n";
xacc::info(ss.str());
ss.str(std::string());
}
xacc::info("Operator pool: " + operatorPool->name() + "\n");
xacc::info("Number of operators in the pool: " + std::to_string(pauliOps.size()) + "\n\n");
// start ADAPT loop
for (int iter = initialIter; iter < _maxIter; iter++){
xacc::info("Iteration: " + std::to_string(iter + 1) + "\n");
xacc::info("Computing [H, A]\n\n");
xacc::info("Printing commutators with absolute value above " + std::to_string(_printThreshold) + "\n");
int maxCommutatorIdx = 0;
double maxCommutator = 0.0;
double gradientNorm = 0.0;
// Loop over non-vanishing commutators and select the one with largest magnitude
for (int operatorIdx = 0; operatorIdx < commutators.size(); operatorIdx++){
// only compute commutators if they aren't zero
int nTermsCommutator = std::dynamic_pointer_cast<PauliOperator>(commutators[operatorIdx])->getTerms().size();
if(nTermsCommutator != 0){
// Print number of instructions for computing <observable>
xacc::info("Number of instructions for commutator calculation: "
+ std::to_string(nTermsCommutator));
// observe the commutators with the updated circuit ansatz
auto grad_vqe = xacc::getAlgorithm(
"vqe", {std::make_pair("observable", commutators[operatorIdx]),
std::make_pair("optimizer", optimizer),
std::make_pair("accelerator", accelerator),
std::make_pair("ansatz", ansatzInstructions)});
auto tmp_buffer = xacc::qalloc(buffer->size());
auto commutatorValue = std::real(grad_vqe->execute(tmp_buffer, x)[0]);
if(abs(commutatorValue) > _printThreshold){
ss << std::setprecision(12) << "[H," << operatorIdx << "] = " << commutatorValue << "\n";
xacc::info(ss.str());
ss.str(std::string());
}
// update maxCommutator
if(abs(commutatorValue) > abs(maxCommutator)){
maxCommutatorIdx = operatorIdx;
maxCommutator = commutatorValue;
}
gradientNorm += commutatorValue * commutatorValue;
}
}
ss << std::setprecision(12) << "Max gradient component: [H, "
<< maxCommutatorIdx << "] = " << maxCommutator << " a.u.\n";
xacc::info(ss.str());
ss.str(std::string());
gradientNorm = std::sqrt(gradientNorm);
ss << std::setprecision(12) << "Norm of gradient vector: " << gradientNorm << " a.u.\n";
xacc::info(ss.str());
ss.str(std::string());
if (gradientNorm < _gradThreshold) { // ADAPT-VQE converged
xacc::info("\nADAPT-VQE converged in " + std::to_string(iter) + " iterations.\n");
ss << std::setprecision(12) << "ADAPT-VQE energy: " << oldEnergy << " a.u.\n";
xacc::info(ss.str());
ss.str(std::string());
ss << "Optimal parameters: \n";
for (auto param : x){
ss << std::setprecision(12) << param << " ";
}
xacc::info(ss.str() + "\n");
ss.str(std::string());
xacc::info("Final ADAPT-VQE circuit\n" + ansatzInstructions->toString());
return;
} else if (iter < _maxIter) { // Add operator and reoptimize
xacc::info("\nVQE optimization of current ansatz.\n\n");
// keep track of growing ansatz
ansatzOps.push_back(maxCommutatorIdx);
// Instruction service for the operator to be added to the ansatz
auto maxCommutatorGate = std::dynamic_pointer_cast<quantum::Circuit>(
xacc::getService<Instruction>("exp_i_theta"));
// Create instruction for new operator
maxCommutatorGate->expand(
{std::make_pair("pauli", pauliOps[maxCommutatorIdx]->toString()),
std::make_pair("param_id", std::string("x") + std::to_string(iter)),
std::make_pair("no-i", true)});
// Add instruction for new operator to the current ansatz
ansatzInstructions->addVariable(std::string("x") + std::to_string(iter));
// Append new instructions to current circuit
for (auto& inst : maxCommutatorGate->getInstructions()){
ansatzInstructions->addInstruction(inst);
}