RcppEigen: Operações Matriciais e Rcpp - Parte 2

Uma das vantagens de se trabalhar com o RcppEigen e Rcpp, além da velocidade, é o fato de que não é necessário implementar todas as funções de interesse em C++. Podemos importar funções já existentes do R para que essas sejam executadas no ambiente Rcpp.

Por exemplo, um problema comum na formação de portfólios é o fato das matrizes de correlações estimadas entre os retornos dos ativos não ser Positivas Definidas. Nesses casos, problemas teóricos e numéricos podem surgir impedindo a obtenção de uma solução ótima para o problema de alocação de portfólios. Higham, N. (2002) propôs uma abordagem numérica para se encontrar a matriz Positiva Definida mais próxima da matriz de ineteresse.

Essa função está disponível no pacote Matrix e é denominada nearPD. Podemos importá-la no ambiente Rcpp da seguinte forma:

// [[Rcpp::export]]
Eigen::MatrixXd nearPDefinite(Eigen::MatrixXd X, int maxit=1e+6, double eigtol = 1e-06, double conv_tol = 1e-07, double posd_tol = 1e-08, bool keepDiagonal=false){
  //Cria o objeto com o Environment.
  Rcpp::Environment G = Rcpp::Environment::global_env();
  //Objeto com as informações do pacote Matrix
  Rcpp::Environment Matrix("package:Matrix");
  //Importa a função de ineteresse
  Function f = Matrix["nearPD"];
  //Obtêm o resultado da função nearPD
  Rcpp::List res = f(X, false, keepDiagonal, true, false, true, false, false,  eigtol, conv_tol, posd_tol, 100, "I", false);
  //Pega o primeiro elemento da lista e transforma para matrix
  Rcpp::NumericMatrix mat = internal::convert_using_rfunction(wrap(res[0]), "as.matrix");
  //Converte de Rcpp::NumericMatrix para Eigen::MatrixXd
  Eigen::MatrixXd Xpd=convertMatrix(mat);
  //Retorna o resultado
  return(Xpd);
}

A função anterior, após compilada, pode ser utilizada no ambiente R de maneira direta pelo nome nearPDefinite. Entretanto, é necessário que a biblioteca Matrix tenha sido habilitada e esteja disponível no environment do R para que seja utilizada. Essa estratégia pode ser para outras funções de outros pacotes também.

RcppEigen: Operações Matriciais e Rcpp - Parte 1

Uma das grandes críticas que os usuários do R fazem é que por vezes algumas operações são muito demoradas no R. Em especial, a execução de loops.

Há, no entanto alguns truques que podem ser utilizados para se aumentar a velocidade dos códigos produzidos em R, quais sejam: usar operações algébricas vetoriais ao invés de loops, utilizar os membros da família apply, utilizar processamento em paralelo , etc.

Entretanto, algumas vezes é necessário programar funções mais complexas e que necessitam do uso de loops clássicos. Quando esse for o caso, uma boa solução é utilizar o pacote RcppEigen que é uma versão em R da biblioteca em C++ Eigen.

Nesse grupo de posts acerca do RcppEigen vamos mostrar como esse pacote pode ser usado e também como novos pacotes para o R podem ser criados. O primeiro passo é construir um projeto de pacote no RStudio:

Definido a criação de um novo pacote o próximo passo é determinar:

Podemos criar um projeto com base em um diretório já existente ou em um novo diretório. Escolhendo um novo diretório temos as seguintes opções:

Escolhemos nesse caso a opção R Package. Em seguida, definimos o local e nome do pacote que desejamos criar:

Automaticamente o R cria alguns arquivos e pastas. Para usar o RcppEigen precisamos de alguns arquivos específicos. Esses arquivos são obtidos executando uma função padrão do pacote:

#Limpa o Workspace
#Habilita o pacote
library(RcppEigen)
#Define o local temporário para criação dos arquivos
setwd("C:\\LocalTemorario")
#Cria os arquivos necessários
RcppEigen.package.skeleton("ArquivosRcppEigen")

Na pasta definida pelo comando setwd um conjunto de arquivos é criado.

