#' AMARETTO_LarsenBased

#' 

#' @param Data 

#' @param Clusters 

#' @param RegulatorData 

#' @param Parameters 

#' @param NrCores 

#' @param random_seeds 

#' @param convergence_cutoff 

#'

#' @import MultiAssayExperiment

#' @import graphics

#' @importFrom Matrix nnzero

#' @importFrom doParallel registerDoParallel

#' 

#' @return result

#' @keywords internal

AMARETTO_LarsenBased <- function(Data, Clusters, RegulatorData, Parameters, NrCores, random_seeds, convergence_cutoff) {

    registerDoParallel(cores = NrCores)

    ptm1 <- proc.time()

    RegulatorData_rownames = rownames(RegulatorData)

    Data_rownames = rownames(Data)

    AutoRegulation = Parameters$AutoRegulation

    RegulatorSign = array(0, length(RegulatorData_rownames))

    Lambda = Parameters$Lambda2

    OneRunStop = Parameters$OneRunStop

    random_seeds = Parameters$random_seeds

    convergence_cutoff = Parameters$convergence_cutoff

    if (AutoRegulation == 1) {

        cat("\tAutoregulation is turned ON.\n")

    } else if (AutoRegulation == 2) {

        cat("\tAutoregulation is turned ON.\n")

    } else {

        cat("\tAutoregulation is turned OFF.\n")

    }

    NrReassignGenes = length(Data_rownames)

    NrReassignGenes_history <- NrReassignGenes

    error_history<-list()

    index=1

    while (NrReassignGenes > convergence_cutoff * length(Data_rownames)) {

        ptm <- proc.time()

        switch(Parameters$Mode, larsen = {

            regulatoryPrograms <- AMARETTO_LearnRegulatoryProgramsLarsen(Data, Clusters, RegulatorData, RegulatorSign, Lambda, AutoRegulation, alpha = Parameters$alpha, pmax = Parameters$pmax, random_seeds)

        })

        error_history[[index]]<-(rowMeans(regulatoryPrograms$error * regulatoryPrograms$error))^0.5

        index<-index+1

        ptm <- proc.time() - ptm

        printf("Elapsed time is %f seconds\n", ptm[3])

        NrClusters = length(unique(Clusters))

        sum = 0

        for (i in 1:NrClusters) {

            sum = sum + Matrix::nnzero(regulatoryPrograms$Beta[i,])

        }

        avg = sum/NrClusters

        printf("Average nr of regulators per module: %f \n", avg)

        PreviousClusters = Clusters

        if (OneRunStop == 1) {

            break

        }

        ptm <- proc.time()

        ReassignGenesToClusters <- AMARETTO_ReassignGenesToClusters(Data, RegulatorData, regulatoryPrograms$Beta, Clusters, AutoRegulation)

        ptm <- proc.time() - ptm

        printf("Elapsed time is %f seconds\n", ptm[3])

        NrReassignGenes = ReassignGenesToClusters$NrReassignGenes

        Clusters = ReassignGenesToClusters$Clusters

        printf("Nr of reassignments is: %i \n", NrReassignGenes)

        NrReassignGenes_history <- c(NrReassignGenes_history, NrReassignGenes)

    }

    ptm1 <- proc.time() - ptm1

    printf("Elapsed time is %f seconds\n", ptm1[3])

    ModuleMembership = as.matrix(PreviousClusters)

    rownames(ModuleMembership) = rownames(Data)

    colnames(ModuleMembership) = c("ModuleNr")

    result <- list(NrModules = length(unique(Clusters)), RegulatoryPrograms = regulatoryPrograms$Beta, AllRegulators = rownames(RegulatorData),

                   AllGenes = rownames(Data), ModuleMembership = ModuleMembership, AutoRegulationReport = regulatoryPrograms$AutoRegulationReport,

                   run_history = list(NrReassignGenes_history = NrReassignGenes_history, error_history = error_history))

    return(result)

}

#' AMARETTO_LearnRegulatoryProgramsLarsen

#' @importFrom foreach foreach

#' @importFrom glmnet cv.glmnet

#' @importFrom stats coef

#' @return result

#' @keywords internal

