Сортировка По Основанию На Месте


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


предварительно

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

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

подробности

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

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

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

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

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

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

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

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

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

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

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

15 185

15 ответов:

Ну, вот простая реализация сортировки 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);
   }
}

Edit:

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

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

причина:

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

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

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

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

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

на основе кода dsimcha, я реализовал более общую версию, которая хорошо вписывается в рамки, которые мы используем (SeqAn). На самом деле, перенос кода был совершенно простым. Только потом я обнаружил, что там are собственно публикации, касающиеся именно этой темы. Самое замечательное: они в основном говорят то же самое, что и вы, ребята. Статья Андерссона и Нильссона о Реализация Radixsort определенно стоит прочитать. Если вы случайно знаете немецкий язык, обязательно также читайте Дэвида Уиза дипломная работа где он реализует общий индекс подстроки. Большая часть диссертации посвящена детальному анализу стоимости построения индекса с учетом вторичной памяти и чрезвычайно больших файлов. Результаты его работы фактически были реализованы в Секане, только не в тех частях, где мне это было нужно.

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

код выполняет более чем в два раза быстрее, чем Introsort для коротких строк. Точка безубыточности находится на длине около 12-13. Тип строки (например, имеет ли она 4, 5 или 16 различных значений) сравнительно неважен. Сортировка > 6,000,000 ДНК считывает из хромосомы 2 генома человека принимает чуть более 2 секунд на моем компьютере. Просто для протокола, это быстро! Особенно учитывая, что я не использую SIMD или любое другое аппаратное ускорение. Кроме того, отчет показывает, что основным узким местом является operator new в строке задания. Он вызывается около 65 000 000 раз-десять раз для каждой строки! Это мертвый поддавки, что swap может быть оптимизирован под эти строки: вместо того, чтобы делать копии, он может просто поменять все знаки. Я не пробовал этого, но я убежден что это будет иметь чертовски большое значение. И еще раз, на случай, если кто-то не слушает: размер радикса почти не влияет на время выполнения - что означает, что я должен наверняка попробуйте реализовать предложение, сделанное FryGuy, Stephan и EvilTeach.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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