R/sim.R

In Blakemassey/ibmr: Inidividual-Based Modeling (IBM) Simulation Framework for R

# ----------------------- SIMULATION FUNCTIONS ---------------------------------
# General functions for simulation construction and analysis of simulation.
#

#' AgingSubModel
#'
#' @usage AgingSubModel(agent_states, step_data, step)
#'
#' @param agent_states agent_states object
#' @param step_data step_data object
#' @param step step object
#'
#' @return agent_states
#' @export
#'
AgingSubModel <- function(agent_states = agent_states,
                          step_data = step_data,
                          step = step){
  agent_states <- agent_states
  step_data <- step_data
  step <- step
  if (nrow(step_data) == 1){
    agent_states$current_age <- CalculateCurrentAge(datetime=int_start(step),
      birth_date=agent_states$birth_date)
  } else {
    current_step_day <- day(step_data[nrow(step_data), "datetime"])
    previous_step_day <- day(step_data[nrow(step_data)-1, "datetime"])
    if(current_step_day != previous_step_day) {
      agent_states$current_age <- CalculateCurrentAge(datetime=int_start(step),
        birth_date=agent_states$birth_date)
    }
  }
  if (is.na(agent_states$start_datetime)) {
    agent_states[["start_datetime"]] <- as.character(int_start(step))
  }
  return(agent_states)
}

#' CalculateCurrentAge
#'
#' Calculates current age of agent at a given step
#'
#' @usage CalculateCurrentAge(datetime, birthdate)
#'
#' @param datetime datetime of step
#' @param birth_date birthdate of the agent
#'
#' @return current_age in days
#' @export
#'
CalculateCurrentAge <- function(datetime,
                                birth_date) {
  current_age <- as.numeric(difftime(datetime, birth_date,
    tz=tz(sim$pars$global$sim_start)),
    units = c("days"))
  return(current_age)
}

#' CalculateReportAge
#'
#' Creates and updates population interval dataframe
#'
#' @usage CalculateReportAge(agent_states, step_data)
#' @param agent_states = agent_states
#' @param step_data = step data dataframe
#'
#' @return report_date in 'report_age_period' time Period
#' @export
#'
#' @examples
#'
CalculateReportAge <- function(agent_states,
                               step_data) {
  birth_date <- lubridate::ymd(agent_states$birth_date)
  datetime <- step_data[nrow(step_data), "datetime"]
  age_period <- sim$pars$global$report_age_period
  if (age_period == "day" || age_period == "days") {
    age <- floor(as.numeric(difftime(datetime, birth_date, units="days")))
  }
  if (age_period == "week" || age_period == "weeks") {
    age <- floor(as.numeric(difftime(datetime, birth_date, units="weeks")))
  }
  if (age_period == "month" || age_period == "months") {
    age <- floor((zoo::as.yearmon(datetime) - zoo::as.yearmon(birth_date))*12)
  }
  if (age_period == "year" || age_period == "years") {
    age <- floor(zoo::as.yearmon(datetime) - zoo::as.yearmon(birth_date))
  }
  return(age)
}

#' CompileAllAgentsStepData
#'
#' Creates a compiled dataframe of all the agent$step_data objects within
#' sim$agents$all
#'
#' @usage CompileAllAgentsStepData(sim)
#'
#' @param sim  a 'sim' list object
#'
#' @return a dataframe of all the step_data objects
#' @export
#'
#' @examples
#'
CompileAllAgentsStepData <- function(sim = sim) {
  all_step_data <- data.frame()
  all <- sim$agents$all
  for (i in 1:length(all)) {
    if (exists("step_data", where=sim$agents$all[[i]])){
      step_data_i <- sim$agents$all[[i]]$step_data
      step_data_i$sex <- sim$agents$all[[i]]$states$sex
      all_step_data <- rbind(all_step_data, step_data_i)
    }
  }
  return(all_step_data)
}

#' CreateAgentsInputClass
#'
#' Creates a data frame of input agents based on class, sex_ration, and min/max
#' age
#'
#' @usage CreateAgentsInputClass(n, class, sex_ratio, age_min, age_max)
#'
#' @param n number of individuals
#' @param class "class" of individuals
#' @param sex_ratio  numeric [1-0], male/female sex ratio
#' @param age_min integer, minimum age (lubridate 'period' is determined later)
#' @param age_max integer, maximum age (lubridate 'period' is determined later)
#' @param start "random" or "regular"
#' @param base RasterLayer, base
#'
#' @return The input data frame with a column "birth_date" containing each
#' agent's birth_date as a POSIXct object.
#' @export
#'
CreateAgentsInputClass <- function(n = 100,
                                   class = "turtle",
                                   sex_ratio = .5,
                                   age_min = 1,
                                   age_max = 10,
                                   start = "random",
                                   base = base){
  id <- 1:n
  sex <- sample(c("male", "female"), n, replace=TRUE, prob = c(sex_ratio,
    1-sex_ratio))
  age <- sample(age_min:age_max, n, replace=TRUE)
  input <- cbind.data.frame(class, sex, age)
  if (start == "random"){
    base_sample <- raster::sampleRandom(x=base, size=n, na.rm=TRUE, xy=TRUE)
  }
  if (start == "regular"){
     base_sample <- raster::sampleRandom(x=base, size=n, na.rm=TRUE, xy=TRUE,
       sp=TRUE)
  }
  start_x <- base_sample[,1]
  start_y <- base_sample[,2]
  input <- cbind.data.frame(id, class, sex, age, start_x, start_y,
    stringsAsFactors = FALSE)
  return(input)
}

