WIP: A big mess after the last maize scenario runs where important information...

WIP: A big mess after the last maize scenario runs where important information was not stored properly on grendel.

Grendel/scenarios/numbers_report.Rmd

@@ -28,6 +28,8 @@ scaler_cutoff <- 7
 
 Document compiled: `r Sys.time()`
 
+# Numbers
+
 ```{r read-results, echo=FALSE}
 current <- path_desktop("goose", date_stamp)
 dir_ls(current, glob = "*.out") %>% 
@@ -133,3 +135,211 @@ numbers_per_day %>%
 ```
 
 
+# Maize
+
+```{r}
+date_stamp <- "2019-01-09"
+```
+
+
+```{r}
+current <- path_desktop("goose", date_stamp)
+dir_ls(current, glob = "*.out") %>% 
+  path_file() %>% 
+  path_ext_remove() %>% 
+  as.character() %>% 
+  str_extract("[[:digit:]]+") -> jobs
+
+# Grain distributions
+grain_dist_file <- glue("{jobs}_maize-availability-grain-distributions_{date_stamp}.fst")
+path(current, jobs, grain_dist_file) %>%
+  fst::read_fst() %>% 
+  as_tibble() -> grain_dists
+
+goose_numbers_file <- glue("{jobs}_maize-availability-goose-numbers_{date_stamp}.fst")
+path(current, jobs, goose_numbers_file) %>%
+  fst::read_fst() %>% 
+  as_tibble() %>% 
+  mutate(day = day - 365L*season,
+         pinkfoot = pf_families + pf_non_breeders,
+         greylag = gl_families + gl_non_breeders,
+         barnacle = bn_families + bn_non_breeders) %>% 
+  filter(param == "GOOSE_BN_STARTNO_SCALER",
+         day == 364) %>% 
+  filter(dplyr::row_number() <= 100) %>% 
+  group_by(value) %>% 
+  summarize(barnacle_tot = mean(barnacle)) -> bn_std
+
+path(current, jobs, goose_numbers_file) %>%
+  fst::read_fst() %>% 
+  as_tibble() %>% 
+  mutate(day = day - 365L*season,
+         pinkfoot = pf_families + pf_non_breeders,
+         greylag = gl_families + gl_non_breeders,
+         barnacle = bn_families + bn_non_breeders) %>% 
+  filter(param == "GOOSE_BN_STARTNO_SCALER",
+         day == 364) %>% 
+  filter(dplyr::between(dplyr::row_number(), 301, 400)) %>% 
+  group_by(value) %>% 
+  summarize(barnacle_tot = mean(barnacle)) -> bn_193
+
+path(current, jobs, goose_numbers_file) %>%
+  fst::read_fst() %>% 
+  as_tibble() %>% 
+  mutate(day = day - 365L*season,
+         pinkfoot = pf_families + pf_non_breeders,
+         greylag = gl_families + gl_non_breeders,
+         barnacle = bn_families + bn_non_breeders) %>% 
+  filter(param == "GOOSE_BN_STARTNO_SCALER",
+         day == 364) %>% 
+  filter(dplyr::between(dplyr::row_number(), 201, 300)) %>% 
+  group_by(value) %>% 
+  summarize(barnacle_tot = mean(barnacle)) -> bn_192
+
+left_join(bn_std, bn_193, by = "value", suffix = c("_std", "_193")) %>% 
+  left_join(bn_192, by = "value", suffix = c("", "_192"))
+  
+  left_join(grain_dists, by = c("season", "param", "value")) %>% 
+  mutate(pinkfoot = pf_families + pf_non_breeders,
+         greylag = gl_families + gl_non_breeders,
+         barnacle = bn_families + bn_non_breeders,
+         day = day - 365L*season) %>% 
+  select(-ends_with("families"), -ends_with("non_breeders")) -> goose_numbers 
+
+
+
+# Grain distributions
+"/Users/au206907/Desktop/goose/2019-01-09/2090033/2090033_WD05_2019-01-08_scenario-results.txt" %>% 
+  read_tsv(col_names = c("season", "day", "pinkfoot", "pinkfoot_dist",
+           "barnacle", "barnacle_dist", "greylag", "greylag_dist",
+           "grain", "maize", "polyref", "param", "value")) -> wd5
+
+wd5 %>% 
+  left_join(wd5_grain_dist, by = c("value", "season")) %>% 
+  filter(grain_dist == "low",
+         barnacle != 0) %>% 
+  mutate(day = day - 365L*season) %>% 
+  select(day, season, value, maize, barnacle) %>% 
+  filter(day < 364 & day > 250) %>% 
+  group_by( value) %>% 
+  summarise(total_maize = sum(maize),
+            total_bn = sum(barnacle)) %>% 
+  arrange(value)
+wd7 %>% 
+  left_join(wd7_grain_dist, by = c("value", "season")) %>% 
+  filter(grain_dist == "low",
+         barnacle != 0) %>% 
+  mutate(day = day - 365L*season) %>% 
+  select(day, season, value, maize, barnacle) %>% 
+  filter(day < 364 & day > 250) %>% 
+  group_by( value) %>% 
+  summarise(total_maize = sum(maize),
+            total_bn = sum(barnacle)) %>% 
+  arrange(value)
+
+wd8 %>% 
+  select(day, season, value, maize) %>% 
+  filter(day == 365+300) %>% 
+  group_by(season, value) %>% 
+  summarise(total_maize = sum(maize))
+
+
+wd5 %>% 
+  select(season, grain, value) %>% 
+  group_by(value, season) %>% 
+  summarise(max_grain = max(grain)) %>% 
+  mutate(grain_dist = case_when(max_grain > 700 ~ "high",
+                                TRUE ~ "low")) %>% 
+  ungroup() %>% 
+  select(value, season, grain_dist) -> wd5_grain_dist
+
+n_seasons <- 10
+
+# Mean numbers per day
+wd8 %>% 
+  left_join(wd8_grain_dist, by = c("value", "season")) %>% 
+  mutate(day = day - 365*season) %>% 
+  select(grain_dist, param, value, season, day, pinkfoot, barnacle, greylag) %>% 
+  group_by(grain_dist, param, value, season, day) %>% 
+  summarise_at(.vars = vars(pinkfoot, barnacle, greylag),
+               .funs = funs(daily = sum)) %>% 
+  group_by(grain_dist, param, value, day) %>% 
+  summarise_at(.vars = vars(ends_with("daily")),
+               .funs = funs(avg = mean,
+                            min, max)) %>% 
+  ungroup() %>% 
+  gather(key = species, value = numbers, -grain_dist, -param, -value, -day) %>% 
+  filter(numbers > 0) -> numbers_per_day
+
+numbers_per_day %>% 
+  filter(value == 1) -> default
+left_join(numbers_per_day, default, 
+          by = c("grain_dist", "param",  "day", "species"),
+          suffix = c("", "_default")) %>% 
+  mutate(standarized = numbers/numbers_default) -> numbers_per_day
+
+numbers_per_day %>% 
+  filter(day == 364,
+         species == "barnacle_daily_avg") %>% 
+  ggplot(aes(value, standarized)) +
+  geom_point() + 
+  facet_grid(~grain_dist)
+  
+
+goose_sp <- "pinkfoot_daily_avg"
+numbers_per_day %>% 
+  filter(param == "GOOSE_BN_STARTNO_SCALER",
+         species == goose_sp) %>% 
+  ggplot(aes(day, standarized)) +
+  geom_line(aes(color = factor(value))) + 
+  scale_color_viridis_d(name = "Barnacle scaler") +
+  hrbrthemes::theme_ipsum_rc(axis_title_size = 12) +
+  labs(title = "Daily average numbers",
+       subtitle = "Pinkfeet numbers as a function of barnacle numbers",
+       caption = "Numbers are average of 10 consecutive seasons",
+       y = "proportion of default numbers") +
+  facet_grid(~grain_dist)
+
+
+
+
+
+goose_numbers_file <- glue("{jobs}_maize-availability-goose-numbers_{date_stamp}.fst")
+path(current, jobs, goose_numbers_file) %>%
+  fst::read_fst() %>% 
+  as_tibble() %>% 
+  left_join(grain_dists, by = c("season", "param", "value")) %>% 
+  mutate(pinkfoot = pf_families + pf_non_breeders,
+         greylag = gl_families + gl_non_breeders,
+         barnacle = bn_families + bn_non_breeders,
+         day = day - 365L*season) %>% 
+  select(-ends_with("families"), -ends_with("non_breeders")) -> goose_numbers 
+
+# Default numbers
+goose_numbers %>% 
+  filter(value == 1) %>% 
+  select(param, day, pinkfoot, greylag, barnacle, grain_dist) %>% 
+  gather(key = species, value = numbers, -param, -day, -grain_dist) %>% 
+  filter(numbers > 0) %>% 
+  group_by(day, species, grain_dist) %>% 
+  summarise(numbers = mean(numbers)) -> default_numbers
+
+# Mean numbers per day
+goose_numbers %>% 
+  select(grain_dist, param, value, season, day, pinkfoot, barnacle, greylag) %>% 
+  group_by(grain_dist, param, value, day) %>% 
+  summarise_at(.vars = vars(pinkfoot, barnacle, greylag),
+               .funs = funs(mean)) %>% 
+  ungroup() %>% 
+  gather(key = species, value = numbers, -grain_dist, -param, -value, -day) %>% 
+  filter(numbers > 0) -> numbers_per_day
+
+# Read the default numbers and standardize
+default_numbers %>% 
+left_join(numbers_per_day, ., 
+          by = c("grain_dist", "day", "species"),
+          suffix = c("", "_default")) %>% 
+  mutate(standarized = numbers/numbers_default) %>% 
+  filter(day <= spring_cutoff,
+         value <= scaler_cutoff) -> numbers_per_day
+```