Refactor nested for loop into smaller testable functions

birds_in_trees.R

### Refactored Birds in Trees code
library(dplyr)
library(mefa4)

# Function to create list of list of args for outer calculation function
#   returns List of list(stand_id, species_list, stand_data)
make_stand_args <- function(stand_id, data){
  splist <- data %>%
    filter(StandID == stand_id) %>%
    dplyr::select(Species) %>%
    distinct()%>%
    pull(Species)

  result <- list(stand_id     = stand_id,
                 species_list = splist)
  return(result)
}

# Function to do the inner calculation
#   calc the mden of single species,
#    given duration and distance intervals of all species, and tau lookup table
calc_mden <- function(species_id, dur_intervals, dis_intervals, tau_ref){

  # duration variable phi for species
  y_dur <- as.matrix(dur_intervals[[species_id]])
  d_dur <- matrix(c(3, 5, 10), nrow(y_dur), 3, byrow=TRUE,
                  dimnames=dimnames(y_dur))

  # distance variable tau
  y_dis <- as.matrix(dis_intervals[[species_id]])
  d_dis <- matrix(c(0.5, 1), nrow(y_dis), 2, byrow=TRUE,
                  dimnames=dimnames(y_dis))

  # grab tau for specific species
  tau <- tau_ref[which(tau_ref$Sp == species_id),3]
  rmax <- apply(d_dis, 1, max)
  q <- (tau^2 / rmax^2) * (1-exp(-(rmax/tau)^2))
  A <- pi * rmax^2
  a_q <- A * q
  y_total <- rowSums(y_dur)
  m_den <- mean(y_total/a_q$tau_best)
  return(m_den)
}

# Outer calculation function
#   return a list of mdens for each stand
#    Each list of mdens would have one entry for each species.
calc_stand_mdens <- function(stand_id, species_list, tau_ref, all_dis, all_dur){

  names(species_list) <- species_list
  species_mdens <- lapply(X=species_list, FUN=calc_mden,
                          dur_intervals=all_dur, dis_intervals=all_dis,
                          tau_ref=tau_ref)

  return(species_mdens)
}

# Wrap do.call because I hate trying to stuff it into lapply directly.
#   Combine dynamic and static args into single arguments list
splat_wrapper <- function(dynamic_args, static_args){
  args <- c(dynamic_args, static_args)
  do.call(what=calc_stand_mdens, args=args)
}

# Input Parameters
raw_data <- source('~/Downloads/SampleBirdData.txt')$value
tau_ref <- source('~/Downloads/tau_ref.txt')$value
stand_list <- source('~/Downloads/sample_standlist.txt')$value

# Correct column is a list of lists issue in source data
qpadbirdData <- raw_data
qpadbirdData$Species <- unlist(raw_data$Species, use.names=F)

# Calculate the crosstabs for distance and duration for all species
# time Intervals
all_dur <- mefa4::Xtab(~ VisitID + Dur + Species, qpadbirdData)
# distance intervals
all_dis <- mefa4::Xtab(~ VisitID + Dis + Species, qpadbirdData)

# Create static args list of the tau_refs and crosstabs
static_args <- list(tau_ref=tau_ref,
                    all_dur=all_dur,
                    all_dis=all_dis)

# Given a list of stands and the source data,
#  make a list of arguments lists for the calculation function
dynamic_args <- lapply(X = stand_list,
                       FUN = make_stand_args,
                       data = qpadbirdData)

# Now run the foo function once for each element (set of args) from stand_args
names(stand_args) <- stand_list
stand_calcs <- lapply(X = dynamic_args,
                      FUN = splat_wrapper,
                      static_args = static_args)

sample_standlist.txt

c("BEF 15", "NUL 10", "BEF 16", "NUL 2", "BEF 17", "NUL 8", "BEF 14", 
"NUL 14", "NUL 1", "NUL 13", "NUL 15", "NUL 6", "BEF 11", "BEF 3", 
"NUL 3", "BEF 2", "BEF 8", "NUL 4", "BEF 9", "BEF 7", "NUL 9", 
"NUL 7", "NUL 5", "BEF 1", "NUL 12", "NUL 11", "BEF 6", "BEF 12"
)

SampleBirdData.txt

