1 Rpk’s devolopment

  • There are many methods to develop an R pacakge

1.1 extension function

1.2 R call C/C++

  • Build foo.dll on Windows or foo.so on MacOS/LinuxOS
R CMD SHLIB foo.c

1.2.1 Examples

  • An example in VectProduce.c:
#include <math.h>
#include <string.h>
#include "Rinternals.h"
#include "R_ext/Rdynload.h"
#include <R.h>
#include <R_ext/Applic.h>

void VectProduce(double *A, int *n1_, double *v, int *n2_, double *outMatrix)
{
    int i,j,n1, n2;
    n1 = n1_[0];
    n2 = n2_[0];
    for (i=0;i<n1;i++){
        for(j=0;j<n2;j++)   outMatrix[j*n1+i]=A[i]*v[j];
    }
}

static const R_CMethodDef cMethods[] = {
  {"VectProduce", (DL_FUNC) &VectProduce, 5},
  NULL
};

void R_init_VectProduce(DllInfo *info)
{
  R_registerRoutines(info,cMethods,NULL,NULL,NULL);
}
  • Or in VectProduce1.c, see Rpk “RidgeVar”:
#include <R.h>
#include <Rinternals.h>  

void VectP(double *A, int n1, double *v, int n2, double *outMatrix)
{
    int i,j;
    for (i=0;i<n1;i++){
        for(j=0;j<n2;j++)   outMatrix[j*n1+i]=A[i]*v[j];
    }
}

SEXP VectProduce(SEXP x, SEXP y)
{
    int nx = length(x), ny = length(y);
    double *A = REAL(x), *B = REAL(y), *Y;
    SEXP outmatrix;

    PROTECT(outmatrix = allocMatrix(REALSXP, nx, ny));  
    Y = REAL(outmatrix);
    VectP(A, nx, B, ny, Y);

    UNPROTECT(1);  
    return(outmatrix);
}
  • or VectProduce2.cpp, see Rpk “IVGC”:
#include <R.h>
#include <Rinternals.h>  
extern "C"{
void VectP(double *A, int n1, double *v, int n2, double *outMatrix)
{
    int i,j;
    for (i=0;i<n1;i++){
        for(j=0;j<n2;j++)   outMatrix[j*n1+i]=A[i]*v[j];
    }
}
}
  • Or return a list in Test_list.c, see Rpk “plvs”:
#include <R.h>
//#include <Rinternals.h>  
#include <Rdefines.h>

SEXP setList() {
   int *p_myint, i;                       // variable declarations
   double *p_double;
   SEXP mydouble, myint, list, list_names;
   char *names[2] = {"integer", "numeric"};

   PROTECT(myint = NEW_INTEGER(5));       // creating an integer vector
   p_myint = INTEGER_POINTER(myint);
   PROTECT(mydouble = NEW_NUMERIC(5));    // ... and a vector of real numbers
   p_double = NUMERIC_POINTER(mydouble);
   for(i = 0; i < 5; i++) {
     p_double[i] = 1/(double)(i + 1);
     p_myint[i] = i + 1;
   }

   // a character string vector of the "names" attribute of the objects in our list
   PROTECT(list_names = allocVector(STRSXP, 2));
   for(i = 0; i < 2; i++)
      SET_STRING_ELT(list_names, i,  mkChar(names[i]));
 
   PROTECT(list = allocVector(VECSXP, 2)); // Creating a list with 2 vector elements
   SET_VECTOR_ELT(list, 0, myint);         // attaching myint vector to list
   SET_VECTOR_ELT(list, 1, mydouble);      // attaching mydouble vector to list
   setAttrib(list, R_NamesSymbol, list_names); //and attaching the vector names
   UNPROTECT(4);
   return (list);
}
  • The corresponding vp.R calling above three foo.dll on Windows (or foo.so on MacOS or LinuxOS):
VectProduce <- function(x, y)
{
    dyn.load("VectProduce.dll") 
    X <- matrix(0, nrow = length(x), ncol = length(y))
    n = length(x)
    p = length(y)
    out <- .C("VectProduce", x, as.integer(n), y, as.integer(p), X)

    return(list( VP = out[[5]] ))
}


