4 IBM Simulation Loop

if any changes are made, need to re-run qmd_to_r_script() at end, in console

Purpose: Define run_simulation(), the IBM simulation loop wrapper. build_habitat() lives in Habitat.qmd; this file focuses solely on the simulation loop, fish initialisation, and result collation.

Requires Functions.R and Habitat.R to be sourced first.

4.1 Design notes

4.1.1 Key modules

Patch choice (movement)
- Optimal movers select habitat patches (move) based on sensed differences in growth potential. Fish compare growth potential in their current habitat to that in the other habitat, where sensed growth potential is either a function of temperature alone or temperature and density (see Configuration options).
- For all habitat sensing options, growth potential in the alternative habitat (i.e., the growth benefit of moving) is penalized by an energetic cost of movement, which is an allometric function of size (relative growth penalty in g/g/d declines with increasing fish weight). Finally, the decision to move based on differences in (movement penalized) growth potential follows a probabilistic (softmax) decision rule, where fish switch with a probability that increases with the sensed growth advantage.
Growth
- After all fish move, fish densities in each patch are calculated, and the temperature adjusted P_Cmax (physiological capacity for consumption given temperature) is penalized by a density-dependent hyperbolic reduction function.
- A Wisconsin bioenergetics model is used to calculate fish growth based on fish weight, temperature, and density (P_Cmax).
- All fish sizes are update accordingly.
Spawning and reproduction
- The decision to spawn (i.e., maturation probability) is function of day of year (potential spawning dates are normally distributed around a mean date), fish size (spawning probability increases with fish weight following a logistic function), and fish condition (probability increases with condition, calculated as current weight / peak weight).
- Fish can only spawn once per year
- When a fish spawns, they experience an energetic cost that is calculated by multiple current body weight by a fixed proportion (e.g., 0.2). In other words, fish lose 20% of their weight if all spawning conditions are met.
- Fecundity (number of eggs per spawner) is an increasing function of weight (log-log linear relationship) with log-scale noise.
- All eggs are subject to a fixed egg-to-fry survival rate, e.g., 0.1.
- Spawning and reproduction occurs instantaneously, in the same time step. I.e., if spawning conditions are met on some day, their offspring are added to the simulation on that same day.
- Offspring inherit their parents patch and life history strategy (e.g., optimal mover or resident)
Survival
- All fish (including new offspring) are subject to weight, temperature, and condition-based mortality sources
  - Weight: daily survival probability increases with fish weight, up to some maximum probability (0.9994) that is tuned to yield an annual survival probability of ~0.8.
    - High temperature exposure: this mortality source represents the acute breakdown of physiological processes at high temperatures. It does not represent chronic effects of high temperatures on bioenergetics. Daily survival probability decreases above some critical high temperature, following a logistic function.
    - Condition/starvation: fish in poor condition (current weight / peak weight) experience reduced daily survival probability, following a logistic function.
- Daily survival probability is calculated as the joint probability of the sources described above, and daily survival is drawn from a binomial distribution. Fish that die are removed from the simulation, but their results are stored.

4.1.2 Configuration options

Stochasticity in movement is accounted for by one of three alternative mechanisms:
- No stochastic movement/entirely deterministic: all movers do the same thing…they all can perfectly sense and compare growth opportunities
- Probabilistic (softmax) rule: the softmax function assigns probabilities to each patch based on the growth difference, with higher probabilities for patches offering greater growth benefits.
- Individual-level movement threshold: each fish has a unique threshold for movement, representing individual variation in boldness/dispersal propensity
Density-dependent food availaility: after fish have selected patches, calculate food availability based on patch density and baseline ration levels. Two alternative approaches:
- Modified Fullerton approach: ration size decreases linearly with fish density, up to some (a priori defined) threshold above which there is no effect of there is no effect of density of ration size.
- Hyperbolic reduction: describes a non-linear decrease in ration size with increasing fish density, where the rate of decrease slows as density increases.
Environment sensing: how omniscient are fish to temperature and density-dependent growth across the habitat patches?
- abiotic_only: fish can only sense temperature-dependent GP, regardless of which habitat is occupied
- density_current: fish sense density-dependent GP in current habitat, and compare to abiotic (temp-based) GP in the alternative habitat
- density_all: fish can sense density-dependent GP in all habitats (this is how Railsback does it: fish can sense habitat and competition conditions in all cells within a specific radius)

4.1.3 Senstive parameters

under construction

What parameters do we think might have very strong effects on model outcomes?

Habitat variables (temperature regimes and maximum P_C_max)
Strength of density-dependence and differences between habitats ()
Egg-to-fry survival rate (could mediate the trade-off between fish size and number)
Strength of size-based dominance ()

Other parameters:

General
- Habitat variables (temperature regimes and food availability in each patch)
- Number of fish per strategy/total number of starting fish
Density-dependence
- Fullerton upper density threshold
- Hyperbolic half saturation density
Movement decision
- Probabilistic/softmax tau (strength of threshold)
Movement cost
- Allometric function b (metabolic exponent) and c (scaling constant)
Movement stochasticity
- Probabilistic/softmax tau (sensitivity to growth difference)
- Individual-level movement threshold, SD in individual boldness

4.1.4 Misc. notes

Currently, we just look-up growth from the precalculated growth array, but this is going to be less accurate for the IBM b/c precalculated growth rates are calculated over 0.25g mass increments, whereas growth in the IBM is continuous. Probably fine for now, but perhaps should borrow code from above/Fullerton IBM to increase accuracy. My guess is that this will dramatically increase run times.
If sense_environment == "density_all", the density used for the alternative patch is the current pre-movement density (i.e., what the fish can observe before anyone decides to move that day). This is the most tractable approximation — computing post-movement density would require solving a joint decision problem.
If using individual-level movement thresholds (i.e., boldness), residents can be given move_threshold = Inf to formally unify all three strategies under one mechanism, rather than branching by strategy label — keeps the code cleaner if we expand to more strategies later.

4.1.5 To do

eventually, should probably collate all parameters into a single control script/page…that way we can see where all the sensitive/arbitrary parameters are
stochasticity in movement: it probably makes sense to combine stochastic movement processes…e.g, probabilistic and individual thresholds, and perhaps even random noise as well
Size-at-maturity could probably be sampled probabilistically, with noise. Preliminary simulation will just set a specific size threshold.
Same for timing of reproduction: rather than all mature fish reproducing on a single day, we could have them spread out over a given period, or make spawn timing size-dependent, with larger fish spawning first. Lots of empirical evidence for this. That would be very interesting because it would give the offspring of the larger/early spawners a competitive advantage over offspring of smaller/later spawners, likely magnifying the effect of size on reproductive success.

4.2 `run_simulation()`

Run the IBM simulation loop over the full habitat_df time series. Returns a named list containing all primary result objects.

Arguments:

habitat_df — output of build_habitat()
params — named list of simulation parameters (see defaults below)
wt_growth — pre-calculated growth array (load from data/wt.growth.array.RData)

Returns:

Element	Description
`ibm_records`	List of per-day data frames (raw loop output)
`ibm_long`	Long-format tibble (one row per fish × day)
`ibm_summary`	Daily summary statistics by strategy
`fish_registry`	One row per fish (founders + all offspring)
`spawn_log`	One row per spawning event
`switches`	Vector of patch-switch counts indexed by pid
`params`	Parameter list passed in (for provenance)
`habitat_df`	Habitat data frame passed in (for provenance)

Code

