Ignore:
Timestamp:
Oct 20, 2006 3:59:21 PM (14 years ago)
Author:
forrest
Message:

for osi methods

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/devel/Cbc/src/CbcHeuristicFPump.cpp

    r444 r463  
    2727   absoluteIncrement_(0.0),
    2828   relativeIncrement_(0.0),
     29   initialWeight_(0.0),
     30   weightFactor_(0.1),
    2931   maximumPasses_(100),
    3032   maximumRetries_(1),
    … …  
    4547   absoluteIncrement_(0.0),
    4648   relativeIncrement_(0.0),
     49   initialWeight_(0.0),
     50   weightFactor_(0.1),
    4751   maximumPasses_(100),
    4852   maximumRetries_(1),
    … …  
    104108    fprintf(fp,"4  heuristicFPump.setRelativeIncrement(%g);\n",relativeIncrement_);
    105109  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicFPump);\n");
     110  if (initialWeight_!=other.initialWeight_)
     111    fprintf(fp,"3  heuristicFPump.setInitialWeight(%g);\n",initialWeight_);
     112  else
     113    fprintf(fp,"4  heuristicFPump.setInitialWeight(%g);\n",initialWeight_);
     114  if (weightFactor_!=other.weightFactor_)
     115    fprintf(fp,"3  heuristicFPump.setWeightFactor(%g);\n",weightFactor_);
     116  else
     117    fprintf(fp,"4  heuristicFPump.setWeightFactor(%g);\n",weightFactor_);
    106118}
    107119
    … …  
    116128  absoluteIncrement_(rhs.absoluteIncrement_),
    117129  relativeIncrement_(rhs.relativeIncrement_),
     130  initialWeight_(rhs.initialWeight_),
     131  weightFactor_(rhs.weightFactor_),
    118132  maximumPasses_(rhs.maximumPasses_),
    119133  maximumRetries_(rhs.maximumRetries_),
    … …  
    154168  cutoff *= direction;
    155169  cutoff = CoinMin(cutoff,solutionValue);
    156   // space for rounded solution
     170  // check plausible and space for rounded solution
    157171  int numberColumns = model_->getNumCols();
     172  int numberIntegers = model_->numberIntegers();
     173  const int * integerVariableOrig = model_->integerVariable();
     174
     175  // 1. initially check 0-1
     176  int i,j;
     177  int general=0;
     178  int * integerVariable = new int[numberIntegers];
     179  const double * lower = model_->solver()->getColLower();
     180  const double * upper = model_->solver()->getColUpper();
     181  bool doGeneral = (accumulate_&4)!=0;
     182  j=0;
     183  for (i=0;i<numberIntegers;i++) {
     184    int iColumn = integerVariableOrig[i];
     185#ifndef NDEBUG
     186    const OsiObject * object = model_->object(i);
     187    const CbcSimpleInteger * integerObject =
     188      dynamic_cast<const  CbcSimpleInteger *> (object);
     189    const OsiSimpleInteger * integerObject2 =
     190      dynamic_cast<const  OsiSimpleInteger *> (object);
     191    assert(integerObject||integerObject2);
     192#endif
     193    if (upper[iColumn]-lower[iColumn]>1.000001) {
     194      general++;
     195      if (doGeneral)
     196        integerVariable[j++]=iColumn;
     197    } else {
     198      integerVariable[j++]=iColumn;
     199    }
     200  }
     201  if (general*3>numberIntegers&&!doGeneral) {
     202    delete [] integerVariable;
     203    return 0;
     204  }
     205  int numberIntegersOrig = numberIntegers;
     206  numberIntegers = j;
     207  if (!doGeneral)
     208    general=0;
    158209  double * newSolution = new double [numberColumns];
    159210  double newSolutionValue=COIN_DBL_MAX;
    … …  
    219270    solver->getDblParam(OsiPrimalTolerance,primalTolerance);
    220271   
    221     int numberIntegers = model_->numberIntegers();
    222     const int * integerVariable = model_->integerVariable();
    223 
    224     // 1. initially check 0-1
    225     int i,j;
    226     int general=0;
    227     for (i=0;i<numberIntegers;i++) {
    228       int iColumn = integerVariable[i];
    229 #ifndef NDEBUG
    230       const OsiObject * object = model_->object(i);
    231       const CbcSimpleInteger * integerObject =
    232         dynamic_cast<const  CbcSimpleInteger *> (object);
    233       const OsiSimpleInteger * integerObject2 =
    234         dynamic_cast<const  OsiSimpleInteger *> (object);
    235       assert(integerObject||integerObject2);
    236 #endif
    237       if (upper[iColumn]-lower[iColumn]>1.000001) {
    238         general++;
    239         break;
    240       }
    241     }
    242     if (general*3>numberIntegers) {
    243       delete solver;
    244       return 0;
    245     }
    246272   
    247273    // 2 space for last rounded solutions
    … …  
    271297    double saveOffset;
    272298    solver->getDblParam(OsiObjOffset,saveOffset);
    273    
     299    // Get amount for original objective
     300    double scaleFactor = 0.0;
     301    for (i=0;i<numberColumns;i++)
     302      scaleFactor += saveObjective[i]*saveObjective[i];
     303    if (scaleFactor)
     304      scaleFactor = (initialWeight_*sqrt((double) numberIntegers))/sqrt(scaleFactor);
    274305    // 5. MAIN WHILE LOOP
    275306    bool newLineNeeded=false;
    … …  
    291322      memcpy(newSolution,solution,numberColumns*sizeof(double));
    292323      int flip;
    293       returnCode = rounds(newSolution,saveObjective,roundExpensive_,downValue_,&flip);
     324      returnCode = rounds(newSolution,saveObjective,numberIntegers,integerVariable,
     325                          roundExpensive_,downValue_,&flip);
    294326      if (returnCode) {
    295327        // SOLUTION IS INTEGER
    … …  
    312344          if (general) {
    313345            int numberLeft=0;
    314             for (i=0;i<numberIntegers;i++) {
    315               int iColumn = integerVariable[i];
     346            for (i=0;i<numberIntegersOrig;i++) {
     347              int iColumn = integerVariableOrig[i];
    316348              double value = floor(newSolution[iColumn]+0.5);
    317               if(solver->isBinary(iColumn)) {
     349              if(solver->isBinary(iColumn)||general) {
    318350                solver->setColLower(iColumn,value);
    319351                solver->setColUpper(iColumn,value);
    … …  
    358390          for (i = 0; i <numberIntegers; i++) {
    359391            int iColumn = integerVariable[i];
    360             if(!solver->isBinary(iColumn))
    361               continue;
    362392            if (newSolution[iColumn]!=b[iColumn]) {
    363393              matched=false;
    … …  
    374404          for (i=0;i<numberIntegers;i++) {
    375405            int iColumn = integerVariable[i];
    376             if(!solver->isBinary(iColumn))
    377               continue;
    378406            double value = max(0.0,CoinDrand48()-0.3);
    379407            double difference = fabs(solution[iColumn]-newSolution[iColumn]);
    380408            if (difference+value>0.5) {
    381               if (newSolution[iColumn]<lower[iColumn]+primalTolerance) newSolution[iColumn] += 1.0;
    382               else if (newSolution[iColumn]>upper[iColumn]-primalTolerance) newSolution[iColumn] -= 1.0;
    383               else abort();
     409              if (newSolution[iColumn]<lower[iColumn]+primalTolerance) {
     410                newSolution[iColumn] += 1.0;
     411              } else if (newSolution[iColumn]>upper[iColumn]-primalTolerance) {
     412                newSolution[iColumn] -= 1.0;
     413              } else {
     414                // general integer
     415                if (difference+value>0.75)
     416                  newSolution[iColumn] += 1.0;
     417                else
     418                  newSolution[iColumn] -= 1.0;
     419              }
    384420            }
    385421          }
    … …  
    393429        // 2. update the objective function based on the new rounded solution
    394430        double offset=0.0;
    395         for (i=0;i<numberIntegers;i++) {
    396           int iColumn = integerVariable[i];
    397           if(!solver->isBinary(iColumn))
     431        double costValue = (1.0-scaleFactor)*solver->getObjSense();
     432       
     433        for (i=0;i<numberColumns;i++) {
     434          // below so we can keep original code and allow for objective
     435          int iColumn = i;
     436          if(!solver->isBinary(iColumn)&&!general)
    398437            continue;
    399           double costValue = 1.0;
    400438          // deal with fixed variables (i.e., upper=lower)
    401439          if (fabs(lower[iColumn]-upper[iColumn]) < primalTolerance) {
    402             if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue);
    403             else                                       solver->setObjCoeff(iColumn,costValue);
     440            //if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue);
     441            //else                                       solver->setObjCoeff(iColumn,costValue);
    404442            continue;
    405443          }
    406444          if (newSolution[iColumn]<lower[iColumn]+primalTolerance) {
    407             solver->setObjCoeff(iColumn,costValue);
     445            solver->setObjCoeff(iColumn,costValue+scaleFactor*saveObjective[iColumn]);
    408446          } else {
    409447            if (newSolution[iColumn]>upper[iColumn]-primalTolerance) {
    410               solver->setObjCoeff(iColumn,-costValue);
     448              solver->setObjCoeff(iColumn,-costValue+scaleFactor*saveObjective[iColumn]);
    411449            } else {
    412450              abort();
    … …  
    416454        }
    417455        solver->setDblParam(OsiObjOffset,-offset);
    418         if (!general&false) {
     456        if (!general&&false) {
    419457          // Solve in two goes - first keep satisfied ones fixed
    420458          double * saveLower = new double [numberIntegers];
    … …  
    444482          memcpy(newSolution,solution,numberColumns*sizeof(double));
    445483          int flip;
    446           returnCode = rounds(newSolution,saveObjective,roundExpensive_,downValue_,&flip);
     484          returnCode = rounds(newSolution,saveObjective,numberIntegers,integerVariable,
     485                              roundExpensive_,downValue_,&flip);
    447486          if (returnCode) {
    448487            // solution - but may not be better
    … …  
    465504          }     
    466505        }
    467         solver->resolve();
    468         assert (solver->isProvenOptimal());
     506        if (!general) {
     507          solver->resolve();
     508          assert (solver->isProvenOptimal());
     509        } else {
     510          int * addStart = new int[2*general+1];
     511          int * addIndex = new int[4*general];
     512          double * addElement = new double[4*general];
     513          double * addLower = new double[2*general];
     514          double * addUpper = new double[2*general];
     515          double * obj = new double[general];
     516          int nAdd=0;
     517          for (i=0;i<numberIntegers;i++) {
     518            int iColumn = integerVariable[i];
     519            if (newSolution[iColumn]>lower[iColumn]+primalTolerance&&
     520                newSolution[iColumn]<upper[iColumn]-primalTolerance) {
     521              obj[nAdd]=1.0;
     522              addLower[nAdd]=0.0;
     523              addUpper[nAdd]=COIN_DBL_MAX;
     524              nAdd++;
     525            }
     526          }
     527          OsiSolverInterface * solver2 = solver;
     528          if (nAdd) {
     529            CoinZeroN(addStart,nAdd+1);
     530            solver2 = solver->clone();
     531            solver2->addCols(nAdd,addStart,NULL,NULL,addLower,addUpper,obj);
     532            // feasible solution
     533            double * sol = new double[nAdd+numberColumns];
     534            memcpy(sol,solution,numberColumns*sizeof(double));
     535            // now rows
     536            int nAdd=0;
     537            int nEl=0;
     538            int nAddRow=0;
     539            for (i=0;i<numberIntegers;i++) {
     540              int iColumn = integerVariable[i];
     541              if (newSolution[iColumn]>lower[iColumn]+primalTolerance&&
     542                  newSolution[iColumn]<upper[iColumn]-primalTolerance) {
     543                addLower[nAddRow]=-newSolution[iColumn];;
     544                addUpper[nAddRow]=COIN_DBL_MAX;
     545                addIndex[nEl] = iColumn;
     546                addElement[nEl++]=-1.0;
     547                addIndex[nEl] = numberColumns+nAdd;
     548                addElement[nEl++]=1.0;
     549                nAddRow++;
     550                addStart[nAddRow]=nEl;
     551                addLower[nAddRow]=newSolution[iColumn];;
     552                addUpper[nAddRow]=COIN_DBL_MAX;
     553                addIndex[nEl] = iColumn;
     554                addElement[nEl++]=1.0;
     555                addIndex[nEl] = numberColumns+nAdd;
     556                addElement[nEl++]=1.0;
     557                nAddRow++;
     558                addStart[nAddRow]=nEl;
     559                sol[nAdd+numberColumns] = fabs(sol[iColumn]-newSolution[iColumn]);
     560                nAdd++;
     561              }
     562            }
     563            solver2->setColSolution(sol);
     564            delete [] sol;
     565            solver2->addRows(nAddRow,addStart,addIndex,addElement,addLower,addUpper);
     566          }
     567          delete [] addStart;
     568          delete [] addIndex;
     569          delete [] addElement;
     570          delete [] addLower;
     571          delete [] addUpper;
     572          delete [] obj;
     573          solver2->resolve();
     574          assert (solver2->isProvenOptimal());
     575          if (nAdd) {
     576            solver->setColSolution(solver2->getColSolution());
     577            delete solver2;
     578          }
     579        }
    469580        if (model_->logLevel())
    470581          printf("\npass %3d: obj. %10.5f --> ", numberPasses+totalNumberPasses,solver->getObjValue());
    … …  
    472583       
    473584      }
     585      // reduce scale factor
     586      scaleFactor *= weightFactor_;
    474587    } // END WHILE
    475588    if (newLineNeeded&&model_->logLevel())
    … …  
    491604      for (i=0;i<numberIntegers;i++) {
    492605        int iColumn=integerVariable[i];
    493         const OsiObject * object = model_->object(i);
    494         double originalLower;
    495         double originalUpper;
    496         getIntegerInformation( object,originalLower, originalUpper);
    497         assert(colLower[iColumn]==originalLower);
    498         newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower));
    499         assert(colUpper[iColumn]==originalUpper);
    500         newSolver->setColUpper(iColumn,CoinMin(colUpper[iColumn],originalUpper));
     606        //const OsiObject * object = model_->object(i);
     607        //double originalLower;
     608        //double originalUpper;
     609        //getIntegerInformation( object,originalLower, originalUpper);
     610        //assert(colLower[iColumn]==originalLower);
     611        //newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower));
     612        newSolver->setColLower(iColumn,colLower[iColumn]);
     613        //assert(colUpper[iColumn]==originalUpper);
     614        //newSolver->setColUpper(iColumn,CoinMin(colUpper[iColumn],originalUpper));
     615        newSolver->setColUpper(iColumn,colUpper[iColumn]);
    501616        if (!usedColumn[iColumn]) {
    502617          double value=lastSolution[iColumn];
    … …  
    598713      cutoff -= CoinMax(CoinMax(gap,absoluteIncrement_),model_->getCutoffIncrement());
    599714      printf("round again with cutoff of %g\n",cutoff);
    600       if (accumulate_<2&&usedColumn)
     715      if ((accumulate_&3)<2&&usedColumn)
    601716        memset(usedColumn,0,numberColumns);
    602717      totalNumberPasses += numberPasses;
    … …  
    608723  delete [] lastSolution;
    609724  delete [] newSolution;
     725  delete [] integerVariable;
    610726  return finalReturnCode;
    611727}
    … …  
    625741CbcHeuristicFPump::rounds(double * solution,
    626742                          const double * objective,
     743                          int numberIntegers, const int * integerVariable,
    627744                          bool roundExpensive, double downValue, int *flip)
    628745{
    … …  
    631748  double primalTolerance ;
    632749  solver->getDblParam(OsiPrimalTolerance,primalTolerance) ;
    633   int numberIntegers = model_->numberIntegers();
    634   const int * integerVariable = model_->integerVariable();
    635750
    636751  int i;
    … …  
    653768  for (i=0;i<numberIntegers;i++) {
    654769    int iColumn = integerVariable[i];
    655     if(!solver->isBinary(iColumn))
    656       continue;
    657 #ifndef NDEBUG
    658       const OsiObject * object = model_->object(i);
    659       const CbcSimpleInteger * integerObject =
    660         dynamic_cast<const  CbcSimpleInteger *> (object);
    661       const OsiSimpleInteger * integerObject2 =
    662         dynamic_cast<const  OsiSimpleInteger *> (object);
    663       assert(integerObject||integerObject2);
    664 #endif
    665770    double value=solution[iColumn];
    666771    double round = floor(value+primalTolerance);
    … …  
    691796  delete [] tmp;
    692797
     798  const double * columnLower = solver->getColLower();
     799  const double * columnUpper = solver->getColUpper();
    693800  if (*flip == 0 && iter != 0) {
    694801    if(model_->logLevel())
    695802      printf(" -- rand = %4d (%4d) ", nnv, nn);
    696      for (i = 0; i < nnv; i++) solution[list[i]] = 1. - solution[list[i]];
     803     for (i = 0; i < nnv; i++) {
     804       // was solution[list[i]] = 1. - solution[list[i]]; but does that work for 7>=x>=6
     805       int index = list[i];
     806       double value = solution[index];
     807       if (value<=1.0) {
     808         solution[index] = 1.0-value;
     809       } else if (value<columnLower[index]+integerTolerance) {
     810         solution[index] = value+1.0;
     811       } else if (value>columnUpper[index]-integerTolerance) {
     812         solution[index] = value-1.0;
     813       } else {
     814         solution[index] = value-1.0;
     815       }
     816     }
    697817     *flip = nnv;
    698818  } else if (model_->logLevel()) {
Note: See TracChangeset for help on using the changeset viewer.