VectProduce1 <- function(x, y)
{
    dyn.load("VectProduce1.dll")  
    out <- .Call("VectProduce", x, y)
    out
}

VectProduce2 <- function(x, y)
{
    dyn.load("VectProduce2.dll")  
    out <- .C("VectProduce", x, y)
    out
}

VectProduce <- function(x, y)
{
    dyn.load("Test_list.dll") 
    out1 <- .Call("out", x, y)
    out1
}

1.2.2 How to develop Rpk?

  • Rstudio
  • on R command-line:
    package.skeleton("myrpk")

1.2.3 How to pack Rpks?

  • on cmd (windows) or terminal (MacOS/Linux):
    R CMD build myrpk
    R CMD Rd2pdf myrpk
    R CMD check myrpk

1.2.4 Where and How to share Rpks?

1.2.5 Example: FactSum

1.3 Call foo.c in Rpks

  • Put all C/C++ files (eg. VectProduce.c) in sub-directory “src” of Rpg
  • Remember to include R.h in the header of foo.c

  • For example, see also Rpk “FactSum”

     #include <math.h>
 #include <string.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <R.h>
 #include <Rinternals.h>

1.4 Call executable files in Rpk

  • executable files (PIexe-unix.exe, or PIexe-linux.exe, or PIexe-windows.exe) in Rpk named “PI”
  • Put all executable files (eg. PIexe.exe) in sub-directory “exec” or “inst/exec” of Rpg
Pi <- function(n=100){
  if(n<1) stop("n must be a non negative integer!")
  progname <- paste("PIexe-",.Platform$OS.type,".exe",sep="")
  progpath <- system.file("exec",progname,package="PI")
  
  if (.Platform$OS.type=="unix") {
      pi_str = system2(progpath,args=as.character(n),stdout=TRUE)
  }
  pi <- paste0("3",pi_str,collapse = "")
  pi
}

1.4.1 Drawbacks:

  • Require different platform’s executable files
  • It is portable on cran: “A source package if possible should not contain binary executable files: they are not portable, and a security risk if they are of the appropriate architecture.”
  • But github is portable.

1.4.2 Advantages:

  • Can call the third C/C++ packages, such as “gmp”.
  • Can hide your secret and important codes

1.4.3 An example can be found in PI

1.5 Rcpp

1.5.1 How to use Rcpp?

1.5.2 An example: tensorMam

  • Include Rcpp.h
//[[Rcpp::depends(RcppEigen)]]
#include <Rcpp.h>
#include <RcppEigen.h>
#include <cmath>
#include <Eigen/Dense>
#include <iostream>
#include <Eigen/SVD>
#include <Eigen/Cholesky>
#include <Eigen/LU>
#include <list>
using namespace Rcpp;
using namespace Eigen;

//----------------------------------------------------------------**
//***----------------------------sequece--------------------------**
// [[Rcpp::export]]
VectorXd SEQ(double min, double max, double by)
{
    int leng = static_cast<int>(floor((max - min) / by + pow(3, -12)) + 1);
    VectorXd C = VectorXd::Constant(leng, 0);
    for (int i = 0; i < leng; i++)
        C[i] = min + by * i;
    return C;
}
//----------------------------------------------------------------**
//***----------------------uppertriangleNorm1 of AA---------------**
double uppertriangleNorm1AA(MatrixXd A)
{
    int i, j, k, p, ii;
    double norm1 = 0, temp1, temp2;
    p = A.cols();

    for (k = 0; k < p; k++) {
        temp2 = 0;
        for (j = 0; j < p; j++) {
            temp1 = 0;
            ii = k < j ? k : j;
            for (i = 0; i <= ii; i++) temp1 += A(i, j) * A(i, k);
            temp2 += fabs(temp1);
        }
        if (temp2 > norm1) norm1 = temp2;
    }
    return norm1;
}

