05-进阶技巧.Rmd

# 进阶技巧 

程序的可读性和执行效率是两个很重要的要求.这一章主要介绍如何提高这两点.

## apply函数族

**apply并不能提高执行效率,只能使代码简洁易读.**

详细的介绍可以在[apply函数族介绍](http://blog.fens.me/r-apply/)查看,其中给了一些简单的例子.下面演示一下稍微复杂的用法.


对矩阵按列进行操作,每列的奇数项求和,偶数项求和:
```{r}
m <- matrix(1:12,4,3)
m
apply(m, 2, tapply,rep(1:2,2), sum)
```


对矩阵按列进行操作,每列前两项求和,后两项求和:

```{r}
apply(m, 2, tapply,rep(1:2,each=2), sum)
```


下面看一个3维情况的例子.




```{r}
a <- array(1:24,dim = c(2,3,4))
a
```
固定1个维度,对另外2个维度求和:

```{r}
apply(a,1,sum)
```
这样看起来不太直观,我们利用aperm函数调整数组的维度顺序:
```{r}
aperm(a,c(2,3,1))
```

固定2个维度,对1个维度求和:

```{r}
apply(a,c(2,3),sum)
```

固定2个维度,对1个维度分组求和:
```{r}
apply(a,c(1,2),tapply,rep(c(-1,-2),2), sum)
```
对第三个维度奇数项、偶数项分别求和(奇数项是-1组,偶数项是-2组).

## 并行

**单次模拟时间越长,重复次数越多,并行得到的提升越明显.**


R中有很多并行包,可以在
[不同并行包比较](https://yulongniu.bionutshell.org/blog/2014/06/25/parallel-package/)查看对比.

这里介绍的foreach是比较友好的一个包.

下面以线性模型估计系数为例,给出示例代码.

sim_single可以在
[单次估计程序]("code/parallel-demo.R")
查看.其中SLP参数用于增加单次模拟的计算量(运行时间).


### 低重复次数,低计算量

```{r,cache=TRUE, cache.path="code/cache/"}
source("code/parallel-demo.R")
library("foreach")
library("doParallel")
beta0 = c(1,-2,3)
N = c(50,100,200)
distribution= c(rnorm,rcauchy)
SIM = 500
tstart=Sys.time()
cl <- makeCluster(detectCores())
registerDoParallel(cl)
result <- foreach(n=N,.combine = rbind) %:%
  foreach(dst=distribution,.combine = c) %:%
  foreach(i=1:SIM,.combine = '+',
          .packages = c("MASS") )%dopar%{
            sim_single(n,beta0,SLP=FALSE,mydist = dst)
          }
stopImplicitCluster()
t.end=Sys.time()
t.end-tstart
result/SIM
```


把申请cluster的命令去掉,%dopar%换成%do%程序就会按串行执行
```{r ,cache=TRUE, cache.path="code/cache/"}
source("code/parallel-demo.R")
library("foreach")
library("doParallel")
beta0 =  c(1,-2,3)
N = c(50,100,200)
distribution= c(rnorm,rcauchy)
SIM = 500
tstart=Sys.time()
#cl <- makeCluster(detectCores())
#registerDoParallel(cl)
result <- foreach(n=N,.combine = rbind) %:%
  foreach(dst=distribution,.combine = c) %:%
  foreach(i=1:SIM,.combine = '+',
          .packages = c("MASS") )%do%{
            sim_single(n,beta0,SLP=FALSE,mydist = dst)
          }
#stopImplicitCluster()
t.end=Sys.time()
t.end-tstart
result/SIM
```


### 高重复次数,低计算量

增加到1000次.

并行:
```{r,cache=TRUE, cache.path="code/cache/"}
source("/Users/wang/Documents/GitHub/Simulation-in-R/code/parallel-demo.R")
library("foreach")
library("doParallel")
beta0 = c(1,-2,3)
N = c(50,100,200)
distribution= c(rnorm,rcauchy)
SIM = 1000
tstart=Sys.time()
cl <- makeCluster(detectCores())
registerDoParallel(cl)
result <- foreach(n=N,.combine = rbind) %:%
  foreach(dst=distribution,.combine = c) %:%
  foreach(i=1:SIM,.combine = '+',
          .packages = c("MASS") )%dopar%{
            sim_single(n,beta0,SLP=FALSE,mydist = dst)
          }
stopImplicitCluster()
t.end=Sys.time()
t.end-tstart
result/SIM
```


串行:
```{r ,cache=TRUE, cache.path="code/cache/"}
source("code/parallel-demo.R")
library("foreach")
library("doParallel")
beta0 =  c(1,-2,3)
N = c(50,100,200)
distribution= c(rnorm,rcauchy)
SIM = 1000
tstart=Sys.time()
#cl <- makeCluster(detectCores())
#registerDoParallel(cl)
result <- foreach(n=N,.combine = rbind) %:%
  foreach(dst=distribution,.combine = c) %:%
  foreach(i=1:SIM,.combine = '+',
          .packages = c("MASS") )%do%{
            sim_single(n,beta0,SLP=FALSE,mydist = dst)
          }
#stopImplicitCluster()
t.end=Sys.time()
t.end-tstart
result/SIM
```

### 低重复次数,高计算量

增加每次模拟的计算量,延长单次模拟的时间.

并行:
```{r ,cache=TRUE, cache.path="code/cache/"}
source("/Users/wang/Documents/GitHub/Simulation-in-R/code/parallel-demo.R")
library("foreach")
library("doParallel")
beta0 = c(1,-2,3)
N = c(50,100,200)
distribution= c(rnorm,rcauchy)
SIM = 500
tstart=Sys.time()
cl <- makeCluster(detectCores())
registerDoParallel(cl)
result <- foreach(n=N,.combine = rbind) %:%
  foreach(dst=distribution,.combine = c) %:%
  foreach(i=1:SIM,.combine = '+',
          .packages = c("MASS") )%dopar%{
            sim_single(n,beta0,SLP=TRUE,mydist = dst)
          }
stopImplicitCluster()
t.end=Sys.time()
t.end-tstart
result/SIM
```

串行:
```{r ,cache=TRUE, cache.path="code/cache/"}
source("code/parallel-demo.R")
library("foreach")
library("doParallel")
beta0 =  c(1,-2,3)
N = c(50,100,200)
distribution= c(rnorm,rcauchy)
SIM = 500
tstart=Sys.time()
#cl <- makeCluster(detectCores())
#registerDoParallel(cl)
result <- foreach(n=N,.combine = rbind) %:%
  foreach(dst=distribution,.combine = c) %:%
  foreach(i=1:SIM,.combine = '+',
          .packages = c("MASS") )%do%{
            sim_single(n,beta0,SLP=TRUE,mydist = dst)
          }
#stopImplicitCluster()
t.end=Sys.time()
t.end-tstart
result/SIM
```


### 模型选择

除了蒙特卡洛模拟之外,另外一个很适合并行的应用是模型选择.下面的例子用AIC对线性模型进行选择:

```{r,cache=TRUE, cache.path="code/cache/"}
AIClm <- function(Y,X,ind){#单次计算AIC,ind为纳入模型的变量下标
  AIC(lm(Y~X[,ind]))
}
gen_ind <- function(bt, index) bt*index#生成候选模型下标

n = 100
p = 18
m0 = p%/%3
p0 = sample(1:p,m0)
beta = rep(0,p)
beta[p0] = 1:2
X <- matrix(rnorm(n*p),n,p)
Y <- X%*%beta + rnorm(n)

##生成所有候选模型的下标
pl <- 1:(2^p - 1)
mt <- matrix(0,2^p-1,p)
for (i in pl) {
  mt[i,]=rev(as.integer(intToBits(i)[1:(p)] ))
}
index = 1:p
lst = t(apply(mt, 1, gen_ind,index=index))

#并行
library("foreach")
library("doParallel")
tstart=Sys.time()
cl <- makeCluster(detectCores())
registerDoParallel(cl)
foreach(i=1:(2^p-1),.combine = c) %dopar%{
  AIClm(Y,X,lst[i,])
}->result
stopImplicitCluster()
t.end=Sys.time()
t.end-tstart
which(lst[which.min(result),]!=0)
which(beta!=0)

#串行
tstart=Sys.time()
cl <- makeCluster(detectCores())
registerDoParallel(cl)
foreach(i=1:(2^p-1),.combine = c) %do%{
  AIClm(Y,X,lst[i,])
}->result
stopImplicitCluster()
t.end=Sys.time()
t.end-tstart
which(lst[which.min(result),]!=0)
which(beta!=0)

```




## Rcpp


Rcpp提供了R与C++的无缝接口,可以很方便的在R中调用编写的C++程序.


[Rcpp文档](https://teuder.github.io/rcpp4everyone_en/)

[如何改写R程序](https://adv-r.hadley.nz/rcpp.html)


[Rcpp已提供的分布函数](https://teuder.github.io/rcpp4everyone_en/220_dpqr_functions.html#list-of-probability-distribution-functions)


可以通过cppFunction直接在R中编写:
```{r}
Rcpp::cppFunction('int add(int x, int y, int z) {
  int sum = x + y + z;
  return sum;
}')
add
add(1, 2, 3)
```
但是这种方式在编写(没有语法高亮)、调试(定位编译报错行号)时都有不便.推荐单独编写cpp文件,通过sourceCpp加载.


下面以矩阵按列求和为例,比较col_mean(自己编写的C++函数)[^1]、colMeans(R base中提供的注重速度的函数)、mean_R(在R中用编写的函数)以及apply的运行速度.



```{r}
Rcpp::sourceCpp('code/Rcpp-demo.cpp')
m = 200
n = 100
X <- matrix(rnorm(m*n),m,n)
col_mean(X) -> l1
mean_R(X) -> l2
all.equal(l1,l2)
colMeans(X) -> l3
all.equal(l1,l3)
apply(X, 2, mean) -> l4
all.equal(l1,l4)
bench::mark(
  col_mean(X),
  colMeans(X),
  mean_R(X),
  apply(X, 2, mean),
  check = FALSE,relative = TRUE
)->results
ggplot2::autoplot(results)
```

对比结果如图所示,其中gc[^2]是一个关于内存使用的指标,越低越好.

apply和在R中用循环编写函数速度差不多,用C++编写的函数明显比其他快,甚至比base库中的函数还要快.不过,当我们增加矩阵的大小时,就会发现不一样的结果:

```{r}
Rcpp::sourceCpp('code/Rcpp-demo.cpp')
m = 2000
n = 1000
X <- matrix(rnorm(m*n),m,n)
col_mean(X) -> l1
mean_R(X) -> l2
all.equal(l1,l2)
colMeans(X) -> l3
all.equal(l1,l3)
apply(X, 2, mean) -> l4
all.equal(l1,l4)
bench::mark(
  col_mean(X),
  colMeans(X),
  mean_R(X),
  apply(X, 2, mean),
  check = FALSE,relative = TRUE
)->results
ggplot2::autoplot(results)
```

可以看到,还是base库中提供的函数速度最快.


### RcppParallel

C++并行库

[官方文档](https://rcppcore.github.io/RcppParallel/index.html)


这里提供一个RcppParallel的例子:
[核估计](https://github.com/Ri0016/kernelCpp)


下面是三个关于线性代数的库,https://gist.github.com/wolfv/ca3ac2b24e1daf70f85eac18ec7b1b8f
这个例子的测试结果表明xtensor最快


### RcppArmadillo

线性代数库
[官方文档](https://cran.r-project.org/web/packages/RcppArmadillo/RcppArmadillo.pdf)


### RcppEigen

线性代数库(更快一点,但是不友好)
[官方文档](https://cran.r-project.org/web/packages/RcppEigen/RcppEigen.pdf)

### xtensor


[GitHub地址](https://github.com/QuantStack/xtensor-r)



[^1]:程序在 [这里]("code/Rcpp-demo.cpp") 查看,其中最后被/\*\*\*R \*/夹住的部分是R程序,每次加载后会自动执行,方便调试.
[^2]:垃圾回收（英语：Garbage Collection，缩写为GC），在计算机科学中是一种自动的存储器管理机制。当一个计算机上的动态存储器不再需要时，就应该予以释放，以让出存储器，这种存储器资源管理，称为垃圾回收。