#' CreateBase
#'
#' Creates a 'base' RasterLayer
#'
#' @usage CreateAgentsInputClass(n, class, sex_ratio, age_min, age_max)
#'
#' @param x_min raster x minimum
#' @param y_min raster y minimum
#' @param x_max raster x maximum
#' @param y_max raster y maximum
#' @param cell_size integer, cell size (meters)
#'
#' @return RasterLayer
#' @export
#'
CreateBase <- function(x_min,
                       y_min,
                       x_max,
                       y_max,
                       cell_size){
  base <- raster::raster(xmn=x_min, xmx=x_max, ymn=y_min, ymx=y_max,
    resolution=cell_size, vals=NULL)
  base[] <- 1
#  base[] <- sample(1:10, ncell(base), replace=TRUE)
  return(base)
}

#' CreateBirthDate
#'
#' Creates a column in the input data frame specifying each agent's birth date
#' baed on its specified age
#'
#' @usage CreateBirthDate(input)
#'
#' @param input a data frame containing the agents used to start the model.
#' Must containing a column specifying "age" as a numeric value.
#'
#' @return The input data frame with a column "birth_date" containing each
#' agent's birth_date as a POSIXct object.
#' @export
#'
CreateBirthDate <- function(sim = sim){
  # Only proceed if there is no birth_date column
  sim = sim
  input = sim$agents$input
  if(!"birth_date" %in% colnames(input)){
    # Loop through each row in the input
    input_age_period <- sim$pars$global$input_age_period
    birth_day <- sim$pars$global$birth_day
    input <- tibble::add_column(input, birth_date = NA)
    for(a in 1:nrow(input)){
      # Is the age_period a year?
      if(input_age_period == "year" || input_age_period == "years") {
        # Determine the first sim_start date after the birth_day
        one_year <- as.period(1, "year")
        s0 <- as.Date(sim$pars$global$sim_start - (one_year*input$age[a]))
        # Set the format of the birth_day
        birth_day_format <- tail(guess_formats(birth_day, orders ="dm"), 1)
        birth_day_format <- paste(birth_day_format,"%Y",sep="")
        # Determine the first birth_day after s0
        s1 <- as.Date(paste(birth_day,year(s0),sep=""), format=birth_day_format)
        if(s0 >= s1) {
          input$birth_date[a] <- as.character(s1)
        } else {
          input$birth_date[a] <- as.character(s1-one_year)
        }
      } else {
        # If age period is not a year
        age_period_unit <- as.period(1, input_age_period)
        input$birth_date[a] <- as.character(sim$pars$global$sim_start -
          (age_period_unit*input$age[a]))
      }
    }
  }
  return(input)
}

#' CreateParsClassConstant
#'
#' Creates a parameter class 'constant'
#'
#' @usage CreateBirthDate(step_cauchy_mu, step_cauchy_rho)
#'
#' @param step_cauchy_mu step cauchy mu
#' @param step_cauchy_rho step cauchy rho
#'
#' @return list
#' @export
#'
CreateParsClassConstant <- function(step_cauchy_mu = 0,
                                    step_cauchy_rho = .5){
  step_cauchy_mu <- step_cauchy_mu
  step_cauchy_rho <- step_cauchy_rho
  fixed <- NamedList(step_cauchy_mu, step_cauchy_rho)
  constant <- NamedList(fixed)
  return(constant)
}

#' CreateParsClassSeason
#'
#' Creates a parameter class 'season'
#'
#' @usage CreateParsClassSeason(step_pareto_scale,
#'   step_pareto_shape, step_max_r)
#'
#' @param step_pareto_scale pareto scale
#' @param step_pareto_shape pareto shape
#' @param step_max_r max radius
#'
#' @return list
#' @export
#'
CreateParsClassSeason <- function(step_pareto_scale,
                                  step_pareto_shape,
                                  step_max_r){
  step_pareto_scale = step_pareto_scale
  step_pareto_shape = step_pareto_shape
  step_max_r = step_max_r
  season <- NamedList(step_pareto_scale, step_pareto_shape, step_max_r)
  return(season)
}

#' CreateParsClassSeason
#'
#' Creates a parameter class 'season'
#'
#' @usage CreateParsClassSeason(input)
#'
#' @param sim_start POSIXct default = as.POSIXct("2017-01-01", tz = "UTC")
#' @param sim_period 'period', default = period(100, "days")
#' @param sim_end POSIXct, default = as.POSIXct("2017-01-01", tz = "UTC")
#' @param rep_period 'period', default = period(1, "days")
#' @param rep_interval character string of dates, default = c("01Jan", "01Jul")
#' @param rep_interval_custom = character string of dates, dafault = NULL
#' @param step_period 'period', default = period(1, "day")
#' @param time_step_period 'period', default = period(6, "hour")
#' @param birth_day = character string of date, default = "01Jan"
#' @param input_age_period = character string of period, default = "years"
#' @param report_age_period = character string of period, default = "months"
#' @param sim_seasons = dataframe, columns = c(start, season), default are
#'   typical sesasons
#'
#' @import lubridate
#'
#' @return list
#' @export
#'
CreateParsGlobal <- function(sim_start = as.POSIXct("2015-01-01", tz = "UTC"),
                             sim_period = period(100, "days"),
                             sim_end = NULL,
                             rep_period = period(1, "days"),
                             rep_interval = c("01Jan", "01Jul"),
                             rep_interval_custom = NULL,
                             step_period = period(1, "day"),
                             time_step_period = period(6, "hour"),
                             birth_day = "01Jan",
                             input_age_period = "years",
                             report_age_period = "months",
                             sim_seasons = NULL){
  sim_start <- sim_start
  sim_period <- sim_period
  sim_end <- sim_end
  rep_period <- rep_period
  rep_interval <- rep_interval
  rep_interval_custom <- rep_interval_custom
  step_period <- step_period
  time_step_period <- time_step_period
  birth_day <- birth_day
  input_age_period = input_age_period
  report_age_period = report_age_period
  if (is.null(sim_seasons)){
    sim_seasons <- data.frame(start = c("20Mar", "21Jun", "23Sep", "22Dec"),
      season = c("spring", "summer", "winter", "fall"), stringsAsFactors=FALSE)
  }
  global <- NamedList(sim_start, sim_period, sim_end, rep_period, rep_interval,
    rep_interval_custom, step_period, time_step_period, sim_seasons, birth_day,
    input_age_period, report_age_period)
  return(global)
}

