Что ж, вот простая реализация сортировки по основанию MSD для ДНК. Он написан на D, потому что это язык, который я использую чаще всего и, следовательно, с меньшей вероятностью сделаю глупые ошибки, но его можно легко перевести на другой язык. Он на месте, но требует, чтобы 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);
   }
}

Очевидно, что это специфично для ДНК, а не является общим, но это должно быть быстро.

Редактировать:

Мне стало любопытно, действительно ли этот код работает, поэтому я протестировал / отладил его, ожидая запуска моего собственного биоинформатического кода. Приведенная выше версия сейчас фактически протестирована и работает. Для 10 миллионов последовательностей по 5 баз в каждой это примерно в 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 сегментов), таким образом вы сделаете меньше ветвей дискового буфера. Это может улучшить производительность, а может и не улучшить, поэтому поэкспериментируйте.

Я никогда не видел сортировки по основанию на месте, и, судя по природе сортировки по основанию, я сомневаюсь, что она намного быстрее, чем сортировка вне места, если временный массив помещается в память.

Причина:

Сортировка выполняет линейное чтение входного массива, но все записи будут почти случайными. Начиная с определенного N и выше, это сводится к пропуску кеш-памяти при каждой записи. Этот промах в кеше замедляет работу вашего алгоритма. Если он на месте или нет, этот эффект не изменится.

Я знаю, что это не ответит на ваш вопрос напрямую, но если сортировка является узким местом, вы можете взглянуть на алгоритмы рядом с сортировкой как на шаг предварительной обработки (wiki-страница в soft-heap может помочь вам начать).

Это могло бы дать очень хороший прирост локальности кеша. Тогда лучше будет справиться сортировка по неуместному основанию из учебника. Записи по-прежнему будут почти случайными, но, по крайней мере, они будут сгруппированы вокруг одних и тех же блоков памяти и, как таковые, увеличат коэффициент попадания в кеш.

Я понятия не имею, сработает ли это на практике.

Кстати: если вы имеете дело только со строками ДНК: вы можете сжать символ на два бита и довольно много упаковать данные. Это сократит потребность в памяти в четыре раза по сравнению с обычным представлением. Адресация становится более сложной, но у ALU вашего процессора в любом случае есть много времени, чтобы потратить его во время всех промахов кеша.

Я собираюсь рискнуть и предложить вам перейти на реализацию heap / heapsort. Это предложение содержит некоторые предположения:

1. Вы контролируете чтение данных
2. Вы можете сделать что-нибудь значимое с отсортированными данными, как только вы «начнете» их сортировать.

Прелесть сортировки кучи / кучи заключается в том, что вы можете создавать кучу во время чтения данных, и вы можете начать получать результаты в тот момент, когда вы создали кучу.

Вернемся назад. Если вам так повезло, что вы можете читать данные асинхронно (то есть вы можете отправить какой-то запрос на чтение и получить уведомление, когда некоторые данные будут готовы), а затем вы можете построить кусок кучи, пока вы ждете следующий блок данных - даже с диска. Часто этот подход может похоронить большую часть стоимости половины вашей сортировки за время, потраченное на получение данных.

Как только вы прочитали данные, первый элемент уже доступен. В зависимости от того, куда вы отправляете данные, это может быть здорово. Если вы отправляете его другому асинхронному считывателю, или какой-либо параллельной модели «событий» или пользовательскому интерфейсу, вы можете отправлять фрагменты и фрагменты по мере продвижения.

