Este é um texto longo. Por favor, tenha paciência comigo. Resumindo, a questão é: existe um algoritmo de classificação de raiz no local viável ?

Preliminares

Eu tenho um grande número de pequenas cadeias de comprimento fixo que usam apenas as letras "A", "C", "G" e "T" (sim, você adivinhou: DNA ) que quero classificar.

No momento, eu uso o std::sortque usa introsort em todas as implementações comuns do STL . Isso funciona muito bem. No entanto, estou convencido de que a classificação por radix se encaixa perfeitamente no meu conjunto de problemas e deve funcionar muito melhor na prática.

Detalhes

Testei essa suposição com uma implementação muito ingênua e, para entradas relativamente pequenas (da ordem de 10.000), isso era verdade (bem, pelo menos mais do que o dobro da velocidade). No entanto, o tempo de execução diminui enormemente quando o tamanho do problema se torna maior ( N > 5.000.000).

O motivo é óbvio: a classificação radix requer a cópia de todos os dados (mais de uma vez na minha implementação ingênua, na verdade). Isso significa que coloquei ~ 4 GiB na memória principal, o que obviamente mata o desempenho. Mesmo se não, não posso me dar ao luxo de usar tanta memória, pois os tamanhos dos problemas realmente se tornam ainda maiores.

Casos de Uso

Idealmente, esse algoritmo deve funcionar com qualquer comprimento de cadeia entre 2 e 100, tanto para o DNA quanto para o DNA5 (que permite um caractere curinga adicional "N"), ou mesmo DNA com códigos de ambiguidade IUPAC (resultando em 16 valores distintos). No entanto, percebo que todos esses casos não podem ser cobertos, por isso estou feliz com qualquer melhoria de velocidade que recebo. O código pode decidir dinamicamente para qual algoritmo enviar.

Pesquisa

Infelizmente, o artigo da Wikipedia sobre classificação radix é inútil. A seção sobre uma variante no local é um lixo completo. A seção NIST-DADS na classificação de raiz é quase inexistente. Existe um artigo que parece promissor, chamado Efficient Adaptive Radix Sorting, que descreve o algoritmo "MSL". Infelizmente, este artigo também é decepcionante.

Em particular, existem as seguintes coisas.

Primeiro, o algoritmo contém vários erros e deixa muita coisa inexplicável. Em particular, ele não detalha a chamada de recursão (simplesmente presumo que ela aumente ou reduz algum ponteiro para calcular os valores atuais de mudança e máscara). Além disso, ele usa as funções dest_groupe dest_addresssem fornecer definições. Não consigo entender como implementá-las com eficiência (ou seja, em O (1); pelo menos dest_addressnão é trivial).

Por último, mas não menos importante, o algoritmo alcança o local, trocando os índices da matriz por elementos dentro da matriz de entrada. Obviamente, isso funciona apenas em matrizes numéricas. Eu preciso usá-lo em cordas. Claro, eu poderia simplesmente estragar a digitação forte e seguir em frente, assumindo que a memória toleraria que eu armazenasse um índice onde não pertence. Mas isso só funciona desde que eu possa espremer minhas strings em 32 bits de memória (assumindo números inteiros de 32 bits). São apenas 16 caracteres (vamos ignorar no momento que 16> log (5.000.000)).

Outro artigo de um dos autores não fornece uma descrição precisa, mas fornece o tempo de execução do MSL como sub-linear, o que é totalmente errado.

Recapitulando : Existe alguma esperança de encontrar uma implementação de referência de trabalho ou pelo menos um bom pseudocódigo / descrição de uma classificação de raiz no local que funcione que funcione em seqüências de DNA?

Bem, aqui está uma implementação simples de uma classificação de radix MSD para DNA. Está escrito em D, porque esse é o idioma que eu mais uso e, portanto, é menos provável que eu cometa erros bobos, mas poderia ser facilmente traduzido para outro idioma. Está no lugar, mas requer 2 * seq.lengthpassagens pela 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, isso é específico do DNA, em vez de geral, mas deve ser rápido.

Editar:

Fiquei curioso para saber se esse código realmente funciona, então testei / depurei enquanto esperava meu código de bioinformática ser executado. A versão acima agora é realmente testada e funciona. Para 10 milhões de seqüências de 5 bases cada, é cerca de 3x mais rápido que um introsort otimizado.