#' CreateRunsList
#'
#' Helper function for RunSimulation(). Creates a list of empty objects, each
#' one named individually as "run_xx" where the number of leading zeroes matches
#' the maximum digit width of the number of runs.
#'
#' @param runs = numeric, number of runs in RunSimulation()
#'
#' @return a list
#' @export
#'
#' @note Where there are < 10 runs, there will have no leading zeros, when
#' there are between 10 and 99 runs, there will be a leading zero for runs
#' numbered < 10.
CreateRunsList <- function(runs = runs) {
  format <- paste0("%0", nchar(runs), "d")
  runs_list <- list()
  for (i in 1:runs) {
    runs_list[[i]] <- NA
    names(runs_list)[[i]] <- paste0("run_", sprintf(format, i))
  }
  return(runs_list)
}

#' CreateReportIntervals
#'
#' Creates summary intervals for the simulation input
#'
#' @usage CreateReportIntervals(sim)
#'
#' @param sim a list that contains a "pars" list which contains a "global" list
#'   that sets the arguments for the Report Interval. See Description section
#'   below.
#'
#' @return list of intervals
#' @export
#'
#' @description
#' sim list that contains a "pars" list which contains a "global" list that
#' contains the following objects:
#'   sim_start = a POSIXct object for simulation's start, required
#'   sim_period = a Period object (from 'lubridate' package), e.g. period(10,
#'     "year"), period(10, "week"), period(10, "week") that sets the length of
#'     time the simulation runs
#'   sim_end = a POSIXct object for simulation's end. Setting this will override
#'     the 'sim_period' parameter. Default = NULL.
#'   rep_period = a Period object that sets the summary interval
#'   rep_interval = vector or dateframe (with date in first column) of the
#'     format %B%d, %b%d, %d%B, or %B%d (see ?strptime for more details on
#'     formatting). Example: c("April01", "May15", "Sep1", "Oct15"). Overrides
#'     'rep_period' parameter.
#'   rep_interval_custom = a vector of POSIXct objects that set the start of the
#'     summary periods. If sim_start precedes the first value, the first value
#'     will go until from start_sim to rep_interval_custom. Last interval will
#'     go from the last date in rep_interval_custom to end of sim_period or
#'     end_sim. Overrides 'rep_period' or 'rep_interval' parameters.
#'
#' @note Either 'sim_period' or 'sim_end' must be set, and either 'rep_period',
#' 'rep_interval', or 'rep_interval_custom' must be set. For 'rep_interval' the
#' dates can be set as numeric for both day and month, but there is a
#' possibility that the day and month will be inverted because dates with both
#' values <12 can be confused (e.g. 02-10 can be either February 10th or October
#' 2nd. Date format: "10Feb" or "02Oct" is preferred.
CreateReportIntervals <- function(sim = sim) {
  sim_start <- sim$pars$global$sim_start
  sim_period <- sim$pars$global$sim_period
  sim_end <- sim$pars$global$sim_end
  rep_period <- sim$pars$global$rep_period
  rep_interval <- sim$pars$global$rep_interval
  rep_interval_custom <- sim$pars$global$rep_interval_custom
  if (is.null(sim_period) && is.null(sim_end)) {
    stop("must provide either 'sim_period' or 'sim_end' parameter value")
  }
  if (!is.null(sim_period))  sim_end <- sim_start + sim_period
  if (!is.null(rep_period) && !is.null(rep_interval)){
    print("'rep_period' used instead of 'rep_interval'")
  }
  if (!is.null(rep_period)) {
    rep_period_unit <- ExtractUnitFromPeriod(rep_period)
    first_int <- as.interval(sim_start, ceiling_date(sim_start+period(1,
      "second"), unit=rep_period_unit))
    rep_intervals <- sim_start
    end_int <- int_end(first_int)
    rep_intervals <- with_tz(append(rep_intervals, end_int), tz(sim_start))
    while (end_int < sim_end) {
      end_int <- end_int + rep_period
      rep_intervals <- with_tz(append(rep_intervals, end_int), tz(sim_start))
    }
    if (tail(rep_intervals, 1) > sim_end) {
      rep_intervals[length(rep_intervals)] <- sim_end
    }
    rep_interval <- NULL
  }
  if (!is.null(rep_interval)){
      if (is.data.frame(rep_interval)) rep_interval <- rep_interval[, 1]
      if (max(names(guess_formats(rep_interval, c("md", "dm")))) == "md") {
        md_format <- max(guess_formats(rep_interval, c("md", "dm")))
        rep_interval <- as.character(as.Date(rep_interval, format = md_format),
          "%d%b")  # date must be in day month order
      }
      dates <- dmy(paste0(rep_interval, 2000)) # creates POSIXct
      rep_interval  <- rep_interval[order(dates)]  # orders rep_interval dates
      first_rep_int <- FindFirstReportInterval(sim_start, rep_interval)
      end_int <- int_end(first_rep_int)
      rep_intervals <- with_tz(append(sim_start, end_int), tz(sim_start))
      while (end_int < sim_end) {
        end_int_jul <- yday(end_int)
        rep_int_jul <- yday(dmy(paste0(rep_interval, year(end_int))))
        rep_int_position <- match(end_int_jul, rep_int_jul)
        if (rep_int_position == length(rep_interval)) {
          next_int <- interval(end_int, dmy(paste0(rep_interval[1],
            (year(end_int)))) + period(1, "year"))
        } else {
          next_int <- interval(end_int, dmy(paste0(rep_interval
            [rep_int_position+1], year(end_int))))
        }
        end_int <- int_end(next_int)
        rep_intervals <- with_tz(append(rep_intervals, end_int), tz(sim_start))
      }
      if (tail(rep_intervals, 1) > sim_end) {
        rep_intervals[length(rep_intervals)] <- sim_end
      }
      rep_interval_custom <- NULL
  }
  if(!is.null(rep_interval_custom)){
    if (is.data.frame(rep_interval)) rep_interval <- rep_interval_custom[,1]
    rep_interval <- rep_interval_custom # needed if rep_interval_custom != df
    rep_intervals <- as.POSIXct(rep_interval, tz=tz(sim_start))
    if (tail(rep_intervals, 1) > sim_end && !is.null(rep_period)) {
      rep_intervals[length(rep_intervals)] <- sim_end
      warning("'sum_inverval_custom' exceeds 'sim_period' time length")
    }
    if (tail(rep_intervals, 1) > sim_end && is.null(rep_period)) {
      rep_intervals[length(rep_intervals)] <- sim_end
      warning("'sum_inverval_custom' exceeds 'sim_end' time length")
    }
    if (tail(rep_intervals, 1) < sim_end) {
      rep_intervals[length(rep_intervals)+1] <- sim_end
    }
  }
  rep_intervals_list <- list()
  for (i in 1:(length(rep_intervals)-1)) {
    interval <- as.interval(rep_intervals[i], rep_intervals[i+1])
    rep_intervals_list[[length(rep_intervals_list)+1]] <- interval
  }
  return(rep_intervals_list)
}

