15 回答

• 3

  你想看看Drs的Large-scale Genome Sequence Processing . 笠原和森下 .

  由四个核苷酸字母A，C，G和T组成的字符串可以特别编码为整数，以便更快地处理 . 基数排序是本书中讨论的众多算法之一;你应该能够根据这个问题调整已接受的答案，并看到一个很大的性能提升 .

  回复于 2024-04-30T12:25:00+08:00
• 1

  “Radix sorting with no extra space”是一篇解决您问题的文章 .

  回复于 2024-04-30T12:25:00+08:00
• 6

  Radix-Sort不是缓存意识，并不是大集合的最快排序算法 . 你可以看看：

  • ti7qsort . ti7qsort是整数的最快排序（可用于小型固定大小的字符串） .

  • Inline QSORT

  • String sorting

  在存储到排序数组之前，您还可以使用压缩并将DNA的每个字母编码为2位 .

  回复于 2024-04-30T12:25:00+08:00
• 55

  dsimcha的MSB基数排序看起来不错，但是Nils更接近问题的核心，观察到缓存局部性是在大问题规模上杀死你的东西 .

  我建议一个非常简单的方法：

  • 根据经验估算基数排序有效的最大尺寸 m .

  • 一次读取 m 元素的块，对它们进行基数排序，并将它们写出来（如果你有足够的内存，则写入内存缓冲区，否则要写入文件），直到你耗尽输入为止 .

  • Mergesort 生成的排序块 .

  Mergesort是我所知道的最适合缓存的排序算法："Read the next item from either array A or B, then write an item to the output buffer."它在磁带驱动器上高效运行 . 它确实需要 2n 空间来排序 n 项目，但我敢打赌，你将看到的大大改进的缓存局部性将使这不重要 - 如果你使用非就地基数排序，你需要额外的空间无论如何 .

  请注意，mergesort最终可以在没有递归的情况下实现，实际上这样做可以清楚地表明真正的线性内存访问模式 .

  回复于 2024-04-30T12:25:00+08:00
• 1

  好吧，这里's a simple implementation of an MSD radix sort for DNA. It'用D写，因为's the language that I use most and therefore am least likely to make silly mistakes in, but it could easily be translated to some other language. It'就位，但需要 2 * seq.length 通过数组 .

  void radixSort(string[] seqs, size_t base = 0) {
      if(seqs.length == 0)
          return;
  
      size_t TPos = seqs.length, APos = 0;
      size_t i = 0;
      while(i < TPos) {
          if(seqs[i][base] == 'A') {
               swap(seqs[i], seqs[APos++]);
               i++;
          }
          else if(seqs[i][base] == 'T') {
              swap(seqs[i], seqs[--TPos]);
          } else i++;
      }
  
      i = APos;
      size_t CPos = APos;
      while(i < TPos) {
          if(seqs[i][base] == 'C') {
              swap(seqs[i], seqs[CPos++]);
          }
          i++;
      }
      if(base < seqs[0].length - 1) {
          radixSort(seqs[0..APos], base + 1);
          radixSort(seqs[APos..CPos], base + 1);
          radixSort(seqs[CPos..TPos], base + 1);
          radixSort(seqs[TPos..seqs.length], base + 1);
     }
  }
  

  显然，这是一种特定的DNA，而不是一般的，但它应该是快速的 .

  编辑：

  我很好奇这段代码是否真的有用，所以我在等待自己的生物信息学代码运行的同时测试/调试了它 . 现在上面的版本实际上已经过测试和运行 . 对于每个5个碱基的1000万个序列，它比优化的introsort快3倍 .

  回复于 2024-04-30T12:25:00+08:00
