Articles by Rcpp Gallery

Performance of the divide-and-conquer SVD algorithm

December 9, 2013 | Rcpp Gallery

The ubiquitous LAPACK library provides several implementations for the singular-value decomposition (SVD). We will illustrate possible speed gains from using the divide-and-conquer method by comparing it to the base case.
<span>#include <RcppArmadillo.h></span>

<span>// [[Rcpp::depends(RcppArmadillo)]]</span>

<span>// [[Rcpp::export]]</span>
<span>arma</span><span>::</span><span>vec</span> <span>baseSVD</span><span>(</span><span>const</span> <span>arma</span><span>::</span><span>mat</span> <span>&</span> <span>X</span><span>)</span> <span>{</span>
    <span>arma</span><span>::</span><span>mat</span> <span>U</span><span>,</span> <span>V</span><span>;</span>
    <span>arma</span><span>::</span><span>vec</span> <span>S</span><span>;</span>
    <span>arma</span><span>::</span><span>svd</span><span>(</span><span>U</span><span>,</span> <span>S</span><span>,</span> <span>V</span><span>,</span> <span>X</span><span>,</span> <span>"standard"</span><span>);</span>
    <span>return</span> <span>S</span><span>;</span>
<span>}</span>

<span>// [[Rcpp::export]]</span>
<span>arma</span><span>::</span><span>vec</span> <span>dcSVD</span><span>(</span><span>const</span> <span>arma</span><span>::</span><span>mat</span> <span>&</span> <span>X</span><span>)</span> <span>{</span>
    <span>arma</span><span>::</span><span>mat</span> <span>U</span><span>,</span> <span>V</span><span>;</span>
    <span>arma</span><span>::</span><span>vec</span> <span>S</span><span>;</span>
    <span>arma</span><span>::</span><span>svd</span><span>(</span><span>U</span><span>,</span> <span>S</span><span>,</span> <span>V</span><span>,</span> <span>X</span><span>,</span> <span>"dc"</span><span>);</span>
    <span>return</span> <span>S</span><span>;</span>
<span>}</span>
Having the two implementations, which differ only in the method argument (added recently in Armadillo 3.930), we are ready to ...
[Read more...]

Converting C code to C++ code: An example from plyr

December 2, 2013 | Rcpp Gallery

The plyr package uses a couple of small C functions to optimise a number of particularly bad bottlenecks. Recently, two functions were converted to C++. This was mostly stimulated by a segmentation fault caused by some large inputs to the split_indices() function: rather than figuring out exactly what was ...
[Read more...]

Munkres’ Assignment Algorithm with RcppArmadillo

September 24, 2013 | Rcpp Gallery

Munkres’ Assignment Algorithm (Munkres (1957), also known as hungarian algorithm) is a well known algorithm in Operations Research solving the problem to optimally assign N jobs to N workers. I needed to solve the Minimal Assignment Problem for a relabeling algorithm in MCMC sampling for finite mixture distributions, where I use ...
[Read more...]

Creating as and wrap for sparse matrices

August 5, 2013 | Rcpp Gallery

An earlier article discussed sparse matrix conversion but stopped short of showing how to create custom as() and wrap() methods or functions. This post starts to close this gap. We will again look at sparse matrices from the Matrix package for R, as well as the SpMat class from Armadillo....
[Read more...]

Gibbs Sampler in C++

August 3, 2013 | Rcpp Gallery

Markov Chain Monte Carlo (MCMC) is a popular simulation method. As it is somewhat demanding, it is also frequently used to benchmark different implementations or algorithms. One particular algorithm has been compared a number of times, starting with an article by Darren Wilkinson, and Darren’s follow–up article which ...
[Read more...]

Faster Multivariate Normal densities with RcppArmadillo and OpenMP

July 13, 2013 | Rcpp Gallery

The Multivariate Normal density function is used frequently in a number of problems. Especially for MCMC problems, fast evaluation is important. Multivariate Normal Likelihoods, Priors and mixtures of Multivariate Normals require numerous evaluations, thus speed of computation is vital. We show a twofold increase in speed by using RcppArmadillo, and ...
[Read more...]

