Ordenación por base in situ

Este es un texto largo. Por favor, tenga paciencia conmigo. En resumen, la pregunta es: ¿ Existe un algoritmo de clasificación de bases 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 ordenar.

Por el momento, utilizo std::sortel que utiliza introsort en todas las implementaciones comunes de STL . Esto funciona bastante bien. Sin embargo, estoy convencido de que la clasificación por base se adapta perfectamente a mi conjunto de problemas y debería funcionar mucho mejor en la práctica.

Detalles

Probé 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: la clasificación por base requiere copiar todos los datos (en realidad, más de una vez en mi implementación ingenua). Esto significa que he puesto ~ 4 GiB en mi memoria principal, lo que obviamente acaba con el rendimiento. Incluso si no fuera así, no puedo permitirme el lujo de utilizar tanta memoria ya que los tamaños de los problemas 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 “N”), o incluso ADN con códigos de ambigüedad IUPAC (lo que da como resultado 16 valores distintos). Sin embargo, me doy cuenta de que no se pueden cubrir todos estos casos, por lo 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 por bases es inútil. La sección sobre una variante in situ es una completa tontería. La sección NIST-DADS sobre clasificación de bases es casi inexistente. Hay un artículo que suena prometedor llamado Efficient Adaptive In-Place Radix Sorting que describe el algoritmo "MSL". Desafortunadamente, este artículo también es decepcionante.

En particular, existen las siguientes cosas.

En primer lugar, el algoritmo contiene varios errores y deja muchas cosas sin explicar. En particular, no detalla la llamada de recursividad (simplemente asumo 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 estar en el lugar intercambiando índices de matriz con elementos dentro de la matriz de entrada. Obviamente, esto sólo funciona en matrices numéricas. Necesito usarlo en cuerdas. Por supuesto, podría simplemente fastidiar la escritura estricta y seguir adelante suponiendo que la memoria tolerará que almacene un índice donde no pertenece. Pero esto sólo funciona siempre que pueda comprimir mis cadenas en 32 bits de memoria (asumiendo números enteros de 32 bits). Son sólo 16 caracteres (ignoremos por el momento que 16 > log(5,000,000)).

Otro artículo de uno de los autores no ofrece ninguna descripción precisa, pero indica que el tiempo de ejecución de MSL es sublineal, lo cual es completamente incorrecto.

En resumen : ¿Existe alguna esperanza de encontrar una implementación de referencia que funcione o al menos un buen pseudocódigo/descripción de un tipo de base que funcione en el lugar y que funcione en cadenas de ADN?

Jan 21 '09 04:01 Konrad Rudolph

Nunca he visto una clasificación de base in situ y, por la naturaleza de la clasificación de base, dudo que sea mucho más rápida que una clasificación fuera de lugar siempre que la matriz temporal quepa en la memoria.

Razón:

La clasificación realiza una lectura lineal en la matriz de entrada, pero todas las escrituras serán casi aleatorias. A partir de un cierto N, esto se reduce a una pérdida de caché por escritura. Esta pérdida 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 cercana como paso previo al procesamiento (la página wiki en el montón suave puede ayudarlo a comenzar).

Eso podría dar un muy buen impulso a la localidad de caché. Una clasificación de base 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 tasa de aciertos de la caché.

Aunque 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 empaquetar una gran cantidad de datos. Esto reducirá el requisito de memoria en un factor cuatro en comparación con una representación ingenua. El direccionamiento se vuelve más complejo, pero la ALU de su CPU tiene mucho tiempo para gastar durante todos los errores de caché de todos modos.

Jan 20 '2009 21:01 Nils Pipenbrinck

Ciertamente puedes eliminar los requisitos de memoria codificando la secuencia en bits. Estás viendo permutaciones, por lo que, para longitud 2, con "ACGT" son 16 estados o 4 bits. Para una longitud 3, son 64 estados, que se pueden codificar en 6 bits. Entonces parece que hay 2 bits para cada letra de la secuencia, o alrededor de 32 bits para 16 caracteres como dijiste.

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

Entonces, para secuencias de longitud 3, se podrían crear 64 depósitos, tal vez de tamaño uint32 o uint64. Inicialícelos a cero. Repita su lista muy grande de secuencias de 3 caracteres y codifíquelas como se indica arriba. Utilice esto como subíndice e incremente ese depósito.
Repita esto hasta que todas sus secuencias hayan sido procesadas.

A continuación, regenera tu lista.

Itere a través de los 64 depósitos en orden, para el recuento encontrado en ese depósito, genere esa cantidad de instancias de la secuencia representada por ese depósito.
cuando se hayan iterado todos los depósitos, tendrá su matriz ordenada.

Una secuencia de 4 agrega 2 bits, por lo que habría 256 depósitos. Una secuencia de 5 agrega 2 bits, por lo que habría 1024 depósitos.

En algún momento, la cantidad de depósitos 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 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 quepan dentro de 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;
}

Jan 20 '2009 22:01 EvilTeach

Si su conjunto de datos es tan grande, entonces creo que lo mejor sería un enfoque de búfer basado en disco:

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 agrupando en una mayor cantidad de depósitos, por ejemplo, si su cadena fuera:

GATTACA

la primera llamada a MSB devolvería el depósito para GATT (256 depósitos en total), de esa manera se crean menos ramas del búfer basado en disco. Esto puede mejorar o no el rendimiento, así que experimente con ello.

Jan 20 '2009 21:01 FryGuy