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

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

Причина:

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

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

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

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

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

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

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

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

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

 namespace seqan { template  inline void prescan(It front, It back, F op, T const& id) { using namespace std; if (front == back) return; typename iterator_traits::value_type accu = *front; *front++ = id; for (; front != back; ++front) { swap(*front, accu); accu = op(accu, *front); } } template  inline void radix_permute(TIter front, TIter back, TSize (& bounds)[RADIX], TSize base) { for (TIter i = front; i != back; ++i) ++bounds[static_cast((*i)[base])]; TSize fronts[RADIX]; std::copy(bounds, bounds + RADIX, fronts); prescan(fronts, fronts + RADIX, std::plus(), 0); std::transform(bounds, bounds + RADIX, fronts, bounds, plus()); 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((*i)[base]); if (current_base <= active_base) ++i; else std::iter_swap(i, front + fronts[current_base]); ++fronts[current_base]; } } template  inline void insertion_sort(TIter front, TIter back, TSize base) { typedef typename Value::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  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(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  inline void radix_sort(TIter front, TIter back) { typedef typename Container::Type TStringSet; typedef typename Value::Type TString; typedef typename Value::Type TChar; typedef typename Size::Type TSize; TSize const RADIX = ValueSize::VALUE; TSize bounds[RADIX]; radix(front, front, back, static_cast(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  #include  #include  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; } 

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

 sort(List elements, int prefix) if (elements.Count < THRESHOLD) return InMemoryRadixSort(elements, prefix) else return DiskBackedRadixSort(elements, prefix) DiskBackedRadixSort(elements, prefix) DiskBackedBuffer[] buckets foreach (element in elements) buckets[element.MSB(prefix)].Add(element); List ret foreach (bucket in buckets) ret.Add(sort(bucket, prefix + 1)) return ret 

Я бы также экспериментировал с группировкой в ​​большее количество ведер, например, если ваша строка была:

 GATTACA 

первый вызов MSB возвратит ведро для GATT (256 полных ковшей), таким образом вы создадите меньше ветвей буфера на основе диска. Это может или не может улучшить производительность, поэтому экспериментируйте с ним.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Для большей производительности вы можете добавить различные типы блоков и больший размер блока. Если блоки имеют одинаковый размер и достаточно большой, вы можете использовать меньшее количество бит для указателей. With a block size of 16 pointers, you already have a byte free in a 32-bit address space. Take a look at the Judy tree documentation for interesting block types. Basically, you add code and engineering time for a space (and runtime) trade-off

You probably want to start with a 256 wide direct radix for the first four characters. That provides a decent space/time tradeoff. In this implementation, you get much less memory overhead than with a simple trie; it is approximately three times smaller (I haven’t measured). O(n) is no problem if the constant is low enough, as you noticed when comparing with the O(n log n) quicksort.

Are you interested in handling doubles? With short sequences, there are going to be. Adapting the blocks to handle counts is tricky, but it can be very space-efficient.