• 1

  我从来没有见过一个就地基数排序，从基数的性质来看，我怀疑只要临时数组适合内存，它比不合适的排序要快得多 .

  Reason:

  排序在输入数组上进行线性读取，但所有写入几乎都是随机的 . 从某个N向上，这可归结为每次写入的高速缓存未命中 . 此缓存未命中会降低您的算法速度 . 如果它到位不会改变这种效果 .

  我知道这不会直接回答你的问题，但如果排序是一个瓶颈你可能想看看 near sorting 算法作为 preprocessing step （软堆上的维基页面可能会让你开始） .

  这可以提供非常好的缓存局部性提升 . 然后，教科书不合适的基数排序将表现得更好 . 写入仍然几乎是随机的，但至少它们将聚集在相同的内存块中，因此增加了缓存命中率 .

  我不知道它是否在实践中有效 .

  顺便说一句：如果你只处理DNA字符串：你可以将一个字符串压缩成两位，并打包你的数据非常多 . 这将使内存需求减少因子4而非naiive表示 . 寻址变得更加复杂，但是你的CPU的ALU有很多时间花在所有的上缓存未命中 .

  回复于 2024-04-30T12:25:00+08:00
• 2

  基于dsimcha 's code, I' ve实现了一个更通用的版本，它非常适合我们使用的框架（SeqAn） . 实际上，移植代码非常简单 . 之后我才发现实际上有关于这个主题的出版物 . 最棒的是：他们基本上和你们说的一样 . Andersson和Nilsson在Implementing Radixsort上发表的一篇论文绝对值得一读 . 如果您碰巧认识德语，请务必阅读David Weese的diploma thesis，其中他实现了通用子字符串索引 . 本文的大部分内容都致力于详细分析构建索引的成本，考虑二级内存和极大文件 . 他的工作成果实际上是在SeqAn中实现的，而不是在我需要它的那些部分 .

  只是为了好玩，这是我写的代码（我不认为任何不使用SeqAn的人会对它有任何用处） . 请注意，它仍然没有考虑更大的基数4.我希望这会对性能产生巨大影响，但遗憾的是我现在没有时间实现这一点 .

  对于短字符串，代码的执行速度是Introsort的两倍多 . 盈亏 balancer 点的长度约为12-13 . 字符串的类型（例如，它是否具有4个，5个或16个不同的值）相对不重要 . 从人类基因组的染色体2中分选> 6,000,000个DNA读数在我的PC上仅需2秒钟 . 只是为了记录，这很快！特别是考虑到我不使用SIMD或任何其他硬件加速 . 此外，valgrind告诉我，字符串赋值中的主要瓶颈是 operator new . 它被调用大约65,000,000次 - 每个字符串十次！这是一个死的赠品 swap 可以针对这些字符串进行优化：它不是复制，而是可以交换所有字符 . 我没有相信它会让人感到有点不同 . 而且，再说一遍，万一有人不听：基数大小对运行时几乎没有影响 - 这意味着我一定要尝试实施FryGuy，Stephan和EvilTeach提出的建议 .

  啊，顺便说一句： cache locality is a noticeable factor ：从1M字符串开始，运行时不再线性增加 . 但是，这可以很容易地修复：我对小子集（<= 20个字符串）使用插入排序 - 而不是随机黑客建议的mergesort . 显然，对于这样的小列表，这比mergesort表现更好（参见我链接的第一篇论文） .

  namespace seqan {
  
  template <typename It, typename F, typename T>
  inline void prescan(It front, It back, F op, T const& id) {
      using namespace std;
      if (front == back) return;
      typename iterator_traits<It>::value_type accu = *front;
      *front++ = id;
      for (; front != back; ++front) {
          swap(*front, accu);
          accu = op(accu, *front);
      }
  }
  
  template <typename TIter, typename TSize, unsigned int RADIX>
  inline void radix_permute(TIter front, TIter back, TSize (& bounds)[RADIX], TSize base) {
      for (TIter i = front; i != back; ++i)
          ++bounds[static_cast<unsigned int>((*i)[base])];
  
      TSize fronts[RADIX];
  
      std::copy(bounds, bounds + RADIX, fronts);
      prescan(fronts, fronts + RADIX, std::plus<TSize>(), 0);
      std::transform(bounds, bounds + RADIX, fronts, bounds, plus<TSize>());
  
      TSize active_base = 0;
  
      for (TIter i = front; i != back; ) {
          if (active_base == RADIX - 1)
              return;
          while (fronts[active_base] >= bounds[active_base])
              if (++active_base == RADIX - 1)
                  return;
          TSize current_base = static_cast<unsigned int>((*i)[base]);
          if (current_base <= active_base)
              ++i;
          else
              std::iter_swap(i, front + fronts[current_base]);
          ++fronts[current_base];
      }
  }
  
  template <typename TIter, typename TSize>
  inline void insertion_sort(TIter front, TIter back, TSize base) {
      typedef typename Value<TIter>::Type T;
      struct {
          TSize base, len;
          bool operator ()(T const& a, T const& b) {
              for (TSize i = base; i < len; ++i)
                  if (a[i] < b[i]) return true;
                  else if (a[i] > b[i]) return false;
              return false;
          }
      } cmp = { base, length(*front) }; // No closures yet. :-(
  
      for (TIter i = front + 1; i != back; ++i) {
          T value = *i;
          TIter j = i;
          for ( ; j != front && cmp(value, *(j - 1)); --j)
              *j = *(j - 1);
          if (j != i)
              *j = value;
      }
  }
  
  template <typename TIter, typename TSize, unsigned int RADIX>
  inline void radix(TIter top, TIter front, TIter back, TSize base, TSize (& parent_bounds)[RADIX], TSize next) {
      if (back - front > 20) {
          TSize bounds[RADIX] = { 0 };
          radix_permute(front, back, bounds, base);
  
          // Sort current bucket recursively by suffix.
          if (base < length(*front) - 1)
              radix(front, front, front + bounds[0], base + 1, bounds, static_cast<TSize>(0));
      }
      else if (back - front > 1)
          insertion_sort(front, back, base);
  
      // Sort next buckets on same level recursively.
      if (next == RADIX - 1) return;
      radix(top, top + parent_bounds[next], top + parent_bounds[next + 1], base, parent_bounds, next + 1);
  }
  
  template <typename TIter>
  inline void radix_sort(TIter front, TIter back) {
      typedef typename Container<TIter>::Type TStringSet;
      typedef typename Value<TStringSet>::Type TString;
      typedef typename Value<TString>::Type TChar;
      typedef typename Size<TStringSet>::Type TSize;
  
      TSize const RADIX = ValueSize<TChar>::VALUE;
      TSize bounds[RADIX];
  
      radix(front, front, back, static_cast<TSize>(0), bounds, RADIX - 1);
  }
  
  } // namespace seqan
  

  回复于 2024-04-30T12:25:00+08:00
