MISolver.cc 23.9 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
 *                                                                           *
7
 *---------------------------------------------------------------------------*
Henrik Zimmer's avatar
Henrik Zimmer committed
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
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
\*===========================================================================*/
Henrik Zimmer's avatar
Henrik Zimmer committed
24

David Bommes's avatar
David Bommes committed
25
#include <CoMISo/Config/config.hh>
26
#include "MISolver.hh"
27
#include "CoMISo/Utils/CoMISoError.hh"
28

29
#if (COMISO_QT_AVAILABLE)
30
#include <CoMISo/QtWidgets/MISolverDialogUI.hh>
31
32
#endif

33
34
35
36
37
#if COMISO_CPLEX_AVAILABLE
#include <ilcplex/ilocplex.h>
ILOSTLBEGIN
#endif

David Bommes's avatar
David Bommes committed
38
#if COMISO_GUROBI_AVAILABLE
39
#include <gurobi_c++.h>
David Bommes's avatar
David Bommes committed
40
41
#endif

42
//#define COMISO_MISOLVER_PERFORMANCE_TEST
43
44
45
46
47
48
#ifdef COMISO_MISOLVER_PERFORMANCE_TEST
#include "SparseQRSolver.hh"
#include "UMFPACKSolver.hh"
#include "EigenLDLTSolver.hh"
#endif

49
50
51
#include <CoMISo/Utils/gmm.hh>

#include <Base/Debug/DebTime.hh>
52
#include <Base/Utils/StopWatch.hh>
David Bommes's avatar
David Bommes committed
53

54
#include <queue>
55
56
#include <float.h>

57
#define ROUND(x) ((x) < 0 ? int((x) - 0.5) : int((x) + 0.5))
David Bommes's avatar
David Bommes committed
58

59
60
61
62
namespace COMISO
{
namespace
{
63

64
65
66
67
68
69
70
class RoundingQueue : protected std::priority_queue<std::pair<double, int>>
{
public:
  RoundingQueue(const double _threshold) : threshold_(_threshold), cur_sum_(0.0)
  {
    c.reserve(128);
  }
71

72
73
74
75
76
  void clear()
  {
    c.clear();
    cur_sum_ = 0.0;
  }
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
  void add(const int _id, const double _rd_val)
  {
    if (empty() || cur_sum_ + _rd_val <= threshold_)
    { // empty? -> always add one element
      emplace(_rd_val, _id);
      cur_sum_ += _rd_val;
      return;
    }
    if (top().first > _rd_val) // replace the last element if worse
    {
      cur_sum_ -= top().first;
      pop();
      emplace(_rd_val, _id);
      cur_sum_ += _rd_val;
    }
  }