run_simulation <- function(habitat_df, params = list(), wt_growth) {

  # ── Parameters ──────────────────────────────────────────────────────────────
  # Configuration options:
  move_stochastic   <- params$move_stochastic   %||% "prob"           # individual movement stochasticity: options = "none", "prob", "indiv"
  food_densdepen    <- params$food_densdepen     %||% "hyperbolic"     # density dependence: options = "fullerton", "hyperbolic"
  sense_environment <- params$sense_environment %||% "density_current" # how omniscient are fish to temperature and density-dependent growth across the habitat patches?
                                                                        #   - "abiotic_only" = Omniscient to abiotic (temp-based) GP only
                                                                        #   - "density_current" = Sense density dependent GP in current habitat, and compare to abiotic GP in the alternative habitat (perhaps the most biologically realistic)
                                                                        #   - "density_all" = Omniscient to density dependent GP in all habitats

  # misc. parameters
  n_fish_per_strat  <- params$n_fish_per_strat  %||% 50        # fish per strategy
  start_wt          <- params$start_wt          %||% 0.5       # initial weight (g)
  wt_sd             <- params$wt_sd             %||% 0.05      # SD of starting weight
  MaxDensity4Growth <- params$MaxDensity4Growth %||% 50        # growth is not depressed further above this density (currently this is arbitrary)
  movecost_c        <- params$movecost_c        %||% 0.1       # movement cost allometric function: `c` = scaling constant
  movecost_b        <- params$movecost_b        %||% 0.3       # movement cost allometric function `b` = decay exponent (higher = steeper drop-off)
  sigma_bold        <- params$sigma_bold        %||% 0.002     # SD of boldness distribution, for dispersal propensity (g/g/d units)
  tau               <- params$tau               %||% 0.003     # sensitivity to growth difference (for probabilistic/softmax variation in movement)
  s_min             <- params$s_min             %||% 0.96      # minimum size-based daily survival probability
  s_w0              <- params$s_w0              %||% 7         # inflection point of size-based survival logistic curve
  s_k               <- params$s_k               %||% 1         # steepness of size-based survival logistic curve
  T1_mort           <- params$T1_mort           %||% 30        # temperature (°C) at which daily survival from thermal stress = 0.1
  T9_mort           <- params$T9_mort           %||% 25.8      # temperature (°C) at which daily survival from thermal stress = 0.9
  K9_starv          <- params$K9_starv          %||% 0.55      # relative condition (W_current/W_peak) at which starvation survival = 0.9
  K1_starv          <- params$K1_starv          %||% 0.45      # relative condition at which starvation survival = 0.1
  age_x0            <- params$age_x0            %||% 14        # logistic inflection point — age (years) at which senescent decline is steepest
  age_k             <- params$age_k             %||% 0.7       # steepness of logistic senescence decline; larger = sharper drop
  age_pmin          <- params$age_pmin          %||% 0.99      # asymptotic daily survival floor for very old fish
  dominance_beta    <- params$dominance_beta    %||% 1         # size-based dominance exponent for effective density (0 = pure scramble, 1 = linear dominance, Inf = pure contest)
  egg_wt            <- params$egg_wt            %||% 0.07      # weight of a single egg (g); energetic cost per offspring
  repro_cost        <- params$repro_cost        %||% 0.2       # energetic cost of reproduction in percentage of body mass
  sigma_fecund      <- params$sigma_fecund      %||% 0.085     # random variation in weight-fecundity relationship
  egg_surv                   <- params$egg_surv                   %||% 0.1    # egg-to-fry (hatching) survival probability ("pass-through" filter for fecundity)
  age_structured_competition <- params$age_structured_competition %||% TRUE   # if TRUE, fry (age-0) compete only with fry; age-1+ compete only with age-1+
  crit_pcmax_lo              <- params$crit_pcmax_lo              %||% 0.30   # pcmax_adjusted_dd at which ~1% critical-period survival is achieved
  crit_pcmax_hi              <- params$crit_pcmax_hi              %||% 0.40   # pcmax_adjusted_dd at which ~99% critical-period survival is achieved (negligible cost)
  crit_period_days           <- params$crit_period_days           %||% 60L    # length of the critical period (days post-hatch)
  sim_seed                   <- params$sim_seed                   %||% 7843

  # Competition grouping vector: used in group_by(across(all_of(eff_grp))) in Steps 2 and 3.
  # When age_structured_competition = TRUE, eff_density is computed within patch × age class.
  # When FALSE, all fish in the same patch compete (original behaviour).
  eff_grp <- if (age_structured_competition) c("patch", "age_class") else "patch"

  # A_warm, A_cold, K_warm, K_cold, S_max_warm, S_max_cold are read from habitat_df each day (see Habitat.qmd)

  # Fullerton has movement cost coded slightly differently, where movement cost is a function of
  # the instantaneous growth rate, distance moved, and a constant...rather than simply a fixed cost
  # as is the case here. But because we only have two patches that are spatially implicit, we can
  # ignore this more complex parameterization...at least for now.

  # ── Local helpers ────────────────────────────────────────────────────────────
  get_patchtemp <- function(t, patch) {
    if (patch == "warm") habitat_df$temp_warm[habitat_df$dayofsim == t]
    else                 habitat_df$temp_cold[habitat_df$dayofsim == t]
  }

  # Lookup sequences — must match the axes of the wt_growth array
  wt_seq_ibm <- seq(0.1,   25,  0.1)    # water temperature (250 values)
  ra_seq_ibm <- seq(0.001, 0.4, 0.001)  # ration (400 values)
  min_wt     <- 0.25                     # absolute minimum fish weight (g)

  # ── Initialise fish ─────────────────────────────────────────────────────────
  set.seed(sim_seed)

  # optimal movers are randomly placed into either warm or cold patches
  fish_pop <- tibble(
    strategy = "optimal_mover",
    patch    = sample(c("warm", "cold"), size = n_fish_per_strat, replace = TRUE),
    weight   = pmax(min_wt, rnorm(n_fish_per_strat, start_wt, wt_sd))
  ) |>
    mutate(
      pid               = row_number(),
      cmax_allometric   = NA_real_,
      pcmax_baseline    = NA_real_,
      pcmax_adjusted    = NA_real_,
      pcmax_adjusted_dd = NA_real_,
      func_temp         = NA_real_,
      ration            = NA_real_,
      # Individual movement threshold — drawn once, fixed for life
      # Positive = reluctant to move (needs a clear advantage)
      # Negative = bold/prone to move (moves even when slightly disadvantaged)
      move_threshold    = rnorm(n_fish_per_strat, mean = 0, sd = sigma_bold),
      peak_weight       = weight,      # track historical max weight for condition factor
      age_days          = 0L,          # age in days, incremented each iteration independent of d
      spawned_this_year = FALSE,       # has this fish already spawned in the current calendar year?
      parent_pid        = NA_integer_, # pid of parent fish (NA for founder cohort)
      cohort            = year(habitat_df$date[1]),  # birth year (cohort)
      prob_surv         = 1,
      survive           = 1
    ) |>
    group_by(patch) |>
    mutate(density = n() / if_else(first(patch) == "warm",
                                   habitat_df$A_warm[1], habitat_df$A_cold[1])) |>
    ungroup()

  # ── Pre-allocate storage ─────────────────────────────────────────────────────
  n_days <- nrow(habitat_df)
  n_fish <- nrow(fish_pop)

  # fish_registry: grows throughout the simulation as offspring are born.
  # Used in place of fish_pop_init in summaries once reproduction is active.
  # One row per fish (founders + all offspring), ordered by pid.
  fish_registry <- fish_pop |>
    select(pid, strategy, parent_pid, cohort) |>
    mutate(birth_dayofsim = 0L)

  # Pre-allocate per-day record list; one data.frame per loop iteration.
  # Avoids large sparse matrices — only alive (fish × day) combos are stored.
  ibm_records <- vector("list", n_days)
  switches    <- integer(n_fish)  # count patch switches per fish (indexed by pid)
  next_pid    <- n_fish + 1L      # next available pid (grows as offspring are born)

  # Spawn log: one row per spawning event, recording parent-offspring relationships
  spawn_log <- tibble(
    parent_pid          = integer(),  # pid of the spawning parent
    dayofsim            = integer(),  # simulation day of spawning
    n_offspring         = integer(),  # number of offspring produced
    weight              = numeric(),  # pre-spawn body weight (g)
    condition           = numeric(),  # pre-spawn relative condition (weight / peak_weight)
    offspring_pid_start = integer(),  # first pid assigned to offspring
    offspring_pid_end   = integer()   # last pid assigned to offspring
  )

  # ── Simulation loop ──────────────────────────────────────────────────────────
  set.seed(sim_seed)

  for (d in seq_len(n_days)) {
    t            <- habitat_df$dayofsim[d]
    current_year <- year(habitat_df$date[d])     # calendar year of this simulation day
    fish_pop$age_days <- fish_pop$age_days + 1L  # increment age for all living fish

    # 1. READ TIME-VARYING HABITAT PARAMETERS ─────────────────────────────────
    A_warm     <- habitat_df$A_warm[d]
    A_cold     <- habitat_df$A_cold[d]
    K_warm     <- habitat_df$K_warm[d]
    K_cold     <- habitat_df$K_cold[d]
    S_max_warm <- habitat_df$S_max_warm[d]
    S_max_cold <- habitat_df$S_max_cold[d]
    pcmax_warm <- habitat_df$pcmax_warm[d]
    pcmax_cold <- habitat_df$pcmax_cold[d]

    # Temperature indices for this day — shared by all fish
    T_warm      <- get_patchtemp(t, "warm")             # get warm patch temp on day t
    T_cold      <- get_patchtemp(t, "cold")             # get cold patch temp on day t
    wt_idx_warm <- which.min(abs(wt_seq_ibm - T_warm))  # water temp index for warm patch on day t
    wt_idx_cold <- which.min(abs(wt_seq_ibm - T_cold))  # water temp index for cold patch on day t

    # If a patch has collapsed to zero area, force all fish out of it before any
    # density calculations (A = 0 would cause 0/0 = NaN in the density step).
    if (A_cold == 0 && any(fish_pop$patch == "cold")) {
      fish_pop$patch[fish_pop$patch == "cold"] <- "warm"
    }
    if (A_warm == 0 && any(fish_pop$patch == "warm")) {
      fish_pop$patch[fish_pop$patch == "warm"] <- "cold"
    }

    # If both patches have zero area, there is no valid habitat: exit immediately.
    # (When A_cold = A_warm = 0, the guard above cycles fish cold→warm→cold, leaving
    # them in cold with A_cold = 0. Any 1-fish age-class group then computes
    # (n-1)/A = 0/0 = NaN, which propagates through eff_density → ration → error.)
    if (A_cold == 0 && A_warm == 0) break

    # Calculate habitat, temperature, and body size dependent rations following convo. with Jonny
    fish_pop <- fish_pop |>
      mutate(
        func_temp        = ifelse(patch == "warm", fncTempDepend(T_warm), fncTempDepend(T_cold)),
        pcmax_baseline   = ifelse(patch == "warm", pcmax_warm, pcmax_cold),
        pcmax_adjusted   = ifelse(func_temp < pcmax_baseline, func_temp, pcmax_baseline),
        cmax_allometric  = fncAllomCmax(weight)
      )

    # 2. CHOOSE PATCH: OPTIMAL MOVERS ─────────────────────────────────────────
    # Mirrors fncGrowthPossible() logic, vectorized across all movers at once.
    # Patch choice senses baseline temperature and ration, but not density.
    mover_rows <- which(fish_pop$strategy == "optimal_mover")  # get row indices for optimal movers

    if (length(mover_rows) > 0) {
      # Mass indices for each mover (clamped to array bounds)
      ma_idx_vec    <- pmax(1L, pmin(4500L, round(fish_pop$weight[mover_rows])))  # get row index for current mass
      ## 2.A.
      in_warm       <- fish_pop$patch[mover_rows] == "warm"  # is the fish in warm patch? T/F
      move_cost_vec <- fncMoveCost_allometric(fish_pop$weight[mover_rows],
                                              c = movecost_c, b = movecost_b)  # calculate size-dependent cost of movement

      # Capture pre-step patch for each mover before any assignments this day.
      # Used by density_all's sequential loop (movement cost direction) and by all
      # modes for switch detection in step 3.
      prev_patch    <- fish_pop$patch[mover_rows]

      # Pre-compute each mover's hypothetical ration in BOTH patches for patch choice sensing.
      # This initial step is done without accounting for the effect of density on p_cmax/ration.
      # pcmax scalars are day-level constants; cmax_allometric is already in fish_pop from step 1.
      pcmax_warm_adj <- min(fncTempDepend(T_warm), pcmax_warm)
      pcmax_cold_adj <- min(fncTempDepend(T_cold), pcmax_cold)
      ra_warm_vec    <- fish_pop$cmax_allometric[mover_rows] * pcmax_warm_adj
      ra_cold_vec    <- fish_pop$cmax_allometric[mover_rows] * pcmax_cold_adj
      ra_idx_warm_vec <- pmax(1L, pmin(400L, map_int(ra_warm_vec, ~which.min(abs(ra_seq_ibm - .x)))))
      ra_idx_cold_vec <- pmax(1L, pmin(400L, map_int(ra_cold_vec, ~which.min(abs(ra_seq_ibm - .x)))))

      if (sense_environment == "abiotic_only") {
        ### 2.A.1. Individuals sense only temperature-based growth potential across all habitats
        g_warm_vec <- wt_growth[cbind(wt_idx_warm, ra_idx_warm_vec, ma_idx_vec)]
        g_cold_vec <- wt_growth[cbind(wt_idx_cold, ra_idx_cold_vec, ma_idx_vec)]
        g_stay     <- ifelse(in_warm, g_warm_vec, g_cold_vec)  # growth if staying in current patch
        g_move_net <- ifelse(in_warm, g_cold_vec - move_cost_vec, g_warm_vec - move_cost_vec)  # growth if moving to alternate patch

      } else if (sense_environment == "density_current") {
        ### 2.A.2. Individuals sense density-dependent GP in current habitat, abiotic GP in alternative habitat
        # Each mover uses its own dominance-weighted effective competitor density (fish per unit
        # area, self excluded) in its current patch.
        # Effective competitor density: grouped by patch × age_class when
        # age_structured_competition = TRUE, by patch only when FALSE.
        eff_density_all      <- fish_pop |>
          mutate(age_class = if_else(cohort == current_year, "age0", "age1plus")) |>
          group_by(across(all_of(eff_grp))) |>
          mutate(
            A_patch     = if_else(first(patch) == "warm", A_warm, A_cold),
            eff_density = (fncEffDensity(weight, beta = dominance_beta) - 1) / A_patch
          ) |>
          ungroup() |>
          pull(eff_density)
        eff_dens_current_vec <- eff_density_all[mover_rows]

        # DD-adjusted ration uses each mover's per-area effective competitor density
        ra_warm_dd_vec     <- fish_pop$cmax_allometric[mover_rows] * pcmax_warm_adj *
                              (K_warm / (K_warm + eff_dens_current_vec))
        ra_cold_dd_vec     <- fish_pop$cmax_allometric[mover_rows] * pcmax_cold_adj *
                              (K_cold / (K_cold + eff_dens_current_vec))
        ra_idx_warm_dd_vec <- pmax(1L, pmin(400L, map_int(ra_warm_dd_vec, ~which.min(abs(ra_seq_ibm - .x)))))
        ra_idx_cold_dd_vec <- pmax(1L, pmin(400L, map_int(ra_cold_dd_vec, ~which.min(abs(ra_seq_ibm - .x)))))

        # Warm fish sense dd-adjusted warm / abiotic cold; cold fish sense abiotic warm / dd-adjusted cold
        g_warm_vec <- wt_growth[cbind(wt_idx_warm,
                                      ifelse(in_warm, ra_idx_warm_dd_vec, ra_idx_warm_vec),
                                      ma_idx_vec)]
        g_cold_vec <- wt_growth[cbind(wt_idx_cold,
                                      ifelse(in_warm, ra_idx_cold_vec, ra_idx_cold_dd_vec),
                                      ma_idx_vec)]
        g_stay     <- ifelse(in_warm, g_warm_vec, g_cold_vec)
        g_move_net <- ifelse(in_warm, g_cold_vec - move_cost_vec, g_warm_vec - move_cost_vec)

      } else if (sense_environment == "density_all") {
        ### 2.A.3. Sequential dominance-ordered patch selection.
        # Fish choose patches in descending body-size order within each competitive age class
        # (largest / most dominant first). Each fish senses density in both patches from the
        # fish already committed this step — not yesterday's pre-movement snapshot. This
        # eliminates the mass-synchronous overcrowding caused by the old pre-movement
        # approximation, and produces biologically realistic distributions: dominant fish
        # claim the thermally optimal patch, subordinates fill in around them.
        #
        # All fish choose anew each day (no incumbency advantage). Movement cost applies
        # only when the chosen patch differs from the fish's patch at the start of this day
        # (captured in prev_patch above).
        #
        # Patch assignments are written directly to fish_pop$patch here; the generic 2.B
        # stochasticity block below is skipped for this sensing mode.

        mover_weights_vec <- fish_pop$weight[mover_rows]
        age_class_all     <- if_else(fish_pop$cohort == current_year, "age0", "age1plus")
        mover_age_class   <- age_class_all[mover_rows]

        # Scalar helper: effective competitor density for one focal fish (weight wi) against a
        # vector of already-placed fish weights (comp_wts), per unit area. Returns Inf when
        # area == 0 (patch inaccessible → growth → 0); returns 0 when no competitors yet.
        eff_dens_scalar <- function(wi, comp_wts, area) {
          if (area == 0)              return(Inf)
          if (length(comp_wts) == 0) return(0)
          sum(ifelse(comp_wts >= wi, 1, (comp_wts / wi)^dominance_beta)) / area
        }

        # Running placement accumulators: weights of fish already assigned this step.
        # Split by age class when age_structured_competition = TRUE (independent pools).
        if (age_structured_competition) {
          placed_warm_age0  <- numeric(0);  placed_cold_age0  <- numeric(0)
          placed_warm_age1p <- numeric(0);  placed_cold_age1p <- numeric(0)
        } else {
          placed_warm_all <- numeric(0);  placed_cold_all <- numeric(0)
        }

        # Processing groups: independent competitive pools, each sorted dominant-first.
        # When age_structured_competition = TRUE, age-0 and age-1+ are separate pools.
        proc_groups <- if (age_structured_competition) {
          list(
            list(idx = which(mover_age_class == "age1plus"), ac = "age1plus"),
            list(idx = which(mover_age_class == "age0"),     ac = "age0")
          )
        } else {
          list(list(idx = seq_along(mover_rows), ac = "all"))
        }

        for (grp in proc_groups) {
          grp_idx <- grp$idx   # positions within mover_rows / mover_weights_vec
          ac      <- grp$ac
          if (length(grp_idx) == 0) next

          # Sort dominant-first within this age class
          grp_idx <- grp_idx[order(mover_weights_vec[grp_idx], decreasing = TRUE)]

          for (ii in grp_idx) {
            fi       <- mover_rows[ii]          # row index in fish_pop
            wi       <- mover_weights_vec[ii]   # fish weight (g)
            cmi      <- fish_pop$cmax_allometric[fi]
            ma_idx_i <- pmax(1L, pmin(4500L, round(wi)))
            prev_i   <- prev_patch[ii]          # patch at start of this day

            # Fetch placement accumulators for this age class
            if (age_structured_competition) {
              pw <- if (ac == "age0") placed_warm_age0 else placed_warm_age1p
              pc <- if (ac == "age0") placed_cold_age0 else placed_cold_age1p
            } else {
              pw <- placed_warm_all
              pc <- placed_cold_all
            }

            # Effective density in each patch given already-placed fish this step
            ed_warm <- eff_dens_scalar(wi, pw, A_warm)
            ed_cold <- eff_dens_scalar(wi, pc, A_cold)

            # Density-adjusted rations in each patch
            ra_warm_i <- cmi * pcmax_warm_adj * (K_warm / (K_warm + ed_warm))
            ra_cold_i <- cmi * pcmax_cold_adj * (K_cold / (K_cold + ed_cold))
            ra_idx_w  <- pmax(1L, pmin(400L, which.min(abs(ra_seq_ibm - ra_warm_i))))
            ra_idx_c  <- pmax(1L, pmin(400L, which.min(abs(ra_seq_ibm - ra_cold_i))))

            # Growth lookup in each patch
            g_warm_i <- wt_growth[wt_idx_warm, ra_idx_w, ma_idx_i]
            g_cold_i <- wt_growth[wt_idx_cold, ra_idx_c, ma_idx_i]

            # Effective growth after movement cost (cost applies only when switching patches).
            # move_cost_vec[ii] is already computed for all movers above.
            move_cost_i <- move_cost_vec[ii]
            gwarm_eff_i <- g_warm_i - if (prev_i == "warm") 0 else move_cost_i
            gcold_eff_i <- g_cold_i - if (prev_i == "cold") 0 else move_cost_i

            # Patch choice — same stochasticity rule as other sensing modes
            chose_warm <- if (move_stochastic == "none") {
              gwarm_eff_i >= gcold_eff_i
            } else if (move_stochastic == "prob") {
              runif(1) < fncMoveSoftmax(gwarm = gwarm_eff_i, gcold = gcold_eff_i, tau = tau)
            } else {  # "indiv": threshold is a per-fish bias on the warm–cold advantage
              (gwarm_eff_i - gcold_eff_i) > fish_pop$move_threshold[fi]
            }

            # Assign patch and update placement accumulator for subsequent fish
            fish_pop$patch[fi] <- if (chose_warm) "warm" else "cold"

            if (age_structured_competition) {
              if (ac == "age0") {
                if (chose_warm) placed_warm_age0  <- c(placed_warm_age0,  wi)
                else            placed_cold_age0  <- c(placed_cold_age0,  wi)
              } else {
                if (chose_warm) placed_warm_age1p <- c(placed_warm_age1p, wi)
                else            placed_cold_age1p <- c(placed_cold_age1p, wi)
              }
            } else {
              if (chose_warm) placed_warm_all <- c(placed_warm_all, wi)
              else            placed_cold_all <- c(placed_cold_all, wi)
            }
          } # end fish loop
        } # end group loop
      }

      ## 2.B. Impose stochasticity in movement
      # density_all handles patch assignment and stochasticity within its sequential loop
      # above (choices written directly to fish_pop$patch). Skip for that sensing mode.
      if (sense_environment != "density_all") {
        if (move_stochastic == "none") {
          ### 2.B.1. Deterministic: All movers do the same thing
          # Update patch based on whether moving is beneficial.
          # 1. If fish is in warm patch and growth potential of moving exceeds growth potential of staying, move to cold.
          # 2. If fish is in cold patch and growth potential of moving exceeds growth potential of staying, move to warm.
          fish_pop$patch[mover_rows] <- case_when(
            in_warm  & g_move_net > g_stay ~ "cold",
            !in_warm & g_move_net > g_stay ~ "warm",
            TRUE                           ~ fish_pop$patch[mover_rows]
          )
        } else if (move_stochastic == "prob") {
          ### 2.B.2. Probabilistic (softmax) rule: Instead of a hard threshold, fish switch with a
          ### probability that increases with the growth advantage. One parameter: τ (sensitivity;
          ### small = nearly deterministic, large = nearly random).
          # Use cost-adjusted net growth rates: staying is free, moving incurs move_cost_vec.
          # g_stay     = growth in current patch (no cost)
          # g_move_net = growth in alternate patch - move_cost_vec
          # Re-map to warm/cold perspective for fncMoveSoftmax:
          gwarm_eff  <- ifelse(in_warm, g_stay, g_move_net)  # effective growth if occupying warm patch
          gcold_eff  <- ifelse(in_warm, g_move_net, g_stay)  # effective growth if occupying cold patch
          p_warm     <- fncMoveSoftmax(gwarm = gwarm_eff, gcold = gcold_eff, tau = tau)  # calculate probability of selecting warm habitat
          choose_warm <- runif(length(mover_rows)) < p_warm  # choose warm if p_warm exceeds a random number between 0 and 1
          fish_pop$patch[mover_rows] <- if_else(choose_warm, "warm", "cold")  # update patch
        } else if (move_stochastic == "indiv") {
          ### 2.B.3. Individual-level movement threshold: each fish has a unique threshold for
          ### movement, representing individual variation in boldness/dispersal propensity
          fish_pop$patch[mover_rows] <- case_when(
            in_warm  & (g_move_net - g_stay) > fish_pop$move_threshold[mover_rows] ~ "cold",
            !in_warm & (g_move_net - g_stay) > fish_pop$move_threshold[mover_rows] ~ "warm",
            TRUE                                                                    ~ prev_patch
          )
        }
      }

    } # end patch choice

    # Recompute func_temp and pcmax_adjusted based on post-choice patch.
    # Required because some fish may have switched patches during step 2;
    # step 1 values reflect the pre-choice patch and would be stale for switchers.
    fish_pop <- fish_pop |>
      mutate(
        func_temp      = ifelse(patch == "warm", fncTempDepend(T_warm), fncTempDepend(T_cold)),
        pcmax_baseline = ifelse(patch == "warm", pcmax_warm, pcmax_cold),
        pcmax_adjusted = ifelse(func_temp < pcmax_baseline, func_temp, pcmax_baseline)
      )

    # 3. UPDATE HABITAT QUALITY (DENSITY-DEPENDENT PCMAX / RATION) ─────────────
    # Re-apply zero-area guard after patch choice: movers may have re-selected a
    # collapsed patch (A = 0) during step 2. Force them back to the surviving patch
    # before density calculations to prevent division-by-zero / NaN rations.
    if (A_cold == 0 && any(fish_pop$patch == "cold")) {
      fish_pop$patch[fish_pop$patch == "cold"] <- "warm"
    }
    if (A_warm == 0 && any(fish_pop$patch == "warm")) {
      fish_pop$patch[fish_pop$patch == "warm"] <- "cold"
    }

    # Tally switches and define `switched` AFTER zero-area guards: fish forced
    # back to their original patch by a guard are not counted as having switched
    # and incur no movement cost. This prevents spurious energy drain in scenarios
    # where one patch has A = 0 and the other patch looks more attractive (fish
    # "choose" the zero-area patch every step but get reversed by the guard).
    switched <- logical(length(mover_rows))  # default FALSE (no movers)
    if (length(mover_rows) > 0) {
      switched <- fish_pop$patch[mover_rows] != prev_patch
      switches[fish_pop$pid[mover_rows]] <- switches[fish_pop$pid[mover_rows]] +
                                            as.integer(switched)
    }

    # Update density (fish per unit area) and dominance-weighted effective competitor density
    # after fish have moved. density counts all fish in the patch regardless of age class.
    # eff_density is computed within the grouping defined by eff_grp:
    #   age_structured_competition = TRUE  → group by patch × age_class (fry vs. age-1+)
    #   age_structured_competition = FALSE → group by patch only (all fish compete together)
    # This reflects ontogenetic habitat segregation in salmonids when enabled. Self is
    # excluded via the (-1) correction from fncEffDensity.
    fish_pop <- fish_pop |>
      mutate(age_class = if_else(cohort == current_year, "age0", "age1plus")) |>
      group_by(patch) |>
      mutate(
        A_patch = if_else(first(patch) == "warm", A_warm, A_cold),
        density = n() / A_patch
      ) |>
      group_by(across(all_of(eff_grp))) |>
      mutate(
        eff_density = (fncEffDensity(weight, beta = dominance_beta) - 1) / A_patch
      ) |>
      ungroup() |>
      select(-A_patch, -age_class)

    if (food_densdepen == "fullerton") {
      ### 3a. Fullerton approach (uses raw patch density; dominance not yet integrated here)
      fdens <- fish_pop$density
      fdens[fdens > MaxDensity4Growth] <- MaxDensity4Growth
      density_effect <- fncRescale((1 - c(fdens, 0.01, MaxDensity4Growth)), c(0.5, 1))
      density_effect <- density_effect[-c(length(density_effect), length(density_effect) - 1)]
      fish_pop <- fish_pop |> mutate(ration = ration * density_effect)
    } else if (food_densdepen == "hyperbolic") {
      ### 3b. Hyperbolic reduction — each fish's P_Cmax is penalised by its own effective density
      # (eff_density, in fish per unit area, self excluded). Large dominant fish face a lower
      # eff_density → less suppression; small subordinate fish face a higher eff_density → more.
      # No -1 correction needed: self is already excluded from eff_density.
      fish_pop <- fish_pop |>
        mutate(
          k                 = if_else(patch == "warm", K_warm, K_cold),
          pcmax_adjusted_dd = pcmax_adjusted * (k / (k + eff_density)),
          ration            = cmax_allometric * pcmax_adjusted_dd
        ) |>
        select(-k)
    }

    # 4. GROW FISH ─────────────────────────────────────────────────────────────
    # Growth lookup for all fish via fncGrowthFish. Growth in g/g/d.
    growth_df <- fncGrowthFish(NA, as.data.frame(fish_pop), T_warm, T_cold)

    # Apply movement cost: deduct from growth rate for fish that actually switched patches.
    # move_cost_vec is in g/g/d and is indexed over mover_rows; switched is the logical
    # subset of those rows that changed patch this step.
    if (length(mover_rows) > 0 && any(switched)) {
      switched_rows <- mover_rows[switched]
      growth_df$growth[switched_rows] <- growth_df$growth[switched_rows] - move_cost_vec[switched]
    }

    # Update weight: W_{t+1} = W_t + (growth_rate * W_t)
    fish_pop$weight      <- fish_pop$weight + (growth_df$growth * fish_pop$weight)
    # Update peak weight: ratchet up only, never down
    fish_pop$peak_weight <- pmax(fish_pop$peak_weight, fish_pop$weight)

    # 5. SPAWNING AND REPRODUCTION ─────────────────────────────────────────────
    # Reset spawned_this_year flag at the calendar year boundary
    if (d > 1) {
      prev_year <- year(habitat_df$date[d - 1])
      if (current_year != prev_year) fish_pop$spawned_this_year <- FALSE
    }

    # Day-of-year and relative condition needed for spawning probability
    doy_d         <- yday(habitat_df$date[d])
    condition_spw <- fish_pop$weight / fish_pop$peak_weight

    # Combined daily spawning probability: size × condition × date
    p_spawn <- fncMaturitySize(fish_pop$weight) *
               fncMaturityCondition(condition_spw) *
               fncMaturityDate(doy_d)

    # Fish that have already spawned this year are ineligible
    p_spawn[fish_pop$spawned_this_year] <- 0

    # Bernoulli trial: which fish spawn today?
    spawns <- as.logical(rbinom(nrow(fish_pop), size = 1, prob = p_spawn))

    # For spawning fish: calculate fecundity, impose energetic cost, mark as spawned
    if (any(spawns)) {
      spawner_idx <- which(spawns)

      # Fecundity: number of eggs as a function of weight, with log-scale noise
      n_eggs      <- round(fncFecundBromage(fish_pop$weight[spawner_idx],
                                            sigma = sigma_fecund, survival = egg_surv))

      # Capture pre-cost weight for logging (spawn decision was made at this weight)
      weight_at_spawn  <- fish_pop$weight[spawner_idx]

      # Energetic cost of reproduction: deduct from parent weight
      # Weight floor at start_wt — spawning cannot reduce a fish below hatch weight
      reproduction_cost <- weight_at_spawn * repro_cost
      fish_pop$weight[spawner_idx]          <- pmax(start_wt, weight_at_spawn - reproduction_cost)

      # Flag these fish as having spawned; they are ineligible to spawn again this year
      fish_pop$spawned_this_year[spawner_idx] <- TRUE

      # Create offspring rows and add to population --------------------------------
      new_fish_list <- vector("list", length(spawner_idx))

      for (i in seq_along(spawner_idx)) {
        si      <- spawner_idx[i]
        parent  <- fish_pop[si, ]
        parent$weight <- weight_at_spawn[i]  # restore pre-cost weight for logging
        n_off   <- n_eggs[i]
        new_pids <- seq(next_pid, next_pid + n_off - 1L)

        # Initial offspring weights drawn from same distribution as founders
        off_wt <- pmax(min_wt, rnorm(n_off, start_wt, wt_sd))

        # Offspring inherit strategy and patch from parent; movement threshold drawn fresh
        off_thresh <- if (parent$strategy == "optimal_mover") {
          rnorm(n_off, mean = 0, sd = sigma_bold)
        } else {
          rep(Inf, n_off)
        }

        new_fish_list[[i]] <- tibble(
          strategy          = parent$strategy,
          patch             = parent$patch,
          weight            = off_wt,
          pid               = new_pids,
          cmax_allometric   = NA_real_,
          pcmax_baseline    = NA_real_,
          pcmax_adjusted    = NA_real_,
          pcmax_adjusted_dd = NA_real_,
          func_temp         = NA_real_,
          ration            = NA_real_,
          move_threshold    = off_thresh,
          peak_weight       = off_wt,
          age_days          = 0L,
          spawned_this_year = FALSE,
          parent_pid        = parent$pid,
          cohort            = current_year,
          prob_surv         = 1,
          survive           = 1
        )

        # Log the spawning event (weight and condition are pre-spawn values)
        spawn_log <- bind_rows(spawn_log, tibble(
          parent_pid          = parent$pid,
          dayofsim            = d,
          n_offspring         = n_off,
          weight              = parent$weight,
          condition           = parent$weight / parent$peak_weight,
          offspring_pid_start = next_pid,
          offspring_pid_end   = next_pid + n_off - 1L
        ))

        next_pid <- next_pid + n_off
      }

      # Bind offspring and add them to the population BEFORE survival.
      # Size-based and density-driven starvation mortality will thin the cohort naturally.
      new_fish <- bind_rows(new_fish_list)
      n_new    <- nrow(new_fish)
      switches <- c(switches, integer(n_new))  # extend switch counter for new fish

      # Extend growth_df with stub rows for offspring so survival indexing stays aligned.
      # WT.actual is set to their patch temperature; growth is NA (they didn't feed today).
      off_wt_actual <- ifelse(new_fish$patch == "warm", T_warm, T_cold)
      growth_df <- bind_rows(growth_df,
                             data.frame(WT.actual = off_wt_actual, growth = NA_real_))

      # Register and add offspring to live population
      fish_registry <- bind_rows(
        fish_registry,
        new_fish |>
          select(pid, strategy, parent_pid, cohort) |>
          mutate(birth_dayofsim = d)
      )
      fish_pop <- bind_rows(fish_pop, new_fish)
    }

    # 6. SURVIVAL ─────────────────────────────────────────────────────────────
    s_max_vec <- ifelse(fish_pop$patch == "warm", S_max_warm, S_max_cold)

    # Size-dependent survival: logistic sigmoid in weight, rising from minprob floor to habitat-specific maxprob
    p_sg      <- fncSurviveSize(fish_pop$weight, minprob = s_min, maxprob = s_max_vec, w0 = s_w0, k = s_k)[[1]]

    # Temperature-dependent survival: multiplied with size/growth probability so both act simultaneously
    p_temp    <- fncSurviveTemp(growth_df$WT.actual, T1 = T1_mort, T9 = T9_mort)

    # Condition-based (starvation) survival: relative condition = current / peak weight
    condition <- fish_pop$weight / fish_pop$peak_weight
    p_starv   <- fncSurviveStarve(condition, K9 = K9_starv, K1 = K1_starv)

    # Age-based survival (senescence): survival probability declines past age_thresh years old.
    p_age <- fncSurviveAge(fish_pop$age_days / 365,
                           x0    = age_x0,
                           k     = age_k,
                           p_min = age_pmin)

    # Critical-period consumption-based survival (Elliott 1989) — disabled.
    # Replaced by the sigmoid shape of fncSurviveSize, which sustains strong size-dependent
    # mortality below w0 (~7g) and achieves the same growth-mediated density dependence
    # without requiring a separate consumption threshold module.
    # p_crit <- fncSurviveConsumption(
    #   pcmax_dd         = fish_pop$pcmax_adjusted_dd,
    #   age_days         = fish_pop$age_days,
    #   crit_pcmax_lo    = crit_pcmax_lo,
    #   crit_pcmax_hi    = crit_pcmax_hi,
    #   crit_period_days = crit_period_days
    # )

    # Calculate combined daily survival rate
    prb.srv   <- pmin(p_sg * p_temp * p_starv * p_age, 1)

    # Minimum weight floor: fish that drop below hatch weight die immediately (backstop)
    prb.srv[fish_pop$weight < start_wt] <- 0

    survivors <- rbinom(nrow(growth_df), size = 1, prob = prb.srv)
    growth_df <- growth_df |>
      mutate(prob_surv = prb.srv, survive = survivors)

    # 7. STORE RESULTS AND REMOVE NON-SURVIVORS ────────────────────────────────
    ibm_records[[d]] <- data.frame(
      pid       = fish_pop$pid,
      dayofsim  = d,
      strategy  = fish_pop$strategy,
      weight    = fish_pop$weight,
      patch     = fish_pop$patch,
      temp      = growth_df$WT.actual,   # includes offspring stub (their patch temp)
      ggd       = growth_df$growth,      # NA for offspring on birth day
      survived  = growth_df$survive,
      condition = condition,
      age       = fish_pop$age_days / 365,
      p_survive = prb.srv                # worth keeping — used in later analyses
    )

    # Remove non-survivors — only living fish carry forward to next iteration
    fish_pop <- fish_pop[growth_df$survive == 1, ]

    # exit loop if all fish have died
    if (nrow(fish_pop) == 0) break
  }

  # ── Collate ibm_long ─────────────────────────────────────────────────────────
  ibm_long <- bind_rows(ibm_records) |>
    left_join(fish_registry |> select(pid, parent_pid, cohort, birth_dayofsim), by = "pid") |>
    left_join(habitat_df |> select(dayofsim, date), by = "dayofsim")

  # ── Standard daily summary by strategy ───────────────────────────────────────
  ibm_summary <- ibm_long |>
    group_by(strategy, dayofsim, date) |>
    summarise(
      mean_weight    = mean(weight,    na.rm = TRUE),
      sd_weight      = sd(weight,      na.rm = TRUE),
      mean_temp      = mean(temp,      na.rm = TRUE),
      sd_temp        = sd(temp,        na.rm = TRUE),
      mean_ggd       = mean(ggd,       na.rm = TRUE),
      sd_ggd         = sd(ggd,         na.rm = TRUE),
      prop_warm      = mean(patch == "warm", na.rm = TRUE),
      n_alive        = sum(survived,   na.rm = TRUE),
      mean_condition = mean(condition, na.rm = TRUE),
      sd_condition   = sd(condition,   na.rm = TRUE),
      .groups        = "drop"
    )

  list(
    ibm_records   = ibm_records,
    ibm_long      = ibm_long,
    ibm_summary   = ibm_summary,
    fish_registry = fish_registry,
    spawn_log     = spawn_log,
    switches      = switches,
    params        = params,
    habitat_df    = habitat_df
  )
}