• 20

  您当然可以通过以位为单位对序列进行编码来降低内存要求 . 您正在考虑排列，因此，对于长度为2，“ACGT”为16个状态或4位 . 对于长度3，这是64个状态，可以以6位编码 . 因此，对于序列中的每个字母，它看起来像2位，或者像你说的那样，对于16个字符，大约32位 .

  如果有办法减少有效“单词”的数量，则可以进一步压缩 .

  因此，对于长度为3的序列，可以创建64个桶，可能是uint32或uint64 . 将它们初始化为零 . 迭代你非常大的3个char序列列表，并按上面的方式对它们进行编码 . 使用此作为下标，并增加该桶 .
  重复此操作，直到处理完所有序列 .

  接下来，重新生成列表 .

  按顺序迭代64个桶，对于在该桶中找到的计数，生成该桶所表示的序列的许多实例 .
  当所有的桶都被迭代后，你就有了你的排序数组 .

  4的序列，加2位，所以会有256个桶 . 5的序列，加2位，所以会有1024个桶 .

  在某些时候，桶的数量将接近你的极限 . 如果从文件中读取序列，而不是将它们保留在内存中，则可以为存储桶提供更多内存 .

  我认为这比在原地进行排序要快，因为铲斗可能适合您的工作组 .

  这是一个显示该技术的黑客

  #include <iostream>
  #include <iomanip>
  
  #include <math.h>
  
  using namespace std;
  
  const int width = 3;
  const int bucketCount = exp(width * log(4)) + 1;
        int *bucket = NULL;
  
  const char charMap[4] = {'A', 'C', 'G', 'T'};
  
  void setup
  (
      void
  )
  {
      bucket = new int[bucketCount];
      memset(bucket, '\0', bucketCount * sizeof(bucket[0]));
  }
  
  void teardown
  (
      void
  )
  {
      delete[] bucket;
  }
  
  void show
  (
      int encoded
  )
  {
      int z;
      int y;
      int j;
      for (z = width - 1; z >= 0; z--)
      {
          int n = 1;
          for (y = 0; y < z; y++)
              n *= 4;
  
          j = encoded % n;
          encoded -= j;
          encoded /= n;
          cout << charMap[encoded];
          encoded = j;
      }
  
      cout << endl;
  }
  
  int main(void)
  {
      // Sort this sequence
      const char *testSequence = "CAGCCCAAAGGGTTTAGACTTGGTGCGCAGCAGTTAAGATTGTTT";
  
      size_t testSequenceLength = strlen(testSequence);
  
      setup();
  
  
      // load the sequences into the buckets
      size_t z;
      for (z = 0; z < testSequenceLength; z += width)
      {
          int encoding = 0;
  
          size_t y;
          for (y = 0; y < width; y++)
          {
              encoding *= 4;
  
              switch (*(testSequence + z + y))
              {
                  case 'A' : encoding += 0; break;
                  case 'C' : encoding += 1; break;
                  case 'G' : encoding += 2; break;
                  case 'T' : encoding += 3; break;
                  default  : abort();
              };
          }
  
          bucket[encoding]++;
      }
  
      /* show the sorted sequences */ 
      for (z = 0; z < bucketCount; z++)
      {
          while (bucket[z] > 0)
          {
              show(z);
              bucket[z]--;
          }
      }
  
      teardown();
  
      return 0;
  }
  

  回复于 2024-04-30T12:25:00+08:00
