Projeto Final
Heitor Victor
16/06/2018
O objetivo deste projeto é otimizar uma função criada para randomizar elementos de uma matriz exparsa, simétrica e com zeros na diagonal prinicipal. Além disso, as matrizes de interesse só possuem valores 1 ou 0. Para entender melhor, suponha que você possua uma matriz desta forma:
\[ \begin{bmatrix} 0 & 1 & 0 & 1 & 0\\ 1 & 0 & 1 & 0 & 1\\ 0 & 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 0 \end{bmatrix} \]
Imagine agora que por algum motivo específico você queira trocar os elementos desta matriz de lugar, mas que ela continue com as propriedades citadas. Este procedimento é muito útil por exemplo, quando a matriz citada é uma matriz de vizinhança e queremos saber se a dependência entre os vizinhos continua caso façamos uma randomização na mesma.
Com esta motivação um código em R foi escrito, contudo, ele deixa a desejar quanto ao desempenho. Daí, através do pacote Rcpp algumas funções que compoem o código em R foram reescritas em C++.
O código possui 2 funções “chave”, que são:
sorteador: Esta função sorteia os elementos a serem trocados de lugar. Possui como entrada 2 vetores, que representam as posições da linha (i) e coluna (j) a serem trocadas.
Funciona da seguinte maneira: Enquanto o valor de \(i \leq j\), o sorteio é refeito. Isto é necessário já que a matriz é simétrica.
taboo: Esta função guarda um número “n” das posições que foram trocadas durante o processo dentro de uma matriz cujo os 2 primeiros elementos são a posição (i,j) de um elemento e os outros dois são as posições (i,j) do outro elemento, ambos sorteados pela função sorteador. A matriz desta função é normalmente chamada na literatura de “Matriz taboo”.
Note que a quantidade de linhas da matriz taboo dada pelo usúario nunca pode ser maior que o número de elementos a serem trocados,caso isso aconteça, o algoritmo entrará em “loop infinito”.
A seguir um exemplo para fixação do modo de uso da função usando a matriz citada anteriormente e seu respectivo código em R.
#Funções e Variaveis Globais
#OBS: Compilar na ordem que está disposto neste script
##########################################################
o.tab = function(e1,e2,tab,n.taboo)
{
flag = -1
for(i in 1:n.taboo)
{
if(anyNA(tab[i,1])==TRUE)
{
break
}
if(((e1[1]==tab[i,1])&&(e1[2]==tab[i,2]))||((e1[1]==tab[i,3])&&(e1[2]==tab[i,4]))||
((e2[1]==tab[i,1])&&(e2[2]==tab[i,2]))||((e2[1]==tab[i,3])&&(e2[2]==tab[i,4])))
{
flag = 0
break
}
else
{
flag = 1
}
}
return (flag)
}
###########################################################
#Sorteador das posicoes da matriz
sorteador = function(In,Jn)
{
j = i = NULL;
aux = numeric(2)
#ij eh a posicao sorteada para a iteração
i = sample(In, 1, replace = T)
j = sample(Jn, 1, replace = T)
while(i<=j)
{
i = sample(In, 1, replace = T)
j = sample(Jn, 1, replace = T)
}
aux = c(i,j)
return(aux)
}
############################################################
smr = function(matriz.original, n = 1000, n.taboo = 5)
{
tab = matrix(nrow = n.taboo ,ncol = 4)
m = matriz.original
m.or = matriz.original
aux1 = aux2 = e1 = e2 = flag = NULL # e1 e e2 sao as posições dos elementos de troca
o = t = ta = 0 #Contador de iterações do while
In = c(1:length(m.or[,1])) #range para o nº de linhas
Jn = c(1:length(m.or[1,])) #range para o nº de colunas
for(i in 1:n)
{
e1 = sorteador(In,Jn)
e2 = sorteador(In,Jn)
while(m.or[e1[1],e1[2]]==m.or[e2[1],e2[2]])
{
e1 = sorteador(In,Jn)
e2 = sorteador(In,Jn)
o = o + 1
}
if(t<n.taboo)
{t = t +1}
else
{t = 0}
tab[t,1] = e1[1]
tab[t,2] = e1[2]
tab[t,3] = e2[1]
tab[t,4] = e2[2]
if(n.taboo>0){
flag = o.tab(e1,e2,tab,n.taboo) #Sinalizador se o indice sorteado pertence a matriz taboo
while(flag < 1||m.or[e1[1],e1[2]]==m.or[e2[1],e2[2]])
{
e1 = sorteador(In,Jn)
e2 = sorteador(In,Jn)
flag = o.tab(e1,e2,tab,n.taboo)
ta = ta + 1
}
}
tab[t,1] = e1[1]
tab[t,2] = e1[2]
tab[t,3] = e2[1]
tab[t,4] = e2[2]
aux1 = m.or[e1[1],e1[2]]
aux2 = m.or[e1[2],e1[1]]
m.or[e1[1],e1[2]] = m.or[e2[1],e2[2]]
m.or[e1[2],e1[1]] = m.or[e2[2],e2[1]]
m.or[e2[1],e2[2]] = aux1
m.or[e2[2],e2[1]] = aux2
}
resu = list(m,m.or,tab,o,ta)
names(resu) = c("Matriz.Original","Matriz.Randomizada","Matriz.Taboo.n.5",
"N.de.vezes.elementos.iguais","Taboo.vezes.entrada")
return(resu)
}
mt <- mt <- cbind(c(0,1,0,1,0),c(1,0,1,0,1),c(0,1,0,0,0),c(1,0,0,0,0),c(0,1,0,0,0))
# matriz teste
smr(# matriz a ser randomizada
matriz.original = mt,
# Número de trocas
n = 1000,
# Número de elementos que irão entrar na matriz taboo
n.taboo = 3)## $Matriz.Original
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0 1 0 1 0
## [2,] 1 0 1 0 1
## [3,] 0 1 0 0 0
## [4,] 1 0 0 0 0
## [5,] 0 1 0 0 0
##
## $Matriz.Randomizada
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0 0 0 1 0
## [2,] 0 0 0 1 1
## [3,] 0 0 0 0 1
## [4,] 1 1 0 0 0
## [5,] 0 1 1 0 0
##
## $Matriz.Taboo.n.5
## [,1] [,2] [,3] [,4]
## [1,] 4 2 2 1
## [2,] 3 1 5 3
## [3,] 4 1 5 4
##
## $N.de.vezes.elementos.iguais
## [1] 1002
##
## $Taboo.vezes.entrada
## [1] 12718Agora a implementação das funções “chave” em C++ utilizando Rcpp e RcppArmadillo.
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <RcppArmadilloExtensions/sample.h>
using namespace Rcpp;
// Matriz taboo
// [[Rcpp::export]]
int tabooC(IntegerVector e1,
IntegerVector e2,
IntegerVector tab1,
IntegerVector tab2,
IntegerVector tab3,
IntegerVector tab4,
int ntaboo)
{
int flag = -1;
int i = 0;
for(i; i <ntaboo; i++)
{
//Rcout << tab[i] << "\n";
if(IntegerMatrix::is_na(tab1[i])==1)
{
break;
}
if(((e1[0]==tab1[i])&&(e1[1]==tab2[i]))||((e1[0]==tab3[i])&&(e1[1]==tab4[i]))||
((e2[0]==tab1[i])&&(e2[1]==tab2[i]))||((e2[0]==tab3[i])&&(e2[1]==tab4[i])))
{
flag = 0;
break;
}
else
{
flag = 1;
}
}
return flag;
}
// Sorteador
// [[Rcpp::export]]
IntegerVector sorteadorC(IntegerVector In,
IntegerVector Jn)
{
int i;
int j;
//ij eh a posicao sorteada para a iteração
i = Rcpp::RcppArmadillo::sample(In,1,1)[0];
j = Rcpp::RcppArmadillo::sample(Jn,1,1)[0];
while(i<=j)
{
i = Rcpp::RcppArmadillo::sample(In,1,1)[0];
j = Rcpp::RcppArmadillo::sample(Jn,1,1)[0];
}
IntegerVector aux;
aux = {i,j};
return aux;
} Agora vamos fazer uma comparação entre as duas funções utilizando o pacote `microbenchmark´.
Primeiro vamos medir o tempo de uma única amostra.
# Funcoes em C
tabooCpre <- function(e1,e2,mt,n){
t1 <- mt[c(1:n),1]
t2 <- mt[c(1:n),2]
t3 <- mt[c(1:n),3]
t4 <- mt[c(1:n),4]
flag <- tabooC(e1,e2,t1,t2,t3,t4,n)
return(flag)
}
smrC = function(matriz.original, n = 1000, n.taboo = 5)
{
tab = matrix(nrow = n.taboo ,ncol = 4)
m = matriz.original
m.or = matriz.original
aux1 = aux2 = e1 = e2 = flag = NULL # e1 e e2 sao as posições dos elementos de troca
o = t = ta = 0 #Contador de iterações do while
In = c(1:length(m.or[,1])) #range para o nº de linhas
Jn = c(1:length(m.or[1,])) #range para o nº de colunas
for(i in 1:n)
{
e1 = sorteadorC(In,Jn)
e2 = sorteadorC(In,Jn)
while(m.or[e1[1],e1[2]]==m.or[e2[1],e2[2]])
{
e1 = sorteadorC(In,Jn)
e2 = sorteadorC(In,Jn)
o = o + 1
}
if(t<n.taboo)
{t = t +1}
else
{t = 0}
tab[t,1] = e1[1]
tab[t,2] = e1[2]
tab[t,3] = e2[1]
tab[t,4] = e2[2]
if(n.taboo>0){
flag = tabooCpre(e1,e2,tab,n.taboo) #Sinalizador se o indice sorteado pertence a matriz taboo
while(flag < 1||m.or[e1[1],e1[2]]==m.or[e2[1],e2[2]])
{
e1 = sorteadorC(In,Jn)
e2 = sorteadorC(In,Jn)
flag = tabooCpre(e1,e2,tab,n.taboo)
ta = ta + 1
}
}
tab[t,1] = e1[1]
tab[t,2] = e1[2]
tab[t,3] = e2[1]
tab[t,4] = e2[2]
aux1 = m.or[e1[1],e1[2]]
aux2 = m.or[e1[2],e1[1]]
m.or[e1[1],e1[2]] = m.or[e2[1],e2[2]]
m.or[e1[2],e1[1]] = m.or[e2[2],e2[1]]
m.or[e2[1],e2[2]] = aux1
m.or[e2[2],e2[1]] = aux2
}
resu = list(m,m.or,tab,o,ta)
names(resu) = c("Matriz.Original","Matriz.Randomizada","Matriz.Taboo.n.5",
"N.de.vezes.elementos.iguais","Taboo.vezes.entrada")
return(resu)
}
library(microbenchmark)
desempenho.mt <- microbenchmark(list = list(smrR = smr(matriz.original = mt,
n.taboo = 3),
smrC = smrC(matriz.original = mt,
n.taboo = 3)),
times = 1)
desempenho.mt## Unit: nanoseconds
## expr min lq mean median uq max neval
## smrR 2161 2161 2161 2161 2161 2161 1
## smrC 675 675 675 675 675 675 1Note que a função do R mostra-se bem inferior a que foi otimizada por C++.
Vamos agora fazer uma amostra de tamanho 10000 e comparar o desempenho das duas funções.
desempenho.mt <- microbenchmark(list = list(smrR = smr(matriz.original = mt,
n.taboo = 3),
smrC = smrC(matriz.original = mt,
n.taboo = 3)),
times = 10000)
desempenho.mt## Unit: nanoseconds
## expr min lq mean median uq max neval
## smrR 9 11 18.3024 13 14 8710 10000
## smrC 9 11 18.6775 13 14 6990 10000autoplot.microbenchmark(desempenho.mt)## Loading required namespace: ggplot2
wilcox.test(desempenho.mt$time~desempenho.mt$expr)##
## Wilcoxon rank sum test with continuity correction
##
## data: desempenho.mt$time by desempenho.mt$expr
## W = 49897000, p-value = 0.7983
## alternative hypothesis: true location shift is not equal to 0Veja que as duas se comportam de forma muito parecidade e que o teste de Wilcox(Distribuição assimétrica) não rejeitou \(H_0\), sendo assim, elas mostram-se com um desempenho igual (estatísticamente falando).
As performances foram muito parecidas, porém, como o processo de cada uma é totalmente aleátorio, essa ainda não é uma comparação tão confiável. Mesmo assim, a função em C++ mostra-se bem superior a em R quando utilizada em outras máquinas.
Sendo assim, a implementação de parte da função em C++ valeu a pena.