Сортировка по месту

Question

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

---

Предварительные

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

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

Подробнее

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

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

Варианты использования

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

Research

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

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

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

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

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

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

Answer 1

Хорошо, вот простая реализация радикальной сортировки 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);
   }
}

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

Edit:

Мне стало любопытно, работает ли этот код на самом деле, поэтому я протестировал / отладил его, ожидая запуска собственного кода биоинформатики. Версия выше теперь на самом деле протестирована и работает. Для 10 миллионов последовательностей по 5 оснований это примерно в 3 раза быстрее, чем оптимизированный интросорт.

Answer 2

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

Причина:

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

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

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

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

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

Answer 3

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

Answer 4

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

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

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

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

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

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

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

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

Answer 5

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

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

Answer 6

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

Answer 7

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

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

Answer 8

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

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

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

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

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

Answer 9

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

Answer 10

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

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