Ignore:
Timestamp:
Jun 14, 2016 10:39:54 AM (6 years ago)
Author:
forrest
Message:

allow heuristics to see if integers are 'optional'

File:
1 edited

Legend:

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

    r2146 r2280  
    216216// Returns 1 if solution, 0 if not
    217217int
    218 CbcHeuristicFPump::solution(double & solutionValue,
     218CbcHeuristicFPump::solutionInternal(double & solutionValue,
    219219                            double * betterSolution)
    220220{
     
    286286    double cutoff;
    287287    model_->solver()->getDblParam(OsiDualObjectiveLimit, cutoff);
     288    double realCutoff=cutoff;
     289    bool secondMajorPass=false;
    288290    double direction = model_->solver()->getObjSense();
    289291    cutoff *= direction;
     
    336338        assert(integerObject || integerObject2);
    337339#endif
     340#ifdef COIN_HAS_CLP
     341        if (!isHeuristicInteger(model_->solver(),iColumn))
     342          continue;
     343#endif
    338344        double value = initialSolution[iColumn];
    339345        double nearest = floor(value + 0.5);
     
    370376        for (i = 0; i < numberIntegers; i++) {
    371377            int iColumn = integerVariableOrig[i];
     378#ifdef COIN_HAS_CLP
     379            if (!isHeuristicInteger(model_->solver(),iColumn))
     380              continue;
     381#endif
    372382            if (upper[iColumn] - lower[iColumn] < 1.000001)
    373383                integerVariable[j++] = iColumn;
     
    476486    double artificialFactor = 0.00001;
    477487    // also try rounding!
    478     double * roundingSolution = new double[numberColumns];
    479     double roundingObjective = solutionValue;
     488    double * roundingSolution = new double[2*numberColumns];
     489    double roundingObjective = realCutoff;
    480490    CbcRounding roundingHeuristic(*model_);
    481491    int dualPass = 0;
     
    807817                randomNumberGenerator_.setSeed(987654321);
    808818            }
     819#ifdef COIN_HAS_CLP
     820            {
     821              OsiClpSolverInterface * clpSolver
     822                = dynamic_cast<OsiClpSolverInterface *> (clonedSolver);
     823              //printf("real cutoff %g fake %g - second pass %c\n",realCutoff,cutoff,
     824              //     secondMajorPass ? 'Y' : 'N');
     825              if (clpSolver && (((accumulate_&16) != 0) ||
     826                                ((accumulate_&8)!=0 && secondMajorPass))) {
     827                // try rounding heuristic
     828                OsiSolverInterface * saveSolver = model_->swapSolver(clonedSolver);
     829                ClpSimplex * simplex = clpSolver->getModelPtr();
     830                double * solverSolution = simplex->primalColumnSolution();
     831                memcpy(solverSolution,solution,numberColumns*sizeof(double));
     832                // Compute using dot product
     833                double newSolutionValue = -saveOffset;
     834                for (  i = 0 ; i < numberColumns ; i++ )
     835                  newSolutionValue += saveObjective[i] * solution[i];
     836                simplex->setObjectiveValue(newSolutionValue);
     837                clpSolver->setObjective(saveObjective);
     838                CbcRounding heuristic1(*model_);
     839                heuristic1.setHeuristicName("rounding in feaspump!");
     840                heuristic1.setWhen(1);
     841                newSolutionValue=realCutoff;
     842                int ifSolR = heuristic1.solution(newSolutionValue,
     843                                                 roundingSolution+numberColumns);
     844                model_->swapSolver(saveSolver);
     845                if (ifSolR&&newSolutionValue<roundingObjective) {
     846                  roundingObjective=newSolutionValue;
     847                  //printf("rounding obj of %g?\n", roundingObjective);
     848                  memcpy(roundingSolution,roundingSolution+numberColumns,
     849                         numberColumns*sizeof(double));
     850                }
     851              }
     852            }
     853#endif
    809854            returnCode = rounds(solver, newSolution,/*saveObjective,*/
    810855                                numberIntegers, integerVariable,
     
    817862            }
    818863            numberPasses++;
     864            if (roundingObjective<realCutoff) {
     865              if (returnCode) {
     866                newSolutionValue = -saveOffset;
     867                for (  i = 0 ; i < numberColumns ; i++ )
     868                    newSolutionValue += saveObjective[i] * newSolution[i];
     869              } else {
     870                newSolutionValue = COIN_DBL_MAX;
     871              }
     872              if (roundingObjective<newSolutionValue&&false) {
     873                returnCode=1;
     874                memcpy(newSolution,roundingSolution,
     875                       numberColumns*sizeof(double));
     876              }
     877            }
    819878            if (returnCode) {
    820879                // SOLUTION IS INTEGER
     
    841900                        for (i = 0; i < numberIntegersOrig; i++) {
    842901                            int iColumn = integerVariableOrig[i];
     902#ifdef COIN_HAS_CLP
     903                            if (!isHeuristicInteger(solver,iColumn))
     904                              continue;
     905#endif
    843906                            double value = floor(newSolution[iColumn] + 0.5);
    844907                            if (solver->isBinary(iColumn)) {
     
    932995                                }
    933996                                newSolutionValue += objective[iColumn] * newSolution[iColumn];
    934                                 if (solver->isInteger(iColumn)) {
     997                                if (isHeuristicInteger(solver,iColumn)) {
    935998                                  double intValue = floor(value+0.5);
    936999                                  if (fabs(value-intValue)>largestAway) {
     
    9931056                                    model_->setBestSolution(saveOldSolution, numberColumns, saveObjectiveValue);
    9941057                                delete [] saveOldSolution;
     1058                                realCutoff = model_->getMinimizationObjValue();
    9951059                            }
    9961060                            if (action == 0 || model_->maximumSecondsReached()) {
     
    11561220                        continue;
    11571221                    // deal with fixed variables (i.e., upper=lower)
    1158                     if (fabs(lower[iColumn] - upper[iColumn]) < primalTolerance || !solver->isInteger(iColumn)) {
     1222                    if (fabs(lower[iColumn] - upper[iColumn]) < primalTolerance || !isHeuristicInteger(solver,iColumn)) {
    11591223                        //if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue);
    11601224                        //else                                       solver->setObjCoeff(iColumn,costValue);
     
    13181382                                for (i = 0; i < numberIntegersOrig; i++) {
    13191383                                    int iColumn = integerVariableOrig[i];
     1384#ifdef COIN_HAS_CLP
     1385                                    if (!isHeuristicInteger(solver,iColumn))
     1386                                      continue;
     1387#endif
    13201388#ifdef RAND_RAND
    13211389                                    if (nChanged > numberChanged)
     
    13411409                                    for (i = 0; i < numberIntegersOrig; i++) {
    13421410                                        int iColumn = integerVariableOrig[i];
     1411#ifdef COIN_HAS_CLP
     1412                                        if (!isHeuristicInteger(solver,iColumn))
     1413                                          continue;
     1414#endif
    13431415                                        if (objective[iColumn] > 0.0) {
    13441416                                            if (basis->getStructStatus(iColumn) ==
     
    13951467                        const double * newSolution = solver->getColSolution();
    13961468                        for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++ ) {
    1397                             if (solver->isInteger(iColumn)) {
     1469                            if (isHeuristicInteger(solver,iColumn)) {
    13981470                                double value = newSolution[iColumn];
    13991471                                double nearest = floor(value + 0.5);
     
    14591531                            for (  i = 0 ; i < numberColumns ; i++ ) {
    14601532                                double value = newSolution[i];
    1461                                 if (solver->isInteger(i)) {
     1533                                if (isHeuristicInteger(solver,i)) {
    14621534                                    double nearest = floor(value + 0.5);
    14631535                                    if (fabs(value - nearest) > 1.0e-6) {
     
    16721744                        const double * newSolution = solver->getColSolution();
    16731745                        for (  i = 0 ; i < numberColumns ; i++ ) {
    1674                             if (solver->isInteger(i)) {
     1746                            if (isHeuristicInteger(solver,i)) {
    16751747                                double value = newSolution[i];
    16761748                                double nearest = floor(value + 0.5);
     
    17401812                    for (i = 0; i < numberIntegersOrig; i++) {
    17411813                        int iColumn = integerVariableOrig[i];
     1814#ifdef COIN_HAS_CLP
     1815                        if (!isHeuristicInteger(solver,iColumn))
     1816                          continue;
     1817#endif
    17421818                        if (objective[iColumn] > 0.0)
    17431819                            closestSolution[i] = 0;
     
    18571933            solutionValue = roundingObjective;
    18581934            newSolutionValue = solutionValue;
     1935            realCutoff=solutionValue-1.0e-5;
    18591936            memcpy(betterSolution, roundingSolution, numberColumns*sizeof(double));
    18601937            solutionFound = true;
     
    19031980                for (i = 0; i < numberIntegersOrig; i++) {
    19041981                    int iColumn = integerVariableOrig[i];
     1982#ifdef COIN_HAS_CLP
     1983                    if (!isHeuristicInteger(newSolver,iColumn))
     1984                      continue;
     1985#endif
    19051986                    //const OsiObject * object = model_->object(i);
    19061987                    //double originalLower;
     
    19342015                if (fixContinuous) {
    19352016                    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    1936                         if (!newSolver->isInteger(iColumn) && usedColumn[iColumn] <= allowedPass) {
     2017                        if (!isHeuristicInteger(newSolver,iColumn) && usedColumn[iColumn] <= allowedPass) {
    19372018                            double value = lastSolution[iColumn];
    19382019                            if (value < colLower[iColumn] + 1.0e-8) {
     
    20422123                    const double * obj = newSolver->getObjCoefficients();
    20432124                    for (int i = 0 ; i < numberColumns ; i++ ) {
    2044                         if (newSolver->isInteger(i)) {
     2125                        if (isHeuristicInteger(newSolver,i)) {
    20452126                            double value = newSolution[i];
    20462127                            double nearest = floor(value + 0.5);
     
    20672148                        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    20682149                            double value = newSolution[iColumn];
    2069                             if (newSolver->isInteger(iColumn)) {
     2150                            if (isHeuristicInteger(newSolver,iColumn)) {
    20702151                                value = floor(newSolution[iColumn] + 0.5);
    20712152                                newSolver->setColLower(iColumn, value);
     
    21502231                                    const double * solution = newSolver->getColSolution();
    21512232                                    for (int  i = 0 ; i < numberColumns ; i++ ) {
    2152                                         if (newSolver->isInteger(i)) {
     2233                                        if (isHeuristicInteger(newSolver,i)) {
    21532234                                            double value = solution[i];
    21542235                                            double nearest = floor(value + 0.5);
     
    21922273                        const double * upper = model_->solver()->getColUpper();
    21932274                        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    2194                             if (newSolver->isInteger(iColumn)) {
     2275                            if (isHeuristicInteger(newSolver,iColumn)) {
    21952276                                double value = floor(newSolution[iColumn] + 0.5);
    21962277                                newSolver->setColLower(iColumn, value);
     
    22532334        if (solutionFound) finalReturnCode = 1;
    22542335        cutoff = CoinMin(cutoff, solutionValue - model_->getCutoffIncrement());
     2336        realCutoff = cutoff;
    22552337        if (numberTries >= maximumRetries_ || !solutionFound || exitAll || cutoff < continuousObjectiveValue + 1.0e-7) {
    22562338            break;
     
    22722354                break;
    22732355            sprintf(pumpPrint, "Round again with cutoff of %g", cutoff);
     2356            secondMajorPass=true;
    22742357            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
    22752358            << pumpPrint
     
    23012384        for (i = 0; i < numberIntegersOrig; i++) {
    23022385            int iColumn = integerVariableOrig[i];
     2386#ifdef COIN_HAS_CLP
     2387            if (!isHeuristicInteger(newSolver,iColumn))
     2388              continue;
     2389#endif
    23032390            int direction = closestSolution[i];
    23042391            closestSolution[i] = iColumn;
     
    23832470
    23842471/**************************END MAIN PROCEDURE ***********************************/
     2472/* If general integers then adds variables to turn into binaries round
     2473   solution
     2474*/
     2475int
     2476CbcHeuristicFPump::solution(double & objectiveValue, double * newSolution)
     2477{
     2478  double * newSolution2=NULL;
     2479  double objective2=COIN_DBL_MAX;
     2480  int returnCode2=0;
     2481  int oddGeneral = (accumulate_&(32|64|128))>>5;
     2482  if (oddGeneral) {
     2483    int maxAround=2;
     2484    bool fixSatisfied=false;
     2485    // clone solver and modify
     2486    OsiSolverInterface * solver = cloneBut(2); // wasmodel_->solver()->clone();
     2487    double cutoff;
     2488    model_->solver()->getDblParam(OsiDualObjectiveLimit, cutoff);
     2489    int numberColumns = model_->getNumCols();
     2490    int numberIntegers = model_->numberIntegers();
     2491    const int * integerVariableOrig = model_->integerVariable();
     2492    int general = 0;
     2493    int nAdd=0;
     2494    const double * lower = solver->getColLower();
     2495    const double * upper = solver->getColUpper();
     2496    const double * initialSolution = solver->getColSolution();
     2497    // we may be being clever so make sure solver lines up wuth model
     2498    for (int i = 0; i < numberColumns; i++)
     2499      solver->setContinuous(i);
     2500    for (int i = 0; i < numberIntegers; i++) {
     2501      int iColumn = integerVariableOrig[i];
     2502#ifdef COIN_HAS_CLP
     2503      if (!isHeuristicInteger(solver,iColumn))
     2504        continue;
     2505#endif
     2506      double value = initialSolution[iColumn];
     2507      double nearest = floor(value + 0.5);
     2508      if (upper[iColumn] - lower[iColumn] > 1.000001) {
     2509        if (fabs(value-nearest)<1.0e-6&&fixSatisfied) {
     2510          solver->setColLower(iColumn,nearest);
     2511          solver->setColUpper(iColumn,nearest);
     2512        } else {
     2513          general++;
     2514          int up=static_cast<int>(upper[iColumn]);
     2515          int lo=static_cast<int>(lower[iColumn]);
     2516          int near=static_cast<int>(nearest);
     2517          up = CoinMin(up,near+maxAround);
     2518          lo = CoinMax(lo,near-maxAround);
     2519          solver->setColLower(iColumn,lo);
     2520          solver->setColUpper(iColumn,up);
     2521          int n=up-lo;
     2522          // 1 - 1, 2,3 - 2, 4567 - 3
     2523          while (n) {
     2524            nAdd++;
     2525            n = n>>1;
     2526          }
     2527        }
     2528      } else {
     2529        solver->setInteger(iColumn);
     2530      }
     2531    }
     2532    if (!general) {
     2533      delete solver;
     2534    }
     2535    if (general) {
     2536      CbcModel * saveModel = model_;
     2537      int * addStart = new int[nAdd+1];
     2538      memset(addStart,0,(nAdd+1)*sizeof(int));
     2539      int * addIndex = new int[general+nAdd];
     2540      double * addElement = new double[general+nAdd];
     2541      double * addLower = new double[nAdd];
     2542      double * addUpper = new double[nAdd];
     2543      for (int i=0;i<nAdd;i++) {
     2544        addLower[i]=0.0;
     2545        addUpper[i]=1.0;
     2546      }
     2547      solver->addCols(nAdd, addStart, NULL, NULL, addLower, addUpper, NULL);
     2548      lower = solver->getColLower();
     2549      upper = solver->getColUpper();
     2550      // now rows
     2551      nAdd = 0;
     2552      int nEl = 0;
     2553      int nAddRow = 0;
     2554      for (int i = 0; i < numberIntegers; i++) {
     2555        int iColumn = integerVariableOrig[i];
     2556#ifdef COIN_HAS_CLP
     2557        if (!isHeuristicInteger(solver,iColumn))
     2558          continue;
     2559#endif
     2560        if (upper[iColumn] - lower[iColumn] > 1.000001) {
     2561          int up=static_cast<int>(upper[iColumn]);
     2562          int lo=static_cast<int>(lower[iColumn]);
     2563          addLower[nAddRow] = lo;
     2564          addUpper[nAddRow] = lo;
     2565          addIndex[nEl] = iColumn;
     2566          addElement[nEl++] = 1.0;
     2567          int n=up-lo;
     2568          // 1 - 1, 2,3 - 2, 4567 - 3
     2569          int el=1;
     2570          while (n) {
     2571            addIndex[nEl] = numberColumns+nAdd;
     2572            addElement[nEl++] = -el;
     2573            nAdd++;
     2574            n = n>>1;
     2575            el = el<<1;
     2576          }
     2577          nAddRow++;
     2578          addStart[nAddRow]=nEl;
     2579        }
     2580      }
     2581      for (int i=0;i<nAdd;i++)
     2582        solver->setInteger(i+numberColumns);
     2583      solver->addRows(nAddRow, addStart, addIndex, addElement, addLower, addUpper);
     2584      delete [] addStart;
     2585      delete [] addIndex;
     2586      delete [] addElement;
     2587      delete [] addLower;
     2588      delete [] addUpper;
     2589      solver->resolve();
     2590      solver->writeMps("test");
     2591      // new CbcModel
     2592      model_=new CbcModel(*solver);
     2593      model_->findIntegers(true);
     2594      // set cutoff
     2595      solver->setDblParam(OsiDualObjectiveLimit, cutoff);
     2596      model_->setCutoff(cutoff);
     2597      newSolution2 = new double [numberColumns+nAdd];
     2598      objective2=objectiveValue;
     2599      returnCode2=solutionInternal(objective2,newSolution2);
     2600      delete solver;
     2601      delete model_;
     2602      model_=saveModel;
     2603    }
     2604  }
     2605  int returnCode=solutionInternal(objectiveValue,newSolution);
     2606  if (returnCode2&&false) {
     2607    int numberColumns=model_->getNumCols();
     2608    memcpy(newSolution,newSolution2,numberColumns*sizeof(double));
     2609  }
     2610  delete [] newSolution2;
     2611  return returnCode;
     2612}
    23852613
    23862614// update model
     
    24332661                for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
    24342662                    int iColumn = column[k];
    2435                     if (elementByRow[k] != 1.0 || !solver->isInteger(iColumn)) {
     2663                    if (elementByRow[k] != 1.0 || !isHeuristicInteger(solver,iColumn)) {
    24362664                        cover = false;
    24372665                        break;
     
    24792707            for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
    24802708                int iColumn = column[k];
    2481                 if (elementByRow[k] != 1.0 || !solver->isInteger(iColumn)) {
     2709                if (elementByRow[k] != 1.0 || !isHeuristicInteger(solver,iColumn)) {
    24822710                    cover = false;
    24832711                    break;
     
    27002928            // First improve by moving continuous ones
    27012929            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    2702                 if (!solver->isInteger(iColumn)) {
     2930                if (!isHeuristicInteger(solver,iColumn)) {
    27032931                    double solValue = solution[iColumn];
    27042932                    double thetaUp = columnUpper[iColumn] - solValue;
Note: See TracChangeset for help on using the changeset viewer.