Sobol Sensitivity Analysis

June 10, 2013 | Rcpp Gallery

Sensitivity analysis is the task of evaluating the sensitivity of a model output Y to input variables (X1,…,Xp). Quite often, it is assumed that this output is related to the input through a known function f :Y= f(X1,…,Xp). Sobol indices are generalizing the coefficient of the coefficient ...
[Read more...]

Dynamic Wrapping and Recursion with Rcpp

April 8, 2013 | Rcpp Gallery

We can leverage small parts of the R’s C API in order to infer the type of objects directly at the run-time of a function call, and use this information to dynamically wrap objects as needed. We’ll also present an example of recursing through a list. To get ...
[Read more...]

Using bigmemory with Rcpp

March 14, 2013 | Rcpp Gallery

The bigmemory package allows users to create matrices that are stored on disk, rather than in RAM. When an element is needed, it is read from the disk and cached in RAM. These objects can be much larger than native R matrices. Objects stored as such larger-than-RAM matrices are defined ...
[Read more...]

Generating a multivariate gaussian distribution using RcppArmadillo

March 12, 2013 | Rcpp Gallery

There are many ways to simulate a multivariate gaussian distribution assuming that you can simulate from independent univariate normal distributions. One of the most popular method is based on the Cholesky decomposition. Let’s see how Rcpp and Armadillo perform on this task.
<span>#include <RcppArmadillo.h></span>
<span>// [[Rcpp::depends(RcppArmadillo)]]</span>

<span>using</span> <span>namespace</span> <span>Rcpp</span><span>;</span>

<span>// [[Rcpp::export]]</span>
<span>arma</span><span>::</span><span>mat</span> <span>mvrnormArma</span><span>(</span><span>int</span> <span>n</span><span>,</span> <span>arma</span><span>::</span><span>vec</span> <span>mu</span><span>,</span> <span>arma</span><span>::</span><span>mat</span> <span>sigma</span><span>)</span> <span>{</span>
   <span>int</span> <span>ncols</span> <span>=</span> <span>sigma</span><span>.</span><span>n_cols</span><span>;</span>
   <span>arma</span><span>::</span><span>mat</span> <span>Y</span> <span>=</span> <span>arma</span><span>::</span><span>randn</span><span>(</span><span>n</span><span>,</span> <span>ncols</span><span>);</span>
   <span>return</span> <span>arma</span><span>::</span><span>repmat</span><span>(</span><span>mu</span><span>,</span> <span>1</span><span>,</span> <span>n</span><span>).</span><span>t</span><span>()</span> <span>+</span> <span>Y</span> <span>*</span> <span>arma</span><span>::</span><span>chol</span><span>(</span><span>sigma</span><span>);</span>
<span>}</span>
The easiest way to perform a ...
[Read more...]

Using Rcpp with Boost.Regex for regular expression

March 1, 2013 | Rcpp Gallery

Gabor asked about Rcpp use with regular expression libraries. This post shows a very simple example, based onone of the Boost.RegEx examples. We need to set linker options. This can be as simple as
Sys.setenv<span>(</span><span>"PKG_LIBS"</span><span>=</span><span>"-lboost_regex"</span><span>)</span>
With that, the following example can be built:
<span>// cf www.boost.org/doc/libs/1_53_0/libs/regex/example/snippets/credit_card_example.cpp</span>

<span>#include <Rcpp.h></span>

<span>#include <string></span>
<span>#include <boost/regex.hpp></span>

<span>bool</span> <span>validate_card_format</span><span>(</span><span>const</span> <span>std</span><span>::</span><span>string</span><span>&</span> <span>s</span><span>)</span> <span>{</span>
   <span>static</span> <span>const</span> <span>boost</span><span>::</span><span>regex</span> <span>e</span><span>(</span><span>"(</span><span>\\</span><span>d{4}[- ]){3}</span><span>\\</span><span>d{4}"</span><span>);</span>
   <span>return</span> <span>boost</span><span>::</span><span>regex_match</span><span>(</span><span>s</span><span>,</span> <span>e</span><span>);</span>
<span>}</span>

