RcppArmadillo Test Notes – A Step-by-Step Tutorial to Create Your Own R Package Based on RcppArmadillo

Test environment:
R version 3.2.2
Rtools 33
Windows 7 64-bit

1. Install the packages “Rcpp”, “RcppArmadillo” and “inline” for R.

2. Download and install Rtools at https://cran.r-project.org/bin/windows/Rtools/, which is needed for building packages for R under Microsoft Windows. When installing, check the option “Edit the System PATH”.

3. Simple test code:

#See https://github.com/petewerner/misc/wiki/RcppArmadillo-cheatsheet for more examples.

rcpp_inc <- '
using namespace Rcpp;
using namespace arma;

m1 <- matrix(1:16, nr=4)
m2 <- matrix(17:32, nr=4)
v1 <- 1:10
v2 <- 11:20

src <- '
mat m1 = as<mat>(m1in);
mat m2 = as<mat>(m2in);
mat out = join_cols(m1, m2);

fn <- cxxfunction(signature(m1in="numeric", m2in="numeric"), src, plugin='RcppArmadillo', rcpp_inc)
res <- fn(m1, m2)
test <- rbind(m1, m2)
all.equal(test, res)

To write an R package based on RcppArmadillo:

4. Use the following R code:


You will get a folder “C:\Users\hp\Documents\testRcppArPackage” (here “testRcppArPackage” is your package name, you can change it to something more meaningful) in your R working directory, which contains all the basic files you need to build an R package.

5. Install the example package:
In windows CMD, “cd” to your R working directory, which contains the folder “testRcppArPackage”, and run the following commands

"D:\Program Files\R\R-3.2.2\bin\x64\Rcmd" check testRcppArPackage
"D:\Program Files\R\R-3.2.2\bin\x64\Rcmd" INSTALL testRcppArPackage

where “D:\Program Files\R\R-3.2.2\bin\x64\Rcmd” is the path of your R installation.
If you get the error that “Error: unexpected symbol in ’’‘~~simple examples’’”, open the file “C:\Users\hp\Documents\testRcppArPackage\man\testRcppArPackage-package.Rd” with your text editor, find the block
~~ simple examples of the most important functions ~~
and delete the content line inside the block and save, so it becomes
}”. Later on you may add example code in that block.
See http://howtomakeanrpackage.pbworks.com/f/How_To_Make_An_R_Package-v1.14-01-11-10.pdf Page 32 for more information.

6. Add your own functions in your package:
Just add your new functions in the file “C:\Users\hp\Documents\testRcppArPackage\src\rcpparma_hello_world.cpp”.
After that, you need to re-generate the “RcppExports.cpp” and “RcppExports.R” files in order to expose your own functions to R (so that you can call your defined functions in R after you load your package). To do this, use the following R code:


See https://cran.r-project.org/web/packages/Rcpp/vignettes/Rcpp-package.pdf for details.

7. Finally, re-build and re-install your package using the same commands in step 5. If you encounter any error during the building process, you can check the file “C:\Users\hp\Documents\testRcppArPackage.Rcheck\00install.out” to find the details of the errors. You can test your package by loading it using


Then you can call your function in R by the same function name defined in the C/C++ source code.

8. In addition, use the command

"D:\Program Files\R\R-3.2.2\bin\x64\Rcmd" build testRcppArPackage

to create a packed package file “testRcppArPackage_1.0.tar.gz” for distributing. You can install your package from this “.tar.gz” file locally using the following R command:

install.packages("testRcppArPackage_1.0.tar.gz", repos = NULL)

However, the above package contains source code only and requires users have compilers (Rtools) to install. To build a binary package that can install directly, use the command

"D:\Program Files\R\R-3.2.2\bin\x64\Rcmd" --build testRcppArPackage_1.0.tar.gz

And to build the binary package for both (Windows) 32-bit and 62-bit,

"D:\Program Files\R\R-3.2.2\bin\x64\Rcmd" --build --compile-both testRcppArPackage_1.0.tar.gz

You’ll get a new binary package file “testRcppArPackage_1.0.zip”, which can be installed directly (for Windows) without compiling. To install,

