npointrun

#!/usr/bin/env Rscript

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

# Rniftilib is not in this location
#if (!suppressPackageStartupMessages(require(Rniftilib))) install.packages("Rniftilib", repos="http://R-Forge.R-project.org")
if (!suppressPackageStartupMessages(require(argparse))) install.packages("argparse")
if (!suppressPackageStartupMessages(require(doParallel))) install.packages("doParallel")
suppressPackageStartupMessages(library(RNifti))
suppressPackageStartupMessages(library(argparse))
suppressPackageStartupMessages(library(doParallel))
suppressPackageStartupMessages(library(neuropointillist))

# Source helper functions

parser <- ArgumentParser()
parser$add_argument("-m", "--mask",  nargs=1, help="Maskfile to process", required=TRUE)
parser$add_argument("-d", "--designmat",  nargs=1, help="Design matrix in RDS format", required=TRUE)
parser$add_argument("-p", "--permutationfile",  nargs=1, help="File name of the permutation output that includes permutation number", required=FALSE)
parser$add_argument("--model", nargs=1, help="R code that defines the voxelwise-model and any initialization", required=TRUE)


args <- parser$parse_args()


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

# reduce to vector and obtain list of nonzero vertices
mask.vector <- as.vector(mask)
mask.vertices <- which(mask.vector > 0)
nvertices <- length(mask.vertices)


###############################################################################
#### 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)

    })
}

if( ! exists('processVoxel')) {
    stop("The model file did not define the processVoxel function")
}


###############################################################################
#### Read in rds data
designmat <- readRDS(args$designmat)
attach(designmat)

voxelfile <- gsub(".nii.gz", ".rds", args$mask)
voxeldat <- readRDS(voxelfile)

stopifnot(dim(voxeldat)[1] == dim(designmat)[1])




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

# if permutation testing, we set the number of the permutation. This can be used as a seed or for other purposes
# we also drop the suffix (must be .nii.gz)
if (!is.null(args$permutationfile)) {
    permutationNumberString <- gsub("[^[:digit:]]", "",  args$permutationfile)
    permutationNumber <- as.numeric(permutationNumberString)

}

cat("Starting sequential job\n")
system.time(results <-sapply(1:nvertices, processVoxel))
prefix <- gsub(".nii.gz$","", args$mask)

# Use the permutation file name provided to modify permutation filenames if we are doing that
if (!is.null(args$permutationfile)) {
    # test to see if more than one name was returned
    # paste the permutation number to the array of names
    if(is.array(results)) {
        rownames(results) <- paste(attributes(results[,1])$names, permutationNumberString, sep=".")
    } else {
        names(results) <- rep(paste(attributes(results)$names[1], permutationNumberString, sep="."), length(results))
    }
}

npointWriteOutputFiles(prefix,results,mask)