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-2018 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  SCIP_CALL( extractFlow(scip, &mcfdata, MAXFLOWVARFLOWROWRATIO, &failed) );
4333  assert(mcfdata.plusflow != NULL);
4334  assert(mcfdata.minusflow != NULL);
4335  assert(mcfdata.colcommodity != NULL);
4336  assert(mcfdata.rowcommodity != NULL);
4337  assert(mcfdata.newcols != NULL);
4338  }
4339 
4340  if( !failed )
4341  {
4342 #ifdef SCIP_DEBUG
4343  printCommodities(scip, &mcfdata);
4344 #endif
4345 
4346  /* 4. identify candidate rows for capacity constraints in the LP
4347  * 5. sort capacity constraint candidates by a ranking on how sure we are that it is indeed a constraint of the desired type
4348  */
4349  SCIP_CALL( extractCapacityRows(scip, &mcfdata) );
4350  assert(mcfdata.capacityrowsigns != NULL);
4351  assert(mcfdata.capacityrowscores != NULL);
4352  assert(mcfdata.capacitycands != NULL);
4353 
4354  if( mcfdata.ncapacitycands == 0 )
4355  failed = TRUE;
4356  }
4357 
4358  if( !failed )
4359  {
4360  /* 6. arc-detection -- identify capacity constraints for the arcs and assign arc ids to columns and capacity constraints */
4361  SCIP_CALL( extractCapacities(scip, &mcfdata) );
4362  assert(mcfdata.colarcid != NULL);
4363  assert(mcfdata.rowarcid != NULL);
4364 
4365  /* 7. node-detection -- assign node ids to flow conservation constraints */
4366  SCIP_CALL( extractNodes(scip, &mcfdata) );
4367  assert(mcfdata.rownodeid != NULL);
4368  assert(mcfdata.colisincident != NULL);
4369  assert(mcfdata.zeroarcarray != NULL);
4370 
4371  /* 8. postprocessing */
4372  /* 8.a if there are still undecided commodity signs, fix them to +1 */
4373  fixCommoditySigns(scip, &mcfdata);
4374 
4375  /* 8.b clean up the network: get rid of commodities without arcs or with at most one node */
4376  SCIP_CALL( cleanupNetwork(scip, &mcfdata) );
4377 
4378  /* 8.c construct incidence function -- assign source and target nodes to capacitated arcs */
4379  SCIP_CALL( identifySourcesTargets(scip, &mcfdata, sepadata, effortlevel) );
4380  assert(mcfdata.arcsources != NULL);
4381  assert(mcfdata.arctargets != NULL);
4382  assert(mcfdata.colsources != NULL);
4383  assert(mcfdata.coltargets != NULL);
4384  assert(mcfdata.firstoutarcs != NULL);
4385  assert(mcfdata.firstinarcs != NULL);
4386  assert(mcfdata.nextoutarcs != NULL);
4387  assert(mcfdata.nextinarcs != NULL);
4388  }
4389 
4390  if( !failed && *effortlevel != MCFEFFORTLEVEL_OFF)
4391  {
4392  int* nodevisited;
4393  int* compnodeid;
4394  int* compnodes;
4395  int* comparcs;
4396  int minnodes;
4397  int v;
4398 
4399  /* 8.d find uncapacitated arcs */
4400  SCIP_CALL( findUncapacitatedArcs(scip, &mcfdata) );
4401 
4402 #ifdef SCIP_DEBUG
4403  printCommodities(scip, &mcfdata);
4404 #endif
4405 
4406  minnodes = MINNODES;
4407 
4408  /* allocate temporary memory for component finding */
4409  SCIP_CALL( SCIPallocBufferArray(scip, &nodevisited, mcfdata.nnodes) );
4410  SCIP_CALL( SCIPallocBufferArray(scip, &compnodes, mcfdata.nnodes) );
4411  SCIP_CALL( SCIPallocBufferArray(scip, &comparcs, mcfdata.narcs) );
4412  BMSclearMemoryArray(nodevisited, mcfdata.nnodes);
4413 
4414  /* allocate temporary memory for v -> compv mapping */
4415  SCIP_CALL( SCIPallocBufferArray(scip, &compnodeid, mcfdata.nnodes) );
4416  for( v = 0; v < mcfdata.nnodes; v++ )
4417  compnodeid[v] = -1;
4418 
4419  /* search components and create a network structure for each of them */
4420  for( v = 0; v < mcfdata.nnodes; v++ )
4421  {
4422  int ncompnodes;
4423  int ncomparcs;
4424 
4425  /* ignore nodes that have been already assigned to a component */
4426  assert(nodevisited[v] == UNKNOWN || nodevisited[v] == VISITED);
4427  if( nodevisited[v] == VISITED )
4428  continue;
4429 
4430  /* identify nodes and arcs of this component */
4431  SCIP_CALL( identifyComponent(scip, &mcfdata, nodevisited, v, compnodes, &ncompnodes, comparcs, &ncomparcs) );
4432  assert(ncompnodes >= 1);
4433  assert(compnodes[0] == v);
4434  assert(nodevisited[v] == VISITED);
4435 
4436  /* ignore network component if it is trivial */
4437  if( ncompnodes >= minnodes && ncomparcs >= MINARCS )
4438  {
4439  SCIP_MCFNETWORK* mcfnetwork;
4440  int i;
4441 
4442  /* make sure that we have enough memory for the new network pointer */
4443  assert(*nmcfnetworks <= MAXNETWORKS);
4444  assert(*nmcfnetworks <= mcfnetworkssize);
4445  if( *nmcfnetworks == mcfnetworkssize )
4446  {
4447  mcfnetworkssize = MAX(2*mcfnetworkssize, *nmcfnetworks+1);
4448  SCIP_CALL( SCIPreallocMemoryArray(scip, mcfnetworks, mcfnetworkssize) );
4449  }
4450  assert(*nmcfnetworks < mcfnetworkssize);
4451 
4452  /* create network data structure */
4453  SCIP_CALL( mcfnetworkCreate(scip, &mcfnetwork) );
4454 
4455  /* fill sparse network structure */
4456  SCIP_CALL( mcfnetworkFill(scip, mcfnetwork, &mcfdata, compnodeid, compnodes, ncompnodes, comparcs, ncomparcs) );
4457 
4458  /* insert in sorted network list */
4459  for( i = *nmcfnetworks; i > 0 && mcfnetwork->nnodes > (*mcfnetworks)[i-1]->nnodes; i-- )
4460  (*mcfnetworks)[i] = (*mcfnetworks)[i-1];
4461  (*mcfnetworks)[i] = mcfnetwork;
4462  (*nmcfnetworks)++;
4463 
4464  /* if we reached the maximal number of networks, update minnodes */
4465  if( *nmcfnetworks >= MAXNETWORKS )
4466  minnodes = MAX(minnodes, (*mcfnetworks)[*nmcfnetworks-1]->nnodes);
4467 
4468  /* if we exceeded the maximal number of networks, delete the last one */
4469  if( *nmcfnetworks > MAXNETWORKS )
4470  {
4471  SCIPdebugMsg(scip, " -> discarded network with %d nodes and %d arcs due to maxnetworks (minnodes=%d)\n",
4472  (*mcfnetworks)[*nmcfnetworks-1]->nnodes, (*mcfnetworks)[*nmcfnetworks-1]->narcs, minnodes);
4473  SCIP_CALL( mcfnetworkFree(scip, &(*mcfnetworks)[*nmcfnetworks-1]) );
4474  (*nmcfnetworks)--;
4475  }
4476  assert(*nmcfnetworks <= MAXNETWORKS);
4477  }
4478  else
4479  {
4480  SCIPdebugMsg(scip, " -> discarded component with %d nodes and %d arcs\n", ncompnodes, ncomparcs);
4481  }
4482  }
4483 
4484  /* free temporary memory */
4485  SCIPfreeBufferArray(scip, &compnodeid);
4486  SCIPfreeBufferArray(scip, &comparcs);
4487  SCIPfreeBufferArray(scip, &compnodes);
4488  SCIPfreeBufferArray(scip, &nodevisited);
4489  }
4490 
4491  /* free memory */
4492  SCIPfreeMemoryArrayNull(scip, &mcfdata.arcsources);
4493  SCIPfreeMemoryArrayNull(scip, &mcfdata.arctargets);
4494  SCIPfreeMemoryArrayNull(scip, &mcfdata.colsources);
4495  SCIPfreeMemoryArrayNull(scip, &mcfdata.coltargets);
4496  SCIPfreeMemoryArrayNull(scip, &mcfdata.firstoutarcs);
4497  SCIPfreeMemoryArrayNull(scip, &mcfdata.firstinarcs);
4498  SCIPfreeMemoryArrayNull(scip, &mcfdata.nextoutarcs);
4499  SCIPfreeMemoryArrayNull(scip, &mcfdata.nextinarcs);
4500  SCIPfreeMemoryArrayNull(scip, &mcfdata.zeroarcarray);
4501  SCIPfreeMemoryArrayNull(scip, &mcfdata.colisincident);
4502  SCIPfreeMemoryArrayNull(scip, &mcfdata.capacityrows);
4503  SCIPfreeMemoryArrayNull(scip, &mcfdata.rownodeid);
4504  SCIPfreeMemoryArrayNull(scip, &mcfdata.rowarcid);
4505  SCIPfreeMemoryArrayNull(scip, &mcfdata.colarcid);
4506  SCIPfreeMemoryArrayNull(scip, &mcfdata.newcols);
4507  SCIPfreeMemoryArrayNull(scip, &mcfdata.rowcommodity);
4508  SCIPfreeMemoryArrayNull(scip, &mcfdata.colcommodity);
4509  SCIPfreeMemoryArrayNull(scip, &mcfdata.commoditysigns);
4510  SCIPfreeMemoryArrayNull(scip, &mcfdata.minusflow);
4511  SCIPfreeMemoryArrayNull(scip, &mcfdata.plusflow);
4512  SCIPfreeMemoryArrayNull(scip, &mcfdata.capacitycands);
4513  SCIPfreeMemoryArrayNull(scip, &mcfdata.flowcands);
4514  SCIPfreeMemoryArrayNull(scip, &mcfdata.capacityrowscores);
4515  SCIPfreeMemoryArrayNull(scip, &mcfdata.capacityrowsigns);
4516  SCIPfreeMemoryArrayNull(scip, &mcfdata.flowrowscores);
4517  SCIPfreeMemoryArrayNull(scip, &mcfdata.flowrowscalars);
4518  SCIPfreeMemoryArrayNull(scip, &mcfdata.flowrowsigns);
4519 
4520  return SCIP_OKAY;
4521 }
4522 #ifdef COUNTNETWORKVARIABLETYPES
4523 /** extracts MCF network structures from the current LP */
4524 static
4525 SCIP_RETCODE printFlowSystemInfo(
4526  SCIP* scip, /**< SCIP data structure */
4527  SCIP_MCFNETWORK** mcfnetworks, /**< array of MCF network structures */
4528  int nmcfnetworks /**< number of MCF networks */
4529  )
4530 {
4531  SCIP_ROW** rows;
4532  SCIP_COL** cols;
4533  SCIP_Bool* colvisited;
4534  int nrows;
4535  int ncols;
4536  int m;
4537  int c;
4538  int a;
4539  int k;
4540  int v;
4541  int nflowrows = 0;
4542  int ncaprows = 0;
4543  int nflowvars = 0;
4544  int nintflowvars = 0;
4545  int nbinflowvars = 0;
4546  int ncontflowvars = 0;
4547  int ncapvars = 0;
4548  int nintcapvars = 0;
4549  int nbincapvars = 0;
4550  int ncontcapvars = 0;
4551 
4552  /* get LP data */
4553  SCIP_CALL( SCIPgetLPRowsData(scip, &rows, &nrows) );
4554  SCIP_CALL( SCIPgetLPColsData(scip, &cols, &ncols) );
4555  SCIP_CALL( SCIPallocBufferArray(scip, &colvisited, ncols) );
4556 
4557  /* get flow variable types */
4558  for(c=0; c < ncols; c++)
4559  colvisited[c]=FALSE;
4560 
4561  MCFdebugMessage("\n\n****** VAR COUNTING ********* \n");
4562 
4563  for(m=0; m < nmcfnetworks; m++)
4564  {
4565  SCIP_MCFNETWORK* mcfnetwork = mcfnetworks[m];
4566 
4567  int narcs = mcfnetwork->narcs;
4568  int nnodes = mcfnetwork->nnodes;
4569  int ncommodities = mcfnetwork->ncommodities;
4570  SCIP_ROW** arccapacityrows = mcfnetwork->arccapacityrows;
4571  SCIP_ROW*** nodeflowrows = mcfnetwork->nodeflowrows;
4572  int* colcommodity = mcfnetwork->colcommodity;
4573 
4574  /* get flow variable types */
4575  for(c=0; c < ncols; c++)
4576  {
4577  SCIP_COL* col;
4578 
4579  if(colcommodity[c] >= 0 && ! colvisited[c])
4580  {
4581  /* this is a flow variable */
4582  nflowvars++;
4583  col = cols[c];
4584  colvisited[c] = TRUE;
4585  switch( SCIPvarGetType(SCIPcolGetVar(col)) )
4586  {
4587  case SCIP_VARTYPE_BINARY:
4588  nbinflowvars++;
4589  break;
4590  case SCIP_VARTYPE_INTEGER:
4591  nintflowvars++;
4592  break;
4593  case SCIP_VARTYPE_IMPLINT:
4594  nintflowvars++;
4595  break;
4597  ncontflowvars++;
4598  break;
4599  default:
4600  SCIPerrorMessage("unknown variable type\n");
4601  SCIPABORT();
4602  return SCIP_INVALIDDATA; /*lint !e527*/
4603  }
4604  }
4605  }
4606  /* get capacity variable types and number of capacity rows*/
4607  for( a = 0; a < narcs; a++ )
4608  {
4609  SCIP_ROW* row;
4610  row = arccapacityrows[a];
4611 
4612  if( row != NULL )
4613  {
4614  SCIP_COL** rowcols;
4615  int rowlen;
4616  int i;
4617 
4618  ncaprows++;
4619  rowcols = SCIProwGetCols(row);
4620  rowlen = SCIProwGetNLPNonz(row);
4621 
4622  for( i = 0; i < rowlen; i++ )
4623  {
4624  c = SCIPcolGetLPPos(rowcols[i]);
4625  assert(0 <= c && c < SCIPgetNLPCols(scip));
4626 
4627  if(colcommodity[c] == -1 && ! colvisited[c] )
4628  {
4629  ncapvars++;
4630  colvisited[c] = TRUE;
4631  switch( SCIPvarGetType(SCIPcolGetVar(rowcols[i]) ) )
4632  {
4633  case SCIP_VARTYPE_BINARY:
4634  nbincapvars++;
4635  break;
4636  case SCIP_VARTYPE_INTEGER:
4637  nintcapvars++;
4638  break;
4639  case SCIP_VARTYPE_IMPLINT:
4640  nintcapvars++;
4641  break;
4643  ncontcapvars++;
4644  break;
4645  default:
4646  SCIPerrorMessage("unknown variable type\n");
4647  SCIPABORT();
4648  return SCIP_INVALIDDATA; /*lint !e527*/
4649  }
4650  }
4651  }
4652  }
4653  }
4654  /* get number of flow rows */
4655  for( k = 0; k < ncommodities; k++ )
4656  {
4657  for( v = 0; v < nnodes; v++ )
4658  {
4659  SCIP_ROW* row;
4660  row = nodeflowrows[v][k];
4661 
4662  if( row != NULL )
4663  nflowrows++;
4664  }
4665  }
4666 
4667  MCFdebugMessage("----- network %i -----\n",m);
4668  MCFdebugMessage(" nof flowrows: %5d\n", nflowrows);
4669  MCFdebugMessage(" nof caprows: %5d\n", ncaprows);
4670  MCFdebugMessage(" nof flowvars: %5d of which [ %d , %d , %d ] are continuous, integer, binary\n",
4671  nflowvars, ncontflowvars, nintflowvars, nbinflowvars);
4672  MCFdebugMessage(" nof capvars: %5d of which [ %d , %d , %d ] are continuous, integer, binary\n",
4673  ncapvars, ncontcapvars, nintcapvars, nbincapvars);
4674  }
4675 
4676  MCFdebugMessage("****** END VAR COUNTING ********* \n\n");
4677 
4678  SCIPfreeBufferArray(scip, &colvisited);
4679 
4680  return SCIP_OKAY;
4681 }
4682 #endif
4683 /*
4684  * Union find methods
4685  * used for generating partitions of node sets and
4686  * for checking connectivity of cut shores
4687  */
4688 
4689 /** initializes a union find data structure by putting each element into its own set */
4690 static
4692  int* representatives, /**< mapping an element v to its representative */
4693  int nelems /**< number of elements in the ground set */
4694  )
4695 {
4696  int v;
4697 
4698  /* we start with each element being in its own set */
4699  for( v = 0; v < nelems; v++ )
4700  representatives[v] = v;
4701 }
4702 
4703 /** applies a union find algorithm to get the representative of v */
4704 static
4706  int* representatives, /**< mapping an element v to its representative */
4707  int v /**< element v to get a representative for */
4708  )
4709 {
4710  assert(representatives != NULL);
4711 
4712  while( v != representatives[v] )
4713  {
4714  representatives[v] = representatives[representatives[v]];
4715  v = representatives[v];
4716  }
4717 
4718  return v;
4719 }
4720 
4721 /** joins two sets in the union find framework */
4722 static
4724  int* representatives, /**< mapping an element v to its representative */
4725  int rep1, /**< representative of first set */
4726  int rep2 /**< representative of second set */
4727  )
4728 {
4729  assert(rep1 != rep2);
4730  assert(representatives[rep1] == rep1);
4731  assert(representatives[rep2] == rep2);
4732 
4733  /* make sure that the smaller representative survives
4734  * -> element 0 is always a representative
4735  */
4736  if( rep1 < rep2 )
4737  representatives[rep2] = rep1;
4738  else
4739  representatives[rep1] = rep2;
4740 }
4741 
4742 /*
4743  * Node pair methods
4744  * used for shrinking the network based on nodepair-weights
4745  * -> creating partition
4746 */
4747 
4748 /** comparison method for weighted nodepairs */
4749 static
4751 {
4752  NODEPAIRENTRY* nodepair1 = (NODEPAIRENTRY*)elem1;
4753  NODEPAIRENTRY* nodepair2 = (NODEPAIRENTRY*)elem2;
4754 
4755  if( nodepair1->weight > nodepair2->weight )
4756  return -1;
4757  else if( nodepair1->weight < nodepair2->weight )
4758  return +1;
4759  else
4760  return 0;
4761 }
4762 
4763 /** NodePair HashTable
4764  * gets the key of the given element */
4765 static
4766 SCIP_DECL_HASHGETKEY(hashGetKeyNodepairs)
4767 {
4768  /*lint --e{715}*/
4769  /* the key is the element itself */
4770  return elem;
4771 }
4772 
4773 /** NodePair HashTable
4774  * returns TRUE iff both keys are equal;
4775  * two nodepairs are equal if both nodes equal
4776  */
4777 static
4778 SCIP_DECL_HASHKEYEQ(hashKeyEqNodepairs)
4779 {
4780 #ifndef NDEBUG
4781  SCIP_MCFNETWORK* mcfnetwork;
4782 #endif
4783  NODEPAIRENTRY* nodepair1;
4784  NODEPAIRENTRY* nodepair2;
4785  int source1;
4786  int source2;
4787  int target1;
4788  int target2;
4789 
4790 #ifndef NDEBUG
4791  mcfnetwork = (SCIP_MCFNETWORK*)userptr;
4792  assert(mcfnetwork != NULL);
4793 #endif
4794 
4795  nodepair1 = (NODEPAIRENTRY*)key1;
4796  nodepair2 = (NODEPAIRENTRY*)key2;
4797 
4798  assert(nodepair1 != NULL);
4799  assert(nodepair2 != NULL);
4800 
4801  source1 = nodepair1->node1;
4802  source2 = nodepair2->node1;
4803  target1 = nodepair1->node2;
4804  target2 = nodepair2->node2;
4805 
4806  assert(source1 >=0 && source1 < mcfnetwork->nnodes);
4807  assert(source2 >=0 && source2 < mcfnetwork->nnodes);
4808  assert(target1 >=0 && target1 < mcfnetwork->nnodes);
4809  assert(target2 >=0 && target2 < mcfnetwork->nnodes);
4810  assert(source1 <= target1);
4811  assert(source2 <= target2);
4812 
4813  return (source1 == source2 && target1 == target2);
4814 }
4815 
4816 /** NodePair HashTable
4817  * returns the hash value of the key */
4818 static
4819 SCIP_DECL_HASHKEYVAL(hashKeyValNodepairs)
4820 {
4821 #ifndef NDEBUG
4822  SCIP_MCFNETWORK* mcfnetwork;
4823 #endif
4824  NODEPAIRENTRY* nodepair;
4825  int source;
4826  int target;
4827  unsigned int hashval;
4828 
4829 #ifndef NDEBUG
4830  mcfnetwork = (SCIP_MCFNETWORK*)userptr;
4831  assert(mcfnetwork != NULL);
4832 #endif
4833 
4834  nodepair = (NODEPAIRENTRY*)key;
4835  assert( nodepair != NULL);
4836 
4837  source = nodepair->node1;
4838  target = nodepair->node2;
4839 
4840  assert(source >=0 && source < mcfnetwork->nnodes);
4841  assert(target >=0 && target < mcfnetwork->nnodes);
4842  assert(source <= target);
4843 
4844  hashval = (source << 16) + target; /*lint !e701*/
4845 
4846  return hashval;
4847 }
4848 
4849 /** creates a priority queue and fills it with the given nodepair entries
4850  *
4851  */
4852 static
4854  SCIP* scip, /**< SCIP data structure */
4855  SCIP_MCFNETWORK* mcfnetwork, /**< MCF network structure */
4856  NODEPAIRQUEUE** nodepairqueue /**< pointer to nodepair priority queue */
4857  )
4858 {
4859  /* For every nodepair that is used in the network (at least one arc exists having this nodepair as endnodes)
4860  * we calculate a weight:
4861  * The weight w_st of a nodepair (s,t) is the minimum of the weights of all s-t and t-s arcs
4862  * The weight w_a of an arc a is calculated as:
4863  * w_a : = s_a + pi_a
4864  * where s_a>=0 is the slack of the capacity constraint and pi_a<=0 its dual.
4865  * The weight of uncapacitated arcs (without capacity constraints) is infinite.
4866  */
4867 #ifdef BETTERWEIGHTFORDEMANDNODES
4868  int ncommodities;
4871  SCIP_Real maxweight;
4872  SCIP_Real minweight;
4873 #endif
4874 
4875 #ifdef TIEBREAKING
4876  int* colcommodity;
4877 #endif
4878 
4879  SCIP_HASHTABLE* hashtable;
4880  NODEPAIRENTRY* nodepairs;
4881 
4882  int hashtablesize;
4883  int a;
4884  int nnodepairs;
4885  int n;
4886 
4887  assert(mcfnetwork != NULL);
4888 
4889 #ifdef BETTERWEIGHTFORDEMANDNODES
4890  ncommodities = mcfnetwork->ncommodities;
4891  nodeflowrows = mcfnetwork->nodeflowrows;
4892  nodeflowscales = mcfnetwork->nodeflowscales;
4893 #endif
4894 
4895 #ifdef TIEBREAKING
4896  colcommodity = mcfnetwork->colcommodity;
4897 #endif
4898 
4899  assert(nodepairqueue != NULL);
4900 
4901  SCIP_CALL( SCIPallocBuffer(scip, nodepairqueue) );
4902 
4903  /* create a hash table for all used node pairs
4904  * hash table is only needed to have unique nodepairs (identify arcs using the same nodepair)
4905  */
4906  hashtablesize = mcfnetwork->narcs;
4907  hashtablesize = MAX(hashtablesize, HASHSIZE_NODEPAIRS);
4908  SCIP_CALL( SCIPhashtableCreate(&hashtable, SCIPblkmem(scip), hashtablesize,
4909  hashGetKeyNodepairs, hashKeyEqNodepairs, hashKeyValNodepairs, (void*) mcfnetwork) );
4910 
4911  /* nodepairs will contain all constructed nodepairs and is used to fill the priority queue */
4912  SCIP_CALL( SCIPallocBufferArray(scip, &(*nodepairqueue)->nodepairs, mcfnetwork->narcs) );
4913 
4914  /* initialize hash table of all used node pairs and fill nodepairs */
4915  nnodepairs = 0;
4916  for( a = 0; a < mcfnetwork->narcs; a++ )
4917  {
4918  NODEPAIRENTRY nodepair;
4919  NODEPAIRENTRY* nodepairptr;
4920  SCIP_ROW* capacityrow;
4921 
4922  capacityrow = mcfnetwork->arccapacityrows[a];
4923 
4924  SCIPdebugMsg(scip, "arc %i = (%i %i)\n", a, mcfnetwork->arcsources[a], mcfnetwork->arctargets[a]);
4925 
4926  /* construct fresh nodepair: smaller node gets node1 in nodeentry */
4927  if( mcfnetwork->arcsources[a] <= mcfnetwork->arctargets[a] )
4928  {
4929  nodepair.node1 = mcfnetwork->arcsources[a];
4930  nodepair.node2 = mcfnetwork->arctargets[a];
4931  }
4932  else
4933  {
4934  nodepair.node2 = mcfnetwork->arcsources[a];
4935  nodepair.node1 = mcfnetwork->arctargets[a];
4936  }
4937 
4938  assert(nodepair.node1 < mcfnetwork->nnodes);
4939  assert(nodepair.node2 < mcfnetwork->nnodes);
4940  if( nodepair.node1 == -1 || nodepair.node2 == -1 )
4941  continue;
4942 
4943  /* construct arc weight of a */
4944  if( capacityrow != NULL )
4945  {
4946  SCIP_Real maxval;
4947  SCIP_Real slack;
4948  SCIP_Real dualsol;
4949  SCIP_Real scale;
4950 #ifdef TIEBREAKING
4951  SCIP_Real totalflow;
4952  SCIP_Real totalcap;
4953  SCIP_COL** rowcols;
4954  int rowlen;
4955  int i;
4956  int c;
4957 #endif
4958 
4959  slack = SCIPgetRowFeasibility(scip, mcfnetwork->arccapacityrows[a]);
4960  slack = MAX(slack, 0.0); /* can only be negative due to numerics */
4961  dualsol = SCIProwGetDualsol(mcfnetwork->arccapacityrows[a]);
4962  maxval = SCIPgetRowMaxCoef(scip, mcfnetwork->arccapacityrows[a]);
4963  scale = ABS(mcfnetwork->arccapacityscales[a])/maxval; /* divide by maxval to normalize rows */
4964  assert(scale > 0.0);
4965 
4966 #ifdef TIEBREAKING
4967  /* get flow on arc for tie breaking */
4968  rowcols = SCIProwGetCols(capacityrow);
4969  rowlen = SCIProwGetNLPNonz(capacityrow);
4970  totalflow = 0.0;
4971  totalcap = 0.0;
4972  SCIPdebugMsg(scip, " row <%s>: \n", SCIProwGetName(capacityrow));
4973 
4974  for( i = 0; i < rowlen; i++ )
4975  {
4976  c = SCIPcolGetLPPos(rowcols[i]);
4977  assert(0 <= c && c < SCIPgetNLPCols(scip));
4978 
4979  SCIPdebugMsg(scip, " col <%s>: %g\n", SCIPvarGetName(SCIPcolGetVar(rowcols[i])), SCIPcolGetPrimsol(rowcols[i]) );
4980  /* sum up flow on arc a*/
4981  if(colcommodity[c] >= 0)
4982  {
4983  SCIPdebugMsg(scip, " flow col <%s>: %g\n", SCIPvarGetName(SCIPcolGetVar(rowcols[i])), REALABS(SCIPcolGetPrimsol(rowcols[i])) );
4984  totalflow += REALABS(SCIPcolGetPrimsol(rowcols[i]));
4985  }
4986  else
4987  {
4988  SCIPdebugMsg(scip, " cap col <%s>: %g\n", SCIPvarGetName(SCIPcolGetVar(rowcols[i])), REALABS(SCIPcolGetPrimsol(rowcols[i])) );
4989  totalcap += REALABS(SCIPcolGetPrimsol(rowcols[i]));
4990  }
4991  }
4992 
4993  SCIPdebugMsg(scip, "cap arc -- slack:%g -- dual:%g -- flow:%g -- cap:%g \n", scale * slack, dualsol/scale, totalflow * scale, totalcap * scale);
4994 #else
4995  SCIPdebugMsg(scip, "cap arc -- slack:%g -- dual:%g1\n", scale * slack, dualsol/scale);
4996 #endif
4997 
4998  /* put the arc weight into a fresh nodepair */
4999  nodepair.weight = scale * slack - ABS(dualsol)/scale;
5000 #ifdef USEFLOWFORTIEBREAKING
5001  if( !SCIPisPositive(scip, nodepair.weight) )
5002  {
5003  nodepair.weight += totalflow * scale;
5004  nodepair.weight = MIN( nodepair.weight, -0.0001);
5005  }
5006 #endif
5007 #ifdef USECAPACITYFORTIEBREAKING
5008  if( !SCIPisPositive(scip, nodepair.weight) )
5009  {
5010  nodepair.weight += totalcap * scale;
5011  nodepair.weight = MIN( nodepair.weight, -0.0001);
5012  }
5013 #endif
5014  }
5015  else
5016  {
5017  /* uncapacitated arc has infinite slack */
5018  SCIPdebugMsg(scip, "uncap arc ... slack infinite\n");
5019  nodepair.weight = SCIPinfinity(scip);
5020  }
5021 
5022  /* check if nodepair already exists in hash-table */
5023  nodepairptr = (NODEPAIRENTRY*)(SCIPhashtableRetrieve(hashtable, (void*) (&nodepair) ));
5024 
5025  /* if nodepair already exists update its weight */
5026  if( nodepairptr != NULL )
5027  {
5028  /* adapt weight */
5029  SCIPdebugMsg(scip, "nodepair known [%d,%d] -- old weight:%g -- new weight:%g\n", nodepair.node1,nodepair.node2,nodepairptr->weight,
5030  MIN(nodepair.weight, nodepairptr->weight));
5031  nodepairptr->weight = MIN(nodepair.weight, nodepairptr->weight);
5032  }
5033  else
5034  {
5035  /* no such nodepair in current hash table: insert into array and hashtable */
5036  nodepairs = (*nodepairqueue)->nodepairs;
5037 
5038  assert(nnodepairs < mcfnetwork->narcs);
5039  nodepairs[nnodepairs] = nodepair;
5040  SCIP_CALL( SCIPhashtableInsert(hashtable, (void*) (&nodepairs[nnodepairs]) ) );
5041 
5042  SCIPdebugMsg(scip, "new nodepair [%d,%d]-- weight:%g\n", nodepair.node1, nodepair.node2, nodepair.weight);
5043 
5044  nnodepairs++;
5045  }
5046  }
5047 
5048  /* free hash table */
5049  SCIPhashtableFree(&hashtable);
5050 
5051  /* we still have to fill the priority queue */
5052 
5053 #ifdef BETTERWEIGHTFORDEMANDNODES
5054  /* initial weights have been calculated
5055  * we modify them now depending on the demand emanating at the endnodes of nodepairs
5056  */
5057 
5058  /* calculate max and min weight */
5059  maxweight = +1; /* we want maxweight to be positive */
5060  minweight = -1; /* we want minweight to be negative */
5061  nodepairs = (*nodepairqueue)->nodepairs;
5062  for( n = 0; n < nnodepairs; n++ )
5063  {
5064  /* maxweight should not be infinity (uncap arcs have infinity weight)*/
5065  if(!SCIPisInfinity(scip,nodepairs[n].weight))
5066  maxweight = MAX(maxweight, nodepairs[n].weight);
5067 
5068  minweight = MIN(minweight, nodepairs[n].weight);
5069  }
5070 
5071  SCIPdebugMsg(scip, "min/max weight:%g / %g\n", minweight, maxweight);
5072 #endif
5073 
5074  /* initialize priority queue */
5075  SCIP_CALL( SCIPpqueueCreate(&(*nodepairqueue)->pqueue, nnodepairs, 2.0, compNodepairs) );
5076 
5077  /* fill priority queue using array nodepairs */
5078  for( n = 0; n < nnodepairs; n++ )
5079  {
5080  int node1 = nodepairs[n].node1;
5081  int node2 = nodepairs[n].node2;
5082 
5083 #ifdef BETTERWEIGHTFORDEMANDNODES
5084  SCIP_Real rhs = 0;
5085  SCIP_Bool hasdemand1 = FALSE;
5086  SCIP_Bool hasdemand2 = FALSE;
5087 
5088  int k; /* commodity */
5089 
5090  SCIPdebugMsg(scip, "nodepair [%d,%d] weight %g\n", node1,node2,nodepairs[n].weight);
5091  /* check both nodes for their demand value in all commodities
5092  * the demand value can be read from the rhs
5093  * of the flowrows
5094  */
5095  /* node1 */
5096  for( k = 0; k < ncommodities; k++ )
5097  {
5098  if( nodeflowrows[node1][k] == NULL )
5099  continue;
5100 
5101  if( nodeflowscales[node1][k] > 0.0 )
5102  rhs = SCIProwGetRhs(nodeflowrows[node1][k]) - SCIProwGetConstant(nodeflowrows[node1][k]);
5103  else
5104  rhs = SCIProwGetLhs(nodeflowrows[node1][k]) - SCIProwGetConstant(nodeflowrows[node1][k]);
5105 
5106  assert( !SCIPisInfinity(scip,ABS(rhs)) );
5107 
5108  if( ! SCIPisZero(scip, rhs) )
5109  {
5110  hasdemand1 = TRUE;
5111  break;
5112  }
5113  }
5114 
5115  /* node2 */
5116  for( k = 0; k < ncommodities; k++ )
5117  {
5118  if( nodeflowrows[node2][k] == NULL )
5119  continue;
5120 
5121  if( nodeflowscales[node2][k] > 0.0 )
5122  rhs = SCIProwGetRhs(nodeflowrows[node2][k]) - SCIProwGetConstant(nodeflowrows[node2][k]);
5123  else
5124  rhs = SCIProwGetLhs(nodeflowrows[node2][k]) - SCIProwGetConstant(nodeflowrows[node2][k]);
5125 
5126  assert(! SCIPisInfinity(scip, ABS(rhs)));
5127 
5128  if( ! SCIPisZero(scip, rhs) )
5129  {
5130  hasdemand2 = TRUE;
5131  break;
5132  }
5133  }
5134 
5135  /* if one of the nodes has no demand increase the score
5136  * (slack arcs are still shrunk first)
5137  *
5138  */
5139  if( SCIPisPositive(scip, nodepairs[n].weight))
5140  {
5141  assert(SCIPisPositive(scip, maxweight));
5142 
5143  if( !hasdemand1 || !hasdemand2 )
5144  nodepairs[n].weight += maxweight;
5145  }
5146  else
5147  {
5148  assert( SCIPisNegative(scip, minweight));
5149 
5150  if( hasdemand1 && hasdemand2)
5151  nodepairs[n].weight += minweight;
5152  }
5153 #endif
5154  SCIPdebugMsg(scip, "nodepair [%d,%d] weight %g\n", node1,node2,nodepairs[n].weight);
5155 
5156  /* fill priority queue */
5157  SCIP_CALL( SCIPpqueueInsert((*nodepairqueue)->pqueue, (void*)&(*nodepairqueue)->nodepairs[n]) );
5158  }
5159 
5160  return SCIP_OKAY;
5161 }
5162 
5163 
5164 /** frees memory of a nodepair queue */
5165 static
5167  SCIP* scip, /**< SCIP data structure */
5168  NODEPAIRQUEUE** nodepairqueue /**< pointer to nodepair priority queue */
5169  )
5170 {
5171  assert(nodepairqueue != NULL);
5172  assert(*nodepairqueue != NULL);
5173 
5174  SCIPpqueueFree(&(*nodepairqueue)->pqueue);
5175  SCIPfreeBufferArray(scip, &(*nodepairqueue)->nodepairs);
5176  SCIPfreeBuffer(scip, nodepairqueue);
5177 }
5178 
5179 
5180 /** returns whether there are any nodepairs left on the queue */
5181 static
5183  NODEPAIRQUEUE* nodepairqueue /**< nodepair priority queue */
5184  )
5185 {
5186  assert(nodepairqueue != NULL);
5187 
5188  return (SCIPpqueueFirst(nodepairqueue->pqueue) == NULL);
5189 }
5190 
5191 
5192 /** removes the top element from the nodepair priority queue, returns nodepair entry or NULL */
5193 static
5195  NODEPAIRQUEUE* nodepairqueue /**< nodepair priority queue */
5196  )
5197 {
5198  assert(nodepairqueue != NULL);
5199 
5200  return (NODEPAIRENTRY*)SCIPpqueueRemove(nodepairqueue->pqueue);
5201 }
5202 
5203 /*
5204  * Node partition methods
5205  * used for creating partitions of nodes
5206  * use union find algorithm
5207 */
5208 
5209 /** returns the representative node in the cluster of the given node */
5210 static
5212  NODEPARTITION* nodepartition, /**< node partition data structure */
5213  int v /**< node id to get representative for */
5214  )
5215 {
5216  return unionfindGetRepresentative(nodepartition->representatives, v);
5217 }
5218 
5219 /** joins two clusters given by their representative nodes */
5220 static
5222  NODEPARTITION* nodepartition, /**< node partition data structure */
5223  int rep1, /**< representative of first cluster */
5224  int rep2 /**< representative of second cluster */
5225  )
5226 {
5227  unionfindJoinSets (nodepartition->representatives, rep1, rep2);
5228 }
5229 
5230 /** partitions nodes into a small number of clusters */
5231 static
5233  SCIP* scip, /**< SCIP data structure */
5234  SCIP_MCFNETWORK* mcfnetwork, /**< MCF network structure */
5235  NODEPARTITION** nodepartition, /**< pointer to node partition data structure */
5236  int nclusters /**< number of clusters to generate */
5237  )
5238 {
5239  NODEPAIRQUEUE* nodepairqueue;
5240  int* clustersize;
5241  int nclustersleft;
5242  int v;
5243  int c;
5244  int pos;
5245 
5246  assert(mcfnetwork != NULL);
5247  assert(nodepartition != NULL);
5248  assert(mcfnetwork->nnodes >= 1);
5249 
5250  /* allocate and initialize memory */
5251  SCIP_CALL( SCIPallocBuffer(scip, nodepartition) );
5252  SCIP_CALL( SCIPallocBufferArray(scip, &(*nodepartition)->representatives, mcfnetwork->nnodes) );
5253  SCIP_CALL( SCIPallocBufferArray(scip, &(*nodepartition)->nodeclusters, mcfnetwork->nnodes) );
5254  SCIP_CALL( SCIPallocBufferArray(scip, &(*nodepartition)->clusternodes, mcfnetwork->nnodes) );
5255  SCIP_CALL( SCIPallocBufferArray(scip, &(*nodepartition)->clusterbegin, nclusters+1) );
5256  (*nodepartition)->nclusters = 0;
5257 
5258  /* we start with each node being in its own cluster */
5259  unionfindInitSets((*nodepartition)->representatives, mcfnetwork->nnodes);
5260 
5261  /* create priority queue for nodepairs */
5262  SCIP_CALL( nodepairqueueCreate(scip, mcfnetwork, &nodepairqueue) );
5263 
5264  /* loop over nodepairs in order of their weights */
5265  nclustersleft = mcfnetwork->nnodes;
5266  while( !nodepairqueueIsEmpty(nodepairqueue) && nclustersleft > nclusters )
5267  {
5268  NODEPAIRENTRY* nodepair;
5269  int node1;
5270  int node2;
5271  SCIP_Real weight;
5272  int node1rep;
5273  int node2rep;
5274 
5275  /* get the next nodepair */
5276  nodepair = nodepairqueueRemove(nodepairqueue);
5277  assert(nodepair != NULL);
5278  node1 = nodepair->node1;
5279  node2 = nodepair->node2;
5280  SCIPdebug( weight = nodepair->weight; )
5281 
5282  assert(node1 >= 0 && node1 < mcfnetwork->nnodes);
5283  assert(node2 >= 0 && node2 < mcfnetwork->nnodes);
5284 
5285  /* identify the representatives of the two nodes */
5286  node1rep = nodepartitionGetRepresentative(*nodepartition, node1);
5287  node2rep = nodepartitionGetRepresentative(*nodepartition, node2);
5288  assert(0 <= node1rep && node1rep < mcfnetwork->nnodes);
5289  assert(0 <= node2rep && node2rep < mcfnetwork->nnodes);
5290 
5291  /* there is nothing to do if the two nodes are already in the same cluster */
5292  if( node1rep == node2rep )
5293  continue;
5294 
5295  /* shrink nodepair by joining the two clusters */
5296  SCIPdebugMsg(scip, "shrinking nodepair (%d,%d) with weight %g: join representatives %d and %d\n",
5297  node1, node2, weight, node1rep, node2rep);
5298  nodepartitionJoin(*nodepartition, node1rep, node2rep);
5299  nclustersleft--;
5300  }
5301 
5302  /* node 0 must be a representative due to our join procedure */
5303  assert((*nodepartition)->representatives[0] == 0);
5304 
5305  /* if there have been too few arcs to shrink the graph to the required number of clusters, join clusters with first cluster
5306  * to create a larger disconnected cluster
5307  */
5308  if( nclustersleft > nclusters )
5309  {
5310  for( v = 1; v < mcfnetwork->nnodes && nclustersleft > nclusters; v++ )
5311  {
5312  int rep;
5313 
5314  rep = nodepartitionGetRepresentative(*nodepartition, v);
5315  if( rep != 0 )
5316  {
5317  nodepartitionJoin(*nodepartition, 0, rep);
5318  nclustersleft--;
5319  }
5320  }
5321  }
5322  assert(nclustersleft <= nclusters);
5323 
5324  /* extract the clusters */
5325  SCIP_CALL( SCIPallocBufferArray(scip, &clustersize, nclusters) );
5326  BMSclearMemoryArray(clustersize, nclusters);
5327  for( v = 0; v < mcfnetwork->nnodes; v++ )
5328  {
5329  int rep;
5330 
5331  /* get cluster of node */
5332  rep = nodepartitionGetRepresentative(*nodepartition, v);
5333  assert(rep <= v); /* due to our joining procedure */
5334  if( rep == v )
5335  {
5336  /* node is its own representative: this is a new cluster */
5337  c = (*nodepartition)->nclusters;
5338  (*nodepartition)->nclusters++;
5339  }
5340  else
5341  c = (*nodepartition)->nodeclusters[rep];
5342  assert(0 <= c && c < nclusters);
5343 
5344  /* assign node to cluster */
5345  (*nodepartition)->nodeclusters[v] = c;
5346  clustersize[c]++;
5347  }
5348 
5349  /* fill the clusterbegin array */
5350  pos = 0;
5351  for( c = 0; c < (*nodepartition)->nclusters; c++ )
5352  {
5353  (*nodepartition)->clusterbegin[c] = pos;
5354  pos += clustersize[c];
5355  }
5356  assert(pos == mcfnetwork->nnodes);
5357  (*nodepartition)->clusterbegin[(*nodepartition)->nclusters] = mcfnetwork->nnodes;
5358 
5359  /* fill the clusternodes array */
5360  BMSclearMemoryArray(clustersize, (*nodepartition)->nclusters);
5361  for( v = 0; v < mcfnetwork->nnodes; v++ )
5362  {
5363  c = (*nodepartition)->nodeclusters[v];
5364  assert(0 <= c && c < (*nodepartition)->nclusters);
5365  pos = (*nodepartition)->clusterbegin[c] + clustersize[c];
5366  assert(pos < (*nodepartition)->clusterbegin[c+1]);
5367  (*nodepartition)->clusternodes[pos] = v;
5368  clustersize[c]++;
5369  }
5370 
5371  /* free temporary memory */
5372  SCIPfreeBufferArray(scip, &clustersize);
5373 
5374  /* free nodepair queue */
5375  nodepairqueueFree(scip, &nodepairqueue);
5376 
5377  return SCIP_OKAY;
5378 }
5379 
5380 /** frees node partition data */
5381 static
5383  SCIP* scip, /**< SCIP data structure */
5384  NODEPARTITION** nodepartition /**< pointer to node partition data structure */
5385  )
5386 {
5387  assert(nodepartition != NULL);
5388  assert(*nodepartition != NULL);
5389 
5390  SCIPfreeBufferArray(scip, &(*nodepartition)->clusterbegin);
5391  SCIPfreeBufferArray(scip, &(*nodepartition)->clusternodes);
5392  SCIPfreeBufferArray(scip, &(*nodepartition)->nodeclusters);
5393  SCIPfreeBufferArray(scip, &(*nodepartition)->representatives);
5394  SCIPfreeBuffer(scip, nodepartition);
5395 }
5396 
5397 /** returns whether given node v is in a cluster that belongs to the partition S */
5398 static
5400  NODEPARTITION* nodepartition, /**< node partition data structure, or NULL */
5401  unsigned int partition, /**< partition of nodes, or node number in single-node partition */
5402  SCIP_Bool inverted, /**< should partition be inverted? */
5403  int v /**< node to check */
5404  )
5405 {
5406  /* if the node does not exist, it is not in the partition
5407  * (and also not in the inverted partition)
5408  */
5409  if( v < 0 )
5410  return FALSE;
5411 
5412  if( nodepartition == NULL )
5413  return ((v == (int)partition) == !inverted);
5414  else
5415  {
5416  int cluster;
5417  unsigned int clusterbit;
5418 
5419  cluster = nodepartition->nodeclusters[v];
5420  assert(0 <= cluster && cluster < nodepartition->nclusters);
5421  clusterbit = (1 << cluster); /*lint !e701*/
5422 
5423  return (((partition & clusterbit) != 0) == !inverted);
5424  }
5425 }
5426 
5427 /** returns whether each of the sets S and T defined by the nodepartition is connected */
5428 static
5430  SCIP* scip, /**< SCIP data structure */
5431  SCIP_MCFNETWORK* mcfnetwork, /**< MCF network structure */
5432  NODEPARTITION* nodepartition, /**< node partition data structure */
5433  unsigned int partition /**< partition of nodes, or node number in single-node partition */
5434  )
5435 {
5436  const int* nodeclusters = nodepartition->nodeclusters;
5437  const int* arcsources = mcfnetwork->arcsources;
5438  const int* arctargets = mcfnetwork->arctargets;
5439  int narcs = mcfnetwork->narcs;
5440  int nclusters;
5441 
5442  int ncomponents;
5443  int a;
5444  int* rep;
5445 
5446  assert(nodepartition->nodeclusters != NULL);
5447  nclusters = nodepartition->nclusters;
5448 
5449  if( SCIPallocBufferArray(scip, &rep, nclusters) != SCIP_OKAY )
5450  return 0;
5451 
5452  /* start with each cluster being isolated */
5453  unionfindInitSets(rep, nclusters);
5454  ncomponents = nclusters;
5455  assert(ncomponents >= 2);
5456 
5457  /* for each arc within S or within T join the connected clusters */
5458  for( a = 0; a < narcs; a++ )
5459  {
5460  int s = arcsources[a];
5461  int t = arctargets[a];
5462 
5463  /* ignore arcs that connect the pseudo node -1 */
5464  if( s == -1 || t == -1 )
5465  continue;
5466 
5467  /* check if arc is within one of the components */
5468  if( nodeInPartition(nodepartition, partition, FALSE, s) == nodeInPartition(nodepartition, partition, FALSE, t) )
5469  {
5470  int cs;
5471  int ct;
5472  int repcs;
5473  int repct;
5474 
5475  assert(0 <= s && s < mcfnetwork->nnodes);
5476  assert(0 <= t && t < mcfnetwork->nnodes);
5477 
5478  /* get clusters of the two nodes */
5479  cs = nodeclusters[s];
5480  ct = nodeclusters[t];
5481  assert(0 <= cs && cs < nclusters);
5482  assert(0 <= ct && ct < nclusters);
5483 
5484  /* nothing to do if we are already in the same cluster */
5485  if( cs == ct )
5486  continue;
5487 
5488  /* get representatives of clusters in the union structure */
5489  repcs = unionfindGetRepresentative (rep, cs);
5490  repct = unionfindGetRepresentative (rep, ct);
5491  if( repcs == repct )
5492  continue;
5493 
5494  /* the arc connects two previously unconnected components of S or T */
5495 
5496  /* check if we already reached two distinct components */
5497  ncomponents--;
5498  if( ncomponents <= 2 )
5499  break;
5500 
5501  /* join the two cluster sets and continue */
5502  unionfindJoinSets (rep, repcs, repct);
5503  }
5504  }
5505 
5506  /* each of the two shores S and T are connected if we were able to shrink the graph
5507  into two components */
5508  assert(ncomponents >= 2);
5509  SCIPfreeBufferArray(scip, &rep);
5510  return (ncomponents == 2);
5511 }
5512 
5513 #ifdef SCIP_DEBUG
5514 static
5515 void nodepartitionPrint(
5516  NODEPARTITION* nodepartition /**< node partition data structure */
5517  )
5518 {
5519  int c;
5520 
5521  for( c = 0; c < nodepartition->nclusters; c++ )
5522  {
5523  int i;
5524 
5525  MCFdebugMessage("cluster %d:", c);
5526  for( i = nodepartition->clusterbegin[c]; i < nodepartition->clusterbegin[c+1]; i++ )
5527  MCFdebugMessage(" %d", nodepartition->clusternodes[i]);
5528  MCFdebugMessage("\n");
5529  }
5530 }
5531 #endif
5532 
5533 #ifdef OUTPUTGRAPH
5534 /** generates a GML file to visualize the network graph and LP solution */
5535 static
5536 SCIP_RETCODE outputGraph(
5537  SCIP* scip, /**< SCIP data structure */
5538  SCIP_MCFNETWORK* mcfnetwork, /**< MCF network structure */
5539  NODEPARTITION* nodepartition, /**< node partition data structure, or NULL */
5540  SCIP_Bool inverted, /**< should partition be inverted? */
5541  unsigned int partition /**< partition of nodes, or node number */
5542  )
5543 {
5544  FILE* file;
5545  char filename[SCIP_MAXSTRLEN];
5546  int v;
5547  int a;
5548 
5549  /* open file */
5550  if( nodepartition == NULL )
5551  (void) SCIPsnprintf(filename, SCIP_MAXSTRLEN, "mcf-node-%d.gml", partition);
5552  else
5553  (void) SCIPsnprintf(filename, SCIP_MAXSTRLEN, "mcf-part-%d.gml", partition);
5554  SCIPinfoMessage(scip, NULL, "creating GML output file <%s>...\n", filename);
5555  file = fopen(filename, "w");
5556  if( file == NULL )
5557  {
5558  SCIPerrorMessage("cannot create GML output file <%s>\n", filename);
5559  return SCIP_FILECREATEERROR;
5560  }
5561 
5562  /* print GML header */
5563  fprintf(file, "graph\n");
5564  fprintf(file, "[\n");
5565  fprintf(file, " hierarchic 1\n");
5566  fprintf(file, " label \"\"\n");
5567  fprintf(file, " directed 1\n");
5568 
5569  /* nodes */
5570  for( v = 0; v < mcfnetwork->nnodes; v++ )
5571  {
5572  char label[SCIP_MAXSTRLEN];
5573  SCIP_Bool inpartition;
5574 
5575  if( mcfnetwork->nodeflowrows[v][0] != NULL )
5576  (void) SCIPsnprintf(label, SCIP_MAXSTRLEN, "%s", SCIProwGetName(mcfnetwork->nodeflowrows[v][0]));
5577  else
5578  (void) SCIPsnprintf(label, SCIP_MAXSTRLEN, "%d", v);
5579  inpartition = nodeInPartition(nodepartition, partition, inverted, v);
5580 
5581  fprintf(file, " node\n");
5582  fprintf(file, " [\n");
5583  fprintf(file, " id %d\n", v);
5584  fprintf(file, " label \"%s\"\n", label);
5585  fprintf(file, " graphics\n");
5586  fprintf(file, " [\n");
5587  fprintf(file, " w 30.0\n");
5588  fprintf(file, " h 30.0\n");
5589  fprintf(file, " type \"ellipse\"\n");
5590  if( inpartition )
5591  fprintf(file, " fill \"#FF0000\"\n");
5592  else
5593  fprintf(file, " fill \"#00FF00\"\n");
5594  fprintf(file, " outline \"#000000\"\n");
5595  fprintf(file, " ]\n");
5596  fprintf(file, " LabelGraphics\n");
5597  fprintf(file, " [\n");
5598  fprintf(file, " text \"%s\"\n", label);
5599  fprintf(file, " fontSize 13\n");
5600  fprintf(file, " fontName \"Dialog\"\n");
5601  fprintf(file, " anchor \"c\"\n");
5602  fprintf(file, " ]\n");
5603  if( inpartition )
5604  fprintf(file, " gid %d\n", mcfnetwork->nnodes+1);
5605  else
5606  fprintf(file, " gid %d\n", mcfnetwork->nnodes+2);
5607  fprintf(file, " ]\n");
5608  }
5609 
5610  /* dummy node for missing arc sources or arc targets */
5611  fprintf(file, " node\n");
5612  fprintf(file, " [\n");
5613  fprintf(file, " id %d\n", mcfnetwork->nnodes);
5614  fprintf(file, " label \"?\"\n");
5615  fprintf(file, " graphics\n");
5616  fprintf(file, " [\n");
5617  fprintf(file, " w 30.0\n");
5618  fprintf(file, " h 30.0\n");
5619  fprintf(file, " type \"ellipse\"\n");
5620  fprintf(file, " fill \"#FFFFFF\"\n");
5621  fprintf(file, " outline \"#000000\"\n");
5622  fprintf(file, " ]\n");
5623  fprintf(file, " LabelGraphics\n");
5624  fprintf(file, " [\n");
5625  fprintf(file, " text \"?\"\n");
5626  fprintf(file, " fontSize 13\n");
5627  fprintf(file, " fontName \"Dialog\"\n");
5628  fprintf(file, " anchor \"c\"\n");
5629  fprintf(file, " ]\n");
5630  fprintf(file, " ]\n");
5631 
5632  /* group node for partition S */
5633  fprintf(file, " node\n");
5634  fprintf(file, " [\n");
5635  fprintf(file, " id %d\n", mcfnetwork->nnodes+1);
5636  fprintf(file, " label \"Partition S\"\n");
5637  fprintf(file, " graphics\n");
5638  fprintf(file, " [\n");
5639  fprintf(file, " type \"roundrectangle\"\n");
5640  fprintf(file, " fill \"#CAECFF84\"\n");
5641  fprintf(file, " outline \"#666699\"\n");
5642  fprintf(file, " outlineStyle \"dotted\"\n");
5643  fprintf(file, " topBorderInset 0\n");
5644  fprintf(file, " bottomBorderInset 0\n");
5645  fprintf(file, " leftBorderInset 0\n");
5646  fprintf(file, " rightBorderInset 0\n");
5647  fprintf(file, " ]\n");
5648  fprintf(file, " LabelGraphics\n");
5649  fprintf(file, " [\n");
5650  fprintf(file, " text \"Partition S\"\n");
5651  fprintf(file, " fill \"#99CCFF\"\n");
5652  fprintf(file, " fontSize 15\n");
5653  fprintf(file, " fontName \"Dialog\"\n");
5654  fprintf(file, " alignment \"right\"\n");
5655  fprintf(file, " autoSizePolicy \"node_width\"\n");
5656  fprintf(file, " anchor \"t\"\n");
5657  fprintf(file, " borderDistance 0.0\n");
5658  fprintf(file, " ]\n");
5659  fprintf(file, " isGroup 1\n");
5660  fprintf(file, " ]\n");
5661 
5662  /* group node for partition T */
5663  fprintf(file, " node\n");
5664  fprintf(file, " [\n");
5665  fprintf(file, " id %d\n", mcfnetwork->nnodes+2);
5666  fprintf(file, " label \"Partition T\"\n");
5667  fprintf(file, " graphics\n");
5668  fprintf(file, " [\n");
5669  fprintf(file, " type \"roundrectangle\"\n");
5670  fprintf(file, " fill \"#CAECFF84\"\n");
5671  fprintf(file, " outline \"#666699\"\n");
5672  fprintf(file, " outlineStyle \"dotted\"\n");
5673  fprintf(file, " topBorderInset 0\n");
5674  fprintf(file, " bottomBorderInset 0\n");
5675  fprintf(file, " leftBorderInset 0\n");
5676  fprintf(file, " rightBorderInset 0\n");
5677  fprintf(file, " ]\n");
5678  fprintf(file, " LabelGraphics\n");
5679  fprintf(file, " [\n");
5680  fprintf(file, " text \"Partition T\"\n");
5681  fprintf(file, " fill \"#99CCFF\"\n");
5682  fprintf(file, " fontSize 15\n");
5683  fprintf(file, " fontName \"Dialog\"\n");
5684  fprintf(file, " alignment \"right\"\n");
5685  fprintf(file, " autoSizePolicy \"node_width\"\n");
5686  fprintf(file, " anchor \"t\"\n");
5687  fprintf(file, " borderDistance 0.0\n");
5688  fprintf(file, " ]\n");
5689  fprintf(file, " isGroup 1\n");
5690  fprintf(file, " ]\n");
5691 
5692  /* arcs */
5693  for( a = 0; a < mcfnetwork->narcs; a++ )
5694  {
5695  SCIP_ROW* row;
5696  SCIP_Real slack;
5697  SCIP_Bool hasfractional;
5698  char label[SCIP_MAXSTRLEN];
5699 
5700  if( mcfnetwork->arccapacityrows[a] != NULL )
5701  (void) SCIPsnprintf(label, SCIP_MAXSTRLEN, "%s", SCIProwGetName(mcfnetwork->arccapacityrows[a]));
5702  else
5703  (void) SCIPsnprintf(label, SCIP_MAXSTRLEN, "%d", a);
5704 
5705  hasfractional = FALSE;
5706  row = mcfnetwork->arccapacityrows[a];
5707  if( row != NULL )
5708  {
5709  SCIP_COL** rowcols;
5710  int rowlen;
5711  int i;
5712 
5713  slack = ABS(mcfnetwork->arccapacityscales[a]) * SCIPgetRowLPFeasibility(scip, row);
5714  rowcols = SCIProwGetCols(row);
5715  rowlen = SCIProwGetNLPNonz(row);
5716  for( i = 0; i < rowlen; i++ )
5717  {
5718  SCIP_VAR* var;
5719 
5720  var = SCIPcolGetVar(rowcols[i]);
5721  if( SCIPvarIsIntegral(var) && !SCIPisFeasIntegral(scip, SCIPvarGetLPSol(var)) )
5722  {
5723  hasfractional = TRUE;
5724  break;
5725  }
5726  }
5727  }
5728  else
5729  slack = SCIPinfinity(scip);
5730 
5731  fprintf(file, " edge\n");
5732  fprintf(file, " [\n");
5733  fprintf(file, " source %d\n", mcfnetwork->arcsources[a] >= 0 ? mcfnetwork->arcsources[a] : mcfnetwork->nnodes);
5734  fprintf(file, " target %d\n", mcfnetwork->arctargets[a] >= 0 ? mcfnetwork->arctargets[a] : mcfnetwork->nnodes);
5735  fprintf(file, " label \"%s\"\n", label);
5736  fprintf(file, " graphics\n");
5737  fprintf(file, " [\n");
5738  if( SCIPisFeasPositive(scip, slack) )
5739  fprintf(file, " fill \"#000000\"\n");
5740  else
5741  fprintf(file, " fill \"#FF0000\"\n");
5742  if( hasfractional )
5743  fprintf(file, " style \"dashed\"\n");
5744  fprintf(file, " width 1\n");
5745  fprintf(file, " targetArrow \"standard\"\n");
5746  fprintf(file, " ]\n");
5747  fprintf(file, " LabelGraphics\n");
5748  fprintf(file, " [\n");
5749  fprintf(file, " text \"%s\"\n", label);
5750  fprintf(file, " ]\n");
5751  fprintf(file, " ]\n");
5752  }
5753 
5754  /* print GML footer */
5755  fprintf(file, "]\n");
5756 
5757  /* close file */
5758  fclose(file);
5759 
5760  return SCIP_OKAY;
5761 }
5762 #endif
5763 
5764 /** adds given cut to LP if violated */
5765 static
5767  SCIP* scip, /**< SCIP data structure */
5768  SCIP_SEPA* sepa, /**< separator */
5769  SCIP_SEPADATA* sepadata, /**< separator data */
5770  SCIP_SOL* sol, /**< the solution that should be separated, or NULL for LP solution */
5771  SCIP_Real* cutcoefs, /**< coefficients of active variables in cut */
5772  SCIP_Real cutrhs, /**< right hand side of cut */
5773  int* cutinds, /**< problem indices of variables with non-zero coefficient */
5774  int cutnnz, /**< number of non-zeros in cut */
5775  SCIP_Bool cutislocal, /**< is the cut only locally valid? */
5776  int cutrank, /**< rank of the cut */
5777  int* ncuts, /**< pointer to count the number of added cuts */
5778  SCIP_Bool* cutoff /**< pointer to store whether a cutoff was found */
5779  )
5780 {
5781  SCIP_ROW* cut;
5782  char cutname[SCIP_MAXSTRLEN];
5783  SCIP_VAR** vars;
5784  int nvars;
5785  int i;
5786 
5787  /* variables for knapsack cover separation */
5788  SCIP_VAR** cutvars;
5789 
5790  assert(scip != NULL);
5791  assert(sepadata != NULL);
5792  assert(cutcoefs != NULL);
5793  assert(ncuts != NULL);
5794  assert(cutoff != NULL);
5795 
5796  /* get active problem variables */
5797  SCIP_CALL( SCIPgetVarsData(scip, &vars, &nvars, NULL, NULL, NULL, NULL) );
5798  assert(nvars == 0 || vars != NULL);
5799 
5800  *cutoff = FALSE;
5801 
5802  SCIP_CALL( SCIPallocBufferArray(scip, &cutvars, cutnnz) );
5803 
5804  for( i = 0; i < cutnnz; ++i )
5805  {
5806  cutvars[i] = vars[cutinds[i]];
5807  }
5808 
5809  /* create the cut */
5810  (void) SCIPsnprintf(cutname, SCIP_MAXSTRLEN, "mcf%d_%d", SCIPgetNLPs(scip), *ncuts);
5811  SCIP_CALL( SCIPcreateEmptyRowSepa(scip, &cut, sepa, cutname, -SCIPinfinity(scip), cutrhs, cutislocal, FALSE,
5812  sepadata->dynamiccuts) );
5813 
5814  SCIP_CALL( SCIPaddVarsToRow(scip, cut, cutnnz, cutvars, cutcoefs) );
5815 
5816  /* set cut rank */
5817  SCIProwChgRank(cut, cutrank);
5818 
5819  /* check efficacy */
5820  assert(SCIPisCutEfficacious(scip, sol, cut));
5821 
5822  SCIPdebugMsg(scip, " -> found MCF cut <%s>: rhs=%f, act=%f eff=%f rank=%d\n",
5823  cutname, cutrhs, SCIPgetRowSolActivity(scip, cut, sol), SCIPgetCutEfficacy(scip, sol, cut), SCIProwGetRank(cut));
5824  /*SCIPdebug( SCIP_CALL(SCIPprintRow(scip, cut, NULL)) );*/
5825 
5826  if( !cutislocal )
5827  {
5828  SCIP_CALL( SCIPaddPoolCut(scip, cut) );
5829  }
5830  else
5831  {
5832  SCIP_CALL( SCIPaddRow(scip, cut, FALSE, cutoff) );
5833  }
5834  (*ncuts)++;
5835 
5836  /* release the row */
5837  SCIP_CALL( SCIPreleaseRow(scip, &cut) );
5838 
5839  if( !(*cutoff) && sepadata->separateknapsack)
5840  {
5841  /* relax cut to knapsack row and separate lifted cover cuts */
5842  SCIP_CALL( SCIPseparateRelaxedKnapsack(scip, NULL, sepa, cutnnz, cutvars, cutcoefs, +1.0, cutrhs, sol, cutoff, ncuts) );
5843  }
5844  SCIPfreeBufferArray(scip, &cutvars);
5845 
5846  return SCIP_OKAY;
5847 }
5848 
5849 /** enumerates cuts between subsets of the clusters
5850  * generates single-node cuts if nodepartition == NULL, otherwise generates cluster cuts
5851  */
5852 static
5854  SCIP* scip, /**< SCIP data structure */
5855  SCIP_SEPA* sepa, /**< separator */
5856  SCIP_SEPADATA* sepadata, /**< separator data */
5857  SCIP_SOL* sol, /**< the solution that should be separated, or NULL for LP solution */
5858  SCIP_Bool allowlocal, /**< should local cuts be allowed */
5859  SCIP_MCFNETWORK* mcfnetwork, /**< MCF network structure */
5860  NODEPARTITION* nodepartition, /**< node partition data structure, or NULL */
5861  int* ncuts, /**< pointer to count the number of added cuts */
5862  SCIP_Bool* cutoff /**< pointer to store whether a cutoff was found */
5863  )
5864 {
5865  SCIP_ROW*** nodeflowrows = mcfnetwork->nodeflowrows;
5866  SCIP_Real** nodeflowscales = mcfnetwork->nodeflowscales;
5868  SCIP_ROW** arccapacityrows = mcfnetwork->arccapacityrows;
5870  int* colcommodity = mcfnetwork->colcommodity;
5871  int* arcsources = mcfnetwork->arcsources;
5872  int* arctargets = mcfnetwork->arctargets;
5873  int nnodes = mcfnetwork->nnodes;
5874  int narcs = mcfnetwork->narcs;
5875  int ncommodities = mcfnetwork->ncommodities;
5876  SCIP_MCFMODELTYPE modeltype = mcfnetwork->modeltype;
5877  MCFEFFORTLEVEL effortlevel = sepadata->effortlevel;
5878  int maxsepacuts;
5879 
5880  int nrows;
5881  int nvars;
5882 
5883  SCIP_AGGRROW* aggrrow;
5884  SCIP_Real* cutcoefs;
5885  SCIP_Real* deltas;
5886  int* cutinds;
5887  int deltassize;
5888  int ndeltas;
5889  SCIP_Real* rowweights;
5890  SCIP_Real* comcutdemands;
5891  SCIP_Real* comdemands;
5892  unsigned int partition;
5893  unsigned int allpartitions;
5894  unsigned int startpartition;
5895  SCIP_Bool useinverted;
5896  SCIP_Bool inverted;
5897  int maxtestdelta;
5898 
5899  int oldncuts = 0; /* to check success of separation for one nodeset */
5900  *cutoff = FALSE;
5901 
5902  assert( effortlevel == MCFEFFORTLEVEL_AGGRESSIVE || effortlevel == MCFEFFORTLEVEL_DEFAULT );
5903  nrows = SCIPgetNLPRows(scip);
5904  nvars = SCIPgetNVars(scip);
5905 
5906  /* get the maximal number of cuts allowed in a separation round */
5907  if( SCIPgetDepth(scip) == 0 )
5908  maxsepacuts = sepadata->maxsepacutsroot;
5909  else
5910  maxsepacuts = sepadata->maxsepacuts;
5911  if( maxsepacuts < 0 )
5912  maxsepacuts = INT_MAX;
5913  else if( effortlevel == MCFEFFORTLEVEL_AGGRESSIVE )
5914  maxsepacuts *= 2;
5915 
5916  /* get the maximal number of deltas to use for cmir separation */
5917  maxtestdelta = sepadata->maxtestdelta;
5918  if( maxtestdelta <= 0 )
5919  maxtestdelta = INT_MAX;
5920  else if( effortlevel == MCFEFFORTLEVEL_AGGRESSIVE )
5921  maxtestdelta *= 2;
5922 
5923  /* Our system has the following form:
5924  * (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
5925  * (2) \sum_{k \in K} f_a^k - c_a x_a <= 0 for all a \in A
5926  *
5927  * The partitioning yields two clusters:
5928  *
5929  * A^+
5930  * cluster S ------> cluster T
5931  * <------
5932  * A^-
5933  *
5934  * Now we look at all commodities in which we have to route flow from T to S:
5935  * K^+ = {k : d^k_S := sum_{v \in S} d_v^k > 0}
5936  *
5937  * Then, the base constraint of the c-MIR cut is the sum of those flow conservation constraints and the
5938  * capacity constraints for arcs A^-:
5939  *
5940  * 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
5941  * + sum_{a \in A^-} sum_{k \in K} f_a^k - c_a x_a <= 0
5942  */
5943 
5944  deltassize = 16;
5945  SCIP_CALL( SCIPallocBufferArray(scip, &deltas, deltassize) );
5946  SCIP_CALL( SCIPallocBufferArray(scip, &rowweights, nrows) );
5947  SCIP_CALL( SCIPallocBufferArray(scip, &comcutdemands, ncommodities) );
5948  SCIP_CALL( SCIPallocBufferArray(scip, &comdemands, ncommodities) );
5949  SCIP_CALL( SCIPallocBufferArray(scip, &cutcoefs, nvars) );
5950  SCIP_CALL( SCIPallocBufferArray(scip, &cutinds, nvars) );
5951 
5952  /* create empty aggregation row */
5953  SCIP_CALL( SCIPaggrRowCreate(scip, &aggrrow) );
5954 
5955  if( nodepartition == NULL )
5956  {
5957  /* loop over all nodes and generate single-node cuts */
5958  startpartition = 0;
5959  allpartitions = (unsigned int) nnodes;
5960  SCIPdebugMsg(scip, "maxtestdelta: %d, maxsepacuts: %d, nnodes: %d \n", maxtestdelta, maxsepacuts, nnodes);
5961  }
5962  else
5963  {
5964  /* loop over all possible partitions of the clusters;
5965  * cluster i is in S iff bit i of 'partition' is 1
5966  */
5967  int nclusters = nodepartition->nclusters;
5968 
5969  assert((unsigned int)nclusters <= 8*sizeof(unsigned int));
5970  SCIPdebugMsg(scip, "maxtestdelta: %d, maxsepacuts: %d, nclusters: %d \n", maxtestdelta, maxsepacuts, nclusters);
5971 
5972  /* We fix the last cluster to belong to partition T. In the undirected case, this is sufficient,
5973  * because S-T is equivalent to T-S. In the directed case, the loop below is conducted two times,
5974  * one with regular S-T and one with inverted partitions.
5975  */
5976  startpartition = 1;
5977  allpartitions = (1 << (nclusters-1)); /*lint !e701*/
5978  }
5979  useinverted = (mcfnetwork->modeltype == SCIP_MCFMODELTYPE_DIRECTED);
5980 
5981  for( partition = startpartition; partition <= allpartitions-1 && !SCIPisStopped(scip) && *ncuts < maxsepacuts && !*cutoff; partition++ )
5982  {
5983  int v;
5984  int a;
5985  int k;
5986  int d;
5987  int nnodesinS;
5988  SCIP_Bool success;
5989  /* variables for flowcutset separation */
5990  SCIP_Real baserhs;
5991  SCIP_Real bestdelta = 1.0;
5992  SCIP_Real bestefficacy;
5993  SCIP_Real f0;
5994 
5995  if( sepadata->checkcutshoreconnectivity )
5996  {
5997  if( nodepartition != NULL && !nodepartitionIsConnected(scip, mcfnetwork, nodepartition, partition ) )
5998  {
5999  /* if S or T are not connected, it is very likely that there is a cut in our cluster partition
6000  that gives dominating inequalities
6001  */
6002  SCIPdebugMsg(scip, "generating cluster cuts for partition 0x%x \n", partition );
6003  SCIPdebugMsg(scip, " -> either shore S or shore T is not connected - skip partition.\n");
6004  continue;
6005  }
6006  }
6007 
6008  for( inverted = FALSE; inverted <= useinverted && !*cutoff; inverted++ )
6009  {
6010  if( nodepartition == NULL )
6011  {
6012  SCIPdebugMsg(scip, "generating single-node cuts for node %u (inverted: %u)\n", partition, inverted);
6013  }
6014  else
6015  {
6016  SCIPdebugMsg(scip, "generating cluster cuts for partition 0x%x (inverted: %u)\n", partition, inverted);
6017  }
6018 
6019 #ifdef OUTPUTGRAPH
6020  SCIP_CALL( outputGraph(scip, mcfnetwork, nodepartition, inverted, partition) );
6021 #endif
6022  /* Check if S and T are connected via union find algorithm.
6023  * We do not check this for single node cuts. First, this can be expensive since
6024  * in this case the graph still contains all nodes. Second, we may not find a dominating cut,
6025  * because we may not have the corresponding dominating node set in our cluster partition.
6026  */
6027 
6028  /* clear memory */
6029  BMSclearMemoryArray(rowweights, nrows);
6030  BMSclearMemoryArray(comcutdemands, ncommodities);
6031  BMSclearMemoryArray(comdemands, ncommodities);
6032 
6033  baserhs = 0.0;
6034  ndeltas = 0;
6035 
6036  /* Identify commodities with positive T -> S demand */
6037  nnodesinS = 0;
6038  for( v = 0; v < nnodes; v++ )
6039  {
6040  /* check if node belongs to S */
6041  if( !nodeInPartition(nodepartition, partition, inverted, v) )
6042  {
6043  /* node does not belong to S */
6044  continue;
6045  }
6046  nnodesinS++;
6047  /* update commodity demand */
6048  for( k = 0; k < ncommodities; k++ )
6049  {
6050  SCIP_Real rhs;
6051 
6052  if( nodeflowrows[v][k] == NULL )
6053  continue;
6054 
6055  if( nodeflowscales[v][k] > 0.0 )
6056  rhs = SCIProwGetRhs(nodeflowrows[v][k]) - SCIProwGetConstant(nodeflowrows[v][k]);
6057  else
6058  rhs = SCIProwGetLhs(nodeflowrows[v][k]) - SCIProwGetConstant(nodeflowrows[v][k]);
6059  if( nodeflowinverted[v][k] )
6060  rhs *= -1.0;
6061 
6062  comcutdemands[k] += rhs * nodeflowscales[v][k];
6063  }
6064  }
6065  assert (1 <= nnodesinS && nnodesinS <= nnodes-1);
6066 
6067  /* ignore cuts with only a single node in S or in T, since these have
6068  * already been tried as single node cuts
6069  */
6070  if( sepadata->separatesinglenodecuts && nodepartition != NULL && (nnodesinS == 1 || nnodesinS == nnodes-1) )
6071  {
6072  SCIPdebugMsg(scip, " -> shore S or T has only one node - skip partition.\n");
6073  break;
6074  }
6075 
6076  /* check if there is at least one useful commodity */
6077  if( modeltype == SCIP_MCFMODELTYPE_DIRECTED )
6078  {
6079  for( k = 0; k < ncommodities; k++ )
6080  {
6081  /* in the directed case, use commodities with positive demand (negative -d_k) */
6082  SCIPdebugMsg(scip, " -> commodity %d: directed cutdemand=%g\n", k, comcutdemands[k]);
6083  if( SCIPisNegative(scip, comcutdemands[k]) )
6084  break;
6085  }
6086  }
6087  else
6088  {
6089  for( k = 0; k < ncommodities; k++ )
6090  {
6091  /* in the undirected case, use commodities with non-zero demand */
6092  SCIPdebugMsg(scip, " -> commodity %d: undirected cutdemand=%g\n", k, comcutdemands[k]);
6093  if( !SCIPisZero(scip, comcutdemands[k]) )
6094  break;
6095  }
6096  }
6097  if( k == ncommodities )
6098  continue;
6099 
6100  /* set weights of capacity rows that go from T to S, i.e., a \in A^- */
6101  for( a = 0; a < narcs; a++ )
6102  {
6103  SCIP_COL** rowcols;
6104  SCIP_Real* rowvals;
6105  SCIP_Real feasibility;
6106  int rowlen;
6107  int r;
6108  int j;
6109 
6110  assert(arcsources[a] < nnodes);
6111  assert(arctargets[a] < nnodes);
6112 
6113  /* check if this is an arc of our cut */
6114  if( modeltype == SCIP_MCFMODELTYPE_DIRECTED )
6115  {
6116  /* in the directed case, check if arc goes from T to S */
6117  if( nodeInPartition(nodepartition, partition, inverted, arcsources[a]) ||
6118  !nodeInPartition(nodepartition, partition, inverted, arctargets[a]) )
6119  {
6120  /* arc has source in S or target in T */
6121  continue;
6122  }
6123  }
6124  else
6125  {
6126  /* in the undirected case, check if the arc has endpoints in S and T */
6127  if( nodeInPartition(nodepartition, partition, inverted, arcsources[a]) &&
6128  nodeInPartition(nodepartition, partition, inverted, arctargets[a]) )
6129  {
6130  /* both endpoints are in S */
6131  continue;
6132  }
6133  if( !nodeInPartition(nodepartition, partition, inverted, arcsources[a]) &&
6134  !nodeInPartition(nodepartition, partition, inverted, arctargets[a]) )
6135  {
6136  /* both endpoints are in T */
6137  continue;
6138  }
6139  }
6140 
6141  /* arc might be uncapacitated */
6142  if( arccapacityrows[a] == NULL )
6143  continue;
6144 
6145  /* use capacity row in c-MIR cut */
6146  r = SCIProwGetLPPos(arccapacityrows[a]);
6147  assert(r < nrows);
6148  if( r == -1 ) /* row might have been removed from LP in the meantime */
6149  continue;
6150  assert(rowweights[r] == 0.0);
6151 
6152  /* if one of the arc nodes is unknown, we only use the capacity row if it does not have slack,
6153  * otherwise, we discard it if the slack is too large
6154  */
6155  feasibility = SCIPgetRowSolFeasibility(scip, arccapacityrows[a], sol);
6156  if( arcsources[a] == -1 || arctargets[a] == -1 )
6157  {
6158  if( SCIPisFeasPositive(scip, feasibility) )
6159  continue;
6160  }
6161  else
6162  {
6163  SCIP_Real maxcoef;
6164 
6165  maxcoef = SCIPgetRowMaxCoef(scip, arccapacityrows[a]);
6166  assert(maxcoef > 0.0);
6167  if( SCIPisFeasGT(scip, feasibility/maxcoef, MAXCAPACITYSLACK) )
6168  continue;
6169  }
6170 
6171  rowweights[r] = arccapacityscales[a];
6172  SCIPdebugMsg(scip, " -> arc %d, r=%d, capacity row <%s>: weight=%g slack=%g dual=%g\n", a, r, SCIProwGetName(arccapacityrows[a]), rowweights[r],
6173  SCIPgetRowFeasibility(scip, arccapacityrows[a]), SCIProwGetDualsol(arccapacityrows[a]));
6174  /*SCIPdebug( SCIP_CALL(SCIPprintRow(scip, arccapacityrows[a], NULL)) );*/
6175 
6176  if( sepadata->separateflowcutset )
6177  {
6178  if( rowweights[r] > 0.0 )
6179  baserhs += rowweights[r] * (