add LPSolve interface

CMakeLists.txt

@@ -472,10 +472,10 @@ endif(NEED_LAPACK AND NOT SUITESPARSE_FOUND)
 set(CMAKE_FIND_LIBRARY_PREFIXES "${TMP_CMAKE_FIND_LIBRARY_PREFIXES}")
 
 
-find_package (lpsolve)
+find_package (LPSolve)
 if (TARGET lpsolve::lpsolve )
   set (COMISO_LPSOLVE_CONFIG_FILE_SETTINGS "#define COMISO_LPSOLVE_AVAILABLE 1" )
-  target_link_libraries     (CoMISo PRIVATE lpsolve::lpsolve)
+  target_link_libraries     (CoMISo PUBLIC lpsolve::lpsolve)
 else ()
   message (STATUS "lpsolve not found!")
   set (COMISO_LPSOLVE_CONFIG_FILE_SETTINGS "#define COMISO_LPSOLVE_AVAILABLE 0" )

Examples/small_linear_problem/main.cc

@@ -33,7 +33,7 @@
 #include <CoMISo/NSolver/CPLEXSolver.hh>
 #include <CoMISo/NSolver/GUROBISolver.hh>
 #include <CoMISo/NSolver/CBCSolver.hh>
-
+#include <CoMISo/NSolver/LPSolveSolver.hh>
 
 // minimize linear problem E = 8*x + 2*y + 3*z subject to x+y+z >= 2 and z-y >= 1 and x, y, z binary
 
@@ -84,34 +84,48 @@ int main(void)
   std::cout << "---------- 3) Get CPLEX and optimize... " << std::endl;
   COMISO::CPLEXSolver csol;
   csol.solve(&lp, constraints, dc);
-#endif
 
   std::cout << "---------- 4) Print solution..." << std::endl;
   for( int i=0; i<n; ++i)
-    std::cerr << "x_" << i << " = " << lp.x()[i] << std::endl;
+    std::cout << "x_" << i << " = " << lp.x()[i] << std::endl;
+#endif
 
   // check if GUROBI solver available in current configuration
 #if( COMISO_GUROBI_AVAILABLE)
     std::cout << "---------- 5) Get GUROBI and optimize... " << std::endl;
     COMISO::GUROBISolver gsol;
     gsol.solve(&lp, constraints, dc);
+
+    std::cout << "---------- 6) Print solution..." << std::endl;
+    for( int i=0; i<n; ++i)
+      std::cout << "x_" << i << " = " << lp.x()[i] << std::endl;
 #endif
 
-  std::cout << "---------- 6) Print solution..." << std::endl;
-  for( int i=0; i<n; ++i)
-    std::cerr << "x_" << i << " = " << lp.x()[i] << std::endl;
 
 
-    // check if CBC solver available in current configuration
+  // check if CBC solver available in current configuration
 #if( COMISO_CBC_AVAILABLE)
-    std::cout << "---------- 6) Get CBC and optimize... " << std::endl;
-    COMISO::CBCSolver cbc_sol;
-    cbc_sol.solve(&lp, constraints, dc);
-#endif
+  std::cout << "---------- 7) Get CBC and optimize... " << std::endl;
+  COMISO::CBCSolver cbc_sol;
+  cbc_sol.solve(&lp, constraints, dc);
 
   std::cout << "---------- 8) Print solution..." << std::endl;
   for( int i=0; i<n; ++i)
-    std::cerr << "x_" << i << " = " << lp.x()[i] << std::endl;
+    std::cout << "x_" << i << " = " << lp.x()[i] << std::endl;
+  std::cout << "Objective is " << lp.eval_f(lp.x().data()) << std::endl;
+#endif
+
+
+  // check if CBC solver available in current configuration
+#if( COMISO_LPSOLVE_AVAILABLE)
+  std::cout << "---------- 9) Get lpsolve and optimize... " << std::endl;
+  COMISO::LPSolveSolver lpsolve_solver;
+  lpsolve_solver.solve(&lp, constraints, dc);
+  std::cout << "---------- 10) Print solution..." << std::endl;
+  for( int i=0; i<n; ++i)
+    std::cout << "x_" << i << " = " << lp.x()[i] << std::endl;
+#endif
+
 
   
   return 0;

NSolver/LPSolveSolver.cc

