Using testParameterOptimality
In BreedingSchemeLanguage: Describe and Simulate Breeding Schemes

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Load the BreedingSchemeLanguage.

library(BreedingSchemeLanguage)

Define the objective function that you will seek to maximize

objFunc <- function(sEnv){
  cross(sEnv, nProgeny=200) # popID 4
  scheme_out <- outputResults(sEnv)
  return(mean(scheme_out[nrow(scheme_out), 1:length(sEnv$sims)]))
}

Find the script, data, and set up objects for results

fileFounderData <- system.file("extdata", "optimizationSchemeData.RData", package = "BreedingSchemeLanguage")
# Take ".RData" off because it gets added
fileFounderData <- substr(fileFounderData, 1, nchar(fileFounderData)-6)
schemeScriptPath <- system.file("extdata", "optimizationScheme.R", package = "BreedingSchemeLanguage")
allowableBudget <- 1000
bestObjective <- -1e6
bestParms <- NULL
allParms <- NULL

Here is the breeding scheme that we want to optimize

Notice that sEnv has to be used explicitly in the breeding scheme

scheme <- readLines(con=schemeScriptPath)
# for (i in 1:10) cat(scheme[i], "\n")

Function to sample parmList at random

sampleParmList <- function(){
  return(list(nRepC0=sample(2, 1), nRepC1=sample(2, 1), nFounder=sample(350, 1) + 50, nProg=sample(350, 1) + 50, yesGeno=runif(1) < 0.5))
}

Sample possible parameter vectors and see how they do

Note that this is a very suboptimal approach to optimization

nIter <- 10
for (i in 1:nIter){
  # Set up a fresh simulation
  if (exists("simEnv")){
    rm(list=names(simEnv), envir=simEnv)
    rm(simEnv)
  }
  simEnv <- defineSpecies(loadData=fileFounderData)
  defineVariances(plotTypeErrVars = c(Standard = 2))
  defineCosts(phenoCost = c(Standard = 1), genoCost = 0.25)

  parmList <- sampleParmList()
  # Get objective value for the scheme
  schemeOut <- testParameterOptimality(schemeFileName = schemeScriptPath, parmList=parmList, objectiveFunc=objFunc, budget=allowableBudget)

  # If over budget, adjust number of progeny, nProg, to keep within budget
  if (is.na(schemeOut$objective)){
    adjust <- allowableBudget / schemeOut$totalCost
    parmList$nFounder <- floor(parmList$nFounder * adjust)
    parmList$nProg <- floor(parmList$nProg * adjust)
    if (exists("simEnv")){
      rm(list=names(simEnv), envir=simEnv)
      rm(simEnv)
    }
    simEnv <- defineSpecies(loadData=fileFounderData)
    defineVariances(plotTypeErrVars = c(Standard = 2))
    defineCosts(phenoCost = c(Standard = 1), genoCost = 0.25)
    schemeOut <- testParameterOptimality(schemeFileName = schemeScriptPath, parmList=parmList, objectiveFunc=objFunc, budget=allowableBudget)
    if (is.na(schemeOut$objective)) stop("Budget allocation problem")
  }

  # If objective better than anything before, save the parmList
  if (schemeOut$objective > bestObjective){
    bestObjective <- schemeOut$objective
    bestParms <- parmList
  }
  allParms <- rbind(allParms, c(i, unlist(parmList), unlist(schemeOut)))
}

Print the output

print(paste("Best Reponse", bestObjective))
print("Parameter List")
print(parmList)
print(head(allParms[order(allParms[,"objective"], decreasing=T),]))