VolVolume.hpp
Go to the documentation of this file.
1 // Copyright (C) 2000, International Business Machines
2 // Corporation and others. All Rights Reserved.
3 #ifndef __VOLUME_HPP__
4 #define __VOLUME_HPP__
5 
6 #include <cfloat>
7 #include <algorithm>
8 #include <cstdio>
9 #include <cmath>
10 
11 #ifndef VOL_DEBUG
12 // When VOL_DEBUG is 1, we check vector indices
13 #define VOL_DEBUG 0
14 #endif
15 
16 template <class T> static inline T
17 VolMax(register const T x, register const T y) {
18  return ((x) > (y)) ? (x) : (y);
19 }
20 
21 template <class T> static inline T
22 VolAbs(register const T x) {
23  return ((x) > 0) ? (x) : -(x);
24 }
25 
26 //############################################################################
27 
28 #if defined(VOL_DEBUG) && (VOL_DEBUG != 0)
29 #define VOL_TEST_INDEX(i, size) \
30 { \
31  if ((i) < 0 || (i) >= (size)) { \
32  printf("bad VOL_?vector index\n"); \
33  abort(); \
34  } \
35 }
36 #define VOL_TEST_SIZE(size) \
37 { \
38  if (s <= 0) { \
39  printf("bad VOL_?vector size\n"); \
40  abort(); \
41  } \
42 }
43 #else
44 #define VOL_TEST_INDEX(i, size)
45 #define VOL_TEST_SIZE(size)
46 #endif
47 
48 //############################################################################
49 
50 class VOL_dvector;
51 class VOL_ivector;
52 class VOL_primal;
53 class VOL_dual;
54 class VOL_swing;
55 class VOL_alpha_factor;
56 class VOL_vh;
57 class VOL_indc;
58 class VOL_user_hooks;
59 class VOL_problem;
60 
61 //############################################################################
62 
66 struct VOL_parms {
68  double lambdainit;
70  double alphainit;
72  double alphamin;
74  double alphafactor;
75 
77  double ubinit;
78 
86  double granularity;
87 
96 
99 
110  int printflag;
112  int printinvl;
114  int heurinvl;
115 
125 
127  int alphaint;
128 
131 };
132 
133 //############################################################################
134 
143 class VOL_dvector {
144 public:
146  double* v;
148  int sz;
149 
150 public:
152  VOL_dvector(const int s) {
153  VOL_TEST_SIZE(s);
154  v = new double[sz = s];
155  }
157  VOL_dvector() : v(0), sz(0) {}
159  VOL_dvector(const VOL_dvector& x) : v(0), sz(0) {
160  sz = x.sz;
161  if (sz > 0) {
162  v = new double[sz];
163  std::copy(x.v, x.v + sz, v);
164  }
165  }
167  ~VOL_dvector() { delete[] v; }
168 
170  inline int size() const {return sz;}
171 
173  inline double& operator[](const int i) {
174  VOL_TEST_INDEX(i, sz);
175  return v[i];
176  }
177 
179  inline double operator[](const int i) const {
180  VOL_TEST_INDEX(i, sz);
181  return v[i];
182  }
183 
186  inline void clear() {
187  delete[] v;
188  v = 0;
189  sz = 0;
190  }
193  inline void cc(const double gamma, const VOL_dvector& w) {
194  if (sz != w.sz) {
195  printf("bad VOL_dvector sizes\n");
196  abort();
197  }
198  double * p_v = v - 1;
199  const double * p_w = w.v - 1;
200  const double * const p_e = v + sz;
201  const double one_gamma = 1.0 - gamma;
202  while ( ++p_v != p_e ){
203  *p_v = one_gamma * (*p_v) + gamma * (*++p_w);
204  }
205  }
206 
209  inline void allocate(const int s) {
210  VOL_TEST_SIZE(s);
211  delete[] v;
212  v = new double[sz = s];
213  }
214 
216  inline void swap(VOL_dvector& w) {
217  std::swap(v, w.v);
218  std::swap(sz, w.sz);
219  }
220 
222  VOL_dvector& operator=(const VOL_dvector& w);
224  VOL_dvector& operator=(const double w);
225 };
226 
227 //-----------------------------------------------------------------------------
237 class VOL_ivector {
238 public:
240  int* v;
242  int sz;
243 public:
245  VOL_ivector(const int s) {
246  VOL_TEST_SIZE(s);
247  v = new int[sz = s];
248  }
250  VOL_ivector() : v(0), sz(0) {}
253  sz = x.sz;
254  if (sz > 0) {
255  v = new int[sz];
256  std::copy(x.v, x.v + sz, v);
257  }
258  }
261  delete [] v;
262  }
263 
265  inline int size() const { return sz; }
267  inline int& operator[](const int i) {
268  VOL_TEST_INDEX(i, sz);
269  return v[i];
270  }
271 
273  inline int operator[](const int i) const {
274  VOL_TEST_INDEX(i, sz);
275  return v[i];
276  }
277 
280  inline void clear() {
281  delete[] v;
282  v = 0;
283  sz = 0;
284  }
285 
288  inline void allocate(const int s) {
289  VOL_TEST_SIZE(s);
290  delete[] v;
291  v = new int[sz = s];
292  }
293 
295  inline void swap(VOL_ivector& w) {
296  std::swap(v, w.v);
297  std::swap(sz, w.sz);
298  }
299 
303  VOL_ivector& operator=(const int w);
304 };
305 
306 //############################################################################
307 // A class describing a primal solution. This class is used only internally
308 class VOL_primal {
309 public:
310  // objective value of this primal solution
311  double value;
312  // the largest of the v[i]'s
313  double viol;
314  // primal solution
316  // v=b-Ax, for the relaxed constraints
318 
319  VOL_primal(const int psize, const int dsize) : x(psize), v(dsize) {}
320  VOL_primal(const VOL_primal& primal) :
321  value(primal.value), viol(primal.viol), x(primal.x), v(primal.v) {}
323  inline VOL_primal& operator=(const VOL_primal& p) {
324  if (this == &p)
325  return *this;
326  value = p.value;
327  viol = p.viol;
328  x = p.x;
329  v = p.v;
330  return *this;
331  }
332 
333  // convex combination. data members in this will be overwritten
334  // convex combination between two primal solutions
335  // x <-- alpha x + (1 - alpha) p.x
336  // v <-- alpha v + (1 - alpha) p.v
337  inline void cc(const double alpha, const VOL_primal& p) {
338  value = alpha * p.value + (1.0 - alpha) * value;
339  x.cc(alpha, p.x);
340  v.cc(alpha, p.v);
341  }
342  // find maximum of v[i]
343  void find_max_viol(const VOL_dvector& dual_lb,
344  const VOL_dvector& dual_ub);
345 };
346 
347 //-----------------------------------------------------------------------------
348 // A class describing a dual solution. This class is used only internally
349 class VOL_dual {
350 public:
351  // lagrangian value
352  double lcost;
353  // reduced costs * (pstar-primal)
354  double xrc;
355  // this information is only printed
356  // dual vector
358 
359  VOL_dual(const int dsize) : u(dsize) { u = 0.0;}
360  VOL_dual(const VOL_dual& dual) :
361  lcost(dual.lcost), xrc(dual.xrc), u(dual.u) {}
363  inline VOL_dual& operator=(const VOL_dual& p) {
364  if (this == &p)
365  return *this;
366  lcost = p.lcost;
367  xrc = p.xrc;
368  u = p.u;
369  return *this;
370  }
371  // dual step
372  void step(const double target, const double lambda,
373  const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
374  const VOL_dvector& v);
375  double ascent(const VOL_dvector& v, const VOL_dvector& last_u) const;
376  void compute_xrc(const VOL_dvector& pstarx, const VOL_dvector& primalx,
377  const VOL_dvector& rc);
378 
379 };
380 
381 
382 //############################################################################
383 /* here we check whether an iteration is green, yellow or red. Also according
384  to this information we decide whether lambda should be changed */
385 class VOL_swing {
386 private:
387  VOL_swing(const VOL_swing&);
388  VOL_swing& operator=(const VOL_swing&);
389 public:
392  int ngs, nrs, nys;
393  int rd;
394 
397  ngs = nrs = nys = 0;
398  }
400 
401  inline void cond(const VOL_dual& dual,
402  const double lcost, const double ascent, const int iter) {
403  double eps = 1.e-3;
404 
405  if (ascent > 0.0 && lcost > dual.lcost + eps) {
406  lastswing = green;
407  lastgreeniter = iter;
408  ++ngs;
409  rd = 0;
410  } else {
411  if (ascent <= 0 && lcost > dual.lcost) {
412  lastswing = yellow;
413  lastyellowiter = iter;
414  ++nys;
415  rd = 0;
416  } else {
417  lastswing = red;
418  lastrediter = iter;
419  ++nrs;
420  rd = 1;
421  }
422  }
423  }
424 
425  inline double
426  lfactor(const VOL_parms& parm, const double lambda, const int iter) {
427  double lambdafactor = 1.0;
428  double eps = 5.e-4;
429  int cons;
430 
431  switch (lastswing) {
432  case green:
433  cons = iter - VolMax(lastyellowiter, lastrediter);
434  if (parm.printflag & 4)
435  printf(" G: Consecutive Gs = %3d\n\n", cons);
436  if (cons >= parm.greentestinvl && lambda < 2.0) {
438  lambdafactor = 2.0;
439  if (parm.printflag & 2)
440  printf("\n ---- increasing lamda to %g ----\n\n",
441  lambda * lambdafactor);
442  }
443  break;
444 
445  case yellow:
446  cons = iter - VolMax(lastgreeniter, lastrediter);
447  if (parm.printflag & 4)
448  printf(" Y: Consecutive Ys = %3d\n\n", cons);
449  if (cons >= parm.yellowtestinvl) {
451  lambdafactor = 1.1;
452  if (parm.printflag & 2)
453  printf("\n **** increasing lamda to %g *****\n\n",
454  lambda * lambdafactor);
455  }
456  break;
457 
458  case red:
459  cons = iter - VolMax(lastgreeniter, lastyellowiter);
460  if (parm.printflag & 4)
461  printf(" R: Consecutive Rs = %3d\n\n", cons);
462  if (cons >= parm.redtestinvl && lambda > eps) {
464  lambdafactor = 0.67;
465  if (parm.printflag & 2)
466  printf("\n **** decreasing lamda to %g *****\n\n",
467  lambda * lambdafactor);
468  }
469  break;
470  }
471  return lambdafactor;
472  }
473 
474  inline void
475  print() {
476  printf("**** G= %i, Y= %i, R= %i ****\n", ngs, nys, nrs);
477  ngs = nrs = nys = 0;
478  }
479 };
480 
481 //############################################################################
482 /* alpha should be decreased if after some number of iterations the objective
483  has increased less that 1% */
485 private:
488 public:
489  double lastvalue;
490 
491  VOL_alpha_factor() {lastvalue = -DBL_MAX;}
493 
494  inline double factor(const VOL_parms& parm, const double lcost,
495  const double alpha) {
496  if (alpha < parm.alphamin)
497  return 1.0;
498  register const double ll = VolAbs(lcost);
499  const double x = ll > 10 ? (lcost-lastvalue)/ll : (lcost-lastvalue);
500  lastvalue = lcost;
501  return (x <= 0.01) ? parm.alphafactor : 1.0;
502  }
503 };
504 
505 //############################################################################
506 /* here we compute the norm of the conjugate direction -hh-, the norm of the
507  subgradient -norm-, the inner product between the subgradient and the
508  last conjugate direction -vh-, and the inner product between the new
509  conjugate direction and the subgradient */
510 class VOL_vh {
511 private:
512  VOL_vh(const VOL_vh&);
513  VOL_vh& operator=(const VOL_vh&);
514 public:
515  double hh;
516  double norm;
517  double vh;
518  double asc;
519 
520  VOL_vh(const double alpha,
521  const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
522  const VOL_dvector& v, const VOL_dvector& vstar,
523  const VOL_dvector& u);
525 };
526 
527 //############################################################################
528 /* here we compute different parameter to be printed. v2 is the square of
529  the norm of the subgradient. vu is the inner product between the dual
530  variables and the subgradient. vabs is the maximum absolute value of
531  the violations of pstar. asc is the inner product between the conjugate
532  direction and the subgradient */
533 class VOL_indc {
534 private:
535  VOL_indc(const VOL_indc&);
536  VOL_indc& operator=(const VOL_indc&);
537 public:
538  double v2;
539  double vu;
540  double vabs;
541  double asc;
542 
543 public:
544  VOL_indc(const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
545  const VOL_primal& primal, const VOL_primal& pstar,
546  const VOL_dual& dual);
548 };
549 
550 //#############################################################################
551 
559 public:
560  virtual ~VOL_user_hooks() {}
561 public:
562  // for all hooks: return value of -1 means that volume should quit
567  virtual int compute_rc(const VOL_dvector& u, VOL_dvector& rc) = 0;
568 
577  virtual int solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
578  double& lcost, VOL_dvector& x, VOL_dvector& v,
579  double& pcost) = 0;
586  virtual int heuristics(const VOL_problem& p,
587  const VOL_dvector& x, double& heur_val) = 0;
588 };
589 
590 //#############################################################################
591 
600 class VOL_problem {
601 private:
602  VOL_problem(const VOL_problem&);
604  void set_default_parm();
605  // ############ INPUT fields ########################
606 public:
610  VOL_problem();
613  VOL_problem(const char *filename);
615  ~VOL_problem();
617 
623  int solve(VOL_user_hooks& hooks, const bool use_preset_dual = false);
625 
626 private:
630  double alpha_;
632  double lambda_;
633  // This union is here for padding (so that data members would be
634  // double-aligned on x86 CPU
635  union {
637  int iter_;
638  double __pad0;
639  };
641 
642 public:
643 
647  double value;
655 
661  int psize;
663  int dsize;
671 
672 public:
676  int iter() const { return iter_; }
678  double alpha() const { return alpha_; }
680  double lambda() const { return lambda_; }
682 
683 private:
687  void read_params(const char* filename);
688 
690  int initialize(const bool use_preset_dual);
691 
693  void print_info(const int iter,
694  const VOL_primal& primal, const VOL_primal& pstar,
695  const VOL_dual& dual);
696 
699  double readjust_target(const double oldtarget, const double lcost) const;
700 
708  double power_heur(const VOL_primal& primal, const VOL_primal& pstar,
709  const VOL_dual& dual) const;
711 };
712 
713 #endif