 Timestamp:
 Apr 16, 2003 11:20:06 AM (17 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/ClpSimplexPrimalQuadratic.cpp
r162 r163 13 13 #include "CoinIndexedVector.hpp" 14 14 #include "CoinWarmStartBasis.hpp" 15 #include "CoinMpsIO.hpp" 15 16 #include "ClpPrimalColumnPivot.hpp" 16 17 #include "ClpMessage.hpp" … … 108 109 last[iPass][jNon]=0; 109 110 } 110 double goodMove=false; 111 // goodMove +1 yes, 0 no, 1 last was bad  just halve gaps, 2 do nothing 112 int goodMove=2; 111 113 char * statusCheck = new char[numberColumns]; 112 114 for (iPass=0;iPass<numberPasses;iPass++) { … … 114 116 double offset=0.0; 115 117 double objValue=objectiveOffset; 118 double lambda=1.0; 119 if (goodMove>=0) { 120 // get best interpolation 121 double coeff0=objectiveOffset,coeff1=0.0,coeff2=0.0; 122 for (iColumn=0;iColumn<numberColumns;iColumn++) { 123 coeff0 += saveObjective[iColumn]*solution[iColumn]; 124 coeff1 += saveObjective[iColumn]*(saveSolution[iColumn]solution[iColumn]); 125 } 126 for (jNon=0;jNon<numberNonLinearColumns;jNon++) { 127 iColumn=listNonLinearColumn[jNon]; 128 double valueI = solution[iColumn]; 129 double valueISave = saveSolution[iColumn]; 130 int j; 131 for (j=columnQuadraticStart[iColumn]; 132 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) { 133 int jColumn = columnQuadratic[j]; 134 double valueJ = solution[jColumn]; 135 double valueJSave = saveSolution[jColumn]; 136 double elementValue = 0.5*quadraticElement[j]; 137 coeff0 += valueI*valueJ*elementValue; 138 coeff1 += (valueI*valueJSave+valueISave*valueJ2.0*valueI*valueJ)*elementValue; 139 coeff2 += (valueISave*valueJSave+valueI*valueJvalueISave*valueJvalueI*valueJSave)*elementValue; 140 } 141 } 142 double lambdaValue; 143 if (fabs(coeff2)<1.0e9) { 144 if (coeff1+coeff2>=0.0) 145 lambda = 0.0; 146 else 147 lambda = 1.0; 148 } else { 149 lambda = (0.5*coeff1)/coeff2; 150 if (lambda>1.0lambda<0.0) { 151 if (coeff1+coeff2>=0.0) 152 lambda = 0.0; 153 else 154 lambda = 1.0; 155 } 156 } 157 lambdaValue = lambda*lambda*coeff2+lambda*coeff1+coeff0; 158 printf("coeffs %g %g %g  lastobj %g\n",coeff0,coeff1,coeff2,lastObjective); 159 printf("obj at saved %g, obj now %g zero deriv at %g  value %g\n", 160 coeff0+coeff1+coeff2,coeff0,lambda,lambdaValue); 161 if (lambda>0.0&&lambda<=1.0) { 162 // update solution 163 for (iColumn=0;iColumn<numberColumns;iColumn++) 164 solution[iColumn] = lambda * saveSolution[iColumn] 165 + (1.0lambda) * solution[iColumn]; 166 if (lambda>0.999) { 167 memcpy(this>dualRowSolution(),savePi,numberRows*sizeof(double)); 168 memcpy(status_,saveStatus,numberRows+numberColumns); 169 } 170 if (lambda>0.99999&&fabs(coeff1+coeff2)>1.0e2) { 171 // tighten all 172 goodMove=1; 173 } 174 } 175 } 116 176 memcpy(objective,saveObjective,numberColumns*sizeof(double)); 117 177 for (iColumn=0;iColumn<numberColumns;iColumn++) … … 173 233 } 174 234 } 175 if (objValue<=lastObjective) 176 goodMove=true; 177 else 178 goodMove=false; 235 // goodMove +1 yes, 0 no, 1 last was bad  just halve gaps, 2 do nothing 179 236 double maxDelta=0.0; 237 if (goodMove>=0) { 238 if (objValue<=lastObjective) 239 goodMove=1; 240 else 241 goodMove=0; 242 } else { 243 maxDelta=1.0e10; 244 } 180 245 double maxGap=0.0; 181 246 for (jNon=0;jNon<numberNonLinearColumns;jNon++) { … … 183 248 maxDelta = max(maxDelta, 184 249 fabs(solution[iColumn]saveSolution[iColumn])); 185 if (goodMove ) {250 if (goodMove>0) { 186 251 if (last[0][jNon]*last[1][jNon]<0) { 187 252 // halve … … 192 257 trust[jNon] *= 1.5; 193 258 } 194 } else if ( trust[jNon]>10.0*deltaTolerance) {259 } else if (goodMove!=2&&trust[jNon]>10.0*deltaTolerance) { 195 260 trust[jNon] *= 0.5; 196 261 } … … 198 263 } 199 264 std::cout<<"largest gap is "<<maxGap<<std::endl; 200 if (goodMove ) {265 if (goodMove>0) { 201 266 double drop = lastObjectiveobjValue; 202 267 std::cout<<"True drop was "<<drop<<std::endl; 203 268 std::cout<<"largest delta is "<<maxDelta<<std::endl; 204 if (maxDelta<deltaTolerance&&drop<1.0e4&&goodMove ) {269 if (maxDelta<deltaTolerance&&drop<1.0e4&&goodMove&&lambda<0.99999) { 205 270 std::cout<<"Exiting"<<std::endl; 206 271 break; 207 272 } 208 } else { 209 //lastObjective += 1.0e3; // to stop exit 210 } 273 } 274 if (!iPass) 275 goodMove=1; 276 targetDrop=0.0; 277 double * r = this>dualColumnSolution(); 211 278 for (jNon=0;jNon<numberNonLinearColumns;jNon++) { 212 279 iColumn=listNonLinearColumn[jNon]; … … 218 285 trueUpper[jNon]); 219 286 } 287 if (iPass) { 288 // get reduced costs 289 this>matrix()>transposeTimes(savePi, 290 this>dualColumnSolution()); 291 for (jNon=0;jNon<numberNonLinearColumns;jNon++) { 292 iColumn=listNonLinearColumn[jNon]; 293 double dj = objective[iColumn]r[iColumn]; 294 r[iColumn]=dj; 295 if (dj<0.0) 296 targetDrop = dj*(columnUpper[iColumn]solution[iColumn]); 297 else 298 targetDrop = dj*(columnLower[iColumn]solution[iColumn]); 299 } 300 } else { 301 memset(r,0,numberColumns*sizeof(double)); 302 } 303 for (jNon=0;jNon<numberNonLinearColumns;jNon++) { 304 iColumn=listNonLinearColumn[jNon]; 305 if (statusCheck[iColumn]=='L'&&r[iColumn]<1.0e4) { 306 columnLower[iColumn]=max(solution[iColumn], 307 trueLower[jNon]); 308 columnUpper[iColumn]=min(solution[iColumn] 309 +trust[jNon], 310 trueUpper[jNon]); 311 } else if (statusCheck[iColumn]=='U'&&r[iColumn]>1.0e4) { 312 columnLower[iColumn]=max(solution[iColumn] 313 trust[jNon], 314 trueLower[jNon]); 315 columnUpper[iColumn]=min(solution[iColumn], 316 trueUpper[jNon]); 317 } else { 318 columnLower[iColumn]=max(solution[iColumn] 319 trust[jNon], 320 trueLower[jNon]); 321 columnUpper[iColumn]=min(solution[iColumn] 322 +trust[jNon], 323 trueUpper[jNon]); 324 } 325 } 220 326 if (goodMove) { 221 327 memcpy(saveSolution,solution,numberColumns*sizeof(double)); … … 223 329 memcpy(saveStatus,status_,numberRows+numberColumns); 224 330 225 targetDrop=0.0;226 if (iPass) {227 // get reduced costs228 this>matrix()>transposeTimes(savePi,229 this>dualColumnSolution());230 double * r = this>dualColumnSolution();231 for (jNon=0;jNon<numberNonLinearColumns;jNon++) {232 iColumn=listNonLinearColumn[jNon];233 double dj = objective[iColumn]r[iColumn];234 r[iColumn]=dj;235 if (dj<0.0)236 targetDrop = dj*(columnUpper[iColumn]solution[iColumn]);237 else238 targetDrop = dj*(columnLower[iColumn]solution[iColumn]);239 }240 }241 331 std::cout<<"Pass  "<<iPass 242 332 <<", target drop is "<<targetDrop … … 251 341 for (jNon=0;jNon<numberNonLinearColumns;jNon++) { 252 342 iColumn=listNonLinearColumn[jNon]; 253 printf("Trust %d %g  solution %d %g obj %g dj %g state %c \n",343 printf("Trust %d %g  solution %d %g obj %g dj %g state %c  bounds %g %g\n", 254 344 jNon,trust[jNon],iColumn,solution[iColumn],objective[iColumn], 255 r[iColumn],statusCheck[iColumn]); 345 r[iColumn],statusCheck[iColumn],columnLower[iColumn], 346 columnUpper[iColumn]); 256 347 } 257 348 } 258 349 setLogLevel(63); 350 this>scaling(false); 259 351 this>primal(1); 352 goodMove=1; 260 353 } else { 261 354 // bad pass  restore solution … … 265 358 memcpy(status_,saveStatus,numberRows+numberColumns); 266 359 iPass; 360 goodMove=1; 267 361 } 268 362 } … … 307 401 setDblParam(ClpObjOffset,objectiveOffsetoffset); 308 402 this>primal(1); 403 // redo values 404 setDblParam(ClpObjOffset,objectiveOffset); 405 objectiveValue_ += optimizationDirection_*offset; 309 406 for (jNon=0;jNon<numberNonLinearColumns;jNon++) { 310 407 iColumn=listNonLinearColumn[jNon]; … … 312 409 columnUpper[iColumn]= trueUpper[jNon]; 313 410 } 411 memcpy(objective,saveObjective,numberColumns*sizeof(double)); 314 412 delete [] saveSolution; 315 413 for (iPass=0;iPass<3;iPass++) … … 375 473 +2*quadratic>getNumElements() 376 474 + 2*numberColumns; 475 // type array 476 // >=0 points to dj or x (whichever this isn't) 477 // 2 row variable 478 // 1 artificial (see previous) 479 int * type = new int[newNumberColumns]; 377 480 double * elements2 = new double[numberElements]; 378 481 int * start2 = new int[newNumberColumns+1]; … … 451 554 for (iRow=0;iRow<numberRows;iRow++) { 452 555 // should look at rows to get correct bounds 453 columnLower2[newNumberColumns]=COIN_DBL_MAX; 454 columnUpper2[newNumberColumns]=COIN_DBL_MAX; 556 if (rowLower[iRow]==rowUpper[iRow]) { 557 columnLower2[newNumberColumns]=COIN_DBL_MAX; 558 columnUpper2[newNumberColumns]=COIN_DBL_MAX; 559 } else if (rowLower[iRow]<1.0e20) { 560 assert(rowUpper[iRow]<1.0e20); 561 columnLower2[newNumberColumns]=COIN_DBL_MAX; 562 columnUpper2[newNumberColumns]=0.0; 563 } else if (rowUpper[iRow]>1.0e20) { 564 columnLower2[newNumberColumns]=0.0; 565 columnUpper2[newNumberColumns]=COIN_DBL_MAX; 566 } else { 567 // can't do ranges just now 568 abort(); 569 } 455 570 int j; 456 571 for (j=rowStart[iRow]; … … 460 575 row2[numberElements++]=column[j]+numberRows; 461 576 } 577 type[newNumberColumns]=2; 462 578 newNumberColumns++; 463 579 start2[newNumberColumns]=numberElements; … … 467 583 for (iColumn=0;iColumn<numberColumns;iColumn++) { 468 584 // dj 469 columnLower2[newNumberColumns]= 0.0;585 columnLower2[newNumberColumns]=COIN_DBL_MAX; 470 586 columnUpper2[newNumberColumns]=COIN_DBL_MAX; 471 587 elements2[numberElements]=1.0; 472 588 row2[numberElements++]=iColumn+numberRows; 589 type[newNumberColumns]=iColumn; 590 type[iColumn]=newNumberColumns; 473 591 newNumberColumns++; 474 592 start2[newNumberColumns]=numberElements; 475 593 // artificial (assuming no bounds) 476 if (getStatus(iColumn)==basicsolution[iColumn]>1.0e7) { 594 if (getStatus(iColumn)==basic 595 (solution[iColumn]>columnLower[iColumn]+1.0e7&& 596 solution[iColumn]<columnUpper[iColumn]1.0e7)) { 597 columnUpper2[newNumberColumns1]=0.0; // fix for now 598 columnLower2[newNumberColumns1]=0.0; // fix for now 477 599 if (fabs(dj[iColumn])>tolerance) { 478 600 columnLower2[newNumberColumns]=0.0; … … 487 609 start2[newNumberColumns]=numberElements; 488 610 } 489 } else if (dj[iColumn]<tolerance) { 490 columnLower2[newNumberColumns]=0.0; 491 columnUpper2[newNumberColumns]=COIN_DBL_MAX; 492 objective2[newNumberColumns]=1.0; 493 elements2[numberElements]=1.0; 494 row2[numberElements++]=iColumn+numberRows; 495 newNumberColumns++; 496 start2[newNumberColumns]=numberElements; 497 } 498 } 499 // Create model 611 } else { 612 if (solution[iColumn]<columnLower[iColumn]+1.0e7) { 613 columnUpper2[iColumn]=columnLower[iColumn]; // fix for now 614 if (dj[iColumn]<tolerance) { 615 columnLower2[newNumberColumns]=0.0; 616 columnUpper2[newNumberColumns]=COIN_DBL_MAX; 617 objective2[newNumberColumns]=1.0; 618 elements2[numberElements]=1.0; 619 row2[numberElements++]=iColumn+numberRows; 620 newNumberColumns++; 621 start2[newNumberColumns]=numberElements; 622 } 623 } else { 624 columnLower2[iColumn]=columnUpper[iColumn]; // fix for now 625 if (dj[iColumn]>tolerance) { 626 columnLower2[newNumberColumns]=0.0; 627 columnUpper2[newNumberColumns]=COIN_DBL_MAX; 628 objective2[newNumberColumns]=1.0; 629 elements2[numberElements]=1.0; 630 row2[numberElements++]=iColumn+numberRows; 631 newNumberColumns++; 632 start2[newNumberColumns]=numberElements; 633 } 634 } 635 } 636 // temp 637 //columnLower2[iColumn]=solution[iColumn]; 638 //columnUpper2[iColumn]=solution[iColumn]; 639 } 640 // Create model 500 641 ClpSimplex model2; 501 642 model2.loadProblem(newNumberColumns,newNumberRows, … … 515 656 for (iColumn=0;iColumn<numberColumns;iColumn++) { 516 657 model2.setStatus(iColumn,getStatus(iColumn)); 517 if (getStatus(iColumn)==basic) { 658 if (getStatus(iColumn)==basic 659 (solution[iColumn]>columnLower[iColumn]+1.0e7&& 660 solution[iColumn]<columnUpper[iColumn]1.0e7)) { 518 661 solution2[newNumberColumns++]=0.0; 519 if (fabs(dj[iColumn])>tolerance) {662 if (fabs(dj[iColumn])>tolerance) 520 663 solution2[newNumberColumns++]=fabs(dj[iColumn]); 521 }522 664 } else if (dj[iColumn]<tolerance) { 523 665 solution2[newNumberColumns++]=0.0; … … 531 673 model2.primalRowSolution()); 532 674 // solve 675 #if 0 676 CoinMpsIO writer; 677 writer.setMpsData(*model2.matrix(), COIN_DBL_MAX, 678 model2.getColLower(), model2.getColUpper(), 679 model2.getObjCoefficients(), 680 (const char*) 0 /*integrality*/, 681 model2.getRowLower(), model2.getRowUpper(), 682 NULL,NULL); 683 writer.writeMps("xx.mps"); 684 #endif 685 model2.scaling(false); 533 686 model2.primal(1); 687 // If objective 0.0 then we can drop artificials 688 // See if any s sub j have wrong sign and/or use djs from infeasibility objective 689 double objectiveOffset; 690 getDblParam(ClpObjOffset,objectiveOffset); 691 double objValue = objectiveOffset; 692 for (iColumn=0;iColumn<numberColumns;iColumn++) 693 objValue += objective[iColumn]*solution2[iColumn]; 694 for (iColumn=0;iColumn<numberColumns;iColumn++) { 695 double valueI = solution2[iColumn]; 696 if (fabs(valueI)>1.0e5) 697 assert(solution2[type[iColumn]]<1.0e7); 698 int j; 699 for (j=columnQuadraticStart[iColumn]; 700 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) { 701 int jColumn = columnQuadratic[j]; 702 double valueJ = solution2[jColumn]; 703 double elementValue = quadraticElement[j]; 704 objValue += 0.5*valueI*valueJ*elementValue; 705 } 706 } 707 printf("Objective value %g\n",objValue); 708 for (iColumn=0;iColumn<newNumberColumns;iColumn++) 709 printf("%d %g\n",iColumn,solution2[iColumn]); 710 delete [] type; 534 711 return 0; 535 712 }
Note: See TracChangeset
for help on using the changeset viewer.