One reason I like using R for data analysis is that R has a great collection of packages that let you easily apply state-of-the-art methods to your problems. But once in a while you find a library that you would like to use that does not have a R wrapper, yet. While the great Rcpp package provides a convenient way to write R extensions in C++, it obviously requires you to write C++ code and to have a compiler installed.

An alternative I found about only recently is the rdyncall package. *rdyncall* provides an improved Foreign Function Interface (FFI), which allows you to dynamically invoke C libraries.

In this blog post I want to give you an example on how to employ *rdyncall* to use

the LWPR libary for *Locally Weighted Projection Regression*.

## Function signatures

The minimal set of functions we need from the LWPR library consists of `lwpr_init_model()`

to initialize a new model, `lwpr_update() `

to train the model and `lwpr_predict()`

to get predictions from the model.

In the ` lwpr.h`

header file you will find the following definitions for these functions:

int lwpr_init_model(LWPR_Model *model, int nIn, int nOut, const char *name); int lwpr_update(LWPR_Model *model, const double *x, const double *y, double *yp, double *max_w); void lwpr_predict(const LWPR_Model *model, const double *x, double cutoff, double *y, double *conf, double *max_w);

In order to use these functions with rdyncall we have to translate the C function definitions into rdyncall’s function signature syntax.

The syntax is `function_name(argument_types)return_type;`

. Most types are respresented by single letters: i = integers, d = double, Z = const char*, … Pointer types are prepended by an asterisk (*), so `*d `

stands for `double*`

in C. Custom types (structs) can also be used and pointers to structs are written as `*`

. The package documentation provides a complete list of available type notations. For our three functions above the signatures for rdyncall look like this:

lwpr_init_model(*<LWPR_Model>iiZ)i; lwpr_update(*<LWPR_Model>*d*d*d*d)i; lwpr_predict(*<LWPR_Model>*dd*d*d*d)v;

## Struct definition

In the function signatures we used the struct *LWPR_Model*. In order to be able to create and modify LWPR_Model “objects” we also have to specify the memory layout of the C struct. The notation for struct definitions is `StructName{field_tyoes}fieldname1 fieldname2 ...;`

.

The C definition of LWPR_Model is:

typedef struct LWPR_Model { int nIn; /**< \brief Number N of input dimensions */ int nInStore; /**< \brief Storage-size of any N-vector, for aligment purposes */ int nOut; /**< \brief Number M of output dimensions */ int n_data; /**< \brief Number of training data the model has seen */ double *mean_x; /**< \brief Mean of all training data the model has seen (Nx1) */ double *var_x; /**< \brief Mean of all training data the model has seen (Nx1) */ char *name; /**< \brief An optional description of the model (Mx1) */ int diag_only; /**< \brief Flag that determines whether distance matrices are handled as diagonal-only */ int meta; /**< \brief Flag that determines wheter 2nd order updates to LWPR_ReceptiveField.M are computed */ double meta_rate; /**< \brief Learning rate for 2nd order updates */ double penalty; /**< \brief Penalty factor used within distance metric updates */ double *init_alpha; /**< \brief Initial learning rate for 2nd order distance metric updates (NxN) */ double *norm_in; /**< \brief Input normalisation (Nx1). Adjust this to the expected variation of your data. */ double *norm_out; /**< \brief Output normalisation. Adjust this to the expected variation of your output data. */ double *init_D; /**< \brief Initial distance metric (NxN). This often requires some tuning (NxN) */ double *init_M; /**< \brief Cholesky factorisation of LWPR_Model.init_D (NxN) */ double w_gen; /**< \brief Threshold that determines the minimum activation before a new RF is created. */ double w_prune; /**< \brief Threshold that determines above which (second highest) activation a RF is pruned. */ double init_lambda; /**< \brief Initial forgetting factor */ double final_lambda; /**< \brief Final forgetting factor */ double tau_lambda; /**< \brief This parameter describes the annealing schedule of the forgetting factor */ double init_S2; /**< \brief Initial value for sufficient statistics LWPR_ReceptiveField.SSs2 */ double add_threshold;/**< \brief Threshold that determines when a new PLS regression axis is added */ LWPR_Kernel kernel; /**< \brief Describes which kernel function is used (Gaussian or BiSquare) */ int update_D; /**< \brief Flag that determines whether distance metric updates are performed (default: 1) */ LWPR_SubModel *sub; /**< \brief Array of SubModels, one for each output dimension. */ struct LWPR_Workspace *ws; /**< \brief Array of Workspaces, one for each thread (cf. LWPR_NUM_THREADS) */ double *storage; /**< \brief Pointer to allocated memory. Do not touch. */ double *xn; /**< \brief Used to hold a normalised input vector (Nx1) */ double *yn; /**< \brief Used to hold a normalised output vector (Nx1) */ } LWPR_Model;

