Commit 26871bc5 authored by Max Lyon's avatar Max Lyon
parents 1a3f03ac aae8fdd0
Pipeline #15445 passed with stages
in 4 minutes and 44 seconds
......@@ -258,7 +258,7 @@ public:
void set_noisy( int _noisy) { noisy_ = _noisy;}
// Get/Set whether the constraint reordering is used (default true)
bool& use_constraint_reordering() { return miso_.use_constraint_reordering(); }
bool use_constraint_reordering = true;
/// Access the MISolver (e.g. to change settings)
COMISO::MISolver& misolver() { return miso_;}
......
......@@ -167,7 +167,7 @@ ConstrainedSolver::solve(
DEB_out_if(_show_timings, 1, "Initital dimension: " << nrows << " x " << ncols
<< ", number of constraints: " << ncons << ", number of integer variables: " << _idx_to_round.size()
<< ", use reordering: " << (use_constraint_reordering() ? "yes" : "no") << "\n")
<< ", use reordering: " << (use_constraint_reordering ? "yes" : "no") << "\n")
// StopWatch for Timings
Base::StopWatch sw, sw2; sw.start(); sw2.start();
......@@ -177,7 +177,7 @@ ConstrainedSolver::solve(
std::vector<int> c_elim(ncons);
// apply sparse gauss elimination to make subsequent _conditions independent
if (use_constraint_reordering())
if (use_constraint_reordering)
make_constraints_independent_reordering(_constraints, _idx_to_round, c_elim);
else
make_constraints_independent(_constraints, _idx_to_round, c_elim);
......
......@@ -9,6 +9,7 @@
//== INCLUDES =================================================================
#include "IterativeSolverT.hh"
#include <Base/Debug/DebOut.hh>
//== NAMESPACES ===============================================================
......@@ -19,8 +20,8 @@ namespace COMISO
template <class RealT>
bool IterativeSolverT<RealT>::gauss_seidel_local(const Matrix& _A, Vector& _x,
const Vector& _rhs, const std::vector<unsigned int>& _idxs,
const int _max_iter, const Real& _tolerance)
const Vector& _rhs, const IndexVector& _idxs, const int _max_iter,
const Real& _tolerance)
{
if (_max_iter == 0)
return false;
......@@ -28,21 +29,20 @@ bool IterativeSolverT<RealT>::gauss_seidel_local(const Matrix& _A, Vector& _x,
typedef typename gmm::linalg_traits<Matrix>::const_sub_col_type ColT;
typedef typename gmm::linalg_traits<ColT>::const_iterator CIter;
// clear old data
i_temp.clear();
q.clear();
updt_vrbl_indcs_.clear();
indx_queue_.clear();
for (unsigned int i = 0; i < _idxs.size(); ++i)
q.push_back(_idxs[i]);
indx_queue_.push_back(_idxs[i]);
int it_count = 0;
while (!q.empty() && it_count < _max_iter)
while (!indx_queue_.empty() && it_count < _max_iter)
{
++it_count;
const auto i = q.front();
q.pop_front();
i_temp.clear();
const auto i = indx_queue_.front();
indx_queue_.pop_front();
indx_temp_.clear();
double res_i = -_rhs[i];
double x_i_new = _rhs[i];
......@@ -56,7 +56,7 @@ bool IterativeSolverT<RealT>::gauss_seidel_local(const Matrix& _A, Vector& _x,
res_i += (*it) * _x[j];
x_i_new -= (*it) * _x[j];
if (j != i)
i_temp.push_back(j);
indx_temp_.push_back(j);
else
diag = *it;
}
......@@ -68,82 +68,13 @@ bool IterativeSolverT<RealT>::gauss_seidel_local(const Matrix& _A, Vector& _x,
if (fabs(res_i * diag) > _tolerance)
{
_x[i] += x_i_new * diag;
for (unsigned int j = 0; j < i_temp.size(); ++j)
q.push_back(i_temp[j]);
updt_vrbl_indcs_.push_back(i);
for (unsigned int j = 0; j < indx_temp_.size(); ++j)
indx_queue_.push_back(indx_temp_[j]);
}
}
return q.empty(); // converged?
}
//-----------------------------------------------------------------------------
template <class RealT>
bool IterativeSolverT<RealT>::gauss_seidel_local2(const Matrix& _A, Vector& _x,
const Vector& _rhs, const std::vector<unsigned int>& _idxs,
const int _max_iter, const Real& _tolerance)
{
typedef typename gmm::linalg_traits<Matrix>::const_sub_col_type ColT;
typedef typename gmm::linalg_traits<ColT>::const_iterator CIter;
double t2 = _tolerance * _tolerance;
// clear old data
i_temp.clear();
s.clear();
for (unsigned int i = 0; i < _idxs.size(); ++i)
s.insert(_idxs[i]);
int it_count = 0;
bool finished = false;
while (!finished && it_count < _max_iter)
{
finished = true;
std::set<int>::iterator s_it = s.begin();
for (; s_it != s.end(); ++s_it)
{
++it_count;
unsigned int cur_i = *s_it;
i_temp.clear();
ColT col = mat_const_col(_A, cur_i);
CIter it = gmm::vect_const_begin(col);
CIter ite = gmm::vect_const_end(col);
double res_i = -_rhs[cur_i];
double x_i_new = _rhs[cur_i];
double diag = 1.0;
for (; it != ite; ++it)
{
res_i += (*it) * _x[it.index()];
x_i_new -= (*it) * _x[it.index()];
if (it.index() != cur_i)
i_temp.push_back(static_cast<int>(it.index()));
else
diag = *it;
}
// compare relative residuum normalized by diagonal entry
if (res_i * res_i / diag > t2)
{
_x[cur_i] += x_i_new / _A(cur_i, cur_i);
for (unsigned int j = 0; j < i_temp.size(); ++j)
{
finished = false;
s.insert(i_temp[j]);
}
}
}
}
// converged?
return finished;
return indx_queue_.empty(); // converged?
}
//-----------------------------------------------------------------------------
......@@ -152,6 +83,7 @@ template <class RealT>
bool IterativeSolverT<RealT>::conjugate_gradient(const Matrix& _A, Vector& _x,
const Vector& _rhs, int& _max_iter, Real& _tolerance)
{
DEB_enter_func;
Real rho, rho_1(0), a;
// initialize vectors
......@@ -179,7 +111,7 @@ bool IterativeSolverT<RealT>::conjugate_gradient(const Matrix& _A, Vector& _x,
while ((not_converged = ((res_norm = vect_norm_rel(r_, d_)) > _tolerance)) &&
cur_iter < _max_iter)
{
// std::cerr << "iter " << cur_iter << " res " << res_norm << std::endl;
DEB_line(11, "iter " << cur_iter << " res " << res_norm);
if (cur_iter != 0)
{
......@@ -200,7 +132,7 @@ bool IterativeSolverT<RealT>::conjugate_gradient(const Matrix& _A, Vector& _x,
_max_iter = cur_iter;
_tolerance = res_norm;
return (!not_converged);
return !not_converged;
}
//-----------------------------------------------------------------------------
......@@ -210,19 +142,11 @@ typename IterativeSolverT<RealT>::Real IterativeSolverT<RealT>::vect_norm_rel(
const Vector& _v, const Vector& _diag) const
{
Real res = 0.0;
for (unsigned int i = 0; i < _v.size(); ++i)
{
res = std::max(fabs(_v[i] * _diag[i]), res);
// Real cur = fabs(_v[i]*_diag[i]);
// if(cur > res)
// res = cur;
}
return res;
}
//-----------------------------------------------------------------------------
//=============================================================================
} // namespace COMISO
......
......@@ -11,7 +11,6 @@
#include <CoMISo/Utils/gmm.hh>
#include <deque>
#include <set>
//== FORWARDDECLARATIONS ======================================================
......@@ -32,19 +31,21 @@ namespace COMISO
template <class RealT> class IterativeSolverT
{
public:
typedef unsigned int uint;
typedef RealT Real;
typedef std::vector<Real> Vector;
typedef std::vector<uint> IndexVector;
typedef gmm::csc_matrix<Real> Matrix;
// local gauss_seidel
// local Gauss-Seidel
bool gauss_seidel_local(const Matrix& _A, Vector& _x, const Vector& _rhs,
const std::vector<unsigned int>& _idxs, const int _max_iter,
const Real& _tolerance);
const IndexVector& _idxs, const int _max_iter, const Real& _tolerance);
// local gauss_seidel
bool gauss_seidel_local2(const Matrix& _A, Vector& _x, const Vector& _rhs,
const std::vector<unsigned int>& _idxs, const int _max_iter,
const Real& _tolerance);
// get the indices of any variables updated during the last local Gauss-Seidel
const IndexVector& updated_variable_indices() const
{
return updt_vrbl_indcs_;
}
// conjugate gradient
bool conjugate_gradient(const Matrix& _A, Vector& _x, const Vector& _rhs,
......@@ -55,16 +56,16 @@ private:
Real vect_norm_rel(const Vector& _v, const Vector& _diag) const;
private:
// helper for conjugate gradient
// context for Conjugate Gradient
Vector p_;
Vector q_;
Vector r_;
Vector d_;
// helper for local gauss seidel
std::vector<unsigned int> i_temp;
std::deque<unsigned int> q;
std::set<int> s;
// context for local Gauss-Seidel
IndexVector indx_temp_;
std::deque<uint> indx_queue_;
IndexVector updt_vrbl_indcs_; // updated variable indices
};
//=============================================================================
......
This diff is collapsed.
......@@ -37,16 +37,7 @@
#include <CoMISo/Config/CoMISoDefines.hh>
#include <CoMISo/Config/config.hh>
#if COMISO_SUITESPARSE_AVAILABLE
#include "CholmodSolver.hh"
#elif COMISO_EIGEN3_AVAILABLE
#include "EigenLDLTSolver.hh"
#else
#print "Warning: MISolver requires Suitesparse or Eigen3 support"
#endif
#include "GMM_Tools.hh"
#include "IterativeSolverT.hh"
#include <vector>
......@@ -59,8 +50,6 @@ class MISolverDialog;
//== CLASS DEFINITION =========================================================
/** \class MISolver MISolver.hh
Mixed-Integer Solver.
......@@ -69,247 +58,199 @@ class MISolverDialog;
rounding of the one variable x_i which is subsequently eliminated from the
system, and the system is solved again ...
*/
class COMISODLLEXPORT MISolver
{
public:
// typedefs
typedef gmm::csc_matrix< double > CSCMatrix;
typedef std::vector< double > Vecd;
typedef std::vector< int > Veci;
typedef std::vector< unsigned int > Vecui;
typedef gmm::csc_matrix<double> CSCMatrix;
typedef std::vector<double> Vecd;
typedef std::vector<int> Veci;
typedef std::vector<unsigned int> Vecui;
enum class RoundingType
{
NONE, // no rounding at all, use with caution
DIRECT, /* simple method that round all integer variables in one pass,
fast, but solutions are far from optimal. */
MULTIPLE, /* greedy method that rounds several variables at once,
controlled by set_multiple_rounding_threshold(). */
GUROBI, // only available if you have the Gurobi solver configured
CPLEX, // only available if you have the CPLEX solver configured
DEFAULT = MULTIPLE // default setting
};
/// default Constructor
MISolver();
/// delete copy constructor
MISolver(const MISolver&) = delete;
/// delete assignment operator
MISolver& operator=(const MISolver& _rhs) = delete;
// destructor
~MISolver();
/// Compute greedy approximation to a mixed integer problem.
/** @param _A symmetric positive semi-definite CSC matrix (Will be \b destroyed!)
* @param _x vector holding solution at the end
/** @param _A symmetric positive semi-definite CSC matrix (Will be \b
* destroyed!)
* @param _x vector holding solution at the end
* @param _rhs right hand side of system Ax=rhs (Will be \b destroyed!)
* @param _to_round vector with variable indices to round to integers
* @param _fixed_order specifies if _to_round indices shall be rounded in the
* given order (\b true) or be greedily selected (\b false)
* */
void solve(
CSCMatrix& _A,
Vecd& _x,
Vecd& _rhs,
Veci& _to_round,
bool _fixed_order = false );
void resolve(
Vecd& _x,
Vecd& _rhs );
* */
void solve(CSCMatrix& _A, Vecd& _x, Vecd& _rhs, const Veci& _to_round);
//! Resolve usig the direct solver
void resolve(Vecd& _x, Vecd& _rhs);
/// Compute greedy approximation to a mixed integer problem.
/** @param _B mx(n+1) matrix with (still non-squared) equations of the energy,
/** @param _B mx(n+1) matrix with (still non-squared) equations of the energy,
* including the right hand side (Will be \b destroyed!)
* @param _x vector holding solution at the end
* @param _x vector holding solution at the end
* @param _to_round vector with variable indices to round to integers
* @param _fixed_order specifies if _to_round indices shall be rounded in the
* given order (\b true) or be greedily selected (\b false)
* */
//template<class CMatrixT>
//void solve(
// CMatrixT& _B,
// Vecd& _x,
// Veci& _to_round,
// bool _fixed_order = false );
* */
// TODO: Missing function??
/// show Qt-Options-Dialog for setting algorithm parameters
/** Requires a Qt Application running and COMISO_GUI to be defined */
void show_options_dialog();
void show_options_dialog() const;
/** @name Get/Set functions for algorithm parameters
* Besides being used by the Qt-Dialog these can also be called explicitly
* to set parameters of the algorithm. */
/*@{*/
/// Shall an initial full solution be computed?
void set_inital_full( bool _b) { initial_full_solution_=_b;}
/// Will an initial full solution be computed?
bool get_inital_full() { return initial_full_solution_;}
/// Shall an full solution be computed if iterative methods did not converged?
void set_iter_full( bool _b) { iter_full_solution_=_b;}
/// Will an full solution be computed if iterative methods did not converged?
bool get_iter_full() { return iter_full_solution_;}
/// Shall a final full solution be computed?
void set_final_full( bool _b) { final_full_solution_=_b;}
/// Will a final full solution be computed?
bool get_final_full() { return final_full_solution_;}
/// Set the solve type
void set_rounding_type(const RoundingType _rt) { rounding_type_ = _rt; }
/// Shall direct (or greedy) rounding be used?
void set_direct_rounding( bool _b) { direct_rounding_=_b;}
/// Will direct rounding be used?
bool get_direct_rounding() { return direct_rounding_;}
/// Get the solve type
RoundingType get_rounding_type() const { return rounding_type_; }
/// Shall no rounding be performed?
void set_no_rounding( bool _b) { no_rounding_=_b;}
/// Will no rounding be performed?
bool get_no_rounding() { return no_rounding_;}
void set_no_rounding() { rounding_type_ = RoundingType::NONE; }
/// Shall direct (or greedy) rounding be used?
void set_direct_rounding() { rounding_type_ = RoundingType::DIRECT; }
/// Shall multiple rounding be performed?
void set_multiple_rounding( bool _b) { multiple_rounding_=_b;}
/// Will multiple rounding be performed?
bool get_multiple_rounding() { return multiple_rounding_;}
void set_multiple_rounding() { rounding_type_ = RoundingType::MULTIPLE; }
/// Shall gurobi solver be used?
void set_gurobi_rounding( bool _b) { gurobi_rounding_=_b;}
/// Will gurobi rounding be performed?
bool get_gurobi_rounding() { return gurobi_rounding_;}
/// Shall Gurobi solver be used?
void set_gurobi_rounding() { rounding_type_ = RoundingType::GUROBI; }
/// Shall cplex solver be used?
void set_cplex_rounding( bool _b) { cplex_rounding_=_b;}
/// Will cplex rounding be performed?
bool get_cplex_rounding() { return cplex_rounding_;}
/// Shall CPLEX solver be used?
void set_cplex_rounding() { rounding_type_ = RoundingType::CPLEX; }
/** @name Get/Set functions for algorithm parameters
* Besides being used by the Qt-Dialog these can also be called explicitly
* to set parameters of the algorithm. */
/*@{*/
/// Shall an initial full solution be computed?
void set_inital_full(bool _b) { initial_full_solution_ = _b; }
/// Will an initial full solution be computed?
bool get_inital_full() const { return initial_full_solution_; }
/// Shall an full solution be computed if iterative methods did not converged?
void set_iter_full(bool _b) { iter_full_solution_ = _b; }
/// Will an full solution be computed if iterative methods did not converged?
bool get_iter_full() const { return iter_full_solution_; }
/// Shall a final full solution be computed?
void set_final_full(bool _b) { final_full_solution_ = _b; }
/// Will a final full solution be computed?
bool get_final_full() const { return final_full_solution_; }
/// Set number of maximum Gauss-Seidel iterations
void set_local_iters( unsigned int _i) { max_local_iters_ = _i;}
void set_local_iters(unsigned int _i) { max_local_iters_ = _i; }
/// Get number of maximum Gauss-Seidel iterations
unsigned int get_local_iters() { return max_local_iters_;}
unsigned int get_local_iters() { return max_local_iters_; }
/// Set error threshold for Gauss-Seidel solver
void set_local_error( double _d) { max_local_error_ = _d;}
void set_local_error(double _d) { max_local_error_ = _d; }
/// Get error threshold for Gauss-Seidel solver
double get_local_error() { return max_local_error_;}
double get_local_error() { return max_local_error_; }
/// Set number of maximum Conjugate Gradient iterations
void set_cg_iters( unsigned int _i) { max_cg_iters_ = _i;}
/// Get number of maximum Conjugate Gradient iterations
unsigned int get_cg_iters() { return max_cg_iters_;}
/// Set number of maximum Conjugate Gradient iterations
void set_cg_iters(unsigned int _i) { max_cg_iters_ = _i; }
/// Get number of maximum Conjugate Gradient iterations
unsigned int get_cg_iters() { return max_cg_iters_; }
/// Set error threshold for Conjugate Gradient
void set_cg_error( double _d) { max_cg_error_ = _d;}
void set_cg_error(double _d) { max_cg_error_ = _d; }
/// Get error threshold for Conjugate Gradient
double get_cg_error() { return max_cg_error_;}
/// Set multiple rounding threshold (upper bound of rounding performed in each iteration)
void set_multiple_rounding_threshold( double _d) { multiple_rounding_threshold_ = _d;}
/// Get multiple rounding threshold (upper bound of rounding performed in each iteration)
double get_multiple_rounding_threshold() { return multiple_rounding_threshold_;}
/// Set noise level of algorithm. 0 - quiet, 1 - more noise, 2 - even more, 100 - all noise
void set_noise( unsigned int _i) { noisy_ = _i;}
/// Get noise level of algorithm
unsigned int get_noise() { return noisy_;}
double get_cg_error() { return max_cg_error_; }
/// Set multiple rounding threshold (upper bound of rounding performed in each
/// iteration)
void set_multiple_rounding_threshold(double _d)
{
multiple_rounding_threshold_ = _d;
}
/// Get multiple rounding threshold (upper bound of rounding performed in
/// each iteration)
double get_multiple_rounding_threshold()
{
return multiple_rounding_threshold_;
}
/// Set time limit for gurobi solver (in seconds)
void set_gurobi_max_time( double _d) { gurobi_max_time_ = _d;}
void set_gurobi_max_time(double _d) { max_time_ = _d; }
/// Get time limit for gurobi solver (in seconds)
double get_gurobi_max_time() { return gurobi_max_time_;}
/// Set output statistics of solver
void set_stats( bool _stats) { stats_ = _stats; }
/// Get output statistics of solver
bool get_stats( ) { return stats_; }
/*@}*/
/// Set/Get use_constraint_reordering for constraint solver (default = true)
bool& use_constraint_reordering() { return use_constraint_reordering_;}
double get_gurobi_max_time() { return max_time_; }
/*@}*/
private:
void solve_no_rounding(
CSCMatrix& _A,
Vecd& _x,
Vecd& _rhs );
void solve_direct_rounding(
CSCMatrix& _A,
Vecd& _x,
Vecd& _rhs,
Veci& _to_round);
void solve_no_rounding(CSCMatrix& _A, Vecd& _x, Vecd& _rhs);
void solve_direct_rounding(
CSCMatrix& _A, Vecd& _x, Vecd& _rhs, const Veci& _to_round);
void solve_multiple_rounding(
CSCMatrix& _A, Vecd& _x, Vecd& _rhs, const Veci& _to_round);
void solve_iterative(
CSCMatrix& _A,
Vecd& _x,
Vecd& _rhs,
Veci& _to_round,
bool _fixed_order );
void solve_gurobi(
CSCMatrix& _A,
Vecd& _x,
Vecd& _rhs,
Veci& _to_round );
inline void solve_cplex(
CSCMatrix& _A,
Vecd& _x,
Vecd& _rhs,
Veci& _to_round );
void update_solution(
CSCMatrix& _A, Vecd& _x, Vecd& _rhs, Veci& _to_round, bool _fixed_order);
void solve_gurobi(CSCMatrix& _A, Vecd& _x, Vecd& _rhs, const Veci& _to_round);
void solve_cplex(CSCMatrix& _A, Vecd& _x, Vecd& _rhs, const Veci& _to_round);
// return true if the solution has been improved only by local iterations
bool update_solution_is_local(
const CSCMatrix& _A, Vecd& _x, const Vecd& _rhs, const Vecui& _neigh_i);
private:
/// Copy constructor (not used)
MISolver(const MISolver& _rhs);
/// Assignment operator (not used)
MISolver& operator=(const MISolver& _rhs);
RoundingType rounding_type_;
// parameters used by the MiSo
bool initial_full_solution_;
bool iter_full_solution_;
bool final_full_solution_;
bool direct_rounding_;
bool no_rounding_;
bool multiple_rounding_;
bool gurobi_rounding_;
bool cplex_rounding_;
double multiple_rounding_threshold_;
unsigned int max_local_iters_;
double max_local_error_;
unsigned int max_cg_iters_;
double max_cg_error_;
double max_full_error_;
unsigned int noisy_;
bool stats_;
// time limit for Gurobi solver (in seconds)
double gurobi_max_time_;
// flag
bool cholmod_step_done_;
// declare direct solver depending on availability
#if COMISO_SUITESPARSE_AVAILABLE
COMISO::CholmodSolver direct_solver_;
#elif COMISO_EIGEN3_AVAILABLE
COMISO::EigenLDLTSolver direct_solver_;
#else
#print "Warning: MISolver requires Suitesparse or Eigen3 support"
#endif
IterativeSolverT<double> siter_;
double multiple_rounding_threshold_; // control solve_multiple_rounding()
double max_time_; // control the time limit for Gurobi and CPLEX (in seconds)
// the actual solver declarations are hidden in the implementation code
class DirectSolver;
class IterativeSolver;
DirectSolver* direct_solver_;