VolVolume.hpp

Go to the documentation of this file.
00001 // $Id: VolVolume.hpp 375 2013-11-22 15:46:16Z tkr $
00002 // Copyright (C) 2000, International Business Machines
00003 // Corporation and others.  All Rights Reserved.
00004 // This code is licensed under the terms of the Eclipse Public License (EPL).
00005 
00006 #ifndef __VOLUME_HPP__
00007 #define __VOLUME_HPP__
00008 
00009 #include <cfloat>
00010 #include <algorithm>
00011 #include <cstdio>
00012 #include <cmath>
00013 #include <cstring>
00014 #include <cstdlib>
00015 
00016 #ifndef VOL_DEBUG
00017 // When VOL_DEBUG is 1, we check vector indices
00018 #define VOL_DEBUG 0
00019 #endif
00020 
00021 template <class T> static inline T
00022 VolMax(register const T x, register const T y) {
00023    return ((x) > (y)) ? (x) : (y);
00024 }
00025 
00026 template <class T> static inline T
00027 VolAbs(register const T x) {
00028    return ((x) > 0) ? (x) : -(x);
00029 }
00030 
00031 //############################################################################
00032 
00033 #if defined(VOL_DEBUG) && (VOL_DEBUG != 0)
00034 #define VOL_TEST_INDEX(i, size)                 \
00035 {                                               \
00036    if ((i) < 0 || (i) >= (size)) {              \
00037       printf("bad VOL_?vector index\n");        \
00038       abort();                                  \
00039    }                                            \
00040 }
00041 #define VOL_TEST_SIZE(size)                     \
00042 {                                               \
00043    if (s <= 0) {                                \
00044       printf("bad VOL_?vector size\n");         \
00045       abort();                                  \
00046    }                                            \
00047 }
00048 #else
00049 #define VOL_TEST_INDEX(i, size)
00050 #define VOL_TEST_SIZE(size)
00051 #endif
00052       
00053 //############################################################################
00054 
00055 class VOL_dvector;
00056 class VOL_ivector;
00057 class VOL_primal;
00058 class VOL_dual;
00059 class VOL_swing;
00060 class VOL_alpha_factor;
00061 class VOL_vh;
00062 class VOL_indc;
00063 class VOL_user_hooks;
00064 class VOL_problem;
00065 
00066 //############################################################################
00067 
00071 struct VOL_parms {
00073    double lambdainit; 
00075    double alphainit; 
00077    double alphamin;
00079    double alphafactor;
00080 
00082    double ubinit;
00083 
00085    double primal_abs_precision;
00087    double gap_abs_precision; 
00089    double gap_rel_precision;  
00091    double granularity;
00092 
00095    double minimum_rel_ascent;
00097    int    ascent_first_check;
00100    int    ascent_check_invl;
00101    
00103    int    maxsgriters; 
00104 
00115    int    printflag; 
00117    int    printinvl; 
00119    int    heurinvl; 
00120 
00123    int greentestinvl; 
00126    int yellowtestinvl; 
00129    int redtestinvl;
00130 
00132    int    alphaint; 
00133 
00135    char* temp_dualfile;
00136 };
00137 
00138 //############################################################################
00139 
00148 class VOL_dvector {
00149 public:
00151    double* v;
00153    int sz;
00154 
00155 public:
00157    VOL_dvector(const int s) {
00158       VOL_TEST_SIZE(s);
00159       v = new double[sz = s];
00160    }
00162    VOL_dvector() : v(0), sz(0) {}
00164    VOL_dvector(const VOL_dvector& x) : v(0), sz(0) {
00165       sz = x.sz;
00166       if (sz > 0) {
00167         v = new double[sz];
00168         std::memcpy(v, x.v, sz * sizeof(double));
00169       }
00170    }
00172    ~VOL_dvector() { delete[] v; }
00173 
00175    inline int size() const {return sz;}
00176 
00178    inline double& operator[](const int i) {
00179       VOL_TEST_INDEX(i, sz);
00180       return v[i];
00181    }
00182 
00184    inline double operator[](const int i) const {
00185       VOL_TEST_INDEX(i, sz);
00186       return v[i];
00187    }
00188 
00191    inline void clear() {
00192       delete[] v;
00193       v = 0;
00194       sz = 0;
00195    }
00198    inline void cc(const double gamma, const VOL_dvector& w) {
00199       if (sz != w.sz) {
00200          printf("bad VOL_dvector sizes\n");
00201          abort();
00202       }
00203       double * p_v = v - 1;
00204       const double * p_w = w.v - 1;
00205       const double * const p_e = v + sz;
00206       const double one_gamma = 1.0 - gamma;
00207       while ( ++p_v != p_e ){
00208          *p_v = one_gamma * (*p_v) + gamma * (*++p_w);
00209       }
00210    }
00211 
00214    inline void allocate(const int s) {
00215       VOL_TEST_SIZE(s);
00216       delete[] v;                 
00217       v = new double[sz = s];
00218    }
00219 
00221    inline void swap(VOL_dvector& w) {
00222       std::swap(v, w.v);
00223       std::swap(sz, w.sz);
00224    }
00225 
00227    VOL_dvector& operator=(const VOL_dvector& w);
00229    VOL_dvector& operator=(const double w);
00230 };
00231 
00232 //-----------------------------------------------------------------------------
00242 class VOL_ivector {
00243 public:
00245    int* v;
00247    int sz;
00248 public:
00250    VOL_ivector(const int s) {
00251       VOL_TEST_SIZE(s);
00252       v = new int[sz = s];
00253    }
00255    VOL_ivector() : v(0), sz(0) {}
00257    VOL_ivector(const VOL_ivector& x) {
00258       sz = x.sz;
00259       if (sz > 0) {
00260          v = new int[sz];
00261          std::memcpy(v, x.v, sz * sizeof(int));
00262       }
00263    }
00265    ~VOL_ivector(){
00266       delete [] v;
00267    }
00268 
00270    inline int size() const { return sz; }
00272    inline int& operator[](const int i) {
00273       VOL_TEST_INDEX(i, sz);
00274       return v[i];
00275    }
00276 
00278    inline int operator[](const int i) const {
00279       VOL_TEST_INDEX(i, sz);
00280       return v[i];
00281    }
00282 
00285    inline void clear() {
00286       delete[] v;
00287       v = 0;
00288       sz = 0;
00289    }
00290 
00293    inline void allocate(const int s) {
00294       VOL_TEST_SIZE(s);
00295       delete[] v;
00296       v = new int[sz = s];
00297    }
00298 
00300    inline void swap(VOL_ivector& w) {
00301       std::swap(v, w.v);
00302       std::swap(sz, w.sz);
00303    }
00304 
00306    VOL_ivector& operator=(const VOL_ivector& v);      
00308    VOL_ivector& operator=(const int w);
00309 };
00310 
00311 //############################################################################
00312 // A class describing a primal solution. This class is used only internally 
00313 class VOL_primal {
00314 public: 
00315    // objective value of this primal solution 
00316    double value;
00317    // the largest of the v[i]'s
00318    double viol;  
00319    // primal solution  
00320    VOL_dvector x;
00321    // v=b-Ax, for the relaxed constraints
00322    VOL_dvector v; 
00323 
00324    VOL_primal(const int psize, const int dsize) : x(psize), v(dsize) {}
00325    VOL_primal(const VOL_primal& primal) :
00326       value(primal.value), viol(primal.viol), x(primal.x), v(primal.v) {}
00327    ~VOL_primal() {}
00328    inline VOL_primal& operator=(const VOL_primal& p) {
00329       if (this == &p) 
00330          return *this;
00331       value = p.value;
00332       viol = p.viol;
00333       x = p.x;
00334       v = p.v;
00335       return *this;
00336    }
00337 
00338    // convex combination. data members in this will be overwritten
00339    // convex combination between two primal solutions
00340    // x <-- alpha x + (1 - alpha) p.x
00341    // v <-- alpha v + (1 - alpha) p.v
00342    inline void cc(const double alpha, const VOL_primal& p) {
00343       value = alpha * p.value + (1.0 - alpha) * value;
00344       x.cc(alpha, p.x);
00345       v.cc(alpha, p.v);
00346    }
00347    // find maximum of v[i]
00348    void find_max_viol(const VOL_dvector& dual_lb, 
00349                       const VOL_dvector& dual_ub);
00350 };
00351 
00352 //-----------------------------------------------------------------------------
00353 // A class describing a dual solution. This class is used only internally 
00354 class VOL_dual {
00355 public:
00356    // lagrangian value
00357    double lcost; 
00358    // reduced costs * (pstar-primal)
00359    double xrc;
00360    // this information is only printed
00361    // dual vector
00362    VOL_dvector u; 
00363 
00364    VOL_dual(const int dsize) : u(dsize) { u = 0.0;}
00365    VOL_dual(const VOL_dual& dual) :
00366       lcost(dual.lcost), xrc(dual.xrc), u(dual.u) {}
00367    ~VOL_dual() {}
00368    inline VOL_dual& operator=(const VOL_dual& p) {
00369       if (this == &p) 
00370          return *this;
00371       lcost = p.lcost;
00372       xrc = p.xrc;
00373       u = p.u;
00374       return *this;
00375    }
00376    // dual step
00377    void   step(const double target, const double lambda,
00378                const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
00379                const VOL_dvector& v);
00380    double ascent(const VOL_dvector& v, const VOL_dvector& last_u) const;
00381    void   compute_xrc(const VOL_dvector& pstarx, const VOL_dvector& primalx,
00382                       const VOL_dvector& rc);
00383 
00384 };
00385 
00386 
00387 //############################################################################
00388 /* here we check whether an iteration is green, yellow or red. Also according
00389    to this information we decide whether lambda should be changed */
00390 class VOL_swing {
00391 private:
00392    VOL_swing(const VOL_swing&);
00393    VOL_swing& operator=(const VOL_swing&);
00394 public:
00395    enum condition {green, yellow, red} lastswing;
00396    int lastgreeniter, lastyellowiter, lastrediter;
00397    int ngs, nrs, nys;
00398    int rd;
00399    
00400    VOL_swing() {
00401       lastgreeniter = lastyellowiter = lastrediter = 0;
00402       ngs = nrs = nys = 0;
00403    }
00404    ~VOL_swing(){}
00405 
00406    inline void cond(const VOL_dual& dual, 
00407                     const double lcost, const double ascent, const int iter) {
00408       double eps = 1.e-3;
00409 
00410       if (ascent > 0.0  &&  lcost > dual.lcost + eps) {
00411          lastswing = green;
00412          lastgreeniter = iter;
00413          ++ngs;
00414          rd = 0;
00415       } else { 
00416          if (ascent <= 0  &&  lcost > dual.lcost) {
00417             lastswing = yellow;
00418             lastyellowiter = iter;
00419             ++nys;
00420             rd = 0;
00421          } else {
00422             lastswing = red;
00423             lastrediter = iter;
00424             ++nrs;
00425             rd = 1;
00426          }
00427       }
00428    }
00429 
00430    inline double
00431    lfactor(const VOL_parms& parm, const double lambda, const int iter) {
00432       double lambdafactor = 1.0;
00433       double eps = 5.e-4;
00434       int cons;
00435 
00436       switch (lastswing) {
00437        case green:
00438          cons = iter - VolMax(lastyellowiter, lastrediter);
00439          if (parm.printflag & 4)
00440             printf("      G: Consecutive Gs = %3d\n\n", cons);
00441          if (cons >= parm.greentestinvl && lambda < 2.0) {
00442             lastgreeniter = lastyellowiter = lastrediter = iter;
00443             lambdafactor = 2.0;
00444             if (parm.printflag & 2)
00445                printf("\n ---- increasing lamda to %g ----\n\n",
00446                       lambda * lambdafactor); 
00447          }
00448          break;
00449       
00450        case yellow:
00451          cons = iter - VolMax(lastgreeniter, lastrediter);
00452          if (parm.printflag & 4)
00453             printf("      Y: Consecutive Ys = %3d\n\n", cons);
00454          if (cons >= parm.yellowtestinvl) {
00455             lastgreeniter = lastyellowiter = lastrediter = iter;
00456             lambdafactor = 1.1;
00457             if (parm.printflag & 2)
00458                printf("\n **** increasing lamda to %g *****\n\n",
00459                       lambda * lambdafactor);
00460          }
00461          break;
00462       
00463        case red:
00464          cons = iter - VolMax(lastgreeniter, lastyellowiter);
00465          if (parm.printflag & 4)
00466             printf("      R: Consecutive Rs = %3d\n\n", cons);
00467          if (cons >= parm.redtestinvl && lambda > eps) {
00468             lastgreeniter = lastyellowiter = lastrediter = iter;
00469             lambdafactor = 0.67;
00470             if (parm.printflag & 2)
00471                printf("\n **** decreasing lamda to %g *****\n\n",
00472                       lambda * lambdafactor);
00473          } 
00474          break;
00475       }
00476       return lambdafactor;
00477    }
00478 
00479    inline void
00480    print() {
00481       printf("**** G= %i, Y= %i, R= %i ****\n", ngs, nys, nrs);
00482       ngs = nrs = nys = 0;  
00483    }
00484 };
00485 
00486 //############################################################################
00487 /* alpha should be decreased if after some number of iterations the objective
00488    has increased less that 1% */
00489 class VOL_alpha_factor {
00490 private:
00491    VOL_alpha_factor(const VOL_alpha_factor&);
00492    VOL_alpha_factor& operator=(const VOL_alpha_factor&);
00493 public:
00494    double lastvalue;
00495 
00496    VOL_alpha_factor() {lastvalue = -DBL_MAX;}
00497    ~VOL_alpha_factor() {}
00498 
00499    inline double factor(const VOL_parms& parm, const double lcost,
00500                         const double alpha) {
00501       if (alpha < parm.alphamin)
00502          return 1.0;
00503       register const double ll = VolAbs(lcost);
00504       const double x = ll > 10 ? (lcost-lastvalue)/ll : (lcost-lastvalue);
00505       lastvalue = lcost;
00506       return (x <= 0.01) ? parm.alphafactor : 1.0;
00507    }
00508 };
00509 
00510 //############################################################################
00511 /* here we compute the norm of the conjugate direction -hh-, the norm of the
00512    subgradient -norm-, the inner product between the subgradient and the 
00513    last conjugate direction -vh-, and the inner product between the new
00514    conjugate direction and the subgradient */
00515 class VOL_vh {
00516 private:
00517    VOL_vh(const VOL_vh&);
00518    VOL_vh& operator=(const VOL_vh&);
00519 public:
00520    double hh;
00521    double norm;
00522    double vh;
00523    double asc;
00524 
00525    VOL_vh(const double alpha,
00526           const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
00527           const VOL_dvector& v, const VOL_dvector& vstar,
00528           const VOL_dvector& u);
00529    ~VOL_vh(){}
00530 };
00531 
00532 //############################################################################
00533 /* here we compute different parameter to be printed. v2 is the square of 
00534    the norm of the subgradient. vu is the inner product between the dual
00535    variables and the subgradient. vabs is the maximum absolute value of
00536    the violations of pstar. asc is the inner product between the conjugate
00537    direction and the subgradient */
00538 class VOL_indc {
00539 private:
00540    VOL_indc(const VOL_indc&);
00541    VOL_indc& operator=(const VOL_indc&);
00542 public:
00543    double v2;
00544    double vu;
00545    double vabs;
00546    double asc;
00547 
00548 public:
00549    VOL_indc(const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
00550             const VOL_primal& primal, const VOL_primal& pstar,
00551             const VOL_dual& dual);
00552    ~VOL_indc() {}
00553 };
00554 
00555 //#############################################################################
00556 
00563 class VOL_user_hooks {
00564 public:
00565    virtual ~VOL_user_hooks() {}
00566 public:
00567    // for all hooks: return value of -1 means that volume should quit
00572    virtual int compute_rc(const VOL_dvector& u, VOL_dvector& rc) = 0;
00573 
00582    virtual int solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
00583                                 double& lcost, VOL_dvector& x, VOL_dvector& v,
00584                                 double& pcost) = 0;
00591    virtual int heuristics(const VOL_problem& p, 
00592                           const VOL_dvector& x, double& heur_val) = 0;
00593 };
00594 
00595 //#############################################################################
00596 
00605 class VOL_problem {
00606 private:
00607    VOL_problem(const VOL_problem&);
00608    VOL_problem& operator=(const VOL_problem&);
00609    void set_default_parm();
00610    // ############ INPUT fields ########################
00611 public: 
00615    VOL_problem();
00618    VOL_problem(const char *filename);
00620    ~VOL_problem();
00622 
00628    int solve(VOL_user_hooks& hooks, const bool use_preset_dual = false);
00630 
00631 private: 
00635    double alpha_; 
00637    double lambda_;
00638    // This union is here for padding (so that data members would be
00639    // double-aligned on x86 CPU
00640    union {
00642       int iter_;
00643       double __pad0;
00644    };
00646 
00647 public:
00648   
00652    double value;
00654    VOL_dvector dsol;
00656    VOL_dvector psol;
00658    VOL_dvector viol;
00660 
00664    VOL_parms parm;
00666    int psize;        
00668    int dsize;      
00671    VOL_dvector dual_lb;
00674    VOL_dvector dual_ub;
00676 
00677 public:
00681    int    iter() const { return iter_; }
00683    double alpha() const { return alpha_; }
00685    double lambda() const { return lambda_; }
00687 
00688 private:
00692    void read_params(const char* filename);
00693 
00695    int initialize(const bool use_preset_dual);
00696 
00698    void print_info(const int iter,
00699                    const VOL_primal& primal, const VOL_primal& pstar,
00700                    const VOL_dual& dual);
00701 
00704    double readjust_target(const double oldtarget, const double lcost) const;
00705 
00713    double power_heur(const VOL_primal& primal, const VOL_primal& pstar,
00714                      const VOL_dual& dual) const;
00716 };
00717 
00718 #endif
 All Classes Files Functions Variables Enumerations Enumerator Friends Defines

Generated on 31 Mar 2015 for Vol by  doxygen 1.6.1