Interfacing one C++ codebase with both Python and R

README.md

Interfacing one C++ codebase with both Python and R

Python

The input numpy array is mapped to the input of a C++ wrapper function using pybind11. Then, core function myfunc with templated input type is executed on the mapped matrix.

To execute the test function and profile memory usage, run

/usr/bin/time -v python test.py

R

The input R matrix is mapped to the input of a C++ wrapper function using RcppEigen. Then, core function myfunc with templated input type is executed on the mapped matrix

To execute the test function and profile memory usage, run

/usr/bin/time -v Rscript test.R

myfunc template code

myfunc is written in template code so that the instantiated type can change depending on how it is compiled. Due to differences in Python vs. R storage order and interface package, the input type cannot be constrained further.

core.hpp

// #include <iostream>
//#include <math.h>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/Sparse>

using namespace std;

template <typename Derived>
void myfunc(Eigen::MatrixBase<Derived> &b, double multiplier)
{
    b *= multiplier;
}

template <typename Derived, typename SparseDerived>
void myfunc_complex(
    Eigen::MatrixBase<Derived> &D, 
    Eigen::SparseMatrixBase<SparseDerived> &Dsp, 
    double multiplier)
{
    D *= multiplier;
    Dsp *= multiplier;
    D += Dsp;
}

double sgn(double val) {
  return (double(0) < val) - (val < double(0));
}

double sthresh(double x, double t ){
  return sgn(x) * max(abs(x)-t, 0.0);
}

void sthreshmat(Eigen::MatrixXd & x,
                double tau,
                const Eigen::MatrixXd & t){

  Eigen::MatrixXd tmp1(x.cols(), x.cols());
  Eigen::MatrixXd tmp2(x.cols(), x.cols());

  tmp1 = x.array().unaryExpr(ptr_fun(sgn));
  tmp2 = (x.cwiseAbs() - tau*t).cwiseMax(0.0);

  x = tmp1.cwiseProduct(tmp2);

  return;
}

inline double shrink(double a, double b) {
  if (b < fabs(a)) {
    if (a > 0) return(a-b);
    else       return(a+b);
  } else {
    return(0.0);
  }
}

// template <typename Derived, typename SparseDerived>
// void myfunc_complex(
//     Eigen::MatrixBase<Derived> &D, 
//     Eigen::SparseMatrixBase<SparseDerived> &Dsp, 

// ista algorithm
template <typename DerivedS, typename SparseDerived, typename DerivedL>
void ccista_core(
    Eigen::MatrixBase<DerivedS> &S, 
    Eigen::Map<Eigen::SparseMatrix<SparseDerived> > &X, 
    Eigen::MatrixBase<DerivedL> &LambdaMat, 
    double lambda2,
    double epstol,
    int maxitr,
    int steptype)
{
    Rcpp::Rcout << "in ccista_core:\n" << S << std::endl;
    return;
    // int p = S.cols();
    // Eigen::DiagonalMatrix<double, Eigen::Dynamic> XdiagM(p);
    // Eigen::SparseMatrix<double, Eigen::ColMajor> Xn;
    // Eigen::SparseMatrix<double, Eigen::ColMajor> Step;
    // Eigen::MatrixXd W = S * X;
    // Eigen::MatrixXd Wn(p, p);
    // Eigen::MatrixXd G(p, p);
    // Eigen::MatrixXd Gn(p, p);
    // Eigen::MatrixXd subg(p, p);
    // Eigen::MatrixXd tmp(p, p);
//    
//    // no member function X.diagonal()
//    double h = - X.diagonal().array().log().sum() + 0.5 * ((X*W).diagonal().sum()) + (lambda2 * pow(X.norm(), 2));
//    double hn = 0;
//    double Qn = 0;
//    double subgnorm, Xnnorm;
//    double tau;
//    double taun = 1.0;
//    double c = 0.5;
//    int itr = 0;
//    int loop = 1;
//    int diagitr = 0;
//    int backitr = 0;
//
//    XdiagM.diagonal() = - X.diagonal();
//    G = XdiagM.inverse();
//    G += 0.5 * (W + W.transpose());
//    if (lambda2 > 0) { G += lambda2 * 2.0 * X; }
//
//    while (loop != 0){
//  
//        tau = taun;
//        diagitr = 0;
//        backitr = 0;
//    
//        while ( 1 ) { // back-tracking line search
//    
//            if (diagitr != 0 || backitr != 0) { tau = tau * c; }
//      
//            tmp = Eigen::MatrixXd(X) - tau*G;
//            sthreshmat(tmp, tau, LambdaMat);
//            Xn = tmp.sparseView();
//      
//            if (Xn.diagonal().minCoeff() < 1e-8 && diagitr < 10) { diagitr += 1; continue; }
//      
//            Step = Xn - X;
//            Wn = S * Xn;
//            Qn = h + Step.cwiseProduct(G).sum() + (1/(2*tau))*pow(Step.norm(),2);
//            hn = - Xn.diagonal().cwiseAbs().array().log().sum() + 0.5 * (Xn.cwiseProduct(Wn).sum());
//      
//            if (lambda2 > 0) { hn += lambda2 * pow(Xn.norm(), 2); }
//      
//            if (hn > Qn) { backitr += 1; } else { break; }
//    
//        }
//    
//        XdiagM.diagonal() = - Xn.diagonal();
//        Gn = XdiagM.inverse();
//        Gn += 0.5 * (Wn + Wn.transpose());
//    
//        if (lambda2 > 0) { Gn += lambda2 * 2 * Eigen::MatrixXd(Xn); }
//    
//        if ( steptype == 0 ) {
//            taun = ( Step * Step ).eval().diagonal().array().sum() / (Step.cwiseProduct( Gn - G ).sum());
//        } else if ( steptype == 1 ) {
//            taun = 1;
//        } else if ( steptype == 2 ){
//            taun = tau;
//        }
//    
//        tmp = Eigen::MatrixXd(Xn).array().unaryExpr(ptr_fun(sgn));
//        tmp = Gn + tmp.cwiseProduct(LambdaMat);
//        subg = Gn;
//        sthreshmat(subg, 1.0, LambdaMat);
//        subg = (Eigen::MatrixXd(Xn).array() != 0).select(tmp, subg);
//        subgnorm = subg.norm();
//        Xnnorm = Xn.norm();
//    
//        //maxd[itr] = subgnorm/Xnnorm;
//        X = Xn;
//        h = hn;
//        G = Gn;
//    
//        itr += 1;
//    
//        loop = int((itr < maxitr) && (subgnorm/Xnnorm > epstol));
//  
//    }

  //nitr = itr;

}

