Changeset 1222


Ignore:
Timestamp:
Sep 4, 2009 6:41:54 AM (11 years ago)
Author:
forrest
Message:

mods

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Cbc/examples/driver4b.cpp

    r1220 r1222  
    142142  bool sameSolution(double objValue,
    143143                    const double * solution);
     144  // creates cut
     145  void createCut(const double * solution);
    144146  // data goes here
    145147  // Best objective found
     
    315317                              const double * solution)
    316318{
    317   if (!solutions_)
    318     return false;
    319319  int numberColumns = model_->getNumCols();
    320320  double same=false;
    321   OsiSolverInterface * solver=model_->solver();
    322   // Optionally check for duplicates
    323   for (int j=0;j<maximumSolutions_;j++) {
    324     if (solutions_[j]) {
    325       double * sol = solutions_[j];
    326       if (fabs(sol[numberColumns]-objValue)<1.0e-5) {
    327         bool thisSame=true;
    328         for (int i=0;i<numberColumns;i++) {
    329           if (fabs(sol[i]-solution[i])>1.0e-6&&
    330               solver->isInteger(i)) {
    331             thisSame=false;
    332             break;
     321  if (solutions_) {
     322    OsiSolverInterface * solver=model_->solver();
     323    // Check for duplicates
     324    for (int j=0;j<maximumSolutions_;j++) {
     325      if (solutions_[j]) {
     326        double * sol = solutions_[j];
     327        if (fabs(sol[numberColumns]-objValue)<1.0e-5) {
     328          bool thisSame=true;
     329          for (int i=0;i<numberColumns;i++) {
     330            if (fabs(sol[i]-solution[i])>1.0e-6&&
     331                solver->isInteger(i)) {
     332              thisSame=false;
     333              break;
     334            }
     335          }
     336          if (thisSame) {
     337            printf("*** Solutions %d and %d same\n",
     338                   numberSolutions_+1,static_cast<int>(sol[numberColumns+1]));
     339            same=true;
    333340          }
    334341        }
    335         if (thisSame) {
    336           printf("*** Solutions %d and %d same\n",
    337                  numberSolutions_+1,static_cast<int>(sol[numberColumns+1]));
    338           same=true;
    339         }
     342      }
     343    }
     344  }
     345  if (!same) {
     346    printf("value of solution is %g\n",objValue);
     347    for (int i=0;i<numberColumns;i++) {
     348      if (fabs(solution[i])>1.0e-8) {
     349        printf("%d %g\n",i,solution[i]);
    340350      }
    341351    }
     
    343353  return same;
    344354}
    345 
    346355CbcEventHandler::CbcAction
    347356MyEventHandler3::event(CbcEvent whichEvent)
     
    353362      return stop; // say finished
    354363#endif
    355       // If preprocessing was done solution will be to processed model
    356       int numberColumns = model_->getNumCols();
    357364      const double * bestSolution = model_->bestSolution();
    358365      if (!bestSolution)
    359366        return noAction;
    360367      double objValue = model_->getObjValue();
    361       printf("value of solution is %g\n",objValue);
    362       for (int i=0;i<numberColumns;i++) {
    363         if (fabs(bestSolution[i])>1.0e-8) {
    364           printf("%d %g\n",i,bestSolution[i]);
    365         }
    366       }
    367368      /* if KILL_SOLUTION not set then should be as driver4 was
    368369         if 1 then all solutions rejected and cuts added
     
    375376#define KILL_SOLUTION 2
    376377      //#define FORCE_CUTS
    377 #ifdef KILL_SOLUTION
    378       if (!sameSolution(objValue,bestSolution))
     378#ifndef KILL_SOLUTION
     379      // If preprocessing was done solution will be to processed model
     380      int numberColumns = model_->getNumCols();
     381      printf("value of solution is %g\n",objValue);
     382      for (int i=0;i<numberColumns;i++) {
     383        if (fabs(bestSolution[i])>1.0e-8) {
     384          printf("%d %g\n",i,bestSolution[i]);
     385        }
     386      }
     387#else
     388      if (!sameSolution(objValue,bestSolution)) {
    379389        saveSolution((whichEvent==solution) ? "solution" : "heuristicSolution",
    380390                       objValue,bestSolution);
     391#if KILL_SOLUTION==2
     392        // Put in cut if we get here
     393        createCut(bestSolution);
     394#endif
     395      }
    381396      double newCutoff = bestSoFar_+0.1*fabs(bestSoFar_);
    382397      // Set cutoff here as may have been changed
     
    402417#if KILL_SOLUTION==1
    403418      // save as won't be done by 'solution'
    404       printf("value of solution is %g\n",objValue);
    405419      bool same = sameSolution(objValue,bestSolution);
    406420      if (!same)
     
    419433        return noAction;
    420434#endif
    421       // Get information on continuous solver
    422       OsiSolverInterface * solver = model_->continuousSolver();
    423       if (solver) {
    424         /* It is up to the user what sort of cuts to apply
    425            (also see FORCE_CUTS lower down)
    426            The classic one here, which just cuts off solution,
    427            does not do a lot - often the cuts are never applied!
    428         */
    429         int numberColumns=solver->getNumCols();
    430         const double * lower = solver->getColLower();
    431         const double * upper = solver->getColUpper();
    432         bool good=true;
    433         for (int iColumn=0;iColumn<numberColumns;iColumn++) {
    434           if (solver->isInteger(iColumn)) {
    435             double value=bestSolution[iColumn];
    436             if (value>lower[iColumn]+0.9&&
    437                 value<upper[iColumn]-0.9) {
    438               good=false;
    439               printf("Can't add cut as general integers\n");
    440               break;
    441             }
    442           }
    443         }
    444         if (good) {
    445           double * cut = new double [numberColumns];
    446           int * which = new int [numberColumns];
    447           // Could make -2.0 or .....
    448           double rhs=-1.0;
    449           int n=0;
    450           double sum=0.0;
    451           for (int iColumn=0;iColumn<numberColumns;iColumn++) {
    452             if (solver->isInteger(iColumn)) {
    453               if (upper[iColumn]>lower[iColumn]) {
    454                 double value=bestSolution[iColumn];
    455                 sum += value;
    456                 if (upper[iColumn]-value<0.1) {
    457                   rhs += upper[iColumn];
    458                   cut[n]=1.0;
    459                   which[n++]=iColumn;
    460                 } else {
    461                   rhs -= lower[iColumn];
    462                   cut[n]=-1.0;
    463                   which[n++]=iColumn;
    464                 }
    465               }
    466             }
    467           }
    468           assert (sum>rhs+0.9);
    469           printf("CUT has %d entries\n",n);
    470           OsiRowCut newCut;
    471           newCut.setLb(-COIN_DBL_MAX);
    472           newCut.setUb(rhs);
    473           newCut.setRow(n,which,cut,false);
    474 #ifdef FORCE_CUTS
    475           /* The next statement  is very optional as it will
    476              always add cut even if not violated.  If this is not done
    477              you can get duplicates from heuristics and strong branching.
    478              If you don't have it you may get some duplicates but it may be
    479              much faster. On one problem I tried I got over 300 cuts each
    480              of which was fully dense.
    481              Even if this is done, you can still get duplicates
    482              close together e.g. strong branching then heuristic
    483              before cut generation entered.*/
    484           newCut.setEffectiveness(COIN_DBL_MAX);
    485 #endif
    486           model_->makeGlobalCut(newCut);
     435      createCut(bestSolution);
    487436#if KILL_SOLUTION==1
    488           return killSolution;
    489 #endif
    490         }
    491       }
     437      return killSolution;
     438#endif
    492439#endif
    493440    } else if (whichEvent==endSearch) {
     
    541488}
    542489
     490// create cut
     491void
     492MyEventHandler3::createCut(const double * solution)
     493{
     494  /* TYPE_CUT 0 is off
     495     1 is row cut
     496     2 is column cut */
     497#define TYPE_CUT 1
     498#if TYPE_CUT
     499  // Get information on continuous solver
     500  OsiSolverInterface * solver = model_->continuousSolver();
     501  if (solver) {
     502    /* It is up to the user what sort of cuts to apply
     503       (also see FORCE_CUTS lower down)
     504       The classic one here, which just cuts off solution,
     505       does not do a lot - often the cuts are never applied!
     506    */
     507    int numberColumns=solver->getNumCols();
     508    const double * lower = solver->getColLower();
     509    const double * upper = solver->getColUpper();
     510    bool good=true;
     511    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     512      if (solver->isInteger(iColumn)) {
     513        double value=solution[iColumn];
     514        if (value>lower[iColumn]+0.9&&
     515            value<upper[iColumn]-0.9) {
     516          good=false;
     517          printf("Can't add cut as general integers\n");
     518          break;
     519        }
     520      }
     521    }
     522    if (good) {
     523#if TYPE_CUT==1
     524      double * cut = new double [numberColumns];
     525      int * which = new int [numberColumns];
     526      /* It is up to the user what sort of cuts to apply
     527         The classic row cut here, which just cuts off solution,
     528         does not do a lot - often the cuts are never applied!
     529      */
     530      // Should be -1.0 for classic cut
     531      double rhs=-6.0;
     532      int n=0;
     533      double sum=0.0;
     534      for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     535        if (solver->isInteger(iColumn)) {
     536          if (upper[iColumn]>lower[iColumn]) {
     537            double value=solution[iColumn];
     538            sum += value;
     539            if (upper[iColumn]-value<0.1) {
     540              rhs += upper[iColumn];
     541              cut[n]=1.0;
     542              which[n++]=iColumn;
     543            } else {
     544              rhs -= lower[iColumn];
     545              cut[n]=-1.0;
     546              which[n++]=iColumn;
     547            }
     548          }
     549        }
     550      }
     551      assert (sum>rhs+0.9);
     552      printf("CUT has %d entries\n",n);
     553      OsiRowCut newCut;
     554      newCut.setLb(-COIN_DBL_MAX);
     555      newCut.setUb(rhs);
     556      newCut.setRow(n,which,cut,false);
     557      delete [] which;
     558      delete [] cut;
     559#else
     560      // Fix a variable other way
     561      int firstAtLb=-1;
     562      int firstAtUb=-1;
     563      for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     564        if (solver->isInteger(iColumn)) {
     565          if (upper[iColumn]>lower[iColumn]) {
     566            double value=solution[iColumn];
     567            if (upper[iColumn]-value<0.1) {
     568              if (firstAtUb<0)
     569                firstAtUb=iColumn;
     570            } else {
     571              if (firstAtLb<0)
     572                firstAtLb=iColumn;
     573            }
     574          }
     575        }
     576      }
     577      OsiColCut newCut;
     578      CoinPackedVector bounds;
     579      if (firstAtUb>=0) {
     580        // fix to lb
     581        bounds.insert(firstAtUb,lower[firstAtUb]);
     582        newCut.setUbs(bounds);
     583      } else {
     584        // fix to ub
     585        bounds.insert(firstAtLb,upper[firstAtLb]);
     586        newCut.setLbs(bounds);
     587      }
     588#endif
     589#ifdef FORCE_CUTS
     590      /* The next statement  is very optional as it will
     591         always add cut even if not violated.  If this is not done
     592         you can get duplicates from heuristics and strong branching.
     593         If you don't have it you may get some duplicates but it may be
     594         much faster. On one problem I tried I got over 300 cuts each
     595         of which was fully dense.
     596         This does not apply to column cuts, which might as well
     597         always be forced on.
     598         
     599         Even if this is done, you can still get duplicates
     600         close together e.g. strong branching then heuristic
     601         before cut generation entered.*/
     602      newCut.setEffectiveness(COIN_DBL_MAX);
     603#endif
     604      model_->makeGlobalCut(newCut);
     605    }
     606  }
     607#endif
     608}
     609
    543610/** Tuning class
    544611 */
Note: See TracChangeset for help on using the changeset viewer.