//----------------------------------------------------------------**
//***--------------------transfer modal of unfoldings-------------**
// [[Rcpp::export]]
MatrixXd TransferModalUnfoldings(MatrixXd S, int d1, int d2, int r1, int r2, int r3)
{
  //From S_(d1) to S_(d2)
  int j;
  MatrixXd S1,S_3;
  if (d1 == 3) {
    if (d2 == 1){
      S1 = S.row(0).transpose();
      S1.resize(r1, r2); 
      for (j = 1; j < r3;j++) {
        S_3 = S.row(j).transpose();
        S_3.resize(r1, r2);
        S1 = cbind_rcpp(S1, S_3);// S3 is r3 *(r2r1) matrix
      }
    }   
    if (d2 == 2) {
      S1 = S.block(0, 0, r3, r1).transpose();
      S1.resize(1,r1*r3);
      for (j = 1; j < r2; j++) {
        S_3 = S.block(0, j*r1, r3, r1).transpose();
        S_3.resize(1,r1*r3);
        S1 = rbind_rcpp(S1, S_3);
      }
    }
  }
  
  if (d1 == 2) {
    if (d2 == 1) {
      S1 = S.block(0, 0, r2, r1).transpose();
      for (j = 1; j < r3; j++) {
        S1 = cbind_rcpp(S1,S.block(0,j*r1,r2,r1).transpose());
      }
    }
    if (d2 == 3) {
      S1 = S.block(0, 0, r2, r1).transpose();
      S1.resize(1, r2*r1);
      for (j = 1; j < r3; j++) {
        S_3 = S.block(0, j*r1, r2, r1).transpose();
        S_3.resize(1, r2*r1);
        S1 = rbind_rcpp(S1, S_3);
      }
    }
    
  }
  
  if (d1 == 1) {
    if (d2 == 2) {
      S1 = S.block(0, 0, r1, r2).transpose();
      for (j = 1; j < r3; j++) {
        S1 = cbind_rcpp(S1, S.block(0, j*r2, r1, r2).transpose());
      }
    }
    if (d2 == 3) {
      S1 = S.block(0, 0, r1, r2);
      S1.resize(1, r1*r2);
      for (j = 1; j < r3; j++) {
        S_3 = S.block(0, j*r2, r1, r2);
        S_3.resize(1, r1*r2);
        S1 = rbind_rcpp(S1, S_3);
      }
    }
  }
  return S1;
}
  • After insert foo.cpp into the sub-directory “src”, we need following steps
    • On R command type: compileAttributes("tensorMam")
    • add using namespace Eigen; in to file RcppExports.cpp, which is in the sub-directory “src”
    • On R command type: system("R CMD check tensorMam"))

1.5.3 Useful functions

  • sourCpp()
    • An example
    sourceCpp("foo.cpp")
  • Rcpp:::SHLIB() to produce an executable file foo.exe
    • On R command
    Rcpp:::SHLIB("foo.cpp")
    • On Shell
    Rscript -e "Rcpp:::SHLIB('foo.cpp')"

2 Matlab Calling C/C++

  • mexFunction()
  • An example VectProduce.c
#include "mex.h" 
void VectProduce(double *A, double *v, int n, double *outMatrix)
{
    int i,j;
    for (i=0;i<n;i++)
    {
        for(j=0;j<n;j++)
        outMatrix[i*n+j]=A[i]*v[j];
    }
}

void mexFunction(
 int nlhs,
 mxArray *plhs[],
 int nrhs,
 const mxArray *prhs[]) 
{
    double *inData, *outData, *Vector;
    int M, dim;
    inData = mxGetPr(prhs[0]);
    M = mxGetM(prhs[0]);
    Vector = mxGetPr(prhs[1]);
    dim = mxGetM(prhs[1]);
    if (M!=dim)
    {
        mexErrMsgTxt("The dimesions of two vectors does not match!");
    }
    plhs[0] = mxCreateDoubleMatrix(M, M, mxREAL);
    outData = mxGetPr(plhs[0]);
    VectProduce(inData, Vector, M, outData);
}
  • Mex command on Matlab

  • Select C/C++ compiler

    mex -setup
