# 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 Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines