Mergesort em C++ usando o pacote Rcpp

A muito tempo atrás eu comecei a ver como dava para usar código em linguagem C do R neste post aqui. Depois disso, no comentário alguém mais esperto que eu falou sobre um tal de pacote Rcpp, que permitia de maneira mais simples, usar código em linguagem C++ do R, depois disso eu fiz mais alguns testes, tentando implementar alguns algorítimos de ordenação, o bubblesort aqui e o insertion sort aqui.

Esse é o tipo de coisa inútil de se fazer, mas é divertido. E esses algorítimos de ordenação tem em qualquer curso introdutório de programação, então são implementações básicas para qualquer linguagem. Depois do mergesort eu ainda tenho o quicksort, mas no mergesort tudo começa a ficar mais complexo.

Por curiosidade, o método padrão na função sort() do R é o Shellsort, criado por Donald Shell mas o que está implementado no R é uma variante do Robert Sedgewick, esse cara da quatro cursos sobre programação gratuitos no Coursera, so olhar aqui.

Mas na função sort do R também temos o Quicksort, mas ele vai ficar pra próxima.

Agora vamos ver o mergesort. Ele é um algorítimo recursivo, ou seja uma função que chama ela mesma, e a ideia básica é, melhor que ordenar todo um vetor, é pegar dois vetores ordenados e então simplesmente intercalar eles. Mas a gente precisa de dois vetores ordenados pra poder usar esse esquema certo?
Pra isso é so a gente quebrar o vetor no meio, e de novo no meio para as duas metades, e de novo para os quatro quartos, e de novo até que vai chegar um momento que so vamos ter mini vetores de um elemento, um vetor de um elemento está necessariamente ordenado, logo agora é so voltar intercalando esses vetores que fica tudo bem.

MergeSort

É até estranho como isso funciona, mas funciona. Bem o legal que diferente do InsertionSort e o BubbleSort, aqui a gente vai precisar de uma função auxiliar, então a gente vai escrever o código em um arquivo, que vamos chamar de merge.cpp, e depois vamos chamar esse arquivo do R usando a função sourceCpp() do pacote Rcpp.

Então o código do mergesort é assim:

#include <Rcpp.h>
using namespace Rcpp;
 
void intercala(int p, int q, int r, NumericVector v,  NumericVector w)
{
  int i, j, k;
   i = p;
   j = q;
   k = 0;
   while (i < q && j < r) {
      if (v[i] < v[j]) {
         w[k] = v[i];
         i++;
      }
      else {
         w[k] = v[j];
         j++;
      }
      k++;
   }
   while (i < q) {
      w[k] = v[i];
      i++;
      k++;
   }
   while (j < r) {
      w[k] = v[j];
      j++;
      k++;
   }
   for (i = p; i < r; i++)
      v[i] = w[i-p];
}
 
void mergesort(int p, int r, NumericVector v, NumericVector aux)
{
   int q;
   if (p < r - 1) {
      q = (p + r) / 2;
      mergesort(p, q, v,aux);
      mergesort(q, r, v,aux);
      intercala(p, q, r, v,aux);
   }
}
 
// [[Rcpp::export]]
NumericVector mergesortC(NumericVector vetor) {
  Rcpp::NumericVector saida = Rcpp::clone(vetor);
  Rcpp::NumericVector aux = Rcpp::clone(vetor);
  int n = saida.size();
  mergesort(0,n,saida,aux);
  return saida;
}

Então, tudo isso está dentro de um arquivo que eu chamei de merge.cpp
Basicamente o código você acha em qualquer esquina da rede mundial de computadores, mas especificamente aqui, temos que importar a biblioteca Rcpp.h, já que estamos usando um objeto específico dela, que é o NumericVector. Outra coisa é que o R manda um ponteiro pra cá, que aponta aonde estão os dados, desse objeto. Se a gente tentar mexer direto nele, não vai dar certo, eu tentei e fiz o R travar com isso heheh, ai perguntando no forum aqui, o cara que criou o pacote me mandou ir estudar, ai lendo o livro dele eu vi que mexer direto nos dados, muitas vezes não é uma boa ideia, dai a gente usa o método clone de NumericVector, a gente clona ele para outro vetor, e então manipula esse novo vetor, pra não mexer na memoria que não devemos.
Outra coisa é que, talvez um ponto fraco do mergeSort, é que para intercalar dois vetores, precisamos de outro lugar para ir guardando eles em ordem, ou seja, diferente do insertionsort, não da para ir arrumando dentro do próprio vetor, por isso precisamos de um vetor auxiliar, aqui eu ja criei um vetor auxiliar chamado aux, para ser usado para esse propósito. Eu até tentei criar vetores dentro da função intercala, mas isso deixou tudo bem lento, desse jeito a velocidade ficou um pouco melhor. Por último, mas bem legal, é que nesse caso, veja que temos três funções, uma que recebe e retorna tudo para o R, e outras duas que fazem o trabalho propriamente dito. Mas quando a gente olha no R, desse jeito a gente só vai ver a mergesortC, isso porque em cima dela tem escrito

