source: pico/stable/1.3/src/milp/pico/BCMip.h @ 5996

Revision 5996, 29.4 KB checked in by wehart, 5 years ago (diff)

Merged revisions 5282-5995 via svnmerge from
 https://software.sandia.gov/svn/public/acro/pico/trunk

........

r5295 | jdsiiro | 2007-12-17 17:48:22 -0700 (Mon, 17 Dec 2007) | 4 lines


Propagating utilib::exception_mngr::exit_fn API change.


See log message for utilib r1429.

........

r5298 | jberry | 2008-02-15 12:49:25 -0700 (Fri, 15 Feb 2008) | 3 lines


Prettified some of the debug output.

........

r5311 | wehart | 2008-02-17 22:30:53 -0700 (Sun, 17 Feb 2008) | 2 lines


Misc change to account for the setup of pyacro.

........

r5332 | wehart | 2008-03-02 13:10:17 -0700 (Sun, 02 Mar 2008) | 2 lines


Update given reor of Python.

........

r5344 | wehart | 2008-03-26 20:32:45 -0700 (Wed, 26 Mar 2008) | 3 lines


Creating the 'ipconvert' command, which doesn't rely on
the PICO MIP solvers.

........

r5345 | wehart | 2008-03-27 09:09:21 -0700 (Thu, 27 Mar 2008) | 3 lines


Rework of ipconvert command line to support --output option
and to enable the specification of multi-file inputs.

........

r5375 | lafisk | 2008-05-23 13:52:36 -0600 (Fri, 23 May 2008) | 2 lines


Add link flags needed when --enable-static-executables

........

r5429 | jeckstei | 2008-06-13 15:09:14 -0600 (Fri, 13 Jun 2008) | 10 lines


Fixed some bugs in paralllel MIP heuristic invocation, although one
big one remains -- parBCMipNode::boundComputationGuts is not getting
called.


In parallel runs, the output summary information about MIP heuristics
is now clean and informative, instead a mishmosh of interleaved
printouts from each individual processor.

........

r5448 | caphill | 2008-06-22 09:08:09 -0600 (Sun, 22 Jun 2008) | 8 lines


Fixing a problem in reading val files in an application created by
gen_milp_app. This code infers the data type. If the first values
for a parameter are integer, and then there are doubles, subsequent
reads of the first "integer" variables fail. This fix stores data
as both integer and double during read-in. There may be a more
elegant solution, but this works.

........

r5458 | caphill | 2008-07-01 08:44:15 -0600 (Tue, 01 Jul 2008) | 4 lines


Fixing a bug in generation of dense solutions. It was causing bombs when running
PICO from ampl.

........

r5459 | caphill | 2008-07-02 05:08:29 -0600 (Wed, 02 Jul 2008) | 8 lines


Bug fix for SOS. When SOS overlap, we should not generally fix them
when we discover a half-infeasible problem during pseudocost
initialization. This will generally affect the sets it overlaps,
especially if it fully determines the set (fixes a 1). If PICO
chooses to (or must) then branch on a collaterally-modified SOS, this
creates chaos.

........

r5471 | caphill | 2008-07-16 15:32:16 -0600 (Wed, 16 Jul 2008) | 8 lines


Fix to reading in .val files for gen_milp_app-generated derivative
applications. The mechanism for splitting pieces from a line
(index tuples and value) was dropping the negative sign from
exponents of numbers input in scientific notation. So every
number that looked like 123e-45 became 123e45, which can be
quite different from the original.

........

r5493 | wehart | 2008-08-04 16:03:30 -0600 (Mon, 04 Aug 2008) | 2 lines


Changes due to reorg of coopr.

........

r5512 | wehart | 2008-08-17 14:53:48 -0600 (Sun, 17 Aug 2008) | 8 lines


A complete reimplementation of the classes used to support PICO
customization in SUCASA. The biggest change is that setupMappedMILP is
no longer supported! Now, the SUCASA script in Coopr is responsible for
managing the creating and reading of the *.map file!


These changes touched some of the core PICO files only to remove
the functionality that was needed for setupMappedMILP.

........

r5513 | wehart | 2008-08-17 21:13:27 -0600 (Sun, 17 Aug 2008) | 2 lines


Fixing various bugs in the sucasa customization classes.

........

r5518 | jwatson | 2008-08-22 14:12:10 -0600 (Fri, 22 Aug 2008) | 2 lines


Some minor structural changes and cleanup to the feasibility pump code base.

........

r5519 | jwatson | 2008-08-22 15:24:45 -0600 (Fri, 22 Aug 2008) | 1 line


Fixed some reporting output in the feasibility pump, and added cycle checking (via hash vectors) for iteration solutions to the core algorithm.

........

r5520 | jwatson | 2008-08-22 23:43:45 -0600 (Fri, 22 Aug 2008) | 1 line


Enhanced cycle detection code in feasibility pump.

........

r5522 | wehart | 2008-08-25 14:29:26 -0600 (Mon, 25 Aug 2008) | 2 lines


Misc fix.

........

r5523 | wehart | 2008-08-25 15:32:43 -0600 (Mon, 25 Aug 2008) | 3 lines


Misc changes to optimize the insertion process when initially building
the maps used for symbols.

........

r5525 | wehart | 2008-08-26 17:27:22 -0600 (Tue, 26 Aug 2008) | 4 lines


Initial rework of SUCASA objects to support initialization from *.val files.


This may change substantially...

........

r5531 | wehart | 2008-08-30 08:03:07 -0600 (Sat, 30 Aug 2008) | 6 lines


An initial rework of the MILPSymbol concept to unify the MILPSymbol and MILPSet logic. This change is needed to support sets of sets.


I'm going to rework this further, but I thought I'd archive this version for
future reference. My next revision of this is going to involve UTILIB Any
objects, which will enable a much more dynamic syntax for set-of-sets.

........

r5553 | wehart | 2008-09-08 00:05:27 -0600 (Mon, 08 Sep 2008) | 7 lines


A major rework of the MILPSymbol functionality. This rework
combines the MILPSymbol and MILPSet classes, and the MILPSet class
has been deleted.


