qaoa: var assignment mode

quantum/plugins/algorithms/qaoa/qaoa.cpp

@@ -9,6 +9,7 @@
  *
  * Contributors:
  *   Thien Nguyen - initial API and implementation
+ *   Milos Prokop - variable assignment mode
  *******************************************************************************/
 
 #include "qaoa.hpp"
@@ -61,6 +62,20 @@ bool QAOA::initialize(const HeterogeneousMap &parameters) {
     m_parameterizedMode = parameters.getString("parameter-scheme");
   }
 
+  // Determine the optimal variable assignment of the QUBO problem
+  // If false, only optimal value is returned
+  m_varAssignmentMode = false;
+  if (parameters.keyExists<bool>("calc-var-assignment")) {
+	  m_varAssignmentMode = parameters.get<bool>("calc-var-assignment");
+  }
+
+  // Number of samples generated to find the optimal problem solution
+  if(m_varAssignmentMode){
+	nbSamples = 1024;
+	if(parameters.keyExists<int>("nbSamples"))
+		nbSamples = parameters.get<int>("nbSamples");
+  }
+
   if (initializeOk) {
     m_costHamObs = parameters.getPointerLike<Observable>("observable");
     m_qpu = parameters.getPointerLike<Accelerator>("accelerator");
@@ -125,6 +140,53 @@ const std::vector<std::string> QAOA::requiredParameters() const {
   return {"accelerator", "optimizer", "observable"};
 }
 
+//evaluate the energy of measurement
+double QAOA::evaluate_assignment(xacc::Observable* const observable, std::string measurement) const{
+
+	double result =	observable->getIdentitySubTerm()->coefficient().real();
+	for(auto &term: observable->getNonIdentitySubTerms()){
+
+		//get the real part, imag expected to be 0
+		double coeff = term->coefficient().real();
+
+		//parse to get hamiltonian terms
+		std::string term_str = term->toString();
+		std::vector<int> qubit_indices;
+
+		char z_coeff[6]; // Number of ciphers in qubit index are expected to fit here.
+		int z_index = 0;
+
+		int term_i = 0;
+
+		for(size_t i = term_str.length()-1; i > 0; --i){
+			char c = term_str[i];
+			if(isdigit(c))
+				z_coeff[5-z_index++] = c;
+			else if(c == 'Z'){
+				int val = 0;
+				for(; z_index>0; --z_index){
+					val += pow(10, z_index-1) * (z_coeff[5-z_index+1] - '0');
+				}
+				qubit_indices.push_back(val);
+				z_index = 0;
+			}
+			else if(c == ' '){
+				if(++term_i == 2)
+					break;
+			}
+			else{
+				break;
+			}
+		}
+		int multiplier = 1;
+		for(auto &qubit_index: qubit_indices)
+			multiplier *= (measurement[qubit_index] == '0') ? 1 : -1;
+		result += multiplier * coeff;
+	}
+	return result;
+}
+
+
 void QAOA::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const {
   const int nbQubits = buffer->size();
   // we need this for ADAPT-QAOA (Daniel)
@@ -278,8 +340,80 @@ void QAOA::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const {
   // 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));
+
+  if(m_varAssignmentMode){
+
+ 	  typedef std::pair<int, double> meas_freq_eval;
+ 	  typedef std::pair<std::string, meas_freq_eval> measurement;
+
+ 	  auto provider = xacc::getIRProvider("quantum");
+ 	  auto evaled = kernel->operator()(result.second);
+ 	  m_qpu->updateShotsNumber(nbSamples);
+
+ 	  //comment below and uncomment commented to use multiple measurement gates
+ 	  for(size_t i=0; i < nbQubits; ++i){
+ 	  	evaled->addInstructions({provider->createInstruction("Measure", i)});
+ 	  }
+
+ 	  //std::vector<size_t> indices;
+ 	  //for(size_t i=0; i < nbQubits; ++i){
+ 	  //   indices.push_back(i);
+ 	  //}
+ 	  //auto meas = provider->createInstruction("Measure", indices);
+ 	  //evaled->addInstructions({meas});
+
+ 	  m_qpu->execute(buffer, evaled);
+ 	  std::vector<measurement> measurements;
+
+ 	  for(std::pair<std::string, int> meas : buffer->getMeasurementCounts()){
+
+ 		  bool found = false;
+ 		  for(auto &instance : measurements)
+ 			  if(instance.first == meas.first){
+ 				  found = true;
+ 				  break;
+ 			  }
+
+ 		  if(!found){
+ 			  measurements.push_back(measurement(meas.first, // bit string
+ 					  	  	  	  	  meas_freq_eval(meas.second, //number of occurences
+ 					  	  	  	  			  evaluate_assignment(m_costHamObs, meas.first))));
+ 		  }
+ 	  }
+
+ 	 //sort by number of occurences
+ 	  if(m_maximize){
+ 		  std::sort(measurements.begin(), measurements.end(),
+ 			  [](const std::pair<std::string, std::pair<int, double>>& a, const std::pair<std::string, std::pair<int,int>>& b) {
+ 				  return a.second.second > b.second.second;
+ 	      });
+ 	  } else {
+ 		  std::sort(measurements.begin(), measurements.end(),
+ 			  [](const std::pair<std::string, std::pair<int, double>>& a, const std::pair<std::string, std::pair<int,int>>& b) {
+ 		  	  	  return a.second.second < b.second.second;
+ 		  });
+ 	  }
+
+ 	  double optimal_value = measurements[0].second.second;
+ 	  std::string opt_config = measurements[0].first;
+
+ 	  int i = 0;
+ 	  int hits;
+ 	  while(i < measurements.size() && abs(measurements[i++].second.second - optimal_value)<10e-1){
+ 		  hits += measurements[i-1].second.first;
+ 	  }
+
+ 	  double hit_rate = hits / double(nbSamples);
+
+ 	  buffer->addExtraInfo("opt-val", ExtraInfo(optimal_value));
+ 	  buffer->addExtraInfo("opt-config", opt_config);
+ 	  buffer->addExtraInfo("hit_rate: ", ExtraInfo(hit_rate));
+ 	  buffer->addExtraInfo("opt-params", ExtraInfo(result.second));
+
+   }else{
+ 	  buffer->addExtraInfo("opt-val", ExtraInfo(finalCost));
+ 	  buffer->addExtraInfo("opt-params", ExtraInfo(result.second));
+   }
 }
 
 std::vector<double>
@@ -368,4 +502,4 @@ QAOA::execute(const std::shared_ptr<AcceleratorBuffer> buffer,
 }
 
 } // namespace algorithm
-} // namespace xacc
+} // namespace xacc