AMARETTO_LearnRegulatoryProgramsLarsen <- function(Data, Clusters, RegulatorData, RegulatorSign, Lambda, AutoRegulation, alpha, pmax, random_seeds) {

    `%dopar%` <- foreach::`%dopar%`

    RegulatorData_rownames = rownames(RegulatorData)

    Data_rownames = rownames(Data)

    trace = 0

    NrFolds = 10

    NrClusters = length(unique(Clusters))

    NrGenes = nrow(Data)

    NrSamples = ncol(Data)

    NrInterpolateSteps = 100

    if (AutoRegulation >= 1) {}

    else if (AutoRegulation == 0) {

        BetaSpecial = list(NrClusters, 1)

        RegulatorPositions = list(NrClusters, 1)

    }

    y_all = mat.or.vec(NrClusters, NrSamples)

    ClusterIDs = unique(Clusters)

    ClusterIDs = sort(ClusterIDs, decreasing = FALSE)

    cnt <- 1:NrClusters

    ptm1 <- proc.time()

    BetaY_all <- foreach(i = 1:NrClusters, .combine = cbind, .init = list(list(), list(), list()), .packages = "glmnet") %dopar% {

            if (length(which(Clusters == ClusterIDs[i])) > 1) {

                y = apply((Data[which(Clusters == ClusterIDs[i]),]), 2, mean)

            } else {

                y = Data[which(Clusters == ClusterIDs[i]),]

            }

            CurrentClusterPositions = which(Clusters %in% ClusterIDs[i])

            nrGenesInClusters = length(CurrentClusterPositions)

            if (AutoRegulation >= 1) {

                X = RegulatorData

            } else if (AutoRegulation == 0) {

                X = RegulatorData[setdiff(RegulatorData_rownames, Data_rownames[CurrentClusterPositions]),]

            }

            suppressWM = function(...) suppressWarnings(suppressMessages(...))

            if(!is.null(random_seeds)){

              seed<-ifelse (length(random_seeds)==2,random_seeds[2],random_seeds[1])

              set.seed(seed)

            }

            fit = suppressWM(cv.glmnet(t(X), y, alpha = alpha, pmax = pmax, lambda = Lambda_Sequence(t(X), y)))

            nonZeroLambdas <- fit$lambda[which(fit$nzero > 0)]

            nonZeroCVMs <- fit$cvm[which(fit$nzero > 0)]

            if (length(which(nonZeroCVMs == min(nonZeroCVMs, na.rm = TRUE))) == 0) {

                warnMessage <- paste0("\nOn cluster ", i, " there were no cv.glm results that gave non-zero coefficients.")

                message(warnMessage)

            }

            bestNonZeroLambda <- nonZeroLambdas[which(nonZeroCVMs == min(nonZeroCVMs, na.rm = TRUE))]

            b_o = coef(fit, s = bestNonZeroLambda)

            b_opt <- c(b_o[2:length(b_o)])

            if (AutoRegulation == 2) {

                CurrentUsedRegulators = RegulatorData_rownames[which(b_opt != 0, arr.ind = TRUE)]

                CurrentClusterMembers = Data_rownames[CurrentClusterPositions]

                nrIterations = 0

                while (length(CurrentClusterMembers[CurrentClusterMembers %in% CurrentUsedRegulators]) != 0) {

                  CurrentClusterMembers = setdiff(CurrentClusterMembers, CurrentUsedRegulators)

                  nrCurrentClusterMembers = length(CurrentClusterMembers)

                  if (nrCurrentClusterMembers > 0) {

                    names = Data_rownames %in% CurrentClusterMembers

                    if (length(which(names == TRUE)) > 1) {

                      y = apply((Data[names, ]), 2, mean)

                    } else {

                      y = Data[names, ]

                    }

                    if(!is.null(random_seeds)){

                      seed<-ifelse (length(random_seeds)==2,random_seeds[2],random_seeds[1])

                      set.seed(seed)

                    }

                    fit = suppressWM(cv.glmnet(t(X), y, alpha = alpha, pmax = pmax, lambda = Lambda_Sequence(t(X), y)))

                    nonZeroLambdas <- fit$lambda[which(fit$nzero > 0)]

                    nonZeroCVMs <- fit$cvm[which(fit$nzero > 0)]

                    if (length(which(nonZeroCVMs == min(nonZeroCVMs, na.rm = TRUE))) == 0) {

                      warnMessage <- paste0("\nOn cluster ", i, " there were no cv.glm results that gave non-zero coefficients during the Autoregulation step.")

                      message(warnMessage)

                    }

                    bestNonZeroLambda <- nonZeroLambdas[which(nonZeroCVMs == min(nonZeroCVMs, na.rm = TRUE))]

                    new_b_o = coef(fit, s = bestNonZeroLambda)

                    new_b_opt <- c(new_b_o[2:length(b_o)])

                    CurrentUsedRegulators = RegulatorData_rownames[which(new_b_opt != 0)]

                    nrIterations = nrIterations + 1

                    b_opt = new_b_opt

                  } else {

                    b_opt = rep(0, length(RegulatorData_rownames))

                  }

                }

                Report <- c(length(CurrentClusterPositions), length(CurrentClusterMembers), nrIterations)

            }

            if (sum(RegulatorSign[which(RegulatorSign != 0)]) > 0) {

                RegulatorCheck = RegulatorSign * t(b_opt)

                WrongRegulators = which(RegulatorCheck < 0)

                if (length(WrongRegulators) == 0) {

                  b_opt[WrongRegulators] = 0

                }

            }

            if (AutoRegulation >= 1) {

            } else {

                BetaSpecial[i] = b_opt

                RegulatorPositions[i] = (RegulatorData_rownames %in% setdiff(RegulatorData_rownames, Data_rownames[CurrentClusterPositions]))

            }

            list(b_opt, y, Report)

        }

    if (AutoRegulation == 0) {

        for (i in 1:NrClusters) {

            Beta[i, RegulatorPositions[i]] = BetaSpecial[i]

        }

    }

    tmpPos = NrClusters + 1

    Beta <- do.call(cbind, BetaY_all[1, 2:tmpPos])

    Beta = t(Beta)

    colnames(Beta) = RegulatorData_rownames

    rownames(Beta) = gsub("result.", "Module_", rownames(Beta))

    y_all <- do.call(cbind, BetaY_all[2, 2:tmpPos])

    y_all = t(y_all)

    rownames(y_all) = gsub("result.", "Module_", rownames(y_all))

    AutoRegulationReport <- do.call(cbind, BetaY_all[3, 2:tmpPos])

    AutoRegulationReport = t(AutoRegulationReport)

    rownames(AutoRegulationReport) = gsub("result.", "Module_", rownames(AutoRegulationReport))

    error = y_all - (Beta %*% RegulatorData)

    result <- list(Beta = Beta, error = error, AutoRegulationReport = AutoRegulationReport)

    return(result)

}

