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