Cgl  0.60.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
OsiTestSolver.hpp
Go to the documentation of this file.
1 // Copyright (C) 2000, International Business Machines
2 // Corporation and others. All Rights Reserved.
3 // This file is licensed under the terms of Eclipse Public License (EPL).
4 
5 // this is a copy of VolVolume (stable/1.1 rev. 233)
6 
7 #ifndef __OSITESTSOLVER_HPP__
8 #define __OSITESTSOLVER_HPP__
9 
10 #include <cfloat>
11 #include <algorithm>
12 #include <cstdio>
13 #include <cmath>
14 #include <cstdlib>
15 
16 #include "CoinFinite.hpp"
17 
18 #ifndef VOL_DEBUG
19 // When VOL_DEBUG is 1, we check vector indices
20 #define VOL_DEBUG 0
21 #endif
22 
23 template <class T> static inline T
24 VolMax(register const T x, register const T y) {
25  return ((x) > (y)) ? (x) : (y);
26 }
27 
28 template <class T> static inline T
29 VolAbs(register const T x) {
30  return ((x) > 0) ? (x) : -(x);
31 }
32 
33 //############################################################################
34 
35 #if defined(VOL_DEBUG) && (VOL_DEBUG != 0)
36 #define VOL_TEST_INDEX(i, size) \
37 { \
38  if ((i) < 0 || (i) >= (size)) { \
39  printf("bad VOL_?vector index\n"); \
40  abort(); \
41  } \
42 }
43 #define VOL_TEST_SIZE(size) \
44 { \
45  if (s <= 0) { \
46  printf("bad VOL_?vector size\n"); \
47  abort(); \
48  } \
49 }
50 #else
51 #define VOL_TEST_INDEX(i, size)
52 #define VOL_TEST_SIZE(size)
53 #endif
54 
55 //############################################################################
56 
57 class VOL_dvector;
58 class VOL_ivector;
59 class VOL_primal;
60 class VOL_dual;
61 class VOL_swing;
62 class VOL_alpha_factor;
63 class VOL_vh;
64 class VOL_indc;
65 class VOL_user_hooks;
66 class VOL_problem;
67 
68 //############################################################################
69 
73 struct VOL_parms {
75  double lambdainit;
77  double alphainit;
79  double alphamin;
81  double alphafactor;
82 
84  double ubinit;
85 
93  double granularity;
94 
103 
106 
117  int printflag;
119  int printinvl;
121  int heurinvl;
122 
132 
134  int alphaint;
135 
138 };
139 
140 //############################################################################
141 
150 class VOL_dvector {
151 public:
153  double* v;
155  int sz;
156 
157 public:
159  VOL_dvector(const int s) {
160  VOL_TEST_SIZE(s);
161  v = new double[sz = s];
162  }
164  VOL_dvector() : v(0), sz(0) {}
166  VOL_dvector(const VOL_dvector& x) : v(0), sz(0) {
167  sz = x.sz;
168  if (sz > 0) {
169  v = new double[sz];
170  std::copy(x.v, x.v + sz, v);
171  }
172  }
174  ~VOL_dvector() { delete[] v; }
175 
177  inline int size() const {return sz;}
178 
180  inline double& operator[](const int i) {
181  VOL_TEST_INDEX(i, sz);
182  return v[i];
183  }
184 
186  inline double operator[](const int i) const {
187  VOL_TEST_INDEX(i, sz);
188  return v[i];
189  }
190 
193  inline void clear() {
194  delete[] v;
195  v = 0;
196  sz = 0;
197  }
200  inline void cc(const double gamma, const VOL_dvector& w) {
201  if (sz != w.sz) {
202  printf("bad VOL_dvector sizes\n");
203  abort();
204  }
205  double * p_v = v - 1;
206  const double * p_w = w.v - 1;
207  const double * const p_e = v + sz;
208  const double one_gamma = 1.0 - gamma;
209  while ( ++p_v != p_e ){
210  *p_v = one_gamma * (*p_v) + gamma * (*++p_w);
211  }
212  }
213 
216  inline void allocate(const int s) {
217  VOL_TEST_SIZE(s);
218  delete[] v;
219  v = new double[sz = s];
220  }
221 
223  inline void swap(VOL_dvector& w) {
224  std::swap(v, w.v);
225  std::swap(sz, w.sz);
226  }
227 
229  VOL_dvector& operator=(const VOL_dvector& w);
231  VOL_dvector& operator=(const double w);
232 };
233 
234 //-----------------------------------------------------------------------------
244 class VOL_ivector {
245 public:
247  int* v;
249  int sz;
250 public:
252  VOL_ivector(const int s) {
253  VOL_TEST_SIZE(s);
254  v = new int[sz = s];
255  }
257  VOL_ivector() : v(0), sz(0) {}
260  sz = x.sz;
261  if (sz > 0) {
262  v = new int[sz];
263  std::copy(x.v, x.v + sz, v);
264  }
265  }
268  delete [] v;
269  }
270 
272  inline int size() const { return sz; }
274  inline int& operator[](const int i) {
275  VOL_TEST_INDEX(i, sz);
276  return v[i];
277  }
278 
280  inline int operator[](const int i) const {
281  VOL_TEST_INDEX(i, sz);
282  return v[i];
283  }
284 
287  inline void clear() {
288  delete[] v;
289  v = 0;
290  sz = 0;
291  }
292 
295  inline void allocate(const int s) {
296  VOL_TEST_SIZE(s);
297  delete[] v;
298  v = new int[sz = s];
299  }
300 
302  inline void swap(VOL_ivector& w) {
303  std::swap(v, w.v);
304  std::swap(sz, w.sz);
305  }
306 
310  VOL_ivector& operator=(const int w);
311 };
312 
313 //############################################################################
314 // A class describing a primal solution. This class is used only internally
315 class VOL_primal {
316 public:
317  // objective value of this primal solution
318  double value;
319  // the largest of the v[i]'s
320  double viol;
321  // primal solution
323  // v=b-Ax, for the relaxed constraints
325 
326  VOL_primal(const int psize, const int dsize) : x(psize), v(dsize) {}
327  VOL_primal(const VOL_primal& primal) :
328  value(primal.value), viol(primal.viol), x(primal.x), v(primal.v) {}
330  inline VOL_primal& operator=(const VOL_primal& p) {
331  if (this == &p)
332  return *this;
333  value = p.value;
334  viol = p.viol;
335  x = p.x;
336  v = p.v;
337  return *this;
338  }
339 
340  // convex combination. data members in this will be overwritten
341  // convex combination between two primal solutions
342  // x <-- alpha x + (1 - alpha) p.x
343  // v <-- alpha v + (1 - alpha) p.v
344  inline void cc(const double alpha, const VOL_primal& p) {
345  value = alpha * p.value + (1.0 - alpha) * value;
346  x.cc(alpha, p.x);
347  v.cc(alpha, p.v);
348  }
349  // find maximum of v[i]
350  void find_max_viol(const VOL_dvector& dual_lb,
351  const VOL_dvector& dual_ub);
352 };
353 
354 //-----------------------------------------------------------------------------
355 // A class describing a dual solution. This class is used only internally
356 class VOL_dual {
357 public:
358  // lagrangian value
359  double lcost;
360  // reduced costs * (pstar-primal)
361  double xrc;
362  // this information is only printed
363  // dual vector
365 
366  VOL_dual(const int dsize) : u(dsize) { u = 0.0;}
367  VOL_dual(const VOL_dual& dual) :
368  lcost(dual.lcost), xrc(dual.xrc), u(dual.u) {}
370  inline VOL_dual& operator=(const VOL_dual& p) {
371  if (this == &p)
372  return *this;
373  lcost = p.lcost;
374  xrc = p.xrc;
375  u = p.u;
376  return *this;
377  }
378  // dual step
379  void step(const double target, const double lambda,
380  const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
381  const VOL_dvector& v);
382  double ascent(const VOL_dvector& v, const VOL_dvector& last_u) const;
383  void compute_xrc(const VOL_dvector& pstarx, const VOL_dvector& primalx,
384  const VOL_dvector& rc);
385 
386 };
387 
388 
389 //############################################################################
390 /* here we check whether an iteration is green, yellow or red. Also according
391  to this information we decide whether lambda should be changed */
392 class VOL_swing {
393 private:
394  VOL_swing(const VOL_swing&);
395  VOL_swing& operator=(const VOL_swing&);
396 public:
399  int ngs, nrs, nys;
400  int rd;
401 
404  ngs = nrs = nys = 0;
405  }
407 
408  inline void cond(const VOL_dual& dual,
409  const double lcost, const double ascent, const int iter) {
410  double eps = 1.e-3;
411 
412  if (ascent > 0.0 && lcost > dual.lcost + eps) {
413  lastswing = green;
414  lastgreeniter = iter;
415  ++ngs;
416  rd = 0;
417  } else {
418  if (ascent <= 0 && lcost > dual.lcost) {
419  lastswing = yellow;
420  lastyellowiter = iter;
421  ++nys;
422  rd = 0;
423  } else {
424  lastswing = red;
425  lastrediter = iter;
426  ++nrs;
427  rd = 1;
428  }
429  }
430  }
431 
432  inline double
433  lfactor(const VOL_parms& parm, const double lambda, const int iter) {
434  double lambdafactor = 1.0;
435  double eps = 5.e-4;
436  int cons;
437 
438  switch (lastswing) {
439  case green:
440  cons = iter - VolMax(lastyellowiter, lastrediter);
441  if (parm.printflag & 4)
442  printf(" G: Consecutive Gs = %3d\n\n", cons);
443  if (cons >= parm.greentestinvl && lambda < 2.0) {
445  lambdafactor = 2.0;
446  if (parm.printflag & 2)
447  printf("\n ---- increasing lamda to %g ----\n\n",
448  lambda * lambdafactor);
449  }
450  break;
451 
452  case yellow:
453  cons = iter - VolMax(lastgreeniter, lastrediter);
454  if (parm.printflag & 4)
455  printf(" Y: Consecutive Ys = %3d\n\n", cons);
456  if (cons >= parm.yellowtestinvl) {
458  lambdafactor = 1.1;
459  if (parm.printflag & 2)
460  printf("\n **** increasing lamda to %g *****\n\n",
461  lambda * lambdafactor);
462  }
463  break;
464 
465  case red:
466  cons = iter - VolMax(lastgreeniter, lastyellowiter);
467  if (parm.printflag & 4)
468  printf(" R: Consecutive Rs = %3d\n\n", cons);
469  if (cons >= parm.redtestinvl && lambda > eps) {
471  lambdafactor = 0.67;
472  if (parm.printflag & 2)
473  printf("\n **** decreasing lamda to %g *****\n\n",
474  lambda * lambdafactor);
475  }
476  break;
477  }
478  return lambdafactor;
479  }
480 
481  inline void
482  print() {
483  printf("**** G= %i, Y= %i, R= %i ****\n", ngs, nys, nrs);
484  ngs = nrs = nys = 0;
485  }
486 };
487 
488 //############################################################################
489 /* alpha should be decreased if after some number of iterations the objective
490  has increased less that 1% */
492 private:
495 public:
496  double lastvalue;
497 
500 
501  inline double factor(const VOL_parms& parm, const double lcost,
502  const double alpha) {
503  if (alpha < parm.alphamin)
504  return 1.0;
505  register const double ll = VolAbs(lcost);
506  const double x = ll > 10 ? (lcost-lastvalue)/ll : (lcost-lastvalue);
507  lastvalue = lcost;
508  return (x <= 0.01) ? parm.alphafactor : 1.0;
509  }
510 };
511 
512 //############################################################################
513 /* here we compute the norm of the conjugate direction -hh-, the norm of the
514  subgradient -norm-, the inner product between the subgradient and the
515  last conjugate direction -vh-, and the inner product between the new
516  conjugate direction and the subgradient */
517 class VOL_vh {
518 private:
519  VOL_vh(const VOL_vh&);
520  VOL_vh& operator=(const VOL_vh&);
521 public:
522  double hh;
523  double norm;
524  double vh;
525  double asc;
526 
527  VOL_vh(const double alpha,
528  const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
529  const VOL_dvector& v, const VOL_dvector& vstar,
530  const VOL_dvector& u);
532 };
533 
534 //############################################################################
535 /* here we compute different parameter to be printed. v2 is the square of
536  the norm of the subgradient. vu is the inner product between the dual
537  variables and the subgradient. vabs is the maximum absolute value of
538  the violations of pstar. asc is the inner product between the conjugate
539  direction and the subgradient */
540 class VOL_indc {
541 private:
542  VOL_indc(const VOL_indc&);
543  VOL_indc& operator=(const VOL_indc&);
544 public:
545  double v2;
546  double vu;
547  double vabs;
548  double asc;
549 
550 public:
551  VOL_indc(const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
552  const VOL_primal& primal, const VOL_primal& pstar,
553  const VOL_dual& dual);
555 };
556 
557 //#############################################################################
558 
566 public:
567  virtual ~VOL_user_hooks() {}
568 public:
569  // for all hooks: return value of -1 means that volume should quit
574  virtual int compute_rc(const VOL_dvector& u, VOL_dvector& rc) = 0;
575 
584  virtual int solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
585  double& lcost, VOL_dvector& x, VOL_dvector& v,
586  double& pcost) = 0;
593  virtual int heuristics(const VOL_problem& p,
594  const VOL_dvector& x, double& heur_val) = 0;
595 };
596 
597 //#############################################################################
598 
607 class VOL_problem {
608 private:
609  VOL_problem(const VOL_problem&);
611  void set_default_parm();
612  // ############ INPUT fields ########################
613 public:
617  VOL_problem();
620  VOL_problem(const char *filename);
622  ~VOL_problem();
624 
630  int solve(VOL_user_hooks& hooks, const bool use_preset_dual = false);
632 
633 private:
637  double alpha_;
639  double lambda_;
640  // This union is here for padding (so that data members would be
641  // double-aligned on x86 CPU
642  union {
644  int iter_;
645  double __pad0;
646  };
648 
649 public:
650 
654  double value;
662 
668  int psize;
670  int dsize;
678 
679 public:
683  int iter() const { return iter_; }
685  double alpha() const { return alpha_; }
687  double lambda() const { return lambda_; }
689 
690 private:
694  void read_params(const char* filename);
695 
697  int initialize(const bool use_preset_dual);
698 
700  void print_info(const int iter,
701  const VOL_primal& primal, const VOL_primal& pstar,
702  const VOL_dual& dual);
703 
706  double readjust_target(const double oldtarget, const double lcost) const;
707 
715  double power_heur(const VOL_primal& primal, const VOL_primal& pstar,
716  const VOL_dual& dual) const;
718 };
719 
720 #endif
static T VolAbs(register const T x)
VOL_alpha_factor & operator=(const VOL_alpha_factor &)
double minimum_rel_ascent
terminate if the relative increase in lcost through ascent_check_invl steps is less than this ...
double lambdainit
initial value of lambda
double norm
double lambda() const
returns the value of lambda
The user hooks should be overridden by the user to provide the problem specific routines for the volu...
int solve(VOL_user_hooks &hooks, const bool use_preset_dual=false)
Solve the problem using the hooks.
int iter() const
returns the iteration number
VOL_dvector psol
final primal solution (OUTPUT)
VOL_primal & operator=(const VOL_primal &p)
VOL_primal(const int psize, const int dsize)
VOL_dvector()
Default constructor creates a vector of size 0.
VOL_dvector dsol
final dual solution (INPUT/OUTPUT)
#define VOL_TEST_SIZE(size)
void swap(VOL_dvector &w)
swaps the vector with w.
VOL_dvector dual_ub
upper bounds for the duals (if 0 length, then filled with +inf) (INPUT)
void cc(const double alpha, const VOL_primal &p)
This class holds every data for the Volume Algorithm and its solve method must be invoked to solve th...
double lcost
double granularity
terminate if best_ub - lcost &lt; granularity
int & operator[](const int i)
Return a reference to the i-th entry.
virtual ~VOL_user_hooks()
int ascent_check_invl
through how many iterations does the relative ascent have to reach a minimum
double ubinit
initial upper bound of the value of an integer solution
VOL_ivector()
Default constructor creates a vector of size 0.
void read_params(const char *filename)
Read in the parameters from the file filename.
~VOL_ivector()
The destructor deletes the data array.
VOL_indc & operator=(const VOL_indc &)
double power_heur(const VOL_primal &primal, const VOL_primal &pstar, const VOL_dual &dual) const
Here we decide the value of alpha1 to be used in the convex combination.
VOL_dvector & operator=(const VOL_dvector &w)
Copy w into the vector.
double alphafactor
when little progress is being done, we multiply alpha by alphafactor
enum VOL_swing::condition lastswing
int maxsgriters
maximum number of iterations
double & operator[](const int i)
Return a reference to the i-th entry.
double ascent(const VOL_dvector &v, const VOL_dvector &last_u) const
void clear()
Delete the content of the vector and replace it with a vector of length 0.
void allocate(const int s)
delete the current vector and allocate space for a vector of size s.
double factor(const VOL_parms &parm, const double lcost, const double alpha)
double asc
double value
final lagrangian value (OUTPUT)
double alphainit
initial value of alpha
void find_max_viol(const VOL_dvector &dual_lb, const VOL_dvector &dual_ub)
VOL_problem()
Default constructor.
double alpha_
value of alpha
VOL_ivector(const VOL_ivector &x)
Copy constructor makes a replica of x.
int dsize
length of dual solution (INPUT)
double operator[](const int i) const
Return the i-th entry.
void allocate(const int s)
delete the current vector and allocate space for a vector of size s.
int sz
The size of the vector.
int printflag
controls the level of printing.
int iter_
iteration number
double lambda_
value of lambda
VOL_swing & operator=(const VOL_swing &)
VOL_problem & operator=(const VOL_problem &)
VOL_ivector & operator=(const VOL_ivector &v)
Copy w into the vector.
double gap_abs_precision
accept if abs gap is less than this
VOL_primal(const VOL_primal &primal)
int alphaint
number of iterations before we check if alpha should be decreased
int size() const
Return the size of the vector.
double * v
The array holding the vector.
VOL_dvector(const VOL_dvector &x)
Copy constructor makes a replica of x.
int yellowtestinvl
how many consecutive yellow iterations are allowed before changing lambda
VOL_vh & operator=(const VOL_vh &)
VOL_indc(const VOL_indc &)
This class contains the parameters controlling the Volume Algorithm.
double hh
VOL_dvector dual_lb
lower bounds for the duals (if 0 length, then filled with -inf) (INPUT)
VOL_dvector u
void clear()
Delete the content of the vector and replace it with a vector of length 0.
int heurinvl
controls how often we run the primal heuristic
VOL_dual(const VOL_dual &dual)
int greentestinvl
how many consecutive green iterations are allowed before changing lambda
int initialize(const bool use_preset_dual)
initializes duals, bounds for the duals, alpha, lambda
double gap_rel_precision
accept if rel gap is less than this
void set_default_parm()
int psize
length of primal solution (INPUT)
double vh
int * v
The array holding the vector.
VOL_dvector v
void swap(VOL_ivector &w)
swaps the vector with w.
int printinvl
controls how often do we print
~VOL_dvector()
The destructor deletes the data array.
char * temp_dualfile
name of file for saving dual solution
~VOL_problem()
Destruct the object.
double alpha() const
returns the value of alpha
const double COIN_DBL_MAX
Definition: CoinFinite.hpp:18
double alphamin
minimum value for alpha
void compute_xrc(const VOL_dvector &pstarx, const VOL_dvector &primalx, const VOL_dvector &rc)
VOL_dvector x
vector of ints.
void cc(const double gamma, const VOL_dvector &w)
Convex combination.
#define VOL_TEST_INDEX(i, size)
virtual int compute_rc(const VOL_dvector &u, VOL_dvector &rc)=0
compute reduced costs
VOL_ivector(const int s)
Construct a vector of size s.
int ascent_first_check
when to check for sufficient relative ascent the first time
VOL_vh(const VOL_vh &)
VOL_parms parm
The parameters controlling the Volume Algorithm (INPUT)
vector of doubles.
int redtestinvl
how many consecutive red iterations are allowed before changing lambda
void step(const double target, const double lambda, const VOL_dvector &dual_lb, const VOL_dvector &dual_ub, const VOL_dvector &v)
int size() const
Return the size of the vector.
int sz
The size of the vector.
double lfactor(const VOL_parms &parm, const double lambda, const int iter)
VOL_dvector(const int s)
Construct a vector of size s.
virtual int solve_subproblem(const VOL_dvector &dual, const VOL_dvector &rc, double &lcost, VOL_dvector &x, VOL_dvector &v, double &pcost)=0
Solve the subproblem for the subgradient step.
VOL_dvector viol
violations (b-Ax) for the relaxed constraints
VOL_dual(const int dsize)
int operator[](const int i) const
Return the i-th entry.
virtual int heuristics(const VOL_problem &p, const VOL_dvector &x, double &heur_val)=0
Starting from the primal vector x, run a heuristic to produce an integer solution.
double primal_abs_precision
accept if max abs viol is less than this
void cond(const VOL_dual &dual, const double lcost, const double ascent, const int iter)
static T VolMax(register const T x, register const T y)
void print_info(const int iter, const VOL_primal &primal, const VOL_primal &pstar, const VOL_dual &dual)
print volume info every parm.printinvl iterations
VOL_dual & operator=(const VOL_dual &p)
double readjust_target(const double oldtarget, const double lcost) const
Checks if lcost is close to the target, if so it increases the target.