Тем не менее, если у вас нет контроля над тем, как данные читаются, и они читаются синхронно, и вы не можете использовать отсортированные данные, пока они не будут полностью записаны, игнорируйте все это. :(

См. Статьи в Википедии:

• Heapsort
• Двоичная куча

Сначала подумайте о коде вашей проблемы. Избавьтесь от строк, замените их двоичным представлением. Используйте первый байт, чтобы указать длину + кодировку. В качестве альтернативы используйте представление фиксированной длины на четырехбайтовой границе. Тогда сортировка по основанию системы счисления становится намного проще. Для поразрядной сортировки самое важное - не иметь обработки исключений в горячей точке внутреннего цикла.

Хорошо, я подумал немного больше о проблеме с четырьмя числами. Для этого вам нужно решение, подобное Джуди дерево. Следующее решение может обрабатывать строки переменной длины; для фиксированной длины просто удалите биты длины, что на самом деле упрощает задачу.

Выделяем блоки по 16 указателей. Младший бит указателей можно использовать повторно, так как ваши блоки всегда будут выровнены. Вам может понадобиться специальный распределитель памяти для него (разбиение большого хранилища на более мелкие блоки). Есть несколько видов блоков:

• Кодирование с 7 битами длины строк переменной длины. По мере их заполнения вы заменяете их:
• Позиция кодирует следующие два символа, у вас есть 16 указателей на следующие блоки, заканчивающиеся на:
• Битовая кодировка последних трех символов строки.

Для каждого типа блока вам необходимо хранить различную информацию в младших разрядах. Поскольку у вас есть строки переменной длины, вам также необходимо хранить конец строки, а последний тип блока может использоваться только для самых длинных строк. 7 битов длины следует заменять на меньшее по мере того, как вы углубляетесь в структуру.

Это дает вам достаточно быстрое и очень эффективное по памяти хранение отсортированных строк. Он будет вести себя как три. Чтобы это работало, обязательно создайте достаточно модульных тестов. Вам нужно охватить все переходы блоков. Вы хотите начать только со вторым типом блоков.

Для еще большей производительности вы можете добавить различные типы блоков и больший размер блока. Если блоки всегда одного и того же размера и достаточно велики, вы можете использовать еще меньше бит для указателей. При размере блока в 16 указателей у вас уже есть свободный байт в 32-битном адресном пространстве. Взгляните на документацию по дереву Джуди, чтобы найти интересные типы блоков. По сути, вы добавляете код и время разработки в обмен на пространство (и время выполнения).

Вероятно, вы захотите начать с 256-кратного прямого основания для первых четырех символов. Это обеспечивает достойный компромисс между пространством и временем. В этой реализации вы получаете гораздо меньше накладных расходов на память, чем с простым деревом; он примерно в три раза меньше (не измерял). O (n) не проблема, если константа достаточно мала, как вы заметили при сравнении с быстрой сортировкой O (n log n).

Вы заинтересованы в работе с дублями? С короткими последовательностями они будут. Адаптировать блоки для обработки подсчетов сложно, но это может быть очень компактно.

Вы, конечно, можете отказаться от требований к памяти, кодируя последовательность в битах. Вы смотрите на перестановки, поэтому для длины 2 с «ACGT» это 16 состояний или 4 бита. Для длины 3 это 64 состояния, которые можно закодировать в 6 бит. Таким образом, это выглядит как 2 бита для каждой буквы в последовательности или примерно 32 бита для 16 символов, как вы сказали.

Если есть способ уменьшить количество допустимых «слов», возможно дальнейшее сжатие.

Таким образом, для последовательностей длиной 3 можно создать 64 сегмента, возможно, размером uint32 или uint64. Инициализируйте их обнулением. Просмотрите ваш очень большой список из 3 последовательностей символов и закодируйте их, как указано выше. Используйте это как нижний индекс и увеличивайте этот сегмент. Повторяйте это, пока все ваши последовательности не будут обработаны.

Затем обновите свой список.

Итерируйте по 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;
}

Вы можете попробовать использовать три. Сортировка данных - это просто итерация по набору данных и его вставка; структура сортируется естественным образом, и вы можете думать о ней как о B-дереве (за исключением того, что вместо сравнения вы всегда используете косвенные ссылки указателя).

Поведение кэширования будет благоприятствовать всем внутренним узлам, поэтому вы, вероятно, не улучшите его; но вы также можете поиграть с коэффициентом ветвления вашего дерева (убедитесь, что каждый узел помещается в одну строку кэша, выделите узлы дерева, аналогичные куче, в виде непрерывного массива, который представляет обход порядка уровней). Поскольку попытки также являются цифровыми структурами (O (k) insert / find / delete для элементов длины k), у вас должна быть конкурентоспособная производительность для сортировки по основанию.

Сортировка по основанию MSB в dsimcha выглядит неплохо, но Нильс ближе к сути проблемы обнаруживает, что локальность кеша убивает вас при больших размерах проблем.

Предлагаю очень простой подход:

1. Эмпирически оцените m наибольшего размера, для которого эффективна сортировка по основанию.
2. Считайте блоки элементов m за один раз, отсортируйте их по основанию и запишите их (в буфер памяти, если у вас достаточно памяти, но в противном случае в файл), пока вы не исчерпаете свой ввод.
3. Сортировка слиянием результирующих отсортированных блоков.

Сортировка слиянием - это наиболее удобный для кеша алгоритм сортировки, о котором я знаю: «Прочтите следующий элемент из массива A или B, затем запишите элемент в выходной буфер». Он эффективно работает на ленточные накопители. Для сортировки элементов 2n требуется пространство n, но я держу пари, что значительно улучшенная локализация кеша, которую вы увидите, сделает это неважным - и если вы использовали сортировку без места, вам нужно это дополнительное пространство в любом случае.

Наконец, обратите внимание, что сортировка слиянием может быть реализована без рекурсии, и, по сути, это проясняет истинный линейный шаблон доступа к памяти.