Eu nunca vi uma classificação radix no local e, pela natureza da classificação radix, duvido que seja muito mais rápido que uma classificação fora do lugar, desde que a matriz temporária caiba na memória.

Razão:

A classificação faz uma leitura linear na matriz de entrada, mas todas as gravações serão quase aleatórias. De um certo N em diante, isso se resume a uma falta de cache por gravação. Essa falta de cache é o que torna seu algoritmo mais lento. Se estiver no lugar ou não, isso não mudará esse efeito.

Sei que isso não responderá diretamente à sua pergunta, mas se a classificação for um gargalo, convém dar uma olhada nos algoritmos de classificação próximos como uma etapa de pré - processamento (a página da wiki na pilha flexível pode ajudar você a começar).

Isso poderia dar um impulso muito bom na localidade do cache. Uma classificação de raiz fora do local do livro-texto terá um desempenho melhor. As gravações ainda serão quase aleatórias, mas pelo menos elas se agruparão em torno dos mesmos pedaços de memória e, como tal, aumentam a taxa de acertos do cache.

Não tenho idéia se funciona na prática.

Btw: Se você estiver lidando apenas com seqüências de DNA: você pode compactar um caractere em dois bits e compactar bastante seus dados. Isso reduzirá o requisito de memória pelo fator quatro, sobre uma representação ingênua. O endereçamento se torna mais complexo, mas a ALU da sua CPU tem muito tempo para gastar durante todas as falhas de cache de qualquer maneira.

Certamente você pode eliminar os requisitos de memória codificando a sequência em bits. Você está procurando permutações, portanto, para o comprimento 2, com "ACGT", que é 16 estados ou 4 bits. Para o comprimento 3, são 64 estados, que podem ser codificados em 6 bits. Portanto, parece 2 bits para cada letra na sequência ou cerca de 32 bits para 16 caracteres, como você disse.

Se houver uma maneira de reduzir o número de 'palavras' válidas, uma compressão adicional pode ser possível.

Portanto, para seqüências de comprimento 3, é possível criar 64 buckets, talvez do tamanho uint32 ou uint64. Inicialize-os para zero. Faça uma iteração na sua lista muito grande de três seqüências de caracteres e codifique-as como acima. Use isso como um subscrito e incremente esse intervalo.
Repita isso até que todas as suas seqüências tenham sido processadas.

Em seguida, gere novamente sua lista.

Faça uma iteração nos 64 depósitos na ordem, para a contagem encontrada nesse depósito, gere muitas instâncias da sequência representada por esse depósito.
quando todos os buckets tiverem sido iterados, você terá sua matriz classificada.

Uma sequência de 4 adiciona 2 bits, para que houvesse 256 buckets. Uma sequência de 5 adiciona 2 bits, para que haja 1024 buckets.

Em algum momento, o número de buckets chegará aos seus limites. Se você ler as seqüências de um arquivo, em vez de mantê-las na memória, mais memória estará disponível para os buckets.

Eu acho que isso seria mais rápido do que fazer o tipo in situ, pois os baldes provavelmente caberão no seu conjunto de trabalho.

Aqui está um hack que mostra a 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;
}

Se seu conjunto de dados for tão grande, eu acho que uma abordagem de buffer baseada em disco seria melhor:

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

Eu também experimentaria agrupar um número maior de buckets, por exemplo, se sua string fosse:

GATTACA

a primeira chamada MSB retornaria o bucket para o GATT (256 total de buckets), dessa forma você cria menos ramificações do buffer baseado em disco. Isso pode ou não melhorar o desempenho, então experimente.

Vou falar sobre um membro e sugerir que você mude para uma implementação heap / heapsort . Esta sugestão vem com algumas suposições:

1. Você controla a leitura dos dados
2. Você pode fazer algo significativo com os dados classificados assim que começar a classificá-los.

A vantagem do heap / heap-sort é que você pode criar o heap enquanto lê os dados e pode começar a obter resultados no momento em que criou o heap.

Vamos voltar. Se você tiver a sorte de poder ler os dados de forma assíncrona (ou seja, poderá postar algum tipo de solicitação de leitura e ser notificado quando alguns dados estiverem prontos), poderá criar um pedaço do heap enquanto aguarda o próximo pedaço de dados a entrar - mesmo a partir do disco. Geralmente, essa abordagem pode ocultar a maior parte do custo de metade da sua classificação, atrás do tempo gasto na obtenção dos dados.

