source: branches/CglWorking101215/src/CglProbing/CglProbingProbe.cpp @ 961

Last change on this file since 961 was 961, checked in by lou, 11 years ago

(nonworking) Rewrite strengthenCoeff. Different case analysis to make earlier
decisions about goodness of cut.

File size: 114.7 KB
Line 
1
2/*
3  Copyright (C) 2010, International Business Machines Corporation and others.
4  All Rights Reserved.
5
6  This code is licensed under the terms of the Eclipse Public License (EPL).
7
8  $Id$
9*/
10
11#if defined(_MSC_VER)
12// Turn off compiler warning about long names
13#  pragma warning(disable:4786)
14#endif
15
16#define PROBING100 0
17
18// #define PRINT_DEBUG
19
20#include "CoinPackedMatrix.hpp"
21#include "CglProbing.hpp"
22#include "CglProbingRowCut.hpp"
23
24// To enable debugging, set CGL_DEBUG in CglProbingDebug.hpp
25#include "CglProbingDebug.hpp"
26
27namespace {
28
29/*
30  Allocating one big block for data arrays is efficient but makes it
31  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.
34*/
35#if CGL_DEBUG == 0
36
37const unsigned int DIratio = sizeof(double)/sizeof(int) ;
38
39#define ONE_ARRAY
40
41#endif
42
43/*
44  Bit constants used to specify bound information.
45
46   0x0: no bounds tightened
47   0x1: u<j> tightened
48   0x2: l<j> tightened
49   0x4: l<j> < -1e10 (-infty)
50   0x8: u<j> >  1e10 (+infty)
51*/
52const unsigned int tightenUpper = 0x1 ;
53const unsigned int tightenLower = 0x2 ;
54const unsigned int infLower = 0x4 ;
55const unsigned int infUpper = 0x8 ;
56
57/*
58  Integer and bit constants used to specify probing.
59
60  The second down pass is used to restore the down probe state when down
61  is feasible but up is not.
62
63  The first set (downProbe, upProbe, downProbe2, oneProbeTooMany) *must*
64  match the iteration number for the corresponding pass in the PROBE
65  loop. The symbolic names make the code a bit easier to read.
66*/
67const unsigned int downProbe = 0 ;
68const unsigned int upProbe = 1 ;
69const unsigned int downProbe2 = 2 ;
70const unsigned int oneProbeTooMany = 3 ;
71
72const unsigned int probeDown = tightenUpper ;
73const unsigned int probeUp = tightenLower ;
74
75const unsigned int downProbeFeas = 0x1 ;
76const unsigned int upProbeFeas = 0x2 ;
77const unsigned int downProbe2Feas = 0x4 ;
78
79}  // end file-local namespace
80
81// ===============================================
82
83
84/*
85  Run through stackC and create disaggregation cuts. See the typeset
86  documentation for the full derivation. Suppose the probe variable is x<p>
87  and we're trying to generate a disaggregation cut for target x<t>.
88
89  For a down probe, we have orig_u<p> reduced to u<p> = down = colUpper[p].
90  It may be that this causes orig_u<t> to be reduced to u<t>. The upper
91  disaggregation cut is
92    -(orig_u<t>-u<t>)x<p> + x<t> <= u<t> - u<p>(orig_u<t>-u<t>)
93
94  It may be that this causes orig_l<t> to be increased to l<t>. The lower
95  disaggregation cut is
96    -(orig_l<t>-l<t>)x<p> + x<t> >= l<t> - u<p>(orig_l<t>-l<t>)
97
98  For an up probe, we have orig_l<p> increased to l<p> = up = colLower[p].
99  It may be that this causes orig_u<t> to be reduced to u<t>. The upper
100  disaggregation cut is
101     (orig_u<t>-u<t>)x<p> + x<t> <= u<t> + l<p>(orig_u<t>-u<t>)
102
103  It may be that this causes orig_l<t> to be increased to l<t>. The lower
104  disaggregation cut is
105     (orig_l<t>-l<t>)x<p> + x<t> >= l<t> + l<p>(orig_l<t>-l<t>)
106
107  Note that stackC[0] = pndx, the index of the probed variable, so we
108  don't want to process stackC[0] in the main loop.
109
110  It must be true that goingToTrueBound == 2 and justReplace == false when
111  this method is called.
112
113  TODO: If you look at the math, the critical thing to get the cut form
114        used here is that the old and new bounds on x<p> be separated by
115        one. E.g., for a down probe, (orig_u<p>-u<p>) = 1.  It looks to
116        me like goingToTrueBound might have been an attempt to guarantee this
117        condition, but the test is wrong, hence it's patched up with tweak
118        that essentially cuts this cut back to binary probing variables.
119        -- lh, 101214 --
120
121  probeDir should be coded using probeDown, probeUp.
122*/
123void disaggCuts (int nstackC, unsigned int probeDir,
124                 double primalTolerance_, double disaggEffectiveness,
125                 const OsiSolverInterface &si,
126                 CglProbingRowCut &rowCut, const int *const stackC,
127                 const double *const colsol,
128                 const double *const colUpper, const double *const colLower,
129                 const double *const saveU, const double *const saveL,
130                 int *const index, double *const element)
131{ 
132  int pndx = stackC[0] ;
133
134# if CGL_DEBUG > 0
135  std::cout
136    << "Entering disaggCuts, probed variable "
137    << si.getColName(pndx) << "(" << pndx << ")." << std::endl ;
138  const OsiRowCutDebugger *debugger = si.getRowCutDebugger() ;
139  // down or up are the only possibilities
140  assert((probeDir&(~(probeDown|probeUp))) == 0) ;
141# endif
142
143  int plusMinus = 0 ;
144  double lu_p = 0.0 ;
145  double origsol_p = colsol[pndx] ;
146  double origdelta_p = 0.0 ;
147/*
148  Set up the quantities that vary depending on whether this is an up or down
149  probe:
150  * lu_p is the new upper (down probe) or lower (up probe) bound on x<p>
151  * plusMinus will be -1 (down probe) or +1 (up probe).
152  * origdelta_p will be x*<p>-u<p> (down probe) or l<p>-x*<p> (up probe)
153*/
154  if (probeDir == probeDown) {
155    lu_p = colUpper[pndx] ;
156    plusMinus = -1 ;
157    origdelta_p = origsol_p-lu_p ;
158  } else {
159    lu_p = colLower[pndx] ;
160    plusMinus = 1 ;
161    origdelta_p = lu_p-origsol_p ;
162  }
163
164  for (int istackC = nstackC-1 ; istackC > 0 ; istackC--) {
165    int icol = stackC[istackC] ;
166    double origu_t = saveU[istackC] ;
167    double origl_t = saveL[istackC] ;
168    double u_t = colUpper[icol] ;
169    double sol_t = colsol[icol] ;
170    double l_t = colLower[icol] ;
171    double boundChange = origu_t-u_t ;
172    double coeff = plusMinus*boundChange ;
173/*
174  Generate the upper disaggregation cut, if it's violated and the coefficients
175  will be reasonable. The effectiveness is the improvement (change) in
176  x<pndx>.
177*/
178    if (boundChange > 0.0 && origu_t < 1.0e10 &&
179        (sol_t > u_t+boundChange*origdelta_p+primalTolerance_)) {
180      OsiRowCut rc ;
181      rc.setLb(-DBL_MAX) ;
182      rc.setUb(u_t+coeff*lu_p) ;
183      index[0] = icol ;
184      element[0] = 1.0 ;
185      index[1] = pndx ;
186      element[1] = coeff ;
187/*
188  Effectiveness is how far x<p> moves. We can save effort by transforming
189  sol_p = lu_p + ((u_t-sol_t)/coeff) to delta_p = ((u_t-sol_t)/coeff).
190*/
191      double delta_p = ((u_t-sol_t)/coeff) ;
192      assert(delta_p > origdelta_p) ;
193      rc.setEffectiveness(delta_p-origdelta_p) ;
194      if (rc.effectiveness() > disaggEffectiveness) {
195        rc.setRow(2,index,element,false) ;
196#       if CGL_DEBUG > 0
197        if (debugger) assert(!debugger->invalidCut(rc)); 
198#       endif
199        rowCut.addCutIfNotDuplicate(rc) ;
200      }
201    }
202/*
203  Generate the upper disaggregation cut.
204*/
205    boundChange = origl_t-l_t ;
206    coeff = plusMinus*boundChange ;
207    if (boundChange < 0.0 && origl_t > -1.0e10 &&
208        (sol_t < l_t+boundChange*origdelta_p-primalTolerance_)) {
209      OsiRowCut rc ;
210      rc.setLb(l_t+coeff*lu_p) ;
211      rc.setUb(DBL_MAX) ;
212      index[0] = icol ;
213      element[0] = 1.0 ;
214      index[1] = pndx ;
215      element[1] = coeff ;
216      double delta_p = ((sol_t-l_t)/coeff) ;
217      assert(delta_p > origdelta_p) ;
218      rc.setEffectiveness(delta_p-origdelta_p) ;
219      if (rc.effectiveness() > disaggEffectiveness) {
220        rc.setRow(2,index,element,false) ;
221#       if CGL_DEBUG > 0
222        if (debugger) assert(!debugger->invalidCut(rc)); 
223#       endif
224        rowCut.addCutIfNotDuplicate(rc) ;
225#if 0
226        printf("%d original bounds %g, %g new Lo %g sol= %g int %d sol= %g\n",icol,origl_t,origu_t,colLower[icol],colsol[icol], pndx, colsol[pndx]) ;
227        printf("-1.0 * x(%d) + %g * y(%d) <= %g\n",
228               icol,boundChange,pndx,rc.ub()) ;
229#endif
230      }
231    }
232  }
233
234# if CGL_DEBUG > 0
235  std::cout
236    << "Leaving disaggCuts, "
237    << rowCut.numberCuts() << " cuts." << std::endl ;
238# endif
239  return ;
240}
241
242
243
244
245/*
246  Walk the column for this variable and propagate a bound change on x<j> to
247  the row bounds for rows referencing x<j>. Recall that the row bounds are
248
249    L<i> = SUM{P}(l<j>) + SUM{M}(u<j>)
250    U<i> = SUM{P}(u<j>) + SUM{M}(l<j>)
251
252  for negative coefficients M and positive coefficients P.
253  The update cases are:
254
255    column bound    a<ij>    row bound
256        l<j>         >0        L<i>
257        l<j>         <0        U<i>
258        u<j>         >0        U<i>
259        u<j>         <0        L<i>
260 
261  Movement is used to determine the column bound that's being tightened. It
262  should be > 0 to tighten l<j>, movement < 0 to tighten u<j>.  Once we've
263  selected the row bound to update, the update is always += value*movement.
264*/
265
266bool updateRowBounds (int j, double movement,
267                      const int *const columnStart,
268                      const int *const columnLength,
269                      const int *const row,
270                      const double *const columnElements, 
271                      const double *const rowUpper,
272                      const double *const rowLower,
273                      int &nstackR, int *const stackR, int *const markR,
274                      double *const minR, double *const maxR,
275                      double *const saveMin, double *const saveMax)
276{
277  bool feasible = true ;
278
279  for (int k = columnStart[j] ; k < columnStart[j]+columnLength[j] ; k++) {
280    int irow = row[k] ;
281    double value = columnElements[k] ;
282    double delta = value*movement ;
283/*
284  markR is initialised by tighten2. A code of -1 indicates this row should be
285  processed; -2 says ignore. Other codes should not happen.
286
287  TODO: This assert, and the disabled continue, go back to the change in
288        tighten2 where rows with infinite bounds are handed arbitrary
289        bounds of 1e60 and labelled suitable for further processing. The
290        assert should be removed and the continue clause reinstated.
291        -- lh, 101127 --
292
293  TODO: Replace the original assert (!= -2) with a more effective one.
294        But don't remove the rest of the code structure. I still think I'll
295        end up reintroducing the `ignore' code.  -- lh, 101203 --
296*/
297#   if CGL_DEBUG > 0
298    if (markR[irow] < -1 || markR[irow]%10000 >= nstackR) {
299      std::cout
300        << "Row " << irow << " marked " << markR[irow]
301        << "; expected between -1 and " << nstackR << "." << std::endl ;
302     
303      dump_row(irow,rowLower[irow],rowUpper[irow],minR[irow],maxR[irow],
304               0,false,false,0,0,0,0,0,0,0) ;
305      std::cout << std::endl ;
306    }
307    assert (markR[irow] >= -1 && markR[irow]%10000 < nstackR) ;
308#   endif
309    if (markR[irow] == -1) {
310      stackR[nstackR] = irow ;
311      markR[irow] = nstackR ;
312      saveMin[nstackR] = minR[irow] ;
313      saveMax[nstackR] = maxR[irow] ;
314      nstackR++ ;
315#   if 0
316    } else if (markR[irow] == -2) {
317      continue ;
318#   endif
319    }
320/*
321  LOU_DEBUG: clear count as row bound will change.
322*/
323    markR[irow] = markR[irow]%10000 ;
324/*
325  Case analysis to tighten the appropriate bound.
326*/
327    if (movement > 0.0) {
328      if (value > 0.0) {
329        if (minR[irow] > -1.0e10)
330          minR[irow] += delta ;
331        if (minR[irow] > rowUpper[irow]+1.0e-5) {
332          feasible = false ;
333          break ;
334        }
335      } else {
336        if (maxR[irow] < 1.0e10)
337          maxR[irow] += delta ;
338        if (maxR[irow] < rowLower[irow]-1.0e-5) {
339          feasible = false ;
340          break ;
341        }
342      }
343    } else {
344      if (value > 0.0) {
345        if (maxR[irow] < 1.0e10)
346          maxR[irow] += delta ;
347        if (maxR[irow] < rowLower[irow]-1.0e-5) {
348          feasible = false ;
349          break ;
350        }
351      } else {
352        if (minR[irow] > -1.0e10)
353          minR[irow] += delta ;
354        if (minR[irow] > rowUpper[irow]+1.0e-5) {
355          feasible = false ;
356          break ;
357        }
358      }
359    }
360  }
361  return (feasible) ;
362}
363
364
365/*
366  jjf: clean up djs and solution
367
368  Calculate domains (u<j>-l<j>) for variables, and potential swings in rows. Do
369  some cleanup of reduced costs.
370
371  If CGL_DEBUG > 0, we'll also do some matrix consistency checks.
372*/
373void groomSoln (double direction, double primalTolerance, double *const djs,
374                const double *const colLower, double *const colsol,
375                const double *const colUpper, double *const columnGap,
376                const CoinPackedMatrix *const rowCopy,
377                const int *const rowStartPos,
378                double *const largestNegativeInRow,
379                double *const largestPositiveInRow)
380{
381  int nRows = rowCopy->getNumRows() ;
382  int nCols = rowCopy->getNumCols() ;
383  const int *rowStart = rowCopy->getVectorStarts() ;
384  const int *column = rowCopy->getIndices() ;
385  const double *rowElements = rowCopy->getElements() ;
386
387# ifndef NDEBUG
388  // Note that rowLength is only used in the assert
389  const int *rowLength = rowCopy->getVectorLengths() ;
390# endif
391
392/*
393  Find the largest potential postive and negative swing in a row, where
394  swing is defined as a<ij>(u<j>-l<j>). Check correctness of rowStartPos if
395  we're debugging.
396*/
397  for (int i = 0 ; i < nRows ; i++) {
398
399    assert (rowStart[i]+rowLength[i] == rowStart[i+1]) ;
400
401    int kk ;
402    double value ;
403
404#if CGL_DEBUG > 0
405    // Check correctness of rowStartPos
406    for (kk = rowStart[i] ; kk < rowStart[i+1] ; kk++) {
407      value = rowElements[kk] ;
408      if (value > 0.0) 
409        break ;
410    }
411    assert (rowStartPos[i] == kk) ;
412#endif
413
414    value = 0.0 ;
415    for (kk = rowStart[i] ; kk < rowStartPos[i] ; kk++) {
416      int iColumn = column[kk] ;
417      double gap = CoinMin(1.0e100,colUpper[iColumn]-colLower[iColumn]) ;
418      value = CoinMin(value,gap*rowElements[kk]) ;
419    }
420    largestNegativeInRow[i] = value ;
421    value = 0.0 ;
422    for ( ; kk < rowStart[i+1] ; kk++) {
423      int iColumn = column[kk] ;
424      double gap = CoinMin(1.0e100,colUpper[iColumn]-colLower[iColumn]) ;
425      value = CoinMax(value,gap*rowElements[kk]) ;
426    }
427    largestPositiveInRow[i] = value ;
428  }
429/*
430  Convert reduced costs to minimisation and clean up any that are on the wrong
431  side of zero.
432
433  TODO: Seems like we could move this ahead of the previous loops and then
434        use columnGap[i] to calculate largest[Negative,Positive]InRow.
435        -- lh, 101125 --
436*/
437  for (int i = 0 ; i < nCols ; ++i) {
438    double djValue = djs[i]*direction ;
439    double gap = colUpper[i]-colLower[i] ;
440    if (gap > 1.0e-8) {
441      if (colsol[i] < colLower[i]+primalTolerance) {
442        colsol[i] = colLower[i] ;
443        djs[i] = CoinMax(0.0,djValue) ;
444      } else if (colsol[i] > colUpper[i]-primalTolerance) {
445        colsol[i] = colUpper[i] ;
446        djs[i] = CoinMin(0.0,djValue) ;
447      } else {
448        djs[i] = 0.0 ;
449      }
450    }
451    columnGap[i] = gap-primalTolerance ;
452  }
453
454  return ;
455}
456
457/*
458  Given that the variable being probed has been discovered to be monotone,
459  this method will save the resulting bound changes for integer variables as
460  column cuts.
461
462  One could do the same for continuous variables, but at present we don't.
463  Presumably we need to test against origColLower, origColUpper (the
464  bounds held in the si) to decide if we have real cuts on the original
465  problem. One could argue that we could test against saveL and saveU and
466  save the trouble of repeatedly recording the same column cut.
467
468  Index and elements are scratch arrays.
469
470  Returns true if any of the bound improvements are actually cuts, false
471  otherwise.
472*/
473bool monotoneActions (double primalTolerance_,
474                      const OsiSolverInterface &si,
475                      OsiCuts &cs,
476                      int nstackC, const int *const stackC,
477                      const char *const intVar,
478                      const double *const colLower,
479                      const double *const colsol,
480                      const double *const colUpper,
481                      int *const index, double *const element)
482{
483  OsiColCut cc ;
484  bool anyColumnCuts = false ;
485  int nTot = 0 ;
486  int nFix = 0 ;
487  const double *origColLower = si.getColLower() ;
488  const double *origColUpper = si.getColUpper() ;
489# if CGL_DEBUG > 0
490  int probedVar = stackC[0] ;
491  std::cout
492    << "    " << si.getColName(probedVar) << " (" << probedVar
493    << ") monotone [" << colLower[probedVar] << ", "
494    << colUpper[probedVar] << "] from " << colsol[probedVar]
495    << "." << std::endl ;
496# endif
497/*
498  Stash the improved column bounds, u<j> first, then l<j>. If any of them
499  actually cut off the current solution, indicate that we have cuts.
500*/
501  bool ifCut = false ;
502  for (int istackC = 0 ; istackC < nstackC ; istackC++) {
503    int icol = stackC[istackC] ;
504    if (intVar[icol]) {
505      if (colUpper[icol] < origColUpper[icol]-1.0e-4) {
506        element[nFix] = colUpper[icol] ;
507        index[nFix++] = icol ;
508        if (colsol[icol] > colUpper[icol]+primalTolerance_) {
509          ifCut = true ;
510          anyColumnCuts = true ;
511        }
512      }
513    }
514  }
515  if (nFix) {
516    nTot = nFix ;
517    cc.setUbs(nFix,index,element) ;
518    nFix = 0 ;
519  }
520  for (int istackC = 0 ; istackC < nstackC ; istackC++) {
521    int icol = stackC[istackC] ;
522    if (intVar[icol]) {
523      if (colLower[icol] > origColLower[icol]+1.0e-4) {
524        element[nFix] = colLower[icol] ;
525        index[nFix++] = icol ;
526        if (colsol[icol] < colLower[icol]-primalTolerance_) {
527          ifCut = true ;
528          anyColumnCuts = true ;
529        }
530      }
531    }
532  }
533  if (nFix) {
534    nTot += nFix ;
535    cc.setLbs(nFix,index,element) ;
536  }
537  if (nTot) {
538    if (ifCut) {
539      cc.setEffectiveness(100.0) ;
540    } else {
541      cc.setEffectiveness(1.0e-5) ;
542    }
543#if CGL_DEBUG > 0
544    checkBounds(si,cc) ;
545#endif
546    cs.insert(cc) ;
547  }
548
549  return (anyColumnCuts) ;
550}
551
552void clearStacks (double primalTolerance_,
553                  int nstackC, int *const stackC, int *const markC,
554                  const double *const colLower, const double *const colUpper,
555                  int nstackR, int *const stackR, int *const markR)
556{
557/*
558  Clear off the bound change markers, except for fixed variables. We'll be
559  moving on to another variable.
560*/
561  for (int istackC = 0 ; istackC < nstackC ; istackC++) {
562    int icol = stackC[istackC] ;
563    if (colUpper[icol]-colLower[icol] > primalTolerance_) {
564      markC[icol] &= ~(tightenLower|tightenUpper) ;
565    } else {
566      markC[icol] = tightenLower|tightenUpper ;
567    }
568  }
569  for (int istackR = 0 ; istackR < nstackR ; istackR++) {
570    int irow = stackR[istackR] ;
571    markR[irow] = -1 ;
572  }
573}
574
575
576
577/*
578  This method examines the rows on stackR looking for redundant rows that
579  consist solely of singleton variables (i.e., variables that appear in just
580  this constraint). It then looks to see if any of the variables in the row
581  can be fixed at bound.
582
583  Consider the case where U<i> < b<i> and all unfixed variables in the row
584  are singletons (occur in no other constraint). Given x<j> is a singleton,
585  the reduced cost is c<j> - ya<j> = c<j> - ye<i> = c<j> - y<i>.  But if
586  U<i> < b<i>, the constraint is loose by definition, hence y<i> = 0 and
587  the reduced cost is c<j>. If a<ij> > 0, x<j> should be u<j> in U<i>. If
588  c<j> < 0, minimising will tend to drive x<j> to u<j>. This cannot violate
589  constraint i (U<i> < b<i>) and (since x<j> is a singleton) will have no
590  effect elsewhere. Hence we can fix x<j> at u<j>.
591
592  Do the case analysis, and what you find is that against U<i> we want
593    a<ij> > 0  ==>  c<j> < 0  ==>  push x<j> to u<j>
594    a<ij> < 0  ==>  c<j> > 0  ==>  push x<j> to l<j>
595  and against L<i> we want
596    a<ij> > 0  ==>  c<j> > 0  ==>  push x<j> to l<j>
597    a<ij> < 0  ==>  c<j> < 0  ==>  push x<j> to u<j>
598
599  Extend it one more time by taking the objective direction (max/min) into
600  account (dir = 1.0 for min, -1.0 for max) and you have
601    against U<i> ==> a<ij>*c<j>*dir < 0
602    against L<i> ==> a<ij>*c<j>*dir > 0
603
604  John's original comment for this code was
605    // also see if singletons can go to good objective
606    // Taken out as should be found elsewhere
607    // and has to be original column length
608  but he reinstated it. Apparently there were cases where fixing the probing
609  variable was required to satisfy the condition that all unfixed variables be
610  singletons. Enough of them to justify reinstating this code.
611*/
612
613void singletonRows (int jProbe, double primalTolerance_,
614                    const OsiSolverInterface &si,
615                    const CoinPackedMatrix *rowCopy,
616                    int *const markC,
617                    int &nstackC, int *const stackC,
618                    double *const saveL, double *const saveU,
619                    double *const colUpper, double *const colLower,
620                    double *const colGap,
621                    const int nstackR, const int *const stackR,
622                    const double *const rowUpper,
623                    const double *const rowLower,
624                    const double *const maxR, const double *const minR)
625{
626/*
627  Unpack a few vectors from the row-major matrix.
628*/
629  const CoinBigIndex *rowStart = rowCopy->getVectorStarts() ;
630  const int *column = rowCopy->getIndices() ;
631  const double *rowElements = rowCopy->getElements() ;
632  const int nCols = rowCopy->getNumCols() ;
633/*
634  `Singleton' must be based on the column length in the original system.
635*/
636  const double *objective = si.getObjCoefficients() ;
637  const int *columnLengths = si.getMatrixByCol()->getVectorLengths() ;
638  const double objSense = si.getObjSense() ;
639/*
640  Open a loop to work through the rows on stackR.
641*/
642  for (int istackR = 0 ; istackR < nstackR ; istackR++) {
643    int i = stackR[istackR] ;
644/*
645  Check the gaps. If the constraint is potentially tight in both directions,
646  there's nothing more to do here.
647*/
648    const double uGap = rowUpper[i]-maxR[i] ;
649    const double lGap = minR[i]-rowLower[i] ;
650    if (uGap < primalTolerance_ && lGap < primalTolerance_) continue ;
651/*
652  Scan the row to see if it meets the `all singletons' condition. Again,
653  if this fails, there's nothing more to be done.
654
655  Note that the original code didn't check the probing variable x<p>,
656  because if you're probing a binary variable it's fixed. But for general
657  integer variables a probe does not in general fix the variable.  So we
658  check all variables.
659
660  We should not be executing this method if we have prima facie infeasibility.
661*/
662    bool canFix = true ;
663    for (int kk = rowStart[i] ; kk < rowStart[i+1] ; kk++) {
664      int j = column[kk] ;
665      assert(colUpper[j]-colLower[j] > -primalTolerance_) ;
666      if (colUpper[j] > colLower[j] && columnLengths[j] != 1) {
667        canFix = false ;
668        break ;
669      }
670    }
671    if (!canFix) continue ;
672/*
673  If we've passed the tests, we can look for variables suitable to drive to
674  bound. Work against U<i> first. We're looking for variables with a<ij> > 0
675  that will be naturally driven to u<j>, or variables with a<ij> < 0 that will
676  be naturally driven to l<j>.
677 
678  Don't overwrite the saved bounds if we've
679  tightened this variable already!
680*/
681    if (uGap > primalTolerance_) {
682      for (int kk = rowStart[i] ; kk < rowStart[i+1] ; kk++) {
683        int j = column[kk] ;
684        const double lj = colLower[j] ;
685        const double uj = colUpper[j] ;
686        if (uj > lj) {
687          double value = rowElements[kk] ;
688          if (objSense*objective[j]*value < 0.0)
689          { assert(jProbe != j) ;
690            if (!(markC[j]&(tightenLower|tightenUpper))) {
691              stackC[nstackC] = j ;
692              saveL[nstackC] = lj ;
693              saveU[nstackC] = uj ;
694            }
695            assert(saveU[nstackC] > saveL[nstackC]) ;
696            if (value > 0.0) {
697              colLower[j] = uj ;
698            } else {
699              colUpper[j] = lj ;
700            }
701            markC[j] |= tightenLower|tightenUpper ;
702            colGap[j] = -primalTolerance_ ;
703            assert(nstackC < nCols) ;
704            nstackC++ ;
705          }
706        }
707      }
708    }
709/*
710  And now the identical code, except that we're working against L<i>, hence
711  the sense of the elibigility test is reversed and we want variables with
712  a<ij> > 0 that will be naturally driven to l<j>, or variables with
713  a<ij> < 0 that will be naturally driven to u<j>.
714*/
715    if (lGap > primalTolerance_) {
716      for (int kk = rowStart[i] ; kk < rowStart[i+1] ; kk++) {
717        int j = column[kk] ;
718        const double lj = colLower[j] ;
719        const double uj = colUpper[j] ;
720        if (uj > lj) {
721          double value = rowElements[kk] ;
722          if (objSense*objective[j]*value > 0.0)
723          { assert(jProbe != j) ;
724            if (!(markC[j]&(tightenLower|tightenUpper))) {
725              stackC[nstackC] = j ;
726              saveL[nstackC] = lj ;
727              saveU[nstackC] = uj ;
728            }
729            assert(saveU[nstackC] > saveL[nstackC]) ;
730            if (value < 0.0) {
731              colLower[j] = uj ;
732            } else {
733              colUpper[j] = lj ;
734            }
735            markC[j] |= tightenLower|tightenUpper ;
736            colGap[j] = -primalTolerance_ ;
737            assert(nstackC < nCols) ;
738            nstackC++ ;
739          }
740        }
741      }
742    }
743  }
744  return ;
745}
746
747
748/*
749  Generate coefficient strengthening cuts.
750 
751  Assume we have a binary probe variable x<p> with a<ip> > 0. Assume a
752  down probe.  Assume that U<i> > b<i> before the probe, but now U'<i> < b<i>
753  (i.e., the constraint is redundant against b<i> after the down probe forces
754  x<p> to 0 and reduces U<i> to U'<i>). We would like to have the following in
755  a strengthened constraint a'<ip>x<p> + (stuff) <= b'<i>:
756    * When x<p> = 0, b'<i> = U'<i>  (the lhs can't be larger than U'<i>)
757    * When x<p> = 1, b'<i> = b<i>   (the original value)
758  Define delta = b<i> - U'<i>, b'<i> = b<i>-delta, and  a'<ip> = a<ip>-delta.
759  When x<p> = 1, the delta term on each side cancels and we're left with the
760  original constraint. When x<p> = 0, the rhs tightens to
761  b'<i> = b<i>+delta<i> = U'<i>.
762
763  For an honest derivation that works for binary and general integer
764  variables, see the typeset documentation. As you'd expect, there are four
765  cases (a<ip> < 0, a<ip> > 0) X (up probe, down probe). Assume that a down
766  probe bounds x<p> as [l<p>, u'<p>] and an up probe bounds x<p> as [l'<p>,
767  u<p>]. Then the summary (omitting subscripts to fit the line) is:
768
769  a<ip>  dir    delta             blow'       a'<ip>      b'
770
771    >0    d    b - U'                         a-delta   b-(u'+1)delta
772    >0    u    blow - L'   blow+(l'-1)delta   a+delta
773    <0    d    b - L'      blow-(u'+1)delta   a-delta
774    <0    u    blow - U'                      a+delta   b+(l'-1)delta
775
776  In the code, negating delta for down probes unifies the up and down math.
777
778  TODO: In the original code, this method will not execute if column cuts
779        are present (i.e., bounds were tightened so that some x*<j> may not
780        be feasible). The current code retains those guards.  Clearly it's
781        a problem if x*<p> isn't feasible --- we can't generate a useful
782        cut over floor(x*<p>), ceil(x*<p>) if the interval is outside the
783        polytope. It's not so clear this is an obstacle for other variables,
784        except that the lhs activity may be inaccurate.
785*/
786
787#define STRENGTHEN_PRINT
788void strengthenCoeff (
789                      int jProbe, unsigned int probeDir,
790                      double primalTolerance_,
791                      bool justReplace, bool canReplace,
792                      double needEffectiveness,
793                      const OsiSolverInterface &si,
794                      CglProbingRowCut &rowCut,
795                      const CoinPackedMatrix *rowCopy,
796                      double *const colUpper, double *const colLower,
797                      const double *const colsol,
798                      const int nstackR, const int *const stackR,
799                      const double *const rowUpper,
800                      const double *const rowLower,
801                      const double *const maxR,
802                      const double *const minR,
803                      const int *const realRows,
804                      double *const element,
805                      int *const index,
806                      CglTreeInfo *const info
807                     )
808{
809# if CGL_DEBUG > 0
810  std::cout
811    << "Entering strengthenCoeff, probing "
812    << si.getColName(jProbe) << "(" << jProbe << ") "
813    << ((probeDir == probeDown)?"down":"up") << "." << std::endl ;
814  const OsiRowCutDebugger *debugger = si.getRowCutDebugger() ;
815# endif
816/*
817  Define magic numbers up front.
818    * Limits on `reasonable' coefficients. We don't want coefficients
819      outside this range.
820    * `Infinity' for the rhs. This'll get fixed later.
821    * Effectiveness in the case that a<ip> = 0 and a'<ip> != 0.
822  All copied from the original code.
823*/
824  const double bigCoeff = 1.0e8 ;
825  const double tinyCoeff = 1.0e-12 ;
826  const double rhsInf = 1.0e20 ;
827  const double aipZeroEffectiveness = 1.0e-7 ;
828/*
829  Unpack a few vectors from the row-major matrix.
830*/
831  const CoinBigIndex *rowStart = rowCopy->getVectorStarts() ;
832  const int *column = rowCopy->getIndices() ;
833  const double *rowElements = rowCopy->getElements() ;
834/*
835  Open up an outer loop to walk stackR and look for interesting constraints.
836*/
837  for (int istackR = 0 ; istackR < nstackR ; istackR++) {
838    int irow = stackR[istackR] ;
839    double bi = rowUpper[irow] ;
840    double blowi = rowLower[irow] ;
841/*
842  We can't get anywhere unless probing has made one end or the other of the
843  constraint redundant. If there's no gap, we can't do anything.
844*/
845    double uGap = bi-maxR[irow] ;
846    double lGap = blowi-minR[irow] ;
847    if (uGap < primalTolerance_ && -lGap < primalTolerance_) continue ;
848    bool isRangeCon = ((blowi  > -rhsInf) && (bi < rhsInf)) ;
849/*
850  We'll need the lhs activity, excluding the probed variable, to evaluate the
851  effectiveness of the cut.  Extract the coefficient of the probe variable
852  while we're here; we need it to determine which case we're working.
853
854  TODO: As noted at the head of the routine, we could correct for column cuts
855        by checking bounds when calculating the sum below, if it was worth the
856        effort.   -- lh, 110113 --
857*/
858    double sum = 0.0 ;
859    double aip = 0.0 ;
860    bool aipNonZero = false ;
861    for (CoinBigIndex kk = rowStart[irow] ; kk < rowStart[irow+1] ; kk++) {
862      int k = column[kk] ;
863      double aik = rowElements[kk] ;
864      if (k == jProbe) {
865        aipNonZero = true ;
866        aip = aik ;
867      } else {
868        sum += aik*colsol[k] ;
869      }
870    }
871/*
872  Now figure out which case we're in and do some setup. At the end, see if the
873  coefficient will have reasonable magnitude. If not, move on to the next
874  iteration.
875*/
876    double delta ;
877    double deltaMul ;
878    bool revisebi = true ;
879    if (probeDir == probeDown) {
880      deltaMul = colUpper[jProbe]+1.0 ; 
881      if (aip >= 0) {
882        delta = -uGap ;
883      } else {
884        delta = -lGap ;
885        revisebi = false ;
886      }
887    } else {
888      deltaMul = colLower[jProbe]-1.0 ; 
889      if (aip >= 0) {
890        delta = lGap ;
891        revisebi = false ;
892      } else {
893        delta = uGap ;
894      }
895    }
896    if (CoinAbs(delta) < primalTolerance_ || CoinAbs(delta) > bigCoeff)
897      continue ;
898/*
899  Now decide if we have something that cuts off the current solution. Augment
900  the lhs activity with the contribution for a'<ip>x*<p>, calculate the new
901  rhs value, and compare.
902  As an alternate measure of effectiveness, consider the gap between the
903  current activity and the revised lhs bound, normalised by the gap between
904  the original rhs and the revised lhs bound.
905
906  We'll also generate the cut if we can strengthen it in place and we're not
907  dealing with a range constraint. (The reason for excluding range constraints
908  is that we might try to alter a<ip> against both b<i> and blow<i>. That way
909  lies madness.)
910
911  TODO: It's not clear to me how the first two effectiveness calculations
912        relate to one another. Think more on this -- lh, 110115 --
913*/
914    double aipPrime = aip+delta ;
915    bool aipPrimeNonZero = true ;
916    if (CoinAbs(aipPrime) < tinyCoeff) {
917      aipPrime = 0.0 ;
918      aipPrimeNonZero = false ;
919    }
920    double xp = colsol[jProbe] ;
921    sum += aipPrime*xp ;
922    double biPrime = DBL_MAX ;
923    double blowiPrime = -DBL_MAX ;
924    double effectiveness = 0.0 ;
925    bool genCut = false ;
926    if (revisebi) {
927      biPrime = rowUpper[irow]+deltaMul*delta ;
928      effectiveness = sum-biPrime ;
929      effectiveness = CoinMax(effectiveness,(sum-maxR[irow])/uGap) ;
930    } else {
931      blowiPrime = rowLower[irow]+deltaMul*delta ;
932      effectiveness = blowiPrime-sum ;
933      effectiveness = CoinMax(effectiveness,(minR[irow]-sum)/lGap) ;
934    }
935    if (!aipNonZero && aipPrimeNonZero)
936      effectiveness = CoinMax(effectiveness,aipZeroEffectiveness) ;
937    if (effectiveness > needEffectiveness) genCut = true ;
938    if (info->strengthenRow && !isRangeCon) genCut = true ;
939/*
940  Are we going to generate the cut? If not, on to the next iteration.
941*/
942    if (!genCut) continue ;
943/*
944  Generate the coefficients. Copy whatever's present and plug in aipPrime at
945  the end.
946*/
947    int n = 0 ;
948    int aip_n = -1 ;
949    for (CoinBigIndex kk = rowStart[irow] ; kk < rowStart[irow+1] ; kk++) {
950      int k = column[kk] ;
951      double aik = rowElements[kk] ;
952      if (k == jProbe)
953        aip_n = k ;
954      index[n] = k ;
955      element[n++] = aik ;
956    }
957    if (aip_n < 0) {
958      aip_n = n ;
959      n++ ;
960    }
961    index[aip_n] = jProbe ;
962    element[aip_n] = aipPrime ;
963/*
964  Fill in a cut structure with the cut.
965*/
966    OsiRowCut rc ;
967    rc.setLb(blowiPrime) ;
968    rc.setUb(biPrime) ;
969    rc.setEffectiveness(effectiveness) ;
970    rc.setRow(n,index,element,false) ;
971#   if CGL_DEBUG > 0
972    if (debugger) assert(!debugger->invalidCut(rc)); 
973#   endif
974#   ifdef STRENGTHEN_PRINT
975    if (canReplace && !isRangeCon) {
976      printf("Strengthen Cut (1):\n") ;
977      dump_row(irow,rc.lb(),rc.ub(),
978               nan(""),nan(""),&si,true,false,4,
979               n,index,element,
980               1.0e-10,colLower,colUpper) ;
981      printf("Original Row:\n") ;
982      int k = rowStart[irow] ;
983      dump_row(irow,rowLower[irow],rowUpper[irow],
984               minR[irow],maxR[irow],&si,true,false,4,
985               rowStart[irow+1]-k,&column[k],&rowElements[k],
986               1.0e-10,colLower,colUpper) ;
987    }
988#   endif
989/*
990  If we're in preprocessing, we might try to simply replace the existing
991  constraint (justReplace = true; preprocessing is the only place this will
992  happen). Otherwise, drop the cut into the cut set.
993
994  realRows comes in as a parameter. This is the translation array created if
995  we modified the constraint system during preprocessing in gutsOfGenerateCuts.
996
997  Effectiveness, if we're strengthening in place, seems to be absolute size of
998  coefficients; smaller is better. (Why?)
999
1000  TODO: It seems to me that we can get here with
1001        justReplace = true and canReplace = false, and this condition is
1002        constant over an iteration of the way loop. In other words, we've done
1003        all the work to generate this cut and we knew before we started that
1004        we would discard it here.  -- lh, 110107 --
1005        Put in an assert to see if we ever do all the work, only to reject
1006        the cut.  -- lh, 110115 --
1007*/
1008      int realRow = (canReplace && !isRangeCon)?(irow):(-1) ;
1009      if (realRows && realRow >= 0)
1010        realRow = realRows[realRow] ;
1011      if (!justReplace) {
1012        rowCut.addCutIfNotDuplicate(rc,realRow) ;
1013      } else if (realRow >= 0) {
1014        double effectiveness = 0.0 ;
1015        for (int i = 0 ; i < n ; i++)
1016          effectiveness += CoinAbs(element[i]) ;
1017        if (!info->strengthenRow[realRow] ||
1018            info->strengthenRow[realRow]->effectiveness() > effectiveness) {
1019          delete info->strengthenRow[realRow] ;
1020          rc.setEffectiveness(effectiveness) ;
1021          info->strengthenRow[realRow] = rc.clone() ;
1022        } else {
1023#       if CGL_DEBUG > 0
1024          std::cout
1025            << "INPLACE: rejected on cut coeff magnitude." << std::endl ;
1026#       endif
1027        }
1028      } else {
1029#       if CGL_DEBUG > 0
1030          std::cout
1031            << "INPLACE: rejected because no real row." << std::endl ;
1032#       endif
1033      }
1034    }
1035
1036# if CGL_DEBUG > 0
1037  std::cout
1038    << "Leaving strengthenCoeff, "
1039    << rowCut.numberCuts() << " cuts." << std::endl ;
1040# endif
1041
1042  return ;
1043}
1044
1045
1046
1047
1048// =========================================================
1049
1050
1051
1052/*
1053  jjf: Does probing and adding cuts
1054
1055  Note that this method is heavily commented and instrumented in CbcAnn.
1056
1057  It would appear that probe() has received more attention that probeClique or
1058  probeSlack. Neither of them has the ONE_ARRAY option.
1059*/
1060
1061/*
1062  We're going to probe integer variables, and we're interested in propagating
1063  the effect of each probe (force a variable up or down).
1064
1065  The bulk of the method consists of three nested loops:
1066    * PASS LOOP: makes multiple passes, probing all variables in lookedAt_
1067    * LOOKEDAT LOOP: iterates over the variables in lookedAt_
1068    * PROBE LOOP: probes down/up/down for each variable.
1069
1070  The body of the probe loop runs a bit over 2600 lines and propagates
1071  the effect of forcing the probed variable down or up.
1072
1073    * The start of the body updates the row bounds of affected rows, then
1074      initialises stackC with the probed variable.
1075
1076    * The middle is a 1,000 line loop that does the bulk of propagation.
1077      At the start there's some dead code (PROBING100) associated with
1078      CglTreeProbingInfo. The purpose remains to be determined. This is
1079      followed by six specialised replications of the same functionality:
1080      walk the column of a variable from stackC, and for each affected
1081      row, try to tighten the bounds of other variables in the row. If
1082      we successfully tighten bounds on a variable, walk the column of
1083      that variable, tighten the bounds of affected rows, and place the
1084      variable on stackC.
1085
1086    * The end of the body deals with the result of the probe. At the end
1087      of each iteration, there's a decision: have we proven infeasibility
1088      (INFEASIBILITY) or monotonicity (MONOTONE), or do we need another
1089      iteration (ITERATE). Monotonicity and iteration are large blocks
1090      in themselves.
1091
1092  There was a fair bit of preamble and postamble around the PASS loop.
1093  The preamble code is protected by PROBING_EXTRA_STUFF and has been disabled
1094  since 2006. It appears to find disaggregation cuts.  The postamble code
1095  purports to find bigM constraints. It's clearly a work in progress. I've
1096  moved both out to tmp files.  -- lh, 101203 --
1097
1098  Comments on the data structures:
1099
1100  markC  0: not fixed
1101         1: tightened upper bound
1102         2: tightened lower bound
1103         3: fixed
1104
1105    or perhaps the above should be interpreted as bit codes, along with the
1106    bits used to indicate infinite bounds:
1107
1108       0x0: no bounds tightened
1109       0x1: u<j> tightened
1110       0x2: l<j> tightened
1111       0x4: l<j> < -1e10 (-infty)
1112       0x8: u<j> >  1e10 (+infty)
1113
1114    Note that we'll only tighten a bound once --- this is enforced during
1115    bounds propagation.
1116
1117  stackC      Records columns to be processed. This record is started anew
1118              each time we probe, pushing a variable in a given direction.
1119              The variable being probed is always entry 0 on stackC.
1120              saveL and saveU are correlated with stackC.
1121
1122              This isn't a stack --- it's a record of every variable
1123              processed. nstackC is the tail of the queue, and istackC is
1124              the head. At the end of a probing pass it's scanned to recover
1125              the actual implications that were generated. It's not clear
1126              to me why it's allocated at 2*ncols. Perhaps, at one time,
1127              there were separate entries for upper and lower bound changes?
1128
1129              No, just as I'm thinking I'm almost done, I notice there's some
1130              fancy footwork in the portion of the code that handles the case
1131              where both the up and down probes were feasible. Seems likely
1132              that we're allowing for each direction to process all variables
1133              at most once.
1134
1135              And as I'm working through the code that handles the case
1136              where down and up probes were both feasible, there in the
1137              move singleton code is an assert that nstackC < nCols.
1138
1139              And by the time I finally reach the end, I'm pretty well
1140              convinced that it's an error if we exceed nCols on stackC.
1141              The only thing to worry about is whether singleton processing
1142              could add a variable that's already been added due to
1143              propagation, and a check of the code says the answer is 'no'.
1144
1145       TODO:  When I rewrite this, allocate stackC as nCols and put in some
1146              asserts to catch any violations. There are already asserts
1147              in place at the point where a variable is placed on stackC
1148              during propagation. (Look for `Is there any change to
1149              propagate'.)   -- lh, 101128 --
1150
1151
1152  stackR      Stack rows to be processed. This stack is started anew each
1153              time we probe pushing a variable in a given direction.
1154
1155  minR, maxR, markR  Initialised externally by tighten2. For row i, minR[i] =
1156                     LB(i), maxR[i] = UB(i) (row lhs lower and upper bounds,
1157                     respectively).  markR is set to -1 if there's at least
1158                     one finite lhs bound, -2 if we shouldn't process this
1159                     row (no finite bounds, or too many coefficients).
1160
1161  Then, as we work, markR[i] will be set to point to the entry in stackR for
1162  row i. saveMin and saveMax are correlated with stackR, and hold the
1163  original values for LB(i) and UB(i) (from minR and maxR) while we're
1164  working. As it turns out, however, the the back pointer from markR is never
1165  used.
1166
1167  largestPositiveInRow  for row i, max{j in P} a<ij>(u<i>-l<j>), where P is
1168                        the set of positive coefficients a<ij>
1169  largestNegativeInRow  for row i, min{j in M} a<ij>(u<i>-l<j>), where M is
1170                        the set of negative coefficients a<ij>
1171
1172  element and index are scratch arrays used to construct cuts.
1173
1174  realRows is the translation array created if we modified the constraint
1175  system during preprocessing in gutsOfGenerateCuts.
1176
1177*/
1178int CglProbing::probe( const OsiSolverInterface &si, 
1179                       OsiCuts &cs, 
1180                       double *colLower, double *colUpper, 
1181                       CoinPackedMatrix *rowCopy,
1182                       CoinPackedMatrix *columnCopy,
1183                       const CoinBigIndex *rowStartPos, const int *realRows, 
1184                       const double *rowLower, const double *rowUpper,
1185                       const char *intVar, double *minR, double *maxR, 
1186                       int *markR, 
1187                       CglTreeInfo *info) const
1188
1189{
1190# if CGL_DEBUG > 0
1191  std::cout << "Entering CglProbing::probe." << std::endl ;
1192# endif
1193/*
1194  PREPARATION
1195
1196  All the code down through `PASS LOOP: HEAD' is preparation. Do all of the
1197  setup that will not change over the nested loops that do the work of probing.
1198  Note that rowCopy may have considerably fewer rows than the original system.
1199*/
1200
1201  int nRows = rowCopy->getNumRows() ;
1202  int nCols = rowCopy->getNumCols() ;
1203  int maxStack = info->inTree ? maxStack_ : maxStackRoot_ ;
1204
1205/*
1206  Fiddle maxStack. Odd sort of function:
1207     1 < maxStack <= 25 => 2*maxStack
1208    26 < maxStack <= 50 => 50
1209    51 < maxStack       => maxStack
1210
1211  TODO: Grepping through the code says there's no way that totalTimesCalled_
1212        can become negative. Perhaps in some derived class? More likely a
1213        magic number for something, or else a quick way to disable this
1214        code.  -- lh, 101021 --
1215*/
1216  if ((totalTimesCalled_%10) == -1) {
1217    int newMax = CoinMin(2*maxStack,50) ;
1218    maxStack = CoinMax(newMax,maxStack) ;
1219  }
1220  totalTimesCalled_++ ;
1221
1222# ifdef ONE_ARRAY
1223  assert((DIratio != 0) && ((DIratio&(DIratio-1)) == 0)) ;
1224  int nSpace = 8*nCols+4*nRows+2*maxStack ;
1225  nSpace += (4*nCols+nRows+maxStack+DIratio-1)>>(DIratio-1) ;
1226  double * colsol = new double[nSpace] ;
1227  double * djs = colsol + nCols ;
1228  double * columnGap = djs + nCols ;
1229  double * saveL = columnGap + nCols ;
1230  double * saveU = saveL + 2*nCols ;
1231  double * saveMin = saveU + 2*nCols ;
1232  double * saveMax = saveMin + nRows ;
1233  double * largestPositiveInRow = saveMax + nRows ;
1234  double * largestNegativeInRow = largestPositiveInRow + nRows ;
1235  double * element = largestNegativeInRow + nRows ;
1236  double * lo0 = element + nCols ;
1237  double * up0 = lo0 + maxStack ;
1238  int * markC = reinterpret_cast<int *> (up0+maxStack) ;
1239  int * stackC = markC + nCols ;
1240  int * stackR = stackC + 2*nCols ;
1241  int * index = stackR + nRows ;
1242  int * stackC0 = index + nCols ;
1243# else
1244  double * colsol = new double[nCols] ;
1245  double * djs = new double[nCols] ;
1246  double * columnGap = new double [nCols] ;
1247  double * saveL = new double [2*nCols] ;
1248  double * saveU = new double [2*nCols] ;
1249  double * saveMin = new double [nRows] ;
1250  double * saveMax = new double [nRows] ;
1251  double * largestPositiveInRow = new double [nRows] ;
1252  double * largestNegativeInRow = new double [nRows] ;
1253  double * element = new double[nCols] ;
1254  double * lo0 = new double[maxStack] ;
1255  double * up0 = new double[maxStack] ;
1256  int * markC = new int [nCols] ;
1257  int * stackC = new int [2*nCols] ;
1258  int * stackR = new int [nRows] ;
1259  int * index = new int[nCols] ;
1260  int * stackC0 = new int[maxStack] ;
1261# endif
1262
1263/*
1264  Create a local container to hold the cuts we generate. At the end we'll
1265  transfer these to the container passed as a parameter.
1266
1267  Typically, rowCopy (nRows) will be smaller than the original system in the
1268  si, but it could be one larger (all original constraints plus the objective).
1269  At the root, allow lots of room; even more if this is the first call at the
1270  root. In the tree, much less.
1271
1272  If all we're doing is strengthening rows in place (option 0x40 indicates
1273  we're in preprocessing, a good place to do this), we can only work with
1274  what we have in rowCopy (and we need a translation array).
1275
1276  TODO: There's a problem with realRows here --- if we didn't delete any
1277        rows during preprocessing in gutsOfGenerateCuts, realRows is not
1278        allocated. So we can't strengthen in place in that case. Not likely
1279        the intent.  -- lh, 101209 --
1280*/
1281  bool justReplace = ((info->options&0x40) != 0) && (realRows != NULL) ;
1282  int nRowsFake = 0 ;
1283  if (justReplace) {
1284    nRowsFake = nRows ;
1285  } else {
1286    int nRowsSafe = CoinMin(nRows,si.getNumRows()) ;
1287    nRowsFake = info->inTree ? nRowsSafe/3 : nRowsSafe*10 ;
1288    if (!info->inTree && info->pass == 0) nRowsFake *= 10 ;
1289  }
1290  CglProbingRowCut rowCut(nRowsFake,!info->inTree) ;
1291
1292
1293  // Unpack matrices
1294  const int * column = rowCopy->getIndices() ;
1295  const CoinBigIndex * rowStart = rowCopy->getVectorStarts() ;
1296  const double * rowElements = rowCopy->getElements() ;
1297
1298  const int * row = columnCopy->getIndices() ;
1299  const CoinBigIndex * columnStart = columnCopy->getVectorStarts() ;
1300  const int * columnLength = columnCopy->getVectorLengths(); 
1301  const double * columnElements = columnCopy->getElements() ;
1302
1303/*
1304  Grab the column solution and reduced costs and groom them.
1305*/
1306  CoinMemcpyN(si.getReducedCost(),nCols,djs) ;
1307  CoinMemcpyN(si.getColSolution(),nCols,colsol) ;
1308  double direction = si.getObjSense() ;
1309  groomSoln(direction,primalTolerance_,djs,colLower,colsol,colUpper,columnGap,
1310            rowCopy,rowStartPos,largestNegativeInRow,largestPositiveInRow) ;
1311/*
1312  Determine objective cutoff (minimisation convention).
1313
1314  TODO: Seems to me that this bit of code is a guaranteed fail for a
1315        maximisation problem with no cutoff. It also repeats in a number
1316        of places.  -- lh, 101125 --
1317*/
1318  double cutoff ;
1319  bool cutoff_available = si.getDblParam(OsiDualObjectiveLimit,cutoff) ;
1320  if (!cutoff_available || usingObjective_ < 0) {
1321    cutoff = si.getInfinity() ;
1322  }
1323  cutoff *= direction ;
1324  if (CoinAbs(cutoff) > 1.0e30)
1325    assert (cutoff > 1.0e30) ;
1326  double current = si.getObjValue() ;
1327  current *= direction ;
1328/*
1329  Scan the variables, noting the ones that are fixed or have a bound too large
1330  to be useful.
1331*/
1332  //int nFix = 0 ;
1333  for (int i = 0 ; i < nCols ; i++) {
1334    if (colUpper[i]-colLower[i] < 1.0e-8) {
1335      markC[i] = tightenLower|tightenUpper ;
1336      //nFix++ ;
1337    } else {
1338      markC[i] = 0 ;
1339      if (colUpper[i] > 1.0e10)
1340        markC[i] |= infUpper ;
1341      if (colLower[i] < -1.0e10)
1342        markC[i] |= infLower ;
1343    }
1344  }
1345  //printf("PROBE %d fixed out of %d\n",nFix,nCols) ;
1346
1347/*
1348  jjf: If we are going to replace coefficient then we don't need to be
1349       effective
1350
1351  Seems like this comment does not agree with CglProbingRowCut.addCuts(),
1352  which checks effectiveness before adding a cut to the OsiCuts collection or
1353  entering it into the strenghtenRow array.
1354
1355  But it does agree with the way coefficient strengthening cuts are handled
1356  down in the end-of-iteration processing for the down/up/down probe loop.
1357
1358  It looks like the idea of strengthenRow is that a cut that's simply
1359  strengthening an existing row i is entered in slot i of strengthenRow. It's
1360  left to the client to deal with it on return. I need to get a better idea of
1361  how justReplace is handled, here and in the client.  -- lh, 101125 --
1362*/
1363  //double needEffectiveness = info->strengthenRow ? -1.0e10 : 1.0e-3 ;
1364  double needEffectiveness = info->strengthenRow ? 1.0e-8 : 1.0e-3 ;
1365  if (justReplace && ((info->pass&1) != 0))
1366    needEffectiveness = -1.0e10 ;
1367
1368  double tolerance = 1.0e1*primalTolerance_ ;
1369
1370  /* for both way coding */
1371  int nstackC0 = -1 ;
1372  int nstackR,nstackC ;
1373/*
1374  So, let's assume that the real business starts here, and the previous block
1375  is irrelevant.
1376*/
1377/*
1378  All the initialisation that occurs here is specific to CglTreeProbingInfo.
1379  Sort of begs the question `Why is this handled as a derived class of
1380  CglTreeInfo?' And why (what?) is CglTreeInfo, in the grand scheme of things?
1381
1382  Some grepping about in Cbc says that we might see a CglTreeProbingInfo
1383  object during processing of the root node. My tentative hypothesis is that
1384  it's a vehicle to pass implications back to cbc, where they will be stashed
1385  in a CglImplication object for further use. It's worth mentioning that the
1386  existence of a CglTreeProbingInfo object depends on an independent symbol
1387  (CLIQUE_ANALYSIS) being defined in CbcModel::branchAndBound. It's also worth
1388  mentioning that the code here will core dump if the info object is not a
1389  CglTreeProbingInfo object.
1390
1391  And finally, it's worth mentioning that this code is disabled --- PROBING100
1392  is defined to 0 at the head of CglProbing --- since 080722.
1393
1394  The trivial implementation CglTreeInfo::initializeFixing() always returns 0,
1395  so in effect this block ensures saveFixingInfo is false.
1396
1397  -- lh, 101126 --
1398*/
1399  bool saveFixingInfo = false ;
1400#if PROBING100
1401  CglTreeProbingInfo *probingInfo = dynamic_cast<CglTreeProbingInfo *> (info) ;
1402  const int *backward = NULL ;
1403  const int *integerVariable = NULL ;
1404  const int *toZero = NULL ;
1405  const int *toOne = NULL ;
1406  const fixEntry *fixEntries = NULL ;
1407#endif
1408  if (info->inTree) {
1409#if PROBING100
1410    backward = probingInfo->backward() ;
1411    integerVariable = probingInfo->integerVariable() ;
1412    toZero = probingInfo->toZero() ;
1413    toOne = probingInfo->toOne() ;
1414    fixEntries = probingInfo->fixEntries() ;
1415#endif
1416  } else {
1417    saveFixingInfo = (info->initializeFixing(&si) > 0) ;
1418  }
1419/*
1420  PASS LOOP: HEAD
1421
1422  Major block #2: Main pass loop.
1423
1424  From CglProbingAnn: anyColumnCuts is set only in the case that we've
1425  fixed a variable by probing (i.e., one of the up or down probe resulted
1426  in infeasibility) and that probe entailed column cuts. Once set, it is
1427  never rescinded. In the reworked code, it's set as the return value of
1428  monotoneActions().
1429*/
1430  bool anyColumnCuts = false ;
1431  int ninfeas = 0 ;
1432  int rowCuts = rowCuts_ ;
1433  double disaggEffectiveness = 1.0e-3 ;
1434  int ipass = 0 ;
1435  int nfixed = -1 ;
1436  int maxPass = info->inTree ? maxPass_ : maxPassRoot_ ;
1437
1438  while (ipass < maxPass && nfixed) {
1439    int iLook ;
1440    ipass++ ;
1441    //printf("pass %d\n",ipass) ;
1442    nfixed = 0 ;
1443/*
1444  We went through a fair bit of trouble in gutsOfGenerateCut to determine
1445  the set of variables to be probed and loaded the indices into lookedAt_;
1446  numberThisTime_ reflects that.
1447
1448  The net effect here is that the root gets special treatment
1449  (maxProbeRoot_) and the first pass at the root gets extra special treatment
1450  (numberThisTime_).
1451
1452  See CbcCutGenerator::generateCuts for the origin of magic number 123.
1453  Comments there indicate it's intended to take effect deep in the tree, but
1454  the code here will also work at the root.
1455
1456  Special cases aside, the net effect of the loops is to walk lookedAt_,
1457  promoting every nth variable (n = cutDown) to the front of lookedAt_. When
1458  the loop finishes, the rest of the variables (which were copied off to
1459  stackC) are appended to lookedAt_. If XXXXXX is not defined, we have an
1460  expensive noop.  -- lh, 101126 --
1461*/
1462    int justFix = (!info->inTree && !info->pass) ? -1 : 0 ;
1463    int maxProbe = info->inTree ? maxProbe_ : maxProbeRoot_ ;
1464    if (justFix < 0)
1465      maxProbe = numberThisTime_ ;
1466    if (maxProbe == 123) {
1467      maxProbe = 0 ;
1468      if (!info->inTree) {
1469        if (!info->pass || numberThisTime_ < -100) {
1470          maxProbe = numberThisTime_ ;
1471        } else {
1472          int cutDown = 4 ;
1473          int offset = info->pass%cutDown ;
1474          int i ;
1475          int k = 0 ;
1476          int kk = offset ;
1477          for (i = 0 ; i < numberThisTime_ ; i++) {
1478            if (!kk) {
1479#define XXXXXX
1480#ifdef XXXXXX
1481              lookedAt_[maxProbe] = lookedAt_[i] ;
1482#endif
1483              maxProbe++ ;
1484              kk = cutDown-1 ;
1485            } else {
1486              stackC[k++] = lookedAt_[i] ;
1487              kk-- ;
1488            }
1489          }
1490#ifdef XXXXXX
1491          memcpy(lookedAt_+maxProbe,stackC,k*sizeof(int)) ;
1492#endif
1493        }
1494      } else {
1495        if (numberThisTime_ < 200) {
1496          maxProbe = numberThisTime_ ;
1497        } else {
1498          int cutDown = CoinMax(numberThisTime_/100,4) ;
1499          int offset = info->pass%cutDown ;
1500          int i ;
1501          int k = 0 ;
1502          int kk = offset ;
1503          for (i = 0 ; i < numberThisTime_ ; i++) {
1504            if (!kk) {
1505#ifdef XXXXXX
1506              lookedAt_[maxProbe] = lookedAt_[i] ;
1507#endif
1508              maxProbe++ ;
1509              kk = cutDown-1 ;
1510            } else {
1511              stackC[k++] = lookedAt_[i] ;
1512              kk-- ;
1513            }
1514          }
1515#ifdef XXXXXX
1516          memcpy(lookedAt_+maxProbe,stackC,k*sizeof(int)) ;
1517#endif
1518        }
1519      }
1520    }  // end maxProbe == 123
1521
1522/*
1523  This looks to be an overall limit on probing. It's decremented every time a
1524  variable is popped off stackC for processing.
1525
1526  TODO: PROBING5 would disable this limit; currently inactive. -- lh, 101127 --
1527*/
1528    int leftTotalStack = maxStack*CoinMax(200,maxProbe) ;
1529#ifdef PROBING5
1530    if (!info->inTree&&!info->pass)
1531      leftTotalStack = 1234567890 ;
1532#endif
1533    //printf("maxStack %d maxPass %d numberThisTime %d info pass %d\n",
1534    //   maxStack,maxPass,numberThisTime_,info->pass) ;
1535/*
1536  LOOKEDAT LOOP: HEAD
1537
1538  Main loop to probe each variable in lookedAt_
1539
1540  We're finally ready to get down to the business of probing. Open a loop to
1541  probe each variable in lookedAt_.
1542*/
1543    for (iLook = 0 ; iLook < numberThisTime_ ; iLook++) {
1544      double solval ;
1545      double down ;
1546      double up ;
1547/*
1548  Too successful? Consider bailing out.
1549
1550  If we're generating row cuts, but haven't fixed any variables or we're in
1551  the tree, break from probing.
1552
1553  JustFix < 0 says this is the first pass at the root; for any other
1554  iteration of the pass loop, it'll be initialised to 0.  If it is the first
1555  pass at the root, turn off row cuts, keep on fixing variables, and stop
1556  at the end of this pass. Otherwise, if we haven't fixed any variables, break.
1557
1558  TODO: Otherwise keep going with no changes? That doesn't seem right. The
1559        logic here does not cover all situations.  -- lh, 101126 --
1560*/
1561      if (rowCut.outOfSpace() || leftTotalStack <= 0) {
1562        if (!justFix && (!nfixed || info->inTree)) {
1563#ifdef COIN_DEVELOP
1564          if (!info->inTree)
1565            printf("Exiting a on pass %d, maxProbe %d\n",
1566                   ipass,maxProbe) ;
1567#endif   
1568          break ;
1569        } else if (justFix <= 0) {
1570          if (!info->inTree) {
1571            rowCuts = 0 ;
1572            justFix = 1 ;
1573            disaggEffectiveness = COIN_DBL_MAX ;
1574            needEffectiveness = COIN_DBL_MAX ;
1575            //maxStack=10 ;
1576            maxPass = 1 ;
1577          } else if (!nfixed) {
1578#ifdef COIN_DEVELOP
1579            printf("Exiting b on pass %d, maxProbe %d\n",
1580                   ipass,maxProbe) ;
1581#endif   
1582            break ;
1583          }
1584        }
1585      }
1586/*
1587  awayFromBound isn't quite accurate --- we'll also set awayFromBound = 0
1588  for a general integer variable at an integer value strictly within bounds.
1589*/
1590      int awayFromBound = 1 ;
1591/*
1592  Have a look at the current variable. We're not interested in processing it
1593  if either or both bounds have been improved (0x02|0x01). If the variable has
1594  become fixed, claim both bounds have improved.
1595
1596  TODO: if x<j> is not an integer variable we have big problems. This should
1597        be an assert. -- lh, 101126 --
1598*/
1599      int j = lookedAt_[iLook] ;
1600      solval = colsol[j] ;
1601      down = floor(solval+tolerance) ;
1602      up = ceil(solval-tolerance) ;
1603      if (columnGap[j] < 0.0) markC[j] = tightenUpper|tightenLower ;
1604      if ((markC[j]&(tightenUpper|tightenLower)) != 0 || !intVar[j]) continue ;
1605/*
1606  `Normalize' variables that are near (or past) their bounds.
1607
1608  In normal circumstances, u<j> - l<j> is at least 1, while tol will be
1609  down around 10e-5 (given a primal tolerance of 10e-6, about as large as
1610  you'd normally see). Then
1611
1612  l<j> < l<j>+tol < l<j>+2tol << u<j>-2tol < u<j>-tol < u<j>
1613
1614  For x<j> > u<j>-tol, (x<j>,u<j>,u<j>) => (u<j>-.1, l<j>, u<j>)
1615  For x<j> < l<j>+tol, (x<j>,l<j>,l<j>) => (l<j>+.1, l<j>, u<j>)
1616  For up == down,      (x<j>,v,v)       => (v+.1,v,v+1)
1617
1618  So we're forcing a spread of 1 between up and down, and moving x<j>
1619  to something far enough (.1) from integer to avoid accuracy problems when
1620  calculating movement. A side effect is that primal infeasibility is
1621  corrected.
1622
1623  TODO: Why this structure? Another approach would be to test using up and
1624        down calculated just above. It works out pretty much the same if
1625        you stick with variables of type double. Moving to int would likely
1626        make it faster.  -- lh, 101127 --
1627
1628  TODO: Things will get seriously wierd if (u<j> - l<j>) <= 3tol. For
1629        u<j>-tol < x<j> < l<j>+2tol, values will change. But the symmetric
1630        case u<j>-2tol < x<j> < l<j>+tol will never change values because
1631        of the order of the if statements. Clearly the code is not designed
1632        to cope with this level of pathology.  -- lh, 101127 --
1633 
1634  TODO: It's worth noting that awayFromBound is never used.  -- lh, 101127 --
1635
1636  TODO: It's worth noting that the values set for solval are never used except
1637        in a debug print.  -- lh, 101213 --
1638
1639  TODO: The final case below corresponds to l<j>+1 <= down == up < u<j>-1,
1640        e.g., an interior integer value. Seems like we'd want to probe +1 and
1641        -1. Presumably the code isn't set up to do this, so we arbitrarily
1642        pick the up probe. Could the code be augmented to do the stronger
1643        probe?  -- lh, 101127 --
1644*/
1645      if (solval >= colUpper[j]-tolerance ||
1646          solval <= colLower[j]+tolerance || up == down) {
1647        awayFromBound = 0 ;
1648        if (solval <= colLower[j]+2.0*tolerance) {
1649          solval = colLower[j]+1.0e-1 ;
1650          down = colLower[j] ;
1651          up = down+1.0 ;
1652        } else if (solval >= colUpper[j]-2.0*tolerance) {
1653          solval = colUpper[j]-1.0e-1 ;
1654          up = colUpper[j] ;
1655          down = up-1 ;
1656        } else {
1657          up = down+1.0 ;
1658          solval = down+1.0e-1 ;
1659        }
1660      }
1661      assert (up <= colUpper[j]) ;
1662      assert (down >= colLower[j]) ;
1663      assert (up > down) ;
1664#     if CGL_DEBUG > 0
1665      const double lj = colLower[j] ;
1666      const double uj = colUpper[j] ;
1667      const bool downIsLower = (CoinAbs(down-colLower[j]) < 1.0e-7) ;
1668      const bool upIsUpper = (CoinAbs(up-colUpper[j]) < 1.0e-7) ;
1669#     endif
1670/*
1671  Set up to probe each variable down (1), up (2), down (1).
1672
1673  The notion is that we'll set a bit (1, 2, 4) in the result record if we're
1674  feasible for a particular way. Way is defined as:
1675    1: we're looking at movement between floor(x*<j>) and x*<j>, u<j>
1676    2: we're looking at movement between ceil(x*<j>) and x*<j>, l<j>
1677  As defined here, movement will be negative for way = 1.
1678
1679  Why do we have a second pass with way = 1? Glad you asked. About 1200
1680  lines farther down, it finally becomes apparent that we execute the
1681  final pass only when the initial down probe was feasible, but the up probe
1682  was not. In a nutshell, we're restoring the state produced by the down
1683  probe, which we'll then nail down in the section of code labelled `keep',
1684  a mere 600 lines farther on.
1685*/
1686      unsigned int iway ;
1687      unsigned int way[] = { probeDown, probeUp, probeDown } ;
1688      unsigned int feasValue[] =
1689          { downProbeFeas, upProbeFeas, downProbe2Feas } ;
1690      unsigned int feasRecord = 0 ;
1691      bool notFeasible = false ;
1692      int istackC = 0 ;
1693      int istackR = 0 ;
1694/*
1695  PROBE LOOP: HEAD
1696
1697  Open a loop to probe up and down for the current variable. Get things
1698  started by placing the variable on the propagation queue (stackC).
1699
1700  As with the previous loops (variables in lookedAt_, passes) this loop
1701  extends to the end of major block #2).
1702*/
1703      for (iway = downProbe ; iway < oneProbeTooMany ; iway++) {
1704        stackC[0] = j ;
1705        markC[j] = way[iway] ;
1706/*
1707  Calculate movement given the current direction. Given the setup above,
1708  movement should be at least +/-1.
1709
1710  TODO: Notice we reset up or down yet again. Really, this shouldn't be
1711        necessary. For that matter, why did we bother with awayFromBound
1712        immediately above? -- lh, 101127 --
1713
1714  TODO: The more I think about this, the less happy I am. Seems to me that
1715        neither colUpper or colLower can change during iterations of this
1716        loop. If we discover monotone, we save bounds and break. If we
1717        don't discover monotone, we restore. Let's test that.
1718        -- lh, 110105 --
1719*/
1720        double solMovement ;
1721        double movement ;
1722        int goingToTrueBound = 0 ;
1723
1724#       if CGL_DEBUG > 0
1725        assert(lj == colLower[j]) ;
1726        assert(uj == colUpper[j]) ;
1727        assert(downIsLower == (CoinAbs(down-colLower[j]) < 1.0e-7)) ;
1728        assert(upIsUpper == (CoinAbs(up-colUpper[j]) < 1.0e-7)) ;
1729#       endif
1730
1731        if (way[iway] == probeDown) {
1732          movement = down-colUpper[j] ;
1733          solMovement = down-colsol[j] ;
1734          assert(movement < -0.99999) ;
1735          if (CoinAbs(down-colLower[j]) < 1.0e-7) {
1736            goingToTrueBound = 2 ;
1737            down = colLower[j] ;
1738          }
1739        } else {
1740          movement = up-colLower[j] ;
1741          solMovement = up-colsol[j] ;
1742          assert(movement > 0.99999) ;
1743          if (CoinAbs(up-colUpper[j]) < 1.0e-7) {
1744            goingToTrueBound = 2 ;
1745            up = colUpper[j] ;
1746          }
1747        }
1748/*
1749  Coding for goingToTrueBound is:
1750    0: We're somewhere in the interior of a general integer.
1751    1: We're heading for one of the bounds on a general integer.
1752    2: We're processing a binary variable (u<j>-l<j> = 1 and l<j> = 0).
1753  You can view this next test as correcting for the fact that the previous
1754  code assumes binary variables are all there is.
1755
1756  TODO: Now that I've figured out the math for the constraints generated
1757        by disaggCuts (see typeset documentation), I see what's wrong
1758        (?) here. The disagg cut that's generated relies on the new probe
1759        bound being distance one from the original bound. But (for example)
1760        that's (colUpper-down) = 1 for a down probe.  Not quite what's
1761        being checked for here. But this next test ensures binary variables,
1762        and then the tests are equivalent. Before I willy-nilly change this, I
1763        need to sort out a couple of other places where goingToTrueBound
1764        controls behaviour.  -- lh, 101214 --
1765
1766  TODO: Apropos the above, for example, the test for coefficient strengthening
1767        cuts includes goingToTrueBound != 0  -- lh, 110105 --
1768*/
1769        if (goingToTrueBound && (colUpper[j]-colLower[j] > 1.5 || colLower[j]))
1770          goingToTrueBound = 1 ;
1771/*
1772  If the user doesn't want disaggregation cuts, pretend we're in the interior
1773  of a general integer variable.
1774
1775  TODO: Why is row strengthening limited to binary variables? If we can
1776        figure out what to do, the limitation seems artificial.
1777        -- lh, 101127 --
1778        The test here also says that we can't have coefficient strengthening
1779        without disaggregation. I don't see any reason for this dominance.
1780        -- lh, 110105 --
1781*/
1782        if ((rowCuts&1) == 0)
1783          goingToTrueBound = 0 ;
1784        bool canReplace = info->strengthenRow && (goingToTrueBound == 2) ;
1785
1786#ifdef PRINT_DEBUG
1787        // Alert for general integer with nontrivial range.
1788        // Also last use of solval
1789        if (CoinAbs(movement)>1.01) {
1790          printf("big %d %g %g %g\n",j,colLower[j],solval,colUpper[j]) ;
1791        }
1792#endif
1793/*
1794  Recall that we adjusted the reduced costs (djs) and current objective
1795  (current) to the minimisation sign convention. objVal will accumulate
1796  objective change due to forced bound changes.
1797*/
1798        double objVal = current ;
1799        if (solMovement*djs[j] > 0.0)
1800          objVal += solMovement*djs[j] ;
1801        nstackC = 1 ;
1802        nstackR = 0 ;
1803        saveL[0] = colLower[j] ;
1804        saveU[0] = colUpper[j] ;
1805        assert (saveU[0] > saveL[0]) ;
1806        notFeasible = false ;
1807/*
1808  Recalculate the upper (down probe) or lower (up probe) bound. You can see
1809  from the debug print that arbitrary integers are an afterthought in this
1810  code.
1811
1812  TODO: And why not just set the bounds to down (colUpper) or up (colLower)?
1813        We set movement = down-colUpper just above, and we've taken a lot of
1814        care to groom up and down. This is pointless effort.  -- lh, 101214 --
1815*/
1816        if (movement < 0.0) {
1817          colUpper[j] += movement ;
1818          colUpper[j] = floor(colUpper[j]+0.5) ;
1819          columnGap[j] = colUpper[j]-colLower[j]-primalTolerance_ ;
1820#ifdef PRINT_DEBUG
1821          printf("** Trying %d down to 0\n",j) ;
1822#endif
1823        } else {
1824          colLower[j] += movement ;
1825          colLower[j] = floor(colLower[j]+0.5) ;
1826          columnGap[j] = colUpper[j]-colLower[j]-primalTolerance_ ;
1827#ifdef PRINT_DEBUG
1828          printf("** Trying %d up to 1\n",j) ;
1829#endif
1830        }
1831/*
1832  Is this variable now fixed?
1833*/
1834        if (CoinAbs(colUpper[j]-colLower[j]) < 1.0e-6)
1835          markC[j] = tightenUpper|tightenLower ;
1836/*
1837  Clear the infinite bound bits (infUpper|infLower) and check again. Bounds
1838  may have improved.
1839*/
1840        markC[j] &= ~(infUpper|infLower) ;
1841        if (colUpper[j] > 1.0e10)
1842          markC[j] |= infUpper ;
1843        if (colLower[j] < -1.0e10)
1844          markC[j] |= infLower ;
1845        istackC = 0 ;
1846/*
1847  Update row bounds to reflect the change in variable bound.
1848
1849  TODO: Replacing code to update row bounds with updateRowBounds(). We'll
1850        see how it looks.  -- lh, 101203 --
1851*/
1852
1853  if (!updateRowBounds(j,movement,columnStart,columnLength,row,columnElements,
1854                       rowUpper,rowLower,nstackR,stackR,markR,
1855                       minR,maxR,saveMin,saveMax)) {
1856    notFeasible = true ;
1857    istackC = 1 ;
1858  }
1859/*
1860  PROBE LOOP: BEGIN PROPAGATION
1861
1862  Row bounds are now adjusted for all rows with a<ij> != 0. Time to consider
1863  the effects. nstackC is incremented each time we add a variable, istackC is
1864  incremented each time we finish processing it.
1865
1866  If we lost feasibility above, istackC = nstackC = 1 and this loop will not
1867  execute.
1868
1869  TODO: Is stackC really a stack? Or a queue? And what is the role of
1870        maxStack? stackC is allocated at 2*ncols, but maxStack defaults to
1871        50. -- lh, 101127 --
1872
1873  TODO: stackC is clearly a queue. Allocation strategy is still unclear.
1874        -- lh, 101128 --
1875
1876  jjf:  could be istackC<maxStack?
1877*/
1878        while (istackC < nstackC && nstackC < maxStack) {
1879          leftTotalStack-- ;
1880          int jway ;
1881          int jcol = stackC[istackC] ;
1882          jway = markC[jcol] ;
1883/*
1884  jjf: If not first and fixed then skip
1885
1886  It's clear that x<j>, the probing variable, can be marked as fixed at this
1887  point (happens just above, when the upper or lower bound is adjusted). And
1888  there's an earlier check that avoids probing a variable that's fixed when
1889  plucked from lookedAt_. But ... why leave the if with an empty body?
1890
1891  TODO: Disabled since r118 050225, but seems like the right thing to do. Is it
1892        somehow guaranteed that fixed variables are not stacked?
1893        -- lh, 101127 --
1894
1895  TODO: Based on what's happening down below, it looks like we can encounter
1896        variables here which have had both bounds tightened and are then
1897        pushed onto istackC to propagate the effect of that. This bit of code
1898        should simply go away.  -- lh, 101127 --
1899*/
1900          if (((jway&(tightenLower|tightenUpper)) ==
1901               (tightenLower|tightenUpper)) &&
1902              istackC) {
1903            //istackC++ ;
1904            //continue ;
1905            //printf("fixed %d on stack\n",jcol) ;
1906          }
1907#if PROBING100
1908/*
1909  This block of code has been disabled since r661 080722. I have no notes on
1910  it from my old annotated CglProbing. Given the use of backward, it's tied to
1911  CglTreeProbingInfo.  -- lh, 101127 --
1912*/
1913          if (backward) {
1914            int jColumn = backward[jcol] ;
1915            if (jColumn>=0) {
1916              int nFix=0 ;
1917              // 0-1 see what else could be fixed
1918              if (jway==tightenUpper) {
1919                // fixed to 0
1920                int j ;
1921                for (j = toZero_[jColumn];j<toOne_[jColumn];j++) {
1922                  int kColumn=fixEntry_[j].sequence ;
1923                  kColumn = integerVariable_[kColumn] ;
1924                  bool fixToOne = fixEntry_[j].oneFixed ;
1925                  if (fixToOne) {
1926                    if (colLower[kColumn]==0.0) {
1927                      if (colUpper[kColumn]==1.0) {
1928                        // See if on list
1929                        if (!(markC[kColumn]&(tightenLower|tightenUpper))) {
1930                          if(nStackC<nCols) {
1931                            stackC[nstackC]=kColumn ;
1932                            saveL[nstackC]=colLower[kColumn] ;
1933                            saveU[nstackC]=colUpper[kColumn] ;
1934                            assert (saveU[nstackC]>saveL[nstackC]) ;
1935                            assert (nstackC<nCols) ;
1936                            nstackC++ ;
1937                            markC[kColumn]|=tightenLower ;
1938                            nFix++ ;
1939                          }
1940                        } else if ((markC[kColumn]&(tightenLower|tightenUpper))==tightenUpper) {
1941                          notFeasible = true ;
1942                        }
1943                      } else {
1944                        // infeasible!
1945                        notFeasible = true ;
1946                      }
1947                    }
1948                  } else {
1949                    if (colUpper[kColumn]==1.0) {
1950                      if (colLower[kColumn]==0.0) {
1951                        // See if on list
1952                        if (!(markC[kColumn]&(tightenLower|tightenUpper))) {
1953                          if(nStackC<nCols) {
1954                            stackC[nstackC]=kColumn ;
1955                            saveL[nstackC]=colLower[kColumn] ;
1956                            saveU[nstackC]=colUpper[kColumn] ;
1957                            assert (saveU[nstackC]>saveL[nstackC]) ;
1958                            assert (nstackC<nCols) ;
1959                            nstackC++ ;
1960                            markC[kColumn]|=tightenUpper ;
1961                            nFix++ ;
1962                          }
1963                        } else if ((markC[kColumn]&(tightenLower|tightenUpper))==tightenLower) {
1964                          notFeasible = true ;
1965                        }
1966                      } else {
1967                        // infeasible!
1968                        notFeasible = true ;
1969                      }
1970                    }
1971                  }
1972                }
1973              } else if (jway==tightenLower) {
1974                int j ;
1975                for ( j=toOne_[jColumn];j<toZero_[jColumn+1];j++) {
1976                  int kColumn=fixEntry_[j].sequence ;
1977                  kColumn = integerVariable_[kColumn] ;
1978                  bool fixToOne = fixEntry_[j].oneFixed ;
1979                  if (fixToOne) {
1980                    if (colLower[kColumn]==0.0) {
1981                      if (colUpper[kColumn]==1.0) {
1982                        // See if on list
1983                        if (!(markC[kColumn]&(tightenLower|tightenUpper))) {
1984                          if(nStackC<nCols) {
1985                            stackC[nstackC]=kColumn ;
1986                            saveL[nstackC]=colLower[kColumn] ;
1987                            saveU[nstackC]=colUpper[kColumn] ;
1988                            assert (saveU[nstackC]>saveL[nstackC]) ;
1989                            assert (nstackC<nCols) ;
1990                            nstackC++ ;
1991                            markC[kColumn]|=tightenLower ;
1992                            nFix++ ;
1993                          }
1994                        } else if ((markC[kColumn]&(tightenLower|tightenUpper))==tightenUpper) {
1995                          notFeasible = true ;
1996                        }
1997                      } else {
1998                        // infeasible!
1999                        notFeasible = true ;
2000                      }
2001                    }
2002                  } else {
2003                    if (colUpper[kColumn]==1.0) {
2004                      if (colLower[kColumn]==0.0) {
2005                        // See if on list
2006                        if (!(markC[kColumn]&(tightenLower|tightenUpper))) {
2007                          if(nStackC<nCols) {
2008                            stackC[nstackC]=kColumn ;
2009                            saveL[nstackC]=colLower[kColumn] ;
2010                            saveU[nstackC]=colUpper[kColumn] ;
2011                            assert (saveU[nstackC]>saveL[nstackC]) ;
2012                            assert (nstackC<nCols) ;
2013                            nstackC++ ;
2014                            markC[kColumn]|=tightenUpper ;
2015                            nFix++ ;
2016                          }
2017                        } else if ((markC[kColumn]&(tightenLower|tightenUpper))==tightenLower) {
2018                          notFeasible = true ;
2019                        }
2020                      } else {
2021                        // infeasible!
2022                        notFeasible = true ;
2023                      }
2024                    }
2025                  }
2026                }
2027              }
2028            }
2029          }
2030#endif  // PROBING100 (disabled)
2031/*
2032  PROBE LOOP: WALK COLUMN
2033
2034  Loop to walk the column of a variable popped off stackC.
2035
2036  We've pulled x<jcol> off stackC. We're going to walk the column and process
2037  each row where a<i,jcol> != 0. Note that we've already updated the row
2038  bounds to reflect the change in bound for x<jcol>. In effect, we've stacked
2039  this variable as a surrogate for stacking the rows whose bounds were
2040  changed by the change to bounds on this variable. Now we're going to examine
2041  those rows and see if any look promising for further propagation.
2042
2043*/
2044          for (int k = columnStart[jcol] ;
2045               k < columnStart[jcol]+columnLength[jcol] ; k++) {
2046            if (notFeasible)
2047              break ;
2048            int irow = row[k] ;
2049            // Row already processed with these bounds.
2050            if (markR[irow]/10000 > 0) continue ;
2051
2052            int rStart = rowStart[irow] ;
2053            int rEnd = rowStartPos[irow] ;
2054            double rowUp = rowUpper[irow] ;
2055            double rowUp2 = 0.0 ;
2056            bool doRowUpN ;
2057            bool doRowUpP ;
2058/*
2059  So, is this row promising?
2060
2061  If our current L<i> (minR) is larger than the original b<i> (rowUp), we're
2062  infeasible.
2063
2064  Otherwise, a bit of linear algebra brings the conclusion that if the
2065  largest negative gap (min{j in N} a<ij>(u<j>-l<j>)) in the row is less than
2066  -(b<i>-L<i>), we aren't going to be able to make any progress shrinking the
2067  domains of variables with negative coefficients.  Similarly, if the largest
2068  positive gap (max{j in P} a<ij>(u<j>-l<j>) is greater than b<i>-L<i>, we
2069  can't shrink the domain of any variable with a positive coefficient.
2070
2071  There are analogous conditions using U<i> and blow<i>.
2072
2073  Note that we don't update the largest negative and positive gaps, so as
2074  propagation grinds on, this simple test becomes less accurate. On the other
2075  hand, the gaps (b<i>-L<i>) and (U<i>-blow<i>) are also shrinking. Still, it's
2076  another justification for the strict limits on propagation.
2077
2078  Summary:
2079
2080  minR = SUM{P}(a<ij>l<j>) + SUM{N}(a<ij>u<j>)
2081  maxR = SUM{P}(a<ij>u<j>) + SUM{N}(a<ij>l<j>)
2082
2083  doRowUpN: a<ij> < 0, minR can violate rowUp => increase l<j> (more negative)
2084  doRowLoN: a<ij> < 0, maxR can violate rowLo => decrease u<j> (more positive)
2085  doRowUpP: a<ij> > 0, minR can violate rowUp => decrease u<j> (less positive)
2086  doRowLoP: a<ij> > 0, maxR can violate rowLo => increase l<j> (less negative)
2087
2088  Note that the bound (l<j>, u<j>) that can be tightened in any given
2089  situation is *not* the bound involved in minR or maxR.
2090
2091  For binary variables, of course, `shrink the domain' corresponds to fixing
2092  the variable.
2093*/
2094
2095            if (rowUp < 1.0e10) {
2096              doRowUpN = true ;
2097              doRowUpP = true ;
2098              rowUp2 = rowUp-minR[irow] ;
2099              if (rowUp2 < -primalTolerance_) {
2100                notFeasible = true ;
2101                break ;
2102              } else {
2103                if (rowUp2+largestNegativeInRow[irow] > 0)
2104                  doRowUpN = false ;
2105                if (rowUp2-largestPositiveInRow[irow] > 0)
2106                  doRowUpP = false ;
2107              }
2108            } else {
2109              doRowUpN = false ;
2110              doRowUpP = false ;
2111              rowUp2 = COIN_DBL_MAX ;
2112            }
2113            double rowLo = rowLower[irow] ;
2114            double rowLo2 = 0.0 ;
2115            bool doRowLoN ;
2116            bool doRowLoP ;
2117            if (rowLo > -1.0e10) {
2118              doRowLoN = true ;
2119              doRowLoP = true ;
2120              rowLo2 = rowLo-maxR[irow] ;
2121              if (rowLo2 > primalTolerance_) {
2122                notFeasible = true ;
2123                break ;
2124              } else {
2125                if (rowLo2-largestNegativeInRow[irow] < 0)
2126                  doRowLoN = false ;
2127                if (rowLo2+largestPositiveInRow[irow] < 0)
2128                  doRowLoP = false ;
2129              }
2130            } else {
2131              doRowLoN = false ;
2132              doRowLoP = false ;
2133              rowLo2 = -COIN_DBL_MAX ;
2134            }
2135            markR[irow] += 10000 ;
2136/*
2137  LOU_DEBUG: see if we're processing again without a bound change.
2138            if (markR[irow]/10000 > 1)
2139              std::cout
2140                << "REDUNDANT: processing " << irow
2141                << " for the " << markR[irow]/10000 << " time." << std::endl ;
2142*/
2143/*
2144  PROBE LOOP: WALK ROW
2145
2146  Given the above analysis, work on the columns with negative coefficients.
2147
2148  doRowUpN: a<ij> < 0, minR can violate rowUp => increase l<j> (more negative)
2149  doRowLoN: a<ij> < 0, maxR can violate rowLo => decrease u<j> (more positive)
2150
2151  TODO: The asserts here should be compiled out unless we're debugging.
2152        -- lh, 101210 --
2153*/
2154            if (doRowUpN && doRowLoN) {
2155              //doRowUpN=doRowLoN=false ;
2156              // Start neg values loop
2157              for (int kk = rStart ; kk < rEnd ; kk++) {
2158                int kcol = column[kk] ;
2159                int markIt = markC[kcol] ;
2160                // Skip columns with fixed variables.
2161                if ((markIt&(tightenLower|tightenUpper)) != (tightenLower|tightenUpper)) {
2162                  double value2 = rowElements[kk] ;
2163                  if (colUpper[kcol] <= 1e10)
2164                    assert ((markIt&infUpper) == 0) ;
2165                  else
2166                    assert ((markIt&infUpper) != 0) ;
2167                  if (colLower[kcol] >= -1e10)
2168                    assert ((markIt&infLower) == 0) ;
2169                  else
2170                    assert ((markIt&infLower) != 0) ;
2171                  assert (value2 < 0.0) ;
2172/*
2173  Not every column will be productive. Can we do anything with this one?
2174*/
2175                  double gap = columnGap[kcol]*value2 ;
2176                  bool doUp = (rowUp2+gap < 0.0) ;
2177                  bool doDown = (rowLo2-gap > 0.0) ;
2178                  if (doUp || doDown) {
2179                    double moveUp = 0.0 ;
2180                    double moveDown = 0.0 ;
2181                    double newUpper = -1.0 ;
2182                    double newLower = 1.0 ;
2183/*
2184  doUp attempts to increase the lower bound. The math looks like this:
2185    a<irow,kcol>x<kcol> + (minR[irow]-a<irow,kcol>u<kcol>) <= b<irow>
2186    value2*x<kcol> - value2*u<kcol> <= (b<irow>-minR[irow]) = rowUp2
2187    x<kcol> - u<kcol> >= rowUp2/value2
2188    l<kcol> = u<kcol> + rowUp2/value2
2189  Hence we need u<kcol> to be finite (0x8 not set). We also refuse to tighten
2190  l<kcol> a second time (0x2 not set).
2191
2192  TODO: So, why don't we allow tightening a second time?  Clearly we might
2193        process a row on some other iteration that imposes a tighter bound.
2194        Is this an optimisation for binary variables, or is there some
2195        deeper issue at work in terms of correctness?  -- lh, 101127 --
2196
2197  TODO: Also, it's clear that this bit of code would like to handle
2198        continuous variables, but it's also clear that it does it incorrectly.
2199        When the new lower bound straddles the existing upper bound, that
2200        looks to me to be sufficient justification to fix the variable. But
2201        the code here does just the opposite: merely tightening the bound is
2202        flagged as fixing the variable, with no consideration that we may have
2203        gone infeasible. But if we prove we straddle the upper bound, all we
2204        claim is we've tightened the upper bound!  -- lh, 101127 --
2205*/
2206                    if (doUp && ((markIt&(tightenLower|infUpper)) == 0)) {
2207                      double dbound = colUpper[kcol]+rowUp2/value2 ;
2208                      if (intVar[kcol]) {
2209                        markIt |= tightenLower ;
2210                        newLower = ceil(dbound-primalTolerance_) ;
2211                      } else {
2212                        newLower = dbound ;
2213                        if (newLower+primalTolerance_ > colUpper[kcol] &&
2214                            newLower-primalTolerance_ <= colUpper[kcol]) {
2215                          newLower = colUpper[kcol] ;
2216                          markIt |= tightenLower ;
2217                          //markIt = (tightenUpper|tightenLower) ;
2218                        } else {
2219                          // avoid problems - fix later ?
2220                          markIt = (tightenUpper|tightenLower) ;
2221                        }
2222                      }
2223                      moveUp = newLower-colLower[kcol] ;
2224                    }
2225/*
2226  The same, but attempt to decrease the upper bound by comparison with the
2227  blow for the row.
2228*/
2229                    if (doDown && ((markIt&(tightenUpper|infLower)) == 0)) {
2230                      double dbound = colLower[kcol] + rowLo2/value2 ;
2231                      if (intVar[kcol]) {
2232                        markIt |= tightenUpper ;
2233                        newUpper = floor(dbound+primalTolerance_) ;
2234                      } else {
2235                        newUpper = dbound ;
2236                        if (newUpper-primalTolerance_ < colLower[kcol] &&
2237                            newUpper+primalTolerance_ >= colLower[kcol]) {
2238                          newUpper = colLower[kcol] ;
2239                          markIt |= tightenUpper ;
2240                          //markIt = (tightenUpper|tightenLower) ;
2241                        } else {
2242                          // avoid problems - fix later ?
2243                          markIt = (tightenUpper|tightenLower) ;
2244                        }
2245                      }
2246                      moveDown = newUpper-colUpper[kcol] ;
2247                    }
2248/*
2249  Is there any change to propagate? Assuming we haven't exceeded the limit
2250  for propagating, queue the variable. (If the variable's already on the list,
2251  that's not a problem.)
2252
2253  Note that markIt is initially loaded from markC and whatever bounds are
2254  changed are then or'd in. So changes accumulate.
2255*/
2256                    if (!moveUp && !moveDown)
2257                      continue ;
2258                    bool onList = ((markC[kcol]&(tightenLower|tightenUpper)) != 0) ;
2259                    if (nstackC < 2*maxStack) {
2260                      markC[kcol] = markIt ;
2261                    }
2262                    if (moveUp && (nstackC < 2*maxStack)) {
2263                      if (!onList) {
2264                        stackC[nstackC] = kcol ;
2265                        saveL[nstackC] = colLower[kcol] ;
2266                        saveU[nstackC] = colUpper[kcol] ;
2267                        assert (saveU[nstackC] > saveL[nstackC]) ;
2268                        assert (nstackC < nCols) ;
2269                        nstackC++ ;
2270                        onList = true ;
2271                      }
2272/*
2273  The first assert here says `If we're minimising, and d<j> < 0, then
2274  x<j> must be NBUB, so if the new l<j> > x*<j> = u<j>, we're infeasible.'
2275
2276  TODO: Why haven't we detected infeasibility? Put another way, why isn't
2277        there an assert(notFeasible)?  -- lh, 101127 --
2278
2279  TODO: Seems like the change in obj should be ((new l<j>) - x*<j>)*d<j>, but
2280        moveUp is (newLower-colLower). Was there some guarantee that x*<j>
2281        equals l<j>?  -- lh, 101127 --
2282  TODO: I think this is another symptom of incomplete conversion for general
2283        integers. For binary variables, it's not an issue.  -- lh, 101203 --
2284*/
2285                      if (newLower > colsol[kcol]) {
2286                        if (djs[kcol] < 0.0) {
2287                          assert (newLower > colUpper[kcol]+primalTolerance_) ;
2288                        } else {
2289                          objVal += moveUp*djs[kcol] ;
2290                        }
2291                      }
2292/*
2293  TODO: Up above we've already set newLower = ceil(dbound-primalTolerance_).
2294        It's highly unlikely that ceil(newLower-1.0e4) will be different.
2295        -- lh, 101127 --
2296  TODO: My first inclination is to claim that newLower should never be worse
2297        than the existing bound. But I'd better extend the benefit of the
2298        doubt for a bit. Some rows will produce stronger bounds than others.
2299        -- lh, 101127 --
2300*/
2301                      if (intVar[kcol]) 
2302                        newLower = CoinMax(colLower[kcol],
2303                                           ceil(newLower-1.0e-4)) ;
2304                      colLower[kcol] = newLower ;
2305                      columnGap[kcol] = colUpper[kcol]-newLower-primalTolerance_ ;
2306/*
2307  Propagate the column bound change.
2308
2309  TODO: Notice that we're not setting a variable here when we discover
2310        infeasibility; rather, we're setting the column bound to a ridiculous
2311        value which will be discovered farther down.  -- lh, 101127 --
2312  TODO: Replacing code to update row bounds with updateRowBounds(). We'll
2313        see how it looks.  -- lh, 101203 --
2314*/
2315                      if (CoinAbs(colUpper[kcol]-colLower[kcol]) < 1.0e-6) {
2316                        markC[kcol] = (tightenLower|tightenUpper) ;
2317                      }
2318                      markC[kcol] &= ~(infLower|infUpper) ;
2319                      if (colUpper[kcol] > 1.0e10)
2320                        markC[kcol] |= infUpper ;
2321                      if (colLower[kcol] < -1.0e10)
2322                        markC[kcol] |= infLower ;
2323
2324                      if (!updateRowBounds(kcol,moveUp,
2325                               columnStart,columnLength,row,columnElements,
2326                               rowUpper,rowLower,nstackR,stackR,markR,
2327                               minR,maxR,saveMin,saveMax)) {
2328                        colLower[kcol] = 1.0e50 ;
2329                      }
2330                    }
2331/*
2332  Repeat the whole business, in the down direction. Stack the variable, if
2333  it's not already there, do the infeasibility check / objective update, and
2334  update minR / maxR for affected constraints.
2335*/
2336                    if (moveDown&&nstackC<2*maxStack) {
2337                      if (!onList) {
2338                        stackC[nstackC]=kcol ;
2339                        saveL[nstackC]=colLower[kcol] ;
2340                        saveU[nstackC]=colUpper[kcol] ;
2341                        assert (saveU[nstackC]>saveL[nstackC]) ;
2342                        assert (nstackC<nCols) ;
2343                        nstackC++ ;
2344                        onList=true ;
2345                      }
2346                      if (newUpper<colsol[kcol]) {
2347                        if (djs[kcol]>0.0) {
2348                          /* should be infeasible */
2349                          assert (colLower[kcol]>newUpper+primalTolerance_) ;
2350                        } else {
2351                          objVal += moveDown*djs[kcol] ;
2352                        }
2353                      }
2354                      if (intVar[kcol])
2355                        newUpper = CoinMin(colUpper[kcol],floor(newUpper+1.0e-4)) ;
2356                      colUpper[kcol]=newUpper ;
2357                      columnGap[kcol] = newUpper-colLower[kcol]-primalTolerance_ ;
2358                      if (CoinAbs(colUpper[kcol]-colLower[kcol])<1.0e-6) {
2359                        markC[kcol] = (tightenLower|tightenUpper); // say fixed
2360                      }
2361                      markC[kcol] &= ~(infLower|infUpper) ;
2362                      if (colUpper[kcol]>1.0e10)
2363                        markC[kcol] |= infUpper ;
2364                      if (colLower[kcol]<-1.0e10)
2365                        markC[kcol] |= infLower ;
2366
2367                      if (!updateRowBounds(kcol,moveDown,
2368                               columnStart,columnLength,row,columnElements,
2369                               rowUpper,rowLower,nstackR,stackR,markR,
2370                               minR,maxR,saveMin,saveMax)) {
2371                        colUpper[kcol] = -1.0e50 ;
2372                      }
2373                    }
2374/*
2375  We've propagated the bound changes for this column. Check to see if we
2376  discovered infeasibility. Abort by claiming we've cleared the propagation
2377  stack.
2378*/
2379                    if (colLower[kcol] > colUpper[kcol]+primalTolerance_) {
2380                      notFeasible = true ;
2381                      k = columnStart[jcol]+columnLength[jcol] ;
2382                      istackC = nstackC+1 ;
2383                      break ;
2384                    }
2385                  }   // end if (doUp || doDown) (productive column)
2386                }  // end column not fixed
2387              } // end loop on negative coefficients of this row
2388            } else if (doRowUpN) {
2389/*
2390  The above block propagated change for a variable with a negative coefficient,
2391  where we could work against both the row upper and lower bounds. Now we're
2392  going to do it again, but specialised for the case where we can only work
2393  against the upper bound.
2394*/
2395              // Start neg values loop
2396              for (int kk = rStart ; kk < rEnd ; kk++) {
2397                int kcol =column[kk] ;
2398                int markIt=markC[kcol] ;
2399                if ((markIt&(tightenLower|tightenUpper))!= (tightenLower|tightenUpper)) {
2400                  double value2=rowElements[kk] ;
2401                  double gap = columnGap[kcol]*value2 ;
2402                  if (!(rowUp2 + gap < 0.0))
2403                    continue ;
2404                  double moveUp=0.0 ;
2405                  double newLower=1.0 ;
2406                  if ((markIt&(tightenLower|infUpper))==0) {
2407                    double dbound = colUpper[kcol]+rowUp2/value2 ;
2408                    if (intVar[kcol]) {
2409                      markIt |= tightenLower ;
2410                      newLower = ceil(dbound-primalTolerance_) ;
2411                    } else {
2412                      newLower=dbound ;
2413                      if (newLower+primalTolerance_>colUpper[kcol]&&
2414                          newLower-primalTolerance_<=colUpper[kcol]) {
2415                        newLower=colUpper[kcol] ;
2416                        markIt |= tightenLower ;
2417                        //markIt= (tightenLower|tightenUpper) ;
2418                      } else {
2419                        // avoid problems - fix later ?
2420                        markIt= (tightenLower|tightenUpper) ;
2421                      }
2422                    }
2423                    moveUp = newLower-colLower[kcol] ;
2424                    if (!moveUp)
2425                      continue ;
2426                    bool onList = ((markC[kcol]&(tightenLower|tightenUpper))!=0) ;
2427                    if (nstackC<2*maxStack) {
2428                      markC[kcol] = markIt ;
2429                    }
2430                    if (moveUp&&nstackC<2*maxStack) {
2431                      if (!onList) {
2432                        stackC[nstackC]=kcol ;
2433                        saveL[nstackC]=colLower[kcol] ;
2434                        saveU[nstackC]=colUpper[kcol] ;
2435                        assert (saveU[nstackC]>saveL[nstackC]) ;
2436                        assert (nstackC<nCols) ;
2437                        nstackC++ ;
2438                        onList=true ;
2439                      }
2440                      if (newLower>colsol[kcol]) {
2441                        if (djs[kcol]<0.0) {
2442                          /* should be infeasible */
2443                          assert (newLower>colUpper[kcol]+primalTolerance_) ;
2444                        } else {
2445                          objVal += moveUp*djs[kcol] ;
2446                        }
2447                      }
2448                      if (intVar[kcol]) 
2449                        newLower = CoinMax(colLower[kcol],ceil(newLower-1.0e-4)) ;
2450                      colLower[kcol]=newLower ;
2451                      columnGap[kcol] = colUpper[kcol]-newLower-primalTolerance_ ;
2452                      if (CoinAbs(colUpper[kcol]-colLower[kcol])<1.0e-6) {
2453                        markC[kcol]= (tightenLower|tightenUpper); // say fixed
2454                      }
2455                      markC[kcol] &= ~(infLower|infUpper) ;
2456                      if (colUpper[kcol]>1.0e10)
2457                        markC[kcol] |= infUpper ;
2458                      if (colLower[kcol]<-1.0e10)
2459                        markC[kcol] |= infLower ;
2460
2461                      if (!updateRowBounds(kcol,moveUp,
2462                               columnStart,columnLength,row,columnElements,
2463                               rowUpper,rowLower,nstackR,stackR,markR,
2464                               minR,maxR,saveMin,saveMax)) {
2465                        colLower[kcol] = 1.0e50 ;
2466                      }
2467                    }
2468                    if (colLower[kcol] > colUpper[kcol]+primalTolerance_) {
2469                      notFeasible = true ;
2470                      k = columnStart[jcol]+columnLength[jcol] ;
2471                      istackC = nstackC+1 ;
2472                      break ;
2473                    }
2474                  }
2475                }
2476              } // end big loop rStart->rPos
2477            } else if (doRowLoN) {
2478/*
2479  And yet again, for the case where we can only work against the lower bound.
2480*/
2481              // Start neg values loop
2482              for (int kk = rStart ; kk < rEnd ; kk++) {
2483                int kcol =column[kk] ;
2484                if ((markC[kcol]&(tightenLower|tightenUpper))!= (tightenLower|tightenUpper)) {
2485                  double moveDown=0.0 ;
2486                  double newUpper=-1.0 ;
2487                  double value2=rowElements[kk] ;
2488                  int markIt=markC[kcol] ;
2489                  assert (value2<0.0) ;
2490                  double gap = columnGap[kcol]*value2 ;
2491                  bool doDown = (rowLo2 -gap > 0.0) ;
2492                  if (doDown && ((markIt&(tightenUpper|infLower)) == 0)) {
2493                    double dbound = colLower[kcol] + rowLo2/value2 ;
2494                    if (intVar[kcol]) {
2495                      markIt |= tightenUpper ;
2496                      newUpper = floor(dbound+primalTolerance_) ;
2497                    } else {
2498                      newUpper=dbound ;
2499                      if (newUpper-primalTolerance_<colLower[kcol]&&
2500                          newUpper+primalTolerance_>=colLower[kcol]) {
2501                        newUpper=colLower[kcol] ;
2502                        markIt |= tightenUpper ;
2503                        //markIt= (tightenLower|tightenUpper) ;
2504                      } else {
2505                        // avoid problems - fix later ?
2506                        markIt= (tightenLower|tightenUpper) ;
2507                      }
2508                    }
2509                    moveDown = newUpper-colUpper[kcol] ;
2510                    if (!moveDown)
2511                      continue ;
2512                    bool onList = ((markC[kcol]&(tightenLower|tightenUpper))!=0) ;
2513                    if (nstackC<2*maxStack) {
2514                      markC[kcol] = markIt ;
2515                    }
2516                    if (moveDown&&nstackC<2*maxStack) {
2517                      if (!onList) {
2518                        stackC[nstackC]=kcol ;
2519                        saveL[nstackC]=colLower[kcol] ;
2520                        saveU[nstackC]=colUpper[kcol] ;
2521                        assert (saveU[nstackC]>saveL[nstackC]) ;
2522                        assert (nstackC<nCols) ;
2523                        nstackC++ ;
2524                        onList=true ;
2525                      }
2526                      if (newUpper<colsol[kcol]) {
2527                        if (djs[kcol]>0.0) {
2528                          /* should be infeasible */
2529                          assert (colLower[kcol]>newUpper+primalTolerance_) ;
2530                        } else {
2531                          objVal += moveDown*djs[kcol] ;
2532                        }
2533                      }
2534                      if (intVar[kcol])
2535                        newUpper = CoinMin(colUpper[kcol],floor(newUpper+1.0e-4)) ;
2536                      colUpper[kcol]=newUpper ;
2537                      columnGap[kcol] = newUpper-colLower[kcol]-primalTolerance_ ;
2538                      if (CoinAbs(colUpper[kcol]-colLower[kcol])<1.0e-6) {
2539                        markC[kcol]= (tightenLower|tightenUpper); // say fixed
2540                      }
2541                      markC[kcol] &= ~(infLower|infUpper) ;
2542                      if (colUpper[kcol]>1.0e10)
2543                        markC[kcol] |= infUpper ;
2544                      if (colLower[kcol]<-1.0e10)
2545                        markC[kcol] |= infLower ;
2546
2547                      if (!updateRowBounds(kcol,moveDown,
2548                               columnStart,columnLength,row,columnElements,
2549                               rowUpper,rowLower,nstackR,stackR,markR,
2550                               minR,maxR,saveMin,saveMax)) {
2551                        colUpper[kcol] = -1.0e50 ;
2552                      }
2553                    }
2554                    if (colLower[kcol] > colUpper[kcol]+primalTolerance_) {
2555                      notFeasible = true ;
2556                      k = columnStart[jcol]+columnLength[jcol] ;
2557                      istackC = nstackC+1 ;
2558                      break ;
2559                    }
2560                  }
2561                }
2562              }  // end big loop rStart->rPos
2563            }
2564/*
2565  We've finished working on the negative coefficients of the row. Advance
2566  rStart and rEnd to cover the positive coefficients and repeat the previous
2567  500 lines.
2568*/
2569            rStart = rEnd ;
2570            rEnd = rowStart[irow+1] ;
2571            if (doRowUpP&&doRowLoP) {
2572              //doRowUpP=doRowLoP=false ;
2573              // Start pos values loop
2574              for (int kk =rStart;kk<rEnd;kk++) {
2575                int kcol=column[kk] ;
2576                int markIt=markC[kcol] ;
2577                if ((markIt&(tightenLower|tightenUpper))!= (tightenLower|tightenUpper)) {
2578                  double value2=rowElements[kk] ;
2579                  assert (value2 > 0.0) ;
2580                  /* positive element */
2581                  double gap = columnGap[kcol]*value2 ;
2582                  bool doDown = (rowLo2 + gap > 0.0) ;
2583                  bool doUp = (rowUp2 - gap < 0.0) ;
2584                  if (doDown||doUp) {
2585                    double moveUp=0.0 ;
2586                    double moveDown=0.0 ;
2587                    double newUpper=-1.0 ;
2588                    double newLower=1.0 ;
2589                    if (doDown && ((markIt&(tightenLower|infUpper)) == 0)) {
2590                      double dbound = colUpper[kcol] + rowLo2/value2 ;
2591                      if (intVar[kcol]) {
2592                        markIt |= tightenLower ;
2593                        newLower = ceil(dbound-primalTolerance_) ;
2594                      } else {
2595                        newLower=dbound ;
2596                        if (newLower+primalTolerance_>colUpper[kcol]&&
2597                            newLower-primalTolerance_<=colUpper[kcol]) {
2598                          newLower=colUpper[kcol] ;
2599                          markIt |= tightenLower ;
2600                          //markIt= (tightenLower|tightenUpper) ;
2601                        } else {
2602                          // avoid problems - fix later ?
2603                          markIt= (tightenLower|tightenUpper) ;
2604                        }
2605                      }
2606                      moveUp = newLower-colLower[kcol] ;
2607                    }
2608                    if (doUp && ((markIt&(tightenUpper|infLower)) == 0)) {
2609                      double dbound = colLower[kcol] + rowUp2/value2 ;
2610                      if (intVar[kcol]) {
2611                        markIt |= tightenUpper ;
2612                        newUpper = floor(dbound+primalTolerance_) ;
2613                      } else {
2614                        newUpper=dbound ;
2615                        if (newUpper-primalTolerance_<colLower[kcol]&&
2616                            newUpper+primalTolerance_>=colLower[kcol]) {
2617                          newUpper=colLower[kcol] ;
2618                          markIt |= tightenUpper ;
2619                          //markIt= (tightenLower|tightenUpper) ;
2620                        } else {
2621                          // avoid problems - fix later ?
2622                          markIt= (tightenLower|tightenUpper) ;
2623                        }
2624                      }
2625                      moveDown = newUpper-colUpper[kcol] ;
2626                    }
2627                    if (!moveUp&&!moveDown)
2628                      continue ;
2629                    bool onList = ((markC[kcol]&(tightenLower|tightenUpper))!=0) ;
2630                    if (nstackC<2*maxStack) {
2631                      markC[kcol] = markIt ;
2632                    }
2633                    if (moveUp&&nstackC<2*maxStack) {
2634                      if (!onList) {
2635                        stackC[nstackC]=kcol ;
2636                        saveL[nstackC]=colLower[kcol] ;
2637                        saveU[nstackC]=colUpper[kcol] ;
2638                        assert (saveU[nstackC]>saveL[nstackC]) ;
2639                        assert (nstackC<nCols) ;
2640                        nstackC++ ;
2641                        onList=true ;
2642                      }
2643                      if (newLower>colsol[kcol]) {
2644                        if (djs[kcol]<0.0) {
2645                          /* should be infeasible */
2646                          assert (newLower>colUpper[kcol]+primalTolerance_) ;
2647                        } else {
2648                          objVal += moveUp*djs[kcol] ;
2649                        }
2650                      }
2651                      if (intVar[kcol]) 
2652                        newLower = CoinMax(colLower[kcol],ceil(newLower-1.0e-4)) ;
2653                      colLower[kcol]=newLower ;
2654                      columnGap[kcol] = colUpper[kcol]-newLower-primalTolerance_ ;
2655                      if (CoinAbs(colUpper[kcol]-colLower[kcol])<1.0e-6) {
2656                        markC[kcol]= (tightenLower|tightenUpper); // say fixed
2657                      }
2658                      markC[kcol] &= ~(infLower|infUpper) ;
2659                      if (colUpper[kcol]>1.0e10)
2660                        markC[kcol] |= infUpper ;
2661                      if (colLower[kcol]<-1.0e10)
2662                        markC[kcol] |= infLower ;
2663
2664                      if (!updateRowBounds(kcol,moveUp,
2665                               columnStart,columnLength,row,columnElements,
2666                               rowUpper,rowLower,nstackR,stackR,markR,
2667                               minR,maxR,saveMin,saveMax)) {
2668                        colLower[kcol] = 1.0e50 ;
2669                      }
2670                    }
2671                    if (moveDown&&nstackC<2*maxStack) {
2672                      if (!onList) {
2673                        stackC[nstackC]=kcol ;
2674                        saveL[nstackC]=colLower[kcol] ;
2675                        saveU[nstackC]=colUpper[kcol] ;
2676                        assert (saveU[nstackC]>saveL[nstackC]) ;
2677                        assert (nstackC<nCols) ;
2678                        nstackC++ ;
2679                        onList=true ;
2680                      }
2681                      if (newUpper<colsol[kcol]) {
2682                        if (djs[kcol]>0.0) {
2683                          /* should be infeasible */
2684                          assert (colLower[kcol]>newUpper+primalTolerance_) ;
2685                        } else {
2686                          objVal += moveDown*djs[kcol] ;
2687                        }
2688                      }
2689                      if (intVar[kcol])
2690                        newUpper = CoinMin(colUpper[kcol],floor(newUpper+1.0e-4)) ;
2691                      colUpper[kcol]=newUpper ;
2692                      columnGap[kcol] = newUpper-colLower[kcol]-primalTolerance_ ;
2693                      if (CoinAbs(colUpper[kcol]-colLower[kcol])<1.0e-6) {
2694                        markC[kcol]= (tightenLower|tightenUpper); // say fixed
2695                      }
2696                      markC[kcol] &= ~(infLower|infUpper) ;
2697                      if (colUpper[kcol]>1.0e10)
2698                        markC[kcol] |= infUpper ;
2699                      if (colLower[kcol]<-1.0e10)
2700                        markC[kcol] |= infLower ;
2701
2702                      if (!updateRowBounds(kcol,moveDown,
2703                               columnStart,columnLength,row,columnElements,
2704                               rowUpper,rowLower,nstackR,stackR,markR,
2705                               minR,maxR,saveMin,saveMax)) {
2706                        colUpper[kcol] = -1.0e50 ;
2707                      }
2708                    }
2709                    if (colLower[kcol]>colUpper[kcol]+primalTolerance_) {
2710                      notFeasible = true; ;
2711                      k=columnStart[jcol]+columnLength[jcol] ;
2712                      istackC=nstackC+1 ;
2713                      break ;
2714                    }
2715                  }
2716                }
2717              } // end big loop rPos->rEnd
2718            } else if (doRowUpP) {
2719              // Start pos values loop
2720              for (int kk =rStart;kk<rEnd;kk++) {
2721                int kcol =column[kk] ;
2722                int markIt=markC[kcol] ;
2723                if ((markIt&(tightenLower|tightenUpper))!= (tightenLower|tightenUpper)) {
2724                  double value2=rowElements[kk] ;
2725                  assert (value2 > 0.0) ;
2726                  /* positive element */
2727                  double gap = columnGap[kcol]*value2 ;
2728                  bool doUp = (rowUp2 - gap < 0.0) ;
2729                  if (doUp && ((markIt&(tightenUpper|infLower)) == 0)) {
2730                    double newUpper=-1.0 ;
2731                    double dbound = colLower[kcol] + rowUp2/value2 ;
2732                    if (intVar[kcol]) {
2733                      markIt |= tightenUpper ;
2734                      newUpper = floor(dbound+primalTolerance_) ;
2735                    } else {
2736                      newUpper=dbound ;
2737                      if (newUpper-primalTolerance_<colLower[kcol]&&
2738                          newUpper+primalTolerance_>=colLower[kcol]) {
2739                        newUpper=colLower[kcol] ;
2740                        markIt |= tightenUpper ;
2741                        //markIt= (tightenLower|tightenUpper) ;
2742                      } else {
2743                        // avoid problems - fix later ?
2744                        markIt= (tightenLower|tightenUpper) ;
2745                      }
2746                    }
2747                    double moveDown = newUpper-colUpper[kcol] ;
2748                    if (!moveDown)
2749                      continue ;
2750                    bool onList = ((markC[kcol]&(tightenLower|tightenUpper))!=0) ;
2751                    if (nstackC<2*maxStack) {
2752                      markC[kcol] = markIt ;
2753                    }
2754                    if (nstackC<2*maxStack) {
2755                      if (!onList) {
2756                        stackC[nstackC]=kcol ;
2757                        saveL[nstackC]=colLower[kcol] ;
2758                        saveU[nstackC]=colUpper[kcol] ;
2759                        assert (saveU[nstackC]>saveL[nstackC]) ;
2760                        assert (nstackC<nCols) ;
2761                        nstackC++ ;
2762                        onList=true ;
2763                      }
2764                      if (newUpper<colsol[kcol]) {
2765                        if (djs[kcol]>0.0) {
2766                          /* should be infeasible */
2767                          assert (colLower[kcol]>newUpper+primalTolerance_) ;
2768                        } else {
2769                          objVal += moveDown*djs[kcol] ;
2770                        }
2771                      }
2772                      if (intVar[kcol])
2773                        newUpper = CoinMin(colUpper[kcol],floor(newUpper+1.0e-4)) ;
2774                      colUpper[kcol]=newUpper ;
2775                      columnGap[kcol] = newUpper-colLower[kcol]-primalTolerance_ ;
2776                      if (CoinAbs(colUpper[kcol]-colLower[kcol])<1.0e-6) {
2777                        markC[kcol]= (tightenLower|tightenUpper); // say fixed
2778                      }
2779                      markC[kcol] &= ~(infLower|infUpper) ;
2780                      if (colUpper[kcol]>1.0e10)
2781                        markC[kcol] |= infUpper ;
2782                      if (colLower[kcol]<-1.0e10)
2783                        markC[kcol] |= infLower ;
2784
2785                      if (!updateRowBounds(kcol,moveDown,
2786                               columnStart,columnLength,row,columnElements,
2787                               rowUpper,rowLower,nstackR,stackR,markR,
2788                               minR,maxR,saveMin,saveMax)) {
2789                        colUpper[kcol] = -1.0e50 ;
2790                      }
2791                    }
2792                    if (colLower[kcol]>colUpper[kcol]+primalTolerance_) {
2793                      notFeasible = true ;
2794                      k=columnStart[jcol]+columnLength[jcol] ;
2795                      istackC=nstackC+1 ;
2796                      break ;
2797                    }
2798                  }
2799                }
2800              } // end big loop rPos->rEnd
2801            } else if (doRowLoP) {
2802              // Start pos values loop
2803              for (int kk =rStart;kk<rEnd;kk++) {
2804                int kcol =column[kk] ;
2805                if ((markC[kcol]&(tightenLower|tightenUpper))!= (tightenLower|tightenUpper)) {
2806                  double value2=rowElements[kk] ;
2807                  int markIt=markC[kcol] ;
2808                  assert (value2 > 0.0) ;
2809                  /* positive element */
2810                  double gap = columnGap[kcol]*value2 ;
2811                  bool doDown = (rowLo2 +gap > 0.0) ;
2812                  if (doDown&&(markIt&(tightenLower|infUpper))==0) {
2813                    double newLower=1.0 ;
2814                    double dbound = colUpper[kcol] + rowLo2/value2 ;
2815                    if (intVar[kcol]) {
2816                      markIt |= tightenLower ;
2817                      newLower = ceil(dbound-primalTolerance_) ;
2818                    } else {
2819                      newLower=dbound ;
2820                      if (newLower+primalTolerance_>colUpper[kcol]&&
2821                          newLower-primalTolerance_<=colUpper[kcol]) {
2822                        newLower=colUpper[kcol] ;
2823                        markIt |= tightenLower ;
2824                        //markIt= (tightenLower|tightenUpper) ;
2825                      } else {
2826                        // avoid problems - fix later ?
2827                        markIt= (tightenLower|tightenUpper) ;
2828                        }
2829                    }
2830                    double moveUp = newLower-colLower[kcol] ;
2831                    if (!moveUp)
2832                      continue ;
2833                    bool onList = ((markC[kcol]&(tightenLower|tightenUpper))!=0) ;
2834                    if (nstackC<2*maxStack) {
2835                      markC[kcol] = markIt ;
2836                    }
2837                    if (nstackC<2*maxStack) {
2838                      if (!onList) {
2839                        stackC[nstackC]=kcol ;
2840                        saveL[nstackC]=colLower[kcol] ;
2841                        saveU[nstackC]=colUpper[kcol] ;
2842                        assert (saveU[nstackC]>saveL[nstackC]) ;
2843                        assert (nstackC<nCols) ;
2844                        nstackC++ ;
2845                        onList=true ;
2846                      }
2847                      if (newLower>colsol[kcol]) {
2848                        if (djs[kcol]<0.0) {
2849                          /* should be infeasible */
2850                          assert (newLower>colUpper[kcol]+primalTolerance_) ;
2851                        } else {
2852                          objVal += moveUp*djs[kcol] ;
2853                        }
2854                      }
2855                      if (intVar[kcol]) 
2856                        newLower = CoinMax(colLower[kcol],ceil(newLower-1.0e-4)) ;
2857                      colLower[kcol]=newLower ;
2858                      columnGap[kcol] = colUpper[kcol]-newLower-primalTolerance_ ;
2859                      if (CoinAbs(colUpper[kcol]-colLower[kcol])<1.0e-6) {
2860                        markC[kcol]= (tightenLower|tightenUpper); // say fixed
2861                      }
2862                      markC[kcol] &= ~(infLower|infUpper) ;
2863                      if (colUpper[kcol]>1.0e10)
2864                        markC[kcol] |= infUpper ;
2865                      if (colLower[kcol]<-1.0e10)
2866                        markC[kcol] |= infLower ;
2867
2868                      if (!updateRowBounds(kcol,moveUp,
2869                               columnStart,columnLength,row,columnElements,
2870                               rowUpper,rowLower,nstackR,stackR,markR,
2871                               minR,maxR,saveMin,saveMax)) {
2872                        colLower[kcol] = 1.0e50 ;
2873                      }
2874                    }
2875                    if (colLower[kcol]>colUpper[kcol]+primalTolerance_) {
2876                      notFeasible = true ;
2877                      k=columnStart[jcol]+columnLength[jcol] ;
2878                      istackC=nstackC+1 ;
2879                      break ;
2880                    }
2881                  }
2882                }
2883              }      // end big loop rPos->rEnd
2884            }    // end processing of positive coefficients of row.
2885          }   // end loop to walk the column of a variable popped off stackC
2886          istackC++ ;
2887        }  // end stackC processing loop
2888/*
2889  PROBE LOOP: END PROPAGATION
2890
2891  End propagation of probe in one direction for a single variable. Hard to
2892  believe 1,000 lines of code can achieve so little.
2893*/
2894/*
2895  PROBE LOOP: INFEASIBILITY
2896
2897  Feasibility check. Primal infeasibility is noted in notFeasible; we still
2898  need to test against the objective bound. If we have shown up and down
2899  infeasibility, we're infeasible, period. Break from the probe loop and
2900  terminate the lookedAt and pass loops by forcing their iteration counts past
2901  the maximum.
2902
2903  If we're not feasible, then we surely can't drive the variable to bound,
2904  so reset goingToTrueBound.
2905
2906  TODO: I'm coming to think forcing goingToTrueBound = 0 is really an overload
2907        to suppress various post-probe activities (cut generation, singleton
2908        motion, etc) that should not be performed when the probe shows
2909        infeasible. It might be that this is not the best approach.
2910        -- lh, 101216 --
2911*/
2912        if (notFeasible || (objVal > cutoff)) {
2913#         if CGL_DEBUG > 1
2914          std::cout << "  Column " << j << " infeasible: " ;
2915          if (notFeasible && (objVal > cutoff))
2916            std::cout << "primal bounds and objective cutoff" ;
2917          else if (notFeasible)
2918            std::cout << "primal bounds" ;
2919          else
2920            std::cout << "objective cutoff" ;
2921          std::cout << "." << std::endl ;
2922#         endif
2923          notFeasible = true ;
2924          if (iway == upProbe && feasRecord == 0) {
2925            ninfeas = 1 ;
2926            j = nCols-1 ;
2927            iLook = numberThisTime_ ;
2928            ipass = maxPass ;
2929            break ;
2930          }
2931          goingToTrueBound = 0 ;
2932        } else {
2933          feasRecord |= feasValue[iway] ; 
2934        }
2935/*
2936  Save the implications? Currently restricted to binary variables. If info
2937  were a CglTreeProbing object, fixes() would save the implication in the
2938  arrays kept there, returning false only when space is exceeded.
2939 
2940  TODO: This would be one of the places that will fail if fixes() and
2941        initializeFixing() are removed from CglTreeInfo.  -- lh, 101127 --
2942*/
2943        if (!notFeasible && saveFixingInfo) {
2944          // save fixing info
2945          assert (j == stackC[0]) ;
2946          int toValue = (way[iway] == probeDown) ? -1 : +1 ;
2947          for (istackC = 1 ; istackC < nstackC ; istackC++) {
2948            int icol = stackC[istackC] ;
2949            // for now back to just 0-1
2950            if (colUpper[icol]-colLower[icol]<1.0e-12 &&
2951                !saveL[istackC] && saveU[istackC]==1.0) {
2952              assert(saveL[istackC] == colLower[icol] ||
2953                     saveU[istackC] == colUpper[icol]) ;
2954              saveFixingInfo =
2955                  info->fixes(j,toValue,icol,
2956                              (colLower[icol] == saveL[istackC])) ;
2957            }
2958          }
2959        }
2960/*
2961  PROBE LOOP: MONOTONE (BREAK AND KEEP)
2962
2963  We can't prove infeasibility. What's left? If we're at the end of the up
2964  probe, we may know enough to prove monotonicity:
2965
2966  * If this is the up probe and we were feasible on this probe but not on
2967    the down pass, we're monotone.
2968  * If this is the second down probe, we were feasible on down, infeasible on
2969    up (monotone) and repeated the down probe to recreate the bounds.
2970
2971  Either way, execute monotoneActions to capture the state and move on to the
2972  next variable. Terminate the probe loop by forcing iway = oneProbeTooMany.
2973
2974  This if-then-else extends to the end of the down/up/down probe loop.
2975
2976  TODO: There's something not right about this `second down pass' business.
2977        Why do we copy stackC to stackC0 on a feasible first pass? Surely that
2978        could be used to restore the state of the first down pass.
2979        -- lh, 101128 --
2980*/
2981        if (iway == downProbe2 ||
2982            (iway == upProbe && feasRecord == upProbeFeas)) {
2983          nfixed++ ;
2984          bool retVal = monotoneActions(primalTolerance_,si,cs,
2985                                        nstackC,stackC,intVar,
2986                                        colLower,colsol,colUpper,
2987                                        index,element) ;
2988          if (retVal) anyColumnCuts = true ;
2989          clearStacks(primalTolerance_,nstackC,stackC,markC,colLower,colUpper,
2990                      nstackR,stackR,markR) ;
2991          break ;
2992        }
2993/*
2994  PROBE LOOP: ITERATE
2995
2996  If we reach here, we may iterate the probe loop.
2997
2998  Cases remaining include
2999
3000  a) End of the first down probe; in this case we will iterate regardless of
3001     feasibility.
3002  b) End of the up probe, up and down feasible; in this case we will not
3003     iterate (we're done, cannot prove monotonicity).
3004  c) End of the up probe, down feasible, up infeasible; in this case we will
3005     iterate to restore the down feasible state.
3006
3007  TODO: Unless I miss my guess, large chunks of this block of code will be
3008        replicated in the code that handles case b), as we generate cuts
3009        related to forcing the variable up.  Further inspection confirms
3010        this. -- lh, 101128 --
3011
3012  TODO: I'm becoming more convinced that the second down probe iteration is
3013        unnecessary. Given that we have stackC0, lo0, and up0, seems like
3014        calling monotoneActions with lo0 = colLower, up0 = colUpper, and
3015        stackC0 = stackC should do the trick.  -- lh, 101216 --
3016
3017  TODO: As I get the control flow sorted out a bit better, I should try to
3018        separate b) (which does not iterate) from a) and c).
3019
3020  The next block deals with the end of the first down probe. If the probe
3021  is feasible, attempt to tighten column bounds with singletonRows, process
3022  stackC to generate disaggregation cuts, and copy the tightened bounds
3023  to stackC0, lo0, and up0 for use after the up probe.  Feasible or not,
3024  we'll need to go on to the up probe, so restore the original bounds.
3025
3026  Note that GoingToTrueBound was set to 0 way back in initial setup unless
3027  the user requested disaggregation cuts (rowCuts&0x01).
3028
3029  I've rearranged the control flow here, moving the call to singletonRows
3030  ahead of the code that captures column bound changes to stackC0, et al.
3031  This makes it possible to use the same code at the end of case b) (down and
3032  up probe both feasible). Also, singletonRows is just another way to tighten
3033  column bounds. In its original position, it was subject to the same
3034  restrictions as coefficient strengthening (coefficient strengthening enabled
3035  and goingToTrueBound for the probing variable). I don't see any reason to
3036  keep them (but it might be good to add a separate control bit precisely for
3037  singletonRows).  -- lh, 110113 --
3038*/
3039        if (iway == downProbe) {
3040          if (!notFeasible) {
3041/*
3042  Attempt to tighten singleton bounds and generate disaggregation cuts.
3043*/
3044            singletonRows(j,primalTolerance_,si,rowCopy,markC,
3045                          nstackC,stackC,saveL,saveU,
3046                          colUpper,colLower,columnGap,
3047                          nstackR,stackR,rowUpper,rowLower,maxR,minR) ;
3048            if (goingToTrueBound == 2 && !justReplace) {
3049              disaggCuts(nstackC,probeDown,primalTolerance_,
3050                         disaggEffectiveness,si,rowCut,stackC,colsol,
3051                         colUpper,colLower,saveU,saveL,index,element) ;
3052            }
3053/*
3054  Copy off the bound changes from the down probe -- we'll need them after the
3055  up probe.
3056*/
3057            nstackC0 = CoinMin(nstackC,maxStack) ;
3058            for (istackC = 0 ; istackC < nstackC0 ; istackC++) {
3059              int icol = stackC[istackC] ;
3060              stackC0[istackC] = icol ;
3061              lo0[istackC] = colLower[icol] ;
3062              up0[istackC] = colUpper[icol] ;
3063            }
3064          } else {
3065            nstackC0 = 0 ;
3066          }
3067/*
3068  Now reset column bounds to their original values in preparation for the
3069  up probe.
3070*/
3071          for (istackC = nstackC-1 ; istackC >= 0 ; istackC--) {
3072            int icol = stackC[istackC] ;
3073            double oldU = saveU[istackC] ;
3074            double oldL = saveL[istackC] ;
3075            colUpper[icol] = oldU ;
3076            colLower[icol] = oldL ;
3077            columnGap[icol] = oldU-oldL-primalTolerance_ ;
3078            markC[icol] = 0 ;
3079            if (oldU > 1.0e10)
3080              markC[icol] |= infUpper ;
3081            if (oldL < -1.0e10)
3082              markC[icol] |= infLower ;
3083          }
3084/*
3085  And now for the rows. Remember, we're still finishing off the first down
3086  probe and readying for the next iteration (up probe). The final three lines
3087  (of about 300) actually restore the row bounds. The rest is cut generation.
3088*/
3089          if ((rowCuts&0x02) != 0 && goingToTrueBound && !anyColumnCuts)
3090            strengthenCoeff(j,iway,primalTolerance_,justReplace,canReplace,
3091                            needEffectiveness,si,rowCut,
3092                            rowCopy,colUpper,colLower,colsol,nstackR,stackR,
3093                            rowUpper,rowLower,maxR,minR,realRows,
3094                            element,index,info) ;
3095/*
3096  Restore row bounds.
3097*/
3098          for (istackR = 0 ; istackR < nstackR ; istackR++) {
3099            int irow = stackR[istackR] ;
3100            minR[irow] = saveMin[istackR] ;
3101            maxR[irow] = saveMax[istackR] ;
3102            markR[irow] = -1 ;
3103          }
3104        }    // end of code to iterate after first down pass
3105/*
3106  PROBE LOOP: DOWN AND UP FEASIBLE
3107
3108  Again, not quite the full story. Cases remaining:
3109  b) End of up probe, up and down feasible; in this case we will not
3110     iterate.
3111  c) End of up probe, down feasible, up infeasible; in this case we will
3112     iterate to restore the down feasible state.
3113
3114  The next block deals with case b). We don't want to iterate the
3115  down/up/down loop, so force iway to oneProbeTooMany. As usual, the move
3116  singleton code consists of two nearly identical blocks, one working off
3117  U(i), the next off L(i).
3118
3119  TODO: The conditions guarding the move singleton code here are different than
3120        those for pass 0: here there's an added guard using strengthenRow.
3121        Just a failure to change all instances? Nor do I see how John's
3122        comment (`just at root') is enforced here, unless it's from outside
3123        via rowCuts or strengthenRow.
3124
3125        The check for any column cuts is not present, so presumably we can
3126        do singleton motion in the presence of column cuts. The check for gap
3127        does not have the high end (gap < 1.0e8) test.
3128
3129        Note also that here the test for cut generation is moved outside the
3130        loop that iterates over rows. Presumably that's because in the
3131        previous case, the loop also restored markR, etc.
3132
3133        Note that the code that handles the bound change is considerably
3134        different than the previous instance in pass 0. We're adding the
3135        changes to stackC rather than stackC0.
3136        -- lh, 101128 --
3137
3138  Moved the call to disaggCuts up from the common code for b) and c), because
3139  c) (up infeasible) sets goingToTrueBound to 0 and we can't execute
3140  disaggCuts. Similarly for strengthenCoeff.
3141*/
3142          else {
3143          if (iway == upProbe &&
3144              feasRecord == (downProbeFeas|upProbeFeas)) {
3145            iway = oneProbeTooMany ;
3146
3147            singletonRows(j,primalTolerance_,si,rowCopy,markC,
3148                          nstackC,stackC,saveL,saveU,
3149                          colUpper,colLower,columnGap,
3150                          nstackR,stackR,rowUpper,rowLower,maxR,minR) ;
3151            if (goingToTrueBound == 2 && !justReplace) {
3152              disaggCuts(nstackC,probeUp,primalTolerance_,
3153                         disaggEffectiveness,si,
3154                         rowCut,stackC,colsol,colUpper,colLower,saveU,saveL,
3155                         index,element) ;
3156            if ((rowCuts&0x02) != 0 && goingToTrueBound && !anyColumnCuts)
3157              strengthenCoeff(j,iway,primalTolerance_,justReplace,canReplace,
3158                              needEffectiveness,si,rowCut,
3159                              rowCopy,colUpper,colLower,colsol,nstackR,stackR,
3160                              rowUpper,rowLower,maxR,minR,realRows,
3161                              element,index,info) ;
3162            }
3163/*
3164  jjf: point back to stack
3165
3166  We're done playing with singletons. Get to the real work here. We don't need
3167  markC any more; coopt it for a cross-reference, var index -> stackC index.
3168  The +1000 offset allows us to distinguish invalid entries (not all variables
3169  are on stackC).
3170*/
3171            for (istackC = nstackC-1 ; istackC >= 0 ; istackC--) {
3172              int icol = stackC[istackC] ;
3173              markC[icol] = istackC+1000 ;
3174            }
3175            OsiColCut cc ;
3176            int nTot = 0, nFix = 0, nInt = 0 ;
3177            bool ifCut = false ;
3178/*
3179  See if we have bounds improvement. Given a lower bound ld<j> from the
3180  down probe, a bound lu<j> from the up probe, and the original bound lo<j>,
3181  the bound we want is max(min(ld<j>,lu<j>),lo<j>). Use a conservative
3182  tolerance.
3183*/
3184            for (istackC = 1 ; istackC < nstackC0 ; istackC++) {
3185              int icol = stackC0[istackC] ;
3186              int istackC1 = markC[icol]-1000 ;
3187              if (istackC1 >= 0) {
3188                if (CoinMin(lo0[istackC],colLower[icol]) >
3189                                              saveL[istackC1]+1.0e-4) {
3190                  saveL[istackC1] = CoinMin(lo0[istackC],colLower[icol]) ;
3191                  if (intVar[icol] /*||!info->inTree*/ ) {
3192                    element[nFix] = saveL[istackC1] ;
3193                    index[nFix++] = icol ;
3194                    nInt++ ;
3195                    if (colsol[icol] < saveL[istackC1]-primalTolerance_)
3196                      ifCut = true ;
3197                  }
3198                  nfixed++ ;
3199                } 
3200              }
3201            }
3202            if (nFix) {
3203              nTot = nFix ;
3204              cc.setLbs(nFix,index,element) ;
3205              nFix = 0 ;
3206            }
3207/*
3208  Repeat for upper bounds. There's an odd bit of asymmetry here; this loop
3209  tries to generate a cut if there's no bound improvement.
3210
3211  TODO: What, pray tell, is a `two cut'?  Binary variables only, it seems.
3212        Judging from the conditions, we're looking for a variable that's
3213        fixed at 0 on a down probe and fixed at 1 on an up probe, or vice
3214        versa.
3215        -- lh, 101128 --
3216*/
3217            for (istackC = 1 ; istackC < nstackC0 ; istackC++) {
3218              int icol = stackC0[istackC] ;
3219              int istackC1 = markC[icol]-1000 ;
3220              if (istackC1 >= 0) {
3221                if (CoinMax(up0[istackC],colUpper[icol]) <
3222                                              saveU[istackC1]-1.0e-4) {
3223                  saveU[istackC1] = CoinMax(up0[istackC],colUpper[icol]) ;
3224                  if (intVar[icol] /*||!info->inTree*/ ) {
3225                    element[nFix] = saveU[istackC1] ;
3226                    index[nFix++] = icol ;
3227                    nInt++ ;
3228                    if (colsol[icol] > saveU[istackC1]+primalTolerance_)
3229                      ifCut = true ;
3230                  }
3231                  nfixed++ ;
3232                } else if (!info->inTree &&
3233                           saveL[0] == 0.0 && saveU[0] == 1.0) {
3234                  // See if can do two cut
3235                  double upperWhenDown = up0[istackC] ;
3236                  double lowerWhenDown = lo0[istackC] ;
3237                  double upperWhenUp = colUpper[icol] ;
3238                  double lowerWhenUp = colLower[icol] ;
3239                  double upperOriginal = saveU[istackC1] ;
3240                  double lowerOriginal = saveL[istackC1] ;
3241                  if (upperWhenDown < lowerOriginal+1.0e-12 &&
3242                      lowerWhenUp > upperOriginal-1.0e-12) {
3243                    OsiRowCut rc ;
3244                    rc.setLb(lowerOriginal) ;
3245                    rc.setUb(lowerOriginal) ;
3246                    rc.setEffectiveness(1.0e-5) ;
3247                    int index[2] ;
3248                    double element[2] ;
3249                    index[0] = j ;
3250                    index[1] = icol ;
3251                    element[0] = -(upperOriginal-lowerOriginal) ;
3252/*
3253  jjf: If zero then - must have been fixed without noticing!
3254
3255  TODO: Uh, fixed without noticing?! That's an algorithm problem.
3256        -- lh, 101128 --
3257*/
3258                    if (CoinAbs(element[0]) > 1.0e-8) {
3259                      element[1] = 1.0 ;
3260                      rc.setRow(2,index,element,false) ;
3261                      cs.insert(rc) ;
3262                    }
3263                  } else if (upperWhenUp < lowerOriginal+1.0e-12 &&
3264                             lowerWhenDown > upperOriginal-1.0e-12) {
3265                    OsiRowCut rc ;
3266                    rc.setLb(upperOriginal) ;
3267                    rc.setUb(upperOriginal) ;
3268                    rc.setEffectiveness(1.0e-5) ;
3269                    int index[2] ;
3270                    double element[2] ;
3271                    index[0] = j ;
3272                    index[1] = icol ;
3273                    element[0] = upperOriginal-lowerOriginal ;
3274                    element[1] = 1.0 ;
3275                    rc.setRow(2,index,element,false) ;
3276                    cs.insert(rc) ;
3277                  } 
3278                }
3279              }
3280            }    // end loop to update upper bounds (w/ two cuts)
3281            if (nFix) {
3282              nTot+=nFix ;
3283              cc.setUbs(nFix,index,element) ;
3284            }
3285            // could tighten continuous as well
3286            if (nInt) {
3287              if (ifCut) {
3288                cc.setEffectiveness(100.0) ;
3289              } else {
3290                cc.setEffectiveness(1.0e-5) ;
3291              }
3292#if CGL_DEBUG > 0
3293              checkBounds(si,cc) ;
3294#endif
3295              cs.insert(cc) ;
3296            }
3297          }    // end of code to handle up and down feasible.
3298/*
3299  PROBE LOOP: DOWN FEASIBLE, UP INFEASIBLE
3300
3301  The only remaining case is down feasible, up infeasible, in which case we
3302  need to reset to the initial state and run the down probe again. Surely
3303  there must be a better way?
3304*/
3305            else {
3306            goingToTrueBound = 0 ;
3307          }
3308/*
3309  PROBE LOOP: RESTORE
3310
3311  And now the code to reset to the initial state.
3312*/
3313          /* restore all */
3314          for (istackC = nstackC-1 ; istackC >= 0 ; istackC--) {
3315            int icol = stackC[istackC] ;
3316            double oldU = saveU[istackC] ;
3317            double oldL = saveL[istackC] ;
3318/*
3319  TODO: This next bit differs in subtle ways from the restore at the end
3320        of the down probe. Here we force markC to show fixed if the bounds are
3321        within 1.0e-4 of one another, a fairly broad tolerance; code for
3322        the pass 0 restore does not restore fixed status.  Also,
3323        we're not distinguishing here between actions for down feasible,
3324        up feasible, vs. down feasible, up infeasible.  -- lh, 101128 --
3325
3326        I still don't see any reason for the restore here to differ from
3327        the restore after the down probe.  -- lh, 110113 --
3328*/
3329            colUpper[icol] = oldU ;
3330            colLower[icol] = oldL ;
3331            columnGap[icol] = oldU-oldL-primalTolerance_ ;
3332            if (oldU > oldL+1.0e-4) {
3333              markC[icol] = 0 ;
3334              if (oldU > 1.0e10)
3335                markC[icol] |= infUpper ;
3336              if (oldL < -1.0e10)
3337                markC[icol] |= infLower ;
3338            } else {
3339              markC[icol] = (tightenUpper|tightenLower) ;
3340            }
3341          }
3342/*
3343  End of column restore. On to rows.
3344*/
3345          for (istackR = 0 ; istackR < nstackR ; istackR++) {
3346            int irow = stackR[istackR] ;
3347            minR[irow] = saveMin[istackR] ;
3348            maxR[irow] = saveMax[istackR] ;
3349            markR[irow] = -1 ;
3350          }
3351        }       // end of PROBE: DOWN AND UP FEASIBLE
3352      }     // PROBE LOOP: END
3353    }    // LOOKEDAT LOOP: END
3354  }   // PASS LOOP: END
3355/*
3356  Free up the space we've been using.
3357*/
3358#ifndef ONE_ARRAY
3359  delete [] stackC0 ;
3360  delete [] lo0 ;
3361  delete [] up0 ;
3362  delete [] columnGap ;
3363  delete [] markC ;
3364  delete [] stackC ;
3365  delete [] stackR ;
3366  delete [] saveL ;
3367  delete [] saveU ;
3368  delete [] saveMin ;
3369  delete [] saveMax ;
3370  delete [] index ;
3371  delete [] element ;
3372  delete [] djs ;
3373  delete [] largestPositiveInRow ;
3374  delete [] largestNegativeInRow ;
3375#endif
3376  delete [] colsol ;
3377/*
3378  jjf:  Add in row cuts
3379
3380  Transfer row cuts from the local container into something that can escape
3381  this method. If we're not in `just replace' mode, simply copy them over to
3382  the container (cs) passed as a parameter using addCuts. If we're in `just
3383  replace' mode, ? need to explore further ?  -- lh, 101125 --
3384
3385  TODO: If there's any possibility of realRows[i] < 0, this code will break
3386        trying to read strengthenRow[]!
3387        -- lh, 101128 --
3388*/
3389  if (!ninfeas) {
3390    if (!justReplace) {
3391      rowCut.addCuts(cs,info->strengthenRow,info->pass) ;
3392    } else {
3393      for (int i = 0 ; i < nRows ; i++) {
3394        int realRow = realRows[i] ;
3395        OsiRowCut *cut = info->strengthenRow[realRow] ;
3396        if (realRow >= 0 && cut) {
3397#ifdef CLP_INVESTIGATE
3398          printf("Row %d, real row %d effectiveness %g\n",
3399                 i,realRow,cut->effectiveness()) ;
3400#endif
3401          cs.insert(cut) ;
3402        }
3403      }
3404    }
3405  }
3406#if 0
3407/*
3408  Debug print.
3409*/
3410  {
3411    int numberRowCutsAfter = cs.sizeRowCuts() ;
3412    int k ;
3413    for (k = 0;k<numberRowCutsAfter;k++) {
3414      OsiRowCut thisCut = cs.rowCut(k) ;
3415      printf("Cut %d is %g <=",k,thisCut.lb()) ;
3416      int n=thisCut.row().getNumElements() ;
3417      const int * column = thisCut.row().getIndices() ;
3418      const double * element = thisCut.row().getElements() ;
3419      assert (n) ;
3420      for (int i=0;i<n;i++) {
3421        printf(" %g*x%d",element[i],column[i]) ;
3422      }
3423      printf(" <= %g\n",thisCut.ub()) ;
3424    }
3425  }
3426#endif
3427
3428# if CGL_DEBUG > 0
3429  std::cout
3430    << "Leaving CglProbing::probe, ninfeas " << ninfeas
3431    << "." << std::endl ;
3432# endif
3433
3434  return (ninfeas) ;
3435}
Note: See TracBrowser for help on using the repository browser.