install.packages("testRcppArPackage_1.0.zip", repos = NULL)

Seamless R and C++ Integration with Rcpp (Use R!)

Posted in programming | Tagged , , | 5 Comments

.m Function to .oct File – Syntax Translation Table

One of the things Octave is superior to Matlab is Octave’s C++ code interface – Oct-Files. With Oct-Files, you can use a similar syntax in C++ as in the m language, which is much more convenient than the Matlab’s .mex extension. When your .m file is slow than you can bear, you’d better try to port it to .oct for a performance boost. However, there is not enough documents about Oct-Files. Here I just want to share a simple table for translating .m function to C++ for dynamic link .oct file.
The following note collects some commonly used syntax based on my own experience.

Function Define:

function out = function_name(args)
%get input
%return value
//To call the function: out = function_name(in);
DEFUN_DLD (function_name, args, nargout, "function description")
//get input length
int nargin = args.length();
//get scalar input
double args0=args(0).scalar_value();
double args1=args(1).scalar_value();
//get matrix input
Matrix y1(args(0).matrix_value());
Matrix y2(args(1).matrix_value());
//return value, use octave_value_list to return an array
//see: https://www.gnu.org/software/octave/doc/interpreter/Character-Strings-in-Oct_002dFiles.html#Character-Strings-in-Oct_002dFiles
return octave_value(args0);

Matrix Define:

myMx = zeros(2,4);
dim_vector dv(2);
dv(0) = 2; dv(1) = 4;
Matrix myMx(dv);//myMx.fill(0.0);


Matrix myMx(2, 4);

Get Matrix Dimensions:

dv = size(myMx);
dim_vector dv=myMx.dims();
//nrows = dv(0);//or myMx.rows();
//ncols = dv(1);//or myMx.cols();

Set/Get Matrix Elements:

myMx(1,2) = 1;
myMx(0,1) = 1;

Matrix Operations:

Mx_c = Mx_a+Mx_b;
Mx_c = Mx_a-Mx_b;
Mx_c = Mx_a * Mx_b;
Mx_c = Mx_a / Mx_b;
Mx_at = Mx_a';
Matrix Mx_c = Mx_a + Mx_b;
Matrix Mx_c = Mx_a - Mx_b;
Matrix Mx_c = Mx_a * Mx_b;
Matrix Mx_c = Mx_a * Mx_b.inverse();
Matrix Mx_at = Mx_a.transpose();

Get Matrix Rows/Columns:

Mx_j = Mx(j,:);
Mx_k = Mx(:,k);
RowVector Mx_j=Mx.row(j-1);
ColumnVector Mx_j=Mx.column(k-1);

Row Sums:

y1sum = sum(y1);
RowVector y1sum = (y1.sum()).row(0);

Matrices Combination:

int rows, columns, columns2;
// initialize row and column values to something
Matrix A(rows,columns);
Matrix B(rows,columns2);
// initialize the contents of A and B to something
A.insert(B, 0, columns);

Call Matlab Function:

stat = 0;
ipv = norcdf(stat);
double stat = 0.0;
octave_value_list ipv =  feval ("normcdf", octave_value(stat), nargout);
feval ("disp", octave_value(ipv(0).scalar_value()), nargout);

Self-Defined Functions:

tr = trace(Mx_a);
sm = sum(sum(Mx_a));
double trace(Matrix mx)
	double value = 0.0;
	dim_vector dv = mx.dims();
	int p = dv(0);
	for(int i=0;i<p;i++)//starts from 0 here!
		value += mx(i,i);
	return value;
double sumall(Matrix mx)
	double value = 0.0;
	dim_vector dv = mx.dims();
	int pr = dv(0);
	int pc = dv(1);
	for (int ir1 = 0 ; ir1 < pr ; ir1++)
		for (int ic1 = 0 ; ic1 < pc ; ic1++)
		return value;
double tr = trace(A);
double sm = sumall(A);

To compile the cpp source file:

mkoctfile octtest.cc;

