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-2024 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
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"
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 */
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 */
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 */
193struct 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 */
217struct 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};
259typedef 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 */
262struct 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};
268typedef struct nodepair NODEPAIRENTRY;
269
270/** nodepair priority queue */
271struct nodepairqueue
272{
273 SCIP_PQUEUE* pqueue; /**< priority queue of elements */
274 NODEPAIRENTRY* nodepairs; /**< elements on the heap */
275};
276typedef struct nodepairqueue NODEPAIRQUEUE;
277
278/** partitioning of the nodes into clusters */
279struct 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};
287typedef 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 */
303static
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 */
329static
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 */
380static
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 */
675static
676void 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 */
721static
722void 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 */
838static
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 */
852static
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 {
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 */
1068static
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 */
1438static
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 */
1462static
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 */
1515static
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 */
1661static
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 */
1726static
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 */
1810static
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 */
1942static
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 */
2059static
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 */
2253static
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 */
2424static
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 */
2509static
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 */
2710static
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 */
2986static
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 */
3003static
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 */
3102static
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 */
3442static
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 */
3760static
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 */
4109static
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 */
4240static
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 */
4549static
4550SCIP_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 {
4613 nbinflowvars++;
4614 break;
4616 nintflowvars++;
4617 break;
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 {
4659 nbincapvars++;
4660 break;
4662 nintcapvars++;
4663 break;
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 */
4715static
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 */
4729static
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 */
4747static
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 */
4774static
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 */
4790static
4791SCIP_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 */
4802static
4803SCIP_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 */
4843static
4844SCIP_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 */
4877static
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;
4894 SCIP_ROW*** nodeflowrows;
4895 SCIP_Real** nodeflowscales;
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 */
5190static
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 */
5206static
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 */
5218static
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 */
5235static
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 */
5245static
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 */
5256static
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 */
5404static
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 */
5421static
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 */
5451static
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);
5533 return (ncomponents == 2);
5534}
5535