+//=============================================================================
+//
+//  CLASS LPSolverSolver - IMPLEMENTATION
+//
+//=============================================================================
+
+//== INCLUDES =================================================================
+
+//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
+#include <CoMISo/Config/config.hh>
+#if COMISO_LPSOLVE_AVAILABLE
+
+//=============================================================================
+#include "LPSolveSolver.hh"
+#include <CoMISo/Utils/CoMISoError.hh>
+
+#include <Base/Debug/DebTime.hh>
+#include <Base/Code/Quality.hh>
+
+#include <lp_lib.h>
+
+#include <stdexcept>
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO {
+
+//== IMPLEMENTATION ==========================================================
+
+namespace {
+
+int getLPSolverRowType(NConstraintInterface* c)
+{
+  switch (c->constraint_type())
+  {
+  case NConstraintInterface::NC_EQUAL: return EQ;
+  case NConstraintInterface::NC_LESS_EQUAL: return LE;
+  case NConstraintInterface::NC_GREATER_EQUAL: return GE;
+  }
+}
+
+std::vector<double> get_linear_energy_coefficients(NProblemInterface* _problem)
+{
+  std::vector<double> zero(_problem->n_unknowns(), 0);
+  std::vector<double> q;
+  q.resize(_problem->n_unknowns());
+  _problem->eval_gradient(zero.data(), q.data());
+  return q;
+}
+
+bool solve_impl(
+    NProblemInterface*                        _problem,
+    const std::vector<NConstraintInterface*>& _constraints,
+    const std::vector<PairIndexVtype>&        _discrete_constraints,
+    bool                                      _unbounded_variables,
+    const double                              _time_limit)
+{
+  lprec* lp = nullptr;
+  int Ncol = _problem->n_unknowns();
+  int ret = 0;
+
+  /* We will build the model row by row */
+  lp = make_lp(0, Ncol);
+  if(lp == nullptr)
+    ret = 1; /* couldn't construct a new model... */
+
+  std::vector<int> col_idxs;
+  col_idxs.reserve(Ncol);
+  std::vector<REAL> row_coefficients;
+  row_coefficients.reserve(Ncol);
+
+  auto reset_vectors = [&col_idxs, &row_coefficients]()
+  {
+    col_idxs.clear();
+    row_coefficients.clear();
+    // LPSolve ignores first entry
+    col_idxs.push_back(0);
+    row_coefficients.push_back(0);
+  };
+
+  if (_unbounded_variables)
+    for (size_t i = 1; i < Ncol+1; ++i)
+      set_unbounded(lp, i);
+
+
+  // setup objective function
+  if(ret == 0)
+  {
+    auto coeffs = get_linear_energy_coefficients(_problem);
+    reset_vectors();
+    for (size_t i = 0; i < coeffs.size(); ++i)
+    {
+      col_idxs.push_back(i+1);
+      row_coefficients.push_back(coeffs[i]);
+    }
+
+    /* set the objective in lpsolve */
+    if(!set_obj_fnex(lp, row_coefficients.size(), row_coefficients.data(), col_idxs.data()))
+      ret = 4;
+  }
+
+
+  // add constraints
+  {
+    set_add_rowmode(lp, TRUE);  /* makes building the model faster if it is done rows by row */
+    std::vector<REAL> all_row_coefficients; // including zeros
+    all_row_coefficients.reserve(Ncol);
+    NConstraintInterface::SVectorNC gc;
+
+    for (auto c : _constraints)
+    {
+      reset_vectors();
+
+      if (!c->is_linear())
+      {
+        DEB_error("LPSolve: non-linear constraints are not supported and thus ignored.");
+        continue;
+      }
+
+      c->eval_gradient(all_row_coefficients.data(), gc);
+      for (NConstraintInterface::SVectorNC::InnerIterator v_it(gc); v_it; ++v_it)
+      {
+        col_idxs.push_back(v_it.index()+1);
+        row_coefficients.push_back(v_it.value());
+      }
+
+      const auto b = c->eval_constraint(all_row_coefficients.data());
+      if(!add_constraintex(lp, row_coefficients.size(), row_coefficients.data(), col_idxs.data(), getLPSolverRowType(c), -b))
+        ret = 3;
+    }
+    set_add_rowmode(lp, FALSE); /* rowmode should be turned off again when done building the model */
+  }
+
+  // setup variable types
+  if (ret == 0)
+  {
+    for (auto c : _discrete_constraints)
+    {
+      switch (c.second)
+      {
+      case Real:
+        // real is default
+        break;
+      case Integer:
+        set_int(lp, c.first+1, true);
+        break;
+      case Binary:
+        set_binary(lp, c.first+1, true);
+        break;
+      }
+    }
+  }
+
+  if(ret == 0)
+  {
+    /* set the object direction to maximize */
+    set_minim(lp);
+
+    /* just out of curioucity, now show the model in lp format on screen */
+    /* this only works if this is a console application. If not, use write_lp and a filename */
+    //write_LP(lp, stdout);
+    /* write_lp(lp, "model.lp"); */
+
+    /* I only want to see important messages on screen while solving */
+    set_verbose(lp, IMPORTANT);
+
+    /* Now let lpsolve calculate a solution */
+    ret = solve(lp);
+    if(ret == OPTIMAL)
+      ret = 0;
+    else
+      ret = 5;
+  }
+
+  if(ret == 0)
+  {
+    /* variable values */
+    get_variables(lp, row_coefficients.data());
+    _problem->store_result(row_coefficients.data());
+    /* we are done now */
+  }
+
+  if(lp != nullptr)
+  {
+    /* clean up such that all used memory by lpsolve is freed */
+    delete_lp(lp);
+  }
+
+  return ret == 0;
+}
+
+}//namespace 
+
+bool LPSolveSolver::solve(
+    NProblemInterface*                        _problem,
+    const std::vector<NConstraintInterface*>& _constraints,
+    const std::vector<PairIndexVtype>&        _discrete_constraints,
+    const double                              _time_limit )
+{
+  DEB_enter_func;
+  bool valid_solution = false;
+  try
+  {
+    valid_solution = solve_impl(_problem, _constraints, _discrete_constraints, unbounded_variables_, _time_limit);
+  }
+  catch (...)
+  {
+    DEB_warning(1, "Caught an error");
+    //COMISO_THROW(UNSPECIFIED_CBC_EXCEPTION);
+  }
+  return valid_solution;
+}
+
+
+//=============================================================================
+} // namespace COMISO
+//=============================================================================
+
+#endif // COMISO_CBC_AVAILABLE
+//=============================================================================
+

