Este es un texto largo. Por favor, tenga paciencia conmigo. En resumen, la pregunta es: ¿existe un algoritmo de clasificación de radix in situ viable ?

Preliminar

Tengo una gran cantidad de pequeñas cadenas de longitud fija que solo usan las letras "A", "C", "G" y "T" (sí, lo has adivinado: ADN ) que quiero clasificar.

Por el momento, uso el std::sortque 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 adapta perfectamente a mi conjunto de problemas 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 era 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: la ordenación de radix requiere copiar todos los datos (más de una vez en mi ingenua implementación, 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 del problema en realidad se vuelven aún mayores.

Casos de uso

Idealmente, este algoritmo debería funcionar con cualquier longitud de cadena entre 2 y 100, tanto para el ADN como para el ADN5 (que permite un carácter comodín adicional "N"), o incluso ADN con códigos de ambigüedad IUPAC (que dan como resultado 16 valores distintos). Sin embargo, me doy cuenta de que todos estos casos no pueden ser cubiertos, así que estoy contento con cualquier mejora de velocidad que obtengo. El código puede decidir dinámicamente a qué algoritmo enviar.

Investigación

Desafortunadamente, el artículo de Wikipedia sobre la clasificación de radix es inútil. La sección sobre una variante en el lugar es basura completa. La sección NIST-DADS en clasificación de radios está al lado de inexistente. Hay un documento prometedor llamado Efficient Adaptive In-Place Radix Sorting que describe el algoritmo "MSL". Desafortunadamente, este documento también es decepcionante.

En particular, hay las siguientes cosas.

Primero, el algoritmo contiene varios errores y deja mucho sin explicar. En particular, no detalla la llamada de recursión (simplemente supongo que incrementa o reduce algún puntero para calcular los valores actuales de desplazamiento y máscara). Además, utiliza las funciones dest_groupy dest_addresssin dar definiciones. No veo cómo implementarlos de manera eficiente (es decir, en O (1); al menos dest_addressno es trivial).

Por último, pero no menos importante, el algoritmo logra el lugar intercambiando í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 simplemente teclear con fuerza y ​​seguir adelante asumiendo que la memoria tolerará que almacene un índice donde no pertenece. Pero esto solo funciona mientras pueda exprimir 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 artículo de uno de los autores no ofrece una descripción precisa, pero da el tiempo de ejecución de MSL como sub-lineal, lo cual es completamente incorrecto.

Para recapitular : ¿Hay alguna esperanza de encontrar una implementación de referencia que funcione o al menos un buen pseudocódigo / descripción de un tipo de radix en el lugar que funcione que funcione en cadenas de ADN?

Bueno, aquí hay una implementación simple de un tipo de radix MSD para ADN. Está escrito en D porque ese es el idioma que más uso 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.lengthpases 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 específico del ADN, en lugar de ser general, pero debería ser rápido.

Editar:

Me dio 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á 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.

Nunca he visto una clasificación de radix in situ, y por la naturaleza de la clasificación de radix dudo que sea mucho más rápido que una clasificación fuera de lugar, siempre y cuando la matriz temporal se ajuste a la memoria.

Razón:

La clasificación realiza una lectura lineal en la matriz de entrada, pero todas las escrituras serán casi aleatorias. Desde un cierto N hacia arriba, esto se reduce a una pérdida de caché por escritura. Este error de caché es lo que ralentiza su algoritmo. Si está en su lugar o no, no cambiará este efecto.

Sé que esto no responderá 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. Una clasificación de radix fuera de lugar de un libro de texto funcionará mejor. Las escrituras seguirán siendo casi aleatorias, pero al menos se agruparán alrededor de los mismos fragmentos de memoria y, como tal, aumentarán la proporción de aciertos de caché.

Sin embargo, no tengo idea si funciona en la práctica.

Por cierto: si solo se trata de cadenas de ADN: puede comprimir un carácter en dos bits y empacar sus datos bastante. Esto reducirá el requisito de memoria por el factor cuatro sobre una representación ingenua. El direccionamiento se vuelve más complejo, pero la ALU de su CPU tiene mucho tiempo para pasar durante todos los fallos de caché de todos modos.

Ciertamente, puede eliminar los requisitos de memoria codificando la secuencia en bits. Estás viendo permutaciones, por lo que, para la longitud 2, con "ACGT" son 16 estados, o 4 bits. Para la longitud 3, son 64 estados, que pueden codificarse en 6 bits. Por lo tanto, parece 2 bits para cada letra en la secuencia, o aproximadamente 32 bits para 16 caracteres como usted dijo.

Si hay una manera de reducir el número de 'palabras' válidas, puede ser posible una mayor compresión.

Entonces, para secuencias de longitud 3, uno podría crear 64 cubos, quizás de tamaño uint32 o uint64. Inicialícelos a cero. Itere a través de su muy extensa lista de 3 secuencias de caracteres y codifíquelas como se indica arriba. Use esto como un subíndice e incremente ese depósito.
Repita esto hasta que todas sus secuencias hayan sido procesadas.

Luego, regenera tu lista.

Itere a través de los 64 cubos en orden, para el recuento encontrado en ese cubo, genere tantas instancias de la secuencia representada por ese cubo.
cuando se han iterado todos los cubos, 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ía 1024 cubos.

En algún momento, el número de cubos se acercará a tus límites. Si lee las secuencias de un archivo, en lugar de guardarlas en la memoria, habrá más memoria disponible para los depósitos.

Creo que esto sería más rápido que hacer la clasificación in situ, ya que es probable que los cubos se ajusten a su 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;
}