#' CreateStepIntervals
#'
#' Create step intervals within a summary interval
#'
#' @usage CreateStepIntervals(rep_interval, step_period)
#'
#' @param rep_interval = an Interval object from lubridate that represents the
#' current summary interval
#' @param step_period = a Period object specifying the length of the step
#' period. The step_period should be of equal or greater length than the
#' time_step_period.
#'
#' @return a list of intervals
#' @export
#'
CreateStepIntervals <- function(rep_interval = rep_interval,
                                step_period = sim$pars$global$step_period) {
  step_period <- step_period
  step_intervals <- list()
  interval_counter <- 1
  current_start <- lubridate::int_start(rep_interval)
  current_end <- (current_start + step_period)
  stop_point <- lubridate::int_end(rep_interval)
  while(current_start < (stop_point)) {
    current_end <- (current_start+step_period)
    step_intervals[[interval_counter]] <- lubridate::interval(current_start,
      current_end)
    interval_counter <- interval_counter + 1
    current_start <- current_start + step_period
    }
  if (lubridate::int_end(step_intervals[[length(step_intervals)]])>stop_point){
    step_intervals[[length(step_intervals)]] <-
      interval(lubridate::int_start(step_intervals[[length(step_intervals)]]),
        stop_point)
    }
  return(step_intervals)
}

#' CreateTimeSteps
#'
#' Create the time steps within a step interval
#'
#' @usage CreateTimeSteps(step_interval)
#'
#' @param step_interval = an Interval object from lubridate that represents the
#' current step interval
#' @param time_step_period = a Period object specifying the length of the time
#' step period. Note that if the length of the current_step_interval is equal to
#' the time_step_period the function will return a POSIXct object with the same
#' date as the end of the current_step_interval.
#'
#'
#' @return a list of intervals
#' @export
#'
CreateTimeSteps <- function(step_interval = step_interval,
                            time_step_period =
                              sim$pars$global$time_step_period) {
  options(lubridate.verbose=FALSE)
  step_interval_end <- lubridate::int_end(step_interval)
  step_interval_start <- lubridate::int_start(step_interval)
  end_int <- lubridate::int_end(lubridate::as.interval(time_step_period,
    step_interval_start))
  steps <- lubridate::with_tz(append(step_interval_start, end_int),
    lubridate::tz(lubridate::int_end(step_interval)))
  while (end_int < step_interval_end) {
    end_int <- end_int + time_step_period
    steps <- lubridate::with_tz(append(steps, end_int),
      lubridate::tz(lubridate::int_end(step_interval)))
  }
  if (tail(steps, 1) > step_interval_end) {
    steps[length(steps)] <- step_interval_end
  }
  time_step_list <- list()
  for (i in 1:(length(steps)-1)) {
    interval <- lubridate::as.interval(steps[i], steps[i+1])
    time_step_list[[length(time_step_list)+1]] <- interval
  }
  return(time_step_list)
}

#' CustomAgentsReportData
#'
#' Creates custom report data
#'
#' @param agent_states = agent_states
#' @param report_data = report_data
#' @param step_data = step_data
#' @param q = counter for reporting interval
#'
#' @return 'report_data' object
#' @export
#'
CustomAgentsReportData <- function(agent_states,
                                   report_data,
                                   step_data,
                                   q) {
    report_data <- report_data
#    report_data[q, "sex"] <- agent_states$sex
#    report_data[q, "move_length"] <- sum(step_data$step_length, na.rm=TRUE)
    return(report_data)
}

#' CustomPopReportData
#'
#' Creates and updates custom pop_report data.
#'
#' @usage CustomPopReportData(report_data, q, agent_states, step_data)
#'
#' @param pop_report = population report
#' @param report_data =  report_data
#' @param q = report_data row for new data
#'
#' @return a population report
#' @export
#'

CustomPopReportData <- function(pop_report,
                                report_data,
                                q) {
  pop_report <- pop_report
#  pop_report[q, "mf_ratio"] <- nrow(subset(report_data, sex=="male")) /
#    nrow(subset(report_data, sex=="female"))
  return(pop_report)
}

