Commit 19dbc365 authored by Martin Marinov's avatar Martin Marinov

merge Solver changes from VCI

parents a0e47ef4 d64f46f6
//=============================================================================
//
// CLASS FiniteElementLogBarrier
//
//=============================================================================
#ifndef COMISO_FINITEELEMENTLOGBARRIER_HH
#define COMISO_FINITEELEMENTLOGBARRIER_HH
//== INCLUDES =================================================================
#include <CoMISo/Config/CoMISoDefines.hh>
#include "NProblemInterface.hh"
//== FORWARDDECLARATIONS ======================================================
//== NAMESPACES ===============================================================
namespace COMISO {
//== CLASS DEFINITION =========================================================
/** \class FiniteElementLogBarrierLowerBound
Implements function of the type f(x) = -c1*log(x-c0) with constants (c0,c1)
A more elaborate description follows.
*/
class FiniteElementLogBarrierLowerBound
{
public:
// define dimensions
const static int NV = 1;
const static int NC = 2;
typedef Eigen::Matrix<size_t,NV,1> VecI;
typedef Eigen::Matrix<double,NV,1> VecV;
typedef Eigen::Matrix<double,NC,1> VecC;
typedef Eigen::Triplet<double> Triplet;
inline double eval_f (const VecV& _x, const VecC& _c) const
{
return -_c[1]*std::log(_x[0]-_c[0]);
}
inline void eval_gradient(const VecV& _x, const VecC& _c, VecV& _g) const
{
_g[0] = -_c[1]/(_x[0]-_c[0]);
}
inline void eval_hessian (const VecV& _x, const VecC& _c, std::vector<Triplet>& _triplets) const
{
_triplets.clear();
_triplets.push_back(Triplet(0,0,_c[1]/std::pow(_x[0]-_c[0],2)));
}
inline double max_feasible_step(const VecV& _x, const VecV& _v, const VecC& _c)
{
if(_v[0] >=0.0)
return DBL_MAX;
else
return 0.999*(_c[0]-_x[0])/_v[0];
}
};
/** \class FiniteElementLogBarrierUpperBound
Implements function of the type f(x) = -c1*log(c0-x) with constants (c0,c1)
A more elaborate description follows.
*/
class FiniteElementLogBarrierUpperBound
{
public:
// define dimensions
const static int NV = 1;
const static int NC = 2;
typedef Eigen::Matrix<size_t,NV,1> VecI;
typedef Eigen::Matrix<double,NV,1> VecV;
typedef Eigen::Matrix<double,NC,1> VecC;
typedef Eigen::Triplet<double> Triplet;
inline double eval_f (const VecV& _x, const VecC& _c) const
{
return -_c[1]*std::log(_c[0]-_x[0]);
}
inline void eval_gradient(const VecV& _x, const VecC& _c, VecV& _g) const
{
_g[0] = _c[1]/(_c[0]-_x[0]);
}
inline void eval_hessian (const VecV& _x, const VecC& _c, std::vector<Triplet>& _triplets) const
{
_triplets.clear();
_triplets.push_back(Triplet(0,0,_c[1]/std::pow(_c[0]-_x[0],2)));
}
inline double max_feasible_step(const VecV& _x, const VecV& _v, const VecC& _c)
{
if(_v[0] <=0.0)
return DBL_MAX;
else
return 0.999*(_c[0]-_x[0])/_v[0];
}
};
//=============================================================================
} // namespace COMISO
//=============================================================================
#endif // COMISO_FINITEELEMENTLOGBARRIER_HH defined
//=============================================================================
......@@ -83,8 +83,8 @@ static void process_gurobi_exception(const GRBException& _exc)
bool
GUROBISolver::
solve(NProblemInterface* _problem,
std::vector<NConstraintInterface*>& _constraints,
std::vector<PairIndexVtype>& _discrete_constraints,
const std::vector<NConstraintInterface *> &_constraints,
const std::vector<PairIndexVtype> &_discrete_constraints,
const double _time_limit)
{
DEB_enter_func;
......@@ -273,8 +273,8 @@ solve(NProblemInterface* _problem,
bool
GUROBISolver::
solve_two_phase(NProblemInterface* _problem, // problem instance
std::vector<NConstraintInterface*>& _constraints, // linear constraints
std::vector<PairIndexVtype>& _discrete_constraints, // discrete constraints
const std::vector<NConstraintInterface*>& _constraints, // linear constraints
const std::vector<PairIndexVtype>& _discrete_constraints, // discrete constraints
const double _time_limit0, // time limit phase 1 in seconds
const double _gap0, // MIP gap phase 1
const double _time_limit1, // time limit phase 2 in seconds
......@@ -288,8 +288,8 @@ solve_two_phase(NProblemInterface* _problem, //
bool
GUROBISolver::
solve_two_phase(NProblemInterface* _problem, // problem instance
std::vector<NConstraintInterface*>& _constraints, // linear constraints
std::vector<PairIndexVtype>& _discrete_constraints, // discrete constraints
const std::vector<NConstraintInterface *> &_constraints, // linear constraints
const std::vector<PairIndexVtype> &_discrete_constraints, // discrete constraints
const double _time_limit0, // time limit phase 1 in seconds
const double _gap0, // MIP gap phase 1
const double _time_limit1, // time limit phase 2 in seconds
......
......@@ -43,14 +43,14 @@ class COMISODLLEXPORT GUROBISolver
{
public:
// ********** 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
const std::vector<NConstraintInterface*>& _constraints, // linear constraints
const 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
const 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);
......@@ -62,16 +62,16 @@ public:
// 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
bool solve_two_phase(NProblemInterface* _problem, // problem instance
std::vector<NConstraintInterface*>& _constraints, // linear constraints
std::vector<PairIndexVtype>& _discrete_constraints, // discrete constraints
const std::vector<NConstraintInterface*>& _constraints, // linear constraints
const 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
bool solve_two_phase(NProblemInterface* _problem, // problem instance
std::vector<NConstraintInterface*>& _constraints, // linear constraints
std::vector<PairIndexVtype>& _discrete_constraints, // discrete constraints
const std::vector<NConstraintInterface*>& _constraints, // linear constraints
const std::vector<PairIndexVtype>& _discrete_constraints, // discrete constraints
const double _time_limit0, // time limit phase 1 in seconds
const double _gap0 , // MIP gap phase 1
const double _time_limit1, // time limit phase 2 in seconds
......
......@@ -101,6 +101,8 @@ public:
// 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");
Ipopt::IpoptApplication& app() {return (*app_); }
......
......@@ -41,6 +41,7 @@
#include <Base/Utils/StopWatch.hh>
#include <Base/Debug/DebOut.hh>
#include <type_traits>
//== NAMESPACES ===============================================================
......@@ -48,6 +49,11 @@ namespace COMISO {
//== IMPLEMENTATION ==========================================================
// cf. issue #3 - gmm5.2 compat
template <typename T>
using linalg_traits = typename gmm::linalg_traits<typename std::remove_const<typename std::remove_reference<T>::type>::type>;
template<class RMatrixT, class CMatrixT, class VectorT, class VectorIT >
void
ConstrainedSolver::solve_const(const RMatrixT& _constraints,
......@@ -286,10 +292,10 @@ ConstrainedSolver::resolve(
gmm::size_type m = gmm::mat_nrows(_B);
gmm::size_type n = gmm::mat_ncols(_B);
typedef typename gmm::linalg_traits<RMatrixT>::const_sub_row_type CRowT;
typedef typename gmm::linalg_traits<RMatrixT>::sub_row_type RowT;
typedef typename gmm::linalg_traits<CRowT>::const_iterator RIter;
typedef typename gmm::linalg_traits<CRowT>::value_type VecValT;
typedef typename linalg_traits<RMatrixT>::const_sub_row_type CRowT;
typedef typename linalg_traits<RMatrixT>::sub_row_type RowT;
typedef typename linalg_traits<CRowT>::const_iterator RIter;
typedef typename linalg_traits<CRowT>::value_type VecValT;
gmm::resize(rhs, n - 1);
gmm::clear(rhs);
......@@ -425,9 +431,9 @@ ConstrainedSolver::make_constraints_independent(
// if not found for integers with value +-1
// and finally take the smallest integer variable
typedef typename gmm::linalg_traits<RMatrixT>::const_sub_row_type CRowT;
typedef typename gmm::linalg_traits<RMatrixT>::sub_row_type RowT;
typedef typename gmm::linalg_traits<CRowT>::const_iterator RIter;
typedef typename linalg_traits<RMatrixT>::const_sub_row_type CRowT;
typedef typename linalg_traits<RMatrixT>::sub_row_type RowT;
typedef typename linalg_traits<CRowT>::const_iterator RIter;
// get current condition row
CRowT row = gmm::mat_const_row(_constraints, i);
......@@ -535,8 +541,8 @@ ConstrainedSolver::make_constraints_independent(
CVector col = constraints_c.col(elim_j);
// iterate over column
typename gmm::linalg_traits<CVector>::const_iterator c_it = gmm::vect_const_begin(col);
typename gmm::linalg_traits<CVector>::const_iterator c_end = gmm::vect_const_end(col);
typename linalg_traits<CVector>::const_iterator c_it = gmm::vect_const_begin(col);
typename linalg_traits<CVector>::const_iterator c_end = gmm::vect_const_end(col);
for (; c_it != c_end; ++c_it)
{
......@@ -635,9 +641,9 @@ ConstrainedSolver::make_constraints_independent_reordering(
// if not found for integers with value +-1
// and finally take the smallest integer variable
typedef typename gmm::linalg_traits<RMatrixT>::const_sub_row_type CRowT;
typedef typename gmm::linalg_traits<RMatrixT>::sub_row_type RowT;
typedef typename gmm::linalg_traits<CRowT>::const_iterator RIter;
typedef typename linalg_traits<RMatrixT>::const_sub_row_type CRowT;
typedef typename linalg_traits<RMatrixT>::sub_row_type RowT;
typedef typename linalg_traits<CRowT>::const_iterator RIter;
// get current condition row
CRowT row = gmm::mat_const_row(_constraints, i);
......@@ -744,8 +750,8 @@ ConstrainedSolver::make_constraints_independent_reordering(
CVector col = constraints_c.col(elim_j);
// iterate over column
typename gmm::linalg_traits<CVector>::const_iterator c_it = gmm::vect_const_begin(col);
typename gmm::linalg_traits<CVector>::const_iterator c_end = gmm::vect_const_end(col);
typename linalg_traits<CVector>::const_iterator c_it = gmm::vect_const_begin(col);
typename linalg_traits<CVector>::const_iterator c_end = gmm::vect_const_end(col);
for (; c_it != c_end; ++c_it)
{
......@@ -828,8 +834,8 @@ ConstrainedSolver::update_constraint_gcd(RMatrixT& _constraints,
return false;
// divide by gcd
typedef typename gmm::linalg_traits<RMatrixT>::const_sub_row_type CRowT;
typedef typename gmm::linalg_traits<CRowT>::const_iterator RIter;
typedef typename linalg_traits<RMatrixT>::const_sub_row_type CRowT;
typedef typename linalg_traits<CRowT>::const_iterator RIter;
// get current constraint row
RIter row_it = gmm::vect_const_begin(gmm::mat_const_row(_constraints, _row_i));
......@@ -882,8 +888,8 @@ ConstrainedSolver::eliminate_constraints(
SVector3T col = _Bcol.col(cur_j);
// iterate over column
typename gmm::linalg_traits<SVector3T>::const_iterator c_it = gmm::vect_const_begin(col);
typename gmm::linalg_traits<SVector3T>::const_iterator c_end = gmm::vect_const_end(col);
typename linalg_traits<SVector3T>::const_iterator c_it = gmm::vect_const_begin(col);
typename linalg_traits<SVector3T>::const_iterator c_end = gmm::vect_const_end(col);
for (; c_it != c_end; ++c_it)
add_row(c_it.index(), -(*c_it) / cur_val, gmm::mat_row(_constraints, i), _Bcol);
......@@ -946,8 +952,8 @@ ConstrainedSolver::eliminate_constraints(
Base::StopWatch sw;
sw.start();
// define iterator on matrix A and on constraints C
typedef typename gmm::linalg_traits<SVector2T>::const_iterator AIter;
typedef typename gmm::linalg_traits<SVector1T>::const_iterator CIter;
typedef typename linalg_traits<SVector2T>::const_iterator AIter;
typedef typename linalg_traits<SVector1T>::const_iterator CIter;
// store variable indices to be eliminated
std::vector<int> elim_varids;
......@@ -1091,7 +1097,7 @@ ConstrainedSolver::add_row(gmm::size_type _row_i,
RowT _row,
MatrixT& _mat)
{
typedef typename gmm::linalg_traits<RowT>::const_iterator RIter;
typedef typename linalg_traits<RowT>::const_iterator RIter;
RIter r_it = gmm::vect_const_begin(_row);
RIter r_end = gmm::vect_const_end(_row);
......@@ -1111,7 +1117,7 @@ ConstrainedSolver::add_row_simultaneously(gmm::size_type _row_i,
RMatrixT& _rmat,
CMatrixT& _cmat)
{
typedef typename gmm::linalg_traits<RowT>::const_iterator RIter;
typedef typename linalg_traits<RowT>::const_iterator RIter;
RIter r_it = gmm::vect_const_begin(_row);
RIter r_end = gmm::vect_const_end(_row);
......@@ -1325,8 +1331,8 @@ ConstrainedSolver::verify_constrained_system(
double _eps)
{
DEB_enter_func;
typedef typename gmm::linalg_traits<RMatrixT>::const_sub_row_type RowT;
typedef typename gmm::linalg_traits<RowT>::const_iterator RIter;
typedef typename linalg_traits<RMatrixT>::const_sub_row_type RowT;
typedef typename linalg_traits<RowT>::const_iterator RIter;
VectorT Ax(_x.size());
gmm::mult(_A, _x, Ax);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment