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

superdesolator

73.ADR之环境之基础

终于到了系统学习这个概念的章节了，环境的东西分4个部分，先讲第一个部分基础的概念：

由于作者在原文用图形结合的方式把有些概念解释的非常清楚，不怕麻烦的同学可以去读原文

总来的说环境就是一些名字的集合，然后它们指向存放在其他地方的对象

技术上说环境（下简称ev)分为frame和parent,frame就是本ev的内容，parent就是本ev的父ev,还记得名字寻值得规则吗？那就是在当前ev找不到就去父中找，当然父中还有父。

最特殊的4个ev是：

1.globalenv()

2.baseenv()

3.emptyenv()

4.environment()

对于一个frame的内容来说，最重要的概念就是reference semantics.（具体的讨论可以看这个帖子)可以不同名字指向相同对象，也可以不同名字指向含有相同值的不同对象（大家可以参考原文的图形）

对于parent来说，我们用parent.env()来看当前ev的父ev，做个实验

x=environment()<br />
y=parent.env(x)<br />
x<br />
y<br />
while(!is.null(y)){<br />
y=parent.env(y)<br />
print(y) }<br />

我们可以看到当前ev的父ev一直往父ev上走，就可以发现其实ev是一个list结构，每个ev都有一个父ev除了所有ev的父ev emptyenv()，只有它没有父ev.事实上这就是一个单向列表结构，我们可以往上找父ev,但是给一个ev我们去没有办法找它的子env,另外从这个parent.env()链，我们看到R中的ev结构如下

交互的工作界面是globalenv()->最近的library()的包的ev->第2近->。。-〉baseenv()->emptyenv()

原文有更好的图形示例，我们还可以用search()看出类似的结果

那么在当前的globalenv()中创建子env呢？如下

e <- new.env()<br />
# the default parent provided by new.env() is environment from<br />
# which it is called - in this case that's the global environment.<br />
parent.env(e)<br />
#> <environment: R_GlobalEnv><br />
ls(e)<br />
#> character(0)<br />
e$a <- 1<br />
e$b <- 2<br />
ls(e)<br />
#> [1] "a" "b"<br />
e$a<br />
#> [1] 1

从中我们也可以看到ls()这个函数的作用，默认ls不列出以.开头的，加上ls(e,all.names=TRUE)就可以了

同理，我们也可以用ls.str()更具体的看，这个时候要注意删除要用rm("a",envir=e)而不能像以前用e$a=NULL这种方法，另外，我们对于给定的名字可以用$,[[,get从一个环境中取值

e$c <- 3<br />
e$c<br />
#> [1] 3<br />
e[["c"]]<br />
#> [1] 3<br />
get("c", envir = e)<br />
#> [1] 3

最后要说的是，我们可以用exists判断ev中存在否，要注意look up的规则，所以我们可以

x <- 10<br />
exists("x", envir = e)<br />
#> [1] TRUE<br />
exists("x", envir = e, inherits = FALSE)<br />
#> [1] FALSE

这就是说，在e中看x存不存在，结果是存在，但我们并没有在e中创建啊，其实它从e的父ev也就是globalenv()中找到的而我说过R的交互界面就是globalenv()环境，我们通过inherits=FALSE来去掉这个行为

最后的最后，比较两个环境

o compare enviroments, you must use identical() not ==:</p>
<p>identical(globalenv(), environment())<br />
#> [1] TRUE<br />
globalenv() == environment()<br />
#> Error: comparison (1) is possible only for atomic and list types

用前者不能用后者

superdesolator

74.ADR之环境之递归和迭代环境

既然73说到环境是一个单向list结构，我也给出了一个自己写的丑陋的迭代看parent.env()的代码，这里我将给出原文作者的程序框架，分别为递归版本和迭代版本

在介绍框架之前，我先列几个函数;

1.exists 73中介绍过，并且注意参数inherits

2.get 以前也介绍过

3.where 这个函数用来判断一个名字在哪个环境中

其实这些函数都要用到遍历环境树

我先列出作者的代码

<br />
where <- function(name, env = parent.frame()) {<br />
  if (identical(env, emptyenv())) {<br />
    # Base case<br />
    stop("Can't find ", name, call. = FALSE)</p>
<p>  } else if (exists(name, envir = env, inherits = FALSE)) {<br />
    # Success case<br />
    env</p>
<p>  } else {<br />
    # Recursive case<br />
    where(name, parent.env(env))</p>
<p>  }<br />
}

虽然我不太懂递归，但这段代码难度也没那么大，大家应该能读懂，3个情况，1是迭代终止条件：递归到了emptyenv()，2是判断成功返回我们需要的，3是没有判断成功，继续递归

由此拿出递归的框架

f <- function(..., env = parent.frame()) {<br />
  if (identical(env, emptyenv())) {<br />
    # base case<br />
  } else if (success) {<br />
    # success case<br />
  } else {<br />
    # recursive case<br />
    f(..., env = parent.env(env))<br />
  }

我们再看看循环迭代的框架：

is_empty <- function(x) identical(x, emptyenv())</p>
<p>f2 <- function(..., env = parent.frame()) {<br />
  while(!is_empty(env)) {<br />
    if (success) {<br />
      # success case<br />
      return()<br />
    }<br />
    # inspect parent<br />
    env <- parent.env(env)<br />
  }</p>
<p>  # base case<br />
}

我在73中的那段代码差不多这个意思，但写的太丑陋了，没有想到循环的终止条件是和emptyenv()判断，大家有兴趣可以改下我在73中那段代码[s:11]

有了这两个框架，我们来依次实现自己版本的exists,get以及补充的任务fget（判断函数是不是function)

我将依次列出这些代码，并且递归和迭代版本交替，有兴趣的同学可以在看之前自己写一遍,最关键的其实是判断成功的条件

<br />
#recursion for exists<br />
exists=function(name,env=parent.frame(),inherits=FALSE){<br />
if(identical(env,emptyenv())){<br />
stop("cant",call.=FALSE)<br />
}else if( name%in%ls(env) ){   #关键的是这里名字在不在ls(e)中，还记得ls(e)吗？<br />
return(TRUE)<br />
}else {<br />
exists(name,env=parent.env(env),inherits=FALSE)<br />
}<br />
}<br />
#iteration for get<br />
is_empty <- function(x) identical(x, emptyenv())<br />
get=function(name,env=parent.frame()){<br />
while(!is_empty(env)){<br />
if(exists(name,envir=env,inherits=FALSE)){<br />
return(env[[name]])<br />
}<br />
env=parent.env(env)<br />
}<br />
}<br />
#recursion for fget<br />
fget=function(name,env=parent.frame(),inherits=FALSE){<br />
if(identical(env,emptyenv()))<br />
{<br />
stop("cant",call. = FALSE)<br />
}else if( exists(name,envir=env,inherits=inherits)&&is.function(get(name,env))){</p>
<p> return(get(name,env))</p>
<p>}else{<br />
fget(name,env=parent.env(env),inherits=inherits)<br />
}<br />
}</p>
<p>

最后提醒大家的是我这里是一下子列出来，如果你也一下子复制到一个R的工作界面里，我们版本的exists会覆盖R中自带的导致get那里不能用自带的失效，对于fget也一样

superdesolator

75.ADR之环境之函数环境

我只能说原文作者太牛了,能把东西讲的那么好,而且图形结合非常到位,可惜发图不方便,大家可以去看原文,这里我就无图说明了;

函数环境分4种：

1.enclosing environment(ee)

函数都有个定义的地方，这个定义的地方所在的环境就是这个函数的ee,ee非常重要，因为R中寻找值的lexical scoping就是用ee，我们可以通过environment(fname)去看一个function的ee（还是可惜没图）

一个函数有且仅有一个ee，其余的environment都可以有0个，1个或多个

2.binding environments(bes)

函数需要个名字，这个名字在哪个ev(environment)，这个函数的be就在哪个，显然我们可以有0个或1个或多个be,下面举个例子来说明

<br />
y<-1<br />
function(x) { x+y }<br />
#OK 这就是说在globalenv()中定义了一个函数，所以这个无名的函数的ee是globalenv(),自然该ee中还有个y名字<br />
e<-new.env()<br />
e$g<-functiion() 1<br />
#OK 这就是一个被binding到e环境中的名字g的一个，这里就是be是一个ev,ee是另一个ev<br />
f<-function() 1<br />
#OK, 这就是我们最常用的情况，ee和be都是一个，那就是globalenv()

ee是属于一个函数的，它决定了函数如何找值，而be决定了我们如何找到函数

值得说的是，就这2个ev形成了包的环境，包package有2个环境，1个是package environment(pe),1个是namespace envvironment(ne)，据我的理解，ne相当于所有functions的ee+be,而pe相当于所有我们可以用的函数的be（读不懂得看原文的图去）,我们可以验证：

environment(sd)<br />
#> <environment: namespace:stats><br />
whpryr::where("sd")<br />
#> <environment: package:stats><br />
#> attr(,"name")<br />
#> [1] "package:stats"<br />
#> attr(,"path")<br />
#> [1] "/usr/lib/R/library/stats"

由于ne是函数定义的ee，所以它需要的parent.env()是一个叫做imports environment的环境，而这个ie的parent.env()是namespace::base这个的parent.env()就是我们的globalenv()了，绕晕了的同学别急，去看原图就懂了

所以stats::sd()用到了var()但如果我们定义在globalenv()中是不会影响到sd的,因为它首先会去namespace::stats中找，而我们用的sd其实是在package::stats中找的，还记得最后一个library进来的包是当前交互界面的parent.env()吗？

OK，下面我们来说第3种：

3.execution environments (exee)

其实R中每个函数执行的时候会产生一个临时的exee,当执行结束就不存在了，对于这样一个函数

h <- function(x) {<br />
a <- 2<br />
x + a<br />
}<br />
y <- h(1)

调用的时候发生的过程是，exe 会包含函数体，x名字（1存储在别的地方，不理解的话请复习环境基础贴），然后就是的a的名字，然后返回3，赋给globalenv()中的名字y的值。然后整个exe e被撤销。说的很抽象，还是得看原文图。。。

有个情况是例外的，就是所谓的闭包工厂函数

<br />
plus <- function(x) {<br />
function(y) x + y<br />
}<br />
plus_one <- plus(1)<br />
identical(parent.env(environment(plus_one)), environment(plus))<br />
#> [1] TRUE

这个时候调用plus(1)产生的exe e对于plusone来说是它的ee，所以这个exe e是不会消失的。而plus_one只是一个被绑定到了globalenv(),也就是它的be是globalenv(),又被绕晕了，没图真难[s:12]

4.最后要说的是calling environments(ce)

又来了个概念，慢慢熬！[s:11]

h <- function() {<br />
x <- 10<br />
function() {<br />
x<br />
}<br />
}<br />
i <- h()<br />
x <- 20<br />
i()

i

对于这段代码，i的函数体定义在了h()的exe e,因此它的e e就是h()的exe e,然后名字i在globalenv(),所以be就是globalenv(),这样i()调用的时候寻值x肯定在ee中找，那就是h()的exe中，那就是x=10,但是x=20对i()意味着什么呢？答案就是i()发生调用的时候，它还有个ce,而x=20是它的ce中的值，这个时候就来了个很重要的概念，就是dynamic scoping!!如果你真的一路读懂了，那么你会“哦，原来如此！“

那么我们怎么来获取这个值呢？想到get就简单了：

f2 <- function() {<br />
x <- 10<br />
function() {<br />
def <- get("x", environment())<br />
cll <- get("x", parent.frame())<br />
list(defined = def, called = cll)<br />
}<br />
}<br />
g2 <- f2()<br />
x <- 20<br />
str(g2())<br />
#> List of 2<br />
#> $ defined: num 10<br />
#> $ called : num 20

我只能说，太神奇了。。。

对于一个调用树来说，我觉得作者那个例子不太好，我自己举了个例子

x<-1<br />
y <- 10<br />
f <- function() {<br />
y<-2<br />
x<-2<br />
h()<br />
}<br />
h <- function() {<br />
x <- 3<br />
x + y<br />
}<br />
f()

这个例子我觉得应该很经典，h()的ee是globalenv(),而它的ce是f()的exe e,按照lexical scoping,它在自己的exe e中找不到值了，它就会去它的ee中找比如y找到10，而不会去它的ce中找，我们当然也可以用parent.frame()来获得它的ce中的值，方法也在前面的代码中展示过了

最后我通过上面的代码展示下综合的找环境的方法

1.environment(f)找f的ee

2.pryr::where("f")找f的be

3.parent.frame()找f的ce，同时如何在一个函数内调用，它也是这个函数的exe e

x<-1<br />
y <- 10<br />
f <- function() {<br />
y<-2<br />
x<-2<br />
h()<br />
}</p>
<p>h <- function() {<br />
print(parent.frame())<br />
print(parent.frame()$x)<br />
print(pryr::where("h"))<br />
print(environment(h))<br />
x <- 3<br />
x + y<br />
}<br />
f()

终于把这段写完了[s:11]

superdesolator

76.ADR之环境之一些扫尾知识

熬过了73-75,我们来点轻松的

1.R中不能以"_"开头命名，也不能用保留字，以及反引号可以用于任何命名

<br />
_abc <- 1<br />
# Error: unexpected input in "_"<br />
if <- 10<br />
#> Error: unexpected assignment in "if <-" #?Reserved<br />
<code>a + b</code> <- 3<br />
<code>:)</code> <- "smile"<br />
<code></code> <- "spaces

2.关于<-,<<-

<-就是赋值到当前环境，而<<-只会修改父环境中对应的值，如果父环境中没有，那么就在globalenv()中创造它！！

这个知识非常重要，作者给出了一个不要这种副作用的代码

rebind <- function(name, value, env = parent.frame()) {<br />
if (identical(env, emptyenv())) {<br />
stop("Can't find ", name, call. = FALSE)<br />
} else if (exists(name, envir = env, inherits = FALSE))<br />
{<br />
assign(name, value, envir = env)<br />
} else {<br />
rebind(name, value, parent.env(env))<br />
}<br />
}<br />
rebind("a", 10)<br />
#> Error: Can't find a<br />
a <- 5<br />
rebind("a", 10)<br />
a<br />
#> [1] 10

这样就是找到了就修改，没找到继续找直到遇到emptyenv(),它不会在globalenv()中创造一个新的！！

3.Explicit environments 这里的解释在73中提到过，还开了个贴，那个帖子里有额外的收获，但不是对reference semantics这个知识的，这个知识的正确意思是：

<br />
modify <- function(x) {<br />
x$a <- 2<br />
invisible()<br />
}<br />
x_l <- list()<br />
x_l$a <- 1<br />
modify(x_l)<br />
x_l$a<br />
#> [1] 1<br />
x_e <- new.env()<br />
x_e$a <- 1<br />
modify(x_e)<br />
x_e$a<br />
#> [1] 2<br />

就说通过函数可以修改，这个意思就是传参的时候用的是类似C++的引用了，而不是copy-on-modify机制了，后者是大家用的时候都指向同一个对象，一旦哪个修改了，就复制一个，然后修改把值再给修改的，其他的不影响，而这里就是任何一个修改，都会改，这里改的原因是传参了，通过参数修改了

4.

x <- 1<br />
e1 <- new.env()<br />
get("x", envir = e1)<br />
#> [1] 1<br />
e2 <- new.env(parent = emptyenv())<br />
get("x", envir = e2)<br />
#> Error: object 'x' not found

前面也说过，get具有inherits=TRUE去找，不想要这个效果就用e2的方法

5.最后，注意3个用途，1是环境避免复制 advoiding copies 2.用环境来修改包 package state 3.hash map （环境本身就类一个hash map)

