Changeset 862


Ignore:
Timestamp:
Jan 22, 2008 4:23:57 PM (12 years ago)
Author:
forrest
Message:

try and make faster

Location:
trunk/Cbc/src
Files:
3 edited

Legend:

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

    r857 r862  
    11941194    }
    11951195  }
     1196
    11961197  solverCharacteristics_->setSolver(solver_);
    11971198  // Set so we can tell we are in initial phase in resolve
     
    25172518            double heurValue = getCutoff() ;
    25182519            int iHeur ;
    2519             for (iHeur = 0 ; iHeur < numberHeuristics_ ; iHeur++)
    2520             { double saveValue = heurValue ;
     2520            for (iHeur = 0 ; iHeur < numberHeuristics_ ; iHeur++) {
     2521              double saveValue = heurValue ;
    25212522              int ifSol = heuristic_[iHeur]->solution(heurValue,newSolution) ;
    25222523              if (ifSol > 0) {
     
    25242525                found = iHeur ;
    25252526                incrementUsed(newSolution);
    2526               }
    2527               else
    2528               if (ifSol < 0)    // just returning an estimate
    2529               { estValue = CoinMin(heurValue,estValue) ;
    2530                 heurValue = saveValue ; } }
    2531             if (found >= 0) {
    2532               lastHeuristic_ = heuristic_[found];
    2533               setBestSolution(CBC_ROUNDING,heurValue,newSolution) ;
    2534             }
     2527                lastHeuristic_ = heuristic_[found];
     2528                setBestSolution(CBC_ROUNDING,heurValue,newSolution) ;
     2529              } else if (ifSol < 0) {
     2530                // just returning an estimate
     2531                estValue = CoinMin(heurValue,estValue) ;
     2532                heurValue = saveValue ;
     2533              }
     2534            }
    25352535            delete [] newSolution ;
    25362536            newNode->setGuessedObjectiveValue(estValue) ;
     
    51085108      numberFixed++ ; } } }
    51095109  numberDJFixed_ += numberFixed;
    5110  
    51115110  return numberFixed; }
    51125111
     
    56365635                printf("Cut generator %d (%s) produced invalid cut (%dth in this go)\n",
    56375636                       i,generator_[i]->cutGeneratorName(),k-numberRowCutsBefore);
     5637                const double *lower = getColLower() ;
     5638                const double *upper = getColUpper() ;
     5639                int numberColumns = solver_->getNumCols();
     5640                for (int i=0;i<numberColumns;i++)
     5641                  printf("%d bounds %g,%g\n",i,lower[i],upper[i]);
    56385642                abort();
    56395643#endif
     
    58975901                printf("Cut generator %d (%s) produced invalid cut (%dth in this go)\n",
    58985902                       i,generator_[i]->cutGeneratorName(),k-numberRowCutsBefore);
     5903                const double *lower = getColLower() ;
     5904                const double *upper = getColUpper() ;
     5905                int numberColumns = solver_->getNumCols();
     5906                for (int i=0;i<numberColumns;i++)
     5907                  printf("%d bounds %g,%g\n",i,lower[i],upper[i]);
    58995908                abort();
    59005909#endif
     
    59905999          found = i ;
    59916000          incrementUsed(newSolution);
     6001          lastHeuristic_ = heuristic_[found];
     6002          setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
    59926003        } else if (ifSol<0) {
    59936004          heuristicValue = saveValue ;
     
    60006011      if (found >= 0) {
    60016012        phase_=4;
    6002         incrementUsed(newSolution);
    6003         lastHeuristic_ = heuristic_[found];
    6004         setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
    60056013        CbcTreeLocal * tree
    60066014          = dynamic_cast<CbcTreeLocal *> (tree_);
     
    63546362  because a generator indicated we're infeasible.
    63556363*/
    6356   if (feasible && solver_->isProvenOptimal())
    6357     reducedCostFix();
     6364  if (feasible && solver_->isProvenOptimal()) 
     6365    reducedCostFix() ;
    63586366  // If at root node do heuristics
    63596367  if (!numberNodes_) {
     
    63936401        found = i ;
    63946402        incrementUsed(newSolution);
     6403        lastHeuristic_ = heuristic_[found];
     6404        setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
    63956405      } else {
    63966406        heuristicValue = saveValue ;
     
    64006410    if (found >= 0) {
    64016411      phase_=4;
    6402       incrementUsed(newSolution);
    6403       lastHeuristic_ = heuristic_[found];
    6404       setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
    64056412    }
    64066413    delete [] newSolution ;
     
    70227029  if (feasible)
    70237030    {
    7024       resolve(solver_) ;
    7025       numberIterations_ += solver_->getIterationCount() ;
    7026       feasible = (solver_->isProvenOptimal() &&
    7027                   !solver_->isDualObjectiveLimitReached()) ;
     7031      bool onOptimalPath=false;
     7032      if ((specialOptions_&1)!=0) {
     7033        const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger() ;
     7034        if (debugger) {
     7035          onOptimalPath=true;
     7036          printf("On optimal path d\n") ;
     7037        }
     7038      }
     7039      int nTightened=0;
     7040#if 1
     7041      if (!currentNode_||(currentNode_->depth()&2)!=0)
     7042        nTightened=tightenBounds();
     7043      if (nTightened) {
     7044        //printf("%d bounds tightened\n",nTightened);
     7045        if ((specialOptions_&1)!=0&&onOptimalPath) {
     7046          const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger() ;
     7047          if (!debugger) {
     7048            // tighten did something???
     7049            solver_->writeMps("infeas4");
     7050            printf("Not on optimalpath aaaa\n");
     7051            abort();
     7052          }
     7053        }
     7054      }
     7055#endif
     7056      if (nTightened>=0) {
     7057        resolve(solver_) ;
     7058        numberIterations_ += solver_->getIterationCount() ;
     7059        feasible = (solver_->isProvenOptimal() &&
     7060                    !solver_->isDualObjectiveLimitReached()) ;
     7061        if ((specialOptions_&1)!=0&&onOptimalPath) {
     7062          if (!solver_->getRowCutDebugger()) {
     7063            // tighten did something???
     7064            solver_->writeMps("infeas4");
     7065            assert (solver_->getRowCutDebugger()) ;
     7066            printf("Not on optimalpath e\n");
     7067            abort();
     7068          }
     7069        }
     7070      } else {
     7071        feasible=false;
     7072      }
    70287073    }
    70297074  if (0&&feasible) {
     
    91479192  setCutoff(saveCutoff);
    91489193  return true;
     9194}
     9195// Tighten bounds - lightweight
     9196int
     9197CbcModel::tightenBounds()
     9198{
     9199  //CoinPackedMatrix matrixByRow(*solver_->getMatrixByRow());
     9200  int numberRows = solver_->getNumRows();
     9201  int numberColumns = solver_->getNumCols();
     9202
     9203  int iRow,iColumn;
     9204
     9205  // Row copy
     9206  //const double * elementByRow = matrixByRow.getElements();
     9207  //const int * column = matrixByRow.getIndices();
     9208  //const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
     9209  //const int * rowLength = matrixByRow.getVectorLengths();
     9210
     9211  const double * columnUpper = solver_->getColUpper();
     9212  const double * columnLower = solver_->getColLower();
     9213  const double * rowUpper = solver_->getRowUpper();
     9214  const double * rowLower = solver_->getRowLower();
     9215
     9216  // Column copy of matrix
     9217  const double * element = solver_->getMatrixByCol()->getElements();
     9218  const int * row = solver_->getMatrixByCol()->getIndices();
     9219  const CoinBigIndex * columnStart = solver_->getMatrixByCol()->getVectorStarts();
     9220  const int * columnLength = solver_->getMatrixByCol()->getVectorLengths();
     9221  const double *objective = solver_->getObjCoefficients() ;
     9222  double direction = solver_->getObjSense();
     9223  double * down = new double [numberRows];
     9224  double * up = new double [numberRows];
     9225  double * sum = new double [numberRows];
     9226  int * type = new int [numberRows];
     9227  CoinZeroN(down,numberRows);
     9228  CoinZeroN(up,numberRows);
     9229  CoinZeroN(sum,numberRows);
     9230  CoinZeroN(type,numberRows);
     9231  double infinity = solver_->getInfinity();
     9232  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     9233    CoinBigIndex start = columnStart[iColumn];
     9234    CoinBigIndex end = start + columnLength[iColumn];
     9235    double lower = columnLower[iColumn];
     9236    double upper = columnUpper[iColumn];
     9237    if (lower==upper) {
     9238      for (CoinBigIndex j=start;j<end;j++) {
     9239        int iRow = row[j];
     9240        double value = element[j];
     9241        sum[iRow]+=2.0*fabs(value*lower);
     9242        if ((type[iRow]&1)==0)
     9243          down[iRow] += value*lower;
     9244        if ((type[iRow]&2)==0)
     9245          up[iRow] += value*lower;
     9246      }
     9247    } else {
     9248      for (CoinBigIndex j=start;j<end;j++) {
     9249        int iRow = row[j];
     9250        double value = element[j];
     9251        if (value>0.0) {
     9252          if ((type[iRow]&1)==0) {
     9253            if (lower!=-infinity) {
     9254              down[iRow] += value*lower;
     9255              sum[iRow]+=fabs(value*lower);
     9256            } else {
     9257              type[iRow] |= 1;
     9258            }
     9259          }
     9260          if ((type[iRow]&2)==0) {
     9261            if (upper!=infinity) {
     9262              up[iRow] += value*upper;
     9263              sum[iRow]+=fabs(value*upper);
     9264            } else {
     9265              type[iRow] |= 2;
     9266            }
     9267          }
     9268        } else {
     9269          if ((type[iRow]&1)==0) {
     9270            if (upper!=infinity) {
     9271              down[iRow] += value*upper;
     9272              sum[iRow]+=fabs(value*upper);
     9273            } else {
     9274              type[iRow] |= 1;
     9275            }
     9276          }
     9277          if ((type[iRow]&2)==0) {
     9278            if (lower!=-infinity) {
     9279              up[iRow] += value*lower;
     9280              sum[iRow]+=fabs(value*lower);
     9281            } else {
     9282              type[iRow] |= 2;
     9283            }
     9284          }
     9285        }
     9286      }
     9287    }
     9288  }
     9289  int nTightened=0;
     9290  double tolerance = 1.0e-6;
     9291  for (iRow=0;iRow<numberRows;iRow++) {
     9292    if ((type[iRow]&1)!=0)
     9293      down[iRow]=-infinity;
     9294    if (down[iRow]>rowUpper[iRow]) {
     9295      if (down[iRow]>rowUpper[iRow]+tolerance+1.0e-8*sum[iRow]) {
     9296        // infeasible
     9297#ifdef COIN_DEVELOP
     9298        printf("infeasible on row %d\n",iRow);
     9299#endif
     9300        nTightened=-1;
     9301        break;
     9302      } else {
     9303        down[iRow]=rowUpper[iRow];
     9304      }
     9305    }
     9306    if ((type[iRow]&2)!=0)
     9307      up[iRow]=infinity;
     9308    if (up[iRow]<rowLower[iRow]) {
     9309      if (up[iRow]<rowLower[iRow]-tolerance-1.0e-8*sum[iRow]) {
     9310        // infeasible
     9311#ifdef COIN_DEVELOP
     9312        printf("infeasible on row %d\n",iRow);
     9313#endif
     9314        nTightened=-1;
     9315        break;
     9316      } else {
     9317        up[iRow]=rowLower[iRow];
     9318      }
     9319    }
     9320  }
     9321  if (nTightened)
     9322    numberColumns=0; // so will skip
     9323  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     9324    CoinBigIndex start = columnStart[iColumn];
     9325    CoinBigIndex end = start + columnLength[iColumn];
     9326    double lower = columnLower[iColumn];
     9327    double upper = columnUpper[iColumn];
     9328    double gap = upper-lower;
     9329    if (!gap)
     9330      continue;
     9331    int canGo=0;
     9332    if (integerInfo_[iColumn]) {
     9333      if (lower!=floor(lower+0.5)) {
     9334#ifdef COIN_DEVELOP
     9335        printf("increasing lower bound on %d from %g to %g\n",iColumn,
     9336               lower,ceil(lower));
     9337#endif
     9338        lower=ceil(lower);
     9339        gap=upper-lower;
     9340        solver_->setColLower(iColumn,lower);
     9341      }
     9342      if (upper!=floor(upper+0.5)) {
     9343#ifdef COIN_DEVELOP
     9344        printf("decreasing upper bound on %d from %g to %g\n",iColumn,
     9345               upper,floor(upper));
     9346#endif
     9347        upper=floor(upper);
     9348        gap=upper-lower;
     9349        solver_->setColUpper(iColumn,upper);
     9350      }
     9351      double newLower=lower;
     9352      double newUpper=upper;
     9353      for (CoinBigIndex j=start;j<end;j++) {
     9354        int iRow = row[j];
     9355        double value = element[j];
     9356        if (value>0.0) {
     9357          if ((type[iRow]&1)==0) {
     9358            // has to be at most something
     9359            if (down[iRow] + value*gap > rowUpper[iRow]+tolerance) {
     9360              double newGap = (rowUpper[iRow]-down[iRow])/value;
     9361              // adjust
     9362              newGap += 1.0e-10*sum[iRow];
     9363              newGap = floor(newGap);
     9364              if (lower+newGap<newUpper)
     9365                newUpper=lower+newGap;
     9366            }
     9367          }
     9368          if (down[iRow]<rowLower[iRow])
     9369            canGo |=1; // can't go down without affecting result
     9370          if ((type[iRow]&2)==0) {
     9371            // has to be at least something
     9372            if (up[iRow] - value*gap < rowLower[iRow]-tolerance) {
     9373              double newGap = (up[iRow]-rowLower[iRow])/value;
     9374              // adjust
     9375              newGap += 1.0e-10*sum[iRow];
     9376              newGap = floor(newGap);
     9377              if (upper-newGap>newLower)
     9378                newLower=upper-newGap;
     9379            }
     9380          }
     9381          if (up[iRow]>rowUpper[iRow])
     9382            canGo |=2; // can't go up without affecting result
     9383        } else {
     9384          if ((type[iRow]&1)==0) {
     9385            // has to be at least something
     9386            if (down[iRow] - value*gap > rowUpper[iRow]+tolerance) {
     9387              double newGap = -(rowUpper[iRow]-down[iRow])/value;
     9388              // adjust
     9389              newGap += 1.0e-10*sum[iRow];
     9390              newGap = floor(newGap);
     9391              if (upper-newGap>newLower)
     9392                newLower=upper-newGap;
     9393            }
     9394          }
     9395          if (up[iRow]>rowUpper[iRow])
     9396            canGo |=1; // can't go down without affecting result
     9397          if ((type[iRow]&2)==0) {
     9398            // has to be at most something
     9399            if (up[iRow] + value*gap < rowLower[iRow]-tolerance) {
     9400              double newGap = -(up[iRow]-rowLower[iRow])/value;
     9401              // adjust
     9402              newGap += 1.0e-10*sum[iRow];
     9403              newGap = floor(newGap);
     9404              if (lower+newGap<newUpper)
     9405                newUpper=lower+newGap;
     9406            }
     9407          }
     9408          if (down[iRow]<rowLower[iRow])
     9409            canGo |=2; // can't go up without affecting result
     9410        }
     9411      }
     9412      if (newUpper<upper||newLower>lower) {
     9413        nTightened++;
     9414        if (newLower>newUpper) {
     9415          // infeasible
     9416#if COIN_DEVELOP>1
     9417          printf("infeasible on column %d\n",iColumn);
     9418#endif
     9419          nTightened=-1;
     9420          break;
     9421        } else {
     9422          solver_->setColLower(iColumn,newLower);
     9423          solver_->setColUpper(iColumn,newUpper);
     9424        }
     9425        for (CoinBigIndex j=start;j<end;j++) {
     9426          int iRow = row[j];
     9427          double value = element[j];
     9428          if (value>0.0) {
     9429            if ((type[iRow]&1)==0) {
     9430              down[iRow] += value*(newLower-lower);
     9431            }
     9432            if ((type[iRow]&2)==0) {
     9433              up[iRow] += value*(newUpper-upper);
     9434            }
     9435          } else {
     9436            if ((type[iRow]&1)==0) {
     9437              down[iRow] += value*(newUpper-upper);
     9438            }
     9439            if ((type[iRow]&2)==0) {
     9440              up[iRow] += value*(newLower-lower);
     9441            }
     9442          }
     9443        }
     9444      } else {
     9445        if (canGo!=3) {
     9446          double objValue = direction*objective[iColumn];
     9447          if (objValue>=0.0&&(canGo&1)==0) {
     9448#if COIN_DEVELOP>1
     9449            printf("dual fix down on column %d\n",iColumn);
     9450#endif
     9451            nTightened++;;
     9452            solver_->setColUpper(iColumn,lower);
     9453          } else if (objValue<=0.0&&(canGo&2)==0) {
     9454#if COIN_DEVELOP>1
     9455            printf("dual fix up on column %d\n",iColumn);
     9456#endif
     9457            nTightened++;;
     9458            solver_->setColLower(iColumn,upper);
     9459          }
     9460        }           
     9461      }
     9462    } else {
     9463      // just do dual tests
     9464      for (CoinBigIndex j=start;j<end;j++) {
     9465        int iRow = row[j];
     9466        double value = element[j];
     9467        if (value>0.0) {
     9468          if (down[iRow]<rowLower[iRow])
     9469            canGo |=1; // can't go down without affecting result
     9470          if (up[iRow]>rowUpper[iRow])
     9471            canGo |=2; // can't go up without affecting result
     9472        } else {
     9473          if (up[iRow]>rowUpper[iRow])
     9474            canGo |=1; // can't go down without affecting result
     9475          if (down[iRow]<rowLower[iRow])
     9476            canGo |=2; // can't go up without affecting result
     9477        }
     9478      }
     9479      if (canGo!=3) {
     9480        double objValue = direction*objective[iColumn];
     9481        if (objValue>=0.0&&(canGo&1)==0) {
     9482#if COIN_DEVELOP>1
     9483          printf("dual fix down on continuous column %d\n",iColumn);
     9484#endif
     9485          nTightened++;;
     9486          solver_->setColUpper(iColumn,lower);
     9487        } else if (objValue<=0.0&&(canGo&2)==0) {
     9488#if COIN_DEVELOP>1
     9489          printf("dual fix up on continuoe column %d\n",iColumn);
     9490#endif
     9491          nTightened++;;
     9492          solver_->setColLower(iColumn,upper);
     9493        }
     9494      }
     9495    }
     9496  }
     9497  delete [] type;
     9498  delete [] down;
     9499  delete [] up;
     9500  delete [] sum;
     9501  return nTightened;
    91499502}
    91509503/*
     
    1086811221        numberSolutions_++;
    1086911222        numberHeuristicSolutions_++;
     11223        lastHeuristic_ = heuristic_[i];
     11224        setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
    1087011225      } else {
    1087111226        heuristicValue = saveValue ;
     
    1087811233    */
    1087911234    if (found >= 0) {
    10880       // For compiler error on terra cluster!
    10881       if (found<numberHeuristics_)
    10882         lastHeuristic_ = heuristic_[found];
    10883       else
    10884         lastHeuristic_ = heuristic_[0];
    10885       setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
    1088611235      CbcTreeLocal * tree
    1088711236        = dynamic_cast<CbcTreeLocal *> (tree_);
     
    1279513144    CbcIntegerBranchingObject * objectI = dynamic_cast<CbcIntegerBranchingObject *> (obj);
    1279613145    //const OsiObject * object2 = obj->orig
     13146#ifndef NDEBUG
    1279713147    const CbcSimpleInteger * object2 = dynamic_cast<const CbcSimpleInteger *> (objectI->object());
    1279813148    assert (object2);
    1279913149    assert (iColumn == object2->columnNumber());
     13150#endif
    1280013151    double bounds[2];
    1280113152    bounds[0]=lower;
     
    1280713158    --nNode;
    1280813159    walkback_[nNode]->applyBounds(iColumn,lower,upper,force);
     13160#if 0
    1280913161    CbcNode * nodeLook = walkback_[nNode]->mutableOwner();
    12810     if (nodeLook&&0) {
     13162    if (nodeLook) {
    1281113163      const OsiBranchingObject * obj = nodeLook->branchingObject();
    1281213164      const CbcIntegerBranchingObject * objectI = dynamic_cast<const CbcIntegerBranchingObject *> (obj);
     
    1281713169      assert (iColumn!=iColumn2);
    1281813170    }
     13171#endif
    1281913172  }
    1282013173}
  • trunk/Cbc/src/CbcModel.hpp

    r854 r862  
    406406  */
    407407  void analyzeObjective();
     408  /// Tighten bounds - lightweight
     409  int tightenBounds();
    408410
    409411
  • trunk/Cbc/src/CbcSolver.cpp

    r857 r862  
    166166#include "OsiClpSolverInterface.hpp"
    167167#include "CbcSolver.hpp"
    168 #define IN_BRANCH_AND_BOUND (0x01000000|262144)
     168//#define IN_BRANCH_AND_BOUND (0x01000000|262144)
     169#define IN_BRANCH_AND_BOUND (0x01000000|262144|128|1024|2048)
     170//#define IN_BRANCH_AND_BOUND (0x01000000|262144|128)
    169171CbcSolver::CbcSolver()
    170172  : babModel_(NULL),
     
    35073509  double statistics_cut_time=0.0;
    35083510  int statistics_nodes=0, statistics_iterations=0;
     3511  int statistics_nrows=0, statistics_ncols=0;
     3512  int statistics_nprocessedrows=0, statistics_nprocessedcols=0;
    35093513  std::string statistics_result;
    35103514  memset(statusUserFunction_,0,numberUserFunctions_*sizeof(int));
     
    51555159              if (logLevel<=1)
    51565160                model_.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
     5161              else
     5162                model_.solver()->setHintParam(OsiDoReducePrint,false,OsiHintTry);
    51575163              {
    51585164                OsiSolverInterface * solver = model_.solver();
     
    51615167                assert (si != NULL);
    51625168                si->getModelPtr()->scaling(doScaling);
     5169                statistics_nrows=si->getNumRows();
     5170                statistics_ncols=si->getNumCols();
     5171                statistics_nprocessedrows=si->getNumRows();
     5172                statistics_nprocessedcols=si->getNumCols();
    51635173                // See if quadratic
    51645174#ifdef COIN_HAS_LINK
     
    57595769                  babModel_->setSecondaryStatus(1);
    57605770                } else {
     5771                  statistics_nprocessedrows=solver2->getNumRows();
     5772                  statistics_nprocessedcols=solver2->getNumCols();
    57615773                  model_.setProblemStatus(-1);
    57625774                  babModel_->setProblemStatus(-1);
     
    88508862                // first header if needed
    88518863                if (state!=2)
    8852                   fputs("Name,result,time,objective,continuous,tightened,cut_time,nodes,iterations\n",fp);
     8864                  fputs("Name,result,time,objective,continuous,tightened,cut_time,nodes,iterations,rows,columns,processed_rows,processed_columns\n",fp);
    88538865                strcpy(buffer,argv[1]);
    88548866                char * slash=buffer;
     
    88578869                    slash=buffer+i+1;
    88588870                }
    8859                 fprintf(fp,"%s,%s,%.2f,%.16g,%g,%g,%.2f,%d,%d\n",
     8871                fprintf(fp,"%s,%s,%.2f,%.16g,%g,%g,%.2f,%d,%d,%d,%d,%d,%d\n",
    88608872                        slash,statistics_result.c_str(),statistics_seconds,statistics_obj,
    88618873                        statistics_continuous,statistics_tighter,statistics_cut_time,statistics_nodes,
    8862                         statistics_iterations);
     8874                        statistics_iterations,statistics_nrows,statistics_ncols,
     8875                        statistics_nprocessedrows,statistics_nprocessedcols);
    88638876                fclose(fp);
    88648877              } else {
Note: See TracChangeset for help on using the changeset viewer.