00001
00002
00003 #if defined(_MSC_VER)
00004
00005 # pragma warning(disable:4786)
00006 #endif
00007 #include <cassert>
00008 #include <cmath>
00009 #include <cfloat>
00010
00011 #include "OsiSolverInterface.hpp"
00012 #include "SbbModel.hpp"
00013 #include "SbbMessage.hpp"
00014 #include "SbbHeuristic.hpp"
00015
00016
00017 SbbHeuristic::SbbHeuristic()
00018 :model_(NULL)
00019 {
00020 }
00021
00022
00023 SbbHeuristic::SbbHeuristic(SbbModel & model)
00024 :
00025 model_(&model)
00026 {
00027 }
00028
00029
00030 SbbHeuristic::~SbbHeuristic ()
00031 {
00032 }
00033
00034
00035 void SbbHeuristic::setModel(SbbModel * model)
00036 {
00037 model_ = model;
00038 }
00039
00040
00041 SbbRounding::SbbRounding()
00042 :SbbHeuristic()
00043 {
00044
00045 }
00046
00047
00048 SbbRounding::SbbRounding(SbbModel & model)
00049 :SbbHeuristic(model)
00050 {
00051
00052 assert(model.solver());
00053 matrix_ = *model.solver()->getMatrixByCol();
00054 matrixByRow_ = *model.solver()->getMatrixByRow();
00055 seed_=1;
00056 }
00057
00058
00059 SbbRounding::~SbbRounding ()
00060 {
00061 }
00062
00063
00064 SbbHeuristic *
00065 SbbRounding::clone() const
00066 {
00067 return new SbbRounding(*this);
00068 }
00069
00070
00071 SbbRounding::SbbRounding(const SbbRounding & rhs)
00072 :
00073 SbbHeuristic(rhs),
00074 matrix_(rhs.matrix_),
00075 matrixByRow_(rhs.matrixByRow_),
00076 seed_(rhs.seed_)
00077 {
00078 }
00079
00080
00081
00082
00083
00084
00085 int
00086 SbbRounding::solution(double & solutionValue,
00087 double * betterSolution)
00088 {
00089
00090 OsiSolverInterface * solver = model_->solver();
00091 const double * lower = solver->getColLower();
00092 const double * upper = solver->getColUpper();
00093 const double * rowLower = solver->getRowLower();
00094 const double * rowUpper = solver->getRowUpper();
00095 const double * solution = solver->getColSolution();
00096 const double * objective = solver->getObjCoefficients();
00097 double integerTolerance = model_->getDblParam(SbbModel::SbbIntegerTolerance);
00098 double primalTolerance;
00099 solver->getDblParam(OsiPrimalTolerance,primalTolerance);
00100
00101 int numberRows = matrix_.getNumRows();
00102
00103 int numberIntegers = model_->numberIntegers();
00104 const int * integerVariable = model_->integerVariable();
00105 int i;
00106 double direction = solver->getObjSense();
00107 double newSolutionValue = direction*solver->getObjValue();
00108 int returnCode = 0;
00109
00110
00111 const double * element = matrix_.getElements();
00112 const int * row = matrix_.getIndices();
00113 const int * columnStart = matrix_.getVectorStarts();
00114 const int * columnLength = matrix_.getVectorLengths();
00115
00116 const double * elementByRow = matrixByRow_.getElements();
00117 const int * column = matrixByRow_.getIndices();
00118 const int * rowStart = matrixByRow_.getVectorStarts();
00119 const int * rowLength = matrixByRow_.getVectorLengths();
00120
00121
00122 int numberColumns = solver->getNumCols();
00123 double * newSolution = new double [numberColumns];
00124 memcpy(newSolution,solution,numberColumns*sizeof(double));
00125
00126 double * rowActivity = new double[numberRows];
00127 memset(rowActivity,0,numberRows*sizeof(double));
00128 for (i=0;i<numberColumns;i++) {
00129 int j;
00130 double value = newSolution[i];
00131 if (value) {
00132 for (j=columnStart[i];
00133 j<columnStart[i]+columnLength[i];j++) {
00134 int iRow=row[j];
00135 rowActivity[iRow] += value*element[j];
00136 }
00137 }
00138 }
00139
00140 for (i=0;i<numberRows;i++) {
00141 if(rowActivity[i]<rowLower[i]) {
00142
00143 rowActivity[i]=rowLower[i];
00144 } else if(rowActivity[i]>rowUpper[i]) {
00145
00146 rowActivity[i]=rowUpper[i];
00147 }
00148 }
00149 for (i=0;i<numberIntegers;i++) {
00150 int iColumn = integerVariable[i];
00151 double value=newSolution[iColumn];
00152 if (fabs(floor(value+0.5)-value)>integerTolerance) {
00153 double below = floor(value);
00154 double newValue=newSolution[iColumn];
00155 double cost = direction * objective[iColumn];
00156 double move;
00157 if (cost>0.0) {
00158
00159 move = 1.0 -(value-below);
00160 } else if (cost<0.0) {
00161
00162 move = below-value;
00163 } else {
00164
00165
00166 move = below-value;
00167 }
00168 newValue += move;
00169 newSolution[iColumn] = newValue;
00170 newSolutionValue += move*cost;
00171 int j;
00172 for (j=columnStart[iColumn];
00173 j<columnStart[iColumn]+columnLength[iColumn];j++) {
00174 int iRow = row[j];
00175 rowActivity[iRow] += move*element[j];
00176 }
00177 }
00178 }
00179
00180 double penalty=0.0;
00181
00182
00183 for (i=0;i<numberRows;i++) {
00184 double value = rowActivity[i];
00185 double thisInfeasibility=0.0;
00186 if (value<rowLower[i]-primalTolerance)
00187 thisInfeasibility = value-rowLower[i];
00188 else if (value>rowUpper[i]+primalTolerance)
00189 thisInfeasibility = value-rowUpper[i];
00190 if (thisInfeasibility) {
00191
00192
00193 double bestCost = 1.0e50;
00194 int k;
00195 int iBest=-1;
00196 double addCost=0.0;
00197 double newValue=0.0;
00198 double changeRowActivity=0.0;
00199 double absInfeasibility = fabs(thisInfeasibility);
00200 for (k=rowStart[i];k<rowStart[i]+rowLength[i];k++) {
00201 int iColumn = column[k];
00202 if (columnLength[iColumn]==1) {
00203 double currentValue = newSolution[iColumn];
00204 double elementValue = elementByRow[k];
00205 double lowerValue = lower[iColumn];
00206 double upperValue = upper[iColumn];
00207 double gap = rowUpper[i]-rowLower[i];
00208 double absElement=fabs(elementValue);
00209 if (thisInfeasibility*elementValue>0.0) {
00210
00211 if ((currentValue-lowerValue)*absElement>=absInfeasibility) {
00212
00213 double distance = absInfeasibility/absElement;
00214 double thisCost = -direction*objective[iColumn]*distance;
00215 if (solver->isInteger(iColumn)) {
00216 distance = ceil(distance-primalTolerance);
00217 assert (currentValue-distance>=lowerValue-primalTolerance);
00218 if (absInfeasibility-distance*absElement< -gap-primalTolerance)
00219 thisCost=1.0e100;
00220 else
00221 thisCost = -direction*objective[iColumn]*distance;
00222 }
00223 if (thisCost<bestCost) {
00224 bestCost=thisCost;
00225 iBest=iColumn;
00226 addCost = thisCost;
00227 newValue = currentValue-distance;
00228 changeRowActivity = -distance*elementValue;
00229 }
00230 }
00231 } else {
00232
00233 if ((upperValue-currentValue)*absElement>=absInfeasibility) {
00234
00235 double distance = absInfeasibility/absElement;
00236 double thisCost = direction*objective[iColumn]*distance;
00237 if (solver->isInteger(iColumn)) {
00238 distance = ceil(distance-1.0e-7);
00239 assert (currentValue-distance<=upperValue+primalTolerance);
00240 if (absInfeasibility-distance*absElement< -gap-primalTolerance)
00241 thisCost=1.0e100;
00242 else
00243 thisCost = direction*objective[iColumn]*distance;
00244 }
00245 if (thisCost<bestCost) {
00246 bestCost=thisCost;
00247 iBest=iColumn;
00248 addCost = thisCost;
00249 newValue = currentValue+distance;
00250 changeRowActivity = distance*elementValue;
00251 }
00252 }
00253 }
00254 }
00255 }
00256 if (iBest>=0) {
00257
00258
00259 newSolution[iBest]=newValue;
00260 thisInfeasibility=0.0;
00261 newSolutionValue += addCost;
00262 rowActivity[i] += changeRowActivity;
00263 }
00264 penalty += fabs(thisInfeasibility);
00265 }
00266 }
00267
00268
00269 if (!penalty) {
00270
00271
00272
00273
00274 double randomNumber = CoinDrand48();
00275 int iPass;
00276 int start[2];
00277 int end[2];
00278 int iRandom = (int) (randomNumber*((double) numberIntegers));
00279 start[0]=iRandom;
00280 end[0]=numberIntegers;
00281 start[1]=0;
00282 end[1]=iRandom;
00283 for (iPass=0;iPass<2;iPass++) {
00284 int i;
00285 for (i=start[iPass];i<end[iPass];i++) {
00286 int iColumn = integerVariable[i];
00287 double value=newSolution[iColumn];
00288 assert (fabs(floor(value+0.5)-value)<integerTolerance);
00289 double cost = direction * objective[iColumn];
00290 double move=0.0;
00291 if (cost>0.0)
00292 move = -1.0;
00293 else if (cost<0.0)
00294 move=1.0;
00295 while (move) {
00296 bool good=true;
00297 double newValue=newSolution[iColumn]+move;
00298 if (newValue<lower[iColumn]-primalTolerance||
00299 newValue>upper[iColumn]+primalTolerance) {
00300 move=0.0;
00301 } else {
00302
00303 int j;
00304 for (j=columnStart[iColumn];
00305 j<columnStart[iColumn]+columnLength[iColumn];j++) {
00306 int iRow = row[j];
00307 double newActivity = rowActivity[iRow] + move*element[j];
00308 if (newActivity<rowLower[iRow]-primalTolerance||
00309 newActivity>rowUpper[iRow]+primalTolerance) {
00310 good=false;
00311 break;
00312 }
00313 }
00314 if (good) {
00315 newSolution[iColumn] = newValue;
00316 newSolutionValue += move*cost;
00317 int j;
00318 for (j=columnStart[iColumn];
00319 j<columnStart[iColumn]+columnLength[iColumn];j++) {
00320 int iRow = row[j];
00321 rowActivity[iRow] += move*element[j];
00322 }
00323 } else {
00324 move=0.0;
00325 }
00326 }
00327 }
00328 }
00329 }
00330 if (newSolutionValue<solutionValue) {
00331
00332 memset(rowActivity,0,numberRows*sizeof(double));
00333 for (i=0;i<numberColumns;i++) {
00334 int j;
00335 double value = newSolution[i];
00336 if (value) {
00337 for (j=columnStart[i];
00338 j<columnStart[i]+columnLength[i];j++) {
00339 int iRow=row[j];
00340 rowActivity[iRow] += value*element[j];
00341 }
00342 }
00343 }
00344
00345 bool feasible=true;
00346 for (i=0;i<numberRows;i++) {
00347 if(rowActivity[i]<rowLower[i]) {
00348 if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
00349 feasible = false;
00350 } else if(rowActivity[i]>rowUpper[i]) {
00351 if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
00352 feasible = false;
00353 }
00354 }
00355 if (feasible) {
00356
00357 memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
00358 solutionValue = newSolutionValue;
00359
00360 returnCode=1;
00361 } else {
00362
00363
00364 }
00365 }
00366 }
00367 delete [] newSolution;
00368 delete [] rowActivity;
00369 return returnCode;
00370 }
00371
00372 void SbbRounding::setModel(SbbModel * model)
00373 {
00374 model_ = model;
00375
00376 assert(model_->solver());
00377 matrix_ = *model_->solver()->getMatrixByCol();
00378 matrixByRow_ = *model_->solver()->getMatrixByRow();
00379 }
00380
00381