Clasificación de radix in situ

https://stackoverflow.com/questions/463105

19-08-2019
|

Pregunta

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 cadenas pequeñas de longitud fija que solo usan las letras & # 8220; A & # 8221 ;, & # 8220; C & # 8221 ;, & # 8220; G & # 8221; y & # 8220; T & # 8221; (sí, lo has adivinado: DNA ) que quiero ordenar.

En este momento, uso std :: sort que usa introsort en todas las implementaciones comunes del STL . Esto funciona bastante bien. Sin embargo, estoy convencido de que clasificación de radios se ajusta 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 clasificación de radix requiere copiar todos los datos (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 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 ADN como para ADN5 (que permite un carácter comodín adicional & # 8220; N & # 8221;), o incluso ADN con IUPAC códigos de ambigüedad (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 clasificación de radios 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 es casi inexistente. Hay un documento prometedor llamado Clasificación eficiente eficiente de radios in situ describe el algoritmo & # 8220; MSL & # 8221 ;. Lamentablemente, este documento también es decepcionante.

En particular, existen las siguientes cosas.

Primero, el algoritmo contiene varios errores y deja mucho sin explicar. En particular, no detalla la llamada de recursión (simplemente asumo que incrementa o reduce algún puntero para calcular los valores actuales de la máscara y el desplazamiento). Además, utiliza 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 situ 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 1

Solución

Bueno, aquí hay una implementación simple de un tipo de matriz 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 que 2 * seq.length pase 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 saber 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.

Otros consejos

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 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á a su pregunta directamente, pero si la clasificación es un cuello de botella, es posible que desee ver los algoritmos de clasificación cercana 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é.

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, entonces, para la longitud 2, con '' ACGT '' eso es 16 estados, o 4 bits. Para la longitud 3, son 64 estados, que pueden codificarse en 6 bits. Por lo tanto, se ven como 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, tal vez 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 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 sus límites. Si lee las secuencias de un archivo, en lugar de guardarlas en la memoria, habrá más memoria disponible para los cubos.

Creo que esto sería más rápido que hacer la clasificación in situ, ya que es probable que los cubos quepan en 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 apuros y sugerirle que cambie a un montón / heapsort implementación. Esta sugerencia viene con algunos supuestos:

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

La belleza del montón / heap-sort es que puedes construir el montón mientras lees los datos, y puedes comenzar a obtener resultados en el momento en que has construido el montón.

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 porción 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 de forma sincrónica, y no tienes uso para los datos ordenados hasta que estén completamente escritos, ignora todo esto. :(

Vea los artículos de Wikipedia:

Heapsort
Montón binario

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 cadena, ¡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 el estallido en ráfaga se basa en intentos, funciona ridículamente bien para los tamaños de alfabeto pequeños utilizados en 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 la forja de origen en http://sourceforge.net/projects / burstsort / , pero no está en su lugar.

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 a Procesamiento de secuencia del genoma a gran escala por los Dres. Kasahara y Morishita.

Las cadenas compuestas de 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 los muchos algoritmos discutidos en el libro; debería poder adaptar la respuesta aceptada a esta pregunta y ver una gran mejora en el rendimiento.

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

Puede 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á 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). Dado que los intentos también son estructuras digitales (O (k) insertar / buscar / eliminar elementos de longitud k), debe tener un rendimiento competitivo para una clasificación de radix.

Me gustaría burstsort una representación empaquetada 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 ordenación más rápido para conjuntos grandes. Puedes mirar:

ti7qsort . ti7qsort es el tipo más rápido 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.

La clasificación de radios 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 problemas de gran tamaño.

Sugiero un enfoque muy simple:

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

Mergesort es el algoritmo de clasificación más amigable con el caché que conozco: "Leer el siguiente elemento de la matriz A o B, luego escribir un elemento en el búfer de salida". Se ejecuta eficientemente en unidades de cinta . Requiere espacio 2n para ordenar los elementos n , pero mi apuesta es que la localidad de caché mejorada que verá hará que eso no sea importante, y si estuviera usando una clasificación de radix no in situ, de todos modos necesitabas ese espacio extra.

Tenga en cuenta finalmente que mergesort puede implementarse sin recurrencia, 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 americana". Se describe aquí: Engineering 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 bins. Luego vuelva a pasar, intercambiando cada elemento en el contenedor correcto. Ahora ordena recursivamente 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.

Bien, pensé un poco más sobre el problema 4-nary. Desea una solución como un Judy tree 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 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 un rendimiento aún mayor, 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 de 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.

Licenciado bajo: CC-BY-SA con atribución

No afiliado a StackOverflow