- Main Page
- Modules
- Namespaces
- Classes
- Files
- Related Pages

- File List
- File Members

# /home/coin/SVN-release/Cgl-0.53.1/Vol/src/VolVolume.hpp

Go to the documentation of this file.

    00001 // Copyright (C) 2000, International Business Machines
    00002 // Corporation and others.  All Rights Reserved.
    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 // When VOL_DEBUG is 1, we check vector indices
    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 // A class describing a primal solution. This class is used only internally
    00308 class VOL_primal {
    00309 public:
    00310    // objective value of this primal solution
    00311    double value;
    00312    // the largest of the v[i]'s
    00313    double viol;
    00314    // primal solution
    00315    VOL_dvector x;
    00316    // v=b-Ax, for the relaxed constraints
    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    // convex combination. data members in this will be overwritten
    00334    // convex combination between two primal solutions
    00335    // x <-- alpha x + (1 - alpha) p.x
    00336    // v <-- alpha v + (1 - alpha) p.v
    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    // find maximum of v[i]
    00343    void find_max_viol(const VOL_dvector& dual_lb,
    00344                       const VOL_dvector& dual_ub);
    00345 };
    00346
    00347 //-----------------------------------------------------------------------------
    00348 // A class describing a dual solution. This class is used only internally
    00349 class VOL_dual {
    00350 public:
    00351    // lagrangian value
    00352    double lcost;
    00353    // reduced costs * (pstar-primal)
    00354    double xrc;
    00355    // this information is only printed
    00356    // dual vector
    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    // dual step
    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 /* here we check whether an iteration is green, yellow or red. Also according
    00384    to this information we decide whether lambda should be changed */
    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 /* alpha should be decreased if after some number of iterations the objective
    00483    has increased less that 1% */
    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 /* here we compute the norm of the conjugate direction -hh-, the norm of the
    00507    subgradient -norm-, the inner product between the subgradient and the
    00508    last conjugate direction -vh-, and the inner product between the new
    00509    conjugate direction and the subgradient */
    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 /* here we compute different parameter to be printed. v2 is the square of
    00529    the norm of the subgradient. vu is the inner product between the dual
    00530    variables and the subgradient. vabs is the maximum absolute value of
    00531    the violations of pstar. asc is the inner product between the conjugate
    00532    direction and the subgradient */
    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    // for all hooks: return value of -1 means that volume should quit
    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    // ############ INPUT fields ########################
    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    // This union is here for padding (so that data members would be
    00634    // double-aligned on x86 CPU
    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