1
2 // #include <cstdio>
3 #include <cassert>
4 #include <climits>
5 #include <set>
6 #include <list>
7 #include <algorithm>
8
9 #include "defs.hh"
10 #include "graph.hh"
11 #include "partition.hh"
12 #include "utils.hh"
13
14 /* use 'and' instead of '&&' */
15 #if _MSC_VER
16 #include <ciso646>
17 #endif
18
19 /*
20 Copyright (c) 2003-2015 Tommi Junttila
21 Released under the GNU Lesser General Public License version 3.
22
23 This file is part of bliss.
24
25 bliss is free software: you can redistribute it and/or modify
26 it under the terms of the GNU Lesser General Public License as published by
27 the Free Software Foundation, version 3 of the License.
28
29 bliss is distributed in the hope that it will be useful,
30 but WITHOUT ANY WARRANTY; without even the implied warranty of
31 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
32 GNU Lesser General Public License for more details.
33
34 You should have received a copy of the GNU Lesser General Public License
35 along with bliss. If not, see <http://www.gnu.org/licenses/>.
36 */
37
38
39 namespace bliss {
40
41 #define _INTERNAL_ERROR() fatal_error("%s:%d: internal error",__FILE__,__LINE__)
42 #define _OUT_OF_MEMORY() fatal_error("%s:%d: out of memory",__FILE__,__LINE__)
43
44 /*-------------------------------------------------------------------------
45 *
46 * Constructor and destructor routines for the abstract graph class
47 *
48 *-------------------------------------------------------------------------*/
49
50
AbstractGraph()51 AbstractGraph::AbstractGraph()
52 {
53 /* Initialize stuff */
54 first_path_labeling = 0;
55 first_path_labeling_inv = 0;
56 best_path_labeling = 0;
57 best_path_labeling_inv = 0;
58 first_path_automorphism = 0;
59 best_path_automorphism = 0;
60 in_search = false;
61 refine_equal_to_first = false;
62
63 /* Default value for using "long prune" */
64 opt_use_long_prune = true;
65 /* Default value for using failure recording */
66 opt_use_failure_recording = true;
67 /* Default value for using component recursion */
68 opt_use_comprec = true;
69
70
71 verbose_level = 0;
72 // verbstr = stdout;
73
74 report_hook = 0;
75 report_user_param = 0;
76 }
77
78
~AbstractGraph()79 AbstractGraph::~AbstractGraph()
80 {
81 if(first_path_labeling) {
82 free(first_path_labeling); first_path_labeling = 0; }
83 if(first_path_labeling_inv) {
84 free(first_path_labeling_inv); first_path_labeling_inv = 0; }
85 if(best_path_labeling) {
86 free(best_path_labeling); best_path_labeling = 0; }
87 if(best_path_labeling_inv) {
88 free(best_path_labeling_inv); best_path_labeling_inv = 0; }
89 if(first_path_automorphism) {
90 free(first_path_automorphism); first_path_automorphism = 0; }
91 if(best_path_automorphism) {
92 free(best_path_automorphism); best_path_automorphism = 0; }
93
94 report_hook = 0;
95 report_user_param = 0;
96 }
97
98
99
100 /*-------------------------------------------------------------------------
101 *
102 * Verbose output management routines
103 *
104 *-------------------------------------------------------------------------*/
105
106 void
set_verbose_level(const unsigned int level)107 AbstractGraph::set_verbose_level(const unsigned int level)
108 {
109 verbose_level = level;
110 }
111
112 void
set_verbose_file(FILE * const fp)113 AbstractGraph::set_verbose_file(FILE* const fp)
114 {
115 // verbstr = fp;
116 }
117
118
119
120 /*-------------------------------------------------------------------------
121 *
122 * Routines for refinement to equitable partition
123 *
124 *-------------------------------------------------------------------------*/
125
126 void
refine_to_equitable()127 AbstractGraph::refine_to_equitable()
128 {
129
130 /* Start refinement from all cells -> push 'em all in the splitting queue */
131 for(Partition::Cell* cell = p.first_cell; cell; cell = cell->next)
132 p.splitting_queue_add(cell);
133
134 do_refine_to_equitable();
135
136 }
137
138 void
refine_to_equitable(Partition::Cell * const unit_cell)139 AbstractGraph::refine_to_equitable(Partition::Cell* const unit_cell)
140 {
141
142 p.splitting_queue_add(unit_cell);
143
144 do_refine_to_equitable();
145 }
146
147
148
149 void
refine_to_equitable(Partition::Cell * const unit_cell1,Partition::Cell * const unit_cell2)150 AbstractGraph::refine_to_equitable(Partition::Cell* const unit_cell1,
151 Partition::Cell* const unit_cell2)
152 {
153
154 p.splitting_queue_add(unit_cell1);
155 p.splitting_queue_add(unit_cell2);
156
157 do_refine_to_equitable();
158 }
159
160
161
162 bool
do_refine_to_equitable()163 AbstractGraph::do_refine_to_equitable()
164 {
165
166 eqref_hash.reset();
167
168 while(!p.splitting_queue_is_empty())
169 {
170 Partition::Cell* const cell = p.splitting_queue_pop();
171
172 if(cell->is_unit())
173 {
174 if(in_search) {
175 const unsigned int index = cell->first;
176 if(first_path_automorphism)
177 {
178 /* Build the (potential) automorphism on-the-fly */
179 first_path_automorphism[first_path_labeling_inv[index]] =
180 p.elements[index];
181 }
182 if(best_path_automorphism)
183 {
184 /* Build the (potential) automorphism on-the-fly */
185 best_path_automorphism[best_path_labeling_inv[index]] =
186 p.elements[index];
187 }
188 }
189 const bool worse = split_neighbourhood_of_unit_cell(cell);
190 if(in_search and worse)
191 goto worse_exit;
192 }
193 else
194 {
195 const bool worse = split_neighbourhood_of_cell(cell);
196 if(in_search and worse)
197 goto worse_exit;
198 }
199 }
200
201 return true;
202
203 worse_exit:
204 /* Clear splitting_queue */
205 p.splitting_queue_clear();
206 return false;
207 }
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224 /*-------------------------------------------------------------------------
225 *
226 * Routines for handling the canonical labeling
227 *
228 *-------------------------------------------------------------------------*/
229
230 /** \internal
231 * Assign the labeling induced by the current partition 'this.p' to
232 * \a labeling.
233 * That is, if the partition is [[2,0],[1]],
234 * then \a labeling will map 0 to 1, 1 to 2, and 2 to 0.
235 */
236 void
update_labeling(unsigned int * const labeling)237 AbstractGraph::update_labeling(unsigned int* const labeling)
238 {
239 const unsigned int N = get_nof_vertices();
240 unsigned int* ep = p.elements;
241 for(unsigned int i = 0; i < N; i++, ep++)
242 labeling[*ep] = i;
243 }
244
245 /** \internal
246 * The same as update_labeling() except that the inverse of the labeling
247 * is also produced and assigned to \a labeling_inv.
248 */
249 void
update_labeling_and_its_inverse(unsigned int * const labeling,unsigned int * const labeling_inv)250 AbstractGraph::update_labeling_and_its_inverse(unsigned int* const labeling,
251 unsigned int* const labeling_inv)
252 {
253 const unsigned int N = get_nof_vertices();
254 unsigned int* ep = p.elements;
255 unsigned int* clip = labeling_inv;
256
257 for(unsigned int i = 0; i < N; ) {
258 labeling[*ep] = i;
259 i++;
260 *clip = *ep;
261 ep++;
262 clip++;
263 }
264 }
265
266
267
268
269
270 /*-------------------------------------------------------------------------
271 *
272 * Routines for handling automorphisms
273 *
274 *-------------------------------------------------------------------------*/
275
276
277 /** \internal
278 * Reset the permutation \a perm to the identity permutation.
279 */
280 void
reset_permutation(unsigned int * perm)281 AbstractGraph::reset_permutation(unsigned int* perm)
282 {
283 const unsigned int N = get_nof_vertices();
284 for(unsigned int i = 0; i < N; i++, perm++)
285 *perm = i;
286 }
287
288 bool
is_automorphism(unsigned int * const perm)289 AbstractGraph::is_automorphism(unsigned int* const perm)
290 {
291 _INTERNAL_ERROR();
292 return false;
293 }
294
295 bool
is_automorphism(const std::vector<unsigned int> & perm) const296 AbstractGraph::is_automorphism(const std::vector<unsigned int>& perm) const
297 {
298 _INTERNAL_ERROR();
299 return false;
300 }
301
302
303
304
305 /*-------------------------------------------------------------------------
306 *
307 * Certificate building
308 *
309 *-------------------------------------------------------------------------*/
310
311 void
cert_add(const unsigned int v1,const unsigned int v2,const unsigned int v3)312 AbstractGraph::cert_add(const unsigned int v1,
313 const unsigned int v2,
314 const unsigned int v3)
315 {
316 if(refine_compare_certificate)
317 {
318 if(refine_equal_to_first)
319 {
320 /* So far equivalent to the first path... */
321 unsigned int index = certificate_current_path.size();
322 if(index >= refine_first_path_subcertificate_end)
323 {
324 refine_equal_to_first = false;
325 }
326 else if(certificate_first_path[index] != v1)
327 {
328 refine_equal_to_first = false;
329 }
330 else if(certificate_first_path[++index] != v2)
331 {
332 refine_equal_to_first = false;
333 }
334 else if(certificate_first_path[++index] != v3)
335 {
336 refine_equal_to_first = false;
337 }
338 if(opt_use_failure_recording and !refine_equal_to_first)
339 {
340 /* We just became different from the first path,
341 * remember the deviation point tree-specific invariant
342 * for the use of failure recording */
343 UintSeqHash h;
344 h.update(v1);
345 h.update(v2);
346 h.update(v3);
347 h.update(index);
348 h.update(eqref_hash.get_value());
349 failure_recording_fp_deviation = h.get_value();
350 }
351 }
352 if(refine_cmp_to_best == 0)
353 {
354 /* So far equivalent to the current best path... */
355 unsigned int index = certificate_current_path.size();
356 if(index >= refine_best_path_subcertificate_end)
357 {
358 refine_cmp_to_best = 1;
359 }
360 else if(v1 > certificate_best_path[index])
361 {
362 refine_cmp_to_best = 1;
363 }
364 else if(v1 < certificate_best_path[index])
365 {
366 refine_cmp_to_best = -1;
367 }
368 else if(v2 > certificate_best_path[++index])
369 {
370 refine_cmp_to_best = 1;
371 }
372 else if(v2 < certificate_best_path[index])
373 {
374 refine_cmp_to_best = -1;
375 }
376 else if(v3 > certificate_best_path[++index])
377 {
378 refine_cmp_to_best = 1;
379 }
380 else if(v3 < certificate_best_path[index])
381 {
382 refine_cmp_to_best = -1;
383 }
384 }
385 if((refine_equal_to_first == false) and
386 (refine_cmp_to_best < 0))
387 return;
388 }
389 /* Update the current path certificate */
390 certificate_current_path.push_back(v1);
391 certificate_current_path.push_back(v2);
392 certificate_current_path.push_back(v3);
393 }
394
395
396 void
cert_add_redundant(const unsigned int v1,const unsigned int v2,const unsigned int v3)397 AbstractGraph::cert_add_redundant(const unsigned int v1,
398 const unsigned int v2,
399 const unsigned int v3)
400 {
401 return cert_add(v1, v2, v3);
402 }
403
404
405
406
407
408
409
410
411
412
413
414 /*-------------------------------------------------------------------------
415 *
416 * Long prune code
417 *
418 *-------------------------------------------------------------------------*/
419
420 void
long_prune_init()421 AbstractGraph::long_prune_init()
422 {
423 const unsigned int N = get_nof_vertices();
424 long_prune_temp.clear();
425 long_prune_temp.resize(N);
426 /* Of how many automorphisms we can store information in
427 the predefined, fixed amount of memory? */
428 const unsigned int nof_fitting_in_max_mem =
429 (long_prune_options_max_mem * 1024 * 1024) / (((N * 2) / 8)+1);
430 long_prune_max_stored_autss = long_prune_options_max_stored_auts;
431 /* Had some problems with g++ in using (a<b)?a:b when constants involved,
432 so had to make this in a stupid way... */
433 if(nof_fitting_in_max_mem < long_prune_options_max_stored_auts)
434 long_prune_max_stored_autss = nof_fitting_in_max_mem;
435
436 long_prune_deallocate();
437 long_prune_fixed.resize(N, 0);
438 long_prune_mcrs.resize(N, 0);
439 long_prune_begin = 0;
440 long_prune_end = 0;
441 }
442
443 void
long_prune_deallocate()444 AbstractGraph::long_prune_deallocate()
445 {
446 while(!long_prune_fixed.empty())
447 {
448 delete long_prune_fixed.back();
449 long_prune_fixed.pop_back();
450 }
451 while(!long_prune_mcrs.empty())
452 {
453 delete long_prune_mcrs.back();
454 long_prune_mcrs.pop_back();
455 }
456 }
457
458 void
long_prune_swap(const unsigned int i,const unsigned int j)459 AbstractGraph::long_prune_swap(const unsigned int i, const unsigned int j)
460 {
461 const unsigned int real_i = i % long_prune_max_stored_autss;
462 const unsigned int real_j = j % long_prune_max_stored_autss;
463 std::vector<bool>* tmp = long_prune_fixed[real_i];
464 long_prune_fixed[real_i] = long_prune_fixed[real_j];
465 long_prune_fixed[real_j] = tmp;
466 tmp = long_prune_mcrs[real_i];
467 long_prune_mcrs[real_i] = long_prune_mcrs[real_j];
468 long_prune_mcrs[real_j] = tmp;
469 }
470
471 std::vector<bool>&
long_prune_allocget_fixed(const unsigned int index)472 AbstractGraph::long_prune_allocget_fixed(const unsigned int index)
473 {
474 const unsigned int i = index % long_prune_max_stored_autss;
475 if(!long_prune_fixed[i])
476 long_prune_fixed[i] = new std::vector<bool>(get_nof_vertices());
477 return *long_prune_fixed[i];
478 }
479
480 std::vector<bool>&
long_prune_get_fixed(const unsigned int index)481 AbstractGraph::long_prune_get_fixed(const unsigned int index)
482 {
483 return *long_prune_fixed[index % long_prune_max_stored_autss];
484 }
485
486 std::vector<bool>&
long_prune_allocget_mcrs(const unsigned int index)487 AbstractGraph::long_prune_allocget_mcrs(const unsigned int index)
488 {
489 const unsigned int i = index % long_prune_max_stored_autss;
490 if(!long_prune_mcrs[i])
491 long_prune_mcrs[i] = new std::vector<bool>(get_nof_vertices());
492 return *long_prune_mcrs[i];
493 }
494
495 std::vector<bool>&
long_prune_get_mcrs(const unsigned int index)496 AbstractGraph::long_prune_get_mcrs(const unsigned int index)
497 {
498 return *long_prune_mcrs[index % long_prune_max_stored_autss];
499 }
500
501 void
long_prune_add_automorphism(const unsigned int * aut)502 AbstractGraph::long_prune_add_automorphism(const unsigned int* aut)
503 {
504 if(long_prune_max_stored_autss == 0)
505 return;
506
507 const unsigned int N = get_nof_vertices();
508
509
510 /* If the buffer of stored auts is full, remove the oldest aut */
511 if(long_prune_end - long_prune_begin == long_prune_max_stored_autss)
512 {
513 long_prune_begin++;
514 }
515 long_prune_end++;
516 std::vector<bool>& fixed = long_prune_allocget_fixed(long_prune_end-1);
517 std::vector<bool>& mcrs = long_prune_allocget_mcrs(long_prune_end-1);
518 /* Mark nodes that are (i) fixed or (ii) minimal orbit representatives
519 * under the automorphism 'aut' */
520 for(unsigned int i = 0; i < N; i++)
521 {
522 fixed[i] = (aut[i] == i);
523 if(long_prune_temp[i] == false)
524 {
525 mcrs[i] = true;
526 unsigned int j = aut[i];
527 while(j != i)
528 {
529 long_prune_temp[j] = true;
530 j = aut[j];
531 }
532 }
533 else
534 {
535 mcrs[i] = false;
536 }
537 /* Clear the temp array on-the-fly... */
538 long_prune_temp[i] = false;
539 }
540
541
542 }
543
544
545
546 /*-------------------------------------------------------------------------
547 *
548 * Routines for handling orbit information
549 *
550 *-------------------------------------------------------------------------*/
551
552 void
update_orbit_information(Orbit & o,const unsigned int * perm)553 AbstractGraph::update_orbit_information(Orbit& o, const unsigned int* perm)
554 {
555 const unsigned int N = get_nof_vertices();
556 for(unsigned int i = 0; i < N; i++)
557 if(perm[i] != i)
558 o.merge_orbits(i, perm[i]);
559 }
560
561
562
563
564
565
566
567
568 /*-------------------------------------------------------------------------
569 *
570 * The actual backtracking search
571 *
572 *-------------------------------------------------------------------------*/
573
574 class TreeNode
575 {
576 //friend class AbstractGraph;
577 public:
578 unsigned int split_cell_first;
579
580 int split_element;
581 static const int SPLIT_START = -1;
582 static const int SPLIT_END = -2;
583
584 Partition::BacktrackPoint partition_bt_point;
585
586 unsigned int certificate_index;
587
588 static const char NO = -1;
589 static const char MAYBE = 0;
590 static const char YES = 1;
591
592 /* First path stuff */
593 bool fp_on;
594 bool fp_cert_equal;
595 char fp_extendable;
596
597 /* Best path stuff */
598 bool in_best_path;
599 int cmp_to_best_path;
600
601 unsigned int failure_recording_ival;
602
603 /* Component recursion related data */
604 unsigned int cr_cep_stack_size;
605 unsigned int cr_cep_index;
606 unsigned int cr_level;
607
608 bool needs_long_prune;
609 unsigned int long_prune_begin;
610 std::set<unsigned int, std::less<unsigned int> > long_prune_redundant;
611
612 UintSeqHash eqref_hash;
613 unsigned int subcertificate_length;
614 };
615
616
617
618
619 typedef struct {
620 unsigned int splitting_element;
621 unsigned int certificate_index;
622 unsigned int subcertificate_length;
623 UintSeqHash eqref_hash;
624 } PathInfo;
625
626
627 void
search(const bool canonical,Stats & stats)628 AbstractGraph::search(const bool canonical, Stats& stats)
629 {
630 const unsigned int N = get_nof_vertices();
631
632 unsigned int all_same_level = UINT_MAX;
633
634 p.graph = this;
635
636 /*
637 * Must be done!
638 */
639 remove_duplicate_edges();
640
641 /*
642 * Reset search statistics
643 */
644 stats.reset();
645 stats.nof_nodes = 1;
646 stats.nof_leaf_nodes = 1;
647
648 /* Free old first path data structures */
649 if(first_path_labeling) {
650 free(first_path_labeling); first_path_labeling = 0; }
651 if(first_path_labeling_inv) {
652 free(first_path_labeling_inv); first_path_labeling_inv = 0; }
653 if(first_path_automorphism) {
654 free(first_path_automorphism); first_path_automorphism = 0; }
655
656 /* Free old best path data structures */
657 if(best_path_labeling) {
658 free(best_path_labeling); best_path_labeling = 0; }
659 if(best_path_labeling_inv) {
660 free(best_path_labeling_inv); best_path_labeling_inv = 0; }
661 if(best_path_automorphism) {
662 free(best_path_automorphism); best_path_automorphism = 0; }
663
664 if(N == 0)
665 {
666 /* Nothing to do, return... */
667 return;
668 }
669
670 /* Initialize the partition ... */
671 p.init(N);
672 /* ... and the component recursion data structures in the partition */
673 if(opt_use_comprec)
674 p.cr_init();
675
676 neighbour_heap.init(N);
677
678 in_search = false;
679 /* Do not compute certificate when building the initial partition */
680 refine_compare_certificate = false;
681 /* The 'eqref_hash' hash value is not computed when building
682 * the initial partition as it is not used for anything at the moment.
683 * This saves some cycles. */
684 compute_eqref_hash = false;
685
686 make_initial_equitable_partition();
687
688 /*
689 * Allocate space for the "first path" and "best path" labelings
690 */
691 if(first_path_labeling) free(first_path_labeling);
692 first_path_labeling = (unsigned int*)calloc(N, sizeof(unsigned int));
693 if(!first_path_labeling) _OUT_OF_MEMORY();
694 if(best_path_labeling) free(best_path_labeling);
695 best_path_labeling = (unsigned int*)calloc(N, sizeof(unsigned int));
696 if(!best_path_labeling) _OUT_OF_MEMORY();
697
698 /*
699 * Is the initial partition discrete?
700 */
701 if(p.is_discrete())
702 {
703 /* Make the best path labeling i.e. the canonical labeling */
704 update_labeling(best_path_labeling);
705 /* Update statistics */
706 stats.nof_leaf_nodes = 1;
707 /* Free component recursion data */
708 if(opt_use_comprec)
709 p.cr_free();
710 return;
711 }
712
713 /*
714 * Allocate the inverses of the "first path" and "best path" labelings
715 */
716 if(first_path_labeling_inv) free(first_path_labeling_inv);
717 first_path_labeling_inv = (unsigned int*)calloc(N, sizeof(unsigned int));
718 if(!first_path_labeling_inv) _OUT_OF_MEMORY();
719 if(best_path_labeling_inv) free(best_path_labeling_inv);
720 best_path_labeling_inv = (unsigned int*)calloc(N, sizeof(unsigned int));
721 if(!best_path_labeling_inv) _OUT_OF_MEMORY();
722
723 /*
724 * Allocate space for the automorphisms
725 */
726 if(first_path_automorphism) free(first_path_automorphism);
727 first_path_automorphism = (unsigned int*)malloc(N * sizeof(unsigned int));
728 if(!first_path_automorphism) _OUT_OF_MEMORY();
729 if(best_path_automorphism) free(best_path_automorphism);
730 best_path_automorphism = (unsigned int*)malloc(N * sizeof(unsigned int));
731 if(!best_path_automorphism) _OUT_OF_MEMORY();
732
733 /*
734 * Initialize orbit information so that all vertices are in their own orbits
735 */
736 first_path_orbits.init(N);
737 best_path_orbits.init(N);
738
739 /*
740 * Initialize certificate memory
741 */
742 initialize_certificate();
743
744 std::vector<TreeNode> search_stack;
745 std::vector<PathInfo> first_path_info;
746 std::vector<PathInfo> best_path_info;
747
748 search_stack.clear();
749
750 /* Initialize "long prune" data structures */
751 if(opt_use_long_prune)
752 long_prune_init();
753
754 /*
755 * Initialize failure recording data structures
756 */
757 typedef std::set<unsigned int, std::less<unsigned int> > FailureRecordingSet;
758 std::vector<FailureRecordingSet> failure_recording_hashes;
759
760 /*
761 * Initialize component recursion data structures
762 */
763 cr_cep_stack.clear();
764 unsigned int cr_cep_index = 0;
765 {
766 /* Inset a sentinel "component end point" */
767 CR_CEP sentinel;
768 sentinel.creation_level = 0;
769 sentinel.discrete_cell_limit = get_nof_vertices();
770 sentinel.next_cr_level = 0;
771 sentinel.next_cep_index = 0;
772 sentinel.first_checked = false;
773 sentinel.best_checked = false;
774 cr_cep_index = 0;
775 cr_cep_stack.push_back(sentinel);
776 }
777 cr_level = 0;
778 if(opt_use_comprec and
779 nucr_find_first_component(cr_level) == true and
780 p.nof_discrete_cells() + cr_component_elements <
781 cr_cep_stack[cr_cep_index].discrete_cell_limit)
782 {
783 cr_level = p.cr_split_level(0, cr_component);
784 CR_CEP cep;
785 cep.creation_level = 0;
786 cep.discrete_cell_limit = p.nof_discrete_cells() + cr_component_elements;
787 cep.next_cr_level = 0;
788 cep.next_cep_index = cr_cep_index;
789 cep.first_checked = false;
790 cep.best_checked = false;
791 cr_cep_index = cr_cep_stack.size();
792 cr_cep_stack.push_back(cep);
793 }
794
795 /*
796 * Build the root node of the search tree
797 */
798 {
799 TreeNode root;
800 Partition::Cell* split_cell = find_next_cell_to_be_splitted(p.first_cell);
801 root.split_cell_first = split_cell->first;
802 root.split_element = TreeNode::SPLIT_START;
803 root.partition_bt_point = p.set_backtrack_point();
804 root.certificate_index = 0;
805 root.fp_on = true;
806 root.fp_cert_equal = true;
807 root.fp_extendable = TreeNode::MAYBE;
808 root.in_best_path = false;
809 root.cmp_to_best_path = 0;
810 root.long_prune_begin = 0;
811 root.needs_long_prune = false;
812
813 root.failure_recording_ival = 0;
814
815 /* Save component recursion info for backtracking */
816 root.cr_level = cr_level;
817 root.cr_cep_stack_size = cr_cep_stack.size();
818 root.cr_cep_index = cr_cep_index;
819 search_stack.push_back(root);
820 }
821
822 /*
823 * Set status and global flags for search related procedures
824 */
825 in_search = true;
826 /* Do not compare certificates during refinement until the first path has been traversed to the leaf */
827 refine_compare_certificate = false;
828
829
830
831
832 /*
833 * The actual backtracking search
834 */
835 while(!search_stack.empty())
836 {
837 TreeNode& current_node = search_stack.back();
838 const unsigned int current_level = (unsigned int)search_stack.size()-1;
839
840
841 if(opt_use_comprec)
842 {
843 CR_CEP& cep = cr_cep_stack[current_node.cr_cep_index];
844 if(cep.first_checked == true and
845 current_node.fp_extendable == TreeNode::MAYBE and
846 !search_stack[cep.creation_level].fp_on)
847 {
848 current_node.fp_extendable = TreeNode::NO;
849 }
850 }
851
852 if(current_node.fp_on)
853 {
854 if(current_node.split_element == TreeNode::SPLIT_END)
855 {
856 search_stack.pop_back();
857 continue;
858 }
859 }
860 else
861 {
862 if(current_node.fp_extendable == TreeNode::YES)
863 {
864 search_stack.pop_back();
865 continue;
866 }
867 if(current_node.split_element == TreeNode::SPLIT_END)
868 {
869 if(opt_use_failure_recording)
870 {
871 TreeNode& parent_node = search_stack[current_level-1];
872 if(parent_node.fp_on)
873 failure_recording_hashes[current_level-1].insert(current_node.failure_recording_ival);
874 }
875 search_stack.pop_back();
876 continue;
877 }
878 if(current_node.fp_extendable == TreeNode::NO and
879 (!canonical or current_node.cmp_to_best_path < 0))
880 {
881 if(opt_use_failure_recording)
882 {
883 TreeNode& parent_node = search_stack[current_level-1];
884 if(parent_node.fp_on)
885 failure_recording_hashes[current_level-1].insert(current_node.failure_recording_ival);
886 }
887 search_stack.pop_back();
888 continue;
889 }
890 }
891
892 /* Restore partition ... */
893 p.goto_backtrack_point(current_node.partition_bt_point);
894 /* ... and re-remember backtracking point */
895 current_node.partition_bt_point = p.set_backtrack_point();
896
897 /* Restore current path certificate */
898 certificate_index = current_node.certificate_index;
899 refine_current_path_certificate_index = current_node.certificate_index;
900 certificate_current_path.resize(certificate_index);
901
902 /* Fetch split cell information */
903 Partition::Cell * const cell =
904 p.get_cell(p.elements[current_node.split_cell_first]);
905
906 /* Restore component recursion information */
907 cr_level = current_node.cr_level;
908 cr_cep_stack.resize(current_node.cr_cep_stack_size);
909 cr_cep_index = current_node.cr_cep_index;
910
911
912 /*
913 * Update long prune redundancy sets
914 */
915 if(opt_use_long_prune and current_level >= 1 and !current_node.fp_on)
916 {
917 unsigned int begin = (current_node.long_prune_begin>long_prune_begin)?current_node.long_prune_begin:long_prune_begin;
918 for(unsigned int i = begin; i < long_prune_end; i++)
919 {
920 const std::vector<bool>& fixed = long_prune_get_fixed(i);
921 #if defined(BLISS_CONSISTENCY_CHECKS)
922 for(unsigned int l = 0; l < search_stack.size()-2; l++)
923 assert(fixed[search_stack[l].split_element]);
924 #endif
925 if(fixed[search_stack[search_stack.size()-1-1].split_element] ==
926 false)
927 {
928 long_prune_swap(begin, i);
929 begin++;
930 current_node.long_prune_begin = begin;
931 continue;
932 }
933 }
934
935 if(current_node.split_element == TreeNode::SPLIT_START)
936 {
937 current_node.needs_long_prune = true;
938 }
939 else if(current_node.needs_long_prune)
940 {
941 current_node.needs_long_prune = false;
942 unsigned int begin = (current_node.long_prune_begin>long_prune_begin)?current_node.long_prune_begin:long_prune_begin;
943 for(unsigned int i = begin; i < long_prune_end; i++)
944 {
945 const std::vector<bool>& fixed = long_prune_get_fixed(i);
946 #if defined(BLISS_CONSISTENCY_CHECKS)
947 for(unsigned int l = 0; l < search_stack.size()-2; l++)
948 assert(fixed[search_stack[l].split_element]);
949 #endif
950 assert(fixed[search_stack[current_level-1].split_element] == true);
951 if(fixed[search_stack[current_level-1].split_element] == false)
952 {
953 long_prune_swap(begin, i);
954 begin++;
955 current_node.long_prune_begin = begin;
956 continue;
957 }
958 const std::vector<bool>& mcrs = long_prune_get_mcrs(i);
959 unsigned int* ep = p.elements + cell->first;
960 for(unsigned int j = cell->length; j > 0; j--, ep++) {
961 if(mcrs[*ep] == false)
962 current_node.long_prune_redundant.insert(*ep);
963 }
964 }
965 }
966 }
967
968
969 /*
970 * Find the next smallest, non-isomorphic element in the cell and
971 * store it in current_node.split_element
972 */
973 {
974 unsigned int next_split_element = UINT_MAX;
975 //unsigned int* next_split_element_pos = 0;
976 unsigned int* ep = p.elements + cell->first;
977 if(current_node.fp_on)
978 {
979 /* Find the next larger splitting element that is
980 * a minimal orbit representative w.r.t. first_path_orbits */
981 for(unsigned int i = cell->length; i > 0; i--, ep++) {
982 if((int)(*ep) > current_node.split_element and
983 *ep < next_split_element and
984 first_path_orbits.is_minimal_representative(*ep)) {
985 next_split_element = *ep;
986 //next_split_element_pos = ep;
987 }
988 }
989 }
990 else if(current_node.in_best_path)
991 {
992 /* Find the next larger splitting element that is
993 * a minimal orbit representative w.r.t. best_path_orbits */
994 for(unsigned int i = cell->length; i > 0; i--, ep++) {
995 if((int)(*ep) > current_node.split_element and
996 *ep < next_split_element and
997 best_path_orbits.is_minimal_representative(*ep) and
998 (!opt_use_long_prune or
999 current_node.long_prune_redundant.find(*ep) ==
1000 current_node.long_prune_redundant.end())) {
1001 next_split_element = *ep;
1002 //next_split_element_pos = ep;
1003 }
1004 }
1005 }
1006 else
1007 {
1008 /* Find the next larger splitting element */
1009 for(unsigned int i = cell->length; i > 0; i--, ep++) {
1010 if((int)(*ep) > current_node.split_element and
1011 *ep < next_split_element and
1012 (!opt_use_long_prune or
1013 current_node.long_prune_redundant.find(*ep) ==
1014 current_node.long_prune_redundant.end())) {
1015 next_split_element = *ep;
1016 //next_split_element_pos = ep;
1017 }
1018 }
1019 }
1020 if(next_split_element == UINT_MAX)
1021 {
1022 /* No more (unexplored children) in the cell */
1023 current_node.split_element = TreeNode::SPLIT_END;
1024 if(current_node.fp_on)
1025 {
1026 /* Update group size */
1027 const unsigned int index = first_path_orbits.orbit_size(first_path_info[search_stack.size()-1].splitting_element);
1028 stats.group_size.multiply(index);
1029 stats.group_size_approx *= (long double)index;
1030 /*
1031 * Update all_same_level
1032 */
1033 if(index == cell->length and all_same_level == current_level+1)
1034 all_same_level = current_level;
1035 }
1036 continue;
1037 }
1038
1039 /* Split on smallest */
1040 current_node.split_element = next_split_element;
1041 }
1042
1043 const unsigned int child_level = current_level+1;
1044 /* Update some statistics */
1045 stats.nof_nodes++;
1046 if(search_stack.size() > stats.max_level)
1047 stats.max_level = search_stack.size();
1048
1049
1050
1051 /* Set flags and indices for the refiner certificate builder */
1052 refine_equal_to_first = current_node.fp_cert_equal;
1053 refine_cmp_to_best = current_node.cmp_to_best_path;
1054 if(!first_path_info.empty())
1055 {
1056 if(refine_equal_to_first)
1057 refine_first_path_subcertificate_end =
1058 first_path_info[search_stack.size()-1].certificate_index +
1059 first_path_info[search_stack.size()-1].subcertificate_length;
1060 if(canonical)
1061 {
1062 if(refine_cmp_to_best == 0)
1063 refine_best_path_subcertificate_end =
1064 best_path_info[search_stack.size()-1].certificate_index +
1065 best_path_info[search_stack.size()-1].subcertificate_length;
1066 }
1067 else
1068 refine_cmp_to_best = -1;
1069 }
1070
1071 const bool was_fp_cert_equal = current_node.fp_cert_equal;
1072
1073 /* Individualize, i.e. split the cell in two, the latter new cell
1074 * will be a unit one containing info.split_element */
1075 Partition::Cell* const new_cell =
1076 p.individualize(cell, current_node.split_element);
1077
1078 /*
1079 * Refine the new partition to equitable
1080 */
1081 if(cell->is_unit())
1082 refine_to_equitable(cell, new_cell);
1083 else
1084 refine_to_equitable(new_cell);
1085
1086
1087
1088
1089 /* Update statistics */
1090 if(p.is_discrete())
1091 stats.nof_leaf_nodes++;
1092
1093
1094 if(!first_path_info.empty())
1095 {
1096 /* We are no longer on the first path */
1097 const unsigned int subcertificate_length =
1098 certificate_current_path.size() - certificate_index;
1099 if(refine_equal_to_first)
1100 {
1101 /* Was equal to the first path so far */
1102 PathInfo& first_pinfo = first_path_info[current_level];
1103 assert(first_pinfo.certificate_index == certificate_index);
1104 if(subcertificate_length != first_pinfo.subcertificate_length)
1105 {
1106 refine_equal_to_first = false;
1107 if(opt_use_failure_recording)
1108 failure_recording_fp_deviation = subcertificate_length;
1109 }
1110 else if(first_pinfo.eqref_hash.cmp(eqref_hash) != 0)
1111 {
1112 refine_equal_to_first = false;
1113 if(opt_use_failure_recording)
1114 failure_recording_fp_deviation = eqref_hash.get_value();
1115 }
1116 }
1117 if(canonical and (refine_cmp_to_best == 0))
1118 {
1119 /* Was equal to the best path so far */
1120 PathInfo& bestp_info = best_path_info[current_level];
1121 assert(bestp_info.certificate_index == certificate_index);
1122 if(subcertificate_length < bestp_info.subcertificate_length)
1123 {
1124 refine_cmp_to_best = -1;
1125 }
1126 else if(subcertificate_length > bestp_info.subcertificate_length)
1127 {
1128 refine_cmp_to_best = 1;
1129 }
1130 else if(bestp_info.eqref_hash.cmp(eqref_hash) > 0)
1131 {
1132 refine_cmp_to_best = -1;
1133 }
1134 else if(bestp_info.eqref_hash.cmp(eqref_hash) < 0)
1135 {
1136 refine_cmp_to_best = 1;
1137 }
1138 }
1139
1140 if(opt_use_failure_recording and
1141 was_fp_cert_equal and
1142 !refine_equal_to_first)
1143 {
1144 UintSeqHash k;
1145 k.update(failure_recording_fp_deviation);
1146 k.update(eqref_hash.get_value());
1147 failure_recording_fp_deviation = k.get_value();
1148
1149 if(current_node.fp_on)
1150 failure_recording_hashes[current_level].insert(failure_recording_fp_deviation);
1151 else
1152 {
1153 for(unsigned int i = current_level; i > 0; i--)
1154 {
1155 if(search_stack[i].fp_on)
1156 break;
1157 const FailureRecordingSet& s = failure_recording_hashes[i];
1158 if(i == current_level and
1159 s.find(failure_recording_fp_deviation) != s.end())
1160 break;
1161 if(s.find(0) != s.end())
1162 break;
1163 search_stack[i].fp_extendable = TreeNode::NO;
1164 }
1165 }
1166 }
1167
1168
1169 /* Check if no longer equal to the first path and,
1170 * if canonical labeling is desired, also worse than the
1171 * current best path */
1172 if(refine_equal_to_first == false and
1173 (!canonical or (refine_cmp_to_best < 0)))
1174 {
1175 /* Yes, backtrack */
1176 stats.nof_bad_nodes++;
1177 if(current_node.fp_cert_equal == true and
1178 current_level+1 > all_same_level)
1179 {
1180 assert(all_same_level >= 1);
1181 for(unsigned int i = all_same_level;
1182 i < search_stack.size();
1183 i++)
1184 {
1185 search_stack[i].fp_extendable = TreeNode::NO;
1186 }
1187 }
1188
1189 continue;
1190 }
1191 }
1192
1193 #if defined(BLISS_VERIFY_EQUITABLEDNESS)
1194 /* The new partition should be equitable */
1195 if(!is_equitable())
1196 fatal_error("consistency check failed - partition after refinement is not equitable");
1197 #endif
1198
1199 /*
1200 * Next level search tree node info
1201 */
1202 TreeNode child_node;
1203
1204 /* No more in the first path */
1205 child_node.fp_on = false;
1206 /* No more in the best path */
1207 child_node.in_best_path = false;
1208
1209 child_node.fp_cert_equal = refine_equal_to_first;
1210 if(current_node.fp_extendable == TreeNode::NO or
1211 (current_node.fp_extendable == TreeNode::MAYBE and
1212 child_node.fp_cert_equal == false))
1213 child_node.fp_extendable = TreeNode::NO;
1214 else
1215 child_node.fp_extendable = TreeNode::MAYBE;
1216 child_node.cmp_to_best_path = refine_cmp_to_best;
1217
1218 child_node.failure_recording_ival = 0;
1219 child_node.cr_cep_stack_size = current_node.cr_cep_stack_size;
1220 child_node.cr_cep_index = current_node.cr_cep_index;
1221 child_node.cr_level = current_node.cr_level;
1222
1223 certificate_index = certificate_current_path.size();
1224
1225 current_node.eqref_hash = eqref_hash;
1226 current_node.subcertificate_length =
1227 certificate_index - current_node.certificate_index;
1228
1229
1230 /*
1231 * The first encountered leaf node at the end of the "first path"?
1232 */
1233 if(p.is_discrete() and first_path_info.empty())
1234 {
1235 //fprintf(stdout, "Level %u: FIRST\n", child_level); fflush(stdout);
1236 stats.nof_canupdates++;
1237 /*
1238 * Update labelings and their inverses
1239 */
1240 update_labeling_and_its_inverse(first_path_labeling,
1241 first_path_labeling_inv);
1242 update_labeling_and_its_inverse(best_path_labeling,
1243 best_path_labeling_inv);
1244 /*
1245 * Reset automorphism array
1246 */
1247 reset_permutation(first_path_automorphism);
1248 reset_permutation(best_path_automorphism);
1249 /*
1250 * Reset orbit information
1251 */
1252 first_path_orbits.reset();
1253 best_path_orbits.reset();
1254 /*
1255 * Reset group size
1256 */
1257 stats.group_size.assign(1);
1258 stats.group_size_approx = 1.0;
1259 /*
1260 * Reset all_same_level
1261 */
1262 all_same_level = child_level;
1263 /*
1264 * Mark the current path to be the first and best one and save it
1265 */
1266 const unsigned int base_size = search_stack.size();
1267 best_path_info.clear();
1268 //fprintf(stdout, " New base is: ");
1269 for(unsigned int i = 0; i < base_size; i++) {
1270 search_stack[i].fp_on = true;
1271 search_stack[i].fp_cert_equal = true;
1272 search_stack[i].fp_extendable = TreeNode::YES;
1273 search_stack[i].in_best_path = true;
1274 search_stack[i].cmp_to_best_path = 0;
1275 PathInfo path_info;
1276 path_info.splitting_element = search_stack[i].split_element;
1277 path_info.certificate_index = search_stack[i].certificate_index;
1278 path_info.eqref_hash = search_stack[i].eqref_hash;
1279 path_info.subcertificate_length = search_stack[i].subcertificate_length;
1280 first_path_info.push_back(path_info);
1281 best_path_info.push_back(path_info);
1282 //fprintf(stdout, "%u ", search_stack[i].split_element);
1283 }
1284 //fprintf(stdout, "\n"); fflush(stdout);
1285 /* Copy certificates */
1286 certificate_first_path = certificate_current_path;
1287 certificate_best_path = certificate_current_path;
1288
1289 /* From now on, compare certificates when refining */
1290 refine_compare_certificate = true;
1291
1292 if(opt_use_failure_recording)
1293 failure_recording_hashes.resize(base_size);
1294
1295 /*
1296 for(unsigned int j = 0; j < search_stack.size(); j++)
1297 fprintf(stderr, "%u ", search_stack[j].split_element);
1298 fprintf(stderr, "\n");
1299 p.print(stderr); fprintf(stderr, "\n");
1300 */
1301
1302 /*
1303 * Backtrack to the previous level
1304 */
1305 continue;
1306 }
1307
1308
1309 if(p.is_discrete() and child_node.fp_cert_equal)
1310 {
1311 /*
1312 * A leaf node that is equal to the first one.
1313 * An automorphism found: aut[i] = elements[first_path_labeling[i]]
1314 */
1315 goto handle_first_path_automorphism;
1316 }
1317
1318
1319 if(!p.is_discrete())
1320 {
1321 Partition::Cell* next_split_cell = 0;
1322 /*
1323 * An internal, non-leaf node
1324 */
1325 if(opt_use_comprec)
1326 {
1327 assert(p.nof_discrete_cells() <=
1328 cr_cep_stack[cr_cep_index].discrete_cell_limit);
1329 assert(cr_level == child_node.cr_level);
1330
1331
1332 if(p.nof_discrete_cells() ==
1333 cr_cep_stack[cr_cep_index].discrete_cell_limit)
1334 {
1335 /* We have reached the end of a component */
1336 assert(cr_cep_index != 0);
1337 CR_CEP& cep = cr_cep_stack[cr_cep_index];
1338
1339 /* First, compare with respect to the first path */
1340 if(first_path_info.empty() or child_node.fp_cert_equal) {
1341 if(cep.first_checked == false)
1342 {
1343 /* First time, go to the next component */
1344 cep.first_checked = true;
1345 }
1346 else
1347 {
1348 assert(!first_path_info.empty());
1349 assert(cep.creation_level < search_stack.size());
1350 TreeNode& old_info = search_stack[cep.creation_level];
1351 /* If the component was found when on the first path,
1352 * handle the found automorphism as the other
1353 * first path automorphisms */
1354 if(old_info.fp_on)
1355 goto handle_first_path_automorphism;
1356 }
1357 }
1358
1359 if(canonical and
1360 !first_path_info.empty() and
1361 child_node.cmp_to_best_path >= 0) {
1362 if(cep.best_checked == false)
1363 {
1364 /* First time, go to the next component */
1365 cep.best_checked = true;
1366 }
1367 else
1368 {
1369 assert(cep.creation_level < search_stack.size());
1370 TreeNode& old_info = search_stack[cep.creation_level];
1371 if(child_node.cmp_to_best_path == 0) {
1372 /* If the component was found when on the best path,
1373 * handle the found automorphism as the other
1374 * best path automorphisms */
1375 if(old_info.in_best_path)
1376 goto handle_best_path_automorphism;
1377 /* Otherwise, we do not remember the automorhism as
1378 * we didn't memorize the path that was invariant
1379 * equal to the best one and passed through the
1380 * component.
1381 * Thus we can only backtrack to the previous level */
1382 child_node.cmp_to_best_path = -1;
1383 if(!child_node.fp_cert_equal)
1384 {
1385 continue;
1386 }
1387 }
1388 else {
1389 assert(child_node.cmp_to_best_path > 0);
1390 if(old_info.in_best_path)
1391 {
1392 stats.nof_canupdates++;
1393 /*
1394 * Update canonical labeling and its inverse
1395 */
1396 for(unsigned int i = 0; i < N; i++) {
1397 if(p.get_cell(p.elements[i])->is_unit()) {
1398 best_path_labeling[p.elements[i]] = i;
1399 best_path_labeling_inv[i] = p.elements[i];
1400 }
1401 }
1402 //update_labeling_and_its_inverse(best_path_labeling, best_path_labeling_inv);
1403 /* Reset best path automorphism */
1404 reset_permutation(best_path_automorphism);
1405 /* Reset best path orbit structure */
1406 best_path_orbits.reset();
1407 /* Mark to be the best one and save prefix */
1408 unsigned int postfix_start = cep.creation_level;
1409 assert(postfix_start < best_path_info.size());
1410 while(p.get_cell(best_path_info[postfix_start].splitting_element)->is_unit()) {
1411 postfix_start++;
1412 assert(postfix_start < best_path_info.size());
1413 }
1414 unsigned int postfix_start_cert = best_path_info[postfix_start].certificate_index;
1415 std::vector<PathInfo> best_path_temp = best_path_info;
1416 best_path_info.clear();
1417 for(unsigned int i = 0; i < search_stack.size(); i++) {
1418 TreeNode& ss_info = search_stack[i];
1419 PathInfo bp_info;
1420 ss_info.cmp_to_best_path = 0;
1421 ss_info.in_best_path = true;
1422 bp_info.splitting_element = ss_info.split_element;
1423 bp_info.certificate_index = ss_info.certificate_index;
1424 bp_info.subcertificate_length = ss_info.subcertificate_length;
1425 bp_info.eqref_hash = ss_info.eqref_hash;
1426 best_path_info.push_back(bp_info);
1427 }
1428 /* Copy the postfix of the previous best path */
1429 for(unsigned int i = postfix_start;
1430 i < best_path_temp.size();
1431 i++)
1432 {
1433 best_path_info.push_back(best_path_temp[i]);
1434 best_path_info[best_path_info.size()-1].certificate_index =
1435 best_path_info[best_path_info.size()-2].certificate_index +
1436 best_path_info[best_path_info.size()-2].subcertificate_length;
1437 }
1438 std::vector<unsigned int> certificate_best_path_old = certificate_best_path;
1439 certificate_best_path = certificate_current_path;
1440 for(unsigned int i = postfix_start_cert; i < certificate_best_path_old.size(); i++)
1441 certificate_best_path.push_back(certificate_best_path_old[i]);
1442 assert(certificate_best_path.size() == best_path_info.back().certificate_index + best_path_info.back().subcertificate_length);
1443 /* Backtrack to the previous level */
1444 continue;
1445 }
1446 }
1447 }
1448 }
1449
1450 /* No backtracking performed, go to next componenet */
1451 cr_level = cep.next_cr_level;
1452 cr_cep_index = cep.next_cep_index;
1453 }
1454
1455 /* Check if the current component has been split into
1456 * new non-uniformity subcomponents */
1457 //if(nucr_find_first_component(cr_level) == true and
1458 // p.nof_discrete_cells() + cr_component_elements <
1459 // cr_cep_stack[cr_cep_index].discrete_cell_limit)
1460 if(nucr_find_first_component(cr_level, cr_component,
1461 cr_component_elements,
1462 next_split_cell) == true and
1463 p.nof_discrete_cells() + cr_component_elements <
1464 cr_cep_stack[cr_cep_index].discrete_cell_limit)
1465 {
1466 const unsigned int next_cr_level =
1467 p.cr_split_level(cr_level, cr_component);
1468 CR_CEP cep;
1469 cep.creation_level = search_stack.size();
1470 cep.discrete_cell_limit =
1471 p.nof_discrete_cells() + cr_component_elements;
1472 cep.next_cr_level = cr_level;
1473 cep.next_cep_index = cr_cep_index;
1474 cep.first_checked = false;
1475 cep.best_checked = false;
1476 cr_cep_index = cr_cep_stack.size();
1477 cr_cep_stack.push_back(cep);
1478 cr_level = next_cr_level;
1479 }
1480 }
1481
1482
1483 /*
1484 * Build the next node info
1485 */
1486 /* Find the next cell to be splitted */
1487 if(!next_split_cell)
1488 next_split_cell = find_next_cell_to_be_splitted(p.get_cell(p.elements[current_node.split_cell_first]));
1489 //Partition::Cell * const next_split_cell = find_next_cell_to_be_splitted(p.get_cell(p.elements[current_node.split_cell_first]));
1490 child_node.split_cell_first = next_split_cell->first;
1491 child_node.split_element = TreeNode::SPLIT_START;
1492 child_node.certificate_index = certificate_index;
1493 child_node.partition_bt_point = p.set_backtrack_point();
1494 child_node.long_prune_redundant.clear();
1495 child_node.long_prune_begin = current_node.long_prune_begin;
1496
1497 /* Save component recursion info for backtracking */
1498 child_node.cr_level = cr_level;
1499 child_node.cr_cep_stack_size = cr_cep_stack.size();
1500 child_node.cr_cep_index = cr_cep_index;
1501
1502 /* Initialize needs_long_prune to prevent a gcc-ubsan warning */
1503 child_node.needs_long_prune = true;
1504
1505 search_stack.push_back(child_node);
1506 continue;
1507 }
1508
1509 /*
1510 * A leaf node not in the first path or equivalent to the first path
1511 */
1512
1513
1514
1515 if(child_node.cmp_to_best_path > 0)
1516 {
1517 /*
1518 * A new, better representative found
1519 */
1520 //fprintf(stdout, "Level %u: NEW BEST\n", child_level); fflush(stdout);
1521 stats.nof_canupdates++;
1522 /*
1523 * Update canonical labeling and its inverse
1524 */
1525 update_labeling_and_its_inverse(best_path_labeling,
1526 best_path_labeling_inv);
1527 /* Reset best path automorphism */
1528 reset_permutation(best_path_automorphism);
1529 /* Reset best path orbit structure */
1530 best_path_orbits.reset();
1531 /*
1532 * Mark the current path to be the best one and save it
1533 */
1534 const unsigned int base_size = search_stack.size();
1535 assert(current_level+1 == base_size);
1536 best_path_info.clear();
1537 for(unsigned int i = 0; i < base_size; i++) {
1538 search_stack[i].cmp_to_best_path = 0;
1539 search_stack[i].in_best_path = true;
1540 PathInfo path_info;
1541 path_info.splitting_element = search_stack[i].split_element;
1542 path_info.certificate_index = search_stack[i].certificate_index;
1543 path_info.subcertificate_length = search_stack[i].subcertificate_length;
1544 path_info.eqref_hash = search_stack[i].eqref_hash;
1545 best_path_info.push_back(path_info);
1546 }
1547 certificate_best_path = certificate_current_path;
1548 /*
1549 * Backtrack to the previous level
1550 */
1551 continue;
1552 }
1553
1554
1555 handle_best_path_automorphism:
1556 /*
1557 *
1558 * Best path automorphism handling
1559 *
1560 */
1561 {
1562
1563 /*
1564 * Equal to the previous best path
1565 */
1566 if(p.is_discrete())
1567 {
1568 #if defined(BLISS_CONSISTENCY_CHECKS)
1569 /* Verify that the automorphism is correctly built */
1570 for(unsigned int i = 0; i < N; i++)
1571 assert(best_path_automorphism[i] ==
1572 p.elements[best_path_labeling[i]]);
1573 #endif
1574 }
1575 else
1576 {
1577 /* An automorphism that was found before the partition was discrete.
1578 * Set the image of all elements in non-disrete cells accordingly */
1579 for(Partition::Cell* c = p.first_nonsingleton_cell; c;
1580 c = c->next_nonsingleton) {
1581 for(unsigned int i = c->first; i < c->first+c->length; i++)
1582 if(p.get_cell(p.elements[best_path_labeling[p.elements[i]]])->is_unit())
1583 best_path_automorphism[p.elements[best_path_labeling[p.elements[i]]]] = p.elements[i];
1584 else
1585 best_path_automorphism[p.elements[i]] = p.elements[i];
1586 }
1587 }
1588
1589 #if defined(BLISS_VERIFY_AUTOMORPHISMS)
1590 /* Verify that it really is an automorphism */
1591 if(!is_automorphism(best_path_automorphism))
1592 fatal_error("Best path automorhism validation check failed");
1593 #endif
1594
1595 unsigned int gca_level_with_first = 0;
1596 for(unsigned int i = search_stack.size(); i > 0; i--) {
1597 if((int)first_path_info[gca_level_with_first].splitting_element !=
1598 search_stack[gca_level_with_first].split_element)
1599 break;
1600 gca_level_with_first++;
1601 }
1602
1603 unsigned int gca_level_with_best = 0;
1604 for(unsigned int i = search_stack.size(); i > 0; i--) {
1605 if((int)best_path_info[gca_level_with_best].splitting_element !=
1606 search_stack[gca_level_with_best].split_element)
1607 break;
1608 gca_level_with_best++;
1609 }
1610
1611 if(opt_use_long_prune)
1612 {
1613 /* Record automorphism */
1614 long_prune_add_automorphism(best_path_automorphism);
1615 }
1616
1617 /*
1618 * Update orbit information
1619 */
1620 update_orbit_information(best_path_orbits, best_path_automorphism);
1621
1622 /*
1623 * Update orbit information
1624 */
1625 const unsigned int nof_old_orbits = first_path_orbits.nof_orbits();
1626 update_orbit_information(first_path_orbits, best_path_automorphism);
1627 if(nof_old_orbits != first_path_orbits.nof_orbits())
1628 {
1629 /* Some orbits were merged */
1630 /* Report automorphism */
1631 if(report_hook)
1632 (*report_hook)(report_user_param,
1633 get_nof_vertices(),
1634 best_path_automorphism);
1635 /* Update statistics */
1636 stats.nof_generators++;
1637 }
1638
1639 /*
1640 * Compute backjumping level
1641 */
1642 unsigned int backjumping_level = current_level+1-1;
1643 if(!first_path_orbits.is_minimal_representative(search_stack[gca_level_with_first].split_element))
1644 {
1645 backjumping_level = gca_level_with_first;
1646 }
1647 else
1648 {
1649 assert(!best_path_orbits.is_minimal_representative(search_stack[gca_level_with_best].split_element));
1650 backjumping_level = gca_level_with_best;
1651 }
1652 /* Backtrack */
1653 search_stack.resize(backjumping_level + 1);
1654 continue;
1655 }
1656
1657
1658 _INTERNAL_ERROR();
1659
1660
1661 handle_first_path_automorphism:
1662 /*
1663 *
1664 * A first-path automorphism: aut[i] = elements[first_path_labeling[i]]
1665 *
1666 */
1667
1668
1669 if(p.is_discrete())
1670 {
1671 #if defined(BLISS_CONSISTENCY_CHECKS)
1672 /* Verify that the complete automorphism is correctly built */
1673 for(unsigned int i = 0; i < N; i++)
1674 assert(first_path_automorphism[i] ==
1675 p.elements[first_path_labeling[i]]);
1676 #endif
1677 }
1678 else
1679 {
1680 /* An automorphism that was found before the partition was discrete.
1681 * Set the image of all elements in non-disrete cells accordingly */
1682 for(Partition::Cell* c = p.first_nonsingleton_cell; c;
1683 c = c->next_nonsingleton) {
1684 for(unsigned int i = c->first; i < c->first+c->length; i++)
1685 if(p.get_cell(p.elements[first_path_labeling[p.elements[i]]])->is_unit())
1686 first_path_automorphism[p.elements[first_path_labeling[p.elements[i]]]] = p.elements[i];
1687 else
1688 first_path_automorphism[p.elements[i]] = p.elements[i];
1689 }
1690 }
1691
1692 #if defined(BLISS_VERIFY_AUTOMORPHISMS)
1693 /* Verify that it really is an automorphism */
1694 if(!is_automorphism(first_path_automorphism))
1695 fatal_error("First path automorphism validation check failed");
1696 #endif
1697
1698 if(opt_use_long_prune)
1699 {
1700 long_prune_add_automorphism(first_path_automorphism);
1701 }
1702
1703 /*
1704 * Update orbit information
1705 */
1706 update_orbit_information(first_path_orbits, first_path_automorphism);
1707
1708 /*
1709 * Compute backjumping level
1710 */
1711 for(unsigned int i = 0; i < search_stack.size(); i++) {
1712 TreeNode& n = search_stack[i];
1713 if(n.fp_on) {
1714 ;
1715 } else {
1716 n.fp_extendable = TreeNode::YES;
1717 }
1718 }
1719
1720 /* Report automorphism by calling the user defined hook function */
1721 if(report_hook)
1722 (*report_hook)(report_user_param,
1723 get_nof_vertices(),
1724 first_path_automorphism);
1725
1726 /* Update statistics */
1727 stats.nof_generators++;
1728 continue;
1729
1730 } /* while(!search_stack.empty()) */
1731
1732
1733
1734
1735 /* Free "long prune" technique memory */
1736 if(opt_use_long_prune)
1737 long_prune_deallocate();
1738
1739 /* Release component recursion data in partition */
1740 if(opt_use_comprec)
1741 p.cr_free();
1742 }
1743
1744
1745
1746
1747 void
find_automorphisms(Stats & stats,void (* hook)(void * user_param,unsigned int n,const unsigned int * aut),void * user_param)1748 AbstractGraph::find_automorphisms(Stats& stats,
1749 void (*hook)(void *user_param,
1750 unsigned int n,
1751 const unsigned int *aut),
1752 void *user_param)
1753 {
1754 report_hook = hook;
1755 report_user_param = user_param;
1756
1757 search(false, stats);
1758
1759 if(first_path_labeling)
1760 {
1761 free(first_path_labeling);
1762 first_path_labeling = 0;
1763 }
1764 if(best_path_labeling)
1765 {
1766 free(best_path_labeling);
1767 best_path_labeling = 0;
1768 }
1769 }
1770
1771
1772 const unsigned int *
canonical_form(Stats & stats,void (* hook)(void * user_param,unsigned int n,const unsigned int * aut),void * user_param)1773 AbstractGraph::canonical_form(Stats& stats,
1774 void (*hook)(void *user_param,
1775 unsigned int n,
1776 const unsigned int *aut),
1777 void *user_param)
1778 {
1779
1780 report_hook = hook;
1781 report_user_param = user_param;
1782
1783 search(true, stats);
1784
1785 return best_path_labeling;
1786 }
1787
1788
1789
1790
1791 /*-------------------------------------------------------------------------
1792 *
1793 * Routines for directed graphs
1794 *
1795 *-------------------------------------------------------------------------*/
1796
Vertex()1797 Digraph::Vertex::Vertex()
1798 {
1799 color = 0;
1800 }
1801
1802
~Vertex()1803 Digraph::Vertex::~Vertex()
1804 {
1805 ;
1806 }
1807
1808
1809 void
add_edge_to(const unsigned int other_vertex)1810 Digraph::Vertex::add_edge_to(const unsigned int other_vertex)
1811 {
1812 edges_out.push_back(other_vertex);
1813 }
1814
1815
1816 void
add_edge_from(const unsigned int other_vertex)1817 Digraph::Vertex::add_edge_from(const unsigned int other_vertex)
1818 {
1819 edges_in.push_back(other_vertex);
1820 }
1821
1822
1823 void
remove_duplicate_edges(std::vector<bool> & tmp)1824 Digraph::Vertex::remove_duplicate_edges(std::vector<bool>& tmp)
1825 {
1826 #if defined(BLISS_CONSISTENCY_CHECKS)
1827 /* Pre-conditions */
1828 for(unsigned int i = 0; i < tmp.size(); i++) assert(tmp[i] == false);
1829 #endif
1830 for(std::vector<unsigned int>::iterator iter = edges_out.begin();
1831 iter != edges_out.end(); )
1832 {
1833 const unsigned int dest_vertex = *iter;
1834 if(tmp[dest_vertex] == true)
1835 {
1836 /* A duplicate edge found! */
1837 iter = edges_out.erase(iter);
1838 }
1839 else
1840 {
1841 /* Not seen earlier, mark as seen */
1842 tmp[dest_vertex] = true;
1843 iter++;
1844 }
1845 }
1846
1847 /* Clear tmp */
1848 for(std::vector<unsigned int>::iterator iter = edges_out.begin();
1849 iter != edges_out.end();
1850 iter++)
1851 {
1852 tmp[*iter] = false;
1853 }
1854
1855 for(std::vector<unsigned int>::iterator iter = edges_in.begin();
1856 iter != edges_in.end(); )
1857 {
1858 const unsigned int dest_vertex = *iter;
1859 if(tmp[dest_vertex] == true)
1860 {
1861 /* A duplicate edge found! */
1862 iter = edges_in.erase(iter);
1863 }
1864 else
1865 {
1866 /* Not seen earlier, mark as seen */
1867 tmp[dest_vertex] = true;
1868 iter++;
1869 }
1870 }
1871
1872 /* Clear tmp */
1873 for(std::vector<unsigned int>::iterator iter = edges_in.begin();
1874 iter != edges_in.end();
1875 iter++)
1876 {
1877 tmp[*iter] = false;
1878 }
1879 #if defined(BLISS_CONSISTENCY_CHECKS)
1880 /* Post-conditions */
1881 for(unsigned int i = 0; i < tmp.size(); i++) assert(tmp[i] == false);
1882 #endif
1883 }
1884
1885
1886 /**
1887 * Sort the edges entering and leaving the vertex according to
1888 * the vertex number of the other edge end.
1889 * Time complexity: O(e log(e)), where e is the number of edges
1890 * entering/leaving the vertex.
1891 */
1892 void
sort_edges()1893 Digraph::Vertex::sort_edges()
1894 {
1895 std::sort(edges_in.begin(), edges_in.end());
1896 std::sort(edges_out.begin(), edges_out.end());
1897 }
1898
1899
1900
1901
1902
1903 /*-------------------------------------------------------------------------
1904 *
1905 * Constructor and destructor for directed graphs
1906 *
1907 *-------------------------------------------------------------------------*/
1908
1909
Digraph(const unsigned int nof_vertices)1910 Digraph::Digraph(const unsigned int nof_vertices)
1911 {
1912 vertices.resize(nof_vertices);
1913 sh = shs_flm;
1914 }
1915
1916
~Digraph()1917 Digraph::~Digraph()
1918 {
1919 ;
1920 }
1921
1922
1923 unsigned int
add_vertex(const unsigned int color)1924 Digraph::add_vertex(const unsigned int color)
1925 {
1926 const unsigned int new_vertex_num = vertices.size();
1927 vertices.resize(new_vertex_num + 1);
1928 vertices.back().color = color;
1929 return new_vertex_num;
1930 }
1931
1932
1933 void
add_edge(const unsigned int vertex1,const unsigned int vertex2)1934 Digraph::add_edge(const unsigned int vertex1, const unsigned int vertex2)
1935 {
1936 assert(vertex1 < get_nof_vertices());
1937 assert(vertex2 < get_nof_vertices());
1938 vertices[vertex1].add_edge_to(vertex2);
1939 vertices[vertex2].add_edge_from(vertex1);
1940 }
1941
1942
1943 void
change_color(const unsigned int vertex,const unsigned int new_color)1944 Digraph::change_color(const unsigned int vertex, const unsigned int new_color)
1945 {
1946 assert(vertex < get_nof_vertices());
1947 vertices[vertex].color = new_color;
1948 }
1949
1950
1951 void
sort_edges()1952 Digraph::sort_edges()
1953 {
1954 for(unsigned int i = 0; i < get_nof_vertices(); i++)
1955 vertices[i].sort_edges();
1956 }
1957
1958
1959 int
cmp(Digraph & other)1960 Digraph::cmp(Digraph& other)
1961 {
1962 /* Compare the numbers of vertices */
1963 if(get_nof_vertices() < other.get_nof_vertices())
1964 return -1;
1965 if(get_nof_vertices() > other.get_nof_vertices())
1966 return 1;
1967 /* Compare vertex colors */
1968 for(unsigned int i = 0; i < get_nof_vertices(); i++)
1969 {
1970 if(vertices[i].color < other.vertices[i].color)
1971 return -1;
1972 if(vertices[i].color > other.vertices[i].color)
1973 return 1;
1974 }
1975 /* Compare vertex degrees */
1976 remove_duplicate_edges();
1977 other.remove_duplicate_edges();
1978 for(unsigned int i = 0; i < get_nof_vertices(); i++)
1979 {
1980 if(vertices[i].nof_edges_in() < other.vertices[i].nof_edges_in())
1981 return -1;
1982 if(vertices[i].nof_edges_in() > other.vertices[i].nof_edges_in())
1983 return 1;
1984 if(vertices[i].nof_edges_out() < other.vertices[i].nof_edges_out())
1985 return -1;
1986 if(vertices[i].nof_edges_out() > other.vertices[i].nof_edges_out())
1987 return 1;
1988 }
1989 /* Compare edges */
1990 for(unsigned int i = 0; i < get_nof_vertices(); i++)
1991 {
1992 Vertex& v1 = vertices[i];
1993 Vertex& v2 = other.vertices[i];
1994 v1.sort_edges();
1995 v2.sort_edges();
1996 std::vector<unsigned int>::const_iterator ei1 = v1.edges_in.begin();
1997 std::vector<unsigned int>::const_iterator ei2 = v2.edges_in.begin();
1998 while(ei1 != v1.edges_in.end())
1999 {
2000 if(*ei1 < *ei2)
2001 return -1;
2002 if(*ei1 > *ei2)
2003 return 1;
2004 ei1++;
2005 ei2++;
2006 }
2007 ei1 = v1.edges_out.begin();
2008 ei2 = v2.edges_out.begin();
2009 while(ei1 != v1.edges_out.end())
2010 {
2011 if(*ei1 < *ei2)
2012 return -1;
2013 if(*ei1 > *ei2)
2014 return 1;
2015 ei1++;
2016 ei2++;
2017 }
2018 }
2019 return 0;
2020 }
2021
2022
2023
2024
2025 Digraph*
permute(const std::vector<unsigned int> & perm) const2026 Digraph::permute(const std::vector<unsigned int>& perm) const
2027 {
2028 Digraph* const g = new Digraph(get_nof_vertices());
2029 for(unsigned int i = 0; i < get_nof_vertices(); i++)
2030 {
2031 const Vertex& v = vertices[i];
2032 g->change_color(perm[i], v.color);
2033 for(std::vector<unsigned int>::const_iterator ei = v.edges_out.begin();
2034 ei != v.edges_out.end();
2035 ei++)
2036 {
2037 g->add_edge(perm[i], perm[*ei]);
2038 }
2039 }
2040 g->sort_edges();
2041 return g;
2042 }
2043
2044
2045 Digraph*
permute(const unsigned int * const perm) const2046 Digraph::permute(const unsigned int* const perm) const
2047 {
2048 Digraph* const g = new Digraph(get_nof_vertices());
2049 for(unsigned int i = 0; i < get_nof_vertices(); i++)
2050 {
2051 const Vertex &v = vertices[i];
2052 g->change_color(perm[i], v.color);
2053 for(std::vector<unsigned int>::const_iterator ei = v.edges_out.begin();
2054 ei != v.edges_out.end();
2055 ei++)
2056 {
2057 g->add_edge(perm[i], perm[*ei]);
2058 }
2059 }
2060 g->sort_edges();
2061 return g;
2062 }
2063
2064
2065
2066
2067
2068 /*-------------------------------------------------------------------------
2069 *
2070 * Print graph in graphviz format
2071 *
2072 *-------------------------------------------------------------------------*/
2073
2074
2075 #if 0
2076 void
2077 Digraph::write_dot(const char* const filename)
2078 {
2079 FILE* const fp = fopen(filename, "w");
2080 if(fp)
2081 {
2082 write_dot(fp);
2083 fclose(fp);
2084 }
2085 }
2086
2087
2088 void
2089 Digraph::write_dot(FILE* const fp)
2090 {
2091 remove_duplicate_edges();
2092
2093 fprintf(fp, "digraph g {\n");
2094
2095 unsigned int vnum = 0;
2096 for(std::vector<Vertex>::const_iterator vi = vertices.begin();
2097 vi != vertices.end();
2098 vi++, vnum++)
2099 {
2100 const Vertex& v = *vi;
2101 fprintf(fp, "v%u [label=\"%u:%u\"];\n", vnum, vnum, v.color);
2102 for(std::vector<unsigned int>::const_iterator ei = v.edges_out.begin();
2103 ei != v.edges_out.end();
2104 ei++)
2105 {
2106 fprintf(fp, "v%u -> v%u\n", vnum, *ei);
2107 }
2108 }
2109
2110 fprintf(fp, "}\n");
2111 }
2112 #endif
2113
2114
2115 void
remove_duplicate_edges()2116 Digraph::remove_duplicate_edges()
2117 {
2118 std::vector<bool> tmp(get_nof_vertices(), false);
2119
2120 for(std::vector<Vertex>::iterator vi = vertices.begin();
2121 vi != vertices.end();
2122 vi++)
2123 {
2124 #if defined(BLISS_EXPENSIVE_CONSISTENCY_CHECKS)
2125 for(unsigned int i = 0; i < tmp.size(); i++) assert(tmp[i] == false);
2126 #endif
2127 (*vi).remove_duplicate_edges(tmp);
2128 }
2129 }
2130
2131
2132
2133
2134
2135 /*-------------------------------------------------------------------------
2136 *
2137 * Get a hash value for the graph.
2138 *
2139 *-------------------------------------------------------------------------*/
2140
2141 unsigned int
get_hash()2142 Digraph::get_hash()
2143 {
2144 remove_duplicate_edges();
2145 sort_edges();
2146
2147 UintSeqHash h;
2148
2149 h.update(get_nof_vertices());
2150
2151 /* Hash the color of each vertex */
2152 for(unsigned int i = 0; i < get_nof_vertices(); i++)
2153 {
2154 h.update(vertices[i].color);
2155 }
2156
2157 /* Hash the edges */
2158 for(unsigned int i = 0; i < get_nof_vertices(); i++)
2159 {
2160 Vertex &v = vertices[i];
2161 for(std::vector<unsigned int>::const_iterator ei = v.edges_out.begin();
2162 ei != v.edges_out.end();
2163 ei++)
2164 {
2165 h.update(i);
2166 h.update(*ei);
2167 }
2168 }
2169
2170 return h.get_value();
2171 }
2172
2173
2174
2175 /*-------------------------------------------------------------------------
2176 *
2177 * Read directed graph in the DIMACS format.
2178 * Returns 0 if an error occurred.
2179 *
2180 *-------------------------------------------------------------------------*/
2181
2182 #if 0
2183 Digraph*
2184 Digraph::read_dimacs(FILE* const fp, FILE* const errstr)
2185 {
2186 Digraph* g = 0;
2187 unsigned int nof_vertices;
2188 unsigned int nof_edges;
2189 unsigned int line_num = 1;
2190
2191 const bool verbose = false;
2192 FILE* const verbstr = stdout;
2193
2194 /* Read comments and the problem definition line */
2195 while(1)
2196 {
2197 int c = getc(fp);
2198 if(c == 'c')
2199 {
2200 /* A comment, ignore the rest of the line */
2201 while((c = getc(fp)) != '\n')
2202 {
2203 if(c == EOF) {
2204 if(errstr)
2205 fprintf(errstr, "error in line %u: not in DIMACS format\n",
2206 line_num);
2207 goto error_exit;
2208 }
2209 }
2210 line_num++;
2211 continue;
2212 }
2213 if(c == 'p')
2214 {
2215 /* The problem definition line */
2216 if(fscanf(fp, " edge %u %u\n", &nof_vertices, &nof_edges) != 2)
2217 {
2218 if(errstr)
2219 fprintf(errstr, "error in line %u: not in DIMACS format\n",
2220 line_num);
2221 goto error_exit;
2222 }
2223 line_num++;
2224 break;
2225 }
2226 if(errstr)
2227 fprintf(errstr, "error in line %u: not in DIMACS format\n", line_num);
2228 goto error_exit;
2229 }
2230
2231 if(nof_vertices <= 0)
2232 {
2233 if(errstr)
2234 fprintf(errstr, "error: no vertices\n");
2235 goto error_exit;
2236 }
2237 if(verbose)
2238 {
2239 fprintf(verbstr, "Instance has %d vertices and %d edges\n",
2240 nof_vertices, nof_edges);
2241 fflush(verbstr);
2242 }
2243
2244 g = new Digraph(nof_vertices);
2245
2246 //
2247 // Read vertex colors
2248 //
2249 if(verbose)
2250 {
2251 fprintf(verbstr, "Reading vertex colors...\n");
2252 fflush(verbstr);
2253 }
2254 while(1)
2255 {
2256 int c = getc(fp);
2257 if(c != 'n')
2258 {
2259 ungetc(c, fp);
2260 break;
2261 }
2262 ungetc(c, fp);
2263 unsigned int vertex;
2264 unsigned int color;
2265 if(fscanf(fp, "n %u %u\n", &vertex, &color) != 2)
2266 {
2267 if(errstr)
2268 fprintf(errstr, "error in line %u: not in DIMACS format\n",
2269 line_num);
2270 goto error_exit;
2271 }
2272 if(!((vertex >= 1) && (vertex <= nof_vertices)))
2273 {
2274 if(errstr)
2275 fprintf(errstr,
2276 "error in line %u: vertex %u not in range [1,...%u]\n",
2277 line_num, vertex, nof_vertices);
2278 goto error_exit;
2279 }
2280 line_num++;
2281 g->change_color(vertex - 1, color);
2282 }
2283 if(verbose)
2284 {
2285 fprintf(verbstr, "Done\n");
2286 fflush(verbstr);
2287 }
2288
2289 //
2290 // Read edges
2291 //
2292 if(verbose)
2293 {
2294 fprintf(verbstr, "Reading edges...\n");
2295 fflush(verbstr);
2296 }
2297 for(unsigned i = 0; i < nof_edges; i++)
2298 {
2299 unsigned int from, to;
2300 if(fscanf(fp, "e %u %u\n", &from, &to) != 2)
2301 {
2302 if(errstr)
2303 fprintf(errstr, "error in line %u: not in DIMACS format\n",
2304 line_num);
2305 goto error_exit;
2306 }
2307 if(not((1 <= from) and (from <= nof_vertices)))
2308 {
2309 if(errstr)
2310 fprintf(errstr,
2311 "error in line %u: vertex %u not in range [1,...%u]\n",
2312 line_num, from, nof_vertices);
2313 goto error_exit;
2314 }
2315 if(not((1 <= to) and (to <= nof_vertices)))
2316 {
2317 if(errstr)
2318 fprintf(errstr,
2319 "error in line %u: vertex %u not in range [1,...%u]\n",
2320 line_num, to, nof_vertices);
2321 goto error_exit;
2322 }
2323 line_num++;
2324 g->add_edge(from-1, to-1);
2325 }
2326 if(verbose)
2327 {
2328 fprintf(verbstr, "Done\n");
2329 fflush(verbstr);
2330 }
2331
2332 return g;
2333
2334 error_exit:
2335 if(g)
2336 delete g;
2337 return 0;
2338 }
2339 #endif
2340
2341
2342
2343 #if 0
2344 void
2345 Digraph::write_dimacs(FILE* const fp)
2346 {
2347 remove_duplicate_edges();
2348 sort_edges();
2349
2350 /* First count the total number of edges */
2351 unsigned int nof_edges = 0;
2352 for(unsigned int i = 0; i < get_nof_vertices(); i++)
2353 {
2354 nof_edges += vertices[i].edges_out.size();
2355 }
2356
2357 /* Output the "header" line */
2358 fprintf(fp, "p edge %u %u\n", get_nof_vertices(), nof_edges);
2359
2360 /* Print the color of each vertex */
2361 for(unsigned int i = 0; i < get_nof_vertices(); i++)
2362 {
2363 Vertex& v = vertices[i];
2364 fprintf(fp, "n %u %u\n", i+1, v.color);
2365 /*
2366 if(v.color != 0)
2367 {
2368 fprintf(fp, "n %u %u\n", i+1, v.color);
2369 }
2370 */
2371 }
2372
2373 /* Print the edges */
2374 for(unsigned int i = 0; i < get_nof_vertices(); i++)
2375 {
2376 Vertex& v = vertices[i];
2377 for(std::vector<unsigned int>::const_iterator ei = v.edges_out.begin();
2378 ei != v.edges_out.end();
2379 ei++)
2380 {
2381 fprintf(fp, "e %u %u\n", i+1, (*ei)+1);
2382 }
2383 }
2384 }
2385 #endif
2386
2387
2388
2389
2390
2391
2392
2393 /*-------------------------------------------------------------------------
2394 *
2395 * Partition independent invariants
2396 *
2397 *-------------------------------------------------------------------------*/
2398
2399 unsigned int
vertex_color_invariant(const Digraph * const g,const unsigned int vnum)2400 Digraph::vertex_color_invariant(const Digraph* const g, const unsigned int vnum)
2401 {
2402 return g->vertices[vnum].color;
2403 }
2404
2405 unsigned int
indegree_invariant(const Digraph * const g,const unsigned int vnum)2406 Digraph::indegree_invariant(const Digraph* const g, const unsigned int vnum)
2407 {
2408 return g->vertices[vnum].nof_edges_in();
2409 }
2410
2411 unsigned int
outdegree_invariant(const Digraph * const g,const unsigned int vnum)2412 Digraph::outdegree_invariant(const Digraph* const g, const unsigned int vnum)
2413 {
2414 return g->vertices[vnum].nof_edges_out();
2415 }
2416
2417 unsigned int
selfloop_invariant(const Digraph * const g,const unsigned int vnum)2418 Digraph::selfloop_invariant(const Digraph* const g, const unsigned int vnum)
2419 {
2420 /* Quite inefficient but luckily not in the critical path */
2421 const Vertex& v = g->vertices[vnum];
2422 for(std::vector<unsigned int>::const_iterator ei = v.edges_out.begin();
2423 ei != v.edges_out.end();
2424 ei++)
2425 {
2426 if(*ei == vnum)
2427 return 1;
2428 }
2429 return 0;
2430 }
2431
2432
2433
2434
2435
2436 /*-------------------------------------------------------------------------
2437 *
2438 * Refine the partition p according to a partition independent invariant
2439 *
2440 *-------------------------------------------------------------------------*/
2441
2442 bool
refine_according_to_invariant(unsigned int (* inv)(const Digraph * const g,const unsigned int v))2443 Digraph::refine_according_to_invariant(unsigned int (*inv)(const Digraph* const g,
2444 const unsigned int v))
2445 {
2446 bool refined = false;
2447
2448 for(Partition::Cell* cell = p.first_nonsingleton_cell; cell; )
2449 {
2450
2451 Partition::Cell* const next_cell = cell->next_nonsingleton;
2452 const unsigned int* ep = p.elements + cell->first;
2453 for(unsigned int i = cell->length; i > 0; i--, ep++)
2454 {
2455 unsigned int ival = inv(this, *ep);
2456 p.invariant_values[*ep] = ival;
2457 if(ival > cell->max_ival) {
2458 cell->max_ival = ival;
2459 cell->max_ival_count = 1;
2460 }
2461 else if(ival == cell->max_ival) {
2462 cell->max_ival_count++;
2463 }
2464 }
2465 Partition::Cell* const last_new_cell = p.zplit_cell(cell, true);
2466 refined |= (last_new_cell != cell);
2467 cell = next_cell;
2468 }
2469
2470 return refined;
2471 }
2472
2473
2474
2475
2476
2477 /*-------------------------------------------------------------------------
2478 *
2479 * Split the neighbourhood of a cell according to the equitable invariant
2480 *
2481 *-------------------------------------------------------------------------*/
2482
2483 bool
split_neighbourhood_of_cell(Partition::Cell * const cell)2484 Digraph::split_neighbourhood_of_cell(Partition::Cell* const cell)
2485 {
2486
2487
2488 const bool was_equal_to_first = refine_equal_to_first;
2489
2490 if(compute_eqref_hash)
2491 {
2492 eqref_hash.update(cell->first);
2493 eqref_hash.update(cell->length);
2494 }
2495
2496 const unsigned int* ep = p.elements + cell->first;
2497 for(unsigned int i = cell->length; i > 0; i--)
2498 {
2499 const Vertex& v = vertices[*ep++];
2500
2501 std::vector<unsigned int>::const_iterator ei = v.edges_out.begin();
2502 for(unsigned int j = v.nof_edges_out(); j != 0; j--)
2503 {
2504 const unsigned int dest_vertex = *ei++;
2505 Partition::Cell* const neighbour_cell = p.get_cell(dest_vertex);
2506 if(neighbour_cell->is_unit())
2507 continue;
2508 const unsigned int ival = ++p.invariant_values[dest_vertex];
2509 if(ival > neighbour_cell->max_ival) {
2510 neighbour_cell->max_ival = ival;
2511 neighbour_cell->max_ival_count = 1;
2512 if(ival == 1)
2513 neighbour_heap.insert(neighbour_cell->first);
2514 }
2515 else if(ival == neighbour_cell->max_ival) {
2516 neighbour_cell->max_ival_count++;
2517 }
2518 }
2519 }
2520
2521 while(!neighbour_heap.is_empty())
2522 {
2523 const unsigned int start = neighbour_heap.remove();
2524 Partition::Cell* const neighbour_cell = p.get_cell(p.elements[start]);
2525
2526 if(compute_eqref_hash)
2527 {
2528 eqref_hash.update(neighbour_cell->first);
2529 eqref_hash.update(neighbour_cell->length);
2530 eqref_hash.update(neighbour_cell->max_ival);
2531 eqref_hash.update(neighbour_cell->max_ival_count);
2532 }
2533
2534
2535 Partition::Cell* const last_new_cell = p.zplit_cell(neighbour_cell, true);
2536
2537 /* Update certificate and hash if needed */
2538 const Partition::Cell* c = neighbour_cell;
2539 while(1)
2540 {
2541 if(in_search)
2542 {
2543 /* Build certificate */
2544 cert_add_redundant(CERT_SPLIT, c->first, c->length);
2545 /* No need to continue? */
2546 if(refine_compare_certificate and
2547 (refine_equal_to_first == false) and
2548 (refine_cmp_to_best < 0))
2549 goto worse_exit;
2550 }
2551 if(compute_eqref_hash)
2552 {
2553 eqref_hash.update(c->first);
2554 eqref_hash.update(c->length);
2555 }
2556 if(c == last_new_cell)
2557 break;
2558 c = c->next;
2559 }
2560 }
2561
2562 if(cell->is_in_splitting_queue())
2563 {
2564 return false;
2565 }
2566
2567
2568 ep = p.elements + cell->first;
2569 for(unsigned int i = cell->length; i > 0; i--)
2570 {
2571 const Vertex& v = vertices[*ep++];
2572
2573 std::vector<unsigned int>::const_iterator ei = v.edges_in.begin();
2574 for(unsigned int j = v.nof_edges_in(); j > 0; j--)
2575 {
2576 const unsigned int dest_vertex = *ei++;
2577 Partition::Cell* const neighbour_cell = p.get_cell(dest_vertex);
2578 if(neighbour_cell->is_unit())
2579 continue;
2580 const unsigned int ival = ++p.invariant_values[dest_vertex];
2581 if(ival > neighbour_cell->max_ival)
2582 {
2583 neighbour_cell->max_ival = ival;
2584 neighbour_cell->max_ival_count = 1;
2585 if(ival == 1)
2586 neighbour_heap.insert(neighbour_cell->first);
2587 }
2588 else if(ival == neighbour_cell->max_ival) {
2589 neighbour_cell->max_ival_count++;
2590 }
2591 }
2592 }
2593
2594 while(!neighbour_heap.is_empty())
2595 {
2596 const unsigned int start = neighbour_heap.remove();
2597 Partition::Cell* const neighbour_cell = p.get_cell(p.elements[start]);
2598
2599 if(compute_eqref_hash)
2600 {
2601 eqref_hash.update(neighbour_cell->first);
2602 eqref_hash.update(neighbour_cell->length);
2603 eqref_hash.update(neighbour_cell->max_ival);
2604 eqref_hash.update(neighbour_cell->max_ival_count);
2605 }
2606
2607 Partition::Cell* const last_new_cell = p.zplit_cell(neighbour_cell, true);
2608
2609 /* Update certificate and hash if needed */
2610 const Partition::Cell* c = neighbour_cell;
2611 while(1)
2612 {
2613 if(in_search)
2614 {
2615 /* Build certificate */
2616 cert_add_redundant(CERT_SPLIT, c->first, c->length);
2617 /* No need to continue? */
2618 if(refine_compare_certificate and
2619 (refine_equal_to_first == false) and
2620 (refine_cmp_to_best < 0))
2621 goto worse_exit;
2622 }
2623 if(compute_eqref_hash)
2624 {
2625 eqref_hash.update(c->first);
2626 eqref_hash.update(c->length);
2627 }
2628 if(c == last_new_cell)
2629 break;
2630 c = c->next;
2631 }
2632 }
2633
2634
2635 if(refine_compare_certificate and
2636 (refine_equal_to_first == false) and
2637 (refine_cmp_to_best < 0))
2638 return true;
2639
2640 return false;
2641
2642 worse_exit:
2643 /* Clear neighbour heap */
2644 UintSeqHash rest;
2645 while(!neighbour_heap.is_empty())
2646 {
2647 const unsigned int start = neighbour_heap.remove();
2648 Partition::Cell* const neighbour_cell = p.get_cell(p.elements[start]);
2649 if(opt_use_failure_recording and was_equal_to_first)
2650 {
2651 rest.update(neighbour_cell->first);
2652 rest.update(neighbour_cell->length);
2653 rest.update(neighbour_cell->max_ival);
2654 rest.update(neighbour_cell->max_ival_count);
2655 }
2656 neighbour_cell->max_ival = 0;
2657 neighbour_cell->max_ival_count = 0;
2658 p.clear_ivs(neighbour_cell);
2659 }
2660 if(opt_use_failure_recording and was_equal_to_first)
2661 {
2662 for(unsigned int i = p.splitting_queue.size(); i > 0; i--)
2663 {
2664 Partition::Cell* const cell = p.splitting_queue.pop_front();
2665 rest.update(cell->first);
2666 rest.update(cell->length);
2667 p.splitting_queue.push_back(cell);
2668 }
2669 rest.update(failure_recording_fp_deviation);
2670 failure_recording_fp_deviation = rest.get_value();
2671 }
2672
2673 return true;
2674 }
2675
2676
2677 bool
split_neighbourhood_of_unit_cell(Partition::Cell * const unit_cell)2678 Digraph::split_neighbourhood_of_unit_cell(Partition::Cell* const unit_cell)
2679 {
2680
2681
2682 const bool was_equal_to_first = refine_equal_to_first;
2683
2684 if(compute_eqref_hash)
2685 {
2686 eqref_hash.update(0x87654321);
2687 eqref_hash.update(unit_cell->first);
2688 eqref_hash.update(1);
2689 }
2690
2691 const Vertex& v = vertices[p.elements[unit_cell->first]];
2692
2693 /*
2694 * Phase 1
2695 * Refine neighbours according to the edges that leave the vertex v
2696 */
2697 std::vector<unsigned int>::const_iterator ei = v.edges_out.begin();
2698 for(unsigned int j = v.nof_edges_out(); j > 0; j--)
2699 {
2700 const unsigned int dest_vertex = *ei++;
2701 Partition::Cell* const neighbour_cell = p.get_cell(dest_vertex);
2702
2703 if(neighbour_cell->is_unit()) {
2704 if(in_search) {
2705 /* Remember neighbour in order to generate certificate */
2706 neighbour_heap.insert(neighbour_cell->first);
2707 }
2708 continue;
2709 }
2710 if(neighbour_cell->max_ival_count == 0)
2711 {
2712 neighbour_heap.insert(neighbour_cell->first);
2713 }
2714 neighbour_cell->max_ival_count++;
2715
2716 unsigned int* const swap_position =
2717 p.elements + neighbour_cell->first + neighbour_cell->length -
2718 neighbour_cell->max_ival_count;
2719 *p.in_pos[dest_vertex] = *swap_position;
2720 p.in_pos[*swap_position] = p.in_pos[dest_vertex];
2721 *swap_position = dest_vertex;
2722 p.in_pos[dest_vertex] = swap_position;
2723 }
2724
2725 while(!neighbour_heap.is_empty())
2726 {
2727 const unsigned int start = neighbour_heap.remove();
2728 Partition::Cell* neighbour_cell = p.get_cell(p.elements[start]);
2729
2730 #if defined(BLISS_CONSISTENCY_CHECKS)
2731 assert(neighbour_cell->first == start);
2732 if(neighbour_cell->is_unit()) {
2733 assert(neighbour_cell->max_ival_count == 0);
2734 } else {
2735 assert(neighbour_cell->max_ival_count > 0);
2736 assert(neighbour_cell->max_ival_count <= neighbour_cell->length);
2737 }
2738 #endif
2739
2740 if(compute_eqref_hash)
2741 {
2742 eqref_hash.update(neighbour_cell->first);
2743 eqref_hash.update(neighbour_cell->length);
2744 eqref_hash.update(neighbour_cell->max_ival_count);
2745 }
2746
2747 if(neighbour_cell->length > 1 and
2748 neighbour_cell->max_ival_count != neighbour_cell->length)
2749 {
2750
2751 Partition::Cell* const new_cell =
2752 p.aux_split_in_two(neighbour_cell,
2753 neighbour_cell->length -
2754 neighbour_cell->max_ival_count);
2755 unsigned int* ep = p.elements + new_cell->first;
2756 unsigned int* const lp = p.elements+new_cell->first+new_cell->length;
2757 while(ep < lp)
2758 {
2759 p.element_to_cell_map[*ep] = new_cell;
2760 ep++;
2761 }
2762 neighbour_cell->max_ival_count = 0;
2763
2764
2765 if(compute_eqref_hash)
2766 {
2767 /* Update hash */
2768 eqref_hash.update(neighbour_cell->first);
2769 eqref_hash.update(neighbour_cell->length);
2770 eqref_hash.update(0);
2771 eqref_hash.update(new_cell->first);
2772 eqref_hash.update(new_cell->length);
2773 eqref_hash.update(1);
2774 }
2775
2776 /* Add cells in splitting_queue */
2777 if(neighbour_cell->is_in_splitting_queue()) {
2778 /* Both cells must be included in splitting_queue in order
2779 to have refinement to equitable partition */
2780 p.splitting_queue_add(new_cell);
2781 } else {
2782 Partition::Cell *min_cell, *max_cell;
2783 if(neighbour_cell->length <= new_cell->length) {
2784 min_cell = neighbour_cell;
2785 max_cell = new_cell;
2786 } else {
2787 min_cell = new_cell;
2788 max_cell = neighbour_cell;
2789 }
2790 /* Put the smaller cell in splitting_queue */
2791 p.splitting_queue_add(min_cell);
2792 if(max_cell->is_unit()) {
2793 /* Put the "larger" cell also in splitting_queue */
2794 p.splitting_queue_add(max_cell);
2795 }
2796 }
2797 /* Update pointer for certificate generation */
2798 neighbour_cell = new_cell;
2799 }
2800 else
2801 {
2802 neighbour_cell->max_ival_count = 0;
2803 }
2804
2805 /*
2806 * Build certificate if required
2807 */
2808 if(in_search)
2809 {
2810 for(unsigned int i = neighbour_cell->first,
2811 j = neighbour_cell->length;
2812 j > 0;
2813 j--, i++)
2814 {
2815 /* Build certificate */
2816 cert_add(CERT_EDGE, unit_cell->first, i);
2817 /* No need to continue? */
2818 if(refine_compare_certificate and
2819 (refine_equal_to_first == false) and
2820 (refine_cmp_to_best < 0))
2821 goto worse_exit;
2822 }
2823 } /* if(in_search) */
2824 } /* while(!neighbour_heap.is_empty()) */
2825
2826 /*
2827 * Phase 2
2828 * Refine neighbours according to the edges that enter the vertex v
2829 */
2830 ei = v.edges_in.begin();
2831 for(unsigned int j = v.nof_edges_in(); j > 0; j--)
2832 {
2833 const unsigned int dest_vertex = *ei++;
2834 Partition::Cell* const neighbour_cell = p.get_cell(dest_vertex);
2835
2836 if(neighbour_cell->is_unit()) {
2837 if(in_search) {
2838 neighbour_heap.insert(neighbour_cell->first);
2839 }
2840 continue;
2841 }
2842 if(neighbour_cell->max_ival_count == 0)
2843 {
2844 neighbour_heap.insert(neighbour_cell->first);
2845 }
2846 neighbour_cell->max_ival_count++;
2847
2848 unsigned int* const swap_position =
2849 p.elements + neighbour_cell->first + neighbour_cell->length -
2850 neighbour_cell->max_ival_count;
2851 *p.in_pos[dest_vertex] = *swap_position;
2852 p.in_pos[*swap_position] = p.in_pos[dest_vertex];
2853 *swap_position = dest_vertex;
2854 p.in_pos[dest_vertex] = swap_position;
2855 }
2856
2857 while(!neighbour_heap.is_empty())
2858 {
2859 const unsigned int start = neighbour_heap.remove();
2860 Partition::Cell* neighbour_cell = p.get_cell(p.elements[start]);
2861
2862 #if defined(BLISS_CONSISTENCY_CHECKS)
2863 assert(neighbour_cell->first == start);
2864 if(neighbour_cell->is_unit()) {
2865 assert(neighbour_cell->max_ival_count == 0);
2866 } else {
2867 assert(neighbour_cell->max_ival_count > 0);
2868 assert(neighbour_cell->max_ival_count <= neighbour_cell->length);
2869 }
2870 #endif
2871
2872 if(compute_eqref_hash)
2873 {
2874 eqref_hash.update(neighbour_cell->first);
2875 eqref_hash.update(neighbour_cell->length);
2876 eqref_hash.update(neighbour_cell->max_ival_count);
2877 }
2878
2879 if(neighbour_cell->length > 1 and
2880 neighbour_cell->max_ival_count != neighbour_cell->length)
2881 {
2882 Partition::Cell* const new_cell =
2883 p.aux_split_in_two(neighbour_cell,
2884 neighbour_cell->length -
2885 neighbour_cell->max_ival_count);
2886 unsigned int* ep = p.elements + new_cell->first;
2887 unsigned int* const lp = p.elements+new_cell->first+new_cell->length;
2888 while(ep < lp) {
2889 p.element_to_cell_map[*ep] = new_cell;
2890 ep++;
2891 }
2892 neighbour_cell->max_ival_count = 0;
2893
2894
2895 if(compute_eqref_hash)
2896 {
2897 eqref_hash.update(neighbour_cell->first);
2898 eqref_hash.update(neighbour_cell->length);
2899 eqref_hash.update(0);
2900 eqref_hash.update(new_cell->first);
2901 eqref_hash.update(new_cell->length);
2902 eqref_hash.update(1);
2903 }
2904
2905 /* Add cells in splitting_queue */
2906 if(neighbour_cell->is_in_splitting_queue()) {
2907 /* Both cells must be included in splitting_queue in order
2908 to have refinement to equitable partition */
2909 p.splitting_queue_add(new_cell);
2910 } else {
2911 Partition::Cell *min_cell, *max_cell;
2912 if(neighbour_cell->length <= new_cell->length) {
2913 min_cell = neighbour_cell;
2914 max_cell = new_cell;
2915 } else {
2916 min_cell = new_cell;
2917 max_cell = neighbour_cell;
2918 }
2919 /* Put the smaller cell in splitting_queue */
2920 p.splitting_queue_add(min_cell);
2921 if(max_cell->is_unit()) {
2922 /* Put the "larger" cell also in splitting_queue */
2923 p.splitting_queue_add(max_cell);
2924 }
2925 }
2926 /* Update pointer for certificate generation */
2927 neighbour_cell = new_cell;
2928 }
2929 else
2930 {
2931 neighbour_cell->max_ival_count = 0;
2932 }
2933
2934 /*
2935 * Build certificate if required
2936 */
2937 if(in_search)
2938 {
2939 for(unsigned int i = neighbour_cell->first,
2940 j = neighbour_cell->length;
2941 j > 0;
2942 j--, i++)
2943 {
2944 /* Build certificate */
2945 cert_add(CERT_EDGE, i, unit_cell->first);
2946 /* No need to continue? */
2947 if(refine_compare_certificate and
2948 (refine_equal_to_first == false) and
2949 (refine_cmp_to_best < 0))
2950 goto worse_exit;
2951 }
2952 } /* if(in_search) */
2953 } /* while(!neighbour_heap.is_empty()) */
2954
2955 if(refine_compare_certificate and
2956 (refine_equal_to_first == false) and
2957 (refine_cmp_to_best < 0))
2958 return true;
2959
2960 return false;
2961
2962 worse_exit:
2963 /* Clear neighbour heap */
2964 UintSeqHash rest;
2965 while(!neighbour_heap.is_empty())
2966 {
2967 const unsigned int start = neighbour_heap.remove();
2968 Partition::Cell* const neighbour_cell = p.get_cell(p.elements[start]);
2969 if(opt_use_failure_recording and was_equal_to_first)
2970 {
2971 rest.update(neighbour_cell->first);
2972 rest.update(neighbour_cell->length);
2973 rest.update(neighbour_cell->max_ival_count);
2974 }
2975 neighbour_cell->max_ival_count = 0;
2976 }
2977 if(opt_use_failure_recording and was_equal_to_first)
2978 {
2979 rest.update(failure_recording_fp_deviation);
2980 failure_recording_fp_deviation = rest.get_value();
2981 }
2982 return true;
2983 }
2984
2985
2986
2987
2988
2989 /*-------------------------------------------------------------------------
2990 *
2991 * Check whether the current partition p is equitable.
2992 * Performance: very slow, use only for debugging purposes.
2993 *
2994 *-------------------------------------------------------------------------*/
2995
2996 bool
is_equitable() const2997 Digraph::is_equitable() const
2998 {
2999 const unsigned int N = get_nof_vertices();
3000 if(N == 0)
3001 return true;
3002
3003 std::vector<unsigned int> first_count = std::vector<unsigned int>(N, 0);
3004 std::vector<unsigned int> other_count = std::vector<unsigned int>(N, 0);
3005
3006 /*
3007 * Check equitabledness w.r.t. outgoing edges
3008 */
3009 for(Partition::Cell* cell = p.first_cell; cell; cell = cell->next)
3010 {
3011 if(cell->is_unit())
3012 continue;
3013
3014 unsigned int* ep = p.elements + cell->first;
3015 const Vertex& first_vertex = vertices[*ep++];
3016
3017 /* Count outgoing edges of the first vertex for cells */
3018 for(std::vector<unsigned int>::const_iterator ei =
3019 first_vertex.edges_out.begin();
3020 ei != first_vertex.edges_out.end();
3021 ei++)
3022 {
3023 first_count[p.get_cell(*ei)->first]++;
3024 }
3025
3026 /* Count and compare outgoing edges of the other vertices */
3027 for(unsigned int i = cell->length; i > 1; i--)
3028 {
3029 const Vertex &vertex = vertices[*ep++];
3030 for(std::vector<unsigned int>::const_iterator ei =
3031 vertex.edges_out.begin();
3032 ei != vertex.edges_out.end();
3033 ei++)
3034 {
3035 other_count[p.get_cell(*ei)->first]++;
3036 }
3037 for(Partition::Cell *cell2 = p.first_cell;
3038 cell2;
3039 cell2 = cell2->next)
3040 {
3041 if(first_count[cell2->first] != other_count[cell2->first])
3042 {
3043 /* Not equitable */
3044 return false;
3045 }
3046 other_count[cell2->first] = 0;
3047 }
3048 }
3049 /* Reset first_count */
3050 for(unsigned int i = 0; i < N; i++)
3051 first_count[i] = 0;
3052 }
3053
3054
3055 /*
3056 * Check equitabledness w.r.t. incoming edges
3057 */
3058 for(Partition::Cell* cell = p.first_cell; cell; cell = cell->next)
3059 {
3060 if(cell->is_unit())
3061 continue;
3062
3063 unsigned int* ep = p.elements + cell->first;
3064 const Vertex& first_vertex = vertices[*ep++];
3065
3066 /* Count incoming edges of the first vertex for cells */
3067 for(std::vector<unsigned int>::const_iterator ei =
3068 first_vertex.edges_in.begin();
3069 ei != first_vertex.edges_in.end();
3070 ei++)
3071 {
3072 first_count[p.get_cell(*ei)->first]++;
3073 }
3074
3075 /* Count and compare incoming edges of the other vertices */
3076 for(unsigned int i = cell->length; i > 1; i--)
3077 {
3078 const Vertex &vertex = vertices[*ep++];
3079 for(std::vector<unsigned int>::const_iterator ei =
3080 vertex.edges_in.begin();
3081 ei != vertex.edges_in.end();
3082 ei++)
3083 {
3084 other_count[p.get_cell(*ei)->first]++;
3085 }
3086 for(Partition::Cell *cell2 = p.first_cell;
3087 cell2;
3088 cell2 = cell2->next)
3089 {
3090 if(first_count[cell2->first] != other_count[cell2->first])
3091 {
3092 /* Not equitable */
3093 return false;
3094 }
3095 other_count[cell2->first] = 0;
3096 }
3097 }
3098 /* Reset first_count */
3099 for(unsigned int i = 0; i < N; i++)
3100 first_count[i] = 0;
3101 }
3102 return true;
3103 }
3104
3105
3106
3107
3108
3109 /*-------------------------------------------------------------------------
3110 *
3111 * Build the initial equitable partition
3112 *
3113 *-------------------------------------------------------------------------*/
3114
3115 void
make_initial_equitable_partition()3116 Digraph::make_initial_equitable_partition()
3117 {
3118 refine_according_to_invariant(&vertex_color_invariant);
3119 p.splitting_queue_clear();
3120 //p.print_signature(stderr); fprintf(stderr, "\n");
3121
3122 refine_according_to_invariant(&selfloop_invariant);
3123 p.splitting_queue_clear();
3124 //p.print_signature(stderr); fprintf(stderr, "\n");
3125
3126 refine_according_to_invariant(&outdegree_invariant);
3127 p.splitting_queue_clear();
3128 //p.print_signature(stderr); fprintf(stderr, "\n");
3129
3130 refine_according_to_invariant(&indegree_invariant);
3131 p.splitting_queue_clear();
3132 //p.print_signature(stderr); fprintf(stderr, "\n");
3133
3134 refine_to_equitable();
3135 //p.print_signature(stderr); fprintf(stderr, "\n");
3136 }
3137
3138
3139
3140
3141
3142 /*-------------------------------------------------------------------------
3143 *
3144 * Find the next cell to be splitted
3145 *
3146 *-------------------------------------------------------------------------*/
3147
3148 Partition::Cell*
find_next_cell_to_be_splitted(Partition::Cell * cell)3149 Digraph::find_next_cell_to_be_splitted(Partition::Cell* cell)
3150 {
3151 switch(sh) {
3152 case shs_f: return sh_first();
3153 case shs_fs: return sh_first_smallest();
3154 case shs_fl: return sh_first_largest();
3155 case shs_fm: return sh_first_max_neighbours();
3156 case shs_fsm: return sh_first_smallest_max_neighbours();
3157 case shs_flm: return sh_first_largest_max_neighbours();
3158 default:
3159 fatal_error("Internal error - unknown splitting heuristics");
3160 return 0;
3161 }
3162 }
3163
3164 /** \internal
3165 * A splitting heuristic.
3166 * Returns the first nonsingleton cell in the current partition.
3167 * The argument \a cell is ignored.
3168 */
3169 Partition::Cell*
sh_first()3170 Digraph::sh_first()
3171 {
3172 Partition::Cell* best_cell = 0;
3173 for(Partition::Cell* cell = p.first_nonsingleton_cell;
3174 cell;
3175 cell = cell->next_nonsingleton)
3176 {
3177 if(opt_use_comprec and p.cr_get_level(cell->first) != cr_level)
3178 continue;
3179 best_cell = cell;
3180 break;
3181 }
3182 return best_cell;
3183 }
3184
3185 /** \internal
3186 * A splitting heuristic.
3187 * Returns the first smallest nonsingleton cell in the current partition.
3188 * The argument \a cell is ignored.
3189 */
3190 Partition::Cell*
sh_first_smallest()3191 Digraph::sh_first_smallest()
3192 {
3193 Partition::Cell* best_cell = 0;
3194 unsigned int best_size = UINT_MAX;
3195 for(Partition::Cell* cell = p.first_nonsingleton_cell;
3196 cell;
3197 cell = cell->next_nonsingleton)
3198 {
3199 if(opt_use_comprec and p.cr_get_level(cell->first) != cr_level)
3200 continue;
3201 if(cell->length < best_size)
3202 {
3203 best_size = cell->length;
3204 best_cell = cell;
3205 }
3206 }
3207 return best_cell;
3208 }
3209
3210 /** \internal
3211 * A splitting heuristic.
3212 * Returns the first largest nonsingleton cell in the current partition.
3213 * The argument \a cell is ignored.
3214 */
3215 Partition::Cell*
sh_first_largest()3216 Digraph::sh_first_largest()
3217 {
3218 Partition::Cell* best_cell = 0;
3219 unsigned int best_size = 0;
3220 for(Partition::Cell* cell = p.first_nonsingleton_cell;
3221 cell;
3222 cell = cell->next_nonsingleton)
3223 {
3224 if(opt_use_comprec and p.cr_get_level(cell->first) != cr_level)
3225 continue;
3226 if(cell->length > best_size)
3227 {
3228 best_size = cell->length;
3229 best_cell = cell;
3230 }
3231 }
3232 return best_cell;
3233 }
3234
3235 /** \internal
3236 * A splitting heuristic.
3237 * Returns the first nonsingleton cell with max number of neighbouring
3238 * nonsingleton cells.
3239 * Assumes that the partition p is equitable.
3240 * Assumes that the max_ival fields of the cells are all 0.
3241 */
3242 Partition::Cell*
sh_first_max_neighbours()3243 Digraph::sh_first_max_neighbours()
3244 {
3245 Partition::Cell* best_cell = 0;
3246 int best_value = -1;
3247 KStack<Partition::Cell*> neighbour_cells_visited;
3248 neighbour_cells_visited.init(get_nof_vertices());
3249 for(Partition::Cell* cell = p.first_nonsingleton_cell;
3250 cell;
3251 cell = cell->next_nonsingleton)
3252 {
3253 if(opt_use_comprec and p.cr_get_level(cell->first) != cr_level)
3254 continue;
3255 int value = 0;
3256 const Vertex &v = vertices[p.elements[cell->first]];
3257 std::vector<unsigned int>::const_iterator ei;
3258 ei = v.edges_in.begin();
3259 for(unsigned int j = v.nof_edges_in(); j > 0; j--)
3260 {
3261 Partition::Cell * const neighbour_cell = p.get_cell(*ei++);
3262 if(neighbour_cell->is_unit())
3263 continue;
3264 neighbour_cell->max_ival++;
3265 if(neighbour_cell->max_ival == 1)
3266 neighbour_cells_visited.push(neighbour_cell);
3267 }
3268 while(!neighbour_cells_visited.is_empty())
3269 {
3270 Partition::Cell* const neighbour_cell = neighbour_cells_visited.pop();
3271 if(neighbour_cell->max_ival != neighbour_cell->length)
3272 value++;
3273 neighbour_cell->max_ival = 0;
3274 }
3275
3276 ei = v.edges_out.begin();
3277 for(unsigned int j = v.nof_edges_out(); j > 0; j--)
3278 {
3279 Partition::Cell * const neighbour_cell = p.get_cell(*ei++);
3280 if(neighbour_cell->is_unit())
3281 continue;
3282 neighbour_cell->max_ival++;
3283 if(neighbour_cell->max_ival == 1)
3284 neighbour_cells_visited.push(neighbour_cell);
3285 }
3286 while(!neighbour_cells_visited.is_empty())
3287 {
3288 Partition::Cell* const neighbour_cell = neighbour_cells_visited.pop();
3289 if(neighbour_cell->max_ival != neighbour_cell->length)
3290 value++;
3291 neighbour_cell->max_ival = 0;
3292 }
3293
3294 if(value > best_value)
3295 {
3296 best_value = value;
3297 best_cell = cell;
3298 }
3299 }
3300 return best_cell;
3301 }
3302
3303 /** \internal
3304 * A splitting heuristic.
3305 * Returns the first smallest nonsingleton cell with max number of neighbouring
3306 * nonsingleton cells.
3307 * Assumes that the partition p is equitable.
3308 * Assumes that the max_ival fields of the cells are all 0.
3309 */
3310 Partition::Cell*
sh_first_smallest_max_neighbours()3311 Digraph::sh_first_smallest_max_neighbours()
3312 {
3313 Partition::Cell* best_cell = 0;
3314 int best_value = -1;
3315 unsigned int best_size = UINT_MAX;
3316 KStack<Partition::Cell*> neighbour_cells_visited;
3317 neighbour_cells_visited.init(get_nof_vertices());
3318 for(Partition::Cell* cell = p.first_nonsingleton_cell;
3319 cell;
3320 cell = cell->next_nonsingleton)
3321 {
3322
3323 if(opt_use_comprec and p.cr_get_level(cell->first) != cr_level)
3324 continue;
3325
3326 int value = 0;
3327 const Vertex& v = vertices[p.elements[cell->first]];
3328 std::vector<unsigned int>::const_iterator ei;
3329
3330 ei = v.edges_in.begin();
3331 for(unsigned int j = v.nof_edges_in(); j > 0; j--)
3332 {
3333 Partition::Cell * const neighbour_cell = p.get_cell(*ei++);
3334 if(neighbour_cell->is_unit())
3335 continue;
3336 neighbour_cell->max_ival++;
3337 if(neighbour_cell->max_ival == 1)
3338 neighbour_cells_visited.push(neighbour_cell);
3339 }
3340 while(!neighbour_cells_visited.is_empty())
3341 {
3342 Partition::Cell * const neighbour_cell = neighbour_cells_visited.pop();
3343 if(neighbour_cell->max_ival != neighbour_cell->length)
3344 value++;
3345 neighbour_cell->max_ival = 0;
3346 }
3347
3348 ei = v.edges_out.begin();
3349 for(unsigned int j = v.nof_edges_out(); j > 0; j--)
3350 {
3351 Partition::Cell * const neighbour_cell = p.get_cell(*ei++);
3352 if(neighbour_cell->is_unit())
3353 continue;
3354 neighbour_cell->max_ival++;
3355 if(neighbour_cell->max_ival == 1)
3356 neighbour_cells_visited.push(neighbour_cell);
3357 }
3358 while(!neighbour_cells_visited.is_empty())
3359 {
3360 Partition::Cell * const neighbour_cell = neighbour_cells_visited.pop();
3361 if(neighbour_cell->max_ival != neighbour_cell->length)
3362 value++;
3363 neighbour_cell->max_ival = 0;
3364 }
3365
3366 if((value > best_value) or
3367 (value == best_value and cell->length < best_size))
3368 {
3369 best_value = value;
3370 best_size = cell->length;
3371 best_cell = cell;
3372 }
3373 }
3374 return best_cell;
3375 }
3376
3377 /** \internal
3378 * A splitting heuristic.
3379 * Returns the first largest nonsingleton cell with max number of neighbouring
3380 * nonsingleton cells.
3381 * Assumes that the partition p is equitable.
3382 * Assumes that the max_ival fields of the cells are all 0.
3383 */
3384 Partition::Cell*
sh_first_largest_max_neighbours()3385 Digraph::sh_first_largest_max_neighbours()
3386 {
3387 Partition::Cell* best_cell = 0;
3388 int best_value = -1;
3389 unsigned int best_size = 0;
3390 KStack<Partition::Cell*> neighbour_cells_visited;
3391 neighbour_cells_visited.init(get_nof_vertices());
3392 for(Partition::Cell* cell = p.first_nonsingleton_cell;
3393 cell;
3394 cell = cell->next_nonsingleton)
3395 {
3396
3397 if(opt_use_comprec and p.cr_get_level(cell->first) != cr_level)
3398 continue;
3399
3400 int value = 0;
3401 const Vertex &v = vertices[p.elements[cell->first]];
3402 std::vector<unsigned int>::const_iterator ei;
3403
3404 ei = v.edges_in.begin();
3405 for(unsigned int j = v.nof_edges_in(); j > 0; j--)
3406 {
3407 Partition::Cell* const neighbour_cell = p.get_cell(*ei++);
3408 if(neighbour_cell->is_unit())
3409 continue;
3410 neighbour_cell->max_ival++;
3411 if(neighbour_cell->max_ival == 1)
3412 neighbour_cells_visited.push(neighbour_cell);
3413 }
3414 while(!neighbour_cells_visited.is_empty())
3415 {
3416 Partition::Cell* const neighbour_cell = neighbour_cells_visited.pop();
3417 if(neighbour_cell->max_ival != neighbour_cell->length)
3418 value++;
3419 neighbour_cell->max_ival = 0;
3420 }
3421
3422 ei = v.edges_out.begin();
3423 for(unsigned int j = v.nof_edges_out(); j > 0; j--)
3424 {
3425 Partition::Cell* const neighbour_cell = p.get_cell(*ei++);
3426 if(neighbour_cell->is_unit())
3427 continue;
3428 neighbour_cell->max_ival++;
3429 if(neighbour_cell->max_ival == 1)
3430 neighbour_cells_visited.push(neighbour_cell);
3431 }
3432 while(!neighbour_cells_visited.is_empty())
3433 {
3434 Partition::Cell* const neighbour_cell = neighbour_cells_visited.pop();
3435 if(neighbour_cell->max_ival != neighbour_cell->length)
3436 value++;
3437 neighbour_cell->max_ival = 0;
3438 }
3439
3440 if((value > best_value) ||
3441 (value == best_value && cell->length > best_size))
3442 {
3443 best_value = value;
3444 best_size = cell->length;
3445 best_cell = cell;
3446 }
3447 }
3448 return best_cell;
3449 }
3450
3451
3452
3453
3454
3455
3456 /*------------------------------------------------------------------------
3457 *
3458 * Initialize the certificate size and memory
3459 *
3460 *-------------------------------------------------------------------------*/
3461
3462 void
initialize_certificate()3463 Digraph::initialize_certificate()
3464 {
3465 certificate_index = 0;
3466 certificate_current_path.clear();
3467 certificate_first_path.clear();
3468 certificate_best_path.clear();
3469 }
3470
3471
3472
3473 /*
3474 * Check whether perm is an automorphism.
3475 * Slow, mainly for debugging and validation purposes.
3476 */
3477 bool
is_automorphism(unsigned int * const perm)3478 Digraph::is_automorphism(unsigned int* const perm)
3479 {
3480 std::set<unsigned int, std::less<unsigned int> > edges1;
3481 std::set<unsigned int, std::less<unsigned int> > edges2;
3482
3483 #if defined(BLISS_CONSISTENCY_CHECKS)
3484 if(!is_permutation(get_nof_vertices(), perm))
3485 _INTERNAL_ERROR();
3486 #endif
3487
3488 for(unsigned int i = 0; i < get_nof_vertices(); i++)
3489 {
3490 Vertex& v1 = vertices[i];
3491 Vertex& v2 = vertices[perm[i]];
3492
3493 edges1.clear();
3494 for(std::vector<unsigned int>::iterator ei = v1.edges_in.begin();
3495 ei != v1.edges_in.end();
3496 ei++)
3497 edges1.insert(perm[*ei]);
3498 edges2.clear();
3499 for(std::vector<unsigned int>::iterator ei = v2.edges_in.begin();
3500 ei != v2.edges_in.end();
3501 ei++)
3502 edges2.insert(*ei);
3503 if(!(edges1 == edges2))
3504 return false;
3505
3506 edges1.clear();
3507 for(std::vector<unsigned int>::iterator ei = v1.edges_out.begin();
3508 ei != v1.edges_out.end();
3509 ei++)
3510 edges1.insert(perm[*ei]);
3511 edges2.clear();
3512 for(std::vector<unsigned int>::iterator ei = v2.edges_out.begin();
3513 ei != v2.edges_out.end();
3514 ei++)
3515 edges2.insert(*ei);
3516 if(!(edges1 == edges2))
3517 return false;
3518 }
3519
3520 return true;
3521 }
3522
3523 bool
is_automorphism(const std::vector<unsigned int> & perm) const3524 Digraph::is_automorphism(const std::vector<unsigned int>& perm) const
3525 {
3526
3527 if(!(perm.size() == get_nof_vertices() and is_permutation(perm)))
3528 return false;
3529
3530 std::set<unsigned int, std::less<unsigned int> > edges1;
3531 std::set<unsigned int, std::less<unsigned int> > edges2;
3532
3533 for(unsigned int i = 0; i < get_nof_vertices(); i++)
3534 {
3535 const Vertex& v1 = vertices[i];
3536 const Vertex& v2 = vertices[perm[i]];
3537
3538 edges1.clear();
3539 for(std::vector<unsigned int>::const_iterator ei = v1.edges_in.begin();
3540 ei != v1.edges_in.end();
3541 ei++)
3542 edges1.insert(perm[*ei]);
3543 edges2.clear();
3544 for(std::vector<unsigned int>::const_iterator ei = v2.edges_in.begin();
3545 ei != v2.edges_in.end();
3546 ei++)
3547 edges2.insert(*ei);
3548 if(!(edges1 == edges2))
3549 return false;
3550
3551 edges1.clear();
3552 for(std::vector<unsigned int>::const_iterator ei = v1.edges_out.begin();
3553 ei != v1.edges_out.end();
3554 ei++)
3555 edges1.insert(perm[*ei]);
3556 edges2.clear();
3557 for(std::vector<unsigned int>::const_iterator ei = v2.edges_out.begin();
3558 ei != v2.edges_out.end();
3559 ei++)
3560 edges2.insert(*ei);
3561 if(!(edges1 == edges2))
3562 return false;
3563 }
3564
3565 return true;
3566 }
3567
3568
3569
3570
3571 bool
nucr_find_first_component(const unsigned int level)3572 Digraph::nucr_find_first_component(const unsigned int level)
3573 {
3574
3575 cr_component.clear();
3576 cr_component_elements = 0;
3577
3578 /* Find first non-discrete cell in the component level */
3579 Partition::Cell* first_cell = p.first_nonsingleton_cell;
3580 while(first_cell)
3581 {
3582 if(p.cr_get_level(first_cell->first) == level)
3583 break;
3584 first_cell = first_cell->next_nonsingleton;
3585 }
3586
3587 /* The component is discrete, return false */
3588 if(!first_cell)
3589 return false;
3590
3591 std::vector<Partition::Cell*> component;
3592 first_cell->max_ival = 1;
3593 component.push_back(first_cell);
3594
3595 for(unsigned int i = 0; i < component.size(); i++)
3596 {
3597 Partition::Cell* const cell = component[i];
3598
3599 const Vertex& v = vertices[p.elements[cell->first]];
3600 std::vector<unsigned int>::const_iterator ei;
3601
3602 ei = v.edges_out.begin();
3603 for(unsigned int j = v.nof_edges_out(); j > 0; j--)
3604 {
3605 const unsigned int neighbour = *ei++;
3606 Partition::Cell* const neighbour_cell = p.get_cell(neighbour);
3607
3608 /* Skip unit neighbours */
3609 if(neighbour_cell->is_unit())
3610 continue;
3611 /* Already marked to be in the same component? */
3612 if(neighbour_cell->max_ival == 1)
3613 continue;
3614 /* Is the neighbour at the same component recursion level? */
3615 if(p.cr_get_level(neighbour_cell->first) != level)
3616 continue;
3617
3618 if(neighbour_cell->max_ival_count == 0)
3619 neighbour_heap.insert(neighbour_cell->first);
3620 neighbour_cell->max_ival_count++;
3621 }
3622 while(!neighbour_heap.is_empty())
3623 {
3624 const unsigned int start = neighbour_heap.remove();
3625 Partition::Cell* const neighbour_cell =
3626 p.get_cell(p.elements[start]);
3627
3628 /* Skip saturated neighbour cells */
3629 if(neighbour_cell->max_ival_count == neighbour_cell->length)
3630 {
3631 neighbour_cell->max_ival_count = 0;
3632 continue;
3633 }
3634 neighbour_cell->max_ival_count = 0;
3635 neighbour_cell->max_ival = 1;
3636 component.push_back(neighbour_cell);
3637 }
3638
3639 ei = v.edges_in.begin();
3640 for(unsigned int j = v.nof_edges_in(); j > 0; j--)
3641 {
3642 const unsigned int neighbour = *ei++;
3643
3644 Partition::Cell* const neighbour_cell = p.get_cell(neighbour);
3645
3646 /* Skip unit neighbours */
3647 if(neighbour_cell->is_unit())
3648 continue;
3649 /* Already marked to be in the same component? */
3650 if(neighbour_cell->max_ival == 1)
3651 continue;
3652 /* Is the neighbour at the same component recursion level? */
3653 if(p.cr_get_level(neighbour_cell->first) != level)
3654 continue;
3655
3656 if(neighbour_cell->max_ival_count == 0)
3657 neighbour_heap.insert(neighbour_cell->first);
3658 neighbour_cell->max_ival_count++;
3659 }
3660 while(!neighbour_heap.is_empty())
3661 {
3662 const unsigned int start = neighbour_heap.remove();
3663 Partition::Cell* const neighbour_cell =
3664 p.get_cell(p.elements[start]);
3665
3666 /* Skip saturated neighbour cells */
3667 if(neighbour_cell->max_ival_count == neighbour_cell->length)
3668 {
3669 neighbour_cell->max_ival_count = 0;
3670 continue;
3671 }
3672 neighbour_cell->max_ival_count = 0;
3673 neighbour_cell->max_ival = 1;
3674 component.push_back(neighbour_cell);
3675 }
3676 }
3677
3678 for(unsigned int i = 0; i < component.size(); i++)
3679 {
3680 Partition::Cell* const cell = component[i];
3681 cell->max_ival = 0;
3682 cr_component.push_back(cell->first);
3683 cr_component_elements += cell->length;
3684 }
3685
3686 return true;
3687 }
3688
3689
3690
3691
3692
3693 bool
nucr_find_first_component(const unsigned int level,std::vector<unsigned int> & component,unsigned int & component_elements,Partition::Cell * & sh_return)3694 Digraph::nucr_find_first_component(const unsigned int level,
3695 std::vector<unsigned int>& component,
3696 unsigned int& component_elements,
3697 Partition::Cell*& sh_return)
3698 {
3699
3700 component.clear();
3701 component_elements = 0;
3702 sh_return = 0;
3703 unsigned int sh_first = 0;
3704 unsigned int sh_size = 0;
3705 unsigned int sh_nuconn = 0;
3706
3707 /* Find first non-discrete cell in the component level */
3708 Partition::Cell* first_cell = p.first_nonsingleton_cell;
3709 while(first_cell)
3710 {
3711 if(p.cr_get_level(first_cell->first) == level)
3712 break;
3713 first_cell = first_cell->next_nonsingleton;
3714 }
3715
3716 if(!first_cell)
3717 {
3718 /* The component is discrete, return false */
3719 return false;
3720 }
3721
3722 std::vector<Partition::Cell*> comp;
3723 KStack<Partition::Cell*> neighbours;
3724 neighbours.init(get_nof_vertices());
3725
3726 first_cell->max_ival = 1;
3727 comp.push_back(first_cell);
3728
3729 for(unsigned int i = 0; i < comp.size(); i++)
3730 {
3731 Partition::Cell* const cell = comp[i];
3732
3733 unsigned int nuconn = 1;
3734
3735 const Vertex& v = vertices[p.elements[cell->first]];
3736 std::vector<unsigned int>::const_iterator ei;
3737
3738 /*| Phase 1: outgoing edges */
3739 ei = v.edges_out.begin();
3740 for(unsigned int j = v.nof_edges_out(); j > 0; j--)
3741 {
3742 const unsigned int neighbour = *ei++;
3743
3744 Partition::Cell* const neighbour_cell = p.get_cell(neighbour);
3745
3746 /* Skip unit neighbours */
3747 if(neighbour_cell->is_unit())
3748 continue;
3749 /* Is the neighbour at the same component recursion level? */
3750 //if(p.cr_get_level(neighbour_cell->first) != level)
3751 // continue;
3752 if(neighbour_cell->max_ival_count == 0)
3753 neighbours.push(neighbour_cell);
3754 neighbour_cell->max_ival_count++;
3755 }
3756 while(!neighbours.is_empty())
3757 {
3758 Partition::Cell* const neighbour_cell = neighbours.pop();
3759 /* Skip saturated neighbour cells */
3760 if(neighbour_cell->max_ival_count == neighbour_cell->length)
3761 {
3762 neighbour_cell->max_ival_count = 0;
3763 continue;
3764 }
3765 nuconn++;
3766 neighbour_cell->max_ival_count = 0;
3767 if(neighbour_cell->max_ival == 0) {
3768 comp.push_back(neighbour_cell);
3769 neighbour_cell->max_ival = 1;
3770 }
3771 }
3772
3773 /*| Phase 2: incoming edges */
3774 ei = v.edges_in.begin();
3775 for(unsigned int j = v.nof_edges_in(); j > 0; j--)
3776 {
3777 const unsigned int neighbour = *ei++;
3778 Partition::Cell* const neighbour_cell = p.get_cell(neighbour);
3779 /*| Skip unit neighbours */
3780 if(neighbour_cell->is_unit())
3781 continue;
3782 /* Is the neighbour at the same component recursion level? */
3783 //if(p.cr_get_level(neighbour_cell->first) != level)
3784 // continue;
3785 if(neighbour_cell->max_ival_count == 0)
3786 neighbours.push(neighbour_cell);
3787 neighbour_cell->max_ival_count++;
3788 }
3789 while(!neighbours.is_empty())
3790 {
3791 Partition::Cell* const neighbour_cell = neighbours.pop();
3792 /* Skip saturated neighbour cells */
3793 if(neighbour_cell->max_ival_count == neighbour_cell->length)
3794 {
3795 neighbour_cell->max_ival_count = 0;
3796 continue;
3797 }
3798 nuconn++;
3799 neighbour_cell->max_ival_count = 0;
3800 if(neighbour_cell->max_ival == 0) {
3801 comp.push_back(neighbour_cell);
3802 neighbour_cell->max_ival = 1;
3803 }
3804 }
3805
3806 /*| Phase 3: splitting heuristics */
3807 switch(sh) {
3808 case shs_f:
3809 if(sh_return == 0 or
3810 cell->first <= sh_first) {
3811 sh_return = cell;
3812 sh_first = cell->first;
3813 }
3814 break;
3815 case shs_fs:
3816 if(sh_return == 0 or
3817 cell->length < sh_size or
3818 (cell->length == sh_size and cell->first <= sh_first)) {
3819 sh_return = cell;
3820 sh_first = cell->first;
3821 sh_size = cell->length;
3822 }
3823 break;
3824 case shs_fl:
3825 if(sh_return == 0 or
3826 cell->length > sh_size or
3827 (cell->length == sh_size and cell->first <= sh_first)) {
3828 sh_return = cell;
3829 sh_first = cell->first;
3830 sh_size = cell->length;
3831 }
3832 break;
3833 case shs_fm:
3834 if(sh_return == 0 or
3835 nuconn > sh_nuconn or
3836 (nuconn == sh_nuconn and cell->first <= sh_first)) {
3837 sh_return = cell;
3838 sh_first = cell->first;
3839 sh_nuconn = nuconn;
3840 }
3841 break;
3842 case shs_fsm:
3843 if(sh_return == 0 or
3844 nuconn > sh_nuconn or
3845 (nuconn == sh_nuconn and
3846 (cell->length < sh_size or
3847 (cell->length == sh_size and cell->first <= sh_first)))) {
3848 sh_return = cell;
3849 sh_first = cell->first;
3850 sh_size = cell->length;
3851 sh_nuconn = nuconn;
3852 }
3853 break;
3854 case shs_flm:
3855 if(sh_return == 0 or
3856 nuconn > sh_nuconn or
3857 (nuconn == sh_nuconn and
3858 (cell->length > sh_size or
3859 (cell->length == sh_size and cell->first <= sh_first)))) {
3860 sh_return = cell;
3861 sh_first = cell->first;
3862 sh_size = cell->length;
3863 sh_nuconn = nuconn;
3864 }
3865 break;
3866 default:
3867 fatal_error("Internal error - unknown splitting heuristics");
3868 return 0;
3869 }
3870 }
3871 assert(sh_return);
3872
3873 for(unsigned int i = 0; i < comp.size(); i++)
3874 {
3875 Partition::Cell* const cell = comp[i];
3876 cell->max_ival = 0;
3877 component.push_back(cell->first);
3878 component_elements += cell->length;
3879 }
3880
3881 return true;
3882 }
3883
3884
3885
3886
3887 /*-------------------------------------------------------------------------
3888 *
3889 * Routines for undirected graphs
3890 *
3891 *-------------------------------------------------------------------------*/
3892
Vertex()3893 Graph::Vertex::Vertex()
3894 {
3895 color = 0;
3896 }
3897
3898
~Vertex()3899 Graph::Vertex::~Vertex()
3900 {
3901 ;
3902 }
3903
3904
3905 void
add_edge(const unsigned int other_vertex)3906 Graph::Vertex::add_edge(const unsigned int other_vertex)
3907 {
3908 edges.push_back(other_vertex);
3909 }
3910
3911
3912 void
remove_duplicate_edges(std::vector<bool> & tmp)3913 Graph::Vertex::remove_duplicate_edges(std::vector<bool>& tmp)
3914 {
3915 #if defined(BLISS_CONSISTENCY_CHECKS)
3916 /* Pre-conditions */
3917 for(unsigned int i = 0; i < tmp.size(); i++) assert(tmp[i] == false);
3918 #endif
3919 for(std::vector<unsigned int>::iterator iter = edges.begin();
3920 iter != edges.end(); )
3921 {
3922 const unsigned int dest_vertex = *iter;
3923 if(tmp[dest_vertex] == true)
3924 {
3925 /* A duplicate edge found! */
3926 iter = edges.erase(iter);
3927 }
3928 else
3929 {
3930 /* Not seen earlier, mark as seen */
3931 tmp[dest_vertex] = true;
3932 iter++;
3933 }
3934 }
3935
3936 /* Clear tmp */
3937 for(std::vector<unsigned int>::iterator iter = edges.begin();
3938 iter != edges.end();
3939 iter++)
3940 {
3941 tmp[*iter] = false;
3942 }
3943 #if defined(BLISS_CONSISTENCY_CHECKS)
3944 /* Post-conditions */
3945 for(unsigned int i = 0; i < tmp.size(); i++) assert(tmp[i] == false);
3946 #endif
3947 }
3948
3949
3950 /**
3951 * Sort the edges leaving the vertex according to
3952 * the vertex number of the other edge end.
3953 * Time complexity: O(e log(e)), where e is the number of edges
3954 * leaving the vertex.
3955 */
3956 void
sort_edges()3957 Graph::Vertex::sort_edges()
3958 {
3959 std::sort(edges.begin(), edges.end());
3960 }
3961
3962
3963
3964 /*-------------------------------------------------------------------------
3965 *
3966 * Constructor and destructor for undirected graphs
3967 *
3968 *-------------------------------------------------------------------------*/
3969
3970
Graph(const unsigned int nof_vertices)3971 Graph::Graph(const unsigned int nof_vertices)
3972 {
3973 vertices.resize(nof_vertices);
3974 sh = shs_flm;
3975 }
3976
3977
~Graph()3978 Graph::~Graph()
3979 {
3980 ;
3981 }
3982
3983
3984 unsigned int
add_vertex(const unsigned int color)3985 Graph::add_vertex(const unsigned int color)
3986 {
3987 const unsigned int vertex_num = vertices.size();
3988 vertices.resize(vertex_num + 1);
3989 vertices.back().color = color;
3990 return vertex_num;
3991 }
3992
3993
3994 void
add_edge(const unsigned int vertex1,const unsigned int vertex2)3995 Graph::add_edge(const unsigned int vertex1, const unsigned int vertex2)
3996 {
3997 //fprintf(stderr, "(%u,%u) ", vertex1, vertex2);
3998 vertices[vertex1].add_edge(vertex2);
3999 vertices[vertex2].add_edge(vertex1);
4000 }
4001
4002
4003 void
change_color(const unsigned int vertex,const unsigned int color)4004 Graph::change_color(const unsigned int vertex, const unsigned int color)
4005 {
4006 vertices[vertex].color = color;
4007 }
4008
4009
4010
4011
4012
4013 /*-------------------------------------------------------------------------
4014 *
4015 * Read graph in the DIMACS format.
4016 * Returns 0 if an error occurred.
4017 *
4018 *-------------------------------------------------------------------------*/
4019
4020 #if 0
4021 Graph*
4022 Graph::read_dimacs(FILE* const fp, FILE* const errstr)
4023 {
4024 Graph *g = 0;
4025 unsigned int nof_vertices;
4026 unsigned int nof_edges;
4027 unsigned int line_num = 1;
4028 int c;
4029
4030 const bool verbose = false;
4031 FILE* const verbstr = stdout;
4032
4033 /* Read comments and the problem definition line */
4034 while(1)
4035 {
4036 c = getc(fp);
4037 if(c == 'c')
4038 {
4039 /* A comment, ignore the rest of the line */
4040 while((c = getc(fp)) != '\n')
4041 {
4042 if(c == EOF)
4043 {
4044 if(errstr)
4045 fprintf(errstr,
4046 "error in line %u: not in DIMACS format\n",
4047 line_num);
4048 goto error_exit;
4049 }
4050 }
4051 line_num++;
4052 continue;
4053 }
4054 if(c == 'p')
4055 {
4056 /* The problem definition line */
4057 if(fscanf(fp, " edge %u %u\n", &nof_vertices, &nof_edges) != 2)
4058 {
4059 if(errstr)
4060 fprintf(errstr, "error in line %u: not in DIMACS format\n",
4061 line_num);
4062 goto error_exit;
4063 }
4064 line_num++;
4065 break;
4066 }
4067 if(errstr)
4068 fprintf(errstr, "error in line %u: not in DIMACS format\n", line_num);
4069 goto error_exit;
4070 }
4071
4072 if(nof_vertices <= 0)
4073 {
4074 if(errstr)
4075 fprintf(errstr, "error: no vertices\n");
4076 goto error_exit;
4077 }
4078 if(verbose)
4079 {
4080 fprintf(verbstr, "Instance has %d vertices and %d edges\n",
4081 nof_vertices, nof_edges);
4082 fflush(verbstr);
4083 }
4084
4085 g = new Graph(nof_vertices);
4086
4087 //
4088 // Read vertex colors
4089 //
4090 if(verbose)
4091 {
4092 fprintf(verbstr, "Reading vertex colors...\n");
4093 fflush(verbstr);
4094 }
4095 while(1)
4096 {
4097 c = getc(fp);
4098 if(c != 'n')
4099 {
4100 ungetc(c, fp);
4101 break;
4102 }
4103 ungetc(c, fp);
4104 unsigned int vertex;
4105 unsigned int color;
4106 if(fscanf(fp, "n %u %u\n", &vertex, &color) != 2)
4107 {
4108 if(errstr)
4109 fprintf(errstr, "error in line %u: not in DIMACS format\n",
4110 line_num);
4111 goto error_exit;
4112 }
4113 if(!((vertex >= 1) && (vertex <= nof_vertices)))
4114 {
4115 if(errstr)
4116 fprintf(errstr,
4117 "error in line %u: vertex %u not in range [1,...,%u]\n",
4118 line_num, vertex, nof_vertices);
4119 goto error_exit;
4120 }
4121 line_num++;
4122 g->change_color(vertex - 1, color);
4123 }
4124 if(verbose)
4125 {
4126 fprintf(verbstr, "Done\n");
4127 fflush(verbstr);
4128 }
4129
4130 //
4131 // Read edges
4132 //
4133 if(verbose)
4134 {
4135 fprintf(verbstr, "Reading edges...\n");
4136 fflush(verbstr);
4137 }
4138 for(unsigned i = 0; i < nof_edges; i++)
4139 {
4140 unsigned int from, to;
4141 if(fscanf(fp, "e %u %u\n", &from, &to) != 2)
4142 {
4143 if(errstr)
4144 fprintf(errstr, "error in line %u: not in DIMACS format\n",
4145 line_num);
4146 goto error_exit;
4147 }
4148 if(!((from >= 1) && (from <= nof_vertices)))
4149 {
4150 if(errstr)
4151 fprintf(errstr,
4152 "error in line %u: vertex %u not in range [1,...,%u]\n",
4153 line_num, from, nof_vertices);
4154 goto error_exit;
4155 }
4156 if(!((to >= 1) && (to <= nof_vertices)))
4157 {
4158 if(errstr)
4159 fprintf(errstr,
4160 "error in line %u: vertex %u not in range [1,...,%u]\n",
4161 line_num, to, nof_vertices);
4162 goto error_exit;
4163 }
4164 line_num++;
4165 g->add_edge(from-1, to-1);
4166 }
4167 if(verbose)
4168 {
4169 fprintf(verbstr, "Done\n");
4170 fflush(verbstr);
4171 }
4172
4173 return g;
4174
4175 error_exit:
4176 if(g)
4177 delete g;
4178 return 0;
4179
4180 }
4181
4182
4183 void
4184 Graph::write_dimacs(FILE* const fp)
4185 {
4186 remove_duplicate_edges();
4187 sort_edges();
4188
4189 /* First count the total number of edges */
4190 unsigned int nof_edges = 0;
4191 for(unsigned int i = 0; i < get_nof_vertices(); i++)
4192 {
4193 Vertex &v = vertices[i];
4194 for(std::vector<unsigned int>::const_iterator ei = v.edges.begin();
4195 ei != v.edges.end();
4196 ei++)
4197 {
4198 const unsigned int dest_i = *ei;
4199 if(dest_i < i)
4200 continue;
4201 nof_edges++;
4202 }
4203 }
4204
4205 /* Output the "header" line */
4206 fprintf(fp, "p edge %u %u\n", get_nof_vertices(), nof_edges);
4207
4208 /* Print the color of each vertex */
4209 for(unsigned int i = 0; i < get_nof_vertices(); i++)
4210 {
4211 Vertex &v = vertices[i];
4212 fprintf(fp, "n %u %u\n", i+1, v.color);
4213 /*
4214 if(v.color != 0)
4215 {
4216 fprintf(fp, "n %u %u\n", i+1, v.color);
4217 }
4218 */
4219 }
4220
4221 /* Print the edges */
4222 for(unsigned int i = 0; i < get_nof_vertices(); i++)
4223 {
4224 Vertex &v = vertices[i];
4225 for(std::vector<unsigned int>::const_iterator ei = v.edges.begin();
4226 ei != v.edges.end();
4227 ei++)
4228 {
4229 const unsigned int dest_i = *ei;
4230 if(dest_i < i)
4231 continue;
4232 fprintf(fp, "e %u %u\n", i+1, dest_i+1);
4233 }
4234 }
4235 }
4236 #endif
4237
4238
4239 void
sort_edges()4240 Graph::sort_edges()
4241 {
4242 for(unsigned int i = 0; i < get_nof_vertices(); i++)
4243 vertices[i].sort_edges();
4244 }
4245
4246
4247 int
cmp(Graph & other)4248 Graph::cmp(Graph& other)
4249 {
4250 /* Compare the numbers of vertices */
4251 if(get_nof_vertices() < other.get_nof_vertices())
4252 return -1;
4253 if(get_nof_vertices() > other.get_nof_vertices())
4254 return 1;
4255 /* Compare vertex colors */
4256 for(unsigned int i = 0; i < get_nof_vertices(); i++)
4257 {
4258 if(vertices[i].color < other.vertices[i].color)
4259 return -1;
4260 if(vertices[i].color > other.vertices[i].color)
4261 return 1;
4262 }
4263 /* Compare vertex degrees */
4264 remove_duplicate_edges();
4265 other.remove_duplicate_edges();
4266 for(unsigned int i = 0; i < get_nof_vertices(); i++)
4267 {
4268 if(vertices[i].nof_edges() < other.vertices[i].nof_edges())
4269 return -1;
4270 if(vertices[i].nof_edges() > other.vertices[i].nof_edges())
4271 return 1;
4272 }
4273 /* Compare edges */
4274 for(unsigned int i = 0; i < get_nof_vertices(); i++)
4275 {
4276 Vertex &v1 = vertices[i];
4277 Vertex &v2 = other.vertices[i];
4278 v1.sort_edges();
4279 v2.sort_edges();
4280 std::vector<unsigned int>::const_iterator ei1 = v1.edges.begin();
4281 std::vector<unsigned int>::const_iterator ei2 = v2.edges.begin();
4282 while(ei1 != v1.edges.end())
4283 {
4284 if(*ei1 < *ei2)
4285 return -1;
4286 if(*ei1 > *ei2)
4287 return 1;
4288 ei1++;
4289 ei2++;
4290 }
4291 }
4292 return 0;
4293 }
4294
4295
4296 Graph*
permute(const std::vector<unsigned int> & perm) const4297 Graph::permute(const std::vector<unsigned int>& perm) const
4298 {
4299 #if defined(BLISS_CONSISTENCY_CHECKS)
4300 #endif
4301
4302 Graph* const g = new Graph(get_nof_vertices());
4303 for(unsigned int i = 0; i < get_nof_vertices(); i++)
4304 {
4305 const Vertex& v = vertices[i];
4306 Vertex& permuted_v = g->vertices[perm[i]];
4307 permuted_v.color = v.color;
4308 for(std::vector<unsigned int>::const_iterator ei = v.edges.begin();
4309 ei != v.edges.end();
4310 ei++)
4311 {
4312 const unsigned int dest_v = *ei;
4313 permuted_v.add_edge(perm[dest_v]);
4314 }
4315 permuted_v.sort_edges();
4316 }
4317 return g;
4318 }
4319
4320 Graph*
permute(const unsigned int * perm) const4321 Graph::permute(const unsigned int* perm) const
4322 {
4323 #if defined(BLISS_CONSISTENCY_CHECKS)
4324 if(!is_permutation(get_nof_vertices(), perm))
4325 _INTERNAL_ERROR();
4326 #endif
4327
4328 Graph* const g = new Graph(get_nof_vertices());
4329 for(unsigned int i = 0; i < get_nof_vertices(); i++)
4330 {
4331 const Vertex& v = vertices[i];
4332 Vertex& permuted_v = g->vertices[perm[i]];
4333 permuted_v.color = v.color;
4334 for(std::vector<unsigned int>::const_iterator ei = v.edges.begin();
4335 ei != v.edges.end();
4336 ei++)
4337 {
4338 const unsigned int dest_v = *ei;
4339 permuted_v.add_edge(perm[dest_v]);
4340 }
4341 permuted_v.sort_edges();
4342 }
4343 return g;
4344 }
4345
4346
4347
4348
4349
4350 /*-------------------------------------------------------------------------
4351 *
4352 * Print graph in graphviz format
4353 *
4354 *-------------------------------------------------------------------------*/
4355
4356
4357 #if 0
4358 void
4359 Graph::write_dot(const char* const filename)
4360 {
4361 FILE *fp = fopen(filename, "w");
4362 if(fp)
4363 {
4364 write_dot(fp);
4365 fclose(fp);
4366 }
4367 }
4368
4369 void
4370 Graph::write_dot(FILE* const fp)
4371 {
4372 remove_duplicate_edges();
4373
4374 fprintf(fp, "graph g {\n");
4375
4376 unsigned int vnum = 0;
4377 for(std::vector<Vertex>::iterator vi = vertices.begin();
4378 vi != vertices.end();
4379 vi++, vnum++)
4380 {
4381 Vertex& v = *vi;
4382 fprintf(fp, "v%u [label=\"%u:%u\"];\n", vnum, vnum, v.color);
4383 for(std::vector<unsigned int>::const_iterator ei = v.edges.begin();
4384 ei != v.edges.end();
4385 ei++)
4386 {
4387 const unsigned int vnum2 = *ei;
4388 if(vnum2 > vnum)
4389 fprintf(fp, "v%u -- v%u\n", vnum, vnum2);
4390 }
4391 }
4392
4393 fprintf(fp, "}\n");
4394 }
4395 #endif
4396
4397
4398
4399
4400
4401
4402 /*-------------------------------------------------------------------------
4403 *
4404 * Get a hash value for the graph.
4405 *
4406 *-------------------------------------------------------------------------*/
4407
4408 unsigned int
get_hash()4409 Graph::get_hash()
4410 {
4411 remove_duplicate_edges();
4412 sort_edges();
4413
4414 UintSeqHash h;
4415
4416 h.update(get_nof_vertices());
4417
4418 /* Hash the color of each vertex */
4419 for(unsigned int i = 0; i < get_nof_vertices(); i++)
4420 {
4421 h.update(vertices[i].color);
4422 }
4423
4424 /* Hash the edges */
4425 for(unsigned int i = 0; i < get_nof_vertices(); i++)
4426 {
4427 Vertex &v = vertices[i];
4428 for(std::vector<unsigned int>::const_iterator ei = v.edges.begin();
4429 ei != v.edges.end();
4430 ei++)
4431 {
4432 const unsigned int dest_i = *ei;
4433 if(dest_i < i)
4434 continue;
4435 h.update(i);
4436 h.update(dest_i);
4437 }
4438 }
4439
4440 return h.get_value();
4441 }
4442
4443
4444
4445
4446
4447 void
remove_duplicate_edges()4448 Graph::remove_duplicate_edges()
4449 {
4450 std::vector<bool> tmp(vertices.size(), false);
4451
4452 for(std::vector<Vertex>::iterator vi = vertices.begin();
4453 vi != vertices.end();
4454 vi++)
4455 {
4456 #if defined(BLISS_EXPENSIVE_CONSISTENCY_CHECKS)
4457 for(unsigned int i = 0; i < tmp.size(); i++) assert(tmp[i] == false);
4458 #endif
4459 (*vi).remove_duplicate_edges(tmp);
4460 }
4461 }
4462
4463
4464
4465
4466
4467 /*-------------------------------------------------------------------------
4468 *
4469 * Partition independent invariants
4470 *
4471 *-------------------------------------------------------------------------*/
4472
4473 /*
4474 * Return the color of the vertex.
4475 * Time complexity: O(1)
4476 */
4477 unsigned int
vertex_color_invariant(const Graph * const g,const unsigned int v)4478 Graph::vertex_color_invariant(const Graph* const g, const unsigned int v)
4479 {
4480 return g->vertices[v].color;
4481 }
4482
4483 /*
4484 * Return the degree of the vertex.
4485 * Time complexity: O(1)
4486 */
4487 unsigned int
degree_invariant(const Graph * const g,const unsigned int v)4488 Graph::degree_invariant(const Graph* const g, const unsigned int v)
4489 {
4490 return g->vertices[v].nof_edges();
4491 }
4492
4493 /*
4494 * Return 1 if the vertex v has a self-loop, 0 otherwise
4495 * Time complexity: O(E_v), where E_v is the number of edges leaving v
4496 */
4497 unsigned int
selfloop_invariant(const Graph * const g,const unsigned int v)4498 Graph::selfloop_invariant(const Graph* const g, const unsigned int v)
4499 {
4500 const Vertex& vertex = g->vertices[v];
4501 for(std::vector<unsigned int>::const_iterator ei = vertex.edges.begin();
4502 ei != vertex.edges.end();
4503 ei++)
4504 {
4505 if(*ei == v)
4506 return 1;
4507 }
4508 return 0;
4509 }
4510
4511
4512
4513
4514
4515
4516 /*-------------------------------------------------------------------------
4517 *
4518 * Refine the partition p according to a partition independent invariant
4519 *
4520 *-------------------------------------------------------------------------*/
4521
4522 bool
refine_according_to_invariant(unsigned int (* inv)(const Graph * const g,const unsigned int v))4523 Graph::refine_according_to_invariant(unsigned int (*inv)(const Graph* const g,
4524 const unsigned int v))
4525 {
4526 bool refined = false;
4527
4528 for(Partition::Cell* cell = p.first_nonsingleton_cell; cell; )
4529 {
4530
4531 Partition::Cell* const next_cell = cell->next_nonsingleton;
4532
4533 const unsigned int* ep = p.elements + cell->first;
4534 for(unsigned int i = cell->length; i > 0; i--, ep++)
4535 {
4536 const unsigned int ival = inv(this, *ep);
4537 p.invariant_values[*ep] = ival;
4538 if(ival > cell->max_ival)
4539 {
4540 cell->max_ival = ival;
4541 cell->max_ival_count = 1;
4542 }
4543 else if(ival == cell->max_ival)
4544 {
4545 cell->max_ival_count++;
4546 }
4547 }
4548 Partition::Cell* const last_new_cell = p.zplit_cell(cell, true);
4549 refined |= (last_new_cell != cell);
4550 cell = next_cell;
4551 }
4552
4553 return refined;
4554 }
4555
4556
4557
4558
4559
4560
4561
4562
4563
4564
4565
4566
4567 /*-------------------------------------------------------------------------
4568 *
4569 * Split the neighbourhood of a cell according to the equitable invariant
4570 *
4571 *-------------------------------------------------------------------------*/
4572
4573 bool
split_neighbourhood_of_cell(Partition::Cell * const cell)4574 Graph::split_neighbourhood_of_cell(Partition::Cell* const cell)
4575 {
4576
4577
4578 const bool was_equal_to_first = refine_equal_to_first;
4579
4580 if(compute_eqref_hash)
4581 {
4582 eqref_hash.update(cell->first);
4583 eqref_hash.update(cell->length);
4584 }
4585
4586 const unsigned int* ep = p.elements + cell->first;
4587 for(unsigned int i = cell->length; i > 0; i--)
4588 {
4589 const Vertex& v = vertices[*ep++];
4590
4591 std::vector<unsigned int>::const_iterator ei = v.edges.begin();
4592 for(unsigned int j = v.nof_edges(); j != 0; j--)
4593 {
4594 const unsigned int dest_vertex = *ei++;
4595 Partition::Cell * const neighbour_cell = p.get_cell(dest_vertex);
4596 if(neighbour_cell->is_unit())
4597 continue;
4598 const unsigned int ival = ++p.invariant_values[dest_vertex];
4599 if(ival > neighbour_cell->max_ival)
4600 {
4601 neighbour_cell->max_ival = ival;
4602 neighbour_cell->max_ival_count = 1;
4603 if(ival == 1) {
4604 neighbour_heap.insert(neighbour_cell->first);
4605 }
4606 }
4607 else if(ival == neighbour_cell->max_ival) {
4608 neighbour_cell->max_ival_count++;
4609 }
4610 }
4611 }
4612
4613 while(!neighbour_heap.is_empty())
4614 {
4615 const unsigned int start = neighbour_heap.remove();
4616 Partition::Cell * const neighbour_cell = p.get_cell(p.elements[start]);
4617
4618 if(compute_eqref_hash)
4619 {
4620 eqref_hash.update(neighbour_cell->first);
4621 eqref_hash.update(neighbour_cell->length);
4622 eqref_hash.update(neighbour_cell->max_ival);
4623 eqref_hash.update(neighbour_cell->max_ival_count);
4624 }
4625
4626
4627 Partition::Cell* const last_new_cell = p.zplit_cell(neighbour_cell, true);
4628
4629 /* Update certificate and hash if needed */
4630 const Partition::Cell* c = neighbour_cell;
4631 while(1)
4632 {
4633 if(in_search)
4634 {
4635 /* Build certificate */
4636 cert_add_redundant(CERT_SPLIT, c->first, c->length);
4637 /* No need to continue? */
4638 if(refine_compare_certificate and
4639 (refine_equal_to_first == false) and
4640 (refine_cmp_to_best < 0))
4641 goto worse_exit;
4642 }
4643 if(compute_eqref_hash)
4644 {
4645 eqref_hash.update(c->first);
4646 eqref_hash.update(c->length);
4647 }
4648 if(c == last_new_cell)
4649 break;
4650 c = c->next;
4651 }
4652 }
4653
4654 if(refine_compare_certificate and
4655 (refine_equal_to_first == false) and
4656 (refine_cmp_to_best < 0))
4657 return true;
4658
4659 return false;
4660
4661 worse_exit:
4662 /* Clear neighbour heap */
4663 UintSeqHash rest;
4664 while(!neighbour_heap.is_empty())
4665 {
4666 const unsigned int start = neighbour_heap.remove();
4667 Partition::Cell * const neighbour_cell = p.get_cell(p.elements[start]);
4668 if(opt_use_failure_recording and was_equal_to_first)
4669 {
4670 rest.update(neighbour_cell->first);
4671 rest.update(neighbour_cell->length);
4672 rest.update(neighbour_cell->max_ival);
4673 rest.update(neighbour_cell->max_ival_count);
4674 }
4675 neighbour_cell->max_ival = 0;
4676 neighbour_cell->max_ival_count = 0;
4677 p.clear_ivs(neighbour_cell);
4678 }
4679 if(opt_use_failure_recording and was_equal_to_first)
4680 {
4681 for(unsigned int i = p.splitting_queue.size(); i > 0; i--)
4682 {
4683 Partition::Cell* const cell = p.splitting_queue.pop_front();
4684 rest.update(cell->first);
4685 rest.update(cell->length);
4686 p.splitting_queue.push_back(cell);
4687 }
4688 rest.update(failure_recording_fp_deviation);
4689 failure_recording_fp_deviation = rest.get_value();
4690 }
4691
4692 return true;
4693 }
4694
4695
4696
4697 bool
split_neighbourhood_of_unit_cell(Partition::Cell * const unit_cell)4698 Graph::split_neighbourhood_of_unit_cell(Partition::Cell* const unit_cell)
4699 {
4700
4701
4702 const bool was_equal_to_first = refine_equal_to_first;
4703
4704 if(compute_eqref_hash)
4705 {
4706 eqref_hash.update(0x87654321);
4707 eqref_hash.update(unit_cell->first);
4708 eqref_hash.update(1);
4709 }
4710
4711 const Vertex& v = vertices[p.elements[unit_cell->first]];
4712
4713 std::vector<unsigned int>::const_iterator ei = v.edges.begin();
4714 for(unsigned int j = v.nof_edges(); j > 0; j--)
4715 {
4716 const unsigned int dest_vertex = *ei++;
4717 Partition::Cell * const neighbour_cell = p.get_cell(dest_vertex);
4718
4719 if(neighbour_cell->is_unit()) {
4720 if(in_search) {
4721 /* Remember neighbour in order to generate certificate */
4722 neighbour_heap.insert(neighbour_cell->first);
4723 }
4724 continue;
4725 }
4726 if(neighbour_cell->max_ival_count == 0)
4727 {
4728 neighbour_heap.insert(neighbour_cell->first);
4729 }
4730 neighbour_cell->max_ival_count++;
4731
4732 unsigned int * const swap_position =
4733 p.elements + neighbour_cell->first + neighbour_cell->length -
4734 neighbour_cell->max_ival_count;
4735 *p.in_pos[dest_vertex] = *swap_position;
4736 p.in_pos[*swap_position] = p.in_pos[dest_vertex];
4737 *swap_position = dest_vertex;
4738 p.in_pos[dest_vertex] = swap_position;
4739 }
4740
4741 while(!neighbour_heap.is_empty())
4742 {
4743 const unsigned int start = neighbour_heap.remove();
4744 Partition::Cell* neighbour_cell = p.get_cell(p.elements[start]);
4745
4746 #if defined(BLISS_CONSISTENCY_CHECKS)
4747 if(neighbour_cell->is_unit()) {
4748 } else {
4749 }
4750 #endif
4751
4752 if(compute_eqref_hash)
4753 {
4754 eqref_hash.update(neighbour_cell->first);
4755 eqref_hash.update(neighbour_cell->length);
4756 eqref_hash.update(neighbour_cell->max_ival_count);
4757 }
4758
4759 if(neighbour_cell->length > 1 and
4760 neighbour_cell->max_ival_count != neighbour_cell->length)
4761 {
4762 Partition::Cell * const new_cell =
4763 p.aux_split_in_two(neighbour_cell,
4764 neighbour_cell->length -
4765 neighbour_cell->max_ival_count);
4766 unsigned int *ep = p.elements + new_cell->first;
4767 unsigned int * const lp = p.elements+new_cell->first+new_cell->length;
4768 while(ep < lp)
4769 {
4770 p.element_to_cell_map[*ep] = new_cell;
4771 ep++;
4772 }
4773 neighbour_cell->max_ival_count = 0;
4774
4775
4776 if(compute_eqref_hash)
4777 {
4778 /* Update hash */
4779 eqref_hash.update(neighbour_cell->first);
4780 eqref_hash.update(neighbour_cell->length);
4781 eqref_hash.update(0);
4782 eqref_hash.update(new_cell->first);
4783 eqref_hash.update(new_cell->length);
4784 eqref_hash.update(1);
4785 }
4786
4787 /* Add cells in splitting_queue */
4788 if(neighbour_cell->is_in_splitting_queue()) {
4789 /* Both cells must be included in splitting_queue in order
4790 to ensure refinement into equitable partition */
4791 p.splitting_queue_add(new_cell);
4792 } else {
4793 Partition::Cell *min_cell, *max_cell;
4794 if(neighbour_cell->length <= new_cell->length) {
4795 min_cell = neighbour_cell;
4796 max_cell = new_cell;
4797 } else {
4798 min_cell = new_cell;
4799 max_cell = neighbour_cell;
4800 }
4801 /* Put the smaller cell in splitting_queue */
4802 p.splitting_queue_add(min_cell);
4803 if(max_cell->is_unit()) {
4804 /* Put the "larger" cell also in splitting_queue */
4805 p.splitting_queue_add(max_cell);
4806 }
4807 }
4808 /* Update pointer for certificate generation */
4809 neighbour_cell = new_cell;
4810 }
4811 else
4812 {
4813 /* neighbour_cell->length == 1 ||
4814 neighbour_cell->max_ival_count == neighbour_cell->length */
4815 neighbour_cell->max_ival_count = 0;
4816 }
4817
4818 /*
4819 * Build certificate if required
4820 */
4821 if(in_search)
4822 {
4823 for(unsigned int i = neighbour_cell->first,
4824 j = neighbour_cell->length;
4825 j > 0;
4826 j--, i++)
4827 {
4828 /* Build certificate */
4829 cert_add(CERT_EDGE, unit_cell->first, i);
4830 /* No need to continue? */
4831 if(refine_compare_certificate and
4832 (refine_equal_to_first == false) and
4833 (refine_cmp_to_best < 0))
4834 goto worse_exit;
4835 }
4836 } /* if(in_search) */
4837 } /* while(!neighbour_heap.is_empty()) */
4838
4839 if(refine_compare_certificate and
4840 (refine_equal_to_first == false) and
4841 (refine_cmp_to_best < 0))
4842 return true;
4843
4844 return false;
4845
4846 worse_exit:
4847 /* Clear neighbour heap */
4848 UintSeqHash rest;
4849 while(!neighbour_heap.is_empty())
4850 {
4851 const unsigned int start = neighbour_heap.remove();
4852 Partition::Cell * const neighbour_cell = p.get_cell(p.elements[start]);
4853 if(opt_use_failure_recording and was_equal_to_first)
4854 {
4855 rest.update(neighbour_cell->first);
4856 rest.update(neighbour_cell->length);
4857 rest.update(neighbour_cell->max_ival_count);
4858 }
4859 neighbour_cell->max_ival_count = 0;
4860 }
4861 if(opt_use_failure_recording and was_equal_to_first)
4862 {
4863 rest.update(failure_recording_fp_deviation);
4864 failure_recording_fp_deviation = rest.get_value();
4865 }
4866 return true;
4867 }
4868
4869
4870
4871
4872
4873
4874
4875
4876
4877 /*-------------------------------------------------------------------------
4878 *
4879 * Check whether the current partition p is equitable.
4880 * Performance: very slow, use only for debugging purposes.
4881 *
4882 *-------------------------------------------------------------------------*/
4883
is_equitable() const4884 bool Graph::is_equitable() const
4885 {
4886 const unsigned int N = get_nof_vertices();
4887 if(N == 0)
4888 return true;
4889
4890 std::vector<unsigned int> first_count = std::vector<unsigned int>(N, 0);
4891 std::vector<unsigned int> other_count = std::vector<unsigned int>(N, 0);
4892
4893 for(Partition::Cell *cell = p.first_cell; cell; cell = cell->next)
4894 {
4895 if(cell->is_unit())
4896 continue;
4897
4898 unsigned int *ep = p.elements + cell->first;
4899 const Vertex &first_vertex = vertices[*ep++];
4900
4901 /* Count how many edges lead from the first vertex to
4902 * the neighbouring cells */
4903 for(std::vector<unsigned int>::const_iterator ei =
4904 first_vertex.edges.begin();
4905 ei != first_vertex.edges.end();
4906 ei++)
4907 {
4908 first_count[p.get_cell(*ei)->first]++;
4909 }
4910
4911 /* Count and compare to the edges of the other vertices */
4912 for(unsigned int i = cell->length; i > 1; i--)
4913 {
4914 const Vertex &vertex = vertices[*ep++];
4915 for(std::vector<unsigned int>::const_iterator ei =
4916 vertex.edges.begin();
4917 ei != vertex.edges.end();
4918 ei++)
4919 {
4920 other_count[p.get_cell(*ei)->first]++;
4921 }
4922 for(Partition::Cell *cell2 = p.first_cell;
4923 cell2;
4924 cell2 = cell2->next)
4925 {
4926 if(first_count[cell2->first] != other_count[cell2->first])
4927 {
4928 /* Not equitable */
4929 return false;
4930 }
4931 other_count[cell2->first] = 0;
4932 }
4933 }
4934 /* Reset first_count */
4935 for(unsigned int i = 0; i < N; i++)
4936 first_count[i] = 0;
4937 }
4938 return true;
4939 }
4940
4941
4942
4943
4944
4945 /*-------------------------------------------------------------------------
4946 *
4947 * Build the initial equitable partition
4948 *
4949 *-------------------------------------------------------------------------*/
4950
make_initial_equitable_partition()4951 void Graph::make_initial_equitable_partition()
4952 {
4953 refine_according_to_invariant(&vertex_color_invariant);
4954 p.splitting_queue_clear();
4955 //p.print_signature(stderr); fprintf(stderr, "\n");
4956
4957 refine_according_to_invariant(&selfloop_invariant);
4958 p.splitting_queue_clear();
4959 //p.print_signature(stderr); fprintf(stderr, "\n");
4960
4961 refine_according_to_invariant(°ree_invariant);
4962 p.splitting_queue_clear();
4963 //p.print_signature(stderr); fprintf(stderr, "\n");
4964
4965 refine_to_equitable();
4966 //p.print_signature(stderr); fprintf(stderr, "\n");
4967
4968 }
4969
4970
4971
4972
4973
4974
4975
4976 /*-------------------------------------------------------------------------
4977 *
4978 * Find the next cell to be splitted
4979 *
4980 *-------------------------------------------------------------------------*/
4981
4982
4983 Partition::Cell*
find_next_cell_to_be_splitted(Partition::Cell * cell)4984 Graph::find_next_cell_to_be_splitted(Partition::Cell* cell)
4985 {
4986 switch(sh) {
4987 case shs_f: return sh_first();
4988 case shs_fs: return sh_first_smallest();
4989 case shs_fl: return sh_first_largest();
4990 case shs_fm: return sh_first_max_neighbours();
4991 case shs_fsm: return sh_first_smallest_max_neighbours();
4992 case shs_flm: return sh_first_largest_max_neighbours();
4993 default:
4994 fatal_error("Internal error - unknown splitting heuristics");
4995 return 0;
4996 }
4997 }
4998
4999 /** \internal
5000 * A splitting heuristic.
5001 * Returns the first nonsingleton cell in the current partition.
5002 */
5003 Partition::Cell*
sh_first()5004 Graph::sh_first()
5005 {
5006 Partition::Cell* best_cell = 0;
5007 for(Partition::Cell* cell = p.first_nonsingleton_cell;
5008 cell;
5009 cell = cell->next_nonsingleton)
5010 {
5011 if(opt_use_comprec and p.cr_get_level(cell->first) != cr_level)
5012 continue;
5013 best_cell = cell;
5014 break;
5015 }
5016 return best_cell;
5017 }
5018
5019 /** \internal
5020 * A splitting heuristic.
5021 * Returns the first smallest nonsingleton cell in the current partition.
5022 */
5023 Partition::Cell*
sh_first_smallest()5024 Graph::sh_first_smallest()
5025 {
5026 Partition::Cell* best_cell = 0;
5027 unsigned int best_size = UINT_MAX;
5028 for(Partition::Cell* cell = p.first_nonsingleton_cell;
5029 cell;
5030 cell = cell->next_nonsingleton)
5031 {
5032 if(opt_use_comprec and p.cr_get_level(cell->first) != cr_level)
5033 continue;
5034 if(cell->length < best_size)
5035 {
5036 best_size = cell->length;
5037 best_cell = cell;
5038 }
5039 }
5040 return best_cell;
5041 }
5042
5043 /** \internal
5044 * A splitting heuristic.
5045 * Returns the first largest nonsingleton cell in the current partition.
5046 */
5047 Partition::Cell*
sh_first_largest()5048 Graph::sh_first_largest()
5049 {
5050 Partition::Cell* best_cell = 0;
5051 unsigned int best_size = 0;
5052 for(Partition::Cell* cell = p.first_nonsingleton_cell;
5053 cell;
5054 cell = cell->next_nonsingleton)
5055 {
5056 if(opt_use_comprec and p.cr_get_level(cell->first) != cr_level)
5057 continue;
5058 if(cell->length > best_size)
5059 {
5060 best_size = cell->length;
5061 best_cell = cell;
5062 }
5063 }
5064 return best_cell;
5065 }
5066
5067 /** \internal
5068 * A splitting heuristic.
5069 * Returns the first nonsingleton cell with max number of neighbouring
5070 * nonsingleton cells.
5071 * Assumes that the partition p is equitable.
5072 * Assumes that the max_ival fields of the cells are all 0.
5073 */
5074 Partition::Cell*
sh_first_max_neighbours()5075 Graph::sh_first_max_neighbours()
5076 {
5077 Partition::Cell* best_cell = 0;
5078 int best_value = -1;
5079 KStack<Partition::Cell*> neighbour_cells_visited;
5080 neighbour_cells_visited.init(get_nof_vertices());
5081 for(Partition::Cell* cell = p.first_nonsingleton_cell;
5082 cell;
5083 cell = cell->next_nonsingleton)
5084 {
5085 if(opt_use_comprec and p.cr_get_level(cell->first) != cr_level)
5086 continue;
5087 const Vertex& v = vertices[p.elements[cell->first]];
5088 std::vector<unsigned int>::const_iterator ei = v.edges.begin();
5089 for(unsigned int j = v.nof_edges(); j > 0; j--)
5090 {
5091 Partition::Cell * const neighbour_cell = p.get_cell(*ei++);
5092 if(neighbour_cell->is_unit())
5093 continue;
5094 neighbour_cell->max_ival++;
5095 if(neighbour_cell->max_ival == 1)
5096 neighbour_cells_visited.push(neighbour_cell);
5097 }
5098 int value = 0;
5099 while(!neighbour_cells_visited.is_empty())
5100 {
5101 Partition::Cell* const neighbour_cell = neighbour_cells_visited.pop();
5102 if(neighbour_cell->max_ival != neighbour_cell->length)
5103 value++;
5104 neighbour_cell->max_ival = 0;
5105 }
5106 if(value > best_value)
5107 {
5108 best_value = value;
5109 best_cell = cell;
5110 }
5111 }
5112 return best_cell;
5113 }
5114
5115 /** \internal
5116 * A splitting heuristic.
5117 * Returns the first smallest nonsingleton cell with max number of neighbouring
5118 * nonsingleton cells.
5119 * Assumes that the partition p is equitable.
5120 * Assumes that the max_ival fields of the cells are all 0.
5121 */
5122 Partition::Cell*
sh_first_smallest_max_neighbours()5123 Graph::sh_first_smallest_max_neighbours()
5124 {
5125 Partition::Cell* best_cell = 0;
5126 int best_value = -1;
5127 unsigned int best_size = UINT_MAX;
5128 KStack<Partition::Cell*> neighbour_cells_visited;
5129 neighbour_cells_visited.init(get_nof_vertices());
5130 for(Partition::Cell* cell = p.first_nonsingleton_cell;
5131 cell;
5132 cell = cell->next_nonsingleton)
5133 {
5134
5135 if(opt_use_comprec and p.cr_get_level(cell->first) != cr_level)
5136 continue;
5137
5138 const Vertex& v = vertices[p.elements[cell->first]];
5139 std::vector<unsigned int>::const_iterator ei = v.edges.begin();
5140 for(unsigned int j = v.nof_edges(); j > 0; j--)
5141 {
5142 Partition::Cell* const neighbour_cell = p.get_cell(*ei++);
5143 if(neighbour_cell->is_unit())
5144 continue;
5145 neighbour_cell->max_ival++;
5146 if(neighbour_cell->max_ival == 1)
5147 neighbour_cells_visited.push(neighbour_cell);
5148 }
5149 int value = 0;
5150 while(!neighbour_cells_visited.is_empty())
5151 {
5152 Partition::Cell* const neighbour_cell = neighbour_cells_visited.pop();
5153 if(neighbour_cell->max_ival != neighbour_cell->length)
5154 value++;
5155 neighbour_cell->max_ival = 0;
5156 }
5157 if((value > best_value) or
5158 (value == best_value and cell->length < best_size))
5159 {
5160 best_value = value;
5161 best_size = cell->length;
5162 best_cell = cell;
5163 }
5164 }
5165 return best_cell;
5166 }
5167
5168 /** \internal
5169 * A splitting heuristic.
5170 * Returns the first largest nonsingleton cell with max number of neighbouring
5171 * nonsingleton cells.
5172 * Assumes that the partition p is equitable.
5173 * Assumes that the max_ival fields of the cells are all 0.
5174 */
5175 Partition::Cell*
sh_first_largest_max_neighbours()5176 Graph::sh_first_largest_max_neighbours()
5177 {
5178 Partition::Cell* best_cell = 0;
5179 int best_value = -1;
5180 unsigned int best_size = 0;
5181 KStack<Partition::Cell*> neighbour_cells_visited;
5182 neighbour_cells_visited.init(get_nof_vertices());
5183 for(Partition::Cell* cell = p.first_nonsingleton_cell;
5184 cell;
5185 cell = cell->next_nonsingleton)
5186 {
5187
5188 if(opt_use_comprec and p.cr_get_level(cell->first) != cr_level)
5189 continue;
5190 const Vertex& v = vertices[p.elements[cell->first]];
5191 std::vector<unsigned int>::const_iterator ei = v.edges.begin();
5192 for(unsigned int j = v.nof_edges(); j > 0; j--)
5193 {
5194 Partition::Cell* const neighbour_cell = p.get_cell(*ei++);
5195 if(neighbour_cell->is_unit())
5196 continue;
5197 neighbour_cell->max_ival++;
5198 if(neighbour_cell->max_ival == 1)
5199 neighbour_cells_visited.push(neighbour_cell);
5200 }
5201 int value = 0;
5202 while(!neighbour_cells_visited.is_empty())
5203 {
5204 Partition::Cell* const neighbour_cell = neighbour_cells_visited.pop();
5205 if(neighbour_cell->max_ival != neighbour_cell->length)
5206 value++;
5207 neighbour_cell->max_ival = 0;
5208 }
5209 if((value > best_value) or
5210 (value == best_value and cell->length > best_size))
5211 {
5212 best_value = value;
5213 best_size = cell->length;
5214 best_cell = cell;
5215 }
5216 }
5217 return best_cell;
5218 }
5219
5220
5221
5222
5223
5224
5225
5226
5227
5228
5229
5230
5231
5232
5233
5234
5235
5236
5237
5238
5239 /*-------------------------------------------------------------------------
5240 *
5241 * Initialize the certificate size and memory
5242 *
5243 *-------------------------------------------------------------------------*/
5244
5245 void
initialize_certificate()5246 Graph::initialize_certificate()
5247 {
5248 certificate_index = 0;
5249 certificate_current_path.clear();
5250 certificate_first_path.clear();
5251 certificate_best_path.clear();
5252 }
5253
5254
5255
5256
5257
5258 /*-------------------------------------------------------------------------
5259 *
5260 * Check whether perm is an automorphism.
5261 * Slow, mainly for debugging and validation purposes.
5262 *
5263 *-------------------------------------------------------------------------*/
5264
5265 bool
is_automorphism(unsigned int * const perm)5266 Graph::is_automorphism(unsigned int* const perm)
5267 {
5268 std::set<unsigned int, std::less<unsigned int> > edges1;
5269 std::set<unsigned int, std::less<unsigned int> > edges2;
5270
5271 #if defined(BLISS_CONSISTENCY_CHECKS)
5272 if(!is_permutation(get_nof_vertices(), perm))
5273 _INTERNAL_ERROR();
5274 #endif
5275
5276 for(unsigned int i = 0; i < get_nof_vertices(); i++)
5277 {
5278 Vertex& v1 = vertices[i];
5279 edges1.clear();
5280 for(std::vector<unsigned int>::iterator ei = v1.edges.begin();
5281 ei != v1.edges.end();
5282 ei++)
5283 edges1.insert(perm[*ei]);
5284
5285 Vertex& v2 = vertices[perm[i]];
5286 edges2.clear();
5287 for(std::vector<unsigned int>::iterator ei = v2.edges.begin();
5288 ei != v2.edges.end();
5289 ei++)
5290 edges2.insert(*ei);
5291
5292 if(!(edges1 == edges2))
5293 return false;
5294 }
5295
5296 return true;
5297 }
5298
5299
5300
5301
5302 bool
is_automorphism(const std::vector<unsigned int> & perm) const5303 Graph::is_automorphism(const std::vector<unsigned int>& perm) const
5304 {
5305
5306 if(!(perm.size() == get_nof_vertices() and is_permutation(perm)))
5307 return false;
5308
5309 std::set<unsigned int, std::less<unsigned int> > edges1;
5310 std::set<unsigned int, std::less<unsigned int> > edges2;
5311
5312 for(unsigned int i = 0; i < get_nof_vertices(); i++)
5313 {
5314 const Vertex& v1 = vertices[i];
5315 edges1.clear();
5316 for(std::vector<unsigned int>::const_iterator ei = v1.edges.begin();
5317 ei != v1.edges.end();
5318 ei++)
5319 edges1.insert(perm[*ei]);
5320
5321 const Vertex& v2 = vertices[perm[i]];
5322 edges2.clear();
5323 for(std::vector<unsigned int>::const_iterator ei = v2.edges.begin();
5324 ei != v2.edges.end();
5325 ei++)
5326 edges2.insert(*ei);
5327
5328 if(!(edges1 == edges2))
5329 return false;
5330 }
5331
5332 return true;
5333 }
5334
5335
5336
5337
5338
5339
5340
5341 bool
nucr_find_first_component(const unsigned int level)5342 Graph::nucr_find_first_component(const unsigned int level)
5343 {
5344
5345 cr_component.clear();
5346 cr_component_elements = 0;
5347
5348 /* Find first non-discrete cell in the component level */
5349 Partition::Cell* first_cell = p.first_nonsingleton_cell;
5350 while(first_cell)
5351 {
5352 if(p.cr_get_level(first_cell->first) == level)
5353 break;
5354 first_cell = first_cell->next_nonsingleton;
5355 }
5356
5357 /* The component is discrete, return false */
5358 if(!first_cell)
5359 return false;
5360
5361 std::vector<Partition::Cell*> component;
5362 first_cell->max_ival = 1;
5363 component.push_back(first_cell);
5364
5365 for(unsigned int i = 0; i < component.size(); i++)
5366 {
5367 Partition::Cell* const cell = component[i];
5368
5369 const Vertex& v = vertices[p.elements[cell->first]];
5370 std::vector<unsigned int>::const_iterator ei = v.edges.begin();
5371 for(unsigned int j = v.nof_edges(); j > 0; j--)
5372 {
5373 const unsigned int neighbour = *ei++;
5374
5375 Partition::Cell* const neighbour_cell = p.get_cell(neighbour);
5376
5377 /* Skip unit neighbours */
5378 if(neighbour_cell->is_unit())
5379 continue;
5380 /* Already marked to be in the same component? */
5381 if(neighbour_cell->max_ival == 1)
5382 continue;
5383 /* Is the neighbour at the same component recursion level? */
5384 if(p.cr_get_level(neighbour_cell->first) != level)
5385 continue;
5386
5387 if(neighbour_cell->max_ival_count == 0)
5388 neighbour_heap.insert(neighbour_cell->first);
5389 neighbour_cell->max_ival_count++;
5390 }
5391 while(!neighbour_heap.is_empty())
5392 {
5393 const unsigned int start = neighbour_heap.remove();
5394 Partition::Cell* const neighbour_cell =
5395 p.get_cell(p.elements[start]);
5396
5397 /* Skip saturated neighbour cells */
5398 if(neighbour_cell->max_ival_count == neighbour_cell->length)
5399 {
5400 neighbour_cell->max_ival_count = 0;
5401 continue;
5402 }
5403 neighbour_cell->max_ival_count = 0;
5404 neighbour_cell->max_ival = 1;
5405 component.push_back(neighbour_cell);
5406 }
5407 }
5408
5409 for(unsigned int i = 0; i < component.size(); i++)
5410 {
5411 Partition::Cell* const cell = component[i];
5412 cell->max_ival = 0;
5413 cr_component.push_back(cell->first);
5414 cr_component_elements += cell->length;
5415 }
5416
5417 return true;
5418 }
5419
5420
5421
5422
5423 bool
nucr_find_first_component(const unsigned int level,std::vector<unsigned int> & component,unsigned int & component_elements,Partition::Cell * & sh_return)5424 Graph::nucr_find_first_component(const unsigned int level,
5425 std::vector<unsigned int>& component,
5426 unsigned int& component_elements,
5427 Partition::Cell*& sh_return)
5428 {
5429
5430 component.clear();
5431 component_elements = 0;
5432 sh_return = 0;
5433 unsigned int sh_first = 0;
5434 unsigned int sh_size = 0;
5435 unsigned int sh_nuconn = 0;
5436
5437 /* Find first non-discrete cell in the component level */
5438 Partition::Cell* first_cell = p.first_nonsingleton_cell;
5439 while(first_cell)
5440 {
5441 if(p.cr_get_level(first_cell->first) == level)
5442 break;
5443 first_cell = first_cell->next_nonsingleton;
5444 }
5445
5446 if(!first_cell)
5447 {
5448 /* The component is discrete, return false */
5449 return false;
5450 }
5451
5452 std::vector<Partition::Cell*> comp;
5453 KStack<Partition::Cell*> neighbours;
5454 neighbours.init(get_nof_vertices());
5455
5456 first_cell->max_ival = 1;
5457 comp.push_back(first_cell);
5458
5459 for(unsigned int i = 0; i < comp.size(); i++)
5460 {
5461 Partition::Cell* const cell = comp[i];
5462
5463 const Vertex& v = vertices[p.elements[cell->first]];
5464 std::vector<unsigned int>::const_iterator ei = v.edges.begin();
5465 for(unsigned int j = v.nof_edges(); j > 0; j--)
5466 {
5467 const unsigned int neighbour = *ei++;
5468
5469 Partition::Cell* const neighbour_cell = p.get_cell(neighbour);
5470
5471 /* Skip unit neighbours */
5472 if(neighbour_cell->is_unit())
5473 continue;
5474 /* Is the neighbour at the same component recursion level? */
5475 //if(p.cr_get_level(neighbour_cell->first) != level)
5476 // continue;
5477 if(neighbour_cell->max_ival_count == 0)
5478 neighbours.push(neighbour_cell);
5479 neighbour_cell->max_ival_count++;
5480 }
5481 unsigned int nuconn = 1;
5482 while(!neighbours.is_empty())
5483 {
5484 Partition::Cell* const neighbour_cell = neighbours.pop();
5485 //neighbours.pop_back();
5486
5487 /* Skip saturated neighbour cells */
5488 if(neighbour_cell->max_ival_count == neighbour_cell->length)
5489 {
5490 neighbour_cell->max_ival_count = 0;
5491 continue;
5492 }
5493 nuconn++;
5494 neighbour_cell->max_ival_count = 0;
5495 if(neighbour_cell->max_ival == 0) {
5496 comp.push_back(neighbour_cell);
5497 neighbour_cell->max_ival = 1;
5498 }
5499 }
5500
5501 switch(sh) {
5502 case shs_f:
5503 if(sh_return == 0 or
5504 cell->first <= sh_first) {
5505 sh_return = cell;
5506 sh_first = cell->first;
5507 }
5508 break;
5509 case shs_fs:
5510 if(sh_return == 0 or
5511 cell->length < sh_size or
5512 (cell->length == sh_size and cell->first <= sh_first)) {
5513 sh_return = cell;
5514 sh_first = cell->first;
5515 sh_size = cell->length;
5516 }
5517 break;
5518 case shs_fl:
5519 if(sh_return == 0 or
5520 cell->length > sh_size or
5521 (cell->length == sh_size and cell->first <= sh_first)) {
5522 sh_return = cell;
5523 sh_first = cell->first;
5524 sh_size = cell->length;
5525 }
5526 break;
5527 case shs_fm:
5528 if(sh_return == 0 or
5529 nuconn > sh_nuconn or
5530 (nuconn == sh_nuconn and cell->first <= sh_first)) {
5531 sh_return = cell;
5532 sh_first = cell->first;
5533 sh_nuconn = nuconn;
5534 }
5535 break;
5536 case shs_fsm:
5537 if(sh_return == 0 or
5538 nuconn > sh_nuconn or
5539 (nuconn == sh_nuconn and
5540 (cell->length < sh_size or
5541 (cell->length == sh_size and cell->first <= sh_first)))) {
5542 sh_return = cell;
5543 sh_first = cell->first;
5544 sh_size = cell->length;
5545 sh_nuconn = nuconn;
5546 }
5547 break;
5548 case shs_flm:
5549 if(sh_return == 0 or
5550 nuconn > sh_nuconn or
5551 (nuconn == sh_nuconn and
5552 (cell->length > sh_size or
5553 (cell->length == sh_size and cell->first <= sh_first)))) {
5554 sh_return = cell;
5555 sh_first = cell->first;
5556 sh_size = cell->length;
5557 sh_nuconn = nuconn;
5558 }
5559 break;
5560 default:
5561 fatal_error("Internal error - unknown splitting heuristics");
5562 return 0;
5563 }
5564 }
5565 assert(sh_return);
5566
5567 for(unsigned int i = 0; i < comp.size(); i++)
5568 {
5569 Partition::Cell* const cell = comp[i];
5570 cell->max_ival = 0;
5571 component.push_back(cell->first);
5572 component_elements += cell->length;
5573 }
5574
5575 return true;
5576 }
5577
5578
5579
5580
5581 }
5582