superdesolator

77.ADR之Functionals编程

本贴我会把3个小节的内容整合在一起，篇幅可能很长很长，由于这是一个比较系统的概念，部分知识要用到前面的帖子中介绍的，例如： lazy function, function enviroments 等等[s:11]

R的核心其实就是一门函数式编程语言，简称FP

一切的一切，从lapply说起，这是一个底层用C写的，用R code意思如下的functional:

<br />
#lapply(x,f,...) 概念上等价于<br />
out <- vector("list", length(x))<br />
for (i in seq_along(x)) {<br />
  out[[i]] <- f(x[[i]], ...)<br />
}

这是基础中的基础，我们用的大部分概念思想源于这个函数，我的总结就是：

1. 一个函数用在一个数据list的不同部分，这里的一个函数也指包装了多个函数的

summary <- function(x) {<br />
  c(mean(x), median(x), sd(x), mad(x), IQR(x))<br />
}<br />
lapply(df, summary)

2.多个函数用在一个数据上，主要思想是将函数分别存储在list里,这就是list of functions的概念

x <- 1:10<br />
funs <- list(<br />
  sum = sum,<br />
  mean = mean,<br />
  median = median<br />
)<br />
lapply(funs, function(f) f(x))

充分理解这两点，其实就可以灵活的使用lapply了，这是一个双向思维，可以这样理解：既可以把函数用于数据，也可以把数据用于函数

3.结合上述2点，加上Anonymous functions 和Closures的概念

前面两点已经刻画了lapply的双向思维，这里来刻画使用的FUN，我们都是使用mean,同样我们也可以自定义一个函数再用它，或者直接使用lapply(x,function(x) {})这样的代码，这就是Anonymous functions了，所以没什么好说的，下面重点说下Closures，比如我们的FUN是power(2)：

其实就是函数返回函数，但这里其实内容很深，但如果你熟悉了我前面介绍的函数环境的概念以及惰性求值的概念，那么这里你会轻松很多，否则我猜你可能云里雾里[s:12]

power <- function(exponent) {<br />
  function(x) {<br />
    x ^ exponent<br />
  }<br />
}<br />
square <- power(2)<br />
square(2)<br />
#> [1] 4<br />
square(4)<br />
#> [1] 16</p>
<p>cube <- power(3)<br />
cube(2)<br />
#> [1] 8<br />
cube(4)<br />
#> [1] 64

这是一个很简单易懂的例子，这个例子里基本上就是解释了什么是闭包Closures.

这里要结合前面深究一下，也就是squre,cube的ee就是power的exe e，我们可以通过如下2种途径看squre,cube的代码

<br />
#不行的方法<br />
square<br />
#> function(x) {<br />
#>     x ^ exponent<br />
#>   }<br />
#> <environment: 0x2b3c6d8><br />
cube<br />
#> function(x) {<br />
#>     x ^ exponent<br />
#>   }<br />
#> <environment: 0x2d8cda0><br />
#方法之一;<br />
as.list(environment(square))<br />
#> $exponent<br />
#> [1] 2<br />
as.list(environment(cube))<br />
#> $exponent<br />
#> [1] 3<br />
#方法之二<br />
library(pryr)<br />
unenclose(square)<br />
#> function (x)<br />
#> {<br />
#>     x^2<br />
#> }<br />
unenclose(cube)<br />
#> function (x)<br />
#> {<br />
#>     x^3<br />
#> }<br />

再用同样的知识，可以实现一个叫mutable state的概念，代码如下：

new_counter <- function() {<br />
  i <- 0<br />
  function() {<br />
    i <<- i + 1<br />
    i<br />
  }<br />
}<br />
counter_one <- new_counter()<br />
counter_two <- new_counter()</p>
<p>counter_one()<br />
#> [1] 1<br />
counter_one()<br />
#> [1] 2<br />
counter_two()<br />
#> [1] 1

其实就是用Closures的ee是工厂函数（产生闭包的函数）的exe e，并且这个exe e不变的原理（前面介绍过）

这个Closures的概念下面还会提到，因为它也是FO的概念的基础，那里还会提到惰性求值的问题.

所以到这里，1,2,3点就概括了基础的东西，这里省略了很多细节的东西了，不然篇幅太长[s:11]

上面基于lapply,下面从它扩展For loop functionals：

sapply,vapply,Map,mclapply等

其实，这些函数系列，初学很烦，感觉很多，分不清，其实lapply是list进list出，而sapply是vector出,出的类型没有vapply严格，它们都是用lapply实现的，只是对于结果做了一些转变，对于Map，它是多个input而已，对于mclapply貌似很厉害，因为lapply的计算顺序不一定，所以可以实现这个mvapply，但是它不是windows下的东西，最后我们用apply,sweep,outer,tapply去实现操作矩阵，数据框，这些函数都是具体的去学的东西，我并不会详细介绍它们的本质就是指定行列之后，就成了1d的了，至于3维以上的也是一定，指定dim,但暂时用不到就不提了，tapply值得看一下，后面有个作业，我自己完成了自己的split版本

<br />
split2=function(x,y){<br />
 tmp=lapply(unique(y),function(x,data=y) {<br />
    which(data==x)<br />
    } )<br />
 results=lapply(tmp,function(data) x[data])<br />
 names(results)=unique(y)<br />
 return(results)<br />
}

大家做的时候可以参考一下哈，然后就是其它的一些函数Predicate functionals：: Filter(), Find() and Position().Mathematical functionals：integrate()，uniroot()，optimise()。大家就自己去看看吧

这部分，最后我要说这样的概念，文章介绍的是sum的制作过程：

首先写个普通的add函数

然后写个处理NA的函数,确定默认值，这两个函数都是处理2个参数的

最后利用Reduce实现到一个向量上，当然其中还有一些控制参数，我自己做练习实现min，pmin,如果你也遇到，我这里可以给你参考：

<br />
rm_na <- function(x, y, identity) {<br />
  if (is.na(x) && is.na(y)) {<br />
    identity<br />
  } else if (is.na(x)) {<br />
    y<br />
  } else {<br />
    x<br />
  }<br />
}<br />
smaller=function(x,y,na.rm=FALSE){<br />
  if (na.rm && (is.na(x) || is.na(y))) rm_na(x, y, Inf) else ifelse(x>y,y,x)   </p>
<p>}<br />
r_smaller <- function(xs, na.rm = TRUE) {<br />
  Reduce(function(x, y) smaller(x, y, na.rm = na.rm), xs,init=Inf)<br />
}<br />
v_smaller <- function(x, y, na.rm = FALSE) {<br />
  stopifnot(length(x) == length(y), is.numeric(x), is.numeric(y))<br />
  if (length(x) == 0) return(numeric())<br />
  simplify2array(Map(function(x, y) smaller(x, y, na.rm = na.rm), x, y))<br />
}

好了，两个部分介绍完了，最后来介绍FO的概念Function operators，总的来说就是在Closures上下文章[s:11]

比如作者举的例子，让一个函数延迟执行:

delay_by <- function(delay, f) {<br />
  function(...) {<br />
    Sys.sleep(delay)<br />
    f(...)<br />
  }<br />
}<br />
system.time(runif(100))<br />
#>    user  system elapsed<br />
#>       0       0       0<br />
system.time(delay_by(0.1, runif)(100))<br />
#>    user  system elapsed<br />
#>   0.000   0.000   0.102

利用这个基本的框架思想，作者开发了非常多的用法思路，确实惊艳，比如Capturing function invocations

怎么做？

其实就是这样：

ignore <- function(...) NULL<br />
tee <- function(f, on_input = ignore, on_output = ignore) {<br />
  function(...) {<br />
    on_input(...)<br />
    output <- f(...)<br />
    on_output(output)<br />
    output<br />
  }<br />
}<br />
g <- function(x) cos(x) - x<br />
show_x <- function(x, ...) cat(sprintf("%+.08f", x), "\n")</p>
<p># The location where the function is evaluated:<br />
zero <- uniroot(tee(g, on_input = show_x), c(-5, 5))

神奇吗？这个...经常出现在这个框架思路中，非常厉害的思想。上面2个基于Closures的FO，作者叫它Behavioural FOs，总共还有Output FOs，Input FOs，Combining FOs,其实思路都一样，只不过功能不同，对于OFO，我只摘出这个例子

failwith <- function(default = NULL, f, quiet = FALSE) {<br />
  force(f)<br />
  function(...) {<br />
    out <- default<br />
    try(out <- f(...), silent = quiet)<br />
    out<br />
  }<br />
}<br />
log("a")<br />
#> Error: non-numeric argument to mathematical function<br />
failwith(NA, log)("a")<br />
#> [1] NA<br />
failwith(NA, log, quiet = TRUE)("a")<br />
#> [1] NA

神奇吗？它给了一个函数失败时提供默认值的功能，其实它做了什么？它什么都没做！就是用一个Closures把...传给了内部的函数。。。但是用处大吗？会很大。。。作者举了个GLM的例子，大概是说你要对一系模型拟合，显然有些会失败，但你不想失败的时候会停下，你只想记录失败的，回头再集中研究失败的，这个failwith就可以很好的完成任务。

对于Input FOs，我提2个东西，一个是partial,一个是多个参数传入函数

前者其实就是固定特定参数，后者是这样的

splat <- function (f) {<br />
  force(f)<br />
  function(args) {<br />
    do.call(f, args)<br />
  }<br />
}<br />
x <- c(NA, runif(100), 1000)<br />
args <- list(<br />
  list(x),<br />
  list(x, na.rm = TRUE),<br />
  list(x, na.rm = TRUE, trim = 0.1)<br />
)<br />
lapply(args, splat(mean))

对于最后一种FO，其实就是

compose <- function(f, g) {<br />
  function(...) f(g(...))<br />
}

没什么说的了

最后，我要提4个特别重要的概念：

1.base::Vectorize()

sample2 <- Vectorize(sample, "size", SIMPLIFY = FALSE)<br />
sample2(1:5, c(1, 1, 3))

2.Memoisation

<br />
library(memoise)<br />
slow_function <- function(x) {<br />
  Sys.sleep(1)<br />
  10<br />
}<br />
system.time(slow_function())<br />
#>    user  system elapsed<br />
#>   0.000   0.001   1.002<br />
system.time(slow_function())<br />
#>    user  system elapsed<br />
#>   0.001   0.000   1.000<br />
fast_function <- memoise(slow_function)<br />
system.time(fast_function())<br />
#>    user  system elapsed<br />
#>   0.001   0.000   1.000<br />
system.time(fast_function())<br />
#>    user  system elapsed<br />
#>       0       0       0<br />

很显然，这个东西是把函数每次的结果存起来了，再调用同样的东西会直接给结果，特别适合斐波那契数列

3.Lazynesss

<br />
funop <- function(f, otherargs) {<br />
  function(...) {<br />
    # maybe do something<br />
    res <- f(...)<br />
    # maybe do something else<br />
    res<br />
  }<br />
}

对于这种框架，如果放入lapply,Map可能会出问题：

funs <- list(mean = mean, sum = sum)<br />
funs_m <- lapply(funs, delay_by, delay = 0.1)</p>
<p>funs_m$mean(1:10)<br />
#> [1] 55

所以要

delay_by <- function(delay, f) {<br />
  force(f)<br />
  function(...) {<br />
    Sys.sleep(delay)<br />
    f(...)<br />
  }<br />
}</p>
<p>funs_m <- lapply(funs, delay_by, delay = 0.1)<br />
funs_m$mean(1:10)

这看上去难理解，如果你理解了函数是个对象的时候，然后结合前面的惰性求值的介绍，也许不那么难理解了

4.最后的最后

remember <- function() {<br />
  memory <- list()<br />
  f <- function(...) {<br />
    # This is inefficient!<br />
    memory <<- append(memory, list(...))<br />
    invisible()<br />
  }</p>
<p>  structure(f, class = "remember")<br />
}<br />
as.list.remember <- function(x, ...) {<br />
  environment(x)$memory<br />
}<br />
print.remember <- function(x, ...) {<br />
  cat("Remembering...\n")<br />
  str(as.list(x))<br />
}<br />
locs <- remember()<br />
vals <- remember()<br />
zero <- uniroot(tee(g, locs, vals), c(-5, 5))<br />
x <- unlist(as.list(locs))<br />
error <- unlist(as.list(vals))<br />
plot(x, type = "b"); abline(h = 0.739, col = "grey50")

这段代码让我赞叹不已，它又让我震惊了，它结合了函数环境+S3方法的本质（把函数转成list都可以。。。而且这种转换得是多巧妙。。。。）+Closures

PS：这部分写的比较混乱，可能同时写3小结的原因，显得是写给我自己的一样，因为我略去了很多，下次我决定分开写了

superdesolator

78.ADR之DEBUG相关

这部分内容我暂时写不出来，水平不够，看不懂原文[s:12]，而且用的很少很少，基本不用，以后能写的时候再写

superdesolator

78.ADR之DEBUG相关

这部分内容我暂时写不出来，水平不够，看不懂原文[s:12]，而且用的很少很少，基本不用，以后能写的时候再写

superdesolator

79.ADR之Metaprogramming之Non-standard evalution

下面简称NSE,这个方式的核心就是substitute函数，大家一定要深刻理解这句话;

It looks at a function argument and instead of seeing the value, it sees the code used to compute the value

举例如下

f <- function(x) {<br />
  substitute(x)<br />
}<br />
f(1:10)<br />
#> 1:10</p>
<p>x <- 10<br />
f(x)<br />
#> x</p>
<p>y <- 13<br />
f(x + y^2)<br />
#> x + y^2

这句话非常重要，它是理解作者后面会说到的Calling from another function中出现问题的关键，我必须在这里强调一下，因为我看了很久Calling from another function这部分，一直没懂，直到我联想到这句话。

substitue之所可以这样捕捉参数，是因为参数是一个特别的对象叫做promise，这种对象包含需要计算（compute）的表达式（expression)和在哪里compute的环境，大家可以暂时不必太深的理解promise，所以substitue只是提取了这个expression的部分.

在深入介绍substitue之前，大家得知道quote,eval这两个函数：

quote就是返回它的输入

quote(1:10)<br />
#> 1:10<br />
quote(x)<br />
#> x<br />
quote(x + y^2)<br />
#> x + y^2

但是要特别注意

quote(x<-1)<br />
#而不是<br />
#quote(x=1)

关于这个原因，请看我的这个学习贴的第49个帖子

eval就是第一个参数是expression,第二个参数是在哪个环境，第三个是这个环境没找到，该去哪个地方找

x <- 10<br />
eval(quote(x))<br />
#> [1] 10<br />
e <- new.env()<br />
e$x <- 20<br />
eval(quote(x), e)

第3个参数后面会提到作用也很大，第2个参数默认是当前环境也可以是自定义list,data.frame或指定环境

eval(quote(x), list(x = 30))<br />
#> [1] 30<br />
eval(quote(x), data.frame(x = 40))<br />
#> [1] 40

懂了基本用法之后，这里又要强调一个概念，就是注意eval的第一个参数是一个expression

sample_df <- data.frame(a = 1:5, b = 5:1, c = c(5, 3, 1, 4, 1))<br />
a <- 10<br />
eval(quote(a), sample_df)<br />
#> [1] 1 2 3 4 5<br />
eval(a, sample_df)<br />
#> [1] 10<br />
eval(quote(b), sample_df)<br />
#> [1] 5 4 3 2 1<br />
eval(b, sample_df)<br />
#> Error: object 'b' not found

这是非常重要的第2个注意的地方，第一个是开头就说的那句话，这2个地方决定后面的escape hatch问题的理解!我是理解了很久很久。。。。。可以看到，是expression的时候才会去第2个参数中找，否则就规则不一样，这个作者给的例子对于新手来说不是十分全面，我说过这是个非常重要的问题，为了研究清楚，我自己又做了一些实验，结果如下：

<br />
#第一个实验<br />
e=new.env()<br />
e$x=20<br />
eval(x)<br />
eval(quote(x))<br />
eval(x,e)<br />
eval(quote(x),e)<br />

从这个结果看，eval(x,e)的时候，并不会去从e中找，而quote(x)则正常

<br />
#第2个实验<br />
a=list(x=10)<br />
x=20<br />
f=function(){<br />
eval(x)<br />
}<br />
f()<br />
a=list(x=10)<br />
x=20<br />
f=function(){<br />
eval(quote(x))<br />
}<br />
f()<br />
a=list(x=10)<br />
x=20<br />
f=function(){<br />
eval(quote(x),a)<br />
}<br />
f()<br />
a=list(x=10)<br />
x=20<br />
f=function(){<br />
eval(x,a)<br />
}<br />
f()<br />

