Tried to figure out weird optim behavior

JGCRI · Jun 19, 2024 · 99137cc · 99137cc
1 parent 909ab50
commit 99137cc
Show file tree

Hide file tree

Showing 5 changed files with 135 additions and 7 deletions.
diff --git 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
@@ -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
@@ -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,28 +53,85 @@ 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
 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)
diff --git 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
@@ -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,