#' AMARETTO_ReassignGenesToClusters

#'

#' @return result

#' @importFrom Matrix nnzero

#' @importFrom stats cor

#' @keywords internal

AMARETTO_ReassignGenesToClusters <- function(Data, RegulatorData, Beta, Clusters, AutoRegulation) {

    `%dopar%` <- foreach::`%dopar%`

    RegulatorData_rownames = rownames(RegulatorData)

    Data_rownames = rownames(Data)

    NrGenes = nrow(Data)

    NrSamples = ncol(Data)

    NrReassignGenes = 0

    X = RegulatorData

    X1 = data.matrix(X)

    ModuleVectors = Beta %*% X1

    GeneNames = rownames(Data)

    ptm1 <- proc.time()

    i <- NULL

    nc <- foreach(i = 1:NrGenes, .combine = c) %dopar%

        {

            OldModule = Clusters[i]

            CurrentGeneVector = Data[i, , drop = FALSE]

            Correlations = cor(t(CurrentGeneVector), t(ModuleVectors))

            corr = data.matrix(Correlations, rownames.force = NA)

            MaxCorrelation = max(corr, na.rm = TRUE)

            MaxPosition = which(signif(corr, digits = 7) == signif(MaxCorrelation, digits = 7))

            MaxPosition = MaxPosition[1]

            if (AutoRegulation > 0) {

                if (MaxPosition != OldModule) {

                  NrReassignGenes = NrReassignGenes + 1

                }

                NewClusters = MaxPosition

            } else {

                if (nnzero(rownames(RegulatorData_rownames) %in% GeneNames[i]) != 0) {

                  if (nnzero(which(which(GeneNames %in% rownames(RegulatorData_rownames)) %in% i) %in% which(Beta[MaxPosition,] != 0)) != 0) {

                    if (MaxPosition != OldModule) {

                      NrReassignGenes = NrReassignGenes + 1

                    }

                    NewClusters = MaxPosition

                  } else {

                    NewClusters = OldModule

                  }

                } else {

                  if (MaxPosition != OldModule) {

                    NrReassignGenes = NrReassignGenes + 1

                  }

                  NewClusters = MaxPosition

                }

            }

        }

    ptm1 <- proc.time() - ptm1

    NrReassignGenes = length(which(nc != Clusters))

    result <- list(NrReassignGenes = NrReassignGenes, Clusters = nc)

    return(result)

}

