R中并行计算-多线程运行

并行计算： 简单来讲，就是同时使用多个计算资源来解决一个计算问题，是提高计算机系统计算速度和处理能力的一种有效手段。

一个问题被分解成为一系列可以并发执行的离散部分；
每个部分可以进一步被分解成为一系列离散指令；
来自每个部分的指令可以在不同的处理器上被同时执行；
需要一个总体的控制/协作机制来负责对不同部分的执行情况进行调度。
而在我们平时的模拟中，在一台电脑或者服务器上，就是将我们的计算任务分散到多个不同的小的核中同时进行处理。

在模拟时什么地方可以用到并行？
并行操作一般适用于重复的操作，比如重复随机按照相同分布生成数据，然后分别同时进行模拟。这里就可以用并行。亦或者我们要做permutation计算p-value等信息，也可以进行并行，因为这种操作是简单的重复即可完成。

但诸如迭代，递归等算法就很难用并行实现，这种都叫串行。因为后一个的对象需要前一个对象的信息，只能先算完前一个，再计算后一个内容。

在进行实际的模拟比较多种方法的优劣时，通常需要重复实验成百上千次，一般可对这里进行并行操作，写在这里的操作是最简单的。但会有个缺点：可能会出现挂服务器跑了半天还没出现结果，但是自己又并不知道运行到哪了的现象。虽然有一些方法可以进行查看（例如snowfall中的sfCat()函数，但是输出的结果是相对来说比较凌乱的，而且有时还会输出不了，具体用法后面会进行介绍），但是还是可能等很久才出一些结果，如果并行某一个地方维度或者代码有些小瑕疵，整段结果都没法进行输出。

所以建议，如果能将并行写到每个算法中间的话，就尽量写到每个具体算法之中（如需要permutation的写到permutation中；如要多次for循环计算统计量以及其它信息的，直接替代for循环），这样后面实际操作时也比较方便。（这样做的缺点是可能导致内存占用过多，从而使并行出错）

多线程运行可在parallel包中实现，或者foreach和doParallel

  package_list <- c("parallel","foreach","doParallel", "dplyr")

使用波士顿的数据集，拟合一个回归模型并计算MSE，循环10,000次。

# data
data(Boston)

# function -  从波士顿数据集的自举样本上计算模型拟合的mse
model.mse <- function(x) {
id <- sample(1:nrow(Boston), 200, replace = T)
mod <- lm(medv ~ ., data = Boston[id,])
mse <- mean((fitted.values(mod) - Boston$medv[id])^2)
return(mse)

}

# data set
#x.list <- lapply(1:2e4, function(x) rnorm(1e3, 0, 200))
x.list <- sapply(1:10000, list)

# detect the number of cores
n.cores <- detectCores()
n.cores

## [1] 12

本地计算机有12核。函数 clusterExport 将数据框架、加载包和其他函数复制到每个集群成员上。这时需要考虑一下并行计算是否真的有好处。如果数据集很大，复制12份并存储在内存中会产生巨大的消耗，而且可能不会加快计算速度。对于这些例子，我们需要将波士顿数据集导出到集群。由于数据集只有0.1Mb，这不会是一个问题。在处理结束时，重要的是要记得用stopCluster关闭集群。

parLapply

使用并行包中的parLapply是将计算并行化的最简单方法，因为 lapply是简单地与parLapply切换，并告诉它集群的设置。首先，我们将使用标准的lapply函数建立一个基线，然后与parLapply进行比较。

# 单核

system.time(a <- lapply(x.list, model.mse))

## 系统用时

## 14.58 0.00 14.66

# 12 核

system.time({
clust <- makeCluster(n.cores)
clusterExport(clust, "Boston")
a <- parLapply(clust, x.list, model.mse)})

## 系统用时
## 0.03 0.02 4.33

stopCluster(clust)

比标准的lapply快得多。另一个简单的功能是mclapply，它工作得非常好，甚至比parLapply更简单，但是Windows机器不支持这个功能，所以这里没有测试。

parSapply

# 单核
system.time(a <- sapply(1:1e4, model.mse))

## 系统用时
## 14.42 0.00 14.45

# 12 核

system.time({
clust <- makeCluster(n.cores)
clusterExport(clust, "Boston")

a <- parSapply(clust, 1:1e4, model.mse)})

## 系统用时
## 0.02 0.05 4.31

stopCluster(clust)

parApply

parApply函数工作原理与上面一样。数据需被转换为一个矩阵.

# 转换为矩阵
x.mat <- matrix(1:1e4, ncol = 1)

# 单核
system.time(a <- apply(x.mat, 1, model.mse))

