Others Calling C/C++

Sat, Mar 14, 2020

Rpk’s devolopment

There are many methods to develop a R-pacakge

extension function

An introduction of Writing R Extensions

R call C/C++

Build foo.dll on Windows or foo.so on MacOS/LinuxOS

R CMD SHLIB foo.c

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 “err":

#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
}

How to develop Rpk?

Rstudio
R command-line:

    package.skeleton("myrpk")

How to pack Rpks?

cmd (windows) or terminal (MacOS/Linux):

    R CMD build myrpk
    R CMD Rd2pdf myrpk
    R CMD check myrpk

cran: https://cran.r-project.org/
github: https://github.com/
Bioconductor: https://www.bioconductor.org/
Omegahat: http://www.omegahat.net/

Example: FactSum

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>

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
}

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.

Advantages:

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

An example can be found in PI

The web-based calculator PI is available too
The R package PI is developed by using GMP
Test on Ubuntu
- Intall R on Ubuntu
- add “deb https://cloud.r-project.org/bin/linux/ubuntu disco-cran35/” to “/etc/apt/sources.list”
```
sudo vi /etc/apt/sources.list
```
- Install R
```
sudo apt-get install r-base-dev
```

Rcpp

It is a R package “Rcpp”
Seamless R and C++ Integration: Introduction

How to use Rcpp?

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”))

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')"
```

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
```

Others Calling C/C++

Rpk’s devolopment

extension function

R call C/C++

Examples

How to develop Rpk?

How to pack Rpks?

Where and How to share Rpks?

Example: FactSum

Call foo.c in Rpks

Call executable files in Rpk

Drawbacks:

Advantages:

An example can be found in PI

Rcpp

How to use Rcpp?

An example: tensorMam

Useful functions

Matlab Calling C/C++