1 /*===========================================================================*\ 2 * * 3 * CoMISo * 4 * Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen * 5 * www.rwth-graphics.de * 6 * * 7 *---------------------------------------------------------------------------* 8 * This file is part of CoMISo. * 9 * * 10 * CoMISo is free software: you can redistribute it and/or modify * 11 * it under the terms of the GNU General Public License as published by * 12 * the Free Software Foundation, either version 3 of the License, or * 13 * (at your option) any later version. * 14 * * 15 * CoMISo is distributed in the hope that it will be useful, * 16 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 18 * GNU General Public License for more details. * 19 * * 20 * You should have received a copy of the GNU General Public License * 21 * along with CoMISo. If not, see <http://www.gnu.org/licenses/>. * 22 * * 23 \*===========================================================================*/ 24 25 26 //============================================================================= 27 // 28 // CLASS MISolver 29 // 30 //============================================================================= 31 32 33 #ifndef COMISO_MISOLVER_HH 34 #define COMISO_MISOLVER_HH 35 36 37 //== INCLUDES ================================================================= 38 #include <CoMISo/Config/CoMISoDefines.hh> 39 #include <CoMISo/Config/config.hh> 40 41 #if COMISO_SUITESPARSE_AVAILABLE 42 #include "CholmodSolver.hh" 43 #elif COMISO_EIGEN3_AVAILABLE 44 #include "EigenLDLTSolver.hh" 45 #else 46 #print "Warning: MISolver requires Suitesparse or Eigen3 support" 47 #endif 48 49 #include "GMM_Tools.hh" 50 #include "IterativeSolverT.hh" 51 52 #include <vector> 53 54 #define ROUND_MI(x) ((x)<0?int((x)-0.5):int((x)+0.5)) 55 56 57 //== FORWARDDECLARATIONS ====================================================== 58 59 60 namespace COMISO { 61 class MISolverDialog; 62 } 63 64 //== NAMESPACES =============================================================== 65 66 namespace COMISO { 67 68 //== CLASS DEFINITION ========================================================= 69 70 71 72 /** \class MISolver MISolver.hh 73 74 Mixed-Integer Solver. 75 Approximates the solution of a (mixed-)integer problem 76 by iteratively computing a continuous(real) minimizer x followed by a 77 rounding of the one variable x_i which is subsequently eliminated from the 78 system, and the system is solved again ... 79 */ 80 81 class COMISODLLEXPORT MISolver 82 { 83 public: 84 85 // typedefs 86 typedef gmm::csc_matrix< double > CSCMatrix; 87 typedef std::vector< double > Vecd; 88 typedef std::vector< int > Veci; 89 typedef std::vector< unsigned int > Vecui; 90 91 // gmm Column and ColumnIterator of CSC matrix 92 typedef gmm::linalg_traits< CSCMatrix >::const_sub_col_type Col; 93 typedef gmm::linalg_traits< Col >::const_iterator ColIter; 94 95 96 /// default Constructor 97 MISolver(); 98 99 /// Destructor ~MISolver()100 ~MISolver() {} 101 102 103 /// Compute greedy approximation to a mixed integer problem. 104 /** @param _A symmetric positive semi-definite CSC matrix (Will be \b destroyed!) 105 * @param _x vector holding solution at the end 106 * @param _rhs right hand side of system Ax=rhs (Will be \b destroyed!) 107 * @param _to_round vector with variable indices to round to integers 108 * @param _fixed_order specifies if _to_round indices shall be rounded in the 109 * given order (\b true) or be greedily selected (\b false) 110 * */ 111 inline void solve( 112 CSCMatrix& _A, 113 Vecd& _x, 114 Vecd& _rhs, 115 Veci& _to_round, 116 bool _fixed_order = false ); 117 118 void resolve( 119 Vecd& _x, 120 Vecd& _rhs ); 121 122 /// Compute greedy approximation to a mixed integer problem. 123 /** @param _B mx(n+1) matrix with (still non-squared) equations of the energy, 124 * including the right hand side (Will be \b destroyed!) 125 * @param _x vector holding solution at the end 126 * @param _to_round vector with variable indices to round to integers 127 * @param _fixed_order specifies if _to_round indices shall be rounded in the 128 * given order (\b true) or be greedily selected (\b false) 129 * */ 130 template<class CMatrixT> 131 void solve( 132 CMatrixT& _B, 133 Vecd& _x, 134 Veci& _to_round, 135 bool _fixed_order = false ); 136 137 138 /// show Qt-Options-Dialog for setting algorithm parameters 139 /** Requires a Qt Application running and COMISO_GUI to be defined */ 140 void show_options_dialog(); 141 142 /** @name Get/Set functions for algorithm parameters 143 * Besides being used by the Qt-Dialog these can also be called explicitly 144 * to set parameters of the algorithm. */ 145 /*@{*/ 146 /// Shall an initial full solution be computed? set_inital_full(bool _b)147 void set_inital_full( bool _b) { initial_full_solution_=_b;} 148 /// Will an initial full solution be computed? get_inital_full()149 bool get_inital_full() { return initial_full_solution_;} 150 151 /// Shall an full solution be computed if iterative methods did not converged? set_iter_full(bool _b)152 void set_iter_full( bool _b) { iter_full_solution_=_b;} 153 /// Will an full solution be computed if iterative methods did not converged? get_iter_full()154 bool get_iter_full() { return iter_full_solution_;} 155 156 /// Shall a final full solution be computed? set_final_full(bool _b)157 void set_final_full( bool _b) { final_full_solution_=_b;} 158 /// Will a final full solution be computed? get_final_full()159 bool get_final_full() { return final_full_solution_;} 160 161 /// Shall direct (or greedy) rounding be used? set_direct_rounding(bool _b)162 void set_direct_rounding( bool _b) { direct_rounding_=_b;} 163 /// Will direct rounding be used? get_direct_rounding()164 bool get_direct_rounding() { return direct_rounding_;} 165 166 /// Shall no rounding be performed? set_no_rounding(bool _b)167 void set_no_rounding( bool _b) { no_rounding_=_b;} 168 /// Will no rounding be performed? get_no_rounding()169 bool get_no_rounding() { return no_rounding_;} 170 171 /// Shall multiple rounding be performed? set_multiple_rounding(bool _b)172 void set_multiple_rounding( bool _b) { multiple_rounding_=_b;} 173 /// Will multiple rounding be performed? get_multiple_rounding()174 bool get_multiple_rounding() { return multiple_rounding_;} 175 176 /// Shall gurobi solver be used? set_gurobi_rounding(bool _b)177 void set_gurobi_rounding( bool _b) { gurobi_rounding_=_b;} 178 /// Will gurobi rounding be performed? get_gurobi_rounding()179 bool get_gurobi_rounding() { return gurobi_rounding_;} 180 181 /// Shall cplex solver be used? set_cplex_rounding(bool _b)182 void set_cplex_rounding( bool _b) { cplex_rounding_=_b;} 183 /// Will cplex rounding be performed? get_cplex_rounding()184 bool get_cplex_rounding() { return cplex_rounding_;} 185 186 /// Set number of maximum Gauss-Seidel iterations set_local_iters(unsigned int _i)187 void set_local_iters( unsigned int _i) { max_local_iters_ = _i;} 188 /// Get number of maximum Gauss-Seidel iterations get_local_iters()189 unsigned int get_local_iters() { return max_local_iters_;} 190 191 /// Set error threshold for Gauss-Seidel solver set_local_error(double _d)192 void set_local_error( double _d) { max_local_error_ = _d;} 193 /// Get error threshold for Gauss-Seidel solver get_local_error()194 double get_local_error() { return max_local_error_;} 195 196 /// Set number of maximum Conjugate Gradient iterations set_cg_iters(unsigned int _i)197 void set_cg_iters( unsigned int _i) { max_cg_iters_ = _i;} 198 /// Get number of maximum Conjugate Gradient iterations get_cg_iters()199 unsigned int get_cg_iters() { return max_cg_iters_;} 200 201 /// Set error threshold for Conjugate Gradient set_cg_error(double _d)202 void set_cg_error( double _d) { max_cg_error_ = _d;} 203 /// Get error threshold for Conjugate Gradient get_cg_error()204 double get_cg_error() { return max_cg_error_;} 205 206 /// Set multiple rounding threshold (upper bound of rounding performed in each iteration) set_multiple_rounding_threshold(double _d)207 void set_multiple_rounding_threshold( double _d) { multiple_rounding_threshold_ = _d;} 208 /// Get multiple rounding threshold (upper bound of rounding performed in each iteration) get_multiple_rounding_threshold()209 double get_multiple_rounding_threshold() { return multiple_rounding_threshold_;} 210 211 /// Set noise level of algorithm. 0 - quiet, 1 - more noise, 2 - even more, 100 - all noise set_noise(unsigned int _i)212 void set_noise( unsigned int _i) { noisy_ = _i;} 213 /// Get noise level of algorithm get_noise()214 unsigned int get_noise() { return noisy_;} 215 216 /// Set time limit for gurobi solver (in seconds) set_gurobi_max_time(double _d)217 void set_gurobi_max_time( double _d) { gurobi_max_time_ = _d;} 218 /// Get time limit for gurobi solver (in seconds) get_gurobi_max_time()219 double get_gurobi_max_time() { return gurobi_max_time_;} 220 221 /// Set output statistics of solver set_stats(bool _stats)222 void set_stats( bool _stats) { stats_ = _stats; } 223 /// Get output statistics of solver get_stats()224 bool get_stats( ) { return stats_; } 225 /*@}*/ 226 227 /// Set/Get use_constraint_reordering for constraint solver (default = true) use_constraint_reordering()228 bool& use_constraint_reordering() { return use_constraint_reordering_;} 229 230 private: 231 232 // find set of variables for simultaneous rounding 233 class RoundingSet 234 { 235 public: 236 typedef std::pair<double,int> PairDI; 237 RoundingSet()238 RoundingSet() : threshold_(0.5), cur_sum_(0.0) {} 239 clear()240 void clear() { rset_.clear(); cur_sum_ = 0.0;} 241 add(int _id,double _rd_val)242 bool add( int _id, double _rd_val) 243 { 244 // empty set? -> always add one element 245 if( rset_.empty() || cur_sum_+_rd_val <= threshold_) 246 { 247 rset_.insert( PairDI(_rd_val,_id) ); 248 cur_sum_ += _rd_val; 249 return true; 250 } 251 else 252 { 253 // move to last element 254 std::set<PairDI>::iterator s_it = rset_.end(); 255 --s_it; 256 257 if( s_it->first > _rd_val) 258 { 259 cur_sum_ -= s_it->first; 260 rset_.erase(s_it); 261 rset_.insert( PairDI(_rd_val,_id) ); 262 cur_sum_ += _rd_val; 263 return true; 264 } 265 } 266 return false; 267 } 268 set_threshold(double _thres)269 void set_threshold( double _thres) { threshold_ = _thres; } 270 get_ids(std::vector<int> & _ids)271 void get_ids( std::vector<int>& _ids ) 272 { 273 _ids.clear(); 274 _ids.reserve( rset_.size()); 275 std::set<PairDI>::iterator s_it = rset_.begin(); 276 for(; s_it != rset_.end(); ++s_it) 277 _ids.push_back( s_it->second); 278 } 279 280 private: 281 282 double threshold_; 283 double cur_sum_; 284 285 std::set<PairDI> rset_; 286 287 std::set<PairDI> test_; 288 }; 289 290 private: 291 292 void solve_no_rounding( 293 CSCMatrix& _A, 294 Vecd& _x, 295 Vecd& _rhs ); 296 297 void solve_direct_rounding( 298 CSCMatrix& _A, 299 Vecd& _x, 300 Vecd& _rhs, 301 Veci& _to_round); 302 303 void solve_multiple_rounding( 304 CSCMatrix& _A, 305 Vecd& _x, 306 Vecd& _rhs, 307 Veci& _to_round ); 308 309 void solve_iterative( 310 CSCMatrix& _A, 311 Vecd& _x, 312 Vecd& _rhs, 313 Veci& _to_round, 314 bool _fixed_order ); 315 316 void solve_gurobi( 317 CSCMatrix& _A, 318 Vecd& _x, 319 Vecd& _rhs, 320 Veci& _to_round ); 321 322 inline void solve_cplex( 323 CSCMatrix& _A, 324 Vecd& _x, 325 Vecd& _rhs, 326 Veci& _to_round ); 327 328 329 void update_solution( 330 CSCMatrix& _A, 331 Vecd& _x, 332 Vecd& _rhs, 333 Vecui& _neigh_i ); 334 335 private: 336 337 /// Copy constructor (not used) 338 MISolver(const MISolver& _rhs); 339 340 /// Assignment operator (not used) 341 MISolver& operator=(const MISolver& _rhs); 342 343 // parameters used by the MiSo 344 bool initial_full_solution_; 345 bool iter_full_solution_; 346 bool final_full_solution_; 347 348 bool direct_rounding_; 349 bool no_rounding_; 350 bool multiple_rounding_; 351 bool gurobi_rounding_; 352 bool cplex_rounding_; 353 354 double multiple_rounding_threshold_; 355 356 unsigned int max_local_iters_; 357 double max_local_error_; 358 unsigned int max_cg_iters_; 359 double max_cg_error_; 360 double max_full_error_; 361 unsigned int noisy_; 362 bool stats_; 363 364 // time limit for gurobi solver (in seconds) 365 double gurobi_max_time_; 366 367 // flag 368 bool cholmod_step_done_; 369 370 // declar direct solver depending on availability 371 #if COMISO_SUITESPARSE_AVAILABLE 372 COMISO::CholmodSolver direct_solver_; 373 #elif COMISO_EIGEN3_AVAILABLE 374 COMISO::EigenLDLTSolver direct_solver_; 375 #else 376 #print "Warning: MISolver requires Suitesparse or Eigen3 support" 377 #endif 378 379 IterativeSolverT<double> siter_; 380 381 // statistics 382 unsigned int n_local_; 383 unsigned int n_cg_; 384 unsigned int n_full_; 385 386 bool use_constraint_reordering_; 387 388 #if(COMISO_QT4_AVAILABLE) 389 friend class COMISO::MISolverDialog; 390 #endif 391 }; 392 393 394 //============================================================================= 395 } // namespace COMISO 396 //============================================================================= 397 #if defined(INCLUDE_TEMPLATES) && !defined(COMISO_MISOLVER_C) 398 #define COMISO_MISOLVER_TEMPLATES 399 #include "MISolverT.cc" 400 #endif 401 //============================================================================= 402 #endif // COMISO_MISOLVER_HH defined 403 //============================================================================= 404 405