#' Simulation Scenario from Bhatnagar et al. (2018+) ggmix paper

#'

#' @description Function that generates data of the different simulation studies

#'   presented in the accompanying paper. This function requires the

#'   \code{popkin} and \code{bnpsd} package to be installed.

#' @param n number of observations to simulate

#' @param p_design number of variables in X_test, i.e., the design matrix

#' @param p_kinship number of variable in X_kinship, i.e., matrix used to

#'   calculate kinship

#' @param k number of intermediate subpopulations.

#' @param s the desired bias coefficient, which specifies sigma indirectly.

#'   Required if sigma is missing

#' @param Fst The desired final FST of the admixed individuals. Required if

#'   sigma is missing

#' @param b0 the true intercept parameter

#' @param eta the true eta parameter, which has to be \code{0 < eta < 1}

#' @param sigma2 the true sigma2 parameter

#' @param nPC number of principal components to include in the design matrix

#'   used for regression adjustment for population structure via principal

#'   components. This matrix is used as the input in a standard lasso regression

#'   routine, where there are no random effects.

#' @param geography the type of geography for simulation the kinship matrix.

#'   "ind" is independent populations where every individuals is actually

#'   unadmixed, "1d" is a 1D geography and "circ" is circular geography.

#'   Default: "ind". See the functions in the \code{bnpsd} for details on how

#'   this data is actually generated.

#' @param percent_causal percentage of \code{p_design} that is causal. must be

#'   \eqn{0 \leq percent_causal \leq 1}. The true regression coefficients are

#'   generated from a standard normal distribution.

#' @param percent_overlap this represents the percentage of causal SNPs that

#'   will also be included in the calculation of the kinship matrix

#' @param train_tune_test the proportion of sample size used for training tuning

#'   parameter selection and testing. default is 60/20/20 split

#' @details The kinship is estimated using the \code{popkin} function from the

#'   \code{popkin} package. This function will multiple that kinship matrix by 2

#'   to give the expected covariance matrix which is subsequently used in the

#'   linear mixed models

#' @return A list with the following elements \describe{\item{ytrain}{simulated

#'   response vector for training set} \item{ytune}{simulated response vector

#'   for tuning parameter selection set} \item{ytest}{simulated response vector

#'   for test set} \item{xtrain}{simulated design matrix for training

#'   set}\item{xtune}{simulated design matrix for tuning parameter selection

#'   set}\item{xtest}{simulated design matrix for testing set}

#'   \item{xtrain_lasso}{simulated design matrix for training set for lasso

#'   model. This is the same as xtrain, but also includes the nPC principal

#'   components} \item{xtune_lasso}{simulated design matrix for tuning parameter

#'   selection set for lasso model. This is the same as xtune, but also includes

#'   the nPC principal components}\item{xtest}{simulated design matrix for

#'   testing set for lasso model. This is the same as xtest, but also includes

#'   the nPC principal components} \item{causal}{character vector of the names

#'   of the causal SNPs} \item{beta}{the vector of true regression coefficients}

#'   \item{kin_train}{2 times the estimated kinship for the training set

#'   individuals} \item{kin_tune_train}{The covariance matrix between the tuning

#'   set and the training set individuals} \item{kin_test_train}{The covariance

#'   matrix between the test set and training set individuals}

#'   \item{Xkinship}{the matrix of SNPs used to estimate the kinship matrix}

#'   \item{not_causal}{character vector of the non-causal SNPs} \item{PC}{the

#'   principal components for population structure adjustment} }

#' @seealso \code{\link[bnpsd]{admix_prop_1d_linear}}

#' @export

#' @examples

#' admixed <- gen_structured_model(n = 100,

#'                                 p_design = 50,

#'                                 p_kinship = 5e2,

#'                                 geography = "1d",

#'                                 percent_causal = 0.10,

#'                                 percent_overlap = "100",

#'                                 k = 5, s = 0.5, Fst = 0.1,

#'                                 b0 = 0, nPC = 10,

#'                                 eta = 0.1, sigma2 = 1,

#'                                 train_tune_test = c(0.8, 0.1, 0.1))

#' names(admixed)