This functionality does not impact core PICO builds, but it is used in
the SUCASA customized PICO solvers.

........

r5568 | wehart | 2008-09-15 14:20:02 -0600 (Mon, 15 Sep 2008) | 2 lines


Hacking out a build problem. For some reason, this builds find with g++ on tofu...?

........

r5572 | jberry | 2008-09-17 11:19:29 -0600 (Wed, 17 Sep 2008) | 3 lines


Fixes relating to the use of initialSolve() versus resolve() in the OsiSolverInterface?; the latter was *not* behaving correctly on a good number of benchmarks in FP iterations >= 2, so I changed it to always call "initialSolve()" and added in warm-start behaviors.

........

r5583 | wehart | 2008-09-18 13:51:27 -0600 (Thu, 18 Sep 2008) | 4 lines


Adding pack/unpack capabilities to MILPSymbol objects.


Deleting old development directories.

........

r5588 | jwatson | 2008-09-20 13:38:04 -0600 (Sat, 20 Sep 2008) | 7 lines


Minor fix involving issues relating to run-time-limit termination. For NSR8K,
the root solve takes longer than an allocated run-time of 10 minutes => negative
run-time limits were showing up in the output of the core FP procedure. I now check
for a positive time limit prior to calling the core procedure, and bail with
a message if it is negative. No functional change - just reporting.

........

r5589 | jwatson | 2008-09-20 21:39:44 -0600 (Sat, 20 Sep 2008) | 5 lines


Added a simple randomized rounding component to FP, which rounds - in the process of forming the x-tilde vector - any
integer variable with a fractional component on [0.33,0.67] to 0 or 1 with equal probability. Eventually will promote
the 0.33/0.67 parameters to FP command-line parameters, but only once I flush out issues with the new scheme.

........

r5590 | jwatson | 2008-09-22 10:13:22 -0600 (Mon, 22 Sep 2008) | 3 lines


Promoted the perturbation interval in FP to a full parameter.

........

r5595 | jwatson | 2008-09-28 21:53:12 -0600 (Sun, 28 Sep 2008) | 3 lines


Commit of Eckstein-Nediak implementation in the new heuristics framework. The version is stripped-down, in particular no parallel and only deals with binary variables. The purpose of this check-in is to assist debug of the generalized FP code. Let's not delete the bbhConfig and mipHeuristic files quite yet.

........

r5597 | jwatson | 2008-09-29 14:37:46 -0600 (Mon, 29 Sep 2008) | 3 lines


Changed the "solve()" invocation in ENHeuristic.cpp to be "resolve()"; this seems to correct much of the weird behavior in the E-N heuristic that I was observing (both here and in the generalized feasibility pump). I haven't a clue why this works, but I'll take it.

........

r5599 | jwatson | 2008-09-29 20:34:30 -0600 (Mon, 29 Sep 2008) | 3 lines


Added run-time bound on Eckstein-Nediak heuristic; parameter is ENHRunTime, quantified in seconds.

........

r5600 | jwatson | 2008-09-29 21:15:25 -0600 (Mon, 29 Sep 2008) | 3 lines


In the Eckstein-Nediak heuristic, fixed problem with imposing a limit on the aggregate number of Gomory cuts.

........

r5602 | jwatson | 2008-10-02 15:37:52 -0600 (Thu, 02 Oct 2008) | 7 lines


Various enhancements of FP output, in an effort to debug the EN heuristic. Some minor fixes to boundary-condition input cases as well.


Complete re-work of the randomization scheme for EN, including reflection of the merit function. It now behaves in that mode as a generalized FP.


I turned off any influence of the objective weight in EN for now, as it is causing *serious* issues with cycling and the like. The relative scaling with the current rule is way off, leading to this and other weird behaviors that I just killed a day mucking with.

........

r5603 | jwatson | 2008-10-03 14:37:36 -0600 (Fri, 03 Oct 2008) | 5 lines


Numerous changes/improvements related to cycle detection and escape from local minima.
Configured the current EH variant to use triangular merit functions, to verify FP-equivalence.
When verified, I will be making the merit function type a command-line argument.

........

r5604 | jwatson | 2008-10-04 16:23:25 -0600 (Sat, 04 Oct 2008) | 4 lines


Added randomized visitation sequence to FP to randomly tie-break among flip candidates whose distances are identical.
Mainly, no reason to do so deterministically, and I have seem cases where the deterministic order leads to artificial cycles.

........

r5605 | jwatson | 2008-10-04 20:08:52 -0600 (Sat, 04 Oct 2008) | 3 lines


When randomized rounding is enabled in FP, apply to construction of initial solution (i.e., the gradient) vector.

........

r5606 | jwatson | 2008-10-05 09:58:36 -0600 (Sun, 05 Oct 2008) | 3 lines


Changed FP randomization behavior to escape local minima. If fractional components are on the interval [0.33, 0.67], always flip - there isn't good gradient information anyway. In all other cases, there is less of an argument to flip, so do so with probability 0.5.

........

r5612 | jwatson | 2008-10-06 14:55:22 -0600 (Mon, 06 Oct 2008) | 3 lines


Minor updates to ENHeuristic to (1) output cumulative # of iterations for statistics gathering and (2) flip integers in the FP-fashion.

........

r5613 | jwatson | 2008-10-07 19:29:43 -0600 (Tue, 07 Oct 2008) | 3 lines


Removed LP relaxation randomization from the FP code, as the performance is rather erratic.

........

r5614 | jwatson | 2008-10-08 09:51:47 -0600 (Wed, 08 Oct 2008) | 3 lines


Added random flip option (FPRandomFlip) to test hypothesis that flipping fractionals doesn't seem to matter all that much... (?!).

........

r5616 | jwatson | 2008-10-08 13:43:10 -0600 (Wed, 08 Oct 2008) | 3 lines


Modification of perturbation mechanism in EN heuristic to break cycling.

........

r5772 | jeckstei | 2008-11-10 15:08:27 -0700 (Mon, 10 Nov 2008) | 4 lines


Changed "lp" to "sol" in various EN-Heuristic debugging statements, so
that PICO will compile with validating and MPI enabled.

