发个学习贴，新手可以进来看看一些问题的理解

superdesolator

回复第119楼的 Ihavenothing：追上你们还早呢[s:11][s:18]

superdesolator

Try 能不能继续发帖，回了站长帖子外面有显示回复，里面没看到

PS：当时找了半天没有显示第7页... [s:12]

superdesolator

82.ADR之Performance code 之 Peformance

我相信接下来的这章的几个小节都会非常有趣，主要是作者介绍关于R的表现方面的知识，包括垃圾回收，包括Rcpp,包括R的C接口这些我们迫不及待想知道的知识[s:11]

这节作者说的知识都是为后面几个小节作预备知识的，但本身营养也蛮丰富的

R语言本身和R语言的一个实现GNU-R（就是目前大家都在用的R）是不同的，这就像C++和GCC C++是不同的，一个是语言本身，一个是语言的一种实现，对于C++这个庞大复杂的语言来说，其实貌似没有什么实现能准确刻画出C++本身，以致于我在学C++的时候总是发现我的GCC（WINDOWS版的）出来的结果和C++ PRIMER或者THINKING IN C++上说的不一样，困惑许久才知道，是实现的问题，对于R来说，也一样，R语言本身的定义比C++要不正式多了，所以作者说R语言本身貌似就存在于GNU-R是如何工作的...

在深入介绍之前，作者介绍了microbenchmarking这个包

library(microbenchmark)</p>
<p>x <- runif(100)<br />
microbenchmark(<br />
  sqrt(x),<br />
  x ^ 0.5<br />
,times=100,unit='ns'<br />
)<br />
#> Unit: nanoseconds<br />
#>     expr    min     lq median     uq    max neval<br />
#>  sqrt(x)  1,600  1,760  1,860  1,960 13,600   100<br />
#>    x^0.5 15,300 15,500 15,600 15,800 62,200   100

这个包主要用来看语言的表现的，比起system.time()要好的多,因为system.time()不是很准确，因此要跑很多次

n <- 1:1e6<br />
system.time(for (i in n) sqrt(x)) / length(n)<br />
system.time(for (i in n) x ^ 0.5) / length(n)

而我们这里的microbenchmark默认是times=100,我们当然也可以自己修改，我们也还可以unit="ns"来替换，应该是默认是ms,我们可以用的有很多例如"eps"看一秒内运行多少次之类的unit="eps"，最后要注意的是，microbenckmark这个函数在计算每一个表达式一次的时候，本身要花掉800ns,所以我们必须从结果中扣除，不方便是吧？我也觉得，为什么要让读者考虑这个事情呢？应该有深层次的原因吧

所以，实际上一次sqrt(长100的向量)的时间是1600ns-800ns=800ns,0.8us.运行1百万次才0.8s,显然这对我们的影响不是很大（我不知道我有没有理解对这个算式啊，作者下面说的又好像不一样，因为有些结果小于800ns，但是花的时间相对大小还是可以看出来的)

简单熟悉了microbenchmark(mbm简写吧)之后，我们就要从几个方便来介绍R为什么慢这个大话题的框架了：

首先从，语言本身的角度来看：

作者从3个角度看

第一，Extreme dynamism

动态的不好的地方是，我们无法预测函数会怎么发生，对于C++这样编译型的语言，输出的类型是什么，输入的类型是什么基本都是知道的，而R则不是：

x <- 0L<br />
for (i in 1:1e6) {<br />
  x <- x + 1<br />
}

我们当然知道x是integer而R则不知道x是什么类型！所以它得花时间去找+的正确method,这就花了时间了，这是语言本身所决定的.另外的例子：

f <- function(x) NULL</p>
<p>s3 <- function(x) UseMethod("s3")<br />
s3.integer <- f</p>
<p>A <- setClass("A", representation(a = "list"))<br />
setGeneric("s4", function(x) standardGeneric("s4"))<br />
setMethod(s4, "A", f)</p>
<p>B <- setRefClass("B", methods = list(rc = f))</p>
<p>a <- A()<br />
b <- B$new()<br />
microbenchmark(<br />
  fun = f(),<br />
  S3 = s3(1L),<br />
  S4 = s4(a),<br />
  RC = b$rc()<br />
)<br />
#> Unit: nanoseconds<br />
#>  expr    min     lq median     uq     max neval<br />
#>   fun    338    535    632    706  10,900   100<br />
#>    S3  4,660  5,550  6,120  6,810  48,400   100<br />
#>    S4 24,100 26,000 27,500 29,200  78,700   100<br />
#>    RC 25,800 27,700 28,900 32,300 915,000   100

这里我们的结果可能是microseconds，不过不重要，这在于mbm函数的unit默认参数是什么

结果主要说明了S3,S4的method dispatch需要花时间所以比较expensive

第二，Name lookup with mutable environments

a <- 1<br />
f <- function() {<br />
  g <- function() {<br />
    print(a)<br />
    assign("a", 2, envir = parent.frame())<br />
    print(a)<br />
    a <- 3<br />
    print(a)<br />
  }<br />
  g()<br />
}<br />
f()<br />
#> [1] 1<br />
#> [1] 2<br />
#> [1] 3

作者举了这个例子，如果你无法知道弄清楚这个例子，赶紧去回头补补环境的知识

作者用这个例子证明了，R每次都要从头开始找名字，也就是第一次找到a了之后，下次遇到a仍然从头找，不是只找一次。

