Ignore:
Timestamp:
Jul 15, 2008 5:56:58 PM (11 years ago)
Author:
forrest
Message:

take out parallel cuts

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Cbc/src/CbcCutGenerator.cpp

    r983 r1006  
    478478      int nEls=0;
    479479      int nCuts= numberRowCutsAfter-numberRowCutsBefore;
     480      // Remove NULL cuts!
     481      int nNull=0;
     482      const double * solution = solver->getColSolution();
     483      bool feasible=true;
     484      for (k = numberRowCutsAfter-1;k>=numberRowCutsBefore;k--) {
     485        const OsiRowCut * thisCut = cs.rowCutPtr(k) ;
     486        double sum=0.0;
     487        if (thisCut->lb()<=thisCut->ub()) {
     488          int n=thisCut->row().getNumElements();
     489          const int * column = thisCut->row().getIndices();
     490          const double * element = thisCut->row().getElements();
     491          if (n<=0) {
     492            // infeasible cut - give up
     493            feasible=false;
     494            break;
     495          }
     496          nEls+= n;
     497          for (int i=0;i<n;i++) {
     498            double value = element[i];
     499            sum += value*solution[column[i]];
     500          }
     501          if (sum>thisCut->ub()) {
     502            sum= sum-thisCut->ub();
     503          } else if (sum<thisCut->lb()) {
     504            sum= thisCut->lb()-sum;
     505          } else {
     506            sum=0.0;
     507            cs.eraseRowCut(k);
     508            nNull++;
     509          }
     510        }
     511      }
     512      //if (nNull)
     513      //printf("%s has %d cuts and %d elements - %d null!\n",generatorName_,
     514      //       nCuts,nEls,nNull);
     515      numberRowCutsAfter = cs.sizeRowCuts() ;
     516      nCuts= numberRowCutsAfter-numberRowCutsBefore;
     517      nEls=0;
    480518      for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
    481         OsiRowCut thisCut = cs.rowCut(k) ;
    482         int n=thisCut.row().getNumElements();
     519        const OsiRowCut * thisCut = cs.rowCutPtr(k) ;
     520        int n=thisCut->row().getNumElements();
    483521        nEls+= n;
    484522      }
     
    486524      //     nCuts,nEls);
    487525      int nElsNow = solver->getMatrixByCol()->getNumElements();
    488       if (nEls*8>nElsNow+8000+10000000) {
     526      int nAdd = model_->parentModel() ? 200 : 10000;
     527      int numberColumns = solver->getNumCols();
     528      int nAdd2 = model_->parentModel() ? 2*numberColumns : 5*numberColumns;
     529      if (/*nEls>CoinMax(nAdd2,nElsNow/8+nAdd)*/nCuts&&feasible) {
    489530        //printf("need to remove cuts\n");
    490531        // just add most effective
    491         int numberColumns = solver->getNumCols();
    492         int nReasonable = CoinMax(5*numberColumns,nElsNow/8);
     532        int nReasonable = CoinMax(nAdd2,nElsNow/8+nAdd);
    493533        int nDelete = nEls - nReasonable;
     534       
    494535        nElsNow = nEls;
    495         const double * solution = solver->getColSolution();
    496536        double * sort = new double [nCuts];
    497537        int * which = new int [nCuts];
     538        // For parallel cuts
     539        double * element2 = new double [numberColumns];
     540        CoinZeroN(element2,numberColumns);
    498541        for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
    499           OsiRowCut thisCut = cs.rowCut(k) ;
     542          const OsiRowCut * thisCut = cs.rowCutPtr(k) ;
    500543          double sum=0.0;
    501           if (thisCut.lb()<=thisCut.ub()) {
    502             int n=thisCut.row().getNumElements();
    503             const int * column = thisCut.row().getIndices();
    504             const double * element = thisCut.row().getElements();
     544          if (thisCut->lb()<=thisCut->ub()) {
     545            int n=thisCut->row().getNumElements();
     546            const int * column = thisCut->row().getIndices();
     547            const double * element = thisCut->row().getElements();
    505548            assert (n);
     549            double norm=0.0;
    506550            for (int i=0;i<n;i++) {
    507551              double value = element[i];
    508552              sum += value*solution[column[i]];
    509             }
    510             if (sum>thisCut.ub()) {
    511               sum= sum-thisCut.ub();
    512             } else if (sum<thisCut.lb()) {
    513               sum= thisCut.lb()-sum;
     553              norm += value*value;
     554            }
     555            if (sum>thisCut->ub()) {
     556              sum= sum-thisCut->ub();
     557            } else if (sum<thisCut->lb()) {
     558              sum= thisCut->lb()-sum;
    514559            } else {
    515               printf("ffffffffffffffff odd\n");
    516560              sum=0.0;
    517561            }
     562            // normalize
     563            sum /= sqrt(norm);
     564            // adjust for length
     565            //sum /= sqrt((double) n);
     566            // randomize
     567            //double randomNumber =
     568            //model_->randomNumberGenerator()->randomDouble();
     569            //sum *= (0.5+randomNumber);
    518570          } else {
    519571            // keep
     
    524576        }
    525577        CoinSort_2(sort,sort+nCuts,which);
     578        // Now see which ones are too similar
     579        int nParallel=0;
     580        for (k = 0;k<nCuts;k++) {
     581          int j=which[k];
     582          const OsiRowCut * thisCut = cs.rowCutPtr(j) ;
     583          if (thisCut->lb()>thisCut->ub())
     584            break; // cut is infeasible
     585          int n=thisCut->row().getNumElements();
     586          const int * column = thisCut->row().getIndices();
     587          const double * element = thisCut->row().getElements();
     588          assert (n);
     589          double norm=0.0;
     590          double lb = thisCut->lb();
     591          double ub = thisCut->ub();
     592          for (int i=0;i<n;i++) {
     593            double value = element[i];
     594            element2[column[i]]=value;
     595            norm += value*value;
     596          }
     597          int kkk = CoinMin(nCuts,k+5);
     598          for (int kk=k+1;kk<kkk;kk++) {
     599            int jj=which[kk];
     600            const OsiRowCut * thisCut2 = cs.rowCutPtr(jj) ;
     601            if (thisCut2->lb()>thisCut2->ub())
     602              break; // cut is infeasible
     603            int nB=thisCut2->row().getNumElements();
     604            const int * columnB = thisCut2->row().getIndices();
     605            const double * elementB = thisCut2->row().getElements();
     606            assert (nB);
     607            double normB=0.0;
     608            double product=0.0;
     609            for (int i=0;i<nB;i++) {
     610              double value = elementB[i];
     611              normB += value*value;
     612              product += value*element2[columnB[i]];
     613            }
     614            if (product>0.0&&product*product>0.99*norm*normB) {
     615              bool parallel=true;
     616              double lbB = thisCut2->lb();
     617              double ubB = thisCut2->ub();
     618              if ((lb<-1.0e20&&lbB>-1.0e20)||
     619                  (lbB<-1.0e20&&lb>-1.0e20))
     620                parallel = false;
     621              double tolerance;
     622              tolerance = CoinMax(fabs(lb),fabs(lbB))+1.0e-6;
     623              if (fabs(lb-lbB)>tolerance)
     624                parallel=false;
     625              if ((ub>1.0e20&&ubB<1.0e20)||
     626                  (ubB>1.0e20&&ub<1.0e20))
     627                parallel = false;
     628              tolerance = CoinMax(fabs(ub),fabs(ubB))+1.0e-6;
     629              if (fabs(ub-ubB)>tolerance)
     630                parallel=false;
     631              if (parallel) {
     632                nParallel++;
     633                sort[k]=0.0;
     634                break;
     635              }
     636            }
     637          }
     638          for (int i=0;i<n;i++) {
     639            element2[column[i]]=0.0;
     640          }
     641        }
     642        delete [] element2;
     643        CoinSort_2(sort,sort+nCuts,which);
    526644        k=0;
    527         while (nDelete>0) {
     645        while (nDelete>0||!sort[k]) {
    528646          int iCut=which[k];
    529           OsiRowCut thisCut = cs.rowCut(iCut) ;
    530           int n=thisCut.row().getNumElements();
     647          const OsiRowCut * thisCut = cs.rowCutPtr(iCut) ;
     648          int n=thisCut->row().getNumElements();
    531649          nDelete-=n;
    532650          k++;
     651          if (k>=nCuts)
     652            break;
    533653        }
    534654        std::sort(which,which+k);
     
    539659        delete [] sort;
    540660        delete [] which;
    541         int numberRowCutsAfter = cs.sizeRowCuts() ;
     661        numberRowCutsAfter = cs.sizeRowCuts() ;
     662#ifdef CLP_INVESTIGATE
    542663        nEls=0;
    543         nCuts= numberRowCutsAfter-numberRowCutsBefore;
     664        int nCuts2= numberRowCutsAfter-numberRowCutsBefore;
    544665        for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
    545           OsiRowCut thisCut = cs.rowCut(k) ;
    546           int n=thisCut.row().getNumElements();
     666          const OsiRowCut * thisCut = cs.rowCutPtr(k) ;
     667          int n=thisCut->row().getNumElements();
    547668          nEls+= n;
    548669        }
    549         printf("%s NOW has %d cuts and %d elements( down from %d els)\n",
    550                generatorName_,
    551                nCuts,nEls,nElsNow);
     670        if (!model_->parentModel()&&nCuts!=nCuts2)
     671          printf("%s NOW has %d cuts and %d elements( down from %d cuts and %d els) - %d parallel\n",
     672                 generatorName_,
     673                 nCuts2,nEls,nCuts,nElsNow,nParallel);
     674#endif
    552675      }
    553676    }
Note: See TracChangeset for help on using the changeset viewer.