Commit 0b0b36e6 authored by Turner's avatar Turner
Browse files

plot update

parent 4dcc2867

CF_trend_map.png

+91.8 KiB (1.02 MiB)
Loading image diff...

R/model cross validation.R

0 → 100644
+298 −0
Original line number Diff line number Diff line 
# 

# 

# flow_CF_models_scored |> 

#   filter(flow_type == "c") |>

#   filter(pearson > 0.2, abs(bias) < 0.0005) ->

#   study_models_and_results

# 

# 

# model_data_ready |> 

#   select(COMPLXID, year, CF, CFPc) ->

#   data_ready_loo_val

# 

# 

# data_ready_loo_val |> 

#   filter(COMPLXID %in% study_models_and_results$COMPLXID) |> 

#   split(~COMPLXID) |> 

#   future_map_dfr(function(CF_data){

#     

# 

#     if(anyNA(CF_data$CFPc)) return(tibble())

#     if(any(is.infinite(CF_data$CFPc))) return(tibble())

#     

#     2000:2019 |> 

#       map_dfr(function(leave_out_yr){

#         

#         CF_data |>

#           filter(year != leave_out_yr) -> preds_in

#         CF_data |>

#           filter(year == leave_out_yr) -> preds_out

#         

#         # linear models  (2-p and 1-p)

#         

#         lm(CF ~ CFPc, preds_in) |> coefficients() -> model_in

#         

# 

#         get_TmodX2_rmse <- function(x){

#           p1 <- x[1]

#           p2 <- x[2]

#           preds_in |>

#             mutate(pred = 1 + CFPc -

#                      (p2 + CFPc ^ p1)^(1/p1)) ->

#             pred_added

#           pull(pred_added, pred) -> CF_pred

#           pull(pred_added, CF) -> target

#           return(

#             sqrt(mean((CF_pred - target)^2, na.rm = T))

#           )

#         }

# 

#         x2sol <- optim(par = c(1,1), fn = get_TmodX2_rmse,

#                        lower = c(1,1), upper = c(50, 20),

#                        method = "L-BFGS-B")

# 

#         preds_out |>

#           mutate(lin = model_in[1] + CFPc * model_in[2],

#                  tur = 1 + CFPc - (x2sol$par[2] + CFPc ^ x2sol$par[1])^(1/x2sol$par[1]))

#         

#         

#         

#       })

#     

#   }) -> cross_val_results

# 

# 

# 

# cross_val_results |> 

#   filter(COMPLXID %in% study_models_and_results$COMPLXID) |> 

#   mutate(rank = rank(CF), .by = COMPLXID) |> 

#   #filter(rank %in% c(4:15)) |> 

#   mutate(lin_error = 100 * abs(lin-CF),

#          tur_error = 100 * abs(tur-CF)) |> 

#   summarise(Linear = mean(lin_error), `New CF model` = mean(tur_error), .by = COMPLXID) |> 

#   pivot_longer(-COMPLXID, values_to = "mae", names_to = "model") |> 

#   ggplot(aes(y= mae, x = model)) + geom_jitter(alpha = 0.5, col = "dodgerblue") +

#   geom_boxplot(fill = NA, outlier.shape = NA) +

#   ylim(0, 20) +

#   theme_classic() + labs(title = "All years", x = NULL, y = "Mean Absolute Error in CF (Perc. Points)") +

#   NULL -> plot_all

# 

# cross_val_results |> 

#   filter(COMPLXID %in% study_models_and_results$COMPLXID) |> 

#   mutate(rank = rank(CF), .by = COMPLXID) |> 

#   filter(rank <= 4) |> 

#   mutate(lin_error = 100 * abs(lin-CF),

#          tur_error = 100 * abs(tur-CF)) |> 

#   summarise(Linear = mean(lin_error), `New CF model` = mean(tur_error), .by = COMPLXID) |> 

#   pivot_longer(-COMPLXID, values_to = "mae", names_to = "model") |> 

#   ggplot(aes(y= mae, x = model)) + geom_jitter(alpha = 0.5, col = "dodgerblue") +

#   geom_boxplot(fill = NA, outlier.shape = NA) +

#   ylim(0, 20) +

#   theme_classic() + labs(title = "Dry years", x = NULL, y = "") +

#   NULL -> plot_dry

# 

# cross_val_results |> 

#   filter(COMPLXID %in% study_models_and_results$COMPLXID) |> 

#   mutate(rank = rank(CF), .by = COMPLXID) |> 

#   filter(rank >= 16) |> 

#   mutate(lin_error = 100 * abs(lin-CF),

#          tur_error = 100 * abs(tur-CF)) |> 

#   summarise(Linear = mean(lin_error), `New CF model` = mean(tur_error), .by = COMPLXID) |> 

#   pivot_longer(-COMPLXID, values_to = "mae", names_to = "model") |> 

#   ggplot(aes(y= mae, x = model)) + geom_jitter(alpha = 0.5, col = "dodgerblue") +

#   geom_boxplot( fill = NA, outlier.shape = NA) +

