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

Razón:

La ordenación hace una lectura lineal en la matriz de entrada, pero todas las escrituras serán casi aleatorias. Desde una cierta N hacia arriba esto se reduce a una falta de caché por escritura. Esta falta de caché es lo que ralentiza tu algoritmo. Si está en su lugar o no no cambiará este efecto.

Sé que esto no responderá a su pregunta directamente, pero si la ordenación es un cuello de botella, es posible que desee echar un vistazo a cerca de los algoritmos de ordenación como un paso de preprocesamiento (la página wiki en el montón de software puede ayudarlo a comenzar).

Eso podría dar un muy buen impulso de localización de caché. Un libro de texto fuera de lugar radix sort funcionará mejor. Las escrituras seguirán siendo casi aleatorias, pero al menos se agruparán alrededor de los mismos trozos de memoria y, como tal aumentar la proporción de aciertos de caché.

No tengo idea si funciona en la práctica, sin embargo.

Por cierto: Si estás tratando solo con cadenas de ADN: Puedes comprimir un char en dos bits y empacar tus datos bastante. Esto reducirá el requisito de memoria por factor cuatro sobre una representación naiive. El direccionamiento se vuelve más complejo, pero la ALU de su CPU tiene mucho tiempo para pasar durante todas las pérdidas de caché de todos modos.

Basado en el código de dsimcha, he implementado una versión más genérica que encaja bien en el framework que usamos (SeqAn). En realidad, portar el código fue completamente sencillo. Solo después descubrí que hay en realidad publicaciones relacionadas con este mismo tema. Lo bueno es que básicamente dicen lo mismo que ustedes. Un artículo de Andersson y Nilsson sobre Implementando Radixsort definitivamente vale la pena leerlo. Si por casualidad sabes alemán, asegúrate de leer también David Weese tesis de diploma donde implementa un índice de subcadenas 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 en realidad se han implementado en SeqAn, solo que no en aquellas 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). Observe que todavía no considera radios mayores 4. Espero que esto tenga un gran impacto en el rendimiento, pero por desgracia, simplemente no tengo el tiempo en este momento para implementar esto.

