Ignore:
Timestamp:
Sep 28, 2007 11:40:51 AM (12 years ago)
Author:
forrest
Message:

add ampl stuff after heuristic

File:
1 edited

Legend:

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

    r789 r809  
    391391  const double *lower = getColLower() ;
    392392  const double *upper = getColUpper() ;
     393  /*
     394    Scan continuous and integer variables to see if continuous
     395    are cover or network with integral rhs.
     396  */
     397  double continuousMultiplier = 1.0;
     398  double * coeffMultiplier=NULL;
     399  {
     400    const double *rowLower = getRowLower() ;
     401    const double *rowUpper = getRowUpper() ;
     402    int numberRows = solver_->getNumRows() ;
     403    double * rhs = new double [numberRows];
     404    memset(rhs,0,numberRows*sizeof(double));
     405    int iColumn;
     406    int numberColumns = solver_->getNumCols() ;
     407    // Column copy of matrix
     408    const double * element = solver_->getMatrixByCol()->getElements();
     409    const int * row = solver_->getMatrixByCol()->getIndices();
     410    const CoinBigIndex * columnStart = solver_->getMatrixByCol()->getVectorStarts();
     411    const int * columnLength = solver_->getMatrixByCol()->getVectorLengths();
     412    for (iColumn=0;iColumn<numberColumns;iColumn++) {
     413      if (upper[iColumn]==lower[iColumn]) {
     414        CoinBigIndex start = columnStart[iColumn];
     415        CoinBigIndex end = start + columnLength[iColumn];
     416        for (CoinBigIndex j=start;j<end;j++) {
     417          int iRow = row[j];
     418          rhs[iRow] += lower[iColumn]*element[j];
     419        }
     420      }
     421    }
     422    int iRow;
     423    for (iRow=0;iRow<numberRows;iRow++) {
     424      if (rowLower[iRow]>-1.0e20&&
     425          fabs(rowLower[iRow]-rhs[iRow]-floor(rowLower[iRow]-rhs[iRow]+0.5))>1.0e-10) {
     426        continuousMultiplier=0.0;
     427        break;
     428      }
     429      if (rowUpper[iRow]<1.0e20&&
     430          fabs(rowUpper[iRow]-rhs[iRow]-floor(rowUpper[iRow]-rhs[iRow]+0.5))>1.0e-10) {
     431        continuousMultiplier=0.0;
     432        break;
     433      }
     434      // set rhs to limiting value
     435      if (rowLower[iRow]!=rowUpper[iRow]) {
     436        if(rowLower[iRow]>-1.0e20) {
     437          if (rowUpper[iRow]<1.0e20) {
     438            // no good
     439            continuousMultiplier=0.0;
     440            break;
     441          } else {
     442            rhs[iRow] = rowLower[iRow]-rhs[iRow];
     443          }
     444        } else {
     445          rhs[iRow] = rowUpper[iRow]-rhs[iRow];
     446        }
     447      } else {
     448        rhs[iRow] = rowUpper[iRow]-rhs[iRow];
     449      }
     450    }
     451    if (continuousMultiplier) {
     452      // 1 network, 2 cover, 4 negative cover
     453      int possible=7;
     454      bool unitRhs=true;
     455      // See which rows could be set cover
     456      for (iColumn=0;iColumn<numberColumns;iColumn++) {
     457        if (upper[iColumn] > lower[iColumn]+1.0e-8) {
     458          CoinBigIndex start = columnStart[iColumn];
     459          CoinBigIndex end = start + columnLength[iColumn];
     460          for (CoinBigIndex j=start;j<end;j++) {
     461            double value = element[j];
     462            if (value!=1.0)
     463              rhs[row[j]]=-0.5;
     464            else if (value!=-1.0)
     465              rhs[row[j]]=-COIN_DBL_MAX;
     466          }
     467        }
     468      }
     469      for (iColumn=0;iColumn<numberColumns;iColumn++) {
     470        if (upper[iColumn] > lower[iColumn]+1.0e-8) {
     471          if (!isInteger(iColumn)) {
     472            CoinBigIndex start = columnStart[iColumn];
     473            CoinBigIndex end = start + columnLength[iColumn];
     474            double rhsValue=0.0;
     475            // 1 all ones, -1 all -1s, 2 all +- 1, 3 no good
     476            int type=0;
     477            for (CoinBigIndex j=start;j<end;j++) {
     478              double value = element[j];
     479              if (fabs(value!=1.0)) {
     480                type=3;
     481                break;
     482              } else if (value==1.0) {
     483                if (!type)
     484                  type=1;
     485                else if (type!=1)
     486                  type=2;
     487              } else {
     488                if (!type)
     489                  type=-1;
     490                else if (type!=-1)
     491                  type=2;
     492              }
     493              int iRow = row[j];
     494              if (rhs[iRow]==-COIN_DBL_MAX) {
     495                type=3;
     496                break;
     497              } else if (rhs[iRow]==-0.5) {
     498                // different values
     499                unitRhs=false;
     500              } else if (rhsValue) {
     501                if (rhsValue!=rhs[iRow])
     502                  unitRhs=false;
     503              } else {
     504                rhsValue=rhs[iRow];
     505              }
     506            }
     507            // if no elements OK
     508            if (type==3) {
     509              // no good
     510              possible=0;
     511              break;
     512            } else if (type==2) {
     513              if (end-start>2) {
     514                // no good
     515                possible=0;
     516                break;
     517              } else {
     518                // only network
     519                possible &= 1;
     520                if (!possible)
     521                  break;
     522              }
     523            } else if (type==1) {
     524              // only cover
     525              possible &= 2;
     526              if (!possible)
     527                break;
     528            } else if (type==-1) {
     529              // only negative cover
     530              possible &= 4;
     531              if (!possible)
     532                break;
     533            }
     534          }
     535        }
     536      }
     537      if ((possible==2||possible==4)&&!unitRhs) {
     538        printf("XXXXXX Continuous all +1 but different rhs\n");
     539        possible=0;
     540      }
     541      // may be all integer
     542      if (possible!=7) {
     543        if (!possible)
     544          continuousMultiplier=0.0;
     545        else if (possible==1)
     546          continuousMultiplier=1.0;
     547        else
     548          continuousMultiplier=0.5;
     549        if (continuousMultiplier)
     550          printf("XXXXXX multiplier of %g\n",continuousMultiplier);
     551        if (continuousMultiplier==0.5) {
     552          coeffMultiplier=new double [numberColumns];
     553          bool allOne=true;
     554          for (iColumn=0;iColumn<numberColumns;iColumn++) {
     555            coeffMultiplier[iColumn]=1.0;
     556            if (upper[iColumn] > lower[iColumn]+1.0e-8) {
     557              if (!isInteger(iColumn)) {
     558                CoinBigIndex start = columnStart[iColumn];
     559                int iRow = row[start];
     560                double value = rhs[iRow];
     561                assert (value>=0.0);
     562                if (value!=0.0&&value!=1.0)
     563                  allOne=false;
     564                coeffMultiplier[iColumn]=0.5*value;
     565              }
     566            }
     567          }
     568          if (allOne) {
     569            // back to old way
     570            delete [] coeffMultiplier;
     571            coeffMultiplier=NULL;
     572          }
     573        }
     574      }
     575    }
     576    delete [] rhs;
     577  }
    393578/*
    394579  Take a first scan to see if there are unfixed continuous variables in the
     
    401586  double maximumCost = 0.0 ;
    402587  trueIncrement=0.0;
    403   bool possibleMultiple = true ;
     588  bool possibleMultiple = continuousMultiplier!=0.0 ;
    404589  int iColumn ;
    405590  int numberColumns = getNumCols() ;
    406   for (iColumn = 0 ; iColumn < numberColumns ; iColumn++)
    407   { if (upper[iColumn] > lower[iColumn]+1.0e-8)
    408     { if (isInteger(iColumn))
    409         maximumCost = CoinMax(maximumCost,fabs(objective[iColumn])) ;
    410       else if (objective[iColumn])
    411         possibleMultiple = false ; } }
     591  if (possibleMultiple) {
     592    for (iColumn = 0 ; iColumn < numberColumns ; iColumn++)
     593      { if (upper[iColumn] > lower[iColumn]+1.0e-8)
     594          { maximumCost = CoinMax(maximumCost,fabs(objective[iColumn])) ; } }
     595  }
    412596  setIntParam(CbcModel::CbcFathomDiscipline,possibleMultiple) ;
    413597/*
     
    428612    for (iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
    429613      if (upper[iColumn] > lower[iColumn]+1.0e-8) {
    430         if (isInteger(iColumn)&&objective[iColumn]) {
    431           double value = fabs(objective[iColumn])*multiplier ;
     614        double objValue = fabs(objective[iColumn]);
     615        if (!isInteger(iColumn)) {
     616          if (!coeffMultiplier)
     617            objValue *= continuousMultiplier;
     618          else
     619            objValue *= coeffMultiplier[iColumn];
     620        }
     621        if (objValue) {
     622          double value = objValue*multiplier ;
    432623          if (value <2.1e9) {
    433624            int nearest = (int) floor(value+0.5) ;
     
    441632          } else {
    442633            // large value - may still be multiple of 1.0
    443             value = fabs(objective[iColumn]);
    444             if (fabs(value-floor(value+0.5)) > 1.0e-8) {
     634            if (fabs(objValue-floor(objValue+0.5)) > 1.0e-8) {
    445635              increment=0;
    446636              break;
     
    452642      }
    453643    }
    454 
     644    delete [] coeffMultiplier;
    455645/*
    456646  If the increment beats the current value for objective change, install it.
Note: See TracChangeset for help on using the changeset viewer.