For dyncall it has to be (manually) transformed to:

LWPR_Model{iiii*d*d*ciidd*d*d*d*d*ddddddddiipp*d*d*d}nIn nInStore nOut n_data mean_x var_x name diag_only meta meta_rate penalty init_alpha norm_in norm_out init_D init_M w_gen w_prune init_lambda final_lambda tau_lambda init_S2 add_threshold kernel update_D sub ws storage xn yn

For fields with pointers to other custom types, such as `LWPR_Workspace`

, we can safely write `p`

for a generic pointer. The enum `LWPR_Kernel`

is represented by an integer (i).

## Putting things together

require("rdyncall") dynbind( "lwpr", # library name; rdyncall will look for liblwpr.so or lwpr.dll on the library path "lwpr_init_model(*<LWPR_Model>iiZ)i; lwpr_update(*<LWPR_Model>*d*d*d*d)i; lwpr_predict(*<LWPR_Model>*dd*d*d*d)v; lwpr_duplicate_mode(*<LWPR_Model>*<LWPR_Model>)i; lwpr_set_init_alpha(*<LWPR_Model>d)i; lwpr_set_init_D(*<LWPR_Model>*dd)i; lwpr_set_init_D_diagonal(*<LWPR_Model>*d)i; lwpr_set_init_D_spherical(*<LWPR_Model>d)i; lwpr_write_xml(*<LWPR_Model>Z)i; lwpr_read_xml(*<LWPR_Model>Z*i)i; ") # create struct type LWPR_Model parseStructInfos("LWPR_Model{iiii*d*d*ciidd*d*d*d*d*ddddddddiipp*d*d*d}nIn nInStore nOut n_data mean_x var_x name diag_only meta meta_rate penalty init_alpha norm_in norm_out init_D init_M w_gen w_prune init_lambda final_lambda tau_lambda init_S2 add_threshold kernel update_D sub ws storage xn yn")

This code will dynamically load the LWPR library and will create a function in the global environment for each given function signature. The `parseStructInfos`

invocation will create a variable LWPR_Model that can be used with `new.struct()`

to create new LWPR_Model objects.

These 20 lines of R code are enough to put everything in place. We are now able to use the LWPR library in R :-).

## The LWPR tutorial code in R

This article would not be complete without an example that actually uses the LWPR library.

So here is the

LWPR tutorial code in R:

testfunc = function(x) { 10 * sin(7.8*log(1+x)) / (1 + 0.1*x**2) } Ntr = 500 Xtr = 10 * runif(Ntr) Ytr = sapply(Xtr, testfunc) + 0.1 * rnorm(Ntr) * Xtr model = new.struct(LWPR_Model) lwpr_init_model(model, 1, 1, "Tutorial") lwpr_set_init_D(model, 20, 0) model$update_D = 1 model$diag_only = TRUE model$penalty = 0.0001 lwpr_set_init_alpha(model, 40) for (i in 1:20) { idx = sample(Ntr,Ntr) for (j in idx) { yp=0 lwpr_update(model, Xtr[j], Ytr[j], yp, NULL) } } Ntest = 500; Xtest = seq(0,10,length.out=Ntest) Ytest = numeric(Ntest) Conf = numeric(Ntest) for (k in 1:Ntest) { yp = 0 conf = 0 lwpr_predict(model, Xtest[k], 0.0, yp, conf, NULL); Ytest[k] = yp Conf[k] = conf } plot(Ytr ~ Xtr, pch=19, cex=0.5, col='red') lines(Xtest , Ytest, col="blue") lines(Xtest , Ytest + Conf, col="cyan") lines(Xtest , Ytest - Conf, col="cyan")

## Conclusion

`rdyncall`

allows you to call C library in R without writing any wrapper code. I find it quicker and easier than the alternatives, but it still requires a little knowledge of C to write the function and struct definitions based on the C header files. In the R Journal Article (PDF) about rdyncall a automatic conversion from C header files to rdyncall specifications is mentioned, but I think it was not released.

## See Also

- Adler, D.: Foreign Library Interface, R Journal Vol.4/1, June 2012 (PDF)
- Complete example: https://gist.github.com/felixr/5564929

Pretty nice, thanks for the post!