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

- File List
- File Members

# /home/coin/SVN-release/Alps-1.2.1/Osi/test/OsiTestSolver.hpp

Go to the documentation of this file.

    00001 // Copyright (C) 2000, International Business Machines
    00002 // Corporation and others.  All Rights Reserved.
    00003
    00004 // this is a copy of VolVolume (stable/1.1 rev. 233)
    00005
    00006 #ifndef __OSITESTSOLVER_HPP__
    00007 #define __OSITESTSOLVER_HPP__
    00008
    00009 #include <cfloat>
    00010 #include <algorithm>
    00011 #include <cstdio>
    00012 #include <cmath>
    00013
    00014 #ifndef VOL_DEBUG
    00015 // When VOL_DEBUG is 1, we check vector indices
    00016 #define VOL_DEBUG 0
    00017 #endif
    00018
    00019 template <class T> static inline T
    00020 VolMax(register const T x, register const T y) {
    00021    return ((x) > (y)) ? (x) : (y);
    00022 }
    00023
    00024 template <class T> static inline T
    00025 VolAbs(register const T x) {
    00026    return ((x) > 0) ? (x) : -(x);
    00027 }
    00028
    00029 //############################################################################
    00030
    00031 #if defined(VOL_DEBUG) && (VOL_DEBUG != 0)
    00032 #define VOL_TEST_INDEX(i, size)                 \
    00033 {                                               \
    00034    if ((i) < 0 || (i) >= (size)) {              \
    00035       printf("bad VOL_?vector index\n");        \
    00036       abort();                                  \
    00037    }                                            \
    00038 }
    00039 #define VOL_TEST_SIZE(size)                     \
    00040 {                                               \
    00041    if (s <= 0) {                                \
    00042       printf("bad VOL_?vector size\n");         \
    00043       abort();                                  \
    00044    }                                            \
    00045 }
    00046 #else
    00047 #define VOL_TEST_INDEX(i, size)
    00048 #define VOL_TEST_SIZE(size)
    00049 #endif
    00050
    00051 //############################################################################
    00052
    00053 class VOL_dvector;
    00054 class VOL_ivector;
    00055 class VOL_primal;
    00056 class VOL_dual;
    00057 class VOL_swing;
    00058 class VOL_alpha_factor;
    00059 class VOL_vh;
    00060 class VOL_indc;
    00061 class VOL_user_hooks;
    00062 class VOL_problem;
    00063
    00064 //############################################################################
    00065
    00069 struct VOL_parms {
    00071    double lambdainit;
    00073    double alphainit;
    00075    double alphamin;
    00077    double alphafactor;
    00078
    00080    double ubinit;
    00081
    00083    double primal_abs_precision;
    00085    double gap_abs_precision;
    00087    double gap_rel_precision;
    00089    double granularity;
    00090
    00093    double minimum_rel_ascent;
    00095    int    ascent_first_check;
    00098    int    ascent_check_invl;
    00099
    00101    int    maxsgriters;
    00102
    00113    int    printflag;
    00115    int    printinvl;
    00117    int    heurinvl;
    00118
    00121    int greentestinvl;
    00124    int yellowtestinvl;
    00127    int redtestinvl;
    00128
    00130    int    alphaint;
    00131
    00133    char* temp_dualfile;
    00134 };
    00135
    00136 //############################################################################
    00137
    00146 class VOL_dvector {
    00147 public:
    00149    double* v;
    00151    int sz;
    00152
    00153 public:
    00155    VOL_dvector(const int s) {
    00156       VOL_TEST_SIZE(s);
    00157       v = new double[sz = s];
    00158    }
    00160    VOL_dvector() : v(0), sz(0) {}
    00162    VOL_dvector(const VOL_dvector& x) : v(0), sz(0) {
    00163       sz = x.sz;
    00164       if (sz > 0) {
    00165          v = new double[sz];
    00166          std::copy(x.v, x.v + sz, v);
    00167       }
    00168    }
    00170    ~VOL_dvector() { delete[] v; }
    00171
    00173    inline int size() const {return sz;}
    00174
    00176    inline double& operator[](const int i) {
    00177       VOL_TEST_INDEX(i, sz);
    00178       return v[i];
    00179    }
    00180
    00182    inline double operator[](const int i) const {
    00183       VOL_TEST_INDEX(i, sz);
    00184       return v[i];
    00185    }
    00186
    00189    inline void clear() {
    00190       delete[] v;
    00191       v = 0;
    00192       sz = 0;
    00193    }
    00196    inline void cc(const double gamma, const VOL_dvector& w) {
    00197       if (sz != w.sz) {
    00198          printf("bad VOL_dvector sizes\n");
    00199          abort();
    00200       }
    00201       double * p_v = v - 1;
    00202       const double * p_w = w.v - 1;
    00203       const double * const p_e = v + sz;
    00204       const double one_gamma = 1.0 - gamma;
    00205       while ( ++p_v != p_e ){
    00206          *p_v = one_gamma * (*p_v) + gamma * (*++p_w);
    00207       }
    00208    }
    00209
    00212    inline void allocate(const int s) {
    00213       VOL_TEST_SIZE(s);
    00214       delete[] v;
    00215       v = new double[sz = s];
    00216    }
    00217
    00219    inline void swap(VOL_dvector& w) {
    00220       std::swap(v, w.v);
    00221       std::swap(sz, w.sz);
    00222    }
    00223
    00225    VOL_dvector& operator=(const VOL_dvector& w);
    00227    VOL_dvector& operator=(const double w);
    00228 };
    00229
    00230 //-----------------------------------------------------------------------------
    00240 class VOL_ivector {
    00241 public:
    00243    int* v;
    00245    int sz;
    00246 public:
    00248    VOL_ivector(const int s) {
    00249       VOL_TEST_SIZE(s);
    00250       v = new int[sz = s];
    00251    }
    00253    VOL_ivector() : v(0), sz(0) {}
    00255    VOL_ivector(const VOL_ivector& x) {
    00256       sz = x.sz;
    00257       if (sz > 0) {
    00258          v = new int[sz];
    00259          std::copy(x.v, x.v + sz, v);
    00260       }
    00261    }
    00263    ~VOL_ivector(){
    00264       delete [] v;
    00265    }
    00266
    00268    inline int size() const { return sz; }
    00270    inline int& operator[](const int i) {
    00271       VOL_TEST_INDEX(i, sz);
    00272       return v[i];
    00273    }
    00274
    00276    inline int operator[](const int i) const {
    00277       VOL_TEST_INDEX(i, sz);
    00278       return v[i];
    00279    }
    00280
    00283    inline void clear() {
    00284       delete[] v;
    00285       v = 0;
    00286       sz = 0;
    00287    }
    00288
    00291    inline void allocate(const int s) {
    00292       VOL_TEST_SIZE(s);
    00293       delete[] v;
    00294       v = new int[sz = s];
    00295    }
    00296
    00298    inline void swap(VOL_ivector& w) {
    00299       std::swap(v, w.v);
    00300       std::swap(sz, w.sz);
    00301    }
    00302
    00304    VOL_ivector& operator=(const VOL_ivector& v);
    00306    VOL_ivector& operator=(const int w);
    00307 };
    00308
    00309 //############################################################################
    00310 // A class describing a primal solution. This class is used only internally
    00311 class VOL_primal {
    00312 public:
    00313    // objective value of this primal solution
    00314    double value;
    00315    // the largest of the v[i]'s
    00316    double viol;
    00317    // primal solution
    00318    VOL_dvector x;
    00319    // v=b-Ax, for the relaxed constraints
    00320    VOL_dvector v;
    00321
    00322    VOL_primal(const int psize, const int dsize) : x(psize), v(dsize) {}
    00323    VOL_primal(const VOL_primal& primal) :
    00324       value(primal.value), viol(primal.viol), x(primal.x), v(primal.v) {}
    00325    ~VOL_primal() {}
    00326    inline VOL_primal& operator=(const VOL_primal& p) {
    00327       if (this == &p)
    00328          return *this;
    00329       value = p.value;
    00330       viol = p.viol;
    00331       x = p.x;
    00332       v = p.v;
    00333       return *this;
    00334    }
    00335
    00336    // convex combination. data members in this will be overwritten
    00337    // convex combination between two primal solutions
    00338    // x <-- alpha x + (1 - alpha) p.x
    00339    // v <-- alpha v + (1 - alpha) p.v
    00340    inline void cc(const double alpha, const VOL_primal& p) {
    00341       value = alpha * p.value + (1.0 - alpha) * value;
    00342       x.cc(alpha, p.x);
    00343       v.cc(alpha, p.v);
    00344    }
    00345    // find maximum of v[i]
    00346    void find_max_viol(const VOL_dvector& dual_lb,
    00347                       const VOL_dvector& dual_ub);
    00348 };
    00349
    00350 //-----------------------------------------------------------------------------
    00351 // A class describing a dual solution. This class is used only internally
    00352 class VOL_dual {
    00353 public:
    00354    // lagrangian value
    00355    double lcost;
    00356    // reduced costs * (pstar-primal)
    00357    double xrc;
    00358    // this information is only printed
    00359    // dual vector
    00360    VOL_dvector u;
    00361
    00362    VOL_dual(const int dsize) : u(dsize) { u = 0.0;}
    00363    VOL_dual(const VOL_dual& dual) :
    00364       lcost(dual.lcost), xrc(dual.xrc), u(dual.u) {}
    00365    ~VOL_dual() {}
    00366    inline VOL_dual& operator=(const VOL_dual& p) {
    00367       if (this == &p)
    00368          return *this;
    00369       lcost = p.lcost;
    00370       xrc = p.xrc;
    00371       u = p.u;
    00372       return *this;
    00373    }
    00374    // dual step
    00375    void   step(const double target, const double lambda,
    00376                const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
    00377                const VOL_dvector& v);
    00378    double ascent(const VOL_dvector& v, const VOL_dvector& last_u) const;
    00379    void   compute_xrc(const VOL_dvector& pstarx, const VOL_dvector& primalx,
    00380                       const VOL_dvector& rc);
    00381
    00382 };
    00383
    00384
    00385 //############################################################################
    00386 /* here we check whether an iteration is green, yellow or red. Also according
    00387    to this information we decide whether lambda should be changed */
    00388 class VOL_swing {
    00389 private:
    00390    VOL_swing(const VOL_swing&);
    00391    VOL_swing& operator=(const VOL_swing&);
    00392 public:
    00393    enum condition {green, yellow, red} lastswing;
    00394    int lastgreeniter, lastyellowiter, lastrediter;
    00395    int ngs, nrs, nys;
    00396    int rd;
    00397
    00398    VOL_swing() {
    00399       lastgreeniter = lastyellowiter = lastrediter = 0;
    00400       ngs = nrs = nys = 0;
    00401    }
    00402    ~VOL_swing(){}
    00403
    00404    inline void cond(const VOL_dual& dual,
    00405                     const double lcost, const double ascent, const int iter) {
    00406       double eps = 1.e-3;
    00407
    00408       if (ascent > 0.0  &&  lcost > dual.lcost + eps) {
    00409          lastswing = green;
    00410          lastgreeniter = iter;
    00411          ++ngs;
    00412          rd = 0;
    00413       } else {
    00414          if (ascent <= 0  &&  lcost > dual.lcost) {
    00415             lastswing = yellow;
    00416             lastyellowiter = iter;
    00417             ++nys;
    00418             rd = 0;
    00419          } else {
    00420             lastswing = red;
    00421             lastrediter = iter;
    00422             ++nrs;
    00423             rd = 1;
    00424          }
    00425       }
    00426    }
    00427
    00428    inline double
    00429    lfactor(const VOL_parms& parm, const double lambda, const int iter) {
    00430       double lambdafactor = 1.0;
    00431       double eps = 5.e-4;
    00432       int cons;
    00433
    00434       switch (lastswing) {
    00435        case green:
    00436          cons = iter - VolMax(lastyellowiter, lastrediter);
    00437          if (parm.printflag & 4)
    00438             printf("      G: Consecutive Gs = %3d\n\n", cons);
    00439          if (cons >= parm.greentestinvl && lambda < 2.0) {
    00440             lastgreeniter = lastyellowiter = lastrediter = iter;
    00441             lambdafactor = 2.0;
    00442             if (parm.printflag & 2)
    00443                printf("\n ---- increasing lamda to %g ----\n\n",
    00444                       lambda * lambdafactor);
    00445          }
    00446          break;
    00447
    00448        case yellow:
    00449          cons = iter - VolMax(lastgreeniter, lastrediter);
    00450          if (parm.printflag & 4)
    00451             printf("      Y: Consecutive Ys = %3d\n\n", cons);
    00452          if (cons >= parm.yellowtestinvl) {
    00453             lastgreeniter = lastyellowiter = lastrediter = iter;
    00454             lambdafactor = 1.1;
    00455             if (parm.printflag & 2)
    00456                printf("\n **** increasing lamda to %g *****\n\n",
    00457                       lambda * lambdafactor);
    00458          }
    00459          break;
    00460
    00461        case red:
    00462          cons = iter - VolMax(lastgreeniter, lastyellowiter);
    00463          if (parm.printflag & 4)
    00464             printf("      R: Consecutive Rs = %3d\n\n", cons);
    00465          if (cons >= parm.redtestinvl && lambda > eps) {
    00466             lastgreeniter = lastyellowiter = lastrediter = iter;
    00467             lambdafactor = 0.67;
    00468             if (parm.printflag & 2)
    00469                printf("\n **** decreasing lamda to %g *****\n\n",
    00470                       lambda * lambdafactor);
    00471          }
    00472          break;
    00473       }
    00474       return lambdafactor;
    00475    }
    00476
    00477    inline void
    00478    print() {
    00479       printf("**** G= %i, Y= %i, R= %i ****\n", ngs, nys, nrs);
    00480       ngs = nrs = nys = 0;
    00481    }
    00482 };
    00483
    00484 //############################################################################
    00485 /* alpha should be decreased if after some number of iterations the objective
    00486    has increased less that 1% */
    00487 class VOL_alpha_factor {
    00488 private:
    00489    VOL_alpha_factor(const VOL_alpha_factor&);
    00490    VOL_alpha_factor& operator=(const VOL_alpha_factor&);
    00491 public:
    00492    double lastvalue;
    00493
    00494    VOL_alpha_factor() {lastvalue = -DBL_MAX;}
    00495    ~VOL_alpha_factor() {}
    00496
    00497    inline double factor(const VOL_parms& parm, const double lcost,
    00498                         const double alpha) {
    00499       if (alpha < parm.alphamin)
    00500          return 1.0;
    00501       register const double ll = VolAbs(lcost);
    00502       const double x = ll > 10 ? (lcost-lastvalue)/ll : (lcost-lastvalue);
    00503       lastvalue = lcost;
    00504       return (x <= 0.01) ? parm.alphafactor : 1.0;
    00505    }
    00506 };
    00507
    00508 //############################################################################
    00509 /* here we compute the norm of the conjugate direction -hh-, the norm of the
    00510    subgradient -norm-, the inner product between the subgradient and the
    00511    last conjugate direction -vh-, and the inner product between the new
    00512    conjugate direction and the subgradient */
    00513 class VOL_vh {
    00514 private:
    00515    VOL_vh(const VOL_vh&);
    00516    VOL_vh& operator=(const VOL_vh&);
    00517 public:
    00518    double hh;
    00519    double norm;
    00520    double vh;
    00521    double asc;
    00522
    00523    VOL_vh(const double alpha,
    00524           const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
    00525           const VOL_dvector& v, const VOL_dvector& vstar,
    00526           const VOL_dvector& u);
    00527    ~VOL_vh(){}
    00528 };
    00529
    00530 //############################################################################
    00531 /* here we compute different parameter to be printed. v2 is the square of
    00532    the norm of the subgradient. vu is the inner product between the dual
    00533    variables and the subgradient. vabs is the maximum absolute value of
    00534    the violations of pstar. asc is the inner product between the conjugate
    00535    direction and the subgradient */
    00536 class VOL_indc {
    00537 private:
    00538    VOL_indc(const VOL_indc&);
    00539    VOL_indc& operator=(const VOL_indc&);
    00540 public:
    00541    double v2;
    00542    double vu;
    00543    double vabs;
    00544    double asc;
    00545
    00546 public:
    00547    VOL_indc(const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
    00548             const VOL_primal& primal, const VOL_primal& pstar,
    00549             const VOL_dual& dual);
    00550    ~VOL_indc() {}
    00551 };
    00552
    00553 //#############################################################################
    00554
    00561 class VOL_user_hooks {
    00562 public:
    00563    virtual ~VOL_user_hooks() {}
    00564 public:
    00565    // for all hooks: return value of -1 means that volume should quit
    00570    virtual int compute_rc(const VOL_dvector& u, VOL_dvector& rc) = 0;
    00571
    00580    virtual int solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
    00581                                 double& lcost, VOL_dvector& x, VOL_dvector& v,
    00582                                 double& pcost) = 0;
    00589    virtual int heuristics(const VOL_problem& p,
    00590                           const VOL_dvector& x, double& heur_val) = 0;
    00591 };
    00592
    00593 //#############################################################################
    00594
    00603 class VOL_problem {
    00604 private:
    00605    VOL_problem(const VOL_problem&);
    00606    VOL_problem& operator=(const VOL_problem&);
    00607    void set_default_parm();
    00608    // ############ INPUT fields ########################
    00609 public:
    00613    VOL_problem();
    00616    VOL_problem(const char *filename);
    00618    ~VOL_problem();
    00620
    00626    int solve(VOL_user_hooks& hooks, const bool use_preset_dual = false);
    00628
    00629 private:
    00633    double alpha_;
    00635    double lambda_;
    00636    // This union is here for padding (so that data members would be
    00637    // double-aligned on x86 CPU
    00638    union {
    00640       int iter_;
    00641       double __pad0;
    00642    };
    00644
    00645 public:
    00646
    00650    double value;
    00652    VOL_dvector dsol;
    00654    VOL_dvector psol;
    00656    VOL_dvector viol;
    00658
    00662    VOL_parms parm;
    00664    int psize;
    00666    int dsize;
    00669    VOL_dvector dual_lb;
    00672    VOL_dvector dual_ub;
    00674
    00675 public:
    00679    int    iter() const { return iter_; }
    00681    double alpha() const { return alpha_; }
    00683    double lambda() const { return lambda_; }
    00685
    00686 private:
    00690    void read_params(const char* filename);
    00691
    00693    int initialize(const bool use_preset_dual);
    00694
    00696    void print_info(const int iter,
    00697                    const VOL_primal& primal, const VOL_primal& pstar,
    00698                    const VOL_dual& dual);
    00699
    00702    double readjust_target(const double oldtarget, const double lcost) const;
    00703
    00711    double power_heur(const VOL_primal& primal, const VOL_primal& pstar,
    00712                      const VOL_dual& dual) const;
    00714 };
    00715
    00716 #endif