NSolver/LPSolveSolver.hh

+// (C) Copyright 2015 by Autodesk, Inc.
+
+//=============================================================================
+//
+//  CLASSLPSolverSolver
+//
+//=============================================================================
+#ifndef COMISO_LPSOLVESOLVER_HH
+#define COMISO_LPSOLVESOLVER_HH
+
+//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
+#include <CoMISo/Config/config.hh>
+
+#if COMISO_LPSOLVE_AVAILABLE
+
+//== INCLUDES =================================================================
+
+#include <CoMISo/Config/CoMISoDefines.hh>
+#include <vector>
+#include <string>
+#include <CoMISo/NSolver/NProblemInterface.hh>
+#include <CoMISo/NSolver/NConstraintInterface.hh>
+#include <CoMISo/NSolver/VariableType.hh>
+
+//== FORWARDDECLARATIONS ======================================================
+
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO {
+
+//== CLASS DEFINITION =========================================================
+
+/**
+    Solver interface for LPSolve.
+*/
+class COMISODLLEXPORT LPSolveSolver
+{
+public:
+  // By default variables have a lower bound of 0.
+  // By choosing unbounded variables LPSolve will double the number of variables
+  // which leads to a more expensive solve.
+  LPSolveSolver(bool _unbounded_variables = false)
+    :
+      unbounded_variables_(_unbounded_variables)
+  {}
+
+  // ********** SOLVE **************** //
+  //! \throws Outcome
+  bool solve(
+    NProblemInterface* _problem, // problem instance
+    const std::vector<NConstraintInterface*>& _constraints, // linear constraints
+    const std::vector<PairIndexVtype>& _discrete_constraints, // discrete constraints
+    const double _time_limit = 60); // time limit in seconds
+
+  //! \throws Outcome
+ bool solve(
+    NProblemInterface* _problem, // problem instance
+    const std::vector<NConstraintInterface*>& _constraints, // linear constraints
+    const double _time_limit = 60) // time limit in seconds
+  {
+    std::vector<PairIndexVtype> dc;
+    return solve(_problem, _constraints, dc, _time_limit);
+  }
+
+private:
+
+ bool unbounded_variables_;
+};
+
+
+//=============================================================================
+} // namespace COMISO
+
+//=============================================================================
+#endif // COMISO_CBC_AVAILABLE
+//=============================================================================
+#endif // COMISO_CBCSolver_HH
+//=============================================================================
+

cmake-library @ 3142dc83

-Subproject commit 7b7646e66c33aaf800fc235c6b18a9fec5706ce3
+Subproject commit 3142dc833a8b5a27c5949e787faa55b2450e63ac