Implemented Scipy optimizer plugin

Signed-off-by: Daniel Claudino <claudinodc@ornl.gov>

quantum/plugins/optimizers/CMakeLists.txt

@@ -1,2 +1,3 @@
 add_subdirectory(nlopt-optimizers)
-add_subdirectory(mlpack)
+add_subdirectory(mlpack)
+add_subdirectory(scipy)

quantum/plugins/optimizers/scipy/CMakeLists.txt

@@ -0,0 +1,58 @@
+message(STATUS "${BoldGreen}Building Scipy Optimizer.${ColorReset}")
+set(LIBRARY_NAME xacc-scipy-optimizer)
+
+file(GLOB
+     SRC
+     scipy_optimizer.cpp)
+
+usfunctiongetresourcesource(TARGET
+                            ${LIBRARY_NAME}
+                            OUT
+                            SRC)
+usfunctiongeneratebundleinit(TARGET
+                             ${LIBRARY_NAME}
+                             OUT
+                             SRC)
+#find_package(pybind11 REQUIRED)
+find_package(Python COMPONENTS Interpreter Development)
+
+add_library(${LIBRARY_NAME} SHARED ${SRC})
+
+if(Python_FOUND)
+  message("FOUND python")
+endif()
+target_include_directories(${LIBRARY_NAME}
+                           PUBLIC . ${CMAKE_SOURCE_DIR}/tpls/pybind11/include ${Python_INCLUDE_DIRS})
+
+                           target_link_libraries(${LIBRARY_NAME} PUBLIC xacc Python::Python)
+
+set(_bundle_name xacc_optimizer_scipy)
+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)

quantum/plugins/optimizers/scipy/manifest.json

@@ -0,0 +1,6 @@
+{
+  "bundle.symbolic_name" : "xacc_optimizer_scipy",
+  "bundle.activator" : true,
+  "bundle.name" : "XACC Scipy Optimizer",
+  "bundle.description" : ""
+}

quantum/plugins/optimizers/scipy/scipy_optimizer.cpp

