Thrill  0.1
prefix_quadrupling.cpp
Go to the documentation of this file.
1 /*******************************************************************************
2  * examples/suffix_sorting/prefix_quadrupling.cpp
3  *
4  * Part of Project Thrill - http://project-thrill.org
5  *
6  * Copyright (C) 2016 Florian Kurpicz <[email protected]>
7  *
8  * All rights reserved. Published under the BSD-2 license in the LICENSE file.
9  ******************************************************************************/
10 
11 #include <examples/suffix_sorting/check_sa.hpp>
12 #include <examples/suffix_sorting/suffix_sorting.hpp>
13 
14 #include <thrill/api/cache.hpp>
15 #include <thrill/api/collapse.hpp>
16 #include <thrill/api/dia.hpp>
17 #include <thrill/api/generate.hpp>
18 #include <thrill/api/max.hpp>
19 #include <thrill/api/merge.hpp>
20 #include <thrill/api/prefix_sum.hpp>
21 #include <thrill/api/print.hpp>
22 #include <thrill/api/size.hpp>
23 #include <thrill/api/sort.hpp>
24 #include <thrill/api/sum.hpp>
25 #include <thrill/api/union.hpp>
26 #include <thrill/api/window.hpp>
27 #include <thrill/api/zip.hpp>
28 #include <thrill/common/logger.hpp>
29 #include <thrill/common/uint_types.hpp>
30 
31 #include <algorithm>
32 #include <limits>
33 #include <random>
34 #include <stdexcept>
35 #include <string>
36 #include <tuple>
37 #include <utility>
38 #include <vector>
39 
40 namespace examples {
41 namespace suffix_sorting {
42 
43 using namespace thrill; // NOLINT
44 using thrill::common::RingBuffer;
45 
46 //! A pair (index, rank)
47 template <typename Index>
48 struct IndexRank {
49  Index index;
50  Index rank;
51 
52  friend std::ostream& operator << (std::ostream& os, const IndexRank& ri) {
53  return os << '(' << ri.index << '|' << ri.rank << ')';
54  }
55 } TLX_ATTRIBUTE_PACKED;
56 
57 //! A pair (index, rank)
58 template <typename Index>
59 struct IndexQuadRank {
60  Index index;
61  Index rank[4];
62 
63  bool operator == (const IndexQuadRank& b) const {
64  return std::tie(rank[0], rank[1], rank[2], rank[3])
65  == std::tie(b.rank[0], b.rank[1], b.rank[2], b.rank[3]);
66  }
67 
68  bool operator < (const IndexQuadRank& b) const {
69  return rank[0] < b.rank[0] ||
70  (rank[0] == b.rank[0] && rank[1] < b.rank[1]) ||
71  (rank[0] == b.rank[0] && rank[1] == b.rank[1]
72  && rank[2] < b.rank[2]) ||
73  (rank[0] == b.rank[0] && rank[1] == b.rank[1]
74  && rank[2] == b.rank[2] && rank[3] < b.rank[3]) ||
75  (rank[0] == b.rank[0] && rank[1] == b.rank[1]
76  && rank[2] == b.rank[2] && rank[3] == b.rank[3]
77  && index > b.index);
78  }
79 
80  friend std::ostream& operator << (std::ostream& os, const IndexQuadRank& ri) {
81  return os << '(' << ri.index << '|' << ri.rank[0] << '|' << ri.rank[1]
82  << '|' << ri.rank[2] << '|' << ri.rank[3] << ')';
83  }
84 } TLX_ATTRIBUTE_PACKED;
85 
86 template <typename AlphabetType, typename Index>
87 struct QuadCharIndex {
88  AlphabetType ch[4];
89  Index index;
90 
91  bool operator == (const QuadCharIndex& b) const {
92  return std::tie(ch[0], ch[1], ch[2], ch[3])
93  == std::tie(b.ch[0], b.ch[1], b.ch[2], b.ch[3]);
94  }
95 
96  bool operator < (const QuadCharIndex& b) const {
97  return std::tie(ch[0], ch[1], ch[2], ch[3])
98  < std::tie(b.ch[0], b.ch[1], b.ch[2], b.ch[3]);
99  }
100 
101  friend std::ostream& operator << (std::ostream& os, const QuadCharIndex& chars) {
102  return os << '[' << chars.index << ": " << chars.ch[0] << ',' << chars.ch[1]
103  << ',' << chars.ch[2] << ',' << chars.ch[3] << ']';
104  }
105 } TLX_ATTRIBUTE_PACKED;
106 
107 enum class Status : uint8_t {
108  UNDECIDED = 0,
109  UNIQUE = 1,
110  FULLY_DISCARDED = 2
111 };
112 
113 //! A triple (index, rank, status)
114 template <typename Index>
115 struct IndexRankStatus {
116  Index index;
117  Index rank;
118  Status status;
119 
120  //! Two IndexRandStatuses are equal iff their ranks are equal.
121  bool operator == (const IndexRankStatus& b) const {
122  return rank == b.rank;
123  }
124 
125  //! A IndexRankStatus is smaller than another iff either its fist rank is
126  //! smaller or if both ranks are equal and the first index is _larger_ than
127  //! the second. The _larger_ is due to suffixes with larger index being
128  //! smaller.
129  bool operator < (const IndexRankStatus& b) const {
130  return rank < b.rank || (rank == b.rank && index > b.index);
131  }
132 
133  friend std::ostream& operator << (std::ostream& os, const IndexRankStatus& irs) {
134  return os << "(i: " << irs.index << "| r: " << irs.rank << "| s: "
135  << static_cast<uint8_t>(irs.status) << ")";
136  }
137 } TLX_ATTRIBUTE_PACKED;
138 
139 //! A triple (index, rank [4], status)
140 template <typename Index>
141 struct IndexQuadRankStatus {
142  Index index;
143  Index rank[4];
144  Status status;
145 
146  friend std::ostream& operator << (std::ostream& os, const IndexQuadRankStatus& iqrs) {
147  return os << "(i: " << iqrs.index << "| r1: " << iqrs.rank[0] << "| r2: "
148  << iqrs.rank[1] << "| r3: " << iqrs.rank[2] << "| r4: "
149  << iqrs.rank[3] << "| s: " << static_cast<uint8_t>(iqrs.status) << ")";
150  }
151 } TLX_ATTRIBUTE_PACKED;
152 
153 //! A triple with index (index, rank_1, rank_2, rank_3)
154 template <typename Index>
155 struct Index3Rank {
156  Index index;
157  Index rank1;
158  Index rank2;
159  Index rank3;
160 
161  friend std::ostream& operator << (std::ostream& os, const Index3Rank& irrr) {
162  return os << "( i: " << irrr.index << "| r1: " << irrr.rank1
163  << "| r2: " << irrr.rank2 << "| r3: " << irrr.rank3 << ")";
164  }
165 } TLX_ATTRIBUTE_PACKED;
166 
167 template <typename Index, typename InputDIA>
168 DIA<Index> PrefixQuadruplingDiscarding(const InputDIA& input_dia, size_t input_size, bool packed) {
169  LOG1 << "Running PrefixQuadruplingDiscarding";
170 
171  using Char = typename InputDIA::ValueType;
172  using IndexRank = suffix_sorting::IndexRank<Index>;
173  using IndexQuadRank = suffix_sorting::IndexQuadRank<Index>;
174  using QuadCharIndex = suffix_sorting::QuadCharIndex<Char, Index>;
175  using IndexQuadRankStatus = suffix_sorting::IndexQuadRankStatus<Index>;
176  using IndexRankStatus = suffix_sorting::IndexRankStatus<Index>;
177  using Index3Rank = suffix_sorting::Index3Rank<Index>;
178 
179  size_t input_bit_size = sizeof(Char) << 3;
180  size_t k_fitting = sizeof(Index) / sizeof(Char);
181 
182  if (packed && k_fitting <= 4) {
183  LOG1 << "Ignoring packing, as it would have no effect."
184  << " More than 4 characters must be packed to make a difference."
185  << " Currently " << k_fitting << " characters are packed.";
186  packed = false;
187  }
188  if (packed && k_fitting < 16) {
189  LOG1 << k_fitting << " bytes can be packed into the space used by one suffix array index.";
190  LOG1 << "This is not enough to skip quadrupling steps. We can decide more suffixes initially, though.";
191  }
192 
193  size_t iteration = 1;
194 
195  if (packed && k_fitting >= 16) {
196  iteration = 0;
197  size_t tmp = k_fitting;
198  while (tmp >>= 1) ++iteration;
199  iteration >>= 1; // Compute log_4(k_fitting)
200  }
201 
202  DIA<IndexRank> names;
203 
204  if (packed) {
205  auto chars_sorted =
206  input_dia
207  .template FlatWindow<IndexRank>(
208  k_fitting,
209  [=](size_t index, const RingBuffer<Char>& rb, auto emit) {
210  size_t result = rb[0];
211  for (size_t i = 1; i < k_fitting; ++i)
212  result = (result << input_bit_size) | rb[i];
213  emit(IndexRank { Index(index), Index(result) });
214  if (index + k_fitting == input_size) {
215  for (size_t i = 1; i < k_fitting; ++i) {
216  result = rb[i];
217  for (size_t j = i + 1; j < k_fitting; ++j)
218  result = (result << input_bit_size) | rb[j];
219  result <<= i * input_bit_size;
220  emit(IndexRank { Index(index + i), Index(result) });
221  }
222  }
223  })
224  .Sort([](const IndexRank& a, const IndexRank& b) {
225  return a.rank < b.rank;
226  });
227 
228  if (debug_print)
229  chars_sorted.Keep().Print("chars_sorted packed");
230 
231  names =
232  chars_sorted
233  .template FlatWindow<IndexRank>(
234  2,
235  [](size_t index, const RingBuffer<IndexRank>& rb, auto emit) {
236  if (index == 0)
237  emit(IndexRank { rb[0].index, Index(1) });
238  emit(IndexRank {
239  rb[1].index, Index(rb[0].rank == rb[1].rank ? 0 : index + 2)
240  });
241  })
242  .PrefixSum([](const IndexRank& a, const IndexRank& b) {
243  return IndexRank {
244  b.index,
245  (a.rank > b.rank ? a.rank : b.rank)
246  };
247  });
248  }
249  else {
250  auto chars_sorted =
251  input_dia
252  .template FlatWindow<QuadCharIndex>(
253  4,
254  [=](size_t index, const RingBuffer<Char>& rb, auto emit) {
255  emit(QuadCharIndex {
256  { rb[0], rb[1], rb[2], rb[3] }, Index(index)
257  });
258  },
259  [=](size_t index, const RingBuffer<Char>& rb, auto emit) {
260  if (index == input_size - 3) {
261  emit(QuadCharIndex {
262  { rb[0], rb[1], rb[2],
263  std::numeric_limits<Char>::lowest() },
264  Index(index)
265  });
266  emit(QuadCharIndex {
267  { rb[1], rb[2],
268  std::numeric_limits<Char>::lowest(),
269  std::numeric_limits<Char>::lowest() },
270  Index(index + 1)
271  });
272  emit(QuadCharIndex {
273  { rb[2],
274  std::numeric_limits<Char>::lowest(),
275  std::numeric_limits<Char>::lowest(),
276  std::numeric_limits<Char>::lowest() },
277  Index(index + 2)
278  });
279  }
280  })
281  .Sort();
282 
283  if (debug_print)
284  chars_sorted.Keep().Print("chars_sorted");
285 
286  names =
287  chars_sorted
288  .template FlatWindow<IndexRank>(
289  2,
290  [](size_t index, const RingBuffer<QuadCharIndex>& rb, auto emit) {
291  if (index == 0)
292  emit(IndexRank { rb[0].index, Index(1) });
293  if (rb[0] == rb[1])
294  emit(IndexRank { rb[1].index, Index(0) });
295  else
296  emit(IndexRank { rb[1].index, Index(index + 2) });
297  })
298  .PrefixSum([](const IndexRank& a, const IndexRank& b) {
299  return IndexRank {
300  b.index,
301  (a.rank > b.rank ? a.rank : b.rank)
302  };
303  });
304  }
305 
306  if (debug_print)
307  names.Keep().Print("names");
308 
309  auto names_unique =
310  names
311  .template FlatWindow<IndexRankStatus>(
312  3,
313  [=](size_t index, const RingBuffer<IndexRank>& rb, auto emit) {
314  if (index == 0) {
315  Status status = rb[0].rank != rb[1].rank ? Status::UNIQUE : Status::UNDECIDED;
316  emit(IndexRankStatus { rb[0].index, rb[0].rank, status });
317  }
318  if (rb[0].rank != rb[1].rank && rb[1].rank != rb[2].rank)
319  emit(IndexRankStatus { rb[1].index, rb[1].rank, Status::UNIQUE });
320  else
321  emit(IndexRankStatus { rb[1].index, rb[1].rank, Status::UNDECIDED });
322  if (index == input_size - 3) {
323  Status status = rb[1].rank != rb[2].rank ? Status::UNIQUE : Status::UNDECIDED;
324  emit(IndexRankStatus { rb[2].index, rb[2].rank, status });
325  }
326  });
327 
328  auto names_unique_sorted =
329  names_unique.Keep()
330  .Sort([iteration](const IndexRankStatus& a, const IndexRankStatus& b) {
331  Index mod_mask = (Index(1) << (iteration << 1)) - 1;
332  Index div_mask = ~mod_mask;
333 
334  if ((a.index & mod_mask) == (b.index & mod_mask))
335  return (a.index & div_mask) < (b.index & div_mask);
336  else
337  return (a.index & mod_mask) < (b.index & mod_mask);
338  });
339 
340  if (debug_print)
341  names_unique_sorted.Keep().Print("Names unique sorted");
342 
343  std::vector<DIA<IndexQuadRank> > fully_discarded;
344 
345  while (true) {
346  size_t next_index = size_t(1) << (iteration << 1);
347  size_t names_size = names_unique_sorted.Keep().Size();
348 
349  ++iteration;
350 
351  auto discarded_names =
352  names_unique_sorted.Keep()
353  .template FlatWindow<IndexQuadRankStatus>(
354  5,
355  [=](size_t index, const RingBuffer<IndexRankStatus>& rb, auto emit) {
356  // Discarded names (we need to change the status since we remove it one step later)
357  if (index == 0) {
358  if (rb[0].status == Status::UNIQUE)
359  emit(IndexQuadRankStatus {
360  rb[0].index,
361  { rb[0].rank, Index(0), Index(0), Index(0) },
362  Status::FULLY_DISCARDED
363  });
364  if (rb[1].status == Status::UNIQUE)
365  emit(IndexQuadRankStatus {
366  rb[1].index,
367  { rb[1].rank, Index(0), Index(0), Index(0) },
368  Status::FULLY_DISCARDED
369  });
370  if (rb[2].status == Status::UNIQUE)
371  emit(IndexQuadRankStatus {
372  rb[2].index,
373  { rb[2].rank, Index(0), Index(0), Index(0) },
374  Status::FULLY_DISCARDED
375  });
376  if (rb[3].status == Status::UNIQUE)
377  emit(IndexQuadRankStatus {
378  rb[3].index,
379  { rb[3].rank, Index(0), Index(0), Index(0) },
380  Status::FULLY_DISCARDED
381  });
382  }
383  if (rb[4].status == Status::UNIQUE) {
384  if (rb[0].status == Status::UNIQUE || rb[1].status == Status::UNIQUE ||
385  rb[2].status == Status::UNIQUE || rb[3].status == Status::UNIQUE)
386  emit(IndexQuadRankStatus {
387  rb[4].index,
388  { rb[4].rank, Index(0), Index(0), Index(0) },
389  Status::FULLY_DISCARDED
390  });
391  else
392  emit(IndexQuadRankStatus {
393  rb[4].index,
394  { rb[4].rank, Index(0), Index(0), Index(0) },
395  Status::UNIQUE
396  });
397  }
398  if (rb[0].status == Status::UNDECIDED) {
399  Index rank1 = (rb[0].index + Index(next_index) == rb[1].index) ? rb[1].rank : Index(0);
400  Index rank2 = (rb[0].index + Index(2 * next_index) == rb[2].index) ? rb[2].rank : Index(0);
401  Index rank3 = (rb[0].index + Index(3 * next_index) == rb[3].index) ? rb[3].rank : Index(0);
402  emit(IndexQuadRankStatus { rb[0].index, { rb[0].rank, rank1, rank2, rank3 }, Status::UNDECIDED });
403  }
404  },
405  [=](size_t index, const RingBuffer<IndexRankStatus>& rb, auto emit) {
406  if (index == 0) {
407  if (rb[0].status == Status::UNIQUE)
408  emit(IndexQuadRankStatus {
409  rb[0].index,
410  { rb[0].rank, Index(0), Index(0), Index(0) },
411  Status::FULLY_DISCARDED
412  });
413  else // (rb[0].status == Status::UNDECIDED)
414  emit(IndexQuadRankStatus {
415  rb[0].index,
416  { rb[0].rank, Index(0), Index(0), Index(0) },
417  Status::UNDECIDED
418  });
419 
420  if (rb.size() > 1 && rb[1].status == Status::UNIQUE)
421  emit(IndexQuadRankStatus {
422  rb[1].index,
423  { rb[1].rank, Index(0), Index(0), Index(0) },
424  Status::FULLY_DISCARDED
425  });
426  else if (rb.size() > 1) // && rb[1].status == Status::UNDECIDED)
427  emit(IndexQuadRankStatus {
428  rb[1].index,
429  { rb[1].rank, Index(0), Index(0), Index(0) },
430  Status::UNDECIDED
431  });
432 
433  if (rb.size() > 2 && rb[2].status == Status::UNIQUE)
434  emit(IndexQuadRankStatus {
435  rb[2].index,
436  { rb[2].rank, Index(0), Index(0), Index(0) },
437  Status::FULLY_DISCARDED
438  });
439  else if (rb.size() > 2) // && rb[2].status == Status::UNDECIDED)
440  emit(IndexQuadRankStatus {
441  rb[2].index,
442  { rb[2].rank, Index(0), Index(0), Index(0) },
443  Status::UNDECIDED
444  });
445 
446  if (rb.size() > 3 && rb[3].status == Status::UNIQUE)
447  emit(IndexQuadRankStatus {
448  rb[3].index,
449  { rb[3].rank, Index(0), Index(0), Index(0) },
450  Status::FULLY_DISCARDED
451  });
452  else if (rb.size() > 3) // && rb[3].status == Status::UNDECIDED)
453  emit(IndexQuadRankStatus {
454  rb[3].index,
455  { rb[3].rank, Index(0), Index(0), Index(0) },
456  Status::UNDECIDED
457  });
458  }
459  if (index == names_size - 4) {
460  Index rank1;
461  Index rank2;
462  Index rank3;
463 
464  if (rb[0].status == Status::UNDECIDED) {
465  rank1 = (rb[0].index + Index(next_index) == rb[1].index) ? rb[1].rank : Index(0);
466  rank2 = (rb[0].index + Index(2 * next_index) == rb[2].index) ? rb[2].rank : Index(0);
467  rank3 = (rb[0].index + Index(3 * next_index) == rb[3].index) ? rb[3].rank : Index(0);
468  emit(IndexQuadRankStatus { rb[0].index, { rb[0].rank, rank1, rank2, rank3 }, Status::UNDECIDED });
469  }
470  if (rb[1].status == Status::UNDECIDED) {
471  rank1 = (rb[1].index + Index(next_index) == rb[2].index) ? rb[2].rank : Index(0);
472  rank2 = (rb[1].index + Index(2 * next_index) == rb[3].index) ? rb[3].rank : Index(0);
473  emit(IndexQuadRankStatus { rb[1].index, { rb[1].rank, rank1, rank2, Index(0) }, Status::UNDECIDED });
474  }
475  if (rb[2].status == Status::UNDECIDED) {
476  rank1 = (rb[2].index + Index(next_index) == rb[3].index) ? rb[3].rank : Index(0);
477  emit(IndexQuadRankStatus { rb[2].index, { rb[2].rank, rank1, Index(0), Index(0) }, Status::UNDECIDED });
478  }
479  if (rb[3].status == Status::UNDECIDED)
480  emit(IndexQuadRankStatus { rb[3].index, { rb[3].rank, Index(0), Index(0), Index(0) }, Status::UNDECIDED });
481  }
482  });
483 
484  auto new_decided =
485  discarded_names.Keep()
486  .Filter([](const IndexQuadRankStatus& iqrs) {
487  return iqrs.status == Status::FULLY_DISCARDED;
488  })
489  .Map([](const IndexQuadRankStatus& iqrs) {
490  return IndexQuadRank {
491  iqrs.index,
492  { iqrs.rank[0], iqrs.rank[1], iqrs.rank[2], iqrs.rank[3] }
493  };
494  });
495 
496  if (debug_print)
497  new_decided.Keep().Print("new_decided");
498 
499  auto partial_discarded =
500  discarded_names.Keep()
501  .Filter([](const IndexQuadRankStatus& iqrs) {
502  return iqrs.status == Status::UNIQUE;
503  })
504  .Map([](const IndexQuadRankStatus& iqrs) {
505  return IndexRankStatus {
506  iqrs.index,
507  iqrs.rank[0],
508  Status::UNIQUE
509  };
510  });
511 
512  if (debug_print)
513  partial_discarded.Keep().Print("partial_discarded");
514 
515  auto undiscarded =
516  discarded_names
517  .Filter([](const IndexQuadRankStatus& iqrs) {
518  return iqrs.status == Status::UNDECIDED;
519  })
520  .Map([](const IndexQuadRankStatus& iqrs) {
521  return IndexQuadRank {
522  iqrs.index,
523  { iqrs.rank[0], iqrs.rank[1], iqrs.rank[2], iqrs.rank[3] }
524  };
525  })
526  .Sort([](const IndexQuadRank& a, const IndexQuadRank& b) {
527  return a < b;
528  });
529 
530  if (debug_print)
531  undiscarded.Keep().Print("undiscarded");
532 
533  fully_discarded.emplace_back(new_decided.Cache());
534 
535  size_t duplicates = undiscarded.Keep().Size();
536 
537  if (input_dia.context().my_rank() == 0) {
538  sLOG1 << "iteration" << iteration - 1
539  << "duplicates" << duplicates;
540  }
541 
542  if (duplicates == 0) {
543  // auto res = Union(fully_discarded);
544  auto sa =
545  Union(fully_discarded)
546  .Sort()
547  .Map([](const IndexQuadRank& iqr) {
548  return iqr.index;
549  });
550  return sa.Collapse();
551  }
552 
553  auto new_ranks =
554  undiscarded
555  .template FlatWindow<Index3Rank>(
556  2,
557  [=](size_t index, const RingBuffer<IndexQuadRank>& rb, auto emit) {
558  if (index == 0) {
559  emit(Index3Rank { rb[0].index, Index(0), Index(0), rb[0].rank[0] });
560  }
561  Index i = rb[0].rank[0] == rb[1].rank[0] ? Index(0) : Index(index + 1);
562  Index r;
563  if (rb[0].rank[0] != rb[1].rank[0])
564  r = Index(index + 1);
565  else
566  r = (rb[0].rank[1] == rb[1].rank[1] &&
567  rb[0].rank[2] == rb[1].rank[2] &&
568  rb[0].rank[3] == rb[1].rank[3]) ? Index(0) : Index(index + 1);
569  emit(Index3Rank { rb[1].index, i, r, rb[1].rank[0] });
570  },
571  [](size_t index, const RingBuffer<IndexQuadRank>& rb, auto emit) {
572  if (index == 0)
573  emit(Index3Rank { rb[0].index, Index(0), Index(0), rb[0].rank[0] });
574  })
575  .PrefixSum([](const Index3Rank& a, const Index3Rank& b) {
576  return Index3Rank {
577  b.index,
578  std::max<Index>(a.rank1, b.rank1),
579  std::max<Index>(a.rank2, b.rank2),
580  b.rank3
581  };
582  })
583  .Map([](const Index3Rank& ir) {
584  return IndexRank { ir.index, ir.rank3 + (ir.rank2 - ir.rank1) };
585  });
586 
587  // new_ranks.Keep().Print("new_ranks");
588 
589  names_unique =
590  new_ranks
591  .template FlatWindow<IndexRankStatus>(
592  3,
593  [=](size_t index, const RingBuffer<IndexRank>& rb, auto emit) {
594  if (index == 0) {
595  Status status = rb[0].rank != rb[1].rank ? Status::UNIQUE : Status::UNDECIDED;
596  emit(IndexRankStatus { rb[0].index, rb[0].rank, status });
597  }
598  if (rb[0].rank != rb[1].rank && rb[1].rank != rb[2].rank)
599  emit(IndexRankStatus { rb[1].index, rb[1].rank, Status::UNIQUE });
600  else
601  emit(IndexRankStatus { rb[1].index, rb[1].rank, Status::UNDECIDED });
602  if (index == duplicates - 3) {
603  Status status = rb[1].rank != rb[2].rank ? Status::UNIQUE : Status::UNDECIDED;
604  emit(IndexRankStatus { rb[2].index, rb[2].rank, status });
605  }
606  },
607  [](size_t index, const RingBuffer<IndexRank>& rb, auto emit) {
608  if (index == 0) { // We know that there are exactly 2 names
609  emit(IndexRankStatus { rb[0].index, rb[0].rank, Status::UNIQUE });
610  emit(IndexRankStatus { rb[1].index, rb[1].rank, Status::UNIQUE });
611  }
612  });
613 
614  if (debug_print)
615  names_unique.Keep().Print("names_unique");
616 
617  names_unique_sorted =
618  names_unique
619  .Union(partial_discarded)
620  .Sort([iteration](const IndexRankStatus& a, const IndexRankStatus& b) {
621  Index mod_mask = (Index(1) << (iteration << 1)) - 1;
622  Index div_mask = ~mod_mask;
623 
624  if ((a.index & mod_mask) == (b.index & mod_mask))
625  return (a.index & div_mask) < (b.index & div_mask);
626  else
627  return (a.index & mod_mask) < (b.index & mod_mask);
628  });
629 
630  if (debug_print)
631  names_unique_sorted.Keep().Print("names_unique_sorted");
632  }
633 }
634 
635 template <typename Index, typename InputDIA>
636 DIA<Index> PrefixQuadrupling(const InputDIA& input_dia, size_t input_size, bool packed) {
637  LOG1 << "Running PrefixQuadrupling";
638 
639  using Char = typename InputDIA::ValueType;
640  using IndexRank = suffix_sorting::IndexRank<Index>;
641  using IndexQuadRank = suffix_sorting::IndexQuadRank<Index>;
642  using QuadCharIndex = suffix_sorting::QuadCharIndex<Char, Index>;
643 
644  size_t input_bit_size = sizeof(Char) << 3;
645  size_t k_fitting = sizeof(Index) / sizeof(Char);
646 
647  if (packed && k_fitting <= 4) {
648  LOG1 << "Ignoring packing, as it would have no effect."
649  << " More than 4 characters must be packed to make a difference."
650  << " Currently " << k_fitting << " characters are packed.";
651  packed = false;
652  }
653  if (packed && k_fitting < 16) {
654  LOG1 << k_fitting << " bytes can be packed into the space used by one suffix array index.";
655  LOG1 << "This is not enough to skip quadrupling steps. We can decide more suffixes initially, though.";
656  }
657 
658  size_t iteration = 1;
659 
660  if (packed && k_fitting >= 16) {
661  iteration = 0;
662  size_t tmp = k_fitting;
663  while (tmp >>= 1) ++iteration;
664  iteration >>= 1; // Compute log_4(k_fitting)
665  }
666 
667  DIA<IndexRank> names;
668 
669  if (packed) {
670  auto chars_sorted =
671  input_dia
672  .template FlatWindow<IndexRank>(
673  k_fitting,
674  [=](size_t index, const RingBuffer<Char>& rb, auto emit) {
675  size_t result = rb[0];
676  for (size_t i = 1; i < k_fitting; ++i)
677  result = (result << input_bit_size) | rb[i];
678  emit(IndexRank { Index(index), Index(result) });
679  if (index + k_fitting == input_size) {
680  for (size_t i = 1; i < k_fitting; ++i) {
681  result = rb[i];
682  for (size_t j = i + 1; j < k_fitting; ++j)
683  result = (result << input_bit_size) | rb[j];
684  result <<= i * input_bit_size;
685  emit(IndexRank { Index(index + i), Index(result) });
686  }
687  }
688  })
689  .Sort([](const IndexRank& a, const IndexRank& b) {
690  return a.rank < b.rank;
691  });
692 
693  if (debug_print)
694  chars_sorted.Keep().Print("chars_sorted packed");
695 
696  names =
697  chars_sorted
698  .template FlatWindow<IndexRank>(
699  2,
700  [](size_t index, const RingBuffer<IndexRank>& rb, auto emit) {
701  if (index == 0)
702  emit(IndexRank { rb[0].index, Index(1) });
703  emit(IndexRank {
704  rb[1].index, Index(rb[0].rank == rb[1].rank ? 0 : index + 2)
705  });
706  });
707  }
708  else {
709  auto chars_sorted =
710  input_dia
711  .template FlatWindow<QuadCharIndex>(
712  4,
713  [=](size_t index, const RingBuffer<Char>& rb, auto emit) {
714  emit(QuadCharIndex {
715  { rb[0], rb[1], rb[2], rb[3] }, Index(index)
716  });
717  },
718  [=](size_t index, const RingBuffer<Char>& rb, auto emit) {
719  if (index == input_size - 3) {
720  emit(QuadCharIndex {
721  { rb[0], rb[1], rb[2],
722  std::numeric_limits<Char>::lowest() },
723  Index(index)
724  });
725  emit(QuadCharIndex {
726  { rb[1], rb[2],
727  std::numeric_limits<Char>::lowest(),
728  std::numeric_limits<Char>::lowest() },
729  Index(index + 1)
730  });
731  emit(QuadCharIndex {
732  { rb[2],
733  std::numeric_limits<Char>::lowest(),
734  std::numeric_limits<Char>::lowest(),
735  std::numeric_limits<Char>::lowest() },
736  Index(index + 2)
737  });
738  }
739  })
740  .Sort();
741 
742  names =
743  chars_sorted
744  .template FlatWindow<IndexRank>(
745  2,
746  [](size_t index, const RingBuffer<QuadCharIndex>& rb, auto emit) {
747  if (index == 0)
748  emit(IndexRank { rb[0].index, Index(1) });
749  if (rb[0] == rb[1])
750  emit(IndexRank { rb[1].index, Index(0) });
751  else
752  emit(IndexRank { rb[1].index, Index(index + 2) });
753  });
754  }
755 
756  auto number_duplicates =
757  names.Keep()
758  .Filter([](const IndexRank& ir) {
759  return ir.rank == Index(0);
760  })
761  .Size();
762 
763  if (number_duplicates == 0) {
764  if (input_dia.context().my_rank() == 0)
765  sLOG1 << "Finished before doubling in loop";
766 
767  auto sa =
768  names
769  .Map([](const IndexRank& ir) {
770  return ir.index;
771  });
772 
773  return sa.Collapse();
774  }
775 
776  names =
777  names
778  .PrefixSum(
779  [](const IndexRank& a, const IndexRank& b) {
780  return IndexRank { b.index, std::max(a.rank, b.rank) };
781  });
782 
783  if (debug_print)
784  names.Keep().Print("names before loop");
785 
786  while (true) {
787  auto names_sorted =
788  names
789  .Sort([iteration](const IndexRank& a, const IndexRank& b) {
790  Index mod_mask = (Index(1) << (iteration << 1)) - 1;
791  Index div_mask = ~mod_mask;
792 
793  if ((a.index & mod_mask) == (b.index & mod_mask))
794  return (a.index & div_mask) < (b.index & div_mask);
795  else
796  return (a.index & mod_mask) < (b.index & mod_mask);
797  });
798 
799  size_t next_index = size_t(1) << (iteration << 1);
800  ++iteration;
801 
802  auto triple_sorted =
803  names_sorted
804  .template FlatWindow<IndexQuadRank>(
805  4,
806  [=](size_t /*index*/, const RingBuffer<IndexRank>& rb, auto emit) {
807  Index rank1 = (rb[0].index + Index(next_index) == rb[1].index) ? rb[1].rank : Index(0);
808  Index rank2 = (rb[0].index + Index(2 * next_index) == rb[2].index) ? rb[2].rank : Index(0);
809  Index rank3 = (rb[0].index + Index(3 * next_index) == rb[3].index) ? rb[3].rank : Index(0);
810  emit(IndexQuadRank { rb[0].index, { rb[0].rank, rank1, rank2, rank3 }
811  });
812  },
813  [=](size_t index, const RingBuffer<IndexRank>& rb, auto emit) {
814  if (index == input_size - 3) {
815  Index rank1 = (rb[0].index + Index(next_index) == rb[1].index) ? rb[1].rank : Index(0);
816  Index rank2 = (rb[0].index + Index(2 * next_index) == rb[2].index) ? rb[2].rank : Index(0);
817  emit(IndexQuadRank { rb[0].index, { rb[0].rank, rank1, rank2, Index(0) }
818  });
819 
820  rank1 = (rb[1].index + Index(next_index) == rb[2].index) ? rb[2].rank : Index(0);
821  emit(IndexQuadRank { rb[1].index, { rb[1].rank, rank1, Index(0), Index(0) }
822  });
823  emit(IndexQuadRank { rb[2].index, { rb[2].rank, Index(0), Index(0), Index(0) }
824  });
825  }
826  })
827  .Sort();
828 
829  names =
830  triple_sorted
831  .template FlatWindow<IndexRank>(
832  2,
833  [](size_t index, const RingBuffer<IndexQuadRank>& rb, auto emit) {
834  if (index == 0) emit(IndexRank { rb[0].index, Index(1) });
835 
836  if (rb[0] == rb[1] && (rb[0].rank[1] != Index(0) || rb[0].rank[2] != Index(0) || rb[0].rank[3] != Index(0)))
837  emit(IndexRank { rb[1].index, Index(0) });
838  else
839  emit(IndexRank { rb[1].index, Index(1) });
840  });
841 
842  number_duplicates =
843  names.Keep()
844  .Filter([](const IndexRank& ir) {
845  return ir.rank == Index(0);
846  })
847  .Size();
848 
849  if (input_dia.context().my_rank() == 0) {
850  sLOG1 << "iteration" << iteration - 1
851  << "duplicates" << number_duplicates;
852  }
853 
854  if (number_duplicates == 0) {
855  auto sa =
856  names
857  .Map([](const IndexRank& ir) {
858  return ir.index;
859  });
860 
861  return sa.Collapse();
862  }
863 
864  names =
865  names
866  .PrefixSum([](const IndexRank& a, const IndexRank& b) {
867  return IndexRank { b.index, std::max(a.rank, b.rank) };
868  });
869 
870  if (debug_print)
871  names.Keep().Print("names");
872  }
873 }
874 
875 template DIA<uint32_t> PrefixQuadrupling<uint32_t>(
876  const DIA<uint8_t>& input_dia, size_t input_size, bool packed);
877 
878 template DIA<uint32_t> PrefixQuadruplingDiscarding<uint32_t>(
879  const DIA<uint8_t>& input_dia, size_t input_size, bool packed);
880 
881 #if !THRILL_ON_TRAVIS
882 
883 template DIA<common::uint40> PrefixQuadrupling<common::uint40>(
884  const DIA<uint8_t>& input_dia, size_t input_size, bool packed);
885 
886 template DIA<common::uint40> PrefixQuadruplingDiscarding<common::uint40>(
887  const DIA<uint8_t>& input_dia, size_t input_size, bool packed);
888 
889 template DIA<uint64_t> PrefixQuadrupling<uint64_t>(
890  const DIA<uint8_t>& input_dia, size_t input_size, bool packed);
891 
892 template DIA<uint64_t> PrefixQuadruplingDiscarding<uint64_t>(
893  const DIA<uint8_t>& input_dia, size_t input_size, bool packed);
894 
895 #endif
896 
897 } // namespace suffix_sorting
898 } // namespace examples
899 
900 /******************************************************************************/
suffix_sorting.hpp
thrill::common::RingBuffer
tlx::RingBuffer< Type, Allocator > RingBuffer
Definition: ring_buffer.hpp:21
examples::suffix_sorting::PrefixQuadruplingDiscarding< uint64_t >
template DIA< uint64_t > PrefixQuadruplingDiscarding< uint64_t >(const DIA< uint8_t > &input_dia, size_t input_size, bool packed)
thrill::api::DIA
DIA is the interface between the user and the Thrill framework.
Definition: dia.hpp:141
max
static uint_pair max()
return an uint_pair instance containing the largest value possible
Definition: uint_types.hpp:226
examples::suffix_sorting::Status::UNDECIDED
thrill::api::Union
auto Union(const FirstDIA &first_dia, const DIAs &... dias)
Union is a LOp, which creates the union of all items from any number of DIAs as a single DIA...
Definition: union.hpp:319
tlx::RingBuffer
A ring (circular) buffer of static (non-growing) size.
Definition: ring_buffer.hpp:36
LOG1
#define LOG1
Definition: logger.hpp:28
examples::suffix_sorting::Status
Status
Definition: prefix_doubling.cpp:130
examples
Definition: bfs.hpp:21
examples::suffix_sorting::Status::UNIQUE
examples::suffix_sorting::IndexRank::index
Index index
Definition: construct_lcp.hpp:37
sLOG1
#define sLOG1
Definition: logger.hpp:38
thrill::api::DIA::PrefixSum
auto PrefixSum(const SumFunction &sum_function=SumFunction(), const ValueType &initial_element=ValueType()) const
PrefixSum is a DOp, which computes the (inclusive) prefix sum of all elements.
Definition: prefix_sum.hpp:134
examples::suffix_sorting::PrefixQuadrupling< uint32_t >
template DIA< uint32_t > PrefixQuadrupling< uint32_t >(const DIA< uint8_t > &input_dia, size_t input_size, bool packed)
prefix_sum.hpp
uint_types.hpp
examples::suffix_sorting::Status::FULLY_DISCARDED
thrill::api::DIA::Print
void Print(const std::string &name=std::string()) const
Print is an Action, which collects all data of the DIA at the worker 0 and prints using ostream seria...
Definition: print.hpp:50
thrill::api::DIA::Sort
auto Sort(const CompareFunction &compare_function=CompareFunction()) const
Sort is a DOp, which sorts a given DIA according to the given compare_function.
Definition: sort.hpp:800
thrill::api::DIA::Map
auto Map(const MapFunction &map_function) const
Map applies map_function : to each item of a DIA and delivers a new DIA contains the returned values...
Definition: dia.hpp:358
operator==
bool operator==(const uint_pair &b) const
equality checking operator
Definition: uint_types.hpp:175
examples::suffix_sorting::operator<<
std::ostream & operator<<(std::ostream &os, const Status &s)
Definition: prefix_doubling.cpp:136
examples::suffix_sorting::PrefixQuadrupling
DIA< Index > PrefixQuadrupling(const InputDIA &input_dia, size_t input_size, bool packed)
Definition: prefix_quadrupling.cpp:636
examples::suffix_sorting::TLX_ATTRIBUTE_PACKED
struct examples::suffix_sorting::IndexRank TLX_ATTRIBUTE_PACKED
examples::suffix_sorting::PrefixQuadruplingDiscarding
DIA< Index > PrefixQuadruplingDiscarding(const InputDIA &input_dia, size_t input_size, bool packed)
Definition: prefix_quadrupling.cpp:168
examples::suffix_sorting::debug_print
bool debug_print
Definition: suffix_sorting.cpp:47
check_sa.hpp
examples::suffix_sorting::PrefixQuadruplingDiscarding< uint32_t >
template DIA< uint32_t > PrefixQuadruplingDiscarding< uint32_t >(const DIA< uint8_t > &input_dia, size_t input_size, bool packed)
examples::suffix_sorting::IndexRank
A pair (index, rank)
Definition: construct_lcp.hpp:36
thrill::api::DIA::Keep
const DIA & Keep(size_t increase=1) const
Mark the referenced DIANode for keeping, which makes children not consume the data when executing...
Definition: dia.hpp:310
operator<
bool operator<(const uint_pair &b) const
less-than comparison operator
Definition: uint_types.hpp:187
examples::suffix_sorting::PrefixQuadrupling< uint64_t >
template DIA< uint64_t > PrefixQuadrupling< uint64_t >(const DIA< uint8_t > &input_dia, size_t input_size, bool packed)