Похоже, вы решили проблему, но для протокола, похоже, что одной из версий работоспособной сортировки по основанию на месте является "Сортировка по американскому флагу". Это описано здесь: Инженерная сортировка по основанию Radix. Общая идея состоит в том, чтобы выполнить 2 прохода для каждого символа - сначала посчитайте, сколько из них у вас есть, чтобы вы могли разделить входной массив на ячейки. Затем пройдите еще раз, переставляя каждый элемент в правильную корзину. Теперь рекурсивно отсортируйте каждую ячейку по позиции следующего символа.

Я бы взрывная сортировка представлял упакованные биты строк. Утверждается, что Burstsort имеет гораздо лучшую локальность, чем сортировка по основанию, уменьшая использование дополнительного пространства с помощью пакетных попыток вместо классических попыток. На оригинальной бумаге есть размеры.

Radix-Sort не учитывает кеширование и не является самым быстрым алгоритмом сортировки для больших наборов. Вы можете посмотреть:

• ti7qsort. ti7qsort - самая быстрая сортировка для целых чисел (может использоваться для строк небольшого фиксированного размера).
• Встроенный QSORT
• Сортировка строк

Вы также можете использовать сжатие и кодировать каждую букву вашей ДНК в 2 бита перед сохранением в массиве сортировки.

С точки зрения производительности вы можете захотеть взглянуть на более общие алгоритмы сортировки сравнения строк.

В настоящее время вы касаетесь каждого элемента каждой струны, но вы можете добиться большего!

В частности, для этого случая очень хорошо подходит пакетная сортировка. В качестве бонуса, поскольку пакетная сортировка основана на попытках, она работает смехотворно хорошо для небольших размеров алфавита, используемых в ДНК / РНК, поскольку вам не нужно встраивать какой-либо тройной узел поиска, хэш или другую схему сжатия узла дерева. три реализации. Попытки также могут быть полезны для вашей конечной цели, подобной массиву суффиксов.

Приличная реализация пакетной сортировки общего назначения доступна в кузнице исходного кода по адресу http://sourceforge.net/projects/burstsort/, но ее нет на месте.

В целях сравнения, реализация C-burstsort, охваченная тестами http://www.cs.mu.oz.au/~rsinha/papers/SinhaRingZobel-2006.pdf, в 4-5 раз быстрее, чем сортировка по принципу быстрой сортировки и сортировки по основанию для некоторых типичных рабочих нагрузок.

Вы захотите взглянуть на Крупномасштабная обработка последовательности генома от Drs. Касахара и Моришита.

Строки, состоящие из четырех нуклеотидных букв A, C, G и T, могут быть специально закодированы в целые числа для более быстрой обработки много. Радиксная сортировка - один из многих алгоритмов, обсуждаемых в книге; вы сможете адаптировать принятый ответ на этот вопрос и увидеть значительное улучшение производительности.

«Сортировка Radix без лишнего места» - это документ, посвященный вашей проблеме.

Хотя принятый ответ полностью отвечает описанию проблемы, я дошел до этого места, тщетно пытаясь найти алгоритм для разделения встроенного массива на N частей. Я сам написал одну, вот она.

Предупреждение: это нестабильный алгоритм разбиения, поэтому для многоуровневого разбиения необходимо переразбивать каждый результирующий раздел, а не весь массив. Преимущество в том, что он встроенный.

Способ, которым это помогает с поставленным вопросом, заключается в том, что вы можете многократно разбивать встроенные разделы на основе буквы строки, а затем сортировать разделы, когда они достаточно малы, с помощью алгоритма по вашему выбору.

  function partitionInPlace(input, partitionFunction, numPartitions, startIndex=0, endIndex=-1) {
    if (endIndex===-1) endIndex=input.length;
    const starts = Array.from({ length: numPartitions + 1 }, () => 0);
    for (let i = startIndex; i < endIndex; i++) {
      const val = input[i];
      const partByte = partitionFunction(val);
      starts[partByte]++;
    }
    let prev = startIndex;
    for (let i = 0; i < numPartitions; i++) {
      const p = prev;
      prev += starts[i];
      starts[i] = p;
    }
    const indexes = [...starts];
    starts[numPartitions] = prev;
  
    let bucket = 0;
    while (bucket < numPartitions) {
      const start = starts[bucket];
      const end = starts[bucket + 1];
      if (end - start < 1) {
        bucket++;
        continue;
      }
      let index = indexes[bucket];
      if (index === end) {
        bucket++;
        continue;
      }
  
      let val = input[index];
      let destBucket = partitionFunction(val);
      if (destBucket === bucket) {
        indexes[bucket] = index + 1;
        continue;
      }
  
      let dest;
      do {
        dest = indexes[destBucket] - 1;
        let destVal;
        let destValBucket = destBucket;
        while (destValBucket === destBucket) {
          dest++;
          destVal = input[dest];
          destValBucket = partitionFunction(destVal);
        }
  
        input[dest] = val;
        indexes[destBucket] = dest + 1;
  
        val = destVal;
        destBucket = destValBucket;
      } while (dest !== index)
    }
    return starts;
  }