Other Things to Note:
You can check the Cpp files folder (e.g.,”D:\Octave\Octave3.6.4_gcc4.6.2\include\octave-3.6.4\octave”) for more Octave .oct APIs.
C++ index starts from 0.
Integer divided by integer is integer in Cpp by default.

Source Code:

Some useful links:
http://wiki.octave.org/Tips_and_tricks#C.2B.2B (Octave C++ Quick Start)



GNU Octave Beginner’s Guide

Gnu Octave Version 3.0.1 Manual: A High-Level Interactive Language For Numerical Computations

Posted in matlab/octave, programming | Leave a comment

Blender for Data Visualization – Updated Code for Blender 2.6


Finally I updated the data visualization code from Blender 2.49 to Blender 2.6. You can find the updated Python source code here:

The Blender Python APIs changed a lot since Blender 2.5, but the biggest problem is lack of documents and tutorials at that time. Now things are much better than before. The changes from Blender 2.5 to 2.6 are not so significant as that from 2.4. Therefore, many resources for Blender 2.5 are also very helpful. Some very useful resources which help me to do the port are listed below:


Posted in data visualization, programming | Leave a comment

Blender for Data Visualization Part 2 – Render Output as Vector Graphics

In my last post, I showed you the basics to use Blender and Python for data visualization. Python is a powerful language which is very suitable for (large) data set processing, we may further explore it later. Blender is one of the most useful tools for 3D stuffs. In this post, I give a simple example of using Blender’s vector rendering function for a better output.

There are several vector render extensions for Blender. One of them I will use here is the “VRM” rendering script, which supports to export Flash’s “SWF” and “SVG” (Scalable Vector Graphics) file formats. Blender 2.49b should contain this extension. You can check that in your Blender’s menu (Render->VRM) or the script folder (“C:\Documents and Settings\Administrator\Application Data\Blender Foundation\Blender\.blender\scripts\py_render_249\vrm.py“). Anyway, I suggest you to download the latest release here:


Just download the snapshot, unzip, copy the “vrm.py” to Blender’s script folder and replace the old one.

Here is the snippet for “SVG” output, a simple modified version of the source code in my earlier post, so only the different parts are listed below:

import sys
sys.path.append( "C:\Documents and Settings\Administrator\Application Data\Blender Foundation\Blender\.blender\scripts\py_render_249" )
import vrm
from vrm import vectorize,ConsoleProgressIndicator,GraphicalProgressIndicator,progress,config
if __name__ == '__main__':
	csv = sys.argv[-1]#get the input file name
	if csv.endswith('.csv'):
		center = barchart(sys.argv[-1])

		#render and save the image

		config.edges['WIDTH'] = 1
		config.edges['COLOR'] = [255,0,0]
		config.edges['SHOW'] = True
		#config.edges['SHOW_HIDDEN'] = False
		#config.edges['STYLE'] = 'MESH'  # MESH or SILHOUETTE
		#config.edges['SHOW_HIDDEN'] = 1
		config.polygons['SHOW'] = False
		config.polygons['SHADING'] = 'FLAT'  # FLAT or TOON
		config.polygons['HSR'] = 'PAINTER'  # PAINTER or NEWELL
		config.polygons['EXPANSION_TRICK'] = True
		config.polygons['TOON_LEVELS'] = 2
		config.output['FORMAT'] = 'SVG'
		#config.output['ANIMATION'] = False
		#config.output['JOIN_OBJECTS'] = True

		vrm.progress = ConsoleProgressIndicator()
		#vrm.progress = GraphicalProgressIndicator()

		print 'input file has no .csv extension'
	print __name__

It simply imports the vrm functions and replaced the “PNG” image rendering part using the “vectorize” function provided by VRM. There are several lines for configuring the render, such as the edges and fill. You can find what those configure parameters mean by reading the source file “vrm.py” or by testing them directly in Blender through the user interface (Blender->Render->VRM), so I won’t explain them here. After execution, you can find the file “test.svg” created in your “C:/” folder; you can open and edit this file using softwares like Inkscape or GIMP, and convert it into “PS” and “PDF” formats.