mex foo.c

  • Another example

  #include <stdlib.h>
  #include <math.h>
  #include <stdio.h>
  #include "mex.h" 

  int UpTriangularInv(double *B, int n, double *A)
  { 
      int i,j,k;
      for(i=0;i<n;i++)
          if(A[i*n+i]==0) return(0);
    for(i=n-1;i>=0;i--)//rows
    {
        if(A[i*n+i]!=1)
            for(j=i;j<n;j++)
                B[j*n+i] = B[j*n+i]/A[i*n+i];
        if(i>0)
        {
            for(j=i;j<n;j++)// columns
                for(k=0;k<i;k++)// rows
                    B[j*n+k] = B[j*n+k] - A[i*n+k]*B[j*n+i];
        }
    }
    return(1);
  }

  int LowTriangularInv(double *B, int n, double *A)
  { 
      int i,j,k;
      for(i=0;i<n;i++)
          if(A[i*n+i]==0) return(0);
    for(i=0;i<n;i++)//rows
    {
        if(A[i*n+i]!=1)
            for(j=0;j<=i;j++)
                B[j*n+i] = B[j*n+i]/A[i*n+i];
        if(i<n-1)
        {
            for(j=0;j<=i;j++)// columns
                for(k=i+1;k<n;k++)// rows
                    B[j*n+k] = B[j*n+k] - A[i*n+k]*B[j*n+i];
        }
    }
    return(1);
  }

  void mexFunction(
 int nlhs,
 mxArray *plhs[],
 int nrhs,
 const mxArray *prhs[]) 
{
    double *inData, *outData, flag, *B, IsUptriang;
    int M, N, i;
    inData = mxGetPr(prhs[0]);
    M = mxGetM(prhs[0]);
    N = mxGetM(prhs[0]);
    if(M!=N) mexErrMsgTxt("Input must be a square matrix.");
    if(nrhs<2)
    {
        IsUptriang =1;
        mexPrintf("Is input a uptriangular matrix!");
    }
    else
        IsUptriang = mxGetScalar(prhs[1]);

    plhs[0] = mxCreateDoubleMatrix(M, M, mxREAL);
    outData = mxGetPr(plhs[0]);
    for(i=0;i<M;i++) outData[i*M+i] = 1;
    if(IsUptriang==1)
        flag = UpTriangularInv(outData, M, inData);
    else
        flag = LowTriangularInv(outData, M, inData);

    if(flag==0) mexErrMsgTxt("Diag of input matrix must be nonzero.");
}
  • Matlab produce an executable file foo.exe

  • mcc function

    mcc -m VectProduce.m

  • Or foo.dll file

    mcc -l VectProduce.m

3 Python calling C/C++

  • There are many methods to call C function from Python. Here, we introduce how to achieve this via the C API.

3.1 The directory tree

-First of all, we create some files and subdirectories. The directory tree looks like the follows.

.
├── README.md
├── __init__.py
├── setup.py
├── LICENSE.txt
└── src
    ├── cfunction.c
    ├── cfunction.h
    └── _cfunction.c

3.2 Prepare the C codes and header files

  • The function below returns a vector that is an array in C, and stores the given numbers. Here, we save this function in a C file named cfunction.c.
#include "cfunction.h"
#include <math.h>
#include <stdlib.h>

void mycfunc(double *x, int *y, double constant1, int constant2, double *output) {
    output[0] = x[0];
    output[1] = y[0] * 1.0;
    output[2] = constant1;
    output[3] = constant2 * 1.0;
}
  • Then, a header file is required, say cfunction.h.
void mycfunc(double *x, int *y, double constant1, int constant2, double *output);

3.3 Other requirements for Wrapping Python pkg

  • To wrap it into Python pkg, a C file named _cfunction.c is required:
#include <Python.h>
#include "cfunction.h"
#include <numpy/arrayobject.h>
#include <stdlib.h>
#include <stdio.h>