test.py

import numpy as np
import cppimport

# cppimport.set_quiet(False)
# cppimport.force_rebuild()
gconcord = cppimport.imp("wrappy")

A = np.ones((5, 10))
m = 3.2

gconcord.myfunc_wrapper(A, m)
print(A[0,0], 'should be equal to', m)

test.R

library(RcppEigen);

Rcpp::sourceCpp('wrapr.cpp')

n = 5
p = 10
m = 3.2

A = matrix(1, n, p)
myfunc_wrapper(A, m)

cat(A[1,1], 'should be equal to', m, '\n')

wrappy.cpp

<%
cfg['compiler_args'] = ['-std=c++11']
cfg['include_dirs'] = ['eigen-3.3.7']
setup_pybind11(cfg)
%>
#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
#include "core.hpp"

namespace py = pybind11;

void myfunc_wrapper(py::EigenDRef<Eigen::MatrixXd> K, double m) {
    myfunc(K, m);
}

PYBIND11_MODULE(wrappy, m) {
    m.def("myfunc_wrapper", &myfunc_wrapper);
}

wrapr.cpp

#include <RcppEigen.h>
#include <iostream>
// #include <Eigen/Core>
// #include <Eigen/Dense>
#include <Eigen/Sparse>
#include "core.hpp"

// [[Rcpp::depends(RcppEigen)]]

// [[Rcpp::export]]
void myfunc_wrapper(Eigen::Map<Eigen::MatrixXd> M, double multiplier){
    myfunc(M, multiplier);
}

// [[Rcpp::export]]
void sparse(
    Eigen::Map<Eigen::MatrixXd> M, 
    Eigen::Map<Eigen::SparseMatrix<double> > Msp,
    double multiplier)
{
    Rcpp::Rcout << "M:\n" << M << std::endl;
    Rcpp::Rcout << "Msp:\n" << Msp << std::endl;
    sparse_core(M, Msp, multiplier);
}

using namespace Rcpp;

namespace impl {

template <int RTYPE>
Vector<RTYPE> ends(const Vector<RTYPE>& x, int n)
{
    n = std::min((R_xlen_t)n, x.size() / 2);
    Vector<RTYPE> res(2 * n);

    std::copy(x.begin(), x.begin() + n, res.begin());
    std::copy(x.end() - n, x.end(), res.begin() + n);

    return res;
}

} // impl

// SEXP ends(SEXP x, int n = 6) {
// [[Rcpp::export]]
SEXP ends(
    Eigen::Map<Eigen::MatrixXd> S, 
    SEXP xx, int n = 6) {
// SEXP ends(Rcpp::NumericMatrix xx, int n = 6) {

    Eigen::Map<Eigen::MatrixXd> x(Rcpp::as<Eigen::Map<Eigen::MatrixXd> >(xx));
    
    // Eigen::Map<Eigen::VectorXd> XS(Rcpp::as<Eigen::Map<Eigen::VectorXd> >(X)); 
    // return true;
    
    // NumericMatrix xx(x);
    Rcout << x << std::endl;
    Rcout << x.rows() << x.cols() << std::endl;
    
    Eigen::MatrixXd asdf = Eigen::MatrixXd::Constant(3, 3, (double)n);
    Rcout << asdf << std::endl;
    
    switch (TYPEOF(xx)) {
        case INTSXP: {
            Rcout << "INTSXP" << std::endl;
            return impl::ends(as<IntegerVector>(xx), n);
        }
        case REALSXP: {
            Rcout << "REALSXP" << std::endl;
            // Eigen::MatrixXd xx = 
            return impl::ends(as<NumericVector>(xx), n);
        }
        case STRSXP: {
            Rcout << "STRSXP" << std::endl;
            return impl::ends(as<CharacterVector>(xx), n);
        }
        case LGLSXP: {
            Rcout << "LGLSXP" << std::endl;
            return impl::ends(as<LogicalVector>(xx), n);
        }
        case CPLXSXP: {
            Rcout << "CPLXSXP" << std::endl;
            return impl::ends(as<ComplexVector>(xx), n);
        }
        default: {
            warning(
                "Invalid SEXPTYPE %d (%s).\n",
                TYPEOF(xx), type2name(xx)
            );
            return R_NilValue;
        }
    }
}