npoint

#!/usr/bin/env Rscript

#### Make sure packages are installed and load them

if (!suppressWarnings(suppressPackageStartupMessages(require(Rniftilib)))) install.packages("http://ascopa.server4you.net/ubuntu/ubuntu/pool/universe/r/r-cran-rniftilib/r-cran-rniftilib_0.0-35.r79.orig.tar.xz", repos=NULL)
if (!suppressWarnings(suppressPackageStartupMessages(require(argparse)))) install.packages("argparse")
if (!suppressWarnings(suppressPackageStartupMessages(require(doParallel)))) install.packages("doParallel")
suppressWarnings(suppressPackageStartupMessages(library(Rniftilib)))
suppressWarnings(suppressPackageStartupMessages(library(argparse)))
suppressWarnings(suppressPackageStartupMessages(library(doParallel)))
suppressWarnings(suppressPackageStartupMessages(library(neuropointillist)))

#### Take in command line arguments

parser <- ArgumentParser(description="This program prepares your MRI data for group-level mixed effects modeling")

parser$add_argument("-m", "--mask", nargs=1, type="character", help="Mask limiting the voxels that will be analyzed", required=TRUE)
parser$add_argument("--set1", nargs=1,  help="List of files at first occasion", required=TRUE)
parser$add_argument("--set2", nargs=1, help="List of files at second occasion")
parser$add_argument("--set3", nargs=1, help="List of files at third occasion")
parser$add_argument("--set4", nargs=1, help="List of files at fourth occasion")
parser$add_argument("--set5", nargs=1, help="List of files at fifth occasion")
parser$add_argument("--setlabels1", nargs=1, help="Covariates for files at first occasion", required=TRUE)
parser$add_argument("--setlabels2", nargs=1, help="Covariates for files at second occasion")
parser$add_argument("--setlabels3", nargs=1, help="Covariates for files at third occasion")
parser$add_argument("--setlabels4", nargs=1,help="Covariates for files at fourth occasion")
parser$add_argument("--setlabels5", nargs=1, help="Covariates for files at fifth occasion")
parser$add_argument("--model", nargs=1, help="R code that defines the voxelwise-model and any initialization", required=TRUE)
parser$add_argument("--covariates", nargs=1, type="character", help="Covariates that will be merged with the design matrix")
parser$add_argument("--output", nargs=1, type="character", help="Output prefix to prepend to output files", required=TRUE)
parser$add_argument("--debugfile", nargs=1, type="character", help="Save voxeldat and designmat objects to this file to develop, test and debug the processVoxel function")
parser$add_argument("-t", "--testvoxel", type="integer", help="Specify a voxel on which the model works to determine output files", default="-1")
parser$add_argument("-p", "--processors", type="integer", help="Run using shared memory with p processors")
parser$add_argument("--sgeN",  type="integer", nargs=1, help="Run using SGE generating N jobs")

if (file.exists("readargs.R")) {
    source("readargs.R")
    if (exists("cmdargs")) { # check that cmdargs is defined
        args <- parser$parse_args(cmdargs)
    } else {
        args <- parser$parse_args()
    }
} else {
    args <- parser$parse_args()
}


###############################################################################

#### Check for mask and read it. It is mandatory and must exist.
maskfile <- args$mask
tryCatch({
    mask <- nifti.image.read(maskfile);
}, error=function(e) {
    cat("Could not read mask file: ", maskfile, "\n")
    stop(e)
})

# save mask dimensions
mask.dims <- dim(mask)

# reduce to vector and obtain list of nonzero vertices and the x,y,x
mask.vector <- as.vector(mask[,,])
mask.vertices <- which(mask.vector > 0)

# assemble the indices to create a reverse lookup table
mask.arrayindices <- data.frame(which(mask[,,] > 0, arr.in=TRUE))
mask.arrayindices$vertex <- mask.vertices

#### Save original arguments for writing out calling info
origargs <- args


#### Do argument checking
args <- npointCheckArguments(args)

#### Are we running in parallel?
if (!is.null(args$processors) || !is.null(args$sgeN)) {
    runningParallel =TRUE
} else {
    runningParallel= FALSE
}


#### A lookup function to convert FSL indices to a vertex number for testing
imagecoordtovertex <- function(x,y,z) {
    # first add 1 to convert from index at zero to index at 1
    nx <- x+1
    ny <- y+1
    nz <- z+1
    row <- mask.arrayindices[which(mask.arrayindices$dim1==nx&
                                   mask.arrayindices$dim2==ny&
                                   mask.arrayindices$dim3==nz),]
    if(is.data.frame(row) && nrow(row)==0) {
        warning("This coordinate is not in the mask; returning 1")
        return(1)
    } else {
        return(as.integer(row.names(row)))
    }
}