#' ExtractUnitFromPeriod
#'
#' Helper function for CreateSumInterval(). Finds unit of Period object.
#'
#' @usage ExtractUnitFromPeriod(period)
#'
#' @param period = period object
#'
#' @return period object's unit as a chararcter
#' @export
#'
ExtractUnitFromPeriod <- function(period) {
  period <- period
  if (period$year > 0) unit <- "year"
  if (period$month > 0) unit <- "month"
  if (period$day == 7 ) unit <- "week"
  if (period$day > 0 && period@day < 7) unit <- "day"
  if (period$day > 7) unit <- "day"
  if (period$hour > 0) unit <- "hour"
  if (period$minute > 0) unit <- "minute"
  if (period$.Data > 0) unit <- "second"
  return(unit)
}

#' FindFirstReportInterval
#'
#' Helper function for CreateSumInterval. Finds first interval based on
#' sim_start datetime and rep_interval dates
#'
#' @usage FindFirstReportInterval(sim_start, rep_interval)
#'
#' @param sim_start  = a POSIXct object for simulation's start, required
#' @param rep_interval = vector or dateframe (with date in first column) of the
#'   format %B%d, %b%d, %d%B, or %B%d (see ?strptime for more details on
#'   formatting). Example: c("April01", "May15", "Sep1", "Oct15"). Required.
#'
#' @import lubridate
#'
#' @return an Interval object
#' @export
FindFirstReportInterval <- function(sim_start,
                                    rep_interval){
  start_year <- as.Date(0, origin=as.Date(floor_date(sim_start, "year")))
  if (length(rep_interval) == 1) {
    first_int <- as.interval(period(1, "year"), dmy(paste0(rep_interval[1],
      year(start_year))))
    if (sim_start %within% first_int) {
      first_rep_int <- interval(sim_start, dmy(paste0(rep_interval[1],
        year(start_year + period(1, "year")))))
    } else {
      first_rep_int <- interval(sim_start, dmy(paste0(rep_interval[1],
        year(start_year))))
    }
  }
  if (length(rep_interval) == 2) {
    first_int <- interval(dmy(paste0("0101", year(start_year))),
      dmy(paste0(rep_interval[1], year(start_year))))
    if ((sim_start + period(1, "second")) %within% first_int) {
      first_rep_int <- interval(sim_start, dmy(paste0(rep_interval[1],
        year(start_year))))
    }
    if (exists("first_rep_int") == FALSE) {
      second_int <- interval(dmy(paste0(rep_interval[1], year(start_year))),
        dmy(paste0(rep_interval[2], year(start_year))))
      if ((sim_start + period(1, "second")) %within% second_int) {
        first_rep_int <- interval(sim_start, dmy(paste0(rep_interval[2],
          year(start_year))))
      } else {
        first_rep_int <- interval(sim_start, dmy(paste0(rep_interval[1],
          year(start_year+period(1, "year")))))
      }
    }
  }
  if (length(rep_interval) > 2) {
    first_int <- interval(dmy(paste0("0101", year(start_year))),
        dmy(paste0(rep_interval[1], year(start_year))))
    if ((sim_start + period(1, "second")) %within% first_int) {
      first_rep_int <- interval(sim_start, dmy(paste0(rep_interval[1],
        year(start_year))))
    }
    for (i in 2:(length(rep_interval)-1)) {
      mid_int <- interval(dmy(paste0(rep_interval[i], year(start_year))),
        dmy(paste0(rep_interval[i+1], year(start_year))))
      if ((sim_start + period(1, "second")) %within% mid_int) {
        first_rep_int <- interval(sim_start, dmy(paste0(rep_interval[i+1],
          year(start_year))))
      }
    }
    for (i in length(rep_interval)) {
      last_int <- interval(dmy(paste0(rep_interval[i], year(start_year))),
        dmy(paste0("0101", (year(start_year)+1) )))
      if ((sim_start + period(1, "second")) %within% last_int) {
        first_rep_int <- interval(sim_start, dmy(paste0(rep_interval[1],
          (year(start_year)+1))))
      }
    }
  }
  return(first_rep_int)
}

#' FindSeasonFromDatetime
#'
#' Determines which season (from sim$pars$global$sim_seasons) a datetime is
#' within
#'
#' @usage FindSeasonFromDatetime(datetime, season)
#'
#' @param datetime = a POSIXct object
#' @param seasons = a dataframe of starting dates ("start" column) and season
#' names ("season" column), usually located at: sim$pars$global$sim_seasons
#'
#' @import lubridate
#'
#' @return  an atomic character object
#' @export
#'
#'
FindSeasonFromDatetime <- function(datetime = datetime,
                                   seasons = sim$pars$global$sim_seasons) {
  start_year <- as.Date(0, origin=as.Date(floor_date(datetime, "year")))
  if (max(names(guess_formats(seasons[, "start"], c("md", "dm")))) == "md") {
    md_format <- max(guess_formats(seasons[, "start"], c("md", "dm")))
    seasons[, "start"] <- as.character(as.Date(seasons[, "start"], format =
      md_format), "%d%b")  # date must be in day month order
  }
  dates <- dmy(paste0(seasons[,"start"], 2000)) # creates POSIXct
  seasons  <- seasons[with(seasons, order(dates)), ]  # orders seasons by dates
  if (nrow(seasons) == 2) {
    first_season <- interval(dmy(paste0(paste0(seasons[1, "start"],
      year(start_year)))), dmy(paste0(seasons[2, "start"], year(start_year)))-1)
    if (datetime %within% first_season) {
      season <- seasons[1, "season"]
    } else {
      season <- seasons[2, "season"]
    }
  }
  if (nrow(seasons) > 2) {
    last_season <- TRUE
    for (i in 1:(nrow(seasons)-1)) {
      mid_season <- interval(dmy(paste0(seasons[i, "start"], year(
        start_year))), dmy(paste0(seasons[i+1, "start"], year(start_year)))-1)
      if (datetime %within% mid_season) {
        season <- seasons[i, "season"]
        last_season = FALSE
      }
    }
    if(last_season == TRUE) season <- seasons[nrow(seasons), "season"]
  }
  return(season)
}