#   ylim(0, 20) +

#   theme_classic() + labs(title = "Wet years", x = NULL, y  = "") +

#   NULL -> plot_wet

# 

# plot_all + plot_dry + plot_wet

# 

# 

# 

# # look at model

# expand.grid(

#   CFp = seq(0.01, to = 2, by = 0.01),

#   phi = c(1.1, 1.5, 2, 3, 5, 10),

#   gamma = 1

# ) |> as_tibble() |> 

#   mutate(CF = 1 + CFp - (gamma + CFp ^ phi)^(1/phi)) |> 

#   ggplot(aes(CFp, CF, group = phi)) +

#   geom_line(col = "hotpink") +

#   annotate("rect", xmin = 0, xmax = 2, ymin = 1, ymax = 1.4,

#            fill = "lightgrey", col = NA) +

#   annotate(geom = "polygon", x = c(0, 0, 1.4), y = c(0, 1.4, 1.4),

#             fill = "lightgrey", col = NA) +

#   geom_abline(slope = 1, linetype = 2) +

#   geom_hline(yintercept = 1, linetype = 2) + 

#   theme_minimal() +

#   annotate("text", label = "ENERGY LIMIT", x = 0.5, y = 0.50, angle = 45,

#            vjust = -0.5, size = 5) +

#   theme_minimal() +

#   annotate("text", label = "CAPACITY LIMIT", x = 1, y = 1, size = 5, vjust = -0.5) +

#   coord_equal() +

#   labs(y = "E / Emax (= CF)", x = "Ep / Emax",

#        title = "    \u03b3 = 1      \u03c6 = {1.1, 1.5, 2, 3, 5, 10}") + ylim(0,1.4) -> phi_flex

# 

# 

# 

# expand.grid(

#   CFp = seq(0.01, to = 2, by = 0.01),

#   gamma = c(1, 1.2, 1.5, 2, 2.5),

#   phi = 3

# ) |> as_tibble() |> 

#   mutate(CF = 1 + CFp - (gamma + CFp ^ phi)^(1/phi)) |> 

#   ggplot(aes(CFp, CF, group = gamma, col = phi)) +

#   geom_line(col = "hotpink") +

#   annotate("rect", xmin = 0, xmax = 2, ymin = 1, ymax = 1.4,

#            fill = "lightgrey", col = NA) +

#   annotate(geom = "polygon", x = c(0, 0, 1.4), y = c(0, 1.4, 1.4),

#            fill = "lightgrey", col = NA) +

#   geom_abline(slope = 1, linetype = 2) +

#   geom_hline(yintercept = 1, linetype = 2) + 

#   theme_minimal() +

#   annotate("text", label = "ENERGY LIMIT", x = 0.5, y = 0.50, angle = 45,

#            vjust = -0.5, size = 5) +

#   theme_minimal() +

#   annotate("text", label = "CAPACITY LIMIT", x = 1, y = 1, size = 5, vjust = -0.5) +

#   #coord_equal() +

#   labs(y = "", x = "Ep / Emax", title =  "    \u03c6 = 3      \u03b3 = {1.0, 1.2, 1.5, 2.0, 2.5}") +

#   ylim(0,1.4) -> gamma_flex

# 

# phi_flex + gamma_flex

# 

# 

# expand.grid(

#   CFp = seq(0.01, to = 2, by = 0.01),

#   phi = seq(0.1, to = 5, by = 0.2)

# ) |> as_tibble() |> 

#   mutate(CF = CFp + (1 - phi) * (1 - exp(-CFp))) |> 

#   ggplot(aes(CFp, CF, group = phi, col = phi)) +

#   geom_line(col = "hotpink") +

#   annotate("rect", xmin = 0, xmax = 2, ymin = 1, ymax = 1.4,

#            fill = "lightgrey", col = NA) +

#   annotate(geom = "polygon", x = c(0, 0, 1.4), y = c(0, 1.4, 1.4),

#            fill = "lightgrey", col = NA) +

#   geom_abline(slope = 1, linetype = 2) +

#   geom_hline(yintercept = 1, linetype = 2) +

#   theme_minimal() +

#   annotate("text", label = "ENERGY LIMIT", x = 0.5, y = 0.50, angle = 45,

#            vjust = -0.5, size = 5) +

#   theme_minimal() +

#   annotate("text", label = "CAPACITY LIMIT", x = 1, y = 1, size = 5, vjust = -0.5) +

#   #coord_equal() +

#   labs(y = "", x = "Ep / Emax", title =  "    \u03c6 = 3      \u03b3 = {1.0, 1.2, 1.5, 2.0, 2.5}") +

#   ylim(0,1.4)

# 



















# 

# 

# 

# 

# 

# get_TmodX2_rmse_ex1 <- function(x){

#   p1 <- x[1]

#   p2 <- x[2]

#   ex1_data |>

#     mutate(pred = 1 + CFPc -

#              (p2 + CFPc ^ p1)^(1/p1)) ->

#     pred_added

#   pull(pred_added, pred) -> CF_pred

