netSI <- function(x,indicator="all", d='HIM', adj.method='cor',
resamp.method="montecarlo", k=3, h=20, n.cores=NULL,save=FALSE,
symm=TRUE, verbose=TRUE, ...){
## Get the function call parameters
## NB the parameters should be evaluated with eval()
Call <- match.call()
##add a check so that an unexisting parameter cannot be passed
id.Call <- match( names(Call),c("x", "indicator", "d", "adj.method" ,
"resamp.method","k","h","n.cores","save","verbose",
"FDR","P","measure","alpha","C","DP","var.thr",
"components","n.nodes","ga","sseed","rho", "use", "symm"), nomatch=0)
if(sum(id.Call[-1]==0)==1){
warning("The parameter '",names(Call)[which(id.Call==0)[2]],"' will be ignored",call.=FALSE)
}
if(sum(id.Call[-1]==0)>1){
msg <- "The following parameters will be ignored:\n"
for(i in which(id.Call==0)[-1]){
msg <- paste(msg,"'",names(Call)[i],"'\n",sep="")
}
warning(msg,call.=FALSE)
}
id.Call <- match(c("x", "indicator", "d", "adj.method","resamp.method","k","h","n.cores"),
names(Call), nomatch=0)
if(id.Call[1]==0){
stop("A dataset must be provided",call.=FALSE)
}else{
## m <- eval(temp, parent.frame())
## x <- eval(Call$x)
if(!(is.matrix(x) | is.data.frame(x))){
stop("x must be a matrix or a data frame", call.=FALSE)
}
}
## Choose the indicators
if(id.Call[2]!=0){
indicator <- eval(Call$indicator)
if(!is.character(indicator))
stop("indicator must be one of 'S','SI','Sw','Sd','all'.", call.=FALSE)
}
INDICATORS <- c('S','SI','Sw','Sd',"all")
indicator <- pmatch(indicator,INDICATORS)
if(is.na(indicator))
stop("invalid indicator", call. =FALSE)
if(indicator == -1)
stop("ambiguous indicator", call. =FALSE)
DISTANCE <- c("HIM","IM","H",
"ipsen","Ipsen","IpsenMikhailov","Ipsen-Mikhailov",
"hamming","Hamming")
d <- pmatch(d, DISTANCE)
if(d==2L | (d>=4 & d<=7L))
d <- 2L
if(d==3L |(d>=8L & d<=9L))
d <- 3L
if(is.na(d))
stop("invalid distance", call. =FALSE)
if(d == -1)
stop("ambiguous distance", call. =FALSE)
d <- c("HIM","IM","H")[d]
##check on save and verbose
if(!is.logical(save))
stop("save must be TRUE or FALSE",call.=FALSE)
if(!is.logical(verbose))
stop("verbose must be TRUE or FALSE", call.=FALSE)
## It can be passed through ...
if (is.null(Call$sseed)){
sseed <- 0
} else {
sseed <- eval(Call$sseed)}
set.seed(sseed)
## Pass parameter gamma to netdist functions
## if(!is.null(Call$ga)){
## ga <- eval(Call$ga)
## } else {
## ga <- Call$ga}
## Pass parameter components to netdist functions
if(!is.null(Call$components)){
components <- eval(Call$components)
if(d=="HIM" & components==TRUE){
warning(paste("components parameter will be ignored. \n",
"The stability indicators will be computed just for",
d, "distance.\n",
"For computing them for Hamming or Ipsen-Mikhailov",
"distance use dist=H or dist=IM"), call.=FALSE)
components <- FALSE
}
}
## n.cores parameter
## n.cores <- eval(Call$n.cores)
## Check availability of cores, otherwise set a default
if(is.null(n.cores)){
if(detectCores()>=2){
n.cores <- detectCores()-1
cl <- makeCluster(n.cores)
warning("The computation has been automatically parallelized", call.=FALSE)
} else {
cl <- NULL
}
} else {
if (n.cores==1){
cl <- NULL
} else {
if (n.cores<detectCores()){
cl <- makeCluster(n.cores)
} else {
if(detectCores()>=2){
n.cores <- detectCores()-1
cl <- makeCluster(n.cores)
warning("The computation has been automatically parallelized", call.=FALSE)
}
}
}
}
## Get the dimension of the input matrix
ddim <- nrow(x)
## Get the resampling indexes
if(verbose==TRUE) cat("computing resampling...\n")
idxs <- resamplingIDX(ddim,resamp.method=resamp.method, k=k, h=h)
## length of the list for optimization purposes
ADJcv <- vector("list",length=length(idxs))
if(verbose==TRUE) cat("computing adjacency matrices...\n")
## Compute the adjacency matrices for each resampling index
if(!is.null(cl)){
## Parallel computation
ADJcv <- parLapply(cl=cl,X=idxs,fun=function(x,DAT,infer.method,...){
ss <- DAT[x,]
tmp <- mat2adj(x=ss, infer.method=infer.method, ...)
return(tmp)
},DAT=x,infer.method=adj.method, ...)
} else {
## One core computation
ADJcv <- lapply(X=idxs,FUN=function(x,DAT,infer.method, ...){
ss <- DAT[x,]
tmp <- mat2adj(x=ss, infer.method=infer.method, ...)
return(tmp)
},DAT=x,infer.method=adj.method, ...)
}
##computing the adjacency matrix on the whole dataset
ADJall <- mat2adj(x=x,infer.method=adj.method, ...)
## Compute the stability indicators
netsi <- list()
if(indicator==1L | indicator==5L){
if(verbose==TRUE) cat("computing stability indicator S...\n")
netsi[["S"]] <- netsiS(ADJall, ADJcv, d=d, cl=cl, ...)
}
if(indicator==3L | indicator==5L){
if(verbose==TRUE) cat("computing stability indicator Sw...\n")
netsi[["Sw"]] <- netsiSw(ADJcv, cl=cl, symm=symm)
}
if(indicator==4L | indicator==5L){
if(verbose==TRUE) cat("computing stability indicator Sd...\n")
netsi[["Sd"]] <- netsiSd(ADJcv, cl=cl)
}
if (!is.null(cl))
stopCluster(cl)
if(indicator==2L | indicator==5L){
if(verbose==TRUE) cat("computing stability indicator SI...\n")
netsi[["SI"]] <- netsiSI(ADJcv, d=d, n.cores=n.cores, ...)
}
results <- list("S"=mean(netsi[["S"]]),
"SI"=mean(netsi[["SI"]]),
"Sw"=apply(netsi[["Sw"]], 1,compute.indicator),
"Sd"=apply(netsi[["Sd"]], 2,compute.indicator)
)
if(save==TRUE){
results$call <- Call
results$ADJ <- ADJall
results$ADJlist <- ADJcv
results$S_boot <- netsi[["S"]]
results$SI_boot <- netsi[["SI"]]
results$Sw_boot <- netsi[["Sw"]]
results$Sd_boot <- netsi[["Sd"]]
}
return(results)
}
## Stability indicator S
netsiS <- function(g, H, d, cl, ...){
DIST <- c("HIM","IM","H")
type <- pmatch(d,DIST)
type <- DIST[type]
if(!is.null(cl)){
s <- parLapply(cl=cl,X=H,fun=function(x,g,type, ...){
res <- netdist(g,x,d=type, n.cores=1, ...)[[type]]
return(res)
}, g=g, type=type, ...)
}else{
s <- lapply(X=H,FUN=function(x,g,type, ...){
res <- netdist(g,x,d=type, n.cores=1, ...)[[type]]
return(res)
},g=g,type=type, ...)
}
return(unlist(s))
}
## Stability indicator SI
netsiSI <- function(H, d, ...){
DIST <- c("HIM","IM","H")
type <- pmatch(d,DIST)
type <- DIST[type]
s <- netdist(H,d=type, ...)[[1]]
return(s[upper.tri(s)])
}
## Degree stability
netsiSd <- function(H,cl){
if (length(H)){
n <- ncol(H[[1]])
} else {
stop("No adjacency matrix computed",call.=FALSE)
}
if (!is.null(cl)){
## Parallel computation
ddin <- parLapply(cl=cl, X=H, rowSums)
ddout <- parLapply(cl=cl, X=H, colSums)
} else {
## One core computation
ddin <- lapply(H, rowSums)
ddout <- lapply(H, colSums)
}
ck <- all(sapply(1:length(ddin),
function(x, ddin, ddout){
all((ddin[[x]] - ddout[[x]]) <= sqrt(.Machine$double.eps))
}, ddin=ddin, ddout=ddout))
## Check if matrix is directed
if (ck){
dd <- matrix(unlist(ddin),ncol=n, byrow=TRUE)
colnames(dd) <- colnames(H[[1]])
} else {
ddinm <- matrix(unlist(ddin),ncol=n, byrow=TRUE)
colnames(ddinm) <- paste(colnames(H[[1]]), "-in", sep="")
ddoutm <- matrix(unlist(ddout),ncol=n, byrow=TRUE)
colnames(ddoutm) <- paste(colnames(H[[1]]), "-out", sep="")
dd <- cbind(ddinm, ddoutm)
}
return(dd)
}
## Edge stability
netsiSw <- function(H,cl, symm=TRUE, ...){
if (length(H))
n <- nrow(H[[1]]) else stop("List of adjacency matrices do not exist")
if (symm){
com <- combn(1:n, 2)
} else {
tmpcom <- expand.grid(1:n, 1:n)
com <- t(tmpcom[tmpcom[,1] != tmpcom[,2],])
dimnames(com) <- list(NULL, NULL)
}
if (!is.null(cl)){
## Parallel computation
tmp <- parLapply(cl, H,
function(x,com){
sapply(1:ncol(com),
function(y, x, allcom){
x[allcom[1,y],allcom[2,y]]},
x=x, allcom=com)
},
com=com)
} else {
## One core computation
tmp <- lapply(H,
function(x,com){
sapply(1:ncol(com),
function(y, x, allcom){
x[allcom[1,y],allcom[2,y]]},
x=x, allcom=com)
},
com=com)
}
## Set up the results in a matrix
m <- matrix(unlist(tmp), ncol=length(H))
rownames(m) <- paste(com[1,],com[2,],sep="-")
colnames(m) <- names(H)
return(m)
}
## Function for the computation of resampling indexes
resamplingIDX <- function(N,resamp.method="montecarlo", k=3, h=20){
## Check of resampling methods
METHODS <- c('montecarlo','LOO','kCV')
resamp.method <- pmatch(resamp.method, METHODS)
if(is.na(resamp.method))
stop("invalid resampling method")
if(resamp.method == -1)
stop("ambiguous resampling method")
## Montecarlo
if (resamp.method==1L)
take <- lapply(1:h,function(x){tmp <- sample(seq(N),floor(N*(1-(1/k))))
return(tmp)})
## Leave One Out
if (resamp.method==2L){
if (k!=1)
warning("h is different than 1 but method is set to LOO (Leave One Out cross-validation schema).\nh will be ignored.")
h <- N
take <- lapply(1:N,function(x,allid){return(allid[which(allid!=x)])},allid=1:N)
}
## K-fold cross-validation
if (resamp.method==3L){
if (k>=N){
stop("Number of fold bigger than samples in the dataset!!!")
}else{
take <- list()
for (H in seq(h)){
tmpcv <- sample(seq(N),N)
s <- split(tmpcv,cut(seq(N),k))
names(s) <- NULL
take <- c(take,lapply(s,function(x,idx){tmp <- idx[!is.element(idx,x)]
return(tmp)},idx=tmpcv))
}
}
}
## return a list with indexes
return(take)
}
compute.indicator <- function(x){
## Compute the indicator value as Range/Mean over resamplings
tmp <- range(x)
mm <- mean(x, na.rm=TRUE)
if (all.equal(mm,0)==TRUE){
rr <- NA
} else {
rr <- (tmp[2] - tmp[1])/mm
}
return(rr)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.