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 grphbase.c
17 * @brief includes several methods for Steiner problem graphs
18 * @author Thorsten Koch
19 * @author Daniel Rehfeldt
20 *
21 * This file contains several basic methods to process Steiner problem graphs.
22 * A graph can not be reduced in terms of edge or node size, but edges can be marked as
23 * EAT_FREE (to not be used anymore) and nodes may have degree one.
24 * The method 'graph_pack()' can be used to build a new graph that discards these nodes and edges.
25 *
26 */
27
28 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
29
30 /*lint -esym(766,stdlib.h) -esym(766,malloc.h) */
31 /*lint -esym(766,string.h) */
32
33 #include "scip/misc.h"
34 #include <stdio.h>
35 #include <stdlib.h>
36 #include <string.h>
37 #include <assert.h>
38 #include "portab.h"
39 #include "misc_stp.h"
40 #include "scip/misc.h"
41 #include "grph.h"
42 #include "heur_tm.h"
43
44 #define STP_DELPSEUDO_MAXGRAD 5
45 #define STP_DELPSEUDO_MAXNEDGES 10
46
47 /*
48 * local functions
49 */
50
51 /** can edge in pseudo-elimination method be cut off? */
52 inline static
cutoffEdge(SCIP * scip,const SCIP_Real * cutoffs,const SCIP_Real * cutoffsrev,const SCIP_Real * ecost,const SCIP_Real * ecostrev,int edgeidx1,int edgeidx2,int cutoffidx)53 SCIP_Bool cutoffEdge(
54 SCIP* scip, /**< SCIP data */
55 const SCIP_Real* cutoffs, /**< cutoff values for each incident edge */
56 const SCIP_Real* cutoffsrev, /**< revere cutoff values (or NULL if undirected) */
57 const SCIP_Real* ecost, /**< edge cost*/
58 const SCIP_Real* ecostrev, /**< reverse edge cost */
59 int edgeidx1, /**< index of first edge to be checked (wrt provided arrays) */
60 int edgeidx2, /**< index of second edge to be checked (wrt provided arrays) */
61 int cutoffidx /**< index for cutoff array */
62 )
63 {
64 SCIP_Real newcost;
65
66 assert(edgeidx1 != edgeidx2);
67
68 if( cutoffs == NULL )
69 return FALSE;
70
71 newcost = ecostrev[edgeidx1] + ecost[edgeidx2];
72
73 if( !SCIPisGT(scip, newcost, cutoffs[cutoffidx]) )
74 return FALSE;
75
76 if( cutoffsrev != NULL )
77 {
78 const SCIP_Real newcostrev = ecost[edgeidx1] + ecostrev[edgeidx2];
79
80 if( !SCIPisGT(scip, newcostrev, cutoffsrev[cutoffidx]) )
81 return FALSE;
82 }
83
84 return TRUE;
85 }
86
87
88 inline static
removeEdge(GRAPH * g,int e)89 void removeEdge(
90 GRAPH* g, /**< the graph */
91 int e /**< the edge to be removed */
92 )
93 {
94 int i;
95 int head;
96 int tail;
97
98 assert(g != NULL);
99 assert(e >= 0);
100 assert(e < g->edges);
101
102 head = g->head[e];
103 tail = g->tail[e];
104
105 if( g->inpbeg[head] == e )
106 g->inpbeg[head] = g->ieat[e];
107 else
108 {
109 if( g->rootedgeprevs != NULL && head == g->source )
110 {
111 i = g->rootedgeprevs[e];
112 assert(g->ieat[i] == e);
113 if( g->ieat[e] >= 0 )
114 g->rootedgeprevs[g->ieat[e]] = i;
115 }
116 else
117 for( i = g->inpbeg[head]; g->ieat[i] != e; i = g->ieat[i] )
118 assert(i >= 0);
119
120 g->ieat[i] = g->ieat[e];
121 }
122 if( g->outbeg[tail] == e )
123 g->outbeg[tail] = g->oeat[e];
124 else
125 {
126 if( g->rootedgeprevs != NULL && tail == g->source )
127 {
128 i = g->rootedgeprevs[e];
129 assert(g->oeat[i] == e);
130 if( g->oeat[e] >= 0 )
131 g->rootedgeprevs[g->oeat[e]] = i;
132 }
133 else
134 for( i = g->outbeg[tail]; g->oeat[i] != e; i = g->oeat[i] )
135 assert(i >= 0);
136
137 g->oeat[i] = g->oeat[e];
138 }
139 }
140
141 /** used by graph_grid_create */
142 static
getNodeNumber(int grid_dim,int shiftcoord,int * ncoords,int * currcoord)143 int getNodeNumber(
144 int grid_dim,
145 int shiftcoord,
146 int* ncoords,
147 int* currcoord
148 )
149 {
150 int number = 0;
151 int tmp;
152 int i;
153 int j;
154 for( i = 0; i < grid_dim; i++ )
155 {
156 tmp = 1;
157 for( j = i + 1; j < grid_dim; j++ )
158 {
159 tmp = tmp * ncoords[j];
160 }
161 if( shiftcoord == i )
162 number += (currcoord[i] + 1) * tmp;
163 else
164 number += currcoord[i] * tmp;
165 }
166 number++;
167 return number;
168 }
169
170 /** used by graph_obstgrid_create */
171 static
compEdgesObst(int coord,int grid_dim,int nobstacles,int * ncoords,int * currcoord,int * edgecosts,int * gridedgecount,int ** coords,int ** gridedges,int ** obst_coords,char * inobstacle)172 void compEdgesObst(
173 int coord,
174 int grid_dim,
175 int nobstacles,
176 int* ncoords,
177 int* currcoord,
178 int* edgecosts,
179 int* gridedgecount,
180 int** coords,
181 int** gridedges,
182 int** obst_coords,
183 char* inobstacle
184 )
185 {
186 char inobst;
187 int i;
188 int j;
189 int z;
190 int x;
191 int y;
192 int node;
193 i = 0;
194 while( i < ncoords[coord] )
195 {
196 currcoord[coord] = i;
197 if( coord < grid_dim - 1 )
198 compEdgesObst(coord + 1, grid_dim, nobstacles, ncoords, currcoord, edgecosts, gridedgecount, coords, gridedges, obst_coords, inobstacle);
199 else
200 {
201 x = coords[0][currcoord[0]];
202 y = coords[1][currcoord[1]];
203 inobst = FALSE;
204 node = getNodeNumber(grid_dim, -1, ncoords, currcoord);
205 for( z = 0; z < nobstacles; z++ )
206 {
207 assert(obst_coords[0][z] < obst_coords[2][z]);
208 assert(obst_coords[1][z] < obst_coords[3][z]);
209 if( x > obst_coords[0][z] && x < obst_coords[2][z] &&
210 y > obst_coords[1][z] && y < obst_coords[3][z] )
211 {
212 inobst = TRUE;
213 inobstacle[node-1] = TRUE;
214 break;
215 }
216 }
217 for( j = 0; j < grid_dim; j++ )
218 {
219 if( currcoord[j] + 1 < ncoords[j] )
220 {
221 if( inobst == FALSE )
222 {
223 gridedges[0][*gridedgecount] = node;
224 gridedges[1][*gridedgecount] = getNodeNumber(grid_dim, j, ncoords, currcoord);
225 edgecosts[*gridedgecount] = coords[j][currcoord[j] + 1] - coords[j][currcoord[j]];
226 (*gridedgecount)++;
227 }
228 }
229 }
230 }
231 i++;
232 }
233 }
234
235 /** used by graph_grid_create */
236 static
compEdges(int coord,int grid_dim,int * ncoords,int * currcoord,int * edgecosts,int * gridedgecount,int ** coords,int ** gridedges)237 void compEdges(
238 int coord,
239 int grid_dim,
240 int* ncoords,
241 int* currcoord,
242 int* edgecosts,
243 int* gridedgecount,
244 int** coords,
245 int** gridedges
246 )
247 {
248 int j;
249 int i = 0;
250 while( i < ncoords[coord] )
251 {
252 currcoord[coord] = i;
253 if( coord < grid_dim - 1 )
254 compEdges(coord + 1, grid_dim, ncoords, currcoord, edgecosts, gridedgecount, coords, gridedges);
255 else
256 {
257 for( j = 0; j < grid_dim; j++ )
258 {
259 if( currcoord[j] + 1 < ncoords[j] )
260 {
261 gridedges[0][*gridedgecount] = getNodeNumber(grid_dim, -1, ncoords, currcoord);
262 gridedges[1][*gridedgecount] = getNodeNumber(grid_dim, j, ncoords, currcoord);
263 edgecosts[*gridedgecount] = coords[j][currcoord[j] + 1] - coords[j][currcoord[j]];
264 (*gridedgecount)++;
265 }
266 }
267 }
268 i++;
269 }
270 }
271
272 /*
273 * global functions
274 */
275
276
277 #if 0
278 /** transforms an MWCSP to an SAP */
279 SCIP_RETCODE graph_MwcsToSap(
280 SCIP* scip, /**< SCIP data structure */
281 GRAPH* graph, /**< the graph */
282 SCIP_Real* maxweights /**< array containing the weight of each node */
283 )
284 {
285 int e;
286 int i;
287 int nnodes;
288 int nterms = 0;
289
290 assert(maxweights != NULL);
291 assert(scip != NULL);
292 assert(graph != NULL);
293 assert(graph->cost != NULL);
294 assert(graph->terms == 0);
295
296 nnodes = graph->knots;
297
298 /* count number of terminals, modify incoming edges for non-terminals */
299 for( i = 0; i < nnodes; i++ )
300 {
301 if( SCIPisLT(scip, maxweights[i], 0.0) )
302 {
303 for( e = graph->inpbeg[i]; e != EAT_LAST; e = graph->ieat[e] )
304 {
305 graph->cost[e] -= maxweights[i];
306 }
307 }
308 else
309 {
310 graph_knot_chg(graph, i, 0);
311 nterms++;
312 }
313 }
314 nterms = 0;
315 for( i = 0; i < nnodes; i++ )
316 {
317 graph->prize[i] = maxweights[i];
318 if( Is_term(graph->term[i]) )
319 {
320 assert(SCIPisGE(scip, maxweights[i], 0.0));
321 nterms++;
322 }
323 else
324 {
325 assert(SCIPisLT(scip, maxweights[i], 0.0));
326 }
327 }
328 assert(nterms == graph->terms);
329 graph->stp_type = STP_MWCSP;
330
331 SCIP_CALL( graph_PcToSap(scip, graph) );
332 assert(graph->stp_type == STP_MWCSP);
333 return SCIP_OKAY;
334 }
335
336
337 /** alters the graph for prize collecting problems */
338 SCIP_RETCODE graph_PcToSap(
339 SCIP* scip, /**< SCIP data structure */
340 GRAPH* graph /**< the graph */
341 )
342 {
343 SCIP_Real* prize;
344 int k;
345 int root;
346 int node;
347 int nnodes;
348 int nterms;
349 int pseudoroot;
350
351 assert(graph != NULL);
352 assert(graph->prize != NULL);
353 assert(graph->knots == graph->ksize);
354 assert(graph->edges == graph->esize);
355
356 prize = graph->prize;
357 nnodes = graph->knots;
358 nterms = graph->terms;
359 graph->norgmodelknots = nnodes;
360 graph->norgmodeledges = graph->edges;
361
362 /* for each terminal, except for the root, one node and three edges (i.e. six arcs) are to be added */
363 SCIP_CALL( graph_resize(scip, graph, (graph->ksize + graph->terms + 2), (graph->esize + graph->terms * 8) , -1) );
364
365 /* create a new nodes */
366 for( k = 0; k < nterms; ++k )
367 graph_knot_add(graph, -1);
368
369 /* new pseudo-root */
370 pseudoroot = graph->knots;
371 graph_knot_add(graph, -1);
372
373 /* new root */
374 root = graph->knots;
375 graph_knot_add(graph, 0);
376
377 nterms = 0;
378 for( k = 0; k < nnodes; ++k )
379 {
380 /* is the kth node a terminal other than the root? */
381 if( Is_term(graph->term[k]) )
382 {
383 /* the copied node */
384 node = nnodes + nterms;
385 nterms++;
386
387 /* switch the terminal property, mark k */
388 graph_knot_chg(graph, k, -2);
389 graph_knot_chg(graph, node, 0);
390 assert(SCIPisGE(scip, prize[k], 0.0));
391
392 /* add one edge going from the root to the 'copied' terminal and one going from the former terminal to its copy */
393 graph_edge_add(scip, graph, root, k, BLOCKED, FARAWAY);
394 graph_edge_add(scip, graph, pseudoroot, node, prize[k], FARAWAY);
395 graph_edge_add(scip, graph, k, node, 0.0, FARAWAY);
396 graph_edge_add(scip, graph, k, pseudoroot, 0.0, FARAWAY);
397 }
398 else if( graph->stp_type != STP_MWCSP )
399 {
400 prize[k] = 0;
401 }
402 }
403 graph->source = root;
404 graph->extended = TRUE;
405 assert((nterms + 1) == graph->terms);
406 if( graph->stp_type != STP_MWCSP )
407 graph->stp_type = STP_PCSPG;
408
409 return SCIP_OKAY;
410 }
411
412
413
414
415 #endif
416
417
418 /** creates a graph out of a given grid */
graph_obstgrid_create(SCIP * scip,GRAPH ** gridgraph,int ** coords,int ** obst_coords,int nterms,int grid_dim,int nobstacles,int scale_order)419 SCIP_RETCODE graph_obstgrid_create(
420 SCIP* scip, /**< SCIP data structure */
421 GRAPH** gridgraph, /**< the (obstacle) grid graph to be constructed */
422 int** coords, /**< coordinates of all points */
423 int** obst_coords, /**< coordinates of obstacles */
424 int nterms, /**< number of terminals */
425 int grid_dim, /**< dimension of the problem */
426 int nobstacles, /**< number of obstacles*/
427 int scale_order /**< scale factor */
428 )
429 {
430 GRAPH* g;
431 GRAPH* gp;
432 double cost;
433 int i;
434 int j;
435 int k;
436 int tmp;
437 int shift;
438 int nnodes;
439 int nedges;
440 double scale_factor;
441 int gridedgecount;
442 int* ncoords;
443 int* currcoord;
444 int* edgecosts;
445 int** termcoords;
446 int** gridedges;
447 char* inobstacle;
448 assert(coords != NULL);
449 assert(grid_dim > 1);
450 assert(nterms > 0);
451 assert(grid_dim == 2);
452 scale_factor = pow(10.0, (double) scale_order);
453
454 /* initialize the terminal-coordinates array */
455 SCIP_CALL( SCIPallocBufferArray(scip, &termcoords, grid_dim) );
456
457 for( i = 0; i < grid_dim; i++ )
458 {
459 SCIP_CALL( SCIPallocBufferArray(scip, &(termcoords[i]), nterms) ); /*lint !e866*/
460 for( j = 0; j < nterms; j++ )
461 termcoords[i][j] = coords[i][j];
462 }
463
464 SCIP_CALL( SCIPallocBufferArray(scip, &ncoords, grid_dim) );
465 SCIP_CALL( SCIPallocBufferArray(scip, &currcoord, grid_dim) );
466
467 /* sort the coordinates and delete multiples */
468 for( i = 0; i < grid_dim; i++ )
469 {
470 ncoords[i] = 1;
471 SCIPsortInt(coords[i], nterms);
472 shift = 0;
473 for( j = 0; j < nterms - 1; j++ )
474 {
475 if( coords[i][j] == coords[i][j + 1] )
476 {
477 shift++;
478 }
479 else
480 {
481 coords[i][j + 1 - shift] = coords[i][j + 1];
482 ncoords[i]++;
483 }
484 }
485 }
486
487 nnodes = 1;
488
489 for( i = 0; i < grid_dim; i++ )
490 nnodes = nnodes * ncoords[i];
491
492 tmp = 0;
493
494 for( i = 0; i < grid_dim; i++ )
495 tmp = tmp + nnodes / ncoords[i];
496
497 nedges = grid_dim * nnodes - tmp;
498 SCIP_CALL( SCIPallocBufferArray(scip, &gridedges, 2) );
499 SCIP_CALL( SCIPallocBufferArray(scip, &edgecosts, nedges) );
500 SCIP_CALL( SCIPallocBufferArray(scip, &(gridedges[0]), nedges) );
501 SCIP_CALL( SCIPallocBufferArray(scip, &(gridedges[1]), nedges) );
502 SCIP_CALL( SCIPallocBufferArray(scip, &(inobstacle), nnodes) );
503 gridedgecount = 0;
504 for( i = 0; i < nnodes; i++ )
505 inobstacle[i] = FALSE;
506 compEdgesObst(0, grid_dim, nobstacles, ncoords, currcoord, edgecosts, &gridedgecount, coords, gridedges, obst_coords, inobstacle);
507 nedges = gridedgecount;
508 /* initialize empty g with allocated slots for nodes and edges */
509 SCIP_CALL( graph_init(scip, gridgraph, nnodes, 2 * nedges, 1) );
510
511 g = *gridgraph;
512 SCIP_CALL( SCIPallocMemoryArray(scip, &(g->grid_ncoords), grid_dim) );
513 for( i = 0; i < grid_dim; i++ )
514 g->grid_ncoords[i] = ncoords[i];
515
516 g->grid_dim = grid_dim;
517 g->grid_coordinates = coords;
518
519 /* add nodes */
520 for( i = 0; i < nnodes; i++ )
521 graph_knot_add(g, -1);
522
523 /* add edges */
524 for( i = 0; i < nedges; i++ )
525 {
526 /* (re) scale edge costs */
527 if( inobstacle[gridedges[1][i] - 1] == FALSE )
528 {
529 cost = ((double) edgecosts[i]) / scale_factor;
530 graph_edge_add(scip, g, gridedges[0][i] - 1, gridedges[1][i] - 1, cost, cost);
531 }
532 }
533
534 /* add terminals */
535 for( i = 0; i < nterms; i++ )
536 {
537 for( j = 0; j < grid_dim; j++ )
538 {
539 for( k = 0; k <= ncoords[j]; k++ )
540 {
541 assert(k != ncoords[j]);
542 if( coords[j][k] == termcoords[j][i] )
543 {
544 currcoord[j] = k;
545 break;
546 }
547 }
548 }
549 /* the position of the (future) terminal */
550 k = getNodeNumber(grid_dim, -1, ncoords, currcoord) - 1;
551
552 if( i == 0 )
553 g->source = k;
554
555 /* make a terminal out of the node */
556 graph_knot_chg(g, k, 0);
557 }
558
559 SCIP_CALL( graph_pack(scip, g, &gp, TRUE) );
560 g = gp;
561 g->stp_type = STP_OARSMT;
562
563 SCIPfreeBufferArray(scip, &inobstacle);
564 SCIPfreeBufferArray(scip, &(gridedges[1]));
565 SCIPfreeBufferArray(scip, &(gridedges[0]));
566 SCIPfreeBufferArray(scip, &edgecosts);
567 SCIPfreeBufferArray(scip, &gridedges);
568 SCIPfreeBufferArray(scip, &currcoord);
569 SCIPfreeBufferArray(scip, &ncoords);
570
571 for( i = grid_dim - 1; i >= 0 ; --i )
572 SCIPfreeBufferArray(scip, &(termcoords[i]));
573
574 SCIPfreeBufferArray(scip, &termcoords);
575
576 return SCIP_OKAY;
577 }
578
579
580
581 /** creates a graph out of a given grid */
graph_grid_create(SCIP * scip,GRAPH ** gridgraph,int ** coords,int nterms,int grid_dim,int scale_order)582 SCIP_RETCODE graph_grid_create(
583 SCIP* scip, /**< SCIP data structure */
584 GRAPH** gridgraph, /**< the grid graph to be constructed */
585 int** coords, /**< coordinates */
586 int nterms, /**< number of terminals*/
587 int grid_dim, /**< problem dimension */
588 int scale_order /**< scale order */
589 )
590 {
591 GRAPH* g;
592 double cost;
593 int i;
594 int j;
595 int k;
596 int tmp;
597 int shift;
598 int nnodes;
599 int nedges;
600 double scale_factor;
601 int gridedgecount;
602 int* ncoords;
603 int* currcoord;
604 int* edgecosts;
605 int** termcoords;
606 int** gridedges;
607 assert(coords != NULL);
608 assert(grid_dim > 1);
609 assert(nterms > 0);
610
611 scale_factor = pow(10.0, (double) scale_order);
612
613 /* initialize the terminal-coordinates array */
614 SCIP_CALL( SCIPallocBufferArray(scip, &termcoords, grid_dim) );
615 for( i = 0; i < grid_dim; i++ )
616 {
617 SCIP_CALL( SCIPallocBufferArray(scip, &(termcoords[i]), nterms) ); /*lint !e866*/
618 for( j = 0; j < nterms; j++ )
619 termcoords[i][j] = coords[i][j];
620 }
621 SCIP_CALL( SCIPallocBufferArray(scip, &ncoords, grid_dim) );
622 SCIP_CALL( SCIPallocBufferArray(scip, &currcoord, grid_dim) );
623
624 /* sort the coordinates and delete multiples */
625 for( i = 0; i < grid_dim; i++ )
626 {
627 ncoords[i] = 1;
628 SCIPsortInt(coords[i], nterms);
629 shift = 0;
630 for( j = 0; j < nterms - 1; j++ )
631 {
632 if( coords[i][j] == coords[i][j + 1] )
633 {
634 shift++;
635 }
636 else
637 {
638 coords[i][j + 1 - shift] = coords[i][j + 1];
639 ncoords[i]++;
640 }
641 }
642 }
643
644 nnodes = 1;
645 for( i = 0; i < grid_dim; i++ )
646 nnodes = nnodes * ncoords[i];
647
648 tmp = 0;
649 for( i = 0; i < grid_dim; i++ )
650 {
651 tmp = tmp + nnodes / ncoords[i];
652 }
653
654 nedges = grid_dim * nnodes - tmp;
655
656 SCIP_CALL( SCIPallocBufferArray(scip, &gridedges, 2) );
657 SCIP_CALL( SCIPallocBufferArray(scip, &edgecosts, nedges) );
658 SCIP_CALL( SCIPallocBufferArray(scip, &(gridedges[0]), nedges) );
659 SCIP_CALL( SCIPallocBufferArray(scip, &(gridedges[1]), nedges) );
660
661 gridedgecount = 0;
662
663 compEdges(0, grid_dim, ncoords, currcoord, edgecosts, &gridedgecount, coords, gridedges);
664
665 /* initialize empty graph with allocated slots for nodes and edges */
666 SCIP_CALL( graph_init(scip, gridgraph, nnodes, 2 * nedges, 1) );
667
668 g = *gridgraph;
669
670 SCIP_CALL( SCIPallocMemoryArray(scip, &(g->grid_ncoords), grid_dim) );
671 for( i = 0; i < grid_dim; i++ )
672 g->grid_ncoords[i] = ncoords[i];
673
674 g->grid_dim = grid_dim;
675 g->grid_coordinates = coords;
676
677 /* add nodes */
678 for( i = 0; i < nnodes; i++ )
679 graph_knot_add(g, -1);
680
681 /* add edges */
682 for( i = 0; i < nedges; i++ )
683 {
684 /* (re) scale edge costs */
685 cost = (double) edgecosts[i] / scale_factor;
686 graph_edge_add(scip, g, gridedges[0][i] - 1, gridedges[1][i] - 1, cost, cost);
687 }
688
689 /* add terminals */
690 for( i = 0; i < nterms; i++ )
691 {
692 for( j = 0; j < grid_dim; j++ )
693 {
694 for( k = 0; k <= ncoords[j]; k++ )
695 {
696 assert(k != ncoords[j]);
697 if( coords[j][k] == termcoords[j][i] )
698 {
699 currcoord[j] = k;
700 break;
701 }
702 }
703 }
704 /* the position of the (future) terminal */
705 k = getNodeNumber(grid_dim, -1, ncoords, currcoord) - 1;
706
707 /* make a terminal out of the node */
708 graph_knot_chg(g, k, 0);
709 }
710
711 g->stp_type = STP_RSMT;
712
713 SCIPfreeBufferArray(scip, &(gridedges[1]));
714 SCIPfreeBufferArray(scip, &(gridedges[0]));
715 SCIPfreeBufferArray(scip, &edgecosts);
716 SCIPfreeBufferArray(scip, &gridedges);
717 SCIPfreeBufferArray(scip, &currcoord);
718 SCIPfreeBufferArray(scip, &ncoords);
719
720 for( i = grid_dim - 1; i >= 0 ; i-- )
721 SCIPfreeBufferArray(scip, &(termcoords[i]));
722
723 SCIPfreeBufferArray(scip, &termcoords);
724
725 return SCIP_OKAY;
726 }
727
728
729 /** computes coordinates of node 'node' */
graph_grid_coordinates(SCIP * scip,int ** coords,int ** nodecoords,int * ncoords,int node,int grid_dim)730 SCIP_RETCODE graph_grid_coordinates(
731 SCIP* scip, /**< SCIP data structure */
732 int** coords, /**< coordinates */
733 int** nodecoords, /**< coordinates of the node (to be computed) */
734 int* ncoords, /**< array with number of coordinate */
735 int node, /**< the node */
736 int grid_dim /**< problem dimension */
737 )
738 {
739 int i;
740 int j;
741 int tmp;
742 int coord;
743 assert(grid_dim > 1);
744 assert(node >= 0);
745 assert(coords != NULL);
746 assert(ncoords != NULL);
747 if( *nodecoords == NULL )
748 SCIP_CALL( SCIPallocMemoryArray(scip, nodecoords, grid_dim) );
749
750 for( i = 0; i < grid_dim; i++ )
751 {
752 tmp = 1;
753 for( j = i; j < grid_dim; j++ )
754 tmp = tmp * ncoords[j];
755
756 coord = node % tmp;
757 tmp = tmp / ncoords[i];
758 coord = coord / tmp;
759 (*nodecoords)[i] = coords[i][coord];
760 }
761 return SCIP_OKAY;
762 }
763
764
765 /** allocates (first and second) and initializes (only second) arrays for PC and MW problems */
graph_pc_init(SCIP * scip,GRAPH * g,int sizeprize,int sizeterm2edge)766 SCIP_RETCODE graph_pc_init(
767 SCIP* scip, /**< SCIP data structure */
768 GRAPH* g, /**< the graph */
769 int sizeprize, /**< size of prize array to allocate (or -1) */
770 int sizeterm2edge /**< size of term2edge array to allocate and initialize to -1 (or -1) */
771 )
772 {
773 assert(scip != NULL);
774 assert(g != NULL);
775
776 if( sizeprize > 0 )
777 {
778 assert(NULL == g->prize);
779 SCIP_CALL( SCIPallocMemoryArray(scip, &(g->prize), sizeprize) );
780 }
781
782 if( sizeterm2edge > 0 )
783 {
784 assert(NULL == g->term2edge);
785 SCIP_CALL( SCIPallocMemoryArray(scip, &(g->term2edge), sizeterm2edge) );
786 for( int i = 0; i < sizeterm2edge; i++ )
787 g->term2edge[i] = -1;
788 }
789
790 return SCIP_OKAY;
791 }
792
793 /** changes graph of PC and MW problems needed for presolving routines */
graph_pc_presolInit(SCIP * scip,GRAPH * g)794 SCIP_RETCODE graph_pc_presolInit(
795 SCIP* scip, /**< SCIP data structure */
796 GRAPH* g /**< the graph */
797 )
798 {
799 int prev;
800 const int root = g->source;
801 const int nedges = g->edges;
802
803 if( g->stp_type == STP_RPCSPG )
804 return SCIP_OKAY;
805
806 assert(scip != NULL && g != NULL);
807 assert(g->rootedgeprevs == NULL);
808 assert(nedges > 0 && g->grad[root] > 0);
809
810 SCIP_CALL( SCIPallocMemoryArray(scip, &(g->rootedgeprevs), nedges) );
811
812 for( int e = 0; e < nedges; e++ )
813 g->rootedgeprevs[e] = -1;
814
815 prev = g->outbeg[root];
816 assert(prev != EAT_LAST);
817
818 for( int e = g->oeat[prev]; e != EAT_LAST; e = g->oeat[e] )
819 {
820 g->rootedgeprevs[e] = prev;
821 prev = e;
822 }
823
824 prev = g->inpbeg[root];
825 assert(prev != EAT_LAST);
826
827 for( int e = g->ieat[prev]; e != EAT_LAST; e = g->ieat[e] )
828 {
829 g->rootedgeprevs[e] = prev;
830 prev = e;
831 }
832
833 return SCIP_OKAY;
834 }
835
836 /** changes graph of PC and MW problems needed after exiting presolving routines */
graph_pc_presolExit(SCIP * scip,GRAPH * g)837 void graph_pc_presolExit(
838 SCIP* scip, /**< SCIP data structure */
839 GRAPH* g /**< the graph */
840 )
841 {
842 assert(scip != NULL && g != NULL);
843
844 if( g->stp_type == STP_RPCSPG )
845 return;
846
847 assert(g->rootedgeprevs != NULL);
848
849 SCIPfreeMemoryArray(scip, &(g->rootedgeprevs));
850 }
851
852 /** checks consistency of term2edge array ONLY for non-extended graphs! */
graph_pc_term2edgeConsistent(const GRAPH * g)853 SCIP_Bool graph_pc_term2edgeConsistent(
854 const GRAPH* g /**< the graph */
855 )
856 {
857 assert(g != NULL);
858 assert(g->term2edge);
859 assert(!g->extended);
860
861 if( g->term2edge[g->source] != -1 )
862 return FALSE;
863
864 for( int i = 0; i < g->knots; i++ )
865 {
866 if( Is_gterm(g->term[i]) && i != g->source && g->term2edge[i] < 0 )
867 {
868 SCIPdebugMessage("term2edge consistency fail1 %d \n", i);
869 return FALSE;
870 }
871
872 if( !Is_gterm(g->term[i]) && g->term2edge[i] != -1 )
873 {
874 SCIPdebugMessage("term2edge consistency fail2 %d \n", i);
875 return FALSE;
876 }
877
878 if( Is_pterm(g->term[i]) && i != g->source )
879 {
880 int k = -1;
881 int e;
882
883 for( e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e] )
884 {
885 k = g->head[e];
886 if( Is_term(g->term[k]) && k != g->source )
887 break;
888 }
889 assert(e != EAT_LAST);
890 assert(k >= 0);
891
892 if( g->term2edge[i] != e )
893 {
894 SCIPdebugMessage("term2edge consistency fail3 %d \n", i);
895 return FALSE;
896 }
897
898 if( g->term2edge[k] != flipedge(e) )
899 {
900 SCIPdebugMessage("term2edge consistency fail4 %d \n", i);
901 return FALSE;
902 }
903 }
904 }
905 return TRUE;
906 }
907
908 /** change property of node to non-terminal */
graph_pc_knot2nonTerm(GRAPH * g,int node)909 void graph_pc_knot2nonTerm(
910 GRAPH* g, /**< the graph */
911 int node /**< node to be changed */
912 )
913 {
914 assert(g != NULL);
915 assert(node >= 0);
916 assert(node < g->knots);
917 assert(g->stp_type == STP_PCSPG || g->stp_type == STP_RPCSPG || g->stp_type == STP_MWCSP || g->stp_type == STP_RMWCSP);
918 assert(g->term2edge);
919
920 if( Is_term(g->term[node]) )
921 g->terms--;
922
923 g->term[node] = -1;
924 g->term2edge[node] = -1;
925 }
926
927 /** updates term2edge array for new graph */
graph_pc_updateTerm2edge(GRAPH * newgraph,const GRAPH * oldgraph,int newtail,int newhead,int oldtail,int oldhead)928 void graph_pc_updateTerm2edge(
929 GRAPH* newgraph, /**< the new graph */
930 const GRAPH* oldgraph, /**< the old graph */
931 int newtail, /**< tail in new graph */
932 int newhead, /**< head in new graph */
933 int oldtail, /**< tail in old graph */
934 int oldhead /**< head in old graph */
935 )
936 {
937 assert(newgraph != NULL);
938 assert(oldgraph != NULL);
939 assert(newgraph->term2edge != NULL);
940 assert(oldgraph->term2edge != NULL);
941 assert(newtail >= 0);
942 assert(newhead >= 0);
943 assert(oldtail >= 0);
944 assert(oldhead >= 0);
945 assert(oldgraph->extended);
946 assert(newgraph->extended);
947
948 assert(newgraph->term2edge != NULL);
949 if( oldgraph->term2edge[oldtail] >= 0 && oldgraph->term2edge[oldhead] >= 0 && oldgraph->term[oldtail] != oldgraph->term[oldhead] )
950 {
951 assert(Is_gterm(newgraph->term[newtail]) && Is_gterm(newgraph->term[newhead]));
952 assert(Is_gterm(oldgraph->term[oldtail]) && Is_gterm(oldgraph->term[oldhead]));
953 assert(oldgraph->source != oldtail && oldgraph->source != oldhead);
954 assert(flipedge(newgraph->edges) == newgraph->edges + 1);
955
956 newgraph->term2edge[newtail] = newgraph->edges;
957 newgraph->term2edge[newhead] = newgraph->edges + 1;
958 }
959
960 assert(-1 == newgraph->term2edge[newgraph->source]);
961 }
962
963 /** mark terminals and switch terminal property to original terminals */
graph_pc_2org(GRAPH * graph)964 void graph_pc_2org(
965 GRAPH* graph /**< the graph */
966 )
967 {
968 int root;
969 int nnodes;
970
971 assert(graph != NULL);
972 assert(graph->extended);
973
974 root = graph->source;
975 nnodes = graph->knots;
976
977 for( int k = 0; k < nnodes; k++ )
978 {
979 graph->mark[k] = (graph->grad[k] > 0);
980
981 if( Is_pterm(graph->term[k]) )
982 {
983 graph_knot_chg(graph, k, 0);
984 }
985 else if( Is_term(graph->term[k]) )
986 {
987 graph->mark[k] = FALSE;
988 if( k != root )
989 graph_knot_chg(graph, k, -2);
990 }
991 }
992
993 if( graph->stp_type == STP_RPCSPG || graph->stp_type == STP_RMWCSP )
994 graph->mark[root] = TRUE;
995
996 graph->extended = FALSE;
997
998 return;
999 }
1000
1001 /** unmark terminals and switch terminal property to transformed terminals */
graph_pc_2trans(GRAPH * graph)1002 void graph_pc_2trans(
1003 GRAPH* graph /**< the graph */
1004 )
1005 {
1006 const int root = graph->source;
1007 const int nnodes = graph->knots;;
1008
1009 assert(graph != NULL);
1010 assert(!(graph->extended));
1011
1012 for( int k = 0; k < nnodes; k++ )
1013 {
1014 graph->mark[k] = (graph->grad[k] > 0);
1015
1016 if( Is_pterm(graph->term[k]) )
1017 graph_knot_chg(graph, k, 0);
1018 else if( Is_term(graph->term[k]) && k != root )
1019 graph_knot_chg(graph, k, -2);
1020 }
1021
1022 graph->extended = TRUE;
1023
1024 return;
1025 }
1026
1027 /** graph_pc_2org if extended */
graph_pc_2orgcheck(GRAPH * graph)1028 void graph_pc_2orgcheck(
1029 GRAPH* graph /**< the graph */
1030 )
1031 {
1032 assert(graph != NULL);
1033
1034 if( !graph->extended )
1035 return;
1036
1037 graph_pc_2org(graph);
1038 }
1039
1040 /** graph_pc_2trans if not extended */
graph_pc_2transcheck(GRAPH * graph)1041 void graph_pc_2transcheck(
1042 GRAPH* graph /**< the graph */
1043 )
1044 {
1045 assert(graph != NULL);
1046
1047 if( graph->extended )
1048 return;
1049
1050 graph_pc_2trans(graph);
1051 }
1052
1053 /* returns sum of positive vertex weights */
graph_pc_getPosPrizeSum(SCIP * scip,const GRAPH * graph)1054 SCIP_Real graph_pc_getPosPrizeSum(
1055 SCIP* scip, /**< SCIP data structure */
1056 const GRAPH* graph /**< the graph */
1057 )
1058 {
1059 SCIP_Real prizesum = 0.0;
1060
1061 assert(scip != NULL);
1062 assert(graph != NULL);
1063 assert(graph->prize != NULL);
1064 assert(!graph->extended);
1065
1066 for( int i = 0; i < graph->knots; i++ )
1067 if( Is_term(graph->term[i]) && i != graph->source && graph->prize[i] < BLOCKED )
1068 prizesum += graph->prize[i];
1069
1070 return prizesum;
1071 }
1072
1073
1074 /** alters the graph for prize collecting problems */
graph_pc_getSap(SCIP * scip,GRAPH * graph,GRAPH ** newgraph,SCIP_Real * offset)1075 SCIP_RETCODE graph_pc_getSap(
1076 SCIP* scip, /**< SCIP data structure */
1077 GRAPH* graph, /**< the graph */
1078 GRAPH** newgraph, /**< the new graph */
1079 SCIP_Real* offset /**< offset */
1080 )
1081 {
1082 SCIP_Real* prize;
1083 SCIP_Real max;
1084 SCIP_Real prizesum;
1085 int k;
1086 int e;
1087 int enext;
1088 int root;
1089 int head;
1090 int nnodes;
1091 int nterms;
1092 int stp_type;
1093 int pseudoroot;
1094
1095 assert(scip != NULL);
1096 assert(graph != NULL);
1097 assert(graph->prize != NULL);
1098 assert(graph->knots == graph->ksize);
1099 assert(graph->edges == graph->esize);
1100
1101 prize = graph->prize;
1102 nnodes = graph->knots;
1103 nterms = graph->terms;
1104 prizesum = 0.0;
1105
1106 stp_type = graph->stp_type;
1107 graph->stp_type = STP_SAP;
1108
1109 /* for each terminal, except for the root, three edges (i.e. six arcs) are to be added */
1110 SCIP_CALL( graph_copy(scip, graph, newgraph) );
1111
1112 graph->stp_type = stp_type;
1113
1114 SCIP_CALL( graph_resize(scip, (*newgraph), ((*newgraph)->ksize + 1), ((*newgraph)->esize + 2 * (nterms - 1)) , -1) );
1115
1116 (*newgraph)->source = graph->source;
1117 root = (*newgraph)->source;
1118
1119 /* new pseudo-root */
1120 pseudoroot = (*newgraph)->knots;
1121 graph_knot_add((*newgraph), -1);
1122
1123 max = 0.0;
1124 for( k = 0; k < nnodes; k++ )
1125 if( Is_pterm(graph->term[k]) )
1126 {
1127 prizesum += prize[k];
1128
1129 if( prize[k] > max )
1130 max = prize[k];
1131 }
1132
1133 prizesum -= max;
1134 *offset -= prizesum;
1135
1136 SCIP_CALL( graph_pc_presolInit(scip, *newgraph) );
1137
1138 e = (*newgraph)->outbeg[root];
1139
1140 while( e != EAT_LAST )
1141 {
1142 enext = (*newgraph)->oeat[e];
1143 head = (*newgraph)->head[e];
1144 if( Is_term((*newgraph)->term[head]) )
1145 {
1146 (void) graph_edge_redirect(scip, (*newgraph), e, pseudoroot, head, graph->cost[e], TRUE);
1147 (*newgraph)->cost[flipedge(e)] = FARAWAY;
1148 assert((*newgraph)->head[e] == head);
1149 assert((*newgraph)->tail[e] == pseudoroot);
1150 }
1151 else
1152 {
1153 (*newgraph)->cost[e] = prizesum;
1154 }
1155
1156 e = enext;
1157 }
1158
1159 graph_pc_presolExit(scip, *newgraph);
1160
1161 for( k = 0; k < nnodes; k++ )
1162 {
1163 /* is the kth node a terminal other than the root? */
1164 if( Is_pterm((*newgraph)->term[k]) )
1165 {
1166 graph_edge_add(scip, (*newgraph), k, pseudoroot, 0.0, FARAWAY);
1167 }
1168 }
1169
1170 return SCIP_OKAY;
1171 }
1172
1173 /** adapts SAP deriving from PCST or MWCS problem with new big M */
graph_pc_adaptSap(SCIP * scip,SCIP_Real bigM,GRAPH * graph,SCIP_Real * offset)1174 void graph_pc_adaptSap(
1175 SCIP* scip, /**< SCIP data structure */
1176 SCIP_Real bigM, /**< new big M value */
1177 GRAPH* graph, /**< the SAP graph */
1178 SCIP_Real* offset /**< the offset */
1179 )
1180 {
1181 SCIP_Real oldbigM;
1182 const int root = graph->source;
1183
1184 assert(bigM > 0.0);
1185 assert(scip != NULL && graph != NULL && offset != NULL);
1186 assert(graph->outbeg[root] >= 0);
1187
1188 oldbigM = graph->cost[graph->outbeg[root]];
1189 assert(oldbigM > 0.0);
1190
1191 *offset += (oldbigM - bigM);
1192
1193 printf("new vs old %f, %f \n", bigM, oldbigM);
1194
1195 for( int e = graph->outbeg[root]; e != EAT_LAST; e = graph->oeat[e] )
1196 {
1197 assert(graph->cost[e] == oldbigM);
1198 graph->cost[e] = bigM;
1199 }
1200 }
1201
1202
1203 /** alters the graph for prize collecting problems and shifts weights to reduce number of terminal */
graph_pc_getSapShift(SCIP * scip,GRAPH * graph,GRAPH ** newgraph,SCIP_Real * offset)1204 SCIP_RETCODE graph_pc_getSapShift(
1205 SCIP* scip, /**< SCIP data structure */
1206 GRAPH* graph, /**< the graph */
1207 GRAPH** newgraph, /**< the new graph */
1208 SCIP_Real* offset /**< offset */
1209 )
1210 {
1211 GRAPH* newg;
1212 SCIP_Real maxp;
1213 SCIP_Real* const prize = graph->prize;
1214 SCIP_Real prizesum;
1215 int e;
1216 int root;
1217 const int nnodes = graph->knots;
1218 const int stp_type = graph->stp_type;
1219 int maxvert;
1220 int pseudoroot;
1221
1222 assert(scip != NULL);
1223 assert(graph != NULL);
1224 assert(graph->prize != NULL);
1225 assert(graph->knots == graph->ksize);
1226 assert(graph->edges == graph->esize);
1227 assert(stp_type == STP_MWCSP || stp_type == STP_PCSPG);
1228 assert(graph->extended);
1229 graph->stp_type = STP_SAP;
1230
1231 /* for each terminal, except for the root, three edges (i.e. six arcs) are to be added */
1232 SCIP_CALL( graph_copy(scip, graph, newgraph) );
1233
1234 graph->stp_type = stp_type;
1235 newg = *newgraph;
1236
1237 /* get max prize and max vertex */
1238 maxvert = -1;
1239 maxp = -1.0;
1240 for( int k = 0; k < nnodes; k++ )
1241 if( Is_pterm(graph->term[k]) && SCIPisGT(scip, graph->prize[k], maxp) )
1242 {
1243 assert(graph->grad[k] > 0);
1244 maxp = graph->prize[k];
1245 maxvert = k;
1246 }
1247
1248 assert(maxvert >= 0);
1249
1250 /* shift the costs */
1251 for( int k = 0; k < nnodes; k++ )
1252 {
1253 newg->mark[k] = (newg->grad[k] > 0);
1254 if( Is_pterm(graph->term[k]) && k != maxvert )
1255 {
1256 SCIP_Real p;
1257
1258 assert(newg->mark[k]);
1259
1260 p = prize[k];
1261 for( e = graph->inpbeg[k]; e != EAT_LAST; e = graph->ieat[e] )
1262 if( SCIPisLT(scip, graph->cost[e], p) && !Is_term(graph->term[graph->tail[e]]) )
1263 break;
1264
1265 /* if there is no incoming arc of lower cost than prize[k], make k a common node */
1266 if( e == EAT_LAST )
1267 {
1268 int e2 = -1;
1269 int term = -1;
1270
1271 newg->term[k] = -1;
1272
1273 for( e = graph->inpbeg[k]; e != EAT_LAST; e = graph->ieat[e] )
1274 {
1275 const int tail = graph->tail[e];
1276
1277 if( Is_term(graph->term[tail]) )
1278 {
1279 if( tail == graph->source )
1280 e2 = e;
1281 else
1282 {
1283 assert(term == -1);
1284 term = tail;
1285 }
1286 }
1287 else
1288 {
1289 newg->cost[e] -= p;
1290 assert(SCIPisGE(scip, newg->cost[e], 0.0));
1291 }
1292 }
1293 prize[k] = 0.0;
1294
1295 (*offset) += p;
1296 assert(e2 != -1);
1297 assert(term != -1);
1298
1299 while( newg->inpbeg[term] != EAT_LAST )
1300 graph_edge_del(scip, newg, newg->inpbeg[term], FALSE);
1301
1302 newg->mark[term] = FALSE;
1303 graph_knot_chg(newg, k, -1);
1304 graph_knot_chg(newg, term, -1);
1305 graph_edge_del(scip, newg, e2, FALSE);
1306 }
1307 }
1308 }
1309
1310 SCIP_CALL( graph_resize(scip, newg, (newg->ksize + 1), (newg->esize + 2 * (newg->terms - 1)) , -1) );
1311
1312 assert(newg->source == graph->source);
1313 root = newg->source;
1314
1315 /* new pseudo-root */
1316 pseudoroot = newg->knots;
1317 graph_knot_add(newg, -1);
1318
1319 prizesum = 0.0;
1320 for( int k = 0; k < nnodes; k++ )
1321 if( Is_pterm(graph->term[k]) )
1322 prizesum += prize[k];
1323
1324 prizesum += 1;
1325
1326 *offset -= prizesum;
1327
1328 /* move edges to terminal from root to pseudo-root */
1329 e = newg->outbeg[root];
1330 while( e != EAT_LAST )
1331 {
1332 const int head = newg->head[e];
1333 const int enext = newg->oeat[e];
1334
1335 if( Is_term(newg->term[head]) )
1336 {
1337 (void) graph_edge_redirect(scip, newg, e, pseudoroot, head, graph->cost[e], TRUE);
1338 newg->cost[flipedge(e)] = FARAWAY;
1339 assert(newg->head[e] == head);
1340 assert(newg->tail[e] == pseudoroot);
1341 }
1342 else
1343 {
1344 newg->cost[e] = prizesum;
1345 }
1346
1347 e = enext;
1348 }
1349
1350 /* add edges from pterminals to pseudo-root */
1351 for( int k = 0; k < nnodes; k++ )
1352 {
1353 /* is the kth node a terminal other than the root? */
1354 if( Is_pterm(newg->term[k]) )
1355 {
1356 assert(newg->mark[k]);
1357 graph_edge_add(scip, newg, k, pseudoroot, 0.0, FARAWAY);
1358 }
1359 }
1360
1361 return SCIP_OKAY;
1362 }
1363
1364 /** alters the graph for prize-collecting problems with given root */
graph_pc_getRsap(SCIP * scip,GRAPH * graph,GRAPH ** newgraph,int * rootcands,int nrootcands,int root)1365 SCIP_RETCODE graph_pc_getRsap(
1366 SCIP* scip, /**< SCIP data structure */
1367 GRAPH* graph, /**< the graph */
1368 GRAPH** newgraph, /**< the new graph */
1369 int* rootcands, /**< array containing all vertices that could be used as root */
1370 int nrootcands, /**< number of all vertices could be used as root */
1371 int root /**< the root of the new SAP */
1372 )
1373 {
1374 GRAPH* p;
1375 int k;
1376 int e;
1377 int e2;
1378 int head;
1379 int aterm;
1380 int proot;
1381 int nnodes;
1382 int stp_type;
1383
1384 assert(graph != NULL);
1385 assert(graph->prize != NULL);
1386 assert(graph->knots == graph->ksize);
1387 assert(graph->edges == graph->esize);
1388
1389 graph_pc_2transcheck(graph);
1390
1391 aterm = -1;
1392 proot = graph->source;
1393 stp_type = graph->stp_type;
1394 graph->stp_type = STP_SAP;
1395
1396 /* copy graph */
1397 SCIP_CALL( graph_copy(scip, graph, newgraph) );
1398
1399 p = *newgraph;
1400
1401 graph->stp_type = stp_type;
1402
1403 assert(Is_pterm(graph->term[root]));
1404
1405 for( e = p->outbeg[root]; e != EAT_LAST; e = p->oeat[e] )
1406 {
1407 head = p->head[e];
1408 if( Is_term(p->term[head]) && head != proot )
1409 {
1410 graph_knot_chg(p, head, -1);
1411 aterm = head;
1412 graph_edge_del(scip, p, e, FALSE);
1413 break;
1414 }
1415 }
1416
1417 assert(aterm != -1);
1418 SCIP_CALL( graph_pc_presolInit(scip, p) );
1419
1420 for( e = graph->outbeg[proot]; e != EAT_LAST; e = graph->oeat[e] )
1421 {
1422 head = graph->head[e];
1423
1424 assert(graph->head[e] == p->head[e]);
1425 assert(graph->tail[e] == p->tail[e]);
1426
1427 if( Is_term(graph->term[head]) && head != aterm )
1428 {
1429 assert(Is_term(p->term[head]));
1430
1431 (void) graph_edge_redirect(scip, p, e, root, head, graph->cost[e], TRUE);
1432 p->cost[flipedge(e)] = FARAWAY;
1433
1434 for( e2 = p->outbeg[head]; e2 != EAT_LAST; e2 = p->oeat[e2] )
1435 if( p->head[e2] != root )
1436 assert(p->term[p->head[e2]] == -2);
1437 }
1438 else
1439 {
1440 graph_edge_del(scip, p, e, FALSE);
1441 }
1442 }
1443
1444 graph_pc_presolExit(scip, p);
1445
1446 assert(p->grad[aterm] == 0);
1447
1448 nnodes = p->knots;
1449 p->source = root;
1450 graph_knot_chg(p, root, 0);
1451
1452 for( k = 0; k < nnodes; k++ )
1453 p->mark[k] = (p->grad[k] > 0);
1454
1455 assert(p->grad[graph->source] == 0);
1456
1457 SCIP_CALL( graph_pc_init(scip, p, nnodes, nnodes) );
1458
1459 assert(graph->term2edge != NULL);
1460
1461 for( k = 0; k < nnodes; k++)
1462 {
1463 p->term2edge[k] = graph->term2edge[k];
1464 if( k < graph->norgmodelknots )
1465 p->prize[k] = graph->prize[k];
1466 else
1467 p->prize[k] = 0.0;
1468 }
1469 p->term2edge[root] = -1;
1470 p->term2edge[aterm] = -1;
1471
1472 if( nrootcands > 0 )
1473 {
1474 SCIP_CALL( graph_pc_presolInit(scip, p) );
1475 for( k = 0; k < nrootcands; k++ )
1476 {
1477 aterm = rootcands[k];
1478 if( aterm == root )
1479 continue;
1480
1481 head = -1;
1482
1483 for( e = p->outbeg[aterm]; e != EAT_LAST; e = p->oeat[e] )
1484 {
1485 head = p->head[e];
1486
1487 if( Is_term(p->term[head]) && p->term2edge[head] >= 0 )
1488 {
1489 assert(p->grad[head] == 2);
1490 assert(head != root);
1491
1492 while( p->outbeg[head] != EAT_LAST )
1493 graph_edge_del(scip, p, p->outbeg[head], FALSE);
1494
1495 graph_knot_chg(p, head, -1);
1496 break;
1497 }
1498 }
1499
1500 assert(head >= 0);
1501
1502 p->term2edge[head] = -1;
1503 p->term2edge[aterm] = -1;
1504
1505 assert(e != EAT_LAST);
1506 graph_knot_chg(p, aterm, 0);
1507 }
1508 graph_pc_presolExit(scip, p);
1509 }
1510
1511 graph_knot_chg(p, proot, -1);
1512 p->prize[root] = 0.0;
1513
1514 return SCIP_OKAY;
1515 }
1516
1517
1518
1519 /** alters the graph for prize collecting problems */
graph_pc_2pc(SCIP * scip,GRAPH * graph)1520 SCIP_RETCODE graph_pc_2pc(
1521 SCIP* scip, /**< SCIP data structure */
1522 GRAPH* graph /**< the graph */
1523 )
1524 {
1525 SCIP_Real* prize;
1526 int root;
1527 const int nnodes = graph->knots;
1528 int nterms;
1529
1530 assert(graph != NULL);
1531 assert(graph->edges == graph->esize);
1532
1533 nterms = graph->terms;
1534 prize = graph->prize;
1535 assert(prize != NULL);
1536 assert(nnodes == graph->ksize);
1537 graph->norgmodeledges = graph->edges;
1538 graph->norgmodelknots = nnodes;
1539
1540 /* for each terminal, except for the root, one node and three edges (i.e. six arcs) are to be added */
1541 SCIP_CALL( graph_resize(scip, graph, (graph->ksize + graph->terms + 1), (graph->esize + graph->terms * 6) , -1) );
1542
1543 /* add new nodes */
1544 for( int k = 0; k < nterms; ++k )
1545 graph_knot_add(graph, -1);
1546
1547 /* new root */
1548 root = graph->knots;
1549 graph_knot_add(graph, 0);
1550
1551 /* allocate and initialize term2edge array */
1552 graph_pc_init(scip, graph, -1, graph->knots);
1553 assert(NULL != graph->term2edge);
1554
1555 nterms = 0;
1556 for( int k = 0; k < nnodes; ++k )
1557 {
1558 /* is the kth node a terminal other than the root? */
1559 if( Is_term(graph->term[k]) )
1560 {
1561 /* the copied node */
1562 const int node = nnodes + nterms;
1563 nterms++;
1564
1565 /* switch the terminal property, mark k */
1566 graph_knot_chg(graph, k, -2);
1567 graph_knot_chg(graph, node, 0);
1568 assert(SCIPisGE(scip, prize[k], 0.0));
1569
1570 /* add one edge going from the root to the 'copied' terminal and one going from the former terminal to its copy */
1571 graph_edge_add(scip, graph, root, k, 0.0, FARAWAY);
1572 graph_edge_add(scip, graph, root, node, prize[k], FARAWAY);
1573
1574 graph->term2edge[k] = graph->edges;
1575 graph->term2edge[node] = graph->edges + 1;
1576 assert(graph->edges + 1 == flipedge(graph->edges));
1577
1578 graph_edge_add(scip, graph, k, node, 0.0, FARAWAY);
1579
1580 assert(graph->head[graph->term2edge[k]] == node);
1581 assert(graph->head[graph->term2edge[node]] == k);
1582 }
1583 else if( graph->stp_type != STP_MWCSP )
1584 {
1585 prize[k] = 0;
1586 }
1587 }
1588 graph->source = root;
1589 graph->extended = TRUE;
1590 assert((nterms + 1) == graph->terms);
1591 if( graph->stp_type != STP_MWCSP )
1592 graph->stp_type = STP_PCSPG;
1593
1594 SCIPdebugMessage("Transformed to PC \n");
1595
1596 return SCIP_OKAY;
1597 }
1598
1599
1600 /** alters the graph for rooted prize collecting problems */
graph_pc_2rpc(SCIP * scip,GRAPH * graph)1601 SCIP_RETCODE graph_pc_2rpc(
1602 SCIP* scip, /**< SCIP data structure */
1603 GRAPH* graph /**< the graph */
1604 )
1605 {
1606 SCIP_Real* prize;
1607 int root;
1608 int node;
1609 int nnodes;
1610 int nterms;
1611
1612 assert(graph != NULL);
1613 assert(graph->edges == graph->esize);
1614
1615 root = graph->source;
1616 nnodes = graph->knots;
1617 nterms = graph->terms;
1618 prize = graph->prize;
1619 graph->norgmodeledges = graph->edges;
1620 graph->norgmodelknots = nnodes;
1621
1622 assert(prize != NULL);
1623 assert(nnodes == graph->ksize);
1624 assert(root >= 0);
1625
1626 /* for each terminal, except for the root, one node and three edges (i.e. six arcs) are to be added */
1627 SCIP_CALL( graph_resize(scip, graph, (graph->ksize + graph->terms), (graph->esize + graph->terms * 4) , -1) );
1628
1629 /* create a new nodes */
1630 for( int k = 0; k < nterms - 1; ++k )
1631 graph_knot_add(graph, -1);
1632
1633 /* allocate and initialize term2edge array */
1634 graph_pc_init(scip, graph, -1, graph->knots);
1635 assert(graph->term2edge != NULL);
1636
1637 nterms = 0;
1638
1639 for( int k = 0; k < nnodes; ++k )
1640 {
1641 /* is the kth node a terminal other than the root? */
1642 if( Is_term(graph->term[k]) && k != root )
1643 {
1644 /* the copied node */
1645 node = nnodes + nterms;
1646 nterms++;
1647 /* switch the terminal property, mark k as former terminal */
1648 graph_knot_chg(graph, k, -2);
1649 graph_knot_chg(graph, node, 0);
1650 assert(SCIPisGE(scip, prize[k], 0.0));
1651
1652 /* add one edge going from the root to the 'copied' terminal and one going from the former terminal to its copy */
1653 graph_edge_add(scip, graph, root, node, prize[k], FARAWAY);
1654
1655 graph->term2edge[k] = graph->edges;
1656 graph->term2edge[node] = graph->edges + 1;
1657 assert(graph->edges + 1 == flipedge(graph->edges));
1658
1659 graph_edge_add(scip, graph, k, node, 0.0, FARAWAY);
1660
1661 assert(graph->head[graph->term2edge[k]] == node);
1662 assert(graph->head[graph->term2edge[node]] == k);
1663 }
1664 else
1665 {
1666 prize[k] = 0.0;
1667 }
1668 }
1669 /* one for the root */
1670 nterms++;
1671
1672 graph->extended = TRUE;
1673 assert(nterms == graph->terms);
1674 graph->stp_type = STP_RPCSPG;
1675
1676 SCIPdebugMessage("Transformed to RPC \n");
1677
1678 return SCIP_OKAY;
1679 }
1680
1681 /** alters the graph for MWCS problems */
graph_pc_2mw(SCIP * scip,GRAPH * graph,SCIP_Real * maxweights)1682 SCIP_RETCODE graph_pc_2mw(
1683 SCIP* scip, /**< SCIP data structure */
1684 GRAPH* graph, /**< the graph */
1685 SCIP_Real* maxweights /**< array containing the weight of each node */
1686 )
1687 {
1688 int nnodes;
1689 int nterms = 0;
1690
1691 assert(maxweights != NULL);
1692 assert(scip != NULL);
1693 assert(graph != NULL);
1694 assert(graph->cost != NULL);
1695 assert(graph->terms == 0);
1696
1697 nnodes = graph->knots;
1698
1699 /* count number of terminals, modify incoming edges for non-terminals */
1700 for( int i = 0; i < nnodes; i++ )
1701 {
1702 if( SCIPisLT(scip, maxweights[i], 0.0) )
1703 {
1704 for( int e = graph->inpbeg[i]; e != EAT_LAST; e = graph->ieat[e] )
1705 graph->cost[e] -= maxweights[i];
1706 }
1707 else if( SCIPisGT(scip, maxweights[i], 0.0) )
1708 {
1709 graph_knot_chg(graph, i, 0);
1710 nterms++;
1711 }
1712 }
1713 nterms = 0;
1714 for( int i = 0; i < nnodes; i++ )
1715 {
1716 graph->prize[i] = maxweights[i];
1717 if( Is_term(graph->term[i]) )
1718 {
1719 assert(SCIPisGE(scip, maxweights[i], 0.0));
1720 nterms++;
1721 }
1722 else
1723 {
1724 assert(SCIPisLE(scip, maxweights[i], 0.0));
1725 }
1726 }
1727 assert(nterms == graph->terms);
1728 graph->stp_type = STP_MWCSP;
1729
1730 SCIP_CALL( graph_pc_2pc(scip, graph) );
1731 assert(graph->stp_type == STP_MWCSP);
1732
1733 SCIPdebugMessage("Transformed to MW \n");
1734
1735 return SCIP_OKAY;
1736 }
1737
1738
1739
1740 /** alters the graph for RMWCS problems */
graph_pc_2rmw(SCIP * scip,GRAPH * graph)1741 SCIP_RETCODE graph_pc_2rmw(
1742 SCIP* scip, /**< SCIP data structure */
1743 GRAPH* graph /**< the graph */
1744 )
1745 {
1746 SCIP_Real* maxweights;
1747 int i;
1748 int root;
1749 int nnodes;
1750 int npterms;
1751 int nrterms;
1752 int maxgrad;
1753
1754 assert(scip != NULL);
1755 assert(graph != NULL);
1756 assert(graph->cost != NULL);
1757
1758 root = -1;
1759 maxgrad = -1;
1760 npterms = 0;
1761 nrterms = 0;
1762 nnodes = graph->knots;
1763 maxweights = graph->prize;
1764
1765 assert(maxweights != NULL);
1766
1767 /* count number of terminals, modify incoming edges for non-terminals */
1768 for( i = 0; i < nnodes; i++ )
1769 {
1770 if( SCIPisLT(scip, maxweights[i], 0.0) )
1771 {
1772 for( int e = graph->inpbeg[i]; e != EAT_LAST; e = graph->ieat[e] )
1773 graph->cost[e] -= maxweights[i];
1774 }
1775 else if( SCIPisGE(scip, maxweights[i], FARAWAY) )
1776 {
1777 assert(Is_term(graph->term[i]));
1778 if( graph->grad[i] > maxgrad )
1779 {
1780 root = i;
1781 maxgrad = graph->grad[i];
1782 }
1783
1784 nrterms++;
1785 }
1786 else if( SCIPisGT(scip, maxweights[i], 0.0) )
1787 {
1788 graph_knot_chg(graph, i, 0);
1789 npterms++;
1790 }
1791 }
1792
1793 assert(root >= 0);
1794 assert(graph->terms == (npterms + nrterms));
1795
1796 graph->norgmodeledges = graph->edges;
1797 graph->norgmodelknots = nnodes;
1798 graph->source = root;
1799
1800 /* for each terminal, except for the root, one node and three edges (i.e. six arcs) are to be added */
1801 SCIP_CALL( graph_resize(scip, graph, (graph->ksize + npterms), (graph->esize + npterms * 4) , -1) );
1802
1803 /* create a new nodes */
1804 for( int k = 0; k < npterms; k++ )
1805 graph_knot_add(graph, -1);
1806
1807 /* allocate and initialize term2edge array */
1808 graph_pc_init(scip, graph, -1, graph->knots);
1809 assert(graph->term2edge != NULL);
1810
1811 i = 0;
1812 for( int k = 0; k < nnodes; ++k )
1813 {
1814 /* is the kth node a terminal other than the root? */
1815 if( Is_term(graph->term[k]) && SCIPisLT(scip, maxweights[k], FARAWAY) )
1816 {
1817 /* the copied node */
1818 const int node = nnodes + i;
1819 i++;
1820
1821 /* switch the terminal property, mark k */
1822 graph_knot_chg(graph, k, -2);
1823 graph_knot_chg(graph, node, 0);
1824 assert(SCIPisGE(scip, maxweights[k], 0.0));
1825
1826 /* add one edge going from the root to the 'copied' terminal and one going from the former terminal to its copy */
1827 graph_edge_add(scip, graph, root, node, maxweights[k], FARAWAY);
1828
1829 graph->term2edge[k] = graph->edges;
1830 graph->term2edge[node] = graph->edges + 1;
1831 assert(graph->edges + 1 == flipedge(graph->edges));
1832
1833 graph_edge_add(scip, graph, k, node, 0.0, FARAWAY);
1834
1835 assert(graph->head[graph->term2edge[k]] == node);
1836 assert(graph->head[graph->term2edge[node]] == k);
1837 }
1838 }
1839
1840 assert(i == npterms);
1841 graph->extended = TRUE;
1842 graph->stp_type = STP_RMWCSP;
1843
1844 SCIPdebugMessage("Transformed to RMW \n");
1845
1846 return SCIP_OKAY;
1847 }
1848
1849 /** transforms MWCSP to RMWCSP if possible */
graph_pc_mw2rmw(SCIP * scip,GRAPH * graph,SCIP_Real prizesum)1850 SCIP_RETCODE graph_pc_mw2rmw(
1851 SCIP* scip, /**< SCIP data structure */
1852 GRAPH* graph, /**< the graph */
1853 SCIP_Real prizesum /**< sum of positive prizes */
1854 )
1855 {
1856 int e;
1857 int p;
1858 int newroot;
1859 int maxgrad;
1860 const int root = graph->source;
1861
1862 assert(scip != NULL);
1863 assert(graph != NULL);
1864 assert(graph->term2edge != NULL);
1865 assert(graph->extended);
1866
1867 newroot = -1;
1868 maxgrad = -1;
1869
1870 e = graph->outbeg[root];
1871 while( e != EAT_LAST )
1872 {
1873 const int enext = graph->oeat[e];
1874 if( SCIPisGE(scip, graph->cost[e], prizesum) )
1875 {
1876 int e2;
1877 const int k = graph->head[e];
1878
1879 assert(Is_term(graph->term[k]));
1880 assert(graph->grad[k] == 2);
1881
1882 for( e2 = graph->outbeg[k]; e2 != EAT_LAST; e2 = graph->oeat[e2] )
1883 if( graph->head[e2] != root )
1884 break;
1885
1886 p = graph->head[e2];
1887 assert(e2 == graph->term2edge[k]);
1888
1889 assert(Is_pterm(graph->term[p]));
1890 assert(SCIPisGE(scip, graph->prize[p], prizesum));
1891
1892 /* delete terminal */
1893 graph_knot_chg(graph, k, -1);
1894 while( graph->outbeg[k] != EAT_LAST )
1895 graph_edge_del(scip, graph, graph->outbeg[k], TRUE);
1896
1897 graph->term2edge[k] = -1;
1898 graph->term2edge[p] = -1;
1899
1900 graph_knot_chg(graph, p, 0);
1901
1902 if( graph->grad[p] > maxgrad )
1903 {
1904 newroot = p;
1905 maxgrad = graph->grad[p];
1906 }
1907 }
1908 e = enext;
1909 }
1910
1911 /* is there a new root? */
1912 if( newroot >= 0 )
1913 {
1914 graph->source = newroot;
1915
1916 e = graph->outbeg[root];
1917 while( e != EAT_LAST )
1918 {
1919 const int enext = graph->oeat[e];
1920 const int k = graph->head[e];
1921 if( Is_term(graph->term[k]) && !SCIPisZero(scip, graph->cost[e]) )
1922 {
1923 (void) graph_edge_redirect(scip, graph, e, newroot, k, graph->cost[e], TRUE);
1924 graph->cost[flipedge(e)] = FARAWAY;
1925 }
1926 e = enext;
1927 }
1928
1929 /* delete old root */
1930 graph_knot_chg(graph, root, -1);
1931 while( graph->outbeg[root] != EAT_LAST )
1932 graph_edge_del(scip, graph, graph->outbeg[root], TRUE);
1933
1934 graph->stp_type = STP_RMWCSP;
1935
1936 }
1937
1938 SCIPdebugMessage("Transformed MW to RMW \n");
1939
1940 return SCIP_OKAY;
1941 }
1942
1943
1944 /** delete a terminal for a (rooted) prize-collecting problem */
graph_pc_deleteTerm(SCIP * scip,GRAPH * g,int i)1945 int graph_pc_deleteTerm(
1946 SCIP* scip, /**< SCIP data structure */
1947 GRAPH* g, /**< graph data structure */
1948 int i /**< index of the terminal */
1949 )
1950 {
1951 int e;
1952 int t;
1953 const int grad = g->grad[i];
1954
1955 assert(g != NULL);
1956 assert(scip != NULL);
1957 assert(Is_term(g->term[i]));
1958
1959 t = UNKNOWN;
1960
1961 /* delete terminal */
1962
1963 assert(g->term2edge[i] != -1);
1964 graph_pc_knot2nonTerm(g, i);
1965 g->mark[i] = FALSE;
1966
1967 while( (e = g->outbeg[i]) != EAT_LAST )
1968 {
1969 const int i1 = g->head[e];
1970
1971 if( Is_pterm(g->term[i1]) && g->source != i1 )
1972 t = g->head[e];
1973 graph_edge_del(scip, g, e, TRUE);
1974 }
1975
1976 assert(g->grad[i] == 0);
1977 assert(t != UNKNOWN);
1978 assert(g->term2edge != NULL);
1979
1980 /* delete artificial terminal */
1981
1982 graph_pc_knot2nonTerm(g, t);
1983
1984 while( g->outbeg[t] != EAT_LAST )
1985 graph_edge_del(scip, g, g->outbeg[t], TRUE);
1986
1987 return grad + 2;
1988 }
1989
1990
1991 /** subtract a given sum from the prize of a terminal */
graph_pc_subtractPrize(SCIP * scip,GRAPH * g,SCIP_Real cost,int i)1992 void graph_pc_subtractPrize(
1993 SCIP* scip, /**< SCIP data structure */
1994 GRAPH* g, /**< the graph */
1995 SCIP_Real cost, /**< cost to be subtracted */
1996 int i /**< the terminal */
1997 )
1998 {
1999 int e;
2000 int j;
2001
2002 assert(scip != NULL);
2003 assert(g != NULL);
2004
2005 if( g->stp_type == STP_RPCSPG && i == g->source )
2006 return;
2007
2008 g->prize[i] -= cost;
2009 for( e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e] )
2010 if( Is_pterm(g->term[g->head[e]]) )
2011 break;
2012
2013 assert(e != EAT_LAST);
2014
2015 j = g->head[e];
2016
2017 assert(j != g->source);
2018 assert(!g->mark[j]);
2019
2020 for( e = g->inpbeg[j]; e != EAT_LAST; e = g->ieat[e] )
2021 if( g->source == g->tail[e] )
2022 break;
2023
2024 assert(e != EAT_LAST);
2025 assert(!g->mark[g->tail[e]] || g->stp_type == STP_RPCSPG);
2026
2027 g->cost[e] -= cost;
2028
2029 assert(g->stp_type == STP_MWCSP || g->stp_type == STP_RMWCSP || SCIPisGE(scip, g->prize[i], 0.0));
2030 assert(SCIPisEQ(scip, g->prize[i], g->cost[e]));
2031 assert(SCIPisGE(scip, g->prize[i], 0.0) || g->stp_type == STP_MWCSP);
2032 }
2033
2034 /** change prize of a terminal */
graph_pc_chgPrize(SCIP * scip,GRAPH * g,SCIP_Real newprize,int i)2035 void graph_pc_chgPrize(
2036 SCIP* scip, /**< SCIP data structure */
2037 GRAPH* g, /**< the graph */
2038 SCIP_Real newprize, /**< prize to be subtracted */
2039 int i /**< the terminal */
2040 )
2041 {
2042 int e;
2043 int j;
2044
2045 assert(scip != NULL);
2046 assert(g != NULL);
2047 assert(newprize > 0.0);
2048
2049 if( g->stp_type == STP_RPCSPG && i == g->source )
2050 return;
2051
2052 g->prize[i] = newprize;
2053 for( e = g->outbeg[i]; e != EAT_LAST; e = g->oeat[e] )
2054 if( Is_pterm(g->term[g->head[e]]) )
2055 break;
2056
2057 assert(e != EAT_LAST);
2058
2059 j = g->head[e];
2060
2061 assert(j != g->source);
2062 assert(!g->mark[j]);
2063
2064 for( e = g->inpbeg[j]; e != EAT_LAST; e = g->ieat[e] )
2065 if( g->source == g->tail[e] )
2066 break;
2067
2068 assert(e != EAT_LAST);
2069 assert(!g->mark[g->tail[e]] || g->stp_type == STP_RPCSPG);
2070
2071 g->cost[e] = newprize;
2072
2073 assert(g->stp_type == STP_MWCSP || g->stp_type == STP_RMWCSP || SCIPisGE(scip, g->prize[i], 0.0));
2074 assert(SCIPisEQ(scip, g->prize[i], g->cost[e]));
2075 assert(SCIPisGE(scip, g->prize[i], 0.0) || g->stp_type == STP_MWCSP);
2076 }
2077
2078 /** contract ancestors of an edge of (rooted) prize-collecting Steiner tree problem or maximum-weight connected subgraph problem */
graph_pc_contractEdgeAncestors(SCIP * scip,GRAPH * g,int t,int s,int ets)2079 SCIP_RETCODE graph_pc_contractEdgeAncestors(
2080 SCIP* scip, /**< SCIP data structure */
2081 GRAPH* g, /**< the graph */
2082 int t, /**< tail node to be contracted (surviving node) */
2083 int s, /**< head node to be contracted */
2084 int ets /**< edge from t to s or -1 */
2085 )
2086 {
2087 assert(g != NULL);
2088 assert(scip != NULL);
2089
2090 if( ets < 0 )
2091 {
2092 for( ets = g->outbeg[t]; ets != EAT_LAST; ets = g->oeat[ets] )
2093 if( g->head[ets] == s )
2094 break;
2095 assert(ets >= 0);
2096 }
2097
2098 SCIP_CALL(SCIPintListNodeAppendCopy(scip, &(g->pcancestors[s]), g->ancestors[ets], NULL));
2099 SCIP_CALL(SCIPintListNodeAppendCopy(scip, &(g->pcancestors[t]), g->ancestors[ets], NULL));
2100
2101 return SCIP_OKAY;
2102 }
2103
2104 /** contract an edge of (rooted) prize-collecting Steiner tree problem or maximum-weight connected subgraph problem */
graph_pc_contractEdge(SCIP * scip,GRAPH * g,int * solnode,int t,int s,int i)2105 SCIP_RETCODE graph_pc_contractEdge(
2106 SCIP* scip, /**< SCIP data structure */
2107 GRAPH* g, /**< the graph */
2108 int* solnode, /**< solution nodes or NULL */
2109 int t, /**< tail node to be contracted (surviving node) */
2110 int s, /**< head node to be contracted */
2111 int i /**< terminal to add offset to */
2112 )
2113 {
2114 int ets;
2115
2116 assert(g != NULL);
2117 assert(scip != NULL);
2118 assert(Is_term(g->term[i]));
2119
2120 /* get edge from t to s */
2121 for( ets = g->outbeg[t]; ets != EAT_LAST; ets = g->oeat[ets] )
2122 if( g->head[ets] == s )
2123 break;
2124
2125 assert(ets != EAT_LAST);
2126
2127 SCIP_CALL( graph_pc_contractEdgeAncestors(scip, g, t, s, ets) );
2128
2129 /* are both endpoints of the edge to be contracted terminals? */
2130 if( Is_term(g->term[t]) && Is_term(g->term[s]) )
2131 {
2132 int e;
2133 int j;
2134
2135 /* get edge from s to its artificial terminal */
2136 for( e = g->outbeg[s]; e != EAT_LAST; e = g->oeat[e] )
2137 if( Is_pterm(g->term[g->head[e]]) )
2138 break;
2139
2140 assert(e != EAT_LAST);
2141 assert(g->pcancestors != NULL);
2142
2143 /* artificial terminal to s */
2144 j = g->head[e];
2145
2146 assert(j != g->source);
2147 assert(!g->mark[j]);
2148 assert(g->term2edge != NULL);
2149
2150 /* delete edge and unmark artificial terminal */
2151 graph_knot_chg(g, j, -1);
2152 graph_edge_del(scip, g, e, TRUE);
2153 g->term2edge[j] = -1;
2154
2155 /* delete remaining incident edge of artificial terminal */
2156 e = g->inpbeg[j];
2157
2158 assert(e != EAT_LAST);
2159 assert(g->source == g->tail[e] || g->source == j);
2160 assert(SCIPisEQ(scip, g->prize[s], g->cost[e]));
2161
2162 graph_pc_subtractPrize(scip, g, g->cost[ets] - g->prize[s], i);
2163 graph_edge_del(scip, g, e, TRUE);
2164
2165 assert(g->inpbeg[j] == EAT_LAST);
2166
2167 /* contract s into t */
2168 SCIP_CALL( graph_knot_contract(scip, g, solnode, t, s) );
2169 g->term2edge[s] = -1;
2170
2171 assert(g->grad[s] == 0);
2172
2173 SCIPdebugMessage("PC contract: %d, %d \n", t, s);
2174 }
2175 else
2176 {
2177 if( g->stp_type != STP_MWCSP && g->stp_type != STP_RMWCSP )
2178 graph_pc_subtractPrize(scip, g, g->cost[ets], i);
2179 else
2180 graph_pc_subtractPrize(scip, g, -(g->prize[s]), i);
2181 SCIP_CALL( graph_knot_contract(scip, g, solnode, t, s) );
2182 }
2183 return SCIP_OKAY;
2184 }
2185
2186
2187 /** is this graph a prize-collecting or maximum-weight variant? */
graph_pc_isPcMw(const GRAPH * g)2188 SCIP_Bool graph_pc_isPcMw(
2189 const GRAPH* g /**< the graph */
2190 )
2191 {
2192 const int type = g->stp_type;
2193 assert(g != NULL);
2194
2195 return (type == STP_PCSPG || type == STP_RPCSPG || type == STP_MWCSP || type == STP_RMWCSP);
2196 }
2197
2198
2199 /** add a vertex */
graph_knot_add(GRAPH * p,int term)2200 void graph_knot_add(
2201 GRAPH* p, /**< the graph */
2202 int term /**< terminal property */
2203 )
2204 {
2205 assert(p != NULL);
2206 assert(p->ksize > p->knots);
2207 assert(term < p->layers);
2208
2209 p->term [p->knots] = term;
2210 p->mark [p->knots] = TRUE;
2211 p->grad [p->knots] = 0;
2212 p->inpbeg[p->knots] = EAT_LAST;
2213 p->outbeg[p->knots] = EAT_LAST;
2214
2215 if( Is_term(term) )
2216 p->terms++;
2217
2218 p->knots++;
2219 }
2220
2221 /** change terminal property of a vertex */
graph_knot_chg(GRAPH * p,int node,int term)2222 void graph_knot_chg(
2223 GRAPH* p, /**< the graph */
2224 int node, /**< node to be changed */
2225 int term /**< terminal property */
2226 )
2227 {
2228 assert(p != NULL);
2229 assert(node >= 0);
2230 assert(node < p->knots);
2231 assert(term < p->layers);
2232
2233 if( term != p->term[node] )
2234 {
2235 if( Is_term(p->term[node]) )
2236 p->terms--;
2237
2238 p->term[node] = term;
2239
2240 if( Is_term(p->term[node]) )
2241 p->terms++;
2242 }
2243 }
2244
2245 /** delete node */
graph_knot_del(SCIP * scip,GRAPH * g,int k,SCIP_Bool freeancestors)2246 void graph_knot_del(
2247 SCIP* scip, /**< SCIP data structure */
2248 GRAPH* g, /**< the graph */
2249 int k, /**< the node */
2250 SCIP_Bool freeancestors /**< free edge ancestors? */
2251 )
2252 {
2253 assert(g != NULL);
2254 assert(k >= 0);
2255 assert(k < g->knots);
2256
2257 while( g->outbeg[k] != EAT_LAST )
2258 graph_edge_del(scip, g, g->outbeg[k], freeancestors);
2259 }
2260
2261 /** pseudo delete node, i.e. reconnect neighbors; maximum degree of 4! */
graph_knot_delPseudo(SCIP * scip,GRAPH * g,const SCIP_Real * edgecosts,const SCIP_Real * cutoffs,const SCIP_Real * cutoffsrev,int vertex,SCIP_Bool * success)2262 SCIP_RETCODE graph_knot_delPseudo(
2263 SCIP* scip, /**< SCIP data structure */
2264 GRAPH* g, /**< the graph */
2265 const SCIP_Real* edgecosts, /**< edge costs for cutoff */
2266 const SCIP_Real* cutoffs, /**< cutoff values for each incident edge */
2267 const SCIP_Real* cutoffsrev, /**< revere cutoff values (or NULL if undirected) */
2268 int vertex, /**< the vertex */
2269 SCIP_Bool* success /**< has node been pseudo-eliminated? */
2270 )
2271 {
2272 IDX* ancestors[STP_DELPSEUDO_MAXGRAD];
2273 IDX* revancestors[STP_DELPSEUDO_MAXGRAD];
2274 SCIP_Real ecost[STP_DELPSEUDO_MAXGRAD];
2275 SCIP_Real ecostrev[STP_DELPSEUDO_MAXGRAD];
2276 SCIP_Real ecostreal[STP_DELPSEUDO_MAXGRAD];
2277 int incedge[STP_DELPSEUDO_MAXGRAD];
2278 int adjvert[STP_DELPSEUDO_MAXGRAD];
2279 int neigbedge[STP_DELPSEUDO_MAXNEDGES];
2280 int edgecount;
2281 int nspareedges;
2282 int replacecount;
2283 const int degree = g->grad[vertex];
2284
2285 assert(scip != NULL);
2286 assert(success != NULL);
2287 assert(g != NULL);
2288 assert(vertex >= 0);
2289 assert(vertex < g->knots);
2290 assert(degree <= STP_DELPSEUDO_MAXGRAD);
2291
2292 #ifndef NDEBUG
2293 {
2294 int sum = 0;
2295 for( int i = 1; i < STP_DELPSEUDO_MAXGRAD; i++ )
2296 sum += i;
2297 assert(sum == STP_DELPSEUDO_MAXNEDGES);
2298 }
2299 #endif
2300
2301 *success = TRUE;
2302
2303 if( degree <= 1 )
2304 return SCIP_OKAY;
2305
2306 nspareedges = degree; /* todo */
2307
2308 edgecount = 0;
2309
2310 for( int i = 0; i < STP_DELPSEUDO_MAXNEDGES; i++ )
2311 neigbedge[i] = -1;
2312
2313 /* save old edges */
2314 for( int e = g->outbeg[vertex]; e != EAT_LAST; e = g->oeat[e] )
2315 {
2316 assert(e >= 0);
2317
2318 incedge[edgecount] = e;
2319 ecostreal[edgecount] = g->cost[e];
2320 ecost[edgecount] = edgecosts[e];
2321 ecostrev[edgecount] = edgecosts[flipedge(e)];
2322
2323 adjvert[edgecount++] = g->head[e];
2324
2325 assert(edgecount <= STP_DELPSEUDO_MAXGRAD);
2326 }
2327
2328 assert(edgecount == degree);
2329 edgecount = 0;
2330 replacecount = 0;
2331
2332 /* check whether there are enough spare edges */
2333 for( int i = 0; i < degree - 1; i++ )
2334 {
2335 const int adjvertex = adjvert[i];
2336 for( int j = i + 1; j < degree; j++ )
2337 {
2338 int e;
2339 const SCIP_Bool cutoff = cutoffEdge(scip, cutoffs, cutoffsrev, ecost, ecostrev, i, j, edgecount);
2340
2341 assert(edgecount < STP_DELPSEUDO_MAXNEDGES);
2342
2343 edgecount++;
2344
2345 /* can edge be discarded? */
2346 if( cutoff )
2347 continue;
2348
2349 /* check whether edge already exists */
2350 for( e = g->outbeg[adjvertex]; e != EAT_LAST; e = g->oeat[e] )
2351 if( g->head[e] == adjvert[j] )
2352 {
2353 assert(e >= 0);
2354 neigbedge[edgecount - 1] = e;
2355 break;
2356 }
2357
2358 if( e != EAT_LAST )
2359 continue;
2360
2361 if( ++replacecount > nspareedges )
2362 {
2363 *success = FALSE;
2364 return SCIP_OKAY;
2365 }
2366 }
2367 }
2368
2369 for( int i = 0; i < degree; i++ )
2370 {
2371 const int e = incedge[i];
2372 ancestors[i] = NULL;
2373 revancestors[i] = NULL;
2374
2375 SCIP_CALL(SCIPintListNodeAppendCopy(scip, &(ancestors[i]), g->ancestors[e], NULL));
2376 SCIP_CALL(SCIPintListNodeAppendCopy(scip, &(revancestors[i]), g->ancestors[flipedge(e)], NULL));
2377 }
2378
2379 /* replace edges */
2380 edgecount = 0;
2381 replacecount = 0;
2382 for( int i = 0; i < degree - 1; i++ )
2383 {
2384 for( int j = i + 1; j < degree; j++ )
2385 {
2386 const SCIP_Bool cutoff = cutoffEdge(scip, cutoffs, cutoffsrev, ecost, ecostrev, i, j, edgecount);
2387
2388 assert(edgecount < STP_DELPSEUDO_MAXNEDGES);
2389
2390 edgecount++;
2391
2392 /* do we need to insert edge at all? */
2393 if( !cutoff )
2394 {
2395 const SCIP_Real newcost = ecostreal[i] + ecostreal[j];
2396 const int oldedge = incedge[(replacecount == nspareedges) ? replacecount - 1 : replacecount];
2397 #ifndef NDEBUG
2398 const int oldtail = g->tail[oldedge];
2399 const int oldhead = g->head[oldedge];
2400 #endif
2401 assert(replacecount <= nspareedges);
2402 assert(replacecount < nspareedges || neigbedge[edgecount - 1] >= 0);
2403
2404 SCIP_CALL( graph_edge_reinsert(scip, g, oldedge, adjvert[i], adjvert[j], newcost, ancestors[i], ancestors[j], revancestors[i], revancestors[j], FALSE));
2405
2406 /* does no edge exist? */
2407 if( neigbedge[edgecount - 1] < 0 )
2408 replacecount++;
2409 #ifndef NDEBUG
2410 else
2411 {
2412 assert(oldtail == g->tail[oldedge]);
2413 assert(oldhead == g->head[oldedge]);
2414 }
2415 #endif
2416 }
2417 }
2418 }
2419
2420 /* delete remaining edges */
2421 graph_knot_del(scip, g, vertex, TRUE);
2422
2423 for( int i = 0; i < degree; i++ )
2424 {
2425 SCIPintListNodeFree(scip, &(ancestors[i]));
2426 SCIPintListNodeFree(scip, &(revancestors[i]));
2427 }
2428
2429 return SCIP_OKAY;
2430 }
2431
2432 /** contract an edge, given by its endpoints */
graph_knot_contract(SCIP * scip,GRAPH * p,int * solnode,int t,int s)2433 SCIP_RETCODE graph_knot_contract(
2434 SCIP* scip, /**< SCIP data structure */
2435 GRAPH* p, /**< graph data structure */
2436 int* solnode, /**< node array to mark whether an node is part of a given solution (CONNECT),
2437 or NULL */
2438 int t, /**< tail node to be contracted */
2439 int s /**< head node to be contracted */
2440 )
2441 {
2442 SCIP_Real* incost = NULL;
2443 SCIP_Real* outcost = NULL;
2444 IDX** ancestors = NULL;
2445 IDX** revancestors = NULL;
2446 int* mark = NULL;
2447 int* edge = NULL;
2448 int* knot = NULL;
2449 int slc = 0;
2450 int i;
2451 int et;
2452 int anti;
2453 int es;
2454 int head;
2455 int tail;
2456 int sgrad;
2457
2458 assert(p != NULL);
2459 assert(t >= 0);
2460 assert(t < p->knots);
2461 assert(s >= 0);
2462 assert(s < p->knots);
2463 assert(s != t);
2464 assert(scip != NULL);
2465 assert(p->grad[s] > 0);
2466 assert(p->grad[t] > 0);
2467 assert(p->layers == 1);
2468
2469 /* save solution */
2470 if( solnode != NULL )
2471 if( solnode[s] == CONNECT )
2472 solnode[t] = CONNECT;
2473
2474 /* change terminal property */
2475 if( Is_term(p->term[s]) )
2476 {
2477 graph_knot_chg(p, t, p->term[s]);
2478 graph_knot_chg(p, s, -1);
2479 }
2480
2481 /* retain root */
2482 if( p->source == s )
2483 p->source = t;
2484
2485 sgrad = p->grad[s];
2486 if( sgrad >= 2 )
2487 {
2488 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &incost, sgrad - 1) );
2489 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &outcost, sgrad - 1) );
2490 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &mark, sgrad - 1) );
2491 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &edge, sgrad - 1) );
2492 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &knot, sgrad - 1) );
2493 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &ancestors, sgrad - 1) );
2494 SCIP_CALL( SCIPallocBlockMemoryArray(scip, &revancestors, sgrad - 1) );
2495 }
2496
2497 /* store edges to be moved/removed */
2498 for( es = p->outbeg[s]; es != EAT_LAST; es = p->oeat[es] )
2499 {
2500 assert(p->tail[es] == s);
2501
2502 if( p->head[es] != t )
2503 {
2504 assert(ancestors != NULL);
2505 assert(revancestors != NULL);
2506 assert(mark != NULL);
2507 assert(incost != NULL);
2508 assert(outcost != NULL);
2509 assert(edge != NULL);
2510 assert(knot != NULL);
2511
2512 ancestors[slc] = NULL;
2513 SCIP_CALL( SCIPintListNodeAppendCopy(scip, &(ancestors[slc]), p->ancestors[es], NULL) );
2514 revancestors[slc] = NULL;
2515 SCIP_CALL( SCIPintListNodeAppendCopy(scip, &(revancestors[slc]), p->ancestors[Edge_anti(es)], NULL) );
2516
2517 mark[slc] = FALSE;
2518 edge[slc] = es;
2519 knot[slc] = p->head[es];
2520 outcost[slc] = p->cost[es];
2521 incost[slc] = p->cost[Edge_anti(es)];
2522 slc++;
2523
2524 assert(slc < sgrad);
2525 }
2526 }
2527
2528 assert(slc == sgrad - 1);
2529
2530 /* traverse edges */
2531 for( i = 0; i < slc; i++ )
2532 {
2533 const int ihead = knot[i];
2534 assert(knot != NULL && outcost != NULL && incost != NULL && mark != NULL);
2535
2536 /* search for an edge out of t with same head as current edge */
2537
2538 if( p->grad[ihead] >= p->grad[t] )
2539 {
2540 for( et = p->outbeg[t]; et >= 0; et = p->oeat[et] )
2541 if( p->head[et] == ihead )
2542 break;
2543 }
2544 else
2545 {
2546 for( et = p->inpbeg[ihead]; et >= 0; et = p->ieat[et] )
2547 if( p->tail[et] == t )
2548 break;
2549 }
2550
2551 /* does such an edge not exist? */
2552 if( et == EAT_LAST )
2553 {
2554 mark[i] = TRUE;
2555 }
2556 else
2557 {
2558 assert(et != EAT_LAST);
2559
2560 /* This is for nodes with edges to s and t.
2561 * Need to adjust the out and in costs of the edge
2562 */
2563 if( p->cost[et] > outcost[i] )
2564 {
2565 SCIPintListNodeFree(scip, &((p->ancestors)[et]));
2566 assert(ancestors != NULL);
2567 SCIP_CALL( SCIPintListNodeAppendCopy(scip, &((p->ancestors)[et]), ancestors[i], NULL) );
2568
2569 p->cost[et] = outcost[i];
2570 }
2571 if( p->cost[Edge_anti(et)] > incost[i] )
2572 {
2573 anti = Edge_anti(et);
2574 SCIPintListNodeFree(scip, &(p->ancestors[anti]));
2575 assert(revancestors != NULL);
2576 assert(incost != NULL);
2577 SCIP_CALL( SCIPintListNodeAppendCopy(scip, &((p->ancestors)[anti]), revancestors[i], NULL) );
2578 p->cost[anti] = incost[i];
2579 }
2580 }
2581 }
2582
2583 /* insert edges */
2584 for( i = 0; i < slc; i++ )
2585 {
2586 assert(mark != NULL);
2587 if( mark[i] )
2588 {
2589 es = p->outbeg[s];
2590
2591 assert(es != EAT_LAST);
2592 assert(ancestors != NULL);
2593 assert(revancestors != NULL);
2594 assert(ancestors[i] != NULL);
2595 assert(revancestors[i] != NULL);
2596 assert(knot != NULL);
2597 assert(outcost != NULL);
2598 assert(incost != NULL);
2599 SCIPintListNodeFree(scip, &(p->ancestors[es]));
2600 SCIP_CALL( SCIPintListNodeAppendCopy(scip, &(p->ancestors[es]), ancestors[i], NULL) );
2601
2602 graph_edge_del(scip, p, es, FALSE);
2603
2604 head = knot[i];
2605 tail = t;
2606
2607 p->grad[head]++;
2608 p->grad[tail]++;
2609
2610 p->cost[es] = outcost[i];
2611 p->tail[es] = tail;
2612 p->head[es] = head;
2613 p->ieat[es] = p->inpbeg[head];
2614 p->oeat[es] = p->outbeg[tail];
2615 p->inpbeg[head] = es;
2616 p->outbeg[tail] = es;
2617
2618 es = Edge_anti(es);
2619 SCIPintListNodeFree(scip, &(p->ancestors[es]));
2620
2621 SCIP_CALL( SCIPintListNodeAppendCopy(scip, &(p->ancestors[es]), revancestors[i], NULL) );
2622
2623 p->cost[es] = incost[i];
2624 p->tail[es] = head;
2625 p->head[es] = tail;
2626 p->ieat[es] = p->inpbeg[tail];
2627 p->oeat[es] = p->outbeg[head];
2628 p->inpbeg[tail] = es;
2629 p->outbeg[head] = es;
2630 }
2631 }
2632
2633 /* delete remaining edges */
2634 while( p->outbeg[s] != EAT_LAST )
2635 {
2636 es = p->outbeg[s];
2637 SCIPintListNodeFree(scip, &(p->ancestors[es]));
2638 SCIPintListNodeFree(scip, &(p->ancestors[Edge_anti(es)]));
2639 graph_edge_del(scip, p, es, FALSE);
2640 }
2641
2642 if( sgrad >= 2 )
2643 {
2644 assert(ancestors != NULL);
2645 assert(revancestors != NULL);
2646 for( i = 0; i < slc; i++ )
2647 {
2648 SCIPintListNodeFree(scip, &(ancestors[i]));
2649 SCIPintListNodeFree(scip, &(revancestors[i]));
2650 }
2651 SCIPfreeBlockMemoryArray(scip, &revancestors, sgrad - 1);
2652 SCIPfreeBlockMemoryArray(scip, &ancestors, sgrad - 1);
2653 SCIPfreeBlockMemoryArray(scip, &knot, sgrad - 1);
2654 SCIPfreeBlockMemoryArray(scip, &edge, sgrad - 1);
2655 SCIPfreeBlockMemoryArray(scip, &mark, sgrad - 1);
2656 SCIPfreeBlockMemoryArray(scip, &outcost, sgrad - 1);
2657 SCIPfreeBlockMemoryArray(scip, &incost, sgrad - 1);
2658 }
2659 assert(p->grad[s] == 0);
2660 assert(p->outbeg[s] == EAT_LAST);
2661 assert(p->inpbeg[s] == EAT_LAST);
2662 return SCIP_OKAY;
2663 }
2664
2665 /** contract endpoint of lower degree into endpoint of higher degree */
graph_knot_contractLowdeg2High(SCIP * scip,GRAPH * g,int * solnode,int t,int s)2666 SCIP_RETCODE graph_knot_contractLowdeg2High(
2667 SCIP* scip, /**< SCIP data structure */
2668 GRAPH* g, /**< graph data structure */
2669 int* solnode, /**< node array to mark whether an node is part of a given solution (CONNECT),
2670 or NULL */
2671 int t, /**< tail node to be contracted */
2672 int s /**< head node to be contracted */
2673 )
2674 {
2675 assert(g != NULL);
2676
2677 if( g->grad[t] >= g->grad[s] )
2678 SCIP_CALL( graph_knot_contract(scip, g, solnode, t, s) );
2679 else
2680 SCIP_CALL( graph_knot_contract(scip, g, solnode, s, t) );
2681
2682 return SCIP_OKAY;
2683 }
2684
graph_edge_redirect(SCIP * scip,GRAPH * g,int eki,int k,int j,SCIP_Real cost,SCIP_Bool forcedelete)2685 int graph_edge_redirect(
2686 SCIP* scip, /**< SCIP data structure */
2687 GRAPH* g, /**< the graph */
2688 int eki, /**< the edge */
2689 int k, /**< new tail */
2690 int j, /**< new head */
2691 SCIP_Real cost, /**< new cost */
2692 SCIP_Bool forcedelete /**< delete edge eki if it is not used? */
2693 )
2694 {
2695 int e;
2696
2697 if( forcedelete )
2698 graph_edge_del(NULL, g, eki, FALSE);
2699
2700 for( e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
2701 {
2702 assert(g->tail[e] == k);
2703
2704 if( g->head[e] == j )
2705 break;
2706 }
2707
2708 /* does edge already exist? */
2709 if( e != EAT_LAST )
2710 {
2711 /* correct cost */
2712 if( SCIPisGT(scip, g->cost[e], cost) )
2713 {
2714 g->cost[e] = cost;
2715 g->cost[Edge_anti(e)] = cost;
2716 }
2717 else
2718 {
2719 e = -1;
2720 }
2721 }
2722 else
2723 {
2724 if( !forcedelete )
2725 graph_edge_del(NULL, g, eki, FALSE);
2726
2727 assert(g->oeat[eki] == EAT_FREE);
2728
2729 e = eki;
2730
2731 g->grad[k]++;
2732 g->grad[j]++;
2733
2734 g->cost[e] = cost;
2735 g->head[e] = j;
2736 g->tail[e] = k;
2737 g->ieat[e] = g->inpbeg[j];
2738 g->oeat[e] = g->outbeg[k];
2739 g->inpbeg[j] = e;
2740 g->outbeg[k] = e;
2741
2742 e = Edge_anti(eki);
2743
2744 g->cost[e] = cost;
2745 g->head[e] = k;
2746 g->tail[e] = j;
2747 g->ieat[e] = g->inpbeg[k];
2748 g->oeat[e] = g->outbeg[j];
2749 g->inpbeg[k] = e;
2750 g->outbeg[j] = e;
2751 return eki;
2752 }
2753 return e;
2754 }
2755
2756 /** reinsert an edge to replace two other edges */
graph_edge_reinsert(SCIP * scip,GRAPH * g,int e1,int k1,int k2,SCIP_Real cost,IDX * ancestors0,IDX * ancestors1,IDX * revancestors0,IDX * revancestors1,SCIP_Bool forcedelete)2757 SCIP_RETCODE graph_edge_reinsert(
2758 SCIP* scip, /**< SCIP data structure */
2759 GRAPH* g, /**< the graph */
2760 int e1, /**< edge to reinsert */
2761 int k1, /**< tail */
2762 int k2, /**< head */
2763 SCIP_Real cost, /**< edge cost */
2764 IDX* ancestors0, /**< ancestors of first edge */
2765 IDX* ancestors1, /**< ancestors of second edge */
2766 IDX* revancestors0, /**< reverse ancestors of first edge */
2767 IDX* revancestors1, /**< reverse ancestors of first edge */
2768 SCIP_Bool forcedelete /**< delete edge e1 if it is not used? */
2769 )
2770 {
2771 /* redirect; store new edge in n1 */
2772 const int n1 = graph_edge_redirect(scip, g, e1, k1, k2, cost, forcedelete);
2773
2774 if( n1 >= 0 )
2775 {
2776 SCIPintListNodeFree(scip, &(g->ancestors[n1]));
2777 SCIPintListNodeFree(scip, &(g->ancestors[Edge_anti(n1)]));
2778
2779 SCIP_CALL( SCIPintListNodeAppendCopy(scip, &(g->ancestors[n1]), revancestors0, NULL) );
2780 SCIP_CALL( SCIPintListNodeAppendCopy(scip, &(g->ancestors[n1]), ancestors1, NULL) );
2781
2782 SCIP_CALL( SCIPintListNodeAppendCopy(scip, &(g->ancestors[Edge_anti(n1)]), ancestors0, NULL) );
2783 SCIP_CALL( SCIPintListNodeAppendCopy(scip, &(g->ancestors[Edge_anti(n1)]), revancestors1, NULL) );
2784 }
2785 return SCIP_OKAY;
2786 }
2787
2788
2789 /** add a new edge to the graph */
graph_edge_add(SCIP * scip,GRAPH * g,int tail,int head,SCIP_Real cost1,SCIP_Real cost2)2790 void graph_edge_add(
2791 SCIP* scip, /**< SCIP data structure */
2792 GRAPH* g, /**< the graph */
2793 int tail, /**< tail of the new edge */
2794 int head, /**< head of the new edge*/
2795 SCIP_Real cost1, /**< tail to head cost */
2796 SCIP_Real cost2 /**< head to tail cost */
2797 )
2798 {
2799 int e;
2800
2801 assert(g != NULL);
2802 assert(SCIPisGE(scip, cost1, 0.0) || SCIPisEQ(scip, cost1, (double) UNKNOWN));
2803 assert(SCIPisGE(scip, cost2, 0.0) || SCIPisEQ(scip, cost2, (double) UNKNOWN));
2804 assert(tail >= 0);
2805 assert(tail < g->knots);
2806 assert(head >= 0);
2807 assert(head < g->knots);
2808
2809 assert(g->esize >= g->edges + 2);
2810
2811 e = g->edges;
2812
2813 g->grad[head]++;
2814 g->grad[tail]++;
2815
2816 if( cost1 != UNKNOWN )
2817 g->cost[e] = cost1;
2818 g->tail[e] = tail;
2819 g->head[e] = head;
2820 g->ieat[e] = g->inpbeg[head];
2821 g->oeat[e] = g->outbeg[tail];
2822 g->inpbeg[head] = e;
2823 g->outbeg[tail] = e;
2824
2825 e++;
2826
2827 if( cost2 != UNKNOWN )
2828 g->cost[e] = cost2;
2829 g->tail[e] = head;
2830 g->head[e] = tail;
2831 g->ieat[e] = g->inpbeg[tail];
2832 g->oeat[e] = g->outbeg[head];
2833 g->inpbeg[tail] = e;
2834 g->outbeg[head] = e;
2835
2836 g->edges += 2;
2837 }
2838
2839
2840 /** delete an edge */
graph_edge_del(SCIP * scip,GRAPH * g,int e,SCIP_Bool freeancestors)2841 void graph_edge_del(
2842 SCIP* scip, /**< SCIP data structure */
2843 GRAPH* g, /**< the graph */
2844 int e, /**< the edge */
2845 SCIP_Bool freeancestors /**< free edge ancestors? */
2846 )
2847 {
2848 assert(g != NULL);
2849 assert(e >= 0);
2850 assert(e < g->edges);
2851
2852 if( freeancestors )
2853 {
2854 assert(scip != NULL);
2855 SCIPintListNodeFree(scip, &((g->ancestors)[e]));
2856 SCIPintListNodeFree(scip, &((g->ancestors)[Edge_anti(e)]));
2857 }
2858
2859 /* delete first arc */
2860 e -= e % 2;
2861 assert(g->head[e] == g->tail[e + 1]);
2862 assert(g->tail[e] == g->head[e + 1]);
2863
2864 g->grad[g->head[e]]--;
2865 g->grad[g->tail[e]]--;
2866
2867 removeEdge(g, e);
2868
2869 assert(g->ieat[e] != EAT_FREE);
2870 assert(g->ieat[e] != EAT_HIDE);
2871 assert(g->oeat[e] != EAT_FREE);
2872 assert(g->oeat[e] != EAT_HIDE);
2873
2874 g->ieat[e] = EAT_FREE;
2875 g->oeat[e] = EAT_FREE;
2876
2877 /* delete second arc */
2878 e++;
2879 removeEdge(g, e);
2880
2881 assert(g->ieat[e] != EAT_FREE);
2882 assert(g->ieat[e] != EAT_HIDE);
2883 assert(g->oeat[e] != EAT_FREE);
2884 assert(g->oeat[e] != EAT_HIDE);
2885
2886 g->ieat[e] = EAT_FREE;
2887 g->oeat[e] = EAT_FREE;
2888 }
2889
2890 /** hide edge */
graph_edge_hide(GRAPH * g,int e)2891 void graph_edge_hide(
2892 GRAPH* g, /**< the graph */
2893 int e /**< the edge */
2894 )
2895 {
2896 assert(g != NULL);
2897 assert(e >= 0);
2898 assert(e < g->edges);
2899
2900 /* Immer mit der ersten von beiden Anfangen
2901 */
2902 e -= e % 2;
2903
2904 assert(g->head[e] == g->tail[e + 1]);
2905 assert(g->tail[e] == g->head[e + 1]);
2906
2907 g->grad[g->head[e]]--;
2908 g->grad[g->tail[e]]--;
2909
2910 removeEdge(g, e);
2911
2912 assert(g->ieat[e] != EAT_FREE);
2913 assert(g->ieat[e] != EAT_HIDE);
2914 assert(g->oeat[e] != EAT_FREE);
2915 assert(g->oeat[e] != EAT_HIDE);
2916
2917 g->ieat[e] = EAT_HIDE;
2918 g->oeat[e] = EAT_HIDE;
2919
2920 e++;
2921
2922 removeEdge(g, e);
2923
2924 assert(g->ieat[e] != EAT_FREE);
2925 assert(g->ieat[e] != EAT_HIDE);
2926 assert(g->oeat[e] != EAT_FREE);
2927 assert(g->oeat[e] != EAT_HIDE);
2928
2929 g->ieat[e] = EAT_HIDE;
2930 g->oeat[e] = EAT_HIDE;
2931 }
2932
2933
2934 /** print edge info */
graph_edge_printInfo(SCIP * scip,const GRAPH * g,int e)2935 void graph_edge_printInfo(
2936 SCIP* scip, /**< SCIP data structure */
2937 const GRAPH* g, /**< the graph */
2938 int e /**< the edge */
2939 )
2940 {
2941 const int t = g->tail[e];
2942 const int h = g->head[e];
2943 printf("e: %d %d->%d (%d->%d) \n", e, t, h, g->term[t], g->term[h]);
2944 }
2945
2946 /** changes solution according to given root */
graph_sol_reroot(SCIP * scip,GRAPH * g,int * result,int newroot)2947 SCIP_RETCODE graph_sol_reroot(
2948 SCIP* scip, /**< SCIP data structure */
2949 GRAPH* g, /**< the graph */
2950 int* result, /**< solution array (CONNECT/UNKNOWN) */
2951 int newroot /**< the new root */
2952 )
2953 {
2954 int* queue;
2955 int* const gmark = g->mark;
2956 int size;
2957 const int nnodes = g->knots;
2958
2959 assert(scip != NULL);
2960 assert(g != NULL);
2961 assert(result != NULL);
2962 assert(Is_term(g->term[newroot]));
2963
2964 if( g->grad[newroot] == 0 )
2965 return SCIP_OKAY;
2966
2967 for( int k = 0; k < nnodes; k++ )
2968 gmark[k] = FALSE;
2969
2970 SCIP_CALL( SCIPallocBufferArray(scip, &queue, nnodes) );
2971
2972 gmark[newroot] = TRUE;
2973 size = 0;
2974 queue[size++] = newroot;
2975
2976 /* BFS loop */
2977 while( size )
2978 {
2979 const int node = queue[--size];
2980
2981 /* traverse outgoing arcs */
2982 for( int a = g->outbeg[node]; a != EAT_LAST; a = g->oeat[a] )
2983 {
2984 const int head = g->head[a];
2985
2986 if( !gmark[head] && (result[a] == CONNECT || result[flipedge(a)] == CONNECT ) )
2987 {
2988 if( result[flipedge(a)] == CONNECT )
2989 {
2990 result[a] = CONNECT;
2991 result[flipedge(a)] = UNKNOWN;
2992 }
2993 gmark[head] = TRUE;
2994 queue[size++] = head;
2995 }
2996 }
2997 }
2998
2999 SCIPfreeBufferArray(scip, &queue);
3000
3001 /* adjust solution if infeasible */
3002 for( int k = 0; k < nnodes; k++ )
3003 {
3004 if( !gmark[k] )
3005 {
3006 for( int a = g->outbeg[k]; a != EAT_LAST; a = g->oeat[a] )
3007 {
3008 result[a] = UNKNOWN;
3009 result[flipedge(a)] = UNKNOWN;
3010 }
3011
3012 /* not yet connected terminal? */
3013 if( Is_term(g->term[k]) )
3014 {
3015 int a;
3016 assert(g->stp_type != STP_SPG);
3017
3018 for( a = g->inpbeg[k]; a != EAT_LAST; a = g->ieat[a] )
3019 {
3020 const int node = g->tail[a];
3021 if( gmark[node] && node != newroot )
3022 {
3023 result[a] = CONNECT;
3024 break;
3025 }
3026 }
3027 if( a == EAT_LAST )
3028 {
3029 for( a = g->inpbeg[k]; a != EAT_LAST; a = g->ieat[a] )
3030 {
3031 const int node = g->tail[a];
3032 if( node == newroot )
3033 {
3034 result[a] = CONNECT;
3035 break;
3036 }
3037 }
3038 }
3039 else
3040 gmark[k] = TRUE;
3041 }
3042 }
3043 }
3044
3045 return SCIP_OKAY;
3046 }
3047
3048
3049 /** checks whether edge(s) of given primal solution have been deleted */
graph_sol_unreduced(SCIP * scip,const GRAPH * graph,const int * result)3050 SCIP_Bool graph_sol_unreduced(
3051 SCIP* scip, /**< SCIP data structure */
3052 const GRAPH* graph, /**< graph data structure */
3053 const int* result /**< solution array, indicating whether an edge is in the solution */
3054 )
3055 {
3056 const int nedges = graph->edges;
3057
3058 assert(scip != NULL);
3059 assert(graph != NULL);
3060 assert(result != NULL);
3061
3062 for( int i = 0; i < nedges; i++ )
3063 if( result[i] == CONNECT && graph->oeat[i] == EAT_FREE )
3064 return FALSE;
3065
3066 return TRUE;
3067 }
3068
3069 /** verifies whether a given primal solution is feasible */
graph_sol_valid(SCIP * scip,const GRAPH * graph,const int * result)3070 SCIP_Bool graph_sol_valid(
3071 SCIP* scip, /**< SCIP data structure */
3072 const GRAPH* graph, /**< graph data structure */
3073 const int* result /**< solution array, indicating whether an edge is in the solution */
3074 )
3075 {
3076 int* queue;
3077 STP_Bool* reached;
3078 int root;
3079 int size;
3080 int nnodes;
3081 int termcount;
3082 SCIP_Bool usepterms;
3083
3084 assert(scip != NULL);
3085 assert(graph != NULL);
3086 assert(result != NULL);
3087
3088 reached = NULL;
3089 nnodes = graph->knots;
3090 root = graph->source;
3091 assert(root >= 0);
3092
3093 SCIP_CALL_ABORT( SCIPallocBufferArray(scip, &reached, nnodes) );
3094 SCIP_CALL_ABORT( SCIPallocBufferArray(scip, &queue, nnodes) );
3095
3096 if( (graph->stp_type == STP_MWCSP || graph->stp_type == STP_PCSPG) && !graph->extended )
3097 usepterms = TRUE;
3098 else
3099 usepterms = FALSE;
3100
3101 assert(reached != NULL);
3102
3103 for( int i = 0; i < nnodes; i++ )
3104 reached[i] = FALSE;
3105
3106 /* BFS until all terminals are reached */
3107
3108 termcount = 1;
3109 size = 0;
3110 reached[root] = TRUE;
3111 queue[size++] = root;
3112
3113 while( size )
3114 {
3115 const int node = queue[--size];
3116
3117 for( int e = graph->outbeg[node]; e != EAT_LAST; e = graph->oeat[e] )
3118 {
3119 if( result[e] == CONNECT )
3120 {
3121 const int i = graph->head[e];
3122
3123 /* cycle? */
3124 if( reached[i] )
3125 {
3126 SCIPfreeBufferArray(scip, &queue);
3127 SCIPfreeBufferArray(scip, &reached);
3128 return FALSE;
3129 }
3130
3131 if( usepterms)
3132 {
3133 if( Is_pterm(graph->term[i]) )
3134 termcount++;
3135 }
3136 else
3137 {
3138 if( Is_term(graph->term[i]) )
3139 termcount++;
3140 }
3141
3142 reached[i] = TRUE;
3143 queue[size++] = i;
3144 }
3145 }
3146 }
3147
3148 #if 0
3149 if(termcount != graph->terms)
3150 {
3151 printf("termcount %d graph->terms %d \n", termcount, graph->terms);
3152 printf("root %d \n", root);
3153
3154 for( int i = 0; i < nnodes && 0; i++ )
3155 {
3156 if( Is_term(graph->term[i]) && !reached[i] )
3157 {
3158 printf("fail %d grad %d\n", i, graph->grad[i]);
3159 for( int e = graph->inpbeg[i]; e != EAT_LAST; e = graph->ieat[e] )
3160 {
3161 printf("tail %d %d \n", graph->tail[e], graph->term[graph->tail[e]]);
3162 }
3163 }
3164 }
3165 }
3166 #endif
3167 SCIPfreeBufferArray(scip, &queue);
3168 SCIPfreeBufferArray(scip, &reached);
3169
3170 return (termcount == graph->terms);
3171 }
3172
3173 /** mark endpoints of edges in given list */
graph_sol_setNodeList(const GRAPH * g,STP_Bool * solnode,IDX * listnode)3174 void graph_sol_setNodeList(
3175 const GRAPH* g, /**< graph data structure */
3176 STP_Bool* solnode, /**< solution nodes array (TRUE/FALSE) */
3177 IDX* listnode /**< edge list */
3178 )
3179 {
3180 int i;
3181 IDX* curr;
3182
3183 assert(g != NULL);
3184 assert(solnode != NULL);
3185
3186 curr = listnode;
3187
3188 while( curr != NULL )
3189 {
3190 i = curr->index;
3191
3192 solnode[g->head[i]] = TRUE;
3193 solnode[g->tail[i]] = TRUE;
3194
3195 curr = curr->parent;
3196 }
3197 }
3198
3199 /** compute solution value for given edge-solution array (CONNECT/UNKNOWN) and offset */
graph_sol_getObj(const SCIP_Real * edgecost,const int * soledge,SCIP_Real offset,int nedges)3200 SCIP_Real graph_sol_getObj(
3201 const SCIP_Real* edgecost,
3202 const int* soledge,
3203 SCIP_Real offset,
3204 int nedges
3205 )
3206 {
3207 SCIP_Real obj = offset;
3208 int e;
3209
3210 for( e = 0; e < nedges; e++ )
3211 if( soledge[e] == CONNECT )
3212 obj += edgecost[e];
3213
3214 return obj;
3215 }
3216
3217 /** get original solution */
graph_sol_getOrg(SCIP * scip,const GRAPH * transgraph,const GRAPH * orggraph,const int * transsoledge,int * orgsoledge)3218 SCIP_RETCODE graph_sol_getOrg(
3219 SCIP* scip, /**< SCIP data structure */
3220 const GRAPH* transgraph, /**< the transformed graph */
3221 const GRAPH* orggraph, /**< the original graph */
3222 const int* transsoledge, /**< solution for transformed problem */
3223 int* orgsoledge /**< new retransformed solution */
3224 )
3225 {
3226 STP_Bool* orgnodearr;
3227 STP_Bool* transnodearr = NULL;
3228
3229 IDX** const ancestors = transgraph->ancestors;
3230
3231 const int transnedges = transgraph->edges;
3232 const int transnnodes = transgraph->knots;
3233 const int orgnnodes = orggraph->knots;
3234 const SCIP_Bool pcmw = graph_pc_isPcMw(transgraph);
3235
3236 assert(transgraph != NULL && orggraph != NULL && transsoledge != NULL && orgsoledge != NULL);
3237 assert(transgraph->ancestors != NULL);
3238 assert(transgraph->stp_type == orggraph->stp_type);
3239
3240 SCIP_CALL( SCIPallocBufferArray(scip, &orgnodearr, orgnnodes) );
3241
3242 if( pcmw )
3243 {
3244 SCIP_CALL( SCIPallocBufferArray(scip, &transnodearr, transnnodes) );
3245
3246 for( int k = 0; k < transnnodes; k++ )
3247 transnodearr[k] = FALSE;
3248
3249 for( int e = 0; e < transnedges; e++ )
3250 if( transsoledge[e] == CONNECT )
3251 {
3252 transnodearr[transgraph->tail[e]] = TRUE;
3253 transnodearr[transgraph->head[e]] = TRUE;
3254 }
3255 }
3256
3257 for( int k = 0; k < orgnnodes; k++ )
3258 orgnodearr[k] = FALSE;
3259
3260 for( int e = 0; e < transnedges; e++ )
3261 if( transsoledge[e] == CONNECT )
3262 graph_sol_setNodeList(orggraph, orgnodearr, ancestors[e]);
3263
3264 /* retransform edges fixed during graph reduction */
3265 graph_sol_setNodeList(orggraph, orgnodearr, transgraph->fixedges);
3266
3267 if( pcmw )
3268 {
3269 SCIP_CALL( graph_sol_markPcancestors(scip, transgraph->pcancestors, orggraph->tail, orggraph->head, orgnnodes,
3270 orgnodearr, NULL, NULL, NULL, NULL ) );
3271 }
3272
3273 for( int e = 0; e < orggraph->edges; e++ )
3274 orgsoledge[e] = UNKNOWN;
3275
3276 /* prune solution (in original graph) */
3277 if( pcmw )
3278 SCIP_CALL( SCIPStpHeurTMPrunePc(scip, orggraph, orggraph->cost, orgsoledge, orgnodearr) );
3279 else
3280 SCIP_CALL( SCIPStpHeurTMPrune(scip, orggraph, orggraph->cost, 0, orgsoledge, orgnodearr) );
3281
3282 SCIPfreeBufferArray(scip, &orgnodearr);
3283 SCIPfreeBufferArrayNull(scip, &transnodearr);
3284
3285 assert(graph_sol_valid(scip, orggraph, orgsoledge));
3286
3287 return SCIP_OKAY;
3288 }
3289
3290
3291 /** mark original solution */
graph_sol_markPcancestors(SCIP * scip,IDX ** pcancestors,const int * tails,const int * heads,int orgnnodes,STP_Bool * solnodemark,STP_Bool * soledgemark,int * solnodequeue,int * nsolnodes,int * nsoledges)3292 SCIP_RETCODE graph_sol_markPcancestors(
3293 SCIP* scip, /**< SCIP data structure */
3294 IDX** pcancestors, /**< the ancestors */
3295 const int* tails, /**< tails array */
3296 const int* heads, /**< heads array */
3297 int orgnnodes, /**< original number of nodes */
3298 STP_Bool* solnodemark, /**< solution nodes mark array */
3299 STP_Bool* soledgemark, /**< solution edges mark array or NULL */
3300 int* solnodequeue, /**< solution nodes queue or NULL */
3301 int* nsolnodes, /**< number of solution nodes or NULL */
3302 int* nsoledges /**< number of solution edges or NULL */
3303 )
3304 {
3305 int* queue;
3306 int nnodes;
3307 int nedges = (nsoledges != NULL)? *nsoledges : 0;
3308 int qstart;
3309 int qend;
3310
3311 assert(scip != NULL && tails != NULL && heads != NULL && pcancestors != NULL && solnodemark != NULL);
3312
3313 if( solnodequeue != NULL )
3314 queue = solnodequeue;
3315 else
3316 SCIP_CALL( SCIPallocBufferArray(scip, &queue, orgnnodes) );
3317
3318 if( nsolnodes == NULL )
3319 {
3320 assert(solnodequeue == NULL);
3321 nnodes = 0;
3322 for( int k = 0; k < orgnnodes; k++ )
3323 if( solnodemark[k] )
3324 queue[nnodes++] = k;
3325 }
3326 else
3327 {
3328 nnodes = *nsolnodes;
3329 assert(solnodequeue != NULL);
3330 }
3331
3332 qstart = 0;
3333 qend = nnodes;
3334
3335 while( qend != qstart )
3336 {
3337 int k = qstart;
3338
3339 assert(qstart < qend);
3340 qstart = qend;
3341
3342 for( ; k < qend; k++ )
3343 {
3344 const int ancestornode = queue[k];
3345
3346 assert(solnodemark[ancestornode]);
3347
3348 for( IDX* curr = pcancestors[ancestornode]; curr != NULL; curr = curr->parent )
3349 {
3350 const int ancestoredge = curr->index;
3351 assert(tails[ancestoredge] < orgnnodes && heads[ancestoredge] < orgnnodes);
3352
3353 if( soledgemark != NULL && !soledgemark[ancestoredge] )
3354 {
3355 soledgemark[ancestoredge] = TRUE;
3356 nedges++;
3357 }
3358 if( !solnodemark[tails[ancestoredge]] )
3359 {
3360 solnodemark[tails[ancestoredge]] = TRUE;
3361 queue[nnodes++] = tails[ancestoredge];
3362 }
3363 if( !solnodemark[heads[ancestoredge]] )
3364 {
3365 solnodemark[heads[ancestoredge]] = TRUE;
3366 queue[nnodes++] = heads[ancestoredge];
3367 }
3368 }
3369 }
3370 qend = nnodes;
3371 }
3372
3373 if( nsolnodes != NULL )
3374 *nsolnodes = nnodes;
3375
3376 if( nsoledges != NULL )
3377 *nsoledges = nedges;
3378
3379 if( solnodequeue == NULL )
3380 SCIPfreeBufferArray(scip, &queue);
3381
3382 return SCIP_OKAY;
3383 }
3384
3385 /** get (real) number of nodes , edges, terminals */
graph_get_NVET(const GRAPH * graph,int * nnodes,int * nedges,int * nterms)3386 void graph_get_NVET(
3387 const GRAPH* graph, /**< the graph */
3388 int* nnodes, /**< number of nodes */
3389 int* nedges, /**< number of edges */
3390 int* nterms /**< number of terminals */
3391 )
3392 {
3393 int v = 0;
3394 int e = 0;
3395 int t = 0;
3396 int vorg;
3397
3398 assert(graph != NULL);
3399
3400 vorg = graph->knots;
3401
3402 for( int k = 0; k < vorg; k++ )
3403 {
3404 if( graph->grad[k] > 0 )
3405 {
3406 v++;
3407 e += graph->grad[k];
3408 if( Is_term(graph->term[k]) )
3409 t++;
3410 }
3411 }
3412
3413 *nnodes = v;
3414 *nedges = e;
3415 *nterms = t;
3416
3417 return;
3418 }
3419
3420 /* get compressed sparse row arrays representing current graph */
graph_get_csr(const GRAPH * g,int * RESTRICT edgearr,int * RESTRICT tailarr,int * RESTRICT start,int * nnewedges)3421 void graph_get_csr(
3422 const GRAPH* g, /**< the graph */
3423 int* RESTRICT edgearr, /**< original edge array [0,...,nedges - 1] */
3424 int* RESTRICT tailarr, /**< tail of csr edge [0,...,nedges - 1] */
3425 int* RESTRICT start, /**< start array [0,...,nnodes] */
3426 int* nnewedges /**< pointer to store number of new edges */
3427 )
3428 {
3429 int i = 0;
3430 const int nnodes = g->knots;
3431
3432 assert(g != NULL);
3433 assert(tailarr != NULL);
3434 assert(edgearr != NULL);
3435 assert(start != NULL);
3436
3437 for( int k = 0; k < nnodes; k++ )
3438 {
3439 start[k] = i;
3440 for( int e = g->inpbeg[k]; e != EAT_LAST; e = g->ieat[e] )
3441 {
3442 edgearr[i] = e;
3443 tailarr[i++] = g->tail[e] + 1;
3444 }
3445 }
3446
3447 *nnewedges = i;
3448 start[nnodes] = i;
3449 }
3450
3451 /* gets edge conflicts */
graph_get_edgeConflicts(SCIP * scip,const GRAPH * g)3452 SCIP_RETCODE graph_get_edgeConflicts(
3453 SCIP* scip, /**< SCIP data structure */
3454 const GRAPH* g /**< the graph */
3455 )
3456 {
3457 int* childcount;
3458 int nconflicts;
3459 const int nedges = g->edges;
3460 const int nedgesorg = g->orgedges;
3461
3462 assert(scip != NULL && g != NULL);
3463 assert(g->ancestors != NULL);
3464 assert(nedgesorg % 2 == 0);
3465
3466 printf("orgedes %d \n", nedgesorg);
3467
3468 SCIP_CALL( SCIPallocBufferArray(scip, &childcount, nedgesorg / 2) );
3469
3470 for( int e = 0; e < nedgesorg / 2; e++ )
3471 childcount[e] = 0;
3472
3473 for( int e = 0; e < nedges; e += 2 )
3474 for( IDX* curr = g->ancestors[e]; curr != NULL; curr = curr->parent )
3475 {
3476 assert(curr->index >= 0 && curr->index / 2 < nedgesorg / 2);
3477 childcount[curr->index / 2]++;
3478 }
3479
3480 nconflicts = 0;
3481
3482 for( int e = 0; e < nedgesorg / 2; e++ )
3483 if( childcount[e] > 1 )
3484 nconflicts++;
3485
3486 printf("nconflicts %d \n", nconflicts);
3487
3488 SCIPfreeBufferArray(scip, &childcount);
3489
3490 return SCIP_OKAY;
3491 }
3492
3493
3494 /** initialize graph */
graph_init(SCIP * scip,GRAPH ** g,int ksize,int esize,int layers)3495 SCIP_RETCODE graph_init(
3496 SCIP* scip, /**< SCIP data structure */
3497 GRAPH** g, /**< new graph */
3498 int ksize, /**< slots for nodes */
3499 int esize, /**< slots for edges */
3500 int layers /**< number of layers (only needed for packing, otherwise 1) */
3501 )
3502 {
3503 GRAPH* p;
3504
3505 assert(ksize > 0);
3506 assert(ksize < INT_MAX);
3507 assert(esize >= 0);
3508 assert(esize < INT_MAX);
3509 assert(layers > 0);
3510 assert(layers < SHRT_MAX);
3511
3512 SCIP_CALL( SCIPallocMemory(scip, g) );
3513 p = *g;
3514 assert(p != NULL);
3515
3516 /* ancestor data for retransformation after reductions */
3517 p->fixedges = NULL;
3518 p->ancestors = NULL;
3519 p->pcancestors = NULL;
3520 p->orgtail = NULL;
3521 p->orghead = NULL;
3522 p->rootedgeprevs = NULL;
3523 p->norgmodelknots = 0;
3524 p->norgmodeledges = 0;
3525 p->ksize = ksize;
3526 p->orgknots = 0;
3527 p->orgedges = 0;
3528 p->knots = 0;
3529 p->terms = 0;
3530 p->orgsource = UNKNOWN;
3531 p->stp_type = UNKNOWN;
3532 p->layers = layers;
3533 p->hoplimit = UNKNOWN;
3534 p->extended = FALSE;
3535 p->source = -1;
3536
3537 SCIP_CALL( SCIPallocMemoryArray(scip, &(p->term), ksize) );
3538 SCIP_CALL( SCIPallocMemoryArray(scip, &(p->mark), ksize) );
3539 SCIP_CALL( SCIPallocMemoryArray(scip, &(p->grad), ksize) );
3540 SCIP_CALL( SCIPallocMemoryArray(scip, &(p->inpbeg), ksize) );
3541 SCIP_CALL( SCIPallocMemoryArray(scip, &(p->outbeg), ksize) );
3542 SCIP_CALL( SCIPallocMemoryArray(scip, &(p->cost), esize) );
3543 SCIP_CALL( SCIPallocMemoryArray(scip, &(p->tail), esize) );
3544 SCIP_CALL( SCIPallocMemoryArray(scip, &(p->head), esize) );
3545 SCIP_CALL( SCIPallocMemoryArray(scip, &(p->ieat), esize) );
3546 SCIP_CALL( SCIPallocMemoryArray(scip, &(p->oeat), esize) );
3547
3548 p->esize = esize;
3549 p->edges = 0;
3550 p->prize = NULL;
3551 p->maxdeg = NULL;
3552 p->grid_coordinates = NULL;
3553 p->grid_ncoords = NULL;
3554 p->mincut_dist = NULL;
3555 p->mincut_head = NULL;
3556 p->mincut_numb = NULL;
3557 p->mincut_prev = NULL;
3558 p->mincut_next = NULL;
3559 p->mincut_temp = NULL;
3560 p->mincut_e = NULL;
3561 p->mincut_x = NULL;
3562 p->mincut_r = NULL;
3563 p->path_heap = NULL;
3564 p->path_state = NULL;
3565 p->term2edge = NULL;
3566
3567 SCIPdebugMessage("Initialized new graph \n");
3568
3569 return SCIP_OKAY;
3570 }
3571
3572 /** initialize data structures required to keep track of reductions */
graph_init_history(SCIP * scip,GRAPH * graph)3573 SCIP_RETCODE graph_init_history(
3574 SCIP* scip, /**< SCIP data structure */
3575 GRAPH* graph /**< graph */
3576 )
3577 {
3578 IDX** ancestors; /* ancestor lists array (over all edges) */
3579 IDX** pcancestors; /* ancestor lists array (over all nodes) */
3580 int* tail; /* tail of all edges */
3581 int* head; /* head of all edges */
3582 int* orgtail; /* (original) tail of all original edges */
3583 int* orghead; /* (original) head of all original edges */
3584 int nedges;
3585 SCIP_Bool pcmw;
3586
3587 assert(scip != NULL);
3588 assert(graph != NULL);
3589
3590 pcmw = graph_pc_isPcMw(graph);
3591
3592 nedges = graph->edges;
3593
3594 SCIP_CALL( SCIPallocMemoryArray(scip, &(graph->orgtail), nedges) );
3595 SCIP_CALL( SCIPallocMemoryArray(scip, &(graph->orghead), nedges) );
3596
3597 tail = graph->tail;
3598 head = graph->head;
3599 orgtail = graph->orgtail;
3600 orghead = graph->orghead;
3601
3602 for( int e = 0; e < nedges; e++ )
3603 {
3604 orgtail[e] = tail[e];
3605 orghead[e] = head[e];
3606 }
3607
3608 if( pcmw )
3609 {
3610 const int nnodes = graph->knots;
3611
3612 SCIP_CALL( SCIPallocMemoryArray(scip, &(graph->pcancestors), nnodes) );
3613
3614 pcancestors = graph->pcancestors;
3615
3616 for( int k = 0; k < nnodes; k++ )
3617 pcancestors[k] = NULL;
3618 }
3619
3620 SCIP_CALL( SCIPallocMemoryArray(scip, &(graph->ancestors), nedges) );
3621
3622 ancestors = graph->ancestors;
3623
3624 for( int e = 0; e < nedges; e++ )
3625 {
3626 SCIP_CALL( SCIPallocBlockMemory(scip, &(ancestors[e])) ); /*lint !e866*/
3627 (ancestors)[e]->index = e;
3628 (ancestors)[e]->parent = NULL;
3629 }
3630
3631 return SCIP_OKAY;
3632 }
3633
3634 /** enlarge given graph */
graph_resize(SCIP * scip,GRAPH * g,int ksize,int esize,int layers)3635 SCIP_RETCODE graph_resize(
3636 SCIP* scip, /**< SCIP data structure */
3637 GRAPH* g, /**< graph to be resized */
3638 int ksize, /**< new node slots */
3639 int esize, /**< new edge slots */
3640 int layers /**< layers (set to -1 by default) */
3641 )
3642 {
3643 assert(scip != NULL);
3644 assert(g != NULL);
3645 assert((ksize < 0) || (ksize >= g->knots));
3646 assert((esize < 0) || (esize >= g->edges));
3647 assert((layers < 0) || (layers >= g->layers));
3648
3649 if( (layers > 0) && (layers != g->layers) )
3650 g->layers = layers;
3651
3652 if( (ksize > 0) && (ksize != g->ksize) )
3653 {
3654 SCIP_CALL( SCIPreallocMemoryArray(scip, &(g->term), ksize) );
3655 SCIP_CALL( SCIPreallocMemoryArray(scip, &(g->mark), ksize) );
3656 SCIP_CALL( SCIPreallocMemoryArray(scip, &(g->grad), ksize) );
3657 SCIP_CALL( SCIPreallocMemoryArray(scip, &(g->inpbeg), ksize) );
3658 SCIP_CALL( SCIPreallocMemoryArray(scip, &(g->outbeg), ksize) );
3659
3660 g->ksize = ksize;
3661 }
3662 if( (esize > 0) && (esize != g->esize) )
3663 {
3664 SCIP_CALL( SCIPreallocMemoryArray(scip, &(g->cost), esize) );
3665 SCIP_CALL( SCIPreallocMemoryArray(scip, &(g->tail), esize) );
3666 SCIP_CALL( SCIPreallocMemoryArray(scip, &(g->head), esize) );
3667 SCIP_CALL( SCIPreallocMemoryArray(scip, &(g->ieat), esize) );
3668 SCIP_CALL( SCIPreallocMemoryArray(scip, &(g->oeat), esize) );
3669
3670 g->esize = esize;
3671 }
3672
3673 return SCIP_OKAY;
3674 }
3675
3676
3677 /** free the graph */
graph_free(SCIP * scip,GRAPH ** graph,SCIP_Bool final)3678 void graph_free(
3679 SCIP* scip, /**< SCIP data structure */
3680 GRAPH** graph, /**< graph to be freed */
3681 SCIP_Bool final /**< delete ancestor data structures? */
3682 )
3683 {
3684 GRAPH* p;
3685
3686 assert(scip != NULL);
3687 assert(graph != NULL);
3688
3689 p = *graph;
3690 assert(p != NULL);
3691
3692 graph_free_history(scip, p);
3693
3694 if( final )
3695 graph_free_historyDeep(scip, p);
3696
3697 if( p->prize != NULL )
3698 {
3699 assert(p->term2edge != NULL);
3700 SCIPfreeMemoryArray(scip, &(p->term2edge));
3701 SCIPfreeMemoryArray(scip, &(p->prize));
3702 }
3703
3704 if( p->stp_type == STP_DCSTP )
3705 {
3706 SCIPfreeMemoryArray(scip, &(p->maxdeg));
3707 }
3708 else if( p->stp_type == STP_RSMT )
3709 {
3710 if( p->grid_coordinates != NULL )
3711 {
3712 assert(p->grid_coordinates != NULL);
3713 for( int i = p->grid_dim - 1; i >= 0; i-- )
3714 SCIPfreeMemoryArray(scip, &(p->grid_coordinates[i]));
3715
3716 SCIPfreeMemoryArray(scip, &(p->grid_coordinates));
3717 }
3718
3719 if( p->grid_ncoords != NULL )
3720 SCIPfreeMemoryArray(scip, &(p->grid_ncoords));
3721 }
3722
3723 SCIPfreeMemoryArray(scip, &(p->oeat));
3724 SCIPfreeMemoryArray(scip, &(p->ieat));
3725 SCIPfreeMemoryArray(scip, &(p->head));
3726 SCIPfreeMemoryArray(scip, &(p->tail));
3727 SCIPfreeMemoryArray(scip, &(p->cost));
3728 SCIPfreeMemoryArray(scip, &(p->outbeg));
3729 SCIPfreeMemoryArray(scip, &(p->inpbeg));
3730 SCIPfreeMemoryArray(scip, &(p->grad));
3731 SCIPfreeMemoryArray(scip, &(p->mark));
3732 SCIPfreeMemoryArray(scip, &(p->term));
3733 SCIPfreeMemoryArrayNull(scip, &(p->rootedgeprevs));
3734
3735 SCIPfreeMemory(scip, graph);
3736 }
3737
3738
3739 /** free the history */
graph_free_history(SCIP * scip,GRAPH * p)3740 void graph_free_history(
3741 SCIP* scip, /**< SCIP data */
3742 GRAPH* p /**< graph data */
3743 )
3744 {
3745 if( p->ancestors != NULL )
3746 {
3747 const int nedges = p->edges;
3748
3749 for( int e = nedges - 1; e >= 0; e-- )
3750 {
3751 IDX* curr = p->ancestors[e];
3752 while( curr != NULL )
3753 {
3754 p->ancestors[e] = curr->parent;
3755 SCIPfreeBlockMemory(scip, &(curr));
3756 curr = p->ancestors[e];
3757 }
3758 }
3759 SCIPfreeMemoryArray(scip, &(p->ancestors));
3760 }
3761 }
3762
3763 /** free the deep history */
graph_free_historyDeep(SCIP * scip,GRAPH * p)3764 void graph_free_historyDeep(
3765 SCIP* scip, /**< SCIP data */
3766 GRAPH* p /**< graph data */
3767 )
3768 {
3769 IDX* curr;
3770
3771 assert(scip != NULL);
3772 assert(p != NULL);
3773 assert(p->path_heap == NULL);
3774 assert(p->path_state == NULL);
3775
3776 if( p->pcancestors != NULL )
3777 {
3778 for( int e = p->norgmodelknots - 1; e >= 0; e-- )
3779 {
3780 curr = p->pcancestors[e];
3781 while( curr != NULL )
3782 {
3783 p->pcancestors[e] = curr->parent;
3784 SCIPfreeBlockMemory(scip, &(curr));
3785 curr = p->pcancestors[e];
3786 }
3787 }
3788 SCIPfreeMemoryArray(scip, &(p->pcancestors));
3789 }
3790
3791 if( p->orgtail != NULL )
3792 {
3793 assert(p->orghead != NULL);
3794
3795 SCIPfreeMemoryArray(scip, &(p->orghead));
3796 SCIPfreeMemoryArray(scip, &(p->orgtail));
3797 }
3798 curr = p->fixedges;
3799 while( curr != NULL )
3800 {
3801 p->fixedges = curr->parent;
3802 SCIPfreeBlockMemory(scip, &(curr));
3803
3804 curr = p->fixedges;
3805 }
3806 }
3807
3808 /** copy the data of the graph */
graph_copy_data(SCIP * scip,const GRAPH * orgraph,GRAPH * copygraph)3809 SCIP_RETCODE graph_copy_data(
3810 SCIP* scip, /**< SCIP data structure */
3811 const GRAPH* orgraph, /**< original graph */
3812 GRAPH* copygraph /**< graph to be copied to */
3813 )
3814 {
3815 GRAPH* g = copygraph;
3816 const GRAPH* p = orgraph;
3817 const int ksize = p->ksize;
3818 const int esize = p->esize;
3819
3820 assert(scip != NULL);
3821 assert(orgraph != NULL);
3822 assert(copygraph != NULL);
3823 assert(ksize == g->ksize && ksize > 0);
3824 assert(esize == g->esize && esize >= 0);
3825
3826 g->norgmodeledges = p->norgmodeledges;
3827 g->norgmodelknots = p->norgmodelknots;
3828 g->knots = p->knots;
3829 g->terms = p->terms;
3830 g->edges = p->edges;
3831 g->source = p->source;
3832 g->orgsource = p->orgsource;
3833 g->orgedges = p->orgedges;
3834 g->orgknots = p->orgknots;
3835 g->grid_dim = p->grid_dim;
3836 g->stp_type = p->stp_type;
3837 g->hoplimit = p->hoplimit;
3838 g->extended = p->extended;
3839 g->term2edge = NULL;
3840 g->prize = NULL;
3841
3842 BMScopyMemoryArray(g->term, p->term, ksize);
3843 BMScopyMemoryArray(g->mark, p->mark, ksize);
3844 BMScopyMemoryArray(g->grad, p->grad, ksize);
3845 BMScopyMemoryArray(g->inpbeg, p->inpbeg, ksize);
3846 BMScopyMemoryArray(g->outbeg, p->outbeg, ksize);
3847 BMScopyMemoryArray(g->cost, p->cost, esize);
3848 BMScopyMemoryArray(g->tail, p->tail, esize);
3849 BMScopyMemoryArray(g->head, p->head, esize);
3850 BMScopyMemoryArray(g->ieat, p->ieat, esize);
3851 BMScopyMemoryArray(g->oeat, p->oeat, esize);
3852
3853 if( g->stp_type == STP_PCSPG || g->stp_type == STP_RPCSPG || g->stp_type == STP_MWCSP || g->stp_type == STP_RMWCSP )
3854 {
3855 SCIP_CALL(SCIPallocMemoryArray(scip, &(g->prize), g->knots));
3856 SCIP_CALL(SCIPallocMemoryArray(scip, &(g->term2edge), g->knots));
3857
3858 for( int k = 0; k < g->knots; k++ )
3859 if( Is_term(p->term[k]) )
3860 g->prize[k] = 0.0;
3861 else
3862 g->prize[k] = p->prize[k];
3863
3864 assert(p->term2edge != NULL);
3865
3866 BMScopyMemoryArray(g->term2edge, p->term2edge, g->knots);
3867 }
3868 else if( g->stp_type == STP_DCSTP )
3869 {
3870 assert(p->maxdeg != NULL);
3871
3872 SCIP_CALL(SCIPallocMemoryArray(scip, &(g->maxdeg), g->knots));
3873
3874 for( int k = 0; k < g->knots; k++ )
3875 g->maxdeg[k] = p->maxdeg[k];
3876 }
3877 else if( p->stp_type == STP_RSMT )
3878 {
3879 assert(p->grid_ncoords != NULL);
3880 assert(p->grid_coordinates != NULL);
3881
3882 SCIP_CALL(SCIPallocMemoryArray(scip, &(g->grid_coordinates), p->grid_dim));
3883
3884 BMScopyMemoryArray(g->grid_coordinates, p->grid_coordinates, p->grid_dim);
3885 for( int k = 0; k < p->grid_dim; k++ )
3886 {
3887 SCIP_CALL(SCIPallocMemoryArray(scip, &(g->grid_coordinates[k]), p->terms)); /*lint !e866*/
3888 BMScopyMemoryArray(g->grid_coordinates[k], p->grid_coordinates[k], p->terms); /*lint !e866*/
3889 }
3890 SCIP_CALL(SCIPallocMemoryArray(scip, &(g->grid_ncoords), p->grid_dim));
3891
3892 BMScopyMemoryArray(g->grid_ncoords, p->grid_ncoords, p->grid_dim);
3893 }
3894 assert(graph_valid(g));
3895
3896 return SCIP_OKAY;
3897 }
3898
3899 /** copy the graph */
graph_copy(SCIP * scip,const GRAPH * orgraph,GRAPH ** copygraph)3900 SCIP_RETCODE graph_copy(
3901 SCIP* scip, /**< SCIP data structure */
3902 const GRAPH* orgraph, /**< original graph */
3903 GRAPH** copygraph /**< graph to be created */
3904 )
3905 {
3906 const GRAPH* p = orgraph;
3907 assert(p != NULL);
3908
3909 SCIP_CALL( graph_init(scip, copygraph, p->ksize, p->esize, p->layers) );
3910
3911 SCIP_CALL( graph_copy_data(scip, orgraph, *copygraph) );
3912
3913 return SCIP_OKAY;
3914 }
3915
graph_show(const GRAPH * p)3916 void graph_show(
3917 const GRAPH* p /**< the graph */
3918 )
3919 {
3920 int i;
3921
3922 assert(p != NULL);
3923
3924 for(i = 0; i < p->knots; i++)
3925 if (p->grad[i] > 0)
3926 (void)printf("Knot %d, term=%d, grad=%d, inpbeg=%d, outbeg=%d\n",
3927 i, p->term[i], p->grad[i], p->inpbeg[i], p->outbeg[i]);
3928
3929 (void)fputc('\n', stdout);
3930
3931 for(i = 0; i < p->edges; i++)
3932 if (p->ieat[i] != EAT_FREE)
3933 (void)printf("Edge %d, cost=%g, tail=%d, head=%d, ieat=%d, oeat=%d\n",
3934 i, p->cost[i], p->tail[i], p->head[i], p->ieat[i], p->oeat[i]);
3935
3936 (void)fputc('\n', stdout);
3937 }
3938
3939
3940 /** reinsert all hidden edges */
graph_uncover(GRAPH * g)3941 void graph_uncover(
3942 GRAPH* g /**< the graph */
3943 )
3944 {/*lint --e{850}*/
3945 int head;
3946 int tail;
3947 int e;
3948
3949 assert(g != NULL);
3950
3951 for( e = 0; e < g->edges; e++ )
3952 {
3953 if( g->ieat[e] == EAT_HIDE )
3954 {
3955 assert(e % 2 == 0);
3956 assert(g->oeat[e] == EAT_HIDE);
3957
3958 head = g->head[e];
3959 tail = g->tail[e];
3960
3961 g->grad[head]++;
3962 g->grad[tail]++;
3963
3964 g->ieat[e] = g->inpbeg[head];
3965 g->oeat[e] = g->outbeg[tail];
3966 g->inpbeg[head] = e;
3967 g->outbeg[tail] = e;
3968
3969 e++;
3970
3971 assert(g->ieat[e] == EAT_HIDE);
3972 assert(g->oeat[e] == EAT_HIDE);
3973 assert(g->head[e] == tail);
3974 assert(g->tail[e] == head);
3975
3976 head = g->head[e];
3977 tail = g->tail[e];
3978 g->ieat[e] = g->inpbeg[head];
3979 g->oeat[e] = g->outbeg[tail];
3980 g->inpbeg[head] = e;
3981 g->outbeg[tail] = e;
3982 }
3983 }
3984 }
3985
3986
3987 /** pack the graph, i.e. build a new graph that discards deleted edges and nodes */
graph_pack(SCIP * scip,GRAPH * graph,GRAPH ** newgraph,SCIP_Bool verbose)3988 SCIP_RETCODE graph_pack(
3989 SCIP* scip, /**< SCIP data structure */
3990 GRAPH* graph, /**< the graph */
3991 GRAPH** newgraph, /**< the new graph */
3992 SCIP_Bool verbose /**< verbose? */
3993 )
3994 {
3995 GRAPH* g;
3996 GRAPH* q;
3997 int* new;
3998 int e;
3999 int oldnnodes;
4000 int oldnedges;
4001 int nnodes;
4002 int nedges;
4003 SCIP_Bool rmw;
4004 SCIP_Bool pcmw;
4005
4006 assert(scip != NULL);
4007 assert(graph != NULL);
4008 assert(graph_valid(graph));
4009
4010 g = graph;
4011 nnodes = 0;
4012 nedges = 0;
4013 oldnnodes = g->knots;
4014 oldnedges = g->edges;
4015 SCIP_CALL( SCIPallocBufferArray(scip, &new, oldnnodes) );
4016
4017 if( verbose )
4018 printf("Reduced graph: ");
4019
4020 /* count nodes */
4021 for( int i = 0; i < oldnnodes; i++ )
4022 {
4023 /* are there incident edges to current node? */
4024 if( g->grad[i] > 0 )
4025 new[i] = nnodes++;
4026 else
4027 new[i] = -1;
4028 }
4029
4030 /* graph vanished? */
4031 if( nnodes == 0 )
4032 {
4033 SCIPfreeBufferArray(scip, &new);
4034 new = NULL;
4035 if( verbose )
4036 printf(" graph vanished!\n");
4037
4038 nnodes = 1;
4039 }
4040
4041 /* count edges */
4042 for( int i = 0; i < oldnedges; i++ )
4043 {
4044 if( g->oeat[i] != EAT_FREE )
4045 {
4046 assert(g->ieat[i] != EAT_FREE);
4047 nedges++;
4048 }
4049 }
4050
4051 assert(nnodes > 1 || nedges == 0);
4052 SCIP_CALL( graph_init(scip, newgraph, nnodes, nedges, g->layers) );
4053 q = *newgraph;
4054 q->norgmodelknots = g->norgmodelknots;
4055 q->norgmodeledges = g->norgmodeledges;
4056 q->orgsource = g->orgsource;
4057 q->orgtail = g->orgtail;
4058 q->orghead = g->orghead;
4059 q->orgknots = g->knots;
4060 q->orgedges = g->edges;
4061 q->stp_type = g->stp_type;
4062 q->maxdeg = g->maxdeg;
4063 q->grid_dim = g->grid_dim;
4064 q->grid_ncoords = g->grid_ncoords;
4065 q->grid_coordinates = g->grid_coordinates;
4066 q->fixedges = g->fixedges;
4067 q->hoplimit = g->hoplimit;
4068 q->extended = g->extended;
4069 q->pcancestors = g->pcancestors;
4070
4071 if( new == NULL )
4072 {
4073 q->ancestors = NULL;
4074 graph_free(scip, &g, FALSE);
4075
4076 if( q->stp_type == STP_RSMT )
4077 {
4078 q->grid_ncoords = NULL;
4079 q->grid_coordinates = NULL;
4080 }
4081
4082 graph_knot_add(q, 0);
4083 q->source = 0;
4084 return SCIP_OKAY;
4085 }
4086
4087 SCIP_CALL( SCIPallocMemoryArray(scip, &(q->ancestors), nedges) );
4088
4089 rmw = g->stp_type == STP_RMWCSP;
4090 pcmw = (g->stp_type == STP_MWCSP || g->stp_type == STP_RPCSPG || g->stp_type == STP_PCSPG || g->stp_type == STP_RMWCSP);
4091 if( pcmw )
4092 SCIP_CALL( graph_pc_init(scip, q, nnodes, nnodes) );
4093
4094 /* add nodes (of positive degree) */
4095 if( rmw )
4096 {
4097 int i;
4098 for( i = 0; i < oldnnodes; i++ )
4099 g->mark[i] = (g->grad[i] > 0);
4100
4101 for( e = g->outbeg[g->source]; e != EAT_LAST; e = g->oeat[e] )
4102 {
4103 if( SCIPisGT(scip, g->cost[e], 0.0) && Is_term(g->term[g->head[e]]) )
4104 {
4105 i = g->head[e];
4106 g->mark[i] = FALSE;
4107 assert(g->grad[i] == 2);
4108 }
4109 }
4110 }
4111
4112 for( int i = 0; i < oldnnodes; i++ )
4113 {
4114 assert(g->term[i] < g->layers);
4115 if( g->grad[i] > 0 )
4116 {
4117 if( pcmw )
4118 {
4119 if( !Is_term(g->term[i]) || (rmw && g->mark[i]) )
4120 q->prize[q->knots] = g->prize[i];
4121 else
4122 q->prize[q->knots] = 0.0;
4123 }
4124 graph_knot_add(q, g->term[i]);
4125 }
4126 }
4127
4128 /* add root */
4129 assert(q->term[new[g->source]] == 0);
4130
4131 q->source = new[g->source];
4132
4133 if( g->stp_type == STP_RPCSPG || g->stp_type == STP_RMWCSP )
4134 q->prize[q->source] = FARAWAY;
4135
4136 /* add edges */
4137 for( int i = 0; i < oldnedges; i += 2 )
4138 {
4139 if( g->ieat[i] == EAT_FREE )
4140 {
4141 assert(g->oeat[i] == EAT_FREE);
4142 assert(g->ieat[i + 1] == EAT_FREE);
4143 assert(g->oeat[i + 1] == EAT_FREE);
4144 SCIPintListNodeFree(scip, &(g->ancestors[i]));
4145 SCIPintListNodeFree(scip, &(g->ancestors[i + 1]));
4146 continue;
4147 }
4148
4149 assert(g->oeat[i] != EAT_FREE);
4150 assert(g->ieat[i + 1] != EAT_FREE);
4151 assert(g->oeat[i + 1] != EAT_FREE);
4152 assert(new[g->tail[i]] >= 0);
4153 assert(new[g->head[i]] >= 0);
4154
4155 e = q->edges;
4156
4157 q->ancestors[e] = NULL;
4158 q->ancestors[e + 1] = NULL;
4159 SCIP_CALL( SCIPintListNodeAppendCopy(scip, &(q->ancestors[e]), g->ancestors[i], NULL) );
4160 SCIP_CALL( SCIPintListNodeAppendCopy(scip, &(q->ancestors[e + 1]), g->ancestors[i + 1], NULL) );
4161
4162 assert(new[g->tail[i]] < nnodes && new[g->head[i]] < nnodes);
4163
4164 if( pcmw )
4165 graph_pc_updateTerm2edge(q, g, new[g->tail[i]], new[g->head[i]], g->tail[i], g->head[i]);
4166
4167 graph_edge_add(scip, q, new[g->tail[i]], new[g->head[i]], g->cost[i], g->cost[Edge_anti(i)]);
4168 }
4169
4170 SCIPfreeBufferArray(scip, &new);
4171
4172 if( g->path_heap != NULL )
4173 graph_path_exit(scip, g);
4174
4175 g->stp_type = UNKNOWN;
4176 graph_free(scip, &g, FALSE);
4177
4178 assert(graph_valid(q));
4179
4180 if( verbose )
4181 printf("Nodes: %d Edges: %d Terminals: %d\n", q->knots, q->edges, q->terms);
4182
4183 return SCIP_OKAY;
4184 }
4185
4186
4187 /** traverse the graph and mark all reached nodes (g->mark[i] has to be FALSE for all i) */
graph_trail(const GRAPH * g,int i)4188 void graph_trail(
4189 const GRAPH* g, /**< the new graph */
4190 int i /**< node to start from */
4191 )
4192 {
4193 int* gmark;
4194
4195 assert(g != NULL);
4196 assert(i >= 0);
4197 assert(i < g->knots);
4198
4199 gmark = g->mark;
4200
4201 if( !gmark[i] )
4202 {
4203 SCIP_QUEUE* queue;
4204 int a;
4205 int head;
4206 int node;
4207 int* pnode;
4208
4209 gmark[i] = TRUE;
4210
4211 if( g->grad[i] == 0 )
4212 return;
4213
4214 SCIP_CALL_ABORT( SCIPqueueCreate(&queue, g->knots, 1.1) );
4215 SCIP_CALL_ABORT( SCIPqueueInsert(queue, &i));
4216
4217 /* BFS loop */
4218 while( !SCIPqueueIsEmpty(queue) )
4219 {
4220 pnode = (SCIPqueueRemove(queue));
4221 node = *pnode;
4222
4223 /* traverse outgoing arcs */
4224 for( a = g->outbeg[node]; a != EAT_LAST; a = g->oeat[a] )
4225 {
4226 head = g->head[a];
4227
4228 if( !gmark[head] )
4229 {
4230 gmark[head] = TRUE;
4231 SCIP_CALL_ABORT(SCIPqueueInsert(queue, &(g->head[a])));
4232 }
4233 }
4234 }
4235 SCIPqueueFree(&queue);
4236 }
4237 }
4238
4239
4240 /** traverse the graph and mark all reached nodes (g->mark[i] has to be FALSE for all i) .... uses an array and should be faster
4241 * than graph_trail, but needs a scip */
graph_trail_arr(SCIP * scip,const GRAPH * g,int i)4242 SCIP_RETCODE graph_trail_arr(
4243 SCIP* scip, /**< scip struct */
4244 const GRAPH* g, /**< the new graph */
4245 int i /**< node to start from */
4246 )
4247 {
4248 int* const gmark = g->mark;
4249
4250 assert(scip != NULL);
4251 assert(g != NULL);
4252 assert(i >= 0);
4253 assert(i < g->knots);
4254
4255 if( !gmark[i] )
4256 {
4257 int* stackarr;
4258 int a;
4259 int head;
4260 int node;
4261 int nnodes;
4262 int stacksize;
4263
4264 gmark[i] = TRUE;
4265
4266 if( g->grad[i] == 0 )
4267 return SCIP_OKAY;
4268
4269 nnodes = g->knots;
4270 stacksize = 0;
4271
4272 SCIP_CALL( SCIPallocBufferArray(scip, &stackarr, nnodes) );
4273
4274 stackarr[stacksize++] = i;
4275
4276 /* DFS loop */
4277 while( stacksize != 0 )
4278 {
4279 node = stackarr[--stacksize];
4280
4281 /* traverse outgoing arcs */
4282 for( a = g->outbeg[node]; a != EAT_LAST; a = g->oeat[a] )
4283 {
4284 head = g->head[a];
4285
4286 if( !gmark[head] )
4287 {
4288 gmark[head] = TRUE;
4289 stackarr[stacksize++] = head;
4290 }
4291 }
4292 }
4293 SCIPfreeBufferArray(scip, &stackarr);
4294 }
4295 return SCIP_OKAY;
4296 }
4297
4298 /** checks whether all terminals are reachable from root */
graph_termsReachable(SCIP * scip,const GRAPH * g,SCIP_Bool * reachable)4299 SCIP_RETCODE graph_termsReachable(
4300 SCIP* scip, /**< scip struct */
4301 const GRAPH* g, /**< the new graph */
4302 SCIP_Bool* reachable /**< are they reachable? */
4303 )
4304 {
4305 const int nnodes = g->knots;
4306
4307 assert(g != NULL);
4308 assert(reachable != NULL);
4309
4310 for( int k = 0; k < nnodes; k++ )
4311 g->mark[k] = FALSE;
4312
4313 *reachable = TRUE;
4314
4315 graph_trail_arr(scip, g, g->source);
4316
4317 for( int k = 0; k < nnodes; k++ )
4318 if( Is_term(g->term[k]) && !g->mark[k] )
4319 {
4320 *reachable = FALSE;
4321 break;
4322 }
4323
4324 return SCIP_OKAY;
4325 }
4326
4327 /** is the given graph valid? */
graph_valid(const GRAPH * g)4328 SCIP_Bool graph_valid(
4329 const GRAPH* g /**< the new graph */
4330 )
4331 {
4332 const char* fehler1 = "*** Graph invalid: Head invalid, Knot %d, Edge %d, Tail=%d, Head=%d\n";
4333 const char* fehler2 = "*** Graph invalid: Tail invalid, Knot %d, Edge %d, Tail=%d, Head=%d\n";
4334 const char* fehler3 = "*** Graph invalid: Source invalid, Layer %d, Source %d, Terminal %d\n";
4335 const char* fehler4 = "*** Graph invalid: FREE invalid, Edge %d/%d\n";
4336 const char* fehler5 = "*** Graph invalid: Anti invalid, Edge %d/%d, Tail=%d/%d, Head=%d/%d\n";
4337 const char* fehler6 = "*** Graph invalid: Knot %d with Grad 0 has Edges\n";
4338 const char* fehler7 = "*** Graph invalid: Knot %d not connected\n";
4339 const char* fehler9 = "*** Graph invalid: Wrong Terminal count, count is %d, should be %d\n";
4340
4341 int k;
4342 int e;
4343 int nterms;
4344 int nnodes;
4345 int nedges;
4346
4347 assert(g != NULL);
4348
4349 nterms = g->terms;
4350 nedges = g->edges;
4351 nnodes = g->knots;
4352
4353 for( k = 0; k < nnodes; k++ )
4354 {
4355 if( Is_term(g->term[k]) )
4356 {
4357 nterms--;
4358 }
4359 for( e = g->inpbeg[k]; e != EAT_LAST; e = g->ieat[e] )
4360 if( g->head[e] != k )
4361 break;
4362
4363 if( e != EAT_LAST )
4364 return((void)fprintf(stderr, fehler1, k, e, g->tail[e], g->head[e]), FALSE);
4365
4366 for( e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
4367 if( g->tail[e] != k )
4368 break;
4369
4370 if( e != EAT_LAST )
4371 return((void)fprintf(stderr, fehler2, k, e, g->tail[e], g->head[e]), FALSE);
4372 }
4373 if( nterms != 0 )
4374 return((void)fprintf(stderr, fehler9, g->terms, g->terms - nterms), FALSE);
4375
4376 if( (g->source < 0 )
4377 || (g->source >= g->knots)
4378 || (g->term[g->source] != 0))
4379 return((void)fprintf(stderr, fehler3,
4380 0, g->source, g->term[g->source]), FALSE);
4381
4382 for( e = 0; e < nedges; e += 2 )
4383 {
4384 if( (g->ieat[e] == EAT_FREE) && (g->oeat[e] == EAT_FREE)
4385 && (g->ieat[e + 1] == EAT_FREE) && (g->oeat[e + 1] == EAT_FREE) )
4386 continue;
4387
4388 if( (g->ieat[e] == EAT_FREE) || (g->oeat[e] == EAT_FREE)
4389 || (g->ieat[e + 1] == EAT_FREE) || (g->oeat[e + 1] == EAT_FREE) )
4390 return((void)fprintf(stderr, fehler4, e, e + 1), FALSE);
4391
4392 if( (g->head[e] != g->tail[e + 1]) || (g->tail[e] != g->head[e + 1]) )
4393 return((void)fprintf(stderr, fehler5,
4394 e, e + 1, g->head[e], g->tail[e + 1],
4395 g->tail[e], g->head[e + 1]), FALSE);
4396 }
4397
4398 for( k = 0; k < nnodes; k++ )
4399 g->mark[k] = FALSE;
4400
4401 graph_trail(g, g->source);
4402
4403 for( k = 0; k < nnodes; k++ )
4404 {
4405 if( (g->grad[k] == 0)
4406 && ((g->inpbeg[k] != EAT_LAST) || (g->outbeg[k] != EAT_LAST)) )
4407 return((void)fprintf(stderr, fehler6, k), FALSE);
4408
4409 if( !g->mark[k] && ((g->grad[k] > 0) || (Is_term(g->term[k])))
4410 && g->stp_type != STP_PCSPG && g->stp_type != STP_MWCSP && g->stp_type != STP_RMWCSP )
4411 return((void)fprintf(stderr, fehler7, k), FALSE);
4412 }
4413
4414 if( (g->stp_type == STP_PCSPG || g->stp_type == STP_MWCSP || g->stp_type == STP_RPCSPG || g->stp_type == STP_RMWCSP) )
4415 {
4416 int npterms = 0;
4417 const int root = g->source;
4418 const SCIP_Bool extended = g->extended;
4419 const SCIP_Bool rooted = (g->stp_type == STP_RPCSPG || g->stp_type == STP_RMWCSP);
4420 nterms = 0;
4421
4422 assert(g->prize != NULL);
4423 assert(g->term2edge != NULL);
4424
4425 for( k = 0; k < nnodes; k++ )
4426 {
4427 if( k == root || (rooted && g->term2edge[k] < 0) )
4428 continue;
4429
4430 if( (extended ? Is_term(g->term[k]) : Is_pterm(g->term[k])) )
4431 {
4432 int e2;
4433 int pterm;
4434 const int term = k;
4435 nterms++;
4436
4437 if( g->grad[k] != 2 )
4438 {
4439 SCIPdebugMessage("terminal degree != 2 for %d \n", k);
4440 return FALSE;
4441 }
4442
4443 for( e = g->inpbeg[term]; e != EAT_LAST; e = g->ieat[e] )
4444 if( g->tail[e] == root )
4445 break;
4446
4447 if( e == EAT_LAST )
4448 {
4449 SCIPdebugMessage("no edge to root for term %d \n", term);
4450 return FALSE;
4451 }
4452
4453 for( e2 = g->outbeg[term]; e2 != EAT_LAST; e2 = g->oeat[e2] )
4454 {
4455 pterm = g->head[e2];
4456 if( (extended ? Is_pterm(g->term[pterm]) : Is_term(g->term[pterm])) && pterm != root )
4457 break;
4458 }
4459
4460 if( e2 == EAT_LAST)
4461 {
4462 SCIPdebugMessage("no terminal for dummy %d \n", g->head[e2]);
4463 return FALSE;
4464 }
4465
4466 assert(pterm != root);
4467
4468 if( e2 != g->term2edge[term] )
4469 {
4470 SCIPdebugMessage("term2edge for node %d faulty \n", term);
4471 return FALSE;
4472 }
4473
4474 if( g->cost[e] != g->prize[pterm] )
4475 {
4476 SCIPdebugMessage("prize mismatch for node %d: \n", k);
4477 return FALSE;
4478 }
4479 }
4480 else if( (extended ? Is_pterm(g->term[k]) : Is_term(g->term[k])) )
4481 {
4482 npterms++;
4483 }
4484 }
4485 if( nterms != npterms || nterms != g->terms - 1 )
4486 {
4487 if( !rooted )
4488 {
4489 SCIPdebugMessage("wrong terminal count \n");
4490 return FALSE;
4491 }
4492 }
4493
4494 for( k = 0; k < nnodes; k++ )
4495 {
4496 g->mark[k] = (g->grad[k] > 0);
4497
4498 if( !extended && (Is_pterm(g->term[k]) || k == root) )
4499 g->mark[k] = FALSE;
4500 }
4501 if( !extended && (g->stp_type == STP_RPCSPG || g->stp_type == STP_RMWCSP) )
4502 g->mark[root] = TRUE;
4503
4504 }
4505 else
4506 {
4507 for( k = 0; k < nnodes; k++ )
4508 g->mark[k] = (g->grad[k] > 0);
4509 }
4510
4511 return TRUE;
4512 }
4513