MISolver.hh 12.1 KB
Newer Older
Henrik Zimmer's avatar
Henrik Zimmer committed
1
2
/*===========================================================================*\
 *                                                                           *
Henrik Zimmer's avatar
Henrik Zimmer committed
3
 *                               CoMISo                                      *
Henrik Zimmer's avatar
Henrik Zimmer committed
4
 *      Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen      *
Henrik Zimmer's avatar
Henrik Zimmer committed
5
 *                           www.rwth-graphics.de                            *
Henrik Zimmer's avatar
Henrik Zimmer committed
6
 *                                                                           *
Henrik Zimmer's avatar
Henrik Zimmer committed
7
8
 *---------------------------------------------------------------------------* 
 *  This file is part of CoMISo.                                             *
Henrik Zimmer's avatar
Henrik Zimmer committed
9
 *                                                                           *
Henrik Zimmer's avatar
Henrik Zimmer committed
10
11
12
13
 *  CoMISo is free software: you can redistribute it and/or modify           *
 *  it under the terms of the GNU General Public License as published by     *
 *  the Free Software Foundation, either version 3 of the License, or        *
 *  (at your option) any later version.                                      *
Henrik Zimmer's avatar
Henrik Zimmer committed
14
 *                                                                           *
Henrik Zimmer's avatar
Henrik Zimmer committed
15
16
17
 *  CoMISo is distributed in the hope that it will be useful,                *
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of           *
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the            *
Henrik Zimmer's avatar
Henrik Zimmer committed
18
 *  GNU General Public License for more details.                             *
Henrik Zimmer's avatar
Henrik Zimmer committed
19
 *                                                                           *
Henrik Zimmer's avatar
Henrik Zimmer committed
20
21
 *  You should have received a copy of the GNU General Public License        *
 *  along with CoMISo.  If not, see <http://www.gnu.org/licenses/>.          *
Henrik Zimmer's avatar
Henrik Zimmer committed
22
23
24
 *                                                                           *
\*===========================================================================*/ 

Henrik Zimmer's avatar
Henrik Zimmer committed
25

26
27
28
29
30
31
32
//=============================================================================
//
//  CLASS MISolver
//
//=============================================================================


33
34
#ifndef COMISO_MISOLVER_HH
#define COMISO_MISOLVER_HH
35
36
37


//== INCLUDES =================================================================
Henrik Zimmer's avatar
Henrik Zimmer committed
38
#include <CoMISo/Config/CoMISoDefines.hh>
David Bommes's avatar
changed    
David Bommes committed
39
#include <CoMISo/Config/config.hh>
40
41
42

#include "GMM_Tools.hh"
#include "CholmodSolver.hh"
David Bommes's avatar
David Bommes committed
43
#include "IterativeSolverT.hh"
44
45
46
47
48
49
50
51
52

#include <vector>

#define ROUND_MI(x) ((x)<0?int((x)-0.5):int((x)+0.5))


//== FORWARDDECLARATIONS ======================================================


53
namespace COMISO {
54
55
56
57
58
class MISolverDialog;
}

//== NAMESPACES ===============================================================

59
namespace COMISO {
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138

//== CLASS DEFINITION =========================================================



/** \class MISolver MISolver.hh 

    Mixed-Integer Solver.
    Approximates the solution of a (mixed-)integer problem
    by iteratively computing a continuous(real) minimizer x followed by a
    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;

  // gmm Column and ColumnIterator of CSC matrix
  typedef gmm::linalg_traits< CSCMatrix >::const_sub_col_type Col;
  typedef gmm::linalg_traits< Col >::const_iterator           ColIter;


  /// default Constructor
  MISolver();

  /// 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 _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 );

  /// Compute greedy approximation to a mixed integer problem.
	/** @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 _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 );

  /// show Qt-Options-Dialog for setting algorithm parameters
  /** Requires a Qt Application running and COMISO_GUI to be defined */
  void show_options_dialog();

  /** @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_;}

David Bommes's avatar
David Bommes committed
139
140
141
142
143
	/// 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_;}

144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
	/// 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_;}

  /// 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_;}

  /// 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_;}

David Bommes's avatar
David Bommes committed
159
160
161
162
163
  /// 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_;}

David Bommes's avatar
David Bommes committed
164
165
166
167
168
169
  /// 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_;}


170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
  /// Set number of maximum Gauss-Seidel iterations
  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_;}

  /// Set error threshold for Gauss-Seidel solver
  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_;}

  /// 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;}
  /// Get error threshold for Conjugate Gradient
  double get_cg_error()           { return max_cg_error_;}

David Bommes's avatar
David Bommes committed
190
191
192
193
  /// 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_;}
194
195
196
197
198
199

  /// 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_;}

David Bommes's avatar
David Bommes committed
200
201
202
203
204
  /// Set time limit for gurobi solver (in seconds)
  void   set_gurobi_max_time( double _d) { gurobi_max_time_ = _d;}
  /// Get time limit for gurobi solver (in seconds)
  double get_gurobi_max_time()          { return gurobi_max_time_;}

205
206
207
208
209
210
  /// Set output statistics of solver
  void set_stats( bool _stats) { stats_ = _stats; }
  /// Get output statistics of solver
  bool get_stats( )            { return stats_; }
	/*@}*/

211
212
213
  /// Set/Get use_constraint_reordering for constraint solver (default = true)
  bool& use_constraint_reordering() { return use_constraint_reordering_;}

David Bommes's avatar
David Bommes committed
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
 private:

