Commit b291d649 authored by Martin Marinov's avatar Martin Marinov

Merge remote-tracking branch 'VCI/master' into marinom/merge-from-VCI

parents 27f5b51a fba5e448
......@@ -76,6 +76,8 @@ solve(NProblemInterface* _problem,
// move to next constraint
++n_constraints;
} else {
std::cerr << "Warning: COMISOSolver received a problem with non-equality constraints!!!" << std::endl;
}
// resize matrix to final number of constraints
......
......@@ -80,7 +80,7 @@ private:
} // namespace COMISO
//=============================================================================
// support std vectors
EIGEN_DEFINE_STL_VECTOR_SPECIALIZATION(COMISO::ConeConstraint);
EIGEN_DEFINE_STL_VECTOR_SPECIALIZATION(COMISO::ConeConstraint)
//=============================================================================
#endif // ACG_CONECONSTRAINT_HH defined
//=============================================================================
......
......@@ -19,7 +19,7 @@
#include <iostream>
#include <vector>
#include <gmm/gmm.h>
#include <CoMISo/Utils/gmm.hh>
#include <CoMISo/Config/CoMISoDefines.hh>
#include <CoMISo/NSolver/NConstraintInterface.hh>
......
......@@ -182,6 +182,7 @@ public:
for(unsigned int j=0; j<NV; ++j)
x_[j] = _x[instances_.index(i,j)];
triplets_.clear();
element_.eval_hessian(x_, instances_.c(i), triplets_);
for(unsigned int j=0; j<triplets_.size(); ++j)
......
This diff is collapsed.
......@@ -46,7 +46,9 @@ public:
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
const double _time_limit = 60, // time limit in seconds
const double _gap = 0.0); // stops when solution with optimality gap
// lower than _gab is reached
bool solve(NProblemInterface* _problem, // problem instance
const std::vector<NConstraintInterface*>& _constraints, // linear constraints
......@@ -89,6 +91,14 @@ public:
const double _time_limit = 60,
const bool _silent = false);
bool solve( NProblemInterface* _problem,
const std::vector<NConstraintInterface*>& _constraints,
const std::vector<NConstraintInterface*>& _lazy_constraints,
const std::vector<PairIndexVtype>& _discrete_constraints, // discrete constraints
const double _time_limit = 60,
const double _gap = 0.0,
const int _lazy_level = 0); // lazy level between 0 and 3
// same as above with additional lazy constraints that are only added iteratively to the problem if not satisfied
bool solve(NProblemInterface* _problem,
......@@ -104,24 +114,19 @@ public:
}
void set_problem_output_path ( const std::string &_problem_output_path);
void set_problem_env_output_path( const std::string &_problem_env_output_path);
void set_solution_input_path ( const std::string &_solution_input_path);
void set_problem_output_path ( const std::string &_problem_output_path);
void set_start_solution_output_path( const std::string &_start_solution_output_path);
void set_problem_env_output_path ( const std::string &_problem_env_output_path);
void set_solution_input_path ( const std::string &_solution_input_path);
protected:
double* P(std::vector<double>& _v)
{
if( !_v.empty())
return ((double*)&_v[0]);
else
return 0;
}
private:
// filenames for exporting/importing gurobi solutions
// if string is empty nothing is imported or exported
std::string problem_output_path_;
std::string start_solution_output_path_;
std::string problem_env_output_path_;
std::string solution_input_path_;
};
......
//=============================================================================
//
// STRUCT IPOPTCallbackParameters
//
//=============================================================================
#ifndef COMISO_IPOPTCALLBACKPARAMETERS_HH
#define COMISO_IPOPTCALLBACKPARAMETERS_HH
#include <IpTNLP.hpp>
//== TYPE DEFINITION CALLBACK PARAMETERS ======================================
namespace COMISO {
struct IPOPTCallbackParameters {
Ipopt::AlgorithmMode mode;
Ipopt::Index iter;
Ipopt::Number obj_value;
Ipopt::Number inf_pr;
Ipopt::Number inf_du;
Ipopt::Number mu;
Ipopt::Number d_norm;
Ipopt::Number regularization_size;
Ipopt::Number alpha_du;
Ipopt::Number alpha_pr;
Ipopt::Index ls_trials;
const Ipopt::IpoptData* ip_data;
Ipopt::IpoptCalculatedQuantities* ip_cq;
};
}
#endif
This diff is collapsed.
......@@ -8,134 +8,144 @@
#ifndef COMISO_IPOPTSOLVER_HH
#define COMISO_IPOPTSOLVER_HH
//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
#include <CoMISo/Config/config.hh>
#if COMISO_IPOPT_AVAILABLE
//== 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 "NProblemIPOPT.hh"
#include "NConstraintInterface.hh"
#include "BoundConstraint.hh"
#include <IpTNLP.hpp>
#include <IpIpoptApplication.hpp>
#include <IpSolveStatistics.hpp>
#include <functional>
#include <vector>
#include <string>
//== FORWARDDECLARATIONS ======================================================
#include <CoMISo/Config/CoMISoDefines.hh>
//== NAMESPACES ===============================================================
namespace COMISO {
//== FORWARDDECLARATIONS ======================================================
class NProblemInterface;
class NConstraintInterface;
struct IPOPTCallbackParameters;
//== CLASS DEFINITION =========================================================
/** \class IPOPTSolver
Solver for Interior Point optimization problems.
/** \class NewtonSolver NewtonSolver.hh <ACG/.../NewtonSolver.hh>
Solves an interior point problem, given an NProblemInterface
instance and optionally a set of constraints as well as "lazy
constraints" via NConstraintInterface.
Brief Description.
A more elaborate description follows.
Lazy constraints are not active while the initial solution to the
problem is computed. After the first solution is found, the lazy
constraints are checked and added to the set of active constraints
if they are violated. This process is then repeated until all
constraints are satisfied OR a maximum number of solution attempts
has been reached. In that case the optimization is started once
more, with all lazy constraints active.
*/
class COMISODLLEXPORT IPOPTSolver
{
public:
/// Default constructor -> set up IpOptApplication
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);
// 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
int solve(NProblemInterface* _problem);
// deprecated interface for backwards compatibility
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");
// app().Options()->SetStringValue("hessian_approximation", "limited-memory");
IPOPTSolver();
~IPOPTSolver();
// *********** OPTIONS **************//
/*!
Set options of the underlying ipopt solver.
For a thorough list and documentation of available options, refer
to: https://www.coin-or.org/Ipopt/documentation/node40.html
*/
void set_ipopt_option(std::string, const int&);
void set_ipopt_option(std::string, const double&);
void set_ipopt_option(std::string, const std::string&);
/*!
Get options of the underlying ipopt solver.
The type of option {int, double, std::string} needs to be passed as
template argument.
*/
template<typename T>
T get_ipopt_option(std::string option);
/*!
Set the maximum number of iterations
*/
void set_max_iterations(const int _max_iterations);
int get_max_iterations() const;
/*! Set the threshold on the lazy inequality constraint to decide
if we are near the constraint boundary.
*/
void set_almost_infeasible_threshold(const double _alm_infsb_thrsh);
double get_almost_infeasible_threshold() const;
/*!
Set the max number of incremental lazy constraint iterations before switching
to the fully constrained problem.
\note The default value is 5.
*/
void set_incremental_lazy_constraint_max_iteration_number
(const int _incr_lazy_cnstr_max_iter_nmbr);
int get_incremental_lazy_constraint_max_iteration_number() const;
/*
Turn on/off solving the fully constraint problem after exhausting the
incremental lazy constraint iterations.
\note The default value of this is true.
*/
void set_enable_all_lazy_contraints(const bool _enbl_all_lzy_cnstr);
bool get_enable_all_lazy_contraints() const;
/*
Set intermediate callback function object. For the definition of
IPOPTCallbackParameters include the IPOPTCallbackParameters.hh
header.
If the callback function returns false, IPOPT will terminate
prematurely with the User_Requested_Stop status.
*/
void set_callback_function
(std::function<bool(const IPOPTCallbackParameters &)>);
Ipopt::IpoptApplication& app() {return (*app_); }
// ********** SOLVE **************** //
//! Solve a problem instance with an optional set of constraints.
//! \throws Outcome
void solve
(NProblemInterface* _problem,
const std::vector<NConstraintInterface*>& _constraints = {});
void set_print_level(const int _level)
{
app().Options()->SetIntegerValue("print_level", _level);
print_level_ = _level;
}
//! 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);
protected:
double* P(std::vector<double>& _v)
{
if( !_v.empty())
return ((double*)&_v[0]);
else
return 0;
}
//! Get the computed solution energy
double energy();
private:
class Impl;
Impl* impl_;
// smart pointer to IpoptApplication to set options etc.
Ipopt::SmartPtr<Ipopt::IpoptApplication> app_;
int print_level_;
// inhibit copy
IPOPTSolver(const IPOPTSolver&);
IPOPTSolver& operator=(const IPOPTSolver&);
};
//=============================================================================
} // namespace COMISO
//=============================================================================
#endif // COMISO_IPOPT_AVAILABLE
//=============================================================================
#endif // ACG_IPOPTSOLVER_HH defined
#endif // COMISO_IPOPTSOLVER_HH defined
//=============================================================================
This diff is collapsed.
//=============================================================================
//
// CLASS IPOPTSolverLean
//
//=============================================================================
#ifndef COMISO_IPOPTLEANSOLVER_HH
#define COMISO_IPOPTLEANSOLVER_HH
//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
#include <CoMISo/Config/config.hh>
#if COMISO_IPOPT_AVAILABLE
//== INCLUDES =================================================================
#include <CoMISo/Config/CoMISoDefines.hh>
#include <vector>
#include <cstddef>
//== NAMESPACES ===============================================================
namespace COMISO {
//== FORWARDDECLARATIONS ======================================================
class NProblemGmmInterface; // deprecated
class NProblemInterface;
class NConstraintInterface;
//== CLASS DEFINITION =========================================================
/** \class NewtonSolver NewtonSolver.hh <ACG/.../NewtonSolver.hh>
Brief Description.
A more elaborate description follows.
*/
class COMISODLLEXPORT IPOPTSolverLean
{
public:
/// Default constructor
IPOPTSolverLean();
/// Destructor
~IPOPTSolverLean();
// *********** OPTIONS **************//
/*!
Set the maximum number of iterations
*/
void set_max_iterations(const int _max_iterations);
int max_iterations() const;
/*!
Set the threshold on the lazy inequality constraint to decide if we are near
the constraint boundary.
*/
void set_almost_infeasible_threshold(const double _alm_infsb_thrsh);
double almost_infeasible_threshold() const;
/*!
Set the max number of incremental lazy constraint iterations before switching
to the fully constrained problem.
\note The default value is 5.
*/
void set_incremental_lazy_constraint_max_iteration_number(const int
_incr_lazy_cnstr_max_iter_nmbr);
int incremental_lazy_constraint_max_iteration_number() const;
/*
Turn on/off solving the fully constraint problem after exhausting the
incremental lazy constraint iterations.
\note The default value of this is true.
*/
void set_enable_all_lazy_contraints(const bool _enbl_all_lzy_cnstr);
bool enable_all_lazy_contraints() const;
// ********** 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);
// for convenience, if no constraints are given
//! \throws Outcome
void solve(NProblemInterface* _problem);
// deprecated interface for backwards compatibility
//! \throws Outcome
void solve(NProblemGmmInterface* _problem,
std::vector<NConstraintInterface*>& _constraints);
//! Get the computed solution energy
double energy();
private:
class Impl;
Impl* impl_;
// inhibit copy
IPOPTSolverLean(const IPOPTSolverLean&);
IPOPTSolverLean& operator=(const IPOPTSolverLean&);
};
//=============================================================================
} // namespace COMISO
//=============================================================================
#endif // COMISO_IPOPT_AVAILABLE
//=============================================================================
#endif // ACG_IPOPTSOLVER_HH defined
//=============================================================================
......@@ -88,7 +88,7 @@ private:
} // namespace COMISO
//=============================================================================
// support std vectors
EIGEN_DEFINE_STL_VECTOR_SPECIALIZATION(COMISO::LinearConstraint);
EIGEN_DEFINE_STL_VECTOR_SPECIALIZATION(COMISO::LinearConstraint)
//=============================================================================
#endif // ACG_LINEARCONSTRAINT_HH defined
//=============================================================================
......
......@@ -13,7 +13,7 @@
#include <CoMISo/Config/CoMISoDefines.hh>
#include <iostream>
#include <gmm/gmm.h>
#include <CoMISo/Utils/gmm.hh>
//== FORWARDDECLARATIONS ======================================================
......
......@@ -12,8 +12,8 @@
//== INCLUDES =================================================================
#include <CoMISo/Config/CoMISoDefines.hh>
#include <CoMISo/Utils/gmm.hh>
#include <iostream>
#include <gmm/gmm.h>
//== FORWARDDECLARATIONS ======================================================
......
......@@ -23,7 +23,7 @@
#include <climits>
#include <CoMISo/Utils/VSToolsT.hh>
#include <gmm/gmm.h>
#include <CoMISo/Utils/gmm.hh>
#include <CoMISo/Config/CoMISoDefines.hh>
......
......@@ -12,7 +12,7 @@
//== INCLUDES =================================================================
#include <CoMISo/Config/CoMISoDefines.hh>
#include <gmm/gmm.h>
#include <CoMISo/Utils/gmm.hh>
#include "NProblemGmmInterface.hh"
#include "LinearConstraintHandlerElimination.hh"
#include "LinearConstraintHandlerPenalty.hh"
......
......@@ -9,7 +9,7 @@
//== INCLUDES =================================================================
#include <CoMISo/Config/CoMISoDefines.hh>
#include <gmm/gmm.h>
#include <CoMISo/Utils/gmm.hh>
#include "NProblemGmmInterface.hh"
#include "LinearConstraintHandlerElimination.hh"
#include "LinearConstraintHandlerPenalty.hh"
......
......@@ -12,7 +12,7 @@
//== INCLUDES =================================================================
#include <Base/Utils/StopWatch.hh>
#include <gmm/gmm.h>
#include <CoMISo/Utils/gmm.hh>
#include "NProblemInterface.hh"
#include <CoMISo/Config/CoMISoDefines.hh>
......
......@@ -11,7 +11,7 @@
//== INCLUDES =================================================================
#include <gmm/gmm.h>
#include <CoMISo/Utils/gmm.hh>
#include <CoMISo/Config/CoMISoDefines.hh>
#include <Base/Debug/DebOut.hh>
......
......@@ -23,6 +23,7 @@ solve(NProblemGmmInterface* _problem)
{
DEB_enter_func;
#if COMISO_SUITESPARSE_AVAILABLE
converged_ = true;
// get problem size
int n = _problem->n_unknowns();
......@@ -98,16 +99,19 @@ solve(NProblemGmmInterface* _problem)
_problem->store_result(P(x));
DEB_line(2, "Newton solver reached max regularization but did not "
"converge");
converged_ = false;
return false;
}
}
}
_problem->store_result(P(x));
DEB_line(2, "Newton Solver did not converge!!! after iterations.");
converged_ = false;
return false;
#else
DEB_warning(1,"NewtonSolver requires not-available CholmodSolver");
converged_ = false;
return false;
#endif
}
......@@ -120,6 +124,7 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
const VectorD& _b)
{
DEB_time_func_def;
converged_ = false;
const double KKT_res_eps = 1e-6;
const int max_KKT_regularization_iters = 40;
......@@ -160,10 +165,11 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
double kkt_res2(0.0);
double constraint_res2(0.0);
int reg_iters(0);
bool fact_ok = true;
do
{
// get Newton search direction by solving LSE
bool fact_ok = factorize(_problem, _A, _b, x, regularize_hessian, regularize_constraints, first_factorization);
fact_ok = factorize(_problem, _A, _b, x, regularize_hessian, regularize_constraints, first_factorization);
first_factorization = false;
if(fact_ok)
......@@ -204,7 +210,7 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
}
++reg_iters;
}
while( (kkt_res2 > KKT_res_eps || constraint_res2 > max_allowed_constraint_violation2) && reg_iters < max_KKT_regularization_iters);
while( (!fact_ok || kkt_res2 > KKT_res_eps || constraint_res2 > max_allowed_constraint_violation2) && reg_iters < max_KKT_regularization_iters);
// no valid step could be found?
if(kkt_res2 > KKT_res_eps || constraint_res2 > max_allowed_constraint_violation2 || reg_iters >= max_KKT_regularization_iters)
......@@ -250,8 +256,11 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
<< ", KKT residual^2 = " << kkt_res2);
// converged?
if(newton_decrement < eps_ || std::abs(t) < eps_ls_)
if(newton_decrement < eps_ || std::abs(t) < eps_ls_ )
{
converged_ = true;
break;
}
++iter;
}
......@@ -273,7 +282,7 @@ bool NewtonSolver::factorize(NProblemInterface* _problem,
{
DEB_enter_func;
const int n = _problem->n_unknowns();
const int m = _A.rows();
const int m = static_cast<int>(_A.rows());
const int nf = n+m;
// get hessian of quadratic problem
......@@ -288,17 +297,17 @@ bool NewtonSolver::factorize(NProblemInterface* _problem,
// add elements of H
for (int k=0; k<H.outerSize(); ++k)
for (SMatrixD::InnerIterator it(H,k); it; ++it)
trips.push_back(Triplet(it.row(),it.col(),it.value()));
trips.push_back(Triplet(static_cast<int>(it.row()),static_cast<int>(it.col()),it.value()));
// add elements of _A
for (int k=0; k<_A.outerSize(); ++k)
for (SMatrixD::InnerIterator it(_A,k); it; ++it)
{
// insert _A block below
trips.push_back(Triplet(it.row()+n,it.col(),it.value()));
trips.push_back(Triplet(static_cast<int>(it.row())+n,static_cast<int>(it.col()),it.value()));
// insert _A^T block right
trips.push_back(Triplet(it.col(),it.row()+n,it.value()));
trips.push_back(Triplet(static_cast<int>(it.col()),static_cast<int>(it.row())+n,it.value()));
}
// regularize constraints
......@@ -397,11 +406,11 @@ bool NewtonSolver::numerical_factorization(SMatrixD& _KKT)
DEB_enter_func;
switch(solver_type_)
{
case LS_EigenLU:
case LS_EigenLU:
lu_solver_.factorize(_KKT);
return (lu_solver_.info() == Eigen::Success);
#if COMISO_SUITESPARSE_AVAILABLE
case LS_Umfpack:
case LS_Umfpack:
umfpack_solver_.factorize(_KKT);
return (umfpack_solver_.info() == Eigen::Success</