---
title: 'Interfacing C and R: The .Call() Interface'
author: "Jan de Leeuw"
date: "Version 004 08/25/2016"
output:
bookdown::gitbook:
config:
toc:
collapse: null
keep_md: yes
split_by: chapter
bookdown::pdf_book:
keep_tex: yes
latex_engine: pdflatex
number_sections: yes
toc: yes
toc_depth: 3
html_document:
toc: yes
toc_depth: '3'
description: null
documentclass: book
fontsize: 12pt
graphics: yes
link-citations: yes
mainfont: Times New Roman
cover-image: graphics/cover.jpg
site: bookdown::bookdown_site
bibliography: dotcall.bib
---
#Introduction
**Note: The directory [gifi.stat.ucla.edu/dotcall](http://gifi.stat.ucla.edu/dotcall) has
the complete Rmd file of this report with all code chunks. It also has pdf and epub versions, and R and C files with the code. The book and the files that go with it are in the public domain. Suggestions for improvement are always welcome. **

The model for using R (@r_core_team_16) in this report is:
* We input the data and store them in objects in R.
* Using `.Call()` we pass the R objects to C subroutines for computation.
* We get the results back as R objects from C.
* We use R again to generate output, including graphics.
This implies, or at least suggests, the following conventions.
* All I/O is handled in R.
* All type checking, and argument coercion is done in R.
* No `.C`, `.external()`, `C++`, and `Rcpp` in this report.
* Our C functions should not have side effects.
I do not use `Rcpp` (@eddelbuettel_13) because in my formative years I was taught to dislike `C++`, and, even after years of therapy, I never got over those negative feelings. This makes the material in this book somewhat old-fashioned, since the excellent and authorative new books by @wickham_15 and @chambers_16 clearly suggest using `Rcpp`. Thus my stubbornness is not related to any shortcomings of `Rcpp`. It merely illustrates that because I am already retired, I have nothing to prove and nothing to lose.
As for coding style in both R and C:
* spaces around binary operators except in x:y,
* braces for groups of statements (even if there is only a single statement in the group) with the closing brace is on its own line,
* no semi-colons, i.e no multiple statements on a single line,
* spaces after commas,
* spaces before opening parentheses, except in function or macro calls,
* no spaces for square brackets indexing vectors, matrices, or lists,
* lower case, possibly with underscores, for variable and function names.
See @wickham_15, chapter 5, for more details.
# Numbers
##Vectors
###Hypotenuse
We first give an example of a function returning a double, i.e. a vector of type numeric and of length one. For this we use the `hypot` function in the C math library, which takes two doubles $x$ and $y$ and returns the double $\sqrt{x^2+y^2}$. Note that
`R.h` includes `math.h`, which declares `hypot`.
```{r hypotC, engine = "c", results ="hide"}
#include
#include
SEXP my_hypot(SEXP x, SEXP y) {
int nprotected = 0;
SEXP z = PROTECT(allocVector(REALSXP, 1));
nprotected++;
REAL(z)[0] = hypot(asReal(x), asReal(y));
UNPROTECT(nprotected);
return z;
}
```
```{r hypotR}
.Call("my_hypot", as.double(3), as.double(4))
```
Some aspects of the code are worth explaining. We always include at least `R.h` and `Rinternals.h`. In our C functions we only pass objects of type **SEXP** and we always return an object of type **SEXP**. If our function does not return anything, we still force it to return `R_NilValue`, which is a particular **SEXP**. But of course functions that do not return anything are not very useful, because they cannot have side effects.
A **SEXP** is a pointer to a C structure of type **SEXPREC**, which stores information about basically all types of R objects. We never access the fields of the **SEXPREC** directly, we use the functions and macros defined in
`Rinternals.h` to work with **SEXP**'s. A **SEXPREC** has a **SEXPTYPE**, which is **REALSXP** for vectors of doubles, **VECSXP** for lists, and so on. We will work with several of the **SEXPTYPE**'s in our later examples.
If we define a **SEXP** in our code, i.e. if we declare it and reserve memory for it, then the memory has to be protected from R's garbage collector. In all our examples we
maintain the counter nProtected that indicates how many objects we have protected, and ultimately how many objects we have to free, using `UNPROTECT()`. We do this even if, as in this example, there is only a single object to protect. The SEXP in the argument list of the function are automatically protected, which is one reason we always treat them as read-only and never modify them.
For a SEXP of SEXPTYPE REALSXP we use the macro `REAL()` to get a pointer to the payload of the **SEXPREC**, i.e. to the first element of the corresponding vector of doubles, while `asReal()` converts the payload to a double, i.e. it dereferences `REAL()`. Note that we do not use any type checking, because we assume this has been taken care of in the R part. There are similar macros and functions for the other **SEXPTYPE**'s.
###Euclidean Norm
More generally, we can write a function that takes a vector $x$ of doubles and returns a double equal to $\sqrt{\sum_{i=1}^n x_i^2}$. This does not really add much to the first example, so just for fun we throw in a macro for the square of a number.
```{r enormC, engine = "c", results ="hide"}
#include
#include
#define square(x) ({ \
typeof(x) _x = (x); \
_x * _x; \
})
SEXP my_enorm(SEXP x) {
R_len_t n = length(x);
int nProtected = 0;
SEXP z = PROTECT(allocVector(REALSXP, 1));
nProtected++;
double *rZ = REAL(z), *rX = REAL(x);
rZ[0] = 0.0;
for (R_len_t i = 0; i < n; i++) {
rZ[0] += square(rX[i]);
}
rZ[0] = sqrt(rZ[0]);
UNPROTECT(nProtected);
return z;
}
```
```{r xlxR}
.Call("my_enorm", as.double(c(1:4)))
```
Note that we use type R_len_t for the number of elements in a standard (not a long) vector. This is currently just int, but it is *future-proof* because it could change to something longer.
###Vectorized f(x) = x * log(x)
```{r xlxC, engine = "c", results ="hide"}
#include
#include
SEXP xlx (SEXP x) {
if (! isReal(x))
error("function argument must be a double\n");
R_len_t n = LENGTH (x), i, nProtected = 0;
SEXP z = PROTECT (allocVector (REALSXP, n)); nProtected++;
double *h = REAL (z), *v = REAL (x);
for (i = 0; i < n; i++)
h[i] = v[i] < 1e-15 ? 0.0 : v[i] * log (v[i]);
UNPROTECT (nProtected);
return (z);
}
```
```{r enormR, fig.align = "center"}
x <- seq(0, 2, length = 101)
plot(x, .Call ("xlx", as.double (x)), type = "l", lwd = 2, col = "RED")
abline (h = 0)
```
###Decode and Encode
```{r aplC, engine = "c", results ="hide"}
#include
#include
SEXP APLDECODE( SEXP cell, SEXP dims )
{
int aux = 1, n = LENGTH (dims), nProtected = 0;
SEXP ind = PROTECT (allocVector( INTSXP, 1 ) ); nProtected++;
int *nind = INTEGER (ind), *ndims = INTEGER (dims), *ncell = INTEGER (cell);
nind[0] = 1;
for( int i = 0; i < n; i++ ) {
nind[0] += aux * ( ncell[i] - 1 );
aux *= ndims[i];
}
UNPROTECT(nProtected);
return (ind);
}
SEXP APLENCODE( SEXP ind, SEXP dims )
{
int n = LENGTH (dims), aux = asInteger (ind), pdim = 1, nProtected = 0;
SEXP cell = PROTECT(allocVector( INTSXP, n ) ); nProtected++;
int *ncell = INTEGER (cell), *ndims = INTEGER (dims);
for ( int i = 0; i < n - 1; i++ )
pdim *= ndims[i];
for ( int i = n - 1; i > 0; i-- ){
ncell[i] = ( aux - 1 ) / pdim;
aux -= pdim * ncell[i];
pdim /= ndims[i - 1];
ncell[i] += 1;
}
ncell[0] = aux;
UNPROTECT( nProtected );
return (cell);
}
```
```{r aplR}
for (i in 1:12)
print (.Call("APLENCODE", as.integer (i), as.integer (c (3, 4))))
for (i in 1:12)
print (.Call("APLENCODE", as.integer (i), as.integer (c (2, 2, 3))))
for (i in 1:12)
print (.Call("APLDECODE", as.integer (.Call("APLENCODE", as.integer (i), as.integer (c (3, 4)))), as.integer (c (3, 4))))
```
## Vectors with NA
If the vector can contain NA's we have to take that into account. The following routine squares the elements of a vector and returns the vector of squares. It skips over the NA's. We also allow for the possibility that some elements are not numbers. If possible, they are converted. If they cannot be converted, R issues a warning and inserts NA's.
```{r squareC, engine = "c", results ="hide"}
#include
#include
SEXP square (SEXP x) {
R_len_t i, n = length (x), nProtected = 0;
SEXP y = PROTECT (allocVector (REALSXP, n)); nProtected++;
double *v = REAL (y), *z;
if (TYPEOF (x) == REALSXP) {
z = REAL (x);
} else {
z = REAL (coerceVector (x, REALSXP));
}
for (i = 0; i < n; i++) {
v[i] = R_NaReal == z[i] ? -1 : z[i] * z[i];
}
UNPROTECT (nProtected);
return (y);
}
```
```{r squareR}
.Call("square", c(1, NA, 2, "3", -1))
.Call("square", c(1, NA, 2, "a", -1))
```
##Matrices and Arrays
###Passing matrix from R to C
###Passing matrix from C to R
##Attributes
#Characters
##From R to C
##From C to R
#Lists
##From R to C
##From C to R
#Functions and Expressions
##Calling an R function from C
```{r normC, engine = "c", results ="hide"}
#include
#include
SEXP my_rnorm (SEXP n) {
R_len_t m = asInteger(n);
int nprotected = 0;
SEXP y = PROTECT(allocVector(REALSXP, m));
nprotected++;
SEXP rnorm_r = PROTECT(install("rnorm"));
nprotected++;
y = eval(lang2(rnorm_r, ScalarInteger(m)), R_GlobalEnv);
UNPROTECT(nprotected);
return y;
}
```
```{r normR}
.Call("my_rnorm", as.integer(5))
```
```{r firstNC, engine = "c", results ="hide"}
#include
#include
SEXP myFirstN(SEXP n) {
R_len_t m = asInteger (n);
int nprotected = 0;
SEXP y = PROTECT(allocVector(REALSXP, m));
nprotected++;
SEXP firstR = PROTECT(install("firstN"));
nprotected++;
y = eval(lang2(firstR, ScalarInteger (m)), R_GlobalEnv);
UNPROTECT (nprotected);
return (y);
}
```
```{r firstNR}
firstN <- function (n) 1:n
.Call("myFirstN", as.integer(5))
```
##Passing R function to C
# Vectors of Characters
#References