1 | | * Testing the sample R program from www |
2 | | * Original R Program (without Rmpi) |
3 | | * Parallel R Program (with Rmpi) |
| 1 | * Testing the sample R program from [here http://ace.acadiau.ca/math/ACMMaC/Rmpi/sample.html] |
| 2 | * Original R Program of 10 loops (without Rmpi) |
| 3 | {{{ |
| 4 | start_time <- Sys.time() |
| 5 | # first make some data |
| 6 | n <- 1000 # number of obs |
| 7 | p <- 30 # number of variables |
| 8 | |
| 9 | x <- matrix(rnorm(n*p),n,p) |
| 10 | beta <- c(rnorm(p/2,0,5),rnorm(p/2,0,.25)) |
| 11 | y <- x %*% beta + rnorm(n,0,20) |
| 12 | thedata <- data.frame(y=y,x=x) |
| 13 | |
| 14 | # summary(lm(y~x)) |
| 15 | |
| 16 | fold <- rep(1:10,length=n) |
| 17 | fold <- sample(fold) |
| 18 | |
| 19 | rssresult <- matrix(0,p,10) |
| 20 | |
| 21 | for (j in 1:10) { |
| 22 | for (i in 1:p) { |
| 23 | templm <- lm(y~.,data=thedata[fold!=j,1:(i+1)]) |
| 24 | yhat <- predict(templm,newdata=thedata[fold==j,1:(i+1)]) |
| 25 | rssresult[i,j] <- sum((yhat-y[fold==j])^2) |
| 26 | } |
| 27 | } |
| 28 | |
| 29 | end_time <- Sys.time() |
| 30 | runTime<-difftime(end_time,start_time,units="secs") |
| 31 | write(runTime,file="single_10_runTime",append=TRUE) |
| 32 | write(rssresult,file="single_10_result") |
| 33 | q(save="no") |
| 34 | |
| 35 | }}} |
| 36 | * Parallel R Program of 10 loops (with Rmpi) |
| 37 | {{{ |
| 38 | # Create Start Timer |
| 39 | start_time<-Sys.time() |
| 40 | |
| 41 | # Initialize MPI |
| 42 | library("Rmpi") |
| 43 | |
| 44 | # Notice we just say "give us all the slaves you've got." |
| 45 | mpi.spawn.Rslaves(nslaves=2) |
| 46 | |
| 47 | if (mpi.comm.size() < 2) { |
| 48 | print("More slave processes are required.") |
| 49 | mpi.quit() |
| 50 | } |
| 51 | |
| 52 | .Last <- function(){ |
| 53 | if (is.loaded("mpi_initialize")){ |
| 54 | if (mpi.comm.size(1) > 0){ |
| 55 | print("Please use mpi.close.Rslaves() to close slaves.") |
| 56 | mpi.close.Rslaves() |
| 57 | } |
| 58 | print("Please use mpi.quit() to quit R") |
| 59 | .Call("mpi_finalize") |
| 60 | } |
| 61 | } |
| 62 | |
| 63 | # Function the slaves will call to perform a validation on the |
| 64 | # fold equal to their slave number. |
| 65 | # Assumes: thedata,fold,foldNumber,p |
| 66 | foldslave <- function() { |
| 67 | # Note the use of the tag for sent messages: |
| 68 | # 1=ready_for_task, 2=done_task, 3=exiting |
| 69 | # Note the use of the tag for received messages: |
| 70 | # 1=task, 2=done_tasks |
| 71 | junk <- 0 |
| 72 | |
| 73 | done <- 0 |
| 74 | while (done != 1) { |
| 75 | # Signal being ready to receive a new task |
| 76 | mpi.send.Robj(junk,0,1) |
| 77 | |
| 78 | # Receive a task |
| 79 | task <- mpi.recv.Robj(mpi.any.source(),mpi.any.tag()) |
| 80 | task_info <- mpi.get.sourcetag() |
| 81 | tag <- task_info[2] |
| 82 | |
| 83 | if (tag == 1) { |
| 84 | foldNumber <- task$foldNumber |
| 85 | |
| 86 | rss <- double(p) |
| 87 | for (i in 1:p) { |
| 88 | # produce a linear model on the first i variables on |
| 89 | # training data |
| 90 | templm <- lm(y~.,data=thedata[fold!=foldNumber,1:(i+1)]) |
| 91 | |
| 92 | # produce predicted data from test data |
| 93 | yhat <- predict(templm,newdata=thedata[fold==foldNumber,1:(i+1)]) |
| 94 | |
| 95 | # get rss of yhat-y |
| 96 | localrssresult <- sum((yhat-thedata[fold==foldNumber,1])^2) |
| 97 | rss[i] <- localrssresult |
| 98 | } |
| 99 | |
| 100 | # Send a results message back to the master |
| 101 | results <- list(result=rss,foldNumber=foldNumber) |
| 102 | mpi.send.Robj(results,0,2) |
| 103 | } |
| 104 | else if (tag == 2) { |
| 105 | done <- 1 |
| 106 | } |
| 107 | # We'll just ignore any unknown messages |
| 108 | } |
| 109 | |
| 110 | mpi.send.Robj(junk,0,3) |
| 111 | } |
| 112 | |
| 113 | # We're in the parent. |
| 114 | # first make some data |
| 115 | n <- 1000 # number of obs |
| 116 | p <- 30 # number of variables |
| 117 | |
| 118 | # Create data as a set of n samples of p independent variables, |
| 119 | # make a "random" beta with higher weights in the front. |
| 120 | # Generate y's as y = beta*x + random |
| 121 | x <- matrix(rnorm(n*p),n,p) |
| 122 | beta <- c(rnorm(p/2,0,5),rnorm(p/2,0,.25)) |
| 123 | y <- x %*% beta + rnorm(n,0,20) |
| 124 | thedata <- data.frame(y=y,x=x) |
| 125 | |
| 126 | fold <- rep(1:10,length=n) |
| 127 | fold <- sample(fold) |
| 128 | |
| 129 | #summary(lm(y~x)) |
| 130 | |
| 131 | # Now, send the data to the slaves |
| 132 | mpi.bcast.Robj2slave(thedata) |
| 133 | mpi.bcast.Robj2slave(fold) |
| 134 | mpi.bcast.Robj2slave(p) |
| 135 | |
| 136 | # Send the function to the slaves |
| 137 | mpi.bcast.Robj2slave(foldslave) |
| 138 | |
| 139 | # Call the function in all the slaves to get them ready to |
| 140 | # undertake tasks |
| 141 | mpi.bcast.cmd(foldslave()) |
| 142 | |
| 143 | |
| 144 | # Create task list |
| 145 | tasks <- vector('list') |
| 146 | for (i in 1:10) { |
| 147 | tasks[[i]] <- list(foldNumber=i) |
| 148 | } |
| 149 | |
| 150 | # Create data structure to store the results |
| 151 | rssresult = matrix(0,p,10) |
| 152 | |
| 153 | junk <- 0 |
| 154 | closed_slaves <- 0 |
| 155 | n_slaves <- mpi.comm.size()-1 |
| 156 | |
| 157 | while (closed_slaves < n_slaves) { |
| 158 | # Receive a message from a slave |
| 159 | message <- mpi.recv.Robj(mpi.any.source(),mpi.any.tag()) |
| 160 | message_info <- mpi.get.sourcetag() |
| 161 | slave_id <- message_info[1] |
| 162 | tag <- message_info[2] |
| 163 | |
| 164 | if (tag == 1) { |
| 165 | # slave is ready for a task. Give it the next task, or tell it tasks |
| 166 | # are done if there are none. |
| 167 | if (length(tasks) > 0) { |
| 168 | # Send a task, and then remove it from the task list |
| 169 | mpi.send.Robj(tasks[[1]], slave_id, 1); |
| 170 | tasks[[1]] <- NULL |
| 171 | } |
| 172 | else { |
| 173 | mpi.send.Robj(junk, slave_id, 2) |
| 174 | } |
| 175 | } |
| 176 | else if (tag == 2) { |
| 177 | # The message contains results. Do something with the results. |
| 178 | # Store them in the data structure |
| 179 | foldNumber <- message$foldNumber |
| 180 | rssresult[,foldNumber] <- message$result |
| 181 | } |
| 182 | else if (tag == 3) { |
| 183 | # A slave has closed down. |
| 184 | closed_slaves <- closed_slaves + 1 |
| 185 | } |
| 186 | } |
| 187 | |
| 188 | |
| 189 | # plot the results |
| 190 | # plot(apply(rssresult,1,mean)) |
| 191 | end_time<-Sys.time() |
| 192 | runTime<-difftime(end_time,start_time,units="secs") |
| 193 | write(runTime, file="parallel_10_runTime",append=TRUE) |
| 194 | write(rssresult, file="parallel_10_result") |
| 195 | mpi.close.Rslaves() |
| 196 | mpi.quit(save="no") |
| 197 | |
| 198 | }}} |
| 199 | |