R高级编程技巧及Rcpp的介绍.ppt

1、R高级编程技巧及Rcpp的介绍,颜林林 北京大学生物信息中心 2011年5月28日,R语言的特点,FREE：开源、免费、自由、灵活, a nrow(a),截至2011年5月28日，CRAN上总计可用软件包数：,3008,遵守各种不同开源协议：,R语言的特点,FREE：开源、免费、自由、灵活, sort(table(gsub( .*, , a,License), decreasing = T),高校、科研机构、企事业单位 金融、通讯、医药、生物、环境,R语言的广泛运用,COS论坛,创办于2006年5月19日 截至2011年5月28日： 整个论坛： 主题：16,000+；回复：109,000+ “

2、S-Plus 发胸牌给p; 发给p; ,生活中的思维方式： 发资料、胸牌等 给参加R大会的人,能用更少的话说清楚，就别罗嗦！, a a * 2 sin(a),R是一门解释型语言, a # 向量运算 b # 循环 b for (i in seq_along(a) bi - sin(ai) ,R是一门解释型语言, system.time(replicate(10000, b b system.time(replicate(10000, for (i in seq_along(a) bi - sin(ai) user system elapsed 6.15 0.02 6.44,把苦力活交给底层，尽量

3、避免使用循环！,如何掌握好向量运算,R语言基本数据类型 vector, matrix, array, data.frame, list, . 元素的提取方法 a1:3, a-2, aa3, ac(x,y), . 数据类型之间的互相转换 as.matrix(), as.vector(), as.character(), . apply系列函数 apply(), sapply(), lapply(), tapply(), . 其他各种函数 sum(), mean(), cumsum(), combn(), .,R语言元素类型,整数（integer） 实数（numeric, double） 复数（c

4、omplex） 字符（character） 逻辑（logical） 原始数据（raw） 因子（factor）,1:10 c(1.1, 3.14, 10) c(1+i, 3-2i) c(a, b, COS) c(TRUE, FALSE) as.raw(48) as.factor(letter1:3),R语言基本数据类型,向量（vector） 列表（list） 矩阵（matrix） 数组（array） 数据框（data.frame）,R语言基本数据类型,向量（vector） 列表（list） 矩阵（matrix） 数组（array） 数据框（data.frame）,一维，同类型元素 一维，不同类型

5、元素 二维，同类型元素 二维或多维，同类型元素 二维，每列元素同类型,学习R函数时应注意,apply系列函数,函数apply, (m apply(m, 1, sum) # 相当于 rowSums(m) 1 16 20 apply(m, 2, sum) # 相当于 colSums(m) 1 3 7 11 15 apply(m, 1:2, function(x) x / 2) ,1 ,2 ,3 ,4 1, 0.5 1.5 2.5 3.5 2, 1.0 2.0 3.0 4.0,函数lapply和sapply, (l lapply(l, sum) # 返回列表 $a 1 15 $b 1 40 sapp

6、ly(l, sum) # 返回向量 a b 15 40,函数sapply和vapply, (X sapply(X, fivenum) ,1 ,2 ,3 1, 1.0 1.0 1 2, 1.5 1.5 2 3, 2.0 2.5 3 4, 2.5 3.5 4 5, 3.0 4.0 5,函数sapply和vapply, (v vapply(X, fivenum, v) ,1 ,2 ,3 Min. 1.0 1.0 1 1st Qu. 1.5 1.5 2 Median 2.0 2.5 3 3rd Qu. 2.5 3.5 4 Max. 3.0 4.0 5,函数tapply, (d c(class(d),

7、class(d$age), class(d$gender) 1 data.frame numeric factor tapply(d$age, d$gender, mean) Female Male 24.0 24.5,函数mapply, mapply(function(x, y) seq_len(x) + y, c(a = 1, b = 2, c = 3), # names from first c(A = 10, B = 0, C = -10) $a 1 11 $b 1 1 2 $c 1 -9 -8 -7,函数rapply, (X rapply(X, function(x) x) a b.

8、c d 3.14159265358979 1 a test,举例：求DNA互补链,DNA由A、C、T、G组成DNA两条链方向相反成对碱基，A与T配对，C与G配对 已知其中一条为“ACTGAAGTGC”求另一条序列,举例：求DNA互补链,# 方法一：循环 revcom1 - function(DNA) result - for (i in nchar(DNA):1) n - substr(DNA, i, i) if (n = A) result - paste(result, T, sep = ) if (n = T) result - paste(result, A, sep = ) if (

9、n = C) result - paste(result, G, sep = ) if (n = G) result - paste(result, C, sep = ) result ,举例：求DNA互补链,# 方法二：*apply() revcom2 - function(DNA) s - strsplit(DNA, )1 s - sapply(s, function(x) if (x = A) return(T) if (x = T) return(A) if (x = G) return(C) if (x = C) return(G) ) paste(rev(s), collapse

10、= ) ,# 方法三：names() revcom3 - function(DNA) tr - c(A,C,G,T) names(tr) - rev(tr) paste(trrev(strsplit(DNA, )1), collapse = ) , tr T G C A A C G T,举例：求DNA互补链,方法比较, library(rbenchmark) benchmark(revcom1(DNA), revcom2(DNA), revcom3(DNA), columns = c(test, replications, elapsed,relative), order = relative

11、,replications = 1000) test replications elapsed relative 3 revcom3(DNA) 1000 0.163 1.000000 2 revcom2(DNA) 1000 0.597 3.662577 1 revcom1(DNA) 1000 0.624 3.828221,R的灵活性,R是什么语言写成的？ C/C+、Fortran、R R的外部扩展 Writing R Extensions,R与各种语言,R + bash / Rscript RDCOM rJava rpy, rpy2 Embeded Rserve、R API、Rcpp,R 与