4.3 Write R file

For sourcing by other scripts.

Code

file.remove("SimulationLoop.R")
quarto::qmd_to_r_script("SimulationLoop.qmd")

--- title: "IBM Simulation Loop" --- ```{r, include = FALSE} library(tidyverse) library(lubridate) ``` ***if any changes are made, need to re-run `qmd_to_r_script()` at end, in console*** **Purpose:** Define `run_simulation()`, the IBM simulation loop wrapper. `build_habitat()` lives in `Habitat.qmd`; this file focuses solely on the simulation loop, fish initialisation, and result collation. Requires `Functions.R` and `Habitat.R` to be sourced first. ## Design notes ### Key modules * **Patch choice (movement)** + Optimal movers select habitat patches (move) based on sensed differences in growth potential. Fish compare growth potential in their current habitat to that in the other habitat, where sensed growth potential is either a function of temperature alone or temperature and density (see *Configuration options*). + For all habitat sensing options, growth potential in the alternative habitat (i.e., the growth benefit of moving) is penalized by an energetic cost of movement, which is an allometric function of size (relative growth penalty in g/g/d declines with increasing fish weight). Finally, the decision to move based on differences in (movement penalized) growth potential follows a probabilistic (softmax) decision rule, where fish switch with a probability that increases with the sensed growth advantage. * **Growth** + After all fish move, fish densities in each patch are calculated, and the temperature adjusted P_Cmax (physiological capacity for consumption given temperature) is penalized by a density-dependent hyperbolic reduction function. + A Wisconsin bioenergetics model is used to calculate fish growth based on fish weight, temperature, and density (P_Cmax). + All fish sizes are update accordingly. * **Spawning and reproduction** + The decision to spawn (i.e., maturation probability) is function of day of year (potential spawning dates are normally distributed around a mean date), fish size (spawning probability increases with fish weight following a logistic function), and fish condition (probability increases with condition, calculated as current weight / peak weight). + Fish can only spawn once per year + When a fish spawns, they experience an energetic cost that is calculated by multiple current body weight by a fixed proportion (e.g., 0.2). In other words, fish lose 20% of their weight if all spawning conditions are met. + Fecundity (number of eggs per spawner) is an increasing function of weight (log-log linear relationship) with log-scale noise. + All eggs are subject to a fixed egg-to-fry survival rate, e.g., 0.1. + Spawning and reproduction occurs instantaneously, in the same time step. I.e., if spawning conditions are met on some day, their offspring are added to the simulation on that same day. + Offspring inherit their parents patch and life history strategy (e.g., optimal mover or resident) * **Survival** + All fish (including new offspring) are subject to weight, temperature, and condition-based mortality sources - Weight: daily survival probability increases with fish weight, up to some maximum probability (0.9994) that is tuned to yield an annual survival probability of ~0.8. - High temperature exposure: this mortality source represents the acute breakdown of physiological processes at high temperatures. It does not represent chronic effects of high temperatures on bioenergetics. Daily survival probability decreases above some critical high temperature, following a logistic function. - Condition/starvation: fish in poor condition (current weight / peak weight) experience reduced daily survival probability, following a logistic function. + Daily survival probability is calculated as the joint probability of the sources described above, and daily survival is drawn from a binomial distribution. Fish that die are removed from the simulation, but their results are stored. ### Configuration options * **Stochasticity in movement** is accounted for by one of three alternative mechanisms: + *No stochastic movement/entirely deterministic:* all movers do the same thing...they all can perfectly sense and compare growth opportunities + *Probabilistic (softmax) rule:* the softmax function assigns probabilities to each patch based on the growth difference, with higher probabilities for patches offering greater growth benefits. + *Individual-level movement threshold:* each fish has a unique threshold for movement, representing individual variation in boldness/dispersal propensity * **Density-dependent food availaility:** after fish have selected patches, calculate food availability based on patch density and baseline ration levels. Two alternative approaches: + *Modified Fullerton approach:* ration size decreases linearly with fish density, up to some (*a priori* defined) threshold above which there is no effect of there is no effect of density of ration size. + *Hyperbolic reduction:* describes a non-linear decrease in ration size with increasing fish density, where the rate of decrease slows as density increases. * **Environment sensing:** how omniscient are fish to temperature and density-dependent growth across the habitat patches? + `abiotic_only`: fish can only sense temperature-dependent GP, regardless of which habitat is occupied + `density_current`: fish sense density-dependent GP in current habitat, and compare to abiotic (temp-based) GP in the alternative habitat + `density_all`: fish can sense density-dependent GP in all habitats (this is how Railsback does it: fish can sense habitat and competition conditions in all cells within a specific radius) ### Senstive parameters *under construction* What parameters do we think might have very strong effects on model outcomes? * Habitat variables (temperature regimes and maximum P_C~max~) * Strength of density-dependence and differences between habitats () * Egg-to-fry survival rate (could mediate the trade-off between fish size and number) * Strength of size-based dominance () Other parameters: * General + Habitat variables (temperature regimes and food availability in each patch) + Number of fish per strategy/total number of starting fish * Density-dependence + Fullerton upper density threshold + Hyperbolic half saturation density * Movement decision + Probabilistic/softmax tau (strength of threshold) * Movement cost + Allometric function b (metabolic exponent) and c (scaling constant) * Movement stochasticity + Probabilistic/softmax tau (sensitivity to growth difference) + Individual-level movement threshold, SD in individual boldness ### Misc. notes * Currently, we just look-up growth from the precalculated growth array, but this is going to be less accurate for the IBM b/c precalculated growth rates are calculated over 0.25g mass increments, whereas growth in the IBM is continuous. Probably fine for now, but perhaps should borrow code from above/Fullerton IBM to increase accuracy. My guess is that this will dramatically increase run times. * If `sense_environment == "density_all"`, the density used for the alternative patch is the current pre-movement density (i.e., what the fish can observe before anyone decides to move that day). This is the most tractable approximation — computing post-movement density would require solving a joint decision problem. * If using individual-level movement thresholds (i.e., boldness), residents can be given `move_threshold = Inf` to formally unify all three strategies under one mechanism, rather than branching by strategy label — keeps the code cleaner if we expand to more strategies later. ### To do * eventually, should probably collate all parameters into a single control script/page...that way we can see where all the sensitive/arbitrary parameters are * stochasticity in movement: it probably makes sense to combine stochastic movement processes...e.g, probabilistic and individual thresholds, and perhaps even random noise as well * Size-at-maturity could probably be sampled probabilistically, with noise. Preliminary simulation will just set a specific size threshold. * Same for timing of reproduction: rather than all mature fish reproducing on a single day, we could have them spread out over a given period, or make spawn timing size-dependent, with larger fish spawning first. Lots of empirical evidence for this. That would be very interesting because it would give the offspring of the larger/early spawners a competitive advantage over offspring of smaller/later spawners, likely magnifying the effect of size on reproductive success. ## `run_simulation()` Run the IBM simulation loop over the full `habitat_df` time series. Returns a named list containing all primary result objects. **Arguments:** * `habitat_df` — output of `build_habitat()` * `params` — named list of simulation parameters (see defaults below) * `wt_growth` — pre-calculated growth array (load from `data/wt.growth.array.RData`) **Returns:** | Element | Description | |---|---| | `ibm_records` | List of per-day data frames (raw loop output) | | `ibm_long` | Long-format tibble (one row per fish × day) | | `ibm_summary` | Daily summary statistics by strategy | | `fish_registry` | One row per fish (founders + all offspring) | | `spawn_log` | One row per spawning event | | `switches` | Vector of patch-switch counts indexed by pid | | `params` | Parameter list passed in (for provenance) | | `habitat_df` | Habitat data frame passed in (for provenance) | ```{r} run_simulation <- function(habitat_df, params = list(), wt_growth) { # ── Parameters ────────────────────────────────────────────────────────────── # Configuration options: move_stochastic <- params$move_stochastic %||% "prob" # individual movement stochasticity: options = "none", "prob", "indiv" food_densdepen <- params$food_densdepen %||% "hyperbolic" # density dependence: options = "fullerton", "hyperbolic" sense_environment <- params$sense_environment %||% "density_current" # how omniscient are fish to temperature and density-dependent growth across the habitat patches? # - "abiotic_only" = Omniscient to abiotic (temp-based) GP only # - "density_current" = Sense density dependent GP in current habitat, and compare to abiotic GP in the alternative habitat (perhaps the most biologically realistic) # - "density_all" = Omniscient to density dependent GP in all habitats # misc. parameters n_fish_per_strat <- params$n_fish_per_strat %||% 50 # fish per strategy start_wt <- params$start_wt %||% 0.5 # initial weight (g) wt_sd <- params$wt_sd %||% 0.05 # SD of starting weight MaxDensity4Growth <- params$MaxDensity4Growth %||% 50 # growth is not depressed further above this density (currently this is arbitrary) movecost_c <- params$movecost_c %||% 0.1 # movement cost allometric function: `c` = scaling constant movecost_b <- params$movecost_b %||% 0.3 # movement cost allometric function `b` = decay exponent (higher = steeper drop-off) sigma_bold <- params$sigma_bold %||% 0.002 # SD of boldness distribution, for dispersal propensity (g/g/d units) tau <- params$tau %||% 0.003 # sensitivity to growth difference (for probabilistic/softmax variation in movement) s_min <- params$s_min %||% 0.96 # minimum size-based daily survival probability s_w0 <- params$s_w0 %||% 7 # inflection point of size-based survival logistic curve s_k <- params$s_k %||% 1 # steepness of size-based survival logistic curve T1_mort <- params$T1_mort %||% 30 # temperature (°C) at which daily survival from thermal stress = 0.1 T9_mort <- params$T9_mort %||% 25.8 # temperature (°C) at which daily survival from thermal stress = 0.9 K9_starv <- params$K9_starv %||% 0.55 # relative condition (W_current/W_peak) at which starvation survival = 0.9 K1_starv <- params$K1_starv %||% 0.45 # relative condition at which starvation survival = 0.1 age_x0 <- params$age_x0 %||% 14 # logistic inflection point — age (years) at which senescent decline is steepest age_k <- params$age_k %||% 0.7 # steepness of logistic senescence decline; larger = sharper drop age_pmin <- params$age_pmin %||% 0.99 # asymptotic daily survival floor for very old fish dominance_beta <- params$dominance_beta %||% 1 # size-based dominance exponent for effective density (0 = pure scramble, 1 = linear dominance, Inf = pure contest) egg_wt <- params$egg_wt %||% 0.07 # weight of a single egg (g); energetic cost per offspring repro_cost <- params$repro_cost %||% 0.2 # energetic cost of reproduction in percentage of body mass sigma_fecund <- params$sigma_fecund %||% 0.085 # random variation in weight-fecundity relationship egg_surv <- params$egg_surv %||% 0.1 # egg-to-fry (hatching) survival probability ("pass-through" filter for fecundity) age_structured_competition <- params$age_structured_competition %||% TRUE # if TRUE, fry (age-0) compete only with fry; age-1+ compete only with age-1+ crit_pcmax_lo <- params$crit_pcmax_lo %||% 0.30 # pcmax_adjusted_dd at which ~1% critical-period survival is achieved crit_pcmax_hi <- params$crit_pcmax_hi %||% 0.40 # pcmax_adjusted_dd at which ~99% critical-period survival is achieved (negligible cost) crit_period_days <- params$crit_period_days %||% 60L # length of the critical period (days post-hatch) sim_seed <- params$sim_seed %||% 7843 # Competition grouping vector: used in group_by(across(all_of(eff_grp))) in Steps 2 and 3. # When age_structured_competition = TRUE, eff_density is computed within patch × age class. # When FALSE, all fish in the same patch compete (original behaviour). eff_grp <- if (age_structured_competition) c("patch", "age_class") else "patch" # A_warm, A_cold, K_warm, K_cold, S_max_warm, S_max_cold are read from habitat_df each day (see Habitat.qmd) # Fullerton has movement cost coded slightly differently, where movement cost is a function of # the instantaneous growth rate, distance moved, and a constant...rather than simply a fixed cost # as is the case here. But because we only have two patches that are spatially implicit, we can # ignore this more complex parameterization...at least for now. # ── Local helpers ──────────────────────────────────────────────────────────── get_patchtemp <- function(t, patch) { if (patch == "warm") habitat_df$temp_warm[habitat_df$dayofsim == t] else habitat_df$temp_cold[habitat_df$dayofsim == t] } # Lookup sequences — must match the axes of the wt_growth array wt_seq_ibm <- seq(0.1, 25, 0.1) # water temperature (250 values) ra_seq_ibm <- seq(0.001, 0.4, 0.001) # ration (400 values) min_wt <- 0.25 # absolute minimum fish weight (g) # ── Initialise fish ───────────────────────────────────────────────────────── set.seed(sim_seed) # optimal movers are randomly placed into either warm or cold patches fish_pop <- tibble( strategy = "optimal_mover", patch = sample(c("warm", "cold"), size = n_fish_per_strat, replace = TRUE), weight = pmax(min_wt, rnorm(n_fish_per_strat, start_wt, wt_sd)) ) |> mutate( pid = row_number(), cmax_allometric = NA_real_, pcmax_baseline = NA_real_, pcmax_adjusted = NA_real_, pcmax_adjusted_dd = NA_real_, func_temp = NA_real_, ration = NA_real_, # Individual movement threshold — drawn once, fixed for life # Positive = reluctant to move (needs a clear advantage) # Negative = bold/prone to move (moves even when slightly disadvantaged) move_threshold = rnorm(n_fish_per_strat, mean = 0, sd = sigma_bold), peak_weight = weight, # track historical max weight for condition factor age_days = 0L, # age in days, incremented each iteration independent of d spawned_this_year = FALSE, # has this fish already spawned in the current calendar year? parent_pid = NA_integer_, # pid of parent fish (NA for founder cohort) cohort = year(habitat_df$date[1]), # birth year (cohort) prob_surv = 1, survive = 1 ) |> group_by(patch) |> mutate(density = n() / if_else(first(patch) == "warm", habitat_df$A_warm[1], habitat_df$A_cold[1])) |> ungroup() # ── Pre-allocate storage ───────────────────────────────────────────────────── n_days <- nrow(habitat_df) n_fish <- nrow(fish_pop) # fish_registry: grows throughout the simulation as offspring are born. # Used in place of fish_pop_init in summaries once reproduction is active. # One row per fish (founders + all offspring), ordered by pid. fish_registry <- fish_pop |> select(pid, strategy, parent_pid, cohort) |> mutate(birth_dayofsim = 0L) # Pre-allocate per-day record list; one data.frame per loop iteration. # Avoids large sparse matrices — only alive (fish × day) combos are stored. ibm_records <- vector("list", n_days) switches <- integer(n_fish) # count patch switches per fish (indexed by pid) next_pid <- n_fish + 1L # next available pid (grows as offspring are born) # Spawn log: one row per spawning event, recording parent-offspring relationships spawn_log <- tibble( parent_pid = integer(), # pid of the spawning parent dayofsim = integer(), # simulation day of spawning n_offspring = integer(), # number of offspring produced weight = numeric(), # pre-spawn body weight (g) condition = numeric(), # pre-spawn relative condition (weight / peak_weight) offspring_pid_start = integer(), # first pid assigned to offspring offspring_pid_end = integer() # last pid assigned to offspring ) # ── Simulation loop ────────────────────────────────────────────────────────── set.seed(sim_seed) for (d in seq_len(n_days)) { t <- habitat_df$dayofsim[d] current_year <- year(habitat_df$date[d]) # calendar year of this simulation day fish_pop$age_days <- fish_pop$age_days + 1L # increment age for all living fish # 1. READ TIME-VARYING HABITAT PARAMETERS ───────────────────────────────── A_warm <- habitat_df$A_warm[d] A_cold <- habitat_df$A_cold[d] K_warm <- habitat_df$K_warm[d] K_cold <- habitat_df$K_cold[d] S_max_warm <- habitat_df$S_max_warm[d] S_max_cold <- habitat_df$S_max_cold[d] pcmax_warm <- habitat_df$pcmax_warm[d] pcmax_cold <- habitat_df$pcmax_cold[d] # Temperature indices for this day — shared by all fish T_warm <- get_patchtemp(t, "warm") # get warm patch temp on day t T_cold <- get_patchtemp(t, "cold") # get cold patch temp on day t wt_idx_warm <- which.min(abs(wt_seq_ibm - T_warm)) # water temp index for warm patch on day t wt_idx_cold <- which.min(abs(wt_seq_ibm - T_cold)) # water temp index for cold patch on day t # If a patch has collapsed to zero area, force all fish out of it before any # density calculations (A = 0 would cause 0/0 = NaN in the density step). if (A_cold == 0 && any(fish_pop$patch == "cold")) { fish_pop$patch[fish_pop$patch == "cold"] <- "warm" } if (A_warm == 0 && any(fish_pop$patch == "warm")) { fish_pop$patch[fish_pop$patch == "warm"] <- "cold" } # If both patches have zero area, there is no valid habitat: exit immediately. # (When A_cold = A_warm = 0, the guard above cycles fish cold→warm→cold, leaving # them in cold with A_cold = 0. Any 1-fish age-class group then computes # (n-1)/A = 0/0 = NaN, which propagates through eff_density → ration → error.) if (A_cold == 0 && A_warm == 0) break # Calculate habitat, temperature, and body size dependent rations following convo. with Jonny fish_pop <- fish_pop |> mutate( func_temp = ifelse(patch == "warm", fncTempDepend(T_warm), fncTempDepend(T_cold)), pcmax_baseline = ifelse(patch == "warm", pcmax_warm, pcmax_cold), pcmax_adjusted = ifelse(func_temp < pcmax_baseline, func_temp, pcmax_baseline), cmax_allometric = fncAllomCmax(weight) ) # 2. CHOOSE PATCH: OPTIMAL MOVERS ───────────────────────────────────────── # Mirrors fncGrowthPossible() logic, vectorized across all movers at once. # Patch choice senses baseline temperature and ration, but not density. mover_rows <- which(fish_pop$strategy == "optimal_mover") # get row indices for optimal movers if (length(mover_rows) > 0) { # Mass indices for each mover (clamped to array bounds) ma_idx_vec <- pmax(1L, pmin(4500L, round(fish_pop$weight[mover_rows]))) # get row index for current mass ## 2.A. in_warm <- fish_pop$patch[mover_rows] == "warm" # is the fish in warm patch? T/F move_cost_vec <- fncMoveCost_allometric(fish_pop$weight[mover_rows], c = movecost_c, b = movecost_b) # calculate size-dependent cost of movement # Capture pre-step patch for each mover before any assignments this day. # Used by density_all's sequential loop (movement cost direction) and by all # modes for switch detection in step 3. prev_patch <- fish_pop$patch[mover_rows] # Pre-compute each mover's hypothetical ration in BOTH patches for patch choice sensing. # This initial step is done without accounting for the effect of density on p_cmax/ration. # pcmax scalars are day-level constants; cmax_allometric is already in fish_pop from step 1. pcmax_warm_adj <- min(fncTempDepend(T_warm), pcmax_warm) pcmax_cold_adj <- min(fncTempDepend(T_cold), pcmax_cold) ra_warm_vec <- fish_pop$cmax_allometric[mover_rows] * pcmax_warm_adj ra_cold_vec <- fish_pop$cmax_allometric[mover_rows] * pcmax_cold_adj ra_idx_warm_vec <- pmax(1L, pmin(400L, map_int(ra_warm_vec, ~which.min(abs(ra_seq_ibm - .x))))) ra_idx_cold_vec <- pmax(1L, pmin(400L, map_int(ra_cold_vec, ~which.min(abs(ra_seq_ibm - .x))))) if (sense_environment == "abiotic_only") { ### 2.A.1. Individuals sense only temperature-based growth potential across all habitats g_warm_vec <- wt_growth[cbind(wt_idx_warm, ra_idx_warm_vec, ma_idx_vec)] g_cold_vec <- wt_growth[cbind(wt_idx_cold, ra_idx_cold_vec, ma_idx_vec)] g_stay <- ifelse(in_warm, g_warm_vec, g_cold_vec) # growth if staying in current patch g_move_net <- ifelse(in_warm, g_cold_vec - move_cost_vec, g_warm_vec - move_cost_vec) # growth if moving to alternate patch } else if (sense_environment == "density_current") { ### 2.A.2. Individuals sense density-dependent GP in current habitat, abiotic GP in alternative habitat # Each mover uses its own dominance-weighted effective competitor density (fish per unit # area, self excluded) in its current patch. # Effective competitor density: grouped by patch × age_class when # age_structured_competition = TRUE, by patch only when FALSE. eff_density_all <- fish_pop |> mutate(age_class = if_else(cohort == current_year, "age0", "age1plus")) |> group_by(across(all_of(eff_grp))) |> mutate( A_patch = if_else(first(patch) == "warm", A_warm, A_cold), eff_density = (fncEffDensity(weight, beta = dominance_beta) - 1) / A_patch ) |> ungroup() |> pull(eff_density) eff_dens_current_vec <- eff_density_all[mover_rows] # DD-adjusted ration uses each mover's per-area effective competitor density ra_warm_dd_vec <- fish_pop$cmax_allometric[mover_rows] * pcmax_warm_adj * (K_warm / (K_warm + eff_dens_current_vec)) ra_cold_dd_vec <- fish_pop$cmax_allometric[mover_rows] * pcmax_cold_adj * (K_cold / (K_cold + eff_dens_current_vec)) ra_idx_warm_dd_vec <- pmax(1L, pmin(400L, map_int(ra_warm_dd_vec, ~which.min(abs(ra_seq_ibm - .x))))) ra_idx_cold_dd_vec <- pmax(1L, pmin(400L, map_int(ra_cold_dd_vec, ~which.min(abs(ra_seq_ibm - .x))))) # Warm fish sense dd-adjusted warm / abiotic cold; cold fish sense abiotic warm / dd-adjusted cold g_warm_vec <- wt_growth[cbind(wt_idx_warm, ifelse(in_warm, ra_idx_warm_dd_vec, ra_idx_warm_vec), ma_idx_vec)] g_cold_vec <- wt_growth[cbind(wt_idx_cold, ifelse(in_warm, ra_idx_cold_vec, ra_idx_cold_dd_vec), ma_idx_vec)] g_stay <- ifelse(in_warm, g_warm_vec, g_cold_vec) g_move_net <- ifelse(in_warm, g_cold_vec - move_cost_vec, g_warm_vec - move_cost_vec) } else if (sense_environment == "density_all") { ### 2.A.3. Sequential dominance-ordered patch selection. # Fish choose patches in descending body-size order within each competitive age class # (largest / most dominant first). Each fish senses density in both patches from the # fish already committed this step — not yesterday's pre-movement snapshot. This # eliminates the mass-synchronous overcrowding caused by the old pre-movement # approximation, and produces biologically realistic distributions: dominant fish # claim the thermally optimal patch, subordinates fill in around them. # # All fish choose anew each day (no incumbency advantage). Movement cost applies # only when the chosen patch differs from the fish's patch at the start of this day # (captured in prev_patch above). # # Patch assignments are written directly to fish_pop$patch here; the generic 2.B # stochasticity block below is skipped for this sensing mode. mover_weights_vec <- fish_pop$weight[mover_rows] age_class_all <- if_else(fish_pop$cohort == current_year, "age0", "age1plus") mover_age_class <- age_class_all[mover_rows] # Scalar helper: effective competitor density for one focal fish (weight wi) against a # vector of already-placed fish weights (comp_wts), per unit area. Returns Inf when # area == 0 (patch inaccessible → growth → 0); returns 0 when no competitors yet. eff_dens_scalar <- function(wi, comp_wts, area) { if (area == 0) return(Inf) if (length(comp_wts) == 0) return(0) sum(ifelse(comp_wts >= wi, 1, (comp_wts / wi)^dominance_beta)) / area } # Running placement accumulators: weights of fish already assigned this step. # Split by age class when age_structured_competition = TRUE (independent pools). if (age_structured_competition) { placed_warm_age0 <- numeric(0); placed_cold_age0 <- numeric(0) placed_warm_age1p <- numeric(0); placed_cold_age1p <- numeric(0) } else { placed_warm_all <- numeric(0); placed_cold_all <- numeric(0) } # Processing groups: independent competitive pools, each sorted dominant-first. # When age_structured_competition = TRUE, age-0 and age-1+ are separate pools. proc_groups <- if (age_structured_competition) { list( list(idx = which(mover_age_class == "age1plus"), ac = "age1plus"), list(idx = which(mover_age_class == "age0"), ac = "age0") ) } else { list(list(idx = seq_along(mover_rows), ac = "all")) } for (grp in proc_groups) { grp_idx <- grp$idx # positions within mover_rows / mover_weights_vec ac <- grp$ac if (length(grp_idx) == 0) next # Sort dominant-first within this age class grp_idx <- grp_idx[order(mover_weights_vec[grp_idx], decreasing = TRUE)] for (ii in grp_idx) { fi <- mover_rows[ii] # row index in fish_pop wi <- mover_weights_vec[ii] # fish weight (g) cmi <- fish_pop$cmax_allometric[fi] ma_idx_i <- pmax(1L, pmin(4500L, round(wi))) prev_i <- prev_patch[ii] # patch at start of this day # Fetch placement accumulators for this age class if (age_structured_competition) { pw <- if (ac == "age0") placed_warm_age0 else placed_warm_age1p pc <- if (ac == "age0") placed_cold_age0 else placed_cold_age1p } else { pw <- placed_warm_all pc <- placed_cold_all } # Effective density in each patch given already-placed fish this step ed_warm <- eff_dens_scalar(wi, pw, A_warm) ed_cold <- eff_dens_scalar(wi, pc, A_cold) # Density-adjusted rations in each patch ra_warm_i <- cmi * pcmax_warm_adj * (K_warm / (K_warm + ed_warm)) ra_cold_i <- cmi * pcmax_cold_adj * (K_cold / (K_cold + ed_cold)) ra_idx_w <- pmax(1L, pmin(400L, which.min(abs(ra_seq_ibm - ra_warm_i)))) ra_idx_c <- pmax(1L, pmin(400L, which.min(abs(ra_seq_ibm - ra_cold_i)))) # Growth lookup in each patch g_warm_i <- wt_growth[wt_idx_warm, ra_idx_w, ma_idx_i] g_cold_i <- wt_growth[wt_idx_cold, ra_idx_c, ma_idx_i] # Effective growth after movement cost (cost applies only when switching patches). # move_cost_vec[ii] is already computed for all movers above. move_cost_i <- move_cost_vec[ii] gwarm_eff_i <- g_warm_i - if (prev_i == "warm") 0 else move_cost_i gcold_eff_i <- g_cold_i - if (prev_i == "cold") 0 else move_cost_i # Patch choice — same stochasticity rule as other sensing modes chose_warm <- if (move_stochastic == "none") { gwarm_eff_i >= gcold_eff_i } else if (move_stochastic == "prob") { runif(1) < fncMoveSoftmax(gwarm = gwarm_eff_i, gcold = gcold_eff_i, tau = tau) } else { # "indiv": threshold is a per-fish bias on the warm–cold advantage (gwarm_eff_i - gcold_eff_i) > fish_pop$move_threshold[fi] } # Assign patch and update placement accumulator for subsequent fish fish_pop$patch[fi] <- if (chose_warm) "warm" else "cold" if (age_structured_competition) { if (ac == "age0") { if (chose_warm) placed_warm_age0 <- c(placed_warm_age0, wi) else placed_cold_age0 <- c(placed_cold_age0, wi) } else { if (chose_warm) placed_warm_age1p <- c(placed_warm_age1p, wi) else placed_cold_age1p <- c(placed_cold_age1p, wi) } } else { if (chose_warm) placed_warm_all <- c(placed_warm_all, wi) else placed_cold_all <- c(placed_cold_all, wi) } } # end fish loop } # end group loop } ## 2.B. Impose stochasticity in movement # density_all handles patch assignment and stochasticity within its sequential loop # above (choices written directly to fish_pop$patch). Skip for that sensing mode. if (sense_environment != "density_all") { if (move_stochastic == "none") { ### 2.B.1. Deterministic: All movers do the same thing # Update patch based on whether moving is beneficial. # 1. If fish is in warm patch and growth potential of moving exceeds growth potential of staying, move to cold. # 2. If fish is in cold patch and growth potential of moving exceeds growth potential of staying, move to warm. fish_pop$patch[mover_rows] <- case_when( in_warm & g_move_net > g_stay ~ "cold", !in_warm & g_move_net > g_stay ~ "warm", TRUE ~ fish_pop$patch[mover_rows] ) } else if (move_stochastic == "prob") { ### 2.B.2. Probabilistic (softmax) rule: Instead of a hard threshold, fish switch with a ### probability that increases with the growth advantage. One parameter: τ (sensitivity; ### small = nearly deterministic, large = nearly random). # Use cost-adjusted net growth rates: staying is free, moving incurs move_cost_vec. # g_stay = growth in current patch (no cost) # g_move_net = growth in alternate patch - move_cost_vec # Re-map to warm/cold perspective for fncMoveSoftmax: gwarm_eff <- ifelse(in_warm, g_stay, g_move_net) # effective growth if occupying warm patch gcold_eff <- ifelse(in_warm, g_move_net, g_stay) # effective growth if occupying cold patch p_warm <- fncMoveSoftmax(gwarm = gwarm_eff, gcold = gcold_eff, tau = tau) # calculate probability of selecting warm habitat choose_warm <- runif(length(mover_rows)) < p_warm # choose warm if p_warm exceeds a random number between 0 and 1 fish_pop$patch[mover_rows] <- if_else(choose_warm, "warm", "cold") # update patch } else if (move_stochastic == "indiv") { ### 2.B.3. Individual-level movement threshold: each fish has a unique threshold for ### movement, representing individual variation in boldness/dispersal propensity fish_pop$patch[mover_rows] <- case_when( in_warm & (g_move_net - g_stay) > fish_pop$move_threshold[mover_rows] ~ "cold", !in_warm & (g_move_net - g_stay) > fish_pop$move_threshold[mover_rows] ~ "warm", TRUE ~ prev_patch ) } } } # end patch choice # Recompute func_temp and pcmax_adjusted based on post-choice patch. # Required because some fish may have switched patches during step 2; # step 1 values reflect the pre-choice patch and would be stale for switchers. fish_pop <- fish_pop |> mutate( func_temp = ifelse(patch == "warm", fncTempDepend(T_warm), fncTempDepend(T_cold)), pcmax_baseline = ifelse(patch == "warm", pcmax_warm, pcmax_cold), pcmax_adjusted = ifelse(func_temp < pcmax_baseline, func_temp, pcmax_baseline) ) # 3. UPDATE HABITAT QUALITY (DENSITY-DEPENDENT PCMAX / RATION) ───────────── # Re-apply zero-area guard after patch choice: movers may have re-selected a # collapsed patch (A = 0) during step 2. Force them back to the surviving patch # before density calculations to prevent division-by-zero / NaN rations. if (A_cold == 0 && any(fish_pop$patch == "cold")) { fish_pop$patch[fish_pop$patch == "cold"] <- "warm" } if (A_warm == 0 && any(fish_pop$patch == "warm")) { fish_pop$patch[fish_pop$patch == "warm"] <- "cold" } # Tally switches and define `switched` AFTER zero-area guards: fish forced # back to their original patch by a guard are not counted as having switched # and incur no movement cost. This prevents spurious energy drain in scenarios # where one patch has A = 0 and the other patch looks more attractive (fish # "choose" the zero-area patch every step but get reversed by the guard). switched <- logical(length(mover_rows)) # default FALSE (no movers) if (length(mover_rows) > 0) { switched <- fish_pop$patch[mover_rows] != prev_patch switches[fish_pop$pid[mover_rows]] <- switches[fish_pop$pid[mover_rows]] + as.integer(switched) } # Update density (fish per unit area) and dominance-weighted effective competitor density # after fish have moved. density counts all fish in the patch regardless of age class. # eff_density is computed within the grouping defined by eff_grp: # age_structured_competition = TRUE → group by patch × age_class (fry vs. age-1+) # age_structured_competition = FALSE → group by patch only (all fish compete together) # This reflects ontogenetic habitat segregation in salmonids when enabled. Self is # excluded via the (-1) correction from fncEffDensity. fish_pop <- fish_pop |> mutate(age_class = if_else(cohort == current_year, "age0", "age1plus")) |> group_by(patch) |> mutate( A_patch = if_else(first(patch) == "warm", A_warm, A_cold), density = n() / A_patch ) |> group_by(across(all_of(eff_grp))) |> mutate( eff_density = (fncEffDensity(weight, beta = dominance_beta) - 1) / A_patch ) |> ungroup() |> select(-A_patch, -age_class) if (food_densdepen == "fullerton") { ### 3a. Fullerton approach (uses raw patch density; dominance not yet integrated here) fdens <- fish_pop$density fdens[fdens > MaxDensity4Growth] <- MaxDensity4Growth density_effect <- fncRescale((1 - c(fdens, 0.01, MaxDensity4Growth)), c(0.5, 1)) density_effect <- density_effect[-c(length(density_effect), length(density_effect) - 1)] fish_pop <- fish_pop |> mutate(ration = ration * density_effect) } else if (food_densdepen == "hyperbolic") { ### 3b. Hyperbolic reduction — each fish's P_Cmax is penalised by its own effective density # (eff_density, in fish per unit area, self excluded). Large dominant fish face a lower # eff_density → less suppression; small subordinate fish face a higher eff_density → more. # No -1 correction needed: self is already excluded from eff_density. fish_pop <- fish_pop |> mutate( k = if_else(patch == "warm", K_warm, K_cold), pcmax_adjusted_dd = pcmax_adjusted * (k / (k + eff_density)), ration = cmax_allometric * pcmax_adjusted_dd ) |> select(-k) } # 4. GROW FISH ───────────────────────────────────────────────────────────── # Growth lookup for all fish via fncGrowthFish. Growth in g/g/d. growth_df <- fncGrowthFish(NA, as.data.frame(fish_pop), T_warm, T_cold) # Apply movement cost: deduct from growth rate for fish that actually switched patches. # move_cost_vec is in g/g/d and is indexed over mover_rows; switched is the logical # subset of those rows that changed patch this step. if (length(mover_rows) > 0 && any(switched)) { switched_rows <- mover_rows[switched] growth_df$growth[switched_rows] <- growth_df$growth[switched_rows] - move_cost_vec[switched] } # Update weight: W_{t+1} = W_t + (growth_rate * W_t) fish_pop$weight <- fish_pop$weight + (growth_df$growth * fish_pop$weight) # Update peak weight: ratchet up only, never down fish_pop$peak_weight <- pmax(fish_pop$peak_weight, fish_pop$weight) # 5. SPAWNING AND REPRODUCTION ───────────────────────────────────────────── # Reset spawned_this_year flag at the calendar year boundary if (d > 1) { prev_year <- year(habitat_df$date[d - 1]) if (current_year != prev_year) fish_pop$spawned_this_year <- FALSE } # Day-of-year and relative condition needed for spawning probability doy_d <- yday(habitat_df$date[d]) condition_spw <- fish_pop$weight / fish_pop$peak_weight # Combined daily spawning probability: size × condition × date p_spawn <- fncMaturitySize(fish_pop$weight) * fncMaturityCondition(condition_spw) * fncMaturityDate(doy_d) # Fish that have already spawned this year are ineligible p_spawn[fish_pop$spawned_this_year] <- 0 # Bernoulli trial: which fish spawn today? spawns <- as.logical(rbinom(nrow(fish_pop), size = 1, prob = p_spawn)) # For spawning fish: calculate fecundity, impose energetic cost, mark as spawned if (any(spawns)) { spawner_idx <- which(spawns) # Fecundity: number of eggs as a function of weight, with log-scale noise n_eggs <- round(fncFecundBromage(fish_pop$weight[spawner_idx], sigma = sigma_fecund, survival = egg_surv)) # Capture pre-cost weight for logging (spawn decision was made at this weight) weight_at_spawn <- fish_pop$weight[spawner_idx] # Energetic cost of reproduction: deduct from parent weight # Weight floor at start_wt — spawning cannot reduce a fish below hatch weight reproduction_cost <- weight_at_spawn * repro_cost fish_pop$weight[spawner_idx] <- pmax(start_wt, weight_at_spawn - reproduction_cost) # Flag these fish as having spawned; they are ineligible to spawn again this year fish_pop$spawned_this_year[spawner_idx] <- TRUE # Create offspring rows and add to population -------------------------------- new_fish_list <- vector("list", length(spawner_idx)) for (i in seq_along(spawner_idx)) { si <- spawner_idx[i] parent <- fish_pop[si, ] parent$weight <- weight_at_spawn[i] # restore pre-cost weight for logging n_off <- n_eggs[i] new_pids <- seq(next_pid, next_pid + n_off - 1L) # Initial offspring weights drawn from same distribution as founders off_wt <- pmax(min_wt, rnorm(n_off, start_wt, wt_sd)) # Offspring inherit strategy and patch from parent; movement threshold drawn fresh off_thresh <- if (parent$strategy == "optimal_mover") { rnorm(n_off, mean = 0, sd = sigma_bold) } else { rep(Inf, n_off) } new_fish_list[[i]] <- tibble( strategy = parent$strategy, patch = parent$patch, weight = off_wt, pid = new_pids, cmax_allometric = NA_real_, pcmax_baseline = NA_real_, pcmax_adjusted = NA_real_, pcmax_adjusted_dd = NA_real_, func_temp = NA_real_, ration = NA_real_, move_threshold = off_thresh, peak_weight = off_wt, age_days = 0L, spawned_this_year = FALSE, parent_pid = parent$pid, cohort = current_year, prob_surv = 1, survive = 1 ) # Log the spawning event (weight and condition are pre-spawn values) spawn_log <- bind_rows(spawn_log, tibble( parent_pid = parent$pid, dayofsim = d, n_offspring = n_off, weight = parent$weight, condition = parent$weight / parent$peak_weight, offspring_pid_start = next_pid, offspring_pid_end = next_pid + n_off - 1L )) next_pid <- next_pid + n_off } # Bind offspring and add them to the population BEFORE survival. # Size-based and density-driven starvation mortality will thin the cohort naturally. new_fish <- bind_rows(new_fish_list) n_new <- nrow(new_fish) switches <- c(switches, integer(n_new)) # extend switch counter for new fish # Extend growth_df with stub rows for offspring so survival indexing stays aligned. # WT.actual is set to their patch temperature; growth is NA (they didn't feed today). off_wt_actual <- ifelse(new_fish$patch == "warm", T_warm, T_cold) growth_df <- bind_rows(growth_df, data.frame(WT.actual = off_wt_actual, growth = NA_real_)) # Register and add offspring to live population fish_registry <- bind_rows( fish_registry, new_fish |> select(pid, strategy, parent_pid, cohort) |> mutate(birth_dayofsim = d) ) fish_pop <- bind_rows(fish_pop, new_fish) } # 6. SURVIVAL ───────────────────────────────────────────────────────────── s_max_vec <- ifelse(fish_pop$patch == "warm", S_max_warm, S_max_cold) # Size-dependent survival: logistic sigmoid in weight, rising from minprob floor to habitat-specific maxprob p_sg <- fncSurviveSize(fish_pop$weight, minprob = s_min, maxprob = s_max_vec, w0 = s_w0, k = s_k)[[1]] # Temperature-dependent survival: multiplied with size/growth probability so both act simultaneously p_temp <- fncSurviveTemp(growth_df$WT.actual, T1 = T1_mort, T9 = T9_mort) # Condition-based (starvation) survival: relative condition = current / peak weight condition <- fish_pop$weight / fish_pop$peak_weight p_starv <- fncSurviveStarve(condition, K9 = K9_starv, K1 = K1_starv) # Age-based survival (senescence): survival probability declines past age_thresh years old. p_age <- fncSurviveAge(fish_pop$age_days / 365, x0 = age_x0, k = age_k, p_min = age_pmin) # Critical-period consumption-based survival (Elliott 1989) — disabled. # Replaced by the sigmoid shape of fncSurviveSize, which sustains strong size-dependent # mortality below w0 (~7g) and achieves the same growth-mediated density dependence # without requiring a separate consumption threshold module. # p_crit <- fncSurviveConsumption( # pcmax_dd = fish_pop$pcmax_adjusted_dd, # age_days = fish_pop$age_days, # crit_pcmax_lo = crit_pcmax_lo, # crit_pcmax_hi = crit_pcmax_hi, # crit_period_days = crit_period_days # ) # Calculate combined daily survival rate prb.srv <- pmin(p_sg * p_temp * p_starv * p_age, 1) # Minimum weight floor: fish that drop below hatch weight die immediately (backstop) prb.srv[fish_pop$weight < start_wt] <- 0 survivors <- rbinom(nrow(growth_df), size = 1, prob = prb.srv) growth_df <- growth_df |> mutate(prob_surv = prb.srv, survive = survivors) # 7. STORE RESULTS AND REMOVE NON-SURVIVORS ──────────────────────────────── ibm_records[[d]] <- data.frame( pid = fish_pop$pid, dayofsim = d, strategy = fish_pop$strategy, weight = fish_pop$weight, patch = fish_pop$patch, temp = growth_df$WT.actual, # includes offspring stub (their patch temp) ggd = growth_df$growth, # NA for offspring on birth day survived = growth_df$survive, condition = condition, age = fish_pop$age_days / 365, p_survive = prb.srv # worth keeping — used in later analyses ) # Remove non-survivors — only living fish carry forward to next iteration fish_pop <- fish_pop[growth_df$survive == 1, ] # exit loop if all fish have died if (nrow(fish_pop) == 0) break } # ── Collate ibm_long ───────────────────────────────────────────────────────── ibm_long <- bind_rows(ibm_records) |> left_join(fish_registry |> select(pid, parent_pid, cohort, birth_dayofsim), by = "pid") |> left_join(habitat_df |> select(dayofsim, date), by = "dayofsim") # ── Standard daily summary by strategy ─────────────────────────────────────── ibm_summary <- ibm_long |> group_by(strategy, dayofsim, date) |> summarise( mean_weight = mean(weight, na.rm = TRUE), sd_weight = sd(weight, na.rm = TRUE), mean_temp = mean(temp, na.rm = TRUE), sd_temp = sd(temp, na.rm = TRUE), mean_ggd = mean(ggd, na.rm = TRUE), sd_ggd = sd(ggd, na.rm = TRUE), prop_warm = mean(patch == "warm", na.rm = TRUE), n_alive = sum(survived, na.rm = TRUE), mean_condition = mean(condition, na.rm = TRUE), sd_condition = sd(condition, na.rm = TRUE), .groups = "drop" ) list( ibm_records = ibm_records, ibm_long = ibm_long, ibm_summary = ibm_summary, fish_registry = fish_registry, spawn_log = spawn_log, switches = switches, params = params, habitat_df = habitat_df ) } ``` ## Write R file For sourcing by other scripts. ```{r} #| eval: false #| purl: false file.remove("SimulationLoop.R") quarto::qmd_to_r_script("SimulationLoop.qmd") ```

4.1 Design notes

4.1.1 Key modules

4.1.2 Configuration options

4.1.3 Senstive parameters

4.1.4 Misc. notes

4.1.5 To do

4.2 run_simulation()

4.3 Write R file

4.2 `run_simulation()`