# 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 #include <cstdlib>
    00015
    00016 #include "CoinFinite.hpp"
    00017
    00018 #ifndef VOL_DEBUG
    00019 // When VOL_DEBUG is 1, we check vector indices
    00020 #define VOL_DEBUG 0
    00021 #endif
    00022
    00023 template <class T> static inline T
    00024 VolMax(register const T x, register const T y) {
    00025    return ((x) > (y)) ? (x) : (y);
    00026 }
    00027
    00028 template <class T> static inline T
    00029 VolAbs(register const T x) {
    00030    return ((x) > 0) ? (x) : -(x);
    00031 }
    00032
    00033 //############################################################################
    00034
    00035 #if defined(VOL_DEBUG) && (VOL_DEBUG != 0)
    00036 #define VOL_TEST_INDEX(i, size)                 \
    00037 {                                               \
    00038    if ((i) < 0 || (i) >= (size)) {              \
    00039       printf("bad VOL_?vector index\n");        \
    00040       abort();                                  \
    00041    }                                            \
    00042 }
    00043 #define VOL_TEST_SIZE(size)                     \
    00044 {                                               \
    00045    if (s <= 0) {                                \
    00046       printf("bad VOL_?vector size\n");         \
    00047       abort();                                  \
    00048    }                                            \
    00049 }
    00050 #else
    00051 #define VOL_TEST_INDEX(i, size)
    00052 #define VOL_TEST_SIZE(size)
    00053 #endif
    00054
    00055 //############################################################################
    00056
    00057 class VOL_dvector;
    00058 class VOL_ivector;
    00059 class VOL_primal;
    00060 class VOL_dual;
    00061 class VOL_swing;
    00062 class VOL_alpha_factor;
    00063 class VOL_vh;
    00064 class VOL_indc;
    00065 class VOL_user_hooks;
    00066 class VOL_problem;
    00067
    00068 //############################################################################
    00069
    00073 struct VOL_parms {
    00075    double lambdainit;
    00077    double alphainit;
    00079    double alphamin;
    00081    double alphafactor;
    00082
    00084    double ubinit;
    00085
    00087    double primal_abs_precision;
    00089    double gap_abs_precision;
    00091    double gap_rel_precision;
    00093    double granularity;
    00094
    00097    double minimum_rel_ascent;
    00099    int    ascent_first_check;
    00102    int    ascent_check_invl;
    00103
    00105    int    maxsgriters;
    00106
    00117    int    printflag;
    00119    int    printinvl;
    00121    int    heurinvl;
    00122
    00125    int greentestinvl;
    00128    int yellowtestinvl;
    00131    int redtestinvl;
    00132
    00134    int    alphaint;
    00135
    00137    char* temp_dualfile;
    00138 };
    00139
    00140 //############################################################################
    00141
    00150 class VOL_dvector {
    00151 public:
    00153    double* v;
    00155    int sz;
    00156
    00157 public:
    00159    VOL_dvector(const int s) {
    00160       VOL_TEST_SIZE(s);
    00161       v = new double[sz = s];
    00162    }
    00164    VOL_dvector() : v(0), sz(0) {}
    00166    VOL_dvector(const VOL_dvector& x) : v(0), sz(0) {
    00167       sz = x.sz;
    00168       if (sz > 0) {
    00169          v = new double[sz];
    00170          std::copy(x.v, x.v + sz, v);
    00171       }
    00172    }
    00174    ~VOL_dvector() { delete[] v; }
    00175
    00177    inline int size() const {return sz;}
    00178
    00180    inline double& operator[](const int i) {
    00181       VOL_TEST_INDEX(i, sz);
    00182       return v[i];
    00183    }
    00184
    00186    inline double operator[](const int i) const {
    00187       VOL_TEST_INDEX(i, sz);
    00188       return v[i];
    00189    }
    00190
    00193    inline void clear() {
    00194       delete[] v;
    00195       v = 0;
    00196       sz = 0;
    00197    }
    00200    inline void cc(const double gamma, const VOL_dvector& w) {
    00201       if (sz != w.sz) {
    00202          printf("bad VOL_dvector sizes\n");
    00203          abort();
    00204       }
    00205       double * p_v = v - 1;
    00206       const double * p_w = w.v - 1;
    00207       const double * const p_e = v + sz;
    00208       const double one_gamma = 1.0 - gamma;
    00209       while ( ++p_v != p_e ){
    00210          *p_v = one_gamma * (*p_v) + gamma * (*++p_w);
    00211       }
    00212    }
    00213
    00216    inline void allocate(const int s) {
    00217       VOL_TEST_SIZE(s);
    00218       delete[] v;
    00219       v = new double[sz = s];
    00220    }
    00221
    00223    inline void swap(VOL_dvector& w) {
    00224       std::swap(v, w.v);
    00225       std::swap(sz, w.sz);
    00226    }
    00227
    00229    VOL_dvector& operator=(const VOL_dvector& w);
    00231    VOL_dvector& operator=(const double w);
    00232 };
    00233
    00234 //-----------------------------------------------------------------------------
    00244 class VOL_ivector {
    00245 public:
    00247    int* v;
    00249    int sz;
    00250 public:
    00252    VOL_ivector(const int s) {
    00253       VOL_TEST_SIZE(s);
    00254       v = new int[sz = s];
    00255    }
    00257    VOL_ivector() : v(0), sz(0) {}
    00259    VOL_ivector(const VOL_ivector& x) {
    00260       sz = x.sz;
    00261       if (sz > 0) {
    00262          v = new int[sz];
    00263          std::copy(x.v, x.v + sz, v);
    00264       }
    00265    }
    00267    ~VOL_ivector(){
    00268       delete [] v;
    00269    }
    00270
    00272    inline int size() const { return sz; }
    00274    inline int& operator[](const int i) {
    00275       VOL_TEST_INDEX(i, sz);
    00276       return v[i];
    00277    }
    00278
    00280    inline int operator[](const int i) const {
    00281       VOL_TEST_INDEX(i, sz);
    00282       return v[i];
    00283    }
    00284
    00287    inline void clear() {
    00288       delete[] v;
    00289       v = 0;
    00290       sz = 0;
    00291    }
    00292
    00295    inline void allocate(const int s) {
    00296       VOL_TEST_SIZE(s);
    00297       delete[] v;
    00298       v = new int[sz = s];
    00299    }
    00300
    00302    inline void swap(VOL_ivector& w) {
    00303       std::swap(v, w.v);
    00304       std::swap(sz, w.sz);
    00305    }
    00306
    00308    VOL_ivector& operator=(const VOL_ivector& v);
    00310    VOL_ivector& operator=(const int w);
    00311 };
    00312
    00313 //############################################################################
    00314 // A class describing a primal solution. This class is used only internally
    00315 class VOL_primal {
    00316 public:
    00317    // objective value of this primal solution
    00318    double value;
    00319    // the largest of the v[i]'s
    00320    double viol;
    00321    // primal solution
    00322    VOL_dvector x;
    00323    // v=b-Ax, for the relaxed constraints
    00324    VOL_dvector v;
    00325
    00326    VOL_primal(const int psize, const int dsize) : x(psize), v(dsize) {}
    00327    VOL_primal(const VOL_primal& primal) :
    00328       value(primal.value), viol(primal.viol), x(primal.x), v(primal.v) {}
    00329    ~VOL_primal() {}
    00330    inline VOL_primal& operator=(const VOL_primal& p) {
    00331       if (this == &p)
    00332          return *this;
    00333       value = p.value;
    00334       viol = p.viol;
    00335       x = p.x;
    00336       v = p.v;
    00337       return *this;
    00338    }
    00339
    00340    // convex combination. data members in this will be overwritten
    00341    // convex combination between two primal solutions
    00342    // x <-- alpha x + (1 - alpha) p.x
    00343    // v <-- alpha v + (1 - alpha) p.v
    00344    inline void cc(const double alpha, const VOL_primal& p) {
    00345       value = alpha * p.value + (1.0 - alpha) * value;
    00346       x.cc(alpha, p.x);
    00347       v.cc(alpha, p.v);
    00348    }
    00349    // find maximum of v[i]
    00350    void find_max_viol(const VOL_dvector& dual_lb,
    00351                       const VOL_dvector& dual_ub);
    00352 };
    00353
    00354 //-----------------------------------------------------------------------------
    00355 // A class describing a dual solution. This class is used only internally
    00356 class VOL_dual {
    00357 public:
    00358    // lagrangian value
    00359    double lcost;
    00360    // reduced costs * (pstar-primal)
    00361    double xrc;
    00362    // this information is only printed
    00363    // dual vector
    00364    VOL_dvector u;
    00365
    00366    VOL_dual(const int dsize) : u(dsize) { u = 0.0;}
    00367    VOL_dual(const VOL_dual& dual) :
    00368       lcost(dual.lcost), xrc(dual.xrc), u(dual.u) {}
    00369    ~VOL_dual() {}
    00370    inline VOL_dual& operator=(const VOL_dual& p) {
    00371       if (this == &p)
    00372          return *this;
    00373       lcost = p.lcost;
    00374       xrc = p.xrc;
    00375       u = p.u;
    00376       return *this;
    00377    }
    00378    // dual step
    00379    void   step(const double target, const double lambda,
    00380                const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
    00381                const VOL_dvector& v);
    00382    double ascent(const VOL_dvector& v, const VOL_dvector& last_u) const;
    00383    void   compute_xrc(const VOL_dvector& pstarx, const VOL_dvector& primalx,
    00384                       const VOL_dvector& rc);
    00385
    00386 };
    00387
    00388
    00389 //############################################################################
    00390 /* here we check whether an iteration is green, yellow or red. Also according
    00391    to this information we decide whether lambda should be changed */
    00392 class VOL_swing {
    00393 private:
    00394    VOL_swing(const VOL_swing&);
    00395    VOL_swing& operator=(const VOL_swing&);
    00396 public:
    00397    enum condition {green, yellow, red} lastswing;
    00398    int lastgreeniter, lastyellowiter, lastrediter;
    00399    int ngs, nrs, nys;
    00400    int rd;
    00401
    00402    VOL_swing() {
    00403       lastgreeniter = lastyellowiter = lastrediter = 0;
    00404       ngs = nrs = nys = 0;
    00405    }
    00406    ~VOL_swing(){}
    00407
    00408    inline void cond(const VOL_dual& dual,
    00409                     const double lcost, const double ascent, const int iter) {
    00410       double eps = 1.e-3;
    00411
    00412       if (ascent > 0.0  &&  lcost > dual.lcost + eps) {
    00413          lastswing = green;
    00414          lastgreeniter = iter;
    00415          ++ngs;
    00416          rd = 0;
    00417       } else {
    00418          if (ascent <= 0  &&  lcost > dual.lcost) {
    00419             lastswing = yellow;
    00420             lastyellowiter = iter;
    00421             ++nys;
    00422             rd = 0;
    00423          } else {
    00424             lastswing = red;
    00425             lastrediter = iter;
    00426             ++nrs;
    00427             rd = 1;
    00428          }
    00429       }
    00430    }
    00431
    00432    inline double
    00433    lfactor(const VOL_parms& parm, const double lambda, const int iter) {
    00434       double lambdafactor = 1.0;
    00435       double eps = 5.e-4;
    00436       int cons;
    00437
    00438       switch (lastswing) {
    00439        case green:
    00440          cons = iter - VolMax(lastyellowiter, lastrediter);
    00441          if (parm.printflag & 4)
    00442             printf("      G: Consecutive Gs = %3d\n\n", cons);
    00443          if (cons >= parm.greentestinvl && lambda < 2.0) {
    00444             lastgreeniter = lastyellowiter = lastrediter = iter;
    00445             lambdafactor = 2.0;
    00446             if (parm.printflag & 2)
    00447                printf("\n ---- increasing lamda to %g ----\n\n",
    00448                       lambda * lambdafactor);
    00449          }
    00450          break;
    00451
    00452        case yellow:
    00453          cons = iter - VolMax(lastgreeniter, lastrediter);
    00454          if (parm.printflag & 4)
    00455             printf("      Y: Consecutive Ys = %3d\n\n", cons);
    00456          if (cons >= parm.yellowtestinvl) {
    00457             lastgreeniter = lastyellowiter = lastrediter = iter;
    00458             lambdafactor = 1.1;
    00459             if (parm.printflag & 2)
    00460                printf("\n **** increasing lamda to %g *****\n\n",
    00461                       lambda * lambdafactor);
    00462          }
    00463          break;
    00464
    00465        case red:
    00466          cons = iter - VolMax(lastgreeniter, lastyellowiter);
    00467          if (parm.printflag & 4)
    00468             printf("      R: Consecutive Rs = %3d\n\n", cons);
    00469          if (cons >= parm.redtestinvl && lambda > eps) {
    00470             lastgreeniter = lastyellowiter = lastrediter = iter;
    00471             lambdafactor = 0.67;
    00472             if (parm.printflag & 2)
    00473                printf("\n **** decreasing lamda to %g *****\n\n",
    00474                       lambda * lambdafactor);
    00475          }
    00476          break;
    00477       }
    00478       return lambdafactor;
    00479    }
    00480
    00481    inline void
    00482    print() {
    00483       printf("**** G= %i, Y= %i, R= %i ****\n", ngs, nys, nrs);
    00484       ngs = nrs = nys = 0;
    00485    }
    00486 };
    00487
    00488 //############################################################################
    00489 /* alpha should be decreased if after some number of iterations the objective
    00490    has increased less that 1% */
    00491 class VOL_alpha_factor {
    00492 private:
    00493    VOL_alpha_factor(const VOL_alpha_factor&);
    00494    VOL_alpha_factor& operator=(const VOL_alpha_factor&);
    00495 public:
    00496    double lastvalue;
    00497
    00498    VOL_alpha_factor() {lastvalue = -COIN_DBL_MAX;}
    00499    ~VOL_alpha_factor() {}
    00500
    00501    inline double factor(const VOL_parms& parm, const double lcost,
    00502                         const double alpha) {
    00503       if (alpha < parm.alphamin)
    00504          return 1.0;
    00505       register const double ll = VolAbs(lcost);
    00506       const double x = ll > 10 ? (lcost-lastvalue)/ll : (lcost-lastvalue);
    00507       lastvalue = lcost;
    00508       return (x <= 0.01) ? parm.alphafactor : 1.0;
    00509    }
    00510 };
    00511
    00512 //############################################################################
    00513 /* here we compute the norm of the conjugate direction -hh-, the norm of the
    00514    subgradient -norm-, the inner product between the subgradient and the
    00515    last conjugate direction -vh-, and the inner product between the new
    00516    conjugate direction and the subgradient */
    00517 class VOL_vh {
    00518 private:
    00519    VOL_vh(const VOL_vh&);
    00520    VOL_vh& operator=(const VOL_vh&);
    00521 public:
    00522    double hh;
    00523    double norm;
    00524    double vh;
    00525    double asc;
    00526
    00527    VOL_vh(const double alpha,
    00528           const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
    00529           const VOL_dvector& v, const VOL_dvector& vstar,
    00530           const VOL_dvector& u);
    00531    ~VOL_vh(){}
    00532 };
    00533
    00534 //############################################################################
    00535 /* here we compute different parameter to be printed. v2 is the square of
    00536    the norm of the subgradient. vu is the inner product between the dual
    00537    variables and the subgradient. vabs is the maximum absolute value of
    00538    the violations of pstar. asc is the inner product between the conjugate
    00539    direction and the subgradient */
    00540 class VOL_indc {
    00541 private:
    00542    VOL_indc(const VOL_indc&);
    00543    VOL_indc& operator=(const VOL_indc&);
    00544 public:
    00545    double v2;
    00546    double vu;
    00547    double vabs;
    00548    double asc;
    00549
    00550 public:
    00551    VOL_indc(const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
    00552             const VOL_primal& primal, const VOL_primal& pstar,
    00553             const VOL_dual& dual);
    00554    ~VOL_indc() {}
    00555 };
    00556
    00557 //#############################################################################
    00558
    00565 class VOL_user_hooks {
    00566 public:
    00567    virtual ~VOL_user_hooks() {}
    00568 public:
    00569    // for all hooks: return value of -1 means that volume should quit
    00574    virtual int compute_rc(const VOL_dvector& u, VOL_dvector& rc) = 0;
    00575
    00584    virtual int solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
    00585                                 double& lcost, VOL_dvector& x, VOL_dvector& v,
    00586                                 double& pcost) = 0;
    00593    virtual int heuristics(const VOL_problem& p,
    00594                           const VOL_dvector& x, double& heur_val) = 0;
    00595 };
    00596
    00597 //#############################################################################
    00598
    00607 class VOL_problem {
    00608 private:
    00609    VOL_problem(const VOL_problem&);
    00610    VOL_problem& operator=(const VOL_problem&);
    00611    void set_default_parm();
    00612    // ############ INPUT fields ########################
    00613 public:
    00617    VOL_problem();
    00620    VOL_problem(const char *filename);
    00622    ~VOL_problem();
    00624
    00630    int solve(VOL_user_hooks& hooks, const bool use_preset_dual = false);
    00632
    00633 private:
    00637    double alpha_;
    00639    double lambda_;
    00640    // This union is here for padding (so that data members would be
    00641    // double-aligned on x86 CPU
    00642    union {
    00644       int iter_;
    00645       double __pad0;
    00646    };
    00648
    00649 public:
    00650
    00654    double value;
    00656    VOL_dvector dsol;
    00658    VOL_dvector psol;
    00660    VOL_dvector viol;
    00662
    00666    VOL_parms parm;
    00668    int psize;
    00670    int dsize;
    00673    VOL_dvector dual_lb;
    00676    VOL_dvector dual_ub;
    00678
    00679 public:
    00683    int    iter() const { return iter_; }
    00685    double alpha() const { return alpha_; }
    00687    double lambda() const { return lambda_; }
    00689
    00690 private:
    00694    void read_params(const char* filename);
    00695
    00697    int initialize(const bool use_preset_dual);
    00698
    00700    void print_info(const int iter,
    00701                    const VOL_primal& primal, const VOL_primal& pstar,
    00702                    const VOL_dual& dual);
    00703
    00706    double readjust_target(const double oldtarget, const double lcost) const;
    00707
    00715    double power_heur(const VOL_primal& primal, const VOL_primal& pstar,
    00716                      const VOL_dual& dual) const;
    00718 };
    00719
    00720 #endif

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines