编写Rcpp包03：为第三方库提供R语言接口

这是编写R包的一个全新尝试。在以往开发R包的过程中，往往是在自动生成的R包框架下开发算法。自动生成的 `RcppExports.cpp` 为R提供了调用算法的接口，其他C++文件是对R接口的具体实现。这事实上是一种“先包后库”的方式。这种方式容易缺乏扩展性，例如如果要将同样的功能移植到Python中、Qt程序中，则往往需要重新编写全部代码。这时就需要一种“先库后包”的开发方式，我们先用C++提供库的核心功能（我们称之为核心库），然后再为其编写R接口，相当于用R语言对核心库进行了一次包装，这个用于包装核心库的包不提供任何计算功能，只对核心库的输入输出进行处理，使之可以实现在特定环境下的功能，如结果的控制台输出、绘图等。下面我们以 [**GWmodel3**](https://github.com/GWmodel-Lab/GWmodel3) 为例，展示如何以这种方式为库 [**libgwmodel**](https://github.com/GWmodel-Lab/libgwodel) 编写R包。

# 包结构

总体来看，新方式编写的R包结构和一般的R包结构一样。但是通常在 `src` 文件夹下面将库整体置入。

> 即使包中使用了Armadillo，也无需使用**RcppArmadillo**包创建框架，使用**Rcpp**包即可。

```diff
GWmodel3
├── DESCRIPTION
├── NAMESPACE
├── R
│   └── RcppExports.R
├── data
├── man
└── src
    ├── Makevars
    ├── Makevars.win
    ├── RcppExports.cpp
+   └── libgwmodel

```

之所以整体置入，是因为可以使用 git submodule 功能进行版本管理。当我们切换到 `src` 文件夹下时，输入以下命令

```bash
git submodule add https://github.com/GWmodel-Lab/libgwmodel.git
```

就完成了 **libgwmodel** 库的置入。

> 在R打包过程中，这个文件夹中的所有文件都会被打包。如果有一些文件或文件夹不想被打包，可以将路径放在 `.Rbuildignore` 文件中。

# 开发

库 **libgwmodel** 提供了一个头文件 `gwmodel.h` 以引入库的功能。在 `RcppExports.cpp` 中，可以像一般C++开发的方式直接引入该库文件。注意不要直接引入 `RcppArmadillo.h` 头文件。

```diff
#include <Rcpp.h>
+ #include <armadillo>
+ #include "gwmodel.h"
+ #include "gwmodelpp/GwmLogger.h"
```

比较麻烦的是，如果不引入 `RcppArmadillo.h` ，那么无法直接将R的 `NumericMatrix`  等对象通过 `as()`函数直接转换为 **Armadillo** 的`mat` 对象。解决的方法是可以仿写，例如

```cpp
inline arma::mat myas(const NumericMatrix& rmat)
{
    return arma::mat(rmat.begin(), rmat.nrow(), rmat.ncol());
}

inline arma::vec myas(const NumericVector& rvec)
{
    return arma::vec(rvec.begin(), rvec.size());
}

inline SEXP mywrap(const arma::mat& amat)
{
    RObject x = wrap(amat.begin(), amat.end());
    x.attr("dim") = Dimension(amat.n_rows, amat.n_cols);
    return x;
}

inline SEXP mywrap(const arma::vec& avec)
{
    RObject x = wrap(avec.begin(), avec.end());
    return x;
}
```

然后在开发中，像一般C++开发调用第三方库的方式开发即可。例如如果我们要写 Basic GWR 预测的函数，就可以这样

```cpp
RcppExport SEXP _GWmodel_gwr_basic_predict(
    SEXP pcoordsSEXP, SEXP xSEXP, SEXP ySEXP, SEXP coordsSEXP,
    SEXP bwSEXP, SEXP adaptiveSEXP, SEXP kernelSEXP,
    SEXP longlatSEXP, SEXP pSEXP, SEXP thetaSEXP,
    SEXP interceptSEXP, SEXP parallel_typeSEXP, SEXP parallel_argSEXP)
{
BEGIN_RCPP
    Rcpp::RObject rcpp_result_gen;
    Rcpp::RNGScope rcpp_rngScope_gen;
    Rcpp::traits::input_parameter< const NumericMatrix& >::type pcoords(pcoordsSEXP);
    Rcpp::traits::input_parameter< const NumericMatrix& >::type x(xSEXP);
    Rcpp::traits::input_parameter< const NumericVector& >::type y(ySEXP);
    Rcpp::traits::input_parameter< const NumericMatrix& >::type coords(coordsSEXP);
    Rcpp::traits::input_parameter< double >::type bw(bwSEXP);
    Rcpp::traits::input_parameter< bool >::type adaptive(adaptiveSEXP);
    Rcpp::traits::input_parameter< size_t >::type kernel(kernelSEXP);
    Rcpp::traits::input_parameter< bool >::type longlat(longlatSEXP);
    Rcpp::traits::input_parameter< double >::type p(pSEXP);
    Rcpp::traits::input_parameter< double >::type theta(thetaSEXP);
    Rcpp::traits::input_parameter< bool >::type intercept(interceptSEXP);
    Rcpp::traits::input_parameter< size_t >::type parallel_type(parallel_typeSEXP);
    Rcpp::traits::input_parameter< IntegerVector >::type parallel_arg(parallel_argSEXP);

    // Logger
    GwmLogger::logger = printer;

    // Convert data types
    arma::mat mpcoords = myas(pcoords);
    arma::mat mx = myas(x);
    arma::vec my = myas(y);
    arma::mat mcoords = myas(coords);
    std::vector<int> vpar_args = as< std::vector<int> >(Rcpp::IntegerVector(parallel_arg));

    // Make Spatial Weight
    CGwmBandwidthWeight bandwidth(bw, adaptive, CGwmBandwidthWeight::KernelFunctionType((size_t)kernel));
    CGwmDistance* distance = nullptr;
    if (longlat)
    {
        distance = new CGwmCRSDistance(true);
    }
    else
    {
        if (p == 2.0 && theta == 0.0)
        {
            distance = new CGwmCRSDistance(false);
        }
        else
        {
            distance = new CGwmMinkwoskiDistance(p, theta);
        }
    }
    CGwmSpatialWeight spatial(&bandwidth, distance);
  
    // Make Algorithm Object
    CGwmGWRBasic algorithm(mx, my, mcoords, spatial, false, intercept);
    switch (ParallelType(size_t(parallel_type)))
    {
    case ParallelType::SerialOnly:
        algorithm.setParallelType(ParallelType::SerialOnly);
        break;
    case ParallelType::OpenMP:
        algorithm.setParallelType(ParallelType::OpenMP);
        algorithm.setOmpThreadNum(vpar_args[0]);
    default:
        algorithm.setParallelType(ParallelType::SerialOnly);
        break;
    }

    // Return Results
    mat betas = algorithm.predict(mpcoords);

    rcpp_result_gen = mywrap(betas);
    return rcpp_result_gen;
END_RCPP
}
```

如果不仿写 `as()` 和 `wrap()` 函数，理论上可以直接引入 `RcppArmadilloAs.h` 和 `RcppArmadilloWrap.h` 以导入 `as()` 和 `wrap()` 函数，但是没有实践过。

# 编译

即使置入了库，但该库不能直接进行编译，需要我们编写 `Makevars` 以保证R可以正确带着库进行编译。如果有类似于**GSL**的依赖库，则需要进一步编写 `configure.ac` 文件。

为了让R知道我们库中有哪些文件需要编译，我们需要将这些文件添加到变量 `OBJECTS` 中，这个变量也要包含了 `RcppExports.cpp` 或其他不在库中的 C++ 文件。例如

```makefile
DEF_DOUBLE_EPS = 1e-8
DEF_M_PI = 3.14159265358979323846
DEF_M_PI_2 = 1.57079632679489661923

LIBGWMODEL_CXXFLAGS = -I./libgwmodel/include -I./libgwmodel/include/gwmodelpp \
	-DM_PI=$(DEF_M_PI) -DM_PI_2=$(DEF_M_PI_2) -DDOUBLE_EPS=$(DEF_DOUBLE_EPS) -DDBL_MAX=__DBL_MAX__ \
	-DARMA_NO_DEBUG

PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) $(LIBGWMODEL_CXXFLAGS) $(GSL_CFLAGS) -DRCPP_USE_GLOBAL_ROSTREAM
PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) $(GSL_LIBS)

OBJECTS_LIBGWMODEL = \ 
	libgwmodel/src/gwmodelpp/CGwmBandwidthSelector.o \ 
	libgwmodel/src/gwmodelpp/CGwmGWRBase.o \
	libgwmodel/src/gwmodelpp/CGwmGWRBasic.o \ 
	libgwmodel/src/gwmodelpp/CGwmSpatialAlgorithm.o \
	libgwmodel/src/gwmodelpp/CGwmSpatialMonoscaleAlgorithm.o \ 
	libgwmodel/src/gwmodelpp/CGwmVariableForwardSelector.o \ 
        ...

OBJECTS_GWMODEL = \
	RcppExports.o

OBJECTS = $(OBJECTS_LIBGWMODEL) $(OBJECTS_GWMODEL)
```

这样R就知道需要编译哪些文件，作为R包的一部分。其中变量 `LIBGWMODEL_CXXFLAGS` 为编译库文件提供了必要的编译参数，使R可以正确编译。其中的一个宏定义 `ARMA_NO_DEBUG` 在库使用了 **Armadillo** 依赖时比较重要，可以避免R在检查中探测到 `std::cout` 和 `std::cerr` 的引用从而违反编写R包的规则。

这种方法的好处是比较简单直接，不容易出问题，在CRAN的编译环境下也可以顺利编译。