Ignore:
File:
1 edited

Legend:

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

    r1883 r2094  
    2626#include "CoinTime.hpp"
    2727#include "CbcEventHandler.hpp"
    28 
     28#ifdef SWITCH_VARIABLES
     29#include "CbcSimpleIntegerDynamicPseudoCost.hpp"
     30#endif
    2931
    3032// Default Constructor
     
    238240    if (!when() || (when() == 1 && model_->phase() != 1))
    239241        return 0; // switched off
     242#ifdef HEURISTIC_INFORM
     243    printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
     244           heuristicName(),numRuns_,numCouldRun_,when_);
     245#endif
    240246    // See if at root node
    241247    bool atRoot = model_->getNodeCount() == 0;
     
    271277            good = false;
    272278        }
    273         if (passNumber != 1 && !good)
     279        if (passNumber > 1 && !good)
    274280            return 0;
    275281    } else {
    276         if (passNumber != 1)
     282        if (passNumber > 1)
    277283            return 0;
    278284    }
     
    309315    const double * upper = model_->solver()->getColUpper();
    310316    bool doGeneral = (accumulate_ & 4) != 0;
     317    int numberUnsatisfied=0;
     318    double sumUnsatisfied=0.0;
     319    const double * initialSolution = model_->solver()->getColSolution();
    311320    j = 0;
    312321/*
     
    327336        assert(integerObject || integerObject2);
    328337#endif
     338        double value = initialSolution[iColumn];
     339        double nearest = floor(value + 0.5);
     340        sumUnsatisfied += fabs(value - nearest);
     341        if (fabs(value - nearest) > 1.0e-6)
     342          numberUnsatisfied++;
    329343        if (upper[iColumn] - lower[iColumn] > 1.000001) {
    330344            general++;
     
    366380        printf("DOing general with %d out of %d\n", general, numberIntegers);
    367381#endif
     382    sprintf(pumpPrint, "Initial state - %d integers unsatisfied sum - %g",
     383            numberUnsatisfied, sumUnsatisfied);
     384    model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
     385      << pumpPrint
     386      << CoinMessageEol;
    368387/*
    369388  This `closest solution' will satisfy integrality, but violate some other
     
    396415    model_->solver()->resolve();
    397416    if (!model_->solver()->isProvenOptimal()) {
    398         // presumably max time or some such
     417
     418        delete[] integerVariable;
     419        delete[] newSolution;
     420        if (closestSolution)
     421            delete[] closestSolution;
     422       
    399423        return 0;
    400424    }
     
    574598        // if cutoff exists then add constraint
    575599        bool useCutoff = (fabs(cutoff) < 1.0e20 && (fakeCutoff_ != COIN_DBL_MAX || numberTries > 1));
     600        bool tryOneClosePass=fakeCutoff_<solver->getObjValue();
    576601        // but there may be a close one
    577602        if (firstCutoff < 2.0*solutionValue && numberTries == 1 && CoinMin(cutoff, fakeCutoff_) < 1.0e20)
    578603            useCutoff = true;
    579         if (useCutoff) {
     604        if (useCutoff || tryOneClosePass) {
    580605            double rhs = CoinMin(cutoff, fakeCutoff_);
     606            if (tryOneClosePass) {
     607              // If way off then .05
     608              if (fakeCutoff_<=-1.0e100) {
     609                // use value as percentage - so 100==0.0, 101==1.0 etc
     610                // probably something like pow I could use but ...
     611                double fraction = 0.0;
     612                while (fakeCutoff_<-1.01e100) {
     613                  fakeCutoff_ *= 0.1;
     614                  fraction += 0.01;
     615                }
     616                rhs = solver->getObjValue()+fraction*fabs(solver->getObjValue());
     617              } else {
     618                rhs = 2.0*solver->getObjValue()-fakeCutoff_; // flip difference
     619              }
     620              fakeCutoff_=COIN_DBL_MAX;
     621            }
    581622            const double * objective = solver->getObjCoefficients();
    582623            int numberColumns = solver->getNumCols();
     
    829870                            if (returnCode && newSolutionValue < saveValue)
    830871                                numberBandBsolutions++;
     872                        } else if (numberColumns>numberIntegersOrig) {
     873                          // relax continuous
     874                          bool takeHint;
     875                          OsiHintStrength strength;
     876                          solver->getHintParam(OsiDoDualInResolve, takeHint, strength);
     877                          //solver->setHintParam(OsiDoReducePrint, false, OsiHintTry);
     878                          solver->setHintParam(OsiDoDualInResolve, false, OsiHintDo);
     879                          //solver->setHintParam(OsiDoScale, false, OsiHintDo);
     880                          solver->resolve();
     881                          solver->setHintParam(OsiDoDualInResolve, takeHint, strength);
     882                          if (solver->isProvenOptimal()) {
     883                            memcpy(newSolution,solver->getColSolution(),
     884                                   numberColumns*sizeof(double));
     885                            newSolutionValue = -saveOffset;
     886                            for (  i = 0 ; i < numberColumns ; i++ ) {
     887                              newSolutionValue += saveObjective[i] * newSolution[i];
     888                            }
     889                            newSolutionValue *= direction;
     890                            sprintf(pumpPrint, "Relaxing continuous gives %g", newSolutionValue);
     891                            //#define DEBUG_BEST
     892#ifdef DEBUG_BEST
     893                            {
     894                              int numberColumns=solver->getNumCols();
     895                              FILE * fp = fopen("solution.data2","wb");
     896                              printf("Solution data on file solution.data2\n");
     897                              size_t numberWritten;
     898                              numberWritten=fwrite(&numberColumns,sizeof(int),1,fp);
     899                              assert (numberWritten==1);
     900                              numberWritten=fwrite(&newSolutionValue,sizeof(double),1,fp);
     901                              assert (numberWritten==1);
     902                              numberWritten=fwrite(newSolution,sizeof(double),numberColumns,fp);
     903                              assert (numberWritten==numberColumns);
     904                              fclose(fp);
     905                              const double * rowLower = solver->getRowLower();
     906                              const double * rowUpper = solver->getRowUpper();
     907                              const double * columnLower = solver->getColLower();
     908                              const double * columnUpper = solver->getColUpper();
     909                              int numberRows = solver->getNumRows() ;
     910                              double *rowActivity = new double[numberRows] ;
     911                              memset(rowActivity, 0, numberRows*sizeof(double)) ;
     912                              const double * element = solver->getMatrixByCol()->getElements();
     913                              const int * row = solver->getMatrixByCol()->getIndices();
     914                              const CoinBigIndex * columnStart = solver->getMatrixByCol()->getVectorStarts();
     915                              const int * columnLength = solver->getMatrixByCol()->getVectorLengths();
     916                              double largestAway=0.0;
     917                              int away=-1;
     918                              double saveOffset;
     919                              solver->getDblParam(OsiObjOffset, saveOffset);
     920                              double newSolutionValue = -saveOffset;
     921                              const double * objective = solver->getObjCoefficients();
     922                              for ( int iColumn=0 ; iColumn<numberColumns ; ++iColumn ) {
     923                                double value=newSolution[iColumn];
     924                                CoinBigIndex start = columnStart[iColumn];
     925                                CoinBigIndex end = start + columnLength[iColumn];
     926                                for (CoinBigIndex j = start; j < end; j++) {
     927                                  int iRow = row[j];
     928                                  if (iRow==1996)
     929                                    printf("fp col %d val %g el %g old y %g\n",
     930                                  iColumn,value,element[j],rowActivity[iRow]);
     931                                  rowActivity[iRow] += value * element[j];
     932                                }
     933                                newSolutionValue += objective[iColumn] * newSolution[iColumn];
     934                                if (solver->isInteger(iColumn)) {
     935                                  double intValue = floor(value+0.5);
     936                                  if (fabs(value-intValue)>largestAway) {
     937                                    largestAway=fabs(value-intValue);
     938                                    away=iColumn;
     939                                  }
     940                                }
     941                              }
     942                              printf("Largest away from int at column %d was %g - obj %g\n",away,
     943                                     largestAway,newSolutionValue);
     944                              double largestInfeasibility=0.0;
     945                              for (int i = 0 ; i < numberRows ; i++) {
     946#if 0 //def CLP_INVESTIGATE
     947                                double inf;
     948                                inf = rowLower[i] - rowActivity[i];
     949                                if (inf > primalTolerance)
     950                                  printf("Row %d inf %g sum %g %g <= %g <= %g\n",
     951                                         i, inf, rowSum[i], rowLower[i], rowActivity[i], rowUpper[i]);
     952                                inf = rowActivity[i] - rowUpper[i];
     953                                if (inf > primalTolerance)
     954                                  printf("Row %d inf %g %g <= %g <= %g\n",
     955                                         i, inf, rowLower[i], rowActivity[i], rowUpper[i]);
     956#endif
     957                                double infeasibility = CoinMax(rowActivity[i]-rowUpper[i],
     958                                                               rowLower[i]-rowActivity[i]);
     959                                if (infeasibility>largestInfeasibility) {
     960                                  largestInfeasibility = infeasibility;
     961                                  printf("Binf of %g on row %d\n",
     962                                         infeasibility,i);
     963                                }
     964                              }
     965                              delete [] rowActivity ;
     966                              printf("Blargest infeasibility is %g - obj %g\n", largestInfeasibility,newSolutionValue);
     967                            }
     968#endif
     969                          } else {
     970                            sprintf(pumpPrint,"Infeasible when relaxing continuous!\n");
     971                          }
     972                          model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
     973                            << pumpPrint
     974                            << CoinMessageEol;
    831975                        }
    832976                    }
     
    9911135                int numberChanged = 0;
    9921136                const double * oldObjective = solver->getObjCoefficients();
     1137                bool fixOnesAtBound=false;
     1138                if (tryOneClosePass&&numberPasses==2) {
     1139                  // take off
     1140                  tryOneClosePass=false;
     1141                  int n=solver->getNumRows()-1;
     1142                  double rhs = solver->getRowUpper()[n];
     1143                  solver->setRowUpper(n,rhs+1.0e15);
     1144                  useRhs+=1.0e15;
     1145                  fixOnesAtBound=true;
     1146                }
    9931147                for (i = 0; i < numberColumns; i++) {
    9941148                    // below so we can keep original code and allow for objective
     
    10091163                    double newValue = 0.0;
    10101164                    if (newSolution[iColumn] < lower[iColumn] + primalTolerance) {
    1011                         newValue = costValue + scaleFactor * saveObjective[iColumn];
     1165                      newValue = costValue + scaleFactor * saveObjective[iColumn];
     1166                      if (fixOnesAtBound)
     1167                        newValue = 100.0*costValue;
    10121168                    } else {
    1013                         if (newSolution[iColumn] > upper[iColumn] - primalTolerance) {
    1014                             newValue = -costValue + scaleFactor * saveObjective[iColumn];
    1015                         }
     1169                      if (newSolution[iColumn] > upper[iColumn] - primalTolerance) {
     1170                        newValue = -costValue + scaleFactor * saveObjective[iColumn];
     1171                      if (fixOnesAtBound)
     1172                        newValue = -100.0*costValue;
     1173                      }
    10161174                    }
    10171175#ifdef RAND_RAND
     
    12361394                    {
    12371395                        const double * newSolution = solver->getColSolution();
    1238                         for (  i = 0 ; i < numberColumns ; i++ ) {
    1239                             if (solver->isInteger(i)) {
    1240                                 double value = newSolution[i];
     1396                        for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++ ) {
     1397                            if (solver->isInteger(iColumn)) {
     1398                                double value = newSolution[iColumn];
    12411399                                double nearest = floor(value + 0.5);
    12421400                                newSumInfeas += fabs(value - nearest);
    1243                                 if (fabs(value - nearest) > 1.0e-6)
     1401                                if (fabs(value - nearest) > 1.0e-6) {
    12441402                                    newNumberInfeas++;
    1245                             }
    1246                             newTrueSolutionValue += saveObjective[i] * newSolution[i];
     1403                                }
     1404                            }
     1405                            newTrueSolutionValue += saveObjective[iColumn] * newSolution[iColumn];
    12471406                        }
    12481407                        newTrueSolutionValue *= direction;
     
    15461705                << pumpPrint
    15471706                << CoinMessageEol;
     1707                CbcEventHandler *eventHandler = model_->getEventHandler() ;
     1708                if (eventHandler) {
     1709                  typedef struct {
     1710                    double newSumInfeas;
     1711                    double trueSolutionValue;
     1712                    double spareDouble[2];
     1713                    OsiSolverInterface * solver;
     1714                    void * sparePointer[2];
     1715                    int numberPasses;
     1716                    int totalNumberPasses;
     1717                    int numberInfeas;
     1718                    int numberIterations;
     1719                    int spareInt[3];
     1720                  } HeurPass;
     1721                  HeurPass temp;
     1722                  temp.solver=solver;
     1723                  temp.newSumInfeas = newSumInfeas;
     1724                  temp.trueSolutionValue = newTrueSolutionValue;
     1725                  temp.numberPasses=numberPasses;
     1726                  temp.totalNumberPasses=totalNumberPasses;
     1727                  temp.numberInfeas=newNumberInfeas;
     1728                  temp.numberIterations=numberIterations;
     1729                  CbcEventHandler::CbcAction status =
     1730                    eventHandler->event(CbcEventHandler::heuristicPass,
     1731                                        &temp);
     1732                  if (status==CbcEventHandler::killSolution) {
     1733                    exitAll = true;
     1734                    break;
     1735                  }
     1736                }
    15481737                if (closestSolution && solver->getObjValue() < closestObjectiveValue) {
    15491738                    int i;
     
    19382127                            if (nFixed*3 > numberColumns*2)
    19392128                                simplex->allSlackBasis(); // may as well go from all slack
     2129                            int logLevel=simplex->logLevel();
     2130                            if (logLevel<=1)
     2131                              simplex->setLogLevel(0);
    19402132                            simplex->primal(1);
     2133                            simplex->setLogLevel(logLevel);
    19412134                            clpSolver->setWarmStart(NULL);
    19422135                        }
     
    23642557    const double * columnUpper = solver->getColUpper();
    23652558    // Check if valid with current solution (allow for 0.99999999s)
     2559    double newSumInfeas = 0.0;
     2560    int newNumberInfeas = 0;
    23662561    for (i = 0; i < numberIntegers; i++) {
    23672562        int iColumn = integerVariable[i];
    23682563        double value = solution[iColumn];
    23692564        double round = floor(value + 0.5);
    2370         if (fabs(value - round) > primalTolerance)
    2371             break;
     2565        if (fabs(value - round) > primalTolerance) {
     2566          newSumInfeas += fabs(value-round);
     2567          newNumberInfeas++;
     2568        }
    23722569    }
    2373     if (i == numberIntegers) {
     2570    if (!newNumberInfeas) {
    23742571        // may be able to use solution even if 0.99999's
    23752572        double * saveLower = CoinCopyOfArray(columnLower, numberColumns);
Note: See TracChangeset for help on using the changeset viewer.