Deve-se copiar os arquivos DESCRIPTION e NAMESPACE e também a pasta src (todas elas dentro da pasta "ArquivosRcppEigen") para a pasta do pacote criado no RStudio (nesse exemplo denominada "NomeDoPacote"). Após essas etapas o pacote estará pronto para receber as funções escritas em RcppEigen, as quais veremos nas próximas partes desse tutorial.

taskscheduleR: Agendando tarefas no R.

A automação de tarefas é importante quando deseja-se executar um mesmo script de maneira regular. Nesse sentido o pacote taskscheduleR pode apoiar no agendamento de scripts feitos em R.

Como exemplo vamos supor que estamos interessados em decidir se devemos comprar ou vender uma ação no dia seguinte baseado no preço de fechamento dos dias anteriores usando um modelo de Máquinas de Suporte Vetorial.

#Limpa o Workspace
rm(list=ls())
#Chama a biblioteca
library(devtools)
library(quantmod)
library(kernlab)
#Instala o pacote taskscheduleR
devtools::install_github("jwijffels/taskscheduleR")
#Chama a biblioteca taskscheduleR
library(taskscheduleR)
#Especifica as datas de interesse
startDate = as.Date("2013-01-01") 
endDate = as.Date("2015-12-31")
#Obtêm os dados do ativo PETR4 e PETR3
getSymbols("PETR4.SA", src="yahoo",from=startDate,to=endDate)
#Faz o Candle Chart
candleChart(PETR4.SA,theme='white', type='candles') 
reChart(major.ticks='months',subset='first 16 weeks') 

O Candle Chart é apresentado a seguir:

Em seguida vamos treinar uma máquina para tentar apontar a direção do resultado futuro, isto é, se a ação irá apresentar um resultado positivo ou negativo:

#Retorno do Preço de Fechamento
ret<-na.omit(ClCl(PETR4.SA))
#Direção do mercado
y<-as.numeric(ifelse(ret<0,0,1))
#Cria a base de dados
df<-data.frame("Y"=c(y[-1],NA),"X"=ret)
#Renomeia as colunas
colnames(df)<-c("Y","X")

Nesse bloco a ideia é tentar prever a direção $Y=sinal(r_{t})$ usando o retorno do período anterior $X=r_{t-1}$. Supondo que o processo de treinamento tenha sido realizado e que tenhamos descobertos $C=0.1$ e $\sigma=16$ para o modelo C-SVC e kernel Gaussiano, a máquina é construída fazendo-se:

#Constrói a máquina
svm <- ksvm(Y~X,data=df, type="C-svc",C=0.1,
            kpar=list(sigma=16),cross=3)
#Faz a previsão
pred<- kernlab::predict(svm,df)
#Compara com os valores reais
table(pred,df$Y)
#Salva a máquina treinada
save.image("Maquina.RData")

Supondo que essa acurácia seja a desejada, agora vamos programar para que o R todo dia baixe a cotação do dia anterior e faça a previsão de qual deveria ser nossa ação no dia atual (vender se a previsão é que o retorno seja positivo ou comprar se a previsão seja do retorno negativo). As opções do pacote taskscheduleR são: ‘ONCE’, ‘MONTHLY’, ‘WEEKLY’, ‘DAILY’, ‘HOURLY’, ‘MINUTE’, ‘ONLOGON’, ‘ONIDLE’. O primeiro passo é montar um script geral que será executado todo dia:

#Chama a biblioteca
library(quantmod)
library(kernlab)
#Lê a máquina criada
load("Maquina.RData")
#Obtêm a cotação atual
getSymbols("PETR4.SA", src="yahoo")
#Último retorno observado
ret<-tail(ClCl(PETR4.SA),1)
#Faz a previsao para o dia
previsao<-kernlab::predict(svm,ret)
#Monta o arquivo de saida
fileConn<-file("Decisao.txt")
writeLines(ifelse(previsao==1,"Vender","Comprar"), fileConn)
close(fileConn)

Esse script deve ser salvo em um arquivo de extensão .R (por exemplo, ScriptPrevisao.R) e, no nosso exemplo, ele deverá ser executado todo dia as 8:00 AM. Após a execução, ele criará um arquivo denominado Decisao.txt que contêm a decisão do dia: Comprar ou Vender. O último passo é programar o computador para executar esse script todo dia no horário determinado, para isso, basta executarmos uma única vez no R:

