/home/coin/SVN-release/Bcp-1.2.1/Applications/Csp/include/KS.hpp

Go to the documentation of this file.
00001 // Copyright (C) 2005, International Business Machines
00002 // Corporation and others.  All Rights Reserved.
00003 #ifndef _KS_HPP
00004 #define _KS_HPP
00005 
00006 #include <utility>
00007 #include <algorithm>
00008 #include <cfloat>
00009 
00010 using std::fill;
00011 
00012 // A goto-less implementation of the Horowitz-Sahni exact solution 
00013 // procedure for solving knapsack problem.
00014 //
00015 // Reference: Martello and Toth, Knapsack Problems, Wiley, 1990, p30-31.
00016 //
00017 // ToDo: Implement a dynamic programming approach for case
00018 //       of knapsacks with integral coefficients
00019 //-------------------------------------------------------------------
00020 
00021 // The knapsack problem is to find:
00022 
00023 // max {sum(j=1,n) p_j*x_j st. sum (j=1,n)w_j*x_j <= c, x binary}
00024 
00025 // Notation:
00026 //     xhat : current solution vector
00027 //     zhat : current solution value = sum (j=1,n) p_j*xhat_j
00028 //     chat : current residual capacity = c - sum (j=1,n) w_j*xhat_j
00029 //     x    : best solution so far, n-vector.
00030 //     z    : value of the best solution so far =  sum (j=1,n) p_j*x_j
00031 
00032 // Input: n, the number of variables; 
00033 //        c, the rhs;
00034 //        p, n-vector of objective func. coefficients;
00035 //        w, n-vector of the row coeff.
00036 
00037 // Output: z, the optimal objective function value;
00038 //         x, the optimal (binary) solution n-vector;
00039 
00040 struct triplet {
00041    double c;
00042    double w;
00043    int i;
00044 };
00045 
00046 inline bool
00047 ratioComp(const triplet& t0, const triplet& t1) {
00048    return (t0.c/t0.w > t1.c/t1.w ?
00049            true : (t0.c/t0.w < t1.c/t1.w ? false : (t0.i > t1.i)));
00050 }
00051 
00052 class Knapsack {
00053 private:
00054    inline void
00055    sortItems(const double* c, const double* w) {
00056       delete perm_;
00057       double* ww = new double[n_+1];
00058       double* cc = new double[n_+1];
00059       perm_ = new int[n_];
00060       triplet* items = new triplet[n_];
00061       for (int i = n_ - 1; i >= 0; --i) {
00062          items[i].c = c[i];
00063          items[i].w = w[i];
00064          items[i].i = i;
00065       }
00066       std::sort(items, items+n_, ratioComp);
00067       for (int i = n_ - 1; i >= 0; --i) {
00068          cc[i] = items[i].c;
00069          ww[i] = items[i].w;
00070          perm_[i] = items[i].i;
00071       }
00072       delete[] items;
00073       delete[] weight_;
00074       delete[] cost_;
00075       ww[n_] = DBL_MAX;
00076       cc[n_] = 0.0;
00077       weight_ = ww;
00078       cost_ = cc;
00079    }
00080 
00082    double cap_;
00084    int n_;
00087    int* perm_;
00090    const double* cost_;
00093    const double* weight_;
00094 
00096    int* x;
00098    double z;
00099 
00100 public:
00101    Knapsack(int n, double cap, const double* c, const double* w) :
00102       cap_(cap), n_(n), perm_(0), cost_(0), weight_(0), x(new int[n])
00103    {
00104       sortItems(c, w);
00105    }
00106    ~Knapsack() {
00107       delete[] perm_;
00108       delete[] x;
00109       delete[] weight_;
00110       delete[] cost_;
00111    }
00112 
00113    void setCosts(const double* c) {
00114       sortItems(c, weight_);
00115    }
00116    void setCapacity(double cap) { cap_ = cap; }
00117 
00119    void optimize(double lb = 0.0, double minimp = 1e-6);
00120 
00121    const int* getBestSol() const { return x; }
00122    const double getBestVal() const { return z; }
00123 };
00124 
00125 //#############################################################################
00126 
00127 void
00128 Knapsack::optimize(double lb, double minimp)
00129 {
00130    z = lb;
00131    fill(x, x+n_, 0);
00132    int * xhat = new int[n_+1];
00133    fill(xhat, xhat+(n_+1), 0);
00134    int j;
00135 
00136    // set up: adding the extra element and
00137    // accounting for the FORTRAN vs C difference in indexing arrays.
00138    double * p = new double[n_+2];
00139    double * w = new double[n_+2];
00140    int ii;
00141    for (ii = 1; ii < n_+1; ii++){
00142       p[ii]=cost_[ii-1];
00143       w[ii]=weight_[ii-1];
00144    }
00145 
00146    // 1. initialize 
00147    double zhat = 0.0;
00148    z = 0.0;
00149    double caphat = cap_;
00150    p[n_+1] = 0.0;
00151    w[n_+1] = DBL_MAX;
00152    j = 1;
00153 
00154    while (true) {
00155       // 2. compute upper bound u
00156       // "find r = min {i: sum k=j,i w_k>chat};"
00157       ii = j;
00158       double wSemiSum = w[j];
00159       double pSemiSum = p[j];
00160       while (wSemiSum <= caphat && ii < n_+2){
00161          ii++;
00162          wSemiSum += w[ii];
00163          pSemiSum += p[ii];
00164       }
00165       if (ii == n_+2) {
00166          printf("Exceeded iterator limit. Aborting...\n");
00167          abort();
00168       }
00169       // r = ii at this point 
00170       wSemiSum -= w[ii];
00171       pSemiSum -= p[ii];
00172       // *FIXME* : why floor
00173       // double u = pSemiSum + floor((caphat - wSemiSum)*p[ii]/w[ii]);
00174       double u = pSemiSum + (caphat - wSemiSum)*p[ii]/w[ii];
00175       
00176       // "if (z >= zhat + u) goto 5: backtrack;"
00177       if (z + minimp < zhat + u) {
00178          do {
00179             // 3. perform a forward step 
00180             while (w[j] <= caphat){
00181                caphat = caphat - w[j];
00182                zhat = zhat + p[j];
00183                xhat[j] = 1;
00184                ++j;
00185             }
00186             if (j <= n_) {
00187                xhat[j] = 0;
00188                ++j;
00189             }
00190          } while (j == n_); 
00191          
00192          // "if (j<n) goto 2: compute_ub;"
00193          if (j < n_)
00194             continue;
00195       
00196          // 4. up date the best solution so far 
00197          if (zhat > z) {
00198             z = zhat;
00199             int k;
00200             for (k = 0; k < n_; ++k){
00201                x[perm_[k]] = xhat[k+1];
00202             }
00203          }
00204          j = n_;
00205          if (xhat[n_] == 1){
00206             caphat = caphat + w[n_];
00207             zhat = zhat - p[n_];
00208             xhat[n_] = 0;
00209          }
00210       }
00211       // 5. backtrack 
00212       // "find i=max{k<j:xhat[k]=1};"
00213       int i = j-1; 
00214       while (!(xhat[i]==1) && i>0) {
00215          i--;
00216       }
00217       
00218       // "if (no such i exists) return;"
00219       if (i == 0) {
00220          delete [] p;
00221          delete [] w;
00222          delete [] xhat;
00223          break;
00224       }
00225 
00226       caphat = caphat + w[i];
00227       zhat = zhat - p[i];
00228       xhat[i] = 0;
00229       j = i+1;
00230       // "goto 2: compute_ub;"
00231    }
00232 }
00233 
00234 #endif

Generated on Thu Jan 15 03:00:58 2009 for coin-Bcp by  doxygen 1.4.7