# run the full algorithm
CoGAPS <- function(data, unc, GStoGenes,
                 outputDir, outputBase="",
                 sep = "\t", isPercentError=FALSE,
                 numPatterns,
                 MaxAtomsA=2^32, alphaA=0.01,
                 MaxAtomsP=2^32, alphaP=0.01,
                 SAIter=1000000000, iter = 500000000, thin=-1,
                 nPerm=500, verbose=TRUE, plot=FALSE, keepChain=FALSE) {


  # decompose the data
  matrixDecomp <- GAPS(data=data, unc=unc,
                       outputDir=outputDir, outputBase=outputBase, 
		       sep=sep, isPercentError=isPercentError,
                       numPatterns=numPatterns, MaxAtomsA=MaxAtomsA, 
		       alphaA=alphaA, MaxAtomsP=MaxAtomsP,
                       alphaP=alphaP, SAIter=SAIter, iter=iter, thin=thin,
                       verbose=verbose, keepChain=keepChain)

  # plot patterns and show heatmap of Anorm
  if (plot) {
    plotGAPS(matrixDecomp$Amean, matrixDecomp$Pmean)
  }

  # compute the gene set scores
  GSP <- calcCoGAPSStat(matrixDecomp$Amean, matrixDecomp$Asd, GStoGenes, nPerm)

  # output gene set results
  if (outputDir != "") {
     outputGSPBase <- paste(outputDir, outputBase, sep="/")
  } else {
     outputGSPBase <- outputBase
  }
  if (outputBase != "") {
    outputGSPBase <- paste(outputGSPBase, "_", sep="")
  }
  
  UpFile <- paste(outputGSPBase, "GSUpStat.txt", sep="")
  DownFile <- paste(outputGSPBase, "GSDownStat.txt", sep="")
  ActEstFile <- paste(outputGSPBase, "GSActEst.txt", sep="")

  write.table(GSP$GSUpreg, file=UpFile, sep=sep)
  write.table(GSP$GSDownreg, file=DownFile, sep=sep)
  write.table(GSP$GSActEst, file=ActEstFile, sep=sep)

  return(list(meanChi2=matrixDecomp$meanChi2,
              D=matrixDecomp$D, Sigma=matrixDecomp$Sigma,
              Amean=matrixDecomp$Amean, Asd=matrixDecomp$Asd,
              Pmean=matrixDecomp$Pmean, Psd=matrixDecomp$Psd,
              meanMock=matrixDecomp$meanMock,
              GSUpreg=GSP$GSUpreg, GSDownreg=GSP$GSDownreg, GSActEst=GSP$GSActEst))
  
  
}