structure(list(StandID = c("BEF 15", "NUL 10", "BEF 16", "NUL 2", 
"BEF 17", "NUL 8", "BEF 17", "BEF 17", "BEF 14", "NUL 14", "NUL 1", 
"NUL 13", "NUL 15", "NUL 6", "NUL 13", "BEF 11", "BEF 14", "NUL 15", 
"BEF 16", "BEF 3", "NUL 14", "BEF 3", "BEF 14", "NUL 3", "BEF 11", 
"BEF 2", "NUL 2", "BEF 16", "BEF 8", "BEF 3", "BEF 8", "NUL 4", 
"BEF 9", "BEF 17", "BEF 7", "NUL 14", "BEF 17", "NUL 10", "NUL 4", 
"NUL 13", "NUL 9", "NUL 10", "BEF 3", "BEF 2", "BEF 3", "NUL 13", 
"NUL 15", "NUL 6", "NUL 4", "NUL 7", "BEF 11", "NUL 10", "NUL 5", 
"NUL 9", "BEF 1", "NUL 12", "NUL 5", "BEF 8", "NUL 7", "NUL 8", 
"BEF 9", "BEF 1", "NUL 2", "NUL 8", "BEF 7", "BEF 8", "BEF 9", 
"NUL 8", "BEF 8", "NUL 15", "NUL 5", "NUL 14", "BEF 9", "NUL 13", 
"NUL 10", "NUL 13", "BEF 14", "NUL 8", "NUL 6", "NUL 10", "BEF 17", 
"BEF 8", "NUL 11", "BEF 6", "BEF 16", "BEF 2", "NUL 7", "NUL 14", 
"BEF 16", "NUL 5", "NUL 1", "BEF 15", "NUL 9", "BEF 12", "BEF 11", 
"BEF 14", "NUL 12", "NUL 5", "NUL 14", "NUL 3"), Species = structure(list(
    Species = c("NOPA", "HETH", "REVI", "RBNU", "BTNW", "MYWA", 
    "OVEN", "BTNW", "WIWR", "BLBW", "BTNW", "HETH", "SCTA", "YBSA", 
    "RBNU", "BLBW", "GCKI", "BTBW", "REVI", "REVI", "REVI", "REVI", 
    "RBNU", "REVI", "OVEN", "OVEN", "HETH", "REVI", "BLBW", "RCKI", 
    "BTBW", "OVEN", "VEER", "BTNW", "BCCH", "REVI", "WIWR", "BCCH", 
    "OVEN", "BLJA", "YBSA", "MYWA", "BTNW", "BCCH", "PAWA", "BTBW", 
    "NOFL", "REVI", "YBSA", "REVI", "CSWA", "SOSP", "BTNW", "REVI", 
    "BTBW", "BOCH", "YBSA", "CAWA", "BTNW", "HETH", "DEJU", "HETH", 
    "WTSP", "YBSA", "BTBW", "BLJA", "VEER", "REVI", "BTBW", "BTNW", 
    "REVI", "BTBW", "BTBW", "YBSA", "AMRE", "RUGR", "BBCU", "RBNU", 
    "BLBW", "REVI", "BTBW", "NAWA", "BTBW", "GCKI", "BCCH", "RBNU", 
    "HETH", "OVEN", "BTNW", "BLBW", "REVI", "REVI", "REVI", "ALFL", 
    "COYE", "HETH", "REVI", "SWTH", "OVEN", "REVI")), row.names = c(NA, 
-100L), class = c("tbl_df", "tbl", "data.frame")), Dur = c("5-10min", 
"0-3min", "0-3min", "5-10min", "5-10min", "5-10min", "0-3min", 
"5-10min", "0-3min", "3-5min", "5-10min", "0-3min", "0-3min", 
"3-5min", "0-3min", "0-3min", "5-10min", "0-3min", "0-3min", 
"5-10min", "3-5min", "5-10min", "0-3min", "0-3min", "0-3min", 
"0-3min", "3-5min", "5-10min", "0-3min", "0-3min", "0-3min", 
"0-3min", "0-3min", "0-3min", "0-3min", "0-3min", "0-3min", "3-5min", 
"0-3min", "0-3min", "3-5min", "0-3min", "0-3min", "0-3min", "3-5min", 
"0-3min", "5-10min", "0-3min", "3-5min", "0-3min", "0-3min", 
"0-3min", "0-3min", "3-5min", "0-3min", "0-3min", "5-10min", 
"0-3min", "0-3min", "0-3min", "0-3min", "0-3min", "0-3min", "5-10min", 
"0-3min", "5-10min", "3-5min", "5-10min", "3-5min", "3-5min", 
"5-10min", "0-3min", "0-3min", "0-3min", "0-3min", "0-3min", 
"3-5min", "5-10min", "3-5min", "0-3min", "0-3min", "0-3min", 
"0-3min", "0-3min", "0-3min", "0-3min", "0-3min", "0-3min", "5-10min", 
"0-3min", "0-3min", "0-3min", "0-3min", "0-3min", "0-3min", "5-10min", 
"0-3min", "0-3min", "0-3min", "0-3min"), Dis = c("0-50m", "50-100m", 
"50-100m", "0-50m", "50-100m", "50-100m", "0-50m", "0-50m", "50-100m", 
"50-100m", "0-50m", "50-100m", "0-50m", "50-100m", "50-100m", 
"0-50m", "50-100m", "50-100m", "50-100m", "50-100m", "50-100m", 
"50-100m", "0-50m", "0-50m", "0-50m", "50-100m", "50-100m", "0-50m", 
"50-100m", "0-50m", "50-100m", "50-100m", "0-50m", "50-100m", 
"50-100m", "50-100m", "50-100m", "50-100m", "50-100m", "50-100m", 
"50-100m", "50-100m", "50-100m", "0-50m", "50-100m", "0-50m", 
"0-50m", "50-100m", "50-100m", "0-50m", "50-100m", "50-100m", 
"50-100m", "50-100m", "50-100m", "0-50m", "0-50m", "0-50m", "0-50m", 
"0-50m", "0-50m", "50-100m", "50-100m", "0-50m", "0-50m", "50-100m", 
"0-50m", "50-100m", "50-100m", "50-100m", "50-100m", "50-100m", 
"0-50m", "0-50m", "0-50m", "50-100m", "0-50m", "0-50m", "50-100m", 
"0-50m", "0-50m", "50-100m", "0-50m", "0-50m", "0-50m", "0-50m", 
"0-50m", "50-100m", "0-50m", "50-100m", "0-50m", "50-100m", "50-100m", 
"50-100m", "50-100m", "0-50m", "0-50m", "50-100m", "0-50m", "50-100m"
), VisitID = c("BEF 15 A 182", "NUL 10 A 154", "BEF 16 B 160", 
"NUL 2 D 172", "BEF 17 D 158", "NUL 8 B 154", "BEF 17 B 182", 
"BEF 17 C 202", "BEF 14 B 158", "NUL 14 B 193", "NUL 1 B 197", 
"NUL 13 D 152", "NUL 15 B 150", "NUL 6 D 173", "NUL 13 A 152", 
"BEF 11 A 159", "BEF 14 A 182", "NUL 15 B 150", "BEF 16 C 160", 
"BEF 3 A 179", "NUL 14 F 152", "BEF 3 D 179", "BEF 14 B 182", 
"NUL 3 D 197", "BEF 11 E 159", "BEF 2 E 162", "NUL 2 C 172", 
"BEF 16 C 180", "BEF 8 A 180", "BEF 3 B 179", "BEF 8 A 180", 
"NUL 4 A 196", "BEF 9 C 160", "BEF 17 F 158", "BEF 7 E 180", 
"NUL 14 D 176", "BEF 17 F 182", "NUL 10 D 194", "NUL 4 A 196", 
"NUL 13 E 193", "NUL 9 C 195", "NUL 10 F 194", "BEF 3 F 179", 
"BEF 2 E 162", "BEF 3 C 162", "NUL 13 C 193", "NUL 15 D 150", 
"NUL 6 C 173", "NUL 4 C 196", "NUL 7 A 154", "BEF 11 C 159", 
"NUL 10 D 154", "NUL 5 B 196", "NUL 9 A 154", "BEF 1 D 202", 
"NUL 12 B 153", "NUL 5 D 196", "BEF 8 F 180", "NUL 7 B 154", 
"NUL 8 C 154", "BEF 9 C 204", "BEF 1 E 179", "NUL 2 A 155", "NUL 8 D 174", 
"BEF 7 D 180", "BEF 8 E 180", "BEF 9 A 180", "NUL 8 A 195", "BEF 8 C 160", 
"NUL 15 E 150", "NUL 5 B 155", "NUL 14 B 176", "BEF 9 A 180", 
"NUL 13 D 193", "NUL 10 B 175", "NUL 13 F 176", "BEF 14 A 158", 
"NUL 8 C 154", "NUL 6 C 155", "NUL 10 E 194", "BEF 17 D 202", 
"BEF 8 F 180", "NUL 11 F 194", "BEF 6 A 178", "BEF 16 C 160", 
"BEF 2 B 162", "NUL 7 B 174", "NUL 14 D 193", "BEF 16 D 180", 
"NUL 5 B 196", "NUL 1 D 172", "BEF 15 C 158", "NUL 9 D 154", 
"BEF 12 F 159", "BEF 11 D 159", "BEF 14 B 201", "NUL 12 E 153", 
"NUL 5 A 155", "NUL 14 C 193", "NUL 3 E 172")), row.names = c(NA, 
-100L), class = c("tbl_df", "tbl", "data.frame"))