#' AMARETTO_CreateModuleData

#'

#' @param AMARETTOinit List output from AMARETTO_Initialize().

#' @param AMARETTOresults List output from AMARETTO_Run()

#'

#' @return result

#' @export

#' @examples

#' data('ProcessedDataLIHC')

#' AMARETTOinit <- AMARETTO_Initialize(ProcessedData = ProcessedDataLIHC,

#'                                     NrModules = 2, VarPercentage = 50)

#' AMARETTOresults <- AMARETTO_Run(AMARETTOinit)

#' AMARETTO_MD <- AMARETTO_CreateModuleData(AMARETTOinit, AMARETTOresults)

AMARETTO_CreateModuleData <- function(AMARETTOinit, AMARETTOresults) {

    ModuleData = matrix(0, AMARETTOresults$NrModules, length(colnames(AMARETTOinit$MA_matrix_Var)))

    rownames(ModuleData) = rownames(AMARETTOresults$AutoRegulationReport)

    colnames(ModuleData) = colnames(AMARETTOinit$MA_matrix_Var)

    for (ModuleNr in 1:AMARETTOresults$NrModules) {

        currentModuleData = AMARETTOinit$MA_matrix_Var[AMARETTOresults$ModuleMembership[, 1] == ModuleNr, ]

        if (length(which(AMARETTOresults$ModuleMembership[, 1] == ModuleNr)) > 1) {

            ModuleData[ModuleNr, ] = colMeans(currentModuleData)

        } else {

            ModuleData[ModuleNr, ] = currentModuleData

        }

    }

    return(ModuleData)

}

#' AMARETTO_CreateRegulatorPrograms

#'

#' @param AMARETTOinit  List output from AMARETTO_Initialize().

#' @param AMARETTOresults List output from AMARETTO_Run()

#'

#' @return result

#' @export

#' @examples

#' data('ProcessedDataLIHC')

#' AMARETTOinit <- AMARETTO_Initialize(ProcessedData = ProcessedDataLIHC,

#'                                     NrModules = 2, VarPercentage = 50)

#' AMARETTOresults <- AMARETTO_Run(AMARETTOinit)

#' AMARETTO_RP <- AMARETTO_CreateRegulatorPrograms(AMARETTOinit,AMARETTOresults)

AMARETTO_CreateRegulatorPrograms <- function(AMARETTOinit, AMARETTOresults) {

    RegulatorProgramData = matrix(0, AMARETTOresults$NrModules, length(colnames(AMARETTOinit$MA_matrix_Var)))

    rownames(RegulatorProgramData) = rownames(AMARETTOresults$AutoRegulationReport)

    colnames(RegulatorProgramData) = colnames(AMARETTOinit$MA_matrix_Var)

    RegulatorNames = rownames(AMARETTOinit$RegulatorData)

    for (ModuleNr in 1:AMARETTOresults$NrModules) {

        currentRegulators = RegulatorNames[which(AMARETTOresults$RegulatoryPrograms[ModuleNr, ] != 0)]

        weights = AMARETTOresults$RegulatoryPrograms[ModuleNr, currentRegulators]

        RegulatorData = AMARETTOinit$RegulatorData[currentRegulators, ]

        RegulatorProgramData[ModuleNr, ] = weights %*% RegulatorData

    }

    return(RegulatorProgramData)

}

#' Lambda_Sequence

#'

#' @return result

#' @keywords internal

Lambda_Sequence <- function(sx, sy) {

    n <- nrow(sx)

    lambda_max <- max(abs(colSums(sx * sy)))/n

    epsilon <- 1e-04

    K <- 100

    lambdaseq <- round(exp(seq(log(lambda_max), log(lambda_max * epsilon), length.out = K)), digits = 10)

    return(lambdaseq)

}