require(mcgraph)
require(huge)
source("huge.wrapper.R")

# Data generation default options
mc.settings <- formals(mcg.graph2data)
hg.settings <- formals(huge)

do.comparison <- function(n=200, nruns=20, nodes=c(40, 80, 400), nw.types=c("random", "scale-free", "hub", "cluster"),
                              pkg=c("huge", "mcgraph"), hg.cls=c("mb", "ct", "glasso"), mcg.cls=c("mcg.lvs", "mcg.ct", "mcg.rpart")) {
    # Save function call
    cll <- match.call()

    # Timestamp of data
    ts <- format(Sys.time(), "%Y-%d-%m-%X")

    # Allocate list for results
    res <- list()
    nms.res <- c("pkg", "nw.type", "classifier", "AUC", "PR.AUC", "BCR", "MCC", "norm_MCC", "nodes", "edges", "components", "seed", "sparsity")
    for (l in seq_along(nms.res)) {
        res[[l]] <- rep(0, length(pkg) * length(nodes) * length(nw.types) * nruns)
    }
    names(res) <- nms.res

    # Pool classifier
    classifier <- c(hg.cls, mcg.cls)

    # Main loops
    i <- 1
    for (nw in nw.types) {
        for (nd in nodes) {
            for (rn in 1:nruns) {
                print(nd)
                data.pkgs = list()
                # Create a 9-digit sampled seed for each iteration based on the numbers 1 to 9 without replacement
                # This gives (n)! unique possibilites for a n-digit integer, i.e. in this case (9)! = 362880
                #dig <- sample(2:9, 1)
                #seed <- strtoi(paste0(sample(1:9, dig), collapse=""))
                # huge data
                #set.seed(seed)
                hg <- huge.generator(n=n, d=nd, graph=nw)
                # Extract relevant objects from huge list object
                data.pkgs$huge <- hg$data
                colnames(data.pkgs$huge) <- mcg.autonames(LETTERS, ncol(data.pkgs$huge))
                # Monte Carlo data
                #set.seed(seed)
                data.pkgs$mcgraph <- t(mcg.graph2data(A=mcg.new(hg$theta)))
                # Estimators of huge
                for (cl in classifier) {
                    for (pk in pkg) {
                        if (cl %in% hg.cls) {
                            est <- huge(data.pkgs[[pk]], method=cl)
                            roc <- mcg.roc(d=est, G=mcg.new(hg$theta), thresholds=rev(est$lambda), fun=huge.wrapper)
                        } else if (cl %in% mcg.cls) {
                            roc <- mcg.roc(d=data.pkgs[[pk]], G=mcg.new(hg$theta), thresholds=seq(0, 0.1, length.out=11), fun=cl)
                        } else {
                            stop("Error. Unknown classifier")
                        }
                        res$pkg[[i]] <- pk
                        res$nw.type[[i]] <- nw
                        res$classifier[[i]] <- if (cl %in% hg.cls) { paste0("huge.", cl) } else { cl }
                        res$AUC[[i]] <- roc$auc
                        res$PR.AUC[[i]] <- roc$pr.auc
                        res$BCR[[i]] <- mean(roc$bcr, na.rm=TRUE)
                        res$MCC[[i]] <- mean(roc$mcc, na.rm=TRUE)
                        res$norm_MCC[[i]] <- mean(roc$mcc.norm, na.rm=TRUE)
                        res$nodes[[i]] <- nd
                        res$edges[[i]] <- length(which(as.matrix(hg$theta) != 0)) / 2
                        res$components[[i]] <- length(unique(mcg.components(mcg.new(hg$theta))))
                        #res$seed[[i]] <- seed
                        res$seed[[i]] <- NA
                        res$sparsity[[i]] <- round(hg$sparsity, 3)
                        i <- i + 1
                    }
                }
            }
        }
    }

    # Convert into data frame
    final <- list(data=data.frame(pkg=res$pkg, nw.type=res$nw.type, classifier=res$classifier, AUC=res$AUC, PR.AUC=res$PR.AUC,
                                  BCR=res$BCR, MCC=res$MCC, norm_MCC=res$norm_MCC, nodes=res$nodes, edges=res$edges,
                                  components=res$components, seed=res$seed, sparsity=res$sparsity),
                  timestamp=ts, call=cll, n=n)

    # If the 'results' directory does not exist and the user is not in it, create it.
    if (dir.exists("results") == FALSE & length(grep(x=getwd(), pattern="results") == 0)) dir.create("results")

    # Change to results directory
    setwd("results")

    # Filename
    fn <- paste0("results_", final$timestamp)

    # Save data as R object and as tab file with date and time
    saveRDS(final, fn)
    write.table(final$data, paste0(fn, ".tab"), sep="\t", row.names=FALSE)

    # Change to inital directory again
    setwd("..")

    return(final)
}

do.comparison(n=200, nodes=c(20, 40, 80, 400), nruns=40)