00001
00002
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
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
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
00137
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
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
00156
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
00170 wSemiSum -= w[ii];
00171 pSemiSum -= p[ii];
00172
00173
00174 double u = pSemiSum + (caphat - wSemiSum)*p[ii]/w[ii];
00175
00176
00177 if (z + minimp < zhat + u) {
00178 do {
00179
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
00193 if (j < n_)
00194 continue;
00195
00196
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
00212
00213 int i = j-1;
00214 while (!(xhat[i]==1) && i>0) {
00215 i--;
00216 }
00217
00218
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
00231 }
00232 }
00233
00234 #endif