<span>const</span> <span>boost</span><span>::</span><span>regex</span> <span>e</span><span>(</span><span>"</span><span>\\</span><span>A(</span><span>\\</span><span>d{3,4})[- ]?(</span><span>\\</span><span>d{4})[- ]?(</span><span>\\</span><span>d{4})[- ]?(</span><span>\\</span><span>d{4})</span><span>\\</span><span>z"</span><span>);</span>
<span>const</span> <span>std</span><span>::</span><span>string</span> <span>machine_format</span><span>(</span><span>"</span><span>\\</span><span>1</span><span>\\</span><span>2</span><span>\\</span><span>3</span><span>\\</span><span>4"</span><span>);</span>
<span>const</span> <span>std</span><span>::</span><span>string</span> <span>human_format</span><span>(</span><span>"</span><span>\\</span><span>1-</span><span>\\</span><span>2-</span><span>\\</span><span>3-</span><span>\\</span><span>4"</span><span>);</span>

<span>std</span><span>::</span><span>string</span> <span>machine_readable_card_number</span><span>(</span><span>const</span> <span>std</span><span>::</span><span>string</span><span>&</span> <span>s</span><span>)</span> <span>{</span>
   <span>return</span> <span>boost</span><span>::</span><span>regex_replace</span><span>(</span><span>s</span><span>,</span> <span>e</span><span>,</span> <span>machine_format</span><span>,</span> <span>boost</span><span>::</span><span>match_default</span> <span>|</span> <span>boost</span><span>::</span><span>format_sed</span><span>);</span>
<span>}</span>

<span>std</span><span>::</span><span>string</span> <span>human_readable_card_number</span><span>(</span><span>const</span> <span>std</span><span>::</span><span>string</span><span>&</span> <span>s</span><span>)</span> <span>{</span>
   <span>return</span> <span>boost</span><span>::</span><span>regex_replace</span><span>(</span><span>s</span><span>,</span> <span>e</span><span>,</span> <span>human_format</span><span>,</span> <span>boost</span><span>::</span><span>match_default</span> <span>|</span> <span>boost</span><span>::</span><span>format_sed</span><span>);</span>
<span>}</span>

<span>// [[Rcpp::export]]</span>
<span>Rcpp</span><span>::</span><span>DataFrame</span> <span>regexDemo</span><span>(</span><span>std</span><span>::</span><span>vector</span><span><</span><span>std</span><span>::</span><span>string</span><span>></span> <span>s</span><span>)</span> <span>{</span>
    <span>int</span> <span>n</span> <span>=</span> <span>s</span><span>.</span><span>size</span><span>();</span>
    
    <span>std</span><span>::</span><span>vector</span><span><</span><span>bool</span><span>></span> <span>valid</span><span>(</span><span>n</span><span>);</span>
    <span>std</span><span>::</span><span>vector</span><span><</span><span>std</span><span>::</span><span>string</span><span>></span> <span>machine</span><span>(</span><span>n</span><span>);</span>
    <span>std</span><span>::</span><span>vector</span><span><</span><span>std</span><span>::</span><span>string</span><span>></span> <span>human</span><span>(</span><span>n</span><span>);</span>
    
    <span>for</span> <span>(</span><span>int</span> <span>i</span><span>=</span><span>0</span><span>;</span> <span>i</span><span><</span><span>n</span><span>;</span> <span>i</span><span>++</span><span>)</span> <span>{</span>
        <span>valid</span><span>[</span><span>i</span><span>]</span>  <span>=</span> <span>validate_card_format</span><span>(</span><span>s</span><span>[</span><span>i</span><span>]);</span>
        <span>machine</span><span>[</span><span>i</span><span>]</span> <span>=</span> <span>machine_readable_card_number</span><span>(</span><span>s</span><span>[</span><span>i</span><span>]);</span>
        <span>human</span><span>[</span><span>i</span><span>]</span> <span>=</span> <span>human_readable_card_number</span><span>(</span><span>s</span><span>[</span><span>i</span><span>]);</span>
    <span>}</span>
    <span>return</span> <span>Rcpp</span><span>::</span><span>DataFrame</span><span>::</span><span>create</span><span>(</span><span>Rcpp</span><span>::</span><span>Named</span><span>(</span><span>"input"</span><span>)</span> <span>=</span> <span>s</span><span>,</span>
                                   <span>Rcpp</span><span>::</span><span>Named</span><span>(</span><span>"valid"</span><span>)</span> <span>=</span> <span>valid</span><span>,</span>
                                   <span>Rcpp</span><span>::</span><span>Named</span><span>(</span><span>"machine"</span><span>)</span> <span>=</span> <span>machine</span><span>,</span>
                                   <span>Rcpp</span><span>::</span><span>Named</span><span>(</span><span>"human"</span><span>)</span> <span>=</span> <span>human</span><span>);</span>
<span>}</span>
We can test the function ...
[Read more...]