quantum/plugins/algorithms/qaoa/qaoa.hpp

@@ -38,9 +38,14 @@ class QAOA : public Algorithm
     std::shared_ptr<CompositeInstruction> externalAnsatz;
     std::shared_ptr<CompositeInstruction> m_single_exec_kernel;
     int m_nbSteps;
+    int nbSamples = 1024;
     std::string m_parameterizedMode;
     bool m_maximize = false;
+    bool m_varAssignmentMode = false;
+    bool m_simplifiedSimulationMode = false;
     CompositeInstruction* m_initial_state;
+
+    double evaluate_assignment(xacc::Observable* const observable, std::string measurement) const;
 };
 } // namespace algorithm
 } // namespace xacc

quantum/plugins/algorithms/qaoa/tests/QAOATester.cpp

@@ -9,6 +9,7 @@
  *
  * Contributors:
  *   Thien Nguyen - initial API and implementation
+ *   Milos Prokop - test evaluate_assignment
  *******************************************************************************/
 #include <random>
 #include <gtest/gtest.h>
@@ -243,6 +244,88 @@ TEST(QAOATester, checkMaxCut) {
   // EXPECT_NEAR((*buffer)["opt-val"].as<double>(), 2.0, 1e-3);
 }
 
+TEST(QAOATester, check_evaluate_assignment) {
+
+	//
+	// Tests the QAOA::evaluate_assignment function. The parameters of QAOA are deliberately very UNperformant such that the resulting state will not be far from initial state
+	// and hence after many shots we expect to get all possible assignments.
+	//
+
+	const int num_qubits = 3;
+	const int num_terms = 7;
+	std::string hamiltonian[num_terms][3] = {{ "-5",   "",   ""},
+											 {  "2", "Z0",   ""},
+											 {  "3", "Z1",   ""},
+											 { "-1", "Z2",   ""},
+											 { "-4", "Z0", "Z1"},
+											 {  "7", "Z1", "Z2"},
+											 {"0.4", "Z0", "Z2"}
+								   	   	    };
+
+	std::string hamiltonian_str = hamiltonian[0][0] + " ";
+	for(size_t i = 1; i < num_terms; ++i){
+
+		hamiltonian_str += "+ ";
+
+		for(size_t i2 = 0; i2 < 3; i2++)
+			hamiltonian_str += hamiltonian[i][i2] + " ";
+
+	}
+
+	auto acc = xacc::getAccelerator("qpp", {std::make_pair("shots", 50)});
+	auto buffer = xacc::qalloc(num_qubits);
+	auto observable = xacc::quantum::getObservable("pauli", hamiltonian_str);
+
+	const int nbParams = observable -> getNonIdentitySubTerms().size() + observable->nBits();;
+	std::vector<double> initialParams;
+
+	// Init random parameters
+	for (int i = 0; i < nbParams; ++i)
+	{
+	  initialParams.emplace_back(0);
+	}
+
+	auto optimizer = xacc::getOptimizer("nlopt",
+	  xacc::HeterogeneousMap {
+		 std::make_pair("initial-parameters", initialParams),
+		 std::make_pair("maximize", false),
+		 std::make_pair("nlopt-maxeval", 1) });
+
+		auto qaoa = xacc::getService<xacc::Algorithm>("QAOA");
+		qaoa->initialize({
+		   std::make_pair("accelerator", acc),
+		   std::make_pair("optimizer", optimizer),
+		   std::make_pair("observable", observable),
+		   std::make_pair("steps", 1),
+		   std::make_pair("calc-var-assignment", true),
+		   std::make_pair("nbSamples", 1)
+		});
+
+		qaoa->execute(buffer);
+		std::string meas = (*buffer)["opt-config"].as<std::string>();
+
+		//calculate the expected result
+		double expected_result = 0.;
+		for(size_t i = 0; i < num_terms; ++i){
+
+			auto term = hamiltonian[i];
+
+			if(term[1] == ""){ //identity term
+				expected_result += std::stod(term[0]);
+			}else if(term[2] == ""){ //single-Z term
+				int qbit = std::stoi(term[1].substr(1));
+				expected_result += (meas[qbit] == '1' ? -1 : 1) * std::stod(term[0]);
+			}else{ //double-Z term
+				int qbit1 = std::stoi(term[1].substr(1));
+				int qbit2 = std::stoi(term[2].substr(1));
+
+				expected_result += (meas[qbit1] == '1' ? -1 : 1) * (meas[qbit2] == '1' ? -1 : 1) * std::stod(term[0]);
+			}
+		}
+
+		EXPECT_NEAR((*buffer)["opt-val"].as<double>(), expected_result, 1e-1);
+}
+
 int main(int argc, char **argv) {
   xacc::Initialize(argc, argv);
   ::testing::InitGoogleTest(&argc, argv);

xacc/accelerator/Accelerator.hpp

@@ -68,6 +68,8 @@ class Accelerator : public Identifiable {
   virtual void updateConfiguration(const HeterogeneousMap &&config) {
     updateConfiguration(config);
   }
+  virtual void updateShotsNumber(int shots) {nbShots = shots;}
+
   virtual const std::vector<std::string> configurationKeys() = 0;
 
   virtual HeterogeneousMap getProperties() { return HeterogeneousMap(); }
@@ -134,6 +136,10 @@ class Accelerator : public Identifiable {
   }
 
   virtual ~Accelerator() {}
+
+protected:
+  int nbShots = -1;
+
 };
 
 template Accelerator* HeterogeneousMap::getPointerLike<Accelerator>(const std::string key) const;