##系统用时
## 14.27 0.00 14.32

# 12 核

system.time({
clust <- makeCluster(n.cores)
clusterExport(clust, "Boston")
a <- parApply(clust, x.mat, 1, model.mse)})

## 系统用时
## 0.00 0.03 4.30

stopCluster(clust)

foreach

foreach函数的工作方式与for循环相似。如果apply函数不合适，你需要使用for循环，foreach应该能完成这个工作。基本上，你可以把通常放在for循环中的东西放在%dopar%操作符之后。这里还有几件事情需要注意。

1. 我们使用doParallel包中的registerDoParallel注册集群
2. 需要用.combined来指定计算后的结果如何组合
3. 需要为cbind或rbind等多参数返回指定.multicombine = TRUE

对于更复杂的处理器，还有一些其他有用的参数，但这些是关键的东西。

# for

system.time({
model.mse.output <- rep(NA, 1e4)

for(k in 1:1e4){
model.mse.output[k] <- model.mse()
}})

## 系统用时
## 14.23 0.00 14.23

# foreach

system.time({
registerDoParallel(n.cores)
foreach(k = 1:1e4, .combine = c) %dopar% model.mse()
})

## 系统用时
## 3.50 0.59 6.53

stopImplicitCluster()

foreach比parXapply慢

实例

以我自己的数据为例，介绍parSapply函数的使用

1.加载包

package_list <- c("DiscriMiner","openxlsx","ggplot2",
                  "parallel","foreach","doParallel",
                  "dplyr")

for(p in package_list){
  if(!suppressWarnings(suppressMessages(require(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))){
    # install.packages(p, repos=site)
    install.packages(p)
    suppressWarnings(suppressMessages(library(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))
  }
}

2.定义函数

我的目的是进行plsda分析，但由于暴露变量和结局变量众多，且需通过定义ncomp寻找最优方程，因此需要进行批量计算。

oplsda_function<-function(j){
  source_dir<-"/home/xuyang/chen"
  #source_dir<-"C:/Users/18896/Desktop/chen"
  setwd(source_dir)
  dat<-openxlsx::read.xlsx("expo meta.xlsx",sheet = 2)
  exposome<-colnames(dat[,2:111])
  outcome<-colnames(dat[,112:250])
  dat2<-na.omit(dat)
  dir.create("oplsda/result_adjusted1")
  setwd("/home/xuyang/chen/oplsda/result_adjusted1")
  #setwd("C:/Users/18896/Desktop/chen/oplsda/newresult_adjusted1")
  result2<-data.frame()
  for (j  in 2:10) {
    result<-data.frame()
    VIP<-data.frame(exposome,rep(1,times=length(exposome)))
    rownames(VIP)<-exposome
    names(VIP)<-c("name","nomeaning")
    for (i in colnames(dat2[,112:250])){
      pls_i<-DiscriMiner::plsDA(dat2[,2:111], dat2[,i], autosel=FALSE, comps=j)
      r2_x_i<-pls_i$R2[j,2]
      r2_y_i<-pls_i$R2[j,4]
      q2_x_i<-pls_i$Q2[j,3]
      error_x_i<-pls_i$error_rate
      result_i<-c(i,r2_x_i,r2_y_i,error_x_i,j)
      result_i<-data.frame(result_i)
      result_i<-as.data.frame(t(result_i))
      result<-rbind(result,result_i)
      VIP_i<-pls_i$VIP[,3]
      # VIP_i<-sort(VIP_i,decreasing = T)
      name<-names(VIP_i)
      VIP_i<-data.frame(name,VIP_i)
      names(VIP_i)<-c("name",i)
      VIP<-dplyr::left_join(VIP,VIP_i,by='name')
    }
    
    names(result)<-c("outcome","R2","Q2","error","ncomp")
    result2<-rbind(result2,result)
    write.csv(result,file = paste("oplsda_ncomp=",j,".csv",sep = ""),row.names = F)
    write.csv(VIP,file = paste("VIP value_ncomp=",j,".csv",sep = ""))
  }
}

3.确定服务器核数 

n.cores <- detectCores()
n.cores

### 64

4.集成运算

clust <- makeCluster(60)
clusterExport(clust, "dat2")
#sapply(2:10, oplsda_function)
parSapply(clust, 2:10, oplsda_function)
stopCluster(clust)

需要注意的是，parSapply命令中，function涉及的函数使用需要在函数名前加上包的名称，如openxlsx::read.xlsx，避免找不到函数的情况。

ref:

https://blog.csdn.net/weixin_41929524/article/details/81707053

Simple Parallel Processing in R | R-bloggers