结果分别是20，20,10,20,这样就很清楚了，我们暂时可以这样理解：eval(x)的时候它按照普通的lexical scoping的规则找变量的值，就好像并不理会，第2个参数设定的环境，永远都是当前的环境，然后再找函数的ee(具体看函数环境的帖子了解什么是ee）,而eval(quote(x),env,enclos)则按我们之前的规则，其中第2个实验的第2个小函数之所以是20，其实也是eval(quote(x),env,enclos)的规则之中，因为默认从当前环境中找，没找到，就从函数的ee中找

这样详细的对比了eval(x),eval(x,env,enclos)，大家应该比较熟悉了

最后这两个函数是相反的

eval(eval(quote(quote(2 + 2))))<br />
#> [1] 4<br />
quote(eval(quote(eval(quote(eval(quote(2 + 2)))))))

大家可以运行第2个例子，结果不一定大家都能说对，我当时做就错了。。。

好了，初步介绍了substitue,quote,eval,没错，初步[s:12]，但愿你还比较清楚...

在更深入之前，我们来看个例子

sample_df <- data.frame(a = 1:5, b = 5:1, c = c(5, 3, 1, 4, 1))</p>
<p>subset(sample_df, a >= 4)<br />

这个subset大家并不陌生吧，我们要做的就是实现自己版本的subset2

其实上面介绍了那么多，我贴下代码，相信大家都会恍然大悟

subset2 <- function(x, condition) {<br />
  condition_call <- substitute(condition)<br />
  r <- eval(condition_call, x)<br />
  x[r, ,drop=F] #如果你还记得subsetting的东西，这里是为了解决只有一列的数据框的时候会出现的问题<br />
}<br />
subset2(sample_df, a >= 4)

让我们来运用前面介绍的substitue,eval的知识来解释这个函数：

condition_call是个普通变量，含有一个substitue的结果，这是一个expression,类似quote(a>=4),第一次eval是一个之前的eval(x)，找到了x，这里还是忘记说了，eval（x),eval(quote(x),env,enclos)的情况还有一种就是x中含有expression,那就先找到这个x，再把它的值也就是一个表达式按第2个eval执行，我们这里就是这个情况，这点也很重要，这是第3个非常重要的点，这样eval(condition_call, x)类似eval(quote(a>=4),x)的结果大家应该可以自己推出来了

前面提到过eval(quote(x),env,enclos)的第3个参数也很有用，下面我们看它用在哪：

y <- 4<br />
x <- 4<br />
condition <- 4<br />
condition_call <- 4</p>
<p>subset2(sample_df, a == 4)<br />
#>   a b c<br />
#> 4 4 2 4<br />
subset2(sample_df, a == y)<br />
#>   a b c<br />
#> 4 4 2 4<br />
subset2(sample_df, a == x)<br />
#>       a  b  c<br />
#> 1     1  5  5<br />
#> 2     2  4  3<br />
#> 3     3  3  1<br />
#> 4     4  2  4<br />
#> 5     5  1  1<br />
#> NA   NA NA NA<br />
#> NA.1 NA NA NA<br />
subset2(sample_df, a == condition)<br />
#> Error: object 'a' not found<br />
subset2(sample_df, a == condition_call)<br />
#> Warning: longer object length is not a multiple of shorter object length<br />
#> [1] a b c<br />
#> <0 rows> (or 0-length row.names)

我们很多时候都会这样的用，因为我们不一定都会用a>=4当条件，更多时候我们希望用个变量就像例子中的那样，可问题出现了，结果和我们想的完全不一样，对于这个问题的理解，关键在于深入理解

eval(exp,env,enclos)这个函数，其中exp的符号首先在env中找，找不到就去它的外部环境中找，对于这里是subset2函数内部环境，找到就不继续往外找了，所以这里的y没找到就继续从外找，x找到了参数x,condition找到了参数，condition_call是局部变量

下面一个个解释：

y=4不用解释了，a==x的结果其实是

a=1:5

sample_df[a==sample_df, ]

而a==sample_df大家实验下结合subsetting中说的logical下标的情况，就会知道

R>a==sample_df<br />
        a     b     c<br />
[1,] TRUE FALSE FALSE<br />
[2,] TRUE FALSE FALSE<br />
[3,] TRUE  TRUE FALSE<br />
[4,] TRUE FALSE  TRUE<br />
[5,] TRUE FALSE FALSE

再取下标的时候只取TRUE的，这样长度就是7，前5个可以正确取出，后2个越界了返回NA，我们之前也说过详细的subset出现问题的各种情况，这里[]越界返回NA只是比较简单的一种情况，不过这里理解够用了

对于a==condition,它会从参数中找到condition的值是a==condition,然后这个时候它无法找到a的值了，

所以a==(a==condition)第一个a是x中找到的，而第2个就得从函数外面找，结果没找到