  // find set of variables for simultaneous rounding
  class RoundingSet
  {
  public:
    typedef std::pair<double,int> PairDI;
    
    RoundingSet() : threshold_(0.5), cur_sum_(0.0) {}

    void clear() { rset_.clear(); cur_sum_ = 0.0;}

    bool add( int _id, double _rd_val)
    {
      // empty set? -> always add one element
      if( rset_.size() == 0 || cur_sum_+_rd_val <= threshold_)
      {
	rset_.insert( PairDI(_rd_val,_id) );
	cur_sum_ += _rd_val;
	return true;
      }
      else
      {
	// move to last element
	std::set<PairDI>::iterator s_it = rset_.end();
	--s_it;

	if( s_it->first > _rd_val)
	{
	  cur_sum_ -= s_it->first;
	  rset_.erase(s_it);
	  rset_.insert( PairDI(_rd_val,_id) );
	  cur_sum_ += _rd_val;
	  return true;
	}
      }
      return false;
    }
    
    void set_threshold( double _thres) { threshold_ = _thres; }

    void get_ids( std::vector<int>& _ids )
    {
      _ids.clear();
      _ids.reserve( rset_.size());
      std::set<PairDI>::iterator s_it = rset_.begin();
      for(; s_it != rset_.end(); ++s_it)
	_ids.push_back( s_it->second);
    }

  private:

    double threshold_;
    double cur_sum_;

    std::set<PairDI> rset_;

    std::set<PairDI> test_;
  };

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_multiple_rounding( 
    CSCMatrix& _A, 
    Vecd&      _x, 
    Vecd&      _rhs, 
    Veci&      _to_round );

  void solve_iterative(
    CSCMatrix& _A, 
    Vecd&      _x, 
    Vecd&      _rhs, 
    Veci&      _to_round,
    bool       _fixed_order );

David Bommes's avatar
David Bommes committed
300
301
302
303
304
305
  void solve_gurobi(
    CSCMatrix& _A,
    Vecd&      _x,
    Vecd&      _rhs,
    Veci&      _to_round );

David Bommes's avatar
David Bommes committed
306
307
308
309
310
311
312

void update_solution( 
    CSCMatrix& _A, 
    Vecd&      _x, 
    Vecd&      _rhs, 
    Vecui&     _neigh_i );

313
314
315
316
317
318
319
320
321
322
private:

  /// Copy constructor (not used)
  MISolver(const MISolver& _rhs);

  /// Assignment operator (not used)
  MISolver& operator=(const MISolver& _rhs);

  // parameters used by the MiSo
  bool initial_full_solution_;
David Bommes's avatar
David Bommes committed
323
  bool iter_full_solution_;
324
325
326
327
  bool final_full_solution_;

  bool direct_rounding_;
  bool no_rounding_;
David Bommes's avatar
David Bommes committed
328
  bool multiple_rounding_;
David Bommes's avatar
David Bommes committed
329
  bool gurobi_rounding_;
David Bommes's avatar
David Bommes committed
330
331

  double multiple_rounding_threshold_;
332
333
334
335
336
337
338
339
340

  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_;

David Bommes's avatar
David Bommes committed
341
342
343
  // time limit for gurobi solver (in seconds)
  double       gurobi_max_time_;

David Bommes's avatar
David Bommes committed
344
345
346
347
348
349
350
351
352
353
354
355
  // flag
  bool         cholmod_step_done_;

  COMISO::CholmodSolver chol_;

  IterativeSolverT<double> siter_;

  // statistics
  unsigned int n_local_;
  unsigned int n_cg_;
  unsigned int n_full_;

356
357
  bool use_constraint_reordering_;

David Bommes's avatar
changed    
David Bommes committed
358
#if(COMISO_QT4_AVAILABLE)
359
  friend class COMISO::MISolverDialog;
David Bommes's avatar
changed    
David Bommes committed
360
#endif
361
362
363
364
};


//=============================================================================
365
} // namespace COMISO
366
//=============================================================================
367
368
#if defined(INCLUDE_TEMPLATES) && !defined(COMISO_MISOLVER_C)
#define COMISO_MISOLVER_TEMPLATES
369
370
371
#include "MISolverT.cc"
#endif
//=============================================================================
372
#endif // COMISO_MISOLVER_HH defined
373
374
//=============================================================================