#Faz o agendamento
library(taskscheduleR)
myscript <- system.file("saida", "ScriptPrevisao.R", package = "taskscheduleR")
#O computador irá executar todo dia as 8:00AM
taskscheduler_create(taskname = "DecisaoPETR4", rscript = myscript, 
                     schedule = "DAILY", starttime = "08:00")

Caso queira deletar o agendamento para que ele não se repita mais, basta fazer:

#Deleta o agendamento
library(taskscheduleR)
taskscheduler_delete(taskname = "DecisaoPETR4")

Introdução ao RCpp - Parte 4.

Frequentemente é necessário inicializar um NumericVector com algum valor de constante fixo, nesse caso podemos usar a função fill(.):

#include 
using namespace Rcpp;
// Função para repetir o valor no vetor
// [[Rcpp::export]]
NumericVector inicializaVetor(double valor, int tamanho) {
  NumericVector x(tamanho);
  x.fill(valor);
  return x;
}

Após compilar o código, a invocação no R é quase imediata:

#Invoca a biblioteca RCpp
library(Rcpp)
#Define o endereço do arquivo Cpp
setwd("C:\\Blog\\Source")
#Compila o código em C++
sourceCpp('Teste.cpp')
inicializaVetor(2,5)

O pacote Rcpp tem muitas funcionalidades. Por meio do Rcpp podemos por exemplo, usar funções já implementadas no R, como:

// [[Rcpp::export]]
NumericVector calculaCurtose(NumericVector x) {
  Rcpp::Environment funcoesPacote("package:e1071");
  Rcpp::Function funcaoCurtose = funcoesPacote["kurtosis"];
  return funcaoCurtose (x);
}

Ou por exemplo, gerar números com a distribuição normal:

// [[Rcpp::export]]
NumericVector geraNormal(){
  Rcpp::Environment stats("package:stats");
  Rcpp::Function rnorm = stats["rnorm"];
  return rnorm(10, Rcpp::Named("sd", 100.0));
}

Comandos esses que após compilados são invocados na forma:

#Invoca a biblioteca RCpp
library(Rcpp)
#Define o endereço do arquivo Cpp
setwd("C:\\Blog\\Source")
#Compila o código em C++
sourceCpp('Teste.cpp')
#Executa as funções:
x<-rnorm(1919)
calculaCurtose(x)
geraNormal()

Entre os objetos mais flexíveis do Rcpp está o objeto GenericVector o qual funciona como uma lista permitindo a inserção de qualquer outro objeto em seus slots. Além do mais, é possível usar os seguintes métodos nesse objeto:

• push_back(x) - Insere o elemento x no final da lista.
• push_front(x) - Insere o elemento x no início da lista.
• insert(i, x) - Insere o elemento x na posição i.
• erase(i) - Apaga o elemento existente na posição i.

Todas essas funções são dinâmicas e não há a necessidade de criar uma nova lista para que os novos elementos sejam inseridos.

Introdução ao RCpp - Parte 3.

Aqui neste post continuaremos com a análise dos possíveis métodos desenvovidos usando Rcpp. O primeiro exemplo trata do uso de matrizes, especificamente, vamos criar uma função para fazer a transposição de uma matriz:

#include 
using namespace Rcpp;
// Função para transposição de matriz
// [[Rcpp::export]]
NumericMatrix transpose(NumericMatrix x) {
  //Obtêm o número de linhas e colunas:
  int nrow = x.nrow(), ncol = x.ncol();
  //Cria a matriz com resultado
  NumericMatrix resultado(ncol,nrow);
  //Faz a transposição:
  for(int i=0;i < nrow; i++){
    for(int j=0;j < ncol; j++){
      resultado(j,i)=x(i,j);
    }
  }
  return resultado;
}

Após compilar o código, a invocação no R é quase imediata:

#Invoca a biblioteca RCpp
library(Rcpp)
#Define o endereço do arquivo Cpp
setwd("C:\\Blog\\Source")
#Compila o código em C++
sourceCpp('Teste.cpp')
#Cria uma matriz de exemplo
mat<-matrix(c(2,8,9,0,0,0,4,3,2), nrow=3,ncol=3)
#Transpoem a matriz
transpose(mat)

Métodos mais sofisticados usando funções implementadas em R e invocadas em C++ também podem ser utilizadas:

//Usando funções do R:
// [[Rcpp::export]]
NumericVector ExemploKernel(NumericVector x, NumericVector y, Function f) {
  NumericVector res = f(x,y);
  return res;
}

Método esse que recebe como argumento dois vetores numéricos e uma função bivariada:

#Cria uma função no R
gaussian<-function(x,y){
  return(sum((x-y)^2))
}
#Gera dois vetores
x<-rnorm(1000)
y<-rnorm(1000)
#Invoca a função em Rcpp:
ExemploKernel(x,y,gaussian)

Introdução ao RCpp - Parte 2.

Dando continuidade ao post de Introdução ao RCpp - Parte 1. Aqui eu mostrarei como os objetos do Rcpp se relacionam com os objetos em R. Os principais objetos do Rcpp são:

1. NumericVector - Objeto representando um vetor numérico no R.
2. CharacterVector - Objeto representando um vetor caracter no R.
3. NumericMatrix - Objeto representando uma matriz numérica no R.
4. CharacterMatrix - Objeto representando uma matriz caracter no R.
5. LogicalVector - Objeto representando um vetor lógico no R.
6. IntegerVector - Objeto representando um vetor inteiro no R.
7. LogicalMatrix - Objeto representando uma matriz lógica no R.
8. IntegerMatrix - Objeto representando uma matriz inteira no R.
9. List - Objeto representando uma lista no R.

Aqui trabalharemos com exemplos de vetores e nos próximos posts o uso de matrizes e listas. Suponha que temos um vetor de retornos e um vetor com os nomes dos ativos representando esses vetores, algo como:

#Invoca a biblioteca RCpp
library(Rcpp)
#Define o endereço do arquivo Cpp
setwd("C:\\Blog\\Source")
#Define a semente para simulação dos dados
set.seed(19393)
#Gera um vetor numérico de retornos:
ret<-rnorm(1000, 0,0.1)
#Gera um vetor caracter de nomes:
nomes<-replicate(1000, paste(sample(LETTERS, 5, replace=TRUE), collapse=""))

Com base nesses dois inputs podemos pensar em algumas funções escritas em Rcpp que utilizem esses insumos. Por exemplo, considere que desejamos obter os 5 menores valores do vetor ret e os respectivos nomes dessas empresas.

//Exemplo de função para obter os k menores valores.
// [[Rcpp::export]]
NumericVector kMenores(int k, NumericVector ret, CharacterVector nomes) {
  //Cria um vetor numérico de tamanho k
  NumericVector vetMenoresRetornos(k);
  //Cria um vetor caracter de tamanho k
  CharacterVector vetMenoresNomes(k);
  //Faz o loop para achar os k menores
  for (int i = 0; i < ret.size(); i++){
    //Verifica qual menor é:
    int iMenor = qualMenor(ret[i],vetMenoresRetornos);
    if(iMenor>-1){
      vetMenoresRetornos[iMenor]=ret[i];
    }
    
  }
  return vetMenoresRetornos;
}

Onde a função qualMenor é criada no topo do mesmo arquivo em C++ da seguinte forma:

#include 
using namespace Rcpp;

//Função para dizer qual posição esse menor é:
int qualMenor(double valor,NumericVector vetMenoresRetornos){
  //Inicializa a variável
  int menor=-1;
  //Procura qual o menor:
  for(int i=0;i< vetMenoresRetornos.size();i++){
    if(valor<=vetMenoresRetornos[i])  return i;
  }
  return(menor);
}

Observe que as linhas 1 e 2 do código anterior são mandatórias e devem constar em todos os arquivo com extensão .cpp, ademais como o código anterior não contêm a linha // [[Rcpp::export]] isso significa que a função é interna e não pode ser invocada pelo R. Finalmente, a invocação no R é feito da seguinte forma:

#Invoca a biblioteca RCpp
library(Rcpp)
#Define o endereço do arquivo Cpp
setwd("C:\\Blog\\Source")
#Compila o código em C++
sourceCpp('Teste.cpp')
#Inova a função compilada
teste<-kMenores(5,ret,nomes)
#Confere com a função do R
head(sort(teste),5)

Repare que em nenhum momento utilizamos o objeto nomes. Para imprimir o nome das empresas com menor retorno no console podemos alterar a função kMenores da seguinte forma:

//Exemplo de função para obter os k menores valores.
// [[Rcpp::export]]
NumericVector kMenores(int k, NumericVector ret, CharacterVector nomes) {
  //Cria um vetor numérico de tamanho k
  NumericVector vetMenoresRetornos(k);
  //Cria um vetor caracter de tamanho k
  CharacterVector vetMenoresNomes(k);
  //Faz o loop para achar os k menores
  for (int i = 0; i < ret.size(); i++){
    //Verifica qual menor é:
    int iMenor = qualMenor(ret[i],vetMenoresRetornos);
    if(iMenor>-1){
      //Guarda os retornos mínimos
      vetMenoresRetornos[iMenor]=ret[i];
      //Guarda os nomes das firmas
      vetMenoresNomes[iMenor]=nomes[i];
    }
  }
  //Antes de retornar os resultados imprimir no console o nome
  for(int j=0;j < k;j++){
    std::cout << "Firma: " << vetMenoresNomes[j] << " Retorno: " << vetMenoresRetornos[j] << std::endl;
  }
  return vetMenoresRetornos;
}

A invocação procede de maneira equivalente a apresentada anteriormente, a diferença é que agora no console do R o nome das firmas e seus retornos simulados são apresentados. Alguns links interessantes para mais detalhes sobre o RCpp:

1. http://adv-r.had.co.nz/Rcpp.html
2. http://gallery.rcpp.org/articles/r-function-from-c++/
3. http://pages.stat.wisc.edu/~gvludwig/Rcpp/

Introdução ao RCpp - Parte 1.

O R possui uma limitação devido a demora de alguns procedimentos realizados nele. O software foi criado para ser ótimo do ponto de vista de processamento computacional quando o usuário utiliza as informações de maneira vetorizada. Quando o processo é realizado de elemento a elemento alguma demora no processamento surge.

Nesse sentido, aqui daremos uma primeira introdução quanto ao uso do pacote RCpp (http://www.rcpp.org/), o qual permite otimizar alguns procedimentos em R usando linguagem C++.

O primeiro passo para desenvolver suas próprias funções em RCpp é ter instalado:

1. R
2. RTools - É preciso que ele esteja instalado na raiz (C:\) e com todas as opções habilitadas.
3. RStudio

Em seguida, abrimos o RStudio e escolhemos a opção File → New File → C++ File:

O RStudio fornecerá um exemplo de como construir funções usando RCpp:

Onde o primeiro bloco de código, da linha 1 até a 24 é código escrito em C++, o que está contido na linha 26 é o código em R usando para a invocação da função em C++. Para compilar a função é necessário clicar no botão Source destacado na figura anterior. Ao clicar nesse botão a função é compilada e fica disponível para ser executada.

É comum no entanto ter dois arquivos: um arquivo com as funções em C++ e outro arquivo com as funções em R, por isso considere o seguinte código para o arquivo com extensão .cpp:

#include 
using namespace Rcpp;
// This is a simple example of exporting a C++ function to R. 
// [[Rcpp::export]]
NumericVector timesTwo(NumericVector x) {
  return x * 2;
}

E no arquivo em R (com extensão .R):

#Invoca a biblioteca RCpp
library(Rcpp)
#Define o endereço do arquivo Cpp
setwd("C:\\Blog\\Source")
#Compila o código em C++
sourceCpp('Teste.cpp')
#Inova a função compilada
timesTwo(42)

Note que nesse exemplo, Teste.cpp é o nome do arquivo com o código em RCpp armazenado na pasta C:\Blog\Source. No próximo post mostrarei como fazer funções mais complexas envolvendo os pacotes do R e seus objetos com o Rcpp.