nhp_output_reports/generate-validation-figures.R at main · The-Strategy-Unit/nhp_output_reports · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
### code to create tables for validation report:
## comparison of projections
## comparison of upper CI with principal projection
## comparison of estimated mitigated activity
purrr::walk(list.files("R", ".R$", , TRUE, TRUE), source)

scheme_code = "XYZ" # add scheme_code for the scenario here to replace XYZ
# If the scheme has site codes already recorded or if all sites are required then set site_codes=NULL, otherwise set sites manually
scenario_name_1 <- "SOC"
scenario_name_2 <- "OBC"

site_codes = NULL
if (is.null(site_codes)) site_codes <- get_sites(scheme_code)
# site_codes = list( # change each element (each can be NULL to mean 'all')
#   ip  = "R0A66",
#   op  = "R0A66",
#   aae = "R0A66"
# )

result_sets = get_nhp_result_sets()

final_report_ndg1 <- result_sets |>
  dplyr::filter(dataset==scheme_code) |>
  dplyr::filter(run_stage=="final_report_ndg1")

final_report_ndg2 <- result_sets |>
  dplyr::filter(dataset==scheme_code) |>
  dplyr::filter(run_stage=="final_report_ndg2")

validation_report_ndg2 <- result_sets |>
  dplyr::filter(dataset==scheme_code) |>
  dplyr::filter(run_stage=="validation_report_ndg2")
validation_report_ndg3 <- result_sets |>
  dplyr::filter(dataset==scheme_code) |>
  dplyr::filter(run_stage=="validation_report_ndg3")
opening_date_scenario <- result_sets |>
  dplyr::filter(dataset==scheme_code) |>
  dplyr::filter(run_stage=="validation_report_ndg2_opening")

selected_results_list <- list(final_report_ndg1,
                              final_report_ndg2,
                              validation_report_ndg2,
                              validation_report_ndg3,
                              opening_date_scenario)


get_final_run_metadata_special <- function(scheme_code, selected_result_set) {

  scheme_results <- selected_result_set |> dplyr::filter(dataset == scheme_code)

  metadata_secondary <- dplyr::filter(scheme_results)
  metadata_primary <- dplyr::filter(scheme_results)

  dplyr::lst(metadata_secondary, metadata_primary)

}

meta <- purrr::map(
  selected_results_list,
  \(result) get_final_run_metadata_special(scheme_code, result)
)

r_final_report_ndg1 <- meta[[1]]$metadata_primary |>
  dplyr::pull(file) |> get_nhp_results(file = _)

r_final_report_ndg2 <- meta[[2]]$metadata_primary |>
  dplyr::pull(file) |>
  get_nhp_results(file = _) #SOC

r_validation_report_ndg2 <- meta[[3]]$metadata_primary |>
  dplyr::pull(file) |> get_nhp_results(file = _) #OBC

r_validation_report_ndg3 <- meta[[4]]$metadata_primary |>
  dplyr::pull(file) |> get_nhp_results(file = _)

r_opening_date_scenario <- meta[[5]]$metadata_primary |>
  dplyr::pull(file) |> get_nhp_results(file = _)

# in CAGR calc, assumes this raises to power of forecast period? Need to account for difference if using opening scenario
# time in years from baseline (23/24) to horizon (41/42 for usual)

get_years <- function(scenario){
  years_to_forecast <- as.numeric(
    scenario[["params"]][["end_year"]] - scenario[["params"]][["start_year"]]
  )
}

fc_period_soc <- get_years(r_final_report_ndg2)
fc_period_obc <- get_years(r_validation_report_ndg2)

get_soc_version <- function(scenario){
  soc_version <- scenario[["params"]][["app_version"]]
  soc_version
}

get_soc_major_version <- function(soc_version){
  soc_major_version <- as.numeric(
    stringr::str_extract(
      soc_version,
      "(?<=v)\\d(?=\\.\\d)" # '4' in e.g. v4.1
    ))
  soc_major_version
}

get_numeric_soc_version <- function(soc_version){
  as.numeric(sub(".", "", soc_version))
}

soc_version <- get_soc_version(r_final_report_ndg2)

soc_major_version <- get_soc_major_version(soc_version)
soc_numeric_version <- get_numeric_soc_version(soc_version)

# get the soc obc data
soc_obc_data <- get_soc_obc(r_final_report_ndg2, r_validation_report_ndg2, site_codes)

# get the soc obc table
soc_obc_table <- get_soc_obc_table(soc_obc_data,soc_numeric_version,scenario_name_1,scenario_name_2)

# get the cagr data
cagr_table <- get_validation_cagr_table(r_final_report_ndg2, r_validation_report_ndg2, site_codes,scenario_name_1,scenario_name_2)

#get the total mitigation table
total_miti_table <- get_total_mitigation_table(r_final_report_ndg2, r_validation_report_ndg2, site_codes,scenario_name_1,scenario_name_2)

# get the mitigation data
tpma_impact_table <- get_tpma_impact_table(r_final_report_ndg2, r_validation_report_ndg2, site_codes,scenario_name_1,scenario_name_2)

# get the p90 table
p90_table <-get_p90_table(soc_obc_data,soc_numeric_version,scenario_name_1,scenario_name_2)

# get soc obc obc_opening table for export to excel
save_soc_obc_open_data <- get_soc_obc_open(r_final_report_ndg2, r_validation_report_ndg2, r_opening_date_scenario, site_codes)

# get the bespoke s curve charts
save_bespoke_ecdf_plots <- get_bespoke_ecdf(r_final_report_ndg2, r_validation_report_ndg2, site_codes)

# get the additional risk table
ecdf_vals <- get_bespoke_ecdf_values(r_final_report_ndg2, r_validation_report_ndg2, site_codes)

# get the additional risk table for opening date scenario
ecdf_vals_open <- get_bespoke_ecdf_values(r_final_report_ndg2, r_opening_date_scenario, site_codes)

# get the details of the model runs featured in these outputs
scenarios_used_details <-  tibble::tibble(
  scheme = c(r_final_report_ndg2[["params"]][["dataset"]],r_validation_report_ndg2[["params"]][["dataset"]]),
  scenario_name = c(scenario_name_1,scenario_name_2),
  scenario = c(r_final_report_ndg2[["params"]][["scenario"]],r_validation_report_ndg2[["params"]][["scenario"]]),
  run_stage = c(meta[[2]][[2]][["run_stage"]],meta[[3]][[2]][["run_stage"]]),
  description = c("Final Report NDG2 Scenario", "Validation Report NDG2 Scenario"),
  baseline = c(glue::glue(r_final_report_ndg2[["params"]][["start_year"]],"/",r_final_report_ndg2[["params"]][["start_year"]]+1),
               glue::glue(r_validation_report_ndg2[["params"]][["start_year"]],"/",r_validation_report_ndg2[["params"]][["start_year"]]+1)),
  horizon = c(glue::glue(r_final_report_ndg2[["params"]][["end_year"]],"/",r_final_report_ndg2[["params"]][["end_year"]]+1),
              glue::glue(r_validation_report_ndg2[["params"]][["end_year"]],"/",r_validation_report_ndg2[["params"]][["end_year"]]+1)),
  model_version = c(r_final_report_ndg2[["params"]][["app_version"]],r_validation_report_ndg2[["params"]][["app_version"]])
)