#' NamedList
#'
#' Creates a list with named objects, and it tries not to replace any
#' already-named arguments.
#'
#' @usage NamedList(...)
#'
#' @param ... = objects to add to list
#'
#' @return  a list
#' @export
#'
#' @note Original function came from Ben Bolker's answer on StackOverflow:
#' "http://stackoverflow.com/questions/16951080"
#'
NamedList <- function(...) {
  list <- list(...)
  str_name <- sapply(substitute(list(...)),deparse)[-1]
  if (is.null(name <- names(list))) name <- str_name
  if (any(no_names <- name == "")) name[no_names] <- str_name[no_names]
    setNames(list, name)
}

#' ReproductionSubModel
#'
#' Reproduction submodel
#'
#' @param agent_states = a 'agent_states' list object
#' @param step_data = list object
#'
#' @return 'agent_states' list object
#' @export
#'
#' @note Just a placeholder for now.
#'
ReproductionSubModel <- function(agent_states = agent_states,
                                 step_data = step_data) {
  dummy = TRUE
  if(dummy){
    agent_states = agent_states
  } else {

  }
  return(agent_states)
}

#' ReturnActiveSeq
#'
#' Create vector of positions of agents that had a step_data$datetime within
#' the current reporting interval
#'
#' @param all = 'all' list
#' @param rep_interval = reporting interval
#'
#' @return 'active_seq' vector
#' @export
#'
ReturnActiveSeq <- function(all = sim$agents$all,
                            rep_interval = rep_interval){
  all <- all
  active_seq <- vector()
  for (i in 1:length(all)) {
    step_data <- all[[i]]$step_data
#    step_data <- sim$agents$all[[i]]$step_data
#    print(paste("The step_data:", step_data))
#    print(paste("The active_seq:", rep_interval))
    if (any(step_data$datetime %within% rep_interval)) {
      active_seq <- alive_seq <- append(active_seq, i)
    }
  }
  return(active_seq)
}

#' ReturnAliveIds
#'
#' Create vector of ids of alive agents
#'
#' @param all = 'all' list
#'
#' @return a vector
#' @export
#'
ReturnAliveIds <- function(all = sim$agents$all){
  all <- all
  alive_ids <- vector()
  for (i in 1:length(all)) {
    agent <- all[[i]]
    if (is.na(agent$states$died)) alive_ids <- append(alive_ids,agent$states$id)
  }
  return(alive_ids)
}

#' ReturnAliveSeq
#'
#' Create vector of alive agents positions
#'
#' @param sim = 'sim' list
#'
#' @return a vector
#' @export
#'
ReturnAliveSeq <- function(sim = sim){
  all <- sim$agents$all
  alive_seq <- vector()
  for (i in 1:length(all)) {
    agent <- all[[i]]
    if (is.na(agent$states$died)) alive_seq <- append(alive_seq, i)
  }
  return(alive_seq)
}

#' ReturnPopReportData
#'
#' Helper function for UpdatePopReport(). Finds intervals with start dates
#' within the report interval.
#'
#' @usage ReturnPopReportData(sim, agents, rep_interval)
#'
#' @param sim = 'sim' list
#' @param agents = agents object
#' @param rep_interval = counter for reporting interval (q)
#'
#' @return period object's unit as a chararcter
#' @export
#'
ReturnPopReportData <- function(sim=sim,
                                agents = agents,
                                rep_interval=rep_interval) {
  sim = sim #
  rep_interval = rep_interval #
  report_data = agents$agents_report
#  print(paste0("Inside PopReportData ", report_data))
#  print(paste0("RepInterval ", rep_interval))
  pop_report_data <- data.frame()
  for (s in 1:nrow(report_data)) {
    if (as.POSIXct(report_data[s, "start"], tz = tz(sim$pars$global$sim_start))
      %within% rep_interval) {
        pop_report_data <- rbind(pop_report_data, report_data[s, ])
    }
  }
  return(pop_report_data)
}

#' ReturnReportStepData
#'
#' Returns agent's step data within the current reporting interval
#'
#' @param agent agent = 'agent' list
#' @param rep_interval = counter for reporting interval (q)
#'
#' @return a step dataframe
#' @export
#'
ReturnReportStepData <- function(agent=agent,
                                 rep_interval = rep_interval) {
  step_data = agent$step_data
  report_step_data <- data.frame()
  for (s in 1:nrow(step_data)) {
    if (step_data[s, "datetime"] %within% rep_interval) report_step_data<-
      rbind(report_step_data, step_data[s, ])
  }
  return(report_step_data)
}

#' SurvivalSubModel
#'
#' Survival submodel
#'
#' @param agent_states = 'agent_states' list object
#' @param step_data = list object
#'
#' @return agent_states
#' @export
#'
#' @note Just a placeholder for now.
SurvivalSubModel <- function(agent_states = agent_states,
                             step_data = step_data) {
  agent_states <- agent_states
  return(agent_states)
}

#' UpdateAgentParsData
#'
#' Creates and updates agents parameters dataframe
#'
#' @usage CreateAgentParsData(sim, init)
#'
#' @param sim = sim list, MUST be included in function's parameter argument
#' needed when: when init == TRUE
#' @param init = logical, whether or not this is the initation step
#'
#' @return 'sim' list
#' @export
#'
UpdateAgentParsData <- function(sim = sim,
                                init = FALSE){
  sim=sim
  if (init == TRUE) {
    input <- sim$agents$input
    all <- sim$agents$all
    for (i in 1:nrow(input)){
      agent <- all[[i]]
      pars_data <- data.frame(id=agent$states$id, rep_int = NA)
      all[[i]] <- append(agent, NamedList(pars_data))
    }
    sim$agents$all <- all
    return(sim)
  } else {
    pars_data <- sim$agents$all
    # placeholder for updates
    sim$agents$all <- pars_data
    return(sim)
  }
}