Depois de ler os dados, o primeiro elemento já estará disponível. Dependendo de onde você está enviando os dados, isso pode ser ótimo. Se você estiver enviando para outro leitor assíncrono, ou algum modelo paralelo de 'evento' ou interface do usuário, poderá enviar trechos e trechos à medida que avança.

Dito isto - se você não tem controle sobre como os dados são lidos e lidos de forma síncrona, e você não usa os dados classificados até que sejam totalmente gravados - ignore tudo isso. :(

Veja os artigos da Wikipedia:

• Heapsort
• Pilha binária

A " classificação Radix sem espaço extra " é um documento que trata do seu problema.

Em termos de desempenho, convém procurar algoritmos mais gerais de classificação de comparação de cadeias.

Atualmente você acaba tocando todos os elementos de cada corda, mas pode fazer melhor!

Em particular, uma classificação de intermitência é um ajuste muito bom para este caso. Como bônus, como o burstsort é baseado em tentativas, ele funciona ridiculamente bem para os pequenos tamanhos de alfabeto usados ​​no DNA / RNA, já que você não precisa criar nenhum tipo de nó de pesquisa ternário, hash ou outro esquema de compressão de nós trie no diretório trie implementação. As tentativas também podem ser úteis para seu objetivo final do tipo sufixo-array.

Uma implementação decente de propósito geral do burstsort está disponível no source forge em http://sourceforge.net/projects/burstsort/ - mas não está no local.

Para fins de comparação, a implementação do C-burstsort coberta em http://www.cs.mu.oz.au/~rsinha/papers/SinhaRingZobel-2006.pdf faz um benchmark 4-5x mais rápido que o quicksort e radix para algumas cargas de trabalho típicas.

Você vai querer dar uma olhada no Processamento de Sequência Genômica em Grande Escala pelos Drs. Kasahara e Morishita.

As seqüências de caracteres compostas pelas quatro letras de nucleotídeo A, C, G e T podem ser especialmente codificadas em números inteiros para um processamento muito mais rápido. A classificação Radix está entre muitos algoritmos discutidos no livro; você poderá adaptar a resposta aceita a esta pergunta e obter uma grande melhoria no desempenho.

Você pode tentar usar um trie . A classificação dos dados é simplesmente iterativa no conjunto de dados e inserida; a estrutura é classificada naturalmente e você pode pensar nela como semelhante a uma árvore B (exceto que, em vez de fazer comparações, você sempre usa indiretos de ponteiro).

O comportamento de armazenamento em cache favorecerá todos os nós internos; portanto, você provavelmente não melhorará isso; mas você também pode mexer com o fator de ramificação de sua tentativa (verifique se todos os nós se encaixam em uma única linha de cache, aloque os nós semelhantes a um heap, como uma matriz contígua que representa um percurso de ordem de nível). Como as tentativas também são estruturas digitais (O (k) inserir / localizar / excluir para elementos de comprimento k), você deve ter desempenho competitivo em uma classificação de raiz.

Eu explodiria uma representação compactada de bits das strings. Alega-se que o Burstsort tem uma localização muito melhor do que as classificações de raiz, mantendo o uso de espaço extra baixo com tentativas de rajada no lugar das tentativas clássicas. O papel original tem medidas.

O Radix-Sort não é consciente do cache e não é o algoritmo de classificação mais rápido para grandes conjuntos. Você pode olhar para:

• ti7qsort . ti7qsort é a ordenação mais rápida para números inteiros (pode ser usada para cadeias de tamanho pequeno e fixo).
• Inline QSORT
• Classificação de cadeias

Você também pode usar a compactação e codificar cada letra do seu DNA em 2 bits antes de armazenar na matriz de classificação.

A classificação de raíz de MSB do dsimcha parece boa, mas Nils se aproxima do cerne do problema com a observação de que a localidade do cache é o que está matando você em grandes tamanhos de problemas.

Sugiro uma abordagem muito simples:

1. Estime empiricamente o maior tamanho mpara o qual uma classificação de base é eficiente.
2. Leia blocos de melementos por vez, classifique-os rapidamente e escreva-os (em um buffer de memória, se você tiver memória suficiente, mas arquivar), até esgotar sua entrada.
3. Mescla os blocos classificados resultantes.

O Mergesort é o algoritmo de classificação mais amigável ao cache que eu conheço: "Leia o próximo item da matriz A ou B e, em seguida, escreva um item no buffer de saída". Ele roda eficientemente em unidades de fita . Ele requer 2nespaço para classificar nitens, mas minha aposta é que a localidade de cache muito aprimorada que você verá tornará isso sem importância - e se você estivesse usando uma classificação de radix não no local, precisaria desse espaço extra.

Finalmente, observe que o mergesort pode ser implementado sem recursão e, de fato, fazê-lo dessa maneira deixa claro o verdadeiro padrão de acesso linear à memória.

Parece que você resolveu o problema, mas, para o registro, parece que uma versão de uma classificação de raiz no local viável é a "Classificação da Bandeira Americana". Está descrito aqui: Engineering Radix Sort . A idéia geral é fazer 2 passes em cada caractere - primeiro conte quantos de cada um você tem, para poder subdividir a matriz de entrada em posições. Em seguida, prossiga novamente, trocando cada elemento na bandeja correta. Agora classifique recursivamente cada posição na próxima posição de caractere.

Primeiro, pense na codificação do seu problema. Livre-se das strings, substitua-as por uma representação binária. Use o primeiro byte para indicar comprimento + codificação. Como alternativa, use uma representação de comprimento fixo em um limite de quatro bytes. Então a classificação do radical se torna muito mais fácil. Para uma classificação de base, o mais importante é não ter tratamento de exceção no ponto quente do loop interno.

OK, pensei um pouco mais sobre o problema quádruplo. Você quer uma solução como uma árvore Judy para isso. A próxima solução pode lidar com cadeias de comprimento variável; para comprimento fixo, remova os bits de comprimento, o que realmente facilita.

Aloque blocos de 16 ponteiros. O bit menos significativo dos ponteiros pode ser reutilizado, pois seus blocos sempre estarão alinhados. Você pode querer um alocador de armazenamento especial para ele (dividindo o armazenamento grande em blocos menores). Existem vários tipos diferentes de blocos:

• Codificação com 7 bits de comprimento de cadeias de comprimento variável. À medida que são preenchidos, você os substitui por:
• A posição codifica os próximos dois caracteres, você tem 16 ponteiros para os próximos blocos, terminando com:
• Codificação de bitmap dos três últimos caracteres de uma sequência.

Para cada tipo de bloco, você precisa armazenar informações diferentes nos LSBs. Como você possui cadeias de comprimento variável, também é necessário armazenar o fim da cadeia, e o último tipo de bloco só pode ser usado para as cadeias mais longas. Os 7 bits de comprimento devem ser substituídos por menos à medida que você se aprofunda na estrutura.

Isso fornece um armazenamento razoavelmente rápido e com muita memória eficiente de seqüências classificadas. Ele se comportará um pouco como um trie . Para fazer isso funcionar, certifique-se de criar testes de unidade suficientes. Você deseja cobertura de todas as transições de bloco. Você deseja começar apenas com o segundo tipo de bloco.

Para obter ainda mais desempenho, convém adicionar diferentes tipos de bloco e um tamanho maior de bloco. Se os blocos sempre tiverem o mesmo tamanho e forem grandes o suficiente, você poderá usar ainda menos bits para os ponteiros. Com um tamanho de bloco de 16 ponteiros, você já tem um byte livre em um espaço de endereço de 32 bits. Veja a documentação da árvore Judy para tipos de blocos interessantes. Basicamente, você adiciona código e tempo de engenharia para uma troca de espaço (e tempo de execução)

Você provavelmente deseja começar com um radix direto de 256 de largura para os quatro primeiros caracteres. Isso fornece uma troca decente de espaço / tempo. Nesta implementação, você obtém muito menos sobrecarga de memória do que com uma simples tentativa; é aproximadamente três vezes menor (não medi). O (n) não é problema se a constante for baixa o suficiente, como você observou ao comparar com o quicksort O (n log n).

Você está interessado em lidar com duplas? Com sequências curtas, haverá. Adaptar os blocos para lidar com contagens é complicado, mas pode ser muito eficiente em termos de espaço.