........

r5779 | wehart | 2008-11-11 01:13:22 -0700 (Tue, 11 Nov 2008) | 2 lines


Update due to rename of stl_auxillary.h

........

r5789 | wehart | 2008-11-12 15:40:21 -0700 (Wed, 12 Nov 2008) | 4 lines


Adding hook for generating *.i files.


Misc fix to get the pico_test code working.

........

r5831 | wehart | 2008-11-19 19:45:31 -0700 (Wed, 19 Nov 2008) | 2 lines


Fixing bug with parallel build.

........

r5835 | wehart | 2008-11-23 20:02:08 -0700 (Sun, 23 Nov 2008) | 2 lines


Update of copyright assertions.

........

r5878 | wehart | 2008-11-24 20:16:34 -0700 (Mon, 24 Nov 2008) | 2 lines


Updates to use new GLPK release.

........

r5884 | wehart | 2008-11-25 23:39:02 -0700 (Tue, 25 Nov 2008) | 4 lines


Adding a test directory of PICO command line drivers.
This uses PyUnit? test mechanisms, leveragin the capabilites
of PyUtilib?.

........

r5885 | wehart | 2008-11-25 23:39:36 -0700 (Tue, 25 Nov 2008) | 2 lines


Adding an NL file used for testing ipconvert.

........

r5886 | wehart | 2008-11-26 13:44:49 -0700 (Wed, 26 Nov 2008) | 6 lines


Update of ipconvert to recognize integer variables.