To export the result as a SWF file, which can be easily embedded online, you need to install “libming” for Python first. You can find libming on its official site: http://www.libming.org/.
For windows users, precompiled installers can be downloaded here:
http://opensourcepack.blogspot.com/2011/03/pming-python-binding-for-libming.html, just choose the right version (2.5/2.6/2.7) for you Python and install. Here is the snippet for SWF output:

		config.edges['WIDTH'] = 2
		config.edges['COLOR'] = [0,0,0]
		config.edges['SHOW'] = False
		#config.edges['SHOW_HIDDEN'] = False
		#config.edges['STYLE'] = 'MESH'  # MESH or SILHOUETTE
		config.polygons['SHOW'] = True
		config.polygons['SHADING'] = 'TOON'  # FLAT or TOON
		config.polygons['HSR'] = 'PAINTER'  # PAINTER or NEWELL
		config.polygons['EXPANSION_TRICK'] = True
		config.polygons['TOON_LEVELS'] = 2
		config.output['FORMAT'] = 'SWF'
		#config.output['ANIMATION'] = False
		#config.output['JOIN_OBJECTS'] = True

		vrm.progress = ConsoleProgressIndicator()
		#vrm.progress = GraphicalProgressIndicator()

You can find the full source code here:


And the result:
Pretty nice, right?

Some useful links:

Blender 2.49 Scripting (Book)

Posted in data visualization, programming | 4 Comments

Central Limit Theorem, Normal Distribution and Entropy

A simple form of The Central Limit Theorem:

For identically independent distributed random vectors Y_{1},Y_{2},\cdots; let \bar{Y}_{n}=\Sigma_{i=1}^{n}Y_{i}/n.

If E||Y_{1}||^{2}\le\infty, \sqrt{n}(\bar{Y}_{n}-EY_{1})\to N(0,CovY_{1}).

When I first learn the central limit theorem, I’m always curious why the sum of a large number of random variables is distributed approximately as a normal/Gaussian distribution. The answer is simple, that’s where normal distribution comes from. Normal distribution is found when we first try to find the limitation of the sum of a large number of random variables. In other words, normal distribution is defined as the limit distribution of a large number of random variables.

Firstly, using the i.i.d. condition, it’s easy to check that E(\bar{Y}_{n}-EY_{1})=0 and Cov(\bar{Y}_{n}-EY_{1})=\frac{CovY_{1}}{n}.

Before we move on, let’s review the basic idea of entropy. The entropy of words s is H(s)=-\Sigma_{\omega\in s}p(\omega)logp(\omega). Entropy is a measure of uncertainty. The less information we have about something, the larger its entropy will be.

Now, ignore the central limit theorem for a while, let’s image what the limit distribution of a large number of random variables can be.

Note that \bar{Y}_{n} is a sample mean. The more we “mean”, the more individual characters will be lost. Mean will hide individual characters, make individuals indistinctive, and hence reduce the information we have. So it’s natural that the more we “mean”, the larger its entropy will be.

OK, that’s enough. From the point above, we can expect that the limit distribution has a large entropy. In other words, the normal distribution should have a large entropy.

The fact is that:

Theorem 8.6.5 ([1], P254)

Let the random vector X\in R^{n} have zero mean and covariance K=EXX^{t}. Then h(X)\leq\frac{1}{2}log(2\pi e)^{n}|K|, with equality iff X\sim N(0,K).

Theorem 8.8.6 ([1], P255, Estimation error and differential entropy)

For any random variable X and estimator \hat{X}, E(X-\hat{X})^{2}\geq\frac{1}{2\pi e}e^{2h(X)}, with equality if and only if X is Gaussian and \hat{X} is the mean of X.

Corollary: Given side information Y and estimation \hat{X}(Y), it follows that E(X-\hat{X}(Y))^{2}\geq\frac{1}{2\pi e}e^{2h(X|Y)}.

8a.6 ([2], P532) N_{p} as a Distribution with Maximum Entropy

The multivariate normal distribution N_{p} has the maximum entropy

(H=-\int P(U)logP(U)dv) subject to the condition that mean and covariance are fixed.

