algorithm - pseudocódigo - En Lugar Radix Sort
radix sort pdf (15)
" Clasificación de radix sin espacio adicional " es un documento que aborda su problema.
Este es un texto largo. Por favor, tenga paciencia conmigo. La pregunta es: ¿hay un algoritmo de ordenación de raíz in situ viable ?
Preliminar
Tengo un gran número de cadenas pequeñas de longitud fija que solo utilizan las letras "A", "C", "G" y "T" (sí, lo has adivinado: DNA ) que quiero ordenar.
Por el momento, uso std::sort que usa introsort en todas las implementaciones comunes de STL . Esto funciona bastante bien Sin embargo, estoy convencido de que la clasificación de radix se ajusta perfectamente a mi problema y debería funcionar mucho mejor en la práctica.
Detalles
He probado esta suposición con una implementación muy ingenua y para entradas relativamente pequeñas (del orden de 10.000) esto fue cierto (bueno, al menos más del doble de rápido). Sin embargo, el tiempo de ejecución se degrada abismalmente cuando el tamaño del problema aumenta ( N > 5,000,000).
La razón es obvia: el ordenamiento de radix requiere copiar los datos completos (más de una vez en mi implementación ingenua, en realidad). Esto significa que he puesto ~ 4 GiB en mi memoria principal, lo que obviamente mata el rendimiento. Incluso si no fuera así, no puedo permitirme usar tanta memoria, ya que los tamaños de los problemas en realidad son aún mayores.
Casos de uso
Idealmente, este algoritmo debería funcionar con cualquier longitud de cadena entre 2 y 100, tanto para ADN como para ADN5 (que permite un carácter comodín adicional "N"), o incluso ADN con códigos de ambigüedad IUPAC (resultando en 16 valores distintos). Sin embargo, me doy cuenta de que no se pueden cubrir todos estos casos, así que estoy contento con cualquier mejora de velocidad que obtenga. El código puede decidir dinámicamente a qué algoritmo enviar.
Investigación
Desafortunadamente, el artículo de Wikipedia sobre clasificación de radix es inútil. La sección sobre una variante in situ es basura completa. La sección NIST-DADS sobre clasificación de radix es casi inexistente. Hay un documento de sonido prometedor llamado Efficient Adaptive in situ Radix Sorting que describe el algoritmo "MSL". Lamentablemente, este documento también es decepcionante.
En particular, hay las siguientes cosas.
Primero, el algoritmo contiene varios errores y deja mucho sin explicación. En particular, no detalla la llamada de recursión (simplemente supongo que incrementa o reduce algún puntero para calcular el cambio actual y los valores de máscara). Además, usa las funciones dest_group y dest_address sin dar definiciones. No veo cómo implementarlos de manera eficiente (es decir, en O (1); al menos dest_address no es trivial).
Por último, pero no menos importante, el algoritmo logra in-place-ness al intercambiar índices de matriz con elementos dentro de la matriz de entrada. Obviamente, esto solo funciona en matrices numéricas. Necesito usarlo en cadenas. Por supuesto, podría atornillar mecanografía fuerte y seguir adelante asumiendo que la memoria tolerará que guarde un índice al que no pertenece. Pero esto solo funciona siempre que pueda comprimir mis cadenas en 32 bits de memoria (suponiendo enteros de 32 bits). Eso es solo 16 caracteres (ignoremos por el momento que 16> log (5,000,000)).
Otro documento de uno de los autores no ofrece ninguna descripción precisa, pero le da al tiempo de ejecución de MSL como sub-lineal, que es totalmente erróneo.
Para recapitular : ¿hay alguna esperanza de encontrar una implementación de referencia de trabajo o, al menos, un buen pseudocódigo / descripción de una ordenación de radix in situ que funcione en cadenas de ADN?
Basado en el código de dsimcha, he implementado una versión más genérica que encaja bien en el marco que usamos (SeqAn). En realidad, portar el código fue completamente sencillo. Solo después encontré que en realidad hay publicaciones sobre este tema. Lo bueno es que básicamente dicen lo mismo que ustedes. Un documento de Andersson y Nilsson sobre Implementando Radixsort definitivamente vale la pena leerlo. Si usted sabe alemán, asegúrese de leer también la tesis de diploma de David Weese donde implementa un índice de subcadena genérico. La mayor parte de la tesis está dedicada a un análisis detallado del costo de construcción del índice, considerando la memoria secundaria y archivos extremadamente grandes. Los resultados de su trabajo se han implementado realmente en SeqAn, pero no en las partes donde lo necesitaba.
Solo por diversión, aquí está el código que he escrito (no creo que nadie que no use SeqAn lo use). Tenga en cuenta que aún no considera radixes mayores 4. Espero que esto tenga un gran impacto en el rendimiento, pero desafortunadamente simplemente no tengo tiempo ahora para implementar esto.
El código realiza más del doble de rápido que Introsort para cadenas cortas. El punto de equilibrio está en una longitud de aproximadamente 12-13. El tipo de cadena (por ejemplo, si tiene 4, 5 o 16 valores diferentes) es comparativamente poco importante. Clasificar> 6,000,000 lecturas de ADN del cromosoma 2 del genoma humano lleva poco más de 2 segundos en mi PC. Solo para el registro, ¡eso es rápido ! Sobre todo teniendo en cuenta que no utilizo SIMD ni ninguna otra aceleración de hardware. Además, valgrind me muestra que el cuello de botella principal es operator new en las asignaciones de cadenas. Se llama alrededor de 65,000,000 de veces, ¡diez veces por cada cadena! Este es un claro ejemplo de que el swap podría optimizarse para estas cadenas: en lugar de hacer copias, podría simplemente intercambiar todos los caracteres. No intenté esto, pero estoy convencido de que sería una gran diferencia. Y, solo para decirlo nuevamente, en caso de que alguien no estuviera escuchando: el tamaño del radix casi no tiene influencia en el tiempo de ejecución , lo que significa que definitivamente debería intentar implementar la sugerencia hecha por FryGuy, Stephan y EvilTeach.
Ah, sí, por cierto: la localidad de caché es un factor notable : a partir de cadenas de 1M, el tiempo de ejecución ya no aumenta linealmente. Sin embargo, esto podría solucionarse con bastante facilidad: utilizo la ordenación por inserción para subconjuntos pequeños (<= 20 cadenas), en lugar de mergesort, como sugiere el hacker aleatorio. Aparentemente, esto funciona incluso mejor que mergesort para listas tan pequeñas (ver el primer documento que he vinculado).
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
Bueno, aquí hay una implementación simple de un tipo de raíz MSD para ADN. Está escrito en D porque ese es el idioma que uso más y, por lo tanto, es menos probable que cometa errores tontos, pero podría traducirse fácilmente a otro idioma. Está en su lugar, pero requiere 2 * seq.length pasa a través de la matriz.
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);
}
}
Obviamente, esto es algo específico para el ADN, en lugar de ser general, pero debe ser rápido.
Editar:
Me entró curiosidad si este código realmente funciona, así que lo probé / depuré mientras esperaba que se ejecutara mi propio código bioinformático. La versión anterior ahora está realmente probada y funciona. Para 10 millones de secuencias de 5 bases cada una, es aproximadamente 3 veces más rápido que un introsort optimizado.
Ciertamente puede eliminar los requisitos de memoria codificando la secuencia en bits. Usted está mirando las permutaciones entonces, para la longitud 2, con "ACGT" eso es 16 estados, o 4 pedacitos. Para la longitud 3, son 64 estados, que pueden codificarse en 6 bits. De modo que parece 2 bits para cada letra de la secuencia, o aproximadamente 32 bits para 16 caracteres como usted dijo.
Si hay una forma de reducir el número de ''palabras'' válidas, es posible que haya una mayor compresión.
Entonces, para secuencias de longitud 3, se podrían crear 64 segmentos, tal vez de tamaño uint32 o uint64. Inicialízalos a cero. Itere a través de su gran lista de 3 secuencias de caracteres y codifíquelas como se indica arriba. Use esto como un subíndice, e incremente ese cubo.
Repite esto hasta que todas tus secuencias hayan sido procesadas.
Luego, regenera tu lista.
Itere a través de los 64 intervalos en orden, para el recuento encontrado en ese segmento, genere tantas instancias de la secuencia representada por ese segmento.
cuando todos los cubos han sido iterados, usted tiene su matriz ordenada.
Una secuencia de 4, agrega 2 bits, por lo que habría 256 cubos. Una secuencia de 5, agrega 2 bits, por lo que habrá 1024 cubos.
En algún momento, la cantidad de cubos se acercará a sus límites. Si lee las secuencias de un archivo, en lugar de mantenerlas en la memoria, habrá más memoria disponible para las cubetas.
Creo que esto sería más rápido que hacer el tipo in situ, ya que es probable que los cubos quepan dentro de tu conjunto de trabajo.
Aquí hay un truco que muestra la técnica
#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;
}
En lo que respecta al rendimiento, es posible que desee ver un algoritmo de clasificación de comparación de cadenas más general.
Actualmente terminas tocando cada elemento de cada cuerda, ¡pero puedes hacerlo mejor!
En particular, un tipo de ráfaga es una muy buena opción para este caso. Como beneficio adicional, dado que burstsort se basa en tries, funciona ridículamente bien para los pequeños tamaños de alfabeto utilizados en DNA / RNA, ya que no es necesario construir ningún tipo de nodo de búsqueda ternario, hash u otro esquema de compresión de nodo trie en el trie implementación. Los intentos también pueden ser útiles para el objetivo final de su matriz de sufijos.
Una implementación decente de propósito general de burstsort está disponible en forge de origen en http://sourceforge.net/projects/burstsort/ , pero no está en su lugar.
Para propósitos de comparación, la implementación de C-burstsort cubierta en http://www.cs.mu.oz.au/~rsinha/papers/SinhaRingZobel-2006.pdf puntos de referencia 4 a 5 veces más rápido que los rangos de clasificación y radix para algunas cargas de trabajo típicas.
La ordenación de radix MSB de dsimcha se ve bien, pero Nils se acerca al corazón del problema con la observación de que la localidad de memoria caché es lo que te está matando a grandes tamaños de problema.
Sugiero un enfoque muy simple:
- Estime empíricamente el tamaño más grande
mpara el cual un tipo de radix es eficiente. - Lea bloques de
melementos a la vez, clasifíquelos y ordénelos (en un búfer de memoria si tiene suficiente memoria, pero si no en el archivo), hasta que agote su entrada. - Fusiona los bloques clasificados resultantes.
Mergesort es el algoritmo de ordenamiento más compatible con caché que conozco: "Lea el siguiente elemento de la matriz A o B, luego escriba un elemento en el búfer de salida". Se ejecuta de manera eficiente en unidades de cinta . Requiere 2n espacio para ordenar n elementos, pero mi apuesta es que la localidad de caché muy mejorada que verá lo dejará sin importancia, y si estuviera usando una ordenación de radix no en el lugar, necesitaba ese espacio adicional de todas formas.
Tenga en cuenta finalmente que mergesort se puede implementar sin recurrencia, y de hecho hacerlo de esta manera deja en claro el verdadero patrón de acceso a la memoria lineal.
Nunca he visto una ordenación de radix en el lugar, y por la naturaleza de la clase de radix dudo que sea mucho más rápido que una ordenación fuera de lugar, siempre que la matriz temporal se ajuste a la memoria.
Razón:
La clasificación hace una lectura lineal en la matriz de entrada, pero todas las escrituras serán casi al azar. Desde cierto N hacia arriba, esto se reduce a una falta de memoria caché por escritura. Esta falta de caché es lo que ralentiza su algoritmo. Si está en su lugar o no, no cambiará este efecto.
Sé que esto no responderá a su pregunta directamente, pero si la clasificación es un cuello de botella es posible que desee echar un vistazo a los algoritmos de clasificación cercanos como un paso de preprocesamiento (la página wiki en el montón dinámico puede ayudarlo a comenzar).
Eso podría dar un impulso de localidad de caché muy agradable. Un ordenamiento de radix fuera de lugar de libro de texto tendrá un mejor rendimiento. Las escrituras seguirán siendo casi aleatorias, pero al menos se agruparán en torno a los mismos fragmentos de memoria y, como tal, aumentarán la proporción de aciertos de la caché.
No tengo idea si funciona en la práctica.
Por cierto: si se trata solo de cadenas de ADN: puedes comprimir un char en dos bits y empacar tus datos bastante. Esto reducirá el requisito de memoria por factor cuatro en lugar de una representación ingenua. El direccionamiento se vuelve más complejo, pero la ALU de su CPU tiene mucho tiempo para pasar durante todas las caídas de todos modos.
Parece que ha resuelto el problema, pero, para que quede constancia, parece que una versión de una clase de base in situ viable es la "Clasificación de bandera estadounidense". Se describe aquí: Ingeniería Radix Sort . La idea general es hacer 2 pases en cada personaje; primero cuente cuántos de cada uno tiene, para poder subdividir la matriz de entrada en bandejas. Luego vuelva a pasar, intercambiando cada elemento en el contenedor correcto. Ahora recursivamente clasifique cada bin en la siguiente posición del personaje.
Primero, piensa en la codificación de tu problema. Deshágase de las cadenas, reemplácelas por una representación binaria. Use el primer byte para indicar la longitud + codificación. Alternativamente, use una representación de longitud fija en un límite de cuatro bytes. Entonces el tipo de radix se vuelve mucho más fácil. Para una ordenación de radix, lo más importante es no tener manejo de excepciones en el punto caliente del bucle interno.
OK, pensé un poco más sobre el problema 4-nary. Desea una solución como un árbol Judy para esto. La siguiente solución puede manejar cadenas de longitud variable; para una longitud fija simplemente elimine los bits de longitud, eso realmente lo hace más fácil.
Asigne bloques de 16 punteros. The least significant bit of the pointers can be reused, as your blocks will always be aligned. You might want a special storage allocator for it (breaking up large storage into smaller blocks). There are a number of different kinds of blocks:
- Encoding with 7 length bits of variable-length strings. As they fill up, you replace them by:
- Position encodes the next two characters, you have 16 pointers to the next blocks, ending with:
- Bitmap encoding of the last three characters of a string.
For each kind of block, you need to store different information in the LSBs. As you have variable length strings you need to store end-of-string too, and the last kind of block can only be used for the longest strings. The 7 length bits should be replaced by less as you get deeper into the structure.
This provides you with a reasonably fast and very memory efficient storage of sorted strings. It will behave somewhat like a trie . To get this working, make sure to build enough unit tests. You want coverage of all block transitions. You want to start with only the second kind of block.
For even more performance, you might want to add different block types and a larger size of block. If the blocks are always the same size and large enough, you can use even fewer bits for the pointers. 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.
Puede intentar usar un trie . Ordenar los datos es simplemente iterar a través del conjunto de datos e insertarlo; la estructura se clasifica de forma natural, y se puede pensar que es similar a un árbol B (excepto que en lugar de hacer comparaciones, siempre se utilizan indirecciones de puntero).
El comportamiento de almacenamiento en caché favorecerá a todos los nodos internos, por lo que probablemente no mejorará en eso; pero también puede jugar con el factor de ramificación de su trie (asegúrese de que cada nodo encaje en una única línea de caché, asigne los nodos similares a un montón, como una matriz contigua que represente un recorrido de nivel). Dado que los intentos también son estructuras digitales (O (k) insertar / encontrar / eliminar para elementos de longitud k), debe tener un rendimiento competitivo a una ordenación de radix.
Querrá echar un vistazo a Procesamiento de secuencia de genoma a gran escala por los Dres. Kasahara y Morishita.
Las cadenas compuestas por las cuatro letras de nucleótido A, C, G y T pueden codificarse especialmente en enteros para un procesamiento mucho más rápido. El género Radix es uno de los muchos algoritmos discutidos en el libro; deberías poder adaptar la respuesta aceptada a esta pregunta y ver una gran mejora en el rendimiento.
Radix-Sort no es consciente de la caché y no es el algoritmo de ordenación más rápido para grandes conjuntos. Puedes mirar:
- ti7qsort . ti7qsort es la clasificación más rápida para enteros (se puede usar para cadenas de tamaño fijo pequeño).
- QSORT en línea
- Clasificación de cuerdas
También puede usar compresión y codificar cada letra de su ADN en 2 bits antes de almacenarla en la matriz de clasificación.
Si su conjunto de datos es tan grande, entonces creo que un enfoque de buffer basado en disco sería lo mejor:
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
También experimentaría la agrupación en un mayor número de segmentos, por ejemplo, si su cadena fuera:
GATTACA
la primera llamada a MSB devolverá el depósito para GATT (256 espacios en total), de esa manera se harán menos ramas del buffer basado en disco. Esto puede o no mejorar el rendimiento, así que experimenta con él.
Voy a dar un paso en falso y sugiero que cambies a una implementación de montón / heapsort . Esta sugerencia viene con algunas suposiciones:
- Usted controla la lectura de los datos
- Puede hacer algo significativo con los datos ordenados tan pronto como ''comience'' a ordenarlos.
La belleza de la ordenación de montón / montón es que puedes construir el montón mientras lees los datos, y puedes comenzar a obtener resultados en el momento en que construyes el montón.
Retrocedamos. Si eres tan afortunado de poder leer los datos de forma asincrónica (es decir, puedes publicar algún tipo de solicitud de lectura y recibir notificaciones cuando algunos datos estén listos), y luego puedes construir una parte del montón mientras esperas el próximo trozo de datos para entrar, incluso desde el disco. A menudo, este enfoque puede enterrar la mayor parte del costo de la mitad de su clasificación detrás del tiempo dedicado a obtener los datos.
Una vez que haya leído los datos, el primer elemento ya está disponible. Dependiendo de dónde envíe los datos, esto puede ser excelente. Si lo está enviando a otro lector asíncrono, o algún modelo de "evento" paralelo, o UI, puede enviar trozos y trozos a medida que avanza.
Dicho eso, si no tiene control sobre cómo se leen los datos, y se los lee de forma sincronizada, y no tiene uso para los datos ordenados hasta que estén completamente escritos, ignore todo esto. :(
Ver los artículos de Wikipedia:
Haría una explosión de una representación de bits empaquetada. Se dice que Burstsort tiene una localidad mucho mejor que los tipos de raíz, manteniendo el uso del espacio extra con intentos de ráfaga en lugar de intentos clásicos. El papel original tiene medidas.