然后,由于the fact that almost every operation is a lexically scoped function call,所以下面的f中包含了+,-,{,(4个函数调用！并且由于它们是在base env里定义的，所以我们要从f()的定义globalenv()中穿越整个search path中的环境去base environment中找到！！！这个问题我猜大家几乎以前是不知道的，同样我也被震惊了，这么简单的+,-,(,{，几乎每个函数都要用的符号要穿越整个search path( 你可以在R中输入search()来查看这个search path，这个知识在前面的环境贴中也介绍了)，就是这么麻烦啊...

f <- function(x, y) {<br />
  (x + y) ^ 2<br />
}<br />
random_env <- function(parent = globalenv()) {<br />
  letter_list <- setNames(as.list(runif(26)), LETTERS)<br />
  list2env(letter_list, envir = new.env(parent = parent))<br />
}<br />
set_env <- function(f, e) {<br />
  environment(f) <- e<br />
  f<br />
}<br />
f2 <- set_env(f, random_env())<br />
f3 <- set_env(f, random_env(environment(f2)))<br />
f4 <- set_env(f, random_env(environment(f3)))</p>
<p>microbenchmark(<br />
  f(1, 2),<br />
  f2(1, 2),<br />
  f3(1, 2),<br />
  f4(1, 2),<br />
  times = 10000<br />
)<br />
#> Unit: nanoseconds<br />
#>      expr   min    lq median    uq       max neval<br />
#>   f(1, 2) 1,010 1,210  1,370 1,650 1,180,000 10000<br />
#>  f2(1, 2) 1,080 1,270  1,440 1,730 1,040,000 10000<br />
#>  f3(1, 2) 1,160 1,370  1,540 1,830    25,800 10000<br />
#>  f4(1, 2) 1,240 1,430  1,610 1,900    61,700 10000

然后作者就实验了一下加一个环境需要增加多少时间，作者不断更改f()的ee（由environment（））的环境，使得它与base环境之间加入更多的环境，结果时间越来越多

这当然也是R语言本身定义的问题，作者建议说用缓存的方法，但又指出这种方法很难实现，因为R里面的修改一个对象是那么的容易，要确保缓存能及时更新不出错很难，又建议加入更多的常量的东西，让R准确的知道+,-,{,(是什么意思，从而不必去找它们的定义，但作者又指出这种方法会使语言不灵活，平衡很重要

第三，Lazy evaluation overhead

我们知道在R中，函数的参数是惰性求值的，为了实现这个惰性求值，R uses a promise object that contains the expression needed to compute the result and the environment in which to perform the computation.

创建这些对象是要开销的，所以参数越多，开销越大

f0 <- function() NULL<br />
f1 <- function(a = 1) NULL<br />
f2 <- function(a = 1, b = 1) NULL<br />
f3 <- function(a = 1, b = 2, c = 3) NULL<br />
f4 <- function(a = 1, b = 2, c = 4, d = 4) NULL<br />
f5 <- function(a = 1, b = 2, c = 4, d = 4, e = 5) NULL<br />
microbenchmark(f0(), f1(), f2(), f3(), f4(), f5(), times = 10000)<br />
#> Unit: nanoseconds<br />
#>  expr min  lq median  uq       max neval<br />
#>  f0() 187 215    247 276     9,070 10000<br />
#>  f1() 237 281    318 356    11,900 10000<br />
#>  f2() 270 319    362 411    23,500 10000<br />
#>  f3() 319 383    428 489    10,600 10000<br />
#>  f4() 363 438    492 562    77,300 10000<br />
#>  f5() 407 494    552 631 1,350,000 10000

而对于大多数其他一些语言，增加参数的开销很少，对于编译型的，参数没用会被警告甚至被自动去除

我们以一个作者的练习来结束语言层面的慢的原因的讨论：

Why is the cost of name lookup less for functions in the base package?

刚开始我不理解，作者不是说了从golbalenv()穿到base env要经过整个search path,然后这里又这样问，我最后才发现，作者问的是那些定义在base包里的函数自己在找值的时候为什么会花的时候少，理解到这里，那还用说吗？自给自足呗，它只在base包本身以及其它很少的包查吧

OK,介绍完R语言本身定义的问题，让我们来看看GNU-R实现的问题，作者顺便吐槽了R core的保守不作为，不愿意修改R的底层

R语言本身的定义就本身没有达到方法最优化的表现，而GNU-R的实现则离的更远

作者上来就给了下面的例子

microbenchmark(<br />
  mtcars[32, 11],<br />
  mtcars$carb[32],<br />
  mtcars[[c(11, 32)]],<br />
  mtcars[[11]][32],<br />
  .subset2(mtcars, 11)[32]<br />
)<br />
#> Unit: nanoseconds<br />
#>                      expr    min     lq median     uq     max neval<br />
#>            mtcars[32, 11] 30,600 31,500 32,100 32,500 109,000   100<br />
#>           mtcars$carb[32] 15,900 16,800 17,200 17,800 471,000   100<br />
#>       mtcars[[c(11, 32)]] 12,900 13,400 13,800 14,300  20,500   100<br />
#>          mtcars[[11]][32] 12,100 13,000 13,400 14,000  20,300   100<br />
#>  .subset2(mtcars, 11)[32]    477    760    867    929  17,500   100

让我觉得感叹的是，即使作者不是专门介绍data.frame,我也从中知道了[[c(1,2)]]和.subset2()的取法，并且这两种取法比最常用的要快，特别是最后一种取法（我以后可能就用它啦[s:11]）

然后作者就比较了ifesle,pmin,pmax在使用中的速度，问题的背景是让一个向量中的值介于a,b之间

squish_ife <- function(x, a, b) {<br />
  ifelse(x <= a, a, ifelse(x >= b, b, x))<br />
}<br />
squish_p <- function(x, a, b) {<br />
  pmax(pmin(x, b), a)<br />
}<br />
squish_in_place <- function(x, a, b) {<br />
  x[x <= a] <- a<br />
  x[x >= b] <- b<br />
  x<br />
}</p>
<p>x <- runif(100, -1.5, 1.5)<br />
microbenchmark(<br />
  squish_ife(x, -1, 1),<br />
  squish_p(x, -1, 1),<br />
  squish_in_place(x, -1, 1)<br />
)<br />
#> Unit: nanoseconds<br />
#>                       expr    min     lq median     uq     max neval<br />
#>       squish_ife(x, -1, 1) 70,500 82,300 90,400 94,700 119,000   100<br />
#>         squish_p(x, -1, 1) 29,200 32,600 34,300 36,600 652,000   100<br />
#>  squish_in_place(x, -1, 1) 10,000 11,200 12,700 14,100  36,100   100<br />

作者解释说ifelse是有名的慢，它会计算它所有的参数，而pmin,pmax看上去很独特应该会快，其实也很慢，因为它可以取任何数量的参数，然后内部要决定用哪个方法

最牛的方法是用C++

#include <Rcpp.h><br />
using namespace Rcpp;</p>
<p>// [[Rcpp::export]]<br />
NumericVector squish_cpp(NumericVector x, double a, double b) {<br />
  int n = x.length();<br />
  NumericVector out(n);</p>
<p>  for (int i = 0; i < n; ++i) {<br />
    double xi = x[i];<br />
    if (xi < a) {<br />
      out[i] = a;<br />
    } else if (xi > b) {<br />
      out[i] = b;<br />
    } else {<br />
      out[i] = xi;<br />
    }<br />
  }</p>
<p>  return out;<br />
}

这里作者提前给出了这个Rcpp,后面的一节专门详细描述

microbenchmark(<br />
  squish_in_place(x, -1, 1),<br />
  squish_cpp(x, -1, 1)<br />
)<br />
#> Unit: nanoseconds<br />
#>                       expr    min     lq median     uq    max neval<br />
#>  squish_in_place(x, -1, 1) 10,600 11,200 11,600 12,000 45,700   100<br />
#>       squish_cpp(x, -1, 1)  4,900  5,400  5,640  5,930 52,600   100

相对比最好的纯R的实现，C++更快

最后作者介绍了一些其它R实现，我就不列出了，值得一提的是改变已有的计算方式可能会有效提高速度

x <- runif(1e6)<br />
y <- runif(1e6)<br />
z <- sample(c(T, F), 1e6, rep = TRUE)</p>
<p>sum((x + y)[z])

由于是R是向量化运算，所以当向量很大的时候，经常读写内存就会减慢速度，x+y,和z都是临时的大向量，它们要被读写内存，减慢速度

#include <Rcpp.h><br />
using namespace Rcpp;</p>
<p>// [[Rcpp::export]]<br />
double cond_sum_cpp(NumericVector x, NumericVector y, LogicalVector z) {<br />
  double sum = 0;<br />
  int n = x.length();</p>
<p>  for(int i = 0; i < n; i++) {<br />
    if (!z[i]) continue;<br />
    sum += x[i] + y[i];<br />
  }</p>
<p>  return sum;<br />
}

作者说如果能变成上面的C++方式会快8倍，因为只一个中间变量sum就可以了,具体效果如下

cond_sum_r <- function(x, y, z) {<br />
  sum((x + y)[z])<br />
}</p>
<p>microbenchmark(<br />
  cond_sum_cpp(x, y, z),<br />
  cond_sum_r(x, y, z),<br />
  unit = "ms"<br />
)<br />
#> Unit: milliseconds<br />
#>                   expr   min    lq median    uq   max neval<br />
#>  cond_sum_cpp(x, y, z)  7.08  7.33    7.9  8.11  8.47   100<br />
#>    cond_sum_r(x, y, z) 27.90 28.70   29.8 30.90 82.40   100

OK，结束！

superdesolator

83.ADR之Performance code 之 Profiling and benchmarking

首先装devtools::install_github("hadley/lineprof")以及shiny包

1.Measuring performance

我们首先要有个profiler,大概意思就是分析器,分析我们代码中各个片段运行所需要的时间，然后找到瓶颈bottlenecks，也就是花时间花的多的.

为了达到各个代码片段计算时间的功能，作者用了自己写的lineprof这个包，并且也同时指出有其他包如：

summaryRprof(), the proftools package and the profr package以及Rprof()可以做这件事情

我们只看lineprof,它的工作原理其实就是执行一段代码，然后每隔几毫秒暂停执行，然后指出哪个函数正在执行，这种分析器叫做sampling or statistical profiler.

假设我们有这段代码

library(lineprof)<br />
f <- function() {<br />
pause(0.1)<br />
g()<br />
h()<br />
}<br />
g <- function() {<br />
pause(0.1)<br />
h()<br />
}<br />
h <- function() {<br />
pause(0.1)<br />
}

作者假设我们有个每0.1秒暂停的分析器，指出哪个函数正在执行，并列出调用函数的树calltree.那么我们会看到

f()<br />
f() > g()<br />
f() > g() > h()<br />
f() > h()

用lineprof的结果是

library(lineprof)<br />
source("profiling-example.R")<br />
l <- lineprof(f())<br />
l<br />
#> time alloc release dups ref src<br />
#> 1 0.074 0.001 0 0 profiling.R#2 f/pause<br />
#> 2 0.143 0.002 0 0 profiling.R#3 f/g<br />
#> 3 0.071 0.000 0 0 profiling.R#4 f/h

特别要注意，lineprof用srcrefs对象来匹配分析，而这种对象只有在代码从硬盘装进内存的时候才会创建，而source("..")就是做了这个事情，所以我们lineprof必须用source进来的代码,至于什么是srcrefs,大家暂时不必理解吧，我也不懂[s:11]

我们只看time，表明各段消耗的时间，ref表明哪行，src就是调用树了

更方便的展示是用shine(l),它是利用shiny包更好的展示结果，可以交互式点击查看，由于这段作者给的是图片，论坛里发图片麻烦，大家其实可以自己试验（要装shiny包）

总之，无论从shiny上看，还是从纯print的输出上看，我们都应该可以找到一些时间花的多的瓶颈代码，然后考虑下面的优化措施,再说各种优化措施之前，还得先提下这种profile不能做的是：

1.无法profile C/C++,以及primitive函数或者byte code compiled code

2.有时无法准确的profile 匿名的函数

3.lazy evaluation会让问题变的不太一样，例如：

i <- function() {<br />
pause(0.1)<br />
10<br />
}<br />
j <- function(x) {<br />
x + 10<br />
}<br />
j(i())

由于lazy,i不会先被执行，然后把结果作为参数传给j，而是等调用j的时候才会执行，所以看上去像是j调用了i,这个应该是特别要注意的

好了，下面我们来介绍各种优化手段，这些手段有些可能比较琐碎难记，但大家经常留意以后可能就熟悉了：

所有优化手段的前提是确保优化了结果和原先是相同的，然后再来看时间，所以综合这两点，我们就要分别用

stopifnot() and all.equal() and microbenckmark，例如：

mean1 <- function(x) mean(x)<br />
mean2 <- function(x) sum(x) / length(x)<br />
x <- runif(100)<br />
stopifnot(all.equal(mean1(x), mean2(x)))<br />
microbenchmark(<br />
mean1(x),<br />
mean2(x)<br />
)

这是所有优化要经历的步骤，首先得确保优化的正确，其次就看优化的速度如何。

第一种优化方式：查查别人怎么解决问题的？

1.上 CRAN task views，看有没有问题相近的

2.看看Rcpp依赖树，看看哪些是用C++解决的问题

3.就是自己上网找，作者推荐了rseek,以及stackoverflow,并提示在sof上加[R]搜索

国内基本就是COS啦[s:11]

第二种优化方式：Do as little as possible

一个函数尽可能的少做一些事情，会变得很快，准确的知道输入输出类型也会变快例如

rowSums(),rowMeans()会比apply(）快，因为它少做了很多事情，然后vapply比sapply快，因为输出类型已知，any(x==10)会比10%in%x快，因为测试相等比测试包含快

作者推荐了两个阅读代码的地方，以此来累积知识面，R-help mailing list 和 stackoverflow，大家赶紧加入这两个地方吧[s:11]

然后作者给出了自己一些个人经验，如下：

read.csv(): specify known columns types with colClasses.
factor(): specify known levels with levels.
cut(): don’t generate labels with labels = FALSE if you don’t need them, or even better, use

findInterval() as mentioned in the “see also” section of the documentation.
unlist(x, use.names = FALSE) is much faster than unlist(x).
interaction(): if you only need combinations that exist in the data, use drop = TRUE

具体大家自己去试验，大概就是read.csv把列的类型确定会快，cut不要输出labels,unlist不要加名字等等

还有比较绝的就是，直接不要method dispath,直接自己来调用正确的方法，如下对S3,S4用findMethod找方法

x <- runif(1e2)<br />
microbenchmark(<br />
mean(x),<br />
mean.default(x)<br />
)<br />
#> Unit: microseconds<br />
#> expr min lq median uq max neval<br />
#> mean(x) 9.94 11.80 12.10 12.50 50.5 100<br />
#> mean.default(x) 2.52 3.16 3.38 3.62 44.6 100

显然这样做很危险，这要求大家非常熟悉自己要输入什么输出什么，如果达到这种境界，就可以修改R中原来的东西，使用自己的快速版本，作者举了个例子，如果一个list内容是长度相等的向量，那可以如下：

quickdf <- function(l) {<br />
class(l) <- "data.frame"<br />
attr(l, "row.names") <- .set_row_names(length(l[[1]]))<br />
l<br />
}<br />
l <- lapply(1:26, function(i) runif(1e3))<br />
names(l) <- letters<br />
microbenchmark(<br />
quickdf(l),<br />
as.data.frame.list(l),<br />
as.data.frame(l)<br />
)<br />
#> Unit: microseconds<br />
#> expr min lq median uq max neval<br />
#> quickdf(l) 25.3 30.6 35.6 42 59.1 100<br />
#> as.data.frame.list(l) 2,180.0 2,270.0 2,380.0 2,490 4,200.0 100<br />
#> as.data.frame(l) 2,200.0 2,300.0 2,410.0 2,550 6,150.0 100

这时间对比！作者逆天了！[s:11]，作者解释as.data.frame()会做很多费事的事情吧，例如把每个东西都变成数据框然后rbind起来，不过，快是快，但得小心

quickdf(list(x = 1, y = 1:2))<br />
#> Warning: corrupt data frame: columns will be truncated or padded with NAs<br />
#> x y<br />
#> 1 1 1

前面说了，必须得长度相同，这里不一致就出问题，而大家可以自己试验as.data.frame是不会出问题的，它会自动采用某种形式，这里其实就是repeat x。

但是我想这仍然阻止不了大家的热情，因为大家自己写程序，如果比较熟悉输入输出，肯定是希望能够有quickdf这种快速版本的，作者就说了，其实是自己花了蛮多时间读源代码然后一条条的去掉源代码的东西拿出最后的结果[s:12]

怎么样？大神也这样做，你是不是有信心了[s:11]我反正有了

好，作者就再给一例，展示了他是如何假定输入是向量，然后一步步修改原diff函数的，我觉得这个“技术”非常实用，篇幅很长我也决定把它复制粘贴下来：

diff1 <- function (x, lag = 1L, differences = 1L) {<br />
ismat <- is.matrix(x)<br />
xlen <- if (ismat) dim(x)[1L] else length(x)<br />
if (length(lag) > 1L || length(differences) > 1L || lag < 1L || differences < 1L)<br />
stop("'lag' and 'differences' must be integers >= 1")<br />
if (lag * differences >= xlen) {<br />
return(x[0L])<br />
}<br />
r <- unclass(x)<br />
i1 <- -seq_len(lag)<br />
if (ismat) {<br />
for (i in seq_len(differences)) {<br />
r <- r[i1, , drop = FALSE] - r[-nrow(r):-(nrow(r) - lag + 1L), ,<br />
drop = FALSE]<br />
}<br />
} else {<br />
for (i in seq_len(differences)) {<br />
r <- r[i1] - r[-length(r):-(length(r) - lag + 1L)]<br />
}<br />
}<br />
class(r) <- oldClass(x)<br />
r<br />
}

假设我们是向量，我们就去掉了matrix部分

diff2 <- function (x, lag = 1L, differences = 1L) {<br />
xlen <- length(x)<br />
if (length(lag) > 1L || length(differences) > 1L || lag < 1L || differences < 1L)<br />
stop("'lag' and 'differences' must be integers >= 1")<br />
if (lag * differences >= xlen) {<br />
return(x[0L])<br />
}<br />
i1 <- -seq_len(lag)<br />
for (i in seq_len(differences)) {<br />
x <- x[i1] - x[-length(x):-(length(x) - lag + 1L)]<br />
}<br />
x<br />
}

然后再假设我们的differences=1L,大家可以自己弄清楚这个参数什么意思

diff3 <- function (x, lag = 1L) {<br />
xlen <- length(x)<br />
if (length(lag) > 1L || lag < 1L)<br />
stop("'lag' must be integer >= 1")<br />
if (lag >= xlen) {<br />
return(x[0L])<br />
}<br />
i1 <- -seq_len(lag)<br />
x[i1] - x[-length(x):-(length(x) - lag + 1L)]<br />
}

最后假设lag=1L,就有了

diff4 <- function (x) {<br />
xlen <- length(x)<br />
if (xlen <= 1) return(x[0L])<br />
x[-1] - x[-xlen]<br />
}

OK,我们来看下速度：

x <- runif(100)<br />
microbenchmark(<br />
diff1(x),<br />
diff2(x),<br />
diff3(x),<br />
diff4(x)<br />
)<br />
#> Unit: microseconds<br />
#> expr min lq median uq max neval<br />
#> diff1(x) 15.40 17.30 17.80 19.70 68.5 100<br />
#> diff2(x) 12.10 13.80 14.60 15.80 39.2 100<br />
#> diff3(x) 10.30 11.60 12.10 12.40 22.8 100<br />
#> diff4(x) 7.49 8.52 9.02 9.44 15.6 100

这只是个例子，速度提高并没有多少，但是设想，我们以后的问题会经常需要一个函数的特定情况，我们能不能不用那些所谓的参数设置，而直接缩简源代码到符合我们的要求的情况，然后用最精简的版本呢？也许那会是速度上的很大的提高，例如quickdf.

所以这个Do as little as possible绝对值得大家反复去体会学习的

作者又顺带提了下，取数据框的下标不如去取数据框每列的下标

sample_rows <- function(df, i) sample.int(nrow(df), i, replace = TRUE)<br />
# Generate a new data frame containing randomly selected rows<br />
boot_cor1 <- function(df, i) {<br />
sub <- df[sample_rows(df, i), , drop = FALSE]<br />
cor(sub$x, sub$y)<br />
}<br />
# Generate new vectors from random rows<br />
boot_cor2 <- function(df, i ) {<br />
idx <- sample_rows(df, i)<br />
cor(df$x[idx], df$y[idx])<br />
}<br />
df <- data.frame(x = runif(100), y = runif(100))<br />
microbenchmark(<br />
boot_cor1(df, 10),<br />
boot_cor2(df, 10)<br />
)<br />
#> Unit: microseconds<br />
#> expr min lq median uq max neval<br />
#> boot_cor1(df, 10) 205 227 281 316 996 100<br />
#> boot_cor2(df, 10) 121 138 161 175 214 100<br />

我觉得这些都是值得记住并且在自己的代码中经常使用

第三种方式：Vectorise

就是用向量化的思维方式，主要就是用C写的函数以及lapply,apply,Vectorise之类的函数，以及用rowSums来替代apply,最后就是向量化取下标，甚至矩阵取下标会很快，缺点就是快的行为不确定，例如查100个不是10个的10X，1000个不是100个的10X，Vectorise这部分作者写的英语我觉得我理解的不太好[s:11]大家仔细去看原文吧

第四种方式：Advoid copies

作者谈到的事情是, 避免复制，循环中初始分配的不够了，于是就搬地方然后再复制过去，这熟悉C/C++的应该比较熟悉

random_string <- function() {<br />
paste(sample(letters, 50, replace = TRUE), collapse = "")<br />
}<br />
strings10 <- replicate(10, random_string())<br />
strings100 <- replicate(100, random_string())<br />
collapse <- function(xs) {<br />
out <- ""<br />
for (x in xs) {<br />
out <- paste0(out, x)<br />
}<br />
out<br />
}<br />
microbenchmark(<br />
loop10 = collapse(strings10),<br />
loop100 = collapse(strings100),<br />
vec10 = paste(strings10, collapse = ""),<br />
vec100 = paste(strings100, collapse = "")<br />
)<br />
#> Unit: microseconds<br />
#> expr min lq median uq max neval<br />
#> loop10 45.7 47.8 52.4 56.1 96.0 100<br />
#> loop100 1,430.0 1,460.0 1,560.0 1,640.0 2,130.0 100<br />
#> vec10 10.6 11.1 11.7 12.9 21.4 100<br />
#> vec100 78.8 79.4 83.0 90.4 124.0 100

这里的out被赋值一个更大的东西的时候，原来的内存装不下，就会再新地方分配一个内存，然后搬去那里，然后就是R中的Modification in place也要复制，这个下贴就会谈到

第五种方式：Byte code compilation

lapply2 <- function(x, f, ...) {<br />
out <- vector("list", length(x))<br />
for (i in seq_along(x)) {<br />
out[[i]] <- f(x[[i]], ...)<br />
}<br />
out<br />
}<br />
lapply2_c <- compiler::cmpfun(lapply2)<br />
x <- list(1:10, letters, c(F, T), NULL)<br />
microbenchmark(<br />
lapply2(x, is.null),<br />
lapply2_c(x, is.null),<br />
lapply(x, is.null)<br />
)<br />
#> Unit: microseconds<br />
#> expr min lq median uq max neval<br />
#> lapply2(x, is.null) 12.00 14.90 15.70 17.30 49.1 100<br />
#> lapply2_c(x, is.null) 7.52 9.15 9.72 10.50 68.3 100<br />
#> lapply(x, is.null) 5.69 6.91 7.41 8.16 16.9 100<br />

对于我们，实现起来很简单，然后速度确实提高了，所以我觉得以后我会经常用这个方法的[s:11]，但要注意很多情况下速度大概只能提高个5%-10%,base R中的函数都是默认byte code compilation, 这个是在R 2.13.0引入的byte code compiler，它可能提高某些代码速度，有些没用，下面就会提到个没用的

作者展示了一个完整的改造t.test的例子，用来总结目前为止的优化手段：

m <- 1000<br />
n <- 50<br />
X <- matrix(rnorm(m * n, mean = 10, sd = 3), nrow = m)<br />
grp <- rep(1:2, each = n / 2)<br />
system.time(for(i in 1:m) t.test(X[i, ] ~ grp)$stat)<br />
#> user system elapsed<br />
#> 1.83 0.00 1.83<br />
system.time(for(i in 1:m) t.test(X[i, grp == 1], X[i, grp == 2])$stat)<br />
#> user system elapsed<br />
#> 0.354 0.000 0.354

首先发现给formula会比较慢，所以改造第2种：

由于for不存储值，所以用apply

compT <- function(x, grp){<br />
t.test(x[grp == 1], x[grp == 2])$stat<br />
}<br />
system.time(t1 <- apply(X, 1, compT, grp = grp))<br />
#> user system elapsed<br />
#> 0.393 0.000 0.394

然后还记得diff4吗？于是读完t.test.default的源码，其中很多是打印p-value,formats之类的（忘说了，问题是只要t统计量，去掉后就如下：（Do as little as possible)

my_t <- function(x, grp) {<br />
t_stat <- function(x) {<br />
m <- mean(x)<br />
n <- length(x)<br />
var <- sum((x - m) ^ 2) / (n - 1)<br />
list(m = m, n = n, var = var)<br />
}<br />
g1 <- t_stat(x[grp == 1])<br />
g2 <- t_stat(x[grp == 2])<br />
se_total <- sqrt(g1$var / g1$n + g2$var / g2$n)<br />
(g1$m - g2$m) / se_total<br />
}<br />
system.time(t2 <- apply(X, 1, my_t, grp = grp))<br />
#> user system elapsed<br />
#> 0.061 0.000 0.060<br />
stopifnot(all.equal(t1, t2))

这样速度就提高了6x,还不够！我们再用rowMeans之类的 ( Vecterise)

rowtstat <- function(X, grp){<br />
t_stat <- function(X) {<br />
m <- rowMeans(X)<br />
n <- ncol(X)<br />
var <- rowSums((X - m) ^ 2) / (n - 1)<br />
list(m = m, n = n, var = var)<br />
}<br />
g1 <- t_stat(X[, grp == 1])<br />
g2 <- t_stat(X[, grp == 2])<br />
se_total <- sqrt(g1$var / g1$n + g2$var / g2$n)<br />
(g1$m - g2$m) / se_total<br />
}<br />
system.time(t3 <- rowtstat(X, grp))<br />
#> user system elapsed<br />
#> 0.003 0.000 0.003<br />
stopifnot(all.equal(t1, t3))<br />

这是什么速度。。。。1000倍的提高！

最后用(byte code compiler):

rowtstat_bc <- compiler::cmpfun(rowtstat)<br />
microbenchmark(<br />
rowtstat(X, grp),<br />
rowtstat_bc(X, grp),<br />
unit = "ms"<br />
)<br />
#> Unit: milliseconds<br />
#> expr min lq median uq max neval<br />
#> rowtstat(X, grp) 2.65 3.18 3.25 3.41 4.80 100<br />
#> rowtstat_bc(X, grp) 2.64 3.17 3.31 3.46 5.95 100

这里却没什么用了

第六种方式：Parallelise

这部分我说不了，不说了，只说一本书 Parallelise R ，这种书得以后才能看[s:11]

OK，结束！

superdesolator

84.ADR之Performance code 之 Memory

首先装下面的包：

install.packages("ggplot2")<br />
install.packages("pryr")<br />
devtools::install_github("hadley/lineprof")

然后进入主题：

1.对象大小

作者用pryr中的object_size()取代R自带的object.size()因为它更好的考虑到了元素共享以及记录了环境.总之，我们用object_size()来衡量一个对象的大小

我们从最简单的integer vector来入手

<br />
sizes <- sapply(0:50, function(n) object_size(seq_len(n)))<br />
plot(0:50, sizes, xlab = "Length", ylab = "Size (bytes)",<br />
type = "s")

图形不贴了，这里最重要的就是我们看到，长度为0的向量占了40B，其实是因为R中任何长度为0的向量都占40B！

object_size(numeric())<br />
#> 40 B<br />
object_size(logical())<br />
#> 40 B<br />
object_size(raw())<br />
#> 40 B<br />
object_size(list())

这40B怎么分配的呢？首先R中任何对象都有4个组成成分(注意我说的对象和向量）

1.Object metadata(4 bytes).这是用来存储base type以及供调试和内存管理的信息

2.Two pointers. 一个指向R在内存中的前一个对象，一个指向后一个，这使得R的内核函数很容易遍历内存中的R对象，这是一个双向链表（2*8 bytes)

3.A pointer to the attributes (8 bytes) 注意这里是一个指向属性的指针，而不是属性本身

以上3个组成部分（前后指针算一个部分吧）是R中任何对象都要占的内存大小28 bytes.而我们的向量又有额外的3个组成部分：

1.The length of the vectors(4 bytes). 显然通过4个字节，R中应该只能创建2^(32-1)个元素的向量，但在R 3.0.0开始以后可以用4个字节创建 2^52个元素的向量，原因大家暂时不必去看了，我也没看

2.The “true” length of the vector (4 bytes). 这个很少用到，有个应用场合是当对象是被用作环境的hash table时，这个时候它表示真实的分配空间，而the length代表已有空间

3.The data(??bytes). 一个空的向量由于没有data这部分就是0 bytes,否则一个数值向量占8 bytes,整型 4bytes,复数 16 bytes.

综上所述, 一个空的向量就是28+8=36 bytes. 还有4 bytes是为了保持地址对齐（熟悉C/C++的应该知道），赋值给非 8 bytes 倍数的地址会很慢，所以一般的CPU会要求内存地址是大小都是8的倍数

深入分析完向量的组成成分之后,我们减去40 bytes,就得到了向量中的data的大小，作者作图更直观的说明了，这里不作图了，只贴原向量

sizes-40<br />
 [1]   0   8   8  16  16  32  32  32  32  48  48  48  48  64<br />
[15]  64  64  64 128 128 128 128 128 128 128 128 128 128 128<br />
[29] 128 128 128 128 128 136 136 144 144 152 152 160 160 168<br />
[43] 168 176 176 184 184 192 192 200 200

我们之前说过integer占4b,又存在内存按8b对齐的问题，所以1个元素不会分4b,会分8b,2个元素也分8b,这都是意料之中，3个元素不会分12，会16，4个也是16，到目前都是对的，但是到向量含有5个integer的时候，应该分5*4，然后对齐应该是24，结果是32b,所以这里引出一个很重要的知识，small vector pool,由于每次R向系统请求内存都是比较expensive的开销，所以R就申请一个大块的内存，然后自己管理，这个大块内存就叫small vector pool,这使得R每次为长度不大（不大于128bytes)的小向量分配内存时不用频繁的向操作系统要，而是从svl中要，否则R会变得很慢.于是R为了效率和简单，就只会创建8,16,32,48,64,128bytes长度的向量，而操作系统很善于分配大块内存，所以对于大于128bytes的，R就会直接向操作系统要内存，并且是8倍数.

作者还提到了关于组成部分能被共享的情况，但是举的例子我运行的结果是没有分享[s:12]

<br />
x <- 1:1e6<br />
object_size(x)<br />
#> 4 MB<br />
y <- list(x, x, x)<br />
object_size(y)<br />
#> 4 MB<br />

我运行的结果是12MB,不知道为什么[s:12](后来问了作者说是，R 3.1.0是对的）然后下面的结果是对的：

x1 <- 1:1e6<br />
y1 <- list(1:1e6, 1:1e6, 1:1e6)<br />
object_size(x1)<br />
#> 4 MB<br />
object_size(y1)<br />
#> 12 MB<br />
object_size(x1, y1)<br />
#> 16 MB<br />
object_size(x1) + object_size(y1) == object_size(x1, y1)<br />
#> [1] TRUE

同样的共享情况也发生在字符串上，作者说R还有个global string pool,不同的字符串只会被存储在一个地方，所以下面的对象大小可能会让你不那么惊讶了：

object_size("banana")<br />
#> 96 B<br />
object_size(rep("banana", 10))<br />
#> 216 B

但是有兴趣的读者还是会好奇到底这个大小怎么来的，比如我[s:11]

虽然我并没有探索完全其中的规律，但我尽可能的以作者给的一个练习来展示一下：

作者让我们比较下下面的两个list

vec <- lapply(0:50, function(i) c("ba", rep("na", i)))<br />
str <- lapply(vec, paste0, collapse = "")

我们先看最简单的情况，就是一个字符串里的字符不断增加,也即第2个list的情况：

再分析之前，我想说下C里面的字符数组都以结尾要占个字节才能形成一个字符串，因为R里面好像也是这个情况

R>a<-"1"<br />
R>object_size(a)<br />
96 B<br />
R>b<-"1234567"<br />
R>object_size(b)<br />
96 B<br />
R>c<-"12345678"<br />
R>object_size(c)<br />
104 B

当含有8个字符的时候就要多内存了，说明还有个“隐形”的字符.

OK，我们来增长这个字符看看什么规律：

vec<-1:150<br />
str<-lapply(vec,function(x) paste(rep("1",x),collapse=""))<br />
str<-c("",str)<br />
R>vapply(str,object_size,numeric(1))-88<br />
  [1]   8   8   8   8   8   8   8   8  16  16  16  16  16  16<br />
 [15]  16  16  32  32  32  32  32  32  32  32  32  32  32  32<br />
 [29]  32  32  32  32  48  48  48  48  48  48  48  48  48  48<br />
 [43]  48  48  48  48  48  48  64  64  64  64  64  64  64  64<br />
 [57]  64  64  64  64  64  64  64  64 128 128 128 128 128 128<br />
 [71] 128 128 128 128 128 128 128 128 128 128 128 128 128 128<br />
 [85] 128 128 128 128 128 128 128 128 128 128 128 128 128 128<br />
 [99] 128 128 128 128 128 128 128 128 128 128 128 128 128 128<br />
[113] 128 128 128 128 128 128 128 128 128 128 128 128 128 128<br />
[127] 128 128 136 136 136 136 136 136 136 136 144 144 144 144<br />
[141] 144 144 144 144 152 152 152 152 152 152 152

首先解释为什么减88b,因为一个空的字符串（可能只含)要占96b,然后1到7个字符占96b,显然1个字符占1b,加上字符，所以96-8=88,这88个b应该是一个字符串对象除了包含data的其它部分的大小，而减去之后，我们发现这个规律正好是我们之前介绍的small vector pool的概念，所以如果对于之前介绍的整数型向量减去40b剩下的就是包含的data的大小，那么这里的一个string,也可以理解为一个减去88b之后的包含data的字符型向量，只不过这里一个元素占1个字节，而integer一个元素占4个.

这样我们就弄清楚了一个string的结构了，也就理解了作者那个练习的第2个list中各个部分了，下面我们再来研究作者说的第一个list.

研究之前，我们先研究包含完全不同字符串的向量的对象大小增长情况，还记得作者说的global string pool么，我们先通过全部取不同字符串来不考虑它

vec<-1:150<br />
str<-lapply(vec,function(x) as.character(seq_len(x)))<br />
res<-vapply(str,object_size,numeric(1))<br />
r<-1:150*8<br />
a<-diff(res-r-40)

结果比较大我不贴了，大家可以复制代码到R中运行看结果.首先我生成了长度从1到150不等的包含不同字符串的向量，并且得到了它们的大小，大家可能觉得毫无规律，由于前面我们说到一个string,其实就是长度为1的string向量，然后每个string对象的data内容是字符型的向量，而我们这里的data内容都限制在了8bytes,所以当我们把一个string向量这样拆分时：40 bytes 用来作为向量 + 8 bytes 用来作为data内容 + 剩余的部分我们未知，于是我们对每个string向量减去40以及总的data的内容的大小（由于是1：150，所以每个string向量的元素个数是1：150，减去的大小自然是1:150*8).然后我们得到了结果a，这表示了两部分内容：

一部分是字符串相同不相同带来的影响，另一部分是，字符串在字符串向量中位置带来的影响

为了去掉第二个部分的影响，我们看下面得代码：

str2<-lapply(2:150,function(x) {<br />
    tmp<-as.character(seq_len(x))<br />
    tmp[length(tmp)]<-as.character(x-1)<br />
    tmp<br />
  })<br />
str2<-c("1",str2)<br />
res2<-vapply(str2,object_size,numeric(1))

这个代码的作用是把每个新增加的不同字符换成出现过的相同的字符，然后我们用res-res2,结果全是48.

这说明了，对于一个string向量，在任何位置新增加一个不同的字符，由于不同字符带来的影响都是48 bytes大小，也就是说，在一个string向量的任何位置加入一个没出现过的字符串给整个string对象带来的大小增加是比出现过的在包含同样data大小的情况下是48bytes.

这点我们也可以从a的结果中看出，几乎都是48bytes.至于为什么前10几个元素会大小不一，我无法解释，但这个大小不一的出现对于相同字符也是一样，加入相同字符也是会出现这种大小不一，所以我前面就说了把它归于位置因素的影响.

总之，通过这个练习我们看到，global string pool 确实影响了string vector的大小，并且当我们把一个string向量分解成：40 bytes 向量基本要求 + string对象包含的data的大小（字符型向量）+ 位置影响的大小 + 是否出现过的影响之后，我们发现是否出现过，如果是，这部分就是0仿佛不需要什么额外的信息一样了，如果不是就是48 bytes.这个结论还是比较漂亮的。

OK，上面讨论的过于细致了，我们可以把上面的研究叫做“global string pool对string share带来的影响",然后我们就继续来回到主题：

2.Memory usage and garbage collection

我们可以使用mem_used()来查看内存使用情况，但这和我们系统报告的不一样，因为它不包含R解释器本身，系统和R都比较lazy,等到需要的时候才会回收内存，内存泄露问题，就是对象之间会有被删掉的对象留下的间隔，R只会计算对象占有的，这些小的碎片不会被利用.

然后就是mem_change

mem_change(x <- 1:1e6)<br />
#> 4.01 MB<br />
# We get that memory back when we delete it<br />
mem_change(rm(x))<br />
#> -4 MB

它会查看内存变化情况，正数表示分配，负数表示回收

但要注意即使什么都不做也会有变化

mem_change(NULL)

因为R会追踪你的操作历史

有了mem_change,我们来介绍garbage collection (or GC for short)，R中利用ref来gc,也就是如果没有名字指向一个对象了，就会被gc掉，但如果还有名字指向就不会被gc

mem_change(x <- 1:1e6)<br />
#> 4 MB<br />
mem_change(y <- x)<br />
#>982 B<br />
# Remove x, no memory freed because y is still pointing to it<br />
mem_change(rm(x))<br />
#> 1.42 kB<br />
# Now nothing points to it and the memory can be freed<br />
mem_change(rm(y))<br />
#>-4 MB

我们可以通过gcinfo(TRUE)看到更具体的信息.作者最后还提到了一个可能的内存泄露，平常都是在函数体分配的会被自动释放，而如果是返回formulas and closures 会保留当时的环境也就是函数的exe e,这个知识在环境那部分介绍过了.

f1 <- function() {<br />
x <- 1:1e6<br />
10<br />
}<br />
mem_change(x <- f1())<br />
#> 1.43 kB<br />
object_size(x)<br />
#> 48 B<br />
f2 <- function() {<br />
x <- 1:1e6<br />
a ~ b<br />
}<br />
mem_change(y <- f2())<br />
#> 4 MB<br />
object_size(y)<br />
#> 4 MB<br />
f3 <- function() {<br />
x <- 1:1e6<br />
function() 10<br />
}<br />
mem_change(z <- f3())<br />
#> 4 MB<br />
object_size(z)<br />
#> 4.01 MB

3.Memory profiling with lineprof

首先我们来熟悉下下面的函数

read_delim <- function(file, header = TRUE, sep = ",") {<br />
# Determine number of fields by reading first line<br />
first <- scan(file, what = character(1), nlines = 1,<br />
sep = sep, quiet = TRUE)<br />
p <- length(first)<br />
# Load all fields as character vectors<br />
all <- scan(file, what = as.list(rep("character", p)),<br />
sep = sep, skip = if (header) 1 else 0, quiet = TRUE)<br />
# Convert from strings to appropriate types (never to factors)<br />
all[] <- lapply(all, type.convert, as.is = TRUE)<br />
# Set column names<br />
if (header) {<br />
names(all) <- first<br />
} else {<br />
names(all) <- paste0("V", seq_along(all))<br />
}<br />
# Convert list into data frame<br />
as.data.frame(all)

这个函数最关键的是弄清楚scan,大家可以?scan然后仔细研究一下，基本就是弄清楚what参数，这里读文件会读入成一个character类型的list,然后又把list每个部分的character转为原来的类型，这用的是type.convert函数，as.is参数是说如果转换不成功就转化为factor,as.is=TRUE丢弃了这个行为.

熟悉了函数之后，我们来profiling:

library(ggplot2)<br />
write.csv(diamonds, "diamonds.csv", row.names = FALSE)<br />
library(lineprof)<br />
source("code/read-delim.R")<br />
prof <- lineprof(read_delim("diamonds.csv"))<br />
shine(prof)

结果的t表明时间，a表明分配，r表明释放（只能说明释放的内存在这条语句之前已经没用了），d表明复制的向量数，对照结果，图不贴了，scan读进来2.5mb,很接近2.8mb在硬盘上，这是因为R不必读逗号，也由于global string pool可以减少大小,然后字符转成原来的类型的时候0.6mb,这里gc还没被触发，最后as.data.frame发生了很多次复制，gc也被触发了.这里我们无法准确的说什么时候gc,要在每次分配内存的时候都gc的话，可以torture=TRUE,我们可以通过gctorture ()来获得更具体的信息，不过这会很慢很慢，会让R变得10-100倍的慢.

4.Modification in place

R中有两个机制，一个是 modify in place (mip), 一个是 modify on copy(moc).（后改：应该是copy on modify，下面moc等价于com)

R根据一个对象的reference来判断这个修改对象的时候采用哪个机制，我们可以C++的引用计数加上R中的名字binding来理解.

当修改一个只有一个名字binding的对象的时候,或者叫做名字指向或者叫做对象只有一个ref的时候，就采用mip,否则就采用moc,这是因为“我和他人共有一个对象，我修改了不能影响他人”这个道理，所以通过一个名字binding去修改对象的时候，会先复制对象到一个新地址然后修改然后把该名字再binding到新地址的对象上，其他不变.

这个时候我们就要用到pryr包中的工具：address(),refs(),tracemem().

我们先看refs():

<br />
> x<-1:10<br />
> refs(x)<br />
[1] 1<br />
> y<-x<br />
> refs(x)<br />
[1] 2<br />
> refs(y)<br />
[1] 2<br />
> z<-x<br />
> refs(x)<br />
[1] 2<br />
> rm(z)<br />
> rm(y)<br />
> refs(x)<br />
[1] 2<br />
> z<br />
Error: object 'z' not found<br />
> y<br />
Error: object 'y' not found<br />
><br />

上面说明了refs结果只有1和2之分，1表示只有1个，2表示多个，因为只要多于一个就会moc,比较奇怪的就是rm之后refs(x)还是2.我们再来对rm做实验：

> a<-1:3<br />
> b<-a<br />
> refs(a)<br />
[1] 2<br />
> refs(b)<br />
[1] 2<br />
> rm(a)<br />
> refs(b)<br />
[1] 2<br />
> refs(a)<br />
[1] 0<br />
> rm(b)<br />
> refs(b)<br />
[1] 0<br />

我的理解是：当一个名字绑定到一个对象的时候，名字就代表了这个对象的地址，所以a,b都代表了内存中为 1:3分配的地址，然后rm(a)，a就不代表地址了，所以refs(a）就是0了，因为它可能查找了内存中对象地址的绑定名字发现没有含有a的，对于b也一样，而最大的迷惑在于rm(a)之后，rm(b),refs(b)之前结果是2.

对于这点：

作者的解释是

x <- 1:5<br />
y <- x<br />
rm(y)<br />
# Should really be 1, because we've deleted y<br />
refs(x)<br />
#> [1] 2

refs() is only an estimate.

就这么多了[s:12] 所以，我无法理解为什么。。。（后求助了下作者，说是：refs(x) is still two because 2 really means 2 or more, so R can't reliably decrement.大概是说R本身就不靠谱，我只是把这个不靠谱展现了出来，别怀疑我的refs函数会给出错的结果！[s:11]好吧，大概就这么个意思）

熟悉了refs之后，要说的就是当refs结果是1就mip,否则就moc.这是个非常好的结论，给我们判断什么时候mip,什么时候moc带来非常方便的判断，很容易记住吧！

然后我们要介绍tracemem()

它是跟踪被复制的对象

x <- 1:10<br />
# Prints the current memory location of the object<br />
tracemem(x)<br />
# [1] "<0x7feeaaa1c6b8>"<br />
x[5] <- 6L<br />
y <- x<br />
# Prints where it has moved from and to<br />
x[5] <- 6L<br />
# tracemem[0x7feeaaa1c6b8 -> 0x7feeaaa1c768]:

第一个x[5]<-6没有出现消息说明了这个时候的x是mip，因为只有一个ref,当y<-x之后就moc了

最后介绍address之前，虽然大家早就猜到address的用途是看对象的地址，我们还是先来介绍一个概念就是引用增加，increment the ref count(itfc)，来看看下面的代码：

<br />
# Touching the object forces an increment<br />
f <- function(x) x<br />
{x <- 1:10; f(x); refs(x)}<br />
#> [1] 2<br />
# Sum is primitive, so no increment<br />
{x <- 1:10; sum(x); refs(x)}<br />
#> [1] 1<br />
# f() and g() never evaluate x, so refs don't increment<br />
f <- function(x) 10<br />
g <- function(x) substitute(x)<br />
{x <- 1:10; f(x); refs(x)}<br />
#> [1] 1<br />
{x <- 1:10; g(x); refs(x)}<br />
#> [1] 1

简单总结就是不是primitive函数且会eval参数的话就会itfc,是的话就不会，这个总结的简单，但真的要深入理解应该会很难，至少作者没有详细解释为什么了，我说下我自己的理解：

primitive函数可能类似C++按引用传递直接传递原对象的ref,导致mip.

非primitive函数且eval后会增加一个ref,但它不同C++按值传递参数会直接灸复制，这里是先增加一个ref,等修改的时候再moc.

一些primitive函数包括 [[<-, [<-, @<-, $<-, attr<-, attributes<-, class<-, dim<-, dimnames<-,

names<-, and levels<-.

这些是直接mip的函数，大家多少都用过吧

OK，这个itcf非常重要，为什么呢？仔细读下面的代码

> x<-1:5<br />
> library(pryr)<br />
> refs(x)<br />
[1] 1<br />
> tracemem(x)<br />
[1] "<0x05bbd0d0>"<br />
> refs(x)<br />
[1] 1<br />
> x[2]<-6L<br />
> refs(x)<br />
[1] 1<br />
> address(x)<br />
[1] "0x5bbd0d0"<br />
> refs(x)<br />
[1] 2<br />
> x[2]<-8L<br />
tracemem[0x05bbd0d0 -> 0x062a9f78]:

首先refs肯定得是个非itcf的函数，然后tracemem实验表明也是，所以在使用[<-这个primitive函数时，itcf不会出现，所以就mip,所以不会打印出tracemem的信息，但是注意，我们使用address之后，就出现问题了！

它是个itcf函数,导致address之后，再修改对象会出现复制！！这个就是我为什么先介绍refs，tracemem,itcf最后介绍address的原因，因为大家会陷入这个陷阱：

<br />
> x<-1:3<br />
> refs(x)<br />
[1] 1<br />
> address(x)<br />
[1] "0x6bf3e10"<br />
> x[2]<-6L<br />
> address(x)<br />
[1] "0x6bf4130"

大家会迷惑为什么x被moc了，而不是mip啊，明明refs是1啊，因为我当时就这么做的...[s:11]（我求助了下作者，他说确实是个dumb bug.原话Yes, that's a dumb bug that I introduced by accident.）

最后说下Rstudio里面试验refs(),基本都是2，因为Rstudio的environment browser会自动增加一个ref.

但下面情况我也没理解：

R>{y<-1:4;refs(y)}<br />
[1] 1<br />
R>refs(y)<br />
[1] 2<br />
R>{a<-1;refs(a)}<br />
[1] 2<br />
R>{aa<-1;refs(aa)}<br />
[1] 2<br />
R>{y<-1:4;refs(y)}<br />
[1] 1<br />
R>{z<-2:3;refs(z)}<br />
[1] 1<br />
R>{zz<-3;refs(zz)}<br />
[1] 2

我反正彻底无语了，因为基本上R environment browser都会自动显示我新创建的，但是一个向量就refs为1，一个值就refs为2. 然后全局赋值不在{}里就

R>zzzz<-4:8<br />
R>refs(zzzz)<br />
[1] 2

无语！！[s:12]

OK，这部分暂时总结到这里，反正有几点是肯定的：

1. 在修改之前的语句是refs()为1就mip,否则就是moc

2. 特别注意会使refs增加的函数基本都是非primitive函数，primitive函数基本都不会增加.

3. R对refs减少好像不靠谱

最后的最后，给个例子，非常有用的例子：

x <- data.frame(matrix(runif(100 * 1e4), ncol = 100))<br />
medians <- vapply(x, median, numeric(1))<br />
for(i in seq_along(medians)) {<br />
x[, i] <- x[, i] - medians[i]<br />
}

这个之所以慢是因为，[<-.data.frame方法不是一个primitive方法，所以每次循环x都会被赋复制

证明如下{code]for(i in 1:5) {

x[, i] <- x[, i] - medians

print(c(address(x), refs(x)))

}[/code]

改进的方法是使用primitive的 [[<-.list方法

y <- as.list(x)<br />
for(i in 1:5) {<br />
y[[i]] <- y[[i]] - medians[i]<br />
print(c(address(y), refs(y)))<br />
}

但是作者又给了个错误的示范，因为address()又增加了y的ref，这样导致即使使用了[[<-.list方法由于y的refs是2，所以还是会moc!!

好吧，这章后面的部分，我是看的比较头疼.... 算是结束了吧，最后ADR还有两章 RCPP和R C API ！

superdesolator

84.ADR之Performance code 之 Rcpp

本来打算看完Rcpp包自带的整个文档再写的，但是还是决定先按照ADR中介绍的写一下，以后应该会深入学习Rcpp这个包的（因为我是C++派的），到时候我会抓出其中的精髓介绍给大家的[s:11]

说起这个Rcpp,貌似主要是Dirk Eddelbuettel and Romain Francois这两个人开发的，也还有其他3个重要的份量级人物，我只提前两个人是因为，我这么一个大菜鸟去了rcppmaillist,问了一个特新手的问题，就是怎么在Rcpp的cpp中调用R自带的函数，RF很亲切，给了我调用函数的方法，但是没告诉我结果的类型转换问题，于是我又问了结果类型转换怎么做，结果DE把我说了一通[s:18]，大概意思是你怎么啥都不会，其实我是理解DE的，每天上来新手总是问重复的问题，他们这些大神确实会烦，哈哈，这些是题外话，我想展示的是，大家研究Rcpp不懂的地方可以去rcpp邮件列表直接去问作者！

好了，我们来看看Hadley wickham怎么介绍Rcpp包的.

首先，Rcpp可以在几个方面很有用：

1.写loops的时候因为R的loop中如果有copy on modify的情况就完了，每次循环都要复制一整个向量，所以用C++的loop，针对性的使用很少的临时变量，避免一大段一大段的向量复制

2.递归的时候，因为递归需要频繁调用函数，而C++调用函数开销比R小的多

3.利用C++ STL中重要的数据结构和算法的时候

Rcpp应该是写了很多类来代表R中的结构，最简单的就是

Scalar的部分

The scalar equivalents of numeric, integer, character, and logical vectors are: double, int, String, and bool

R中的scalar其实就是长度为1的vector

Vector的部分

NumericVector, IntegerVector,CharacterVector, and LogicalVector.

比较顾名思义了，对于矩阵的部分;

NumericMatrix, IntegerMatrix, CharacterMatrix, and LogicalMatrix

有了这些类之后，这些类已经写好了很多类方法了，我们就可以调用它们的方法，这些类很多操作符也已经重载好了，给个例子：

cppFunction('NumericVector rowSumsC(NumericMatrix x) {<br />
int nrow = x.nrow(), ncol = x.ncol();<br />
NumericVector out(nrow);<br />
for (int i = 0; i < nrow; i++) {<br />
double total = 0;<br />
for (int j = 0; j < ncol; j++) {<br />
total += x(i, j);<br />
}<br />
out[i] = total;<br />
}<br />
return out;<br />
}')<br />
set.seed(1014)<br />
x <- matrix(sample(100), 10)<br />
rowSums(x)<br />
#> [1] 458 558 488 458 536 537 488 491 508 528<br />
rowSumsC(x)<br />
#> [1] 458 558 488 458 536 537 488 491 508 528

对于这个例子，对于NumericMatrix这个类，.ncol(),.nrow()都是写好的类方法作用是我们熟知的R中的作用，然后还展示了NumericVector的构造函数，其中x(i,j)用来取矩阵的元素，我觉得是重载了（）操作符，out，可以预见NumericVector底层用了数组或者是c++的vector吧

这样一来我们基本了解Rcpp的框架了：对应R中的数据结构的类，类方法对象方法，重载操作符.

我们之前展示的是cppFunction包裹起来，这样不方便，所以我们通常这样写cpp文件;

#include <Rcpp.h><br />
using namespace Rcpp;<br />
// [[Rcpp::export]]<br />
double meanC(NumericVector x) {<br />
int n = x.size();<br />
double total = 0;<br />
for(int i = 0; i < n; ++i) {<br />
total += x[i];<br />
}<br />
return total / n;<br />
}<br />
/*** R<br />
library(microbenchmark)<br />
x <- runif(1e5)<br />
microbenchmark(<br />
mean(x),<br />
meanC(x)<br />
)<br />
*/

注意其中的#include <Rcpp.h>

using namespace Rcpp; 必须得在cpp文件中

// [[Rcpp::export]] 必须用在每个要给用户使用的函数前面

最后特别注意的是/*** R */这个也很有用，我们可以在其中写R的测试代码，然后sourceCpp()的时候这些代码被自动运行然后打印出来，要注意R和***有个空格

下面举个复杂点的例子

int f4(Function pred, List x) {<br />
int n = x.size();<br />
for(int i = 0; i < n; ++i) {<br />
LogicalVector res = pred(x[i]);<br />
if (res[0]) return i + 1;<br />
}<br />
return 0;<br />
}

这个例子大家可以想想对应R中的什么功能

好了，我们来做个练习,实现：

diff(). Start by assuming lag 1, and then generalise for lag n.

下面是我的解答，大家可以给出自己的版本

#include <Rcpp.h><br />
using namespace Rcpp;<br />
// [[Rcpp::export]]<br />
NumericVector diffC(NumericVector x ,int lag=1) {<br />
    int n=x.size();<br />
    if(lag>n/2) stop("wrong");<br />
    if(n>1){<br />
    NumericVector out(n-lag);<br />
    for(int i=0;i<n-lag;++i) {<br />
      out[i]=x[i+lag]-x[i];<br />
    }<br />
    return out;<br />
    } else {<br />
      return x;<br />
    }<br />
}<br />

接着我们看看attributes怎么用：

#include <Rcpp.h><br />
using namespace Rcpp;<br />
// [[Rcpp::export]]<br />
NumericVector attribs() {<br />
NumericVector out = NumericVector::create(1, 2, 3);<br />
out.names() = CharacterVector::create("a", "b", "c");<br />
out.attr("my-attr") = "my-value";<br />
out.attr("class") = "my-class";<br />
return out;<br />
}

我们可以用属于类的方法create来通过scalar来创建vector.由于R中每个对象都有attributes，我们可以通过.attr()来查询或修改，当然.names()是名字属性的一个别名，然后class也是一种属性，这些基础的知识大家得知道哈

For S4 objects, .slot() plays a similar role to .attr(). 这个我们就不说了[s:11]因为我看不懂

接着我们来看看List和DataFrame怎么用：

#include <Rcpp.h><br />
using namespace Rcpp;<br />
// [[Rcpp::export]]<br />
double mpe(List mod) {<br />
if (!mod.inherits("lm")) stop("Input must be a linear model");<br />
NumericVector resid = as<NumericVector>(mod["residuals"]);<br />
NumericVector fitted = as<NumericVector>(mod["fitted.values"]);<br />
int n = resid.size();<br />
double err = 0;<br />
for(int i = 0; i < n; ++i) {<br />
err += resid[i] / (fitted[i] + resid[i]);<br />
}<br />
return err / n;<br />
}<br />
mod <- lm(mpg ~ wt, data = mtcars)<br />
mpe(mod)<br />
#> [1] -0.0154

很多R返回的结果都是一个list,然后加一些class属性，我们针对这个最通常情况的Rcpp写法就是，参数为List,然后第一个语句就用.inherits()判断是不是设定的某个class，不是就用stop停止，是的话，就针对我们的需要提取list特定的部分然后用as转成我们需要的类型，再分析，注意这里的as是个模板类型

这个例子很好的展示了怎么分析class为lm的对象，这个流程用来分析s3对象是很重要的

接着我们来看看Function怎么用：

#include <Rcpp.h><br />
using namespace Rcpp;<br />
// [[Rcpp::export]]<br />
RObject callWithOne(Function f) {<br />
return f(1);<br />
}<br />
callWithOne(function(x) x + 1)<br />
#> [1] 2<br />
callWithOne(paste)<br />
#> [1] "1"

这个例子展示了最基本的如何从C++调用R的函数，对于按位置传参基本没什么，但是对于命名参数的传参，Rcpp

给了_[""]的写法：

RObject ff(Function f){<br />
  NumericVector a=NumericVector::create(1,2,3);<br />
  return f(a,_["lag"]=2);<br />
}

如果我们ff(mean)就相当于mean(a,lag=2).

RObject是可以捕捉所有类型，因为我们不知道调用的函数会出现什么样的结果

我们还可以返回一个List

#include <Rcpp.h><br />
using namespace Rcpp;<br />
// [[Rcpp::export]]<br />
List lapply1(List input, Function f) {<br />
int n = input.size();<br />
List out(n);<br />
for(int i = 0; i < n; i++) {<br />
out[i] = f(input[i]);<br />
}<br />
return out;<br />
}

如果看过之前的帖子，或ADR这本书读下来，我们知道这是一个lapply的一个CPP版本

当然我们还有其他类型

There are also classes for many more specialised language objects: Environment, ComplexVector,

RawVector, DottedPair, Language, Promise, Symbol, WeakReference, and so on.

这就需要大家自己去读Rcpp那200多页的pdf了，我还没读呢，我觉得近期是肯定得读的，有可能让我们更了解R的结构

开头我提到了我去rcppmaillist问了一个问题，就是怎么调用R的函数，基本上是这样

？？meanC(NumericVector x) {<br />
    Function mean=Environment("package::base")["mean"];</p>
<p>    return mean(x);<br />
 }<br />

这样就是去base包里取出mean然后赋值给Function mean但是这个代码的问题是结果不知道什么类型，mean(x)返回的是SEXP,我让函数返回RObject都出错... 这个问题等我看完Rcpp再来解决吧，但是这个例子展示了如何调用R里面的函数，通过用包名构造evironment对象，然后提取对应的R函数，这个方法很重要，所以我提前展示这个例子

下面，我们来讨论缺失值怎么处理

#include <Rcpp.h><br />
using namespace Rcpp;<br />
// [[Rcpp::export]]<br />
List scalar_missings() {<br />
int int_s = NA_INTEGER;<br />
String chr_s = NA_STRING;<br />
bool lgl_s = NA_LOGICAL;<br />
double num_s = NA_REAL;<br />
return List::create(int_s, chr_s, lgl_s, num_s);<br />
}<br />
str(scalar_missings())<br />
#> List of 4<br />
#> $ : int NA<br />
#> $ : chr NA<br />
#> $ : logi TRUE<br />
#> $ : num NA

我们来一个个分析每个类型的缺失：

对于Integers，缺失值被当作最小的整数储存，R中设定了它们一些行为，而C++不会知道，所以evalCpp('NA_INTEGER + 1') 会给出-2147483647这个结果

所以为了得到正常的效果，我们要创建长度为1的IntegerVector

<br />
IntegerVector a=IntegerVector::create(NA_INTEGER)

这样就可以了

对于Doubles,作者说R中的NA是IEEE浮点数NaN的一种特殊形式（C++中NAN），它的表现如下：

表达式涉及到NAN出现FALSE：

evalCpp("NAN == 1")<br />
#> [1] FALSE

但要注意和逻辑值结合的时候

evalCpp("NAN && TRUE")<br />
#> [1] TRUE<br />
evalCpp("NAN || FALSE")<br />
#> [1] TRUE

最后，在数值的上下文中

evalCpp("NAN + 1")<br />
#> [1] NaN

对于Strings,由于String是Rcpp写的类，所以知道怎么对付缺失值

对于Boolean，C++的bool只有false,true而R是FALSE,TRUE,NA，所以要注意如果一个长度为1逻辑向量中含有缺失值就会被转为TRUE

bool f() {<br />
   LogicalVector a=LogicalVector::create(NA_LOGICAL);<br />
   return a;<br />
}

最开始的scalar_missings也展示了同样的效果

为了判断一个向量中的一个值是不是缺失，用类方法is_na():

#include <Rcpp.h><br />
using namespace Rcpp;<br />
// [[Rcpp::export]]<br />
LogicalVector is_naC(NumericVector x) {<br />
int n = x.size();<br />
LogicalVector out(n);<br />
for (int i = 0; i < n; ++i) {<br />
out[i] = NumericVector::is_na(x[i]);<br />
}<br />
return out;<br />
}<br />
is_naC(c(NA, 5.4, 3.2, NA))<br />
#> [1] TRUE FALSE FALSE TRUE

或者用一个语法糖

#include <Rcpp.h><br />
using namespace Rcpp;<br />
// [[Rcpp::export]]<br />
LogicalVector is_naC2(NumericVector x) {<br />
return is_na(x);<br />
}<br />
is_naC2(c(NA, 5.4, 3.2, NA))<br />
#> [1] TRUE FALSE FALSE TRUE

接下来，我们就做个练习：

Rewrite cumsum() and diff() so they can handle missing values. Note that these functions have

slightly more complicated behaviour

我们写cumsum:

NumericVector cumsumC(NumericVector x,bool narm=true,bool reserve=true){<br />
  LogicalVector pos=!is_na(x);<br />
 if(narm){<br />
   NumericVector narmx=x[pos];<br />
   int n=narmx.size();<br />
   NumericVector out(n);<br />
   out[0]=narmx[0];<br />
   for(int i=1;i<n;++i){<br />
     out[i]=out[i-1]+narmx[i];<br />
   }<br />
     if(reserve){<br />
      NumericVector z=clone(x);<br />
      z[pos]=out;<br />
      return z;<br />
     } else {<br />
      return out;<br />
     }<br />
  } else {<br />
    if(pos[0]){<br />
     int m=0;<br />
     int flag=0;<br />
    for(int i=0;i<x.size();++i){<br />
      if(!pos[i]) {<br />
       m=i-1;<br />
       flag=1;<br />
        break;<br />
     }}<br />
    if(!flag) m=x.size()-1;<br />
    NumericVector y(m+1);<br />
    for(int i=0;i<m+1;++i){<br />
      y[i]=x[i];<br />
    }<br />
    NumericVector tmp=cumsumC(y);<br />
    NumericVector out(x.size());<br />
    for(int i=0;i<x.size();++i){<br />
      if(i<=m) out[i]=tmp[i];<br />
      else  out[i]=NA_REAL;<br />
     }<br />
    return out;<br />
    } else {<br />
      NumericVector z=clone(x);<br />
      z[pos]=NA_REAL;<br />
      return z;    }<br />
    }<br />
}

这段代码我觉得主要是由于不熟练，所以写的很长，其次是貌似没有一个vector[beg:end]这样的一个重载，所以比较麻烦(当然，也许Rcpp中有介绍，可我目前不知道）

程序的功能就是narm控制去不去除NA，一旦不去除，就只算到第1个非NA的数，这部分用正常的cumsum逻辑，其余为NA，一旦去除，而不保留，就是直接提取非NA数据按正常的cumsum算，结果也直接显示，但是保留的话就是把结果对应到对应位置上去，其余NA还是NA，大家可以自己测试测试

下面来介绍一下Rcpp sugar

主要分为4类

arithmetic and logical operators

logical summary functions

vector views

other useful functions

第一类，arithmetic and logical operators

其实很多基本的算术与逻辑操作符都被向量化了+ *, -, /, pow, <, <=, >, >=, ==, !=, !.

pdistR <- function(x, ys) {<br />
sqrt((x - ys) ^ 2)<br />
}<br />
#include <Rcpp.h><br />
using namespace Rcpp;<br />
// [[Rcpp::export]]<br />
NumericVector pdistC2(double x, NumericVector ys) {<br />
return sqrt(pow((x - ys), 2));<br />
}

这其实就是C++的重载操作符，让这些操作符针对特定类型进行特定行为，这里就是实现了向量化

Logical summary functions

<br />
any_naR <- function(x) any(is.na(x))<br />
#include <Rcpp.h><br />
using namespace Rcpp;<br />
// [[Rcpp::export]]<br />
bool any_naC(NumericVector x) {<br />
return is_true(any(is_na(x)));<br />
}<br />

这里Rcpp的any返回一个可以被is_true,is_false,is_na转化成bool类型的对象

Vector views

head(), tail(), rep_each(), rep_len(), rev(),

seq_along(), and seq_len()

好处就是R的版本会发生很多次复制，Rcpp版本则不会所以效率很高

Other useful functions （我就复制粘贴了啊）

Math functions: abs(), acos(), asin(), atan(), beta(), ceil(), ceiling(), choose(), cos(), cosh(),

digamma(), exp(), expm1(), factorial(), floor(), gamma(), lbeta(), lchoose(), lfactorial(),

lgamma(), log(), log10(), log1p(), pentagamma(), psigamma(), round(), signif(), sin(), sinh(),

sqrt(), tan(), tanh(), tetragamma(), trigamma(), trunc().

Scalar summaries: mean(), min(), max(), sum(), sd(), and (for vectors) var().

Vector summaries: cumsum(), diff(), pmin(), and pmax().

Finding values: match(), self_match(), which_max(), which_min().

Dealing with duplicates: duplicated(), unique().

d/q/p/r for all standard distributions.

Finally, noNA(x) asserts that the vector x does not contain any missing values, and allows optimisation of some mathematical operations.

可以看到，Rcpp sugar用处很大，也可以预见，将来有更多的更方便的Rcpp sugar出现

最后，我们来看看激动人心的STL

我觉得这部分太重要了，特别是对熟悉C++的同学们，作者基本从迭代器的角度出发

#include <Rcpp.h><br />
using namespace Rcpp;<br />
// [[Rcpp::export]]<br />
double sum3(NumericVector x) {<br />
double total = 0;<br />
NumericVector::iterator it;<br />
for(it = x.begin(); it != x.end(); ++it) {<br />
total += *it;<br />
}<br />
return total;<br />
}

这是最基本的例子

然后就是开始使用CPP的资源了

#include <numeric><br />
#include <Rcpp.h><br />
using namespace Rcpp;<br />
// [[Rcpp::export]]<br />
double sum4(NumericVector x) {<br />
return std::accumulate(x.begin(), x.end(), 0.0);<br />
}

这里我们就使用了numeric当中的一个算法，简单的实现了自己的sum,我们已经可以预见接下来要介绍的是有多方便了[s:11]

算法的部分，algorithm含有很多基于迭代器的有用的算法（我还记得TICPP上的那段话，正是有了迭代器才使得算法才能泛型化）

下面的代码展示了使用algorithm的部分算法来实现R中的findInterval

#include <algorithm><br />
#include <Rcpp.h><br />
using namespace Rcpp;<br />
// [[Rcpp::export]]<br />
IntegerVector findInterval2(NumericVector x, NumericVector breaks) {<br />
IntegerVector out(x.size());<br />
NumericVector::iterator it, pos;<br />
IntegerVector::iterator out_it;<br />
for(it = x.begin(), out_it = out.begin(); it != x.end();<br />
++it, ++out_it) {<br />
pos = std::upper_bound(breaks.begin(), breaks.end(), *it);<br />
*out_it = std::distance(breaks.begin(), pos);<br />
}<br />
return out;<br />
}

数据结构的部分：

The STL provides a large set of data structures: array, bitset, list, forward_list, map, multimap,

multiset, priority_queue, queue, dequeue, set, stack, unordered_map, unordered_set,

unordered_multimap, unordered_multiset, and vector

下面我们介绍3种比较常用的：vector, the unordered_set, and the unordered_map.

vectors

我们看看如何用vector来实现rle:

<br />
#include <Rcpp.h><br />
using namespace Rcpp;<br />
// [[Rcpp::export]]<br />
List rleC(NumericVector x) {<br />
std::vector<int> lengths;<br />
std::vector<double> values;<br />
// Initialise first value<br />
int i = 0;<br />
double prev = x[0];<br />
values.push_back(prev);<br />
lengths.push_back(1);<br />
NumericVector::iterator it;<br />
for(it = x.begin() + 1; it != x.end(); ++it) {<br />
if (prev == *it) {<br />
lengths[i]++;<br />
} else {<br />
values.push_back(*it);<br />
lengths.push_back(1);<br />
i++;<br />
prev = *it;<br />
}<br />
}<br />
return List::create(<br />
_["lengths"] = lengths,<br />
_["values"] = values<br />
);<br />
}

具体不作多解释了

然后就是sets的部分，来实现一个duplicated()的功能

// [[Rcpp::plugins(cpp11)]]<br />
#include <Rcpp.h><br />
#include <unordered_set><br />
using namespace Rcpp;<br />
// [[Rcpp::export]]<br />
LogicalVector duplicatedC(IntegerVector x) {<br />
std::unordered_set<int> seen;<br />
int n = x.size();<br />
LogicalVector out(n);<br />
for (int i = 0; i < n; ++i) {<br />
out[i] = !seen.insert(x[i]).second;<br />
}<br />
return out;<br />
}

值得注意的是unordered_set在C++ 11中，所以得加个// [[Rcpp::plugins(cpp11)]]，然后就是.insert().second这个返回的是插入的值是不是新的，不是新的就相当于重了.first是返回的指向元素的迭代器

最后就是Map来是实现下table

#include <Rcpp.h><br />
using namespace Rcpp;<br />
// [[Rcpp::export]]<br />
std::map<double, int> tableC(NumericVector x) {<br />
std::map<double, int> counts;<br />
int n = x.size();<br />
for (int i = 0; i < n; i++) {<br />
counts[x[i]]++;<br />
}<br />
return counts;<br />
}

我们来做个练习，我们使用CPP的资源来实现which.max

#include <Rcpp.h><br />
#include <algorithm><br />
using namespace Rcpp;<br />
// [[Rcpp::export]]<br />
 int  whichmax(NumericVector x ) {<br />
   double val=*std::max_element<NumericVector::iterator>(x.begin(),x.end());<br />
   int n=x.size();<br />
   for(int i=0;i<n;++i){<br />
     if(val==x[i])<br />
      return i+1;<br />
    }<br />
 }<br />

这个是我自己写的，大家可以参考

最后作者给了个case-study:

模拟gibbs sampler

gibbs_r <- function(N, thin) {<br />
mat <- matrix(nrow = N, ncol = 2)<br />
x <- y <- 0<br />
for (i in 1:N) {<br />
for (j in 1:thin) {<br />
x <- rgamma(1, 3, y * y + 4)<br />
y <- rnorm(1, 1 / (x + 1), 1 / sqrt(2 * (x + 1)))<br />
}<br />
mat[i, ] <- c(x, y)<br />
}<br />
mat<br />
}<br />
#include <Rcpp.h><br />
using namespace Rcpp;<br />
// [[Rcpp::export]]<br />
NumericMatrix gibbs_cpp(int N, int thin) {<br />
NumericMatrix mat(N, 2);<br />
double x = 0, y = 0;<br />
for(int i = 0; i < N; i++) {<br />
for(int j = 0; j < thin; j++) {<br />
x = rgamma(1, 3, 1 / (y * y + 4))[0];<br />
y = rnorm(1, 1 / (x + 1), 1 / sqrt(2 * (x + 1)))[0];<br />
}<br />
mat(i, 0) = x;<br />
mat(i, 1) = y;<br />
}<br />
return(mat);<br />
}

对于从R到CPP的版本，这里主要就是Rcpp的矩阵构造以及rgamma等返回一个向量得通过取下标转成一个scalar

然后用矩阵的()操作符

最后作者简单对比了下R的循环版本，向量化版本和C++的循环版本

首先是循环版本

vacc1a <- function(age, female, ily) {<br />
p <- 0.25 + 0.3 * 1 / (1 - exp(0.04 * age)) + 0.1 * ily<br />
p <- p * if (female) 1.25 else 0.75<br />
p <- max(0, p)<br />
p <- min(1, p)<br />
p<br />
}<br />
vacc1 <- function(age, female, ily) {<br />
n <- length(age)<br />
out <- numeric(n)<br />
for (i in seq_len(n)) {<br />
out[i] <- vacc1a(age[i], female[i], ily[i])<br />
}<br />
out<br />
}

其次是向量化版本

vacc2 <- function(age, female, ily) {<br />
p <- 0.25 + 0.3 * 1 / (1 - exp(0.04 * age)) + 0.1 * ily<br />
p <- p * ifelse(female, 1.25, 0.75)<br />
p <- pmax(0, p)<br />
p <- pmin(1, p)<br />
p<br />
}

最后是CPP循环版本

#include <Rcpp.h><br />
using namespace Rcpp;<br />
double vacc3a(double age, bool female, bool ily){<br />
double p = 0.25 + 0.3 * 1 / (1 - exp(0.04 * age)) + 0.1 * ily;<br />
p = p * (female ? 1.25 : 0.75);<br />
p = std::max(p, 0.0);<br />
p = std::min(p, 1.0);<br />
return p;<br />
}<br />
// [[Rcpp::export]]<br />
NumericVector vacc3(NumericVector age, LogicalVector female,<br />
LogicalVector ily) {<br />
int n = age.size();<br />
NumericVector out(n);<br />
for(int i = 0; i < n; ++i) {<br />
out[i] = vacc3a(age[i], female[i], ily[i]);<br />
}<br />
return out;<br />
}

我们来看看benchmark

<br />
n <- 1000<br />
age <- rnorm(n, mean = 50, sd = 10)<br />
female <- sample(c(T, F), n, rep = TRUE)<br />
ily <- sample(c(T, F), n, prob = c(0.8, 0.2), rep = TRUE)<br />
stopifnot(<br />
all.equal(vacc1(age, female, ily), vacc2(age, female, ily)),<br />
all.equal(vacc1(age, female, ily), vacc3(age, female, ily))<br />
)<br />
microbenchmark(<br />
vacc1 = vacc1(age, female, ily),<br />
vacc2 = vacc2(age, female, ily),<br />
vacc3 = vacc3(age, female, ily)<br />
)<br />
#> Unit: microseconds<br />
#> expr min lq median uq max neval<br />
#> vacc1 7,160.0 7,460 7,590.0 7,960.0 11,700.0 100<br />
#> vacc2 352.0 362 404.0 422.0 758.0 100<br />
#> vacc3 54.2 56 63.5 69.3 79.6 100<br />

相信这绝对可以震撼大家，向量化本来就可以有很大提高了，但是CPP还能有10X的提高！

最后的最后，作者写了如何在R包里使用Rcpp,这个以后我会写，目前我还不会写R包更别提看懂这部分，作者也不忘给我们指明了解更多的方向：Rcpp包专门介绍的pdf 以及 C++的经典著作

OK，结束！这几天比较“忙”，战线拖的有点长[s:11]

yanlinlin82

回复第125楼的 superdesolator：随着修炼ADR，superdesolator 现在已然进阶成为R之大牛了。

superdesolator

回复第126楼的 yanlinlin82：哈哈，大牛别黑我啦，我前几天去Rcpp的maillist上问了个问题，结果被作者之一Dirk Eddelbuettel给批评的体无完肤啊[s:11]，大概意思是：你怎么啥都不懂？[s:18] 不过被大牛黑的感觉还蛮好

superdesolator

回复第110楼的 superdesolator：这里针对这个环境知识，有个讨论

superdesolator

回复第128楼的 superdesolator：针对环境知识，还有个

关于binding environment 名字的讨论因为说到名字，其实就是涉及到绑定环境的知识了，这个帖子深入讨论了一些问题

superdesolator

85.search path的一些探究

见我写在了这里的帖子

PS：ADR一直欠着R的C interface 没写，一直感到自己功力不够，写不出尽可能好点的指引贴，这几天还是决定先按着ADR上介绍的加上自己探究的一些先写着，等以后学完 R internals manul后有机会写个详细的

superdesolator

明明写了一些文字，怕丢失先发出来了，怎么不见了？

superdesolator

86.ADR之Performance code之 R's C interface (一）

浪费感情啊，写了一部分，先发表了，结果没有了，这部分我就不打算写了[s:15]

主要从两部分来介绍：

一是结合整本ADR，加上自己阅读Rinternals.h的理解

二是 ADR的R's C interface 部分

一.Rinternals.h解读

我们几乎顺着Rinternals.h从上往下漂流

1.首先解释 extern "C",（解释了半天结果丢失了[s:12]，这次简短解释下了）

当C++调用C写的模块的时候，C的文件需要加extern,然后C++的头文件需要加extern "C",这样编译C++模块的时候就会把extern "C" {}部分按C的方式编译，然后连接C模块的时候，就可以找到对应的函数了

更具体的解释看这里吧

2.接下来干的事情是对长向量的支持,无论什么情况, R_len_t都是个int,而R_xlen_t在定义了长向量支持后，是个ptrdiff_t,而这个ptrdiff_t本身又是个long int,然后没定义长向量支持，R_xlen_t就是个int,也就是说，我们用R_xlen_t就可以了，不必关心内部的情况,这部分还要注意的就是

typedef struct { R_xlen_t lv_length, lv_truelength; } R_long_vec_hdr_t;

在ADR的memory部分说过，向量的组成有包含的元素长度以及真实长度，所以上面这个结构体就是这两部分的组合了

3.然后我们来看看最基本的类型：

typedef unsigned int SEXPTYPE;</p>
<p>#define NILSXP	     0	  /* nil = NULL */<br />
#define SYMSXP	     1	  /* symbols */<br />
#define LISTSXP	     2	  /* lists of dotted pairs */<br />
#define CLOSXP	     3	  /* closures */<br />
#define ENVSXP	     4	  /* environments */<br />
#define PROMSXP	     5	  /* promises: [un]evaluated closure arguments */<br />
#define LANGSXP	     6	  /* language constructs (special lists) */<br />
#define SPECIALSXP   7	  /* special forms */<br />
#define BUILTINSXP   8	  /* builtin non-special forms */<br />
#define CHARSXP	     9	  /* "scalar" string type (internal only)*/<br />
#define LGLSXP	    10	  /* logical vectors */<br />
/* 11 and 12 were factors and ordered factors in the 1990s */<br />
#define INTSXP	    13	  /* integer vectors */<br />
#define REALSXP	    14	  /* real variables */<br />
#define CPLXSXP	    15	  /* complex variables */<br />
#define STRSXP	    16	  /* string vectors */<br />
#define DOTSXP	    17	  /* dot-dot-dot object */<br />
#define ANYSXP	    18	  /* make "any" args work.<br />
			     Used in specifying types for symbol<br />
			     registration to mean anything is okay  */<br />
#define VECSXP	    19	  /* generic vectors */<br />
#define EXPRSXP	    20	  /* expressions vectors */<br />
#define BCODESXP    21    /* byte code */<br />
#define EXTPTRSXP   22    /* external pointer */<br />
#define WEAKREFSXP  23    /* weak reference */<br />
#define RAWSXP      24    /* raw bytes */<br />
#define S4SXP       25    /* S4, non-vector */</p>
<p>/* used for detecting PROTECT issues in memory.c */<br />
#define NEWSXP      30    /* fresh node creaed in new page */<br />
#define FREESXP     31    /* node released by GC */</p>
<p>#define FUNSXP      99    /* Closure or Builtin or Special */<br />

很显然我们不可能全都理解，但是至少有些常见的R-level的类型可以找到对应的C的类型，例如INTSXP,VECSXP等，但要把向量和标量的区别分清楚.有了这个最基本的定义后，我们来看看，我们怎么把类型和对应的类型描述的字符放在一起，这个表在util.c中

const static struct {<br />
  const char * const str;<br />
  const int type;<br />
}<br />
TypeTable[] = {<br />
{ "NULL",  	NILSXP	   },  /* real types */<br />
{ "symbol",		SYMSXP	   },<br />
{ "pairlist",	LISTSXP	   },<br />
{ "closure",	CLOSXP	   },<br />
{ "environment",	ENVSXP	   },<br />
{ "promise",	PROMSXP	   },<br />
{ "language",	LANGSXP	   },<br />
{ "special",	SPECIALSXP },<br />
{ "builtin",	BUILTINSXP },<br />
{ "char",		CHARSXP	   },<br />
{ "logical",	LGLSXP	   },<br />
{ "integer",	INTSXP	   },<br />
{ "double",		REALSXP	   }, /*-  "real", for R <= 0.61.x */<br />
{ "complex",	CPLXSXP	   },<br />
{ "character",	STRSXP	   },<br />
{ "...",		DOTSXP	   },<br />
{ "any",		ANYSXP	   },<br />
{ "expression",	EXPRSXP	   },<br />
{ "list",		VECSXP	   },<br />
{ "externalptr",	EXTPTRSXP  },<br />
{ "bytecode",	BCODESXP   },<br />
{ "weakref",	WEAKREFSXP },<br />
{ "raw",		RAWSXP },<br />
{ "S4",		S4SXP },<br />
/* aliases : */<br />
{ "numeric",	REALSXP	   },<br />
{ "name",		SYMSXP	   },</p>
<p>{ (char *)NULL,	-1	   }<br />
};

Rinternals.h还从enum的角度定义了一次上面的类型，不管怎么样，大体意思不变的

4.然后就是解释USE_RINTERNALS这个宏定义的作用了，

这里另外开了个帖子基本作用就是为了分开同名宏定义和函数

5.最最重要的部分，就是下面的几个结构体的认识：

struct sxpinfo_struct {<br />
    SEXPTYPE type      :  5;/* ==> (FUNSXP == 99) %% 2^5 == 3 == CLOSXP<br />
			     * -> warning: `type' is narrower than values<br />
			     *              of its type<br />
			     * when SEXPTYPE was an enum */<br />
    unsigned int obj   :  1;<br />
    unsigned int named :  2;<br />
    unsigned int gp    : 16;<br />
    unsigned int mark  :  1;<br />
    unsigned int debug :  1;<br />
    unsigned int trace :  1;  /* functions and memory tracing */<br />
    unsigned int spare :  1;  /* currently unused */<br />
    unsigned int gcgen :  1;  /* old generation number */<br />
    unsigned int gccls :  3;  /* node class */<br />
}; /*		    Tot: 32 */

这个结构体做的事情是用1个int通过位域来携带很多信息，这里面很多东西暂时不必关心，我们可以关注下

SEXPTYPE type : 5;这意味着用5个bit来表示类型，2^5=32，所以最多可以表示32种类型，而大家可以数下SEXPTYPE那里，只有29个类型，完全足够了.

我们再先来看看下面的代码

#define SEXPREC_HEADER \<br />
    struct sxpinfo_struct sxpinfo; \<br />
    struct SEXPREC *attrib; \<br />
    struct SEXPREC *gengc_next_node, *gengc_prev_node

我们看到SEXPREC结构体中将会定义sxpinfo_struct类型的字段sxpinfo,OK,我们再继续先看：

#define TYPEOF(x)	((x)->sxpinfo.type)

大家秒懂了吧[s:11]，原来TYPEOF是个宏，干的事情就是去结构体的某个位域数数，然后去SEXPTYPE那里对号入座而已（我略去了那些宏就是数字这层意思，为了好理解）

然后就是一些具体的结构体，我选一个举个例子：

struct closxp_struct {<br />
    struct SEXPREC *formals;<br />
    struct SEXPREC *body;<br />
    struct SEXPREC *env;<br />
};

这就是一个closure结构体，熟悉函数的应该立马认出formals,body, environment(这个有点出入）

OK，接下来就是描述一个R对象的结构体了

typedef struct SEXPREC {<br />
    SEXPREC_HEADER;<br />
    union {<br />
	struct primsxp_struct primsxp;<br />
	struct symsxp_struct symsxp;<br />
	struct listsxp_struct listsxp;<br />
	struct envsxp_struct envsxp;<br />
	struct closxp_struct closxp;<br />
	struct promsxp_struct promsxp;<br />
    } u;<br />
} SEXPREC, *SEXP;

这就是我们最常见的SEXP这个东西的源头了，它就是上面结构体的指针，其中的SEXPREC_HEADER;已经在前面给出了，那么我们不惜代价的人为展开：

typedef struct SEXPREC {<br />
    struct sxpinfo_struct sxpinfo;<br />
    struct SEXPREC *attrib;<br />
    struct SEXPREC *gengc_next_node, *gengc_prev_node;<br />
    union {<br />
	struct primsxp_struct primsxp;<br />
	struct symsxp_struct symsxp;<br />
	struct listsxp_struct listsxp;<br />
	struct envsxp_struct envsxp;<br />
	struct closxp_struct closxp;<br />
	struct promsxp_struct promsxp;<br />
    } u;<br />
} SEXPREC, *SEXP;

我们回忆memory那里说过的R对象的大小分析的时候提到的，我们来一一对号，sxpinfo就是信息4个字节，attrib就是属性指针，gengc_next_node, gengc_prev_node就是指向前后2个对象（或叫node节点）的指针，这些指针ADR作者告诉我们是8字节一个（我自己没看出为什么8字节一个指针,查看了R internals manul中只说了一个节点的32位平台上的总大小是28字节,然后union占3个words,应该是12个字节，3个指针占3个words,12个字节，加上sxpinfo4个字节，这才是28个字节吧？我会问问ADR作者，PS：后记，作者认为我是正确的[s:11]，在32位机上），后面就是具体的数据结构了

再接下来就是向量意义上的对象（节点）的结构体：

typedef struct VECTOR_SEXPREC {<br />
    SEXPREC_HEADER;<br />
    struct vecsxp_struct vecsxp;<br />
} VECTOR_SEXPREC, *VECSEXP;</p>
<p>typedef union { VECTOR_SEXPREC s; double align; } SEXPREC_ALIGN;<br />

我们来看看vecsxp_struct怎么定义的

struct vecsxp_struct {<br />
    R_len_t	length;<br />
    R_len_t	truelength;<br />
};

回忆起memory那部分说的向量的组成部分，发现ADR作者说的也是吻合的.

这部分最后要注意的是一个地址对齐的问题，之前也提到过，这里的方法通过union,当然不熟悉的也不必担心，暂时知道就可以了

6.解释do {} while(0)的意义

这其实就是个一个C语言中宏定义的技巧，这里顺带提下

7.各种宏或者函数定义

说完上面的几点之后，接下来就几乎都是各种实用的access了，例如：

不支持长向量特征的时候

# define LENGTH(x)	(((VECSEXP) (x))->vecsxp.length)

一看就觉得这个函数肯定是lengthof()的底层吧,支持长向量的时候比较复杂，这里不提了

又如非常重要的一个宏：

#define DATAPTR(x)	(((SEXPREC_ALIGN *) (x)) + 1)

解释是这样的：

Under the generational allocator the data for vector nodes comes immediately after the node tructure, so the data address is a known offset from the node SEXP

再回忆起memory中的向量组成，跟在truelength后面就是data了,那么 ((SEXPREC_ALIGN *) (x)) + 1 就是正好跳过了一个向量节点的部分指向了该向量节点的第一个数据地址（要注意向量节点SEXPREC_ALIGN的含义，以及指针+1的深刻理解）

以上面为依据，我们再看看这个拿到数据指针之后的用法

#define CHAR(x)		((const char *) DATAPTR(x))<br />
#define LOGICAL(x)	((int *) DATAPTR(x))<br />
#define INTEGER(x)	((int *) DATAPTR(x))<br />
#define RAW(x)		((Rbyte *) DATAPTR(x))<br />
#define COMPLEX(x)	((Rcomplex *) DATAPTR(x))<br />
#define REAL(x)		((double *) DATAPTR(x))<br />
#define STRING_ELT(x,i)	((SEXP *) DATAPTR(x))[i]<br />
#define VECTOR_ELT(x,i)	((SEXP *) DATAPTR(x))[i]<br />
#define STRING_PTR(x)	((SEXP *) DATAPTR(x))<br />
#define VECTOR_PTR(x)	((SEXP *) DATAPTR(x))

这些宏定义在第二部分，介绍ADR作者说的R's C interface的还会用到

再接下去就是List Access Macros,例如

#define TAG(e)		((e)->u.listsxp.tagval)<br />
#define CAR(e)		((e)->u.listsxp.carval)<br />
#define CDR(e)		((e)->u.listsxp.cdrval)<br />
#define CAAR(e)		CAR(CAR(e))<br />
#define CDAR(e)		CDR(CAR(e))

要注意这里的list和R level的list不一样，R level的list是上面介绍的向量节点SEXPREC_ALIGN的结构

再来 Closure Access Macros

define FORMALS(x)	((x)->u.closxp.formals)<br />
#define BODY(x)		((x)->u.closxp.body)<br />
#define CLOENV(x)	((x)->u.closxp.env)

不用说很熟悉了

Symbol Access Macros，Environment Access Macros 就不说了

说完了宏定义，下面就有一大堆类似的函数定义了，我拿几个第二部分用到的

void SET_STRING_ELT(SEXP x, R_xlen_t i, SEXP v);<br />
SEXP SET_VECTOR_ELT(SEXP x, R_xlen_t i, SEXP v);<br />
#define CONS(a, b)	cons((a), (b))		/* data lists */<br />
#define LCONS(a, b)	lcons((a), (b))		/* language lists */<br />
SEXP Rf_cons(SEXP, SEXP);<br />
EXP (TAG)(SEXP e);<br />
SEXP (CAR)(SEXP e);<br />
SEXP (CDR)(SEXP e);<br />
SEXP (CAAR)(SEXP e);<br />
#define PROTECT(s)	Rf_protect(s)<br />
#define UNPROTECT(n)	Rf_unprotect(n)<br />
#define UNPROTECT_PTR(s)	Rf_unprotect_ptr(s)<br />
typedef int PROTECT_INDEX;<br />
#define PROTECT_WITH_INDEX(x,i) R_ProtectWithIndex(x,i)<br />
SEXP Rf_allocVector(SEXPTYPE, R_xlen_t);<br />
SEXP Rf_duplicate(SEXP);<br />
 SEXP Rf_mkChar(const char *);<br />
int Rf_ncols(SEXP);<br />
int Rf_nrows(SEXP); #这两个很熟悉吧<br />

8.再来就是一些符号的定义

LibExtern SEXP	R_GlobalEnv;	    /* The "global" environment */</p>
<p>LibExtern SEXP  R_EmptyEnv;<br />
LibExtern SEXP	R_MissingArg;	    /* Missing argument marker */<br />
#ifdef __MAIN__<br />
attribute_hidden<br />
#else<br />
extern<br />
#endif<br />
SEXP	R_RestartToken;     /* Marker for restarted function calls */</p>
<p>/* Symbol Table Shortcuts */<br />
LibExtern SEXP	R_Bracket2Symbol;   /* "[[" */<br />
LibExtern SEXP	R_BracketSymbol;    /* "[" */<br />
LibExtern SEXP	R_BraceSymbol;<br />

这部分不说了，在相当长的往下拉的部分都看不懂，这让我觉得要想真的从C的层面了解R（当然Rinernals.h只是一个头文件啊还有很多其他头文件涉及到R.h,Rconfig.h什么之类的就更不懂了，而且还有很多个实现文件，要想能把这么多东西拼在一起，R level上都得再熟悉个几年，等pure R的东西都几乎知道了，才能算真正开始研究这些C API，包括改变R，不过这对我来说至少还有4年左右的距离吧[s:12]，或者我也不会那么底层，能熟练R level就已经够了[s:11]）

9.在经历第8部分大段大段完全不知道说什么之后，在整个Rinternals.h头文件的最后终于看懂了宏的名字的取名，以及一些逻辑判断宏定义

#define acopy_string		Rf_acopy_string<br />
#define alloc3DArray            Rf_alloc3DArray<br />
#define allocArray		Rf_allocArray<br />
undef isNull<br />
#define isNull(s)	(TYPEOF(s) == NILSXP)<br />
#undef isSymbol<br />
#define isSymbol(s)	(TYPEOF(s) == SYMSXP)<br />
#undef isLogical<br />
#define isLogical(s)	(TYPEOF(s) == LGLSXP)

OK，到此为止对Rinternals.h的初步解读结束了！

PS：第二部分准备写在后面的帖子中

Ihavenothing

回复第132楼的 superdesolator：

指针的大小，一般来说在32位系统上是4个字节，64位上是8个字节。

superdesolator

回复第133楼的 Ihavenothing：嗯嗯，Hadley wickham应该说错了，他在说一个节点大小为28 bytes的情况下（显然就是指在32位了），说了一个指针是8 bytes, 按照我的计算和常理应该是4 bytes

信息4 bytes +属性指针4 bytes +前后2个指针 2*4 bytes + 一个union的大小(3个指针）3*4 bytes=28 bytes.

作者以为是

4+8+2*8=28 bytes,很显然他没有计算union，我发了这个计算公式给他，他说

"Yes,I think that's right on 32-bit."

我就郁闷了，如果他在64-bit上说是8 bytes,那么节点大小应该是52 bytes. 天知道他仔细看了我说的没有[s:12]

我还有次跟他说，他的pryr的包中的refs函数有时莫名其妙的返回2，结果他说“i will look into it",结果到现在都没有就这个问题回复我

superdesolator

87.ADR之Performance code之 R's C interface（二）

我们来讲讲ADR作者介绍的部分

首先作者上来就说，建议使用Rcpp[s:12],其实我也觉得真的要用C/C++来写R的函数的时候我会选Rcpp，不过作为了解，还是值得介绍的

1.R中调用C语言函数

通过.Call()的方式如下

// In C ----------------------------------------<br />
#include <R.h><br />
#include <Rinternals.h><br />
SEXP add(SEXP a, SEXP b) {<br />
SEXP result = PROTECT(allocVector(REALSXP, 1));<br />
REAL(result)[0] = asReal(a) + asReal(b);<br />
UNPROTECT(1);<br />
return result;<br />
}<br />
# In R ----------------------------------------<br />
add <- function(a, b) {<br />
.Call("add", a, b)<br />
}

但是通过使用inline包，我们可以直接合并两步

add <- cfunction(c(a = "integer", b = "integer"), "<br />
SEXP result = PROTECT(allocVector(REALSXP, 1));<br />
REAL(result)[0] = asReal(a) + asReal(b);<br />
UNPROTECT(1);<br />
return result;<br />
")<br />
add(1, 5)<br />
#> [1] 6

2.C的数据结构

其实我第一部分已经写的很仔细了，来看看作者的解释

REALSXP: numeric vector<br />
INTSXP: integer vector<br />
LGLSXP: logical vector<br />
STRSXP: character vector<br />
VECSXP: list<br />
CLOSXP: function (closure)<br />
ENVSXP: environment<br />
CPLXSXP: complex vectors<br />
LISTSXP: “pair” lists. At the R level, you only need to care about the distinction lists and pairlists for<br />
function arguments, but internally they are used in many more places<br />
DOTSXP: ‘…’<br />
SYMSXP: names/symbols<br />
NILSXP: NULL<br />
LANGSXP: language constructs<br />
CHARSXP: “scalar” strings<br />
PROMSXP: promises, lazily evaluated function arguments<br />
EXPRSXP: expressions

然后作者的pryr包里提供了sexp_type()来查看

library(pryr)<br />
sexp_type(10L)<br />
#> [1] "INTSXP"<br />
sexp_type("a")<br />
#> [1] "STRSXP"<br />
sexp_type(T)<br />
#> [1] "LGLSXP"<br />
sexp_type(list(a = 1))<br />
#> [1] "VECSXP"<br />
sexp_type(pairlist(a = 1))<br />
#> [1] "LISTSXP

不过,library(pryr)

sexp_type(10L)

#> [1] "INTSXP"

sexp_type("a")

#> [1] "STRSXP"

sexp_type(T)

#> [1] "LGLSXP"

sexp_type(list(a = 1))

#> [1] "VECSXP"

sexp_type(pairlist(a = 1))

#> [1] "LISTSXP [/code]

不过，要注意的是，这里作者说的类型，我在第一部分也说了，都是数字啊！只是用来区别不同的东西，产生类型的概念，真正的有结构体的类型，我们在第一部分的第5小部分全都列出来了

3.创建修改向量

其实很多R level的概念，例如环境，promise参数等,要想弄清楚底层和R level对应很难，因为光作者介绍的基本类型的向量，list等就已经够新手琢磨很久了[s:12]，所以下面看上去很难的东西，说到低，在R的level上就是基本类型向量的操作,根本没有涉及到例如“修改函数体，获得函数参数，NSE,函数式编程，环境”等等的概念

3.1创建向量与垃圾收集

dummy <- cfunction(body = '<br />
SEXP dbls = PROTECT(allocVector(REALSXP, 4));<br />
SEXP lgls = PROTECT(allocVector(LGLSXP, 4));<br />
SEXP ints = PROTECT(allocVector(INTSXP, 4));<br />
SEXP vec = PROTECT(allocVector(VECSXP, 3));<br />
SET_VECTOR_ELT(vec, 0, dbls);<br />
SET_VECTOR_ELT(vec, 1, lgls);<br />
SET_VECTOR_ELT(vec, 2, ints);<br />
UNPROTECT(4);<br />
return vec;<br />
')<br />
dummy()<br />
#> [[1]]<br />
#> [1] 3.458460e-323 6.950929e-310 0.000000e+00 6.950929e-310<br />
#><br />
#> [[2]]<br />
#> [1] TRUE TRUE TRUE TRUE<br />
#><br />
#> [[3]]<br />
#> [1] 187 48 58 142<br />

我们得使用PROTECT来告诉R我们这个东西在用不要垃圾回收了，UNPROTECT(4)就是释放最近保护的4个东西

结合第一部分说的，SET_VECTOR_ELT是个vector Access Functions,再次提醒一定要区分vector和list在C底层的概念，C的底层来说，R的list是个vector,R的pairlist才是个list[s:11]

UNPROTECT_PTR() ，PROTECT_WITH_INDEX() 这两个没看懂，不提了

下面介绍的是memset,allocVector()在输出分配之前不会清理，所以得自己来，至于asInteger,INTEGER()什么意思，第一部分也介绍过了，简单来说这个函数接受的是个长度为1的向量要转成一个int 标量，然后INTEGER是获得数据指针.allocVector的作用，我觉得就是把适当根据输入类型加上一些节点来组成对象，不必深入啦

zeroes <- cfunction(c(n_ = "integer"), '<br />
int n = asInteger(n_);<br />
SEXP out = PROTECT(allocVector(INTSXP, n));<br />
memset(INTEGER(out), 0, n * sizeof(int));<br />
UNPROTECT(1);<br />
return out;<br />
')<br />
zeroes(10);<br />
#> [1] 0 0 0 0 0 0 0 0 0 0

3.2缺失和无穷值

INTSXP: NA_INTEGER<br />
LGLSXP: NA_LOGICAL<br />
STRSXP: NA_STRING<br />
#对于浮点数<br />
#Use ISNA(), ISNAN(), and !R_FINITE() macros to check for missing, NaN, or non-finite<br />
#values. Use the constants NA_REAL, R_NaN, R_PosInf, and R_NegInf to set those values

下面用上面介绍过的知识来写个is.na

is_na <- cfunction(c(x = "ANY"), '<br />
int n = length(x);<br />
SEXP out = PROTECT(allocVector(LGLSXP, n));<br />
for (int i = 0; i < n; i++) {<br />
switch(TYPEOF(x)) {<br />
case LGLSXP:<br />
LOGICAL(out)[i] = (LOGICAL(x)[i] == NA_LOGICAL);<br />
break;<br />
case INTSXP:<br />
LOGICAL(out)[i] = (INTEGER(x)[i] == NA_INTEGER);<br />
break;<br />
case REALSXP:<br />
LOGICAL(out)[i] = ISNA(REAL(x)[i]);<br />
break;<br />
case STRSXP:<br />
LOGICAL(out)[i] = (STRING_ELT(x, i) == NA_STRING);<br />
break;<br />
default:<br />
LOGICAL(out)[i] = NA_LOGICAL;<br />
}<br />
}<br />
UNPROTECT(1);<br />
return out;<br />
')<br />
is_na(c(NA, 1L))<br />
#> [1] TRUE FALSE<br />
is_na(c(NA, 1))<br />
#> [1] TRUE FALSE<br />
is_na(c(NA, "a"))<br />
#> [1] TRUE FALSE<br />
is_na(c(NA, TRUE))<br />
#> [1] TRUE FALSE

R的base::is.na()会对NA和NaN返回TRUE，但上面只对NA

3.3获得向量数据

add_two <- cfunction(c(x = "numeric"), "<br />
int n = length(x);<br />
double *px, *pout;<br />
SEXP out = PROTECT(allocVector(REALSXP, n));<br />
px = REAL(x);<br />
pout = REAL(out);<br />
for (int i = 0; i < n; i++) {<br />
pout[i] = px[i] + 2;<br />
}<br />
UNPROTECT(1);<br />
return out;<br />
")<br />
add_two(as.numeric(1:10))<br />
#> [1] 3 4 5 6 7 8 9 10 11 12

如果深刻理解了我第一部分写的，再加上对allocVector的理解，那么这里就是取出创建的节点向量的数据的指针，然后利用指针操作

4.字符型向量和lists（R的角度）

其实早在memory那章我就探究了R的字符型向量的规律，其实STRSXP的每个元素都是一个CHARSXP，我们用STRING_ELT(x, i)来获得CHARSXP,用CHAR(STRING_ELT(x, i))来获得真实的const char* string，然后用SET_STRING_ELT(x, i, value)来设定值，用mkChar() 来创建一个CHARSXP来插入一个STRSXP向量，用mkString() 来创建一个长度为1的向量，你晕了没[s:11]

举个例子，综合说明一下

<br />
abc <- cfunction(NULL, '<br />
SEXP out = PROTECT(allocVector(STRSXP, 3));<br />
SET_STRING_ELT(out, 0, mkChar("a"));<br />
SET_STRING_ELT(out, 1, mkChar("b"));<br />
SET_STRING_ELT(out, 2, mkChar("c"));<br />
UNPROTECT(1);<br />
return out;<br />
')<br />
abc()<br />
#> [1] "a" "b" "c"

这就使得纯用C写处理R的字符型向量很麻烦，对于lists来说，每个元素类型也不一样，如果用C处理要用到大量的switch也是非常麻烦，所以作者建议大家用Rcpp

5.会修改输入，要注意的问题

add_three <- cfunction(c(x = "numeric"), '<br />
REAL(x)[0] = REAL(x)[0] + 3;<br />
return x;<br />
')<br />
x <- 1<br />
y <- x<br />
add_three(x)<br />
#> [1] 4<br />
x<br />
#> [1] 4<br />
y<br />
#> [1] 4<br />

这就是C底层写的函数和R level写的不同之处了，在R level在函数中我们往往想绕过com的机制，想修改原来的对象，所以经常用subsitute,eval配合各种环境来做好，而C底层，其实就是修改到原来的对象，从C的角度看这个事情，叫做R’s lazy copy-onmodify semantics，别被绕晕，我觉得就是说从C的角度看，不会自动发生com而已

要做正确得这样

add_four <- cfunction(c(x = "numeric"), '<br />
SEXP x_copy = PROTECT(duplicate(x));<br />
REAL(x_copy)[0] = REAL(x_copy)[0] + 4;<br />
UNPROTECT(1);<br />
return x_copy;<br />
')<br />
x <- 1<br />
y <- x<br />
add_four(x)<br />
#> [1] 5<br />
x<br />
#> [1] 1<br />
y<br />
#> [1] 1

所以单纯从一个C写的函数来看，R的解释器在运行的时候可不会帮你com机制，除非是你用R写的函数（个人理解）

6.强制转换

回到3.1就已经接触这个强制转换了

asLogical(x): INTSXP -> int  #也许错了，在跟作者沟通中，会更新结果<br />
asInteger(x): INTSXP -> int<br />
asReal(x): REALSXP -> double<br />
CHAR(asChar(x)): STRSXP -> const char*<br />
ScalarLogical(x): int -> LGLSXP<br />
ScalarInteger(x): int -> INTSXP<br />
ScalarReal(x): double -> REALSXP<br />
mkString(x): const char* -> STRSXP

7.长向量

从R 3.0.0开始支持长向量，感兴趣的可以自己去找相关资料

8.Pairlists（R的角度）

再次以R的角度（C的角度）方式总结下：lists (VECSXPs), pairlists (LISTSXPs)

这两个概念是不同的，第一部分的第7小部分介绍了CAR,CDR等

car <- cfunction(c(x = "ANY"), 'return CAR(x);')<br />
cdr <- cfunction(c(x = "ANY"), 'return CDR(x);')<br />
cadr <- cfunction(c(x = "ANY"), 'return CADR(x);')<br />
x <- quote(f(a = 1, b = 2))<br />
# The first element<br />
car(x)<br />
#> f<br />
# Second and third elements<br />
cdr(x)<br />
#> $a<br />
#> [1] 1<br />
#><br />
#> $b<br />
#> [1] 2<br />
# Second element<br />
car(cdr(x))<br />
#> [1] 1<br />
cadr(x)<br />
#> [1] 1

我们可以使用下面的模板来游走pairlists:

count <- cfunction(c(x = "ANY"), '<br />
SEXP el, nxt;<br />
int i = 0;<br />
for(nxt = x; nxt != R_NilValue; el = CAR(nxt), nxt = CDR(nxt)) {<br />
i++;<br />
}<br />
return ScalarInteger(i);<br />
')<br />
count(quote(f(a, b, c)))<br />
#> [1] 4<br />
count(quote(f()))<br />
#> [1] 1

我们也可以自己通过CONS来创建新pairlists,或LCONS来创建新的call，这些在R level上应该蛮简单的，但是下面是用C写的就复杂多了

new_call <- cfunction(NULL, '<br />
SEXP REALSXP_10 = PROTECT(ScalarReal(10));<br />
SEXP REALSXP_5 = PROTECT(ScalarReal(5));<br />
SEXP out = PROTECT(LCONS(install("+"), LCONS(<br />
REALSXP_10, LCONS(<br />
REALSXP_5, R_NilValue #注意这里的结尾<br />
)<br />
)));<br />
UNPROTECT(3);<br />
return out;<br />
')<br />
gctorture(TRUE)<br />
new_call()<br />
#> 10 + 5<br />
gctorture(FALSE)

install用来创建一个symbol,我们还可以用TAG() and SET_TAG()来获得symbol,按我的理解应该是参数的名称

同时，属性也是个pairlists(这是新概念，在R的level上没有接触过的),不过它已经有setAttrib() and getAttrib()这两个写好的帮助函数了

set_attr <- cfunction(c(obj = "SEXP", attr = "SEXP", value = "SEXP"), '<br />
const char* attr_s = CHAR(asChar(attr));<br />
duplicate(obj);<br />
setAttrib(obj, install(attr_s), value);<br />
return obj;<br />
')<br />
x <- 1:10<br />
set_attr(x, "a", 1)<br />
#> [1] 1 2 3 4 5 6 7 8 9 10<br />
#> attr(,"a")<br />
#> [1] 1

这里做的事情就是，通过CHAR（asChar）来从STRSXP->const char *，然后利用install来获得symbol,最后借助帮助函数来完成，要注意的是，setAttrib() and getAttrib()每次都要在属性pairlists上做个线性搜索

另外，classgets(), namesgets(),dimgets(), and dimnamesgets() are the internal versions of the default methods of class<-, names<-, dim<-, and dimnames<-.

这部分的知识非常好，特别是说明了属性其实是个pairlists,所以在R的level上，我们说我们只要知道参数是pairlists就够了，其实在C的底层很多其它东西都是用pairlists实现的

9.输入验证

这个概念其实就是让最底层的用C做，其它的用pure R 来做

add_ <- cfunction(signature(a = "integer", b = "integer"), "<br />
SEXP result = PROTECT(allocVector(REALSXP, 1));<br />
REAL(result)[0] = asReal(a) + asReal(b);<br />
UNPROTECT(1);<br />
return result;<br />
")<br />
add <- function(a, b) {<br />
stopifnot(is.numeric(a), is.numeric(b))<br />
stopifnot(length(a) == 1, length(b) == 1)<br />
add_(a, b)<br />
}

当然有的时候，用C做验证会更方便

10.寻找C源代码

作者给出一套方法,再src/main/names.c中找C函数名字，然后

pryr::show_c_source(.Internal(tabulate(bin, nbins)))

它利用了GitHub代码搜索功能

然后针对.Internal() and .Primitive().它们和.Call()还是不一样的，先看下参数形式

SEXP attribute_hidden do_tabulate(SEXP call, SEXP op, SEXP args,<br />
SEXP rho)

然后看看作者对参数的解释

SEXP call: the complete call to the function. CAR(call) gives the name of the function (as a symbol);

CDR(call) gives the arguments.

SEXP op: an “offset pointer”. This is used when multiple R functions use the same C function. For

example do_logic() implements &, |, and !. show_c_source() prints this out for you.

SEXP args: a pairlist containing the unevaluated arguments to the function.

SEXP rho: the environment in which the call was executed.

最后作者结合PROTECT,input validation,把.Internal() and .Primitive() 转成了.Call()

OK,随着这部分结束，整本ADR全部结束！！ [s:11]

superdesolator

88. R navigation tools

搬砖贴，点这里看原文

该贴主要说明了怎么在R中导航，我依据这个帖子提取它的要点以及加上其他资料，总结下导航的途径：

1.R-bloggers

2.Task views

3.stackoverflow

4.r-maillist

5.Rsitesearch() 或者点这里去搜索

6.Rseek

7.R-wiki （6.7由于google的问题暂时不能用）

8.Gmane's mirrored R-lists

9. Rdocumentation.org

10.sos package

11.github

12. help.search() or ??

13.apropos() (12.13 只搜索search list）

14.ls() (例如 ls("package:utils", pat="zip") ）

15.getAnywhere() (注意它并没有NSE)

16.::: （例如 stats:::predict.poly）

17.find （find("split", simple=FALSE)）

18.sessionInfo()

19.dump (dump("x", file=stdout()))

20.其他常见的 :

str,methods,args,formals,file.choose,getwd,grep,%in% 等等

superdesolator

89.PKD之介绍

从这个贴开始,我们来介绍怎么制作R包

首先，准备工具：

1.R版本大于3.1.0,Rstudio最新，Rtools 3.1.0相关的最新的，并且把旧版本的在PATH中的路径去掉（针对windows,MAC用XCode,linux用r-base-dev之类的开发工具包)

2.下载相关包install.packages(c("devtools", "roxygen2", "testthat", "knitr"))

3.保持devtools最新devtools::install_github("hadley/devtools")

4.Check下面代码返回TRUE一切就OK

library(devtools)<br />
has_devel()

然后先介绍一些最常见的组成部分：

R/: where your R code lives in .R files.

DESCRIPTION: metadata about the package

man/: function documentation

vignettes/: long-form documentation which show how to combine multiple parts of your package to

solve real problems.

NAMESPACE: ensures that your package plays nicely with others.

tests/: stores unit tests that ensure that your package is operating as designed.

data/: sample datasets (or other R objects)

src/: compiled C, C++ and fortran source code

superdesolator

90.PKD之R code

最简单的R包有4个组成部分：

1. An RStudio project file, pkgname.Rproj.

2. An R/ directory.

3. A basic DESCRIPTION file.

4. A basic NAMESPACE file

我们通过devtools::create("path/to/package/pkgname"来创建，我们所有的R代码就要放入R/文件中

下面我们来演示一下，当中会用到dir.create,file.copy,我们要做的事情是：假设我们当前目录有2个.r文件，叫pkd.r,pkd2.r,然后我们想在当前wd中创建一个新目录pkd,并在其中创建一个新包newpkg,最后把那2个.r文件放入R/中，我们的流程(workflow)如下：

dir.create("pkd")<br />
path<-getwd()<br />
rfiles<-c("pkd.r","pkd2.r")<br />
rfilespath<-paste0(path,"/",rfiles)<br />
devtools::create("./pkd/newpkg")<br />
Rdirpath<-paste0(path,"/pkd/newpkg/R")<br />
sapply(rfilespath,file.copy,to=Rdirpath)<br />

create()会创建一个.Rproj文件,如果你通过其它方式（我也不知道）创建了一个包，可以这样加

devtools::use_rstudio("path/to/package"）

.Rproj只是文本文件，里面有一些选项，我们最好通过双击.Rproj文件打开新包工作界面后点右上角的选项project options选择.

用Rstudio project的好处是：

1.每个project是独立的

2.F2（在函数名字的地方按F2就会自动找到定义的地方），Ctrl+.找函数定义功能（注意一定是Rstudio open a project的时候才可以用）

3.Alt+shift+K 所有快捷键

快捷键的功能得慢慢熟悉

还有两个非常重要的组合键：

1.devtools::load_all(),ctrl+shift+L 它会把source直接装入内存，并会保存所有打开的文件

2.Build&reload键，ctrl+shift+b，它比较复杂，并只能用在Rstudio,它会先装包，然后再重起R，最后再library包

具体的区别，在下面的帖子会有详细的解释

下面我们介绍一个避免side effect的概念：

1.不要使用library(),require(),把需要的东西描述在DESCRIPTION文件

2.不要修改全局的options()或者graphics par().把这些功能写入函数让用户自己选择

3.不要写文件到disk,用data/缓冲

主要原因就是：第一，不要背着用户做一些行为，第二，library()如果是在load_all()的时候实验用到了，却没放入包就失效了.（我的理解，作者这段原话我还是贴出来吧：

[quote]There are two reasons to avoid side-effects. The first is pragmatic: these funtions will work while you’re developing a package locally with load_all(), but they won’t work when you release your package. That’s because your R code is run once when the package is built, not every time library() is called. The second is principled: you shouldn’t change global state behind your users’ backs[/quote]）

但是有的时候还是得用side-effects,一般用.onLoad() and .onAttach()，它们的通常一些用法如下：

1.动态导入DLL，这个后面会在NAMESPACE提到，而不用.onLoad()

2.显示重要的信息，例如包被导入的时候的提示信息：

.onAttach <- function(libname, pkgname) {<br />
packageStartupMessage("Welcome to my package")<br />
}

我们可以把上述代码放入pkd.r中，然后ctrl+shift+b,会看到library包的时候出现了提示信息

3.和其他语言交流:[quote]For example, if you use RJava to talk to a .jar file, you

need to call rJava::.jpackage(). To make C++ classes available as reference classes in R with RCpp

modules, you call Rcpp::loadRcppModules([/quote]

4.To register vignette engines with tools::vignetteEngine().

5.设置options()，看看作者的方法：

.onLoad <- function(libname, pkgname) {<br />
op <- options()<br />
op.devtools <- list(<br />
devtools.path = "~/R-dev",<br />
devtools.install.args = "",<br />
devtools.name = "Your name goes here",<br />
devtools.desc.author = '"First Last <first.last@example.com> [aut, cre]"',<br />
devtools.desc.license = "What license is it under?",<br />
devtools.desc.suggests = NULL,<br />
devtools.desc = list()<br />
)<br />
toset <- !(names(op.devtools) %in% names(op))<br />
if(any(toset)) options(op.devtools[toset])<br />
invisible()<br />
}

另外一种side effect是定义S4，方法和泛型的时候，load顺序很重要，这个事情交给了DESRIPTION文件中Collate字段负责，后面会提到

最后，.onLoad,.onAttach的区别后面也会提到，一般会把它们集中放入某个约定文件，zzz.r (理解能力有限，附带原文,当然它也提到了.onUnload：

[quote]If you use .onLoad(), consider using .onUnload() to clean up any side effects. By convention, .onLoad()

and friends are usually saved in a file called zzz.R.[/quote]

NB:如果包要提交到CRAN（这是我的目标啊，我要完成一个包提交到CRAN上！！[s:11]），必须得用ASCII，但仍然可以包含unicode:

x <- "This is a bullet •"<br />
y <- "This is a bullet \u2022"<br />
identical(x, y)<br />
#> [1] TRUE<br />
cat(stringi::stri_escape_unicode(x))<br />
#> This is a bullet \u2022

OK，最后的一部分，我们来理解理解R包的类型以及相关的概念：

R包的类型，source,bundle,binary,installed 四种类型

1.source就是我们develop的时候的样子

2.bundle就是我们devtools::build()之后的.tar.gz,如果我们把它解压，我们会得到和source不太一样的地方：

1)vignettes被创建好，而不是原始的markdown,latex

2)src中的临时文件不在了

3).Rbuildignore中的文件全都不在

3.如果一个人没有R的development工具，我们就要给binary包，我们通过devtools:build(binary=TRUE)来编译（Mac以.tgz结尾,windows,zip)

如果我们解压，则与source的不同是：

1)R目录没有.r文件了,换成3个解析过的函数的有效形式[quote]there are three files that store the parsed functions

in an efficient format. This is basically the result of loading all the R code and then saving the

functions with save(), but with a little extra metadata to make things as fast as possible[/quote]

2)存在Meta目录包含了很多Rds文件（下面又是我翻译无能)

[quote]These contain cached metadata about the package,

like what topics the help files cover and parsed versions of the DESCRIPTION files. (You can use

readRDS() to see exactly what’s in those files). These files make package loading faster by caching

costly computations[/quote]

3)html文件夹

4)src变为libs，其中包含的是编译好的代码，分32bit,64bit

5)inst目录中的内容会被拿到最上层去

作者给了个图来描述区别，非常形象，可惜发图乏力[s:12]

installed packages只是那些被解压到package library的binary package.

最正常的流程应该是source->bundle,bundle->binary,binary->installed.（没有r development工具的，给个binary版本就能解压到正确的目录，而不能从source或者bundle->binary,PS：顺带提下，这就是1年半以前COS里一位叫做helend的大神给我编译过一次RMysql的binary包，我用了一年半的时间，现在才懂什么意思[s:11])

知道包的类型区别之后，我们再来看看安装包的时候，发生了什么，作者也给了个图，但那个图可能不好理解：

1.install.packages(), download binary ->installed

2.install.packages( ,type="source") download bundle->installed

3.R CMD install binary or bundle or source ->installed

4.install() 作者说是R CMD install的wrapper,但只给了source->installed的箭头

5.build() R CMD build的wrapper, source--- build() ---->bundle----- R CMD install--------->installed

6.install_github() download---github----> source---build()------>bundle------R CMD install--------->installed

另外，install_url(),install_gitorious(), install_bitbucket(）也是和github差不多。

综合来说，build()负责把source变成bundle,install.xxxx负责download各种类型的包，然后有些经过build,有些直接变成installed.例如，install.packages(,type="source")从bundle变成installed,而install_github首先用build从source变成bundle生成vignettes，再从bundle->installed.

了解完上面复杂的区别之后，掌握load_all(),Build&reload的区别也就很简单了：

load_all()直接source->in memory

Build&reload source->installed->in memory

library() installed->in memory

最后我们来看看什么是library:

它是包的集合，就是一个installed packages的目录，我们可以通过.libPaths()来看目录，当我们要求library一个包的时候，R会查找.libPaths()中的每个path，library()一个不存在的包会出现error,而require()会打印结果返回FALSE,但仍要记住，制作一个包的时候不要用这两个，后面会介绍用DESRIPTION

这也就是每个project可以独立依赖包而不影响其它project，就是.libPaths不同就可以了

superdesolator

91.PKD之Package metadata

我们来看看DESCRIPTION文件，我们用devtools:create出来的newpkg中，就有一个最小的文件，这种文件看上去如下

Package: mypackage<br />
Title: What the package does (one line)<br />
Version: 0.1<br />
Authors@R: "First Last <first.last@example.com> [aut, cre]"<br />
Description: What the package does (one paragraph)<br />
Depends: R (>= 3.1.0)<br />
License: What license is it under?<br />
LazyData: true

它是一种DCF文件，全名Debian control format,特点就是一个filed + : +value,多行value的话，第二行开始每行要缩进，一般4个spaces，例如：

Description: The description of a package is usually long,<br />
    spanning multiple lines. The second and subsequent lines<br />
    should be indented, usually with four spaces.

再深入之前，得加个作者的注释

[quote](If you’re creating a lot of packages, you can set global options devtools.desc.author , devtools.desc.license, devtools.desc.suggests, and devtools.desc to modify the defaults. See package?devtools for more details.)[/quote]

上面虽然介绍了最基本的一些部分，但是下面要先介绍最基本的DCF文件没有包括的两个非常重要的field:

Imports and Suggests.我们先介绍这两个字段，再一个个介绍其他的：

1.Imports and Suggests

它们看上去是这样：

Imports:<br />
    dplyr,<br />
    ggvis<br />
Suggests:<br />
    dplyr,<br />
    ggvis,

Imports表明我们的newpkg需要什么样的包，这些包如果使用newpkg的用户没有，那么在装newpkg的时候就会把这个字段下的包也同时装了，这点我想大家早有体会，然后我们在R代码里如果需要用户这些包的东西，最好使用::来引用，（当然后面作者说NAMESPACE会提到更好的解决方法），之所以要::引用是因为，装newpkg的时候帮你自动装上这些Imports包而不会帮你自动library它们，无论如何，这些包是我们的newpkg经常需要用到的依赖包

而Suggests中提供的包，则并不是必须的，它们通常会用在引用数据集，或被用作了测试，或者被用来生成vignettes,或者是我们的包里面只用到它的一个函数，这里其实我个人理解是，假设我们的包里有100个函数，然后只有1个函数要用到，而且这个函数也不是经常会被用到的，例如pryr中的object_size（并不是说它不会被经常用到，只是举例）,那么我们可以不把它放在Imports字段，虽然放了对用户来说是万无一失，但是性价比很不高，因为这个百分之一的函数被用户用到的可能性不高，用户没有必要去为了这个极小可能用到的函数而去装整个包

那这个时候，我们得采用作者推荐的模板：

# You need the suggested package for this function<br />
my_fun <- function(a, b) {<br />
  if (!requireNamespace("pkg", quietly = TRUE)) {<br />
    stop("Pkg needed for this function to work. Please install it.",<br />
      call. = FALSE)<br />
  }<br />
}</p>
<p># There's a fallback method if the package isn't available<br />
my_fun <- function(a, b) {<br />
  if (requireNamespace("pkg", quietly = TRUE)) {<br />
    pkg::f()<br />
  } else {<br />
    g()<br />
  }<br />
}<br />

其实就是万一当用户用到的时候，我们可以提示一个信息，让用户再去装,我们最简单的添加这两个字段的方法是 devtools::use_package()：

devtools::use_package("dplyr") # Defaults to imports<br />
#> Adding dplyr to Imports<br />
#> Refer to functions with dplyr::fun()<br />
devtools::use_package("dplyr", "Suggests")<br />
#> Adding dplyr to Suggests<br />
#> Use requireNamespace("dplyr", quietly = TRUE) to test if package is<br />
#>  installed, then use dplyr::fun() to refer to functions.

当然，我们可以手工添加，最后，作者建议，Imports and Suggests 字段的依赖包，最好是每行一个包，然后按字典序排（这可以从上面的函数创建的结果看到，也可以去作者的ggplot2 github source中看到,注意cran上下的那个source其实是bundle，它的DESCRIPTION中这两个字段貌似在一行）

2.Versioning

简单来说吧，形如：

Imports:<br />
    ggvis (>= 0.2),<br />
    dplyr (>= 0.3.0.1)<br />
Suggests:<br />
    MASS (>= 7.3.0)

因为R不能同时装入不同版本的同一个包，所以最好是大于等于某个版本，而不是精确的=

3.Other dependencies

Depends: devtools::create()帮我们自动弄好了R的版本需求，至于R包，作者建议用Imports,然后说明了在R 2.14.0之前貌似只能用Depends，在R 3.1.1以及之前通常得depends:methods（如果要用S4），这个BUG在R 3.2.0被

修复，并且使用

LinkingTo:

如果需要用到链接编译C代码（具体会在后面的帖子中提到）

Enhances:

哪些包被你这个包加强了，例如你为定义在那些包的类写了类方法，作者建议不熟悉的先不管

最后作者提到SystemRequirements，但这个字段不会被自动检测到，只能给个快速草考，我们最好在README中给个详细的介绍

4.Naming your package

我就不介绍作者说的了，因为取名字嘛，发挥个人想象力呗，但是要列出一个必须的条件

[quote]There’s only one formal requirement: the package name can only consist of letters, numbers and . (and it must start with a letter and cannot end with a period)[/quote]

这也就意味着不能用-或_

5.The title and description fields

Title字段只是一行描述，它应该只是一段普通文本（不加标记），最好65个字符内

Description字段更详细点，可以多个句子，但是只能一个段落，如果要分行，每行80个字符内，行缩进4个spaces

大家可以自己去看看ggplot2的这两个字段，这两个字段在当你想要提交包到CRAN上的时候非常重要，因为它会被展示出来，最后要还是要提醒，把更具体的内容放入README.md（后面也会专门提到）

6.Who wrote the package?

使用Authors@R,这个字段的展示功能很丰富，最基本如下

Authors@R: c(<br />
    person("Hadley", "Wickham", email = "hadley@rstudio.com", role = "cre"),<br />
    person("Winston", "Chang", email = "winston@rstudio.com", role = "aut"))

这个字段是包含可以被执行的R代码，一般使用person这个函数，它对具有first,last name的姓名工作的很好，然后加个email地址，最后得介绍下role参数，引用作者：

[quote]cre: the package maintainer (creator), the person you should bother if you have problems.

aut: full authors who have contributed much to the package.

ctb: people who have made smaller contributions, like patches.

cph: copyright holder. This is used if copyright is held by someone other than the author, typically a company (their employer).

(The full list of roles is extremely comprehensive. Should your package have a woodcutter (“wdc”), lyricist (“lyr”) or costume designer (“cst”), rest comfortably that you can correctly describe their role in creating your package.)[/quote]

最后作者还提到当封装C库的时候的一些问题（会在后面提到）以及最好附带一些其他的可获得资源的途径例如：

[quote]URL: http://yihui.name/knitr/

BugReports: https://github.com/yihui/knitr/issues[/quote]

作者也同时建议使用Authors@R来替代Maintained以及Author，另外在CRAN上的email也要确保不是一个maillist而是供个人用的，同时保证不要有任何的过滤规则

7.license

作者介绍的，我也不好翻译，主要有三种吧：

[quote]MIT (v. similar: to BSD 2 and 3 clause licenses): this is a simple and permissive license. It lets people use your code and freely distribute subject to only one restriction: the license must always be distributed with the code.

The MIT license is a “template”, so if you use it, you need License: MIT + file LICENSE, and LICENSE file that looks like this:

YEAR: <Year or years when changes have been made>

COPYRIGHT HOLDER: <Name of the copyright holder>

GPL-2 or GPL-3: these are “copy-left” licenses, which means that any one who distributes your code in a bundle must license the whole bundle in a GPL-compatible way. Additionally anyone who distributes modified versions of your code (derivative works) must also make the source code avaialble. GPL-3 is a little stricter than GPL-2, closing some older loopholes.

CC0: It relinquishes all your rights on the code and data so that it can be freely used by any one for any purpose. This is sometimes called putting it in the public domain, although that term is not well-defined, and not meaningful in all countries.[/quote]

最后一种常见的是给data packages,因为至少在美国，data是没有版权的（应该没理解错吧），原话如下

[quote]This license is most appropriate for data packages. Data, at least in the US, is not copyrightable anyway, so you’re not really giving up much. This license just makes it clear[/quote]

这个概念第一次深入接触肯定会很多不熟悉，但是我觉得慢慢来吧，作者也给出了更多的学习的地方

choosealicense.com

https://tldrlegal.com/

“Writing R Extensions” section on licensing.

我们自己在用别人的代码的时候，得注意服从它们的license,后面会提到针对C source code

最后如果想要在CRAN上发包，得符合这些license

8.最后其他的组成部分

Collate 控制R文件的source顺序

Version 非常重要

LazyData 非常重要以至于被包含在了minimal descriptio

这些会在后面说的，更多的完整的列表可以在"The DESCRIPTION file” section of the [R extensions manual][description]如果我们又想用自己的fields又要放到CRAN上，请注意用一个正确的英文单词并且不存在在那个完整的列表上的名字