gen_structured_model <- function(n, p_design, p_kinship, k, s, Fst, b0, nPC = 10,

                                 eta, sigma2, geography = c("ind", "1d", "circ"),

                                 percent_causal, percent_overlap, train_tune_test = c(0.6, 0.2, 0.2)) {

  # p_design:

  # p_kinship:

  # k:	N

  # s:

  # F: The length-k vector of inbreeding coefficients (or FST's) of the intermediate subpopulations,

  # up to a scaling factor (which cancels out in calculations). Required if sigma is missing

  # Fst: T

  # browser()

  # define population structure

    )

  }

  # if (!requireNamespace("simtrait", quietly = TRUE)) {

  #   stop(strwrap("Package \"simtrait\" needed to simulate data.

  #                Please install it."),

  #        call. = FALSE

  #   )

  # }

    )

  }

  ########################

  ### Generate Kinship ###

  ########################

  # FF <- 1:k # subpopulation FST vector, up to a scalar

  # s <- 0.5 # desired bias coefficient

  # Fst <- 0.1 # desired FST for the admixed individuals

    # Q <- obj$Q

    # FF <- obj$F

    # rescaled inbreeding vector for intermediate subpopulations

    # get pop structure parameters of the admixed individuals

    # here’s the labels (for simplicity, list all individuals of S1 first, then S2, then S3)

    # data dimensions infered from labs:

    # desired admixture matrix ("is" stands for "Independent Subpopulations")

    # number of subpopulations "k_subpops"

    # desired admixture matrix

    # subpopulation FST vector, unnormalized so far

    # normalized to have the desired FST

    # NOTE fst is a function in the `popkin` package

    # verify FST for the intermediate subpopulations

    # fst(inbr_subpops)

    #> [1] 0.2

    # get coancestry of the admixed individuals

    # before getting FST for individuals, weigh then inversely proportional to subpop sizes

    # obj <- bnpsd::admix_prop_1d_circular(n_ind = n, k_subpops = k, s = s, F = FF, Fst = Fst)

    # Q <- obj$Q

    # FF <- obj$F

    # admixture proportions from *circular* 1D geography

    )

    # get pop structure parameters of the admixed individuals

  }

  #####################

  ### GENERATE SNPS ###

  #####################

  # browser()

    # this contains all SNPs (X_{Design}:X_{kinship})

    # out <- bnpsd::rbnpsd(Q = Q, F = FF, m = total_snps_to_simulate)

    # draw all random Allele Freqs (AFs) and genotypes

    # reuse the previous inbr_subpops, admix_proportions

      # NOTE by default p_subpops and p_ind are not returned, but here we will ask for them

    )

    ###################

    ### CAUSAL LOCI ###

    ###################

    # browser()

    # Snps used for kinship

    # all causal snps are in kinship matrix

      # browser()

      # compute marginal allele frequencies

      # select random SNPs! this performs the magic...

      # also runs additional checks

      # draw random SNP coefficients for selected loci

      # causal_coeffs <- stats::rnorm(ncausal, 0, 1)

      # subset data to consider causal loci only

    }

    # Xdesign_causal <- Xall[,causal_indexes,drop=F] # the subset of causal data (keep as a matrix even if m_causal == 1)

    # now estimate kinship using popkin

    # PhiHat <- popkin::popkin(Xkinship, lociOnCols = TRUE)

    # this contains all SNPs (X_{Testing}:X_{kinship})

    # out <- bnpsd::rbnpsd(Q = Q, F = FF, m = total_snps_to_simulate)

    # draw all random Allele Freqs (AFs) and genotypes

    # reuse the previous inbr_subpops, admix_proportions

      # NOTE by default p_subpops and p_ind are not returned, but here we will ask for them

    )

    # Snps used for kinship

    # length(snps_design)

    # setdiff(cnames, snps_kinships) %>% length()

      # compute marginal allele frequencies

      # select random SNPs! this performs the magic...

      # also runs additional checks

      # draw random SNP coefficients for selected loci

      # causal_coeffs <- stats::rnorm(ncausal, 0, 1)

      # causal <- sample(snps_design, ncausal, replace = FALSE)

    }

    # now estimate kinship using popkin

    # PhiHat <- popkin::popkin(Xkinship, lociOnCols = TRUE)

  }

  }

  # y <- trait

  }

  # y <- mu + sigma * matrix(rnorm(nsim * n), n, nsim)

  # y <- b0 + mu + t(P) + t(E)

  # y <- MASS::mvrnorm(1, mu = mu, Sigma = eta * sigma2 * kin + (1 - eta) * sigma2 * diag(n))

  # partition the data into train/tune/test

  ))

  # g %>% table

  # res = split(admixed$x, g)

  # xtrain <- Xdesign[ind,,drop=FALSE]

  # xtune <- Xdesign[-ind,,drop=FALSE]

              # used in manuscript to generate simulation study figures

  ))

}

"%ni%" <- Negate("%in%")

# internal function

# taken verbatim from glmnet package.

# it used to be exported by glmnet, but is no longer exported.

lambda.interp <- function (lambda, s) {

  }

  else {

  }

}

# internal function

# taken verbatim from glmnet package.

# it used to be exported by glmnet, but is no longer exported.

nonzeroCoef <- function (beta, bystep = FALSE) {

    else {

    }

  }

  else {

      }

      else {

      }

    }

  }

}

# internal function

# taken verbatim from https://github.com/OchoaLab/simtrait/blob/master/R/select_loci.R

select_loci <- function(maf, m_causal, maf_cut = 0.05) {

  # check for missing parameters

  # data dimensions

  # other checks

  # select random loci!

  # we might not want to pick extremely rare alleles, so set MAF thresholds

}