Это длинный текст. Пожалуйста, потерпите меня. Вопрос сводится к следующему: существует ли работоспособный алгоритм сортировки по основанию ?

предварительный

У меня есть огромное количество маленьких строк фиксированной длины, которые используют только буквы «A», «C», «G» и «T» (да, вы уже догадались: ДНК ), которые я хочу отсортировать.

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

подробности

Я проверил это предположение с очень наивной реализацией и для относительно небольших входных данных (порядка 10000) это было правдой (ну, по крайней мере, более чем в два раза быстрее). Тем не менее, время выполнения ужасно ухудшается, когда размер проблемы становится больше ( N > 5 000 000).

Причина очевидна: радикальная сортировка требует копирования всех данных (на самом деле, не раз в моей наивной реализации). Это означает, что я поместил ~ 4 ГиБ в мою основную память, что, очевидно, снижает производительность. Даже если это не так, я не могу позволить себе использовать столько памяти, поскольку размеры проблем на самом деле становятся еще больше.

Случаи использования

В идеале этот алгоритм должен работать с любой длиной строки от 2 до 100, для ДНК, а также для ДНК5 (которая допускает дополнительный подстановочный знак «N»), или даже для ДНК с кодами неоднозначности IUPAC (в результате получается 16 различных значений). Тем не менее, я понимаю, что все эти случаи не могут быть покрыты, поэтому я доволен любым улучшением скорости, которое я получаю. Код может динамически решать, какой алгоритм отправить.

Исследовательская работа

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

В частности, есть следующие вещи.

Во-первых, алгоритм содержит несколько ошибок и оставляет много необъяснимых. В частности, он не детализирует рекурсивный вызов (я просто предполагаю, что он увеличивает или уменьшает некоторый указатель для вычисления текущих значений сдвига и маски). Также он использует функции dest_groupи dest_addressне дает определения. Я не вижу, как реализовать это эффективно (то есть в O (1); по крайней мере dest_address, не тривиально).

Наконец, что не менее важно, алгоритм достигает на месте путем замены индексов массива элементами внутри входного массива. Это очевидно работает только на числовых массивах. Мне нужно использовать его на строки. Конечно, я мог бы просто прокрутить строгую типизацию и продолжить, предполагая, что память будет терпеть мое хранение индекса там, где он не принадлежит. Но это работает только до тех пор, пока я могу втиснуть свои строки в 32-битную память (при условии, что 32-битные целые числа). Это всего 16 символов (на данный момент давайте игнорируем 16> log (5 000 000)).

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

Напомним : есть ли надежда найти работающую эталонную реализацию или, по крайней мере, хороший псевдокод / ​​описание работающей радикальной сортировки на месте, которая работает со строками ДНК?

Ну, вот простая реализация радикальной сортировки 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 раза быстрее, чем оптимизированный интросорт.

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

Причина:

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

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

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

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

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

Конечно, вы можете отбросить требования к памяти, кодируя последовательность в битах. Вы смотрите на перестановки, поэтому для длины 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;
}

Если ваш набор данных такой большой, то я думаю, что дисковый буферный подход будет лучше:

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

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

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

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

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

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

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

Смотрите статьи в Википедии:

• Пирамидальная сортировка
• Двоичная куча

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

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

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

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

Достойная универсальная реализация burstsort доступна в source forge по адресу http://sourceforge.net/projects/burstsort/, но она не на месте.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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