Pergunta

Este é um texto longo. Por favor, tenha paciência comigo. Resumia, a pergunta é:? Existe uma viável no local radix sort algoritmo

Preliminar

Eu tenho um número enorme de pequena de comprimento fixo cordas que usam apenas as letras “A”, “C”, “G” e “T” (sim, você adivinhou -lo:. DNA ) que eu quero tipo

No momento, eu uso std::sort que usa introsort em todas as implementações comuns do STL . Isso funciona muito bem. No entanto, estou convencido de que radix encaixa tipo meu conjunto problema perfeitamente e deve funcionar muito melhor na prática.

Detalhes

Eu testei esta suposição com uma implementação muito ingênua e relativamente pequenos entradas (na ordem de 10.000) isso era verdade (bem, pelo menos duas vezes mais rápido). No entanto, em tempo de execução degrada abismalmente quando o tamanho do problema se torna maior ( N > 5.000.000).

A razão é óbvia: Radix tipo requer copiar os dados inteiros (mais de uma vez na minha implementação ingênua, na verdade). Isso significa que eu coloquei ~ 4 GiB em minha principal memória que, obviamente, mata desempenho. Mesmo que isso não aconteceu, eu não posso dar ao luxo de usar essa quantidade de memória já que o problema tamanhos realmente tornar-se ainda maior.

Casos de Uso

Idealmente, esse algoritmo deve funcionar com qualquer comprimento de cadeia de entre 2 e 100, por DNA, bem como DNA5 (que permite que um caractere curinga adicional “N”), ou DNA, mesmo com IUPAC ambigüidade códigos (resultando em 16 valores distintos). No entanto, percebo que todos estes casos não pode ser coberto, por isso estou feliz com qualquer melhoria de velocidade que recebo. O código pode decidir dinamicamente qual algoritmo para expedição para.

Research

Infelizmente, o Wikipedia sobre radix sort é inútil. A seção sobre uma variante no local é um lixo completo. A seção NIST-DADS em radix sort fica ao lado inexistente. Há um papel prometendo-som chamado Eficiente Adaptive In-Place Radix Sorting que descreve o algoritmo “MSL”. Infelizmente, este trabalho, também, é decepcionante.

Em particular, existem as seguintes coisas.

Primeiro, o algoritmo contém vários erros e deixa muito inexplicada. Em particular, ele não detalha a chamada recursão (I simplesmente assumir que ele aumenta ou reduz alguns ponteiro para calcular a mudança atual e mascarar valores). Além disso, ele usa o funções dest_group e dest_address sem dar definições. Eu não consigo ver como implementar las de forma eficiente (isto é, em O (1), pelo menos dest_address não é trivial)

.

Por último, mas não menos importante, os alcança algoritmo in-place-ness trocando índices de array com elementos dentro da matriz de entrada. Isto, obviamente, só funciona em matrizes numéricas. Eu preciso usá-lo nas cordas. Claro, eu poderia apenas parafuso tipagem forte e ir em frente assumindo que a memória vai tolerar meu armazenar um índice onde ele não pertence. Mas isso só funciona, desde que eu posso espremer minhas cordas em 32 bits de memória (supondo que 32 bit inteiros). Isso é apenas 16 caracteres (vamos ignorar por um momento que 16> log (5.000.000)).

Outro papel por um dos autores dá nenhuma descrição precisa em tudo, mas giruntime do ves MSL como sub-linear que é totalmente errado.

Para recapitular :? Existe alguma esperança de encontrar uma implementação de referência de trabalho ou pelo menos uma boa pseudocódigo / descrição de uma raiz de trabalho no local tipo que funciona em cadeias de DNA

Foi útil?

Solução

Bem, aqui está uma implementação simples de um radix MSD tipo de DNA. Está escrito em D porque essa é a linguagem que eu uso mais e, portanto, sou menos propensos a cometer erros bobos em, mas pode ser facilmente traduzido para outra língua. É no local, mas requer 2 * seq.length passa através da 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, este é um tipo de específico para DNA, ao contrário de ser geral, mas deve ser rápido.

Editar:

