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

https://stackoverflow.com/questions/463105

19-08-2019
|

Вопрос

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

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

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

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

Подробные сведения

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

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

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

В идеале, этот алгоритм должен работать с любой длиной строки от 2 до 100, как для ДНК, так и для DNA5 (что допускает дополнительный подстановочный знак “N”), или даже с ДНК с ИЮПАК коды двусмысленности (в результате получается 16 различных значений).Тем не менее, я понимаю, что все эти случаи не могут быть охвачены, поэтому я доволен любым улучшением скорости, которое я получаю.Код может динамически решать, к какому алгоритму обращаться.

Исследования

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

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

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

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

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

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

Решение

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

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

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

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

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

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

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

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

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

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

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

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

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

Вы наверняка захотите взглянуть на Крупномасштабная Обработка Последовательностей генома автор : доктор Дж .Касахара и Моришита.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Лицензировано под: CC-BY-SA с атрибуция

Не связан с StackOverflow