# runs matrix decomposition with MCMC using rjags
 RunGAPS <- function(data, unc, numPatterns,
                     MaxAtomsA, alphaA, lambdaA, 
                     MaxAtomsP, alphaP, lambdaP, 
                     SAIter, iter, thin=1) {
 
 
   # get the data needed for MCMC
   matrixDecompData <- GetDataListGAPS(data, unc, numPatterns,
                                       MaxAtomsA, alphaA, lambdaA,
                                       MaxAtomsP, alphaP, lambdaP,
                                       SAIter, thin)
   matrixDecompInits <- list(A = matrix(
                               data=0, matrixDecompData$numGenes, numPatterns),
                             P = matrix(
                               data=0, numPatterns, matrixDecompData$numSamples))
 
   # create the model file used by this decomposition
   FID <- sprintf('%.10f', thin - floor(thin))
   ModelFile <- paste("GAPS", FID, "bug", sep=".")
 
   cat("var\n", file=ModelFile)
   cat("\tA[numGenes, numPatterns], ADim[numGenes, numPatterns],\n",
       file=ModelFile, append=TRUE)
   cat("\tP[numPatterns, numSamples], PDim[numPatterns, numSamples],\n",
       file=ModelFile, append=TRUE)
   cat("\tD[numGenes, numSamples], S[numGenes, numSamples];\n",
       file=ModelFile, append=TRUE)
   cat("\nmodel {\n", file=ModelFile, append=TRUE)
   cat("A[,] ~ datomicgibbs(ADim, MaxAtomsA, alphaA, lambdaA, SAIterA, fileID)\n", file=ModelFile, append=TRUE)
   cat("P[,] ~ datomicgibbs(PDim, MaxAtomsP, alphaP, lambdaP, SAIterP, fileID)\n", file=ModelFile, append=TRUE)
   cat("D[,] ~ gapsnorm(A,P,S)\n", file=ModelFile, append=TRUE)
   cat("}\n", file=ModelFile, append=TRUE) 
 
   # set the model being used
   matrixDecompModel <- jags.model(file=ModelFile,
                                   data=matrixDecompData,
                                   inits=matrixDecompInits,
                                   n.adapt=SAIter)
   # update and specify which nodes are monitored
   matrixDecompOutput <- jags.samples(model=matrixDecompModel,
                                      variable.names=c('A','P'),  
                                      n.iter=iter,
                                      thin=iter)
 
   file.remove(ModelFile)
   
   return(matrixDecompOutput)
   
 
 }