VolVolume.hpp

Go to the documentation of this file.

// $Id: VolVolume.hpp 375 2013-11-22 15:46:16Z tkr $
// Copyright (C) 2000, International Business Machines
// Corporation and others.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).

#ifndef __VOLUME_HPP__
#define __VOLUME_HPP__

#include <cfloat>
#include <algorithm>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <cstdlib>

#ifndef VOL_DEBUG
// When VOL_DEBUG is 1, we check vector indices
#define VOL_DEBUG 0
#endif

template <class T> static inline T
VolMax(register const T x, register const T y) {
   return ((x) > (y)) ? (x) : (y);
}

template <class T> static inline T
VolAbs(register const T x) {
   return ((x) > 0) ? (x) : -(x);
}

//############################################################################

#if defined(VOL_DEBUG) && (VOL_DEBUG != 0)
#define VOL_TEST_INDEX(i, size)                 \
{                                               \
   if ((i) < 0 || (i) >= (size)) {              \
      printf("bad VOL_?vector index\n");        \
      abort();                                  \
   }                                            \
}
#define VOL_TEST_SIZE(size)                     \
{                                               \
   if (s <= 0) {                                \
      printf("bad VOL_?vector size\n");         \
      abort();                                  \
   }                                            \
}
#else
#define VOL_TEST_INDEX(i, size)
#define VOL_TEST_SIZE(size)
#endif
      
//############################################################################

class VOL_dvector;
class VOL_ivector;
class VOL_primal;
class VOL_dual;
class VOL_swing;
class VOL_alpha_factor;
class VOL_vh;
class VOL_indc;
class VOL_user_hooks;
class VOL_problem;

//############################################################################

struct VOL_parms {
   double lambdainit; 
   double alphainit; 
   double alphamin;
   double alphafactor;

   double ubinit;

   double primal_abs_precision;
   double gap_abs_precision; 
   double gap_rel_precision;  
   double granularity;

   double minimum_rel_ascent;
   int    ascent_first_check;
   int    ascent_check_invl;
   
   int    maxsgriters; 

   int    printflag; 
   int    printinvl; 
   int    heurinvl; 

   int greentestinvl; 
   int yellowtestinvl; 
   int redtestinvl;

   int    alphaint; 

   char* temp_dualfile;
};

//############################################################################

class VOL_dvector {
public:
   double* v;
   int sz;

public:
   VOL_dvector(const int s) {
      VOL_TEST_SIZE(s);
      v = new double[sz = s];
   }
   VOL_dvector() : v(0), sz(0) {}
   VOL_dvector(const VOL_dvector& x) : v(0), sz(0) {
      sz = x.sz;
      if (sz > 0) {
        v = new double[sz];
        std::memcpy(v, x.v, sz * sizeof(double));
      }
   }
   ~VOL_dvector() { delete[] v; }

   inline int size() const {return sz;}

   inline double& operator[](const int i) {
      VOL_TEST_INDEX(i, sz);
      return v[i];
   }

   inline double operator[](const int i) const {
      VOL_TEST_INDEX(i, sz);
      return v[i];
   }

   inline void clear() {
      delete[] v;
      v = 0;
      sz = 0;
   }
   inline void cc(const double gamma, const VOL_dvector& w) {
      if (sz != w.sz) {
         printf("bad VOL_dvector sizes\n");
         abort();
      }
      double * p_v = v - 1;
      const double * p_w = w.v - 1;
      const double * const p_e = v + sz;
      const double one_gamma = 1.0 - gamma;
      while ( ++p_v != p_e ){
         *p_v = one_gamma * (*p_v) + gamma * (*++p_w);
      }
   }

   inline void allocate(const int s) {
      VOL_TEST_SIZE(s);
      delete[] v;                 
      v = new double[sz = s];
   }

   inline void swap(VOL_dvector& w) {
      std::swap(v, w.v);
      std::swap(sz, w.sz);
   }

   VOL_dvector& operator=(const VOL_dvector& w);
   VOL_dvector& operator=(const double w);
};

//-----------------------------------------------------------------------------
class VOL_ivector {
public:
   int* v;
   int sz;
public:
   VOL_ivector(const int s) {
      VOL_TEST_SIZE(s);
      v = new int[sz = s];
   }
   VOL_ivector() : v(0), sz(0) {}
   VOL_ivector(const VOL_ivector& x) {
      sz = x.sz;
      if (sz > 0) {
         v = new int[sz];
         std::memcpy(v, x.v, sz * sizeof(int));
      }
   }
   ~VOL_ivector(){
      delete [] v;
   }

   inline int size() const { return sz; }
   inline int& operator[](const int i) {
      VOL_TEST_INDEX(i, sz);
      return v[i];
   }

   inline int operator[](const int i) const {
      VOL_TEST_INDEX(i, sz);
      return v[i];
   }

   inline void clear() {
      delete[] v;
      v = 0;
      sz = 0;
   }

   inline void allocate(const int s) {
      VOL_TEST_SIZE(s);
      delete[] v;
      v = new int[sz = s];
   }

