-
Skill scores
-
-
Mean skill scores
-
cvRM %>%
- dcast(season ~ lambda, value.var = c('R2', 'RE', 'CE')) %>%
- roundDT(2)
-
## season R2_0.0 R2_1.0 RE_0.0 RE_1.0 CE_0.0 CE_1.0
-## 1: NJ 0.49 0.46 0.53 0.52 0.45 0.43
-## 2: JO 0.44 0.43 0.38 0.39 0.23 0.28
-## 3: WY 0.46 0.46 0.49 0.49 0.40 0.39
+
+
Skill scores as a function of λ
+
+
Mean skill scores and Wilcoxon test
+
cvMean <- cv[, lapply(.SD, mean), .SDcols = c('R2', 'RE', 'CE'), keyby = .(lambda, season)]
+cvMeanLong <- cvMean %>% melt(id.var = c('lambda', 'season'), variable.name = 'metric')
+
+cv0 <- cvLong[lambda == 0]
+
+cvp <- cvLong[,
+ .(p.value = {
+ s <- .BY$season
+ m <- .BY$metric
+ wilcox.test(value, cv0[season == s & metric == m, value])$p.value
+ }),
+ by = .(lambda, season, metric)]
+cvMeanMerged <- merge(cvMeanLong, cvp, by = c('lambda', 'season', 'metric')) %>%
+ merge(cvMeanLong[lambda == 0, .(season, metric, v0 = value)], by = c('season', 'metric'))
+
Figure 3
+
ggplot(cvMeanMerged) +
+ geom_hline(aes(yintercept = value), cvMeanLong[lambda == 0], colour = 'gray') +
+ geom_linerange(
+ aes(lambda, ymin = value, ymax = v0, colour = season, alpha = p.value < 0.1),
+ size = 0.2) +
+ geom_point(
+ aes(lambda, value, colour = season, alpha = p.value < 0.1),
+ size = 1) +
+ geom_point(
+ aes(lambda, value, colour = season),
+ cvMeanLong[lambda == 0],
+ shape = 'triangle',
+ size = 2,
+ alpha = 1) +
+ facet_grid(
+ vars(metric), vars(season),
+ labeller = labeller(.rows = metricLab, .cols = ssnLabBrief),
+ scales = 'free') +
+ scale_alpha_manual(
+ name = 'Significantly different from \u03bb = 0',
+ values = c(0.4, 1)) +
+ scale_colour_manual(values = dark2) +
+ scale_x_continuous(breaks = seq(0, 3, 0.2), labels = skip_label(5)) +
+ guides(color = FALSE) +
+ labs(x= '\u03bb', y = 'Metric value [-]') +
+ panel_border('black', 0.2) +
+ theme(
+ strip.text.y = element_text(angle = 0),
+ legend.position = 'none')
+
Skill distribution
-
Figure S12
-
DT <- cv %>%
- melt(
- id.vars = c('rep', 'season', 'lambda'),
- measure.vars = c('R2', 'RE', 'CE'),
- variable.name = 'metric',
- value.name = 'value')
-hline <- DT[metric != 'R2',
- .(y = 0),
- by = .(season, metric)]
-ggplot(DT) +
- geom_hline(
- aes(yintercept = y),
- hline,
- size = 0.3,
- colour = 'firebrick') +
- geom_boxplot(
- aes(lambda, value, fill = season),
- width = 0.5,
- alpha = 0.5,
- size = 0.3,
- outlier.size = 0.8,
- show.legend = FALSE) +
- facet_grid(
- vars(metric), vars(season),
- labeller = labeller(
- metric = as_labeller(c(R2 = 'R\u00b2', RE = 'RE', CE = 'CE')),
- season = as_labeller(ssnLabBrief)),
- scales = 'free_y') +
- scale_x_discrete(labels = lambdaLab) +
- scale_fill_manual(values = dark2) +
- labs(x = NULL, y = 'Metric value [-]') +
- theme(
- panel.spacing.y = unit(1.5, 'lines'),
- strip.background = element_rect('gray90', NA),
- axis.line.x = element_blank(),
- axis.ticks.x = element_blank())
-
-
Statistically skillful? We calculate the probability of RE > 0 and of CE > 0.
-
cv[, .(pRE = .SD[RE > 0, .N] / .N,
- pCE = .SD[CE > 0, .N] / .N), keyby = .(lambda, season)]
-
## lambda season pRE pCE
-## 1: 0.0 NJ 1.00 0.92
-## 2: 0.0 JO 0.88 0.72
-## 3: 0.0 WY 1.00 0.92
-## 4: 1.0 NJ 1.00 0.92
-## 5: 1.0 JO 1.00 0.94
-## 6: 1.0 WY 1.00 0.92
-
Model 0’s wet season is not statistically skillful.
+
Figure S11
+
ggplot(cvLong[lambda %in% round(seq(0, 3, 0.2), 1)]) +
+ geom_boxplot(
+ aes(lambda, value, group = lambda, fill = season, alpha = I(lambda / 3)),
+ size = 0.2,
+ outlier.size = 0.4, outlier.alpha = 1) +
+ facet_grid(
+ vars(metric), vars(season),
+ labeller = labeller(.rows = metricLab, .cols = ssnLabBrief),
+ scales = 'free') +
+ scale_x_continuous(breaks = seq(0, 3, 0.2), labels = skip_label(5)) +
+ scale_fill_manual(values = dark2) +
+ panel_border('black', 0.2) +
+ labs(x = '\u03bb', y = 'Metric value [-]') +
+ theme(
+ strip.text.y = element_text(angle = 0),
+ legend.position = 'none')
+
Trajectories
-
# Make strings to merge to plot legend
-lambdaLegendDT <- cvRM[, roundDT(.SD, type = 'char')
- ][, .(lab = paste(lambdaLab(.BY$lambda), R2, RE, CE, sep = ' ')),
- by = .(season = as.character(season), lambda)]
-lambdaLegend <- function(s) {
- function(x) lambdaLegendDT[season == s & lambda %in% x, lab]
-}
-
-# Select Model 1, apply low pass filter
-recFinal <- rec[lambda == '1.0']
-recFinal[, season := factor(season, seasons)]
-recFinal[, Qlpf := dplR::pass.filt(Q, 20, 'low')]
-
Figure 3
-
p <- lapply(seasons, function(s) {
- ggplot(rec[season == s & year %in% 1922:2003]) +
- geom_line(aes(year, Q, colour = lambda, linetype = lambda)) +
- geom_line(aes(year, Qa), p1tar[s], colour = 'gray', size = 0.2) +
- scale_linetype_discrete(
- name = c(' R\u00b2 RE CE '),
- labels = lambdaLegend(s)) +
- scale_colour_manual(
- name = c(' R\u00b2 RE CE '),
- values = c('steelblue', 'darkorange'),
- labels = lambdaLegend(s)) +
- scale_x_continuous(
- breaks = seq(1920, 2000, 10),
- labels = skip_label(2)) +
- labs(x = NULL,
- y = if (s == 'JO') 'Q [million m\u00b3]' else '') +
- facet_wrap(
- vars(season),
- strip.position = 'right',
- labeller = as_labeller(ssnLabBrief)) +
- theme(
- legend.title = element_text(hjust = 1, family = 'Lucida Sans Unicode'),
- legend.text = element_text(family = 'Lucida Sans Unicode'),
- axis.text.x = if (s != 'WY') element_blank(),
- axis.ticks.x = if (s != 'WY') element_blank(),
- axis.line.x = if (s != 'WY') element_blank(),
- strip.background = element_rect('gray95', NA),
- legend.box.margin = unit(c(0, 0, 0, -0.5), 'cm'),
- legend.key.width = unit(1, 'cm'))
-
-})
-
-trajPlotInst <- wrap_plots(p) + plot_layout(ncol = 1)
-
-p <- ggplot(recFinal) +
- geom_hline(
- aes(yintercept = Qm, colour = 'Long-term mean'),
- recFinal[, .(Qm = mean(Q)), by = season],
- size = 0.2) +
- geom_line(aes(year, Q, colour = 'Reconstruction (\u03bb = 1)'), size = 0.2) +
- geom_line(aes(year, Qlpf, colour = '20-year low-pass filter'), size = 0.3) +
- facet_wrap(
- vars(season),
- scales = 'free_y',
- labeller = as_labeller(ssnLabBrief),
- strip.position = 'right',
- ncol = 1) +
- scale_x_continuous(
- breaks = seq(1750, 2000, 10),
- labels = skip_label(5),
- sec.axis = dup_axis()) +
- scale_colour_manual(
- name = NULL,
- values = c('gray', 'black', 'firebrick'),
- breaks = c('Reconstruction (\u03bb = 1)', '20-year low-pass filter', 'Long-term mean'),
- ) +
- labs(x = NULL, y = 'Q [million m\u00b3]') +
- theme(
- plot.margin = margin(t = 20),
- strip.background = element_rect('gray95', NA),
- legend.position = 'top',
- legend.key.width = unit(1.2, 'cm'),
- legend.text = element_text(family = 'Lucida Sans Unicode'))
-
-trajPlotFull <- wrap_plots(p)
-
-layout <- c(
- patchwork::area(1, 1, 30, 19),
- patchwork::area(31, 1, 60, 20)
-)
-
-wrap_plots(trajPlotInst) + wrap_plots(trajPlotFull) +
- plot_layout(design = layout)
-
+
# Make strings to merge to plot legend
+lambdaLegendDT <- cvMean[, season := as.character(season)
+ ][, roundDT(.SD, digits = c(1, 2, 2, 2, 2), type = 'char')
+ ][, .(lab = paste(lambdaLab(.BY$lambda), R2, RE, CE, sep = ' ')), by = .(season, lambda)]
+
+
+lambdaLegend <- function(s) {
+ function(x) lambdaLegendDT[season == s & lambda %in% x, lab]
+}
+
Figure 4
+
# Select Model 1, apply low pass filter
+recFinal <- rec[lambda == 1.2]
+recFinal[, season := factor(season, seasons)]
+recFinal[, Qlpf := dplR::pass.filt(Q, 20, 'low')]
+
+p <- lapply(seasons, function(s) {
+ ggplot(rec[season == s & year %in% 1922:2003 & lambda %in% c(0, 1.2)]) +
+ geom_line(aes(year, Q, colour = to_char_1dp(lambda), linetype = to_char_1dp(lambda))) +
+ geom_line(aes(year, Qa), p1tar[s], colour = 'gray', size = 0.2) +
+ scale_linetype_discrete(
+ name = c(' R\u00b2 RE CE '),
+ labels = lambdaLegend(s)) +
+ scale_colour_manual(
+ name = c(' R\u00b2 RE CE '),
+ values = c('steelblue', 'darkorange'),
+ labels = lambdaLegend(s)) +
+ scale_x_continuous(
+ breaks = seq(1920, 2000, 10),
+ labels = skip_label(2)) +
+ labs(x = NULL,
+ y = if (s == 'JO') 'Q [million m\u00b3]' else '') +
+ facet_wrap(
+ vars(season),
+ strip.position = 'right',
+ labeller = as_labeller(ssnLabBrief)) +
+ theme(
+ legend.title = element_text(hjust = 1, family = 'Lucida Sans Unicode'),
+ legend.text = element_text(family = 'Lucida Sans Unicode'),
+ axis.text.x = if (s != 'WY') element_blank(),
+ axis.ticks.x = if (s != 'WY') element_blank(),
+ axis.line.x = if (s != 'WY') element_blank(),
+ strip.background = element_rect('gray95', NA),
+ legend.box.margin = unit(c(0, 0, 0, -0.5), 'cm'),
+ legend.key.width = unit(1, 'cm'))
+
+})
+
+trajPlotInst <- wrap_plots(p) + plot_layout(ncol = 1)
+
+p <- ggplot(recFinal) +
+ geom_hline(
+ aes(yintercept = Qm, colour = 'Long-term mean'),
+ recFinal[, .(Qm = mean(Q)), by = season],
+ size = 0.2) +
+ geom_line(aes(year, Q, colour = 'Reconstruction (\u03bb = 1.2)'), size = 0.2) +
+ geom_line(aes(year, Qlpf, colour = '20-year low-pass filter'), size = 0.3) +
+ facet_wrap(
+ vars(season),
+ scales = 'free_y',
+ labeller = as_labeller(ssnLabBrief),
+ strip.position = 'right',
+ ncol = 1) +
+ scale_x_continuous(
+ breaks = seq(1750, 2000, 10),
+ labels = skip_label(5),
+ sec.axis = dup_axis()) +
+ scale_colour_manual(
+ name = NULL,
+ values = c('gray', 'black', 'firebrick'),
+ breaks = c('Reconstruction (\u03bb = 1.2)', '20-year low-pass filter', 'Long-term mean'),
+ ) +
+ labs(x = NULL, y = 'Q [million m\u00b3]') +
+ theme(
+ plot.margin = margin(t = 20),
+ strip.background = element_rect('gray95', NA),
+ legend.position = 'top',
+ legend.key.width = unit(1.2, 'cm'),
+ legend.text = element_text(family = 'Lucida Sans Unicode'))
+
+trajPlotFull <- wrap_plots(p)
+
+layout <- c(
+ patchwork::area(1, 1, 30, 19),
+ patchwork::area(31, 1, 60, 20)
+)
+
+wrap_plots(trajPlotInst) + wrap_plots(trajPlotFull) +
+ plot_layout(design = layout)
+
Magnitude of adjustments
-
Figure S13
-
# Amount of adjustment
-adj <- rec[, dcast(.SD, season + year ~ lambda, value.var = 'Q')
- ][, adj := `1.0` - `0.0`]
-# Density of adjustment amount
-adjDens <- adj[, {
- s <- .BY$season
- bw <- fcase(s == 'NJ', 30,
- s == 'JO', 50,
- s == 'WY', 9)
- d <- density(adj, bw = bw, from = min(adj), to = max(adj))
- data.table(x = d$x, y = d$y)
-}, by = season]
-
-trajPlotFull2 <- ggplot(rec) +
- geom_line(aes(year, Q, colour = lambda, linetype = lambda), size = 0.4) +
- scale_linetype_discrete(
- name = NULL,
- label = lambdaLab) +
- scale_colour_manual(
- name = NULL,
- values = c('steelblue', 'darkorange'),
- labels = lambdaLab) +
- labs(x = 'Year',
- y = 'Q [Mm\u00b3]') +
- facet_wrap(vars(season), labeller = as_labeller(ssnLabBrief), ncol = 1, scales = 'free_y') +
- theme(
- legend.text = element_text(family = 'Lucida Sans Unicode'),
- legend.position = 'top',
- legend.key.width = unit(1.5, 'cm'),
- strip.background = element_rect('gray95', NA))
-
-adjDensPlot <- ggplot(adjDens) +
- geom_line(aes(x, y), size = 0.2) +
- facet_wrap(vars(season), scales = 'free', ncol = 1, labeller = as_labeller(ssnLabBrief)) +
- theme(strip.background = element_rect('gray95', NA)) +
- labs(x = 'Adjustment [Mm\u00b3]', y = 'Density')
-
-trajPlotFull2 + adjDensPlot +
- plot_layout(widths = c(3, 1))
-
+
Figure S12
+
# Amount of adjustment
+adj <- rec[lambda %in% c(0, 1.2), dcast(.SD, season + year ~ lambda, value.var = 'Q')
+ ][, adj := `1.2` - `0`]
+# Density of adjustment amount
+adjDens <- adj[, {
+ s <- .BY$season
+ bw <- fcase(s == 'NJ', 30,
+ s == 'JO', 60,
+ s == 'WY', 90)
+ d <- density(adj, bw = bw, from = min(adj), to = max(adj))
+ data.table(x = d$x, y = d$y)
+}, by = season]
+
+trajPlotFull2 <- ggplot(rec[lambda %in% c(0, 1.2)]) +
+ geom_line(aes(year, Q, colour = to_char_1dp(lambda), linetype = to_char_1dp(lambda)), size = 0.4) +
+ scale_linetype_discrete(
+ name = NULL,
+ label = lambdaLab) +
+ scale_colour_manual(
+ name = NULL,
+ values = c('steelblue', 'darkorange'),
+ labels = lambdaLab) +
+ labs(x = 'Year',
+ y = 'Q [Mm\u00b3]') +
+ facet_wrap(vars(season), labeller = as_labeller(ssnLabBrief), ncol = 1, scales = 'free_y') +
+ theme(
+ legend.text = element_text(family = 'Lucida Sans Unicode'),
+ legend.position = 'top',
+ legend.key.width = unit(1.5, 'cm'),
+ strip.background = element_rect('gray95', NA))
+
+adjDensPlot <- ggplot(adjDens) +
+ geom_line(aes(x, y), size = 0.2) +
+ facet_wrap(vars(season), scales = 'free', ncol = 1, labeller = as_labeller(ssnLabBrief)) +
+ theme(strip.background = element_rect('gray95', NA)) +
+ labs(x = 'Adjustment [Mm\u00b3]', y = 'Density')
+
+trajPlotFull2 + adjDensPlot +
+ plot_layout(widths = c(3, 1))
+
+
+
+
Flow ratio
+
Figure 5
+
p1tarWide <- dcast(p1tar, year ~ season, value.var = 'Qa')
+p1tarWide[, rQ := NJ / WY]
+
+deltaQ <- rec[, dcast(.SD, lambda + year ~ season, value.var = 'Q')
+ ][, dQ := NJ + JO - WY
+ ][, rQ := NJ / WY]
+densEst <- deltaQ[, {
+ d <- density(dQ, cut = 0)
+ data.table(x = d$x, y = d$y)
+}, by = lambda]
+Qmean <- p1tar[season == 'WY', mean(Qa)]
+
ggplot(deltaQ[lambda == 1.2], aes(WY, rQ * 100)) +
+ geom_hline(yintercept = 50, colour = 'gray', size = 0.2) +
+ geom_point(aes(colour = 'Reconstructed'), size = 1) +
+ geom_point(aes(colour = 'Observed'), p1tarWide) +
+ geom_text(
+ aes(label = year),
+ p1tarWide[year %in% c(1971, 1973, 1975)],
+ size = 3,
+ nudge_x = 150) +
+ geom_text(
+ aes(label = year),
+ deltaQ[lambda == 1.2 & year == 1815],
+ size = 3,
+ nudge_x = 150) +
+ scale_colour_manual(
+ name = NULL,
+ values = pal) +
+ scale_y_continuous(breaks = seq(20, 60, 5)) +
+ scale_x_continuous(breaks = seq(1000, 4000, 500)) +
+ labs(
+ x = 'Annual streamflow [million m\u00b3]',
+ y = 'Fraction of dry season flow to annual flow [%]') +
+ theme(
+ legend.margin = margin(5, 5, 5, 5),
+ legend.background = element_rect(NA, 'black', 0.2),
+ legend.position = c(0.65, 0.9))
+
Selected inputs
-
Extract selected inputs.
-
mbPool[, prox := fifelse(site %in% crnMeta$site, 'RW', 'OX')]
-mbPool[, site.display := paste0(siteLab(site), ' ', prox, ' (', lag, ')')]
-mbPool[, selected := '\u2713']
-
-tickMat <- dcast(mbPool, lambda + prox + site.display ~ season, value.var = 'selected') %>%
- melt(
- id.var = c('lambda', 'prox', 'site.display'),
- variable.name = 'season',
- value.name = 'selected')
-tickMat[is.na(selected), selected := ' ']
-tickMat[, site.display := factor(site.display,
- levels = sort(unique(site.display), decreasing = TRUE))]
-
-count <- tickMat[selected != ' ', .N, by = .(season, lambda)]
-
Figure 4
-
yLabCol <- function(x) {
- sapply(x, function(xx) {
- ind <- grep('RW', xx)
- if (length(ind) == 0) xx <- paste0(ind, '<span style="color:firebrick">', xx, '</span>')
- xx
- })
-}
-ssnLab2Lines <- function(l) {
- n <- count[lambda == l, N]
- function(s) {
- out <- fcase(
- s == 'NJ', 'Dry<br>season',
- s == 'JO', 'Wet<br>season',
- default = 'Water<br>year')
- paste0(out, '<br>**', n, '**')
- }
-}
-
-p1 <- ggplot(tickMat[lambda == '0.0']) +
- geom_text(
- aes(season, site.display, label = selected, colour = prox),
- family = 'Arial Unicode MS',
- show.legend = FALSE,
- vjust = 0.25) +
- facet_wrap(vars(lambda), labeller = as_labeller(lambdaLab)) +
- scale_x_discrete(labels = ssnLab2Lines('0.0')) +
- scale_y_discrete(labels = yLabCol, drop = FALSE) +
- scale_colour_manual(values = c('firebrick', 'black')) +
- labs(x = NULL, y = NULL) +
- cowplot::panel_border('black', 0.2) +
- theme(
- axis.text.y = element_markdown(),
- axis.text.x = element_markdown(lineheight = 1.5),
- strip.text = element_text(family = 'Lucida Sans Unicode', face = 'bold'),
- axis.ticks = element_blank())
-p2 <- ggplot(tickMat[lambda == '1.0']) +
- geom_text(
- aes(season, site.display, label = selected, colour = prox),
- family = 'Arial Unicode MS',
- show.legend = FALSE,
- vjust = 0.25) +
- facet_wrap(vars(lambda), labeller = as_labeller(lambdaLab)) +
- scale_x_discrete(labels = ssnLab2Lines('1.0')) +
- scale_y_discrete(labels = yLabCol, drop = FALSE) +
- scale_colour_manual(values = c('firebrick', 'black')) +
- labs(x = NULL, y = NULL) +
- cowplot::panel_border('black', 0.2) +
- theme(
- axis.text.y = element_blank(),
- axis.text.x = element_markdown(lineheight = 1.5),
- strip.text = element_text(family = 'Lucida Sans Unicode', face = 'bold'),
- strip.background = element_blank(),
- axis.ticks = element_blank())
-
-p1 + p2
-
+
Extract selected inputs
+
weightDT <- best[lambda %in% c(0, 1.2), {
+ DT <- poolDT[c(ga[[1]]@solution == 1)]
+ DT[, prox := fifelse(site %in% crnMeta$site, 'RW', 'OX')]
+ DT[, site.display := paste0(siteLab(site), ' ', prox, ' (', lag, ')')]
+ DT[, weight := weight]
+ DT[]
+}, by = lambda]
+
+weightDT[, site.display :=
+ factor(site.display,
+ levels = sort(unique(site.display), decreasing = TRUE))]
+
+weightDT[, sameSign := fifelse(weight * rho < 0, 'x', ' ')]
+
+count <-weightDT[, .N, by = .(season, lambda)]
+
Figure 6
+
yLabCol <- function(x) {
+ sapply(x, function(xx) {
+ ind <- grep('RW', xx)
+ if (length(ind) == 0) xx <- paste0(ind, '<span style="color:firebrick">', xx, '</span>')
+ xx
+ })
+}
+ssnLab2Lines <- function(l) {
+ n <- count[lambda == l, N]
+ function(s) {
+ out <- fcase(
+ s == 'NJ', 'Dry<br>season',
+ s == 'JO', 'Wet<br>season',
+ default = 'Water<br>year')
+ paste0(out, '<br>**', n, '**')
+ }
+}
+
+pl <- lapply(c(0, 1.2), function(l) {
+ ggplot(weightDT[lambda == l]) +
+ geom_tile(aes(season, site.display, fill = weight)) +
+ geom_text(
+ aes(season, site.display, label = sameSign),
+ vjust = 0.25,
+ colour = 'gray30') +
+ facet_wrap(vars(lambda), labeller = as_labeller(lambdaLab)) +
+ scale_x_discrete(labels = ssnLab2Lines(l), expand = c(0, 0)) +
+ scale_y_discrete(labels = yLabCol, drop = FALSE) +
+ scale_fill_distiller(
+ palette = 'RdBu',
+ direction = 1,
+ name = 'Weight',
+ limits = abs_range(weightDT$weight)) +
+ labs(x = NULL, y = NULL) +
+ cowplot::panel_border('black', 0.2) +
+ theme(
+ legend.key.height = unit(1, 'cm'),
+ axis.text.y = element_markdown(),
+ axis.text.x = element_markdown(lineheight = 1.5),
+ strip.text = element_text(family = 'Lucida Sans Unicode', face = 'bold'),
+ axis.ticks = element_blank())
+})
+
+for (i in 2:length(pl)) pl[[i]] <- pl[[i]] + theme(axis.text.y = element_blank())
+
+wrap_plots(pl, nrow = 1, guides = 'collect')
+
Annual mass balance
-
Calculate mass difference.
-
deltaQ <- rec[, dcast(.SD, year + lambda ~ season, value.var = 'Q')
- ][, dQ := NJ + JO - WY
- ][, rQ := NJ / WY]
-densEst <- deltaQ[, {
- bw <- if (.BY$lambda == '1.0') 50 else 80
- ct <- if (.BY$lambda == '1.0') 1.5 else 1
- d <- density(dQ, bw = bw, cut = ct)
- data.table(x = d$x, y = d$y)
-}, by = lambda]
-
Probability of ΔQ within 10% mean annual flow (190 Mm3)? Range?
-
Qmean <- p1merge[season == 'WY', mean(Qa)]
-deltaQ[, .(p = .SD[abs(dQ) > Qmean / 10, .N] / .N * 100,
- min = min(dQ),
- max = max(dQ)),
- by = lambda] %>%
- roundDT(0)
-
## lambda p min max
-## 1: 0.0 28 -640 600
-## 2: 1.0 13 -289 426
-
Figure 5
-
cols <- c('#984ea3', '#4daf4a')
-
-p1 <- ggplot(deltaQ) +
- geom_line(aes(year, WY, linetype = 'Water year', colour = 'Water year'), size = 0.4) +
- geom_line(aes(year, JO + NJ, linetype = 'Dry + Wet', colour = 'Dry + Wet'), size = 0.4) +
- facet_wrap(
- vars(lambda),
- ncol = 1,
- labeller = as_labeller(lambdaLab)) +
- scale_colour_manual(
- name = NULL,
- values = cols) +
- scale_linetype_manual(
- name = NULL,
- values = c(2, 1)) +
- labs(x = 'Year', y = 'Q [million m\u00b3]') +
- theme(
- legend.position = 'none',
- strip.text = element_text(family = 'Lucida Sans Unicode'),
- strip.background = element_rect('gray90', NA))
-
-p2 <- ggplot(deltaQ) +
- stat_density(
- aes(y = WY, colour = 'Water year', linetype = 'Water year'),
- bw = 160,
- position = 'identity',
- geom = 'line') +
- stat_density(
- aes(y = JO + NJ, colour = 'Dry + Wet', linetype = 'Dry + Wet'),
- bw = 160,
- position = 'identity',
- geom = 'line') +
- scale_x_continuous(breaks = scales::pretty_breaks(3), expand = c(0.01, 0)) +
- scale_colour_manual(
- name = NULL,
- values = cols) +
- scale_linetype_manual(
- name = NULL,
- values = c(2, 1)) +
- facet_wrap(
- vars(lambda),
- ncol = 1,
- labeller = as_labeller(lambdaLab)) +
- labs(x = 'Density', y = NULL) +
- theme(
- legend.key.width = unit(0.95, 'cm'),
- axis.text.y = element_blank(),
- axis.ticks.y = element_blank(),
- strip.text = element_text(family = 'Lucida Sans Unicode'),
- strip.background = element_rect('gray90', NA))
-
-p3 <- ggplot(deltaQ) +
- geom_line(aes(year, dQ, colour = lambda, linetype = lambda), size = 0.4) +
- scale_y_continuous(breaks = seq(-800, 800, 200), limits = abs_range(densEst$x)) +
- scale_colour_manual(
- name = NULL,
- values = c('steelblue', 'darkorange'),
- labels = lambdaLab) +
- scale_linetype_manual(
- name = NULL,
- labels = lambdaLab,
- values = c(1, 2)) +
- labs(x = 'Year', y = '\u0394Q [million m\u00b3]') +
- theme(legend.position = 'none')
-
-p4 <- ggplot(densEst) +
- geom_area(
- aes(x, y),
- densEst[x %between% c(-Qmean / 10, Qmean / 10) & lambda == 'l=1'],
- fill = 'gray95') +
- geom_line(aes(x, y, colour = lambda, linetype = lambda)) +
- scale_colour_manual(
- name = NULL,
- values = c('steelblue', 'darkorange'),
- labels = lambdaLab) +
- scale_x_continuous(breaks = seq(-800, 800, 200), limits = abs_range(densEst$x)) +
- scale_y_continuous(expand = c(0.01, 0)) +
- scale_linetype_manual(
- name = NULL,
- labels = lambdaLab,
- values = c(1, 2)) +
- coord_flip() +
- labs(x = NULL, y = 'Density') +
- theme(
- legend.key.width = unit(0.95, 'cm'),
- legend.text = element_text(family = 'Lucida Sans Unicode'),
- axis.text.y = element_blank(),
- axis.ticks.y = element_blank())
-
-p1 + p2 + p3 + p4 +
- plot_layout(ncol = 2, widths = c(2, 1), heights = c(2, 1))
-
-
-
-
Seasonal mass balance
-
cols1 <- RColorBrewer::brewer.pal(9, 'Set1')[c(3, 4)]
-cols2 <- RColorBrewer::brewer.pal(9, 'Set1')[c(2, 5)]
-
-p1 <- ggplot(deltaQ) +
- geom_line(aes(year, WY - JO, linetype = 'WY - WS', colour = 'WY - WS'), size = 0.4) +
- geom_line(aes(year, NJ, linetype = 'DS', colour = 'DS'), size = 0.4) +
- facet_wrap(
- vars(lambda),
- ncol = 1,
- labeller = as_labeller(lambdaLab)) +
- scale_colour_manual(
- name = NULL,
- breaks = c('WY - WS', 'DS'),
- values = cols1) +
- scale_linetype_manual(
- name = NULL,
- breaks = c('WY - WS', 'DS'),
- values = c(1, 2)) +
- labs(x = 'Year', y = 'Q [million m\u00b3]') +
- theme(
- legend.position = 'none',
- strip.text = element_text(family = 'Lucida Sans Unicode'),
- strip.background = element_rect('gray90', NA))
-
-p2 <- ggplot(deltaQ) +
- stat_density(
- aes(y = NJ, colour = 'DS', linetype = 'DS'),
- bw = 100,
- position = 'identity',
- geom = 'line') +
- stat_density(
- aes(y = WY - JO, colour = 'WY - WS', linetype = 'WY - WS'),
- bw = 100,
- position = 'identity',
- geom = 'line') +
- scale_x_continuous(breaks = scales::pretty_breaks(2), expand = c(0.01, 0)) +
- scale_colour_manual(
- name = NULL,
- breaks = c('WY - WS', 'DS'),
- values = cols1) +
- scale_linetype_manual(
- name = NULL,
- breaks = c('WY - WS', 'DS'),
- values = c(1, 2)) +
- facet_wrap(
- vars(lambda),
- ncol = 1,
- labeller = as_labeller(lambdaLab)) +
- labs(x = 'Density', y = NULL) +
- theme(
- legend.key.width = unit(0.95, 'cm'),
- axis.text.y = element_blank(),
- axis.ticks.y = element_blank(),
- strip.text = element_text(family = 'Lucida Sans Unicode'),
- strip.background = element_rect('gray90', NA))
-
-p3 <- ggplot(deltaQ) +
- geom_line(aes(year, WY - NJ, linetype = 'WY - DS', colour = 'WY - DS'), size = 0.4) +
- geom_line(aes(year, JO, linetype = 'WS', colour = 'WS'), size = 0.4) +
- facet_wrap(
- vars(lambda),
- ncol = 1,
- labeller = as_labeller(lambdaLab)) +
- scale_colour_manual(
- name = NULL,
- breaks = c('WY - DS', 'WS'),
- values = cols2) +
- scale_linetype_manual(
- name = NULL,
- breaks = c('WY - DS', 'WS'),
- values = c(1, 2)) +
- labs(x = 'Year', y = 'Q [million m\u00b3]') +
- theme(
- legend.position = 'none',
- strip.text = element_text(family = 'Lucida Sans Unicode'),
- strip.background = element_rect('gray90', NA))
-
-p4 <- ggplot(deltaQ) +
- stat_density(
- aes(y = JO, colour = 'WS', linetype = 'WS'),
- bw = 170,
- position = 'identity',
- geom = 'line') +
- stat_density(
- aes(y = WY - NJ, colour = 'WY - DS', linetype = 'WY - DS'),
- bw = 170,
- position = 'identity',
- geom = 'line') +
- scale_x_continuous(breaks = scales::pretty_breaks(3), expand = c(0.01, 0)) +
- scale_colour_manual(
- name = NULL,
- breaks = c('WY - DS', 'WS'),
- values = cols2) +
- scale_linetype_manual(
- name = NULL,
- breaks = c('WY - DS', 'WS'),
- values = c(1, 2)) +
- facet_wrap(
- vars(lambda),
- ncol = 1,
- labeller = as_labeller(lambdaLab)) +
- labs(x = 'Density', y = NULL) +
- theme(
- legend.key.width = unit(0.95, 'cm'),
- axis.text.y = element_blank(),
- axis.ticks.y = element_blank(),
- strip.text = element_text(family = 'Lucida Sans Unicode'),
- strip.background = element_rect('gray90', NA))
-
-p1 + p2 + p3 + p4 +
- plot_layout(ncol = 2, widths = c(2, 1))
-
-
-
-
Flow ratio
-
Figure S14
-
p1mergeWide <- dcast(p1merge, year ~ season, value.var = 'Qa')
-p1mergeWide[, rQ := NJ / WY]
-
-rqDensEstRec <- deltaQ[, {
- bw <- if (.BY$lambda == '1.0') 0.031 else 0.033
- d <- density(rQ, bw = bw)
- data.table(x = d$x, y = d$y)
-}, by = lambda]
-
-rqDensEstRecShort <- deltaQ[year %in% 1922:2005, {
- bw <- if (.BY$lambda == '1.0') 0.031 else 0.033
- d <- density(rQ, bw = bw)
- data.table(x = d$x, y = d$y)
-}, by = lambda]
-
-rqDensEstInst <- p1mergeWide[, {
- d <- density(rQ, bw = 0.045)
- data.table(x = d$x, y = d$y)
- }]
-
ggplot(deltaQ[lambda == '1.0'], aes(WY, rQ * 100)) +
- geom_hline(yintercept = 50, colour = 'gray', size = 0.2) +
- geom_point(aes(colour = 'Reconstructed'), size = 1) +
- geom_point(aes(colour = 'Observed'), p1mergeWide) +
- geom_text(
- aes(label = year),
- rbind(
- p1mergeWide[year %in% c(1971, 1973, 1975)],
- p1mergeWide[order(rQ)][.N]),
- size = 3,
- nudge_x = 120) +
- geom_text(
- aes(label = year),
- deltaQ[lambda == '1.0'][order(rQ)][.N],
- size = 3,
- nudge_x = 150) +
- geom_text(
- aes(label = year),
- deltaQ[lambda == '1.0' & year %in% c(1815, 1819, 1879)],
- size = 3,
- nudge_x = 120) +
- scale_colour_manual(
- name = NULL,
- values = pal) +
- scale_y_continuous(breaks = seq(20, 60, 5)) +
- scale_x_continuous(breaks = seq(1000, 4000, 500)) +
- labs(
- x = 'Annual streamflow [million m\u00b3]',
- y = 'Fraction of dry season flow to annual flow [%]') +
- theme(
- legend.margin = margin(5, 5, 5, 5),
- legend.background = element_rect(NA, 'black', 0.2),
- legend.position = c(0.75, 0.9))
-
+
Figure 7
+
densEst <- deltaQ[, {
+ d <- density(dQ, cut = 0, bw = if (lambda == 0) 120 else 100)
+ data.table(x = d$x, y = d$y)
+}, by = lambda]
+
+densEst[, yScaled := y / max(y) * 3]
+densEst[, l2 := to_char_1dp(lambda)]
+
dQlims <- abs_range(densEst[lambda %in% c(0, 1.2), x])
+p1 <- ggplot(deltaQ[lambda %in% c(0, 1.2)]) +
+ geom_line(aes(year, dQ, colour = as.character(lambda), linetype = as.character(lambda)), size = 0.4) +
+ scale_y_continuous(
+ breaks = seq(-800, 800, 200),
+ limits = dQlims) +
+ scale_colour_manual(
+ name = NULL,
+ values = c('steelblue', 'darkorange'),
+ labels = lambdaLab) +
+ scale_linetype_manual(
+ name = NULL,
+ labels = lambdaLab,
+ values = c(1, 2)) +
+ labs(x = 'Year', y = '\u0394Q [million m\u00b3]') +
+ theme(legend.position = 'none')
+
+p2 <- ggplot(densEst[lambda %in% c(0, 1.2)]) +
+ geom_area(
+ aes(x, y),
+ densEst[x %between% c(-Qmean / 10, Qmean / 10) & lambda == 1.2],
+ fill = 'gray95') +
+ geom_line(aes(x, y, colour = as.character(lambda), linetype = as.character(lambda))) +
+ scale_colour_manual(
+ name = NULL,
+ values = c('steelblue', 'darkorange'),
+ labels = lambdaLab) +
+ scale_x_continuous(
+ breaks = seq(-1000, 1000, 200),
+ limits = dQlims) +
+ scale_y_continuous(expand = expansion(c(0, 0.1), 0), labels = skip_label(2)) +
+ scale_linetype_manual(
+ name = NULL,
+ labels = lambdaLab,
+ values = c(1, 2)) +
+ coord_flip() +
+ labs(x = NULL, y = 'Density') +
+ theme(
+ legend.key.width = unit(0.95, 'cm'),
+ legend.text = element_text(family = 'Lucida Sans Unicode'),
+ axis.text.y = element_blank(),
+ axis.ticks.y = element_blank())
+
+lambdas <- round(seq(0, 3, 0.2), 1)
+
+p3 <- ggplot(densEst[lambda %in% lambdas]) +
+ geom_ribbon(
+ aes(x,
+ ymin = lambda * 3,
+ ymax = yScaled + lambda * 3,
+ group = rev(l2),
+ fill = lambda),
+ colour = 'gray30',
+ size = 0.2) +
+ scale_fill_distiller(
+ name = '\u03bb',
+ palette = 'Oranges', direction = 1) +
+ scale_y_continuous(expand = c(0, 0)) +
+ labs(x = '\u0394Q [million m\u00b3]') +
+ theme(
+ panel.grid.major.x = element_line('gray'),
+ axis.text.y = element_blank(),
+ axis.line.y = element_blank(),
+ axis.ticks.y = element_blank())
+
+p4 <- ggplot(
+ deltaQ[lambda %in% lambdas, .(p = .SD[abs(dQ) > Qmean / 10, .N] / .N * 100,
+ min = min(dQ),
+ max = max(dQ)),
+ by = lambda]) +
+ geom_line(aes(lambda, p)) +
+ geom_point(aes(lambda, p)) +
+ scale_x_continuous(breaks = seq(0, 3, 0.2), labels = skip_label(5)) +
+ scale_y_continuous(limits = c(0, 50)) +
+ labs(x = '\u03bb', y = '\u0394Q within 10% mean flow [%]')
+
+p1 + p2 + p3 + p4 +
+ plot_layout(ncol = 2, widths = c(2, 1), guides = 'collect') +
+ plot_annotation(tag_levels = 'a', tag_suffix = ')')
+
diff --git a/results/Xrwdo.RDS b/results/Xrwdo.RDS
new file mode 100644
index 0000000..86a982d
Binary files /dev/null and b/results/Xrwdo.RDS differ
diff --git a/results/p1merge.csv b/results/p1merge.csv
deleted file mode 100644
index d84f7b1..0000000
--- a/results/p1merge.csv
+++ /dev/null
@@ -1,253 +0,0 @@
-season,year,Qa
-NJ,1922,576.02
-NJ,1923,629.17
-NJ,1924,583.9
-NJ,1925,680.41
-NJ,1926,564.32
-NJ,1927,1099.87
-NJ,1928,675.23
-NJ,1929,559.97
-NJ,1930,575.6
-NJ,1931,431.91
-NJ,1932,335.75
-NJ,1933,649.25
-NJ,1934,632.19
-NJ,1935,600.91
-NJ,1936,844.14
-NJ,1937,570.67
-NJ,1938,1031.37
-NJ,1939,974.85
-NJ,1940,1102.03
-NJ,1941,784.93
-NJ,1942,858.72
-NJ,1943,818.48
-NJ,1944,1018.82
-NJ,1945,1100.03
-NJ,1946,768.53
-NJ,1947,669.69
-NJ,1948,710.89
-NJ,1949,722.97
-NJ,1950,1151.52
-NJ,1951,1087.87
-NJ,1952,803.09
-NJ,1953,904.67
-NJ,1954,998.78
-NJ,1955,647.74
-NJ,1956,589.73
-NJ,1957,570.13
-NJ,1958,366.6
-NJ,1959,343.2
-NJ,1960,357.25
-NJ,1961,589.05
-NJ,1962,599.46
-NJ,1963,370.64
-NJ,1964,998.21
-NJ,1965,649.59
-NJ,1966,805.49
-NJ,1967,566.15
-NJ,1968,548.26
-NJ,1969,487.99
-NJ,1970,833.78
-NJ,1971,989.71
-NJ,1972,818.75
-NJ,1973,1038.21
-NJ,1974,922.9
-NJ,1975,1116.74
-NJ,1976,967.8
-NJ,1977,642.93
-NJ,1978,781.72
-NJ,1979,672.02
-NJ,1980,412.34
-NJ,1981,659.53
-NJ,1982,790.69
-NJ,1983,289.26
-NJ,1984,699.65
-NJ,1985,382.25
-NJ,1986,882.975757593226
-NJ,1987,605.868822678699
-NJ,1988,1062.56701953332
-NJ,1989,561.555787190882
-NJ,1990,451.738375126651
-NJ,1991,482.452743073149
-NJ,1992,363.600814194708
-NJ,1993,687.596765387529
-NJ,1994,916.122501379285
-NJ,1995,440.146579201838
-NJ,1996,483.192606887584
-NJ,1997,407.944509603207
-NJ,1998,463.098168185596
-NJ,1999,527.75356590599
-NJ,2000,711.609266083842
-NJ,2001,688.893034146715
-NJ,2002,609.198229295548
-NJ,2003,928.235241088233
-NJ,2004,483.712152118832
-NJ,2005,537.983028861537
-JO,1922,808.19
-JO,1923,773.27
-JO,1924,1109.73
-JO,1925,886.13
-JO,1926,1051.41
-JO,1927,1168.57
-JO,1928,885.34
-JO,1929,1331.43
-JO,1930,1039.22
-JO,1931,433.72
-JO,1932,1046
-JO,1933,1751.42
-JO,1934,1091.41
-JO,1935,1228.09
-JO,1936,864.44
-JO,1937,1304.73
-JO,1938,1810.43
-JO,1939,1763.69
-JO,1940,1077.93
-JO,1941,1032.49
-JO,1942,1568.95
-JO,1943,1673.96
-JO,1944,883.36
-JO,1945,1535.93
-JO,1946,938.92
-JO,1947,1277.87
-JO,1948,1140.04
-JO,1949,1504.72
-JO,1950,1448.66
-JO,1951,1513.46
-JO,1952,1506.29
-JO,1953,1421.2
-JO,1954,854.22
-JO,1955,1231.11
-JO,1956,1631.23
-JO,1957,976.41
-JO,1958,763.52
-JO,1959,1247
-JO,1960,780.05
-JO,1961,1353.2
-JO,1962,827.7
-JO,1963,1125.19
-JO,1964,1269.56
-JO,1965,923.46
-JO,1966,1003
-JO,1967,1267.16
-JO,1968,961.64
-JO,1969,1212.62
-JO,1970,2272.85
-JO,1971,2842
-JO,1972,1117.64
-JO,1973,3280.17
-JO,1974,1352.51
-JO,1975,2570.49
-JO,1976,975.47
-JO,1977,1137.93
-JO,1978,2161.9
-JO,1979,638.54
-JO,1980,1148.43
-JO,1981,1379.1
-JO,1982,656.49
-JO,1983,718.09
-JO,1984,951.79
-JO,1985,665.59
-JO,1986,964.347155490639
-JO,1987,1218.76824001899
-JO,1988,876.495049808683
-JO,1989,648.165362955593
-JO,1990,599.323723433828
-JO,1991,895.693120347209
-JO,1992,762.713083753879
-JO,1993,854.948586830613
-JO,1994,1686.14043147657
-JO,1995,1486.00845350929
-JO,1996,1173.95841626334
-JO,1997,1310.4754800794
-JO,1998,457.261314552663
-JO,1999,810.498640417454
-JO,2000,523.660164351221
-JO,2001,1209.18159109615
-JO,2002,1161.21876703603
-JO,2003,925.165553819053
-JO,2004,1304.22752854299
-JO,2005,1850.2131561004
-WY,1922,1384.21
-WY,1923,1402.44
-WY,1924,1693.63
-WY,1925,1566.54
-WY,1926,1615.73
-WY,1927,2268.44
-WY,1928,1560.57
-WY,1929,1891.4
-WY,1930,1614.82
-WY,1931,865.63
-WY,1932,1381.75
-WY,1933,2400.67
-WY,1934,1723.6
-WY,1935,1829
-WY,1936,1708.58
-WY,1937,1875.4
-WY,1938,2841.8
-WY,1939,2738.54
-WY,1940,2179.96
-WY,1941,1817.42
-WY,1942,2427.67
-WY,1943,2492.44
-WY,1944,1902.18
-WY,1945,2635.96
-WY,1946,1707.45
-WY,1947,1947.56
-WY,1948,1850.93
-WY,1949,2227.69
-WY,1950,2600.18
-WY,1951,2601.33
-WY,1952,2309.38
-WY,1953,2325.87
-WY,1954,1853
-WY,1955,1878.85
-WY,1956,2220.96
-WY,1957,1546.54
-WY,1958,1130.12
-WY,1959,1590.2
-WY,1960,1137.3
-WY,1961,1942.25
-WY,1962,1427.16
-WY,1963,1495.83
-WY,1964,2267.77
-WY,1965,1573.05
-WY,1966,1808.49
-WY,1967,1833.31
-WY,1968,1509.9
-WY,1969,1700.61
-WY,1970,3106.63
-WY,1971,3831.71
-WY,1972,1936.39
-WY,1973,4318.38
-WY,1974,2275.41
-WY,1975,3687.23
-WY,1976,1943.27
-WY,1977,1780.86
-WY,1978,2943.62
-WY,1979,1310.56
-WY,1980,1560.77
-WY,1981,2038.63
-WY,1982,1447.18
-WY,1983,1007.35
-WY,1984,1651.44
-WY,1985,1047.84
-WY,1986,1847.32291308387
-WY,1987,1824.63706269769
-WY,1988,1939.062069342
-WY,1989,1209.72115014647
-WY,1990,1051.06209856048
-WY,1991,1378.14586342036
-WY,1992,1126.31389794859
-WY,1993,1542.54535221814
-WY,1994,2602.26293285585
-WY,1995,1926.15503271112
-WY,1996,1657.15102315093
-WY,1997,1718.41998968261
-WY,1998,920.359482738259
-WY,1999,1338.25220632344
-WY,2000,1235.26943043506
-WY,2001,1898.07462524286
-WY,2002,1770.41699633158
-WY,2003,1853.40079490729
-WY,2004,1787.93968066182
-WY,2005,2388.19618496193
diff --git a/results/p1tar.RDS b/results/p1tar.RDS
new file mode 100644
index 0000000..9b9211a
Binary files /dev/null and b/results/p1tar.RDS differ
diff --git a/results/poolDT.RDS b/results/poolDT.RDS
new file mode 100644
index 0000000..605e4f6
Binary files /dev/null and b/results/poolDT.RDS differ
diff --git a/results/pop500_gen500_lambda0.RDS b/results/pop500_gen500_lambda0.RDS
deleted file mode 100644
index a1d7bb3..0000000
Binary files a/results/pop500_gen500_lambda0.RDS and /dev/null differ
diff --git a/results/pop500_gen500_lambda00_s100.RDS b/results/pop500_gen500_lambda00_s100.RDS
new file mode 100644
index 0000000..31d21a5
Binary files /dev/null and b/results/pop500_gen500_lambda00_s100.RDS differ
diff --git a/results/pop500_gen500_lambda00_s150.RDS b/results/pop500_gen500_lambda00_s150.RDS
new file mode 100644
index 0000000..3e6647b
Binary files /dev/null and b/results/pop500_gen500_lambda00_s150.RDS differ
diff --git a/results/pop500_gen500_lambda00_s200.RDS b/results/pop500_gen500_lambda00_s200.RDS
new file mode 100644
index 0000000..e4cd3a5
Binary files /dev/null and b/results/pop500_gen500_lambda00_s200.RDS differ
diff --git a/results/pop500_gen500_lambda00_s250.RDS b/results/pop500_gen500_lambda00_s250.RDS
new file mode 100644
index 0000000..447a15e
Binary files /dev/null and b/results/pop500_gen500_lambda00_s250.RDS differ
diff --git a/results/pop500_gen500_lambda00_s300.RDS b/results/pop500_gen500_lambda00_s300.RDS
new file mode 100644
index 0000000..3d5c378
Binary files /dev/null and b/results/pop500_gen500_lambda00_s300.RDS differ
diff --git a/results/pop500_gen500_lambda00_s350.RDS b/results/pop500_gen500_lambda00_s350.RDS
new file mode 100644
index 0000000..9117fd6
Binary files /dev/null and b/results/pop500_gen500_lambda00_s350.RDS differ
diff --git a/results/pop500_gen500_lambda00_s400.RDS b/results/pop500_gen500_lambda00_s400.RDS
new file mode 100644
index 0000000..04fc1aa
Binary files /dev/null and b/results/pop500_gen500_lambda00_s400.RDS differ
diff --git a/results/pop500_gen500_lambda00_s450.RDS b/results/pop500_gen500_lambda00_s450.RDS
new file mode 100644
index 0000000..fc8636e
Binary files /dev/null and b/results/pop500_gen500_lambda00_s450.RDS differ
diff --git a/results/pop500_gen500_lambda02_s100.RDS b/results/pop500_gen500_lambda02_s100.RDS
new file mode 100644
index 0000000..e12f006
Binary files /dev/null and b/results/pop500_gen500_lambda02_s100.RDS differ
diff --git a/results/pop500_gen500_lambda02_s150.RDS b/results/pop500_gen500_lambda02_s150.RDS
new file mode 100644
index 0000000..f4ecebe
Binary files /dev/null and b/results/pop500_gen500_lambda02_s150.RDS differ
diff --git a/results/pop500_gen500_lambda02_s200.RDS b/results/pop500_gen500_lambda02_s200.RDS
new file mode 100644
index 0000000..40fbabb
Binary files /dev/null and b/results/pop500_gen500_lambda02_s200.RDS differ
diff --git a/results/pop500_gen500_lambda02_s250.RDS b/results/pop500_gen500_lambda02_s250.RDS
new file mode 100644
index 0000000..f7eb923
Binary files /dev/null and b/results/pop500_gen500_lambda02_s250.RDS differ
diff --git a/results/pop500_gen500_lambda02_s300.RDS b/results/pop500_gen500_lambda02_s300.RDS
new file mode 100644
index 0000000..c0d2db1
Binary files /dev/null and b/results/pop500_gen500_lambda02_s300.RDS differ
diff --git a/results/pop500_gen500_lambda02_s350.RDS b/results/pop500_gen500_lambda02_s350.RDS
new file mode 100644
index 0000000..bc7ec44
Binary files /dev/null and b/results/pop500_gen500_lambda02_s350.RDS differ
diff --git a/results/pop500_gen500_lambda02_s400.RDS b/results/pop500_gen500_lambda02_s400.RDS
new file mode 100644
index 0000000..756dc3e
Binary files /dev/null and b/results/pop500_gen500_lambda02_s400.RDS differ
diff --git a/results/pop500_gen500_lambda02_s450.RDS b/results/pop500_gen500_lambda02_s450.RDS
new file mode 100644
index 0000000..3b8348f
Binary files /dev/null and b/results/pop500_gen500_lambda02_s450.RDS differ
diff --git a/results/pop500_gen500_lambda04_s100.RDS b/results/pop500_gen500_lambda04_s100.RDS
new file mode 100644
index 0000000..65a7033
Binary files /dev/null and b/results/pop500_gen500_lambda04_s100.RDS differ
diff --git a/results/pop500_gen500_lambda04_s150.RDS b/results/pop500_gen500_lambda04_s150.RDS
new file mode 100644
index 0000000..12ef701
Binary files /dev/null and b/results/pop500_gen500_lambda04_s150.RDS differ
diff --git a/results/pop500_gen500_lambda04_s200.RDS b/results/pop500_gen500_lambda04_s200.RDS
new file mode 100644
index 0000000..fefdcc9
Binary files /dev/null and b/results/pop500_gen500_lambda04_s200.RDS differ
diff --git a/results/pop500_gen500_lambda04_s250.RDS b/results/pop500_gen500_lambda04_s250.RDS
new file mode 100644
index 0000000..ebb693e
Binary files /dev/null and b/results/pop500_gen500_lambda04_s250.RDS differ
diff --git a/results/pop500_gen500_lambda04_s300.RDS b/results/pop500_gen500_lambda04_s300.RDS
new file mode 100644
index 0000000..b756cd7
Binary files /dev/null and b/results/pop500_gen500_lambda04_s300.RDS differ
diff --git a/results/pop500_gen500_lambda04_s350.RDS b/results/pop500_gen500_lambda04_s350.RDS
new file mode 100644
index 0000000..33cf2cb
Binary files /dev/null and b/results/pop500_gen500_lambda04_s350.RDS differ
diff --git a/results/pop500_gen500_lambda04_s400.RDS b/results/pop500_gen500_lambda04_s400.RDS
new file mode 100644
index 0000000..b7f682c
Binary files /dev/null and b/results/pop500_gen500_lambda04_s400.RDS differ
diff --git a/results/pop500_gen500_lambda04_s450.RDS b/results/pop500_gen500_lambda04_s450.RDS
new file mode 100644
index 0000000..7f78e2d
Binary files /dev/null and b/results/pop500_gen500_lambda04_s450.RDS differ
diff --git a/results/pop500_gen500_lambda06_s100.RDS b/results/pop500_gen500_lambda06_s100.RDS
new file mode 100644
index 0000000..feabac3
Binary files /dev/null and b/results/pop500_gen500_lambda06_s100.RDS differ
diff --git a/results/pop500_gen500_lambda06_s150.RDS b/results/pop500_gen500_lambda06_s150.RDS
new file mode 100644
index 0000000..6d1cf0f
Binary files /dev/null and b/results/pop500_gen500_lambda06_s150.RDS differ
diff --git a/results/pop500_gen500_lambda06_s200.RDS b/results/pop500_gen500_lambda06_s200.RDS
new file mode 100644
index 0000000..3453445
Binary files /dev/null and b/results/pop500_gen500_lambda06_s200.RDS differ
diff --git a/results/pop500_gen500_lambda06_s250.RDS b/results/pop500_gen500_lambda06_s250.RDS
new file mode 100644
index 0000000..68ea455
Binary files /dev/null and b/results/pop500_gen500_lambda06_s250.RDS differ
diff --git a/results/pop500_gen500_lambda06_s300.RDS b/results/pop500_gen500_lambda06_s300.RDS
new file mode 100644
index 0000000..1479937
Binary files /dev/null and b/results/pop500_gen500_lambda06_s300.RDS differ
diff --git a/results/pop500_gen500_lambda06_s350.RDS b/results/pop500_gen500_lambda06_s350.RDS
new file mode 100644
index 0000000..8de0605
Binary files /dev/null and b/results/pop500_gen500_lambda06_s350.RDS differ
diff --git a/results/pop500_gen500_lambda06_s400.RDS b/results/pop500_gen500_lambda06_s400.RDS
new file mode 100644
index 0000000..63dbd36
Binary files /dev/null and b/results/pop500_gen500_lambda06_s400.RDS differ
diff --git a/results/pop500_gen500_lambda06_s450.RDS b/results/pop500_gen500_lambda06_s450.RDS
new file mode 100644
index 0000000..4ef98a6
Binary files /dev/null and b/results/pop500_gen500_lambda06_s450.RDS differ
diff --git a/results/pop500_gen500_lambda08_s100.RDS b/results/pop500_gen500_lambda08_s100.RDS
new file mode 100644
index 0000000..b8eff53
Binary files /dev/null and b/results/pop500_gen500_lambda08_s100.RDS differ
diff --git a/results/pop500_gen500_lambda08_s150.RDS b/results/pop500_gen500_lambda08_s150.RDS
new file mode 100644
index 0000000..2e42877
Binary files /dev/null and b/results/pop500_gen500_lambda08_s150.RDS differ
diff --git a/results/pop500_gen500_lambda08_s200.RDS b/results/pop500_gen500_lambda08_s200.RDS
new file mode 100644
index 0000000..1198b52
Binary files /dev/null and b/results/pop500_gen500_lambda08_s200.RDS differ
diff --git a/results/pop500_gen500_lambda08_s250.RDS b/results/pop500_gen500_lambda08_s250.RDS
new file mode 100644
index 0000000..61a523a
Binary files /dev/null and b/results/pop500_gen500_lambda08_s250.RDS differ
diff --git a/results/pop500_gen500_lambda08_s300.RDS b/results/pop500_gen500_lambda08_s300.RDS
new file mode 100644
index 0000000..3d392b8
Binary files /dev/null and b/results/pop500_gen500_lambda08_s300.RDS differ
diff --git a/results/pop500_gen500_lambda08_s350.RDS b/results/pop500_gen500_lambda08_s350.RDS
new file mode 100644
index 0000000..9c79df6
Binary files /dev/null and b/results/pop500_gen500_lambda08_s350.RDS differ
diff --git a/results/pop500_gen500_lambda08_s400.RDS b/results/pop500_gen500_lambda08_s400.RDS
new file mode 100644
index 0000000..4576ea3
Binary files /dev/null and b/results/pop500_gen500_lambda08_s400.RDS differ
diff --git a/results/pop500_gen500_lambda08_s450.RDS b/results/pop500_gen500_lambda08_s450.RDS
new file mode 100644
index 0000000..b6eb728
Binary files /dev/null and b/results/pop500_gen500_lambda08_s450.RDS differ
diff --git a/results/pop500_gen500_lambda1.RDS b/results/pop500_gen500_lambda1.RDS
deleted file mode 100644
index cf18891..0000000
Binary files a/results/pop500_gen500_lambda1.RDS and /dev/null differ
diff --git a/results/pop500_gen500_lambda10_s100.RDS b/results/pop500_gen500_lambda10_s100.RDS
new file mode 100644
index 0000000..5405cff
Binary files /dev/null and b/results/pop500_gen500_lambda10_s100.RDS differ
diff --git a/results/pop500_gen500_lambda10_s150.RDS b/results/pop500_gen500_lambda10_s150.RDS
new file mode 100644
index 0000000..1bba16a
Binary files /dev/null and b/results/pop500_gen500_lambda10_s150.RDS differ
diff --git a/results/pop500_gen500_lambda10_s200.RDS b/results/pop500_gen500_lambda10_s200.RDS
new file mode 100644
index 0000000..79d4cc0
Binary files /dev/null and b/results/pop500_gen500_lambda10_s200.RDS differ
diff --git a/results/pop500_gen500_lambda10_s250.RDS b/results/pop500_gen500_lambda10_s250.RDS
new file mode 100644
index 0000000..0c60162
Binary files /dev/null and b/results/pop500_gen500_lambda10_s250.RDS differ
diff --git a/results/pop500_gen500_lambda10_s300.RDS b/results/pop500_gen500_lambda10_s300.RDS
new file mode 100644
index 0000000..4f72ea2
Binary files /dev/null and b/results/pop500_gen500_lambda10_s300.RDS differ
diff --git a/results/pop500_gen500_lambda10_s350.RDS b/results/pop500_gen500_lambda10_s350.RDS
new file mode 100644
index 0000000..603d500
Binary files /dev/null and b/results/pop500_gen500_lambda10_s350.RDS differ
diff --git a/results/pop500_gen500_lambda10_s400.RDS b/results/pop500_gen500_lambda10_s400.RDS
new file mode 100644
index 0000000..9b8e4ac
Binary files /dev/null and b/results/pop500_gen500_lambda10_s400.RDS differ
diff --git a/results/pop500_gen500_lambda10_s450.RDS b/results/pop500_gen500_lambda10_s450.RDS
new file mode 100644
index 0000000..51e9892
Binary files /dev/null and b/results/pop500_gen500_lambda10_s450.RDS differ
diff --git a/results/pop500_gen500_lambda12_s100.RDS b/results/pop500_gen500_lambda12_s100.RDS
new file mode 100644
index 0000000..790e0d4
Binary files /dev/null and b/results/pop500_gen500_lambda12_s100.RDS differ
diff --git a/results/pop500_gen500_lambda12_s150.RDS b/results/pop500_gen500_lambda12_s150.RDS
new file mode 100644
index 0000000..97433b0
Binary files /dev/null and b/results/pop500_gen500_lambda12_s150.RDS differ
diff --git a/results/pop500_gen500_lambda12_s200.RDS b/results/pop500_gen500_lambda12_s200.RDS
new file mode 100644
index 0000000..9abbfb4
Binary files /dev/null and b/results/pop500_gen500_lambda12_s200.RDS differ
diff --git a/results/pop500_gen500_lambda12_s250.RDS b/results/pop500_gen500_lambda12_s250.RDS
new file mode 100644
index 0000000..ca03f32
Binary files /dev/null and b/results/pop500_gen500_lambda12_s250.RDS differ
diff --git a/results/pop500_gen500_lambda12_s300.RDS b/results/pop500_gen500_lambda12_s300.RDS
new file mode 100644
index 0000000..7161bf3
Binary files /dev/null and b/results/pop500_gen500_lambda12_s300.RDS differ
diff --git a/results/pop500_gen500_lambda12_s350.RDS b/results/pop500_gen500_lambda12_s350.RDS
new file mode 100644
index 0000000..993bde5
Binary files /dev/null and b/results/pop500_gen500_lambda12_s350.RDS differ
diff --git a/results/pop500_gen500_lambda12_s400.RDS b/results/pop500_gen500_lambda12_s400.RDS
new file mode 100644
index 0000000..d943967
Binary files /dev/null and b/results/pop500_gen500_lambda12_s400.RDS differ
diff --git a/results/pop500_gen500_lambda12_s450.RDS b/results/pop500_gen500_lambda12_s450.RDS
new file mode 100644
index 0000000..11a3575
Binary files /dev/null and b/results/pop500_gen500_lambda12_s450.RDS differ
diff --git a/results/pop500_gen500_lambda14_s100.RDS b/results/pop500_gen500_lambda14_s100.RDS
new file mode 100644
index 0000000..33e38ca
Binary files /dev/null and b/results/pop500_gen500_lambda14_s100.RDS differ
diff --git a/results/pop500_gen500_lambda14_s150.RDS b/results/pop500_gen500_lambda14_s150.RDS
new file mode 100644
index 0000000..abeabb6
Binary files /dev/null and b/results/pop500_gen500_lambda14_s150.RDS differ
diff --git a/results/pop500_gen500_lambda14_s200.RDS b/results/pop500_gen500_lambda14_s200.RDS
new file mode 100644
index 0000000..d5146bb
Binary files /dev/null and b/results/pop500_gen500_lambda14_s200.RDS differ
diff --git a/results/pop500_gen500_lambda14_s250.RDS b/results/pop500_gen500_lambda14_s250.RDS
new file mode 100644
index 0000000..d9077b5
Binary files /dev/null and b/results/pop500_gen500_lambda14_s250.RDS differ
diff --git a/results/pop500_gen500_lambda14_s300.RDS b/results/pop500_gen500_lambda14_s300.RDS
new file mode 100644
index 0000000..b05f12b
Binary files /dev/null and b/results/pop500_gen500_lambda14_s300.RDS differ
diff --git a/results/pop500_gen500_lambda14_s350.RDS b/results/pop500_gen500_lambda14_s350.RDS
new file mode 100644
index 0000000..621e650
Binary files /dev/null and b/results/pop500_gen500_lambda14_s350.RDS differ
diff --git a/results/pop500_gen500_lambda14_s400.RDS b/results/pop500_gen500_lambda14_s400.RDS
new file mode 100644
index 0000000..c29dd53
Binary files /dev/null and b/results/pop500_gen500_lambda14_s400.RDS differ
diff --git a/results/pop500_gen500_lambda14_s450.RDS b/results/pop500_gen500_lambda14_s450.RDS
new file mode 100644
index 0000000..7989b22
Binary files /dev/null and b/results/pop500_gen500_lambda14_s450.RDS differ
diff --git a/results/pop500_gen500_lambda16_s100.RDS b/results/pop500_gen500_lambda16_s100.RDS
new file mode 100644
index 0000000..f3dfbb2
Binary files /dev/null and b/results/pop500_gen500_lambda16_s100.RDS differ
diff --git a/results/pop500_gen500_lambda16_s150.RDS b/results/pop500_gen500_lambda16_s150.RDS
new file mode 100644
index 0000000..a4bd4d7
Binary files /dev/null and b/results/pop500_gen500_lambda16_s150.RDS differ
diff --git a/results/pop500_gen500_lambda16_s200.RDS b/results/pop500_gen500_lambda16_s200.RDS
new file mode 100644
index 0000000..78d7618
Binary files /dev/null and b/results/pop500_gen500_lambda16_s200.RDS differ
diff --git a/results/pop500_gen500_lambda16_s250.RDS b/results/pop500_gen500_lambda16_s250.RDS
new file mode 100644
index 0000000..135196d
Binary files /dev/null and b/results/pop500_gen500_lambda16_s250.RDS differ
diff --git a/results/pop500_gen500_lambda16_s300.RDS b/results/pop500_gen500_lambda16_s300.RDS
new file mode 100644
index 0000000..29f1dc3
Binary files /dev/null and b/results/pop500_gen500_lambda16_s300.RDS differ
diff --git a/results/pop500_gen500_lambda16_s350.RDS b/results/pop500_gen500_lambda16_s350.RDS
new file mode 100644
index 0000000..30f7120
Binary files /dev/null and b/results/pop500_gen500_lambda16_s350.RDS differ
diff --git a/results/pop500_gen500_lambda16_s400.RDS b/results/pop500_gen500_lambda16_s400.RDS
new file mode 100644
index 0000000..3fdca82
Binary files /dev/null and b/results/pop500_gen500_lambda16_s400.RDS differ
diff --git a/results/pop500_gen500_lambda16_s450.RDS b/results/pop500_gen500_lambda16_s450.RDS
new file mode 100644
index 0000000..11947b3
Binary files /dev/null and b/results/pop500_gen500_lambda16_s450.RDS differ
diff --git a/results/pop500_gen500_lambda18_s100.RDS b/results/pop500_gen500_lambda18_s100.RDS
new file mode 100644
index 0000000..5102378
Binary files /dev/null and b/results/pop500_gen500_lambda18_s100.RDS differ
diff --git a/results/pop500_gen500_lambda18_s150.RDS b/results/pop500_gen500_lambda18_s150.RDS
new file mode 100644
index 0000000..55637af
Binary files /dev/null and b/results/pop500_gen500_lambda18_s150.RDS differ
diff --git a/results/pop500_gen500_lambda18_s200.RDS b/results/pop500_gen500_lambda18_s200.RDS
new file mode 100644
index 0000000..7ac0574
Binary files /dev/null and b/results/pop500_gen500_lambda18_s200.RDS differ
diff --git a/results/pop500_gen500_lambda18_s250.RDS b/results/pop500_gen500_lambda18_s250.RDS
new file mode 100644
index 0000000..4d8a447
Binary files /dev/null and b/results/pop500_gen500_lambda18_s250.RDS differ
diff --git a/results/pop500_gen500_lambda18_s300.RDS b/results/pop500_gen500_lambda18_s300.RDS
new file mode 100644
index 0000000..8c3518d
Binary files /dev/null and b/results/pop500_gen500_lambda18_s300.RDS differ
diff --git a/results/pop500_gen500_lambda18_s350.RDS b/results/pop500_gen500_lambda18_s350.RDS
new file mode 100644
index 0000000..38ff9e3
Binary files /dev/null and b/results/pop500_gen500_lambda18_s350.RDS differ
diff --git a/results/pop500_gen500_lambda18_s400.RDS b/results/pop500_gen500_lambda18_s400.RDS
new file mode 100644
index 0000000..fbe3c15
Binary files /dev/null and b/results/pop500_gen500_lambda18_s400.RDS differ
diff --git a/results/pop500_gen500_lambda18_s450.RDS b/results/pop500_gen500_lambda18_s450.RDS
new file mode 100644
index 0000000..24dcfd6
Binary files /dev/null and b/results/pop500_gen500_lambda18_s450.RDS differ
diff --git a/results/pop500_gen500_lambda20_s100.RDS b/results/pop500_gen500_lambda20_s100.RDS
new file mode 100644
index 0000000..847e465
Binary files /dev/null and b/results/pop500_gen500_lambda20_s100.RDS differ
diff --git a/results/pop500_gen500_lambda20_s150.RDS b/results/pop500_gen500_lambda20_s150.RDS
new file mode 100644
index 0000000..e746029
Binary files /dev/null and b/results/pop500_gen500_lambda20_s150.RDS differ
diff --git a/results/pop500_gen500_lambda20_s200.RDS b/results/pop500_gen500_lambda20_s200.RDS
new file mode 100644
index 0000000..05e7cf4
Binary files /dev/null and b/results/pop500_gen500_lambda20_s200.RDS differ
diff --git a/results/pop500_gen500_lambda20_s250.RDS b/results/pop500_gen500_lambda20_s250.RDS
new file mode 100644
index 0000000..f68cdaa
Binary files /dev/null and b/results/pop500_gen500_lambda20_s250.RDS differ
diff --git a/results/pop500_gen500_lambda20_s300.RDS b/results/pop500_gen500_lambda20_s300.RDS
new file mode 100644
index 0000000..0f55ab0
Binary files /dev/null and b/results/pop500_gen500_lambda20_s300.RDS differ
diff --git a/results/pop500_gen500_lambda20_s350.RDS b/results/pop500_gen500_lambda20_s350.RDS
new file mode 100644
index 0000000..7c9d41b
Binary files /dev/null and b/results/pop500_gen500_lambda20_s350.RDS differ
diff --git a/results/pop500_gen500_lambda20_s400.RDS b/results/pop500_gen500_lambda20_s400.RDS
new file mode 100644
index 0000000..5ba60d6
Binary files /dev/null and b/results/pop500_gen500_lambda20_s400.RDS differ
diff --git a/results/pop500_gen500_lambda20_s450.RDS b/results/pop500_gen500_lambda20_s450.RDS
new file mode 100644
index 0000000..5fd15f0
Binary files /dev/null and b/results/pop500_gen500_lambda20_s450.RDS differ
diff --git a/results/pop500_gen500_lambda22_s100.RDS b/results/pop500_gen500_lambda22_s100.RDS
new file mode 100644
index 0000000..dd03384
Binary files /dev/null and b/results/pop500_gen500_lambda22_s100.RDS differ
diff --git a/results/pop500_gen500_lambda22_s150.RDS b/results/pop500_gen500_lambda22_s150.RDS
new file mode 100644
index 0000000..fc7bc44
Binary files /dev/null and b/results/pop500_gen500_lambda22_s150.RDS differ
diff --git a/results/pop500_gen500_lambda22_s200.RDS b/results/pop500_gen500_lambda22_s200.RDS
new file mode 100644
index 0000000..662ee41
Binary files /dev/null and b/results/pop500_gen500_lambda22_s200.RDS differ
diff --git a/results/pop500_gen500_lambda22_s250.RDS b/results/pop500_gen500_lambda22_s250.RDS
new file mode 100644
index 0000000..0d3e804
Binary files /dev/null and b/results/pop500_gen500_lambda22_s250.RDS differ
diff --git a/results/pop500_gen500_lambda22_s300.RDS b/results/pop500_gen500_lambda22_s300.RDS
new file mode 100644
index 0000000..ee3f342
Binary files /dev/null and b/results/pop500_gen500_lambda22_s300.RDS differ
diff --git a/results/pop500_gen500_lambda22_s350.RDS b/results/pop500_gen500_lambda22_s350.RDS
new file mode 100644
index 0000000..bd49e22
Binary files /dev/null and b/results/pop500_gen500_lambda22_s350.RDS differ
diff --git a/results/pop500_gen500_lambda22_s400.RDS b/results/pop500_gen500_lambda22_s400.RDS
new file mode 100644
index 0000000..4f05036
Binary files /dev/null and b/results/pop500_gen500_lambda22_s400.RDS differ
diff --git a/results/pop500_gen500_lambda22_s450.RDS b/results/pop500_gen500_lambda22_s450.RDS
new file mode 100644
index 0000000..63a4ead
Binary files /dev/null and b/results/pop500_gen500_lambda22_s450.RDS differ
diff --git a/results/pop500_gen500_lambda24_s100.RDS b/results/pop500_gen500_lambda24_s100.RDS
new file mode 100644
index 0000000..3cb6323
Binary files /dev/null and b/results/pop500_gen500_lambda24_s100.RDS differ
diff --git a/results/pop500_gen500_lambda24_s150.RDS b/results/pop500_gen500_lambda24_s150.RDS
new file mode 100644
index 0000000..4884826
Binary files /dev/null and b/results/pop500_gen500_lambda24_s150.RDS differ
diff --git a/results/pop500_gen500_lambda24_s200.RDS b/results/pop500_gen500_lambda24_s200.RDS
new file mode 100644
index 0000000..2c477d1
Binary files /dev/null and b/results/pop500_gen500_lambda24_s200.RDS differ
diff --git a/results/pop500_gen500_lambda24_s250.RDS b/results/pop500_gen500_lambda24_s250.RDS
new file mode 100644
index 0000000..a635dd3
Binary files /dev/null and b/results/pop500_gen500_lambda24_s250.RDS differ
diff --git a/results/pop500_gen500_lambda24_s300.RDS b/results/pop500_gen500_lambda24_s300.RDS
new file mode 100644
index 0000000..3766a95
Binary files /dev/null and b/results/pop500_gen500_lambda24_s300.RDS differ
diff --git a/results/pop500_gen500_lambda24_s350.RDS b/results/pop500_gen500_lambda24_s350.RDS
new file mode 100644
index 0000000..8bd86c7
Binary files /dev/null and b/results/pop500_gen500_lambda24_s350.RDS differ
diff --git a/results/pop500_gen500_lambda24_s400.RDS b/results/pop500_gen500_lambda24_s400.RDS
new file mode 100644
index 0000000..f2a61c8
Binary files /dev/null and b/results/pop500_gen500_lambda24_s400.RDS differ
diff --git a/results/pop500_gen500_lambda24_s450.RDS b/results/pop500_gen500_lambda24_s450.RDS
new file mode 100644
index 0000000..c9672e4
Binary files /dev/null and b/results/pop500_gen500_lambda24_s450.RDS differ
diff --git a/results/pop500_gen500_lambda26_s100.RDS b/results/pop500_gen500_lambda26_s100.RDS
new file mode 100644
index 0000000..ad055bf
Binary files /dev/null and b/results/pop500_gen500_lambda26_s100.RDS differ
diff --git a/results/pop500_gen500_lambda26_s150.RDS b/results/pop500_gen500_lambda26_s150.RDS
new file mode 100644
index 0000000..addb811
Binary files /dev/null and b/results/pop500_gen500_lambda26_s150.RDS differ
diff --git a/results/pop500_gen500_lambda26_s200.RDS b/results/pop500_gen500_lambda26_s200.RDS
new file mode 100644
index 0000000..d5b210c
Binary files /dev/null and b/results/pop500_gen500_lambda26_s200.RDS differ
diff --git a/results/pop500_gen500_lambda26_s250.RDS b/results/pop500_gen500_lambda26_s250.RDS
new file mode 100644
index 0000000..bdbd1b2
Binary files /dev/null and b/results/pop500_gen500_lambda26_s250.RDS differ
diff --git a/results/pop500_gen500_lambda26_s300.RDS b/results/pop500_gen500_lambda26_s300.RDS
new file mode 100644
index 0000000..5a09b02
Binary files /dev/null and b/results/pop500_gen500_lambda26_s300.RDS differ
diff --git a/results/pop500_gen500_lambda26_s350.RDS b/results/pop500_gen500_lambda26_s350.RDS
new file mode 100644
index 0000000..fc8a9f7
Binary files /dev/null and b/results/pop500_gen500_lambda26_s350.RDS differ
diff --git a/results/pop500_gen500_lambda26_s400.RDS b/results/pop500_gen500_lambda26_s400.RDS
new file mode 100644
index 0000000..96066e0
Binary files /dev/null and b/results/pop500_gen500_lambda26_s400.RDS differ
diff --git a/results/pop500_gen500_lambda26_s450.RDS b/results/pop500_gen500_lambda26_s450.RDS
new file mode 100644
index 0000000..9d62ffb
Binary files /dev/null and b/results/pop500_gen500_lambda26_s450.RDS differ
diff --git a/results/pop500_gen500_lambda28_s100.RDS b/results/pop500_gen500_lambda28_s100.RDS
new file mode 100644
index 0000000..3a55906
Binary files /dev/null and b/results/pop500_gen500_lambda28_s100.RDS differ
diff --git a/results/pop500_gen500_lambda28_s150.RDS b/results/pop500_gen500_lambda28_s150.RDS
new file mode 100644
index 0000000..9764960
Binary files /dev/null and b/results/pop500_gen500_lambda28_s150.RDS differ
diff --git a/results/pop500_gen500_lambda28_s200.RDS b/results/pop500_gen500_lambda28_s200.RDS
new file mode 100644
index 0000000..90d9318
Binary files /dev/null and b/results/pop500_gen500_lambda28_s200.RDS differ
diff --git a/results/pop500_gen500_lambda28_s250.RDS b/results/pop500_gen500_lambda28_s250.RDS
new file mode 100644
index 0000000..b5d7a0e
Binary files /dev/null and b/results/pop500_gen500_lambda28_s250.RDS differ
diff --git a/results/pop500_gen500_lambda28_s300.RDS b/results/pop500_gen500_lambda28_s300.RDS
new file mode 100644
index 0000000..08d0bff
Binary files /dev/null and b/results/pop500_gen500_lambda28_s300.RDS differ
diff --git a/results/pop500_gen500_lambda28_s350.RDS b/results/pop500_gen500_lambda28_s350.RDS
new file mode 100644
index 0000000..ec38ceb
Binary files /dev/null and b/results/pop500_gen500_lambda28_s350.RDS differ
diff --git a/results/pop500_gen500_lambda28_s400.RDS b/results/pop500_gen500_lambda28_s400.RDS
new file mode 100644
index 0000000..6ad5d0c
Binary files /dev/null and b/results/pop500_gen500_lambda28_s400.RDS differ
diff --git a/results/pop500_gen500_lambda28_s450.RDS b/results/pop500_gen500_lambda28_s450.RDS
new file mode 100644
index 0000000..faced9e
Binary files /dev/null and b/results/pop500_gen500_lambda28_s450.RDS differ
diff --git a/results/pop500_gen500_lambda30_s100.RDS b/results/pop500_gen500_lambda30_s100.RDS
new file mode 100644
index 0000000..0ef1a6b
Binary files /dev/null and b/results/pop500_gen500_lambda30_s100.RDS differ
diff --git a/results/pop500_gen500_lambda30_s150.RDS b/results/pop500_gen500_lambda30_s150.RDS
new file mode 100644
index 0000000..77ffe9c
Binary files /dev/null and b/results/pop500_gen500_lambda30_s150.RDS differ
diff --git a/results/pop500_gen500_lambda30_s200.RDS b/results/pop500_gen500_lambda30_s200.RDS
new file mode 100644
index 0000000..6b2de82
Binary files /dev/null and b/results/pop500_gen500_lambda30_s200.RDS differ
diff --git a/results/pop500_gen500_lambda30_s250.RDS b/results/pop500_gen500_lambda30_s250.RDS
new file mode 100644
index 0000000..5d54016
Binary files /dev/null and b/results/pop500_gen500_lambda30_s250.RDS differ
diff --git a/results/pop500_gen500_lambda30_s300.RDS b/results/pop500_gen500_lambda30_s300.RDS
new file mode 100644
index 0000000..b0f31c1
Binary files /dev/null and b/results/pop500_gen500_lambda30_s300.RDS differ
diff --git a/results/pop500_gen500_lambda30_s350.RDS b/results/pop500_gen500_lambda30_s350.RDS
new file mode 100644
index 0000000..4e88c49
Binary files /dev/null and b/results/pop500_gen500_lambda30_s350.RDS differ
diff --git a/results/pop500_gen500_lambda30_s400.RDS b/results/pop500_gen500_lambda30_s400.RDS
new file mode 100644
index 0000000..2dfa48a
Binary files /dev/null and b/results/pop500_gen500_lambda30_s400.RDS differ
diff --git a/results/pop500_gen500_lambda30_s450.RDS b/results/pop500_gen500_lambda30_s450.RDS
new file mode 100644
index 0000000..2850a6c
Binary files /dev/null and b/results/pop500_gen500_lambda30_s450.RDS differ
diff --git a/site-selection-GA.R b/site-selection-GA.R
index 4f3b3a7..5fbebf5 100644
--- a/site-selection-GA.R
+++ b/site-selection-GA.R
@@ -1,145 +1,47 @@
-lambda <- 0
+args <- commandArgs(trailingOnly = TRUE)
+lambdas <- seq(0, 2, 0.2)
+lambda <- lambdas[as.integer(args)]
+
pop <- 500
gen <- 500
-source('R/init.R')
-source('R/correlation_functions.R')
+library(magrittr)
+library(glue)
+library(data.table)
+library(mbr)
library(GA)
library(mbr)
-if (!require(memoise)) install.packages('memoise')
-
-if (!require(mbr)) {
- if (!require(remotes)) install.packages('remotes')
- remotes::install_github('ntthung/mbr')
-}
-
-# Streamflow ------------------------------
-p1merge <- fread('results/p1merge.csv', key = 'season')
-
-## Hinkley D and transformation type
-p1HD <- p1merge[, .(D = round(hinkley(Qa), 2),
- Dlog = round(hinkley(log(Qa)), 2)),
- keyby = season
- ][, trans := fifelse(abs(D) > abs(Dlog), 'log', 'none')]
-
-# Proxies -------------------------------
-
-oxi <- fread('data/oxi.csv', key = 'site')
-crn <- fread('data/crn.csv', key = 'site')
-
-crnMeta <- fread('data/crnMeta.csv', key = 'site')
-oxiMeta <- fread('data/oxiMeta.csv', key = 'site')
-
-# In-filling
-crnWideFull <- crn[year %in% 1748:2005, dcast(.SD, year ~ site, value.var = 'rwi')]
-crnMatFull <- as.matrix(crnWideFull[, -'year'])
-crnMatFilled <- imputePCA(crnMatFull, ncp = 19)$completeObs
-crnWideFilled <- as.data.table(crnMatFilled)
-
-oxiWideFull <- oxi[year %in% 1748:2005, dcast(.SD, year ~ site, value.var = 'do18')]
-oxiMatFull <- as.matrix(oxiWideFull[, -'year'])
-oxiMatFilled <- imputePCA(oxiMatFull, ncp = 3)$completeObs
-oxiWideFilled <- as.data.table(oxiMatFilled)
-
-# Input matrices
-allLags <- -2:2 # Core years: 1750:2003, and shift back or forth
-Xrw <- do.call(cbind, lapply(allLags, function(l) {
- x <- crnMatFilled[3:256 - l, ]
- colnames(x) <- paste0(colnames(x), l)
- x
-}))
-Xdo <- do.call(cbind, lapply(allLags, function(l) {
- x <- oxiMatFilled[3:256 - l, ]
- colnames(x) <- paste0(colnames(x), l)
- x
-}))
-Xrwdo <- cbind(Xrw, Xdo)
+source('R/input_selection_functions.R')
+p1tar <- readRDS('results/p1tar.RDS')
+Xrwdo <- readRDS('results/Xrwdo.RDS')
+poolDT <- readRDS('results/poolDT.RDS')
-# Cross-validation -------------------------
set.seed(24)
cvFolds <- make_Z(1922:2003, nRuns = 50, frac = 0.25, contiguous = TRUE)
-# Labelling ---------------------------------
-oxiLab <- function(s) oxiMeta[s, display]
-crnLab <- function(s) crnMeta[s, display]
+lambdaChar <- if (lambda < 1) paste0('0', lambda * 10) else as.character(lambda * 10)
+
+lapply(seq(100, 400), function(s) {
+
+ siteOptim <- ga(
+ type = 'binary',
+ fitness = memoise::memoise(cv_site_selection),
+ pool = poolDT,
+ Xpool = Xrwdo,
+ instQ = p1tar,
+ cv.folds = cvFolds,
+ start.year = 1750,
+ lambda = lambda,
+ log.trans = 1:3,
+ popSize = pop,
+ maxiter = gen,
+ run = min(c(gen, 100)),
+ parallel = TRUE,
+ monitor = FALSE,
+ nBits = nrow(poolDT),
+ seed = s)
+
+ saveRDS(siteOptim, glue('results/pop{pop}_gen{gen}_lambda{lambdaChar}_s{s}.RDS'))
-lambdaLab <- function(x) sapply(x, function(v) paste0('\u03bb = ', substr(v, 3, 5)))
-ssnLabBrief <- function(x) seasonClass[x, type]
-names(seasons) <- seasons <- c('NJ', 'JO', 'WY')
-
-seasonClass <- fread('data/season_class.csv', key = 'season')
-
-# Streamflow: 1922 to 2003
-# Targets: 1922 to 2003
-QDry <- p1merge['NJ'][1:82]
-QWet <- p1merge['JO'][1:82]
-QAnn <- p1merge['WY'][1:82]
-p1tar <- rbind(QDry, QWet, QAnn)
-p1tar[, season := factor(season, seasons)]
-setkey(p1tar, season)
-
-QMat <- log(as.matrix(dcast(p1tar, year ~ season, value.var = 'Qa')[, -'year']))
-cor_lag <- function(X, Y, l, Xyears, Yyears, alpha) {
- ind <- which(Xyears %in% Yyears) - l
- cor_boot(cbind(X[ind, ], Y),
- 1:ncol(X),
- 1:ncol(Y) + ncol(X),
- c('site', 'season'),
- alpha = alpha)
-}
-set.seed(42)
-rhoCrn <- lapplyrbind(allLags, function(l) {
- DT <- cor_lag(crnMatFilled, QMat, l, 1748:2005, 1922:2003, 0.05)
- DT[, lag := l]
- DT[, alpha := 0.05][]
})
-rhoOxi <- lapplyrbind(allLags, function(l) {
- DT <- cor_lag(oxiMatFilled, QMat, l, 1748:2005, 1922:2003, 0.05)
- DT[, lag := l]
- DT[, alpha := 0.05][]
-})
-
-dt <- CJ(s = factor(seasons, seasons), l = allLags)
-dt[, sl := paste0(s, '(', l, ')')]
-
-rhoCrn[, season_lag := paste0(season, '(', lag, ')')]
-rhoCrn[, season_lag := factor(season_lag, levels = dt$sl)]
-rhoCrn[, season := factor(season, levels = seasons)]
-rhoCrn[, lag := factor(lag, levels = allLags)]
-rhoOxi[, season_lag := paste0(season, '(', lag, ')')]
-rhoOxi[, season_lag := factor(season_lag, levels = dt$sl)]
-rhoOxi[, season := factor(season, levels = seasons)]
-rhoOxi[, lag := factor(lag, levels = allLags)]
-
-poolDT <- rbind(rhoCrn[{signif} & abs(rho0) >= 0.2, .(site, lag), by = season],
- rhoOxi[{signif} & abs(rho0) >= 0.2, .(site, lag), by = season]
- )[, site := paste0(site, lag)]
-setkey(poolDT, season)
-
-Xused <- Xrwdo[, poolDT[, unique(site)]]
-
-lambdaChar <- fcase(
- lambda == 0, '00',
- lambda == 1, '10')
-
-siteOptim <- ga(
- type = 'binary',
- fitness = memoise::memoise(mbr::cv_site_selection),
- pool = poolDT,
- seasons = seasons,
- Xpool = Xused,
- instQ = p1tar,
- cv.folds = cvFolds,
- start.year = 1750,
- lambda = lambda,
- log.trans = 1:3,
- num.targets = 3,
- popSize = pop,
- maxiter = gen,
- run = min(c(gen, 100)),
- parallel = TRUE,
- monitor = FALSE,
- nBits = nrow(poolDT))
-
-saveRDS(siteOptim, glue('results/pop{pop}_gen{gen}_lambda{lambdaChar}_refactored.RDS'))
