Calling Apophenia library from R

I struggled to get the code (markov.c) on p.129 of the book “Modeling with Data” to work.

It turns out that the code needs to be modified to the following (the author of the book tells me that the current version of Apophenia no longer uses spaces as delimiters as a default option). Also, you may want to download the current version of the Apophenia library (available here), as there was a related bug that’s been fixed.

#include <apop.h>
int main(){
    apop_data *t = apop_text_to_data("data-markov", .has_row_names='n', 
		.has_col_names='n', .delimiters=" ");
    apop_data *out  = apop_dot(t, t);
    apop_data_show(out);
}

However, before discovering the source of the issue, I tried to determine if it was an issue with reading the data or an issue with the subsequent code. As such, I tried bringing in the data using R then sending it over to C.

The following requires the Apophenia library (available here) and the GSL library too.

Here is the C code (I saved this in markov2.c).

#include <apop.h>

void dot_prod(double *t, int *n, int *k, double *out) {
	// Create a gsl_matrix view on the data, then make this into
	// an apop_data structure that can be passed to apop_dot()
	gsl_matrix_view t_gsl = gsl_matrix_view_array (t, *n, *k);
	apop_data *t_apop = apop_matrix_to_data(&t_gsl.matrix);
	apop_data *out_apop = apop_dot(t_apop, t_apop, 0, 0);

	// Populate matrix passed back to R
	for (int i=0; i < *n; i++) {
		for (int j=0; j < *n; j++) {
			out[i * *n + j] = apop_data_get(out_apop, i, j);
		}
	}
}

Here is the R code:

# Read in data
setwd("~/Dropbox/academic/mwd2")
tt <- read.delim("data-markov", header=FALSE, as.is=TRUE, sep=" ")
tt <- as.matrix(tt)

# Compile and load C function dot_prod()
system('PKG_LIBS="-lapophenia -lgsl" R CMD SHLIB markov2.c')
dyn.load("markov2.so")

# Function that calculates the product XX of a square matrix X.
# Wrapper for C function loaded above.
dot_prod <- function(x) {
	n <- dim(x)[1]
	k <- dim(x)[2]
	if (n==k) {
	  z <- .C("dot_prod", as.double(x), n, k, results=double(n^2))
	  return(matrix(z$results, nrow=n))
	} else {
		return(NA)
	}
}

print(dot_prod(tt))

This confirmed the issue to be one in reading in the data.

Advertisements
This entry was posted in R, Uncategorized. Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s