Changeset 759


Ignore:
Timestamp:
Aug 17, 2007 10:15:27 AM (12 years ago)
Author:
forrest
Message:

start of work on new vub heuristic

Location:
trunk/Cbc/src
Files:
7 edited

Legend:

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

    r747 r759  
    141141    CglPreProcess process;
    142142    /* Do not try and produce equality cliques and
    143        do up to 5 passes */
     143       do up to 2 passes */
    144144    if (logLevel<=1)
    145145      process.messageHandler()->setLogLevel(0);
    146     OsiSolverInterface * solver2= process.preProcess(*solver);
     146    OsiSolverInterface * solver2= process.preProcessNonDefault(*solver,false,2);
    147147    if (!solver2) {
    148148      if (logLevel>1)
     
    153153      double before = solver->getNumRows()+solver->getNumCols();
    154154      double after = solver2->getNumRows()+solver2->getNumCols();
    155       if (logLevel>1)
    156         printf("before %d rows %d columns, after %d rows %d columns\n",
    157                solver->getNumRows(),solver->getNumCols(),
    158                solver2->getNumRows(),solver2->getNumCols());
    159       if (after>fractionSmall_*before&&after>300)
    160         return 0;
     155      char generalPrint[100];
     156      if (after>fractionSmall_*before&&after>300) {
     157        sprintf(generalPrint,"Full problem %d rows %d columns, reduced to %d rows %d columns - too large",
     158                solver->getNumRows(),solver->getNumCols(),
     159                solver2->getNumRows(),solver2->getNumCols());
     160        model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
     161          << generalPrint
     162          <<CoinMessageEol;
     163        return -1;
     164      } else {
     165        sprintf(generalPrint,"Full problem %d rows %d columns, reduced to %d rows %d columns",
     166                solver->getNumRows(),solver->getNumCols(),
     167                solver2->getNumRows(),solver2->getNumCols());
     168        model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
     169          << generalPrint
     170          <<CoinMessageEol;
     171      }
    161172      solver2->resolve();
    162173      CbcModel model(*solver2);
  • trunk/Cbc/src/CbcHeuristicFPump.cpp

    r758 r759  
    478478              returnCode = smallBranchAndBound(solver,numberNodes_,newSolution,newSolutionValue,
    479479                                               solutionValue,"CbcHeuristicFpump");
     480              if (returnCode<0)
     481                returnCode=0; // returned on size
    480482              if ((returnCode&2)!=0) {
    481483                // could add cut
     
    979981      returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,newSolutionValue,
    980982                                       cutoff,"CbcHeuristicLocalAfterFPump");
     983      if (returnCode<0)
     984        returnCode=0; // returned on size - could try changing
    981985      if ((returnCode&2)!=0) {
    982986        // could add cut
     
    11071111    int returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,newSolutionValue,
    11081112                                     newSolutionValue,"CbcHeuristicLocalAfterFPump");
     1113    if (returnCode<0)
     1114      returnCode=0; // returned on size
    11091115    if ((returnCode&2)!=0) {
    11101116      // could add cut
  • trunk/Cbc/src/CbcHeuristicGreedy.cpp

    r687 r759  
    731731    int returnCode = smallBranchAndBound(newSolver,200,newSolution,newSolutionValue,
    732732                                         solutionValue,"CbcHeuristicGreedy");
     733    if (returnCode<0)
     734      returnCode=0; // returned on size
    733735    if ((returnCode&2)!=0) {
    734736      // could add cut
  • trunk/Cbc/src/CbcHeuristicLocal.cpp

    r738 r759  
    156156  int returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,objectiveValue,
    157157                                         objectiveValue,"CbcHeuristicLocal");
     158  if (returnCode<0)
     159    returnCode=0; // returned on size
    158160  if ((returnCode&2)!=0) {
    159161    // could add cut
  • trunk/Cbc/src/CbcHeuristicRINS.cpp

    r640 r759  
    203203      returnCode = smallBranchAndBound(newSolver,numberNodes_,betterSolution,solutionValue,
    204204                                         model_->getCutoff(),"CbcHeuristicRINS");
     205      if (returnCode<0)
     206        returnCode=0; // returned on size
    205207      if ((returnCode&1)!=0)
    206208        numberSuccesses_++;
  • trunk/Cbc/src/CbcMain.cpp

    r388 r759  
    836836            }
    837837          } else if (valid==1) {
    838             abort();
     838            std::cout<<" is illegal for double parameter "<<parameters[iParam].name()<<" value remains "<<
     839              parameters[iParam].doubleValue()<<std::endl;
    839840          } else {
    840841            switch(type) {
     
    860861            parameters[iParam].setIntParameter(*model,value);
    861862          } else if (valid==1) {
    862             abort();
     863            std::cout<<" is illegal for integer parameter "<<parameters[iParam].name()<<" value remains "<<
     864              parameters[iParam].intValue()<<std::endl;
    863865          } else {
    864866            std::cout<<parameters[iParam].name()<<" has value "<<
  • trunk/Cbc/src/CbcSolver.cpp

    r758 r759  
    487487  }
    488488}
     489#if 1
     490/*
     491  On input
     492  doAction - 0 just fix in original and return NULL
     493             1 return fixed non-presolved solver
     494             2 return fixed presolved solver
     495             3 do heuristics and set best solution
     496             4 do BAB and just set best solution
     497             -2 cleanup afterwards if using 2
     498  On ouput - number fixed
     499*/
     500static OsiClpSolverInterface * fixVubs(CbcModel & model, int numberLevels,
     501                                       int & doAction, CoinMessageHandler * generalMessageHandler)
     502{
     503  OsiSolverInterface * solver = model.solver()->clone();
     504  OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
     505  ClpSimplex * lpSolver = clpSolver->getModelPtr();
     506  // Tighten bounds
     507  lpSolver->tightenPrimalBounds(0.0,11,true);
     508  char generalPrint[200];
     509  //const double *objective = solver->getObjCoefficients() ;
     510  double *columnLower = lpSolver->columnLower() ;
     511  double *columnUpper = lpSolver->columnUpper() ;
     512  int numberColumns = solver->getNumCols() ;
     513  int numberRows = solver->getNumRows();
     514  int iRow,iColumn;
     515
     516  // Row copy
     517  CoinPackedMatrix matrixByRow(*solver->getMatrixByRow());
     518  const double * elementByRow = matrixByRow.getElements();
     519  const int * column = matrixByRow.getIndices();
     520  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
     521  const int * rowLength = matrixByRow.getVectorLengths();
     522
     523  // Column copy
     524  CoinPackedMatrix  matrixByCol(*solver->getMatrixByCol());
     525  //const double * element = matrixByCol.getElements();
     526  const int * row = matrixByCol.getIndices();
     527  const CoinBigIndex * columnStart = matrixByCol.getVectorStarts();
     528  const int * columnLength = matrixByCol.getVectorLengths();
     529
     530  const double * rowLower = solver->getRowLower();
     531  const double * rowUpper = solver->getRowUpper();
     532
     533  // Get maximum size of VUB tree
     534  // otherColumn is one fixed to 0 if this one zero
     535  int nEl = matrixByCol.getNumElements();
     536  CoinBigIndex * fixColumn = new CoinBigIndex [numberColumns+1];
     537  int * otherColumn = new int [nEl];
     538  int * fix = new int[numberColumns];
     539  char * mark = new char [numberColumns];
     540  memset(mark,0,numberColumns);
     541  int numberInteger=0;
     542  int numberOther=0;
     543  fixColumn[0]=0;
     544  double large = lpSolver->largeValue(); // treat bounds > this as infinite
     545#ifndef NDEBUG
     546  double large2= 1.0e10*large;
     547#endif
     548  double tolerance = lpSolver->primalTolerance();
     549  int * check = new int[numberRows];
     550  for (iRow=0;iRow<numberRows;iRow++)
     551    check[iRow]=-1;
     552  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     553    fix[iColumn]=-1;
     554    if (columnUpper[iColumn] > columnLower[iColumn]+1.0e-8) {
     555      if (solver->isInteger(iColumn))
     556        numberInteger++;
     557      if (columnLower[iColumn]==0.0) {
     558        bool infeasible=false;
     559        fix[iColumn]=0;
     560        // fake upper bound
     561        double saveUpper = columnUpper[iColumn];
     562        columnUpper[iColumn]=0.0;
     563        for (CoinBigIndex i=columnStart[iColumn];
     564             i<columnStart[iColumn]+columnLength[iColumn];i++) {
     565          iRow = row[i];
     566          if (rowLower[iRow]<-1.0e6&&rowUpper[iRow]>1.0e6)
     567            continue; // unlikely
     568          //==
     569          // possible row
     570          int infiniteUpper = 0;
     571          int infiniteLower = 0;
     572          double maximumUp = 0.0;
     573          double maximumDown = 0.0;
     574          double newBound;
     575          CoinBigIndex rStart = rowStart[iRow];
     576          CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
     577          CoinBigIndex j;
     578          int kColumn;
     579          // Compute possible lower and upper ranges
     580          for (j = rStart; j < rEnd; ++j) {
     581            double value=elementByRow[j];
     582            kColumn = column[j];
     583            if (value > 0.0) {
     584              if (columnUpper[kColumn] >= large) {
     585                ++infiniteUpper;
     586              } else {
     587                maximumUp += columnUpper[kColumn] * value;
     588              }
     589              if (columnLower[kColumn] <= -large) {
     590                ++infiniteLower;
     591              } else {
     592                maximumDown += columnLower[kColumn] * value;
     593              }
     594            } else if (value<0.0) {
     595              if (columnUpper[kColumn] >= large) {
     596                ++infiniteLower;
     597              } else {
     598                maximumDown += columnUpper[kColumn] * value;
     599              }
     600              if (columnLower[kColumn] <= -large) {
     601                ++infiniteUpper;
     602              } else {
     603                maximumUp += columnLower[kColumn] * value;
     604              }
     605            }
     606          }
     607          // Build in a margin of error
     608          maximumUp += 1.0e-8*fabs(maximumUp);
     609          maximumDown -= 1.0e-8*fabs(maximumDown);
     610          double maxUp = maximumUp+infiniteUpper*1.0e31;
     611          double maxDown = maximumDown-infiniteLower*1.0e31;
     612          if (maxUp <= rowUpper[iRow] + tolerance &&
     613              maxDown >= rowLower[iRow] - tolerance) {
     614            //printf("Redundant row in vubs %d\n",iRow);
     615          } else {
     616            if (maxUp < rowLower[iRow] -100.0*tolerance ||
     617                maxDown > rowUpper[iRow]+100.0*tolerance) {
     618              infeasible=true;
     619              break;
     620            }
     621            double lower = rowLower[iRow];
     622            double upper = rowUpper[iRow];
     623            for (j = rStart; j < rEnd; ++j) {
     624              double value=elementByRow[j];
     625              kColumn = column[j];
     626              double nowLower = columnLower[kColumn];
     627              double nowUpper = columnUpper[kColumn];
     628              if (value > 0.0) {
     629                // positive value
     630                if (lower>-large) {
     631                  if (!infiniteUpper) {
     632                    assert(nowUpper < large2);
     633                    newBound = nowUpper +
     634                      (lower - maximumUp) / value;
     635                    // relax if original was large
     636                    if (fabs(maximumUp)>1.0e8)
     637                      newBound -= 1.0e-12*fabs(maximumUp);
     638                  } else if (infiniteUpper==1&&nowUpper>large) {
     639                    newBound = (lower -maximumUp) / value;
     640                    // relax if original was large
     641                    if (fabs(maximumUp)>1.0e8)
     642                      newBound -= 1.0e-12*fabs(maximumUp);
     643                  } else {
     644                    newBound = -COIN_DBL_MAX;
     645                  }
     646                  if (newBound > nowLower + 1.0e-12&&newBound>-large) {
     647                    // Tighten the lower bound
     648                    // check infeasible (relaxed)
     649                    if (nowUpper < newBound) {
     650                      if (nowUpper - newBound <
     651                          -100.0*tolerance) {
     652                        infeasible=true;
     653                        break;
     654                      }
     655                    }
     656                  }
     657                }
     658                if (upper <large) {
     659                  if (!infiniteLower) {
     660                    assert(nowLower >- large2);
     661                    newBound = nowLower +
     662                      (upper - maximumDown) / value;
     663                    // relax if original was large
     664                    if (fabs(maximumDown)>1.0e8)
     665                    newBound += 1.0e-12*fabs(maximumDown);
     666                  } else if (infiniteLower==1&&nowLower<-large) {
     667                    newBound =   (upper - maximumDown) / value;
     668                    // relax if original was large
     669                    if (fabs(maximumDown)>1.0e8)
     670                      newBound += 1.0e-12*fabs(maximumDown);
     671                  } else {
     672                    newBound = COIN_DBL_MAX;
     673                  }
     674                  if (newBound < nowUpper - 1.0e-12&&newBound<large) {
     675                    // Tighten the upper bound
     676                    // check infeasible (relaxed)
     677                    if (nowLower > newBound) {
     678                      if (newBound - nowLower <
     679                          -100.0*tolerance) {
     680                        infeasible=true;
     681                        break;
     682                      } else {
     683                        newBound=nowLower;
     684                      }
     685                    }
     686                    if (!newBound||(solver->isInteger(kColumn)&&newBound<0.999)) {
     687                      // fix to zero
     688                      if (!mark[kColumn]) {
     689                        otherColumn[numberOther++]=kColumn;
     690                        mark[kColumn]=1;
     691                        if (check[iRow]==-1)
     692                          check[iRow]=iColumn;
     693                        else
     694                          assert(check[iRow]==iColumn);
     695                      }
     696                    }
     697                  }
     698                }
     699              } else {
     700                // negative value
     701                if (lower>-large) {
     702                  if (!infiniteUpper) {
     703                    assert(nowLower < large2);
     704                    newBound = nowLower +
     705                      (lower - maximumUp) / value;
     706                    // relax if original was large
     707                    if (fabs(maximumUp)>1.0e8)
     708                      newBound += 1.0e-12*fabs(maximumUp);
     709                  } else if (infiniteUpper==1&&nowLower<-large) {
     710                    newBound = (lower -maximumUp) / value;
     711                    // relax if original was large
     712                    if (fabs(maximumUp)>1.0e8)
     713                      newBound += 1.0e-12*fabs(maximumUp);
     714                  } else {
     715                    newBound = COIN_DBL_MAX;
     716                  }
     717                  if (newBound < nowUpper - 1.0e-12&&newBound<large) {
     718                    // Tighten the upper bound
     719                    // check infeasible (relaxed)
     720                    if (nowLower > newBound) {
     721                      if (newBound - nowLower <
     722                          -100.0*tolerance) {
     723                        infeasible=true;
     724                        break;
     725                      } else {
     726                        newBound=nowLower;
     727                      }
     728                    }
     729                    if (!newBound||(solver->isInteger(kColumn)&&newBound<0.999)) {
     730                      // fix to zero
     731                      if (!mark[kColumn]) {
     732                        otherColumn[numberOther++]=kColumn;
     733                        mark[kColumn]=1;
     734                        if (check[iRow]==-1)
     735                          check[iRow]=iColumn;
     736                        else
     737                          assert(check[iRow]==iColumn);
     738                      }
     739                    }
     740                  }
     741                }
     742                if (upper <large) {
     743                  if (!infiniteLower) {
     744                    assert(nowUpper < large2);
     745                    newBound = nowUpper +
     746                      (upper - maximumDown) / value;
     747                    // relax if original was large
     748                    if (fabs(maximumDown)>1.0e8)
     749                      newBound -= 1.0e-12*fabs(maximumDown);
     750                  } else if (infiniteLower==1&&nowUpper>large) {
     751                    newBound =   (upper - maximumDown) / value;
     752                    // relax if original was large
     753                    if (fabs(maximumDown)>1.0e8)
     754                      newBound -= 1.0e-12*fabs(maximumDown);
     755                  } else {
     756                    newBound = -COIN_DBL_MAX;
     757                  }
     758                  if (newBound > nowLower + 1.0e-12&&newBound>-large) {
     759                    // Tighten the lower bound
     760                    // check infeasible (relaxed)
     761                    if (nowUpper < newBound) {
     762                      if (nowUpper - newBound <
     763                          -100.0*tolerance) {
     764                        infeasible=true;
     765                        break;
     766                      }
     767                    }
     768                  }
     769                }
     770              }
     771            }
     772          }
     773        }
     774        for (int i=fixColumn[iColumn];i<numberOther;i++)
     775          mark[otherColumn[i]]=0;
     776        // reset bound unless infeasible
     777        if (!infeasible||!solver->isInteger(iColumn))
     778          columnUpper[iColumn]=saveUpper;
     779        else if (solver->isInteger(iColumn))
     780          columnLower[iColumn]=1.0;
     781      }
     782    }
     783    fixColumn[iColumn+1]=numberOther;
     784  }
     785  delete [] check;
     786  delete [] mark;
     787  // Now do reverse way
     788  int * counts = new int [numberColumns];
     789  CoinZeroN(counts,numberColumns);
     790  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     791    for (int i=fixColumn[iColumn];i<fixColumn[iColumn+1];i++)
     792      counts[otherColumn[i]]++;
     793  }
     794  numberOther=0;
     795  CoinBigIndex * fixColumn2 = new CoinBigIndex [numberColumns+1];
     796  int * otherColumn2 = new int [fixColumn[numberColumns]];
     797  fixColumn2[0]=0;
     798  for ( iColumn=0;iColumn<numberColumns;iColumn++) {
     799    numberOther += counts[iColumn];
     800    counts[iColumn]=0;
     801    fixColumn2[iColumn+1]=numberOther;
     802  }
     803  for ( iColumn=0;iColumn<numberColumns;iColumn++) {
     804    for (int i=fixColumn[iColumn];i<fixColumn[iColumn+1];i++) {
     805      int jColumn=otherColumn[i];
     806      CoinBigIndex put = fixColumn2[jColumn] + counts[jColumn];
     807      counts[jColumn]++;
     808      otherColumn2[put]=iColumn;
     809    }
     810  }
     811  // get top layer i.e. those which are not fixed by any other
     812  int kLayer=0;
     813  while (true) {
     814    int numberLayered=0;
     815    int numberInteger=0;
     816    for ( iColumn=0;iColumn<numberColumns;iColumn++) {
     817      if (fix[iColumn]==kLayer) {
     818        for (int i=fixColumn2[iColumn];i<fixColumn2[iColumn+1];i++) {
     819          int jColumn=otherColumn2[i];
     820          if (fix[jColumn]==kLayer) {
     821            fix[iColumn]=kLayer+1;
     822          }
     823        }
     824      }
     825      if (fix[iColumn]==kLayer) {
     826        numberLayered++;
     827        if (solver->isInteger(iColumn))
     828          numberInteger++;
     829      }
     830    }
     831    if (numberLayered) {
     832      printf("%d (%d integer) at priority %d\n",numberLayered,numberInteger,kLayer);
     833      kLayer++;
     834    } else {
     835      break;
     836    }
     837  }
     838  delete [] counts;
     839  delete [] fixColumn;
     840  delete [] otherColumn;
     841  delete [] otherColumn2;
     842  delete [] fixColumn2;
     843  // Now do fixing
     844  delete [] fix;
     845  delete solver;
     846  return NULL;
     847}
     848#endif
    489849#ifdef COIN_HAS_ASL
    490850/*  Returns OsiSolverInterface (User should delete)
     
    14531813    anyToDo=true;
    14541814    CbcHeuristicFPump heuristic4(*model);
    1455     heuristic4.setFractionSmall(0.6);
     1815    heuristic4.setFractionSmall(0.5);
     1816    double dextra3 = parameters[whichParam(DEXTRA3,numberParameters,parameters)].doubleValue();
     1817    if (dextra3)
     1818      heuristic4.setFractionSmall(dextra3);
    14561819    heuristic4.setMaximumPasses(parameters[whichParam(FPUMPITS,numberParameters,parameters)].intValue());
    14571820    int pumpTune=parameters[whichParam(FPUMPTUNE,numberParameters,parameters)].intValue();
     
    23002663            }
    23012664          } else if (valid==1) {
    2302             abort();
     2665            std::cout<<" is illegal for double parameter "<<parameters[iParam].name()<<" value remains "<<
     2666              parameters[iParam].doubleValue()<<std::endl;
    23032667          } else {
    23042668            std::cout<<parameters[iParam].name()<<" has value "<<
     
    23402704              else if (parameters[iParam].type()==CUTPASSINTREE)
    23412705                cutPassInTree = value;
    2342               else if (parameters[iParam].type()==FPUMPITS)
    2343                 { parameters[iParam].setIntValue(value);}
    23442706              else if (parameters[iParam].type()==STRONGBRANCHING||
    23452707                       parameters[iParam].type()==NUMBERBEFORE)
     
    23482710            }
    23492711          } else if (valid==1) {
    2350             abort();
     2712            std::cout<<" is illegal for integer parameter "<<parameters[iParam].name()<<" value remains "<<
     2713              parameters[iParam].intValue()<<std::endl;
    23512714          } else {
    23522715            std::cout<<parameters[iParam].name()<<" has value "<<
     
    30563419          case DOHEURISTIC:
    30573420            if (goodModel) {
     3421              int vubAction = parameters[whichParam(VUBTRY,numberParameters,parameters)].intValue();
     3422              if (vubAction!=-1) {
     3423                // look at vubs
     3424                OsiClpSolverInterface * newSolver = fixVubs(model,1,vubAction,generalMessageHandler);
     3425                assert (!newSolver);
     3426              }
    30583427              // Actually do heuristics
    30593428              doHeuristics(&model,2);
     
    49845353              // For best solution
    49855354              double * bestSolution = NULL;
     5355              printf ("bab best %g %x\n",babModel->getMinimizationObjValue(),babModel->bestSolution());
    49865356              if (babModel->getMinimizationObjValue()<1.0e50&&type==BAB) {
    49875357                // post process
     
    50205390                }
    50215391                checkSOS(babModel, babModel->solver());
     5392              } else if (model.bestSolution()&&type==BAB&&model.getMinimizationObjValue()&&preProcess) {
     5393                sprintf(generalPrint,"Restoring heuristic best solution of %g",model.getMinimizationObjValue());
     5394                generalMessageHandler->message(CLP_GENERAL,generalMessages)
     5395                  << generalPrint
     5396                  <<CoinMessageEol;
     5397                int n = saveSolver->getNumCols();
     5398                bestSolution = new double [n];
     5399                // Put solution now back in saveSolver
     5400                babModel->assignSolver(saveSolver);
     5401                saveSolver->setColSolution(model.bestSolution());
     5402                memcpy(bestSolution,babModel->solver()->getColSolution(),n*sizeof(double));
     5403                // and put back in very original solver
     5404                {
     5405                  ClpSimplex * original = originalSolver->getModelPtr();
     5406                  double * lower = original->columnLower();
     5407                  double * upper = original->columnUpper();
     5408                  double * solution = original->primalColumnSolution();
     5409                  int n = original->numberColumns();
     5410                  //assert (!n||n==babModel->solver()->getNumCols());
     5411                  for (int i=0;i<n;i++) {
     5412                    solution[i]=bestSolution[i];
     5413                    if (originalSolver->isInteger(i)) {
     5414                      lower[i]=solution[i];
     5415                      upper[i]=solution[i];
     5416                    }
     5417                  }
     5418                }
    50225419              }
    50235420              if (type==STRENGTHEN&&strengthenedModel)
Note: See TracChangeset for help on using the changeset viewer.