Commit 31037bdd authored by David Bommes's avatar David Bommes
Browse files

added Eigen3 support

git-svn-id: http://www.openflipper.org/svnrepo/CoMISo/trunk@67 1355f012-dd97-4b2f-ae87-10fa9f823a57
parent 75338fad
...@@ -97,6 +97,15 @@ else () ...@@ -97,6 +97,15 @@ else ()
set (COMISO_IPOPT_CONFIG_FILE_SETTINGS "#define COMISO_IPOPT_AVAILABLE 0" ) set (COMISO_IPOPT_CONFIG_FILE_SETTINGS "#define COMISO_IPOPT_AVAILABLE 0" )
endif () endif ()
find_package (Eigen3)
if (Eigen3_FOUND )
set (COMISO_Eigen3_CONFIG_FILE_SETTINGS "#define COMISO_Eigen3_AVAILABLE 1" )
list( APPEND COMISO_INCLUDE_DIRECTORIES ${Eigen3_INCLUDE_DIR} )
else ()
message (STATUS "Eigen3 not found!")
set (COMISO_Eigen3_CONFIG_FILE_SETTINGS "#define COMISO_Eigen3_AVAILABLE 0" )
endif ()
if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/factored_solver/CMakeLists.txt" ) if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/factored_solver/CMakeLists.txt" )
add_subdirectory (Examples/factored_solver) add_subdirectory (Examples/factored_solver)
endif() endif()
...@@ -113,6 +122,9 @@ endif() ...@@ -113,6 +122,9 @@ endif()
if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_factored_example/CMakeLists.txt" ) if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_factored_example/CMakeLists.txt" )
add_subdirectory (Examples/small_factored_example) add_subdirectory (Examples/small_factored_example)
endif() endif()
if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/super_sparse_matrix/CMakeLists.txt" )
add_subdirectory (Examples/super_sparse_matrix)
endif()
include_directories ( include_directories (
.. ..
......
...@@ -24,7 +24,7 @@ namespace COMISO { ...@@ -24,7 +24,7 @@ namespace COMISO {
int int
IPOPTSolver:: IPOPTSolver::
solve(NProblemGmmInterface* _problem, std::vector<NConstraintGmmInterface*>& _constraints) solve(NProblemGmmInterface* _problem, std::vector<NConstraintInterface*>& _constraints)
{ {
//---------------------------------------------------------------------------- //----------------------------------------------------------------------------
// 1. Create an instance of IPOPT NLP // 1. Create an instance of IPOPT NLP
...@@ -38,7 +38,7 @@ solve(NProblemGmmInterface* _problem, std::vector<NConstraintGmmInterface*>& _co ...@@ -38,7 +38,7 @@ solve(NProblemGmmInterface* _problem, std::vector<NConstraintGmmInterface*>& _co
Ipopt::SmartPtr<Ipopt::IpoptApplication> app = IpoptApplicationFactory(); Ipopt::SmartPtr<Ipopt::IpoptApplication> app = IpoptApplicationFactory();
app->Options()->SetStringValue("linear_solver", "ma57"); app->Options()->SetStringValue("linear_solver", "ma57");
// app->Options()->SetStringValue("derivative_test", "second-order"); // app->Options()->SetStringValue("derivative_test", "second-order");
// app->Options()->SetIntegerValue("print_level", 0); // app->Options()->SetIntegerValue("print_level", 0);
// Initialize the IpoptApplication and process the options // Initialize the IpoptApplication and process the options
...@@ -89,11 +89,12 @@ bool NProblemIPOPT::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g, ...@@ -89,11 +89,12 @@ bool NProblemIPOPT::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
// ToDo: perturb x // ToDo: perturb x
// nonzeros in the jacobian of C_ and the hessian of the lagrangian // nonzeros in the jacobian of C_ and the hessian of the lagrangian
SVectorNP g; SMatrixNP HP;
SMatrixNP H; SVectorNC g;
problem_->eval_hessian(&(x[0]), H); SMatrixNC H;
problem_->eval_hessian(&(x[0]), HP);
nnz_jac_g = 0; nnz_jac_g = 0;
nnz_h_lag = gmm::nnz(H); nnz_h_lag = 0;
// clear old data // clear old data
jac_g_iRow_.clear(); jac_g_iRow_.clear();
...@@ -105,7 +106,7 @@ bool NProblemIPOPT::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g, ...@@ -105,7 +106,7 @@ bool NProblemIPOPT::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
// iterate over rows // iterate over rows
for( int i=0; i<n; ++i) for( int i=0; i<n; ++i)
{ {
SVectorNP& ri = H.row(i); SVectorNP& ri = HP.row(i);
SVectorNP_citer v_it = gmm::vect_const_begin(ri); SVectorNP_citer v_it = gmm::vect_const_begin(ri);
SVectorNP_citer v_end = gmm::vect_const_end (ri); SVectorNP_citer v_end = gmm::vect_const_end (ri);
...@@ -117,6 +118,7 @@ bool NProblemIPOPT::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g, ...@@ -117,6 +118,7 @@ bool NProblemIPOPT::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
{ {
h_lag_iRow_.push_back(i); h_lag_iRow_.push_back(i);
h_lag_jCol_.push_back(v_it.index()); h_lag_jCol_.push_back(v_it.index());
++nnz_h_lag;
} }
} }
} }
...@@ -127,36 +129,26 @@ bool NProblemIPOPT::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g, ...@@ -127,36 +129,26 @@ bool NProblemIPOPT::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
{ {
constraints_[i]->eval_gradient(&(x[0]),g); constraints_[i]->eval_gradient(&(x[0]),g);
constraints_[i]->eval_hessian (&(x[0]),H); constraints_[i]->eval_hessian (&(x[0]),H);
nnz_jac_g += gmm::nnz(g);
nnz_h_lag += gmm::nnz(H);
// build table // iterate over sparse vector
SVectorNP_citer v_it = gmm::vect_const_begin(g); SVectorNC::InnerIterator v_it(g);
SVectorNP_citer v_end = gmm::vect_const_end (g); for(; v_it; ++v_it)
for(; v_it != v_end; ++v_it)
{ {
jac_g_iRow_.push_back(i); jac_g_iRow_.push_back(i);
jac_g_jCol_.push_back(v_it.index()); jac_g_jCol_.push_back(v_it.index());
++nnz_jac_g;
} }
for( int i=0; i<n; ++i) // iterate over superSparseMatrix
{ SMatrixNC::iterator m_it = H.begin();
SVectorNP& ri = H.row(i); SMatrixNC::iterator m_end = H.end();
for(; m_it != m_end; ++m_it)
v_it = gmm::vect_const_begin(ri); if( m_it.row() >= m_it.col())
v_end = gmm::vect_const_end (ri);
for(; v_it != v_end; ++v_it)
{ {
// store lower triangular part only h_lag_iRow_.push_back(m_it.row());
if( i >= (int)v_it.index()) h_lag_jCol_.push_back(m_it.col());
{ ++nnz_h_lag;
h_lag_iRow_.push_back(i);
h_lag_jCol_.push_back(v_it.index());
}
} }
}
} }
// store for error checking... // store for error checking...
...@@ -192,10 +184,10 @@ bool NProblemIPOPT::get_bounds_info(Index n, Number* x_l, Number* x_u, ...@@ -192,10 +184,10 @@ bool NProblemIPOPT::get_bounds_info(Index n, Number* x_l, Number* x_u,
// enum ConstraintType {NC_EQUAL, NC_LESS_EQUAL, NC_GREATER_EQUAL}; // enum ConstraintType {NC_EQUAL, NC_LESS_EQUAL, NC_GREATER_EQUAL};
switch(constraints_[i]->constraint_type()) switch(constraints_[i]->constraint_type())
{ {
case NConstraintGmmInterface::NC_EQUAL : g_u[i] = 0.0 ; g_l[i] = 0.0 ; break; case NConstraintInterface::NC_EQUAL : g_u[i] = 0.0 ; g_l[i] = 0.0 ; break;
case NConstraintGmmInterface::NC_LESS_EQUAL : g_u[i] = 0.0 ; g_l[i] = -1.0e19; break; case NConstraintInterface::NC_LESS_EQUAL : g_u[i] = 0.0 ; g_l[i] = -1.0e19; break;
case NConstraintGmmInterface::NC_GREATER_EQUAL : g_u[i] = 1.0e19; g_l[i] = 0.0 ; break; case NConstraintInterface::NC_GREATER_EQUAL : g_u[i] = 1.0e19; g_l[i] = 0.0 ; break;
default : g_u[i] = 1.0e19; g_l[i] = -1.0e19; break; default : g_u[i] = 1.0e19; g_l[i] = -1.0e19; break;
} }
} }
...@@ -273,20 +265,18 @@ bool NProblemIPOPT::eval_jac_g(Index n, const Number* x, bool new_x, ...@@ -273,20 +265,18 @@ bool NProblemIPOPT::eval_jac_g(Index n, const Number* x, bool new_x,
// return the structure of the jacobian of the constraints // return the structure of the jacobian of the constraints
// global index // global index
int gi = 0; int gi = 0;
SVectorNP g; SVectorNC g;
for( int i=0; i<m; ++i) for( int i=0; i<m; ++i)
{ {
constraints_[i]->eval_gradient(x, g); constraints_[i]->eval_gradient(x, g);
// iterate over non-zero values SVectorNC::InnerIterator v_it(g);
SVectorNP_citer it = gmm::vect_const_begin(g);
SVectorNP_citer ite = gmm::vect_const_end(g);
for (; it != ite; ++it) for( ; v_it; ++v_it)
{ {
if(gi < nele_jac) if(gi < nele_jac)
values[gi] = *it; values[gi] = v_it.value();
++gi; ++gi;
} }
} }
...@@ -322,12 +312,11 @@ bool NProblemIPOPT::eval_h(Index n, const Number* x, bool new_x, ...@@ -322,12 +312,11 @@ bool NProblemIPOPT::eval_h(Index n, const Number* x, bool new_x,
int gi = 0; int gi = 0;
// get hessian of problem // get hessian of problem
SMatrixNP H; problem_->eval_hessian(x, HP_);
problem_->eval_hessian(x, H);
for( int i=0; i<n; ++i) for( int i=0; i<n; ++i)
{ {
SVectorNP& ri = H.row(i); SVectorNP& ri = HP_.row(i);
SVectorNP_citer v_it = gmm::vect_const_begin(ri); SVectorNP_citer v_it = gmm::vect_const_begin(ri);
SVectorNP_citer v_end = gmm::vect_const_end (ri); SVectorNP_citer v_end = gmm::vect_const_end (ri);
...@@ -347,24 +336,20 @@ bool NProblemIPOPT::eval_h(Index n, const Number* x, bool new_x, ...@@ -347,24 +336,20 @@ bool NProblemIPOPT::eval_h(Index n, const Number* x, bool new_x,
// Hessians of Constraints // Hessians of Constraints
for(unsigned int j=0; j<constraints_.size(); ++j) for(unsigned int j=0; j<constraints_.size(); ++j)
{ {
SMatrixNC H;
constraints_[j]->eval_hessian(x, H); constraints_[j]->eval_hessian(x, H);
for( int i=0; i<n; ++i) SMatrixNC::iterator m_it = H.begin();
{ SMatrixNC::iterator m_end = H.end();
SVectorNP& ri = H.row(i);
SVectorNP_citer v_it = gmm::vect_const_begin(ri); for(; m_it != m_end; ++m_it)
SVectorNP_citer v_end = gmm::vect_const_end (ri); {
// store lower triangular part only
for(; v_it != v_end; ++v_it) if( m_it.row() >= m_it.col())
{ {
// store lower triangular part only if( gi < nele_hess)
if( i >= (int)v_it.index()) values[gi] = lambda[j]*(*m_it);
{ ++gi;
if( gi < nele_hess)
values[gi] = lambda[j]*(*v_it);
++gi;
}
} }
} }
} }
...@@ -398,4 +383,4 @@ void NProblemIPOPT::finalize_solution(SolverReturn status, ...@@ -398,4 +383,4 @@ void NProblemIPOPT::finalize_solution(SolverReturn status,
} // namespace COMISO } // namespace COMISO
//============================================================================= //=============================================================================
#endif // COMISO_IPOPT_AVAILABLE #endif // COMISO_IPOPT_AVAILABLE
//============================================================================= //=============================================================================
\ No newline at end of file
...@@ -19,7 +19,7 @@ ...@@ -19,7 +19,7 @@
#include <cstddef> #include <cstddef>
#include <gmm/gmm.h> #include <gmm/gmm.h>
#include "NProblemGmmInterface.hh" #include "NProblemGmmInterface.hh"
#include "NConstraintGmmInterface.hh" #include "NConstraintInterface.hh"
#include <IpTNLP.hpp> #include <IpTNLP.hpp>
#include <IpIpoptApplication.hpp> #include <IpIpoptApplication.hpp>
#include <IpSolveStatistics.hpp> #include <IpSolveStatistics.hpp>
...@@ -78,7 +78,7 @@ public: ...@@ -78,7 +78,7 @@ public:
// }; // };
//------------------------------------------------------ //------------------------------------------------------
int solve(NProblemGmmInterface* _problem, std::vector<NConstraintGmmInterface*>& _constraints); int solve(NProblemGmmInterface* _problem, std::vector<NConstraintInterface*>& _constraints);
protected: protected:
double* P(std::vector<double>& _v) double* P(std::vector<double>& _v)
...@@ -108,8 +108,11 @@ public: ...@@ -108,8 +108,11 @@ public:
typedef Ipopt::IpoptData IpoptData; typedef Ipopt::IpoptData IpoptData;
typedef Ipopt::IpoptCalculatedQuantities IpoptCalculatedQuantities; typedef Ipopt::IpoptCalculatedQuantities IpoptCalculatedQuantities;
typedef NConstraintGmmInterface::SVectorNP SVectorNP; // sparse matrix and vector types
typedef NConstraintGmmInterface::SMatrixNP SMatrixNP; 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< double* > VectorPT;
typedef gmm::array1D_reference< const double* > VectorPTC; typedef gmm::array1D_reference< const double* > VectorPTC;
...@@ -121,7 +124,7 @@ public: ...@@ -121,7 +124,7 @@ public:
typedef typename gmm::linalg_traits<SVectorNP>::iterator SVectorNP_iter; typedef typename gmm::linalg_traits<SVectorNP>::iterator SVectorNP_iter;
/** default constructor */ /** default constructor */
NProblemIPOPT(NProblemGmmInterface* _problem, std::vector<NConstraintGmmInterface*>& _constraints) NProblemIPOPT(NProblemGmmInterface* _problem, std::vector<NConstraintInterface*>& _constraints)
: problem_(_problem), constraints_(_constraints) {} : problem_(_problem), constraints_(_constraints) {}
/** default destructor */ /** default destructor */
...@@ -205,7 +208,7 @@ private: ...@@ -205,7 +208,7 @@ private:
// pointer to problem instance // pointer to problem instance
NProblemGmmInterface* problem_; NProblemGmmInterface* problem_;
// reference to constraints vector // reference to constraints vector
std::vector<NConstraintGmmInterface*>& constraints_; std::vector<NConstraintInterface*>& constraints_;
int nnz_jac_g_; int nnz_jac_g_;
int nnz_h_lag_; int nnz_h_lag_;
...@@ -215,6 +218,9 @@ private: ...@@ -215,6 +218,9 @@ private:
std::vector<Index> jac_g_jCol_; std::vector<Index> jac_g_jCol_;
std::vector<Index> h_lag_iRow_; std::vector<Index> h_lag_iRow_;
std::vector<Index> h_lag_jCol_; std::vector<Index> h_lag_jCol_;
// Sparse Matrix of problem (don't initialize every time!!!)
SMatrixNP HP_;
}; };
......
...@@ -11,8 +11,7 @@ ...@@ -11,8 +11,7 @@
//== INCLUDES ================================================================= //== INCLUDES =================================================================
#include <gmm/gmm.h> #include "NConstraintInterface.hh"
#include "NConstraintGmmInterface.hh"
//== FORWARDDECLARATIONS ====================================================== //== FORWARDDECLARATIONS ======================================================
...@@ -30,13 +29,9 @@ namespace COMISO { ...@@ -30,13 +29,9 @@ namespace COMISO {
A more elaborate description follows. A more elaborate description follows.
*/ */
class LinearConstraint : public NConstraintGmmInterface class LinearConstraint : public NConstraintInterface
{ {
public: public:
// ToDo: appropriate Vector/MatrixType ???
typedef gmm::wsvector<double> SVectorNP;
typedef gmm::row_matrix< SVectorNP > SMatrixNP;
// use c-arrays as vectors for gmm // use c-arrays as vectors for gmm
typedef gmm::array1D_reference<double*> VectorPT; typedef gmm::array1D_reference<double*> VectorPT;
...@@ -45,31 +40,13 @@ public: ...@@ -45,31 +40,13 @@ public:
// enum ConstraintType {NC_EQUAL, NC_LESS_EQUAL, NC_GREATER_EQUAL}; // enum ConstraintType {NC_EQUAL, NC_LESS_EQUAL, NC_GREATER_EQUAL};
/// Default constructor /// Default constructor
LinearConstraint(const ConstraintType _type = NC_EQUAL) : NConstraintGmmInterface(_type) LinearConstraint(const ConstraintType _type = NC_EQUAL) : NConstraintInterface(_type)
{} {}
// linear equation of the form -> coeffs_^T * (x,1) =_type= 0 // linear equation of the form -> coeffs_^T * (x,1) =_type= 0
LinearConstraint(const SVectorNP& _coeffs_plus_const, const ConstraintType _type = NC_EQUAL) : NConstraintGmmInterface(_type) LinearConstraint(const SVectorNC& _coeffs, const double _b, const ConstraintType _type = NC_EQUAL) : NConstraintInterface(_type)
{
int n = gmm::vect_size(_coeffs_plus_const)-1;
gmm::resize(coeffs_, n+1);
gmm::copy(_coeffs_plus_const, coeffs_);
if( n >= 0)
b_ = coeffs_[n];
else
b_ = 0.0;
gmm::resize(coeffs_, n);
}
// linear equation of the form -> coeffs_^T * (x,1) =_type= 0
LinearConstraint(const SVectorNP& _coeffs, const double _b, const ConstraintType _type = NC_EQUAL) : NConstraintGmmInterface(_type)
{ {
int n = gmm::vect_size(_coeffs); coeffs_ = _coeffs;
gmm::resize(coeffs_, n);
gmm::copy(_coeffs, coeffs_);
b_ = _b; b_ = _b;
} }
...@@ -78,24 +55,29 @@ public: ...@@ -78,24 +55,29 @@ public:
virtual int n_unknowns() virtual int n_unknowns()
{ {
return gmm::vect_size(coeffs_); return coeffs_.innerSize();
} }
virtual double eval_constraint ( const double* _x ) virtual double eval_constraint ( const double* _x )
{ {
return (gmm::vect_sp(coeffs_, VectorPT((double*)_x, gmm::vect_size(coeffs_))) + b_); double v = b_;
SVectorNC::InnerIterator c_it(coeffs_);
for(; c_it; ++c_it)
v += c_it.value()*_x[c_it.index()];
return v;
} }
virtual void eval_gradient( const double* _x, SVectorNP& _g ) virtual void eval_gradient( const double* _x, SVectorNC& _g )
{ {
gmm::resize(_g, gmm::vect_size(coeffs_)); _g = coeffs_;
gmm::copy (coeffs_, _g);
} }
virtual void eval_hessian ( const double* _x, SMatrixNP& _h ) virtual void eval_hessian ( const double* _x, SMatrixNC& _h )
{ {
gmm::resize(_h, gmm::vect_size(coeffs_), gmm::vect_size(coeffs_)); _h.clear();
gmm::clear(_h); _h.resize(coeffs_.innerSize(), coeffs_.innerSize());
} }
// inherited from base // inherited from base
...@@ -104,7 +86,7 @@ public: ...@@ -104,7 +86,7 @@ public:
private: private:
// linear equation of the form -> coeffs_^T * x + b_ // linear equation of the form -> coeffs_^T * x + b_
SVectorNP coeffs_; SVectorNC coeffs_;
double b_; double b_;
}; };
......
//=============================================================================
//
// CLASS NConstraintInterface
//
//=============================================================================
#ifndef COMISO_NCONSTRAINTINTERFACE_HH
#define COMISO_NCONSTRAINTINTERFACE_HH
//== INCLUDES =================================================================
#include "SuperSparseMatrixT.hh"
#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
#include <Eigen/Sparse>
//== FORWARDDECLARATIONS ======================================================
//== NAMESPACES ===============================================================
namespace COMISO {
//== CLASS DEFINITION =========================================================
/** \class NProblemGmmInterface NProblemGmmInterface.hh <ACG/.../NPRoblemGmmInterface.hh>
Brief Description.
A more elaborate description follows.
*/
class NConstraintInterface
{
public:
// define Sparse Datatypes
typedef Eigen::SparseVector<double> SVectorNC;
typedef SuperSparseMatrixT<double> SMatrixNC;
// different types of constraints
enum ConstraintType {NC_EQUAL, NC_LESS_EQUAL, NC_GREATER_EQUAL};
/// Default constructor
NConstraintInterface(const ConstraintType _type = NC_EQUAL) : type_(_type) {}
/// Destructor
~NConstraintInterface() {}
virtual int n_unknowns ( ) = 0;
virtual double eval_constraint ( const double* _x ) = 0;
virtual void eval_gradient ( const double* _x, SVectorNC& _g ) = 0;
virtual void eval_hessian ( const double* _x, SMatrixNC& _h ) = 0;
virtual ConstraintType constraint_type ( ) { return type_; }
virtual bool is_satisfied ( const double* _x, double _eps = 1e-6 )
{
switch( type_)
{
case NC_EQUAL : return (fabs(eval_constraint(_x)) <= _eps); break;
case NC_LESS_EQUAL : return ( eval_constraint(_x) <= _eps); break;
case NC_GREATER_EQUAL: return ( eval_constraint(_x) >= -_eps); break;
}
return false;
}
virtual double gradient_update_factor( const double* _x, double _eps = 1e-6 )
{
double val = eval_constraint(_x);
bool upper_bound_ok = ( val <= _eps);
bool lower_bound_ok = ( val >= -_eps);
if(upper_bound_ok)
{
if(lower_bound_ok || type_ == NC_LESS_EQUAL) return 0.0;
else return 1.0;
}
else
{
if(lower_bound_ok && type_ == NC_GREATER_EQUAL) return 0.0;
else return -1.0;
}
}
private:
// constraint type
ConstraintType type_;
};
//=============================================================================