R: Как распределить многострочный график с решеткой в ​​R 3.2.1?

Question

Я новичок в программировании R и хотел бы знать, как работать параллельно plot на 12 объектов trellis сделано с lattice пакетом.

В принципе, после того, как много шагов предварительной обработки, у меня есть следующие команды:

plot(adhd_plot, split = c(1,1,4,3)) #plot adhd trellis object at 1,1 in a grid of 4 by 3 i.e 4 COLUMNS x 3 ROWS 
plot(bpd_plot, split = c(2,1,4,3), newpage = F) #plot bpd trellis object in 2nd Column in a grid of 4colx3row 
plot(bmi_plot, split = c(3,1,4,3), newpage = F) 
plot(dbp_plot, split = c(4,1,4,3), newpage = F) 
plot(height_plot, split = c(1,2,4,3), newpage = F) 
plot(hdl_plot, split = c(2,2,4,3), newpage = F) 
plot(ldl_plot, split = c(3,2,4,3), newpage = F) 
plot(ra_plot, split = c(4,2,4,3), newpage = F) 
plot(sbp_plot, split = c(1,3,4,3), newpage = F) 
plot(scz_plot, split = c(2,3,4,3), newpage = F) 
plot(tc_plot, split = c(3,3,4,3), newpage = F) 
plot(tg_plot, split = c(4,3,4,3), newpage = F) 

Проблема в том, что в то время как эти команды работают, они принимают очень долго (> 4 часа) на Mac OSX для получения такой цифры, как:

Поскольку мой Mac имеет 8 ядер, мне показалось, что я должен попытаться разделить команду графика на разные ядра, чтобы ускорить построение графика.

После поиска по другим вопросам распараллеливания, я нашел doParallel пакет и думал, что я мог бы потенциально реализовать функцию parLapply в нем, как выглядит следующим образом:

library(doParallel) 
detectCores() 
cl <- makeCluster(6) #6 out of 8 cores 
registerdoParallel(cl) 
parLapply(cl, list_of_all_trellis_objects, plot) 

Однако, я не знаю, как использовать параметр split в приведенной выше команде parLapply для размещения участков в разных местах сетки.

Мне обязательно нужны 12 участков, размещенных отдельно и не накладываемых друг на друга, так как это сделать?

Благодарим вас за рассмотрение моего запроса, и я с нетерпением жду ваших намеков и решений.

Answer 1

Как было предложено в комментариях, нет возможности писать для построения устройства параллельно.

Некоторые обходные пути, чтобы ускорить рисование отдельных участков:

1. Снизить количество точек в QQ сюжета, см:

   https://stats.stackexchange.com/questions/35220/removing-extraneous-points-near-the-centre-of-a-qq-plot

2. Загрузка данных быстрее, применяя эти советы :

   http://cbio.ensmp.fr/~thocking/reading-large-text-files-into-R.html

3. Вы можете попытаться рисовать/сохранять несколько графиков параллельно (где каждый сюжет использует методы из пунктов 1 и 2), но запись на диск может стать значительным узким местом.

Edit:

Вот грубый код рисовать быстрый Qq-сюжет:

https://github.com/vforget/fastqq

код ниже:

find_conf_intervals = function(row){ 
    i = row[1] 
    len = row[2] 
    if (i < 10000 | i %% 100 == 0){ 
    return(c(-log10(qbeta(0.95,i,len-i+1)), -log10(qbeta(0.05,i,len-i+1)))) 
    } else { # Speed up 
    return(c(NA,NA)) 
    } 
} 

confidence.intervals <- function(e){ 
    xspace = 0.078 
    print("1") 
    ci = apply(cbind(1:length(e), rep(length(e),length(e))), MARGIN=1, FUN=find_conf_intervals) 
    print("2") 
    bks = append(seq(10000,length(e),100),length(e)+1) 
    print("3") 
    for (i in 1:(length(bks)-1)){ 
    ci[1, bks[i]:(bks[i+1]-1)] = ci[1, bks[i]] 
    ci[2, bks[i]:(bks[i+1]-1)] = ci[2, bks[i]] 
    } 
    colnames(ci) = names(e) 
    ## Extrapolate to make plotting prettier (doesn't affect intepretation at data points) 
    slopes = c((ci[1,1] - ci[1,2])/(e[1] - e[2]), (ci[2,1] - ci[2,2])/(e[1] - e[2])) 
    print("4") 
    extrap_x = append(e[1]+xspace,e) ## extrapolate slightly for plotting purposes only 
    extrap_y = cbind(c(ci[1,1] + slopes[1]*xspace, ci[2,1] + slopes[2]*xspace), ci) 
    print("5") 
    polygon(c(extrap_x, rev(extrap_x)), c(extrap_y[1,], rev(extrap_y[2,])), 
      col = "grey81", border = "grey81") 
} 

quant.subsample <- function(y, m=100, e=1) { 
    ## m: size of a systematic sample 
    ## e: number of extreme values at either end to use 
    x <- sort(y) 
    n <- length(x) 
    quants <- (1 + sin(1:m/(m+1) * pi - pi/2))/2 
    sort(c(x[1:e], quantile(x, probs=quants), x[(n+1-e):n])) 
    ## Returns m + 2*e sorted values from the EDF of y 
} 

get.points <- function(pv) { 
    suppressWarnings(as.numeric(pv)) 
    names(d) = names(pv) 
    d = d[!is.na(d)] 
    d = d[d>0 & d<1] 
    d = d[order(d,decreasing=F)] 
    y = -log10(d) 
    x = -log10(ppoints(length(d))) 
    m <- 0.001 * length(x) 
    e <- floor(0.0005 * length(x)) 
    return(list(x=quant.subsample(x, m, e), y=quant.subsample(y, m, e))) 
} 

fqq <- function(x, y, ...) { 
    plot(0, 
     col=FALSE, 
     xlim=range(x), 
     ylim=range(y), 
     xlab=expression(Expected~~-log[10](italic(p))), 
     ylab=expression(Observed~~-log[10](italic(p))), 
     ...) 
    abline(0,1,col=2) 
    points(x,y, ...) 
} 

args <- commandArgs(trailingOnly = TRUE) 
pv.f = args[1] 
qq.f = args[2] 
nrows = as.numeric(args[3]) 
message(Sys.time()) 
message("READING") 
d <- read.table(pv.f, header=TRUE, sep=" ", nrows=nrows, colClasses=c("numeric")) 
message(Sys.time()) 
message("LAMBDA") 
chisq <- qchisq(1-d$P_VAL,1) 
lambda = median(chisq)/qchisq(0.5,1) 
message(Sys.time()) 
message("PLOTING") 
p <- get.points(d$P_VAL) 
png(file=qq.f) 
fqq(p$x, p$y, main=paste(pv.f, lambda, sep="\n"), cex.axis=1.5, cex.lab=1.5) 
dev.off() 
message(Sys.time())