fiquei curioso se este código realmente funciona, então eu testei / depurado-lo enquanto espera para o meu próprio código de bioinformática para ser executado. A versão anterior agora é realmente testado e funciona. Por 10 milhões de sequências de 5 bases cada, é sobre 3x mais rápido que um introsort otimizado.

Outras dicas

Eu nunca vi uma raiz no local espécie, e da natureza do radix-sort eu duvido que ele é muito mais rápido do que um fora de lugar espécie, desde que os ataques matriz temporária para a memória.

Motivo:

A classificação faz uma leitura linear na matriz de entrada, mas todas as gravações serão quase aleatória. De um certo N para cima isto resume-se a um erro de cache por gravação. Este erro de cache é o que retarda o seu algoritmo. Se é no lugar ou não, não vai mudar este efeito.

Eu sei que isso não vai responder a sua pergunta diretamente, mas se a classificação é um gargalo que você pode querer ter um olhar para perto triagem algoritmos como uma etapa de pré-processamento ( o wiki páginas sobre o soft-heap pode começar).

Isso poderia dar um impulso muito agradável de cache localidade. Um livro-texto fora-de-lugar radix sort, então, um melhor desempenho. As gravações ainda será quase aleatória, mas pelo menos eles vão agrupar em torno dos mesmos pedaços de memória e, como tal, aumentar o cache taxa de acertos.

Eu não tenho idéia se ele funciona na prática embora.

BTW: Se você está lidando apenas com cordas de DNA: Você pode compactar um char em dois pedaços e embalar seus dados bastante. Isso irá reduzir o requisito de memória pelo fator de quatro sobre uma representação naiive. Dirigindo-se torna mais complexo, mas a ALU de sua CPU tem muito tempo para gastar durante todo o esconderijo-acidentes de qualquer maneira.

Você pode certamente deixar cair os requisitos de memória ao codificar a seqüência em bits. Você está navegando na permutações assim, para o comprimento 2, com "ACGT" que é 16 estados, ou 4 bits. Para o comprimento de 3, que é 64 estados, os quais podem ser codificados em 6 bits. Portanto, parece que 2 bits para cada letra na seqüência, ou cerca de 32 bits para 16 caracteres como você disse.

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

Assim, para seqüências de comprimento 3, pode-se criar 64 baldes, talvez uint32 porte, ou uint64. Inicializar-los para zero. Iterate através de sua lista muito grande de 3 seqüências de char, e codificá-los como acima. Use isso como um subscrito e incremento que balde.
Repita este procedimento até todas as suas sequências foram processadas.

Em seguida, regenerar sua lista.

percorrer os baldes 64 em ordem, para a contagem encontrados em que balde, gerar que muitas instâncias da sequência representada por esse balde.
quando todos os baldes foram iterado, você tem sua matriz classificada.

Uma sequência de quatro, adiciona 2 bits, de modo que não haveria 256 baldes. Uma seqüência de 5, adiciona 2 bits, por isso não seria de 1024 baldes.

Em algum momento o número de baldes vai aproximar seus limites. Se você ler as seqüências de um arquivo, em vez de mantê-los na memória, mais memória estaria disponível para baldes.

Eu acho que isso seria mais rápido do que fazendo o tipo in situ como os baldes são susceptíveis de caber dentro de seu conjunto de trabalho.

Aqui está um corte 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 o seu conjunto de dados é tão grande, então eu acho que uma abordagem tampão baseado 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

Também gostaria de experimentar agrupamento em um maior número de baldes, por exemplo, se a seqüência foi: 

GATTACA

a primeira chamada MSB iria devolver o balde para GATT (256 baldes totais), de que maneira você faz menos galhos do buffer baseado em disco. Isto pode ou não pode melhorar o desempenho, então experimento com ele.

Vou sair em um membro e sugiro que você mudar para uma pilha / heapsort implementação. 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 você 'Iniciar' ficando classificados.

A beleza da pilha / pilha-sort é que você pode construir a pilha enquanto você lê os dados, e você pode começar a obter resultados no momento em que você construiu o heap.