• 3

  如果您的数据集太大，那么我认为基于磁盘的缓冲区方法最好：

  sort(List<string> elements, int prefix)
      if (elements.Count < THRESHOLD)
           return InMemoryRadixSort(elements, prefix)
      else
           return DiskBackedRadixSort(elements, prefix)
  
  DiskBackedRadixSort(elements, prefix)
      DiskBackedBuffer<string>[] buckets
      foreach (element in elements)
          buckets[element.MSB(prefix)].Add(element);
  
      List<string> ret
      foreach (bucket in buckets)
          ret.Add(sort(bucket, prefix + 1))
  
      return ret
  

  我还会尝试分组到更多数量的存储桶，例如，如果您的字符串是：

  GATTACA
  

  第一个MSB调用将返回GATT的桶（总共256个桶），这样就可以减少基于磁盘的缓冲区的分支 . 这可能会也可能不会提高性能，所以请尝试一下 .

  回复于 2024-04-30T12:25:00+08:00
• 3

  我打算走出去，建议你切换到堆/ heapsort实现 . 这个建议有一些假设：

  • 您可以控制数据的读取

  • 只要'start'对其进行排序，您就可以对排序后的数据执行有意义的操作 .

  堆/堆排序的优点在于，您可以在读取数据时构建堆，并且可以在构建堆时开始获取结果 .

  我们退一步吧 . 如果你是如此幸运，你可以异步读取数据（也就是你可以发布某种读取请求并在某些数据准备就绪时得到通知），然后您可以在等待下一个数据块进入时构建一大块堆 - 甚至可以从磁盘中获取 . 通常，这种方法可以将一半分类的大部分成本埋在获取数据所花费的时间之后 .

  读取数据后，第一个元素已经可用 . 根据您发送数据的位置，这可能很棒 . 如果要将其发送到另一个异步读取器或某个并行“事件”模型或UI，则可以随时发送块和块 .

  也就是说 - 如果您无法控制数据的读取方式，并且它是同步读取的，并且在完全写出之前您没有使用已排序的数据 - 请忽略所有这些 . :(

  查看维基百科文章：

  • Heapsort

  • Binary heap

  回复于 2024-04-30T12:25:00+08:00
• 3

  在性能方面，您可能希望查看更通用的字符串比较排序算法 .

  目前你最终触及每个字符串的每个元素，但你可以做得更好！

  特别是，burst sort非常适合这种情况 . 作为奖励，因为burstsort基于尝试，它对于DNA / RNA中使用的小字母大小非常有效，因为您不需要构建任何类型的三元搜索节点，散列或其他节点压缩方案 . 实施 . 尝试对于类似后缀数组的最终目标也许有用 .

  一个体面的burstsort通用实现可以在http://sourceforge.net/projects/burstsort/的源伪造上获得 - 但它不是就地的 .

  出于比较目的，http://www.cs.mu.oz.au/~rsinha/papers/SinhaRingZobel-2006.pdf基准测试中的C-burstsort实现比一些典型工作负载的快速排序和基数排序快4-5倍 .

  回复于 2024-04-30T12:25:00+08:00
• 6

  您可以尝试使用trie . 对数据进行排序只需迭代数据集并插入即可;结构是自然排序的，您可以将其视为类似于B树（除了进行比较，您总是使用指针间接） .

  缓存行为将有利于所有内部节点，因此您可能无法改进;但你也可以调整你的trie的分支因子（确保每个节点都适合单个缓存行，分配类似于堆的trie节点，作为表示级别顺序遍历的连续数组） . 由于尝试也是数字结构（对于长度为k的元素，O（k）插入/查找/删除），您应该具有基数排序的竞争性能 .

  回复于 2024-04-30T12:25:00+08:00
• 4

  我会burstsort字符串的打包位表示 . 声称Burstsort比基数排序具有更好的局部性，通过爆发尝试代替经典尝试来保持额外的空间使用 . 原始论文有测量 .

  回复于 2024-04-30T12:25:00+08:00
• 19

  看起来你在这里描述的是've solved the problem, but for the record, it appears that one version of a workable in-place radix sort is the 869032 . It'：Engineering Radix Sort . 一般的想法是对每个字符进行2次传递 - 首先计算每个字符的数量，这样就可以将输入数组细分为二进制数 . 然后再次检查，将每个元素交换到正确的bin中 . 现在递归地对下一个字符位置的每个bin进行排序 .

  回复于 2024-04-30T12:25:00+08:00
• 8

  首先，考虑一下您的问题的编码 . 摆脱了字符串，用二进制表示替换它们 . 使用第一个字节表示长度编码 . 或者，在四字节边界处使用固定长度表示 . 然后基数排序变得更容易 . 对于基数排序，最重要的是不要在内循环的热点处进行异常处理 .

  好吧，我想了解更多有关4-nary问题的内容 . 你想要一个像Judy tree这样的解决方案 . 下一个解决方案可以处理变长字符串;对于固定长度，只需删除长度位，这实际上使它更容易 .

  分配16个指针的块 . 指针的最低位可以重复使用，因为您的块将始终对齐 . 您可能需要一个特殊的存储分配器（将大存储分成更小的块） . 有许多不同类型的块：

  • 使用7个长度位的可变长度字符串进行编码 . 当他们填满时，你将它们替换为：

  • 位置编码接下来的两个字符，你有16个指向下一个块的指针，结束于：

  • 字符串最后三个字符的位图编码 .

  对于每种块，您需要在LSB中存储不同的信息 . 由于你有可变长度的字符串，你也需要存储字符串结尾，最后一种块只能用于最长的字符串 . 当你深入到结构中时，7个长度的位应该被更少的替换 .

  这为您提供了排序字符串的相当快速且非常高效的内存存储 . 它的行为有点像trie . 要使其正常工作，请确保构建足够的单元测试 . 您想要覆盖所有块转换 . 您只想从第二种块开始 .

  为了获得更高的性能，您可能希望添加不同的块类型和更大的块大小 . 如果块总是大小相同且足够大，则可以使用更少的位作为指针 . 块大小为16指针时，您已经在32位地址空间中释放了一个字节 . 查看Judy树文档，了解有趣的块类型 . 基本上，您为空间（和运行时）权衡添加代码和工程时间

  您可能希望从前四个字符的256宽直接基数开始 . 这提供了一个合适的空间/时间权衡 . 在这个实现中，与简单的trie相比，你获得的内存开销要少得多;它大约小三倍（我没有测量过） . 如果常量足够低，则O（n）没有问题，正如您在与O（n log n）快速排序进行比较时所注意到的那样 .

  你有兴趣处理双打吗？有了短序列，就会有 . 调整块以处理计数是棘手的，但它可以非常节省空间 .

  回复于 2024-04-30T12:25:00+08:00