Commit 5bc17d5c authored by Martin Marinov's avatar Martin Marinov
Browse files

Merged the CoMISo changes from ReForm into the official git repo.

parent a9b5509a
Pipeline #302 failed with stage
// (C) Copyright 2014 by Autodesk, Inc.
//
// The information contained herein is confidential, proprietary
// to Autodesk, Inc., and considered a trade secret as defined
// in section 499C of the penal code of the State of California.
// Use of this information by anyone other than authorized
// employees of Autodesk, Inc. is granted only under a written
// non-disclosure agreement, expressly prescribing the scope
// and manner of such use.
#ifndef GMMTYPES_HH_INCLUDED
#define GMMTYPES_HH_INCLUDED
#include <gmm/gmm_matrix.h>
namespace COMISO_GMM
{
// Matrix Types
typedef gmm::col_matrix< gmm::wsvector<double> > WSColMatrix;
typedef gmm::row_matrix< gmm::wsvector<double> > WSRowMatrix;
typedef gmm::col_matrix< gmm::rsvector<double> > RSColMatrix;
typedef gmm::row_matrix< gmm::rsvector<double> > RSRowMatrix;
typedef gmm::csc_matrix<double> CSCMatrix;
}//namespace COMISO_GMM
#endif//GMMTYPES_HH_INCLUDED
// (C) Copyright 2014 by Autodesk, Inc.
//
// The information contained herein is confidential, proprietary
// to Autodesk, Inc., and considered a trade secret as defined
// in section 499C of the penal code of the State of California.
// Use of this information by anyone other than authorized
// employees of Autodesk, Inc. is granted only under a written
// non-disclosure agreement, expressly prescribing the scope
// and manner of such use.
#ifndef STDTYPES_HH_INCLUDED
#define STDTYPES_HH_INCLUDED
#include <vector>
namespace COMISO_STD
{
// Vector Types
typedef std::vector<double> DoubleVector;
typedef std::vector<int> IntVector;
typedef std::vector<unsigned int> UIntVector;
}//namespace COMISO_GMM
#endif//STDTYPES_HH_INCLUDED
......@@ -12,6 +12,9 @@
//=============================================================================
#include "CBCSolver.hh"
#include <CoMISo/Utils/CoMISoError.hh>
#include <Base/Debug/DebTime.hh>
// For Branch and bound
#include "OsiSolverInterface.hpp"
......@@ -50,6 +53,8 @@ namespace COMISO {
//== IMPLEMENTATION ==========================================================
namespace {
// These are some "tweaks" for the CbcModel, copied from some sample Cbc code
// I have no idea what any of these do!
void model_tweak(CbcModel& model)
......@@ -127,6 +132,7 @@ void model_tweak(CbcModel& model)
//model.addHeuristic(&heuristic2);
}
#define TRACE_CBT(DESCR, EXPR) DEB_line(7, DESCR << ": " << EXPR)
#define P(X) ((X).data())
bool solve_impl(
......@@ -136,10 +142,11 @@ bool solve_impl(
const double _time_limit
)
{
if(!_problem->constant_hessian())
std::cerr << "Warning: CBCSolver received a problem with non-constant hessian!" << std::endl;
if(!_problem->constant_gradient())
std::cerr << "Warning: CBCSolver received a problem with non-constant gradient!" << std::endl;
DEB_enter_func;
DEB_warning_if(!_problem->constant_hessian(), 1,
"CBCSolver received a problem with non-constant hessian!");
DEB_warning_if(!_problem->constant_gradient(), 1,
"CBCSolver received a problem with non-constant gradient!");
const int n_rows = _constraints.size(); // Constraints #
const int n_cols = _problem->n_unknowns(); // Unknowns #
......@@ -162,8 +169,8 @@ bool solve_impl(
for (size_t i = 0; i < _constraints.size(); ++i)
{
if(!_constraints[i]->is_linear())
std::cerr << "Warning: constraint " << i << " is non-linear" << std::endl;
DEB_error_if(!_constraints[i]->is_linear(), "constraint " << i <<
" is non-linear");
NConstraintInterface::SVectorNC gc;
_constraints[i]->eval_gradient(P(x), gc);
......@@ -189,8 +196,8 @@ bool solve_impl(
break;
}
// TRACE_CBT("Constraint " << i << " is of type " <<
// _constraints[i]->constraint_type() << "; RHS val", b);
TRACE_CBT("Constraint " << i << " is of type " <<
_constraints[i]->constraint_type() << "; RHS val", b);
}
//----------------------------------------------
......@@ -199,15 +206,15 @@ bool solve_impl(
std::vector<double> objective(n_cols);
_problem->eval_gradient(P(x), P(objective));
// TRACE_CBT("Objective linear term", objective);
TRACE_CBT("Objective linear term", objective);
const double c = _problem->eval_f(P(x));
// ICGB: Where does objective constant term go in CBC?
// MCM: I could not find this either: It is possible that the entire model
// can be translated to accomodate the constant (if != 0). Ask DB!
// DEB_error_if(c > FLT_EPSILON, "Ignoring a non-zero constant objective term: "
// << c);
// TRACE_CBT("Objective constant term", c);
DEB_error_if(c > FLT_EPSILON, "Ignoring a non-zero constant objective term: "
<< c);
TRACE_CBT("Objective constant term", c);
// CBC Problem initialize
OsiClpSolverInterface si;
......@@ -257,7 +264,7 @@ bool solve_impl(
model.branchAndBound();
const double* solution = model.bestSolution();
// THROW_OUTCOME_if(solution == nullptr, UNSPECIFIED_CBC_EXCEPTION);
//COMISO_THROW_if(solution == nullptr, UNSPECIFIED_CBC_EXCEPTION);
// store if solution available
if(solution != 0)
......@@ -271,6 +278,8 @@ bool solve_impl(
#undef P
}//namespace
bool CBCSolver::solve(
NProblemInterface* _problem,
const std::vector<NConstraintInterface*>& _constraints,
......@@ -278,6 +287,7 @@ bool CBCSolver::solve(
const double _time_limit
)
{
DEB_enter_func;
bool valid_solution = false;
try
{
......@@ -285,8 +295,8 @@ bool CBCSolver::solve(
}
catch (CoinError& ce)
{
std::cerr << "CoinError code = " << ce.message() << std::endl;
// THROW_OUTCOME(UNSPECIFIED_CBC_EXCEPTION);
DEB_warning(1, "CoinError code = " << ce.message() << "]\n");
//COMISO_THROW(UNSPECIFIED_CBC_EXCEPTION);
}
return valid_solution;
}
......
......@@ -4,10 +4,18 @@
//
//=============================================================================
#define COMISO_CPLEXSOLVER_C
//== INCLUDES =================================================================
//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
#include "CPLEXSolver.hh"
#include "LinearConstraint.hh"
#include "BoundConstraint.hh"
#include "ConeConstraint.hh"
#include <Base/Debug/DebTime.hh>
#include <Base/Utils/StopWatch.hh>
#if COMISO_CPLEX_AVAILABLE
//=============================================================================
......@@ -22,6 +30,1278 @@ namespace COMISO {
// ********** SOLVE **************** //
bool
CPLEXSolver::
solve2(NProblemInterface* _problem,
std::vector<NConstraintInterface*>& _constraints,
std::vector<PairIndexVtype>& _discrete_constraints,
const double _time_limit,
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
//----------------------------------------------
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;
}
if(!_silent)
{
std::cerr << "#unknowns : " << _problem->n_unknowns() << std::endl;
std::cerr << "#CPLEX variables: " << _problem->n_unknowns() << 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);
for(unsigned int i=0; i<_constraints.size(); ++i)
{
if(!_constraints[i]->is_linear())
std::cerr << "Warning: CPLEXSolver received a problem with non-linear constraints!!!" << std::endl;
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)
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;
}
}
// 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)
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_limit);
{ // 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();
// if (solution_input_path_.empty())
// {
// if (!problem_env_output_path_.empty())
// {
// std::cout << "Writing problem's environment into file \"" << problem_env_output_path_ << "\"." << std::endl;
// model.getEnv().writeParams(problem_env_output_path_);
// }
// if (!problem_output_path_.empty())
// {
// std::cout << "Writing problem into file \"" << problem_output_path_ << "\"." << std::endl;
// GurobiHelper::outputModelToMpsGz(model, problem_output_path_);
// }
//
// model.optimize();
// }
// else
// {
// std::cout << "Reading solution from file \"" << solution_input_path_ << "\"." << std::endl;
// }
//
//----------------------------------------------
// 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;
}
//-----------------------------------------------------------------------------
bool
CPLEXSolver::
solve(NProblemInterface* _problem,
std::vector<NConstraintInterface*>& _constraints,
std::vector<PairIndexVtype>& _discrete_constraints,
const double _time_limit,
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())
{