Changeset 38 for trunk/CbcFathomDynamicProgramming.cpp
 Timestamp:
 Dec 10, 2004 2:53:37 PM (15 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/CbcFathomDynamicProgramming.cpp
r36 r38 15 15 #include "CoinHelperFunctions.hpp" 16 16 #include "CoinPackedMatrix.hpp" 17 #include "CoinSort.hpp" 17 18 // Default Constructor 18 19 CbcFathomDynamicProgramming::CbcFathomDynamicProgramming() … … 25 26 lookup_(NULL), 26 27 numberActive_(0), 28 maximumSizeAllowed_(1000000), 27 29 startBit_(NULL), 28 30 numberBits_(NULL), 29 rhs_(NULL) 31 rhs_(NULL), 32 numberNonOne_(0) 30 33 { 31 34 … … 40 43 lookup_(NULL), 41 44 numberActive_(0), 45 maximumSizeAllowed_(1000000), 42 46 startBit_(NULL), 43 47 numberBits_(NULL), 44 rhs_(NULL) 48 rhs_(NULL), 49 numberNonOne_(0) 45 50 { 46 51 type_=gutsOfCheckPossible(); … … 67 72 id_ = NULL; 68 73 lookup_ = NULL; 69 numberActive_ = 0;70 74 startBit_ = NULL; 71 75 numberBits_ = NULL; … … 90 94 lookup_(NULL), 91 95 numberActive_(rhs.numberActive_), 96 maximumSizeAllowed_(rhs.maximumSizeAllowed_), 92 97 startBit_(NULL), 93 98 numberBits_(NULL), 94 rhs_(NULL) 99 rhs_(NULL), 100 numberNonOne_(rhs.numberNonOne_) 95 101 { 96 102 if (size_) { … … 273 279 } 274 280 delete [] rhs; 281 if (allowableSize&&size_>allowableSize) { 282 printf("Too large  need %d entries x 16 bytes\n",size_); 283 return 1; // too big 284 } 275 285 if (n01==numberColumns&&!nbadcoeff) 276 286 return 0; // easiest … … 290 300 { 291 301 int returnCode=0; 292 int type=gutsOfCheckPossible( 1000000);302 int type=gutsOfCheckPossible(maximumSizeAllowed_); 293 303 if (type>=0) { 294 304 bool gotSolution=false; … … 309 319 310 320 int numberColumns = solver>getNumCols(); 321 int * indices = new int [numberActive_]; 322 double offset; 323 solver>getDblParam(OsiObjOffset,offset); 324 double fixedObj=offset; 325 int i; 311 326 // may be possible 312 // for test  just if all 1313 327 if (type==0) { 314 int i; 315 316 int * indices = new int [numberActive_]; 317 double offset; 318 solver>getDblParam(OsiObjOffset,offset); 319 double fixedObj=offset; 328 // rhs 1 and coefficients 1 320 329 for (i=0;i<numberColumns;i++) { 321 330 int j; … … 343 352 } 344 353 } 345 int needed=0; 354 returnCode=1; 355 } else { 356 // more complex 357 int * coefficients = new int[numberActive_]; 358 // If not too many general rhs then we can be more efficient 359 numberNonOne_=0; 360 for (i=0;i<numberActive_;i++) { 361 if (rhs_[i]!=1) 362 numberNonOne_++; 363 } 364 if (numberNonOne_*2<numberActive_) { 365 // put rhs >1 every second 366 int * permute = new int[numberActive_]; 367 int * temp = new int[numberActive_]; 368 // try different ways 369 int k=0; 370 for (i=0;i<numberRows;i++) { 371 int newRow = lookup_[i]; 372 if (newRow>=0&&rhs_[newRow]>1) { 373 permute[newRow]=k; 374 k +=2; 375 } 376 } 377 // adjust so k points to last 378 k = 2; 379 // and now rest 380 int k1=1; 381 for (i=0;i<numberRows;i++) { 382 int newRow = lookup_[i]; 383 if (newRow>=0&&rhs_[newRow]==1) { 384 permute[newRow]=k1; 385 k1++; 386 if (k1<=k) 387 k1++; 388 } 389 } 390 for (i=0;i<numberActive_;i++) { 391 int put = permute[i]; 392 temp[put]=rhs_[i]; 393 } 394 memcpy(rhs_,temp,numberActive_*sizeof(int)); 395 for (i=0;i<numberActive_;i++) { 396 int put = permute[i]; 397 temp[put]=numberBits_[i]; 398 } 399 memcpy(numberBits_,temp,numberActive_*sizeof(int)); 400 k=0; 401 for (i=0;i<numberActive_;i++) { 402 startBit_[i]=k; 403 k+= numberBits_[i]; 404 } 405 for (i=0;i<numberRows;i++) { 406 int newRow = lookup_[i]; 407 if (newRow>=0) 408 lookup_[i]=permute[newRow]; 409 } 410 delete [] permute; 411 delete [] temp; 412 // mark new method 413 type=2; 414 } 415 for (i=0;i<numberColumns;i++) { 416 int j; 417 double lowerValue = lower[i]; 418 assert (lowerValue==floor(lowerValue)); 419 double cost = direction * objective[i]; 420 fixedObj += lowerValue*cost; 421 if (lowerValue==upper[i]) 422 continue; 423 int n=0; 424 double gap=upper[i]lowerValue; 425 for (j=columnStart[i]; 426 j<columnStart[i]+columnLength[i];j++) { 427 int iRow=row[j]; 428 double value = element[j]; 429 int iValue = (int) value; 430 int newRow = lookup_[iRow]; 431 if (newRow<0iValue>rhs_[newRow]) { 432 n=0; 433 break; //can't use 434 } else { 435 coefficients[n]=iValue; 436 indices[n++]=newRow; 437 if (gap*iValue>rhs_[newRow]) { 438 gap = rhs_[newRow]/iValue; 439 } 440 } 441 } 442 int nTotal = (int) gap; 443 if (n) { 444 if (type==1) { 445 for (int k=1;k<=nTotal;k++) { 446 addOneColumn1(i,n,indices,coefficients,cost); 447 } 448 } else { 449 // may be able to take sort out with new method 450 CoinSort_2(indices,indices+n,coefficients); 451 for (int k=1;k<=nTotal;k++) { 452 addOneColumn1A(i,n,indices,coefficients,cost); 453 } 454 } 455 } 456 } 457 delete [] coefficients; 458 returnCode=1; 459 } 460 int needed=0; 461 double bestValue=COIN_DBL_MAX; 462 int iBest=1; 463 if (type==0) { 346 464 int numberActive=0; 347 465 for (i=0;i<numberRows;i++) { … … 350 468 if (rowLower[i]==rowUpper[i]) { 351 469 needed += 1<<numberActive; 352 } 353 numberActive++; 354 } 355 } 356 double bestValue=COIN_DBL_MAX; 357 int iBest=1; 470 numberActive++; 471 } 472 } 473 } 358 474 for (i=0;i<size_;i++) { 359 475 if ((i&needed)==needed) { … … 365 481 } 366 482 } 367 returnCode=1; 368 if (bestValue<COIN_DBL_MAX) { 369 bestValue += fixedObj; 370 printf("Can get solution of %g\n",bestValue); 371 if (bestValue<model_>getMinimizationObjValue()) { 372 // set up solution 373 betterSolution = new double[numberColumns]; 374 memcpy(betterSolution,lower,numberColumns*sizeof(double)); 375 while (iBest>0) { 376 int iColumn = id_[iBest]; 377 assert (iColumn>=0); 378 betterSolution[iColumn]++; 379 assert (betterSolution[iColumn]<=upper[iColumn]); 380 iBest = back_[iBest]; 381 } 382 } else { 383 } 384 } 385 delete [] indices; 386 } 483 } else { 484 int * lower = new int[numberActive_]; 485 for (i=0;i<numberRows;i++) { 486 int newRow = lookup_[i]; 487 if (newRow>=0) { 488 int gap=(int) (rowUpper[i]CoinMax(0.0,rowLower[i])); 489 lower[newRow]=rhs_[newRow]gap; 490 int numberBits = numberBits_[newRow]; 491 int startBit = startBit_[newRow]; 492 if (numberBits==1&&!gap) { 493 needed = 1<<startBit; 494 } 495 } 496 } 497 for (i=0;i<size_;i++) { 498 if ((i&needed)==needed) { 499 // this one may do 500 bool good = true; 501 for (int kk=0;kk<numberActive_;kk++) { 502 int numberBits = numberBits_[kk]; 503 int startBit = startBit_[kk]; 504 int size = 1<<numberBits; 505 int start = 1<<startBit; 506 int mask = start*(size1); 507 int level=(i&mask)>>startBit; 508 if (level<lower[kk]) { 509 good=false; 510 break; 511 } 512 } 513 if (good&&cost_[i]<bestValue) { 514 bestValue=cost_[i]; 515 iBest=i; 516 } 517 } 518 } 519 delete [] lower; 520 } 521 if (bestValue<COIN_DBL_MAX) { 522 bestValue += fixedObj; 523 printf("Can get solution of %g\n",bestValue); 524 if (bestValue<model_>getMinimizationObjValue()) { 525 // set up solution 526 betterSolution = new double[numberColumns]; 527 memcpy(betterSolution,lower,numberColumns*sizeof(double)); 528 while (iBest>0) { 529 int iColumn = id_[iBest]; 530 assert (iColumn>=0); 531 betterSolution[iColumn]++; 532 assert (betterSolution[iColumn]<=upper[iColumn]); 533 iBest = back_[iBest]; 534 } 535 } 536 } 537 delete [] indices; 387 538 gutsOfDelete(); 388 539 if (gotSolution) { … … 436 587 mask = 1<<iRow; 437 588 } 438 i= 0;589 i=size_1mask; 439 590 bool touched = false; 440 while (i <size_) {591 while (i>=0) { 441 592 int kMask = i&mask; 442 593 if (kMask==0) { … … 453 604 } 454 605 } 455 i ++;606 i; 456 607 } else { 457 608 // we can skip some 458 int k=i; 459 int iBit=0; 460 k &= ~1; 461 while ((k&kMask)!=0) { 462 iBit++; 463 k &= ~(1<<iBit); 464 } 465 // onto next 466 k += 1<<(iBit+1); 609 int k=(i&~mask)1; 467 610 #ifdef CBC_DEBUG 468 for (int j=i +1;j<k;j++) {611 for (int j=i1;j>k;j) { 469 612 int jMask = j&mask; 470 613 assert (jMask!=0); … … 476 619 return touched; 477 620 } 621 /* Adds one attempt of one column of type 1, 622 returns true if was used in making any changes. 623 At present the user has to call it once for each possible value 624 */ 625 bool 626 CbcFathomDynamicProgramming::addOneColumn1(int id,int numberElements, const int * rows, 627 const int * coefficients, double cost) 628 { 629 /* build up masks. 630 a) mask for 1 rhs 631 b) mask for addition 632 c) mask so adding will overflow 633 d) individual masks 634 */ 635 int mask1=0; 636 int maskAdd=0; 637 int mask2=0; 638 int i; 639 int n2=0; 640 int mask[40]; 641 int nextbit[40]; 642 int adjust[40]; 643 assert (numberElements<=40); 644 for (i=0;i<numberElements;i++) { 645 int iRow=rows[i]; 646 int numberBits = numberBits_[iRow]; 647 int startBit = startBit_[iRow]; 648 if (numberBits==1) { 649 mask1 = 1<<startBit; 650 maskAdd = 1<<startBit; 651 mask2 = 1<<startBit; 652 } else { 653 int value=coefficients[i]; 654 int size = 1<<numberBits; 655 int start = 1<<startBit; 656 assert (value<size); 657 maskAdd = start*value; 658 int gap = sizerhs_[iRow]1; 659 assert (gap>=0); 660 int hi2=rhs_[iRow]value; 661 if (hi2<size1) 662 hi2++; 663 adjust[n2] = start*hi2; 664 mask2 += start*gap; 665 nextbit[n2]=startBit+numberBits; 666 mask[n2++] = start*(size1); 667 } 668 } 669 i=size_1maskAdd; 670 bool touched = false; 671 while (i>=0) { 672 int kMask = i&mask1; 673 if (kMask==0) { 674 bool good=true; 675 for (int kk=n21;kk>=0;kk) { 676 int iMask = mask[kk]; 677 int jMask = iMask&mask2; 678 int kkMask = iMask&i; 679 kkMask += jMask; 680 if (kkMask>iMask) { 681 // we can skip some 682 int k=(i&~iMask); 683 k = adjust[kk]; 684 #ifdef CBC_DEBUG 685 for (int j=i1;j>k;j) { 686 int jMask = j&mask1; 687 if (jMask==0) { 688 bool good=true; 689 for (int kk=n21;kk>=0;kk) { 690 int iMask = mask[kk]; 691 int jMask = iMask&mask2; 692 int kkMask = iMask&i; 693 kkMask += jMask; 694 if (kkMask>iMask) { 695 good=false; 696 break; 697 } 698 } 699 assert (!good); 700 } 701 } 702 #endif 703 i=k; 704 good=false; 705 break; 706 } 707 } 708 if (good) { 709 double thisCost = cost_[i]; 710 if (thisCost!=COIN_DBL_MAX) { 711 // possible 712 double newCost=thisCost+cost; 713 int next = i + maskAdd; 714 if (cost_[next]>newCost) { 715 cost_[next]=newCost; 716 back_[next]=i; 717 id_[next]=id; 718 touched=true; 719 } 720 } 721 } 722 i; 723 } else { 724 // we can skip some 725 // we can skip some 726 int k=(i&~mask1)1; 727 #ifdef CBC_DEBUG 728 for (int j=i1;j>k;j) { 729 int jMask = j&mask; 730 assert (jMask!=0); 731 } 732 #endif 733 i=k; 734 } 735 } 736 return touched; 737 } 738 /* Adds one attempt of one column of type 1, 739 returns true if was used in making any changes. 740 At present the user has to call it once for each possible value 741 This version is when there are enough 1 rhs to do faster 742 */ 743 bool 744 CbcFathomDynamicProgramming::addOneColumn1A(int id,int numberElements, const int * rows, 745 const int * coefficients, double cost) 746 { 747 /* build up masks. 748 a) mask for 1 rhs 749 b) mask for addition 750 c) mask so adding will overflow 751 d) mask for non 1 rhs 752 */ 753 int maskA=0; 754 int maskAdd=0; 755 int maskC=0; 756 int maskD=0; 757 int i; 758 for (i=0;i<numberElements;i++) { 759 int iRow=rows[i]; 760 int numberBits = numberBits_[iRow]; 761 int startBit = startBit_[iRow]; 762 if (numberBits==1) { 763 maskA = 1<<startBit; 764 maskAdd = 1<<startBit; 765 } else { 766 int value=coefficients[i]; 767 int size = 1<<numberBits; 768 int start = 1<<startBit; 769 assert (value<size); 770 maskAdd = start*value; 771 int gap = sizerhs_[iRow]1; 772 assert (gap>=0); 773 maskC = start*gap; 774 maskD = start*(size1); 775 } 776 } 777 int maskDiff = maskDmaskC; 778 i=size_1maskAdd; 779 bool touched = false; 780 while (i>=0) { 781 int kMask = i&maskA; 782 if (kMask==0) { 783 int added = i & maskD; // just bits belonging to non 1 rhs 784 added += maskC; // will overflow mask if bad 785 added &= (~maskD); 786 if (added == 0) { 787 double thisCost = cost_[i]; 788 if (thisCost!=COIN_DBL_MAX) { 789 // possible 790 double newCost=thisCost+cost; 791 int next = i + maskAdd; 792 if (cost_[next]>newCost) { 793 cost_[next]=newCost; 794 back_[next]=i; 795 id_[next]=id; 796 touched=true; 797 } 798 } 799 i; 800 } else { 801 // we can skip some 802 int k = i & ~ maskD; // clear all 803 // Put back enough  but only below where we are 804 int kk=(numberNonOne_<<1)2; 805 assert (rhs_[kk]>1); 806 int iMask=0; 807 for(;kk>=0;kk=2) { 808 iMask = 1<<startBit_[kk+1]; 809 if ((added&iMask)!=0) { 810 iMask; 811 break; 812 } 813 } 814 assert (kk>=0); 815 iMask &= maskDiff; 816 k = iMask; 817 assert (k<i); 818 i=k; 819 } 820 } else { 821 // we can skip some 822 // we can skip some 823 int k=(i&~maskA)1; 824 #ifdef CBC_DEBUG 825 for (int j=i1;j>k;j) { 826 int jMask = j&mask; 827 assert (jMask!=0); 828 } 829 #endif 830 i=k; 831 } 832 } 833 return touched; 834 } 478 835 // update model 479 836 void CbcFathomDynamicProgramming::setModel(CbcModel * model)
Note: See TracChangeset
for help on using the changeset viewer.