Scippy

SCIP

Solving Constraint Integer Programs

sepa_mcf.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2017 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not email to scip@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file sepa_mcf.c
17  * @brief multi-commodity-flow network cut separator
18  * @author Tobias Achterberg
19  * @author Christian Raack
20  *
21  * We try to identify a multi-commodity flow structure in the LP relaxation of the
22  * following type:
23  *
24  * (1) sum_{a in delta^+(v)} f_a^k - sum_{a in delta^-(v)} f_a^k <= -d_v^k for all v in V and k in K
25  * (2) sum_{k in K} f_a^k - c_a x_a <= 0 for all a in A
26  *
27  * Constraints (1) are flow conservation constraints, which say that for each commodity k and node v the
28  * outflow (delta^+(v)) minus the inflow (delta^-(v)) of a node v must not exceed the negative of the demand of
29  * node v in commodity k. To say it the other way around, inflow minus outflow must be at least equal to the demand.
30  * Constraints (2) are the arc capacity constraints, which say that the sum of all flow over an arc a must not
31  * exceed its capacity c_a x_a, with x being a binary or integer variable.
32  * c_a x_a does not need to be a single product of a capacity and an integer variable; we also accept general scalar
33  * products.
34  */
35 
36 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
37 
38 
39 /* #define COUNTNETWORKVARIABLETYPES */
40 /* #define MCF_DEBUG */
41 
42 /* algorithmic defines in testing phase*/
43 /* #define USEFLOWFORTIEBREAKING */
44 /* #define USECAPACITYFORTIEBREAKING */
45 /* #define TIEBREAKING */
46 #define BETTERWEIGHTFORDEMANDNODES
47 
48 #include <assert.h>
49 #include <string.h>
50 
51 #include "scip/sepa_mcf.h"
52 #include "scip/cons_knapsack.h"
53 #include "scip/pub_misc.h"
54 
55 
56 #define SEPA_NAME "mcf"
57 #define SEPA_DESC "multi-commodity-flow network cut separator"
58 #define SEPA_PRIORITY -10000
59 #define SEPA_FREQ 0
60 #define SEPA_MAXBOUNDDIST 0.0
61 #define SEPA_USESSUBSCIP FALSE /**< does the separator use a secondary SCIP instance? */
62 #define SEPA_DELAY FALSE /**< should separation method be delayed, if other separators found cuts? */
63 
64 /* changeable parameters*/
65 #define DEFAULT_NCLUSTERS 5 /**< number of clusters to generate in the shrunken network */
66 #define DEFAULT_MAXWEIGHTRANGE 1e+06 /**< maximal valid range max(|weights|)/min(|weights|) of row weights for CMIR */
67 #define DEFAULT_MAXTESTDELTA 20 /**< maximal number of different deltas to try (-1: unlimited) for CMIR */
68 #define DEFAULT_TRYNEGSCALING FALSE /**< should negative values also be tested in scaling? for CMIR */
69 #define DEFAULT_FIXINTEGRALRHS TRUE /**< should an additional variable be complemented if f0 = 0? for CMIR */
70 #define DEFAULT_DYNAMICCUTS TRUE /**< should generated cuts be removed from the LP if they are no longer tight? */
71 #define DEFAULT_MODELTYPE 0 /**< model type of network (0: auto, 1:directed, 2:undirected) */
72 #define DEFAULT_MAXSEPACUTS 100 /**< maximal number of cuts separated per separation round (-1: unlimited) */
73 #define DEFAULT_MAXSEPACUTSROOT 200 /**< maximal number of cuts separated per separation round in root node (-1: unlimited) */
74 #define DEFAULT_MAXINCONSISTENCYRATIO 0.02 /**< maximum inconsistency ratio (inconsistencies/(arcs*commodities)) at all */
75 #define DEFAULT_MAXARCINCONSISTENCYRATIO 0.5 /**< maximum inconsistency ratio for arcs not to be deleted */
76 #define DEFAULT_CHECKCUTSHORECONNECTIVITY TRUE /**< should we only separate if the cuts shores are connected */
77 #define DEFAULT_SEPARATESINGLENODECUTS TRUE /**< should we separate inequalities based on single node cuts */
78 #define DEFAULT_SEPARATEFLOWCUTSET TRUE /**< should we separate flowcutset inequalities */
79 #define DEFAULT_SEPARATEKNAPSACK TRUE /**< should we separate knapsack cover inequalities */
80 
81 /* fixed parameters */
82 #define BOUNDSWITCH 0.5
83 #define USEVBDS TRUE
84 #define ALLOWLOCAL TRUE
85 #define MINFRAC 0.05
86 #define MAXFRAC 0.999
87 
88 #define MAXCOLS 2000000 /**< maximum number of columns */
89 #define MAXAGGRLEN(nvars) (0.1*(nvars)+1000) /**< maximal length of base inequality */
90 #define MINCOLROWRATIO 0.01 /**< minimum ratio cols/rows for the separator to be switched on */
91 #define MAXCOLROWRATIO 100.0 /**< maximum ratio cols/rows for the separator to be switched on */
92 #define MAXFLOWVARFLOWROWRATIO 100.0 /**< maximum ratio flowvars/flowrows for the separator to be switched on */
93 #define MAXARCNODERATIO 100.0 /**< maximum ratio arcs/nodes for the separator to be switched on */
94 #define MAXNETWORKS 4 /**< maximum number of networks to consider */
95 #define MAXFLOWCANDDENSITY 0.1 /**< maximum density of rows to be accepted as flow candidates*/
96 #define MINCOMNODESFRACTION 0.5 /**< minimal size of commodity relative to largest commodity to keep it in the network */
97 #define MINNODES 3 /**< minimal number of nodes in network to keep it for separation */
98 #define MINARCS 3 /**< minimal number of arcs in network to keep it for separation */
99 #define MAXCAPACITYSLACK 0.1 /**< maximal slack of weighted capacity constraints to use in aggregation */
100 #define UNCAPACITATEDARCSTRESHOLD 0.8 /**< threshold for the percentage of commodities an uncapacitated arc should appear in */
101 #define HASHSIZE_NODEPAIRS 500 /**< minimal size of hash table for nodepairs */
102 
103 /* #define OUTPUTGRAPH should a .gml graph of the network be generated for debugging purposes? */
104 
105 
106 #ifdef SCIP_DEBUG
107 #ifndef MCF_DEBUG
108 #define MCF_DEBUG
109 #endif
110 #endif
111 
112 #ifdef MCF_DEBUG
113 #define MCFdebugMessage printf
114 #else
115 /*lint --e{506}*/
116 /*lint --e{681}*/
117 #define MCFdebugMessage while(FALSE) printf
118 #endif
119 
120 /*
121  * Data structures
122  */
123 
124 /** model type of the network */
126 {
127  SCIP_MCFMODELTYPE_AUTO = 0, /**< model type should be detected automatically */
128  SCIP_MCFMODELTYPE_DIRECTED = 1, /**< directed network where each arc has its own capacity */
129  SCIP_MCFMODELTYPE_UNDIRECTED = 2 /**< directed network where anti-parallel arcs share the capacity */
130 };
132 
133 /** effort level for separation */
135 {
136  MCFEFFORTLEVEL_OFF = 0, /**< no separation */
137  MCFEFFORTLEVEL_DEFAULT = 1, /**< conservative separation */
138  MCFEFFORTLEVEL_AGGRESSIVE = 2 /**< aggressive separation */
139 };
141 
142 /** extracted multi-commodity-flow network */
144 {
145  SCIP_ROW*** nodeflowrows; /**< nodeflowrows[v][k]: flow conservation constraint for node v and
146  * commodity k; NULL if this node does not exist in the commodity */
147  SCIP_Real** nodeflowscales; /**< scaling factors to convert nodeflowrows[v][k] into a +/-1 <= row */
148  SCIP_Bool** nodeflowinverted; /**< does nodeflowrows[v][k] have to be inverted to fit the network structure? */
149  SCIP_ROW** arccapacityrows; /**< arccapacity[a]: capacity constraint on arc a;
150  * NULL if uncapacitated */
151  SCIP_Real* arccapacityscales; /**< scaling factors to convert arccapacity[a] into a <= row with
152  * positive entries for the flow variables */
153  int* arcsources; /**< source node ids of arcs */
154  int* arctargets; /**< target node ids of arcs */
155  int* colcommodity; /**< commodity number of each column, or -1 */
156  int nnodes; /**< number of nodes in the graph */
157  int narcs; /**< number of arcs in the graph */
158  int nuncapacitatedarcs; /**< number of uncapacitated arcs in the graph */
159  int ncommodities; /**< number of commodities */
160  SCIP_MCFMODELTYPE modeltype; /**< detected model type of the network */
161 };
163 
164 /** separator data */
165 struct SCIP_SepaData
166 {
167  SCIP_MCFNETWORK** mcfnetworks; /**< array of multi-commodity-flow network structures */
168  int nmcfnetworks; /**< number of multi-commodity-flow networks (-1: extraction not yet done) */
169  int nclusters; /**< number of clusters to generate in the shrunken network -- default separation */
170  SCIP_Real maxweightrange; /**< maximal valid range max(|weights|)/min(|weights|) of row weights */
171  int maxtestdelta; /**< maximal number of different deltas to try (-1: unlimited) -- default separation */
172  SCIP_Bool trynegscaling; /**< should negative values also be tested in scaling? */
173  SCIP_Bool fixintegralrhs; /**< should an additional variable be complemented if f0 = 0? */
174  SCIP_Bool dynamiccuts; /**< should generated cuts be removed from the LP if they are no longer tight? */
175  int modeltype; /**< model type of the network */
176  int maxsepacuts; /**< maximal number of cmir cuts separated per separation round */
177  int maxsepacutsroot; /**< maximal number of cmir cuts separated per separation round in root node -- default separation */
178  SCIP_Real maxinconsistencyratio; /**< maximum inconsistency ratio (inconsistencies/(arcs*commodities)) for separation at all*/
179  SCIP_Real maxarcinconsistencyratio; /**< maximum inconsistency ratio for arcs not to be deleted */
180  SCIP_Bool checkcutshoreconnectivity;/**< should we only separate if the cuts shores are connected */
181  SCIP_Bool separatesinglenodecuts; /**< should we separate inequalities based on single node cuts ? */
182  SCIP_Bool separateflowcutset; /**< should we separate flowcutset inequalities on the network cuts ? */
183  SCIP_Bool separateknapsack; /**< should we separate knapsack cover inequalities on the network cuts ? */
184  SCIP_Bool lastroundsuccess; /**< did we find a cut in the last round? */
185  MCFEFFORTLEVEL effortlevel; /**< effort level of separation (off / aggressive / default) */
186 };
187 
188 /** internal MCF extraction data to pass to subroutines */
189 struct mcfdata
190 {
191  unsigned char* flowrowsigns; /**< potential or actual sides of rows to be used as flow conservation constraint */
192  SCIP_Real* flowrowscalars; /**< scalar of rows to transform into +/-1 coefficients */
193  SCIP_Real* flowrowscores; /**< score value indicating how sure we are that this is indeed a flow conservation constraint */
194  unsigned char* capacityrowsigns; /**< potential or actual sides of rows to be used as capacity constraint */
195  SCIP_Real* capacityrowscores; /**< score value indicating how sure we are that this is indeed a capacity constraint */
196  int* flowcands; /**< list of row indices that are candidates for flow conservation constraints */
197  int nflowcands; /**< number of elements in flow candidate list */
198  int* capacitycands; /**< list of row indices that are candidates for capacity constraints */
199  int ncapacitycands; /**< number of elements in capacity candidate list */
200  SCIP_Bool* plusflow; /**< is column c member of a flow row with coefficient +1? */
201  SCIP_Bool* minusflow; /**< is column c member of a flow row with coefficient -1? */
202  int ncommodities; /**< number of commodities */
203  int nemptycommodities; /**< number of commodities that have been discarded but still counted in 'ncommodities' */
204  int* commoditysigns; /**< +1: regular, -1: all arcs have opposite direction; 0: undecided */
205  int commoditysignssize; /**< size of commoditysigns array */
206  int* colcommodity; /**< commodity number of each column, or -1 */
207  int* rowcommodity; /**< commodity number of each row, or -1 */
208  int* colarcid; /**< arc id of each flow column, or -1 */
209  int* rowarcid; /**< arc id of each capacity row, or -1 */
210  int* rownodeid; /**< node id of each flow conservation row, or -1 */
211  int arcarraysize; /**< size of arrays indexed by arcs */
212  int* arcsources; /**< source node ids of arcs */
213  int* arctargets; /**< target node ids of arcs */
214  int* colsources; /**< source node for flow columns, -1 for non existing, -2 for unknown */
215  int* coltargets; /**< target node for flow columns, -1 for non existing, -2 for unknown */
216  int* firstoutarcs; /**< for each node the first arc id for which the node is the source node */
217  int* firstinarcs; /**< for each node the first arc id for which the node is the target node */
218  int* nextoutarcs; /**< for each arc the next outgoing arc in the adjacency list */
219  int* nextinarcs; /**< for each arc the next outgoing arc in the adjacency list */
220  int* newcols; /**< columns of current commodity that have to be inspected for incident flow conservation rows */
221  int nnewcols; /**< number of newcols */
222  int narcs; /**< number of arcs in the extracted graph */
223  int nnodes; /**< number of nodes in the extracted graph */
224  SCIP_Real ninconsistencies; /**< number of inconsistencies between the commodity graphs */
225  SCIP_ROW** capacityrows; /**< capacity row for each arc */
226  int capacityrowssize; /**< size of array */
227  SCIP_Bool* colisincident; /**< temporary memory for column collection */
228  int* zeroarcarray; /**< temporary array of zeros */
229  SCIP_MCFMODELTYPE modeltype; /**< model type that is used for this network extraction */
230 };
231 typedef struct mcfdata MCFDATA; /**< internal MCF extraction data to pass to subroutines */
232 
233 /** data structure to put a nodepair on the heap -- it holds node1 <= node2 */
234 struct nodepair
235 {
236  int node1; /**< index of the first node */
237  int node2; /**< index of the second node */
238  SCIP_Real weight; /**< weight of the arc in the separation problem */
239 };
240 typedef struct nodepair NODEPAIRENTRY;
241 
242 /** nodepair priority queue */
243 struct nodepairqueue
244 {
245  SCIP_PQUEUE* pqueue; /**< priority queue of elements */
246  NODEPAIRENTRY* nodepairs; /**< elements on the heap */
247 };
248 typedef struct nodepairqueue NODEPAIRQUEUE;
249 
250 /** partitioning of the nodes into clusters */
251 struct nodepartition
252 {
253  int* representatives; /**< mapping of node ids to their representatives within their cluster */
254  int* nodeclusters; /**< cluster for each node id */
255  int* clusternodes; /**< node ids sorted by cluster */
256  int* clusterbegin; /**< first entry in clusternodes for each cluster (size: nclusters+1) */
257  int nclusters; /**< number of clusters */
258 };
259 typedef struct nodepartition NODEPARTITION;
260 
261 /*
262  * Local methods
263  */
264 
265 #define LHSPOSSIBLE 1u /**< we may use the constraint as lhs <= a*x */
266 #define RHSPOSSIBLE 2u /**< we may use the constraint as a*x <= rhs */
267 #define LHSASSIGNED 4u /**< we have chosen to use the constraint as lhs <= a*x */
268 #define RHSASSIGNED 8u /**< we have chosen to use the constraint as a*x <= rhs */
269 #define INVERTED 16u /**< we need to invert the row */
270 #define DISCARDED 32u /**< we have chosen to not use the constraint */
271 #define UNDIRECTED 64u /**< the capacity candidate has two flow variables for a commodity */
272 
273 
274 /** creates an empty MCF network data structure */
275 static
277  SCIP* scip, /**< SCIP data structure */
278  SCIP_MCFNETWORK** mcfnetwork /**< MCF network structure */
279  )
280 {
281  assert(mcfnetwork != NULL);
282 
283  SCIP_CALL( SCIPallocMemory(scip, mcfnetwork) );
284  (*mcfnetwork)->nodeflowrows = NULL;
285  (*mcfnetwork)->nodeflowscales = NULL;
286  (*mcfnetwork)->nodeflowinverted = NULL;
287  (*mcfnetwork)->arccapacityrows = NULL;
288  (*mcfnetwork)->arccapacityscales = NULL;
289  (*mcfnetwork)->arcsources = NULL;
290  (*mcfnetwork)->arctargets = NULL;
291  (*mcfnetwork)->colcommodity = NULL;
292  (*mcfnetwork)->nnodes = 0;
293  (*mcfnetwork)->nuncapacitatedarcs = 0;
294  (*mcfnetwork)->narcs = 0;
295  (*mcfnetwork)->ncommodities = 0;
296 
297  return SCIP_OKAY;
298 }
299 
300 /** frees MCF network data structure */
301 static
303  SCIP* scip, /**< SCIP data structure */
304  SCIP_MCFNETWORK** mcfnetwork /**< MCF network structure */
305  )
306 {
307  assert(mcfnetwork != NULL);
308 
309  if( *mcfnetwork != NULL )
310  {
311  int v;
312  int a;
313 
314  for( v = 0; v < (*mcfnetwork)->nnodes; v++ )
315  {
316  int k;
317 
318  for( k = 0; k < (*mcfnetwork)->ncommodities; k++ )
319  {
320  if( (*mcfnetwork)->nodeflowrows[v][k] != NULL )
321  {
322  SCIP_CALL( SCIPreleaseRow(scip, &(*mcfnetwork)->nodeflowrows[v][k]) );
323  }
324  }
325  SCIPfreeMemoryArrayNull(scip, &(*mcfnetwork)->nodeflowrows[v]);
326  SCIPfreeMemoryArrayNull(scip, &(*mcfnetwork)->nodeflowscales[v]);
327  SCIPfreeMemoryArrayNull(scip, &(*mcfnetwork)->nodeflowinverted[v]);
328  }
329  for( a = 0; a < (*mcfnetwork)->narcs; a++ )
330  {
331  if( (*mcfnetwork)->arccapacityrows[a] != NULL )
332  {
333  SCIP_CALL( SCIPreleaseRow(scip, &(*mcfnetwork)->arccapacityrows[a]) );
334  }
335  }
336  SCIPfreeMemoryArrayNull(scip, &(*mcfnetwork)->nodeflowrows);
337  SCIPfreeMemoryArrayNull(scip, &(*mcfnetwork)->nodeflowscales);
338  SCIPfreeMemoryArrayNull(scip, &(*mcfnetwork)->nodeflowinverted);
339  SCIPfreeMemoryArrayNull(scip, &(*mcfnetwork)->arccapacityrows);
340  SCIPfreeMemoryArrayNull(scip, &(*mcfnetwork)->arccapacityscales);
341  SCIPfreeMemoryArrayNull(scip, &(*mcfnetwork)->arcsources);
342  SCIPfreeMemoryArrayNull(scip, &(*mcfnetwork)->arctargets);
343  SCIPfreeMemoryArrayNull(scip, &(*mcfnetwork)->colcommodity);
344 
345  SCIPfreeMemory(scip, mcfnetwork);
346  }
347 
348  return SCIP_OKAY;
349 }
350 
351 /** fills the MCF network structure with the MCF data */
352 static
354  SCIP* scip, /**< SCIP data structure */
355  SCIP_MCFNETWORK* mcfnetwork, /**< MCF network structure */
356  MCFDATA* mcfdata, /**< internal MCF extraction data to pass to subroutines */
357  int* compnodeid, /**< temporary storage for v -> compv mapping; must be set to -1 for all v */
358  int* compnodes, /**< array of node ids of the component */
359  int ncompnodes, /**< number of nodes in the component */
360  int* comparcs, /**< array of arc ids of the component */
361  int ncomparcs /**< number of arcs in the component */
362  )
363 {
364  unsigned char* flowrowsigns = mcfdata->flowrowsigns;
365  SCIP_Real* flowrowscalars = mcfdata->flowrowscalars;
366  unsigned char* capacityrowsigns = mcfdata->capacityrowsigns;
367  int* flowcands = mcfdata->flowcands;
368  int nflowcands = mcfdata->nflowcands;
369  int ncommodities = mcfdata->ncommodities;
370  int* commoditysigns = mcfdata->commoditysigns;
371  int* colcommodity = mcfdata->colcommodity;
372  int* rowcommodity = mcfdata->rowcommodity;
373  int* rownodeid = mcfdata->rownodeid;
374  SCIP_ROW** capacityrows = mcfdata->capacityrows;
375  SCIP_MCFMODELTYPE modeltype = mcfdata->modeltype;
376 
377  SCIP_Real* comdemands;
378  SCIP_ROW** rows;
379  SCIP_COL** cols;
380  int nrows;
381  int ncols;
382  int* compcommodity;
383  int ncompcommodities;
384  int v;
385  int a;
386  int k;
387  int i;
388  int c;
389 
390  assert(mcfnetwork != NULL);
391  assert(modeltype != SCIP_MCFMODELTYPE_AUTO);
392  assert(2 <= ncompnodes && ncompnodes <= mcfdata->nnodes);
393  assert(1 <= ncomparcs && ncomparcs <= mcfdata->narcs);
394  assert(ncommodities > 0);
395 
396 #ifndef NDEBUG
397  /* v -> compv mapping must be all -1 */
398  for( v = 0; v < mcfdata->nnodes; v++ )
399  assert(compnodeid[v] == -1);
400 #endif
401 
402  /* allocate temporary memory */
403  SCIP_CALL( SCIPallocBufferArray(scip, &comdemands, ncommodities) );
404  SCIP_CALL( SCIPallocBufferArray(scip, &compcommodity, ncommodities) );
405 
406  /* initialize demand array */
407  BMSclearMemoryArray(comdemands, ncommodities);
408 
409  /* initialize k -> compk mapping */
410  for( k = 0; k < ncommodities; k++ )
411  compcommodity[k] = -1;
412 
413  /* get LP rows and cols data */
414  SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
415  SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
416 
417  /* generate v -> compv mapping */
418  for( i = 0; i < ncompnodes; i++ )
419  {
420  v = compnodes[i];
421  assert(0 <= v && v < mcfdata->nnodes);
422  compnodeid[v] = i;
423  }
424 
425  /* generate k -> compk mapping */
426  ncompcommodities = 0;
427  for( i = 0; i < nflowcands; i++ )
428  {
429  int r;
430  int rv;
431 
432  r = flowcands[i];
433  assert(0 <= r && r < nrows);
434 
435  rv = rownodeid[r];
436  if( rv >= 0 && compnodeid[rv] >= 0 )
437  {
438  k = rowcommodity[r];
439  assert(0 <= k && k < ncommodities);
440  if( compcommodity[k] == -1 )
441  {
442  compcommodity[k] = ncompcommodities;
443  ncompcommodities++;
444  }
445  }
446  }
447 
448  /** @todo model type and flow type may be different for each component */
449  /* record model and flow type */
450  mcfnetwork->modeltype = modeltype;
451 
452  /* record network size */
453  mcfnetwork->nnodes = ncompnodes;
454  mcfnetwork->narcs = ncomparcs;
455  mcfnetwork->nuncapacitatedarcs = 0;
456  mcfnetwork->ncommodities = ncompcommodities;
457 
458  /* allocate memory for arrays and initialize with default values */
459  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfnetwork->nodeflowrows, mcfnetwork->nnodes) );
460  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfnetwork->nodeflowscales, mcfnetwork->nnodes) );
461  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfnetwork->nodeflowinverted, mcfnetwork->nnodes) );
462  for( v = 0; v < mcfnetwork->nnodes; v++ )
463  {
464  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfnetwork->nodeflowrows[v], mcfnetwork->ncommodities) ); /*lint !e866*/
465  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfnetwork->nodeflowscales[v], mcfnetwork->ncommodities) ); /*lint !e866*/
466  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfnetwork->nodeflowinverted[v], mcfnetwork->ncommodities) ); /*lint !e866*/
467  for( k = 0; k < mcfnetwork->ncommodities; k++ )
468  {
469  mcfnetwork->nodeflowrows[v][k] = NULL;
470  mcfnetwork->nodeflowscales[v][k] = 0.0;
471  mcfnetwork->nodeflowinverted[v][k] = FALSE;
472  }
473  }
474 
475  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfnetwork->arccapacityrows, mcfnetwork->narcs) );
476  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfnetwork->arccapacityscales, mcfnetwork->narcs) );
477  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfnetwork->arcsources, mcfnetwork->narcs) );
478  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfnetwork->arctargets, mcfnetwork->narcs) );
479  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfnetwork->colcommodity, ncols) );
480  for( a = 0; a < mcfnetwork->narcs; a++ )
481  {
482  mcfnetwork->arccapacityrows[a] = NULL;
483  mcfnetwork->arccapacityscales[a] = 0.0;
484  mcfnetwork->arcsources[a] = -1;
485  mcfnetwork->arctargets[a] = -1;
486  }
487  BMSclearMemoryArray(mcfnetwork->colcommodity, mcfnetwork->ncommodities);
488 
489  /* fill in existing node data */
490  for( i = 0; i < nflowcands; i++ )
491  {
492  int r;
493  int rv;
494 
495  r = flowcands[i];
496  assert(0 <= r && r < nrows);
497 
498  rv = rownodeid[r];
499  if( rv >= 0 && compnodeid[rv] >= 0 )
500  {
501  SCIP_Real scale;
502  int rk;
503 
504  v = compnodeid[rv];
505  rk = rowcommodity[r];
506  assert(v < mcfnetwork->nnodes);
507  assert(0 <= rk && rk < ncommodities);
508  assert((flowrowsigns[r] & (LHSASSIGNED | RHSASSIGNED)) != 0);
509 
510  k = compcommodity[rk];
511  assert(0 <= k && k < mcfnetwork->ncommodities);
512 
513  /* fill in node -> row assignment */
514  SCIP_CALL( SCIPcaptureRow(scip, rows[r]) );
515  mcfnetwork->nodeflowrows[v][k] = rows[r];
516  scale = flowrowscalars[r];
517  if( (flowrowsigns[r] & LHSASSIGNED) != 0 )
518  scale *= -1.0;
519  if( commoditysigns[rk] == -1 )
520  scale *= -1.0;
521  mcfnetwork->nodeflowscales[v][k] = scale;
522  mcfnetwork->nodeflowinverted[v][k] = ((flowrowsigns[r] & INVERTED) != 0);
523  }
524  }
525 
526  /* fill in existing arc data */
527  for( a = 0; a < mcfnetwork->narcs; a++ )
528  {
529  SCIP_ROW* capacityrow;
530  SCIP_COL** rowcols;
531  SCIP_Real* rowvals;
532  int rowlen;
533  int globala;
534  int r;
535  int j;
536 
537  globala = comparcs[a];
538  capacityrow = capacityrows[globala];
539 
540  mcfnetwork->arccapacityscales[a] = 1.0;
541 
542  /* If arc is capacitated */
543  if( capacityrow != NULL)
544  {
545  r = SCIProwGetLPPos(capacityrow);
546  assert(0 <= r && r < nrows);
547  assert((capacityrowsigns[r] & (LHSASSIGNED | RHSASSIGNED)) != 0);
548  assert((capacityrowsigns[r] & INVERTED) == 0);
549  assert(mcfdata->rowarcid[r] == globala);
550 
551  SCIP_CALL( SCIPcaptureRow(scip, capacityrow) );
552  mcfnetwork->arccapacityrows[a] = capacityrow;
553 
554  /* Scale constraint such that the coefficients for the flow variables are normalized in such a way that coefficients in
555  * multiple capacity constraints that belong to the same commodity are (hopefully) equal.
556  * This is needed for binary flow variable models in which the demand of each commodity is stored as the coefficient in
557  * the capacity constraints. Since it may happen (e.g., due to presolve) that different capacity constraints are scaled
558  * differently, we need to find scaling factors to make the demand coefficients of each commodity equal.
559  * To do this, we scale the first capacity constraint with +1 and then store the coefficients of the flow variables
560  * as target demands for the commodities. Then, we scale the other constraints in such a way that these demands are hit, if possible.
561  * Note that for continuous flow variable models, the coefficients in the capacity constraints are usually +1.0.
562  * This is achieved automatically by our scaling routine.
563  */
564  rowcols = SCIProwGetCols(capacityrow);
565  rowvals = SCIProwGetVals(capacityrow);
566  rowlen = SCIProwGetNLPNonz(capacityrow);
567  for( j = 0; j < rowlen; j++ )
568  {
569  c = SCIPcolGetLPPos(rowcols[j]);
570  assert(0 <= c && c < SCIPgetNLPCols(scip));
571  k = colcommodity[c];
572  if( k >= 0 )
573  {
574  if( comdemands[k] != 0.0 )
575  {
576  /* update the scaling factor */
577  mcfnetwork->arccapacityscales[a] = comdemands[k]/rowvals[j];
578  break;
579  }
580  }
581  }
582 
583  /* use negative scaling if we use the left hand side, use positive scaling if we use the right hand side */
584  mcfnetwork->arccapacityscales[a] = ABS(mcfnetwork->arccapacityscales[a]);
585  if( (capacityrowsigns[r] & LHSASSIGNED) != 0 )
586  mcfnetwork->arccapacityscales[a] *= -1.0;
587 
588  /* record the commodity demands */
589  for( j = 0; j < rowlen; j++ )
590  {
591  c = SCIPcolGetLPPos(rowcols[j]);
592  assert(0 <= c && c < SCIPgetNLPCols(scip));
593  k = colcommodity[c];
594  if( k >= 0 && comdemands[k] == 0.0 )
595  comdemands[k] = mcfnetwork->arccapacityscales[a] * rowvals[j];
596  }
597  }
598  else
599  {
600  /* arc is uncapacitated */
601  mcfnetwork->arccapacityrows[a] = NULL;
602  mcfnetwork->nuncapacitatedarcs++;
603  }
604 
605  /* copy the source/target node assignment */
606  if( mcfdata->arcsources[globala] >= 0 )
607  {
608  assert(mcfdata->arcsources[globala] < mcfdata->nnodes);
609  assert(0 <= compnodeid[mcfdata->arcsources[globala]] && compnodeid[mcfdata->arcsources[globala]] < mcfnetwork->nnodes);
610  mcfnetwork->arcsources[a] = compnodeid[mcfdata->arcsources[globala]];
611  }
612  if( mcfdata->arctargets[globala] >= 0 )
613  {
614  assert(mcfdata->arctargets[globala] < mcfdata->nnodes);
615  assert(0 <= compnodeid[mcfdata->arctargets[globala]] && compnodeid[mcfdata->arctargets[globala]] < mcfnetwork->nnodes);
616  mcfnetwork->arctargets[a] = compnodeid[mcfdata->arctargets[globala]];
617  }
618  }
619 
620  /* translate colcommodity array */
621  for( c = 0; c < ncols; c++ )
622  {
623  k = colcommodity[c];
624  if( k >= 0 )
625  mcfnetwork->colcommodity[c] = compcommodity[k];
626  else
627  mcfnetwork->colcommodity[c] = -1;
628  }
629 
630  /* reset v -> compv mapping */
631  for( i = 0; i < ncompnodes; i++ )
632  {
633  assert(0 <= compnodes[i] && compnodes[i] < mcfdata->nnodes);
634  assert(compnodeid[compnodes[i]] == i);
635  compnodeid[compnodes[i]] = -1;
636  }
637 
638  /* free temporary memory */
639  SCIPfreeBufferArray(scip, &compcommodity);
640  SCIPfreeBufferArray(scip, &comdemands);
641 
642  return SCIP_OKAY;
643 }
644 
645 #ifdef SCIP_DEBUG
646 /** displays the MCF network */
647 static
648 void mcfnetworkPrint(
649  SCIP_MCFNETWORK* mcfnetwork /**< MCF network structure */
650  )
651 {
652  if( mcfnetwork == NULL )
653  MCFdebugMessage("MCF network is empty\n");
654  else
655  {
656  int v;
657  int a;
658 
659  for( v = 0; v < mcfnetwork->nnodes; v++ )
660  {
661  int k;
662 
663  MCFdebugMessage("node %2d:\n", v);
664  for( k = 0; k < mcfnetwork->ncommodities; k++ )
665  {
666  MCFdebugMessage(" commodity %2d: ", k);
667  if( mcfnetwork->nodeflowrows[v][k] != NULL )
668  {
669  MCFdebugMessage("<%s> [%+g] [inv:%u]\n", SCIProwGetName(mcfnetwork->nodeflowrows[v][k]),
670  mcfnetwork->nodeflowscales[v][k], mcfnetwork->nodeflowinverted[v][k]);
671  /*SCIP_CALL( SCIProwPrint(mcfnetwork->nodeflowrows[v][k], NULL) );*/
672  }
673  else
674  MCFdebugMessage("-\n");
675  }
676  }
677 
678  for( a = 0; a < mcfnetwork->narcs; a++ )
679  {
680  MCFdebugMessage("arc %2d [%2d -> %2d]: ", a, mcfnetwork->arcsources[a], mcfnetwork->arctargets[a]);
681  if( mcfnetwork->arccapacityrows[a] != NULL )
682  {
683  MCFdebugMessage("<%s> [%+g]\n", SCIProwGetName(mcfnetwork->arccapacityrows[a]), mcfnetwork->arccapacityscales[a]);
684  /*SCIProwPrint(mcfnetwork->arccapacityrows[a], NULL);*/
685  }
686  else
687  MCFdebugMessage("-\n");
688  }
689  }
690 }
691 
692 /** displays commodities and its members */
693 static
694 void printCommodities(
695  SCIP* scip, /**< SCIP data structure */
696  MCFDATA* mcfdata /**< internal MCF extraction data to pass to subroutines */
697  )
698 {
699  unsigned char* flowrowsigns = mcfdata->flowrowsigns;
700  unsigned char* capacityrowsigns = mcfdata->capacityrowsigns;
701  int ncommodities = mcfdata->ncommodities;
702  int* commoditysigns = mcfdata->commoditysigns;
703  int* colcommodity = mcfdata->colcommodity;
704  int* rowcommodity = mcfdata->rowcommodity;
705  int* colarcid = mcfdata->colarcid;
706  int* rownodeid = mcfdata->rownodeid;
707  int nnodes = mcfdata->nnodes;
708  SCIP_ROW** capacityrows = mcfdata->capacityrows;
709 
710  SCIP_COL** cols;
711  SCIP_ROW** rows;
712  int ncols;
713  int nrows;
714  int k;
715  int c;
716  int r;
717  int a;
718 
719  cols = SCIPgetLPCols(scip);
720  ncols = SCIPgetNLPCols(scip);
721  rows = SCIPgetLPRows(scip);
722  nrows = SCIPgetNLPRows(scip);
723 
724  for( k = 0; k < ncommodities; k++ )
725  {
726  MCFdebugMessage("commodity %d (sign: %+d):\n", k, commoditysigns[k]);
727 
728  for( c = 0; c < ncols; c++ )
729  {
730  if( colcommodity[c] == k )
731  MCFdebugMessage(" col <%s>: arc %d\n", SCIPvarGetName(SCIPcolGetVar(cols[c])), colarcid != NULL ? colarcid[c] : -1);
732  }
733  for( r = 0; r < nrows; r++ )
734  {
735  if( rowcommodity[r] == k )
736  MCFdebugMessage(" row <%s>: node %d [sign:%+d, inv:%+d]\n", SCIProwGetName(rows[r]), rownodeid != NULL ? rownodeid[r] : -1,
737  (flowrowsigns[r] & RHSASSIGNED) != 0 ? +1 : -1,
738  (flowrowsigns[r] & INVERTED) != 0 ? -1 : +1);
739  }
740  MCFdebugMessage("\n");
741  }
742 
743  if( rownodeid != NULL )
744  {
745  int v;
746 
747  for( v = 0; v < nnodes; v++ )
748  {
749  MCFdebugMessage("node %d:\n", v);
750  for( r = 0; r < nrows; r++ )
751  {
752  if( rownodeid[r] == v )
753  MCFdebugMessage(" row <%s> [sign:%+d, inv:%+d]\n", SCIProwGetName(rows[r]),
754  (flowrowsigns[r] & RHSASSIGNED) != 0 ? +1 : -1,
755  (flowrowsigns[r] & INVERTED) != 0 ? -1 : +1);
756  }
757  MCFdebugMessage("\n");
758  }
759  }
760 
761  assert(capacityrows != NULL || mcfdata->narcs == 0);
762 
763  MCFdebugMessage("capacities:\n");
764  for( a = 0; a < mcfdata->narcs; a++ )
765  {
766  MCFdebugMessage(" arc %d: ", a);
767  if( capacityrows[a] != NULL ) /*lint !e613*/
768  {
769  r = SCIProwGetLPPos(capacityrows[a]); /*lint !e613*/
770  assert(0 <= r && r < nrows);
771  if( (capacityrowsigns[r] & LHSASSIGNED) != 0 )
772  MCFdebugMessage(" row <%s> [sign:-1]\n", SCIProwGetName(rows[r]));
773  else if( (capacityrowsigns[r] & RHSASSIGNED) != 0 )
774  MCFdebugMessage(" row <%s> [sign:+1]\n", SCIProwGetName(rows[r]));
775  }
776  else
777  MCFdebugMessage(" -\n");
778  }
779  MCFdebugMessage("\n");
780 
781  assert(colcommodity != NULL || ncols == 0);
782 
783  MCFdebugMessage("unused columns:\n");
784  for( c = 0; c < ncols; c++ )
785  {
786  if( colcommodity[c] == -1 ) /*lint !e613*/
787  {
788  SCIP_VAR* var = SCIPcolGetVar(cols[c]);
789  MCFdebugMessage(" col <%s> [%g,%g]\n", SCIPvarGetName(var), SCIPvarGetLbGlobal(var), SCIPvarGetUbGlobal(var));
790  }
791  }
792  MCFdebugMessage("\n");
793 
794  MCFdebugMessage("unused rows:\n");
795  for( r = 0; r < nrows; r++ )
796  {
797  assert(rowcommodity != NULL);
798 
799  if( rowcommodity[r] == -1 && (capacityrowsigns == NULL || (capacityrowsigns[r] & (LHSASSIGNED | RHSASSIGNED)) == 0) )
800  {
801  MCFdebugMessage(" row <%s>\n", SCIProwGetName(rows[r]));
802  /*SCIPdebug( SCIP_CALL(SCIPprintRow(scip, rows[r], NULL)) );*/
803  }
804  }
805  MCFdebugMessage("\n");
806 }
807 #endif
808 
809 /** comparator method for flow and capacity row candidates */
810 static
812 {
813  SCIP_Real* rowscores = (SCIP_Real*)dataptr;
814 
815  if( rowscores[ind2] < rowscores[ind1] )
816  return -1;
817  else if( rowscores[ind2] > rowscores[ind1] )
818  return +1;
819  else
820  return 0;
821 }
822 
823 /** extracts flow conservation from the LP */
824 static
826  SCIP* scip, /**< SCIP data structure */
827  MCFDATA* mcfdata /**< internal MCF extraction data to pass to subroutines */
828  )
829 {
830  unsigned char* flowrowsigns;
831  SCIP_Real* flowrowscalars;
832  SCIP_Real* flowrowscores;
833  int* flowcands;
834 
835  SCIP_ROW** rows;
836  int nrows;
837  int ncols;
838  int r;
839 
840  SCIP_Real maxdualflow;
841 
842  SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
843  ncols = SCIPgetNLPCols(scip);
844 
845  /* allocate temporary memory for extraction data */
846  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->flowrowsigns, nrows) ); /*lint !e685*/
847  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->flowrowscalars, nrows) );
848  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->flowrowscores, nrows) );
849  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->flowcands, nrows) );
850  flowrowsigns = mcfdata->flowrowsigns;
851  flowrowscalars = mcfdata->flowrowscalars;
852  flowrowscores = mcfdata->flowrowscores;
853  flowcands = mcfdata->flowcands;
854 
855  assert(mcfdata->nflowcands == 0);
856 
857  maxdualflow = 0.0;
858  for( r = 0; r < nrows; r++ )
859  {
860  SCIP_ROW* row;
861  SCIP_COL** rowcols;
862  SCIP_Real* rowvals;
863  SCIP_Real rowlhs;
864  SCIP_Real rowrhs;
865  int rowlen;
866  int nbinvars;
867  int nintvars;
868  int nimplintvars;
869  int ncontvars;
870  SCIP_Real coef;
871  SCIP_Bool hasposcoef;
872  SCIP_Bool hasnegcoef;
873  SCIP_Real absdualsol;
874  int i;
875 
876  row = rows[r];
877  assert(SCIProwGetLPPos(row) == r);
878 
879  /* get dual solution, if available */
880  absdualsol = SCIProwGetDualsol(row);
881  if( absdualsol == SCIP_INVALID ) /*lint !e777*/
882  absdualsol = 0.0;
883  absdualsol = ABS(absdualsol);
884 
885  flowrowsigns[r] = 0;
886  flowrowscalars[r] = 0.0;
887  flowrowscores[r] = 0.0;
888 
889  /* ignore modifiable rows */
890  if( SCIProwIsModifiable(row) )
891  continue;
892 
893  /* ignore empty rows */
894  rowlen = SCIProwGetNLPNonz(row);
895  if( rowlen == 0 )
896  continue;
897 
898  /* No dense rows please */
899  if( rowlen > MAXFLOWCANDDENSITY * ncols )
900  continue;
901 
902  rowcols = SCIProwGetCols(row);
903  rowvals = SCIProwGetVals(row);
904  rowlhs = SCIProwGetLhs(row) - SCIProwGetConstant(row);
905  rowrhs = SCIProwGetRhs(row) - SCIProwGetConstant(row);
906 
907  /* identify flow conservation constraints */
908  coef = ABS(rowvals[0]);
909  hasposcoef = FALSE;
910  hasnegcoef = FALSE;
911  nbinvars = 0;
912  nintvars = 0;
913  nimplintvars = 0;
914  ncontvars = 0;
915  for( i = 0; i < rowlen; i++ )
916  {
917  SCIP_Real absval = ABS(rowvals[i]);
918  if( !SCIPisEQ(scip, absval, coef) )
919  break;
920 
921  hasposcoef = hasposcoef || (rowvals[i] > 0.0);
922  hasnegcoef = hasnegcoef || (rowvals[i] < 0.0);
923  switch( SCIPvarGetType(SCIPcolGetVar(rowcols[i])) )
924  {
925  case SCIP_VARTYPE_BINARY:
926  nbinvars++;
927  break;
929  nintvars++;
930  break;
932  nimplintvars++;
933  break;
935  ncontvars++;
936  break;
937  default:
938  SCIPerrorMessage("unknown variable type\n");
939  SCIPABORT();
940  return SCIP_INVALIDDATA; /*lint !e527*/
941  }
942  }
943  if( i == rowlen )
944  {
945  /* Flow conservation constraints should always be a*x <= -d.
946  * If lhs and rhs are finite, both sides are still valid candidates.
947  */
948  if( !SCIPisInfinity(scip, -rowlhs) )
949  flowrowsigns[r] |= LHSPOSSIBLE;
950  if( !SCIPisInfinity(scip, rowrhs) )
951  flowrowsigns[r] |= RHSPOSSIBLE;
952  flowrowscalars[r] = 1.0/coef;
953  flowcands[mcfdata->nflowcands] = r;
954  mcfdata->nflowcands++;
955  }
956 
957  /* calculate flow row score */
958  if( (flowrowsigns[r] & (LHSPOSSIBLE | RHSPOSSIBLE)) != 0 )
959  {
960  /* row does not need to be scaled: score +1000 */
961  if( SCIPisEQ(scip, flowrowscalars[r], 1.0) )
962  flowrowscores[r] += 1000.0;
963 
964  /* row has positive and negative coefficients: score +500 */
965  if( hasposcoef && hasnegcoef )
966  flowrowscores[r] += 500.0;
967 
968  /* all variables are of the same type:
969  * continuous: score +1000
970  * integer: score +500
971  * binary: score +100
972  */
973  if( ncontvars == rowlen )
974  flowrowscores[r] += 1000.0;
975  else if( nintvars + nimplintvars == rowlen )
976  flowrowscores[r] += 500.0;
977  else if( nbinvars == rowlen )
978  flowrowscores[r] += 100.0;
979 
980  /* the longer the row, the earlier we want to process it: score +10*len/(len+10) */
981  /* value is in [1,10) */
982  flowrowscores[r] += 10.0*rowlen/(rowlen+10.0);
983 
984  /* row is an equation: score +50, tie-breaking */
985  if( (flowrowsigns[r] & (LHSPOSSIBLE | RHSPOSSIBLE)) == (LHSPOSSIBLE | RHSPOSSIBLE) )
986  flowrowscores[r] += 50.0;
987 
988  assert(flowrowscores[r] > 0.0);
989 
990  /* update maximum dual solution value for additional score tie breaking */
991  maxdualflow = MAX(maxdualflow, absdualsol);
992 
993  /** @todo go through list of several model types, depending on the current model type throw away invalid constraints
994  * instead of assigning a low score
995  */
996  }
997  }
998 
999  /* apply additional score tie breaking using the dual solutions */
1000  if( SCIPisPositive(scip, maxdualflow) )
1001  {
1002  int i;
1003 
1004  for( i = 0; i < mcfdata->nflowcands; i++ )
1005  {
1006  SCIP_Real dualsol;
1007 
1008  r = flowcands[i];
1009  assert(0 <= r && r < nrows);
1010  dualsol = SCIProwGetDualsol(rows[r]);
1011  if( dualsol == SCIP_INVALID ) /*lint !e777*/
1012  dualsol = 0.0;
1013  else if( flowrowsigns[r] == (LHSPOSSIBLE | RHSPOSSIBLE) )
1014  dualsol = ABS(dualsol);
1015  else if( flowrowsigns[r] == RHSPOSSIBLE )
1016  dualsol = -dualsol;
1017  assert(maxdualflow > 0.0); /*for flexelint*/
1018  flowrowscores[r] += dualsol/maxdualflow + 1.0;
1019  assert(flowrowscores[r] > 0.0);
1020  }
1021  }
1022 
1023  /* sort candidates by score */
1024  SCIPsortInd(mcfdata->flowcands, compCands, (void*)flowrowscores, mcfdata->nflowcands);
1025 
1026  MCFdebugMessage("flow conservation candidates [%d]\n", mcfdata->nflowcands);
1027 #ifdef SCIP_DEBUG
1028  for( r = 0; r < mcfdata->nflowcands; r++ )
1029  {
1030  /*SCIPdebug( SCIP_CALL(SCIPprintRow(scip, rows[mcfdata->flowcands[r]], NULL)) );*/
1031  SCIPdebugMsg(scip, "%4d [score: %2g]: %s\n", mcfdata->flowcands[r], flowrowscores[mcfdata->flowcands[r]],
1032  SCIProwGetName(rows[mcfdata->flowcands[r]]));
1033  }
1034 #endif
1035 
1036  return SCIP_OKAY;
1037 }
1038 
1039 /** extracts capacity rows from the LP */
1040 static
1042  SCIP* scip, /**< SCIP data structure */
1043  MCFDATA* mcfdata /**< internal MCF extraction data to pass to subroutines */
1044  )
1045 {
1046  unsigned char* flowrowsigns = mcfdata->flowrowsigns;
1047  int* colcommodity = mcfdata->colcommodity;
1048  int ncommodities = mcfdata->ncommodities;
1049  int nactivecommodities = mcfdata->ncommodities - mcfdata->nemptycommodities;
1050  SCIP_MCFMODELTYPE modeltype = mcfdata->modeltype;
1051 
1052  unsigned char* capacityrowsigns;
1053  SCIP_Real* capacityrowscores;
1054  int* capacitycands;
1055 
1056  SCIP_ROW** rows;
1057  int nrows;
1058  int r;
1059 
1060  SCIP_Real maxdualcapacity;
1061  int maxcolspercommoditylimit;
1062  int *ncolspercommodity;
1063  int *maxcolspercommodity;
1064  SCIP_Real directedcandsscore;
1065  SCIP_Real undirectedcandsscore;
1066 
1067  SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
1068 
1069  /* allocate temporary memory for extraction data */
1070  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->capacityrowsigns, nrows) ); /*lint !e685*/
1071  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->capacityrowscores, nrows) );
1072  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->capacitycands, nrows) );
1073  capacityrowsigns = mcfdata->capacityrowsigns;
1074  capacityrowscores = mcfdata->capacityrowscores;
1075  capacitycands = mcfdata->capacitycands;
1076 
1077  assert(mcfdata->ncapacitycands == 0);
1078 
1079  /* allocate temporary memory for model type identification */
1080  SCIP_CALL( SCIPallocBufferArray(scip, &ncolspercommodity, ncommodities) );
1081  SCIP_CALL( SCIPallocBufferArray(scip, &maxcolspercommodity, nrows) );
1082 
1083  /* identify model type and set the maximal number of flow variables per capacity constraint and commodity */
1084  switch( modeltype )
1085  {
1087  maxcolspercommoditylimit = 2; /* will be set to 1 later if we detect that the network is directed */
1088  break;
1090  maxcolspercommoditylimit = 1;
1091  break;
1093  maxcolspercommoditylimit = 2;
1094  break;
1095  default:
1096  SCIPerrorMessage("invalid parameter value %d for model type\n", modeltype);
1097  return SCIP_INVALIDDATA;
1098  }
1099 
1100  maxdualcapacity = 0.0;
1101  directedcandsscore = 0.0;
1102  undirectedcandsscore = 0.0;
1103  for( r = 0; r < nrows; r++ )
1104  {
1105  SCIP_ROW* row;
1106  SCIP_COL** rowcols;
1107  SCIP_Real* rowvals;
1108  SCIP_Real rowlhs;
1109  SCIP_Real rowrhs;
1110  int rowlen;
1111  int nposflowcoefs;
1112  int nnegflowcoefs;
1113  int nposcapacitycoefs;
1114  int nnegcapacitycoefs;
1115  int nbadcoefs;
1116  int ncoveredcommodities;
1117  SCIP_Real sameflowcoef;
1118  SCIP_Real sameabsflowcoef;
1119  SCIP_Real maxabscapacitycoef;
1120  SCIP_Real absdualsol;
1121  unsigned char rowsign;
1122  int i;
1123 
1124  row = rows[r];
1125  assert(SCIProwGetLPPos(row) == r);
1126 
1127  capacityrowsigns[r] = 0;
1128  capacityrowscores[r] = 0.0;
1129 
1130  /* ignore modifiable rows */
1131  if( SCIProwIsModifiable(row) )
1132  continue;
1133 
1134  /* ignore empty rows */
1135  rowlen = SCIProwGetNLPNonz(row);
1136  if( rowlen == 0 )
1137  continue;
1138 
1139  /* ignore rows that have already been used as flow conservation constraints */
1140  if( (flowrowsigns[r] & (LHSASSIGNED | RHSASSIGNED)) != 0 )
1141  continue;
1142 
1143  /* get dual solution, if available */
1144  absdualsol = SCIProwGetDualsol(row);
1145  if( absdualsol == SCIP_INVALID ) /*lint !e777*/
1146  absdualsol = 0.0;
1147  absdualsol = ABS(absdualsol);
1148 
1149  rowcols = SCIProwGetCols(row);
1150  rowvals = SCIProwGetVals(row);
1151  rowlhs = SCIProwGetLhs(row) - SCIProwGetConstant(row);
1152  rowrhs = SCIProwGetRhs(row) - SCIProwGetConstant(row);
1153 
1154  /* reset commodity counting array */
1155  BMSclearMemoryArray(ncolspercommodity, ncommodities);
1156  maxcolspercommodity[r] = 0;
1157 
1158  /* identify capacity constraints */
1159  nposflowcoefs = 0;
1160  nnegflowcoefs = 0;
1161  nposcapacitycoefs = 0;
1162  nnegcapacitycoefs = 0;
1163  nbadcoefs = 0;
1164  ncoveredcommodities = 0;
1165  sameflowcoef = 0.0;
1166  sameabsflowcoef = 0.0;
1167  maxabscapacitycoef = 0.0;
1168 
1169  rowsign = 0;
1170  if( !SCIPisInfinity(scip, -rowlhs) )
1171  rowsign |= LHSPOSSIBLE;
1172  if( !SCIPisInfinity(scip, rowrhs) )
1173  rowsign |= RHSPOSSIBLE;
1174  for( i = 0; i < rowlen; i++ )
1175  {
1176  int c;
1177  int k;
1178 
1179  c = SCIPcolGetLPPos(rowcols[i]);
1180  assert(0 <= c && c < SCIPgetNLPCols(scip));
1181 
1182  /* check if this is a flow variable */
1183  k = colcommodity[c];
1184  assert(-1 <= k && k < ncommodities);
1185  if( k >= 0 )
1186  {
1187  SCIP_Real abscoef;
1188 
1189  abscoef = ABS(rowvals[i]);
1190  if( sameflowcoef == 0.0 )
1191  sameflowcoef = rowvals[i];
1192  else if( !SCIPisEQ(scip, sameflowcoef, rowvals[i]) )
1193  sameflowcoef = SCIP_REAL_MAX;
1194  if( sameabsflowcoef == 0.0 )
1195  sameabsflowcoef = abscoef;
1196  else if( !SCIPisEQ(scip, sameabsflowcoef, abscoef) )
1197  sameabsflowcoef = SCIP_REAL_MAX;
1198 
1199  if( rowvals[i] > 0.0 )
1200  nposflowcoefs++;
1201  else
1202  nnegflowcoefs++;
1203 
1204  /* count number of covered commodities in capacity candidate */
1205  if( ncolspercommodity[k] == 0 )
1206  ncoveredcommodities++;
1207  ncolspercommodity[k]++;
1208  maxcolspercommodity[r] = MAX(maxcolspercommodity[r], ncolspercommodity[k]);
1209 
1210  if( ncolspercommodity[k] >= 2 )
1211  capacityrowsigns[r] |= UNDIRECTED;
1212  }
1213  else
1214 /* if( SCIPvarGetType(SCIPcolGetVar(rowcols[i])) != SCIP_VARTYPE_CONTINUOUS ) */
1215  {
1216  SCIP_Real abscoef;
1217 
1218  /* save maximal capacity coef*/
1219  abscoef = ABS(rowvals[i]);
1220  if( abscoef > maxabscapacitycoef )
1221  maxabscapacitycoef = abscoef;
1222 
1223  /* a variable which is not a flow variable can be used as capacity variable */
1224  if( rowvals[i] > 0.0 )
1225  nposcapacitycoefs++;
1226  else
1227  nnegcapacitycoefs++;
1228 
1229  /* a continuous variable is considered to be not so nice*/
1231  nbadcoefs++;
1232  }
1233  }
1234 
1235  /* check if this is a valid capacity constraint */
1236  /* it has at least one flow variable */
1237  if( rowsign != 0 && nposflowcoefs + nnegflowcoefs > 0 )
1238  {
1239  SCIP_Real commodityexcessratio;
1240 
1241  capacityrowsigns[r] |= rowsign;
1242  capacitycands[mcfdata->ncapacitycands] = r;
1243  mcfdata->ncapacitycands++;
1244 
1245  /* calculate capacity row score */
1246  capacityrowscores[r] = 1.0;
1247 
1248  /* calculate mean commodity excess: in the (un)directed case there should be exactly */
1249  /* one (two) flow variable per commodity. in this case commodityexcessratio = 0 */
1250  assert(ncoveredcommodities > 0);
1251  commodityexcessratio =
1252  ABS((nposflowcoefs + nnegflowcoefs)/(SCIP_Real)ncoveredcommodities - maxcolspercommoditylimit);
1253 
1254  capacityrowscores[r] += 1000.0 * MAX(0.0, 2.0 - commodityexcessratio);
1255 
1256  /* row has at most 'maxcolspercommoditylimit' columns per commodity: score +1000 */
1257 /* if( maxcolspercommodity[r] <= maxcolspercommoditylimit )
1258  capacityrowscores[r] += 1000.0;*/
1259 
1260  /* row is of type f - c*x <= b: score +1000 */
1261  if( (capacityrowsigns[r] & RHSPOSSIBLE) != 0 && nnegflowcoefs == 0 && nposcapacitycoefs == 0 && nnegcapacitycoefs > 0 )
1262  capacityrowscores[r] += 1000.0;
1263  if( (capacityrowsigns[r] & LHSPOSSIBLE) != 0 && nposflowcoefs == 0 && nposcapacitycoefs > 0 && nnegcapacitycoefs == 0 )
1264  capacityrowscores[r] += 1000.0;
1265 
1266  /* row has no continuous variables that are not flow variables: score +1000 */
1267 /* if( nbadcoefs == 0 )
1268  capacityrowscores[r] += 1000.0;*/
1269 
1270  /* almost all commodities are covered: score +2000*ncoveredcommodities/(nactivecommodities+3)
1271  * use slightly increased denominator in order to not increase score too much for very few commodities
1272  */
1273  assert(nactivecommodities + 3 > 0);
1274  capacityrowscores[r] += 2000.0 * ncoveredcommodities/(SCIP_Real)(nactivecommodities + 3);
1275 
1276  /* all coefficients of flow variables are +1 or all are -1: score +500 */
1277  if( SCIPisEQ(scip, ABS(sameflowcoef), 1.0) )
1278  capacityrowscores[r] += 500.0;
1279 
1280  /* all coefficients of flow variables are equal: score +250 */
1281  if( sameflowcoef != 0.0 && sameflowcoef != SCIP_REAL_MAX ) /*lint !e777*/
1282  capacityrowscores[r] += 250.0;
1283 
1284  /* all coefficients of flow variables are +1 or -1: score +100 */
1285  if( SCIPisEQ(scip, sameabsflowcoef, 1.0) )
1286  capacityrowscores[r] += 100.0;
1287 
1288  /* there is at least one capacity variable with coefficient not equal to +/-1: score +100 */
1289  if( maxabscapacitycoef > 0.0 && !SCIPisEQ(scip, maxabscapacitycoef, 1.0) )
1290  capacityrowscores[r] += 100.0;
1291 
1292  /* flow coefficients are mostly of the same sign: score +20*max(npos,nneg)/(npos+nneg) */
1293  capacityrowscores[r] += 20.0 * MAX(nposflowcoefs, nnegflowcoefs)/MAX(1.0,(SCIP_Real)(nposflowcoefs + nnegflowcoefs));
1294 
1295  /* capacity coefficients are mostly of the same sign: score +10*max(npos,nneg)/(npos+nneg+1) */
1296  capacityrowscores[r] += 10.0 * MAX(nposcapacitycoefs, nnegcapacitycoefs)/(SCIP_Real)(nposcapacitycoefs+nnegcapacitycoefs+1.0);
1297 
1298  /* row is a <= row with non-negative right hand side: score +10 */
1299  if( (capacityrowsigns[r] & RHSPOSSIBLE) != 0 && !SCIPisNegative(scip, rowrhs) )
1300  capacityrowscores[r] += 10.0;
1301 
1302  /* row is an inequality: score +10 */
1303  if( SCIPisInfinity(scip, -rowlhs) != SCIPisInfinity(scip, rowrhs) )
1304  capacityrowscores[r] += 10.0;
1305 
1306  assert(capacityrowscores[r] > 0.0);
1307  SCIPdebugMsg(scip, "row <%s>: maxcolspercommodity=%d capacityrowsign=%d nposflowcoefs=%d nnegflowcoefs=%d nposcapacitycoefs=%d nnegcapacitycoefs=%d nbadcoefs=%d nactivecommodities=%d sameflowcoef=%g -> score=%g\n",
1308  SCIProwGetName(row), maxcolspercommodity[r], capacityrowsigns[r], nposflowcoefs, nnegflowcoefs, nposcapacitycoefs, nnegcapacitycoefs, nbadcoefs, nactivecommodities, sameflowcoef, capacityrowscores[r]);
1309 
1310  /* update maximum dual solution value for additional score tie breaking */
1311  maxdualcapacity = MAX(maxdualcapacity, absdualsol);
1312 
1313  /* if the model type should be detected automatically, count the number of directed and undirected capacity candidates */
1314  if( modeltype == SCIP_MCFMODELTYPE_AUTO )
1315  {
1316  assert(maxcolspercommoditylimit == 2);
1317  if( (capacityrowsigns[r] & UNDIRECTED) != 0 )
1318  undirectedcandsscore += capacityrowscores[r];
1319  else
1320  directedcandsscore += capacityrowscores[r];
1321  }
1322  }
1323  else
1324  {
1325  SCIPdebugMsg(scip, "row <%s>: rowsign = %d nposflowcoefs = %d nnegflowcoefs = %d -> discard\n",
1326  SCIProwGetName(row), rowsign, nposflowcoefs, nnegflowcoefs);
1327  }
1328  }
1329 
1330  /* if the model type should be detected automatically, decide it by a majority vote */
1331  if( modeltype == SCIP_MCFMODELTYPE_AUTO )
1332  {
1333  if( directedcandsscore > undirectedcandsscore )
1334  modeltype = SCIP_MCFMODELTYPE_DIRECTED;
1335  else
1336  modeltype = SCIP_MCFMODELTYPE_UNDIRECTED;
1337 
1338  MCFdebugMessage("detected model type: %s (%g directed score, %g undirected score)\n",
1339  modeltype == SCIP_MCFMODELTYPE_DIRECTED ? "directed" : "undirected", directedcandsscore, undirectedcandsscore);
1340 
1341  if( modeltype == SCIP_MCFMODELTYPE_DIRECTED )
1342  {
1343  int i;
1344 
1345  /* discard all undirected arcs */
1346  for( i = 0; i < mcfdata->ncapacitycands; i++ )
1347  {
1348  r = capacitycands[i];
1349  assert(0 <= r && r < nrows);
1350  if( (capacityrowsigns[r] & UNDIRECTED) != 0 )
1351  {
1352  /* reduce the score of the undirected row in the directed model */
1353  if( maxcolspercommodity[r] <= maxcolspercommoditylimit )
1354  capacityrowscores[r] -= 1000.0;
1355  }
1356  }
1357  }
1358 
1359  /* record the detected model type */
1360  mcfdata->modeltype = modeltype;
1361  }
1362 
1363  /* apply additional score tie breaking using the dual solutions */
1364  if( SCIPisPositive(scip, maxdualcapacity) )
1365  {
1366  int i;
1367 
1368  for( i = 0; i < mcfdata->ncapacitycands; i++ )
1369  {
1370  SCIP_Real dualsol;
1371 
1372  r = capacitycands[i];
1373  assert(0 <= r && r < nrows);
1374  dualsol = SCIProwGetDualsol(rows[r]);
1375  if( dualsol == SCIP_INVALID ) /*lint !e777*/
1376  dualsol = 0.0;
1377  else if( capacityrowsigns[r] == (LHSPOSSIBLE | RHSPOSSIBLE) )
1378  dualsol = ABS(dualsol);
1379  else if( capacityrowsigns[r] == RHSPOSSIBLE )
1380  dualsol = -dualsol;
1381  assert(maxdualcapacity > 0.0); /*for flexelint*/
1382  capacityrowscores[r] += MAX(dualsol, 0.0)/maxdualcapacity;
1383  assert(capacityrowscores[r] > 0.0);
1384  }
1385  }
1386 
1387  /* sort candidates by score */
1388  SCIPsortInd(mcfdata->capacitycands, compCands, (void*)capacityrowscores, mcfdata->ncapacitycands);
1389 
1390  MCFdebugMessage("capacity candidates [%d]\n", mcfdata->ncapacitycands);
1391 #ifdef SCIP_DEBUG
1392  for( r = 0; r < mcfdata->ncapacitycands; r++ )
1393  {
1394  SCIPdebugMsg(scip, "row %4d [score: %2g]: %s\n", mcfdata->capacitycands[r],
1395  capacityrowscores[mcfdata->capacitycands[r]], SCIProwGetName(rows[mcfdata->capacitycands[r]]));
1396  /*SCIPdebug( SCIP_CALL(SCIPprintRow(scip, rows[mcfdata->capacitycands[r]], NULL)) );*/
1397  }
1398 #endif
1399 
1400  /* free temporary memory */
1401  SCIPfreeBufferArray(scip, &maxcolspercommodity);
1402  SCIPfreeBufferArray(scip, &ncolspercommodity);
1403 
1404  return SCIP_OKAY;
1405 }
1406 
1407 /** creates a new commodity */
1408 static
1410  SCIP* scip, /**< SCIP data structure */
1411  MCFDATA* mcfdata /**< internal MCF extraction data to pass to subroutines */
1412  )
1413 {
1414  /* get memory for commoditysigns array */
1415  assert(mcfdata->ncommodities <= mcfdata->commoditysignssize);
1416  if( mcfdata->ncommodities == mcfdata->commoditysignssize )
1417  {
1418  mcfdata->commoditysignssize = MAX(2*mcfdata->commoditysignssize, mcfdata->ncommodities+1);
1419  SCIP_CALL( SCIPreallocMemoryArray(scip, &mcfdata->commoditysigns, mcfdata->commoditysignssize) );
1420  }
1421  assert(mcfdata->ncommodities < mcfdata->commoditysignssize);
1422 
1423  /* create commodity */
1424  SCIPdebugMsg(scip, "**** creating new commodity %d ****\n", mcfdata->ncommodities);
1425  mcfdata->commoditysigns[mcfdata->ncommodities] = 0;
1426  mcfdata->ncommodities++;
1427 
1428  return SCIP_OKAY;
1429 }
1430 
1431 /** creates a new arc */
1432 static
1434  SCIP* scip, /**< SCIP data structure */
1435  MCFDATA* mcfdata, /**< internal MCF extraction data to pass to subroutines */
1436  int source, /**< source of new arc */
1437  int target, /**< target of new arc */
1438  int* newarcid /**< pointer to store id of new arc */
1439  )
1440 {
1441  assert(source != target );
1442  assert(0 <= source && source < mcfdata->nnodes);
1443  assert(0 <= target && target < mcfdata->nnodes);
1444  assert(newarcid != NULL);
1445 
1446  *newarcid = mcfdata->narcs;
1447 
1448  /* get memory for arrays indexed by arcs */
1449  assert(mcfdata->narcs <= mcfdata->arcarraysize);
1450  if( mcfdata->narcs == mcfdata->arcarraysize )
1451  {
1452  mcfdata->arcarraysize = MAX(2*mcfdata->arcarraysize, mcfdata->narcs+1);
1453  SCIP_CALL( SCIPreallocMemoryArray(scip, &mcfdata->arcsources, mcfdata->arcarraysize) );
1454  SCIP_CALL( SCIPreallocMemoryArray(scip, &mcfdata->arctargets, mcfdata->arcarraysize) );
1455  SCIP_CALL( SCIPreallocMemoryArray(scip, &mcfdata->nextinarcs, mcfdata->arcarraysize) );
1456  SCIP_CALL( SCIPreallocMemoryArray(scip, &mcfdata->nextoutarcs, mcfdata->arcarraysize) );
1457  }
1458  assert(mcfdata->narcs < mcfdata->arcarraysize);
1459 
1460  /* capacityrows is a special case since it is used earlier */
1461  if( mcfdata->capacityrowssize < mcfdata->arcarraysize )
1462  {
1463  mcfdata->capacityrowssize = mcfdata->arcarraysize;
1464  SCIP_CALL( SCIPreallocMemoryArray(scip, &mcfdata->capacityrows, mcfdata->capacityrowssize) );
1465  }
1466  assert(mcfdata->narcs < mcfdata->capacityrowssize);
1467 
1468  /* create new arc */
1469  SCIPdebugMsg(scip, "**** creating new arc %d: %d -> %d ****\n", mcfdata->narcs, source, target);
1470 
1471  mcfdata->arcsources[*newarcid] = source;
1472  mcfdata->arctargets[*newarcid] = target;
1473  mcfdata->nextoutarcs[*newarcid] = mcfdata->firstoutarcs[source];
1474  mcfdata->firstoutarcs[source] = *newarcid;
1475  mcfdata->nextinarcs[*newarcid] = mcfdata->firstinarcs[target];
1476  mcfdata->firstinarcs[target] = *newarcid;
1477  mcfdata->capacityrows[*newarcid] = NULL;
1478 
1479  mcfdata->narcs++;
1480 
1481  return SCIP_OKAY;
1482 }
1483 
1484 /** adds the given flow row and all involved columns to the current commodity */
1485 static
1487  SCIP* scip, /**< SCIP data structure */
1488  MCFDATA* mcfdata, /**< internal MCF extraction data to pass to subroutines */
1489  SCIP_ROW* row, /**< flow row to add to current commodity */
1490  unsigned char rowsign, /**< possible flow row signs to use */
1491  int* comcolids, /**< array of column indices of columns in commodity */
1492  int* ncomcolids /**< pointer to number of columns in commodity */
1493  )
1494 {
1495  unsigned char* flowrowsigns = mcfdata->flowrowsigns;
1496  SCIP_Bool* plusflow = mcfdata->plusflow;
1497  SCIP_Bool* minusflow = mcfdata->minusflow;
1498  int ncommodities = mcfdata->ncommodities;
1499  int* commoditysigns = mcfdata->commoditysigns;
1500  int* colcommodity = mcfdata->colcommodity;
1501  int* rowcommodity = mcfdata->rowcommodity;
1502  int* newcols = mcfdata->newcols;
1503 
1504  SCIP_COL** rowcols;
1505  SCIP_Real* rowvals;
1506  int rowlen;
1507  int rowscale;
1508  SCIP_Bool invertrow;
1509  int r;
1510  int k;
1511  int i;
1512 
1513  assert(comcolids != NULL);
1514  assert(ncomcolids != NULL);
1515 
1516  k = ncommodities-1;
1517  assert(k >= 0);
1518 
1519  r = SCIProwGetLPPos(row);
1520  assert(r >= 0);
1521 
1522  /* check if row has to be inverted */
1523  invertrow = ((rowsign & INVERTED) != 0);
1524  rowsign &= ~INVERTED;
1525 
1526  assert(rowcommodity[r] == -1);
1527  assert((flowrowsigns[r] | rowsign) == flowrowsigns[r]);
1528  assert((rowsign & (LHSPOSSIBLE | RHSPOSSIBLE)) == rowsign);
1529  assert(rowsign != 0);
1530 
1531  /* if the row is only usable as flow row in one direction, we cannot change the sign
1532  * of the whole commodity anymore
1533  */
1534  if( (flowrowsigns[r] & (LHSPOSSIBLE | RHSPOSSIBLE)) != (LHSPOSSIBLE | RHSPOSSIBLE) )
1535  commoditysigns[k] = +1; /* we cannot switch directions */
1536 
1537  /* decide the sign (direction) of the row */
1538  if( rowsign == LHSPOSSIBLE )
1539  rowsign = LHSASSIGNED;
1540  else if( rowsign == RHSPOSSIBLE )
1541  rowsign = RHSASSIGNED;
1542  else
1543  {
1544  SCIP_Real dualsol = SCIProwGetDualsol(row);
1545 
1546  assert(rowsign == (LHSPOSSIBLE | RHSPOSSIBLE));
1547 
1548  /* if we have a valid non-zero dual solution, choose the side which is tight */
1549  if( !SCIPisZero(scip, dualsol) && dualsol != SCIP_INVALID ) /*lint !e777*/
1550  {
1551  if( dualsol > 0.0 )
1552  rowsign = LHSASSIGNED;
1553  else
1554  rowsign = RHSASSIGNED;
1555  }
1556  else
1557  {
1558  SCIP_Real rowlhs = SCIProwGetLhs(row) - SCIProwGetConstant(row);
1559  SCIP_Real rowrhs = SCIProwGetRhs(row) - SCIProwGetConstant(row);
1560 
1561  /* choose row sign such that we get a*x <= -d with d non-negative */
1562  if( rowrhs < 0.0 )
1563  rowsign = RHSASSIGNED;
1564  else if( rowlhs > 0.0 )
1565  rowsign = LHSASSIGNED;
1566  else
1567  rowsign = RHSASSIGNED; /* if we are still undecided, choose rhs */
1568  }
1569  }
1570  if( rowsign == RHSASSIGNED )
1571  rowscale = +1;
1572  else
1573  rowscale = -1;
1574 
1575  /* reintroduce inverted flag */
1576  if( invertrow )
1577  {
1578  rowsign |= INVERTED;
1579  rowscale *= -1;
1580  }
1581  flowrowsigns[r] |= rowsign;
1582 
1583  SCIPdebugMsg(scip, "adding flow row %d <%s> with sign %+d%s to commodity %d [score:%g]\n",
1584  r, SCIProwGetName(row), rowscale, (rowsign & INVERTED) != 0 ? " (inverted)" : "",
1585  k, mcfdata->flowrowscores[r]);
1586  /*SCIPdebug( SCIP_CALL(SCIPprintRow(scip, row, NULL)) );*/
1587 
1588  /* add row to commodity */
1589  rowcommodity[r] = k;
1590  rowcols = SCIProwGetCols(row);
1591  rowvals = SCIProwGetVals(row);
1592  rowlen = SCIProwGetNLPNonz(row);
1593  for( i = 0; i < rowlen; i++ )
1594  {
1595  SCIP_Real val;
1596  int c;
1597 
1598  c = SCIPcolGetLPPos(rowcols[i]);
1599  assert(0 <= c && c < SCIPgetNLPCols(scip));
1600 
1601  /* assign column to commodity */
1602  if( colcommodity[c] == -1 )
1603  {
1604  assert(!plusflow[c]);
1605  assert(!minusflow[c]);
1606  assert(mcfdata->nnewcols < SCIPgetNLPCols(scip));
1607  colcommodity[c] = k;
1608  newcols[mcfdata->nnewcols] = c;
1609  mcfdata->nnewcols++;
1610  comcolids[*ncomcolids] = c;
1611  (*ncomcolids)++;
1612  }
1613  assert(colcommodity[c] == k);
1614 
1615  /* update plusflow/minusflow */
1616  val = rowscale * rowvals[i];
1617  if( val > 0.0 )
1618  {
1619  assert(!plusflow[c]);
1620  plusflow[c] = TRUE;
1621  }
1622  else
1623  {
1624  assert(!minusflow[c]);
1625  minusflow[c] = TRUE;
1626  }
1627  }
1628 }
1629 
1630 /* inverts the lhs/rhs assignment of all rows in the given commodity */
1631 static
1633  SCIP* scip, /**< SCIP data structure */
1634  MCFDATA* mcfdata, /**< internal MCF extraction data to pass to subroutines */
1635  int k, /**< commodity that the flow row should enter */
1636  SCIP_ROW** comrows, /**< flow rows in commodity k */
1637  int ncomrows, /**< number of flow rows (number of nodes) in commodity k */
1638  int* comcolids, /**< column indices of columns in commodity k */
1639  int ncomcolids /**< number of columns in commodity k */
1640  )
1641 {
1642  unsigned char* flowrowsigns = mcfdata->flowrowsigns;
1643  SCIP_Bool* plusflow = mcfdata->plusflow;
1644  SCIP_Bool* minusflow = mcfdata->minusflow;
1645 
1646  int i;
1647 
1648  assert(mcfdata->commoditysigns[k] == 0);
1649  assert(comrows != NULL || ncomrows == 0);
1650  assert(comcolids != NULL);
1651 
1652  /* switch assignments of rows */
1653  for( i = 0; i < ncomrows; i++ )
1654  {
1655  SCIP_ROW* row;
1656  int r;
1657  unsigned char rowsign;
1658 
1659  assert(comrows != NULL);
1660  row = comrows[i];
1661  assert( row != NULL );
1662  r = SCIProwGetLPPos(row);
1663  assert(0 <= r && r < SCIPgetNLPRows(scip));
1664  assert(mcfdata->rowcommodity[r] == k);
1665  assert(!SCIPisInfinity(scip, -SCIProwGetLhs(row)));
1666  assert(!SCIPisInfinity(scip, SCIProwGetRhs(row)));
1667 
1668  rowsign = flowrowsigns[r];
1669  assert((rowsign & (LHSASSIGNED | RHSASSIGNED)) != 0);
1670  assert((rowsign & INVERTED) == 0);
1671 
1672  flowrowsigns[r] &= ~(LHSASSIGNED | RHSASSIGNED);
1673  if( (rowsign & LHSASSIGNED) != 0 )
1674  flowrowsigns[r] |= RHSASSIGNED;
1675  else
1676  flowrowsigns[r] |= LHSASSIGNED;
1677  }
1678 
1679  /* switch plus/minusflow of columns of the given commodity */
1680  for( i = 0; i < ncomcolids; i++ )
1681  {
1682  int c;
1683  SCIP_Bool tmp;
1684 
1685  c = comcolids[i];
1686  assert(0 <= c && c < SCIPgetNLPCols(scip));
1687  assert(mcfdata->colcommodity[c] == k);
1688 
1689  tmp = plusflow[c];
1690  plusflow[c] = minusflow[c];
1691  minusflow[c] = tmp;
1692  }
1693 }
1694 
1695 /** deletes a commodity and removes the flow rows again from the system */
1696 static
1698  SCIP* scip, /**< SCIP data structure */
1699  MCFDATA* mcfdata, /**< internal MCF extraction data to pass to subroutines */
1700  int k, /**< commodity to delete */
1701  SCIP_ROW** comrows, /**< flow rows of the commodity */
1702  int nrows, /**< number of flow rows in the commodity */
1703  int* ndelflowrows, /**< pointer to store number of flow rows in deleted commodity */
1704  int* ndelflowvars /**< pointer to store number of flow vars in deleted commodity */
1705  )
1706 {
1707  unsigned char* flowrowsigns = mcfdata->flowrowsigns;
1708  SCIP_Bool* plusflow = mcfdata->plusflow;
1709  SCIP_Bool* minusflow = mcfdata->minusflow;
1710  int ncommodities = mcfdata->ncommodities;
1711  int* colcommodity = mcfdata->colcommodity;
1712  int* rowcommodity = mcfdata->rowcommodity;
1713 
1714  int n;
1715 
1716  assert(0 <= k && k < ncommodities);
1717 
1718  assert( ndelflowrows != NULL );
1719  assert( ndelflowvars != NULL );
1720 
1721  SCIPdebugMsg(scip, "deleting commodity %d (%d total commodities) with %d flow rows\n", k, ncommodities, nrows);
1722 
1723  *ndelflowrows = 0;
1724  *ndelflowvars = 0;
1725 
1726  for( n = 0; n < nrows; n++ )
1727  {
1728  SCIP_ROW* row;
1729  SCIP_COL** rowcols;
1730  int rowlen;
1731  int r;
1732  int i;
1733 
1734  row = comrows[n];
1735  r = SCIProwGetLPPos(row);
1736  assert(r >= 0);
1737  assert(rowcommodity[r] == k);
1738  assert((flowrowsigns[r] & (LHSASSIGNED | RHSASSIGNED)) != 0);
1739 
1740  SCIPdebugMsg(scip, " -> removing row <%s> from commodity\n", SCIProwGetName(row));
1741 
1742  /* remove the lhs/rhs assignment and the inverted flag */
1743  flowrowsigns[r] &= ~(LHSASSIGNED | RHSASSIGNED | INVERTED);
1744 
1745  /* remove row from commodity */
1746  rowcommodity[r] = -1;
1747  rowcols = SCIProwGetCols(row);
1748  rowlen = SCIProwGetNLPNonz(row);
1749  for( i = 0; i < rowlen; i++ )
1750  {
1751  int c;
1752 
1753  c = SCIPcolGetLPPos(rowcols[i]);
1754  assert(0 <= c && c < SCIPgetNLPCols(scip));
1755 
1756  /* remove column from commodity */
1757  assert(colcommodity[c] == k || colcommodity[c] == -1);
1758  if(colcommodity[c] == k)
1759  (*ndelflowvars)++;
1760  colcommodity[c] = -1;
1761 
1762  /* reset plusflow/minusflow */
1763  plusflow[c] = FALSE;
1764  minusflow[c] = FALSE;
1765  }
1766 
1767  (*ndelflowrows)++;
1768  }
1769 
1770  /* get rid of commodity if it is the last one; otherwise, just leave it
1771  * as an empty commodity which will be discarded later
1772  */
1773  if( k == ncommodities-1 )
1774  mcfdata->ncommodities--;
1775  else
1776  mcfdata->nemptycommodities++;
1777 }
1778 
1779 /** checks whether the given row fits into the given commodity and returns the possible flow row signs */
1780 static
1782  SCIP* scip, /**< SCIP data structure */
1783  MCFDATA* mcfdata, /**< internal MCF extraction data to pass to subroutines */
1784  SCIP_ROW* row, /**< flow row to check */
1785  int k, /**< commodity that the flow row should enter */
1786  unsigned char* rowsign, /**< pointer to store the possible flow row signs */
1787  SCIP_Bool* invertcommodity /**< pointer to store whether the commodity has to be inverted to accommodate the row */
1788  )
1789 {
1790  unsigned char* flowrowsigns = mcfdata->flowrowsigns;
1791  SCIP_Bool* plusflow = mcfdata->plusflow;
1792  SCIP_Bool* minusflow = mcfdata->minusflow;
1793  int* colcommodity = mcfdata->colcommodity;
1794  int* rowcommodity = mcfdata->rowcommodity;
1795  int* commoditysigns = mcfdata->commoditysigns;
1796 
1797  SCIP_COL** rowcols;
1798  SCIP_Real* rowvals;
1799  int rowlen;
1800  unsigned char flowrowsign;
1801  unsigned char invflowrowsign;
1802  int r;
1803  int j;
1804 
1805  assert(invertcommodity != NULL);
1806 
1807  *rowsign = 0;
1808  *invertcommodity = FALSE;
1809 
1810  r = SCIProwGetLPPos(row);
1811  assert(0 <= r && r < SCIPgetNLPRows(scip));
1812 
1813  /* ignore rows that are already used */
1814  if( rowcommodity[r] != -1 )
1815  return;
1816 
1817  /* check if row is an available flow row */
1818  flowrowsign = flowrowsigns[r];
1819  assert((flowrowsign & (LHSPOSSIBLE | RHSPOSSIBLE | DISCARDED)) == flowrowsign);
1820  if( (flowrowsign & DISCARDED) != 0 )
1821  return;
1822  if( (flowrowsign & (LHSPOSSIBLE | RHSPOSSIBLE)) == 0 )
1823  return;
1824  invflowrowsign = flowrowsign;
1825 
1826  /* check whether the row fits w.r.t. the signs of the coefficients */
1827  rowcols = SCIProwGetCols(row);
1828  rowvals = SCIProwGetVals(row);
1829  rowlen = SCIProwGetNLPNonz(row);
1830  for( j = 0; j < rowlen && (flowrowsign != 0 || invflowrowsign != 0); j++ )
1831  {
1832  int rowc;
1833 
1834  rowc = SCIPcolGetLPPos(rowcols[j]);
1835  assert(0 <= rowc && rowc < SCIPgetNLPCols(scip));
1836 
1837  /* check if column already belongs to the same commodity */
1838  if( colcommodity[rowc] == k )
1839  {
1840  /* column only fits if it is not yet present with the same sign */
1841  if( plusflow[rowc] )
1842  {
1843  /* column must not be included with positive sign */
1844  if( rowvals[j] > 0.0 )
1845  {
1846  flowrowsign &= ~RHSPOSSIBLE;
1847  invflowrowsign &= ~LHSPOSSIBLE;
1848  }
1849  else
1850  {
1851  flowrowsign &= ~LHSPOSSIBLE;
1852  invflowrowsign &= ~RHSPOSSIBLE;
1853  }
1854  }
1855  if( minusflow[rowc] )
1856  {
1857  /* column must not be included with negative sign */
1858  if( rowvals[j] > 0.0 )
1859  {
1860  flowrowsign &= ~LHSPOSSIBLE;
1861  invflowrowsign &= ~RHSPOSSIBLE;
1862  }
1863  else
1864  {
1865  flowrowsign &= ~RHSPOSSIBLE;
1866  invflowrowsign &= ~LHSPOSSIBLE;
1867  }
1868  }
1869  }
1870  else if( colcommodity[rowc] != -1 )
1871  {
1872  /* column does not fit if it already belongs to a different commodity */
1873  flowrowsign = 0;
1874  invflowrowsign = 0;
1875  }
1876  }
1877 
1878  if( flowrowsign != 0 )
1879  {
1880  /* flow row fits without inverting anything */
1881  *rowsign = flowrowsign;
1882  *invertcommodity = FALSE;
1883  }
1884  else if( invflowrowsign != 0 )
1885  {
1886  /* this must be an inequality */
1887  assert((flowrowsigns[r] & (LHSPOSSIBLE | RHSPOSSIBLE)) != (LHSPOSSIBLE | RHSPOSSIBLE));
1888 
1889  /* flow row fits only if row or commodity is inverted */
1890  if( commoditysigns == NULL || commoditysigns[k] == 0 )
1891  {
1892  /* commodity can be inverted */
1893  *rowsign = invflowrowsign;
1894  *invertcommodity = TRUE;
1895  }
1896  else
1897  {
1898  /* row has to be inverted */
1899  *rowsign = (invflowrowsign | INVERTED);
1900  *invertcommodity = FALSE;
1901  }
1902  }
1903  else
1904  {
1905  /* we can discard the row, since it can also not be member of a different commodity */
1906  SCIPdebugMsg(scip, " -> discard flow row %d <%s>, comoditysign=%d\n", r, SCIProwGetName(row), commoditysigns[k]);
1907  flowrowsigns[r] |= DISCARDED;
1908  }
1909 }
1910 
1911 /** returns a flow conservation row that fits into the current commodity, or NULL */
1912 static
1914  SCIP* scip, /**< SCIP data structure */
1915  MCFDATA* mcfdata, /**< internal MCF extraction data to pass to subroutines */
1916  SCIP_ROW** nextrow, /**< pointer to store next row */
1917  unsigned char* nextrowsign, /**< pointer to store possible signs of next row */
1918  SCIP_Bool* nextinvertcommodity /**< pointer to store whether current commodity has to be inverted to accommodate the next row */
1919  )
1920 {
1921  SCIP_Real* flowrowscores = mcfdata->flowrowscores;
1922  SCIP_Bool* plusflow = mcfdata->plusflow;
1923  SCIP_Bool* minusflow = mcfdata->minusflow;
1924  int* newcols = mcfdata->newcols;
1925  int ncommodities = mcfdata->ncommodities;
1926 
1927  SCIP_COL** cols;
1928  int k;
1929 
1930  assert(nextrow != NULL);
1931  assert(nextrowsign != NULL);
1932 
1933  *nextrow = NULL;
1934  *nextrowsign = 0;
1935  *nextinvertcommodity = FALSE;
1936 
1937  k = ncommodities-1;
1938 
1939  cols = SCIPgetLPCols(scip);
1940  assert(cols != NULL);
1941 
1942  /* check if there are any columns left in the commodity that have not yet been inspected for incident flow rows */
1943  while( mcfdata->nnewcols > 0 )
1944  {
1945  SCIP_COL* col;
1946  SCIP_ROW** colrows;
1947  int collen;
1948  SCIP_ROW* bestrow;
1949  unsigned char bestrowsign;
1950  SCIP_Bool bestinvertcommodity;
1951  SCIP_Real bestscore;
1952  int c;
1953  int i;
1954 
1955  /* pop next new column from stack */
1956  c = newcols[mcfdata->nnewcols-1];
1957  mcfdata->nnewcols--;
1958  assert(0 <= c && c < SCIPgetNLPCols(scip));
1959 
1960  /* check if this columns already as both signs */
1961  assert(plusflow[c] || minusflow[c]);
1962  if( plusflow[c] && minusflow[c] )
1963  continue;
1964 
1965  /* check whether column is incident to a valid flow row that fits into the current commodity */
1966  bestrow = NULL;
1967  bestrowsign = 0;
1968  bestinvertcommodity = FALSE;
1969  bestscore = 0.0;
1970  col = cols[c];
1971  colrows = SCIPcolGetRows(col);
1972  collen = SCIPcolGetNLPNonz(col);
1973  for( i = 0; i < collen; i++ )
1974  {
1975  SCIP_ROW* row;
1976  unsigned char flowrowsign;
1977  SCIP_Bool invertcommodity;
1978 
1979  row = colrows[i];
1980 
1981  /* check if row fits into the current commodity */
1982  getFlowrowFit(scip, mcfdata, row, k, &flowrowsign, &invertcommodity);
1983 
1984  /* do we have a winner? */
1985  if( flowrowsign != 0 )
1986  {
1987  int r;
1988  SCIP_Real score;
1989 
1990  r = SCIProwGetLPPos(row);
1991  assert(0 <= r && r < SCIPgetNLPRows(scip));
1992  score = flowrowscores[r];
1993  assert(score > 0.0);
1994 
1995  /* If we have to invert the row, this will lead to a negative slack variable in the MIR cut,
1996  * which needs to be substituted in the end. We like to avoid this and therefore reduce the
1997  * score.
1998  */
1999  if( (flowrowsign & INVERTED) != 0 )
2000  score *= 0.75;
2001 
2002  if( score > bestscore )
2003  {
2004  bestrow = row;
2005  bestrowsign = flowrowsign;
2006  bestinvertcommodity = invertcommodity;
2007  bestscore = score;
2008  }
2009  }
2010  }
2011 
2012  /* if there was a valid row for this column, pick the best one
2013  * Note: This is not the overall best row, only the one for the first column that has a valid row.
2014  * However, picking the overall best row seems to be too expensive
2015  */
2016  if( bestrow != NULL )
2017  {
2018  assert(bestscore > 0.0);
2019  assert(bestrowsign != 0);
2020  *nextrow = bestrow;
2021  *nextrowsign = bestrowsign;
2022  *nextinvertcommodity = bestinvertcommodity;
2023  break;
2024  }
2025  }
2026 }
2027 
2028 /** extracts flow conservation rows and puts them into commodities */
2029 static
2031  SCIP* scip, /**< SCIP data structure */
2032  MCFDATA* mcfdata, /**< internal MCF extraction data to pass to subroutines */
2033  SCIP_Real maxflowvarflowrowratio, /**< maximum ratio of flowvars and flowrows */
2034  SCIP_Bool* failed /**< pointer to store whether flowrowflowvarratio exceeded */
2035  )
2036 {
2037  int* flowcands = mcfdata->flowcands;
2038 
2039  SCIP_Bool* plusflow;
2040  SCIP_Bool* minusflow;
2041  int* colcommodity;
2042  int* rowcommodity;
2043 
2044  SCIP_ROW** comrows;
2045  int* ncomnodes;
2046  int* comcolids;
2047  int ncomcolids;
2048  SCIP_ROW** rows;
2049  int nrows;
2050  int ncols;
2051  int maxnnodes;
2052  int nflowrows;
2053  int nflowvars;
2054  int i;
2055  int c;
2056  int r;
2057  int k;
2058 
2059  /* get LP data */
2060  rows = SCIPgetLPRows(scip);
2061  nrows = SCIPgetNLPRows(scip);
2062  ncols = SCIPgetNLPCols(scip);
2063 
2064  assert(failed != NULL);
2065  assert(!*failed);
2066 
2067  /* allocate memory */
2068  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->plusflow, ncols) );
2069  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->minusflow, ncols) );
2070  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->colcommodity, ncols) );
2071  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->rowcommodity, nrows) );
2072  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->newcols, ncols) );
2073  plusflow = mcfdata->plusflow;
2074  minusflow = mcfdata->minusflow;
2075  colcommodity = mcfdata->colcommodity;
2076  rowcommodity = mcfdata->rowcommodity;
2077 
2078  /* allocate temporary memory */
2079  SCIP_CALL( SCIPallocBufferArray(scip, &comrows, nrows) );
2080  SCIP_CALL( SCIPallocBufferArray(scip, &ncomnodes, nrows) );
2081  SCIP_CALL( SCIPallocBufferArray(scip, &comcolids, ncols) );
2082 
2083  /* 3. Extract network structure of flow conservation constraints:
2084  * (a) Initialize plusflow[c] = minusflow[c] = FALSE for all columns c and other local data.
2085  */
2086  BMSclearMemoryArray(plusflow, ncols);
2087  BMSclearMemoryArray(minusflow, ncols);
2088  for( c = 0; c < ncols; c++ )
2089  colcommodity[c] = -1;
2090  for( r = 0; r < nrows; r++ )
2091  rowcommodity[r] = -1;
2092 
2093  assert(flowcands != NULL || mcfdata->nflowcands == 0);
2094 
2095  /* (b) As long as there are flow conservation candidates left:
2096  * (i) Create new commodity and use first flow conservation constraint as new row.
2097  * (ii) Add new row to commodity, update pluscom/minuscom accordingly.
2098  * (iii) For the newly added columns search for an incident flow conservation constraint. Pick the one of highest ranking.
2099  * (iv) If found, set new row to this row and goto (ii).
2100  */
2101  maxnnodes = 0;
2102  nflowrows = 0;
2103  nflowvars = 0;
2104  for( i = 0; i < mcfdata->nflowcands; i++ )
2105  {
2106  SCIP_ROW* newrow;
2107  unsigned char newrowsign;
2108  SCIP_Bool newinvertcommodity;
2109  int nnodes;
2110 
2111  assert(flowcands != NULL);
2112  r = flowcands[i];
2113  assert(0 <= r && r < nrows);
2114  newrow = rows[r];
2115 
2116  /* check if row fits into a new commodity */
2117  getFlowrowFit(scip, mcfdata, newrow, mcfdata->ncommodities, &newrowsign, &newinvertcommodity);
2118  if( newrowsign == 0 )
2119  continue;
2120  assert(!newinvertcommodity);
2121  assert((newrowsign & INVERTED) == 0);
2122 
2123  /* start new commodity */
2124  SCIPdebugMsg(scip, " -------------------start new commodity %d--------------------- \n", mcfdata->ncommodities );
2125  SCIP_CALL( createNewCommodity(scip, mcfdata) );
2126  nnodes = 0;
2127  ncomcolids = 0;
2128 
2129  /* fill commodity with flow conservation constraints */
2130  do
2131  {
2132  /* if next flow row demands an inverting of the commodity, do it now */
2133  if( newinvertcommodity )
2134  invertCommodity(scip, mcfdata, mcfdata->ncommodities-1, comrows, nnodes, comcolids, ncomcolids);
2135 
2136  /* add new row to commodity */
2137  SCIPdebugMsg(scip, " -> add flow row <%s> \n", SCIProwGetName(newrow));
2138  addFlowrowToCommodity(scip, mcfdata, newrow, newrowsign, comcolids, &ncomcolids);
2139  comrows[nnodes] = newrow;
2140  nnodes++;
2141  nflowrows++;
2142 
2143  /* get next row to add */
2144  getNextFlowrow(scip, mcfdata, &newrow, &newrowsign, &newinvertcommodity);
2145  }
2146  while( newrow != NULL );
2147 
2148  ncomnodes[mcfdata->ncommodities-1] = nnodes;
2149  maxnnodes = MAX(maxnnodes, nnodes);
2150  nflowvars += ncomcolids;
2151  SCIPdebugMsg(scip, " -> finished commodity %d: identified %d nodes, maxnnodes=%d\n", mcfdata->ncommodities-1, nnodes, maxnnodes);
2152 
2153  /* if the commodity has too few nodes, or if it has much fewer nodes than the largest commodity, discard it */
2154  if( nnodes < MINNODES || nnodes < MINCOMNODESFRACTION * maxnnodes )
2155  {
2156  int ndelflowrows;
2157  int ndelflowvars;
2158 
2159  deleteCommodity(scip, mcfdata, mcfdata->ncommodities-1, comrows, nnodes, &ndelflowrows, &ndelflowvars);
2160  nflowrows -= ndelflowrows;
2161  nflowvars -= ndelflowvars;
2162  assert(nflowvars >= 0);
2163  assert(nflowrows >= 0);
2164  }
2165  }
2166  /* final cleanup of small commodities */
2167  for( k = 0; k < mcfdata->ncommodities; k++ )
2168  {
2169  assert(ncomnodes[k] >= MINNODES);
2170 
2171  /* if the commodity has much fewer nodes than the largest commodity, discard it */
2172  if( ncomnodes[k] < MINCOMNODESFRACTION * maxnnodes )
2173  {
2174  int nnodes;
2175  int ndelflowrows;
2176  int ndelflowvars;
2177 
2178  nnodes = 0;
2179  for( i = 0; i < mcfdata->nflowcands; i++ )
2180  {
2181  assert(flowcands != NULL);
2182  r = flowcands[i];
2183  if( rowcommodity[r] == k )
2184  {
2185  comrows[nnodes] = rows[r];
2186  nnodes++;
2187 #ifdef NDEBUG
2188  if( nnodes == ncomnodes[k] )
2189  break;
2190 #endif
2191  }
2192  }
2193  assert(nnodes == ncomnodes[k]);
2194  deleteCommodity(scip, mcfdata, k, comrows, nnodes, &ndelflowrows, &ndelflowvars);
2195  nflowrows -= ndelflowrows;
2196  nflowvars -= ndelflowvars;
2197  assert(nflowvars >= 0);
2198  assert(nflowrows >= 0);
2199  }
2200  }
2201 
2202  /* free temporary memory */
2203  SCIPfreeBufferArray(scip, &comcolids);
2204  SCIPfreeBufferArray(scip, &ncomnodes);
2205  SCIPfreeBufferArray(scip, &comrows);
2206 
2207  MCFdebugMessage("identified %d commodities (%d empty) with a maximum of %d nodes and %d flowrows, %d flowvars \n",
2208  mcfdata->ncommodities, mcfdata->nemptycommodities, maxnnodes, nflowrows, nflowvars);
2209 
2210  assert(nflowvars >= 0);
2211  assert(nflowrows >= 0);
2212 
2213  /* do not allow flow system exceeding the flowvarflowrowratio (average node degree)*/
2214  if( nflowrows == 0)
2215  *failed = TRUE;
2216  else if( (SCIP_Real)nflowvars / (SCIP_Real)nflowrows > maxflowvarflowrowratio )
2217  *failed = TRUE;
2218 
2219  return SCIP_OKAY;
2220 }
2221 
2222 /** Arc-Detection -- identifies capacity constraints for the arcs and assigns arc ids to columns and capacity constraints */
2223 static
2225  SCIP* scip, /**< SCIP data structure */
2226  MCFDATA* mcfdata /**< internal MCF extraction data to pass to subroutines */
2227  )
2228 {
2229 
2230  unsigned char* capacityrowsigns = mcfdata->capacityrowsigns;
2231  int* colcommodity = mcfdata->colcommodity;
2232 #ifndef NDEBUG
2233  unsigned char* flowrowsigns = mcfdata->capacityrowsigns;
2234  int* rowcommodity = mcfdata->rowcommodity;
2235 #endif
2236 
2237  int* colarcid;
2238  int* rowarcid;
2239 
2240  SCIP_ROW** rows;
2241  SCIP_COL** cols;
2242  int nrows;
2243  int ncols;
2244 
2245  int r;
2246  int c;
2247  int i;
2248 
2249 #ifndef NDEBUG
2250  SCIP_Real* capacityrowscores = mcfdata->capacityrowscores;
2251 #endif
2252  int *capacitycands = mcfdata->capacitycands;
2253  int ncapacitycands = mcfdata->ncapacitycands;
2254 
2255  assert(mcfdata->narcs == 0);
2256  assert(capacitycands != NULL || ncapacitycands == 0);
2257 
2258  /* get LP data */
2259  SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
2260  SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
2261 
2262  /* allocate temporary memory for extraction data */
2263  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->colarcid, ncols) );
2264  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->rowarcid, nrows) );
2265  colarcid = mcfdata->colarcid;
2266  rowarcid = mcfdata->rowarcid;
2267 
2268  /* initialize arcid arrays */
2269  for( c = 0; c < ncols; c++ )
2270  colarcid[c] = -1;
2271  for( r = 0; r < nrows; r++ )
2272  rowarcid[r] = -1;
2273 
2274  /* -> loop through the list of capacity cands in non-increasing score order */
2275  for( i = 0; i < ncapacitycands; i++ )
2276  {
2277 
2278  SCIP_ROW* capacityrow;
2279  SCIP_COL** rowcols;
2280  int rowlen;
2281  int nassignedflowvars;
2282  int nunassignedflowvars;
2283  int k;
2284 
2285  assert(capacitycands != NULL);
2286  r = capacitycands[i];
2287  assert(0 <= r && r < nrows );
2288  capacityrow = rows[r];
2289 
2290  /* row must be a capacity candidate */
2291  assert((capacityrowsigns[r] & (LHSPOSSIBLE | RHSPOSSIBLE)) != 0);
2292  assert((capacityrowsigns[r] & DISCARDED) == 0);
2293  assert(capacityrowscores[r] > 0.0);
2294 
2295  /* row must not be already assigned */
2296  assert((capacityrowsigns[r] & (LHSASSIGNED | RHSASSIGNED)) == 0);
2297  assert(rowarcid[r] == -1);
2298 
2299  /* row should not be a flow conservation constraint */
2300  assert( rowcommodity[r] == -1 );
2301  assert( (flowrowsigns[r] & (LHSASSIGNED | RHSASSIGNED)) == 0 );
2302 
2303  /* count the number of already assigned and not yet assigned flow variables */
2304  rowcols = SCIProwGetCols(capacityrow);
2305  rowlen = SCIProwGetNLPNonz(capacityrow);
2306  nassignedflowvars = 0;
2307  nunassignedflowvars = 0;
2308  for( k = 0; k < rowlen; k++ )
2309  {
2310  c = SCIPcolGetLPPos(rowcols[k]);
2311  assert(0 <= c && c < ncols);
2312 
2313  assert(-1 <= colcommodity[c] && colcommodity[c] < mcfdata->ncommodities);
2314  assert(-1 <= colarcid[c] && colarcid[c] < mcfdata->narcs);
2315 
2316  /* ignore columns that are not flow variables */
2317  if( colcommodity[c] == -1 )
2318  continue;
2319 
2320  /* check if column is already assigned to an arc */
2321  if( colarcid[c] >= 0 )
2322  nassignedflowvars++;
2323  else
2324  nunassignedflowvars++;
2325  }
2326 
2327  /* Ignore row if all of its flow variables have already been assigned to some other arc.
2328  * Only accept the row as capacity constraint if at least 1/3 of its flow vars are
2329  * not yet assigned to some other arc.
2330  */
2331  if( nunassignedflowvars == 0 || nassignedflowvars >= nunassignedflowvars * 2 )
2332  {
2333  SCIPdebugMsg(scip, "discarding capacity candidate row %d <%s> [score:%g]: %d assigned flowvars, %d unassigned flowvars\n",
2334  r, SCIProwGetName(capacityrow), mcfdata->capacityrowscores[r], nassignedflowvars, nunassignedflowvars);
2335  capacityrowsigns[r] |= DISCARDED;
2336  continue;
2337  }
2338 
2339  /* create new arc -- store capacity row */
2340  assert(mcfdata->narcs <= mcfdata->capacityrowssize);
2341  if( mcfdata->narcs == mcfdata->capacityrowssize )
2342  {
2343  mcfdata->capacityrowssize = MAX(2*mcfdata->capacityrowssize, mcfdata->narcs+1);
2344  SCIP_CALL( SCIPreallocMemoryArray(scip, &mcfdata->capacityrows, mcfdata->capacityrowssize) );
2345  }
2346  assert(mcfdata->narcs < mcfdata->capacityrowssize);
2347  assert(mcfdata->narcs < nrows);
2348 
2349  mcfdata->capacityrows[mcfdata->narcs] = capacityrow;
2350 
2351  /* assign the capacity row to a new arc id */
2352  assert(0 <= r && r < nrows);
2353  rowarcid[r] = mcfdata->narcs;
2354 
2355  /* decide which sign to use */
2356  if( (capacityrowsigns[r] & RHSPOSSIBLE) != 0 )
2357  capacityrowsigns[r] |= RHSASSIGNED;
2358  else
2359  {
2360  assert((capacityrowsigns[r] & LHSPOSSIBLE) != 0);
2361  capacityrowsigns[r] |= LHSASSIGNED;
2362  }
2363 
2364  SCIPdebugMsg(scip, "assigning capacity row %d <%s> with sign %+d to arc %d [score:%g]: %d assigned flowvars, %d unassigned flowvars\n",
2365  r, SCIProwGetName(capacityrow), (capacityrowsigns[r] & RHSASSIGNED) != 0 ? +1 : -1, mcfdata->narcs,
2366  mcfdata->capacityrowscores[r], nassignedflowvars, nunassignedflowvars);
2367 
2368  /* assign all involved flow variables to the new arc id */
2369  for( k = 0; k < rowlen; k++ )
2370  {
2371  int rowc = SCIPcolGetLPPos(rowcols[k]);
2372  assert(0 <= rowc && rowc < ncols);
2373 
2374  /* due to aggregations in preprocessing it may happen that a flow variable appears in multiple capacity constraints;
2375  * in this case, assign it to the first that has been found
2376  */
2377  if( colcommodity[rowc] >= 0 && colarcid[rowc] == -1 )
2378  colarcid[rowc] = mcfdata->narcs;
2379  }
2380 
2381  /* increase number of arcs */
2382  mcfdata->narcs++;
2383  }
2384  return SCIP_OKAY;
2385 } /* END extractcapacities */
2386 
2387 
2388 /** collects all flow columns of all commodities (except the one of the base row) that are incident to the node described by the given flow row */
2389 static
2391  SCIP* scip, /**< SCIP data structure */
2392  MCFDATA* mcfdata, /**< internal MCF extraction data to pass to subroutines */
2393  SCIP_ROW* flowrow, /**< flow conservation constraint that defines the node */
2394  int basecommodity /**< commodity of the base row */
2395  )
2396 {
2397  int* colcommodity = mcfdata->colcommodity;
2398  int* colarcid = mcfdata->colarcid;
2399  int* newcols = mcfdata->newcols;
2400  SCIP_ROW** capacityrows = mcfdata->capacityrows;
2401  SCIP_Bool* colisincident = mcfdata->colisincident;
2402 
2403  SCIP_COL** rowcols;
2404  int rowlen;
2405  int i;
2406 
2407 #ifndef NDEBUG
2408  /* check that the marker array is correctly initialized */
2409  for( i = 0; i < SCIPgetNLPCols(scip); i++ )
2410  assert(!colisincident[i]);
2411 #endif
2412 
2413  /* loop through all flow columns in the flow conservation constraint */
2414  rowcols = SCIProwGetCols(flowrow);
2415  rowlen = SCIProwGetNLPNonz(flowrow);
2416  mcfdata->nnewcols = 0;
2417 
2418  for( i = 0; i < rowlen; i++ )
2419  {
2420  SCIP_COL** capacityrowcols;
2421  int capacityrowlen;
2422  int arcid;
2423  int c;
2424  int j;
2425 
2426 
2427  c = SCIPcolGetLPPos(rowcols[i]);
2428  assert(0 <= c && c < SCIPgetNLPCols(scip));
2429 
2430  /* get arc id of the column in the flow conservation constraint */
2431  arcid = colarcid[c];
2432  if( arcid == -1 )
2433  continue;
2434  assert(arcid < mcfdata->narcs);
2435 
2436  /* collect flow variables in the capacity constraint of this arc */
2437  assert(capacityrows[arcid] != NULL);
2438  capacityrowcols = SCIProwGetCols(capacityrows[arcid]);
2439  capacityrowlen = SCIProwGetNLPNonz(capacityrows[arcid]);
2440 
2441  for( j = 0; j < capacityrowlen; j++ )
2442  {
2443  int caprowc;
2444 
2445  caprowc = SCIPcolGetLPPos(capacityrowcols[j]);
2446  assert(0 <= caprowc && caprowc < SCIPgetNLPCols(scip));
2447 
2448  /* ignore columns that do not belong to a commodity, i.e., are not flow variables */
2449  if( colcommodity[caprowc] == -1 )
2450  {
2451  assert(colarcid[caprowc] == -1);
2452  continue;
2453  }
2454  assert(colarcid[caprowc] <= arcid); /* colarcid < arcid if column belongs to multiple arcs, for example, due to an aggregation in presolving */
2455 
2456  /* ignore columns in the same commodity as the base row */
2457  if( colcommodity[caprowc] == basecommodity )
2458  continue;
2459 
2460  /* if not already done, collect the column */
2461  if( !colisincident[caprowc] )
2462  {
2463  assert(mcfdata->nnewcols < SCIPgetNLPCols(scip));
2464  colisincident[caprowc] = TRUE;
2465  newcols[mcfdata->nnewcols] = caprowc;
2466  mcfdata->nnewcols++;
2467  }
2468  }
2469  }
2470 }
2471 
2472 /** compares given row against a base node flow row and calculates a similarity score;
2473  * score is 0.0 if the rows are incompatible
2474  */
2475 static
2477  SCIP* scip, /**< SCIP data structure */
2478  MCFDATA* mcfdata, /**< internal MCF extraction data to pass to subroutines */
2479  int baserowlen, /**< length of base node flow row */
2480  int* basearcpattern, /**< arc pattern of base node flow row */
2481  int basenposuncap, /**< number of uncapacitated vars in base node flow row with positive coeff*/
2482  int basenneguncap, /**< number of uncapacitated vars in base node flow row with negative coeff*/
2483  SCIP_ROW* row, /**< row to compare against base node flow row */
2484  SCIP_Real* score, /**< pointer to store the similarity score */
2485  SCIP_Bool* invertcommodity /**< pointer to store whether the arcs in the commodity of the row have
2486  * to be inverted for the row to be compatible to the base row */
2487  )
2488 {
2489  unsigned char* flowrowsigns = mcfdata->flowrowsigns;
2490  int* commoditysigns = mcfdata->commoditysigns;
2491  int narcs = mcfdata->narcs;
2492  int* rowcommodity = mcfdata->rowcommodity;
2493  int* colarcid = mcfdata->colarcid;
2494  int* arcpattern = mcfdata->zeroarcarray;
2495  SCIP_MCFMODELTYPE modeltype = mcfdata->modeltype;
2496 
2497  SCIP_COL** rowcols;
2498  SCIP_Real* rowvals;
2499  int nposuncap;
2500  int nneguncap;
2501  int ncols;
2502  int rowlen;
2503  int rowcom;
2504  int rowcomsign;
2505  SCIP_Bool incompatible;
2506  SCIP_Real overlap;
2507  int* overlappingarcs;
2508  int noverlappingarcs;
2509  int r;
2510  int i;
2511 
2512  *score = 0.0;
2513  *invertcommodity = FALSE;
2514 
2515 #ifndef NDEBUG
2516  for( i = 0; i < narcs; i++ )
2517  assert(arcpattern[i] == 0);
2518 #endif
2519 
2520  /* allocate temporary memory */
2521  SCIP_CALL( SCIPallocBufferArray(scip, &overlappingarcs, narcs) );
2522 
2523  r = SCIProwGetLPPos(row);
2524  assert(0 <= r && r < SCIPgetNLPRows(scip));
2525  assert((flowrowsigns[r] & (LHSASSIGNED | RHSASSIGNED)) != 0);
2526  rowcom = rowcommodity[r];
2527  assert(0 <= rowcom && rowcom < mcfdata->ncommodities);
2528  rowcomsign = commoditysigns[rowcom];
2529  assert(-1 <= rowcomsign && rowcomsign <= +1);
2530 
2531  rowcols = SCIProwGetCols(row);
2532  rowvals = SCIProwGetVals(row);
2533  rowlen = SCIProwGetNLPNonz(row);
2534  incompatible = FALSE;
2535  noverlappingarcs = 0;
2536  nposuncap=0;
2537  nneguncap=0;
2538  ncols = SCIPgetNLPCols(scip);
2539  for( i = 0; i < rowlen; i++ )
2540  {
2541  int c;
2542  int arcid;
2543  int valsign;
2544 
2545  c = SCIPcolGetLPPos(rowcols[i]);
2546  assert(0 <= c && c < SCIPgetNLPCols(scip));
2547 
2548  /* get the sign of the coefficient in the flow conservation constraint */
2549  valsign = (rowvals[i] > 0.0 ? +1 : -1);
2550  if( (flowrowsigns[r] & LHSASSIGNED) != 0 )
2551  valsign *= -1;
2552  if( (flowrowsigns[r] & INVERTED) != 0 )
2553  valsign *= -1;
2554 
2555  arcid = colarcid[c];
2556  if( arcid == -1 )
2557  {
2558  if( valsign > 0 )
2559  nposuncap++;
2560  else
2561  nneguncap++;
2562  continue;
2563  }
2564  assert(arcid < narcs);
2565 
2566  /* check if this arc is also member of the base row */
2567  if( basearcpattern[arcid] != 0 )
2568  {
2569  /* check if the sign of the arc matches in the directed case */
2570  if( modeltype == SCIP_MCFMODELTYPE_DIRECTED )
2571  {
2572  int validcomsign;
2573 
2574  if( ( valsign * basearcpattern[arcid] ) > 0 )
2575  validcomsign = +1;
2576  else
2577  validcomsign = -1;
2578 
2579  if( rowcomsign == 0 )
2580  {
2581  /* the first entry decides whether we have to invert the commodity */
2582  rowcomsign = validcomsign;
2583  }
2584  else if( rowcomsign != validcomsign )
2585  {
2586  /* the signs do not fit: this is incompatible */
2587  incompatible = TRUE;
2588  break;
2589  }
2590  }
2591  else
2592  {
2593  /* in the undirected case, we ignore the sign of the coefficient */
2594  valsign = +1;
2595  }
2596 
2597  /* store overlapping arc pattern */
2598  if( arcpattern[arcid] == 0 )
2599  {
2600  overlappingarcs[noverlappingarcs] = arcid;
2601  noverlappingarcs++;
2602  }
2603  arcpattern[arcid] += valsign;
2604  }
2605  }
2606 
2607  /* calculate the weighted overlap and reset the zeroarcarray */
2608  overlap = 0.0;
2609  for( i = 0; i < noverlappingarcs; i++ )
2610  {
2611  SCIP_Real basenum;
2612  SCIP_Real arcnum;
2613  int arcid;
2614 
2615  arcid = overlappingarcs[i];
2616  assert(0 <= arcid && arcid < narcs);
2617  assert(modeltype == SCIP_MCFMODELTYPE_UNDIRECTED || rowcomsign * basearcpattern[arcid] * arcpattern[arcid] > 0);
2618 
2619  basenum = ABS(basearcpattern[arcid]);
2620  arcnum = ABS(arcpattern[arcid]);
2621  assert(basenum != 0.0);
2622  assert(arcnum != 0.0);
2623 
2624  if( basenum > arcnum )
2625  overlap += arcnum/basenum;
2626  else
2627  overlap += basenum/arcnum;
2628 
2629  arcpattern[arcid] = 0;
2630  }
2631 
2632 /* calculate the score: maximize overlap and use minimal number of non-overlapping entries as tie breaker */
2633  if( !incompatible && overlap > 0.0 )
2634  {
2635  /* flow variables with arc-id */
2636  int rowarcs = rowlen - nposuncap - nneguncap;
2637  int baserowarcs = baserowlen - basenposuncap - basenneguncap;
2638 
2639  assert(overlap <= (SCIP_Real) rowlen);
2640  assert(overlap <= (SCIP_Real) baserowlen);
2641  assert(noverlappingarcs >= 1);
2642 
2643  *invertcommodity = (rowcomsign == -1);
2644 
2645  /* only one overlapping arc is very dangerous,
2646  since this can also be the other end node of the arc */
2647  if( noverlappingarcs >= 2 )
2648  *score += 1000.0;
2649 
2650  assert(rowarcs >= 0 && baserowarcs >= 0 );
2651  /* in the ideal undirected case there are two flow variables with the same arc-id */
2652  if( modeltype == SCIP_MCFMODELTYPE_DIRECTED )
2653  *score = overlap - (rowarcs + baserowarcs - 2.0 * overlap)/(2.0 * ncols + 1.0);
2654  else
2655  *score = overlap - (rowarcs + baserowarcs - 4.0 * overlap)/(2.0 * ncols + 1.0);
2656 
2657  /* Also use number of uncapacitated flowvars (variables without arcid) as tie-breaker */
2658  if(*invertcommodity)
2659  *score += 1.0 - (ABS(nneguncap - basenposuncap) + ABS(nposuncap - basenneguncap))/(2.0 * ncols + 1.0);
2660  else
2661  *score += 1.0 - (ABS(nposuncap - basenposuncap) + ABS(nneguncap - basenneguncap))/(2.0 * ncols + 1.0);
2662 
2663  *score = MAX(*score, 1e-6); /* score may get negative due to many columns having the same arcid */
2664 
2665  }
2666 
2667  SCIPdebugMsg(scip, " -> node similarity: row <%s>: incompatible=%u overlap=%g rowlen=%d baserowlen=%d score=%g\n",
2668  SCIProwGetName(row), incompatible, overlap, rowlen, baserowlen, *score);
2669 
2670  /* free temporary memory */
2671  SCIPfreeBufferArray(scip, &overlappingarcs);
2672 
2673  return SCIP_OKAY;
2674 }
2675 
2676 /** assigns node ids to flow conservation constraints */
2677 static
2679  SCIP* scip, /**< SCIP data structure */
2680  MCFDATA* mcfdata /**< internal MCF extraction data to pass to subroutines */
2681  )
2682 {
2683  unsigned char* flowrowsigns = mcfdata->flowrowsigns;
2684  int ncommodities = mcfdata->ncommodities;
2685  int* commoditysigns = mcfdata->commoditysigns;
2686  int narcs = mcfdata->narcs;
2687  int* flowcands = mcfdata->flowcands;
2688  int nflowcands = mcfdata->nflowcands;
2689  int* rowcommodity = mcfdata->rowcommodity;
2690  int* colarcid = mcfdata->colarcid;
2691  int* newcols = mcfdata->newcols;
2692  SCIP_MCFMODELTYPE modeltype = mcfdata->modeltype;
2693  int* rownodeid;
2694  SCIP_Bool* colisincident;
2695  SCIP_Bool* rowprocessed;
2696 
2697  SCIP_ROW** rows;
2698  SCIP_COL** cols;
2699  int nrows;
2700  int ncols;
2701 
2702  int* arcpattern;
2703  int nposuncap;
2704  int nneguncap;
2705  SCIP_ROW** bestflowrows;
2706  SCIP_Real* bestscores;
2707  SCIP_Bool* bestinverted;
2708  int r;
2709  int c;
2710  int n;
2711 
2712  assert(mcfdata->nnodes == 0);
2713  assert(modeltype != SCIP_MCFMODELTYPE_AUTO);
2714 
2715  /* get LP data */
2716  SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
2717  SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
2718 
2719  /* allocate temporary memory */
2720  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->rownodeid, nrows) );
2721  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->colisincident, ncols) );
2722  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->zeroarcarray, narcs) );
2723  BMSclearMemoryArray(mcfdata->zeroarcarray, narcs);
2724  rownodeid = mcfdata->rownodeid;
2725  colisincident = mcfdata->colisincident;
2726 
2727  /* allocate temporary local memory */
2728  SCIP_CALL( SCIPallocBufferArray(scip, &arcpattern, narcs) );
2729  SCIP_CALL( SCIPallocBufferArray(scip, &bestflowrows, ncommodities) );
2730  SCIP_CALL( SCIPallocBufferArray(scip, &bestscores, ncommodities) );
2731  SCIP_CALL( SCIPallocBufferArray(scip, &bestinverted, ncommodities) );
2732  SCIP_CALL( SCIPallocBufferArray(scip, &rowprocessed, nrows) );
2733 
2734  /* initialize temporary memory */
2735  for( r = 0; r < nrows; r++ )
2736  rownodeid[r] = -1;
2737  for( c = 0; c < ncols; c++ )
2738  colisincident[c] = FALSE;
2739 
2740  assert(flowcands != NULL || nflowcands == 0);
2741 
2742  /* process all flow conservation constraints that have been used */
2743  for( n = 0; n < nflowcands; n++ )
2744  {
2745  SCIP_COL** rowcols;
2746  SCIP_Real* rowvals;
2747  int rowlen;
2748  int rowscale;
2749  int basecommodity;
2750  int i;
2751 
2752  assert(flowcands != NULL);
2753  r = flowcands[n];
2754  assert(0 <= r && r < nrows);
2755 
2756  /* ignore rows that are not used as flow conservation constraint */
2757  basecommodity = rowcommodity[r];
2758  if( basecommodity == -1 )
2759  continue;
2760  assert((flowrowsigns[r] & (LHSASSIGNED | RHSASSIGNED)) != 0);
2761  assert(mcfdata->rowarcid[r] == -1);
2762 
2763  /* skip rows that are already assigned to a node */
2764  if( rownodeid[r] >= 0 )
2765  continue;
2766 
2767  /* assign row to new node id */
2768  SCIPdebugMsg(scip, "assigning row %d <%s> of commodity %d to node %d [score: %g]\n",
2769  r, SCIProwGetName(rows[r]), basecommodity, mcfdata->nnodes, mcfdata->flowrowscores[r]);
2770  rownodeid[r] = mcfdata->nnodes;
2771 
2772  /* increase number of nodes */
2773  mcfdata->nnodes++;
2774 
2775  /* For single commodity models we are done --
2776  * no matching flow rows need to be found
2777  */
2778  if(ncommodities == 1)
2779  continue;
2780 
2781  /* get the arc pattern of the flow row */
2782  BMSclearMemoryArray(arcpattern, narcs);
2783  nposuncap=0;
2784  nneguncap=0;
2785 
2786  rowcols = SCIProwGetCols(rows[r]);
2787  rowvals = SCIProwGetVals(rows[r]);
2788  rowlen = SCIProwGetNLPNonz(rows[r]);
2789  if( (flowrowsigns[r] & RHSASSIGNED) != 0 )
2790  rowscale = +1;
2791  else
2792  rowscale = -1;
2793  if( (flowrowsigns[r] & INVERTED) != 0 )
2794  rowscale *= -1;
2795  if( commoditysigns[basecommodity] == -1 )
2796  rowscale *= -1;
2797 
2798  for( i = 0; i < rowlen; i++ )
2799  {
2800  int arcid;
2801 
2802  c = SCIPcolGetLPPos(rowcols[i]);
2803  assert(0 <= c && c < ncols);
2804  arcid = colarcid[c];
2805  if( arcid >= 0 )
2806  {
2807  /* due to presolving we may have multiple flow variables of the same arc in the row */
2808  if( modeltype == SCIP_MCFMODELTYPE_UNDIRECTED || rowscale * rowvals[i] > 0.0 )
2809  arcpattern[arcid]++;
2810  else
2811  arcpattern[arcid]--;
2812  }
2813  /* we also count variables that have no arc -- these have no capacity constraint --> uncapacitated */
2814  else
2815  {
2816  if( modeltype == SCIP_MCFMODELTYPE_UNDIRECTED || rowscale * rowvals[i] > 0.0 )
2817  nposuncap++;
2818  else
2819  nneguncap++;
2820  }
2821  }
2822 
2823  /* initialize arrays to store best flow rows */
2824  for( i = 0; i < ncommodities; i++ )
2825  {
2826  bestflowrows[i] = NULL;
2827  bestscores[i] = 0.0;
2828  bestinverted[i] = FALSE;
2829  }
2830 
2831  /* collect columns that are member of incident arc capacity constraints */
2832  collectIncidentFlowCols(scip, mcfdata, rows[r], basecommodity);
2833 
2834  /* initialize rowprocessed array */
2835  BMSclearMemoryArray(rowprocessed, nrows);
2836 
2837  /* identify flow conservation constraints in other commodities that match this node;
2838  * search for flow rows in the column vectors of the incident columns
2839  */
2840  for( i = 0; i < mcfdata->nnewcols; i++ )
2841  {
2842  SCIP_ROW** colrows;
2843  int collen;
2844  int j;
2845 
2846  c = newcols[i];
2847  assert(0 <= c && c < ncols);
2848  assert(mcfdata->colcommodity[c] >= 0);
2849  assert(mcfdata->colcommodity[c] != basecommodity);
2850 
2851  /* clean up the marker array */
2852  assert(colisincident[c]);
2853  colisincident[c] = FALSE;
2854 
2855  /* scan column vector for flow conservation constraints */
2856  colrows = SCIPcolGetRows(cols[c]);
2857  collen = SCIPcolGetNLPNonz(cols[c]);
2858 
2859  for( j = 0; j < collen; j++ )
2860  {
2861  int colr;
2862  int rowcom;
2863  SCIP_Real score;
2864  SCIP_Bool invertcommodity;
2865 
2866  colr = SCIProwGetLPPos(colrows[j]);
2867  assert(0 <= colr && colr < nrows);
2868 
2869  /* ignore rows that have already been processed */
2870  if( rowprocessed[colr] )
2871  continue;
2872  rowprocessed[colr] = TRUE;
2873 
2874  /* ignore rows that are not flow conservation constraints in the network */
2875  rowcom = rowcommodity[colr];
2876  assert(rowcom != basecommodity);
2877  if( rowcom == -1 )
2878  continue;
2879 
2880  assert(rowcom == mcfdata->colcommodity[c]);
2881  assert((flowrowsigns[colr] & (LHSASSIGNED | RHSASSIGNED)) != 0);
2882  assert(mcfdata->rowarcid[colr] == -1);
2883 
2884  /* ignore rows that are already assigned to a node */
2885  if( rownodeid[colr] >= 0 )
2886  continue;
2887 
2888  /* compare row against arc pattern and calculate score */
2889  SCIP_CALL( getNodeSimilarityScore(scip, mcfdata, rowlen, arcpattern,
2890  nposuncap, nneguncap, colrows[j], &score, &invertcommodity) );
2891  assert( !SCIPisNegative(scip, score) );
2892 
2893  if( score > bestscores[rowcom] )
2894  {
2895  bestflowrows[rowcom] = colrows[j];
2896  bestscores[rowcom] = score;
2897  bestinverted[rowcom] = invertcommodity;
2898  }
2899  }
2900  }
2901  assert(bestflowrows[basecommodity] == NULL);
2902 
2903  /* for each commodity, pick the best flow conservation constraint to define this node */
2904  for( i = 0; i < ncommodities; i++ )
2905  {
2906  int comr;
2907 
2908  if( bestflowrows[i] == NULL )
2909  continue;
2910 
2911  comr = SCIProwGetLPPos(bestflowrows[i]);
2912  assert(0 <= comr && comr < nrows);
2913  assert(rowcommodity[comr] == i);
2914  assert((flowrowsigns[comr] & (LHSASSIGNED | RHSASSIGNED)) != 0);
2915  assert(rownodeid[comr] == -1);
2916  assert(mcfdata->nnodes >= 1);
2917  /* assign flow row to current node */
2918  SCIPdebugMsg(scip, " -> assigning row %d <%s> of commodity %d to node %d [invert:%u]\n",
2919  comr, SCIProwGetName(rows[comr]), i, mcfdata->nnodes-1, bestinverted[i]);
2920  rownodeid[comr] = mcfdata->nnodes-1;
2921 
2922  /* fix the direction of the arcs of the commodity */
2923  if( bestinverted[i] )
2924  {
2925  assert(commoditysigns[i] != +1);
2926  commoditysigns[i] = -1;
2927  }
2928  else
2929  {
2930  assert(commoditysigns[i] != -1);
2931  commoditysigns[i] = +1;
2932  }
2933  }
2934  }
2935 
2936  /* free local temporary memory */
2937 
2938  SCIPfreeBufferArray(scip, &rowprocessed);
2939  SCIPfreeBufferArray(scip, &bestinverted);
2940  SCIPfreeBufferArray(scip, &bestscores);
2941  SCIPfreeBufferArray(scip, &bestflowrows);
2942  SCIPfreeBufferArray(scip, &arcpattern);
2943 
2944  return SCIP_OKAY;
2945 }
2946 
2947 /** if there are still undecided commodity signs, fix them to +1 */
2948 static
2950  SCIP* scip, /**< SCIP data structure */
2951  MCFDATA* mcfdata /**< internal MCF extraction data to pass to subroutines */
2952  )
2953 {
2954  int* commoditysigns = mcfdata->commoditysigns;
2955  int k;
2956 
2957  for( k = 0; k < mcfdata->ncommodities; k++ )
2958  {
2959  if( commoditysigns[k] == 0 )
2960  commoditysigns[k] = +1;
2961  }
2962 }
2963 
2964 
2965 /** identifies the (at most) two nodes which contain the given flow variable */
2966 static
2968  SCIP* scip, /**< SCIP data structure */
2969  MCFDATA* mcfdata, /**< internal MCF extraction data to pass to subroutines */
2970  SCIP_COL* col, /**< flow column */
2971  int* sourcenode, /**< pointer to store the source node of the flow column */
2972  int* targetnode /**< pointer to store the target node of the flow column */
2973  )
2974 {
2975  unsigned char* flowrowsigns = mcfdata->flowrowsigns;
2976  int* commoditysigns = mcfdata->commoditysigns;
2977  int* rowcommodity = mcfdata->rowcommodity;
2978  int* rownodeid = mcfdata->rownodeid;
2979  int* colsources = mcfdata->colsources;
2980  int* coltargets = mcfdata->coltargets;
2981 
2982  SCIP_ROW** colrows;
2983  SCIP_Real* colvals;
2984  int collen;
2985  int c;
2986  int i;
2987 
2988  assert(sourcenode != NULL);
2989  assert(targetnode != NULL);
2990  assert(colsources != NULL);
2991  assert(coltargets != NULL);
2992 
2993  c = SCIPcolGetLPPos(col);
2994  assert(0 <= c && c < SCIPgetNLPCols(scip));
2995 
2996  /* check if we have this column already in cache */
2997  if( colsources[c] >= -1 )
2998  {
2999  assert(coltargets[c] >= -1);
3000  *sourcenode = colsources[c];
3001  *targetnode = coltargets[c];
3002  }
3003  else
3004  {
3005  *sourcenode = -1;
3006  *targetnode = -1;
3007 
3008  /* search for flow conservation rows in the column vector */
3009  colrows = SCIPcolGetRows(col);
3010  colvals = SCIPcolGetVals(col);
3011  collen = SCIPcolGetNLPNonz(col);
3012  for( i = 0; i < collen; i++ )
3013  {
3014  int r;
3015 
3016  r = SCIProwGetLPPos(colrows[i]);
3017  assert(0 <= r && r < SCIPgetNLPRows(scip));
3018 
3019  if( rownodeid[r] >= 0 )
3020  {
3021  int v;
3022  int k;
3023  int scale;
3024 
3025  v = rownodeid[r];
3026  k = rowcommodity[r];
3027  assert(0 <= v && v < mcfdata->nnodes);
3028  assert(0 <= k && k < mcfdata->ncommodities);
3029  assert((flowrowsigns[r] & (LHSASSIGNED | RHSASSIGNED)) != 0);
3030 
3031  /* check whether the flow row is inverted */
3032  scale = +1;
3033  if( (flowrowsigns[r] & LHSASSIGNED) != 0 )
3034  scale *= -1;
3035  if( (flowrowsigns[r] & INVERTED) != 0 )
3036  scale *= -1;
3037  if( commoditysigns[k] == -1 )
3038  scale *= -1;
3039 
3040  /* decide whether this node is source or target */
3041  if( ( scale * colvals[i] ) > 0.0 )
3042  {
3043  assert(*sourcenode == -1);
3044  *sourcenode = v;
3045  if( *targetnode >= 0 )
3046  break;
3047  }
3048  else
3049  {
3050  assert(*targetnode == -1);
3051  *targetnode = v;
3052  if( *sourcenode >= 0 )
3053  break;
3054  }
3055  }
3056  }
3057 
3058  /* cache result for future use */
3059  colsources[c] = *sourcenode;
3060  coltargets[c] = *targetnode;
3061  }
3062 }
3063 
3064 /** find uncapacitated arcs for flow columns that have no associated arc yet */
3065 static
3067  SCIP* scip, /**< SCIP data structure */
3068  MCFDATA* mcfdata /**< internal MCF extraction data to pass to subroutines */
3069  )
3070 {
3071  int* flowcands = mcfdata->flowcands;
3072  int nflowcands = mcfdata->nflowcands;
3073 #ifndef NDEBUG
3074  unsigned char* flowrowsigns = mcfdata->flowrowsigns;
3075  int* colcommodity = mcfdata->colcommodity;
3076  int* rowcommodity = mcfdata->rowcommodity;
3077 #endif
3078  int* rownodeid = mcfdata->rownodeid;
3079  int* colarcid = mcfdata->colarcid;
3080  int nnodes = mcfdata->nnodes;
3081  int ncommodities = mcfdata->ncommodities;
3082  SCIP_MCFMODELTYPE modeltype = mcfdata->modeltype;
3083 
3084  SCIP_ROW** rows;
3085  SCIP_COL** cols;
3086  int ncols;
3087 
3088  int* sortedflowcands;
3089  int* sortedflowcandnodeid;
3090  int* sourcecount;
3091  int* targetcount;
3092  int* adjnodes;
3093  int nadjnodes;
3094  int* inccols;
3095  int ninccols;
3096  int arcsthreshold;
3097 
3098  int v;
3099  int n;
3100 
3101  /* there should have been a cleanup already */
3102  assert(mcfdata->nemptycommodities == 0);
3103  assert(ncommodities >= 0);
3104  assert(modeltype == SCIP_MCFMODELTYPE_UNDIRECTED || modeltype == SCIP_MCFMODELTYPE_DIRECTED);
3105 
3106  /* avoid trivial cases */
3107  if( ncommodities == 0 || nflowcands == 0 || nnodes == 0 )
3108  return SCIP_OKAY;
3109 
3110  SCIPdebugMsg(scip, "finding uncapacitated arcs\n");
3111 
3112  /* get LP data */
3113  rows = SCIPgetLPRows(scip);
3114  cols = SCIPgetLPCols(scip);
3115  ncols = SCIPgetNLPCols(scip);
3116  assert(rows != NULL);
3117  assert(cols != NULL || ncols == 0);
3118 
3119  /* allocate temporary memory */
3120  SCIP_CALL( SCIPallocBufferArray(scip, &sortedflowcands, nflowcands) );
3121  SCIP_CALL( SCIPallocBufferArray(scip, &sortedflowcandnodeid, nflowcands) );
3122  SCIP_CALL( SCIPallocBufferArray(scip, &sourcecount, nnodes) );
3123  SCIP_CALL( SCIPallocBufferArray(scip, &targetcount, nnodes) );
3124  SCIP_CALL( SCIPallocBufferArray(scip, &adjnodes, nnodes) );
3125  SCIP_CALL( SCIPallocBufferArray(scip, &inccols, ncols) );
3126 
3127  /* copy flowcands and initialize sortedflowcandnodeid arrays */
3128  for( n = 0; n < nflowcands; n++ )
3129  {
3130  sortedflowcands[n] = flowcands[n];
3131  sortedflowcandnodeid[n] = rownodeid[flowcands[n]];
3132  }
3133 
3134  /* sort flow candidates by node id */
3135  SCIPsortIntInt(sortedflowcandnodeid, sortedflowcands, nflowcands);
3136  assert(sortedflowcandnodeid[0] <= 0);
3137  assert(sortedflowcandnodeid[nflowcands-1] == nnodes-1);
3138 
3139  /* initialize sourcecount and targetcount arrays */
3140  for( v = 0; v < nnodes; v++ )
3141  {
3142  sourcecount[v] = 0;
3143  targetcount[v] = 0;
3144  }
3145  nadjnodes = 0;
3146  ninccols = 0;
3147 
3148  /* we only accept an arc if at least this many flow variables give rise to this arc */
3149  arcsthreshold = (int) SCIPceil(scip, (SCIP_Real) ncommodities * UNCAPACITATEDARCSTRESHOLD );
3150 
3151  /* in the undirected case, there are two variables per commodity in each capacity row */
3152  if( modeltype == SCIP_MCFMODELTYPE_UNDIRECTED )
3153  arcsthreshold *= 2;
3154 
3155  /* skip unused flow candidates */
3156  for( n = 0; n < nflowcands; n++ )
3157  {
3158  if( sortedflowcandnodeid[n] >= 0 )
3159  break;
3160  assert(0 <= sortedflowcands[n] && sortedflowcands[n] < SCIPgetNLPRows(scip));
3161  assert(rowcommodity[sortedflowcands[n]] == -1);
3162  }
3163  assert(n < nflowcands);
3164  assert(sortedflowcandnodeid[n] == 0);
3165 
3166  /* for each node s, count for each other node t the number of flow variables that are not yet assigned
3167  * to an arc and that give rise to an (s,t) arc or an (t,s) arc
3168  */
3169  for( v = 0; n < nflowcands; v++ ) /*lint !e440*/ /* for flexelint: n is used as abort criterion for loop */
3170  {
3171  int l;
3172 
3173  assert(v < nnodes);
3174  assert(0 <= sortedflowcands[n] && sortedflowcands[n] < SCIPgetNLPRows(scip));
3175  assert(rowcommodity[sortedflowcands[n]] >= 0);
3176  assert(rownodeid[sortedflowcands[n]] == sortedflowcandnodeid[n]);
3177  assert(sortedflowcandnodeid[n] == v); /* we must have at least one row per node */
3178  assert(nadjnodes == 0);
3179  assert(ninccols == 0);
3180 
3181  SCIPdebugMsg(scip, " node %d starts with flowcand %d: <%s>\n", v, n, SCIProwGetName(rows[sortedflowcands[n]]));
3182 
3183  /* process all flow rows that belong to node v */
3184  for( ; n < nflowcands && sortedflowcandnodeid[n] == v; n++ )
3185  {
3186  SCIP_COL** rowcols;
3187  int rowlen;
3188  int r;
3189  int i;
3190 
3191  r = sortedflowcands[n];
3192  assert((flowrowsigns[r] & (LHSASSIGNED | RHSASSIGNED)) != 0);
3193  assert(mcfdata->rowarcid[r] == -1);
3194 
3195  /* update sourcecount and targetcount for all flow columns in the row that are not yet assigned to an arc */
3196  rowcols = SCIProwGetCols(rows[r]);
3197  rowlen = SCIProwGetNLPNonz(rows[r]);
3198  for( i = 0; i < rowlen; i++ )
3199  {
3200  SCIP_COL* col;
3201  int arcid;
3202  int c;
3203  int s;
3204  int t;
3205 
3206  col = rowcols[i];
3207  c = SCIPcolGetLPPos(col);
3208  assert(0 <= c && c < SCIPgetNLPCols(scip));
3209  arcid = colarcid[c];
3210  assert(-2 <= arcid && arcid < mcfdata->narcs);
3211  assert(rowcommodity[r] == colcommodity[c]);
3212 
3213  if( arcid == -2 )
3214  {
3215  /* This is the second time we see this column, and we were unable to assign an arc
3216  * to this column at the first time. So, this time we can ignore it. Just reset the
3217  * temporary arcid -2 to -1.
3218  */
3219  colarcid[c] = -1;
3220  }
3221  else if( arcid == -1 )
3222  {
3223  int u;
3224 
3225  /* identify the (at most) two nodes which contain this flow variable */
3226  getIncidentNodes(scip, mcfdata, col, &s, &t);
3227 
3228  SCIPdebugMsg(scip, " col <%s> [%g,%g] (s,t):(%i,%i)\n", SCIPvarGetName(SCIPcolGetVar(col)),
3230 
3231  assert(-1 <= s && s < nnodes);
3232  assert(-1 <= t && t < nnodes);
3233  assert(s == v || t == v);
3234  assert(s != t);
3235 
3236  /* in the undirected case, always use s as other node */
3237  if( modeltype == SCIP_MCFMODELTYPE_UNDIRECTED && s == v )
3238  {
3239  s = t;
3240  t = v;
3241  }
3242 
3243  /* if there is no other node than v, ignore column */
3244  if( s < 0 || t < 0 )
3245  continue;
3246 
3247  /* remember column in incidence list
3248  * Note: each column can be collected at most once for node v, because each column can appear in at most one
3249  * commodity, and in each commodity a column can have at most one +1 and one -1 entry. One of the two +/-1 entries
3250  * is already used for v.
3251  */
3252  assert(ninccols < ncols);
3253  inccols[ninccols] = c;
3254  ninccols++;
3255 
3256  /* update source or target count */
3257  if( s != v )
3258  {
3259  sourcecount[s]++;
3260  u = s;
3261  }
3262  else
3263  {
3264  targetcount[t]++;
3265  u = t;
3266  }
3267 
3268  /* if other node has been seen the first time, store it in adjlist for sparse access of count arrays */
3269  if( sourcecount[u] + targetcount[u] == 1 )
3270  {
3271  assert(nadjnodes < nnodes);
3272  adjnodes[nadjnodes] = u;
3273  nadjnodes++;
3274  }
3275  }
3276  }
3277  }
3278 
3279  /* check if we want to add uncapacitated arcs s -> v or v -> t */
3280  for( l = 0; l < nadjnodes; l++ )
3281  {
3282  int u;
3283 
3284  u = adjnodes[l];
3285  assert(0 <= u && u < nnodes);
3286  assert(sourcecount[u] > 0 || targetcount[u] > 0);
3287  assert(modeltype != SCIP_MCFMODELTYPE_UNDIRECTED || targetcount[u] == 0);
3288  assert(ninccols >= 0);
3289 
3290  /* add arcs u -> v */
3291  if( sourcecount[u] >= arcsthreshold )
3292  {
3293  int arcid;
3294  int m;
3295 
3296  /* create new arc */
3297  SCIP_CALL( createNewArc(scip, mcfdata, u, v, &arcid) );
3298  SCIPdebugMsg(scip, " -> new arc: <%i> = (%i,%i)\n", arcid, u, v);
3299 
3300  /* assign arcid to all involved columns */
3301  for( m = 0; m < ninccols; m++ )
3302  {
3303  int c;
3304  int s;
3305  int t;
3306 
3307  c = inccols[m];
3308  assert(0 <= c && c < ncols);
3309 
3310  assert(cols != NULL);
3311  getIncidentNodes(scip, mcfdata, cols[c], &s, &t);
3312  assert(s == v || t == v);
3313 
3314  if( s == u || (modeltype == SCIP_MCFMODELTYPE_UNDIRECTED && t == u) )
3315  {
3316  SCIPdebugMsg(scip, " -> assign arcid:%i to column <%s>\n", arcid, SCIPvarGetName(SCIPcolGetVar(cols[c])));
3317  colarcid[c] = arcid;
3318 
3319  /* remove column from incidence array */
3320  inccols[m] = inccols[ninccols-1];
3321  ninccols--;
3322  m--;
3323  }
3324  } /*lint --e{850}*/
3325  }
3326 
3327  /* add arcs v -> u */
3328  if( targetcount[u] >= arcsthreshold )
3329  {
3330  int arcid;
3331  int m;
3332 
3333  /* create new arc */
3334  SCIP_CALL( createNewArc(scip, mcfdata, v, u, &arcid) );
3335  SCIPdebugMsg(scip, " -> new arc: <%i> = (%i,%i)\n", arcid, v, u);
3336 
3337  /* assign arcid to all involved columns */
3338  for( m = 0; m < ninccols; m++ )
3339  {
3340  int c;
3341  int s;
3342  int t;
3343 
3344  c = inccols[m];
3345  assert(0 <= c && c < ncols);
3346 
3347  assert(cols != NULL);
3348  getIncidentNodes(scip, mcfdata, cols[c], &s, &t);
3349  assert(s == v || t == v);
3350 
3351  if( t == u )
3352  {
3353  SCIPdebugMsg(scip, " -> assign arcid:%i to column <%s>\n", arcid, SCIPvarGetName(SCIPcolGetVar(cols[c])));
3354  colarcid[c] = arcid;
3355 
3356  /* remove column from incidence array */
3357  inccols[m] = inccols[ninccols-1];
3358  ninccols--;
3359  m--;
3360  }
3361  } /*lint --e{850}*/
3362  }
3363  }
3364 
3365  /* reset sourcecount and targetcount arrays */
3366  for( l = 0; l < nadjnodes; l++ )
3367  {
3368  sourcecount[l] = 0;
3369  targetcount[l] = 0;
3370  }
3371  nadjnodes = 0;
3372 
3373  /* mark the incident columns that could not be assigned to a new arc such that we do not inspect them again */
3374  for( l = 0; l < ninccols; l++ )
3375  {
3376  assert(colarcid[inccols[l]] == -1);
3377  colarcid[inccols[l]] = -2;
3378  }
3379  ninccols = 0;
3380  }
3381  assert(n == nflowcands);
3382  assert(v == nnodes);
3383 
3384 #ifdef SCIP_DEBUG
3385  /* eventually, we must have reset all temporary colarcid[c] = -2 settings to -1 */
3386  for( n = 0; n < ncols; n++ )
3387  assert(colarcid[n] >= -1);
3388 #endif
3389 
3390  /* free temporary memory */
3391  SCIPfreeBufferArray(scip, &inccols);
3392  SCIPfreeBufferArray(scip, &adjnodes);
3393  SCIPfreeBufferArray(scip, &targetcount);
3394  SCIPfreeBufferArray(scip, &sourcecount);
3395  SCIPfreeBufferArray(scip, &sortedflowcandnodeid);
3396  SCIPfreeBufferArray(scip, &sortedflowcands);
3397 
3398  MCFdebugMessage("network after finding uncapacitated arcs has %d nodes, %d arcs, and %d commodities\n",
3399  mcfdata->nnodes, mcfdata->narcs, mcfdata->ncommodities);
3400 
3401  return SCIP_OKAY;
3402 }
3403 
3404 /** cleans up the network: gets rid of commodities without arcs or with at most one node */
3405 static
3407  SCIP* scip, /**< SCIP data structure */
3408  MCFDATA* mcfdata /**< internal MCF extraction data to pass to subroutines */
3409  )
3410 {
3411  int* flowcands = mcfdata->flowcands;
3412  int nflowcands = mcfdata->nflowcands;
3413  int* colcommodity = mcfdata->colcommodity;
3414  int* rowcommodity = mcfdata->rowcommodity;
3415  int* colarcid = mcfdata->colarcid;
3416  int* rowarcid = mcfdata->rowarcid;
3417  int* rownodeid = mcfdata->rownodeid;
3418  int ncommodities = mcfdata->ncommodities;
3419  int* commoditysigns = mcfdata->commoditysigns;
3420  int narcs = mcfdata->narcs;
3421  int nnodes = mcfdata->nnodes;
3422  SCIP_ROW** capacityrows = mcfdata->capacityrows;
3423 
3424  SCIP_ROW** rows;
3425  int nrows;
3426  int ncols;
3427 
3428  int* nnodespercom;
3429  int* narcspercom;
3430  SCIP_Bool* arcisincom;
3431  int* perm;
3432  int permsize;
3433  int maxnnodes;
3434  int nnodesthreshold;
3435  int newncommodities;
3436 
3437  int i;
3438  int a;
3439  int k;
3440 
3441  MCFdebugMessage("network before cleanup has %d nodes, %d arcs, and %d commodities\n", nnodes, narcs, ncommodities);
3442 
3443  /* get LP data */
3444  SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
3445  ncols = SCIPgetNLPCols(scip);
3446 
3447  /* allocate temporary memory */
3448  permsize = ncommodities;
3449  permsize = MAX(permsize, narcs);
3450  permsize = MAX(permsize, nnodes);
3451  SCIP_CALL( SCIPallocBufferArray(scip, &nnodespercom, ncommodities) );
3452  SCIP_CALL( SCIPallocBufferArray(scip, &narcspercom, ncommodities) );
3453  SCIP_CALL( SCIPallocBufferArray(scip, &arcisincom, ncommodities) );
3454  SCIP_CALL( SCIPallocBufferArray(scip, &perm, permsize) );
3455  BMSclearMemoryArray(nnodespercom, ncommodities);
3456  BMSclearMemoryArray(narcspercom, ncommodities);
3457 
3458  /** @todo remove nodes without any incoming and outgoing arcs */
3459 
3460  assert(flowcands != NULL || nflowcands == 0);
3461 
3462  /* count the number of nodes in each commodity */
3463  for( i = 0; i < nflowcands; i++ )
3464  {
3465  int r;
3466 
3467  assert(flowcands != NULL);
3468  r = flowcands[i];
3469  assert(0 <= r && r < nrows);
3470  assert((rownodeid[r] >= 0) == (rowcommodity[r] >= 0));
3471  if( rowcommodity[r] >= 0 )
3472  {
3473  assert(rowcommodity[r] < ncommodities);
3474  nnodespercom[rowcommodity[r]]++;
3475  }
3476  }
3477 
3478  assert(capacityrows != NULL || narcs == 0);
3479 
3480  /* count the number of arcs in each commodity */
3481  for( a = 0; a < narcs; a++ )
3482  {
3483  SCIP_COL** rowcols;
3484  int rowlen;
3485  int r;
3486  int j;
3487 
3488  assert(capacityrows != NULL);
3489  r = SCIProwGetLPPos(capacityrows[a]);
3490  assert(0 <= r && r < nrows);
3491  assert(rowarcid[r] == a);
3492 
3493  /* identify commodities which are touched by this arc capacity constraint */
3494  BMSclearMemoryArray(arcisincom, ncommodities);
3495  rowcols = SCIProwGetCols(rows[r]);
3496  rowlen = SCIProwGetNLPNonz(rows[r]);
3497  for( j = 0; j < rowlen; j++ )
3498  {
3499  int c;
3500 
3501  c = SCIPcolGetLPPos(rowcols[j]);
3502  assert(0 <= c && c < ncols);
3503  if( colcommodity[c] >= 0 && colarcid[c] == a )
3504  {
3505  assert(colcommodity[c] < ncommodities);
3506  arcisincom[colcommodity[c]] = TRUE;
3507  }
3508  }
3509 
3510  /* increase arc counters of touched commodities */
3511  for( k = 0; k < ncommodities; k++ )
3512  {
3513  if( arcisincom[k] )
3514  narcspercom[k]++;
3515  }
3516  }
3517 
3518  /* calculate maximal number of nodes per commodity */
3519  maxnnodes = 0;
3520  for( k = 0; k < ncommodities; k++ )
3521  maxnnodes = MAX(maxnnodes, nnodespercom[k]);
3522 
3523  /* we want to keep only commodities that have at least a certain size relative
3524  * to the largest commodity
3525  */
3526 
3527  nnodesthreshold = (int)(MINCOMNODESFRACTION * maxnnodes);
3528  nnodesthreshold = MAX(nnodesthreshold, MINNODES);
3529  SCIPdebugMsg(scip, " -> node threshold: %d\n", nnodesthreshold);
3530 
3531  /* discard trivial commodities */
3532  newncommodities = 0;
3533  for( k = 0; k < ncommodities; k++ )
3534  {
3535  SCIPdebugMsg(scip, " -> commodity %d: %d nodes, %d arcs\n", k, nnodespercom[k], narcspercom[k]);
3536 
3537  /* only keep commodities of a certain size that have at least one arc */
3538  if( nnodespercom[k] >= nnodesthreshold && narcspercom[k] >= 1 )
3539  {
3540  assert(newncommodities <= k);
3541  perm[k] = newncommodities;
3542  commoditysigns[newncommodities] = commoditysigns[k];
3543  newncommodities++;
3544  }
3545  else
3546  perm[k] = -1;
3547  }
3548 
3549  if( newncommodities < ncommodities )
3550  {
3551  SCIP_Bool* arcisused;
3552  SCIP_Bool* nodeisused;
3553  int newnarcs;
3554  int newnnodes;
3555  int c;
3556  int v;
3557 
3558  SCIPdebugMsg(scip, " -> discarding %d of %d commodities\n", ncommodities - newncommodities, ncommodities);
3559 
3560  SCIP_CALL( SCIPallocBufferArray(scip, &arcisused, narcs) );
3561  SCIP_CALL( SCIPallocBufferArray(scip, &nodeisused, nnodes) );
3562 
3563  /* update data structures to new commodity ids */
3564  BMSclearMemoryArray(arcisused, narcs);
3565  BMSclearMemoryArray(nodeisused, nnodes);
3566  for( c = 0; c < ncols; c++ )
3567  {
3568  if( colcommodity[c] >= 0 )
3569  {
3570  assert(-1 <= colarcid[c] && colarcid[c] < narcs);
3571  assert(colcommodity[c] < mcfdata->ncommodities);
3572  colcommodity[c] = perm[colcommodity[c]];
3573  assert(colcommodity[c] < newncommodities);
3574  if( colcommodity[c] == -1 )
3575  {
3576  /* we are lazy and do not update plusflow and minusflow */
3577  colarcid[c] = -1;
3578  }
3579  else if( colarcid[c] >= 0 )
3580  arcisused[colarcid[c]] = TRUE;
3581  }
3582  }
3583  for( i = 0; i < nflowcands; i++ )
3584  {
3585  int r;
3586 
3587  assert(flowcands != NULL);
3588  r = flowcands[i];
3589  assert(0 <= r && r < nrows);
3590  assert((rownodeid[r] >= 0) == (rowcommodity[r] >= 0));
3591  if( rowcommodity[r] >= 0 )
3592  {
3593  assert(0 <= rownodeid[r] && rownodeid[r] < nnodes);
3594  assert(rowcommodity[r] < mcfdata->ncommodities);
3595  rowcommodity[r] = perm[rowcommodity[r]];
3596  assert(rowcommodity[r] < newncommodities);
3597  if( rowcommodity[r] == -1 )
3598  {
3599  /* we are lazy and do not update flowrowsigns */
3600  rownodeid[r] = -1;
3601  }
3602  else
3603  nodeisused[rownodeid[r]] = TRUE;
3604  }
3605  }
3606 
3607  mcfdata->ncommodities = newncommodities;
3608  ncommodities = newncommodities;
3609 
3610  /* discard unused arcs */
3611  newnarcs = 0;
3612  for( a = 0; a < narcs; a++ )
3613  {
3614  int r;
3615 
3616  assert(capacityrows != NULL);
3617 
3618  if( arcisused[a] )
3619  {
3620  assert(newnarcs <= a);
3621  perm[a] = newnarcs;
3622  capacityrows[newnarcs] = capacityrows[a];
3623  newnarcs++;
3624  }
3625  else
3626  {
3627  /* we are lazy and do not update capacityrowsigns */
3628  perm[a] = -1;
3629  }
3630  r = SCIProwGetLPPos(capacityrows[a]);
3631  assert(0 <= r && r < nrows);
3632  assert(rowarcid[r] == a);
3633  rowarcid[r] = perm[a];
3634  }
3635 
3636  /* update remaining data structures to new arc ids */
3637  if( newnarcs < narcs )
3638  {
3639  SCIPdebugMsg(scip, " -> discarding %d of %d arcs\n", narcs - newnarcs, narcs);
3640 
3641  for( c = 0; c < ncols; c++ )
3642  {
3643  if( colarcid[c] >= 0 )
3644  {
3645  colarcid[c] = perm[colarcid[c]];
3646  assert(colarcid[c] >= 0); /* otherwise colarcid[c] was set to -1 in the colcommodity update */
3647  }
3648  }
3649  mcfdata->narcs = newnarcs;
3650  narcs = newnarcs;
3651  }
3652 #ifndef NDEBUG
3653  for( a = 0; a < narcs; a++ )
3654  {
3655  int r;
3656  assert(capacityrows != NULL);
3657  r = SCIProwGetLPPos(capacityrows[a]);
3658  assert(0 <= r && r < nrows);
3659  assert(rowarcid[r] == a);
3660  }
3661 #endif
3662 
3663  /* discard unused nodes */
3664  newnnodes = 0;
3665  for( v = 0; v < nnodes; v++ )
3666  {
3667  if( nodeisused[v] )
3668  {
3669  assert(newnnodes <= v);
3670  perm[v] = newnnodes;
3671  newnnodes++;
3672  }
3673  else
3674  perm[v] = -1;
3675  }
3676 
3677  /* update data structures to new node ids */
3678  if( newnnodes < nnodes )
3679  {
3680  SCIPdebugMsg(scip, " -> discarding %d of %d nodes\n", nnodes - newnnodes, nnodes);
3681 
3682  for( i = 0; i < nflowcands; i++ )
3683  {
3684  int r;
3685 
3686  assert(flowcands != NULL);
3687  r = flowcands[i];
3688  assert(0 <= r && r < nrows);
3689  assert((rownodeid[r] >= 0) == (rowcommodity[r] >= 0));
3690  if( rowcommodity[r] >= 0 )
3691  {
3692  assert(rowcommodity[r] < ncommodities);
3693  rownodeid[r] = perm[rownodeid[r]];
3694  assert(rownodeid[r] >= 0); /* otherwise we would have deleted the commodity in the rowcommodity update above */
3695  }
3696  }
3697  mcfdata->nnodes = newnnodes;
3698 #ifdef MCF_DEBUG
3699  nnodes = newnnodes;
3700 #endif
3701  }
3702 
3703  /* free temporary memory */
3704  SCIPfreeBufferArray(scip, &nodeisused);
3705  SCIPfreeBufferArray(scip, &arcisused);
3706  }
3707 
3708  /* empty commodities have been removed here */
3709  mcfdata->nemptycommodities = 0;
3710 
3711  /* free temporary memory */
3712  SCIPfreeBufferArray(scip, &perm);
3713  SCIPfreeBufferArray(scip, &arcisincom);
3714  SCIPfreeBufferArray(scip, &narcspercom);
3715  SCIPfreeBufferArray(scip, &nnodespercom);
3716 
3717  MCFdebugMessage("network after cleanup has %d nodes, %d arcs, and %d commodities\n", nnodes, narcs, ncommodities);
3718 
3719  return SCIP_OKAY;
3720 }
3721 
3722 /** for each arc identifies a source and target node */
3723 static
3725  SCIP* scip, /**< SCIP data structure */
3726  MCFDATA* mcfdata, /**< internal MCF extraction data to pass to subroutines */
3727  SCIP_SEPADATA* sepadata, /**< separator data */
3728  MCFEFFORTLEVEL* effortlevel /**< pointer to store effort level of separation */
3729  )
3730 {
3731  int* colarcid = mcfdata->colarcid;
3732  int* colcommodity = mcfdata->colcommodity;
3733  int narcs = mcfdata->narcs;
3734  int nnodes = mcfdata->nnodes;
3735  int ncommodities = mcfdata->ncommodities;
3736  SCIP_ROW** capacityrows = mcfdata->capacityrows;
3737  SCIP_MCFMODELTYPE modeltype = mcfdata->modeltype;
3738  SCIP_Real maxinconsistencyratio = sepadata->maxinconsistencyratio;
3739  SCIP_Real maxarcinconsistencyratio = sepadata->maxarcinconsistencyratio;
3740  int* arcsources;
3741  int* arctargets;
3742  int* colsources;
3743  int* coltargets;
3744  int* firstoutarcs;
3745  int* firstinarcs;
3746  int* nextoutarcs;
3747  int* nextinarcs;
3748 
3749  SCIP_Real *sourcenodecnt;
3750  SCIP_Real *targetnodecnt;
3751  int *flowvarspercom;
3752  int *comtouched;
3753  int *touchednodes;
3754  int ntouchednodes;
3755 
3756  int ncols;
3757  SCIP_Real maxninconsistencies;
3758 
3759  int c;
3760  int v;
3761  int a;
3762 
3763  /* initialize effort level of separation */
3764  assert(effortlevel != NULL);
3765  *effortlevel = MCFEFFORTLEVEL_DEFAULT;
3766 
3767  ncols = SCIPgetNLPCols(scip);
3768 
3769  /* allocate memory in mcfdata */
3770  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->arcsources, narcs) );
3771  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->arctargets, narcs) );
3772  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->colsources, ncols) );
3773  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->coltargets, ncols) );
3774  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->firstoutarcs, nnodes) );
3775  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->firstinarcs, nnodes) );
3776  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->nextoutarcs, narcs) );
3777  SCIP_CALL( SCIPallocMemoryArray(scip, &mcfdata->nextinarcs, narcs) );
3778  arcsources = mcfdata->arcsources;
3779  arctargets = mcfdata->arctargets;
3780  colsources = mcfdata->colsources;
3781  coltargets = mcfdata->coltargets;
3782  firstoutarcs = mcfdata->firstoutarcs;
3783  firstinarcs = mcfdata->firstinarcs;
3784  nextoutarcs = mcfdata->nextoutarcs;
3785  nextinarcs = mcfdata->nextinarcs;
3786 
3787  mcfdata->arcarraysize = narcs;
3788 
3789  /* initialize colsources and coltargets */
3790  for( c = 0; c < ncols; c++ )
3791  {
3792  colsources[c] = -2;
3793  coltargets[c] = -2;
3794  }
3795 
3796  /* initialize adjacency lists */
3797  for( v = 0; v < nnodes; v++ )
3798  {
3799  firstoutarcs[v] = -1;
3800  firstinarcs[v] = -1;
3801  }
3802  for( a = 0; a < narcs; a++ )
3803  {
3804  nextoutarcs[a] = -1;
3805  nextinarcs[a] = -1;
3806  }
3807 
3808  /* allocate temporary memory for source and target node identification */
3809  SCIP_CALL( SCIPallocBufferArray(scip, &sourcenodecnt, nnodes) );
3810  SCIP_CALL( SCIPallocBufferArray(scip, &targetnodecnt, nnodes) );
3811  SCIP_CALL( SCIPallocBufferArray(scip, &flowvarspercom, ncommodities) );
3812  SCIP_CALL( SCIPallocBufferArray(scip, &comtouched, ncommodities) );
3813  SCIP_CALL( SCIPallocBufferArray(scip, &touchednodes, nnodes) );
3814 
3815  BMSclearMemoryArray(sourcenodecnt, nnodes);
3816  BMSclearMemoryArray(targetnodecnt, nnodes);
3817 
3818  mcfdata->ninconsistencies = 0.0;
3819  maxninconsistencies = maxinconsistencyratio * (SCIP_Real)narcs;
3820 
3821  /* search for source and target nodes */
3822  for( a = 0; a < narcs; a++ )
3823  {
3824  SCIP_COL** rowcols;
3825  int rowlen;
3826  int bestsourcev;
3827  int besttargetv;
3828  SCIP_Real bestsourcecnt;
3829  SCIP_Real besttargetcnt;
3830  SCIP_Real totalsourcecnt;
3831  SCIP_Real totaltargetcnt;
3832  SCIP_Real totalnodecnt;
3833  SCIP_Real nsourceinconsistencies;
3834  SCIP_Real ntargetinconsistencies;
3835  int ntouchedcoms;
3836  int i;
3837 #ifndef NDEBUG
3838  int r;
3839 
3840  r = SCIProwGetLPPos(capacityrows[a]);
3841 #endif
3842  assert(0 <= r && r < SCIPgetNLPRows(scip));
3843  assert((mcfdata->capacityrowsigns[r] & (LHSASSIGNED | RHSASSIGNED)) != 0);
3844  assert(mcfdata->rowarcid[r] == a);
3845 
3846 #ifndef NDEBUG
3847  for( i = 0; i < nnodes; i++ )
3848  {
3849  assert(sourcenodecnt[i] == 0);
3850  assert(targetnodecnt[i] == 0);
3851  }
3852 #endif
3853 
3854  rowcols = SCIProwGetCols(capacityrows[a]);
3855  rowlen = SCIProwGetNLPNonz(capacityrows[a]);
3856 
3857  /* count number of flow variables per commodity */
3858  BMSclearMemoryArray(flowvarspercom, ncommodities);
3859  BMSclearMemoryArray(comtouched, ncommodities);
3860  ntouchedcoms = 0;
3861  for( i = 0; i < rowlen; i++ )
3862  {
3863  c = SCIPcolGetLPPos(rowcols[i]);
3864  assert(0 <= c && c < SCIPgetNLPCols(scip));
3865  if( colarcid[c] >= 0 )
3866  {
3867  int k = colcommodity[c];
3868  assert (0 <= k && k < ncommodities);
3869  flowvarspercom[k]++;
3870  if( !comtouched[k] )
3871  {
3872  ntouchedcoms++;
3873  comtouched[k] = TRUE;
3874  }
3875  }
3876  }
3877 
3878  /* if the row does not have any flow variable, it is not a capacity constraint */
3879  if( ntouchedcoms == 0 )
3880  {
3881  capacityrows[a] = NULL;
3882  arcsources[a] = -1;
3883  arctargets[a] = -1;
3884  continue;
3885  }
3886 
3887  /* check the flow variables of the capacity row for flow conservation constraints */
3888  ntouchednodes = 0;
3889  totalsourcecnt = 0.0;
3890  totaltargetcnt = 0.0;
3891  totalnodecnt = 0.0;
3892  for( i = 0; i < rowlen; i++ )
3893  {
3894  c = SCIPcolGetLPPos(rowcols[i]);
3895  assert(0 <= c && c < SCIPgetNLPCols(scip));
3896  if( colarcid[c] >= 0 )
3897  {
3898  int k = colcommodity[c];
3899  int sourcev;
3900  int targetv;
3901  SCIP_Real weight;
3902 
3903  assert (0 <= k && k < ncommodities);
3904  assert (comtouched[k]);
3905  assert (flowvarspercom[k] >= 1);
3906 
3907  /* identify the (at most) two nodes which contain this flow variable */
3908  getIncidentNodes(scip, mcfdata, rowcols[i], &sourcev, &targetv);
3909 
3910  /* count the nodes */
3911  weight = 1.0/flowvarspercom[k];
3912  if( sourcev >= 0 )
3913  {
3914  if( sourcenodecnt[sourcev] == 0.0 && targetnodecnt[sourcev] == 0.0 )
3915  {
3916  touchednodes[ntouchednodes] = sourcev;
3917  ntouchednodes++;
3918  }
3919  sourcenodecnt[sourcev] += weight;
3920  totalsourcecnt += weight;
3921  }
3922  if( targetv >= 0 )
3923  {
3924  if( sourcenodecnt[targetv] == 0.0 && targetnodecnt[targetv] == 0.0 )
3925  {
3926  touchednodes[ntouchednodes] = targetv;
3927  ntouchednodes++;
3928  }
3929  targetnodecnt[targetv] += weight;
3930  totaltargetcnt += weight;
3931  }
3932  if( sourcev >= 0 || targetv >= 0 )
3933  totalnodecnt += weight;
3934  }
3935  }
3936 
3937  /* perform a majority vote on source and target node */
3938  bestsourcev = -1;
3939  besttargetv = -1;
3940  bestsourcecnt = 0.0;
3941  besttargetcnt = 0.0;
3942  for( i = 0; i < ntouchednodes; i++ )
3943  {
3944  v = touchednodes[i];
3945  assert(0 <= v && v < nnodes);
3946 
3947  if( modeltype == SCIP_MCFMODELTYPE_DIRECTED )
3948  {
3949  /* in the directed model, we distinguish between source and target */
3950  if( sourcenodecnt[v] >= targetnodecnt[v] )
3951  {
3952  if( sourcenodecnt[v] > bestsourcecnt )
3953  {
3954  bestsourcev = v;
3955  bestsourcecnt = sourcenodecnt[v];
3956  }
3957  }
3958  else
3959  {
3960  if( targetnodecnt[v] > besttargetcnt )
3961  {
3962  besttargetv = v;
3963  besttargetcnt = targetnodecnt[v];
3964  }
3965  }
3966  }
3967  else
3968  {
3969  SCIP_Real nodecnt = sourcenodecnt[v] + targetnodecnt[v];
3970 
3971  /* in the undirected model, we use source for the maximum and target for the second largest number of total hits */
3972  assert( modeltype == SCIP_MCFMODELTYPE_UNDIRECTED );
3973  if( nodecnt > bestsourcecnt )
3974  {
3975  besttargetv = bestsourcev;
3976  besttargetcnt = bestsourcecnt;
3977  bestsourcev = v;
3978  bestsourcecnt = nodecnt;
3979  }
3980  else if( nodecnt > besttargetcnt )
3981  {
3982  besttargetv = v;
3983  besttargetcnt = nodecnt;
3984  }
3985  }
3986 
3987  /* clear the nodecnt arrays */
3988  sourcenodecnt[v] = 0;
3989  targetnodecnt[v] = 0;
3990  }
3991 
3992  /* check inconsistency of arcs */
3993  if( modeltype == SCIP_MCFMODELTYPE_UNDIRECTED )
3994  {
3995  totalsourcecnt = totalnodecnt;
3996  totaltargetcnt = totalnodecnt;
3997  }
3998  assert(SCIPisGE(scip,totalsourcecnt,bestsourcecnt));
3999  assert(SCIPisGE(scip,totaltargetcnt,besttargetcnt));
4000  nsourceinconsistencies = (totalsourcecnt - bestsourcecnt)/ntouchedcoms;
4001  ntargetinconsistencies = (totaltargetcnt - besttargetcnt)/ntouchedcoms;
4002 
4003  /* delete arcs that have to large inconsistency */
4004  if( nsourceinconsistencies > maxarcinconsistencyratio )
4005  {
4006  /* delete source assignment */
4007  bestsourcev = -1;
4008  }
4009 
4010  if( ntargetinconsistencies > maxarcinconsistencyratio )
4011  {
4012  /* delete target assignment */
4013  besttargetv = -1;
4014  }
4015 
4016  /* assign the incident nodes */
4017  assert(bestsourcev == -1 || bestsourcev != besttargetv);
4018  arcsources[a] = bestsourcev;
4019  arctargets[a] = besttargetv;
4020  SCIPdebugMsg(scip, "arc %d: %d -> %d (len=%d, sourcecnt=%g/%g, targetcnt=%g/%g, %g/%g inconsistencies)\n",
4021  a, bestsourcev, besttargetv, rowlen,
4022  bestsourcecnt, totalsourcecnt, besttargetcnt, totaltargetcnt,
4023  nsourceinconsistencies, ntargetinconsistencies);
4024 
4025  /* update adjacency lists */
4026  if( bestsourcev != -1 )
4027  {
4028  nextoutarcs[a] = firstoutarcs[bestsourcev];
4029  firstoutarcs[bestsourcev] = a;
4030  }
4031  if( besttargetv != -1 )
4032  {
4033  nextinarcs[a] = firstinarcs[besttargetv];
4034  firstinarcs[besttargetv] = a;
4035  }
4036 
4037  /* update the number of inconsistencies */
4038  mcfdata->ninconsistencies += 0.5*(nsourceinconsistencies + ntargetinconsistencies);
4039 
4040  if( mcfdata->ninconsistencies > maxninconsistencies )
4041  {
4042  MCFdebugMessage(" -> reached maximal number of inconsistencies: %g > %g\n",
4043  mcfdata->ninconsistencies, maxninconsistencies);
4044  break;
4045  }
4046  }
4047 
4048  /**@todo should we also use an aggressive parameter setting -- this should be done here */
4049  if( mcfdata->ninconsistencies <= maxninconsistencies && narcs > 0 && ncommodities > 0 )
4050  *effortlevel = MCFEFFORTLEVEL_DEFAULT;
4051  else
4052  *effortlevel = MCFEFFORTLEVEL_OFF;
4053 
4054  MCFdebugMessage("extracted network has %g inconsistencies (ratio %g) -> separating with effort %d\n",
4055  mcfdata->ninconsistencies, mcfdata->ninconsistencies/(SCIP_Real)narcs, *effortlevel);
4056 
4057  /* free temporary memory */
4058  SCIPfreeBufferArray(scip, &touchednodes);
4059  SCIPfreeBufferArray(scip, &comtouched);
4060  SCIPfreeBufferArray(scip, &flowvarspercom);
4061  SCIPfreeBufferArray(scip, &targetnodecnt);
4062  SCIPfreeBufferArray(scip, &sourcenodecnt);
4063 
4064  return SCIP_OKAY;
4065 }
4066 
4067 #define UNKNOWN 0 /**< node has not yet been seen */
4068 #define ONSTACK 1 /**< node is currently on the processing stack */
4069 #define VISITED 2 /**< node has been visited and assigned to some component */
4070 
4071 /** returns lists of nodes and arcs in the connected component of the given startv */
4072 static
4074  SCIP* scip, /**< SCIP data structure */
4075  MCFDATA* mcfdata, /**< internal MCF extraction data to pass to subroutines */
4076  int* nodevisited, /**< array to mark visited nodes */
4077  int startv, /**< node for which the connected component should be generated */
4078  int* compnodes, /**< array to store node ids of the component */
4079  int* ncompnodes, /**< pointer to store the number of nodes in the component */
4080  int* comparcs, /**< array to store arc ids of the component */
4081  int* ncomparcs /**< pointer to store the number of arcs in the component */
4082  )
4083 {
4084  int* arcsources = mcfdata->arcsources;
4085  int* arctargets = mcfdata->arctargets;
4086  int* firstoutarcs = mcfdata->firstoutarcs;
4087  int* firstinarcs = mcfdata->firstinarcs;
4088  int* nextoutarcs = mcfdata->nextoutarcs;
4089  int* nextinarcs = mcfdata->nextinarcs;
4090  int nnodes = mcfdata->nnodes;
4091 
4092  int* stacknodes;
4093  int nstacknodes;
4094 
4095  assert(nodevisited != NULL);
4096  assert(0 <= startv && startv < nnodes);
4097  assert(nodevisited[startv] == UNKNOWN);
4098  assert(compnodes != NULL);
4099  assert(ncompnodes != NULL);
4100  assert(comparcs != NULL);
4101  assert(ncomparcs != NULL);
4102 
4103  *ncompnodes = 0;
4104  *ncomparcs = 0;
4105 
4106  /* allocate temporary memory for node stack */
4107  SCIP_CALL( SCIPallocBufferArray(scip, &stacknodes, nnodes) );
4108 
4109  /* put startv on stack */
4110  stacknodes[0] = startv;
4111  nstacknodes = 1;
4112  nodevisited[startv] = ONSTACK;
4113 
4114  /* perform depth-first search */
4115  while( nstacknodes > 0 )
4116  {
4117  int v;
4118  int a;
4119 
4120  assert(firstoutarcs != NULL);
4121  assert(firstinarcs != NULL);
4122  assert(nextoutarcs != NULL);
4123  assert(nextinarcs != NULL);
4124 
4125  /* pop first element from stack */
4126  v = stacknodes[nstacknodes-1];
4127  nstacknodes--;
4128  assert(0 <= v && v < nnodes);
4129  assert(nodevisited[v] == ONSTACK);
4130  nodevisited[v] = VISITED;
4131 
4132  /* put node into component */
4133  assert(*ncompnodes < nnodes);
4134  compnodes[*ncompnodes] = v;
4135  (*ncompnodes)++;
4136 
4137  /* go through the list of outgoing arcs */
4138  for( a = firstoutarcs[v]; a != -1; a = nextoutarcs[a] )
4139  {
4140  int targetv;
4141 
4142  assert(0 <= a && a < mcfdata->narcs);
4143  assert(arctargets != NULL);
4144 
4145  targetv = arctargets[a];
4146 
4147  /* check if we have already visited the target node */
4148  if( targetv != -1 && nodevisited[targetv] == VISITED )
4149  continue;
4150 
4151  /* put arc to component */
4152  assert(*ncomparcs < mcfdata->narcs);
4153  comparcs[*ncomparcs] = a;
4154  (*ncomparcs)++;
4155 
4156  /* push target node to stack */
4157  if( targetv != -1 && nodevisited[targetv] == UNKNOWN )
4158  {
4159  assert(nstacknodes < nnodes);
4160  stacknodes[nstacknodes] = targetv;
4161  nstacknodes++;
4162  nodevisited[targetv] = ONSTACK;
4163  }
4164  }
4165 
4166  /* go through the list of ingoing arcs */
4167  for( a = firstinarcs[v]; a != -1; a = nextinarcs[a] )
4168  {
4169  int sourcev;
4170 
4171  assert(0 <= a && a < mcfdata->narcs);
4172  assert(arcsources != NULL);
4173 
4174  sourcev = arcsources[a];
4175 
4176  /* check if we have already seen the source node */
4177  if( sourcev != -1 && nodevisited[sourcev] == VISITED )
4178  continue;
4179 
4180  /* put arc to component */
4181  assert(*ncomparcs < mcfdata->narcs);
4182  comparcs[*ncomparcs] = a;
4183  (*ncomparcs)++;
4184 
4185  /* push source node to stack */
4186  if( sourcev != -1 && nodevisited[sourcev] == UNKNOWN )
4187  {
4188  assert(nstacknodes < nnodes);
4189  stacknodes[nstacknodes] = sourcev;
4190  nstacknodes++;
4191  nodevisited[sourcev] = ONSTACK;
4192  }
4193  }
4194  }
4195 
4196  /* free temporary memory */
4197  SCIPfreeBufferArray(scip, &stacknodes);
4198 
4199  return SCIP_OKAY;
4200 }
4201 
4202 /** extracts MCF network structures from the current LP */
4203 static
4205  SCIP* scip, /**< SCIP data structure */
4206  SCIP_SEPADATA* sepadata, /**< separator data */
4207  SCIP_MCFNETWORK*** mcfnetworks, /**< pointer to store array of MCF network structures */
4208  int* nmcfnetworks, /**< pointer to store number of MCF networks */
4209  MCFEFFORTLEVEL* effortlevel /**< effort level of separation */
4210  )
4211 {
4212  MCFDATA mcfdata;
4213 
4214  SCIP_MCFMODELTYPE modeltype = (SCIP_MCFMODELTYPE) sepadata->modeltype;
4215 
4216  SCIP_Bool failed;
4217 
4218  SCIP_ROW** rows;
4219  SCIP_COL** cols;
4220  int nrows;
4221  int ncols;
4222  int mcfnetworkssize;
4223 
4224  assert(mcfnetworks != NULL);
4225  assert(nmcfnetworks != NULL);
4226  assert(effortlevel != NULL);
4227 
4228  failed = FALSE;
4229  *effortlevel = MCFEFFORTLEVEL_OFF;
4230  *mcfnetworks = NULL;
4231  *nmcfnetworks = 0;
4232  mcfnetworkssize = 0;
4233 
4234  /* Algorithm to identify multi-commodity-flow network with capacity constraints
4235  *
4236  * 1. Identify candidate rows for flow conservation constraints in the LP.
4237  * 2. Sort flow conservation candidates by a ranking on how sure we are that it is indeed a constraint of the desired type.
4238  * 3. Extract network structure of flow conservation constraints:
4239  * (a) Initialize plusflow[c] = minusflow[c] = FALSE for all columns c and other local data.
4240  * (b) As long as there are flow conservation candidates left:
4241  * (i) Create new commodity and use first flow conservation constraint as new row.
4242  * (ii) Add new row to commodity, update pluscom/minuscom accordingly.
4243  * (iii) For the newly added columns search for an incident flow conservation constraint. Pick the one of highest ranking.
4244  * Reflect row or commodity if necessary (multiply with -1)
4245  * (iv) If found, set new row to this row and goto (ii).
4246  * (v) If only very few flow rows have been used, discard the commodity immediately.
4247  * 4. Identify candidate rows for capacity constraints in the LP.
4248  * 5. Sort capacity constraint candidates by a ranking on how sure we are that it is indeed a constraint of the desired type.
4249  * 6. Identify capacity constraints for the arcs and assign arc ids to columns and capacity constraints.
4250  * 7. Assign node ids to flow conservation constraints.
4251  * 8. PostProcessing
4252  * a if there are still undecided commodity signs, fix them to +1
4253  * b clean up the network: get rid of commodities without arcs or with at most one node
4254  * c assign source and target nodes to capacitated arc
4255  * d find uncapacitated arcs
4256  */
4257 
4258  /* get LP data */
4259  SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
4260  SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
4261 
4262  /* initialize local extraction data */
4263  mcfdata.flowrowsigns = NULL;
4264  mcfdata.flowrowscalars = NULL;
4265  mcfdata.flowrowscores = NULL;
4266  mcfdata.capacityrowsigns = NULL;
4267  mcfdata.capacityrowscores = NULL;
4268  mcfdata.flowcands = NULL;
4269  mcfdata.nflowcands = 0;
4270  mcfdata.capacitycands = NULL;
4271  mcfdata.ncapacitycands = 0;
4272  mcfdata.plusflow = NULL;
4273  mcfdata.minusflow = NULL;
4274  mcfdata.ncommodities = 0;
4275  mcfdata.nemptycommodities = 0;
4276  mcfdata.commoditysigns = NULL;
4277  mcfdata.commoditysignssize = 0;
4278  mcfdata.colcommodity = NULL;
4279  mcfdata.rowcommodity = NULL;
4280  mcfdata.colarcid = NULL;
4281  mcfdata.rowarcid = NULL;
4282  mcfdata.rownodeid = NULL;
4283  mcfdata.arcarraysize = 0;
4284  mcfdata.arcsources = NULL;
4285  mcfdata.arctargets = NULL;
4286  mcfdata.colsources = NULL;
4287  mcfdata.coltargets = NULL;
4288  mcfdata.firstoutarcs = NULL;
4289  mcfdata.firstinarcs = NULL;
4290  mcfdata.nextoutarcs = NULL;
4291  mcfdata.nextinarcs = NULL;
4292  mcfdata.newcols = NULL;
4293  mcfdata.nnewcols = 0;
4294  mcfdata.narcs = 0;
4295  mcfdata.nnodes = 0;
4296  mcfdata.ninconsistencies = 0.0;
4297  mcfdata.capacityrows = NULL;
4298  mcfdata.capacityrowssize = 0;
4299  mcfdata.colisincident = NULL;
4300  mcfdata.zeroarcarray = NULL;
4301  mcfdata.modeltype = modeltype;
4302 
4303  /* 1. identify candidate rows for flow conservation constraints in the LP
4304  * 2. Sort flow conservation candidates by a ranking on how sure we are that it is indeed a constraint of the desired type
4305  */
4306  SCIP_CALL( extractFlowRows(scip, &mcfdata) );
4307  assert(mcfdata.flowrowsigns != NULL);
4308  assert(mcfdata.flowrowscalars != NULL);
4309  assert(mcfdata.flowrowscores != NULL);
4310  assert(mcfdata.flowcands != NULL);
4311 
4312  if( mcfdata.nflowcands == 0 )
4313  failed = TRUE;
4314 
4315  if( !failed )
4316  {
4317  /* 3. extract network structure of flow conservation constraints. */
4318  SCIP_CALL( extractFlow(scip, &mcfdata, MAXFLOWVARFLOWROWRATIO, &failed) );
4319  assert(mcfdata.plusflow != NULL);
4320  assert(mcfdata.minusflow != NULL);
4321  assert(mcfdata.colcommodity != NULL);
4322  assert(mcfdata.rowcommodity != NULL);
4323  assert(mcfdata.newcols != NULL);
4324  }
4325 
4326  if( !failed )
4327  {
4328 #ifdef SCIP_DEBUG
4329  printCommodities(scip, &mcfdata);
4330 #endif
4331 
4332  /* 4. identify candidate rows for capacity constraints in the LP
4333  * 5. sort capacity constraint candidates by a ranking on how sure we are that it is indeed a constraint of the desired type
4334  */
4335  SCIP_CALL( extractCapacityRows(scip, &mcfdata) );
4336  assert(mcfdata.capacityrowsigns != NULL);
4337  assert(mcfdata.capacityrowscores != NULL);
4338  assert(mcfdata.capacitycands != NULL);
4339 
4340  if( mcfdata.ncapacitycands == 0 )
4341  failed = TRUE;
4342  }
4343 
4344  if( !failed )
4345  {
4346  /* 6. arc-detection -- identify capacity constraints for the arcs and assign arc ids to columns and capacity constraints */
4347  SCIP_CALL( extractCapacities(scip, &mcfdata) );
4348  assert(mcfdata.colarcid != NULL);
4349  assert(mcfdata.rowarcid != NULL);
4350 
4351  /* 7. node-detection -- assign node ids to flow conservation constraints */
4352  SCIP_CALL( extractNodes(scip, &mcfdata) );
4353  assert(mcfdata.rownodeid != NULL);
4354  assert(mcfdata.colisincident != NULL);
4355  assert(mcfdata.zeroarcarray != NULL);
4356 
4357  /* 8. postprocessing */
4358  /* 8.a if there are still undecided commodity signs, fix them to +1 */
4359  fixCommoditySigns(scip, &mcfdata);
4360 
4361  /* 8.b clean up the network: get rid of commodities without arcs or with at most one node */
4362  SCIP_CALL( cleanupNetwork(scip, &mcfdata) );
4363 
4364  /* 8.c construct incidence function -- assign source and target nodes to capacitated arcs */
4365  SCIP_CALL( identifySourcesTargets(scip, &mcfdata, sepadata, effortlevel) );
4366  assert(mcfdata.arcsources != NULL);
4367  assert(mcfdata.arctargets != NULL);
4368  assert(mcfdata.colsources != NULL);
4369  assert(mcfdata.coltargets != NULL);
4370  assert(mcfdata.firstoutarcs != NULL);
4371  assert(mcfdata.firstinarcs != NULL);
4372  assert(mcfdata.nextoutarcs != NULL);
4373  assert(mcfdata.nextinarcs != NULL);
4374  }
4375 
4376  if( !failed && *effortlevel != MCFEFFORTLEVEL_OFF)
4377  {
4378  int* nodevisited;
4379  int* compnodeid;
4380  int* compnodes;
4381  int* comparcs;
4382  int minnodes;
4383  int v;
4384 
4385  /* 8.d find uncapacitated arcs */
4386  SCIP_CALL( findUncapacitatedArcs(scip, &mcfdata) );
4387 
4388 #ifdef SCIP_DEBUG
4389  printCommodities(scip, &mcfdata);
4390 #endif
4391 
4392  minnodes = MINNODES;
4393 
4394  /* allocate temporary memory for component finding */
4395  SCIP_CALL( SCIPallocBufferArray(scip, &nodevisited, mcfdata.nnodes) );
4396  SCIP_CALL( SCIPallocBufferArray(scip, &compnodes, mcfdata.nnodes) );
4397  SCIP_CALL( SCIPallocBufferArray(scip, &comparcs, mcfdata.narcs) );
4398  BMSclearMemoryArray(nodevisited, mcfdata.nnodes);
4399 
4400  /* allocate temporary memory for v -> compv mapping */
4401  SCIP_CALL( SCIPallocBufferArray(scip, &compnodeid, mcfdata.nnodes) );
4402  for( v = 0; v < mcfdata.nnodes; v++ )
4403  compnodeid[v] = -1;
4404 
4405  /* search components and create a network structure for each of them */
4406  for( v = 0; v < mcfdata.nnodes; v++ )
4407  {
4408  int ncompnodes;
4409  int ncomparcs;
4410 
4411  /* ignore nodes that have been already assigned to a component */
4412  assert(nodevisited[v] == UNKNOWN || nodevisited[v] == VISITED);
4413  if( nodevisited[v] == VISITED )
4414  continue;
4415 
4416  /* identify nodes and arcs of this component */
4417  SCIP_CALL( identifyComponent(scip, &mcfdata, nodevisited, v, compnodes, &ncompnodes, comparcs, &ncomparcs) );
4418  assert(ncompnodes >= 1);
4419  assert(compnodes[0] == v);
4420  assert(nodevisited[v] == VISITED);
4421 
4422  /* ignore network component if it is trivial */
4423  if( ncompnodes >= minnodes && ncomparcs >= MINARCS )
4424  {
4425  SCIP_MCFNETWORK* mcfnetwork;
4426  int i;
4427 
4428  /* make sure that we have enough memory for the new network pointer */
4429  assert(*nmcfnetworks <= MAXNETWORKS);
4430  assert(*nmcfnetworks <= mcfnetworkssize);
4431  if( *nmcfnetworks == mcfnetworkssize )
4432  {
4433  mcfnetworkssize = MAX(2*mcfnetworkssize, *nmcfnetworks+1);
4434  SCIP_CALL( SCIPreallocMemoryArray(scip, mcfnetworks, mcfnetworkssize) );
4435  }
4436  assert(*nmcfnetworks < mcfnetworkssize);
4437 
4438  /* create network data structure */
4439  SCIP_CALL( mcfnetworkCreate(scip, &mcfnetwork) );
4440 
4441  /* fill sparse network structure */
4442  SCIP_CALL( mcfnetworkFill(scip, mcfnetwork, &mcfdata, compnodeid, compnodes, ncompnodes, comparcs, ncomparcs) );
4443 
4444  /* insert in sorted network list */
4445  for( i = *nmcfnetworks; i > 0 && mcfnetwork->nnodes > (*mcfnetworks)[i-1]->nnodes; i-- )
4446  (*mcfnetworks)[i] = (*mcfnetworks)[i-1];
4447  (*mcfnetworks)[i] = mcfnetwork;
4448  (*nmcfnetworks)++;
4449 
4450  /* if we reached the maximal number of networks, update minnodes */
4451  if( *nmcfnetworks >= MAXNETWORKS )
4452  minnodes = MAX(minnodes, (*mcfnetworks)[*nmcfnetworks-1]->nnodes);
4453 
4454  /* if we exceeded the maximal number of networks, delete the last one */
4455  if( *nmcfnetworks > MAXNETWORKS )
4456  {
4457  SCIPdebugMsg(scip, " -> discarded network with %d nodes and %d arcs due to maxnetworks (minnodes=%d)\n",
4458  (*mcfnetworks)[*nmcfnetworks-1]->nnodes, (*mcfnetworks)[*nmcfnetworks-1]->narcs, minnodes);
4459  SCIP_CALL( mcfnetworkFree(scip, &(*mcfnetworks)[*nmcfnetworks-1]) );
4460  (*nmcfnetworks)--;
4461  }
4462  assert(*nmcfnetworks <= MAXNETWORKS);
4463  }
4464  else
4465  {
4466  SCIPdebugMsg(scip, " -> discarded component with %d nodes and %d arcs\n", ncompnodes, ncomparcs);
4467  }
4468 
4469  }
4470 
4471  /* free temporary memory */
4472  SCIPfreeBufferArray(scip, &compnodeid);
4473  SCIPfreeBufferArray(scip, &comparcs);
4474  SCIPfreeBufferArray(scip, &compnodes);
4475  SCIPfreeBufferArray(scip, &nodevisited);
4476  }
4477 
4478  /* free memory */
4479  SCIPfreeMemoryArrayNull(scip, &mcfdata.arcsources);
4480  SCIPfreeMemoryArrayNull(scip, &mcfdata.arctargets);
4481  SCIPfreeMemoryArrayNull(scip, &mcfdata.colsources);
4482  SCIPfreeMemoryArrayNull(scip, &mcfdata.coltargets);
4483  SCIPfreeMemoryArrayNull(scip, &mcfdata.firstoutarcs);
4484  SCIPfreeMemoryArrayNull(scip, &mcfdata.firstinarcs);
4485  SCIPfreeMemoryArrayNull(scip, &mcfdata.nextoutarcs);
4486  SCIPfreeMemoryArrayNull(scip, &mcfdata.nextinarcs);
4487  SCIPfreeMemoryArrayNull(scip, &mcfdata.zeroarcarray);
4488  SCIPfreeMemoryArrayNull(scip, &mcfdata.colisincident);
4489  SCIPfreeMemoryArrayNull(scip, &mcfdata.capacityrows);
4490  SCIPfreeMemoryArrayNull(scip, &mcfdata.rownodeid);
4491  SCIPfreeMemoryArrayNull(scip, &mcfdata.rowarcid);
4492  SCIPfreeMemoryArrayNull(scip, &mcfdata.colarcid);
4493  SCIPfreeMemoryArrayNull(scip, &mcfdata.newcols);
4494  SCIPfreeMemoryArrayNull(scip, &mcfdata.rowcommodity);
4495  SCIPfreeMemoryArrayNull(scip, &mcfdata.colcommodity);
4496  SCIPfreeMemoryArrayNull(scip, &mcfdata.commoditysigns);
4497  SCIPfreeMemoryArrayNull(scip, &mcfdata.minusflow);
4498  SCIPfreeMemoryArrayNull(scip, &mcfdata.plusflow);
4499  SCIPfreeMemoryArrayNull(scip, &mcfdata.capacitycands);
4500  SCIPfreeMemoryArrayNull(scip, &mcfdata.flowcands);
4501  SCIPfreeMemoryArrayNull(scip, &mcfdata.capacityrowscores);
4502  SCIPfreeMemoryArrayNull(scip, &mcfdata.capacityrowsigns);
4503  SCIPfreeMemoryArrayNull(scip, &mcfdata.flowrowscores);
4504  SCIPfreeMemoryArrayNull(scip, &mcfdata.flowrowscalars);
4505  SCIPfreeMemoryArrayNull(scip, &mcfdata.flowrowsigns);
4506 
4507  return SCIP_OKAY;
4508 }
4509 #ifdef COUNTNETWORKVARIABLETYPES
4510 /** extracts MCF network structures from the current LP */
4511 static
4512 SCIP_RETCODE printFlowSystemInfo(
4513  SCIP* scip, /**< SCIP data structure */
4514  SCIP_MCFNETWORK** mcfnetworks, /**< array of MCF network structures */
4515  int nmcfnetworks /**< number of MCF networks */
4516  )
4517 {
4518  SCIP_ROW** rows;
4519  SCIP_COL** cols;
4520  SCIP_Bool* colvisited;
4521  int nrows;
4522  int ncols;
4523  int m;
4524  int c;
4525  int a;
4526  int k;
4527  int v;
4528  int nflowrows = 0;
4529  int ncaprows = 0;
4530  int nflowvars = 0;
4531  int nintflowvars = 0;
4532  int nbinflowvars = 0;
4533  int ncontflowvars = 0;
4534  int ncapvars = 0;
4535  int nintcapvars = 0;
4536  int nbincapvars = 0;
4537  int ncontcapvars = 0;
4538 
4539  /* get LP data */
4540  SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
4541  SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
4542  SCIP_CALL( SCIPallocBufferArray(scip, &colvisited, ncols) );
4543  /* get flow variable types */
4544  for(c=0; c < ncols; c++)
4545  colvisited[c]=FALSE;
4546 
4547 
4548  MCFdebugMessage("\n\n****** VAR COUNTING ********* \n");
4549 
4550  for(m=0; m < nmcfnetworks; m++)
4551  {
4552  SCIP_MCFNETWORK* mcfnetwork = mcfnetworks[m];
4553 
4554  int narcs = mcfnetwork->narcs;
4555  int nnodes = mcfnetwork->nnodes;
4556  int ncommodities = mcfnetwork->ncommodities;
4557  SCIP_ROW** arccapacityrows = mcfnetwork->arccapacityrows;
4558  SCIP_ROW*** nodeflowrows = mcfnetwork->nodeflowrows;
4559  int* colcommodity = mcfnetwork->colcommodity;
4560 
4561  /* get flow variable types */
4562  for(c=0; c < ncols; c++)
4563  {
4564  SCIP_COL* col;
4565 
4566  if(colcommodity[c] >= 0 && ! colvisited[c])
4567  {
4568  /* this is a flow variable */
4569  nflowvars++;
4570  col = cols[c];
4571  colvisited[c] = TRUE;
4572  switch( SCIPvarGetType(SCIPcolGetVar(col)) )
4573  {
4574  case SCIP_VARTYPE_BINARY:
4575  nbinflowvars++;
4576  break;
4577  case SCIP_VARTYPE_INTEGER:
4578  nintflowvars++;
4579  break;
4580  case SCIP_VARTYPE_IMPLINT:
4581  nintflowvars++;
4582  break;
4584  ncontflowvars++;
4585  break;
4586  default:
4587  SCIPerrorMessage("unknown variable type\n");
4588  SCIPABORT();
4589  return SCIP_INVALIDDATA; /*lint !e527*/
4590  }
4591  }
4592  }
4593  /* get capacity variable types and number of capacity rows*/
4594  for( a = 0; a < narcs; a++ )
4595  {
4596  SCIP_ROW* row;
4597  row = arccapacityrows[a];
4598 
4599  if( row != NULL )
4600  {
4601  SCIP_COL** rowcols;
4602  int rowlen;
4603  int i;
4604 
4605  ncaprows++;
4606  rowcols = SCIProwGetCols(row);
4607  rowlen = SCIProwGetNLPNonz(row);
4608 
4609  for( i = 0; i < rowlen; i++ )
4610  {
4611  c = SCIPcolGetLPPos(rowcols[i]);
4612  assert(0 <= c && c < SCIPgetNLPCols(scip));
4613 
4614  if(colcommodity[c] == -1 && ! colvisited[c] )
4615  {
4616  ncapvars++;
4617  colvisited[c] = TRUE;
4618  switch( SCIPvarGetType(SCIPcolGetVar(rowcols[i]) ) )
4619  {
4620  case SCIP_VARTYPE_BINARY:
4621  nbincapvars++;
4622  break;
4623  case SCIP_VARTYPE_INTEGER:
4624  nintcapvars++;
4625  break;
4626  case SCIP_VARTYPE_IMPLINT:
4627  nintcapvars++;
4628  break;
4630  ncontcapvars++;
4631  break;
4632  default:
4633  SCIPerrorMessage("unknown variable type\n");
4634  SCIPABORT();
4635  return SCIP_INVALIDDATA; /*lint !e527*/
4636  }
4637  }
4638  }
4639  }
4640  }
4641  /* get number of flow rows */
4642  for( k = 0; k < ncommodities; k++ )
4643  {
4644  for( v = 0; v < nnodes; v++ )
4645  {
4646  SCIP_ROW* row;
4647  row = nodeflowrows[v][k];
4648 
4649  if( row != NULL )
4650  nflowrows++;
4651  }
4652  }
4653 
4654  MCFdebugMessage("----- network %i -----\n",m);
4655  MCFdebugMessage(" nof flowrows: %5d\n", nflowrows);
4656  MCFdebugMessage(" nof caprows: %5d\n", ncaprows);
4657  MCFdebugMessage(" nof flowvars: %5d of which [ %d , %d , %d ] are continuous, integer, binary\n",
4658  nflowvars, ncontflowvars, nintflowvars, nbinflowvars);
4659  MCFdebugMessage(" nof capvars: %5d of which [ %d , %d , %d ] are continuous, integer, binary\n",
4660  ncapvars, ncontcapvars, nintcapvars, nbincapvars);
4661  }
4662 
4663 
4664  MCFdebugMessage("****** END VAR COUNTING ********* \n\n");
4665 
4666  SCIPfreeBufferArray(scip, &colvisited);
4667 
4668  return SCIP_OKAY;
4669 }
4670 #endif
4671 /*
4672  * Union find methods
4673  * used for generating partitions of node sets and
4674  * for checking connectivity of cut shores
4675  */
4676 
4677 /** initializes a union find data structure by putting each element into its own set */
4678 static
4680  int* representatives, /**< mapping an element v to its representative */
4681  int nelems /**< number of elements in the ground set */
4682  )
4683 {
4684  int v;
4685 
4686  /* we start with each element being in its own set */
4687  for( v = 0; v < nelems; v++ )
4688  representatives[v] = v;
4689 }
4690 
4691 /** applies a union find algorithm to get the representative of v */
4692 static
4694  int* representatives, /**< mapping an element v to its representative */
4695  int v /**< element v to get a representative for */
4696  )
4697 {
4698  assert(representatives != NULL);
4699 
4700  while( v != representatives[v] )
4701  {
4702  representatives[v] = representatives[representatives[v]];
4703  v = representatives[v];
4704  }
4705 
4706  return v;
4707 }
4708 
4709 /** joins two sets in the union find framework */
4710 static
4712  int* representatives, /**< mapping an element v to its representative */
4713  int rep1, /**< representative of first set */
4714  int rep2 /**< representative of second set */
4715  )
4716 {
4717  assert(rep1 != rep2);
4718  assert(representatives[rep1] == rep1);
4719  assert(representatives[rep2] == rep2);
4720 
4721  /* make sure that the smaller representative survives
4722  * -> element 0 is always a representative
4723  */
4724  if( rep1 < rep2 )
4725  representatives[rep2] = rep1;
4726  else
4727  representatives[rep1] = rep2;
4728 }
4729 
4730 /*
4731  * Node pair methods
4732  * used for shrinking the network based on nodepair-weights
4733  * -> creating partition
4734 */
4735 
4736 /** comparison method for weighted nodepairs */
4737 static
4739 {
4740  NODEPAIRENTRY* nodepair1 = (NODEPAIRENTRY*)elem1;
4741  NODEPAIRENTRY* nodepair2 = (NODEPAIRENTRY*)elem2;
4742 
4743  if( nodepair1->weight > nodepair2->weight )
4744  return -1;
4745  else if( nodepair1->weight < nodepair2->weight )
4746  return +1;
4747  else
4748  return 0;
4749 }
4750 
4751 /** NodePair HashTable
4752  * gets the key of the given element */
4753 static
4754 SCIP_DECL_HASHGETKEY(hashGetKeyNodepairs)
4755 {
4756  /*lint --e{715}*/
4757  /* the key is the element itself */
4758  return elem;
4759 }
4760 
4761 /** NodePair HashTable
4762  * returns TRUE iff both keys are equal;
4763  * two nodepairs are equal if both nodes equal
4764  */
4765 static
4766 SCIP_DECL_HASHKEYEQ(hashKeyEqNodepairs)
4767 {
4768 #ifndef NDEBUG
4769  SCIP_MCFNETWORK* mcfnetwork;
4770 #endif
4771  NODEPAIRENTRY* nodepair1;
4772  NODEPAIRENTRY* nodepair2;
4773  int source1;
4774  int source2;
4775  int target1;
4776  int target2;
4777 
4778 #ifndef NDEBUG
4779  mcfnetwork = (SCIP_MCFNETWORK*)userptr;
4780  assert(mcfnetwork != NULL);
4781 #endif
4782 
4783  nodepair1 = (NODEPAIRENTRY*)key1;
4784  nodepair2 = (NODEPAIRENTRY*)key2;
4785 
4786  assert(nodepair1 != NULL);
4787  assert(nodepair2 != NULL);
4788 
4789  source1 = nodepair1->node1;
4790  source2 = nodepair2->node1;
4791  target1 = nodepair1->node2;
4792  target2 = nodepair2->node2;
4793 
4794  assert(source1 >=0 && source1 < mcfnetwork->nnodes);
4795  assert(source2 >=0 && source2 < mcfnetwork->nnodes);
4796  assert(target1 >=0 && target1 < mcfnetwork->nnodes);
4797  assert(target2 >=0 && target2 < mcfnetwork->nnodes);
4798  assert(source1 <= target1);
4799  assert(source2 <= target2);
4800 
4801  return (source1 == source2 && target1 == target2);
4802 }
4803 
4804 /** NodePair HashTable
4805  * returns the hash value of the key */
4806 static
4807 SCIP_DECL_HASHKEYVAL(hashKeyValNodepairs)
4808 {
4809 #ifndef NDEBUG
4810  SCIP_MCFNETWORK* mcfnetwork;
4811 #endif
4812  NODEPAIRENTRY* nodepair;
4813  int source;
4814  int target;
4815  unsigned int hashval;
4816 
4817 #ifndef NDEBUG
4818  mcfnetwork = (SCIP_MCFNETWORK*)userptr;
4819  assert(mcfnetwork != NULL);
4820 #endif
4821 
4822  nodepair = (NODEPAIRENTRY*)key;
4823  assert( nodepair != NULL);
4824 
4825  source = nodepair->node1;
4826  target = nodepair->node2;
4827 
4828  assert(source >=0 && source < mcfnetwork->nnodes);
4829  assert(target >=0 && target < mcfnetwork->nnodes);
4830  assert(source <= target);
4831 
4832  hashval = (source << 16) + target; /*lint !e701*/
4833 
4834  return hashval;
4835 }
4836 
4837 /** creates a priority queue and fills it with the given nodepair entries
4838  *
4839  */
4840 static
4842  SCIP* scip, /**< SCIP data structure */
4843  SCIP_MCFNETWORK* mcfnetwork, /**< MCF network structure */
4844  NODEPAIRQUEUE** nodepairqueue /**< pointer to nodepair priority queue */
4845  )
4846 {
4847  /* For every nodepair that is used in the network (at least one arc exists having this nodepair as endnodes)
4848  * we calculate a weight:
4849  * The weight w_st of a nodepair (s,t) is the minimum of the weights of all s-t and t-s arcs
4850  * The weight w_a of an arc a is calculated as:
4851  * w_a : = s_a + pi_a
4852  * where s_a>=0 is the slack of the capacity constraint and pi_a<=0 its dual.
4853  * The weight of uncapacitated arcs (without capacity constraints) is infinite.
4854  */
4855 #ifdef BETTERWEIGHTFORDEMANDNODES
4856  int ncommodities;
4859  SCIP_Real maxweight;
4860  SCIP_Real minweight;
4861 #endif
4862 
4863 #ifdef TIEBREAKING
4864  int* colcommodity;
4865 #endif
4866 
4867 
4868  SCIP_HASHTABLE* hashtable;
4869  NODEPAIRENTRY* nodepairs;
4870 
4871  int hashtablesize;
4872  int a;
4873  int nnodepairs;
4874  int n;
4875 
4876  assert(mcfnetwork != NULL);
4877 
4878 #ifdef BETTERWEIGHTFORDEMANDNODES
4879  ncommodities = mcfnetwork->ncommodities;
4880  nodeflowrows = mcfnetwork->nodeflowrows;
4881  nodeflowscales = mcfnetwork->nodeflowscales;
4882 #endif
4883 
4884 #ifdef TIEBREAKING
4885  colcommodity = mcfnetwork->colcommodity;
4886 #endif
4887 
4888  assert(nodepairqueue != NULL);
4889 
4890  SCIP_CALL( SCIPallocBuffer(scip, nodepairqueue) );
4891 
4892  /* create a hash table for all used node pairs
4893  * hash table is only needed to have unique nodepairs (identify arcs using the same nodepair)
4894  */
4895  hashtablesize = mcfnetwork->narcs;
4896  hashtablesize = MAX(hashtablesize, HASHSIZE_NODEPAIRS);
4897  SCIP_CALL( SCIPhashtableCreate(&hashtable, SCIPblkmem(scip), hashtablesize,
4898  hashGetKeyNodepairs, hashKeyEqNodepairs, hashKeyValNodepairs, (void*) mcfnetwork) );
4899 
4900  /* nodepairs will contain all constructed nodepairs and is used to fill the priority queue */
4901  SCIP_CALL( SCIPallocBufferArray(scip, &(*nodepairqueue)->nodepairs, mcfnetwork->narcs) );
4902 
4903  /* initialize hash table of all used node pairs and fill nodepairs */
4904  nnodepairs = 0;
4905  for( a = 0; a < mcfnetwork->narcs; a++ )
4906  {
4907  NODEPAIRENTRY nodepair;
4908  NODEPAIRENTRY* nodepairptr;
4909  SCIP_ROW* capacityrow;
4910 
4911  capacityrow = mcfnetwork->arccapacityrows[a];
4912 
4913  SCIPdebugMsg(scip, "arc %i = (%i %i)\n", a, mcfnetwork->arcsources[a], mcfnetwork->arctargets[a]);
4914 
4915  /* construct fresh nodepair: smaller node gets node1 in nodeentry */
4916  if( mcfnetwork->arcsources[a] <= mcfnetwork->arctargets[a] )
4917  {
4918  nodepair.node1 = mcfnetwork->arcsources[a];
4919  nodepair.node2 = mcfnetwork->arctargets[a];
4920  }
4921  else
4922  {
4923  nodepair.node2 = mcfnetwork->arcsources[a];
4924  nodepair.node1 = mcfnetwork->arctargets[a];
4925  }
4926 
4927  assert(nodepair.node1 < mcfnetwork->nnodes);
4928  assert(nodepair.node2 < mcfnetwork->nnodes);
4929  if( nodepair.node1 == -1 || nodepair.node2 == -1 )
4930  continue;
4931 
4932  /* construct arc weight of a */
4933  if( capacityrow != NULL )
4934  {
4935  SCIP_Real maxval;
4936  SCIP_Real slack;
4937  SCIP_Real dualsol;
4938  SCIP_Real scale;
4939 #ifdef TIEBREAKING
4940  SCIP_Real totalflow;
4941  SCIP_Real totalcap;
4942  SCIP_COL** rowcols;
4943  int rowlen;
4944  int i;
4945  int c;
4946 #endif
4947 
4948  slack = SCIPgetRowFeasibility(scip, mcfnetwork->arccapacityrows[a]);
4949  slack = MAX(slack, 0.0); /* can only be negative due to numerics */
4950  dualsol = SCIProwGetDualsol(mcfnetwork->arccapacityrows[a]);
4951  maxval = SCIPgetRowMaxCoef(scip, mcfnetwork->arccapacityrows[a]);
4952  scale = ABS(mcfnetwork->arccapacityscales[a])/maxval; /* divide by maxval to normalize rows */
4953  assert(scale > 0.0);
4954 
4955 #ifdef TIEBREAKING
4956  /* get flow on arc for tie breaking */
4957  rowcols = SCIProwGetCols(capacityrow);
4958  rowlen = SCIProwGetNLPNonz(capacityrow);
4959  totalflow = 0.0;
4960  totalcap = 0.0;
4961  SCIPdebugMsg(scip, " row <%s>: \n", SCIProwGetName(capacityrow));
4962 
4963  for( i = 0; i < rowlen; i++ )
4964  {
4965  c = SCIPcolGetLPPos(rowcols[i]);
4966  assert(0 <= c && c < SCIPgetNLPCols(scip));
4967 
4968  SCIPdebugMsg(scip, " col <%s>: %g\n", SCIPvarGetName(SCIPcolGetVar(rowcols[i])), SCIPcolGetPrimsol(rowcols[i]) );
4969  /* sum up flow on arc a*/
4970  if(colcommodity[c] >= 0)
4971  {
4972  SCIPdebugMsg(scip, " flow col <%s>: %g\n", SCIPvarGetName(SCIPcolGetVar(rowcols[i])), REALABS(SCIPcolGetPrimsol(rowcols[i])) );
4973  totalflow += REALABS(SCIPcolGetPrimsol(rowcols[i]));
4974  }
4975  else
4976  {
4977  SCIPdebugMsg(scip, " cap col <%s>: %g\n", SCIPvarGetName(SCIPcolGetVar(rowcols[i])), REALABS(SCIPcolGetPrimsol(rowcols[i])) );
4978  totalcap += REALABS(SCIPcolGetPrimsol(rowcols[i]));
4979  }
4980  }
4981 
4982  SCIPdebugMsg(scip, "cap arc -- slack:%g -- dual:%g -- flow:%g -- cap:%g \n", scale * slack, dualsol/scale, totalflow * scale, totalcap * scale);
4983 #else
4984  SCIPdebugMsg(scip, "cap arc -- slack:%g -- dual:%g1\n", scale * slack, dualsol/scale);
4985 #endif
4986 
4987  /* put the arc weight into a fresh nodepair */
4988  nodepair.weight = scale * slack - ABS(dualsol)/scale;
4989 #ifdef USEFLOWFORTIEBREAKING
4990  if( !SCIPisPositive(scip, nodepair.weight) )
4991  {
4992  nodepair.weight += totalflow * scale;
4993  nodepair.weight = MIN( nodepair.weight, -0.0001);
4994  }
4995 #endif
4996 #ifdef USECAPACITYFORTIEBREAKING
4997  if( !SCIPisPositive(scip, nodepair.weight) )
4998  {
4999  nodepair.weight += totalcap * scale;
5000  nodepair.weight = MIN( nodepair.weight, -0.0001);
5001  }
5002 #endif
5003 
5004  }
5005  else
5006  {
5007  /* uncapacitated arc has infinite slack */
5008  SCIPdebugMsg(scip, "uncap arc ... slack infinite\n");
5009  nodepair.weight = SCIPinfinity(scip);
5010  }
5011 
5012  /* check if nodepair already exists in hash-table */
5013  nodepairptr = (NODEPAIRENTRY*)(SCIPhashtableRetrieve(hashtable, (void*) (&nodepair) ));
5014 
5015  /* if nodepair already exists update its weight */
5016  if( nodepairptr != NULL )
5017  {
5018  /* adapt weight */
5019  SCIPdebugMsg(scip, "nodepair known [%d,%d] -- old weight:%g -- new weight:%g\n", nodepair.node1,nodepair.node2,nodepairptr->weight,
5020  MIN(nodepair.weight, nodepairptr->weight));
5021  nodepairptr->weight = MIN(nodepair.weight, nodepairptr->weight);
5022  }
5023  else
5024  {
5025  /* no such nodepair in current hash table: insert into array and hashtable */
5026  nodepairs = (*nodepairqueue)->nodepairs;
5027 
5028  assert(nnodepairs < mcfnetwork->narcs);
5029  nodepairs[nnodepairs] = nodepair;
5030  SCIP_CALL( SCIPhashtableInsert(hashtable, (void*) (&nodepairs[nnodepairs]) ) );
5031 
5032  SCIPdebugMsg(scip, "new nodepair [%d,%d]-- weight:%g\n", nodepair.node1, nodepair.node2, nodepair.weight);
5033 
5034  nnodepairs++;
5035  }
5036  }
5037 
5038  /* free hash table */
5039  SCIPhashtableFree(&hashtable);
5040 
5041  /* we still have to fill the priority queue */
5042 
5043 #ifdef BETTERWEIGHTFORDEMANDNODES
5044  /* initial weights have been calculated
5045  * we modify them now depending on the demand emanating at the endnodes of nodepairs
5046  */
5047 
5048  /* calculate max and min weight */
5049  maxweight = +1; /* we want maxweight to be positive */
5050  minweight = -1; /* we want minweight to be negative */
5051  nodepairs = (*nodepairqueue)->nodepairs;
5052  for( n = 0; n < nnodepairs; n++ )
5053  {
5054  /* maxweight should not be infinity (uncap arcs have infinity weight)*/
5055  if(!SCIPisInfinity(scip,nodepairs[n].weight))
5056  maxweight = MAX(maxweight, nodepairs[n].weight);
5057 
5058  minweight = MIN(minweight, nodepairs[n].weight);
5059  }
5060 
5061  SCIPdebugMsg(scip, "min/max weight:%g / %g\n", minweight, maxweight);
5062 #endif
5063 
5064  /* initialize priority queue */
5065  SCIP_CALL( SCIPpqueueCreate(&(*nodepairqueue)->pqueue, nnodepairs, 2.0, compNodepairs) );
5066 
5067  /* fill priority queue using array nodepairs */
5068  for( n = 0; n < nnodepairs; n++ )
5069  {
5070  int node1 = nodepairs[n].node1;
5071  int node2 = nodepairs[n].node2;
5072 
5073 #ifdef BETTERWEIGHTFORDEMANDNODES
5074  SCIP_Real rhs = 0;
5075  SCIP_Bool hasdemand1 = FALSE;
5076  SCIP_Bool hasdemand2 = FALSE;
5077 
5078  int k; /* commodity */
5079 
5080  SCIPdebugMsg(scip, "nodepair [%d,%d] weight %g\n", node1,node2,nodepairs[n].weight);
5081  /* check both nodes for their demand value in all commodities
5082  * the demand value can be read from the rhs
5083  * of the flowrows
5084  */
5085  /* node1 */
5086  for( k = 0; k < ncommodities; k++ )
5087  {
5088  if( nodeflowrows[node1][k] == NULL )
5089  continue;
5090 
5091  if( nodeflowscales[node1][k] > 0.0 )
5092  rhs = SCIProwGetRhs(nodeflowrows[node1][k]) - SCIProwGetConstant(nodeflowrows[node1][k]);
5093  else
5094  rhs = SCIProwGetLhs(nodeflowrows[node1][k]) - SCIProwGetConstant(nodeflowrows[node1][k]);
5095 
5096  assert( !SCIPisInfinity(scip,ABS(rhs)) );
5097 
5098  if( ! SCIPisZero(scip, rhs) )
5099  {
5100  hasdemand1 = TRUE;
5101  break;
5102  }
5103 
5104  }
5105  /* node2 */
5106  for( k = 0; k < ncommodities; k++ )
5107  {
5108  if( nodeflowrows[node2][k] == NULL )
5109  continue;
5110 
5111  if( nodeflowscales[node2][k] > 0.0 )
5112  rhs = SCIProwGetRhs(nodeflowrows[node2][k]) - SCIProwGetConstant(nodeflowrows[node2][k]);
5113  else
5114  rhs = SCIProwGetLhs(nodeflowrows[node2][k]) - SCIProwGetConstant(nodeflowrows[node2][k]);
5115 
5116  assert(! SCIPisInfinity(scip, ABS(rhs)));
5117 
5118  if( ! SCIPisZero(scip, rhs) )
5119  {
5120  hasdemand2 = TRUE;
5121  break;
5122  }
5123  }
5124 
5125  /* if one of the nodes has no demand increase the score
5126  * (slack arcs are still shrunk first)
5127  *
5128  */
5129  if( SCIPisPositive(scip, nodepairs[n].weight))
5130  {
5131  assert(SCIPisPositive(scip, maxweight));
5132 
5133  if( !hasdemand1 || !hasdemand2 )
5134  nodepairs[n].weight += maxweight;
5135  }
5136  else
5137  {
5138  assert( SCIPisNegative(scip, minweight));
5139 
5140  if( hasdemand1 && hasdemand2)
5141  nodepairs[n].weight += minweight;
5142  }
5143 #endif
5144  SCIPdebugMsg(scip, "nodepair [%d,%d] weight %g\n", node1,node2,nodepairs[n].weight);
5145 
5146  /* fill priority queue */
5147  SCIP_CALL( SCIPpqueueInsert((*nodepairqueue)->pqueue, (void*)&(*nodepairqueue)->nodepairs[n]) );
5148  }
5149 
5150  return SCIP_OKAY;
5151 }
5152 
5153 
5154 /** frees memory of a nodepair queue */
5155 static
5157  SCIP* scip, /**< SCIP data structure */
5158  NODEPAIRQUEUE** nodepairqueue /**< pointer to nodepair priority queue */
5159  )
5160 {
5161  assert(nodepairqueue != NULL);
5162  assert(*nodepairqueue != NULL);
5163 
5164  SCIPpqueueFree(&(*nodepairqueue)->pqueue);
5165  SCIPfreeBufferArray(scip, &(*nodepairqueue)->nodepairs);
5166  SCIPfreeBuffer(scip, nodepairqueue);
5167 }
5168 
5169 
5170 /** returns whether there are any nodepairs left on the queue */
5171 static
5173  NODEPAIRQUEUE* nodepairqueue /**< nodepair priority queue */
5174  )
5175 {
5176  assert(nodepairqueue != NULL);
5177 
5178  return (SCIPpqueueFirst(nodepairqueue->pqueue) == NULL);
5179 }
5180 
5181 
5182 /** removes the top element from the nodepair priority queue, returns nodepair entry or NULL */
5183 static
5185  NODEPAIRQUEUE* nodepairqueue /**< nodepair priority queue */
5186  )
5187 {
5188 
5189  assert(nodepairqueue != NULL);
5190 
5191  return (NODEPAIRENTRY*)SCIPpqueueRemove(nodepairqueue->pqueue);
5192 }
5193 
5194 /*
5195  * Node partition methods
5196  * used for creating partitions of nodes
5197  * use union find algorithm
5198 */
5199 
5200 /** returns the representative node in the cluster of the given node */
5201 static
5203  NODEPARTITION* nodepartition, /**< node partition data structure */
5204  int v /**< node id to get representative for */
5205  )
5206 {
5207  return unionfindGetRepresentative(nodepartition->representatives, v);
5208 }
5209 
5210 /** joins two clusters given by their representative nodes */
5211 static
5213  NODEPARTITION* nodepartition, /**< node partition data structure */
5214  int rep1, /**< representative of first cluster */
5215  int rep2 /**< representative of second cluster */
5216  )
5217 {
5218  unionfindJoinSets (nodepartition->representatives, rep1, rep2);
5219 }
5220 
5221 /** partitions nodes into a small number of clusters */
5222 static
5224  SCIP* scip, /**< SCIP data structure */
5225  SCIP_MCFNETWORK* mcfnetwork, /**< MCF network structure */
5226  NODEPARTITION** nodepartition, /**< pointer to node partition data structure */
5227  int nclusters /**< number of clusters to generate */
5228  )
5229 {
5230  NODEPAIRQUEUE* nodepairqueue;
5231  int* clustersize;
5232  int nclustersleft;
5233  int v;
5234  int c;
5235  int pos;
5236 
5237  assert(mcfnetwork != NULL);
5238  assert(nodepartition != NULL);
5239  assert(mcfnetwork->nnodes >= 1);
5240 
5241  /* allocate and initialize memory */
5242  SCIP_CALL( SCIPallocBuffer(scip, nodepartition) );
5243  SCIP_CALL( SCIPallocBufferArray(scip, &(*nodepartition)->representatives, mcfnetwork->nnodes) );
5244  SCIP_CALL( SCIPallocBufferArray(scip, &(*nodepartition)->nodeclusters, mcfnetwork->nnodes) );
5245  SCIP_CALL( SCIPallocBufferArray(scip, &(*nodepartition)->clusternodes, mcfnetwork->nnodes) );
5246  SCIP_CALL( SCIPallocBufferArray(scip, &(*nodepartition)->clusterbegin, nclusters+1) );
5247  (*nodepartition)->nclusters = 0;
5248 
5249  /* we start with each node being in its own cluster */
5250  unionfindInitSets((*nodepartition)->representatives, mcfnetwork->nnodes);
5251 
5252  /* create priority queue for nodepairs */
5253  SCIP_CALL( nodepairqueueCreate(scip, mcfnetwork, &nodepairqueue) );
5254 
5255  /* loop over nodepairs in order of their weights */
5256  nclustersleft = mcfnetwork->nnodes;
5257  while( !nodepairqueueIsEmpty(nodepairqueue) && nclustersleft > nclusters )
5258  {
5259  NODEPAIRENTRY* nodepair;
5260  int node1;
5261  int node2;
5262  SCIP_Real weight;
5263  int node1rep;
5264  int node2rep;
5265 
5266  /* get the next nodepair */
5267  nodepair = nodepairqueueRemove(nodepairqueue);
5268  assert(nodepair != NULL);
5269  node1 = nodepair->node1;
5270  node2 = nodepair->node2;
5271  SCIPdebug( weight = nodepair->weight; )
5272 
5273  assert(node1 >= 0 && node1 < mcfnetwork->nnodes);
5274  assert(node2 >= 0 && node2 < mcfnetwork->nnodes);
5275 
5276  /* identify the representatives of the two nodes */
5277  node1rep = nodepartitionGetRepresentative(*nodepartition, node1);
5278  node2rep = nodepartitionGetRepresentative(*nodepartition, node2);
5279  assert(0 <= node1rep && node1rep < mcfnetwork->nnodes);
5280  assert(0 <= node2rep && node2rep < mcfnetwork->nnodes);
5281 
5282  /* there is nothing to do if the two nodes are already in the same cluster */
5283  if( node1rep == node2rep )
5284  continue;
5285 
5286  /* shrink nodepair by joining the two clusters */
5287  SCIPdebugMsg(scip, "shrinking nodepair (%d,%d) with weight %g: join representatives %d and %d\n",
5288  node1, node2, weight, node1rep, node2rep);
5289  nodepartitionJoin(*nodepartition, node1rep, node2rep);
5290  nclustersleft--;
5291  }
5292 
5293  /* node 0 must be a representative due to our join procedure */
5294  assert((*nodepartition)->representatives[0] == 0);
5295 
5296  /* if there have been too few arcs to shrink the graph to the required number of clusters, join clusters with first cluster
5297  * to create a larger disconnected cluster
5298  */
5299  if( nclustersleft > nclusters )
5300  {
5301  for( v = 1; v < mcfnetwork->nnodes && nclustersleft > nclusters; v++ )
5302  {
5303  int rep;
5304 
5305  rep = nodepartitionGetRepresentative(*nodepartition, v);
5306  if( rep != 0 )
5307  {
5308  nodepartitionJoin(*nodepartition, 0, rep);
5309  nclustersleft--;
5310  }
5311  }
5312  }
5313  assert(nclustersleft <= nclusters);
5314 
5315  /* extract the clusters */
5316  SCIP_CALL( SCIPallocBufferArray(scip, &clustersize, nclusters) );
5317  BMSclearMemoryArray(clustersize, nclusters);
5318  for( v = 0; v < mcfnetwork->nnodes; v++ )
5319  {
5320  int rep;
5321 
5322  /* get cluster of node */
5323  rep = nodepartitionGetRepresentative(*nodepartition, v);
5324  assert(rep <= v); /* due to our joining procedure */
5325  if( rep == v )
5326  {
5327  /* node is its own representative: this is a new cluster */
5328  c = (*nodepartition)->nclusters;
5329  (*nodepartition)->nclusters++;
5330  }
5331  else
5332  c = (*nodepartition)->nodeclusters[rep];
5333  assert(0 <= c && c < nclusters);
5334 
5335  /* assign node to cluster */
5336  (*nodepartition)->nodeclusters[v] = c;
5337  clustersize[c]++;
5338  }
5339 
5340  /* fill the clusterbegin array */
5341  pos = 0;
5342  for( c = 0; c < (*nodepartition)->nclusters; c++ )
5343  {
5344  (*nodepartition)->clusterbegin[c] = pos;
5345  pos += clustersize[c];
5346  }
5347  assert(pos == mcfnetwork->nnodes);
5348  (*nodepartition)->clusterbegin[(*nodepartition)->nclusters] = mcfnetwork->nnodes;
5349 
5350  /* fill the clusternodes array */
5351  BMSclearMemoryArray(clustersize, (*nodepartition)->nclusters);
5352  for( v = 0; v < mcfnetwork->nnodes; v++ )
5353  {
5354  c = (*nodepartition)->nodeclusters[v];
5355  assert(0 <= c && c < (*nodepartition)->nclusters);
5356  pos = (*nodepartition)->clusterbegin[c] + clustersize[c];
5357  assert(pos < (*nodepartition)->clusterbegin[c+1]);
5358  (*nodepartition)->clusternodes[pos] = v;
5359  clustersize[c]++;
5360  }
5361 
5362  /* free temporary memory */
5363  SCIPfreeBufferArray(scip, &clustersize);
5364 
5365  /* free nodepair queue */
5366  nodepairqueueFree(scip, &nodepairqueue);
5367 
5368  return SCIP_OKAY;
5369 }
5370 
5371 /** frees node partition data */
5372 static
5374  SCIP* scip, /**< SCIP data structure */
5375  NODEPARTITION** nodepartition /**< pointer to node partition data structure */
5376  )
5377 {
5378  assert(nodepartition != NULL);
5379  assert(*nodepartition != NULL);
5380 
5381  SCIPfreeBufferArray(scip, &(*nodepartition)->clusterbegin);
5382  SCIPfreeBufferArray(scip, &(*nodepartition)->clusternodes);
5383  SCIPfreeBufferArray(scip, &(*nodepartition)->nodeclusters);
5384  SCIPfreeBufferArray(scip, &(*nodepartition)->representatives);
5385  SCIPfreeBuffer(scip, nodepartition);
5386 }
5387 
5388 /** returns whether given node v is in a cluster that belongs to the partition S */
5389 static
5391  NODEPARTITION* nodepartition, /**< node partition data structure, or NULL */
5392  unsigned int partition, /**< partition of nodes, or node number in single-node partition */
5393  SCIP_Bool inverted, /**< should partition be inverted? */
5394  int v /**< node to check */
5395  )
5396 {
5397  /* if the node does not exist, it is not in the partition
5398  * (and also not in the inverted partition)
5399  */
5400  if( v < 0 )
5401  return FALSE;
5402 
5403  if( nodepartition == NULL )
5404  return ((v == (int)partition) == !inverted);
5405  else
5406  {
5407  int cluster;
5408  unsigned int clusterbit;
5409 
5410  cluster = nodepartition->nodeclusters[v];
5411  assert(0 <= cluster && cluster < nodepartition->nclusters);
5412  clusterbit = (1 << cluster); /*lint !e701*/
5413 
5414  return (((partition & clusterbit) != 0) == !inverted);
5415  }
5416 }
5417 
5418 /** returns whether each of the sets S and T defined by the nodepartition is connected */
5419 static
5421  SCIP* scip, /**< SCIP data structure */
5422  SCIP_MCFNETWORK* mcfnetwork, /**< MCF network structure */
5423  NODEPARTITION* nodepartition, /**< node partition data structure */
5424  unsigned int partition /**< partition of nodes, or node number in single-node partition */
5425  )
5426 {
5427  const int* nodeclusters = nodepartition->nodeclusters;
5428  const int* arcsources = mcfnetwork->arcsources;
5429  const int* arctargets = mcfnetwork->arctargets;
5430  int narcs = mcfnetwork->narcs;
5431  int nclusters;
5432 
5433  int ncomponents;
5434  int a;
5435  int* rep;
5436 
5437  assert(nodepartition->nodeclusters != NULL);
5438  nclusters = nodepartition->nclusters;
5439 
5440  if( SCIPallocBufferArray(scip, &rep, nclusters) != SCIP_OKAY )
5441  return 0;
5442 
5443  /* start with each cluster being isolated */
5444  unionfindInitSets(rep, nclusters);
5445  ncomponents = nclusters;
5446  assert(ncomponents >= 2);
5447 
5448  /* for each arc within S or within T join the connected clusters */
5449  for( a = 0; a < narcs; a++ )
5450  {
5451  int s = arcsources[a];
5452  int t = arctargets[a];
5453 
5454  /* ignore arcs that connect the pseudo node -1 */
5455  if( s == -1 || t == -1 )
5456  continue;
5457 
5458  /* check if arc is within one of the components */
5459  if( nodeInPartition(nodepartition, partition, FALSE, s) == nodeInPartition(nodepartition, partition, FALSE, t) )
5460  {
5461  int cs;
5462  int ct;
5463  int repcs;
5464  int repct;
5465 
5466  assert(0 <= s && s < mcfnetwork->nnodes);
5467  assert(0 <= t && t < mcfnetwork->nnodes);
5468 
5469  /* get clusters of the two nodes */
5470  cs = nodeclusters[s];
5471  ct = nodeclusters[t];
5472  assert(0 <= cs && cs < nclusters);
5473  assert(0 <= ct && ct < nclusters);
5474 
5475  /* nothing to do if we are already in the same cluster */
5476  if( cs == ct )
5477  continue;
5478 
5479  /* get representatives of clusters in the union structure */
5480  repcs = unionfindGetRepresentative (rep, cs);
5481  repct = unionfindGetRepresentative (rep, ct);
5482  if( repcs == repct )
5483  continue;
5484 
5485  /* the arc connects two previously unconnected components of S or T */
5486 
5487  /* check if we already reached two distinct components */
5488  ncomponents--;
5489  if( ncomponents <= 2 )
5490  break;
5491 
5492  /* join the two cluster sets and continue */
5493  unionfindJoinSets (rep, repcs, repct);
5494  }
5495  }
5496 
5497  /* each of the two shores S and T are connected if we were able to shrink the graph
5498  into two components */
5499  assert(ncomponents >= 2);
5500  SCIPfreeBufferArray(scip, &rep);
5501  return (ncomponents == 2);
5502 }
5503 
5504 #ifdef SCIP_DEBUG
5505 static
5506 void nodepartitionPrint(
5507  NODEPARTITION* nodepartition /**< node partition data structure */
5508  )
5509 {
5510  int c;
5511 
5512  for( c = 0; c < nodepartition->nclusters; c++ )
5513  {
5514  int i;
5515 
5516  MCFdebugMessage("cluster %d:", c);
5517  for( i = nodepartition->clusterbegin[c]; i < nodepartition->clusterbegin[c+1]; i++ )
5518  MCFdebugMessage(" %d", nodepartition->clusternodes[i]);
5519  MCFdebugMessage("\n");
5520  }
5521 }
5522 #endif
5523 
5524 #ifdef OUTPUTGRAPH
5525 /** generates a GML file to visualize the network graph and LP solution */
5526 static
5527 SCIP_RETCODE outputGraph(
5528  SCIP* scip, /**< SCIP data structure */
5529  SCIP_MCFNETWORK* mcfnetwork, /**< MCF network structure */
5530  NODEPARTITION* nodepartition, /**< node partition data structure, or NULL */
5531  SCIP_Bool inverted, /**< should partition be inverted? */
5532  unsigned int partition /**< partition of nodes, or node number */
5533  )
5534 {
5535  FILE* file;
5536  char filename[SCIP_MAXSTRLEN];
5537  int v;
5538  int a;
5539 
5540  /* open file */
5541  if( nodepartition == NULL )
5542  (void) SCIPsnprintf(filename, SCIP_MAXSTRLEN, "mcf-node-%d.gml", partition);
5543  else
5544  (void) SCIPsnprintf(filename, SCIP_MAXSTRLEN, "mcf-part-%d.gml", partition);
5545  SCIPinfoMessage(scip, NULL, "creating GML output file <%s>...\n", filename);
5546  file = fopen(filename, "w");
5547  if( file == NULL )
5548  {
5549  SCIPerrorMessage("cannot create GML output file <%s>\n", filename);
5550  return SCIP_FILECREATEERROR;
5551  }
5552 
5553  /* print GML header */
5554  fprintf(file, "graph\n");
5555  fprintf(file, "[\n");
5556  fprintf(file, " hierarchic 1\n");
5557  fprintf(file, " label \"\"\n");
5558  fprintf(file, " directed 1\n");
5559 
5560  /* nodes */
5561  for( v = 0; v < mcfnetwork->nnodes; v++ )
5562  {
5563  char label[SCIP_MAXSTRLEN];
5564  SCIP_Bool inpartition;
5565 
5566  if( mcfnetwork->nodeflowrows[v][0] != NULL )
5567  (void) SCIPsnprintf(label, SCIP_MAXSTRLEN, "%s", SCIProwGetName(mcfnetwork->nodeflowrows[v][0]));
5568  else
5569  (void) SCIPsnprintf(label, SCIP_MAXSTRLEN, "%d", v);
5570  inpartition = nodeInPartition(nodepartition, partition, inverted, v);
5571 
5572  fprintf(file, " node\n");
5573  fprintf(file, " [\n");
5574  fprintf(file, " id %d\n", v);
5575  fprintf(file, " label \"%s\"\n", label);
5576  fprintf(file, " graphics\n");
5577  fprintf(file, " [\n");
5578  fprintf(file, " w 30.0\n");
5579  fprintf(file, " h 30.0\n");
5580  fprintf(file, " type \"ellipse\"\n");
5581  if( inpartition )
5582  fprintf(file, " fill \"#FF0000\"\n");
5583  else
5584  fprintf(file, " fill \"#00FF00\"\n");
5585  fprintf(file, " outline \"#000000\"\n");
5586  fprintf(file, " ]\n");
5587  fprintf(file, " LabelGraphics\n");
5588  fprintf(file, " [\n");
5589  fprintf(file, " text \"%s\"\n", label);
5590  fprintf(file, " fontSize 13\n");
5591  fprintf(file, " fontName \"Dialog\"\n");
5592  fprintf(file, " anchor \"c\"\n");
5593  fprintf(file, " ]\n");
5594  if( inpartition )
5595  fprintf(file, " gid %d\n", mcfnetwork->nnodes+1);
5596  else
5597  fprintf(file, " gid %d\n", mcfnetwork->nnodes+2);
5598  fprintf(file, " ]\n");
5599  }
5600 
5601  /* dummy node for missing arc sources or arc targets */
5602  fprintf(file, " node\n");
5603  fprintf(file, " [\n");
5604  fprintf(file, " id %d\n", mcfnetwork->nnodes);
5605  fprintf(file, " label \"?\"\n");
5606  fprintf(file, " graphics\n");
5607  fprintf(file, " [\n");
5608  fprintf(file, " w 30.0\n");
5609  fprintf(file, " h 30.0\n");
5610  fprintf(file, " type \"ellipse\"\n");
5611  fprintf(file, " fill \"#FFFFFF\"\n");
5612  fprintf(file, " outline \"#000000\"\n");
5613  fprintf(file, " ]\n");
5614  fprintf(file, " LabelGraphics\n");
5615  fprintf(file, " [\n");
5616  fprintf(file, " text \"?\"\n");
5617  fprintf(file, " fontSize 13\n");
5618  fprintf(file, " fontName \"Dialog\"\n");
5619  fprintf(file, " anchor \"c\"\n");
5620  fprintf(file, " ]\n");
5621  fprintf(file, " ]\n");
5622 
5623  /* group node for partition S */
5624  fprintf(file, " node\n");
5625  fprintf(file, " [\n");
5626  fprintf(file, " id %d\n", mcfnetwork->nnodes+1);
5627  fprintf(file, " label \"Partition S\"\n");
5628  fprintf(file, " graphics\n");
5629  fprintf(file, " [\n");
5630  fprintf(file, " type \"roundrectangle\"\n");
5631  fprintf(file, " fill \"#CAECFF84\"\n");
5632  fprintf(file, " outline \"#666699\"\n");
5633  fprintf(file, " outlineStyle \"dotted\"\n");
5634  fprintf(file, " topBorderInset 0\n");
5635  fprintf(file, " bottomBorderInset 0\n");
5636  fprintf(file, " leftBorderInset 0\n");
5637  fprintf(file, " rightBorderInset 0\n");
5638  fprintf(file, " ]\n");
5639  fprintf(file, " LabelGraphics\n");
5640  fprintf(file, " [\n");
5641  fprintf(file, " text \"Partition S\"\n");
5642  fprintf(file, " fill \"#99CCFF\"\n");
5643  fprintf(file, " fontSize 15\n");
5644  fprintf(file, " fontName \"Dialog\"\n");
5645  fprintf(file, " alignment \"right\"\n");
5646  fprintf(file, " autoSizePolicy \"node_width\"\n");
5647  fprintf(file, " anchor \"t\"\n");
5648  fprintf(file, " borderDistance 0.0\n");
5649  fprintf(file, " ]\n");
5650  fprintf(file, " isGroup 1\n");
5651  fprintf(file, " ]\n");
5652 
5653  /* group node for partition T */
5654  fprintf(file, " node\n");
5655  fprintf(file, " [\n");
5656  fprintf(file, " id %d\n", mcfnetwork->nnodes+2);
5657  fprintf(file, " label \"Partition T\"\n");
5658  fprintf(file, " graphics\n");
5659  fprintf(file, " [\n");
5660  fprintf(file, " type \"roundrectangle\"\n");
5661  fprintf(file, " fill \"#CAECFF84\"\n");
5662  fprintf(file, " outline \"#666699\"\n");
5663  fprintf(file, " outlineStyle \"dotted\"\n");
5664  fprintf(file, " topBorderInset 0\n");
5665  fprintf(file, " bottomBorderInset 0\n");
5666  fprintf(file, " leftBorderInset 0\n");
5667  fprintf(file, " rightBorderInset 0\n");
5668  fprintf(file, " ]\n");
5669  fprintf(file, " LabelGraphics\n");
5670  fprintf(file, " [\n");
5671  fprintf(file, " text \"Partition T\"\n");
5672  fprintf(file, " fill \"#99CCFF\"\n");
5673  fprintf(file, " fontSize 15\n");
5674  fprintf(file, " fontName \"Dialog\"\n");
5675  fprintf(file, " alignment \"right\"\n");
5676  fprintf(file, " autoSizePolicy \"node_width\"\n");
5677  fprintf(file, " anchor \"t\"\n");
5678  fprintf(file, " borderDistance 0.0\n");
5679  fprintf(file, " ]\n");
5680  fprintf(file, " isGroup 1\n");
5681  fprintf(file, " ]\n");
5682 
5683  /* arcs */
5684  for( a = 0; a < mcfnetwork->narcs; a++ )
5685  {
5686  SCIP_ROW* row;
5687  SCIP_Real slack;
5688  SCIP_Bool hasfractional;
5689  char label[SCIP_MAXSTRLEN];
5690 
5691  if( mcfnetwork->arccapacityrows[a] != NULL )
5692  (void) SCIPsnprintf(label, SCIP_MAXSTRLEN, "%s", SCIProwGetName(mcfnetwork->arccapacityrows[a]));
5693  else
5694  (void) SCIPsnprintf(label, SCIP_MAXSTRLEN, "%d", a);
5695 
5696  hasfractional = FALSE;
5697  row = mcfnetwork->arccapacityrows[a];
5698  if( row != NULL )
5699  {
5700  SCIP_COL** rowcols;
5701  int rowlen;
5702  int i;
5703 
5704  slack = ABS(mcfnetwork->arccapacityscales[a]) * SCIPgetRowLPFeasibility(scip, row);
5705  rowcols = SCIProwGetCols(row);
5706  rowlen = SCIProwGetNLPNonz(row);
5707  for( i = 0; i < rowlen; i++ )
5708  {
5709  SCIP_VAR* var;
5710 
5711  var = SCIPcolGetVar(rowcols[i]);
5712  if( SCIPvarIsIntegral(var) && !SCIPisFeasIntegral(scip, SCIPvarGetLPSol(var)) )
5713  {
5714  hasfractional = TRUE;
5715  break;
5716  }
5717  }
5718  }
5719  else
5720  slack = SCIPinfinity(scip);
5721 
5722  fprintf(file, " edge\n");
5723  fprintf(file, " [\n");
5724  fprintf(file, " source %d\n", mcfnetwork->arcsources[a] >= 0 ? mcfnetwork->arcsources[a] : mcfnetwork->nnodes);
5725  fprintf(file, " target %d\n", mcfnetwork->arctargets[a] >= 0 ? mcfnetwork->arctargets[a] : mcfnetwork->nnodes);
5726  fprintf(file, " label \"%s\"\n", label);
5727  fprintf(file, " graphics\n");
5728  fprintf(file, " [\n");
5729  if( SCIPisFeasPositive(scip, slack) )
5730  fprintf(file, " fill \"#000000\"\n");
5731  else
5732  fprintf(file, " fill \"#FF0000\"\n");
5733  if( hasfractional )
5734  fprintf(file, " style \"dashed\"\n");
5735  fprintf(file, " width 1\n");
5736  fprintf(file, " targetArrow \"standard\"\n");
5737  fprintf(file, " ]\n");
5738  fprintf(file, " LabelGraphics\n");
5739  fprintf(file, " [\n");
5740  fprintf(file, " text \"%s\"\n", label);
5741  fprintf(file, " ]\n");
5742  fprintf(file, " ]\n");
5743  }
5744 
5745  /* print GML footer */
5746  fprintf(file, "]\n");
5747 
5748  /* close file */
5749  fclose(file);
5750 
5751  return SCIP_OKAY;
5752 }
5753 #endif
5754 
5755 /** adds given cut to LP if violated */
5756 static
5758  SCIP* scip, /**< SCIP data structure */
5759  SCIP_SEPA* sepa, /**< separator */
5760  SCIP_SEPADATA* sepadata, /**< separator data */
5761  SCIP_SOL* sol, /**< the solution that should be separated, or NULL for LP solution */
5762  SCIP_Real* cutcoefs, /**< coefficients of active variables in cut */
5763  SCIP_Real cutrhs, /**< right hand side of cut */
5764  SCIP_Bool cutislocal, /**< is the cut only locally valid? */
5765  int cutrank, /**< rank of the cut */
5766  int* ncuts, /**< pointer to count the number of added cuts */
5767  SCIP_Bool* cutoff /**< pointer to store whether a cutoff was found */
5768  )
5769 {
5770  SCIP_ROW* cut;
5771  char cutname[SCIP_MAXSTRLEN];
5772  SCIP_VAR** vars;
5773  int nvars;
5774  int v;
5775 
5776 /* variables for knapsack cover separation */
5777  SCIP_VAR** cutvars = NULL;
5778  SCIP_Real* cutvals = NULL;
5779  int ncutvars;
5780 
5781  assert(scip != NULL);
5782  assert(sepadata != NULL);
5783  assert(cutcoefs != NULL);
5784  assert(ncuts != NULL);
5785  assert(cutoff != NULL);
5786 
5787  /* get active problem variables */
5788  SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
5789  assert(nvars == 0 || vars != NULL);
5790 
5791  ncutvars = 0;
5792  *cutoff = FALSE;
5793 
5794  if( sepadata->separateknapsack )
5795  {
5796  /* allocate temporary memory */
5797  SCIP_CALL( SCIPallocBufferArray(scip, &cutvars, nvars) );
5798  SCIP_CALL( SCIPallocBufferArray(scip, &cutvals, nvars) );
5799  }
5800 
5801  /* create the cut */
5802  (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "mcf%d_%d", SCIPgetNLPs(scip), *ncuts);
5803  SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &cut, sepa, cutname, -SCIPinfinity(scip), cutrhs,
5804  cutislocal, FALSE, sepadata->dynamiccuts) );
5805 
5806  /* add coefficients */
5807  SCIP_CALL( SCIPcacheRowExtensions(scip, cut) );
5808  for( v = 0; v < nvars; v++ )
5809  {
5810  if( SCIPisZero(scip, cutcoefs[v]) )
5811  continue;
5812 
5813  SCIP_CALL( SCIPaddVarToRow(scip, cut, vars[v], cutcoefs[v]) );
5814 
5815  if( sepadata->separateknapsack )
5816  {
5817  assert(cutvars != NULL && cutvals != NULL);
5818  cutvars[ncutvars] = vars[v];
5819  cutvals[ncutvars] = cutcoefs[v];
5820  ncutvars++;
5821  }
5822 
5823  }
5824  SCIP_CALL( SCIPflushRowExtensions(scip, cut) );
5825 
5826  /* set cut rank */
5827  SCIProwChgRank(cut, cutrank);
5828 
5829  /* check efficacy */
5830  if( SCIPisCutEfficacious(scip, sol, cut) )
5831  {
5832  SCIPdebugMsg(scip, " -> found MCF cut <%s>: rhs=%f, act=%f eff=%f rank=%d\n",
5833  cutname, cutrhs, SCIPgetRowSolActivity(scip, cut, sol), SCIPgetCutEfficacy(scip, sol, cut), SCIProwGetRank(cut));
5834  /*SCIPdebug( SCIP_CALL(SCIPprintRow(scip, cut, NULL)) );*/
5835  SCIP_CALL( SCIPaddCut(scip, sol, cut, FALSE, cutoff) );
5836 
5837  if( !(*cutoff) && !cutislocal )
5838  {
5839  SCIP_CALL( SCIPaddPoolCut(scip, cut) );
5840  }
5841  (*ncuts)++;
5842  }
5843 
5844  /* release the row */
5845  SCIP_CALL( SCIPreleaseRow(scip, &cut) );
5846 
5847  if( !(*cutoff) && sepadata->separateknapsack)
5848  {
5849  /* relax cut to knapsack row and separate lifted cover cuts */
5850  SCIP_CALL( SCIPseparateRelaxedKnapsack(scip, NULL, sepa, ncutvars, cutvars, cutvals, +1.0, cutrhs, sol, cutoff, ncuts) );
5851 
5852  /* free temporary memory */
5853  SCIPfreeBufferArray(scip, &cutvals);
5854  SCIPfreeBufferArray(scip, &cutvars);
5855  }
5856 
5857  return SCIP_OKAY;
5858 }
5859 
5860 /** enumerates cuts between subsets of the clusters
5861  * generates single-node cuts if nodepartition == NULL, otherwise generates cluster cuts
5862  */
5863 static
5865  SCIP* scip, /**< SCIP data structure */
5866  SCIP_SEPA* sepa, /**< separator */
5867  SCIP_SEPADATA* sepadata, /**< separator data */
5868  SCIP_SOL* sol, /**< the solution that should be separated, or NULL for LP solution */
5869  SCIP_MCFNETWORK* mcfnetwork, /**< MCF network structure */
5870  NODEPARTITION* nodepartition, /**< node partition data structure, or NULL */
5871  int* ncuts, /**< pointer to count the number of added cuts */
5872  SCIP_Bool* cutoff /**< pointer to store whether a cutoff was found */
5873  )
5874 {
5875  SCIP_ROW*** nodeflowrows = mcfnetwork->nodeflowrows;
5876  SCIP_Real** nodeflowscales = mcfnetwork->nodeflowscales;
5878  SCIP_ROW** arccapacityrows = mcfnetwork->arccapacityrows;
5880  int* colcommodity = mcfnetwork->colcommodity;
5881  int* arcsources = mcfnetwork->arcsources;
5882  int* arctargets = mcfnetwork->arctargets;
5883  int nnodes = mcfnetwork->nnodes;
5884  int narcs = mcfnetwork->narcs;
5885  int ncommodities = mcfnetwork->ncommodities;
5886  SCIP_MCFMODELTYPE modeltype = mcfnetwork->modeltype;
5887  MCFEFFORTLEVEL effortlevel = sepadata->effortlevel;
5888  int maxsepacuts;
5889 
5890  int nrows;
5891  int nvars;
5892 
5893  SCIP_Real* cutcoefs;
5894  SCIP_Real* deltas;
5895  int deltassize;
5896  int ndeltas;
5897  SCIP_Real* rowweights;
5898  SCIP_Real* comcutdemands;
5899  SCIP_Real* comdemands;
5900  unsigned int partition;
5901  unsigned int allpartitions;
5902  unsigned int startpartition;
5903  SCIP_Bool useinverted;
5904  SCIP_Bool inverted;
5905  int maxtestdelta;
5906 
5907  int oldncuts = 0; /* to check success of separation for one nodeset */
5908  *cutoff = FALSE;
5909 
5910  assert( effortlevel == MCFEFFORTLEVEL_AGGRESSIVE || effortlevel == MCFEFFORTLEVEL_DEFAULT );
5911  nrows = SCIPgetNLPRows(scip);
5912  nvars = SCIPgetNVars(scip);
5913 
5914  /* get the maximal number of cuts allowed in a separation round */
5915  if( SCIPgetDepth(scip) == 0 )
5916  maxsepacuts = sepadata->maxsepacutsroot;
5917  else
5918  maxsepacuts = sepadata->maxsepacuts;
5919  if( maxsepacuts < 0 )
5920  maxsepacuts = INT_MAX;
5921  else if( effortlevel == MCFEFFORTLEVEL_AGGRESSIVE )
5922  maxsepacuts *= 2;
5923 
5924  /* get the maximal number of deltas to use for cmir separation */
5925  maxtestdelta = sepadata->maxtestdelta;
5926  if( maxtestdelta <= 0 )
5927  maxtestdelta = INT_MAX;
5928  else if( effortlevel == MCFEFFORTLEVEL_AGGRESSIVE )
5929  maxtestdelta *= 2;
5930 
5931 
5932  /* Our system has the following form:
5933  * (1) \sum_{a \in \delta^+(v)} f_a^k - \sum_{a \in \delta^-(v)} f_a^k <= -d_v^k for all v \in V and k \in K
5934  * (2) \sum_{k \in K} f_a^k - c_a x_a <= 0 for all a \in A
5935  *
5936  * The partitioning yields two clusters:
5937  *
5938  * A^+
5939  * cluster S ------> cluster T
5940  * <------
5941  * A^-
5942  *
5943  * Now we look at all commodities in which we have to route flow from T to S:
5944  * K^+ = {k : d^k_S := sum_{v \in S} d_v^k > 0}
5945  *
5946  * Then, the base constraint of the c-MIR cut is the sum of those flow conservation constraints and the
5947  * capacity constraints for arcs A^-:
5948  *
5949  * sum_{k \in K^+} sum_{v \in S} (sum_{a \in \delta^+(v)} f_a^k - sum_{a \in \delta^-(v)} f_a^k) <= sum_{k \in K^+} sum_{v \in S} -d_v^k
5950  * + sum_{a \in A^-} sum_{k \in K} f_a^k - c_a x_a <= 0
5951  */
5952 
5953  deltassize = 16;
5954  SCIP_CALL( SCIPallocBufferArray(scip, &deltas, deltassize) );
5955  SCIP_CALL( SCIPallocBufferArray(scip, &rowweights, nrows) );
5956  SCIP_CALL( SCIPallocBufferArray(scip, &comcutdemands, ncommodities) );
5957  SCIP_CALL( SCIPallocBufferArray(scip, &comdemands, ncommodities) );
5958  SCIP_CALL( SCIPallocBufferArray(scip, &cutcoefs, nvars) );
5959 
5960  if( nodepartition == NULL )
5961  {
5962  /* loop over all nodes and generate single-node cuts */
5963  startpartition = 0;
5964  allpartitions = (unsigned int) nnodes;
5965  SCIPdebugMsg(scip, "maxtestdelta: %d, maxsepacuts: %d, nnodes: %d \n", maxtestdelta, maxsepacuts, nnodes);
5966  }
5967  else
5968  {
5969  /* loop over all possible partitions of the clusters;
5970  * cluster i is in S iff bit i of 'partition' is 1
5971  */
5972  int nclusters = nodepartition->nclusters;
5973 
5974  assert((unsigned int)nclusters <= 8*sizeof(unsigned int));
5975  SCIPdebugMsg(scip, "maxtestdelta: %d, maxsepacuts: %d, nclusters: %d \n", maxtestdelta, maxsepacuts, nclusters);
5976 
5977  /* We fix the last cluster to belong to partition T. In the undirected case, this is sufficient,
5978  * because S-T is equivalent to T-S. In the directed case, the loop below is conducted two times,
5979  * one with regular S-T and one with inverted partitions.
5980  */
5981  startpartition = 1;
5982  allpartitions = (1 << (nclusters-1)); /*lint !e701*/
5983  }
5984  useinverted = (mcfnetwork->modeltype == SCIP_MCFMODELTYPE_DIRECTED);
5985 
5986  for( partition = startpartition; partition <= allpartitions-1 && !SCIPisStopped(scip) && *ncuts < maxsepacuts && !*cutoff; partition++ )
5987  {
5988  int v;
5989  int a;
5990  int k;
5991  int d;
5992  int nnodesinS;
5993  /* variables for flowcutset separation */
5994  SCIP_Real baserhs;
5995  SCIP_Real bestdelta = 1.0;
5996  SCIP_Real bestrelviolation;
5997  SCIP_Real f0;
5998 
5999 
6000  if( sepadata->checkcutshoreconnectivity )
6001  {
6002  if( nodepartition != NULL && !nodepartitionIsConnected(scip, mcfnetwork, nodepartition, partition ) )
6003  {
6004  /* if S or T are not connected, it is very likely that there is a cut in our cluster partition
6005  that gives dominating inequalities
6006  */
6007  SCIPdebugMsg(scip, "generating cluster cuts for partition 0x%x \n", partition );
6008  SCIPdebugMsg(scip, " -> either shore S or shore T is not connected - skip partition.\n");
6009  continue;
6010  }
6011  }
6012 
6013  for( inverted = FALSE; inverted <= useinverted && !*cutoff; inverted++ )
6014  {
6015 
6016  if( nodepartition == NULL )
6017  {
6018  SCIPdebugMsg(scip, "generating single-node cuts for node %u (inverted: %u)\n", partition, inverted);
6019  }
6020  else
6021  {
6022  SCIPdebugMsg(scip, "generating cluster cuts for partition 0x%x (inverted: %u)\n", partition, inverted);
6023  }
6024 
6025 #ifdef OUTPUTGRAPH
6026  SCIP_CALL( outputGraph(scip, mcfnetwork, nodepartition, inverted, partition) );
6027 #endif
6028  /* Check if S and T are connected via union find algorithm.
6029  * We do not check this for single node cuts. First, this can be expensive since
6030  * in this case the graph still contains all nodes. Second, we may not find a dominating cut,
6031  * because we may not have the corresponding dominating node set in our cluster partition.
6032  */
6033 
6034  /* clear memory */
6035  BMSclearMemoryArray(rowweights, nrows);
6036  BMSclearMemoryArray(comcutdemands, ncommodities);
6037  BMSclearMemoryArray(comdemands, ncommodities);
6038 
6039  baserhs = 0.0;
6040  ndeltas = 0;
6041 
6042  /* Identify commodities with positive T -> S demand */
6043  nnodesinS = 0;
6044  for( v = 0; v < nnodes; v++ )
6045  {
6046  /* check if node belongs to S */
6047  if( !nodeInPartition(nodepartition, partition, inverted, v) )
6048  {
6049  /* node does not belong to S */
6050  continue;
6051  }
6052  nnodesinS++;
6053  /* update commodity demand */
6054  for( k = 0; k < ncommodities; k++ )
6055  {
6056  SCIP_Real rhs;
6057 
6058  if( nodeflowrows[v][k] == NULL )
6059  continue;
6060 
6061  if( nodeflowscales[v][k] > 0.0 )
6062  rhs = SCIProwGetRhs(nodeflowrows[v][k]) - SCIProwGetConstant(nodeflowrows[v][k]);
6063  else
6064  rhs = SCIProwGetLhs(nodeflowrows[v][k]) - SCIProwGetConstant(nodeflowrows[v][k]);
6065  if( nodeflowinverted[v][k] )
6066  rhs *= -1.0;
6067 
6068  comcutdemands[k] += rhs * nodeflowscales[v][k];
6069  }
6070  }
6071  assert (1 <= nnodesinS && nnodesinS <= nnodes-1);
6072 
6073  /* ignore cuts with only a single node in S or in T, since these have
6074  * already been tried as single node cuts
6075  */
6076  if( sepadata->separatesinglenodecuts && nodepartition != NULL && (nnodesinS == 1 || nnodesinS == nnodes-1) )
6077  {
6078  SCIPdebugMsg(scip, " -> shore S or T has only one node - skip partition.\n");
6079  break;
6080  }
6081 
6082  /* check if there is at least one useful commodity */
6083  if( modeltype == SCIP_MCFMODELTYPE_DIRECTED )
6084  {
6085  for( k = 0; k < ncommodities; k++ )
6086  {
6087  /* in the directed case, use commodities with positive demand (negative -d_k) */
6088  SCIPdebugMsg(scip, " -> commodity %d: directed cutdemand=%g\n", k, comcutdemands[k]);
6089  if( SCIPisNegative(scip, comcutdemands[k]) )
6090  break;
6091  }
6092  }
6093  else
6094  {
6095  for( k = 0; k < ncommodities; k++ )
6096  {
6097  /* in the undirected case, use commodities with non-zero demand */
6098  SCIPdebugMsg(scip, " -> commodity %d: undirected cutdemand=%g\n", k, comcutdemands[k]);
6099  if( !SCIPisZero(scip, comcutdemands[k]) )
6100  break;
6101  }
6102  }
6103  if( k == ncommodities )
6104  continue;
6105 
6106  /* set weights of capacity rows that go from T to S, i.e., a \in A^- */
6107  for( a = 0; a < narcs; a++ )
6108  {
6109  SCIP_COL** rowcols;
6110  SCIP_Real* rowvals;
6111  SCIP_Real feasibility;
6112  int rowlen;
6113  int r;
6114  int j;
6115 
6116  assert(arcsources[a] < nnodes);
6117  assert(arctargets[a] < nnodes);
6118 
6119  /* check if this is an arc of our cut */
6120  if( modeltype == SCIP_MCFMODELTYPE_DIRECTED )
6121  {
6122  /* in the directed case, check if arc goes from T to S */
6123  if( nodeInPartition(nodepartition, partition, inverted, arcsources[a]) ||
6124  !nodeInPartition(nodepartition, partition, inverted, arctargets[a]) )
6125  {
6126  /* arc has source in S or target in T */
6127  continue;
6128  }
6129  }
6130  else
6131  {
6132  /* in the undirected case, check if the arc has endpoints in S and T */
6133  if( nodeInPartition(nodepartition, partition, inverted, arcsources[a]) &&
6134  nodeInPartition(nodepartition, partition, inverted, arctargets[a]) )
6135  {
6136  /* both endpoints are in S */
6137  continue;
6138  }
6139  if( !nodeInPartition(nodepartition, partition, inverted, arcsources[a]) &&
6140  !nodeInPartition(nodepartition, partition, inverted, arctargets[a]) )
6141  {
6142  /* both endpoints are in T */
6143  continue;
6144  }
6145  }
6146 
6147  /* arc might be uncapacitated */
6148  if( arccapacityrows[a] == NULL )
6149  continue;
6150 
6151  /* use capacity row in c-MIR cut */
6152  r = SCIProwGetLPPos(arccapacityrows[a]);
6153  assert(r < nrows);
6154  if( r == -1 ) /* row might have been removed from LP in the meantime */
6155  continue;
6156  assert(rowweights[r] == 0.0);
6157 
6158  /* if one of the arc nodes is unknown, we only use the capacity row if it does not have slack,
6159  * otherwise, we discard it if the slack is too large
6160  */
6161  feasibility = SCIPgetRowSolFeasibility(scip, arccapacityrows[a], sol);
6162  if( arcsources[a] == -1 || arctargets[a] == -1 )
6163  {
6164  if( SCIPisFeasPositive(scip, feasibility) )
6165  continue;
6166  }
6167  else
6168  {
6169  SCIP_Real maxcoef;
6170 
6171  maxcoef = SCIPgetRowMaxCoef(scip, arccapacityrows[a]);
6172  assert(maxcoef > 0.0);
6173  if( SCIPisFeasGT(scip, feasibility/maxcoef, MAXCAPACITYSLACK) )
6174  continue;
6175  }
6176 
6177  rowweights[r] = arccapacityscales[a];
6178  SCIPdebugMsg(scip, " -> arc %d, r=%d, capacity row <%s>: weight=%g slack=%g dual=%g\n", a, r, SCIProwGetName(arccapacityrows[a]), rowweights[r],
6179  SCIPgetRowFeasibility(scip, arccapacityrows[a]), SCIProwGetDualsol(arccapacityrows[a]));
6180  /*SCIPdebug( SCIP_CALL(SCIPprintRow(scip, arccapacityrows[a], NULL)) );*/
6181 
6182  if( sepadata->separateflowcutset )
6183  {
6184  if( rowweights[r] > 0.0 )
6185  baserhs += rowweights[r] * (SCIProwGetRhs(arccapacityrows[a]) - SCIProwGetConstant(arccapacityrows[a]));
6186  else
6187  baserhs += rowweights[r] * (SCIProwGetLhs(arccapacityrows[a]) - SCIProwGetConstant(arccapacityrows[a]));
6188  }
6189 
6190  /* extract useful deltas for c-MIR scaling and update the demand value for commodities (in binary flow model) */
6191  rowcols = SCIProwGetCols(arccapacityrows[a]);
6192  rowvals = SCIProwGetVals(arccapacityrows[a]);
6193  rowlen = SCIProwGetNLPNonz(arccapacityrows[a]);
6194  for( j = 0; j < rowlen; j++ )
6195  {
6196  SCIP_Real coef;
6197  int c;
6198 
6199  coef = rowvals[j] * arccapacityscales[a];
6200  coef = ABS(coef);
6201 
6202  /* update commodity demands */
6203  c = SCIPcolGetLPPos(rowcols[j]);
6204  assert(0 <= c && c < SCIPgetNLPCols(scip));
6205  k = colcommodity[c];
6206  if( k >= 0 )
6207  comdemands[k] = coef;
6208 
6209  /* insert coefficients of integer variables into deltas array */
6210  /* coefficients should not be too small and not be too big */
6211  if( !SCIPisZero(scip, 1.0 / coef) && !SCIPisFeasZero(scip, coef) && SCIPcolIsIntegral(rowcols[j]) )
6212  {
6213  SCIP_Bool exists;
6214  int left;
6215  int right;
6216 
6217  /* binary search if we already know this coefficient */
6218  exists = FALSE;
6219  left = 0;
6220  right = ndeltas-1;
6221  while( left <= right )
6222  {
6223  int mid = (left+right)/2;
6224  /* take deltas that are not too close */
6225  if( REALABS( deltas[mid] / coef - 1.0 ) < 1e-03 ) /*lint !e771*/
6226  {
6227  exists = TRUE;
6228  break;
6229  }
6230  else if( coef < deltas[mid] )
6231  right = mid-1;
6232  else
6233  left = mid+1;
6234  }
6235 
6236  /* insert new candidate value */
6237  if( !exists )
6238  {
6239  assert(right == left-1);
6240  assert(ndeltas <= deltassize);
6241  if( ndeltas == deltassize )
6242  {
6243  deltassize *= 2;
6244  SCIP_CALL( SCIPreallocBufferArray(scip, &deltas, deltassize) );
6245  }
6246  if( left < ndeltas )
6247  {
6248  for( d = ndeltas; d > left; d-- )
6249  deltas[d] = deltas[d-1];
6250  }
6251  deltas[left] = coef;
6252  ndeltas++;
6253  SCIPdebugMsg(scip, " -> new capacity %g considered as delta\n", coef);
6254  }
6255  }
6256  }
6257  }
6258 
6259  /* set weights of node flow conservation constraints in c-MIR aggregation */
6260  for( v = 0; v < nnodes; v++ )
6261  {
6262  /* aggregate flow conservation constraints of the 'active' commodities */
6263  for( k = 0; k < ncommodities; k++ )
6264  {
6265  SCIP_Real scale;
6266  int r;
6267 
6268  /* if commodity was not hit by the capacity constraints of the cut in the graph, ignore the commodity */
6269  if( comdemands[k] == 0.0 )
6270  continue;
6271 
6272  scale = comdemands[k];
6273  if( modeltype == SCIP_MCFMODELTYPE_DIRECTED )
6274  {
6275  /* in the directed case, use rows of commodities with positive demand (negative -d_k) */
6276  if( !SCIPisNegative(scip, comcutdemands[k]) )
6277  continue;
6278 
6279  /* check if node belongs to S */
6280  if( !nodeInPartition(nodepartition, partition, inverted, v) )
6281  {
6282  /* node belongs to T */
6283  continue;
6284  }
6285  }
6286  else
6287  {
6288  /* in the undirected case, use rows of commodities with non-zero demand */
6289  if( SCIPisZero(scip, comcutdemands[k]) )
6290  continue;
6291 
6292  /* If the demand (-d_k) is negative (i.e., points into the wrong direction), we use the flow
6293  * in the opposite direction, i.e., sum over all nodes in T instead of S.
6294  */
6295  if( comcutdemands[k] > 0.0 )
6296  {
6297  /* check if node belongs to T */
6298  if( nodeInPartition(nodepartition, partition, inverted, v) )
6299  {
6300  /* node b