  void get_ids(std::vector<int>& _ids)
  {
    _ids.clear();
    _ids.reserve(size());
    std::sort_heap(c.begin(), c.end());
    for (auto s_it = c.begin(); s_it != c.end(); ++s_it)
      _ids.push_back(s_it->second);
  }

private:
  const double threshold_;
  double cur_sum_;
};

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

} // namespace
114
115

// Constructor
116
MISolver::MISolver()
117
118
119
{
  // default parameters
  initial_full_solution_ = true;
120
121
  iter_full_solution_ = true;
  final_full_solution_ = true;
122

123
124
  direct_rounding_ = false;
  no_rounding_ = false;
David Bommes's avatar
David Bommes committed
125
  multiple_rounding_ = true;
126
127
  gurobi_rounding_ = false;
  cplex_rounding_ = false;
128

David Bommes's avatar
David Bommes committed
129
  max_local_iters_ = 100000;
130
  max_local_error_ = 1e-3;
131
132
  max_cg_iters_ = 50;
  max_cg_error_ = 1e-3;
David Bommes's avatar
David Bommes committed
133
134
135

  multiple_rounding_threshold_ = 0.5;

David Bommes's avatar
David Bommes committed
136
137
  gurobi_max_time_ = 60;

138
139
  noisy_ = 0;
  stats_ = true;
140
141

  use_constraint_reordering_ = true;
142
143
}

144
145
//-----------------------------------------------------------------------------

146
147
void MISolver::solve(
    CSCMatrix& _A, Vecd& _x, Vecd& _rhs, Veci& _to_round, bool _fixed_order)
148
{
149
  DEB_enter_func;
150

151
152
153
  DEB_out(6, "# integer    variables: "
                 << _to_round.size() << "\n# continuous variables: "
                 << _x.size() - _to_round.size() << "\n");
154

155
  // nothing to solve?
156
157
  if (gmm::mat_ncols(_A) == 0 || gmm::mat_nrows(_A) == 0)
    return;
158

159
  if (gurobi_rounding_)
160
    solve_gurobi(_A, _x, _rhs, _to_round);
161
  else if (cplex_rounding_)
162
    solve_cplex(_A, _x, _rhs, _to_round);
163
164
165
166
167
168
  else if (no_rounding_ || _to_round.size() == 0)
    solve_no_rounding(_A, _x, _rhs);
  else if (direct_rounding_)
    solve_direct_rounding(_A, _x, _rhs, _to_round);
  else if (multiple_rounding_)
    solve_multiple_rounding(_A, _x, _rhs, _to_round);
169
  else
170
    solve_iterative(_A, _x, _rhs, _to_round, _fixed_order);
171
172
173
174
175
}

//-----------------------------------------------------------------------------

// inline function...
176
void MISolver::solve_cplex(CSCMatrix& _A, Vecd& _x, Vecd& _rhs, Veci& _to_round)
177
178
{
  DEB_enter_func;
179
  DEB_out(2, "gurobi_max_time_: " << gurobi_max_time_ << "\n");
180
181
182
183
184

#if COMISO_CPLEX_AVAILABLE

  // get round-indices in set
  std::set<int> to_round;
185
  for (unsigned int i = 0; i < _to_round.size(); ++i)
186
187
    to_round.insert(_to_round[i]);

188
189
  try
  {
190
191
192
193
194
195
196
197
198
199
200

    IloEnv env_;
    IloModel model(env_);

    // set time limite
    //    model.getEnv().set(GRB_DoubleParam_TimeLimit, gurobi_max_time_);

    unsigned int n = _rhs.size();

    // 1. allocate variables
    std::vector<IloNumVar> vars;
201
202
203
    for (unsigned int i = 0; i < n; ++i)
      if (to_round.count(i))
        vars.push_back(IloNumVar(env_, -IloIntMax, IloIntMax, IloNumVar::Int));
204
      else
205
206
        vars.push_back(
            IloNumVar(env_, -IloInfinity, IloInfinity, IloNumVar::Float));
207
208
209
210
211
212

    // 2. setup_energy

    // build objective function from linear system E = x^tAx - 2x^t*rhs
    IloExpr objective(env_);

213
214
    for (unsigned int i = 0; i < _A.nc; ++i)
      for (unsigned int j = _A.jc[i]; j < _A.jc[i + 1]; ++j)
215
      {
216
        objective += _A.pr[j] * vars[_A.ir[j]] * vars[i];
217
      }
218
219
    for (unsigned int i = 0; i < n; ++i)
      objective -= 2 * _rhs[i] * vars[i];
220
221
222
223
224
225
226
227
228
229

    // ToDo: objective correction!!!

    //    _A.jc[c+1]
    //    _A.pr[write]
    //    _A.ir[write]
    //    _A.nc
    //    _A.nr

    // minimize
230
    model.add(IloMinimize(env_, objective));
231
232
233
234
235

    // 4. solve
    IloCplex cplex(model);
    cplex.setParam(IloCplex::TiLim, gurobi_max_time_);

Max Lyon's avatar
Max Lyon committed
236
237
238
239
240
241
242
243
244
#ifdef 0
    // set parameters comparable to CoMISo
    {
      cplex.setParam(IloCplex::MIPSearch  , 1);  // Traditional Branch-and-Cut
      cplex.setParam(IloCplex::NodeSel    , 0);  // Depth-First
      cplex.setParam(IloCplex::VarSel     , -1); // closest to integer
      cplex.setParam(IloCplex::MIPEmphasis, 1);  // concentrate on feasibility
    }
#endif
245
246
247
248
249

    cplex.solve();

    // 5. store result
    _x.resize(n);
250
    for (unsigned int i = 0; i < n; ++i)
251
252
      _x[i] = cplex.getValue(vars[i]);

253
    DEB_out(2, "CPLEX objective: " << cplex.getObjValue() << "\n");
254
255
256
  }
  catch (IloException& e)
  {
257
    PROGRESS_RESUME_ABORT; // resume a processed abort request
258
    DEB_warning(2, "CPLEX Concert exception caught: " << e.getMessage())
259
260
261
  }
  catch (...)
  {
262
    PROGRESS_RESUME_ABORT; // resume a processed abort request
263
    DEB_warning(1, "CPLEX Unknown exception caught")
264
265
266
267
268
269
270
  }

#else
  DEB_out(1, "CPLEX solver is not available, please install it...\n")
#endif
}

271
272
//-----------------------------------------------------------------------------

273
void MISolver::solve_no_rounding(CSCMatrix& _A, Vecd& _x, Vecd& _rhs)
David Bommes's avatar
David Bommes committed
274
{
275
276
  COMISO_THROW_if(
      !direct_solver_.calc_system_gmm(_A), UNSPECIFIED_EIGEN_FAILURE);
277
  COMISO_THROW_if(!direct_solver_.solve(_x, _rhs), UNSPECIFIED_EIGEN_FAILURE);
David Bommes's avatar
David Bommes committed
278
279
280
281
}

//-----------------------------------------------------------------------------

282
void MISolver::resolve(Vecd& _x, Vecd& _rhs) { direct_solver_.solve(_x, _rhs); }
283
284
285

//-----------------------------------------------------------------------------

286
287
void MISolver::solve_direct_rounding(
    CSCMatrix& _A, Vecd& _x, Vecd& _rhs, Veci& _to_round)
David Bommes's avatar
David Bommes committed
288
{
Max Lyon's avatar
Max Lyon committed
289
290
  DEB_enter_func;
  Veci to_round(_to_round);
David Bommes's avatar
David Bommes committed
291
292
293
294
  // copy to round vector and make it unique
  std::sort(to_round.begin(), to_round.end());
  Veci::iterator last_unique;
  last_unique = std::unique(to_round.begin(), to_round.end());
Max Lyon's avatar
Max Lyon committed
295
  size_t r = (size_t)(last_unique - to_round.begin());
296
  to_round.resize(r);
David Bommes's avatar
David Bommes committed
297
298
299

  // initalize old indices
  Veci old_idx(_rhs.size());
300
  for (unsigned int i = 0; i < old_idx.size(); ++i)
David Bommes's avatar
David Bommes committed
301
    old_idx[i] = i;
302
303
  direct_solver_.calc_system_gmm(_A);
  direct_solver_.solve(_x, _rhs);
David Bommes's avatar
David Bommes committed
304

305
#ifdef COMISO_MISOLVER_PERFORMANCE_TEST
306
307
  // check solver performance (only for testing!!!)
  {
308
    Base::StopWatch sw;
309
310

    // performance comparison code
311
#if (COMISO_SUITESPARSE_SPQR_AVAILABLE)
312
313
314
315
    {
      sw.start();
      COMISO::SparseQRSolver spqr;
      spqr.calc_system_gmm(_A);
316
      std::cerr << "SparseQR factor took: " << sw.stop() / 1000.0 << "s\n";
317
318
      Vecd x2(_x);
      sw.start();
319
320
      spqr.solve(x2, _rhs);
      std::cerr << "SparseQR solve took: " << sw.stop() / 1000.0 << "s\n";
321
      Vecd res(_x);
322
323
324
      gmm::add(_x, gmm::scaled(x2, -1.0), res);
      std::cerr << "DIFFERENCE IN RESULT: " << gmm::vect_norm2(res)
                << std::endl;
325
    }
David Bommes's avatar
changed    
David Bommes committed
326
#endif
327
328

    // performance comparison code
329
#if (COMISO_SUITESPARSE_AVAILABLE)
330
331
332
333
    {
      sw.start();
      COMISO::UMFPACKSolver umf;
      umf.calc_system_gmm(_A);
334
      std::cerr << "UMFPack factor took: " << sw.stop() / 1000.0 << "s\n";
335
336
      Vecd x3(_x);
      sw.start();
337
338
      umf.solve(x3, _rhs);
      std::cerr << "UMFPack solve took: " << sw.stop() / 1000.0 << "s\n";
339
      Vecd res2(_x);
340
341
342
      gmm::add(_x, gmm::scaled(x3, -1.0), res2);
      std::cerr << "UMFPACK DIFFERENCE IN RESULT: " << gmm::vect_norm2(res2)
                << std::endl;
343
344
345
346
347
348
349
    }

    // performance comparison code
    {
      sw.start();
      COMISO::CholmodSolver chol;
      chol.calc_system_gmm(_A);
350
      std::cerr << "Choldmod factor took: " << sw.stop() / 1000.0 << "s\n";
351
352
      Vecd x4(_x);
      sw.start();
353
354
      chol.solve(x4, _rhs);
      std::cerr << "Choldmod solve took: " << sw.stop() / 1000.0 << "s\n";
355
      Vecd res(_x);
356
357
358
      gmm::add(_x, gmm::scaled(x4, -1.0), res);
      std::cerr << "DIFFERENCE IN RESULT: " << gmm::vect_norm2(res)
                << std::endl;
359
    }
360
361
#endif

362
#if (COMISO_EIGEN3_AVAILABLE)
363
364
365
366
367
    // performance comparison code
    {
      sw.start();
      COMISO::EigenLDLTSolver ldlt;
      ldlt.calc_system_gmm(_A);
368
369
      DEB_out(2, "Eigen LDLT factor took: " << sw.stop() / 1000.0 << "s\n")
          Vecd x5(_x);
370
      sw.start();
371
372
373
374
375
      ldlt.solve(x5, _rhs);
      DEB_out(2, "Eigen LDLT solve took: " << sw.stop() / 1000.0 << "s\n")
          Vecd res(_x);
      gmm::add(_x, gmm::scaled(x5, -1.0), res);
      DEB_warning(2, "DIFFERENCE IN RESULT: " << gmm::vect_norm2(res))
376
377
    }
#endif
378
  }
379
#endif
David Bommes's avatar
David Bommes committed
380
381
382

  // round and eliminate variables
  Vecui elim_i;
383
384
  Vecd elim_v;
  for (unsigned int i = 0; i < to_round.size(); ++i)
David Bommes's avatar
David Bommes committed
385
386
387
388
389
390
391
392
  {
    _x[to_round[i]] = ROUND(_x[to_round[i]]);
    elim_i.push_back(to_round[i]);
    elim_v.push_back(_x[to_round[i]]);
    // update old idx
    old_idx[to_round[i]] = -1;
  }

393
394
  Veci::iterator new_end = std::remove(old_idx.begin(), old_idx.end(), -1);
  old_idx.resize(new_end - old_idx.begin());
David Bommes's avatar
David Bommes committed
395
396
  // eliminate vars from linear system
  Vecd xr(_x);
397
  COMISO_GMM::eliminate_csc_vars2(elim_i, elim_v, _A, xr, _rhs);
David Bommes's avatar
David Bommes committed
398

399
  // std::cerr << "size A: " << gmm::mat_nrows(_A) << " " << gmm::mat_ncols(_A)
David Bommes's avatar
David Bommes committed
400
401
402
403
404
  // 	    << std::endl;
  // std::cerr << "size x  : " << xr.size() << std::endl;
  // std::cerr << "size rhs: " << _rhs.size() << std::endl;

  // final full solution
405
  if (gmm::mat_ncols(_A) > 0)
David Bommes's avatar
David Bommes committed
406
  {
407
408
    //    direct_solver_.update_system_gmm(_A);
    direct_solver_.calc_system_gmm(_A);
409
    direct_solver_.solve(xr, _rhs);
David Bommes's avatar
David Bommes committed
410
411
412
  }

  // store solution values to result vector
413
414
  for (unsigned int i = 0; i < old_idx.size(); ++i)
    _x[old_idx[i]] = xr[i];
David Bommes's avatar
David Bommes committed
415
416
417
418
}

//-----------------------------------------------------------------------------

419
420
void MISolver::solve_iterative(
    CSCMatrix& _A, Vecd& _x, Vecd& _rhs, Veci& _to_round, bool _fixed_order)
David Bommes's avatar
David Bommes committed
421
{
422
  DEB_enter_func;
David Bommes's avatar
David Bommes committed
423
  // StopWatch
424
  Base::StopWatch sw;
David Bommes's avatar
David Bommes committed
425
426
  double time_search_next_integer = 0;

427
  // some statistics
David Bommes's avatar
David Bommes committed
428
  n_local_ = 0;
429
430
  n_cg_ = 0;
  n_full_ = 0;
David Bommes's avatar
David Bommes committed
431
432
433

  // reset cholmod step flag
  cholmod_step_done_ = false;
434
435
436

  Veci to_round(_to_round);
  // if the order is not fixed, uniquify the indices
437
  if (!_fixed_order)
438
439
440
441
442
  {
    // copy to round vector and make it unique
    std::sort(to_round.begin(), to_round.end());
    Veci::iterator last_unique;
    last_unique = std::unique(to_round.begin(), to_round.end());
Max Lyon's avatar
Max Lyon committed
443
    size_t r = (size_t)(last_unique - to_round.begin());
444
    to_round.resize(r);
445
446
447
448
  }

  // initalize old indices
  Veci old_idx(_rhs.size());
449
  for (unsigned int i = 0; i < old_idx.size(); ++i)
450
451
    old_idx[i] = i;

452
  if (initial_full_solution_)
453
  {
454
455
    DEB_out_if(noisy_ > 2, 2, "initial full solution\n")
        direct_solver_.calc_system_gmm(_A);
456
    direct_solver_.solve(_x, _rhs);
457

David Bommes's avatar
David Bommes committed
458
    cholmod_step_done_ = true;
459

David Bommes's avatar
David Bommes committed
460
461
    ++n_full_;
  }
462
463
464
465
466
467
468
469

  // neighbors for local optimization
  Vecui neigh_i;

  // Vector for reduced solution
  Vecd xr(_x);

  // loop until solution computed
470
  for (unsigned int i = 0; i < to_round.size(); ++i)
471
  {
472
473
474
475
    DEB_out_if(noisy_ > 0, 2,
        "Integer DOF's left: " << to_round.size() - (i + 1)
                               << " ") DEB_out_if(noisy_ > 1, 1,
        "residuum_norm: " << COMISO_GMM::residuum_norm(_A, xr, _rhs) << "\n");
476
477
478
479

    // index to eliminate
    unsigned int i_best = 0;

David Bommes's avatar
David Bommes committed
480
481
482
    // position in round vector
    unsigned int tr_best = 0;

483
    if (_fixed_order) // if order is fixed, simply get next index from to_round
484
485
486
    {
      i_best = to_round[i];
    }
487
    else // else search for best rounding candidate
488
    {
David Bommes's avatar
David Bommes committed
489
      sw.start();
490
      // find index yielding smallest rounding error
491
492
      double r_best = FLT_MAX;
      for (unsigned int j = 0; j < to_round.size(); ++j)
493
      {
494
        if (to_round[j] != -1)
495
496
        {
          int cur_idx = to_round[j];
497
498
          double rnd_error = fabs(ROUND(xr[cur_idx]) - xr[cur_idx]);
          if (rnd_error < r_best)
499
          {
500
501
502
            i_best = cur_idx;
            r_best = rnd_error;
            tr_best = j;
503
504
505
          }
        }
      }
David Bommes's avatar
David Bommes committed
506
      time_search_next_integer += sw.stop();
507
508
509
510
    }

    // store rounded value
    double rnd_x = ROUND(xr[i_best]);
511
    _x[old_idx[i_best]] = rnd_x;
512
513
514
515

    // compute neighbors
    neigh_i.clear();
    Col col = gmm::mat_const_col(_A, i_best);
516
517
518
519
520
    ColIter it = gmm::vect_const_begin(col);
    ColIter ite = gmm::vect_const_end(col);
    for (; it != ite; ++it)
    {
      if (it.index() != i_best)
Andreas Fabri's avatar
Andreas Fabri committed
521
        neigh_i.push_back(static_cast<int>(it.index()));
522
    }
523
    // eliminate var
David Bommes's avatar
David Bommes committed
524
525
    COMISO_GMM::fix_var_csc_symmetric(i_best, rnd_x, _A, xr, _rhs);
    to_round[tr_best] = -1;
526

David Bommes's avatar
David Bommes committed
527
528
    // 3-stage update of solution w.r.t. roundings
    // local GS / CG / SparseCholesky
529
    update_solution(_A, xr, _rhs, neigh_i);
David Bommes's avatar
David Bommes committed
530
  }
531

David Bommes's avatar
David Bommes committed
532
  // final full solution?
533
  if (final_full_solution_)
David Bommes's avatar
David Bommes committed
534
  {
535
    DEB_line_if(noisy_ > 2, 2, "final full solution");
536

537
    if (gmm::mat_ncols(_A) > 0)
538
    {
539
540
      if (cholmod_step_done_)
        direct_solver_.update_system_gmm(_A);
David Bommes's avatar
David Bommes committed
541
      else
542
        direct_solver_.calc_system_gmm(_A);
543

544
      direct_solver_.solve(xr, _rhs);
David Bommes's avatar
David Bommes committed
545
      ++n_full_;
546
    }
David Bommes's avatar
David Bommes committed
547
548
549
  }

  // store solution values to result vector
550
551
  for (unsigned int i = 0; i < old_idx.size(); ++i)
    _x[old_idx[i]] = xr[i];
David Bommes's avatar
David Bommes committed
552
553

  // output statistics
554
555
556
557
558
559
560
561
  DEB_out_if(stats_, 2,
      " *** Statistics of MiSo Solver ***"
          << "\n Number of CG    iterations  = " << n_cg_
          << "\n Number of LOCAL iterations  = " << n_local_
          << "\n Number of FULL  iterations  = " << n_full_
          << "\n Number of ROUNDING          = " << _to_round.size()
          << "\n time searching next integer = "
          << time_search_next_integer / 1000.0 << "s\n\n");
David Bommes's avatar
David Bommes committed
562
}
563

David Bommes's avatar
David Bommes committed
564
565
//-----------------------------------------------------------------------------

566
567
void MISolver::update_solution(
    const CSCMatrix& _A, Vecd& _x, const Vecd& _rhs, const Vecui& _neigh_i)
David Bommes's avatar
David Bommes committed
568
{
569
570
  DEB_enter_func;
  bool converged = false; // set to not converged
David Bommes's avatar
David Bommes committed
571

572
  if (max_local_iters_ > 0) // compute new solution
David Bommes's avatar
David Bommes committed
573
  {
574
    DEB_out(6, "use local iteration ");
David Bommes's avatar
David Bommes committed
575

576
    int n_its = max_local_iters_;
David Bommes's avatar
David Bommes committed
577
    double tolerance = max_local_error_;
578
579
    converged =
        siter_.gauss_seidel_local(_A, _x, _rhs, _neigh_i, n_its, tolerance);
David Bommes's avatar
David Bommes committed
580
581
582
583
584

    ++n_local_;
  }

  // conjugate gradient
585
  if (!converged && max_cg_iters_ > 0)
David Bommes's avatar
David Bommes committed
586
  {
587
    DEB_out(6, ", cg ");
David Bommes's avatar
David Bommes committed
588
589
590

    int max_cg_iters = max_cg_iters_;
    double tolerance = max_cg_error_;
591
592
    converged =
        siter_.conjugate_gradient(_A, _x, _rhs, max_cg_iters, tolerance);
David Bommes's avatar
David Bommes committed
593

594
595
596
    DEB_out(6, "( converged " << converged << " "
                              << " iters " << max_cg_iters << " "
                              << " res_norm " << tolerance << "\n");
David Bommes's avatar
David Bommes committed
597
598
599
    ++n_cg_;
  }

600
  if (!converged && iter_full_solution_ && gmm::mat_ncols(_A) > 0)
David Bommes's avatar
David Bommes committed
601
  {
602
603
604
605
    DEB_out(6, ", full ");
    if (cholmod_step_done_)
      direct_solver_.update_system_gmm(_A);
    else
David Bommes's avatar
David Bommes committed
606
    {
607
608
      direct_solver_.calc_system_gmm(_A);
      cholmod_step_done_ = true;
609
    }
610
611
612
613
    // const_cast<> to work around broken const correctness in Eigen
    direct_solver_.solve(_x, const_cast<Vecd&>(_rhs)); 

    ++n_full_;
David Bommes's avatar
David Bommes committed
614
615
  }

616
  DEB_line(6, "");
David Bommes's avatar
David Bommes committed
617
618
619
}

//-----------------------------------------------------------------------------
620

621
622
void MISolver::solve_multiple_rounding(
    CSCMatrix& _A, Vecd& _x, Vecd& _rhs, const Veci& _to_round)
David Bommes's avatar
David Bommes committed
623
{
624
  DEB_time_func_def;
David Bommes's avatar
David Bommes committed
625
626
  // some statistics
  n_local_ = 0;
627
628
  n_cg_ = 0;
  n_full_ = 0;
David Bommes's avatar
David Bommes committed
629
630
631
632
633

  // reset cholmod step flag
  cholmod_step_done_ = false;

  Veci to_round(_to_round);
634
635
636
637
638
  { // copy to round vector and make it unique
    std::sort(to_round.begin(), to_round.end());
    const auto last_unique = std::unique(to_round.begin(), to_round.end());
    to_round.resize(last_unique - to_round.begin());
  }
David Bommes's avatar
David Bommes committed
639

640
  if (initial_full_solution_)
David Bommes's avatar
David Bommes committed
641
  {
642
643
644
645
646
647
    DEB_out_if(noisy_ > 2, 2, "initial full solution\n");
    // TODO: we can throw more specific outcomes in the body of the
    // functions below
    COMISO_THROW_if(
        !direct_solver_.calc_system_gmm(_A), UNSPECIFIED_EIGEN_FAILURE);
    COMISO_THROW_if(!direct_solver_.solve(_x, _rhs), UNSPECIFIED_EIGEN_FAILURE);
David Bommes's avatar
David Bommes committed
648
649
650
651
652
653

    cholmod_step_done_ = true;

    ++n_full_;
  }

654
655
656
657
  //  data context for the main rounding cycle
  Vecui neigh_i;            // neighbors for local optimization
  std::vector<int> tr_best; // best indices to round
  RoundingQueue rndg_queue(multiple_rounding_threshold_); // rounding queue
David Bommes's avatar
David Bommes committed
658

659
660
661
662
#ifdef DEB_ON
  Base::StopWatch sw; // additional timer for the search next integer time
  double time_search_next_integer = 0;
#endif // DEB_ON
David Bommes's avatar
David Bommes committed
663
664

  // loop until solution computed
665
  while (!to_round.empty())
David Bommes's avatar
David Bommes committed
666
  {
667
668
    DEB_line(11, "Integer DOF's left: " << to_round.size());
    DEB_line(11, "residuum_norm: " << COMISO_GMM::residuum_norm(_A, _x, _rhs));
David Bommes's avatar
David Bommes committed
669

670
671
    DEB_only(sw.start());
    rndg_queue.clear();
672

David Bommes's avatar
David Bommes committed
673
    // find index yielding smallest rounding error
674
    for (unsigned int j = 0; j < to_round.size(); ++j)
David Bommes's avatar
David Bommes committed
675
    {
676
677
678
      const auto xr_j = _x[to_round[j]];
      const auto rnd_error = fabs(ROUND(xr_j) - xr_j);
      rndg_queue.add(j, rnd_error);
679
    }
680
    rndg_queue.get_ids(tr_best);
681

682
    DEB_only(time_search_next_integer += sw.stop());
David Bommes's avatar
David Bommes committed
683
684

    // nothing more to do?
685
    if (tr_best.empty())
David Bommes's avatar
David Bommes committed
686
687
      break;

688
    DEB_line(11, "rounding " << tr_best.size() << " variables simultaneously");
David Bommes's avatar
David Bommes committed
689
690
691
692

    // clear neigh for local update
    neigh_i.clear();

693
    for (const auto tr_indx : tr_best)
David Bommes's avatar
David Bommes committed
694
    {
695
696
      const unsigned i_cur = static_cast<unsigned>(to_round[tr_indx]);
      const double rnd_x = ROUND(_x[i_cur]); // store rounded value
David Bommes's avatar
David Bommes committed
697
698

      // compute neighbors
699
700
701
702
703
704
705
706
707
708
      const Col col = gmm::mat_const_col(_A, i_cur);
      for (auto it = gmm::vect_const_begin(col), ite = gmm::vect_const_end(col);
           it != ite; ++it)
      {
        if (it.index() != (unsigned int)i_cur)
          neigh_i.push_back(static_cast<int>(it.index()));
      }
      COMISO_GMM::fix_var_csc_symmetric(
          i_cur, rnd_x, _A, _x, _rhs); // eliminate
      to_round[tr_indx] = -1;
David Bommes's avatar
David Bommes committed
709
710
711
712
    }

    // 3-stage update of solution w.r.t. roundings
    // local GS / CG / SparseCholesky
713
    update_solution(_A, _x, _rhs, neigh_i);
714

715
716
    // remove all rounded indices from the to_round[] array
    for (const auto tr_indx : tr_best)
717
    {
718
719
720
721
722
723
      while (!to_round.empty() && to_round.back() < 0)
        to_round.pop_back(); // remove entries that are already rounded
      if (tr_indx >= to_round.size())
        continue; // already removed
      std::swap(to_round[tr_indx], to_round.back());
      to_round.pop_back();
724
725
726
    }
  }

727
728
  // final full solution?
  if (final_full_solution_ && gmm::mat_ncols(_A) > 0)
729
  {
730
731
732
733
734
735
736
737
    DEB_line(3, "final full solution");
    if (cholmod_step_done_)
      direct_solver_.update_system_gmm(_A);
    else
      direct_solver_.calc_system_gmm(_A);

    direct_solver_.solve(_x, _rhs);
    ++n_full_;
738
739
  }

740
741
742
743
744
745
746
747
  DEB_out_if(stats_, 2, // output statistics
      " *** Statistics of MiSo Solver ***"
          << "\n Number of CG    iterations  = " << n_cg_
          << "\n Number of LOCAL iterations  = " << n_local_
          << "\n Number of FULL  iterations  = " << n_full_
          << "\n Number of ROUNDING          = " << _to_round.size()
          << "\n time searching next integer = "
          << time_search_next_integer / 1000.0 << "s\n\n");
748
749
}

David Bommes's avatar
David Bommes committed
750
751
//-----------------------------------------------------------------------------

752
753
void MISolver::solve_gurobi(
    CSCMatrix& _A, Vecd& _x, Vecd& _rhs, Veci& _to_round)
David Bommes's avatar
David Bommes committed
754
{
755
  DEB_enter_func;
David Bommes's avatar
David Bommes committed
756
757
758
759
#if COMISO_GUROBI_AVAILABLE

  // get round-indices in set
  std::set<int> to_round;
760
  for (unsigned int i = 0; i < _to_round.size(); ++i)
David Bommes's avatar
David Bommes committed
761
762
    to_round.insert(_to_round[i]);

763
764
  try
  {
David Bommes's avatar
David Bommes committed
765
766
767
768
769
770
771
772
773
774
775
    GRBEnv env = GRBEnv();

    GRBModel model = GRBModel(env);

    // set time limite
    model.getEnv().set(GRB_DoubleParam_TimeLimit, gurobi_max_time_);

    unsigned int n = _rhs.size();

    // 1. allocate variables
    std::vector<GRBVar> vars;
776
777
778
779
    for (unsigned int i = 0; i < n; ++i)
      if (to_round.count(i))
        vars.push_back(
            model.addVar(-GRB_INFINITY, GRB_INFINITY, 0.0, GRB_INTEGER));
David Bommes's avatar
David Bommes committed
780
      else
781
782
        vars.push_back(
            model.addVar(-GRB_INFINITY, GRB_INFINITY, 0.0, GRB_CONTINUOUS));
David Bommes's avatar
David Bommes committed
783
784
785
786
787
788
789
790
791

    // Integrate new variables
    model.update();

    // 2. setup_energy

    // build objective function from linear system E = x^tAx - 2x^t*rhs
    GRBQuadExpr objective;

792
793
    for (unsigned int i = 0; i < _A.nc; ++i)
      for (unsigned int j = _A.jc[i]; j < _A.jc[i + 1]; ++j)
David Bommes's avatar
David Bommes committed
794
      {
795
        objective += _A.pr[j] * vars[_A.ir[j]] * vars[i];
David Bommes's avatar
David Bommes committed
796
      }
797
798
    for (unsigned int i = 0; i < n; ++i)
      objective -= 2 * _rhs[i] * vars[i];
David Bommes's avatar
David Bommes committed
799

800
801
802
803
804
    //    _A.jc[c+1]
    //    _A.pr[write]
    //    _A.ir[write]
    //    _A.nc
    //    _A.nr
David Bommes's avatar
David Bommes committed
805
806
807
808
809
810
811
812
813
814

    // minimize
    model.set(GRB_IntAttr_ModelSense, 1);
    model.setObjective(objective);

    // 4. solve
    model.optimize();

    // 5. store result
    _x.resize(n);
815
    for (unsigned int i = 0; i < n; ++i)
David Bommes's avatar
David Bommes committed
816
817
      _x[i] = vars[i].get(GRB_DoubleAttr_X);

818
819
    DEB_out(
        2, "GUROBI objective: " << model.get(GRB_DoubleAttr_ObjVal) << "\n");
David Bommes's avatar
David Bommes committed
820
  }
821
  catch (GRBException& e)
David Bommes's avatar
David Bommes committed
822
  {
823
    PROGRESS_RESUME_ABORT; // resume a processed abort request
824
825
    DEB_warning(2,
        "Error code = " << e.getErrorCode() << "[" << e.getMessage() << "]\n");
David Bommes's avatar
David Bommes committed
826
  }
827
  catch (...)
David Bommes's avatar
David Bommes committed
828
  {
829
    PROGRESS_RESUME_ABORT; // resume a processed abort request
830
    DEB_warning(1, "Exception during optimization");
David Bommes's avatar
David Bommes committed
831
832
833
  }

#else
834
  DEB_out(1, "GUROBI solver is not available, please install it...\n");
David Bommes's avatar
David Bommes committed
835
836
837
#endif
}

838
839
//----------------------------------------------------------------------------

840
void MISolver::show_options_dialog()
841
{
842
  DEB_enter_func;
843
#if (COMISO_QT4_AVAILABLE)
844
845
  MISolverDialog* pd = new MISolverDialog(*this);
  pd->exec();
David Bommes's avatar
changed    
David Bommes committed
846
#else
847
  DEB_warning(1, "Qt not available to show solver dialog!!!");
848
849
850
#endif
}

851
// end namespace COMISO
852
} // namespace COMISO