41 Rcpp糖
在C++中,向量和矩阵的运算通常需要逐个元素进行, 或者调用相应的函数。 Rcpp通过C++的表达式模板(expression template)功能, 可以在C++中写出像R中对向量和矩阵运算那样的表达式。 这称为Rcpp糖(sugar)。
R中的很多函数如sin
等是向量化的,
Rcpp糖也提供了这样的功能。
Rcpp糖提供了一些向量化的函数如ifelse
, sapply
等。
比如,两个向量相加可以直接写成x + y
而不是用循环或迭代器(iterator)逐元素计算;
若x
是一个NumericVector,
用sin(x)
可以返回由x
每个元素的正弦值组成的NumericVector。
Rcpp糖不仅简化了程序, 还提高了运行效率。
41.1 简单示例
比如,函数 \[\begin{aligned} f(x,y) = \begin{cases} x^2 & x < y, \\ -y^2 & x \geq y \end{cases} \end{aligned}\]
如下的程序可以在C++中定义一个向量化的版本:
sourceCpp(code='
#include <Rcpp.h>
using namespace Rcpp;
//[[Rcpp::export]]
NumericVector f11(NumericVector x, NumericVector y){
return ifelse(x < y, x*x, -(y*y));
}
')
f11(c(1, 3), c(4,2))
## [1] 1 -4
上面简单例子中,x < y
是向量化的,
x * x
, y * y
, -(y * y)
都是向量化的。
ifelse
也是向量化的。
41.2 向量化的运算符
41.2.1 向量化的四则运算
Rcpp糖向量化了+, -, *, /
。
设x, y
都是NumericVector, 则
x + y
, x - y
, x * y
, x * y
将返回对应元素进行四则运算后的NumericVector变量。
向量与标量运算,
如 x + 2.0
, 2.0 - x
, x * 2.0
, 2.0 / x
将返回标量与向量每个元素进行四则运算后的NumericVector变量。
还可以进行混合四则运算,如:
NumericVector res = x * y + y / 2.0;
NumericVector res = x * (y - 2.0);
NumericVector res = x / (y * y);
参加四则运算的或者都是同基本类型的向量而且长度相等, 或者一边是向量,一边是同类型的标量。
注意:对向量整体赋一个相同的值, 不能简单地写成如x=0;
这样的赋值,
需要用循环或STL的fill算法, 如std::fill(x.begin(), x.end(), 0);
。
41.2.2 向量化的二元逻辑运算
Rcpp糖扩展了两个元素的比较到两个等长向量的比较,
以及一个向量与一个同类型标量的比较,
结果是一个同长度的LogicalVector。
比较运算符包括<, >, <=, >=, ==, !=
。
也可以使用嵌套的表达式, 比如,设x, y是两个NumericVector,可以写:
41.3 用Rcpp访问数学函数
在C和C++代码中可以调用R中的函数, 但是格式比较复杂, 而且原来不支持向量化。 Rcpp糖则使得从C++中调用R函数变得和在R调用函数格式类似。
R源程序中提供了许多数学和统计相关的函数, 这些函数可以编译成一个独立的函数库, 供其它程序链接使用,函数内容在R的Rmath.h头文件中有详细列表。
Rcpp糖把R中的数学函数在C++中向量化了, 输入一个向量, 结果是对应元素为函数值的向量。 自变量类型可以是数值型或整数型。 这些数学函数包括abs, exp, floor, ceil, pow等。 在Rcpp中支持的C++源程序中调用这些函数, 最简单的一种方法是直接在Rcpp 名字空间中使用和R中相同的函数名和调用方法, 类似sqrt这样的函数允许输入一个向量, 对向量的每个元素计算。
比如,下面的程序输入一个向量,计算相应的标准正态分布函数值:
sourceCpp(code='
#include <Rcpp.h>
using namespace Rcpp;
//[[Rcpp::export]]
NumericVector f10(NumericVector x){
NumericVector y(x.size());
y = pnorm(x);
return y;
}
')
f10(c(0, 1.96, 2.58))
## [1] 0.5000000 0.9750021 0.9950600
如果需要对单个的数值计算,
可以使用Rmath.h中定义的带有Rf_
的版本,
如:
这里用::
使用了缺省的名字空间。
注意所有自变量都不能省略。
Rcpp还提供了一个R
名字空间,
可以用不带Rf_
前缀的函数,
但是自变量也不能省略。如
在Rcpp的R
名字空间中有许多的数学和统计相关的函数,
各函数的自变量参见Rcpp的Rmath.h文件。
R中提供了许多分布的分布密度(或概率质量函数)、分布函数、分位数函数, 分布密度函数和概率质量函数命名类似dxxxx, 分布函数命名类似pxxxx, 分位数函数命名类似qxxxx。 在Rcpp糖支持下, 这些函数可以直接用类似R中的格式调用。 如
sourceCpp(code='
#include <Rcpp.h>
using namespace Rcpp;
//[[Rcpp::export]]
List f14(){
NumericVector x =
NumericVector::create(0, 1.96, 2.58);
NumericVector p =
NumericVector::create(0.95, 0.975, 0.995);
NumericVector y1 = dnorm(x, 0.0, 1.0);
NumericVector y2 = pnorm(x, 0.0, 1.0);
NumericVector y3 = qnorm(p, 0.0, 1.0);
return List::create(Named("y1")=y1,
Named("y2")=y2, Named("y3")=y3);
}
')
f14()
## $y1
## [1] 0.39894228 0.05844094 0.01430511
##
## $y2
## [1] 0.5000000 0.9750021 0.9950600
##
## $y3
## [1] 1.644854 1.959964 2.575829
##
R中的rxxxx类的函数可以产生各种分布的随机数向量, 随机数向量与当前种子有关。 Rcpp属性会自动地维护随机数发生器的状态使其与R同步。 如
41.4 返回单一逻辑值的函数
在R中, any()
和all()
对一个逻辑向量分别判断是否有任何真值,
以及所有元素为真值。
Rcpp糖在C++中也提供了这样的any()
和all()
函数。
如
sourceCpp(code='
#include <Rcpp.h>
using namespace Rcpp;
//[[Rcpp::export]]
List f12(){
IntegerVector x = seq_len(1000);
LogicalVector res1 = any( x*x < 3 );
LogicalVector res2 = all( x*x < 3 );
return List::create(Named("any")=res1,
Named("all")=res2);
}
')
f12()
## $any
## [1] TRUE
##
## $all
## [1] FALSE
any()
和all()
不是直接返回逻辑值,
而是返回一个类对象,
该类定义了is_true
, is_false
, is_na
方法与向SEXP转换的运算符。
此种结果可以保存到LogicalVector中, 但不能赋值到bool类型,
因为可能有缺失值。
逻辑型糖表达式结果可以用 is_true
, is_false
, is_na
来判断结果。
如
在求any()
结果时,一旦遇到一个真值结果就为真值,
即使后面有缺失值也没有关系。
在求all()
结果时,一旦遇到一个假值结果就为假值,
即使后面有缺失值也没有关系。
41.5 返回糖表达式的函数
is_na
以任何糖表达式为输入,
输出一个逻辑类型的元素个数相同的糖表达式。
结果每个元素当输入中对应元素缺失时为TRUE, 否则为FALSE。
如
IntegerVector x = IntegerVector::create(0, 1, NA_INTEGER, 3);
LogicalVector res1 is_na(x);
LogicalVector res2 = all( is_na(x) );
if( is_true( any( ! is_na(x) ) ) ) ...
seq_along
输入一个向量,
输出一个元素为该向量的各个下标值的糖表达式。
如
注意上述程序不会计算x*x*x*x*x
的值而只利用其结果的元素个数。
所以两次调用的计算量是一样的。
seq_len
自变量为个数,
返回一个元素值为正数的元素个数等于自变量值的糖表达式,
第\(i\)个元素等于\(i\)。
常可与sapply
, lapply
配合使用。
如
sourceCpp(code='
#include <Rcpp.h>
using namespace Rcpp;
//[[Rcpp::export]]
List f13(){
IntegerVector x =seq_len(3);
List y = lapply( seq_len(3), seq_len);
return List::create(Named("x")=x,
Named("y")=y);
}
')
f13()
## $x
## [1] 1 2 3
##
## $y
## $y[[1]]
## [1] 1
##
## $y[[2]]
## [1] 1 2
##
## $y[[3]]
## [1] 1 2 3
##
pmin
和pmax
自变量为两个等长的向量或糖表达式,
返回长度相同的结果,结果元素只等于对应元素的最小值或最大值。
自变量也可以取一个向量一个标量, 向量的每个元素与标量比较得到最小值或最大值。 如
ifelse
函数有三个自变量,
第一自变量是一个逻辑型向量值的糖表达式,
第二和第三自变量或者是两个同类型并与第一自变量等长的糖表达式,
或者其中一个是同类型标量,
结果仍为向量型糖表达式,
第\(i\)元素当第一自变量第\(i\)元素为真时取第一自变量的第\(i\)元素,
当第一自变量第\(i\)元素为假时取第二自变量的第\(i\)元素,
当第一自变量第\(i\)元素为缺失值时取缺失值。
如
IntegerVector x, y;
IntegerVector res1 = ifelse( x < y, x, (x+y)*y );
IntegerVector res2 = ifelse( x > y, x, 2 );
sapply
函数第一自变量是一个向量或列表, 第二自变量是一个函数。
返回值类型在编译时从函数的结果类型导出。
第二自变量可以是任意的C++函数,
比如, 可以是如下的重载的模板化函数:
下面的例子中sapply
第二自变量使用了functor,
这是一种能够产生函数的函数:
template <typename T>
struct square : std::unaray_function<T,T> {
T operator() (const T& x) {
return x * x;
}
}
sapply( seq_len(4), square<int>() );
lapply
函数与sapply
函数基本相同,
只不过lapply
函数总是返回列表,
列表在Rcpp中为List或GenericVector, 在R API中类型为VECSXP。
sign
函数输入一个数值型或整型表达式,
返回各元素值在\(\{1, 0, -1, \text{NA} \}\)中取值的糖表达式,
可以保存到IntegerVector中。
结果各元素的取值表示输入中对应元素为正、零、负和缺失。
如
diff
函数自变量为数值型或整数型向量或有这样结果的糖表达式,
输出后一元素减去前一元素的结果, 结果长度比输入长度少一。
如
41.6 R与Rcpp不同语法示例
RCpp糖通过一些现代的C++技术, 支持部分的向量化运算, 比如,NumericVector与标量相加, 两个等长NumericVector相加, 对NumericVector计算如绝对值这样的数学函数值等。 但是,C++毕竟和R有很大差别, 要区分Rcpp能做的和不能做的。
例如,NumericVector为了把向量所有元素都改成同一值,
不能直接用等于赋值。
可以用std::fill(y.begin(), y.end(), 99.99)
这样的做法。