@@ -0,0 +1,162 @@
+#include "scipy_optimizer.hpp"
+#include "Optimizer.hpp"
+#include "pybind11/stl.h"
+#include "pybind11/numpy.h"
+#include "xacc.hpp"
+#include "xacc_plugin.hpp"
+
+namespace py = pybind11;
+
+namespace xacc {
+
+const std::string ScipyOptimizer::get_algorithm() const {
+  std::string optimizerAlgo = "COBYLA";
+  if (options.stringExists("algorithm")) {
+    optimizerAlgo = options.getString("algorithm");
+  }
+  if (options.stringExists("scipy-optimizer")) {
+    optimizerAlgo = options.getString("scipy-optimizer");
+  }
+  return optimizerAlgo;
+}
+
+const bool ScipyOptimizer::isGradientBased() const {
+
+  std::string optimizerAlgo = "cobyla";
+  if (options.stringExists("algorithm")) {
+    optimizerAlgo = options.getString("algorithm");
+  }
+  if (options.stringExists("scipy-optimizer")) {
+    optimizerAlgo = options.getString("scipy-optimizer");
+  }
+
+  if (options.stringExists("optimizer")) {
+    optimizerAlgo = options.getString("optimizer");
+  }
+
+  if (optimizerAlgo == "bfgs") {
+    return true;
+  } else {
+    return false;
+  }
+}
+
+OptResult ScipyOptimizer::optimize(OptFunction& function) {
+
+  bool maximize = false;
+  if (options.keyExists<bool>("maximize")) {
+      xacc::info("Turning on maximize!");
+    maximize = options.get<bool>("maximize");
+  }
+
+  std::string algo = "COBYLA";
+  if (options.stringExists("algorithm")) {
+    algo = options.getString("algorithm");
+  }
+  if (options.stringExists("scipy-optimizer")) {
+    algo = options.getString("scipy-optimizer");
+  }
+  if (options.stringExists("optimizer")) {
+    algo = options.getString("optimizer");
+  }
+
+    if (algo == "cobyla" || algo == "COBYLA") {
+      algo ="COBYLA";
+    } else if (algo == "nelder-mead" || algo == "Nelder-Mead") {
+      algo = "Nelder-Mead";
+    } else if (algo == "bfgs" || algo == "BFGS" || algo == "l-bfgs") {
+      algo = "BFGS";
+    } else {
+      xacc::XACCLogger::instance()->error("Invalid optimizer at this time: " +
+                                          algo);
+    }
+
+  double tol = 1e-6;
+  if (options.keyExists<double>("ftol")) {
+    tol = options.get<double>("ftol");
+    xacc::info("[Scipy] function tolerance set to " + std::to_string(tol));
+  }
+  if (options.keyExists<double>("scipy-ftol")) {
+    tol = options.get<double>("ftol");
+    xacc::info("[Scipy] function tolerance set to " + std::to_string(tol));
+  }
+
+  int maxeval = 1000;
+  if (options.keyExists<int>("maxeval")) {
+    maxeval = options.get<int>("maxeval");
+    xacc::info("[Scipy] max function evaluations set to " +
+               std::to_string(maxeval));
+  }
+  if (options.keyExists<int>("scipy-maxeval")) {
+    maxeval = options.get<int>("maxeval");
+    xacc::info("[Scipy] max function evaluations set to " +
+               std::to_string(maxeval));
+  }
+
+  std::vector<double> x(function.dimensions());
+  if (options.keyExists<std::vector<double>>("initial-parameters")) {
+    x = options.get_with_throw<std::vector<double>>("initial-parameters");
+  } else if (options.keyExists<std::vector<int>>("initial-parameters")) {
+    auto tmpx = options.get<std::vector<int>>("initial-parameters");
+    x = std::vector<double>(tmpx.begin(), tmpx.end());
+  }
+
+  // here the python stuff starts
+  py::list pyInitialParams;
+  for (const auto &param : x) {
+      pyInitialParams.append(param);
+  }
+
+  if (isGradientBased()) std::cout << algo << "\n";
+
+  // wrap the objective function in this lambda
+  // scipy passes a numpy array to this function, hence the py::array_t type
+  py::object pyObjFunction = py::cpp_function([&function](const py::array_t<double>& pyParams) {
+    std::vector<double> params(pyParams.size());
+    std::memcpy(params.data(), pyParams.data(), pyParams.size() * sizeof(double));
+    return function(std::move(params));
+  });
+
+  // call this for gradient-based optimization
+  py::object pyObjFunctionWithGrad = py::cpp_function([&function](const py::array_t<double>& pyParams) {
+    std::vector<double> params(pyParams.size());
+    std::memcpy(params.data(), pyParams.data(), pyParams.size() * sizeof(double));
+
+    std::vector<double> grad(params.size());
+    double result = function(params, grad);
+    py::array_t<double> pyGrad(grad.size());
+    std::memcpy(pyGrad.mutable_data(), grad.data(), grad.size() * sizeof(double));
+
+    return py::make_tuple(result, pyGrad);
+  });
+
+  py::module scipy_optimize = py::module::import("scipy.optimize");
+
+  // error handling helps here to see if it's coming from C++ or python
+  try {
+
+    py::object result = scipy_optimize.attr("minimize")(
+        isGradientBased() ? pyObjFunctionWithGrad : pyObjFunction,
+        pyInitialParams,
+        py::arg("args") = py::tuple(),
+        py::arg("method") = algo,
+        py::arg("tol") = tol,
+        py::arg("jac") = (isGradientBased() ? true : false)
+    );
+
+    std::vector<double> optimizedParams = result.attr("x").cast<std::vector<double>>();
+    double optimalValue = result.attr("fun").cast<double>();
+
+    return {optimalValue, optimizedParams};
+    } catch (const py::error_already_set& e) {
+      std::cerr << "Python error: " << e.what() << std::endl;
+      throw;
+    } catch (const std::exception& e) {
+      std::cerr << "Error: " << e.what() << std::endl;
+      throw;
+    }
+}
+} // namespace xacc
+
+// Register the plugin with XACC
+REGISTER_OPTIMIZER(xacc::ScipyOptimizer)

quantum/plugins/optimizers/scipy/scipy_optimizer.hpp

@@ -0,0 +1,26 @@
+#ifndef SCIPYOPTIMIZER_HPP
+#define SCIPYOPTIMIZER_HPP
+
+#include <xacc.hpp>
+#include <xacc_service.hpp>
+#include <Optimizer.hpp>
+#include <pybind11/embed.h>
+
+namespace xacc {
+
+class ScipyOptimizer : public xacc::Optimizer {
+public:
+
+  ScipyOptimizer() = default;
+  ~ScipyOptimizer() = default;
+
+  const std::string name() const override { return "scipy"; }
+  const std::string description() const override { return ""; }
+
+  OptResult optimize(OptFunction &function) override;
+  const bool isGradientBased() const override;
+  virtual const std::string get_algorithm() const override;
+};
+
+} // namespace xacc
+#endif // SCIPYOPTIMIZER_HPP