Fast factor generation with Rcpp

February 27, 2013 | Rcpp Gallery

Recall that factors are really just integer vectors with ‘levels’, i.e., character labels that get mapped to each integer in the vector. How can we take an arbitrary character, integer, numeric, or logical vector and coerce it to a factor with Rcpp? It’s actually quite easy with Rcpp ...
[Read more...]

Sorting Numeric Vectors in C++ and R

January 31, 2013 | Rcpp Gallery

Consider the problem to sort all elements of the given vector in ascending order. We can simply use the function std::sort from the C++ STL.
<span>#include <Rcpp.h></span>
<span>using</span> <span>namespace</span> <span>Rcpp</span><span>;</span>

<span>// [[Rcpp::export]]</span>
<span>NumericVector</span> <span>stl_sort</span><span>(</span><span>NumericVector</span> <span>x</span><span>)</span> <span>{</span>
   <span>NumericVector</span> <span>y</span> <span>=</span> <span>clone</span><span>(</span><span>x</span><span>);</span>
   <span>std</span><span>::</span><span>sort</span><span>(</span><span>y</span><span>.</span><span>begin</span><span>(),</span> <span>y</span><span>.</span><span>end</span><span>());</span>
   <span>return</span> <span>y</span><span>;</span>
<span>}</span>
library<span>(</span>rbenchmark<span>)</span>
set.seed<span>(</span><span>123</span><span>)</span>
z <span><-</span> rnorm<span>(</span><span>100000</span><span>)</span>
x <span><-</span> rnorm<span>(</span><span>100</span><span>)</span>

<span># check that stl_sort is the same as sort</span>
stopifnot<span>(</span>all.equal<span>(</span>stl_sort<span>(</span>x<span>),</span> sort<span>(</span>x<span>)))</span>

<span># benchmark stl_sort and sort</span>
benchmark<span>(</span>stl_sort<span>(</span>z<span>),</span> sort<span>(</span>z<span>),</span> order<span>=</span><span>"relative"</span><span>)[,</span><span>1</span><span>:</span><span>4</span><span>]</span>
test replications elapsed relative 1 stl_sort(z) 100 0.632 1.000 2 sort(z) 100 1.164 1.842 Consider the problem of sorting the first n elements of a given vector. ...
[Read more...]

Using Boost via the new BH package

January 31, 2013 | Rcpp Gallery

Earlier today the new BH package arrived on CRAN. Over the years, Jay Emerson, Michael Kane and I had numerous discussions about a basic Boost infrastructure package providing Boost headers for other CRAN packages. JJ and Romain chipped in as well, and Jay finally took the lead by first creating ...
[Read more...]

Using Boost’s foreach macro

January 30, 2013 | Rcpp Gallery

Boost provides a macro, BOOST_FOREACH, that allows us to easily iterate over elements in a container, similar to what we might do in R with sapply. In particular, it frees us from having to deal with iterators as we do with std::for_each and std::transform. The macro ...
[Read more...]
1 2 3 4 5 6

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)