Scippy

SCIP

Solving Constraint Integer Programs

extreduce_dbg.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2022 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scip.zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file extreduce_dbg.c
17  * @brief extended reductions for Steiner tree problems
18  * @author Daniel Rehfeldt
19  *
20  * This file implements extended reduction debugging routines for several Steiner problems.
21  *
22  * A list of all interface methods can be found in reduce.h.
23  *
24  */
25 
26 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
27 // #define SCIP_DEBUG
28 // #define STP_DEBUG_EXTPC // use if special sds for PC are deactivated
29 
30 #include <stdio.h>
31 #include <stdlib.h>
32 #include <assert.h>
33 #include "graph.h"
34 #include "portab.h"
35 #include <limits.h>
36 #include "extreduce.h"
37 
38 
39 /** get SD MST weight
40  * NOTE: might deviate because only getSd is used...maybe only use double in the code? */
41 static
43  SCIP* scip, /**< SCIP */
44  const GRAPH* graph, /**< graph data structure */
45  int nnodes, /**< number of nodes for MST computation */
46  const int nodes[], /**< nodes (from graph) for MST computation */
47  SCIP_Bool computeLower, /**< compute lower bound? (otherwise upper) */
48  EXTDATA* extdata /**< extension data */
49  )
50 {
51  SCIP_Real mstweight;
52  CGRAPH* cgraph;
53  CMST* cmst;
54  SCIP_Real* adjcosts;
55 
58 
59  adjcosts = cgraph->adjedgecosts;
60 
61  assert(adjcosts);
62  assert(nnodes > 0 && nnodes <= STP_EXTTREE_MAXNLEAVES_GUARD);
63 
64  //printf(" \n");
65 
66  /* build the MST graph */
67  for( int i = 0; i < nnodes; i++ )
68  cgraph_node_append(cgraph, i);
69 
70  for( int i = 0; i < nnodes; i++ )
71  {
72  const int startnode = nodes[i];
73 
74  for( int j = 0; j < nnodes; j++ )
75  {
76  SCIP_Real specialDist;
77 
78  if( i == j )
79  {
80  specialDist = FARAWAY;
81  }
82  else
83  {
84  const int endnode = nodes[j];
85 
86  if( computeLower )
87  {
88  specialDist = extreduce_extGetSdDouble(scip, graph, startnode, endnode, extdata);
89  }
90  else
91  {
92  const SCIP_Real sd1 = extreduce_extGetSd(scip, graph, startnode, endnode, extdata);
93  const SCIP_Real sd2 = extreduce_extGetSd(scip, graph, endnode, startnode, extdata);
94 
95  if( EQ(sd1, -1.0) || EQ(sd2, -1.0) )
96  {
97  specialDist = -1.0;
98  }
99  else
100  {
101  specialDist = MAX(sd1, sd2);
102  }
103  }
104 
105  if( specialDist <= 0.0 )
106  {
107  assert(EQ(specialDist, -1.0));
108  specialDist = BLOCKED;
109  }
110 
111  // printf("sd %d->%d: %f \n", startnode, endnode, specialDist);
112  }
113 
114  adjcosts[j] = specialDist;
115  }
116 
117  cgraph_node_applyMinAdjCosts(cgraph, i, i);
118  }
119 
120  /* compute the MST */
121  cmst_computeMst(cgraph, 0, cmst);
122 
123  mstweight = cmst->mstobj;
124 
125  cmst_free(scip, &cmst);
126  cgraph_free(scip, &cgraph);
127 
128  assert(GE(mstweight, 0.0));
129 
130  if( GE(mstweight, BLOCKED) )
131  {
132  mstweight = FARAWAY;
133  }
134 
135  return mstweight;
136 }
137 
138 
139 /** gets MST weight for SD MST spanning all leaves and extension vertex
140  * NOTE: only for debugging! very slow!
141  * */
142 static
144  SCIP* scip, /**< SCIP */
145  const GRAPH* graph, /**< graph data structure */
146  int extvert, /**< extended vertex */
147  SCIP_Bool computeLower, /**< compute lower bound? (otherwise upper) */
148  EXTDATA* extdata /**< extension data */
149 )
150 {
151  SCIP_Real mstweight;
152  int* extleaves;
153  const int* const leaves = extdata->tree_leaves;
154  const int nleaves = extdata->tree_nleaves;
155 
156  assert(extvert >= 0 && extvert < graph->knots);
157 
159 
160  assert(nleaves <= STP_EXTTREE_MAXNLEAVES_GUARD);
161 
162  for( int i = 0; i < nleaves; ++i )
163  {
164  const int leaf = leaves[i];
165  assert(extvert != leaf);
166 
167  extleaves[i] = leaf;
168  }
169 
170  extleaves[nleaves] = extvert;
171 
172  mstweight = sdmstGetWeight(scip, graph, nleaves + 1, extleaves, computeLower, extdata);
173 
174  SCIPfreeBufferArray(scip, &extleaves);
175 
176  return mstweight;
177 }
178 
179 
180 /** gets nodes of top levelbase MST */
181 static
183  int extnode, /**< extension nodes */
184  int nnodes, /**< number of nodes */
185  int nodes[], /**< nodes (from graph) to be filled in */
186  const EXTDATA* extdata /**< extension data */
187 )
188 {
189  const int* const leaves = extdata->tree_leaves;
190  int j = 0;
191  const int nleaves = extdata->tree_nleaves;
192 
193  assert(nnodes == nleaves - 1);
194 
195  SCIPdebugMessage("mstTopLevel nodes: \n");
196 
197  for( int i = 0; i < nleaves; ++i )
198  {
199  const int leaf = leaves[i];
200 
201  if( leaf == extnode )
202  continue;
203 
204  SCIPdebugMessage(" %d \n", leaf);
205 
206  nodes[j] = leaf;
207 
208  j++;
209  }
210 
211  assert(j == nnodes);
212 }
213 
214 
215 /** is the weight of the top levelbase MST valid? */
216 static
218  SCIP* scip, /**< SCIP */
219  const GRAPH* graph, /**< graph data structure */
220  int base_nnodes, /**< number of nodes */
221  int base_nodes[], /**< MST nodes (from graph) */
222  const CSR* mst_base, /**< the stored MST */
223  EXTDATA* extdata /**< extension data */
224 )
225 {
226  const SCIP_Real mstweight_org = reduce_dcmstGetWeight(scip, mst_base);
227  const SCIP_Real mstweight_upper = sdmstGetWeight(scip, graph, base_nnodes, base_nodes, FALSE, extdata);
228  SCIP_Bool isValid = TRUE;
229 
230  assert(base_nnodes == mst_base->nnodes);
231 
232  if( GT(mstweight_org, mstweight_upper) )
233  {
234  SCIPdebugMessage("top baselevel MST weight violates upper bound: %f > %f \n",
235  mstweight_org, mstweight_upper);
236 
237  isValid = FALSE;
238  }
239  else
240  {
241  const SCIP_Real mstweight_lower = sdmstGetWeight(scip, graph, base_nnodes, base_nodes, TRUE, extdata);
242 
243  if( LT(mstweight_org, mstweight_lower) )
244  {
245  SCIPdebugMessage("top baselevel MST weight violates lower bound: %f < %f \n",
246  mstweight_org, mstweight_lower);
247 
248  isValid = FALSE;
249  }
250  }
251 
252  return isValid;
253 }
254 
255 
256 
257 /** returns shortest distance between given vertices */
258 static
260  SCIP* scip, /**< SCIP */
261  const GRAPH* g, /**< graph data structure */
262  int vertexBlocked, /**< forbidden vertex */
263  const DISTDATA* distdata, /**< distance data */
264  int vertex1, /**< first vertex */
265  int vertex2 /**< second vertex */
266  )
267 {
268  SCIP_Real* dist;
269  DHEAP* dheap = NULL;
270  int* state;
271  const SCIP_Real* const pc_costshifts = distdata->dijkdata->node_bias;
272  const int nnodes = g->knots;
273  const SCIP_Bool isPc = graph_pc_isPc(g);
274  SCIP_Real distance = FARAWAY;
275 
276  assert(!isPc || pc_costshifts);
277 
278  SCIP_CALL_ABORT( graph_heap_create(scip, nnodes, NULL, NULL, &dheap) );
279  SCIP_CALL_ABORT( SCIPallocMemoryArray(scip, &dist, nnodes) );
280 
281  state = dheap->position;
282 
283  for( int k = 0; k < nnodes; k++ )
284  {
285  dist[k] = FARAWAY;
286  assert(state[k] == UNKNOWN);
287  }
288 
289  dist[vertex1] = 0.0;
290  graph_heap_correct(vertex1, 0.0, dheap);
291 
292  assert(dheap->size == 1);
293 
294  /* main loop */
295  while( dheap->size > 0 )
296  {
297  /* get nearest labeled node */
298  const int k = graph_heap_deleteMinReturnNode(dheap);
299  const SCIP_Real k_dist = dist[k];
300 
301  if( isPc && k != vertex1 )
302  {
303  dist[k] += pc_costshifts[k];
304  }
305 
306  if( k == vertex2 )
307  {
308  distance = k_dist;
309  break;
310  }
311 
312  /* correct adjacent nodes */
313  for( int e = g->outbeg[k]; e != EAT_LAST; e = g->oeat[e] )
314  {
315  const int m = g->head[e];
316  assert(g->mark[m]);
317 
318  if( m == vertexBlocked )
319  continue;
320 
321  if( state[m] != CONNECT )
322  {
323  const SCIP_Real distnew = isPc ?
324  k_dist + g->cost[e] - pc_costshifts[m]
325  : k_dist + g->cost[e];
326 
327  if( distnew < dist[m] )
328  {
329  dist[m] = distnew;
330  graph_heap_correct(m, distnew, dheap);
331  }
332  }
333  }
334  }
335 
336  SCIPfreeMemoryArray(scip, &dist);
337  graph_heap_free(scip, TRUE, TRUE, &dheap);
338 
339  return distance;
340 }
341 
342 
343 /** Helper.
344  * Gives maximum cost among all close nodes. */
345 static inline
347  int vertex, /**< vertex for which to get the maximum cost*/
348  const DISTDATA* distdata /**< distance data */
349  )
350 {
351  const SCIP_Real* const distances = distdata->closenodes_distances;
352  const RANGE* const range = distdata->closenodes_range;
353  const int range_start = range[vertex].start;
354  const int range_end = range[vertex].end;
355  SCIP_Real maxcost = -FARAWAY;
356 
357  for( int i = range_start; i < range_end; ++i )
358  {
359  const SCIP_Real dist = distances[i];
360 
361  if( dist > maxcost )
362  maxcost = dist;
363  }
364 
365  assert(GE(maxcost, 0.0) || range_start == range_end);
366 
367  return maxcost;
368 }
369 
370 
371 /** Gets close nodes and corresponding distances.
372  * NOTE: needs to correspond to 'distDataComputeCloseNodes' in 'extreduce_dbg.c' */
373 static
375  SCIP* scip, /**< SCIP */
376  const GRAPH* g, /**< graph data structure */
377  const DISTDATA* distdata, /**< distance data */
378  int startvertex, /**< start vertex */
379  int* closenodes_indices, /**< indices of close nodes */
380  SCIP_Real* closenodes_dists, /**< distances of close nodes */
381  int* nclosenodes /**< number of added close nodes */
382  )
383 {
384  SCIP_Real* dist;
385  DHEAP* dheap = NULL;
386  int* state;
387  DCSR* const dcsr = g->dcsr_storage;
388  const RANGE* const RESTRICT range_csr = dcsr->range;
389  const int* const RESTRICT head_csr = dcsr->head;
390  const SCIP_Real* const RESTRICT cost_csr = dcsr->cost;
391  const SCIP_Real* const pc_costshifts = distdata->dijkdata->node_bias;
392  const int nnodes = g->knots;
393  const SCIP_Real closenodes_maxcost = distCloseNodesGetMaxCost(startvertex, distdata);
394  int clodenode_count;
395  const SCIP_Bool isPc = graph_pc_isPc(g);
396 
397  assert(dcsr && g && distdata);
398  assert(startvertex >= 0 && startvertex < g->knots);
399  assert(!isPc || pc_costshifts);
400 
401  SCIP_CALL( graph_heap_create(scip, nnodes, NULL, NULL, &dheap) );
402  SCIP_CALL( SCIPallocMemoryArray(scip, &dist, nnodes) );
403 
404  state = dheap->position;
405 
406  for( int k = 0; k < nnodes; k++ )
407  {
408  dist[k] = FARAWAY;
409  assert(state[k] == UNKNOWN);
410  }
411 
412  clodenode_count = 0;
413  dist[startvertex] = 0.0;
414  graph_heap_correct(startvertex, 0.0, dheap);
415 
416  assert(dheap->size == 1);
417 
418  /* main loop */
419  while( dheap->size > 0 )
420  {
421  /* get nearest labeled node */
422  const int k = graph_heap_deleteMinReturnNode(dheap);
423  const int k_start = range_csr[k].start;
424  const int k_end = range_csr[k].end;
425  const SCIP_Real k_dist = dist[k];
426 
427  if( k != startvertex )
428  {
429  if( isPc )
430  dist[k] += pc_costshifts[k];
431 
432  if( GT(k_dist, closenodes_maxcost) )
433  break;
434 
435  closenodes_indices[clodenode_count] = k;
436  closenodes_dists[clodenode_count] = dist[k];
437 
438  clodenode_count++;
439  }
440 
441  /* correct adjacent nodes */
442  for( int e = k_start; e < k_end; e++ )
443  {
444  const int m = head_csr[e];
445  assert(g->mark[m]);
446 
447  if( state[m] != CONNECT )
448  {
449  const SCIP_Real distnew = isPc ?
450  k_dist + cost_csr[e] - pc_costshifts[m]
451  : k_dist + cost_csr[e];
452 
453  if( distnew < dist[m] )
454  {
455  dist[m] = distnew;
456  graph_heap_correct(m, distnew, dheap);
457  }
458  }
459  }
460  }
461 
462  SCIPfreeMemoryArray(scip, &dist);
463  graph_heap_free(scip, TRUE, TRUE, &dheap);
464 
465  *nclosenodes = clodenode_count;
466 
467  return SCIP_OKAY;
468 }
469 
470 
471 
472 /* prints info about lost node */
473 static
475  const GRAPH* g, /**< graph data structure */
476  const DISTDATA* distdata, /**< distance data */
477  int vertex_base, /**< vertex for check */
478  int index_lost, /**< vertex that is lost */
479  const int* closenodes_indices, /**< indices of newly found close nodes */
480  const SCIP_Real* closenodes_dists, /**< distances of newly found close nodes */
481  int nclosenodes /**< number of newly found close nodes */
482 )
483 {
484 #ifdef SCIP_DEBUG
485  RANGE* const org_range = distdata->closenodes_range;
486  int* const org_indices = distdata->closenodes_indices;
487  SCIP_Real* const org_dists = distdata->closenodes_distances;
488  const int org_start = org_range[vertex_base].start;
489  const int org_end = org_range[vertex_base].end;
490  const int vertex_lost = org_indices[index_lost];
491 
492  assert(org_start <= index_lost && index_lost < org_end);
493 
494  SCIPdebugMessage("could not find vertex_base %d in new close nodes \n", vertex_lost);
495  printf("new nodes: \n");
496 
497  for( int k = 0; k < nclosenodes; ++k )
498  {
499  printf("(idx=%d dist=%f) ", k, closenodes_dists[k]);
500  graph_knot_printInfo(g, closenodes_indices[k]);
501  }
502 
503  printf("original nodes: \n");
504 
505  for( int k = org_start; k < org_end; ++k )
506  {
507  printf("(dist=%f) ", org_dists[k]);
508  graph_knot_printInfo(g, org_indices[k]);
509  }
510 
511 #ifndef NDEBUG
512  {
513  int index = index_lost;
514  int vertex = vertex_lost;
515  const int* const prededges = distdata->closenodes_prededges;
516 
517  assert(prededges);
518 
519  printf("path from lost node=%d to base=%d \n", vertex_lost, vertex_base);
520 
521  while( vertex != vertex_base )
522  {
523  const int prededge = prededges[index];
524 
525  assert(prededge >= 0);
526 
527  graph_edge_printInfo(g, prededge);
528 
529  if( distdata->dijkdata->node_bias )
530  printf("vertex %d shift=%f \n", vertex, distdata->dijkdata->node_bias[vertex]);
531 
532  assert(vertex == g->head[prededge] || vertex == g->tail[prededge]);
533 
534  if( vertex == g->head[prededge] )
535  vertex = g->tail[prededge];
536  else
537  vertex = g->head[prededge];
538 
539  if( vertex != vertex_base )
540  {
541  /* get the index of the next vertex */
542  for( index = org_start; index != org_end; index++ )
543  {
544  if( org_indices[index] == vertex )
545  break;
546  }
547 
548  assert(index != org_end);
549  }
550  }
551  }
552 #endif
553 #endif
554 }
555 
556 /* Is each original close node to 'vertex' included in the 'closenodes_indices' array?
557  * And if so, with correct costs? */
558 static
560  SCIP* scip, /**< SCIP */
561  const GRAPH* g, /**< graph data structure */
562  const DISTDATA* distdata, /**< distance data */
563  int vertex, /**< vertex for check */
564  const int* closenodes_indices, /**< indices of newly found close nodes */
565  const SCIP_Real* closenodes_dists, /**< distances of newly found close nodes */
566  int nclosenodes /**< number of newly found close nodes */
567 )
568 {
569  RANGE* const org_range = distdata->closenodes_range;
570  int* const org_indices = distdata->closenodes_indices;
571  SCIP_Real* const org_dists = distdata->closenodes_distances;
572  int* newnodeindex;
573  const int nnodes = graph_get_nNodes(g);
574  const int org_start = org_range[vertex].start;
575  const int org_end = org_range[vertex].end;
576  const int org_nclosenodes = org_end - org_start;
577  SCIP_Bool isIncluded = TRUE;
578 
579  if( nclosenodes < org_nclosenodes )
580  {
581  SCIPdebugMessage("too few new closenodes! %d < %d \n", nclosenodes, org_nclosenodes);
582  return FALSE;
583  }
584 
585  SCIP_CALL_ABORT( SCIPallocCleanBufferArray(scip, &newnodeindex, nnodes) );
586 
587  /* mark the indices of new close nodes */
588  for( int j = 0; j < nclosenodes; j++ )
589  {
590  const int cnode = closenodes_indices[j];
591 
592  assert(newnodeindex[cnode] == 0);
593 
594  newnodeindex[cnode] = j + 1;
595  }
596 
597  /* the actual checks */
598  for( int i = org_start; i < org_end; ++i )
599  {
600  const SCIP_Real org_dist = org_dists[i];
601  const int org_node = org_indices[i];
602  const int new_index = newnodeindex[org_node] - 1;
603 #ifndef NDEBUG
604  const int prededge = distdata->closenodes_prededges[i];
605 
606  assert(new_index >= -1);
607  assert(distdata->closenodes_prededges);
608  assert(org_node == g->tail[prededge] || org_node == g->head[prededge]);
609 #endif
610 
611  if( new_index == -1 )
612  {
613  distCloseNodesPrintLostNodeInfo(g, distdata, vertex, i, closenodes_indices, closenodes_dists, nclosenodes);
614 
615  isIncluded = FALSE;
616  break;
617  }
618 
619  assert(org_node == closenodes_indices[new_index]);
620 
621  if( !EQ(closenodes_dists[new_index], org_dist) )
622  {
623  SCIPdebugMessage("wrong distances: %f != %f \n", org_dists[i], closenodes_dists[new_index]);
624 
625  isIncluded = FALSE;
626  break;
627  }
628  }
629 
630  /* unmark the new close nodes */
631  for( int j = 0; j < nclosenodes; j++ )
632  {
633  const int cnode = closenodes_indices[j];
634 
635  assert(newnodeindex[cnode] == j + 1);
636 
637  newnodeindex[cnode] = 0;
638  }
639 
640  SCIPfreeCleanBufferArray(scip, &newnodeindex);
641 
642  return isIncluded;
643 }
644 
645 
646 /** computes counters for degrees and edges */
647 static
649  const GRAPH* graph, /**< graph data structure */
650  const EXTDATA* extdata, /**< extension data */
651  int edgecount[], /**< edge count */
652  int degreecount[], /**< degree count */
653  SCIP_Bool* treeIsFlawed /**< is tree flawed?) */
654 )
655 {
656  const int* const tree_edges = extdata->tree_edges;
657  const int* const tree_deg = extdata->tree_deg;
658  const int tree_nedges = extdata->tree_nedges;
659 
660  for( int i = 0; i < tree_nedges; i++ )
661  {
662  const int e = tree_edges[i];
663  const int head = graph->head[e];
664  const int tail = graph->tail[e];
665 
666  assert(e >= 0 && e < graph->edges);
667 
668  if( edgecount[e] > 0 || edgecount[flipedge(e)] > 0 )
669  {
670  printf("tree_nedges %d \n", tree_nedges);
671  printf("FLAW: double edge \n");
672  graph_edge_printInfo(graph, e);
673  *treeIsFlawed = TRUE;
674  }
675 
676  if( tree_deg[tail] <= 0 )
677  {
678  printf("FLAW: non-positive degree for %d (%d) \n", tail, tree_deg[tail]);
679  *treeIsFlawed = TRUE;
680  }
681 
682  if( tree_deg[head] <= 0 )
683  {
684  printf("FLAW: non-positive degree for %d (%d) \n", head, tree_deg[head]);
685  *treeIsFlawed = TRUE;
686  }
687 
688  degreecount[tail]++;
689  degreecount[head]++;
690 
691  edgecount[e]++;
692  }
693 }
694 
695 
696 /** reset counters */
697 static
699  const GRAPH* graph, /**< graph data structure */
700  const EXTDATA* extdata, /**< extension data */
701  int edgecount[], /**< edge count */
702  int degreecount[] /**< degree count */
703 )
704 {
705  const int* const tree_edges = extdata->tree_edges;
706  const int nnodes = graph_get_nNodes(graph);
707  const int nedges = graph_get_nEdges(graph);
708  const int tree_nedges = extdata->tree_nedges;
709 
710  for( int i = 0; i < tree_nedges; i++ )
711  {
712  const int e = tree_edges[i];
713  const int head = graph->head[e];
714  const int tail = graph->tail[e];
715 
716  edgecount[e] = 0;
717  degreecount[tail] = 0;
718  degreecount[head] = 0;
719  }
720 
721  for( int i = 0; i < nnodes; i++ )
722  {
723  assert(degreecount[i] == 0);
724  }
725 
726  for( int i = 0; i < nedges; i++ )
727  {
728  assert(edgecount[i] == 0);
729  }
730 }
731 
732 
733 /** are the degrees flawed? */
734 static
736  const GRAPH* graph, /**< graph data structure */
737  const EXTDATA* extdata, /**< extension data */
738  const int degreecount[] /**< degree count */
739 )
740 {
741  const int* const tree_edges = extdata->tree_edges;
742  const int* const tree_deg = extdata->tree_deg;
743  const int tree_nedges = extdata->tree_nedges;
744 
745  /* degree check */
746  for( int i = 0; i < tree_nedges; i++ )
747  {
748  const int e = tree_edges[i];
749  const int head = graph->head[e];
750 
751  if( degreecount[head] != tree_deg[head] )
752  {
753  printf("FLAW: wrong degree \n");
754  return TRUE;
755  }
756  }
757 
758  return FALSE;
759 }
760 
761 
762 /** flaw with the tree leaves? */
763 static
765  const GRAPH* graph, /**< graph data structure */
766  const EXTDATA* extdata, /**< extension data */
767  const int degreecount[] /**< degree count */
768 )
769 {
770  const int* const tree_edges = extdata->tree_edges;
771  const int* const tree_leaves = extdata->tree_leaves;
772  const int tree_nedges = extdata->tree_nedges;
773  const int nleaves = extdata->tree_nleaves;
774  int leavescount = 1; /* for tail of initial edge */
775 
776  assert(nleaves >= 1);
777 
778  for( int i = 0; i < tree_nedges; i++ )
779  {
780  const int e = tree_edges[i];
781  const int head = graph->head[e];
782 
783  if( degreecount[head] == 1 )
784  leavescount++;
785  }
786 
787  if( leavescount != nleaves )
788  {
789  printf("FLAW wrong leaves count %d != %d \n", leavescount, nleaves);
790  return TRUE;
791  }
792 
793  for( int i = 0; i < nleaves; i++ )
794  {
795  const int leaf = tree_leaves[i];
796 
797  if( degreecount[leaf] != 1 )
798  {
799  printf("FLAW wrong leaf %d degree %d != %d \n", leaf, degreecount[leaf], 1);
800  return TRUE;
801  }
802  }
803 
804  return FALSE;
805 }
806 
807 
808 /** flaw with the inner nodes? */
809 static
811  const GRAPH* graph, /**< graph data structure */
812  const EXTDATA* extdata, /**< extension data */
813  const int degreecount[] /**< degree count */
814 )
815 {
816  const int* const tree_edges = extdata->tree_edges;
817  const int* const tree_innerNodes = extdata->tree_innerNodes;
818  const int tree_nedges = extdata->tree_nedges;
819  const int ninnerNodes = extdata->tree_ninnerNodes;
820  int innercount = 0;
821 
822  for( int i = 0; i < tree_nedges; i++ )
823  {
824  const int e = tree_edges[i];
825  const int head = graph->head[e];
826 
827  if( degreecount[head] > 1 )
828  innercount++;
829  }
830 
831  if( innercount != ninnerNodes )
832  {
833  printf("FLAW wrong inner count %d != %d \n", innercount, ninnerNodes);
834  return TRUE;
835  }
836 
837  for( int i = 0; i < ninnerNodes; i++ )
838  {
839  const int node = tree_innerNodes[i];
840 
841  if( degreecount[node] <= 1 )
842  {
843  printf("FLAW wrong inner node %d degree %d <= %d \n", node, degreecount[node], 1);
844  return TRUE;
845  }
846  }
847 
848  return FALSE;
849 }
850 
851 
852 /** are distance values of tree flawed? */
853 static
855  const GRAPH* graph, /**< graph data structure */
856  const EXTDATA* extdata /**< extension data */
857 )
858 {
859  const int nnodes = graph_get_nNodes(graph);
860  const SCIP_Bool isPc = graph_pc_isPcMw(graph);
861 
862  for( int i = 0; i < nnodes; i++ )
863  {
864  if( isPc && extdata->pcdata->pcSdToNode[i] >= -0.5 )
865  {
866  printf("FLAW wrong pcSdToNode entry[%d]=%f \n", i, extdata->pcdata->pcSdToNode[i]);
867  return TRUE;
868  }
869  if( extdata->tree_bottleneckDistNode[i] >= -0.5 )
870  {
871  printf("FLAW wrong tree_bottleneckDistNode entry[%d]=%f \n", i, extdata->tree_bottleneckDistNode[i]);
872  return TRUE;
873  }
874  }
875 
876  return FALSE;
877 }
878 
879 
880 /** cleans extension data extensively for debugging */
882  const GRAPH* graph, /**< graph data structure */
883  EXTDATA* extdata /**< extension data */
884 )
885 {
886  const int nnodes = graph_get_nNodes(graph);
887  const int maxstacksize = extreduce_getMaxStackSize();
888  const int maxncomponents = extreduce_getMaxStackNcomponents(graph);
889  int* const extstack_data = extdata->extstack_data;
890  int* const extstack_start = extdata->extstack_start;
891  int* const extstack_state = extdata->extstack_state;
892  int* const tree_edges = extdata->tree_edges;
893  int* const tree_leaves = extdata->tree_leaves;
894  int* const tree_innerNodes = extdata->tree_innerNodes;
895  int* const tree_parentNode = extdata->tree_parentNode;
896  SCIP_Real* const tree_parentEdgeCost = extdata->tree_parentEdgeCost;
897  SCIP_Real* const redcost_treenodeswaps = extdata->reddata->redcost_treenodeswaps;
898  const int redcost_nlevels = extdata->reddata->redcost_nlevels;
899 
900  for( int i = 0; i < maxstacksize; ++i )
901  {
902  extstack_data[i] = INT_MIN;
903  }
904 
905  for( int i = 0; i < maxncomponents + 1; ++i )
906  {
907  extstack_start[i] = INT_MIN;
908  extstack_state[i] = INT_MIN;
909  }
910 
911  for( int i = 0; i < nnodes; ++i )
912  {
913  tree_edges[i] = INT_MIN;
914  tree_leaves[i] = INT_MIN;
915  tree_innerNodes[i] = INT_MIN;
916  tree_parentNode[i] = INT_MIN;
917  }
918 
919  for( int i = 0; i < nnodes; ++i )
920  {
921  tree_edges[i] = INT_MIN;
922  tree_leaves[i] = INT_MIN;
923  }
924 
925 
926  for( int i = 0; i < nnodes; ++i )
927  {
928  tree_parentEdgeCost[i] = -FARAWAY;
929  }
930 
931  for( int i = 0; i < nnodes * redcost_nlevels; i++ )
932  {
933  redcost_treenodeswaps[i] = -FARAWAY;
934  }
935 }
936 
937 
938 /** is current tree flawed? */
940  SCIP* scip, /**< SCIP */
941  const GRAPH* graph, /**< graph data structure */
942  const EXTDATA* extdata /**< extension data */
943 )
944 {
945  int* edgecount;
946  int* degreecount;
947  const int nnodes = graph_get_nNodes(graph);
948  const int nedges = graph_get_nEdges(graph);
949  SCIP_Bool flawed = FALSE;
950 
951  SCIP_CALL_ABORT( SCIPallocCleanBufferArray(scip, &edgecount, nedges) );
952  SCIP_CALL_ABORT( SCIPallocCleanBufferArray(scip, &degreecount, nnodes) );
953 
954  treeGetCounters(graph, extdata, edgecount, degreecount, &flawed);
955 
956  if( flawed )
957  {
958  flawed = TRUE;
959  }
960  else if( treeDegreesAreFlawed(graph, extdata, degreecount) )
961  {
962  flawed = TRUE;
963  }
964  else if( treeLeavesAreFlawed(graph, extdata, degreecount) )
965  {
966  flawed = TRUE;
967  }
968  else if( treeInnerNodesAreFlawed(graph, extdata, degreecount) )
969  {
970  flawed = TRUE;
971  }
972  else if( treeDistsAreFlawed(graph, extdata) )
973  {
974  flawed = TRUE;
975  }
976 
977  treeResetCounters(graph, extdata, edgecount, degreecount);
978 
979  SCIPfreeCleanBufferArray(scip, &degreecount);
980  SCIPfreeCleanBufferArray(scip, &edgecount);
981 
982  return flawed;
983 }
984 
985 
986 /** is current tree completely hashed? */
988  const GRAPH* graph, /**< graph data structure */
989  const EXTDATA* extdata /**< extension data */
990 )
991 {
992  const REDDATA* const reddata = extdata->reddata;
993 
994  for( int i = 0; i < extdata->tree_nedges; i++ )
995  {
996  const int edge = extdata->tree_edges[i];
997  const int nAncestors = graph_edge_nPseudoAncestors(graph, edge);
998 
999  assert(nAncestors >= 0);
1000 
1001  if( nAncestors == 0 )
1002  continue;
1003 
1005  return FALSE;
1006  }
1007 
1008  return TRUE;
1009 }
1010 
1011 
1012 /** prints the leaves of the tree */
1014  const EXTDATA* extdata /**< extension data */
1015 )
1016 {
1017  const int nleaves = extdata->tree_nleaves;
1018  const int* const leaves = extdata->tree_leaves;
1019 
1020  printf("tree leaves: \n");
1021 
1022  for( int i = 0; i < nleaves; ++i )
1023  {
1024  const int leaf = leaves[i];
1025  printf("%d: leaf=%d \n", i, leaf);
1026  }
1027 }
1028 
1029 
1030 /** prints the current stack */
1032  const GRAPH* graph, /**< graph data structure */
1033  const EXTDATA* extdata /**< extension data */
1034 )
1035 {
1036  const int* const extstack_data = extdata->extstack_data;
1037  const int* const extstack_start = extdata->extstack_start;
1038  const int stackpos = extdata->extstack_ncomponents - 1;
1039 
1040  for( int j = 0; j <= stackpos; j++ )
1041  {
1042  if( extdata->extstack_state[j] == EXT_STATE_NONE )
1043  printf("pos=%d state=NONE \n", j);
1044  else if( extdata->extstack_state[j] == EXT_STATE_EXPANDED )
1045  printf("pos=%d state=EXPANDED \n", j);
1046  else
1047  {
1048  assert(extdata->extstack_state[j] == EXT_STATE_MARKED);
1049 
1050  printf("pos=%d state=MARKED \n", j);
1051  }
1052 
1053  /* check all leaves of current component */
1054  for( int i = extstack_start[j]; i < extstack_start[j + 1]; i++ )
1055  {
1056  const int edge = extstack_data[i];
1057 
1058  if( edge == EXT_EDGE_WRAPPED )
1059  {
1060  printf(" EDGE_WRAPPED\n");
1061  continue;
1062  }
1063 
1064  assert(edge >= 0 && edge < graph->edges);
1065  printf(" ");
1066  graph_edge_printInfo(graph, edge);
1067  }
1068  }
1069 }
1070 
1071 
1072 /** Prints top horizontal level */
1074  const EXTDATA* extdata /**< extension data */
1075 )
1076 {
1077  const REDDATA* const reddata = extdata->reddata;
1078  const MLDISTS* const sds_horizontal = reddata->sds_horizontal;
1079  const int levelsize = extreduce_mldistsTopLevelNSlots(sds_horizontal);
1080  const int* const levelids = extreduce_mldistsTopLevelBases(sds_horizontal);
1081 
1082  printf("current (horizontal) top level: \n");
1083 
1084  for( int i = 0; i < levelsize; i++ )
1085  {
1086  printf("%d ", levelids[i]);
1087  }
1088 
1089  printf("\n");
1090 }
1091 
1092 
1093 /** is the node in the current top component of the stack? */
1095  const GRAPH* graph, /**< graph data structure */
1096  const EXTDATA* extdata, /**< extension data */
1097  int node /**< the node */
1098  )
1099 {
1100  const int stackpos = extStackGetPosition(extdata);
1101  const int* const stack_start = extdata->extstack_start;
1102  const int* const stack_data = extdata->extstack_data;
1103 
1104  assert(graph);
1105  assert(stack_start && stack_data);
1106  assert(node >= 0 && node < graph->knots);
1107  assert(extdata->extstack_state[stackpos] == EXT_STATE_EXPANDED || extdata->extstack_state[stackpos] == EXT_STATE_MARKED);
1108 
1109  for( int i = stack_start[stackpos]; i < stack_start[stackpos + 1]; i++ )
1110  {
1111  const int topnode = graph->head[stack_data[i]];
1112 
1113  assert(topnode >= 0 && topnode < graph->knots);
1114 
1115  if( topnode == node )
1116  return TRUE;
1117  }
1118 
1119  return FALSE;
1120 }
1121 
1122 
1123 /** Are the close-nodes still valid?
1124  * NOTE: expensive method, just designed for debugging! */
1126  SCIP* scip, /**< SCIP */
1127  const GRAPH* g, /**< graph data structure */
1128  const DISTDATA* distdata /**< distance data */
1129 )
1130 {
1131  SCIP_Real* closenodes_dists;
1132  int* closenodes_indices;
1133  int nclosenodes = 0;
1134  const int nnodes = graph_get_nNodes(g);
1135  SCIP_Bool isValid = TRUE;
1136 
1137  SCIP_CALL_ABORT( SCIPallocMemoryArray(scip, &closenodes_indices, nnodes) );
1138  SCIP_CALL_ABORT( SCIPallocMemoryArray(scip, &closenodes_dists, nnodes) );
1139 
1140  assert(scip && distdata);
1141  assert(distdata->pathroot_isdirty);
1142 
1143  for( int i = 0; i < nnodes; i++ )
1144  {
1145  if( distdata->pathroot_isdirty[i] )
1146  continue;
1147 
1148  SCIP_CALL_ABORT( distCloseNodesCompute(scip, g, distdata, i,
1149  closenodes_indices, closenodes_dists, &nclosenodes) );
1150 
1151  if( !distCloseNodesIncluded(scip, g, distdata, i, closenodes_indices, closenodes_dists, nclosenodes) )
1152  {
1153 #ifdef SCIP_DEBUG
1154  SCIPdebugMessage("corrupted node: ");
1155  graph_knot_printInfo(g, i);
1156 #endif
1157 
1158  isValid = FALSE;
1159  break;
1160  }
1161  }
1162 
1163  SCIPfreeMemoryArray(scip, &closenodes_dists);
1164  SCIPfreeMemoryArray(scip, &closenodes_indices);
1165 
1166  return isValid;
1167 }
1168 
1169 
1170 /** Computes actual distance between two nodes.
1171  * NOTE: expensive method, just designed for debugging! */
1173  SCIP* scip, /**< SCIP */
1174  const GRAPH* g, /**< graph data structure */
1175  int vertexBlocked, /**< forbidden vertex */
1176  const DISTDATA* distdata, /**< distance data */
1177  int vertex1, /**< first vertex */
1178  int vertex2 /**< second vertex */
1179 )
1180 {
1181  SCIP_Real dist;
1182 
1183  assert(scip && g && distdata);
1184  assert(graph_knot_isInRange(g, vertex1));
1185  assert(graph_knot_isInRange(g, vertex2));
1186  assert(vertex1 != vertex2);
1187  assert(vertexBlocked != vertex1 && vertexBlocked != vertex2);
1188 
1189  dist = distGetRestricted(scip, g, vertexBlocked, distdata, vertex1, vertex2);
1190 
1191  return dist;
1192 }
1193 
1194 
1195 /** debug initialization */
1197  int* extedgesstart, /**< array */
1198  int* extedges /**< array */
1199 )
1200 {
1201  assert(extedgesstart && extedges);
1202 
1203  for( int i = 0; i < STP_EXT_MAXGRAD; i++ )
1204  extedgesstart[i] = -1;
1205 
1206  for( int i = 0; i < STP_EXT_MAXGRAD * STP_EXT_MAXGRAD; i++ )
1207  extedges[i] = -1;
1208 }
1209 
1210 
1211 /** check whether vertical SDs are up to date for given leaf of component */
1213  SCIP* scip, /**< SCIP */
1214  const GRAPH* graph, /**< graph data structure */
1215  int compsize, /**< size of component */
1216  int nleaves_ancestors, /**< number of leaves to ancestors */
1217  int topleaf, /**< component leaf to check for */
1218  EXTDATA* extdata /**< extension data */
1219  )
1220 {
1221  const MLDISTS* const sds_vertical = extdata->reddata->sds_vertical;
1222  const int* const leaves = extdata->tree_leaves;
1223  const int nleaves = extdata->tree_nleaves;
1224  const SCIP_Real* const adjedgecosts = extreduce_mldistsTopTargetDists(sds_vertical, topleaf);
1225  const SCIP_Bool atInitialComp = extIsAtInitialComp(extdata);
1226  const int nleaves_old = atInitialComp ? 1 : (nleaves - compsize);
1227  SCIP_Bool isInSync = TRUE;
1228 #ifndef NDEBUG
1229  const int* const adjids = extreduce_mldistsTopTargetIds(sds_vertical, topleaf);
1230 #endif
1231 
1232  assert(adjedgecosts);
1233  assert(nleaves_old == nleaves_ancestors);
1234  assert(nleaves_old > 0 && nleaves_old < nleaves);
1235 
1236 #ifndef STP_DEBUG_EXTPC
1237  if( graph_pc_isPc(graph) )
1238  return TRUE;
1239 #endif
1240 
1241  /* get the SDs to the ancestor (lower) leafs and compare */
1242  for( int j = 0; j < nleaves_old; j++ )
1243  {
1244  const int leaf = leaves[j];
1245  const SCIP_Real sd_old = adjedgecosts[j];
1246  const SCIP_Real specialDist_new = extreduce_extGetSd(scip, graph, topleaf, leaf, extdata);
1247  const SCIP_Real sd_new = (specialDist_new >= -0.5) ? specialDist_new : FARAWAY;
1248 
1249  assert(!extreduce_nodeIsInStackTop(graph, extdata, leaf));
1250  assert(extdata->tree_deg[leaf] == 1);
1251  assert(leaf != topleaf);
1252  assert(adjids[j] == leaf);
1253 
1254  if( !EQ(sd_old, sd_new) )
1255  {
1256  SCIPdebugMessage("vertical SDs are wrong! %f!=%f \n", sd_old, sd_new);
1257 
1258  isInSync = FALSE;
1259  break;
1260  }
1261  }
1262 
1263 
1264  return isInSync;
1265 }
1266 
1267 
1268 /** check whether horizontal SDs are up to date for given leaf of component */
1270  SCIP* scip, /**< SCIP */
1271  const GRAPH* graph, /**< graph data structure */
1272  int topleaf, /**< component leaf to check for */
1273  EXTDATA* extdata /**< extension data */
1274  )
1275 {
1276  const MLDISTS* const sds_horizontal = extdata->reddata->sds_horizontal;
1277  const int* const extstack_data = extdata->extstack_data;
1278  const int* const ghead = graph->head;
1279  const int stackpos = extStackGetPosition(extdata);
1280  const int stackstart = extStackGetTopOutEdgesStart(extdata, stackpos);
1281  const int stackend = extStackGetTopOutEdgesEnd(extdata, stackpos);
1282  SCIP_Bool isInSync = TRUE;
1283 #ifndef NDEBUG
1284  SCIP_Bool hitTopLeaf = FALSE;
1285 #endif
1286 
1287 #ifndef STP_DEBUG_EXTPC
1288  if( graph_pc_isPc(graph) )
1289  return TRUE;
1290 #endif
1291 
1292  for( int i = stackstart, j = 0; i < stackend; i++, j++ )
1293  {
1294  const int edge2sibling = extstack_data[i];
1295  const int sibling = ghead[edge2sibling];
1296 
1297  assert(extreduce_nodeIsInStackTop(graph, extdata, sibling));
1298  assert(extdata->tree_deg[sibling] == 1);
1299 
1300  if( sibling == topleaf )
1301  {
1302 #ifndef NDEBUG
1303  hitTopLeaf = TRUE;
1304 #endif
1305  continue;
1306  }
1307  else
1308  {
1309  const SCIP_Real sd_old = extreduce_mldistsTopTargetDist(sds_horizontal, topleaf, sibling);
1310  const SCIP_Real sd_new = extreduce_extGetSdProperDouble(scip, graph, topleaf, sibling, extdata);
1311 
1312  if( !EQ(sd_old, sd_new) )
1313  {
1314  SCIPdebugMessage("vertical SDs are wrong! %f!=%f \n", sd_old, sd_new);
1315 
1316  isInSync = FALSE;
1317  break;
1318  }
1319  }
1320  }
1321 
1322  assert(hitTopLeaf || !isInSync);
1323 
1324  return isInSync;
1325 }
1326 
1327 
1328 /** are sds from top component leaf corresponding to current tree? */
1330  SCIP* scip, /**< SCIP */
1331  const GRAPH* graph, /**< graph data structure */
1332  const SCIP_Real sds[], /**< SDs from top leaf */
1333  int topleaf, /**< component leaf to check for */
1334  EXTDATA* extdata /**< extension data */
1335  )
1336 {
1337  const int* const leaves = extdata->tree_leaves;
1338  const int nleaves = extdata->tree_nleaves;
1339  SCIP_Bool isInSync = TRUE;
1340 
1341 #ifndef STP_DEBUG_EXTPC
1342  if( graph_pc_isPc(graph) )
1343  return TRUE;
1344 #endif
1345 
1346  for( int j = 0; j < nleaves; j++ )
1347  {
1348  const int leaf = leaves[j];
1349 
1350  if( leaf != topleaf )
1351  {
1352  const SCIP_Real sd_double = extreduce_extGetSdProperDouble(scip, graph, topleaf, leaf, extdata);
1353  const SCIP_Real sd = extreduce_extGetSdProper(scip, graph, topleaf, leaf, extdata);
1354 
1355  if( !EQ(sds[j], sd_double) && !EQ(sds[j], sd) )
1356  {
1357  SCIPdebugMessage("SD from %d to %d not correct! \n", topleaf, leaf);
1358  SCIPdebugMessage("new sds (double, single): %f, %f ... old sd: %f \n", sd_double, sd, sds[j]);
1359 
1360  isInSync = FALSE;
1361  break;
1362  }
1363  }
1364  else
1365  {
1366  if( !EQ(sds[j], FARAWAY) )
1367  {
1368  SCIPdebugMessage("SD to topleaf not FARAWAY! (but %f) \n", sds[j]);
1369 
1370  isInSync = FALSE;
1371  break;
1372  }
1373  }
1374  }
1375 
1376  return isInSync;
1377 }
1378 
1379 
1380 /** is the top level base MST objective in sync with the current tree? */
1382  SCIP* scip, /**< SCIP */
1383  const GRAPH* graph, /**< graph data structure */
1384  int extnode, /**< node from which the level was extended */
1385  EXTDATA* extdata /**< extension data */
1386 )
1387 {
1388  SCIP_Bool isValid = TRUE;
1389 
1390  assert(extnode >= 0 && extnode < graph->knots);
1391  assert(extdata->tree_deg[extnode] == 1);
1392 
1393  if( !graph_pc_isPcMw(graph) )
1394  {
1395  CSR mst_base;
1396  const REDDATA* const reddata = extdata->reddata;
1397  const CSRDEPO* const msts_levelbase = reddata->msts_levelbase;
1398  int* nodes;
1399  const int nleaves = extdata->tree_nleaves;
1400  const int nnodes = nleaves - 1;
1401 
1402  assert(nnodes >= 1);
1403 
1404  SCIP_CALL_ABORT( SCIPallocMemoryArray(scip, &nodes, nnodes) );
1405 
1406  graph_csrdepo_getTopCSR(msts_levelbase, &mst_base);
1407 
1408  assert(nnodes == mst_base.nnodes);
1409  assert(reduce_dcmstMstIsValid(scip, &mst_base));
1410 
1411  /* get the nodes for the new MST computations */
1412  mstTopLevelBaseGetNodes(extnode, nnodes, nodes, extdata);
1413 
1414  if( !mstTopLevelBaseValidWeight(scip, graph, nnodes, nodes, &mst_base, extdata) )
1415  {
1416  isValid = FALSE;
1417  }
1418 
1419  SCIPfreeMemoryArray(scip, &nodes);
1420  }
1421 
1422  return isValid;
1423 }
1424 
1425 
1426 /** is the objective of the top MST extension valid for the tree? */
1428  SCIP* scip, /**< SCIP */
1429  const GRAPH* graph, /**< graph data structure */
1430  int extvert, /**< extended vertex */
1431  SCIP_Real extobj, /**< objective of extension */
1432  EXTDATA* extdata /**< extension data */
1433 )
1434 {
1435  SCIP_Bool isValid = TRUE;
1436 
1437  /** does currently not work because of weird SD computation */
1438  if( !graph_pc_isPcMw(graph) )
1439  {
1440  const SCIP_Real mstweight_upper = sdmstGetExtWeight(scip, graph, extvert, FALSE, extdata);
1441 
1442  if( GT(extobj, mstweight_upper) )
1443  {
1444  SCIPdebugMessage("top extension MST weight violates upper bound: %f > %f \n",
1445  extobj, mstweight_upper);
1446 
1447  isValid = FALSE;
1448  }
1449  else
1450  {
1451  const SCIP_Real mstweight_lower = sdmstGetExtWeight(scip, graph, extvert, TRUE, extdata);
1452 
1453  if( LT(extobj, mstweight_lower) )
1454  {
1455  SCIPdebugMessage("top extension MST weight violates lower bound: %f < %f \n",
1456  extobj, mstweight_lower);
1457 
1458  isValid = FALSE;
1459  }
1460  }
1461  }
1462 
1463  return isValid;
1464 }
1465 
1466 
1467 /** is the objective of the top MST sync with the tree? */
1469  SCIP* scip, /**< SCIP */
1470  const GRAPH* graph, /**< graph data structure */
1471  SCIP_Real compobj, /**< alleged objective of component */
1472  EXTDATA* extdata /**< extension data */
1473 )
1474 {
1475  SCIP_Bool isValid = TRUE;
1476 
1477  if( !graph_pc_isPcMw(graph) )
1478  {
1479  const int nleaves = extdata->tree_nleaves;
1480  const int* const leaves = extdata->tree_leaves;
1481  const SCIP_Real mstweight_upper = sdmstGetWeight(scip, graph, nleaves, leaves, FALSE, extdata);
1482  const SCIP_Real mstweight_lower = sdmstGetWeight(scip, graph, nleaves, leaves, TRUE, extdata);
1483 
1484  if( GT(compobj, mstweight_upper) )
1485  {
1486  SCIPdebugMessage("top MST weight violates upper bound: %f > %f \n",
1487  compobj, mstweight_upper);
1488 
1489  isValid = FALSE;
1490  }
1491  else if( LT(compobj, mstweight_lower) )
1492  {
1493  SCIPdebugMessage("top MST weight violates lower bound: %f < %f \n",
1494  compobj, mstweight_lower);
1495 
1496  isValid = FALSE;
1497  }
1498  }
1499 
1500  return isValid;
1501 }
1502 
1503 
1504 /** is the top MST sync with the tree? (does not check objective) */
1506  SCIP* scip, /**< SCIP */
1507  const GRAPH* graph, /**< graph data structure */
1508  EXTDATA* extdata /**< extension data */
1509 )
1510 {
1511  CSR topmst;
1512  const REDDATA* const reddata = extdata->reddata;
1513  const CSRDEPO* const msts_comp = reddata->msts_comp;
1514  const int* const leaves = extdata->tree_leaves;
1515  const int nleaves = extdata->tree_nleaves;
1516 
1517  graph_csrdepo_getTopCSR(msts_comp, &topmst);
1518 
1519  if( topmst.nnodes != nleaves )
1520  {
1521  SCIPdebugMessage("wrong top mst size nodes=%d leaves=%d \n", topmst.nnodes, nleaves);
1522 
1523  return FALSE;
1524  }
1525 
1526  if( graph_pc_isPcMw(graph) )
1527  {
1528  return TRUE;
1529  }
1530 
1531  /* make sure that the edge costs of the MST correspond to the SDs in the tree */
1532  for( int i = 0; i < nleaves; i++ )
1533  {
1534  const int startleaf = leaves[i];
1535 
1536  for( int j = topmst.start[i]; j != topmst.start[i + 1]; j++ )
1537  {
1538  const int head = topmst.head[j];
1539  const SCIP_Real edgecost = topmst.cost[j];
1540  int endleaf;
1541  SCIP_Real sd_lower;
1542  SCIP_Real sd_upper;
1543 
1544  assert(head >= 0 && head < nleaves);
1545 
1546  endleaf = leaves[head];
1547  sd_lower = extreduce_extGetSdProperDouble(scip, graph, startleaf, endleaf, extdata);
1548  sd_upper = MAX(extreduce_extGetSdProper(scip, graph, startleaf, endleaf, extdata),
1549  extreduce_extGetSdProper(scip, graph, endleaf, startleaf, extdata));
1550 
1551  if( LT(edgecost, sd_lower) )
1552  {
1553  SCIPdebugMessage("MST edge %d->%d violates lower bound: %f < %f \n",
1554  startleaf, endleaf , edgecost, sd_lower);
1555 
1556  return FALSE;
1557  }
1558  else if( GT(edgecost, sd_upper) )
1559  {
1560  SCIPdebugMessage("MST edge %d->%d violates upper bound: %f > %f \n",
1561  startleaf, endleaf , edgecost, sd_upper);
1562 
1563  return FALSE;
1564  }
1565  }
1566  }
1567 
1568  return TRUE;
1569 }
1570 
1571 
1572 /** are the internal data MST structures in sync. with each other? */
1574  const EXTDATA* extdata /**< extension data */
1575 )
1576 {
1577  const REDDATA* const reddata = extdata->reddata;
1578  const CSRDEPO* const msts_comp = reddata->msts_comp;
1579  const CSRDEPO* const msts_levelbase = reddata->msts_levelbase;
1580  const MLDISTS* const sds_vertical = reddata->sds_vertical;
1581  const MLDISTS* const sds_horizontal = reddata->sds_horizontal;
1582 
1583  assert(!extreduce_mldistsEmptySlotExists(sds_horizontal));
1584  assert(!extreduce_mldistsEmptySlotExists(sds_vertical));
1585 
1586  if( extreduce_mldistsNlevels(sds_horizontal) != extreduce_mldistsNlevels(sds_vertical) )
1587  {
1588  SCIPdebugMessage("MLDISTS out of sync %d!=%d \n", extreduce_mldistsNlevels(sds_horizontal),
1589  extreduce_mldistsNlevels(sds_vertical));
1590 
1591  return FALSE;
1592  }
1593 
1594  if( graph_csrdepo_getNcsrs(msts_comp) != graph_csrdepo_getNcsrs(msts_levelbase) )
1595  {
1596  SCIPdebugMessage("CSR depos out of sync %d!=%d \n", graph_csrdepo_getNcsrs(msts_comp),
1597  graph_csrdepo_getNcsrs(msts_levelbase));
1598 
1599  return FALSE;
1600  }
1601 
1602  if (extreduce_mldistsTopLevelNSlots(sds_horizontal) != extreduce_mldistsTopLevelNSlots(sds_vertical) )
1603  {
1604  SCIPdebugMessage("sds out of sync %d!=%d \n", extreduce_mldistsTopLevelNSlots(sds_horizontal),
1605  extreduce_mldistsTopLevelNSlots(sds_vertical));
1606 
1607  return FALSE;
1608  }
1609 
1610  return TRUE;
1611 }
1612 
1613 
1614 /** is that complete current stack hashed? */
1616  const GRAPH* graph, /**< graph data structure */
1617  const EXTDATA* extdata /**< extension data */
1618 )
1619 {
1620  const int* const extstack_data = extdata->extstack_data;
1621  const int* const extstack_start = extdata->extstack_start;
1622  const int* const hasharr = extdata->reddata->pseudoancestor_mark;
1623  const int stackpos = extStackGetPosition(extdata);
1624 
1625  for( int i = extstack_start[stackpos]; i < extstack_start[stackpos + 1]; i++ )
1626  {
1627  const int edge = extstack_data[i];
1628 
1629  if( graph_edge_nPseudoAncestors(graph, edge) == 0 )
1630  continue;
1631 
1632  if( !graph_pseudoAncestors_edgeIsHashed(graph->pseudoancestors, edge, hasharr) )
1633  {
1634  SCIPdebugMessage("stacktop edge %d is not hashed! \n", edge);
1635 
1636  return FALSE;
1637  }
1638  }
1639 
1640  return TRUE;
1641 }
1642 
1643 
1644 /** returns size of component on the stack */
1646  const EXTDATA* extdata, /**< extension data */
1647  int stackpos /**< position on the stack */
1648 )
1649 {
1650  const int* const stack_start = extdata->extstack_start;
1651  int size = stack_start[stackpos + 1] - stack_start[stackpos];
1652  const SCIP_Bool atInitialStar = extInitialCompIsStar(extdata) && (stackpos == 0);
1653  const SCIP_Bool atInitialGenStar = extInitialCompIsGenStar(extdata) && (stackpos == 0);
1654 
1655  if( atInitialStar )
1656  {
1657  size -= 1;
1658  }
1659  else if( atInitialGenStar )
1660  {
1661  size -= 2;
1662  }
1663 
1664  assert(extdata->extstack_state[stackpos] != EXT_STATE_NONE);
1665  assert(size > 0 && size < STP_EXT_MAXGRAD);
1666 
1667  return size;
1668 }
int graph_get_nEdges(const GRAPH *)
Definition: graph_stats.c:548
int graph_edge_nPseudoAncestors(const GRAPH *, int)
SCIP_RETCODE cmst_init(SCIP *scip, CMST **cmst, int maxnnodes)
int extreduce_extStackCompNOutedges(const EXTDATA *extdata, int stackpos)
static SCIP_Bool treeDegreesAreFlawed(const GRAPH *graph, const EXTDATA *extdata, const int degreecount[])
void extreduce_extendInitDebug(int *extedgesstart, int *extedges)
static SCIP_Bool extInitialCompIsGenStar(const EXTDATA *extdata)
SCIP_Bool extreduce_mstTopCompExtObjValid(SCIP *scip, const GRAPH *graph, int extvert, SCIP_Real extobj, EXTDATA *extdata)
SCIP_Bool extreduce_treeIsHashed(const GRAPH *graph, const EXTDATA *extdata)
int *RESTRICT head
Definition: graphdefs.h:224
#define EXT_EDGE_WRAPPED
Definition: extreducedefs.h:49
int * head
Definition: graphdefs.h:141
SCIP_Bool extreduce_distCloseNodesAreValid(SCIP *scip, const GRAPH *g, const DISTDATA *distdata)
int start
Definition: graphdefs.h:152
#define SCIPfreeMemoryArray(scip, ptr)
Definition: scip_mem.h:71
static SCIP_Bool treeInnerNodesAreFlawed(const GRAPH *graph, const EXTDATA *extdata, const int degreecount[])
void extreduce_extdataCleanArraysDbg(const GRAPH *graph, EXTDATA *extdata)
SCIP_Bool extreduce_mldistsEmptySlotExists(const MLDISTS *)
#define SCIPallocMemoryArray(scip, ptr, num)
Definition: scip_mem.h:55
SCIP_Bool extreduce_mstTopLevelBaseObjValid(SCIP *scip, const GRAPH *graph, int extnode, EXTDATA *extdata)
int extreduce_getMaxStackNcomponents(const GRAPH *)
SCIP_Bool extreduce_stackTopIsHashed(const GRAPH *graph, const EXTDATA *extdata)
void graph_knot_printInfo(const GRAPH *, int)
Definition: graph_stats.c:151
void cgraph_node_append(CGRAPH *cgraph, int nodeid)
#define FALSE
Definition: def.h:87
int *const pseudoancestor_mark
static SCIP_Bool mstTopLevelBaseValidWeight(SCIP *scip, const GRAPH *graph, int base_nnodes, int base_nodes[], const CSR *mst_base, EXTDATA *extdata)
SCIP_Real * cost
Definition: graphdefs.h:142
#define TRUE
Definition: def.h:86
int extreduce_mldistsNlevels(const MLDISTS *)
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
includes various files containing graph methods used for Steiner tree problems
int *const extstack_data
int *const tree_edges
SCIP_Real * cost
Definition: graphdefs.h:164
SCIP_Real * node_bias
Definition: graphdefs.h:319
static int extStackGetTopOutEdgesStart(const EXTDATA *extdata, int stackpos)
void graph_csrdepo_getTopCSR(const CSRDEPO *, CSR *)
Definition: graph_util.c:684
#define EXT_STATE_EXPANDED
Definition: extreducedefs.h:53
#define EAT_LAST
Definition: graphdefs.h:58
static SCIP_Bool extInitialCompIsStar(const EXTDATA *extdata)
SCIP_Real extreduce_extGetSdProperDouble(SCIP *, const GRAPH *, int, int, EXTDATA *)
#define SCIPdebugMessage
Definition: pub_message.h:87
#define FARAWAY
Definition: graphdefs.h:89
static int extStackGetTopOutEdgesEnd(const EXTDATA *extdata, int stackpos)
static SCIP_Bool treeLeavesAreFlawed(const GRAPH *graph, const EXTDATA *extdata, const int degreecount[])
int * closenodes_prededges
Definition: extreducedefs.h:86
static int extStackGetPosition(const EXTDATA *extdata)
SCIP_Bool extreduce_mstTopCompInSync(SCIP *scip, const GRAPH *graph, EXTDATA *extdata)
SCIP_Bool extreduce_mstInternalsInSync(const EXTDATA *extdata)
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip_mem.h:127
RANGE * closenodes_range
Definition: extreducedefs.h:84
int * closenodes_indices
Definition: extreducedefs.h:85
int * start
Definition: graphdefs.h:140
REDDATA *const reddata
SCIP_Bool extreduce_sdshorizontalInSync(SCIP *scip, const GRAPH *graph, int topleaf, EXTDATA *extdata)
void extreduce_printTopLevel(const EXTDATA *extdata)
SCIP_Bool extreduce_nodeIsInStackTop(const GRAPH *graph, const EXTDATA *extdata, int node)
int *RESTRICT mark
Definition: graphdefs.h:198
int extreduce_mldistsTopLevelNSlots(const MLDISTS *)
#define SCIPallocCleanBufferArray(scip, ptr, num)
Definition: scip_mem.h:133
int * position
Definition: graphdefs.h:305
int *const tree_leaves
int *RESTRICT oeat
Definition: graphdefs.h:231
This file implements extended reduction techniques for several Steiner problems.
#define GE(a, b)
Definition: portab.h:85
#define EXT_STATE_MARKED
Definition: extreducedefs.h:54
SCIP_Bool extreduce_mstTopCompObjValid(SCIP *scip, const GRAPH *graph, SCIP_Real compobj, EXTDATA *extdata)
SCIP_Real *const pcSdToNode
static void treeGetCounters(const GRAPH *graph, const EXTDATA *extdata, int edgecount[], int degreecount[], SCIP_Bool *treeIsFlawed)
void extreduce_printStack(const GRAPH *graph, const EXTDATA *extdata)
int *const tree_parentNode
static SCIP_Real sdmstGetWeight(SCIP *scip, const GRAPH *graph, int nnodes, const int nodes[], SCIP_Bool computeLower, EXTDATA *extdata)
Definition: extreduce_dbg.c:42
PSEUDOANS * pseudoancestors
Definition: graphdefs.h:236
int *const tree_innerNodes
static SCIP_Bool extIsAtInitialComp(const EXTDATA *extdata)
static SCIP_Real sdmstGetExtWeight(SCIP *scip, const GRAPH *graph, int extvert, SCIP_Bool computeLower, EXTDATA *extdata)
static SCIP_Real distGetRestricted(SCIP *scip, const GRAPH *g, int vertexBlocked, const DISTDATA *distdata, int vertex1, int vertex2)
SCIP_Real reduce_dcmstGetWeight(SCIP *, const CSR *)
Definition: reduce_util.c:1573
void cmst_free(SCIP *scip, CMST **cmst)
SCIP_Real extreduce_extGetSdProper(SCIP *, const GRAPH *, int, int, EXTDATA *)
SCIP_Real * closenodes_distances
Definition: extreducedefs.h:87
#define NULL
Definition: lpi_spx1.cpp:155
SCIP_Real extreduce_extGetSdDouble(SCIP *, const GRAPH *, int, int, EXTDATA *)
int *const extstack_state
SCIP_Real extreduce_extGetSd(SCIP *, const GRAPH *, int, int, EXTDATA *)
#define EQ(a, b)
Definition: portab.h:79
int knots
Definition: graphdefs.h:190
#define SCIP_CALL(x)
Definition: def.h:384
#define LT(a, b)
Definition: portab.h:81
SCIP_Real extreduce_mldistsTopTargetDist(const MLDISTS *, int, int)
void cgraph_free(SCIP *scip, CGRAPH **cgraph)
CSRDEPO *const msts_comp
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:115
#define STP_EXTTREE_MAXNLEAVES_GUARD
Definition: extreducedefs.h:48
void graph_heap_free(SCIP *, SCIP_Bool, SCIP_Bool, DHEAP **)
Definition: graph_util.c:1034
#define GT(a, b)
Definition: portab.h:84
SCIP_Bool reduce_dcmstMstIsValid(SCIP *, const CSR *)
Definition: reduce_util.c:1650
#define SCIP_Bool
Definition: def.h:84
const SCIP_Real * extreduce_mldistsTopTargetDists(const MLDISTS *, int)
SCIP_Bool extreduce_sdsTopInSync(SCIP *scip, const GRAPH *graph, const SCIP_Real sds[], int topleaf, EXTDATA *extdata)
SCIP_RETCODE cgraph_init(SCIP *scip, CGRAPH **cgraph, int maxnnodes)
SCIP_Bool extreduce_treeIsFlawed(SCIP *scip, const GRAPH *graph, const EXTDATA *extdata)
static SCIP_Real distCloseNodesGetMaxCost(int vertex, const DISTDATA *distdata)
int *RESTRICT tail
Definition: graphdefs.h:223
CSRDEPO *const msts_levelbase
int *const extstack_start
#define flipedge(edge)
Definition: graphdefs.h:84
#define MAX(x, y)
Definition: tclique_def.h:83
SCIP_Real *const redcost_treenodeswaps
void extreduce_printLeaves(const EXTDATA *extdata)
SCIP_Real *const tree_bottleneckDistNode
PCDATA *const pcdata
void graph_edge_printInfo(const GRAPH *, int)
Definition: graph_stats.c:133
SCIP_Bool graph_pc_isPc(const GRAPH *)
#define CONNECT
Definition: graphdefs.h:87
Portable definitions.
#define EXT_STATE_NONE
Definition: extreducedefs.h:52
SCIP_Real * adjedgecosts
Definition: completegraph.h:46
SCIP_Bool extreduce_sdsverticalInSync(SCIP *scip, const GRAPH *graph, int compsize, int nleaves_ancestors, int topleaf, EXTDATA *extdata)
SCIP_Real mstobj
Definition: completegraph.h:74
SCIP_Real * cost
Definition: graphdefs.h:221
const int redcost_nlevels
int extreduce_getMaxStackSize(void)
static SCIP_Bool treeDistsAreFlawed(const GRAPH *graph, const EXTDATA *extdata)
#define STP_EXT_MAXGRAD
Definition: extreducedefs.h:42
SCIP_RETCODE graph_heap_create(SCIP *, int, int *, DENTRY *, DHEAP **)
Definition: graph_util.c:992
SCIP_Bool graph_pc_isPcMw(const GRAPH *)
#define SCIP_Real
Definition: def.h:177
#define BLOCKED
Definition: graphdefs.h:90
int graph_heap_deleteMinReturnNode(DHEAP *)
Definition: graph_util.c:1076
#define SCIPfreeCleanBufferArray(scip, ptr)
Definition: scip_mem.h:137
MLDISTS *const sds_vertical
int *RESTRICT outbeg
Definition: graphdefs.h:204
static SCIP_RETCODE distCloseNodesCompute(SCIP *scip, const GRAPH *g, const DISTDATA *distdata, int startvertex, int *closenodes_indices, SCIP_Real *closenodes_dists, int *nclosenodes)
int graph_get_nNodes(const GRAPH *)
Definition: graph_stats.c:536
SCIP_Bool graph_knot_isInRange(const GRAPH *, int)
Definition: graph_node.c:52
MLDISTS *const sds_horizontal
#define UNKNOWN
Definition: sepa_mcf.c:4095
const int * extreduce_mldistsTopLevelBases(const MLDISTS *)
int *const tree_deg
#define nnodes
Definition: gastrans.c:65
void cgraph_node_applyMinAdjCosts(CGRAPH *cgraph, int nodepos, int nodeid)
static SCIP_Bool distCloseNodesIncluded(SCIP *scip, const GRAPH *g, const DISTDATA *distdata, int vertex, const int *closenodes_indices, const SCIP_Real *closenodes_dists, int nclosenodes)
DCSR * dcsr_storage
Definition: graphdefs.h:271
#define SCIP_CALL_ABORT(x)
Definition: def.h:363
static void distCloseNodesPrintLostNodeInfo(const GRAPH *g, const DISTDATA *distdata, int vertex_base, int index_lost, const int *closenodes_indices, const SCIP_Real *closenodes_dists, int nclosenodes)
SCIP_Bool graph_pseudoAncestors_edgeIsHashed(const PSEUDOANS *, int, const int *)
const int * extreduce_mldistsTopTargetIds(const MLDISTS *, int)
SCIP_Bool * pathroot_isdirty
Definition: extreducedefs.h:91
static void treeResetCounters(const GRAPH *graph, const EXTDATA *extdata, int edgecount[], int degreecount[])
static void mstTopLevelBaseGetNodes(int extnode, int nnodes, int nodes[], const EXTDATA *extdata)
void graph_heap_correct(int, SCIP_Real, DHEAP *)
Definition: graph_util.c:1166
void cmst_computeMst(const CGRAPH *cgraph, int mstroot, CMST *cmst)
SCIP_Real extreduce_distComputeRestrictedDist(SCIP *scip, const GRAPH *g, int vertexBlocked, const DISTDATA *distdata, int vertex1, int vertex2)
int graph_csrdepo_getNcsrs(const CSRDEPO *)
Definition: graph_util.c:741
SCIP_Real *const tree_parentEdgeCost