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-2018 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 scip.zib.de. */
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( SCIPcreateEmptyRowCons(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( SCIPcreateEmptyRowCons(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( SCIPcreateEmptyRowCons(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( SCIPcreateEmptyRowCons(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( SCIPcreateEmptyRowCons(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  SCIP_CALL( SCIPaddVarLocksType(scip, vars[v], SCIP_LOCKTYPE_MODEL, 1, 1) );
1535 
1536  return SCIP_OKAY;
1537 }
1538 
1539 /** constraint copying method of constraint handler */
1540 static
1541 SCIP_DECL_CONSCOPY(consCopyStp)
1542 { /*lint --e{715}*/
1543  const char* consname;
1544  SCIP_PROBDATA* probdata;
1545  GRAPH* graph;
1546 
1547  probdata = SCIPgetProbData(scip);
1548  assert(probdata != NULL);
1549 
1550  graph = SCIPprobdataGetGraph(probdata);
1551  assert(graph != NULL);
1552 
1553  consname = SCIPconsGetName(sourcecons);
1554 
1555  /* creates and captures a and constraint */
1556  SCIP_CALL( SCIPcreateConsStp(scip, cons, consname, graph) );
1557 
1558  *valid = TRUE;
1559 
1560  return SCIP_OKAY;
1561 }
1562 
1563 
1564 /**@} */
1565 
1566 /**@name Interface methods
1567  *
1568  * @{
1569  */
1570 
1571 /** creates the handler for stp constraints and includes it in SCIP */
1573  SCIP* scip /**< SCIP data structure */
1574  )
1575 {
1576  SCIP_CONSHDLRDATA* conshdlrdata;
1577  SCIP_CONSHDLR* conshdlr;
1578 
1579  /* create stp constraint handler data */
1580  SCIP_CALL( SCIPallocMemory(scip, &conshdlrdata) );
1581 
1582  conshdlr = NULL;
1583  /* include constraint handler */
1586  consEnfolpStp, consEnfopsStp, consCheckStp, consLockStp,
1587  conshdlrdata) );
1588  assert(conshdlr != NULL);
1589 
1590  SCIP_CALL( SCIPsetConshdlrCopy(scip, conshdlr, conshdlrCopyStp, consCopyStp) );
1591  SCIP_CALL( SCIPsetConshdlrDelete(scip, conshdlr, consDeleteStp) );
1592  SCIP_CALL( SCIPsetConshdlrTrans(scip, conshdlr, consTransStp) );
1593  SCIP_CALL( SCIPsetConshdlrInitsol(scip, conshdlr, consInitsolStp) );
1594  SCIP_CALL( SCIPsetConshdlrInitlp(scip, conshdlr, consInitlpStp) );
1595  SCIP_CALL( SCIPsetConshdlrProp(scip, conshdlr, consPropStp, CONSHDLR_PROPFREQ, CONSHDLR_DELAYPROP,
1597  SCIP_CALL( SCIPsetConshdlrSepa(scip, conshdlr, consSepalpStp, NULL, CONSHDLR_SEPAFREQ,
1599  SCIP_CALL( SCIPsetConshdlrFree(scip, conshdlr, consFreeStp) );
1600 
1601  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/stp/backcut", "Try Back-Cuts",
1602  &conshdlrdata->backcut, TRUE, DEFAULT_BACKCUT, NULL, NULL) );
1603  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/stp/creepflow", "Use Creep-Flow",
1604  &conshdlrdata->creepflow, TRUE, DEFAULT_CREEPFLOW, NULL, NULL) );
1605  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/stp/disjunctcut", "Only disjunct Cuts",
1606  &conshdlrdata->disjunctcut, TRUE, DEFAULT_DISJUNCTCUT, NULL, NULL) );
1607  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/stp/nestedcut", "Try Nested-Cuts",
1608  &conshdlrdata->nestedcut, TRUE, DEFAULT_NESTEDCUT, NULL, NULL) );
1609  SCIP_CALL( SCIPaddBoolParam(scip, "constraints/stp/flowsep", "Try Flow-Cuts",
1610  &conshdlrdata->flowsep, TRUE, DEFAULT_FLOWSEP, NULL, NULL) );
1611  SCIP_CALL( SCIPaddIntParam(scip,
1612  "constraints/"CONSHDLR_NAME"/maxrounds",
1613  "maximal number of separation rounds per node (-1: unlimited)",
1614  &conshdlrdata->maxrounds, FALSE, DEFAULT_MAXROUNDS, -1, INT_MAX, NULL, NULL) );
1615  SCIP_CALL( SCIPaddIntParam(scip,
1616  "constraints/"CONSHDLR_NAME"/maxroundsroot",
1617  "maximal number of separation rounds per node in the root node (-1: unlimited)",
1618  &conshdlrdata->maxroundsroot, FALSE, DEFAULT_MAXROUNDSROOT, -1, INT_MAX, NULL, NULL) );
1619  SCIP_CALL( SCIPaddIntParam(scip,
1620  "constraints/"CONSHDLR_NAME"/maxsepacuts",
1621  "maximal number of cuts separated per separation round",
1622  &conshdlrdata->maxsepacuts, FALSE, DEFAULT_MAXSEPACUTS, 0, INT_MAX, NULL, NULL) );
1623  SCIP_CALL( SCIPaddIntParam(scip,
1624  "constraints/"CONSHDLR_NAME"/maxsepacutsroot",
1625  "maximal number of cuts separated per separation round in the root node",
1626  &conshdlrdata->maxsepacutsroot, FALSE, DEFAULT_MAXSEPACUTSROOT, 0, INT_MAX, NULL, NULL) );
1627 
1628 
1629  return SCIP_OKAY;
1630 }
1631 
1632 /** creates and captures a stp constraint */
1634  SCIP* scip, /**< SCIP data structure */
1635  SCIP_CONS** cons, /**< pointer to hold the created constraint */
1636  const char* name, /**< name of constraint */
1637  GRAPH* graph /**< graph data structure */
1638  )
1639 {
1640  SCIP_CONSHDLR* conshdlr;
1641  SCIP_CONSDATA* consdata;
1642 
1643  /* find the stp constraint handler */
1644  conshdlr = SCIPfindConshdlr(scip, CONSHDLR_NAME);
1645  if( conshdlr == NULL )
1646  {
1647  SCIPerrorMessage("stp constraint handler not found\n");
1648  return SCIP_PLUGINNOTFOUND;
1649  }
1650 
1651  SCIP_CALL( SCIPallocBlockMemory(scip, &consdata) );
1652 
1653  consdata->graph = graph;
1654 
1655  /* create constraint */
1656  SCIP_CALL( SCIPcreateCons(scip, cons, name, conshdlr, consdata, FALSE, TRUE, TRUE, TRUE, TRUE,
1657  FALSE, FALSE, FALSE, FALSE, FALSE) );
1658 
1659  return SCIP_OKAY;
1660 }
1661 
1662 /** sets graph */
1664  SCIP* scip, /**< SCIP data structure */
1665  const GRAPH* g /**< graph data structure */
1666  )
1667 {
1668  SCIP_CONSDATA* consdata;
1669  SCIP_CONSHDLR* conshdlr;
1670 
1671  conshdlr = SCIPfindConshdlr(scip, "stp");
1672  assert(conshdlr != NULL);
1673  assert(SCIPconshdlrGetNConss(conshdlr) > 0);
1674 
1675  consdata = SCIPconsGetData(SCIPconshdlrGetConss(conshdlr)[0]);
1676 
1677  assert(consdata != NULL);
1678 
1679  consdata->graph = SCIPprobdataGetGraph2(scip);
1680  assert(consdata->graph != NULL);
1681 }
1682 
1683 /* dual ascent heuristic */
1685  SCIP* scip, /**< SCIP data structure */
1686  const GRAPH* g, /**< graph data structure */
1687  SCIP_Real* RESTRICT redcost, /**< array to store reduced costs or NULL */
1688  SCIP_Real* RESTRICT nodearrreal, /**< real vertices array for internal computations or NULL */
1689  SCIP_Real* objval, /**< pointer to store objective value */
1690  SCIP_Bool addcuts, /**< should dual ascent add Steiner cuts? */
1691  SCIP_Bool ascendandprune, /**< should the ascent-and-prune heuristic be executed? */
1692  GNODE** gnodearrterms, /**< gnode terminals array for internal computations or NULL */
1693  const int* result, /**< solution array or NULL */
1694  int* RESTRICT edgearrint, /**< int edges array for internal computations or NULL */
1695  int* RESTRICT nodearrint, /**< int vertices array for internal computations or NULL */
1696  int root, /**< the root */
1697  SCIP_Bool is_pseudoroot, /**< is the root a pseudo root? */
1698  SCIP_Real damaxdeviation, /**< maximum deviation for dual-ascent ( -1.0 for default) */
1699  STP_Bool* RESTRICT nodearrchar /**< STP_Bool vertices array for internal computations or NULL */
1700  )
1701 {
1702  SCIP_CONSHDLR* conshdlr = NULL;
1703  SCIP_PQUEUE* pqueue;
1704  SCIP_VAR** vars;
1705  SCIP_Real* RESTRICT rescap;
1706  GNODE** gnodearr;
1707  int* RESTRICT edgearr;
1708  int* RESTRICT tailarr;
1709  int* RESTRICT start;
1710  int* RESTRICT stackarr;
1711  int* RESTRICT cutverts;
1712  int* RESTRICT unsatarcs;
1713  int* RESTRICT unsattails;
1714  int* RESTRICT gmark;
1715  int* RESTRICT active;
1716  SCIP_Real dualobj;
1717  SCIP_Real currscore;
1718  const SCIP_Real maxdeviation = (damaxdeviation > 0.0) ? damaxdeviation : DEFAULT_DAMAXDEVIATION;
1719  const int nnodes = g->knots;
1720  const int nterms = g->terms;
1721  const int nedges = g->edges;
1722  int ncsredges;
1723  int norgcutverts;
1724  int stacklength;
1725  int augmentingcomponent;
1726  const SCIP_Bool addconss = (SCIPgetStage(scip) < SCIP_STAGE_INITSOLVE);
1727 
1728  /* should currently not be activated */
1729  assert(addconss || !addcuts);
1730  assert(g != NULL);
1731  assert(scip != NULL);
1732  assert(objval != NULL);
1733  assert(Is_term(g->term[root]));
1734  assert(maxdeviation >= DA_MAXDEVIATION_LOWER && maxdeviation <= DA_MAXDEVIATION_UPPER);
1735  assert(damaxdeviation == -1.0 || damaxdeviation > 0.0);
1736 
1737  if( nnodes == 1 )
1738  return SCIP_OKAY;
1739 
1740  if( addcuts )
1741  {
1742  vars = SCIPprobdataGetVars(scip);
1743  assert(vars != NULL);
1744 
1745  if( !addconss )
1746  {
1747  conshdlr = SCIPfindConshdlr(scip, "stp");
1748  assert(conshdlr != NULL);
1749  }
1750  }
1751  else
1752  {
1753  vars = NULL;
1754  }
1755 
1756  /* if specified root is not a terminal, take default root */
1757  if( !Is_term(g->term[root]) )
1758  root = g->source;
1759 
1760 #ifdef BITFIELDSARRAY
1761  u_int32_t* bitarr;
1762  SCIP_CALL( SCIPallocBufferArray(scip, &bitarr, nedges / ARRLENGTH + 1) );
1763 #endif
1764 
1765  stacklength = 0;
1766 
1767  SCIP_CALL( SCIPallocBufferArray(scip, &unsattails, nedges) );
1768 
1769  if( redcost == NULL )
1770  SCIP_CALL( SCIPallocBufferArray(scip, &rescap, nedges) );
1771  else
1772  rescap = redcost;
1773 
1774  if( nodearrint == NULL )
1775  SCIP_CALL( SCIPallocBufferArray(scip, &cutverts, nnodes) );
1776  else
1777  cutverts = nodearrint;
1778 
1779  if( edgearrint == NULL )
1780  SCIP_CALL( SCIPallocBufferArray(scip, &unsatarcs, nedges) );
1781  else
1782  unsatarcs = edgearrint;
1783 
1784  if( gnodearrterms == NULL )
1785  {
1786  SCIP_CALL( SCIPallocBufferArray(scip, &gnodearr, nterms - 1) );
1787  for( int i = 0; i < nterms - 1; i++ )
1788  SCIP_CALL( SCIPallocBlockMemory(scip, &gnodearr[i]) ); /*lint !e866*/
1789  }
1790  else
1791  {
1792  gnodearr = gnodearrterms;
1793  }
1794 
1795  SCIP_CALL( SCIPpqueueCreate(&pqueue, nterms, 2.0, GNODECmpByDist) );
1796 
1798  SCIP_CALL( SCIPallocMemoryArray(scip, &edgearr, nedges) );
1799  SCIP_CALL( SCIPallocMemoryArray(scip, &tailarr, nedges) );
1800  SCIP_CALL( SCIPallocMemoryArray(scip, &start, nnodes + 1) );
1801  SCIP_CALL( SCIPallocMemoryArray(scip, &gmark, nnodes + 1) );
1802  SCIP_CALL( SCIPallocMemoryArray(scip, &stackarr, nnodes) );
1803 
1804  /* fill auxiliary adjacent vertex/edges arrays */
1805  graph_get_csr(g, edgearr, tailarr, start, &ncsredges);
1806 
1807  /* initialize priority queue and res. capacity */
1808  SCIP_CALL( dualascent_init(scip, g, start, edgearr, root, is_pseudoroot, ncsredges, gmark, active, pqueue,
1809  gnodearr, rescap, &dualobj, &augmentingcomponent) );
1810 
1811  /* mark whether an arc is satisfied (has capacity 0) */
1812  for( int i = 0; i < ncsredges; i++ )
1813  {
1814 #ifdef BITFIELDSARRAY
1815  if( SCIPisZero(scip, rescap[i]) )
1816  SetBit(bitarr, i);
1817  else
1818  CleanBit(bitarr, i);
1819 #else
1820  if( rescap[i] == 0.0 )
1821  {
1822  if( active[tailarr[i] - 1] == 0 )
1823  tailarr[i] = 0;
1824  else
1825  tailarr[i] *= -1;
1826  }
1827 #endif
1828  }
1829 
1830  norgcutverts = 0;
1831 
1832  /* (main) dual ascent loop */
1833  while( SCIPpqueueNElems(pqueue) > 0 && !SCIPisStopped(scip) )
1834  {
1835  /* get active vertex of minimum score */
1836  GNODE* const gnodeact = (GNODE*) SCIPpqueueRemove(pqueue);
1837  const SCIP_Real prio1 = gnodeact->dist;
1838  const SCIP_Real prio2 = (SCIPpqueueNElems(pqueue) > 0) ? ((GNODE*) SCIPpqueueFirst(pqueue))->dist : FARAWAY;
1839  const int v = gnodeact->number;
1840  SCIP_Real degsum = g->grad[v];
1841  int ncutverts = 0;
1842  int nunsatarcs = 0;
1843 
1844  SCIP_Bool firstrun = TRUE;
1845 
1846  SCIPdebugMessage("DA: START WITH v %d prio1 %f prio2 %f \n", v, prio1, prio2);
1847 
1848  /* perform augmentation as long as priority of root component does not exceed max deviation */
1849  for( ; ; )
1850  {
1851  assert(stacklength == 0);
1852 
1853  /* 1. step: BFS from v (or connected component) on saturated, incoming arcs */
1854 
1855  if( firstrun )
1856  {
1857  firstrun = FALSE;
1858  gmark[v + 1] = TRUE;
1859  cutverts[ncutverts++] = v;
1860  assert(stacklength < nnodes);
1861  stackarr[stacklength++] = v;
1862  }
1863  /* not in first processing of root component: */
1864  else
1865  {
1866  for( int i = norgcutverts; i < ncutverts; i++ )
1867  {
1868  const int s = cutverts[i];
1869 
1870  assert(gmark[s + 1]);
1871  assert(active[s] != 0);
1872  assert(stacklength < nnodes);
1873 
1874  stackarr[stacklength++] = s;
1875  }
1876  }
1877 #ifdef DFS
1878  while( stacklength )
1879  {
1880  const int node = stackarr[--stacklength];
1881 #else
1882  for( int n = 0; n < stacklength; n++ )
1883  {
1884  int end;
1885 
1886  assert(n < nnodes);
1887  node = stackarr[n];
1888 #endif
1889 
1890  /* traverse incoming arcs */
1891  for( int i = start[node], end = start[node + 1]; i != end; i++ )
1892  {
1893  int tail = tailarr[i];
1894 
1895  /* zero reduced-cost arc? */
1896  if( tail <= 0 )
1897  {
1898  tail *= -1;
1899  if( !gmark[tail] )
1900  {
1901  /* if an active vertex has been hit (other than v), break */
1902  if( 0 == tail )
1903  {
1904  const int realtail = g->tail[edgearr[i]];
1905 
1906  /* v should not be processed */
1907  if( realtail == v )
1908  continue;
1909 
1910  /* is realtail active or does realtail lead to an active vertex other than v? */
1911  if( is_active(active, realtail, v) )
1912  {
1913  active[v] = realtail + 1;
1914  stacklength = 0;
1915  goto ENDOFLOOP;
1916  }
1917 
1918  tail = realtail + 1;
1919 
1920  /* have we processed tail already? */
1921  if( gmark[tail] )
1922  continue;
1923  }
1924 
1925  assert(tail > 0);
1926 
1927  gmark[tail] = TRUE;
1928  tail--;
1929  cutverts[ncutverts++] = tail;
1930  degsum += g->grad[tail];
1931 
1932  assert(stacklength < nnodes);
1933  stackarr[stacklength++] = tail;
1934  } /* marked */
1935  } /* zero reduced-cost arc */
1936  else if( !gmark[tail] )
1937  {
1938  unsattails[nunsatarcs] = tail;
1939  unsatarcs[nunsatarcs++] = i;
1940  }
1941  }
1942  }
1943 #ifndef DFS
1944  stacklength = 0;
1945 #endif
1946  currscore = degsum - (ncutverts - 1);
1947 
1948  /* guiding solution provided? */
1949  if( result != NULL )
1950  {
1951  int nsolarcs = 0;
1952  for( int i = 0; i < nunsatarcs; i++ )
1953  {
1954  const int a = unsatarcs[i];
1955 
1956  assert(tailarr[a] > 0);
1957 
1958  if( !(gmark[tailarr[a]]) )
1959  {
1960  if( result[edgearr[a]] == CONNECT )
1961  nsolarcs++;
1962  }
1963  }
1964 
1965  assert(nsolarcs > 0);
1966  assert(currscore <= nedges);
1967 
1968  if( nsolarcs > 1 )
1969  currscore += (SCIP_Real) ((nsolarcs - 1) * (g->knots * 2.0));
1970  }
1971  else
1972  {
1973  assert(SCIPisGE(scip, currscore, prio1));
1974  }
1975 
1976  SCIPdebugMessage("DA: deviation %f \n", (currscore - prio1) / prio1);
1977  SCIPdebugMessage("DA: currscore %f prio1 %f prio2 %f \n", currscore, prio1, prio2);
1978 
1979  /* augmentation criteria met? */
1980  if( ((currscore - prio1) / prio1) <= maxdeviation || currscore <= prio2 )
1981  {
1982  SCIP_CONS* cons = NULL;
1983  SCIP_ROW* row = NULL;
1984 
1985  int shift = 0;
1986  SCIP_Real min = FARAWAY;
1987  SCIP_Bool isactive = FALSE;
1988 
1989  /* 2. step: get minimum residual capacity among cut-arcs */
1990 
1991  /* adjust array of unsatisfied arcs */
1992 
1993  for( int i = 0; i < nunsatarcs; i++ )
1994  {
1995  const int tail = unsattails[i];
1996 
1997  if( gmark[tail] )
1998  {
1999  shift++;
2000  }
2001  else
2002  {
2003  const int a = unsatarcs[i];
2004 
2005  assert(tailarr[a] > 0);
2006  assert(rescap[a] > 0);
2007 
2008  if( rescap[a] < min )
2009  min = rescap[a];
2010  if( shift )
2011  {
2012  unsattails[i - shift] = tail;
2013  unsatarcs[i - shift] = a;
2014  }
2015  }
2016  }
2017 
2018  assert(SCIPisLT(scip, min, FARAWAY));
2019  nunsatarcs -= shift;
2020 
2021  norgcutverts = ncutverts;
2022 
2023  /* 3. step: perform augmentation */
2024 
2025  /* create constraints/cuts ? */
2026  if( addcuts )
2027  {
2028  if( addconss )
2029  {
2030  SCIP_CALL( SCIPcreateConsLinear(scip, &cons, "da", 0, NULL, NULL,
2031  1.0, SCIPinfinity(scip), TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, TRUE, FALSE) );
2032  }
2033  else
2034  {
2035  SCIP_CALL( SCIPcreateEmptyRowCons(scip, &row, conshdlr, "da", 1.0,
2036  SCIPinfinity(scip), FALSE, FALSE, TRUE) );
2037 
2038  SCIP_CALL( SCIPcacheRowExtensions(scip, row) );
2039  }
2040  }
2041 
2042  shift = 0;
2043 
2044  /* update (dual) objective */
2045  dualobj += min;
2046 
2047  for( int i = 0; i < nunsatarcs; i++ )
2048  {
2049  const int a = unsatarcs[i];
2050  assert(a >= 0);
2051 
2052  if( addcuts )
2053  {
2054  assert(vars != NULL);
2055 
2056  if( addconss )
2057  SCIP_CALL( SCIPaddCoefLinear(scip, cons, vars[edgearr[a]], 1.0) );
2058  else
2059  SCIP_CALL( SCIPaddVarToRow(scip, row, vars[edgearr[a]], 1.0) );
2060  }
2061  rescap[a] -= min;
2062 
2063  assert(SCIPisGE(scip, rescap[a], 0.0));
2064 
2065  if( rescap[a] <= DA_EPS )
2066  {
2067  int tail = unsattails[i];
2068 
2069  rescap[a] = 0.0;
2070 
2071  assert(tail > 0);
2072  assert(tailarr[a] > 0);
2073 
2074  tailarr[a] *= -1;
2075 
2076  if( active[tail - 1] >= 0 && is_active(active, tail - 1, v) )
2077  {
2078  assert(tail - 1 != v);
2079  tailarr[a] = 0;
2080  if( !isactive )
2081  {
2082  isactive = TRUE;
2083  active[v] = tail;
2084  }
2085  }
2086 
2087 
2088  if( !(gmark[tail]) )
2089  {
2090  assert(tail != 0);
2091 
2092  gmark[tail] = TRUE;
2093  tail--;
2094  degsum += g->grad[tail];
2095  cutverts[ncutverts++] = tail;
2096  }
2097 
2098  shift++;
2099  }
2100  else if( shift )
2101  {
2102  unsattails[i - shift] = unsattails[i];
2103  unsatarcs[i - shift] = a;
2104  }
2105  }
2106 
2107  if( addcuts )
2108  {
2109  if( addconss )
2110  {
2111  SCIP_CALL( SCIPaddCons(scip, cons) );
2112  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
2113  }
2114  else
2115  {
2116  SCIP_Bool infeasible;
2117 
2118  SCIP_CALL( SCIPflushRowExtensions(scip, row) );
2119  SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
2120  SCIP_CALL( SCIPreleaseRow(scip, &row) );
2121 
2122  assert(!infeasible);
2123  }
2124  }
2125 
2126  if( isactive )
2127  {
2128  stacklength = 0;
2129  goto ENDOFLOOP;
2130  }
2131  nunsatarcs -= shift;
2132 
2133  continue;
2134  }
2135  else
2136  {
2137  SCIP_Bool insert = TRUE;
2138 
2139  if( is_pseudoroot )
2140  {
2141  int i = start[v];
2142  const int end = start[v + 1];
2143 
2144  assert(end - i == 2);
2145 
2146  for( ; i != end; i++ )
2147  if( rescap[i] != 0.0 )
2148  break;
2149 
2150  if( i == end )
2151  {
2152  if( augmentingcomponent == -1 )
2153  augmentingcomponent = v;
2154 
2155  if( augmentingcomponent != v )
2156  insert = FALSE;
2157  }
2158  }
2159 
2160  if( insert )
2161  {
2162  /* reinsert active vertex */
2163  gnodeact->dist = currscore;
2164  SCIP_CALL( SCIPpqueueInsert(pqueue, gnodeact) );
2165  }
2166  }
2167 
2168  ENDOFLOOP:
2169 
2170  for( int i = 0; i < ncutverts; i++ )
2171  gmark[cutverts[i] + 1] = FALSE;
2172 
2173  for( int i = 0; i < nnodes + 1; i++ )
2174  {
2175  assert(!gmark[i]);
2176  }
2177 
2178  break;
2179  } /* augmentation loop */
2180  } /* dual ascent loop */
2181 
2182  SCIPdebugMessage("DA: dualglobal: %f \n", dualobj);
2183  *objval = dualobj;
2184 
2185  for( int i = ncsredges; i < nedges; i++ )
2186  {
2187  edgearr[i] = i;
2188  rescap[i] = g->cost[i];
2189  }
2190 
2191  /* re-extend rescap array */
2192  for( int i = 0; i < ncsredges; i++ )
2193  {
2194  if( edgearr[i] != i )
2195  {
2196  SCIP_Real bufferedval = rescap[i];
2197  int a = i;
2198 
2199  rescap[i] = g->cost[i];
2200  while( edgearr[a] != a )
2201  {
2202  const int shift = edgearr[a];
2203  const SCIP_Real min = rescap[shift];
2204 
2205  rescap[shift] = bufferedval;
2206  bufferedval = min;
2207  edgearr[a] = a;
2208  a = shift;
2209  }
2210  }
2211  }
2212 
2213 #ifdef BITFIELDSARRAY
2214  SCIPfreeBufferArray(scip, &bitarr);
2215 #endif
2216 
2217  SCIPfreeMemoryArray(scip, &stackarr);
2218  SCIPfreeMemoryArray(scip, &gmark);
2219  SCIPfreeMemoryArray(scip, &start);
2220  SCIPfreeMemoryArray(scip, &tailarr);
2221  SCIPfreeMemoryArray(scip, &edgearr);
2222  SCIPfreeMemoryArray(scip, &active);
2223 
2224  SCIPpqueueFree(&pqueue);
2225 
2226  if( gnodearrterms == NULL )
2227  {
2228  for( int i = nterms - 2; i >= 0; i-- )
2229  SCIPfreeBlockMemory(scip, &gnodearr[i]);
2230  SCIPfreeBufferArray(scip, &gnodearr);
2231  }
2232 
2233  /* call Ascend-And-Prune? */
2234  if( ascendandprune )
2235  {
2236  SCIP_Bool success;
2237  STP_Bool* RESTRICT mynodearrchar = NULL;
2238 
2239  if( nodearrchar == NULL )
2240  SCIP_CALL( SCIPallocBufferArray(scip, &mynodearrchar, nnodes) );
2241  else
2242  mynodearrchar = nodearrchar;
2243 
2244  SCIP_CALL( SCIPStpHeurAscendPruneRun(scip, NULL, g, rescap, unsatarcs, cutverts, root, mynodearrchar, &success, TRUE) );
2245 
2246  if( nodearrchar == NULL )
2247  SCIPfreeBufferArray(scip, &mynodearrchar);
2248  }
2249 
2250  if( edgearrint == NULL )
2251  SCIPfreeBufferArray(scip, &unsatarcs);
2252 
2253  if( nodearrint == NULL )
2254  SCIPfreeBufferArray(scip, &cutverts);
2255 
2256  if( redcost == NULL )
2257  SCIPfreeBufferArray(scip, &rescap);
2258 
2259  SCIPfreeBufferArray(scip, &unsattails);
2260 
2261  return SCIP_OKAY;
2262 }
2263 
2264 /** dual ascent heuristic for PCSPG and MWCSP */
2266  SCIP* scip, /**< SCIP data structure */
2267  GRAPH* g, /**< graph data structure */
2268  SCIP_Real* redcost, /**< array to store reduced costs or NULL */
2269  SCIP_Real* objval, /**< pointer to store objective value */
2270  SCIP_Bool addcuts, /**< should dual-ascent add Steiner cuts? */
2271  SCIP_Bool ascendandprune, /**< perform ascend-and-prune and add solution? */
2272  int nruns /**< number of dual ascent runs */
2273  )
2274 {
2275  SCIP_CONSHDLR* conshdlr = NULL;
2276  SCIP_PQUEUE* pqueue;
2277  SCIP_VAR** vars;
2278  GRAPH* transgraph;
2279  SCIP_Real min;
2280  SCIP_Real prio1;
2281  SCIP_Real offset;
2282  SCIP_Real dualobj;
2283  SCIP_Real currscore;
2284  SCIP_Real maxdeviation;
2285  SCIP_Real* rescap;
2286  GNODE* gnodeact;
2287  GNODE** gnodearr;
2288  int s;
2289  int i;
2290  int k;
2291  int v;
2292  int a;
2293  int tail;
2294  int pnode;
2295  int shift;
2296  int root;
2297  int nnodes;
2298  int nterms;
2299  int nedges;
2300  int degsum;
2301  int ncutverts;
2302  int pseudoroot;
2303  int nunsatarcs;
2304  int stacklength;
2305  int norgcutverts;
2306  int* cutverts;
2307  int* stackarr;
2308  STP_Bool* origedge;
2309  int* unsatarcs;
2310  STP_Bool firstrun;
2311  STP_Bool* sat;
2312  STP_Bool* active;
2313  const SCIP_Bool addconss = (SCIPgetStage(scip) < SCIP_STAGE_INITSOLVE);
2314 
2315  /* should currently not be activated */
2316  assert(addconss || !addcuts);
2317 
2318  assert(g != NULL);
2319  assert(scip != NULL);
2320  assert(nruns >= 0);
2321  assert(objval != NULL);
2322 
2323  if( g->knots == 1 )
2324  return SCIP_OKAY;
2325 
2326  if( addcuts )
2327  {
2328  vars = SCIPprobdataGetVars(scip);
2329  assert(vars != NULL);
2330  if( !addconss )
2331  {
2332  conshdlr = SCIPfindConshdlr(scip, "stp");
2333  assert(conshdlr != NULL);
2334  }
2335  }
2336  else
2337  {
2338  vars = NULL;
2339  }
2340 
2341  root = g->source;
2342  degsum = 0;
2343  offset = 0.0;
2344  dualobj = 0.0;
2345 
2346  ncutverts = 0;
2347  norgcutverts = 0;
2348  maxdeviation = DEFAULT_DAMAXDEVIATION;
2349 
2350  SCIP_CALL( graph_pc_getSap(scip, g, &transgraph, &offset) );
2351 
2352  nnodes = transgraph->knots;
2353  nedges = transgraph->edges;
2354  nterms = transgraph->terms;
2355  pseudoroot = nnodes - 1;
2356 
2357  if( redcost == NULL )
2358  {
2359  SCIP_CALL( SCIPallocBufferArray(scip, &rescap, nedges) );
2360  }
2361  else
2362  {
2363  rescap = redcost;
2364  }
2365 
2366  stacklength = 0;
2367  SCIP_CALL( SCIPallocBufferArray(scip, &stackarr, nnodes) );
2368  SCIP_CALL( SCIPallocBufferArray(scip, &sat, nedges) );
2370  SCIP_CALL( SCIPallocBufferArray(scip, &cutverts, nnodes) );
2371  SCIP_CALL( SCIPallocBufferArray(scip, &gnodearr, nterms - 1) );
2372  SCIP_CALL( SCIPallocBufferArray(scip, &unsatarcs, nedges) );
2373  SCIP_CALL( SCIPallocBufferArray(scip, &origedge, nedges) );
2374 
2375  for( i = 0; i < nedges; i++ )
2376  if( !Is_term(transgraph->term[transgraph->tail[i]]) && transgraph->head[i] == pseudoroot )
2377  origedge[i] = FALSE;
2378  else if( transgraph->tail[i] == pseudoroot && !Is_term(transgraph->term[transgraph->head[i]]) )
2379  origedge[i] = FALSE;
2380  else
2381  origedge[i] = TRUE;
2382 
2383  for( i = 0; i < nterms - 1; i++ )
2384  {
2385  SCIP_CALL( SCIPallocBuffer(scip, &gnodearr[i]) ); /*lint !e866*/
2386  }
2387 
2388  SCIP_CALL( SCIPpqueueCreate( &pqueue, nnodes, 2.0, GNODECmpByDist) );
2389 
2390  k = 0;
2391  /* mark terminals as active, add all except root to pqueue */
2392  for( i = 0; i < nnodes; i++ )
2393  {
2394  if( Is_term(transgraph->term[i]) )
2395  {
2396  active[i] = TRUE;
2397  assert(transgraph->grad[i] > 0);
2398  if( i != root )
2399  {
2400  gnodearr[k]->number = i;
2401  gnodearr[k]->dist = transgraph->grad[i];
2402 
2403  for( a = transgraph->inpbeg[i]; a != EAT_LAST; a = transgraph->ieat[a] )
2404  if( SCIPisEQ(scip, transgraph->cost[a], 0.0) )
2405  break;
2406 
2407  if( a != EAT_LAST )
2408  gnodearr[k]->dist += transgraph->grad[transgraph->tail[a]] - 1;
2409 
2410  assert(gnodearr[k]->dist > 0);
2411 
2412  SCIP_CALL( SCIPpqueueInsert(pqueue, gnodearr[k++]) );
2413  }
2414  }
2415  else
2416  {
2417  active[i] = FALSE;
2418  }
2419  transgraph->mark[i] = FALSE;
2420  }
2421 
2422  for( i = 0; i < nedges; i++ )
2423  {
2424  rescap[i] = transgraph->cost[i];
2425  if( SCIPisZero(scip, rescap[i]) )
2426  sat[i] = TRUE;
2427  else
2428  sat[i] = FALSE;
2429  }
2430 
2431  /* dual ascent loop */
2432  while( SCIPpqueueNElems(pqueue) > 0 && !SCIPisStopped(scip) )
2433  {
2434  /* get active vertex of minimum score */
2435  gnodeact = (GNODE*) SCIPpqueueRemove(pqueue);
2436 
2437  v = gnodeact->number;
2438  prio1 = gnodeact->dist;
2439 
2440  firstrun = TRUE;
2441  nunsatarcs = 0;
2442 
2443  /* perform augmentation as long as ... */
2444  for( ; ; )
2445  {
2446  assert(stacklength == 0);
2447  /* 1. step: BFS from v (or connected component) on saturated, incoming arcs */
2448 
2449  if( firstrun )
2450  {
2451  degsum = transgraph->grad[v];
2452  ncutverts = 0;
2453  firstrun = FALSE;
2454  nunsatarcs = 0;
2455  transgraph->mark[v] = TRUE;
2456  cutverts[ncutverts++] = v;
2457  stackarr[stacklength++] = v;
2458  }
2459  /* not in first processing of root component: */
2460  else
2461  {
2462  for( i = norgcutverts; i < ncutverts; i++ )
2463  {
2464  s = cutverts[i];
2465  assert(transgraph->mark[s]);
2466  if( active[s] )
2467  {
2468  active[v] = FALSE;
2469  stacklength = 0;
2470  goto ENDOFLOOP;
2471  }
2472 
2473  stackarr[stacklength++] = s;
2474  }
2475  }
2476 
2477  while( stacklength )
2478  {
2479  pnode = stackarr[--stacklength];
2480 
2481  /* traverse incoming arcs */
2482  for( a = transgraph->inpbeg[pnode]; a != EAT_LAST; a = transgraph->ieat[a] )
2483  {
2484  tail = transgraph->tail[a];
2485  if( sat[a] )
2486  {
2487  if( !transgraph->mark[tail] )
2488  {
2489  /* if an active vertex has been hit, break */
2490  if( active[tail] )
2491  {
2492  active[v] = FALSE;
2493  stacklength = 0;
2494  goto ENDOFLOOP;
2495  }
2496 
2497  degsum += transgraph->grad[tail];
2498  transgraph->mark[tail] = TRUE;
2499  cutverts[ncutverts++] = tail;
2500  stackarr[stacklength++] = tail;
2501  }
2502  }
2503  else if( !transgraph->mark[tail] )
2504  {
2505  unsatarcs[nunsatarcs++] = a;
2506  }
2507  }
2508  }
2509 
2510  currscore = degsum - (ncutverts - 1);
2511 
2512  assert(SCIPisGE(scip, currscore, prio1));
2513 
2514  /* augmentation criteria met? */
2515  if( SCIPisLE(scip, (currscore - prio1) / prio1, maxdeviation) || (SCIPpqueueNElems(pqueue) == 0) )
2516  {
2517  SCIP_Bool in = FALSE;
2518  SCIP_ROW* row;
2519  SCIP_CONS* cons = NULL;
2520 
2521  /* 2. pass: get minimum residual capacity among cut-arcs */
2522 
2523  /* adjust array of unsatisfied arcs */
2524  min = FARAWAY;
2525  shift = 0;
2526 
2527  for( i = 0; i < nunsatarcs; i++ )
2528  {
2529  a = unsatarcs[i];
2530  if( transgraph->mark[transgraph->tail[a]] )
2531  {
2532  shift++;
2533  }
2534  else
2535  {
2536 
2537  assert(!sat[a]);
2538  if( SCIPisLT(scip, rescap[a], min) )
2539  min = rescap[a];
2540  if( shift != 0 )
2541  unsatarcs[i - shift] = a;
2542  }
2543  }
2544 
2545  assert(SCIPisLT(scip, min, FARAWAY));
2546  nunsatarcs -= shift;
2547 
2548  if( nunsatarcs > 0)
2549  assert(!transgraph->mark[transgraph->tail[unsatarcs[nunsatarcs-1]]]);
2550 
2551  norgcutverts = ncutverts;
2552 
2553 
2554  /* 3. pass: perform augmentation */
2555 
2556 
2557  /* create constraint/row */
2558 
2559  if( addcuts )
2560  {
2561  if( addconss )
2562  {
2563  SCIP_CALL( SCIPcreateConsLinear(scip, &cons, "da", 0, NULL, NULL,
2564  1.0, SCIPinfinity(scip), TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, TRUE, TRUE, FALSE) );
2565  }
2566  else
2567  {
2568  SCIP_CALL(SCIPcreateEmptyRowCons(scip, &row, conshdlr, "da", 1.0, SCIPinfinity(scip), FALSE, FALSE, TRUE));
2569  SCIP_CALL(SCIPcacheRowExtensions(scip, row));
2570  }
2571  }
2572 
2573  dualobj += min;
2574  for( i = 0; i < nunsatarcs; i++ )
2575  {
2576  a = unsatarcs[i];
2577  if( a == -1 )
2578  continue;
2579 
2580  if( addcuts && origedge[a] )
2581  {
2582  assert(vars != NULL);
2583  assert(cons != NULL);
2584 
2585  if( g->tail[a] == root && g->head[a] == v )
2586  in = TRUE;
2587 
2588  if( addconss )
2589  SCIP_CALL( SCIPaddCoefLinear(scip, cons, vars[a], 1.0) );
2590  else
2591  SCIP_CALL( SCIPaddVarToRow(scip, row, vars[a], 1.0) );
2592  }
2593  rescap[a] -= min;
2594 
2595  assert(SCIPisGE(scip, rescap[a], 0.0));
2596 
2597  if( SCIPisEQ(scip, rescap[a], 0.0) )
2598  {
2599  sat[a] = TRUE;
2600  if( !(transgraph->mark[transgraph->tail[a]]) )
2601  {
2602  tail = transgraph->tail[a];
2603  degsum += transgraph->grad[tail];
2604  transgraph->mark[tail] = TRUE;
2605  cutverts[ncutverts++] = tail;
2606  }
2607  }
2608  }
2609 
2610  if( addcuts )
2611  {
2612  assert(vars != NULL);
2613 
2614  if( !in )
2615  {
2616  for( i = g->outbeg[root]; i != EAT_LAST; i = g->oeat[i] )
2617  if( g->head[i] == v )
2618  {
2619  if( addconss )
2620  SCIP_CALL( SCIPaddCoefLinear(scip, cons, vars[i], 1.0) );
2621  else
2622  SCIP_CALL( SCIPaddVarToRow(scip, row, vars[i], 1.0) );
2623  }
2624  }
2625 
2626  if( addconss )
2627  {
2628  SCIP_CALL( SCIPaddCons(scip, cons) );
2629  SCIP_CALL( SCIPreleaseCons(scip, &cons) );
2630  }
2631  else
2632  {
2633  SCIP_Bool infeasible;
2634  assert(row != NULL);
2635 
2636  SCIP_CALL( SCIPflushRowExtensions(scip, row) );
2637  SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
2638  SCIP_CALL( SCIPreleaseRow(scip, &row) );
2639 
2640  assert(!infeasible);
2641  }
2642  }
2643 
2644  continue;
2645  }
2646  else
2647  {
2648  /* reinsert active vertex */
2649  gnodeact->dist = currscore;
2650  SCIP_CALL( SCIPpqueueInsert(pqueue, gnodeact) );
2651  }
2652 
2653  ENDOFLOOP:
2654 
2655  for( i = 0; i < ncutverts; i++ )
2656  transgraph->mark[cutverts[i]] = FALSE;
2657 
2658  break;
2659  } /* augmentation loop */
2660  } /* dual ascent loop */
2661 
2662 
2663  *objval = dualobj + offset;
2664  SCIPdebugMessage("DA: dualglobal: %f \n", *objval + SCIPprobdataGetOffset(scip));
2665 
2666  /* call dual Ascend-And-Prune? */
2667  if( ascendandprune )
2668  {
2669  SCIP_Bool success;
2670  SCIP_CALL( SCIPStpHeurAscendPruneRun(scip, NULL, g, rescap, unsatarcs, cutverts, -1, active, &success, TRUE));
2671  }
2672 
2673  /* free memory */
2674  SCIPpqueueFree(&pqueue);
2675 
2676  for( i = nterms - 2; i >= 0; i-- )
2677  SCIPfreeBuffer(scip, &gnodearr[i]);
2678 
2679  SCIPfreeBufferArray(scip, &origedge);
2680  SCIPfreeBufferArray(scip, &unsatarcs);
2681  SCIPfreeBufferArray(scip, &cutverts);
2682  SCIPfreeBufferArray(scip, &gnodearr);
2683  SCIPfreeBufferArray(scip, &active);
2684  SCIPfreeBufferArray(scip, &sat);
2685  SCIPfreeBufferArray(scip, &stackarr);
2686 
2687  if( redcost == NULL )
2688  SCIPfreeBufferArray(scip, &rescap);
2689 
2690  graph_free(scip, &transgraph, TRUE);
2691 
2692  return SCIP_OKAY;
2693 
2694 }
2695 
2696 /**@} */
void SCIPconshdlrSetData(SCIP_CONSHDLR *conshdlr, SCIP_CONSHDLRDATA *conshdlrdata)
Definition: cons.c:4221
int SCIPpqueueNElems(SCIP_PQUEUE *pqueue)
Definition: misc.c:1263
static int graph_next_term(const GRAPH *g, int terms, int *term, const int *w, const SCIP_Bool firstrun)
Definition: cons_stp.c:277
#define CONSHDLR_PROPFREQ
Definition: cons_stp.c:76
SCIP_RETCODE SCIPsetConshdlrDelete(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_DECL_CONSDELETE((*consdelete)))
Definition: scip_cons.c:640
static volatile int nterms
Definition: interrupt.c:37
int *RESTRICT mincut_e
Definition: grph.h:113
#define DEFAULT_CREEPFLOW
Definition: cons_stp.c:93
#define NULL
Definition: def.h:239
static SCIP_DECL_CONSINITSOL(consInitsolStp)
Definition: cons_stp.c:1245
SCIP_RETCODE SCIPcacheRowExtensions(SCIP *scip, SCIP_ROW *row)
Definition: scip_lp.c:1547
#define DEFAULT_NESTEDCUT
Definition: cons_stp.c:95
int *RESTRICT head
Definition: grph.h:96
SCIP_Bool SCIPisFeasEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
SCIP_NODE * SCIPgetCurrentNode(SCIP *scip)
Definition: scip_tree.c:158
SCIP_STAGE SCIPgetStage(SCIP *scip)
Definition: scip_general.c:412
Definition: grph.h:57
SCIP_Bool SCIPconsIsDynamic(SCIP_CONS *cons)
Definition: cons.c:8335
SCIP_RETCODE SCIPsetConshdlrTrans(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_DECL_CONSTRANS((*constrans)))
Definition: scip_cons.c:663
int source
Definition: grph.h:67
SCIP_Bool SCIPisFeasLT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
static SCIP_DECL_CONSENFOLP(consEnfolpStp)
Definition: cons_stp.c:1401
static SCIP_DECL_CONSSEPALP(consSepalpStp)
Definition: cons_stp.c:1323
SCIP_CONSHDLR * SCIPfindConshdlr(SCIP *scip, const char *name)
Definition: scip_cons.c:954
SCIP_RETCODE SCIPflushRowExtensions(SCIP *scip, SCIP_ROW *row)
Definition: scip_lp.c:1570
Constraint handler for Steiner problems.
#define SCIPfreeMemoryArray(scip, ptr)
Definition: scip_mem.h:88
#define DEFAULT_BACKCUT
Definition: cons_stp.c:92
int terms
Definition: grph.h:64
SCIP_RETCODE SCIPaddVarToRow(SCIP *scip, SCIP_ROW *row, SCIP_VAR *var, SCIP_Real val)
Definition: scip_lp.c:1602
SCIP_VAR ** SCIPprobdataGetVars(SCIP *scip)
#define SCIPallocMemoryArray(scip, ptr, num)
Definition: scip_mem.h:72
SCIP_Bool SCIPisGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
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:2266
int *RESTRICT maxdeg
Definition: grph.h:78
#define EAT_LAST
Definition: grph.h:31
#define RESTRICT
Definition: def.h:251
SCIP_Bool SCIPisFeasNegative(SCIP *scip, SCIP_Real val)
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_Bool SCIPisFeasGE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
reduction and dual-cost based primal heuristic for Steiner problems
SCIP_CONS ** SCIPconshdlrGetConss(SCIP_CONSHDLR *conshdlr)
Definition: cons.c:4563
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:65
int *RESTRICT inpbeg
Definition: grph.h:74
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:243
#define FLOW_FACTOR
Definition: cons_stp.c:108
void SCIPStpBranchruleInitNodeState(const GRAPH *g, int *nodestate)
Definition: branch_stp.c:683
SCIP_Real SCIPinfinity(SCIP *scip)
void graph_free(SCIP *, GRAPH **, SCIP_Bool)
Definition: grphbase.c:3674
Problem data for stp problem.
#define TRUE
Definition: def.h:64
#define SCIPdebug(x)
Definition: pub_message.h:74
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:53
SCIP_Bool SCIPconsIsStickingAtNode(SCIP_CONS *cons)
Definition: cons.c:8355
void * SCIPpqueueFirst(SCIP_PQUEUE *pqueue)
Definition: misc.c:1249
int SCIPprobdataGetNVars(SCIP *scip)
static GRAPHNODE ** active
#define SCIPfreeBlockMemory(scip, ptr)
Definition: scip_mem.h:114
void SCIPpqueueFree(SCIP_PQUEUE **pqueue)
Definition: misc.c:1160
#define SCIPdebugMessage
Definition: pub_message.h:77
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:297
#define DEFAULT_MAXROUNDS
Definition: cons_stp.c:84
SCIP_Bool SCIPisEQ(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
int SCIPnodeGetDepth(SCIP_NODE *node)
Definition: tree.c:7354
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip_mem.h:142
int number
Definition: misc_stp.h:44
#define SCIPallocBlockMemory(scip, ptr)
Definition: scip_mem.h:97
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
SCIP_Bool SCIPconsIsRemovable(SCIP_CONS *cons)
Definition: cons.c:8345
static SCIP_DECL_CONSLOCK(consLockStp)
Definition: cons_stp.c:1522
SCIP_RETCODE SCIPsetConshdlrInitlp(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_DECL_CONSINITLP((*consinitlp)))
Definition: scip_cons.c:686
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:155
void graph_knot_chg(GRAPH *, int, int)
Definition: grphbase.c:2218
GRAPH * SCIPprobdataGetGraph(SCIP_PROBDATA *probdata)
SCIP_RETCODE SCIPaddCoefLinear(SCIP *scip, SCIP_CONS *cons, SCIP_VAR *var, SCIP_Real val)
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:1685
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:1011
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)
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_RETCODE SCIPaddVarLocksType(SCIP *scip, SCIP_VAR *var, SCIP_LOCKTYPE locktype, int nlocksdown, int nlocksup)
Definition: scip_var.c:4199
#define DA_EPS
Definition: cons_stp.c:101
SCIP_RETCODE SCIPsetConshdlrInitsol(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_DECL_CONSINITSOL((*consinitsol)))
Definition: scip_cons.c:506
SCIP_VAR * w
Definition: circlepacking.c:58
int *RESTRICT oeat
Definition: grph.h:103
SCIP_Bool SCIPisCutEfficacious(SCIP *scip, SCIP_SOL *sol, SCIP_ROW *cut)
Definition: scip_cut.c:161
#define CONNECT
Definition: grph.h:154
int *RESTRICT mincut_dist
Definition: grph.h:106
SCIP_RETCODE SCIPsetConshdlrCopy(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_DECL_CONSHDLRCOPY((*conshdlrcopy)), SCIP_DECL_CONSCOPY((*conscopy)))
Definition: scip_cons.c:409
#define CONSHDLR_DESC
Definition: cons_stp.c:71
int stp_type
Definition: grph.h:127
#define SCIPerrorMessage
Definition: pub_message.h:45
const char * SCIPconshdlrGetName(SCIP_CONSHDLR *conshdlr)
Definition: cons.c:4191
SCIP_RETCODE SCIPaddCons(SCIP *scip, SCIP_CONS *cons)
Definition: scip_prob.c:2822
SCIP_Bool SCIPisLT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
unsigned char STP_Bool
Definition: grph.h:52
SCIP_Real SCIPprobdataGetOffset(SCIP *scip)
#define CONSHDLR_SEPAPRIORITY
Definition: cons_stp.c:72
#define SCIPallocBuffer(scip, ptr)
Definition: scip_mem.h:128
#define STP_DCSTP
Definition: grph.h:43
#define SCIPfreeBufferArrayNull(scip, ptr)
Definition: scip_mem.h:143
SCIP_Real dist
Definition: misc_stp.h:45
SCIP_RETCODE SCIPincludeConshdlrStp(SCIP *scip)
Definition: cons_stp.c:1573
const char * SCIPconsGetName(SCIP_CONS *cons)
Definition: cons.c:8076
#define DEFAULT_MAXROUNDSROOT
Definition: cons_stp.c:85
int *RESTRICT grad
Definition: grph.h:73
SCIP_Bool SCIPconsIsPropagated(SCIP_CONS *cons)
Definition: cons.c:8295
#define DA_MAXDEVIATION_UPPER
Definition: cons_stp.c:100
SCIP_RETCODE SCIPsetConshdlrFree(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_DECL_CONSFREE((*consfree)))
Definition: scip_cons.c:434
SCIP_CONSHDLRDATA * SCIPconshdlrGetData(SCIP_CONSHDLR *conshdlr)
Definition: cons.c:4211
int *RESTRICT mincut_head_inact
Definition: grph.h:108
internal miscellaneous methods
int knots
Definition: grph.h:62
#define SCIP_CALL(x)
Definition: def.h:351
SCIP_Bool SCIPisFeasGT(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
#define CONSHDLR_DELAYPROP
Definition: cons_stp.c:81
SCIP_Bool SCIPconsIsLocal(SCIP_CONS *cons)
Definition: cons.c:8315
SCIP_RETCODE SCIPaddRow(SCIP *scip, SCIP_ROW *row, SCIP_Bool forcecut, SCIP_Bool *infeasible)
Definition: scip_cut.c:294
struct SCIP_ConsData SCIP_CONSDATA
Definition: type_cons.h:51
void * SCIPpqueueRemove(SCIP_PQUEUE *pqueue)
Definition: misc.c:1208
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 SCIPconshdlrGetNConss(SCIP_CONSHDLR *conshdlr)
Definition: cons.c:4593
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
#define STP_SPG
Definition: grph.h:38
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:130
#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:1542
#define SCIP_Bool
Definition: def.h:62
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 SCIPcreateEmptyRowCons(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:1336
int *RESTRICT tail
Definition: grph.h:95
int SCIPgetDepth(SCIP *scip)
Definition: scip_tree.c:715
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_RETCODE SCIPaddPoolCut(SCIP *scip, SCIP_ROW *row)
Definition: scip_cut.c:405
SCIP_Bool SCIPconsIsChecked(SCIP_CONS *cons)
Definition: cons.c:8275
SCIP_Bool SCIPconsIsInitial(SCIP_CONS *cons)
Definition: cons.c:8245
SCIP_RETCODE SCIPStpValidateSol(SCIP *, const GRAPH *, const double *, SCIP_Bool *)
Definition: validate.c:136
int *RESTRICT term
Definition: grph.h:68
#define CONSHDLR_EAGERFREQ
Definition: cons_stp.c:77
#define BMScopyMemoryArray(ptr, source, num)
Definition: memory.h:116
#define CREEP_VALUE
Definition: cons_stp.c:109
Constraint handler for linear constraints in their most general form, .
#define DA_MAXDEVIATION_LOWER
Definition: cons_stp.c:99
#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:1634
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:86
Steiner vertex branching rule.
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)
SCIP_RETCODE SCIPreleaseRow(SCIP *scip, SCIP_ROW **row)
Definition: scip_lp.c:1474
#define SCIPfreeBuffer(scip, ptr)
Definition: scip_mem.h:140
#define DEFAULT_MAXSEPACUTSROOT
Definition: cons_stp.c:87
#define Is_term(a)
Definition: grph.h:168
#define MAX(x, y)
Definition: def.h:208
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
SCIP_RETCODE SCIPpqueueInsert(SCIP_PQUEUE *pqueue, void *elem)
Definition: misc.c:1181
SCIP_CONSDATA * SCIPconsGetData(SCIP_CONS *cons)
Definition: cons.c:8106
SCIP_PROBDATA * SCIPgetProbData(SCIP *scip)
Definition: scip_prob.c:1020
#define CONSHDLR_DELAYSEPA
Definition: cons_stp.c:80
SCIP_RETCODE SCIPpqueueCreate(SCIP_PQUEUE **pqueue, int initsize, SCIP_Real sizefac, SCIP_DECL_SORTPTRCOMP((*ptrcomp)))
Definition: misc.c:1135
SCIP_Real * cost
Definition: grph.h:94
#define CONSHDLR_SEPAFREQ
Definition: cons_stp.c:75
SCIP_RETCODE SCIPreleaseCons(SCIP *scip, SCIP_CONS **cons)
Definition: scip_cons.c:1187
SCIP_VAR * a
Definition: circlepacking.c:57
#define SCIP_Real
Definition: def.h:150
SCIP_Bool SCIPconsIsModifiable(SCIP_CONS *cons)
Definition: cons.c:8325
SCIP_Bool SCIPisStopped(SCIP *scip)
Definition: scip_general.c:739
#define CONSHDLR_PROP_TIMING
Definition: cons_stp.c:90
SCIP_Bool SCIPconsIsEnforced(SCIP_CONS *cons)
Definition: cons.c:8265
SCIP_Real * SCIPprobdataGetXval(SCIP *scip, SCIP_SOL *sol)
int *RESTRICT outbeg
Definition: grph.h:76
SCIP_Bool SCIPconsIsSeparated(SCIP_CONS *cons)
Definition: cons.c:8255
SCIP_RETCODE SCIPprintRow(SCIP *scip, SCIP_ROW *row, FILE *file)
Definition: scip_lp.c:2094
#define CONSHDLR_ENFOPRIORITY
Definition: cons_stp.c:73
int edges
Definition: grph.h:92
SCIP_RETCODE graph_pc_getSap(SCIP *, GRAPH *, GRAPH **, SCIP_Real *)
Definition: grphbase.c:1075
#define flipedge(edge)
Definition: grph.h:150
#define Q_NULL
Definition: cons_stp.c:63
#define SCIPallocMemory(scip, ptr)
Definition: scip_mem.h:70
SCIP_Bool SCIPisZero(SCIP *scip, SCIP_Real val)
SCIP_Bool SCIPisLE(SCIP *scip, SCIP_Real val1, SCIP_Real val2)
struct SCIP_ConshdlrData SCIP_CONSHDLRDATA
Definition: type_cons.h:50
int *RESTRICT mincut_r
Definition: grph.h:115
SCIP_Real SCIPvarGetUbLocal(SCIP_VAR *var)
Definition: var.c:17409
#define nnodes
Definition: gastrans.c:65
#define CONSHDLR_NEEDSCONS
Definition: cons_stp.c:82
static SCIP_DECL_CONSCHECK(consCheckStp)
Definition: cons_stp.c:1451
GRAPH * SCIPprobdataGetGraph2(SCIP *scip)
static SCIP_DECL_CONSPROP(consPropStp)
Definition: cons_stp.c:1473
SCIP callable library.
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:129
SCIP_RETCODE SCIPsetConshdlrProp(SCIP *scip, SCIP_CONSHDLR *conshdlr, SCIP_DECL_CONSPROP((*consprop)), int propfreq, SCIP_Bool delayprop, SCIP_PROPTIMING proptiming)
Definition: scip_cons.c:343
void SCIPStpConshdlrSetGraph(SCIP *scip, const GRAPH *g)
Definition: cons_stp.c:1664