12、C+,R（解释型语言） 无需编译 不需要其它文件 弱类型语言 灵活性好 运行速度较慢 其他特性 向量运算 大量统计和绘图函数,C+（编译型语言） 需要事先编译 需要头文件和库文件 强类型语言 灵活性相对较差 运行速度快 其他特性 面向对象 模板与泛型编程,从C到C+,字符串、数组等类型 内存的管理 模板和泛型编程 STL (Standard Template Library),Rcpp相关历史,RQuantLib RcppTemplate Rserve rcppbind / RAbstraction / RObjects / CXXR R API Rcpp / RInside / RProto

13、Buf,Rcpp的运用,在R中调用C+模块（Rcpp） 用C+实现计算功能 通过C+调用其他软件库 在C+中使用R（RInside） 向量运算（STL） 使用R的统计函数 使用R的绘图函数,举例：Rcpp和inline带来的速度,f - function(n, x=1) for (i in 1:n) x=1/(1+x) g - function(n, x=1) for (i in 1:n) x=(1/(1+x) h - function(n, x=1) for (i in 1:n) x=(1+x)(-1) j - function(n, x=1) for (i in 1:n) x=1/1+x

14、k - function(n, x=1) for (i in 1:n) x=1/1+x library(rbenchmark) N - 1e5 benchmark(f(N,1), g(N,1), h(N,1), j(N,1), k(N,1), columns=c(test, replications, elapsed, relative), order=relative, replications=10),举例：Rcpp和inline带来的速度,f - function(n, x=1) for (i in 1:n) x=1/(1+x) g - function(n, x=1) for (i i

15、n 1:n) x=(1/(1+x) h - function(n, x=1) for (i in 1:n) x=(1+x)(-1) j - function(n, x=1) for (i in 1:n) x=1/1+x k - function(n, x=1) for (i in 1:n) x=1/1+x test replications elapsed relative 5 k(N, 1) 10 6.103 1.000000 1 f(N, 1) 10 6.426 1.052925 4 j(N, 1) 10 6.835 1.119941 2 g(N, 1) 10 7.339 1.202523

16、 3 h(N, 1) 10 9.226 1.511716,举例：Rcpp和inline带来的速度,library(inline) src (ns); double x = as(xs); for (int i=0; in; i+) x=1/(1+x); return wrap(x); l - cxxfunction(signature(ns=integer, xs=numeric), body=src, plugin=Rcpp),举例：Rcpp和inline带来的速度, benchmark( f(N,1), g(N,1), h(N,1), j(N,1), k(N,1), l(N,1), col

17、umns=c(test, replications, elapsed, relative), order=relative, eplications=10) test replications elapsed relative 6 l(N, 1) 10 0.027 1.0000 5 k(N, 1) 10 6.190 229.2593 1 f(N, 1) 10 6.379 236.2593 4 j(N, 1) 10 6.808 252.1481 2 g(N, 1) 10 7.267 269.1481 3 h(N, 1) 10 9.184 340.1481,举例：Rcpp API,#include RcppExport SEXP foo(SEXP a, SEXP b) Rcpp:NumericVector xa(a); Rcpp:NumericVector xb(b); int n = std:min(a.size(),