1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2021 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15
16 /**@file cons_stp.c
17 * @brief Constraint handler for Steiner problems
18 * @author Gerald Gamrath
19 * @author Daniel Rehfeldt
20 * @author Michael Winkler
21 *
22 * This file checks solutions for feasibility and separates violated model constraints. For more details see \ref STP_CONS page.
23 *
24 * @page STP_CONS Separating violated constraints
25 *
26 * In this file a constraint handler checking solutions for feasibility and separating violated model constraints is implemented.
27 * The separation problem for the cut inequalities described in \ref STP_PROBLEM can be solved by a max-flow algorithm in
28 * polynomial time. Regarding the variable values of a given LP solution as capacities on the edges, one can check for each
29 * \f$ t \in T \setminus \{r\} \f$, with \f$ r \f$ being the root, whether the minimal \f$ (r, t) \f$-cut is less than one. In this case,
30 * a violated cut inequality has been found, otherwise none exists. In order to calculate such a minimal cut an adaptation of Hao
31 * and Orlin's preflow-push algorithm (see A Faster Algorithm for Finding the Minimum Cut in a Directed Graph) is used. Furthermore, the file implements a dual ascent heuristic, based on a concept described
32 * in "A dual ascent approach for Steiner tree problems on a directed graph." by R. Wong.
33 *
34 */
35
36 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
37
38 #include <assert.h>
39 #include <string.h>
40
41 #include "cons_stp.h"
42 #include "probdata_stp.h"
43 #include "grph.h"
44 #include "heur_prune.h"
45 #include "heur_ascendprune.h"
46 #include "portab.h"
47 #include "branch_stp.h"
48
49 #include "scip/scip.h"
50 #include "scip/misc.h"
51 #include "scip/cons_linear.h"
52 #include <time.h>
53 #if 0
54 #ifdef WITH_UG
55 #define ADDCUTSTOPOOL 1
56 #else
57 #define ADDCUTSTOPOOL 0
58 #endif
59 #endif
60
61 #define ADDCUTSTOPOOL 0
62
63 #define Q_NULL -1 /* NULL element of queue/list */
64
65 /**@name Constraint handler properties
66 *
67 * @{
68 */
69
70 #define CONSHDLR_NAME "stp"
71 #define CONSHDLR_DESC "steiner tree constraint handler"
72 #define CONSHDLR_SEPAPRIORITY 0 /**< priority of the constraint handler for separation */
73 #define CONSHDLR_ENFOPRIORITY 0 /**< priority of the constraint handler for constraint enforcing */
74 #define CONSHDLR_CHECKPRIORITY 9999999 /**< priority of the constraint handler for checking feasibility */
75 #define CONSHDLR_SEPAFREQ 1 /**< frequency for separating cuts; zero means to separate only in the root node */
76 #define CONSHDLR_PROPFREQ 0 /**< frequency for propagating domains; zero means only preprocessing propagation */
77 #define CONSHDLR_EAGERFREQ 1 /**< frequency for using all instead of only the useful constraints in separation,
78 * propagation and enforcement, -1 for no eager evaluations, 0 for first only */
79 #define CONSHDLR_DELAYSEPA FALSE /**< should separation method be delayed, if other separators found cuts? */
80 #define CONSHDLR_DELAYPROP FALSE /**< should propagation method be delayed, if other propagators found reductions? */
81 #define CONSHDLR_NEEDSCONS TRUE /**< should the constraint handler be skipped, if no constraints are available? */
82
83 #define DEFAULT_MAXROUNDS 10 /**< maximal number of separation rounds per node (-1: unlimited) */
84 #define DEFAULT_MAXROUNDSROOT -1 /**< maximal number of separation rounds in the root node (-1: unlimited) */
85 #define DEFAULT_MAXSEPACUTS INT_MAX /**< maximal number of cuts separated per separation round */
86 #define DEFAULT_MAXSEPACUTSROOT INT_MAX /**< maximal number of cuts separated per separation round in the root node */
87
88
89 #define CONSHDLR_PROP_TIMING SCIP_PROPTIMING_BEFORELP
90
91 #define DEFAULT_BACKCUT FALSE /**< Try Back-Cuts FALSE*/
92 #define DEFAULT_CREEPFLOW TRUE /**< Use Creep-Flow */
93 #define DEFAULT_DISJUNCTCUT FALSE /**< Only disjunct Cuts FALSE */
94 #define DEFAULT_NESTEDCUT FALSE /**< Try Nested-Cuts FALSE*/
95 #define DEFAULT_FLOWSEP TRUE /**< Try Flow-Cuts */
96
97 #define DEFAULT_DAMAXDEVIATION 0.25 /**< max deviation for dual ascent */
98 #define DA_MAXDEVIATION_LOWER 0.01 /**< lower bound for max deviation for dual ascent */
99 #define DA_MAXDEVIATION_UPPER 0.9 /**< upper bound for max deviation for dual ascent */
100 #define DA_EPS (5.0 * 1e-7)
101
102 /* *
103 #define FLOW_FACTOR 100000
104 #define CREEP_VALUE 1 this is the original value todo check what is better
105 */
106
107 #define FLOW_FACTOR 1000000
108 #define CREEP_VALUE 10
109
110 /* do depth-first search */
111 #define DFS
112
113
114 #ifdef BITFIELDSARRAY
115 #define ARRLENGTH 32
116 #define SetBit(Arr, pos) ( Arr[(pos/ARRLENGTH)] |= (1 << (pos%ARRLENGTH)) )
117 #define CleanBit(Arr, pos) ( Arr[(pos/ARRLENGTH)] &= ~(1 << (pos%ARRLENGTH)) )
118 #define BitTrue(Arr, pos) ( Arr[(pos/ARRLENGTH)] & (1 << (pos%ARRLENGTH)) )
119 #endif
120
121
122 /**@} */
123
124 /*
125 * Data structures
126 */
127
128 /** @brief Constraint data for \ref cons_stp.c "Stp" constraints */
129 struct SCIP_ConsData
130 {
131 GRAPH* graph; /**< graph data structure */
132 };
133
134 /** @brief Constraint handler data for \ref cons_stp.c "Stp" constraint handler */
135 struct SCIP_ConshdlrData
136 {
137 SCIP_Bool backcut; /**< should backcuts be applied? */
138 SCIP_Bool creepflow; /**< should creepflow cuts be applied? */
139 SCIP_Bool disjunctcut; /**< should disjunction cuts be applied? */
140 SCIP_Bool nestedcut; /**< should nested cuts be applied? */
141 SCIP_Bool flowsep; /**< should flow separation be applied? */
142 int maxrounds; /**< maximal number of separation rounds per node (-1: unlimited) */
143 int maxroundsroot; /**< maximal number of separation rounds in the root node (-1: unlimited) */
144 int maxsepacuts; /**< maximal number of cuts separated per separation round */
145 int maxsepacutsroot; /**< maximal number of cuts separated per separation round in the root node */
146 };
147
148
149 /**@name Local methods
150 *
151 * @{
152 */
153
154 /** returns whether node realtail is active or leads to active node other than dfsbase */
155 static
is_active(const int * active,int realtail,int dfsbase)156 SCIP_Bool is_active(
157 const int* active, /**< active nodes array */
158 int realtail, /**< vertex to start from */
159 int dfsbase /**< DFS source vertex */
160 )
161 {
162 int curr;
163
164 for( curr = active[realtail]; curr != 0 && curr != dfsbase + 1; curr = active[curr - 1] )
165 {
166 assert(curr >= 0);
167 }
168
169 return (curr == 0);
170 }
171
172
173
174 /** add a cut */
175 static
cut_add(SCIP * scip,SCIP_CONSHDLR * conshdlr,const GRAPH * g,const SCIP_Real * xval,int * capa,const int updatecapa,int * ncuts,SCIP_Bool local,SCIP_Bool * success)176 SCIP_RETCODE cut_add(
177 SCIP* scip, /**< SCIP data structure */
178 SCIP_CONSHDLR* conshdlr, /**< constraint handler */
179 const GRAPH* g, /**< graph data structure */
180 const SCIP_Real* xval, /**< edge values */
181 int* capa, /**< edges capacities (scaled) */
182 const int updatecapa, /**< update capacities? */
183 int* ncuts, /**< pointer to store number of cuts */
184 SCIP_Bool local, /**< is the cut local? */
185 SCIP_Bool* success /**< pointer to store whether add cut be added */
186 )
187 {
188 SCIP_ROW* row;
189 SCIP_VAR** vars = SCIPprobdataGetVars(scip);
190 SCIP_Real sum = 0.0;
191 SCIP_Bool inccapa = FALSE;
192 unsigned int i;
193 int* gmark = g->mark;
194 int* ghead = g->head;
195 int* gtail = g->tail;
196 unsigned int nedges = (unsigned int) g->edges;
197
198 assert(g->knots > 0);
199
200 (*success) = FALSE;
201
202 assert(g != NULL);
203 assert(scip != NULL);
204
205 SCIP_CALL( SCIPcreateEmptyRowConshdlr(scip, &row, conshdlr, "twocut", 1.0, SCIPinfinity(scip), local, FALSE, TRUE) );
206
207 SCIP_CALL( SCIPcacheRowExtensions(scip, row) );
208
209 assert(gmark[g->source]);
210
211 for( i = 0; i < nedges; i++ )
212 {
213 if( gmark[gtail[i]] && !gmark[ghead[i]] ) // todo bool array?
214 {
215 if( updatecapa )
216 {
217 if( capa[i] < FLOW_FACTOR )
218 inccapa = TRUE;
219
220 SCIPdebugMessage("set capa[%d] from %6d to %6d\n", i, capa[i], FLOW_FACTOR);
221 capa[i] = FLOW_FACTOR;
222
223 if( !inccapa )
224 {
225 SCIP_CALL(SCIPflushRowExtensions(scip, row));
226 SCIP_CALL(SCIPreleaseRow(scip, &row));
227 return SCIP_OKAY;
228 }
229 }
230
231 if( xval != NULL )
232 {
233 sum += xval[i];
234
235 if( SCIPisFeasGE(scip, sum, 1.0) )
236 {
237 SCIP_CALL(SCIPflushRowExtensions(scip, row));
238 SCIP_CALL(SCIPreleaseRow(scip, &row));
239 return SCIP_OKAY;
240 }
241 }
242 SCIP_CALL(SCIPaddVarToRow(scip, row, vars[i], 1.0));
243 }
244 }
245
246 assert(sum < 1.0);
247
248 SCIP_CALL( SCIPflushRowExtensions(scip, row) );
249
250 /* checks whether cut is sufficiently violated */
251 if( SCIPisCutEfficacious(scip, NULL, row) )
252 {
253 SCIP_Bool infeasible;
254
255 SCIPdebug( SCIP_CALL( SCIPprintRow(scip, row, NULL) ) );
256
257 SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
258
259 #if ADDCUTSTOPOOL
260 /* if at root node, add cut to pool */
261 if( !infeasible )
262 SCIP_CALL( SCIPaddPoolCut(scip, row) );
263 #endif
264 (*ncuts)++;
265 (*success) = TRUE;
266 }
267
268 SCIP_CALL( SCIPreleaseRow(scip, &row) );
269
270 return SCIP_OKAY;
271 }
272
273
274
275 static
graph_next_term(const GRAPH * g,int terms,int * term,const int * w,const SCIP_Bool firstrun)276 int graph_next_term(
277 const GRAPH* g, /**< graph data structure */
278 int terms, /**< number of terminals */
279 int* term, /**< terminal array */
280 const int* w, /**< awake level */
281 const SCIP_Bool firstrun /**< first run? */
282 )
283 {
284 int i;
285 int k;
286 int t;
287 int wmax;
288 int mindist = g->knots + 1;
289
290 assert(term != NULL);
291
292 if( firstrun ) // todo randomize?
293 {
294 assert(w[term[terms - 1]] == 0);
295 return term[terms - 1];
296 }
297
298 k = -1;
299
300 for( i = 0; (i < terms); i++ )
301 {
302 assert(w[term[i]] >= 0);
303
304 if( w[term[i]] == 0 )
305 {
306 assert(g->mincut_dist[term[i]] < g->knots + 1);
307
308 if( g->mincut_dist[term[i]] < mindist )
309 {
310 k = i;
311 mindist = g->mincut_dist[term[i]];
312 }
313 }
314 }
315
316 if( k == -1 )
317 {
318 wmax = 0;
319
320 for( i = 0; (i < terms); i++ )
321 {
322 if( w[term[i]] > wmax )
323 {
324 k = i;
325 wmax = w[term[i]];
326 mindist = g->mincut_dist[term[i]];
327 }
328 else if( w[term[i]] == wmax && g->mincut_dist[term[i]] < mindist )
329 {
330 assert(wmax != 0);
331
332 k = i;
333 mindist = g->mincut_dist[term[i]];
334 }
335 }
336 }
337
338 assert(k >= 0);
339 assert(k < terms);
340
341 t = term[k];
342 term[k] = term[terms - 1];
343
344 return t;
345 }
346
347 static
set_capacity(const GRAPH * g,const SCIP_Bool creep_flow,const int flip,int * capa,const SCIP_Real * xval)348 void set_capacity(
349 const GRAPH* g, /**< graph data structure */
350 const SCIP_Bool creep_flow, /**< creep flow? */
351 const int flip, /**< reverse the flow? */
352 int* capa, /**< edges capacities (scaled) */
353 const SCIP_Real* xval /**< edge values */
354 )
355 {
356 int k;
357 int krev;
358 int nedges = g->edges;
359
360 assert(g != NULL);
361 assert(xval != NULL);
362
363 for( k = 0; k < nedges; k += 2 )
364 {
365 krev = k + 1;
366 if( !flip )
367 {
368 capa[k] = (int)(xval[k ]
369 * FLOW_FACTOR + 0.5);
370 capa[krev] = (int)(xval[krev]
371 * FLOW_FACTOR + 0.5);
372 }
373 else
374 {
375 capa[k] = (int)(xval[krev]
376 * FLOW_FACTOR + 0.5);
377 capa[krev] = (int)(xval[k ]
378 * FLOW_FACTOR + 0.5);
379 }
380
381 if( creep_flow )
382 {
383 capa[k] += CREEP_VALUE;
384 capa[krev] += CREEP_VALUE;
385 }
386 }
387 }
388
389 /** separate */
390 static
sep_flow(SCIP * scip,SCIP_CONSHDLR * conshdlr,SCIP_CONSHDLRDATA * conshdlrdata,SCIP_CONSDATA * consdata,int maxcuts,int * ncuts)391 SCIP_RETCODE sep_flow(
392 SCIP* scip, /**< SCIP data structure */
393 SCIP_CONSHDLR* conshdlr, /**< constraint handler */
394 SCIP_CONSHDLRDATA* conshdlrdata, /**< constraint handler data */
395 SCIP_CONSDATA* consdata, /**< constraint data */
396 int maxcuts, /**< maximal number of cuts */
397 int* ncuts /**< pointer to store number of cuts */
398 )
399 {
400 GRAPH* g;
401 SCIP_VAR** vars;
402 SCIP_ROW* row = NULL;
403 SCIP_Real* xval;
404 SCIP_Real sum;
405 int i;
406 int k;
407 int j;
408 int ind;
409 int layer;
410 int count = 0;
411 unsigned int flowsep;
412
413 assert(scip != NULL);
414 assert(conshdlr != NULL);
415 assert(conshdlrdata != NULL);
416
417 vars = SCIPprobdataGetVars(scip);
418 flowsep = conshdlrdata->flowsep;
419
420 /* get the graph */
421 g = consdata->graph;
422 assert(g != NULL);
423
424 xval = SCIPprobdataGetXval(scip, NULL);
425 assert(xval != NULL);
426
427 for(i = 0; i < g->knots; i++)
428 {
429 for(layer = 0; layer < g->layers; layer++)
430 {
431 /* continue at root */
432 if( i == g->source )
433 continue;
434
435 /* at terminal: input sum == 1
436 * basically a cut (starcut))
437 */
438 if( g->term[i] == layer )
439 {
440 sum = 0.0;
441
442 for( k = g->inpbeg[i]; k != EAT_LAST; k = g->ieat[k] )
443 {
444 ind = layer * g->edges + k;
445 sum += (xval != NULL) ? xval[ind] : 0.0;
446 }
447
448 if( !SCIPisFeasEQ(scip, sum, 1.0) )
449 {
450 SCIP_Bool infeasible;
451
452 SCIP_CALL( SCIPcreateEmptyRowConshdlr(scip, &row, conshdlr, "term", 1.0,
453 1.0, FALSE, FALSE, TRUE) );
454
455 SCIP_CALL( SCIPcacheRowExtensions(scip, row) );
456
457 for(k = g->inpbeg[i]; k != EAT_LAST; k = g->ieat[k])
458 {
459 ind = layer * g->edges + k;
460
461 SCIP_CALL( SCIPaddVarToRow(scip, row, vars[ind], 1.0) );
462 }
463
464 SCIP_CALL( SCIPflushRowExtensions(scip, row) );
465
466 SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
467
468 #if ADDCUTSTOPOOL
469 /* add cut to pool */
470 if( !infeasible )
471 SCIP_CALL( SCIPaddPoolCut(scip, row) );
472 #endif
473
474 count++;
475
476 SCIP_CALL( SCIPreleaseRow(scip, &row) );
477
478 if( *ncuts + count >= maxcuts )
479 goto TERMINATE;
480 }
481 }
482
483 /* flow cuts disabled? */
484 if( !flowsep )
485 continue;
486
487 /* the value of each outgoing edge needs to be smaller than the sum of the ingoing edges */
488 for( j = g->outbeg[i]; j != EAT_LAST; j = g->oeat[j] )
489 {
490 ind = layer * g->edges + j;
491 sum = (xval != NULL) ? -xval[ind] : -1.0;
492
493 for( k = g->inpbeg[i]; k != EAT_LAST; k = g->ieat[k] )
494 {
495 ind = layer * g->edges + k;
496 sum += (xval != NULL) ? xval[ind] : 0.0;
497 }
498 if( SCIPisFeasNegative(scip, sum) )
499 {
500 SCIP_Bool infeasible;
501
502 SCIP_CALL( SCIPcreateEmptyRowConshdlr(scip, &row, conshdlr, "flow", 0.0, SCIPinfinity(scip),
503 FALSE, FALSE, TRUE) );
504
505 SCIP_CALL( SCIPcacheRowExtensions(scip, row) );
506
507 ind = layer * g->edges + j;
508
509 SCIP_CALL( SCIPaddVarToRow(scip, row, vars[ind], -1.0) );
510
511 for( k = g->inpbeg[i]; k != EAT_LAST; k = g->ieat[k] )
512 {
513 ind = layer * g->edges + k;
514
515 SCIP_CALL( SCIPaddVarToRow(scip, row, vars[ind], 1.0) );
516 }
517
518 SCIP_CALL( SCIPflushRowExtensions(scip, row) );
519
520 SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
521
522 #if ADDCUTSTOPOOL
523 /* add cut to pool */
524 if( !infeasible )
525 SCIP_CALL( SCIPaddPoolCut(scip, row) );
526 #endif
527
528 count++;
529
530 SCIP_CALL( SCIPreleaseRow(scip, &row) );
531
532 if( *ncuts + count >= maxcuts )
533 goto TERMINATE;
534 }
535 }
536
537 /* consider only non terminals */
538 if( g->term[i] == layer )
539 continue;
540
541 /* input of a vertex has to be <= 1.0 */
542 sum = 0.0;
543
544 for( k = g->inpbeg[i]; k != EAT_LAST; k = g->ieat[k] )
545 {
546 ind = layer * g->edges + k;
547 sum += (xval != NULL) ? xval[ind] : 1.0;
548 }
549 if( SCIPisFeasGT(scip, sum, 1.0) )
550 {
551 SCIP_Bool infeasible;
552
553 SCIP_CALL( SCIPcreateEmptyRowConshdlr(scip, &row, conshdlr, "infl", -SCIPinfinity(scip),
554 1.0, FALSE, FALSE, TRUE) );
555
556 SCIP_CALL( SCIPcacheRowExtensions(scip, row) );
557
558 for( k = g->inpbeg[i]; k != EAT_LAST; k = g->ieat[k] )
559 {
560 ind = layer * g->edges + k;
561
562 SCIP_CALL( SCIPaddVarToRow(scip, row, vars[ind], 1.0) );
563 }
564
565 SCIP_CALL( SCIPflushRowExtensions(scip, row) );
566
567 SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
568
569 #if ADDCUTSTOPOOL
570 /* if at root node, add cut to pool */
571 if( !infeasible )
572 SCIP_CALL( SCIPaddPoolCut(scip, row) );
573 #endif
574
575 count++;
576
577 SCIP_CALL( SCIPreleaseRow(scip, &row) );
578
579 if( *ncuts + count >= maxcuts )
580 goto TERMINATE;
581 }
582
583 /* incoming flow <= outgoing flow */
584 sum = 0.0;
585
586 for( k = g->inpbeg[i]; k != EAT_LAST; k = g->ieat[k] )
587 {
588 ind = layer * g->edges + k;
589 sum -= (xval != NULL) ? xval[ind] : 1.0;
590 }
591 for( k = g->outbeg[i]; k != EAT_LAST; k = g->oeat[k] )
592 {
593 ind = layer * g->edges + k;
594 sum += (xval != NULL) ? xval[ind] : 0.0;
595 }
596 if( SCIPisFeasNegative(scip, sum) )
597 {
598 SCIP_Bool infeasible;
599
600 SCIP_CALL( SCIPcreateEmptyRowConshdlr(scip, &row, conshdlr, "bala", 0.0,
601 (g->terms == 2) ? 0.0 : SCIPinfinity(scip), FALSE, FALSE, TRUE) );
602
603 SCIP_CALL( SCIPcacheRowExtensions(scip, row) );
604
605 for( k = g->inpbeg[i]; k != EAT_LAST; k = g->ieat[k] )
606 {
607 ind = layer * g->edges + k;
608
609 SCIP_CALL( SCIPaddVarToRow(scip, row, vars[ind], -1.0) );
610 }
611 for( k = g->outbeg[i]; k != EAT_LAST; k = g->oeat[k] )
612 {
613 ind = layer * g->edges + k;
614
615 SCIP_CALL( SCIPaddVarToRow(scip, row, vars[ind], 1.0) );
616 }
617
618 SCIP_CALL( SCIPflushRowExtensions(scip, row) );
619
620 SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
621
622 #if ADDCUTSTOPOOL
623 /* if at root node, add cut to pool */
624 if( !infeasible )
625 SCIP_CALL( SCIPaddPoolCut(scip, row) );
626 #endif
627
628 count++;
629
630 SCIP_CALL( SCIPreleaseRow(scip, &row) );
631
632 if( *ncuts + count >= maxcuts )
633 goto TERMINATE;
634 }
635 }
636 }
637
638 TERMINATE:
639 SCIPdebugMessage("In/Out Separator: %d Inequalities added\n", count);
640
641 *ncuts += count;
642
643 return SCIP_OKAY;
644 }
645
646 /** separate 2-cuts */
647 static
sep_2cut(SCIP * scip,SCIP_CONSHDLR * conshdlr,SCIP_CONSHDLRDATA * conshdlrdata,SCIP_CONSDATA * consdata,const int * termorg,int maxcuts,int * ncuts)648 SCIP_RETCODE sep_2cut(
649 SCIP* scip, /**< SCIP data structure */
650 SCIP_CONSHDLR* conshdlr, /**< constraint handler */
651 SCIP_CONSHDLRDATA* conshdlrdata, /**< constraint handler data */
652 SCIP_CONSDATA* consdata, /**< constraint data */
653 const int* termorg, /**< original terminals or NULL */
654 int maxcuts, /**< maximal number of cuts */
655 int* ncuts /**< pointer to store number of cuts */
656 )
657 {
658 #if 0
659 const SCIP_Bool nested_cut = conshdlrdata->nestedcut;
660 const SCIP_Bool back_cut = conshdlrdata->backcut;
661 const SCIP_Bool creep_flow = conshdlrdata->creepflow;
662 const SCIP_Bool disjunct_cut = conshdlrdata->disjunctcut;
663 #endif
664 /* we do not longer support any other flow as they slow everything down and are of little use anyway todo remove user parameter */
665 const SCIP_Bool flowsep = conshdlrdata->flowsep;
666 const SCIP_Bool nested_cut = FALSE;
667 const SCIP_Bool creep_flow = TRUE;
668 const SCIP_Bool disjunct_cut = FALSE;
669 const SCIP_Bool intree = (SCIPgetDepth(scip) > 0);
670
671 SCIP_VAR** vars;
672 GRAPH* g;
673 SCIP_Real* xval;
674 int* w;
675 int* capa;
676 int* term;
677 int* start;
678 int* excess;
679 int* rootcut;
680 int* edgearr;
681 int* headarr;
682 int* residual;
683 int* edgecurr;
684 int* headactive;
685 int* edgeflipped;
686 int* headinactive;
687 int i;
688 int k;
689 int e;
690 int root;
691 int head;
692 int count;
693 int terms;
694 int nedges;
695 int nnodes;
696 int newnnodes;
697 int newnedges;
698 int rootcutsize;
699 SCIP_Bool rerun;
700 SCIP_Bool addedcut;
701
702 assert(scip != NULL);
703 assert(conshdlr != NULL);
704 assert(conshdlrdata != NULL);
705
706 g = consdata->graph;
707 assert(g != NULL);
708
709 root = g->source;
710 excess = g->mincut_e;
711 nedges = g->edges;
712 nnodes = g->knots;
713 addedcut = FALSE;
714 residual = g->mincut_r;
715 edgecurr = g->mincut_numb;
716 headactive = g->mincut_head;
717 headinactive = g->mincut_head_inact;
718
719 assert(residual != NULL);
720 assert(edgecurr != NULL);
721 assert(headactive != NULL);
722 assert(headinactive != NULL);
723
724 xval = SCIPprobdataGetXval(scip, NULL);
725 assert(xval != NULL);
726
727 assert(creep_flow == TRUE);
728 assert(nested_cut == FALSE);
729 assert(disjunct_cut == FALSE);
730
731 /* for 2-terminal nets no cuts are necessary if flows are given */
732 if( flowsep && (g->terms == 2) )
733 return SCIP_OKAY;
734
735 SCIP_CALL( SCIPallocBufferArray(scip, &capa, nedges) );
736 SCIP_CALL( SCIPallocBufferArray(scip, &w, nnodes) );
737 SCIP_CALL( SCIPallocBufferArray(scip, &term, g->terms) );
738 SCIP_CALL( SCIPallocBufferArray(scip, &edgearr, nedges) );
739 SCIP_CALL( SCIPallocBufferArray(scip, &headarr, nedges) );
740 SCIP_CALL( SCIPallocBufferArray(scip, &edgeflipped, nedges) );
741 SCIP_CALL( SCIPallocBufferArray(scip, &start, nnodes + 1) );
742 SCIP_CALL( SCIPallocBufferArray(scip, &rootcut, nnodes + 1) );
743
744 #if 0
745 clock_t startt, endt;
746 double cpu_time_used;
747 startt = clock();
748 #endif
749
750 vars = SCIPprobdataGetVars(scip);
751 assert(vars != NULL);
752
753 assert(nedges >= nnodes);
754
755 for( k = 0; k < nnodes; k++ )
756 {
757 w[k] = 0;
758 excess[k] = 0;
759 }
760
761 for( e = 0; e < nedges; e += 2 )
762 {
763 const int erev = e + 1;
764
765 if( intree && SCIPvarGetUbLocal(vars[e]) < 0.5 && SCIPvarGetUbLocal(vars[erev]) < 0.5 )
766 {
767 capa[e] = 0;
768 capa[erev] = 0;
769 residual[e] = 0;
770 residual[erev] = 0;
771
772 headarr[e] = 1;
773 headarr[erev] = 1;
774 }
775 else
776 {
777 capa[e] = (int)(xval[e] * FLOW_FACTOR + 0.5) + CREEP_VALUE;
778 capa[erev] = (int)(xval[erev] * FLOW_FACTOR + 0.5) + CREEP_VALUE;
779 residual[e] = capa[e];
780 residual[erev] = capa[erev];
781
782 headarr[e] = SCIPisFeasLT(scip, xval[e], 1.0) ? 1 : 0;
783 headarr[erev] = SCIPisFeasLT(scip, xval[erev], 1.0) ? 1 : 0;
784 }
785 edgearr[e] = -1;
786 edgearr[erev] = -1;
787 }
788
789 /*
790 * bfs along 0 edges from the root
791 * */
792
793 w[root] = 1;
794 rootcutsize = 0;
795 rootcut[rootcutsize++] = root;
796
797 /* bfs loop */
798 for( i = 0; i < rootcutsize; i++ )
799 {
800 assert(rootcutsize <= nnodes);
801
802 k = rootcut[i];
803
804 assert(k < nnodes);
805
806 /* traverse outgoing arcs */
807 for( e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
808 {
809 head = g->head[e];
810
811 /* head not been added to root cut yet? */
812 if( w[head] == 0 )
813 {
814 if( headarr[e] == 0 )
815 {
816 w[head] = 1;
817 rootcut[rootcutsize++] = head;
818 }
819 else
820 {
821 /* push as much as possible out of perpetually dormant nodes (possibly to other dormant nodes) */
822 assert(w[head] == 0);
823 #if 1 /* for debug */
824 residual[e] = 0;
825 #endif
826 excess[head] += capa[e];
827 }
828 }
829 }
830 }
831
832 i = 0;
833 terms = 0;
834 newnnodes = 0;
835
836 /* fill auxiliary adjacent vertex/edges arrays and get useable terms */
837 for( k = 0; k < nnodes; k++ )
838 {
839 headactive[k] = Q_NULL;
840 headinactive[k] = Q_NULL;
841
842 start[k] = i;
843
844 /* non-dormant node? */
845 if( w[k] == 0 )
846 {
847 newnnodes++;
848 edgecurr[k] = i;
849 for( e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
850 {
851 if( w[g->head[e]] == 0 && capa[e] != 0 )
852 {
853 edgearr[e] = i;
854 residual[i] = capa[e];
855 headarr[i++] = g->head[e];
856 }
857 }
858
859 /* unreachable node? */
860 if( edgecurr[k] == i )
861 {
862 w[k] = 1;
863 newnnodes--;
864 }
865 else if( Is_term(g->term[k]) )
866 {
867 term[terms++] = k;
868 }
869 }
870 else
871 {
872 edgecurr[k] = -1;
873 }
874 }
875
876 newnedges = i;
877 start[nnodes] = i;
878
879 /* initialize edgeflipped */
880 for( e = nedges - 1; e >= 0; e-- )
881 {
882 if( edgearr[e] >= 0 )
883 {
884 i = edgearr[e];
885 edgeflipped[i] = edgearr[flipedge(e)];
886 }
887 }
888
889 SCIPdebugMessage("Cut Pretest: %d eliminations\n", g->terms - terms - 1);
890
891 #if 0
892 endt = clock();
893 cpu_time_used = ((double) (endt - startt)) / CLOCKS_PER_SEC;
894 startt = clock();
895 #endif
896
897 count = 0;
898 rerun = FALSE;
899
900 while( terms > 0 )
901 {
902 if( ((unsigned) terms) % 32 == 0 && SCIPisStopped(scip) )
903 break;
904
905 /* look for reachable terminal */
906 i = graph_next_term(g, terms, term, w, !rerun);
907
908 terms--;
909
910 assert(g->term[i] == 0);
911 assert(g->source != i);
912
913 if( nested_cut && !disjunct_cut )
914 set_capacity(g, creep_flow, 0, capa, xval);
915
916 do
917 {
918 #if 0
919 /* write flow problem in extended dimacs format */
920 FILE *fptr;
921
922 fptr = fopen("flow", "w+");
923 assert(fptr != NULL);
924
925 fprintf(fptr, "p max %d %d \n", nnodes, nedges);
926 fprintf(fptr, "n %d s \n", g->source + 1);
927 fprintf(fptr, "n %d t \n", i + 1);
928
929 for( k = 0; k < nnodes; k++ )
930 {
931 for( e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
932 {
933 fprintf(fptr, "a %d %d %d \n", k + 1, g->head[e] + 1, capa[e]);
934 }
935 }
936
937 fprintf(fptr, "x\n");
938
939 fclose(fptr);
940 #endif
941 // declare cuts on branched-on (artificial) terminals as local
942 const SCIP_Bool localcut = (termorg != NULL && termorg[i] != g->term[i]);
943
944 /* non-trivial cut? */
945 if( w[i] != 1 )
946 {
947 graph_mincut_exec(g, root, i, nnodes, newnedges, rootcutsize, rootcut, capa, w, start, edgeflipped, headarr, rerun);
948
949 /* cut */
950 for( k = nnodes - 1; k >= 0; k-- )
951 g->mark[k] = (w[k] != 0);
952
953 assert(g->mark[root]);
954 }
955 else
956 {
957 SCIP_Real flowsum = 0.0;
958
959 assert(rerun);
960
961 for( e = g->inpbeg[i]; e != EAT_LAST; e = g->ieat[e] )
962 flowsum += xval[e];
963
964 if( SCIPisFeasGE(scip, flowsum, 1.0) )
965 continue;
966
967 for( k = nnodes - 1; k >= 0; k-- )
968 g->mark[k] = TRUE;
969
970 g->mark[i] = FALSE;
971 }
972
973 rerun = TRUE;
974
975 SCIP_CALL( cut_add(scip, conshdlr, g, xval, capa, nested_cut || disjunct_cut, ncuts, localcut, &addedcut) );
976 if( addedcut )
977 {
978 count++;
979
980 if( *ncuts >= maxcuts )
981 goto TERMINATE;
982 }
983 else
984 break;
985 }
986 while( nested_cut ); /* Nested Cut is CONSTANT ! */
987 } /* while terms > 0 */
988
989
990 #if 0
991 endt = clock();
992 cpu_time_used = ((double) (endt - startt)) / CLOCKS_PER_SEC;
993 #endif
994
995 #if 0
996 /*
997 * back cuts currently not supported
998 * */
999 /* back cuts enabled? */
1000 if( back_cut )
1001 {
1002 for( k = 0; k < nnodes; k++ )
1003 w[k] = 0;
1004
1005 if( !nested_cut || disjunct_cut )
1006 set_capacity(g, creep_flow, 1, capa, xval);
1007
1008 terms = tsave;
1009
1010 while( terms > 0 )
1011 {
1012 /* look for reachable terminal */
1013 i = graph_next_term(g, terms, term, w, TRUE);
1014
1015 terms--;
1016
1017 assert(g->term[i] == 0);
1018 assert(g->source != i);
1019
1020 if( nested_cut && !disjunct_cut )
1021 set_capacity(g, creep_flow, 1, capa, xval);
1022
1023 rerun = FALSE;
1024
1025 do
1026 {
1027 graph_mincut_exec(g, i, g->source, nedges, capa, w, start, edgearr, headarr, rerun);
1028
1029 rerun = TRUE;
1030
1031 for( k = 0; k < nnodes; k++ )
1032 {
1033 g->mark[k] = (w[k] != 0) ? 0 : 1; // todo not the other way around??
1034 w[k] = 0;
1035 }
1036
1037 SCIP_CALL( cut_add(scip, conshdlr, g, xval, capa, nested_cut || disjunct_cut, ncuts, &addedcut) );
1038 if( addedcut )
1039 {
1040 count++;
1041
1042 if( *ncuts >= maxcuts )
1043 goto TERMINATE;
1044 }
1045 else
1046 break;
1047 #if 0
1048 if (nested_cut || disjunct_cut)
1049 for(k = p->beg[p->rcnt - 1]; k < p->nzcnt; k++)
1050 capa[p->ind[k] % nedges
1051 + (((p->ind[k] % nedges) % 2)
1052 ? -1 : 1)] = FLOW_FACTOR;
1053 #endif
1054 }
1055 while( nested_cut ); /* Nested Cut is CONSTANT todo why not only one round? seems to make no sense whatsoever */
1056
1057 rerun = FALSE;
1058 }
1059 }
1060 #endif
1061 TERMINATE:
1062 SCIPfreeBufferArray(scip, &rootcut);
1063 SCIPfreeBufferArray(scip, &start);
1064 SCIPfreeBufferArray(scip, &edgeflipped);
1065 SCIPfreeBufferArray(scip, &headarr);
1066 SCIPfreeBufferArray(scip, &edgearr);
1067
1068 SCIPfreeBufferArray(scip, &term);
1069 SCIPfreeBufferArray(scip, &w);
1070
1071 SCIPfreeBufferArray(scip, &capa);
1072
1073 SCIPdebugMessage("2-cut Separator: %d Inequalities added\n", count);
1074
1075 return SCIP_OKAY;
1076 }
1077
1078
1079 static
dualascent_init(SCIP * scip,const GRAPH * g,const int * RESTRICT start,const int * RESTRICT edgearr,int root,SCIP_Bool is_pseudoroot,int ncsredges,int * gmark,int * RESTRICT active,SCIP_PQUEUE * pqueue,GNODE ** gnodearr,SCIP_Real * RESTRICT rescap,SCIP_Real * dualobj,int * augmentingcomponent)1080 SCIP_RETCODE dualascent_init(
1081 SCIP* scip, /**< SCIP */
1082 const GRAPH* g, /**< graph data structure */
1083 const int* RESTRICT start, /**< CSR start array [0,...,nnodes] */
1084 const int* RESTRICT edgearr, /**< CSR ancestor edge array */
1085 int root, /**< the root */
1086 SCIP_Bool is_pseudoroot, /**< is the root a pseudo root? */
1087 int ncsredges, /**< number of CSR edges */
1088 int* gmark, /**< array for marking nodes */
1089 int* RESTRICT active, /**< active vertices mark */
1090 SCIP_PQUEUE* pqueue, /**< priority queue */
1091 GNODE** gnodearr, /**< array containing terminal nodes*/
1092 SCIP_Real* RESTRICT rescap, /**< residual capacity */
1093 SCIP_Real* dualobj, /**< dual objective */
1094 int* augmentingcomponent /**< augmenting component */
1095 )
1096 {
1097 const int nnodes = g->knots;
1098 *dualobj = 0.0;
1099 *augmentingcomponent = -1;
1100
1101 for( int i = 0; i < ncsredges; i++ )
1102 rescap[i] = g->cost[edgearr[i]];
1103
1104 /* mark terminals as active, add all except root to pqueue */
1105 for( int i = 0, k = 0; i < nnodes; i++ )
1106 {
1107 if( Is_term(g->term[i]) )
1108 {
1109 active[i] = 0;
1110 assert(g->grad[i] > 0);
1111 if( i != root )
1112 {
1113 SCIP_Real warmstart = FALSE;
1114 gnodearr[k]->number = i;
1115 gnodearr[k]->dist = g->grad[i];
1116
1117 /* for variants with dummy terminals */
1118 if( g->grad[i] == 2 )
1119 {
1120 int a;
1121
1122 for( a = g->inpbeg[i]; a != EAT_LAST; a = g->ieat[a] )
1123 if( g->cost[a] == 0.0 )
1124 break;
1125
1126 if( a != EAT_LAST )
1127 {
1128 const int tail = g->tail[a];
1129 gnodearr[k]->dist += g->grad[tail] - 1;
1130
1131 if( is_pseudoroot )
1132 {
1133 SCIP_Bool zeroedge = FALSE;
1134 for( a = g->inpbeg[tail]; a != EAT_LAST; a = g->ieat[a] )
1135 if( g->cost[a] == 0.0 )
1136 {
1137 zeroedge = TRUE;
1138 gnodearr[k]->dist += g->grad[g->tail[a]] - 1;
1139 }
1140
1141 /* warmstart possible? */
1142 if( !zeroedge )
1143 {
1144 int j;
1145 int end;
1146 int prizearc;
1147 SCIP_Real prize;
1148
1149 if( rescap[start[i]] == 0.0 )
1150 prizearc = start[i] + 1;
1151 else
1152 prizearc = start[i];
1153
1154 prize = rescap[prizearc];
1155 assert(prize > 0.0);
1156
1157 for( j = start[tail], end = start[tail + 1]; j != end; j++ )
1158 if( rescap[j] < prize )
1159 break;
1160
1161 if( j == end )
1162 {
1163 warmstart = TRUE;
1164 *dualobj += prize;
1165 rescap[prizearc] = 0.0;
1166 for( j = start[tail], end = start[tail + 1]; j != end; j++ )
1167 rescap[j] -= prize;
1168 }
1169 }
1170 }
1171 }
1172
1173 assert(gnodearr[k]->dist > 0);
1174 }
1175 if( !warmstart )
1176 SCIP_CALL(SCIPpqueueInsert(pqueue, gnodearr[k]));
1177 else if( *augmentingcomponent == -1 )
1178 {
1179 SCIP_CALL(SCIPpqueueInsert(pqueue, gnodearr[k]));
1180 *augmentingcomponent = i;
1181 }
1182 k++;
1183 }
1184 }
1185 else
1186 {
1187 active[i] = -1;
1188 }
1189 }
1190
1191 for( int i = 0; i < nnodes + 1; i++ )
1192 gmark[i] = FALSE;
1193
1194 return SCIP_OKAY;
1195 }
1196
1197 /**@} */
1198
1199
1200 /**@name Callback methods
1201 *
1202 * @{
1203 */
1204
1205 /** copy method for constraint handler plugins (called when SCIP copies plugins) */
1206 static
SCIP_DECL_CONSHDLRCOPY(conshdlrCopyStp)1207 SCIP_DECL_CONSHDLRCOPY(conshdlrCopyStp)
1208 { /*lint --e{715}*/
1209 assert(scip != NULL);
1210 assert(conshdlr != NULL);
1211 assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
1212
1213 /* call inclusion method of constraint handler */
1214 SCIP_CALL( SCIPincludeConshdlrStp(scip) );
1215
1216 *valid = TRUE;
1217
1218 return SCIP_OKAY;
1219 }
1220
1221 /** destructor of constraint handler to free constraint handler data (called when SCIP is exiting) */
1222 static
SCIP_DECL_CONSFREE(consFreeStp)1223 SCIP_DECL_CONSFREE(consFreeStp)
1224 { /*lint --e{715}*/
1225 SCIP_CONSHDLRDATA* conshdlrdata;
1226
1227 assert(scip != NULL);
1228 assert(conshdlr != NULL);
1229 assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
1230
1231 /* free constraint handler data */
1232 conshdlrdata = SCIPconshdlrGetData(conshdlr);
1233 assert(conshdlrdata != NULL);
1234
1235 SCIPfreeMemory(scip, &conshdlrdata);
1236
1237 SCIPconshdlrSetData(conshdlr, NULL);
1238
1239 return SCIP_OKAY;
1240 }
1241
1242 /** solving process initialization method of constraint handler (called when branch and bound process is about to begin) */
1243 static
SCIP_DECL_CONSINITSOL(consInitsolStp)1244 SCIP_DECL_CONSINITSOL(consInitsolStp)
1245 { /*lint --e{715}*/
1246 #ifdef WITH_UG
1247 SCIPStpConshdlrSetGraph(scip, SCIPprobdataGetGraph2(scip));
1248 #endif
1249 return SCIP_OKAY;
1250 }
1251
1252 /** frees specific constraint data */
1253 static
SCIP_DECL_CONSDELETE(consDeleteStp)1254 SCIP_DECL_CONSDELETE(consDeleteStp)
1255 { /*lint --e{715}*/
1256 assert(conshdlr != NULL);
1257 assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
1258 assert(consdata != NULL);
1259 assert(*consdata != NULL);
1260
1261 SCIPfreeBlockMemory(scip, consdata);
1262
1263 return SCIP_OKAY;
1264 }
1265
1266 /** transforms constraint data into data belonging to the transformed problem */
1267 static
SCIP_DECL_CONSTRANS(consTransStp)1268 SCIP_DECL_CONSTRANS(consTransStp)
1269 { /*lint --e{715}*/
1270 SCIP_CONSDATA* sourcedata;
1271 SCIP_CONSDATA* targetdata;
1272
1273 assert(conshdlr != NULL);
1274 assert(strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0);
1275 assert(SCIPgetStage(scip) == SCIP_STAGE_TRANSFORMING);
1276 assert(sourcecons != NULL);
1277 assert(targetcons != NULL);
1278
1279 sourcedata = SCIPconsGetData(sourcecons);
1280 assert(sourcedata != NULL);
1281
1282 /* create constraint data for target constraint */
1283 SCIP_CALL( SCIPallocBlockMemory(scip, &targetdata) );
1284
1285 targetdata->graph = sourcedata->graph;
1286
1287 /* create target constraint */
1288 SCIP_CALL( SCIPcreateCons(scip, targetcons, SCIPconsGetName(sourcecons), conshdlr, targetdata,
1289 SCIPconsIsInitial(sourcecons), SCIPconsIsSeparated(sourcecons), SCIPconsIsEnforced(sourcecons),
1290 SCIPconsIsChecked(sourcecons), SCIPconsIsPropagated(sourcecons),
1291 SCIPconsIsLocal(sourcecons), SCIPconsIsModifiable(sourcecons),
1292 SCIPconsIsDynamic(sourcecons), SCIPconsIsRemovable(sourcecons), SCIPconsIsStickingAtNode(sourcecons)) );
1293
1294 return SCIP_OKAY;
1295 }
1296
1297 #if 1
1298 /** LP initialization method of constraint handler (called before the initial LP relaxation at a node is solved) */
1299 static
SCIP_DECL_CONSINITLP(consInitlpStp)1300 SCIP_DECL_CONSINITLP(consInitlpStp)
1301 { /*lint --e{715}*/
1302 #if 0
1303 SCIP_PROBDATA* probdata;
1304 GRAPH* graph;
1305
1306 SCIP_Real lpobjval;
1307
1308 probdata = SCIPgetProbData(scip);
1309 assert(probdata != NULL);
1310
1311 graph = SCIPprobdataGetGraph(probdata);
1312 assert(graph != NULL);
1313
1314 SCIP_CALL( SCIPdualAscentPcStp(scip, graph, NULL, &lpobjval, TRUE, 1) );
1315 #endif
1316
1317 return SCIP_OKAY;
1318 }
1319 #endif
1320 /** separation method of constraint handler for LP solutions */
1321 static
SCIP_DECL_CONSSEPALP(consSepalpStp)1322 SCIP_DECL_CONSSEPALP(consSepalpStp)
1323 { /*lint --e{715}*/
1324 SCIP_CONSHDLRDATA* conshdlrdata;
1325 SCIP_CONSDATA* consdata;
1326 GRAPH* g;
1327 int* termorg = NULL;
1328 int* nodestatenew = NULL;
1329 int maxcuts;
1330 int ncuts = 0;
1331 const SCIP_Bool atrootnode = (SCIPnodeGetDepth(SCIPgetCurrentNode(scip)) == 0);
1332 #ifndef NDEBUG
1333 int nterms;
1334 #endif
1335
1336 *result = SCIP_DIDNOTRUN;
1337
1338 conshdlrdata = SCIPconshdlrGetData(conshdlr);
1339 assert(conshdlrdata != NULL);
1340
1341 maxcuts = atrootnode ? conshdlrdata->maxsepacutsroot : conshdlrdata->maxsepacuts;
1342
1343 assert(nconss == 1);
1344 consdata = SCIPconsGetData(conss[0]);
1345
1346 assert(consdata != NULL);
1347
1348 g = consdata->graph;
1349 assert(g != NULL);
1350
1351 #ifndef NDEBUG
1352 nterms = g->terms;
1353 #endif
1354
1355 SCIP_CALL( sep_flow(scip, conshdlr, conshdlrdata, consdata, maxcuts, &ncuts) );
1356
1357 /* change graph according to branch-and-bound terminal changes */
1358 if( !atrootnode && g->stp_type == STP_SPG )
1359 {
1360 const int nnodes = g->knots;
1361
1362 SCIP_CALL(SCIPallocBufferArray(scip, &nodestatenew, nnodes));
1363 SCIP_CALL(SCIPallocBufferArray(scip, &termorg, nnodes));
1364 BMScopyMemoryArray(termorg, g->term, nnodes);
1365
1366 SCIPStpBranchruleInitNodeState(g, nodestatenew);
1367 SCIP_CALL( SCIPStpBranchruleApplyVertexChgs(scip, nodestatenew, NULL) );
1368
1369 for( int k = 0; k < nnodes; k++ )
1370 if( nodestatenew[k] == BRANCH_STP_VERTEX_TERM && !Is_term(g->term[k]) )
1371 graph_knot_chg(g, k, 0);
1372 }
1373
1374 SCIP_CALL( sep_2cut(scip, conshdlr, conshdlrdata, consdata, termorg, maxcuts, &ncuts) );
1375
1376 if( ncuts > 0 )
1377 *result = SCIP_SEPARATED;
1378
1379 /* restore graph */
1380 if( !atrootnode && g->stp_type == STP_SPG )
1381 {
1382 for( int k = 0; k < g->knots; k++ )
1383 if( g->term[k] != termorg[k] )
1384 graph_knot_chg(g, k, termorg[k]);
1385 }
1386
1387 #ifndef NDEBUG
1388 assert(g->terms == nterms);
1389 #endif
1390
1391 SCIPfreeBufferArrayNull(scip, &termorg);
1392 SCIPfreeBufferArrayNull(scip, &nodestatenew);
1393
1394 return SCIP_OKAY;
1395 }
1396
1397
1398 /** constraint enforcing method of constraint handler for LP solutions */
1399 static
SCIP_DECL_CONSENFOLP(consEnfolpStp)1400 SCIP_DECL_CONSENFOLP(consEnfolpStp)
1401 { /*lint --e{715}*/
1402 SCIP_Bool feasible;
1403 SCIP_CONSDATA* consdata;
1404 int i;
1405
1406 for( i = 0; i < nconss; i++ )
1407 {
1408 consdata = SCIPconsGetData(conss[i]);
1409
1410 SCIP_CALL( SCIPStpValidateSol(scip, consdata->graph, SCIPprobdataGetXval(scip, NULL), &feasible) );
1411
1412 if( !feasible )
1413 {
1414 *result = SCIP_INFEASIBLE;
1415 return SCIP_OKAY;
1416 }
1417 }
1418 *result = SCIP_FEASIBLE;
1419
1420 return SCIP_OKAY;
1421 }
1422
1423 /** constraint enforcing method of constraint handler for pseudo solutions */
1424 static
SCIP_DECL_CONSENFOPS(consEnfopsStp)1425 SCIP_DECL_CONSENFOPS(consEnfopsStp)
1426 { /*lint --e{715}*/
1427 SCIP_Bool feasible;
1428
1429 assert(nconss == 1);
1430
1431 for( int i = 0; i < nconss; i++ )
1432 {
1433 const SCIP_CONSDATA* consdata = SCIPconsGetData(conss[i]);
1434
1435 SCIP_CALL( SCIPStpValidateSol(scip, consdata->graph, SCIPprobdataGetXval(scip, NULL), &feasible) );
1436
1437 if( !feasible )
1438 {
1439 *result = SCIP_INFEASIBLE;
1440 return SCIP_OKAY;
1441 }
1442 }
1443 *result = SCIP_FEASIBLE;
1444
1445 return SCIP_OKAY;
1446 }
1447
1448 /** feasibility check method of constraint handler for integral solutions */
1449 static
SCIP_DECL_CONSCHECK(consCheckStp)1450 SCIP_DECL_CONSCHECK(consCheckStp)
1451 { /*lint --e{715}*/
1452 const GRAPH* g = SCIPprobdataGetGraph2(scip);
1453 SCIP_Bool feasible;
1454
1455 assert(g != NULL);
1456
1457 SCIP_CALL(SCIPStpValidateSol(scip, g, SCIPprobdataGetXval(scip, sol), &feasible));
1458
1459 if( !feasible )
1460 {
1461 *result = SCIP_INFEASIBLE;
1462 return SCIP_OKAY;
1463 }
1464
1465 *result = SCIP_FEASIBLE;
1466
1467 return SCIP_OKAY;
1468 }
1469
1470 /** domain propagation method of constraint handler */
1471 static
SCIP_DECL_CONSPROP(consPropStp)1472 SCIP_DECL_CONSPROP(consPropStp)
1473 { /*lint --e{715}*/
1474 SCIP_PROBDATA* probdata;
1475 GRAPH* graph;
1476
1477 probdata = SCIPgetProbData(scip);
1478 assert(probdata != NULL);
1479
1480 graph = SCIPprobdataGetGraph(probdata);
1481 assert(graph != NULL);
1482
1483 /* for degree constrained model, check whether problem is infeasible */
1484 if( graph->stp_type == STP_DCSTP )
1485 {
1486 int k;
1487 int nnodes;
1488 int degsum;
1489 int* maxdegs;
1490
1491 nnodes = graph->knots;
1492 maxdegs = graph->maxdeg;
1493
1494 assert(maxdegs != NULL);
1495
1496 degsum = 0;
1497 for( k = 0; k < nnodes; k++ )
1498 {
1499 if( Is_term(graph->term[k]) )
1500 {
1501 assert(maxdegs[k] > 0);
1502 degsum += maxdegs[k] - 1;
1503 }
1504 else
1505 {
1506 assert(maxdegs[k] >= 0);
1507 degsum += MAX(maxdegs[k] - 2, 0);
1508 }
1509 }
1510
1511 if( degsum < graph->terms - 2 )
1512 *result = SCIP_CUTOFF;
1513 else
1514 *result = SCIP_DIDNOTFIND;
1515 }
1516 return SCIP_OKAY;
1517 }
1518
1519 /** variable rounding lock method of constraint handler */
1520 static
SCIP_DECL_CONSLOCK(consLockStp)1521 SCIP_DECL_CONSLOCK(consLockStp)
1522 { /*lint --e{715}*/
1523 SCIP_VAR** vars;
1524 int nvars;
1525 int v;
1526
1527 assert(scip != NULL);
1528 assert(cons != NULL);
1529
1530 vars = SCIPprobdataGetVars(scip);
1531 nvars = SCIPprobdataGetNVars(scip);
1532
1533 for( v = 0; v < nvars; ++v )
1534 {
1535 SCIP_CALL( SCIPaddVarLocksType(scip, vars[v], locktype, nlockspos, nlocksneg) );
1536 }
1537
1538 return SCIP_OKAY;
1539 }
1540
1541 /** constraint copying method of constraint handler */
1542 static
SCIP_DECL_CONSCOPY(consCopyStp)1543 SCIP_DECL_CONSCOPY(consCopyStp)
1544 { /*lint --e{715}*/
1545 const char* consname;
1546 SCIP_PROBDATA* probdata;
1547 GRAPH* graph;
1548
1549 probdata = SCIPgetProbData(scip);
1550 assert(probdata != NULL);
1551
1552 graph = SCIPprobdataGetGraph(probdata);
1553 assert(graph != NULL);
1554
1555 consname = SCIPconsGetName(sourcecons);
1556
1557 /* creates and captures a and constraint */
1558 SCIP_CALL( SCIPcreateConsStp(scip, cons, consname, graph) );
1559
1560 *valid = TRUE;
1561
1562 return SCIP_OKAY;
1563 }
1564
1565
1566 /**@} */
1567
1568 /**@name Interface methods
1569 *
1570 * @{
1571 */
1572
1573 /** creates the handler for stp constraints and includes it in SCIP */
SCIPincludeConshdlrStp(SCIP * scip)1574 SCIP_RETCODE SCIPincludeConshdlrStp(
1575 SCIP* scip /**< SCIP data structure */
1576 )
1577 {
1578 SCIP_CONSHDLRDATA* conshdlrdata;
1579 SCIP_CONSHDLR* conshdlr;
1580
1581 /* create stp constraint handler data */
1582 SCIP_CALL( SCIPallocMemory(scip, &conshdlrdata) );
1583
1584 conshdlr = NULL;
1585 /* include constraint handler */
1586 SCIP_CALL( SCIPincludeConshdlrBasic(scip, &conshdlr, CONSHDLR_NAME, CONSHDLR_DESC,
1587 CONSHDLR_ENFOPRIORITY, CONSHDLR_CHECKPRIORITY, CONSHDLR_EAGERFREQ, CONSHDLR_NEEDSCONS,
1588 consEnfolpStp, consEnfopsStp, consCheckStp, consLockStp,
1589 conshdlrdata) );
1590 assert(conshdlr != NULL);
1591
1592 SCIP_CALL( SCIPsetConshdlrCopy(scip, conshdlr, conshdlrCopyStp, consCopyStp) );
1593 SCIP_CALL( SCIPsetConshdlrDelete(scip, conshdlr, consDeleteStp) );
1594 SCIP_CALL( SCIPsetConshdlrTrans(scip, conshdlr, consTransStp) );
1595 SCIP_CALL( SCIPsetConshdlrInitsol(scip, conshdlr, consInitsolStp) );
1596 SCIP_CALL( SCIPsetConshdlrInitlp(scip, conshdlr, consInitlpStp) );
1597 SCIP_CALL( SCIPsetConshdlrProp(scip, conshdlr, consPropStp, CONSHDLR_PROPFREQ, CONSHDLR_DELAYPROP,
1598 CONSHDLR_PROP_TIMING) );
1599 SCIP_CALL( SCIPsetConshdlrSepa(scip, conshdlr, consSepalpStp, NULL, CONSHDLR_SEPAFREQ,
1600 CONSHDLR_SEPAPRIORITY, CONSHDLR_DELAYSEPA) );
1601 SCIP_CALL( SCIPsetConshdlrFree(scip, conshdlr, consFreeStp) );
1602
1603 SCIP_CALL( SCIPaddBoolParam(scip, "constraints/stp/backcut", "Try Back-Cuts",
1604 &conshdlrdata->backcut, TRUE, DEFAULT_BACKCUT, NULL, NULL) );
1605 SCIP_CALL( SCIPaddBoolParam(scip, "constraints/stp/creepflow", "Use Creep-Flow",
1606 &conshdlrdata->creepflow, TRUE, DEFAULT_CREEPFLOW, NULL, NULL) );
1607 SCIP_CALL( SCIPaddBoolParam(scip, "constraints/stp/disjunctcut", "Only disjunct Cuts",
1608 &conshdlrdata->disjunctcut, TRUE, DEFAULT_DISJUNCTCUT, NULL, NULL) );
1609 SCIP_CALL( SCIPaddBoolParam(scip, "constraints/stp/nestedcut", "Try Nested-Cuts",
1610 &conshdlrdata->nestedcut, TRUE, DEFAULT_NESTEDCUT, NULL, NULL) );
1611 SCIP_CALL( SCIPaddBoolParam(scip, "constraints/stp/flowsep", "Try Flow-Cuts",
1612 &conshdlrdata->flowsep, TRUE, DEFAULT_FLOWSEP, NULL, NULL) );
1613 SCIP_CALL( SCIPaddIntParam(scip,
1614 "constraints/"CONSHDLR_NAME"/maxrounds",
1615 "maximal number of separation rounds per node (-1: unlimited)",
1616 &conshdlrdata->maxrounds, FALSE, DEFAULT_MAXROUNDS, -1, INT_MAX, NULL, NULL) );
1617 SCIP_CALL( SCIPaddIntParam(scip,
1618 "constraints/"CONSHDLR_NAME"/maxroundsroot",
1619 "maximal number of separation rounds per node in the root node (-1: unlimited)",
1620 &conshdlrdata->maxroundsroot, FALSE, DEFAULT_MAXROUNDSROOT, -1, INT_MAX, NULL, NULL) );
1621 SCIP_CALL( SCIPaddIntParam(scip,
1622 "constraints/"CONSHDLR_NAME"/maxsepacuts",
1623 "maximal number of cuts separated per separation round",
1624 &conshdlrdata->maxsepacuts, FALSE, DEFAULT_MAXSEPACUTS, 0, INT_MAX, NULL, NULL) );
1625 SCIP_CALL( SCIPaddIntParam(scip,
1626 "constraints/"CONSHDLR_NAME"/maxsepacutsroot",
1627 "maximal number of cuts separated per separation round in the root node",
1628 &conshdlrdata->maxsepacutsroot, FALSE, DEFAULT_MAXSEPACUTSROOT, 0, INT_MAX, NULL, NULL) );
1629
1630
1631 return SCIP_OKAY;
1632 }
1633
1634 /** creates and captures a stp constraint */
SCIPcreateConsStp(SCIP * scip,SCIP_CONS ** cons,const char * name,GRAPH * graph)1635 SCIP_RETCODE SCIPcreateConsStp(
1636 SCIP* scip, /**< SCIP data structure */
1637 SCIP_CONS** cons, /**< pointer to hold the created constraint */
1638 const char* name, /**< name of constraint */
1639 GRAPH* graph /**< graph data structure */
1640 )
1641 {
1642 SCIP_CONSHDLR* conshdlr;
1643 SCIP_CONSDATA* consdata;
1644
1645 /* find the stp constraint handler */
1646 conshdlr = SCIPfindConshdlr(scip, CONSHDLR_NAME);
1647 if( conshdlr == NULL )
1648 {
1649 SCIPerrorMessage("stp constraint handler not found\n");
1650 return SCIP_PLUGINNOTFOUND;
1651 }
1652
1653 SCIP_CALL( SCIPallocBlockMemory(scip, &consdata) );
1654
1655 consdata->graph = graph;
1656
1657 /* create constraint */
1658 SCIP_CALL( SCIPcreateCons(scip, cons, name, conshdlr, consdata, FALSE, TRUE, TRUE, TRUE, TRUE,
1659 FALSE, FALSE, FALSE, FALSE, FALSE) );
1660
1661 return SCIP_OKAY;
1662 }
1663
1664 /** sets graph */
SCIPStpConshdlrSetGraph(SCIP * scip,const GRAPH * g)1665 void SCIPStpConshdlrSetGraph(
1666 SCIP* scip, /**< SCIP data structure */
1667 const GRAPH* g /**< graph data structure */
1668 )
1669 {
1670 SCIP_CONSDATA* consdata;
1671 SCIP_CONSHDLR* conshdlr;
1672
1673 conshdlr = SCIPfindConshdlr(scip, "stp");
1674 assert(conshdlr != NULL);
1675 assert(SCIPconshdlrGetNConss(conshdlr) > 0);
1676
1677 consdata = SCIPconsGetData(SCIPconshdlrGetConss(conshdlr)[0]);
1678
1679 assert(consdata != NULL);
1680
1681 consdata->graph = SCIPprobdataGetGraph2(scip);
1682 assert(consdata->graph != NULL);
1683 }
1684
1685 /* dual ascent heuristic */
SCIPStpDualAscent(SCIP * scip,const GRAPH * g,SCIP_Real * RESTRICT redcost,SCIP_Real * RESTRICT nodearrreal,SCIP_Real * objval,SCIP_Bool addcuts,SCIP_Bool ascendandprune,GNODE ** gnodearrterms,const int * result,int * RESTRICT edgearrint,int * RESTRICT nodearrint,int root,SCIP_Bool is_pseudoroot,SCIP_Real damaxdeviation,STP_Bool * RESTRICT nodearrchar)1686 SCIP_RETCODE SCIPStpDualAscent(
1687 SCIP* scip, /**< SCIP data structure */
1688 const GRAPH* g, /**< graph data structure */
1689 SCIP_Real* RESTRICT redcost, /**< array to store reduced costs or NULL */
1690 SCIP_Real* RESTRICT nodearrreal, /**< real vertices array for internal computations or NULL */
1691 SCIP_Real* objval, /**< pointer to store objective value */
1692 SCIP_Bool addcuts, /**< should dual ascent add Steiner cuts? */
1693 SCIP_Bool ascendandprune, /**< should the ascent-and-prune heuristic be executed? */
1694 GNODE** gnodearrterms, /**< gnode terminals array for internal computations or NULL */
1695 const int* result, /**< solution array or NULL */
1696 int* RESTRICT edgearrint, /**< int edges array for internal computations or NULL */
1697 int* RESTRICT nodearrint, /**< int vertices array for internal computations or NULL */
1698 int root, /**< the root */
1699 SCIP_Bool is_pseudoroot, /**< is the root a pseudo root? */
1700 SCIP_Real damaxdeviation, /**< maximum deviation for dual-ascent ( -1.0 for default) */
1701 STP_Bool* RESTRICT nodearrchar /**< STP_Bool vertices array for internal computations or NULL */
1702 )
1703 {
1704 SCIP_CONSHDLR* conshdlr = NULL;
1705 SCIP_PQUEUE* pqueue;
1706 SCIP_VAR** vars;
1707 SCIP_Real* RESTRICT rescap;
1708 GNODE** gnodearr;
1709 int* RESTRICT edgearr;
1710 int* RESTRICT tailarr;
1711 int* RESTRICT start;
1712 int* RESTRICT stackarr;
1713 int* RESTRICT cutverts;
1714 int* RESTRICT unsatarcs;
1715 int* RESTRICT unsattails;
1716 int* RESTRICT gmark;
1717 int* RESTRICT active;
1718 SCIP_Real dualobj;
1719 SCIP_Real currscore;
1720 const SCIP_Real maxdeviation = (damaxdeviation > 0.0) ? damaxdeviation : DEFAULT_DAMAXDEVIATION;
1721 const int nnodes = g->knots;
1722 const int nterms = g->terms;
1723 const int nedges = g->edges;
1724 int ncsredges;
1725 int norgcutverts;
1726 int stacklength;
1727 int augmentingcomponent;
1728 const SCIP_Bool addconss = (SCIPgetStage(scip) < SCIP_STAGE_INITSOLVE);
1729
1730 /* should currently not be activated */
1731 assert(addconss || !addcuts);
1732 assert(g != NULL);
1733 assert(scip != NULL);
1734 assert(objval != NULL);
1735 assert(Is_term(g->term[root]));
1736 assert(maxdeviation >= DA_MAXDEVIATION_LOWER && maxdeviation <= DA_MAXDEVIATION_UPPER);
1737 assert(damaxdeviation == -1.0 || damaxdeviation > 0.0);
1738
1739 if( nnodes == 1 )
1740 return SCIP_OKAY;
1741
1742 if( addcuts )
1743 {
1744 vars = SCIPprobdataGetVars(scip);
1745 assert(vars != NULL);
1746
1747 if( !addconss )
1748 {
1749 conshdlr = SCIPfindConshdlr(scip, "stp");
1750 assert(conshdlr != NULL);
1751 }
1752 }
1753 else
1754 {
1755 vars = NULL;
1756 }
1757
1758 /* if specified root is not a terminal, take default root */
1759 if( !Is_term(g->term[root]) )
1760 root = g->source;
1761
1762 #ifdef BITFIELDSARRAY
1763 u_int32_t* bitarr;
1764 SCIP_CALL( SCIPallocBufferArray(scip, &bitarr, nedges / ARRLENGTH + 1) );
1765 #endif
1766
1767 stacklength = 0;
1768
1769 SCIP_CALL( SCIPallocBufferArray(scip, &unsattails, nedges) );
1770
1771 if( redcost == NULL )
1772 SCIP_CALL( SCIPallocBufferArray(scip, &rescap, nedges) );
1773 else
1774 rescap = redcost;
1775
1776 if( nodearrint == NULL )
1777 SCIP_CALL( SCIPallocBufferArray(scip, &cutverts, nnodes) );
1778 else
1779 cutverts = nodearrint;
1780
1781 if( edgearrint == NULL )
1782 SCIP_CALL( SCIPallocBufferArray(scip, &unsatarcs, nedges) );
1783 else
1784 unsatarcs = edgearrint;
1785
1786 if( gnodearrterms == NULL )
1787 {
1788 SCIP_CALL( SCIPallocBufferArray(scip, &gnodearr, nterms - 1) );
1789 for( int i = 0; i < nterms - 1; i++ )
1790 SCIP_CALL( SCIPallocBlockMemory(scip, &gnodearr[i]) ); /*lint !e866*/
1791 }
1792 else
1793 {
1794 gnodearr = gnodearrterms;
1795 }
1796
1797 SCIP_CALL( SCIPpqueueCreate(&pqueue, nterms, 2.0, GNODECmpByDist, NULL) );
1798
1799 SCIP_CALL( SCIPallocMemoryArray(scip, &active, nnodes) );
1800 SCIP_CALL( SCIPallocMemoryArray(scip, &edgearr, nedges) );
1801 SCIP_CALL( SCIPallocMemoryArray(scip, &tailarr, nedges) );
1802 SCIP_CALL( SCIPallocMemoryArray(scip, &start, nnodes + 1) );
1803 SCIP_CALL( SCIPallocMemoryArray(scip, &gmark, nnodes + 1) );
1804 SCIP_CALL( SCIPallocMemoryArray(scip, &stackarr, nnodes) );
1805
1806 /* fill auxiliary adjacent vertex/edges arrays */
1807 graph_get_csr(g, edgearr, tailarr, start, &ncsredges);
1808
1809 /* initialize priority queue and res. capacity */
1810 SCIP_CALL( dualascent_init(scip, g, start, edgearr, root, is_pseudoroot, ncsredges, gmark, active, pqueue,
1811 gnodearr, rescap, &dualobj, &augmentingcomponent) );
1812
1813 /* mark whether an arc is satisfied (has capacity 0) */
1814 for( int i = 0; i < ncsredges; i++ )
1815 {
1816 #ifdef BITFIELDSARRAY
1817 if( SCIPisZero(scip, rescap[i]) )
1818 SetBit(bitarr, i);
1819 else
1820 CleanBit(bitarr, i);
1821 #else
1822 if( rescap[i] == 0.0 )
1823 {
1824 if( active[tailarr[i] - 1] == 0 )
1825 tailarr[i] = 0;
1826 else
1827 tailarr[i] *= -1;
1828 }
1829 #endif
1830 }
1831
1832 norgcutverts = 0;
1833
1834 /* (main) dual ascent loop */
1835 while( SCIPpqueueNElems(pqueue) > 0 && !SCIPisStopped(scip) )
1836 {
1837 /* get active vertex of minimum score */
1838 GNODE* const gnodeact = (GNODE*) SCIPpqueueRemove(pqueue);
1839 const SCIP_Real prio1 = gnodeact->dist;
1840 const SCIP_Real prio2 = (SCIPpqueueNElems(pqueue) > 0) ? ((GNODE*) SCIPpqueueFirst(pqueue))->dist : FARAWAY;
1841 const int v = gnodeact->number;
1842 SCIP_Real degsum = g->grad[v];
1843 int ncutverts = 0;
1844 int nunsatarcs = 0;
1845
1846 SCIP_Bool firstrun = TRUE;
1847
1848 SCIPdebugMessage("DA: START WITH v %d prio1 %f prio2 %f \n", v, prio1, prio2);
1849
1850 /* perform augmentation as long as priority of root component does not exceed max deviation */
1851 for( ; ; )
1852 {
1853 assert(stacklength == 0);
1854
1855 /* 1. step: BFS from v (or connected component) on saturated, incoming arcs */
1856
1857 if( firstrun )
1858 {
1859 firstrun = FALSE;
1860 gmark[v + 1] = TRUE;
1861 cutverts[ncutverts++] = v;
1862 assert(stacklength < nnodes);
1863 stackarr[stacklength++] = v;
1864 }
1865 /* not in first processing of root component: */
1866 else
1867 {
1868 for( int i = norgcutverts; i < ncutverts; i++ )
1869 {
1870 const int s = cutverts[i];
1871
1872 assert(gmark[s + 1]);
1873 assert(active[s] != 0);
1874 assert(stacklength < nnodes);
1875
1876 stackarr[stacklength++] = s;
1877 }
1878 }
1879 #ifdef DFS
1880 while( stacklength )
1881 {
1882 const int node = stackarr[--stacklength];
1883 #else
1884 for( int n = 0; n < stacklength; n++ )
1885 {
1886 int end;
1887
1888 assert(n < nnodes);
1889 node = stackarr[n];
1890 #endif
1891
1892 /* traverse incoming arcs */
1893 for( int i = start[node], end = start[node + 1]; i != end; i++ )
1894 {
1895 int tail = tailarr[i];
1896
1897 /* zero reduced-cost arc? */
1898 if( tail <= 0 )
1899 {
1900 tail *= -1;
1901 if( !gmark[tail] )
1902 {
1903 /* if an active vertex has been hit (other than v), break */
1904 if( 0 == tail )
1905 {
1906 const int realtail = g->tail[edgearr[i]];
1907
1908 /* v should not be processed */
1909 if( realtail == v )
1910 continue;
1911
1912 /* is realtail active or does realtail lead to an active vertex other than v? */
1913 if( is_active(active, realtail, v) )
1914 {
1915 active[v] = realtail + 1;
1916 stacklength = 0;
1917 goto ENDOFLOOP;
1918 }
1919
1920 tail = realtail + 1;
1921
1922 /* have we processed tail already? */
1923 if( gmark[tail] )
1924 continue;
1925 }
1926
1927 assert(tail > 0);
1928
1929 gmark[tail] = TRUE;
1930 tail--;
1931 cutverts[ncutverts++] = tail;
1932 degsum += g->grad[tail];
1933
1934 assert(stacklength < nnodes);
1935 stackarr[stacklength++] = tail;
1936 } /* marked */
1937 } /* zero reduced-cost arc */
1938 else if( !gmark[tail] )
1939 {
1940 unsattails[nunsatarcs] = tail;
1941 unsatarcs[nunsatarcs++] = i;
1942 }
1943 }
1944 }
1945 #ifndef DFS
1946 stacklength = 0;
1947 #endif
1948 currscore = degsum - (ncutverts - 1);
1949
1950 /* guiding solution provided? */
1951 if( result != NULL )
1952 {
1953 int nsolarcs = 0;
1954 for( int i = 0; i < nunsatarcs; i++ )
1955 {
1956 const int a = unsatarcs[i];
1957
1958 assert(tailarr[a] > 0);
1959
1960 if( !(gmark[tailarr[a]]) )
1961 {
1962 if( result[edgearr[a]] == CONNECT )
1963 nsolarcs++;
1964 }
1965 }
1966
1967 assert(nsolarcs > 0);
1968 assert(currscore <= nedges);
1969
1970 if( nsolarcs > 1 )
1971 currscore += (SCIP_Real) ((nsolarcs - 1) * (g->knots * 2.0));
1972 }
1973 else
1974 {
1975 assert(SCIPisGE(scip, currscore, prio1));
1976 }
1977
1978 SCIPdebugMessage("DA: deviation %f \n", (currscore - prio1) / prio1);
1979 SCIPdebugMessage("DA: currscore %f prio1 %f prio2 %f \n", currscore, prio1, prio2);
1980
1981 /* augmentation criteria met? */
1982 if( ((currscore - prio1) / prio1) <= maxdeviation || currscore <= prio2 )
1983 {
1984 SCIP_CONS* cons = NULL;
1985 SCIP_ROW* row = NULL;
1986
1987 int shift = 0;
1988 SCIP_Real min = FARAWAY;
1989 SCIP_Bool isactive = FALSE;
1990
1991 /* 2. step: get minimum residual capacity among cut-arcs */
1992
1993 /* adjust array of unsatisfied arcs */
1994
1995 for( int i = 0; i < nunsatarcs; i++ )
1996 {
1997 const int tail = unsattails[i];
1998
1999 if( gmark[tail] )
2000 {
2001 shift++;
2002 }
2003 else
2004 {
2005 const int a = unsatarcs[i];
2006
2007 assert(tailarr[a] > 0);
2008 assert(rescap[a] > 0);
2009
2010 if( rescap[a] < min )
2011 min = rescap[a];
2012 if( shift )
2013 {
2014 unsattails[i - shift] = tail;
2015 unsatarcs[i - shift] = a;
2016 }
2017 }
2018 }
2019
2020 assert(SCIPisLT(scip, min, FARAWAY));
2021 nunsatarcs -= shift;
2022
2023 norgcutverts = ncutverts;
2024
2025 /* 3. step: perform augmentation */
2026
2027 /* create constraints/cuts ? */
2028 if( addcuts )
2029 {
2030 if( addconss )
2031 {
2032 SCIP_CALL( SCIPcreateConsLinear(scip, &cons, "da", 0, NULL, NULL,
2033 1.0, SCIPinfinity(scip), TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, FALSE) );
2034 }
2035 else
2036 {
2037 SCIP_CALL( SCIPcreateEmptyRowConshdlr(scip, &row, conshdlr, "da", 1.0,
2038 SCIPinfinity(scip), FALSE, FALSE, TRUE) );
2039
2040 SCIP_CALL( SCIPcacheRowExtensions(scip, row) );
2041 }
2042 }
2043
2044 shift = 0;
2045
2046 /* update (dual) objective */
2047 dualobj += min;
2048
2049 for( int i = 0; i < nunsatarcs; i++ )
2050 {
2051 const int a = unsatarcs[i];
2052 assert(a >= 0);
2053
2054 if( addcuts )
2055 {
2056 assert(vars != NULL);
2057
2058 if( addconss )
2059 SCIP_CALL( SCIPaddCoefLinear(scip, cons, vars[edgearr[a]], 1.0) );
2060 else
2061 SCIP_CALL( SCIPaddVarToRow(scip, row, vars[edgearr[a]], 1.0) );
2062 }
2063 rescap[a] -= min;
2064
2065 assert(SCIPisGE(scip, rescap[a], 0.0));
2066
2067 if( rescap[a] <= DA_EPS )
2068 {
2069 int tail = unsattails[i];
2070
2071 rescap[a] = 0.0;
2072
2073 assert(tail > 0);
2074 assert(tailarr[a] > 0);
2075
2076 tailarr[a] *= -1;
2077
2078 if( active[tail - 1] >= 0 && is_active(active, tail - 1, v) )
2079 {
2080 assert(tail - 1 != v);
2081 tailarr[a] = 0;
2082 if( !isactive )
2083 {
2084 isactive = TRUE;
2085 active[v] = tail;
2086 }
2087 }
2088
2089
2090 if( !(gmark[tail]) )
2091 {
2092 assert(tail != 0);
2093
2094 gmark[tail] = TRUE;
2095 tail--;
2096 degsum += g->grad[tail];
2097 cutverts[ncutverts++] = tail;
2098 }
2099
2100 shift++;
2101 }
2102 else if( shift )
2103 {
2104 unsattails[i - shift] = unsattails[i];
2105 unsatarcs[i - shift] = a;
2106 }
2107 }
2108
2109 if( addcuts )
2110 {
2111 if( addconss )
2112 {
2113 SCIP_CALL( SCIPaddCons(scip, cons) );
2114 SCIP_CALL( SCIPreleaseCons(scip, &cons) );
2115 }
2116 else
2117 {
2118 SCIP_Bool infeasible;
2119
2120 SCIP_CALL( SCIPflushRowExtensions(scip, row) );
2121 SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
2122 SCIP_CALL( SCIPreleaseRow(scip, &row) );
2123
2124 assert(!infeasible);
2125 }
2126 }
2127
2128 if( isactive )
2129 {
2130 stacklength = 0;
2131 goto ENDOFLOOP;
2132 }
2133 nunsatarcs -= shift;
2134
2135 continue;
2136 }
2137 else
2138 {
2139 SCIP_Bool insert = TRUE;
2140
2141 if( is_pseudoroot )
2142 {
2143 int i = start[v];
2144 const int end = start[v + 1];
2145
2146 assert(end - i == 2);
2147
2148 for( ; i != end; i++ )
2149 if( rescap[i] != 0.0 )
2150 break;
2151
2152 if( i == end )
2153 {
2154 if( augmentingcomponent == -1 )
2155 augmentingcomponent = v;
2156
2157 if( augmentingcomponent != v )
2158 insert = FALSE;
2159 }
2160 }
2161
2162 if( insert )
2163 {
2164 /* reinsert active vertex */
2165 gnodeact->dist = currscore;
2166 SCIP_CALL( SCIPpqueueInsert(pqueue, gnodeact) );
2167 }
2168 }
2169
2170 ENDOFLOOP:
2171
2172 for( int i = 0; i < ncutverts; i++ )
2173 gmark[cutverts[i] + 1] = FALSE;
2174
2175 for( int i = 0; i < nnodes + 1; i++ )
2176 {
2177 assert(!gmark[i]);
2178 }
2179
2180 break;
2181 } /* augmentation loop */
2182 } /* dual ascent loop */
2183
2184 SCIPdebugMessage("DA: dualglobal: %f \n", dualobj);
2185 *objval = dualobj;
2186
2187 for( int i = ncsredges; i < nedges; i++ )
2188 {
2189 edgearr[i] = i;
2190 rescap[i] = g->cost[i];
2191 }
2192
2193 /* re-extend rescap array */
2194 for( int i = 0; i < ncsredges; i++ )
2195 {
2196 if( edgearr[i] != i )
2197 {
2198 SCIP_Real bufferedval = rescap[i];
2199 int a = i;
2200
2201 rescap[i] = g->cost[i];
2202 while( edgearr[a] != a )
2203 {
2204 const int shift = edgearr[a];
2205 const SCIP_Real min = rescap[shift];
2206
2207 rescap[shift] = bufferedval;
2208 bufferedval = min;
2209 edgearr[a] = a;
2210 a = shift;
2211 }
2212 }
2213 }
2214
2215 #ifdef BITFIELDSARRAY
2216 SCIPfreeBufferArray(scip, &bitarr);
2217 #endif
2218
2219 SCIPfreeMemoryArray(scip, &stackarr);
2220 SCIPfreeMemoryArray(scip, &gmark);
2221 SCIPfreeMemoryArray(scip, &start);
2222 SCIPfreeMemoryArray(scip, &tailarr);
2223 SCIPfreeMemoryArray(scip, &edgearr);
2224 SCIPfreeMemoryArray(scip, &active);
2225
2226 SCIPpqueueFree(&pqueue);
2227
2228 if( gnodearrterms == NULL )
2229 {
2230 for( int i = nterms - 2; i >= 0; i-- )
2231 SCIPfreeBlockMemory(scip, &gnodearr[i]);
2232 SCIPfreeBufferArray(scip, &gnodearr);
2233 }
2234
2235 /* call Ascend-And-Prune? */
2236 if( ascendandprune )
2237 {
2238 SCIP_Bool success;
2239 STP_Bool* RESTRICT mynodearrchar = NULL;
2240
2241 if( nodearrchar == NULL )
2242 SCIP_CALL( SCIPallocBufferArray(scip, &mynodearrchar, nnodes) );
2243 else
2244 mynodearrchar = nodearrchar;
2245
2246 SCIP_CALL( SCIPStpHeurAscendPruneRun(scip, NULL, g, rescap, unsatarcs, cutverts, root, mynodearrchar, &success, TRUE) );
2247
2248 if( nodearrchar == NULL )
2249 SCIPfreeBufferArray(scip, &mynodearrchar);
2250 }
2251
2252 if( edgearrint == NULL )
2253 SCIPfreeBufferArray(scip, &unsatarcs);
2254
2255 if( nodearrint == NULL )
2256 SCIPfreeBufferArray(scip, &cutverts);
2257
2258 if( redcost == NULL )
2259 SCIPfreeBufferArray(scip, &rescap);
2260
2261 SCIPfreeBufferArray(scip, &unsattails);
2262
2263 return SCIP_OKAY;
2264 }
2265
2266 /** dual ascent heuristic for PCSPG and MWCSP */
2267 SCIP_RETCODE SCIPStpDualAscentPcMw(
2268 SCIP* scip, /**< SCIP data structure */
2269 GRAPH* g, /**< graph data structure */
2270 SCIP_Real* redcost, /**< array to store reduced costs or NULL */
2271 SCIP_Real* objval, /**< pointer to store objective value */
2272 SCIP_Bool addcuts, /**< should dual-ascent add Steiner cuts? */
2273 SCIP_Bool ascendandprune, /**< perform ascend-and-prune and add solution? */
2274 int nruns /**< number of dual ascent runs */
2275 )
2276 {
2277 SCIP_CONSHDLR* conshdlr = NULL;
2278 SCIP_PQUEUE* pqueue;
2279 SCIP_VAR** vars;
2280 GRAPH* transgraph;
2281 SCIP_Real min;
2282 SCIP_Real prio1;
2283 SCIP_Real offset;
2284 SCIP_Real dualobj;
2285 SCIP_Real currscore;
2286 SCIP_Real maxdeviation;
2287 SCIP_Real* rescap;
2288 GNODE* gnodeact;
2289 GNODE** gnodearr;
2290 int s;
2291 int i;
2292 int k;
2293 int v;
2294 int a;
2295 int tail;
2296 int pnode;
2297 int shift;
2298 int root;
2299 int nnodes;
2300 int nterms;
2301 int nedges;
2302 int degsum;
2303 int ncutverts;
2304 int pseudoroot;
2305 int nunsatarcs;
2306 int stacklength;
2307 int norgcutverts;
2308 int* cutverts;
2309 int* stackarr;
2310 STP_Bool* origedge;
2311 int* unsatarcs;
2312 STP_Bool firstrun;
2313 STP_Bool* sat;
2314 STP_Bool* active;
2315 const SCIP_Bool addconss = (SCIPgetStage(scip) < SCIP_STAGE_INITSOLVE);
2316
2317 /* should currently not be activated */
2318 assert(addconss || !addcuts);
2319
2320 assert(g != NULL);
2321 assert(scip != NULL);
2322 assert(nruns >= 0);
2323 assert(objval != NULL);
2324
2325 if( g->knots == 1 )
2326 return SCIP_OKAY;
2327
2328 if( addcuts )
2329 {
2330 vars = SCIPprobdataGetVars(scip);
2331 assert(vars != NULL);
2332 if( !addconss )
2333 {
2334 conshdlr = SCIPfindConshdlr(scip, "stp");
2335 assert(conshdlr != NULL);
2336 }
2337 }
2338 else
2339 {
2340 vars = NULL;
2341 }
2342
2343 root = g->source;
2344 degsum = 0;
2345 offset = 0.0;
2346 dualobj = 0.0;
2347
2348 ncutverts = 0;
2349 norgcutverts = 0;
2350 maxdeviation = DEFAULT_DAMAXDEVIATION;
2351
2352 SCIP_CALL( graph_pc_getSap(scip, g, &transgraph, &offset) );
2353
2354 nnodes = transgraph->knots;
2355 nedges = transgraph->edges;
2356 nterms = transgraph->terms;
2357 pseudoroot = nnodes - 1;
2358
2359 if( redcost == NULL )
2360 {
2361 SCIP_CALL( SCIPallocBufferArray(scip, &rescap, nedges) );
2362 }
2363 else
2364 {
2365 rescap = redcost;
2366 }
2367
2368 stacklength = 0;
2369 SCIP_CALL( SCIPallocBufferArray(scip, &stackarr, nnodes) );
2370 SCIP_CALL( SCIPallocBufferArray(scip, &sat, nedges) );
2371 SCIP_CALL( SCIPallocBufferArray(scip, &active, nnodes) );
2372 SCIP_CALL( SCIPallocBufferArray(scip, &cutverts, nnodes) );
2373 SCIP_CALL( SCIPallocBufferArray(scip, &gnodearr, nterms - 1) );
2374 SCIP_CALL( SCIPallocBufferArray(scip, &unsatarcs, nedges) );
2375 SCIP_CALL( SCIPallocBufferArray(scip, &origedge, nedges) );
2376
2377 for( i = 0; i < nedges; i++ )
2378 if( !Is_term(transgraph->term[transgraph->tail[i]]) && transgraph->head[i] == pseudoroot )
2379 origedge[i] = FALSE;
2380 else if( transgraph->tail[i] == pseudoroot && !Is_term(transgraph->term[transgraph->head[i]]) )
2381 origedge[i] = FALSE;
2382 else
2383 origedge[i] = TRUE;
2384
2385 for( i = 0; i < nterms - 1; i++ )
2386 {
2387 SCIP_CALL( SCIPallocBuffer(scip, &gnodearr[i]) ); /*lint !e866*/
2388 }
2389
2390 SCIP_CALL( SCIPpqueueCreate( &pqueue, nnodes, 2.0, GNODECmpByDist, NULL) );
2391
2392 k = 0;
2393 /* mark terminals as active, add all except root to pqueue */
2394 for( i = 0; i < nnodes; i++ )
2395 {
2396 if( Is_term(transgraph->term[i]) )
2397 {
2398 active[i] = TRUE;
2399 assert(transgraph->grad[i] > 0);
2400 if( i != root )
2401 {
2402 gnodearr[k]->number = i;
2403 gnodearr[k]->dist = transgraph->grad[i];
2404
2405 for( a = transgraph->inpbeg[i]; a != EAT_LAST; a = transgraph->ieat[a] )
2406 if( SCIPisEQ(scip, transgraph->cost[a], 0.0) )
2407 break;
2408
2409 if( a != EAT_LAST )
2410 gnodearr[k]->dist += transgraph->grad[transgraph->tail[a]] - 1;
2411
2412 assert(gnodearr[k]->dist > 0);
2413
2414 SCIP_CALL( SCIPpqueueInsert(pqueue, gnodearr[k++]) );
2415 }
2416 }
2417 else
2418 {
2419 active[i] = FALSE;
2420 }
2421 transgraph->mark[i] = FALSE;
2422 }
2423
2424 for( i = 0; i < nedges; i++ )
2425 {
2426 rescap[i] = transgraph->cost[i];
2427 if( SCIPisZero(scip, rescap[i]) )
2428 sat[i] = TRUE;
2429 else
2430 sat[i] = FALSE;
2431 }
2432
2433 /* dual ascent loop */
2434 while( SCIPpqueueNElems(pqueue) > 0 && !SCIPisStopped(scip) )
2435 {
2436 /* get active vertex of minimum score */
2437 gnodeact = (GNODE*) SCIPpqueueRemove(pqueue);
2438
2439 v = gnodeact->number;
2440 prio1 = gnodeact->dist;
2441
2442 firstrun = TRUE;
2443 nunsatarcs = 0;
2444
2445 /* perform augmentation as long as ... */
2446 for( ; ; )
2447 {
2448 assert(stacklength == 0);
2449 /* 1. step: BFS from v (or connected component) on saturated, incoming arcs */
2450
2451 if( firstrun )
2452 {
2453 degsum = transgraph->grad[v];
2454 ncutverts = 0;
2455 firstrun = FALSE;
2456 nunsatarcs = 0;
2457 transgraph->mark[v] = TRUE;
2458 cutverts[ncutverts++] = v;
2459 stackarr[stacklength++] = v;
2460 }
2461 /* not in first processing of root component: */
2462 else
2463 {
2464 for( i = norgcutverts; i < ncutverts; i++ )
2465 {
2466 s = cutverts[i];
2467 assert(transgraph->mark[s]);
2468 if( active[s] )
2469 {
2470 active[v] = FALSE;
2471 stacklength = 0;
2472 goto ENDOFLOOP;
2473 }
2474
2475 stackarr[stacklength++] = s;
2476 }
2477 }
2478
2479 while( stacklength )
2480 {
2481 pnode = stackarr[--stacklength];
2482
2483 /* traverse incoming arcs */
2484 for( a = transgraph->inpbeg[pnode]; a != EAT_LAST; a = transgraph->ieat[a] )
2485 {
2486 tail = transgraph->tail[a];
2487 if( sat[a] )
2488 {
2489 if( !transgraph->mark[tail] )
2490 {
2491 /* if an active vertex has been hit, break */
2492 if( active[tail] )
2493 {
2494 active[v] = FALSE;
2495 stacklength = 0;
2496 goto ENDOFLOOP;
2497 }
2498
2499 degsum += transgraph->grad[tail];
2500 transgraph->mark[tail] = TRUE;
2501 cutverts[ncutverts++] = tail;
2502 stackarr[stacklength++] = tail;
2503 }
2504 }
2505 else if( !transgraph->mark[tail] )
2506 {
2507 unsatarcs[nunsatarcs++] = a;
2508 }
2509 }
2510 }
2511
2512 currscore = degsum - (ncutverts - 1);
2513
2514 assert(SCIPisGE(scip, currscore, prio1));
2515
2516 /* augmentation criteria met? */
2517 if( SCIPisLE(scip, (currscore - prio1) / prio1, maxdeviation) || (SCIPpqueueNElems(pqueue) == 0) )
2518 {
2519 SCIP_Bool in = FALSE;
2520 SCIP_ROW* row;
2521 SCIP_CONS* cons = NULL;
2522
2523 /* 2. pass: get minimum residual capacity among cut-arcs */
2524
2525 /* adjust array of unsatisfied arcs */
2526 min = FARAWAY;
2527 shift = 0;
2528
2529 for( i = 0; i < nunsatarcs; i++ )
2530 {
2531 a = unsatarcs[i];
2532 if( transgraph->mark[transgraph->tail[a]] )
2533 {
2534 shift++;
2535 }
2536 else
2537 {
2538
2539 assert(!sat[a]);
2540 if( SCIPisLT(scip, rescap[a], min) )
2541 min = rescap[a];
2542 if( shift != 0 )
2543 unsatarcs[i - shift] = a;
2544 }
2545 }
2546
2547 assert(SCIPisLT(scip, min, FARAWAY));
2548 nunsatarcs -= shift;
2549
2550 if( nunsatarcs > 0)
2551 assert(!transgraph->mark[transgraph->tail[unsatarcs[nunsatarcs-1]]]);
2552
2553 norgcutverts = ncutverts;
2554
2555
2556 /* 3. pass: perform augmentation */
2557
2558
2559 /* create constraint/row */
2560
2561 if( addcuts )
2562 {
2563 if( addconss )
2564 {
2565 SCIP_CALL( SCIPcreateConsLinear(scip, &cons, "da", 0, NULL, NULL,
2566 1.0, SCIPinfinity(scip), TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, TRUE, TRUE, FALSE) );
2567 }
2568 else
2569 {
2570 SCIP_CALL(SCIPcreateEmptyRowConshdlr(scip, &row, conshdlr, "da", 1.0, SCIPinfinity(scip), FALSE, FALSE, TRUE));
2571 SCIP_CALL(SCIPcacheRowExtensions(scip, row));
2572 }
2573 }
2574
2575 dualobj += min;
2576 for( i = 0; i < nunsatarcs; i++ )
2577 {
2578 a = unsatarcs[i];
2579 if( a == -1 )
2580 continue;
2581
2582 if( addcuts && origedge[a] )
2583 {
2584 assert(vars != NULL);
2585 assert(cons != NULL);
2586
2587 if( g->tail[a] == root && g->head[a] == v )
2588 in = TRUE;
2589
2590 if( addconss )
2591 SCIP_CALL( SCIPaddCoefLinear(scip, cons, vars[a], 1.0) );
2592 else
2593 SCIP_CALL( SCIPaddVarToRow(scip, row, vars[a], 1.0) );
2594 }
2595 rescap[a] -= min;
2596
2597 assert(SCIPisGE(scip, rescap[a], 0.0));
2598
2599 if( SCIPisEQ(scip, rescap[a], 0.0) )
2600 {
2601 sat[a] = TRUE;
2602 if( !(transgraph->mark[transgraph->tail[a]]) )
2603 {
2604 tail = transgraph->tail[a];
2605 degsum += transgraph->grad[tail];
2606 transgraph->mark[tail] = TRUE;
2607 cutverts[ncutverts++] = tail;
2608 }
2609 }
2610 }
2611
2612 if( addcuts )
2613 {
2614 assert(vars != NULL);
2615
2616 if( !in )
2617 {
2618 for( i = g->outbeg[root]; i != EAT_LAST; i = g->oeat[i] )
2619 if( g->head[i] == v )
2620 {
2621 if( addconss )
2622 SCIP_CALL( SCIPaddCoefLinear(scip, cons, vars[i], 1.0) );
2623 else
2624 SCIP_CALL( SCIPaddVarToRow(scip, row, vars[i], 1.0) );
2625 }
2626 }
2627
2628 if( addconss )
2629 {
2630 SCIP_CALL( SCIPaddCons(scip, cons) );
2631 SCIP_CALL( SCIPreleaseCons(scip, &cons) );
2632 }
2633 else
2634 {
2635 SCIP_Bool infeasible;
2636 assert(row != NULL);
2637
2638 SCIP_CALL( SCIPflushRowExtensions(scip, row) );
2639 SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
2640 SCIP_CALL( SCIPreleaseRow(scip, &row) );
2641
2642 assert(!infeasible);
2643 }
2644 }
2645
2646 continue;
2647 }
2648 else
2649 {
2650 /* reinsert active vertex */
2651 gnodeact->dist = currscore;
2652 SCIP_CALL( SCIPpqueueInsert(pqueue, gnodeact) );
2653 }
2654
2655 ENDOFLOOP:
2656
2657 for( i = 0; i < ncutverts; i++ )
2658 transgraph->mark[cutverts[i]] = FALSE;
2659
2660 break;
2661 } /* augmentation loop */
2662 } /* dual ascent loop */
2663
2664
2665 *objval = dualobj + offset;
2666 SCIPdebugMessage("DA: dualglobal: %f \n", *objval + SCIPprobdataGetOffset(scip));
2667
2668 /* call dual Ascend-And-Prune? */
2669 if( ascendandprune )
2670 {
2671 SCIP_Bool success;
2672 SCIP_CALL( SCIPStpHeurAscendPruneRun(scip, NULL, g, rescap, unsatarcs, cutverts, -1, active, &success, TRUE));
2673 }
2674
2675 /* free memory */
2676 SCIPpqueueFree(&pqueue);
2677
2678 for( i = nterms - 2; i >= 0; i-- )
2679 SCIPfreeBuffer(scip, &gnodearr[i]);
2680
2681 SCIPfreeBufferArray(scip, &origedge);
2682 SCIPfreeBufferArray(scip, &unsatarcs);
2683 SCIPfreeBufferArray(scip, &cutverts);
2684 SCIPfreeBufferArray(scip, &gnodearr);
2685 SCIPfreeBufferArray(scip, &active);
2686 SCIPfreeBufferArray(scip, &sat);
2687 SCIPfreeBufferArray(scip, &stackarr);
2688
2689 if( redcost == NULL )
2690 SCIPfreeBufferArray(scip, &rescap);
2691
2692 graph_free(scip, &transgraph, TRUE);
2693
2694 return SCIP_OKAY;
2695
2696 }
2697
2698 /**@} */
2699