El código funciona más del doble de rápido que Introsort para cadenas cortas. El punto de equilibrio tiene 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 de lecturas de ADN del cromosoma 2 del genoma humano toma poco más de 2 segundos en mi PC. Solo para que conste, eso es rápido ! Especialmente teniendo en cuenta que no uso SIMD o cualquier otra aceleración de hardware. Además, valgrind me muestra que el principal cuello de botella es operator new en las asignaciones de cadenas. Se llama alrededor de 65.000.000 veces-diez veces por cada cadena! Esto es un claro indicio de que swap podría optimizarse para estas cadenas: en lugar de hacer copias, podría simplemente intercambiar todos los caracteres. No lo intenté, pero estoy convencido de que haría una gran diferencia. Y, por decirlo de nuevo, en caso de que alguien no estuviera escuchando: el tamaño radix casi no tiene influencia en el tiempo de ejecución, lo que significa que debería definitivamente 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 arreglarse fácilmente: uso la ordenación por inserción para subconjuntos pequeños (

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

Ciertamente puede eliminar los requisitos de memoria codificando la secuencia en bits. Usted está mirando permutaciones así, para la longitud 2, con "ACGT" que es 16 estados, o 4 bits. Para la longitud 3, eso es 64 estados, que se pueden codificar en 6 bits. Así que se ve como 2 bits para cada letra en la secuencia, o alrededor de 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.

Así que para secuencias de longitud 3, uno podría crear 64 cubos, tal vez del tamaño uint32 o uint64. Inicializarlos a cero. Itere a través de su lista muy muy grande de 3 secuencias de caracteres, y codifíquelas como arriba. Use esto como un subíndice, e incremente ese cubo.
Repita esto hasta que todas sus secuencias hayan sido procesadas.

A continuación, regenera tu lista.

Iterar a través de los 64 cubos en orden, para el recuento encontrado en ese cubo, generar que muchas instancias de la secuencia representada por que cubo.
cuando todos los cubos se han iterado, 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 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 los cubos son es probable que quepa dentro de su equipo 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 creo 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 agrupando en un mayor número de cubos, por ejemplo, si su cadena fuera:

GATTACA

La primera llamada a MSB devolvería el bucket para GATT (256 buckets totales), de esa manera se hacen menos ramas del buffer basado en disco. Esto puede o no mejorar el rendimiento, así que experimenta con él.

Voy a arriesgarme y sugerirle que cambie a una implementación heap/ heapsort. Esta sugerencia viene con algunas suposiciones:

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 de heap / heap-sort es que puede construir el montón mientras lee los datos, y puede comenzar a obtener resultados en el momento en que ha construido el pila.

Una vez que haya leído los datos, el primer elemento ya está disponible. Dependiendo de dónde esté enviando los datos, esto puede ser genial. Si lo está enviando a otro lector asíncrono, o a algún modelo de 'evento' paralelo, o interfaz de usuario, puede enviar trozos y trozos a medida que avanza.

Dicho esto, si no tiene control sobre cómo se leen los datos, y se leen sincrónicamente, y no tiene uso para los datos ordenados hasta que se escriben completamente, ignore todo esto. :(

Ver la Wikipedia artículos:

• Heapsort
• Montón binario

En cuanto al rendimiento, es posible que desee mirar un algoritmo de clasificación de comparación de cadenas más general.

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

En particular, un burst sort es una muy buena opción para este caso. Como ventaja adicional, ya que burstsort se basa en tries, funciona ridículamente bien para los pequeños tamaños de alfabeto utilizados en ADN / ARN, ya que no necesita construir ningún tipo de nodo de búsqueda ternaria, hash u otro nodo trie esquema de compresión en la implementación de trie. Los intentos también pueden ser útiles para su objetivo final de matriz de sufijos.

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

Para fines de comparación, la implementación de C-burstsort abarcada en http://www.cs.mu.oz.au/~rsinha/papers / SinhaRingZobel-2006.pdf puntos de referencia 4-5 veces más rápidos que quicksort y radix ordena algunas cargas de trabajo típicas.

Puede intentar usar un trie. Ordenar los datos es simplemente iterar a través del conjunto de datos e insertarlo; la estructura se ordena naturalmente, y se puede pensar en ella como similar a un árbol B (excepto que en lugar de hacer comparaciones, siempre utiliza indirectas 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 se ajuste a un una sola línea de caché, asignar nodos 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) insert/find/delete para elementos de longitud k), debería tener un rendimiento competitivo para una clasificación radix.

Me gustaría burstsort una representación de bits empaquetados de las cadenas. Se afirma que Burstsort tiene una ubicación mucho mejor que radix sorts, manteniendo el uso de espacio extra bajo con intentos de ráfaga en lugar de intentos clásicos. El papel original tiene medidas.

Usted querrá echar un vistazo a Procesamiento de Secuencias Genómicas a Gran Escala por los doctores Kasahara y Morishita.

Las cadenas compuestas por las cuatro letras de nucleótidos A, C, G y T pueden codificarse especialmente en enteros para un procesamiento mucho más rápido. Radix sort es uno de los muchos algoritmos discutidos en el libro; debería ser capaz de adaptar la respuesta aceptada a esta pregunta y ver una gran mejora en el rendimiento.

"Radix sorting with no extra space " es un documento que aborda su problema.

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

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

También puede utilizar la compresión y codificar cada letra de su ADN en 2 bits antes de almacenar en la matriz de ordenación.

La ordenación de MSB radix 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 tamaños de problemas.

Sugiero un enfoque muy simple:

1. Calcule empíricamente el tamaño más grande m para el cual una clasificación radix es eficiente.
2. Leer bloques de elementos m a la vez, radix ordenarlos, y escribirlos (a un búfer de memoria si tiene suficiente memoria, pero de lo contrario a archivo), hasta que agota tu entrada.
3. Mergesort los bloques ordenados resultantes.

Mergesort es el algoritmo de clasificación más amigable con la 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 accionadors de cinta. Requiere 2n espacio para ordenar los elementos n, pero mi apuesta es que la localidad de caché mucho mejor que verá hará que eso carezca de importancia -- y si estuviera utilizando una ordenación radix no en el lugar, necesitabas ese espacio extra de todos modos.

Tenga en cuenta finalmente que mergesort se puede implementar sin recursión, y de hecho hacerlo de esta manera deja claro el verdadero patrón de acceso a la memoria lineal.

Parece que ha resuelto el problema, pero para el registro, parece que una versión de una ordenación radix en el lugar viable es la "Ordenación de bandera americana". Se describe aquí: Engineering Radix Sort. La idea general es hacer 2 pases en cada carácter-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 la bandeja correcta. Ahora ordena recursivamente cada bin en la siguiente posición de carácter.

Primero, piensa en la codificación de tu problema. Deshazte de las cadenas, reemplázalas por una representación binaria. Utilice el primer byte para indicar longitud + codificación. Alternativamente, use una representación de longitud fija en un límite de cuatro bytes. Entonces el orden radix se vuelve mucho más fácil. Para una clasificación 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. Usted quiere una solución como un Judy árbol para esto. La siguiente solución puede manejar cadenas de longitud variable; para una longitud fija, simplemente elimine los bits de longitud, lo que en realidad lo hace más fácil.

Asigna bloques de 16 punteros. La parte menos significativa 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 (dividir el almacenamiento grande en bloques más pequeños). Hay un número de diferentes tipos de bloques:

• Codificación con 7 bits de longitud de cuerdas de longitud variable. A medida que se llenan, los reemplazas por:
• Position codifica los siguientes dos caracteres, tienes 16 punteros a los siguientes bloques, terminando con:
• Codificación de mapa de bits de los tres últimos caracteres de una cadena.

Para cada tipo de bloque, necesita almacenar información diferente en los LSBs. 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 debe ser reemplazado por menos a medida que se profundiza en la estructura.

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

Para obtener aún más rendimiento, es posible que desee agregar diferentes tipos de bloques y un tamaño de bloque más grande. Si los bloques siempre son 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. Echa un vistazo a la documentación de Judy Tree para encontrar tipos de bloques interesantes. Básicamente, se agrega código y tiempo de ingeniería para un espacio (y tiempo de ejecución) trade-off

Probablemente quieras empezar con una raíz directa de 256 de 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ña (no he medido). O(n) no es un problema si la constante es lo suficientemente baja, como notó al comparar con el O (n log n) quicksort.

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