对于最后一个a==condition_call我暂时没理解结果，我开了个[a href=http://cos.name/cn/topic/157299]讨论帖，大家可以看看有没有什么解答[/a]

分析了这么多，解决方案就是一种，就是让eval(exp,x,enclos)中的exp中的符号在env中找不到的时候去一个指定的地方找，而不是从eval所处的环境中，这个参数就是enclos,如果再找不到，就从enclos的parent.env中找。。。

所以我们指定一个ce环境

subset2 <- function(x, condition) {<br />
  condition_call <- substitute(condition)<br />
  r <- eval(condition_call, x, parent.frame())<br />
  x[r, ]<br />
}</p>
<p>x <- 4<br />
subset2(sample_df, a == x)<br />
#>   a b c<br />
#> 4 4 2 4

这样，就不必纠结从参数，从局部变量中找到的结果了。。。。

还有个方式就是

subset2a <- function(x, condition) {<br />
  condition_call <- substitute(condition)<br />
  env <- list2env(x, parent = parent.frame())<br />
  r <- eval(condition_call, env)<br />
  x[r, ]<br />
}</p>
<p>x <- 5<br />
subset2a(sample_df, a == x)

所以这里又一个总结就出来了，这个总结非常重要：

对于这个函数eval(exp,env,enclos)，如果env不是一个real env而是list,df没有父环境的，我们exp一旦没找到值，如果没有指定enclos参数，那就会从当前环境中往外找，一旦指定就从那里开始找，而对于env是一个real env,就不会去看enclos，直接从env的父环境中找

上面的总结非常重要，基本上是对eval工作所有情况的覆盖

好了，下面我们来看几个例子（ADR的几个习题）来回顾下：

第一个例子，为什么不能在subset2用quote替代substitute？

subset2 <- function(x, condition) {<br />
  condition_call <- quote(condition)<br />
  r <- eval(condition_call, x)<br />
  x[r, ]<br />
}

因为，eval(quote(condition),x)怎么找condtion? 先从x中找，而x是一个df,又没有指定enclos所以从subset2内部找condition,就找到了参数的condtion,所以显然结果又是找不到'a'!!

第二个例子，我们想subset(mtcars, , -cyl) ， subset(mtcars, , disp:drat)这样来取值

select <- function(df, vars) {<br />
  vars <- substitute(vars)<br />
  var_pos <- setNames(as.list(seq_along(df)), names(df))<br />
  pos <- eval(vars, var_pos)<br />
  df[, pos, drop = FALSE]<br />
}<br />
select(mtcars, -cyl)

我只能说这段代码太赞了，它用list存储表达式的列位置，然后eval出来，然后根据这个结果取子集

显然这个功能也是因为 df[ , -1] ,df[ ,2:3]之类的可以实现，plyr这个包好像有专门这样做的函数，那里做的

好像高能更高级

第三个例子是实现transform

df <- data.frame(x = 1:5)<br />
transform(df, x2 = x * x, x3 = x2 * x)<br />
plyr::mutate(df, x2 = x * x, x3 = x2 * x)

我就不贴transform,mutateR代码了，大家可以直接在交互界面查看，然后仔细研究，大概思路是

产生x2,x3然后对比原先的df，重名就替换，不重名就加入

至于mutate可以利用刚刚创建的列，是因为它是通过一个循环来修改df,每循环一次改一次df,transform则是一次性eval好再找，所以不能利用刚创建的列

好了，下面就来讲最难理解的部分了，Calling from another function，当然如果你充分理解了我之前提到的几个要注意的点，也许没那么难理解

subset2 <- function(x, condition) {<br />
  condition_call <- substitute(condition)<br />
  r <- eval(condition_call, x, parent.frame())<br />
  x[r, ]<br />
}</p>
<p>scramble <- function(x) x[sample(nrow(x)), ]</p>
<p>subscramble <- function(x, condition) {<br />
  scramble(subset2(x, condition))<br />
}

这是一个随机打乱顺序抽取提取下标的子集的函数，2个函数我们应该非常熟悉，合在一起也没什么，不过就是不工作！

为什么。。。我试着读作者写的原因，理解了几天，终于大概懂了。。。[s:12]

作者用debug技术弄的，我说过我不懂，我很少debug....所以我直接说原因：

在这一段subscramble <- function(x, condition) {

scramble(subset2(x, condition))

}中，subset2里的condition_call <- substitute(condition)结果是什么？？这个非常重要，如果没有外层函数，直接用条件a>=4那就是a>=4这个表达式，而这里呢？我们会想外层函数把condition的值传给内层啊，所以也是a>=4这个表达式啊！！难道不是？确实不是，回到我开这个帖子的时候说的那句话：

It looks at a function argument and instead of seeing the value, it sees the code used to compute the value

就是说，它只输出用来产生参数值的那个表达式！所以这里condition的值是a>=4，但是condition用来产生a>=4,a>=4用来产生结果，所以condition_call <- substitute(condition)结果是condition表达式.

理解这个非常关键，然后就进入了eval(condition_call, x, parent.frame())，这个函数我们前面弄的很清楚了，

我们要从x中找condition这个符号的值，显然没找到，那我们就去enclos指定的parent.frame()中找，而它的parent.frame()是什么呢？是subscramble函数，所以从哪里找到condition呢？从subscramble的参数中找到，OK！

值是a>=4，所以我们就要eval(a>=4)这个东西，但是，从subset2往外直到subscramble的ee都没有叫a的对象，所以失败！

所以这个问题涉及到的，对eval(exp,env,enclos)的深刻理解，这个之前总结过了，以及对substitute的深刻理解，开头那句话非常关键，加起来就得出结论了，与此同时，如果我

a <- 4<br />
subscramble(sample_df, a == 4)<br />
#>   a b c<br />
#> 4 4 2 4<br />
#> 5 5 1 1<br />
#> 1 1 5 5<br />
#> 2 2 4 3<br />
#> 3 3 3 1</p>
<p>a <- c(1, 1, 4, 4, 4, 4)<br />
subscramble(sample_df, a >= 4)<br />
#>     a  b  c<br />
#> NA NA NA NA<br />
#> 4   4  2  4<br />
#> 5   5  1  1<br />
#> 3   3  3  1

为什么会出现这些结果也就好理解多了，它们都是在global.env()找到，这个也是subscramble的ee中找到，所以用来a>=4，第一个返回logical取下标会自动扩展，第二个越界

最大的问题理解了！我们来看解决方案

subset2_q <- function(x, condition) {<br />
  r <- eval(condition, x, parent.frame())<br />
  x[r, ]<br />
}<br />
subset2 <- function(x, condition) {<br />
  subset2_q(x, substitute(condition))<br />
}</p>
<p>subscramble <- function(x, condition) {<br />
  condition <- substitute(condition)<br />
  scramble(subset2_q(x, condition))<br />
}<br />

这样就确保传给subset_q的是一个表达式，而subset_q直接默认用表达式，而不用NSE的方式

所以escape hatch的意思是：

As a developer, you should always provide an escape hatch: an alternative version of the function that uses standard evaluation

这是原话，大概是提供一个非NSE版本即SE版本

好了，最大的问题解决了，我们再来深入substitute的全面理解

a <- 1<br />
b <- 2<br />
substitute(a + b + z)<br />
#> a + b + z<br />
f <- function() {<br />
  a <- 1<br />
  b <- 2<br />
  substitute(a + b + z)<br />
}<br />
f()<br />
#> 1 + 2 + z

可以看到在全局和函数内部行为的不同！作者用pryr ::subs()替代使得在全局也一样，具体就不介绍了，大家可以当substitute用，只不过全局行为也一致了

我们来看下substitute的规则：

1.an ordinary variable, it’s replaced by the value of the variable.

2.a promise (a function argument), it’s replaced by the expression associated with the promise.

3...., it’s replaced by the contents of ...

这些规则是substitute的核心东西：

x <- quote(mpg)<br />
y <- quote(disp)<br />
subs(xyplot(x ~ y, data = mtcars))<br />
xyplot3 <- function(x, y, ...) {<br />
  substitute(xyplot(x ~ y, ...))<br />
}<br />
xyplot3(mpg, disp, data = mtcars, col = "red", aspect = "xy")<br />
#> xyplot(mpg ~ disp, data = mtcars, col = "red", aspect = "xy")

最后我们要注意substitue,

x <- quote(a + b)<br />
substitute(x, list(a = 1, b = 2))<br />
#> x

这和eval不一样，所以为了实现这个功能我们这样

substitute_q <- function(x, env) {<br />
  call <- substitute(substitute(y, env), list(y = x))<br />
  eval(call)<br />
}</p>
<p>x <- quote(a + b)<br />
substitute_q(x, list(a = 1, b = 2))<br />
#> 1 + 2

这里其实就是类似 eval(quote(quote(1+2)))的思想，只不过substitute更复杂，外层让y变成quote(a+b),内层让它变成quote(quote(1+2))，要注意，substitute本身也返回表达式

所以这里又来个小总结，substitute和quote的区别，quote就是返回输入，而substitute遵循上面3个规则，以及注意substitute(exp,env)和eval(exp,env)的区别，对于第一个参数都是exp的情况行为比较一致，但如果是变量，就要注意上面介绍的前者需要一个转换，就是substitute_q

最后来个神奇的东西

subs(a + b, list("+" = quote(f)))<br />
#> f(a, b)<br />
subs(a + b, list("+" = quote(<code>*</code>)))<br />
#> a * b

最后的最后来一些扫尾的知识，NSE的缺点是不referentially transparent.，也就是对于NSE实现的f(x),f(y),f(10)的结果可能不一样，当然x,y的值都为10，这点我们前面已经很清楚了，适合用的地方就是

library(ggplot2)<br />
library("ggplot2")<br />
ggplot2 <- "plyr"<br />
library(ggplot2)

这里主要用了

g <- function(x) deparse(substitute(x))<br />
g(1:10)<br />
#> [1] "1:10"<br />
g(x)<br />
#> [1] "x"<br />
g(x + y^2)<br />
#> [1] "x + y^2"

还有就是

x <- 10<br />
y <- "a"<br />
df <- data.frame(x, y)<br />
names(df)<br />
#> [1] "x" "y"

如果你仔细读了前面的，我觉得这个应该可以轻松实现[s:11]

PS;写了3个多小时终于写完这个非常复杂，精细的NSE了

思路总结:

1.quote vs substitute

2.substitute的3个规则等

3.substitute(exp,env）vs eval(exp ,env)

4.eval(x) vs eval(exp,env,enclos)

5.案例问题，特别是函数调用NSE函数的问题，大部分都是exp到底是什么，怎么来的，eval(exp,env,enclos)的寻值规则是什么

不过Metaprogramming的东西还有两章呢，下章就是expression，继续加油！

superdesolator

80.ADR之Metaprogramming之expressions

前面介绍了NSE，也提到了expression,quote()返回的就是一个expr(简写下同），到底什么是expr呢，expr就是一个可以被R执行的动作，而quote()就可以捕捉这个动作，expr也被叫做abstract synatax tree(ast)抽象语法树，因为它也代表了R代码的结构层次树，所以你既可以把它理解为静态的代码结构层次树，也可以理解为动态的可以被R执行的动作.我们可以用pryr::ast来看，注意这个ast函数非常有用！后面经常用到！

<br />
#我们通过ast来看代码的结构树<br />
ast(y<-x*10)<br />
#或者用quote把这个动作捕捉下来<br />
quote(y<-x*10)<br />

而正如很多时候R语言的不一致性带来的困扰，expression()却不是返回一个expr,而是一个list of exprs,关于quote()和expression()更多不同请看这个帖子.

正如我所说，ast是非常好的东西，我们可以用它看清楚很多东西：

ast(if(x>1) x else 1/x)<br />
#或者也许你想看看 if else if else 好的<br />
ast(if (x>3) 1 else if(x>2) 2 else 3)<br />
#感觉到这个函数牛x的你 也许会试验while<br />
ast( while( i>1) {a<-a+1})<br />
#同样你也可以看清楚for,其实一切控制在R中都是函数<br />

我们可以用ast干什么呢？既然都可以看到代码的结构了，还有什么不能做啊，比如看清楚if else的原理之后，结合if是个函数，基础再好点，写出下面不难

#如果前面你基础好，下面对你应该很熟悉<br />
"if"(3>1,2,3)

那么我也想写个for的函数，OK，用ast看下

<br />
ast( for( i in 1:3) {print(4)})<br />

我们发现函数name是for,第一个参数是i,第二个是一个1:3返回的向量，第三个是循环体

那我们写个"for"

"for"(i,1:3,print(4))

怎么样？ast神奇不！

熟悉了ast之后，我们再回到主线，讨论expr,通常一个expr可能有4种形式：constants,names,calls,pairlists.

下面我们结合ast来先大概看下这4个部分到底都是怎么样的:

1.constants

constants就是长度为1的atomic vectors，ast原样返回

ast("a")<br />
ast(1L)

2.names

就是通常对象的名字，ast返回名字前面加个backtik（反引号）

ast(mean)

3.calls

代表了调用一个函数的动作，它是递归的

ast(f(g(),h(1,a)))

4.pairlists

这是一个R的遗留的东西，很多地方被list取代了，除了函数参数中,ast会用[]标出

ast(function(x=1,y) x)

大致了解了要说的内容之后，先介绍几个有趣的知识，再来详细深入names,calls,pairlists(constants就u需要了吧）

第一个有趣的知识是

ast(x+y%+%z)<br />
ast(x^y%+%z)

这样我们就可以看出%+%之类的函数的优先权了

第二个有趣的知识见这里

就是substitute会插入calltree引起的一些问题

好，下面深入

1.names

对于string,我们可以用as.name转成名字,对于name我们用is.name判断,稍微注意下a a这样的名字，用backtick反引一下

identical(quote(name),as.name("name"))

as.name的用途就在于当函数传的是字符的时候，我们需要处理字符，否则我们直接用quote更省事

假如我们要实现自己版本的get,assign,我们怎么做呢？我们用到的工具有as.name,substitute,eval，至于环境的问题，我们不作深入

get2=function(x){<br />
  eval(as.name(x),parent.frame())<br />
}<br />
assign2=function(x,y){<br />
x<-as.name(x)<br />
tmp<-substitute(x<-y)<br />
eval(tmp,parent.frame())<br />
}<br />
#下面实验<br />
R>a=1<br />
R>get2("a")<br />
[1] 1<br />
R>yyyy<br />
Error: object 'yyyy' not found<br />
R>assign2("yyyy",2:4)<br />
R>yyyy<br />
[1] 2 3 4<br />

因为我们必须处理string,所以使用as.name,我们无法使用quote,还要注意的是assign2这个函数，其实我想了蛮久才写出来，最关键是对substitute本身以及3个规则的理解，用x它就会在当前环境中找，并用as.name(x)的结果替换，这点非常重要因为直接

assign2=function(x,y){<br />
tmp<-substitute(as.name(x)<-y)<br />
eval(tmp,parent.frame())<br />
}

是不行的！我觉得虽然代码很短，但是道理蛮深的，首先必须quote(a<-4)这个形式才能eval,而quote本身不能把a<-y，用4替换y,把这个形式变出来，因为它只会原样返回它的输入，所以得用substitute,而要完成类似substitute(a<-4)，而不是substitute(as.name("a")<-4)因为我们不能as.name("a")<-4，因为as.name("name")和quote(name)效果一样，所以as.name("a")<-4相当于quote(name)<-4和我们原意不一样，所以我们用x=as.name("a")，然后substitute(x<-4)这样x就被as.name("a")的结果a替换掉了

对于names,除了上面最基本的应用，有个比较特别的名字，empty name,它用来代表missing arguments,但它不能被绑定到一个名字

f<-function(x) 10<br />
formals(f)$x<br />
is.name(formals(f)$x)<br />
as.character(formals(f)$x)<br />
a<-formals(f)$x

我们可以用formals(f)去修改f的参数

g<-function(x=20,y=x) {x+y}<br />
formals(g)

对于这个函数我们看到g的参数情况，如果我们想把y设置为没有默认的，那就相当于给个empty name给它，方式为

formals(g)$y<-quote(expr =)<br />
g

quote(expr =)是显示创建empty name的方法，这样的话，没有默认值到有默认值，有默认值到没有默认值，修改默认值都可以做了，一个函数的参数就可以被我们改动了，甚至我们可以增减：

formals(g)$z<-quote(expr=)<br />
g<br />
#function (x = 20, y = x, z)<br />
#{<br />
 #   x + y<br />
#}<br />
formals(g)$x<-NULL<br />
g<br />
#function (y = x, z)<br />
#{<br />
 #   x + y<br />
#}

2.calls

这部分是重点，一个call很像一个list,它有length,[[,[方法，并且是递归的，因为参数中可能也是一个call等

通常一个call的第一个部分是函数名或者另一个call

x<-quote(mean(1:10))<br />
x[[1]]<br />
x[[2]]<br />
x<-quote(add(10)(20))<br />
x[[1]]

其它部分就是参数了，它们可以通过位置或者名字取到，也可以通过$<-,[[<-来修改，这部分我就不写了，大家可以增加或者减少（NULL），自己试验下，和names中介绍的修改参数差不多

当然我们也可以通过[方法，x[-1]来修改，然后as.list(x[-1])这样就把参数显示转换为一个列表了

另外需要注意的是，函数通过位置修改要特别注意，这个位置到底是什么参数，因为R的关键字参数的概念，我们可以用f(x=3,4)之类的形式调用一个函数，这和f(4,x=3)效果是一样的，但是要注意这两个call的参数list是不同的，我们得小心的对付通过位置取参数的方式，如果你觉得很麻烦，作者介绍了standardise_call()，它内部调用match.call()来完成参数的标准化，具体看下例子

m1<-quote(read.delim("data.txt",sep="|"))<br />
m2<-quote(read.delim(s="|","data.txt"))<br />
standardise_call(m1)<br />
standardise_call(m2)

我们再看下standardise_call的源码

standardise_call<br />
function (call, env = parent.frame())<br />
{<br />
    stopifnot(is.call(call))<br />
    f <- eval(call[[1]], env)<br />
    if (is.primitive(f))<br />
        return(call)<br />
    match.call(f, call)<br />
}<br />
<environment: namespace:pryr>

基本上就是取得函数名字的值，也就是函数对象，然后进行match.call(f,call)这个match.call()就帮我们标准化了，包括参数部分匹配等等问题全都自己解决了，至于match.call()的源码是什么，暂时先不管了

当然我们quote()可以返回一个call,我们也可以自己创建call,方法有两个：

第一，利用call()，它的第一个参数是函数名字，其它都是exprs代表的参数：

call("mean",quote(1:10),na.rm=TRUE)

这里有个问题，其它参数应该是表达式，所以下面是不同的：

<br />
a<-call("mean",1:10)<br />
b<-call("mean",quote(1:10))<br />
identical(a,b)<br />

关键就在于参数是不是expression了，第一个不是，第二个是，那么为什么第一个会形成一个call呢？不是说其它参数应该是expression吗？而na.rm=TRUE又是什么？看这里的讨论

第二，利用as.call（）将一个list变成一个call,list第一个参数要求是一个name或者一个call,其他的是参数

as.call(list(quote(mean),quote(1:10)))<br />
as.call(list(quote(add(10)),20))

知道了基本的call内容之后，我们来看看一些实例，这些实例都是ADR这部分的题目中选出的，基本上是我做出来供大家参考的，不会的没办法呈现给大家啦[s:11]

第一个问题：实现自己的do.call

do.call2=function(i,...){<br />
  tmp= as.call(list(as.name(i),...))<br />
  eval(tmp,parent.frame())<br />
}

第二个问题：实现如下concat()函数

concat(quote(f),a=1,b=quote(mean(a)))<br />
#f(a=1,b=mean(a))

concat<-function(fname,...){<br />
  as.call(list(fname,...))<br />
}

第三个问题：利用list()传参，如下

<br />
make_call(quote(mean),list(quote(x),na.rm=TRUE))<br />
#mean(x,na.rm=TRUE)<br />
make_call(quote(mean),quote(x),na.rm=TRUE)<br />
#mean(x,na.rm=TRUE)

我自己把这个问题升级了，就是第一种方式和第二种混合也可以也就是说

make_call(quote(mean),list(quote(x),na.rm=TRUE),a=1,4,list(c=3))<br />
#mean(x, na.rm = TRUE, a = 1, 4, 3)

为了实现这个功能，我是花了蛮久的时间代码如下：

<br />
make_call<-function(fname,...){<br />
  tmp=list(...)<br />
  nametmp=names(tmp)<br />
  tmp=lapply(tmp,function(x){<br />
    if(!is.list(x))<br />
      list(x)<br />
    else<br />
      x<br />
  }<br />
  )<br />
  for( i in seq_along(tmp)) {<br />
    names(tmp[[i]])[1]=nametmp[i]<br />
  }<br />
  transtmp=Reduce('c',tmp)<br />
  as.call(c(fname,transtmp))</p>
<p>}

实现上述代码着实花了我不少功夫，最大的技巧是利用c这个函数去合并list,这个细节非常有用，这在我ADR系列的第一个帖子中就有提到，这也是transtmp=Reduce('c',tmp)代码的由来，它遍历包含list的list,tmp,然后把那些list给合并在一起，然后怎么得到tmp呢，就是对参数list(...)这样得到的是含list和其他结构的list,利用lapply把非list转化为list就可以了，事情的框架基本就是这两步，其余的代码干的事情是，光靠上面两步会出现一个a=4这种参数会被弄成f(...,4,)名字a会丢，原因就是在lapply完之后，名字并没有给大list的对应位置的等待合并的小list,那个循环就是用来给这样的名字的

通过这3个问题，我们都会发现基本上解决思路都是化为list(quote(fname),....)然后转成call的形式！

在继续深入call之前，我们来轻松一下，看下mode这个函数的实现

function (x)<br />
{<br />
    if (is.expression(x))<br />
        return("expression")<br />
    if (is.call(x))<br />
        return(switch(deparse(x[[1L]])[1L], <code>(</code> = "(", "call"))<br />
    if (is.name(x))<br />
        "name"<br />
    else switch(tx <- typeof(x), double = , integer = "numeric",<br />
        closure = , builtin = , special = "function", tx)<br />
}<br />
<bytecode: 0x082049d8><br />
<environment: namespace:base>

class和typeof的关系前面的帖子已经说清楚了，当时提到了mode,现在正式面对这个问题了，我们这里的概念，name,typeof出来是symbol,call，typeof出来是language,不过不用管了，看以看到mode是基于is.call,is.name来判断的，包括expression,然后就是switch(typeof(x),....)如果是double,integer返回numeric,如果是closure,builtin,special返回function,其它就是不动的返回

这里的closure就是我们写的函数，builtin是primitive函数，special是不计算参数的函数，它们都会返回"function"

typeof(expression)<br />
typeof(quote)

这样的mode有的时候对于初学者很有用，试想一开始typeof(x)的时候就出现closure,builtin,special,谁能想到它们都是function呢？好了，关于这方面知识，知道到这里暂时够了，以后深入再说吧[s:11]

好，下面来继续说call,我们要讨论的是捕捉当前的call.也就是在函数内部捕捉调用这个函数的表达式：

f <- function(abc = 1, def = 2, ghi = 3) {<br />
  list(sys = sys.call(), match = match.call())<br />
}<br />
f(d = 2, 2)<br />
#> $sys<br />
#> f(d = 2, 2)<br />
#><br />
#> $match<br />
#> f(abc = 2, def = 2)

match.call()经常用于建模的函数，用于更新模型：

mod <- lm(mpg ~ wt, data = mtcars)<br />
update(mod, formula = . ~ . + cyl)<br />
#><br />
#> Call:<br />
#> lm(formula = mpg ~ wt + cyl, data = mtcars)<br />
#><br />
#> Coefficients:<br />
#> (Intercept)           wt          cyl<br />
#>       39.69        -3.19        -1.51

update就用来更新了模型，这是怎么做到的？我们来写出来

<br />
update_call <- function (object, formula., ...) {<br />
  call <- object$call</p>
<p>  # Use update.formula to deal with formulas like . ~ .<br />
  if (!missing(formula.)) {<br />
    call$formula <- update.formula(formula(object), formula.)<br />
  }</p>
<p>  modify_call(call, dots(...))<br />
}<br />
update_model <- function(object, formula., ...) {<br />
  call <- update_call(object, formula., ...)<br />
  eval(call, parent.frame())<br />
}<br />
update_model(mod, formula = . ~ . + cyl)<br />
#><br />
#> Call:<br />
#> lm(formula = mpg ~ wt + cyl, data = mtcars)<br />
#><br />
#> Coefficients:<br />
#> (Intercept)           wt          cyl<br />
#>       39.69        -3.19        -1.51

其实它的本质就是提取object的call，然后利用call的性质，用call$formula去修改formula这个参数，其中又利用了update.formula,这个函数是用C写的，作者一道练习是让我们用R写出来，我只能说，天啊，我还不了解formula这个类型，怎么能写[s:12]？

call的formula参数更新好后，就需要更新来自...的参数了，这个时候我们需要用到dots()来把...转成输入内容的list,然后又调用modify_call去修改，基于篇幅(呵呵[s:11]已经相当长了）我就不列出这两个函数的源码了，大家可以自己pryr::dots,主要的思路就是用dots获得...的输入内容的list，然后在modify_call里面，然后取出list的名字，在call本身的参数名字里写个循环，如果list的名字出现就修改call,用基本的$<-修改，modify_call还是值得一读的

这样一来的话，新的call的包括主要的formula在内的所有参数都被更新了，然后在update_model里eval一下就好了

最后值得注意一下的是，env的问题，我们自己的是在parent.frame里执行，而R原来的版本update会在globalenv()中执行，所以对下面的例子：

f <- function() {<br />
  n <- 3<br />
  lm(mpg ~ poly(wt, n), data = mtcars)<br />
}<br />
mod <- f()<br />
update(mod, data = mtcars)<br />
#> Error: object 'n' not found<br />
update_model <- function(object, formula., ...) {<br />
  call <- update_call(object, formula., ...)<br />
  eval(call, environment(formula(object)))<br />
}

就好了，本来的update是找不到n的，因为n在f的exe e中，而environment(formula())就记录了这个环境,formula我觉得后面会遇到，到时候会详细介绍的

3.pairlists

前面说过,pairlists是R过去的遗产，目前只保留在函数参数里，它们的表现和list很像，但底层是用链表而不是向量，原话如下：

Pairlists are a holdover from R’s past. They behave identically to lists, but have a different internal representation (as a linked list rather than a vector). Pairlists have been replaced by lists everywhere except in function arguments.

然后下面我们要做非常有趣的事情，创造函数，主要利用"function",给它参数，函数体然后在env，eval就可以了

make_function <- function(args, body, env = parent.frame()) {<br />
  args <- as.pairlist(args)</p>
<p>  eval(call("function", args, body), env)<br />
}

通常我们给的参数是alist()它不计算它的参数，alist(x=a)是list(x=quote(a))的等效简写形式，配合alist我们

add <- make_function(alist(a = 1, b = 2), quote(a + b))<br />
add(1)<br />
#> [1] 3<br />
add(1, 2)<br />
#> [1] 3</p>
<p># To have an argument with no default, you need an explicit =<br />
make_function(alist(a = , b = a), quote(a + b))<br />
#> function (a, b = a)<br />
#> a + b<br />
# To take <code>...</code> as an argument put it on the LHS of =<br />
make_function(alist(a = , b = , ... =), quote(a + b))<br />
#> function (a, b, ...)<br />
#> a + b

可以看到传特定参数的时候要注意的都写在上面了，它在闭包的应用也很方便，因为它可以更好的看函数体

adder <- function(x) {<br />
  make_function(alist(y =), substitute({x + y}), parent.frame())<br />
}<br />
adder(10)<br />
#> function (y)<br />
#> {<br />
#>     10 + y<br />
#> }

这里稍微提下alist(a)和alist(a=)是有区别的，一个是某个位置是参数的值为a的name，一个是名字为a的参数的值为emptyname.

我们来看下make_function的应用之一partial函数,大家可以pryr::partial,这里不复制了，仔细阅读完源码，我们发现，它做的事情是得到一个指定参数以及加入...的call表达式，然后创建一个函数，参数是...，源码中用的是

list(...=expr=),我觉得也可以用我们这里提到的alist(...=),然后函数体就是那个call表达式，在加入环境以及一些细节的控制参数，一个partical就得到了，这个函数的思路非常值得学习，下面再看下unenclose，它的作用是在函数的ee环境中找到函数定义需要的值，它的源码我也不贴了，大家pryr::unenclose就可以了，在读这个函数之前，需要回忆下substitute_q的作用，它的作用就是传入的是quote，然后里面用双层substitute来产生结果，大家可以看下这个函数的源码，回到unenclose,其实就简单了，利用environment(f)返回f的ee,然后利用formals(f)返回参数，然后substitute_q(body(f),env)就得到替换值过的函数了，这个源码也值得推荐

最后我们来看下make_function的一个应用

curve(sin(exp(4 * x)), n = 1000)

这里的x是一个pronoun,这种函数叫anaphoric，在其他语言，Perl之类有应用.我们来看实现：

curve2 <- function(expr, xlim = c(0, 1), n = 100, env = parent.frame()) {<br />
  f <- make_function(alist(x = ), substitute(expr), env)</p>
<p>  x <- seq(xlim[1], xlim[2], length = n)<br />
  y <- f(x)</p>
<p>  plot(x, y, type = "l", ylab = deparse(substitute(expr)))<br />
}

其实很简单就是利用make_function,和x,以及函数表达式创造了一个函数，作者还列出了一种方法，利用substitute加env来找x的值

curve3 <- function(expr, xlim = c(0, 1), n = 100,<br />
                   env = parent.frame()) {<br />
  env2 <- new.env(parent = env)<br />
  env2$x <- seq(xlim[1], xlim[2], length = n)</p>
<p>  y <- eval(substitute(expr), env2)<br />
  plot(env2$x, y, type = "l", ylab = deparse(substitute(expr)))<br />
}

我们前面提到过deparse的作用，parse则相反的功能，它把字符串变成expression

exp <- parse(text = c("<br />
  x <- 4<br />
  x<br />
  5<br />
"))<br />
length(exp)<br />
#> [1] 3<br />
typeof(exp)<br />
#> [1] "expression"</p>
<p>exp[[1]]<br />
#> x <- 4<br />
exp[[2]]

OK,paste返回一个expression里面包含各种具体的expression类型call,name,...

能把字符串处理成这样，那显然我们想到了source,从字符串文件得到并运行源码，我们来实现自己的source

<br />
simple_source <- function(file, envir = new.env()) {<br />
  stopifnot(file.exists(file))<br />
  stopifnot(is.environment(envir))</p>
<p>  lines <- readLines(file, warn = FALSE)<br />
  exprs <- parse(text = lines)</p>
<p>  n <- length(exprs)<br />
  if (n == 0L) return(invisible())</p>
<p>  for (i in seq_len(n - 1)) {<br />
    eval(exprs[i], envir)<br />
  }<br />
  invisible(eval(exprs[n], envir))<br />
}<br />

这段代码最主要的就是exprs<-parse(text=lines),然后遍历exprs，在一个新环境中eval它们，最后的隐式返回大家可以不必太在意

最后的最后，我们来介绍遍历AST的内容：

基本的都介绍了，我们可以用substitute,modify_call之类的函数以及相应的知识来修改call,但是如果要更复杂的功能，我们必须学着遍历ast,codetools包里有个findGlobals()可以寻找一个函数的依赖情况

f<-function(x=a,y=2){<br />
  x+y<br />
}<br />
codetools::findGlobals(f)

对于我们要做的遍历ast树，基本上框架如下

recurse_call <- function(x) {<br />
  if (is.atomic(x)) {<br />
    # Return a value<br />
  } else if (is.name(x)) {<br />
    # Return a value<br />
  } else if (is.call(x)) {<br />
    # Call recurse_call recursively<br />
  } else if (is.pairlist(x)) {<br />
    # Call recurse_call recursively<br />
  } else {<br />
    # User supplied incorrect input<br />
    stop("Don't know how to handle type ", typeof(x),<br />
      call. = FALSE)<br />
  }<br />
}

这部分挺难的，我们来仔细理解，从简单的事情做，首先第一个任务是寻找一个expression中是否存在T,F之类的逻辑简写

ast(TRUE)<br />
ast(T)

通过这个观察发现，如果是name有可能，如果是atomic不可能含有T,F所以代码如下

logical_abbr <- function(x) {<br />
  if (is.atomic(x)) {<br />
    FALSE<br />
  } else if (is.name(x)) {<br />
    identical(x, quote(T)) || identical(x, quote(F))<br />
  } else if (is.call(x) || is.pairlist(x)) {<br />
    for (i in seq_along(x)) {<br />
      if (logical_abbr(x[[i]])) return(TRUE)<br />
    }<br />
  } else {<br />
    stop("Don't know how to handle type ", typeof(x),<br />
      call. = FALSE)<br />
  }<br />
}<br />
logical_abbr(quote(mean(x, na.rm = T)))<br />
#> [1] TRUE

看上去挺简单，但是要真的理解我觉得还是蛮难的，关键在于x是call,pairlist那里的循环递归的理解，循环第一层次，遇到一个call深入遍历进去，一有T,F就返回，没有就继续第一层，如果深入一个call的时候又遇到一个call就又深入，然后遍历玩这个call,返回上个call的其余部分遍历，没有就跳出call的遍历，继续第一层遍历，一次类推。。。你被绕晕了没。。。[s:11]

第二个任务是找出所有被<-创建的变量,作者修改了这个函数好几次，每次带入一个问题，不断解决得到最后的版本，大家可以去原书看，太长了就不贴了

大家可以对着原书，看我的思路，思路是：第一遍当是call的时候判断name是不是'<-'是就返回call第2个参数，否则就继续遍历这个call,这样做的问题是输出结果不漂亮，于是第二遍都返回character,这样的结果会出现重复的名字，以及<-含有<-的情况就没有被考虑，于是第三遍用了unique，以及把继续遍历call的结构放到判断<-的外面，就是即使是<-也要遍历它，这些遍历可以进行的关键是，lapply一个call的时候，第一个已经是name了不会是call了所以不会死循环（我开始读的时候就纠结这里），最后针对names(a)<-1会出现问题，作者使用is.name()来确保返回的不是一个call的函数名

很复杂吧，慢慢熬吧[s:11]

作者还利用上面的框架实现了自己的bquote2，基本上就是is.call为TRUE的时候看是不是.()函数是的话eval第2个参数，不是就继续递归a<-lapply(x,bquote2,env)并要as.call(a),对于pairlist则直接递归，具体的大家自己去研究吧

PS:当我做这个练习的时候发现了作者的一个BUG

f<-function(x=TRUE) {<br />
+ g(x+T)<br />
+ }<br />
#try let this work<br />
logical_abbr(f)

我的解决方法

logical_abbr(call("function",formals(f),body(f)))

本来应该是对的，结果老错，我就以为是我自己写错了，结果后来发现竟然是作者些的logical_abbr函数有问题！

logical_abbr <- function(x) {<br />
  if (is.atomic(x)) {<br />
    FALSE<br />
  } else if (is.name(x)) {<br />
    identical(x, quote(T)) || identical(x, quote(F))<br />
  } else if (is.call(x) || is.pairlist(x)) {<br />
    for (i in seq_along(x)) {<br />
      if (logical_abbr(x[[i]])) return(TRUE)<br />
    }<br />
    FALSE   #应该在这里加个FALSE ！！<br />
  } else {<br />
    stop("Don't know how to handle type ", typeof(x),<br />
         call. = FALSE)<br />
  }<br />
}

superdesolator

81.ADR之Metaprogramming之Domain specific languages(DSLS)

对于这章来说,纯学R的人应该不太会用到,所以我不会写的很详细，但读读总不错，如果你要用到对应的知识，具体的去看原书吧

这一章刻画了DSLS,也就是把R的code转成其它语言，作者举了个dplyr包的，将Rcode转成SQL的函数translate_sql()

然后作者实现了转成两种语言的任务！

下面我们就分别一步步介绍作者怎么转的：

第一，HTML

关于HTML tag的意思，实体等意思我就不介绍了，直接介绍作者思路：

作者首先创建一个html类的构造函数，然后定义S3的print方法，大家可以不必细推怎么print出来

html <- function(x) structure(x, class = "html")<br />
print.html <- function(x, ...) {<br />
  out <- paste0("<HTML> ", x)<br />
  cat(paste(strwrap(out), collapse = "\n"), "\n", sep = "")<br />
}

那么接下来作者首先做的事情就是把<,>之类的符号转成实体定义&lt,&gt

escape <- function(x) UseMethod("escape")<br />
escape.html <- function(x) x<br />
escape.character <- function(x) {<br />
  x <- gsub("&", "&", x)<br />
  x <- gsub("<", "<", x)<br />
  x <- gsub(">", ">", x)</p>
<p>  html(x)<br />
}<br />
escape.list <- function(x) {<br />
  lapply(x, escape)<br />
}

这个思路就是如果generic函数输入是html类就直接返回，如果是字符就替换，如果是list就对每个部分escape一下

这样的话，就不必担心实体定义的问题了，这个技术大概叫escape，还记得前面的escape hatch吗？这类思想就是帮你直接处理好一些琐碎的问题

下面就是基本的tag函数了，就是说，我输入

p("some text",b("some bold text"),class="mypara")

就要输出

<br />
<p class="mypara"> some text <b> some bold text </b> </p>

这种HTML语言的形式，依次类推更复杂的情况

因为tags有很多很多，我们不可能全部写出来，所以方法就是尝试写一个tag比如p,然后用closure的思想包装一下

对于p tag我们看到它的输入有位置参数有命名参数，位置参数都是内容，命名参数都是标签的属性，所以我们需要对一个list提取命名的和未命名的

<br />
named <- function(x) {<br />
  if (is.null(names(x))) return(NULL)<br />
  x[names(x) != ""]<br />
}<br />
unnamed <- function(x) {<br />
  if (is.null(names(x))) return(x)<br />
  x[names(x) == ""]<br />
}<br />

我们还需要用到一个html_attributes，其实它就是把这些命名的组合了起来，作者也没有列出这个函数的细节

<br />
p <- function(...) {<br />
  args <- list(...)<br />
  attribs <- html_attributes(named(args))<br />
  children <- unlist(escape(unnamed(args)))<br />
  html(paste0(<br />
    "<p", attribs, ">",<br />
    paste(children, collapse = ""),<br />
    "</p>"<br />
  ))<br />
}

其实很简单，就是把命名参数给提取出来然后组合一下放入中，然后未命名的都要esape一下，注意这里named函数返回的是list所以要unlist,然后作为内容放入之间。这样一个P就做好了，这么轻松的做完，还是因为作者的框架思路太好了

那么从这个P标签扩展到其它类似标签的工厂函数也就好写了：

<br />
tag <- function(tag) {<br />
  force(tag)<br />
  function(...) {<br />
    args <- list(...)<br />
    attribs <- html_attributes(named(args))<br />
    children <- unlist(escape(unnamed(args)))</p>
<p>    html(paste0(<br />
      "<", tag, attribs, ">",<br />
      paste(children, collapse = ""),<br />
      "</", tag, ">"<br />
    ))<br />
  }<br />
}<br />

这样的话根据输入的字符就可以获得对应的函数

p <- tag("p")<br />
b <- tag("b")<br />
i <- tag("i")<br />
p("Some text.", b("Some bold text"), i("Some italic text"),<br />
  class = "mypara")<br />
#> <HTML> <p class = 'mypara'>Some text.<b>Some bold text</b><i>Some<br />
#> italic text</i></p>

如果稍微熟悉点HTML应该知道类似这种标签，它不能有未命名的参数，所以对这类标签，我们这样处理：

<br />
void_tag <- function(tag) {<br />
  force(tag)<br />
  function(...) {<br />
    args <- list(...)<br />
    if (length(unnamed(args)) > 0) {<br />
      stop("Tag ", tag, " can not have children", call. = FALSE)<br />
    }<br />
    attribs <- html_attributes(named(args))</p>
<p>    html(paste0("<", tag, attribs, " />"))<br />
  }<br />
}<br />

这种标签更简单，好了，我们下面牺牲篇幅列出所有HTML标签，让大家感受下

tags <- c("a", "abbr", "address", "article", "aside", "audio", "b",<br />
  "bdi", "bdo", "blockquote", "body", "button", "canvas", "caption",<br />
  "cite", "code", "colgroup", "data", "datalist", "dd", "del",<br />
  "details", "dfn", "div", "dl", "dt", "em", "eventsource",<br />
  "fieldset", "figcaption", "figure", "footer", "form", "h1", "h2",<br />
  "h3", "h4", "h5", "h6", "head", "header", "hgroup", "html", "i",<br />
  "iframe", "ins", "kbd", "label", "legend", "li", "mark", "map",<br />
  "menu", "meter", "nav", "noscript", "object", "ol", "optgroup",<br />
  "option", "output", "p", "pre", "progress", "q", "ruby", "rp",<br />
  "rt", "s", "samp", "script", "section", "select", "small", "span",<br />
  "strong", "style", "sub", "summary", "sup", "table", "tbody",<br />
  "td", "textarea", "tfoot", "th", "thead", "time", "title", "tr",<br />
  "u", "ul", "var", "video")</p>
<p>void_tags <- c("area", "base", "br", "col", "command", "embed",<br />
  "hr", "img", "input", "keygen", "link", "meta", "param", "source",<br />
  "track", "wbr")

到这里，我们其实基本的已经够用了，但是不方便，每次都要提前弄好p..之类的tags很不方便，于是结合NSE的方式，我们可以这样（我自己想不到，只有看了作者写的才恍然大悟）

<br />
with_html <- function(code) {<br />
  eval(substitute(code), tag_fs)<br />
}<br />
with_html(body(<br />
  h1("A heading", id = "first"),<br />
  p("Some text &", b("some bold text.")),<br />
  img(src = "myimg.png", width = 100, height = 100)<br />
))<br />
#> <HTML> <body><h1 id = 'first'>A heading</h1><p>Some text<br />
#> &<b>some bold text.</b></p><img src = 'myimg.png' width =<br />
#> '100' height = '100' /></body><br />

这样的话,我们立刻明白了需要构造出这样的tag_fs list,它的组成部分名字应该是p,b然后内容对应为tag("p"),tag("b")这样的闭包函数,而这应该不难

<br />
tag_fs <- c(<br />
  setNames(lapply(tags, tag), tags),<br />
  setNames(lapply(void_tags, void_tag), void_tags)<br />
)<br />

这其实很简单吧，就是对一个名字列表应用tag之类的函数，然后命名为对应的名字，最后用c结合起来，还记得c这个可以组合list的函数吗？

基本上到这里，我们就可以使用with_html来写HTML了！！其实我觉得最关键的还是NSE的eval(substitute(code),tag_fs)的使用，简直太经典了，自己内部就做了p -> tag("p")的转换，这个思想我觉得以后还会遇到，并且R本身的with应该也是用的这个思想

由于这部分内容我们大概懂就好，暂时不需要太深入（或者说还没有能力深入[s:11]），所以先跳过作者给的任何练习吧，下面也一样

第二，LaTeX

下面我们要把R code转成LaTex语句，分几种情况，这里假设大家有一定的LaTex基础了

首先是已知符号：

greek <- c(<br />
  "alpha", "theta", "tau", "beta", "vartheta", "pi", "upsilon",<br />
  "gamma", "gamma", "varpi", "phi", "delta", "kappa", "rho",<br />
  "varphi", "epsilon", "lambda", "varrho", "chi", "varepsilon",<br />
  "mu", "sigma", "psi", "zeta", "nu", "varsigma", "omega", "eta",<br />
  "xi", "Gamma", "Lambda", "Sigma", "Psi", "Delta", "Xi", "Upsilon",<br />
  "Omega", "Theta", "Pi", "Phi")<br />
greek_list <- setNames(paste0("\\", greek), greek)<br />
greek_env <- list2env(as.list(greek_list), parent = emptyenv())

注意这里，前面的tag_fs那是用了lapply再用c所以结果是个list,而这里结果先是个vector所以得转成list再转成env

这样的话，我们就可以了

to_math <- function(x) {<br />
  expr <- substitute(x)<br />
  eval(expr, latex_env(expr))<br />
}

其次，对于未知的符号，我们再一次walk with ast treh

all_names <- function(x) {<br />
  if (is.atomic(x)) {<br />
    character()<br />
  } else if (is.name(x)) {<br />
    as.character(x)<br />
  } else if (is.call(x) || is.pairlist(x)) {<br />
    children <- lapply(x[-1], all_names)<br />
    unique(unlist(children))<br />
  } else {<br />
    stop("Don't know how to handle type ", typeof(x), call. = FALSE)<br />
  }<br />
}</p>
<p>all_names(quote(x + y + f(a, b, c, 10)))<br />
#> [1] "x" "y" "a" "b" "c"

这里有一个小问题就是f(a=1,b,c,10)中的a不会被提取，也许这是对的，因为它只是个参数名称而已，具体不深入了，然后就要为这个做个env

latex_env <- function(expr) {<br />
  names <- all_names(expr)<br />
  symbol_list <- setNames(as.list(names), names)<br />
  symbol_env <- list2env(symbol_list)</p>
<p>  symbol_env<br />
}

但是显然会出现问题，因为如果对于已经的名字它也只会返回对应字符串，所以我们的解决方法是让symbol_env成为greek_env的父环境，这样首先在子环境中，找不到才去父环境，为了解决这个问题，我们得自己写个函数，R本身没有

clone_env <- function(env, parent = parent.env(env)) {<br />
  list2env(as.list(env), parent = parent)<br />
}

其实很简单，就是把子环境转成list,然后设置parent,然后转会env，我们把这个放回

latex_env <- function(expr) {<br />
  # Unknown symbols<br />
  names <- all_names(expr)<br />
  symbol_list <- setNames(as.list(names), names)<br />
  symbol_env <- list2env(symbol_list)</p>
<p>  # Known symbols<br />
  clone_env(greek_env, symbol_env)<br />
}

这样就OK了

再次，对于已知函数

unary_op <- function(left, right) {<br />
  force(left)<br />
  force(right)<br />
  function(e1) {<br />
    paste0(left, e1, right)<br />
  }<br />
}</p>
<p>binary_op <- function(sep) {<br />
  force(sep)<br />
  function(e1, e2) {<br />
    paste0(e1, sep, e2)<br />
  }<br />
}<br />
# Binary operators<br />
f_env <- new.env(parent = emptyenv())<br />
f_env$"+" <- binary_op(" + ")<br />
f_env$"-" <- binary_op(" - ")<br />
f_env$"*" <- binary_op(" * ")<br />
f_env$"/" <- binary_op(" / ")<br />
f_env$"^" <- binary_op("^")<br />
f_env$"[" <- binary_op("_")</p>
<p># Grouping<br />
f_env$"{" <- unary_op("\\left{ ", " \\right}")<br />
f_env$"(" <- unary_op("\\left( ", " \\right)")<br />
f_env$paste <- paste</p>
<p># Other math functions<br />
f_env$sqrt <- unary_op("\\sqrt{", "}")<br />
f_env$sin <- unary_op("\\sin(", ")")<br />
f_env$log <- unary_op("\\log(", ")")<br />
f_env$abs <- unary_op("\\left| ", "\\right| ")<br />
f_env$frac <- function(a, b) {<br />
  paste0("\\frac{", a, "}{", b, "}")<br />
}</p>
<p># Labelling<br />
f_env$hat <- unary_op("\\hat{", "}")<br />
f_env$tilde <- unary_op("\\tilde{", "}")

为了使得sin(sin)可以工作，f_env应该是最后个被查找的env,也就是要使得它为symbol_env的父环境

latex_env <- function(expr) {<br />
  # Known functions<br />
  f_env</p>
<p>  # Default symbols<br />
  names <- all_names(expr)<br />
  symbol_list <- setNames(as.list(names), names)<br />
  symbol_env <- list2env(symbol_list, parent = f_env)</p>
<p>  # Known symbols<br />
  greek_env <- clone_env(greek_env, parent = symbol_env)<br />
}<br />

最后，我们来对于未知函数，又要来次 walk with ast tree[s:11]

all_calls <- function(x) {<br />
  if (is.atomic(x) || is.name(x)) {<br />
    character()<br />
  } else if (is.call(x)) {<br />
    fname <- as.character(x[[1]])<br />
    children <- lapply(x[-1], all_calls)<br />
    unique(c(fname, unlist(children)))<br />
  } else if (is.pairlist(x)) {<br />
    unique(unlist(lapply(x[-1], all_calls), use.names = FALSE))<br />
  } else {<br />
    stop("Don't know how to handle type ", typeof(x), call. = FALSE)<br />
  }<br />
}</p>
<p>all_calls(quote(f(g + b, c, d(a))))<br />
#> [1] "f" "+" "d"

这段代码我也没细读，大概是说提取不重复的函数名字，我们需要为这些函数特别写个函数

unknown_op <- function(op) {<br />
  force(op)<br />
  function(...) {<br />
    contents <- paste(..., collapse = ", ")<br />
    paste0("\\mathrm{", op, "}(", contents, ")")<br />
  }<br />
}

好了，最后的最后，我们要让这个未命名的函数环境成为命名的函数环境的父环境

latex_env <- function(expr) {<br />
  calls <- all_calls(expr)<br />
  call_list <- setNames(lapply(calls, unknown_op), calls)<br />
  call_env <- list2env(call_list)</p>
<p>  # Known functions<br />
  f_env <- clone_env(f_env, call_env)</p>
<p>  # Default symbols<br />
  symbols <- all_names(expr)<br />
  symbol_list <- setNames(as.list(symbols), symbols)<br />
  symbol_env <- list2env(symbol_list, parent = f_env)</p>
<p>  # Known symbols<br />
  greek_env <- clone_env(greek_env, parent = symbol_env)<br />
}<br />

OK，结束！[s:11]

Ihavenothing

越来越深了啊！

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：这里针对这个环境知识，有个讨论