   inline void swap(VOL_ivector& w) {
      std::swap(v, w.v);
      std::swap(sz, w.sz);
   }

   VOL_ivector& operator=(const VOL_ivector& v);      
   VOL_ivector& operator=(const int w);
};

//############################################################################
// A class describing a primal solution. This class is used only internally 
class VOL_primal {
public: 
   // objective value of this primal solution 
   double value;
   // the largest of the v[i]'s
   double viol;  
   // primal solution  
   VOL_dvector x;
   // v=b-Ax, for the relaxed constraints
   VOL_dvector v; 

   VOL_primal(const int psize, const int dsize) : x(psize), v(dsize) {}
   VOL_primal(const VOL_primal& primal) :
      value(primal.value), viol(primal.viol), x(primal.x), v(primal.v) {}
   ~VOL_primal() {}
   inline VOL_primal& operator=(const VOL_primal& p) {
      if (this == &p) 
         return *this;
      value = p.value;
      viol = p.viol;
      x = p.x;
      v = p.v;
      return *this;
   }

   // convex combination. data members in this will be overwritten
   // convex combination between two primal solutions
   // x <-- alpha x + (1 - alpha) p.x
   // v <-- alpha v + (1 - alpha) p.v
   inline void cc(const double alpha, const VOL_primal& p) {
      value = alpha * p.value + (1.0 - alpha) * value;
      x.cc(alpha, p.x);
      v.cc(alpha, p.v);
   }
   // find maximum of v[i]
   void find_max_viol(const VOL_dvector& dual_lb, 
                      const VOL_dvector& dual_ub);
};

//-----------------------------------------------------------------------------
// A class describing a dual solution. This class is used only internally 
class VOL_dual {
public:
   // lagrangian value
   double lcost; 
   // reduced costs * (pstar-primal)
   double xrc;
   // this information is only printed
   // dual vector
   VOL_dvector u; 

   VOL_dual(const int dsize) : u(dsize) { u = 0.0;}
   VOL_dual(const VOL_dual& dual) :
      lcost(dual.lcost), xrc(dual.xrc), u(dual.u) {}
   ~VOL_dual() {}
   inline VOL_dual& operator=(const VOL_dual& p) {
      if (this == &p) 
         return *this;
      lcost = p.lcost;
      xrc = p.xrc;
      u = p.u;
      return *this;
   }
   // dual step
   void   step(const double target, const double lambda,
               const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
               const VOL_dvector& v);
   double ascent(const VOL_dvector& v, const VOL_dvector& last_u) const;
   void   compute_xrc(const VOL_dvector& pstarx, const VOL_dvector& primalx,
                      const VOL_dvector& rc);

};


//############################################################################
/* here we check whether an iteration is green, yellow or red. Also according
   to this information we decide whether lambda should be changed */
class VOL_swing {
private:
   VOL_swing(const VOL_swing&);
   VOL_swing& operator=(const VOL_swing&);
public:
   enum condition {green, yellow, red} lastswing;
   int lastgreeniter, lastyellowiter, lastrediter;
   int ngs, nrs, nys;
   int rd;
   
   VOL_swing() {
      lastgreeniter = lastyellowiter = lastrediter = 0;
      ngs = nrs = nys = 0;
   }
   ~VOL_swing(){}

   inline void cond(const VOL_dual& dual, 
                    const double lcost, const double ascent, const int iter) {
      double eps = 1.e-3;

      if (ascent > 0.0  &&  lcost > dual.lcost + eps) {
         lastswing = green;
         lastgreeniter = iter;
         ++ngs;
         rd = 0;
      } else { 
         if (ascent <= 0  &&  lcost > dual.lcost) {
            lastswing = yellow;
            lastyellowiter = iter;
            ++nys;
            rd = 0;
         } else {
            lastswing = red;
            lastrediter = iter;
            ++nrs;
            rd = 1;
         }
      }
   }

   inline double
   lfactor(const VOL_parms& parm, const double lambda, const int iter) {
      double lambdafactor = 1.0;
      double eps = 5.e-4;
      int cons;

      switch (lastswing) {
       case green:
         cons = iter - VolMax(lastyellowiter, lastrediter);
         if (parm.printflag & 4)
            printf("      G: Consecutive Gs = %3d\n\n", cons);
         if (cons >= parm.greentestinvl && lambda < 2.0) {
            lastgreeniter = lastyellowiter = lastrediter = iter;
            lambdafactor = 2.0;
            if (parm.printflag & 2)
               printf("\n ---- increasing lamda to %g ----\n\n",
                      lambda * lambdafactor); 
         }
         break;
      
       case yellow:
         cons = iter - VolMax(lastgreeniter, lastrediter);
         if (parm.printflag & 4)
            printf("      Y: Consecutive Ys = %3d\n\n", cons);
         if (cons >= parm.yellowtestinvl) {
            lastgreeniter = lastyellowiter = lastrediter = iter;
            lambdafactor = 1.1;
            if (parm.printflag & 2)
               printf("\n **** increasing lamda to %g *****\n\n",
                      lambda * lambdafactor);
         }
         break;
      
       case red:
         cons = iter - VolMax(lastgreeniter, lastyellowiter);
         if (parm.printflag & 4)
            printf("      R: Consecutive Rs = %3d\n\n", cons);
         if (cons >= parm.redtestinvl && lambda > eps) {
            lastgreeniter = lastyellowiter = lastrediter = iter;
            lambdafactor = 0.67;
            if (parm.printflag & 2)
               printf("\n **** decreasing lamda to %g *****\n\n",
                      lambda * lambdafactor);
         } 
         break;
      }
      return lambdafactor;
   }

