Coverage Report

Created: 2025-10-10 06:10

next uncovered line (L), next uncovered region (R), next uncovered branch (B)
/src/igraph/src/isomorphism/bliss/partition.cc
Line
Count
Source
1
#include <new>
2
#include <cassert>
3
4
#include "graph.hh"
5
#include "partition.hh"
6
7
#include "igraph_decls.h"
8
9
/* Allow using 'and' instead of '&&' with MSVC */
10
#if _MSC_VER
11
#include <ciso646>
12
#endif
13
14
/*
15
  Copyright (c) 2003-2021 Tommi Junttila
16
  Released under the GNU Lesser General Public License version 3.
17
18
  This file is part of bliss.
19
20
  bliss is free software: you can redistribute it and/or modify
21
  it under the terms of the GNU Lesser General Public License as published by
22
  the Free Software Foundation, version 3 of the License.
23
24
  bliss is distributed in the hope that it will be useful,
25
  but WITHOUT ANY WARRANTY; without even the implied warranty of
26
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
27
  GNU Lesser General Public License for more details.
28
29
  You should have received a copy of the GNU Lesser General Public License
30
  along with bliss.  If not, see <http://www.gnu.org/licenses/>.
31
*/
32
33
namespace bliss {
34
35
Partition::Partition()
36
4.32k
{
37
4.32k
  N = 0;
38
4.32k
  elements = 0;
39
4.32k
  in_pos = 0;
40
4.32k
  invariant_values = 0;
41
4.32k
  cells = 0;
42
4.32k
  free_cells = 0;
43
4.32k
  element_to_cell_map = 0;
44
4.32k
  graph = 0;
45
4.32k
  discrete_cell_count = 0;
46
  /* Initialize a distribution count sorting array. */
47
1.11M
  for(unsigned int i = 0; i < 256; i++)
48
1.10M
    dcs_count[i] = 0;
49
50
4.32k
  cr_enabled = false;
51
4.32k
  cr_cells = 0;
52
4.32k
  cr_levels = 0;
53
4.32k
}
54
55
56
57
Partition::~Partition()
58
4.32k
{
59
4.32k
  delete[] elements; elements = nullptr;
60
4.32k
  delete[] cells; cells = nullptr;
61
4.32k
  delete[] element_to_cell_map; element_to_cell_map = nullptr;
62
4.32k
  delete[] in_pos; in_pos = nullptr;
63
4.32k
  delete[] invariant_values; invariant_values = nullptr;
64
4.32k
  N = 0;
65
4.32k
}
66
67
68
69
void Partition::init(const unsigned int M)
70
4.31k
{
71
4.31k
  assert(M > 0);
72
4.31k
  N = M;
73
74
4.31k
  delete[] elements;
75
4.31k
  elements = new unsigned int[N];
76
182k
  for(unsigned int i = 0; i < N; i++)
77
178k
    elements[i] = i;
78
79
4.31k
  delete[] in_pos;
80
4.31k
  in_pos = new unsigned int*[N];
81
182k
  for(unsigned int i = 0; i < N; i++)
82
178k
    in_pos[i] = elements + i;
83
84
4.31k
  delete[] invariant_values;
85
4.31k
  invariant_values = new unsigned int[N];
86
182k
  for(unsigned int i = 0; i < N; i++)
87
178k
    invariant_values[i] = 0;
88
89
4.31k
  delete[] cells;
90
4.31k
  cells = new Cell[N];
91
92
4.31k
  cells[0].first = 0;
93
4.31k
  cells[0].length = N;
94
4.31k
  cells[0].max_ival = 0;
95
4.31k
  cells[0].max_ival_count = 0;
96
4.31k
  cells[0].in_splitting_queue = false;
97
4.31k
  cells[0].in_neighbour_heap = false;
98
4.31k
  cells[0].prev = 0;
99
4.31k
  cells[0].next = 0;
100
4.31k
  cells[0].next_nonsingleton = 0;
101
4.31k
  cells[0].prev_nonsingleton = 0;
102
4.31k
  cells[0].split_level = 0;
103
4.31k
  first_cell = &cells[0];
104
4.31k
  if(N == 1)
105
2
    {
106
2
      first_nonsingleton_cell = 0;
107
2
      discrete_cell_count = 1;
108
2
    }
109
4.31k
  else
110
4.31k
    {
111
4.31k
      first_nonsingleton_cell = &cells[0];
112
4.31k
      discrete_cell_count = 0;
113
4.31k
    }
114
115
178k
  for(unsigned int i = 1; i < N; i++)
116
174k
    {
117
174k
      cells[i].first = 0;
118
174k
      cells[i].length = 0;
119
174k
      cells[i].max_ival = 0;
120
174k
      cells[i].max_ival_count = 0;
121
174k
      cells[i].in_splitting_queue = false;
122
174k
      cells[i].in_neighbour_heap = false;
123
174k
      cells[i].prev = 0;
124
174k
      cells[i].next = (i < N-1)?&cells[i+1]:0;
125
174k
      cells[i].next_nonsingleton = 0;
126
174k
      cells[i].prev_nonsingleton = 0;
127
174k
    }
128
4.31k
  if(N > 1)
129
4.31k
    free_cells = &cells[1];
130
2
  else
131
2
    free_cells = 0;
132
133
4.31k
  delete[] element_to_cell_map;
134
4.31k
  element_to_cell_map = new Cell*[N];
135
182k
  for(unsigned int i = 0; i < N; i++)
136
178k
    element_to_cell_map[i] = first_cell;
137
138
4.31k
  splitting_queue.init(N);
139
4.31k
  refinement_stack.init(N);
140
141
  /* Reset the main backtracking stack */
142
4.31k
  bt_stack.clear();
143
4.31k
}
144
145
146
147
148
149
150
Partition::BacktrackPoint
151
Partition::set_backtrack_point()
152
4.63M
{
153
4.63M
  BacktrackInfo info;
154
4.63M
  info.refinement_stack_size = refinement_stack.size();
155
4.63M
  if(cr_enabled)
156
4.63M
    info.cr_backtrack_point = cr_get_backtrack_point();
157
4.63M
  BacktrackPoint p = bt_stack.size();
158
4.63M
  bt_stack.push_back(info);
159
4.63M
  return p;
160
4.63M
}
161
162
163
164
void
165
Partition::goto_backtrack_point(BacktrackPoint p)
166
2.44M
{
167
2.44M
  assert(p < bt_stack.size());
168
2.44M
  BacktrackInfo info = bt_stack[p];
169
2.44M
  bt_stack.resize(p);
170
171
2.44M
  if(cr_enabled)
172
2.44M
    cr_goto_backtrack_point(info.cr_backtrack_point);
173
174
2.44M
  const unsigned int dest_refinement_stack_size = info.refinement_stack_size;
175
176
2.44M
  assert(refinement_stack.size() >= dest_refinement_stack_size);
177
4.93M
  while(refinement_stack.size() > dest_refinement_stack_size)
178
2.49M
    {
179
2.49M
      RefInfo i = refinement_stack.pop();
180
2.49M
      const unsigned int first = i.split_cell_first;
181
2.49M
      Cell* cell = get_cell(elements[first]);
182
183
2.49M
      if(cell->first != first)
184
2.20M
        {
185
2.20M
          assert(cell->first < first);
186
2.20M
          assert(cell->split_level <= dest_refinement_stack_size);
187
2.20M
          goto done;
188
2.20M
        }
189
2.49M
      assert(cell->split_level > dest_refinement_stack_size);
190
191
621k
      while(cell->split_level > dest_refinement_stack_size)
192
336k
        {
193
336k
          assert(cell->prev);
194
336k
          cell = cell->prev;
195
336k
        }
196
2.78M
      while(cell->next and
197
2.76M
            cell->next->split_level > dest_refinement_stack_size)
198
2.49M
        {
199
          /* Merge next cell */
200
2.49M
          Cell* const next_cell = cell->next;
201
2.49M
          if(cell->length == 1)
202
158k
            discrete_cell_count--;
203
2.49M
          if(next_cell->length == 1)
204
2.47M
            discrete_cell_count--;
205
          /* Update element_to_cell_map values of elements added in cell */
206
2.49M
          unsigned int* ep = elements + next_cell->first;
207
2.49M
          unsigned int* const lp = ep + next_cell->length;
208
5.02M
          for( ; ep < lp; ep++)
209
2.53M
            element_to_cell_map[*ep] = cell;
210
          /* Update cell parameters */
211
2.49M
          cell->length += next_cell->length;
212
2.49M
          if(next_cell->next)
213
2.48M
            next_cell->next->prev = cell;
214
2.49M
          cell->next = next_cell->next;
215
          /* (Pseudo)free next_cell */
216
2.49M
          next_cell->first = 0;
217
2.49M
          next_cell->length = 0;
218
2.49M
          next_cell->prev = 0;
219
2.49M
          next_cell->next = free_cells;
220
2.49M
          free_cells = next_cell;
221
2.49M
        }
222
223
2.49M
    done:
224
2.49M
      if(i.prev_nonsingleton_first >= 0)
225
66.4k
        {
226
66.4k
          Cell* const prev_cell = get_cell(elements[i.prev_nonsingleton_first]);
227
66.4k
          assert(prev_cell->length > 1);
228
66.4k
          cell->prev_nonsingleton = prev_cell;
229
66.4k
          prev_cell->next_nonsingleton = cell;
230
66.4k
        }
231
2.42M
      else
232
2.42M
        {
233
          //assert(cell->prev_nonsingleton == 0);
234
2.42M
          cell->prev_nonsingleton = 0;
235
2.42M
          first_nonsingleton_cell = cell;
236
2.42M
        }
237
238
2.49M
      if(i.next_nonsingleton_first >= 0)
239
2.04M
        {
240
2.04M
          Cell* const next_cell = get_cell(elements[i.next_nonsingleton_first]);
241
2.04M
          assert(next_cell->length > 1);
242
2.04M
          cell->next_nonsingleton = next_cell;
243
2.04M
          next_cell->prev_nonsingleton = cell;
244
2.04M
        }
245
452k
      else
246
452k
        {
247
          //assert(cell->next_nonsingleton == 0);
248
452k
          cell->next_nonsingleton = 0;
249
452k
        }
250
2.49M
    }
251
252
2.44M
}
253
254
255
256
Partition::Cell*
257
Partition::individualize(Partition::Cell * const cell,
258
                         const unsigned int element)
259
2.32M
{
260
2.32M
  assert(!cell->is_unit());
261
262
2.32M
  unsigned int * const pos = in_pos[element];
263
2.32M
  assert((unsigned int)(pos - elements) >= cell->first);
264
2.32M
  assert((unsigned int)(pos - elements) < cell->first + cell->length);
265
2.32M
  assert(*pos == element);
266
267
2.32M
  const unsigned int last = cell->first + cell->length - 1;
268
2.32M
  *pos = elements[last];
269
2.32M
  in_pos[*pos] = pos;
270
2.32M
  elements[last] = element;
271
2.32M
  in_pos[element] = elements + last;
272
273
2.32M
  Partition::Cell * const new_cell = aux_split_in_two(cell, cell->length-1);
274
2.32M
  assert(elements[new_cell->first] == element);
275
2.32M
  element_to_cell_map[element] = new_cell;
276
277
2.32M
  return new_cell;
278
2.32M
}
279
280
281
282
Partition::Cell*
283
Partition::aux_split_in_two(Partition::Cell* const cell,
284
                            const unsigned int first_half_size)
285
2.49M
{
286
2.49M
  RefInfo i;
287
288
2.49M
  assert(0 < first_half_size && first_half_size < cell->length);
289
290
  /* (Pseudo)allocate new cell */
291
2.49M
  Cell * const new_cell = free_cells;
292
2.49M
  assert(new_cell != 0);
293
2.49M
  free_cells = new_cell->next;
294
  /* Update new cell parameters */
295
2.49M
  new_cell->first = cell->first + first_half_size;
296
2.49M
  new_cell->length = cell->length - first_half_size;
297
2.49M
  new_cell->next = cell->next;
298
2.49M
  if(new_cell->next)
299
2.48M
    new_cell->next->prev = new_cell;
300
2.49M
  new_cell->prev = cell;
301
2.49M
  new_cell->split_level = refinement_stack.size()+1;
302
  /* Update old, splitted cell parameters */
303
2.49M
  cell->length = first_half_size;
304
2.49M
  cell->next = new_cell;
305
  /* CR */
306
2.49M
  if(cr_enabled)
307
2.49M
    cr_create_at_level_trailed(new_cell->first, cr_get_level(cell->first));
308
309
  /* Add cell in refinement_stack for backtracking */
310
2.49M
  i.split_cell_first = new_cell->first;
311
2.49M
  if(cell->prev_nonsingleton)
312
75.3k
    i.prev_nonsingleton_first = cell->prev_nonsingleton->first;
313
2.42M
  else
314
2.42M
    i.prev_nonsingleton_first = -1;
315
2.49M
  if(cell->next_nonsingleton)
316
2.03M
    i.next_nonsingleton_first = cell->next_nonsingleton->first;
317
461k
  else
318
461k
    i.next_nonsingleton_first = -1;
319
2.49M
  refinement_stack.push(i);
320
321
  /* Modify nonsingleton cell list */
322
2.49M
  if(new_cell->length > 1)
323
38.6k
    {
324
38.6k
      new_cell->prev_nonsingleton = cell;
325
38.6k
      new_cell->next_nonsingleton = cell->next_nonsingleton;
326
38.6k
      if(new_cell->next_nonsingleton)
327
27.0k
        new_cell->next_nonsingleton->prev_nonsingleton = new_cell;
328
38.6k
      cell->next_nonsingleton = new_cell;
329
38.6k
    }
330
2.45M
  else
331
2.45M
    {
332
2.45M
      new_cell->next_nonsingleton = 0;
333
2.45M
      new_cell->prev_nonsingleton = 0;
334
2.45M
      discrete_cell_count++;
335
2.45M
    }
336
337
2.49M
  if(cell->is_unit())
338
183k
    {
339
183k
      if(cell->prev_nonsingleton)
340
20.8k
        cell->prev_nonsingleton->next_nonsingleton = cell->next_nonsingleton;
341
162k
      else
342
162k
        first_nonsingleton_cell = cell->next_nonsingleton;
343
183k
      if(cell->next_nonsingleton)
344
143k
        cell->next_nonsingleton->prev_nonsingleton = cell->prev_nonsingleton;
345
183k
      cell->next_nonsingleton = 0;
346
183k
      cell->prev_nonsingleton = 0;
347
183k
      discrete_cell_count++;
348
183k
    }
349
350
2.49M
  return new_cell;
351
2.49M
}
352
353
354
355
void
356
Partition::splitting_queue_add(Cell* const cell)
357
2.71M
{
358
2.71M
  static const unsigned int smallish_cell_threshold = 1;
359
2.71M
  assert(!cell->in_splitting_queue);
360
2.71M
  cell->in_splitting_queue = true;
361
2.71M
  if(cell->length <= smallish_cell_threshold)
362
2.64M
    splitting_queue.push_front(cell);
363
66.2k
  else
364
66.2k
    splitting_queue.push_back(cell);
365
2.71M
}
366
367
368
369
void
370
Partition::splitting_queue_clear()
371
24.4k
{
372
53.8k
  while(!splitting_queue_is_empty())
373
29.3k
    splitting_queue_pop();
374
24.4k
}
375
376
377
378
379
380
/*
381
 * Assumes that the invariant values are NOT the same
382
 * and that the cell contains more than one element
383
 */
384
Partition::Cell*
385
Partition::sort_and_split_cell1(Partition::Cell* const cell)
386
21.1k
{
387
#if defined(BLISS_EXPENSIVE_CONSISTENCY_CHECKS)
388
  assert(cell->length > 1);
389
  assert(cell->first + cell->length <= N);
390
  unsigned int nof_0_found = 0;
391
  unsigned int nof_1_found = 0;
392
  for(unsigned int i = cell->first; i < cell->first + cell->length; i++)
393
    {
394
      const unsigned int ival = invariant_values[elements[i]];
395
      assert(ival == 0 or ival == 1);
396
      if(ival == 0) nof_0_found++;
397
      else nof_1_found++;
398
    }
399
  assert(nof_0_found > 0);
400
  assert(nof_1_found > 0);
401
  assert(nof_1_found == cell->max_ival_count);
402
  assert(nof_0_found + nof_1_found == cell->length);
403
  assert(cell->max_ival == 1);
404
#endif
405
406
407
  /* (Pseudo)allocate new cell */
408
21.1k
  Cell* const new_cell = free_cells;
409
21.1k
  assert(new_cell != 0);
410
21.1k
  free_cells = new_cell->next;
411
412
21.1k
#define NEW_SORT1
413
21.1k
#ifdef NEW_SORT1
414
21.1k
      unsigned int *ep0 = elements + cell->first;
415
21.1k
      unsigned int *ep1 = ep0 + cell->length - cell->max_ival_count;
416
21.1k
      if(cell->max_ival_count > cell->length / 2)
417
2.40k
        {
418
          /* There are more ones than zeros, only move zeros */
419
2.40k
          unsigned int * const end = ep0 + cell->length;
420
24.1k
          while(ep1 < end)
421
21.7k
            {
422
28.4k
              while(invariant_values[*ep1] == 0)
423
6.64k
                {
424
6.64k
                  const unsigned int tmp = *ep1;
425
6.64k
                  *ep1 = *ep0;
426
6.64k
                  *ep0 = tmp;
427
6.64k
                  in_pos[tmp] = ep0;
428
6.64k
                  in_pos[*ep1] = ep1;
429
6.64k
                  ep0++;
430
6.64k
                }
431
21.7k
              element_to_cell_map[*ep1] = new_cell;
432
21.7k
              invariant_values[*ep1] = 0;
433
21.7k
              ep1++;
434
21.7k
            }
435
2.40k
        }
436
18.7k
      else
437
18.7k
        {
438
          /* There are more zeros than ones, only move ones */
439
18.7k
          unsigned int * const end = ep1;
440
321k
          while(ep0 < end)
441
302k
            {
442
346k
              while(invariant_values[*ep0] != 0)
443
43.7k
                {
444
43.7k
                  const unsigned int tmp = *ep0;
445
43.7k
                  *ep0 = *ep1;
446
43.7k
                  *ep1 = tmp;
447
43.7k
                  in_pos[tmp] = ep1;
448
43.7k
                  in_pos[*ep0] = ep0;
449
43.7k
                  ep1++;
450
43.7k
                }
451
302k
              ep0++;
452
302k
            }
453
18.7k
          ep1 = end;
454
72.8k
          while(ep1 < elements + cell->first + cell->length)
455
54.0k
            {
456
54.0k
              element_to_cell_map[*ep1] = new_cell;
457
54.0k
              invariant_values[*ep1] = 0;
458
54.0k
              ep1++;
459
54.0k
            }
460
18.7k
        }
461
  /* Update new cell parameters */
462
21.1k
  new_cell->first = cell->first + cell->length - cell->max_ival_count;
463
21.1k
  new_cell->length = cell->length - (new_cell->first - cell->first);
464
21.1k
  new_cell->next = cell->next;
465
21.1k
  if(new_cell->next)
466
18.8k
    new_cell->next->prev = new_cell;
467
21.1k
  new_cell->prev = cell;
468
21.1k
  new_cell->split_level = refinement_stack.size()+1;
469
  /* Update old, splitted cell parameters */
470
21.1k
  cell->length = new_cell->first - cell->first;
471
21.1k
  cell->next = new_cell;
472
  /* CR */
473
21.1k
  if(cr_enabled)
474
21.1k
    cr_create_at_level_trailed(new_cell->first, cr_get_level(cell->first));
475
476
#else
477
  /* Sort vertices in the cell according to the invariant values */
478
  unsigned int *ep0 = elements + cell->first;
479
  unsigned int *ep1 = ep0 + cell->length;
480
  while(ep1 > ep0)
481
    {
482
      const unsigned int element = *ep0;
483
      const unsigned int ival = invariant_values[element];
484
      invariant_values[element] = 0;
485
      assert(ival <= 1);
486
      assert(element_to_cell_map[element] == cell);
487
      assert(in_pos[element] == ep0);
488
      if(ival == 0)
489
        {
490
          ep0++;
491
        }
492
      else
493
        {
494
          ep1--;
495
          *ep0 = *ep1;
496
          *ep1 = element;
497
          element_to_cell_map[element] = new_cell;
498
          in_pos[element] = ep1;
499
          in_pos[*ep0] = ep0;
500
        }
501
    }
502
503
  assert(ep1 != elements + cell->first);
504
  assert(ep0 != elements + cell->first + cell->length);
505
506
  /* Update new cell parameters */
507
  new_cell->first = ep1 - elements;
508
  new_cell->length = cell->length - (new_cell->first - cell->first);
509
  new_cell->next = cell->next;
510
  if(new_cell->next)
511
    new_cell->next->prev = new_cell;
512
  new_cell->prev = cell;
513
  new_cell->split_level = cell->split_level;
514
  /* Update old, splitted cell parameters */
515
  cell->length = new_cell->first - cell->first;
516
  cell->next = new_cell;
517
  cell->split_level = refinement_stack.size()+1;
518
  /* CR */
519
  if(cr_enabled)
520
    cr_create_at_level_trailed(new_cell->first, cr_get_level(cell->first));
521
522
#endif /* ifdef NEW_SORT1*/
523
524
  /* Add cell in refinement stack for backtracking */
525
21.1k
  {
526
21.1k
    RefInfo i;
527
21.1k
    i.split_cell_first = new_cell->first;
528
21.1k
    if(cell->prev_nonsingleton)
529
7.05k
      i.prev_nonsingleton_first = cell->prev_nonsingleton->first;
530
14.1k
    else
531
14.1k
      i.prev_nonsingleton_first = -1;
532
21.1k
    if(cell->next_nonsingleton)
533
17.8k
      i.next_nonsingleton_first = cell->next_nonsingleton->first;
534
3.34k
    else
535
3.34k
      i.next_nonsingleton_first = -1;
536
    /* Modify nonsingleton cell list */
537
21.1k
    if(new_cell->length > 1)
538
19.0k
      {
539
19.0k
        new_cell->prev_nonsingleton = cell;
540
19.0k
        new_cell->next_nonsingleton = cell->next_nonsingleton;
541
19.0k
        if(new_cell->next_nonsingleton)
542
16.5k
          new_cell->next_nonsingleton->prev_nonsingleton = new_cell;
543
19.0k
        cell->next_nonsingleton = new_cell;
544
19.0k
      }
545
2.09k
    else
546
2.09k
      {
547
2.09k
        new_cell->next_nonsingleton = 0;
548
2.09k
        new_cell->prev_nonsingleton = 0;
549
2.09k
        discrete_cell_count++;
550
2.09k
      }
551
21.1k
    if(cell->is_unit())
552
1.95k
      {
553
1.95k
        if(cell->prev_nonsingleton)
554
1.22k
          cell->prev_nonsingleton->next_nonsingleton = cell->next_nonsingleton;
555
731
        else
556
731
          first_nonsingleton_cell = cell->next_nonsingleton;
557
1.95k
        if(cell->next_nonsingleton)
558
1.69k
          cell->next_nonsingleton->prev_nonsingleton = cell->prev_nonsingleton;
559
1.95k
        cell->next_nonsingleton = 0;
560
1.95k
        cell->prev_nonsingleton = 0;
561
1.95k
        discrete_cell_count++;
562
1.95k
      }
563
21.1k
    refinement_stack.push(i);
564
21.1k
  }
565
566
567
  /* Add cells in splitting queue */
568
21.1k
  assert(!new_cell->in_splitting_queue);
569
21.1k
  if(cell->in_splitting_queue) {
570
    /* Both cells must be included in splitting_queue in order to have
571
       refinement to equitable partition */
572
2.40k
    splitting_queue_add(new_cell);
573
18.7k
  } else {
574
18.7k
    Cell *min_cell, *max_cell;
575
18.7k
    if(cell->length <= new_cell->length) {
576
4.15k
      min_cell = cell;
577
4.15k
      max_cell = new_cell;
578
14.6k
    } else {
579
14.6k
      min_cell = new_cell;
580
14.6k
      max_cell = cell;
581
14.6k
    }
582
    /* Put the smaller cell in splitting_queue */
583
18.7k
    splitting_queue_add(min_cell);
584
18.7k
    if(max_cell->is_unit()) {
585
      /* Put the "larger" cell also in splitting_queue */
586
335
      splitting_queue_add(max_cell);
587
335
    }
588
18.7k
  }
589
590
591
21.1k
  return new_cell;
592
21.1k
}
593
594
595
596
597
598
/**
599
 * An auxiliary function for distribution count sorting.
600
 * Build start array so that
601
 * dcs_start[0] = 0 and dcs_start[i+1] = dcs_start[i] + dcs_count[i].
602
 */
603
void
604
Partition::dcs_cumulate_count(const unsigned int max)
605
12.6k
{
606
12.6k
  assert(max <= 255);
607
12.6k
  unsigned int* count_p = dcs_count;
608
12.6k
  unsigned int* start_p = dcs_start;
609
12.6k
  unsigned int sum = 0;
610
57.2k
  for(unsigned int i = max+1; i > 0; i--)
611
44.5k
    {
612
44.5k
      *start_p = sum;
613
44.5k
      start_p++;
614
44.5k
      sum += *count_p;
615
44.5k
      count_p++;
616
44.5k
    }
617
12.6k
}
618
619
620
/**
621
 * Distribution count sorting of cells with invariant values less than 256.
622
 */
623
Partition::Cell*
624
Partition::sort_and_split_cell255(Partition::Cell* const cell,
625
                                  const unsigned int max_ival)
626
12.6k
{
627
12.6k
  assert(max_ival <= 255);
628
629
12.6k
  if(cell->is_unit())
630
0
    {
631
      /* Reset invariant value */
632
0
      invariant_values[elements[cell->first]] = 0;
633
0
      return cell;
634
0
    }
635
636
#ifdef BLISS_CONSISTENCY_CHECKS
637
  for(unsigned int i = 0; i < 256; i++)
638
    assert(dcs_count[i] == 0);
639
#endif
640
641
  /*
642
   * Compute the distribution of invariant values to the count array
643
   */
644
12.6k
  {
645
12.6k
    const unsigned int *ep = elements + cell->first;
646
12.6k
    assert(element_to_cell_map[*ep] == cell);
647
12.6k
    const unsigned int ival = invariant_values[*ep];
648
12.6k
    assert(ival <= 255);
649
12.6k
    dcs_count[ival]++;
650
12.6k
    ep++;
651
#if defined(BLISS_CONSISTENCY_CHECKS)
652
    bool equal_invariant_values = true;
653
#endif
654
270k
    for(unsigned int i = cell->length - 1; i != 0; i--)
655
258k
      {
656
258k
        assert(element_to_cell_map[*ep] == cell);
657
258k
        const unsigned int ival2 = invariant_values[*ep];
658
258k
        assert(ival2 <= 255);
659
258k
        assert(ival2 <= max_ival);
660
258k
        dcs_count[ival2]++;
661
#if defined(BLISS_CONSISTENCY_CHECKS)
662
        if(ival2 != ival) {
663
          equal_invariant_values = false;
664
        }
665
#endif
666
258k
        ep++;
667
258k
      }
668
#if defined(BLISS_CONSISTENCY_CHECKS)
669
    assert(!equal_invariant_values);
670
    if(equal_invariant_values) {
671
      assert(dcs_count[ival] == cell->length);
672
      dcs_count[ival] = 0;
673
      clear_ivs(cell);
674
      return cell;
675
    }
676
#endif
677
12.6k
  }
678
679
  /* Build start array */
680
12.6k
  dcs_cumulate_count(max_ival);
681
682
  //assert(dcs_start[255] + dcs_count[255] == cell->length);
683
12.6k
  assert(dcs_start[max_ival] + dcs_count[max_ival] == cell->length);
684
685
  /* Do the sorting */
686
57.2k
  for(unsigned int i = 0; i <= max_ival; i++)
687
44.5k
    {
688
44.5k
      unsigned int *ep = elements + cell->first + dcs_start[i];
689
245k
      for(unsigned int j = dcs_count[i]; j > 0; j--)
690
201k
        {
691
270k
          while(true)
692
270k
            {
693
270k
              const unsigned int element = *ep;
694
270k
              const unsigned int ival = invariant_values[element];
695
270k
              if(ival == i)
696
201k
                break;
697
270k
              assert(ival > i);
698
69.6k
              assert(dcs_count[ival] > 0);
699
69.6k
              *ep = elements[cell->first + dcs_start[ival]];
700
69.6k
              elements[cell->first + dcs_start[ival]] = element;
701
69.6k
              dcs_start[ival]++;
702
69.6k
              dcs_count[ival]--;
703
69.6k
            }
704
201k
          ep++;
705
201k
        }
706
44.5k
      dcs_count[i] = 0;
707
44.5k
    }
708
709
#if defined(BLISS_CONSISTENCY_CHECKS)
710
  for(unsigned int i = 0; i < 256; i++)
711
    assert(dcs_count[i] == 0);
712
#endif
713
714
  /* split cell */
715
12.6k
  Cell* const new_cell = split_cell(cell);
716
12.6k
  assert(new_cell != cell);
717
12.6k
  return new_cell;
718
12.6k
}
719
720
721
722
/*
723
 * Sort the elements in a cell according to their invariant values.
724
 * The invariant values are not cleared.
725
 * Warning: the in_pos array is left in incorrect state.
726
 */
727
bool
728
Partition::shellsort_cell(Partition::Cell* const cell)
729
0
{
730
0
  unsigned int h;
731
0
  unsigned int* ep;
732
733
  //assert(cell->first + cell->length <= N);
734
735
0
  if(cell->is_unit())
736
0
    return false;
737
738
  /* Check whether all the elements have the same invariant value */
739
0
  bool equal_invariant_values = true;
740
0
  {
741
0
    ep = elements + cell->first;
742
0
    const unsigned int ival = invariant_values[*ep];
743
0
    assert(element_to_cell_map[*ep] == cell);
744
0
    ep++;
745
0
    for(unsigned int i = cell->length - 1; i > 0; i--)
746
0
      {
747
0
        assert(element_to_cell_map[*ep] == cell);
748
0
        if(invariant_values[*ep] != ival) {
749
0
          equal_invariant_values = false;
750
0
          break;
751
0
        }
752
0
        ep++;
753
0
      }
754
0
  }
755
0
  if(equal_invariant_values)
756
0
    return false;
757
758
0
  ep = elements + cell->first;
759
760
0
  for(h = 1; h <= cell->length/9; h = 3*h + 1)
761
0
    ;
762
0
  for( ; h > 0; h = h/3) {
763
0
    for(unsigned int i = h; i < cell->length; i++) {
764
0
      const unsigned int element = ep[i];
765
0
      const unsigned int ival = invariant_values[element];
766
0
      unsigned int j = i;
767
0
      while(j >= h and invariant_values[ep[j-h]] > ival) {
768
0
        ep[j] = ep[j-h];
769
0
        j -= h;
770
0
      }
771
0
      ep[j] = element;
772
0
    }
773
0
  }
774
0
  return true;
775
0
}
776
777
778
779
void
780
Partition::clear_ivs(Cell* const cell)
781
33.9k
{
782
33.9k
  unsigned int* ep = elements + cell->first;
783
156k
  for(unsigned int i = cell->length; i > 0; i--, ep++)
784
123k
    invariant_values[*ep] = 0;
785
33.9k
}
786
787
788
/*
789
 * Assumes that the elements in the cell are sorted according to their
790
 * invariant values.
791
 */
792
Partition::Cell*
793
Partition::split_cell(Partition::Cell* const original_cell)
794
12.6k
{
795
12.6k
  Cell* cell = original_cell;
796
12.6k
  const bool original_cell_was_in_splitting_queue =
797
12.6k
    original_cell->in_splitting_queue;
798
12.6k
  Cell* largest_new_cell = 0;
799
800
32.2k
  while(true)
801
32.2k
    {
802
32.2k
      unsigned int* ep = elements + cell->first;
803
32.2k
      const unsigned int* const lp = ep + cell->length;
804
32.2k
      const unsigned int ival = invariant_values[*ep];
805
32.2k
      invariant_values[*ep] = 0;
806
32.2k
      element_to_cell_map[*ep] = cell;
807
32.2k
      in_pos[*ep] = ep;
808
32.2k
      ep++;
809
270k
      while(ep < lp)
810
258k
        {
811
258k
          const unsigned int e = *ep;
812
258k
          if(invariant_values[e] != ival)
813
19.5k
            break;
814
238k
          invariant_values[e] = 0;
815
238k
          in_pos[e] = ep;
816
238k
          ep++;
817
238k
          element_to_cell_map[e] = cell;
818
238k
        }
819
32.2k
      if(ep == lp)
820
12.6k
        break;
821
822
19.5k
      Cell* const new_cell = aux_split_in_two(cell,
823
19.5k
                                              (ep - elements) - cell->first);
824
825
19.5k
      if(graph and graph->compute_eqref_hash)
826
0
        {
827
0
          graph->eqref_hash.update(new_cell->first);
828
0
          graph->eqref_hash.update(new_cell->length);
829
0
          graph->eqref_hash.update(ival);
830
0
        }
831
832
      /* Add cells in splitting_queue */
833
19.5k
      assert(!new_cell->is_in_splitting_queue());
834
19.5k
      if(original_cell_was_in_splitting_queue)
835
616
        {
836
          /* In this case, all new cells are inserted in splitting_queue */
837
616
          assert(cell->is_in_splitting_queue());
838
616
          splitting_queue_add(new_cell);
839
616
        }
840
18.9k
      else
841
18.9k
        {
842
          /* Otherwise, we can omit one new cell from splitting_queue */
843
18.9k
          assert(!cell->is_in_splitting_queue());
844
18.9k
          if(largest_new_cell == 0) {
845
12.2k
            largest_new_cell = cell;
846
12.2k
          } else {
847
6.70k
            assert(!largest_new_cell->is_in_splitting_queue());
848
6.70k
            if(cell->length > largest_new_cell->length) {
849
2.22k
              splitting_queue_add(largest_new_cell);
850
2.22k
              largest_new_cell = cell;
851
4.47k
            } else {
852
4.47k
              splitting_queue_add(cell);
853
4.47k
            }
854
6.70k
          }
855
18.9k
        }
856
      /* Process the rest of the cell */
857
19.5k
      cell = new_cell;
858
19.5k
    }
859
860
861
12.6k
  if(original_cell == cell) {
862
    /* All the elements in cell had the same invariant value */
863
0
    return cell;
864
0
  }
865
866
  /* Add cells in splitting_queue */
867
12.6k
  if(!original_cell_was_in_splitting_queue)
868
12.2k
    {
869
      /* Also consider the last new cell */
870
12.2k
      assert(largest_new_cell);
871
12.2k
      if(cell->length > largest_new_cell->length)
872
846
        {
873
846
          splitting_queue_add(largest_new_cell);
874
846
          largest_new_cell = cell;
875
846
        }
876
11.3k
      else
877
11.3k
        {
878
11.3k
          splitting_queue_add(cell);
879
11.3k
        }
880
12.2k
      if(largest_new_cell->is_unit())
881
545
        {
882
          /* Needed in certificate computation */
883
545
          splitting_queue_add(largest_new_cell);
884
545
        }
885
12.2k
    }
886
887
12.6k
  return cell;
888
12.6k
}
889
890
891
Partition::Cell*
892
Partition::zplit_cell(Partition::Cell* const cell,
893
                      const bool max_ival_info_ok)
894
74.2k
{
895
74.2k
  assert(cell != 0);
896
897
74.2k
  Cell* last_new_cell = cell;
898
899
74.2k
  if(!max_ival_info_ok)
900
0
    {
901
      /* Compute max_ival info */
902
0
      assert(cell->max_ival == 0);
903
0
      assert(cell->max_ival_count == 0);
904
0
      unsigned int *ep = elements + cell->first;
905
0
      for(unsigned int i = cell->length; i > 0; i--, ep++)
906
0
        {
907
0
          const unsigned int ival = invariant_values[*ep];
908
0
          if(ival > cell->max_ival)
909
0
            {
910
0
              cell->max_ival = ival;
911
0
              cell->max_ival_count = 1;
912
0
            }
913
0
          else if(ival == cell->max_ival)
914
0
            {
915
0
              cell->max_ival_count++;
916
0
            }
917
0
        }
918
0
    }
919
920
#ifdef BLISS_CONSISTENCY_CHECKS
921
  /* Verify max_ival info */
922
  {
923
    unsigned int nof_zeros = 0;
924
    unsigned int max_ival = 0;
925
    unsigned int max_ival_count = 0;
926
    unsigned int *ep = elements + cell->first;
927
    for(unsigned int i = cell->length; i > 0; i--, ep++)
928
      {
929
        const unsigned int ival = invariant_values[*ep];
930
        if(ival == 0)
931
          nof_zeros++;
932
        if(ival > max_ival)
933
          {
934
            max_ival = ival;
935
            max_ival_count = 1;
936
          }
937
        else if(ival == max_ival)
938
          max_ival_count++;
939
      }
940
    assert(max_ival == cell->max_ival);
941
    assert(max_ival_count == cell->max_ival_count);
942
  }
943
#endif
944
945
  /* max_ival info has been computed */
946
947
74.2k
  if(cell->max_ival_count == cell->length)
948
40.3k
    {
949
      /* All invariant values are the same, clear 'em */
950
40.3k
      if(cell->max_ival > 0)
951
30.9k
        clear_ivs(cell);
952
40.3k
    }
953
33.8k
  else
954
33.8k
    {
955
      /* All invariant values are not the same */
956
33.8k
      if(cell->max_ival == 1)
957
21.1k
        {
958
          /* Specialized splitting for cells with binary invariant values */
959
21.1k
          last_new_cell = sort_and_split_cell1(cell);
960
21.1k
        }
961
12.6k
      else if(cell->max_ival < 256)
962
12.6k
        {
963
          /* Specialized splitting for cells with invariant values < 256 */
964
12.6k
          last_new_cell = sort_and_split_cell255(cell, cell->max_ival);
965
12.6k
        }
966
0
      else
967
0
        {
968
          /* Generic sorting and splitting */
969
0
          const bool sorted = shellsort_cell(cell);
970
0
          assert(sorted);
971
0
          IGRAPH_UNUSED(sorted);
972
0
          last_new_cell = split_cell(cell);
973
0
        }
974
33.8k
    }
975
74.2k
  cell->max_ival = 0;
976
74.2k
  cell->max_ival_count = 0;
977
74.2k
  return last_new_cell;
978
74.2k
}
979
980
981
982
/*
983
 *
984
 * Component recursion specific code
985
 *
986
 */
987
void
988
Partition::cr_init()
989
4.31k
{
990
4.31k
  assert(bt_stack.empty());
991
992
4.31k
  cr_enabled = true;
993
994
4.31k
  delete[] cr_cells;
995
4.31k
  cr_cells = new CRCell[N];
996
997
4.31k
  delete[] cr_levels;
998
4.31k
  cr_levels = new CRCell*[N];
999
1000
182k
  for(unsigned int i = 0; i < N; i++) {
1001
178k
    cr_levels[i] = 0;
1002
178k
    cr_cells[i].level = UINT_MAX;
1003
178k
    cr_cells[i].next = 0;
1004
178k
    cr_cells[i].prev_next_ptr = 0;
1005
178k
  }
1006
1007
8.62k
  for(const Cell *cell = first_cell; cell; cell = cell->next)
1008
4.31k
    cr_create_at_level_trailed(cell->first, 0);
1009
1010
4.31k
  cr_max_level = 0;
1011
4.31k
}
1012
1013
1014
void
1015
Partition::cr_free()
1016
4.31k
{
1017
4.31k
  delete[] cr_cells; cr_cells = nullptr;
1018
4.31k
  delete[] cr_levels; cr_levels = nullptr;
1019
1020
4.31k
  cr_created_trail.clear();
1021
4.31k
  cr_splitted_level_trail.clear();
1022
4.31k
  cr_bt_info.clear();
1023
4.31k
  cr_max_level = 0;
1024
1025
4.31k
  cr_enabled = false;
1026
4.31k
}
1027
1028
1029
unsigned int
1030
Partition::cr_split_level(const unsigned int level,
1031
                          const std::vector<unsigned int>& splitted_cells)
1032
18.5k
{
1033
18.5k
  assert(cr_enabled);
1034
18.5k
  assert(level <= cr_max_level);
1035
18.5k
  cr_levels[++cr_max_level] = 0;
1036
18.5k
  cr_splitted_level_trail.push_back(level);
1037
1038
45.8k
  for(unsigned int i = 0; i < splitted_cells.size(); i++)
1039
27.2k
    {
1040
27.2k
      const unsigned int cell_index = splitted_cells[i];
1041
27.2k
      assert(cell_index < N);
1042
27.2k
      CRCell& cr_cell = cr_cells[cell_index];
1043
27.2k
      assert(cr_cell.level == level);
1044
27.2k
      cr_cell.detach();
1045
27.2k
      cr_create_at_level(cell_index, cr_max_level);
1046
27.2k
    }
1047
1048
18.5k
  return cr_max_level;
1049
18.5k
}
1050
1051
1052
unsigned int
1053
Partition::cr_get_backtrack_point()
1054
4.63M
{
1055
4.63M
  assert(cr_enabled);
1056
4.63M
  CR_BTInfo info;
1057
4.63M
  info.created_trail_index = cr_created_trail.size();
1058
4.63M
  info.splitted_level_trail_index = cr_splitted_level_trail.size();
1059
4.63M
  cr_bt_info.push_back(info);
1060
4.63M
  return cr_bt_info.size()-1;
1061
4.63M
}
1062
1063
1064
void
1065
Partition::cr_goto_backtrack_point(const unsigned int btpoint)
1066
2.44M
{
1067
2.44M
  assert(cr_enabled);
1068
2.44M
  assert(btpoint < cr_bt_info.size());
1069
4.93M
  while(cr_created_trail.size() > cr_bt_info[btpoint].created_trail_index)
1070
2.49M
    {
1071
2.49M
      const unsigned int cell_index = cr_created_trail.back();
1072
2.49M
      cr_created_trail.pop_back();
1073
2.49M
      CRCell& cr_cell = cr_cells[cell_index];
1074
2.49M
      assert(cr_cell.level != UINT_MAX);
1075
2.49M
      assert(cr_cell.prev_next_ptr);
1076
2.49M
      cr_cell.detach();
1077
2.49M
    }
1078
1079
2.46M
  while(cr_splitted_level_trail.size() >
1080
2.46M
        cr_bt_info[btpoint].splitted_level_trail_index)
1081
15.6k
    {
1082
15.6k
      const unsigned int dest_level = cr_splitted_level_trail.back();
1083
15.6k
      cr_splitted_level_trail.pop_back();
1084
15.6k
      assert(cr_max_level > 0);
1085
15.6k
      assert(dest_level < cr_max_level);
1086
38.3k
      while(cr_levels[cr_max_level]) {
1087
22.7k
        CRCell *cr_cell = cr_levels[cr_max_level];
1088
22.7k
        cr_cell->detach();
1089
22.7k
        cr_create_at_level(cr_cell - cr_cells, dest_level);
1090
22.7k
      }
1091
15.6k
      cr_max_level--;
1092
15.6k
    }
1093
2.44M
  cr_bt_info.resize(btpoint);
1094
2.44M
}
1095
1096
1097
void
1098
Partition::cr_create_at_level(const unsigned int cell_index,
1099
                              const unsigned int level)
1100
2.57M
{
1101
2.57M
  assert(cr_enabled);
1102
2.57M
  assert(cell_index < N);
1103
2.57M
  assert(level < N);
1104
2.57M
  CRCell& cr_cell = cr_cells[cell_index];
1105
2.57M
  assert(cr_cell.level == UINT_MAX);
1106
2.57M
  assert(cr_cell.next == 0);
1107
2.57M
  assert(cr_cell.prev_next_ptr == 0);
1108
2.57M
  if(cr_levels[level])
1109
2.53M
    cr_levels[level]->prev_next_ptr = &(cr_cell.next);
1110
2.57M
  cr_cell.next = cr_levels[level];
1111
2.57M
  cr_levels[level] = &cr_cell;
1112
2.57M
  cr_cell.prev_next_ptr = &cr_levels[level];
1113
2.57M
  cr_cell.level = level;
1114
2.57M
}
1115
1116
1117
void
1118
Partition::cr_create_at_level_trailed(const unsigned int cell_index,
1119
                                      const unsigned int level)
1120
2.52M
{
1121
2.52M
  assert(cr_enabled);
1122
2.52M
  cr_create_at_level(cell_index, level);
1123
2.52M
  cr_created_trail.push_back(cell_index);
1124
2.52M
}
1125
1126
1127
} // namespace bliss