For given mean and covariance, the Gaussian distribution is the maximum entropy distribution. “It gives the lowest log likelihood to its members on average. That means Gaussian distribution is the safest assumption when the true distribution is unknown” ([3]). That also explains why many people tend to “abuse” Gaussian assumption so much.


[1] Elements Of Information Theory, Second Edition, Thomas M. Cover, Thomas M. Cover, 2006

[2] Linear Statistical Inference and its Applications, Second Edition, C. Radhakfushna Rao, 2002

[3] A Short Introduction to Model Selection, Kolmogorov Complexity and Minimum Description Length (MDL), Volker Nannen, 2003

Posted in asymptotic, entropy, probability | 3 Comments

Blender for Data Visualization

There are many tools and languages for visualizing data, such as R (1
, 2, 3), Processing, Flash (1, 2, 3). And the powerful open source 3D editor Blender is of course one of them. The advantage of using Blender is you can use the high-level scripting language Python and the handy Blender API.

Here is a simple example showing how to draw a bar-chart from a “.csv” file using Blender (Source code in Python, GPL, by Michel Anders, adapted from the “Chapter02/README.txt” of the source code package for the book “Blender 2.49 Scripting“):

__author__    = "Michel Anders (varkenvarken)"
__version__   = "1.0 2009/07/28"
__copyright__ = "copyright 2009,2010 Michel J. Anders."
__url__       = ["author's site, http://www.swineworld.org"]

#    This program is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    GNU General Public License for more details.
#    You should have received a copy of the GNU General Public License
#    along with this program.  If not, see <http://www.gnu.org/licenses/>.

from Blender import Scene,Mesh,Window,Camera,Lamp,Text3d,World
from csv import DictReader
import sys

from math import pi

def label(text,position,orientation='z'):
	draw a label
	if orientation=='x':
	elif orientation=='z':
	print 'label %s at %s along %s' %(text,position,orientation)

def bar(value,position,shape='Cube'):
	draw a 3D bar
	print 'bar'
	me = Mesh.Primitives.Cube(1.0)
	print 'mesh added'
	sc = Scene.GetCurrent()
	print 'got scene'
	ob = sc.objects.new(me,'Mesh')
	print 'object added'
	print 'location set'
	print '%s-bar height %s at %s' %(shape,value,position)

def barchart(filename):
	draw the bar-chart and return the center position of the chart
	csv = open(filename)
	data = DictReader(csv)#standard csv Python module to read the data

	xlabel = data.fieldnames[0]#column headers
	print xlabel
	rows = [d for d in data]#other data
	print rows
	#find the extremes of the data
	maximum = max([float(r[n]) for n in data.fieldnames[1:] for r in rows])
	minimum = min([float(r[n]) for n in data.fieldnames[1:] for r in rows])
	print maximum
	print minimum

	for x,row in enumerate(rows):
	    label(row[xlabel],(x,10,0))#draw row label
	    for y,ylabel in enumerate(data.fieldnames[1:]):
	        bar(10.0*(float(row[ylabel])-minimum)/maximum,(x,0,y+1))#draw the bar for each column of row x
	x = lastx+1
	for y,ylabel in enumerate(data.fieldnames[1:]):
	    label(ylabel,(x,0,y+0.5),'x')#draw column label


	return (lastx/2.0,5.0,0.0)#return the center position of the chart

def removeobjects():
	sc = Scene.GetCurrent()
	for ob in sc.objects:

def addcamera(center):
	sc = Scene.GetCurrent()
	ca = Camera.New('persp','Camera')
	ob = sc.objects.new(ca)

def addlamp(loc=(0.0,0.0,10.0)):
	sc = Scene.GetCurrent()
	la = Lamp.New('Lamp')
	ob = sc.objects.new(la)

if __name__ == '__main__':
	csv = sys.argv[-1]#get the input file name
	if csv.endswith('.csv'):
		center = barchart(sys.argv[-1])
		#render and save the image

		print 'input file has no .csv extension'
	print __name__

The result:

Download the source code: https://bzstat.googlecode.com/svn/trunk/DataVisualization/Blender/

