00001
00002
00003 #ifndef __VOLUME_HPP__
00004 #define __VOLUME_HPP__
00005
00006 #include <cfloat>
00007 #include <algorithm>
00008 #include <cstdio>
00009 #include <cmath>
00010
00011 #ifndef VOL_DEBUG
00012
00013 #define VOL_DEBUG 0
00014 #endif
00015
00016 template <class T> static inline T
00017 VolMax(register const T x, register const T y) {
00018 return ((x) > (y)) ? (x) : (y);
00019 }
00020
00021 template <class T> static inline T
00022 VolAbs(register const T x) {
00023 return ((x) > 0) ? (x) : -(x);
00024 }
00025
00026
00027
00028 #if defined(VOL_DEBUG) && (VOL_DEBUG != 0)
00029 #define VOL_TEST_INDEX(i, size) \
00030 { \
00031 if ((i) < 0 || (i) >= (size)) { \
00032 printf("bad VOL_?vector index\n"); \
00033 abort(); \
00034 } \
00035 }
00036 #define VOL_TEST_SIZE(size) \
00037 { \
00038 if (s <= 0) { \
00039 printf("bad VOL_?vector size\n"); \
00040 abort(); \
00041 } \
00042 }
00043 #else
00044 #define VOL_TEST_INDEX(i, size)
00045 #define VOL_TEST_SIZE(size)
00046 #endif
00047
00048
00049
00050 class VOL_dvector;
00051 class VOL_ivector;
00052 class VOL_primal;
00053 class VOL_dual;
00054 class VOL_swing;
00055 class VOL_alpha_factor;
00056 class VOL_vh;
00057 class VOL_indc;
00058 class VOL_user_hooks;
00059 class VOL_problem;
00060
00061
00062
00066 struct VOL_parms {
00068 double lambdainit;
00070 double alphainit;
00072 double alphamin;
00074 double alphafactor;
00075
00077 double ubinit;
00078
00080 double primal_abs_precision;
00082 double gap_abs_precision;
00084 double gap_rel_precision;
00086 double granularity;
00087
00090 double minimum_rel_ascent;
00092 int ascent_first_check;
00095 int ascent_check_invl;
00096
00098 int maxsgriters;
00099
00110 int printflag;
00112 int printinvl;
00114 int heurinvl;
00115
00118 int greentestinvl;
00121 int yellowtestinvl;
00124 int redtestinvl;
00125
00127 int alphaint;
00128
00130 char* temp_dualfile;
00131 };
00132
00133
00134
00143 class VOL_dvector {
00144 public:
00146 double* v;
00148 int sz;
00149
00150 public:
00152 VOL_dvector(const int s) {
00153 VOL_TEST_SIZE(s);
00154 v = new double[sz = s];
00155 }
00157 VOL_dvector() : v(0), sz(0) {}
00159 VOL_dvector(const VOL_dvector& x) : v(0), sz(0) {
00160 sz = x.sz;
00161 if (sz > 0) {
00162 v = new double[sz];
00163 std::copy(x.v, x.v + sz, v);
00164 }
00165 }
00167 ~VOL_dvector() { delete[] v; }
00168
00170 inline int size() const {return sz;}
00171
00173 inline double& operator[](const int i) {
00174 VOL_TEST_INDEX(i, sz);
00175 return v[i];
00176 }
00177
00179 inline double operator[](const int i) const {
00180 VOL_TEST_INDEX(i, sz);
00181 return v[i];
00182 }
00183
00186 inline void clear() {
00187 delete[] v;
00188 v = 0;
00189 sz = 0;
00190 }
00193 inline void cc(const double gamma, const VOL_dvector& w) {
00194 if (sz != w.sz) {
00195 printf("bad VOL_dvector sizes\n");
00196 abort();
00197 }
00198 double * p_v = v - 1;
00199 const double * p_w = w.v - 1;
00200 const double * const p_e = v + sz;
00201 const double one_gamma = 1.0 - gamma;
00202 while ( ++p_v != p_e ){
00203 *p_v = one_gamma * (*p_v) + gamma * (*++p_w);
00204 }
00205 }
00206
00209 inline void allocate(const int s) {
00210 VOL_TEST_SIZE(s);
00211 delete[] v;
00212 v = new double[sz = s];
00213 }
00214
00216 inline void swap(VOL_dvector& w) {
00217 std::swap(v, w.v);
00218 std::swap(sz, w.sz);
00219 }
00220
00222 VOL_dvector& operator=(const VOL_dvector& w);
00224 VOL_dvector& operator=(const double w);
00225 };
00226
00227
00237 class VOL_ivector {
00238 public:
00240 int* v;
00242 int sz;
00243 public:
00245 VOL_ivector(const int s) {
00246 VOL_TEST_SIZE(s);
00247 v = new int[sz = s];
00248 }
00250 VOL_ivector() : v(0), sz(0) {}
00252 VOL_ivector(const VOL_ivector& x) {
00253 sz = x.sz;
00254 if (sz > 0) {
00255 v = new int[sz];
00256 std::copy(x.v, x.v + sz, v);
00257 }
00258 }
00260 ~VOL_ivector(){
00261 delete [] v;
00262 }
00263
00265 inline int size() const { return sz; }
00267 inline int& operator[](const int i) {
00268 VOL_TEST_INDEX(i, sz);
00269 return v[i];
00270 }
00271
00273 inline int operator[](const int i) const {
00274 VOL_TEST_INDEX(i, sz);
00275 return v[i];
00276 }
00277
00280 inline void clear() {
00281 delete[] v;
00282 v = 0;
00283 sz = 0;
00284 }
00285
00288 inline void allocate(const int s) {
00289 VOL_TEST_SIZE(s);
00290 delete[] v;
00291 v = new int[sz = s];
00292 }
00293
00295 inline void swap(VOL_ivector& w) {
00296 std::swap(v, w.v);
00297 std::swap(sz, w.sz);
00298 }
00299
00301 VOL_ivector& operator=(const VOL_ivector& v);
00303 VOL_ivector& operator=(const int w);
00304 };
00305
00306
00307
00308 class VOL_primal {
00309 public:
00310
00311 double value;
00312
00313 double viol;
00314
00315 VOL_dvector x;
00316
00317 VOL_dvector v;
00318
00319 VOL_primal(const int psize, const int dsize) : x(psize), v(dsize) {}
00320 VOL_primal(const VOL_primal& primal) :
00321 value(primal.value), viol(primal.viol), x(primal.x), v(primal.v) {}
00322 ~VOL_primal() {}
00323 inline VOL_primal& operator=(const VOL_primal& p) {
00324 if (this == &p)
00325 return *this;
00326 value = p.value;
00327 viol = p.viol;
00328 x = p.x;
00329 v = p.v;
00330 return *this;
00331 }
00332
00333
00334
00335
00336
00337 inline void cc(const double alpha, const VOL_primal& p) {
00338 value = alpha * p.value + (1.0 - alpha) * value;
00339 x.cc(alpha, p.x);
00340 v.cc(alpha, p.v);
00341 }
00342
00343 void find_max_viol(const VOL_dvector& dual_lb,
00344 const VOL_dvector& dual_ub);
00345 };
00346
00347
00348
00349 class VOL_dual {
00350 public:
00351
00352 double lcost;
00353
00354 double xrc;
00355
00356
00357 VOL_dvector u;
00358
00359 VOL_dual(const int dsize) : u(dsize) { u = 0.0;}
00360 VOL_dual(const VOL_dual& dual) :
00361 lcost(dual.lcost), xrc(dual.xrc), u(dual.u) {}
00362 ~VOL_dual() {}
00363 inline VOL_dual& operator=(const VOL_dual& p) {
00364 if (this == &p)
00365 return *this;
00366 lcost = p.lcost;
00367 xrc = p.xrc;
00368 u = p.u;
00369 return *this;
00370 }
00371
00372 void step(const double target, const double lambda,
00373 const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
00374 const VOL_dvector& v);
00375 double ascent(const VOL_dvector& v, const VOL_dvector& last_u) const;
00376 void compute_xrc(const VOL_dvector& pstarx, const VOL_dvector& primalx,
00377 const VOL_dvector& rc);
00378
00379 };
00380
00381
00382
00383
00384
00385 class VOL_swing {
00386 private:
00387 VOL_swing(const VOL_swing&);
00388 VOL_swing& operator=(const VOL_swing&);
00389 public:
00390 enum condition {green, yellow, red} lastswing;
00391 int lastgreeniter, lastyellowiter, lastrediter;
00392 int ngs, nrs, nys;
00393 int rd;
00394
00395 VOL_swing() {
00396 lastgreeniter = lastyellowiter = lastrediter = 0;
00397 ngs = nrs = nys = 0;
00398 }
00399 ~VOL_swing(){}
00400
00401 inline void cond(const VOL_dual& dual,
00402 const double lcost, const double ascent, const int iter) {
00403 double eps = 1.e-3;
00404
00405 if (ascent > 0.0 && lcost > dual.lcost + eps) {
00406 lastswing = green;
00407 lastgreeniter = iter;
00408 ++ngs;
00409 rd = 0;
00410 } else {
00411 if (ascent <= 0 && lcost > dual.lcost) {
00412 lastswing = yellow;
00413 lastyellowiter = iter;
00414 ++nys;
00415 rd = 0;
00416 } else {
00417 lastswing = red;
00418 lastrediter = iter;
00419 ++nrs;
00420 rd = 1;
00421 }
00422 }
00423 }
00424
00425 inline double
00426 lfactor(const VOL_parms& parm, const double lambda, const int iter) {
00427 double lambdafactor = 1.0;
00428 double eps = 5.e-4;
00429 int cons;
00430
00431 switch (lastswing) {
00432 case green:
00433 cons = iter - VolMax(lastyellowiter, lastrediter);
00434 if (parm.printflag & 4)
00435 printf(" G: Consecutive Gs = %3d\n\n", cons);
00436 if (cons >= parm.greentestinvl && lambda < 2.0) {
00437 lastgreeniter = lastyellowiter = lastrediter = iter;
00438 lambdafactor = 2.0;
00439 if (parm.printflag & 2)
00440 printf("\n ---- increasing lamda to %g ----\n\n",
00441 lambda * lambdafactor);
00442 }
00443 break;
00444
00445 case yellow:
00446 cons = iter - VolMax(lastgreeniter, lastrediter);
00447 if (parm.printflag & 4)
00448 printf(" Y: Consecutive Ys = %3d\n\n", cons);
00449 if (cons >= parm.yellowtestinvl) {
00450 lastgreeniter = lastyellowiter = lastrediter = iter;
00451 lambdafactor = 1.1;
00452 if (parm.printflag & 2)
00453 printf("\n **** increasing lamda to %g *****\n\n",
00454 lambda * lambdafactor);
00455 }
00456 break;
00457
00458 case red:
00459 cons = iter - VolMax(lastgreeniter, lastyellowiter);
00460 if (parm.printflag & 4)
00461 printf(" R: Consecutive Rs = %3d\n\n", cons);
00462 if (cons >= parm.redtestinvl && lambda > eps) {
00463 lastgreeniter = lastyellowiter = lastrediter = iter;
00464 lambdafactor = 0.67;
00465 if (parm.printflag & 2)
00466 printf("\n **** decreasing lamda to %g *****\n\n",
00467 lambda * lambdafactor);
00468 }
00469 break;
00470 }
00471 return lambdafactor;
00472 }
00473
00474 inline void
00475 print() {
00476 printf("**** G= %i, Y= %i, R= %i ****\n", ngs, nys, nrs);
00477 ngs = nrs = nys = 0;
00478 }
00479 };
00480
00481
00482
00483
00484 class VOL_alpha_factor {
00485 private:
00486 VOL_alpha_factor(const VOL_alpha_factor&);
00487 VOL_alpha_factor& operator=(const VOL_alpha_factor&);
00488 public:
00489 double lastvalue;
00490
00491 VOL_alpha_factor() {lastvalue = -DBL_MAX;}
00492 ~VOL_alpha_factor() {}
00493
00494 inline double factor(const VOL_parms& parm, const double lcost,
00495 const double alpha) {
00496 if (alpha < parm.alphamin)
00497 return 1.0;
00498 register const double ll = VolAbs(lcost);
00499 const double x = ll > 10 ? (lcost-lastvalue)/ll : (lcost-lastvalue);
00500 lastvalue = lcost;
00501 return (x <= 0.01) ? parm.alphafactor : 1.0;
00502 }
00503 };
00504
00505
00506
00507
00508
00509
00510 class VOL_vh {
00511 private:
00512 VOL_vh(const VOL_vh&);
00513 VOL_vh& operator=(const VOL_vh&);
00514 public:
00515 double hh;
00516 double norm;
00517 double vh;
00518 double asc;
00519
00520 VOL_vh(const double alpha,
00521 const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
00522 const VOL_dvector& v, const VOL_dvector& vstar,
00523 const VOL_dvector& u);
00524 ~VOL_vh(){}
00525 };
00526
00527
00528
00529
00530
00531
00532
00533 class VOL_indc {
00534 private:
00535 VOL_indc(const VOL_indc&);
00536 VOL_indc& operator=(const VOL_indc&);
00537 public:
00538 double v2;
00539 double vu;
00540 double vabs;
00541 double asc;
00542
00543 public:
00544 VOL_indc(const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
00545 const VOL_primal& primal, const VOL_primal& pstar,
00546 const VOL_dual& dual);
00547 ~VOL_indc() {}
00548 };
00549
00550
00551
00558 class VOL_user_hooks {
00559 public:
00560 virtual ~VOL_user_hooks() {}
00561 public:
00562
00567 virtual int compute_rc(const VOL_dvector& u, VOL_dvector& rc) = 0;
00568
00577 virtual int solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
00578 double& lcost, VOL_dvector& x, VOL_dvector& v,
00579 double& pcost) = 0;
00586 virtual int heuristics(const VOL_problem& p,
00587 const VOL_dvector& x, double& heur_val) = 0;
00588 };
00589
00590
00591
00600 class VOL_problem {
00601 private:
00602 VOL_problem(const VOL_problem&);
00603 VOL_problem& operator=(const VOL_problem&);
00604 void set_default_parm();
00605
00606 public:
00610 VOL_problem();
00613 VOL_problem(const char *filename);
00615 ~VOL_problem();
00617
00623 int solve(VOL_user_hooks& hooks, const bool use_preset_dual = false);
00625
00626 private:
00630 double alpha_;
00632 double lambda_;
00633
00634
00635 union {
00637 int iter_;
00638 double __pad0;
00639 };
00641
00642 public:
00643
00647 double value;
00649 VOL_dvector dsol;
00651 VOL_dvector psol;
00653 VOL_dvector viol;
00655
00659 VOL_parms parm;
00661 int psize;
00663 int dsize;
00666 VOL_dvector dual_lb;
00669 VOL_dvector dual_ub;
00671
00672 public:
00676 int iter() const { return iter_; }
00678 double alpha() const { return alpha_; }
00680 double lambda() const { return lambda_; }
00682
00683 private:
00687 void read_params(const char* filename);
00688
00690 int initialize(const bool use_preset_dual);
00691
00693 void print_info(const int iter,
00694 const VOL_primal& primal, const VOL_primal& pstar,
00695 const VOL_dual& dual);
00696
00699 double readjust_target(const double oldtarget, const double lcost) const;
00700
00708 double power_heur(const VOL_primal& primal, const VOL_primal& pstar,
00709 const VOL_dual& dual) const;
00711 };
00712
00713 #endif