tau_ref.txt

structure(list(Sp = c("RTHU", "BLBW", "GCKI", "BOCH", "BTBW", 
"BRCR", "MAWA", "CMWA", "BAWW", "AMRE", "BHVI", "NOPA", "CSWA", 
"CAWA", "BTNW", "MOWA", "WIWR", "REVI", "YBFL", "RBGR", "LEFL", 
"PAWA", "BLPW", "DOWO", "YBSA", "OVEN", "DEJU", "CEDW", "SCTA", 
"VEER", "SOSP", "MYWA", "EAWP", "BCCH", "WBNU", "SWSP", "INBU", 
"BEKI", "NAWA", "HAWO", "ALFL", "GCFL", "AMGO", "WTSP", "SWTH", 
"LISP", "PUFI", "RCKI", "AMRO", "BLJA", "RBNU", "COYE", "RUGR", 
"NOWA", "HETH", "BBCU", "NOFL", "PIWO"), tau = c(0.298999, 0.470045, 
0.388122, 0.474593, 0.452795, 0.464176, 0.607747, 0.5223915, 
0.615874, 0.6017003, 0.560901, 0.563834, 0.747081, 0.625857, 
0.595276, 0.777684, 0.590151, 0.69453, 0.633242, 0.931138, 0.769924, 
0.631699821, 0.631699821, 0.6087145, 0.802057, 0.798386, 0.618904, 
0.7378935, 0.676319, 0.918361, 0.603962, 0.681626, 0.837621, 
0.8056039, 0.743488, 0.932883, 0.702112, 0.818148, 0.689438, 
0.782298, 0.978404, 0.914948, 0.836709, 1.169883, 0.7810633, 
0.777799, 0.794372, 0.898033, 0.929143, 1.074654, 0.932979, 0.770593, 
1.10043, 0.723556, 1.27119, 1.88797, 1.37911, 1.45686), tau_best = c(0.2110982, 
0.364099, 0.379388, 0.409157, 0.447038, 0.449505, 0.452729, 0.460238, 
0.461154, 0.4969422, 0.513969, 0.535236, 0.535671, 0.536698, 
0.537636, 0.537947, 0.545593, 0.550853, 0.58209, 0.58465, 0.589032, 
0.595738, 0.595738, 0.601205, 0.608763, 0.612224, 0.621086, 0.62365, 
0.631127, 0.666315, 0.66755, 0.677063, 0.706598, 0.7102644, 0.713479, 
0.719708, 0.747621, 0.751698, 0.753327, 0.7644505, 0.7648203, 
0.767414, 0.780721, 0.792044, 0.7953055, 0.811148, 0.83448, 0.879949, 
0.918185, 0.918586, 0.93674, 0.942806, 0.955041, 1.09869, 1.22448, 
1.25339, 1.41392, 1.43079)), class = c("tbl_df", "tbl", "data.frame"
), row.names = c(NA, -58L))

