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