Commit 2fe2db99 authored by Martin Marinov's avatar Martin Marinov
Browse files

Merging latest CoMISo updates from RWTH.

[git-p4: depot-paths = "//ReForm/ReForm/main/CoMISo/": change = 12969]
parent 69a442bd
// Build time dependencies for CoMiso
#define COMISO_QT4_AVAILABLE 0
#define COMISO_QT_AVAILABLE 0
#define COMISO_BOOST_AVAILABLE 0
#define COMISO_BLAS_AVAILABLE 1
#define COMISO_GMM_AVAILABLE 1
......@@ -16,11 +16,13 @@
#define COMISO_METIS_AVAILABLE 1
#define COMISO_CGAL_AVAILABLE 1
#define COMISO_TAUCS_AVAILABLE 0
#define COMISO_CBC_AVAILABLE 0
#define COMISO_GUROBI_AVAILABLE 0
#define COMISO_CPLEX_AVAILABLE 0
#define COMISO_DOCLOUD_AVAILABLE 1
#define COMISO_ARPACK_AVAILABLE 0
#define COMISO_CPLEX_AVAILABLE 0
#define COMISO_EIGEN3_AVAILABLE 1
#define COMISO_DCO_AVAILABLE 0
#define COMISO_CBC_AVAILABLE 0
#define COMISO_CLP_AVAILABLE 0
#define COMISO_CGL_AVAILABLE 0
#define COMISO_COINUTILS_AVAILABLE 0
#define COMISO_DOCLOUD_AVAILABLE 1
\ No newline at end of file
......@@ -41,6 +41,7 @@
// Heuristics
#include "CbcHeuristic.hpp"
#include "CbcCompareDepth.hpp"
#include <stdexcept>
#include <stdio.h>
......@@ -257,12 +258,16 @@ void solve_impl(
// Pass the OsiSolver with the problem to be solved to CbcModel
CbcModel model(si);
model.solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry);
model.setMaximumSolutions(4);
TRACE_CBT("CbcModel::getMaximumSolutions()", model.getMaximumSolutions());
//model.setMaximumSeconds(_time_limit);
TRACE_CBT("CbcModel::getMaximumSeconds()", model.getMaximumSeconds());
// model.setMaximumSolutions(4);
// TRACE_CBT("CbcModel::getMaximumSolutions()", model.getMaximumSolutions());
model.setMaximumSeconds(_time_limit);
// TRACE_CBT("CbcModel::getMaximumSeconds()", model.getMaximumSeconds());
model_tweak(model);
CbcCompareDepth compare_depth;
model.setNodeComparison(&compare_depth);
model.branchAndBound();
auto solution = model.bestSolution();
......
......@@ -18,7 +18,6 @@
#include "NProblemInterface.hh"
#include "NConstraintInterface.hh"
#include "VariableType.hh"
#include "GurobiHelper.hh"
//== FORWARDDECLARATIONS ======================================================
class GRBModel;
......@@ -61,6 +60,7 @@ public:
std::vector<PairIndexVtype> dc;
return solve(_problem, _constraints, dc, _time_limit);
}
};
......
......@@ -9,6 +9,7 @@ SET(my_headers
${CMAKE_CURRENT_SOURCE_DIR}/DOCloudConfig.hh
${CMAKE_CURRENT_SOURCE_DIR}/DOCloudJob.hh
${CMAKE_CURRENT_SOURCE_DIR}/DOCloudSolver.hh
${CMAKE_CURRENT_SOURCE_DIR}/FiniteElementProblem.hh
${CMAKE_CURRENT_SOURCE_DIR}/GurobiHelper.hh
${CMAKE_CURRENT_SOURCE_DIR}/GUROBISolver.hh
${CMAKE_CURRENT_SOURCE_DIR}/IPOPTSolver.hh
......@@ -24,6 +25,7 @@ SET(my_headers
${CMAKE_CURRENT_SOURCE_DIR}/NProblemGmmInterface.hh
${CMAKE_CURRENT_SOURCE_DIR}/NProblemInterface.hh
${CMAKE_CURRENT_SOURCE_DIR}/NProblemInterfaceADOLC.hh
${CMAKE_CURRENT_SOURCE_DIR}/NProblemInterfaceDCO.hh
${CMAKE_CURRENT_SOURCE_DIR}/NPTiming.hh
${CMAKE_CURRENT_SOURCE_DIR}/QuadraticProblem.hh
${CMAKE_CURRENT_SOURCE_DIR}/SuperSparseMatrixT.hh
......@@ -52,6 +54,7 @@ SET(my_sources
${CMAKE_CURRENT_SOURCE_DIR}/DOCloudCache.cc
${CMAKE_CURRENT_SOURCE_DIR}/DOCloudJob.cc
${CMAKE_CURRENT_SOURCE_DIR}/DOCloudSolver.cc
${CMAKE_CURRENT_SOURCE_DIR}/FiniteElementProblem.cc
${CMAKE_CURRENT_SOURCE_DIR}/GurobiHelper.cc
${CMAKE_CURRENT_SOURCE_DIR}/GUROBISolver.cc
${CMAKE_CURRENT_SOURCE_DIR}/IPOPTSolver.cc
......@@ -59,6 +62,7 @@ SET(my_sources
${CMAKE_CURRENT_SOURCE_DIR}/LinearConstraint.cc
${CMAKE_CURRENT_SOURCE_DIR}/LinearConstraintHandlerElimination.cc
${CMAKE_CURRENT_SOURCE_DIR}/LinearConstraintHandlerPenalty.cc
${CMAKE_CURRENT_SOURCE_DIR}/LinearProblem.cc
${CMAKE_CURRENT_SOURCE_DIR}/NewtonSolver.cc
${CMAKE_CURRENT_SOURCE_DIR}/NPDerivativeChecker.cc
${CMAKE_CURRENT_SOURCE_DIR}/NPLinearConstraints.cc
......
......@@ -4,6 +4,8 @@
//
//=============================================================================
#define COMISO_CPLEXSOLVER_C
//== INCLUDES =================================================================
//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
......@@ -15,13 +17,20 @@
#if COMISO_CPLEX_AVAILABLE
//=============================================================================
#include <stdexcept>
//== NAMESPACES ===============================================================
namespace COMISO {
//== IMPLEMENTATION ==========================================================
// ********** SOLVE **************** //
bool
CPLEXSolver::
......@@ -267,13 +276,9 @@ CPLEXSolver::
solve(NProblemInterface* _problem,
std::vector<NConstraintInterface*>& _constraints,
std::vector<PairIndexVtype>& _discrete_constraints,
const double /*_time_limit*/,
const double _time_limit,
const bool _silent)
{
//std::cout << "_time_limit: " << _time_limit << std::endl;
MY_PAUSE;
try
{
//----------------------------------------------
......@@ -464,9 +469,8 @@ solve(NProblemInterface* _problem,
if(!_silent)
std::cerr << "cplex -> generate model...\n";
IloCplex cplex(model);
cplex.setParam(IloCplex::IntSolLim, 2);
//cplex.setParam(IloCplex::TiLim, _time_limit);
cplex.setParam(IloCplex::TiLim, _time_limit);
{ // hack
// 0 [CPX_NODESEL_DFS] Depth-first search
// 1 [CPX_NODESEL_BESTBOUND] Best-bound search
......@@ -481,10 +485,6 @@ solve(NProblemInterface* _problem,
if(_silent)
cplex.setOut(env_.getNullStream());
const char model_filename[] = "C:\\temp\\cplex_model.sav";
cout << "Saving the CPLEX model to " << model_filename << std::endl;
cplex.exportModel(model_filename);
IloBool solution_found = cplex.solve();
......@@ -566,14 +566,300 @@ solve(NProblemInterface* _problem,
}
//-----------------------------------------------------------------------------
bool
CPLEXSolver::
solve_two_phase(NProblemInterface* _problem,
std::vector<NConstraintInterface*>& _constraints,
std::vector<PairIndexVtype>& _discrete_constraints,
const double _time_limit0,
const double _gap0,
const double _time_limit1,
const double _gap1,
const bool _silent)
{
try
{
//----------------------------------------------
// 0. set up environment
//----------------------------------------------
if(!_silent)
std::cerr << "cplex -> get environment...\n";
IloEnv env_;
if(!_silent)
std::cerr << "cplex -> get model...\n";
IloModel model(env_);
// model.getEnv().set(GRB_DoubleParam_TimeLimit, _time_limit);
//----------------------------------------------
// 1. allocate variables and initialize limits
//----------------------------------------------
if(!_silent)
std::cerr << "cplex -> allocate variables...\n";
// determine variable types: 0->real, 1->integer, 2->bool
std::vector<char> vtypes (_problem->n_unknowns(),0);
for(unsigned int i=0; i<_discrete_constraints.size(); ++i)
if(_discrete_constraints[i].first < vtypes.size())
{
switch(_discrete_constraints[i].second)
{
case Integer: vtypes [_discrete_constraints[i].first] = 1;
break;
case Binary : vtypes[_discrete_constraints[i].first] = 2;
break;
default : break;
}
}
else
std::cerr << "ERROR: requested a discrete variable which is above the total number of variables"
<< _discrete_constraints[i].first << " vs " << vtypes.size() << std::endl;
// CPLEX variables
std::vector<IloNumVar> vars; vars.reserve(_problem->n_unknowns());
// first all
for( int i=0; i<_problem->n_unknowns(); ++i)
switch(vtypes[i])
{
case 0 : vars.push_back( IloNumVar(env_,-IloInfinity, IloInfinity, IloNumVar::Float) ); break;
case 1 : vars.push_back( IloNumVar(env_, -IloIntMax, IloIntMax, IloNumVar::Int) ); break;
case 2 : vars.push_back( IloNumVar(env_, 0, 1, IloNumVar::Bool) ); break;
default: std::cerr << "ERROR: unhandled varibalbe type in CPLEXSolver!!!" << std::endl;
}
if(!_silent)
{
std::cerr << "#unknowns : " << _problem->n_unknowns() << std::endl;
std::cerr << "#CPLEX variables: " << vars.size() << std::endl;
}
// Integrate new variables
// model.update();
//----------------------------------------------
// 2. setup constraints
//----------------------------------------------
if(!_silent)
std::cerr << "cplex -> setup constraints...\n";
// get zero vector
std::vector<double> x(_problem->n_unknowns(), 0.0);
// handle constraints depending on their tyep
for(unsigned int i=0; i<_constraints.size(); ++i)
{
// is linear constraint?
LinearConstraint* lin_ptr = dynamic_cast<LinearConstraint*>(_constraints[i]);
if(lin_ptr)
{
IloExpr lin_expr(env_);
NConstraintInterface::SVectorNC gc;
_constraints[i]->eval_gradient(P(x), gc);
NConstraintInterface::SVectorNC::InnerIterator v_it(gc);
for(; v_it; ++v_it)
{
assert(v_it.index() >= 0 && v_it.index() < _problem->n_unknowns());
lin_expr += vars[v_it.index()]*v_it.value();
}
double b = _constraints[i]->eval_constraint(P(x));
switch(_constraints[i]->constraint_type())
{
case NConstraintInterface::NC_EQUAL : model.add(lin_expr + b == 0); break;
case NConstraintInterface::NC_LESS_EQUAL : model.add(lin_expr + b <= 0); break;
case NConstraintInterface::NC_GREATER_EQUAL : model.add(lin_expr + b >= 0); break;
}
}
else
{
BoundConstraint* bnd_ptr = dynamic_cast<BoundConstraint*>(_constraints[i]);
if(bnd_ptr)
{
assert( int(bnd_ptr->idx()) < _problem->n_unknowns());
switch(bnd_ptr->constraint_type())
{
case NConstraintInterface::NC_EQUAL : vars[bnd_ptr->idx()].setBounds(bnd_ptr->bound(), bnd_ptr->bound()); break;
case NConstraintInterface::NC_LESS_EQUAL : vars[bnd_ptr->idx()].setUB(bnd_ptr->bound()); break;
case NConstraintInterface::NC_GREATER_EQUAL : vars[bnd_ptr->idx()].setLB(bnd_ptr->bound()); break;
}
}
else
{
ConeConstraint* cone_ptr = dynamic_cast<ConeConstraint*>(_constraints[i]);
if(cone_ptr)
{
IloExpr soc_lhs(env_);
IloExpr soc_rhs(env_);
assert(cone_ptr->i() >= 0 && cone_ptr->i() <= _problem->n_unknowns());
soc_rhs= 0.5*cone_ptr->c()*vars[cone_ptr->i()]*vars[cone_ptr->i()];
NConstraintInterface::SMatrixNC::iterator q_it = cone_ptr->Q().begin();
for(; q_it != cone_ptr->Q().end(); ++q_it)
{
assert( int(q_it.col()) < _problem->n_unknowns());
assert( int(q_it.row()) < _problem->n_unknowns());
soc_lhs += 0.5*(*q_it)*vars[q_it.col()]*vars[q_it.row()];
}
model.add(soc_lhs <= soc_rhs);
}
else
std::cerr << "Warning: CPLEXSolver received a constraint of unknow type!!! -> skipping it" << std::endl;
}
}
}
// model.update();
//----------------------------------------------
// 3. setup energy
//----------------------------------------------
if(!_silent)
std::cerr << "cplex -> setup energy...\n";
if(!_problem->constant_hessian())
std::cerr << "Warning: CPLEXSolver received a problem with non-constant hessian!!!" << std::endl;
// GRBQuadExpr objective;
IloExpr objective(env_);
// add quadratic part
NProblemInterface::SMatrixNP H;
_problem->eval_hessian(P(x), H);
for( int i=0; i<H.outerSize(); ++i)
for (NProblemInterface::SMatrixNP::InnerIterator it(H,i); it; ++it)
{
assert(it.row() >= 0 && it.row() < _problem->n_unknowns());
assert(it.col() >= 0 && it.col() < _problem->n_unknowns());
objective += 0.5*it.value()*vars[it.row()]*vars[it.col()];
}
// add linear part
std::vector<double> g(_problem->n_unknowns());
_problem->eval_gradient(P(x), P(g));
for(unsigned int i=0; i<g.size(); ++i)
{
objective += g[i]*vars[i];
}
// add constant part
objective += _problem->eval_f(P(x));
model.add(IloMinimize(env_,objective));
// model.set(GRB_IntAttr_ModelSense, 1);
// model.setObjective(objective);
// model.update();
//----------------------------------------------
// 4. solve problem
//----------------------------------------------
if(!_silent)
std::cerr << "cplex -> generate model...\n";
IloCplex cplex(model);
cplex.setParam(IloCplex::TiLim, _time_limit0);
cplex.setParam(IloCplex::EpGap, _gap0);
{ // hack
// 0 [CPX_NODESEL_DFS] Depth-first search
// 1 [CPX_NODESEL_BESTBOUND] Best-bound search
// 2 [CPX_NODESEL_BESTEST] Best-estimate search
// 3 [CPX_NODESEL_BESTEST_ALT] Alternative best-estimate search
// cplex.setParam(IloCplex::NodeSel , 0);
}
if(!_silent)
std::cerr << "cplex -> solve...\n";
// silent mode?
if(_silent)
cplex.setOut(env_.getNullStream());
IloBool solution_found = cplex.solve();
// phase 2?
if(_time_limit1 > 0 && (!solution_found || cplex.getMIPRelativeGap() > _gap1))
{
cplex.setParam(IloCplex::TiLim, _time_limit1);
cplex.setParam(IloCplex::EpGap, _gap1);
solution_found = cplex.solve();
}
//----------------------------------------------
// 5. store result
//----------------------------------------------
if(solution_found != IloFalse)
{
if(!_silent)
std::cerr << "cplex -> store result...\n";
// store computed result
for(unsigned int i=0; i<vars.size(); ++i)
x[i] = cplex.getValue(vars[i]);
_problem->store_result(P(x));
}
/*
if (solution_input_path_.empty())
{
// store computed result
for(unsigned int i=0; i<vars.size(); ++i)
x[i] = vars[i].get(GRB_DoubleAttr_X);
}
*/
// else
// {
// std::cout << "Loading stored solution from \"" << solution_input_path_ << "\"." << std::endl;
// // store loaded result
// const size_t oldSize = x.size();
// x.clear();
// GurobiHelper::readSolutionVectorFromSOL(x, solution_input_path_);
// if (oldSize != x.size()) {
// std::cerr << "oldSize != x.size() <=> " << oldSize << " != " << x.size() << std::endl;
// throw std::runtime_error("Loaded solution vector doesn't have expected dimension.");
// }
// }
//
// _problem->store_result(P(x));
//
// // ObjVal is only available if the optimize was called.
// if (solution_input_path_.empty())
// std::cout << "GUROBI Objective: " << model.get(GRB_DoubleAttr_ObjVal) << std::endl;
return solution_found;
}
catch (IloException& e)
{
cerr << "Concert exception caught: " << e << endl;
return false;
}
catch (...)
{
cerr << "Unknown exception caught" << endl;
return false;
}
return false;
}
//-----------------------------------------------------------------------------
void
CPLEXSolver::
test()
{
MY_PAUSE;
try
{
std::cerr << "cplex2 -> get environment...\n";
......@@ -638,13 +924,12 @@ solve(NProblemInterface* _problem,
const double _time_limit,
const bool _silent)
{
// std::cerr << "Warning: CPLEXSolver does not support lazy constraints yet -> solve with all constraints instead" << std::endl;
// std::vector<NConstraintInterface*> C;
// std::copy(_constraints.begin(),_constraints.end(),std::back_inserter(C));
// return solve(_problem, C, _time_limit, _silent);
DEB_time_func_def;
StopWatch sw; sw.start();
bool feasible_point_found = false;
int cur_pass = 0;
......@@ -900,6 +1185,8 @@ solve(NProblemInterface* _problem,
}
}
const double overall_time = sw.stop()/1000.0;
//----------------------------------------------------------------------------
// 4. output statistics
//----------------------------------------------------------------------------
......@@ -914,6 +1201,7 @@ solve(NProblemInterface* _problem,
// }
std::cerr <<"############# CPLEX with lazy constraints statistics ###############" << std::endl;
std::cerr << "overall time: " << overall_time << "s" << std::endl;
std::cerr << "#passes : " << cur_pass << "( of " << _max_passes << ")" << std::endl;
for(unsigned int i=0; i<n_inf.size(); ++i)
std::cerr << "pass " << i << " induced " << n_inf[i] << " infeasible and " << n_almost_inf[i] << " almost infeasible" << std::endl;
......@@ -943,7 +1231,6 @@ void
CPLEXSolver::
add_constraint_to_model( NConstraintInterface* _constraint, std::vector<IloNumVar>& _vars, IloModel& _model, IloEnv& _env, double* _x)
{
MY_PAUSE;
// is linear constraint?
LinearConstraint* lin_ptr = dynamic_cast<LinearConstraint*>(_constraint);
......@@ -956,7 +1243,7 @@ add_constraint_to_model( NConstraintInterface* _constraint, std::vector<IloNumVa
NConstraintInterface::SVectorNC::InnerIterator v_it(gc);
for(; v_it; ++v_it)
{
assert(v_it.index() >= 0 && v_it.index() < _vars.size());
assert(v_it.index() >= 0 && v_it.index() < (int)_vars.size());
lin_expr += _vars[v_it.index()]*v_it.value();
}
......@@ -975,7 +1262,7 @@ add_constraint_to_model( NConstraintInterface* _constraint, std::vector<IloNumVa
BoundConstraint* bnd_ptr = dynamic_cast<BoundConstraint*>(_constraint);
if(bnd_ptr)
{
assert( int(bnd_ptr->idx()) < _vars.size());
assert( int(bnd_ptr->idx()) < (int)_vars.size());
switch(bnd_ptr->constraint_type())
{
......@@ -992,15 +1279,15 @@ add_constraint_to_model( NConstraintInterface* _constraint, std::vector<IloNumVa
IloExpr soc_lhs(_env);
IloExpr soc_rhs(_env);
assert(cone_ptr->i() >= 0 && cone_ptr->i() <= _vars.size());
assert(cone_ptr->i() >= 0 && cone_ptr->i() <= (int)_vars.size());
soc_rhs= 0.5*cone_ptr->c()*_vars[cone_ptr->i()]*_vars[cone_ptr->i()];
NConstraintInterface::SMatrixNC::iterator q_it = cone_ptr->Q().begin();
for(; q_it != cone_ptr->Q().end(); ++q_it)
{
assert( int(q_it.col()) < _vars.size());
assert( int(q_it.row()) < _vars.size());
assert( int(q_it.col()) < (int)_vars.size());
assert( int(q_it.row()) < (int)_vars.size());
soc_lhs += 0.5*(*q_it)*_vars[q_it.col()]*_vars[q_it.row()];
}
......@@ -1013,41 +1300,6 @@ add_constraint_to_model( NConstraintInterface* _constraint, std::vector<IloNumVa
}
}
//-----------------------------------------------------------------------------
//void
//CPLEXSolver::
//set_problem_output_path( const std::string &_problem_output_path)
//{
// problem_output_path_ = _problem_output_path;
//}
//
//
////-----------------------------------------------------------------------------
//
//
//void
//CPLEXSolver::
//set_problem_env_output_path( const std::string &_problem_env_output_path)
//{
// problem_env_output_path_ = _problem_env_output_path;
//}
//
//
////-----------------------------------------------------------------------------
//
//
//void
//CPLEXSolver::
//set_solution_input_path(const std::string &_solution_input_path)
//{
// solution_input_path_ = _solution_input_path;
//}
//== IMPLEMENTATION ==========================================================
......@@ -1094,4 +1346,4 @@ CPLEXSolver()
} // namespace COMISO
//=============================================================================
#endif // COMISO_CPLEX_AVAILABLE
//=============================================================================
//=============================================================================
\ No newline at end of file
......@@ -66,6 +66,20 @@ public:
const bool _silent = false) // time limit in seconds
{ std::vector<PairIndexVtype> dc; return solve(_problem, _constraints, dc, _time_limit, _silent);}
// 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);
// optimization with additional lazy constraints that are only added iteratively to the problem if not satisfied
bool solve(NProblemInterface* _problem,
......
#include "FiniteElementProblem.hh"
namespace COMISO {
// FiniteElementProblem