Changeset 1572


Ignore:
Timestamp:
May 11, 2020 1:07:02 PM (2 months ago)
Author:
forrest
Message:

try and improve preprocessing

Location:
trunk/src
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/CglPreProcess/CglPreProcess.cpp

    r1566 r1572  
    10801080  printf("saving %s from %s - %d row, %d columns\n",
    10811081    name, where, solver->getNumRows(), solver->getNumCols());
     1082  if (mpsNumber>20) {
     1083    printf("Not saving\n");
     1084    return;
     1085  }
    10821086  solver->writeMpsNative(name, NULL, NULL, 0, 1, 0);
    10831087#if CGL_WRITEMPS > 1
     
    11231127#define writeDebugMps(x, y, z)
    11241128#endif
     1129#define USE_CGL_RATIONAL 1
     1130#if USE_CGL_RATIONAL>0
     1131#include "CoinRational.hpp"
     1132static long computeGcd(long a, long b) {
     1133  // This is the standard Euclidean algorithm for gcd
     1134  long remainder = 1;
     1135  // Make sure a<=b (will always remain so)
     1136  if (a > b) {
     1137    // Swap a and b
     1138    long temp = a;
     1139    a = b;
     1140    b = temp;
     1141  }
     1142  // If zero then gcd is nonzero
     1143  if (!a) {
     1144    if (b) {
     1145      return b;
     1146    }
     1147    else {
     1148      printf("### WARNING: CglGMI::computeGcd() given two zeroes!\n");
     1149      exit(1);
     1150    }
     1151  }
     1152  while (remainder) {
     1153    remainder = b % a;
     1154    b = a;
     1155    a = remainder;
     1156  }
     1157  return b;
     1158} /* computeGcd */
     1159static bool scaleRowIntegral(double* rowElem, int rowNz)
     1160{
     1161  long gcd, lcm;
     1162  double maxdelta = 1.0e-13;
     1163  double maxscale = 1000;
     1164  long maxdnom = 1000;
     1165  //long numerator = 0, denominator = 0;
     1166  // Initialize gcd and lcm
     1167  CoinRational r = CoinRational(rowElem[0], maxdelta, maxdnom);
     1168  if (r.getNumerator() != 0){
     1169    gcd = labs(r.getNumerator());
     1170    lcm = r.getDenominator();
     1171  } else {
     1172    return false;
     1173  }
     1174  for (int i = 1; i < rowNz; ++i) {
     1175    if (rowElem[i]) {
     1176      r = CoinRational(rowElem[i], maxdelta, maxdnom);
     1177      if (r.getNumerator() != 0){
     1178        gcd = computeGcd(gcd, r.getNumerator());
     1179        lcm *= r.getDenominator()/(computeGcd(lcm,r.getDenominator()));
     1180      } else {
     1181        return false;
     1182      }
     1183    }
     1184  }
     1185  double scale = ((double)lcm)/((double)gcd);
     1186  if (fabs(scale) > maxscale) {
     1187    return false;
     1188  }
     1189  scale = fabs(scale);
     1190  // Looks like we have a good scaling factor; scale and return;
     1191  for (int i = 0; i < rowNz; ++i) {
     1192    double value = rowElem[i]*scale;
     1193    rowElem[i] = floor(value+0.5);
     1194    assert (fabs(rowElem[i]-value)<1.0e-9);
     1195  }
     1196  return true;
     1197} /* scaleRowIntegral */
     1198#endif
    11251199OsiSolverInterface *
    11261200CglPreProcess::preProcessNonDefault(OsiSolverInterface &model,
     
    11291203{
    11301204  double ppstart = getCurrentCPUTime();
    1131 #ifdef CGL_WRITEMPS
     1205#if 0 //def CGL_WRITEMPS
    11321206  bool rcdActive = true;
    11331207  std::string modelName;
     
    12231297  delete startModel_;
    12241298  startModel_ = originalModel_->clone();
    1225   CoinPackedMatrix matrixByRow(*originalModel_->getMatrixByRow());
    12261299  int numberRows = originalModel_->getNumRows();
     1300  int numberColumns = originalModel_->getNumCols();
     1301  int *rows = new int[numberRows];
     1302  double *element = new double[numberRows];
     1303  // probably okay with rowType_ but for now ..
     1304  CoinPackedMatrix matrixByRow(*originalModel_->getMutableMatrixByRow());
     1305  if (!rowType_ && true) {
     1306    // clean matrix
     1307    double *elementByRow = matrixByRow.getMutableElements();
     1308    const int *column = matrixByRow.getIndices();
     1309    const CoinBigIndex *rowStart = matrixByRow.getVectorStarts();
     1310    const int *rowLength = matrixByRow.getVectorLengths();
     1311   
     1312    const double *rowLower = originalModel_->getRowLower();
     1313    const double *rowUpper = originalModel_->getRowUpper();
     1314    double * temp = new double [2*numberColumns+4];
     1315    double * tempSave = temp+numberColumns+2;
     1316    int nChanged = 0;
     1317    int nScaled = 0;
     1318    for (int iRow=0;iRow<numberRows;iRow++) {
     1319      int n = 0;
     1320      // make majority positive
     1321      CoinBigIndex start = rowStart[iRow];
     1322      CoinBigIndex end = start+rowLength[iRow];
     1323      double multiplier = 1.0;
     1324      for (CoinBigIndex j=start;j<end;j++) {
     1325        if (elementByRow[j]<0)
     1326          n++;
     1327      }
     1328      if (2*n>end-start)
     1329        multiplier = -1.0;
     1330      int nInRow = end-start;
     1331      n = nInRow;
     1332      for (int i=0;i<n;i++)
     1333        temp[i] = multiplier*elementByRow[i+start];
     1334      double lower = rowLower[iRow];
     1335      double upper = rowUpper[iRow];
     1336      if (lower>-1.0e20)
     1337        temp[n++] = multiplier*lower;
     1338      if (upper<1.0e20)
     1339        temp[n++] = multiplier*upper;
     1340      memcpy(tempSave,temp,n*sizeof(double));
     1341      if (scaleRowIntegral(temp, n)) {
     1342        // double check
     1343        double largestError = 0.0;
     1344        double mult = temp[0]/elementByRow[start];
     1345        if (fabs(mult-floor(mult+0.1))<1.0e-12)
     1346          mult = floor(mult+0.1);
     1347        for (int i=0;i<n;i++) {
     1348          double value = mult*tempSave[i];
     1349          if (value) {
     1350            double vint = floor(value+0.01);
     1351            largestError = CoinMax(largestError,fabs(value-vint));
     1352            assert (fabs(vint)>0.9);
     1353          }
     1354        }
     1355        if (largestError<1.0e-9) {
     1356          multiplier = mult;
     1357        }
     1358      }
     1359      element[iRow] = multiplier;
     1360      if (multiplier!=1.0) {
     1361        nChanged++;
     1362        if (fabs(multiplier)!=1.0)
     1363          nScaled++;
     1364        memcpy(elementByRow+start,temp,(end-start)*sizeof(double));
     1365        if (multiplier<0.0) {
     1366          double tempV = lower;
     1367          lower = -upper;
     1368          upper = -tempV;
     1369        }
     1370        if (lower>-1.0e20)
     1371          lower *= fabs(multiplier);
     1372        if (upper<1.0e20)
     1373          upper *= fabs(multiplier);
     1374        originalModel_->setRowLower(iRow,lower);
     1375        originalModel_->setRowUpper(iRow,upper);
     1376      }
     1377    }
     1378    if (nChanged) {
     1379      CoinPackedMatrix * matrixByColumn =
     1380        originalModel_->getMutableMatrixByCol();
     1381      const int *row = matrixByColumn->getIndices();
     1382      const CoinBigIndex *columnStart = matrixByColumn->getVectorStarts();
     1383      const int *columnLength = matrixByColumn->getVectorLengths();
     1384      double *columnElements = matrixByColumn->getMutableElements();
     1385      double largestDelta = 0.0;
     1386      for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     1387        for (CoinBigIndex j=columnStart[iColumn];
     1388             j<columnStart[iColumn]+columnLength[iColumn];j++) {
     1389          int iRow = row[j];
     1390          double value = columnElements[j];
     1391          if (fabs(element[iRow])!=1.0) {
     1392            value *= element[iRow];
     1393            double vint = floor(value+0.01);
     1394            largestDelta = CoinMax(largestDelta,fabs(value-vint));
     1395            assert (largestDelta<1.0e-9);
     1396            columnElements[j] = vint;
     1397            assert (fabs(vint)>0.9);
     1398          } else if (element[iRow]==-1.0) {
     1399            columnElements[j] = -value;
     1400          }
     1401        }
     1402      }
     1403#if CBC_USEFUL_PRINTING > 1
     1404      printf("%d rows changed %d scaled - largest error %g\n",
     1405             nChanged,nScaled,largestDelta);
     1406#endif
     1407      //writeDebugMps(originalModel_, "scaled", 0);
     1408      delete startModel_;
     1409      startModel_ = originalModel_->clone();
     1410    }
     1411    delete [] temp;
     1412  }
    12271413  if (rowType_)
    12281414    assert(numberRowType_ == numberRows);
    1229   int numberColumns = originalModel_->getNumCols();
    12301415  //int originalNumberColumns=numberColumns;
    12311416  int minimumLength = 5;
     
    12501435  // We want to add columns
    12511436  int numberSlacks = 0;
    1252   int *rows = new int[numberRows];
    1253   double *element = new double[numberRows];
    12541437
    12551438  int iRow;
     
    19502133  numberRows = startModel_->getNumRows();
    19512134  numberColumns = startModel_->getNumCols();
     2135  // Column copy
     2136#ifdef CBC_PREPROCESS_EXPERIMENT
    19522137  //CoinPackedMatrix * matrixByColumn = const_cast<CoinPackedMatrix *>(startModel_->getMatrixByCol());
    1953   // Column copy
    1954   //const int * row = matrixByColumn->getIndices();
    1955   //const CoinBigIndex * columnStart = matrixByColumn->getVectorStarts();
    1956   //const int * columnLength = startModel_->getMatrixByCol()->getVectorLengths();
    1957   //const double * columnElements = matrixByColumn->getElements();
     2138  CoinPackedMatrix * matrixByColumn = startModel_->getMutableMatrixByCol();
     2139  const int * row = matrixByColumn->getIndices();
     2140  const CoinBigIndex * columnStart = matrixByColumn->getVectorStarts();
     2141  const int * columnLength = startModel_->getMatrixByCol()->getVectorLengths();
     2142  const double * columnElements = matrixByColumn->getElements();
     2143  // look for costed slacks on equality rows first
     2144  int nCosted = 0;
     2145  int nCostedI = 0;
     2146#define CAN_MOVE_LARGE_OBJ
     2147#ifndef CAN_MOVE_LARGE_OBJ
     2148  double goodRatio=0.001;
     2149#else
     2150  double goodRatio=0.0;
     2151#endif
     2152#endif
    19582153  double *obj = CoinCopyOfArray(startModel_->getObjCoefficients(), numberColumns);
     2154#ifdef CBC_PREPROCESS_EXPERIMENT
     2155#ifdef PRINT_STUFF
     2156  {
     2157    int * counts = new int[numberRows];
     2158    memset(counts,0,numberRows*sizeof(int));
     2159    int nn=0;
     2160    int nz=0;
     2161    int neq=0;
     2162    for (int i=0;i<numberColumns;i++) {
     2163      if (columnLength[i]==1) {
     2164        int iRow = row[columnStart[i]];
     2165        counts[iRow]++;
     2166        if (rowLower[iRow]==rowUpper[iRow])
     2167          neq++;
     2168        if (obj[i])
     2169          nn++;
     2170        else
     2171          nz++;
     2172      }
     2173    }
     2174    int cc[10];
     2175    memset(cc,0,sizeof(cc));
     2176    for (int i=0;i<numberRows;i++) {
     2177      if (rowLower[i]==rowUpper[i]) {
     2178        int k = CoinMin(counts[i],9);
     2179        cc[k]++;
     2180      }
     2181    }
     2182    for (int i=1;i<10;i++)
     2183      if (cc[i])
     2184        printf("(%d E rows have %d singletons) ",cc[i],i);
     2185    printf("\n");
     2186    memset(cc,0,sizeof(cc));
     2187    for (int i=0;i<numberRows;i++) {
     2188      if (rowLower[i]!=rowUpper[i]) {
     2189        int k = CoinMin(counts[i],9);
     2190        cc[k]++;
     2191      }
     2192    }
     2193    for (int i=1;i<10;i++)
     2194      if (cc[i])
     2195        printf("(%d L/G rows have %d singletons) ",cc[i],i);
     2196    printf("\n");
     2197    delete [] counts;
     2198    printf("%d single with cost %d zero cost %d equality rows\n",
     2199           nn,nz,neq);
     2200  }
     2201#endif
     2202#endif
    19592203  double offset;
    19602204  int numberMoved = 0;
    19612205  startModel_->getDblParam(OsiObjOffset, offset);
     2206#ifdef CBC_PREPROCESS_EXPERIMENT
     2207  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     2208    if (columnLength[iColumn]==1 && obj[iColumn]) {
     2209      double testObj = fabs(obj[iColumn])*goodRatio;
     2210      int iRow = row[columnStart[iColumn]];
     2211      if (rowLower[iRow]==rowUpper[iRow]) {
     2212        double rhs = rowLower[iRow];
     2213        bool canDo = false;
     2214        if (startModel_->isInteger(iColumn)) {
     2215          // are all integer
     2216          bool allInteger=true;
     2217          double element = columnElements[columnStart[iColumn]];
     2218          if ((element-floor(element)))
     2219            allInteger=false;
     2220          for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
     2221            int jColumn = column[j];
     2222            if (!startModel_->isInteger(jColumn)) {
     2223              allInteger=false;
     2224              break;
     2225            } else if (fabs(obj[jColumn])<testObj) {
     2226              // scaling
     2227              allInteger=false;
     2228              break;
     2229            }
     2230          }
     2231          if (allInteger) {
     2232            //printf("can take off cost on integer %d and make slack\n",iColumn);
     2233            nCostedI++;
     2234            canDo = true;
     2235            marked[iColumn]=1;
     2236            //printf("zz %d\n",iColumn);
     2237          }
     2238        } else {
     2239          bool goodScaling=true;
     2240          for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
     2241            int jColumn = column[j];
     2242            if (fabs(obj[jColumn])<testObj) {
     2243              // scaling
     2244              goodScaling=false;
     2245              break;
     2246            }
     2247          }
     2248          if (goodScaling) {
     2249            //printf("can take off cost on %d and make slack\n",iColumn);
     2250            nCosted++;
     2251            marked[iColumn]=1;
     2252            canDo = true;
     2253            //printf("zz %d\n",iColumn);
     2254          }
     2255        }
     2256        if (canDo) {
     2257          // make cost zero
     2258          double valueSlack = columnElements[columnStart[iColumn]];
     2259          double objValue = obj[iColumn];
     2260          double subtract = objValue/valueSlack;
     2261          offset -= subtract * rhs;
     2262          startModel_->setContinuous(iColumn); // to allow change
     2263          for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
     2264            int jColumn = column[j];
     2265            if (iColumn != jColumn) {
     2266              double value = elementByRow[j];
     2267              obj[jColumn] -= subtract*value;;
     2268            } else {
     2269              obj[iColumn] = 0.0;
     2270            }
     2271          }
     2272        }
     2273      }
     2274    }
     2275  }
     2276#ifdef PRINT_STUFF
     2277  if (nCosted||nCostedI)
     2278    printf("%d integer costed slacks and %d continuous\n",nCostedI,nCosted);
     2279  {
     2280    int nn=0;
     2281    int nz=0;
     2282    for (int i=0;i<numberColumns;i++) {
     2283      if (columnLength[i]==1) {
     2284        if (obj[i])
     2285          nn++;
     2286        else
     2287          nz++;
     2288      }
     2289    }
     2290    printf("%d single with cost %d zero cost\n",nn,nz);
     2291  }
     2292#endif
     2293#endif
    19622294  // This is not a vital loop so be careful
    19632295#ifndef SKIP_MOVE_COSTS_TO_INTEGERS
    19642296  for (iRow = 0; iRow < numberRows; iRow++) {
    1965     //int slack = -1;
    19662297    int nPlus = 0;
    19672298    int nMinus = 0;
     
    19742305    if (rhs != rowUpper[iRow])
    19752306      continue;
    1976     //int multiple=0;
    1977     //int iSlack=-1;
    19782307    int numberContinuous = 0;
    19792308    for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
    19802309      int iColumn = column[j];
     2310#ifdef CBC_PREPROCESS_EXPERIMENT
     2311      if (marked[iColumn]) {
     2312        nPlus = -numberColumns;
     2313        nMinus = -numberColumns;
     2314        break;
     2315      }
     2316#endif
    19812317      double value = elementByRow[j];
    19822318      if (upper[iColumn] > lower[iColumn]) {
     
    21242460    // allow transfer of costs
    21252461    //presolveActions |= 4; can be slow
     2462#ifdef CBC_PREPROCESS_EXPERIMENT
     2463    if ((tuning&8192)!=0)
     2464      presolveActions |= 0x40 | 0x80; // more for dupcol and doubleton
     2465#endif
    21262466    // If trying for SOS don't allow some transfers
    21272467    if (makeEquality == 2 || makeEquality == 3)
     
    22152555      << CoinMessageEol;
    22162556    return NULL;
     2557  } else {
     2558    //printf("skipping tightenPrimalBounds\n");
    22172559  }
    22182560  OsiSolverInterface *returnModel = NULL;
     
    25862928    }
    25872929  }
    2588 #if CGL_WRITEMPS
    2589   const OsiRowCutDebugger *debugger = returnModel->getRowCutDebugger();
    2590   if (debugger)
    2591     printf("Contains optimal before tighten\n");
    2592 #endif
    25932930  if (returnModel) {
    25942931    if (returnModel->getNumRows()) {
    25952932      // tighten bounds
    2596       int infeas = tightenPrimalBounds(*returnModel);
     2933#if CGL_WRITEMPS
     2934      const OsiRowCutDebugger *debugger = returnModel->getRowCutDebugger();
     2935      if (debugger)
     2936        printf("Contains optimal before tighten\n");
     2937#endif
     2938      int infeas = tightenPrimalBounds(*returnModel,true);
    25972939#if CGL_WRITEMPS
    25982940      if (debugger)
    25992941        assert(returnModel->getRowCutDebugger());
     2942      writeDebugMps(returnModel, "afterTighten", NULL);
    26002943#endif
    26012944      if (infeas) {
     2945        handler_->message(CGL_INFEASIBLE, messages_)
     2946          << CoinMessageEol;
    26022947        delete returnModel;
    26032948        for (int iPass = 0; iPass < numberSolvers_; iPass++) {
     
    30723417  return returnModel;
    30733418}
     3419static int
     3420tighten(double *colLower, double * colUpper,
     3421        const int *column, const double *rowElements,
     3422        const CoinBigIndex *rowStart,
     3423        const CoinBigIndex * rowStartPos,const int * rowLength,
     3424        double *rowLower, double *rowUpper,
     3425        int nRows,int nCols,char * intVar,int maxpass,
     3426        double tolerance);
    30743427
    30753428/* Tightens primal bounds to make dual and branch and cutfaster.  Unless
    30763429   fixed, bounds are slightly looser than they could be.
    30773430   Returns non-zero if problem infeasible
    3078    Fudge for branch and bound - put bounds on columns of factor *
    3079    largest value (at continuous) - should improve stability
    3080    in branch and bound on infeasible branches (0.0 is off)
    30813431*/
    3082 int CglPreProcess::tightenPrimalBounds(OsiSolverInterface &model, double factor)
     3432int
     3433CglPreProcess::tightenPrimalBounds(OsiSolverInterface &model,
     3434                                   bool tightenRowBounds)
    30833435{
    3084 
     3436  //return 0;
    30853437  // Get a row copy in standard format
    30863438  CoinPackedMatrix copy = *model.getMatrixByRow();
    30873439  // get matrix data pointers
    3088   const int *column = copy.getIndices();
     3440  int *column = copy.getMutableIndices();
    30893441  const CoinBigIndex *rowStart = copy.getVectorStarts();
    30903442  const int *rowLength = copy.getVectorLengths();
     
    31013453  const double *colLower = model.getColLower();
    31023454  const double *colUpper = model.getColUpper();
     3455#ifdef CBC_PREPROCESS_EXPERIMENT
     3456  const double *objective = model.getObjCoefficients();
     3457#endif
     3458  double direction = model.getObjSense();
     3459  const int *columnLength = model.getMatrixByCol()->getVectorLengths();
    31033460  // New and saved column bounds
    31043461  double *newLower = new double[numberColumns];
     
    31103467  double *columnUpper = new double[numberColumns];
    31113468  memcpy(columnUpper, colUpper, numberColumns * sizeof(double));
    3112 
    31133469  int iRow, iColumn;
    31143470
    3115   // If wanted - tighten column bounds using solution
    3116   if (factor) {
    3117     /*
    3118       Callers need to ensure that the solution is fresh
    3119     */
    3120     const double *solution = model.getColSolution();
    3121     double largest = 0.0;
    3122     if (factor > 0.0) {
    3123       assert(factor > 1.0);
    3124       for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    3125         if (columnUpper[iColumn] - columnLower[iColumn] > tolerance) {
    3126           largest = CoinMax(largest, fabs(solution[iColumn]));
    3127         }
    3128       }
    3129       largest *= factor;
    3130     } else {
    3131       // absolute
    3132       largest = -factor;
    3133     }
    3134     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    3135       if (columnUpper[iColumn] - columnLower[iColumn] > tolerance) {
    3136         newUpper[iColumn] = CoinMin(columnUpper[iColumn], largest);
    3137         newLower[iColumn] = CoinMax(columnLower[iColumn], -largest);
    3138       }
    3139     }
    3140   }
    31413471  int numberRows = model.getNumRows();
    31423472  const double *rowLower = model.getRowLower();
    31433473  const double *rowUpper = model.getRowUpper();
     3474#if 1
     3475  // TEMP COPIES
     3476  {
     3477    double * cLower = new double [2*numberRows+2*numberColumns];
     3478    double * cUpper = cLower+numberColumns;
     3479    double * rLower = cUpper+numberColumns;
     3480    double * rUpper = rLower+numberRows;
     3481    memcpy(cLower, colLower, numberColumns * sizeof(double));
     3482    memcpy(cUpper, colUpper, numberColumns * sizeof(double));
     3483    memcpy(rLower, rowLower, numberRows * sizeof(double));
     3484    memcpy(rUpper, rowUpper, numberRows * sizeof(double));
     3485    char * intVar = new char [numberColumns];
     3486    for (int i=0;i<numberColumns;i++) {
     3487      if (model.isInteger(i))
     3488        intVar[i]=1;
     3489      else
     3490        intVar[i]=0;
     3491    }
     3492    CoinBigIndex * rowStartPos = new CoinBigIndex[numberRows];
     3493    // sort rows and set up rowStartPos
     3494    for (int iRow=0;iRow<numberRows;iRow++) {
     3495      CoinBigIndex rStart = rowStart[iRow];
     3496      CoinBigIndex rEnd = rStart + rowLength[iRow];
     3497      CoinSort_2(element+rStart,element+rEnd,column+rStart);
     3498      CoinBigIndex j;
     3499      for (j=rStart;j<rEnd;j++) {
     3500        if (element[j]>0)
     3501          break;
     3502      }
     3503      rowStartPos[iRow]=j;
     3504    }
     3505    int nInfeas = tighten(cLower, cUpper, column, element,
     3506                          rowStart, rowStartPos, rowLength, rLower, rUpper,
     3507                          numberRows, numberColumns, intVar,20,tolerance);
     3508    // see if any free rows
     3509    for (int iRow=0;iRow<numberRows;iRow++) {
     3510      if (rLower[iRow]<-1.0e20&&rUpper[iRow]>1.0e20) {
     3511        model.setRowLower(iRow,-COIN_DBL_MAX);
     3512        model.setRowUpper(iRow,COIN_DBL_MAX);
     3513      }
     3514    }   
     3515    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     3516      if (cLower[iColumn]>newLower[iColumn]) {
     3517        newLower[iColumn] = cLower[iColumn];
     3518        columnLower[iColumn] = cLower[iColumn];
     3519        model.setColLower(iColumn,cLower[iColumn]);
     3520      }
     3521      if (cUpper[iColumn]<newUpper[iColumn]) {
     3522        newUpper[iColumn] = cUpper[iColumn];
     3523        columnUpper[iColumn] = cUpper[iColumn];
     3524        model.setColUpper(iColumn,cUpper[iColumn]);
     3525      }
     3526    }   
     3527    delete [] rowStartPos;
     3528    delete [] intVar;
     3529    delete [] cLower;
     3530  }
     3531#endif
    31443532#ifndef NDEBUG
    31453533  double large2 = 1.0e10 * large;
     
    31723560        CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
    31733561        CoinBigIndex j;
     3562#ifdef CBC_PREPROCESS_EXPERIMENT
     3563        bool allInteger = true;
     3564#endif
    31743565        // Compute possible lower and upper ranges
    31753566
     
    31773568          double value = element[j];
    31783569          iColumn = column[j];
     3570#ifdef CBC_PREPROCESS_EXPERIMENT
     3571          if (!model.isInteger(iColumn) || value != floor(value+0.5))
     3572            allInteger = false;
     3573#endif
    31793574          if (value > 0.0) {
    31803575            if (newUpper[iColumn] >= large) {
     
    32173612          double lower = rowLower[iRow];
    32183613          double upper = rowUpper[iRow];
     3614#ifdef CBC_PREPROCESS_EXPERIMENT
     3615          if (tightenRowBounds && upper>lower) {
     3616            double lowerNew = lower;
     3617            double upperNew = upper;
     3618            maxUp = CoinMax(maxUp,lower);
     3619            maxDown = CoinMin(maxDown,upper);
     3620            if (!infiniteLower && maxDown > lower + 1.0e-6)
     3621              lowerNew = CoinMax(maxDown-1.0e-6,lower);
     3622            if (!infiniteUpper && maxUp < upper - 1.0e-6)
     3623              upperNew = CoinMin(maxUp+1.0e-6,upper);
     3624            if (lowerNew > upperNew) {
     3625              printf("BAD bounds of %g (%g) and %g (%g) on row %d\n",
     3626                     lowerNew,lower,upperNew,upper,iRow);
     3627              lowerNew = upperNew;
     3628            }
     3629            if (lower!=lowerNew||upper!=upperNew) {
     3630              if (allInteger) {
     3631                lower = floor(lower + 0.5);
     3632                upper = floor(upper + 0.5);
     3633              } else {
     3634                lower = lowerNew;
     3635                upper = upperNew;
     3636              }
     3637              if (lower!=rowLower[iRow]||upper!=rowUpper[iRow])
     3638#ifdef LOTS_OF_PRINTING
     3639                printf("On row %d bounds %g,%g -> %g,%g\n",
     3640                       iRow,rowLower[iRow],rowUpper[iRow],
     3641                       lower,upper);
     3642#endif
     3643              model.setRowLower(iRow,lower);
     3644              model.setRowUpper(iRow,upper);
     3645            }
     3646          }
     3647#endif
    32193648          for (j = rStart; j < rEnd; ++j) {
    32203649            double value = element[j];
     
    32223651            double nowLower = newLower[iColumn];
    32233652            double nowUpper = newUpper[iColumn];
     3653            bool anyChange = false;
    32243654            if (value > 0.0) {
    32253655              // positive value
     
    32423672                  // Tighten the lower bound
    32433673                  newLower[iColumn] = newBound;
    3244                   numberChanged++;
     3674                  anyChange = true;
    32453675                  // check infeasible (relaxed)
    32463676                  if (nowUpper - newBound < -100.0 * tolerance) {
     
    32773707                  // Tighten the upper bound
    32783708                  newUpper[iColumn] = newBound;
    3279                   numberChanged++;
     3709                  anyChange = true;
    32803710                  // check infeasible (relaxed)
    32813711                  if (newBound - nowLower < -100.0 * tolerance) {
     
    33143744                  // Tighten the upper bound
    33153745                  newUpper[iColumn] = newBound;
    3316                   numberChanged++;
     3746                  anyChange = true;
    33173747                  // check infeasible (relaxed)
    33183748                  if (newBound - nowLower < -100.0 * tolerance) {
     
    33493779                  // Tighten the lower bound
    33503780                  newLower[iColumn] = newBound;
    3351                   numberChanged++;
     3781                  anyChange = true;
    33523782                  // check infeasible (relaxed)
    33533783                  if (nowUpper - newBound < -100.0 * tolerance) {
     
    33673797              }
    33683798            }
     3799            if (anyChange) {
     3800              numberChanged++;
     3801            } else if (columnLength[iColumn] == 1) {
     3802#ifdef CBC_PREPROCESS_EXPERIMENT
     3803              // may be able to do better
     3804              // should be picked up elsewhere if no objective
     3805              if (objective[iColumn]) {
     3806                double newBound;
     3807                if (direction*objective[iColumn]>0.0) {
     3808                  // want objective as low as possible so reduce upper bound
     3809                  if (value < 0.0) {
     3810                    double gap = maxUp-upper;
     3811                    newBound = newLower[iColumn] - gap/value;
     3812                  } else {
     3813                    double gap = lower-maxDown;
     3814                    newBound = newLower[iColumn] + gap/value;
     3815                  }
     3816                  if (newBound<newUpper[iColumn]-1.0e-7) {
     3817                    if (model.isInteger(iColumn)) {
     3818                      newBound = ceil(newBound);
     3819                      newBound = CoinMax(newLower[iColumn],newBound);
     3820                    } else {
     3821                      newBound = CoinMax(newLower[iColumn],newBound+1.0e-7);
     3822                    }
     3823                  }
     3824                  if (newBound<newUpper[iColumn]-1.0e-7) {
     3825                    numberChanged++;
     3826#ifdef LOTS_OF_PRINTING
     3827                    printf("singleton %d obj %g %g <= %g rlo %g rup %g maxd %g maxu %g el %g\n",
     3828                           iColumn,objective[iColumn],newLower[iColumn],newUpper[iColumn],
     3829                           lower,upper,maxDown,maxUp,value);
     3830                    printf("upperbound changed to %g\n",newBound);
     3831#endif
     3832                    newUpper[iColumn] = newBound;
     3833                    anyChange = true;
     3834                    // check infeasible (relaxed)
     3835                    if (newBound - nowLower < -100.0 * tolerance) {
     3836                      numberInfeasible++;
     3837                    }
     3838                    // adjust
     3839                    double now;
     3840                    if (nowUpper > large) {
     3841                      now = 0.0;
     3842                      infiniteUpper--;
     3843                    } else {
     3844                      now = nowUpper;
     3845                    }
     3846                    maximumUp += (newBound - now) * value;
     3847                    maxUp += (newBound - now) * value;
     3848                    nowUpper = newBound;
     3849                  }
     3850                } else {
     3851                  // want objective as low as possible so increase lower bound
     3852                  if (value > 0.0) { // ?
     3853                    double gap = maxUp-upper;
     3854                    newBound = newUpper[iColumn] - gap/value;
     3855                  } else {
     3856                    double gap = lower-maxDown;
     3857                    newBound = newUpper[iColumn] + gap/value;
     3858                  }
     3859                  if (newBound>newLower[iColumn]+1.0e-7) {
     3860                    if (model.isInteger(iColumn)) {
     3861                      newBound = ceil(newBound);
     3862                      newBound = CoinMin(newUpper[iColumn],newBound);
     3863                    } else {
     3864                      newBound = CoinMin(newUpper[iColumn],newBound+1.0e-7);
     3865                    }
     3866                  }
     3867                  if (newBound>newLower[iColumn]+1.0e-7) {
     3868                    numberChanged++;
     3869#ifdef LOTS_OF_PRINTING
     3870                    printf("singleton %d obj %g %g <= %g rlo %g rup %g maxd %g maxu %g el %g\n",
     3871                           iColumn,objective[iColumn],newLower[iColumn],newUpper[iColumn],
     3872                           lower,upper,maxDown,maxUp,value);
     3873                    printf("lowerbound changed to %g\n",newBound);
     3874#endif
     3875                    newLower[iColumn] = newBound;
     3876                    anyChange = true;
     3877                    // check infeasible (relaxed)
     3878                    if (nowUpper - newBound < -100.0 * tolerance) {
     3879                      numberInfeasible++;
     3880                    }
     3881                    // adjust
     3882                    double now;
     3883                    if (nowLower < -large) {
     3884                      now = 0.0;
     3885                      infiniteLower--;
     3886                    } else {
     3887                      now = nowLower;
     3888                    }
     3889                    maximumDown += (newBound - now) * value;
     3890                    maxDown += (newBound - now) * value;
     3891                    nowLower = newBound;
     3892                  }
     3893                }
     3894              }
     3895#endif
     3896            }
    33693897          }
    33703898        }
     
    35794107  }
    35804108  if (numberInfeasible) {
    3581     handler_->message(CGL_INFEASIBLE, messages_)
    3582       << CoinMessageEol;
    35834109    // restore column bounds
    35844110    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     
    35924118  delete[] columnUpper;
    35934119  return (numberInfeasible);
     4120}
     4121#define CGL_REASONABLE_INTEGER_BOUND 1.23456789e10
     4122// This tightens column bounds (and can declare infeasibility)
     4123// It may also declare rows to be redundant
     4124static int
     4125tighten(double *colLower, double * colUpper,
     4126        const int *column, const double *rowElements,
     4127        const CoinBigIndex *rowStart,
     4128        const CoinBigIndex * rowStartPos,const int * rowLength,
     4129        double *rowLower, double *rowUpper,
     4130        int nRows,int nCols,char * intVar,int maxpass,
     4131        double tolerance)
     4132{
     4133  int i, j;
     4134  CoinBigIndex k, krs, kre;
     4135  int dolrows;
     4136  int iflagu, iflagl;
     4137  int ntotal=0,nchange=1,jpass=0;
     4138  double dmaxup, dmaxdown, dbound;
     4139  int ninfeas=0;
     4140  // Later try with cliques????
     4141  assert (rowStartPos);
     4142  while(nchange) {
     4143    nchange = 0;
     4144    if (jpass==maxpass) break;
     4145    jpass++;
     4146    dolrows = (jpass & 1) == 1;
     4147   
     4148    for (i = 0; i < nRows; ++i) {
     4149      if (rowLower[i]>-1.0e20||rowUpper[i]<1.0e20) {
     4150        int iflagu = 0;
     4151        int iflagl = 0;
     4152        double dmaxup = 0.0;
     4153        double dmaxdown = 0.0;
     4154        CoinBigIndex krs = rowStart[i];
     4155        CoinBigIndex krs2 = rowStartPos[i];
     4156        CoinBigIndex kre = rowStart[i]+rowLength[i];
     4157       
     4158        /* ------------------------------------------------------------*/
     4159        /* Compute L(i) and U(i) */
     4160        /* ------------------------------------------------------------*/
     4161        for (k = krs; k < krs2; ++k) {
     4162          double value=rowElements[k];
     4163          int j = column[k];
     4164          if (colUpper[j] < 1.0e12)
     4165            dmaxdown += colUpper[j] * value;
     4166          else
     4167            ++iflagl;
     4168          if (colLower[j] > -1.0e12)
     4169            dmaxup += colLower[j] * value;
     4170          else
     4171            ++iflagu;
     4172        }
     4173        for (k = krs2; k < kre; ++k) {
     4174          double value=rowElements[k];
     4175          int j = column[k];
     4176          if (colUpper[j] < 1.0e12)
     4177            dmaxup += colUpper[j] * value;
     4178          else
     4179            ++iflagu;
     4180          if (colLower[j] > -1.0e12)
     4181            dmaxdown += colLower[j] * value;
     4182          else
     4183            ++iflagl;
     4184        }
     4185        if (iflagu)
     4186          dmaxup=1.0e31;
     4187        if (iflagl)
     4188          dmaxdown=-1.0e31;
     4189        if (dmaxup <= rowUpper[i] + tolerance && dmaxdown >= rowLower[i] - tolerance) {
     4190          /*
     4191           * The sum of the column maxs is at most the row ub, and
     4192           * the sum of the column mins is at least the row lb;
     4193           * this row says nothing at all.
     4194           * I suspect that this corresponds to
     4195           * an implied column singleton in the paper (viii, on p. 325),
     4196           * where the singleton in question is the row slack.
     4197           */
     4198          ++nchange;
     4199          rowLower[i]=-COIN_DBL_MAX;
     4200          rowUpper[i]=COIN_DBL_MAX;
     4201        } else {
     4202          if (dmaxup < rowLower[i] -tolerance || dmaxdown > rowUpper[i]+tolerance) {
     4203            ninfeas++;
     4204            break;
     4205          }
     4206          /*        Finite U(i) */
     4207          /* -------------------------------------------------------------*/
     4208          /* below is deliberate mistake (previously was by chance) */
     4209          /*        never do both */
     4210          if (iflagu == 0 && rowLower[i] > 0.0 && iflagl == 0 && rowUpper[i] < 1e15) {
     4211            if (dolrows) {
     4212              iflagu = 1;
     4213            } else {
     4214              iflagl = 1;
     4215            }
     4216          }
     4217          if (iflagu == 0 && rowLower[i] > -1e15) {
     4218            for (k = krs; k < kre; ++k) {
     4219              double value=rowElements[k];
     4220              j = column[k];
     4221              if (value > 0.0) {
     4222                if (colUpper[j] < 1.0e12) {
     4223                  dbound = colUpper[j] + (rowLower[i] - dmaxup) / value;
     4224                  if (dbound > colLower[j] + 1.0e-8) {
     4225                    /* we can tighten the lower bound */
     4226                    /* the paper mentions this as a possibility on p. 227 */
     4227                    colLower[j] = dbound;
     4228                    ++nchange;
     4229                   
     4230                    /* this may have fixed the variable */
     4231                    /* I believe that this roughly corresponds to a
     4232                     * forcing constraint in the paper (p. 226).
     4233                     * If there is a forcing constraint (with respect
     4234                     * to the original, untightened bounds), then in this
     4235                     * loop we will go through all the columns and fix
     4236                     * each of them to their implied bound, rather than
     4237                     * determining that the row as a whole is forced
     4238                     * and just fixing them without doing computation for
     4239                     * each column (as in the paper).
     4240                     * By doing it this way, we can tighten bounds and
     4241                     * get forcing constraints for free.
     4242                     */
     4243                    if (colUpper[j] - colLower[j] <= tolerance) {
     4244                      /* --------------------------------------------------*/
     4245                      /*                check if infeasible !!!!! */
     4246                      /* --------------------------------------------------*/
     4247                      if (colUpper[j] - colLower[j] < -100.0*tolerance) {
     4248                        ninfeas++;
     4249                      }
     4250                    }
     4251                  }
     4252                }
     4253              } else {
     4254                if (colLower[j] > -1.0e12) {
     4255                  dbound = colLower[j] + (rowLower[i] - dmaxup) / value;
     4256                  if (dbound < colUpper[j] - 1.0e-8) {
     4257                    colUpper[j] = dbound;
     4258                    ++nchange;
     4259                    if (colUpper[j] - colLower[j] <= tolerance) {
     4260                      /* --------------------------------------------------*/
     4261                      /*                check if infeasible !!!!! */
     4262                      /* --------------------------------------------------*/
     4263                      if (colUpper[j] - colLower[j] < -100.0*tolerance) {
     4264                        ninfeas++;
     4265                      }
     4266                    }
     4267                  }
     4268                }
     4269              }
     4270            }
     4271          }
     4272         
     4273          /* ----------------------------------------------------------------*/
     4274          /*        Finite L(i) */
     4275          /* ----------------------------------------------------------------*/
     4276          if (iflagl == 0 && rowUpper[i] < 1e15) {
     4277            for (k = krs; k < kre; ++k) {
     4278              double value=rowElements[k];
     4279              j = column[k];
     4280              if (value < 0.0) {
     4281                if (colUpper[j] < 1.0e12) {
     4282                  dbound = colUpper[j] + (rowUpper[i] - dmaxdown) / value;
     4283                  if (dbound > colLower[j] + 1.0e-8) {
     4284                    colLower[j] = dbound;
     4285                    ++nchange;
     4286                    if (! (colUpper[j] - colLower[j] > tolerance)) {
     4287                      /* --------------------------------------------------*/
     4288                      /*                check if infeasible !!!!! */
     4289                      /* --------------------------------------------------*/
     4290                      if (colUpper[j] - colLower[j] < -100.0*tolerance) {
     4291                        ninfeas++;
     4292                      }
     4293                    }
     4294                  }
     4295                }
     4296              } else {
     4297                if (colLower[j] > -1.0e12) {
     4298                  dbound = colLower[j] + (rowUpper[i] - dmaxdown) / value;
     4299                  if (dbound < colUpper[j] - 1.0e-8) {
     4300                    colUpper[j] = dbound;
     4301                    ++nchange;
     4302                    if (! (colUpper[j] - colLower[j] > tolerance)) {
     4303                      /* --------------------------------------------------*/
     4304                      /*                check if infeasible !!!!! */
     4305                      /* --------------------------------------------------*/
     4306                      if (colUpper[j] - colLower[j] < -100.0*tolerance) {
     4307                        ninfeas++;
     4308                      }
     4309                    }
     4310                  }
     4311                }
     4312              }
     4313            }
     4314          }
     4315        }
     4316      }
     4317    }
     4318    for (j = 0; j < nCols; ++j) {
     4319      if (intVar[j]) {
     4320        if (colUpper[j]-colLower[j]>1.0e-8) {
     4321          if (floor(colUpper[j]+1.0e-4)<colUpper[j])
     4322            nchange++;
     4323          // clean up anyway
     4324          colUpper[j]=floor(colUpper[j]+1.0e-4);
     4325          if (ceil(colLower[j]-1.0e-4)>colLower[j])
     4326            nchange++;
     4327          // clean up anyway
     4328          colLower[j]=ceil(colLower[j]-1.0e-4);
     4329          if (colUpper[j]<colLower[j]) {
     4330            /*printf("infeasible\n");*/
     4331            ninfeas++;
     4332          }
     4333        } else {
     4334          // clean
     4335          colUpper[j]=floor(colUpper[j]+1.0e-4);
     4336          colLower[j]=ceil(colLower[j]-1.0e-4);
     4337          if (colUpper[j]<colLower[j]) {
     4338            /*printf("infeasible\n");*/
     4339            ninfeas++;
     4340          }
     4341        }
     4342      }
     4343    }
     4344    if (ninfeas) break;
     4345  }
     4346  return (ninfeas);
    35944347}
    35954348
  • trunk/src/CglPreProcess/CglPreProcess.hpp

    r1526 r1572  
    7171      fixed or integral, bounds are slightly looser than they could be.
    7272      Returns non-zero if problem infeasible
    73       Fudge for branch and bound - put bounds on columns of factor *
    74       largest value (at continuous) - should improve stability
    75       in branch and bound on infeasible branches (0.0 is off)
    76   */
    77   int tightenPrimalBounds(OsiSolverInterface &model, double factor = 0.0);
     73  */
     74  int tightenPrimalBounds(OsiSolverInterface &model,
     75                          bool tightenRowBounds=false);
    7876  /** Fix some of problem - returning new problem.
    7977      Uses reduced costs.
  • trunk/src/CglProbing/CglProbing.cpp

    r1567 r1572  
    16851685  bool feasible=true;
    16861686  if (!info->inTree&&!info->pass) {
     1687#if CBC_PREPROCESS_EXPERIMENT>1
     1688    int nRows = si.getNumRows();
     1689    double * columnLo = new double [2*nRows+2*nCols];
     1690    double * columnUp = columnLo + nCols;
     1691    double * rowLo = columnUp + nCols;
     1692    double * rowUp = rowLo + nRows;
     1693    // temp
     1694    memcpy(columnLo,si.getColLower(),nCols*sizeof(double));
     1695    memcpy(columnUp,si.getColUpper(),nCols*sizeof(double));
     1696    memcpy(rowLo,si.getRowLower(),nRows*sizeof(double));
     1697    memcpy(rowUp,si.getRowUpper(),nRows*sizeof(double));
     1698    // Get a row copy in standard format
     1699    CoinPackedMatrix rowCopy = *si.getMatrixByRow();
     1700    int columnsTightened;
     1701    int rowsTightened;
     1702    int elementsChanged;
     1703    si.tightPrimalBounds(rowLo,rowUp,columnLo,columnUp,
     1704                         columnsTightened,rowsTightened,
     1705                         elementsChanged,rowCopy);
     1706    int nChanged = 0;
     1707    CoinPackedVector lbs;
     1708    CoinPackedVector ubs;
     1709    for (int iColumn = 0;iColumn<nCols;iColumn++) {
     1710      if (columnLo[iColumn]>colLower[iColumn]) {
     1711        nChanged++;
     1712        colLower[iColumn] = columnLo[iColumn];
     1713        if (intVar[iColumn])
     1714          lbs.insert(iColumn,columnLo[iColumn]);
     1715      }
     1716      if (columnUp[iColumn]<colUpper[iColumn]) {
     1717        nChanged++;
     1718        colUpper[iColumn] = columnUp[iColumn];
     1719        if (intVar[iColumn])
     1720          ubs.insert(iColumn,columnUp[iColumn]);
     1721      }
     1722    }
     1723#ifdef LOTS_OF_PRINTING
     1724    printf("%d bounds changed, %d lower int, %d upper int\n",
     1725           nChanged,lbs.getNumElements(),ubs.getNumElements());
     1726#endif
     1727    if (lbs.getNumElements()) {
     1728      OsiColCut cc;
     1729      cc.setLbs(lbs);
     1730      cs.insert(cc);
     1731    }
     1732    if (ubs.getNumElements()) {
     1733      OsiColCut cc;
     1734      cc.setUbs(ubs);
     1735      cs.insert(cc);
     1736    }
     1737    if (rowCopy_) {
     1738      //printf("Can't tighten row bounds\n");
     1739    } else {
     1740      nChanged=0;
     1741      for (int iRow = 0;iRow<nRows;iRow++) {
     1742        if (rowLo[iRow]>rowLower[iRow]) {
     1743          nChanged++;
     1744          rowLower[iRow] = rowLo[iRow];
     1745        }
     1746        if (rowUp[iRow]<rowUpper[iRow]) {
     1747          nChanged++;
     1748          rowUpper[iRow] = rowUp[iRow];
     1749        }
     1750      }
     1751#ifdef LOTS_OF_PRINTING
     1752      printf("%d row bounds changed\n",nChanged);
     1753#endif
     1754    }
     1755#endif
    16871756    // make more integer
    16881757    feasible = analyze(&si,intVar,colLower,colUpper);
Note: See TracChangeset for help on using the changeset viewer.