Commit 187238d8 authored by Mccaskey, Alex's avatar Mccaskey, Alex
Browse files

adding uccsd ir generator


Signed-off-by: Mccaskey, Alex's avatarAlex McCaskey <mccaskeyaj@ornl.gov>
parent 11b667fa
Pipeline #66310 failed with stage
in 6 minutes and 34 seconds
import xacc
opstr = '''(0.174073,0) Z2 Z3 +
(0.1202,0) Z1 Z3 +
(0.165607,0) Z1 Z2 +
(0.165607,0) Z0 Z3 +
(0.1202,0) Z0 Z2 +
(-0.0454063,0) Y0 Y1 X2 X3 +
(-0.220041,0) Z3 +
(-0.106477,0) +
(0.17028,0) Z0 +
(-0.220041,0) Z2 +
(0.17028,0) Z1 +
(-0.0454063,0) X0 X1 Y2 Y3 +
(0.0454063,0) X0 Y1 Y2 X3 +
(0.168336,0) Z0 Z1 +
(0.0454063,0) Y0 X1 X2 Y3'''
op = xacc.getObservable('pauli',opstr)
qpu = xacc.getAccelerator('tnqvm')
@xacc.qpu(algo='vqe', accelerator=qpu, observable=op)
def ansatz_vqe(buffer, *args):
uccsd(nelectrons=2,nqubits=4)
buffer = xacc.qalloc(4)
ansatz_vqe(buffer, 0,0)
\ No newline at end of file
......@@ -69,6 +69,14 @@ FermionOperator::FermionOperator(Operators operators,
FermionOperator::FermionOperator(Operators operators, double coeff)
: FermionOperator(operators, std::complex<double>(coeff, 0)) {}
FermionOperator::FermionOperator(Operators operators, double coeff, std::string var)
{
terms.emplace(std::piecewise_construct,
std::forward_as_tuple(FermionTerm::id(operators)),
std::forward_as_tuple(std::complex<double>(coeff,0.0), operators, var));
}
void FermionOperator::clear() { terms.clear(); }
std::vector<std::shared_ptr<Function>>
......@@ -81,7 +89,7 @@ const std::string FermionOperator::toString() {
std::stringstream s;
for (auto &kv : terms) {
std::complex<double> c = std::get<0>(kv.second);
s << c << " ";
s << c << " " << std::get<2>(kv.second) << " ";
Operators ops = std::get<1>(kv.second);
std::vector<int> creations, annhilations;
......@@ -157,7 +165,6 @@ const int FermionOperator::nBits() {
FermionOperator &
FermionOperator::operator+=(const FermionOperator &v) noexcept {
FermionOperator vv = v;
// std::cout << "adding " << toString() << " and " << vv.toString() << "\n";
for (auto &kv : v.terms) {
auto termId = kv.first;
......
......@@ -22,7 +22,7 @@ using Operator = std::pair<int, bool>;
using Operators = std::vector<Operator>;
// using SiteMap = std::map<Operators, std::complex<double>>;
using FermionTermTuple = std::tuple<std::complex<double>, Operators>;
using FermionTermTuple = std::tuple<std::complex<double>, Operators, std::string>;
class FermionTerm : public FermionTermTuple,
public tao::operators::commutative_multipliable<FermionTerm>,
public tao::operators::equality_comparable<FermionTerm> {
......@@ -31,31 +31,43 @@ public:
FermionTerm() {
std::get<0>(*this) = std::complex<double>(1.0, 0.0);
std::get<1>(*this) = {};
std::get<2>(*this) = "";
}
FermionTerm(const FermionTerm &t) {
std::get<0>(*this) = std::get<0>(t);
std::get<1>(*this) = std::get<1>(t);
std::get<2>(*this) = std::get<2>(t);
}
FermionTerm(std::complex<double> c) {
std::get<0>(*this) = c;
std::get<1>(*this) = {};
std::get<2>(*this) = "";
}
FermionTerm(double c) {
std::get<0>(*this) = std::complex<double>(c,0.0);
std::get<1>(*this) = {};
std::get<2>(*this) = "";
}
FermionTerm(std::complex<double> c, Operators ops) {
std::get<0>(*this) = c;
std::get<1>(*this) = ops;
std::get<2>(*this) = "";
}
FermionTerm(std::complex<double> c, Operators ops, std::string var) {
std::get<0>(*this) = c;
std::get<1>(*this) = ops;
std::get<2>(*this) = var;
}
FermionTerm(Operators ops) {
std::get<0>(*this) = std::complex<double>(1.0, 0.0);
std::get<1>(*this) = ops;
std::get<2>(*this) = "";
}
static const std::string id(const Operators &ops) {
......@@ -118,6 +130,7 @@ public:
Operators &ops() { return std::get<1>(*this); }
std::complex<double> &coeff() { return std::get<0>(*this); }
std::string var() {return std::get<2>(*this);}
FermionTerm &operator*=(const FermionTerm &v) noexcept;
......@@ -150,6 +163,7 @@ public:
FermionOperator(Operators operators);
FermionOperator(Operators operators, std::complex<double> coeff);
FermionOperator(Operators operators, double coeff);
FermionOperator(Operators operators, double coeff, std::string var);
std::vector<std::shared_ptr<Function>>
observe(std::shared_ptr<Function> function) override;
......
......@@ -25,10 +25,11 @@ JW::transform(std::shared_ptr<Observable> observable) {
// Loop over all Fermionic terms...
for (auto &kv : terms) {
auto var = kv.second.var();
auto coeff = kv.second.coeff();
auto operators = kv.second.ops();
PauliOperator current(coeff);
PauliOperator current(coeff, var);
for (auto &kv2 : operators) {
std::map<int, std::string> zs;
auto isCreation = kv2.second;
......
......@@ -3,6 +3,7 @@ set(LIBRARY_NAME xacc-generators)
file(GLOB SRC hwe/hwe.cpp
range/range.cpp
exp/exp.cpp
uccsd/uccsd.cpp
GeneratorsActivator.cpp)
usfunctiongetresourcesource(TARGET ${LIBRARY_NAME} OUT SRC)
......@@ -12,7 +13,7 @@ add_library(${LIBRARY_NAME} SHARED ${SRC})
target_include_directories(
${LIBRARY_NAME}
PUBLIC . hwe exp range)
PUBLIC . hwe exp range uccsd)
target_link_libraries(${LIBRARY_NAME} PUBLIC xacc PRIVATE xacc-pauli xacc-fermion)
......@@ -46,6 +47,7 @@ if(XACC_BUILD_TESTS)
add_subdirectory(hwe/tests)
add_subdirectory(exp/tests)
add_subdirectory(range/tests)
add_subdirectory(uccsd/tests)
endif()
install(TARGETS ${LIBRARY_NAME} DESTINATION ${CMAKE_INSTALL_PREFIX}/plugins)
#include "hwe.hpp"
#include "exp.hpp"
#include "range.hpp"
#include "uccsd.hpp"
#include "cppmicroservices/BundleActivator.h"
#include "cppmicroservices/BundleContext.h"
......@@ -16,9 +17,11 @@ public:
auto hwe = std::make_shared<xacc::generators::HWE>();
auto expit = std::make_shared<xacc::generators::Exp>();
auto r = std::make_shared<xacc::generators::Range>();
auto u = std::make_shared<xacc::generators::UCCSD>();
context.RegisterService<xacc::IRGenerator>(hwe);
context.RegisterService<xacc::IRGenerator>(expit);
context.RegisterService<xacc::IRGenerator>(r);
context.RegisterService<xacc::IRGenerator>(u);
}
......
/*
* CommutingSetGenerator.hpp
*
* Created on: Aug 4, 2017
* Author: aqw
*/
#ifndef VQE_TRANSFORMATION_COMMUTINGSETGENERATOR_HPP_
#define VQE_TRANSFORMATION_COMMUTINGSETGENERATOR_HPP_
#include "PauliOperator.hpp"
#include <Eigen/Core>
#include <numeric>
using namespace xacc::quantum;
namespace xacc {
namespace vqe {
class CommutingSetGenerator {
private:
std::pair<Eigen::VectorXi, Eigen::VectorXi> bv(Term& op, int nQubits) {
Eigen::VectorXi vx = Eigen::VectorXi::Zero(nQubits);
Eigen::VectorXi vz = Eigen::VectorXi::Zero(nQubits);
for (auto term : op.ops()) {
if (term.second == "X") {
vx(term.first) += 1;
} else if (term.second == "Z") {
vz(term.first) += 1;
} else if (term.second == "Y") {
vx(term.first) += 1;
vz(term.first) += 1;
}
}
return std::make_pair(vx, vz);
};
int bv_commutator(Term& term1, Term& term2, int nQubits) {
auto pair1 = bv(term1, nQubits);
auto pair2 = bv(term2, nQubits);
auto scalar = pair1.first.dot(pair2.second) + pair1.second.dot(pair2.first);
return scalar % 2;
};
public:
std::vector<std::vector<Term>> getCommutingSet(
PauliOperator& composite, int n_qubits) {
std::vector<std::vector<Term>> commuting_ops;
std::vector<Term> allTerms;
for (auto& kv : composite.getTerms()) {
allTerms.push_back(kv.second);
}
for (int i = 0; i < allTerms.size(); i++) {
auto t_i = allTerms[i];
if (i == 0) {
commuting_ops.push_back({t_i});
} else {
auto comm_ticker = 0;
for (int j = 0; j < commuting_ops.size(); j++) {
auto j_op_list = commuting_ops[j];
int sum = 0;
int innerCounter = 0;
for (auto j_op : j_op_list) {
auto t_jopPtr = allTerms[innerCounter];
sum += bv_commutator(t_i, t_jopPtr,
n_qubits);
innerCounter++;
}
if (sum == 0) {
commuting_ops[j].push_back(t_i);
comm_ticker += 1;
break;
}
}
if (comm_ticker == 0) {
commuting_ops.push_back({t_i});
}
}
}
return commuting_ops;
}
};
}
}
#endif
add_xacc_test(UCCSD)
target_link_libraries(UCCSDTester xacc)
\ No newline at end of file
/***********************************************************************************
* 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
*
**********************************************************************************/
#include <gtest/gtest.h>
// #include "UCCSD.hpp"
#include "Instruction.hpp"
#include "XACC.hpp"
using namespace xacc;
TEST(UCCSDTester,checkUCCSD) {
xacc::Initialize();
auto uccsd = xacc::getIRGenerator("uccsd");
auto f = uccsd->generate({{"nelectrons", 2}, {"nqubits", 4}});
std::cout << f->toString() << "\n";
xacc::Finalize();
}
int main(int argc, char** argv) {
::testing::InitGoogleTest(&argc, argv);
return RUN_ALL_TESTS();
}
#include "uccsd.hpp"
#include "xacc_service.hpp"
#include "FermionOperator.hpp"
#include "PauliOperator.hpp"
#include "XACC.hpp"
#include "ObservableTransform.hpp"
using namespace xacc::quantum;
namespace xacc {
namespace generators {
std::shared_ptr<Function>
UCCSD::generate(std::map<std::string, InstructionParameter> &parameters) {
if (!parameters.count("nelectrons")) {
xacc::error("Invalid mapping of parameters for UCCSD generator, missing "
"nelectrons key.");
}
if (!parameters.count("nqubits")) {
xacc::error("Invalid mapping of parameters for UCCSD generator, missing "
"nqubits key.");
}
std::vector<xacc::InstructionParameter> variables;
if (!parameters.empty()) {
options = parameters;
}
int nQubits = options["nqubits"].as<int>();
int nElectrons = options["nelectrons"].as<int>();
// xacc::info("UCCSD Generator (nqubits,nelectrons) = " +
// std::to_string(nQubits)+", " + std::to_string(nElectrons) +".");
// Compute the number of parameters
auto _nOccupied = (int)std::ceil(nElectrons / 2.0);
auto _nVirtual = nQubits / 2 - _nOccupied;
auto nSingle = _nOccupied * _nVirtual;
auto nDouble = nSingle * (nSingle + 1) / 2;
auto _nParameters = nSingle + nDouble;
std::vector<std::string> params;
if (variables.empty()) {
for (int i = 0; i < _nParameters; i++) {
params.push_back("theta" + std::to_string(i));
variables.push_back(InstructionParameter("theta" + std::to_string(i)));
}
} else {
for (int i = 0; i < _nParameters; i++) {
params.push_back(variables[i].as<std::string>());
}
}
auto slice = [](const std::vector<std::string> &v, int start = 0,
int end = -1) {
int oldlen = v.size();
int newlen;
if (end == -1 or end >= oldlen) {
newlen = oldlen - start;
} else {
newlen = end - start;
}
std::vector<std::string> nv(newlen);
for (int i = 0; i < newlen; i++) {
nv[i] = v[start + i];
}
return nv;
};
auto singleParams = slice(params, 0, nSingle);
auto doubleParams1 = slice(params, nSingle, 2 * nSingle);
auto doubleParams2 = slice(params, 2 * nSingle);
std::vector<std::function<int(int)>> fs{[](int i) { return 2 * i; },
[](int i) { return 2 * i + 1; }};
using OpType = std::vector<std::pair<int, bool>>;
int count = 0;
std::stringstream fOpSS;
FermionOperator myOp;
for (int i = 0; i < _nVirtual; i++) {
for (int j = 0; j < _nOccupied; j++) {
auto vs = _nOccupied + i;
auto os = j;
for (int s = 0; s < 2; s++) {
auto ti = fs[s];
auto oi = fs[1 - s];
auto vt = ti(vs);
auto vo = oi(vs);
auto ot = ti(os);
auto oo = oi(os);
OpType op1{{vt, 1}, {ot, 0}}, op2{{ot, 1}, {vt, 0}};
FermionOperator op(op1, 1.0, singleParams[count]);
FermionOperator opp(op2, -1., singleParams[count]);
OpType op3{{vt, 1}, {ot, 0}, {vo, 1}, {oo, 0}},
op4{{oo, 1}, {vo, 0}, {ot, 1}, {vt, 0}};
FermionOperator oppp(op3, -1., doubleParams1[count]);
FermionOperator opppp(op4, 1., doubleParams1[count]);
myOp += op + opp + oppp + opppp;
}
count++;
}
}
count = 0;
// routine for converting amplitudes for use by UCCSD
std::vector<std::tuple<int, int>> tupleVec;
for (int i = 0; i < _nVirtual; i++) {
for (int j = 0; j < _nOccupied; j++) {
tupleVec.push_back(std::make_tuple(i, j));
}
}
// Combination lambda used to determine indices
auto Combination = [=](std::vector<std::tuple<int, int>> t) {
std::vector<std::tuple<int, int, int, int>> comboVec;
for (int i = 0; i < t.size(); i++) {
for (int j = i + 1; j < t.size(); j++) {
const std::tuple<int, int, int, int> newTuple =
std::tuple_cat(t[i], t[j]);
comboVec.push_back(newTuple);
}
}
return comboVec;
};
auto combineVec = Combination(tupleVec);
for (auto i : combineVec) {
auto p = std::get<0>(i);
auto q = std::get<1>(i);
auto r = std::get<2>(i);
auto s = std::get<3>(i);
auto vs1 = _nOccupied + p;
auto os1 = q;
auto vs2 = _nOccupied + r;
auto os2 = s;
for (int sa = 0; sa < 2; sa++) {
for (int sb = 0; sb < 2; sb++) {
auto ia = fs[sa];
auto ib = fs[sb];
auto v1a = ia(vs1);
auto o1a = ia(os1);
auto v2b = ib(vs2);
auto o2b = ib(os2);
OpType op5{{v1a, 1}, {o1a, 0}, {v2b, 1}, {o2b, 0}},
op6{{o2b, 1}, {o1a, 0}, {v2b, 1}, {o2b, 0}};
FermionOperator op(op5, -1., doubleParams2[count]);
FermionOperator op2(op6, 1., doubleParams2[count]);
myOp += op + op2;
}
}
count++;
}
// std::cout << "HELLO: " << myOp.toString() << "\n";
auto jw = xacc::getService<ObservableTransform>("jw");
auto compositeResult =
jw->transform(std::shared_ptr<Observable>(&myOp, [](Observable *) {}));
// std::cout << "COMP:\n" << compositeResult->toString() << "\n";
// Create the Spin Hamiltonian
auto transformedIR =
std::dynamic_pointer_cast<PauliOperator>(compositeResult)->toXACCIR();
std::unordered_map<std::string, Term> terms =
std::dynamic_pointer_cast<PauliOperator>(compositeResult)->getTerms();
// CommutingSetGenerator gen;
auto commutingSets =
std::dynamic_pointer_cast<PauliOperator>(compositeResult)
->getTerms(); // gen.getCommutingSet(compositeResult, nQubits);
auto pi = 3.14159265358979323; // boost::math::constants::pi<double>();
auto gateRegistry = xacc::getService<IRProvider>("quantum");
auto uccsdGateFunction =
gateRegistry->createFunction("uccsdPrep", {}, variables);
// Perform Trotterization...
for (auto inst : commutingSets) {
// for (auto inst : s) {
Term spinInst = inst.second;
// Get the individual pauli terms
auto termsMap = std::get<2>(spinInst);
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});
}
}
// The largest qubit index is on the last term
int largestQbitIdx = terms[terms.size() - 1].first;
auto tempFunction = gateRegistry->createFunction("temp", {}, {});
for (int i = 0; i < terms.size(); i++) {
auto qbitIdx = terms[i].first;
auto gateName = terms[i].second;
if (i < terms.size() - 1) {
auto cnot = gateRegistry->createInstruction(
"CNOT", std::vector<int>{qbitIdx, terms[i + 1].first});
tempFunction->addInstruction(cnot);
}
if (gateName == "X") {
auto hadamard =
gateRegistry->createInstruction("H", std::vector<int>{qbitIdx});
tempFunction->insertInstruction(0, hadamard);
} else if (gateName == "Y") {
auto rx =
gateRegistry->createInstruction("Rx", std::vector<int>{qbitIdx});
InstructionParameter p(pi / 2.0);
rx->setParameter(0, p);
tempFunction->insertInstruction(0, rx);
}
// Add the Rotation for the last term
if (i == terms.size() - 1) {
// FIXME DONT FORGET DIVIDE BY 2
std::stringstream ss;
ss << 2 * std::imag(std::get<0>(spinInst)) << " * "
<< std::get<1>(spinInst);
auto rz =
gateRegistry->createInstruction("Rz", std::vector<int>{qbitIdx});
InstructionParameter p(ss.str());