###############################################################################
#### Calculate the number of data sets
numberdatasets <- sum(!is.null(args$set1),
                          !is.null(args$set2),
                          !is.null(args$set3),
                          !is.null(args$set4),
                          !is.null(args$set5))

###############################################################################
#### Read in all the data sets

cat("Reading", numberdatasets, "data sets.\n")
data <- npointReadDataSets(args,numberdatasets,mask.vertices);
voxeldat <- data$voxeldat
designmat <-data$designmat
rm(data)
gc()

###############################################################################
#### Read in covariates if specified and merge with other covariates specified
#### on the command line
if (!is.null(args$covariates)) {
    designmat <- npointMergeDesignmatWithCovariates(designmat,args$covariates,dim(voxeldat)[1])
}

###############################################################################
### If debugging file is specified, save out design matrix and voxel matrix
### to this file
if(!is.null(args$debugfile)) {
    dir <- dirname(args$output)
    if (!dir.exists(dir)) {
        dir.create(dir, recursive=TRUE)
    }
    # add the prefix to the debug file name
    debugfilename <- paste(dir, args$debugfile,sep="/")
    save(designmat,voxeldat, imagecoordtovertex, mask.arrayindices, file=debugfilename)
}

###############################################################################
#### read model code
if (!is.null(args$model)) {
    modelfile <- args$model
    if (!file.exists(modelfile)) {
        stop("model file ", modelfile, " does not exist!")
    }
    result <- tryCatch({
        source(modelfile)
    }, error=function(e) {
        cat("There were errors in the model file ", modelfile, "\n")
        stop(e)

    })
}

# check to see that processVoxel is defined
if(!exists('processVoxel')) {
    stop("The model file did not define the processVoxel function")
}


###############################################################################
#### Do the parallel processing

nvertices <- length(mask.vertices)
if (runningParallel) {
    if (!is.null(args$processors)) {
        attach(designmat) # we attach to the designmat
        cl <- makeCluster(args$processors, type="FORK")
        cat("Exporting data to cluster.\n")
        clusterExport(cl,varlist=c("voxeldat"))
        cat("Starting parallel job using", args$processors, "cores.\n")
        cat("Use top to make sure that no threads are using more than 100% of the CPU.\n")
        system.time(results <-parSapply(cl,1:nvertices, processVoxel))
        stopCluster(cl)
        npointWriteOutputFiles(args$output,results,mask)
        npointWriteCallingInfo(origargs)
    } else { # we are using SGE
        # split up the data into chunks and write out scripts to process
        if (args$sgeN > nvertices) {
            stop("Number of SGE jobs requested is greater than the number of vertices")
        } else {
          cat("no. of vertices", nvertices, "\n")
          #sgeN * size is now larger than the # of vertices
  	  size <- ceiling(nvertices/args$sgeN)
          cat("size", trunc(nvertices/args$sgeN), "\n")
          njobs <- npointSplitDataSize(size,voxeldat,args$output,mask)
          cat("no. of jobs", njobs ,"\n")
          # save the design matrix
          designmatname <- paste(args$output, "designmat.rds", sep="")
          makefilename <- paste(dirname(args$output), "/Makefile", sep="")
          masterscript.local <- paste(dirname(args$output), "/runme.local", sep="")
          masterscript.sge <- paste(dirname(args$output), "/runme.sge", sep="")
          jobscript <- paste(dirname(args$output), "/sgejob.bash", sep="")
          saveRDS(designmat,designmatname)
          attach(designmat) # attach to the designmat
                        # test one voxel to obtain names for return vals
          if (args$testvoxel < 0) {
              args$testvoxel <- trunc(dim(voxeldat)[2]/2)
                                                 }
          tryCatch({ out <- processVoxel(args$testvoxel)},error=function(e) {
              message("error testing the model on a random voxel to determine output filenames")
              message("try providing the -t option to give a specific voxel for testing")
              message(e)
          stop("Exiting before generating makefile.")})
          npointWriteMakefile(basename(args$output), names(out), paste(getwd(), "/",args$model,sep=""), basename(designmatname), makefilename, masterscript.local)
          npointWriteSGEsubmitscript(basename(args$output), names(out), paste(getwd(), "/",args$model,sep=""), basename(designmatname), masterscript.sge,jobscript, njobs)
          npointWriteCallingInfo(origargs)
      }
    }
    } else {

    cat("Starting sequential job\n")
    cat("You might want to check whether your model is multithreaded\n")
    cat("because your code might run faster if you limit the number of threads\n")
    system.time(results <-sapply(1:nvertices, processVoxel))
    npointWriteOutputFiles(args$output,results,mask)
    npointWriteCallingInfo(origargs)
}