Commit 078354d6 authored by Martin Marinov's avatar Martin Marinov
Browse files

Merged official CoMISo compatibility changes from the research branch.

[git-p4: depot-paths = "//ReForm/ReForm/main/CoMISo/": change = 13619]
parent 4b3f7b30
......@@ -11,7 +11,8 @@
#define COMISO_HSL_AVAILABLE 0
#define COMISO_PETSC_AVAILABLE 0
#define COMISO_TAO_AVAILABLE 0
#define COMISO_IPOPT_AVAILABLE 1
#define COMISO_IPOPT_AVAILABLE 0
#define COMISO_IPOPTLEAN_AVAILABLE 1
#define COMISO_MUMPS_AVAILABLE 1
#define COMISO_METIS_AVAILABLE 1
#define COMISO_CGAL_AVAILABLE 1
......
......@@ -44,16 +44,9 @@
#include "CbcCompareDepth.hpp"
#include <stdexcept>
#include <stdio.h>
DEB_module("CBCSolver")
#define CBC_INFINITY COIN_DBL_MAX
#ifndef _MSC_VER
#define sprintf_s snprintf
#endif
//== NAMESPACES ===============================================================
namespace COMISO {
......@@ -142,7 +135,7 @@ void model_tweak(CbcModel& model)
#define TRACE_CBT(DESCR, EXPR) DEB_line(7, DESCR << ": " << EXPR)
#define P(X) ((X).data())
void solve_impl(
bool solve_impl(
NProblemInterface* _problem,
const std::vector<NConstraintInterface*>& _constraints,
const std::vector<PairIndexVtype>& _discrete_constraints,
......@@ -253,12 +246,8 @@ void solve_impl(
// TODO: make this accessible through the DEB system instead
volatile static bool dump_problem = false; // change on run-time if necessary
if (dump_problem)
{
static int n_mps_dumps = 0;
char filename[64];
sprintf_s(filename, sizeof(filename), "CBC_problem_dump_%04i", n_mps_dumps++);
si.writeMps(filename); //output problem as .MPS
}
si.writeMps("CBC_problem_dump"); //output problem as .MPS
// Pass the OsiSolver with the problem to be solved to CbcModel
CbcModel model(si);
model.solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry);
......@@ -274,17 +263,24 @@ void solve_impl(
model.branchAndBound();
auto solution = model.bestSolution();
COMISO_THROW_if(solution == nullptr, UNSPECIFIED_CBC_EXCEPTION);
_problem->store_result(solution);
const double* solution = model.bestSolution();
//COMISO_THROW_if(solution == nullptr, UNSPECIFIED_CBC_EXCEPTION);
// store if solution available
if(solution != 0)
{
_problem->store_result(solution);
return true;
}
else
return false;
}
#undef TRACE_CBT
#undef P
}//namespace
void CBCSolver::solve(
bool CBCSolver::solve(
NProblemInterface* _problem,
const std::vector<NConstraintInterface*>& _constraints,
const std::vector<PairIndexVtype>& _discrete_constraints,
......@@ -292,15 +288,17 @@ void CBCSolver::solve(
)
{
DEB_enter_func;
bool valid_solution = false;
try
{
solve_impl(_problem, _constraints, _discrete_constraints, _time_limit);
valid_solution = solve_impl(_problem, _constraints, _discrete_constraints, _time_limit);
}
catch (CoinError& ce)
{
DEB_warning(1, "CoinError code = " << ce.message() << "]\n");
COMISO_THROW(UNSPECIFIED_CBC_EXCEPTION);
//COMISO_THROW(UNSPECIFIED_CBC_EXCEPTION);
}
return valid_solution;
}
......
//=============================================================================
//
// CLASS CPCSolver
// CLASS CBCSolver
//
//=============================================================================
#ifndef COMISO_CBCSolver_HH
......@@ -45,14 +45,14 @@ public:
// ********** SOLVE **************** //
//! \throws Outcome
void solve(
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
void solve(
bool solve(
NProblemInterface* _problem, // problem instance
const std::vector<NConstraintInterface*>& _constraints, // linear constraints
const double _time_limit = 60) // time limit in seconds
......@@ -61,6 +61,19 @@ public:
return solve(_problem, _constraints, dc, _time_limit);
}
// // same as previous but more control over stopping criteria
// // the optimizer runs in two phases
// // phase 1) stops if solution with a MIP gap lower than _gap0 is found or _time_limit0 is reached
// // phase 2) starts only if in phase 1 no solution with a MIP gap lower than _gap1 was found and
// // uses _gap1 and _time_limit2 as parameters
// inline bool solve_two_phase(NProblemInterface* _problem, // problem instance
// std::vector<NConstraintInterface*>& _constraints, // linear constraints
// std::vector<PairIndexVtype>& _discrete_constraints, // discrete constraints
// const double _time_limit0 = 60, // time limit phase 1 in seconds
// const double _gap0 = 0.001, // MIP gap phase 1
// const double _time_limit1 = 120, // time limit phase 2 in seconds
// const double _gap1 = 0.2, // MIP gap phase 2
// const bool _silent = false);
};
......
......@@ -12,7 +12,7 @@ SET(my_headers
${CMAKE_CURRENT_SOURCE_DIR}/FiniteElementProblem.hh
${CMAKE_CURRENT_SOURCE_DIR}/GurobiHelper.hh
${CMAKE_CURRENT_SOURCE_DIR}/GUROBISolver.hh
${CMAKE_CURRENT_SOURCE_DIR}/IPOPTSolver.hh
${CMAKE_CURRENT_SOURCE_DIR}/IPOPTSolverLean.hh
${CMAKE_CURRENT_SOURCE_DIR}/LeastSquaresProblem.hh
${CMAKE_CURRENT_SOURCE_DIR}/LinearConstraint.hh
${CMAKE_CURRENT_SOURCE_DIR}/LinearConstraintHandlerElimination.hh
......@@ -57,7 +57,7 @@ SET(my_sources
${CMAKE_CURRENT_SOURCE_DIR}/FiniteElementProblem.cc
${CMAKE_CURRENT_SOURCE_DIR}/GurobiHelper.cc
${CMAKE_CURRENT_SOURCE_DIR}/GUROBISolver.cc
${CMAKE_CURRENT_SOURCE_DIR}/IPOPTSolver.cc
${CMAKE_CURRENT_SOURCE_DIR}/IPOPTSolverLean.cc
${CMAKE_CURRENT_SOURCE_DIR}/LeastSquaresProblem.cc
${CMAKE_CURRENT_SOURCE_DIR}/LinearConstraint.cc
${CMAKE_CURRENT_SOURCE_DIR}/LinearConstraintHandlerElimination.cc
......
......@@ -22,10 +22,6 @@
#include <Base/Debug/DebTime.hh>
#include <stdexcept>
DEB_module("Gurobi")
//== NAMESPACES ===============================================================
namespace COMISO {
......@@ -90,12 +86,12 @@ static void process_gurobi_exception(const GRBException& _exc)
}
}
void
bool
GUROBISolver::
solve(NProblemInterface* _problem,
std::vector<NConstraintInterface*>& _constraints,
std::vector<PairIndexVtype>& _discrete_constraints,
const double _time_limit)
std::vector<NConstraintInterface*>& _constraints,
std::vector<PairIndexVtype>& _discrete_constraints,
const double _time_limit)
{
DEB_enter_func;
try
......@@ -163,9 +159,9 @@ solve(NProblemInterface* _problem,
switch(_constraints[i]->constraint_type())
{
case NConstraintInterface::NC_EQUAL : model.addConstr(lin_expr + b == 0); break;
case NConstraintInterface::NC_LESS_EQUAL : model.addConstr(lin_expr + b <= 0); break;
case NConstraintInterface::NC_GREATER_EQUAL : model.addConstr(lin_expr + b >= 0); break;
case NConstraintInterface::NC_EQUAL : model.addConstr(lin_expr + b == 0); break;
case NConstraintInterface::NC_LESS_EQUAL : model.addConstr(lin_expr + b <= 0); break;
case NConstraintInterface::NC_GREATER_EQUAL : model.addConstr(lin_expr + b >= 0); break;
}
}
model.update();
......@@ -259,17 +255,28 @@ solve(NProblemInterface* _problem,
// ObjVal is only available if the optimize was called.
DEB_out_if(solution_input_path_.empty(), 2,
"GUROBI Objective: " << model.get(GRB_DoubleAttr_ObjVal) << "\n");
return true;
}
catch (const GRBException& e)
catch(GRBException& e)
{
process_gurobi_exception(e);
//process_gurobi_exception(e);
std::cout << "Error code = " << e.getErrorCode() << std::endl;
std::cout << e.getMessage() << std::endl;
return false;
}
catch(...)
{
std::cout << "Exception during optimization" << std::endl;
return false;
}
return false;
}
//-----------------------------------------------------------------------------
void
bool
GUROBISolver::
solve_two_phase(NProblemInterface* _problem, // problem instance
std::vector<NConstraintInterface*>& _constraints, // linear constraints
......@@ -284,7 +291,7 @@ solve_two_phase(NProblemInterface* _problem, //
}
void
bool
GUROBISolver::
solve_two_phase(NProblemInterface* _problem, // problem instance
std::vector<NConstraintInterface*>& _constraints, // linear constraints
......@@ -295,7 +302,6 @@ solve_two_phase(NProblemInterface* _problem, //
const double _gap1, // MIP gap phase 2
double& _final_gap) //return final gap
{
DEB_enter_func;
try
{
//----------------------------------------------
......@@ -468,18 +474,29 @@ solve_two_phase(NProblemInterface* _problem, //
// ObjVal is only available if the optimize was called.
DEB_line_if(solution_input_path_.empty(), 2,
"GUROBI Objective: " << model.get(GRB_DoubleAttr_ObjVal));
return true;
}
catch (const GRBException& e)
catch(GRBException& e)
{
process_gurobi_exception(e);
//process_gurobi_exception(e);
std::cout << "Error code = " << e.getErrorCode() << std::endl;
std::cout << e.getMessage() << std::endl;
return false;
}
catch(...)
{
std::cout << "Exception during optimization" << std::endl;
return false;
}
return false;
}
//-----------------------------------------------------------------------------
void
bool
GUROBISolver::
solve(NProblemInterface* _problem,
const std::vector<NConstraintInterface*>& _constraints,
......@@ -695,6 +712,7 @@ solve(NProblemInterface* _problem,
x[i] = vars[i].get(GRB_DoubleAttr_X);
}
const double overall_time = sw.stop()/1000.0;
//----------------------------------------------------------------------------
// 4. output statistics
......@@ -713,15 +731,33 @@ solve(NProblemInterface* _problem,
//----------------------------------------------
// 5. store result
//----------------------------------------------
COMISO_THROW_TODO_if(model.get(GRB_IntAttr_Status) != GRB_OPTIMAL,
"Gurobi solution not optimal");
_problem->store_result(P(x));
//COMISO_THROW_TODO_if(model.get(GRB_IntAttr_Status) != GRB_OPTIMAL,
// "Gurobi solution not optimal");
if(model.get(GRB_IntAttr_Status) != GRB_OPTIMAL)
return false;
else
{
_problem->store_result(P(x));
return true;
}
}
catch(GRBException& e)
{
//process_gurobi_exception(e);
std::cout << "Error code = " << e.getErrorCode() << std::endl;
std::cout << e.getMessage() << std::endl;
return false;
}
catch (const GRBException& e)
catch(...)
{
process_gurobi_exception(e);
std::cout << "Exception during optimization" << std::endl;
return false;
}
return false;
// //----------------------------------------------
// // 4. iteratively solve problem
// //----------------------------------------------
......
......@@ -50,16 +50,12 @@ public:
~GUROBISolver() {}
// ********** SOLVE **************** //
//! \throws Outcome
void solve(
NProblemInterface* _problem, // problem instance
bool solve(NProblemInterface* _problem, // problem instance
std::vector<NConstraintInterface*>& _constraints, // linear constraints
std::vector<PairIndexVtype>& _discrete_constraints, // discrete constraints
const double _time_limit = 60); // time limit in seconds
//! \throws Outcome
void solve(
NProblemInterface* _problem, // problem instance
bool solve(NProblemInterface* _problem, // problem instance
std::vector<NConstraintInterface*>& _constraints, // linear constraints
const double _time_limit = 60) // time limit in seconds
{
......@@ -72,8 +68,7 @@ public:
// phase 1) stops if solution with a MIP gap lower than _gap0 is found or _time_limit0 is reached
// phase 2) starts only if in phase 1 no solution with a MIP gap lower than _gap1 was found and
// uses _gap1 and _time_limit2 as parameters
void solve_two_phase(
NProblemInterface* _problem, // problem instance
bool solve_two_phase(NProblemInterface* _problem, // problem instance
std::vector<NConstraintInterface*>& _constraints, // linear constraints
std::vector<PairIndexVtype>& _discrete_constraints, // discrete constraints
const double _time_limit0 = 60, // time limit phase 1 in seconds
......@@ -81,8 +76,7 @@ public:
const double _time_limit1 = 120, // time limit phase 2 in seconds
const double _gap1 = 0.2); // MIP gap phase 2
void solve_two_phase(
NProblemInterface* _problem, // problem instance
bool solve_two_phase(NProblemInterface* _problem, // problem instance
std::vector<NConstraintInterface*>& _constraints, // linear constraints
std::vector<PairIndexVtype>& _discrete_constraints, // discrete constraints
const double _time_limit0, // time limit phase 1 in seconds
......@@ -93,7 +87,7 @@ public:
// optimization with additional lazy constraints that are only added iteratively to the problem if not satisfied
void solve(NProblemInterface* _problem,
bool solve(NProblemInterface* _problem,
const std::vector<NConstraintInterface*>& _constraints,
const std::vector<NConstraintInterface*>& _lazy_constraints,
std::vector<PairIndexVtype>& _discrete_constraints, // discrete constraints
......@@ -104,7 +98,7 @@ public:
// same as above with additional lazy constraints that are only added iteratively to the problem if not satisfied
void solve(NProblemInterface* _problem,
bool solve(NProblemInterface* _problem,
const std::vector<NConstraintInterface*>& _constraints,
const std::vector<NConstraintInterface*>& _lazy_constraints,
const double _almost_infeasible = 0.5,
......
This diff is collapsed.
......@@ -16,21 +16,27 @@
//== INCLUDES =================================================================
#include <CoMISo/Config/CoMISoDefines.hh>
#include <CoMISo/Utils/StopWatch.hh>
#include <vector>
#include <cstddef>
#include <gmm/gmm.h>
#include "NProblemGmmInterface.hh"
#include "NProblemInterface.hh"
#include "NConstraintInterface.hh"
#include "BoundConstraint.hh"
#include <IpTNLP.hpp>
#include <IpIpoptApplication.hpp>
#include <IpSolveStatistics.hpp>
//== FORWARDDECLARATIONS ======================================================
//== NAMESPACES ===============================================================
namespace COMISO {
//== FORWARDDECLARATIONS ======================================================
class NProblemGmmInterface; // deprecated
class NProblemInterface;
class NConstraintInterface;
//== CLASS DEFINITION =========================================================
/** \class NewtonSolver NewtonSolver.hh <ACG/.../NewtonSolver.hh>
......@@ -41,45 +47,366 @@ class NConstraintInterface;
class COMISODLLEXPORT IPOPTSolver
{
public:
/// Default constructor
/// Default constructor -> set up IpOptApplication
IPOPTSolver();
/// Destructor
~IPOPTSolver();
~IPOPTSolver() {}
// ********** SOLVE **************** //
//! \throws Outcome
void solve(NProblemInterface* _problem,
const std::vector<NConstraintInterface*>& _constraints);
//! Same as above with additional lazy constraints that are only added iteratively to the problem if not satisfied
//! \throws Outcome
void solve(NProblemInterface* _problem,
const std::vector<NConstraintInterface*>& _constraints,
const std::vector<NConstraintInterface*>& _lazy_constraints,
const double _almost_infeasible = 0.5,
const int _max_passes = 5 );
// solve -> returns ipopt status code
//------------------------------------------------------
// enum ApplicationReturnStatus
// {
// Solve_Succeeded=0,
// Solved_To_Acceptable_Level=1,
// Infeasible_Problem_Detected=2,
// Search_Direction_Becomes_Too_Small=3,
// Diverging_Iterates=4,
// User_Requested_Stop=5,
// Feasible_Point_Found=6,
//
// Maximum_Iterations_Exceeded=-1,
// Restoration_Failed=-2,
// Error_In_Step_Computation=-3,
// Maximum_CpuTime_Exceeded=-4,
// Not_Enough_Degrees_Of_Freedom=-10,
// Invalid_Problem_Definition=-11,
// Invalid_Option=-12,
// Invalid_Number_Detected=-13,
//
// Unrecoverable_Exception=-100,
// NonIpopt_Exception_Thrown=-101,
// Insufficient_Memory=-102,
// Internal_Error=-199
// };
//------------------------------------------------------
int solve(NProblemInterface* _problem, const std::vector<NConstraintInterface*>& _constraints);
// same as above with additional lazy constraints that are only added iteratively to the problem if not satisfied
int solve(NProblemInterface* _problem,
const std::vector<NConstraintInterface*>& _constraints,
const std::vector<NConstraintInterface*>& _lazy_constraints,
const double _almost_infeasible = 0.5,
const int _max_passes = 5 );
// for convenience, if no constraints are given
//! \throws Outcome
void solve(NProblemInterface* _problem);
int solve(NProblemInterface* _problem);
// deprecated interface for backwards compatibility
//! \throws Outcome
void solve(NProblemGmmInterface* _problem,
std::vector<NConstraintInterface*>& _constraints);
int solve(NProblemGmmInterface* _problem, std::vector<NConstraintInterface*>& _constraints);
// ********* CONFIGURATION ********************* //
// access the ipopt-application (for setting parameters etc.)
// examples: app().Options()->SetIntegerValue("max_iter", 100);
// app().Options()->SetStringValue("derivative_test", "second-order");
Ipopt::IpoptApplication& app() {return (*app_); }
void set_print_level(const int _level)
{
app().Options()->SetIntegerValue("print_level", _level);
print_level_ = _level;
}
protected:
double* P(std::vector<double>& _v)
{
if( !_v.empty())
return ((double*)&_v[0]);
else
return 0;
}
private:
// smart pointer to IpoptApplication to set options etc.
Ipopt::SmartPtr<Ipopt::IpoptApplication> app_;
int print_level_;
};
//== CLASS DEFINITION PROBLEM INSTANCE=========================================================
class NProblemIPOPT : public Ipopt::TNLP
{
public:
// Ipopt Types
typedef Ipopt::Number Number;
typedef Ipopt::Index Index;
typedef Ipopt::SolverReturn SolverReturn;
typedef Ipopt::IpoptData IpoptData;
typedef Ipopt::IpoptCalculatedQuantities IpoptCalculatedQuantities;
// sparse matrix and vector types
typedef NConstraintInterface::SVectorNC SVectorNC;
typedef NConstraintInterface::SMatrixNC SMatrixNC;
typedef NProblemInterface::SMatrixNP SMatrixNP;
/** default constructor */
NProblemIPOPT(NProblemInterface* _problem, const std::vector<NConstraintInterface*>& _constraints, const bool _hessian_approximation = false)
: problem_(_problem), store_solution_(false), hessian_approximation_(_hessian_approximation)
{
// initialize solution
x_.resize(_problem->n_unknowns());
_problem->initial_x(P(x_));
split_constraints(_constraints);
analyze_special_properties(_problem, _constraints);
}
/** default destructor */
virtual ~NProblemIPOPT() {};
/**@name Overloaded from TNLP */
//@{
/** Method to return some info about the nlp */
virtual bool get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
Index& nnz_h_lag, IndexStyleEnum& index_style);
/** Method to return the bounds for my problem */
virtual bool get_bounds_info(Index n, Number* x_l, Number* x_u,
Index m, Number* g_l, Number* g_u);
/** Method to return the starting point for the algorithm */
virtual bool get_starting_point(Index n, bool init_x, Number* x,
bool init_z, Number* z_L, Number* z_U,
Index m, bool init_lambda,
Number* lambda);
/** Method to return the objective value */
virtual bool eval_f(Index n, const Number* x, bool new_x, Number& obj_value);
/** Method to return the gradient of the objective */
virtual bool eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f);
/** Method to return the constraint residuals */
virtual bool eval_g(Index n, const Number* x, bool new_x, Index m, Number* g);
/** Method to return:
* 1) The structure of the jacobian (if "values" is NULL)
* 2) The values of the jacobian (if "values" is not NULL)
*/
virtual bool eval_jac_g(Index n, const Number* x, bool new_x,
Index m, Index nele_jac, Index* iRow, Index *jCol,
Number* values);
/** Method to return:
* 1) The structure of the hessian of the lagrangian (if "values" is NULL)
* 2) The values of the hessian of the lagrangian (if "values" is not NULL)
*/
virtual bool eval_h(Index n, const Number* x, bool new_x,
Number obj_factor, Index m, const Number* lambda,
bool new_lambda, Index nele_hess, Index* iRow,
Index* jCol, Number* values);
//@}
/** @name Solution Methods */
//@{
/** This method is called when the algorithm is complete so the TNLP can store/write the solution */
virtual void finalize_solution(SolverReturn status,
Index n, const Number* x, const Number* z_L, const Number* z_U,
Index m, const Number* g, const Number* lambda,
Number obj_value,
const IpoptData* ip_data,
IpoptCalculatedQuantities* ip_cq);
//@}
// special properties of problem
bool hessian_constant() const;
bool jac_c_constant() const;
bool jac_d_constant() const;
bool& store_solution() {return store_solution_; }
std::vector<double>& solution() {return x_;}
//! Get the computed solution energy