#' UpdateAgentsReport
#'
#' Creates and updates agents interval data dataframe
#'
#' @usage UpdateAgentsReport(agents, init)
#'
#' @param sim = 'sim' list
#' @param rep_interval = rep_interval
#' @param step_intervals = step_intervals
#'
#' @return an 'agents' list
#' @export
#'
#' @note rep_data, init are internal to RunSimualtion()
UpdateAgentsReport <- function(sim = sim,
                               rep_interval = rep_interval,
                               step_intervals = step_intervals,
                               custom_agent_report = FALSE) {
  agents <- sim$agents
  alive_seq  <- ReturnActiveSeq(all=sim$agents$all, rep_interval=rep_interval)
  for (p in alive_seq) {
    if(!("report_data" %in% names(agents$all[[p]]))) {
      report_data <- data.frame()
    } else {
      report_data <- agents$all[[p]]$report_data
    }
    agent <- agents$all[[p]]
    q <- nrow(report_data) + 1
    agent_states <- agent$states
    step_data <- ReturnReportStepData(agent, rep_interval)
    pars_data <- agent$pars
    report_data[q, "start"] <- as.character(int_start(step_intervals[[1]]))
    report_data[q, "end"] <- as.character(int_end(step_intervals[[length(
      step_intervals)]]))
    report_data[q, "id"] <- agent_states$id
    report_data[q, "age"] <- CalculateReportAge(agent_states, step_data)
    if (custom_agent_report == TRUE) {
      report_data <- CustomAgentsReportData(agent_states, report_data,
        step_data, q)
    }
    agents$all[[p]][["report_data"]] <- report_data
  }
  if(!("agents_report" %in% names(agents))) {
    agents_report <- data.frame()
    print("Created agents_report")
  } else {
    agents_report <- sim$agents$agents_report
  }
  active_seq <- ReturnActiveSeq(all=sim$agents$all, rep_interval=rep_interval)
  for (r in active_seq) {
    report_data <- agents$all[[r]]$report_data
    agents_report  <- rbind(agents_report, report_data[nrow(report_data), ])
  }
  print("Updated agents_report")
  agents[["agents_report"]] <- agents_report
   return(agents)
}

#' UpdateAgentStates
#'
#' Creates and updates agents list from the sim object
#'
#' @param agent_states a list from sim$agents$all$agent$states
#' @param sim sim list, MUST be included in function's parameter arguments
#' needed when: when init == TRUE
#' @param init logical, whether or not this is the initation step
#'
#' @return an 'agents' list
#' @export
#'
UpdateAgentStates <- function(agent_states = NULL,
                              sim = sim,
                              init = FALSE) {
  if (init == TRUE) {
    input <- sim$agents$input
    input <- CreateBirthDate(sim)
    input_columns <- colnames(input)
    na_columns <- c("start_datetime", "died")
    all <- list()
    for (i in 1:nrow(input)) {
      states <- list()
      for (j in input_columns) states <- append(states, input[i, j])
      for (k in 1:length(na_columns)) states <- append(states, NA)
      states <- setNames(states, c(input_columns, na_columns))
      agent <- NamedList(states)
      all <- append(all, NamedList(agent))
    }
    sim$agents <- append(sim$agents, NamedList(all))
    return(sim)
  } else {
    agent_states <- agent_states
    return(agent_states)
  }
}

#' UpdateAgentStepData
#'
#' Creates and updates agents step data database
#'
#' @usage UpdateAgentStepData(agents, sim, init)
#'
#' @param step_data = 'step data' dataframe
#' @param sim = 'sim' object
#' @param init = logical, whether or not this is the initation step
#'
#' @return an 'agents' list
#' @export
#'
UpdateAgentStepData <- function(step_data = NULL,
                                sim = sim,
                                init = FALSE) {
  if (init == TRUE) {
    sim_start <- sim$pars$global$sim_start
    all <- sim$agents$all
    for (i in 1:length(all)) {
      agent <- all[[i]]
      step_data <- data.frame(id=agent$states$id, datetime=sim_start,
        x=agent$states$start_x, y=agent$states$start_y, exp_angle=NA,
        abs_angle=NA)
      agent  <-  append(agent, NamedList(step_data))
      all[[i]] <- agent
    }
    sim$agents$all <- all
    return(sim)
  } else {
    step_data <- step_data
    return(step_data)
  }
}



#' UpdatePopReport
#'
#' Creates and updates population interval dataframe
#'
#' @usage UpdatePopReportData(sim, rep_interval, step_intervals)
#'
#' @param sim  'sim' list
#' @param rep_interval rep_interval
#' @param step_intervals  step_intervals
#'
#' @return a 'pop_report' dataframe
#' @export
#'
UpdatePopReport <- function(sim = sim,
                            rep_interval=rep_interval,
                            step_intervals=step_intervals,
                            custom_pop_report = FALSE) {
  agents <- sim$agents
  rep_interval = rep_interval
  if(!("pop_report" %in% names(agents))) {
    pop_report <- data.frame()
    print("Created pop_report")
  } else {
    pop_report <- sim$agents$pop_report
  }
  q <- nrow(pop_report) + 1
  report_data <- ReturnPopReportData(sim=sim, agents=agents, rep_interval=rep_interval)
  pop_report[q, "start"] <- as.character(int_start(step_intervals[[1]]))
  pop_report[q, "end"] <- as.character(int_end(step_intervals[[length(
    step_intervals)]]))
  pop_report[q, "total_n"] <- length(unique(report_data$id))
  if (custom_pop_report == TRUE) {
    pop_report <- CustomPopReportData(pop_report, report_data, q)
  }
  agents[["pop_report"]] <- pop_report
  return(agents)
}

