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

- File List
- File Members

# /home/coin/svn-release/Alps-1.4.6/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 // This file is licensed under the terms of Eclipse Public License (EPL).
    00004
    00005 // this is a copy of VolVolume (stable/1.1 rev. 233)
    00006
    00007 #ifndef __OSITESTSOLVER_HPP__
    00008 #define __OSITESTSOLVER_HPP__
    00009
    00010 #include <cfloat>
    00011 #include <algorithm>
    00012 #include <cstdio>
    00013 #include <cmath>
    00014
    00015 #include "CoinFinite.hpp"
    00016
    00017 #ifndef VOL_DEBUG
    00018 // When VOL_DEBUG is 1, we check vector indices
    00019 #define VOL_DEBUG 0
    00020 #endif
    00021
    00022 template <class T> static inline T
    00023 VolMax(register const T x, register const T y) {
    00024    return ((x) > (y)) ? (x) : (y);
    00025 }
    00026
    00027 template <class T> static inline T
    00028 VolAbs(register const T x) {
    00029    return ((x) > 0) ? (x) : -(x);
    00030 }
    00031
    00032 //############################################################################
    00033
    00034 #if defined(VOL_DEBUG) && (VOL_DEBUG != 0)
    00035 #define VOL_TEST_INDEX(i, size)                 \
    00036 {                                               \
    00037    if ((i) < 0 || (i) >= (size)) {              \
    00038       printf("bad VOL_?vector index\n");        \
    00039       abort();                                  \
    00040    }                                            \
    00041 }
    00042 #define VOL_TEST_SIZE(size)                     \
    00043 {                                               \
    00044    if (s <= 0) {                                \
    00045       printf("bad VOL_?vector size\n");         \
    00046       abort();                                  \
    00047    }                                            \
    00048 }
    00049 #else
    00050 #define VOL_TEST_INDEX(i, size)
    00051 #define VOL_TEST_SIZE(size)
    00052 #endif
    00053
    00054 //############################################################################
    00055
    00056 class VOL_dvector;
    00057 class VOL_ivector;
    00058 class VOL_primal;
    00059 class VOL_dual;
    00060 class VOL_swing;
    00061 class VOL_alpha_factor;
    00062 class VOL_vh;
    00063 class VOL_indc;
    00064 class VOL_user_hooks;
    00065 class VOL_problem;
    00066
    00067 //############################################################################
    00068
    00072 struct VOL_parms {
    00074    double lambdainit;
    00076    double alphainit;
    00078    double alphamin;
    00080    double alphafactor;
    00081
    00083    double ubinit;
    00084
    00086    double primal_abs_precision;
    00088    double gap_abs_precision;
    00090    double gap_rel_precision;
    00092    double granularity;
    00093
    00096    double minimum_rel_ascent;
    00098    int    ascent_first_check;
    00101    int    ascent_check_invl;
    00102
    00104    int    maxsgriters;
    00105
    00116    int    printflag;
    00118    int    printinvl;
    00120    int    heurinvl;
    00121
    00124    int greentestinvl;
    00127    int yellowtestinvl;
    00130    int redtestinvl;
    00131
    00133    int    alphaint;
    00134
    00136    char* temp_dualfile;
    00137 };
    00138
    00139 //############################################################################
    00140
    00149 class VOL_dvector {
    00150 public:
    00152    double* v;
    00154    int sz;
    00155
    00156 public:
    00158    VOL_dvector(const int s) {
    00159       VOL_TEST_SIZE(s);
    00160       v = new double[sz = s];
    00161    }
    00163    VOL_dvector() : v(0), sz(0) {}
    00165    VOL_dvector(const VOL_dvector& x) : v(0), sz(0) {
    00166       sz = x.sz;
    00167       if (sz > 0) {
    00168          v = new double[sz];
    00169          std::copy(x.v, x.v + sz, v);
    00170       }
    00171    }
    00173    ~VOL_dvector() { delete[] v; }
    00174
    00176    inline int size() const {return sz;}
    00177
    00179    inline double& operator[](const int i) {
    00180       VOL_TEST_INDEX(i, sz);
    00181       return v[i];
    00182    }
    00183
    00185    inline double operator[](const int i) const {
    00186       VOL_TEST_INDEX(i, sz);
    00187       return v[i];
    00188    }
    00189
    00192    inline void clear() {
    00193       delete[] v;
    00194       v = 0;
    00195       sz = 0;
    00196    }
    00199    inline void cc(const double gamma, const VOL_dvector& w) {
    00200       if (sz != w.sz) {
    00201          printf("bad VOL_dvector sizes\n");
    00202          abort();
    00203       }
    00204       double * p_v = v - 1;
    00205       const double * p_w = w.v - 1;
    00206       const double * const p_e = v + sz;
    00207       const double one_gamma = 1.0 - gamma;
    00208       while ( ++p_v != p_e ){
    00209          *p_v = one_gamma * (*p_v) + gamma * (*++p_w);
    00210       }
    00211    }
    00212
    00215    inline void allocate(const int s) {
    00216       VOL_TEST_SIZE(s);
    00217       delete[] v;
    00218       v = new double[sz = s];
    00219    }
    00220
    00222    inline void swap(VOL_dvector& w) {
    00223       std::swap(v, w.v);
    00224       std::swap(sz, w.sz);
    00225    }
    00226
    00228    VOL_dvector& operator=(const VOL_dvector& w);
    00230    VOL_dvector& operator=(const double w);
    00231 };
    00232
    00233 //-----------------------------------------------------------------------------
    00243 class VOL_ivector {
    00244 public:
    00246    int* v;
    00248    int sz;
    00249 public:
    00251    VOL_ivector(const int s) {
    00252       VOL_TEST_SIZE(s);
    00253       v = new int[sz = s];
    00254    }
    00256    VOL_ivector() : v(0), sz(0) {}
    00258    VOL_ivector(const VOL_ivector& x) {
    00259       sz = x.sz;
    00260       if (sz > 0) {
    00261          v = new int[sz];
    00262          std::copy(x.v, x.v + sz, v);
    00263       }
    00264    }
    00266    ~VOL_ivector(){
    00267       delete [] v;
    00268    }
    00269
    00271    inline int size() const { return sz; }
    00273    inline int& operator[](const int i) {
    00274       VOL_TEST_INDEX(i, sz);
    00275       return v[i];
    00276    }
    00277
    00279    inline int operator[](const int i) const {
    00280       VOL_TEST_INDEX(i, sz);
    00281       return v[i];
    00282    }
    00283
    00286    inline void clear() {
    00287       delete[] v;
    00288       v = 0;
    00289       sz = 0;
    00290    }
    00291
    00294    inline void allocate(const int s) {
    00295       VOL_TEST_SIZE(s);
    00296       delete[] v;
    00297       v = new int[sz = s];
    00298    }
    00299
    00301    inline void swap(VOL_ivector& w) {
    00302       std::swap(v, w.v);
    00303       std::swap(sz, w.sz);
    00304    }
    00305
    00307    VOL_ivector& operator=(const VOL_ivector& v);
    00309    VOL_ivector& operator=(const int w);
    00310 };
    00311
    00312 //############################################################################
    00313 // A class describing a primal solution. This class is used only internally
    00314 class VOL_primal {
    00315 public:
    00316    // objective value of this primal solution
    00317    double value;
    00318    // the largest of the v[i]'s
    00319    double viol;
    00320    // primal solution
    00321    VOL_dvector x;
    00322    // v=b-Ax, for the relaxed constraints
    00323    VOL_dvector v;
    00324
    00325    VOL_primal(const int psize, const int dsize) : x(psize), v(dsize) {}
    00326    VOL_primal(const VOL_primal& primal) :
    00327       value(primal.value), viol(primal.viol), x(primal.x), v(primal.v) {}
    00328    ~VOL_primal() {}
    00329    inline VOL_primal& operator=(const VOL_primal& p) {
    00330       if (this == &p)
    00331          return *this;
    00332       value = p.value;
    00333       viol = p.viol;
    00334       x = p.x;
    00335       v = p.v;
    00336       return *this;
    00337    }
    00338
    00339    // convex combination. data members in this will be overwritten
    00340    // convex combination between two primal solutions
    00341    // x <-- alpha x + (1 - alpha) p.x
    00342    // v <-- alpha v + (1 - alpha) p.v
    00343    inline void cc(const double alpha, const VOL_primal& p) {
    00344       value = alpha * p.value + (1.0 - alpha) * value;
    00345       x.cc(alpha, p.x);
    00346       v.cc(alpha, p.v);
    00347    }
    00348    // find maximum of v[i]
    00349    void find_max_viol(const VOL_dvector& dual_lb,
    00350                       const VOL_dvector& dual_ub);
    00351 };
    00352
    00353 //-----------------------------------------------------------------------------
    00354 // A class describing a dual solution. This class is used only internally
    00355 class VOL_dual {
    00356 public:
    00357    // lagrangian value
    00358    double lcost;
    00359    // reduced costs * (pstar-primal)
    00360    double xrc;
    00361    // this information is only printed
    00362    // dual vector
    00363    VOL_dvector u;
    00364
    00365    VOL_dual(const int dsize) : u(dsize) { u = 0.0;}
    00366    VOL_dual(const VOL_dual& dual) :
    00367       lcost(dual.lcost), xrc(dual.xrc), u(dual.u) {}
    00368    ~VOL_dual() {}
    00369    inline VOL_dual& operator=(const VOL_dual& p) {
    00370       if (this == &p)
    00371          return *this;
    00372       lcost = p.lcost;
    00373       xrc = p.xrc;
    00374       u = p.u;
    00375       return *this;
    00376    }
    00377    // dual step
    00378    void   step(const double target, const double lambda,
    00379                const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
    00380                const VOL_dvector& v);
    00381    double ascent(const VOL_dvector& v, const VOL_dvector& last_u) const;
    00382    void   compute_xrc(const VOL_dvector& pstarx, const VOL_dvector& primalx,
    00383                       const VOL_dvector& rc);
    00384
    00385 };
    00386
    00387
    00388 //############################################################################
    00389 /* here we check whether an iteration is green, yellow or red. Also according
    00390    to this information we decide whether lambda should be changed */
    00391 class VOL_swing {
    00392 private:
    00393    VOL_swing(const VOL_swing&);
    00394    VOL_swing& operator=(const VOL_swing&);
    00395 public:
    00396    enum condition {green, yellow, red} lastswing;
    00397    int lastgreeniter, lastyellowiter, lastrediter;
    00398    int ngs, nrs, nys;
    00399    int rd;
    00400
    00401    VOL_swing() {
    00402       lastgreeniter = lastyellowiter = lastrediter = 0;
    00403       ngs = nrs = nys = 0;
    00404    }
    00405    ~VOL_swing(){}
    00406
    00407    inline void cond(const VOL_dual& dual,
    00408                     const double lcost, const double ascent, const int iter) {
    00409       double eps = 1.e-3;
    00410
    00411       if (ascent > 0.0  &&  lcost > dual.lcost + eps) {
    00412          lastswing = green;
    00413          lastgreeniter = iter;
    00414          ++ngs;
    00415          rd = 0;
    00416       } else {
    00417          if (ascent <= 0  &&  lcost > dual.lcost) {
    00418             lastswing = yellow;
    00419             lastyellowiter = iter;
    00420             ++nys;
    00421             rd = 0;
    00422          } else {
    00423             lastswing = red;
    00424             lastrediter = iter;
    00425             ++nrs;
    00426             rd = 1;
    00427          }
    00428       }
    00429    }
    00430
    00431    inline double
    00432    lfactor(const VOL_parms& parm, const double lambda, const int iter) {
    00433       double lambdafactor = 1.0;
    00434       double eps = 5.e-4;
    00435       int cons;
    00436
    00437       switch (lastswing) {
    00438        case green:
    00439          cons = iter - VolMax(lastyellowiter, lastrediter);
    00440          if (parm.printflag & 4)
    00441             printf("      G: Consecutive Gs = %3d\n\n", cons);
    00442          if (cons >= parm.greentestinvl && lambda < 2.0) {
    00443             lastgreeniter = lastyellowiter = lastrediter = iter;
    00444             lambdafactor = 2.0;
    00445             if (parm.printflag & 2)
    00446                printf("\n ---- increasing lamda to %g ----\n\n",
    00447                       lambda * lambdafactor);
    00448          }
    00449          break;
    00450
    00451        case yellow:
    00452          cons = iter - VolMax(lastgreeniter, lastrediter);
    00453          if (parm.printflag & 4)
    00454             printf("      Y: Consecutive Ys = %3d\n\n", cons);
    00455          if (cons >= parm.yellowtestinvl) {
    00456             lastgreeniter = lastyellowiter = lastrediter = iter;
    00457             lambdafactor = 1.1;
    00458             if (parm.printflag & 2)
    00459                printf("\n **** increasing lamda to %g *****\n\n",
    00460                       lambda * lambdafactor);
    00461          }
    00462          break;
    00463
    00464        case red:
    00465          cons = iter - VolMax(lastgreeniter, lastyellowiter);
    00466          if (parm.printflag & 4)
    00467             printf("      R: Consecutive Rs = %3d\n\n", cons);
    00468          if (cons >= parm.redtestinvl && lambda > eps) {
    00469             lastgreeniter = lastyellowiter = lastrediter = iter;
    00470             lambdafactor = 0.67;
    00471             if (parm.printflag & 2)
    00472                printf("\n **** decreasing lamda to %g *****\n\n",
    00473                       lambda * lambdafactor);
    00474          }
    00475          break;
    00476       }
    00477       return lambdafactor;
    00478    }
    00479
    00480    inline void
    00481    print() {
    00482       printf("**** G= %i, Y= %i, R= %i ****\n", ngs, nys, nrs);
    00483       ngs = nrs = nys = 0;
    00484    }
    00485 };
    00486
    00487 //############################################################################
    00488 /* alpha should be decreased if after some number of iterations the objective
    00489    has increased less that 1% */
    00490 class VOL_alpha_factor {
    00491 private:
    00492    VOL_alpha_factor(const VOL_alpha_factor&);
    00493    VOL_alpha_factor& operator=(const VOL_alpha_factor&);
    00494 public:
    00495    double lastvalue;
    00496
    00497    VOL_alpha_factor() {lastvalue = -COIN_DBL_MAX;}
    00498    ~VOL_alpha_factor() {}
    00499
    00500    inline double factor(const VOL_parms& parm, const double lcost,
    00501                         const double alpha) {
    00502       if (alpha < parm.alphamin)
    00503          return 1.0;
    00504       register const double ll = VolAbs(lcost);
    00505       const double x = ll > 10 ? (lcost-lastvalue)/ll : (lcost-lastvalue);
    00506       lastvalue = lcost;
    00507       return (x <= 0.01) ? parm.alphafactor : 1.0;
    00508    }
    00509 };
    00510
    00511 //############################################################################
    00512 /* here we compute the norm of the conjugate direction -hh-, the norm of the
    00513    subgradient -norm-, the inner product between the subgradient and the
    00514    last conjugate direction -vh-, and the inner product between the new
    00515    conjugate direction and the subgradient */
    00516 class VOL_vh {
    00517 private:
    00518    VOL_vh(const VOL_vh&);
    00519    VOL_vh& operator=(const VOL_vh&);
    00520 public:
    00521    double hh;
    00522    double norm;
    00523    double vh;
    00524    double asc;
    00525
    00526    VOL_vh(const double alpha,
    00527           const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
    00528           const VOL_dvector& v, const VOL_dvector& vstar,
    00529           const VOL_dvector& u);
    00530    ~VOL_vh(){}
    00531 };
    00532
    00533 //############################################################################
    00534 /* here we compute different parameter to be printed. v2 is the square of
    00535    the norm of the subgradient. vu is the inner product between the dual
    00536    variables and the subgradient. vabs is the maximum absolute value of
    00537    the violations of pstar. asc is the inner product between the conjugate
    00538    direction and the subgradient */
    00539 class VOL_indc {
    00540 private:
    00541    VOL_indc(const VOL_indc&);
    00542    VOL_indc& operator=(const VOL_indc&);
    00543 public:
    00544    double v2;
    00545    double vu;
    00546    double vabs;
    00547    double asc;
    00548
    00549 public:
    00550    VOL_indc(const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
    00551             const VOL_primal& primal, const VOL_primal& pstar,
    00552             const VOL_dual& dual);
    00553    ~VOL_indc() {}
    00554 };
    00555
    00556 //#############################################################################
    00557
    00564 class VOL_user_hooks {
    00565 public:
    00566    virtual ~VOL_user_hooks() {}
    00567 public:
    00568    // for all hooks: return value of -1 means that volume should quit
    00573    virtual int compute_rc(const VOL_dvector& u, VOL_dvector& rc) = 0;
    00574
    00583    virtual int solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
    00584                                 double& lcost, VOL_dvector& x, VOL_dvector& v,
    00585                                 double& pcost) = 0;
    00592    virtual int heuristics(const VOL_problem& p,
    00593                           const VOL_dvector& x, double& heur_val) = 0;
    00594 };
    00595
    00596 //#############################################################################
    00597
    00606 class VOL_problem {
    00607 private:
    00608    VOL_problem(const VOL_problem&);
    00609    VOL_problem& operator=(const VOL_problem&);
    00610    void set_default_parm();
    00611    // ############ INPUT fields ########################
    00612 public:
    00616    VOL_problem();
    00619    VOL_problem(const char *filename);
    00621    ~VOL_problem();
    00623
    00629    int solve(VOL_user_hooks& hooks, const bool use_preset_dual = false);
    00631
    00632 private:
    00636    double alpha_;
    00638    double lambda_;
    00639    // This union is here for padding (so that data members would be
    00640    // double-aligned on x86 CPU
    00641    union {
    00643       int iter_;
    00644       double __pad0;
    00645    };
    00647
    00648 public:
    00649
    00653    double value;
    00655    VOL_dvector dsol;
    00657    VOL_dvector psol;
    00659    VOL_dvector viol;
    00661
    00665    VOL_parms parm;
    00667    int psize;
    00669    int dsize;
    00672    VOL_dvector dual_lb;
    00675    VOL_dvector dual_ub;
    00677
    00678 public:
    00682    int    iter() const { return iter_; }
    00684    double alpha() const { return alpha_; }
    00686    double lambda() const { return lambda_; }
    00688
    00689 private:
    00693    void read_params(const char* filename);
    00694
    00696    int initialize(const bool use_preset_dual);
    00697
    00699    void print_info(const int iter,
    00700                    const VOL_primal& primal, const VOL_primal& pstar,
    00701                    const VOL_dual& dual);
    00702
    00705    double readjust_target(const double oldtarget, const double lcost) const;
    00706
    00714    double power_heur(const VOL_primal& primal, const VOL_primal& pstar,
    00715                      const VOL_dual& dual) const;
    00717 };
    00718
    00719 #endif