@grosscol it doesn't look like it is seing the tau_ref argument. I made a sample dataset that you can try to run. is there a way to upload it to this github?

I think I fixed that issue by adding tau_ref to the make_stand_args() function

make_stand_args <- function(stand_id, tau_ref, data){
  splist <- data %>%
    filter(StandID == stand_id) %>%
    dplyr::select(Species) %>%
    distinct()
  
  spdata <- data %>%
    filter(StandID == stand_id)
  
  result <- list(stand_id     = stand_id,
                 species_list = splist,
                 stand_data   = spdata,
                 tau_ref      = tau_ref)
  return(result)
}

Good catch. The invariant argument would need to be added to splat_wrapper, and the static args added to the dynamic ones. Derp.

Adding it to the args list works. I usually keep those args lists for things that change with each call invocation, but in this case I think your way is simpler.

As for uploading a sample data set, you could upload the output of dput for a sample of the data as your own gist. The r-help-readme channel probably has a bunch of tips about that.

ok I just uploaded a sample. I am getting an error that says this when you show traceback. it doesn't like the fact that all_dur is a list but if I get stand_calc_mdens to return all_dur and then run it seperately and use as.matrix(all_dur[[species_id]] manually it works fine

Error in h(simpleError(msg, call)) : 
error in evaluating the argument 'x' in selecting a method for function 'as.matrix': invalid subscript type 'list'
9.
h(simpleError(msg, call))
8.
.handleSimpleError(function (cond) 
.Internal(C_tryCatchHelper(addr, 1L, cond)), "invalid subscript type 'list'", 
base::quote(dur_intervals[[species_id]]))
7.
as.matrix(dur_intervals[[species_id]])
6.
FUN(X[[i]], ...)
5.
lapply(X = species_list, FUN = calc_mden, dur_intervals = all_dur, 
dis_intervals = all_dis, tau_ref = tau_ref)
4.
(function (stand_id, species_list, stand_data, tau_ref) 
{
all_dur <- Xtab(~VisitID + Dur + Species$Species, stand_data)
all_dis <- Xtab(~VisitID + Dis + Species$Species, stand_data) ...
3.
do.call(what = calc_stand_mdens, args = args)
2.
FUN(X[[i]], ...)
1.
lapply(X = stand_args, FUN = splat_wrapper)