   inline void
   print() {
      printf("**** G= %i, Y= %i, R= %i ****\n", ngs, nys, nrs);
      ngs = nrs = nys = 0;  
   }
};

//############################################################################
/* alpha should be decreased if after some number of iterations the objective
   has increased less that 1% */
class VOL_alpha_factor {
private:
   VOL_alpha_factor(const VOL_alpha_factor&);
   VOL_alpha_factor& operator=(const VOL_alpha_factor&);
public:
   double lastvalue;

   VOL_alpha_factor() {lastvalue = -DBL_MAX;}
   ~VOL_alpha_factor() {}

   inline double factor(const VOL_parms& parm, const double lcost,
                        const double alpha) {
      if (alpha < parm.alphamin)
         return 1.0;
      register const double ll = VolAbs(lcost);
      const double x = ll > 10 ? (lcost-lastvalue)/ll : (lcost-lastvalue);
      lastvalue = lcost;
      return (x <= 0.01) ? parm.alphafactor : 1.0;
   }
};

//############################################################################
/* here we compute the norm of the conjugate direction -hh-, the norm of the
   subgradient -norm-, the inner product between the subgradient and the 
   last conjugate direction -vh-, and the inner product between the new
   conjugate direction and the subgradient */
class VOL_vh {
private:
   VOL_vh(const VOL_vh&);
   VOL_vh& operator=(const VOL_vh&);
public:
   double hh;
   double norm;
   double vh;
   double asc;

   VOL_vh(const double alpha,
          const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
          const VOL_dvector& v, const VOL_dvector& vstar,
          const VOL_dvector& u);
   ~VOL_vh(){}
};

//############################################################################
/* here we compute different parameter to be printed. v2 is the square of 
   the norm of the subgradient. vu is the inner product between the dual
   variables and the subgradient. vabs is the maximum absolute value of
   the violations of pstar. asc is the inner product between the conjugate
   direction and the subgradient */
class VOL_indc {
private:
   VOL_indc(const VOL_indc&);
   VOL_indc& operator=(const VOL_indc&);
public:
   double v2;
   double vu;
   double vabs;
   double asc;

public:
   VOL_indc(const VOL_dvector& dual_lb, const VOL_dvector& dual_ub,
            const VOL_primal& primal, const VOL_primal& pstar,
            const VOL_dual& dual);
   ~VOL_indc() {}
};

//#############################################################################

class VOL_user_hooks {
public:
   virtual ~VOL_user_hooks() {}
public:
   // for all hooks: return value of -1 means that volume should quit
   virtual int compute_rc(const VOL_dvector& u, VOL_dvector& rc) = 0;

   virtual int solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
                                double& lcost, VOL_dvector& x, VOL_dvector& v,
                                double& pcost) = 0;
   virtual int heuristics(const VOL_problem& p, 
                          const VOL_dvector& x, double& heur_val) = 0;
};

//#############################################################################

class VOL_problem {
private:
   VOL_problem(const VOL_problem&);
   VOL_problem& operator=(const VOL_problem&);
   void set_default_parm();
   // ############ INPUT fields ########################
public: 
   VOL_problem();
   VOL_problem(const char *filename);
   ~VOL_problem();

   int solve(VOL_user_hooks& hooks, const bool use_preset_dual = false);

private: 
   double alpha_; 
   double lambda_;
   // This union is here for padding (so that data members would be
   // double-aligned on x86 CPU
   union {
      int iter_;
      double __pad0;
   };

public:
  
   double value;
   VOL_dvector dsol;
   VOL_dvector psol;
   VOL_dvector viol;

   VOL_parms parm;
   int psize;        
   int dsize;      
   VOL_dvector dual_lb;
   VOL_dvector dual_ub;

public:
   int    iter() const { return iter_; }
   double alpha() const { return alpha_; }
   double lambda() const { return lambda_; }

private:
   void read_params(const char* filename);

   int initialize(const bool use_preset_dual);

   void print_info(const int iter,
                   const VOL_primal& primal, const VOL_primal& pstar,
                   const VOL_dual& dual);

   double readjust_target(const double oldtarget, const double lcost) const;

   double power_heur(const VOL_primal& primal, const VOL_primal& pstar,
                     const VOL_dual& dual) const;
};

#endif