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.

### Like this:

Like Loading...

*Related*