1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2021 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15
16 /**@file reduce_bnd.c
17 * @brief bound based reductions for Steiner tree problems
18 * @author Daniel Rehfeldt
19 *
20 * This file implements bound-based reduction techniques for several Steiner problems.
21 * All tests can be found in "A Generic Approach to Solving the Steiner Tree Problem and Variants" by Daniel Rehfeldt.
22 *
23 * A list of all interface methods can be found in grph.h.
24 *
25 */
26
27 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
28 #include <stdio.h>
29 #include <stdlib.h>
30 #include <string.h>
31 #include <assert.h>
32 #include "scip/scip.h"
33 #include "grph.h"
34 #include "heur_tm.h"
35 #include "heur_ascendprune.h"
36 #include "cons_stp.h"
37 #include "heur_local.h"
38 #include "misc_stp.h"
39 #include "prop_stp.h"
40 #include "probdata_stp.h"
41 #include "heur_rec.h"
42 #include "heur_slackprune.h"
43
44 #define DEFAULT_HEURRUNS 100 /**< number of runs of constructive heuristic */
45 #define DEFAULT_DARUNS 7 /**< number of runs for dual ascent heuristic */
46 #define DEFAULT_NMAXROOTS 8 /**< max number of roots to use for new graph in dual ascent heuristic */
47 #define PERTUBATION_RATIO 0.05 /**< pertubation ratio for dual-ascent primal bound computation */
48 #define PERTUBATION_RATIO_PC 0.005 /**< pertubation ratio for dual-ascent primal bound computation */
49 #define SOLPOOL_SIZE 20 /**< size of presolving solution pool */
50 #define STP_RED_MINBNDTERMS 750
51 #define STP_DABD_MAXDEGREE 5
52 #define STP_DABD_MAXDNEDGES 10
53 #define STP_DAEX_MAXDFSDEPTH 6
54 #define STP_DAEX_MINDFSDEPTH 4
55 #define STP_DAEX_MAXGRAD 8
56 #define STP_DAEX_MINGRAD 6
57 #define STP_DAEX_EDGELIMIT 50000
58 #define DAMAXDEVIATION_RANDOM_LOWER 0.15 /**< random upper bound for max deviation for dual ascent */
59 #define DAMAXDEVIATION_RANDOM_UPPER 0.30 /**< random upper bound for max deviation for dual ascent */
60 #define DAMAXDEVIATION_FAST 0.75
61
62 #define EXEDGE_FREE 0
63 #define EXEDGE_FIXED 1
64 #define EXEDGE_KILLED 2
65
66
67 /** mark ancestors of given edge */
68 static
markAncestorsConflict(const GRAPH * graph,int edge,int * ancestormark)69 SCIP_Bool markAncestorsConflict(
70 const GRAPH* graph, /**< graph data structure */
71 int edge, /**< edge to use */
72 int* ancestormark /**< ancestor mark array */
73 )
74 {
75 int count = 0;
76 assert(edge >= 0);
77
78 for( IDX* curr = graph->ancestors[edge]; curr != NULL && count <= STP_DAEX_MAXGRAD; curr = curr->parent, count++ )
79 {
80 const unsigned idx = ((unsigned) curr->index) / 2;
81
82 assert(curr->index >= 0 && idx < (unsigned) (MAX(graph->edges, graph->orgedges) / 2));
83
84 if( ancestormark[idx] )
85 return TRUE;
86
87 ancestormark[idx] = 1;
88 }
89
90 return FALSE;
91 }
92
93
94 /** mark ancestors of given edge */
95 static
markAncestors(const GRAPH * graph,int edge,int * ancestormark)96 void markAncestors(
97 const GRAPH* graph, /**< graph data structure */
98 int edge, /**< edge to use */
99 int* ancestormark /**< ancestor mark array */
100 )
101 {
102 int count = 0;
103 assert(edge >= 0);
104
105 for( IDX* curr = graph->ancestors[edge]; curr != NULL && count <= STP_DAEX_MAXGRAD; curr = curr->parent, count++ )
106 {
107 assert(curr->index >= 0 && curr->index / 2 < (MAX(graph->edges, graph->orgedges) / 2));
108 assert((ancestormark[((unsigned) curr->index) / 2]) == 0);
109
110 ancestormark[((unsigned) curr->index) / 2] = 1;
111 }
112 }
113
114 /** unmark ancestors of given edge */
115 static
unmarkAncestorsConflict(const GRAPH * graph,int edge,int * ancestormark)116 void unmarkAncestorsConflict(
117 const GRAPH* graph, /**< graph data structure */
118 int edge, /**< edge to use */
119 int* ancestormark /**< ancestor mark array */
120 )
121 {
122 int count = 0;
123 assert(edge >= 0);
124
125 for( IDX* curr = graph->ancestors[edge]; curr != NULL && count <= STP_DAEX_MAXGRAD; curr = curr->parent, count++ )
126 {
127 assert(curr->index >= 0 && curr->index / 2 < (MAX(graph->edges, graph->orgedges) / 2));
128 ancestormark[((unsigned) curr->index) / 2] = 0;
129 }
130 }
131
132 /** unmark ancestors of given edge */
133 static
unmarkAncestors(const GRAPH * graph,int edge,int * ancestormark)134 void unmarkAncestors(
135 const GRAPH* graph, /**< graph data structure */
136 int edge, /**< edge to use */
137 int* ancestormark /**< ancestor mark array */
138 )
139 {
140 int count = 0;
141 assert(edge >= 0);
142
143 for( IDX* curr = graph->ancestors[edge]; curr != NULL && count <= STP_DAEX_MAXGRAD; curr = curr->parent, count++ )
144 {
145 const unsigned idx = ((unsigned) curr->index) / 2;
146
147 assert(curr->index >= 0 && curr->index / 2 < (MAX(graph->edges, graph->orgedges) / 2));
148 assert(ancestormark[idx] == 1);
149
150 ancestormark[idx] = 0;
151 }
152 }
153
154
155 /** finalize subtree computations (clean up, update global bound) */
156 static
finalizeSubtree(const GRAPH * graph,const int * edgeends,const int * treeedges,int dfsdepth,SCIP_Bool stopped,SCIP_Real localbound,SCIP_Real * globalbound,SCIP_Bool * ruleout,int * nodepos,int * ancestormark)157 void finalizeSubtree(
158 const GRAPH* graph, /**< graph data structure */
159 const int* edgeends, /**< heads or tail of edge */
160 const int* treeedges, /**< tree edges */
161 int dfsdepth, /**< dfs depth */
162 SCIP_Bool stopped, /**< has extension been stopped? */
163 SCIP_Real localbound, /**< local bound */
164 SCIP_Real* globalbound, /**< points to global bound */
165 SCIP_Bool* ruleout, /**< rule out? */
166 int* nodepos, /**< node position in tree */
167 int* ancestormark /**< ancestor mark array */
168 )
169 {
170
171 if( localbound > *globalbound )
172 *globalbound = localbound;
173
174 if( !stopped )
175 *ruleout = TRUE;
176
177 if( !(*ruleout) )
178 {
179 for( int etree = 0; etree < dfsdepth; etree++ )
180 {
181 const int node = edgeends[treeedges[etree]];
182 assert(nodepos[node]);
183 nodepos[node] = 0;
184 unmarkAncestors(graph, treeedges[etree], ancestormark);
185 }
186 }
187 else
188 {
189 assert(dfsdepth == 0);
190 }
191 }
192
193 /** should we truncate subtree? */
194 static
truncateSubtree(const GRAPH * graph,SCIP_Real extendedcost,int root,int currhead,int maxgrad,int dfsdepth,int maxdfsdepth,SCIP_Real * minbound,SCIP_Bool * stopped)195 SCIP_Bool truncateSubtree(
196 const GRAPH* graph, /**< graph data structure */
197 SCIP_Real extendedcost, /**< cost of subtree */
198 int root, /**< DA root */
199 int currhead, /**< latest node */
200 int maxgrad, /**< maximum allowed degree */
201 int dfsdepth, /**< dfs depth */
202 int maxdfsdepth, /**< max dfs depth */
203 SCIP_Real* minbound, /**< bound */
204 SCIP_Bool* stopped /**< real truncation? */
205 )
206 {
207 if( Is_term(graph->term[currhead]) || graph->grad[currhead] > maxgrad || dfsdepth >= maxdfsdepth )
208 {
209 /* real truncation? */
210 if( root != currhead )
211 {
212 *stopped = TRUE;
213 if( extendedcost < *minbound )
214 *minbound = extendedcost;
215 }
216 return TRUE;
217 }
218 return FALSE;
219 }
220 #if 0
221
222 /** extend subtree and return new nadded_edges */
223 static
224 int reduceAndExtendSubtree(
225 SCIP* scip, /**< SCIP */
226 const GRAPH* graph, /**< graph data structure */
227 const int* graph_start, /**< start */
228 const int* graph_next, /**< next */
229 const int* graph_endp, /**< next */
230 int currnode, /**< latest node */
231 int nadded_edges, /**< added edges */
232 int* nodepos, /**< node position in tree */
233 int* edgestack, /**< stack */
234 unsigned char* extensionmark /**< mark array for extension */
235 )
236 {
237 int n = nadded_edges;
238 const int start = n;
239 SCIP_Bool cleanup = FALSE;
240
241 /* extend from currnode */
242 for( int e = graph_start[currnode]; e != EAT_LAST; e = graph_next[e] )
243 if( !nodepos[graph_endp[e]] )
244 {
245 /* check whether edge needs to be added */
246
247 const int vertex_exnew = graph_endp[e];
248
249 assert((n - start) <= STP_DAEX_MAXGRAD);
250 extensionmark[n - start] = EXEDGE_FREE;
251
252 for( int e2 = graph->outbeg[vertex_exnew]; e2 != EAT_LAST; e2 = graph->oeat[e2] )
253 {
254 const int head = graph->head[e2];
255
256 /* have we already extended to head? */
257 if( nodepos[head] < 0 )
258 {
259 const int stackposition = (-nodepos[head]) - 1;
260 const int edge_exold = edgestack[stackposition];
261 const SCIP_Real cost_exnew = graph->cost[e];
262 const SCIP_Real cost_exold = graph->cost[edge_exold];
263 const SCIP_Real cost_alt = graph->cost[e2];
264
265 assert(edge_exold >= 0);
266 assert(e != e2 && e != edge_exold && edge_exold != e2);
267
268 if( (SCIPisLT(scip, cost_alt, cost_exold) || SCIPisLT(scip, cost_alt, cost_exnew)) )
269 {
270 /* can we eliminate new edge? */
271 if( cost_exold <= cost_exnew && extensionmark[n - start] == EXEDGE_FREE )
272 {
273 assert(stackposition >= start);
274
275 extensionmark[stackposition - start] = EXEDGE_FIXED;
276 extensionmark[n - start] = EXEDGE_KILLED;
277 break;
278 }
279
280 /* can we eliminate old edge? */
281 if( cost_exold >= cost_exnew && extensionmark[stackposition - start] == EXEDGE_FREE )
282 {
283 extensionmark[stackposition - start] = EXEDGE_KILLED;
284 extensionmark[n - start] = EXEDGE_FIXED;
285
286 assert(head == graph->head[edge_exold]);
287 nodepos[head] = 0;
288
289 cleanup = TRUE;
290 }
291 }
292 }
293 }
294
295 if( extensionmark[n - start] != EXEDGE_KILLED )
296 {
297 edgestack[n++] = e;
298
299 assert(nodepos[vertex_exnew] == 0);
300 nodepos[vertex_exnew] = -n;
301 }
302 }
303
304 for( int i = start; i < n; i++ )
305 {
306 const int edge = edgestack[i];
307
308 assert(edge >= 0);
309 assert(nodepos[graph->head[edge]] <= 0);
310 assert(nodepos[graph->head[edge]] < 0 || (extensionmark[i - start] == EXEDGE_KILLED));
311
312 nodepos[graph->head[edge]] = 0;
313 }
314
315 if( cleanup )
316 {
317 int newn = start;
318 for( int i = start; i < n; i++ )
319 if( extensionmark[i - start] != EXEDGE_KILLED )
320 edgestack[newn++] = edgestack[i];
321
322
323 assert(newn < n);
324 n = newn;
325 }
326
327
328 return n;
329 }
330 #endif
331
332 /** remove latest edge from subtree and returns new tree cost */
333 static
shortenSubtree(SCIP * scip,const GRAPH * graph,const SCIP_Real * redcost,const int * treeedges,SCIP_Real treecostold,SCIP_Real treecostoffset,int dfsdepth,int lastnode,int * nodepos,int * ancestormark,unsigned int * nstallrounds)334 SCIP_Real shortenSubtree(
335 SCIP* scip, /**< SCIP */
336 const GRAPH* graph, /**< graph data structure */
337 const SCIP_Real* redcost, /**< reduced costs */
338 const int* treeedges, /**< edges of tree */
339 SCIP_Real treecostold, /**< old tree cost */
340 SCIP_Real treecostoffset, /**< tree cost offset */
341 int dfsdepth, /**< DFS depth */
342 int lastnode, /**< last node */
343 int* nodepos, /**< array to mark node position*/
344 int* ancestormark, /**< ancestor mark array */
345 unsigned int* nstallrounds /**< rounds since latest update */
346 )
347 {
348 const int lastedge = treeedges[dfsdepth - 1];
349
350 assert(dfsdepth >= 1 && lastnode >= 0);
351 assert(SCIPisGE(scip, treecostold, 0.0) && treecostoffset >= 0.0);
352
353 nodepos[lastnode] = 0;
354 unmarkAncestors(graph, lastedge, ancestormark);
355
356 /* recompute tree cost? */
357 if( (*nstallrounds)++ >= 9 )
358 {
359 SCIP_Real treecost = treecostoffset;
360
361 *nstallrounds = 0;
362
363 for( int i = 0; i < dfsdepth - 1; i++ )
364 treecost += redcost[treeedges[i]];
365
366 #if 0 // todo deleteeme
367 if( !SCIPisEQ(scip, treecost, treecostold - redcost[lastedge]) )
368 {
369 printf("%.15f %.15f \n", treecost, treecostold - redcost[lastedge]);
370 assert(0);
371 exit(1);
372 }
373 #endif
374
375 assert(SCIPisEQ(scip, treecost, treecostold - redcost[lastedge]));
376 }
377
378 return (treecostold - redcost[lastedge]);
379 }
380
381 /** extend subtree and return new nadded_edges */
382 static
extendSubtreeHead(SCIP * scip,const GRAPH * graph,const SCIP_Real * redcost,int curredge,int currhead,int dfsdepth,int nadded_edges,SCIP_Real * treecost,SCIP_Real * treecosts,int * nodepos,int * treeedges,int * edgestack,int * ancestormark)383 int extendSubtreeHead(
384 SCIP* scip, /**< SCIP */
385 const GRAPH* graph, /**< graph data structure */
386 const SCIP_Real* redcost, /**< reduced costs */
387 int curredge, /**< added edges */
388 int currhead, /**< latest node */
389 int dfsdepth, /**< DFS depth*/
390 int nadded_edges, /**< added edges */
391 SCIP_Real* treecost, /**< pointer to treecost */
392 SCIP_Real* treecosts, /**< edge costs of tree */
393 int* nodepos, /**< node position in tree */
394 int* treeedges, /**< edges of tree */
395 int* edgestack, /**< stack */
396 int* ancestormark /**< ancestor mark array */
397 )
398 {
399 int n = nadded_edges + 1;
400 nodepos[currhead] = dfsdepth + 1;
401 *treecost += redcost[curredge];
402 treeedges[dfsdepth] = curredge;
403 treecosts[dfsdepth] = graph->cost[curredge];
404
405 markAncestors(graph, curredge, ancestormark);
406
407 for( int e = graph->outbeg[currhead]; e != EAT_LAST; e = graph->oeat[e] )
408 if( !nodepos[graph->head[e]] )
409 edgestack[n++] = e;
410
411 return n;
412 }
413
414 /** extend subtree and return new nadded_edges */
415 static
extendSubtreeTail(const GRAPH * graph,const SCIP_Real * redcost,int curredge,int currtail,int dfsdepth,int nadded_edges,SCIP_Real * treecost,SCIP_Real * treecosts,int * nodepos,int * treeedges,int * edgestack,int * ancestormark)416 int extendSubtreeTail(
417 const GRAPH* graph, /**< graph data structure */
418 const SCIP_Real* redcost, /**< reduced costs */
419 int curredge, /**< added edges */
420 int currtail, /**< latest node */
421 int dfsdepth, /**< DFS depth*/
422 int nadded_edges, /**< added edges */
423 SCIP_Real* treecost, /**< pointer to treecost */
424 SCIP_Real* treecosts, /**< edge costs of tree */
425 int* nodepos, /**< node position in tree */
426 int* treeedges, /**< edges of tree */
427 int* edgestack, /**< stack */
428 int* ancestormark /**< ancestor mark array */
429 )
430 {
431 int n = nadded_edges + 1;
432 nodepos[currtail] = dfsdepth + 1;
433 *treecost += redcost[curredge];
434 treeedges[dfsdepth] = curredge;
435 treecosts[dfsdepth] = graph->cost[curredge];
436
437 markAncestors(graph, curredge, ancestormark);
438
439 for( int e = graph->inpbeg[currtail]; e != EAT_LAST; e = graph->ieat[e] )
440 if( !nodepos[graph->tail[e]] )
441 edgestack[n++] = e;
442
443 return n;
444 }
445
446
447 /** small helper */
448 static
updateEqArrays(int edge,unsigned int * eqstack,int * eqstack_size,SCIP_Bool * eqmark)449 void updateEqArrays(
450 int edge, /**< the edge */
451 unsigned int* eqstack, /**< stores edges that were used for equality comparison */
452 int* eqstack_size, /**< pointer to size of eqstack */
453 SCIP_Bool* eqmark /**< marks edges that were used for equality comparison */
454 )
455 {
456 const unsigned int halfedge = ((unsigned int) edge) / 2;
457
458 assert(edge >= 0);
459
460 if( !eqmark[halfedge] )
461 {
462 eqmark[halfedge] = TRUE;
463 eqstack[*eqstack_size] = halfedge;
464
465 *eqstack_size = *eqstack_size + 1;
466 }
467 }
468
469
470 /** compare with tree bottleneck */
471 static
bottleneckRuleOut(SCIP * scip,const GRAPH * graph,const SCIP_Real * treecosts,const SCIP_Real * basebottlenecks,const int * nodepos,const int * treeedges,SCIP_Real orgedgecost,SCIP_Real extedgecost,int orgedge,int extedge,int dfsdepth,SCIP_Bool allow_eq,unsigned int * eqstack,int * eqstack_size,SCIP_Bool * eqmark)472 SCIP_Bool bottleneckRuleOut(
473 SCIP* scip, /**< SCIP */
474 const GRAPH* graph, /**< graph data structure */
475 const SCIP_Real* treecosts, /**< costs of edges of the tree */
476 const SCIP_Real* basebottlenecks, /**< bottleneck costs for innode and basenode */
477 const int* nodepos, /**< node position in tree */
478 const int* treeedges, /**< edges in tree (corresponding to treecosts) */
479 SCIP_Real orgedgecost, /**< cost of original edge */
480 SCIP_Real extedgecost, /**< cost of short-cut edge */
481 int orgedge, /**< original edge */
482 int extedge, /**< short-cut edge */
483 int dfsdepth, /**< dfs depth */
484 SCIP_Bool allow_eq, /**< test for equality? */
485 unsigned int* eqstack, /**< stores edges that were used for equality comparison */
486 int* eqstack_size, /**< pointer to size of eqstack */
487 SCIP_Bool* eqmark /**< marks edges that were used for equality comparison */
488 )
489 {
490 int start;
491
492 assert(!allow_eq ||(eqstack != NULL && eqstack_size != NULL && eqmark != NULL));
493
494 if( allow_eq )
495 {
496 if( SCIPisLE(scip, extedgecost, orgedgecost) )
497 {
498 if( SCIPisEQ(scip, extedgecost, orgedgecost) )
499 updateEqArrays(orgedge, eqstack, eqstack_size, eqmark);
500
501 return TRUE;
502 }
503 }
504 else if( SCIPisLT(scip, extedgecost, orgedgecost) )
505 return TRUE;
506
507 if( nodepos[graph->tail[extedge]] > graph->knots )
508 {
509 start = nodepos[graph->tail[extedge]] - 1 - graph->knots;
510 assert(start >= 0 && start <= 3);
511
512 if( start <= 2 )
513 {
514 assert(basebottlenecks[start] > 0);
515
516 if( SCIPisLT(scip, extedgecost, basebottlenecks[start]) )
517 return TRUE;
518 }
519 start = 0;
520 }
521 else
522 {
523 start = nodepos[graph->tail[extedge]]; /* not -1! We save the incoming bottlenecks */
524 assert(start >= 1 && start <= dfsdepth);
525 assert(start < dfsdepth || graph->tail[orgedge] == graph->tail[extedge]);
526
527 }
528
529 for( int i = start; i < dfsdepth; i++ )
530 {
531 assert(treecosts[i] >= 0.0);
532
533 if( allow_eq )
534 {
535 if( SCIPisLE(scip, extedgecost, treecosts[i]) )
536 {
537 if( SCIPisEQ(scip, extedgecost, treecosts[i]) )
538 updateEqArrays(treeedges[i], eqstack, eqstack_size, eqmark);
539
540 return TRUE;
541 }
542 }
543 else if( SCIPisLT(scip, extedgecost, treecosts[i]) )
544 return TRUE;
545 }
546
547 return FALSE;
548 }
549
550 /** can subtree be ruled out? */
551 static
ruleOutSubtree(SCIP * scip,const GRAPH * graph,const SCIP_Real * treecosts,const SCIP_Real * basebottlenecks,const int * nodepos,const int * treeedges,SCIP_Real cutoff,SCIP_Real extendedcost,int dfsdepth,int curredge,SCIP_Bool allowequality,const int * ancestormark,unsigned int * eqstack,int * eqstack_size,SCIP_Bool * eqmark)552 SCIP_Bool ruleOutSubtree(
553 SCIP* scip, /**< SCIP */
554 const GRAPH* graph, /**< graph data structure */
555 const SCIP_Real* treecosts, /**< costs of edges of the tree */
556 const SCIP_Real* basebottlenecks, /**< bottleneck costs for innode and basenode */
557 const int* nodepos, /**< node position in tree */
558 const int* treeedges, /**< edges in tree (corresponding to treecosts) */
559 SCIP_Real cutoff, /**< cut-off bound */
560 SCIP_Real extendedcost, /**< cost of subtree */
561 int dfsdepth, /**< dfs depth */
562 int curredge, /**< latest edge */
563 SCIP_Bool allowequality, /**< rule out also in case of equality */
564 const int* ancestormark, /**< ancestor mark array */
565 unsigned int* eqstack, /**< stores edges that were used for equality comparison */
566 int* eqstack_size, /**< pointer to size of eqstack */
567 SCIP_Bool* eqmark /**< marks edges that were used for equality comparison */
568 )
569 {
570 if( allowequality ? (extendedcost >= cutoff) : SCIPisGT(scip, extendedcost, cutoff) )
571 {
572 return TRUE;
573 }
574 else
575 {
576 SCIP_Real currcost;
577 int currhead;
578
579 if( ancestormark != NULL )
580 {
581 int count = 0;
582 for( IDX* curr = graph->ancestors[curredge]; curr != NULL && count <= STP_DAEX_MAXGRAD; curr = curr->parent, count++ )
583 {
584 assert(curr->index >= 0 && curr->index / 2 < (MAX(graph->edges, graph->orgedges) / 2));
585 if( ancestormark[((unsigned) curr->index) / 2] )
586 return TRUE;
587 }
588 }
589
590 currcost = graph->cost[curredge];
591 currhead = graph->head[curredge];
592
593 for( int e = graph->inpbeg[currhead]; e != EAT_LAST; e = graph->ieat[e] )
594 {
595 int tail;
596 SCIP_Real ecost;
597
598 if( e == curredge )
599 continue;
600
601 assert(flipedge(e) != curredge);
602
603 tail = graph->tail[e];
604 ecost = graph->cost[e];
605
606 if( nodepos[tail] )
607 {
608 const SCIP_Bool allow_eq = (eqmark != NULL && eqmark[((unsigned) e) / 2] == FALSE);
609
610 if( bottleneckRuleOut(scip, graph, treecosts, basebottlenecks, nodepos, treeedges, currcost,
611 ecost, curredge, e, dfsdepth, allow_eq, eqstack, eqstack_size, eqmark) )
612 return TRUE;
613 }
614 else
615 {
616 const SCIP_Bool is_term = Is_term(graph->term[tail]);
617
618 for( int e2 = graph->outbeg[tail], count = 0; e2 != EAT_LAST && count < STP_DAEX_MAXGRAD;
619 e2 = graph->oeat[e2], count++ )
620 {
621 int head2;
622
623 if( e2 == e )
624 continue;
625
626 assert(flipedge(e2) != e && e2 != curredge && e2 != flipedge(curredge) );
627
628 head2 = graph->head[e2];
629 assert(head2 != tail);
630
631 if( nodepos[head2] )
632 {
633 const SCIP_Real extcost = is_term ? MAX(ecost, graph->cost[e2]) : (ecost + graph->cost[e2]);
634 const SCIP_Bool allow_eq = (eqmark != NULL && eqmark[((unsigned) e) / 2] == FALSE && eqmark[((unsigned) e2) / 2] == FALSE);
635
636 assert(head2 != currhead);
637
638 if( bottleneckRuleOut(scip, graph, treecosts, basebottlenecks,
639 nodepos, treeedges, currcost, extcost, curredge, flipedge(e2), dfsdepth, allow_eq, eqstack, eqstack_size, eqmark) )
640 {
641 return TRUE;
642 }
643 }
644 }
645 }
646 }
647 }
648
649 return FALSE;
650 }
651
652 /** updates node bounds for reduced cost fixings */
653 static
updateNodeFixingBounds(SCIP_Real * fixingbounds,const GRAPH * graph,const SCIP_Real * pathdist,const PATH * vnoi,SCIP_Real lpobjval,SCIP_Bool initialize)654 void updateNodeFixingBounds(
655 SCIP_Real* fixingbounds, /**< fixing bounds */
656 const GRAPH* graph, /**< graph data structure */
657 const SCIP_Real* pathdist, /**< shortest path distances */
658 const PATH* vnoi, /**< Voronoi paths */
659 SCIP_Real lpobjval, /**< LP objective */
660 SCIP_Bool initialize /**< initialize fixing bounds? */
661 )
662 {
663 const int nnodes = graph->knots;
664
665 assert(graph->stp_type == STP_SPG || graph->stp_type == STP_RSMT || !graph->extended);
666
667 if( initialize )
668 for( int k = 0; k < nnodes; k++ )
669 fixingbounds[k] = -FARAWAY;
670
671 for( int k = 0; k < nnodes; k++ )
672 {
673 SCIP_Real fixbnd;
674
675 if( !graph->mark[k] )
676 continue;
677
678 assert(!Is_pterm(graph->term[k]));
679
680 fixbnd = pathdist[k] + vnoi[k].dist + lpobjval;
681
682 if( fixbnd > fixingbounds[k] )
683 fixingbounds[k] = fixbnd;
684 }
685 }
686
687 /** updates node bounds for reduced cost replacement */
688 static
updateNodeReplaceBounds(SCIP * scip,SCIP_Real * replacebounds,const GRAPH * graph,const SCIP_Real * cost,const SCIP_Real * pathdist,const PATH * vnoi,const int * vbase,int * nodearrint,SCIP_Real lpobjval,SCIP_Real upperbound,int root,SCIP_Bool initialize,SCIP_Bool extendedsearch)689 SCIP_RETCODE updateNodeReplaceBounds(
690 SCIP* scip, /**< SCIP */
691 SCIP_Real* replacebounds, /**< replacement bounds */
692 const GRAPH* graph, /**< graph data structure */
693 const SCIP_Real* cost, /**< reduced costs */
694 const SCIP_Real* pathdist, /**< shortest path distances */
695 const PATH* vnoi, /**< Voronoi paths */
696 const int* vbase, /**< bases to Voronoi paths */
697 int* nodearrint, /**< for internal computations */
698 SCIP_Real lpobjval, /**< LP objective */
699 SCIP_Real upperbound, /**< upper bound */
700 int root, /**< DA root */
701 SCIP_Bool initialize, /**< initialize fixing bounds? */
702 SCIP_Bool extendedsearch /**< perform extended searching? */
703 )
704 {
705 unsigned int* eqstack = NULL;
706 SCIP_Bool* eqmark = NULL;
707 int outedges[STP_DABD_MAXDEGREE];
708 const int nnodes = graph->knots;
709 const SCIP_Real cutoff = upperbound - lpobjval;
710 const int halfnedges = graph->edges / 2;
711
712 assert(!SCIPisNegative(scip, cutoff));
713 assert(graph->stp_type == STP_SPG || graph->stp_type == STP_RSMT || !graph->extended);
714
715 if( initialize )
716 for( int k = 0; k < nnodes; k++ )
717 replacebounds[k] = -FARAWAY;
718
719 if( extendedsearch )
720 {
721 SCIP_CALL( SCIPallocCleanBufferArray(scip, &eqmark, halfnedges) );
722 SCIP_CALL( SCIPallocBufferArray(scip, &eqstack, halfnedges) );
723 }
724
725 for( int node = 0; node < nnodes; node++ )
726 {
727 const int degree = graph->grad[node];
728
729 if( degree >= 3 && !Is_gterm(graph->term[node]) )
730 {
731 SCIP_Real fixbnd;
732
733 /* bound already good enough? */
734 if( SCIPisLT(scip, upperbound, replacebounds[node]) )
735 continue;
736
737 fixbnd = pathdist[node] + vnoi[node].dist + vnoi[node + nnodes].dist + lpobjval;
738
739 assert(!Is_pterm(graph->term[node]));
740
741 /* Y-test for small degrees */
742 if( degree <= STP_DABD_MAXDEGREE && !SCIPisLT(scip, upperbound, fixbnd) )
743 {
744 int eqstack_size = 0;
745 int edgecount = 0;
746
747 fixbnd = FARAWAY;
748
749 for( int e = graph->outbeg[node]; e != EAT_LAST; e = graph->oeat[e] )
750 outedges[edgecount++] = e;
751
752 assert(edgecount == degree);
753
754 /* compute cost for each root and save minimum */
755 for( int i = 0; i < degree; i++ )
756 {
757 const int tmproot = graph->head[outedges[i]];
758 const int rootedge = flipedge(outedges[i]);
759
760 /* loop over all combinations of Y with root tmproot */
761 for( int j = i + 1; j <= i + degree - 2; j++ )
762 {
763 for( int k = j + 1; k <= i + degree - 1; k++ )
764 {
765 const int outedge1 = outedges[j % degree];
766 const int outedge2 = outedges[k % degree];
767 const int leaf1 = graph->head[outedge1];
768 const int leaf2 = graph->head[outedge2];
769
770 SCIP_Real tmpcostY;
771
772 assert(leaf1 != leaf2 && tmproot != leaf1 && tmproot != leaf2);
773 assert(vbase[leaf1] >= 0 || vnoi[leaf1].dist >= FARAWAY);
774 assert(vbase[leaf2] >= 0 || vnoi[leaf2].dist >= FARAWAY);
775
776 if( leaf1 == root || leaf2 == root )
777 continue;
778
779 tmpcostY = pathdist[tmproot] + cost[rootedge] + cost[outedge1] + cost[outedge2];
780
781 if( vbase[leaf1] != vbase[leaf2] )
782 {
783 tmpcostY += vnoi[leaf1].dist + vnoi[leaf2].dist;
784 }
785 else
786 {
787 /* also covers the case that leaf is a terminal */
788 const SCIP_Real leaf1far = vnoi[leaf1 + nnodes].dist;
789 const SCIP_Real leaf2far = vnoi[leaf2 + nnodes].dist;
790
791 assert(vbase[leaf1 + nnodes] >= 0 || leaf1far == FARAWAY);
792 assert(vbase[leaf2 + nnodes] >= 0 || leaf2far == FARAWAY);
793
794 tmpcostY += MIN(leaf1far + vnoi[leaf2].dist, vnoi[leaf1].dist + leaf2far);
795 }
796
797 if( tmpcostY < fixbnd )
798 {
799 if( extendedsearch && SCIPisLE(scip, tmpcostY, cutoff) )
800 {
801 int tree3outedges[2];
802 SCIP_Bool ruleout;
803 #ifndef NDEBUG
804 const SCIP_Real tmpcostYorg = tmpcostY;
805 #endif
806 tree3outedges[0] = outedge1;
807 tree3outedges[1] = outedge2;
808
809 SCIP_CALL( reduce_check3Tree(scip, graph, root, cost, pathdist, vnoi, vbase, cutoff, tree3outedges, rootedge, nodearrint,
810 &tmpcostY, &ruleout, eqstack, &eqstack_size, eqmark) );
811
812 if( ruleout )
813 tmpcostY = FARAWAY;
814
815 #ifndef NDEBUG
816 assert(tmpcostY >= tmpcostYorg);
817 #endif
818 }
819
820 if( tmpcostY < fixbnd )
821 fixbnd = tmpcostY;
822 }
823 }
824 } /* Y loop */
825 } /* root loop */
826
827 fixbnd += lpobjval;
828
829 assert(SCIPisGE(scip, fixbnd, pathdist[node] + vnoi[node].dist + vnoi[node + nnodes].dist + lpobjval)
830 || fixbnd >= FARAWAY);
831
832 if( extendedsearch )
833 {
834 for( int i = 0; i < eqstack_size; i++ )
835 eqmark[eqstack[i]] = FALSE;
836
837 for( int i = 0; i < halfnedges; i++ )
838 assert(eqmark[i] == FALSE);
839 }
840 }
841
842 if( fixbnd > replacebounds[node] )
843 replacebounds[node] = fixbnd;
844 }
845 }
846 if( extendedsearch )
847 {
848 assert(eqstack != NULL && eqmark != NULL);
849
850 SCIPfreeBufferArray(scip, &eqstack);
851 SCIPfreeCleanBufferArray(scip, &eqmark);
852 }
853
854 return SCIP_OKAY;
855 }
856
857 /** updates edge fixing bounds for reduced cost fixings */
858 static
updateEdgeFixingBounds(SCIP_Real * fixingbounds,const GRAPH * graph,const SCIP_Real * cost,const SCIP_Real * pathdist,const PATH * vnoi,SCIP_Real lpobjval,int extnedges,SCIP_Bool initialize,SCIP_Bool undirected)859 void updateEdgeFixingBounds(
860 SCIP_Real* fixingbounds, /**< fixing bounds */
861 const GRAPH* graph, /**< graph data structure */
862 const SCIP_Real* cost, /**< reduced costs */
863 const SCIP_Real* pathdist, /**< shortest path distances */
864 const PATH* vnoi, /**< Voronoi paths */
865 SCIP_Real lpobjval, /**< LP objective */
866 int extnedges, /**< number of edges for extended problem */
867 SCIP_Bool initialize, /**< initialize fixing bounds? */
868 SCIP_Bool undirected /**< only consider undirected edges */
869 )
870 {
871 const int nnodes = graph->knots;
872
873 assert(graph->stp_type == STP_SPG || graph->stp_type == STP_RSMT || !graph->extended);
874 assert(extnedges > 0);
875
876 if( initialize )
877 for( int e = 0; e < extnedges; e++ )
878 fixingbounds[e] = -FARAWAY;
879
880 for( int k = 0; k < nnodes; k++ )
881 {
882 if( !graph->mark[k] )
883 continue;
884
885 assert(!Is_pterm(graph->term[k]));
886
887 if( undirected )
888 {
889 for( int e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
890 {
891 const int head = graph->head[e];
892
893 if( graph->mark[head] )
894 {
895 const int erev = flipedge(e);
896 const SCIP_Real fixbnd = pathdist[k] + cost[e] + vnoi[head].dist + lpobjval;
897 const SCIP_Real fixbndrev = pathdist[head] + cost[erev] + vnoi[k].dist + lpobjval;
898 const SCIP_Real minbnd = MIN(fixbnd, fixbndrev);
899
900 assert(fixingbounds[e] == fixingbounds[erev]);
901
902 if( minbnd > fixingbounds[e] )
903 {
904 fixingbounds[e] = minbnd;
905 fixingbounds[erev] = minbnd;
906 }
907 }
908 }
909 }
910 else
911 {
912 for( int e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
913 {
914 const int head = graph->head[e];
915
916 if( graph->mark[head] )
917 {
918 const SCIP_Real fixbnd = pathdist[k] + cost[e] + vnoi[head].dist + lpobjval;
919
920 if( fixbnd > fixingbounds[e] )
921 fixingbounds[e] = fixbnd;
922 }
923 }
924 }
925 }
926 }
927
928 /** eliminate nodes by using fixing-bounds and reduced costs */
929 static
reduceWithNodeFixingBounds(SCIP * scip,GRAPH * graph,GRAPH * transgraph,const SCIP_Real * fixingbounds,SCIP_Real upperbound)930 int reduceWithNodeFixingBounds(
931 SCIP* scip, /**< SCIP data structure */
932 GRAPH* graph, /**< graph data structure */
933 GRAPH* transgraph, /**< graph data structure or NULL */
934 const SCIP_Real* fixingbounds, /**< fixing bounds */
935 SCIP_Real upperbound /**< best upperbound */
936 )
937 {
938 int nfixed = 0;
939 int nnodes = graph->knots;
940
941 assert(graph->stp_type == STP_SPG || graph->stp_type == STP_RSMT || !graph->extended);
942
943 for( int k = 0; k < nnodes; k++ )
944 {
945 if( !graph->mark[k] || Is_term(graph->term[k]) )
946 continue;
947
948 assert(!Is_pterm(graph->term[k]));
949
950 if( SCIPisLT(scip, upperbound, fixingbounds[k]) )
951 {
952 SCIPdebugMessage("delete knot %d %f < %f %d\n", k, upperbound, fixingbounds[k], graph->grad[k]);
953
954 graph_knot_del(scip, graph, k, TRUE);
955 if( transgraph != NULL )
956 graph_knot_del(scip, transgraph, k, FALSE);
957 }
958 }
959 return nfixed;
960 }
961
962 static
reduceWithNodeReplaceBounds(SCIP * scip,GRAPH * graph,const PATH * vnoi,const SCIP_Real * pathdist,const SCIP_Real * cost,const SCIP_Real * replacebounds,int * nodetouched,SCIP_Real lpobjval,SCIP_Real upperbound)963 int reduceWithNodeReplaceBounds(
964 SCIP* scip, /**< SCIP data structure */
965 GRAPH* graph, /**< graph data structure */
966 const PATH* vnoi, /**< Voronoi data structure */
967 const SCIP_Real* pathdist, /**< distance array from shortest path calculations */
968 const SCIP_Real* cost, /**< dual ascent costs */
969 const SCIP_Real* replacebounds, /**< replacement bounds */
970 int* nodetouched, /**< node array for internal computations */
971 SCIP_Real lpobjval, /**< lower bound corresponding to pathdist and vnoi */
972 SCIP_Real upperbound /**< best upperbound */
973 )
974 {
975 int adjvert[STP_DABD_MAXDEGREE];
976 SCIP_Real cutoffs[STP_DABD_MAXDNEDGES];
977 SCIP_Real cutoffsrev[STP_DABD_MAXDNEDGES];
978 int nfixed = 0;
979 const int nnodes = graph->knots;
980
981 for( int k = 0; k < nnodes; k++ )
982 nodetouched[k] = 0;
983
984 /* main loop */
985 for( int degree = 3; degree <= STP_DABD_MAXDEGREE; degree++ )
986 {
987 for( int k = 0; k < nnodes; k++ )
988 {
989 if( degree != graph->grad[k] || Is_gterm(graph->term[k]) )
990 continue;
991
992 if( SCIPisLT(scip, upperbound, replacebounds[k]) && nodetouched[k] == 0 )
993 {
994 int edgecount = 0;
995 SCIP_Bool success;
996
997 for( int e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
998 nodetouched[graph->head[e]]++;
999
1000 /* fill cutoff */
1001
1002 for( int e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
1003 adjvert[edgecount++] = graph->head[e];
1004
1005 assert(edgecount == degree);
1006
1007 edgecount = 0;
1008 for( int i = 0; i < degree - 1; i++ )
1009 {
1010 const int vert = adjvert[i];
1011 for( int i2 = i + 1; i2 < degree; i2++ )
1012 {
1013 const int vert2 = adjvert[i2];
1014
1015 assert(edgecount < STP_DABD_MAXDNEDGES);
1016
1017 cutoffs[edgecount] = upperbound - lpobjval - (pathdist[vert] + vnoi[vert2].dist);
1018 cutoffsrev[edgecount] = upperbound - lpobjval - (pathdist[vert2] + vnoi[vert].dist);
1019
1020 edgecount++;
1021 }
1022 }
1023
1024 assert(edgecount > 0);
1025
1026 /* try to eliminate */
1027
1028 SCIP_CALL_ABORT(graph_knot_delPseudo(scip, graph, cost, cutoffs, cutoffsrev, k, &success));
1029
1030 if( success )
1031 nfixed++;
1032 else
1033 {
1034 assert(graph->grad[k] == degree);
1035
1036 for( int e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
1037 {
1038 nodetouched[graph->head[e]]--;
1039 assert(nodetouched[graph->head[e]] >= 0);
1040 }
1041 }
1042 }
1043 }
1044 }
1045
1046 return nfixed;
1047 }
1048 /** eliminate edges by using fixing-bounds and reduced costs */
1049 static
reduceWithEdgeFixingBounds(SCIP * scip,GRAPH * graph,GRAPH * transgraph,const SCIP_Real * fixingbounds,const int * result,SCIP_Real upperbound)1050 int reduceWithEdgeFixingBounds(
1051 SCIP* scip, /**< SCIP data structure */
1052 GRAPH* graph, /**< graph data structure */
1053 GRAPH* transgraph, /**< graph data structure or NULL */
1054 const SCIP_Real* fixingbounds, /**< fixing bounds */
1055 const int* result, /**< solution */
1056 SCIP_Real upperbound /**< best upperbound */
1057 )
1058 {
1059 int nfixed = 0;
1060 int nnodes = graph->knots;
1061 const SCIP_Bool solgiven = (result != NULL);
1062
1063 assert(graph->stp_type == STP_SPG || graph->stp_type == STP_RSMT || !graph->extended);
1064 assert(!solgiven || graph_sol_valid(scip, graph, result));
1065
1066 for( int k = 0; k < nnodes; k++ )
1067 {
1068 int e;
1069 if( !graph->mark[k] )
1070 continue;
1071
1072 assert(!Is_pterm(graph->term[k]));
1073
1074 e = graph->outbeg[k];
1075 while( e != EAT_LAST )
1076 {
1077 const int enext = graph->oeat[e];
1078
1079 if( graph->mark[graph->head[e]] )
1080 {
1081 const int erev = flipedge(e);
1082 SCIP_Bool delete;
1083
1084 if( !solgiven || result[e] == CONNECT || result[erev] == CONNECT )
1085 delete = (SCIPisLT(scip, upperbound, fixingbounds[e]) && SCIPisLT(scip, upperbound, fixingbounds[erev]));
1086 else
1087 delete = (upperbound <= fixingbounds[e] && upperbound <= fixingbounds[erev]);
1088
1089 if( delete )
1090 {
1091 SCIPdebugMessage("delete edge %d \n", e);
1092
1093 graph_edge_del(scip, graph, e, TRUE);
1094 if( transgraph != NULL )
1095 graph_edge_del(scip, transgraph, e, FALSE);
1096 nfixed++;
1097 }
1098 }
1099
1100 e = enext;
1101 }
1102 }
1103
1104 return nfixed;
1105 }
1106
1107
1108 /** find roots for PC and MW during DA reduction */
1109 static
findDaRoots(SCIP * scip,GRAPH * graph,const GRAPH * transgraph,const SCIP_Real * cost,const SCIP_Real * bestcost,SCIP_Real lpobjval,SCIP_Real bestlpobjval,SCIP_Real upperbound,SCIP_Real prizesum,SCIP_Bool rerun,SCIP_Bool probrooted,PATH * vnoi,int * state,int * pathedge,int * vbase,STP_Bool * nodearrchar,int * roots,int * rootscount)1110 void findDaRoots(
1111 SCIP* scip, /**< SCIP data structure */
1112 GRAPH* graph, /**< graph data structure */
1113 const GRAPH* transgraph, /**< graph data structure */
1114 const SCIP_Real* cost, /**< da reduced cost */
1115 const SCIP_Real* bestcost, /**< best incumbent da reduced cost */
1116 SCIP_Real lpobjval, /**< da lower bound */
1117 SCIP_Real bestlpobjval, /**< best da lower bound */
1118 SCIP_Real upperbound, /**< da upper bound */
1119 SCIP_Real prizesum, /**< sum of positive prizes */
1120 SCIP_Bool rerun, /**< not the first run? */
1121 SCIP_Bool probrooted, /**< is transgraph a rooted RMW or RPC? */
1122 PATH* vnoi, /**< SP array */
1123 int* state, /**< array */
1124 int* pathedge, /**< array */
1125 int* vbase, /**< array */
1126 STP_Bool* nodearrchar, /**< bool array */
1127 int* roots, /**< root array */
1128 int* rootscount /**< number of roots */
1129 )
1130 {
1131 const int root = graph->source;
1132 const int nnodes = graph->knots;
1133 const SCIP_Bool graphextended = graph->extended;
1134 int nroots = *rootscount;
1135
1136 assert(transgraph->extended);
1137
1138 if( graphextended )
1139 graph_pc_2org(graph);
1140
1141 assert(rerun || nroots == 0);
1142
1143 /*
1144 * get possible roots
1145 */
1146
1147 BMSclearMemoryArray(nodearrchar, nnodes);
1148
1149 if( rerun )
1150 for( int i = 0; i < nroots; i++ )
1151 nodearrchar[roots[i]] = TRUE;
1152
1153 if( probrooted )
1154 {
1155 const int transroot = transgraph->source;
1156
1157 assert(transgraph->term2edge != NULL);
1158
1159 for( int e = transgraph->outbeg[transroot]; e != EAT_LAST; e = transgraph->oeat[e] )
1160 {
1161 const int pseudoterm = transgraph->head[e];
1162
1163 if( Is_term(transgraph->term[pseudoterm]) && transgraph->term2edge[pseudoterm] >= 0 )
1164 {
1165 int realterm;
1166 assert(transgraph->tail[transgraph->term2edge[pseudoterm]] == pseudoterm);
1167
1168 realterm = transgraph->head[transgraph->term2edge[pseudoterm]];
1169 assert(Is_pterm(transgraph->term[realterm]));
1170
1171 if( rerun && nodearrchar[realterm] )
1172 continue;
1173
1174 if( SCIPisGT(scip, cost[e], upperbound - lpobjval) || SCIPisGT(scip, bestcost[e], upperbound - bestlpobjval) )
1175 {
1176 assert(SCIPisPositive(scip, graph->prize[realterm]));
1177 assert(nroots < graph->terms);
1178
1179 roots[nroots++] = realterm;
1180 nodearrchar[realterm] = TRUE;
1181 }
1182 }
1183 }
1184 }
1185 else
1186 {
1187 for( int e = graph->outbeg[root]; e != EAT_LAST; e = graph->oeat[e] )
1188 {
1189 const int pseudoterm = graph->head[e];
1190
1191 if( Is_pterm(graph->term[pseudoterm]) )
1192 {
1193 int realterm = -1;
1194 int e3 = -1;
1195
1196 assert(graph->grad[pseudoterm] == 2);
1197
1198 for( int e2 = transgraph->inpbeg[pseudoterm]; e2 != EAT_LAST; e2 = transgraph->ieat[e2] )
1199 {
1200 if( transgraph->cost[e2] == 0.0 )
1201 realterm = transgraph->tail[e2];
1202 else
1203 e3 = e2;
1204 }
1205
1206 if( rerun && nodearrchar[realterm] )
1207 continue;
1208
1209 assert(realterm >= 0);
1210 assert(e3 >= 0);
1211 assert(realterm != root);
1212 assert(SCIPisEQ(scip, graph->cost[e], graph->cost[e3]));
1213 assert(graph->mark[realterm]);
1214
1215 if( SCIPisGT(scip, cost[e3], upperbound - lpobjval) || SCIPisGT(scip, bestcost[e3], upperbound - bestlpobjval) )
1216 {
1217 assert(SCIPisPositive(scip, graph->prize[realterm]));
1218 assert(nroots < graph->terms);
1219
1220 roots[nroots++] = realterm;
1221 nodearrchar[realterm] = TRUE;
1222 }
1223 }
1224 }
1225 }
1226
1227 /* could more roots be found? */
1228 if( nroots > *rootscount && graph->terms > 2 )
1229 {
1230 /*
1231 * try to find additional roots by shortest path computations
1232 */
1233
1234 SCIP_Real maxprize = 0.0;
1235 int qstart = *rootscount;
1236
1237 for( int i = 0; i < nnodes; i++ )
1238 {
1239 state[i] = UNKNOWN;
1240 vnoi[i].dist = FARAWAY;
1241 vnoi[i].edge = UNKNOWN;
1242
1243 if( Is_term(graph->term[i]) && !nodearrchar[i] && graph->mark[i] && maxprize < graph->prize[i] && graph->prize[i] < prizesum )
1244 maxprize = graph->prize[i];
1245 }
1246
1247 for( int rounds = 0; rounds < 10 && qstart != nroots; rounds++ )
1248 {
1249 const int oldnroots = nroots;
1250
1251 SCIPdebugMessage("root-finding round %d \n", rounds);
1252
1253 for( int i = qstart; i < nroots; i++ )
1254 {
1255 int nlabeled;
1256 const int term = roots[i];
1257 assert(nodearrchar[term]);
1258 assert(Is_term(graph->term[term]));
1259 assert(SCIPisPositive(scip, graph->prize[term]));
1260
1261 /* look for terminals in neighborhood of term */
1262 graph_sdPaths(scip, graph, vnoi, graph->cost, maxprize, pathedge, state, vbase, &nlabeled, term, term, 200);
1263
1264 /* check labelled nodes */
1265 for( int k = 0; k < nlabeled; k++ )
1266 {
1267 const int l = vbase[k];
1268
1269 assert(graph->mark[l]);
1270 assert(state[l] != UNKNOWN);
1271
1272 if( Is_term(graph->term[l]) && !nodearrchar[l] && vnoi[l].dist <= graph->prize[l] )
1273 {
1274 roots[nroots++] = l;
1275 nodearrchar[l] = TRUE;
1276
1277 SCIPdebugMessage("added new root %d dist: %f, prize: %f \n", l, vnoi[l].dist, graph->prize[l]);
1278 }
1279
1280 /* restore data structures */
1281 state[l] = UNKNOWN;
1282 vnoi[l].dist = FARAWAY;
1283 vnoi[l].edge = UNKNOWN;
1284 }
1285 }
1286 #ifndef NDEBUG
1287 for( int k = 0; k < nnodes; k++ )
1288 {
1289 assert(state[k] == UNKNOWN);
1290 assert(vnoi[k].dist == FARAWAY);
1291 assert(vnoi[k].edge == UNKNOWN);
1292 }
1293 #endif
1294
1295 qstart = oldnroots;
1296 }
1297 SCIPdebugMessage("number of terminals found by DA: %d \n", nroots);
1298 }
1299
1300 if( rerun )
1301 SCIPdebugMessage("new roots in rerun %d \n", nroots - *rootscount);
1302
1303 *rootscount = nroots;
1304
1305 if( graphextended )
1306 graph_pc_2trans(graph);
1307 }
1308
1309
1310 /** set prize of marked terminals to blocked */
1311 static
markPcMwRoots(SCIP * scip,const int * roots,int nrootsold,int nrootsnew,SCIP_Real prizesum,GRAPH * graph,SCIP_Bool * userec,STPSOLPOOL ** solpool)1312 void markPcMwRoots(
1313 SCIP* scip, /**< SCIP data structure */
1314 const int* roots, /**< root array */
1315 int nrootsold, /**< old number of roots */
1316 int nrootsnew, /**< new number of roots */
1317 SCIP_Real prizesum, /**< sum of positive prizes */
1318 GRAPH* graph, /**< graph data structure */
1319 SCIP_Bool* userec, /**< recombination? */
1320 STPSOLPOOL** solpool /**< solution pool */
1321 )
1322 {
1323 const int root = graph->source;
1324 const SCIP_Bool graphextended = graph->extended;
1325
1326 if( graphextended )
1327 graph_pc_2org(graph);
1328
1329 if( *userec && *solpool != NULL )
1330 {
1331 *userec = FALSE;
1332 SCIPStpHeurRecFreePool(scip, solpool);
1333 *solpool = NULL;
1334 }
1335
1336 assert(graph != NULL);
1337 assert(roots != NULL);
1338 assert(!graph->extended);
1339 assert(nrootsnew >= 0 && nrootsold >= 0);
1340
1341 for( int i = nrootsold; i < nrootsnew; i++ )
1342 {
1343 int e;
1344 const int term = roots[i];
1345 const int pterm = graph->head[graph->term2edge[term]];
1346
1347 assert(Is_term(graph->term[term]));
1348 assert(SCIPisPositive(scip, graph->prize[term]));
1349 assert(pterm >= 0);
1350 assert(Is_pterm(graph->term[pterm]));
1351
1352 for( e = graph->inpbeg[pterm]; e != EAT_LAST; e = graph->ieat[e] )
1353 if( root == graph->tail[e] )
1354 break;
1355
1356 assert(e != EAT_LAST);
1357 assert(SCIPisEQ(scip, graph->prize[term], graph->cost[e]));
1358
1359 graph->prize[term] = prizesum;
1360 graph->cost[e] = prizesum;
1361 }
1362
1363 if( graphextended )
1364 graph_pc_2trans(graph);
1365 }
1366
1367 /** are the dual costs still valid, i.e. are there zero cost paths from the root to all terminals? */
1368 static
dualCostIsValid(SCIP * scip,GRAPH * transgraph,const SCIP_Real * cost,int * nodearrint,STP_Bool * nodearrbool)1369 SCIP_Bool dualCostIsValid(
1370 SCIP* scip, /**< SCIP data structure */
1371 GRAPH* transgraph, /**< graph data structure */
1372 const SCIP_Real* cost, /**< dual ascent costs */
1373 int* nodearrint, /**< int array */
1374 STP_Bool* nodearrbool /**< bool array */
1375 )
1376 {
1377 int* const queue = nodearrint;
1378 STP_Bool* visited = nodearrbool;
1379 int qsize;
1380 const int root = transgraph->source;
1381 const int nnodes = transgraph->knots;
1382
1383 /*
1384 * construct new graph corresponding to zero cost paths from the root to all terminals
1385 */
1386
1387 BMSclearMemoryArray(visited, nnodes);
1388
1389 qsize = 0;
1390 visited[root] = TRUE;
1391 queue[qsize++] = root;
1392
1393 /* DFS */
1394 while( qsize )
1395 {
1396 const int k = queue[--qsize];
1397
1398 /* traverse outgoing arcs */
1399 for( int a = transgraph->outbeg[k]; a != EAT_LAST; a = transgraph->oeat[a] )
1400 {
1401 const int head = transgraph->head[a];
1402
1403 if( !visited[head] && SCIPisZero(scip, cost[a]) )
1404 {
1405 visited[head] = TRUE;
1406 queue[qsize++] = head;
1407 }
1408 }
1409 assert(qsize <= nnodes);
1410 }
1411
1412 for( int k = 0; k < nnodes; k++ )
1413 if( Is_term(transgraph->term[k]) && !visited[k] )
1414 {
1415 return FALSE;
1416 }
1417
1418 return TRUE;
1419 }
1420
1421
1422 /** pertubate edge costs for PCMW dual-ascent */
1423 static
pertubateEdgeCosts(SCIP * scip,const GRAPH * graph,GRAPH * transgraph,const int * result,STP_Bool * nodearrchar,int randomize)1424 void pertubateEdgeCosts(
1425 SCIP* scip,
1426 const GRAPH* graph,
1427 GRAPH* transgraph,
1428 const int* result,
1429 STP_Bool* nodearrchar,
1430 int randomize
1431 )
1432 {
1433 int e;
1434 const int root = graph->source;
1435 const int newroot = transgraph->source;
1436 const int nnodes = graph->knots;
1437 const int nedges = graph->edges;
1438
1439 BMSclearMemoryArray(nodearrchar, nnodes);
1440
1441 /* mark all vertices visited in regular graph */
1442 for( e = 0; e < nedges; e++ )
1443 if( result[e] == CONNECT && graph->tail[e] != root )
1444 nodearrchar[graph->head[e]] = TRUE;
1445 srand((unsigned)graph->terms);
1446
1447 if( graph->stp_type != STP_MWCSP )
1448 {
1449 SCIP_Real pratio = PERTUBATION_RATIO_PC;
1450
1451 for( int k = 0; k < nnodes; k++ )
1452 {
1453 assert(Is_gterm(graph->term[k]) == Is_gterm(transgraph->term[k]) || transgraph->grad[k] == 0);
1454
1455 if( randomize > 8 )
1456 pratio = ((SCIP_Real)(rand() % 10)) / (100.0) - 1.0 / 100.0;
1457 else if( randomize > 6 )
1458 pratio = ((SCIP_Real)(rand() % 10)) / (200.0);
1459 else if( randomize > 4 )
1460 pratio = ((SCIP_Real)(rand() % 10)) / (300.0);
1461 else if( randomize > 0 )
1462 pratio = ((SCIP_Real)(rand() % 10)) / 1000.0;
1463 else
1464 pratio = PERTUBATION_RATIO_PC + ((SCIP_Real)(rand() % 10)) / 1000.0;
1465
1466 assert(SCIPisPositive(scip, 1.0 - pratio));
1467 assert(SCIPisPositive(scip, 1.0 + pratio));
1468
1469 if( !Is_gterm(graph->term[k]) )
1470 {
1471 for( e = transgraph->inpbeg[k]; e != EAT_LAST; e = transgraph->ieat[e] )
1472 {
1473 assert(transgraph->tail[e] != root);
1474
1475 if( result[e] == CONNECT || result[flipedge(e)] == CONNECT )
1476 transgraph->cost[e] *= 1.0 - pratio;
1477 else
1478 transgraph->cost[e] *= 1.0 + pratio;
1479 }
1480 }
1481 else if( Is_term(transgraph->term[k]) && k != root && k != newroot )
1482 {
1483 assert(transgraph->grad[k] == 2);
1484
1485 for( e = transgraph->inpbeg[k]; e != EAT_LAST; e = transgraph->ieat[e] )
1486 if( SCIPisPositive(scip, transgraph->cost[e]) )
1487 {
1488 assert(!Is_pterm(transgraph->term[transgraph->tail[e]]));
1489 assert(transgraph->tail[e] != root);
1490 assert(result[flipedge(e)] != CONNECT);
1491
1492 if( result[e] == CONNECT )
1493 transgraph->cost[e] *= 1.0 - pratio;
1494 else
1495 transgraph->cost[e] *= 1.0 + pratio;
1496
1497 assert(SCIPisPositive(scip, transgraph->cost[e]));
1498 }
1499 }
1500 }
1501
1502 return;
1503 }
1504
1505 for( int k = 0; k < nnodes; k++ )
1506 {
1507 SCIP_Real pratio = PERTUBATION_RATIO;
1508
1509 assert(Is_gterm(graph->term[k]) == Is_gterm(transgraph->term[k]));
1510
1511 if( randomize > 8 )
1512 pratio = ((SCIP_Real)(rand() % 10)) / (50.0) - 1.0 / 10.0;
1513 else if( randomize > 6 )
1514 pratio = ((SCIP_Real)(rand() % 10)) / (20.0);
1515 else if( randomize > 4 )
1516 pratio = ((SCIP_Real)(rand() % 10)) / (30.0);
1517 else if( randomize > 0 )
1518 pratio = ((SCIP_Real)(rand() % 10)) / 100.0;
1519 else
1520 pratio = PERTUBATION_RATIO + ((SCIP_Real)(rand() % 10)) / 200.0;
1521
1522 if( !Is_gterm(graph->term[k]) )
1523 {
1524 if( nodearrchar[k] )
1525 {
1526 for( e = transgraph->inpbeg[k]; e != EAT_LAST; e = transgraph->ieat[e] )
1527 transgraph->cost[e] *= 1.0 - pratio;
1528 }
1529 else
1530 {
1531 for( e = transgraph->inpbeg[k]; e != EAT_LAST; e = transgraph->ieat[e] )
1532 transgraph->cost[e] *= 1.0 + pratio;
1533 }
1534 }
1535 else if( Is_term(transgraph->term[k]) && k != root && k != newroot )
1536 {
1537 assert(transgraph->grad[k] == 2);
1538
1539 if( nodearrchar[k] )
1540 {
1541 for( e = transgraph->inpbeg[k]; e != EAT_LAST; e = transgraph->ieat[e] )
1542 if( SCIPisPositive(scip, transgraph->cost[e]) )
1543 {
1544 assert(!Is_pterm(transgraph->term[transgraph->tail[e]]));
1545
1546 transgraph->cost[e] *= 1.0 + pratio;
1547 }
1548 }
1549 else
1550 {
1551 for( e = transgraph->inpbeg[k]; e != EAT_LAST; e = transgraph->ieat[e] )
1552 if( SCIPisPositive(scip, transgraph->cost[e]) )
1553 {
1554 assert(!Is_pterm(transgraph->term[transgraph->tail[e]]));
1555 transgraph->cost[e] *= 1.0 - pratio;
1556 }
1557 }
1558 }
1559 }
1560 }
1561
1562 /** order roots */
1563 static
orderDaRoots(SCIP * scip,const GRAPH * graph,int * terms,int nterms,SCIP_Bool randomize,SCIP_RANDNUMGEN * randnumgen)1564 SCIP_RETCODE orderDaRoots(
1565 SCIP* scip, /**< SCIP data structure */
1566 const GRAPH* graph, /**< graph structure */
1567 int* terms, /**< sol int array corresponding to upper bound */
1568 int nterms, /**< number of terminals */
1569 SCIP_Bool randomize, /**< randomize */
1570 SCIP_RANDNUMGEN* randnumgen /**< random number generator */
1571 )
1572 {
1573 int* termdegs;
1574 int maxdeg = 0;
1575
1576 assert(terms != NULL);
1577 assert(nterms > 0);
1578
1579 SCIP_CALL( SCIPallocBufferArray(scip, &termdegs, nterms) );
1580
1581 for( int i = 0; i < nterms; i++ )
1582 {
1583 const int grad = graph->grad[terms[i]];
1584 assert(terms[i] >= 0);
1585
1586 termdegs[i] = -grad;
1587
1588 if( grad > maxdeg )
1589 maxdeg = termdegs[i];
1590 }
1591
1592 if( randomize )
1593 for( int i = 0; i < nterms; i++ )
1594 termdegs[i] -= SCIPrandomGetInt(randnumgen, 0, maxdeg);
1595
1596 SCIPsortIntInt(termdegs, terms, nterms);
1597
1598 SCIPfreeBufferArray(scip, &termdegs);
1599
1600 return SCIP_OKAY;
1601 }
1602
1603
1604 /** compute primal solution during dual-ascent routine for PCSTP or MWCSP */
1605 static
computeDaSolPcMw(SCIP * scip,GRAPH * graph,STPSOLPOOL * pool,PATH * vnoi,const SCIP_Real * cost,SCIP_Real * pathdist,SCIP_Real * upperbound,int * result1,int * result2,int * vbase,int * pathedge,STP_Bool * nodearrchar,SCIP_Bool * apsol)1606 SCIP_RETCODE computeDaSolPcMw(
1607 SCIP* scip, /**< SCIP data structure */
1608 GRAPH* graph, /**< graph data structure */
1609 STPSOLPOOL* pool, /**< solution pool */
1610 PATH* vnoi, /**< Voronoi data structure */
1611 const SCIP_Real* cost, /**< dual ascent costs */
1612 SCIP_Real* pathdist, /**< distance array from shortest path calculations */
1613 SCIP_Real* upperbound, /**< upperbound pointer */
1614 int* result1, /**< sol int array corresponding to upper bound */
1615 int* result2, /**< sol int array */
1616 int* vbase, /**< int array */
1617 int* pathedge, /**< int array */
1618 STP_Bool* nodearrchar, /**< node array storing solution vertices */
1619 SCIP_Bool* apsol /**< ascend-prune sol? */
1620 )
1621 {
1622 SCIP_Real ub;
1623 const int nedges = graph->edges;
1624 SCIP_Bool success;
1625
1626 /* compute new solution and store it in result2 */
1627
1628 SCIP_CALL( SCIPStpHeurAscendPruneRun(scip, NULL, graph, cost, result2, vbase, -1, nodearrchar, &success, FALSE) );
1629 assert(success);
1630
1631 SCIP_CALL( SCIPStpHeurLocalRun(scip, graph, graph->cost, result2) );
1632
1633 assert(graph_sol_valid(scip, graph, result2));
1634
1635 /* compute objective value */
1636 ub = graph_sol_getObj(graph->cost, result2, 0.0, nedges);
1637 SCIPdebugMessage("DA: first new sol value in computeDaSolPcMw: %f ... old value: %f \n", ub, *upperbound);
1638
1639 /* try recombination? */
1640 if( pool != NULL )
1641 {
1642 SCIPdebugMessage("ub %f vs best sol %f\n", ub, pool->sols[0]->obj);
1643 SCIP_CALL( SCIPStpHeurRecAddToPool(scip, ub, result2, pool, &success) );
1644
1645 #ifdef SCIP_DEBUG
1646 for( int i = 0; i < pool->size; i++ )
1647 {
1648 printf(" %f ", pool->sols[i]->obj);
1649 }
1650 printf("\n ");
1651 #endif
1652
1653 if( success && pool->size >= 2 )
1654 {
1655 /* get index of just added solution */
1656 int solindex = pool->maxindex;
1657
1658 SCIP_Bool solfound;
1659
1660 SCIPdebugMessage("POOLSIZE %d \n", pool->size);
1661
1662 SCIP_CALL( SCIPStpHeurRecRun(scip, pool, NULL, NULL, graph, NULL, &solindex, 3, pool->size, FALSE, &solfound) );
1663
1664 if( solfound )
1665 {
1666 STPSOL* sol = SCIPStpHeurRecSolfromIdx(pool, solindex);
1667
1668 assert(pool->size >= 2);
1669 assert(sol != NULL);
1670
1671 SCIPdebugMessage("DA: rec found better solution with obj %f vs %f \n", sol->obj, ub);
1672
1673 if( SCIPisLE(scip, sol->obj, ub) )
1674 {
1675 BMScopyMemoryArray(result2, sol->soledges, nedges);
1676
1677 SCIP_CALL( SCIPStpHeurLocalRun(scip, graph, graph->cost, result2) );
1678 ub = graph_sol_getObj(graph->cost, result2, 0.0, nedges);
1679
1680 if( SCIPisLT(scip, ub, sol->obj) )
1681 SCIP_CALL( SCIPStpHeurRecAddToPool(scip, ub, result2, pool, &success) );
1682 }
1683 }
1684 }
1685 }
1686
1687 if( SCIPisLE(scip, ub, *upperbound) )
1688 {
1689 SCIPdebugMessage("DA: improved incumbent %f vs %f, return \n", ub, *upperbound);
1690
1691 *apsol = TRUE;
1692 *upperbound = ub;
1693 BMScopyMemoryArray(result1, result2, nedges);
1694 }
1695
1696 if( graph->stp_type != STP_MWCSP || !(*apsol) )
1697 return SCIP_OKAY;
1698
1699 #if 1
1700 SCIP_CALL( SCIPStpHeurRecExclude(scip, graph, result1, result2, pathedge, nodearrchar, &success) );
1701
1702 if( success )
1703 {
1704 BMScopyMemoryArray(result1, result2, nedges);
1705 *upperbound = graph_sol_getObj(graph->cost, result1, 0.0, nedges);
1706 SCIPdebugMessage("DA: afterLastExclusion %f \n", *upperbound);
1707 }
1708 #endif
1709
1710 return SCIP_OKAY;
1711 }
1712
1713
1714 /** try to improve both dual and primal solution */
1715 static
computePertubedSol(SCIP * scip,GRAPH * graph,GRAPH * transgraph,STPSOLPOOL * pool,PATH * vnoi,GNODE ** gnodearr,SCIP_Real * cost,SCIP_Real * costrev,SCIP_Real * bestcost,SCIP_Real * pathdist,int * state,int * vbase,int * pathedge,int * result,int * result2,int * transresult,STP_Bool * nodearrchar,SCIP_Real * upperbound,SCIP_Real * lpobjval,SCIP_Real * bestlpobjval,SCIP_Real * minpathcost,SCIP_Bool * apsol,SCIP_Real offset,int extnedges,int pertubation)1716 SCIP_RETCODE computePertubedSol(
1717 SCIP* scip,
1718 GRAPH* graph,
1719 GRAPH* transgraph,
1720 STPSOLPOOL* pool,
1721 PATH* vnoi,
1722 GNODE** gnodearr,
1723 SCIP_Real* cost,
1724 SCIP_Real* costrev,
1725 SCIP_Real* bestcost,
1726 SCIP_Real* pathdist,
1727 int* state,
1728 int* vbase,
1729 int* pathedge,
1730 int* result,
1731 int* result2,
1732 int* transresult,
1733 STP_Bool* nodearrchar,
1734 SCIP_Real* upperbound,
1735 SCIP_Real* lpobjval,
1736 SCIP_Real* bestlpobjval,
1737 SCIP_Real* minpathcost,
1738 SCIP_Bool* apsol,
1739 SCIP_Real offset,
1740 int extnedges,
1741 int pertubation
1742 )
1743 {
1744 SCIP_Real lb;
1745 int e;
1746 const int root = graph->source;
1747 const int nedges = graph->edges;
1748 const int transnedges = transgraph->edges;
1749
1750 graph_pc_2transcheck(graph);
1751
1752 /* pertubate the reduced cost? */
1753 if( graph->stp_type == STP_MWCSP )
1754 {
1755 SCIP_Real* transcost;
1756
1757 SCIP_CALL( SCIPallocBufferArray(scip, &transcost, transnedges) );
1758 BMScopyMemoryArray(transcost, transgraph->cost, transnedges);
1759
1760 /* result contains no valid solution?*/
1761 if( !(*apsol) )
1762 {
1763 /* compute new solution */
1764 SCIP_Real bound;
1765 SCIP_Bool success;
1766 SCIP_CALL( SCIPStpHeurAscendPruneRun(scip, NULL, graph, bestcost, result2, vbase, -1, nodearrchar, &success, FALSE) );
1767 assert(success);
1768
1769 SCIP_CALL( SCIPStpHeurLocalRun(scip, graph, graph->cost, result2) );
1770
1771 assert(graph_sol_valid(scip, graph, result2));
1772
1773 bound = graph_sol_getObj(graph->cost, result2, 0.0, nedges);
1774
1775 if( SCIPisLE(scip, bound, *upperbound) )
1776 {
1777 *upperbound = bound;
1778 *apsol = TRUE;
1779 BMScopyMemoryArray(result, result2, nedges);
1780 }
1781 pertubateEdgeCosts(scip, graph, transgraph, result2, nodearrchar, pertubation);
1782 }
1783 else
1784 {
1785 pertubateEdgeCosts(scip, graph, transgraph, result, nodearrchar, pertubation);
1786 }
1787
1788 /* todo use result as guiding solution? */
1789 SCIP_CALL( SCIPStpDualAscent(scip, transgraph, cost, pathdist, &lb, FALSE, FALSE, gnodearr, NULL, transresult, state, root, TRUE, -1.0, nodearrchar) );
1790
1791 BMScopyMemoryArray(transgraph->cost, transcost, transnedges);
1792
1793 SCIPfreeBufferArray(scip, &transcost);
1794 }
1795
1796 SCIP_CALL( computeDaSolPcMw(scip, graph, pool, vnoi, cost, pathdist, upperbound, result, result2, vbase, pathedge, nodearrchar, apsol) );
1797
1798 /* does result not contain a valid solution? */
1799 if( !(*apsol) )
1800 BMScopyMemoryArray(result, result2, nedges);
1801
1802 SCIPdebugMessage("DA: pertubated sol value %f \n", *upperbound);
1803
1804 /* compute guiding solution */
1805
1806 BMScopyMemoryArray(transresult, result, nedges);
1807
1808 for( e = nedges; e < transnedges; e++ )
1809 transresult[e] = UNKNOWN;
1810
1811 /* non-trivial solution */
1812 for( e = graph->outbeg[root]; e != EAT_LAST; e = graph->oeat[e] )
1813 if( Is_term(graph->term[graph->head[e]]) && result[e] != CONNECT )
1814 break;
1815
1816 if( e != EAT_LAST)
1817 {
1818 int k;
1819 for( e = graph->outbeg[root]; e != EAT_LAST; e = graph->oeat[e] )
1820 if( !Is_term(graph->term[graph->head[e]]) && result[e] == CONNECT )
1821 break;
1822
1823 assert(e != EAT_LAST);
1824 k = graph->head[e];
1825
1826 for( e = transgraph->outbeg[k]; e != EAT_LAST; e = transgraph->oeat[e] )
1827 {
1828 if( transgraph->head[e] == graph->knots )
1829 {
1830 transresult[e] = CONNECT;
1831 break;
1832 }
1833 }
1834 }
1835
1836 assert(!(*apsol) || graph_sol_valid(scip, graph, result));
1837 assert(graph_valid(transgraph));
1838 assert(root == transgraph->source);
1839
1840 SCIP_CALL( SCIPStpDualAscent(scip, transgraph, cost, pathdist, &lb, FALSE, FALSE, gnodearr, transresult, NULL, state, root, TRUE, -1.0, nodearrchar) );
1841
1842 lb += offset;
1843 *lpobjval = lb;
1844
1845 SCIP_CALL( computeDaSolPcMw(scip, graph, pool, vnoi, cost, pathdist, upperbound, result, result2, vbase, pathedge, nodearrchar, apsol) );
1846
1847 assert(!(*apsol) || graph_sol_valid(scip, graph, result));
1848
1849 if( SCIPisGE(scip, lb, *bestlpobjval) )
1850 {
1851 *bestlpobjval = lb;
1852 BMScopyMemoryArray(bestcost, cost, extnedges);
1853 }
1854
1855 /* the required reduced path cost to be surpassed */
1856 *minpathcost = *upperbound - *lpobjval;
1857
1858 return SCIP_OKAY;
1859 }
1860
1861 /** compute Voronoi region for dual-ascent elimination for PC/MW */
1862 static
computeTransVoronoi(SCIP * scip,GRAPH * transgraph,PATH * vnoi,const SCIP_Real * cost,SCIP_Real * costrev,SCIP_Real * pathdist,int * vbase,int * pathedge)1863 void computeTransVoronoi(
1864 SCIP* scip,
1865 GRAPH* transgraph,
1866 PATH* vnoi,
1867 const SCIP_Real* cost,
1868 SCIP_Real* costrev,
1869 SCIP_Real* pathdist,
1870 int* vbase,
1871 int* pathedge
1872 )
1873 {
1874 const int root = transgraph->source;
1875 const int transnnodes = transgraph->knots;
1876 const int transnedges = transgraph->edges;
1877
1878 for( int k = 0; k < transnnodes; k++ )
1879 transgraph->mark[k] = (transgraph->grad[k] > 0);
1880
1881 /* distance from root to all nodes */
1882 graph_path_execX(scip, transgraph, root, cost, pathdist, pathedge);
1883
1884 for( int k = 0; k < transnnodes; k++ )
1885 if( Is_term(transgraph->term[k]) )
1886 assert(SCIPisEQ(scip, pathdist[k], 0.0));
1887
1888 for( int e = 0; e < transnedges; e++ )
1889 costrev[e] = cost[flipedge(e)];
1890
1891 /* no paths should go back to the root */
1892 for( int e = transgraph->outbeg[root]; e != EAT_LAST; e = transgraph->oeat[e] )
1893 costrev[e] = FARAWAY;
1894
1895 /* build Voronoi diagram wrt incoming arcs */
1896 graph_voronoiTerms(scip, transgraph, costrev, vnoi, vbase, transgraph->path_heap, transgraph->path_state);
1897 }
1898
1899
1900 /** reduce SPG graph based on information from dual ascent and given upper bound */
1901 static
reduceSPG(SCIP * scip,GRAPH * graph,SCIP_Real * offset,STP_Bool * marked,STP_Bool * nodearrchar,const PATH * vnoi,const SCIP_Real * cost,const SCIP_Real * pathdist,const int * result,SCIP_Real minpathcost,int root,SCIP_Bool solgiven)1902 int reduceSPG(
1903 SCIP* scip, /**< SCIP data structure */
1904 GRAPH* graph, /**< graph data structure */
1905 SCIP_Real* offset, /**< offset pointer */
1906 STP_Bool* marked, /**< edge array to mark which (directed) edge can be removed */
1907 STP_Bool* nodearrchar, /**< node array storing solution vertices */
1908 const PATH* vnoi, /**< Voronoi data structure */
1909 const SCIP_Real* cost, /**< dual ascent costs */
1910 const SCIP_Real* pathdist, /**< distance array from shortest path calculations */
1911 const int* result, /**< sol int array */
1912 SCIP_Real minpathcost, /**< the required reduced path cost to be surpassed */
1913 int root, /**< the root */
1914 SCIP_Bool solgiven /**< is sol given? */
1915 )
1916 {
1917 int nfixed = 0;
1918 const int nedges = graph->edges;
1919 const int nnodes = graph->knots;
1920 const SCIP_Bool rpc = (graph->stp_type == STP_RPCSPG);
1921 const SCIP_Bool keepsol = (solgiven && SCIPisZero(scip, minpathcost));
1922
1923 if( solgiven )
1924 {
1925 for( int k = 0; k < nnodes; k++ )
1926 nodearrchar[k] = FALSE;
1927 for( int e = 0; e < nedges; e++ )
1928 {
1929 if( result[e] == CONNECT )
1930 {
1931 nodearrchar[graph->tail[e]] = TRUE;
1932 nodearrchar[graph->head[e]] = TRUE;
1933 }
1934 }
1935 }
1936
1937 /* RPC? If yes, restore original graph */
1938 if( rpc )
1939 graph_pc_2org(graph);
1940
1941 for( int k = 0; k < nnodes; k++ )
1942 {
1943 SCIP_Real redcost;
1944 int e;
1945
1946 if( !graph->mark[k] )
1947 continue;
1948
1949 assert(!Is_pterm(graph->term[k]));
1950
1951 redcost = pathdist[k] + vnoi[k].dist;
1952
1953 if( rpc && Is_term(graph->term[k]) && SCIPisGT(scip, redcost, minpathcost) && k != root )
1954 {
1955 (*offset) += graph->prize[k];
1956 nfixed += graph_pc_deleteTerm(scip, graph, k);
1957 continue;
1958 }
1959
1960 if( !Is_term(graph->term[k]) && !keepsol &&
1961 (SCIPisGT(scip, redcost, minpathcost) || (solgiven && SCIPisEQ(scip, redcost, minpathcost) && !nodearrchar[k])) )
1962 {
1963 nfixed += graph->grad[k];
1964 graph_knot_del(scip, graph, k, TRUE);
1965 continue;
1966 }
1967
1968 e = graph->outbeg[k];
1969 while( e != EAT_LAST )
1970 {
1971 const int head = graph->head[e];
1972 const int enext = graph->oeat[e];
1973
1974 /* for rpc no artificial terminal arcs should be deleted */
1975 if( (rpc && !graph->mark[head])
1976 || (keepsol && (result[e] == CONNECT || result[flipedge(e)] == CONNECT)) )
1977 {
1978 e = enext;
1979 continue;
1980 }
1981
1982 redcost = pathdist[k] + cost[e] + vnoi[head].dist;
1983
1984 if( SCIPisGT(scip, redcost, minpathcost)
1985 || (solgiven && SCIPisEQ(scip, redcost, minpathcost) && result[e] != CONNECT && result[flipedge(e)] != CONNECT) )
1986 {
1987 if( marked[flipedge(e)] )
1988 {
1989 graph_edge_del(scip, graph, e, TRUE);
1990 nfixed++;
1991 }
1992 else
1993 {
1994 marked[e] = TRUE;
1995 }
1996 }
1997 e = enext;
1998 }
1999 }
2000
2001 for( int k = 0; k < nnodes; k++ )
2002 if( graph->grad[k] == 0 && k != root )
2003 graph->mark[k] = FALSE;
2004
2005 return nfixed;
2006 }
2007
2008 /** reduce PCSTP or MWCS graph based on information from dual ascent and given upper bound */
2009 static
reducePcMw(SCIP * scip,GRAPH * graph,GRAPH * transgraph,const PATH * vnoi,const SCIP_Real * cost,const SCIP_Real * pathdist,SCIP_Real minpathcost,const int * result,STP_Bool * marked,STP_Bool * nodearrchar,SCIP_Bool solgiven)2010 int reducePcMw(
2011 SCIP* scip, /**< SCIP data structure */
2012 GRAPH* graph, /**< graph data structure */
2013 GRAPH* transgraph, /**< graph data structure */
2014 const PATH* vnoi, /**< Voronoi data structure */
2015 const SCIP_Real* cost, /**< dual ascent costs */
2016 const SCIP_Real* pathdist, /**< distance array from shortest path calculations */
2017 SCIP_Real minpathcost, /**< the required reduced path cost to be surpassed */
2018 const int* result, /**< sol int array */
2019 STP_Bool* marked, /**< edge array to mark which (directed) edge can be removed */
2020 STP_Bool* nodearrchar, /**< node array storing solution vertices */
2021 SCIP_Bool solgiven /**< is sol given? */
2022 )
2023 {
2024 int nnodes;
2025 int nfixed;
2026 SCIP_Real tmpcost;
2027 SCIP_Bool keepsol = FALSE;
2028
2029 assert(scip != NULL);
2030 assert(graph != NULL);
2031 assert(pathdist != NULL);
2032 assert(result != NULL);
2033 assert(cost != NULL);
2034 assert(vnoi != NULL);
2035
2036 assert(SCIPisGE(scip, minpathcost, 0.0));
2037 assert(!solgiven || graph_sol_valid(scip, graph, result));
2038
2039 if( minpathcost < 0.0 )
2040 minpathcost = 0.0;
2041
2042 nfixed = 0;
2043 nnodes = graph->knots;
2044
2045 graph_pc_2orgcheck(graph);
2046
2047 if( solgiven )
2048 {
2049 const int nedges = graph->edges;
2050
2051 for( int k = 0; k < nnodes; k++ )
2052 nodearrchar[k] = FALSE;
2053
2054 for( int e = 0; e < nedges; e++ )
2055 {
2056 if( result[e] == CONNECT )
2057 {
2058 assert(graph->oeat[e] != EAT_FREE);
2059
2060 nodearrchar[graph->head[e]] = TRUE;
2061 nodearrchar[graph->tail[e]] = TRUE;
2062 }
2063 }
2064 if( SCIPisZero(scip, minpathcost) )
2065 keepsol = TRUE;
2066 }
2067
2068 /* try to eliminate vertices and edges */
2069 for( int k = 0; k < nnodes; k++ )
2070 {
2071 if( !graph->mark[k] )
2072 continue;
2073
2074 assert(!Is_pterm(graph->term[k]));
2075
2076 if( Is_term(graph->term[k]) )
2077 {
2078 int e = graph->outbeg[k];
2079 while( e != EAT_LAST )
2080 {
2081 const int etmp = graph->oeat[e];
2082 tmpcost = pathdist[k] + cost[e] + vnoi[graph->head[e]].dist;
2083
2084 if( graph->mark[graph->head[e]] &&
2085 ((SCIPisGT(scip, tmpcost, minpathcost) && (!keepsol || (result[e] != CONNECT && result[flipedge(e)] != CONNECT))) ||
2086 (solgiven && tmpcost >= minpathcost && result[e] != CONNECT && result[flipedge(e)] != CONNECT)) )
2087 {
2088 if( marked[flipedge(e)] )
2089 {
2090 assert(!Is_pterm(graph->term[graph->head[e]]));
2091
2092 graph_edge_del(scip, graph, e, TRUE);
2093 graph_edge_del(scip, transgraph, e, FALSE);
2094 nfixed++;
2095 }
2096 else
2097 {
2098 marked[e] = TRUE;
2099 }
2100 }
2101 e = etmp;
2102 }
2103 continue;
2104 }
2105
2106 tmpcost = pathdist[k] + vnoi[k].dist;
2107
2108 if( (SCIPisGT(scip, tmpcost, minpathcost) && !keepsol) ||
2109 (solgiven && tmpcost >= minpathcost && !nodearrchar[k]))
2110 {
2111 while( transgraph->outbeg[k] != EAT_LAST )
2112 {
2113 const int e = transgraph->outbeg[k];
2114
2115 graph_edge_del(scip, transgraph, e, FALSE);
2116 graph_edge_del(scip, graph, e, TRUE);
2117 nfixed++;
2118 }
2119 assert(graph->outbeg[k] == EAT_LAST);
2120 }
2121 else
2122 {
2123 int e = transgraph->outbeg[k];
2124 while( e != EAT_LAST )
2125 {
2126 const int etmp = transgraph->oeat[e];
2127 tmpcost = pathdist[k] + cost[e] + vnoi[transgraph->head[e]].dist;
2128
2129 if( (SCIPisGT(scip, tmpcost, minpathcost) && (!keepsol || (result[e] != CONNECT && result[flipedge(e)] != CONNECT))) ||
2130 (solgiven && tmpcost >= minpathcost && result[e] != CONNECT && result[flipedge(e)] != CONNECT) )
2131 {
2132 if( marked[flipedge(e)] )
2133 {
2134 graph_edge_del(scip, graph, e, TRUE);
2135 graph_edge_del(scip, transgraph, e, FALSE);
2136
2137 nfixed++;
2138 }
2139 else
2140 {
2141 marked[e] = TRUE;
2142 }
2143 }
2144 e = etmp;
2145 }
2146 }
2147 }
2148 SCIPdebugMessage("DA: eliminations %d \n", nfixed);
2149
2150 for( int k = 0; k < nnodes; k++ )
2151 if( graph->grad[k] == 0 && k != graph->source )
2152 graph->mark[k] = FALSE;
2153
2154 return nfixed;
2155 }
2156
2157 /** try to run reduction method for best known reduced costs (if they are valid) */
2158 static
reducePcMwTryBest(SCIP * scip,GRAPH * graph,GRAPH * transgraph,PATH * vnoi,const SCIP_Real * cost,SCIP_Real * costrev,SCIP_Real * bestcost,SCIP_Real * pathdist,SCIP_Real * upperbound,SCIP_Real * lpobjval,SCIP_Real * bestlpobjval,SCIP_Real * minpathcost,SCIP_Real oldupperbound,const int * result,int * vbase,int * state,int * pathedge,STP_Bool * marked,STP_Bool * nodearrchar,SCIP_Bool * solgiven,int extnedges)2159 int reducePcMwTryBest(
2160 SCIP* scip, /**< SCIP data structure */
2161 GRAPH* graph, /**< graph data structure */
2162 GRAPH* transgraph, /**< graph data structure */
2163 PATH* vnoi, /**< Voronoi data structure */
2164 const SCIP_Real* cost, /**< dual ascent costs */
2165 SCIP_Real* costrev, /**< SCIP_Real array */
2166 SCIP_Real* bestcost, /**< best dual ascent costs */
2167 SCIP_Real* pathdist, /**< distance array from shortest path calculations */
2168 SCIP_Real* upperbound, /**< upper bound */
2169 SCIP_Real* lpobjval, /**< reduced cost value */
2170 SCIP_Real* bestlpobjval, /**< best reduced cost value */
2171 SCIP_Real* minpathcost, /**< the required reduced path cost to be surpassed */
2172 SCIP_Real oldupperbound, /**< old upper bound */
2173 const int* result, /**< sol int array */
2174 int* vbase, /**< array for Voronoi bases */
2175 int* state, /**< int 4 * nnodes array for internal computations */
2176 int* pathedge, /**< array for predecessor edge on a path */
2177 STP_Bool* marked, /**< edge array to mark which (directed) edge can be removed */
2178 STP_Bool* nodearrchar, /**< node array storing solution vertices */
2179 SCIP_Bool* solgiven, /**< is sol given? */
2180 int extnedges /**< number of edges for transgraph */
2181 )
2182 {
2183 const SCIP_Bool dualisvalid = dualCostIsValid(scip, transgraph, bestcost, state, nodearrchar);
2184
2185 if( !dualisvalid )
2186 {
2187 *bestlpobjval = *lpobjval;
2188 BMScopyMemoryArray(bestcost, cost, extnedges);
2189 }
2190
2191 /* has upperbound and changed, and is best reduced cost valid and different from cost? */
2192 if( SCIPisGT(scip, *bestlpobjval, *lpobjval) && SCIPisLT(scip, *upperbound, oldupperbound) )
2193 {
2194 *solgiven = *solgiven && graph_sol_unreduced(scip, graph, result);
2195
2196 assert(!(*solgiven) || graph_sol_valid(scip, graph, result));
2197
2198 *minpathcost = *upperbound - *bestlpobjval;
2199 assert(SCIPisGE(scip, *minpathcost, 0.0));
2200
2201 computeTransVoronoi(scip, transgraph, vnoi, bestcost, costrev, pathdist, vbase, pathedge);
2202
2203 return reducePcMw(scip, graph, transgraph, vnoi, bestcost, pathdist, *minpathcost, result, marked, nodearrchar, *solgiven);
2204 }
2205 return 0;
2206 }
2207
2208
2209 /** check (directed) edge */
2210 static
reduceCheckEdge(SCIP * scip,const GRAPH * graph,int root,const SCIP_Real * redcost,const SCIP_Real * pathdist,const PATH * vnoi,SCIP_Real cutoff,int edge,SCIP_Bool equality,int * edgestack,SCIP_Bool * deletable,unsigned int * eqstack,int * eqstack_size,SCIP_Bool * eqmark)2211 SCIP_RETCODE reduceCheckEdge(
2212 SCIP* scip, /**< SCIP data structure */
2213 const GRAPH* graph, /**< graph data structure */
2214 int root, /**< graph root from dual ascent */
2215 const SCIP_Real* redcost, /**< reduced costs */
2216 const SCIP_Real* pathdist, /**< shortest path distances */
2217 const PATH* vnoi, /**< Voronoi paths */
2218 SCIP_Real cutoff, /**< cutoff value */
2219 int edge, /**< (directed) edge to be checked */
2220 SCIP_Bool equality, /**< allow equality? */
2221 int* edgestack, /**< array of size nodes for internal computations */
2222 SCIP_Bool* deletable, /**< is edge deletable? */
2223 unsigned int* eqstack, /**< stores edges that were used for equality comparison */
2224 int* eqstack_size, /**< pointer to size of eqstack */
2225 SCIP_Bool* eqmark /**< marks edges that were used for equality comparison */
2226 )
2227 {
2228 const int head = graph->head[edge];
2229 const int tail = graph->tail[edge];
2230 const int maxgrad = (graph->edges > STP_DAEX_EDGELIMIT) ? STP_DAEX_MINGRAD : STP_DAEX_MAXGRAD;
2231 SCIP_Real edgebound = redcost[edge] + pathdist[tail] + vnoi[head].dist;
2232
2233 assert(scip != NULL);
2234 assert(graph != NULL);
2235 assert(redcost != NULL);
2236 assert(pathdist != NULL);
2237 assert(vnoi != NULL);
2238 assert(edge >= 0 && edge < graph->edges);
2239
2240 *deletable = FALSE;
2241
2242 /* trivial rule-out? */
2243 if( SCIPisGT(scip, edgebound, cutoff) || (equality && SCIPisEQ(scip, edgebound, cutoff)) || head == root )
2244 {
2245 *deletable = TRUE;
2246 }
2247 else if( (graph->grad[tail] <= maxgrad && !Is_term(graph->term[tail]))
2248 || (graph->grad[head] <= maxgrad && !Is_term(graph->term[head])) )
2249 {
2250 SCIP_Real treecosts[STP_DAEX_MAXDFSDEPTH + 1];
2251 int treeedges[STP_DAEX_MAXDFSDEPTH + 1];
2252 SCIP_Real basebottlenecks[3];
2253 int* nodepos;
2254 int* ancestormark;
2255
2256 const int nnodes = graph->knots;
2257 const int maxdfsdepth = (graph->edges > STP_DAEX_EDGELIMIT) ? STP_DAEX_MINDFSDEPTH : STP_DAEX_MAXDFSDEPTH;
2258
2259 /* allocate clean arrays */
2260 SCIP_CALL( SCIPallocCleanBufferArray(scip, &nodepos, nnodes) );
2261 SCIP_CALL( SCIPallocCleanBufferArray(scip, &ancestormark, (MAX(graph->edges, graph->orgedges) / 2)) );
2262
2263 basebottlenecks[0] = graph->cost[edge];
2264 #ifndef NEBUG
2265 basebottlenecks[1] = -1.0;
2266 basebottlenecks[2] = -1.0;
2267
2268 for( int i = 0; i < STP_DAEX_MAXDFSDEPTH + 1; i++ )
2269 {
2270 treecosts[i] = -1.0;
2271 treeedges[i] = -1;
2272 }
2273 #endif
2274
2275 markAncestors(graph, edge, ancestormark);
2276
2277 /* check whether to extend from head */
2278 if( !Is_term(graph->term[head]) && graph->grad[head] <= maxgrad )
2279 {
2280 const SCIP_Real treecostoffset = redcost[edge] + pathdist[tail];
2281 SCIP_Real minbound = FARAWAY;
2282 SCIP_Real treecost = treecostoffset;
2283 unsigned int rounds = 0;
2284 int dfsdepth = 0;
2285 int nadded_edges = 0;
2286 SCIP_Bool stopped = FALSE;
2287
2288 nodepos[tail] = nnodes + 1;
2289 nodepos[head] = nnodes + 4;
2290
2291 for( int e = graph->outbeg[head]; e != EAT_LAST; e = graph->oeat[e] )
2292 if( graph->head[e] != tail )
2293 edgestack[nadded_edges++] = e;
2294
2295 /* limited DFS */
2296 while( nadded_edges > 0 )
2297 {
2298 const int curredge = edgestack[--nadded_edges];
2299 const int currhead = graph->head[curredge];
2300
2301 /* subtree already processed? */
2302 if( nodepos[currhead] )
2303 {
2304 assert(dfsdepth >= 1 && treeedges[dfsdepth - 1] == curredge);
2305
2306 treecost = shortenSubtree(scip, graph, redcost, treeedges, treecost, treecostoffset, dfsdepth,
2307 currhead, nodepos, ancestormark, &rounds);
2308
2309 dfsdepth--;
2310 }
2311 else
2312 {
2313 const SCIP_Real extendedcost = treecost + redcost[curredge] + vnoi[currhead].dist;
2314
2315 if( ruleOutSubtree(scip, graph, treecosts, basebottlenecks, nodepos, treeedges, cutoff, extendedcost, dfsdepth, curredge, equality,
2316 ancestormark, eqstack, eqstack_size, eqmark) )
2317 continue;
2318
2319 if( truncateSubtree(graph, extendedcost, root, currhead, maxgrad, dfsdepth, maxdfsdepth, &minbound, &stopped) )
2320 {
2321 /* stopped and no further improvement of bound possible? */
2322 if( stopped && minbound <= edgebound )
2323 break;
2324
2325 continue;
2326 }
2327
2328 nadded_edges = extendSubtreeHead(scip, graph, redcost, curredge, currhead, dfsdepth, nadded_edges, &treecost, treecosts, nodepos,
2329 treeedges, edgestack, ancestormark);
2330 dfsdepth++;
2331 }
2332 } /* DFS loop */
2333
2334 assert(stopped || minbound == FARAWAY);
2335 assert(SCIPisGE(scip, minbound, redcost[edge] + pathdist[tail] + vnoi[head].dist));
2336
2337 finalizeSubtree(graph, graph->head, treeedges, dfsdepth, stopped, minbound, &edgebound, deletable, nodepos, ancestormark);
2338 } /* extend from head */
2339
2340 /* check whether to extend from tail */
2341 if( !(*deletable) && !Is_term(graph->term[tail]) && graph->grad[tail] <= maxgrad )
2342 {
2343 const SCIP_Real treecostoffset = edgebound - pathdist[tail];
2344 SCIP_Real minbound = FARAWAY;
2345 SCIP_Real treecost = treecostoffset;
2346 int dfsdepth = 0;
2347 int nadded_edges = 0;
2348 unsigned int rounds = 0;
2349 SCIP_Bool stopped = FALSE;
2350
2351 nodepos[head] = nnodes + 1;
2352 nodepos[tail] = nnodes + 4;
2353
2354 for( int e = graph->inpbeg[tail]; e != EAT_LAST; e = graph->ieat[e] )
2355 if( graph->tail[e] != head )
2356 edgestack[nadded_edges++] = e;
2357
2358 /* limited DFS */
2359 while( nadded_edges > 0 )
2360 {
2361 const int curredge = edgestack[--nadded_edges];
2362 const int currtail = graph->tail[curredge];
2363
2364 /* subtree already processed? */
2365 if( nodepos[currtail] )
2366 {
2367 assert(dfsdepth >= 1 && treeedges[dfsdepth - 1] == curredge);
2368
2369 treecost = shortenSubtree(scip, graph, redcost, treeedges, treecost, treecostoffset, dfsdepth,
2370 currtail, nodepos, ancestormark, &rounds);
2371
2372 dfsdepth--;
2373 }
2374 else
2375 {
2376 const SCIP_Real extendedcost = treecost + redcost[curredge] + pathdist[currtail];
2377
2378 if( ruleOutSubtree(scip, graph, treecosts, basebottlenecks, nodepos, treeedges, cutoff, extendedcost, dfsdepth, flipedge((unsigned)curredge), equality,
2379 ancestormark, eqstack, eqstack_size, eqmark) )
2380 continue;
2381
2382 if( truncateSubtree(graph, extendedcost, -1, currtail, maxgrad, dfsdepth, maxdfsdepth, &minbound, &stopped) )
2383 {
2384 if( stopped )
2385 break;
2386 continue;
2387 }
2388
2389 nadded_edges = extendSubtreeTail(graph, redcost, curredge, currtail, dfsdepth, nadded_edges, &treecost, treecosts, nodepos, treeedges,
2390 edgestack, ancestormark);
2391 dfsdepth++;
2392 }
2393 } /* DFS loop */
2394
2395 assert(stopped || minbound == FARAWAY);
2396 assert(SCIPisGE(scip, minbound, redcost[edge] + pathdist[tail] + vnoi[head].dist));
2397
2398 finalizeSubtree(graph, graph->tail, treeedges, dfsdepth, stopped, minbound, &edgebound, deletable, nodepos, ancestormark);
2399 } /* extend from tail */
2400
2401 /* clean arrays */
2402 nodepos[head] = 0;
2403 nodepos[tail] = 0;
2404 unmarkAncestors(graph, edge, ancestormark);
2405
2406 /* free memory */
2407 for( int i = 0; i < nnodes; i++ )
2408 assert(nodepos[i] == 0);
2409 for( int i = 0; i < MAX(graph->edges, graph->orgedges) / 2; i++ )
2410 assert(ancestormark[i] == 0);
2411
2412 SCIPfreeCleanBufferArray(scip, &ancestormark);
2413 SCIPfreeCleanBufferArray(scip, &nodepos);
2414 }
2415
2416 return SCIP_OKAY;
2417 }
2418
2419
2420 /** todo don't use this function, but check for conflicts in misc_stp.c and handle them */
reduce_deleteConflictEdges(SCIP * scip,GRAPH * g)2421 SCIP_RETCODE reduce_deleteConflictEdges(
2422 SCIP* scip, /**< SCIP data structure */
2423 GRAPH* g /**< the graph */
2424 )
2425 {
2426 int* edgemark;
2427 const int nedges = g->edges;
2428 const int nancestors = MAX(g->edges, g->orgedges) / 2;
2429
2430 assert(scip != NULL && g != NULL);
2431 assert(g->ancestors != NULL);
2432
2433 SCIP_CALL( SCIPallocBufferArray(scip, &edgemark, nancestors) );
2434
2435 for( int e = 0; e < nancestors; e++ )
2436 edgemark[e] = FALSE;
2437
2438 for( int e = 0; e < nedges; e += 2 )
2439 {
2440 SCIP_Bool conflict;
2441
2442 if( g->oeat[e] == EAT_FREE )
2443 continue;
2444
2445 conflict = markAncestorsConflict(g, e, edgemark);
2446
2447 unmarkAncestorsConflict(g, e, edgemark);
2448 if( conflict )
2449 {
2450 conflict = markAncestorsConflict(g, e + 1, edgemark);
2451 unmarkAncestorsConflict(g, e + 1, edgemark);
2452
2453 if( conflict )
2454 graph_edge_del(scip, g, e, TRUE);
2455 }
2456 }
2457
2458 SCIPfreeBufferArray(scip, &edgemark);
2459
2460 return SCIP_OKAY;
2461 }
2462
2463 /** check 3-tree */
reduce_check3Tree(SCIP * scip,const GRAPH * graph,int root,const SCIP_Real * redcost,const SCIP_Real * pathdist,const PATH * vnoi,const int * vbase,SCIP_Real cutoff,const int * outedges,int inedge,int * edgestack,SCIP_Real * treebound,SCIP_Bool * ruleout,unsigned int * eqstack,int * eqstack_size,SCIP_Bool * eqmark)2464 SCIP_RETCODE reduce_check3Tree(
2465 SCIP* scip, /**< SCIP data structure */
2466 const GRAPH* graph, /**< graph data structure */
2467 int root, /**< graph root from dual ascent */
2468 const SCIP_Real* redcost, /**< reduced costs */
2469 const SCIP_Real* pathdist, /**< shortest path distances */
2470 const PATH* vnoi, /**< Voronoi paths */
2471 const int* vbase, /**< bases to Voronoi paths */
2472 SCIP_Real cutoff, /**< cutoff value */
2473 const int* outedges, /**< two outgoing edge */
2474 int inedge, /**< incoming edge */
2475 int* edgestack, /**< array of size nodes for internal computations */
2476 SCIP_Real* treebound, /**< to store a lower bound for the tree */
2477 SCIP_Bool* ruleout, /**< could tree be ruled out? */
2478 unsigned int* eqstack, /**< stores edges that were used for equality comparison */
2479 int* eqstack_size, /**< pointer to size of eqstack */
2480 SCIP_Bool* eqmark /**< marks edges that were used for equality comparison */
2481 )
2482 {
2483 #ifndef NDEBUG
2484 const SCIP_Real orgtreebound = *treebound;
2485 #endif
2486 assert(scip != NULL);
2487 assert(graph != NULL);
2488 assert(redcost != NULL);
2489 assert(ruleout != NULL);
2490 assert(pathdist != NULL);
2491 assert(vnoi != NULL);
2492 assert(vbase != NULL);
2493 assert(outedges != NULL);
2494 assert(inedge >= 0 && outedges[0] >= 0 && outedges[1] >= 0);
2495 assert(inedge < graph->edges && outedges[0] < graph->edges && outedges[1] < graph->edges);
2496
2497 /* trivial rule-out? */
2498 if( SCIPisGT(scip, *treebound, cutoff) )
2499 {
2500 *ruleout = TRUE;
2501 }
2502 else
2503 {
2504 SCIP_Real treecosts[STP_DAEX_MAXDFSDEPTH + 1];
2505 int treeedges[STP_DAEX_MAXDFSDEPTH + 1];
2506 SCIP_Real basebottlenecks[3];
2507 const SCIP_Real basecost = redcost[inedge] + redcost[outedges[0]] + redcost[outedges[1]] + pathdist[graph->tail[inedge]];
2508 int* nodepos;
2509 int* ancestormark;
2510 const int nnodes = graph->knots;
2511 const int innode = graph->tail[inedge];
2512 const int basenode = graph->head[inedge];
2513 const int maxdfsdepth = (graph->edges > STP_DAEX_EDGELIMIT) ? STP_DAEX_MINDFSDEPTH : STP_DAEX_MAXDFSDEPTH;
2514 const int maxgrad = (graph->edges > STP_DAEX_EDGELIMIT) ? STP_DAEX_MINGRAD : STP_DAEX_MAXGRAD;
2515
2516 assert(basenode == graph->tail[outedges[0]] && basenode == graph->tail[outedges[1]]);
2517
2518 /* allocate clean array */
2519 SCIP_CALL(SCIPallocCleanBufferArray(scip, &nodepos, nnodes));
2520 SCIP_CALL( SCIPallocCleanBufferArray(scip, &ancestormark, (MAX(graph->edges, graph->orgedges) / 2)) );
2521
2522 nodepos[basenode] = nnodes + 1; /* basenode */
2523 nodepos[innode] = nnodes + 2; /* innode */
2524
2525 #ifndef NDEBUG
2526 for( int i = 0; i < STP_DAEX_MAXDFSDEPTH + 1; i++ )
2527 {
2528 treecosts[i] = -1.0;
2529 treeedges[i] = -1;
2530 }
2531 #endif
2532
2533 *ruleout = FALSE;
2534
2535 /* can we rule out subgraph beforehand? */
2536 if( markAncestorsConflict(graph, inedge, ancestormark)
2537 || markAncestorsConflict(graph, outedges[0], ancestormark)
2538 || markAncestorsConflict(graph, outedges[1], ancestormark) )
2539 {
2540 *ruleout = TRUE;
2541 }
2542
2543 /* main loop: extend tree and try to rule it out */
2544 for( int i = 0; i < 2 && !(*ruleout); i++ )
2545 {
2546 SCIP_Real minbound = FARAWAY;
2547 SCIP_Real treecost = basecost;
2548 const int startedge = outedges[i];
2549 const int costartedge = outedges[(i + 1) % 2];
2550 const int startnode = graph->head[startedge];
2551 const int costartnode = graph->head[costartedge];
2552 int dfsdepth = 0;
2553 int nadded_edges = 0;
2554 SCIP_Bool stopped = FALSE;
2555 unsigned int rounds = 0;
2556
2557 /* try to extend the tree from startnode */
2558
2559 assert(startnode != root && costartnode != root);
2560
2561 /* for basenode */
2562 basebottlenecks[0] = graph->cost[startedge];
2563
2564 /* for innode */
2565 basebottlenecks[1] = MAX(graph->cost[startedge], graph->cost[inedge]);
2566
2567 /* for costartnode */
2568 basebottlenecks[2] = MAX(graph->cost[startedge], graph->cost[costartedge]);
2569
2570 nodepos[costartnode] = nnodes + 3;
2571 nodepos[startnode] = nnodes + 4;
2572
2573 /* can we rule out entire subtree already? */
2574 if( ruleOutSubtree(scip, graph, treecosts, basebottlenecks, nodepos, treeedges, FARAWAY, 0.0, 0, startedge, FALSE,
2575 NULL, NULL, NULL, NULL) )
2576 {
2577 *ruleout = TRUE;
2578 break;
2579 }
2580
2581 /* cannot extend over terminals or vertices of high degree */
2582 if( Is_term(graph->term[startnode]) || graph->grad[startnode] > maxgrad )
2583 continue;
2584
2585 assert(nodepos[startnode]);
2586
2587 for( int e = graph->outbeg[startnode]; e != EAT_LAST; e = graph->oeat[e] )
2588 if( !nodepos[graph->head[e]] )
2589 edgestack[nadded_edges++] = e;
2590
2591 /* limited DFS starting from startnode */
2592 while ( nadded_edges > 0 )
2593 {
2594 const int curredge = edgestack[--nadded_edges];
2595 const int currhead = graph->head[curredge];
2596
2597 /* subtree already processed? */
2598 if( nodepos[currhead] )
2599 {
2600 assert(dfsdepth >= 1 && treeedges[dfsdepth - 1] == curredge);
2601
2602 treecost = shortenSubtree(scip, graph, redcost, treeedges, treecost, basecost, dfsdepth,
2603 currhead, nodepos, ancestormark, &rounds);
2604
2605 dfsdepth--;
2606 }
2607 else
2608 {
2609 SCIP_Real extendedcost = treecost + redcost[curredge];
2610
2611 if( vbase[costartnode] == vbase[currhead] )
2612 {
2613 /* also covers the case that leaf is a terminal */
2614 const SCIP_Real costartfar = vnoi[costartnode + nnodes].dist;
2615 const SCIP_Real currentfar = vnoi[currhead + nnodes].dist;
2616
2617 assert(vbase[costartnode + nnodes] >= 0 || costartfar == FARAWAY);
2618 assert(vbase[currhead + nnodes] >= 0 || currentfar == FARAWAY);
2619 extendedcost += MIN(costartfar + vnoi[currhead].dist, vnoi[costartnode].dist + currentfar);
2620 }
2621 else
2622 extendedcost += vnoi[costartnode].dist + vnoi[currhead].dist;
2623
2624 /* can we rule out subtree? */
2625 if( ruleOutSubtree(scip, graph, treecosts, basebottlenecks, nodepos, treeedges, cutoff, extendedcost, dfsdepth, curredge, FALSE,
2626 ancestormark, eqstack, eqstack_size, eqmark) )
2627 continue;
2628
2629 /* do we need to stop extension? */
2630 if( truncateSubtree(graph, extendedcost, root, currhead, maxgrad, dfsdepth, maxdfsdepth, &minbound, &stopped) )
2631 {
2632 if( stopped && minbound <= (*treebound) )
2633 break;
2634
2635 continue;
2636 }
2637
2638 nadded_edges = extendSubtreeHead(scip, graph, redcost, curredge, currhead, dfsdepth, nadded_edges, &treecost, treecosts, nodepos,
2639 treeedges, edgestack, ancestormark);
2640 dfsdepth++;
2641 }
2642 } /* DFS loop */
2643
2644 assert(stopped || minbound == FARAWAY);
2645 #ifndef NDEBUG
2646 assert(SCIPisGE(scip, minbound, orgtreebound));
2647 #endif
2648
2649 finalizeSubtree(graph, graph->head, treeedges, dfsdepth, stopped, minbound, treebound, ruleout, nodepos, ancestormark);
2650 } /* main loop for outgoing edges */
2651
2652 #ifndef NDEBUG
2653 assert(*treebound >= orgtreebound);
2654 #endif
2655
2656 if( !(*ruleout) && !Is_term(graph->term[innode]) && graph->grad[innode] <= maxgrad )
2657 {
2658 /* move down the incoming edge */
2659 const SCIP_Real treecostoffset = *treebound - pathdist[innode];
2660 SCIP_Real minbound = FARAWAY;
2661 SCIP_Real treecost = treecostoffset;
2662 int dfsdepth = 0;
2663 int nadded_edges = 0;
2664 SCIP_Bool stopped = FALSE;
2665 unsigned int rounds = 0;
2666
2667 assert(treecost >= 0.0);
2668 assert(nodepos[innode] == nnodes + 2);
2669 assert(nodepos[basenode] == nnodes + 1);
2670
2671 nodepos[graph->head[outedges[0]]] = nnodes + 2;
2672 nodepos[graph->head[outedges[1]]] = nnodes + 3;
2673 nodepos[innode] = nnodes + 4;
2674
2675 /* for basenode */
2676 basebottlenecks[0] = graph->cost[inedge];
2677
2678 /* for outedges[0] */
2679 basebottlenecks[1] = MAX(graph->cost[outedges[0]], graph->cost[inedge]);
2680
2681 /* for outedges[1] */
2682 basebottlenecks[2] = MAX(graph->cost[outedges[1]], graph->cost[inedge]);
2683
2684 /* can we rule out entire subtree already? */
2685 if( ruleOutSubtree(scip, graph, treecosts, basebottlenecks, nodepos, treeedges, FARAWAY, 0.0, 0, flipedge(inedge), FALSE,
2686 NULL, NULL, NULL, NULL) )
2687 {
2688 *ruleout = TRUE;
2689 }
2690 else
2691 {
2692 for( int e = graph->inpbeg[innode]; e != EAT_LAST; e = graph->ieat[e] )
2693 if( !nodepos[graph->tail[e]] )
2694 edgestack[nadded_edges++] = e;
2695 }
2696
2697 /* limited DFS starting from inedge */
2698 while( nadded_edges > 0 )
2699 {
2700 const int curredge = edgestack[--nadded_edges];
2701 const int currtail = graph->tail[curredge];
2702
2703 /* subtree already processed? */
2704 if( nodepos[currtail] )
2705 {
2706 assert(dfsdepth >= 1 && treeedges[dfsdepth - 1] == curredge);
2707
2708 treecost = shortenSubtree(scip, graph, redcost, treeedges, treecost, treecostoffset, dfsdepth,
2709 currtail, nodepos, ancestormark, &rounds);
2710
2711 dfsdepth--;
2712 }
2713 else
2714 {
2715 const SCIP_Real extendedcost = treecost + redcost[curredge] + pathdist[currtail];
2716
2717 if( ruleOutSubtree(scip, graph, treecosts, basebottlenecks, nodepos, treeedges, cutoff, extendedcost, dfsdepth, flipedge((unsigned)curredge), FALSE,
2718 ancestormark, eqstack, eqstack_size, eqmark) )
2719 continue;
2720
2721 if( truncateSubtree(graph, extendedcost, -1, currtail, maxgrad, dfsdepth, maxdfsdepth, &minbound, &stopped) )
2722 {
2723 if( stopped && minbound <= (*treebound) )
2724 break;
2725
2726 continue;
2727 }
2728
2729 nadded_edges = extendSubtreeTail(graph, redcost, curredge, currtail, dfsdepth, nadded_edges, &treecost, treecosts, nodepos,
2730 treeedges, edgestack, ancestormark);
2731 dfsdepth++;
2732 }
2733 } /* DFS loop */
2734
2735 #ifndef NDEBUG
2736 assert(stopped || minbound == FARAWAY);
2737 assert(SCIPisGE(scip, minbound, orgtreebound));
2738 #endif
2739
2740 finalizeSubtree(graph, graph->tail, treeedges, dfsdepth, stopped, minbound, treebound, ruleout, nodepos, ancestormark);
2741 } /* inedge */
2742
2743 /* clean nodepos array */
2744 nodepos[basenode] = 0;
2745 nodepos[innode] = 0;
2746
2747 unmarkAncestorsConflict(graph, inedge, ancestormark);
2748
2749 for( int i = 0; i < 2; i++ )
2750 {
2751 nodepos[graph->head[outedges[i]]] = 0;
2752 unmarkAncestorsConflict(graph, outedges[i], ancestormark);
2753 }
2754
2755 #ifndef NDEBUG
2756 assert(*treebound >= orgtreebound);
2757
2758 for( int i = 0; i < nnodes; i++ )
2759 assert(nodepos[i] == 0);
2760
2761 for( int i = 0; i < MAX(graph->edges, graph->orgedges) / 2; i++ )
2762 assert(ancestormark[i] == 0);
2763 #endif
2764
2765 SCIPfreeCleanBufferArray(scip, &ancestormark);
2766 SCIPfreeCleanBufferArray(scip, &nodepos);
2767 }
2768
2769 return SCIP_OKAY;
2770 }
2771
2772
2773 /** reduce SPG graph based on reduced cost information and given upper bound */
reduce_extendedEdge(SCIP * scip,GRAPH * graph,const PATH * vnoi,const SCIP_Real * cost,const SCIP_Real * pathdist,const int * result,SCIP_Real minpathcost,int root,int * nodearr,STP_Bool * marked)2774 int reduce_extendedEdge(
2775 SCIP* scip, /**< SCIP data structure */
2776 GRAPH* graph, /**< graph data structure */
2777 const PATH* vnoi, /**< Voronoi data structure */
2778 const SCIP_Real* cost, /**< dual ascent costs */
2779 const SCIP_Real* pathdist, /**< distance array from shortest path calculations */
2780 const int* result, /**< sol int array */
2781 SCIP_Real minpathcost, /**< the required reduced path cost to be surpassed */
2782 int root, /**< the root */
2783 int* nodearr, /**< for internal stuff */
2784 STP_Bool* marked /**< edge array to mark which (directed) edge can be removed */
2785 )
2786 {
2787 unsigned int* eqstack;
2788 SCIP_Bool* eqmark;
2789 int nfixed = 0;
2790 const int nedges = graph->edges;
2791 const int nnodes = graph->knots;
2792 const int halfnedges = graph->edges / 2;
2793 const SCIP_Bool rpc = (graph->stp_type == STP_RPCSPG);
2794
2795 if( rpc )
2796 graph_pc_2orgcheck(graph);
2797
2798 if( SCIPisZero(scip, minpathcost) )
2799 return 0;
2800
2801 SCIP_CALL( SCIPallocCleanBufferArray(scip, &eqmark, halfnedges) );
2802 SCIP_CALL( SCIPallocBufferArray(scip, &eqstack, halfnedges) );
2803
2804 /* main loop */
2805 for( int e = 0; e < nedges; e += 2 )
2806 {
2807 if( graph->oeat[e] != EAT_FREE )
2808 {
2809 const int erev = e + 1;
2810 int eqstack_size = 0;
2811 SCIP_Bool deletable = TRUE;
2812 const SCIP_Bool allowequality = (result != NULL && result[e] != CONNECT && result[erev] != CONNECT);
2813
2814 assert(graph->oeat[erev] != EAT_FREE);
2815
2816 if( SCIPisZero(scip, cost[e]) || SCIPisZero(scip, cost[erev]) )
2817 continue;
2818
2819 if( marked == NULL || !marked[e] )
2820 {
2821 SCIP_CALL_ABORT(reduceCheckEdge(scip, graph, root, cost, pathdist, vnoi, minpathcost, e, allowequality, nodearr,
2822 &deletable, eqstack, &eqstack_size, eqmark));
2823
2824 if( marked != NULL && deletable )
2825 marked[e] = TRUE;
2826 }
2827
2828 if( (marked == NULL || !marked[erev]) && deletable )
2829 {
2830 SCIP_CALL_ABORT(reduceCheckEdge(scip, graph, root, cost, pathdist, vnoi, minpathcost, erev, allowequality,
2831 nodearr, &deletable, eqstack, &eqstack_size, eqmark));
2832
2833 if( marked != NULL && deletable )
2834 marked[erev] = TRUE;
2835 }
2836
2837 for( int i = 0; i < eqstack_size; i++ )
2838 eqmark[eqstack[i]] = FALSE;
2839 for( int i = 0; i < halfnedges; i++ )
2840 assert(eqmark[i] == FALSE);
2841
2842 if( deletable )
2843 {
2844 graph_edge_del(scip, graph, e, TRUE);
2845 nfixed++;
2846 }
2847 }
2848 }
2849
2850 SCIPfreeBufferArray(scip, &eqstack);
2851 SCIPfreeCleanBufferArray(scip, &eqmark);
2852
2853 for( int k = 0; k < nnodes; k++ )
2854 if( graph->grad[k] == 0 && k != root )
2855 graph->mark[k] = FALSE;
2856
2857 return nfixed;
2858 }
2859
2860 /** dual ascent based reductions */
reduce_da(SCIP * scip,GRAPH * graph,PATH * vnoi,GNODE ** gnodearr,SCIP_Real * cost,SCIP_Real * costrev,SCIP_Real * pathdist,SCIP_Real * ub,SCIP_Real * offset,int * edgearrint,int * vbase,int * state,int * pathedge,int * nodearrint,int * heursources,STP_Bool * nodearrchar,int * nelims,int prevrounds,SCIP_RANDNUMGEN * randnumgen,SCIP_Bool userec,SCIP_Bool extended)2861 SCIP_RETCODE reduce_da(
2862 SCIP* scip, /**< SCIP data structure */
2863 GRAPH* graph, /**< graph data structure */
2864 PATH* vnoi, /**< Voronoi data structure */
2865 GNODE** gnodearr, /**< GNODE* terminals array for internal computations or NULL */
2866 SCIP_Real* cost, /**< edge costs */
2867 SCIP_Real* costrev, /**< reverse edge costs */
2868 SCIP_Real* pathdist, /**< distance array for shortest path calculations */
2869 SCIP_Real* ub, /**< pointer to provide upper bound and return upper bound found during ascent and prune (if better) */
2870 SCIP_Real* offset, /**< pointer to store offset */
2871 int* edgearrint, /**< int edges array for internal computations or NULL */
2872 int* vbase, /**< array for Voronoi bases */
2873 int* state, /**< int 4 * nnodes array for internal computations */
2874 int* pathedge, /**< array for predecessor edge on a path */
2875 int* nodearrint, /**< int nnodes array for internal computations */
2876 int* heursources, /**< array to store starting points of TM heuristic */
2877 STP_Bool* nodearrchar, /**< STP_Bool node array for internal computations */
2878 int* nelims, /**< pointer to store number of reduced edges */
2879 int prevrounds, /**< number of reduction rounds that have been performed already */
2880 SCIP_RANDNUMGEN* randnumgen, /**< random number generator */
2881 SCIP_Bool userec, /**< use recombination heuristic? */
2882 SCIP_Bool extended /**< use extended tests? */
2883 )
2884 {
2885 STPSOLPOOL* pool;
2886
2887 SCIP_Real* edgefixingbounds;
2888 SCIP_Real* nodefixingbounds;
2889 SCIP_Real* nodereplacebounds;
2890 SCIP_Real lpobjval;
2891 SCIP_Real upperbound;
2892 SCIP_Real minpathcost;
2893 SCIP_Real damaxdeviation;
2894 SCIP_Bool success;
2895 const SCIP_Bool rpc = (graph->stp_type == STP_RPCSPG);
2896 const SCIP_Bool directed = (graph->stp_type == STP_SAP || graph->stp_type == STP_NWSPG);
2897 int* terms;
2898 int* result;
2899 int* starts;
2900 int k;
2901 int runs;
2902 int nruns;
2903 int nterms;
2904 const int nedges = graph->edges;
2905 const int nnodes = graph->knots;
2906 int nfixed;
2907 int best_start;
2908 STP_Bool* marked;
2909
2910 assert(ub != NULL);
2911 assert(scip != NULL);
2912 assert(graph != NULL);
2913 assert(nelims != NULL);
2914 assert(nodearrint != NULL);
2915
2916 nfixed = 0;
2917 minpathcost = 0.0;
2918
2919 if( graph->terms <= 1 )
2920 return SCIP_OKAY;
2921
2922 if( SCIPisGE(scip, *ub, 0.0) )
2923 upperbound = *ub;
2924 else
2925 upperbound = FARAWAY;
2926
2927 /* allocate memory */
2928 SCIP_CALL( SCIPallocBufferArray(scip, &result, nedges) );
2929 SCIP_CALL( SCIPallocBufferArray(scip, &marked, nedges) );
2930 SCIP_CALL( SCIPallocBufferArray(scip, &edgefixingbounds, nedges) );
2931 SCIP_CALL( SCIPallocBufferArray(scip, &nodefixingbounds, nnodes) );
2932 SCIP_CALL( SCIPallocBufferArray(scip, &nodereplacebounds, nnodes) );
2933
2934 if( !rpc && !directed )
2935 SCIP_CALL( SCIPallocBufferArray(scip, &terms, graph->terms) );
2936 else
2937 terms = NULL;
2938
2939 /*
2940 * 1. step: compute upper bound
2941 */
2942 if( userec )
2943 SCIP_CALL( SCIPStpHeurRecInitPool(scip, &pool, nedges, SOLPOOL_SIZE) );
2944 else
2945 pool = NULL;
2946
2947 if( !rpc && extended )
2948 SCIP_CALL( reduce_deleteConflictEdges(scip, graph) );
2949
2950 /* initialize */
2951 k = 0;
2952 nterms = 0;
2953 for( int i = 0; i < nnodes; i++ )
2954 {
2955 if( !rpc )
2956 graph->mark[i] = (graph->grad[i] > 0);
2957
2958 if( graph->mark[i] )
2959 {
2960 k++;
2961 if( Is_term(graph->term[i]) )
2962 {
2963 if( !rpc && !directed )
2964 terms[nterms] = i;
2965 nterms++;
2966 }
2967 }
2968 }
2969
2970 /* not more than two terminals? */
2971 if( nterms <= 2 )
2972 goto TERMINATE;
2973
2974 /* number of runs should not exceed number of connected vertices */
2975 if( directed )
2976 runs = MIN(k, DEFAULT_HEURRUNS);
2977 else
2978 runs = MIN(k, DEFAULT_HEURRUNS / 5);
2979
2980 /* neither PC, MW, RPC nor HC? */
2981 if( !rpc && heursources != NULL )
2982 {
2983 /* choose starting points for TM heuristic */
2984 starts = heursources;
2985 SCIPStpHeurTMCompStarts(graph, starts, &runs);
2986 }
2987 else
2988 {
2989 starts = NULL;
2990 }
2991
2992 for( int e = 0; e < nedges; e++ )
2993 {
2994 cost[e] = graph->cost[e];
2995 costrev[e] = graph->cost[flipedge(e)];
2996 result[e] = UNKNOWN;
2997 }
2998
2999 if( directed || rpc )
3000 {
3001 SCIP_Real ubnew;
3002
3003 SCIP_CALL( SCIPStpHeurTMRun(scip, NULL, graph, starts, &best_start, result, runs, graph->source, cost, costrev, NULL, NULL, 0.0, &success, FALSE) );
3004 assert(success);
3005
3006 /* calculate objective value of solution */
3007 ubnew = graph_sol_getObj(graph->cost, result, 0.0, nedges);
3008
3009 if( SCIPisLT(scip, ubnew, upperbound) )
3010 upperbound = ubnew;
3011 }
3012
3013 /*
3014 * 2. step: repeatedly compute lower bound and reduced costs and try to eliminate
3015 */
3016
3017
3018 damaxdeviation = -1.0;
3019
3020 /* if not RPC, select roots for dual ascent */
3021 if( !rpc && !directed )
3022 {
3023 assert(terms != NULL);
3024
3025 nruns = MIN(graph->terms, DEFAULT_DARUNS);
3026 SCIP_CALL( orderDaRoots(scip, graph, terms, graph->terms, (prevrounds > 0), randnumgen) );
3027
3028 if( prevrounds > 0 )
3029 damaxdeviation = SCIPrandomGetReal(randnumgen, DAMAXDEVIATION_RANDOM_LOWER, DAMAXDEVIATION_RANDOM_UPPER);
3030 }
3031 else
3032 {
3033 nruns = 1;
3034 if( rpc )
3035 graph_pc_2transcheck(graph);
3036 }
3037
3038 assert(nruns > 0);
3039
3040 for( int outerrounds = 0; outerrounds < 2; outerrounds++ )
3041 {
3042 /* main reduction loop */
3043 for( int run = 0; run < nruns; run++ )
3044 {
3045 int root = graph->source;
3046 SCIP_Bool apsol = FALSE;
3047 SCIP_Bool usesol = (run > 1);
3048
3049 if( rpc )
3050 usesol = TRUE;
3051
3052 /* graph vanished? */
3053 if( graph->grad[graph->source] == 0 )
3054 break;
3055
3056 if( !rpc && !directed )
3057 {
3058 assert(terms != NULL);
3059 root = terms[run];
3060 usesol = usesol && (SCIPrandomGetInt(randnumgen, 0, 2) < 2) && graph->stp_type != STP_RSMT;
3061 }
3062
3063 if( usesol )
3064 {
3065 /* solution might not be valid anymore */
3066 if( !graph_sol_valid(scip, graph, result) )
3067 {
3068 SCIPdebugMessage("solution not valid; run normal dual-ascent \n");
3069 SCIP_CALL(
3070 SCIPStpDualAscent(scip, graph, cost, pathdist, &lpobjval, FALSE, FALSE, gnodearr, NULL, edgearrint, state, root, FALSE, damaxdeviation, nodearrchar));
3071 }
3072 else
3073 {
3074 if( !rpc )
3075 {
3076 SCIPdebugMessage("reroot solution and run guided dual-ascent \n");
3077 SCIP_CALL(graph_sol_reroot(scip, graph, result, root));
3078 #ifndef NDEBUG
3079 {
3080 const int realroot = graph->source;
3081 graph->source = root;
3082 assert(graph_sol_valid(scip, graph, result));
3083 graph->source = realroot;
3084 }
3085 #endif
3086 }
3087
3088 SCIP_CALL(
3089 SCIPStpDualAscent(scip, graph, cost, pathdist, &lpobjval, FALSE, FALSE, gnodearr, result, edgearrint, state, root, FALSE, damaxdeviation, nodearrchar));
3090 }
3091 }
3092 else
3093 {
3094 SCIPdebugMessage("no rerooting, run normal dual-ascent \n");
3095 SCIP_CALL(
3096 SCIPStpDualAscent(scip, graph, cost, pathdist, &lpobjval, FALSE, FALSE, gnodearr, NULL, edgearrint, state, root, FALSE, damaxdeviation, nodearrchar));
3097 }
3098
3099 /* perform ascent and prune */
3100 if( !directed )
3101 {
3102 SCIP_Real ubnew;
3103 SCIP_Bool soladded = FALSE;
3104
3105 SCIP_CALL(SCIPStpHeurAscendPruneRun(scip, NULL, graph, cost, result, nodearrint, root, nodearrchar, &success, FALSE));
3106 assert(success && graph_sol_valid(scip, graph, result));
3107
3108 ubnew = graph_sol_getObj(graph->cost, result, 0.0, nedges);
3109
3110 if( userec && !rpc )
3111 {
3112 SCIPdebugMessage("obj before local %f \n", ubnew);
3113
3114 SCIP_CALL(SCIPStpHeurLocalRun(scip, graph, graph->cost, result));
3115 ubnew = graph_sol_getObj(graph->cost, result, 0.0, nedges);
3116
3117 SCIPdebugMessage("obj after local %f \n", ubnew);
3118 assert(graph_sol_valid(scip, graph, result));
3119
3120 SCIP_CALL(SCIPStpHeurRecAddToPool(scip, ubnew, result, pool, &soladded));
3121 }
3122
3123 if( userec && soladded && pool->size >= 2 && ubnew < upperbound )
3124 {
3125 /* get index of just added solution */
3126 int solindex = pool->maxindex;
3127 SCIP_Bool solfound;
3128
3129 SCIPdebugMessage("POOLSIZE %d \n", pool->size);
3130
3131 SCIP_CALL(SCIPStpHeurRecRun(scip, pool, NULL, NULL, graph, NULL, &solindex, 1, pool->size, FALSE, &solfound));
3132
3133 if( solfound )
3134 {
3135 const STPSOL* const sol = SCIPStpHeurRecSolfromIdx(pool, solindex);
3136
3137 assert(sol != NULL);
3138
3139 SCIPdebugMessage("DA: rec found better solution with obj %f vs %f \n", sol->obj, ubnew);
3140
3141 if( sol->obj < ubnew )
3142 {
3143 BMScopyMemoryArray(result, sol->soledges, nedges);
3144
3145 SCIP_CALL(SCIPStpHeurLocalRun(scip, graph, graph->cost, result));
3146 ubnew = graph_sol_getObj(graph->cost, result, 0.0, nedges);
3147
3148 assert(SCIPisLE(scip, ubnew, sol->obj));
3149
3150 if( ubnew < sol->obj )
3151 SCIP_CALL(SCIPStpHeurRecAddToPool(scip, ubnew, result, pool, &solfound));
3152 }
3153 }
3154 }
3155
3156 if( SCIPisLE(scip, ubnew, upperbound) )
3157 {
3158 apsol = TRUE;
3159 upperbound = ubnew;
3160 }
3161 }
3162
3163 /* the required reduced path cost to be surpassed */
3164 minpathcost = upperbound - lpobjval;
3165 assert(SCIPisGE(scip, minpathcost, 0.0));
3166
3167 SCIPdebugMessage("upper: %f lower: %f \n", upperbound, lpobjval);
3168
3169 for( k = 0; k < nnodes; k++ )
3170 graph->mark[k] = (graph->grad[k] > 0);
3171
3172 /* distance from root to all nodes */
3173 graph_path_execX(scip, graph, root, cost, pathdist, pathedge);
3174
3175 for( unsigned e = 0; e < (unsigned) nedges; e++ )
3176 {
3177 marked[e] = FALSE;
3178 costrev[e] = cost[flipedge(e)];
3179 }
3180
3181 /* no paths should go back to the root */
3182 for( int e = graph->outbeg[root]; e != EAT_LAST; e = graph->oeat[e] )
3183 costrev[e] = FARAWAY;
3184
3185 /* build Voronoi diagram */
3186 if( directed || rpc )
3187 graph_voronoiTerms(scip, graph, costrev, vnoi, vbase, graph->path_heap, state);
3188 else
3189 {
3190 graph_get4nextTerms(scip, graph, costrev, costrev, vnoi, vbase, graph->path_heap, state);
3191
3192 #ifndef NDEBUG
3193 for( int i = 0; i < nnodes; i++ )
3194 if( !Is_term(graph->term[i]) )
3195 {
3196 assert(vbase[i] != root || vnoi[i].dist >= FARAWAY);
3197 assert(vbase[i + nnodes] != root || vnoi[i + nnodes].dist >= FARAWAY);
3198 }
3199 else
3200 assert(vbase[i] == i);
3201 #endif
3202 updateNodeFixingBounds(nodefixingbounds, graph, pathdist, vnoi, lpobjval, (run == 0));
3203 updateEdgeFixingBounds(edgefixingbounds, graph, cost, pathdist, vnoi, lpobjval, nedges, (run == 0), TRUE);
3204 }
3205
3206 nfixed += reduceSPG(scip, graph, offset, marked, nodearrchar, vnoi, cost, pathdist, result, minpathcost, root, apsol);
3207
3208 if( !directed && !rpc )
3209 {
3210 if( !SCIPisZero(scip, minpathcost) )
3211 {
3212 nfixed += reduceWithNodeFixingBounds(scip, graph, NULL, nodefixingbounds, upperbound);
3213 nfixed += reduceWithEdgeFixingBounds(scip, graph, NULL, edgefixingbounds, (apsol ? result : NULL), upperbound);
3214 }
3215
3216 if( extended )
3217 {
3218 nfixed += reduce_extendedEdge(scip, graph, vnoi, cost, pathdist, (apsol ? result : NULL), minpathcost, root, nodearrint, marked);
3219 }
3220
3221 if( !SCIPisZero(scip, minpathcost) )
3222 SCIP_CALL(
3223 updateNodeReplaceBounds(scip, nodereplacebounds, graph, cost, pathdist, vnoi, vbase, nodearrint, lpobjval, upperbound, root,
3224 (run == 0), extended));
3225 }
3226
3227 if( nfixed > 0 )
3228 SCIP_CALL(level0(scip, graph));
3229 assert(graph_valid(graph));
3230
3231 if( !rpc )
3232 {
3233 for( k = 0; k < nnodes; k++ )
3234 graph->mark[k] = graph->grad[k] > 0;
3235 }
3236 else
3237 {
3238 graph_pc_2trans(graph);
3239 graph_pc_2org(graph);
3240 }
3241 } /* root loop */
3242
3243 if( !directed && !rpc && !SCIPisZero(scip, minpathcost) )
3244 nfixed += reduceWithNodeReplaceBounds(scip, graph, vnoi, pathdist, cost, nodereplacebounds, nodearrint, lpobjval, upperbound);
3245
3246 if( nfixed == 0 || !userec )
3247 {
3248 break;
3249 }
3250 else if( userec )
3251 {
3252 damaxdeviation = SCIPrandomGetReal(randnumgen, DAMAXDEVIATION_RANDOM_LOWER, DAMAXDEVIATION_RANDOM_UPPER);
3253 }
3254 }
3255
3256 TERMINATE:
3257 *nelims = nfixed;
3258
3259 if( SCIPisLT(scip, upperbound, *ub) || SCIPisLT(scip, *ub, 0.0) )
3260 *ub = upperbound;
3261
3262 if( userec )
3263 SCIPStpHeurRecFreePool(scip, &pool);
3264
3265 SCIPfreeBufferArrayNull(scip, &terms);
3266 SCIPfreeBufferArray(scip, &nodereplacebounds);
3267 SCIPfreeBufferArray(scip, &nodefixingbounds);
3268 SCIPfreeBufferArray(scip, &edgefixingbounds);
3269 SCIPfreeBufferArray(scip, &marked);
3270 SCIPfreeBufferArray(scip, &result);
3271
3272 assert(graph_valid(graph));
3273
3274 return SCIP_OKAY;
3275 }
3276
3277
3278 /** dual ascent reduction for slack-and-prune heuristic */
reduce_daSlackPrune(SCIP * scip,SCIP_VAR ** vars,GRAPH * graph,PATH * vnoi,GNODE ** gnodearr,SCIP_Real * cost,SCIP_Real * costrev,SCIP_Real * pathdist,SCIP_Real * upperbound,int * edgearrint,int * edgearrint2,int * vbase,int * pathedge,int * state,int * solnode,STP_Bool * nodearrchar,STP_Bool * edgearrchar,int * nelims,int minelims,SCIP_Bool solgiven)3279 SCIP_RETCODE reduce_daSlackPrune(
3280 SCIP* scip, /**< SCIP data structure */
3281 SCIP_VAR** vars, /**< problem variables or NULL */
3282 GRAPH* graph, /**< graph data structure */
3283 PATH* vnoi, /**< Voronoi data structure */
3284 GNODE** gnodearr, /**< GNODE* terminals array for internal computations or NULL */
3285 SCIP_Real* cost, /**< array to store reduced costs */
3286 SCIP_Real* costrev, /**< reverse edge costs */
3287 SCIP_Real* pathdist, /**< distance array for shortest path calculations */
3288 SCIP_Real* upperbound, /**< pointer to store new upper bound */
3289 int* edgearrint, /**< int edges array to store solution value */
3290 int* edgearrint2, /**< int edges array for internal computations */
3291 int* vbase, /**< array for Voronoi bases */
3292 int* pathedge, /**< array for predecessor edge on a path */
3293 int* state, /**< int 4 * nnodes array for internal computations */
3294 int* solnode, /**< array of nodes of current solution that is not to be destroyed */
3295 STP_Bool* nodearrchar, /**< STP_Bool node array for internal computations */
3296 STP_Bool* edgearrchar, /**< STP_Bool edge array for internal computations */
3297 int* nelims, /**< pointer to store number of reduced edges */
3298 int minelims, /**< minimum number of edges to eliminate */
3299 SCIP_Bool solgiven /**< solution provided? */
3300 )
3301 {
3302 IDX** ancestors;
3303 IDX** revancestors;
3304 SCIP_Real obj;
3305 SCIP_Real tmpcost;
3306 SCIP_Real lpobjval;
3307 SCIP_Real objprune;
3308 SCIP_Real minpathcost;
3309 SCIP_Real* sd;
3310 SCIP_Real* ecost;
3311 SCIP_Bool rpc;
3312 SCIP_Bool success;
3313 SCIP_Bool eliminate;
3314
3315 int* grad;
3316 int* adjvert;
3317 int* incedge;
3318 int* reinsert;
3319 int i;
3320 int k;
3321 int e;
3322 int e2;
3323 int e3;
3324 int etmp;
3325 int root;
3326 int head;
3327 int nterms;
3328 int nedges;
3329 int nnodes;
3330 int nfixed;
3331 int redrounds;
3332 STP_Bool* marked;
3333
3334 assert(scip != NULL);
3335 assert(cost != NULL);
3336 assert(graph != NULL);
3337 assert(nelims != NULL);
3338 assert(solnode != NULL);
3339 assert(costrev != NULL);
3340 assert(pathedge != NULL);
3341 assert(upperbound != NULL);
3342 assert(edgearrint != NULL);
3343 assert(edgearrint2 != NULL);
3344 assert(nodearrchar != NULL);
3345 assert(edgearrchar != NULL);
3346
3347 /* 1. step: initialize */
3348
3349 rpc = (graph->stp_type == STP_RPCSPG);
3350 grad = graph->grad;
3351 root = graph->source;
3352 nfixed = 0;
3353 nedges = graph->edges;
3354 nnodes = graph->knots;
3355
3356 /* graph vanished? */
3357 if( grad[graph->source] == 0 )
3358 return SCIP_OKAY;
3359
3360 marked = edgearrchar;
3361
3362 /* allocate length-4 buffer memory */
3363 SCIP_CALL( SCIPallocBufferArray(scip, &sd, 4) );
3364 SCIP_CALL( SCIPallocBufferArray(scip, &ancestors, 4) );
3365 SCIP_CALL( SCIPallocBufferArray(scip, &revancestors, 4) );
3366 SCIP_CALL( SCIPallocBufferArray(scip, &ecost, 4) );
3367 SCIP_CALL( SCIPallocBufferArray(scip, &adjvert, 4) );
3368 SCIP_CALL( SCIPallocBufferArray(scip, &reinsert, 4) );
3369 SCIP_CALL( SCIPallocBufferArray(scip, &incedge, 4) );
3370
3371 for( i = 0; i < 4; i++ )
3372 {
3373 sd[i] = 0.0;
3374 ancestors[i] = NULL;
3375 revancestors[i] = NULL;
3376 }
3377
3378 k = 0;
3379 nterms = 0;
3380 for( i = 0; i < nnodes; i++ )
3381 {
3382 if( !rpc )
3383 graph->mark[i] = (grad[i] > 0);
3384 if( graph->mark[i] )
3385 {
3386 k++;
3387 if( Is_term(graph->term[i]) )
3388 nterms++;
3389 }
3390 }
3391
3392 /* not more than two terminals? */
3393 if( nterms <= 2 )
3394 goto TERMINATE;
3395
3396 /* 2. step: - if not provided, compute lower bound and reduced costs
3397 * - try to eliminate edges and nodes */
3398
3399 for( i = 0; i < nnodes; i++ )
3400 if( Is_term(graph->term[i]) )
3401 assert(grad[i] > 0);
3402
3403 if( rpc )
3404 {
3405 graph_pc_2trans(graph);
3406 }
3407
3408
3409 if( !solgiven )
3410 {
3411 /* try to build MST on solnode nodes */
3412 for( i = 0; i < nnodes; i++ )
3413 graph->mark[i] = (solnode[i] == CONNECT);
3414
3415 for( e = 0; e < nedges; e++ )
3416 edgearrint[e] = UNKNOWN;
3417
3418 graph_path_exec(scip, graph, MST_MODE, root, graph->cost, vnoi);
3419
3420 for( i = 0; i < nnodes; i++ )
3421 {
3422 e = vnoi[i].edge;
3423 if( e >= 0 )
3424 {
3425 edgearrint[e] = CONNECT;
3426 }
3427 else if( Is_term(graph->term[i]) && i != root )
3428 {
3429 break;
3430 }
3431 }
3432
3433 if( i == nnodes )
3434 {
3435 int l;
3436 int count;
3437
3438 do
3439 {
3440 count = 0;
3441
3442 for( l = nnodes - 1; l >= 0; --l )
3443 {
3444 if( (solnode[l] != CONNECT) || Is_term(graph->term[l]) )
3445 continue;
3446
3447 for( e = graph->outbeg[l]; e != EAT_LAST; e = graph->oeat[e] )
3448 if( edgearrint[e] == CONNECT )
3449 break;
3450
3451 if( e == EAT_LAST )
3452 {
3453 /* there has to be exactly one incoming edge */
3454 for( e = graph->inpbeg[l]; e != EAT_LAST; e = graph->ieat[e] )
3455 {
3456 if( edgearrint[e] == CONNECT )
3457 {
3458 edgearrint[e] = UNKNOWN;
3459 solnode[l] = UNKNOWN;
3460 count++;
3461 break;
3462 }
3463 }
3464 }
3465 }
3466 }
3467 while( count > 0 );
3468 }
3469 }
3470
3471
3472 if( solgiven || i == nnodes )
3473 {
3474 obj = graph_sol_getObj(graph->cost, edgearrint, 0.0, nedges);
3475
3476 SCIP_CALL( SCIPStpDualAscent(scip, graph, cost, pathdist, &lpobjval, FALSE, FALSE, gnodearr, edgearrint, edgearrint2, state, root, FALSE, -1.0, nodearrchar) );
3477 }
3478 else
3479 {
3480 obj = FARAWAY;
3481 SCIP_CALL( SCIPStpDualAscent(scip, graph, cost, pathdist, &lpobjval, FALSE, FALSE, gnodearr, NULL, edgearrint2, state, root, FALSE, -1.0, nodearrchar) );
3482 }
3483
3484 #if 0
3485 SCIP_QUEUE* queue;
3486 SCIP_CALL( SCIPqueueCreate(&queue, nnodes, 2.0) );
3487
3488 GRAPH* prunegraph = graph;
3489 int* mark = graph->mark;
3490 int* pnode;
3491 int a;
3492 for( k = 0; k < nnodes; k++ )
3493 mark[k] = FALSE;
3494
3495
3496 /* BFS from root along incoming arcs of zero cost */
3497
3498 mark[prunegraph->source] = TRUE;
3499
3500 SCIP_CALL( SCIPqueueInsert(queue, &(prunegraph->source)) );
3501
3502 while( !SCIPqueueIsEmpty(queue) )
3503 {
3504 pnode = (SCIPqueueRemove(queue));
3505 k = *pnode;
3506
3507 /* traverse outgoing arcs */
3508 for( a = prunegraph->outbeg[k]; a != EAT_LAST; a = prunegraph->oeat[a] )
3509 {
3510 head = prunegraph->head[a];
3511
3512 if( SCIPisEQ(scip, cost[a], 0.0) )
3513 {
3514 /* vertex not labeled yet? */
3515 if( !mark[head] )
3516 {
3517 mark[head] = TRUE;
3518 SCIP_CALL( SCIPqueueInsert(queue, &(prunegraph->head[a])) );
3519 }
3520 }
3521 }
3522 }
3523 SCIPqueueFree(&queue);
3524 for( k = 0; k < nnodes; k++ ) // TODO
3525 if( Is_term(prunegraph->term[k]) && !mark[k] )
3526 printf("in bnd FAIL %d not marked, but terminal, \n", k);
3527 #endif
3528
3529 SCIP_CALL( SCIPStpHeurAscendPruneRun(scip, NULL, graph, cost, edgearrint2, vbase, root, nodearrchar, &success, FALSE) );
3530
3531 objprune = graph_sol_getObj(graph->cost, edgearrint2, 0.0, nedges);
3532
3533 assert(success);
3534
3535 if( success && SCIPisLT(scip, objprune, obj ) )
3536 {
3537
3538 for( i = 0; i < nnodes; i++ )
3539 solnode[i] = UNKNOWN;
3540
3541 for( e = 0; e < nedges; e++ )
3542 {
3543 edgearrint[e] = edgearrint2[e];
3544 if( edgearrint[e] == CONNECT )
3545 {
3546 solnode[graph->tail[e]] = CONNECT;
3547 solnode[graph->head[e]] = CONNECT;
3548 }
3549 }
3550 }
3551
3552 obj = 0.0;
3553
3554 for( e = 0; e < nedges; e++ )
3555 {
3556 if( edgearrint[e] == CONNECT )
3557 obj += graph->cost[e];
3558
3559 marked[e] = FALSE;
3560 costrev[e] = cost[flipedge(e)];
3561 }
3562
3563 *upperbound = obj;
3564
3565 for( k = 0; k < nnodes; k++ )
3566 graph->mark[k] = (grad[k] > 0);
3567
3568 /* distance from root to all nodes */
3569 graph_path_execX(scip, graph, root, cost, pathdist, pathedge);
3570
3571 /* no paths should go back to the root */
3572 for( e = graph->outbeg[root]; e != EAT_LAST; e = graph->oeat[e] )
3573 costrev[e] = FARAWAY;
3574
3575 /* build Voronoi diagram */
3576 graph_get4nextTerms(scip, graph, costrev, costrev, vnoi, vbase, graph->path_heap, state);
3577
3578 for( k = 0; k < nnodes; k++ )
3579 if( !Is_term(graph->term[k]) )
3580 assert(vbase[k + nnodes] != root );
3581
3582 /* RPC? If yes, restore original graph */
3583 if( rpc )
3584 {
3585 graph_pc_2org(graph);
3586 graph->mark[root] = FALSE;
3587 }
3588
3589 for( e = 0; e < nedges; e++ )
3590 costrev[e] = -1.0;
3591
3592 for( redrounds = 0; redrounds < 3; redrounds++ )
3593 {
3594 if( redrounds == 0 )
3595 {
3596 eliminate = FALSE;
3597 minpathcost = FARAWAY;
3598 }
3599 else if( redrounds == 1 )
3600 {
3601 assert(minelims > 0);
3602 assert(2 * minelims < nedges);
3603
3604 eliminate = TRUE;
3605 SCIPsortReal(costrev, nedges);
3606
3607 /* the required reduced path cost to be surpassed */
3608 minpathcost = costrev[nedges - 2 * minelims];
3609
3610 if( SCIPisLE(scip, minpathcost, 0.0) )
3611 minpathcost = 2 * SCIPepsilon(scip);
3612
3613 k = nedges - 2 * minelims;
3614
3615 /* try to first eliminate edges with higher gap */
3616 for( e = nedges - 1; e > k && e >= 2; e = e - 2 )
3617 {
3618 if( SCIPisLE(scip, costrev[e - 2], minpathcost) )
3619 break;
3620 }
3621
3622 if( SCIPisGT(scip, costrev[e], minpathcost) )
3623 minpathcost = costrev[e];
3624
3625 if( SCIPisLE(scip, minpathcost, 0.0) )
3626 minpathcost = 2 * SCIPepsilon(scip);
3627
3628 for( e = 0; e < nedges; e++ )
3629 marked[e] = FALSE;
3630 }
3631 else
3632 {
3633 eliminate = TRUE;
3634
3635 /* the required reduced path cost to be surpassed */
3636 minpathcost = costrev[nedges - 2 * minelims];
3637
3638 if( SCIPisLE(scip, minpathcost, 0.0) )
3639 minpathcost = 2 * SCIPepsilon(scip);
3640
3641 for( e = 0; e < nedges; e++ )
3642 marked[e] = FALSE;
3643 }
3644
3645 for( k = 0; k < nnodes; k++ )
3646 {
3647 if( grad[k] <= 0 )
3648 continue;
3649
3650 if( nfixed > minelims )
3651 break;
3652
3653 if( !Is_term(graph->term[k]) && (!eliminate || pathdist[k] + vnoi[k].dist >= minpathcost) && solnode[k] != CONNECT )
3654 {
3655 if( !eliminate )
3656 {
3657 tmpcost = pathdist[k] + vnoi[k].dist;
3658
3659 for( e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
3660 {
3661 if( SCIPisGT(scip, tmpcost, costrev[e]) )
3662 costrev[e] = tmpcost;
3663
3664 e2 = flipedge(e);
3665
3666 if( SCIPisGT(scip, tmpcost, costrev[e2]) )
3667 costrev[e2] = tmpcost;
3668 }
3669
3670 continue;
3671 }
3672 nfixed += grad[k];
3673
3674 while( graph->outbeg[k] != EAT_LAST )
3675 graph_edge_del(scip, graph, graph->outbeg[k], TRUE);
3676 }
3677 else
3678 {
3679 e = graph->outbeg[k];
3680 while( e != EAT_LAST )
3681 {
3682 etmp = graph->oeat[e];
3683 head = graph->head[e];
3684
3685 /* for rpc no artificial terminal arcs should be deleted; in general: delete no solution edges */
3686 if( (rpc && !graph->mark[head])
3687 || (edgearrint[e] == CONNECT) || (edgearrint[flipedge(e)] == CONNECT) )
3688 {
3689 e = etmp;
3690 continue;
3691 }
3692
3693 tmpcost = pathdist[k] + cost[e] + vnoi[head].dist;
3694
3695 if( (!eliminate) || tmpcost >= minpathcost )
3696 {
3697 if( marked[flipedge(e)] )
3698 {
3699 if( eliminate )
3700 {
3701 graph_edge_del(scip, graph, e, TRUE);
3702 nfixed++;
3703 }
3704 else
3705 {
3706 if( SCIPisGT(scip, tmpcost, costrev[e]) )
3707 costrev[e] = tmpcost;
3708
3709 if( SCIPisGT(scip, tmpcost, costrev[flipedge(e)]) )
3710 costrev[flipedge(e)] = tmpcost;
3711 }
3712 }
3713 else
3714 {
3715 marked[e] = TRUE;
3716 }
3717 }
3718 e = etmp;
3719 }
3720 }
3721 }
3722
3723 for( k = 0; k < nnodes; k++ )
3724 {
3725 if( nfixed > minelims )
3726 break;
3727
3728 if( !graph->mark[k] || Is_term(graph->term[k]) || solnode[k] == CONNECT )
3729 continue;
3730
3731 if( grad[k] == 3 )
3732 {
3733 tmpcost = pathdist[k] + vnoi[k].dist + vnoi[k + nnodes].dist;
3734
3735 if( !eliminate || tmpcost >= minpathcost )
3736 {
3737 e = graph->outbeg[k];
3738 assert(graph->oeat[e] != EAT_LAST);
3739 e2 = graph->oeat[e];
3740 assert(graph->oeat[e2] != EAT_LAST);
3741 e3 = graph->oeat[e2];
3742 assert(graph->oeat[e3] == EAT_LAST);
3743
3744 if( SCIPisLE(scip, cost[e], 0.0) || SCIPisLE(scip, cost[e2], 0.0) || SCIPisLE(scip, cost[e3], 0.0) )
3745 continue;
3746
3747 if( !eliminate )
3748 {
3749 for( e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
3750 {
3751 if( SCIPisGT(scip, tmpcost, costrev[e]) )
3752 costrev[e] = tmpcost;
3753
3754 e2 = flipedge(e);
3755
3756 if( SCIPisGT(scip, tmpcost, costrev[e2]) )
3757 costrev[e2] = tmpcost;
3758 }
3759
3760 continue;
3761 }
3762 nfixed += 3;
3763 while( graph->outbeg[k] != EAT_LAST )
3764 graph_edge_del(scip, graph, graph->outbeg[k], TRUE);
3765 }
3766 }
3767 }
3768 }
3769 SCIPdebugMessage("deleted by da: %d \n", nfixed );
3770
3771 if( rpc )
3772 graph_pc_2trans(graph);
3773 assert(graph->mark[root]);
3774
3775 TERMINATE:
3776 *nelims = nfixed;
3777
3778 /* free memory */
3779 for( k = 0; k < 4; k++ )
3780 {
3781 SCIPintListNodeFree(scip, &(ancestors[k]));
3782 SCIPintListNodeFree(scip, &(revancestors[k]));
3783 }
3784
3785 SCIPfreeBufferArray(scip, &incedge);
3786 SCIPfreeBufferArray(scip, &reinsert);
3787 SCIPfreeBufferArray(scip, &adjvert);
3788 SCIPfreeBufferArray(scip, &ecost);
3789 SCIPfreeBufferArray(scip, &revancestors);
3790 SCIPfreeBufferArray(scip, &ancestors);
3791 SCIPfreeBufferArray(scip, &sd);
3792
3793 if( edgearrchar == NULL )
3794 SCIPfreeBufferArray(scip, &marked);
3795
3796 return SCIP_OKAY;
3797 }
3798
3799
3800 /** dual ascent based reductions for PCSPG and MWCSP */
reduce_daPcMw(SCIP * scip,GRAPH * graph,PATH * vnoi,GNODE ** gnodearr,SCIP_Real * cost,SCIP_Real * costrev,SCIP_Real * pathdist,int * vbase,int * pathedge,int * edgearrint,int * state,STP_Bool * nodearrchar,int * nelims,SCIP_Bool solbasedda,SCIP_Bool varyroot,SCIP_Bool markroots,SCIP_Bool userec,SCIP_Bool fastmode,SCIP_RANDNUMGEN * randnumgen,SCIP_Real prizesum)3801 SCIP_RETCODE reduce_daPcMw(
3802 SCIP* scip, /**< SCIP data structure */
3803 GRAPH* graph, /**< graph data structure */
3804 PATH* vnoi, /**< Voronoi data structure array */
3805 GNODE** gnodearr, /**< GNODE* terminals array for internal computations or NULL */
3806 SCIP_Real* cost, /**< reduced edge costs */
3807 SCIP_Real* costrev, /**< reduced reverse edge costs */
3808 SCIP_Real* pathdist, /**< distance array for shortest path calculations */
3809 int* vbase, /**< Voronoi base array */
3810 int* pathedge, /**< shortest path incoming edge array for shortest path calculations */
3811 int* edgearrint, /**< int edges array for internal computations or NULL */
3812 int* state, /**< int 4 * vertices array */
3813 STP_Bool* nodearrchar, /**< STP_Bool node array for internal computations */
3814 int* nelims, /**< pointer to store number of reduced edges */
3815 SCIP_Bool solbasedda, /**< rerun Da based on best primal solution */
3816 SCIP_Bool varyroot, /**< vary root for DA if possible */
3817 SCIP_Bool markroots, /**< should terminals proven to be part of an opt. sol. be marked as such? */
3818 SCIP_Bool userec, /**< use recombination heuristic? */
3819 SCIP_Bool fastmode, /**< run heuristic in fast mode? */
3820 SCIP_RANDNUMGEN* randnumgen, /**< random number generator */
3821 SCIP_Real prizesum /**< sum of positive prizes */
3822 )
3823 {
3824 STPSOLPOOL* pool;
3825 GRAPH* transgraph;
3826 SCIP_Real* bestcost;
3827 SCIP_Real* edgefixingbounds;
3828 SCIP_Real* nodefixingbounds;
3829 SCIP_Real ub;
3830 SCIP_Real offset;
3831 SCIP_Real lpobjval;
3832 SCIP_Real bestlpobjval;
3833 SCIP_Real upperbound;
3834 SCIP_Real minpathcost;
3835 SCIP_Real damaxdeviation;
3836 int* roots;
3837 int* result;
3838 int* result2;
3839 int* transresult;
3840 STP_Bool* marked;
3841 int nroots;
3842 int nfixed;
3843 int nusedroots;
3844 const int root = graph->source;
3845 const int nedges = graph->edges;
3846 const int extnedges = nedges + 2 * (graph->terms - 1);
3847 const int nnodes = graph->knots;
3848 SCIP_Bool apsol;
3849 SCIP_Bool success;
3850
3851 assert(scip != NULL);
3852 assert(graph != NULL);
3853 assert(nelims != NULL);
3854 assert(nodearrchar != NULL);
3855
3856 /* not more than one terminal? */
3857 if( graph->terms <= 1 )
3858 return SCIP_OKAY;
3859
3860 nroots = 0;
3861 nfixed = 0;
3862 varyroot = varyroot && userec;
3863
3864 /* allocate memory */
3865 if( edgearrint == NULL )
3866 SCIP_CALL( SCIPallocBufferArray(scip, &result, nedges) );
3867 else
3868 result = edgearrint;
3869
3870 SCIP_CALL( SCIPallocBufferArray(scip, &transresult, extnedges) );
3871 SCIP_CALL( SCIPallocBufferArray(scip, &marked, extnedges) );
3872 SCIP_CALL( SCIPallocBufferArray(scip, &result2, nedges) );
3873 SCIP_CALL( SCIPallocBufferArray(scip, &bestcost, extnedges) );
3874 SCIP_CALL( SCIPallocBufferArray(scip, &edgefixingbounds, extnedges) );
3875 SCIP_CALL( SCIPallocBufferArray(scip, &nodefixingbounds, nnodes + 1) );
3876
3877 if( userec )
3878 SCIP_CALL( SCIPStpHeurRecInitPool(scip, &pool, nedges, SOLPOOL_SIZE) );
3879 else
3880 pool = NULL;
3881
3882 #ifndef NDEBUG
3883 {
3884 int nterms = 0;
3885 for( int i = 0; i < nnodes; i++ )
3886 if( graph->mark[i] )
3887 if( Is_term(graph->term[i]) )
3888 nterms++;
3889
3890 assert(nterms == (graph->terms - ((graph->stp_type != STP_RPCSPG)? 1 : 0)));
3891 }
3892 #endif
3893
3894 /*
3895 * 1. step: compute lower bound and reduced costs
3896 */
3897
3898 graph_pc_2trans(graph);
3899 offset = 0.0;
3900
3901 /* transform the problem to a real SAP */
3902 SCIP_CALL( graph_pc_getSap(scip, graph, &transgraph, &offset) );
3903
3904 /* initialize data structures for shortest paths */
3905 SCIP_CALL( graph_path_init(scip, transgraph) );
3906
3907 damaxdeviation = fastmode ? DAMAXDEVIATION_FAST : -1.0;
3908
3909 SCIP_CALL( SCIPStpDualAscent(scip, transgraph, cost, pathdist, &lpobjval, FALSE, FALSE,
3910 gnodearr, NULL, transresult, state, root, TRUE, damaxdeviation, nodearrchar) );
3911
3912 lpobjval += offset;
3913 bestlpobjval = lpobjval;
3914 BMScopyMemoryArray(bestcost, cost, extnedges);
3915
3916 /* compute first primal solution */
3917 upperbound = FARAWAY;
3918 apsol = FALSE;
3919 SCIP_CALL( computeDaSolPcMw(scip, graph, NULL, vnoi, cost, pathdist, &upperbound, result, result2, vbase, pathedge, nodearrchar, &apsol) );
3920
3921 assert(apsol && upperbound < FARAWAY);
3922 assert(graph_sol_valid(scip, graph, result));
3923
3924 /* the required reduced path cost to be surpassed */
3925 minpathcost = upperbound - lpobjval;
3926 if( userec)
3927 SCIPdebugMessage("DA first minpathcost %f \n", minpathcost);
3928
3929 /* initialize data structures for transgraph */
3930 SCIP_CALL( graph_init_history(scip, transgraph) );
3931 computeTransVoronoi(scip, transgraph, vnoi, cost, costrev, pathdist, vbase, pathedge);
3932
3933 /*
3934 * 2. step: try to eliminate
3935 */
3936
3937 /* restore original graph */
3938 graph_pc_2org(graph);
3939
3940 for( int e = 0; e < extnedges; e++ )
3941 marked[e] = FALSE;
3942
3943 /* try to reduce the graph */
3944 assert(dualCostIsValid(scip, transgraph, cost, state, nodearrchar));
3945
3946 updateEdgeFixingBounds(edgefixingbounds, graph, cost, pathdist, vnoi, lpobjval, extnedges, TRUE, FALSE);
3947 updateNodeFixingBounds(nodefixingbounds, graph, pathdist, vnoi, lpobjval, TRUE);
3948
3949 nfixed += reducePcMw(scip, graph, transgraph, vnoi, cost, pathdist, minpathcost, result, marked, nodearrchar, TRUE);
3950
3951 /* edges from result might have been deleted! */
3952 apsol = apsol && graph_sol_unreduced(scip, graph, result);
3953 assert(!apsol || graph_sol_valid(scip, graph, result));
3954
3955 if( userec )
3956 SCIPdebugMessage("DA: 1. NFIXED %d \n", nfixed);
3957
3958 /* rerun dual ascent? */
3959 if( solbasedda && graph->terms > 2 && SCIPisGT(scip, minpathcost, 0.0) )
3960 {
3961 const SCIP_Real oldupperbound = upperbound;
3962
3963 /* with recombination? */
3964 if( userec && graph->stp_type != STP_MWCSP )
3965 {
3966 /* compute second solution and add to pool */
3967 SCIP_CALL( SCIPStpHeurTMRun(scip, NULL, graph, NULL, NULL, result2, DEFAULT_HEURRUNS / 5, root, graph->cost, graph->cost, NULL, NULL, 0.0, &success, FALSE) );
3968 assert(success);
3969
3970 SCIP_CALL( SCIPStpHeurLocalRun(scip, graph, graph->cost, result2) );
3971 ub = graph_sol_getObj(graph->cost, result2, 0.0, nedges);
3972
3973 SCIP_CALL( SCIPStpHeurRecAddToPool(scip, ub, result2, pool, &success) );
3974 SCIPdebugMessage("added initial TM sol to pool? %d , ub %f \n", success, ub);
3975 }
3976
3977 /* try to improve both dual and primal bound */
3978 SCIP_CALL( computePertubedSol(scip, graph, transgraph, pool, vnoi, gnodearr, cost, costrev, bestcost, pathdist, state, vbase, pathedge, result, result2,
3979 transresult, nodearrchar, &upperbound, &lpobjval, &bestlpobjval, &minpathcost, &apsol, offset, extnedges, 0) );
3980
3981 assert(graph_sol_valid(scip, graph, result));
3982 assert(!apsol || SCIPisEQ(scip, graph_sol_getObj(graph->cost, result, 0.0, nedges), upperbound));
3983
3984 graph_pc_2orgcheck(graph);
3985
3986 assert(dualCostIsValid(scip, transgraph, cost, state, nodearrchar));
3987 computeTransVoronoi(scip, transgraph, vnoi, cost, costrev, pathdist, vbase, pathedge);
3988 updateEdgeFixingBounds(edgefixingbounds, graph, cost, pathdist, vnoi, lpobjval, extnedges, FALSE, FALSE);
3989 updateNodeFixingBounds(nodefixingbounds, graph, pathdist, vnoi, lpobjval, FALSE);
3990
3991 nfixed += reducePcMw(scip, graph, transgraph, vnoi, cost, pathdist, minpathcost, result, marked, nodearrchar, apsol);
3992
3993 nfixed += reducePcMwTryBest(scip, graph, transgraph, vnoi, cost, costrev, bestcost, pathdist, &upperbound,
3994 &lpobjval, &bestlpobjval, &minpathcost, oldupperbound, result, vbase, state, pathedge, marked, nodearrchar, &apsol, extnedges);
3995
3996 nfixed += reduceWithEdgeFixingBounds(scip, graph, transgraph, edgefixingbounds, NULL, upperbound);
3997 nfixed += reduceWithNodeFixingBounds(scip, graph, transgraph, nodefixingbounds, upperbound);
3998
3999 if( userec )
4000 SCIPdebugMessage("eliminations after sol based run2 with best dual sol %d bestlb %f newlb %f\n", nfixed, bestlpobjval, lpobjval);
4001 }
4002
4003 /* pertubation runs for MWCSP */
4004 if( varyroot && graph->stp_type == STP_MWCSP )
4005 {
4006 for( int run = 0; run < DEFAULT_NMAXROOTS && graph->terms > STP_RED_MINBNDTERMS; run++ )
4007 {
4008 SCIP_Real oldupperbound = upperbound;
4009
4010 graph_pc_2trans(graph);
4011
4012 apsol = apsol && graph_sol_unreduced(scip, graph, result);
4013 assert(!apsol || graph_sol_valid(scip, graph, result));
4014
4015 assert(SCIPisEQ(scip, upperbound, graph_sol_getObj(graph->cost, result, 0.0, nedges)));
4016
4017 /* try to improve both dual and primal bound */
4018 SCIP_CALL( computePertubedSol(scip, graph, transgraph, pool, vnoi, gnodearr, cost, costrev, bestcost, pathdist, state, vbase, pathedge, result, result2,
4019 transresult, nodearrchar, &upperbound, &lpobjval, &bestlpobjval, &minpathcost, &apsol, offset, extnedges, run) );
4020
4021 SCIPdebugMessage("DA: pertubated run %d ub: %f \n", run, upperbound);
4022 SCIPdebugMessage("DA: pertubated run %d minpathcost: %f \n", run, upperbound - lpobjval);
4023
4024 computeTransVoronoi(scip, transgraph, vnoi, cost, costrev, pathdist, vbase, pathedge);
4025 updateEdgeFixingBounds(edgefixingbounds, graph, cost, pathdist, vnoi, lpobjval, extnedges, FALSE, FALSE);
4026 updateNodeFixingBounds(nodefixingbounds, graph, pathdist, vnoi, lpobjval, FALSE);
4027
4028 nfixed += reducePcMw(scip, graph, transgraph, vnoi, cost, pathdist, minpathcost, result, marked, nodearrchar, apsol);
4029
4030 nfixed += reducePcMwTryBest(scip, graph, transgraph, vnoi, cost, costrev, bestcost, pathdist, &upperbound,
4031 &lpobjval, &bestlpobjval, &minpathcost, oldupperbound, result, vbase, state, pathedge, marked, nodearrchar, &apsol, extnedges);
4032
4033 nfixed += reduceWithEdgeFixingBounds(scip, graph, transgraph, edgefixingbounds, NULL, upperbound);
4034 nfixed += reduceWithNodeFixingBounds(scip, graph, transgraph, nodefixingbounds, upperbound);
4035
4036 SCIPdebugMessage("DA: pertubated run %d NFIXED %d \n", run, nfixed);
4037 }
4038 }
4039
4040 if( graph->stp_type == STP_MWCSP && graph->terms < STP_RED_MINBNDTERMS )
4041 varyroot = FALSE;
4042
4043 /* change or mark roots? */
4044 if( varyroot || markroots )
4045 {
4046 SCIP_CALL( SCIPallocBufferArray(scip, &roots, graph->terms) );
4047
4048 findDaRoots(scip, graph, transgraph, cost, bestcost, lpobjval, bestlpobjval, upperbound, prizesum, FALSE, FALSE,
4049 vnoi, state, pathedge, vbase, nodearrchar, roots, &nroots);
4050
4051 /* should prize of terminals be changed? */
4052 if( nroots > 0 && markroots )
4053 markPcMwRoots(scip, roots, 0, nroots, prizesum, graph, &userec, &pool);
4054
4055 if( nroots > 0 && varyroot )
4056 SCIP_CALL( orderDaRoots(scip, graph, roots, nroots, TRUE, randnumgen) );
4057 }
4058 else
4059 {
4060 roots = NULL;
4061 }
4062
4063 if( varyroot )
4064 nusedroots = MIN(DEFAULT_NMAXROOTS, nroots);
4065 else
4066 nusedroots = -1;
4067
4068 graph_path_exit(scip, transgraph);
4069 graph_free(scip, &transgraph, TRUE);
4070
4071 /* loop and change root for dual ascent run */
4072 for( int run = 0; run < nusedroots; run++ )
4073 {
4074 const int tmproot = roots[run];
4075 int transnnodes;
4076 int transnedges;
4077
4078 assert(nroots > 0);
4079 assert(roots != NULL);
4080 assert(Is_term(graph->term[tmproot]));
4081
4082 if( graph->terms <= 2 )
4083 break;
4084
4085 SCIP_CALL( graph_pc_getRsap(scip, graph, &transgraph, roots, nroots, tmproot) );
4086
4087 assert(graph_valid(transgraph));
4088
4089 transnnodes = transgraph->knots;
4090 transnedges = transgraph->edges;
4091
4092 for( int k = 0; k < transnnodes; k++ )
4093 transgraph->mark[k] = (transgraph->grad[k] > 0);
4094
4095 /* init data structures for shortest paths and history */
4096 SCIP_CALL( graph_path_init(scip, transgraph) );
4097 SCIP_CALL( graph_init_history(scip, transgraph ) );
4098
4099 transgraph->stp_type = STP_SAP;
4100
4101 if( apsol && run > 1 )
4102 {
4103 BMScopyMemoryArray(transresult, result, graph->edges);
4104 SCIP_CALL(graph_sol_reroot(scip, transgraph, transresult, tmproot));
4105 SCIP_CALL( SCIPStpDualAscent(scip, transgraph, cost, pathdist, &lpobjval, FALSE, FALSE,
4106 gnodearr, transresult, result2, state, tmproot, FALSE, -1.0, nodearrchar));
4107 }
4108 else
4109 {
4110 SCIP_CALL( SCIPStpDualAscent(scip, transgraph, cost, pathdist, &lpobjval, FALSE, FALSE,
4111 gnodearr, NULL, transresult, state, tmproot, FALSE, -1.0, nodearrchar));
4112 }
4113
4114 assert(graph_valid(transgraph));
4115
4116 for( int e = graph->outbeg[tmproot]; e != EAT_LAST; e = graph->oeat[e] )
4117 {
4118 const int k = graph->head[e];
4119 if( Is_term(graph->term[k]) )
4120 {
4121 if( k == root )
4122 cost[flipedge(e)] = 0.0;
4123 else
4124 cost[e] = 0.0;
4125 }
4126 }
4127
4128 for( int k = 0; k < nnodes; k++ )
4129 {
4130 if( Is_pterm(graph->term[k]) && graph->prize[k] >= prizesum )
4131 {
4132 for( int e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
4133 {
4134 const int head = graph->head[e];
4135 if( Is_term(graph->term[head]) && head != root )
4136 cost[e] = 0.0;
4137 }
4138 }
4139 }
4140
4141 apsol = FALSE;
4142 SCIP_CALL( computeDaSolPcMw(scip, graph, pool, vnoi, cost, pathdist, &upperbound, result, result2, vbase, pathedge, nodearrchar, &apsol) );
4143
4144 SCIPdebugMessage("ROOTRUNS upperbound %f \n", upperbound);
4145 SCIPdebugMessage("ROOTRUNS best sol in pool %f \n", pool->sols[0]->obj);
4146
4147 for( int k = 0; k < transnnodes; k++ )
4148 transgraph->mark[k] = (transgraph->grad[k] > 0);
4149
4150 /* the required reduced path cost to be surpassed */
4151 minpathcost = upperbound - lpobjval;
4152
4153 if( markroots )
4154 {
4155 const int oldnroots = nroots;
4156 findDaRoots(scip, graph, transgraph, cost, cost, lpobjval, lpobjval, upperbound, prizesum, TRUE, TRUE,
4157 vnoi, state, pathedge, vbase, nodearrchar, roots, &nroots);
4158
4159 /* should prize of terminals be changed? */
4160 if( nroots > oldnroots )
4161 markPcMwRoots(scip, roots, oldnroots, nroots, prizesum, graph, &userec, &pool);
4162 }
4163
4164 SCIPdebugMessage("ROOTRUNS: minpathcost %f \n", minpathcost);
4165 SCIPdebugMessage("lb: %f ub: %f \n", lpobjval, upperbound);
4166
4167 /* distance from root to all nodes */
4168 graph_path_execX(scip, transgraph, tmproot, cost, pathdist, pathedge);
4169
4170 for( int e = 0; e < transnedges; e++ )
4171 costrev[e] = cost[flipedge(e)];
4172
4173 /* no paths should go back to the root */
4174 for( int e = transgraph->outbeg[tmproot]; e != EAT_LAST; e = transgraph->oeat[e] )
4175 costrev[e] = FARAWAY;
4176
4177 /* build Voronoi diagram */
4178 graph_voronoiTerms(scip, transgraph, costrev, vnoi, vbase, transgraph->path_heap, state);
4179 graph_get2next(scip, transgraph, costrev, costrev, vnoi, vbase, transgraph->path_heap, state);
4180
4181 /* restore original graph */
4182 graph_pc_2org(graph);
4183
4184 assert(graph->mark[tmproot]);
4185 graph->mark[tmproot] = FALSE;
4186
4187 /* at first run switch to undirected case */
4188 if( run == 0 )
4189 for( int e = 0; e < extnedges; e++ )
4190 edgefixingbounds[e] = MIN(edgefixingbounds[e], edgefixingbounds[flipedge(e)]);
4191
4192 updateEdgeFixingBounds(edgefixingbounds, graph, cost, pathdist, vnoi, lpobjval, extnedges, FALSE, TRUE);
4193 updateNodeFixingBounds(nodefixingbounds, graph, pathdist, vnoi, lpobjval, FALSE);
4194
4195 for( int e = 0; e < transnedges; e++ )
4196 marked[e] = FALSE;
4197
4198 /* try to eliminate vertices and edges */
4199 nfixed += reducePcMw(scip, graph, transgraph, vnoi, cost, pathdist, minpathcost, result, marked, nodearrchar, apsol);
4200 nfixed += reduceWithEdgeFixingBounds(scip, graph, transgraph, edgefixingbounds, NULL, upperbound);
4201 nfixed += reduceWithNodeFixingBounds(scip, graph, transgraph, nodefixingbounds, upperbound);
4202
4203 apsol = apsol && graph_sol_unreduced(scip, graph, result);
4204 assert(!apsol || graph_sol_valid(scip, graph, result));
4205 SCIPdebugMessage("FIXED with changed root %d \n\n", nfixed);
4206
4207 graph->mark[tmproot] = TRUE;
4208
4209 transgraph->stp_type = STP_RPCSPG;
4210 graph_path_exit(scip, transgraph);
4211 graph_free(scip, &transgraph, TRUE);
4212 }
4213
4214 *nelims = nfixed;
4215
4216 /* free memory */
4217 if( pool != NULL )
4218 SCIPStpHeurRecFreePool(scip, &pool);
4219 SCIPfreeBufferArrayNull(scip, &roots);
4220 SCIPfreeBufferArray(scip, &nodefixingbounds);
4221 SCIPfreeBufferArray(scip, &edgefixingbounds);
4222 SCIPfreeBufferArray(scip, &bestcost);
4223 SCIPfreeBufferArray(scip, &result2);
4224 SCIPfreeBufferArray(scip, &marked);
4225 SCIPfreeBufferArray(scip, &transresult);
4226
4227 if( edgearrint == NULL )
4228 SCIPfreeBufferArray(scip, &result);
4229
4230 assert(graph_valid(graph));
4231
4232 return SCIP_OKAY;
4233 }
4234
4235
4236 /** dual ascent based heuristic reductions for MWCSP */
reduce_daSlackPruneMw(SCIP * scip,GRAPH * graph,PATH * vnoi,GNODE ** gnodearr,SCIP_Real * cost,SCIP_Real * costrev,SCIP_Real * pathdist,int * vbase,int * pathedge,int * soledge,int * state,int * solnode,STP_Bool * nodearrchar,int * nelims,int minelims,SCIP_Bool solgiven)4237 SCIP_RETCODE reduce_daSlackPruneMw(
4238 SCIP* scip, /**< SCIP data structure */
4239 GRAPH* graph, /**< graph data structure */
4240 PATH* vnoi, /**< Voronoi data structure array */
4241 GNODE** gnodearr, /**< GNODE* terminals array for internal computations or NULL */
4242 SCIP_Real* cost, /**< edge costs */
4243 SCIP_Real* costrev, /**< reverse edge costs */
4244 SCIP_Real* pathdist, /**< distance array for shortest path calculations */
4245 int* vbase, /**< Voronoi base array */
4246 int* pathedge, /**< shortest path incoming edge array for shortest path calculations */
4247 int* soledge, /**< edge solution array (CONNECT/UNKNOWN) or NULL; needs to contain solution if solgiven == TRUE */
4248 int* state, /**< int 4 * vertices array */
4249 int* solnode, /**< array of nodes of current solution that is not to be destroyed */
4250 STP_Bool* nodearrchar, /**< STP_Bool node array for internal computations */
4251 int* nelims, /**< pointer to store number of reduced edges */
4252 int minelims, /**< minimum number of edges to eliminate */
4253 SCIP_Bool solgiven /**< solution provided? */
4254 )
4255 {
4256 SCIP_HEURDATA* tmheurdata;
4257 IDX** ancestors;
4258 IDX** revancestors;
4259 GRAPH* transgraph;
4260 SCIP_Real* sd;
4261 SCIP_Real* ecost;
4262 SCIP_Real ub;
4263 SCIP_Real offset;
4264 SCIP_Bool success;
4265 SCIP_Real tmpcost;
4266 SCIP_Real lpobjval;
4267 SCIP_Real hopfactor;
4268 SCIP_Real upperbound;
4269 SCIP_Real minpathcost;
4270 SCIP_Bool eliminate;
4271 int* result;
4272 int* adjvert;
4273 int* incedge;
4274 int* reinsert;
4275 int* transresult;
4276 int i;
4277 int k;
4278 int e;
4279 int e2;
4280 int etmp;
4281 int runs;
4282 int root;
4283 int nterms;
4284 int nedges;
4285 int nnodes;
4286 int nfixed;
4287 int redrounds;
4288 int best_start;
4289 int transnnodes;
4290 int transnedges;
4291 STP_Bool* marked;
4292
4293 assert(scip != NULL);
4294 assert(graph != NULL);
4295 assert(nelims != NULL);
4296 assert(nodearrchar != NULL);
4297
4298 root = graph->source;
4299 nfixed = 0;
4300 nnodes = graph->knots;
4301 nedges = graph->edges;
4302 success = FALSE;
4303 hopfactor = DEFAULT_HOPFACTOR;
4304
4305 /* not more than two terminals? */
4306 if( graph->terms <= 3 )
4307 return SCIP_OKAY;
4308
4309 /* allocate memory */
4310 if( soledge == NULL )
4311 {
4312 SCIP_CALL( SCIPallocBufferArray(scip, &result, nedges) );
4313 }
4314 else
4315 {
4316 result = soledge;
4317 }
4318
4319 SCIP_CALL( SCIPallocBufferArray(scip, &transresult, nedges + 2 * (graph->terms - 1)) );
4320 SCIP_CALL( SCIPallocBufferArray(scip, &marked, nedges + 2 * (graph->terms - 1)) );
4321
4322 /* allocate length-4 buffer memory */
4323 SCIP_CALL( SCIPallocBufferArray(scip, &sd, 4) );
4324 SCIP_CALL( SCIPallocBufferArray(scip, &ancestors, 4) );
4325 SCIP_CALL( SCIPallocBufferArray(scip, &revancestors, 4) );
4326 SCIP_CALL( SCIPallocBufferArray(scip, &ecost, 4) );
4327 SCIP_CALL( SCIPallocBufferArray(scip, &adjvert, 4) );
4328 SCIP_CALL( SCIPallocBufferArray(scip, &reinsert, 4) );
4329 SCIP_CALL( SCIPallocBufferArray(scip, &incedge, 4) );
4330
4331 for( i = 0; i < 4; i++ )
4332 {
4333 sd[i] = 0.0;
4334 ancestors[i] = NULL;
4335 revancestors[i] = NULL;
4336 }
4337
4338 /* 1. step: compute upper bound */
4339
4340 runs = 0;
4341 nterms = 0;
4342 for( i = 0; i < nnodes; i++ )
4343 {
4344 if( graph->grad[i] > 0 )
4345 {
4346 runs++;
4347 if( Is_term(graph->term[i]) )
4348 nterms++;
4349 }
4350 }
4351
4352 assert(nterms == (graph->terms - ((graph->stp_type != STP_RPCSPG)? 1 : 0)));
4353
4354 offset = 0.0;
4355
4356 /* transform the problem to a real SAP */
4357 SCIP_CALL( graph_pc_getSap(scip, graph, &transgraph, &offset) );
4358
4359 /* initialize data structures for shortest paths and history */
4360 SCIP_CALL( graph_path_init(scip, transgraph) );
4361 SCIP_CALL( graph_init_history(scip, transgraph ) );
4362 transnnodes = transgraph->knots;
4363 transnedges = transgraph->edges;
4364
4365 if( !solgiven )
4366 {
4367 for( k = 0; k < nnodes; k++ )
4368 nodearrchar[k] = (solnode[k] == CONNECT);
4369
4370 for( e = 0; e < nedges; e++ )
4371 {
4372 cost[e] = graph->cost[e];
4373 costrev[e] = graph->cost[flipedge(e)];
4374 result[e] = UNKNOWN;
4375 }
4376
4377 /* build trivial solution */
4378 SCIP_CALL( SCIPStpHeurTMPrunePc(scip, graph, graph->cost, result, nodearrchar) );
4379 }
4380 else
4381 {
4382 for( e = 0; e < nedges; e++ )
4383 {
4384 cost[e] = graph->cost[e];
4385 costrev[e] = graph->cost[flipedge(e)];
4386 }
4387 }
4388
4389 upperbound = graph_sol_getObj(graph->cost, result, 0.0, nedges);
4390
4391 /* number of runs should not exceed number of connected vertices */
4392 runs = MIN(runs, DEFAULT_HEURRUNS);
4393
4394 /* get TM heuristic data */
4395 assert(SCIPfindHeur(scip, "TM") != NULL);
4396 tmheurdata = SCIPheurGetData(SCIPfindHeur(scip, "TM"));
4397
4398
4399 /* compute Steiner tree to obtain upper bound */
4400 SCIP_CALL( SCIPStpHeurTMRun(scip, tmheurdata, graph, NULL, &best_start, transresult, runs, root, cost, costrev, &hopfactor, NULL, 0.0, &success, FALSE) );
4401
4402 /* feasible solution found? */
4403 if( success )
4404 {
4405 /* calculate objective value of solution */
4406 ub = graph_sol_getObj(graph->cost, transresult, 0.0, nedges);
4407
4408 SCIPdebugMessage("TM upperbound in reduce_daSlackPruneMw %f \n", ub + SCIPprobdataGetOffset(scip));
4409
4410 if( SCIPisLE(scip, ub, upperbound) )
4411 {
4412 upperbound = ub;
4413 for( e = 0; e < nedges; e++ )
4414 result[e] = transresult[e];
4415 }
4416 }
4417
4418 /* compute lower bound and reduced costs todo use SCIPdualAscentStpSol */
4419 SCIP_CALL( SCIPStpDualAscent(scip, transgraph, cost, pathdist, &lpobjval, FALSE, FALSE, gnodearr, NULL, transresult, state, root, FALSE, -1.0, nodearrchar) );
4420
4421 SCIP_CALL( SCIPStpHeurAscendPruneRun(scip, NULL, graph, cost, transresult, vbase, -1, nodearrchar, &success, FALSE) );
4422
4423 assert(success);
4424 assert(graph_sol_valid(scip, graph, transresult));
4425
4426 if( success )
4427 {
4428 ub = graph_sol_getObj(graph->cost, transresult, 0.0, nedges);
4429
4430 SCIPdebugMessage("AP upperbound in reduce_daSlackPruneMw %f \n", ub + SCIPprobdataGetOffset(scip));
4431
4432 if( SCIPisLE(scip, ub, upperbound) )
4433 for( e = 0; e < nedges; e++ )
4434 result[e] = transresult[e];
4435 }
4436
4437 /*
4438 * 2. step: try to eliminate
4439 * */
4440
4441 for( e = 0; e < transnedges; e++ )
4442 {
4443 costrev[e] = cost[flipedge(e)];
4444 marked[e] = FALSE;
4445 }
4446
4447 for( k = 0; k < transnnodes; k++ )
4448 transgraph->mark[k] = (transgraph->grad[k] > 0);
4449
4450 /* distance from root to all nodes */
4451 graph_path_execX(scip, transgraph, root, cost, pathdist, pathedge);
4452
4453 for( i = 0; i < transnnodes; i++ )
4454 if( Is_term(transgraph->term[i]) )
4455 assert(SCIPisEQ(scip, pathdist[i], 0.0));
4456
4457 /* no paths should go back to the root */
4458 for( e = transgraph->outbeg[root]; e != EAT_LAST; e = transgraph->oeat[e] )
4459 costrev[e] = FARAWAY;
4460
4461 /* build Voronoi diagram */
4462 graph_voronoiTerms(scip, transgraph, costrev, vnoi, vbase, transgraph->path_heap, transgraph->path_state);
4463
4464 /* restore original graph */
4465 graph_pc_2org(graph);
4466
4467 for( k = 0; k < nnodes; k++ )
4468 solnode[k] = UNKNOWN;
4469
4470 for( e = 0; e < nedges; e++ )
4471 {
4472 costrev[e] = -1.0;
4473 if( result[e] == CONNECT )
4474 {
4475 solnode[graph->head[e]] = CONNECT;
4476 solnode[graph->tail[e]] = CONNECT;
4477 }
4478 }
4479
4480 for( redrounds = 0; redrounds < 3; redrounds++ )
4481 {
4482 if( redrounds == 0 )
4483 {
4484 eliminate = FALSE;
4485 minpathcost = FARAWAY;
4486 }
4487 else if( redrounds == 1 )
4488 {
4489 assert(minelims > 0);
4490 assert(2 * minelims < nedges);
4491
4492 eliminate = TRUE;
4493 SCIPsortReal(costrev, nedges);
4494
4495 /* the required reduced path cost to be surpassed */
4496 minpathcost = costrev[nedges - 2 * minelims];
4497
4498 if( SCIPisLE(scip, minpathcost, 0.0) )
4499 minpathcost = 2 * SCIPepsilon(scip);
4500
4501 k = nedges - 2 * minelims;
4502
4503 /* try to first eliminate edges with higher gap */
4504 for( e = nedges - 1; e > k && e >= 2; e = e - 2 )
4505 {
4506 if( SCIPisLE(scip, costrev[e - 2], minpathcost) )
4507 break;
4508 }
4509
4510 if( SCIPisGT(scip, costrev[e], minpathcost) )
4511 minpathcost = costrev[e];
4512
4513 if( SCIPisLE(scip, minpathcost, 0.0) )
4514 minpathcost = 2 * SCIPepsilon(scip);
4515
4516
4517 for( e = 0; e < nedges; e++ )
4518 marked[e] = FALSE;
4519 }
4520 else
4521 {
4522 eliminate = TRUE;
4523
4524 /* the required reduced path cost to be surpassed */
4525 minpathcost = costrev[nedges - 2 * minelims];
4526
4527 if( SCIPisLE(scip, minpathcost, 0.0) )
4528 minpathcost = 2 * SCIPepsilon(scip);
4529
4530 for( e = 0; e < nedges; e++ )
4531 marked[e] = FALSE;
4532 }
4533
4534 /* try to eliminate vertices and edges */
4535 for( k = 0; k < nnodes; k++ )
4536 {
4537 if( !graph->mark[k] )
4538 continue;
4539
4540 if( nfixed > minelims )
4541 break;
4542
4543
4544 if( Is_term(graph->term[k]) )
4545 continue;
4546
4547 tmpcost = pathdist[k] + vnoi[k].dist;
4548
4549 if( SCIPisGT(scip, tmpcost, minpathcost) ||
4550 (SCIPisGE(scip, tmpcost, minpathcost) && solnode[k] != CONNECT) || (!eliminate && solnode[k] != CONNECT) )
4551 {
4552 if( !eliminate )
4553 {
4554 for( e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
4555 {
4556 if( SCIPisGT(scip, tmpcost, costrev[e]) )
4557 costrev[e] = tmpcost;
4558
4559 e2 = flipedge(e);
4560
4561 if( SCIPisGT(scip, tmpcost, costrev[e2]) )
4562 costrev[e2] = tmpcost;
4563 }
4564
4565 continue;
4566 }
4567
4568 while( transgraph->outbeg[k] != EAT_LAST )
4569 {
4570 e = transgraph->outbeg[k];
4571 graph_edge_del(scip, transgraph, e, FALSE);
4572 graph_edge_del(scip, graph, e, TRUE);
4573 nfixed++;
4574 }
4575 assert(graph->outbeg[k] == EAT_LAST);
4576 }
4577 else
4578 {
4579 e = transgraph->outbeg[k];
4580 while( e != EAT_LAST )
4581 {
4582 etmp = transgraph->oeat[e];
4583 tmpcost = pathdist[k] + cost[e] + vnoi[transgraph->head[e]].dist;
4584
4585 if( SCIPisGT(scip, tmpcost, minpathcost) ||
4586 ( (SCIPisGE(scip, tmpcost, minpathcost) || !eliminate) && result[e] != CONNECT && result[flipedge(e)] != CONNECT) )
4587 {
4588 if( marked[flipedge(e)] )
4589 {
4590 if( eliminate )
4591 {
4592 graph_edge_del(scip, graph, e, TRUE);
4593 graph_edge_del(scip, transgraph, e, FALSE);
4594 nfixed++;
4595 }
4596 else
4597 {
4598 if( SCIPisGT(scip, tmpcost, costrev[e]) )
4599 costrev[e] = tmpcost;
4600
4601 if( SCIPisGT(scip, tmpcost, costrev[flipedge(e)]) )
4602 costrev[flipedge(e)] = tmpcost;
4603 }
4604 }
4605 else
4606 {
4607 marked[e] = TRUE;
4608 }
4609 }
4610 e = etmp;
4611 }
4612 }
4613 }
4614 }
4615 assert(root == transgraph->source);
4616
4617 SCIPdebugMessage("fixed by da: %d \n", nfixed );
4618
4619 graph_path_exit(scip, transgraph);
4620 graph_free(scip, &transgraph, TRUE);
4621
4622 /* restore transformed graph */
4623 graph_pc_2trans(graph);
4624
4625 *nelims = nfixed;
4626
4627 /* free memory */
4628 for( k = 0; k < 4; k++ )
4629 {
4630 SCIPintListNodeFree(scip, &(ancestors[k]));
4631 SCIPintListNodeFree(scip, &(revancestors[k]));
4632 }
4633
4634 SCIPfreeBufferArray(scip, &incedge);
4635 SCIPfreeBufferArray(scip, &reinsert);
4636 SCIPfreeBufferArray(scip, &adjvert);
4637 SCIPfreeBufferArray(scip, &ecost);
4638 SCIPfreeBufferArray(scip, &revancestors);
4639 SCIPfreeBufferArray(scip, &ancestors);
4640 SCIPfreeBufferArray(scip, &sd);
4641 SCIPfreeBufferArray(scip, &marked);
4642 SCIPfreeBufferArray(scip, &transresult);
4643
4644 if( soledge == NULL )
4645 SCIPfreeBufferArray(scip, &result);
4646
4647 return SCIP_OKAY;
4648 }
4649
4650
4651
4652 /** bound-based reductions for the (R)PCSTP, the MWCSP and the STP */
reduce_bound(SCIP * scip,GRAPH * graph,PATH * vnoi,SCIP_Real * cost,SCIP_Real * prize,SCIP_Real * radius,SCIP_Real * costrev,SCIP_Real * offset,SCIP_Real * upperbound,int * heap,int * state,int * vbase,int * nelims)4653 SCIP_RETCODE reduce_bound(
4654 SCIP* scip, /**< SCIP data structure */
4655 GRAPH* graph, /**< graph data structure */
4656 PATH* vnoi, /**< Voronoi data structure */
4657 SCIP_Real* cost, /**< edge cost array */
4658 SCIP_Real* prize, /**< prize (nodes) array */
4659 SCIP_Real* radius, /**< radius array */
4660 SCIP_Real* costrev, /**< reversed edge cost array */
4661 SCIP_Real* offset, /**< pointer to the offset */
4662 SCIP_Real* upperbound, /**< pointer to an upper bound */
4663 int* heap, /**< heap array */
4664 int* state, /**< array to store state of a node during Voronoi computation*/
4665 int* vbase, /**< Voronoi base to each node */
4666 int* nelims /**< pointer to store number of eliminated edges */
4667 )
4668 {
4669 SCIP_HEURDATA* tmheurdata;
4670 GRAPH* adjgraph;
4671 PATH* mst;
4672 SCIP_Real r;
4673 SCIP_Real obj;
4674 SCIP_Real max;
4675 SCIP_Real radii;
4676 SCIP_Real bound;
4677 SCIP_Real tmpcost;
4678 SCIP_Real mstobj;
4679 SCIP_Real mstobj2;
4680 SCIP_Real maxcost;
4681 SCIP_Real radiim2;
4682 SCIP_Real radiim3;
4683 int* perm;
4684 int* result;
4685 int* starts;
4686 int e;
4687 int k;
4688 int l;
4689 int head;
4690 int tail;
4691 int runs;
4692 int root;
4693 int etemp;
4694 int nterms;
4695 int nnodes;
4696 int nedges;
4697 int best_start;
4698 STP_Bool* stnode;
4699 SCIP_Bool ub;
4700 SCIP_Bool pc;
4701 SCIP_Bool mw;
4702 SCIP_Bool pcmw;
4703 SCIP_Bool success;
4704
4705 assert(scip != NULL);
4706 assert(graph != NULL);
4707 assert(vnoi != NULL);
4708 assert(cost != NULL);
4709 assert(radius != NULL);
4710 assert(costrev != NULL);
4711 assert(heap != NULL);
4712 assert(state != NULL);
4713 assert(vbase != NULL);
4714 assert(nelims != NULL);
4715 assert(graph->source >= 0);
4716 assert(upperbound != NULL);
4717
4718 mst = NULL;
4719 obj = DEFAULT_HOPFACTOR;
4720 perm = NULL;
4721 root = graph->source;
4722 nedges = graph->edges;
4723 nnodes = graph->knots;
4724 mstobj = 0.0;
4725 *nelims = 0;
4726 mstobj2 = 0.0;
4727 best_start = 0;
4728 ub = SCIPisGT(scip, *upperbound, 0.0);
4729 mw = (graph->stp_type == STP_MWCSP);
4730 pc = (graph->stp_type == STP_RPCSPG) || (graph->stp_type == STP_PCSPG);
4731 pcmw = pc || mw;
4732 success = TRUE;
4733
4734 /* no upper bound provided? */
4735 if( !ub )
4736 {
4737 /* allocate memory */
4738 SCIP_CALL( SCIPallocBufferArray(scip, &result, nedges) );
4739 SCIP_CALL( SCIPallocBufferArray(scip, &stnode, nnodes) );
4740 }
4741 else
4742 {
4743 result = NULL;
4744 stnode = NULL;
4745 }
4746
4747 /* initialize */
4748 e = 0;
4749 nterms = 0;
4750 for( k = 0; k < nnodes; k++ )
4751 {
4752 if( !ub )
4753 {
4754 assert(stnode != NULL);
4755 stnode[k] = FALSE;
4756 }
4757 if( !pcmw )
4758 graph->mark[k] = (graph->grad[k] > 0);
4759 if( graph->mark[k] )
4760 {
4761 e++;
4762 if( Is_term(graph->term[k]) )
4763 nterms++;
4764 }
4765 }
4766
4767 /* not more than two terminals? */
4768 if( nterms <= 2 || (pcmw && nterms <= 3) )
4769 {
4770 /* free memory and return */
4771 SCIPfreeBufferArrayNull(scip, &stnode);
4772 SCIPfreeBufferArrayNull(scip, &result);
4773 return SCIP_OKAY;
4774 }
4775
4776 assert(nterms == (graph->terms - ((graph->stp_type == STP_PCSPG || mw)? 1 : 0)));
4777
4778 runs = MIN(e, DEFAULT_HEURRUNS);
4779
4780 /* neither PC, MW, RPC nor HC? */
4781 if( !pcmw && graph->stp_type != STP_DHCSTP )
4782 {
4783 /* choose starting points for TM heuristic */
4784
4785 SCIP_CALL( SCIPallocBufferArray(scip, &starts, nnodes) );
4786
4787 SCIPStpHeurTMCompStarts(graph, starts, &runs);
4788 }
4789 else
4790 {
4791 starts = NULL;
4792 }
4793
4794 /* initialise cost and costrev array */
4795 maxcost = 0.0;
4796 for( e = 0; e < nedges; e++ )
4797 {
4798 if( !ub )
4799 {
4800 assert(result != NULL);
4801 result[e] = UNKNOWN;
4802 }
4803 cost[e] = graph->cost[e];
4804 costrev[e] = graph->cost[flipedge(e)];
4805
4806 assert(SCIPisGE(scip, cost[e], 0.0));
4807
4808 if( graph->stp_type == STP_DHCSTP && SCIPisLT(scip, graph->cost[e], BLOCKED) && SCIPisGT(scip, graph->cost[e], maxcost) )
4809 maxcost = graph->cost[e];
4810 }
4811
4812 /* init auxiliary graph */
4813 if( !mw )
4814 {
4815 SCIP_CALL( graph_init(scip, &adjgraph, nterms, MIN(nedges, (nterms - 1) * nterms), 1) );
4816 }
4817 else
4818 {
4819 adjgraph = NULL;
4820 }
4821
4822 /* build voronoi regions, concomitantly building adjgraph and computing radii values*/
4823 SCIP_CALL( graph_voronoiWithRadius(scip, graph, adjgraph, vnoi, radius, cost, costrev, vbase, heap, state) );
4824
4825 /* get 2nd next terminals to all non-terminal nodes */
4826 graph_get2next(scip, graph, cost, costrev, vnoi, vbase, heap, state);
4827
4828 /* get 3th next terminals to all non-terminal nodes */
4829 graph_get3next(scip, graph, cost, costrev, vnoi, vbase, heap, state);
4830
4831 /* for (rooted) prize collecting get 4th next terminals to all non-terminal nodes */
4832 if( pc )
4833 graph_get4next(scip, graph, cost, costrev, vnoi, vbase, heap, state);
4834
4835 /* no MWCS problem? */
4836 if( !mw )
4837 {
4838 assert(adjgraph != NULL);
4839 graph_knot_chg(adjgraph, 0, 0);
4840 adjgraph->source = 0;
4841
4842 /* compute MST on adjgraph */
4843 SCIP_CALL( SCIPallocBufferArray(scip, &mst, nterms) );
4844 SCIP_CALL( graph_path_init(scip, adjgraph) );
4845 graph_path_exec(scip, adjgraph, MST_MODE, 0, adjgraph->cost, mst);
4846
4847 max = -1.0;
4848 r = -1.0;
4849 for( k = 1; k < nterms; k++ )
4850 {
4851 assert(adjgraph->path_state[k] == CONNECT);
4852 e = mst[k].edge;
4853 assert(e >= 0);
4854 tmpcost = adjgraph->cost[e];
4855 mstobj += tmpcost;
4856 if( SCIPisGT(scip, tmpcost, max) )
4857 max = tmpcost;
4858 else if( SCIPisGT(scip, tmpcost, r) )
4859 r = tmpcost;
4860 }
4861 mstobj -= max;
4862 mstobj2 = mstobj - r;
4863 }
4864
4865 /* for (rooted) prize collecting an maximum weight problems: correct radius values */
4866 if( graph->stp_type == STP_RPCSPG )
4867 {
4868 assert(graph->mark[graph->source]);
4869 for( k = 0; k < nnodes; k++ )
4870 {
4871 if( !Is_term(graph->term[k]) || !graph->mark[k] )
4872 continue;
4873
4874 if( SCIPisGT(scip, radius[k], prize[k]) && k != graph->source )
4875 radius[k] = prize[k];
4876 }
4877 }
4878 else if( graph->stp_type == STP_PCSPG )
4879 {
4880 for( k = 0; k < nnodes; k++ )
4881 {
4882 if( !graph->mark[k] )
4883 continue;
4884
4885 if( Is_term(graph->term[k]) && SCIPisGT(scip, radius[k], prize[k]) )
4886 radius[k] = prize[k];
4887 }
4888 }
4889 else if( graph->stp_type == STP_MWCSP )
4890 {
4891 max = 0.0;
4892 for( k = 0; k < nnodes; k++ )
4893 {
4894 if( !Is_term(graph->term[k]) || !graph->mark[k] )
4895 continue;
4896
4897 assert(SCIPisGE(scip, prize[k], 0.0));
4898 if( SCIPisGT(scip, prize[k], max) )
4899 max = prize[k];
4900
4901 if( SCIPisGE(scip, radius[k], 0.0 ) )
4902 radius[k] = 0.0;
4903 else
4904 radius[k] = -radius[k];
4905 }
4906 }
4907
4908 /* sort radius values */
4909 if( pc )
4910 {
4911 SCIP_CALL( SCIPallocBufferArray(scip, &perm, nnodes) );
4912 for( k = 0; k < nnodes; k++ )
4913 perm[k] = k;
4914 SCIPsortRealInt(radius, perm, nnodes);
4915 }
4916 else
4917 {
4918 SCIPsortReal(radius, nnodes);
4919 }
4920
4921 radiim2 = 0.0;
4922
4923 /* sum all but two radius values of highest/lowest value */
4924 if( mw )
4925 {
4926 for( k = 2; k < nterms; k++ )
4927 {
4928 assert( SCIPisGT(scip, FARAWAY, radius[k]) );
4929 radiim2 += radius[k];
4930 }
4931 }
4932 else
4933 {
4934 for( k = 0; k < nterms - 2; k++ )
4935 {
4936 assert( SCIPisGT(scip, FARAWAY, radius[k]) );
4937 radiim2 += radius[k];
4938 }
4939 }
4940 radii = radiim2 + radius[nterms - 2] + radius[nterms - 1];
4941 #if 1
4942 if( nterms >= 3 )
4943 radiim3 = radiim2 - radius[nterms - 3];
4944 else
4945 radiim3 = 0;
4946 #endif
4947 /* get TM heuristic data */
4948 assert(SCIPfindHeur(scip, "TM") != NULL);
4949 tmheurdata = SCIPheurGetData(SCIPfindHeur(scip, "TM"));
4950
4951 /* no upper bound available? */
4952 if( !ub )
4953 {
4954 assert(stnode != NULL);
4955 assert(result != NULL);
4956
4957 /* PC or RPC? Then restore transformed graph */
4958 if( pcmw )
4959 graph_pc_2trans(graph);
4960
4961 SCIP_CALL( SCIPStpHeurTMRun(scip, tmheurdata, graph, starts, &best_start, result, runs, root, cost, costrev, &obj, NULL, maxcost, &success, FALSE) );
4962
4963 /* PC or RPC? Then restore oringinal graph */
4964 if( pcmw )
4965 graph_pc_2org(graph);
4966
4967 /* no feasible solution found? */
4968 if( !success )
4969 {
4970 /* free memory and return */
4971 if( !mw )
4972 {
4973 graph_path_exit(scip, adjgraph);
4974 graph_free(scip, &adjgraph, TRUE);
4975 }
4976 SCIPfreeBufferArrayNull(scip, &mst);
4977 SCIPfreeBufferArrayNull(scip, &starts);
4978 SCIPfreeBufferArray(scip, &result);
4979 SCIPfreeBufferArray(scip, &stnode);
4980 return SCIP_OKAY;
4981 }
4982
4983 /* calculate objective value of solution */
4984 obj = 0.0;
4985 for( e = 0; e < nedges; e++ )
4986 {
4987 if( result[e] == CONNECT )
4988 {
4989 head = graph->head[e];
4990 if( mw )
4991 {
4992 if( graph->mark[head] )
4993 {
4994 assert(stnode[head] == FALSE);
4995 obj += prize[head];
4996 }
4997 }
4998 else
4999 {
5000 obj += graph->cost[e];
5001 stnode[head] = TRUE;
5002 stnode[graph->tail[e]] = TRUE;
5003 }
5004 }
5005 }
5006 *upperbound = obj;
5007 }
5008 else
5009 {
5010 obj = *upperbound;
5011 assert(SCIPisGE(scip, obj, 0.0));
5012 }
5013
5014 if( SCIPisGT(scip, radiim2, mstobj) )
5015 bound = radiim2;
5016 else
5017 bound = mstobj;
5018
5019 /* traverse all node, try to eliminate each node or incident edges */
5020 for( k = 0; k < nnodes; k++ )
5021 {
5022 if( (!graph->mark[k] && (pcmw)) || graph->grad[k] == 0 )
5023 continue;
5024
5025 if( mw )
5026 {
5027 tmpcost = -vnoi[k].dist - vnoi[k + nnodes].dist + bound - graph->prize[k];
5028
5029 if( !Is_term(graph->term[k]) && (SCIPisLT(scip, tmpcost, obj)) )
5030 {
5031 while( graph->outbeg[k] != EAT_LAST )
5032 {
5033 (*nelims)++;
5034 graph_edge_del(scip, graph, graph->outbeg[k], TRUE);
5035 }
5036 }
5037 else
5038 {
5039 e = graph->outbeg[k];
5040 while( e != EAT_LAST )
5041 {
5042 etemp = graph->oeat[e];
5043 tail = graph->tail[e];
5044 head = graph->head[e];
5045 if( !graph->mark[head] )
5046 {
5047 e = etemp;
5048 continue;
5049 }
5050 tmpcost = bound - graph->prize[k];
5051 if( vbase[tail] != vbase[head] )
5052 {
5053 tmpcost -= vnoi[head].dist + vnoi[tail].dist;
5054 }
5055 else
5056 {
5057 if( SCIPisGT(scip, vnoi[tail].dist + vnoi[head + nnodes].dist, vnoi[tail + nnodes].dist + vnoi[head].dist) )
5058 tmpcost -= vnoi[tail + nnodes].dist + vnoi[head].dist;
5059 else
5060 tmpcost -= vnoi[tail].dist + vnoi[head + nnodes].dist;
5061 }
5062 if( (SCIPisLT(scip, tmpcost, obj)) )
5063 {
5064 (*nelims)++;
5065 graph_edge_del(scip, graph, e, TRUE);
5066 }
5067 e = etemp;
5068 }
5069 }
5070 continue;
5071 }
5072
5073 tmpcost = vnoi[k].dist + vnoi[k + nnodes].dist + bound;
5074
5075 /* can node k be deleted? */
5076 if( !Is_term(graph->term[k]) && (SCIPisGT(scip, tmpcost, obj)
5077 || (((stnode != NULL) ? !stnode[k] : 0) && SCIPisGE(scip, tmpcost, obj))) )
5078 {
5079 /* delete all incident edges */
5080 while( graph->outbeg[k] != EAT_LAST )
5081 {
5082 e = graph->outbeg[k];
5083 (*nelims)++;
5084 assert(!pc || graph->tail[e] != root);
5085 assert(!pc || graph->mark[graph->head[e]]);
5086 assert(!Is_pterm(graph->term[graph->head[e]]));
5087 assert(!Is_pterm(graph->term[graph->tail[e]]));
5088
5089 graph_edge_del(scip, graph, e, TRUE);
5090 }
5091 }
5092 else
5093 {
5094 e = graph->outbeg[k];
5095 while( e != EAT_LAST )
5096 {
5097 etemp = graph->oeat[e];
5098 tail = graph->tail[e];
5099 head = graph->head[e];
5100 tmpcost = graph->cost[e] + bound;
5101
5102 if( vbase[tail] != vbase[head] )
5103 {
5104 tmpcost += vnoi[head].dist + vnoi[tail].dist;
5105 }
5106 else
5107 {
5108 if( SCIPisGT(scip, vnoi[tail].dist + vnoi[head + nnodes].dist, vnoi[tail + nnodes].dist + vnoi[head].dist) )
5109 tmpcost += vnoi[tail + nnodes].dist + vnoi[head].dist;
5110 else
5111 tmpcost += vnoi[tail].dist + vnoi[head + nnodes].dist;
5112 assert(SCIPisGE(scip, tmpcost, vnoi[head].dist + vnoi[tail].dist + graph->cost[e] + bound));
5113 }
5114 /* can edge e or arc e be deleted? */
5115 if( (SCIPisGT(scip, tmpcost, obj) || (((result != NULL) ? (result[e] != CONNECT) : 0) && result[flipedge(e)] != CONNECT && SCIPisGE(scip, tmpcost, obj)))
5116 && SCIPisLT(scip, graph->cost[e], FARAWAY) && (!(pc || mw) || graph->mark[head]) )
5117 {
5118 if( graph->stp_type == STP_DHCSTP && SCIPisGT(scip, graph->cost[e], graph->cost[flipedge(e)]) )
5119 {
5120 graph->cost[e] = FARAWAY;
5121 (*nelims)++;
5122 }
5123 else
5124 {
5125 assert(!Is_pterm(graph->term[head]));
5126 assert(!Is_pterm(graph->term[tail]));
5127 graph_edge_del(scip, graph, e, TRUE);
5128 (*nelims)++;
5129 }
5130 }
5131 e = etemp;
5132 }
5133 }
5134 }
5135 #if 1
5136 /* traverse all node, try to eliminate 3 degree nodes */
5137 if( !mw )
5138 {
5139 for( k = 0; k < nnodes; k++ )
5140 {
5141 if( (!graph->mark[k] && pc) || graph->grad[k] == 0 )
5142 continue;
5143
5144 if( (graph->grad[k] == 3 || graph->grad[k] == 4) && !Is_term(graph->term[k]) )
5145 {
5146 tmpcost = vnoi[k].dist + vnoi[k + nnodes].dist + vnoi[k + 2 * nnodes].dist + radiim3;
5147 if( SCIPisGT(scip, tmpcost, obj) )
5148 {
5149 SCIP_CALL( graph_knot_delPseudo(scip, graph, graph->cost, NULL, NULL, k, &success) );
5150
5151 assert(graph->grad[k] == 0 || (graph->grad[k] == 4 && !success));
5152 }
5153 }
5154 }
5155 }
5156 #endif
5157
5158 /* for(R)PC: try to eliminate terminals */
5159 if( pc )
5160 {
5161 SCIP_CALL( graph_get4nextTTerms(scip, graph, cost, vnoi, vbase, heap, state) );
5162
5163 for( k = 0; k < nnodes; k++ )
5164 {
5165 /* is k a terminal other than the root? */
5166 if( Is_term(graph->term[k]) && graph->mark[k] && graph->grad[k] > 2 && k != graph->source )
5167 {
5168 assert(perm != NULL);
5169 for( l = 0; l < nterms; l++ )
5170 if( perm[l] == k )
5171 break;
5172
5173 tmpcost = vnoi[k].dist + radii - radius[l];
5174
5175 if( l == nterms - 1 )
5176 tmpcost -= radius[nterms - 2];
5177 else
5178 tmpcost -= radius[nterms - 1];
5179
5180
5181 if( SCIPisGT(scip, tmpcost, obj) )
5182 {
5183 SCIPdebugMessage("alternative bnd elimination!!! \n\n");
5184 (*nelims) += graph_pc_deleteTerm(scip, graph, k);
5185 (*offset) += graph->prize[k];
5186 }
5187 else
5188 {
5189 tmpcost += vnoi[k + nnodes].dist;
5190 if( l >= nterms - 2 )
5191 tmpcost -= radius[nterms - 3];
5192 else
5193 tmpcost -= radius[nterms - 2];
5194 if( SCIPisGT(scip, tmpcost, obj) || SCIPisGT(scip, mstobj2 + vnoi[k].dist + vnoi[k + nnodes].dist, obj) )
5195 {
5196 for( e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
5197 if( graph->mark[graph->head[e]] && SCIPisLT(scip, graph->cost[e], graph->prize[k]) )
5198 break;
5199
5200 if( e == EAT_LAST )
5201 {
5202 SCIPdebugMessage("second elimination!!! prize: %f \n\n", graph->prize[k]);
5203 (*offset) += graph->prize[k];
5204 (*nelims) += graph_pc_deleteTerm(scip, graph, k);
5205 }
5206 }
5207 }
5208 }
5209 }
5210 }
5211
5212 SCIPdebugMessage("nelims (edges) in bound reduce: %d,\n", *nelims);
5213
5214 /* free adjgraph */
5215 if( !mw )
5216 {
5217 graph_path_exit(scip, adjgraph);
5218 graph_free(scip, &adjgraph, TRUE);
5219 }
5220
5221 /* free memory*/
5222 SCIPfreeBufferArrayNull(scip, &perm);
5223 SCIPfreeBufferArrayNull(scip, &mst);
5224 SCIPfreeBufferArrayNull(scip, &starts);
5225 SCIPfreeBufferArrayNull(scip, &stnode);
5226 SCIPfreeBufferArrayNull(scip, &result);
5227
5228 assert(graph_valid(graph));
5229
5230 return SCIP_OKAY;
5231 }
5232
5233
5234
5235
5236 /** Bound-based reduction method for the MWCSP .
5237 * The reduction method tries to eliminate nodes and vertices
5238 * by checking whether an upper bound for each solution that contains them
5239 * is smaller than the best known solution value.
5240 * The essence of the approach is a decomposition of the graph such that this upper bound
5241 * is minimized.
5242 * */
reduce_boundMw(SCIP * scip,GRAPH * graph,PATH * vnoi,PATH * path,SCIP_Real * cost,SCIP_Real * radius,SCIP_Real * costrev,SCIP_Real * offset,int * heap,int * state,int * vbase,int * result,int * nelims)5243 SCIP_RETCODE reduce_boundMw(
5244 SCIP* scip, /**< SCIP data structure */
5245 GRAPH* graph, /**< graph data structure */
5246 PATH* vnoi, /**< Voronoi data structure (size 3 * nnodes) */
5247 PATH* path, /**< shortest path data structure (size nnodes) */
5248 SCIP_Real* cost, /**< edge cost array */
5249 SCIP_Real* radius, /**< radius array */
5250 SCIP_Real* costrev, /**< reversed edge cost array */
5251 SCIP_Real* offset, /**< pointer to the offset */
5252 int* heap, /**< heap array */
5253 int* state, /**< array to store state of a node during Voronoi computation*/
5254 int* vbase, /**< Voronoi base to each node */
5255 int* result, /**< solution array or NULL */
5256 int* nelims /**< pointer to store number of eliminated edges */
5257 )
5258 {
5259 PATH* mst;
5260 SCIP_Real* prize;
5261 SCIP_Real obj;
5262 SCIP_Real bound;
5263 SCIP_Real tmpcost;
5264 SCIP_Real radiim2;
5265 int e;
5266 int k;
5267 int head;
5268 int nterms;
5269 int nnodes;
5270 int nedges;
5271
5272 assert(scip != NULL);
5273 assert(graph != NULL);
5274 assert(vnoi != NULL);
5275 assert(path != NULL);
5276 assert(cost != NULL);
5277 assert(radius != NULL);
5278 assert(costrev != NULL);
5279 assert(heap != NULL);
5280 assert(state != NULL);
5281 assert(vbase != NULL);
5282 assert(nelims != NULL);
5283 assert(graph->source >= 0);
5284
5285 mst = NULL;
5286 prize = graph->prize;
5287 nedges = graph->edges;
5288 nnodes = graph->knots;
5289 nterms = graph->terms - 1;
5290 *nelims = 0;
5291
5292 assert(prize != NULL);
5293
5294 /* not more than two nodes of positive weight? */
5295 if( nterms <= 2 )
5296 /* return */
5297 return SCIP_OKAY;
5298
5299 /* initialize cost and costrev array */
5300 for( e = 0; e < nedges; e++ )
5301 {
5302 cost[e] = graph->cost[e];
5303 costrev[e] = graph->cost[flipedge(e)];
5304
5305 assert(SCIPisGE(scip, cost[e], 0.0));
5306 }
5307
5308 /* compute decomposition of graph and radius values */
5309 graph_voronoiWithRadiusMw(scip, graph, path, cost, radius, vbase, heap, state);
5310
5311 /* sum all radius values, exclude two radius values of lowest value */
5312 for( k = 0; k < nnodes; k++ )
5313 {
5314 if( !Is_term(graph->term[k]) || !graph->mark[k] )
5315 continue;
5316
5317 assert(vbase[k] == k);
5318 assert(SCIPisGE(scip, prize[k], 0.0));
5319
5320 if( SCIPisGE(scip, radius[k], FARAWAY) )
5321 radius[k] = 0.0;
5322 else
5323 {
5324 if( SCIPisGE(scip, radius[k], prize[k] ) )
5325 radius[k] = 0.0;
5326 else
5327 radius[k] = prize[k] - radius[k];
5328 }
5329 }
5330
5331 for( k = 0; k < nnodes; k++ )
5332 {
5333 if( !graph->mark[k] || graph->grad[k] == 0 || Is_term(graph->term[k]) )
5334 continue;
5335 }
5336
5337 /* build Voronoi regions */
5338 graph_voronoiMw(scip, graph, costrev, vnoi, vbase, heap, state);
5339
5340 /* get 2nd next positive node to all non-positive nodes */
5341 graph_get2next(scip, graph, cost, costrev, vnoi, vbase, heap, state);
5342
5343 for( k = 0; k < nnodes; k++ )
5344 {
5345 if( !Is_term(graph->term[k]) || !graph->mark[k] )
5346 continue;
5347
5348 assert(vbase[k] == k);
5349 }
5350
5351 SCIPsortReal(radius, nnodes);
5352
5353 radiim2 = 0.0;
5354
5355 for( k = 2; k < nterms; k++ )
5356 {
5357 assert( SCIPisGT(scip, FARAWAY, radius[k]) );
5358 radiim2 += radius[k];
5359 }
5360
5361 /* solution available? */
5362 if( result != NULL)
5363 {
5364 /* calculate objective value of solution */
5365 obj = 0.0;
5366 for( e = 0; e < nedges; e++ )
5367 {
5368 if( result[e] == CONNECT )
5369 {
5370 head = graph->head[e];
5371
5372 if( graph->mark[head] )
5373 obj += prize[head];
5374 }
5375 }
5376 }
5377 else
5378 {
5379 obj = 0.0;
5380 }
5381
5382 bound = radiim2;
5383
5384 /* traverse all nodes, try to eliminate each non-positive node */
5385 for( k = 0; k < nnodes; k++ )
5386 {
5387 if( !graph->mark[k] || graph->grad[k] == 0 || Is_term(graph->term[k]) )
5388 continue;
5389
5390 assert(SCIPisLE(scip, graph->prize[k], 0.0));
5391
5392 tmpcost = -vnoi[k].dist - vnoi[k + nnodes].dist + bound + graph->prize[k];
5393
5394 if( (SCIPisLT(scip, tmpcost, obj)) )
5395 {
5396 while( graph->outbeg[k] != EAT_LAST )
5397 {
5398 (*nelims)++;
5399 graph_edge_del(scip, graph, graph->outbeg[k], TRUE);
5400 }
5401 }
5402 }
5403
5404 SCIPdebugMessage("nelims (edges) in MWCSP bound reduce: %d,\n", *nelims);
5405
5406 /* free memory*/
5407 SCIPfreeBufferArrayNull(scip, &mst);
5408
5409 return SCIP_OKAY;
5410 }
5411
5412
5413
5414 /** bound-based reductions for the (R)PCSTP, the MWCSP and the STP; used by prune heuristic */
reduce_boundPrune(SCIP * scip,GRAPH * graph,PATH * vnoi,SCIP_Real * cost,SCIP_Real * radius,SCIP_Real * costrev,SCIP_Real * offset,int * heap,int * state,int * vbase,const int * solnode,const int * soledge,int * nelims,int minelims)5415 SCIP_RETCODE reduce_boundPrune(
5416 SCIP* scip, /**< SCIP data structure */
5417 GRAPH* graph, /**< graph data structure */
5418 PATH* vnoi, /**< Voronoi data structure */
5419 SCIP_Real* cost, /**< edge cost array */
5420 SCIP_Real* radius, /**< radius array */
5421 SCIP_Real* costrev, /**< reversed edge cost array */
5422 SCIP_Real* offset, /**< pointer to the offset */
5423 int* heap, /**< heap array */
5424 int* state, /**< array to store state of a node during Voronoi computation */
5425 int* vbase, /**< Voronoi base to each node */
5426 const int* solnode, /**< array of nodes of current solution that is not to be destroyed */
5427 const int* soledge, /**< array of edges of current solution that is not to be destroyed */
5428 int* nelims, /**< pointer to store number of eliminated edges */
5429 int minelims /**< minimum number of edges to be eliminated */
5430 )
5431 {
5432 SCIP_Real* const prize = graph->prize;
5433 SCIP_Real obj;
5434 SCIP_Real max;
5435 SCIP_Real bound;
5436 SCIP_Real tmpcost;
5437 SCIP_Real mstobj;
5438 SCIP_Real maxcost;
5439 SCIP_Real radiim2;
5440 SCIP_Real radiim3;
5441 const int root = graph->source;
5442 const int nnodes = graph->knots;
5443 const int nedges = graph->edges;
5444 const SCIP_Bool mw = (graph->stp_type == STP_MWCSP);
5445 const SCIP_Bool pc = (graph->stp_type == STP_RPCSPG) || (graph->stp_type == STP_PCSPG);
5446 const SCIP_Bool pcmw = (pc || mw);
5447 const int nterms = graph->terms - ( (mw || graph->stp_type == STP_PCSPG) ? 1 : 0 );
5448 SCIP_Bool eliminate;
5449
5450 assert(scip != NULL);
5451 assert(vnoi != NULL);
5452 assert(cost != NULL);
5453 assert(heap != NULL);
5454 assert(graph != NULL);
5455 assert(state != NULL);
5456 assert(vbase != NULL);
5457 assert(nelims != NULL);
5458 assert(radius != NULL);
5459 assert(costrev != NULL);
5460 assert(solnode != NULL);
5461 assert(soledge != NULL);
5462 assert(graph->source >= 0);
5463 assert(graph_valid(graph));
5464
5465 mstobj = 0.0;
5466 *nelims = 0;
5467
5468 if( nterms <= 2 )
5469 return SCIP_OKAY;
5470
5471 assert( graph->stp_type != STP_RPCSPG || Is_term(graph->term[root]) );
5472
5473 if( !pcmw )
5474 for( int k = 0; k < nnodes; k++ )
5475 graph->mark[k] = (graph->grad[k] > 0);
5476
5477 /* initialize cost and costrev array */
5478 maxcost = 0.0;
5479 for( int e = 0; e < nedges; e++ )
5480 {
5481 cost[e] = graph->cost[e];
5482 costrev[e] = graph->cost[flipedge(e)];
5483
5484 assert(SCIPisGE(scip, cost[e], 0.0));
5485
5486 if( graph->stp_type == STP_DHCSTP && SCIPisLT(scip, graph->cost[e], BLOCKED) && SCIPisGT(scip, graph->cost[e], maxcost) )
5487 maxcost = graph->cost[e];
5488 }
5489
5490 /* no MWCSP? */
5491 if( !mw )
5492 {
5493 GRAPH* adjgraph;
5494 PATH* mst;
5495
5496 SCIP_CALL( graph_init(scip, &adjgraph, nterms, MIN(nedges, (nterms - 1) * nterms), 1) );
5497
5498 /* build Voronoi regions, concomitantly building adjgraph and computing radii values*/
5499 SCIP_CALL( graph_voronoiWithRadius(scip, graph, adjgraph, vnoi, radius, cost, costrev, vbase, heap, state) );
5500
5501 /* get 2nd next terminals to all non-terminal nodes */
5502 graph_get2next(scip, graph, cost, costrev, vnoi, vbase, heap, state);
5503
5504 /* get 3th next terminals to all non-terminal nodes */
5505 graph_get3next(scip, graph, cost, costrev, vnoi, vbase, heap, state);
5506
5507 assert(adjgraph != NULL);
5508 graph_knot_chg(adjgraph, 0, 0);
5509 adjgraph->source = 0;
5510 assert(graph_valid(adjgraph));
5511
5512 /* compute MST on adjgraph */
5513 SCIP_CALL( SCIPallocBufferArray(scip, &mst, nterms) );
5514 SCIP_CALL( graph_path_init(scip, adjgraph) );
5515 graph_path_exec(scip, adjgraph, MST_MODE, 0, adjgraph->cost, mst);
5516
5517 max = -1.0;
5518 for( int k = 1; k < nterms; k++ )
5519 {
5520 const int e = mst[k].edge;
5521 assert(adjgraph->path_state[k] == CONNECT);
5522 assert(e >= 0);
5523 tmpcost = adjgraph->cost[e];
5524 mstobj += tmpcost;
5525 if( SCIPisGT(scip, tmpcost, max) )
5526 max = tmpcost;
5527 }
5528 mstobj -= max;
5529
5530 /* free memory*/
5531 SCIPfreeBufferArray(scip, &mst);
5532 graph_path_exit(scip, adjgraph);
5533 graph_free(scip, &adjgraph, TRUE);
5534
5535 }
5536 else
5537 {
5538
5539 /* build Voronoi regions */
5540 graph_voronoiMw(scip, graph, costrev, vnoi, vbase, heap, state);
5541
5542 /* get 2nd next positive node to all non-positive nodes */
5543 graph_get2next(scip, graph, cost, costrev, vnoi, vbase, heap, state);
5544 }
5545
5546
5547 /* for (rooted) prize collecting an maximum weight problems: correct radius values */
5548 if( graph->stp_type == STP_RPCSPG )
5549 {
5550 assert(graph->mark[graph->source]);
5551 for( int k = 0; k < nnodes; k++ )
5552 {
5553 if( !Is_term(graph->term[k]) || !graph->mark[k] )
5554 continue;
5555
5556 if( SCIPisGT(scip, radius[k], prize[k]) && k != graph->source )
5557 radius[k] = prize[k];
5558 }
5559 }
5560 else if( graph->stp_type == STP_PCSPG )
5561 {
5562 for( int k = 0; k < nnodes; k++ )
5563 {
5564 if( !graph->mark[k] )
5565 continue;
5566
5567 if( Is_term(graph->term[k]) && SCIPisGT(scip, radius[k], prize[k]) )
5568 radius[k] = prize[k];
5569 }
5570 }
5571
5572 /* sort radius values */
5573 if( !mw )
5574 SCIPsortReal(radius, nnodes);
5575
5576 radiim2 = 0.0;
5577
5578 if( mw )
5579 {
5580 for( int e = 0; e < nedges; e++ )
5581 costrev[e] = FARAWAY;
5582 bound = 0.0;
5583 radiim3 = 0.0;
5584 }
5585 else
5586 {
5587 for( int e = 0; e < nedges; e++ )
5588 costrev[e] = -1.0;
5589
5590 /* sum all but two radius values of highest/lowest value */
5591 for( int k = 0; k < nterms - 2; k++ )
5592 {
5593 assert( SCIPisGT(scip, FARAWAY, radius[k]) );
5594 radiim2 += radius[k];
5595 }
5596 if( nterms >= 3 )
5597 radiim3 = radiim2 - radius[nterms - 3];
5598 else
5599 radiim3 = 0.0;
5600
5601 if( SCIPisGT(scip, radiim2, mstobj) )
5602 bound = radiim2;
5603 else
5604 bound = mstobj;
5605 }
5606
5607 for( int redrounds = 0; redrounds < 3; redrounds++ )
5608 {
5609 int nrealelims = MIN(2 * minelims, nedges - 1);
5610
5611 if( redrounds == 0 )
5612 {
5613 eliminate = FALSE;
5614 obj = FARAWAY;
5615 }
5616 else if( redrounds == 1 )
5617 {
5618 assert(minelims > 0);
5619 assert(2 * minelims < nedges);
5620 eliminate = TRUE;
5621 SCIPsortReal(costrev, nedges);
5622
5623 if( mw )
5624 {
5625 obj = costrev[nrealelims];
5626 }
5627 else
5628 {
5629 obj = costrev[nedges - nrealelims];
5630
5631 if( SCIPisLT(scip, obj, 0.0) )
5632 obj = 0.0;
5633 }
5634 }
5635 else
5636 {
5637 if( mw )
5638 {
5639 obj = costrev[nrealelims] + 2 * SCIPepsilon(scip);
5640 }
5641 else
5642 {
5643 obj = costrev[nedges - nrealelims] - 2 * SCIPepsilon(scip);
5644
5645 if( SCIPisLT(scip, obj, 0.0) )
5646 obj = 0.0;
5647 }
5648 eliminate = TRUE;
5649 }
5650
5651 if( mw )
5652 {
5653 for( int k = 0; k < nnodes; k++ )
5654 {
5655 if( (*nelims) >= minelims )
5656 break;
5657
5658 if( root == k )
5659 continue;
5660
5661 if( !graph->mark[k] || graph->grad[k] == 0 || Is_gterm(graph->term[k] ) )
5662 continue;
5663
5664 tmpcost = -vnoi[k].dist - vnoi[k + nnodes].dist + graph->prize[k];
5665
5666 if( (!eliminate || SCIPisLT(scip, tmpcost, obj)) && solnode[k] != CONNECT )
5667 {
5668 if( eliminate )
5669 {
5670 while (graph->outbeg[k] != EAT_LAST)
5671 {
5672 (*nelims)++;
5673 graph_edge_del(scip, graph, graph->outbeg[k], TRUE);
5674 }
5675 }
5676 else
5677 {
5678 for( int e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
5679 {
5680 if( SCIPisLT(scip, tmpcost, costrev[e]) )
5681 {
5682 costrev[e] = tmpcost;
5683 costrev[flipedge(e)] = tmpcost;
5684 }
5685 }
5686 }
5687 }
5688 else if( solnode[k] == CONNECT )
5689 {
5690 int e = graph->outbeg[k];
5691
5692 while( e != EAT_LAST )
5693 {
5694 const int etemp = graph->oeat[e];
5695 const int tail = graph->tail[e];
5696 const int head = graph->head[e];
5697
5698 if( tail > head || soledge[e] == CONNECT || soledge[flipedge(e)] == CONNECT )
5699 {
5700 e = etemp;
5701 continue;
5702 }
5703
5704 tmpcost = graph->prize[head] + graph->prize[tail];
5705
5706 if( vbase[tail] != vbase[head] )
5707 {
5708 tmpcost -= vnoi[head].dist + vnoi[tail].dist;
5709 }
5710 else
5711 {
5712 if( SCIPisGT(scip, -vnoi[tail].dist -vnoi[head + nnodes].dist, -vnoi[tail + nnodes].dist -vnoi[head].dist) )
5713 tmpcost -= vnoi[tail].dist + vnoi[head + nnodes].dist;
5714 else
5715 tmpcost -= vnoi[tail + nnodes].dist + vnoi[head].dist;
5716 }
5717 /* can edge e or arc e be deleted? */
5718 if( (!eliminate || SCIPisLT(scip, tmpcost, obj))
5719 && SCIPisLT(scip, graph->cost[e], FARAWAY) && (graph->mark[head]) )
5720 {
5721 assert(!Is_pterm(graph->term[head]));
5722 assert(!Is_pterm(graph->term[tail]));
5723
5724 if( eliminate )
5725 {
5726 graph_edge_del(scip, graph, e, TRUE);
5727 (*nelims)++;
5728 }
5729 else if( SCIPisLT(scip, tmpcost, costrev[e]) )
5730 {
5731 costrev[e] = tmpcost;
5732 costrev[flipedge(e)] = tmpcost;
5733 }
5734
5735 }
5736 e = etemp;
5737 }
5738 }
5739 }
5740 }
5741 /* no MWCSP */
5742 else
5743 {
5744 /* traverse all nodes, try to eliminate each node or incident edges */
5745 for( int k = 0; k < nnodes; k++ )
5746 {
5747 if( (*nelims) >= minelims )
5748 break;
5749 if( root == k )
5750 continue;
5751
5752 if( (!graph->mark[k] && (pcmw)) || graph->grad[k] == 0 )
5753 continue;
5754
5755 tmpcost = vnoi[k].dist + vnoi[k + nnodes].dist + bound;
5756
5757 /* can node k be deleted? */
5758 if( !Is_term(graph->term[k]) && (!eliminate || SCIPisGT(scip, tmpcost, obj)) && solnode[k] != CONNECT )
5759 {
5760 /* delete all incident edges */
5761 if( eliminate )
5762 {
5763 while( graph->outbeg[k] != EAT_LAST )
5764 {
5765 const int e = graph->outbeg[k];
5766 (*nelims)++;
5767
5768 assert(!pc || graph->tail[e] != graph->source);
5769 assert(!pc || graph->mark[graph->head[e]]);
5770 assert(!Is_pterm(graph->term[graph->head[e]]));
5771 assert(!Is_pterm(graph->term[graph->tail[e]]));
5772
5773 graph_edge_del(scip, graph, e, TRUE);
5774 }
5775 }
5776 else
5777 {
5778 for( int e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
5779 {
5780 if( SCIPisGT(scip, tmpcost, costrev[e]) )
5781 {
5782 costrev[e] = tmpcost;
5783 costrev[flipedge(e)] = tmpcost;
5784 }
5785 }
5786 }
5787 }
5788 else
5789 {
5790 int e = graph->outbeg[k];
5791 while( e != EAT_LAST )
5792 {
5793 const int etemp = graph->oeat[e];
5794 const int tail = graph->tail[e];
5795 const int head = graph->head[e];
5796
5797 if( tail > head )
5798 {
5799 e = etemp;
5800 continue;
5801 }
5802
5803 if( soledge[e] == CONNECT || soledge[flipedge(e)] == CONNECT )
5804 {
5805 e = etemp;
5806 continue;
5807 }
5808
5809 tmpcost = graph->cost[e] + bound;
5810
5811 if( vbase[tail] != vbase[head] )
5812 {
5813 tmpcost += vnoi[head].dist + vnoi[tail].dist;
5814 }
5815 else
5816 {
5817 if( SCIPisGT(scip, vnoi[tail].dist + vnoi[head + nnodes].dist, vnoi[tail + nnodes].dist + vnoi[head].dist) )
5818 tmpcost += vnoi[tail + nnodes].dist + vnoi[head].dist;
5819 else
5820 tmpcost += vnoi[tail].dist + vnoi[head + nnodes].dist;
5821
5822 assert(SCIPisGE(scip, tmpcost, vnoi[head].dist + vnoi[tail].dist + graph->cost[e] + bound));
5823 }
5824 /* can edge e or arc e be deleted? */
5825 if( (!eliminate || SCIPisGT(scip, tmpcost, obj))
5826 && SCIPisLT(scip, graph->cost[e], FARAWAY) && (!(pc) || graph->mark[head]) )
5827 {
5828 assert(!Is_pterm(graph->term[head]));
5829 assert(!Is_pterm(graph->term[tail]));
5830
5831 if( eliminate )
5832 {
5833 graph_edge_del(scip, graph, e, TRUE);
5834 (*nelims)++;
5835 }
5836 else if( SCIPisGT(scip, tmpcost, costrev[e]) )
5837 {
5838 costrev[e] = tmpcost;
5839 costrev[flipedge(e)] = tmpcost;
5840 }
5841 }
5842 e = etemp;
5843 }
5844 }
5845 }
5846 #if 1
5847 /* traverse all nodes, try to eliminate 3,4 degree nodes */
5848 for( int k = 0; k < nnodes; k++ )
5849 {
5850 if( (*nelims) >= minelims )
5851 break;
5852
5853 if( (!graph->mark[k] && pc) || graph->grad[k] <= 0 )
5854 continue;
5855
5856 if( solnode[k] == CONNECT )
5857 continue;
5858
5859 if( !eliminate )
5860 {
5861 if( graph->grad[k] >= 3 && graph->grad[k] <= 4 && !Is_term(graph->term[k]) )
5862 {
5863 tmpcost = vnoi[k].dist + vnoi[k + nnodes].dist + vnoi[k + 2 * nnodes].dist + radiim3;
5864 for( int e = graph->outbeg[k]; e != EAT_LAST; e = graph->oeat[e] )
5865 {
5866 if( SCIPisGT(scip, tmpcost, costrev[e]) )
5867 {
5868 costrev[e] = tmpcost;
5869 costrev[flipedge(e)] = tmpcost;
5870 }
5871 }
5872 }
5873 continue;
5874 }
5875
5876 if( graph->grad[k] >= 3 && graph->grad[k] <= 4 && !Is_term(graph->term[k]) )
5877 {
5878 tmpcost = vnoi[k].dist + vnoi[k + nnodes].dist + vnoi[k + 2 * nnodes].dist + radiim3;
5879 if( SCIPisGT(scip, tmpcost, obj) )
5880 {
5881 SCIP_Bool success;
5882
5883 SCIP_CALL( graph_knot_delPseudo(scip, graph, graph->cost, NULL, NULL, k, &success) );
5884
5885 assert(graph->grad[k] == 0 || (graph->grad[k] == 4 && !success));
5886
5887 if( success )
5888 (*nelims)++;
5889 }
5890 }
5891 }
5892 } /* no MWCS */
5893 #endif
5894 } /* redrounds */
5895
5896 SCIPdebugMessage("nelims (edges) in bound reduce: %d,\n", *nelims);
5897
5898 return SCIP_OKAY;
5899 }
5900
5901
5902
5903 /** bound-based reduction test for the HCDSTP */
reduce_boundHop(SCIP * scip,GRAPH * graph,PATH * vnoi,SCIP_Real * cost,SCIP_Real * radius,SCIP_Real * costrev,int * heap,int * state,int * vbase,int * nelims)5904 SCIP_RETCODE reduce_boundHop(
5905 SCIP* scip,
5906 GRAPH* graph,
5907 PATH* vnoi,
5908 SCIP_Real* cost,
5909 SCIP_Real* radius,
5910 SCIP_Real* costrev,
5911 int* heap,
5912 int* state,
5913 int* vbase,
5914 int* nelims
5915 )
5916 {
5917 SCIP_Real max;
5918 SCIP_Real tmpcost;
5919 SCIP_Real bound;
5920 SCIP_Real mstobj;
5921 SCIP_Real radiim2;
5922
5923 GRAPH* adjgraph;
5924 PATH* mst;
5925 int e;
5926 int k;
5927 int tail;
5928 int head;
5929 int etemp;
5930 int nnodes;
5931 int nedges;
5932 int nterms;
5933 SCIP_Real hoplimit;
5934
5935 assert(scip != NULL);
5936 assert(graph != NULL);
5937 assert(vnoi != NULL);
5938 assert(cost != NULL);
5939 assert(radius != NULL);
5940 assert(costrev != NULL);
5941 assert(heap != NULL);
5942 assert(state != NULL);
5943 assert(vbase != NULL);
5944 assert(nelims != NULL);
5945
5946 *nelims = 0;
5947 nterms = 0;
5948 nedges = graph->edges;
5949 nnodes = graph->knots;
5950
5951 for( k = 0; k < nnodes; k++ )
5952 {
5953 graph->mark[k] = (graph->grad[k] > 0);
5954 if( graph->mark[k] && Is_term(graph->term[k]) )
5955 nterms++;
5956 }
5957
5958 for( e = 0; e < nedges; e++ )
5959 {
5960 if( SCIPisLT(scip, graph->cost[e], FARAWAY) )
5961 cost[e] = 1.0;
5962 else
5963 cost[e] = FARAWAY;
5964 if( SCIPisLT(scip, graph->cost[flipedge(e)], FARAWAY) )
5965 costrev[e] = 1.0;
5966 else
5967 costrev[e] = FARAWAY;
5968 }
5969
5970 /* init auxiliary graph */
5971 SCIP_CALL( graph_init(scip, &adjgraph, nterms, MIN(nedges, (nterms - 1) * nterms), 1) );
5972
5973 SCIP_CALL( graph_voronoiWithRadius(scip, graph, adjgraph, vnoi, radius, cost, costrev, vbase, heap, state) );
5974
5975 /* get 2nd next terminals to all nodes */
5976 graph_get2next(scip, graph, cost, costrev, vnoi, vbase, heap, state);
5977
5978 /* compute MST on adjgraph */
5979 graph_knot_chg(adjgraph, 0, 0);
5980 adjgraph->source = 0;
5981 assert(graph_valid(adjgraph));
5982 SCIP_CALL( SCIPallocBufferArray(scip, &mst, nterms) );
5983 SCIP_CALL( graph_path_init(scip, adjgraph) );
5984 graph_path_exec(scip, adjgraph, MST_MODE, 0, adjgraph->cost, mst);
5985
5986 max = -1;
5987 assert(mst[0].edge == -1);
5988 mstobj = 0.0;
5989
5990 /* compute MST cost ...*/
5991 for( k = 1; k < nterms; k++ )
5992 {
5993 e = mst[k].edge;
5994 assert(adjgraph->path_state[k] == CONNECT);
5995 assert(e >= 0);
5996 tmpcost = adjgraph->cost[e];
5997 mstobj += tmpcost;
5998 if( SCIPisGT(scip, tmpcost, max) )
5999 max = tmpcost;
6000 }
6001 /* ...minus longest edge */
6002 mstobj -= max;
6003
6004 SCIPsortReal(radius, nnodes);
6005 radiim2 = 0.0;
6006
6007 for( e = 0; e < nterms - 2; e++ )
6008 {
6009 assert( SCIPisGT(scip, FARAWAY, radius[e]) );
6010 radiim2 += radius[e];
6011 }
6012
6013 hoplimit = (SCIP_Real) graph->hoplimit;
6014
6015 if( SCIPisGT(scip, radiim2, mstobj) )
6016 bound = radiim2;
6017 else
6018 bound = mstobj;
6019
6020 /* traverse all node, try to eliminate first the node and then all incident edges */
6021 for( k = 0; k < nnodes; k++ )
6022 {
6023 /* can node k be deleted? */
6024 if( !Is_term(graph->term[k]) && SCIPisGT(scip, vnoi[k].dist + vnoi[k + nnodes].dist + bound, hoplimit) )
6025 {
6026 e = graph->outbeg[k];
6027
6028 /* delete incident edges */
6029 while( e != EAT_LAST )
6030 {
6031 assert(e >= 0);
6032 (*nelims)++;
6033 etemp = graph->oeat[e];
6034 graph_edge_del(scip, graph, e, TRUE);
6035 e = etemp;
6036 }
6037 }
6038 else
6039 {
6040 e = graph->outbeg[k];
6041 while( e != EAT_LAST )
6042 {
6043 assert(e >= 0);
6044 tail = graph->tail[e];
6045 head = graph->head[e];
6046 tmpcost = 1.0 + bound;
6047 if( vbase[tail] != vbase[head] )
6048 {
6049 tmpcost += vnoi[head].dist + vnoi[tail].dist;
6050 }
6051 else
6052 {
6053 if( SCIPisGT(scip, vnoi[tail].dist + vnoi[head + nnodes].dist, vnoi[tail + nnodes].dist + vnoi[head].dist) )
6054 tmpcost += vnoi[tail + nnodes].dist + vnoi[head].dist;
6055 else
6056 tmpcost += vnoi[tail].dist + vnoi[head + nnodes].dist;
6057 assert(SCIPisGE(scip, tmpcost, vnoi[head].dist + vnoi[tail].dist + 1.0 + bound));
6058 }
6059
6060 /* can edge e (i.e. both arc e and its reverse arc) or arc e be deleted? */
6061 if( SCIPisGT(scip, tmpcost, hoplimit) && SCIPisLT(scip, graph->cost[e], FARAWAY) )
6062 {
6063 etemp = graph->oeat[e];
6064 if( SCIPisGT(scip, graph->cost[e], graph->cost[flipedge(e)]) )
6065 {
6066 graph->cost[e] = FARAWAY;
6067 (*nelims)++;
6068 }
6069 else
6070 {
6071 graph_edge_del(scip, graph, e, TRUE);
6072 (*nelims)++;
6073 }
6074 e = etemp;
6075 }
6076 else
6077 {
6078 e = graph->oeat[e];
6079 }
6080 }
6081 }
6082 }
6083
6084 SCIPdebugMessage("nelimsX (edges) in hop bound reduce: %d,\n", *nelims);
6085
6086 /* free adjgraph */
6087 graph_path_exit(scip, adjgraph);
6088 graph_free(scip, &adjgraph, TRUE);
6089
6090 /* free memory*/
6091 SCIPfreeBufferArray(scip, &mst);
6092 assert(graph_valid(graph));
6093
6094 return SCIP_OKAY;
6095 }
6096
6097 /** hop bound-based reduction test for the HCDSTP */
reduce_boundHopR(SCIP * scip,GRAPH * graph,PATH * vnoi,SCIP_Real * cost,SCIP_Real * costrev,SCIP_Real * pathdist,int * heap,int * state,int * vbase,int * nelims,int * pathedge)6098 SCIP_RETCODE reduce_boundHopR(
6099 SCIP* scip,
6100 GRAPH* graph,
6101 PATH* vnoi,
6102 SCIP_Real* cost,
6103 SCIP_Real* costrev,
6104 SCIP_Real* pathdist,
6105 int* heap,
6106 int* state,
6107 int* vbase,
6108 int* nelims,
6109 int* pathedge
6110 )
6111 {
6112 SCIP_Real tmpcost;
6113 int e;
6114 int k;
6115 int root;
6116 int head;
6117 int etemp;
6118 int bound;
6119 int nnodes;
6120 int nedges;
6121 int nterms;
6122 int hoplimit;
6123
6124 assert(scip != NULL);
6125 assert(graph != NULL);
6126 assert(vnoi != NULL);
6127 assert(cost != NULL);
6128 assert(costrev != NULL);
6129 assert(heap != NULL);
6130 assert(state != NULL);
6131 assert(vbase != NULL);
6132 assert(nelims != NULL);
6133
6134 *nelims = 0;
6135 nterms = 0;
6136 root = graph->source;
6137 nedges = graph->edges;
6138 nnodes = graph->knots;
6139 hoplimit = graph->hoplimit;
6140
6141 for( k = 0; k < nnodes; k++ )
6142 {
6143 graph->mark[k] = (graph->grad[k] > 0);
6144 if( graph->mark[k] && Is_term(graph->term[k]) )
6145 nterms++;
6146 }
6147 bound = nterms - 2;
6148 for( e = 0; e < nedges; e++ )
6149 {
6150 if( SCIPisLT(scip, graph->cost[e], FARAWAY) )
6151 cost[e] = 1.0;
6152 else
6153 cost[e] = graph->cost[e];
6154 if( SCIPisLT(scip, graph->cost[flipedge(e)], FARAWAY) )
6155 costrev[e] = 1.0;
6156 else
6157 costrev[e] = graph->cost[flipedge(e)];
6158 }
6159
6160 /* distance from root to all nodes */
6161 graph_path_execX(scip, graph, root, cost, pathdist, pathedge);
6162
6163 /* no paths should go back to the root */
6164 for( e = graph->outbeg[root]; e != EAT_LAST; e = graph->oeat[e] )
6165 costrev[e] = FARAWAY;
6166
6167 /* build voronoi diagram */
6168 graph_voronoiTerms(scip, graph, costrev, vnoi, vbase, heap, state);
6169
6170 /* traverse all node, try to eliminate first the node and then all incident edges */
6171 for( k = 0; k < nnodes; k++ )
6172 {
6173 /* can node k be deleted? */
6174 if( !Is_term(graph->term[k]) && SCIPisGT(scip, vnoi[k].dist + pathdist[k] + (double) bound, (double) hoplimit) )
6175 {
6176 e = graph->outbeg[k];
6177
6178 /* delete incident edges */
6179 while( e != EAT_LAST )
6180 {
6181 assert(e >= 0);
6182 (*nelims)++;
6183 etemp = graph->oeat[e];
6184 graph_edge_del(scip, graph, e, TRUE);
6185 e = etemp;
6186 }
6187 }
6188 else
6189 {
6190 e = graph->outbeg[k];
6191 while( e != EAT_LAST )
6192 {
6193 assert(e >= 0);
6194 head = graph->head[e];
6195 tmpcost = pathdist[k] + 1.0 + vnoi[head].dist + bound;
6196
6197 etemp = graph->oeat[e];
6198 /* can edge e (i.e. both arc e and its reverse arc) or arc e be deleted? */
6199 if( SCIPisGT(scip, tmpcost, (double) hoplimit) && SCIPisLT(scip, graph->cost[e], FARAWAY) )
6200 {
6201
6202 if( SCIPisGT(scip, FARAWAY, graph->cost[flipedge(e)]) )
6203 {
6204 graph->cost[e] = FARAWAY;
6205 (*nelims)++;
6206 }
6207 else
6208 {
6209 graph_edge_del(scip, graph, e, TRUE);
6210 (*nelims)++;
6211 }
6212 }
6213 e = etemp;
6214 }
6215 }
6216 }
6217
6218 SCIPdebugMessage("eliminated (edges) in hcr bound reduce: %d,\n", *nelims);
6219
6220 assert(graph_valid(graph));
6221
6222 return SCIP_OKAY;
6223 }
6224
6225 /* reduction method for HCSTP */
reduce_boundHopRc(SCIP * scip,GRAPH * graph,PATH * vnoi,SCIP_Real * cost,SCIP_Real * costrev,SCIP_Real * pathdist,SCIP_Real objval,int * heap,int * state,int * vbase,int * nelims,int * pathedge,SCIP_Bool fix)6226 SCIP_RETCODE reduce_boundHopRc(
6227 SCIP* scip,
6228 GRAPH* graph,
6229 PATH* vnoi,
6230 SCIP_Real* cost,
6231 SCIP_Real* costrev,
6232 SCIP_Real* pathdist,
6233 SCIP_Real objval,
6234 int* heap,
6235 int* state,
6236 int* vbase,
6237 int* nelims,
6238 int* pathedge,
6239 SCIP_Bool fix
6240 )
6241 {
6242 SCIP_VAR** vars;
6243 SCIP_HEURDATA* tmheurdata;
6244 SCIP_Real min;
6245 SCIP_Real bound;
6246 SCIP_Real maxmin;
6247 SCIP_Real maxcost;
6248 SCIP_Real tmpcost;
6249 SCIP_Real hopfactor;
6250 int* result;
6251 int e;
6252 int k;
6253 int root;
6254 int head;
6255 int etemp;
6256 int nnodes;
6257 int nedges;
6258 int best_start;
6259 SCIP_Bool success;
6260
6261 assert(scip != NULL);
6262 assert(graph != NULL);
6263 assert(vnoi != NULL);
6264 assert(cost != NULL);
6265 assert(costrev != NULL);
6266 assert(heap != NULL);
6267 assert(state != NULL);
6268 assert(vbase != NULL);
6269 assert(nelims != NULL);
6270
6271 hopfactor = DEFAULT_HOPFACTOR;
6272 bound = 0.0;
6273 *nelims = 0;
6274 best_start = 0;
6275 success = TRUE;
6276 vars = NULL;
6277 root = graph->source;
6278 nedges = graph->edges;
6279 nnodes = graph->knots;
6280
6281 SCIP_CALL( SCIPallocBufferArray(scip, &result, nedges) );
6282 maxcost = 0.0;
6283 if( fix )
6284 {
6285 vars = SCIPprobdataGetVars(scip);
6286 assert(vars != NULL);
6287 for( e = 0; e < nedges; e += 2 )
6288 {
6289 result[e] = UNKNOWN;
6290 result[e + 1] = UNKNOWN;
6291
6292 if( SCIPvarGetUbLocal(vars[e + 1]) < 0.5 )
6293 {
6294 costrev[e] = BLOCKED;
6295 }
6296 else
6297 {
6298 costrev[e] = graph->cost[e + 1];
6299
6300 if( SCIPisGT(scip, costrev[e], maxcost) && SCIPisLT(scip, costrev[e], BLOCKED) )
6301 maxcost = costrev[e];
6302 }
6303 cost[e + 1] = costrev[e];
6304 if( SCIPvarGetUbLocal(vars[e]) < 0.5 )
6305 {
6306 costrev[e + 1] = BLOCKED;
6307 }
6308 else
6309 {
6310 costrev[e + 1] = graph->cost[e];
6311
6312 if( SCIPisGT(scip, graph->cost[e], maxcost) && SCIPisLT(scip, costrev[e + 1], BLOCKED) )
6313 maxcost = graph->cost[e];
6314 }
6315 cost[e] = costrev[e + 1];
6316 }
6317 }
6318 else
6319 {
6320 for( e = 0; e < nedges; e++ )
6321 {
6322 result[e] = UNKNOWN;
6323 cost[e] = graph->cost[e];
6324 costrev[e] = graph->cost[flipedge(e)];
6325 if( SCIPisGT(scip, graph->cost[e], maxcost) )
6326 maxcost = graph->cost[e];
6327 }
6328 }
6329
6330 maxmin = -1.0;
6331 for( k = 0; k < nnodes; k++ )
6332 {
6333 graph->mark[k] = (graph->grad[k] > 0);
6334 if( graph->mark[k] && Is_term(graph->term[k]) )
6335 {
6336 if( k != root )
6337 {
6338 min = FARAWAY;
6339 for( e = graph->inpbeg[k]; e != EAT_LAST; e = graph->ieat[e] )
6340 if( SCIPisLT(scip, cost[e], min) )
6341 min = cost[e];
6342 assert(SCIPisGT(scip, BLOCKED, min));
6343 if( SCIPisGT(scip, min, maxmin) )
6344 maxmin = min;
6345 bound += min;
6346 }
6347 }
6348 }
6349 bound -= maxmin;
6350
6351
6352 /* distance from root to all nodes */
6353 graph_path_execX(scip, graph, root, cost, pathdist, pathedge);
6354
6355 /* no paths should go back to the root */
6356 for( e = graph->outbeg[root]; e != EAT_LAST; e = graph->oeat[e] )
6357 costrev[e] = FARAWAY;
6358
6359 /* build voronoi diagram */
6360 graph_voronoiTerms(scip, graph, costrev, vnoi, vbase, heap, state);
6361
6362 if( SCIPisLT(scip, objval, 0.0) )
6363 {
6364 /* get TM heuristic data */
6365 assert(SCIPfindHeur(scip, "TM") != NULL);
6366 tmheurdata = SCIPheurGetData(SCIPfindHeur(scip, "TM"));
6367
6368 /* compute UB */
6369 SCIP_CALL( SCIPStpHeurTMRun(scip, tmheurdata, graph, NULL, &best_start, result, 50, root, cost, costrev, &hopfactor, NULL, maxcost, &success, FALSE) );
6370
6371 objval = 0.0;
6372 for( e = 0; e < nedges; e++ )
6373 if( result[e] == CONNECT )
6374 objval += graph->cost[e];
6375 }
6376 else
6377 {
6378 /* objval = objval - fixed; */
6379 objval = SCIPgetCutoffbound(scip);
6380 assert(SCIPisGT(scip, objval, 0.0));
6381 }
6382
6383 /* traverse all node, try to eliminate first the node and then all incident edges */
6384 for( k = 0; k < nnodes; k++ )
6385 {
6386 if( Is_term(graph->term[k]) )
6387 continue;
6388 /* can node k be deleted? */
6389 if( SCIPisGT(scip, vnoi[k].dist + pathdist[k] + bound, objval) )
6390 {
6391 e = graph->outbeg[k];
6392
6393 /* delete incident edges */
6394 while( e != EAT_LAST )
6395 {
6396 assert(e >= 0);
6397
6398 etemp = graph->oeat[e];
6399 if( fix )
6400 {
6401 assert(vars != NULL);
6402 /* try to fix edge */
6403 SCIP_CALL( fixedgevar(scip, vars[e], nelims) );
6404
6405 /* try to fix reversed edge */
6406 SCIP_CALL( fixedgevar(scip, vars[flipedge(e)], nelims) );
6407 }
6408 else
6409 {
6410 graph_edge_del(scip, graph, e, TRUE);
6411 (*nelims)++;
6412 }
6413 e = etemp;
6414 }
6415 }
6416 else
6417 {
6418 e = graph->outbeg[k];
6419 while( e != EAT_LAST )
6420 {
6421 assert(e >= 0);
6422 head = graph->head[e];
6423 tmpcost = pathdist[k] + graph->cost[e] + vnoi[head].dist + bound;
6424
6425 etemp = graph->oeat[e];
6426 /* can edge e (i.e. both arc e and its reverse arc) or arc e be deleted? */
6427 if( SCIPisGT(scip, tmpcost, objval) && SCIPisLT(scip, graph->cost[e], FARAWAY) )
6428 {
6429 if( fix )
6430 {
6431 assert(vars != NULL);
6432
6433 /* try to fix edge */
6434 SCIP_CALL( fixedgevar(scip, vars[e], nelims) );
6435 }
6436 else
6437 {
6438 if( SCIPisGT(scip, FARAWAY, graph->cost[flipedge(e)]) )
6439 {
6440 graph->cost[e] = FARAWAY;
6441 (*nelims)++;
6442 }
6443 else
6444 {
6445 graph_edge_del(scip, graph, e, TRUE);
6446 (*nelims)++;
6447 }
6448 }
6449 }
6450 e = etemp;
6451 }
6452 }
6453 }
6454
6455 SCIPdebugMessage("CCC eliminated (edges) in hcrc bound reduce: %d,\n", *nelims);
6456 /* free memory */
6457 SCIPfreeBufferArray(scip, &result);
6458
6459 assert(graph_valid(graph));
6460
6461 return SCIP_OKAY;
6462 }
6463