Commit 578d2d21 authored by Martin Marinov's avatar Martin Marinov
Browse files

Added an Outcome-based exception machinery for unrecoverable error reporting.

[git-p4: depot-paths = "//ReForm/ReForm/main/CoMISo/": change = 10406]
parent b58c436f
......@@ -16,6 +16,7 @@
#include "GUROBISolver.hh"
#include <CoMISo/Utils/StopWatch.hh>
#include <Base/OutcomeUtils.hh>
#include <stdexcept>
//== NAMESPACES ===============================================================
......@@ -46,18 +47,12 @@ GUROBISolver()
//#define TRACE_GUROBI_INPUT(DESCR, EXPR)
// ********** SOLVE **************** //
bool
GUROBISolver::
solve(NProblemInterface* _problem,
std::vector<NConstraintInterface*>& _constraints,
std::vector<PairIndexVtype>& _discrete_constraints,
const double /*_time_limit*/)
void GUROBISolver::solve(
NProblemInterface* _problem,
std::vector<NConstraintInterface*>& _constraints,
std::vector<PairIndexVtype>& _discrete_constraints,
const double /*_time_limit*/)
{
//std::cout << "_time_limit: " << _time_limit << std::endl;
MY_PAUSE;
//std::cout.precision(std::numeric_limits< double >::digits10);
try
{
//----------------------------------------------
......@@ -69,8 +64,8 @@ solve(NProblemInterface* _problem,
//model.getEnv().set(GRB_DoubleParam_TimeLimit, _time_limit);
// stop when a solution is found
model.getEnv().set(GRB_IntParam_SolutionLimit, 2);
// stop when a solution is found
model.getEnv().set(GRB_IntParam_SolutionLimit, 2);
//----------------------------------------------
// 1. allocate variables
......@@ -110,7 +105,7 @@ solve(NProblemInterface* _problem,
for(unsigned int i=0; i<_constraints.size(); ++i)
{
TRACE_GUROBI_INPUT("Constraint index: ", i);
TRACE_GUROBI_INPUT("Constraint index: ", i);
if(!_constraints[i]->is_linear())
std::cerr << "Warning: GUROBISolver received a problem with non-linear constraints!!!" << std::endl;
......@@ -118,26 +113,26 @@ solve(NProblemInterface* _problem,
NConstraintInterface::SVectorNC gc;
_constraints[i]->eval_gradient(P(x), gc);
TRACE_GUROBI_INPUT("Constraint type: ",
(int)_constraints[i]->constraint_type());
TRACE_GUROBI_INPUT("Constraint type: ",
(int)_constraints[i]->constraint_type());
NConstraintInterface::SVectorNC::InnerIterator v_it(gc);
for(; v_it; ++v_it)
{
// lin_expr += v_it.value()*vars[v_it.index()];
{
// lin_expr += v_it.value()*vars[v_it.index()];
lin_expr += vars[v_it.index()]*_QNT(v_it.value());
TRACE_GUROBI_INPUT("Constraint linear coeff", _QNT(v_it.value()));
TRACE_GUROBI_INPUT("Constraint linear var idx", v_it.index());
}
TRACE_GUROBI_INPUT("Constraint linear coeff", _QNT(v_it.value()));
TRACE_GUROBI_INPUT("Constraint linear var idx", v_it.index());
}
double b = _QNT(_constraints[i]->eval_constraint(P(x)));
TRACE_GUROBI_INPUT("Constraint B", b);
TRACE_GUROBI_INPUT("Constraint B", b);
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();
......@@ -155,29 +150,29 @@ solve(NProblemInterface* _problem,
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*_QNT(it.value())*vars[it.row()]*vars[it.col()];
TRACE_GUROBI_INPUT("Objective quadratic term coeff", _QNT(it.value()));
TRACE_GUROBI_INPUT("Objective quadratic term var i idx", it.row());
TRACE_GUROBI_INPUT("Objective quadratic term var j idx", it.col());
}
}
TRACE_GUROBI_INPUT("Objective quadratic term coeff", _QNT(it.value()));
TRACE_GUROBI_INPUT("Objective quadratic term var i idx", it.row());
TRACE_GUROBI_INPUT("Objective quadratic term var j idx", 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 += _QNT(g[i])*vars[i];
TRACE_GUROBI_INPUT("Objective linear term coeff", _QNT(g[i]));
TRACE_GUROBI_INPUT("Objective linear term var idx", i);
}
TRACE_GUROBI_INPUT("Objective linear term coeff", _QNT(g[i]));
TRACE_GUROBI_INPUT("Objective linear term var idx", i);
}
// add constant part
const double c = _QNT(_problem->eval_f(P(x)));
const double c = _QNT(_problem->eval_f(P(x)));
objective += c;
TRACE_GUROBI_INPUT("Objective constant term coeff", c);
TRACE_GUROBI_INPUT("Objective constant term coeff", c);
model.set(GRB_IntAttr_ModelSense, 1);
model.setObjective(objective);
......@@ -206,7 +201,7 @@ solve(NProblemInterface* _problem,
}
else
{
std::cout << "Reading solution from file \"" << solution_input_path_ << "\"." << std::endl;
std::cout << "Reading solution from file \"" << solution_input_path_ << "\"." << std::endl;
}
//----------------------------------------------
......@@ -236,30 +231,27 @@ solve(NProblemInterface* _problem,
// 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 true;
std::cout << "GUROBI Objective: " << model.get(GRB_DoubleAttr_ObjVal) << std::endl;
}
catch(GRBException& e)
catch (GRBException& 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;
}
MY_PAUSE;
return false;
// NOTE: we could propagate e.getMessage() either using std::exception, or a specialized Reform exception type
switch ( e.getErrorCode() )
{
// TODO: Find where these are defined and use the correct name for 10010
case 10010: THROW_OUTCOME(GUROBI_LICENCE_MODEL_TOO_LARGE);
default: THROW_OUTCOME(UNSPECIFIED_GUROBI_EXCEPTION);
}
}
}
//-----------------------------------------------------------------------------
#if 0
bool
GUROBISolver::
solve(NProblemInterface* _problem,
......@@ -507,7 +499,7 @@ solve(NProblemInterface* _problem,
return false;
}
catch(...)
{
{
std::cout << "Exception during optimization" << std::endl;
return false;
}
......@@ -666,6 +658,7 @@ solve(NProblemInterface* _problem,
// return (solution_found != IloFalse);
}
#endif//0
//-----------------------------------------------------------------------------
......
......@@ -23,8 +23,6 @@
#include "VariableType.hh"
#include "GurobiHelper.hh"
//#include <gurobi_c++.h>
//== FORWARDDECLARATIONS ======================================================
class GRBModel;
class GRBVar;
......@@ -54,42 +52,49 @@ public:
~GUROBISolver() {}
// ********** SOLVE **************** //
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
bool solve(NProblemInterface* _problem, // problem instance
std::vector<NConstraintInterface*>& _constraints, // linear constraints
const double _time_limit = 60 ) // time limit in seconds
//! \throws Outcome
void 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
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);
std::vector<PairIndexVtype> dc;
return solve(_problem, _constraints, dc, _time_limit);
}
#if 0
// optimization with additional lazy constraints that are only added iteratively to the problem if not satisfied
bool solve(NProblemInterface* _problem,
const std::vector<NConstraintInterface*>& _constraints,
const std::vector<NConstraintInterface*>& _lazy_constraints,
std::vector<PairIndexVtype>& _discrete_constraints, // discrete constraints
const double _almost_infeasible = 0.5,
const int _max_passes = 5,
const double _time_limit = 60,
const bool _silent = false);
const std::vector<NConstraintInterface*>& _constraints,
const std::vector<NConstraintInterface*>& _lazy_constraints,
std::vector<PairIndexVtype>& _discrete_constraints, // discrete constraints
const double _almost_infeasible = 0.5,
const int _max_passes = 5,
const double _time_limit = 60,
const bool _silent = false);
// same as above with additional lazy constraints that are only added iteratively to the problem if not satisfied
bool 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,
const double _time_limit = 60,
const bool _silent = false)
const std::vector<NConstraintInterface*>& _constraints,
const std::vector<NConstraintInterface*>& _lazy_constraints,
const double _almost_infeasible = 0.5,
const int _max_passes = 5,
const double _time_limit = 60,
const bool _silent = false)
{
std::vector<PairIndexVtype> dc; return solve(_problem, _constraints, _lazy_constraints, dc, _almost_infeasible, _max_passes, _time_limit, _silent);
std::vector<PairIndexVtype> dc;
return solve(_problem, _constraints, _lazy_constraints, dc, _almost_infeasible, _max_passes, _time_limit, _silent);
}
#endif//0
void set_problem_output_path ( const std::string &_problem_output_path);
void set_problem_env_output_path( const std::string &_problem_env_output_path);
......
This diff is collapsed.
......@@ -24,9 +24,6 @@
#include "NProblemInterface.hh"
#include "NConstraintInterface.hh"
#include "BoundConstraint.hh"
#include <IpTNLP.hpp>
#include <IpIpoptApplication.hpp>
#include <IpSolveStatistics.hpp>
//== FORWARDDECLARATIONS ======================================================
......@@ -36,7 +33,6 @@ namespace COMISO {
//== CLASS DEFINITION =========================================================
/** \class NewtonSolver NewtonSolver.hh <ACG/.../NewtonSolver.hh>
......@@ -48,191 +44,38 @@ class COMISODLLEXPORT IPOPTSolver
{
public:
/// Default constructor -> set up IpOptApplication
/// Default constructor
IPOPTSolver();
/// Destructor
~IPOPTSolver() {}
~IPOPTSolver();
// ********** SOLVE **************** //
// 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);
//! \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
int solve(NProblemInterface* _problem,
//! \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 );
// for convenience, if no constraints are given
int solve(NProblemInterface* _problem);
//! \throws Outcome
void solve(NProblemInterface* _problem);
// deprecated interface for backwards compatibility
int solve(NProblemGmmInterface* _problem, std::vector<NConstraintInterface*>& _constraints);
//! \throws Outcome
void 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_); }
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_;
};
//== 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)
: problem_(_problem), store_solution_(false) { 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_;}
private:
/**@name Methods to block default compiler methods.
* The compiler automatically generates the following three methods.
* Since the default compiler implementation is generally not what
* you want (for all but the most simple classes), we usually
* put the declarations of these methods in the private section
* and never implement them. This prevents the compiler from
* implementing an incorrect "default" behavior without us
* knowing. (See Scott Meyers book, "Effective C++")
*
*/
//@{
// MyNLP();
NProblemIPOPT(const NProblemIPOPT&);
NProblemIPOPT& operator=(const NProblemIPOPT&);
//@}
// split user-provided constraints into general-constraints and bound-constraints
void split_constraints(const std::vector<NConstraintInterface*>& _constraints);
// determine if hessian_constant, jac_c_constant or jac_d_constant
void analyze_special_properties(const NProblemInterface* _problem, const std::vector<NConstraintInterface*>& _constraints);
//Ipopt::IpoptApplication& app() {return (*app_); }
protected:
......@@ -245,150 +88,8 @@ protected:
}
private:
// pointer to problem instance
NProblemInterface* problem_;
// reference to constraints vector
std::vector<NConstraintInterface*> constraints_;
std::vector<BoundConstraint*> bound_constraints_;
bool hessian_constant_;
bool jac_c_constant_;
bool jac_d_constant_;
bool store_solution_;
std::vector<double> x_;
};
//== CLASS DEFINITION PROBLEM INSTANCE=========================================================
class NProblemGmmIPOPT : 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 gmm::wsvector<double> SVectorNP;
typedef NProblemGmmInterface::SMatrixNP SMatrixNP;
typedef gmm::array1D_reference< double* > VectorPT;
typedef gmm::array1D_reference< const double* > VectorPTC;
typedef gmm::array1D_reference< Index* > VectorPTi;
typedef gmm::array1D_reference< const Index* > VectorPTCi;
typedef gmm::linalg_traits<SVectorNP>::const_iterator SVectorNP_citer;
typedef gmm::linalg_traits<SVectorNP>::iterator SVectorNP_iter;
/** default constructor */
NProblemGmmIPOPT(NProblemGmmInterface* _problem, std::vector<NConstraintInterface*>& _constraints)
: problem_(_problem), constraints_(_constraints), nnz_jac_g_(0), nnz_h_lag_(0)
{}
/** default destructor */
virtual ~NProblemGmmIPOPT() {};
/**@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);
//@}
private:
/**@name Methods to block default compiler methods.
* The compiler automatically generates the following three methods.
* Since the default compiler implementation is generally not what
* you want (for all but the most simple classes), we usually