diff --git a/results/default_stats.txt b/results/default_stats.txt index e7c577a..4f10642 100644 --- a/results/default_stats.txt +++ b/results/default_stats.txt @@ -4,11 +4,25 @@ Unnamed Hector core NA q10_rh 1.76 (unitless) Unnamed Hector core NA diff 2.38 cm2/s Mean of T, CO2 MSEs: 2.33 +Mean of T, CO2 NMSEs: 0.00204 +Mean of T, CO2 NMSEs (with unc): 0.00091 -CO2 MSE: 4.58 -T MSE: 0.0771 +Mean of smooth T, CO2 MSEs (k = 3 ): 2.32 +Mean of smooth T, CO2 NMSEs (k = 3 ): 0.00201 +Mean of smooth T, CO2 MSEs (k = 5 ): 2.32 +Mean of smooth T, CO2 NMSEs (k = 5 ): 0.00201 +Mean of smooth T, CO2 MSEs (k = 10 ): 2.32 +Mean of smooth T, CO2 NMSEs (k = 10 ): 0.00198 + +CO2 MSE: 4.58 +T MSE: 0.0771 +T MSE with unc: 0.0343 RMSE: 2.14 +CO2 NMSE: 2.74e-07 +T NMSE: 0.00409 +T NMSE with unc: 0.00182 + ***Key Metrics*** TCRE: 1.51 TCR: 1.77 diff --git a/results/unc_nmse_big_box.txt b/results/unc_nmse_big_box.txt index 60200cf..c3a7cb8 100644 --- a/results/unc_nmse_big_box.txt +++ b/results/unc_nmse_big_box.txt @@ -40,3 +40,54 @@ ssp370 2081 2100 2.7 ssp585 2021 2040 0.805 ssp585 2041 2060 1.61 ssp585 2081 2100 3.58 +parameters values +beta 1.08393761240082 +q10_rh 3.52 +diff 2 + +Objective Function Value: 0.000642 +Counts: 18 +Counts: 18 + +CO2 MSE: 233 +T MSE: 0.0592 +RMSE: 15.3 +T MSE accounting for unc: 0.0239 + +***Key Metrics*** +TCRE: 1.51 +TCR: 1.82 + +***Historical Warming and ERF*** +GSAT Warming: 0.538 +Ocean Heat Content Change: 347 +Total Aerosol ERF: -1.24 +WMGHG ERF: 3.31 +Methane ERF: 0.536 + +***Future Warming*** +scenario start end GSAT +ssp119 2021 2040 0.651 +ssp119 2041 2060 0.754 +ssp119 2081 2100 0.601 +ssp126 2021 2040 0.666 +ssp126 2041 2060 0.933 +ssp126 2081 2100 0.916 +ssp245 2021 2040 0.672 +ssp245 2041 2060 1.14 +ssp245 2081 2100 1.73 +ssp370 2021 2040 0.682 +ssp370 2041 2060 1.28 +ssp370 2081 2100 2.7 +ssp585 2021 2040 0.805 +ssp585 2041 2060 1.61 +ssp585 2081 2100 3.58 +parameters values +beta 1.0838541250268 +q10_rh 3.52 +diff 2 + +Objective Function Value: 0.000642 +Counts: 20 +Counts: 20 +Convergence: 0 diff --git a/scripts/default_stats.R b/scripts/default_stats.R index 87e45b5..e5f2473 100644 --- a/scripts/default_stats.R +++ b/scripts/default_stats.R @@ -6,6 +6,7 @@ # Importing libraries library(hector) +library(zoo) # Setting up file paths COMP_DATA_DIR <- file.path(here::here(), "comparison_data") @@ -28,8 +29,8 @@ source(file.path(SCRIPTS_DIR, "major_functions.R")) ### Getting observational data ### -co2_data <- get_co2_data(CO2_PATH) -temp_data <- get_temp_data(TEMP_PATH) +co2_data <- get_co2_data(CO2_PATH, include_unc = T) +temp_data <- get_temp_data(TEMP_PATH, include_unc = T) obs_data <- rbind(co2_data, temp_data) @@ -52,12 +53,42 @@ hector_data <- run_hector(ini_file = INI_FILE, params = NULL, vals = best_pars, yrs = 1750:2014, - vars = c(GMST(), CONCENTRATIONS_CO2())) + vars = c(GMST(), CONCENTRATIONS_CO2()), + include_unc = T) +# Outputting overall MSEs write_metric("Mean of T, CO2 MSEs:", mean_T_CO2_mse(obs_data, hector_data), OUTPUT) +write_metric("Mean of T, CO2 NMSEs:", + mean_T_CO2_nmse(obs_data, hector_data), + OUTPUT) +write_metric("Mean of T, CO2 NMSEs (with unc):", + mean_T_CO2_nmse_unc(obs_data, hector_data), + OUTPUT) +write("", OUTPUT, append = TRUE) + +# Outputting smoothed MSEs +for (k in c(3, 5, 10)) { + smoothed_temps <- rollmean(temp_data$value, + k = k, + align = "center", + fill = NA) + smoothed_temp_data <- temp_data + smoothed_temp_data$variable <- "Smooth T" + smoothed_temp_data$value <- smoothed_temps + + new_data <- rbind(co2_data, temp_data, smoothed_temp_data) + + write_metric(paste("Mean of smooth T, CO2 MSEs (k =", k, "):"), + smooth_T_CO2_mse(new_data, hector_data), + OUTPUT) + + write_metric(paste("Mean of smooth T, CO2 NMSEs (k =", k, "):"), + smooth_T_CO2_nmse(new_data, hector_data), + OUTPUT) +} write("", OUTPUT, append = TRUE) # Outputting individual MSEs @@ -65,15 +96,42 @@ T_mse <- get_var_mse(obs_data = obs_data, hector_data = hector_data, var = GMST(), yrs = 1850:2014) +T_mse_unc <- get_var_mse_unc(obs_data = obs_data, + hector_data = hector_data, + var = GMST(), + yrs = c(1850:2014), + mse_fn = mse_unc) CO2_mse <- get_var_mse(obs_data = obs_data, hector_data = hector_data, var = CONCENTRATIONS_CO2(), yrs = c(1750, 1850:2014)) -write_metric("CO2 MSE:", CO2_mse, OUTPUT) -write_metric("T MSE: ", T_mse, OUTPUT) +# Getting NMSEs +T_nmse <- get_var_mse(obs_data = obs_data, + hector_data = hector_data, + var = GMST(), + yrs = 1850:2014, + mse_fn = nmse) +T_nmse_unc <- get_var_mse_unc(obs_data = obs_data, + hector_data = hector_data, + var = GMST(), + yrs = c(1850:2014), + mse_fn = nmse_unc) +CO2_nmse <- get_var_mse(obs_data = obs_data, + hector_data = hector_data, + var = CONCENTRATIONS_CO2(), + yrs = c(1750, 1850:2014), + mse_fn = nmse) + +write_metric("CO2 MSE: ", CO2_mse, OUTPUT) +write_metric("T MSE: ", T_mse, OUTPUT) +write_metric("T MSE with unc:", T_mse_unc, OUTPUT) write_metric("RMSE: ", sqrt(mean(CO2_mse, T_mse)), OUTPUT) # not 100% sure this is how we want to calculate this write("", OUTPUT, append = TRUE) +write_metric("CO2 NMSE:", CO2_nmse, OUTPUT) +write_metric("T NMSE: ", T_nmse, OUTPUT) +write_metric("T NMSE with unc:", T_nmse_unc, OUTPUT) +write("", OUTPUT, append = TRUE) ### Outputting table metrics ### calc_table_metrics(NULL, NULL, OUTPUT) \ No newline at end of file diff --git a/scripts/major_functions.R b/scripts/major_functions.R index f412e17..12c413c 100644 --- a/scripts/major_functions.R +++ b/scripts/major_functions.R @@ -400,6 +400,10 @@ run_optim <- function(obs_data, ini_file, params, par = NULL, sd = NULL, # Output counts write_metric("Counts:", optim_output$counts, output_file) + + # Output convergence + write_metric("Convergence: ", optim_output$convergence, output_file) + write(paste("Messages:", optim_output$message), output_file, append = TRUE) write("", output_file, append = TRUE) return(best_pars) diff --git a/scripts/unc_nmse.R b/scripts/unc_nmse.R index 09f176f..9f2db2e 100644 --- a/scripts/unc_nmse.R +++ b/scripts/unc_nmse.R @@ -36,6 +36,7 @@ obs_data <- rbind(co2_data, temp_data) best_pars <- run_optim(obs_data = obs_data, ini_file = INI_FILE, params = PARAMS, + par = c(0, 2.2 - 0.44*3, 2.6), lower = c(0, 2.2 - 0.44 * 3, 2.3 - 0.1 * 3), upper = c(0.5 + 0.232 * 3, 2.2 + 0.44 * 3, 2.3 + 0.1 * 3), yrs = 1750:2014,