passo para trás Vamos. Se você está tão feliz que você pode ler os dados de forma assíncrona (isto é, você pode enviar algum tipo de pedido de leitura e ser notificado quando alguns dados estão prontos), e então você pode construir um pedaço da pilha enquanto você está esperando para o próximo bloco de dados para entrar - mesmo a partir do disco. Muitas vezes, esta abordagem pode enterrar a maior parte do custo de metade da sua ordenação por trás o tempo gasto recebendo os dados.

Depois de ter os dados lidos, o primeiro elemento já está disponível. Dependendo de onde você está enviando os dados, isso pode ser grande. Se você está enviando-a para um outro leitor assíncrona, ou algum modelo paralelo 'evento', ou UI, você pode enviar pedaços e pedaços como você vai.

Dito isso - se você não tem controle sobre como os dados são lidos, e ele é lido de forma síncrona, e você não tem nenhum uso para os dados ordenados até que seja inteiramente escrito - ignorar tudo isso. : (

Ver artigos da Wikipédia:

Em termos de performance que você pode querer olhar para uma seqüência de comparação mais geral algoritmos de ordenação.

Atualmente você acaba tocando cada elemento de cada corda, mas pode fazer melhor!

Em particular, um burstsort é um ajuste muito bom para este caso. Como um bônus, desde burstsort é baseada em tentativas, ele funciona ridiculamente bem para o pequeno alfabeto tamanhos utilizados no DNA / RNA, desde que você não precisa para construir qualquer tipo de nó ternário pesquisa, hash ou outro esquema de compressão nó trie para o implementação trie. As tentativas pode ser útil para o seu objetivo sufixo-array-like final, bem.

A implementação de propósito geral decente de burstsort está disponível na forja fonte em http://sourceforge.net/projects / burstsort / -. mas não é no local

Para fins de comparação, a implementação C-burstsort coberto a http://www.cs.mu.oz.au/~rsinha/papers/SinhaRingZobel-2006.pdf benchmarks 4-5x mais rápido do que quicksort e Radix sortes para algumas cargas de trabalho típicos.

Você vai querer dar uma olhada grande escala Genoma seqüência de processamento pelos Drs. Kasahara e Morishita.

cordas compostas por as quatro letras de nucleótidos A, C, G, e T podem ser codificados especialmente em números inteiros para muito processamento mais rápido. Radix espécie está entre muitos algoritmos discutidos no livro; você deve ser capaz de se adaptar a resposta aceita a esta questão e ver uma melhoria grande desempenho.

" Radix classificação sem espaço extra" é um documento abordando o seu problema.

Você pode tentar usar um trie . A classificação dos dados é simplesmente a iteração através do conjunto de dados e inseri-lo; a estrutura é naturalmente ordenada, e você pode pensar nisso como semelhante a um B-Tree (exceto em vez de fazer comparações, você sempre uso indireções ponteiro).

Cache comportamento favorecerá todos os nós internos, então você provavelmente não vai melhorar a situação; mas você pode mexer com o fator de ramificação de sua trie bem (garantir que cada nó se encaixa em uma única linha de cache, alocar trie nodos semelhante a uma pilha, como uma matriz contígua que representa um percurso nível ordem). Desde tentativas também são estruturas digitais (O (k) insert / encontrar / delete para elementos de comprimento k), você deve ter desempenho competitivo a um radix tipo.

burstsort uma representação das cordas bits embalado. Burstsort é reivindicado ter muito melhor localidade de Radix Sort, mantendo o baixo uso de espaço extra com tentativas de ruptura no lugar de tentativas clássicas. O documento original tem medições.

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

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

MSB do dsimcha Radix tipo parece bom, mas Nils fica mais perto do coração do problema com a observação de que o cache local é o que está matando você em tamanhos grandes problemas.

Eu sugiro uma abordagem muito simples:

  1. Empiricamente estimar o maior m tamanho para o qual uma espécie radix é eficiente.
  2. Leia blocos de elementos m de cada vez, radix sort-los e escrevê-los (para um buffer de memória, se você tem memória suficiente, mas de outra forma para o arquivo), até esgotar a sua entrada.
  3. Mergesort os blocos classificados resultantes.

Mergesort é o algoritmo de classificação mais amigável-cache Estou ciente: "Leia o próximo item a partir de qualquer conjunto A ou B, em seguida, escrever um artigo para o buffer de saída." Ele é executado de forma eficiente em unidades de fita . Ele requer espaço 2n a itens tipo n, mas minha aposta é que o muito melhorada de cache localidade você verá vai fazer isso sem importância - e se você estivesse usando uma raiz não-in-place tipo, você precisava esse espaço extra de qualquer maneira.

Por favor, note finalmente que mergesort pode ser implementado sem recursão, e de fato fazê-lo desta maneira deixa claro o verdadeiro padrão de acesso de memória linear.

Parece que você já resolveu o problema, mas para o registro, parece que uma versão de um radix viável no local espécie é a "bandeira americana Sort". É descrito aqui: Engenharia Radix Sort . A ideia geral é a de fazer 2 passes em cada personagem - primeiro contar quantos de cada você tem, para que possa subdividir a matriz de entrada em caixas. Em seguida, passar de novo, trocando cada elemento para o lixo correta. Agora recursivamente tipo cada bin na próxima posição de caractere.

Primeiro, pense sobre a codificação do seu problema. Se livrar das cordas, substituí-los por uma representação binária. Usar o primeiro byte para indicar + codificação de comprimento. Alternativamente, usam uma representação de comprimento fixado em um limite de quatro bytes. Em seguida, a radix sort torna-se muito mais fácil. Para um radix sort, a coisa mais importante é não ter manipulação de exceção no ponto quente do circuito interno.

OK, eu pensei um pouco mais sobre o problema 4-nary. Você quer uma solução como um Judy árvore para isso. A solução seguinte pode lidar com cadeias de comprimento variável; para comprimento fixo basta remover os bits de comprimento, que realmente faz com que seja mais fácil.

Alocar blocos de 16 ponteiros. O bit menos significativo dos ponteiros podem ser reutilizados, como seus blocos ficarão sempre alinhados. Você pode querer um alocador de armazenamento especial para ele (quebrando grande armazenamento em blocos menores). Há um número de diferentes tipos de blocos:

  • Encoding com 7 bits de comprimento de cadeias de comprimento variável. Como eles encher-se, você substituí-los por:
  • Posição codifica os próximos dois personagens, você tem 16 ponteiros para os próximos blocos, terminando com:
  • Bitmap codificação dos três últimos caracteres de uma string.

Para cada tipo de bloco, você precisa armazenar informações diferentes nas LSB. Como você tem cadeias de comprimento variável que você precisa para armazenar fim-de-corda também, e último tipo de bloco só pode ser usado para as cordas mais longas. Os 7 bits de comprimento deve ser substituído pelo menos como você se aprofundar na estrutura.

Isto proporciona-lhe uma memória de armazenamento de strings ordenadas razoavelmente rápido e muito eficiente. Ele irá se comportar um pouco como um trie . Para começar este trabalho, certifique-se de construir testes de unidade suficientes. Você quer que a cobertura de todas as transições de bloco. Você quer começar com apenas o segundo tipo de bloco.

Para ainda mais o desempenho, você pode querer adicionar diferentes tipos de blocos e um tamanho maior do bloco. Se os blocos são sempre do mesmo tamanho e grande o suficiente, você pode usar até mesmo menos bits para os ponteiros. Com um tamanho de bloco de 16 pontos, você já tem um livre byte em um espaço de endereçamento de 32 bits. Dê uma olhada na documentação árvore Judy para tipos de blocos interessantes. Basicamente, você adicionar código e engenharia tempo para um espaço (e tempo de execução) trade-off

Você provavelmente vai querer começar com um 256 radix ampla direta para os primeiros quatro caracteres. Que fornece uma decente tradeoff espaço / tempo. Nesta implementação, você recebe muito menos sobrecarga de memória do que com um simples trie; é cerca de três vezes menor (I não medido). O (n) não é problema se a constante é baixo o suficiente, como você notou quando se compara com o O (n log n) quicksort.

Você está interessado em manipulação de duplas? Com seqüências curtas, não vão ser. Adaptar os blocos para contagem alça é complicado, mas pode ser muito eficiente de espaço.

Licenciado em: CC-BY-SA com atribuição
Não afiliado a StackOverflow
scroll top