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