NOTE: when converting to an MPS file, this conversion
translates the problem into a minimization problem. Many MILP solvers
cannot understand the OBJSENSE field in MPS files... :(

........

r5887 | wehart | 2008-11-26 14:51:32 -0700 (Wed, 26 Nov 2008) | 5 lines


Renaming ipconvert to pico_convert. Cindy and I realized that we
would also like to leverage glpsol and/or ampl to do conversions, but
that's not directly possible. I'm going to create an ipconvert script
in Coopr to support that more general functionality.

........

r5902 | wehart | 2008-12-01 21:14:01 -0700 (Mon, 01 Dec 2008) | 3 lines


Updates for UTILIB commit for enum labels.
Portability fix for gcc 4.3.1

........

r5920 | caphill | 2008-12-03 00:57:23 -0700 (Wed, 03 Dec 2008) | 5 lines


Outputs a solution file problem.sol.txt when we are solving only the
root LP. It prints the objective value, the (non-zero, within a hard-coded
tolerance of 1e-16) primal solution and dual solution in the form
name = value

........

r5931 | wehart | 2008-12-03 22:01:19 -0700 (Wed, 03 Dec 2008) | 2 lines


Adding diagnostic output when debug=1

........

r5932 | wehart | 2008-12-03 22:18:21 -0700 (Wed, 03 Dec 2008) | 2 lines


Printing out the optimization sense.

........

r5936 | wehart | 2008-12-04 10:59:07 -0700 (Thu, 04 Dec 2008) | 3 lines


If the IP file defines the name of the problem, then
override the default name, which is derived from the filename.

........

r5938 | wehart | 2008-12-04 15:30:27 -0700 (Thu, 04 Dec 2008) | 5 lines


When reading a *.nl file, initialize the problem name with
the filename.


Fix to milp.cpp to strip off *.nl or *.lp file suffixes.

........

r5955 | wehart | 2008-12-05 12:19:53 -0700 (Fri, 05 Dec 2008) | 2 lines


Fixes for building with validation.

........

r5981 | wehart | 2008-12-06 18:39:42 -0700 (Sat, 06 Dec 2008) | 8 lines


Setting the sense of optimization after loading the LP.


When an LP maximization problem is loaded, the sense was not being set
here.


NOTE: this may not be the right place to set the sense, but I don't
see a better place in the MILP object.

........

r5982 | wehart | 2008-12-06 19:24:45 -0700 (Sat, 06 Dec 2008) | 5 lines


Rework of pico_convert to use the utilib::OptionParser?.


Misc update of Makefile to link with tinyxml (which this
OptionParser? needs).

........

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
Line 
1/*  _________________________________________________________________________
2 *
3 *  Acro: A Common Repository for Optimizers
4 *  Copyright (c) 2008 Sandia Corporation.
5 *  This software is distributed under the BSD License.
6 *  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
7 *  the U.S. Government retains certain rights in this software.
8 *  For more information, see the README.txt file in the top Acro directory.
9 *  _________________________________________________________________________
10 */
11
12/**
13 * \file BCMip.h
14 *
15 * Classes for a branch-and-cut search: BCMip, BCMipProblem, BCMipNode
16 */
17
18#ifndef pico_BCMip_h
19#define pico_BCMip_h
20
21#include <acro_config.h>
22#include <utilib/HashedSet.h>
23#include <utilib/_math.h>
24#include <utilib/seconds.h>
25#include <pico/cutFinder.h>
26#include <pico/milpNode.h>
27#include <pico/milp.h>
28#include <pico/MILPProblem.h>
29#include <pico/BCMipBase.h>
30#include <pico/picoRowCut.h>
31#include <pico/PicoLPInterface.h>
32#include <pico/lpSetupObjCuts.h>
33
34namespace pico {
35
36  class SNLKnapsack;
37
38using utilib::HashedSet;
39
40// An index into the set of cut finders.  For now there's
41// always a pool scanner and it's always first.
42#define POOLSCANNERTYPE 0
43
44class cutFinder;
45
46// Do we need to define debug as a global parameter here?
47
48class BCMip;
49class BCMipProblem;
50class BCMipNode;
51
52// Used to manage the cutting schedule
53
54class cutQObject
55{
56friend class BCMip;
57protected:
58// If this ever changes, be sure to propagate effects to parBCMip::(un)packCutQueue
59  BasicArray<picoRowCut *>new_cuts;
60  double time_to_find;
61  cutFinder *finder_used;
62  double finder_readiness;
63  bool found_infeasibility;
64  bool first_run_this_finder_this_SP;
65public:
66
67  // This is an int rather than a size_type because it's passed as an argument
68  // (by reference) to generateCuts(...).  We want that method, which could be
69  // used by other PICO users, to have more generic types whenever possible.
70  int num_cuts;
71  cutQObject()
72    {
73    time_to_find = 0.0;
74    finder_used = NULL;
75    finder_readiness = 0.0;
76    found_infeasibility = false;
77    first_run_this_finder_this_SP = false;
78    };
79
80BasicArray<picoRowCut *>& getCutContainer()
81  {return(new_cuts);};
82
83double getTime()
84  {return(time_to_find);};
85
86 double getAverageTime()
87   {return(time_to_find/num_cuts);};
88
89cutFinder *getFinderUsed()
90  {return(finder_used);};
91
92double getFinderReadiness()
93  {return(finder_readiness);};
94
95bool getInfeasibility()
96  {return(found_infeasibility);};
97
98bool getFirstRun()
99  {return(first_run_this_finder_this_SP);}
100
101void setCutQObject(double new_time, cutFinder *new_finder,
102                   double new_readiness, bool new_infeas,
103                   bool firstSolveForSP)
104  {
105  time_to_find = new_time;
106  finder_used = new_finder;
107  finder_readiness = new_readiness;
108  found_infeasibility = new_infeas;
109  first_run_this_finder_this_SP = firstSolveForSP;
110  }
111
112 void removeCut(size_type to_remove)
113   {
114   delete new_cuts[to_remove];   
115   new_cuts[to_remove] = new_cuts[--num_cuts];
116   }
117
118 void removeCut(picoRowCut *to_remove)
119   {
120     for (int i = 0; i < num_cuts; i++)
121       {
122         if (new_cuts[i] == to_remove)
123           {
124           removeCut(i);
125           return;
126           }
127       }
128     EXCEPTION_MNGR(std::runtime_error, "Trying to remove a cut from a queue that doesn't contain it\n");
129   }
130
131 // Like removeCut, but no deletion
132
133 void releaseCut(size_type to_remove)
134   {
135   new_cuts[to_remove] = new_cuts[--num_cuts];
136   }
137
138 void releaseCut(picoRowCut *to_remove)
139   {
140     for (int i = 0; i < num_cuts; i++)
141       {
142         if (new_cuts[i] == to_remove)
143           {
144           releaseCut(i);
145           return;
146           }
147       }
148     EXCEPTION_MNGR(std::runtime_error, "Trying to release a cut from a queue that doesn't contain it\n");
149   }
150
151
152 // Replace cut, don't delete the old one
153
154 void replaceCut(size_type to_replace, picoRowCut *new_cut)
155   {
156   new_cuts[to_replace] = new_cut;
157   }
158
159 void replaceCut(picoRowCut *to_replace, picoRowCut *new_cut)
160   {
161     for (int i = 0; i < num_cuts; i++)
162       {
163         if (new_cuts[i] == to_replace)
164           {
165           replaceCut(i, new_cut);
166           return;
167           }
168       }
169     EXCEPTION_MNGR(std::runtime_error, "Trying to replace a cut in a queue that doesn't contain it\n");
170   }
171
172
173void printCuts()
174  {
175  for (int i=0; i < num_cuts; i++)
176    (new_cuts[i])->print();
177  }
178
179};
180
181
182class BCMipProblem: public MILPProblem, public BCMipBase
183{
184 public:
185  // If this returns true, when a BCMipNode is made active, the LP
186  // solver state is restored to exactly what it was when the current
187  // LP solution was created.  For example, CGL cut generators get
188  // their problem information this way.  The knapsack cut generator
189  // at least uses the matrix directly.  If there is a derived class
190  // that uses only cut generators that do not use the
191  // OsiSolverInterface for problem data (e.g. applications within
192  // PICO that use the BCMipProblem-derived class), then setting this
193  // to false, means we can leave old (globally valid) cuts in the LP
194  // if we want.  Probably means a more efficient swap.
195
196  // Can this be static in derived classes?
197   virtual bool resetNodeState()
198     {return(true);};
199
200  //// Make an empty MIP, for parallel initialization
201  BCMipProblem()
202    {};
203};
204
205class BCMip: virtual public MILP, virtual public BCMipBase, public BCMipParams
206{
207
208  // For tracking performance
209public:
210  DoubleVector BCMIPTimingInfo;
211
212protected:
213
214  // The SNL knapsack cut finder.  This will probably move into a finder manager once that's
215  // implemented
216
217  SNLKnapsack *SNLKnapsack_finder;
218
219  // Gives cuts serial numbers similar to subproblem serial numbers.  Mostly
220  // used for debugging.
221  int cutCounter;
222
223  // TODO: make sure this is reasonably efficient.  May be necessary to use
224  // picoRowCut here rather than picoRowCut * for correct operation of the hash table.
225
226  // Cuts currently loaded into the LP
227  HashedSet<picoRowCut,false> loadedCuts;
228
229  // Globally valid, possibly useful cuts *not* loaded
230
231  HashedSet<picoRowCut,false> cutPool;
232
233  // Maximum number of cuts loaded at any given time (unless all cuts are
234  // binding)
235
236  int loadedPoolMaxSize;
237
238  // Usually generate this many cuts before incorporating them
239
240  int cutBlockSize;
241
242  // Number of subproblems for which all cut Finders have been run.
243  // This is used to control the cutting process.  Early on, when
244  // BCMip hasn't learned anything about the power of the various
245  // cutFinders on this problem, run all finders once on a subproblem
246  // before starting the Finder competition.
247
248  int numFullFindersBounds;
249
250        ///
251  // These can be applied optionally
252  BasicArray<cutFinder*> cutGenerators;
253  // These must be applied continuously until they return no more cuts
254  // before we can consider the LP relaxation complete.
255  BasicArray<cutFinder*> separationAlgorithms;
256
257  void BCMIPInit();
258
259  BasicArray<cutQObject> cutQueue;
260  int num_queue_objects;
261  int queue_total;
262
263  // Quality-tracking information associated with each cut
264  // finder in cutGenerators
265
266  BasicArray<finderQualityObject> cutQuality;
267
268  // Used for scheduling cut generators
269
270  DoubleVector finderSchedBias;
271  DoubleVector lastRunTime;
272 
273  // Estimate of the amount of time needed to solve the
274  // first LP for a subproblem.
275  double sp_first_solve_time_est;
276
277  // If bit i is on, then disable the cutFinder of type i during set
278  // up in preprocessing.  There may be a better way to do this, but
279  // we can't prevent the cutter anouncing its availability (that's
280  // done in static initialization). That would require editing the
281  // cut finder class and we'd like the users to be able to disable
282  // with a call in the driver.  We can't disable after the array of
283  // cut generators is made, because that's done after the search
284  // starts (and that's a monolithic process in the driver).  So we
285  // set a bit before the search starts and don't enter this cut
286  // finder into the array of cut finders.  We'll use a pointer so we
287  // can ditch this after it's used.
288  BitArray *disableBuiltInCuts;
289
290public:
291
292  ///
293  BCMipProblem *BCMIP;
294
295  // If we know that we have stopped the computation early
296  // (gracefully, no abort), then we will need to make the
297  // final output reflect that.  Return true if there is a
298  // reason for an early termination and set the abortReason
299  // char * to say why.  Called when the computation is over.
300  // We override the MILP class method, since BCMip controls
301  // the onlyRootHeuristic parameter.  This is a reason for
302  // stopping early.
303
304  virtual bool knownStopEarly()
305    {
306      if (!onlyRootHeuristic)
307        return (false);
308      if (!abortReason)
309        abortReason = "Only doing Root Heuristic";
310    }
311
312  // Users can inform the system of their own domain-specific cut
313  // finders.  If isSeparator is true, then this is treated as a
314  // separation algorithm.  It is applied repeatedly after every LP
315  // solution until it returns no more cuts.  If isSeparator is false,
316  // then this is treated as a cut finder that can be applied
317  // optionally.  This method returns the cut type, which the user
318  // should store locally and return in the cut finder's cutType()
319  // method.  Default behavior is to turn off all built-in (generic)
320  // cut finders if users provide a new (optional) cut finder, but to
321  // leave them on for separation algorithms.  If they users want to
322  // change this, they will have to turn on/off cutfinders using.
323
324  int declareCutFinder(cutFinder *newFinder, bool isSeparator);
325
326  bool separationAlgorithmsExist()
327    {
328      return (separationAlgorithms.size() > 0);
329    }
330
331  // Apply separation algorithms.  If one has a cut to add, this adds it
332  // and returns true.  If all separation algorithms return no cuts, return false.
333
334  bool addSepCut(BCMipNode *curr_node);
335
336
337  // Highest quality (weighted by readiness) among all finders just
338  // after the last resolve
339
340  double bestFinderQuality;
341
342  // finders with effective quality (i.e. quality weighted by
343  // readiness) at least this high can compete to run during a
344  // competitive phase.
345  double finderThreshold;
346
347  void updateFinderQualityThreshold(BCMipNode *curr_node);
348
349  virtual bool finderCompetitive(size_type finderID, BCMipNode *curr_node)
350    {
351      return(cutGenerators[finderID]->readiness(curr_node) * cutQuality[finderID].AvgQuality >= finderThreshold);
352    }
353
354  ///
355  bool setup(int& argc, char**& argv, BCMipProblem& _MIP)
356        {
357          setBCMipProblem(_MIP);
358          return MILP::setup(argc,argv,_MIP);
359        }
360
361  // If a cut finder runs, it must run for a positive time.  This is
362  // the smallest unit of time the timer can distinguish from zero.
363  double min_runtime;
364
365  // Used for debugging.  Print information used to schedule cut finders and branching
366
367  void printCutSchedulingInfo(BCMipNode *curr_node);
368
369  // Used for debugging and tuning.
370
371  void printCutFinderQualStats(bool fullInfo = false)
372    {
373      size_type numFinders = cutGenerators.size();
374      ucout << "Current cut finder quality statistics:\n";
375        for (size_type i = 0; i < numFinders; i++)
376          {
377            if (cutGenerators[i] == NULL)  continue;
378            ucout << cutGenerators[i]->cutName() << ": ";
379            cutQuality[i].printFinderQualInfo(fullInfo);
380          }
381    }
382
383  // Number of rows in the original (core) problem.  Inactive subproblems only
384  // store slack variables for these rows in the basis.
385
386  int numCoreRows;
387
388  // For debugging.  Indices of solutions in the solutionsToWatch array (from milp)
389  // contained by the currently active problem, and the number of such solutions.
390  // Used to check for cuts that cut off known feasible (optimal) solutions.
391
392  IntVector containedWatchedPointsInd;
393  int numContainedWatchedPoints;
394
395  // This used to set the MIP pointer too, but that's called elsewhere
396  // (through the setup methods
397
398  void setBCMipProblem(BCMipProblem& inputMIP)
399        {
400        BCMIP = &inputMIP;
401        }
402
403  // Note:  Does NOT assume ownership of inputMIP
404  // This no longer attempts to initialize the hash table names, since
405  // currently the hash tables in utilib ignore this name.  If that changes,
406  // initialize it with a direct hash-table method from here.
407  void MILPInit(bool initVB=true)
408    {
409    if (initVB)
410       MILP::MILPInit();
411    BCMIPInit();
412    }
413
414  /* This has gone away at the MILP level at least for now
415  //For parallel versions that have I/O capabilities
416  BCMip(int argc,char** argv):
417    MILP(argc, argv), loadedCuts("LoadedCuts"), cutPool("CutPool")
418    {
419    MCMIPInit();
420    }; */
421
422  BCMip();
423
424  /// Destructor
425  ~BCMip();
426
427  void reset(bool resetVB = true);
428
429  virtual pebbl::branchSub* blankSub();
430
431  int getCutSerial()
432    {return(cutCounter++);}
433
434  int poolSize()
435    {
436    return (cutPool.size());
437    }
438
439  // For debugging.
440
441  void printCutDataStructures();
442
443  // Update the cut pools to reflect this new set of active cuts
444  // (and all previous ones deleted) At least for now, this is called
445// when a subproblem is made active.  Each of these new cuts should
446// have their reference count decremented whether they're moving from the
447// cut pool or staying in the loaded pool (previous problem incremented this
448// count before deactivating.  The two sets of cuts are persistent and regular.
449
450  void replaceLoadedCuts(BasicArray<picoRowCut *> &new_cuts1,BasicArray<picoRowCut *> &new_cuts2);
451
452  // Update the loaded cut pool to reflect addition of these cuts
453  // Also set the loaded flag in each of these cuts.
454
455  void addLoadedCuts(BasicArray<picoRowCut *> &new_cuts, int num_cuts);
456
457  // Adds cuts to the LP and updates the pools.   If adding these cuts puts us over
458  // the loaded-cut limit, remove nonbinding cuts in order of current scaled dual.
459
460  void addCuts(BasicArray<picoRowCut *> &new_cuts, int num_cuts);
461
462  // Remove a cut from the loaded-cuts pool.  Move it to the cut
463  // pool if it's poolable, and otherwise, if its reference count is
464  // zero, delete it
465
466  void unloadCut(picoRowCut *to_unload);
467
468  // Add a cut to the pool.  If there is already a cut in the pool
469  // parallel to the new one, keep only one in the pool, make sure
470  // both have the strongest upper bound (rhs), and, if either has
471  // a zero reference count, delete it.
472
473  void addToCutPool(picoRowCut *to_add);
474
475  // As above for a set of cuts
476  void unloadCuts(BasicArray<picoRowCut *> &to_unload);
477
478
479// Put all the binding cuts into cutlist.  Put the index (into cutlist)
480// of all those with degenerate basic slack variables into degenerate_list.
481// Binding cuts have zero slack.  BasicArrays
482// are sized exactly afterward (assumed, eg, by the PICO cut
483// management piece in the LP interface).  We no longer assume the basis
484// in the LP interface is valid (store degeneracy info with the cuts
485// immediately after the solve).  Pseudocost initialization, eg, when
486// determining branching readiness, made basis invalid sometimes.
487// Also record this basis information (slackness) for persistent cuts.
488
489
490  void getBindingRows(BasicArray<picoRowCut*>& persistent_cuts,
491                      BasicArray<int>& slack_persistent_list,
492                      BasicArray<picoRowCut*>& cutlist,
493                      BasicArray<int>& degenerate_list);
494
495  // Called after a resolve during cutting.  Age nonbinding cuts.
496  // Remove if they're old enough.  Return true if cuts were removed
497  // and false otherwise.
498
499  bool ageCuts();
500
501  // Add the cuts in the queue into the LP and loaded pool.  Return
502  // the number added. If multiple finders have found the same cut,
503  // removes the redundancy and gives credit to the finder that found
504  // the tightest cut (if the bounds aren't the same).  For ties,
505  // credit goes to the fastest finder.  This assumes that the only
506  // finder that can return cuts with a positive reference count is the
507  // cut pool scanner.  That is, we assume all
508// other finders are making new cuts.  If a finder generates a cut redundant with the cut pool and
509// wins (stronger or faster), then we use the strongest cut, and the correct finder gets credit.
510
511  int incorporateCuts();
512
513  // Removes cut redundancy from a single cut finder (certainly an issue with CGL
514  // finders, or any finder that doesn't police itself).  This is separate from the
515  // redundancy elimination from BCMip::incorporateCuts because there's no need to
516  // assign credit.  Just keep the best.  Run this as we go along to avoid keeping a
517  // lot of redundant cuts around (a minor efficiency improvement).
518
519void removeOneFinderCutRedundancy(int queue_index);
520
521// Clean up the cut management objects if a cut finder determines a subproblem is infeasible
522
523void cleanUpAfterInfeasibility(double SPBound, bool newSubproblem);
524
525void trimLoadedPool(int num_to_remove);
526
527  // Estimate of the magnitude of the (optimal) objective function
528
529  double solMagnitudeEst();
530
531  // Takes an LP solution x.  Scans the cut pool.  Removes all cuts
532  // that are violated by x and returns them.
533
534  void scanCutPool(BasicArray<picoRowCut *>& violated_cuts, int &num_cuts, DoubleVector &x);
535
536  // Should we run all cut Finders on this subproblem?  Currently yes if this is one of the first
537  // suproblems we've bound.  If this method changes, consider effect on postRunAllFinders().
538
539  bool shouldRunAllFinders()
540    {return(numFullFindersBounds < tryAllFindersLimit);}
541
542  // Keep track of the number of times we've run all finders (early), so we'll know when to stop.
543  // If this is the last such round, set cut qualities for the competition-only phase.  If the
544  // pass might have terminated early because of infeasibility detection, don't count the pass.
545  void postRunAllFinders(bool maybeNotAllRun);
546
547  // Run all finders that have positive readiness
548
549  void runAllFinders(BCMipNode *curr_problem,
550                     BitArray &finders_run,
551                     bool &problemInfeasible);
552
553  // Assumes the cutFinder passed in has positive readiness
554  void runFinder(cutFinder *curr_finder, BCMipNode *curr_problem, bool &problemInfeasible);
555
556  // Select a finder and run it
557  cutFinder  *runOneFinder(BCMipNode *curr_problem, bool &problemInfeasible);
558
559    // Return true if and only if any of the cut finders' readiness
560    // measure beats the branching readiness
561  bool readinessCheck(BCMipNode *curr_problem);
562
563  bool canStillCut(BCMipNode *curr_problem,
564                   PicoLPInterface::problemState &this_state);
565
566  cutFinder *getCutFinder(int finderType)
567        {return(cutGenerators[finderType]);};
568
569  int numCutFinders()
570        { return((int)cutGenerators.size()); }
571
572  // For each set of cuts in the cut queue, update the
573  // finder quality measures.
574
575  void updateCutFinderQuality(DoubleVector &x, double time_to_add, double boundMovement,
576                              double branchingThreshold);
577
578  void recomputeCutFinderQualityAverages(finderQualityObject *qualObject, double new_qual,
579                                         double newTime, double branchingThreshold, bool newSubproblem);
580
581  // Update quality measures after a cut finder has determined a subproblem is infeasible
582
583  void updateCutFinderQualityFromInfeasible(finderQualityObject *qualObject, double this_time,
584                                            double SPBound, bool newSubproblem);
585
586  // Upon starting a new subproblem (never cut before), set the cut
587  // quality the scheduler.  Usually this is the average quality for
588  // this finder on new subproblems, but if that's very low compared
589  // to the branching quality, add a little so finders still have a
590  // chance to run occasionally.
591
592  void newSPSetCutQual(double branchQual);
593
594  // Decrease all cut schedule bias values by the min, to insure numerical
595  // stability (e.g. no wrap)
596  void resetSchedBias();
597
598  virtual void preprocess();
599
600  double spFirstSolveTimeEst()
601    {return(sp_first_solve_time_est);};
602
603  void updateSPFirstSolveTimeEst(double new_time, bool initializing = false);
604
605  // All cut queue objects should have no cuts (handled before this call)
606  void resetCutQueue()
607    {
608    num_queue_objects = 0;
609    queue_total = 0;
610    };
611
612  void removeCutQObject(size_type to_remove);
613
614  bool cutQueueBig()
615    {
616    if (queue_total >= cutBlockSize)
617      return(true);
618    return(false);
619    }
620  virtual void printAllStatistics(std::ostream& stream = std::cout);
621
622// This should be called from the driver before search starts.  Indicate a
623// built-in cut finder shouldn't be added to the list of cut finders in
624// BCMip preprocessing.
625
626void disableBuiltInCutFinder(cutFinder::cutTypeID thisFinder)
627  {disableBuiltInCuts->set(thisFinder);};
628
629  void enumCheck() { };   // Don't barf in reset if enumerating
630
631};
632
633
634// Note: don't have node classes share parameter classes (like BCMipBase)
635// with branching classes and problem classes (parameters would need to be
636// reset on construction of each node).
637class BCMipNode: virtual public MILPNode, virtual public BCMipBase
638{
639  friend class lpSetupObjCuts;
640 
641protected:
642
643  bool separationAlgsExist;
644
645  // Mostly for safety
646  bool active;
647
648  DoubleVector full_solution;
649
650  // These are passed into routines in the LP interface class (set from there)
651
652  // List of binding cuts (needed, eg, to get this LP solution)
653  BasicArray<picoRowCut *> binding;
654  // The indices (into binding) of the cuts with basic slack variables
655  IntVector basic_slack_cuts;
656
657  // Cuts explicitly enforced as long as the node is alive.  For example, these
658  // are used for branching on constraints.
659
660  BasicArray<picoRowCut *> persistent_cuts;
661
662  // The indices (into persistent_cuts) of the persistent cuts with
663  // basic slack variables.  Since these are always enforced, they can
664  // go slack and will therefore have basic slacks more often.  Thus
665  // this will be a higher percentage of the persistent cuts than the
666  // percentage of binding in basic_slack_cuts, but there probably
667  // won't be a lot of persistent cuts.  If that changes (there are a
668  // lot of persistent cuts), consider representing this as a
669  // BasisArray.
670
671  IntVector persistent_basic_slacks;
672
673  // Bit i is 1 iff cut finder i has been invoked since the last
674  // time this LP was solved
675
676  BitArray finderHistory;
677
678  // Bit i is 1 iff cut finder i has been invoked at least once on
679  // this subproblem
680
681  BitArray finderSPHistory;
682
683  int NumLPSolves;
684 
685  // The last cut finder used for this subproblem
686  cutFinder *lastFinder;
687
688  // A measure of the quality of branching
689  double branchQuality;
690
691  // How ready are we to branch? (between 0 and 1)
692
693  double branchReadiness;
694
695public:
696
697  // Is this the first time this finder has been run on this
698  // subproblem?  Because this bit is set *after* the run and
699  // corresponding quality updates, this will be 0 if and only if the
700  // finder was not run before on this subproblem.
701
702  bool firstRun(int finderType)
703    {
704      return(finderSPHistory(finderType)==0);
705    }
706
707  bool useSavedSolutionNow()
708    {
709      return(MILPNode::useSavedSolutionNow() && NumLPSolves == 0);
710    }
711
712// Add these cuts to the array of persistent cuts and mark
713// them as persistent.  We assume all cuts are either persistent or
714// regular in all subproblems to which they apply.
715// We'll likely only add one at a time, but we'll use this more
716// general form for now.
717
718  void addPersistentCuts(BasicArray<picoRowCut *> newCuts);
719
720  // The current LP relaxation solution
721  DoubleVector &x(){return(full_solution);};
722
723  virtual int initialBasisSize();
724
725  void makeActiveLP();
726
727  // If we're checking for cutting off watched points, this records which
728  // watched points are contained in this subproblem
729
730  void recordWatchedPoints();
731
732  void noLongerActiveLP();
733
734  PicoLPInterface::problemState boundComputationGuts();
735
736  virtual void beforeCutsHeuristics();
737
738  // Enumeration stuff, in some cases overriding the standard MILPNode
739  // versions that just throw exceptions.
740
741  int terminalSplit();
742
743  void terminalChildSetup(int whichChild,BCMipNode* parent);
744
745  IntVector tSplitVars;  // Extra node data used for splitting terminal
746                         // nodes.  Will usually just have zero size
747                         // and not bother anybody.
748
749  // Other overrides
750 
751  void setBasis()
752    {
753      if (!basisValid)
754        return;
755      // For activating problems.  Cuts must already be in the LP.
756      if (!active && (binding.size() > 0 || persistent_cuts.size() > 0))
757        {
758        lp->setBasisWithCuts(persistent_cuts,
759                           persistent_basic_slacks,
760                           binding,
761                           basic_slack_cuts,
762                           basis);
763        }
764    // For active problems with full basis stored locally
765    else
766      MILPNode::setBasis();
767    }
768
769  void setBounds()
770    {
771    if (!IwasLast())
772      MILPNode::setBounds();
773    }
774   
775  void getSolution()
776    {
777    MILPNode::getSolution(); // for integer vars
778    lp->getColSolution(full_solution.data()); 
779    }
780
781  void getBasis(BasisArray& basisBuffer)
782    {
783      // TODO: This could be inefficient if it causes a lot of
784      // allocation/deallocation.  but we have to make sure the size()
785      // function is correct (used heavily in the PicoLPInterface).
786      // Check that at least resizes that don't require a new
787      // allocation (e.g. ones that move only within the last
788      // allocated byte) don't do allocation.
789      size_type basisSize = lp->getNumCols() + lp->getNumRows();
790      if (basisSize > basisBuffer.size())
791        basisBuffer.resize(basisSize);
792      lp->getBasis(basisBuffer);
793      basisValid = lp->basisAvailable();
794    }
795
796  // Pseudocost initialization uses the LP solver.  Make sure the LP
797  // is back in the state it was after solving the node initially.
798  // Otherwise, the cut generators could see the wrong point
799  // (solution).  Also, if this is the root, use average pseudocost
800  // resolve time to set the initial subproblem solve time estimate.
801
802  virtual void postPseudoCostInit(IntVector &initVars,
803                                  bool dontRemoveAllBranchChoices,
804                                  int &numBranches);
805
806  PicoLPInterface::problemState resolveLP(bool &continue_cutting, double &resolve_time);
807
808  // Override pseudocost initialization solve call to record timings for
809  // initialization of subproblem-first-time-solution time at the root
810  PicoLPInterface::problemState PCInitSolve();
811
812  // Add cuts from the cut queues to the current problem and resolve the LP.
813
814  PicoLPInterface::problemState incorporateCuts(bool &continue_cutting);
815
816
817  // Used by cut finders to determine readiness
818
819  // True iff the LP has been resolved since the last
820  // time this finder was invoked for this subproblem.
821  // Vacuously true if this finder has never been invoked.
822
823  bool LPResolved(int finderType)
824    {return(finderHistory(finderType)==0);};
825
826  // True iff another cut finder has been called since
827  // the last use of this cut finder for this subproblem
828  bool otherFinderCalled(int finderType);
829
830  // True iff this is a new subproblem for which no cutting has taken
831  // place yet.  (other than the addition of cuts from the pool, for example).
832
833  bool newSubproblem()
834    {return (NumLPSolves <= 1);};
835
836  // True if this problem has had at most one set of cuts
837  // incorporated.  We determine using number of LP solves.  The first
838  // is the bounding after branching.  The second is after the first
839  // cut incorporation.  This case is special since we track cut
840  // finder quality based on how well cuts work on a "new" subproblem
841
842  bool firstCutSet()
843    {return(NumLPSolves <= 2);}
844
845  // Record the pseudocost effects from a branch only on the LP solve
846  // immediately after the branch (not on subsequent resolves for cut
847  // incorporation)
848  bool updatePCFromBoundingSolve()
849    {return(NumLPSolves == 0 && MILPNode::updatePCFromBoundingSolve());};
850
851  // Controls whether bounds all bounds are loaded into the LP, or
852  // just the few that have changed from the parent.
853
854  bool isWarmStart();
855
856  // Controls whether all the cuts in the LP are replaced with those
857  // from the subproblem, or are just held over because the the last
858  // LP was for the same subproblem or its parent.  Typically this is
859  // when either this node or its parent was the last to use the LP.
860  // But inheritance from parents is turned off for special
861  // enumeration children, because they can be weird.
862
863  bool replaceCutsToActivate();
864
865  // MILPNode required the state be bounded, but we will want to
866  // recognize this as a possible reason for stopping cutting (before
867  // we declare the state bounded).
868  bool candidateSolution()
869  {
870    if (nFrac == 0)
871      return true;
872    return false;
873  }
874
875
876  // Update the measure of the benefit of branching at this point.
877  // This can change anytime the LP is solved.
878
879  bool computeBranchQuality();
880
881  // Update the measure of readiness to branch.  Should probably be
882  // called with the above update of branching quality
883
884  void computeBranchReadiness();
885
886  double branchingReadyMeasure();
887
888  virtual branchSub *makeChild(int whichChild);
889  void childInit(BCMipNode *child, int whichChild, bool initVB = true);
890
891  virtual BCMip *bcGlobal() {return globalPtr;};
892  MILP *mGlobal() {return (MILP *)globalPtr;};
893  pebbl::branching *bGlobal() const {return (pebbl::branching *)globalPtr;};
894
895  void BCMipNodeFromBCMip(BCMip* master, bool initVB=true);  // For making the root
896  void BCMipNodeAsChildOf(BCMipNode* parent, bool initVB=true); // For making a child
897  // For use with derived classes that will call constructors for MILPNode and branchSub
898  BCMipNode()
899    { }
900
901  ~BCMipNode();
902
903  void BCMipNodeInit();
904
905  virtual void printMILPNode();
906
907  // Stuff for setting up heuristic LP's.  Override the similar methods
908  // in MILPNode.
909
910  // Make an object that will set up an LP later
911
912  lpSetupObj* makeLPSetupObj()
913    {
914      return new lpSetupObjCuts(this);
915    };
916
917  // Set up an LP for the Eckstein-Nediak heuristic code in .cpp file)
918 
919  void setupLP(__LPSolver_Interface* sol);
920
921  // End of heuristic LP setup stuff.
922
923 protected:
924  BCMip *globalPtr;
925};
926
927
928
929// This is used by the handlers (search engines) to create the root.
930// As an inline, it will generate a call to BCMipNode with the BCMip (branching
931// class) pointer as an argument.
932
933inline pebbl::branchSub* BCMip::blankSub()
934{
935MEMDEBUG_START_NEW("BCMipNode")
936BCMipNode *tmp = new BCMipNode;
937#ifdef ACRO_VALIDATING
938 if (tmp == NULL)
939   EXCEPTION_MNGR(std::runtime_error, "Failed to allocate a new BCMipNode\n");
940#endif
941tmp->BCMipNodeFromBCMip((BCMip *)this);
942MEMDEBUG_END_NEW("BCMipNode")
943return (pebbl::branchSub *)tmp;
944};
945
946} // namespace pico
947
948#endif
Note: See TracBrowser for help on using the repository browser.