And read the notes on how to run the script: https://bzstat.googlecode.com/svn/trunk/DataVisualization/Blender/ReadMe.txt


Blender 2.49 Scripting (Book)

Blender 2.49 Python API: http://www.blender.org/documentation/249PythonDoc/


Posted in data visualization, programming | 4 Comments

Using SDL and OpenGL in R

SDL(Simple DirectMedia Layer) is “A cross-platform multimedia library designed to provide fast access to the graphics framebuffer and audio device”. OpenGL is “The Industry Standard for High Performance Graphics”.

With rdyncall package, which provides "a cross-platform framework for dynamic binding of C libraries using a flexible Foreign Function Interface", we can use SDL and OpenGL in R!

Testing Notes:

1. Install rdyncall package from CRAN: http://cran.r-project.org/web/packages/rdyncall/index.html

2. Download required DLLs for the test from here: http://dyncall.org/bindings/rdyncall/rdyncall-demo-dlls.zip

Unzip it and copy the DLL files to the system path location (e.g. C:\windows\system32).

3. Test Code (Copied from “rdyncall/demo/SDL.R”)
In R:

# Package: rdyncall
# File: demo/SDL.R
# Description: 3D Rotating Cube Demo using SDL,OpenGL and GLU. (dynport demo)


# Globals.

surface <- NULL

# Init.

init <- function()
  err <- SDL_Init(SDL_INIT_VIDEO)
  if (err != 0) error("SDL_Init failed")
  surface <<- SDL_SetVideoMode(512,512,32,SDL_DOUBLEBUF+SDL_OPENGL)

# GL Display Lists

makeCubeDisplaylist <- function()
  vertices <- as.double(c(
  -1, 1,-1,
   1, 1,-1,
  -1,-1, 1,
   1,-1, 1,
  -1, 1, 1,
   1, 1, 1

  colors <- as.raw( col2rgb( rainbow(8) ) )

  triangleIndices <- as.integer(c(
    0, 2, 1,
    2, 3, 1,
    1, 3, 7,
    1, 7, 5,
    4, 5, 7,
    4, 7, 6,
    6, 2, 0,
    6, 0, 4,
    2, 7, 3,
    2, 6, 7,
    4, 0, 5,
    0, 1, 5

  glVertexPointer(3, GL_DOUBLE, 0, vertices )

  glColorPointer(3, GL_UNSIGNED_BYTE, 0, colors )

  displaylistId <- glGenLists(1)
  glNewList( displaylistId, GL_COMPILE )
  glDrawElements(GL_TRIANGLES, 36, GL_UNSIGNED_INT, triangleIndices)



# Mainloop.

mainloop <- function()
  displaylistId <- makeCubeDisplaylist()
  evt <- new.struct(SDL_Event)
  blink <- 0
  tbase <- SDL_GetTicks()
  quit <- FALSE
    tnow <- SDL_GetTicks()
    tdemo <- ( tnow - tbase ) / 1000


    aspect <- 512/512
    gluPerspective(60, aspect, 3, 1000)

    glRotated(sin(tdemo)*60.0, 0, 1, 0);
    glRotated(cos(tdemo)*90.0, 1, 0, 0);




    SDL_WM_SetCaption(paste("time:", tdemo),NULL)
    blink <- blink + 0.01
    while (blink > 1) blink <- blink - 1
    while( SDL_PollEvent(evt) != 0 )
      if ( evt$type == SDL_QUIT ) quit <- TRUE
      else if (evt$type == SDL_MOUSEBUTTONDOWN )
        button <- evt$button
        cat("button ",button$button," at ",button$x,",",button$y,"\n")
    glerr <- glGetError()
    if (glerr != 0)
      cat("GL Error:", gluErrorString(glerr) )
      quit <- 1
  glDeleteLists(displaylistId, 1)

cleanup <- function()

run <- function()


You will see a window showing a rotating cube:


SDL: http://www.libsdl.org/

OpenGL: http://www.opengl.org/

Dyncall: http://dyncall.org/index.shtml

Rdyncall Package: http://cran.r-project.org/web/packages/rdyncall/index.html

Posted in programming, r | 1 Comment