# Author: Kosuke Hamazaki hamazaki@ut-biomet.org
# 2020 The University of Tokyo
#
# Description:
# Definition of simBsCma class
#' R6 Class Representing a Breeding Scheme
#'
#' @description
#' simBsCma object store specific information of simulation results of breeding scheme optimized by CMA-ES.
#'
# @details
# Details: simBsCma object store specific information of simulation results of breeding scheme optimized by CMA-ES.
#'
#' @export
#' @import R6
simBsCma <- R6::R6Class(
"simBsCma",
public = list(
#' @field simBsName [character] Name of this simulation of breeding schemes
simBsName = NULL,
#' @field bsInfoInit [bsInfo] breeding scheme info
#' (see:\link[myBreedSimulatR]{bsInfo})
bsInfoInit = NULL,
#' @field breederInfoInit [breederInfo] breeder info
#' (see:\link[myBreedSimulatR]{breederInfo})
breederInfoInit = NULL,
#' @field lociEffMethod [character] How to compute loci effects (use true QTL effects ("true") or estimated marker effects ("estimated"))
lociEffMethod = NULL,
#' @field trainingPopType [character] Training population consists of all populations created in the breeding scheme for `trainingPopType = "all"`,
#' or consists of the latest population in breeding scheme for `trainingPopType = "latest"`.
trainingPopType = NULL,
#' @field trainingPopInit [character / numeric] Training population names or No.s (not generations!!) for initial population
trainingPopInit = NULL,
#' @field trainingIndNamesInit [character] Names of training individuals for initia; population
trainingIndNamesInit = NULL,
#' @field methodMLRInit [character] Methods for estimating marker effects for initial population.
#' The following methods are offered:
#'
#' "Ridge", "LASSO", "ElasticNet", "RR-BLUP", "GBLUP", "BayesA" (uni-trait),
#' "BayesB" (uni-trait), "BayesC" (uni-trait), "BRR", "BL" (uni-trait), "SpikeSlab" (multi-trait)
methodMLRInit = NULL,
#' @field multiTraitInit [logical] Use multiple-trait model for estimation of marker effects or not for initial population
multiTraitInit = NULL,
#' @field nIterSimulation [numeric] Number of iterations for this simulation setting
nIterSimulation = NULL,
#' @field nGenerationProceed [numeric] Number of generations to be proceeded
nGenerationProceed = NULL,
#' @field nGenerationProceedSimulation [numeric] Number of generations to be proceeded in simulation for optimizing policy
nGenerationProceedSimulation = NULL,
#' @field setGoalAsFinalGeneration [logical] Set goal for simulation of breeding scheme as the real final generation
setGoalAsFinalGeneration = NULL,
#' @field performOptimization [logical] Perform optimization of policy in each generation or not
performOptimization = NULL,
#' @field useFirstOptimizedValue [logical] Perform optimization of policy only for the initial population and use the optimized hyperprameters permanently.
useFirstOptimizedValue = NULL,
#' @field performRobustOptimization [logical] Whether or not performing robust optimization.
performRobustOptimization = NULL,
#' @field lowerQuantile [logical] Lower quantile for the robust optimization (value at risk)
lowerQuantile = NULL,
#' @field sameAcrossGeneration [logical] Use same hyper prameter across generations or not
sameAcrossGeneration = NULL,
#' @field propSelOpt [logical] Optimize proportion of selected individuals or not
propSelOpt = NULL,
#' @field hMin [numeric] Lower bound of hyper parameters
hMin = NULL,
#' @field hMax [numeric] Upper bound of hyper parameters
hMax = NULL,
#' @field nTotalIterForOneOptimization [numeric] Number of total iterations that can be assigned for one optimization process
nTotalIterForOneOptimization = NULL,
#' @field nIterSimulationPerEvaluation [numeric] Number of simulations per one evaluation of the hyperparameter set of your interest
nIterSimulationPerEvaluation = NULL,
#' @field nIterSimulationForOneMrkEffect [numeric] Number of simulations per one evaluation of the hyperparameter set of your interest for one set of marker effects
nIterSimulationForOneMrkEffect = NULL,
#' @field nIterMrkEffectForRobustOptimization [numeric] Number of sets of marker effects utilized for the robust optimization
nIterMrkEffectForRobustOptimization = NULL,
#' @field nIterOptimization [numeric] Number of iterations required for one optimization
nIterOptimization = NULL,
#' @field lambda [numeric] Number of offspring in CMA-ES. Must be greater than or equal to mu.
lambda = NULL,
#' @field mu [numeric] Population size in CMA-ES.
mu = NULL,
#' @field saveObjNameBase [character] Base name of `cmaES` object to be saved
saveObjNameBase = NULL,
#' @field whenToSaveObj [numeric] When (how many iterations) to save `cmaES` object in CMA-ES
whenToSaveObj = NULL,
#' @field nTopEvalForOpt [numeric] Number of individuals to be evaluated when evaluating population max or population min for optimization of hyperparameters
nTopEvalForOpt = NULL,
#' @field rewardWeightVec [numeric] When returning reward function, `rewardWeightVec` will be multiplied by estimated GVs for each generation to evaluate the method.
#' If you want to apply discounted method, you can achieve by `rewardWeightVec = sapply(1:nGenerationProceedSimulation, function(genProceedNo) gamma ^ (genProceedNo - 1))` where `gamma` is discounted rate.
#' Or if you want to evaluate just the final generation, you can achieve by `rewardWeightVec = c(rep(0, nGenerationProceedSimulation - 1), 1)`.
rewardWeightVec = NULL,
#' @field digitsEval [numeric] When you evaluate each hyperparameter, you can round the average of evaluates (`eval`) with `round(eval, digitsEval)`.
digitsEval = NULL,
#' @field nRefreshMemoryEvery [numeric] Every `nRefreshMemoryEvery` iterations, we refresh memory used for simulations by `gc(reset = TRUE)`.
nRefreshMemoryEvery = NULL,
#' @field updateBreederInfo [logical] Update breederInfo or not for each generation
updateBreederInfo = NULL,
#' @field phenotypingInds [logical] Phenotyping individuals or not for each generation
phenotypingInds = NULL,
#' @field nRepForPhenoInit [numeric] Number of replications to be phenotyped for initial population
nRepForPhenoInit = NULL,
#' @field nRepForPheno [numeric] Number of replications to be phenotyped for each generation
nRepForPheno = NULL,
#' @field updateModels [logical] Whether or not updating the model for each generation
#' (If `phenotypingInds = FALSE` for generation of your interest, `updateModels` will be automatically `FALSE` for that generation.)
updateModels = NULL,
#' @field updateBreederInfoSimulation [logical] Update breederInfo or not for each generation in simulations for optimization
updateBreederInfoSimulation = NULL,
#' @field phenotypingIndsSimulation [logical] Phenotyping individuals or not for each generation in simulations for optimization
phenotypingIndsSimulation = NULL,
#' @field nRepForPhenoSimulation [numeric] Number of replications to be phenotyped for each generation in simulations for optimization
nRepForPhenoSimulation = NULL,
#' @field updateModelsSimulation [logical] Whether or not updating the model for each generation in simulations for optimization
#' (If `phenotypingInds = FALSE` for generation of your interest, `updateModels` will be automatically `FALSE` for that generation.)
updateModelsSimulation = NULL,
#' @field methodMLR [character] Methods for estimating marker effects.
#' The following methods are offered:
#'
#' "Ridge", "LASSO", "ElasticNet", "RR-BLUP", "GBLUP", "BayesA" (uni-trait),
#' "BayesB" (uni-trait), "BayesC" (uni-trait), "BRR", "BL" (uni-trait), "SpikeSlab" (multi-trait)
methodMLR = NULL,
#' @field multiTrait [logical] Use multiple-trait model for estimation of marker effects or not
multiTrait = NULL,
#' @field nSelectionWaysVec [numeric] Number of selection ways for each generation
nSelectionWaysVec = NULL,
#' @field selectionMethodList [list] (list of) Selection method for each generation
selectionMethodList = NULL,
#' @field traitNoSelList [list] (list of list of) Number of trait No of your interest for selection for each generation
traitNoSelList = NULL,
#' @field blockSplitMethod [character] How to determine the markers belonging to each block when computing OHV.
#' You can define the number of markers in each block (`nMrkInBlock`) or the minimum length of each segment (`minimumSegmentLength`).
blockSplitMethod = NULL,
#' @field nMrkInBlock [numeric] Number of markers in each block. This will be used for the computation of OHV.
nMrkInBlock = NULL,
#' @field minimumSegmentLength [numeric] Minimum length of each segment [cM]. This will be used for the computation of OHV.
minimumSegmentLength = NULL,
#' @field nSelInitOPVList [list] (list of) Number of selected candiates for first screening before selecting parent candidates by OPV for each generation
nSelInitOPVList = NULL,
#' @field nIterOPV [numeric] Number of iterations for computation of OPV
nIterOPV = NULL,
#' @field nProgeniesEMBVVec [numeric] Number of progenies of double haploids produced when computing EMBV for each generation
nProgeniesEMBVVec = NULL,
#' @field nIterEMBV [numeric] Number of iterations to estimate EMBV
nIterEMBV = NULL,
#' @field nCoresEMBV [numeric] Number of cores used for EMBV estimation
nCoresEMBV = NULL,
#' @field clusteringForSelList [list] (list of) Apply clustering results of marker genotype to selection or not for each generation
clusteringForSelList = NULL,
#' @field nClusterList [list] (list of) Number of clusters for each generation
nClusterList = NULL,
#' @field nTopClusterList [list] (list of) Number of top clusters used for selection for each generation
nTopClusterList = NULL,
#' @field nTopEachList [list] (list of) Number of selected individuals in each cluster for each generation
nTopEachList = NULL,
#' @field nSelList [list] (list of) Number of selection candidates for each generation
nSelList = NULL,
#' @field multiTraitsEvalMethodList [list] (list of) When evaluating multiple traits, you can choose how to evaluate these traits simultaneously.
#' One method is to take a weighted sum of these traits, and the other is to compute a product of these traits (adjusted to be positive values in advance.)
multiTraitsEvalMethodList = NULL,
#' @field hSelList [list] (list of) Hyperparameter which determines which trait is weighted when selecting parent candidates for multiple traits for each generation
hSelList = NULL,
#' @field matingMethodVec [character] Mating method for each generation
matingMethodVec = NULL,
#' @field allocateMethodVec [character] Allocation method for each generation
allocateMethodVec = NULL,
#' @field weightedAllocationMethodList [list] (list of) Which selection index will be used for weighted resource allocation for each generation
weightedAllocationMethodList = NULL,
#' @field traitNoRAList [list] (list of) Trait No of your interest for resource allocation for each generation
traitNoRAList = NULL,
#' @field minimumUnitAllocateVec [numeric] (vector of) Minimum number of allocated progenies for each pair.
minimumUnitAllocateVec = NULL,
#' @field includeGVPVec [logical] Whether or not to consider genetic variance of progenies of each pair when determining the number of progenies per each pair for each generation
includeGVPVec = NULL,
#' @field nNextPopVec [numeric] Number of progenies for the next generation for each generation
nNextPopVec = NULL,
#' @field targetGenGVPVec [numeric] Vector of target generation of selfing when evaluating the genetic variance of progenies (F#)
targetGenGVPVec = NULL,
#' @field performOCSVec [logical] Vector of whether or not performing OCS before allocation
performOCSVec = NULL,
#' @field targetGenOCSVec [numeric] Vector of target generation in OCS
targetGenOCSVec = NULL,
#' @field HeStarRatioVec [numeric] Vector of ratio of He0 to HeStar
HeStarRatioVec = NULL,
#' @field degreeOCSVec [numeric] Vector of how generation effect is emphasized when computing He(t)
degreeOCSVec = NULL,
#' @field includeGVPOCSVec [logical] Vector of whether or not to consider genetic variance of progenies of each pair when evaluating each pair in OCS
includeGVPOCSVec = NULL,
#' @field hOCSList [list] List of weight for each trait/criterion to evaluate each mating pair
hOCSList = NULL,
#' @field weightedAllocationMethodOCSList [list] List of which selection index will be used to evaluate mating pair in OCS
weightedAllocationMethodOCSList = NULL,
#' @field traitNoOCSList [list] List of trait No to evaluate mating pair in OCS
traitNoOCSList = NULL,
#' @field nCrossesOCSVec [numeric] Vector of number of crosses that will be selected in OCS
nCrossesOCSVec = NULL,
#' @field nameMethod [character] Method for naming individuals
nameMethod = NULL,
#' @field nCores [numeric] Number of cores used for simulations of breeding scheme
nCores = NULL,
#' @field nCoresPerOptimization [numeric] Number of cores used for simulations of breeding scheme when optimizing hyperparameters
nCoresPerOptimization = NULL,
#' @field overWriteRes [logical] Overwrite simulation results when the targeted results already exists
overWriteRes = NULL,
#' @field showProgress [logical] Show progress bar or not
showProgress = NULL,
#' @field returnMethod [character] Which type of results will be returned (saved) in the object
returnMethod = NULL,
#' @field saveAllResAt [character] If NULL, we won't save the simulation results. Else, we will save bsInfo and breederInfo for each iteration in the directory defined by `saveAllResAt` path.
saveAllResAt = NULL,
#' @field evaluateGVMethod [character] Which type of GV (true GV or estimated GV) will be used for evaluation of the simulation results.
evaluateGVMethod = NULL,
#' @field nTopEval [numeric] Number of individuals to be evaluated when evaluating population max or population min
nTopEval = NULL,
#' @field traitNoEval [numeric] Trait No of your interest for evaluation of the method
traitNoEval = NULL,
#' @field hEval [numeric] Hyperparameter which determines which BV is emphasized when evaluating simulation results of each method
hEval = NULL,
#' @field summaryAllResAt [character] If NULL, we will summary the simulation results in `self$trueGVMatList` and `self$estimatedGVMatList`.
#' Else, we will summary the simulation results in `summaryAllResAt` path.
summaryAllResAt = NULL,
#' @field verbose [logical] Display info (optional)
verbose = NULL,
#' @field lociEffectsInit [matrix] Marker and QTL effects used for crossInfo object for initial population
lociEffectsInit = NULL,
#' @field propSelOrNot [numeric] Which ones correspond to proportion of selected individuals
propSelOrNot = NULL,
#' @field hOriLens [numeric] Length of original weighting parameter vector for resource allocation
hOriLens = NULL,
#' @field propSelLens [numeric] Length of proportion of selected individuals vector
propSelLens = NULL,
#' @field hLens [numeric] Length of hyperparameter vector
hLens = NULL,
#' @field hStart [numeric] Initial values of hyper parameters
hStart = NULL,
#' @field solnInit [list] List of solution of StoSOO for initial state
solnInit = NULL,
#' @field hVecOptsList [list] List of optimized hyperparameters in each generation
hVecOptsList = list(),
#' @field optimalHyperParamMatsList [list] List of optimized hyper parameters given finite numbers of budget for optimization
optimalHyperParamMatsList = list(),
#' @field simBsRes [list] Simulation results of each method
simBsRes = list(),
#' @field trueGVMatInit [matrix] A true GV matrix for initial population
trueGVMatInit = NULL,
#' @field trueGVMatList [list] A list of true GV matrix of simulation results
trueGVMatList = list(),
#' @field trueGVSummaryArray [array] An array of summary statistics of true GVs for each population & each iteration
trueGVSummaryArray = NULL,
#' @field estimatedGVMatInit [matrix] An estimated GV matrix for initial population
estimatedGVMatInit = NULL,
#' @field estimatedGVMatList [list] A list of estimated GV matrix of simulation results
estimatedGVMatList = list(),
#' @field estimatedGVSummaryArray [array] An array of summary statistics of estimated GVs for each population & each iteration
estimatedGVSummaryArray = NULL,
#' @description Create a new simBsCma object.
#' @param simBsName [character] Name of this simulation of breeding schemes
#' @param bsInfoInit [bsInfo] breeding scheme info
#' (see:\link[myBreedSimulatR]{bsInfo})
#' @param breederInfoInit [breederInfo] breeder info
#' (see:\link[myBreedSimulatR]{breederInfo})
#' @param lociEffMethod [character] How to compute loci effects (use true QTL effects ("true") or estimated marker effects ("estimated"))
#' @param trainingPopType [character] Training population consists of all populations created in the breeding scheme for `trainingPopType = "all"`,
#' or consists of the latest population in breeding scheme for `trainingPopType = "latest"`.
#' @param trainingPopInit [character / numeric] Training population names or No.s (not generations!!) for initial population
#' @param trainingIndNamesInit [character] Names of training individuals for initia; population
#' @param methodMLRInit [character] Methods for estimating marker effects for initial population.
#' The following methods are offered:
#'
#' "Ridge", "LASSO", "ElasticNet", "RR-BLUP", "GBLUP", "BayesA" (uni-trait),
#' "BayesB" (uni-trait), "BayesC" (uni-trait), "BRR", "BL" (uni-trait), "SpikeSlab" (multi-trait)
#' @param multiTraitInit [logical] Use multiple-trait model for estimation of marker effects or not for initial population
#' @param nIterSimulation [numeric] Number of iterations for this simulation setting
#' @param nGenerationProceed [numeric] Number of generations to be proceeded
#' @param nGenerationProceedSimulation [numeric] Number of generations to be proceeded in simulation for optimizing policy
#' @param setGoalAsFinalGeneration [logical] Set goal for simulation of breeding scheme as the real final generation
#' @param performOptimization [logical] Perform optimization of policy in each generation or not
#' @param useFirstOptimizedValue [logical] Perform optimization of policy only for the initial population and use the optimized hyperprameters permanently.
#' @param performRobustOptimization [logical] Whether or not performing robust optimization.
#' @param lowerQuantile [logical] Lower quantile for the robust optimization (value at risk)
#' @param sameAcrossGeneration [logical] Use same hyper parameter across generations or not
#' @param propSelOpt [logical] Optimize proportion of selected individuals or not
#' @param hMin [numeric] Lower bound of hyper parameters
#' @param hMax [numeric] Upper bound of hyper parameters
#' @param nTotalIterForOneOptimization [numeric] Number of total iterations that can be assigned for one optimization process
#' @param nIterSimulationPerEvaluation [numeric] Number of simulations per one evaluation of the hyperparameter set of your interest
#' @param nIterSimulationForOneMrkEffect [numeric] Number of simulations per one evaluation of the hyperparameter set of your interest for one set of marker effects
#' @param nIterMrkEffectForRobustOptimization [numeric] Number of sets of marker effects utilized for the robust optimization
#' @param nIterOptimization [numeric] Number of iterations required for one optimization
#' @param lambda [numeric] Number of offspring in CMA-ES. Must be greater than or equal to mu.
#' @param mu [numeric] Population size in CMA-ES.
#' @param saveObjNameBase [character] Base name of the Obj to be saved
#' @param whenToSaveObj [numeric] When (how many iterations) to save the Obj in StoSOO
#' @param nTopEvalForOpt [numeric] Number of individuals to be evaluated when evaluating population max or population min for optimization of hyperparameters
#' @param rewardWeightVec [numeric] When returning reward function, `rewardWeightVec` will be multiplied by estimated GVs for each generation to evaluate the method.
#' If you want to apply discounted method, you can achieve by `rewardWeightVec = sapply(1:nGenerationProceedSimulation, function(genProceedNo) gamma ^ (genProceedNo - 1))` where `gamma` is discounted rate.
#' Or if you want to evaluate just the final generation, you can achieve by `rewardWeightVec = c(rep(0, nGenerationProceedSimulation - 1), 1)`.
#' @param digitsEval [numeric] When you evaluate each hyperparameter, you can round the average of evaluates (`eval`) with `round(eval, digitsEval)`.
#' @param nRefreshMemoryEvery [numeric] Every `nRefreshMemoryEvery` iterations, we refresh memory used for simulations by `gc(reset = TRUE)`.
#' @param updateBreederInfo [logical] Update breederInfo or not for each generation
#' @param phenotypingInds [logical] Phenotyping individuals or not for each generation
#' @param nRepForPhenoInit [numeric] Number of replications to be phenotyped for initial population
#' @param nRepForPheno [numeric] Number of replications to be phenotyped for each generation
#' @param updateModels [logical] Whether or not updating the model for each generation
#' (If `phenotypingInds = FALSE` for generation of your interest, `updateModels` will be automatically `FALSE` for that generation.)
#' @param updateBreederInfoSimulation [logical] Update breederInfo or not for each generation in simulations for optimization
#' @param phenotypingIndsSimulation [logical] Phenotyping individuals or not for each generation in simulations for optimization
#' @param nRepForPhenoSimulation [numeric] Number of replications to be phenotyped for each generation in simulations for optimization
#' @param updateModelsSimulation [logical] Whether or not updating the model for each generation in simulations for optimization
#' (If `phenotypingInds = FALSE` for generation of your interest, `updateModels` will be automatically `FALSE` for that generation.)
#' @param performRobustOptimization [logical] Whether or not performing robust optimization.
#' @param methodMLR [character] Methods for estimating marker effects.
#' The following methods are offered:
#'
#' "Ridge", "LASSO", "ElasticNet", "RR-BLUP", "GBLUP", "BayesA" (uni-trait),
#' "BayesB" (uni-trait), "BayesC" (uni-trait), "BRR", "BL" (uni-trait), "SpikeSlab" (multi-trait)
#' @param multiTrait [logical] Use multiple-trait model for estimation of marker effects or not
#' @param nSelectionWaysVec [numeric] Number of selection ways for each generation
#' @param selectionMethodList [list] (list of) Selection method for each generation
#' @param traitNoSelList [list] (list of list of) Number of trait No of your interest for selection for each generation
#' @param blockSplitMethod [character] How to determine the markers belonging to each block when computing OHV.
#' You can define the number of markers in each block (`nMrkInBlock`) or the minimum length of each segment (`minimumSegmentLength`).
#' @param nMrkInBlock [numeric] Number of markers in each block. This will be used for the computation of OHV.
#' @param minimumSegmentLength [numeric] Minimum length of each segment [cM]. This will be used for the computation of OHV.
#' @param nSelInitOPVList [list] (list of) Number of selected candiates for first screening before selecting parent candidates by OPV for each generation
#' @param nIterOPV [numeric] Number of iterations for computation of OPV
#' @param nProgeniesEMBVVec [numeric] Number of progenies of double haploids produced when computing EMBV for each generation
#' @param nIterEMBV [numeric] Number of iterations to estimate EMBV
#' @param nCoresEMBV [numeric] Number of cores used for EMBV estimation
#' @param clusteringForSelList [list] (list of) Apply clustering results of marker genotype to selection or not for each generation
#' @param nClusterList [list] (list of) Number of clusters for each generation
#' @param nTopClusterList [list] (list of) Number of top clusters used for selection for each generation
#' @param nTopEachList [list] (list of) Number of selected individuals in each cluster for each generation
#' @param nSelList [list] (list of) Number of selection candidates for each generation
#' @param multiTraitsEvalMethodList [list] (list of) When evaluating multiple traits, you can choose how to evaluate these traits simultaneously.
#' One method is to take a weighted sum of these traits, and the other is to compute a product of these traits (adjusted to be positive values in advance.)
#' @param hSelList [list] (list of) Hyperparameter which determines which trait is weighted when selecting parent candidates for multiple traits for each generation
#' @param matingMethodVec [character] Mating method for each generation
#' @param allocateMethodVec [character] Allocation method for each generation
#' @param weightedAllocationMethodList [list] (list of) Which selection index will be used for weighted resource allocation for each generation
#' @param traitNoRAList [list] (list of) Trait No of your interest for resource allocation for each generation
#' @param minimumUnitAllocateVec [numeric] (vector of) Minimum number of allocated progenies for each pair.
#' @param includeGVPVec [logical] Whether or not to consider genetic variance of progenies of each pair when determining the number of progenies per each pair for each generation
#' @param targetGenGVPVec [numeric] Vector of target generation of selfing when evaluating the genetic variance of progenies (F#)
#' @param nNextPopVec [numeric] Vector of number of progenies for the next generation for each generation
#' @param performOCSVec [logical] Vector of whether or not performing OCS before allocation
#' @param targetGenOCSVec [numeric] Vector of target generation in OCS
#' @param HeStarRatioVec [numeric] Vector of ratio of He0 to HeStar
#' @param degreeOCSVec [numeric] Vector of how generation effect is emphasized when computing He(t)
#' @param includeGVPOCSVec [logical] Vector of whether or not to consider genetic variance of progenies of each pair when evaluating each pair in OCS
#' @param hOCSList [list] List of weight for each trait/criterion to evaluate each mating pair
#' @param weightedAllocationMethodOCSList [list] List of which selection index will be used to evaluate mating pair in OCS
#' @param traitNoOCSList [list] List of trait No to evaluate mating pair in OCS
#' @param nCrossesOCSVec [numeric] Vector of number of crosses that will be selected in OCS
#' @param nameMethod [character] Method for naming individuals
#' @param nCores [numeric] Number of cores used for simulations of breeding scheme
#' @param nCoresPerOptimization [numeric] Number of cores used for simulations of breeding scheme when optimizing hyperparameters
#' @param overWriteRes [logical] Overwrite simulation results when the targeted results already exists
#' @param showProgress [logical] Show progress bar or not
#' @param returnMethod [character] Which type of results will be returned (saved) in the object
#' @param saveAllResAt [character] If NULL, we won't save the simulation results. Else, we will save bsInfo and breederInfo for each iteration in the directory defined by `saveAllResAt` path.
#' @param evaluateGVMethod [character] Which type of GV (true GV or estimated GV) will be used for evaluation of the simulation results
#' @param nTopEval [numeric] Number of individuals to be evaluated when evaluating population max or population min#'
#' @param traitNoEval [numeric] Trait No of your interest for evaluation of the method
#' @param hEval [numeric] Hyperparameter which determines which BV is emphasized when evaluating simulation results of each method
#' @param summaryAllResAt [character] If NULL, we will summary the simulation results in `self$trueGVMatList` and `self$estimatedGVMatList`.
#' Else, we will summary the simulation results in `summaryAllResAt` path.
#' @param verbose [logical] Display info (optional)
#' @return A new `simBs` object.
#' @examples
#' ### create simulation information
#' mySimInfo <- simInfo$new(simName = "Simulation Example",
#' simGeno = TRUE,
#' simPheno = TRUE,
#' #' nSimGeno = 1,
#' #' nSimPheno = 3,
#' #' nCoreMax = 4,
#' #' nCorePerGeno = 1,
#' #' nCorePerPheno = 3,
#' saveDataFileName = NULL)
#
#' ### create specie information
#' mySpec <- specie$new(nChr = 3,
#' lChr = c(100, 150, 200),
#' specName = "Example 1",
#' ploidy = 2,
#' mutRate = 10^-8,
#' recombRate = 10^-6,
#' chrNames = c("C1", "C2", "C3"),
#' nLoci = 100,
#' recombRateOneVal = FALSE,
#' effPopSize = 100,
#' simInfo = mySimInfo,
#' verbose = TRUE)
#
#' ### create lociInfo object
#' myLoci <- lociInfo$new(genoMap = NULL, specie = mySpec)
#' plot(myLoci, alpha = 0.1)
#
#' ### create traitInfo object
#' myTrait <- traitInfo$new(lociInfo = myLoci,
#' nMarkers = 80,
#' nTraits = 3,
#' nQTLs = c(4, 8, 3),
#' actionTypeEpiSimple = TRUE,
#' qtlOverlap = TRUE,
#' nOverlap = c(2, 3, 0),
#' effCor = 0.1,
#' propDomi = 0.2,
#' interactionMean = c(4, 1, 2))
#' myTrait$plot(alphaMarker = 0.1)
#
#
#
#' ### create bsInfo object
#' myBsInfo <- bsInfo$new(simInfo = mySimInfo,
#' specie = mySpec,
#' lociInfo = myLoci,
#' traitInfo = myTrait,
#' geno = NULL,
#' haplo = NULL,
#' founderIsInitPop = TRUE,
#' popNameBase = "Population",
#' initIndNames = NULL,
#' verbose = TRUE)
#
#' ### create cross information object
#' for (i in 1:10) {
#' myCrossInfo <- crossInfo$new(parentPopulation = myBsInfo$populations[[myBsInfo$generation]],
#' method = "randomMate",
#' nNextPop = 100)
#
#' myBsInfo$nextGeneration(crossInfo = myCrossInfo)
#' }
#' geno <- myBsInfo$overGeneration()$genoMat
#' myBsInfo$print()
#' myBsInfo$plot(plotTarget = "trueAGV",
#' targetTrait = 1:3,
#' targetPopulation = 1:11,
#' plotType = "jitter")
#'
initialize = function(simBsName = "Undefined",
bsInfoInit,
breederInfoInit = NULL,
lociEffMethod = NULL,
trainingPopType = NULL,
trainingPopInit = NULL,
trainingIndNamesInit = NULL,
methodMLRInit = NULL,
multiTraitInit = FALSE,
nIterSimulation = NULL,
nGenerationProceed = NULL,
nGenerationProceedSimulation = NULL,
setGoalAsFinalGeneration = TRUE,
performOptimization = NULL,
useFirstOptimizedValue = NULL,
performRobustOptimization = FALSE,
lowerQuantile = NULL,
sameAcrossGeneration = TRUE,
propSelOpt = TRUE,
hMin = NULL,
hMax = NULL,
nTotalIterForOneOptimization = NULL,
nIterSimulationPerEvaluation = NULL,
nIterSimulationForOneMrkEffect = NULL,
nIterMrkEffectForRobustOptimization = NULL,
nIterOptimization = NULL,
lambda = NULL,
mu = NULL,
saveObjNameBase = NULL,
whenToSaveObj = NA,
nTopEvalForOpt = NULL,
rewardWeightVec = NULL,
digitsEval = NULL,
nRefreshMemoryEvery = NULL,
updateBreederInfo = TRUE,
phenotypingInds = FALSE,
nRepForPhenoInit = NULL,
nRepForPheno = NULL,
updateModels = FALSE,
updateBreederInfoSimulation = TRUE,
phenotypingIndsSimulation = FALSE,
nRepForPhenoSimulation = NULL,
updateModelsSimulation = FALSE,
methodMLR = NULL,
multiTrait = FALSE,
nSelectionWaysVec = NULL,
selectionMethodList = NULL,
traitNoSelList = NULL,
blockSplitMethod = NULL,
nMrkInBlock = NULL,
minimumSegmentLength = NULL,
nSelInitOPVList = NULL,
nIterOPV = NULL,
nProgeniesEMBVVec = NULL,
nIterEMBV = NULL,
nCoresEMBV = NULL,
clusteringForSelList = FALSE,
nClusterList = NULL,
nTopClusterList = NULL,
nTopEachList = NULL,
nSelList = NULL,
multiTraitsEvalMethodList = NULL,
hSelList = NULL,
matingMethodVec = NULL,
allocateMethodVec = NULL,
weightedAllocationMethodList = NULL,
traitNoRAList = NULL,
minimumUnitAllocateVec = NULL,
includeGVPVec = FALSE,
targetGenGVPVec = NULL,
nNextPopVec = NULL,
performOCSVec = NULL,
targetGenOCSVec = NULL,
HeStarRatioVec = NULL,
degreeOCSVec = NULL,
includeGVPOCSVec = FALSE,
hOCSList = NULL,
weightedAllocationMethodOCSList = NULL,
traitNoOCSList = NULL,
nCrossesOCSVec = NULL,
nameMethod = "pairBase",
nCores = NULL,
nCoresPerOptimization = NULL,
overWriteRes = FALSE,
showProgress = TRUE,
returnMethod = "all",
saveAllResAt = NULL,
evaluateGVMethod = "true",
nTopEval = NULL,
traitNoEval = NULL,
hEval = NULL,
summaryAllResAt = NULL,
verbose = TRUE) {
# define some methods
lociEffMethodsOffered <- c("true", "estimated")
trainingPopTypesOffered <- c("all", "latest")
supportedMethodsMLR <- c("Ridge", "LASSO", "ElasticNet", "RR-BLUP", "GBLUP",
"BayesA", "BayesB", "BayesC", "BRR", "BL", "SpikeSlab")
supportedMethodsGlmnet <- c("Ridge", "LASSO", "ElasticNet")
supportedMethodsBGLR <- c("BayesA", "BayesB", "BayesC", "BRR", "BL", "SpikeSlab")
selectionMethodsOffered <- c("nonSelection", "selectBV", "selectWBV", "selectOHV",
"selectEMBV", "selectOPV", "userSI", "userSpecific")
selectionMethodsWithSelection <- c("selectBV", "selectWBV", "selectOHV",
"selectEMBV", "selectOPV", "userSI", "userSpecific")
selectionMethodsWithMrkEff <- c("selectBV", "selectWBV", "selectOHV", "selectEMBV", "selectOPV")
matingMethodsOffered <- c("randomMate", "roundRobin", "diallel", "diallelWithSelfing",
"selfing", "maxGenDist", "makeDH", "nonCross", "userSpecific")
allocateMethodsOffered <- c("equalAllocation", "weightedAllocation", "userSpecific")
blockSplitMethodsOffered <- c("nMrkInBlock", "minimumSegmentLength")
multiTraitsEvalMethodsOffered <- c("sum", "prod")
nameMethodsOffered <- c("pairBase", "individualBase")
returnMethodsOffered <- c("all", "summary", "max", "mean", "median", "min", "var")
# simBsName
if (is.null(simBsName)) {
simBsName <- "Undefined"
}
stopifnot(is.character(simBsName))
# bsInfoInit class
if (class(bsInfoInit)[1] != "bsInfo") {
stop(paste('class(bsInfoInit)[1] != "bsInfo"\n"bsInfo" must be a',
'bsInfo object see: ?bsInfo'))
}
# define some variables
nIndNow <- bsInfoInit$populations[[length(bsInfoInit$populations)]]$nInd
nTraits <- bsInfoInit$traitInfo$nTraits
# lociEffMethod
if (!is.null(lociEffMethod)) {
if (!(lociEffMethod %in% lociEffMethodsOffered)) {
stop(paste0("We only offer the following methods for naming individuals: ",
paste(lociEffMethodsOffered, collapse = "; ")))
}
} else {
lociEffMethod <- "estimated"
message(paste0("`lociEffMethod` is not specified. We substitute `lociEffMethod = ",
lociEffMethod,"` instead."))
}
# nIterSimulation
if (!is.null(nIterSimulation)) {
stopifnot(is.numeric(nIterSimulation))
nIterSimulation <- floor(nIterSimulation)
stopifnot(nIterSimulation >= 1)
} else {
nIterSimulation <- 1
message(paste0("`nIterSimulation` is not specified. We substitute `nIterSimulation = ",
nIterSimulation,"` instead."))
}
# nGenerationProceed
if (!is.null(nGenerationProceed)) {
stopifnot(is.numeric(nGenerationProceed))
nGenerationProceed <- floor(nGenerationProceed)
stopifnot(nGenerationProceed >= 1)
} else {
nGenerationProceed <- 1
message(paste0("`nGenerationProceed` is not specified. We substitute `nGenerationProceed = ",
nGenerationProceed,"` instead."))
}
# nGenerationProceedSimulation
if (!is.null(nGenerationProceedSimulation)) {
stopifnot(is.numeric(nGenerationProceedSimulation))
nGenerationProceedSimulation <- floor(nGenerationProceedSimulation)
stopifnot(nGenerationProceedSimulation >= 1)
} else {
nGenerationProceedSimulation <- nGenerationProceed
message(paste0("`nGenerationProceedSimulation` is not specified. We substitute `nGenerationProceedSimulation = ",
nGenerationProceedSimulation,"` instead."))
}
# setGoalAsFinalGeneration
stopifnot(is.logical(setGoalAsFinalGeneration))
# performOptimization
if (!is.null(performOptimization)) {
stopifnot(is.logical(performOptimization))
} else {
performOptimization <- rep(FALSE, nGenerationProceed)
performOptimization[1] <- TRUE
message(paste0("`performOptimization` is not specified. Instead, we substitute `performOptimization = c(",
paste(performOptimization, collapse = ", "), ")`"))
}
if (!(length(performOptimization) %in% c(1, nGenerationProceed))) {
stop(paste("length(performOptimization) must be equal to 1 or equal to nGenerationProceed."))
} else if (length(performOptimization) == 1) {
performOptimization <- rep(performOptimization, nGenerationProceed)
}
names(performOptimization) <- 1:nGenerationProceed
# useFirstOptimizedValue
if (!is.null(useFirstOptimizedValue)){
stopifnot(is.logical(useFirstOptimizedValue))
if (useFirstOptimizedValue) {
if (nGenerationProceed != nGenerationProceedSimulation) {
useFirstOptimizedValue <- FALSE
}
if (any(performOptimization[-1])) {
useFirstOptimizedValue <- FALSE
}
if (!useFirstOptimizedValue) {
message(paste0("You cannot set `useFirstOptimizedValue` as `TRUE` with this `nGenerationProceed`, `nGenerationProceedSimulation`, and `performOptimization`. We substitute `useFirstOptimizedValue = ",
useFirstOptimizedValue,"` instead."))
}
}
} else {
if ((nGenerationProceed == nGenerationProceedSimulation) | (!all(performOptimization[-1]))) {
useFirstOptimizedValue <- TRUE
} else {
useFirstOptimizedValue <- FALSE
}
message(paste0("`useFirstOptimizedValue` is not specified. We substitute `useFirstOptimizedValue = ",
useFirstOptimizedValue,"` instead."))
}
# performRobustOptimization
if (performRobustOptimization & (lociEffMethod == "true")) {
performRobustOptimization <- FALSE
message(paste0("When `lociEffMethod == 'true'`, we don't need to perform robust optimization. We substitute `performRobustOptimization = ",
performRobustOptimization,"` instead."))
}
# lowerQuantile
if (!is.null(lowerQuantile)) {
stopifnot(is.numeric(lowerQuantile))
stopifnot(lowerQuantile > 0)
stopifnot(lowerQuantile < 1)
} else {
lowerQuantile <- 0.05
message(paste0("`lowerQuantile` is not specified. We substitute `lowerQuantile = ",
lowerQuantile,"` instead."))
}
# nTotalIterForOneOptimization
if (!is.null(nTotalIterForOneOptimization)) {
stopifnot(is.numeric(nTotalIterForOneOptimization))
nTotalIterForOneOptimization <- floor(nTotalIterForOneOptimization)
stopifnot(nTotalIterForOneOptimization >= 1)
} else {
nTotalIterForOneOptimization <- 1e04
message(paste0("`nTotalIterForOneOptimization` is not specified. We substitute `nTotalIterForOneOptimization = ",
nTotalIterForOneOptimization,"` instead."))
}
# nIterSimulationPerEvaluation
if (!is.null(nIterSimulationPerEvaluation)) {
stopifnot(is.numeric(nIterSimulationPerEvaluation))
nIterSimulationPerEvaluation <- floor(nIterSimulationPerEvaluation)
stopifnot(nIterSimulationPerEvaluation >= 1)
} else {
nIterSimulationPerEvaluation <- 10
message(paste0("`nIterSimulationPerEvaluation` is not specified. We substitute `nIterSimulationPerEvaluation = ",
nIterSimulationPerEvaluation,"` instead."))
}
# nIterSimulationForOneMrkEffect
if (!is.null(nIterSimulationForOneMrkEffect)) {
stopifnot(is.numeric(nIterSimulationForOneMrkEffect))
nIterSimulationForOneMrkEffect <- floor(nIterSimulationForOneMrkEffect)
stopifnot(nIterSimulationForOneMrkEffect >= 1)
} else {
if (performRobustOptimization) {
nIterSimulationForOneMrkEffect <- 1
} else {
nIterSimulationForOneMrkEffect <- nIterSimulationPerEvaluation
}
message(paste0("`nIterSimulationForOneMrkEffect` is not specified. We substitute `nIterSimulationForOneMrkEffect = ",
nIterSimulationForOneMrkEffect,"` instead."))
}
# nIterOptimization
if (!is.null(nIterMrkEffectForRobustOptimization)) {
stopifnot(is.numeric(nIterMrkEffectForRobustOptimization))
nIterMrkEffectForRobustOptimization <- floor(nIterMrkEffectForRobustOptimization)
stopifnot(nIterMrkEffectForRobustOptimization >= 1)
} else {
if (performRobustOptimization) {
nIterMrkEffectForRobustOptimization <- round(nIterSimulationPerEvaluation / nIterSimulationForOneMrkEffect)
} else {
nIterMrkEffectForRobustOptimization <- 1
}
message(paste0("`nIterMrkEffectForRobustOptimization` is not specified. We substitute `nIterMrkEffectForRobustOptimization = ",
nIterMrkEffectForRobustOptimization,"` instead."))
}
if (nIterSimulationPerEvaluation != nIterSimulationForOneMrkEffect * nIterMrkEffectForRobustOptimization) {
nIterSimulationPerEvaluation <- nIterSimulationForOneMrkEffect * nIterMrkEffectForRobustOptimization
message((paste0("We substitute `nIterSimulationPerEvaluation = ",
nIterSimulationPerEvaluation,"` instead.")))
}
# lambda
if (!is.null(lambda)) {
stopifnot(is.numeric(lambda))
} else {
lambda <- 4 + floor(3 * log(paramLen))
message(paste0("You do not specify `lambda`. We set `lambda = ",
lambda, "`."))
}
# mu
if (!is.null(mu)) {
stopifnot(is.numeric(mu))
} else {
mu <- floor(lambda / 2)
message(paste0("You do not specify `mu`. We set `mu = ",
mu, "`."))
}
# nIterOptimization
if (!is.null(nIterOptimization)) {
stopifnot(is.numeric(nIterOptimization))
nIterOptimization <- floor(nIterOptimization)
stopifnot(nIterOptimization >= 1)
} else {
nIterOptimization <- round(nTotalIterForOneOptimization / (nIterSimulationPerEvaluation * lambda))
message(paste0("`nIterOptimization` is not specified. We substitute `nIterOptimization = ",
nIterOptimization,"` instead."))
}
if (nTotalIterForOneOptimization != nIterSimulationPerEvaluation * nIterOptimization * lambda) {
nTotalIterForOneOptimization <- nIterSimulationPerEvaluation * nIterOptimization * lambda
message((paste0("We substitute `nTotalIterForOneOptimization = ",
nTotalIterForOneOptimization,"` instead.")))
}
# saveObjNameBase
if (is.null(saveObjNameBase)) {
whenToSaveObj <- NULL
}
# whenToSaveObj
if (!is.null(whenToSaveObj)) {
if (!all(is.na(whenToSaveObj))) {
stopifnot(is.numeric(whenToSaveObj))
whenToSaveObj <- whenToSaveObj[!is.na(whenToSaveObj)]
whenToSaveObj <- floor(whenToSaveObj)
stopifnot(all(whenToSaveObj >= 1))
stopifnot(all(whenToSaveObj <= nIterOptimization))
} else {
whenToSaveObj <- nIterOptimization
message(paste0("You do not specify `whenToSaveObj`. We set `whenToSaveObj = ",
whenToSaveObj, "`."))
}
}
# rewardWeightVec
if (!is.null(rewardWeightVec)) {
stopifnot(is.numeric(rewardWeightVec))
stopifnot(all(rewardWeightVec >= 0))
} else {
rewardWeightVec <- c(rep(0, nGenerationProceedSimulation - 1), 1)
message(paste0("`rewardWeightVec` is not specified. Instead, we substitute `rewardWeightVec = c(",
paste(rewardWeightVec, collapse = ", "), ")`."))
}
if (!(length(rewardWeightVec) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(rewardWeightVec) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(rewardWeightVec) == 1) {
rewardWeightVec <- rep(rewardWeightVec, nGenerationProceedSimulation)
}
names(rewardWeightVec) <- 1:nGenerationProceedSimulation
# digitsEval
if (!is.null(digitsEval)) {
stopifnot(is.numeric(digitsEval))
digitsEval <- floor(digitsEval)
} else {
digitsEval <- 3
message(paste0("`digitsEval` is not specified. We substitute `digitsEval = ",
digitsEval,"` instead."))
}
# nRefreshMemoryEvery
if (!is.null(nRefreshMemoryEvery)) {
stopifnot(is.numeric(nRefreshMemoryEvery))
nRefreshMemoryEvery <- floor(nRefreshMemoryEvery)
} else {
nRefreshMemoryEvery <- 2
message(paste0("`nRefreshMemoryEvery` is not specified. We substitute `nRefreshMemoryEvery = ",
nRefreshMemoryEvery,"` instead."))
}
# updateBreederInfo
stopifnot(is.logical(updateBreederInfo))
if (!(length(updateBreederInfo) %in% c(1, nGenerationProceed))) {
stop(paste("length(updateBreederInfo) must be equal to 1 or equal to nGenerationProceed."))
} else if (length(updateBreederInfo) == 1) {
updateBreederInfo <- rep(updateBreederInfo, nGenerationProceed)
}
names(updateBreederInfo) <- 1:nGenerationProceed
# phenotypingInds
stopifnot(is.logical(phenotypingInds))
if (!(length(phenotypingInds) %in% c(1, nGenerationProceed))) {
stop(paste("length(phenotypingInds) must be equal to 1 or equal to nGenerationProceed."))
} else if (length(phenotypingInds) == 1) {
phenotypingInds <- rep(phenotypingInds, nGenerationProceed)
}
names(phenotypingInds) <- 1:nGenerationProceed
phenotypingInds[!updateBreederInfo] <- FALSE
# nRepForPhenoInit
if (!is.null(nRepForPhenoInit)) {
stopifnot(is.numeric(nRepForPhenoInit))
nRepForPhenoInit <- floor(nRepForPhenoInit)
stopifnot(nRepForPhenoInit >= 1)
} else {
nRepForPhenoInit <- 1
message(paste0("`nRepForPhenoInit` is not specified. We substitute `nRepForPhenoInit = ",
nRepForPhenoInit,"` instead."))
}
# nRepForPheno
if (!is.null(nRepForPheno)) {
stopifnot(is.numeric(nRepForPheno))
nRepForPheno <- floor(nRepForPheno)
stopifnot(all(nRepForPheno >= 0))
} else {
nRepForPheno <- 1
message(paste0("`nRepForPheno` is not specified. We substitute `nRepForPheno = ",
nRepForPheno,"` instead."))
}
if (!(length(nRepForPheno) %in% c(1, nGenerationProceed))) {
stop(paste("length(nRepForPheno) must be equal to 1 or equal to nGenerationProceed."))
} else if (length(nRepForPheno) == 1) {
nRepForPheno <- rep(nRepForPheno, nGenerationProceed)
}
names(nRepForPheno) <- 1:nGenerationProceed
nRepForPheno[!phenotypingInds] <- 0
# updateBreederInfoSimulation
stopifnot(is.logical(updateBreederInfoSimulation))
if (!(length(updateBreederInfoSimulation) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(updateBreederInfoSimulation) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(updateBreederInfoSimulation) == 1) {
updateBreederInfoSimulation <- rep(updateBreederInfoSimulation, nGenerationProceedSimulation)
}
names(updateBreederInfoSimulation) <- 1:nGenerationProceedSimulation
# phenotypingIndsSimulation
stopifnot(is.logical(phenotypingIndsSimulation))
if (!(length(phenotypingIndsSimulation) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(phenotypingIndsSimulation) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(phenotypingIndsSimulation) == 1) {
phenotypingIndsSimulation <- rep(phenotypingIndsSimulation, nGenerationProceedSimulation)
}
names(phenotypingIndsSimulation) <- 1:nGenerationProceedSimulation
phenotypingIndsSimulation[!updateBreederInfoSimulation] <- FALSE
# nRepForPhenoSimulation
if (!is.null(nRepForPhenoSimulation)) {
stopifnot(is.numeric(nRepForPhenoSimulation))
nRepForPhenoSimulation <- floor(nRepForPhenoSimulation)
stopifnot(all(nRepForPhenoSimulation >= 0))
} else {
nRepForPhenoSimulation <- 1
message(paste0("`nRepForPhenoSimulation` is not specified. We substitute `nRepForPhenoSimulation = ",
nRepForPhenoSimulation,"` instead."))
}
if (!(length(nRepForPhenoSimulation) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(nRepForPhenoSimulation) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(nRepForPhenoSimulation) == 1) {
nRepForPhenoSimulation <- rep(nRepForPhenoSimulation, nGenerationProceedSimulation)
}
names(nRepForPhenoSimulation) <- 1:nGenerationProceedSimulation
nRepForPhenoSimulation[!phenotypingIndsSimulation] <- 0
# trainingPopType
if (!is.null(trainingPopType)) {
if (!(trainingPopType %in% trainingPopTypesOffered)) {
stop(paste0("We only offer the following methods for naming individuals: ",
paste(trainingPopTypesOffered, collapse = "; ")))
}
} else {
trainingPopType <- "latest"
message(paste0("`trainingPopType` is not specified. We substitute `trainingPopType = ",
trainingPopType,"` instead."))
}
# multiTraitInit
stopifnot(is.logical(multiTraitInit))
# methodMLRInit
if (!is.null(methodMLRInit)) {
if (!(methodMLRInit %in% supportedMethodsMLR)) {
stop(paste0("We only offer the following methods for naming individuals: ",
paste(supportedMethodsMLR, collapse = "; ")))
}
} else {
if (multiTraitInit) {
methodMLRInit <- "SpikeSlab"
} else {
methodMLRInit <- "BayesB"
}
message(paste0("`methodMLRInit` is not specified. We substitute `methodMLRInit = ",
methodMLRInit,"` instead."))
}
if (methodMLRInit %in% supportedMethodsBGLR) {
if (!multiTraitInit) {
if (methodMLRInit == "SpikeSlab") {
message(paste0("For uni-trait model, `methodMLR = 'SpikeSlab'` is not offered.\n",
"We use `methodMLR = 'BayesC'` instead."))
methodMLRInit <- "BayesC"
}
} else {
if (methodMLRInit %in% c("BL", "BayesB", "BayesC")) {
message(paste0("For multi-trait model, `methodMLR = 'BL'`, `methodMLR = 'BayesB'`, `methodMLR = 'BayesC'` are not offered.\n",
"We use `methodMLR = 'SpikeSlab'` instead. This is equivalent to BayesC model."))
methodMLRInit <- "SpikeSlab"
} else if (methodMLRInit == "BayesA") {
message(paste0("For multi-trait model, `methodMLR = 'BayesA'` is not offered.\n",
"We use `methodMLR = 'BRR'` instead."))
methodMLRInit <- "BRR"
}
}
}
# updateModels
stopifnot(is.logical(updateModels))
if (!(length(updateModels) %in% c(1, nGenerationProceed))) {
stop(paste("length(updateModels) must be equal to 1 or equal to nGenerationProceed."))
} else if (length(updateModels) == 1) {
updateModels <- rep(updateModels, nGenerationProceed)
}
names(updateModels) <- 1:nGenerationProceed
updateModels[!phenotypingInds] <- FALSE
# updateModelsSimulation
stopifnot(is.logical(updateModelsSimulation))
if (!(length(updateModelsSimulation) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(updateModelsSimulation) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(updateModelsSimulation) == 1) {
updateModelsSimulation <- rep(updateModelsSimulation, nGenerationProceedSimulation)
}
names(updateModelsSimulation) <- 1:nGenerationProceedSimulation
updateModelsSimulation[!phenotypingIndsSimulation] <- FALSE
# multiTrait
stopifnot(is.logical(multiTrait))
# methodMLR
if (!is.null(methodMLR)) {
if (!(methodMLR %in% supportedMethodsMLR)) {
stop(paste0("We only offer the following methods for naming individuals: ",
paste(supportedMethodsMLR, collapse = "; ")))
}
} else {
methodMLR <- "LASSO"
message(paste0("`methodMLR` is not specified. We substitute `methodMLR = ",
methodMLR,"` instead."))
}
if (methodMLR %in% supportedMethodsBGLR) {
if (!multiTrait) {
if (methodMLR == "SpikeSlab") {
message(paste0("For uni-trait model, `methodMLR = 'SpikeSlab'` is not offered.\n",
"We use `methodMLR = 'BayesC'` instead."))
methodMLR <- "BayesC"
}
} else {
if (methodMLR %in% c("BL", "BayesB", "BayesC")) {
message(paste0("For multi-trait model, `methodMLR = 'BL'`, `methodMLR = 'BayesB'`, `methodMLR = 'BayesC'` are not offered.\n",
"We use `methodMLR = 'SpikeSlab'` instead. This is equivalent to BayesC model."))
methodMLR <- "SpikeSlab"
} else if (methodMLR == "BayesA") {
message(paste0("For multi-trait model, `methodMLR = 'BayesA'` is not offered.\n",
"We use `methodMLR = 'BRR'` instead."))
methodMLR <- "BRR"
}
}
}
# breederInfoInit class
if (!is.null(breederInfoInit)) {
if (class(breederInfoInit)[1] != "breederInfo") {
stop(paste('class(breederInfoInit)[1] != "breederInfo"\n"breederInfo" must be a',
'breederInfo object see: ?breederInfo'))
}
} else {
breederInfoInit <- myBreedSimulatR::breederInfo$new(breederName = "Undefined",
bsInfo = bsInfoInit,
mrkNames = NULL,
initGenotyping = TRUE,
initGenotypedIndNames = NULL,
mafThres = 0.05,
heteroThres = 1,
calculateGRMBase = TRUE,
methodsGRMBase = "addNOIA",
calcEpistasisBase = FALSE,
verbose = verbose)
breederInfoInit$phenotyper(bsInfo = bsInfoInit,
generationOfInterest = NULL,
nRep = nRepForPhenoInit,
multiTraitsAsEnvs = FALSE,
phenotypedIndNames = NULL,
estimateGV = TRUE,
estimatedGVMethod = "lme4")
}
# nSelectionWaysVec
if (!is.null(nSelectionWaysVec)) {
stopifnot(is.numeric(nSelectionWaysVec))
nSelectionWaysVec <- floor(nSelectionWaysVec)
stopifnot(all(nSelectionWaysVec >= 1))
} else {
nSelectionWaysVec <- 1
message(paste0("`nSelectionWaysVec` is not specified. We substitute `nSelectionWaysVec = ",
nSelectionWaysVec,"` instead."))
}
if (!(length(nSelectionWaysVec) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(nSelectionWaysVec) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(nSelectionWaysVec) == 1) {
nSelectionWaysVec <- rep(nSelectionWaysVec, nGenerationProceedSimulation)
}
names(nSelectionWaysVec) <- 1:nGenerationProceedSimulation
# selectionMethodList
if (!is.null(selectionMethodList)) {
if (!is.list(selectionMethodList)) {
selectionMethodList <- list(selectionMethodList)
}
} else {
selectionMethodList <- "nonSelection"
message(paste0("`selectionMethodList` is not specified. We substitute `selectionMethodList = list(",
selectionMethodList,")` instead."))
selectionMethodList <- sapply(X = nSelectionWaysVec,
FUN = function(nSelectionWays) {
rep(selectionMethodList, nSelectionWays)
}, simplify = FALSE)
}
if (!(length(selectionMethodList) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(selectionMethodList) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(selectionMethodList) == 1) {
selectionMethodList <- rep(selectionMethodList, nGenerationProceedSimulation)
}
names(selectionMethodList) <- 1:nGenerationProceedSimulation
stopifnot(all(unlist(lapply(selectionMethodList, length)) == nSelectionWaysVec))
stopifnot(all(unlist(lapply(selectionMethodList, function(x) all(x %in% selectionMethodsOffered)))))
whereSelectionList <- lapply(selectionMethodList, function(x) x %in% selectionMethodsWithSelection)
# traitNoSelList
if (!is.null(traitNoSelList)) {
if (!is.list(traitNoSelList)) {
traitNoSelList <- sapply(nSelectionWaysVec,
function(nSelectionWays) {
traitNoSelListNow <- rep(list(traitNoSelList), nSelectionWays)
return(traitNoSelListNow)
}, simplify = FALSE)
} else if (!is.list(traitNoSelList[[1]])) {
traitNoSelList <- sapply(nSelectionWaysVec,
function(nSelectionWays) {
traitNoSelListNow <- rep(traitNoSelList, nSelectionWays)
return(traitNoSelListNow)
}, simplify = FALSE)
}
} else {
traitNoSelList <- 1
message(paste0("`traitNoSelList` is not specified. We substitute `traitNoSelList = list(list(",
traitNoSelList,"))` instead."))
traitNoSelList <- sapply(nSelectionWaysVec,
function(nSelectionWays) {
traitNoSelListNow <- rep(list(traitNoSelList), nSelectionWays)
return(traitNoSelListNow)
}, simplify = FALSE)
}
if (!(length(traitNoSelList) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(traitNoSelList) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(traitNoSelList) == 1) {
traitNoSelList <- rep(traitNoSelList, nGenerationProceedSimulation)
}
stopifnot(all(unlist(lapply(traitNoSelList, is.list))))
stopifnot(all(unlist(lapply(traitNoSelList, function(traitNoSel) all(unlist(lapply(traitNoSel, is.numeric)))))))
stopifnot(all(sapply(traitNoSelList, function(traitNoSel) all(unlist(lapply(traitNoSel, function(x) all(x >= 1)))))))
stopifnot(all(sapply(traitNoSelList, function(traitNoSel) all(unlist(lapply(traitNoSel, function(x) all(x <= nTraits)))))))
stopifnot(all(unlist(lapply(traitNoRAList, length)) == nSelectionWaysVec))
names(traitNoSelList) <- 1:nGenerationProceedSimulation
# blockSplitMethod
if (!is.null(blockSplitMethod)) {
if (!(blockSplitMethod %in% blockSplitMethodsOffered)) {
stop(paste0("We only offer the following block splitting methods: ",
paste(blockSplitMethodsOffered, collapse = "; ")))
}
} else {
blockSplitMethod <- "minimumSegmentLength"
message("You do not specify the block splitting method. We will define blocks by the minimum length of segements.")
}
# nMrkInBlock
if (!is.null(nMrkInBlock)) {
stopifnot(is.numeric(nMrkInBlock))
nMrkInBlock <- floor(nMrkInBlock)
stopifnot(nMrkInBlock >= 1)
stopifnot(nMrkInBlock <= min(bsInfoInit$specie$nLoci))
} else {
nMrkInBlock <- min(bsInfoInit$specie$nLoci) %/% 10
if (any(unlist(lapply(selectionMethodList, function(selectionMethod) c("selectOHV", "selectOPV") %in% selectionMethod)))) {
message(paste0("`nMrkInBlock` is not specified even though you choose 'selectOHV' / 'selectOPV' method. We substitute `nMrkInBlock = ", nMrkInBlock,"` instead."))
}
}
# minimumSegmentLength
if (!is.null(minimumSegmentLength)) {
stopifnot(is.numeric(minimumSegmentLength))
minimumSegmentLength <- floor(minimumSegmentLength)
stopifnot(minimumSegmentLength >= 1)
stopifnot(minimumSegmentLength <= min(bsInfoInit$specie$lChr))
} else {
minimumSegmentLength <- 6.25
if (any(unlist(lapply(selectionMethodList, function(selectionMethod) c("selectOHV", "selectOPV") %in% selectionMethod)))) {
message(paste0("`nMrkInBlock` is not specified even though you choose 'selectOHV' / 'selectOPV' method. We substitute `nMrkInBlock = ", nMrkInBlock,"` instead."))
}
}
# nIterOPV
if (!is.null(nIterOPV)) {
stopifnot(is.numeric(nIterOPV))
nIterOPV <- floor(nIterOPV)
stopifnot(nIterOPV >= 1)
} else {
nIterOPV <- 5000
if (any(unlist(lapply(selectionMethodList, function(selectionMethod) "selectOPV" %in% selectionMethod)))) {
message(paste0("`nIterOPV` is not specified even though you choose 'selectOPV' method. We substitute `nIterOPV = ", nIterOPV,"` instead."))
}
}
# clusteringForSelList
if (!is.null(clusteringForSelList)) {
if (!is.list(clusteringForSelList)) {
clusteringForSelList <- list(clusteringForSelList)
}
} else {
clusteringForSelList <- FALSE
message(paste0("`clusteringForSelList` is not specified. We substitute `clusteringForSelList = list(",
clusteringForSelList,")` instead."))
clusteringForSelList <- sapply(X = nSelectionWaysVec,
FUN = function(nSelectionWays) {
rep(clusteringForSelList, nSelectionWays)
}, simplify = FALSE)
}
if (!(length(clusteringForSelList) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(clusteringForSelList) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(clusteringForSelList) == 1) {
clusteringForSelList <- rep(clusteringForSelList, nGenerationProceedSimulation)
}
names(clusteringForSelList) <- 1:nGenerationProceedSimulation
stopifnot(all(unlist(lapply(clusteringForSelList, length)) == nSelectionWaysVec))
stopifnot(all(unlist(lapply(clusteringForSelList, is.logical))))
# nSelList
if (!is.null(nSelList)) {
if (!is.list(nSelList)) {
nSelList <- list(nSelList)
}
} else {
nSelList <- nIndNow %/% 10
message(paste0("`nSelList` is not specified. We substitute `nSelList = list(",
nSelList,")` instead."))
nSelList <- sapply(X = nSelectionWaysVec,
FUN = function(nSelectionWays) {
rep(nSelList, nSelectionWays)
}, simplify = FALSE)
}
if (!(length(nSelList) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(nSelList) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(nSelList) == 1) {
nSelList <- rep(nSelList, nGenerationProceedSimulation)
}
stopifnot(all(unlist(lapply(nSelList, length)) == nSelectionWaysVec))
stopifnot(all(unlist(lapply(nSelList, is.numeric))))
nSelList <- sapply(X = 1:nGenerationProceedSimulation,
FUN = function(generationProceedNo) {
nSel <- nSelList[[generationProceedNo]]
whereSelection <- whereSelectionList[[generationProceedNo]]
nSel[!whereSelection] <- nIndNow
return(nSel)
}, simplify = FALSE)
names(nSelList) <- 1:nGenerationProceedSimulation
# nSelInitOPVList
if (!is.null(nSelInitOPVList)) {
if (!is.list(nSelInitOPVList)) {
nSelInitOPVList <- list(nSelInitOPVList)
}
} else {
nSelInitOPVList <- nIndNow %/% 2
if (any(unlist(lapply(selectionMethodList, function(selectionMethod) "selectOPV" %in% selectionMethod)))) {
message(paste0("`nSelInitOPVList` is not specified. We substitute `nSelInitOPVList = list(",
nSelInitOPVList,")` instead."))
}
nSelInitOPVList <- sapply(X = nSelectionWaysVec,
FUN = function(nSelectionWays) {
rep(nSelInitOPVList, nSelectionWays)
}, simplify = FALSE)
}
if (!(length(nSelInitOPVList) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(nSelInitOPVList) must be equal to 1 or equal to nGenerationProceedSimulation"))
} else if (length(nSelInitOPVList) == 1) {
nSelInitOPVList <- rep(nSelInitOPVList, nGenerationProceedSimulation)
}
stopifnot(all(unlist(lapply(nSelInitOPVList, length)) == nSelectionWaysVec))
stopifnot(all(unlist(lapply(nSelInitOPVList, is.numeric))))
stopifnot(all(unlist(lapply(nSelInitOPVList, is.numeric))))
whereNSelSatisfy <- unlist(mapply(FUN = function(nSelInitOPV, nSel) {
all(nSelInitOPV >= nSel)
},
nSelInitOPVList,
nSelList))
if (any(!whereNSelSatisfy)) {
if (any(lapply(selectionMethodList, function(selectionMethod) "selectOPV" %in% selectionMethod))) {
message("`nSelInitOPV` should be larger than `nSel`. We substitute `nSelInitOPV` by `nSel` when `nSelInitOPV` is smaller than `nSel`.")
}
nSelInitOPVList[whereNSelSatisfy] <- sapply(X = (1:nGenerationProceedSimulation)[whereNSelSatisfy],
FUN = function(generationProceedNo) {
nSelInitOPV <- nSelInitOPVList[generationProceedNo]
nSel <- nSelList[generationProceedNo]
nSelInitOPV[nSelInitOPV < nSel] <- nSel[nSelInitOPV < nSel]
return(nSelInitOPV)
}, simplify = FALSE)
}
names(nSelInitOPVList) <- 1:nGenerationProceedSimulation
# nClusterList
if (!is.null(nClusterList)) {
if (!is.list(nClusterList)) {
nClusterList <- list(nClusterList)
}
} else {
nClusterList <- 5
message(paste0("`nClusterList` is not specified. We substitute `nClusterList = list(",
nClusterList,")` instead."))
nClusterList <- sapply(X = nSelectionWaysVec,
FUN = function(nSelectionWays) {
rep(nClusterList, nSelectionWays)
}, simplify = FALSE)
}
if (!(length(nClusterList) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(nClusterList) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(nClusterList) == 1) {
nClusterList <- rep(nClusterList, nGenerationProceedSimulation)
}
stopifnot(all(unlist(lapply(nClusterList, length)) == nSelectionWaysVec))
stopifnot(all(unlist(lapply(nClusterList, is.numeric))))
nClusterList <- sapply(X = 1:nGenerationProceedSimulation,
FUN = function(generationProceedNo) {
nCluster <- nClusterList[[generationProceedNo]]
clusteringForSel <- clusteringForSelList[[generationProceedNo]]
nCluster[!clusteringForSel] <- 1
return(nCluster)
}, simplify = FALSE)
names(nClusterList) <- 1:nGenerationProceedSimulation
# nTopClusterList
if (!is.null(nTopClusterList)) {
if (!is.list(nTopClusterList)) {
nTopClusterList <- list(nTopClusterList)
}
} else {
nTopClusterList <- 5
message(paste0("`nTopClusterList` is not specified. We substitute `nTopClusterList = list(",
nTopClusterList,")` instead."))
nTopClusterList <- sapply(X = nSelectionWaysVec,
FUN = function(nSelectionWays) {
rep(nTopClusterList, nSelectionWays)
}, simplify = FALSE)
}
if (!(length(nTopClusterList) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(nTopClusterList) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(nTopClusterList) == 1) {
nTopClusterList <- rep(nTopClusterList, nGenerationProceedSimulation)
}
stopifnot(all(unlist(lapply(nTopClusterList, length)) == nSelectionWaysVec))
stopifnot(all(unlist(lapply(nTopClusterList, is.numeric))))
nTopClusterList <- sapply(X = 1:nGenerationProceedSimulation,
FUN = function(generationProceedNo) {
nTopCluster <- nTopClusterList[[generationProceedNo]]
clusteringForSel <- clusteringForSelList[[generationProceedNo]]
nTopCluster[!clusteringForSel] <- 1
return(nTopCluster)
}, simplify = FALSE)
names(nTopClusterList) <- 1:nGenerationProceedSimulation
# nTopEachList
if (!is.null(nTopEachList)) {
if (!is.list(nTopEachList)) {
nTopEachList <- list(nTopEachList)
}
} else {
nTopEachList <- sapply(X = 1:nGenerationProceedSimulation,
FUN = function(generationProceedNo) {
nTopEach <- nTopEachList[[generationProceedNo]]
nSel <- nSelList[[generationProceedNo]]
nTopCluster <- nTopClusterList[[generationProceedNo]]
nTopEach <- nSel %/% nTopCluster
return(nTopEach)
}, simplify = FALSE)
message(paste0("`nTopEachList` is not specified. We substitute `nTopEachList = list(",
nTopEachList,")` instead."))
}
if (!(length(nTopEachList) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(nTopEachList) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(nTopEachList) == 1) {
nTopEachList <- rep(nTopEachList, nGenerationProceedSimulation)
}
stopifnot(all(unlist(lapply(nTopEachList, length)) == nSelectionWaysVec))
stopifnot(all(unlist(lapply(nTopEachList, is.numeric))))
names(nTopEachList) <- 1:nGenerationProceedSimulation
# multiTraitsEvalMethodList
if (!is.null(multiTraitsEvalMethodList)) {
if (!is.list(multiTraitsEvalMethodList)) {
multiTraitsEvalMethodList <- list(multiTraitsEvalMethodList)
}
} else {
multiTraitsEvalMethodList <- "sum"
message(paste0("`multiTraitsEvalMethodList` is not specified. We substitute `multiTraitsEvalMethodList = list(",
multiTraitsEvalMethodList,")` instead."))
multiTraitsEvalMethodList <- sapply(X = nSelectionWaysVec,
FUN = function(nSelectionWays) {
rep(multiTraitsEvalMethodList, nSelectionWays)
}, simplify = FALSE)
}
if (!(length(multiTraitsEvalMethodList) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(multiTraitsEvalMethodList) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(multiTraitsEvalMethodList) == 1) {
multiTraitsEvalMethodList <- rep(multiTraitsEvalMethodList, nGenerationProceedSimulation)
}
names(multiTraitsEvalMethodList) <- 1:nGenerationProceedSimulation
stopifnot(all(unlist(lapply(multiTraitsEvalMethodList, length)) == nSelectionWaysVec))
stopifnot(all(unlist(lapply(multiTraitsEvalMethodList, function(x) all(x %in% multiTraitsEvalMethodsOffered)))))
# hSelList
if (!is.null(hSelList)) {
if (!is.list(hSelList)) {
hSelList <- sapply(nSelectionWaysVec,
function(nSelectionWays) {
hSelListNow <- rep(list(hSelList), nSelectionWays)
return(hSelListNow)
}, simplify = FALSE)
} else if (!is.list(hSelList[[1]])) {
hSelList <- sapply(nSelectionWaysVec,
function(nSelectionWays) {
hSelListNow <- rep(hSelList, nSelectionWays)
return(hSelListNow)
}, simplify = FALSE)
}
} else {
hSelList <- 1
message(paste0("`hSelList` is not specified. We substitute `hSelList = list(list(",
hSelList,"))` instead."))
hSelList <- sapply(1:nGenerationProceedSimulation,
function(generationProceedNo) {
hSelListNow <- sapply(X = traitNoSelList[[generationProceedNo]],
FUN = function(traitNoSelNow) {
rep(hSelList, length(traitNoSelNow))
}, simplify = FALSE)
return(hSelListNow)
}, simplify = FALSE)
}
if (!(length(hSelList) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(hSelList) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(hSelList) == 1) {
hSelList <- rep(hSelList, nGenerationProceedSimulation)
}
stopifnot(all(unlist(lapply(hSelList, is.list))))
stopifnot(all(unlist(lapply(hSelList, function(hSel) all(unlist(lapply(hSel, is.numeric)))))))
stopifnot(all(sapply(hSelList, function(hSel) all(unlist(lapply(hSel, function(x) all(x >= 0)))))))
stopifnot(all(unlist(lapply(hSelList, length)) == nSelectionWaysVec))
names(hSelList) <- 1:nGenerationProceedSimulation
stopifnot(all(unlist(lapply(X = hSelList, FUN = function(x) lapply(x, length))) ==
unlist(lapply(X = traitNoSelList, FUN = function(x) lapply(x, length)))))
# matingMethodVec
if (!is.null(matingMethodVec)) {
stopifnot(is.character(matingMethodVec))
stopifnot(all(matingMethodVec %in% matingMethodsOffered))
} else {
matingMethodVec <- "randomMate"
message(paste0("`matingMethodVec` is not specified. We substitute `matingMethodVec = ",
matingMethodVec,"` instead."))
}
if (!(length(matingMethodVec) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(matingMethodVec) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(matingMethodVec) == 1) {
matingMethodVec <- rep(matingMethodVec, nGenerationProceedSimulation)
}
names(matingMethodVec) <- 1:nGenerationProceedSimulation
# allocateMethodVec
if (!is.null(allocateMethodVec)) {
stopifnot(is.character(allocateMethodVec))
stopifnot(all(allocateMethodVec %in% allocateMethodsOffered))
} else {
allocateMethodVec <- "equalAllocation"
message(paste0("`allocateMethodVec` is not specified. We substitute `allocateMethodVec = ",
allocateMethodVec,"` instead."))
}
if (!(length(allocateMethodVec) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(allocateMethodVec) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(allocateMethodVec) == 1) {
allocateMethodVec <- rep(allocateMethodVec, nGenerationProceedSimulation)
}
names(allocateMethodVec) <- 1:nGenerationProceedSimulation
# weightedAllocationMethodList
if (!is.null(weightedAllocationMethodList)) {
if (!is.list(weightedAllocationMethodList)) {
weightedAllocationMethodList <- list(weightedAllocationMethodList)
}
} else {
weightedAllocationMethodList <- "selectBV"
message(paste0("`weightedAllocationMethodList` is not specified. We substitute `weightedAllocationMethodList = list(",
weightedAllocationMethodList,")` instead."))
weightedAllocationMethodList <- list(weightedAllocationMethodList)
}
if (!(length(weightedAllocationMethodList) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(weightedAllocationMethodList) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(weightedAllocationMethodList) == 1) {
weightedAllocationMethodList <- rep(weightedAllocationMethodList, nGenerationProceedSimulation)
}
weightedAllocationMethodList <- lapply(weightedAllocationMethodList, function(x) x[x %in% selectionMethodsWithSelection])
names(weightedAllocationMethodList) <- 1:nGenerationProceedSimulation
stopifnot(all(unlist(lapply(weightedAllocationMethodList, is.character))))
stopifnot(all(unlist(lapply(weightedAllocationMethodList, length)) >= 1))
# traitNoRAList
if (!is.null(traitNoRAList)) {
if (!is.list(traitNoRAList)) {
traitNoRAList <- list(traitNoRAList)
}
} else {
traitNoRAList <- 1
message(paste0("`traitNoRAList` is not specified. We substitute `traitNoRAList = list(",
traitNoRAList,")` instead."))
traitNoRAList <- list(traitNoRAList)
}
if (!(length(traitNoRAList) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(traitNoRAList) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(traitNoRAList) == 1) {
traitNoRAList <- rep(traitNoRAList, nGenerationProceedSimulation)
}
stopifnot(all(unlist(lapply(traitNoRAList, is.numeric))))
stopifnot(all(sapply(traitNoRAList, function(traitNoRA) all(traitNoRA >= 1))))
names(traitNoRAList) <- 1:nGenerationProceedSimulation
# minimumUnitAllocateVec
if (!is.null(minimumUnitAllocateVec)) {
stopifnot(is.numeric(minimumUnitAllocateVec))
minimumUnitAllocateVec <- floor(minimumUnitAllocateVec)
stopifnot(all(minimumUnitAllocateVec >= 1))
} else {
minimumUnitAllocateVec <- 1
message(paste0("`minimumUnitAllocateVec` is not specified. We substitute `minimumUnitAllocateVec = ",
minimumUnitAllocateVec,"` instead."))
}
if (!(length(minimumUnitAllocateVec) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(minimumUnitAllocateVec) must be equal to 1 or equal to nGenerationProceedSimulation"))
} else if (length(minimumUnitAllocateVec) == 1) {
minimumUnitAllocateVec <- rep(minimumUnitAllocateVec, nGenerationProceedSimulation)
}
names(minimumUnitAllocateVec) <- 1:nGenerationProceedSimulation
# includeGVPVec
stopifnot(is.logical(includeGVPVec))
if (!(length(includeGVPVec) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(includeGVPVec) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(includeGVPVec) == 1) {
includeGVPVec <- rep(includeGVPVec, nGenerationProceedSimulation)
}
names(includeGVPVec) <- 1:nGenerationProceedSimulation
# targetGenGVPVec
if (!is.null(targetGenGVPVec)) {
stopifnot(is.numeric(targetGenGVPVec))
targetGenGVPVec <- floor(targetGenGVPVec)
stopifnot(all(targetGenGVPVec >= 0))
} else {
targetGenGVPVec <- 0
message(paste0("`targetGenGVPVec` is not specified. We substitute `targetGenGVPVec = ",
targetGenGVPVec,"` instead."))
}
if (!(length(targetGenGVPVec) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(targetGenGVPVec) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(targetGenGVPVec) == 1) {
targetGenGVPVec <- rep(targetGenGVPVec, nGenerationProceedSimulation)
}
names(targetGenGVPVec) <- 1:nGenerationProceedSimulation
# performOCSVec
stopifnot(is.logical(performOCSVec))
if (!(length(performOCSVec) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(performOCSVec) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(performOCSVec) == 1) {
performOCSVec <- rep(performOCSVec, nGenerationProceedSimulation)
}
names(performOCSVec) <- 1:nGenerationProceedSimulation
# targetGenOCSVec
if (!is.null(targetGenOCSVec)) {
stopifnot(is.numeric(targetGenOCSVec))
targetGenOCSVec <- floor(targetGenOCSVec)
stopifnot(all(targetGenOCSVec >= 1))
} else {
targetGenOCSVec <- nGenerationProceedSimulation
message(paste0("`targetGenOCSVec` is not specified. We substitute `targetGenOCSVec = ",
targetGenOCSVec,"` instead."))
}
if (!(length(targetGenOCSVec) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(targetGenOCSVec) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(targetGenOCSVec) == 1) {
targetGenOCSVec <- rep(targetGenOCSVec, nGenerationProceedSimulation)
}
names(targetGenOCSVec) <- 1:nGenerationProceedSimulation
# HeStarRatioVec
if (!is.null(HeStarRatioVec)) {
stopifnot(is.numeric(HeStarRatioVec))
stopifnot(all(HeStarRatioVec > 0))
} else {
HeStarRatioVec <- 0.1
message(paste0("`HeStarRatioVec` is not specified. We substitute `HeStarRatioVec = ",
HeStarRatioVec,"` instead."))
}
if (!(length(HeStarRatioVec) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(HeStarRatioVec) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(HeStarRatioVec) == 1) {
HeStarRatioVec <- rep(HeStarRatioVec, nGenerationProceedSimulation)
}
names(HeStarRatioVec) <- 1:nGenerationProceedSimulation
# degreeOCSVec
if (!is.null(degreeOCSVec)) {
stopifnot(is.numeric(degreeOCSVec))
stopifnot(all(degreeOCSVec > 0))
} else {
degreeOCSVec <- 1
message(paste0("`degreeOCSVec` is not specified. We substitute `degreeOCSVec = ",
degreeOCSVec,"` instead."))
}
if (!(length(degreeOCSVec) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(degreeOCSVec) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(degreeOCSVec) == 1) {
degreeOCSVec <- rep(degreeOCSVec, nGenerationProceedSimulation)
}
names(degreeOCSVec) <- 1:nGenerationProceedSimulation
# includeGVPOCSVec
stopifnot(is.logical(includeGVPOCSVec))
if (!(length(includeGVPOCSVec) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(includeGVPOCSVec) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(includeGVPOCSVec) == 1) {
includeGVPOCSVec <- rep(includeGVPOCSVec, nGenerationProceedSimulation)
}
names(includeGVPOCSVec) <- 1:nGenerationProceedSimulation
# hOCSList
if (!is.null(hOCSList)) {
if (!is.list(hOCSList)) {
hOCSList <- list(hOCSList)
}
} else {
hOCSList <- 1
message(paste0("`hList` is not specified. We substitute `hList = list(",
hList,")` instead."))
hOCSList <- sapply(X = 1:nGenerationProceedSimulation,
FUN = function(generationProceedNo) {
hOCSLenNow <- length(traitNoOCSList[[generationProceedNo]]) *
(length(weightedAllocationMethodOCSList[[generationProceedNo]]) +
includeGVPOCSVec[generationProceedNo])
return(rep(hOCSList, hLenOCSNow))
}, simplify = FALSE)
}
if (!(length(hOCSList) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(hOCSList) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(hOCSList) == 1) {
hOCSList <- rep(hOCSList, nGenerationProceedSimulation)
}
stopifnot(all(unlist(lapply(hOCSList, is.numeric))))
# stopifnot(all(sapply(hOCSList, function(h) all(h >= 0))))
stopifnot(all(sapply(hOCSList, function(h) all(h <= 10))))
stopifnot(all(unlist(lapply(hOCSList, length)) == (unlist(lapply(traitNoOCSList, length)) *
(unlist(lapply(weightedAllocationMethodOCSList, length)) + includeGVPOCSVec))))
names(hOCSList) <- 1:nGenerationProceedSimulation
# weightedAllocationMethodOCSList
if (!is.null(weightedAllocationMethodOCSList)) {
if (!is.list(weightedAllocationMethodOCSList)) {
weightedAllocationMethodOCSList <- list(weightedAllocationMethodOCSList)
}
} else {
weightedAllocationMethodOCSList <- "selectBV"
message(paste0("`weightedAllocationMethodOCSList` is not specified. We substitute `weightedAllocationMethodOCSList = list(",
weightedAllocationMethodOCSList,")` instead."))
weightedAllocationMethodOCSList <- list(weightedAllocationMethodOCSList)
}
if (!(length(weightedAllocationMethodOCSList) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(weightedAllocationMethodOCSList) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(weightedAllocationMethodOCSList) == 1) {
weightedAllocationMethodOCSList <- rep(weightedAllocationMethodOCSList, nGenerationProceedSimulation)
}
weightedAllocationMethodOCSList <- lapply(weightedAllocationMethodOCSList, function(x) x[x %in% selectionMethodsWithSelection])
names(weightedAllocationMethodOCSList) <- 1:nGenerationProceedSimulation
stopifnot(all(unlist(lapply(weightedAllocationMethodOCSList, is.character))))
stopifnot(all(unlist(lapply(weightedAllocationMethodOCSList, length)) >= 1))
# traitNoOCSList
if (!is.null(traitNoOCSList)) {
if (!is.list(traitNoOCSList)) {
traitNoOCSList <- list(traitNoOCSList)
}
} else {
traitNoOCSList <- 1
message(paste0("`traitNoOCSList` is not specified. We substitute `traitNoRAList = list(",
traitNoOCSList,")` instead."))
traitNoOCSList <- list(traitNoOCSList)
}
if (!(length(traitNoOCSList) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(traitNoOCSList) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(traitNoOCSList) == 1) {
traitNoOCSList <- rep(traitNoOCSList, nGenerationProceedSimulation)
}
stopifnot(all(unlist(lapply(traitNoOCSList, is.numeric))))
stopifnot(all(sapply(traitNoOCSList, function(traitNoOCS) all(traitNoOCS >= 1))))
names(traitNoOCSList) <- 1:nGenerationProceedSimulation
# nCrossesOCSVec
if (!is.null(nCrossesOCSVec)) {
stopifnot(is.numeric(nCrossesOCSVec))
nCrossesOCSVec <- floor(nCrossesOCSVec)
stopifnot(all(nCrossesOCSVec >= 1))
} else {
nCrossesOCSVec <- 0
message(paste0("`nCrossesOCSVec` is not specified. We substitute `nCrossesOCSVec = ",
nCrossesOCSVec,"` instead."))
}
if (!(length(nCrossesOCSVec) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(nCrossesOCSVec) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(nCrossesOCSVec) == 1) {
nCrossesOCSVec <- rep(nCrossesOCSVec, nGenerationProceedSimulation)
}
names(nCrossesOCSVec) <- 1:nGenerationProceedSimulation
# propSelOpt
stopifnot(is.logical(propSelOpt))
# hLens
propSelLens <- propSelOpt * nSelectionWaysVec
hOriLens <- sapply(X = 1:nGenerationProceedSimulation,
FUN = function(genNow) {
length(traitNoRAList[[genNow]]) * (includeGVPVec[genNow] + length(weightedAllocationMethodList[[genNow]]))
}, simplify = TRUE)
hLens <- propSelLens + hOriLens
propSelOrNot <- unlist(sapply(X = 1:nGenerationProceedSimulation,
FUN = function(genNow) {
c(rep(TRUE, propSelLens[genNow]),
rep(FALSE, hOriLens[genNow]))
}, simplify = FALSE))
hNamesList <- sapply(X = 1:nGenerationProceedSimulation,
FUN = function(genNow) {
selectNames <- paste0("Selection-", 1:nSelectionWaysVec[genNow])
hNamesOri <- paste(bsInfoInit$traitInfo$traitNames[traitNoRAList[[genNow]]],
c(weightedAllocationMethodList[[genNow]], "genVarProgeny"), sep = "-")[1:hLens[genNow]]
if (propSelOpt) {
return(c(selectNames, hNamesOri))
} else {
return(hNamesOri)
}
}, simplify = FALSE)
names(hNamesList) <- paste0("Generation_", 1:nGenerationProceedSimulation)
hVecNamesAll <- paste0(rep(names(hNamesList), hLens), "-", unlist(hNamesList, use.names = FALSE))
# sameAcrossGeneration
stopifnot(is.logical(sameAcrossGeneration))
if (!((length(unique(traitNoRAList)) == 1) & (length(unique(weightedAllocationMethodList)) == 1))) {
sameAcrossGeneration <- FALSE
message(paste0("You cannot set `sameAcrossGeneration` as `TRUE` with this `traitNoRAList` and `weightedAllocationMethodList`. We substitute `sameAcrossGeneration = ",
sameAcrossGeneration,"` instead."))
}
hVecLen <- ifelse(sameAcrossGeneration, max(hLens), sum(hLens))
if (sameAcrossGeneration) {
hVecNames <- hNamesList[[which.max(hLens)]]
} else {
hVecNames <- hVecNamesAll
}
# hMin
if (!is.null(hMin)) {
stopifnot(is.numeric(hMin))
# stopifnot(hMin >= 0)
} else {
hMin <- 0
message(paste0("`hMin` is not specified. We substitute `hMin = ",
hMin,"` instead."))
}
if (!(length(hMin) %in% c(1, hVecLen))) {
stop(paste0("length(hMin) must be equal to 1 or equal to ", hVecLen, "."))
} else if (length(hMin) == 1) {
hMin <- rep(hMin, hVecLen)
}
if (propSelOpt) {
hMin[propSelOrNot] <- 0
}
# hMax
if (!is.null(hMax)) {
stopifnot(is.numeric(hMax))
stopifnot(hMax > 0)
} else {
hMax <- 2
message(paste0("`hMax` is not specified. We substitute `hMax = ",
hMax,"` instead."))
}
if (!(length(hMax) %in% c(1, hVecLen))) {
stop(paste0("length(hMax) must be equal to 1 or equal to ", hVecLen, "."))
} else if (length(hMax) == 1) {
hMax <- rep(hMax, hVecLen)
}
if (propSelOpt) {
hMax[propSelOrNot] <- 1
}
# hStart
hStart <- (hMin + hMax) / 2
names(hMin) <- names(hMax) <- names(hStart) <- hVecNames
# nNextPopVec
if (!is.null(nNextPopVec)) {
stopifnot(is.numeric(nNextPopVec))
nNextPopVec <- floor(nNextPopVec)
stopifnot(all(nNextPopVec >= 1))
} else {
nNextPopVec <- nIndNow
message(paste0("`nNextPopVec` is not specified. We substitute `nNextPopVec = ",
nNextPopVec,"` instead."))
}
if (!(length(nNextPopVec) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(nNextPopVec) must be equal to 1 or equal to nGenerationProceedSimulation."))
} else if (length(nNextPopVec) == 1) {
nNextPopVec <- rep(nNextPopVec, nGenerationProceedSimulation)
}
names(nNextPopVec) <- 1:nGenerationProceedSimulation
# nProgeniesEMBV
if (!is.null(nProgeniesEMBVVec)) {
stopifnot(is.numeric(nProgeniesEMBVVec))
nProgeniesEMBVVec <- floor(nProgeniesEMBVVec)
stopifnot(all(nProgeniesEMBVVec >= 1))
} else {
nProgeniesEMBVVec <- round(2 * nNextPopVec / min(unlist(lapply(nSelList, sum)), nNextPopVec))
if (any(unlist(lapply(selectionMethodList, function(selectionMethod) "selectEMBV" %in% selectionMethod)))) {
message(paste0("`nProgeniesEMBVVec` is not specified. We substitute `nProgeniesEMBVVec = c(",
paste(nProgeniesEMBVVec, collapse = ", "),
")` instead."))
}
}
if (!(length(nProgeniesEMBVVec) %in% c(1, nGenerationProceedSimulation))) {
stop(paste("length(nProgeniesEMBVVec) must be equal to 1 or equal to nGenerationProceedSimulation"))
} else if (length(nProgeniesEMBVVec) == 1) {
nProgeniesEMBVVec <- rep(nProgeniesEMBVVec, nGenerationProceedSimulation)
}
names(nProgeniesEMBVVec) <- 1:nGenerationProceedSimulation
# nIterEMBV
if (!is.null(nIterEMBV)) {
stopifnot(is.numeric(nIterEMBV))
nIterEMBV <- floor(nIterEMBV)
stopifnot(nIterEMBV >= 1)
} else {
nIterEMBV <- 10
if (any(unlist(lapply(selectionMethodList, function(selectionMethod) "selectEMBV" %in% selectionMethod)))) {
message(paste0("`nIterEMBV` is not specified. We substitute `nIterEMBV = ", nIterEMBV,"` instead."))
}
}
# nCoresEMBV
if (!is.null(nCoresEMBV)) {
stopifnot(is.numeric(nCoresEMBV))
nCoresEMBV <- floor(nCoresEMBV)
stopifnot(nCoresEMBV >= 1)
} else {
nCoresEMBV <- 1
if (any(unlist(lapply(selectionMethodList, function(selectionMethod) "selectEMBV" %in% selectionMethod)))) {
message(paste0("`nCoresEMBV` is not specified. We substitute `nCoresEMBV = ",
nCoresEMBV,"` instead."))
}
}
if (nCoresEMBV >= parallel::detectCores()) {
warning("You are going to assign the number of cores larger than that of your PC to `nCoresEMBV` ! Is it OK ?")
}
# nameMethod
if (!is.null(nameMethod)) {
if (!(nameMethod %in% nameMethodsOffered)) {
stop(paste0("We only offer the following methods for naming individuals: ",
paste(nameMethodsOffered, collapse = "; ")))
}
}
# returnMethod
if (!is.null(returnMethod)) {
if (!all(returnMethod %in% returnMethodsOffered)) {
returnMethod <- returnMethod[returnMethod %in% returnMethodsOffered]
message(paste0("We only offer the following methods for returining the simulation results: ",
paste(returnMethodsOffered, collapse = "; ")))
stopifnot(length(returnMethod) >= 1)
}
} else {
returnMethod <- "summary"
}
# evaluateGVMethod
if (!is.null(evaluateGVMethod)) {
if (!(evaluateGVMethod %in% lociEffMethodsOffered)) {
stop(paste0("We only offer the following methods for evaluating the simulation results: ",
paste(lociEffMethodsOffered, collapse = "; ")))
}
} else {
evaluateGVMethod <- "true"
}
# nTopEvalForOpt
if (!is.null(nTopEvalForOpt)) {
stopifnot(is.numeric(nTopEvalForOpt))
nTopEvalForOpt <- floor(nTopEvalForOpt)
stopifnot(nTopEvalForOpt >= 1)
stopifnot(nTopEvalForOpt <= min(nNextPopVec))
} else {
nTopEvalForOpt <- 1
message(paste0("`nTopEvalForOpt` is not specified. We substitute `nTopEvalForOpt = ",
nTopEvalForOpt,"` instead."))
}
# nTopEval
if (!is.null(nTopEval)) {
stopifnot(is.numeric(nTopEval))
nTopEval <- floor(nTopEval)
stopifnot(nTopEval >= 1)
stopifnot(nTopEval <= min(nNextPopVec))
} else {
nTopEval <- 1
message(paste0("`nTopEval` is not specified. We substitute `nTopEval = ",
nTopEval,"` instead."))
}
# traitNoEval
if (!is.null(traitNoEval)) {
stopifnot(is.numeric(traitNoEval))
traitNoEval <- floor(traitNoEval)
stopifnot(all(traitNoEval >= 1))
stopifnot(all(traitNoEval <= nTraits))
} else {
traitNoEval <- 1
message(paste0("`traitNoEval` is not specified. We substitute `traitNoEval = ", traitNoEval,"` instead."))
}
# hEval
if (!is.null(hEval)) {
stopifnot(is.numeric(hEval))
stopifnot(all(hEval >= 0))
} else {
hEval <- 0.1
message(paste0("`hEval` is not specified. We substitute `hEval = ", hEval, "` instead."))
}
hEvalLen <- length(traitNoEval)
if (!(length(hEval) %in% c(1, hEvalLen))) {
stop(paste("length(hEval) must be equal to 1 or equal to `length(traitNoEval)`"))
} else if (length(hEval) == 1) {
hEval <- rep(hEval, hEvalLen)
}
# nCores
if (!is.null(nCores)) {
stopifnot(is.numeric(nCores))
nCores <- floor(nCores)
stopifnot(nCores >= 1)
} else {
nCores <- 1
message(paste0("`nCores` is not specified. We substitute `nCores = ",
nCores,"` instead."))
}
if (nCores >= parallel::detectCores()) {
warning("You are going to assign the number of cores larger than that of your PC to `nCores` ! Is it OK ?")
}
# nCoresPerOptimization
if (!is.null(nCoresPerOptimization)) {
stopifnot(is.numeric(nCoresPerOptimization))
nCoresPerOptimization <- floor(nCoresPerOptimization)
stopifnot(nCoresPerOptimization >= 1)
} else {
nCoresPerOptimization <- 1
message(paste0("`nCoresPerOptimization` is not specified. We substitute `nCoresPerOptimization = ",
nCoresPerOptimization,"` instead."))
}
if (nCoresPerOptimization >= parallel::detectCores()) {
warning("You are going to assign the number of cores larger than that of your PC to `nCoresPerOptimization` ! Is it OK ?")
}
estimatedGVInitExist <- !is.null(breederInfoInit$estimatedGVByMLRInfo[[names(bsInfoInit$populations[length(bsInfoInit$populations)])]])
if (!estimatedGVInitExist) {
# trainingPopInit
if (is.null(trainingPopInit)) {
if (trainingPopType == "all") {
trainingPopInit <- 1:length(breederInfoInit$populationsFB)
} else {
trainingPopInit <- length(breederInfoInit$populationsFB)
}
} else {
stopifnot(is.numeric(trainingPopInit))
trainingPopInit <- floor(trainingPopInit)
if (!all(trainingPopInit %in% (1:bsInfoInit$generation))) {
message(paste0("The following population cannot be used for initial training population; ",
paste(trainingPopInit[!(trainingPopInit %in% (1:bsInfoInit$generation))], collapse = ", ")))
trainingPopInit <- trainingPopInit[trainingPopInit %in% (1:bsInfoInit$generation)]
}
}
# trainingIndNamesInit
if (!is.null(trainingIndNamesInit)) {
stopifnot(is.character(trainingIndNamesInit))
stopifnot(length(trainingIndNamesInit) >= 1)
}
}
# saveAllResAt
if (!is.null(saveAllResAt)) {
stopifnot(is.character(saveAllResAt))
tailCheck <- TRUE
while (tailCheck) {
strLength <- stringr::str_length(string = saveAllResAt)
tailCheck <- stringr::str_sub(string = saveAllResAt, start = strLength, end = strLength) == "/"
if (tailCheck) {
saveAllResAt <- stringr::str_sub(string = saveAllResAt, start = 1, end = strLength - 1)
}
}
if (!dir.exists(paths = saveAllResAt)) {
dir.create(path = saveAllResAt)
}
}
# summaryAllResAt
if (!is.null(summaryAllResAt)) {
stopifnot(is.character(summaryAllResAt))
tailCheck <- TRUE
while (tailCheck) {
strLength <- stringr::str_length(string = summaryAllResAt)
tailCheck <- stringr::str_sub(string = summaryAllResAt, start = strLength, end = strLength) == "/"
if (tailCheck) {
summaryAllResAt <- stringr::str_sub(string = summaryAllResAt, start = 1, end = strLength - 1)
}
}
if (!dir.exists(paths = summaryAllResAt)) {
message(paste0("There is no directory named '", summaryAllResAt, "'. We will summarize the simulation results inside the `simBsCma` object."))
summaryAllResAt <- NULL
}
}
# Save arguments in `self`
self$simBsName <- simBsName
self$bsInfoInit <- bsInfoInit
self$breederInfoInit <- breederInfoInit
self$lociEffMethod <- lociEffMethod
self$trainingPopType <- trainingPopType
self$trainingPopInit <- trainingPopInit
self$trainingIndNamesInit <- trainingIndNamesInit
self$methodMLRInit <- methodMLRInit
self$multiTraitInit <- multiTraitInit
self$nIterSimulation <- nIterSimulation
self$nGenerationProceed <- nGenerationProceed
self$nGenerationProceedSimulation <- nGenerationProceedSimulation
self$setGoalAsFinalGeneration <- setGoalAsFinalGeneration
self$performOptimization <- performOptimization
self$useFirstOptimizedValue <- useFirstOptimizedValue
self$performRobustOptimization <- performRobustOptimization
self$sameAcrossGeneration <- sameAcrossGeneration
self$lowerQuantile <- lowerQuantile
self$hMin <- hMin
self$hMax <- hMax
self$nTotalIterForOneOptimization <- nTotalIterForOneOptimization
self$nIterSimulationPerEvaluation <- nIterSimulationPerEvaluation
self$nIterSimulationForOneMrkEffect <- nIterSimulationForOneMrkEffect
self$nIterMrkEffectForRobustOptimization <- nIterMrkEffectForRobustOptimization
self$nIterOptimization <- nIterOptimization
self$lambda <- lambda
self$mu <- mu
self$saveObjNameBase <- saveObjNameBase
self$whenToSaveObj <- whenToSaveObj
self$nTopEvalForOpt <- nTopEvalForOpt
self$rewardWeightVec <- rewardWeightVec
self$digitsEval <- digitsEval
self$nRefreshMemoryEvery <- nRefreshMemoryEvery
self$propSelOpt <- propSelOpt
self$propSelOrNot <- propSelOrNot
self$hOriLens <- hOriLens
self$propSelLens <- propSelLens
self$hLens <- hLens
self$hStart <- hStart
self$updateBreederInfo <- updateBreederInfo
self$phenotypingInds <- phenotypingInds
self$nRepForPhenoInit <- nRepForPhenoInit
self$nRepForPheno <- nRepForPheno
self$updateModels <- updateModels
self$updateBreederInfoSimulation <- updateBreederInfoSimulation
self$phenotypingIndsSimulation <- phenotypingIndsSimulation
self$nRepForPhenoSimulation <- nRepForPhenoSimulation
self$updateModelsSimulation <- updateModelsSimulation
self$methodMLR <- methodMLR
self$multiTrait <- multiTrait
self$nSelectionWaysVec <- nSelectionWaysVec
self$selectionMethodList <- selectionMethodList
self$traitNoSelList <- traitNoSelList
self$blockSplitMethod <- blockSplitMethod
self$nMrkInBlock <- nMrkInBlock
self$minimumSegmentLength <- minimumSegmentLength
self$nSelInitOPVList <- nSelInitOPVList
self$nIterOPV <- nIterOPV
self$nProgeniesEMBVVec <- nProgeniesEMBVVec
self$nIterEMBV <- nIterEMBV
self$nCoresEMBV <- nCoresEMBV
self$clusteringForSelList <- clusteringForSelList
self$nClusterList <- nClusterList
self$nTopClusterList <- nTopClusterList
self$nTopEachList <- nTopEachList
self$nSelList <- nSelList
self$multiTraitsEvalMethodList <- multiTraitsEvalMethodList
self$hSelList <- hSelList
self$matingMethodVec <- matingMethodVec
self$allocateMethodVec <- allocateMethodVec
self$weightedAllocationMethodList <- weightedAllocationMethodList
self$traitNoRAList <- traitNoRAList
self$minimumUnitAllocateVec <- minimumUnitAllocateVec
self$includeGVPVec <- includeGVPVec
self$nNextPopVec <- nNextPopVec
self$targetGenGVPVec <- targetGenGVPVec
self$performOCSVec <- performOCSVec
self$targetGenOCSVec <- targetGenOCSVec
self$HeStarRatioVec <- HeStarRatioVec
self$degreeOCSVec <- degreeOCSVec
self$includeGVPOCSVec <- includeGVPOCSVec
self$hOCSList <- hOCSList
self$weightedAllocationMethodOCSList <- weightedAllocationMethodOCSList
self$traitNoOCSList <- traitNoOCSList
self$nCrossesOCSVec <- nCrossesOCSVec
self$nameMethod <- nameMethod
self$nCores <- nCores
self$nCoresPerOptimization <- nCoresPerOptimization
self$overWriteRes <- overWriteRes
self$showProgress <- showProgress
self$returnMethod <- returnMethod
self$saveAllResAt <- saveAllResAt
self$evaluateGVMethod <- evaluateGVMethod
self$nTopEval <- nTopEval
self$traitNoEval <- traitNoEval
self$hEval <- hEval
self$summaryAllResAt <- summaryAllResAt
self$verbose <- verbose
# trueGVMatList
trueGVMatInit <- list(bsInfoInit$populations[[length(bsInfoInit$populations)]]$trueGVMat)
names(trueGVMatInit) <- bsInfoInit$populations[[length(bsInfoInit$populations)]]$name
trueGVMatInitList <- rep(list(trueGVMatInit), nIterSimulation)
names(trueGVMatInitList) <- paste0("Iteration_", 1:nIterSimulation)
self$trueGVMatInit <- trueGVMatInit[[1]]
self$trueGVMatList <- trueGVMatInitList
# estimatedGVMatList
# if (!estimatedGVInitExist) {
# self$computeLociEffInit()
# breederInfoInit$estimateGVByMLR(trainingPop = trainingPopInit,
# trainingIndNames = trainingIndNamesInit,
# testingPop = length(breederInfoInit$populationsFB),
# methodMLR = methodMLRInit,
# multiTrait = multiTraitInit,
# alpha = 0.5,
# nIter = 12000,
# burnIn = 3000,
# thin = 5,
# bayesian = TRUE)
# }
#
# estimatedGVMatInit <- list(breederInfoInit$estimatedGVByMLRInfo[[names(bsInfoInit$populations[length(bsInfoInit$populations)])]]$testingEstimatedGVByMLR)
self$computeLociEffInit()
lociEffects <- self$lociEffectsInit
genoMatNow <- bsInfoInit$populations[[length(bsInfoInit$populations)]]$genoMat
genoMatWithIntNow <- cbind(Intercept = rep(1, nrow(genoMatNow)),
genoMatNow)
estimatedGVMatInit <- list(genoMatWithIntNow[, rownames(lociEffects)] %*% lociEffects)
names(estimatedGVMatInit) <- bsInfoInit$populations[[length(bsInfoInit$populations)]]$name
estimatedGVMatInitList <- rep(list(estimatedGVMatInit), nIterSimulation)
names(estimatedGVMatInitList) <- paste0("Iteration_", 1:nIterSimulation)
self$estimatedGVMatInit <- estimatedGVMatInit[[1]]
self$estimatedGVMatList <- estimatedGVMatInitList
},
#' @description
#' estimate/extract marker/QTL effects information
computeLociEffInit = function() {
# Read arguments from `self`
bsInfoInit <- self$bsInfoInit
breederInfoInit <- self$breederInfoInit
lociEffMethod <- self$lociEffMethod
trainingPopType <- self$trainingPopType
trainingPopInit <- self$trainingPopInit
trainingIndNamesInit <- self$trainingIndNamesInit
methodMLRInit <- self$methodMLRInit
multiTraitInit <- self$multiTraitInit
verbose <- self$verbose
if (lociEffMethod == "true") {
lociEffectsInit <- bsInfoInit$lociEffects
} else if (lociEffMethod == "estimated") {
trainingPopName <- names(breederInfoInit$populationsFB)[trainingPopInit]
estimatedMrkEffName <- paste0(trainingPopName[length(trainingPopName)], "_", methodMLRInit)
if (is.null(breederInfoInit$estimatedMrkEffInfo[[estimatedMrkEffName]])) {
breederInfoInit$estimateMrkEff(trainingPop = trainingPopInit,
trainingIndNames = trainingIndNamesInit,
methodMLR = methodMLRInit,
multiTrait = multiTraitInit,
alpha = 0.5,
nIter = 12000,
burnIn = 3000,
thin = 5,
bayesian = TRUE)
}
lociEffectsInit <- breederInfoInit$lociEffects(bsInfo = bsInfoInit,
trainingPop = trainingPopInit,
methodMLR = methodMLRInit)
}
self$lociEffectsInit <- lociEffectsInit
},
#' @description
#' start simulation of breeding scheme
startSimulation = function() {
# Read arguments from `self`
simBsName <- self$simBsName
bsInfoInit <- self$bsInfoInit
breederInfoInit <- self$breederInfoInit
lociEffMethod <- self$lociEffMethod
methodMLRInit <- self$methodMLRInit
multiTraitInit <- self$multiTraitInit
nIterSimulation <- self$nIterSimulation
nGenerationProceed <- self$nGenerationProceed
nGenerationProceedSimulation <- self$nGenerationProceedSimulation
setGoalAsFinalGeneration <- self$setGoalAsFinalGeneration
performOptimization <- self$performOptimization
useFirstOptimizedValue <- self$useFirstOptimizedValue
performRobustOptimization <- self$performRobustOptimization
lowerQuantile <- self$lowerQuantile
sameAcrossGeneration <- self$sameAcrossGeneration
hMin <- self$hMin
hMax <- self$hMax
nTotalIterForOneOptimization <- self$nTotalIterForOneOptimization
nIterSimulationPerEvaluation <- self$nIterSimulationPerEvaluation
nIterSimulationForOneMrkEffect <- self$nIterSimulationForOneMrkEffect
nIterMrkEffectForRobustOptimization <- self$nIterMrkEffectForRobustOptimization
nIterOptimization <- self$nIterOptimization
lambda <- self$lambda
mu <- self$mu
nTopEvalForOpt <- self$nTopEvalForOpt
rewardWeightVec <- self$rewardWeightVec
digitsEval <- self$digitsEval
nRefreshMemoryEvery <- self$nRefreshMemoryEvery
propSelOpt <- self$propSelOpt
propSelOrNot <- self$propSelOrNot
hOriLens <- self$hOriLens
propSelLens <- self$propSelLens
hLens <- self$hLens
hStart <- self$hStart
updateBreederInfo <- self$updateBreederInfo
phenotypingInds <- self$phenotypingInds
nRepForPhenoInit <- self$nRepForPhenoInit
nRepForPheno <- self$nRepForPheno
updateModels <- self$updateModels
updateBreederInfoSimulation <- self$updateBreederInfoSimulation
phenotypingIndsSimulation <- self$phenotypingIndsSimulation
nRepForPhenoSimulation <- self$nRepForPhenoSimulation
updateModelsSimulation <- self$updateModelsSimulation
methodMLR <- self$methodMLR
multiTrait <- self$multiTrait
nSelectionWaysVec <- self$nSelectionWaysVec
selectionMethodList <- self$selectionMethodList
traitNoSelList <- self$traitNoSelList
blockSplitMethod <- self$blockSplitMethod
nMrkInBlock <- self$nMrkInBlock
minimumSegmentLength <- self$minimumSegmentLength
nSelInitOPVList <- self$nSelInitOPVList
nIterOPV <- self$nIterOPV
nProgeniesEMBVVec <- self$nProgeniesEMBVVec
nIterEMBV <- self$nIterEMBV
nCoresEMBV <- self$nCoresEMBV
clusteringForSelList <- self$clusteringForSelList
nClusterList <- self$nClusterList
nTopClusterList <- self$nTopClusterList
nTopEachList <- self$nTopEachList
nSelList <- self$nSelList
multiTraitsEvalMethodList <- self$multiTraitsEvalMethodList
hSelList <- self$hSelList
matingMethodVec <- self$matingMethodVec
allocateMethodVec <- self$allocateMethodVec
weightedAllocationMethodList <- self$weightedAllocationMethodList
traitNoRAList <- self$traitNoRAList
minimumUnitAllocateVec <- self$minimumUnitAllocateVec
includeGVPVec <- self$includeGVPVec
nNextPopVec <- self$nNextPopVec
nameMethod <- self$nameMethod
nCores <- self$nCores
overWriteRes <- self$overWriteRes
showProgress <- self$showProgress
returnMethod <- self$returnMethod
saveAllResAt <- self$saveAllResAt
evaluateGVMethod <- self$evaluateGVMethod
nTopEval <- self$nTopEval
traitNoEval <- self$traitNoEval
hEval <- self$hEval
summaryAllResAt <- self$summaryAllResAt
verbose <- self$verbose
saveObjNameBase <- self$saveObjNameBase
populationNameInit <- names(bsInfoInit$populations[length(bsInfoInit$populations)])
iterNames <- paste0("Iteration_", 1:nIterSimulation)
hVecOptsList <- self$hVecOptsList
if (!is.null(saveAllResAt)) {
saveAllResAtSplit <- stringr::str_split(string = list.files(saveAllResAt),
pattern = "_")
saveAllResAtSplitLast <- lapply(X = saveAllResAtSplit,
FUN = function(saveAllResAtSplitVec) {
return(saveAllResAtSplitVec[length(saveAllResAtSplitVec)])
})
saveAllNumeric <- unique(sort(as.numeric(stringr::str_remove(saveAllResAtSplitLast, ".rds"))))
}
if (is.null(self$lociEffectsInit)) {
self$computeLociEffInit()
}
lociEffectsInit <- self$lociEffectsInit
if (is.null(self$simBsRes[[simBsName]])) {
self$simBsRes[[simBsName]] <- list()
if ("all" %in% returnMethod) {
self$simBsRes[[simBsName]]$all <- list()
}
if ("max" %in% returnMethod) {
self$simBsRes[[simBsName]]$max <- c()
}
if ("mean" %in% returnMethod) {
self$simBsRes[[simBsName]]$mean <- c()
}
if ("median" %in% returnMethod) {
self$simBsRes[[simBsName]]$median <- c()
}
if ("min" %in% returnMethod) {
self$simBsRes[[simBsName]]$min <- c()
}
if ("var" %in% returnMethod) {
self$simBsRes[[simBsName]]$var <- c()
}
}
if (verbose) {
if (useFirstOptimizedValue) {
print("Perform optimization of hyperparameters once.")
} else {
print(paste0("Iteration: ", "1-", nIterSimulation, ", Generation: ", 1,
"; Perform optimization of hyperparameters."))
}
}
if (setGoalAsFinalGeneration) {
nGenerationProceedSimulationNow <- min(nGenerationProceed, nGenerationProceedSimulation)
} else {
nGenerationProceedSimulationNow <- nGenerationProceedSimulation
}
hVecLenNow <- sum(hLens[1:nGenerationProceedSimulationNow])
# soln <- OOR::StoSOO(par = hStart[1:hVecLenNow], fn = private$maximizeFunc,
# nGenerationProceedSimulation = nGenerationProceedSimulation,
# lower = hMin[1:hVecLenNow], upper = hMax[1:hVecLenNow],
# nb_iter = nIterOptimization,
# control = list(type = "sto", verbose = showProgress, max = TRUE))
# self$solnInit <- soln
# hVecOpt <- soln$par
if (is.null(saveObjNameBase)) {
saveObjNameBaseInit <- NULL
fileNameSolResInit <- ""
} else {
saveObjNameBaseInit <- paste0(saveObjNameBase, "_Initial")
fileNameSolResInit <- paste0(saveObjNameBaseInit, "_solution_results.rds")
}
if (!file.exists(fileNameSolResInit)) {
if (!is.null(saveObjNameBase)) {
savedObjFilesAll <- list.files(dirname(saveObjNameBase))
saveObjNameInitial <- paste0(saveObjNameBase, "_Initial")
wherePastObjFiles <- grep(pattern = basename(saveObjNameInitial),
x = savedObjFilesAll)
savedObjFilesNow <- savedObjFilesAll[wherePastObjFiles]
savedCurrentObjFile <- savedObjFilesNow[
grep(pattern = ".*_CMA-ES_object.rds",
x = savedObjFilesNow)
]
if (length(savedCurrentObjFile) == 1) {
cmaESNow <- readRDS(file = paste0(dirname(saveObjNameBase), "/",
savedCurrentObjFile))
} else {
cmaESNow <- NULL
}
} else {
cmaESNow <- NULL
}
if (is.null(cmaESNow)) {
cmaESNow <- myBreedSimulatR::cmaES$new(parameter = hStart[1:hVecLenNow],
optimizeFunc = private$maximizeFunc,
nGenerationProceedSimulation = nGenerationProceedSimulation,
lowerBound = hMin[1:hVecLenNow],
upperBound = hMax[1:hVecLenNow],
nIterOptimization = nIterOptimization,
lambda = lambda,
mu = mu,
maximize = TRUE,
saveObjNameBase = saveObjNameBaseInit,
whenToSaveObj = self$whenToSaveObj,
verbose = showProgress)
}
cmaESNow$startOptimization()
saveRDS(object = cmaESNow,
file = cmaESNow$fileNameSaveObj)
optimalHyperParamMat <- cmaESNow$optimalHyperParamMat
hVecOpt <- cmaESNow$optimalParameter
solnInit <- list(value = cmaESNow$optimalValue,
par = hVecOpt)
rm(cmaESNow)
gc(reset = TRUE); gc(reset = TRUE)
solnRes <- c(solnInit,
list(optimalHyperParamMat = optimalHyperParamMat))
if (!is.null(saveObjNameBase)) {
saveRDS(object = solnRes, file = fileNameSolResInit)
}
} else {
solnRes <- readRDS(file = fileNameSolResInit)
solnInit <- solnRes[1:2]
hVecOpt <- solnRes$par
optimalHyperParamMat <- solnRes$optimalHyperParamMat
}
self$solnInit <- solnInit
self$optimalHyperParamMatsList[["Initial"]] <- optimalHyperParamMat
hVecOptsList[["Initial"]] <- hVecOpt
if (sameAcrossGeneration) {
hListOpt <- sapply(X = hLens[1:nGenerationProceedSimulationNow],
FUN = function(hLen) {
hVecOpt[1:hLen]
}, simplify = FALSE)
} else {
hListOpt <- split(x = hVecOpt[1:hVecLenNow],
f = rep(1:nGenerationProceedSimulationNow,
hLens[1:nGenerationProceedSimulationNow]))
}
if (useFirstOptimizedValue) {
# save
if (verbose) {
print("Perform simulation based on optimized hyperparameters.")
}
simBsCma <- myBreedSimulatR::simBs$new(simBsName = simBsName,
bsInfoInit = bsInfoInit,
breederInfoInit = breederInfoInit,
lociEffMethod = lociEffMethod,
trainingPopType = trainingPopType,
trainingPopInit = trainingPopInit,
trainingIndNamesInit = trainingIndNamesInit,
methodMLRInit = methodMLRInit,
multiTraitInit = multiTraitInit,
samplingMrkEffInit = FALSE,
seedMrkEffSamplingInit = NA,
nIterSimulation = nIterSimulation,
nGenerationProceed = nGenerationProceed,
nRefreshMemoryEvery = nRefreshMemoryEvery,
updateBreederInfo = updateBreederInfo,
phenotypingInds = phenotypingInds,
nRepForPhenoInit = nRepForPhenoInit,
nRepForPheno = nRepForPheno,
updateModels = updateModels,
methodMLR = methodMLR,
multiTrait = multiTrait,
nSelectionWaysVec = nSelectionWaysVec,
selectionMethodList = selectionMethodList,
traitNoSelList = traitNoSelList,
blockSplitMethod = blockSplitMethod,
nMrkInBlock = nMrkInBlock,
minimumSegmentLength = minimumSegmentLength,
nSelInitOPVList = nSelInitOPVList,
nIterOPV = nIterOPV,
nProgeniesEMBVVec = nProgeniesEMBVVec,
nIterEMBV = nIterEMBV,
nCoresEMBV = nCoresEMBV,
clusteringForSelList = clusteringForSelList,
nClusterList = nClusterList,
nTopClusterList = nTopClusterList,
nTopEachList = nTopEachList,
nSelList = nSelList,
multiTraitsEvalMethodList = multiTraitsEvalMethodList,
hSelList = hSelList,
matingMethodVec = matingMethodVec,
allocateMethodVec = allocateMethodVec,
weightedAllocationMethodList = weightedAllocationMethodList,
includeGVPVec = includeGVPVec,
targetGenGVPVec = targetGenGVPVec,
performOCSVec = performOCSVec,
targetGenOCSVec = targetGenOCSVec,
HeStarRatioVec = HeStarRatioVec,
degreeOCSVec = degreeOCSVec,
includeGVPOCSVec = includeGVPOCSVec,
hOCSList = hOCSList,
weightedAllocationMethodOCSList = weightedAllocationMethodOCSList,
traitNoOCSList = traitNoOCSList,
nCrossesOCSVec = nCrossesOCSVec,
traitNoRAList = traitNoRAList,
hList = hListOpt,
minimumUnitAllocateVec = minimumUnitAllocateVec,
nNextPopVec = nNextPopVec,
nameMethod = nameMethod,
nCores = nCoresPerOptimization,
overWriteRes = overWriteRes,
showProgress = showProgress,
returnMethod = returnMethod,
saveAllResAt = saveAllResAt,
evaluateGVMethod = evaluateGVMethod,
nTopEval = nTopEval,
traitNoEval = traitNoEval,
hEval = hEval,
summaryAllResAt = summaryAllResAt,
verbose = verbose)
simBsCma$startSimulation()
self$simBsRes[[simBsName]] <- simBsCma$simBsRes[[simBsName]]
self$trueGVMatList <- simBsCma$trueGVMatList
self$estimatedGVMatList <- simBsCma$estimatedGVMatList
} else {
if (nCores == 1) {
conductSimulations <- sapply(X = iterNames,
FUN = function(iterName) {
if (!((is.null(self$simBsRes[[simBsName]]$all[[iterName]])) &
(length(self$trueGVMatList[[iterName]]) <= 1))) {
if (overWriteRes) {
conductSimulation <- TRUE
} else {
conductSimulation <- FALSE
}
} else {
if ("all" %in% returnMethod) {
conductSimulation <- TRUE
} else {
if (any(sapply(returnMethod, function(x) is.null(self$simBsRes[[simBsName]][[x]][iterName])))) {
conductSimulation <- TRUE
} else {
if (overWriteRes) {
conductSimulation <- TRUE
} else {
conductSimulation <- FALSE
}
}
}
}
return(conductSimulation)
})
if (!is.null(saveAllResAt)) {
if (!overWriteRes) {
conductSimulations[saveAllNumeric] <- FALSE
}
}
simulationCounts <- 0
# iterNo <- 1
for (iterNo in 1:nIterSimulation) {
iterName <- iterNames[iterNo]
if (conductSimulations[iterNo]) {
simulationCounts <- simulationCounts + 1
bsInfo <- bsInfoInit$clone(deep = FALSE)
breederInfo <- breederInfoInit$clone(deep = FALSE)
lociEffects <- lociEffectsInit
hCount <- 0
# trueGVMatList
if (is.null(self$trueGVMatList[[iterName]])) {
self$trueGVMatList[[iterName]][[populationNameInit]] <- self$trueGVMatInit
}
# estimatedGVMatList
if (is.null(self$estimatedGVMatList[[iterName]])) {
self$estimatedGVMatList[[iterName]][[populationNameInit]] <- self$estimatedGVMatInit
}
simRes <- private$performOneSimulationOpt(iterNo = iterNo)
# simRes <- private$performOneSimulationOptTryError(iterNo = iterNo)
if ("all" %in% returnMethod) {
self$simBsRes[[simBsName]]$all[[iterName]] <- simRes$all
}
if ("max" %in% returnMethod) {
self$simBsRes[[simBsName]]$max[iterName] <- simRes$max
}
if ("mean" %in% returnMethod) {
self$simBsRes[[simBsName]]$mean[iterName] <- simRes$mean
}
if ("median" %in% returnMethod) {
self$simBsRes[[simBsName]]$median[iterName] <- simRes$median
}
if ("min" %in% returnMethod) {
self$simBsRes[[simBsName]]$min[iterName] <- simRes$min
}
if ("var" %in% returnMethod) {
self$simBsRes[[simBsName]]$var[iterName] <- simRes$var
}
self$trueGVMatList[[iterName]] <- c(self$trueGVMatList[[iterName]],
simRes$trueGVMatList)
self$estimatedGVMatList[[iterName]] <- c(self$estimatedGVMatList[[iterName]],
simRes$estimatedGVMatList)
if (simulationCounts %% nRefreshMemoryEvery == 0) {
rm(simRes); rm(bsInfo); rm(breederInfo)
gc(reset = TRUE); gc(reset = TRUE)
}
}
}
if (showProgress) {
cat("\n")
}
} else {
conductSimulations <- sapply(X = iterNames,
FUN = function(iterName) {
if (!((is.null(self$simBsRes[[simBsName]]$all[[iterName]])) &
(length(self$trueGVMatList[[iterName]]) <= 1))) {
if (overWriteRes) {
conductSimulation <- TRUE
} else {
conductSimulation <- FALSE
}
} else {
if ("all" %in% returnMethod) {
conductSimulation <- TRUE
} else {
if (any(sapply(returnMethod, function(x) is.null(self$simBsRes[[simBsName]][[x]][iterName])))) {
conductSimulation <- TRUE
} else {
if (overWriteRes) {
conductSimulation <- TRUE
} else {
conductSimulation <- FALSE
}
}
}
}
return(conductSimulation)
})
if (!is.null(saveAllResAt)) {
if (!overWriteRes) {
conductSimulations[saveAllNumeric] <- FALSE
}
}
if (any(conductSimulations)) {
if (showProgress) {
simResAll <- pbmcapply::pbmclapply(X = (1:nIterSimulation)[conductSimulations],
FUN = private$performOneSimulationOptTryError,
mc.cores = nCores)
} else {
simResAll <- parallel::mclapply(X = (1:nIterSimulation)[conductSimulations],
FUN = private$performOneSimulationOptTryError,
mc.cores = nCores)
}
names(simResAll) <- iterNames[conductSimulations]
if ("all" %in% returnMethod) {
self$simBsRes[[simBsName]]$all[iterNames[conductSimulations]] <- lapply(simResAll, function(x) x$all)
}
if ("max" %in% returnMethod) {
self$simBsRes[[simBsName]]$max[iterNames[conductSimulations]] <- unlist(lapply(simResAll, function(x) x$max))
}
if ("mean" %in% returnMethod) {
self$simBsRes[[simBsName]]$mean[iterNames[conductSimulations]] <- unlist(lapply(simResAll, function(x) x$mean))
}
if ("median" %in% returnMethod) {
self$simBsRes[[simBsName]]$median[iterNames[conductSimulations]] <- unlist(lapply(simResAll, function(x) x$median))
}
if ("min" %in% returnMethod) {
self$simBsRes[[simBsName]]$min[iterNames[conductSimulations]] <- unlist(lapply(simResAll, function(x) x$min))
}
if ("var" %in% returnMethod) {
self$simBsRes[[simBsName]]$var[iterNames[conductSimulations]] <- unlist(lapply(simResAll, function(x) x$var))
}
trueGVMatListNow <- lapply(simResAll, function(x) x$trueGVMatList)
self$trueGVMatList[iterNames[conductSimulations]] <- sapply(iterNames[conductSimulations],
function(iterName) {
c(self$trueGVMatList[[iterName]],
trueGVMatListNow[[iterName]])
}, simplify = FALSE)
estimatedGVMatListNow <- lapply(simResAll, function(x) x$estimatedGVMatList)
self$estimatedGVMatList[iterNames[conductSimulations]] <- sapply(iterNames[conductSimulations],
function(iterName) {
c(self$estimatedGVMatList[[iterName]],
estimatedGVMatListNow[[iterName]])
}, simplify = FALSE)
}
}
}
},
#' @description
#' start simulation of breeding scheme
summaryResults = function() {
# Read arguments from `self`
simBsName <- self$simBsName
bsInfoInit <- self$bsInfoInit
showProgress <- self$showProgress
evaluateGVMethod <- self$evaluateGVMethod
nTopEval <- self$nTopEval
simBsRes <- self$simBsRes
summaryAllResAt <- self$summaryAllResAt
nIterSimulation <- self$nIterSimulation
iterNames <- paste0("Iteration_", 1:nIterSimulation)
if (!is.null(summaryAllResAt)) {
fileNameBsInfoRes0 <- here::here(summaryAllResAt,
paste0(simBsName, "_bsInfo_"))
fileNameBreederInfoRes0 <- here::here(summaryAllResAt,
paste0(simBsName, "_breederInfo_"))
if (showProgress) {
listOfGVMatList <- pbmcapply::pbmclapply(X = iterNames,
FUN = private$extractGVMatList,
mc.cores = self$nCores)
} else {
listOfGVMatList <- parallel::mclapply(X = iterNames,
FUN = private$extractGVMatList,
mc.cores = self$nCores)
}
if (!is.null(listOfGVMatList$warning)) {
listOfGVMatList <- listOfGVMatList$value
}
trueGVMatList <- lapply(listOfGVMatList, function(x) x$trueGVMatEachList)
names(trueGVMatList) <- iterNames
trueGVMatListNonNULL <- which(!unlist(lapply(trueGVMatList, is.null)))
trueGVMatList <- trueGVMatList[trueGVMatListNonNULL]
self$trueGVMatList <- trueGVMatList
estimatedGVMatList <- lapply(listOfGVMatList, function(x) x$estimatedGVMatEachList)
names(estimatedGVMatList) <- iterNames
estimatedGVMatListNonNULL <- which(!unlist(lapply(estimatedGVMatList, is.null)))
estimatedGVMatList <- estimatedGVMatList[estimatedGVMatListNonNULL]
self$estimatedGVMatList <- estimatedGVMatList
} else {
trueGVMatList <- self$trueGVMatList
estimatedGVMatList <- self$estimatedGVMatList
if ((is.null(self$simBsRes[[simBsName]]$all)) &
(length(self$trueGVMatList[[1]]) == 1)) {
stop("Please start simulation with `returnMethod = 'summary'`. You do not have simulation results.")
}
}
trueGVSummaryArrayList <- lapply(X = trueGVMatList,
FUN = private$extractTrueSummaryRes)
trueGVSummaryArray <- do.call(what = abind::abind,
args = trueGVSummaryArrayList)
dimnames(trueGVSummaryArray)[c(1, 3, 4)] <- list(Index = c("max", "mean", "median", "min", "var"),
Population = names(trueGVMatList[[1]]),
Iteration = names(trueGVMatList))
self$trueGVSummaryArray <- trueGVSummaryArray
estimatedGVSummaryArrayList <- lapply(X = estimatedGVMatList,
FUN = private$extractEstimatedSummaryRes)
estimatedGVSummaryArray <- do.call(what = abind::abind,
args = estimatedGVSummaryArrayList)
dimnames(estimatedGVSummaryArray)[c(1, 3, 4)] <- list(Index = c("max", "mean", "median", "min", "var"),
Population = names(estimatedGVMatList[[1]]),
Iteration = names(estimatedGVMatList))
self$estimatedGVSummaryArray <- estimatedGVSummaryArray
},
#' @description
#' Display information about the object
print = function() {
cat(paste0(
"Name of Simulation Settings: ", self$simBsName, "\n",
"Number of Simulations: ", self$nIterSimulation, "\n",
"Number of Generations to be Proceeded: ", self$nGenerationProceed, "\n",
"Number of Individuals for Next Generation: \n"
))
print(self$nNextPopVec)
},
#' @description Draw figures for visualization of simulation results for summary statistics
#' @param targetTrait [numeric] Target trait. character is OK, but numeric vector
#' corresponding to target traits is preferred. It should be a vector with length 1.
#' @param targetPopulation [numeric] Target populations. character is OK, but numeric vector
#' corresponding to target traits is preferred.
#' @param plotType [character] We offer "box", "violin", "lines", "density"
#' to draw figures for simulation results.
#' @param plotTargetDensity [character] When you choose "density" for `plotType`,
#' you should select which summary statistics will be plotted. It should be a vector with length 1.
#' @param plotGVMethod [character] Which type of GV (true GV or estimated GV) will be used for plotting the simulation results
#' @param adjust [numeric] the bandwidth used is actually adjust*bw. This makes it easy to specify values like ‘half the default’ bandwidth.
#' (see: `adjust` in \link[stats]{density})
#'
plot = function(targetTrait = 1,
targetPopulation = NULL,
plotType = "box",
plotTargetDensity = "max",
plotGVMethod = "true",
adjust = 1e-05) {
lociEffMethodsOffered <- c("true", "estimated")
# targetTrait
if (is.numeric(targetTrait)) {
targetTraitName <- self$bsInfoInit$traitInfo$traitNames[targetTrait]
} else if (is.character(targetTrait)) {
targetTraitName <- targetTrait
} else {
stop("`targetTraitName` must be `numeric` or `character`!")
}
stopifnot(length(targetTraitName) == 1)
# plotType
plotTypeOffered <- c("box", "violin", "lines", "density")
stopifnot(length(plotType) == 1)
stopifnot(plotType %in% plotTypeOffered)
# plotGVMethod
if (!is.null(plotGVMethod)) {
if (!(plotGVMethod %in% lociEffMethodsOffered)) {
stop(paste0("We only offer the following methods for evaluating the simulation results: ",
paste(lociEffMethodsOffered, collapse = "; ")))
}
} else {
plotGVMethod <- "true"
}
if (plotGVMethod == "true") {
if (is.null(self$trueGVSummaryArray)) {
self$summaryResults()
}
trueGVSummaryArray <- self$trueGVSummaryArray
} else {
if (is.null(self$estimatedGVSummaryArray)) {
self$summaryResults()
}
trueGVSummaryArray <- self$estimatedGVSummaryArray
}
# targetPopulation
if (is.null(targetPopulation)) {
targetPopulation <- 1:dim(trueGVSummaryArray)[3]
}
targetPopulation <- targetPopulation[targetPopulation %in% (1:dim(trueGVSummaryArray)[3])]
trueGVSummaryArray <- trueGVSummaryArray[ , , targetPopulation, , drop = FALSE]
dimSummary <- dim(trueGVSummaryArray)
dimnamesSummary <- dimnames(trueGVSummaryArray)
trueGVSummaryDf <- data.frame(SummaryStatistics = rep(dimnamesSummary[[1]], prod(dimSummary[-1])),
Trait = rep(rep(dimnamesSummary[[2]], each = dimSummary[1]),
prod(dimSummary[3:4])),
Population = rep(rep(dimnamesSummary[[3]], each = prod(dimSummary[1:2])),
prod(dimSummary[4])),
Iteration = rep(dimnamesSummary[[4]], each = prod(dimSummary[1:3])),
Value = c(trueGVSummaryArray))
trueGVSummaryDf$SummaryStatistics <- factor(trueGVSummaryDf$SummaryStatistics, levels = dimnamesSummary[[1]])
trueGVSummaryDf$Trait <- factor(trueGVSummaryDf$Trait, levels = dimnamesSummary[[2]])
trueGVSummaryDf$Population <- factor(trueGVSummaryDf$Population, levels = dimnamesSummary[[3]])
trueGVSummaryDf$Iteration <- factor(trueGVSummaryDf$Iteration, levels = dimnamesSummary[[4]])
if (plotType %in% c("box", "violin")) {
# trueGVSummaryDfTarget <- trueGVSummaryDf[trueGVSummaryDf$SummaryStatistics %in% plotTarget, ]
trueGVSummaryDfTarget <- trueGVSummaryDf[trueGVSummaryDf$Trait %in% targetTraitName, ]
trueGVSummaryDfTarget$Value <- round(trueGVSummaryDfTarget$Value, 3)
plt <- plotly::plot_ly(
data = trueGVSummaryDfTarget,
x = ~ Population,
y = ~ Value,
split = ~ SummaryStatistics,
type = plotType,
hoverinfo = "text",
# boxpoints = "all",
# jitter = 0.3,
# pointpos = -1.8,
text = paste0(apply(trueGVSummaryDfTarget, 1, function(l) {
paste(names(l), ":", l, collapse = "\n")
}))
) %>%
plotly::layout(title = list(text = targetTraitName),
# yaxis = list(title = list(text = plotTarget))
yaxis = list(title = list(text = paste0(plotGVMethod, " GV")))
)
if (plotType == "box") {
plt <- plt %>% plotly::layout(boxmode = "group")
} else if (plotType == "violin") {
plt <- plt %>% plotly::layout(violinmode = "group")
}
} else if (plotType %in% c("lines")) {
trueGVSummaryMeanArray <- apply(X = trueGVSummaryArray,
MARGIN = 1:3, FUN = mean)
dimSummaryMean <- dim(trueGVSummaryMeanArray)
dimnamesSummaryMean <- dimnames(trueGVSummaryMeanArray)
trueGVSummaryMeanDf <- data.frame(SummaryStatistics = rep(dimnamesSummaryMean[[1]], prod(dimSummaryMean[-1])),
Trait = rep(rep(dimnamesSummaryMean[[2]], each = dimSummaryMean[1]),
prod(dimSummaryMean[3])),
Population = rep(dimnamesSummaryMean[[3]], each = prod(dimSummaryMean[1:2])),
Value = c(trueGVSummaryMeanArray))
trueGVSummaryMeanDf$SummaryStatistics <- factor(trueGVSummaryMeanDf$SummaryStatistics, levels = dimnamesSummaryMean[[1]])
trueGVSummaryMeanDf$Trait <- factor(trueGVSummaryMeanDf$Trait, levels = dimnamesSummaryMean[[2]])
trueGVSummaryMeanDf$Population <- factor(trueGVSummaryMeanDf$Population, levels = dimnamesSummaryMean[[3]])
trueGVSummaryMeanDfTarget <- trueGVSummaryMeanDf[trueGVSummaryMeanDf$Trait %in% targetTraitName, ]
trueGVSummaryMeanDfTarget$Value <- round(trueGVSummaryMeanDfTarget$Value, 3)
plt <- plot_ly(
data = trueGVSummaryMeanDfTarget,
x = ~ Population,
y = ~ Value,
split = ~ SummaryStatistics,
# line = list(color = colorVec[1],
# width = widthVec[1],
# dash = dashVec[1]),
type = "scatter",
mode = "markers+lines",
hoverinfo = "text",
text = paste0(apply(trueGVSummaryMeanDfTarget, 1, function(l) {
paste(names(l), ":", l, collapse = "\n")
}))
) %>%
plotly::layout(title = list(text = targetTraitName),
yaxis = list(title = list(text = paste0(plotGVMethod, " GV"))))
} else if (plotType == "density") {
trueGVSummaryDfTarget <- trueGVSummaryDf[trueGVSummaryDf$Trait %in% targetTraitName, ]
trueGVSummaryDfTargetSS <- trueGVSummaryDfTarget[trueGVSummaryDfTarget$SummaryStatistics %in% plotTargetDensity, ]
densityValueDfList <- lapply(X = dimnames(trueGVSummaryArray)[[3]],
FUN = function(popName) {
trueGVSummaryDfTargetSSEachPop <- trueGVSummaryDfTargetSS[trueGVSummaryDfTargetSS$Population %in% popName, ]
densityResEachPop <- density(x = sort(trueGVSummaryDfTargetSSEachPop$Value), adjust = adjust)
x <- densityResEachPop$x
x <- c(min(x), x)
y <- cumsum(densityResEachPop$y / sum(densityResEachPop$y))
y <- c(0, y)
return(data.frame(x, y))
})
densityValueDf <- do.call(what = rbind,
args = densityValueDfList)
densityValueDf$Population <- rep(dimnames(trueGVSummaryArray)[[3]], unlist(lapply(densityValueDfList, nrow)))
densityValueDf$Population <- factor(densityValueDf$Population, levels = dimnames(trueGVSummaryArray)[[3]])
plt <- plot_ly(
data = densityValueDf,
x = ~ x,
y = ~ y,
split = ~ Population,
# line = list(color = colorVec[1],
# width = widthVec[1],
# dash = dashVec[1]),
type = "scatter",
mode = "lines"
# name = nameVec[1]
) %>%
plotly::layout(title = list(text = paste0(targetTraitName, "-", plotTargetDensity)),
xaxis = list(title = list(text = paste0(plotGVMethod, " GV"))))
}
return(plt)
}
),
private = list(
# @description marker and QTL effects used for crossInfo object
#
# @param ind1 [individual class] parent 1
# @param ind2 [individual class] parent 2
# @param names [character] names of the descendants
# @param n [numeric] number of descendants
lociEffects = function(bsInfo,
breederInfo,
alpha = 0.5,
nIter = 5000,
burnIn = 1000,
thin = 5,
bayesian = FALSE,
gcReset = FALSE) {
# Read arguments from `self`
lociEffMethod <- self$lociEffMethod
trainingPopType <- self$trainingPopType
methodMLR <- self$methodMLR
multiTrait <- self$multiTrait
verbose <- self$verbose
if (lociEffMethod == "true") {
lociEffects <- bsInfo$lociEffects
} else if (lociEffMethod == "estimated") {
trainingPopName <- names(breederInfo$populationsFB)
infoName <- paste0(trainingPopName[length(trainingPopName)], "_", methodMLR)
if (trainingPopType == "all") {
trainingPop <- 1:length(breederInfo$populationsFB)
} else {
trainingPop <- length(breederInfo$populationsFB)
}
if (is.null(breederInfo$estimatedMrkEffInfo[[infoName]])) {
breederInfo$estimateMrkEff(trainingPop = trainingPop,
methodMLR = methodMLR,
multiTrait = multiTrait,
alpha = alpha,
nIter = nIter,
burnIn = burnIn,
thin = thin,
bayesian = bayesian)
}
lociEffects <- breederInfo$lociEffects(bsInfo = bsInfo,
trainingPop = trainingPop,
methodMLR = methodMLR)
if (gcReset) {
rm(bsInfo); rm(breederInfo)
gc(reset = TRUE); gc(reset = TRUE)
}
}
return(lociEffects)
},
# @description marker and QTL effects used for crossInfo object
# @param hVec [numeric] hyperparameter to be optimized
maximizeFunc = function(hVec, nGenerationProceedSimulation,
bsInfo = NULL, breederInfo = NULL) {
performRobustOptimization <- self$performRobustOptimization
if (self$sameAcrossGeneration) {
if (self$propSelOpt) {
propSelList <- sapply(X = self$propSelLens[1:nGenerationProceedSimulation],
FUN = function(propSelLen) {
hVec[self$propSelOrNot][1:propSelLen]
}, simplify = FALSE)
nSelList <- mapply(FUN = function(propSel, nInd) {
round(nInd * propSel)
},
propSelList,
c(self$bsInfoInit$populations[[length(bsInfoInit$populations)]]$nInd,
self$nNextPopVec),
SIMPLIFY = FALSE)
hList <- sapply(X = self$hLensOri[1:nGenerationProceedSimulation],
FUN = function(hLen) {
hVec[!self$propSelOrNot][1:hLen]
}, simplify = FALSE)
} else {
propSelList <- NULL
nSelList <- self$nSelList
hList <- sapply(X = self$hLens[1:nGenerationProceedSimulation],
FUN = function(hLen) {
hVec[1:hLen]
}, simplify = FALSE)
}
} else {
hVecLenNow <- sum(self$hLens[1:nGenerationProceedSimulation])
if (self$propSelOpt) {
propSelList <- split(x = hVec[1:hVecLenNow][self$propSelOrNot],
f = rep(1:nGenerationProceedSimulation,
self$propSelLens[1:nGenerationProceedSimulation]))
nSelList <- mapply(FUN = function(propSel, nInd) {
round(nInd * propSel)
},
propSelList,
c(self$bsInfoInit$populations[[length(bsInfoInit$populations)]]$nInd,
self$nNextPopVec),
SIMPLIFY = FALSE)
hList <- split(x = hVec[1:hVecLenNow][!self$propSelOrNot],
f = rep(1:nGenerationProceedSimulation,
self$hLensOri[1:nGenerationProceedSimulation]))
} else {
propSelList <- NULL
nSelList <- self$nSelList
hList <- split(x = hVec[1:hVecLenNow],
f = rep(1:nGenerationProceedSimulation,
self$hLens[1:nGenerationProceedSimulation]))
}
}
rewardWeightVec <- self$rewardWeightVec[(self$nGenerationProceedSimulation - nGenerationProceedSimulation + 1):self$nGenerationProceedSimulation]
rewardWeightVec <- rewardWeightVec / sum(rewardWeightVec)
nonZeroWeight <- rewardWeightVec != 0
rewardOnlyLast <- (!any(nonZeroWeight[1:(nGenerationProceedSimulation - 1)])) & nonZeroWeight[nGenerationProceedSimulation]
if (rewardOnlyLast) {
returnMethod <- "max"
} else {
returnMethod <- "summary"
}
if (performRobustOptimization) {
samplingMrkEffInit <- TRUE
nCoresForOneMrkEff <- 1
} else {
samplingMrkEffInit <- FALSE
nCoresForOneMrkEff <- self$nCoresPerOptimization
}
if (is.null(bsInfo)) {
bsInfo <- self$bsInfoInit$clone(deep = FALSE)
}
if (is.null(breederInfo)) {
breederInfo <- self$breederInfoInit$clone(deep = FALSE)
}
maximizeFuncForOneMrkEffect <- function(iterNoForMrkEffects) {
simBsNow <- myBreedSimulatR::simBs$new(simBsName = self$simBsName,
bsInfoInit = bsInfo,
breederInfoInit = breederInfo,
lociEffMethod = self$lociEffMethod,
trainingPopType = self$trainingPopType,
trainingPopInit = self$trainingPopInit,
trainingIndNamesInit = self$trainingIndNamesInit,
methodMLRInit = self$methodMLRInit,
multiTraitInit = self$multiTraitInit,
samplingMrkEffInit = samplingMrkEffInit,
seedMrkEffSamplingInit = NA,
nIterSimulation = self$nIterSimulationForOneMrkEffect,
nGenerationProceed = nGenerationProceedSimulation,
nRefreshMemoryEvery = self$nRefreshMemoryEvery,
updateBreederInfo = self$updateBreederInfoSimulation[1:nGenerationProceedSimulation],
phenotypingInds = self$phenotypingIndsSimulation[1:nGenerationProceedSimulation],
nRepForPhenoInit = self$nRepForPhenoInit,
nRepForPheno = self$nRepForPhenoSimulation[1:nGenerationProceedSimulation],
updateModels = self$updateModelsSimulation[1:nGenerationProceedSimulation],
methodMLR = self$methodMLR,
multiTrait = self$multiTrait,
nSelectionWaysVec = self$nSelectionWaysVec[1:nGenerationProceedSimulation],
selectionMethodList = self$selectionMethodList[1:nGenerationProceedSimulation],
traitNoSelList = self$traitNoSelList[1:nGenerationProceedSimulation],
blockSplitMethod = self$blockSplitMethod,
nMrkInBlock = self$nMrkInBlock,
minimumSegmentLength = self$minimumSegmentLength,
nSelInitOPVList = self$nSelInitOPVList[1:nGenerationProceedSimulation],
nIterOPV = self$nIterOPV,
nProgeniesEMBVVec = self$nProgeniesEMBVVec[1:nGenerationProceedSimulation],
nIterEMBV = self$nIterEMBV,
nCoresEMBV = self$nCoresEMBV,
clusteringForSelList = self$clusteringForSelList[1:nGenerationProceedSimulation],
nClusterList = self$nClusterList[1:nGenerationProceedSimulation],
nTopClusterList = self$nTopClusterList[1:nGenerationProceedSimulation],
nTopEachList = self$nTopEachList[1:nGenerationProceedSimulation],
nSelList = nSelList[1:nGenerationProceedSimulation],
multiTraitsEvalMethodList = self$multiTraitsEvalMethodList[1:nGenerationProceedSimulation],
hSelList = self$hSelList[1:nGenerationProceedSimulation],
matingMethodVec = self$matingMethodVec[1:nGenerationProceedSimulation],
allocateMethodVec = self$allocateMethodVec[1:nGenerationProceedSimulation],
weightedAllocationMethodList = self$weightedAllocationMethodList[1:nGenerationProceedSimulation],
includeGVPVec = self$includeGVPVec[1:nGenerationProceedSimulation],
targetGenGVPVec = targetGenGVPVec[1:nGenerationProceedSimulation],
performOCSVec = performOCSVec[1:nGenerationProceedSimulation],
targetGenOCSVec = targetGenOCSVec[1:nGenerationProceedSimulation],
HeStarRatioVec = HeStarRatioVec[1:nGenerationProceedSimulation],
degreeOCSVec = degreeOCSVec[1:nGenerationProceedSimulation],
includeGVPOCSVec = includeGVPOCSVec[1:nGenerationProceedSimulation],
hOCSList = hOCSList[1:nGenerationProceedSimulation],
weightedAllocationMethodOCSList = weightedAllocationMethodOCSList[1:nGenerationProceedSimulation],
traitNoOCSList = traitNoOCSList[1:nGenerationProceedSimulation],
nCrossesOCSVec = nCrossesOCSVec[1:nGenerationProceedSimulation],
traitNoRAList = self$traitNoRAList[1:nGenerationProceedSimulation],
hList = hList[1:nGenerationProceedSimulation],
minimumUnitAllocateVec = self$minimumUnitAllocateVec[1:nGenerationProceedSimulation],
nNextPopVec = self$nNextPopVec[1:nGenerationProceedSimulation],
nameMethod = self$nameMethod,
nCores = nCoresForOneMrkEff,
overWriteRes = self$overWriteRes,
showProgress = FALSE,
returnMethod = returnMethod,
saveAllResAt = NULL,
evaluateGVMethod = "estimated",
nTopEval = self$nTopEvalForOpt,
traitNoEval = self$traitNoEval,
hEval = self$hEval,
summaryAllResAt = NULL,
verbose = FALSE)
simBsNow$startSimulation()
if (rewardOnlyLast) {
maxEvalForOneMrkEffect <- mean(simBsNow$simBsRes[[simBsNow$simBsName]]$max, na.rm = TRUE)
} else {
rewardEachIter <- lapply(X = simBsNow$estimatedGVMatList,
FUN = function(estimatedGVMatEachIter) {
estimatedGVMatMax <- try(do.call(what = rbind,
args = lapply(X = estimatedGVMatEachIter,
FUN = function(estimatedGVMatEachPop) {
apply(estimatedGVMatEachPop, 2, max)
})
), silent = TRUE)
if ("try-error" %in% class(estimatedGVMatMax)) {
reward <- NA
} else {
estimatedGVMaxEachPop <- try(estimatedGVMatMax[, self$traitNoEval, drop = FALSE] %*% as.matrix(self$hEval),
silent = TRUE)
if ("try-error" %in% class(estimatedGVMaxEachPop)) {
reward <- NA
} else {
reward <- try(sum(estimatedGVMaxEachPop[-1, 1] * rewardWeightVec),
silent = TRUE)
if ("try-error" %in% class(reward)) {
reward <- NA
}
}
}
return(reward)
})
maxEvalForOneMrkEffect <- mean(unlist(rewardEachIter), na.rm = TRUE)
}
return(maxEvalForOneMrkEffect)
}
if (performRobustOptimization) {
maximizeFuncForOneMrkEffectTryError <- function(iterNoForMrkEffects) {
performMaximizeFuncForOneMrkEffect <- TRUE
while (performMaximizeFuncForOneMrkEffect) {
maxEvalForOneMrkEffect <- try(maximizeFuncForOneMrkEffect(iterNoForMrkEffects = iterNoForMrkEffects),
silent = TRUE)
performMaximizeFuncForOneMrkEffect <- "try-error" %in% class(maxEvalForOneMrkEffect)
}
return(maxEvalForOneMrkEffect)
}
if (self$nCoresPerOptimization == 1) {
if (self$showProgress) {
maxEvalsForMrkEffects <- unlist(
pbapply::pblapply(X = 1:self$nIterMrkEffectForRobustOptimization,
FUN = maximizeFuncForOneMrkEffectTryError)
)
} else {
maxEvalsForMrkEffects <- unlist(
lapply(X = 1:self$nIterMrkEffectForRobustOptimization,
FUN = maximizeFuncForOneMrkEffectTryError)
)
}
} else {
if (self$showProgress) {
maxEvalsForMrkEffects <- unlist(
pbmcapply::pbmclapply(X = 1:self$nIterMrkEffectForRobustOptimization,
FUN = maximizeFuncForOneMrkEffectTryError,
mc.cores = self$nCoresPerOptimization)
)
} else {
maxEvalsForMrkEffects <- unlist(
parallel::mclapply(X = 1:self$nIterMrkEffectForRobustOptimization,
FUN = maximizeFuncForOneMrkEffectTryError,
mc.cores = self$nCoresPerOptimization)
)
}
}
maxEval <- quantile(x = maxEvalsForMrkEffects, probs = self$lowerQuantile)
} else {
maxEval <- maximizeFuncForOneMrkEffect()
}
if (!is.na(maxEval)) {
maxEval <- round(maxEval,
digits = self$digitsEval)
} else {
maxEval <- -Inf
message("Computation failed: maybe internal error occured !")
}
return(maxEval)
},
# @description Proceed optimized breeding scheme (one simulation)
#
# @param iterNo [numeric] Iteration No.
performOneSimulationOpt = function(iterNo) {
# Read arguments from `self`
simBsName <- self$simBsName
bsInfoInit <- self$bsInfoInit
breederInfoInit <- self$breederInfoInit
lociEffMethod <- self$lociEffMethod
methodMLRInit <- self$methodMLRInit
multiTraitInit <- self$multiTraitInit
nIterSimulation <- self$nIterSimulation
nGenerationProceed <- self$nGenerationProceed
nGenerationProceedSimulation <- self$nGenerationProceedSimulation
setGoalAsFinalGeneration <- self$setGoalAsFinalGeneration
performOptimization <- self$performOptimization
useFirstOptimizedValue <- self$useFirstOptimizedValue
sameAcrossGeneration <- self$sameAcrossGeneration
hMin <- self$hMin
hMax <- self$hMax
nTotalIterForOneOptimization <- self$nTotalIterForOneOptimization
nIterSimulationPerEvaluation <- self$nIterSimulationPerEvaluation
nIterOptimization <- self$nIterOptimization
lambda <- self$lambda
mu <- self$mu
nTopEvalForOpt <- self$nTopEvalForOpt
digitsEval <- self$digitsEval
nRefreshMemoryEvery <- self$nRefreshMemoryEvery
hLens <- self$hLens
hStart <- self$hStart
updateBreederInfo <- self$updateBreederInfo
phenotypingInds <- self$phenotypingInds
nRepForPhenoInit <- self$nRepForPhenoInit
nRepForPheno <- self$nRepForPheno
updateModels <- self$updateModels
updateBreederInfoSimulation <- self$updateBreederInfoSimulation
phenotypingIndsSimulation <- self$phenotypingIndsSimulation
nRepForPhenoSimulation <- self$nRepForPhenoSimulation
updateModelsSimulation <- self$updateModelsSimulation
methodMLR <- self$methodMLR
multiTrait <- self$multiTrait
nSelectionWaysVec <- self$nSelectionWaysVec
selectionMethodList <- self$selectionMethodList
traitNoSelList <- self$traitNoSelList
blockSplitMethod <- self$blockSplitMethod
nMrkInBlock <- self$nMrkInBlock
minimumSegmentLength <- self$minimumSegmentLength
nSelInitOPVList <- self$nSelInitOPVList
nIterOPV <- self$nIterOPV
nProgeniesEMBVVec <- self$nProgeniesEMBVVec
nIterEMBV <- self$nIterEMBV
nCoresEMBV <- self$nCoresEMBV
clusteringForSelList <- self$clusteringForSelList
nClusterList <- self$nClusterList
nTopClusterList <- self$nTopClusterList
nTopEachList <- self$nTopEachList
nSelList <- self$nSelList
multiTraitsEvalMethodList <- self$multiTraitsEvalMethodList
hSelList <- self$hSelList
matingMethodVec <- self$matingMethodVec
allocateMethodVec <- self$allocateMethodVec
weightedAllocationMethodList <- self$weightedAllocationMethodList
traitNoRAList <- self$traitNoRAList
minimumUnitAllocateVec <- self$minimumUnitAllocateVec
includeGVPVec <- self$includeGVPVec
nNextPopVec <- self$nNextPopVec
nameMethod <- self$nameMethod
nCores <- self$nCores
overWriteRes <- self$overWriteRes
showProgress <- self$showProgress
returnMethod <- self$returnMethod
saveAllResAt <- self$saveAllResAt
evaluateGVMethod <- self$evaluateGVMethod
nTopEval <- self$nTopEval
traitNoEval <- self$traitNoEval
hEval <- self$hEval
summaryAllResAt <- self$summaryAllResAt
verbose <- self$verbose
saveObjNameBase <- self$saveObjNameBase
populationNameInit <- names(bsInfoInit$populations[length(bsInfoInit$populations)])
iterNames <- paste0("Iteration_", 1:nIterSimulation)
hVecOptsList <- self$hVecOptsList
lociEffectsInit <- self$lociEffectsInit
if (setGoalAsFinalGeneration) {
nGenerationProceedSimulationNow <- min(nGenerationProceed, nGenerationProceedSimulation)
} else {
nGenerationProceedSimulationNow <- nGenerationProceedSimulation
}
hVecLenNow <- sum(hLens[1:nGenerationProceedSimulationNow])
hVecOpt <- self$solnInit$par
hVecOptsList[["Initial"]] <- hVecOpt
optimalHyperParamMat <- self$optimalHyperParamMatsList[["Initial"]]
if (sameAcrossGeneration) {
if (propSelOpt) {
propSelListOpt <- sapply(X = propSelLens[1:nGenerationProceedSimulation],
FUN = function(propSelLen) {
hVecOpt[propSelOrNot][1:propSelLen]
}, simplify = FALSE)
nSelListOpt <- mapply(FUN = function(propSel, nInd) {
round(nInd * propSel)
},
propSelListOpt,
c(bsInfoInit$populations[[length(bsInfoInit$populations)]]$nInd,
nNextPopVec),
SIMPLIFY = FALSE)
hListOpt <- sapply(X = hLensOri[1:nGenerationProceedSimulation],
FUN = function(hLen) {
hVecOpt[!propSelOrNot][1:hLen]
}, simplify = FALSE)
} else {
propSelListOpt <- NULL
nSelListOpt <- nSelList
hListOpt <- sapply(X = hLens[1:nGenerationProceedSimulation],
FUN = function(hLen) {
hVecOpt[1:hLen]
}, simplify = FALSE)
}
} else {
hVecLenNow <- sum(hLens[1:nGenerationProceedSimulation])
if (propSelOpt) {
propSelListOpt <- split(x = hVecOpt[1:hVecLenNow][propSelOrNot],
f = rep(1:nGenerationProceedSimulation,
propSelLens[1:nGenerationProceedSimulation]))
nSelListOpt <- mapply(FUN = function(propSel, nInd) {
round(nInd * propSel)
},
propSelListOpt,
c(bsInfoInit$populations[[length(bsInfoInit$populations)]]$nInd,
nNextPopVec),
SIMPLIFY = FALSE)
hListOpt <- split(x = hVecOpt[1:hVecLenNow][!propSelOrNot],
f = rep(1:nGenerationProceedSimulation,
hLensOri[1:nGenerationProceedSimulation]))
} else {
propSelListOpt <- NULL
nSelListOpt <- nSelList
hListOpt <- split(x = hVecOpt[1:hVecLenNow],
f = rep(1:nGenerationProceedSimulation,
hLens[1:nGenerationProceedSimulation]))
}
}
iterName <- iterNames[iterNo]
simRes <- list()
simRes$trueGVMatList <- list()
simRes$estimatedGVMatList <- list()
if (is.null(hVecOptsList[[iterName]])) {
hVecOptsList[[iterName]] <- list()
}
bsInfo <- bsInfoInit$clone(deep = FALSE)
breederInfo <- breederInfoInit$clone(deep = FALSE)
lociEffects <- lociEffectsInit
hCount <- 0
# trueGVMatList
if (is.null(self$trueGVMatList[[iterName]])) {
simRes$trueGVMatList[[populationNameInit]] <- self$trueGVMatInit
}
# estimatedGVMatList
if (is.null(self$estimatedGVMatList[[iterName]])) {
simRes$estimatedGVMatList[[populationNameInit]] <- self$estimatedGVMatInit
}
if (!is.null(saveObjNameBase)) {
savedObjFilesAll <- list.files(dirname(saveObjNameBase))
saveObjNameIter <- paste0(saveObjNameBase, "_", iterName)
saveObjNameIterBase <- basename(saveObjNameIter)
fileNameSimResIter <- paste0(saveObjNameIter, "_simulation_results.rds")
fileNameBreederInfoIters <- savedObjFilesAll[
grep(pattern = paste0(saveObjNameIterBase, ".*_breederInfo.rds"),
x = savedObjFilesAll)
]
genSavedBreederInfoIters <- readr::parse_number(
stringr::str_extract(
string = fileNameBreederInfoIters,
pattern = "Generation_.*_breederInfo.rds"
)
)
if (length(genSavedBreederInfoIters) == 0) {
genSavedBreederInfoIterMax <- 0
} else {
genSavedBreederInfoIterMax <- max(genSavedBreederInfoIters)
}
} else {
fileNameSimResIter <- NULL
genSavedBreederInfoIterMax <- 0
}
for (genProceedNo in 1:nGenerationProceed) {
if (is.null(saveObjNameBase)) {
saveObjNameIterGen <- NULL
fileNameSolResIterGen <- ""
} else {
saveObjNameIterGen <- paste0(saveObjNameBase, "_", iterName, "_Generation_", genProceedNo)
fileNameSolResIterGen <- paste0(saveObjNameIterGen, "_solution_results.rds")
fileNameBsInfoIterGen <- paste0(saveObjNameIterGen, "_bsInfo.rds")
fileNameBreederInfoIterGen <- paste0(saveObjNameIterGen, "_breederInfo.rds")
}
if ((genProceedNo >= 2) & (performOptimization[genProceedNo])) {
if (verbose) {
print(paste0("Iteration: ", iterNo, ", Generation: ", genProceedNo,
"; Perform optimization of hyperparameters."))
}
if (setGoalAsFinalGeneration) {
nGenerationProceedSimulationNow <- min(nGenerationProceed - genProceedNo + 1,
nGenerationProceedSimulation)
} else {
nGenerationProceedSimulationNow <- nGenerationProceedSimulation
}
hVecLenNow <- sum(hLens[1:nGenerationProceedSimulationNow])
if (!file.exists(fileNameSolResIterGen)) {
if (!is.null(saveObjNameBase)) {
saveObjNameIterGenBase <- basename(saveObjNameIterGen)
wherePastObjFiles <- grep(pattern = saveObjNameIterGenBase,
x = savedObjFilesAll)
if (length(wherePastObjFiles) >= 1) {
savedObjFilesNow <- savedObjFilesAll[wherePastObjFiles]
savedCurrentObjFile <- savedObjFilesNow[
grep(pattern = ".*_CMA-ES_object.rds",
x = savedObjFilesNow)
]
cmaESNow <- readRDS(file = paste0(dirname(saveObjNameBase), "/",
savedCurrentObjFile))
} else {
cmaESNow <- NULL
}
} else {
cmaESNow <- NULL
}
if (is.null(cmaESNow)) {
cmaESNow <- myBreedSimulatR::cmaES$new(parameter = hStart[1:hVecLenNow],
optimizeFunc = private$maximizeFunc,
nGenerationProceedSimulation = nGenerationProceedSimulation,
lowerBound = hMin[1:hVecLenNow],
upperBound = hMax[1:hVecLenNow],
nIterOptimization = nIterOptimization,
lambda = lambda,
mu = mu,
maximize = TRUE,
saveObjNameBase = saveObjNameBaseInit,
whenToSaveObj = self$whenToSaveObj,
verbose = showProgress)
}
cmaESNow$startOptimization()
saveRDS(object = cmaESNow,
file = cmaESNow$fileNameSaveObj)
optimalHyperParamMat <- cmaESNow$optimalHyperParamMat
hVecOpt <- cmaESNow$optimalParameter
solnInit <- list(value = cmaESNow$optimalValue,
par = hVecOpt)
rm(cmaESNow)
gc(reset = TRUE); gc(reset = TRUE)
solnRes <- c(solnInit,
list(optimalHyperParamMat = optimalHyperParamMat))
if (!is.null(saveObjNameBase)) {
saveRDS(object = solnRes, file = fileNameSolResInit)
}
} else {
solnRes <- readRDS(file = fileNameSolResIterGen)
optimalHyperParamMat <- solnRes$optimalHyperParamMat
hVecOpt <- solnRes$par
}
if (sameAcrossGeneration) {
if (propSelOpt) {
propSelListOpt <- sapply(X = propSelLens[1:nGenerationProceedSimulation],
FUN = function(propSelLen) {
hVecOpt[propSelOrNot][1:propSelLen]
}, simplify = FALSE)
nSelListOpt <- mapply(FUN = function(propSel, nInd) {
round(nInd * propSel)
},
propSelListOpt,
c(bsInfoInit$populations[[length(bsInfoInit$populations)]]$nInd,
nNextPopVec),
SIMPLIFY = FALSE)
hListOpt <- sapply(X = hLensOri[1:nGenerationProceedSimulation],
FUN = function(hLen) {
hVecOpt[!propSelOrNot][1:hLen]
}, simplify = FALSE)
} else {
propSelListOpt <- NULL
nSelListOpt <- nSelList
hListOpt <- sapply(X = hLens[1:nGenerationProceedSimulation],
FUN = function(hLen) {
hVecOpt[1:hLen]
}, simplify = FALSE)
}
} else {
hVecLenNow <- sum(hLens[1:nGenerationProceedSimulation])
if (propSelOpt) {
propSelListOpt <- split(x = hVecOpt[1:hVecLenNow][propSelOrNot],
f = rep(1:nGenerationProceedSimulation,
propSelLens[1:nGenerationProceedSimulation]))
nSelListOpt <- mapply(FUN = function(propSel, nInd) {
round(nInd * propSel)
},
propSelListOpt,
c(bsInfoInit$populations[[length(bsInfoInit$populations)]]$nInd,
nNextPopVec),
SIMPLIFY = FALSE)
hListOpt <- split(x = hVecOpt[1:hVecLenNow][!propSelOrNot],
f = rep(1:nGenerationProceedSimulation,
hLensOri[1:nGenerationProceedSimulation]))
} else {
propSelListOpt <- NULL
nSelListOpt <- nSelList
hListOpt <- split(x = hVecOpt[1:hVecLenNow],
f = rep(1:nGenerationProceedSimulation,
hLens[1:nGenerationProceedSimulation]))
}
}
self$optimalHyperParamMatsList[[iterName]][[paste0("Generation_", genProceedNo)]] <- optimalHyperParamMat
hVecOptsList[[iterName]][[paste0("Generation_", genProceedNo)]] <- hVecOpt
hCount <- 1
} else {
self$optimalHyperParamMatsList[[iterName]][[paste0("Generation_", genProceedNo)]] <- optimalHyperParamMat
hVecOptsList[[iterName]][[paste0("Generation_", genProceedNo)]] <- hVecOpt
hCount <- hCount + 1
}
self$hVecOptsList <- hVecOptsList
if (!file.exists(fileNameSimResIter)) {
if (genProceedNo > genSavedBreederInfoIterMax) {
crossInfoNow <- myBreedSimulatR::crossInfo$new(parentPopulation = bsInfo$populations[[length(bsInfo$populations)]],
nSelectionWays = nSelectionWaysVec[hCount],
selectionMethod = selectionMethodList[[hCount]],
traitNoSel = traitNoSelList[[hCount]],
userSI = NULL,
lociEffects = lociEffects,
blockSplitMethod = blockSplitMethod,
nMrkInBlock = nMrkInBlock,
minimumSegmentLength = minimumSegmentLength,
nSelInitOPV = nSelInitOPVList[[hCount]],
nIterOPV = nIterOPV,
nProgeniesEMBV = nProgeniesEMBVVec[hCount],
nIterEMBV = nIterEMBV,
nCoresEMBV = nCoresEMBV,
clusteringForSel = clusteringForSelList[[hCount]],
nCluster = nClusterList[[hCount]],
nTopCluster = nTopClusterList[[hCount]],
nTopEach = nTopEachList[[hCount]],
nSel = nSelListOpt[[hCount]],
multiTraitsEvalMethod = multiTraitsEvalMethodList[[hCount]],
hSel = hSelList[[hCount]],
matingMethod = matingMethodVec[hCount],
allocateMethod = allocateMethodVec[hCount],
weightedAllocationMethod = weightedAllocationMethodList[[hCount]],
nProgenies = NULL,
traitNoRA = traitNoRAList[[hCount]],
h = hListOpt[[hCount]],
minimumUnitAllocate = minimumUnitAllocateVec[hCount],
includeGVP = includeGVPVec[hCount],
targetGenGVP = targetGenGVPVec[hCount],
nNextPop = nNextPopVec[hCount],
performOCS = performOCSVec[hCount],
targetGenOCS = targetGenOCSVec[hCount],
HeStarRatio = HeStarRatioVec[hCount],
degreeOCS = degreeOCSVec[hCount],
includeGVPOCS = includeGVPOCSVec[hCount],
hOCS = hOCSList[[hCount]],
weightedAllocationMethodOCS = weightedAllocationMethodOCSList[[hCount]],
traitNoOCS = traitNoOCSList[[hCount]],
nCrossesOCS = nCrossesOCSVec[hCount],
nPairs = NULL,
nameMethod = nameMethod,
indNames = NULL,
seedSimRM = NA,
seedSimMC = NA,
selCands = NULL,
crosses = NULL,
verbose = verbose)
bsInfo$nextGeneration(crossInfo = crossInfoNow)
if (updateBreederInfo[genProceedNo]) {
breederInfo$getNewPopulation(bsInfo = bsInfo,
generationNew = NULL,
genotyping = TRUE,
genotypedIndNames = NULL)
if (phenotypingInds[genProceedNo]) {
breederInfo$phenotyper(bsInfo = bsInfo,
generationOfInterest = NULL,
estimateGV = TRUE,
estimatedGVMethod = "lme4",
nRep = nRepForPheno[genProceedNo])
breederInfo$populationsFB[[genProceedNo]]$lmerResList <- NULL
gc(reset = TRUE); gc(reset = TRUE)
if (updateModels[genProceedNo]) {
if (performOptimization[genProceedNo + 1] & (!is.null(saveObjNameBase))) {
saveRDS(object = bsInfo, file = fileNameBsInfoIterGen)
saveRDS(object = breederInfo, file = fileNameBreederInfoIterGen)
}
lociEffects <- private$lociEffects(bsInfo = bsInfo$clone(deep = FALSE),
breederInfo = breederInfo$clone(deep = FALSE),
gcReset = TRUE)
}
}
}
} else if (genProceedNo == genSavedBreederInfoIterMax) {
if (updateModels[genProceedNo]) {
if (performOptimization[genProceedNo + 1] & (!is.null(saveObjNameBase))) {
bsInfo <- readRDS(file = fileNameBsInfoIterGen)
breederInfo <- readRDS(file = fileNameBreederInfoIterGen)
}
lociEffects <- private$lociEffects(bsInfo = bsInfo$clone(deep = FALSE),
breederInfo = breederInfo$clone(deep = FALSE),
gcReset = TRUE)
}
}
if (any(c("all", "summary") %in% returnMethod)) {
populationNameNow <- names(bsInfo$populations)[length(bsInfo$populations)]
trueGVMat <- bsInfo$populations[[length(bsInfo$populations)]]$trueGVMat
simRes$trueGVMatList[[populationNameNow]] <- trueGVMat
# if (breederInfo$generation < bsInfo$generation) {
# breederInfo$getNewPopulation(bsInfo = bsInfo,
# generationNew = bsInfo$generation,
# genotyping = TRUE,
# genotypedIndNames = NULL)
# }
#
# if (is.null(breederInfo$estimatedGVByMLRInfo[[names(bsInfo$populations[length(bsInfo$populations)])]])) {
# breederInfo$estimateGVByMLR(trainingPop = self$trainingPopInit,
# trainingIndNames = self$trainingIndNamesInit,
# testingPop = length(breederInfo$populationsFB),
# methodMLR = self$methodMLRInit,
# multiTrait = self$multiTraitInit,
# alpha = 0.5,
# nIter = 12000,
# burnIn = 3000,
# thin = 5,
# bayesian = TRUE)
# }
# estimatedGVMat <- breederInfo$estimatedGVByMLRInfo[[names(bsInfo$populations[length(bsInfo$populations)])]]$testingEstimatedGVByMLR
genoMatNow <- bsInfo$populations[[length(bsInfo$populations)]]$genoMat
genoMatWithIntNow <- cbind(Intercept = rep(1, nrow(genoMatNow)),
genoMatNow)
estimatedGVMat <- genoMatWithIntNow[, rownames(lociEffects)] %*% lociEffects
simRes$estimatedGVMatList[[populationNameNow]] <- estimatedGVMat
}
}
}
if (!file.exists(fileNameSimResIter)) {
if (!is.null(saveAllResAt)) {
fileNameBsInfoRes <- here::here(saveAllResAt,
paste0(simBsName, "_bsInfo_", iterName, ".rds"))
fileNameBreederInfoRes <- here::here(saveAllResAt,
paste0(simBsName, "_breederInfo_", iterName, ".rds"))
saveRDS(object = bsInfo, file = fileNameBsInfoRes)
saveRDS(object = breederInfo, file = fileNameBreederInfoRes)
}
if ("all" %in% returnMethod) {
simRes$all <- list(bsInfo = bsInfo,
breederInfo = breederInfo)
}
if (any(returnMethod %in% c("summary", "max", "mean", "median", "min", "var"))) {
populationNameNow <- names(bsInfo$populations)[length(bsInfo$populations)]
if (any(c("all", "summary") %in% returnMethod)) {
trueGVMat <- simRes$trueGVMatList[[populationNameNow]]
estimatedGVMat <- simRes$estimatedGVMatList[[populationNameNow]]
} else {
trueGVMat <- bsInfo$populations[[length(bsInfo$populations)]]$trueGVMat
simRes$trueGVMatList[[populationNameNow]] <- trueGVMat
# if (breederInfo$generation < bsInfo$generation) {
# breederInfo$getNewPopulation(bsInfo = bsInfo,
# generationNew = bsInfo$generation,
# genotyping = TRUE,
# genotypedIndNames = NULL)
# }
#
# if (is.null(breederInfo$estimatedGVByMLRInfo[[names(bsInfo$populations[length(bsInfo$populations)])]])) {
# breederInfo$estimateGVByMLR(trainingPop = self$trainingPopInit,
# trainingIndNames = self$trainingIndNamesInit,
# testingPop = length(breederInfo$populationsFB),
# methodMLR = self$methodMLRInit,
# multiTrait = self$multiTraitInit,
# alpha = 0.5,
# nIter = 12000,
# burnIn = 3000,
# thin = 5,
# bayesian = TRUE)
# }
# estimatedGVMat <- breederInfo$estimatedGVByMLRInfo[[names(bsInfo$populations[length(bsInfo$populations)])]]$testingEstimatedGVByMLR
genoMatNow <- bsInfo$populations[[length(bsInfo$populations)]]$genoMat
genoMatWithIntNow <- cbind(Intercept = rep(1, nrow(genoMatNow)),
genoMatNow)
estimatedGVMat <- genoMatWithIntNow[, rownames(lociEffects)] %*% lociEffects
simRes$estimatedGVMatList[[populationNameNow]] <- estimatedGVMat
}
if (evaluateGVMethod == "true") {
trueGVMatNow <- trueGVMat
} else {
trueGVMatNow <- estimatedGVMat
}
# trueGVMatScaled <- apply(X = trueGVMatNow, MARGIN = 2,
# FUN = function(trueGV) {
# return(scale(x = trueGV, center = TRUE,
# scale = as.logical(sd(trueGV))))
# })
trueGVMatInit <- self$trueGVMatInit
trueGVMatScaled <- do.call(what = cbind,
args = sapply(X = 1:ncol(trueGVMatNow),
FUN = function(traitNo) {
trueGVMean <- mean(trueGVMatInit[, traitNo])
trueGVSd <- sd(trueGVMatInit[, traitNo])
if (trueGVSd != 0) {
trueGVScaled <- (trueGVMatNow[, traitNo] - trueGVMean) / trueGVSd
} else {
trueGVScaled <- (trueGVMatNow[, traitNo] - trueGVMean)
}
return(trueGVScaled)
}, simplify = FALSE))
rownames(trueGVMatScaled) <- rownames(trueGVMatNow)
colnames(trueGVMatScaled) <- colnames(trueGVMatNow)
trueEvals <- (trueGVMatScaled[, traitNoEval, drop = FALSE] %*% hEval)[, 1]
if ("max" %in% returnMethod) {
simRes$max <- mean(x = sort(x = trueEvals, decreasing = TRUE)[1:nTopEval])
}
if ("mean" %in% returnMethod) {
simRes$mean <- mean(x = trueEvals)
}
if ("median" %in% returnMethod) {
simRes$median <- median(x = trueEvals)
}
if ("min" %in% returnMethod) {
simRes$min <- mean(x = sort(x = trueEvals, decreasing = FALSE)[1:nTopEval])
}
if ("var" %in% returnMethod) {
simRes$var <- var(x = trueEvals)
}
}
if (iterNo %% nRefreshMemoryEvery == 0) {
rm(trueGVMat); rm(estimatedGVMat); rm(bsInfo); rm(breederInfo)
if (!all(returnMethod == "all")) {
rm(trueGVMatNow); rm(trueGVMatScaled)
}
gc(reset = TRUE); gc(reset = TRUE)
}
saveRDS(object = simRes, file = fileNameSimResIter)
} else {
simRes <- readRDS(file = fileNameSimResIter)
}
return(simRes)
},
# @description Proceed optimized breeding scheme with try-error (one simulation)
#
# @param iterNo [numeric] Iteration No.
performOneSimulationOptTryError = function(iterNo) {
simRes <- try(private$performOneSimulationOpt(iterNo = iterNo),
silent = TRUE)
nIterSimulation <- self$nIterSimulation
iterNames <- paste0("Iteration_", 1:nIterSimulation)
iterName <- iterNames[iterNo]
bsInfoInit <- self$bsInfoInit
populationNameInit <- names(bsInfoInit$populations[length(bsInfoInit$populations)])
if ("try-error" %in% class(simRes)) {
simRes <- list()
simRes$trueGVMatList <- list()
simRes$estimatedGVMatList <- list()
# trueGVMatList
if (is.null(self$trueGVMatList[[iterName]])) {
simRes$trueGVMatList[[populationNameInit]] <- self$trueGVMatInit
}
# estimatedGVMatList
if (is.null(self$estimatedGVMatList[[iterName]])) {
simRes$estimatedGVMatList[[populationNameInit]] <- self$estimatedGVMatInit
}
}
return(simRes)
},
# @description Extract GV matrix as a list from bsInfo & breederInfo objects
#
# @param iterName [character] Iteration Name
extractGVMatList = function(iterName) {
# Read arguments from `self`
simBsName <- self$simBsName
bsInfoInit <- self$bsInfoInit
showProgress <- self$showProgress
evaluateGVMethod <- self$evaluateGVMethod
nTopEval <- self$nTopEval
simBsRes <- self$simBsRes
summaryAllResAt <- self$summaryAllResAt
nIterSimulation <- self$nIterSimulation
iterNames <- paste0("Iteration_", 1:nIterSimulation)
fileNameBsInfoRes0 <- here::here(summaryAllResAt,
paste0(simBsName, "_bsInfo_"))
fileNameBreederInfoRes0 <- here::here(summaryAllResAt,
paste0(simBsName, "_breederInfo_"))
fileNameBsInfoRes <- paste0(fileNameBsInfoRes0, iterName, ".rds")
bsInfoEach <- try(readRDS(file = fileNameBsInfoRes), silent = TRUE)
if (!("try-error" %in% class(bsInfoEach))) {
trueGVMatEachList <- lapply(X = bsInfoEach$populations,
FUN = function(eachPop) {
trueGVMatEachPop <- eachPop$trueGVMat
return(trueGVMatEachPop)
})
} else {
trueGVMatEachList <- NULL
}
# fileNameBreederInfoRes <- paste0(fileNameBreederInfoRes0, iterName, ".rds")
# breederInfoEach <- try(readRDS(file = fileNameBreederInfoRes), silent = TRUE)
# if (!("try-error" %in% class(breederInfoEach))) {
# if (breederInfoEach$generation < bsInfoEach$generation) {
# for (generationAdd in (breederInfoEach$generation + 1):bsInfoEach$generation) {
# breederInfoEach$getNewPopulation(bsInfo = bsInfoEach,
# generationNew = generationAdd,
# genotyping = estimated,
# genotypedIndNames = NULL)
# }
# }
#
# for (generationNow in 1:length(breederInfoEach$populationsFB)) {
# eachPop <- breederInfoEach$populationsFB[[generationNow]]
# if (is.null(breederInfoEach$estimatedGVByMLRInfo[[eachPop$name]])) {
# breederInfoEach$estimateGVByMLR(trainingPop = self$trainingPopInit,
# trainingIndNames = self$trainingIndNamesInit,
# testingPop = generationNow,
# methodMLR = self$methodMLRInit,
# multiTrait = self$multiTraitInit,
# alpha = 0.5,
# nIter = 12000,
# burnIn = 3000,
# thin = 5,
# bayesian = TRUE)
# }
# }
# estimatedGVMatEachList <- lapply(X = breederInfoEach$populationsFB,
# FUN = function(eachPop) {
# estimatedGVMatEachPop <- breederInfoEach$estimatedGVByMLRInfo[[eachPop$name]]$testingEstimatedGVByMLR
#
# return(estimatedGVMatEachPop)
# })
# } else {
# estimatedGVMatEachList <- NULL
# }
estimatedGVMatEachList <- lapply(X = bsInfoEach$populations,
FUN = function(eachPop) {
genoMatNow <- eachPop$genoMat
genoMatWithIntNow <- cbind(Intercept = rep(1, nrow(genoMatNow)),
genoMatNow)
estimatedGVMatEachPop <- genoMatWithIntNow[, rownames(lociEffects)] %*% lociEffects
return(estimatedGVMatEachPop)
})
return(list(trueGVMatEachList = trueGVMatEachList,
estimatedGVMatEachList = estimatedGVMatEachList))
},
# @description Extract summary results from true GV matrix
#
# @param trueGVMatListEach [list] Each list of true GV matrix
extractTrueSummaryRes = function(trueGVMatListEach) {
nTopEval <- self$nTopEval
trueGVSummaryArrayEachList <- lapply(X = trueGVMatListEach,
FUN = function(trueGVMatEachPop) {
trueGVSummaryEachPop <- apply(X = trueGVMatEachPop,
MARGIN = 2,
FUN = function(trueGVMatEachPopEachTrait) {
trueGVSummaryEachPopEachTrait <- c(mean(x = sort(x = trueGVMatEachPopEachTrait, decreasing = TRUE)[1:nTopEval]),
mean(trueGVMatEachPopEachTrait),
median(trueGVMatEachPopEachTrait),
mean(x = sort(x = trueGVMatEachPopEachTrait, decreasing = FALSE)[1:nTopEval]),
var(trueGVMatEachPopEachTrait))
return(trueGVSummaryEachPopEachTrait)
})
trueGVSummaryArrayEachPop <- array(data = trueGVSummaryEachPop,
dim = c(dim(trueGVSummaryEachPop), 1),
dimnames = c(dimnames(trueGVSummaryEachPop),
list(Population = "")))
return(trueGVSummaryArrayEachPop)
})
trueGVSummaryArrayEach <- do.call(what = abind::abind,
args = trueGVSummaryArrayEachList)
trueGVSummaryArrayEach <- array(data = trueGVSummaryArrayEach,
dim = c(dim(trueGVSummaryArrayEach), 1),
dimnames = c(dimnames(trueGVSummaryArrayEach),
list(Iteration = "")))
return(trueGVSummaryArrayEach)
},
# @description Extract summary results from estimated GV matrix
#
# @param estimatedGVMatListEach [list] Each list of estimated GV matrix
extractEstimatedSummaryRes = function(estimatedGVMatListEach) {
nTopEval <- self$nTopEval
estimatedGVSummaryArrayEachList <- lapply(X = estimatedGVMatListEach,
FUN = function(estimatedGVMatEachPop) {
estimatedGVSummaryEachPop <- apply(X = estimatedGVMatEachPop,
MARGIN = 2,
FUN = function(estimatedGVMatEachPopEachTrait) {
estimatedGVSummaryEachPopEachTrait <- c(mean(x = sort(x = estimatedGVMatEachPopEachTrait, decreasing = TRUE)[1:nTopEval]),
mean(estimatedGVMatEachPopEachTrait),
median(estimatedGVMatEachPopEachTrait),
mean(x = sort(x = estimatedGVMatEachPopEachTrait, decreasing = FALSE)[1:nTopEval]),
var(estimatedGVMatEachPopEachTrait))
return(estimatedGVSummaryEachPopEachTrait)
})
estimatedGVSummaryArrayEachPop <- array(data = estimatedGVSummaryEachPop,
dim = c(dim(estimatedGVSummaryEachPop), 1),
dimnames = c(dimnames(estimatedGVSummaryEachPop),
list(Population = "")))
return(estimatedGVSummaryArrayEachPop)
})
estimatedGVSummaryArrayEach <- do.call(what = abind::abind,
args = estimatedGVSummaryArrayEachList)
estimatedGVSummaryArrayEach <- array(data = estimatedGVSummaryArrayEach,
dim = c(dim(estimatedGVSummaryArrayEach), 1),
dimnames = c(dimnames(estimatedGVSummaryArrayEach),
list(Iteration = "")))
return(estimatedGVSummaryArrayEach)
}
)
)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.