Changeset 2385 for trunk/Clp/src/ClpNonLinearCost.cpp
 Timestamp:
 Jan 6, 2019 2:43:06 PM (4 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Clp/src/ClpNonLinearCost.cpp
r2235 r2385 21 21 // Default Constructor 22 22 // 23 ClpNonLinearCost::ClpNonLinearCost () : 24 changeCost_(0.0), 25 feasibleCost_(0.0), 26 infeasibilityWeight_(1.0), 27 largestInfeasibility_(0.0), 28 sumInfeasibilities_(0.0), 29 averageTheta_(0.0), 30 numberRows_(0), 31 numberColumns_(0), 32 start_(NULL), 33 whichRange_(NULL), 34 offset_(NULL), 35 lower_(NULL), 36 cost_(NULL), 37 model_(NULL), 38 infeasible_(NULL), 39 numberInfeasibilities_(1), 40 status_(NULL), 41 bound_(NULL), 42 cost2_(NULL), 43 method_(1), 44 convex_(true), 45 bothWays_(false) 46 { 47 23 ClpNonLinearCost::ClpNonLinearCost() 24 : changeCost_(0.0) 25 , feasibleCost_(0.0) 26 , infeasibilityWeight_(1.0) 27 , largestInfeasibility_(0.0) 28 , sumInfeasibilities_(0.0) 29 , averageTheta_(0.0) 30 , numberRows_(0) 31 , numberColumns_(0) 32 , start_(NULL) 33 , whichRange_(NULL) 34 , offset_(NULL) 35 , lower_(NULL) 36 , cost_(NULL) 37 , model_(NULL) 38 , infeasible_(NULL) 39 , numberInfeasibilities_(1) 40 , status_(NULL) 41 , bound_(NULL) 42 , cost2_(NULL) 43 , method_(1) 44 , convex_(true) 45 , bothWays_(false) 46 { 48 47 } 49 48 //#define VALIDATE 50 49 #ifdef VALIDATE 51 static double * 52 static double * 50 static double *saveLowerV = NULL; 51 static double *saveUpperV = NULL; 53 52 #ifdef NDEBUG 54 53 Validate sgould not be set if no debug … … 59 58 later may do dual analysis and even finding duplicate columns 60 59 */ 61 ClpNonLinearCost::ClpNonLinearCost ( ClpSimplex * model, int method) 62 { 63 method = 2; 64 model_ = model; 65 numberRows_ = model_>numberRows(); 66 //if (numberRows_==402) { 67 //model_>setLogLevel(63); 68 //model_>setMaximumIterations(30000); 69 //} 70 numberColumns_ = model_>numberColumns(); 71 // If gub then we need this extra 72 int numberExtra = model_>numberExtraRows(); 73 if (numberExtra) 74 method = 1; 75 int numberTotal1 = numberRows_ + numberColumns_; 76 int numberTotal = numberTotal1 + numberExtra; 77 convex_ = true; 78 bothWays_ = false; 79 method_ = method; 80 numberInfeasibilities_ = 0; 81 changeCost_ = 0.0; 82 feasibleCost_ = 0.0; 83 infeasibilityWeight_ = 1.0; 84 double * cost = model_>costRegion(); 85 // check if all 0 86 int iSequence; 87 bool allZero = true; 88 for (iSequence = 0; iSequence < numberTotal1; iSequence++) { 89 if (cost[iSequence]) { 90 allZero = false; 91 break; 92 } 93 } 94 if (allZero&&model_>clpMatrix()>type()<15) 95 model_>setInfeasibilityCost(1.0); 96 double infeasibilityCost = model_>infeasibilityCost(); 97 sumInfeasibilities_ = 0.0; 98 averageTheta_ = 0.0; 99 largestInfeasibility_ = 0.0; 100 // All arrays NULL to start 101 status_ = NULL; 102 bound_ = NULL; 103 cost2_ = NULL; 104 start_ = NULL; 105 whichRange_ = NULL; 106 offset_ = NULL; 107 lower_ = NULL; 108 cost_ = NULL; 109 infeasible_ = NULL; 110 111 double * upper = model_>upperRegion(); 112 double * lower = model_>lowerRegion(); 113 114 // See how we are storing things 115 bool always4 = (model_>clpMatrix()> 116 generalExpanded(model_, 10, iSequence) != 0); 117 if (always4) 118 method_ = 1; 119 if (CLP_METHOD1) { 120 start_ = new int [numberTotal+1]; 121 whichRange_ = new int [numberTotal]; 122 offset_ = new int [numberTotal]; 123 memset(offset_, 0, numberTotal * sizeof(int)); 124 125 126 // First see how much space we need 127 int put = 0; 128 129 // For quadratic we need inf,0,0,+inf 130 for (iSequence = 0; iSequence < numberTotal1; iSequence++) { 131 if (!always4) { 132 if (lower[iSequence] > COIN_DBL_MAX) 133 put++; 134 if (upper[iSequence] < COIN_DBL_MAX) 135 put++; 136 put += 2; 137 } else { 138 put += 4; 139 } 140 } 141 142 // and for extra 143 put += 4 * numberExtra; 60 ClpNonLinearCost::ClpNonLinearCost(ClpSimplex * model, int method) 61 { 62 method = 2; 63 model_ = model; 64 numberRows_ = model_>numberRows(); 65 //if (numberRows_==402) { 66 //model_>setLogLevel(63); 67 //model_>setMaximumIterations(30000); 68 //} 69 numberColumns_ = model_>numberColumns(); 70 // If gub then we need this extra 71 int numberExtra = model_>numberExtraRows(); 72 if (numberExtra) 73 method = 1; 74 int numberTotal1 = numberRows_ + numberColumns_; 75 int numberTotal = numberTotal1 + numberExtra; 76 convex_ = true; 77 bothWays_ = false; 78 method_ = method; 79 numberInfeasibilities_ = 0; 80 changeCost_ = 0.0; 81 feasibleCost_ = 0.0; 82 infeasibilityWeight_ = 1.0; 83 double *cost = model_>costRegion(); 84 // check if all 0 85 int iSequence; 86 bool allZero = true; 87 for (iSequence = 0; iSequence < numberTotal1; iSequence++) { 88 if (cost[iSequence]) { 89 allZero = false; 90 break; 91 } 92 } 93 if (allZero && model_>clpMatrix()>type() < 15) 94 model_>setInfeasibilityCost(1.0); 95 double infeasibilityCost = model_>infeasibilityCost(); 96 sumInfeasibilities_ = 0.0; 97 averageTheta_ = 0.0; 98 largestInfeasibility_ = 0.0; 99 // All arrays NULL to start 100 status_ = NULL; 101 bound_ = NULL; 102 cost2_ = NULL; 103 start_ = NULL; 104 whichRange_ = NULL; 105 offset_ = NULL; 106 lower_ = NULL; 107 cost_ = NULL; 108 infeasible_ = NULL; 109 110 double *upper = model_>upperRegion(); 111 double *lower = model_>lowerRegion(); 112 113 // See how we are storing things 114 bool always4 = (model_>clpMatrix()>generalExpanded(model_, 10, iSequence) != 0); 115 if (always4) 116 method_ = 1; 117 if (CLP_METHOD1) { 118 start_ = new int[numberTotal + 1]; 119 whichRange_ = new int[numberTotal]; 120 offset_ = new int[numberTotal]; 121 memset(offset_, 0, numberTotal * sizeof(int)); 122 123 // First see how much space we need 124 int put = 0; 125 126 // For quadratic we need inf,0,0,+inf 127 for (iSequence = 0; iSequence < numberTotal1; iSequence++) { 128 if (!always4) { 129 if (lower[iSequence] > COIN_DBL_MAX) 130 put++; 131 if (upper[iSequence] < COIN_DBL_MAX) 132 put++; 133 put += 2; 134 } else { 135 put += 4; 136 } 137 } 138 139 // and for extra 140 put += 4 * numberExtra; 144 141 #ifndef NDEBUG 145 146 #endif 147 lower_ = new double[put];148 cost_ = new double[put];149 infeasible_ = new unsigned int[(put+31)>>5];150 memset(infeasible_, 0, ((put + 31) >> 5)*sizeof(unsigned int));151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 start_[iSequence+1] = put;187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 start_[iSequence+1] = put;202 203 assert(put <= kPut);204 142 int kPut = put; 143 #endif 144 lower_ = new double[put]; 145 cost_ = new double[put]; 146 infeasible_ = new unsigned int[(put + 31) >> 5]; 147 memset(infeasible_, 0, ((put + 31) >> 5) * sizeof(unsigned int)); 148 149 put = 0; 150 151 start_[0] = 0; 152 153 for (iSequence = 0; iSequence < numberTotal1; iSequence++) { 154 if (!always4) { 155 if (lower[iSequence] > COIN_DBL_MAX) { 156 lower_[put] = COIN_DBL_MAX; 157 setInfeasible(put, true); 158 cost_[put++] = cost[iSequence]  infeasibilityCost; 159 } 160 whichRange_[iSequence] = put; 161 lower_[put] = lower[iSequence]; 162 cost_[put++] = cost[iSequence]; 163 lower_[put] = upper[iSequence]; 164 cost_[put++] = cost[iSequence] + infeasibilityCost; 165 if (upper[iSequence] < COIN_DBL_MAX) { 166 lower_[put] = COIN_DBL_MAX; 167 setInfeasible(put  1, true); 168 cost_[put++] = 1.0e50; 169 } 170 } else { 171 lower_[put] = COIN_DBL_MAX; 172 setInfeasible(put, true); 173 cost_[put++] = cost[iSequence]  infeasibilityCost; 174 whichRange_[iSequence] = put; 175 lower_[put] = lower[iSequence]; 176 cost_[put++] = cost[iSequence]; 177 lower_[put] = upper[iSequence]; 178 cost_[put++] = cost[iSequence] + infeasibilityCost; 179 lower_[put] = COIN_DBL_MAX; 180 setInfeasible(put  1, true); 181 cost_[put++] = 1.0e50; 182 } 183 start_[iSequence + 1] = put; 184 } 185 for (; iSequence < numberTotal; iSequence++) { 186 187 lower_[put] = COIN_DBL_MAX; 188 setInfeasible(put, true); 189 put++; 190 whichRange_[iSequence] = put; 191 lower_[put] = 0.0; 192 cost_[put++] = 0.0; 193 lower_[put] = 0.0; 194 cost_[put++] = 0.0; 195 lower_[put] = COIN_DBL_MAX; 196 setInfeasible(put  1, true); 197 cost_[put++] = 1.0e50; 198 start_[iSequence + 1] = put; 199 } 200 assert(put <= kPut); 201 } 205 202 #ifdef FAST_CLPNON 206 // See how we are storing things 207 CoinAssert (model_>clpMatrix()> 208 generalExpanded(model_, 10, iSequence) == 0); 209 #endif 210 if (CLP_METHOD2) { 211 assert (!numberExtra); 212 bound_ = new double[numberTotal]; 213 cost2_ = new double[numberTotal]; 214 status_ = new unsigned char[numberTotal]; 203 // See how we are storing things 204 CoinAssert(model_>clpMatrix()>generalExpanded(model_, 10, iSequence) == 0); 205 #endif 206 if (CLP_METHOD2) { 207 assert(!numberExtra); 208 bound_ = new double[numberTotal]; 209 cost2_ = new double[numberTotal]; 210 status_ = new unsigned char[numberTotal]; 215 211 #ifdef VALIDATE 216 delete[] saveLowerV;217 218 delete[] saveUpperV;219 220 #endif 221 222 223 224 225 226 212 delete[] saveLowerV; 213 saveLowerV = CoinCopyOfArray(model_>lowerRegion(), numberTotal); 214 delete[] saveUpperV; 215 saveUpperV = CoinCopyOfArray(model_>upperRegion(), numberTotal); 216 #endif 217 for (iSequence = 0; iSequence < numberTotal; iSequence++) { 218 bound_[iSequence] = 0.0; 219 cost2_[iSequence] = cost[iSequence]; 220 setInitialStatus(status_[iSequence]); 221 } 222 } 227 223 } 228 224 #if 0 … … 281 277 #endif 282 278 // Refresh one assuming regions OK 283 void 284 ClpNonLinearCost::refresh(int iSequence) 285 { 286 double infeasibilityCost = model_>infeasibilityCost(); 287 double primalTolerance = model_>currentPrimalTolerance(); 288 double * cost = model_>costRegion(); 289 double * upper = model_>upperRegion(); 290 double * lower = model_>lowerRegion(); 291 double * solution = model_>solutionRegion(); 292 cost2_[iSequence] = cost[iSequence]; 293 double value = solution[iSequence]; 294 double lowerValue = lower[iSequence]; 295 double upperValue = upper[iSequence]; 296 if (value  upperValue <= primalTolerance) { 297 if (value  lowerValue >= primalTolerance) { 298 // feasible 299 status_[iSequence] = static_cast<unsigned char>(CLP_FEASIBLE  (CLP_SAME << 4)); 300 bound_[iSequence] = 0.0; 301 } else { 302 // below 303 cost[iSequence] = infeasibilityCost; 304 status_[iSequence] = static_cast<unsigned char>(CLP_BELOW_LOWER  (CLP_SAME << 4)); 305 bound_[iSequence] = upperValue; 306 upper[iSequence] = lowerValue; 307 lower[iSequence] = COIN_DBL_MAX; 308 } 309 } else { 310 // above 311 cost[iSequence] += infeasibilityCost; 312 status_[iSequence] = static_cast<unsigned char>(CLP_ABOVE_UPPER  (CLP_SAME << 4)); 313 bound_[iSequence] = lowerValue; 314 lower[iSequence] = upperValue; 315 upper[iSequence] = COIN_DBL_MAX; 316 } 317 279 void ClpNonLinearCost::refresh(int iSequence) 280 { 281 double infeasibilityCost = model_>infeasibilityCost(); 282 double primalTolerance = model_>currentPrimalTolerance(); 283 double *cost = model_>costRegion(); 284 double *upper = model_>upperRegion(); 285 double *lower = model_>lowerRegion(); 286 double *solution = model_>solutionRegion(); 287 cost2_[iSequence] = cost[iSequence]; 288 double value = solution[iSequence]; 289 double lowerValue = lower[iSequence]; 290 double upperValue = upper[iSequence]; 291 if (value  upperValue <= primalTolerance) { 292 if (value  lowerValue >= primalTolerance) { 293 // feasible 294 status_[iSequence] = static_cast< unsigned char >(CLP_FEASIBLE  (CLP_SAME << 4)); 295 bound_[iSequence] = 0.0; 296 } else { 297 // below 298 cost[iSequence] = infeasibilityCost; 299 status_[iSequence] = static_cast< unsigned char >(CLP_BELOW_LOWER  (CLP_SAME << 4)); 300 bound_[iSequence] = upperValue; 301 upper[iSequence] = lowerValue; 302 lower[iSequence] = COIN_DBL_MAX; 303 } 304 } else { 305 // above 306 cost[iSequence] += infeasibilityCost; 307 status_[iSequence] = static_cast< unsigned char >(CLP_ABOVE_UPPER  (CLP_SAME << 4)); 308 bound_[iSequence] = lowerValue; 309 lower[iSequence] = upperValue; 310 upper[iSequence] = COIN_DBL_MAX; 311 } 318 312 } 319 313 // Refreshes costs always makes row costs zero 320 void 321 ClpNonLinearCost::refreshCosts(const double * columnCosts) 322 { 323 double * cost = model_>costRegion(); 324 // zero row costs 325 memset(cost + numberColumns_, 0, numberRows_ * sizeof(double)); 326 // copy column costs 327 CoinMemcpyN(columnCosts, numberColumns_, cost); 328 if ((method_ & 1) != 0) { 329 for (int iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) { 330 int start = start_[iSequence]; 331 int end = start_[iSequence+1]  1; 332 double thisFeasibleCost = cost[iSequence]; 333 if (infeasible(start)) { 334 cost_[start] = thisFeasibleCost  infeasibilityWeight_; 335 cost_[start+1] = thisFeasibleCost; 336 } else { 337 cost_[start] = thisFeasibleCost; 338 } 339 if (infeasible(end  1)) { 340 cost_[end1] = thisFeasibleCost + infeasibilityWeight_; 341 } 342 } 343 } 344 if (CLP_METHOD2) { 345 for (int iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) { 346 cost2_[iSequence] = cost[iSequence]; 347 } 348 } 349 } 350 ClpNonLinearCost::ClpNonLinearCost(ClpSimplex * model, const int * starts, 351 const double * lowerNon, const double * costNon) 314 void ClpNonLinearCost::refreshCosts(const double *columnCosts) 315 { 316 double *cost = model_>costRegion(); 317 // zero row costs 318 memset(cost + numberColumns_, 0, numberRows_ * sizeof(double)); 319 // copy column costs 320 CoinMemcpyN(columnCosts, numberColumns_, cost); 321 if ((method_ & 1) != 0) { 322 for (int iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) { 323 int start = start_[iSequence]; 324 int end = start_[iSequence + 1]  1; 325 double thisFeasibleCost = cost[iSequence]; 326 if (infeasible(start)) { 327 cost_[start] = thisFeasibleCost  infeasibilityWeight_; 328 cost_[start + 1] = thisFeasibleCost; 329 } else { 330 cost_[start] = thisFeasibleCost; 331 } 332 if (infeasible(end  1)) { 333 cost_[end  1] = thisFeasibleCost + infeasibilityWeight_; 334 } 335 } 336 } 337 if (CLP_METHOD2) { 338 for (int iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) { 339 cost2_[iSequence] = cost[iSequence]; 340 } 341 } 342 } 343 ClpNonLinearCost::ClpNonLinearCost(ClpSimplex *model, const int *starts, 344 const double *lowerNon, const double *costNon) 352 345 { 353 346 #ifndef FAST_CLPNON 354 // what about scaling?  only try without it initially 355 assert(!model>scalingFlag()); 356 model_ = model; 357 numberRows_ = model_>numberRows(); 358 numberColumns_ = model_>numberColumns(); 359 int numberTotal = numberRows_ + numberColumns_; 360 convex_ = true; 361 bothWays_ = true; 362 start_ = new int [numberTotal+1]; 363 whichRange_ = new int [numberTotal]; 364 offset_ = new int [numberTotal]; 365 memset(offset_, 0, numberTotal * sizeof(int)); 366 367 double whichWay = model_>optimizationDirection(); 368 COIN_DETAIL_PRINT(printf("Direction %g\n", whichWay)); 369 370 numberInfeasibilities_ = 0; 371 changeCost_ = 0.0; 372 feasibleCost_ = 0.0; 373 double infeasibilityCost = model_>infeasibilityCost(); 374 infeasibilityWeight_ = infeasibilityCost;; 375 largestInfeasibility_ = 0.0; 376 sumInfeasibilities_ = 0.0; 377 378 int iSequence; 379 assert (!model_>rowObjective()); 380 double * cost = model_>objective(); 381 382 // First see how much space we need 383 // Also set up feasible bounds 384 int put = starts[numberColumns_]; 385 386 double * columnUpper = model_>columnUpper(); 387 double * columnLower = model_>columnLower(); 388 for (iSequence = 0; iSequence < numberColumns_; iSequence++) { 389 if (columnLower[iSequence] > 1.0e20) 390 put++; 391 if (columnUpper[iSequence] < 1.0e20) 392 put++; 393 } 394 395 double * rowUpper = model_>rowUpper(); 396 double * rowLower = model_>rowLower(); 397 for (iSequence = 0; iSequence < numberRows_; iSequence++) { 398 if (rowLower[iSequence] > 1.0e20) 399 put++; 400 if (rowUpper[iSequence] < 1.0e20) 401 put++; 402 put += 2; 403 } 404 lower_ = new double [put]; 405 cost_ = new double [put]; 406 infeasible_ = new unsigned int[(put+31)>>5]; 407 memset(infeasible_, 0, ((put + 31) >> 5)*sizeof(unsigned int)); 408 409 // now fill in 410 put = 0; 411 412 start_[0] = 0; 413 for (iSequence = 0; iSequence < numberTotal; iSequence++) { 414 lower_[put] = COIN_DBL_MAX; 415 whichRange_[iSequence] = put + 1; 416 double thisCost; 417 double lowerValue; 418 double upperValue; 419 if (iSequence >= numberColumns_) { 420 // rows 421 lowerValue = rowLower[iSequencenumberColumns_]; 422 upperValue = rowUpper[iSequencenumberColumns_]; 423 if (lowerValue > 1.0e30) { 424 setInfeasible(put, true); 425 cost_[put++] = infeasibilityCost; 426 lower_[put] = lowerValue; 427 } 428 cost_[put++] = 0.0; 429 thisCost = 0.0; 430 } else { 431 // columns  move costs and see if convex 432 lowerValue = columnLower[iSequence]; 433 upperValue = columnUpper[iSequence]; 434 if (lowerValue > 1.0e30) { 435 setInfeasible(put, true); 436 cost_[put++] = whichWay * cost[iSequence]  infeasibilityCost; 437 lower_[put] = lowerValue; 438 } 439 int iIndex = starts[iSequence]; 440 int end = starts[iSequence+1]; 441 assert (fabs(columnLower[iSequence]  lowerNon[iIndex]) < 1.0e8); 442 thisCost = COIN_DBL_MAX; 443 for (; iIndex < end; iIndex++) { 444 if (lowerNon[iIndex] < columnUpper[iSequence]  1.0e8) { 445 lower_[put] = lowerNon[iIndex]; 446 cost_[put++] = whichWay * costNon[iIndex]; 447 // check convexity 448 if (whichWay * costNon[iIndex] < thisCost  1.0e12) 449 convex_ = false; 450 thisCost = whichWay * costNon[iIndex]; 451 } else { 452 break; 453 } 454 } 455 } 456 lower_[put] = upperValue; 457 setInfeasible(put, true); 458 cost_[put++] = thisCost + infeasibilityCost; 459 if (upperValue < 1.0e20) { 460 lower_[put] = COIN_DBL_MAX; 461 cost_[put++] = 1.0e50; 462 } 463 int iFirst = start_[iSequence]; 464 if (lower_[iFirst] != COIN_DBL_MAX) { 465 setInfeasible(iFirst, true); 466 whichRange_[iSequence] = iFirst + 1; 467 } else { 468 whichRange_[iSequence] = iFirst; 469 } 470 start_[iSequence+1] = put; 471 } 472 // can't handle nonconvex at present 473 assert(convex_); 474 status_ = NULL; 475 bound_ = NULL; 476 cost2_ = NULL; 477 method_ = 1; 347 // what about scaling?  only try without it initially 348 assert(!model>scalingFlag()); 349 model_ = model; 350 numberRows_ = model_>numberRows(); 351 numberColumns_ = model_>numberColumns(); 352 int numberTotal = numberRows_ + numberColumns_; 353 convex_ = true; 354 bothWays_ = true; 355 start_ = new int[numberTotal + 1]; 356 whichRange_ = new int[numberTotal]; 357 offset_ = new int[numberTotal]; 358 memset(offset_, 0, numberTotal * sizeof(int)); 359 360 double whichWay = model_>optimizationDirection(); 361 COIN_DETAIL_PRINT(printf("Direction %g\n", whichWay)); 362 363 numberInfeasibilities_ = 0; 364 changeCost_ = 0.0; 365 feasibleCost_ = 0.0; 366 double infeasibilityCost = model_>infeasibilityCost(); 367 infeasibilityWeight_ = infeasibilityCost; 368 ; 369 largestInfeasibility_ = 0.0; 370 sumInfeasibilities_ = 0.0; 371 372 int iSequence; 373 assert(!model_>rowObjective()); 374 double *cost = model_>objective(); 375 376 // First see how much space we need 377 // Also set up feasible bounds 378 int put = starts[numberColumns_]; 379 380 double *columnUpper = model_>columnUpper(); 381 double *columnLower = model_>columnLower(); 382 for (iSequence = 0; iSequence < numberColumns_; iSequence++) { 383 if (columnLower[iSequence] > 1.0e20) 384 put++; 385 if (columnUpper[iSequence] < 1.0e20) 386 put++; 387 } 388 389 double *rowUpper = model_>rowUpper(); 390 double *rowLower = model_>rowLower(); 391 for (iSequence = 0; iSequence < numberRows_; iSequence++) { 392 if (rowLower[iSequence] > 1.0e20) 393 put++; 394 if (rowUpper[iSequence] < 1.0e20) 395 put++; 396 put += 2; 397 } 398 lower_ = new double[put]; 399 cost_ = new double[put]; 400 infeasible_ = new unsigned int[(put + 31) >> 5]; 401 memset(infeasible_, 0, ((put + 31) >> 5) * sizeof(unsigned int)); 402 403 // now fill in 404 put = 0; 405 406 start_[0] = 0; 407 for (iSequence = 0; iSequence < numberTotal; iSequence++) { 408 lower_[put] = COIN_DBL_MAX; 409 whichRange_[iSequence] = put + 1; 410 double thisCost; 411 double lowerValue; 412 double upperValue; 413 if (iSequence >= numberColumns_) { 414 // rows 415 lowerValue = rowLower[iSequence  numberColumns_]; 416 upperValue = rowUpper[iSequence  numberColumns_]; 417 if (lowerValue > 1.0e30) { 418 setInfeasible(put, true); 419 cost_[put++] = infeasibilityCost; 420 lower_[put] = lowerValue; 421 } 422 cost_[put++] = 0.0; 423 thisCost = 0.0; 424 } else { 425 // columns  move costs and see if convex 426 lowerValue = columnLower[iSequence]; 427 upperValue = columnUpper[iSequence]; 428 if (lowerValue > 1.0e30) { 429 setInfeasible(put, true); 430 cost_[put++] = whichWay * cost[iSequence]  infeasibilityCost; 431 lower_[put] = lowerValue; 432 } 433 int iIndex = starts[iSequence]; 434 int end = starts[iSequence + 1]; 435 assert(fabs(columnLower[iSequence]  lowerNon[iIndex]) < 1.0e8); 436 thisCost = COIN_DBL_MAX; 437 for (; iIndex < end; iIndex++) { 438 if (lowerNon[iIndex] < columnUpper[iSequence]  1.0e8) { 439 lower_[put] = lowerNon[iIndex]; 440 cost_[put++] = whichWay * costNon[iIndex]; 441 // check convexity 442 if (whichWay * costNon[iIndex] < thisCost  1.0e12) 443 convex_ = false; 444 thisCost = whichWay * costNon[iIndex]; 445 } else { 446 break; 447 } 448 } 449 } 450 lower_[put] = upperValue; 451 setInfeasible(put, true); 452 cost_[put++] = thisCost + infeasibilityCost; 453 if (upperValue < 1.0e20) { 454 lower_[put] = COIN_DBL_MAX; 455 cost_[put++] = 1.0e50; 456 } 457 int iFirst = start_[iSequence]; 458 if (lower_[iFirst] != COIN_DBL_MAX) { 459 setInfeasible(iFirst, true); 460 whichRange_[iSequence] = iFirst + 1; 461 } else { 462 whichRange_[iSequence] = iFirst; 463 } 464 start_[iSequence + 1] = put; 465 } 466 // can't handle nonconvex at present 467 assert(convex_); 468 status_ = NULL; 469 bound_ = NULL; 470 cost2_ = NULL; 471 method_ = 1; 478 472 #else 479 480 473 printf("recompile ClpNonLinearCost without FAST_CLPNON\n"); 474 abort(); 481 475 #endif 482 476 } … … 484 478 // Copy constructor 485 479 // 486 ClpNonLinearCost::ClpNonLinearCost (const ClpNonLinearCost & rhs) :487 changeCost_(0.0),488 feasibleCost_(0.0),489 infeasibilityWeight_(1.0),490 largestInfeasibility_(0.0),491 sumInfeasibilities_(0.0),492 averageTheta_(0.0),493 numberRows_(rhs.numberRows_),494 numberColumns_(rhs.numberColumns_),495 start_(NULL),496 whichRange_(NULL),497 offset_(NULL),498 lower_(NULL),499 cost_(NULL),500 model_(NULL),501 infeasible_(NULL),502 numberInfeasibilities_(1),503 status_(NULL),504 bound_(NULL),505 cost2_(NULL),506 method_(rhs.method_),507 convex_(true),508 509 { 510 511 512 513 514 515 516 517 518 519 520 521 522 start_ = new int [numberTotal+1];523 524 whichRange_ = new int[numberTotal];525 526 offset_ = new int[numberTotal];527 528 529 lower_ = new double[numberEntries];530 531 cost_ = new double[numberEntries];532 533 infeasible_ = new unsigned int[(numberEntries+31)>>5];534 535 536 537 538 539 540 541 480 ClpNonLinearCost::ClpNonLinearCost(const ClpNonLinearCost &rhs) 481 : changeCost_(0.0) 482 , feasibleCost_(0.0) 483 , infeasibilityWeight_(1.0) 484 , largestInfeasibility_(0.0) 485 , sumInfeasibilities_(0.0) 486 , averageTheta_(0.0) 487 , numberRows_(rhs.numberRows_) 488 , numberColumns_(rhs.numberColumns_) 489 , start_(NULL) 490 , whichRange_(NULL) 491 , offset_(NULL) 492 , lower_(NULL) 493 , cost_(NULL) 494 , model_(NULL) 495 , infeasible_(NULL) 496 , numberInfeasibilities_(1) 497 , status_(NULL) 498 , bound_(NULL) 499 , cost2_(NULL) 500 , method_(rhs.method_) 501 , convex_(true) 502 , bothWays_(rhs.bothWays_) 503 { 504 if (numberRows_) { 505 int numberTotal = numberRows_ + numberColumns_; 506 model_ = rhs.model_; 507 numberInfeasibilities_ = rhs.numberInfeasibilities_; 508 changeCost_ = rhs.changeCost_; 509 feasibleCost_ = rhs.feasibleCost_; 510 infeasibilityWeight_ = rhs.infeasibilityWeight_; 511 largestInfeasibility_ = rhs.largestInfeasibility_; 512 sumInfeasibilities_ = rhs.sumInfeasibilities_; 513 averageTheta_ = rhs.averageTheta_; 514 convex_ = rhs.convex_; 515 if (CLP_METHOD1) { 516 start_ = new int[numberTotal + 1]; 517 CoinMemcpyN(rhs.start_, (numberTotal + 1), start_); 518 whichRange_ = new int[numberTotal]; 519 CoinMemcpyN(rhs.whichRange_, numberTotal, whichRange_); 520 offset_ = new int[numberTotal]; 521 CoinMemcpyN(rhs.offset_, numberTotal, offset_); 522 int numberEntries = start_[numberTotal]; 523 lower_ = new double[numberEntries]; 524 CoinMemcpyN(rhs.lower_, numberEntries, lower_); 525 cost_ = new double[numberEntries]; 526 CoinMemcpyN(rhs.cost_, numberEntries, cost_); 527 infeasible_ = new unsigned int[(numberEntries + 31) >> 5]; 528 CoinMemcpyN(rhs.infeasible_, ((numberEntries + 31) >> 5), infeasible_); 529 } 530 if (CLP_METHOD2) { 531 bound_ = CoinCopyOfArray(rhs.bound_, numberTotal); 532 cost2_ = CoinCopyOfArray(rhs.cost2_, numberTotal); 533 status_ = CoinCopyOfArray(rhs.status_, numberTotal); 534 } 535 } 542 536 } 543 537 … … 545 539 // Destructor 546 540 // 547 ClpNonLinearCost::~ClpNonLinearCost 548 { 549 delete[] start_;550 delete[] whichRange_;551 delete[] offset_;552 delete[] lower_;553 delete[] cost_;554 delete[] infeasible_;555 delete[] status_;556 delete[] bound_;557 delete[] cost2_;541 ClpNonLinearCost::~ClpNonLinearCost() 542 { 543 delete[] start_; 544 delete[] whichRange_; 545 delete[] offset_; 546 delete[] lower_; 547 delete[] cost_; 548 delete[] infeasible_; 549 delete[] status_; 550 delete[] bound_; 551 delete[] cost2_; 558 552 } 559 553 … … 562 556 // 563 557 ClpNonLinearCost & 564 ClpNonLinearCost::operator=(const ClpNonLinearCost &rhs)565 { 566 567 568 569 delete[] start_;570 delete[] whichRange_;571 delete[] offset_;572 delete[] lower_;573 delete[] cost_;574 delete[] infeasible_;575 delete[] status_;576 delete[] bound_;577 delete[] cost2_;578 579 580 581 582 583 584 585 586 587 588 589 590 start_ = new int [numberTotal+1];591 592 whichRange_ = new int[numberTotal];593 594 offset_ = new int[numberTotal];595 596 597 lower_ = new double[numberEntries];598 599 cost_ = new double[numberEntries];600 601 infeasible_ = new unsigned int[(numberEntries+31)>>5];602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 558 ClpNonLinearCost::operator=(const ClpNonLinearCost &rhs) 559 { 560 if (this != &rhs) { 561 numberRows_ = rhs.numberRows_; 562 numberColumns_ = rhs.numberColumns_; 563 delete[] start_; 564 delete[] whichRange_; 565 delete[] offset_; 566 delete[] lower_; 567 delete[] cost_; 568 delete[] infeasible_; 569 delete[] status_; 570 delete[] bound_; 571 delete[] cost2_; 572 start_ = NULL; 573 whichRange_ = NULL; 574 lower_ = NULL; 575 cost_ = NULL; 576 infeasible_ = NULL; 577 status_ = NULL; 578 bound_ = NULL; 579 cost2_ = NULL; 580 method_ = rhs.method_; 581 if (numberRows_) { 582 int numberTotal = numberRows_ + numberColumns_; 583 if (CLP_METHOD1) { 584 start_ = new int[numberTotal + 1]; 585 CoinMemcpyN(rhs.start_, (numberTotal + 1), start_); 586 whichRange_ = new int[numberTotal]; 587 CoinMemcpyN(rhs.whichRange_, numberTotal, whichRange_); 588 offset_ = new int[numberTotal]; 589 CoinMemcpyN(rhs.offset_, numberTotal, offset_); 590 int numberEntries = start_[numberTotal]; 591 lower_ = new double[numberEntries]; 592 CoinMemcpyN(rhs.lower_, numberEntries, lower_); 593 cost_ = new double[numberEntries]; 594 CoinMemcpyN(rhs.cost_, numberEntries, cost_); 595 infeasible_ = new unsigned int[(numberEntries + 31) >> 5]; 596 CoinMemcpyN(rhs.infeasible_, ((numberEntries + 31) >> 5), infeasible_); 597 } 598 if (CLP_METHOD2) { 599 bound_ = CoinCopyOfArray(rhs.bound_, numberTotal); 600 cost2_ = CoinCopyOfArray(rhs.cost2_, numberTotal); 601 status_ = CoinCopyOfArray(rhs.status_, numberTotal); 602 } 603 } 604 model_ = rhs.model_; 605 numberInfeasibilities_ = rhs.numberInfeasibilities_; 606 changeCost_ = rhs.changeCost_; 607 feasibleCost_ = rhs.feasibleCost_; 608 infeasibilityWeight_ = rhs.infeasibilityWeight_; 609 largestInfeasibility_ = rhs.largestInfeasibility_; 610 sumInfeasibilities_ = rhs.sumInfeasibilities_; 611 averageTheta_ = rhs.averageTheta_; 612 convex_ = rhs.convex_; 613 bothWays_ = rhs.bothWays_; 614 } 615 return *this; 622 616 } 623 617 // Changes infeasible costs and computes number and cost of infeas … … 625 619 // We will also need a 2 bit per variable array for some 626 620 // purpose which will come to me later 627 void 628 ClpNonLinearCost::checkInfeasibilities(double oldTolerance) 629 { 630 numberInfeasibilities_ = 0; 631 double infeasibilityCost = model_>infeasibilityCost(); 632 changeCost_ = 0.0; 633 largestInfeasibility_ = 0.0; 634 sumInfeasibilities_ = 0.0; 635 double primalTolerance = model_>currentPrimalTolerance(); 636 int iSequence; 637 double * solution = model_>solutionRegion(); 638 double * upper = model_>upperRegion(); 639 double * lower = model_>lowerRegion(); 640 double * cost = model_>costRegion(); 641 bool toNearest = oldTolerance <= 0.0; 642 feasibleCost_ = 0.0; 643 //bool checkCosts = (infeasibilityWeight_ != infeasibilityCost); 644 infeasibilityWeight_ = infeasibilityCost; 645 int numberTotal = numberColumns_ + numberRows_; 646 //#define NONLIN_DEBUG 621 void ClpNonLinearCost::checkInfeasibilities(double oldTolerance) 622 { 623 numberInfeasibilities_ = 0; 624 double infeasibilityCost = model_>infeasibilityCost(); 625 changeCost_ = 0.0; 626 largestInfeasibility_ = 0.0; 627 sumInfeasibilities_ = 0.0; 628 double primalTolerance = model_>currentPrimalTolerance(); 629 int iSequence; 630 double *solution = model_>solutionRegion(); 631 double *upper = model_>upperRegion(); 632 double *lower = model_>lowerRegion(); 633 double *cost = model_>costRegion(); 634 bool toNearest = oldTolerance <= 0.0; 635 feasibleCost_ = 0.0; 636 //bool checkCosts = (infeasibilityWeight_ != infeasibilityCost); 637 infeasibilityWeight_ = infeasibilityCost; 638 int numberTotal = numberColumns_ + numberRows_; 639 //#define NONLIN_DEBUG 647 640 #ifdef NONLIN_DEBUG 648 double *saveSolution = NULL;649 double *saveLower = NULL;650 double *saveUpper = NULL;651 unsigned char *saveStatus = NULL;652 653 654 655 656 657 658 641 double *saveSolution = NULL; 642 double *saveLower = NULL; 643 double *saveUpper = NULL; 644 unsigned char *saveStatus = NULL; 645 if (method_ == 3) { 646 // Save solution as we will be checking 647 saveSolution = CoinCopyOfArray(solution, numberTotal); 648 saveLower = CoinCopyOfArray(lower, numberTotal); 649 saveUpper = CoinCopyOfArray(upper, numberTotal); 650 saveStatus = CoinCopyOfArray(model_>statusArray(), numberTotal); 651 } 659 652 #else 660 assert(method_ != 3);661 #endif 662 663 664 665 666 667 668 669 670 671 int end = start_[iSequence+1]  1;672 673 674 675 676 thisFeasibleCost = cost_[start+1];677 678 679 680 thisFeasibleCost = cost_[end2];681 cost_[end1] = thisFeasibleCost + infeasibilityCost;682 683 684 if (value < lower_[iRange+1] + primalTolerance) {685 686 if (value >= lower_[iRange+1]  primalTolerance && infeasible(iRange) && iRange == start)687 688 689 690 691 692 693 694 upperValue = lower_[iRange+1];695 696 697 698 699 700 701 702 703 #if PRINT_DETAIL7 >1704 705 iSequence,solution[iSequence],706 lowerValue,upperValue);707 #endif 708 switch(status) {709 710 711 712 713 714 715 716 717 718 lowerValue = lower_[iRange+1];719 720 653 assert(method_ != 3); 654 #endif 655 if (CLP_METHOD1) { 656 // nonbasic should be at a valid bound 657 for (iSequence = 0; iSequence < numberTotal; iSequence++) { 658 double lowerValue; 659 double upperValue; 660 double value = solution[iSequence]; 661 int iRange; 662 // get correct place 663 int start = start_[iSequence]; 664 int end = start_[iSequence + 1]  1; 665 // correct costs for this infeasibility weight 666 // If free then true cost will be first 667 double thisFeasibleCost = cost_[start]; 668 if (infeasible(start)) { 669 thisFeasibleCost = cost_[start + 1]; 670 cost_[start] = thisFeasibleCost  infeasibilityCost; 671 } 672 if (infeasible(end  1)) { 673 thisFeasibleCost = cost_[end  2]; 674 cost_[end  1] = thisFeasibleCost + infeasibilityCost; 675 } 676 for (iRange = start; iRange < end; iRange++) { 677 if (value < lower_[iRange + 1] + primalTolerance) { 678 // put in better range if infeasible 679 if (value >= lower_[iRange + 1]  primalTolerance && infeasible(iRange) && iRange == start) 680 iRange++; 681 whichRange_[iSequence] = iRange; 682 break; 683 } 684 } 685 assert(iRange < end); 686 lowerValue = lower_[iRange]; 687 upperValue = lower_[iRange + 1]; 688 ClpSimplex::Status status = model_>getStatus(iSequence); 689 if (upperValue == lowerValue && status != ClpSimplex::isFixed) { 690 if (status != ClpSimplex::basic) { 691 model_>setStatus(iSequence, ClpSimplex::isFixed); 692 status = ClpSimplex::isFixed; 693 } 694 } 695 //#define PRINT_DETAIL7 2 696 #if PRINT_DETAIL7 > 1 697 printf("NL %d sol %g bounds %g %g\n", 698 iSequence, solution[iSequence], 699 lowerValue, upperValue); 700 #endif 701 switch (status) { 702 703 case ClpSimplex::basic: 704 case ClpSimplex::superBasic: 705 // iRange is in correct place 706 // slot in here 707 if (infeasible(iRange)) { 708 if (lower_[iRange] < 1.0e50) { 709 //cost_[iRange] = cost_[iRange+1]infeasibilityCost; 710 // possibly below 711 lowerValue = lower_[iRange + 1]; 712 if (value  lowerValue < primalTolerance) { 713 value = lowerValue  value  primalTolerance; 721 714 #ifndef NDEBUG 722 if(value > 1.0e15)723 724 iSequence, lowerValue, solution[iSequence], lower_[iRange+2]);715 if (value > 1.0e15) 716 printf("nonlincostb %d %g %g %g\n", 717 iSequence, lowerValue, solution[iSequence], lower_[iRange + 2]); 725 718 #endif 726 719 #if PRINT_DETAIL7 727 printf("**NL %d sol %g below %g\n", 728 iSequence,solution[iSequence], 729 lowerValue); 730 #endif 731 sumInfeasibilities_ += value; 732 largestInfeasibility_ = CoinMax(largestInfeasibility_, value); 733 changeCost_ = lowerValue * 734 (cost_[iRange]  cost[iSequence]); 735 numberInfeasibilities_++; 736 } 737 } else { 738 //cost_[iRange] = cost_[iRange1]+infeasibilityCost; 739 // possibly above 740 upperValue = lower_[iRange]; 741 if (value  upperValue > primalTolerance) { 742 value = value  upperValue  primalTolerance; 720 printf("**NL %d sol %g below %g\n", 721 iSequence, solution[iSequence], 722 lowerValue); 723 #endif 724 sumInfeasibilities_ += value; 725 largestInfeasibility_ = CoinMax(largestInfeasibility_, value); 726 changeCost_ = lowerValue * (cost_[iRange]  cost[iSequence]); 727 numberInfeasibilities_++; 728 } 729 } else { 730 //cost_[iRange] = cost_[iRange1]+infeasibilityCost; 731 // possibly above 732 upperValue = lower_[iRange]; 733 if (value  upperValue > primalTolerance) { 734 value = value  upperValue  primalTolerance; 743 735 #ifndef NDEBUG 744 if(value > 1.0e15)745 746 iSequence, lower_[iRange1], solution[iSequence], upperValue);736 if (value > 1.0e15) 737 printf("nonlincostu %d %g %g %g\n", 738 iSequence, lower_[iRange  1], solution[iSequence], upperValue); 747 739 #endif 748 740 #if PRINT_DETAIL7 749 printf("**NL %d sol %g above %g\n", 750 iSequence,solution[iSequence], 751 upperValue); 752 #endif 753 sumInfeasibilities_ += value; 754 largestInfeasibility_ = CoinMax(largestInfeasibility_, value); 755 changeCost_ = upperValue * 756 (cost_[iRange]  cost[iSequence]); 757 numberInfeasibilities_++; 758 } 759 } 760 } 761 //lower[iSequence] = lower_[iRange]; 762 //upper[iSequence] = lower_[iRange+1]; 763 //cost[iSequence] = cost_[iRange]; 764 break; 765 case ClpSimplex::isFree: 766 //if (toNearest) 767 //solution[iSequence] = 0.0; 768 break; 769 case ClpSimplex::atUpperBound: 770 if (!toNearest) { 771 // With increasing tolerances  we may be at wrong place 772 if (fabs(value  upperValue) > oldTolerance * 1.0001) { 773 if (fabs(value  lowerValue) <= oldTolerance * 1.0001) { 774 if (fabs(value  lowerValue) > primalTolerance) 775 solution[iSequence] = lowerValue; 776 model_>setStatus(iSequence, ClpSimplex::atLowerBound); 777 } else { 778 model_>setStatus(iSequence, ClpSimplex::superBasic); 779 } 780 } else if (fabs(value  upperValue) > primalTolerance) { 781 solution[iSequence] = upperValue; 782 } 783 } else { 784 // Set to nearest and make at upper bound 785 int kRange; 786 iRange = 1; 787 double nearest = COIN_DBL_MAX; 788 for (kRange = start; kRange < end; kRange++) { 789 if (fabs(lower_[kRange]  value) < nearest) { 790 nearest = fabs(lower_[kRange]  value); 791 iRange = kRange; 792 } 793 } 794 assert (iRange >= 0); 795 iRange; 796 whichRange_[iSequence] = iRange; 797 solution[iSequence] = lower_[iRange+1]; 798 } 799 break; 800 case ClpSimplex::atLowerBound: 801 if (!toNearest) { 802 // With increasing tolerances  we may be at wrong place 803 // below stops compiler error with gcc 3.2!!! 804 if (iSequence == 119) 805 printf("ZZ %g %g %g %g\n", lowerValue, value, upperValue, oldTolerance); 806 if (fabs(value  lowerValue) > oldTolerance * 1.0001) { 807 if (fabs(value  upperValue) <= oldTolerance * 1.0001) { 808 if (fabs(value  upperValue) > primalTolerance) 809 solution[iSequence] = upperValue; 810 model_>setStatus(iSequence, ClpSimplex::atUpperBound); 811 } else { 812 model_>setStatus(iSequence, ClpSimplex::superBasic); 813 } 814 } else if (fabs(value  lowerValue) > primalTolerance) { 815 solution[iSequence] = lowerValue; 816 } 817 } else { 818 // Set to nearest and make at lower bound 819 int kRange; 820 iRange = 1; 821 double nearest = COIN_DBL_MAX; 822 for (kRange = start; kRange < end; kRange++) { 823 if (fabs(lower_[kRange]  value) < nearest) { 824 nearest = fabs(lower_[kRange]  value); 825 iRange = kRange; 826 } 827 } 828 assert (iRange >= 0); 829 whichRange_[iSequence] = iRange; 830 solution[iSequence] = lower_[iRange]; 831 } 832 break; 833 case ClpSimplex::isFixed: 834 if (toNearest) { 835 // Set to true fixed 836 for (iRange = start; iRange < end; iRange++) { 837 if (lower_[iRange] == lower_[iRange+1]) 838 break; 839 } 840 if (iRange == end) { 841 // Odd  but make sensible 842 // Set to nearest and make at lower bound 843 int kRange; 844 iRange = 1; 845 double nearest = COIN_DBL_MAX; 846 for (kRange = start; kRange < end; kRange++) { 847 if (fabs(lower_[kRange]  value) < nearest) { 848 nearest = fabs(lower_[kRange]  value); 849 iRange = kRange; 850 } 851 } 852 assert (iRange >= 0); 853 whichRange_[iSequence] = iRange; 854 if (lower_[iRange] != lower_[iRange+1]) 855 model_>setStatus(iSequence, ClpSimplex::atLowerBound); 856 else 857 model_>setStatus(iSequence, ClpSimplex::atUpperBound); 858 } 859 solution[iSequence] = lower_[iRange]; 860 } 861 break; 862 } 863 lower[iSequence] = lower_[iRange]; 864 upper[iSequence] = lower_[iRange+1]; 865 cost[iSequence] = cost_[iRange]; 866 feasibleCost_ += thisFeasibleCost * solution[iSequence]; 867 //assert (iRange==whichRange_[iSequence]); 741 printf("**NL %d sol %g above %g\n", 742 iSequence, solution[iSequence], 743 upperValue); 744 #endif 745 sumInfeasibilities_ += value; 746 largestInfeasibility_ = CoinMax(largestInfeasibility_, value); 747 changeCost_ = upperValue * (cost_[iRange]  cost[iSequence]); 748 numberInfeasibilities_++; 749 } 868 750 } 869 } 751 } 752 //lower[iSequence] = lower_[iRange]; 753 //upper[iSequence] = lower_[iRange+1]; 754 //cost[iSequence] = cost_[iRange]; 755 break; 756 case ClpSimplex::isFree: 757 //if (toNearest) 758 //solution[iSequence] = 0.0; 759 break; 760 case ClpSimplex::atUpperBound: 761 if (!toNearest) { 762 // With increasing tolerances  we may be at wrong place 763 if (fabs(value  upperValue) > oldTolerance * 1.0001) { 764 if (fabs(value  lowerValue) <= oldTolerance * 1.0001) { 765 if (fabs(value  lowerValue) > primalTolerance) 766 solution[iSequence] = lowerValue; 767 model_>setStatus(iSequence, ClpSimplex::atLowerBound); 768 } else { 769 model_>setStatus(iSequence, ClpSimplex::superBasic); 770 } 771 } else if (fabs(value  upperValue) > primalTolerance) { 772 solution[iSequence] = upperValue; 773 } 774 } else { 775 // Set to nearest and make at upper bound 776 int kRange; 777 iRange = 1; 778 double nearest = COIN_DBL_MAX; 779 for (kRange = start; kRange < end; kRange++) { 780 if (fabs(lower_[kRange]  value) < nearest) { 781 nearest = fabs(lower_[kRange]  value); 782 iRange = kRange; 783 } 784 } 785 assert(iRange >= 0); 786 iRange; 787 whichRange_[iSequence] = iRange; 788 solution[iSequence] = lower_[iRange + 1]; 789 } 790 break; 791 case ClpSimplex::atLowerBound: 792 if (!toNearest) { 793 // With increasing tolerances  we may be at wrong place 794 // below stops compiler error with gcc 3.2!!! 795 if (iSequence == 119) 796 printf("ZZ %g %g %g %g\n", lowerValue, value, upperValue, oldTolerance); 797 if (fabs(value  lowerValue) > oldTolerance * 1.0001) { 798 if (fabs(value  upperValue) <= oldTolerance * 1.0001) { 799 if (fabs(value  upperValue) > primalTolerance) 800 solution[iSequence] = upperValue; 801 model_>setStatus(iSequence, ClpSimplex::atUpperBound); 802 } else { 803 model_>setStatus(iSequence, ClpSimplex::superBasic); 804 } 805 } else if (fabs(value  lowerValue) > primalTolerance) { 806 solution[iSequence] = lowerValue; 807 } 808 } else { 809 // Set to nearest and make at lower bound 810 int kRange; 811 iRange = 1; 812 double nearest = COIN_DBL_MAX; 813 for (kRange = start; kRange < end; kRange++) { 814 if (fabs(lower_[kRange]  value) < nearest) { 815 nearest = fabs(lower_[kRange]  value); 816 iRange = kRange; 817 } 818 } 819 assert(iRange >= 0); 820 whichRange_[iSequence] = iRange; 821 solution[iSequence] = lower_[iRange]; 822 } 823 break; 824 case ClpSimplex::isFixed: 825 if (toNearest) { 826 // Set to true fixed 827 for (iRange = start; iRange < end; iRange++) { 828 if (lower_[iRange] == lower_[iRange + 1]) 829 break; 830 } 831 if (iRange == end) { 832 // Odd  but make sensible 833 // Set to nearest and make at lower bound 834 int kRange; 835 iRange = 1; 836 double nearest = COIN_DBL_MAX; 837 for (kRange = start; kRange < end; kRange++) { 838 if (fabs(lower_[kRange]  value) < nearest) { 839 nearest = fabs(lower_[kRange]  value); 840 iRange = kRange; 841 } 842 } 843 assert(iRange >= 0); 844 whichRange_[iSequence] = iRange; 845 if (lower_[iRange] != lower_[iRange + 1]) 846 model_>setStatus(iSequence, ClpSimplex::atLowerBound); 847 else 848 model_>setStatus(iSequence, ClpSimplex::atUpperBound); 849 } 850 solution[iSequence] = lower_[iRange]; 851 } 852 break; 853 } 854 lower[iSequence] = lower_[iRange]; 855 upper[iSequence] = lower_[iRange + 1]; 856 cost[iSequence] = cost_[iRange]; 857 feasibleCost_ += thisFeasibleCost * solution[iSequence]; 858 //assert (iRange==whichRange_[iSequence]); 859 } 860 } 870 861 #ifdef NONLIN_DEBUG 871 double saveCost = feasibleCost_; 872 if (method_ == 3) { 873 feasibleCost_ = 0.0; 874 // Put back solution as we will be checking 875 unsigned char * statusA = model_>statusArray(); 876 for (iSequence = 0; iSequence < numberTotal; iSequence++) { 877 double value = solution[iSequence]; 878 solution[iSequence] = saveSolution[iSequence]; 879 saveSolution[iSequence] = value; 880 value = lower[iSequence]; 881 lower[iSequence] = saveLower[iSequence]; 882 saveLower[iSequence] = value; 883 value = upper[iSequence]; 884 upper[iSequence] = saveUpper[iSequence]; 885 saveUpper[iSequence] = value; 886 unsigned char value2 = statusA[iSequence]; 887 statusA[iSequence] = saveStatus[iSequence]; 888 saveStatus[iSequence] = value2; 862 double saveCost = feasibleCost_; 863 if (method_ == 3) { 864 feasibleCost_ = 0.0; 865 // Put back solution as we will be checking 866 unsigned char *statusA = model_>statusArray(); 867 for (iSequence = 0; iSequence < numberTotal; iSequence++) { 868 double value = solution[iSequence]; 869 solution[iSequence] = saveSolution[iSequence]; 870 saveSolution[iSequence] = value; 871 value = lower[iSequence]; 872 lower[iSequence] = saveLower[iSequence]; 873 saveLower[iSequence] = value; 874 value = upper[iSequence]; 875 upper[iSequence] = saveUpper[iSequence]; 876 saveUpper[iSequence] = value; 877 unsigned char value2 = statusA[iSequence]; 878 statusA[iSequence] = saveStatus[iSequence]; 879 saveStatus[iSequence] = value2; 880 } 881 } 882 #endif 883 if (CLP_METHOD2) { 884 //#define CLP_NON_JUST_BASIC 885 #ifndef CLP_NON_JUST_BASIC 886 // nonbasic should be at a valid bound 887 for (iSequence = 0; iSequence < numberTotal; iSequence++) { 888 #else 889 const int *pivotVariable = model_>pivotVariable(); 890 for (int i = 0; i < numberRows_; i++) { 891 int iSequence = pivotVariable[i]; 892 #endif 893 double value = solution[iSequence]; 894 unsigned char iStatus = status_[iSequence]; 895 assert(currentStatus(iStatus) == CLP_SAME); 896 double lowerValue = lower[iSequence]; 897 double upperValue = upper[iSequence]; 898 double costValue = cost2_[iSequence]; 899 double trueCost = costValue; 900 int iWhere = originalStatus(iStatus); 901 if (iWhere == CLP_BELOW_LOWER) { 902 lowerValue = upperValue; 903 upperValue = bound_[iSequence]; 904 costValue = infeasibilityCost; 905 } else if (iWhere == CLP_ABOVE_UPPER) { 906 upperValue = lowerValue; 907 lowerValue = bound_[iSequence]; 908 costValue += infeasibilityCost; 909 } 910 // get correct place 911 int newWhere = CLP_FEASIBLE; 912 ClpSimplex::Status status = model_>getStatus(iSequence); 913 if (upperValue == lowerValue && status != ClpSimplex::isFixed) { 914 if (status != ClpSimplex::basic) { 915 model_>setStatus(iSequence, ClpSimplex::isFixed); 916 status = ClpSimplex::isFixed; 917 } 918 } 919 switch (status) { 920 921 case ClpSimplex::basic: 922 case ClpSimplex::superBasic: 923 if (value  upperValue <= primalTolerance) { 924 if (value  lowerValue >= primalTolerance) { 925 // feasible 926 //newWhere=CLP_FEASIBLE; 927 } else { 928 // below 929 newWhere = CLP_BELOW_LOWER; 930 assert(fabs(lowerValue) < 1.0e100); 931 double infeasibility = lowerValue  value  primalTolerance; 932 sumInfeasibilities_ += infeasibility; 933 largestInfeasibility_ = CoinMax(largestInfeasibility_, infeasibility); 934 costValue = trueCost  infeasibilityCost; 935 changeCost_ = lowerValue * (costValue  cost[iSequence]); 936 numberInfeasibilities_++; 889 937 } 890 } 891 #endif 892 if (CLP_METHOD2) { 893 //#define CLP_NON_JUST_BASIC 894 #ifndef CLP_NON_JUST_BASIC 895 // nonbasic should be at a valid bound 896 for (iSequence = 0; iSequence < numberTotal; iSequence++) { 897 #else 898 const int * pivotVariable = model_>pivotVariable(); 899 for (int i=0;i<numberRows_;i++) { 900 int iSequence = pivotVariable[i]; 901 #endif 902 double value = solution[iSequence]; 903 unsigned char iStatus = status_[iSequence]; 904 assert (currentStatus(iStatus) == CLP_SAME); 905 double lowerValue = lower[iSequence]; 906 double upperValue = upper[iSequence]; 907 double costValue = cost2_[iSequence]; 908 double trueCost = costValue; 909 int iWhere = originalStatus(iStatus); 910 if (iWhere == CLP_BELOW_LOWER) { 911 lowerValue = upperValue; 912 upperValue = bound_[iSequence]; 913 costValue = infeasibilityCost; 914 } else if (iWhere == CLP_ABOVE_UPPER) { 915 upperValue = lowerValue; 916 lowerValue = bound_[iSequence]; 917 costValue += infeasibilityCost; 918 } 919 // get correct place 920 int newWhere = CLP_FEASIBLE; 921 ClpSimplex::Status status = model_>getStatus(iSequence); 922 if (upperValue == lowerValue && status != ClpSimplex::isFixed) { 923 if (status != ClpSimplex::basic) { 924 model_>setStatus(iSequence, ClpSimplex::isFixed); 925 status = ClpSimplex::isFixed; 926 } 927 } 928 switch(status) { 929 930 case ClpSimplex::basic: 931 case ClpSimplex::superBasic: 932 if (value  upperValue <= primalTolerance) { 933 if (value  lowerValue >= primalTolerance) { 934 // feasible 935 //newWhere=CLP_FEASIBLE; 936 } else { 937 // below 938 newWhere = CLP_BELOW_LOWER; 939 assert (fabs(lowerValue) < 1.0e100); 940 double infeasibility = lowerValue  value  primalTolerance; 941 sumInfeasibilities_ += infeasibility; 942 largestInfeasibility_ = CoinMax(largestInfeasibility_, infeasibility); 943 costValue = trueCost  infeasibilityCost; 944 changeCost_ = lowerValue * (costValue  cost[iSequence]); 945 numberInfeasibilities_++; 946 } 947 } else { 948 // above 949 newWhere = CLP_ABOVE_UPPER; 950 double infeasibility = value  upperValue  primalTolerance; 951 sumInfeasibilities_ += infeasibility; 952 largestInfeasibility_ = CoinMax(largestInfeasibility_, infeasibility); 953 costValue = trueCost + infeasibilityCost; 954 changeCost_ = upperValue * (costValue  cost[iSequence]); 955 numberInfeasibilities_++; 956 } 957 break; 958 case ClpSimplex::isFree: 959 break; 960 case ClpSimplex::atUpperBound: 961 if (!toNearest) { 962 // With increasing tolerances  we may be at wrong place 963 if (fabs(value  upperValue) > oldTolerance * 1.0001) { 964 if (fabs(value  lowerValue) <= oldTolerance * 1.0001) { 965 if (fabs(value  lowerValue) > primalTolerance) { 966 solution[iSequence] = lowerValue; 967 value = lowerValue; 968 } 969 model_>setStatus(iSequence, ClpSimplex::atLowerBound); 970 } else { 971 if (value < upperValue) { 972 if (value > lowerValue) { 973 model_>setStatus(iSequence, ClpSimplex::superBasic); 974 } else { 975 // set to lower bound as infeasible 976 solution[iSequence] = lowerValue; 977 value = lowerValue; 978 model_>setStatus(iSequence, ClpSimplex::atLowerBound); 979 } 980 } else { 981 // set to upper bound as infeasible 982 solution[iSequence] = upperValue; 983 value = upperValue; 984 } 985 } 986 } else if (fabs(value  upperValue) > primalTolerance) { 987 solution[iSequence] = upperValue; 988 value = upperValue; 989 } 990 } else { 991 // Set to nearest and make at bound 992 if (fabs(value  lowerValue) < fabs(value  upperValue)) { 993 solution[iSequence] = lowerValue; 994 value = lowerValue; 995 model_>setStatus(iSequence, ClpSimplex::atLowerBound); 996 } else { 997 solution[iSequence] = upperValue; 998 value = upperValue; 999 } 1000 } 1001 break; 1002 case ClpSimplex::atLowerBound: 1003 if (!toNearest) { 1004 // With increasing tolerances  we may be at wrong place 1005 if (fabs(value  lowerValue) > oldTolerance * 1.0001) { 1006 if (fabs(value  upperValue) <= oldTolerance * 1.0001) { 1007 if (fabs(value  upperValue) > primalTolerance) { 1008 solution[iSequence] = upperValue; 1009 value = upperValue; 1010 } 1011 model_>setStatus(iSequence, ClpSimplex::atUpperBound); 1012 } else { 1013 if (value < upperValue) { 1014 if (value > lowerValue) { 1015 model_>setStatus(iSequence, ClpSimplex::superBasic); 1016 } else { 1017 // set to lower bound as infeasible 1018 solution[iSequence] = lowerValue; 1019 value = lowerValue; 1020 } 1021 } else { 1022 // set to upper bound as infeasible 1023 solution[iSequence] = upperValue; 1024 value = upperValue; 1025 model_>setStatus(iSequence, ClpSimplex::atUpperBound); 1026 } 1027 } 1028 } else if (fabs(value  lowerValue) > primalTolerance) { 1029 solution[iSequence] = lowerValue; 1030 value = lowerValue; 1031 } 1032 } else { 1033 // Set to nearest and make at bound 1034 if (fabs(value  lowerValue) < fabs(value  upperValue)) { 1035 solution[iSequence] = lowerValue; 1036 value = lowerValue; 1037 } else { 1038 solution[iSequence] = upperValue; 1039 value = upperValue; 1040 model_>setStatus(iSequence, ClpSimplex::atUpperBound); 1041 } 1042 } 1043 break; 1044 case ClpSimplex::isFixed: 1045 solution[iSequence] = lowerValue; 1046 value = lowerValue; 1047 break; 1048 } 938 } else { 939 // above 940 newWhere = CLP_ABOVE_UPPER; 941 double infeasibility = value  upperValue  primalTolerance; 942 sumInfeasibilities_ += infeasibility; 943 largestInfeasibility_ = CoinMax(largestInfeasibility_, infeasibility); 944 costValue = trueCost + infeasibilityCost; 945 changeCost_ = upperValue * (costValue  cost[iSequence]); 946 numberInfeasibilities_++; 947 } 948 break; 949 case ClpSimplex::isFree: 950 break; 951 case ClpSimplex::atUpperBound: 952 if (!toNearest) { 953 // With increasing tolerances  we may be at wrong place 954 if (fabs(value  upperValue) > oldTolerance * 1.0001) { 955 if (fabs(value  lowerValue) <= oldTolerance * 1.0001) { 956 if (fabs(value  lowerValue) > primalTolerance) { 957 solution[iSequence] = lowerValue; 958 value = lowerValue; 959 } 960 model_>setStatus(iSequence, ClpSimplex::atLowerBound); 961 } else { 962 if (value < upperValue) { 963 if (value > lowerValue) { 964 model_>setStatus(iSequence, ClpSimplex::superBasic); 965 } else { 966 // set to lower bound as infeasible 967 solution[iSequence] = lowerValue; 968 value = lowerValue; 969 model_>setStatus(iSequence, ClpSimplex::atLowerBound); 970 } 971 } else { 972 // set to upper bound as infeasible 973 solution[iSequence] = upperValue; 974 value = upperValue; 975 } 976 } 977 } else if (fabs(value  upperValue) > primalTolerance) { 978 solution[iSequence] = upperValue; 979 value = upperValue; 980 } 981 } else { 982 // Set to nearest and make at bound 983 if (fabs(value  lowerValue) < fabs(value  upperValue)) { 984 solution[iSequence] = lowerValue; 985 value = lowerValue; 986 model_>setStatus(iSequence, ClpSimplex::atLowerBound); 987 } else { 988 solution[iSequence] = upperValue; 989 value = upperValue; 990 } 991 } 992 break; 993 case ClpSimplex::atLowerBound: 994 if (!toNearest) { 995 // With increasing tolerances  we may be at wrong place 996 if (fabs(value  lowerValue) > oldTolerance * 1.0001) { 997 if (fabs(value  upperValue) <= oldTolerance * 1.0001) { 998 if (fabs(value  upperValue) > primalTolerance) { 999 solution[iSequence] = upperValue; 1000 value = upperValue; 1001 } 1002 model_>setStatus(iSequence, ClpSimplex::atUpperBound); 1003 } else { 1004 if (value < upperValue) { 1005 if (value > lowerValue) { 1006 model_>setStatus(iSequence, ClpSimplex::superBasic); 1007 } else { 1008 // set to lower bound as infeasible 1009 solution[iSequence] = lowerValue; 1010 value = lowerValue; 1011 } 1012 } else { 1013 // set to upper bound as infeasible 1014 solution[iSequence] = upperValue; 1015 value = upperValue; 1016 model_>setStatus(iSequence, ClpSimplex::atUpperBound); 1017 } 1018 } 1019 } else if (fabs(value  lowerValue) > primalTolerance) { 1020 solution[iSequence] = lowerValue; 1021 value = lowerValue; 1022 } 1023 } else { 1024 // Set to nearest and make at bound 1025 if (fabs(value  lowerValue) < fabs(value  upperValue)) { 1026 solution[iSequence] = lowerValue; 1027 value = lowerValue; 1028 } else { 1029 solution[iSequence] = upperValue; 1030 value = upperValue; 1031 model_>setStatus(iSequence, ClpSimplex::atUpperBound); 1032 } 1033 } 1034 break; 1035 case ClpSimplex::isFixed: 1036 solution[iSequence] = lowerValue; 1037 value = lowerValue; 1038 break; 1039 } 1049 1040 #ifdef NONLIN_DEBUG 1050 1051 1052 1053 1054 1055 #endif 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 assert (lowerValue<=upperValue);1074 1075 1076 1041 double lo = saveLower[iSequence]; 1042 double up = saveUpper[iSequence]; 1043 double cc = cost[iSequence]; 1044 unsigned char ss = saveStatus[iSequence]; 1045 unsigned char snow = model_>statusArray()[iSequence]; 1046 #endif 1047 if (iWhere != newWhere) { 1048 setOriginalStatus(status_[iSequence], newWhere); 1049 if (newWhere == CLP_BELOW_LOWER) { 1050 bound_[iSequence] = upperValue; 1051 upperValue = lowerValue; 1052 lowerValue = COIN_DBL_MAX; 1053 costValue = trueCost  infeasibilityCost; 1054 } else if (newWhere == CLP_ABOVE_UPPER) { 1055 bound_[iSequence] = lowerValue; 1056 lowerValue = upperValue; 1057 upperValue = COIN_DBL_MAX; 1058 costValue = trueCost + infeasibilityCost; 1059 } else { 1060 costValue = trueCost; 1061 } 1062 lower[iSequence] = lowerValue; 1063 upper[iSequence] = upperValue; 1064 assert(lowerValue <= upperValue); 1065 } 1066 // always do as other things may change 1067 cost[iSequence] = costValue; 1077 1068 #ifdef NONLIN_DEBUG 1078 1079 assert(ss == snow);1080 assert(cc == cost[iSequence]);1081 assert(lo == lower[iSequence]);1082 assert(up == upper[iSequence]);1083 assert(value == saveSolution[iSequence]);1084 1085 #endif 1086 1087 1088 1069 if (method_ == 3) { 1070 assert(ss == snow); 1071 assert(cc == cost[iSequence]); 1072 assert(lo == lower[iSequence]); 1073 assert(up == upper[iSequence]); 1074 assert(value == saveSolution[iSequence]); 1075 } 1076 #endif 1077 feasibleCost_ += trueCost * value; 1078 } 1079 } 1089 1080 #ifdef NONLIN_DEBUG 1090 1091 assert(fabs(saveCost  feasibleCost_) < 0.001 * (1.0 + CoinMax(fabs(saveCost), fabs(feasibleCost_))));1092 delete[] saveSolution;1093 delete[] saveLower;1094 delete[] saveUpper;1095 delete[] saveStatus;1096 #endif 1097 1081 if (method_ == 3) 1082 assert(fabs(saveCost  feasibleCost_) < 0.001 * (1.0 + CoinMax(fabs(saveCost), fabs(feasibleCost_)))); 1083 delete[] saveSolution; 1084 delete[] saveLower; 1085 delete[] saveUpper; 1086 delete[] saveStatus; 1087 #endif 1088 //feasibleCost_ /= (model_>rhsScale()*model_>objScale()); 1098 1089 } 1099 1090 // Puts feasible bounds into lower and upper 1100 void 1101 ClpNonLinearCost::feasibleBounds() 1102 { 1103 if (CLP_METHOD2) { 1104 int iSequence; 1105 double * upper = model_>upperRegion(); 1106 double * lower = model_>lowerRegion(); 1107 double * cost = model_>costRegion(); 1108 int numberTotal = numberColumns_ + numberRows_; 1109 for (iSequence = 0; iSequence < numberTotal; iSequence++) { 1110 unsigned char iStatus = status_[iSequence]; 1111 assert (currentStatus(iStatus) == CLP_SAME); 1112 double lowerValue = lower[iSequence]; 1113 double upperValue = upper[iSequence]; 1114 double costValue = cost2_[iSequence]; 1115 int iWhere = originalStatus(iStatus); 1116 if (iWhere == CLP_BELOW_LOWER) { 1117 lowerValue = upperValue; 1118 upperValue = bound_[iSequence]; 1119 assert (fabs(lowerValue) < 1.0e100); 1120 } else if (iWhere == CLP_ABOVE_UPPER) { 1121 upperValue = lowerValue; 1122 lowerValue = bound_[iSequence]; 1123 } 1124 setOriginalStatus(status_[iSequence], CLP_FEASIBLE); 1125 lower[iSequence] = lowerValue; 1126 upper[iSequence] = upperValue; 1127 cost[iSequence] = costValue; 1128 } 1129 } 1091 void ClpNonLinearCost::feasibleBounds() 1092 { 1093 if (CLP_METHOD2) { 1094 int iSequence; 1095 double *upper = model_>upperRegion(); 1096 double *lower = model_>lowerRegion(); 1097 double *cost = model_>costRegion(); 1098 int numberTotal = numberColumns_ + numberRows_; 1099 for (iSequence = 0; iSequence < numberTotal; iSequence++) { 1100 unsigned char iStatus = status_[iSequence]; 1101 assert(currentStatus(iStatus) == CLP_SAME); 1102 double lowerValue = lower[iSequence]; 1103 double upperValue = upper[iSequence]; 1104 double costValue = cost2_[iSequence]; 1105 int iWhere = originalStatus(iStatus); 1106 if (iWhere == CLP_BELOW_LOWER) { 1107 lowerValue = upperValue; 1108 upperValue = bound_[iSequence]; 1109 assert(fabs(lowerValue) < 1.0e100); 1110 } else if (iWhere == CLP_ABOVE_UPPER) { 1111 upperValue = lowerValue; 1112 lowerValue = bound_[iSequence]; 1113 } 1114 setOriginalStatus(status_[iSequence], CLP_FEASIBLE); 1115 lower[iSequence] = lowerValue; 1116 upper[iSequence] = upperValue; 1117 cost[iSequence] = costValue; 1118 } 1119 } 1130 1120 } 1131 1121 /* Goes through one bound for each variable. … … 1133 1123 The indices are row indices and need converting to sequences 1134 1124 */ 1135 void 1136 ClpNonLinearCost::goThru(int numberInArray, double multiplier, 1137 const int * index, const double * array, 1138 double * rhs) 1139 { 1140 assert (model_ != NULL); 1141 abort(); 1142 const int * pivotVariable = model_>pivotVariable(); 1143 if (CLP_METHOD1) { 1144 for (int i = 0; i < numberInArray; i++) { 1145 int iRow = index[i]; 1146 int iSequence = pivotVariable[iRow]; 1147 double alpha = multiplier * array[iRow]; 1148 // get where in bound sequence 1149 int iRange = whichRange_[iSequence]; 1150 iRange += offset_[iSequence]; //add temporary bias 1151 double value = model_>solution(iSequence); 1152 if (alpha > 0.0) { 1153 // down one 1154 iRange; 1155 assert(iRange >= start_[iSequence]); 1156 rhs[iRow] = value  lower_[iRange]; 1157 } else { 1158 // up one 1159 iRange++; 1160 assert(iRange < start_[iSequence+1]  1); 1161 rhs[iRow] = lower_[iRange+1]  value; 1162 } 1163 offset_[iSequence] = iRange  whichRange_[iSequence]; 1164 } 1165 } 1125 void ClpNonLinearCost::goThru(int numberInArray, double multiplier, 1126 const int *index, const double *array, 1127 double *rhs) 1128 { 1129 assert(model_ != NULL); 1130 abort(); 1131 const int *pivotVariable = model_>pivotVariable(); 1132 if (CLP_METHOD1) { 1133 for (int i = 0; i < numberInArray; i++) { 1134 int iRow = index[i]; 1135 int iSequence = pivotVariable[iRow]; 1136 double alpha = multiplier * array[iRow]; 1137 // get where in bound sequence 1138 int iRange = whichRange_[iSequence]; 1139 iRange += offset_[iSequence]; //add temporary bias 1140 double value = model_>solution(iSequence); 1141 if (alpha > 0.0) { 1142 // down one 1143 iRange; 1144 assert(iRange >= start_[iSequence]); 1145 rhs[iRow] = value  lower_[iRange]; 1146 } else { 1147 // up one 1148 iRange++; 1149 assert(iRange < start_[iSequence + 1]  1); 1150 rhs[iRow] = lower_[iRange + 1]  value; 1151 } 1152 offset_[iSequence] = iRange  whichRange_[iSequence]; 1153 } 1154 } 1166 1155 #ifdef NONLIN_DEBUG 1167 double *saveRhs = NULL;1168 1169 1170 1171 1172 #endif 1173 1174 const double *solution = model_>solutionRegion();1175 const double *upper = model_>upperRegion();1176 const double *lower = model_>lowerRegion();1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 } else if(iWhere == CLP_BELOW_LOWER) {1197 assert(alpha < 0);1198 1199 1200 1201 1202 assert(iWhere == CLP_ABOVE_UPPER);1203 1204 1205 1206 1156 double *saveRhs = NULL; 1157 if (method_ == 3) { 1158 int numberRows = model_>numberRows(); 1159 saveRhs = CoinCopyOfArray(rhs, numberRows); 1160 } 1161 #endif 1162 if (CLP_METHOD2) { 1163 const double *solution = model_>solutionRegion(); 1164 const double *upper = model_>upperRegion(); 1165 const double *lower = model_>lowerRegion(); 1166 for (int i = 0; i < numberInArray; i++) { 1167 int iRow = index[i]; 1168 int iSequence = pivotVariable[iRow]; 1169 double alpha = multiplier * array[iRow]; 1170 double value = solution[iSequence]; 1171 unsigned char iStatus = status_[iSequence]; 1172 int iWhere = currentStatus(iStatus); 1173 if (iWhere == CLP_SAME) 1174 iWhere = originalStatus(iStatus); 1175 if (iWhere == CLP_FEASIBLE) { 1176 if (alpha > 0.0) { 1177 // going below 1178 iWhere = CLP_BELOW_LOWER; 1179 rhs[iRow] = value  lower[iSequence]; 1180 } else { 1181 // going above 1182 iWhere = CLP_ABOVE_UPPER; 1183 rhs[iRow] = upper[iSequence]  value; 1184 } 1185 } else if (iWhere == CLP_BELOW_LOWER) { 1186 assert(alpha < 0); 1187 // going feasible 1188 iWhere = CLP_FEASIBLE; 1189 rhs[iRow] = upper[iSequence]  value; 1190 } else { 1191 assert(iWhere == CLP_ABOVE_UPPER); 1192 // going feasible 1193 iWhere = CLP_FEASIBLE; 1194 rhs[iRow] = value  lower[iSequence]; 1195 } 1207 1196 #ifdef NONLIN_DEBUG 1208 1209 assert(rhs[iRow] == saveRhs[iRow]);1210 #endif 1211 1212 1213 1197 if (method_ == 3) 1198 assert(rhs[iRow] == saveRhs[iRow]); 1199 #endif 1200 setCurrentStatus(status_[iSequence], iWhere); 1201 } 1202 } 1214 1203 #ifdef NONLIN_DEBUG 1215 delete[] saveRhs;1204 delete[] saveRhs; 1216 1205 #endif 1217 1206 } 1218 1207 /* Takes off last iteration (i.e. offsets closer to 0) 1219 1208 */ 1220 void 1221 ClpNonLinearCost::goBack(int numberInArray, const int * index, 1222 double * rhs) 1223 { 1224 assert (model_ != NULL); 1225 abort(); 1226 const int * pivotVariable = model_>pivotVariable(); 1227 if (CLP_METHOD1) { 1228 for (int i = 0; i < numberInArray; i++) { 1229 int iRow = index[i]; 1230 int iSequence = pivotVariable[iRow]; 1231 // get where in bound sequence 1232 int iRange = whichRange_[iSequence]; 1233 // get closer to original 1234 if (offset_[iSequence] > 0) { 1235 offset_[iSequence]; 1236 assert (offset_[iSequence] >= 0); 1237 iRange += offset_[iSequence]; //add temporary bias 1238 double value = model_>solution(iSequence); 1239 // up one 1240 assert(iRange < start_[iSequence+1]  1); 1241 rhs[iRow] = lower_[iRange+1]  value; // was earlier lower_[iRange] 1242 } else { 1243 offset_[iSequence]++; 1244 assert (offset_[iSequence] <= 0); 1245 iRange += offset_[iSequence]; //add temporary bias 1246 double value = model_>solution(iSequence); 1247 // down one 1248 assert(iRange >= start_[iSequence]); 1249 rhs[iRow] = value  lower_[iRange]; // was earlier lower_[iRange+1] 1250 } 1251 } 1252 } 1209 void ClpNonLinearCost::goBack(int numberInArray, const int *index, 1210 double *rhs) 1211 { 1212 assert(model_ != NULL); 1213 abort(); 1214 const int *pivotVariable = model_>pivotVariable(); 1215 if (CLP_METHOD1) { 1216 for (int i = 0; i < numberInArray; i++) { 1217 int iRow = index[i]; 1218 int iSequence = pivotVariable[iRow]; 1219 // get where in bound sequence 1220 int iRange = whichRange_[iSequence]; 1221 // get closer to original 1222 if (offset_[iSequence] > 0) { 1223 offset_[iSequence]; 1224 assert(offset_[iSequence] >= 0); 1225 iRange += offset_[iSequence]; //add temporary bias 1226 double value = model_>solution(iSequence); 1227 // up one 1228 assert(iRange < start_[iSequence + 1]  1); 1229 rhs[iRow] = lower_[iRange + 1]  value; // was earlier lower_[iRange] 1230 } else { 1231 offset_[iSequence]++; 1232 assert(offset_[iSequence] <= 0); 1233 iRange += offset_[iSequence]; //add temporary bias 1234 double value = model_>solution(iSequence); 1235 // down one 1236 assert(iRange >= start_[iSequence]); 1237 rhs[iRow] = value  lower_[iRange]; // was earlier lower_[iRange+1] 1238 } 1239 } 1240 } 1253 1241 #ifdef NONLIN_DEBUG 1254 double *saveRhs = NULL;1255 1256 1257 1258 1259 #endif 1260 1261 const double *solution = model_>solutionRegion();1262 const double *upper = model_>upperRegion();1263 const double *lower = model_>lowerRegion();1264 1265 1266 1267 1268 1269 1270 1271 assert(iWhere != CLP_SAME);1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 } else if(iWhere == CLP_BELOW_LOWER) {1283 1284 1285 1286 1287 1288 1289 1290 1242 double *saveRhs = NULL; 1243 if (method_ == 3) { 1244 int numberRows = model_>numberRows(); 1245 saveRhs = CoinCopyOfArray(rhs, numberRows); 1246 } 1247 #endif 1248 if (CLP_METHOD2) { 1249 const double *solution = model_>solutionRegion(); 1250 const double *upper = model_>upperRegion(); 1251 const double *lower = model_>lowerRegion(); 1252 for (int i = 0; i < numberInArray; i++) { 1253 int iRow = index[i]; 1254 int iSequence = pivotVariable[iRow]; 1255 double value = solution[iSequence]; 1256 unsigned char iStatus = status_[iSequence]; 1257 int iWhere = currentStatus(iStatus); 1258 int original = originalStatus(iStatus); 1259 assert(iWhere != CLP_SAME); 1260 if (iWhere == CLP_FEASIBLE) { 1261 if (original == CLP_BELOW_LOWER) { 1262 // going below 1263 iWhere = CLP_BELOW_LOWER; 1264 rhs[iRow] = lower[iSequence]  value; 1265 } else { 1266 // going above 1267 iWhere = CLP_ABOVE_UPPER; 1268 rhs[iRow] = value  upper[iSequence]; 1269 } 1270 } else if (iWhere == CLP_BELOW_LOWER) { 1271 // going feasible 1272 iWhere = CLP_FEASIBLE; 1273 rhs[iRow] = value  upper[iSequence]; 1274 } else { 1275 // going feasible 1276 iWhere = CLP_FEASIBLE; 1277 rhs[iRow] = lower[iSequence]  value; 1278 } 1291 1279 #ifdef NONLIN_DEBUG 1292 1293 assert(rhs[iRow] == saveRhs[iRow]);1294 #endif 1295 1296 1297 1280 if (method_ == 3) 1281 assert(rhs[iRow] == saveRhs[iRow]); 1282 #endif 1283 setCurrentStatus(status_[iSequence], iWhere); 1284 } 1285 } 1298 1286 #ifdef NONLIN_DEBUG 1299 delete [] saveRhs; 1300 #endif 1301 } 1302 void 1303 ClpNonLinearCost::goBackAll(const CoinIndexedVector * update) 1304 { 1305 assert (model_ != NULL); 1306 const int * pivotVariable = model_>pivotVariable(); 1307 int number = update>getNumElements(); 1308 const int * index = update>getIndices(); 1309 if (CLP_METHOD1) { 1310 for (int i = 0; i < number; i++) { 1311 int iRow = index[i]; 1312 int iSequence = pivotVariable[iRow]; 1313 offset_[iSequence] = 0; 1314 } 1287 delete[] saveRhs; 1288 #endif 1289 } 1290 void ClpNonLinearCost::goBackAll(const CoinIndexedVector *update) 1291 { 1292 assert(model_ != NULL); 1293 const int *pivotVariable = model_>pivotVariable(); 1294 int number = update>getNumElements(); 1295 const int *index = update>getIndices(); 1296 if (CLP_METHOD1) { 1297 for (int i = 0; i < number; i++) { 1298 int iRow = index[i]; 1299 int iSequence = pivotVariable[iRow]; 1300 offset_[iSequence] = 0; 1301 } 1315 1302 #ifdef CLP_DEBUG 1316 for (i = 0; i < numberRows_ + numberColumns_; i++) 1317 assert(!offset_[i]); 1318 #endif 1319 } 1320 if (CLP_METHOD2) { 1321 for (int i = 0; i < number; i++) { 1322 int iRow = index[i]; 1323 int iSequence = pivotVariable[iRow]; 1324 setSameStatus(status_[iSequence]); 1325 } 1326 } 1327 } 1328 void 1329 ClpNonLinearCost::checkInfeasibilities(int numberInArray, const int * index) 1330 { 1331 assert (model_ != NULL); 1332 double primalTolerance = model_>currentPrimalTolerance(); 1333 const int * pivotVariable = model_>pivotVariable(); 1334 if (CLP_METHOD1) { 1335 for (int i = 0; i < numberInArray; i++) { 1336 int iRow = index[i]; 1337 int iSequence = pivotVariable[iRow]; 1338 // get where in bound sequence 1339 int iRange; 1340 int currentRange = whichRange_[iSequence]; 1341 double value = model_>solution(iSequence); 1342 int start = start_[iSequence]; 1343 int end = start_[iSequence+1]  1; 1344 for (iRange = start; iRange < end; iRange++) { 1345 if (value < lower_[iRange+1] + primalTolerance) { 1346 // put in better range 1347 if (value >= lower_[iRange+1]  primalTolerance && infeasible(iRange) && iRange == start) 1348 iRange++; 1349 break; 1350 } 1351 } 1352 assert(iRange < end); 1353 assert(model_>getStatus(iSequence) == ClpSimplex::basic); 1354 double & lower = model_>lowerAddress(iSequence); 1355 double & upper = model_>upperAddress(iSequence); 1356 double & cost = model_>costAddress(iSequence); 1357 whichRange_[iSequence] = iRange; 1358 if (iRange != currentRange) { 1359 if (infeasible(iRange)) 1360 numberInfeasibilities_++; 1361 if (infeasible(currentRange)) 1362 numberInfeasibilities_; 1363 } 1364 lower = lower_[iRange]; 1365 upper = lower_[iRange+1]; 1366 cost = cost_[iRange]; 1367 } 1368 } 1369 if (CLP_METHOD2) { 1370 double * solution = model_>solutionRegion(); 1371 double * upper = model_>upperRegion(); 1372 double * lower = model_>lowerRegion(); 1373 double * cost = model_>costRegion(); 1374 for (int i = 0; i < numberInArray; i++) { 1375 int iRow = index[i]; 1376 int iSequence = pivotVariable[iRow]; 1377 double value = solution[iSequence]; 1378 unsigned char iStatus = status_[iSequence]; 1379 assert (currentStatus(iStatus) == CLP_SAME); 1380 double lowerValue = lower[iSequence]; 1381 double upperValue = upper[iSequence]; 1382 double costValue = cost2_[iSequence]; 1383 int iWhere = originalStatus(iStatus); 1384 if (iWhere == CLP_BELOW_LOWER) { 1385 lowerValue = upperValue; 1386 upperValue = bound_[iSequence]; 1387 numberInfeasibilities_; 1388 assert (fabs(lowerValue) < 1.0e100); 1389 } else if (iWhere == CLP_ABOVE_UPPER) { 1390 upperValue = lowerValue; 1391 lowerValue = bound_[iSequence]; 1392 numberInfeasibilities_; 1393 } 1394 // get correct place 1395 int newWhere = CLP_FEASIBLE; 1396 if (value  upperValue <= primalTolerance) { 1397 if (value  lowerValue >= primalTolerance) { 1398 // feasible 1399 //newWhere=CLP_FEASIBLE; 1400 } else { 1401 // below 1402 newWhere = CLP_BELOW_LOWER; 1403 assert (fabs(lowerValue) < 1.0e100); 1404 costValue = infeasibilityWeight_; 1405 numberInfeasibilities_++; 1406 } 1407 } else { 1408 // above 1409 newWhere = CLP_ABOVE_UPPER; 1410 costValue += infeasibilityWeight_; 1411 numberInfeasibilities_++; 1412 } 1413 if (iWhere != newWhere) { 1414 setOriginalStatus(status_[iSequence], newWhere); 1415 if (newWhere == CLP_BELOW_LOWER) { 1416 bound_[iSequence] = upperValue; 1417 upperValue = lowerValue; 1418 lowerValue = COIN_DBL_MAX; 1419 } else if (newWhere == CLP_ABOVE_UPPER) { 1420 bound_[iSequence] = lowerValue; 1421 lowerValue = upperValue; 1422 upperValue = COIN_DBL_MAX; 1423 } 1424 lower[iSequence] = lowerValue; 1425 upper[iSequence] = upperValue; 1426 cost[iSequence] = costValue; 1427 } 1428 } 1429 } 1303 for (i = 0; i < numberRows_ + numberColumns_; i++) 1304 assert(!offset_[i]); 1305 #endif 1306 } 1307 if (CLP_METHOD2) { 1308 for (int i = 0; i < number; i++) { 1309 int iRow = index[i]; 1310 int iSequence = pivotVariable[iRow]; 1311 setSameStatus(status_[iSequence]); 1312 } 1313 } 1314 } 1315 void ClpNonLinearCost::checkInfeasibilities(int numberInArray, const int *index) 1316 { 1317 assert(model_ != NULL); 1318 double primalTolerance = model_>currentPrimalTolerance(); 1319 const int *pivotVariable = model_>pivotVariable(); 1320 if (CLP_METHOD1) { 1321 for (int i = 0; i < numberInArray; i++) { 1322 int iRow = index[i]; 1323 int iSequence = pivotVariable[iRow]; 1324 // get where in bound sequence 1325 int iRange; 1326 int currentRange = whichRange_[iSequence]; 1327 double value = model_>solution(iSequence); 1328 int start = start_[iSequence]; 1329 int end = start_[iSequence + 1]  1; 1330 for (iRange = start; iRange < end; iRange++) { 1331 if (value < lower_[iRange + 1] + primalTolerance) { 1332 // put in better range 1333 if (value >= lower_[iRange + 1]  primalTolerance && infeasible(iRange) && iRange == start) 1334 iRange++; 1335 break; 1336 } 1337 } 1338 assert(iRange < end); 1339 assert(model_>getStatus(iSequence) == ClpSimplex::basic); 1340 double &lower = model_>lowerAddress(iSequence); 1341 double &upper = model_>upperAddress(iSequence); 1342 double &cost = model_>costAddress(iSequence); 1343 whichRange_[iSequence] = iRange; 1344 if (iRange != currentRange) { 1345 if (infeasible(iRange)) 1346 numberInfeasibilities_++; 1347 if (infeasible(currentRange)) 1348 numberInfeasibilities_; 1349 } 1350 lower = lower_[iRange]; 1351 upper = lower_[iRange + 1]; 1352 cost = cost_[iRange]; 1353 } 1354 } 1355 if (CLP_METHOD2) { 1356 double *solution = model_>solutionRegion(); 1357 double *upper = model_>upperRegion(); 1358 double *lower = model_>lowerRegion(); 1359 double *cost = model_>costRegion(); 1360 for (int i = 0; i < numberInArray; i++) { 1361 int iRow = index[i]; 1362 int iSequence = pivotVariable[iRow]; 1363 double value = solution[iSequence]; 1364 unsigned char iStatus = status_[iSequence]; 1365 assert(currentStatus(iStatus) == CLP_SAME); 1366 double lowerValue = lower[iSequence]; 1367 double upperValue = upper[iSequence]; 1368 double costValue = cost2_[iSequence]; 1369 int iWhere = originalStatus(iStatus); 1370 if (iWhere == CLP_BELOW_LOWER) { 1371 lowerValue = upperValue; 1372 upperValue = bound_[iSequence]; 1373 numberInfeasibilities_; 1374 assert(fabs(lowerValue) < 1.0e100); 1375 } else if (iWhere == CLP_ABOVE_UPPER) { 1376 upperValue = lowerValue; 1377 lowerValue = bound_[iSequence]; 1378 numberInfeasibilities_; 1379 } 1380 // get correct place 1381 int newWhere = CLP_FEASIBLE; 1382 if (value  upperValue <= primalTolerance) { 1383 if (value  lowerValue >= primalTolerance) { 1384 // feasible 1385 //newWhere=CLP_FEASIBLE; 1386 } else { 1387 // below 1388 newWhere = CLP_BELOW_LOWER; 1389 assert(fabs(lowerValue) < 1.0e100); 1390 costValue = infeasibilityWeight_; 1391 numberInfeasibilities_++; 1392 } 1393 } else { 1394 // above 1395 newWhere = CLP_ABOVE_UPPER; 1396 costValue += infeasibilityWeight_; 1397 numberInfeasibilities_++; 1398 } 1399 if (iWhere != newWhere) { 1400 setOriginalStatus(status_[iSequence], newWhere); 1401 if (newWhere == CLP_BELOW_LOWER) { 1402 bound_[iSequence] = upperValue; 1403 upperValue = lowerValue; 1404 lowerValue = COIN_DBL_MAX; 1405 } else if (newWhere == CLP_ABOVE_UPPER) { 1406 bound_[iSequence] = lowerValue; 1407 lowerValue = upperValue; 1408 upperValue = COIN_DBL_MAX; 1409 } 1410 lower[iSequence] = lowerValue; 1411 upper[iSequence] = upperValue; 1412 cost[iSequence] = costValue; 1413 } 1414 } 1415 } 1430 1416 } 1431 1417 /* Puts back correct infeasible costs for each variable … … 1435 1421 changed costs will be stored as normal CoinIndexedVector 1436 1422 */ 1437 void 1438 ClpNonLinearCost::checkChanged(int numberInArray, CoinIndexedVector * update) 1439 { 1440 assert (model_ != NULL); 1441 double primalTolerance = model_>currentPrimalTolerance(); 1442 const int * pivotVariable = model_>pivotVariable(); 1443 int number = 0; 1444 int * index = update>getIndices(); 1445 double * work = update>denseVector(); 1446 if (CLP_METHOD1) { 1447 for (int i = 0; i < numberInArray; i++) { 1448 int iRow = index[i]; 1449 int iSequence = pivotVariable[iRow]; 1450 // get where in bound sequence 1451 int iRange; 1452 double value = model_>solution(iSequence); 1453 int start = start_[iSequence]; 1454 int end = start_[iSequence+1]  1; 1455 for (iRange = start; iRange < end; iRange++) { 1456 if (value < lower_[iRange+1] + primalTolerance) { 1457 // put in better range 1458 if (value >= lower_[iRange+1]  primalTolerance && infeasible(iRange) && iRange == start) 1459 iRange++; 1460 break; 1461 } 1462 } 1463 assert(iRange < end); 1464 assert(model_>getStatus(iSequence) == ClpSimplex::basic); 1465 int jRange = whichRange_[iSequence]; 1466 if (iRange != jRange) { 1467 // changed 1468 work[iRow] = cost_[jRange]  cost_[iRange]; 1469 index[number++] = iRow; 1470 double & lower = model_>lowerAddress(iSequence); 1471 double & upper = model_>upperAddress(iSequence); 1472 double & cost = model_>costAddress(iSequence); 1473 whichRange_[iSequence] = iRange; 1474 if (infeasible(iRange)) 1475 numberInfeasibilities_++; 1476 if (infeasible(jRange)) 1477 numberInfeasibilities_; 1478 lower = lower_[iRange]; 1479 upper = lower_[iRange+1]; 1480 cost = cost_[iRange]; 1481 } 1482 } 1483 } 1484 if (CLP_METHOD2) { 1485 double * solution = model_>solutionRegion(); 1486 double * upper = model_>upperRegion(); 1487 double * lower = model_>lowerRegion(); 1488 double * cost = model_>costRegion(); 1489 for (int i = 0; i < numberInArray; i++) { 1490 int iRow = index[i]; 1491 int iSequence = pivotVariable[iRow]; 1492 double value = solution[iSequence]; 1493 unsigned char iStatus = status_[iSequence]; 1494 assert (currentStatus(iStatus) == CLP_SAME); 1495 double lowerValue = lower[iSequence]; 1496 double upperValue = upper[iSequence]; 1497 double costValue = cost2_[iSequence]; 1498 int iWhere = originalStatus(iStatus); 1499 if (iWhere == CLP_BELOW_LOWER) { 1500 lowerValue = upperValue; 1501 upperValue = bound_[iSequence]; 1502 numberInfeasibilities_; 1503 assert (fabs(lowerValue) < 1.0e100); 1504 } else if (iWhere == CLP_ABOVE_UPPER) { 1505 upperValue = lowerValue; 1506 lowerValue = bound_[iSequence]; 1507 numberInfeasibilities_; 1508 } 1509 // get correct place 1510 int newWhere = CLP_FEASIBLE; 1511 if (value  upperValue <= primalTolerance) { 1512 if (value  lowerValue >= primalTolerance) { 1513 // feasible 1514 //newWhere=CLP_FEASIBLE; 1515 } else { 1516 // below 1517 newWhere = CLP_BELOW_LOWER; 1518 costValue = infeasibilityWeight_; 1519 numberInfeasibilities_++; 1520 assert (fabs(lowerValue) < 1.0e100); 1521 } 1522 } else { 1523 // above 1524 newWhere = CLP_ABOVE_UPPER; 1525 costValue += infeasibilityWeight_; 1526 numberInfeasibilities_++; 1527 } 1528 if (iWhere != newWhere) { 1529 work[iRow] = cost[iSequence]  costValue; 1530 index[number++] = iRow; 1531 setOriginalStatus(status_[iSequence], newWhere); 1532 if (newWhere == CLP_BELOW_LOWER) { 1533 bound_[iSequence] = upperValue; 1534 upperValue = lowerValue; 1535 lowerValue = COIN_DBL_MAX; 1536 } else if (newWhere == CLP_ABOVE_UPPER) { 1537 bound_[iSequence] = lowerValue; 1538 lowerValue = upperValue; 1539 upperValue = COIN_DBL_MAX; 1540 } 1541 lower[iSequence] = lowerValue; 1542 upper[iSequence] = upperValue; 1543 cost[iSequence] = costValue; 1544 } 1545 } 1546 } 1547 update>setNumElements(number); 1423 void ClpNonLinearCost::checkChanged(int numberInArray, CoinIndexedVector *update) 1424 { 1425 assert(model_ != NULL); 1426 double primalTolerance = model_>currentPrimalTolerance(); 1427 const int *pivotVariable = model_>pivotVariable(); 1428 int number = 0; 1429 int *index = update>getIndices(); 1430 double *work = update>denseVector(); 1431 if (CLP_METHOD1) { 1432 for (int i = 0; i < numberInArray; i++) { 1433 int iRow = index[i]; 1434 int iSequence = pivotVariable[iRow]; 1435 // get where in bound sequence 1436 int iRange; 1437 double value = model_>solution(iSequence); 1438 int start = start_[iSequence]; 1439 int end = start_[iSequence + 1]  1; 1440 for (iRange = start; iRange < end; iRange++) { 1441 if (value < lower_[iRange + 1] + primalTolerance) { 1442 // put in better range 1443 if (value >= lower_[iRange + 1]  primalTolerance && infeasible(iRange) && iRange == start) 1444 iRange++; 1445 break; 1446 } 1447 } 1448 assert(iRange < end); 1449 assert(model_>getStatus(iSequence) == ClpSimplex::basic); 1450 int jRange = whichRange_[iSequence]; 1451 if (iRange != jRange) { 1452 // changed 1453 work[iRow] = cost_[jRange]  cost_[iRange]; 1454 index[number++] = iRow; 1455 double &lower = model_>lowerAddress(iSequence); 1456 double &upper = model_>upperAddress(iSequence); 1457 double &cost = model_>costAddress(iSequence); 1458 whichRange_[iSequence] = iRange; 1459 if (infeasible(iRange)) 1460 numberInfeasibilities_++; 1461 if (infeasible(jRange)) 1462 numberInfeasibilities_; 1463 lower = lower_[iRange]; 1464 upper = lower_[iRange + 1]; 1465 cost = cost_[iRange]; 1466 } 1467 } 1468 } 1469 if (CLP_METHOD2) { 1470 double *solution = model_>solutionRegion(); 1471 double *upper = model_>upperRegion(); 1472 double *lower = model_>lowerRegion(); 1473 double *cost = model_>costRegion(); 1474 for (int i = 0; i < numberInArray; i++) { 1475 int iRow = index[i]; 1476 int iSequence = pivotVariable[iRow]; 1477 double value = solution[iSequence]; 1478 unsigned char iStatus = status_[iSequence]; 1479 assert(currentStatus(iStatus) == CLP_SAME); 1480 double lowerValue = lower[iSequence]; 1481 double upperValue = upper[iSequence]; 1482 double costValue = cost2_[iSequence]; 1483 int iWhere = originalStatus(iStatus); 1484 if (iWhere == CLP_BELOW_LOWER) { 1485 lowerValue = upperValue; 1486 upperValue = bound_[iSequence]; 1487 numberInfeasibilities_; 1488 assert(fabs(lowerValue) < 1.0e100); 1489 } else if (iWhere == CLP_ABOVE_UPPER) { 1490 upperValue = lowerValue; 1491 lowerValue = bound_[iSequence]; 1492 numberInfeasibilities_; 1493 } 1494 // get correct place 1495 int newWhere = CLP_FEASIBLE; 1496 if (value  upperValue <= primalTolerance) { 1497 if (value  lowerValue >= primalTolerance) { 1498 // feasible 1499 //newWhere=CLP_FEASIBLE; 1500 } else { 1501 // below 1502 newWhere = CLP_BELOW_LOWER; 1503 costValue = infeasibilityWeight_; 1504 numberInfeasibilities_++; 1505 assert(fabs(lowerValue) < 1.0e100); 1506 } 1507 } else { 1508 // above 1509 newWhere = CLP_ABOVE_UPPER; 1510 costValue += infeasibilityWeight_; 1511 numberInfeasibilities_++; 1512 } 1513 if (iWhere != newWhere) { 1514 work[iRow] = cost[iSequence]  costValue; 1515 index[number++] = iRow; 1516 setOriginalStatus(status_[iSequence], newWhere); 1517 if (newWhere == CLP_BELOW_LOWER) { 1518 bound_[iSequence] = upperValue; 1519 upperValue = lowerValue; 1520 lowerValue = COIN_DBL_MAX; 1521 } else if (newWhere == CLP_ABOVE_UPPER) { 1522 bound_[iSequence] = lowerValue; 1523 lowerValue = upperValue; 1524 upperValue = COIN_DBL_MAX; 1525 } 1526 lower[iSequence] = lowerValue; 1527 upper[iSequence] = upperValue; 1528 cost[iSequence] = costValue; 1529 } 1530 } 1531 } 1532 update>setNumElements(number); 1548 1533 } 1549 1534 /* Sets bounds and cost for one variable  returns change in cost*/ … … 1551 1536 ClpNonLinearCost::setOne(int iSequence, double value) 1552 1537 { 1553 assert (model_ != NULL); 1554 double primalTolerance = model_>currentPrimalTolerance(); 1555 // difference in cost 1556 double difference = 0.0; 1557 if (CLP_METHOD1) { 1558 // get where in bound sequence 1559 int iRange; 1560 int currentRange = whichRange_[iSequence]; 1561 int start = start_[iSequence]; 1562 int end = start_[iSequence+1]  1; 1563 if (!bothWays_) { 1564 // If fixed try and get feasible 1565 if (lower_[start+1] == lower_[start+2] && fabs(value  lower_[start+1]) < 1.001 * primalTolerance) { 1566 iRange = start + 1; 1567 } else { 1568 for (iRange = start; iRange < end; iRange++) { 1569 if (value <= lower_[iRange+1] + primalTolerance) { 1570 // put in better range 1571 if (value >= lower_[iRange+1]  primalTolerance && infeasible(iRange) && iRange == start) 1572 iRange++; 1573 break; 1574 } 1575 } 1576 } 1577 } else { 1578 // leave in current if possible 1579 iRange = whichRange_[iSequence]; 1580 if (value < lower_[iRange]  primalTolerance  value > lower_[iRange+1] + primalTolerance) { 1581 for (iRange = start; iRange < end; iRange++) { 1582 if (value < lower_[iRange+1] + primalTolerance) { 1583 // put in better range 1584 if (value >= lower_[iRange+1]  primalTolerance && infeasible(iRange) && iRange == start) 1585 iRange++; 1586 break; 1587 } 1588 } 1589 } 1538 assert(model_ != NULL); 1539 double primalTolerance = model_>currentPrimalTolerance(); 1540 // difference in cost 1541 double difference = 0.0; 1542 if (CLP_METHOD1) { 1543 // get where in bound sequence 1544 int iRange; 1545 int currentRange = whichRange_[iSequence]; 1546 int start = start_[iSequence]; 1547 int end = start_[iSequence + 1]  1; 1548 if (!bothWays_) { 1549 // If fixed try and get feasible 1550 if (lower_[start + 1] == lower_[start + 2] && fabs(value  lower_[start + 1]) < 1.001 * primalTolerance) { 1551 iRange = start + 1; 1552 } else { 1553 for (iRange = start; iRange < end; iRange++) { 1554 if (value <= lower_[iRange + 1] + primalTolerance) { 1555 // put in better range 1556 if (value >= lower_[iRange + 1]  primalTolerance && infeasible(iRange) && iRange == start) 1557 iRange++; 1558 break; 1590 1559 } 1591 assert(iRange < end); 1592 whichRange_[iSequence] = iRange; 1593 if (iRange != currentRange) { 1594 if (infeasible(iRange)) 1595 numberInfeasibilities_++; 1596 if (infeasible(currentRange)) 1597 numberInfeasibilities_; 1560 } 1561 } 1562 } else { 1563 // leave in current if possible 1564 iRange = whichRange_[iSequence]; 1565 if (value < lower_[iRange]  primalTolerance  value > lower_[iRange + 1] + primalTolerance) { 1566 for (iRange = start; iRange < end; iRange++) { 1567 if (value < lower_[iRange + 1] + primalTolerance) { 1568 // put in better range 1569 if (value >= lower_[iRange + 1]  primalTolerance && infeasible(iRange) && iRange == start) 1570 iRange++; 1571 break; 1598 1572 } 1599 double & lower = model_>lowerAddress(iSequence); 1600 double & upper = model_>upperAddress(iSequence); 1601 double & cost = model_>costAddress(iSequence); 1602 lower = lower_[iRange]; 1603 upper = lower_[iRange+1]; 1604 ClpSimplex::Status status = model_>getStatus(iSequence); 1605 if (upper == lower) { 1606 if (status != ClpSimplex::basic) { 1607 model_>setStatus(iSequence, ClpSimplex::isFixed); 1608 status = ClpSimplex::basic; // so will skip 1609 } 1610 } 1611 switch(status) { 1612 1613 case ClpSimplex::basic: 1614 case ClpSimplex::superBasic: 1615 case ClpSimplex::isFree: 1616 break; 1617 case ClpSimplex::atUpperBound: 1618 case ClpSimplex::atLowerBound: 1619 case ClpSimplex::isFixed: 1620 // set correctly 1621 if (fabs(value  lower) <= primalTolerance * 1.001) { 1622 model_>setStatus(iSequence, ClpSimplex::atLowerBound); 1623 } else if (fabs(value  upper) <= primalTolerance * 1.001) { 1624 model_>setStatus(iSequence, ClpSimplex::atUpperBound); 1625 } else { 1626 // set superBasic 1627 model_>setStatus(iSequence, ClpSimplex::superBasic); 1628 } 1629 break; 1630 } 1631 difference = cost  cost_[iRange]; 1632 cost = cost_[iRange]; 1633 } 1634 if (CLP_METHOD2) { 1635 double * upper = model_>upperRegion(); 1636 double * lower = model_>lowerRegion(); 1637 double * cost = model_>costRegion(); 1638 unsigned char iStatus = status_[iSequence]; 1639 assert (currentStatus(iStatus) == CLP_SAME); 1640 double lowerValue = lower[iSequence]; 1641 double upperValue = upper[iSequence]; 1642 double costValue = cost2_[iSequence]; 1643 int iWhere = originalStatus(iStatus); 1573 } 1574 } 1575 } 1576 assert(iRange < end); 1577 whichRange_[iSequence] = iRange; 1578 if (iRange != currentRange) { 1579 if (infeasible(iRange)) 1580 numberInfeasibilities_++; 1581 if (infeasible(currentRange)) 1582 numberInfeasibilities_; 1583 } 1584 double &lower = model_>lowerAddress(iSequence); 1585 double &upper = model_>upperAddress(iSequence); 1586 double &cost = model_>costAddress(iSequence); 1587 lower = lower_[iRange]; 1588 upper = lower_[iRange + 1]; 1589 ClpSimplex::Status status = model_>getStatus(iSequence); 1590 if (upper == lower) { 1591 if (status != ClpSimplex::basic) { 1592 model_>setStatus(iSequence, ClpSimplex::isFixed); 1593 status = ClpSimplex::basic; // so will skip 1594 } 1595 } 1596 switch (status) { 1597 1598 case ClpSimplex::basic: 1599 case ClpSimplex::superBasic: 1600 case ClpSimplex::isFree: 1601 break; 1602 case ClpSimplex::atUpperBound: 1603 case ClpSimplex::atLowerBound: 1604 case ClpSimplex::isFixed: 1605 // set correctly 1606 if (fabs(value  lower) <= primalTolerance * 1.001) { 1607 model_>setStatus(iSequence, ClpSimplex::atLowerBound); 1608 } else if (fabs(value  upper) <= primalTolerance * 1.001) { 1609 model_>setStatus(iSequence, ClpSimplex::atUpperBound); 1610 } else { 1611 // set superBasic 1612 model_>setStatus(iSequence, ClpSimplex::superBasic); 1613 } 1614 break; 1615 } 1616 difference = cost  cost_[iRange]; 1617 cost = cost_[iRange]; 1618 } 1619 if (CLP_METHOD2) { 1620 double *upper = model_>upperRegion(); 1621 double *lower = model_>lowerRegion(); 1622 double *cost = model_>costRegion(); 1623 unsigned char iStatus = status_[iSequence]; 1624 assert(currentStatus(iStatus) == CLP_SAME); 1625 double lowerValue = lower[iSequence]; 1626 double upperValue = upper[iSequence]; 1627 double costValue = cost2_[iSequence]; 1628 int iWhere = originalStatus(iStatus); 1644 1629 #undef CLP_USER_DRIVEN 1645 1630 #ifdef CLP_USER_DRIVEN 1646 1647 info.type=3;1648 info.sequence=iSequence;1649 info.value=value;1650 info.change=0;1651 info.lower=lowerValue;1652 info.upper=upperValue;1653 info.cost=costValue;1654 #endif 1655 1656 1657 1658 1659 assert(fabs(lowerValue) < 1.0e100);1660 1661 1662 1663 1664 1665 1666 1667 1668 1669 1670 1671 1672 1673 1674 1675 1676 assert(fabs(lowerValue) < 1.0e100);1677 1678 1679 1680 1681 1682 1683 1684 1685 1686 1687 1688 1689 1690 1691 1692 1693 1694 1695 1696 1697 1698 1699 1700 1701 1702 1703 1704 1705 1706 1707 switch(status) {1708 1709 1710 1711 1712 1713 1714 1715 1716 1717 1718 1719 1720 1721 1722 1723 1724 1725 1726 1631 clpUserStruct info; 1632 info.type = 3; 1633 info.sequence = iSequence; 1634 info.value = value; 1635 info.change = 0; 1636 info.lower = lowerValue; 1637 info.upper = upperValue; 1638 info.cost = costValue; 1639 #endif 1640 if (iWhere == CLP_BELOW_LOWER) { 1641 lowerValue = upperValue; 1642 upperValue = bound_[iSequence]; 1643 numberInfeasibilities_; 1644 assert(fabs(lowerValue) < 1.0e100); 1645 } else if (iWhere == CLP_ABOVE_UPPER) { 1646 upperValue = lowerValue; 1647 lowerValue = bound_[iSequence]; 1648 numberInfeasibilities_; 1649 } 1650 // get correct place 1651 int newWhere = CLP_FEASIBLE; 1652 if (value  upperValue <= primalTolerance) { 1653 if (value  lowerValue >= primalTolerance) { 1654 // feasible 1655 //newWhere=CLP_FEASIBLE; 1656 } else { 1657 // below 1658 newWhere = CLP_BELOW_LOWER; 1659 costValue = infeasibilityWeight_; 1660 numberInfeasibilities_++; 1661 assert(fabs(lowerValue) < 1.0e100); 1662 } 1663 } else { 1664 // above 1665 newWhere = CLP_ABOVE_UPPER; 1666 costValue += infeasibilityWeight_; 1667 numberInfeasibilities_++; 1668 } 1669 if (iWhere != newWhere) { 1670 difference = cost[iSequence]  costValue; 1671 setOriginalStatus(status_[iSequence], newWhere); 1672 if (newWhere == CLP_BELOW_LOWER) { 1673 bound_[iSequence] = upperValue; 1674 upperValue = lowerValue; 1675 lowerValue = COIN_DBL_MAX; 1676 } else if (newWhere == CLP_ABOVE_UPPER) { 1677 bound_[iSequence] = lowerValue; 1678 lowerValue = upperValue; 1679 upperValue = COIN_DBL_MAX; 1680 } 1681 lower[iSequence] = lowerValue; 1682 upper[iSequence] = upperValue; 1683 cost[iSequence] = costValue; 1684 } 1685 ClpSimplex::Status status = model_>getStatus(iSequence); 1686 if (upperValue == lowerValue) { 1687 if (status != ClpSimplex::basic) { 1688 model_>setStatus(iSequence, ClpSimplex::isFixed); 1689 status = ClpSimplex::basic; // so will skip 1690 } 1691 } 1692 switch (status) { 1693 1694 case ClpSimplex::basic: 1695 case ClpSimplex::superBasic: 1696 case ClpSimplex::isFree: 1697 break; 1698 case ClpSimplex::atUpperBound: 1699 case ClpSimplex::atLowerBound: 1700 case ClpSimplex::isFixed: 1701 // set correctly 1702 if (fabs(value  lowerValue) <= primalTolerance * 1.001) { 1703 model_>setStatus(iSequence, ClpSimplex::atLowerBound); 1704 } else if (fabs(value  upperValue) <= primalTolerance * 1.001) { 1705 model_>setStatus(iSequence, ClpSimplex::atUpperBound); 1706 } else { 1707 // set superBasic 1708 model_>setStatus(iSequence, ClpSimplex::superBasic); 1709 } 1710 break; 1711 } 1727 1712 #ifdef CLP_USER_DRIVEN 1728 1729 1730 1731 1732 #endif 1733 1734 1735 1713 model_>eventHandler()>eventWithInfo(ClpEventHandler::pivotRow, 1714 &info); 1715 if (info.change) { 1716 } 1717 #endif 1718 } 1719 changeCost_ += value * difference; 1720 return difference; 1736 1721 } 1737 1722 /* Sets bounds and infeasible cost and true cost for one variable 1738 1723 This is for gub and column generation etc */ 1739 void 1740 ClpNonLinearCost::setOne(int sequence, double solutionValue, double lowerValue, double upperValue, 1741 double costValue) 1742 { 1743 if (CLP_METHOD1) { 1744 int iRange = 1; 1745 int start = start_[sequence]; 1746 double infeasibilityCost = model_>infeasibilityCost(); 1747 cost_[start] = costValue  infeasibilityCost; 1748 lower_[start+1] = lowerValue; 1749 cost_[start+1] = costValue; 1750 lower_[start+2] = upperValue; 1751 cost_[start+2] = costValue + infeasibilityCost; 1752 double primalTolerance = model_>currentPrimalTolerance(); 1753 if (solutionValue  lowerValue >= primalTolerance) { 1754 if (solutionValue  upperValue <= primalTolerance) { 1755 iRange = start + 1; 1756 } else { 1757 iRange = start + 2; 1758 } 1759 } else { 1760 iRange = start; 1761 } 1762 model_>costRegion()[sequence] = cost_[iRange]; 1763 whichRange_[sequence] = iRange; 1764 } 1765 if (CLP_METHOD2) { 1766 bound_[sequence]=0.0; 1767 cost2_[sequence]=costValue; 1768 setInitialStatus(status_[sequence]); 1769 } 1770 1724 void ClpNonLinearCost::setOne(int sequence, double solutionValue, double lowerValue, double upperValue, 1725 double costValue) 1726 { 1727 if (CLP_METHOD1) { 1728 int iRange = 1; 1729 int start = start_[sequence]; 1730 double infeasibilityCost = model_>infeasibilityCost(); 1731 cost_[start] = costValue  infeasibilityCost; 1732 lower_[start + 1] = lowerValue; 1733 cost_[start + 1] = costValue; 1734 lower_[start + 2] = upperValue; 1735 cost_[start + 2] = costValue + infeasibilityCost; 1736 double primalTolerance = model_>currentPrimalTolerance(); 1737 if (solutionValue  lowerValue >= primalTolerance) { 1738 if (solutionValue  upperValue <= primalTolerance) { 1739 iRange = start + 1; 1740 } else { 1741 iRange = start + 2; 1742 } 1743 } else { 1744 iRange = start; 1745 } 1746 model_>costRegion()[sequence] = cost_[iRange]; 1747 whichRange_[sequence] = iRange; 1748 } 1749 if (CLP_METHOD2) { 1750 bound_[sequence] = 0.0; 1751 cost2_[sequence] = costValue; 1752 setInitialStatus(status_[sequence]); 1753 } 1771 1754 } 1772 1755 /* Sets bounds and cost for outgoing variable 1773 1756 may change value 1774 1757 Returns direction */ 1775 int 1776 ClpNonLinearCost::setOneOutgoing(int iSequence, double & value) 1777 { 1778 assert (model_ != NULL); 1779 double primalTolerance = model_>currentPrimalTolerance(); 1780 // difference in cost 1781 double difference = 0.0; 1782 int direction = 0; 1758 int ClpNonLinearCost::setOneOutgoing(int iSequence, double &value) 1759 { 1760 assert(model_ != NULL); 1761 double primalTolerance = model_>currentPrimalTolerance(); 1762 // difference in cost 1763 double difference = 0.0; 1764 int direction = 0; 1783 1765 #ifdef CLP_USER_DRIVEN 1784 double saveLower = model_>lowerRegion()[iSequence]; 1785 double saveUpper = model_>upperRegion()[iSequence]; 1786 double saveCost = model_>costRegion()[iSequence]; 1787 #endif 1788 if (CLP_METHOD1) { 1789 // get where in bound sequence 1790 int iRange; 1791 int currentRange = whichRange_[iSequence]; 1792 int start = start_[iSequence]; 1793 int end = start_[iSequence+1]  1; 1794 // Set perceived direction out 1795 if (value <= lower_[currentRange] + 1.001 * primalTolerance) { 1796 direction = 1; 1797 } else if (value >= lower_[currentRange+1]  1.001 * primalTolerance) { 1798 direction = 1; 1799 } else { 1800 // odd 1801 direction = 0; 1766 double saveLower = model_>lowerRegion()[iSequence]; 1767 double saveUpper = model_>upperRegion()[iSequence]; 1768 double saveCost = model_>costRegion()[iSequence]; 1769 #endif 1770 if (CLP_METHOD1) { 1771 // get where in bound sequence 1772 int iRange; 1773 int currentRange = whichRange_[iSequence]; 1774 int start = start_[iSequence]; 1775 int end = start_[iSequence + 1]  1; 1776 // Set perceived direction out 1777 if (value <= lower_[currentRange] + 1.001 * primalTolerance) { 1778 direction = 1; 1779 } else if (value >= lower_[currentRange + 1]  1.001 * primalTolerance) { 1780 direction = 1; 1781 } else { 1782 // odd 1783 direction = 0; 1784 } 1785 // If fixed try and get feasible 1786 if (lower_[start + 1] == lower_[start + 2] && fabs(value  lower_[start + 1]) < 1.001 * primalTolerance) { 1787 iRange = start + 1; 1788 } else { 1789 // See if exact 1790 for (iRange = start; iRange < end; iRange++) { 1791 if (value == lower_[iRange + 1]) { 1792 // put in better range 1793 if (infeasible(iRange) && iRange == start) 1794 iRange++; 1795 break; 1796 } 1797 } 1798 if (iRange == end) { 1799 // not exact 1800 for (iRange = start; iRange < end; iRange++) { 1801 if (value <= lower_[iRange + 1] + primalTolerance) { 1802 // put in better range 1803 if (value >= lower_[iRange + 1]  primalTolerance && infeasible(iRange) && iRange == start) 1804 iRange++; 1805 break; 1802 1806 } 1803 // If fixed try and get feasible 1804 if (lower_[start+1] == lower_[start+2] && fabs(value  lower_[start+1]) < 1.001 * primalTolerance) { 1805 iRange = start + 1; 1806 } else { 1807 // See if exact 1808 for (iRange = start; iRange < end; iRange++) { 1809 if (value == lower_[iRange+1]) { 1810 // put in better range 1811 if (infeasible(iRange) && iRange == start) 1812 iRange++; 1813 break; 1814 } 1815 } 1816 if (iRange == end) { 1817 // not exact 1818 for (iRange = start; iRange < end; iRange++) { 1819 if (value <= lower_[iRange+1] + primalTolerance) { 1820 // put in better range 1821 if (value >= lower_[iRange+1]  primalTolerance && infeasible(iRange) && iRange == start) 1822 iRange++; 1823 break; 1824 } 1825 } 1826 } 1827 } 1828 assert(iRange < end); 1829 whichRange_[iSequence] = iRange; 1830 if (iRange != currentRange) { 1831 if (infeasible(iRange)) 1832 numberInfeasibilities_++; 1833 if (infeasible(currentRange)) 1834 numberInfeasibilities_; 1835 } 1836 double & lower = model_>lowerAddress(iSequence); 1837 double & upper = model_>upperAddress(iSequence); 1838 double & cost = model_>costAddress(iSequence); 1839 lower = lower_[iRange]; 1840 upper = lower_[iRange+1]; 1841 if (upper == lower) { 1842 value = upper; 1843 } else { 1844 // set correctly 1845 if (fabs(value  lower) <= primalTolerance * 1.001) { 1846 value = CoinMin(value, lower + primalTolerance); 1847 } else if (fabs(value  upper) <= primalTolerance * 1.001) { 1848 value = CoinMax(value, upper  primalTolerance); 1849 } else { 1850 //printf("*** variable wandered off bound %g %g %g!\n", 1851 // lower,value,upper); 1852 if (value  lower <= upper  value) 1853 value = lower + primalTolerance; 1854 else 1855 value = upper  primalTolerance; 1856 } 1857 } 1858 difference = cost  cost_[iRange]; 1859 cost = cost_[iRange]; 1860 } 1861 if (CLP_METHOD2) { 1862 double * upper = model_>upperRegion(); 1863 double * lower = model_>lowerRegion(); 1864 double * cost = model_>costRegion(); 1865 unsigned char iStatus = status_[iSequence]; 1866 assert (currentStatus(iStatus) == CLP_SAME); 1867 double lowerValue = lower[iSequence]; 1868 double upperValue = upper[iSequence]; 1869 double costValue = cost2_[iSequence]; 1870 // Set perceived direction out 1871 if (value <= lowerValue + 1.001 * primalTolerance) { 1872 direction = 1; 1873 } else if (value >= upperValue  1.001 * primalTolerance) { 1874 direction = 1; 1875 } else { 1876 // odd 1877 direction = 0; 1878 } 1879 int iWhere = originalStatus(iStatus); 1880 if (iWhere == CLP_BELOW_LOWER) { 1881 lowerValue = upperValue; 1882 upperValue = bound_[iSequence]; 1883 numberInfeasibilities_; 1884 assert (fabs(lowerValue) < 1.0e100); 1885 } else if (iWhere == CLP_ABOVE_UPPER) { 1886 upperValue = lowerValue; 1887 lowerValue = bound_[iSequence]; 1888 numberInfeasibilities_; 1889 } 1890 // get correct place 1891 // If fixed give benefit of doubt 1892 if (lowerValue == upperValue) 1893 value = lowerValue; 1894 int newWhere = CLP_FEASIBLE; 1895 if (value  upperValue <= primalTolerance) { 1896 if (value  lowerValue >= primalTolerance) { 1897 // feasible 1898 //newWhere=CLP_FEASIBLE; 1899 } else { 1900 // below 1901 newWhere = CLP_BELOW_LOWER; 1902 costValue = infeasibilityWeight_; 1903 numberInfeasibilities_++; 1904 assert (fabs(lowerValue) < 1.0e100); 1905 } 1906 } else { 1907 // above 1908 newWhere = CLP_ABOVE_UPPER; 1909 costValue += infeasibilityWeight_; 1910 numberInfeasibilities_++; 1911 } 1912 if (iWhere != newWhere) { 1913 difference = cost[iSequence]  costValue; 1914 setOriginalStatus(status_[iSequence], newWhere); 1915 if (newWhere == CLP_BELOW_LOWER) { 1916 bound_[iSequence] = upperValue; 1917 upper[iSequence] = lowerValue; 1918 lower[iSequence] = COIN_DBL_MAX; 1919 } else if (newWhere == CLP_ABOVE_UPPER) { 1920 bound_[iSequence] = lowerValue; 1921 lower[iSequence] = upperValue; 1922 upper[iSequence] = COIN_DBL_MAX; 1923 } else { 1924 lower[iSequence] = lowerValue; 1925 upper[iSequence] = upperValue; 1926 } 1927 cost[iSequence] = costValue; 1928 } 1929 // set correctly 1930 if (fabs(value  lowerValue) <= primalTolerance * 1.001) { 1931 value = CoinMin(value, lowerValue + primalTolerance); 1932 } else if (fabs(value  upperValue) <= primalTolerance * 1.001) { 1933 value = CoinMax(value, upperValue  primalTolerance); 1934 } else { 1935 //printf("*** variable wandered off bound %g %g %g!\n", 1936 // lowerValue,value,upperValue); 1937 if (value  lowerValue <= upperValue  value) 1938 value = lowerValue + primalTolerance; 1939 else 1940 value = upperValue  primalTolerance; 1941 } 1942 } 1943 changeCost_ += value * difference; 1807 } 1808 } 1809 } 1810 assert(iRange < end); 1811 whichRange_[iSequence] = iRange; 1812 if (iRange != currentRange) { 1813 if (infeasible(iRange)) 1814 numberInfeasibilities_++; 1815 if (infeasible(currentRange)) 1816 numberInfeasibilities_; 1817 } 1818 double &lower = model_>lowerAddress(iSequence); 1819 double &upper = model_>upperAddress(iSequence); 1820 double &cost = model_>costAddress(iSequence); 1821 lower = lower_[iRange]; 1822 upper = lower_[iRange + 1]; 1823 if (upper == lower) { 1824 value = upper; 1825 } else { 1826 // set correctly 1827 if (fabs(value  lower) <= primalTolerance * 1.001) { 1828 value = CoinMin(value, lower + primalTolerance); 1829 } else if (fabs(value  upper) <= primalTolerance * 1.001) { 1830 value = CoinMax(value, upper  primalTolerance); 1831 } else { 1832 //printf("*** variable wandered off bound %g %g %g!\n", 1833 // lower,value,upper); 1834 if (value  lower <= upper  value) 1835 value = lower + primalTolerance; 1836 else 1837 value = upper  primalTolerance; 1838 } 1839 } 1840 difference = cost  cost_[iRange]; 1841 cost = cost_[iRange]; 1842 } 1843 if (CLP_METHOD2) { 1844 double *upper = model_>upperRegion(); 1845 double *lower = model_>lowerRegion(); 1846 double *cost = model_>costRegion(); 1847 unsigned char iStatus = status_[iSequence]; 1848 assert(currentStatus(iStatus) == CLP_SAME); 1849 double lowerValue = lower[iSequence]; 1850 double upperValue = upper[iSequence]; 1851 double costValue = cost2_[iSequence]; 1852 // Set perceived direction out 1853 if (value <= lowerValue + 1.001 * primalTolerance) { 1854 direction = 1; 1855 } else if (value >= upperValue  1.001 * primalTolerance) { 1856 direction = 1; 1857 } else { 1858 // odd 1859 direction = 0; 1860 } 1861 int iWhere = originalStatus(iStatus); 1862 if (iWhere == CLP_BELOW_LOWER) { 1863 lowerValue = upperValue; 1864 upperValue = bound_[iSequence]; 1865 numberInfeasibilities_; 1866 assert(fabs(lowerValue) < 1.0e100); 1867 } else if (iWhere == CLP_ABOVE_UPPER) { 1868 upperValue = lowerValue; 1869 lowerValue = bound_[iSequence]; 1870 numberInfeasibilities_; 1871 } 1872 // get correct place 1873 // If fixed give benefit of doubt 1874 if (lowerValue == upperValue) 1875 value = lowerValue; 1876 int newWhere = CLP_FEASIBLE; 1877 if (value  upperValue <= primalTolerance) { 1878 if (value  lowerValue >= primalTolerance) { 1879 // feasible 1880 //newWhere=CLP_FEASIBLE; 1881 } else { 1882 // below 1883 newWhere = CLP_BELOW_LOWER; 1884 costValue = infeasibilityWeight_; 1885 numberInfeasibilities_++; 1886 assert(fabs(lowerValue) < 1.0e100); 1887 } 1888 } else { 1889 // above 1890 newWhere = CLP_ABOVE_UPPER; 1891 costValue += infeasibilityWeight_; 1892 numberInfeasibilities_++; 1893 } 1894 if (iWhere != newWhere) { 1895 difference = cost[iSequence]  costValue; 1896 setOriginalStatus(status_[iSequence], newWhere); 1897 if (newWhere == CLP_BELOW_LOWER) { 1898 bound_[iSequence] = upperValue; 1899 upper[iSequence] = lowerValue; 1900 lower[iSequence] = COIN_DBL_MAX; 1901 } else if (newWhere == CLP_ABOVE_UPPER) { 1902 bound_[iSequence] = lowerValue; 1903 lower[iSequence] = upperValue; 1904 upper[iSequence] = COIN_DBL_MAX; 1905 } else { 1906 lower[iSequence] = lowerValue; 1907 upper[iSequence] = upperValue; 1908 } 1909 cost[iSequence] = costValue; 1910 } 1911 // set correctly 1912 if (fabs(value  lowerValue) <= primalTolerance * 1.001) { 1913 value = CoinMin(value, lowerValue + primalTolerance); 1914 } else if (fabs(value  upperValue) <= primalTolerance * 1.001) { 1915 value = CoinMax(value, upperValue  primalTolerance); 1916 } else { 1917 //printf("*** variable wandered off bound %g %g %g!\n", 1918 // lowerValue,value,upperValue); 1919 if (value  lowerValue <= upperValue  value) 1920 value = lowerValue + primalTolerance; 1921 else 1922 value = upperValue  primalTolerance; 1923 } 1924 } 1925 changeCost_ += value * difference; 1944 1926 #ifdef CLP_USER_DRIVEN 1945 1946 if (differencesaveCost!=model_>costRegion()[iSequence]) {1947 1948 iSequence,iSequencemodel_>numberColumns(),1949 saveLower,value,saveUpper,saveCost,1950 1951 1952 1953 1954 #endif 1955 1927 //if (iSequence>=model_>numberColumns) 1928 if (difference  saveCost != model_>costRegion()[iSequence]) { 1929 printf("Out %d (%d) %g <= %g <= %g cost %g x=>x %g < %g cost %g\n", 1930 iSequence, iSequence  model_>numberColumns(), 1931 saveLower, value, saveUpper, saveCost, 1932 model_>lowerRegion()[iSequence], 1933 model_>upperRegion()[iSequence], 1934 model_>costRegion()[iSequence]); 1935 } 1936 #endif 1937 return direction; 1956 1938 } 1957 1939 // Returns nearest bound … … 1959 1941 ClpNonLinearCost::nearest(int iSequence, double solutionValue) 1960 1942 { 1961 assert(model_ != NULL);1962 1963 1964 1965 1966 1967 int end = start_[iSequence+1];1968 1969 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 const double *upper = model_>upperRegion();1981 const double *lower = model_>lowerRegion();1982 1983 1984 1985 1986 1987 1988 assert(fabs(lowerValue) < 1.0e100);1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1943 assert(model_ != NULL); 1944 double nearest = 0.0; 1945 if (CLP_METHOD1) { 1946 // get where in bound sequence 1947 int iRange; 1948 int start = start_[iSequence]; 1949 int end = start_[iSequence + 1]; 1950 int jRange = 1; 1951 nearest = COIN_DBL_MAX; 1952 for (iRange = start; iRange < end; iRange++) { 1953 if (fabs(solutionValue  lower_[iRange]) < nearest) { 1954 jRange = iRange; 1955 nearest = fabs(solutionValue  lower_[iRange]); 1956 } 1957 } 1958 assert(jRange < end); 1959 nearest = lower_[jRange]; 1960 } 1961 if (CLP_METHOD2) { 1962 const double *upper = model_>upperRegion(); 1963 const double *lower = model_>lowerRegion(); 1964 double lowerValue = lower[iSequence]; 1965 double upperValue = upper[iSequence]; 1966 int iWhere = originalStatus(status_[iSequence]); 1967 if (iWhere == CLP_BELOW_LOWER) { 1968 lowerValue = upperValue; 1969 upperValue = bound_[iSequence]; 1970 assert(fabs(lowerValue) < 1.0e100); 1971 } else if (iWhere == CLP_ABOVE_UPPER) { 1972 upperValue = lowerValue; 1973 lowerValue = bound_[iSequence]; 1974 } 1975 if (fabs(solutionValue  lowerValue) < fabs(solutionValue  upperValue)) 1976 nearest = lowerValue; 1977 else 1978 nearest = upperValue; 1979 } 1980 return nearest; 1999 1981 } 2000 1982 /// Feasible cost with offset and direction (i.e. for reporting) … … 2002 1984 ClpNonLinearCost::feasibleReportCost() const 2003 1985 { 2004 double value; 2005 model_>getDblParam(ClpObjOffset, value); 2006 return (feasibleCost_ + model_>objectiveAsObject()>nonlinearOffset()) * model_>optimizationDirection() / 2007 (model_>objectiveScale() * model_>rhsScale())  value; 1986 double value; 1987 model_>getDblParam(ClpObjOffset, value); 1988 return (feasibleCost_ + model_>objectiveAsObject()>nonlinearOffset()) * model_>optimizationDirection() / (model_>objectiveScale() * model_>rhsScale())  value; 2008 1989 } 2009 1990 // Get rid of real costs (just for moment) 2010 void 2011 ClpNonLinearCost::zapCosts() 2012 { 2013 int iSequence; 2014 double infeasibilityCost = model_>infeasibilityCost(); 2015 // zero out all costs 2016 int numberTotal = numberColumns_ + numberRows_; 2017 if (CLP_METHOD1) { 2018 int n = start_[numberTotal]; 2019 memset(cost_, 0, n * sizeof(double)); 2020 for (iSequence = 0; iSequence < numberTotal; iSequence++) { 2021 int start = start_[iSequence]; 2022 int end = start_[iSequence+1]  1; 2023 // correct costs for this infeasibility weight 2024 if (infeasible(start)) { 2025 cost_[start] = infeasibilityCost; 2026 } 2027 if (infeasible(end  1)) { 2028 cost_[end1] = infeasibilityCost; 2029 } 2030 } 2031 } 2032 if (CLP_METHOD2) { 2033 } 1991 void ClpNonLinearCost::zapCosts() 1992 { 1993 int iSequence; 1994 double infeasibilityCost = model_>infeasibilityCost(); 1995 // zero out all costs 1996 int numberTotal = numberColumns_ + numberRows_; 1997 if (CLP_METHOD1) { 1998 int n = start_[numberTotal]; 1999 memset(cost_, 0, n * sizeof(double)); 2000 for (iSequence = 0; iSequence < numberTotal; iSequence++) { 2001 int start = start_[iSequence]; 2002 int end = start_[iSequence + 1]  1; 2003 // correct costs for this infeasibility weight 2004 if (infeasible(start)) { 2005 cost_[start] = infeasibilityCost; 2006 } 2007 if (infeasible(end  1)) { 2008 cost_[end  1] = infeasibilityCost; 2009 } 2010 } 2011 } 2012 if (CLP_METHOD2) { 2013 } 2034 2014 } 2035 2015 #ifdef VALIDATE 2036 2016 // For debug 2037 void 2038 ClpNonLinearCost::validate() 2039 { 2040 double primalTolerance = model_>currentPrimalTolerance(); 2041 int iSequence; 2042 const double * solution = model_>solutionRegion(); 2043 const double * upper = model_>upperRegion(); 2044 const double * lower = model_>lowerRegion(); 2045 const double * cost = model_>costRegion(); 2046 double infeasibilityCost = model_>infeasibilityCost(); 2047 int numberTotal = numberRows_ + numberColumns_; 2048 int numberInfeasibilities = 0; 2049 double sumInfeasibilities = 0.0; 2050 2051 for (iSequence = 0; iSequence < numberTotal; iSequence++) { 2052 double value = solution[iSequence]; 2053 int iStatus = status_[iSequence]; 2054 assert (currentStatus(iStatus) == CLP_SAME); 2055 double lowerValue = lower[iSequence]; 2056 double upperValue = upper[iSequence]; 2057 double costValue = cost2_[iSequence]; 2058 int iWhere = originalStatus(iStatus); 2059 if (iWhere == CLP_BELOW_LOWER) { 2060 lowerValue = upperValue; 2061 upperValue = bound_[iSequence]; 2062 assert (fabs(lowerValue) < 1.0e100); 2063 costValue = infeasibilityCost; 2064 assert (value <= lowerValue  primalTolerance); 2065 numberInfeasibilities++; 2066 sumInfeasibilities += lowerValue  value  primalTolerance; 2067 assert (model_>getStatus(iSequence) == ClpSimplex::basic); 2068 } else if (iWhere == CLP_ABOVE_UPPER) { 2069 upperValue = lowerValue; 2070 lowerValue = bound_[iSequence]; 2071 costValue += infeasibilityCost; 2072 assert (value >= upperValue + primalTolerance); 2073 numberInfeasibilities++; 2074 sumInfeasibilities += value  upperValue  primalTolerance; 2075 assert (model_>getStatus(iSequence) == ClpSimplex::basic); 2076 } else { 2077 assert (value >= lowerValue  primalTolerance && value <= upperValue + primalTolerance); 2078 } 2079 assert (lowerValue == saveLowerV[iSequence]); 2080 assert (upperValue == saveUpperV[iSequence]); 2081 assert (costValue == cost[iSequence]); 2082 } 2083 if (numberInfeasibilities) 2084 printf("JJ %d infeasibilities summing to %g\n", 2085 numberInfeasibilities, sumInfeasibilities); 2086 } 2087 #endif 2017 void ClpNonLinearCost::validate() 2018 { 2019 double primalTolerance = model_>currentPrimalTolerance(); 2020 int iSequence; 2021 const double *solution = model_>solutionRegion(); 2022 const double *upper = model_>upperRegion(); 2023 const double *lower = model_>lowerRegion(); 2024 const double *cost = model_>costRegion(); 2025 double infeasibilityCost = model_>infeasibilityCost(); 2026 int numberTotal = numberRows_ + numberColumns_; 2027 int numberInfeasibilities = 0; 2028 double sumInfeasibilities = 0.0; 2029 2030 for (iSequence = 0; iSequence < numberTotal; iSequence++) { 2031 double value = solution[iSequence]; 2032 int iStatus = status_[iSequence]; 2033 assert(currentStatus(iStatus) == CLP_SAME); 2034 double lowerValue = lower[iSequence]; 2035 double upperValue = upper[iSequence]; 2036 double costValue = cost2_[iSequence]; 2037 int iWhere = originalStatus(iStatus); 2038 if (iWhere == CLP_BELOW_LOWER) { 2039 lowerValue = upperValue; 2040 upperValue = bound_[iSequence]; 2041 assert(fabs(lowerValue) < 1.0e100); 2042 costValue = infeasibilityCost; 2043 assert(value <= lowerValue  primalTolerance); 2044 numberInfeasibilities++; 2045 sumInfeasibilities += lowerValue  value  primalTolerance; 2046 assert(model_>getStatus(iSequence) == ClpSimplex::basic); 2047 } else if (iWhere == CLP_ABOVE_UPPER) { 2048 upperValue = lowerValue; 2049 lowerValue = bound_[iSequence]; 2050 costValue += infeasibilityCost; 2051 assert(value >= upperValue + primalTolerance); 2052 numberInfeasibilities++; 2053 sumInfeasibilities += value  upperValue  primalTolerance; 2054 assert(model_>getStatus(iSequence) == ClpSimplex::basic); 2055 } else { 2056 assert(value >= lowerValue  primalTolerance && value <= upperValue + primalTolerance); 2057 } 2058 assert(lowerValue == saveLowerV[iSequence]); 2059 assert(upperValue == saveUpperV[iSequence]); 2060 assert(costValue == cost[iSequence]); 2061 } 2062 if (numberInfeasibilities) 2063 printf("JJ %d infeasibilities summing to %g\n", 2064 numberInfeasibilities, sumInfeasibilities); 2065 } 2066 #endif 2067 2068 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2 2069 */
Note: See TracChangeset
for help on using the changeset viewer.