Si su conjunto de datos es tan grande, entonces pensaría que un enfoque de búfer basado en disco sería 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 cubos, por ejemplo, si su cadena era:

GATTACA

la primera llamada de MSB devolvería el depósito para GATT (256 depósitos totales), de esa forma se crean menos ramas del almacenamiento intermedio basado en disco. Esto puede o no mejorar el rendimiento, así que experimente con él.

Voy a salir en una extremidad y sugieren que cambie a un montón / heapsort aplicación. Esta sugerencia viene con algunos supuestos:

1. Usted controla la lectura de los datos.
2. Puede hacer algo significativo con los datos ordenados tan pronto como 'comience' a ordenarlos.

La belleza del montón / clasificación de montón es que puede construir el montón mientras lee los datos, y puede comenzar a obtener resultados en el momento en que lo haya hecho.

Retrocedamos Si tiene la suerte de poder leer los datos de forma asincrónica (es decir, puede publicar algún tipo de solicitud de lectura y recibir una notificación cuando haya algunos datos listos), y luego puede construir una parte del montón mientras espera el siguiente fragmento de datos que ingresará, 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 paralelo de 'evento', o UI, puede enviar fragmentos y fragmentos a medida que avanza.

Dicho esto, si no tienes control sobre cómo se leen los datos, y si se leen sincrónicamente, y no tienes uso para los datos ordenados hasta que se escriban por completo, ignora todo esto. :(

Ver los artículos de Wikipedia:

• Heapsort
• Montón binario

" Clasificación de radix sin espacio adicional " es un documento que aborda su problema.

En cuanto al rendimiento, es posible que desee ver algoritmos de clasificación de comparación de cadenas más generales.

Actualmente terminas tocando cada elemento de cada cuerda, ¡pero puedes hacerlo mejor!

En particular, un estallido es una muy buena opción para este caso. Como beneficio adicional, dado que la ráfaga de ráfagas se basa en intentos, funciona ridículamente bien para los pequeños tamaños de alfabeto utilizados en el ADN / ARN, 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 su objetivo final tipo matriz de sufijo.

Una implementación decente de propósito general de burstsort está disponible en source forge en http://sourceforge.net/projects/burstsort/ , pero no está implementada.

Para fines de comparación, la implementación de C-burstsort se cubre en http://www.cs.mu.oz.au/~rsinha/papers/SinhaRingZobel-2006.pdf puntos de referencia 4-5 veces más rápidos que los de clasificación rápida y tipo de raíz para algunas cargas de trabajo típicas.

Querrá echar un vistazo al procesamiento de secuencia del genoma a gran escala por los Dres. Kasahara y Morishita.

Las cadenas compuestas por las cuatro letras de nucleótidos A, C, G y T pueden codificarse especialmente en números enteros para un procesamiento mucho más rápido. La clasificación por radix se encuentra entre muchos algoritmos discutidos en el libro; debe poder adaptar la respuesta aceptada a esta pregunta y ver una gran mejora en el rendimiento.

Puedes intentar usar un trie . Ordenar los datos es simplemente iterar a través del conjunto de datos e insertarlos; la estructura está naturalmente ordenada, y puede pensar que es similar a un B-Tree (excepto que en lugar de hacer comparaciones, siempre usa 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 sola línea de caché, asigne nodos de trie similares a un montón, como una matriz contigua que representa un recorrido de orden de nivel). Como los intentos también son estructuras digitales (O (k) insert / find / delete para elementos de longitud k), debe tener un rendimiento competitivo para una clasificación de radix.

Yo clasificaría una representación en bits de las cadenas. Se dice que Burstsort tiene una localidad mucho mejor que los tipos de radix, manteniendo el uso de espacio adicional bajo con intentos de ráfaga en lugar de intentos clásicos. El papel original tiene medidas.

Radix-Sort no es consciente de la memoria caché y no es el algoritmo de clasificación más rápido para conjuntos grandes. Puedes mirar:

• ti7qsort . ti7qsort es la clasificación más rápida para enteros (se puede usar para cadenas de tamaño pequeño y fijo).
• QSORT en línea
• Clasificación de cadenas

También puede usar la compresión y codificar cada letra de su ADN en 2 bits antes de almacenarla en la matriz de clasificación.

El tipo de raíz MSB de dsimcha se ve bien, pero Nils se acerca al corazón del problema con la observación de que la localidad de caché es lo que te está matando en grandes problemas.

Sugiero un enfoque muy simple:

1. Estima empíricamente el tamaño más grande mpara el cual una clasificación de radix es eficiente.
2. Lea bloques de melementos a la vez, clasifíquelos en radix y escríbalos (en un búfer de memoria si tiene suficiente memoria, pero de lo contrario para archivar), hasta agotar su entrada.
3. Combina los bloques ordenados resultantes.

Mergesort es el algoritmo de clasificación más amigable con el 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 eficientemente en unidades de cinta . Requiere 2nespacio para ordenar nelementos, pero mi apuesta es que la localidad de caché mejorada que verá hará que eso no sea importante, y si estaba usando una clasificación de radix no in situ, de todos modos necesitaba ese espacio adicional.

Finalmente, tenga en cuenta que mergesort puede implementarse sin recursividad, y de hecho al hacerlo de esta manera deja en claro el verdadero patrón de acceso lineal a la memoria.

Parece que ha resuelto el problema, pero para el registro, parece que una versión de una clasificación de radix in situ viable es la "Clasificación de la 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 que pueda subdividir la matriz de entrada en contenedores. Luego vuelva a pasar, intercambiando cada elemento en el contenedor correcto. Ahora recursivamente ordena cada bin en la siguiente posición de personaje.

Primero, piense en la codificación de su problema. Deshágase de las cadenas, reemplácelas por una representación binaria. Use el primer byte para indicar longitud + codificación. Alternativamente, use una representación de longitud fija en un límite de cuatro bytes. Entonces la clasificación de radix se vuelve mucho más fácil. Para una clasificación de radix, lo más importante es no tener un manejo de excepciones en el punto caliente del bucle interno.

OK, pensé un poco más sobre el problema 4-nary. Quieres una solución como un árbol Judy para esto. La siguiente solución puede manejar cadenas de longitud variable; para una longitud fija, simplemente quite los bits de longitud, eso en realidad lo hace más fácil

Asignar bloques de 16 punteros. El bit menos significativo de los punteros se puede reutilizar, ya que sus bloques siempre estarán alineados. Es posible que desee un asignador de almacenamiento especial para él (dividiendo el almacenamiento grande en bloques más pequeños). Hay varios tipos diferentes de bloques:

• Codificación con 7 bits de longitud de cadenas de longitud variable. A medida que se llenan, los reemplaza por:
• La posición codifica los siguientes dos caracteres, tiene 16 punteros a los siguientes bloques, que terminan con:
• Codificación de mapa de bits de los últimos tres caracteres de una cadena.

Para cada tipo de bloque, debe almacenar información diferente en los LSB. Como tiene cadenas de longitud variable, también necesita almacenar el final de la cadena, y el último tipo de bloque solo se puede usar para las cadenas más largas. Los 7 bits de longitud deben reemplazarse por menos a medida que profundiza en la estructura.

Esto le proporciona un almacenamiento razonablemente rápido y muy eficiente en la memoria de cadenas ordenadas. Se comportará como un trie . Para que esto funcione, asegúrese de construir suficientes pruebas unitarias. Desea cobertura de todas las transiciones de bloque. Desea comenzar solo con el segundo tipo de bloque.

Para obtener aún más rendimiento, es posible que desee agregar diferentes tipos de bloque y un tamaño de bloque más grande. Si los bloques son siempre del mismo tamaño y lo suficientemente grandes, puede usar incluso menos bits para los punteros. Con un tamaño de bloque de 16 punteros, ya tiene un byte libre en un espacio de direcciones de 32 bits. Eche un vistazo a la documentación del árbol Judy para ver tipos de bloques interesantes. Básicamente, agrega código y tiempo de ingeniería para un intercambio de espacio (y tiempo de ejecución)

Probablemente quiera comenzar con una raíz directa de 256 ancho para los primeros cuatro caracteres. Eso proporciona una compensación decente espacio / tiempo. En esta implementación, obtienes mucha menos sobrecarga de memoria que con un simple trie; Es aproximadamente tres veces más pequeño (no lo he medido). O (n) no es un problema si la constante es lo suficientemente baja, como se notó al comparar con la clasificación rápida O (n log n).

¿Estás interesado en manejar dobles? Con secuencias cortas, habrá. Adaptar los bloques para manejar conteos es complicado, pero puede ser muy eficiente en cuanto al espacio.