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