static PyObject *py_mycfunc(PyObject *self, PyObject *args) {
    
    PyObject *py_x, *py_y;
    double number_a;
    int number_b;

    if (!PyArg_ParseTuple(args,"OOdi", &py_x, &py_y, &number_a, &number_b)) {
        return NULL;
    }

    // Interpret the input objects as numpy arrays with contiguous storage. 
    PyObject *py_x_array = PyArray_FROM_OTF(py_x, NPY_DOUBLE, NPY_IN_ARRAY);
    PyObject *py_y_array = PyArray_FROM_OTF(py_y, NPY_INT, NPY_IN_ARRAY);

    // If that did not work, throw an exception. 
    if (py_x_array == NULL || py_y_array== NULL ) {
        Py_XDECREF(py_x_array);
        Py_XDECREF(py_y_array);
        return NULL;
    }

    // Dimensions.
    int nx, px, ny;
    nx = (int)PyArray_DIM(py_x_array, 0);
    px = (int)PyArray_DIM(py_x_array, 1);
    ny = (int)PyArray_DIM(py_y_array, 0);

    if ( nx !=  ny ) {
        PyErr_SetString(PyExc_RuntimeError, "Dimensions don't match!!");
        return NULL;
    }

    // Get pointers to the data as C-types.
    double *x = (double*) PyArray_DATA(py_x_array);
    int    *y = (int*)    PyArray_DATA(py_y_array);

    // Maybe you need to create an output array.
    double *c_output;
    c_output = (double*)malloc(sizeof(double)*ny);

    // Call the external C function.
    mycfunc( x, y, number_a, number_b, c_output );

    // Create a numpy array and return it.
    npy_intp dims[1] = {ny};
    PyObject *ret = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, c_output);

    // Clean up. 
    free(c_output);
    Py_DECREF(py_x_array);
    Py_DECREF(py_y_array);

    return ret;
}

/* Documentations */
static char module_docs[] = "Some thing about your module.";

static char cfunction_docs[] = "Some thing about your C function.";


/* Module method table */
static PyMethodDef MydemoMethods[] = {
  {"mycfunc",  py_mycfunc, METH_VARARGS, cfunction_docs},
  /* If your extension module has many export functions, you need write interfaces for each of them. */
  /* And then add them to Module method table. For example: */
  /* {"mycfunc2",  py_mycfunc2, METH_VARARGS, cfunction2_docs}, */
  { NULL, NULL, 0, NULL}
};

/* Module structure */
static struct PyModuleDef demomodule = {
    PyModuleDef_HEAD_INIT,
    "Demo",           /* name of module */
    module_docs,        /* Doc string (may be NULL) */
    -1,                 /* Size of per-interpreter state or -1 */
    MydemoMethods       /* Method table */
};

/* Module initialization function */
PyMODINIT_FUNC PyInit_Demo(void) {
    PyObject *object = PyModule_Create(&demomodule);
    if(object == NULL) {
        return NULL;
    }
    import_array();
    return object;
}
  • Create an empty python file named __init__.py.
  • And create a python file named setup.py:
from setuptools import setup, Extension, find_packages
import numpy

extensions = Extension('Demo',  # name of our extension module
                        sources = ['src/cfunction.c', 'src/_cfunction.c'] )

setup(
    name='mydemo', # A name show on Pypi.
    version='0.1', 
    ext_modules = [extensions],
    license='MIT',
    description='An example python package',
    long_description=open('README.md').read(),
    install_requires=[],
    url='https://github.com/',
    author='someone',
    author_email='myemail@example.com',
    include_dirs = [numpy.get_include()]
)

3.4 Compile pkg

  • Type following codes on command line to compile pkg
# cd to your folder
python setup.py build   
python setup.py sdist   // Create a source distribution.
python setup.py install

3.5 Test the kpg

  • At the end, it will be done by running the following python test file:
import numpy as np
import Demo

n  = 4
p  = 5
X  = []
for _ in range(n):
    X.append(np.random.normal(size= p))
x = np.asarray(X)

y = np.random.choice(10, p)

Demo.mycfunc(x, y, 3.14, 2020)  // Trigger the error we set.

y = np.random.choice(10, n)

Demo.mycfunc(x, y, 3.14, 2020)
  • More details can be found mydemo.

4 Julia calling C/C++

  • It is illustrated later