// [[Rcpp::export]]

Isso faz com que a gente só veja ela, e não as funções auxiliares.

Concluído a construção da função em C++, agora é simples no R

> library(Rcpp) > sourceCpp("merge.cpp") > mergesortC(sample(1:10)) [1] 1 2 3 4 5 6 7 8 9 10

Basicamente, só damos um sourceCpp no arquivo com o código, podemos usar o argumento verbose=TRUE se quisermos acompanhar a compilação e ta tudo pronto para usar, e podemos ver que o mergesortC funciona certinho.

Agora a gente pode até medir o tempo que ele precisa para organizar, por exemplo, 1000 números inteiros, usando a função system.time.

> system.time(mergesortC(sample(1:1000))) usuário sistema decorrido 0.000 0.000 0.001

Alias, para ter uma noção de porque as pessoas continuaram inventando algorítimos de ordenação, a gente pode fazer uma comparação, entre o mergesort e o insertionsort por exemplo, e ver qual deles é o mais rápido, para organizar um vetor de um tamanho qualquer, é so a gente ficar bagunçando um vetor, com sample, e registar o tempo que ele leva para organizar, e repetir isso varias vezes.

dados<-data.frame(Tempo=rep(NA,1000),Algoritimo=rep(c("InsertionSort","MergeSort"),each=500)) for(i in 1:1000){ if(i<=500){ dados[i,1]<-system.time(insertionsortC(sample(1:10000)))[3] }else{ dados[i,1]<-system.time(mergesortC(sample(1:10000)))[3] } print(i) }

A gente cria então um dataframe para receber os dados, e faz o experimento, pegando o tempo total gasto no processo, o insertionsortC eu peguei do post passado.

O resultado disso é o seguinte.

01

É, o mergesort pode ser mais complicado de entender, e implementar, mas é muito mais rápido.

> aggregate(dados$Tempo, list(dados$Algoritimo), mean) Group.1 x 1 InsertionSort 0.043494 2 MergeSort 0.005362 > aggregate(dados$Tempo, list(dados$Algoritimo), sd) Group.1 x 1 InsertionSort 0.0044236978 2 MergeSort 0.0004933993

Ele é muitas vezes mais rápido, e varia muito menos no tempo de execução, não é a toa que o sort do R não usa o InsertionSort, mas ele também não usa o mergesort, ou seja, tem como melhorar ainda, mas essa vai ficar para a próxima porque ja está na hora de dormir.

Bem o script está aqui embaixo além do repositório recologia no github. Esse foi o primeiro código no Rcpp que envolvia mais de uma função, e no final muitas coisas envolvem varias funções, então acho que assim vamos abrindo portas para possibilidades legais. Outra coisa é que esse tipo de estratégia para resolver problemas, de partir o problema em partes menores, e resolver essas partes antes de resolver um problema maior, pode parecer algo longe para biologia, mas a gente vai pela mesma estrada em problemas como o de alinhamento de sequências. E vamos dormir.

library(Rcpp)
 
sourceCpp("merge.cpp",verbose=T)
 
mergesortC(sample(1:10))
 
system.time(mergesortC(sample(1:1000)))
 
cppFunction("
    NumericVector insertionsortC(NumericVector vetor) {
        int n = vetor.size();
 
        double aux;
        int i , j;
 
        for(i=1;i<n;i++) {
            aux=vetor[i];
            j=i-1;
            while(j>=0 && vetor[j]>aux) {
                vetor[j+1]=vetor[j];
                j=j-1;
                }
            vetor[j+1]=aux;
            }
        return vetor;
        }
")
 
dados<-data.frame(Tempo=rep(NA,1000),Algoritimo=rep(c("InsertionSort","MergeSort"),each=500))
 
for(i in 1:1000){
    if(i<=500){
        dados[i,1]<-system.time(insertionsortC(sample(1:10000)))[3]
    }else{
        dados[i,1]<-system.time(mergesortC(sample(1:10000)))[3]
    }
    print(i)
}
 
#Figura 1
boxplot(Tempo~Algoritimo,data=dados)
 
aggregate(dados$Tempo, list(dados$Algoritimo), mean)
aggregate(dados$Tempo, list(dados$Algoritimo), sd)