#   pull(pred_added, CF) -> target

#   return(

#     sqrt(mean((CF_pred - target)^2, na.rm = T))

#   )

# }

# get_TmodX2_rmse_ex2 <- function(x){

#   p1 <- x[1]

#   p2 <- x[2]

#   ex2_data |>

#     mutate(pred = 1 + CFPc -

#              (p2 + CFPc ^ p1)^(1/p1)) ->

#     pred_added

#   pull(pred_added, pred) -> CF_pred

#   pull(pred_added, CF) -> target

#   return(

#     sqrt(mean((CF_pred - target)^2, na.rm = T))

#   )

# }

# 

# 

# data_ready_loo_val |> 

#   filter(COMPLXID == "hc0215") |> 

#   filter(year %in% 2000:2019) -> ex1_data

# 

# data_ready_loo_val |> 

#   filter(COMPLXID == "hc0046") |> 

#   filter(year %in% 2000:2019) -> ex2_data

# 

# lm(CF ~ CFPc, ex1_data) -> lm_params_ex1

# 

# cfmod_params_ex1 <- optim(par = c(1,1), fn = get_TmodX2_rmse_ex1,

#                lower = c(1,1), upper = c(20, 5),

#                method = "L-BFGS-B")

# 

# tibble(CFPc = seq(0, to = 1.1, 0.01)) |> 

#   mutate(CF = lm_params_ex1$coefficients[1] + lm_params_ex1$coefficients[2] * CFPc) ->

#   lin_ex1_data

# 

# tibble(CFPc = seq(0, to = 1.1, 0.01)) |> 

#   mutate(CF = 1 + CFPc - (cfmod_params_ex1$par[2] + CFPc ^ cfmod_params_ex1$par[1])^(1/cfmod_params_ex1$par[1])) ->

#   CFmod_ex1_data

# 

# lm(CF ~ CFPc, ex2_data) -> lm_params_ex2

# 

# cfmod_params_ex2 <- optim(par = c(1,1), fn = get_TmodX2_rmse_ex2,

#                           lower = c(1,1), upper = c(20, 5),

#                           method = "L-BFGS-B")

# 

# tibble(CFPc = seq(0, to = 1.1, 0.01)) |> 

#   mutate(CF = lm_params_ex2$coefficients[1] + lm_params_ex2$coefficients[2] * CFPc) ->

#   lin_ex2_data

# 

# tibble(CFPc = seq(0, to = 1.1, 0.01)) |> 

#   mutate(CF = 1 + CFPc - (cfmod_params_ex2$par[2] + CFPc ^ cfmod_params_ex2$par[1])^(1/cfmod_params_ex2$par[1])) ->

#   CFmod_ex2_data

# 

# 

# ex1_data |> 

#   ggplot(aes(CFPc, CF)) +

#   annotate("rect", xmin = 0, xmax = 1.1, ymin = 1, ymax = 1.1,

#            fill = "lightgrey", col = NA) +

#   annotate(geom = "polygon", x = c(0, 0, 1.1), y = c(0, 1.1, 1.1),

#            fill = "lightgrey", col = NA) +

#   geom_abline(slope = 1, linetype = 2) +

#   geom_hline(yintercept = 1, linetype = 2) + 

#   geom_point(col = "#0070B9") +

#   theme_minimal() +

#   geom_line(data = lin_ex1_data, col = "black", linetype = 3) +

#   geom_line(data = CFmod_ex1_data, col = "#0070B9") +

#   labs(title = "Chief Joseph, WA", subtitle = "Fitted parameters: \u03c6 = 2.24, \u03b3 = 1.04",

#        x = "Capacity Potential (Ep / Emax)", y = "Capacity Factor (E / Emax)") +

#   ylim(0,1.1) + xlim(0,1.1) +

# 

#   ex2_data |> 

#   ggplot(aes(CFPc, CF)) +

#   annotate("rect", xmin = 0, xmax = 1.1, ymin = 1, ymax = 1.1,

#            fill = "lightgrey", col = NA) +

#   annotate(geom = "polygon", x = c(0, 0, 1.1), y = c(0, 1.1, 1.1),

#            fill = "lightgrey", col = NA) +

#   geom_abline(slope = 1, linetype = 2) +

#   geom_hline(yintercept = 1, linetype = 2) + 

#   geom_point(col = "#0070B9") +

#   theme_minimal() +

#   geom_line(data = lin_ex2_data, col = "black", linetype = 3) +

#   geom_line(data = CFmod_ex2_data, col = "#0070B9") +

#   labs(title = "Shasta, CA", subtitle = "Fitted parameters: \u03c6 = 1.85, \u03b3 = 1.02",

#        x = "Capacity Potential (Ep / Emax)", y = "") +

#   ylim(0,1.1) + xlim(0,1.1) +

#   NULL

# 

#

Rplots.pdf

deleted100644 → 0
−3.58 KiB

File deleted.

View file

_targets/meta/meta

+21 −21

File changed.

Preview size limit exceeded, changes collapsed.