Scippy

SCIP

Solving Constraint Integer Programs

cons_stp.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2020 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file 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)
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
110 /* do depth-first search */
111 #define DFS
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
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
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
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
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
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
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) );
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
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
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 */
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
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
1244 SCIP_DECL_CONSINITSOL(consInitsolStp)
1245 { /*lint --e{715}*/
1246 #ifdef WITH_UG
1248 #endif
1249  return SCIP_OKAY;
1250 }
1251 
1252 /** frees specific constraint data */
1253 static
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
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
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
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
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
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
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
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
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
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 */
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 */
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,
1599  SCIP_CALL( SCIPsetConshdlrSepa(scip, conshdlr, consSepalpStp, NULL, CONSHDLR_SEPAFREQ,
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 */
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 */
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 */
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 
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 */
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) );
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 /**@} */
SCIP_Bool SCIPisEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
void * SCIPpqueueRemove(SCIP_PQUEUE *pqueue)
Definition: misc.c:1433
void SCIPconshdlrSetData(SCIP_CONSHDLR *conshdlr, SCIP_CONSHDLRDATA *conshdlrdata)
Definition: cons.c:4209
static int graph_next_term(const GRAPH *g, int terms, int *term, const int *w, const SCIP_Bool firstrun)
Definition: cons_stp.c:277
SCIP_RETCODE SCIPcreateConsLinear(SCIP *scip, SCIP_CONS **cons, const char *name, int nvars, SCIP_VAR **vars, SCIP_Real *vals, SCIP_Real lhs, SCIP_Real rhs, SCIP_Bool initial, SCIP_Bool separate, SCIP_Bool enforce, SCIP_Bool check, SCIP_Bool propagate, SCIP_Bool local, SCIP_Bool modifiable, SCIP_Bool dynamic, SCIP_Bool removable, SCIP_Bool stickingatnode)
#define CONSHDLR_PROPFREQ
Definition: cons_stp.c:76
static volatile int nterms
Definition: interrupt.c:38
int *RESTRICT mincut_e
Definition: grph.h:113
#define DEFAULT_CREEPFLOW
Definition: cons_stp.c:93
static SCIP_DECL_CONSINITSOL(consInitsolStp)
Definition: cons_stp.c:1245
SCIP_Bool SCIPisStopped(SCIP *scip)
Definition: scip_general.c:687
void SCIPpqueueFree(SCIP_PQUEUE **pqueue)
Definition: misc.c:1262
SCIP_RETCODE SCIPsetConshdlrTrans(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_DECL_CONSTRANS((*constrans)))
Definition: scip_cons.c:586
SCIP_CONSHDLR * SCIPfindConshdlr(SCIP *scip, const char *name)
Definition: scip_cons.c:877
#define DEFAULT_NESTEDCUT
Definition: cons_stp.c:95
int *RESTRICT head
Definition: grph.h:96
Definition: grph.h:57
int source
Definition: grph.h:67
SCIP_Bool SCIPisGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
void * SCIPpqueueFirst(SCIP_PQUEUE *pqueue)
Definition: misc.c:1453
SCIP_RETCODE SCIPcreateCons(SCIP *scip, SCIP_CONS **cons, const char *name, SCIP_CONSHDLR *conshdlr, SCIP_CONSDATA *consdata, SCIP_Bool initial, SCIP_Bool separate, SCIP_Bool enforce, SCIP_Bool check, SCIP_Bool propagate, SCIP_Bool local, SCIP_Bool modifiable, SCIP_Bool dynamic, SCIP_Bool removable, SCIP_Bool stickingatnode)
Definition: scip_cons.c:934
static SCIP_DECL_CONSENFOLP(consEnfolpStp)
Definition: cons_stp.c:1401
SCIP_Bool SCIPconsIsLocal(SCIP_CONS *cons)
Definition: cons.c:8328
static SCIP_DECL_CONSSEPALP(consSepalpStp)
Definition: cons_stp.c:1323
Constraint handler for Steiner problems.
SCIP_Bool SCIPconsIsChecked(SCIP_CONS *cons)
Definition: cons.c:8288
#define SCIPfreeMemoryArray(scip, ptr)
Definition: scip_mem.h:69
SCIP_RETCODE SCIPincludeConshdlrBasic(SCIP *scip, SCIP_CONSHDLR **conshdlrptr, const char *name, const char *desc, int enfopriority, int chckpriority, int eagerfreq, SCIP_Bool needscons, SCIP_DECL_CONSENFOLP((*consenfolp)), SCIP_DECL_CONSENFOPS((*consenfops)), SCIP_DECL_CONSCHECK((*conscheck)), SCIP_DECL_CONSLOCK((*conslock)), SCIP_CONSHDLRDATA *conshdlrdata)
Definition: scip_cons.c:166
SCIP_RETCODE SCIPaddVarLocksType(SCIP *scip, SCIP_VAR *var, SCIP_LOCKTYPE locktype, int nlocksdown, int nlocksup)
Definition: scip_var.c:4263
#define DEFAULT_BACKCUT
Definition: cons_stp.c:92
int terms
Definition: grph.h:64
SCIP_VAR ** SCIPprobdataGetVars(SCIP *scip)
#define SCIPallocMemoryArray(scip, ptr, num)
Definition: scip_mem.h:53
SCIP_RETCODE SCIPprintRow(SCIP *scip, SCIP_ROW *row, FILE *file)
Definition: scip_lp.c:2152
SCIP_RETCODE SCIPaddCoefLinear(SCIP *scip, SCIP_CONS *cons, SCIP_VAR *var, SCIP_Real val)
SCIP_RETCODE SCIPStpDualAscentPcMw(SCIP *scip, GRAPH *g, SCIP_Real *redcost, SCIP_Real *objval, SCIP_Bool addcuts, SCIP_Bool ascendandprune, int nruns)
Definition: cons_stp.c:2268
int *RESTRICT maxdeg
Definition: grph.h:78
#define EAT_LAST
Definition: grph.h:31
int SCIPconshdlrGetNConss(SCIP_CONSHDLR *conshdlr)
Definition: cons.c:4594
static SCIP_RETCODE sep_flow(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_CONSHDLRDATA *conshdlrdata, SCIP_CONSDATA *consdata, int maxcuts, int *ncuts)
Definition: cons_stp.c:392
SCIP_NODE * SCIPgetCurrentNode(SCIP *scip)
Definition: scip_tree.c:81
reduction and dual-cost based primal heuristic for Steiner problems
void graph_mincut_exec(const GRAPH *, const int, const int, const int, const int, const int, const int *, const int *, int *RESTRICT, const int *, const int *, const int *, const SCIP_Bool)
#define FALSE
Definition: def.h:73
int *RESTRICT inpbeg
Definition: grph.h:74
#define FLOW_FACTOR
Definition: cons_stp.c:108
void SCIPStpBranchruleInitNodeState(const GRAPH *g, int *nodestate)
Definition: branch_stp.c:683
SCIP_EXPORT int SCIPnodeGetDepth(SCIP_NODE *node)
Definition: tree.c:7430
void graph_free(SCIP *, GRAPH **, SCIP_Bool)
Definition: grphbase.c:3678
SCIP_Bool SCIPisFeasNegative(SCIP *scip, SCIP_Real val)
Problem data for stp problem.
#define TRUE
Definition: def.h:72
#define SCIPdebug(x)
Definition: pub_message.h:84
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
int SCIPpqueueNElems(SCIP_PQUEUE *pqueue)
Definition: misc.c:1467
SCIP_Bool SCIPconsIsInitial(SCIP_CONS *cons)
Definition: cons.c:8258
SCIP_PROBDATA * SCIPgetProbData(SCIP *scip)
Definition: scip_prob.c:962
int SCIPprobdataGetNVars(SCIP *scip)
SCIP_RETCODE SCIPsetConshdlrDelete(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_DECL_CONSDELETE((*consdelete)))
Definition: scip_cons.c:563
static GRAPHNODE ** active
SCIP_RETCODE SCIPaddBoolParam(SCIP *scip, const char *name, const char *desc, SCIP_Bool *valueptr, SCIP_Bool isadvanced, SCIP_Bool defaultvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition: scip_param.c:48
#define SCIPfreeBlockMemory(scip, ptr)
Definition: scip_mem.h:95
#define SCIPdebugMessage
Definition: pub_message.h:87
SCIP_RETCODE SCIPflushRowExtensions(SCIP *scip, SCIP_ROW *row)
Definition: scip_lp.c:1604
#define DEFAULT_MAXROUNDS
Definition: cons_stp.c:84
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip_mem.h:123
int number
Definition: misc_stp.h:43
#define SCIPallocBlockMemory(scip, ptr)
Definition: scip_mem.h:78
static SCIP_DECL_CONSFREE(consFreeStp)
Definition: cons_stp.c:1224
#define DEFAULT_DISJUNCTCUT
Definition: cons_stp.c:94
#define CONSHDLR_CHECKPRIORITY
Definition: cons_stp.c:74
static SCIP_DECL_CONSLOCK(consLockStp)
Definition: cons_stp.c:1522
void graph_knot_chg(GRAPH *, int, int)
Definition: grphbase.c:2222
GRAPH * SCIPprobdataGetGraph(SCIP_PROBDATA *probdata)
SCIP_RETCODE 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)
Definition: cons_stp.c:1687
SCIP_RETCODE SCIPStpHeurAscendPruneRun(SCIP *scip, SCIP_HEUR *heur, const GRAPH *g, const SCIP_Real *redcosts, int *edgearrint, int *nodearrint, int root, STP_Bool *nodearrchar, SCIP_Bool *solfound, SCIP_Bool addsol)
SCIP_RETCODE SCIPsetConshdlrFree(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_DECL_CONSFREE((*consfree)))
Definition: scip_cons.c:357
static SCIP_Bool is_active(const int *active, int realtail, int dfsbase)
Definition: cons_stp.c:157
static SCIP_DECL_CONSHDLRCOPY(conshdlrCopyStp)
Definition: cons_stp.c:1208
int *RESTRICT mark
Definition: grph.h:70
SCIP_Bool SCIPisLE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
#define DA_EPS
Definition: cons_stp.c:101
SCIP_VAR * w
Definition: circlepacking.c:58
int *RESTRICT oeat
Definition: grph.h:103
#define CONNECT
Definition: grph.h:154
SCIP_RETCODE SCIPcacheRowExtensions(SCIP *scip, SCIP_ROW *row)
Definition: scip_lp.c:1581
int *RESTRICT mincut_dist
Definition: grph.h:106
SCIP_RETCODE SCIPcreateEmptyRowConshdlr(SCIP *scip, SCIP_ROW **row, SCIP_CONSHDLR *conshdlr, const char *name, SCIP_Real lhs, SCIP_Real rhs, SCIP_Bool local, SCIP_Bool modifiable, SCIP_Bool removable)
Definition: scip_lp.c:1337
#define CONSHDLR_DESC
Definition: cons_stp.c:71
int stp_type
Definition: grph.h:127
#define SCIPerrorMessage
Definition: pub_message.h:55
SCIP_RETCODE SCIPreleaseRow(SCIP *scip, SCIP_ROW **row)
Definition: scip_lp.c:1508
unsigned char STP_Bool
Definition: grph.h:52
SCIP_Bool SCIPisFeasEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Real SCIPprobdataGetOffset(SCIP *scip)
#define CONSHDLR_SEPAPRIORITY
Definition: cons_stp.c:72
#define SCIPallocBuffer(scip, ptr)
Definition: scip_mem.h:109
SCIP_RETCODE SCIPsetConshdlrInitsol(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_DECL_CONSINITSOL((*consinitsol)))
Definition: scip_cons.c:429
#define STP_DCSTP
Definition: grph.h:43
#define SCIPfreeBufferArrayNull(scip, ptr)
Definition: scip_mem.h:124
SCIP_Real dist
Definition: misc_stp.h:44
SCIP_RETCODE SCIPincludeConshdlrStp(SCIP *scip)
Definition: cons_stp.c:1575
#define DEFAULT_MAXROUNDSROOT
Definition: cons_stp.c:85
int *RESTRICT grad
Definition: grph.h:73
SCIP_Bool SCIPisZero(SCIP *scip, SCIP_Real val)
#define DA_MAXDEVIATION_UPPER
Definition: cons_stp.c:100
int *RESTRICT mincut_head_inact
Definition: grph.h:108
internal miscellaneous methods
#define NULL
Definition: lpi_spx1.cpp:155
int knots
Definition: grph.h:62
SCIP_CONSDATA * SCIPconsGetData(SCIP_CONS *cons)
Definition: cons.c:8119
#define SCIP_CALL(x)
Definition: def.h:364
#define CONSHDLR_DELAYPROP
Definition: cons_stp.c:81
struct SCIP_ConsData SCIP_CONSDATA
Definition: type_cons.h:56
static SCIP_RETCODE 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)
Definition: cons_stp.c:177
reduction-based primal heuristic for Steiner problems
#define FARAWAY
Definition: grph.h:156
int *RESTRICT mincut_head
Definition: grph.h:107
#define CONSHDLR_NAME
Definition: cons_stp.c:70
int GNODECmpByDist(void *first_arg, void *second_arg)
static SCIP_DECL_CONSENFOPS(consEnfopsStp)
Definition: cons_stp.c:1426
static SCIP_RETCODE 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)
Definition: cons_stp.c:1081
SCIP_Bool SCIPconsIsRemovable(SCIP_CONS *cons)
Definition: cons.c:8358
#define STP_SPG
Definition: grph.h:38
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:111
SCIP_Real SCIPinfinity(SCIP *scip)
SCIP_RETCODE SCIPsetConshdlrSepa(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_DECL_CONSSEPALP((*conssepalp)), SCIP_DECL_CONSSEPASOL((*conssepasol)), int sepafreq, int sepapriority, SCIP_Bool delaysepa)
Definition: scip_cons.c:220
#define DEFAULT_MAXSEPACUTS
Definition: cons_stp.c:86
static void set_capacity(const GRAPH *g, const SCIP_Bool creep_flow, const int flip, int *capa, const SCIP_Real *xval)
Definition: cons_stp.c:349
static SCIP_DECL_CONSCOPY(consCopyStp)
Definition: cons_stp.c:1544
int SCIPgetDepth(SCIP *scip)
Definition: scip_tree.c:638
#define SCIP_Bool
Definition: def.h:70
int *RESTRICT ieat
Definition: grph.h:102
static SCIP_DECL_CONSTRANS(consTransStp)
Definition: cons_stp.c:1269
#define DEFAULT_DAMAXDEVIATION
Definition: cons_stp.c:98
SCIP_RETCODE SCIPsetConshdlrCopy(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_DECL_CONSHDLRCOPY((*conshdlrcopy)), SCIP_DECL_CONSCOPY((*conscopy)))
Definition: scip_cons.c:332
int *RESTRICT tail
Definition: grph.h:95
static SCIP_RETCODE sep_2cut(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_CONSHDLRDATA *conshdlrdata, SCIP_CONSDATA *consdata, const int *termorg, int maxcuts, int *ncuts)
Definition: cons_stp.c:649
static SCIP_DECL_CONSDELETE(consDeleteStp)
Definition: cons_stp.c:1255
#define BRANCH_STP_VERTEX_TERM
Definition: branch_stp.h:38
SCIP_CONSHDLRDATA * SCIPconshdlrGetData(SCIP_CONSHDLR *conshdlr)
Definition: cons.c:4199
#define MAX(x, y)
Definition: tclique_def.h:83
SCIP_RETCODE SCIPStpValidateSol(SCIP *, const GRAPH *, const double *, SCIP_Bool *)
Definition: validate.c:136
SCIP_Bool SCIPisCutEfficacious(SCIP *scip, SCIP_SOL *sol, SCIP_ROW *cut)
Definition: scip_cut.c:88
int *RESTRICT term
Definition: grph.h:68
#define CONSHDLR_EAGERFREQ
Definition: cons_stp.c:77
#define BMScopyMemoryArray(ptr, source, num)
Definition: memory.h:126
SCIP_Bool SCIPconsIsSeparated(SCIP_CONS *cons)
Definition: cons.c:8268
#define CREEP_VALUE
Definition: cons_stp.c:109
Constraint handler for linear constraints in their most general form, .
SCIP_Bool SCIPisFeasGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
#define DA_MAXDEVIATION_LOWER
Definition: cons_stp.c:99
SCIP_Bool SCIPconsIsDynamic(SCIP_CONS *cons)
Definition: cons.c:8348
#define DEFAULT_FLOWSEP
Definition: cons_stp.c:96
SCIP_RETCODE SCIPcreateConsStp(SCIP *scip, SCIP_CONS **cons, const char *name, GRAPH *graph)
Definition: cons_stp.c:1636
includes various files containing graph methods used for Steiner tree problems
Portable defintions.
int *RESTRICT mincut_numb
Definition: grph.h:109
int layers
Definition: grph.h:65
#define SCIPfreeMemory(scip, ptr)
Definition: scip_mem.h:67
SCIP_EXPORT SCIP_Real SCIPvarGetUbLocal(SCIP_VAR *var)
Definition: var.c:17728
Steiner vertex branching rule.
SCIP_RETCODE SCIPaddRow(SCIP *scip, SCIP_ROW *row, SCIP_Bool forcecut, SCIP_Bool *infeasible)
Definition: scip_cut.c:221
#define SCIPfreeBuffer(scip, ptr)
Definition: scip_mem.h:121
#define DEFAULT_MAXSEPACUTSROOT
Definition: cons_stp.c:87
#define Is_term(a)
Definition: grph.h:168
SCIP_RETCODE SCIPpqueueCreate(SCIP_PQUEUE **pqueue, int initsize, SCIP_Real sizefac, SCIP_DECL_SORTPTRCOMP((*ptrcomp)), SCIP_DECL_PQUEUEELEMCHGPOS((*elemchgpos)))
Definition: misc.c:1235
SCIP_Bool SCIPisFeasGT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
struct SCIP_ProbData SCIP_PROBDATA
Definition: type_prob.h:44
void graph_get_csr(const GRAPH *, int *RESTRICT, int *RESTRICT, int *RESTRICT, int *)
static SCIP_DECL_CONSINITLP(consInitlpStp)
Definition: cons_stp.c:1301
SCIP_RETCODE SCIPStpBranchruleApplyVertexChgs(SCIP *scip, int *vertexchgs, GRAPH *graph)
Definition: branch_stp.c:708
#define CONSHDLR_DELAYSEPA
Definition: cons_stp.c:80
SCIP_Bool SCIPconsIsEnforced(SCIP_CONS *cons)
Definition: cons.c:8278
SCIP_Real * cost
Definition: grph.h:94
#define CONSHDLR_SEPAFREQ
Definition: cons_stp.c:75
SCIP_Bool SCIPisFeasLT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_RETCODE SCIPaddVarToRow(SCIP *scip, SCIP_ROW *row, SCIP_VAR *var, SCIP_Real val)
Definition: scip_lp.c:1641
SCIP_CONS ** SCIPconshdlrGetConss(SCIP_CONSHDLR *conshdlr)
Definition: cons.c:4551
SCIP_VAR * a
Definition: circlepacking.c:57
SCIP_RETCODE SCIPaddIntParam(SCIP *scip, const char *name, const char *desc, int *valueptr, SCIP_Bool isadvanced, int defaultvalue, int minvalue, int maxvalue, SCIP_DECL_PARAMCHGD((*paramchgd)), SCIP_PARAMDATA *paramdata)
Definition: scip_param.c:74
const char * SCIPconsGetName(SCIP_CONS *cons)
Definition: cons.c:8089
#define SCIP_Real
Definition: def.h:163
#define CONSHDLR_PROP_TIMING
Definition: cons_stp.c:90
SCIP_Bool SCIPisLT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_Real * SCIPprobdataGetXval(SCIP *scip, SCIP_SOL *sol)
int *RESTRICT outbeg
Definition: grph.h:76
#define CONSHDLR_ENFOPRIORITY
Definition: cons_stp.c:73
SCIP_Bool SCIPconsIsStickingAtNode(SCIP_CONS *cons)
Definition: cons.c:8368
int edges
Definition: grph.h:92
SCIP_RETCODE graph_pc_getSap(SCIP *, GRAPH *, GRAPH **, SCIP_Real *)
Definition: grphbase.c:1075
const char * SCIPconshdlrGetName(SCIP_CONSHDLR *conshdlr)
Definition: cons.c:4179
#define flipedge(edge)
Definition: grph.h:150
#define Q_NULL
Definition: cons_stp.c:63
#define SCIPallocMemory(scip, ptr)
Definition: scip_mem.h:51
SCIP_RETCODE SCIPaddCons(SCIP *scip, SCIP_CONS *cons)
Definition: scip_prob.c:2764
SCIP_RETCODE SCIPaddPoolCut(SCIP *scip, SCIP_ROW *row)
Definition: scip_cut.c:332
struct SCIP_ConshdlrData SCIP_CONSHDLRDATA
Definition: type_cons.h:55
int *RESTRICT mincut_r
Definition: grph.h:115
#define nnodes
Definition: gastrans.c:65
#define CONSHDLR_NEEDSCONS
Definition: cons_stp.c:82
static SCIP_DECL_CONSCHECK(consCheckStp)
Definition: cons_stp.c:1451
SCIP_RETCODE SCIPreleaseCons(SCIP *scip, SCIP_CONS **cons)
Definition: scip_cons.c:1110
SCIP_STAGE SCIPgetStage(SCIP *scip)
Definition: scip_general.c:356
SCIP_RETCODE SCIPsetConshdlrProp(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_DECL_CONSPROP((*consprop)), int propfreq, SCIP_Bool delayprop, SCIP_PROPTIMING proptiming)
Definition: scip_cons.c:266
GRAPH * SCIPprobdataGetGraph2(SCIP *scip)
SCIP_Bool SCIPconsIsPropagated(SCIP_CONS *cons)
Definition: cons.c:8308
static SCIP_DECL_CONSPROP(consPropStp)
Definition: cons_stp.c:1473
SCIP_RETCODE SCIPsetConshdlrInitlp(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_DECL_CONSINITLP((*consinitlp)))
Definition: scip_cons.c:609
SCIP callable library.
SCIP_Bool SCIPconsIsModifiable(SCIP_CONS *cons)
Definition: cons.c:8338
SCIP_RETCODE SCIPpqueueInsert(SCIP_PQUEUE *pqueue, void *elem)
Definition: misc.c:1334
void SCIPStpConshdlrSetGraph(SCIP *scip, const GRAPH *g)
Definition: cons_stp.c:1666