MISolver.hh 11.2 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>
39
40
41

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

#include <vector>

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


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


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

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

58
namespace COMISO {
59
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

//== 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
138
139
140
141
142
	/// 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_;}

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

163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
  /// 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
183
184
185
186
  /// 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_;}
187
188
189
190
191
192
193
194
195
196
197
198

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

  /// Set output statistics of solver
  void set_stats( bool _stats) { stats_ = _stats; }
  /// Get output statistics of solver
  bool get_stats( )            { return stats_; }
	/*@}*/

David Bommes's avatar
David Bommes committed
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
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
 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 );


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

292
293
294
295
296
297
298
299
300
301
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
302
  bool iter_full_solution_;
303
304
305
306
  bool final_full_solution_;

  bool direct_rounding_;
  bool no_rounding_;
David Bommes's avatar
David Bommes committed
307
308
309
  bool multiple_rounding_;

  double multiple_rounding_threshold_;
310
311
312
313
314
315
316
317
318

  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
319
320
321
322
323
324
325
326
327
328
329
330
  // flag
  bool         cholmod_step_done_;

  COMISO::CholmodSolver chol_;

  IterativeSolverT<double> siter_;

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

331
  friend class COMISO::MISolverDialog;
332
333
334
335
};


//=============================================================================
336
} // namespace COMISO
337
//=============================================================================
338
339
#if defined(INCLUDE_TEMPLATES) && !defined(COMISO_MISOLVER_C)
#define COMISO_MISOLVER_TEMPLATES
340
341
342
#include "MISolverT.cc"
#endif
//=============================================================================
343
#endif // COMISO_MISOLVER_HH defined
344
345
//=============================================================================