Changeset 964


Ignore:
Timestamp:
Feb 10, 2011 10:33:31 PM (11 years ago)
Author:
lou
Message:

Two-cuts pulled out as implicationsCuts. Initial construction and grooming of
row copy pulled to setupRowCopy, groomModel. Rewrote CglProbingRowCut?, better
structure, excised OsiRowCut2 (folded whichRow into OsiRowCut?). Continue to
disentangle control variables. Bug fixes. Does not die on any miplib3 problem.
Slooow on some.

Location:
branches/CglWorking101215/src/CglProbing
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • branches/CglWorking101215/src/CglProbing/CglProbing.cpp

    r963 r964  
    2424#include "CglProbingDebug.hpp"
    2525
    26 /*
    27   PROBING_EXTRA_STUFF has been false since June 2006, so anything guarded by
    28   this symbol is likely to be extremely dusty.
    29 */
    30 //#define PROBING_EXTRA_STUFF true
    31 #define PROBING_EXTRA_STUFF false
    32 // Seems to be unused  -- lh, 101202 --
    33 // #define FIXED_ALLOWANCE 10
    34 #define SIZE_ROW_MULT 4
    35 #define SIZE_ROW_ADD 2000
    3626
    3727typedef struct {double infeasibility;int sequence;} double_int_pair;
     
    4939#endif
    5040
     41namespace {
     42/*
     43  This method will attempt to identify variables that will naturally take on
     44  integer values but were not originally identified as such, and mark them as
     45  integer. Upper and lower bounds are forced to integrality.
     46 
     47  Executed once, at first cut generation pass at the root.
     48
     49  Which begs various questions:
     50    * Why not execute this independently, before any cut generation occurs?
     51    * Why not execute repeatedly, as the problem is simplified due to fixed
     52      variables?
     53
     54  The expression abs(x - floor(x+.5)) evaluates to zero when x is integer and
     55  is nonzero otherwise. Farther down in the code, there's an odd-looking test
     56  that checks abs(1/aij - floor(1/aij + .5)). What we're looking for is
     57  coefficients aij = 1/k, k integer. This guarantees that b/aij is integer
     58  as long as b is integer.
     59*/
     60
     61bool analyze (const OsiSolverInterface *solverX, char *intVar,
     62              double *lower, double *upper)
     63{
     64  const double e20Inf = 1.0e20 ;
     65  const int changedToInt = 77 ;
     66
     67  /*
     68    I don't see why we clone the solver here. There are no modifications to
     69    anything except the parameters intVar, lower, and upper.
     70  */
     71  OsiSolverInterface *solver = solverX->clone() ;
     72  const double *objective = solver->getObjCoefficients() ;
     73  int numberColumns = solver->getNumCols() ;
     74  int numberRows = solver->getNumRows() ;
     75  double direction = solver->getObjSense() ;
     76  int iRow,iColumn ;
     77
     78  /*
     79    Nor do I understand why we make copies of the matrix. Two, in fact! Surely
     80    we could ask the OSI for these. Aaaaah, but then we might be triggering
     81    changes inside the OSI? Shouldn't matter, getMatrixBy* are const.
     82  */
     83  // Row copy
     84  CoinPackedMatrix matrixByRow(*solver->getMatrixByRow()) ;
     85  const double *elementByRow = matrixByRow.getElements() ;
     86  const int *column = matrixByRow.getIndices() ;
     87  const CoinBigIndex *rowStart = matrixByRow.getVectorStarts() ;
     88  const int *rowLength = matrixByRow.getVectorLengths() ;
     89
     90  // Column copy
     91  CoinPackedMatrix matrixByCol(*solver->getMatrixByCol()) ;
     92  const double *element = matrixByCol.getElements() ;
     93  const int *row = matrixByCol.getIndices() ;
     94  const CoinBigIndex *columnStart = matrixByCol.getVectorStarts() ;
     95  const int *columnLength = matrixByCol.getVectorLengths() ;
     96
     97  const double * rowLower = solver->getRowLower() ;
     98  const double * rowUpper = solver->getRowUpper() ;
     99
     100  char *ignore = new char [numberRows] ;
     101  int *which = new int[numberRows] ;
     102  double *changeRhs = new double[numberRows] ;
     103  memset(changeRhs,0,numberRows*sizeof(double)) ;
     104  memset(ignore,0,numberRows) ;
     105
     106  int numberChanged = 0 ;
     107  bool finished = false ;
     108/*
     109  Open main analysis loop. Keep going until nothing changes.
     110*/
     111  while (!finished) {
     112    int saveNumberChanged = numberChanged ;
     113/*
     114  Open loop to scan each constraint. With luck, we can conclude that some
     115  variable must be integer. With less luck, the constraint will not prevent the
     116  variable from being integer.
     117*/
     118    for (iRow = 0 ; iRow < numberRows ; iRow++) {
     119      int numberContinuous = 0 ;
     120      double value1 = 0.0,value2 = 0.0 ;
     121      bool allIntegerCoeff = true ;
     122      double sumFixed = 0.0 ;
     123      int jColumn1 = -1,jColumn2 = -1 ;
     124/*
     125  Scan the coefficients of the constraint and accumulate some information:
     126    * Contribution of fixed variables.
     127    * Count of continuous variables, and column index and coefficients for
     128      the first and last continuous variable.
     129    * A boolean, allIntegerCoeff, true if all coefficients of unfixed
     130      integer variables are integer.
     131*/
     132      for (CoinBigIndex j = rowStart[iRow] ;
     133           j < rowStart[iRow]+rowLength[iRow] ; j++) {
     134        int jColumn = column[j] ;
     135        double value = elementByRow[j] ;
     136        if (upper[jColumn] > lower[jColumn]+1.0e-8) {
     137          if (!intVar[jColumn]) {
     138            if (numberContinuous == 0) {
     139              jColumn1 = jColumn ;
     140              value1 = value ;
     141            } else {
     142              jColumn2 = jColumn ;
     143              value2 = value ;
     144            }
     145            numberContinuous++ ;
     146          } else {
     147            if (fabs(value-floor(value+0.5)) > 1.0e-12)
     148              allIntegerCoeff = false ;
     149          }
     150        } else {
     151          sumFixed += lower[jColumn]*value ;
     152        }
     153      }
     154/*
     155  See if the row bounds are integer after adjusting for the contribution of
     156  fixed variables. Serendipitous cancellation is possible here (if unlikely).
     157  The fixed contribution could change a non-integer bound to an integer.
     158*/
     159      double low = rowLower[iRow] ;
     160      if (low > -e20Inf) {
     161        low -= sumFixed ;
     162        if (fabs(low-floor(low+0.5)) > 1.0e-12)
     163          allIntegerCoeff = false ;
     164      }
     165      double up = rowUpper[iRow] ;
     166      if (up<e20Inf) {
     167        up -= sumFixed ;
     168        if (fabs(up-floor(up+0.5)) > 1.0e-12)
     169          allIntegerCoeff = false ;
     170      }
     171/*
     172  To make progress, it must be true that the coefficients of all unfixed
     173  integer variables are integer and the row bounds are integer.
     174*/
     175      if (!allIntegerCoeff) continue ;
     176/*
     177  Case 1: We have a single continuous variable to consider.
     178
     179  If a<ij> is of the form 1/k, k integer, then x<j> can be integer. Put another
     180  way, this constraint will not prevent integrality and we can ignore it in
     181  further processing.  An equality forces integrality.
     182
     183  So, for an equality with one continuous variable, shouldn't we mark the
     184  variable as `cannot be declared integer' so we don't waste time in subsequent
     185  passes?
     186*/
     187      if (numberContinuous == 1) {
     188        if (low == up) {
     189          if (fabs(value1) > 1.0e-3) {
     190            value1 = 1.0/value1 ;
     191            if (fabs(value1-floor(value1+0.5)) < 1.0e-12) {
     192              numberChanged++ ;
     193              intVar[jColumn1] = changedToInt ;
     194            }
     195          }
     196        } else {
     197          if (fabs(value1) > 1.0e-3) {
     198            value1 = 1.0/value1 ;
     199            if (fabs(value1-floor(value1+0.5)) < 1.0e-12) {
     200              ignore[iRow] = 1 ;
     201            }
     202          }
     203        }
     204      } else
     205/*
     206  Case 2: Two continuous variables.
     207
     208  John's comment is `need general theory - for now just look at 2 cases'. I've
     209  worked out the general theory. See the paper documentation. The code needs
     210  to be rewritten. What's here is overly restrictive for case 2.1 and wrong for
     211  case 2.2.
     212
     213  2.1) Constraint is (integer lhs) + x<1> - x<2> = (integer rhs), x<1> and x<2>
     214     have lower bounds of zero and do not appear in any other constraint.
     215     The overall contribution to the objective, (c<1>+c<2>), is unfavourable.
     216     Consider any given value of (rhs-lhs) = (integer). Some x<i> will head
     217     towards infinity, the other towards zero, modulo finite bounds and
     218     feasibility. See the paper documentation.
     219
     220  2.2) Constraint as above but two coefficients in each column. The second
     221     coefficients feed into G/L row(s) which will try and minimize. Also
     222     consider that the second coefficients may be in the same row. Code as
     223     written is incorrect.
     224
     225  Case 2.2) has been disabled from the start (with a note that says `take out
     226  until fixed').
     227*/
     228      if (numberContinuous == 2) {
     229        if (low == up) {
     230          if (fabs(value1) == 1.0 && value1*value2 == -1.0 &&
     231              !lower[jColumn1] && !lower[jColumn2] &&
     232              columnLength[jColumn1] == 1 && columnLength[jColumn2] == 1) {
     233            int n = 0 ;
     234            int i ;
     235            double objChange =
     236                direction*(objective[jColumn1]+objective[jColumn2]) ;
     237            double bound = CoinMin(upper[jColumn1],upper[jColumn2]) ;
     238            bound = CoinMin(bound,e20Inf) ;
     239            // Since column lengths are 1, there is no second coeff.
     240            for (i = columnStart[jColumn1] ;
     241                 i < columnStart[jColumn1]+columnLength[jColumn1] ; i++) {
     242              int jRow = row[i] ;
     243              double value = element[i] ;
     244              if (jRow != iRow) {
     245                which[n++] = jRow ;
     246                changeRhs[jRow] = value ;
     247              }
     248            }
     249            for (i = columnStart[jColumn2] ;
     250                 i < columnStart[jColumn2]+columnLength[jColumn2] ; i++) {
     251              int jRow = row[i] ;
     252              double value = element[i] ;
     253              if (jRow != iRow) {
     254                if (!changeRhs[jRow]) {
     255                  which[n++] = jRow ;
     256                  changeRhs[jRow] = value ;
     257                } else {
     258                  changeRhs[jRow] += value ;
     259                }
     260              }
     261            }
     262            if (objChange >= 0.0) {
     263              bool good = true ;
     264              // Since n = 0, this loop never executes
     265              for (i = 0 ; i < n ; i++) {
     266                int jRow = which[i] ;
     267                double value = changeRhs[jRow] ;
     268                if (value) {
     269                  value *= bound ;
     270                  if (rowLength[jRow] == 1) {
     271                    if (value>0.0) {
     272                      double rhs = rowLower[jRow] ;
     273                      if (rhs>0.0) {
     274                        double ratio =rhs/value ;
     275                        if (fabs(ratio-floor(ratio+0.5))>1.0e-12)
     276                          good = false ;
     277                      }
     278                    } else {
     279                      double rhs = rowUpper[jRow] ;
     280                      if (rhs<0.0) {
     281                        double ratio =rhs/value ;
     282                        if (fabs(ratio-floor(ratio+0.5))>1.0e-12)
     283                          good = false ;
     284                      }
     285                    }
     286                  } else if (rowLength[jRow] == 2) {
     287                    if (value>0.0) {
     288                      if (rowLower[jRow]>-e20Inf)
     289                        good = false ;
     290                    } else {
     291                      if (rowUpper[jRow]<e20Inf)
     292                        good = false ;
     293                    }
     294                  } else {
     295                    good = false ;
     296                  }
     297                }
     298              }
     299
     300              if (good) {
     301                // both can be integer
     302                numberChanged++ ;
     303                intVar[jColumn1] = changedToInt ;
     304                numberChanged++ ;
     305                intVar[jColumn2] = changedToInt ;
     306              }
     307            }
     308            // clear (another noop since n is zero)
     309            for (i = 0 ; i < n ; i++) {
     310              changeRhs[which[i]] = 0.0 ;
     311            }
     312          }
     313        }
     314      }
     315    } // end loop to scan each constraint
     316/*
     317  Check for variables x<j> with integer bounds l<j> and u<j> and no
     318  constraint that prevents x<j> from being integer.
     319*/
     320    for (iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
     321      if (upper[iColumn] > lower[iColumn]+1.0e-8 && !intVar[iColumn]) {
     322        double value ;
     323        value = upper[iColumn] ;
     324        if (value < e20Inf && fabs(value-floor(value+0.5)) > 1.0e-12)
     325          continue ;
     326        value = lower[iColumn] ;
     327        if (value > -e20Inf && fabs(value-floor(value+0.5)) > 1.0e-12)
     328          continue ;
     329        bool integer = true ;
     330        for (CoinBigIndex j = columnStart[iColumn] ;
     331             j < columnStart[iColumn]+columnLength[iColumn] ; j++) {
     332          int iRow = row[j] ;
     333          if (!ignore[iRow]) {
     334            integer = false ;
     335            break ;
     336          }
     337        }
     338        if (integer) {
     339          numberChanged++ ;
     340          intVar[iColumn] = changedToInt ;
     341        }
     342      }
     343    }
     344    finished = (numberChanged == saveNumberChanged) ;
     345  } // end main loop: while !finished
     346/*
     347  Nothing changed on the last pass. Time to clean up the new integer
     348  variables.  Force bounds to integer values. Feasibility can be lost here.
     349  The codes assigned here are curious choices. If the variable is fixed at
     350  zero, declare it binary (1), but if fixed at some other value, declare it
     351  continuous (0). Otherwise, declare it general integer (2).
     352*/
     353  bool feasible = true ;
     354  for (iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
     355    if (intVar[iColumn] == changedToInt) {
     356      if (upper[iColumn] > e20Inf) {
     357        upper[iColumn] = e20Inf ;
     358      } else {
     359        upper[iColumn] = floor(upper[iColumn]+1.0e-5) ;
     360      }
     361      if (lower[iColumn] < -e20Inf) {
     362        lower[iColumn] = -e20Inf ;
     363      } else {
     364        lower[iColumn] = ceil(lower[iColumn]-1.0e-5) ;
     365        if (lower[iColumn] > upper[iColumn])
     366          feasible = false ;
     367      }
     368      if (lower[iColumn] == 0.0 && upper[iColumn] == 1.0)
     369        intVar[iColumn] = 1 ;
     370      else if (lower[iColumn] == upper[iColumn])
     371        intVar[iColumn] = 0 ;
     372      else
     373        intVar[iColumn] = 2 ;
     374    }
     375  }
     376/*
     377  Clean up and return.
     378*/
     379  delete [] which ;
     380  delete [] changeRhs ;
     381  delete [] ignore ;
     382  delete solver ;
     383
     384# if CGLPROBING_DEBUG > 0
     385  std::cout
     386    << "CglProbing::analyze: " << numberChanged
     387    << " variables could be made integer." << std::endl ;
     388# endif
     389
     390  return feasible ;
     391}
     392
     393/*
     394  Set up the working row copy.
     395 
     396  If we have a snapshot handy, most of the setup is done -- we just make
     397  a copy of the snapshot.  Otherwise, we need to create a row-major copy
     398  of the constraint matrix and the rhs and rhslow vectors.  If the client
     399  has specified mode 0 (work from snapshot), we silently correct it.
     400*/
     401
     402CoinPackedMatrix *setupRowCopy (int mode, bool useObj,
     403                                double cutoff, double offset,
     404                                const CoinPackedMatrix *rowCopy_,
     405                                int numberRows_, int numberColumns_,
     406                                const double *const rowLower_,
     407                                const double *const rowUpper_,
     408                                const OsiSolverInterface &si,
     409                                double *rowLower, double *rowUpper)
     410{
     411# if CGL_DEBUG > 0
     412  std::cout << "Entering setupRowCopy, mode " << mode << ", " ;
     413  if (!rowCopy_) std::cout << "no " ;
     414  std::cout << "row copy." << std::endl ;
     415# endif
     416
     417  int nCols = si.getNumCols() ;
     418  int nRows = -1 ;
     419  CoinPackedMatrix *rowCopy = NULL ;
     420/*
     421  Make the working copy from the model in the si. If we're using the
     422  objective, add in the objective coefficients.
     423*/
     424  if (!rowCopy_) {
     425
     426    nRows = si.getNumRows();
     427
     428    CoinMemcpyN(si.getRowLower(),nRows,rowLower) ;
     429    CoinMemcpyN(si.getRowUpper(),nRows,rowUpper) ;
     430
     431    if (!useObj) {
     432      rowCopy = new CoinPackedMatrix(*si.getMatrixByRow()) ;
     433    } else {
     434      rowCopy = new CoinPackedMatrix(*si.getMatrixByRow(),1,nCols,false) ;
     435      const double *objective = si.getObjCoefficients() ;
     436      const double objmul = si.getObjSense() ;
     437      int *columns = new int[nCols] ;
     438      double *elements = new double[nCols] ;
     439      int kk = 0 ;
     440      for (int j = 0 ; j < nCols ; j++) {
     441        if (objective[j]) {
     442          elements[kk] = objmul*objective[j] ;
     443          columns[kk++] = j ;
     444        }
     445      }
     446      rowCopy->appendRow(kk,columns,elements) ;
     447      delete[] columns ;
     448      delete[] elements ;
     449      rowLower[nRows] = -DBL_MAX ;
     450      rowUpper[nRows] = cutoff+offset ;
     451      nRows++ ;
     452    }
     453  } else {
     454/*
     455  Make the working copy from the snapshot. It must be true that the snapshot
     456  has the same number of columns as the constraint system in the solver. The
     457  asserts are sanity checks:
     458    * the column count could go bad because the user has changed the model
     459      outside of CglProbing;
     460    * rowLower and rowUpper are allocated to the number of rows in the model
     461      plus 1 (in case we want the objective), so we need to be sure there's
     462      enough space;
     463    * the final check of row counts is simple internal consistency.
     464
     465  TODO: Another good reason to hive off snapshot as a separate class.
     466        -- lh, 100917 --
     467*/
     468    assert(nCols == numberColumns_) ;
     469    assert(nRows <= si.getNumRows()+1) ;
     470    assert(rowCopy_->getNumRows() == numberRows_) ;
     471
     472    nRows = numberRows_ ;
     473    rowCopy = new CoinPackedMatrix(*rowCopy_) ;
     474    CoinMemcpyN(rowLower_,nRows,rowLower) ;
     475    CoinMemcpyN(rowUpper_,nRows,rowUpper) ;
     476    if (useObj) {
     477      rowLower[nRows-1] = -DBL_MAX ;
     478      rowUpper[nRows-1] = cutoff+offset ;
     479    }
     480  }
     481# if CGL_DEBUG > 0
     482  if (rowCopy) {
     483    std::cout << "  verifying matrix." << std::endl ;
     484    rowCopy->verifyMtx(2) ;
     485  }
     486  std::cout << "Leaving setupRowCopy." << std::endl ;
     487# endif
     488  return (rowCopy) ;
     489}
     490
     491
     492/*
     493  Remove rows with too many elements and rows that are deemed useless for
     494  further probing.  Try to strengthen rows in place, if allowed. Deal with
     495  fixed variables. We're going to do this in place by editing the internal
     496  data structures of rowCopy. realRows will become a cross-reference from
     497  row indices in the groomed model back to the original indices.
     498
     499  Strengthening in place is highly specialised. We're looking for constraints
     500  that satisfy:
     501    * All variables binary
     502    * All coefficients negative (N), except for a single positive (t), and
     503      b<i> > 0. In this case, we can convert
     504        a<it>x<t> + SUM{N}a<ik>x<k> <= b<i>
     505      to
     506        (a<it>-b<i>)x<t> + SUM{N}a<ik>x<k> <= 0
     507    * All coefficients positive (P), except for a single negative (t), and
     508      blow<i> < 0. In this case, we can convert
     509        a<it>x<t> + SUM{P}a<ik>x<k> >= blow<i>
     510      to
     511        (a<it>-blow<i>)x<t> + SUM{P}a<ik>x<k> >= 0
     512  See the typeset documentation for the full derivation.
     513
     514  The original code here went to some trouble to avoid deleting the objective
     515  for density, but it would delete it for an ineffective bound. But we've
     516  gone to some trouble back in gutsOfGenerateCuts to use the objective only
     517  if there's a reasonable cutoff. Let's not second guess the decision here.
     518
     519  While we're here, we will sort the coefficients of each row so that all
     520  negative coefficients precede all positive coefficients. Record the break
     521  points in rowStartPos.
     522
     523  An alternate code block for deleting constraints on bounds would delete
     524  any constraint where an infinite (1.0e3 :-) bound on a variable would
     525  lead to an infinite row LB or UB. This is a recipe for deleting the
     526  entire model, when run prior to an initial round of bound propagation. It's
     527  been disabled since r675 (2008). I've chopped it entirely. Something useful
     528  could be done along this line, but only after an initial round of bound
     529  propagation. Arguably, an initial round of bound propagation would help all
     530  of CglProbing and should be done in any event.
     531
     532  Returns true if we're good to proceed, false if the groomed model is not
     533  worth pursuing.
     534*/
     535
     536bool groomModel (
     537                 bool useObj, int maxRowLen,
     538                 const OsiSolverInterface &si, const char *const intVar,
     539                 CoinPackedMatrix *rowCopy,
     540                 double *const rowLower, double *const rowUpper,
     541                 const double *const colLower, const double *const colUpper,
     542                 int *&realRows, CoinBigIndex *&rowStartPos,
     543                 const CglTreeInfo *info
     544                )
     545{
     546# if CGL_DEBUG > 0
     547  std::cout << "Entering groomModel." << std::endl ;
     548# endif
     549
     550  realRows = 0 ;
     551  rowStartPos = 0 ;
     552
     553  const double weakRhs = 1.0e20 ;
     554  const double strengthenRhs = 1.0e10 ;
     555/*
     556  In the first half (row selection), only elements is modified. In the second
     557  half (fixed variable removal), all are subject to change.
     558*/
     559  CoinBigIndex *rowStart = rowCopy->getMutableVectorStarts() ;
     560  int *rowLength = rowCopy->getMutableVectorLengths();
     561  double *elements = rowCopy->getMutableElements() ;
     562  int *column = rowCopy->getMutableIndices() ;
     563/*
     564  nRealRows excludes the objective, if we're using it.
     565*/
     566  int nRows = rowCopy->getNumRows() ;
     567  const int nRealRows = (useObj)?(nRows-1):(nRows) ;
     568  const int nCols = rowCopy->getNumCols() ;
     569  CoinBigIndex nElements = rowCopy->getNumElements() ;
     570
     571/*
     572  Coefficient strengthening is allowed only on the first call of the
     573  generator, and only if the user has provided an array to record the
     574  strengthened rows.
     575*/
     576  bool allowStrengthen = (info->strengthenRow && !info->pass) ;
     577  int *which = new int [nRealRows] ;
     578
     579  int nDense = 0 ;
     580  int nWeak = 0 ;
     581/*
     582  Preliminary assessment for dense rows and free rows.  Count the number of
     583  coefficients that we'll remove. We're looking for rows that are too long or
     584  have row bounds so weak as to be ineffective.  If the total coefficients
     585  in dense rows amounts to less than 1/10 of the matrix, just deal with it
     586  (maxRowLen set to nCols won't exclude anything).
     587
     588  At the end of this markup, which[i] holds the row length, or (nCols+1) for
     589  rows we're going to delete for weak rhs.
     590*/
     591  for (int i = 0 ; i < nRealRows ; i++) {
     592    const int leni = rowLength[i] ;
     593    if ((rowLower[i] < -weakRhs) && (rowUpper[i] > weakRhs)) {
     594      nWeak += leni ;
     595      which[i] = nCols+1 ;
     596    } else {
     597      if (leni > maxRowLen) nDense += leni ;
     598      which[i] = leni ;
     599    }
     600  }
     601  if (nDense*10 < nElements) maxRowLen = nCols ;
     602
     603/*
     604  Walk the rows again. Rows that should be deleted will come up positive
     605  when we subtract maxRowLen and can be summarily dismissed.  Note that we're
     606  converting which[] as we go, so that at the end it will contain the
     607  indices of the constraints to be deleted in a block at the beginning.
     608*/
     609  int nDelete = 0 ;
     610  int nKeep = 0 ;
     611  for (int i = 0 ; i < nRealRows ; i++) {
     612    if (which[i]-maxRowLen > 0) {
     613      which[nDelete++] = i ;
     614      continue ;
     615    }
     616/*
     617  Process the row further only if we can strengthen rows in place. Avoid
     618  range constraints to avoid complications. See the comments at the head
     619  for details.
     620 
     621  Shift the rhs values here --- it has to be done for any constraint that we
     622  keep.
     623*/
     624    const double bi = rowUpper[i] ;
     625    const double blowi = rowLower[i] ;
     626    rowLower[nKeep] = blowi ;
     627    rowUpper[nKeep] = bi ;
     628    nKeep++ ;
     629
     630    if (!allowStrengthen || (blowi > -weakRhs && bi < weakRhs)) continue ;
     631/*
     632  Assess the row. If we find anything other than binary variables, we're
     633  immediately done.
     634
     635  We're lumping a<ij> = 0 in with a<ij> < 0, but arguably we shouldn't see
     636  coefficients of zero. Introducing a whole lot of checks into this loop is
     637  counterproductive; that's why we don't abort as soon as nPlus and nMinus
     638  both exceed 1.
     639*/
     640    int nPlus = 0 ;
     641    int nMinus = 0 ;
     642    const CoinBigIndex rstarti = rowStart[i] ;
     643    const int leni = rowLength[i] ;
     644    const CoinBigIndex rendi = rstarti+leni ;
     645    for (CoinBigIndex jj = rstarti ; jj < rendi ; jj++) {
     646      int j = column[jj] ;
     647      if (intVar[j] == 1) {
     648        double value = elements[jj] ;
     649        if (value > 0.0) {
     650          nPlus++ ;
     651        } else {
     652          nMinus++ ;
     653        }
     654      } else {
     655        nPlus = 2 ;
     656        nMinus = 2 ;
     657        break ;
     658      }
     659    }
     660/*
     661  If there's one positive coefficient, find it and reduce it, given a
     662  reasonable value for b<i>. If there's one negative coefficient, find it
     663  and increase it, given a reasonable value for blow<i>. Note that only one of
     664  these will execute, even if nPlus = nMinus = 1, because only one of b<i> or
     665  blow<i> will have a reasonable value. We just don't know which one.
     666
     667  Be careful with rhs values here! Remember that shift up above.
     668*/
     669    double effectiveness = 0.0 ;
     670    if (nPlus == 1 && bi > 0.0 && bi < strengthenRhs) {
     671      for (CoinBigIndex jj = rstarti ; jj < rendi ; jj++) {
     672        double value = elements[jj] ;
     673        if (value > 0.0) {
     674          elements[jj] -= bi ;
     675#         if CGL_DEBUG > 1
     676          std::cout
     677            << "      row " << i << ": decrease a<" << i << ","
     678            << column[jj] << "> from " << value << " to "
     679            << elements[jj] << "." << std::endl ;
     680#         endif
     681        }
     682        effectiveness += fabs(value) ;
     683      }
     684      rowUpper[nKeep-1] = 0.0 ;
     685    }
     686    if (nMinus == 1 && blowi < 0.0 && blowi > -strengthenRhs) {
     687      for (CoinBigIndex jj = rstarti ; jj < rendi ; jj++) {
     688        double value = elements[jj] ;
     689        if (value < 0.0) {
     690          elements[jj] -= blowi ;
     691#         if CGL_DEBUG > 1
     692          std::cout
     693            << "      row " << i << ": increase a<" << i << ","
     694            << column[jj] << "> from " << value << " to "
     695            << elements[jj] << "." << std::endl ;
     696#         endif
     697        }
     698        effectiveness += fabs(value) ;
     699      }
     700      rowLower[nKeep-1] = 0.0 ;
     701    }
     702/*
     703  If we actually strengthened a coefficient, stash the strengthened row in
     704  strengthenRow for return to the caller.
     705*/
     706    if (effectiveness) {
     707      OsiRowCut *rc = new OsiRowCut() ;
     708      rc->setLb(rowLower[nKeep-1]) ;
     709      rc->setUb(rowUpper[nKeep-1]) ;
     710      rc->setRow(leni,column+rstarti,elements+rstarti,false) ;
     711      rc->setEffectiveness(effectiveness) ;
     712      assert (!info->strengthenRow[i]) ;
     713      info->strengthenRow[i] = rc ;
     714    }
     715  }
     716/*
     717  Deal with the objective, if we're using it. In particular, if we've deleted
     718  constraints, we need to shift the rhs entries for the objective.
     719*/
     720  if (useObj && nDelete)
     721  { rowLower[nKeep] = rowLower[nRealRows] ;
     722    rowUpper[nKeep] = rowUpper[nRealRows] ;
     723  }
     724/*
     725  Did we decide to delete anything (but not everything) in the processing
     726  loop?  If so, deal with it.  We need to physically delete the rows, of
     727  course. There may also be bookkeeping, if we're trying to strengthen rows
     728  in place.  In particular, strengthenRow will be passed back to the caller
     729  and is indexed using the original row index. We need a cross-reference,
     730  realRows, to keep track.
     731
     732  realRows is a tradeoff -- it's bigger than necessary, by nDelete entries,
     733  but nDelete will typically be small and this approach to building the
     734  xref saves some work (we need a working array of size nRows, and we want
     735  to retain the block of indices at the front of which).
     736
     737  But keep in mind that we can't strengthen the objective in place!
     738*/
     739  if (nDelete > 0 && nDelete < nRealRows) {
     740    if (info->strengthenRow) {
     741      realRows = new int [nRows] ;
     742      CoinZeroN(realRows,nRows) ;
     743      for (int k = 0 ; k < nDelete ; k++) {
     744        realRows[which[k]] = -1 ;
     745      }
     746      if (useObj) realRows[nRealRows] = -1 ;
     747      int inew = 0 ;
     748      for (int iold = 0 ; iold < nRows ; iold++) {
     749        if (!realRows[iold]) {
     750          realRows[inew++] = iold ;
     751        }
     752      }
     753      assert(inew == nKeep) ;
     754    }
     755    rowCopy->deleteRows(nDelete,which) ;
     756    nRows = rowCopy->getNumRows() ;
     757    assert((useObj && (nRows == nKeep+1)) || (!useObj && (nRows == nKeep))) ;
     758  }
     759  delete[] which ;
     760/*
     761  Do we have any real constraints left to work with? (The objective doesn't
     762  count here.) If not, clean up and go home.
     763*/
     764  if (!nKeep) {
     765#   if CGL_DEBUG > 0
     766      std::cout
     767        << "  All rows unsuitable for probing! Cleaning up." << std::endl ;
     768#   endif
     769    delete[] realRows ;
     770    return (false) ;
     771  }
     772/*
     773  So far, so good. We have something that looks like a reasonable collection
     774  of rows. Time to deal with fixed variables. We're going to walk the rows,
     775  compressing out the coefficients of the fixed variables, and we'll sort
     776  the remaining coefficients so that all negative coefficients precede
     777  all positive coefficients.  There will be no gaps when we're done.
     778
     779  NOTE that columns are not physically removed.
     780
     781  We need to refresh the mutable pointers; deleteRows above may have rendered
     782  the previous pointers invalid.
     783
     784  TODO: This code duplicates a block in snapshot(), except that the code in
     785        snapshot doesn't substitute for the values of fixed variables. The
     786        two blocks should probably be pulled out into a method.
     787*/
     788  elements = rowCopy->getMutableElements() ;
     789  column = rowCopy->getMutableIndices() ;
     790  rowStart = rowCopy->getMutableVectorStarts() ;
     791  rowLength = rowCopy->getMutableVectorLengths();
     792
     793  CoinBigIndex newSize = 0 ;
     794  int *columnPos = new int [nCols] ;
     795  double *elementsPos = new double [nCols] ;
     796  rowStartPos = new CoinBigIndex [nRows] ;
     797  for (int i = 0 ; i < nRows ; i++) {
     798    double offset = 0.0 ;
     799    const CoinBigIndex rstarti = rowStart[i] ;
     800    const CoinBigIndex rendi = rstarti+rowLength[i] ;
     801    rowStart[i] = newSize ;
     802    const CoinBigIndex saveNewStarti = newSize ;
     803    int nPos = 0 ;
     804    for (CoinBigIndex jj = rstarti ; jj < rendi ; jj++) {
     805      int j = column[jj] ;
     806      double value = elements[jj] ;
     807      if (colUpper[j] > colLower[j]) {
     808        if (value < 0.0) {
     809          elements[newSize] = value ;
     810          column[newSize++] = j ;
     811        } else {
     812          elementsPos[nPos] = value ;
     813          columnPos[nPos++] = j ;
     814        }
     815      } else {
     816        offset += colUpper[j]*value ;
     817      }
     818    }
     819    rowStartPos[i] = newSize ;
     820    for (int k = 0 ; k < nPos ; k++) {
     821      elements[newSize] = elementsPos[k] ;
     822      column[newSize++] = columnPos[k] ;
     823    }
     824    rowLength[i] = newSize-saveNewStarti ;
     825    if (offset) {
     826      if (rowLower[i] > -weakRhs)
     827        rowLower[i] -= offset ;
     828      if (rowUpper[i] < weakRhs)
     829        rowUpper[i] -= offset ;
     830    }
     831  }
     832  delete [] columnPos ;
     833  delete [] elementsPos ;
     834  rowStart[nRows] = newSize ;
     835  rowCopy->setNumElements(newSize) ;
     836
     837  return (true) ;
     838}
     839
     840}  // end file local namespace
     841
     842
    51843//-------------------------------------------------------------------
    52844// Generate disaggregation cuts
     
    57849
    58850# if CGL_DEBUG > 0
    59   std::cout << "Entering CglProbing::generateCuts." << std::endl ;
     851  std::cout
     852    << "Entering CglProbing::generateCuts, matrix " << si.getNumRows()
     853    << " x " << si.getNumCols() << "." << std::endl ;
    60854/*
    61855  Note that the debugger must be activated before generateCuts is called.
     
    85879/*
    86880  Setup. Create arrays that will hold row and column bounds while we're
     881  working. Note that the row bound arrays are allocated with one extra
     882  position, in case we want to incorporate the objective as a constraint while
    87883  working.
    88884
     
    111907  CglTreeInfo info = info2 ;
    112908/*
    113   Do the work.
     909  Do the work. There's no guarantee about the state of the row and column
     910  bounds arrays if we come back infeasible.
    114911*/
    115912  int ninfeas =
    116913        gutsOfGenerateCuts(si,cs,rowLower,rowUpper,colLower,colUpper,&info) ;
     914# if CGL_DEBUG > 0
     915  { int errs = 0 ;
     916    const CoinPackedMatrix *mtx = si.getMatrixByRow() ;
     917    errs = mtx->verifyMtx(2) ;
     918    if (errs > 0) {
     919      std::cout
     920        << "    generateCuts: verifyMtx failed, "
     921        << errs << " errors." << std::endl ;
     922      assert (false) ;
     923    } else {
     924      std::cout << "    generateCuts: matrix verified." << std::endl ;
     925    }
     926  }
     927# endif
    117928/*
    118929  Infeasibility is indicated by a stylized infeasible column cut.
     
    168979{
    169980# if CGL_DEBUG > 0
    170   std::cout << "Entering CglProbing::generateCutsAndModify." << std::endl ;
     981  std::cout
     982    << "Entering CglProbing::generateCutsAndModify, matrix " << si.getNumRows()
     983    << " x " << si.getNumCols() << "." << std::endl ;
    171984  const OsiRowCutDebugger *debugger = si.getRowCutDebugger() ;
    172985  if (debugger) {
     
    2231036  Do the work.
    2241037*/
     1038# if CGL_DEBUG > 0
     1039  { int errs = 0 ;
     1040    const CoinPackedMatrix *mtx = si.getMatrixByRow() ;
     1041    errs = mtx->verifyMtx(2) ;
     1042    if (errs > 0) {
     1043      std::cout
     1044        << "    generateCutsAndModify: verifyMtx failed pre-cutgen, "
     1045        << errs << " errors." << std::endl ;
     1046      assert (errs == 0) ;
     1047    }
     1048  }
     1049# endif
     1050/*
     1051  Do the work. There's no guarantees about the state of the row and column
     1052  bounds arrays if we come back infeasible.
     1053*/
    2251054  int ninfeas = gutsOfGenerateCuts(si,cs,
    2261055                                   rowLower,rowUpper,colLower,colUpper,info) ;
     1056# if CGL_DEBUG > 0
     1057  { int errs = 0 ;
     1058    const CoinPackedMatrix *mtx = si.getMatrixByRow() ;
     1059    errs = mtx->verifyMtx(2) ;
     1060    if (errs > 0) {
     1061      std::cout
     1062        << "    generateCutsAndModify: verifyMtx failed post-cutgen, "
     1063        << errs << " errors." << std::endl ;
     1064      assert (errs == 0) ;
     1065    }
     1066  }
     1067# endif
    2271068/*
    2281069  Infeasibility is indicated by a stylized infeasible column cut.
     
    2901131}
    2911132
    292 bool analyze(const OsiSolverInterface * solverX, char * intVar,
    293              double * lower, double * upper)
    294 {
    295   OsiSolverInterface * solver = solverX->clone();
    296   const double *objective = solver->getObjCoefficients() ;
    297   int numberColumns = solver->getNumCols() ;
    298   int numberRows = solver->getNumRows();
    299   double direction = solver->getObjSense();
    300   int iRow,iColumn;
    301 
    302   // Row copy
    303   CoinPackedMatrix matrixByRow(*solver->getMatrixByRow());
    304   const double * elementByRow = matrixByRow.getElements();
    305   const int * column = matrixByRow.getIndices();
    306   const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
    307   const int * rowLength = matrixByRow.getVectorLengths();
    308 
    309   // Column copy
    310   CoinPackedMatrix  matrixByCol(*solver->getMatrixByCol());
    311   const double * element = matrixByCol.getElements();
    312   const int * row = matrixByCol.getIndices();
    313   const CoinBigIndex * columnStart = matrixByCol.getVectorStarts();
    314   const int * columnLength = matrixByCol.getVectorLengths();
    315 
    316   const double * rowLower = solver->getRowLower();
    317   const double * rowUpper = solver->getRowUpper();
    318 
    319   char * ignore = new char [numberRows];
    320   int * which = new int[numberRows];
    321   double * changeRhs = new double[numberRows];
    322   memset(changeRhs,0,numberRows*sizeof(double));
    323   memset(ignore,0,numberRows);
    324   int numberChanged=0;
    325   bool finished=false;
    326   while (!finished) {
    327     int saveNumberChanged = numberChanged;
    328     for (iRow=0;iRow<numberRows;iRow++) {
    329       int numberContinuous=0;
    330       double value1=0.0,value2=0.0;
    331       bool allIntegerCoeff=true;
    332       double sumFixed=0.0;
    333       int jColumn1=-1,jColumn2=-1;
    334       for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
    335         int jColumn = column[j];
    336         double value = elementByRow[j];
    337         if (upper[jColumn] > lower[jColumn]+1.0e-8) {
    338           if (!intVar[jColumn]) {
    339             if (numberContinuous==0) {
    340               jColumn1=jColumn;
    341               value1=value;
    342             } else {
    343               jColumn2=jColumn;
    344               value2=value;
    345             }
    346             numberContinuous++;
    347           } else {
    348             if (fabs(value-floor(value+0.5))>1.0e-12)
    349               allIntegerCoeff=false;
    350           }
    351         } else {
    352           sumFixed += lower[jColumn]*value;
    353         }
    354       }
    355       double low = rowLower[iRow];
    356       if (low>-1.0e20) {
    357         low -= sumFixed;
    358         if (fabs(low-floor(low+0.5))>1.0e-12)
    359           allIntegerCoeff=false;
    360       }
    361       double up = rowUpper[iRow];
    362       if (up<1.0e20) {
    363         up -= sumFixed;
    364         if (fabs(up-floor(up+0.5))>1.0e-12)
    365           allIntegerCoeff=false;
    366       }
    367       if (!allIntegerCoeff)
    368         continue; // can't do
    369       if (numberContinuous==1) {
    370         // see if really integer
    371         // This does not allow for complicated cases
    372         if (low==up) {
    373           if (fabs(value1)>1.0e-3) {
    374             value1 = 1.0/value1;
    375             if (fabs(value1-floor(value1+0.5))<1.0e-12) {
    376               // integer
    377               numberChanged++;
    378               intVar[jColumn1]=77;
    379             }
    380           }
    381         } else {
    382           if (fabs(value1)>1.0e-3) {
    383             value1 = 1.0/value1;
    384             if (fabs(value1-floor(value1+0.5))<1.0e-12) {
    385               // This constraint will not stop it being integer
    386               ignore[iRow]=1;
    387             }
    388           }
    389         }
    390       } else if (numberContinuous==2) {
    391         if (low==up) {
    392           /* need general theory - for now just look at 2 cases -
    393              1 - +- 1 one in column and just costs i.e. matching objective
    394              2 - +- 1 two in column but feeds into G/L row which will try and minimize
    395              (take out 2 for now - until fixed)
    396           */
    397           if (fabs(value1)==1.0&&value1*value2==-1.0&&!lower[jColumn1]
    398               &&!lower[jColumn2]&&columnLength[jColumn1]==1&&columnLength[jColumn2]==1) {
    399             int n=0;
    400             int i;
    401             double objChange=direction*(objective[jColumn1]+objective[jColumn2]);
    402             double bound = CoinMin(upper[jColumn1],upper[jColumn2]);
    403             bound = CoinMin(bound,1.0e20);
    404             for ( i=columnStart[jColumn1];i<columnStart[jColumn1]+columnLength[jColumn1];i++) {
    405               int jRow = row[i];
    406               double value = element[i];
    407               if (jRow!=iRow) {
    408                 which[n++]=jRow;
    409                 changeRhs[jRow]=value;
    410               }
    411             }
    412             for ( i=columnStart[jColumn2];i<columnStart[jColumn2]+columnLength[jColumn2];i++) {
    413               int jRow = row[i];
    414               double value = element[i];
    415               if (jRow!=iRow) {
    416                 if (!changeRhs[jRow]) {
    417                   which[n++]=jRow;
    418                   changeRhs[jRow]=value;
    419                 } else {
    420                   changeRhs[jRow]+=value;
    421                 }
    422               }
    423             }
    424             if (objChange>=0.0) {
    425               // see if all rows OK
    426               bool good=true;
    427               for (i=0;i<n;i++) {
    428                 int jRow = which[i];
    429                 double value = changeRhs[jRow];
    430                 if (value) {
    431                   value *= bound;
    432                   if (rowLength[jRow]==1) {
    433                     if (value>0.0) {
    434                       double rhs = rowLower[jRow];
    435                       if (rhs>0.0) {
    436                         double ratio =rhs/value;
    437                         if (fabs(ratio-floor(ratio+0.5))>1.0e-12)
    438                           good=false;
    439                       }
    440                     } else {
    441                       double rhs = rowUpper[jRow];
    442                       if (rhs<0.0) {
    443                         double ratio =rhs/value;
    444                         if (fabs(ratio-floor(ratio+0.5))>1.0e-12)
    445                           good=false;
    446                       }
    447                     }
    448                   } else if (rowLength[jRow]==2) {
    449                     if (value>0.0) {
    450                       if (rowLower[jRow]>-1.0e20)
    451                         good=false;
    452                     } else {
    453                       if (rowUpper[jRow]<1.0e20)
    454                         good=false;
    455                     }
    456                   } else {
    457                     good=false;
    458                   }
    459                 }
    460               }
    461               if (good) {
    462                 // both can be integer
    463                 numberChanged++;
    464                 intVar[jColumn1]=77;
    465                 numberChanged++;
    466                 intVar[jColumn2]=77;
    467               }
    468             }
    469             // clear
    470             for (i=0;i<n;i++) {
    471               changeRhs[which[i]]=0.0;
    472             }
    473           }
    474         }
    475       }
    476     }
    477     for (iColumn=0;iColumn<numberColumns;iColumn++) {
    478       if (upper[iColumn] > lower[iColumn]+1.0e-8&&!intVar[iColumn]) {
    479         double value;
    480         value = upper[iColumn];
    481         if (value<1.0e20&&fabs(value-floor(value+0.5))>1.0e-12)
    482           continue;
    483         value = lower[iColumn];
    484         if (value>-1.0e20&&fabs(value-floor(value+0.5))>1.0e-12)
    485           continue;
    486         bool integer=true;
    487         for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
    488           int iRow = row[j];
    489           if (!ignore[iRow]) {
    490             integer=false;
    491             break;
    492           }
    493         }
    494         if (integer) {
    495           // integer
    496           numberChanged++;
    497           intVar[iColumn]=77;
    498         }
    499       }
    500     }
    501     finished = numberChanged==saveNumberChanged;
    502   }
    503   bool feasible=true;
    504   for (iColumn=0;iColumn<numberColumns;iColumn++) {
    505     if (intVar[iColumn]==77) {
    506       if (upper[iColumn]>1.0e20) {
    507         upper[iColumn] = 1.0e20;
    508       } else {
    509         upper[iColumn] = floor(upper[iColumn]+1.0e-5);
    510       }
    511       if (lower[iColumn]<-1.0e20) {
    512         lower[iColumn] = -1.0e20;
    513       } else {
    514         lower[iColumn] = ceil(lower[iColumn]-1.0e-5);
    515         if (lower[iColumn]>upper[iColumn])
    516           feasible=false;
    517       }
    518       if (lower[iColumn]==0.0&&upper[iColumn]==1.0)
    519         intVar[iColumn]=1;
    520       else if (lower[iColumn]==upper[iColumn])
    521         intVar[iColumn]=0;
    522       else
    523         intVar[iColumn]=2;
    524     }
    525   }
    526   delete [] which;
    527   delete [] changeRhs;
    528   delete [] ignore;
    529   //if (numberChanged)
    530   //printf("%d variables could be made integer\n",numberChanged);
    531   delete solver;
    532   return feasible;
    533 }
    534 
    5351133
    5361134/*
     
    5381136        have an extra slot that can be used for bounds on the objective,
    5391137        in the event that gutsOfGenerateCuts decides to put it in the matrix.
     1138
     1139  rowLower, rowUpper, colLower, and colUpper should be allocated by the caller
     1140  and will be filled in by gutsOfGenerateCuts.
    5401141*/
    5411142int CglProbing::gutsOfGenerateCuts (const OsiSolverInterface &si,
     
    5481149# if CGL_DEBUG > 0
    5491150  std::cout
    550     << "Entering CglProbing::gutsOfGenerateCuts, mode "
    551     << mode_ << std::endl ;
     1151    << "Entering CglProbing::gutsOfGenerateCuts, mode " << mode_ << ", " ;
     1152  if (!info->inTree) std::cout << "not " ;
     1153  std::cout << "in tree, pass " << info->pass << "." << std::endl ;
    5521154# endif
    5531155
    554   int nRows ;
    555  
    556   CoinPackedMatrix *rowCopy = NULL ;
    5571156  int numberRowCutsBefore = cs.sizeRowCuts() ;
    5581157
    559 /*
    560   Establish a cutoff -- nontrivial, if possible. The client can forbid use of
    561   the objective by setting usingObjective_ to -1.
    562 */
    563   double cutoff ;
    564   bool cutoff_available = si.getDblParam(OsiDualObjectiveLimit,cutoff) ;
    565   if (!cutoff_available || usingObjective_ < 0) {
    566     cutoff = si.getInfinity() ;
     1158  int ninfeas = 0 ;
     1159  bool feasible = true ;
     1160/*
     1161  If we don't have a snapshot, force mode 0 to mode 1.
     1162
     1163  TODO: Ah ... maybe a `You're confused' message and an error return would be
     1164        more appropriate? Seems like a serious algorithm error at a higher
     1165        level.  -- lh, 100917 --
     1166
     1167        Or maybe the client is within the documentation. After all, John's
     1168        comment for mode 0 says `only available with snapshot; otherwise as
     1169        mode 1'.  -- lh, 101021 --
     1170
     1171        So, maybe I should try to figure out if requesting mode 0 results in
     1172        creating a snapshot somewhere along the way.  -- lh, 110204 --
     1173*/
     1174  int mode = mode_ ;
     1175  if (!rowCopy_ && mode == 0) {
     1176    mode = 1 ;
     1177#   if CGL_DEBUG > 0
     1178    std::cout << "  forcing mode 0 to mode 1; no snapshot." << std::endl ;
     1179#   endif
    5671180  }
    568   cutoff *= si.getObjSense() ;
    569 /*
    570   This looks like a sure failure if we're maximising and grabbed the solver's
    571   infinity.  In any event, why is large and positive ok but large and
    572   negative is not? --lh, 100917 --
    573 */
    574   if (fabs(cutoff) > 1.0e30)
    575     assert (cutoff > 1.0e30) ;
    576 
    577   int mode = mode_ ;
    578   int nCols = si.getNumCols();
     1181
     1182  int maxProbe = info->inTree ? maxProbe_ : maxProbeRoot_ ;
     1183  int maxRowLen = info->inTree ? maxElements_ : maxElementsRoot_ ;
     1184  // Forcing probing to consider really dense constraints wasn't good,
     1185  // even at the root.
     1186  // if (!info->inTree && !info->pass) maxRowLen = nCols ;
     1187
    5791188/*
    5801189  Get a modifiable array of variable types. Force reevaluation of the
     
    5821191  bounds.
    5831192*/
     1193  int nCols = si.getNumCols();
    5841194  const char *intVarOriginal = si.getColType(true) ;
    5851195  char *intVar = CoinCopyOfArray(intVarOriginal,nCols) ;
     
    5891199  CoinMemcpyN(si.getColLower(),nCols,colLower) ;
    5901200  CoinMemcpyN(si.getColUpper(),nCols,colUpper) ;
    591   const double * colsol = si.getColSolution() ;
     1201  const double *colsol = si.getColSolution() ;
    5921202/*
    5931203  Stage 1: Preliminary processing.
     
    6171227    }
    6181228  }
    619   bool feasible = true ;
    6201229/*
    6211230  If we're at the root of the search tree, scan the constraint system to see
     
    6231232  Convert them to integer type.
    6241233
    625   TODO: If we lose feasibility here, we should clean up and bail. This is not
    626         handled well in subsequent blocks of code.  -- lh, 101021 --
     1234  If we lose feasibility here, cut our losses and bail out while it's still
     1235  uncomplicated.
     1236
     1237  TODO: Really, analyze should return a count instead of a boolean.
     1238        -- lh, 110204 --
    6271239*/
    6281240  if (!info->inTree && !info->pass) {
    6291241    feasible = analyze(&si,intVar,colLower,colUpper) ;
     1242#   if CGL_DEBUG > 1
     1243    std::cout
     1244     << "  Analyze completed; " << ((feasible)?"feasible":"infeasible")
     1245     << "." << std::endl ;
     1246#   endif
     1247    if (!feasible) {
     1248      delete[] intVar ;
     1249      return (1) ;
     1250    }
    6301251  }
    631 # if CGL_DEBUG > 1
    632   std::cout
    633    << "  Analyze completed; " << ((feasible)?"feasible":"infeasible")
    634    << "." << std::endl ;
    635 # endif
    636 /*
    637   Attempt to tighten bounds based on reduced costs. Consider movement by +/-1
    638   and use that as a filter. For integer variables, this will fix the variable
    639   at the current bound. For continuous variables, we can restrict the gap
    640   between lower and upper bound to less than 1.0.
    641 
    642   jjf: Should be in CbcCutGenerator and check if basic.
    643 
    644   lh: (100917) Not sure why that's relevant, except that reduced cost is
    645       zero by definition for basic variables. Why not just test for d<j> = 0?
    646  
    647   lh: (100917) Seems like a prime candidate for a subroutine (inlined, perhaps).
    648       And that's a really good explanation for the redundant code. This
    649       probably *was* a subroutine at one time.
    650 
    651   This block disabled as of 100917 (PROBING_EXTRA_STUFF == false). Probably
    652   disabled further back than that.
    653 
    654   lh: (101021) And I'm asking myself `Why why?' in the print statements below.
    655       Apparently John was of the opinion that anything that could be fixed here
    656       should have been fixed elsewhere. Did this functionality actually move
    657       elsewhere, so that those print statements are just telltales in case
    658       something slipped by?
    659 */
    660   if (feasible && PROBING_EXTRA_STUFF) {
    661     const double * djs = si.getReducedCost() ;
    662 
    663     // Begin redundant code block -- we did this just above!
    664     const double * colsol = si.getColSolution() ;
    665     double direction = si.getObjSense() ;
    666     double cutoff ;
    667     bool cutoff_available = si.getDblParam(OsiDualObjectiveLimit,cutoff) ;
    668     if (!cutoff_available || usingObjective_ < 0) {
    669       cutoff = si.getInfinity() ;
    670     }
    671     cutoff *= direction ;
    672     if (fabs(cutoff) > 1.0e30)
    673       assert (cutoff > 1.0e30) ;
    674     // End redundant code block
    675 
    676     double current = si.getObjValue() ;
    677     current *= direction ;
    678     double gap = CoinMax(cutoff-current,1.0e-1) ;
    679     for (int i = 0; i < nCols; ++i) {
    680       double djValue = djs[i]*direction ;
    681       // Consider only variables that are not fixed.
    682       if (colUpper[i]-colLower[i] > 1.0e-8) {
    683         // If we're currently at lower bound ...
    684         if (colsol[i] < colLower[i]+primalTolerance_) {
    685           // ... and +1 will drive us over the cutoff ...
    686           if (djValue > gap) {
    687             // ... and we have to move at least +1 ...
    688             if (si.isInteger(i)) {
    689               // ... fix at lower bound.
    690               printf("why int %d not fixed at lb\n",i) ;
    691               colUpper[i] = colLower[i] ;
    692             } else {
    693               // ... if a fractional move is possible, calculate the max
    694               // movement and crank down the upper bound.
    695               double newUpper = colLower[i] + gap/djValue ;
    696               if (newUpper < colUpper[i]) {
    697                 //printf("%d ub from %g to %g\n",i,colUpper[i],newUpper) ;
    698                 colUpper[i] = CoinMax(newUpper,colLower[i]+1.0e-5) ;
    699               }
    700             }
    701           }
    702         } else if (colsol[i] > colUpper[i]-primalTolerance_) {
    703           if (-djValue > gap) {
    704             if (si.isInteger(i)) {
    705               printf("why int %d not fixed at ub\n",i) ;
    706               colLower[i] = colUpper[i] ;
    707             } else {
    708               double newLower = colUpper[i] + gap/djValue ;
    709               if (newLower>colLower[i]) {
    710                 //printf("%d lb from %g to %g\n",i,colLower[i],newLower) ;
    711                 colLower[i]= CoinMin(newLower,colUpper[i]-1.0e-5) ;
    712               }
    713             }
    714           }
    715         }
    716       }
    717     }
    718   }  // End cost bound tightening using reduced costs.
    719 
    720   int ninfeas = 0 ;
    721   // Set up maxes
    722   int maxProbe = info->inTree ? maxProbe_ : maxProbeRoot_ ;
    723   int maxElements = info->inTree ? maxElements_ : maxElementsRoot_ ;
    724 
    725   // Forcing probing to consider really dense constraints wasn't good,
    726   // even at the root.
    727   // if (!info->inTree && !info->pass) maxElements = nCols ;
    728 
    729   // Get objective offset
     1252/*
     1253  Get the objective offset and try to establish a nontrivial cutoff. The
     1254  client can forbid use of the objective by setting usingObjective_ to -1.
     1255  If the si claims a nontrivial cutoff, check that it's reasonable.
     1256
     1257  TODO: Check through the code and make sure we're using cutoff and offset
     1258        correctly. Should add the offset in one place and keep the value
     1259        for later use.  -- lh, 110205 --
     1260
     1261  TODO: Down in probeClique, there's a clause that forces the cutoff to
     1262        infty for mode 0.
     1263*/
    7301264  double offset ;
    7311265  si.getDblParam(OsiObjOffset,offset) ;
    732 # if CGL_DEBUG > 3
    733   // print it just once
    734   if (offset && !info->inTree && !info->pass)
    735     printf("CglProbing obj offset %g\n",offset) ;
     1266  double cutoff ;
     1267  bool cutoff_available = si.getDblParam(OsiDualObjectiveLimit,cutoff) ;
     1268  if (!cutoff_available) {
     1269    cutoff = si.getInfinity() ;
     1270  } else {
     1271    cutoff *= si.getObjSense() ;
     1272    if (cutoff > 1.0e30) cutoff_available = false ;
     1273  }
     1274  bool useObj = (cutoff_available && (usingObjective_ >= 1)) ;
     1275  bool useCutoff = (cutoff_available && (usingObjective_ >= 0)) ;
     1276# if CGL_DEBUG > 1
     1277  std::cout << "  cutoff " ;
     1278  if (!cutoff_available)
     1279    std::cout << "not available" ;
     1280  else
     1281    std::cout << cutoff ;
     1282  std::cout << ", offset " << offset << "." << std::endl ;
     1283  if (useCutoff)
     1284    std::cout << "  using cutoff" ;
     1285  else
     1286    std::cout << "  not using cutoff" ;
     1287  if (useObj)
     1288    std::cout << ", using objective as constraint" ;
     1289  else
     1290    std::cout << ", not using objective as constraint" ;
     1291  std::cout << "." << std::endl ;
    7361292# endif
    737 /*
    738   If we have a snapshot handy, most of the setup is done. Otherwise, we need to
    739   create row- and column-major copies of the constraint matrix. If the client
    740   has specified mode 0 (work from snapshot), it is confused, but we silently
    741   correct it.
    742 
    743   TODO: Ah ... maybe a `You're confused' message and an error return would be
    744         more appropriate? Seems like a serious algorithm error at a higher
    745         level.  -- lh, 100917 --
    746 
    747         Or maybe the client is within the documentation. After all, John's
    748         comment for mode 0 says `only available with snapshot; otherwise as
    749         mode 1'.  -- lh, 101021 --
    750 
    751   TODO: Note that the column-major copy is not created until much further
    752         down, after we've settled on the set of rows to process and applied
    753         various techniques for row strengthening.   -- lh, 101021 --
    754 */
    755   if (!rowCopy_) {
    756     nRows = si.getNumRows();
    757     if (mode == 0)
    758       mode = 1 ;
    759 /*
    760   Make the row copy and create modifiable row bound arrays. If we have a
    761   nontrivial cutoff and are not forbidden from using it, add in the objective
    762   coefficients.
    763 
    764   TODO: For the objective, we make an entry in rowLower and rowUpper and
    765         install our bounds. That implies that the arrays passed in must have
    766         sufficient capacity.  -- lh, 100917 --
    767 
    768   TODO: Really, there should be a boolean for using the objective, so we don't
    769         keep reevaluating this condition. And why do we generate a boolean,
    770         maximise, instead of using sense*objcoeff, as is common elsewhere?
    771         -- lh, 100917 --
    772 
    773   TODO: No column copy?  -- lh, 100917 --
    774         It's made far down below, at the end of the row strengthening block,
    775         after we've decided on the rows to keep.
    776 */
    777     if (cutoff < 1.0e30 && usingObjective_ > 0) {
    778       rowCopy = new CoinPackedMatrix(*si.getMatrixByRow(),1,nCols,false) ;
    779     } else {
    780       rowCopy = new CoinPackedMatrix(*si.getMatrixByRow()) ;
    781     }
    782     if (cutoff < 1.0e30 && usingObjective_ > 0) {
    783       int *columns = new int[nCols] ;
    784       double *elements = new double[nCols] ;
    785       int n = 0 ;
    786       const double *objective = si.getObjCoefficients() ;
    787       bool maximize = (si.getObjSense() == -1) ;
    788       for (int i = 0 ; i < nCols ; i++) {
    789         if (objective[i]) {
    790           elements[n] = (maximize) ? -objective[i] : objective[i] ;
    791           columns[n++] = i ;
    792         }
    793       }
    794       rowCopy->appendRow(n,columns,elements) ;
    795       delete [] columns ;
    796       delete [] elements ;
    797       CoinMemcpyN(si.getRowLower(),nRows,rowLower) ;
    798       CoinMemcpyN(si.getRowUpper(),nRows,rowUpper) ;
    799       rowLower[nRows] = -DBL_MAX ;
    800       rowUpper[nRows] = cutoff+offset ;
    801       nRows++ ;
    802     } else {
    803       CoinMemcpyN(si.getRowLower(),nRows,rowLower) ;
    804       CoinMemcpyN(si.getRowUpper(),nRows,rowUpper) ;
    805     }
     1293
     1294/* EXTRA_PROBING_STUFF block (tighten using reduced cost) was here. */
     1295
     1296/*
     1297  Create the working row copy from the si's model or from a snapshot (only if
     1298  the user specified mode 0 and a snapshot exists).
     1299*/
     1300  CoinPackedMatrix *rowCopy = setupRowCopy(mode,useObj,cutoff,
     1301                                           offset,rowCopy_,
     1302                                           numberRows_,numberColumns_,
     1303                                           rowLower_,rowUpper_,
     1304                                           si,rowLower,rowUpper) ;
     1305  int nRows = rowCopy->getNumRows() ;
     1306/*
     1307  Groom the rowCopy: delete useless rows, substitute for fixed variables, sort
     1308  coefficients so that negative precede positive, and try some specialised
     1309  coefficient strengthening on the rows that remain. If groomModel comes
     1310  back false, it's an indication that probing isn't worth pursuing. Clean
     1311  up and go home while it's relatively easy.
     1312
     1313  realRows is the cross-reference between rows in rowCopy and rows of the
     1314  original model. rowStartPos is the first positive coefficient in a row.
     1315*/
     1316  int *realRows = 0 ;
     1317  int *rowStartPos = 0 ;
     1318  feasible = groomModel(useObj,maxRowLen,si,intVar,rowCopy,rowLower,rowUpper,
     1319                        colLower,colUpper,realRows,rowStartPos,info) ;
     1320  nRows = rowCopy->getNumRows() ;
     1321  if (!feasible) {
     1322    delete[] intVar ;
     1323    delete[] realRows ;
     1324    delete[] rowStartPos ;
     1325    delete rowCopy ;
     1326    return (0) ;
    8061327  }
    8071328/*
    808   If we're starting from a snapshot, make a copy for further processing. Make
    809   sure we adjust rowLower and rowUpper to accommodate the objective (the
    810   coefficients are already in).
    811 
    812   TODO: Another good reason to hive off snapshot as a separate class.
    813         -- lh, 100917 --
    814 
    815   TODO: Ummmm ... what if the objective isn't already in? I see, I'm checking
    816         usingObjective_, which (one hopes) is correctly set from before. But
    817         it's publically accessible. This should be an internal state thing.
    818         -- lh, 100917 --
    819  
    820   TODO: Note that we don't make a working copy of the column copy held by the
    821         snapshot.  -- lh, 101021 --
    822 
    823   TODO: Why do we need to avoid modifying the snapshot? Could we get by without
    824         making yet another copy?  -- lh, 101125 --
    825 */
    826   else {
    827     nRows = numberRows_ ;
    828     // Sanity -- our cached copy must have the same number of columns!
    829     assert(nCols == numberColumns_) ;
    830     rowCopy = new CoinPackedMatrix(*rowCopy_) ;
    831     // More sanity.
    832     assert (rowCopy_->getNumRows() == numberRows_) ;
    833     rowLower = new double[nRows] ;
    834     rowUpper = new double[nRows] ;
    835     CoinMemcpyN(rowLower_,nRows,rowLower) ;
    836     CoinMemcpyN(rowUpper_,nRows,rowUpper) ;
    837     if (usingObjective_ > 0) {
    838       rowLower[nRows-1] = -DBL_MAX ;
    839       rowUpper[nRows-1] = cutoff+offset ;
    840     }
    841   }
    842 /*
    843   Remove rows with too many elements and rows that are deemed useless for
    844   further probing.  Try to strengthen rows in place, if allowed. Deal with
    845   fixed variables. We're going to do this in place by editing the internal
    846   data structures of rowCopy.
    847 
    848   TODO: At first glance, it looks like nRealRows is just so we can avoid
    849         tossing the objective (that was John's comment). But then why not
    850         just set the upper bound to nRows-1 or nRows-2, depending on
    851         usingObjective_?  -- lh, 100917 --
    852 
    853   TODO: This block of code looks like an excellent candidate for a subroutine.
    854         Maybe it once was, given that it's wrapped in its own code block.
    855         -- lh, 100923 --
    856 
    857   TODO: OUTRUBBISH is disabled, so either it didn't improve performance or
    858         it's buggy. -- lh, 100923
    859 */
    860   // int i ;
    861   CoinBigIndex * rowStartPos = NULL ;
    862   int * realRows = NULL ;
    863   // Begin row strengthen block
    864   {
    865     //#define OUTRUBBISH
    866 
    867     int *rowLength = rowCopy->getMutableVectorLengths();
    868     double * elements = rowCopy->getMutableElements() ;
    869     int * column = rowCopy->getMutableIndices() ;
    870     CoinBigIndex *rowStart = rowCopy->getMutableVectorStarts() ;
    871 #ifdef OUTRUBBISH
    872     double large = 1.0e3 ;
    873 #endif
    874     int nDelete = 0 ;
    875     int nKeep = 0 ;
    876     int *which = new int[nRows] ;
    877     int nElements = rowCopy->getNumElements() ;
    878     int nTotalOut = 0 ;
    879     int nRealRows = si.getNumRows() ;
    880 /*
    881   Preliminary assessment for dense rows and free rows.  Count the number of
    882   coefficients that we'll remove. We're looking for rows that are too long or
    883   have row bounds so weak as to be ineffective. Don't toss the objective
    884   (index larger than nRealRows.
    885 */
    886     for (int i = 0 ; i < nRows ; i++) {
    887       if (rowLength[i] > maxElements ||
    888           (rowLower[i] < -1.0e20 && rowUpper[i] > 1.0e20)) {
    889         if (i < nRealRows)
    890           nTotalOut += rowLength[i] ;
    891       }
    892     }
    893 /*
    894   If the total coefficients in dense rows amounts to less than 1/10 of the
    895   matrix, just deal with it (maxElements set to nCols won't exclude anything).
    896 */
    897     if (nTotalOut*10 < nElements)
    898       maxElements = nCols ;
    899 #ifdef OUTRUBBISH
    900     int nExtraDel = 0 ;
    901 #endif
    902 /*
    903   Main processing loop(1).
    904  
    905   Tossing the junk will significantly reduce the constraint system, so let's
    906   get to it.
    907 
    908   Examine each row for length and effectiveness.  We can summarily toss free
    909   rows and rows that are too dense.
    910 
    911   TODO: We carefully exclude the objective from deletion for density, but note
    912         that 1e30 is the break for adding it and 1e20 is the break for removing
    913         it. Admittedly, between 1e20 and 1e30 there's not likely to be much
    914         effect, but it'd be nice to be consistent.  -- lh, 100923 --
    915 */
    916     for (int i = 0 ; i < nRows ; i++) {
    917       if ((rowLength[i] > maxElements && i < nRealRows) ||
    918           (rowLower[i] < -1.0e20 && rowUpper[i] > 1.0e20)) {
    919         which[nDelete++] = i ;
    920       } else {
    921 /*
    922   This row passes the coefficient limit and has at least one finite bound.
    923   But is it useful for bound propagation and probing?
    924 
    925   TODO: This is lhs bound evaluation. Should be done up front.
    926         -- lh, 100917 --
    927 
    928   TODO: The OUTRUBBISH and !OUTRUBBISH blocks are in no way mutually exclusive,
    929         even though that's the configuration here. It's easy to conceive of a
    930         blend in which the OUTRUBBISH block runs to remove `weak' constraints
    931         at some threshold. Constraints that pass are then considered for
    932         strengthening.  -- lh, 100923 --
    933 */
    934 #ifdef OUTRUBBISH
    935 /*
    936   jjf: out all rows with infinite plus and minus
    937 
    938   It's actually a bit more extreme than that; large is set to 1.0e3 above,
    939   which is a fair bit short of infinity. Very aggressive removal of constraints
    940   perceived to be weak. -- lh, 100923 --
    941 */
    942         int nPlus = (rowUpper[i] > -large) ? 0 : 1 ;
    943         int nMinus = (rowLower[i] < large) ? 0 : 1 ;
    944         CoinBigIndex start = rowStart[i] ;
    945         CoinBigIndex end = start + rowLength[i] ;
    946         for (CoinBigIndex j = start ; j < end ; j++) {
    947           int iColumn = column[j] ;
    948           if (colUpper[iColumn] > large) {
    949             if (elements[j] > 0.0)
    950               nPlus++ ;
    951             else
    952               nMinus++ ;
    953           }
    954           if (colLower[iColumn] < -large) {
    955             if (elements[j] < 0.0)
    956               nPlus++ ;
    957             else
    958               nMinus++ ;
    959           }
    960         }
    961         if (!nPlus||!nMinus) {
    962           rowLower[nKeep] = rowLower[i] ;
    963           rowUpper[nKeep] = rowUpper[i] ;
    964           nKeep++ ;
    965         } else {
    966           nExtraDel++ ;
    967           which[nDelete++] = i ;
    968         }
    969 #else  // OUTRUBBISH
    970 /*
    971   The (way) less aggressive approach. We won't actually delete anything
    972   here.  We're only looking at rows that can be strengthened in place.
    973   Call this Strengthen(1). See notes for details.
    974 
    975   Process this row only if info says we can strengthen rows in place; only on
    976   the first pass for this node; only if the row has at least one infinite
    977   bound (lower or upper)
    978 */
    979         if (info->strengthenRow && !info->pass &&
    980             (rowLower[i] < -1.0e20 || rowUpper[i] > 1.0e20)) {
    981 /*
    982   Assess the row. If we find anything other than binary variables, we're
    983   immediately done.
    984 
    985   TODO: Really, the loop should abort as soon as nPlus > 1 & nMinus > 1.
    986         Presumably it doesn't because the test would take more time than
    987         just finishing the processing?
    988         -- lh, 100923 --
    989 
    990   TODO: coefficients of zero are lumped into nMinus, but arguably a coeff
    991         of zero shouldn't exist in a sparse matrix structure.
    992         -- lh, 100923 --
    993 
    994   TODO: Is the code not consistent about labelling binary variables? This is
    995         precisely the condition used to set the code for binary in intVar.
    996         -- lh, 100923 --
    997 
    998   TODO: Doesn't allow for nPlus == 1 and nMinus == 1. Maybe that's a special
    999         case that gets caught elsewhere? And why increment effectiveness
    1000         rather than simply assign to it?   -- lh, 100923 --
    1001 */
    1002           int nPlus = 0 ;
    1003           int nMinus = 0 ;
    1004           for (CoinBigIndex j = rowStart[i] ; j < rowLength[i] ; j++) {
    1005             int jColumn = column[j] ;
    1006             if (intVar[jColumn] &&
    1007                 colLower[jColumn] == 0.0 && colUpper[jColumn] == 1.0) {
    1008               double value = elements[j] ;
    1009               if (value > 0.0) {
    1010                 nPlus++ ;
    1011               } else {
    1012                 nMinus++ ;
    1013               }
    1014             } else {
    1015               nPlus = 2 ;
    1016               nMinus = 2 ;
    1017               break ;
    1018             }
    1019           }
    1020 /*
    1021   Look for the single positive coefficient in the row and reduce it. (Similarly
    1022   for negative coefficient in the next block.
    1023 
    1024   TODO: So why don't we just break from the loop after we've tweaked the
    1025         coefficient?  There's no real loop here to accumulate effectiveness.
    1026         -- lh, 100923 --
    1027 */
    1028           double effectiveness = 0.0 ;
    1029           if (nPlus == 1 && rowUpper[i] > 0.0 && rowUpper[i] < 1.0e10) {
    1030             // can make element smaller
    1031             for (CoinBigIndex j = rowStart[i] ; j < rowStart[i+1] ; j++) {
    1032               double value = elements[j] ;
    1033               if (value > 0.0) {
    1034                 elements[j] -= rowUpper[i] ;
    1035                 //printf("pass %d row %d plus el from %g to %g\n",info->pass,
    1036                 //     i,elements[j]+rowUpper[i],elements[j]) ;
    1037               }
    1038               effectiveness += fabs(elements[j]) ;
    1039             }
    1040             rowUpper[i] = 0.0 ;
    1041           } else
    1042           if (nMinus == 1 && rowLower[i] < 0.0 && rowLower[i] > -1.0e10) {
    1043             // can make element smaller in magnitude
    1044             for (CoinBigIndex j = rowStart[i] ; j < rowStart[i+1] ; j++) {
    1045               double value = elements[j] ;
    1046               if (value < 0.0) {
    1047                 elements[j] -= rowLower[i] ;
    1048                 //printf("pass %d row %d minus el from %g to %g\n",info->pass,
    1049                 //     i,elements[j]+rowLower[i],elements[j]) ;
    1050               }
    1051               effectiveness += fabs(elements[j]) ;
    1052             }
    1053             rowLower[i]=0.0 ;
    1054           }
    1055 /*
    1056   If the coefficient is nonzero, remember the result. We need to remember the
    1057   result because (as with everything done here) we're working with a copy.
    1058   This business of strengthening the coefficient must be propagated back out
    1059   to the client, who gets to make the decision about whether to use it.
    1060 
    1061   TODO: Why do we create a block-local object, then create a clone for later
    1062         use? Why not just use new?  -- lh, 100923 --
    1063 */
    1064           if (effectiveness) {
    1065             OsiRowCut rc ;
    1066             rc.setLb(rowLower[i]) ;
    1067             rc.setUb(rowUpper[i]) ;
    1068             int start = rowStart[i] ;
    1069             rc.setRow(rowLength[i],column+start,elements+start,false) ;
    1070             rc.setEffectiveness(effectiveness) ;
    1071             assert (!info->strengthenRow[i]) ;
    1072             info->strengthenRow[i] = rc.clone() ;
    1073           }
    1074         }
    1075         rowLower[nKeep] = rowLower[i] ;
    1076         rowUpper[nKeep] = rowUpper[i] ;
    1077         nKeep++ ;
    1078 #endif   // OUTRUBBISH
    1079       }
    1080     }  // End main processing loop(1)
    1081 /*
    1082   Did we delete anything in processing loop(1)? If so, deal with it. We need
    1083   to delete the rows, of course. There may also be bookkeeping.  In
    1084   particular, strengthened constraints are held in strengthenRow by the
    1085   original row index. We need a cross-reference, realRows, to keep track.
    1086 
    1087   TODO: Huh? What's the deal with the objective here? Why does it need
    1088         special treatment in realRows?  -- lh, 100923 --
    1089 
    1090   TODO: If we have no constraints remaining (nKeep == 0), then surely we could
    1091         skip this bit of work.  -- lh, 101021 --
    1092 
    1093   TODO: Because realRows is not allocated unless we delete rows from the
    1094         system, we will not be able to strengthen in place. See the matching
    1095         TODO in probe().  -- lh, 101209 --
    1096 */
    1097     if (nDelete) {
    1098 #ifdef OUTRUBBISH
    1099       if (nExtraDel) {
    1100         printf("%d rows deleted (extra %d)\n",nDelete,nExtraDel) ;
    1101       }
    1102 #else
    1103 #endif
    1104       if (info->strengthenRow) {
    1105         realRows = new int [nRows] ;
    1106         CoinZeroN(realRows,nRows) ;
    1107         for (int i = 0 ; i < nDelete ; i++)
    1108           realRows[which[i]] = -1 ;
    1109         int k = 0 ;
    1110         for (int i = 0 ; i < nRows ; i++) {
    1111           if (!realRows[i]) {
    1112             if (i < nRealRows)
    1113               realRows[k++] = i ; // keep
    1114             else
    1115               realRows[k++] = -1 ; // objective - discard
    1116           }
    1117         }
    1118       }
    1119       rowCopy->deleteRows(nDelete,which) ;
    1120       nRows = nKeep ;
    1121     }
    1122     delete [] which ;
    1123 /*
    1124   Do we have any constraints left to work with? If not, clean up and go home.
    1125 */
    1126     if (!nRows) {
    1127 #ifdef COIN_DEVELOP
    1128       printf("All rows too long for probing\n") ;
    1129 #endif
    1130       delete rowCopy ;
    1131       // Read `if we're working from a snapshot'
    1132       if (rowCopy_) {
    1133         delete [] rowLower ;
    1134         delete [] rowUpper ;
    1135       }
    1136       delete [] intVar ;
    1137 /*
    1138   jjf: and put back unreasonable bounds on integer variables
    1139 
    1140   TODO: Do we need to be so careful here? Why not just do a block copy?
    1141         -- lh, 100923 --
    1142 */
    1143       const double * trueLower = si.getColLower() ;
    1144       const double * trueUpper = si.getColUpper() ;
    1145       for (int i = 0 ; i < nCols ; i++) {
    1146         if (intVarOriginal[i] == 2) {
    1147           if (colUpper[i] == CGL_REASONABLE_INTEGER_BOUND)
    1148             colUpper[i] = trueUpper[i] ;
    1149           if (colLower[i] == -CGL_REASONABLE_INTEGER_BOUND)
    1150             colLower[i] = trueLower[i] ;
    1151         }
    1152       }
    1153       delete [] realRows ;
    1154       return 0 ;
    1155     }
    1156 /*
    1157   We're still inside the row strengthen block ...
    1158 
    1159   jjf: Out elements for fixed columns and sort
    1160 
    1161   In a few more words: Process the rows, adjusting rowLower and rowUpper for
    1162   fixed variables; the columns are physically removed. As we rewrite the
    1163   coefficients, sort them so that negative coefficients precede positive
    1164   coefficients.
    1165 
    1166   We need to refresh the mutable pointers; deleteRows above may have rendered
    1167   the previous pointers invalid.
    1168 
    1169   TODO: This code duplicates a block in snapshot(), except that the code in
    1170         snapshot doesn't substitute for the values of fixed variables. The
    1171         two blocks should probably be pulled out into a method.
    1172 */
    1173     elements = rowCopy->getMutableElements() ;
    1174     column = rowCopy->getMutableIndices() ;
    1175     rowStart = rowCopy->getMutableVectorStarts() ;
    1176     rowLength = rowCopy->getMutableVectorLengths();
    1177 
    1178 #if 0
    1179     // For debugging?
    1180     int nFixed = 0 ;
    1181     for (i = 0 ; i < nCols ; i++) {
    1182       if (colUpper[i] == colLower[i])
    1183         nFixed++ ;
    1184     }
    1185     printf("%d columns fixed\n",nFixed) ;
    1186 #endif
    1187 
    1188     CoinBigIndex newSize = 0 ;
    1189     int * column2 = new int[nCols] ;
    1190     double * elements2 = new double[nCols] ;
    1191     rowStartPos = new CoinBigIndex [nRows] ;
    1192     for (int i = 0 ; i < nRows ; i++) {
    1193       double offset = 0.0 ;
    1194       CoinBigIndex start = rowStart[i] ;
    1195       CoinBigIndex end = start+rowLength[i] ;
    1196       rowStart[i] = newSize ;
    1197       CoinBigIndex save = newSize ;
    1198       int nOther = 0 ;
    1199       for (CoinBigIndex j = start ; j < end ; j++) {
    1200         int iColumn = column[j] ;
    1201         if (colUpper[iColumn] > colLower[iColumn]) {
    1202           double value = elements[j] ;
    1203           if (value < 0.0) {
    1204             elements[newSize] = value ;
    1205             column[newSize++] = iColumn ;
    1206           } else {
    1207             elements2[nOther] = value ;
    1208             column2[nOther++] = iColumn ;
    1209           }
    1210         } else {
    1211           offset += colUpper[iColumn]*elements[j] ;
    1212         }
    1213       }
    1214       rowStartPos[i] = newSize ;
    1215       for (int k = 0 ; k < nOther ; k++) {
    1216         elements[newSize] = elements2[k] ;
    1217         column[newSize++] = column2[k] ;
    1218       }
    1219       rowLength[i] = newSize-save ;
    1220       if (offset) {
    1221         if (rowLower[i] > -1.0e20)
    1222           rowLower[i] -= offset ;
    1223         if (rowUpper[i] < 1.0e20)
    1224           rowUpper[i] -= offset ;
    1225       }
    1226     }
    1227     delete [] column2 ;
    1228     delete [] elements2 ;
    1229     rowStart[nRows] = newSize ;
    1230     rowCopy->setNumElements(newSize) ;
    1231   }
    1232   // End row strengthen block
    1233 
    1234 /*
    1235   Here, finally, we're making a column-major copy, after we've winnowed the
    1236   rows for the set we want to work with.
    1237 
    1238   TODO: Still, I'm asking myself, ``Why does the snapshot keep a colum-major
     1329  Make a column-major copy, after we've winnowed the rows to the set we want
     1330  to work with.
     1331
     1332  TODO: I'm asking myself, ``Why does the snapshot keep a colum-major
    12391333        copy?  -- lh, 101021 --
    12401334*/
     
    14671561  processing are tagged with -1 in markR
    14681562*/
    1469 #if 0
    1470 /*
    1471   TODO: This bit of code seems to have migrated up to the head of Stage 1,
    1472         main processing loop (1), where it's used to limit the constraints
    1473         installed in rowCopy. Presumably this is the origin of the abort().
    1474         We should never see a dense row in the working matrix.
    1475 */
    1476         // Only look at short rows
    1477         for (i=0;i<nRows;i++) {
    1478           if (rowLength[i]>maxElements)
    1479             abort(); //markR[i]=-2 ;
    1480         }
    1481 #endif
    14821563/*
    14831564  TODO: Eh? Sort by index for repeatability? Seems pointless, and it is
     
    14921573          ninfeas = probe(si,cs,colLower,colUpper,rowCopy,columnCopy,
    14931574                          rowStartPos,realRows,rowLower,rowUpper,
    1494                           intVar,minR,maxR,markR,info) ;
     1575                          intVar,minR,maxR,markR,info,
     1576                          useObj,useCutoff,cutoff) ;
    14951577        } else {
    14961578          ninfeas = probeCliques(si,debugger,cs,
    14971579                                 colLower,colUpper,rowCopy,columnCopy,
    14981580                                 realRows,rowLower,rowUpper,
    1499                                  intVar,minR,maxR,markR,info) ;
     1581                                 intVar,minR,maxR,markR,info,
     1582                                 useObj,useCutoff,cutoff) ;
    15001583        }
    15011584      } // end maxProbe > 0
     
    15031586  } // end mode != 0
    15041587/*
    1505   End of block for mode = 1 or 2; begin mode = 0. This works with the snapshot
    1506   and does not do bound propagation.
     1588  End of block for mode = 1 or 2; begin mode = 0.
     1589 
     1590  Mode 0 works with the snapshot and does not do bound propagation.
    15071591
    15081592  TODO: In fact, it's unclear to me just what we're doing here. The activity
     
    15831667    OsiCuts csNew ;
    15841668    if (rowCuts_) {
    1585 #if 0
    1586 /*
    1587   As for mode = 1 or 2; this bit of code migrated up to the head of the Stage 1
    1588   main processing loop.
    1589 */
    1590       // Only look at short rows
    1591       for (i=0;i<nRows;i++) {
    1592         if (rowLength[i]>maxElements)
    1593           abort(); //markR[i]=-2 ;
    1594       }
    1595 #endif
    15961669      ninfeas = probeCliques(si,debugger,csNew,colLower,colUpper,
    15971670                             rowCopy,columnCopy,
    15981671                             realRows,rowLower,rowUpper,
    1599                              intVar,minR,maxR,markR,info) ;
     1672                             intVar,minR,maxR,markR,info,
     1673                             useObj,useCutoff,cutoff) ;
    16001674    }
    16011675/*
     
    23132387      (*oneCut)->print() ;
    23142388# endif
    2315   } else {
     2389  } else if (si.getRowCutDebuggerAlways()) {
    23162390    std::cout << "  !! Not on optimal path." << std::endl ;
    23172391  }
     
    23222396
    23232397
     2398/*
     2399  Parameters useObj, useCutoff, and cutoff are unused as of 110208. Added to
     2400  probe(), added here for symmetry.
     2401*/
    23242402
    23252403// Does probing and adding cuts
     
    23362414                       char * intVar, double * minR, double * maxR,
    23372415                       int * markR,
    2338                        CglTreeInfo * info) const
     2416                       CglTreeInfo * info,
     2417                       bool useObj, bool useCutoff, double cutoff) const
    23392418{
    23402419  // Set up maxes
     
    24322511
    24332512  int ipass=0,nfixed=-1;
    2434 
     2513/*
     2514  Chop now that useCutoff and cutoff come in as parameters.
    24352515  double cutoff;
    24362516  bool cutoff_available = si.getDblParam(OsiDualObjectiveLimit,cutoff);
     
    24412521  if (fabs(cutoff)>1.0e30)
    24422522    assert (cutoff>1.0e30);
     2523*/
    24432524  double current = si.getObjValue();
    24442525  // make irrelevant if mode is 0
     
    33983479                      if (realRows&&realRow>0)
    33993480                        realRow=realRows[realRow];
    3400                       rowCut.addCutIfNotDuplicate(rc,realRow);
     3481                      rc.setWhichRow(realRow) ;
     3482                      rowCut.addCutIfNotDuplicate(rc);
    34013483                      }
    34023484                    }
     
    34773559                      if (realRows&&realRow>0)
    34783560                        realRow=realRows[realRow];
    3479                       rowCut.addCutIfNotDuplicate(rc,realRow);
     3561                      rc.setWhichRow(realRow) ;
     3562                      rowCut.addCutIfNotDuplicate(rc);
    34803563                      }
    34813564                    }
     
    37243807                        //printf("c point to row %d\n",irow);
    37253808#ifdef STRENGTHEN_PRINT
    3726                       if (rowLower[irow]<-1.0e20) {
    3727                         printf("7Cut %g <= ",rc.lb());
    3728                         int k;
    3729                         for ( k=0;k<n;k++) {
    3730                           int iColumn = index[k];
    3731                           printf("%g*",element[k]);
    3732                           if (si.isInteger(iColumn))
    3733                             printf("i%d ",iColumn);
    3734                           else
    3735                             printf("x%d ",iColumn);
     3809                        if (rowLower[irow]<-1.0e20) {
     3810                          printf("7Cut %g <= ",rc.lb());
     3811                          int k;
     3812                          for ( k=0;k<n;k++) {
     3813                            int iColumn = index[k];
     3814                            printf("%g*",element[k]);
     3815                            if (si.isInteger(iColumn))
     3816                              printf("i%d ",iColumn);
     3817                            else
     3818                              printf("x%d ",iColumn);
     3819                          }
     3820                          printf("<= %g\n",rc.ub());
     3821                          printf("Row %g <= ",rowLower[irow]);
     3822                          for (k=rowStart[irow];k<rowStart[irow]+rowLength[irow];k++) {
     3823                            int iColumn = column[k];
     3824                            printf("%g*",rowElements[k]);
     3825                            if (si.isInteger(iColumn))
     3826                              printf("i%d ",iColumn);
     3827                            else
     3828                              printf("x%d ",iColumn);
     3829                          }
     3830                          printf("<= %g\n",rowUpper[irow]);
    37363831                        }
    3737                         printf("<= %g\n",rc.ub());
    3738                         printf("Row %g <= ",rowLower[irow]);
    3739                         for (k=rowStart[irow];k<rowStart[irow]+rowLength[irow];k++) {
    3740                           int iColumn = column[k];
    3741                           printf("%g*",rowElements[k]);
    3742                           if (si.isInteger(iColumn))
    3743                             printf("i%d ",iColumn);
    3744                           else
    3745                             printf("x%d ",iColumn);
    3746                         }
    3747                         printf("<= %g\n",rowUpper[irow]);
    3748                       }
    3749 #endif
    3750                       int realRow = (rowLower[irow]<-1.0e20) ? irow : -1;
    3751                       if (realRows&&realRow>0)
    3752                         realRow=realRows[realRow];
    3753                         rowCut.addCutIfNotDuplicate(rc,realRow);
     3832  #endif
     3833                        int realRow = (rowLower[irow]<-1.0e20) ? irow : -1;
     3834                        if (realRows&&realRow>0)
     3835                          realRow=realRows[realRow];
     3836                        rc.setWhichRow(realRow) ;
     3837                        rowCut.addCutIfNotDuplicate(rc);
    37543838                      }
    37553839                    }
     
    38033887                        //printf("d point to row %d\n",irow);
    38043888#ifdef STRENGTHEN_PRINT
    3805                       if (rowUpper[irow]>1.0e20) {
    3806                         printf("8Cut %g <= ",rc.lb());
    3807                         int k;
    3808                         for ( k=0;k<n;k++) {
    3809                           int iColumn = index[k];
    3810                           printf("%g*",element[k]);
    3811                           if (si.isInteger(iColumn))
    3812                             printf("i%d ",iColumn);
    3813                           else
    3814                             printf("x%d ",iColumn);
     3889                        if (rowUpper[irow]>1.0e20) {
     3890                          printf("8Cut %g <= ",rc.lb());
     3891                          int k;
     3892                          for ( k=0;k<n;k++) {
     3893                            int iColumn = index[k];
     3894                            printf("%g*",element[k]);
     3895                            if (si.isInteger(iColumn))
     3896                              printf("i%d ",iColumn);
     3897                            else
     3898                              printf("x%d ",iColumn);
     3899                          }
     3900                          printf("<= %g\n",rc.ub());
     3901                          printf("Row %g <= ",rowLower[irow]);
     3902                          for (k=rowStart[irow];k<rowStart[irow]+rowLength[irow];k++) {
     3903                            int iColumn = column[k];
     3904                            printf("%g*",rowElements[k]);
     3905                            if (si.isInteger(iColumn))
     3906                              printf("i%d ",iColumn);
     3907                            else
     3908                              printf("x%d ",iColumn);
     3909                          }
     3910                          printf("<= %g\n",rowUpper[irow]);
    38153911                        }
    3816                         printf("<= %g\n",rc.ub());
    3817                         printf("Row %g <= ",rowLower[irow]);
    3818                         for (k=rowStart[irow];k<rowStart[irow]+rowLength[irow];k++) {
    3819                           int iColumn = column[k];
    3820                           printf("%g*",rowElements[k]);
    3821                           if (si.isInteger(iColumn))
    3822                             printf("i%d ",iColumn);
    3823                           else
    3824                             printf("x%d ",iColumn);
    3825                         }
    3826                         printf("<= %g\n",rowUpper[irow]);
    3827                       }
    3828 #endif
    3829                       int realRow = (rowUpper[irow]>1.0e20) ? irow : -1;
    3830                       if (realRows&&realRow>0)
    3831                         realRow=realRows[realRow];
    3832                         rowCut.addCutIfNotDuplicate(rc,realRow);
     3912  #endif
     3913                        int realRow = (rowUpper[irow]>1.0e20) ? irow : -1;
     3914                        if (realRows&&realRow>0)
     3915                          realRow=realRows[realRow];
     3916                        rc.setWhichRow(realRow) ;
     3917                        rowCut.addCutIfNotDuplicate(rc);
    38333918                      }
    38343919                    }
     
    38863971                          CglTreeInfo * info) const
    38873972{
     3973  // Check if this is ever used.
     3974  assert(false) ;
    38883975  if (!numberCliques_)
    38893976    return 0;
     
    47064793                      }
    47074794#endif
    4708                       rowCut.addCutIfNotDuplicate(rc,irow);
     4795                      rc.setWhichRow(irow) ;
     4796                      rowCut.addCutIfNotDuplicate(rc);
    47094797                    }
    47104798                  }
     
    47824870                      }
    47834871#endif
    4784                       rowCut.addCutIfNotDuplicate(rc,irow);
     4872                      rc.setWhichRow(irow) ;
     4873                      rowCut.addCutIfNotDuplicate(rc);
    47854874                    }
    47864875                  }
     
    49485037                      }
    49495038#endif
    4950                         rowCut.addCutIfNotDuplicate(rc,irow);
     5039                        rc.setWhichRow(irow) ;
     5040                        rowCut.addCutIfNotDuplicate(rc);
    49515041                      }
    49525042                    }
     
    50245114                      }
    50255115#endif
    5026                         rowCut.addCutIfNotDuplicate(rc,irow);
     5116                        rc.setWhichRow(irow) ;
     5117                        rowCut.addCutIfNotDuplicate(rc);
    50275118                      }
    50285119                    }
  • branches/CglWorking101215/src/CglProbing/CglProbing.hpp

    r963 r964  
    1414      Can be used for building cliques
    1515
    16       ?!?? What is this?
     16      ?!?? What is this?  -- lh, 10???? --
    1717   */
    1818  typedef struct {
     
    4747     returning the information. This is a big TODO.
    4848
     49     \todo
     50     Actually, it *is* done in the code here, labelled as two-cuts, in
     51     conjunction with CglPreProcess. It's just that John never bothered
     52     to change the comment. Update this once I'm a bit further along in
     53     the rewrite. -- lh, 110201 --
     54
    4955  3) x -> 1 causes a constraint to become slack by some amount c.  We can
    5056     tighten the constraint ax + .... <= b (where a may be zero) to
    5157     (a+c)x + .... <= b.  If this cut is violated then it is generated.
    5258
    53   4) Similarly, we can generate implied disaggregation cuts
     59  4) Similarly, we can generate implied disaggregation cuts.
    5460
    5561  Note the following differences from cuts in OSL:
     
    5864  b) The "chaining" in 2) was implemented.
    5965  c) Row cuts modified the original constraint rather than adding a cut.
    60   b) This code can cope with general integer variables.
     66  d) This code can cope with general integer variables. Sort of. Sometimes
     67     it trails off into a disabled dead end.
     68
     69  \todo
     70  c) is also untrue for the existing code, in the sense that CglPreProcess
     71  will replace original rows with coefficient strengthening cuts. Again,
     72  the comment was never updated. Update once I'm a bit further along in
     73  the rewrite. -- lh, 110201 --
    6174
    6275  The generated cuts are inserted into an OsiCut object. Both row cuts and
     
    348361            const int *realRow, const double *rowLower, const double *rowUpper,
    349362            const char *intVar, double *minR, double *maxR, int *markR,
    350             CglTreeInfo *info) const;
     363            CglTreeInfo *info,
     364            bool useObj, bool useCutoff, double cutoff) const;
    351365  /// Does probing and adding cuts (with cliques)
    352366  int probeCliques( const OsiSolverInterface & si,
     
    357371             double * rowLower, double * rowUpper,
    358372             char * intVar, double * minR, double * maxR, int * markR,
    359              CglTreeInfo * info) const;
     373             CglTreeInfo * info,
     374             bool useObj, bool useCutoff, double cutoff) const;
    360375  /// Does probing and adding cuts for clique slacks
    361376  int probeSlacks( const OsiSolverInterface & si,
     
    508523
    509524  /*! \brief Row cuts flags
     525
     526    Broadly, the possibilities are:
     527    - generate row cuts that will be returned in an OsiCuts collection;
     528    - generate row cuts that can replace (strengthen) rows in the original
     529      constraint system, returned through a CglTreeInfo structure;
     530    - generate column cuts
     531    If the CglTreeInfo structure contains a strengthenRow vector, and you
     532    specify 0x4 here, no row cuts are returned in the cut collection, but rows
     533    are strengthened if possible.
    510534
    511535    - 0x0: do not generate row cuts
  • branches/CglWorking101215/src/CglProbing/CglProbingDebug.cpp

    r963 r964  
    7575
    7676/*
    77   Utility to scan a packed matrix for damage. Not exhaustive, but occasionally
    78   useful. Assumes column-major.
    79 */
    80 
    81 int verifyMtx (const CoinPackedMatrix *const mtx)
    82 
    83 {
    84   const double smallCoeff = 1.0e-50 ;
    85   const double largeCoeff = 1.0e50 ;
    86 
    87   int n = mtx->getNumCols() ;
    88   int m = mtx->getNumRows() ;
    89   int e = mtx->getNumElements() ;
    90 
    91   int majLen, minLen ;
    92   std::string majName, minName ;
    93 
    94   if (mtx->isColOrdered()) {
    95     majLen = n ;
    96     majName = "col" ;
    97     minLen = m ;
    98     minName = "row" ;
    99   } else {
    100     majLen = m ;
    101     majName = "row" ;
    102     minLen = n ;
    103     minName = "col" ;
    104   }
    105 
    106   double extraGap = mtx->getExtraGap() ;
    107   double extraMajor = mtx->getExtraMajor() ;
    108 
    109   std::cout
    110     << " Matrix is " << ((mtx->isColOrdered())?"column":"row") << "-major, "
    111     << m << " rows X " << n << " cols; ex maj " << extraMajor
    112     << ", ex gap " << extraGap << "." << std::endl ;
    113   if (mtx->hasGaps())
    114     std::cout << "  matrix has gaps." << std::endl ;
    115 
    116   const CoinBigIndex *const majStarts = mtx->getVectorStarts() ;
    117   const int *const majLens = mtx->getVectorLengths() ;
    118   const int *const minInds = mtx->getIndices() ;
    119   const double *const coeffs = mtx->getElements() ;
    120 /*
    121   Scan the major dimension and make sure that the minor dimension indices are
    122   reasonable.
    123 */
    124   int errs = 0 ;
    125 
    126   for (int majndx = 0 ; majndx < majLen ; majndx++) {
    127     CoinBigIndex majStart = majStarts[majndx] ;
    128     CoinBigIndex majStartp1 = majStarts[majndx+1] ;
    129     int majLen = majLens[majndx] ;
    130 
    131     if (majLen != (majStartp1-majStart)) {
    132       std::cout
    133         << "  " << majName << majndx
    134         << ": (end+1)-start = " << majStartp1 << " - " << majStart
    135         << " = " << majStartp1-majStart << " != len = " << majLen
    136         << "." << std::endl ;
    137       // errs++ ;
    138     }
    139 
    140     for (CoinBigIndex ii = majStart ;  ii < majLen ; ii++) {
    141       int minndx = minInds[ii] ;
    142       if (minndx < 0 || minndx >= minLen) {
    143         std::cout
    144           << "  " << majName << " " << majndx << ": "
    145           << minName << " index " << ii << " is " << minndx
    146           << ", should be between 0 and " << minLen-1 << "." << std::endl ;
    147         errs++ ;
    148       }
    149       double aij = coeffs[ii] ;
    150       if (CoinIsnan(aij) ||
    151           CoinAbs(aij) > largeCoeff || CoinAbs(aij) < smallCoeff) {
    152         std::cout
    153           << "  a<" << majndx << "," << minndx << "> = " << aij
    154           << " appears bogus." << std::endl ;
    155         errs++ ;
    156       }
    157     }
    158   }
    159   if (errs > 0) {
    160     std::cout << "  Detected " << errs << " errors in matrix." << std::endl ;
    161   }
    162 
    163   return (errs) ;
    164 }
    165 
    166 /*
    16777  Debugging utility to print row information. Far too many parameters, but
    16878  doing it right requires too much code to replicate at multiple locations.
  • branches/CglWorking101215/src/CglProbing/CglProbingDebug.hpp

    r963 r964  
    3636
    3737namespace CglProbingDebug {
    38 
    39 /*! \brief Verify integrity of a column-major matrix.
    40 
    41   Scans the index and coefficient arrays looking for obviously bogus values.
    42   Can be handy when trying to pin down where a matrix is getting corrupted.
    43 */
    44 int verifyMtx (const CoinPackedMatrix *const mtx) ;
    4538
    4639/*! \brief Check column bounds against optimal solution
  • branches/CglWorking101215/src/CglProbing/CglProbingProbe.cpp

    r963 r964  
    3030  Allocating one big block for data arrays is efficient but makes it
    3131  difficult for debugging tools to detect accesses outside the boundary of an
    32   array. The assert checks that a double is larger than an int and that the
    33   ratio is a power of two.
     32  array.
    3433*/
    3534#if CGL_DEBUG == 0
     
    8079
    8180// ===============================================
     81
     82/*
     83  This method looks for implication cuts: implications that can be generated
     84  given a binary probing variable x<p>. For an arbitrary binary variable x<j>,
     85  we're looking for things like the following:
     86    1) x<p> -> 0 ==> x<j> -> 0 and x<p> -> 1 ==> x<j> -> 0
     87       In this case, we can fix x<j> to 0. Similarly if were driven to 1.
     88    2) x<p> -> 0 ==> x<j> -> 0 and x<p> -> 1 ==> x<j> -> 1
     89       In this case, we can write x<j> = x<p>. If the sense is reversed, we
     90       get x<j> = (1-x<p>).
     91  Form 2) are actual cuts; form 1) just amounts to a bound change.
     92
     93  You can extend this to arbitrary bounds on x<j>. Let ld<j> be the lower
     94  bound for the down probe and ud<j> be the upper bound for the down probe.
     95  Similarly, lu<j> and uu<j> are lower and upper bounds from the up probe, and
     96  l<j> and u<j> are the original bounds.
     97    1) The equivalent is to look for max(min(ld<j>,lu<j>),l<j>) and
     98       min(max(ud<j>,uu<j>),u<j>).
     99    2) The equivalent constraints are x<j> = (u<j>-l<j>)x<p> + l<j> and
     100       x<j> = -(u<j>-l<j>)x<p> + u<j>
     101  As with binary variables, form 2) are actual cuts; form 1) just amounts to
     102  a bound change.
     103
     104  1) is applied to continuous variables as well as integer variables,
     105  and the result is captured in vec_l and vec_u. Normally, these
     106  parameters will be saveL and saveU, hence the result will propagate
     107  out from CglProbing::genCutsAndModify when the `original' bounds are
     108  restored. Only integer variables are recorded as column cuts.
     109
     110  The method uses jProbe to decide if it should generate form 2) cuts. If
     111  jProbe >= 0, it's assumed to be a binary variable and form 2) cuts are
     112  attempted.
     113
     114  Bounds held in vec_lu and vec_uu are assumed to be indexed by column
     115  number.  Bounds held in vec_l, and vec_u are assumed to be correlated
     116  with entries on stackC. Bounds held in vec_ld and vec_ud are assumed to
     117  be correlated with entries on stackC0.
     118
     119  In particular, stackC[0] and stackC0[0] are the probing variable.
     120
     121  MarkC, element, and index are used as work arrays. It's relatively rare that
     122  any given probe will change bounds on half the variables, so we're going to
     123  bet on that. We'll store lower bounds from [0] of index and element and
     124  upper bounds from [nCols-1].
     125
     126  TODO: Check that the above comment is correct --- do improvements to
     127        continuous bounds really propagate back to the caller?
     128        -- lh, 110201 --
     129*/
     130
     131void implicationCuts (int jProbe, double primalTol, int nCols,
     132                      OsiCuts &cs, const OsiSolverInterface &si,
     133                      const char *const intVar,
     134                      const double *const colsol,
     135                      int nC, const int *const stackC, int *const markC,
     136                      const double *const vec_uu, const double *const vec_lu,
     137                      double *const vec_u, double *const vec_l,
     138                      int nC0, const int *const stackC0,
     139                      const double *const vec_ud, const double *const vec_ld,
     140                      double *const element, int *const index
     141                     )
     142{
     143# if CGL_DEBUG > 0
     144  std::cout
     145    << "Entering implicationCuts, probe " << jProbe << "." << std::endl ;
     146  bool hitTheWall = false ;
     147  int origRowCutCnt = cs.sizeRowCuts() ;
     148  int origColCutCnt = cs.sizeColCuts() ;
     149/*
     150  If we have a binary probe variable, it should not be fixed.
     151*/
     152  if (jProbe >= 0) assert(vec_l[0] == 0.0 && vec_u[0] == 1.0) ;
     153# endif
     154/*
     155  Construct a cross-reference so that we can correlate entries on stackC0 and
     156  stackC. markC is no longer of use, so use it. Keep in mind that stackC0 and
     157  stackC will not necessarily hold the same variables; the 1000 offset
     158  guarantees we can recognise the entries made here.
     159*/
     160  for (int iC = nC-1 ; iC >= 0 ; iC--) {
     161    int j = stackC[iC] ;
     162    markC[j] = iC+1000 ;
     163  }
     164/*
     165  See if we have bounds improvement. Given a lower bound ld<j> from the
     166  down probe, a bound lu<j> from the up probe, and the original bound lo<j>,
     167  we can surely impose the lesser of ld<j> and lu<j> if this is more than the
     168  original lo<j>.  In math, max(min(ld<j>,lu<j>),lo<j>). Use a conservative
     169  tolerance. Stash the improved bound in saveL so that it'll be restored
     170  when we restore the `original' (pre-probing) bounds. For integer variables,
     171  add a column cut.
     172*/
     173  int nTot = 0 ;
     174  int nInt = 0 ;
     175  int ilb = 0 ;
     176  int iub = nCols-1 ;
     177  bool ifCut = false ;
     178  const double minChange = 1.0e-4 ;
     179  int implIndices[2] ;
     180  double implCoeffs[2] ;
     181/*
     182  Scan stackC0 looking for variables that had bound improvement on both the up
     183  and down probe. If the variable was changed by only one probe, move on.
     184*/
     185  for (int iC0 = 1 ; iC0 < nC0 ; iC0++) {
     186    int j = stackC0[iC0] ;
     187    int iC = markC[j]-1000 ;
     188    if (iC < 0) continue ;
     189    double udj = vec_ud[iC0] ;
     190    double ldj = vec_ld[iC0] ;
     191    double uuj = vec_uu[j] ;
     192    double luj = vec_lu[j] ;
     193    double &uj = vec_u[iC] ;
     194    double &lj = vec_l[iC] ;
     195/*
     196  Obtain the consensus lower bound.
     197*/
     198    double bestProbelj = CoinMin(ldj,luj) ;
     199    if (bestProbelj > lj+minChange) {
     200#     if CGL_DEBUG > 0
     201      std::cout << "      consensus lb for " ;
     202      if (intVar[j]) std::cout << "(int) " ;
     203      std::cout
     204        << " x<" << j << "> improved from "
     205        << lj << " to " << bestProbelj << "." << std::endl ;
     206#     endif
     207      lj = bestProbelj ;
     208      nTot++ ;
     209      if (intVar[j] && nInt < nCols) {
     210        element[ilb] = bestProbelj ;
     211        index[ilb++] = j ;
     212        nInt++ ;
     213        if (colsol[j] < bestProbelj-primalTol)
     214          ifCut = true ;
     215      }
     216#     if CGL_DEBUG > 0
     217      else {
     218        if (intVar[j] && nInt >= nCols && hitTheWall == false) {
     219          std::cout
     220            << "CglProbing::implicationCuts: (lb) exceeded column cut space!."
     221            << std::endl ;
     222          hitTheWall = true ;
     223        }
     224      }
     225#     endif
     226    }
     227/*
     228  Obtain the consensus upper bound.
     229*/
     230    double bestProbeuj = CoinMax(udj,uuj) ;
     231    if (bestProbeuj < uj-minChange) {
     232#     if CGL_DEBUG > 0
     233      std::cout << "      consensus ub for " ;
     234      if (intVar[j]) std::cout << "(int) " ;
     235      std::cout
     236        << " x<" << j << "> improved from "
     237        << uj << " to " << bestProbeuj << "." << std::endl ;
     238#     endif
     239      uj = bestProbeuj ;
     240      nTot++ ;
     241      if (intVar[j] && nInt < nCols) {
     242        element[iub] = bestProbeuj ;
     243        index[iub--] = j ;
     244        nInt++ ;
     245        if (colsol[j] > bestProbeuj+primalTol)
     246          ifCut = true ;
     247      }
     248#     if CGL_DEBUG > 0
     249      else {
     250        if (intVar[j] && nInt >= nCols && hitTheWall == false) {
     251          std::cout
     252            << "CglProbing::implicationCuts: (ub) exceeded column cut space!."
     253            << std::endl ;
     254          hitTheWall = true ;
     255        }
     256      }
     257#     endif
     258    }
     259/*
     260  Now look for honest implication cuts (form 2) in the introductory comments.
     261  Recall that jProbe >= 0 is taken as an indication that the probe variable
     262  is binary and these cuts should be attempted.  It's possible (miracles
     263  happen) that the above code has managed to fix x<j> by pulling together
     264  the bounds from the down and up probe (but we should not lose feasibility
     265  here). In the case that x<j> is fixed, move on.
     266*/
     267    if (jProbe >= 0) {
     268      assert(uj >= lj) ;
     269      if (CoinAbs(uj-lj) < primalTol) continue ;
     270/*
     271  Check for x<j> = (u<j>-l<j>)x<p> + l<j>. In words, check that the down probe
     272  forced x<j> to l<j> and the up probe forced x<j> to u<j>.
     273*/
     274      if (udj < lj+primalTol && luj > uj-primalTol) {
     275        OsiRowCut rc ;
     276        rc.setLb(lj) ;
     277        rc.setUb(lj) ;
     278        rc.setEffectiveness(1.0e-5) ;
     279        implIndices[0] = jProbe ;
     280        implIndices[1] = j ;
     281        implCoeffs[0] = -(uj-lj) ;
     282        implCoeffs[1] = 1.0 ;
     283        rc.setRow(2,implIndices,implCoeffs,false) ;
     284#       if CGL_DEBUG > 0
     285        std::cout << "      implication:" ;
     286        rc.print() ;
     287#       endif
     288        cs.insert(rc) ;
     289      } else if (uuj < lj+primalTol && ldj > uj-primalTol) {
     290/*
     291  Repeat for x<j> = -(u<j>-l<j>)x<p> + u<j>
     292*/
     293        OsiRowCut rc ;
     294        rc.setLb(uj) ;
     295        rc.setUb(uj) ;
     296        rc.setEffectiveness(1.0e-5) ;
     297        implIndices[0] = jProbe ;
     298        implIndices[1] = j ;
     299        implCoeffs[0] = uj-lj ;
     300        implCoeffs[1] = 1.0 ;
     301        rc.setRow(2,implIndices,implCoeffs,false) ;
     302#       if CGL_DEBUG > 0
     303        std::cout << "      inverse implication:" ;
     304        rc.print() ;
     305#       endif
     306        cs.insert(rc) ;
     307      }
     308    }
     309  }
     310/*
     311  If we have column cuts, load them into the cut collection passed in as a
     312  parameter.
     313*/
     314  if (nInt > 0) {
     315    OsiColCut cc ;
     316    if (ilb > 0) {
     317      cc.setLbs(ilb,index,element) ;
     318    }
     319    if (iub < nCols-1) {
     320      cc.setUbs((nCols-1)-iub,&index[iub+1],&element[iub+1]) ;
     321    }
     322    if (ifCut) {
     323      cc.setEffectiveness(100.0) ;
     324    } else {
     325      cc.setEffectiveness(1.0e-5) ;
     326    }
     327#   if CGL_DEBUG > 0
     328    std::cout
     329      << "    tightened " << ilb << " lb, " << (nCols-1)-iub << " ub" ;
     330    if (ifCut)
     331      std::cout << "; cut." << std::endl ;
     332    else
     333      std::cout << "; no cut." << std::endl ;
     334    CglProbingDebug::checkBounds(si,cc) ;
     335#   endif
     336    cs.insert(cc) ;
     337  }
     338# if CGL_DEBUG > 0
     339  std::cout
     340    << "Leaving implicationCuts, " << cs.sizeRowCuts()-origRowCutCnt
     341    << " row cuts, " << cs.sizeColCuts()-origColCutCnt << " col cuts."
     342    << std::endl ;
     343# endif
     344
     345  return ;
     346}
    82347
    83348
     
    142407    << "Entering disaggCuts, "
    143408    << ((probeDir == probeDown)?"down":"up") << " probe on "
    144     << si.getColName(pndx) << "<" << pndx << ">." << std::endl ;
     409    << "x<" << pndx << ">." << std::endl ;
    145410  const OsiRowCutDebugger *debugger = si.getRowCutDebugger() ;
    146411  int cutsAtStart = rowCut.numberCuts() ;
     
    351616    double delta = value*movement ;
    352617/*
    353   markR is initialised by tighten2. A code of -1 indicates this row should be
    354   processed; -2 says ignore. Other codes should not happen.
     618  markR is initialised by calcRowBounds. A code of -1 indicates this row
     619  should be processed; -2 says ignore. Other codes should not happen.
    355620
    356621  TODO: This assert, and the disabled continue, go back to the change in
    357         tighten2 where rows with infinite bounds are handed arbitrary
     622        calcRowBounds where rows with infinite bounds are handed arbitrary
    358623        bounds of 1e60 and labelled suitable for further processing. The
    359624        assert should be removed and the continue clause reinstated.
     
    460725                double *const largestPositiveInRow)
    461726{
     727# if CGL_DEBUG > 0
     728  std::cout << "Entering groomSoln." << std::endl ;
     729  std::cout << "  verifying matrix." ;
     730  rowCopy->verifyMtx(3) ;
     731# endif
    462732  int nRows = rowCopy->getNumRows() ;
    463733  int nCols = rowCopy->getNumCols() ;
     
    483753      if (value > 0.0)
    484754        break ;
     755    }
     756    if (rowStartPos[i] != kk) {
     757      std::cout
     758        << "  Row " << i << ": " << kkstart << ".." << kkend << ", len "
     759        << rowLength[i] << ", pos at " << rowStartPos[i] << "." << std::endl ;
    485760    }
    486761    assert (rowStartPos[i] == kk) ;
     
    526801    columnGap[i] = gap-primalTolerance ;
    527802  }
    528 
     803# if CGL_DEBUG > 0
     804  std::cout << "Leaving groomSoln." << std::endl ;
     805# endif
    529806  return ;
    530807}
     
    565842  int probedVar = stackC[0] ;
    566843  std::cout << "Entering monotoneActions." ;
    567 # if CGL_DEBUG > 1
     844# if CGL_DEBUG > 0
    568845  std::cout
    569846    << std::endl
    570     << "    " << si.getColName(probedVar) << " (" << probedVar
     847    << "    " << " x(" << probedVar
    571848    << ") monotone [" << colLower[probedVar] << ", "
    572849    << colUpper[probedVar] << "] from " << colsol[probedVar]
     
    589866          anyColumnCuts++ ;
    590867        }
     868#       if CGL_DEBUG > 0
     869        std::cout
     870          << "    " << " x(" << icol
     871          << ") ub " << origColUpper[icol] << " ==> " << colUpper[icol] ;
     872        if (ifCut) std::cout << " (cut)" ;
     873        std::cout << "." << std::endl ;
     874#       endif
    591875      }
    592876    }
     
    607891          anyColumnCuts++ ;
    608892        }
     893#       if CGL_DEBUG > 0
     894        std::cout
     895          << "    " << " x(" << icol
     896          << ") lb " << origColLower[icol] << " ==> " << colLower[icol] ;
     897        if (ifCut) std::cout << " (cut)" ;
     898        std::cout << "." << std::endl ;
     899#       endif
    609900      }
    610901    }
     
    628919# if CGL_DEBUG > 0
    629920  std::cout << "  Leaving monotoneActions" ;
     921  if (nTot > 0)
     922    std::cout << ", " << nTot << " monotone" ;
    630923  if (anyColumnCuts > 0)
    631924    std::cout << ", " << anyColumnCuts << " column cuts" ;
     
    687980    against L<i> ==> a<ij>*c<j>*dir > 0
    688981
     982  Note that we cannot process the objective here --- that'd be nonsense.
     983  It's a minor pain to exclude it.
     984
    689985  John's original comment for this code was
    690986    // also see if singletons can go to good objective
     
    696992*/
    697993
    698 void singletonRows (int jProbe, double primalTolerance_,
     994void singletonRows (int jProbe, double primalTolerance_, bool useObj,
    699995                    const OsiSolverInterface &si,
    700996                    const CoinPackedMatrix *rowCopy,
     
    7211017  const double *rowElements = rowCopy->getElements() ;
    7221018  const int nCols = rowCopy->getNumCols() ;
     1019  const int nRows = rowCopy->getNumRows() ;
    7231020/*
    7241021  `Singleton' must be based on the column length in the original system.
     
    7281025  const double objSense = si.getObjSense() ;
    7291026/*
    730   Open a loop to work through the rows on stackR.
     1027  Open a loop to work through the rows on stackR. Don't process the objective
     1028  (if present, it's the last row).
    7311029*/
    7321030  for (int istackR = 0 ; istackR < nstackR ; istackR++) {
    7331031    int i = stackR[istackR] ;
     1032    if (useObj && i == (nRows-1)) continue ;
    7341033/*
    7351034  Check the gaps. If the constraint is potentially tight in both directions,
     
    7791078  be naturally driven to l<j>.
    7801079 
    781   Don't overwrite the saved bounds if we've
    782   tightened this variable already!
     1080  Don't overwrite the saved bounds if we've tightened this variable already!
    7831081*/
    7841082    if (uGap > primalTolerance_) {
     
    7901088          double value = rowElements[kk] ;
    7911089          if (objSense*objective[j]*value < 0.0)
    792           { assert(jProbe != j) ;
    793             if (!(markC[j]&(tightenLower|tightenUpper))) {
     1090          { if (!(markC[j]&(tightenLower|tightenUpper))) {
     1091              assert(nstackC < nCols) ;
    7941092              stackC[nstackC] = j ;
    7951093              saveL[nstackC] = lj ;
    7961094              saveU[nstackC] = uj ;
     1095              nstackC++ ;
    7971096            }
    798             assert(saveU[nstackC] > saveL[nstackC]) ;
    7991097            if (value > 0.0) {
    8001098              colLower[j] = uj ;
     
    8041102            markC[j] |= tightenLower|tightenUpper ;
    8051103            colGap[j] = -primalTolerance_ ;
    806             assert(nstackC < nCols) ;
    807             nstackC++ ;
    8081104#           if CGL_DEBUG > 1
    8091105            std::cout
    8101106              << "    row " << i << " uGap " << rowUpper[i] << "-"
    8111107              << maxR[i] << " = " << uGap
    812               << " " << si.getColName(j) << "(" << j
     1108              << " " << " x(" << j
    8131109              << ") c = " << objective[j] << ", a = " << value ;
    8141110            if (colLower[j] == uj)
    815               std::cout << " l " << saveL[nstackC-1] << " -> " << uj ;
     1111              std::cout << " l " << lj << " -> " << uj ;
    8161112            else
    817               std::cout << " u " << saveU[nstackC-1] << " -> " << lj ;
     1113              std::cout << " u " << uj << " -> " << lj ;
    8181114            std::cout << std::endl ;
    8191115#           endif
     
    8361132          double value = rowElements[kk] ;
    8371133          if (objSense*objective[j]*value > 0.0)
    838           { assert(jProbe != j) ;
    839             if (!(markC[j]&(tightenLower|tightenUpper))) {
     1134          { if (!(markC[j]&(tightenLower|tightenUpper))) {
     1135              assert(nstackC < nCols) ;
    8401136              stackC[nstackC] = j ;
    8411137              saveL[nstackC] = lj ;
    8421138              saveU[nstackC] = uj ;
     1139              nstackC++ ;
    8431140            }
    844             assert(saveU[nstackC] > saveL[nstackC]) ;
    8451141            if (value < 0.0) {
    8461142              colLower[j] = uj ;
     
    8501146            markC[j] |= tightenLower|tightenUpper ;
    8511147            colGap[j] = -primalTolerance_ ;
    852             assert(nstackC < nCols) ;
    853             nstackC++ ;
    8541148#           if CGL_DEBUG > 1
    8551149            std::cout
    8561150              << "    row " << i << " lGap " << minR[i] << "-"
    8571151              << rowLower[i] << " = " << lGap
    858               << " " << si.getColName(j) << "(" << j
     1152              << " " << " x(" << j
    8591153              << ") c = " << objective[j] << ", a = " << value ;
    8601154            if (colLower[j] == uj)
    861               std::cout << " l " << saveL[nstackC-1] << " -> " << uj ;
     1155              std::cout << " l " << lj << " -> " << uj ;
    8621156            else
    863               std::cout << " u " << saveU[nstackC-1] << " -> " << lj ;
     1157              std::cout << " u " << uj << " -> " << lj ;
    8641158            std::cout << std::endl ;
    8651159#           endif
     
    9221216                      int jProbe, unsigned int probeDir,
    9231217                      double primalTolerance_,
    924                       bool justReplace, bool canReplace,
     1218                      bool justReplace,
    9251219                      double needEffectiveness,
    9261220                      const OsiSolverInterface &si,
     
    9431237  std::cout
    9441238    << "Entering strengthenCoeff, probing "
    945     << si.getColName(jProbe) << "(" << jProbe << ") "
     1239    << "x<" << jProbe << "> "
    9461240    << ((probeDir == probeDown)?"down":"up") << "." << std::endl ;
    9471241  const OsiRowCutDebugger *debugger = si.getRowCutDebugger() ;
    9481242  int newCuts = 0 ;
     1243  int cutsAtStart = rowCut.numberCuts() ;
    9491244# endif
    9501245/*
     
    9771272  We can't get anywhere unless probing has made one end or the other of the
    9781273  constraint redundant. If there's no gap, we can't do anything.
     1274
     1275  Range constraints pose the danger that we'll want to assign two different
     1276  values to a<ip>. If we're in `justReplace' mode, this is a nonstarter.
    9791277*/
    9801278    double uGap = bi-maxR[irow] ;
     
    9821280    if (uGap < primalTolerance_ && -lGap < primalTolerance_) continue ;
    9831281    bool isRangeCon = ((blowi  > -rhsInf) && (bi < rhsInf)) ;
     1282    if (isRangeCon && justReplace) continue ;
    9841283/*
    9851284  We'll need the lhs activity, excluding the probed variable, to evaluate the
     
    10391338  the original rhs and the revised lhs bound.
    10401339
    1041   We'll also generate the cut if we can strengthen it in place and we're not
    1042   dealing with a range constraint. (The reason for excluding range constraints
    1043   is that we might try to alter a<ip> against both b<i> and blow<i>. That way
    1044   lies madness.)
     1340  We'll also generate the cut if we can strengthen it in place.
    10451341
    10461342  TODO: It's not clear to me how the first two effectiveness calculations
     
    10711367      effectiveness = CoinMax(effectiveness,aipZeroEffectiveness) ;
    10721368    if (effectiveness > needEffectiveness) genCut = true ;
    1073     if (info->strengthenRow && !isRangeCon) genCut = true ;
     1369    if (info->strengthenRow) genCut = true ;
    10741370/*
    10751371  Are we going to generate the cut? If not, on to the next iteration.
     
    10771373    if (!genCut) continue ;
    10781374/*
    1079   Generate the coefficients. Copy whatever's present and plug in aipPrime at
    1080   the end.
     1375  Generate the coefficients. Copy all except a<ip>. Add a<ip>' at the end, if
     1376  it's nonzero. (a<ip>' can be zero and still have a perfectly valid cut.)
    10811377*/
    10821378    int n = 0 ;
    1083     int aip_n = -1 ;
    10841379    for (CoinBigIndex kk = kkstart ; kk < kkend ; kk++) {
    10851380      int k = column[kk] ;
    10861381      double aik = rowElements[kk] ;
    1087       if (k == jProbe)
    1088         aip_n = n ;
    1089       index[n] = k ;
    1090       element[n++] = aik ;
     1382      if (k != jProbe) {
     1383        index[n] = k ;
     1384        element[n++] = aik ;
     1385      }
    10911386    }
    1092     if (aip_n < 0) {
    1093       aip_n = n ;
    1094       n++ ;
     1387    if (aipPrimeNonZero) {
     1388      index[n] = jProbe ;
     1389      element[n++] = aipPrime ;
    10951390    }
    1096     index[aip_n] = jProbe ;
    1097     element[aip_n] = aipPrime ;
    10981391/*
    10991392  Fill in a cut structure with the cut.
     
    11191412/*
    11201413  If we're in preprocessing, we might try to simply replace the existing
    1121   constraint (justReplace = true; preprocessing is the only place this will
    1122   happen). Otherwise, drop the cut into the cut set.
     1414  constraint (justReplace = true). Otherwise, drop the cut into the cut set.
    11231415
    11241416  realRows comes in as a parameter. This is the translation array created if
     
    11261418
    11271419  Effectiveness, if we're strengthening in place, seems to be absolute size of
    1128   coefficients; smaller is better. (Why?)
    1129 
    1130   TODO: It seems to me that we can get here with
    1131         justReplace = true and canReplace = false, and this condition is
    1132         constant over an iteration of the way loop. In other words, we've done
    1133         all the work to generate this cut and we knew before we started that
    1134         we would discard it here.  -- lh, 110107 --
    1135         Put in an assert to see if we ever do all the work, only to reject
    1136         the cut.  -- lh, 110115 --
    1137 */
    1138       int realRow = (canReplace && !isRangeCon)?(irow):(-1) ;
    1139       if (realRows && realRow >= 0)
     1420  coefficients; smaller is better. (Why? Proxy for fewer coefficients?)
     1421*/
     1422      int realRow = irow ;
     1423      if (realRows)
    11401424        realRow = realRows[realRow] ;
     1425      assert(realRow >= 0) ;
     1426      rc.setWhichRow(realRow) ;
    11411427      if (!justReplace) {
    11421428#       if CGL_DEBUG > 0
    11431429          std::cout
    1144             << "    Adding cut, real row " << realRow << "." << std::endl ;
     1430            << "    strengthen coeff cut on real row "
     1431            << realRow << " added to cut collection." << std::endl ;
    11451432          newCuts++ ;
    11461433#       endif
    1147         rowCut.addCutIfNotDuplicate(rc,realRow) ;
    1148       } else if (realRow >= 0) {
     1434        rowCut.addCutIfNotDuplicate(rc) ;
     1435      } else {
    11491436        double effectiveness = 0.0 ;
    11501437        for (int i = 0 ; i < n ; i++)
     
    11561443#         if CGL_DEBUG > 0
    11571444          std::cout
    1158             << "    Cut will replace row " << realRow << "." << std::endl ;
     1445            << "    strengthen coeff on real row " << realRow
     1446            << " will replace row, eff " << effectiveness << "."
     1447            << std::endl ;
    11591448          newCuts++ ;
    11601449#         endif
    11611450          info->strengthenRow[realRow] = rc.clone() ;
    11621451        } else {
    1163 #       if CGL_DEBUG > 0
    1164           std::cout
    1165             << "INPLACE: rejected on cut coeff magnitude." << std::endl ;
    1166 #       endif
     1452#         if CGL_DEBUG > 0
     1453            std::cout
     1454              << "    strengthen coeff cut on real row " << realRow
     1455              << " rejected; eff " << effectiveness << " >  eff "
     1456              << info->strengthenRow[realRow]->effectiveness()
     1457              << " of incumbent." << std::endl ;
     1458#         endif
    11671459        }
    1168       } else {
    1169 #       if CGL_DEBUG > 0
    1170           std::cout
    1171             << "INPLACE: rejected because no real row." << std::endl ;
    1172 #       endif
    11731460      }
    11741461    }
    11751462
    11761463# if CGL_DEBUG > 0
     1464  int cutsAtEnd = rowCut.numberCuts() ;
     1465  int outplaceCuts = cutsAtEnd-cutsAtStart ;
     1466  int inplaceCuts = newCuts-outplaceCuts ;
    11771467  std::cout << "Leaving strengthenCoeff" ;
    11781468  if (newCuts > 0)
    1179     std::cout << ", " << newCuts << " cuts" ;
     1469    std::cout << ", " << newCuts << " cuts, " << inplaceCuts << " in place" ;
    11801470  std::cout << "." << std::endl ;
    11811471# endif
     
    12941584              time we probe pushing a variable in a given direction.
    12951585
    1296   minR, maxR, markR  Initialised externally by tighten2. For row i, minR[i] =
    1297                      LB(i), maxR[i] = UB(i) (row lhs lower and upper bounds,
    1298                      respectively).  markR is set to -1 if there's at least
    1299                      one finite lhs bound, -2 if we shouldn't process this
    1300                      row (no finite bounds, or too many coefficients).
     1586  minR, maxR, markR  Initialised externally by calcRowBounds. For row i,
     1587                     minR[i] = LB(i), maxR[i] = UB(i) (row lhs lower and
     1588                     upper bounds, respectively).  markR is set to -1 if
     1589                     there's at least one finite lhs bound, -2 if we shouldn't
     1590                     process this row (no finite bounds, or too many
     1591                     coefficients).
    13011592
    13021593  Then, as we work, markR[i] will be set to point to the entry in stackR for
     
    13261617                       const char *intVar, double *minR, double *maxR,
    13271618                       int *markR,
    1328                        CglTreeInfo *info) const
     1619                       CglTreeInfo *info,
     1620                       bool useObj, bool useCutoff, double cutoff) const
    13291621
    13301622{
    13311623# if CGL_DEBUG > 0
    13321624  std::cout << "Entering CglProbing::probe." << std::endl ;
     1625  int numberRowCutsBefore = cs.sizeRowCuts() ;
     1626  int numberColCutsBefore = cs.sizeColCuts() ;
    13331627# endif
    13341628/*
     
    14111705  root. In the tree, much less.
    14121706
    1413   If all we're doing is strengthening rows in place (option 0x40 indicates
    1414   we're in preprocessing, a good place to do this), we can only work with
    1415   what we have in rowCopy (and we need a translation array).
    1416 
    1417   TODO: There's a problem with realRows here --- if we didn't delete any
    1418         rows during preprocessing in gutsOfGenerateCuts, realRows is not
    1419         allocated. So we can't strengthen in place in that case. Not likely
    1420         the intent.  -- lh, 101209 --
    1421 */
    1422   bool justReplace = ((info->options&0x40) != 0) && (realRows != NULL) ;
    1423   int nRowsFake = 0 ;
     1707  If what we're doing is strengthening rows in place (option 0x40 indicates
     1708  we're in preprocessing, a good place to do this), we can generate at most
     1709  one cut per original row (nRows) and we need the vector strengthenRow
     1710  to pass the cut back to the caller.
     1711
     1712  TODO: I'm asking myself `Why do we need a cut container at all, if we're
     1713        in `just replace' mode?'. The cuts are kept in strengthenRow. But
     1714        at the end we copy them over into cs anyhow. This can probably be
     1715        rationalised.  -- lh, 110209 --
     1716*/
     1717  bool justReplace =
     1718      ((info->options&0x40) != 0) && (info->strengthenRow != NULL) ;
     1719  int cutCapacity = 0 ;
    14241720  if (justReplace) {
    1425     nRowsFake = nRows ;
     1721    cutCapacity = nRows ;
    14261722  } else {
    14271723    int nRowsSafe = CoinMin(nRows,si.getNumRows()) ;
    1428     nRowsFake = info->inTree ? nRowsSafe/3 : nRowsSafe*10 ;
    1429     if (!info->inTree && info->pass == 0) nRowsFake *= 10 ;
     1724    cutCapacity = info->inTree ? nRowsSafe/3 : nRowsSafe*10 ;
     1725    if (!info->inTree && info->pass == 0) cutCapacity *= 10 ;
    14301726  }
    1431   CglProbingRowCut rowCut(nRowsFake,!info->inTree) ;
     1727  CglProbingRowCut rowCut(cutCapacity,!info->inTree) ;
    14321728
    14331729
     
    14521748            rowCopy,rowStartPos,largestNegativeInRow,largestPositiveInRow) ;
    14531749/*
    1454   Determine objective cutoff (minimisation convention).
    1455 
    1456   TODO: Seems to me that this bit of code is a guaranteed fail for a
    1457         maximisation problem with no cutoff. It also repeats in a number
    1458         of places.  -- lh, 101125 --
    1459 */
    1460   double cutoff ;
    1461   bool cutoff_available = si.getDblParam(OsiDualObjectiveLimit,cutoff) ;
    1462   if (!cutoff_available || usingObjective_ < 0) {
    1463     cutoff = si.getInfinity() ;
    1464   }
    1465   cutoff *= direction ;
    1466   if (CoinAbs(cutoff) > 1.0e30)
    1467     assert (cutoff > 1.0e30) ;
    1468   double current = si.getObjValue() ;
    1469   current *= direction ;
    1470 /*
    14711750  Scan the variables, noting the ones that are fixed or have a bound too large
    14721751  to be useful.
    14731752*/
    1474   //int nFix = 0 ;
    14751753  for (int i = 0 ; i < nCols ; i++) {
    14761754    if (colUpper[i]-colLower[i] < 1.0e-8) {
    14771755      markC[i] = tightenLower|tightenUpper ;
    1478       //nFix++ ;
    14791756    } else {
    14801757      markC[i] = 0 ;
     
    14851762    }
    14861763  }
    1487   //printf("PROBE %d fixed out of %d\n",nFix,nCols) ;
    1488 
    14891764/*
    14901765  jjf: If we are going to replace coefficient then we don't need to be
     
    14931768  Seems like this comment does not agree with CglProbingRowCut.addCuts(),
    14941769  which checks effectiveness before adding a cut to the OsiCuts collection or
    1495   entering it into the strenghtenRow array.
     1770  entering it into the strengthenRow array.
    14961771
    14971772  But it does agree with the way coefficient strengthening cuts are handled
    14981773  down in the end-of-iteration processing for the down/up/down probe loop.
    14991774
    1500   It looks like the idea of strengthenRow is that a cut that's simply
    1501   strengthening an existing row i is entered in slot i of strengthenRow. It's
    1502   left to the client to deal with it on return. I need to get a better idea of
    1503   how justReplace is handled, here and in the client.  -- lh, 101125 --
     1775  The idea of strengthenRow is that a cut that's simply strengthening an
     1776  existing row i is entered in slot i of strengthenRow. It's left to the
     1777  client to deal with it on return. I need to get a better idea of how
     1778  justReplace is handled, here and in the client.  -- lh, 101125 --
     1779
    15041780*/
    15051781  //double needEffectiveness = info->strengthenRow ? -1.0e10 : 1.0e-3 ;
     
    15131789  int nstackC0 = -1 ;
    15141790  int nstackR,nstackC ;
    1515 /*
    1516   So, let's assume that the real business starts here, and the previous block
    1517   is irrelevant.
    1518 */
    15191791/*
    15201792  All the initialisation that occurs here is specific to CglTreeProbingInfo.
     
    15591831    saveFixingInfo = (info->initializeFixing(&si) > 0) ;
    15601832  }
     1833# if CGL_DEBUG > 0
     1834  std::cout
     1835    << "    prior to probe loop, justReplace "
     1836    << ((justReplace)?"true":"false") << ", required effectiveness "
     1837    << needEffectiveness << "." << std::endl ;
     1838# endif
    15611839/*
    15621840  PASS LOOP: HEAD
     
    15831861    nfixed = 0 ;
    15841862#   if CGL_DEBUG > 0
    1585     std::cout << "  probe: starting pass " << ipass << "." << std::endl ;
     1863    std::cout
     1864      << "  probe loop: starting pass " << ipass << "." << std::endl ;
    15861865#   endif
    15871866/*
     
    17492028#     if CGL_DEBUG > 1
    17502029      std::cout
    1751         << "  probe: looking at " << si.getColName(j) << " (" << j
     2030        << "  probe: looking at " << " x(" << j
    17522031        << ") = " << solval << ", l = " << colLower[j]
    17532032        << ", u = " << colUpper[j] << "." << std::endl ;
     
    19152194        need to sort out a couple of other places where goingToTrueBound
    19162195        controls behaviour.  -- lh, 101214 --
    1917 
    1918   TODO: Apropos the above, for example, the test for coefficient strengthening
    1919         cuts includes goingToTrueBound != 0  -- lh, 110105 --
    19202196*/
    19212197        if (goingToTrueBound && (colUpper[j]-colLower[j] > 1.5 || colLower[j]))
     
    19252201  of a general integer variable.
    19262202
    1927   TODO: Why is row strengthening limited to binary variables? If we can
    1928         figure out what to do, the limitation seems artificial.
    1929         -- lh, 101127 --
    1930         The test here also says that we can't have coefficient strengthening
     2203  TODO: The test here also says that we can't have coefficient strengthening
    19312204        without disaggregation. I don't see any reason for this dominance.
    19322205        -- lh, 110105 --
     
    19342207        if ((rowCuts&1) == 0)
    19352208          goingToTrueBound = 0 ;
    1936         bool canReplace = info->strengthenRow && (goingToTrueBound == 2) ;
    19372209
    19382210#       if CGL_DEBUG > 1
     
    19542226  objective change due to forced bound changes.
    19552227*/
    1956         double objVal = current ;
     2228        double objVal = si.getObjValue() ;
    19572229        if (solMovement*djs[j] > 0.0)
    19582230          objVal += solMovement*djs[j] ;
     
    23392611                        << "        " << nstackC-1
    23402612                        << " col " << kcol << " " << colsol[kcol]
    2341                         << " l " << saveL[nstackC]
     2613                        << " l " << saveL[nstackC-1]
    23422614                        << " -> " << newLower << std::endl ;
    23432615#                     endif
     
    23982670                        << "        " << nstackC-1
    23992671                        << " col " << kcol << " " << colsol[kcol]
    2400                         << " u " << saveU[nstackC]
     2672                        << " u " << saveU[nstackC-1]
    24012673                        << " -> " << newUpper << std::endl ;
    24022674#                     endif
     
    24992771                        << "        " << nstackC-1
    25002772                        << " col " << kcol << " " << colsol[kcol]
    2501                         << " l " << saveL[nstackC]
     2773                        << " l " << saveL[nstackC-1]
    25022774                        << " -> " << newLower << std::endl ;
    25032775#                     endif
     
    25922864                        << "        " << nstackC-1
    25932865                        << " col " << kcol << " " << colsol[kcol]
    2594                         << " u " << saveU[nstackC]
     2866                        << " u " << saveU[nstackC-1]
    25952867                        << " -> " << newUpper << std::endl ;
    25962868#                     endif
     
    27142986                        << "        " << nstackC-1
    27152987                        << " col " << kcol << " " << colsol[kcol]
    2716                         << " l " << saveL[nstackC]
     2988                        << " l " << saveL[nstackC-1]
    27172989                        << " -> " << newLower << std::endl ;
    27182990#                     endif
     
    27593031                        << "        " << nstackC-1
    27603032                        << " col " << kcol << " " << colsol[kcol]
    2761                         << " u " << saveU[nstackC]
     3033                        << " u " << saveU[nstackC-1]
    27623034                        << " -> " << newUpper << std::endl ;
    27633035#                     endif
     
    28493121                        << "        " << nstackC-1
    28503122                        << " col " << kcol << " " << colsol[kcol]
    2851                         << " u " << saveU[nstackC]
     3123                        << " u " << saveU[nstackC-1]
    28523124                        << " -> " << newUpper << std::endl ;
    28533125#                     endif
     
    29393211                        << "        " << nstackC-1
    29403212                        << " col " << kcol << " " << colsol[kcol]
    2941                         << " l " << saveL[nstackC]
     3213                        << " l " << saveL[nstackC-1]
    29423214                        << " -> " << newLower << std::endl ;
    29433215#                     endif
     
    30003272        -- lh, 101216 --
    30013273*/
    3002         if (notFeasible || (objVal > cutoff)) {
     3274        if (notFeasible || (useCutoff && (objVal > cutoff))) {
    30033275#         if CGL_DEBUG > 1
    3004           std::cout << "  Column " << j << " infeasible: " ;
    3005           if (notFeasible && (objVal > cutoff))
     3276          std::cout
     3277            << "    " << ((way[iway] == probeDown)?"down":"up")
     3278            << " probe for x<" << j << "> infeasible on " ;
     3279          if (notFeasible && (useCutoff && (objVal > cutoff)))
    30063280            std::cout << "primal bounds and objective cutoff" ;
    30073281          else if (notFeasible)
     
    31073381  strengthening.
    31083382 
    3109   Note that GoingToTrueBound was set to 0 way back in initial setup unless
     3383  Note that goingToTrueBound was set to 0 way back in initial setup unless
    31103384  the user requested disaggregation cuts (rowCuts&0x01).
    31113385*/
    31123386        if (iway == downIter) {
    31133387          if (!notFeasible) {
    3114             singletonRows(j,primalTolerance_,si,rowCopy,markC,
     3388            singletonRows(j,primalTolerance_,useObj,si,rowCopy,markC,
    31153389                          nstackC,stackC,saveL,saveU,
    31163390                          colUpper,colLower,columnGap,
     
    31223396            if ((rowCuts&0x02) != 0 && goingToTrueBound && !anyColumnCuts)
    31233397              strengthenCoeff(j,probeDown,primalTolerance_,justReplace,
    3124                               canReplace,needEffectiveness,si,rowCut,
     3398                              needEffectiveness,si,rowCut,
    31253399                              rowCopy,colUpper,colLower,colsol,nstackR,stackR,
    31263400                              rowUpper,rowLower,maxR,minR,realRows,
     
    31853459            iway = oneIterTooMany ;
    31863460
    3187             singletonRows(j,primalTolerance_,si,rowCopy,markC,
     3461            singletonRows(j,primalTolerance_,useObj,si,rowCopy,markC,
    31883462                          nstackC,stackC,saveL,saveU,
    31893463                          colUpper,colLower,columnGap,
     
    31963470            if ((rowCuts&0x02) != 0 && goingToTrueBound && !anyColumnCuts)
    31973471              strengthenCoeff(j,probeUp,primalTolerance_,justReplace,
    3198                               canReplace,needEffectiveness,si,rowCut,
     3472                              needEffectiveness,si,rowCut,
    31993473                              rowCopy,colUpper,colLower,colsol,nstackR,stackR,
    32003474                              rowUpper,rowLower,maxR,minR,realRows,
    32013475                              element,index,info) ;
    3202 /*
    3203   jjf: point back to stack
    3204 
    3205   We're done playing with singletons. Get to the real work here. We don't need
    3206   markC any more; coopt it for a cross-reference, var index -> stackC index.
    3207   The +1000 offset allows us to distinguish invalid entries (not all variables
    3208   are on stackC).
    3209 */
    3210             for (istackC = nstackC-1 ; istackC >= 0 ; istackC--) {
    3211               int icol = stackC[istackC] ;
    3212               markC[icol] = istackC+1000 ;
    3213             }
    3214             OsiColCut cc ;
    3215             int nTot = 0, nFix = 0, nInt = 0 ;
    3216             bool ifCut = false ;
    3217 /*
    3218   See if we have bounds improvement. Given a lower bound ld<j> from the
    3219   down probe, a bound lu<j> from the up probe, and the original bound lo<j>,
    3220   the bound we want is max(min(ld<j>,lu<j>),lo<j>). Use a conservative
    3221   tolerance.
    3222 */
    3223             for (istackC = 1 ; istackC < nstackC0 ; istackC++) {
    3224               int icol = stackC0[istackC] ;
    3225               int istackC1 = markC[icol]-1000 ;
    3226               if (istackC1 >= 0) {
    3227                 if (CoinMin(lo0[istackC],colLower[icol]) >
    3228                                               saveL[istackC1]+1.0e-4) {
    3229                   saveL[istackC1] = CoinMin(lo0[istackC],colLower[icol]) ;
    3230                   if (intVar[icol] /*||!info->inTree*/ ) {
    3231                     element[nFix] = saveL[istackC1] ;
    3232                     index[nFix++] = icol ;
    3233                     nInt++ ;
    3234                     if (colsol[icol] < saveL[istackC1]-primalTolerance_)
    3235                       ifCut = true ;
    3236                   }
    3237                   nfixed++ ;
    3238                 }
    3239               }
    3240             }
    3241             if (nFix) {
    3242               nTot = nFix ;
    3243               cc.setLbs(nFix,index,element) ;
    3244               nFix = 0 ;
    3245             }
    3246 /*
    3247   Repeat for upper bounds. There's an odd bit of asymmetry here; this loop
    3248   tries to generate a cut if there's no bound improvement.
    3249 
    3250   TODO: What, pray tell, is a `two cut'?  Binary variables only, it seems.
    3251         Judging from the conditions, we're looking for a variable that's
    3252         fixed at 0 on a down probe and fixed at 1 on an up probe, or vice
    3253         versa.
    3254         -- lh, 101128 --
    3255 */
    3256             for (istackC = 1 ; istackC < nstackC0 ; istackC++) {
    3257               int icol = stackC0[istackC] ;
    3258               int istackC1 = markC[icol]-1000 ;
    3259               if (istackC1 >= 0) {
    3260                 if (CoinMax(up0[istackC],colUpper[icol]) <
    3261                                               saveU[istackC1]-1.0e-4) {
    3262                   saveU[istackC1] = CoinMax(up0[istackC],colUpper[icol]) ;
    3263                   if (intVar[icol] /*||!info->inTree*/ ) {
    3264                     element[nFix] = saveU[istackC1] ;
    3265                     index[nFix++] = icol ;
    3266                     nInt++ ;
    3267                     if (colsol[icol] > saveU[istackC1]+primalTolerance_)
    3268                       ifCut = true ;
    3269                   }
    3270                   nfixed++ ;
    3271                 } else if (!info->inTree &&
    3272                            saveL[0] == 0.0 && saveU[0] == 1.0) {
    3273                   // See if can do two cut
    3274                   double upperWhenDown = up0[istackC] ;
    3275                   double lowerWhenDown = lo0[istackC] ;
    3276                   double upperWhenUp = colUpper[icol] ;
    3277                   double lowerWhenUp = colLower[icol] ;
    3278                   double upperOriginal = saveU[istackC1] ;
    3279                   double lowerOriginal = saveL[istackC1] ;
    3280                   if (upperWhenDown < lowerOriginal+1.0e-12 &&
    3281                       lowerWhenUp > upperOriginal-1.0e-12) {
    3282                     OsiRowCut rc ;
    3283                     rc.setLb(lowerOriginal) ;
    3284                     rc.setUb(lowerOriginal) ;
    3285                     rc.setEffectiveness(1.0e-5) ;
    3286                     int index[2] ;
    3287                     double element[2] ;
    3288                     index[0] = j ;
    3289                     index[1] = icol ;
    3290                     element[0] = -(upperOriginal-lowerOriginal) ;
    3291 /*
    3292   jjf: If zero then - must have been fixed without noticing!
    3293 
    3294   TODO: Uh, fixed without noticing?! That's an algorithm problem.
    3295         -- lh, 101128 --
    3296 */
    3297                     if (CoinAbs(element[0]) > 1.0e-8) {
    3298                       element[1] = 1.0 ;
    3299                       rc.setRow(2,index,element,false) ;
    3300                       cs.insert(rc) ;
    3301                     }
    3302                   } else if (upperWhenUp < lowerOriginal+1.0e-12 &&
    3303                              lowerWhenDown > upperOriginal-1.0e-12) {
    3304                     OsiRowCut rc ;
    3305                     rc.setLb(upperOriginal) ;
    3306                     rc.setUb(upperOriginal) ;
    3307                     rc.setEffectiveness(1.0e-5) ;
    3308                     int index[2] ;
    3309                     double element[2] ;
    3310                     index[0] = j ;
    3311                     index[1] = icol ;
    3312                     element[0] = upperOriginal-lowerOriginal ;
    3313                     element[1] = 1.0 ;
    3314                     rc.setRow(2,index,element,false) ;
    3315                     cs.insert(rc) ;
    3316                   }
    3317                 }
    3318               }
    3319             }    // end loop to update upper bounds (w/ two cuts)
    3320             if (nFix) {
    3321               nTot+=nFix ;
    3322               cc.setUbs(nFix,index,element) ;
    3323             }
    3324             // could tighten continuous as well
    3325             if (nInt) {
    3326               if (ifCut) {
    3327                 cc.setEffectiveness(100.0) ;
    3328               } else {
    3329                 cc.setEffectiveness(1.0e-5) ;
    3330               }
    3331 #if CGL_DEBUG > 0
    3332               CglProbingDebug::checkBounds(si,cc) ;
    3333 #endif
    3334               cs.insert(cc) ;
    3335             }
     3476            int jProbe = j ;
     3477            if (info->inTree || !(saveL[0] == 0.0 && saveU[0] == 1.0))
     3478              jProbe = -1 ;
     3479            implicationCuts(jProbe,primalTolerance_,nCols,cs,si,intVar,colsol,
     3480                            nstackC,stackC,markC,colUpper,colLower,
     3481                            saveU,saveL,nstackC0,stackC0,up0,lo0,
     3482                            element,index) ;
    33363483          }    // end of code to handle up and down feasible.
    33373484/*
     
    34253572  this method. If we're not in `just replace' mode, simply copy them over to
    34263573  the container (cs) passed as a parameter using addCuts. If we're in `just
    3427   replace' mode, ? need to explore further ?  -- lh, 101125 --
    3428 
    3429   TODO: If there's any possibility of realRows[i] < 0, this code will break
    3430         trying to read strengthenRow[]!
    3431         -- lh, 101128 --
     3574  replace' mode, they will be carried back to the caller in info.strengthenRow.
    34323575*/
    34333576  if (!ninfeas) {
     
    34353578      rowCut.addCuts(cs,info->strengthenRow,info->pass) ;
    34363579    } else {
     3580      assert(info->strengthenRow != 0) ;
    34373581      for (int i = 0 ; i < nRows ; i++) {
    3438         int realRow = realRows[i] ;
     3582        int realRow = (realRows)?realRows[i]:i ;
     3583        assert(realRow >= 0 && realRow < si.getNumRows()) ;
    34393584        OsiRowCut *cut = info->strengthenRow[realRow] ;
    34403585        if (realRow >= 0 && cut) {
    3441 #ifdef CLP_INVESTIGATE
    3442           printf("Row %d, real row %d effectiveness %g\n",
    3443                  i,realRow,cut->effectiveness()) ;
    3444 #endif
     3586#         if CGL_DEBUG > 1
     3587          std::cout
     3588            << "    row " << realRow << " (local " << i
     3589            << ") strengthened, effectiveness " << cut->effectiveness()
     3590            << "." << std::endl ;
     3591#         endif
    34453592          cs.insert(cut) ;
    34463593        }
     
    34483595    }
    34493596  }
    3450 #if 0
    3451 /*
    3452   Debug print.
    3453 */
    3454   {
    3455     int numberRowCutsAfter = cs.sizeRowCuts() ;
    3456     int k ;
    3457     for (k = 0;k<numberRowCutsAfter;k++) {
     3597# if CGL_DEBUG > 0
     3598  int numberRowCutsAfter = cs.sizeRowCuts() ;
     3599  int numberColCutsAfter = cs.sizeColCuts() ;
     3600
     3601  std::cout
     3602    << "Leaving CglProbing::probe, ninfeas " << ninfeas << ", "
     3603    << (numberRowCutsAfter-numberRowCutsBefore) << " row cuts, "
     3604    << (numberColCutsAfter-numberColCutsBefore) << " col cuts."
     3605    << std::endl ;
     3606 
     3607  if ((numberRowCutsAfter-numberRowCutsBefore) > 0) {
     3608    std::cout << "  row cuts:" << std::endl ;
     3609    for (int k = numberRowCutsBefore ; k < numberRowCutsAfter ; k++) {
    34583610      OsiRowCut thisCut = cs.rowCut(k) ;
    3459       printf("Cut %d is %g <=",k,thisCut.lb()) ;
    3460       int n=thisCut.row().getNumElements() ;
    3461       const int * column = thisCut.row().getIndices() ;
    3462       const double * element = thisCut.row().getElements() ;
    3463       assert (n) ;
    3464       for (int i=0;i<n;i++) {
    3465         printf(" %g*x%d",element[i],column[i]) ;
    3466       }
    3467       printf(" <= %g\n",thisCut.ub()) ;
     3611      thisCut.print() ;
    34683612    }
    34693613  }
    3470 #endif
    3471 
    3472 # if CGL_DEBUG > 0
    3473   std::cout
    3474     << "Leaving CglProbing::probe, ninfeas " << ninfeas
    3475     << "." << std::endl ;
     3614  if ((numberColCutsAfter-numberColCutsBefore) > 0) {
     3615    std::cout << "  col cuts:" << std::endl ;
     3616    for (int k = numberColCutsBefore ; k < numberColCutsAfter ; k++) {
     3617      OsiColCut thisCut = cs.colCut(k) ;
     3618      thisCut.print() ;
     3619    }
     3620  }
    34763621# endif
    34773622
  • branches/CglWorking101215/src/CglProbing/CglProbingRowCut.cpp

    r928 r964  
    1717#include "CglProbingRowCut.hpp"
    1818
    19 #define SIZE_ROW_MULT 4
    20 #define SIZE_ROW_ADD 2000
    21 
    22 /*
    23   Hash functions. Why is there no CoinUtils hash package? Apparently it was
    24   easier to just copy.  -- lh, 101202 --
     19#include "CglProbingDebug.hpp"
     20
     21/* \file CglProbingRowCut.cpp
     22   \brief Implementation of the CglProbingRowCut class.
     23*/
     24
     25/*
     26  Hash functions.
     27 
     28  TODO: Why is there no CoinUtils hash package?  -- lh, 101202 --
    2529*/
    2630namespace {
    2731
    28 double multiplier[] = { 1.23456789e2, -9.87654321 } ;
    29 
    30 int hashCut (const OsiRowCut2 & x, int size)
    31 {
    32   int xN =x.row().getNumElements();
    33   double xLb = x.lb();
    34   double xUb = x.ub();
    35   const int * xIndices = x.row().getIndices();
    36   const double * xElements = x.row().getElements();
    37   unsigned int hashValue;
    38   double value=1.0;
    39   if (xLb>-1.0e10)
    40     value += xLb*multiplier[0];
    41   if (xUb<1.0e10)
    42     value += xUb*multiplier[1];
    43   for( int j=0;j<xN;j++) {
    44     int xColumn = xIndices[j];
    45     double xValue = xElements[j];
    46     int k=(j&1);
    47     value += (j+1)*multiplier[k]*(xColumn+1)*xValue;
    48   }
    49   // should be compile time but too lazy for now
    50   if (sizeof(value)>sizeof(hashValue)) {
    51     assert (sizeof(value)==2*sizeof(hashValue));
    52     union { double d; int i[2]; } xx;
    53     xx.d = value;
    54     hashValue = (xx.i[0] + xx.i[1]);
     32const double multiplier[] = { 1.23456789e2, -9.87654321 } ;
     33
     34/*
     35  Calculate a hash value for the cut, modulo the size of the hash table.
     36  Note that the hash value is sensitive to coefficient order.
     37*/
     38int hashCut (const OsiRowCut &x, int size)
     39{
     40  int xN = x.row().getNumElements() ;
     41  double xLb = x.lb() ;
     42  double xUb = x.ub() ;
     43  const int *xIndices = x.row().getIndices() ;
     44  const double *xElements = x.row().getElements() ;
     45  double value = 1.0 ;
     46  if (xLb > -1.0e10)
     47    value += xLb*multiplier[0] ;
     48  if (xUb < 1.0e10)
     49    value += xUb*multiplier[1] ;
     50  // This calculation is order-sensitive
     51  for (int j = 0 ; j < xN ; j++) {
     52    int xColumn = xIndices[j] ;
     53    double xValue = xElements[j] ;
     54    int k = (j&1) ;
     55    value += (j+1)*multiplier[k]*(xColumn+1)*xValue ;
     56  }
     57/*
     58  It really would be nice if sizeof where legal in the preprocessor. But a
     59  good optimizing compiler shouldn't have any trouble.
     60*/
     61  unsigned int hashValue ;
     62  if (sizeof(double) == sizeof(int)) {
     63    union { double d; int i; } xx ;
     64    xx.d = value ;
     65    hashValue = xx.i ;
     66  } else if (sizeof(double) == 2*sizeof(int)) {
     67    union { double d; int i[2]; } xx ;
     68    xx.d = value ;
     69    hashValue = (xx.i[0] + xx.i[1]) ;
    5570  } else {
    56     assert (sizeof(value)==sizeof(hashValue));
    57     union { double d; unsigned int i[2]; } xx;
    58     xx.d = value;
    59     hashValue = xx.i[0];
    60   }
    61   return hashValue%(size);
    62 }
    63 
    64 bool same (const OsiRowCut2 & x, const OsiRowCut2 & y)
    65 {
    66   int xN =x.row().getNumElements();
    67   int yN =y.row().getNumElements();
    68   bool identical=false;
    69   if (xN==yN) {
    70     double xLb = x.lb();
    71     double xUb = x.ub();
    72     double yLb = y.lb();
    73     double yUb = y.ub();
    74     if (fabs(xLb-yLb)<1.0e-8&&fabs(xUb-yUb)<1.0e-8) {
    75       const int * xIndices = x.row().getIndices();
    76       const double * xElements = x.row().getElements();
    77       const int * yIndices = y.row().getIndices();
    78       const double * yElements = y.row().getElements();
    79       int j;
    80       for( j=0;j<xN;j++) {
    81         if (xIndices[j]!=yIndices[j])
    82           break;
    83         if (fabs(xElements[j]-yElements[j])>1.0e-12)
    84           break;
    85       }
    86       identical =  (j==xN);
    87     }
    88   }
    89   return identical;
     71    std::cout << "Hey! Time to implement another case, eh?" ;
     72    assert(false) ;
     73  }
     74
     75  return (hashValue%size) ;
     76
     77}
     78
     79/*
     80  Compare two cuts to see if they're the same. Invoked when there's a hash
     81  collision. Coefficient order is critical --- the coefficients are compared
     82  in the order given.
     83
     84  Arguably we could use OsiRowCut::operator==, but it's not a toleranced
     85  comparison, and that can be a problem. Also, we don't check efficiency here.
     86*/
     87bool same (const OsiRowCut &x, const OsiRowCut &y)
     88{
     89  int xN = x.row().getNumElements() ;
     90  int yN = y.row().getNumElements() ;
     91  bool identical = false ;
     92  if (xN == yN) {
     93    double xLb = x.lb() ;
     94    double xUb = x.ub() ;
     95    double yLb = y.lb() ;
     96    double yUb = y.ub() ;
     97    if (fabs(xLb-yLb) < 1.0e-8 && fabs(xUb-yUb) < 1.0e-8) {
     98      const int * xIndices = x.row().getIndices() ;
     99      const double * xElements = x.row().getElements() ;
     100      const int * yIndices = y.row().getIndices() ;
     101      const double * yElements = y.row().getElements() ;
     102      int j ;
     103      for (j = 0 ; j < xN ; j++) {
     104        if (xIndices[j] != yIndices[j]) break ;
     105        if (fabs(xElements[j]-yElements[j]) > 1.0e-12) break ;
     106      }
     107      identical = (j == xN) ;
     108    }
     109  }
     110  return (identical) ;
    90111}
    91112
    92113} // end file-local namespace
    93114
    94 
     115/*
     116  Note that while nRows is used to size the collection, the connection won't
     117  be obvious to the casual observer. Expanding the capacity later is
     118  expensive, so don't scrimp on the initial allocation.
     119*/
    95120CglProbingRowCut::CglProbingRowCut (int nRows, bool initialPass)
    96121{
    97   numberCuts_ = 0 ;
     122  // Calculate the maximum size of the collection
    98123  if (nRows < 500) {
    99124    maxSize_ = SIZE_ROW_MULT*nRows+SIZE_ROW_ADD ;
    100   } else if (nRows<5000) {
     125  } else if (nRows < 5000) {
    101126    maxSize_ = (SIZE_ROW_MULT*nRows+SIZE_ROW_ADD)>>1 ;
    102   } else if (nRows<10000) {
     127  } else if (nRows < 10000) {
    103128    maxSize_ = (SIZE_ROW_MULT*(nRows>>1)+SIZE_ROW_ADD)>>1 ;
    104129  } else {
    105130    maxSize_ = (SIZE_ROW_MULT*CoinMin(nRows,100000)+SIZE_ROW_ADD)>>2 ;
    106131  }
     132
     133  // Calculate the initial size of the cut collection
    107134  size_ = (maxSize_>>3)+10 ;
    108 /*
    109   TODO: When this is the initial pass, cut the capacity by half? Seems
    110         counterintuitive. For that matter, I don't expect the constructor
    111         to be making this decision for me.  -- lh, 100923 --
    112 */
    113135  if (initialPass)
    114136    size_ = size_>>1 ;
    115   if (size_<1000)
    116     hashSize_=4*size_ ;
     137  // Set up the vector of cuts
     138  nRows_ = nRows ;
     139  rowCut_ = new  OsiRowCut* [size_] ;
     140  numberCuts_ = 0 ;
     141
     142  // Set up the hash table
     143  if (size_ < 1000)
     144    hashSize_ = 4*size_ ;
    117145  else
    118     hashSize_=2*size_ ;
    119   nRows_ = nRows ;
    120   rowCut_ = new  OsiRowCut2 * [size_] ;
    121   hash_ = new HashLink[hashSize_] ;
    122   for (int i=0;i<hashSize_;i++) {
    123     hash_[i].index=-1 ;
    124     hash_[i].next=-1 ;
    125   }
    126   numberCuts_=0 ;
    127   lastHash_=-1 ;
    128 }
    129 
    130 
     146    hashSize_ = 2*size_ ;
     147  hash_ = new HashLink [hashSize_] ;
     148  for (int i = 0 ; i < hashSize_ ; i++) {
     149    hash_[i].index = -1 ;
     150    hash_[i].next = -1 ;
     151  }
     152  lastHash_ = -1 ;
     153}
     154
     155/*
     156  Destructor. Delete any cuts, then remove the cut vector and hash table.
     157*/
    131158CglProbingRowCut::~CglProbingRowCut ()
    132159{
    133   for (int i=0;i<numberCuts_;i++)
    134     delete rowCut_[i] ;
     160  for (int i = 0 ; i < numberCuts_ ; i++) delete rowCut_[i] ;
    135161  delete [] rowCut_ ;
    136162  delete [] hash_ ;
    137163}
    138164
    139 // Return 0 if added, 1 if not, -1 if not added because of space
    140 int CglProbingRowCut::addCutIfNotDuplicate(OsiRowCut & cut,int whichRow)
    141 {
    142   if (numberCuts_==size_&&numberCuts_<maxSize_) {
     165/*
     166  If the cut is not in the table, and there is room for the cut, returns
     167  the index of the appropriate empty hash table entry (index field = -1).
     168  If the cut is already in the table, returns the entry for the cut (index
     169  field >= 0).  If there's no room, returns -1
     170*/
     171int CglProbingRowCut::findHashTableEntry (const OsiRowCut &newCut)
     172{
     173/*
     174  Hash the cut and hope we hit an empty entry. If so, we're done.
     175*/
     176  int newNdx = hashCut(newCut,hashSize_) ;
     177  if (hash_[newNdx].index == -1) return (newNdx) ;
     178/*
     179  Collision. Walk the chain of entries starting from here, looking for a
     180  duplicate cut or the end of the chain.
     181*/
     182  int dupNdx = -1 ;
     183  for (int ndx = newNdx ; ndx >= 0 ; ndx = hash_[ndx].next) {
     184    int cutNdx = hash_[ndx].index ;
     185    if (same(newCut,*rowCut_[cutNdx])) {
     186      dupNdx = ndx ;
     187      break ;
     188    }
     189  }
     190  if (dupNdx >= 0) return (dupNdx) ;
     191/*
     192  We didn't find a duplicate, so we need to add to the chain. Start a linear
     193  search for an empty entry from the point where we ended the last search.
     194*/
     195  assert(lastHash_ < hashSize_) ;
     196  int startSearch = lastHash_ ;
     197  for (lastHash_++ ;
     198       lastHash_ < hashSize_ && hash_[lastHash_].index >= 0 ;
     199       lastHash_++) ;
     200  if (lastHash_ >= hashSize_) {
     201    for (lastHash_ = 0 ;
     202         lastHash_ < startSearch && hash_[lastHash_].index >= 0 ;
     203         lastHash_++) ;
     204  }
     205  if (hash_[lastHash_].index >= 0) {
     206    std::cout
     207      << "  CglProbingRowCut::insertInHashTable: no room in hash table!"
     208      << std::endl ;
     209    return (-1) ;
     210  }
     211  assert(hash_[lastHash_].index == -1) ;
     212/*
     213  Insert the new link at the beginning of the chain
     214*/
     215  hash_[lastHash_].next = hash_[newNdx].next ;
     216  hash_[newNdx].next = lastHash_ ;
     217
     218  return (lastHash_) ;
     219}
     220
     221
     222/*
     223  Add a cut only if it's not a duplicate of an existing cut.
     224 
     225  Return 0 if added, 1 if not, -1 if not added because of space
     226*/
     227int CglProbingRowCut::addCutIfNotDuplicate(OsiRowCut &cut)
     228{
     229/*
     230  Time to expand the cut storage vector and hash table? We have to repopulate
     231  the hash table, so this is a fairly expensive operation. If we're at the
     232  max, return -1.
     233*/
     234  if (numberCuts_ >= size_) {
     235    if (numberCuts_ >= maxSize_) return (-1) ;
    143236    size_ = CoinMin(2*size_+100,maxSize_) ;
    144     if (size_<1000)
    145       hashSize_=4*size_ ;
     237    if (size_ < 1000)
     238      hashSize_ = 4*size_ ;
    146239    else
    147       hashSize_=2*size_ ;
    148 #ifdef COIN_DEVELOP
    149     printf("increaing size from %d to %d (hash size %d, maxsize %d)\n",
    150            numberCuts_,size_,hashSize_,maxSize_) ;
    151 #endif
    152     OsiRowCut2 ** temp = new  OsiRowCut2 * [size_] ;
    153     delete [] hash_ ;
     240      hashSize_ = 2*size_ ;
     241#   if CGL_DEBUG > 0
     242    std::cout
     243      << "  addCutIfNotDuplicate: increasing size from " << numberCuts_
     244      << " to " << size_ << " (max " << maxSize_ << ", hash " << hashSize_
     245      << "." << std::endl ;
     246#   endif
     247    OsiRowCut **oldRowCut = rowCut_ ;
     248    rowCut_ = new  OsiRowCut* [size_] ;
     249    HashLink *oldHash = hash_ ;
    154250    hash_ = new HashLink[hashSize_] ;
    155     for (int i=0;i<hashSize_;i++) {
    156       hash_[i].index=-1 ;
    157       hash_[i].next=-1 ;
    158     }
    159     for (int i=0;i<numberCuts_;i++) {
    160       temp[i]=rowCut_[i] ;
    161       int ipos = hashCut(*temp[i],hashSize_) ;
    162       int found = -1 ;
    163       int jpos=ipos ;
    164       while ( true ) {
    165         int j1 = hash_[ipos].index ;
    166        
    167         if ( j1 >= 0 ) {
    168           if ( !same(*temp[i],*temp[j1]) ) {
    169             int k = hash_[ipos].next ;
    170             if ( k != -1 )
    171               ipos = k ;
    172             else
    173               break ;
    174           } else {
    175             found = j1 ;
    176             break ;
    177           }
    178         } else {
    179           break ;
    180         }
    181       }
    182       if (found<0) {
    183         assert (hash_[ipos].next==-1) ;
    184         if (ipos==jpos) {
    185           // first
    186           hash_[ipos].index=i ;
    187         } else {
    188           // find next space
    189           while ( true ) {
    190             ++lastHash_ ;
    191             assert (lastHash_<hashSize_) ;
    192             if ( hash_[lastHash_].index == -1 )
    193               break ;
    194           }
    195           hash_[ipos].next = lastHash_ ;
    196           hash_[lastHash_].index = i ;
    197         }
    198       }
    199     }
    200     delete [] rowCut_ ;
    201     rowCut_ = temp ;
    202   }
    203   if (numberCuts_<size_) {
    204     double newLb = cut.lb() ;
    205     double newUb = cut.ub() ;
    206     CoinPackedVector vector = cut.row() ;
    207     int numberElements =vector.getNumElements() ;
    208     int * newIndices = vector.getIndices() ;
    209     double * newElements = vector.getElements() ;
    210     CoinSort_2(newIndices,newIndices+numberElements,newElements) ;
    211     int i ;
    212     bool bad=false ;
    213     for (i=0;i<numberElements;i++) {
    214       double value = fabs(newElements[i]) ;
    215       if (value<1.0e-12||value>1.0e12)
    216         bad=true ;
    217     }
    218     if (bad)
    219       return 1 ;
    220     OsiRowCut2 newCut(whichRow) ;
    221     newCut.setLb(newLb) ;
    222     newCut.setUb(newUb) ;
    223     newCut.setRow(vector) ;
    224     int ipos = hashCut(newCut,hashSize_) ;
    225     int found = -1 ;
    226     int jpos=ipos ;
    227     while ( true ) {
    228       int j1 = hash_[ipos].index ;
    229 
    230       if ( j1 >= 0 ) {
    231         if ( !same(newCut,*rowCut_[j1]) ) {
    232           int k = hash_[ipos].next ;
    233           if ( k != -1 )
    234             ipos = k ;
    235           else
    236             break ;
    237         } else {
    238           found = j1 ;
    239           break ;
    240         }
    241       } else {
     251    for (int i = 0 ; i < hashSize_ ; i++) {
     252      hash_[i].index = -1 ;
     253      hash_[i].next = -1 ;
     254    }
     255/*
     256  Repopulate the hash table. A return of -1 indicates no space in the hash
     257  table, which really shouldn't happen here.
     258*/
     259    bool failed = false ;
     260    for (int i = 0 ; i < numberCuts_ ; i++) {
     261      rowCut_[i] = oldRowCut[i] ;
     262      int hashNdx = findHashTableEntry(*rowCut_[i]) ;
     263      if (hashNdx < 0) {
     264        failed = true ;
    242265        break ;
    243266      }
    244     }
    245     if (found<0) {
    246       assert (hash_[ipos].next==-1) ;
    247       if (ipos==jpos) {
    248         // first
    249         hash_[ipos].index=numberCuts_ ;
    250       } else {
    251         // find next space
    252         while ( true ) {
    253           ++lastHash_ ;
    254           assert (lastHash_<hashSize_) ;
    255           if ( hash_[lastHash_].index == -1 )
    256             break ;
    257         }
    258         hash_[ipos].next = lastHash_ ;
    259         hash_[lastHash_].index = numberCuts_ ;
    260       }
    261       OsiRowCut2 * newCutPtr = new OsiRowCut2(whichRow) ;
    262       newCutPtr->setLb(newLb) ;
    263       newCutPtr->setUb(newUb) ;
    264       newCutPtr->setRow(vector) ;
    265       rowCut_[numberCuts_++]=newCutPtr ;
    266       return 0 ;
    267     } else {
    268       return 1 ;
    269     }
    270   } else {
    271     return -1 ;
    272   }
    273 }
    274 
     267      if (hash_[hashNdx].index < 0) {
     268        hash_[hashNdx].index = i ;
     269      }
     270#     if CGL_DEBUG > 0
     271      else {
     272        std::cout
     273          << "    addCutIfNotDuplicate: encountered duplicate cut while "
     274          << "repopulating hash table after expansion." << std::endl ;
     275      }
     276#     endif
     277    }
     278    if (failed) {
     279      delete [] rowCut_ ;
     280      rowCut_ = oldRowCut ;
     281      delete [] hash_ ;
     282      hash_ = oldHash ;
     283      return (-1) ;
     284    }
     285    delete [] oldHash ;
     286    delete [] oldRowCut ;
     287  }
     288/*
     289  There is now space to add the cut, or we would have already returned.
     290*/
     291  CoinPackedVector &cutVec = cut.mutableRow() ;
     292# if 1 // CGL_DEBUG > 0
     293  int bad = 0 ;
     294  int coeffCnt = cutVec.getNumElements() ;
     295  const double *constCoeffs = cutVec.getElements() ;
     296  for (int i = 0 ; i < coeffCnt ; i++) {
     297    double coeff = fabs(constCoeffs[i]) ;
     298    if (CoinAbs(coeff) < 1.0e-12 || CoinAbs(coeff) > 1.0e12) bad++ ;
     299  }
     300  if (bad) {
     301    std::cout
     302      << "    addCutIfNotDuplicate: BAD: " << bad << " coefficients in cut."
     303      << std::endl ;
     304    assert(false) ;
     305  }
     306# endif
     307/*
     308  Find the hash table entry that matches the cut. We need to sort the indices
     309  to ensure a standard order, because the hash function is order-sensitive.
     310  The exact order isn't important.
     311*/
     312  cutVec.sortIncrIndex() ;
     313  int hashNdx = findHashTableEntry(cut) ;
     314  if (hash_[hashNdx].index < 0)
     315  { OsiRowCut *newCut = cut.clone() ;
     316    rowCut_[numberCuts_] = newCut ;
     317    hash_[hashNdx].index = numberCuts_ ;
     318    numberCuts_++ ;
     319  }
     320
     321  return (0) ;
     322}
     323
     324/*
     325  Transfer up to nRows_ cuts out of this CglProbingRowCut structure
     326  into the standard OsiCuts structure given as a parameter. If the total
     327  number of cuts is less than nRows_, take them all, otherwise, just take
     328  the best. nRows_ is a nominal limit which can be exceeded if the cuts
     329  straddling the limit all have the same effectiveness.
     330
     331  If whichRows is supplied, entries which do not already associate a cut with
     332  the row will be filled in, if a cut is transferred that identifies itself
     333  with the row.
     334
     335  When we're taking all cuts, scan from first or last cut on alternate passes,
     336  so that cuts generated later in the pass have a chance of being associated
     337  with a row. When we sort for the best cuts, use the best.
     338*/
    275339void CglProbingRowCut::addCuts (OsiCuts &cs, OsiRowCut **whichRow, int iPass)
    276340{
    277   int numberCuts=cs.sizeRowCuts() ;
     341  int numberCuts = cs.sizeRowCuts() ;
    278342  int i ;
    279   if (numberCuts_<nRows_) {
    280     if ((iPass&1)==1) {
    281       for (i=0;i<numberCuts_;i++) {
     343/*
     344  Take all, working from first or last on alternate passes so that later cuts
     345  have a chance to claim places in whichRow.
     346*/
     347  if (numberCuts_ < nRows_) {
     348    if ((iPass&1) == 1) {
     349      for (i = 0 ; i < numberCuts_ ; i++) {
    282350        cs.insert(*rowCut_[i]) ;
    283351        if (whichRow) {
    284           int iRow= rowCut_[i]->whichRow() ;
    285           if (iRow>=0&&!whichRow[iRow])
    286             whichRow[iRow]=cs.rowCutPtr(numberCuts); ;
     352          int iRow = rowCut_[i]->whichRow() ;
     353          if (iRow >= 0 && !whichRow[iRow])
     354            whichRow[iRow] = cs.rowCutPtr(numberCuts) ;
    287355        }
    288356        numberCuts++ ;
    289357      }
    290358    } else {
    291       for (i=numberCuts_-1;i>=0;i--) {
     359      for (i = numberCuts_-1 ; i >= 0 ; i--) {
    292360        cs.insert(*rowCut_[i]) ;
    293361        if (whichRow) {
    294           int iRow= rowCut_[i]->whichRow() ;
    295           if (iRow>=0&&!whichRow[iRow])
    296             whichRow[iRow]=cs.rowCutPtr(numberCuts); ;
     362          int iRow = rowCut_[i]->whichRow() ;
     363          if (iRow >= 0 && !whichRow[iRow])
     364            whichRow[iRow] = cs.rowCutPtr(numberCuts) ;
    297365        }
    298366        numberCuts++ ;
     
    300368    }
    301369  } else {
    302     // just best
    303     double * effectiveness = new double[numberCuts_] ;
    304     int iCut=0 ;
    305     for (i=0;i<numberCuts_;i++) {
    306       double value = -rowCut_[i]->effectiveness() ;
     370/*
     371  Transfer just the best cuts. Load effectiveness values into a vector, sort,
     372  select the threshold, then scan the cuts and take the ones over the
     373  threshold. Cuts that are associated with rows get star billing.
     374*/
     375    double *effectiveness = new double[numberCuts_] ;
     376    int iCut = 0 ;
     377    for (iCut = 0 ; iCut < numberCuts_ ; iCut++) {
     378      double value = -rowCut_[iCut]->effectiveness() ;
    307379      if (whichRow) {
    308         int iRow= rowCut_[i]->whichRow() ;
    309         if (iRow>=0)
     380        int iRow = rowCut_[iCut]->whichRow() ;
     381        if (iRow >= 0)
    310382          value -= 1.0e10 ;
    311383      }
    312       effectiveness[iCut++]=value ;
     384      effectiveness[iCut] = value ;
    313385    }
    314386    std::sort(effectiveness,effectiveness+numberCuts_) ;
    315     double threshold = -1.0e20 ;
    316     if (iCut>nRows_)
    317       threshold = effectiveness[nRows_] ;
    318     for ( i=0;i<numberCuts_;i++) {
    319       if (rowCut_[i]->effectiveness()>threshold) {
     387    double threshold = effectiveness[nRows_-1] ;
     388    for (i = 0 ; i < numberCuts_ ; i++) {
     389      if (rowCut_[i]->effectiveness() > threshold) {
    320390        cs.insert(*rowCut_[i]) ;
    321391        if (whichRow) {
    322           int iRow= rowCut_[i]->whichRow() ;
    323           if (iRow>=0&&!whichRow[iRow])
    324             whichRow[iRow]=cs.rowCutPtr(numberCuts); ;
     392          int iRow = rowCut_[i]->whichRow() ;
     393          if (iRow >= 0 && !whichRow[iRow])
     394            whichRow[iRow] = cs.rowCutPtr(numberCuts) ;
    325395        }
    326396        numberCuts++ ;
     
    332402  { delete rowCut_[i] ;
    333403    rowCut_[i] = 0 ; }
    334   numberCuts_=0 ;
    335 }
    336 
     404  numberCuts_ = 0 ;
     405}
     406
  • branches/CglWorking101215/src/CglProbing/CglProbingRowCut.hpp

    r928 r964  
    1111#define CglProbingRowCut_H
    1212
    13 /*
    14   TODO: What does this provide that's better / distinct from OsiCuts?
    15         There's a cryptic comment down in gutsOfGenerateCuts, `keep cuts out
    16         of cs until end so we can find duplicates quickly', but that doesn't
    17         quite justify a new class.  -- lh, 100923 --
    18 
    19   TODO: This should be a separate file. It's about 200 lines.
    20         -- lh, 100923 --
    21 
    22   TODO: This should be renamed so that it follows the usual conventions for a
    23         class name.  -- lh, 100923 --
    24 
    25   Reading a bit more, it would appear that this class has two reasons for
    26   existence: OsiCuts is more general, handling row and column cuts, and it
    27   uses OsiRowCut, instead of OsiRowCut2 used here. OsiRowCut2 adds a
    28   reference (simple integer index) to the original constraint. Remains to be
    29   seen whether this really justifies a new class or whether it'd be better to
    30   add the reference to the original constraint in OsiRowCut.  -- lh, 101125 --
    31 
    32   Separated and renamed CglProbingRowCut  -- lh, 101202 --
    33 */
    34 
    3513#include <OsiCuts.hpp>
    3614
     15/*
     16  \brief Augmented row cut class to accumulate probing cuts.
     17
     18  A collection of OsiRowCut objects, specialised for use in CglProbing.
     19  It supports fast identification of duplicate cuts and some bookkeeping
     20  that's useful in CglProbing.
     21
     22  The value given to the constructor for \p nRows will be used to to set
     23  the capacity of the collection, based on a complicated calculation that
     24  tries to come up with something that's unlikely to overflow during the
     25  probe. Expansion is costly because the hash table must be expanded and
     26  repopulated.  The actual capacity will have no obvious relation to \p
     27  nRows. Check the code for the constructor if this concerns you.
     28*/
    3729class CglProbingRowCut {
    3830
    3931  public:
    4032
     33  /*!  \brief Default constructor
     34
     35    The value given for \p nRows will be used to set #nRows_. It's also
     36    used in a magic calculation to set the capacity of the collection. See
     37    the notes for the class for details.  Setting \p initialPass to true
     38    will reduce the size of the hash table by half, in the context of a
     39    similar magic calculation.
     40  */
    4141  CglProbingRowCut(int nRows, bool initialPass) ;
    4242
     43  /// Destructor
    4344  ~CglProbingRowCut() ;
    4445
    45   inline OsiRowCut2 *cut (int i) const { return rowCut_[i] ; }
     46  inline OsiRowCut *cut (int i) const { return rowCut_[i] ; }
    4647
     48  /// The number of cuts in the collection.
    4749  inline int numberCuts() const { return numberCuts_ ; }
    4850
     51  /// Check if space remains to add more cuts.
    4952  inline bool outOfSpace() const { return (maxSize_ == numberCuts_) ; }
    5053
    51   // Return 0 if added, 1 if not, -1 if not added because of space
    52   int addCutIfNotDuplicate(OsiRowCut &cut, int whichRow = -1) ;
     54  /*! \brief Add a cut if it's not a duplicate
     55 
     56    Add a cut if it's not a duplicate of a cut already in the collection.
    5357
     58    Return value:
     59     - -1 if the cut is not added because there's no space,
     60     -  0 if the cut is added,
     61     -  1 if the cut is not added because it's a duplicate
     62  */
     63  int addCutIfNotDuplicate(OsiRowCut &cut) ;
     64
     65  /*! \brief Transfer cuts into the given OsiCuts collection, updating
     66             whichRow.
     67
     68    Transfer cuts from this collection into \p cs, deleting all cuts from this
     69    collection when finished. Nominally, up to #nRows_ cuts are transferred
     70    out but this limit can be exceeded if the cuts straddling the limit have
     71    equal effectiveness.
     72
     73    During the transfer, if a cut has an associated row index, \p whichRow is
     74    checked. If the entry for the row has no associated cut, a pointer to the
     75    cut is stored. The \p iPass parameter gives some randomness in lieu of
     76    sorting when the number of cuts in the container is less than the transfer
     77    limit.
     78  */
    5479  void addCuts(OsiCuts &cs, OsiRowCut **whichRow, int iPass) ;
    5580
    5681  private:
    5782
    58   typedef struct { int index, next ; } HashLink;
     83  /*! \brief find an entry for a cut in the hash table
    5984
    60   OsiRowCut2 **rowCut_ ;
    61   /// Hash table
     85    If the cut is already in the hash table, returns the index of the entry.
     86    If the cut is not in the hash table, returns the index of an appropriate
     87    empty entry. If there's no space, returns -1.
     88  */
     89    int findHashTableEntry (const OsiRowCut &newCut) ;
     90
     91  /*! \brief hash table entry
     92 
     93    index is the cut index, next is the next hash table entry for this
     94    hash value
     95  */
     96
     97  typedef struct { int index ;
     98                   int next ; } HashLink ;
     99
     100  /// Fudge factor used by constructor to calculate #maxSize_
     101  static const int SIZE_ROW_MULT = 4 ;
     102  /// Another fudge factor used by constructor to calculate #maxSize_
     103  static const int SIZE_ROW_ADD = 2000 ;
     104
     105  /// The number of cuts in the collection
     106  int numberCuts_ ;
     107  /// The maximum number of cuts that can be held in the collection.
     108  int maxSize_ ;
     109  /// The initial size of the cut storage vector
     110  int size_ ;
     111  /// The cut storage vector
     112  OsiRowCut **rowCut_ ;
     113  /// Nominal number of cuts transferred out by #addCuts
     114  int nRows_ ;
     115
     116  /// Size of hash table
     117  int hashSize_ ;
     118  /// Hash table for fast duplicate check
    62119  HashLink *hash_ ;
    63   int size_ ;
    64   int maxSize_ ;
    65   int hashSize_ ;
    66   int nRows_ ;
    67   int numberCuts_ ;
     120  /*! \brief Most recently filled hash table entry
     121 
     122    And the starting point to search for a new empty entry when there's a hash
     123    collision.
     124  */
    68125  int lastHash_ ;
    69126
  • branches/CglWorking101215/src/CglProbing/CglProbingTest.cpp

    r957 r964  
    130130    siP->setIntParam(OsiNameDiscipline,1) ;
    131131
    132     std::string localDir = "/export/Lou2/Split/Data/Miplib3/" ;
     132    std::string localDir = "/devel/Coin/Split/Data/Miplib3/" ;
     133    // std::string localDir = "/export/Lou2/Split/Data/Miplib3/" ;
    133134    // std::string probName = "p0033" ;
    134135    // std::string fn = mpsDir+probName ;
     
    210211    siP->setIntParam(OsiNameDiscipline,1) ;
    211212
    212     std::string localDir = "/export/Lou2/Split/Data/Sample/" ;
     213    std::string localDir = "/devel/Coin/Split/Data/Sample/" ;
    213214    std::string probName = "p0033" ;
    214215    std::string fn = localDir+probName ;
  • branches/CglWorking101215/src/CglProbing/CglProbingTighten.cpp

    r963 r964  
    255255
    256256  TODO: There doesn't seem to be any good reason to pass the guts of the row
    257         and column copies as parameters. It should be easy enough to unpack
    258         them.  -- lh, 100924 --
     257        copy as parameters. It should be easy enough to unpack them.
     258        -- lh, 100924 --
    259259
    260260  TODO: minR and maxR will have artificial values of +/- 1e60 installed as
     
    274274*/
    275275
    276 void
    277 CglProbing::calcRowBounds(double *colLower, double * colUpper,
    278                              const int *column, const double *rowElements,
    279                              const CoinBigIndex *rowStart,
    280                              const int * rowLength,
    281                              double *rowLower, double *rowUpper,
    282                              double * minR, double * maxR, int * markR,
    283                              int nRows) const
     276void CglProbing::calcRowBounds (double *colLower, double * colUpper,
     277                                const int *column, const double *rowElements,
     278                                const CoinBigIndex *rowStart,
     279                                const int * rowLength,
     280                                double *rowLower, double *rowUpper,
     281                                double * minR, double * maxR, int * markR,
     282                                int nRows) const
    284283{
    285284  int i, j, k, kre ;
     
    370369#     endif
    371370    }
    372 #   if CGL_DEBUG > 1
    373     dump_row(i,rowLower[i],rowUpper[i],minR[i],maxR[i],0,false,false,0,
    374              rowLength[i],&column[krs],&rowElements[krs],
    375              primalTolerance_,colLower,colUpper) ;
     371#   if CGL_DEBUG > 2
     372    // Just in case we executed the else and didn't set these ...
     373    krs = rowStart[i] ;
     374    kre = rowStart[i]+rowLength[i] ;
     375    CglProbingDebug::dump_row(i,rowLower[i],rowUpper[i],minR[i],maxR[i],
     376                              0,false,false,0,
     377                              rowLength[i],&column[krs],&rowElements[krs],
     378                              primalTolerance_,colLower,colUpper) ;
    376379#   endif
    377380  }
Note: See TracChangeset for help on using the changeset viewer.