################################################### ### chunk number 1: set width ################################################### options(width=60) options(continue=" ") .getmaxmem <- function() { ml <- system("free", intern=TRUE)[2] m <- sub(" .*$", "", sub("^[^1-9]*", "", ml)) format(as.numeric(m) / 1024 / 1024, digits=4) } .getnoproc <- function() { pc <- system("cat /proc/cpuinfo | grep processor | tail -1", intern=TRUE) as.numeric(sub("^[^1-9]*", "", pc)) + 1 } .getproctype <- function() { mn <- system("cat /proc/cpuinfo | grep 'model name' | head -1", intern=TRUE) sub("^.*:[ ]*", "", gsub("[\t ]+", " ", mn)) } ################################################### ### chunk number 2: load saved data ################################################### load("ISA_internals.Rdata") ################################################### ### chunk number 3: download data matrix ################################################### GEO <- "GSE18858" GEOfile <- paste(sep="", GEO, "_series_matrix.txt.gz") GEOurl <- paste(sep="", "ftp://ftp.ncbi.nih.gov/pub/geo/DATA/SeriesMatrix/", GEO, "/", GEOfile) if (!file.exists(GEOfile)) { download.file(GEOurl, GEOfile) } ################################################### ### chunk number 4: read in the data ################################################### data <- read.table(gzfile(GEOfile), comment.char="!", header=TRUE, row.names=1) data <- as.matrix(data) ################################################### ### chunk number 5: dimensions ################################################### dim(data) ################################################### ### chunk number 6: load the isa2 package ################################################### library(isa2) ################################################### ### chunk number 7: create a function that measures ISA running time ################################################### mesISA <- function(E, thr.row, thr.col, no.seeds) { t1 <- system.time({ NE <- isa.normalize(E) }) t2 <- system.time({ seeds <- generate.seeds(length=nrow(E), count=no.seeds) }) t3 <- system.time({ mods <- isa.iterate(NE, row.seeds=seeds, thr.row=thr.row, thr.col=thr.col) }) t4 <- system.time({ mods2 <- isa.unique(NE, mods) }) cbind(normalization=t1, seeds.generation=t2, isa.iteration=t3, module.merge=t4, full=t1+t2+t3+t4) } ################################################### ### chunk number 8: test the masurement function ################################################### mydata <- data[ sample(nrow(data), 1000), sample(ncol(data), 30)] mesISA(mydata, 3, 3, 100) ################################################### ### chunk number 9: running time in the function of rows set up things ################################################### row.thresholds <- seq(1,3,by=0.5) col.thresholds <- 2 ################################################### ### chunk number 10: running time in the function of rows set up things 2 ################################################### no.rows <- seq(5000, min(40000, nrow(data)), by=5000) do.row.thr <- function(thr, rep=5) { res <- lapply(no.rows, function(x) { lapply(1:rep, function(xxx) { mydata <- data[sample(nrow(data), x),] mesISA(mydata, thr, col.thresholds, 100) }) }) res } ################################################### ### chunk number 11: calculate by rows eval=FALSE ################################################### ## by.rows <- lapply(row.thresholds, do.row.thr) ################################################### ### chunk number 12: function to plot with error bars ################################################### myplot <- function(x, y, sd, xlim=range(x), ylim=c(min(y-sd),max(y+sd)), xlab="", ylab="", ...) { plot(NA, type="n", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab) xmin <- par("usr")[1] xr <- par("usr")[2]-xmin bw <- xr/200 segments(x,y-sd,x,y+sd) segments(x-bw,y-sd,x+bw,y-sd) segments(x-bw,y+sd,x+bw,y+sd) points(x, y, ...) } ################################################### ### chunk number 13: function to extract mean and sd ################################################### get.mean <- function(xx) { sapply(xx, function(x) mean(sapply(x, function(y) sum(y[1:2,5])))) } get.sd <- function(xx) { sapply(xx, function(x) sd(sapply(x, function(y) sum(y[1:2,5])))) } ################################################### ### chunk number 14: create the plot eval=FALSE ################################################### ## layout(rbind(1:2,3:4,5:6)) ## for (i in 1:length(row.thresholds)) { ## par(mar=c(5,4,1,1)+0.1) ## y <- get.mean(by.rows[[i]]) ## s <- get.sd(by.rows[[i]]) ## myplot(no.rows, y, s, ## type="b", pch=20, xlab="# of rows", ylab="running time [s]") ## rt <- row.thresholds[i] ## text(min(no.rows), max(y+s), ## substitute(Theta[r]==rt, list(rt=rt)), ## adj=c(0,1), cex=1.3) ## } ################################################### ### chunk number 15: by-rows-plot ################################################### layout(rbind(1:2,3:4,5:6)) for (i in 1:length(row.thresholds)) { par(mar=c(5,4,1,1)+0.1) y <- get.mean(by.rows[[i]]) s <- get.sd(by.rows[[i]]) myplot(no.rows, y, s, type="b", pch=20, xlab="# of rows", ylab="running time [s]") rt <- row.thresholds[i] text(min(no.rows), max(y+s), substitute(Theta[r]==rt, list(rt=rt)), adj=c(0,1), cex=1.3) } ################################################### ### chunk number 16: change number of columns set up things ################################################### row.thresholds2 <- 2 col.thresholds2 <- seq(1,3,by=0.5) ################################################### ### chunk number 17: change number of columns set up things ################################################### no.cols <- seq(30, ncol(data), by=30) do.col.thr <- function(thr, rep=5) { res <- lapply(no.cols, function(x) { lapply(1:rep, function(xxx) { mydata <- data[, sample(ncol(data), x)] mesISA(mydata, row.thresholds2, thr, 100) }) }) res } ################################################### ### chunk number 18: calculate by cols eval=FALSE ################################################### ## by.cols <- lapply(col.thresholds2, do.col.thr) ################################################### ### chunk number 19: make by cols plot eval=FALSE ################################################### ## layout(rbind(1:2,3:4,5:6)) ## for (i in 1:length(col.thresholds2)) { ## par(mar=c(5,4,1,1)+0.1) ## y <- get.mean(by.cols[[i]]) ## s <- get.sd(by.cols[[i]]) ## myplot(no.cols, y, s, ## type="b", pch=20, xlab="# of cols", ylab="running time [s]") ## rt <- col.thresholds2[i] ## text(min(no.cols), max(y+s), ## substitute(Theta[c]==rt, list(rt=rt)), ## adj=c(0,1), cex=1.3) ## } ################################################### ### chunk number 20: by-cols-plot ################################################### layout(rbind(1:2,3:4,5:6)) for (i in 1:length(col.thresholds2)) { par(mar=c(5,4,1,1)+0.1) y <- get.mean(by.cols[[i]]) s <- get.sd(by.cols[[i]]) myplot(no.cols, y, s, type="b", pch=20, xlab="# of cols", ylab="running time [s]") rt <- col.thresholds2[i] text(min(no.cols), max(y+s), substitute(Theta[c]==rt, list(rt=rt)), adj=c(0,1), cex=1.3) } ################################################### ### chunk number 21: no seeds parameters ################################################### thr.comb <- list( c(2,2), c(1,3), c(3,1) ) no.seeds <- seq(50, 400, by=50) ################################################### ### chunk number 22: no seeds function ################################################### do.no.seeds <- function(thr, rep=5) { nr <- min(20000, nrow(data)) nc <- min(100, ncol(data)) res <- lapply(no.seeds, function(x) { lapply(1:rep, function(xxx) { mydata <- data[sample(nrow(data), nr), sample(ncol(data), nc)] mesISA(mydata, thr[1], thr[2], x) }) }) res } ################################################### ### chunk number 23: calculate no needs eval=FALSE ################################################### ## by.no.seeds <- lapply(thr.comb, do.no.seeds) ################################################### ### chunk number 24: plot by no seeds eval=FALSE ################################################### ## layout(rbind(1:2,3:4)) ## for (i in 1:length(thr.comb)) { ## par(mar=c(5,4,1,1)+0.1) ## y <- get.mean(by.no.seeds[[i]]) ## s <- get.sd(by.no.seeds[[i]]) ## myplot(no.seeds, y, s, ## type="b", pch=20, xlab="# of seeds", ylab="running time [s]") ## th <- thr.comb[[i]] ## text(min(no.seeds), max(y+s), ## substitute(paste(Theta[r]==r1, ", ", Theta[c]==r2), ## list(r1=th[1], r2=th[2])), ## adj=c(0,1), cex=1.3) ## } ################################################### ### chunk number 25: by-seeds-plot ################################################### layout(rbind(1:2,3:4)) for (i in 1:length(thr.comb)) { par(mar=c(5,4,1,1)+0.1) y <- get.mean(by.no.seeds[[i]]) s <- get.sd(by.no.seeds[[i]]) myplot(no.seeds, y, s, type="b", pch=20, xlab="# of seeds", ylab="running time [s]") th <- thr.comb[[i]] text(min(no.seeds), max(y+s), substitute(paste(Theta[r]==r1, ", ", Theta[c]==r2), list(r1=th[1], r2=th[2])), adj=c(0,1), cex=1.3) } ################################################### ### chunk number 26: sessioninfo ################################################### toLatex(sessionInfo())