#' UpdateSpatial
#'
#' Creates and updates population interval dataframe
#'
#' @usage UpdateSpatial(spatial, init)
#'
#' @param sim = sim list, MUST be included in function's parameter argument
#' needed when: when init == TRUE
#' @param init = logical, whether or not this is the initation step
#'
#' @return an 'agents' list
#' @export
#'
UpdateSpatial <- function(sim = sim,
                          init = TRUE) {
  sim <- sim
  if (init == TRUE) {
    spatial <- sim$spatial
    sim$spatial <- spatial
    return(sim)
  } else {
#   spatial_timer <- UpdateSpatialTimer(spatial_timer)
#   if (spatial_timer == timer_number) UpdateSpatial(); rm(spatial_timer)
#   spatial <- append()
    spatial <- sim$spatial
    sim$spatial <- spatial
    return(sim)
  }
}

#' WriteSimList
#'
#' Writes "sim" list to a directory,
#'
#' @usage WriteSimList(write, run, sim, output_dir, components)
#'
#' @param write  = logical, whether or not to write the sim list to a file.
#' Default is TRUE.
#' @param run = name of run. Default uses the name from the runs object to
#' ensure that the proper number of leading zeros is used. If default is not
#' use, the name will simply be: "sim_(run).RData"
#' @param sim = a 'sim' list object
#' @param output_dir = output directory, default is working directory
#' @param components = components of 'sim' to write. Options include: "agents",
#' "pars", and "spatial". Default is all.
#'
#' @return Writes a file to the output_dir
#' @export
#'
#' @note NOT FINISHED - Code not implemented for writing subset of components
WriteSimList <- function(write = TRUE,
                         run = names(runs[j]),
                         sim = sim,
                         output_dir = getwd(),
                         components = "all") {
  sim = sim
  if (write == TRUE) {
    if (components == "agents") {
      file_path = file.path(output_dir, paste0("sim_",run,".RData"))
      print(paste0("Writing: ", file_path, " - agents only"))
      save(sim$agents, file = file_path)
    }
    if (components == "all") {
      file_path = file.path(output_dir, paste0("sim_",run,".RData"))
      print(paste0("Writing: ", file_path))
      save(sim, file = file_path)
    }
  }
}

#------------------------------------------------------------------------------#
################################ OLD CODE ######################################
#------------------------------------------------------------------------------#

# BehaviorSubModelBAEA <- function(sim = sim,
#                                  agent_states = agent_states,
#                                  step_data = step_data,
#                                  step = step){
#   sex <- agent_states$sex
#   beta <- as.matrix(sim$pars$classes[[sex]]$constant$fixed$behavior_betas)
#   step_row <- which(step_data$datetime == step)
#   current_behavior <- as.numeric(step_data[step_row, "behavior"])
#   current_time_prop <- as.numeric(step_data[step_row, "time_proportion"])
#   next_time_prop <- as.numeric(step_data[step_row + 1, "time_proportion"])
#   step_data <- as.data.frame(step_data)
#   gamma <- diag(5)
#   g <- beta[1, ]  #  g = state transition probabilities intercepts
#   g <- g +
#     beta[2, ] * cos(2*pi * (step_data[step_row, "julian"]/365)) +
#     beta[3, ] * sin(2*pi * (step_data[step_row, "julian"]/365)) +
#     beta[4, ] * cos(2*pi * (step_data[step_row, "time_proportion"])) +
#     beta[5, ] * sin(2*pi * (step_data[step_row, "time_proportion"]))
#   exp(g)
#   gamma[!gamma] <- exp(g) # Beta values in non-diagonal locations in matrix
#   gamma2 <- t(gamma) # probabilities for state transitions are now in rows
#   gamma3 <- gamma2/apply(gamma2, 1, sum) # rows sum to 1
#   if(current_time_prop <= .5 & current_behavior != 5){
#     gamma3[, 5] <- 0
#     gamma3 <- gamma3/apply(gamma3, 1, sum)
#     gamma3[, 5] <- 0
#   }
#   # next control is new - it forces agent to leave roost after .4 time prop
#   if(current_time_prop >= .4 & current_time_prop <= .5 & current_behavior == 5){
#     gamma3[, 5] <- 0
#     gamma3 <- gamma3/apply(gamma3, 1, sum)
#     gamma3[, 5] <- 0
#   }
#   if(current_time_prop > .5 & current_behavior == 5){
#     gamma3[, 1:4] <- 0
#     gamma3[, 5] <- 1
#     #gamma3 <- gamma3/apply(gamma3, 1, sum)
#   }
#   if(current_behavior == 1){ # prevents 1->5 (Cruise to Roost)
#     gamma3[, 5] <- 0
#     gamma3 <- gamma3/apply(gamma3, 1, sum)
#     gamma3[, 5] <- 0
#   }
#   if(current_behavior == 5){ # prevents 5->1 (Roost to Cruise)
#     gamma3[, 1] <- 0
#     gamma3 <- gamma3/apply(gamma3, 1, sum)
#     gamma3[, 1] <- 0
#   }
#   # trans prob. given current behavior
#   next_behavior <- sample(1:5, size = 1, prob = gamma3[current_behavior, ])
#
#   if(step_row != nrow(step_data)) { # prevent the creation of an "extra" step
#     if(next_time_prop == 1){ # select Nest or Roost only for last location of day
#       if (next_behavior %in% c(3,5)){
#         step_data[step_row + 1, "behavior"] <- next_behavior
#       } else { # forces selection of Nest or Roost
#         gamma3[, c(1,2,4)] <- 0
#         gamma3 <- gamma3/apply(gamma3, 1, sum)
#         gamma3[, c(1,2,4)] <- 0
#         next_behavior <- sample(1:5, size = 1, prob = gamma3[current_behavior,])
#         step_data[step_row + 1, "behavior"] <- next_behavior
#       }
#     }
#     if(current_time_prop == 1){
#       step_data[step_row + 1, "behavior"] <- current_behavior
#     }
#     step_data[step_row + 1, "behavior"] <- next_behavior
#   }
#   return(step_data)
# }