quantum/plugins/optimizers/scipy/tests/CMakeLists.txt

@@ -0,0 +1,4 @@
+include_directories(${CMAKE_SOURCE_DIR}/tpls/pybind11/include ${Python_INCLUDE_DIRS})
+
+add_xacc_test(ScipyOptimizer)
+target_link_libraries(ScipyOptimizerTester xacc Python::Python)

quantum/plugins/optimizers/scipy/tests/ScipyOptimizerTester.cpp

@@ -0,0 +1,103 @@
+#include <gtest/gtest.h>
+#include "xacc.hpp"
+#include "xacc_service.hpp"
+#include <dlfcn.h>
+#include <pybind11/embed.h>
+
+using namespace xacc;
+namespace py = pybind11;
+
+TEST(ScipyOptimizerTester, checkSimple) {
+
+
+  auto optimizer =
+      xacc::getService<Optimizer>("scipy"); // NLOptimizer optimizer;
+
+  OptFunction f([](const std::vector<double> &x, std::vector<double>& g) { return x[0] * x[0] + 5; },
+                1);
+
+  EXPECT_EQ(optimizer->name(), "scipy");
+  EXPECT_EQ(1, f.dimensions());
+
+
+  optimizer->setOptions(HeterogeneousMap{std::make_pair("maxeval", 20)});
+
+  auto result = optimizer->optimize(f);
+  EXPECT_NEAR(5.0, result.first, 1.0e-6);
+  EXPECT_NEAR(result.second[0], 0.0, 1.0e-6);
+
+}
+
+TEST(ScipyOptimizerTester, checkGradient) {
+
+  auto optimizer =
+      xacc::getService<Optimizer>("scipy");
+
+  OptFunction f(
+      [](const std::vector<double> &x, std::vector<double> &grad) {
+        if (!grad.empty()) {
+          std::cout << "GRAD\n";
+          grad[0] = 2. * x[0];
+        }
+        auto xx = x[0] * x[0] + 5;
+        std::cout << xx << "\n";
+        return xx;
+      },
+      1);
+
+  EXPECT_EQ(1, f.dimensions());
+
+
+  optimizer->setOptions(
+      HeterogeneousMap{std::make_pair("maxeval", 20),std::make_pair("initial-parameters", std::vector<double>{1.0}),
+                       std::make_pair("optimizer", "bfgs")});
+
+  auto result = optimizer->optimize(f);
+
+  EXPECT_NEAR(result.first, 5.0, 1e-4);
+  EXPECT_NEAR(result.second[0], 0.0, 1e-4);
+}
+
+TEST(ScipyOptimizerTester, checkGradientRosenbrock) {
+
+  auto optimizer =
+      xacc::getService<Optimizer>("scipy");
+
+  OptFunction f(
+      [](const std::vector<double> &x, std::vector<double> &grad) {
+        if (!grad.empty()) {
+        //   std::cout << "GRAD\n";
+          grad[0] = -2 * (1 - x[0]) + 400 * (std::pow(x[0], 3) - x[1] * x[0]);
+          grad[1] = 200 * (x[1] - std::pow(x[0],2));
+        }
+        auto xx = 100 * std::pow(x[1] - std::pow(x[0], 2), 2) + std::pow(1 - x[0], 2);
+        std::cout << xx << ", " << x << ", " << grad << "\n";
+
+        return xx;
+      },
+      2);
+
+  EXPECT_EQ(2, f.dimensions());
+
+  optimizer->setOptions(
+      HeterogeneousMap{std::make_pair("maxeval", 200),
+                       std::make_pair("optimizer", "bfgs")});
+
+  auto result = optimizer->optimize(f);
+
+  EXPECT_NEAR(result.first, 0.0, 1e-4);
+  EXPECT_NEAR(result.second[0], 1.0, 1e-4);
+  EXPECT_NEAR(result.second[1], 1.0, 1e-4);
+
+}
+
+int main(int argc, char **argv) {
+  dlopen("libpython3.8.so", RTLD_LAZY | RTLD_GLOBAL);
+  xacc::Initialize(argc, argv);
+  py::initialize_interpreter();
+  ::testing::InitGoogleTest(&argc, argv);
+  auto ret = RUN_ALL_TESTS();
+  py::finalize_interpreter();
+  xacc::Finalize();
+  return ret;
+}