VolVolume.hpp

Go to the documentation of this file.
00001 // $Id: VolVolume.hpp 321 2011-07-31 16:05:18Z stefan $
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 
00015 #ifndef VOL_DEBUG
00016 // When VOL_DEBUG is 1, we check vector indices
00017 #define VOL_DEBUG 0
00018 #endif
00019 
00020 template <class T> static inline T
00021 VolMax(register const T x, register const T y) {
00022    return ((x) > (y)) ? (x) : (y);
00023 }
00024 
00025 template <class T> static inline T
00026 VolAbs(register const T x) {
00027    return ((x) > 0) ? (x) : -(x);
00028 }
00029 
00030 //############################################################################
00031 
00032 #if defined(VOL_DEBUG) && (VOL_DEBUG != 0)
00033 #define VOL_TEST_INDEX(i, size)                 \
00034 {                                               \
00035    if ((i) < 0 || (i) >= (size)) {              \
00036       printf("bad VOL_?vector index\n");        \
00037       abort();                                  \
00038    }                                            \
00039 }
00040 #define VOL_TEST_SIZE(size)                     \
00041 {                                               \
00042    if (s <= 0) {                                \
00043       printf("bad VOL_?vector size\n");         \
00044       abort();                                  \
00045    }                                            \
00046 }
00047 #else
00048 #define VOL_TEST_INDEX(i, size)
00049 #define VOL_TEST_SIZE(size)
00050 #endif
00051       
00052 //############################################################################
00053 
00054 class VOL_dvector;
00055 class VOL_ivector;
00056 class VOL_primal;
00057 class VOL_dual;
00058 class VOL_swing;
00059 class VOL_alpha_factor;
00060 class VOL_vh;
00061 class VOL_indc;
00062 class VOL_user_hooks;
00063 class VOL_problem;
00064 
00065 //############################################################################
00066 
00070 struct VOL_parms {
00072    double lambdainit; 
00074    double alphainit; 
00076    double alphamin;
00078    double alphafactor;
00079 
00081    double ubinit;
00082 
00084    double primal_abs_precision;
00086    double gap_abs_precision; 
00088    double gap_rel_precision;  
00090    double granularity;
00091 
00094    double minimum_rel_ascent;
00096    int    ascent_first_check;
00099    int    ascent_check_invl;
00100    
00102    int    maxsgriters; 
00103 
00114    int    printflag; 
00116    int    printinvl; 
00118    int    heurinvl; 
00119 
00122    int greentestinvl; 
00125    int yellowtestinvl; 
00128    int redtestinvl;
00129 
00131    int    alphaint; 
00132 
00134    char* temp_dualfile;
00135 };
00136 
00137 //############################################################################
00138 
00147 class VOL_dvector {
00148 public:
00150    double* v;
00152    int sz;
00153 
00154 public:
00156    VOL_dvector(const int s) {
00157       VOL_TEST_SIZE(s);
00158       v = new double[sz = s];
00159    }
00161    VOL_dvector() : v(0), sz(0) {}
00163    VOL_dvector(const VOL_dvector& x) : v(0), sz(0) {
00164       sz = x.sz;
00165       if (sz > 0) {
00166         v = new double[sz];
00167         std::memcpy(v, x.v, sz * sizeof(double));
00168       }
00169    }
00171    ~VOL_dvector() { delete[] v; }
00172 
00174    inline int size() const {return sz;}
00175 
00177    inline double& operator[](const int i) {
00178       VOL_TEST_INDEX(i, sz);
00179       return v[i];
00180    }
00181 
00183    inline double operator[](const int i) const {
00184       VOL_TEST_INDEX(i, sz);
00185       return v[i];
00186    }
00187 
00190    inline void clear() {
00191       delete[] v;
00192       v = 0;
00193       sz = 0;
00194    }
00197    inline void cc(const double gamma, const VOL_dvector& w) {
00198       if (sz != w.sz) {
00199          printf("bad VOL_dvector sizes\n");
00200          abort();
00201       }
00202       double * p_v = v - 1;
00203       const double * p_w = w.v - 1;
00204       const double * const p_e = v + sz;
00205       const double one_gamma = 1.0 - gamma;
00206       while ( ++p_v != p_e ){
00207          *p_v = one_gamma * (*p_v) + gamma * (*++p_w);
00208       }
00209    }
00210 
00213    inline void allocate(const int s) {
00214       VOL_TEST_SIZE(s);
00215       delete[] v;                 
00216       v = new double[sz = s];
00217    }
00218 
00220    inline void swap(VOL_dvector& w) {
00221       std::swap(v, w.v);
00222       std::swap(sz, w.sz);
00223    }
00224 
00226    VOL_dvector& operator=(const VOL_dvector& w);
00228    VOL_dvector& operator=(const double w);
00229 };
00230 
00231 //-----------------------------------------------------------------------------
00241 class VOL_ivector {
00242 public:
00244    int* v;
00246    int sz;
00247 public:
00249    VOL_ivector(const int s) {
00250       VOL_TEST_SIZE(s);
00251       v = new int[sz = s];
00252    }
00254    VOL_ivector() : v(0), sz(0) {}
00256    VOL_ivector(const VOL_ivector& x) {
00257       sz = x.sz;
00258       if (sz > 0) {
00259          v = new int[sz];
00260          std::memcpy(v, x.v, sz * sizeof(int));
00261       }
00262    }
00264    ~VOL_ivector(){
00265       delete [] v;
00266    }
00267 
00269    inline int size() const { return sz; }
00271    inline int& operator[](const int i) {
00272       VOL_TEST_INDEX(i, sz);
00273       return v[i];
00274    }
00275 
00277    inline int operator[](const int i) const {
00278       VOL_TEST_INDEX(i, sz);
00279       return v[i];
00280    }
00281 
00284    inline void clear() {
00285       delete[] v;
00286       v = 0;
00287       sz = 0;
00288    }
00289 
00292    inline void allocate(const int s) {
00293       VOL_TEST_SIZE(s);
00294       delete[] v;
00295       v = new int[sz = s];
00296    }
00297 
00299    inline void swap(VOL_ivector& w) {
00300       std::swap(v, w.v);
00301       std::swap(sz, w.sz);
00302    }
00303 
00305    VOL_ivector& operator=(const VOL_ivector& v);      
00307    VOL_ivector& operator=(const int w);
00308 };
00309 
00310 //############################################################################
00311 // A class describing a primal solution. This class is used only internally 
00312 class VOL_primal {
00313 public: 
00314    // objective value of this primal solution 
00315    double value;
00316    // the largest of the v[i]'s
00317    double viol;  
00318    // primal solution  
00319    VOL_dvector x;
00320    // v=b-Ax, for the relaxed constraints
00321    VOL_dvector v; 
00322 
00323    VOL_primal(const int psize, const int dsize) : x(psize), v(dsize) {}
00324    VOL_primal(const VOL_primal& primal) :
00325       value(primal.value), viol(primal.viol), x(primal.x), v(primal.v) {}
00326    ~VOL_primal() {}
00327    inline VOL_primal& operator=(const VOL_primal& p) {
00328       if (this == &p) 
00329          return *this;
00330       value = p.value;
00331       viol = p.viol;
00332       x = p.x;
00333       v = p.v;
00334       return *this;
00335    }
00336 
00337    // convex combination. data members in this will be overwritten
00338    // convex combination between two primal solutions
00339    // x <-- alpha x + (1 - alpha) p.x
00340    // v <-- alpha v + (1 - alpha) p.v
00341    inline void cc(const double alpha, const VOL_primal& p) {
00342       value = alpha * p.value + (1.0 - alpha) * value;
00343       x.cc(alpha, p.x);
00344       v.cc(alpha, p.v);
00345    }
00346    // find maximum of v[i]
00347    void find_max_viol(const VOL_dvector& dual_lb, 
00348                       const VOL_dvector& dual_ub);
00349 };
00350 
00351 //-----------------------------------------------------------------------------
00352 // A class describing a dual solution. This class is used only internally 
00353 class VOL_dual {
00354 public:
00355    // lagrangian value
00356    double lcost; 
00357    // reduced costs * (pstar-primal)
00358    double xrc;
00359    // this information is only printed
00360    // dual vector
00361    VOL_dvector u; 
00362 
00363    VOL_dual(const int dsize) : u(dsize) { u = 0.0;}
00364    VOL_dual(const VOL_dual& dual) :
00365       lcost(dual.lcost), xrc(dual.xrc), u(dual.u) {}
00366    ~VOL_dual() {}
00367    inline VOL_dual& operator=(const VOL_dual& p) {
00368       if (this == &p) 
00369          return *this;
00370       lcost = p.lcost;
00371       xrc = p.xrc;
00372       u = p.u;
00373       return *this;
00374    }
00375    // dual step
00376    void   step(const double target, const double lambda,
00377                const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
00378                const VOL_dvector& v);
00379    double ascent(const VOL_dvector& v, const VOL_dvector& last_u) const;
00380    void   compute_xrc(const VOL_dvector& pstarx, const VOL_dvector& primalx,
00381                       const VOL_dvector& rc);
00382 
00383 };
00384 
00385 
00386 //############################################################################
00387 /* here we check whether an iteration is green, yellow or red. Also according
00388    to this information we decide whether lambda should be changed */
00389 class VOL_swing {
00390 private:
00391    VOL_swing(const VOL_swing&);
00392    VOL_swing& operator=(const VOL_swing&);
00393 public:
00394    enum condition {green, yellow, red} lastswing;
00395    int lastgreeniter, lastyellowiter, lastrediter;
00396    int ngs, nrs, nys;
00397    int rd;
00398    
00399    VOL_swing() {
00400       lastgreeniter = lastyellowiter = lastrediter = 0;
00401       ngs = nrs = nys = 0;
00402    }
00403    ~VOL_swing(){}
00404 
00405    inline void cond(const VOL_dual& dual, 
00406                     const double lcost, const double ascent, const int iter) {
00407       double eps = 1.e-3;
00408 
00409       if (ascent > 0.0  &&  lcost > dual.lcost + eps) {
00410          lastswing = green;
00411          lastgreeniter = iter;
00412          ++ngs;
00413          rd = 0;
00414       } else { 
00415          if (ascent <= 0  &&  lcost > dual.lcost) {
00416             lastswing = yellow;
00417             lastyellowiter = iter;
00418             ++nys;
00419             rd = 0;
00420          } else {
00421             lastswing = red;
00422             lastrediter = iter;
00423             ++nrs;
00424             rd = 1;
00425          }
00426       }
00427    }
00428 
00429    inline double
00430    lfactor(const VOL_parms& parm, const double lambda, const int iter) {
00431       double lambdafactor = 1.0;
00432       double eps = 5.e-4;
00433       int cons;
00434 
00435       switch (lastswing) {
00436        case green:
00437          cons = iter - VolMax(lastyellowiter, lastrediter);
00438          if (parm.printflag & 4)
00439             printf("      G: Consecutive Gs = %3d\n\n", cons);
00440          if (cons >= parm.greentestinvl && lambda < 2.0) {
00441             lastgreeniter = lastyellowiter = lastrediter = iter;
00442             lambdafactor = 2.0;
00443             if (parm.printflag & 2)
00444                printf("\n ---- increasing lamda to %g ----\n\n",
00445                       lambda * lambdafactor); 
00446          }
00447          break;
00448       
00449        case yellow:
00450          cons = iter - VolMax(lastgreeniter, lastrediter);
00451          if (parm.printflag & 4)
00452             printf("      Y: Consecutive Ys = %3d\n\n", cons);
00453          if (cons >= parm.yellowtestinvl) {
00454             lastgreeniter = lastyellowiter = lastrediter = iter;
00455             lambdafactor = 1.1;
00456             if (parm.printflag & 2)
00457                printf("\n **** increasing lamda to %g *****\n\n",
00458                       lambda * lambdafactor);
00459          }
00460          break;
00461       
00462        case red:
00463          cons = iter - VolMax(lastgreeniter, lastyellowiter);
00464          if (parm.printflag & 4)
00465             printf("      R: Consecutive Rs = %3d\n\n", cons);
00466          if (cons >= parm.redtestinvl && lambda > eps) {
00467             lastgreeniter = lastyellowiter = lastrediter = iter;
00468             lambdafactor = 0.67;
00469             if (parm.printflag & 2)
00470                printf("\n **** decreasing lamda to %g *****\n\n",
00471                       lambda * lambdafactor);
00472          } 
00473          break;
00474       }
00475       return lambdafactor;
00476    }
00477 
00478    inline void
00479    print() {
00480       printf("**** G= %i, Y= %i, R= %i ****\n", ngs, nys, nrs);
00481       ngs = nrs = nys = 0;  
00482    }
00483 };
00484 
00485 //############################################################################
00486 /* alpha should be decreased if after some number of iterations the objective
00487    has increased less that 1% */
00488 class VOL_alpha_factor {
00489 private:
00490    VOL_alpha_factor(const VOL_alpha_factor&);
00491    VOL_alpha_factor& operator=(const VOL_alpha_factor&);
00492 public:
00493    double lastvalue;
00494 
00495    VOL_alpha_factor() {lastvalue = -DBL_MAX;}
00496    ~VOL_alpha_factor() {}
00497 
00498    inline double factor(const VOL_parms& parm, const double lcost,
00499                         const double alpha) {
00500       if (alpha < parm.alphamin)
00501          return 1.0;
00502       register const double ll = VolAbs(lcost);
00503       const double x = ll > 10 ? (lcost-lastvalue)/ll : (lcost-lastvalue);
00504       lastvalue = lcost;
00505       return (x <= 0.01) ? parm.alphafactor : 1.0;
00506    }
00507 };
00508 
00509 //############################################################################
00510 /* here we compute the norm of the conjugate direction -hh-, the norm of the
00511    subgradient -norm-, the inner product between the subgradient and the 
00512    last conjugate direction -vh-, and the inner product between the new
00513    conjugate direction and the subgradient */
00514 class VOL_vh {
00515 private:
00516    VOL_vh(const VOL_vh&);
00517    VOL_vh& operator=(const VOL_vh&);
00518 public:
00519    double hh;
00520    double norm;
00521    double vh;
00522    double asc;
00523 
00524    VOL_vh(const double alpha,
00525           const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
00526           const VOL_dvector& v, const VOL_dvector& vstar,
00527           const VOL_dvector& u);
00528    ~VOL_vh(){}
00529 };
00530 
00531 //############################################################################
00532 /* here we compute different parameter to be printed. v2 is the square of 
00533    the norm of the subgradient. vu is the inner product between the dual
00534    variables and the subgradient. vabs is the maximum absolute value of
00535    the violations of pstar. asc is the inner product between the conjugate
00536    direction and the subgradient */
00537 class VOL_indc {
00538 private:
00539    VOL_indc(const VOL_indc&);
00540    VOL_indc& operator=(const VOL_indc&);
00541 public:
00542    double v2;
00543    double vu;
00544    double vabs;
00545    double asc;
00546 
00547 public:
00548    VOL_indc(const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
00549             const VOL_primal& primal, const VOL_primal& pstar,
00550             const VOL_dual& dual);
00551    ~VOL_indc() {}
00552 };
00553 
00554 //#############################################################################
00555 
00562 class VOL_user_hooks {
00563 public:
00564    virtual ~VOL_user_hooks() {}
00565 public:
00566    // for all hooks: return value of -1 means that volume should quit
00571    virtual int compute_rc(const VOL_dvector& u, VOL_dvector& rc) = 0;
00572 
00581    virtual int solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
00582                                 double& lcost, VOL_dvector& x, VOL_dvector& v,
00583                                 double& pcost) = 0;
00590    virtual int heuristics(const VOL_problem& p, 
00591                           const VOL_dvector& x, double& heur_val) = 0;
00592 };
00593 
00594 //#############################################################################
00595 
00604 class VOL_problem {
00605 private:
00606    VOL_problem(const VOL_problem&);
00607    VOL_problem& operator=(const VOL_problem&);
00608    void set_default_parm();
00609    // ############ INPUT fields ########################
00610 public: 
00614    VOL_problem();
00617    VOL_problem(const char *filename);
00619    ~VOL_problem();
00621 
00627    int solve(VOL_user_hooks& hooks, const bool use_preset_dual = false);
00629 
00630 private: 
00634    double alpha_; 
00636    double lambda_;
00637    // This union is here for padding (so that data members would be
00638    // double-aligned on x86 CPU
00639    union {
00641       int iter_;
00642       double __pad0;
00643    };
00645 
00646 public:
00647   
00651    double value;
00653    VOL_dvector dsol;
00655    VOL_dvector psol;
00657    VOL_dvector viol;
00659 
00663    VOL_parms parm;
00665    int psize;        
00667    int dsize;      
00670    VOL_dvector dual_lb;
00673    VOL_dvector dual_ub;
00675 
00676 public:
00680    int    iter() const { return iter_; }
00682    double alpha() const { return alpha_; }
00684    double lambda() const { return lambda_; }
00686 
00687 private:
00691    void read_params(const char* filename);
00692 
00694    int initialize(const bool use_preset_dual);
00695 
00697    void print_info(const int iter,
00698                    const VOL_primal& primal, const VOL_primal& pstar,
00699                    const VOL_dual& dual);
00700 
00703    double readjust_target(const double oldtarget, const double lcost) const;
00704 
00712    double power_heur(const VOL_primal& primal, const VOL_primal& pstar,
00713                      const VOL_dual& dual) const;
00715 };
00716 
00717 #endif

Generated on Sun Oct 16 03:19:09 2011 for Vol by  doxygen 1.4.7