Merge branch 'main' into more_plots

JGCRI · Jul 11, 2024 · 50211c3 · 50211c3
2 parents 12f1925 + 7b9ca28
commit 50211c3
Show file tree

Hide file tree

Showing 12 changed files with 658 additions and 13 deletions.
diff --git a/results/14_alpha_ohc_unc_nmse.txt b/results/14_alpha_ohc_unc_nmse.txt
@@ -0,0 +1,46 @@
+parameters	values
+beta	0.732
+q10_rh	1.76
+diff	1.042
+alpha	0.61285349058506
+
+Objective Function Value: 5.85e-05
+Counts: 13
+Counts: 13
+Convergence:  0
+Messages: CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH
+
+CO2 MSE: 27.1
+T MSE:   0.0134
+RMSE:    5.2
+T MSE accounting for unc: 0.00301
+OHC MSE accounting for unc: 29.6
+
+***Key Metrics***
+TCRE: 1.55
+TCR:  2
+
+***Historical Warming and ERF***
+GSAT Warming:              1.04
+Ocean Heat Content Change: 451
+Total Aerosol ERF:         -1.24
+WMGHG ERF:                 3.69
+Methane ERF:               0.546
+
+***Future Warming***
+scenario	start	end	GSAT
+ssp119	2021	2040	0.678
+ssp119	2041	2060	0.721
+ssp119	2081	2100	0.389
+ssp126	2021	2040	0.71
+ssp126	2041	2060	0.947
+ssp126	2081	2100	0.787
+ssp245	2021	2040	0.758
+ssp245	2041	2060	1.27
+ssp245	2081	2100	1.84
+ssp370	2021	2040	0.8
+ssp370	2041	2060	1.51
+ssp370	2081	2100	3.08
+ssp585	2021	2040	0.88
+ssp585	2041	2060	1.77
+ssp585	2081	2100	3.91
diff --git a/results/14_alpha_ohc_unc_nmse_big_box.txt b/results/14_alpha_ohc_unc_nmse_big_box.txt
@@ -0,0 +1,46 @@
+parameters	values
+beta	0.903540335704792
+q10_rh	0.88
+diff	0.806
+alpha	0.459527640832312
+
+Objective Function Value: 4.63e-05
+Counts: 31
+Counts: 31
+Convergence:  0
+Messages: CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH
+
+CO2 MSE: 119
+T MSE:   0.0125
+RMSE:    10.9
+T MSE accounting for unc: 0.00248
+OHC MSE accounting for unc: 0.529
+
+***Key Metrics***
+TCRE: 1.24
+TCR:  2.07
+
+***Historical Warming and ERF***
+GSAT Warming:              1.08
+Ocean Heat Content Change: 413
+Total Aerosol ERF:         -1.24
+WMGHG ERF:                 3.39
+Methane ERF:               0.546
+
+***Future Warming***
+scenario	start	end	GSAT
+ssp119	2021	2040	0.549
+ssp119	2041	2060	0.464
+ssp119	2081	2100	0.0238
+ssp126	2021	2040	0.588
+ssp126	2041	2060	0.702
+ssp126	2081	2100	0.384
+ssp245	2021	2040	0.655
+ssp245	2041	2060	1.06
+ssp245	2081	2100	1.42
+ssp370	2021	2040	0.712
+ssp370	2041	2060	1.34
+ssp370	2081	2100	2.7
+ssp585	2021	2040	0.77
+ssp585	2041	2060	1.57
+ssp585	2081	2100	3.48
diff --git a/results/default_stats.txt b/results/default_stats.txt
@@ -18,10 +18,11 @@ Mean of smooth T, CO2 NMSEs (k = 5 ): 0.000846
 Mean of smooth T, CO2 MSEs (k = 10 ): 2.3
 Mean of smooth T, CO2 NMSEs (k = 10 ): 0.000836
 
-CO2 MSE:         4.58
-T MSE:           0.0358
-T MSE with unc:  0.0118
-OHC MSE with unc: 71.4
+
+CO2 MSE:        4.58
+T MSE:          0.0358
+T MSE w/ unc:   0.0118
+OHC MSE w/ unc: 71.4
 RMSE:    2.14
 
 CO2 NMSE: 2.74e-07

diff --git a/results/only_CO2.txt b/results/only_CO2.txt
@@ -0,0 +1,47 @@
+parameters	values
+beta	0.345007750594335
+q10_rh	2.51138174655301
+diff	1.14341924905439
+S	4.28574861619889
+alpha	2.63322316338636
+
+Objective Function Value: 2.62
+Counts: 163
+Counts: 163
+Convergence:  0
+Messages: CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH
+
+CO2 MSE: 2.62
+T MSE:   0.76
+RMSE:    1.62
+T MSE accounting for unc: 0.587
+OHC MSE accounting for unc: 46500
+
+***Key Metrics***
+TCRE: 2.45
+TCR:  2.48
+
+***Historical Warming and ERF***
+GSAT Warming:              -0.575
+Ocean Heat Content Change: -29.9
+Total Aerosol ERF:         -1.24
+WMGHG ERF:                 3.92
+Methane ERF:               0.532
+
+***Future Warming***
+scenario	start	end	GSAT
+ssp119	2021	2040	1.58
+ssp119	2041	2060	2.34
+ssp119	2081	2100	2.55
+ssp126	2021	2040	1.53
+ssp126	2041	2060	2.53
+ssp126	2081	2100	3.17
+ssp245	2021	2040	1.35
+ssp245	2041	2060	2.39
+ssp245	2081	2100	4.03
+ssp370	2021	2040	1.21
+ssp370	2041	2060	2.2
+ssp370	2081	2100	4.58
+ssp585	2021	2040	1.65
+ssp585	2041	2060	3.18
+ssp585	2081	2100	6.55
diff --git a/results/only_OHC.txt b/results/only_OHC.txt
@@ -0,0 +1,47 @@
+parameters	values
+beta	0.537736488380451
+q10_rh	1.76
+diff	1.26924979875569
+S	2.96192583593624
+alpha	1.07684594297377
+
+Objective Function Value: 0
+Counts: 50
+Counts: 50
+Convergence:  52
+Messages: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
+
+CO2 MSE: 3.9
+T MSE:   0.0437
+RMSE:    1.98
+T MSE accounting for unc: 0.0154
+OHC MSE accounting for unc: 0
+
+***Key Metrics***
+TCRE: 1.59
+TCR:  1.93
+
+***Historical Warming and ERF***
+GSAT Warming:              0.741
+Ocean Heat Content Change: 390
+Total Aerosol ERF:         -1.24
+WMGHG ERF:                 3.88
+Methane ERF:               0.54
+
+***Future Warming***
+scenario	start	end	GSAT
+ssp119	2021	2040	0.818
+ssp119	2041	2060	0.988
+ssp119	2081	2100	0.763
+ssp126	2021	2040	0.832
+ssp126	2041	2060	1.2
+ssp126	2081	2100	1.18
+ssp245	2021	2040	0.831
+ssp245	2041	2060	1.41
+ssp245	2081	2100	2.14
+ssp370	2021	2040	0.837
+ssp370	2041	2060	1.56
+ssp370	2081	2100	3.18
+ssp585	2021	2040	0.98
+ssp585	2041	2060	1.93
+ssp585	2081	2100	4.14
diff --git a/results/only_T.txt b/results/only_T.txt
@@ -0,0 +1,47 @@
+parameters	values
+beta	0.651707169999549
+q10_rh	1.76
+diff	1.042
+S	2.50476947896062
+alpha	0.440119271611972
+
+Objective Function Value: 0.00247
+Counts: 83
+Counts: 83
+Convergence:  52
+Messages: ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
+
+CO2 MSE: 8.45
+T MSE:   0.0127
+RMSE:    2.91
+T MSE accounting for unc: 0.00247
+OHC MSE accounting for unc: 71.1
+
+***Key Metrics***
+TCRE: 1.32
+TCR:  1.76
+
+***Historical Warming and ERF***
+GSAT Warming:              1.07
+Ocean Heat Content Change: 448
+Total Aerosol ERF:         -1.24
+WMGHG ERF:                 3.79
+Methane ERF:               0.547
+
+***Future Warming***
+scenario	start	end	GSAT
+ssp119	2021	2040	0.55
+ssp119	2041	2060	0.537
+ssp119	2081	2100	0.184
+ssp126	2021	2040	0.585
+ssp126	2041	2060	0.749
+ssp126	2081	2100	0.542
+ssp245	2021	2040	0.646
+ssp245	2041	2060	1.07
+ssp245	2081	2100	1.52
+ssp370	2021	2040	0.698
+ssp370	2041	2060	1.32
+ssp370	2081	2100	2.69
+ssp585	2021	2040	0.747
+ssp585	2041	2060	1.52
+ssp585	2081	2100	3.38
diff --git a/scripts/alpha_fix_S_ohc_unc_nmse.R b/scripts/alpha_fix_S_ohc_unc_nmse.R
@@ -0,0 +1,88 @@
+# Script for experiment 14
+# Script to use normalized T, OHC and CO2 MSEs while accounting for T, OHC unc
+# Also includes alpha as param to optimize over (but ECS is fixed at 3)
+# Uses OHC range from Matilda table rather than manuscript script
+# Author: Peter Scully
+# Date: 6/25/24
+
+### Constants and Imports ###
+
+# Importing libraries
+library(hector)
+
+# Setting up file paths
+COMP_DATA_DIR <- file.path(here::here(), "comparison_data")
+SCRIPTS_DIR <- file.path(here::here(), "scripts")
+RESULTS_DIR <- file.path(here::here(), "results")
+
+CO2_PATH <- file.path(COMP_DATA_DIR,
+                      "Supplementary_Table_UoM_GHGConcentrations-1-1-0_annualmeans_v23March2017.csv")
+TEMP_PATH <-
+  file.path(COMP_DATA_DIR,
+            "HadCRUT.5.0.2.0.analysis.summary_series.global.annual.csv")
+
+OHC_PATH <- file.path(COMP_DATA_DIR, "OHC_ensemble_Kuhlbrodt_etal_2022.csv")
+
+
+INI_FILE <- system.file("input/hector_ssp245.ini", package = "hector")
+PARAMS <- c(BETA(), Q10_RH(), DIFFUSIVITY(), AERO_SCALE())
+
+OUTPUT <- file.path(RESULTS_DIR, "14_alpha_ohc_unc_nmse_big_box.txt")
+
+
+source(file.path(SCRIPTS_DIR, "major_functions.R"))
+
+### Getting observational data ###
+co2_data <- get_co2_data(CO2_PATH, include_unc = TRUE)
+temp_data <- get_temp_data(TEMP_PATH, include_unc = TRUE)
+ohc_data <- get_ohc_data(OHC_PATH, include_unc = T)
+obs_data <- rbind(co2_data, temp_data, ohc_data)
+
+### Calling optim ###
+best_pars <- run_optim(obs_data = obs_data,
+                       ini_file = INI_FILE,
+                       params = PARAMS,
+                       lower = c(0,               2.2 - 0.44 * 3, 1.16 - 0.118 * 3, 0),
+                       upper = c(0.5 + 0.232 * 3, 2.2 + 0.44 * 3, 1.16 + 0.118 * 3, 3),
+                       yrs = 1750:2014,
+                       vars = c(GMST(), CONCENTRATIONS_CO2(), HEAT_FLUX()),
+                       error_fn = mean_T_CO2_OHC_nmse_unc,
+                       include_unc = T,
+                       method = "L-BFGS-B",
+                       output_file = OUTPUT)
+
+### Outputting individual MSEs ###
+hector_data <- run_hector(ini_file = INI_FILE, 
+                          params = PARAMS, 
+                          vals = best_pars, 
+                          yrs = 1750:2014, 
+                          vars = c(GMST(), CONCENTRATIONS_CO2(), HEAT_FLUX()))
+
+T_mse <- get_var_mse(obs_data = obs_data, 
+                     hector_data = hector_data, 
+                     var = GMST(), 
+                     yrs = 1850:2014)
+CO2_mse <- get_var_mse(obs_data = obs_data, 
+                       hector_data = hector_data, 
+                       var = CONCENTRATIONS_CO2(), 
+                       yrs = c(1750, 1850:2014))
+T_mse_unc <- get_var_mse_unc(obs_data = obs_data,
+                             hector_data = hector_data,
+                             var = GMST(),
+                             yrs = 1850:2014,
+                             mse_fn = mse_unc)
+OHC_mse_unc <- get_var_mse_unc(obs_data = obs_data, 
+                                hector_data = hector_data, 
+                                var = "OHC", 
+                                yrs = 1957:2014,
+                                mse_fn = mse_unc)
+
+write_metric("CO2 MSE:", CO2_mse, OUTPUT)
+write_metric("T MSE:  ", T_mse, OUTPUT)
+write_metric("RMSE:   ", sqrt(mean(CO2_mse, T_mse)), OUTPUT) # not 100% sure this is how we want to calculate this
+write_metric("T MSE accounting for unc:", T_mse_unc, OUTPUT)
+write_metric("OHC MSE accounting for unc:", OHC_mse_unc, OUTPUT)
+write("", OUTPUT, append = TRUE)
+
+### Outputting table metrics ###
+calc_table_metrics(PARAMS, best_pars, OUTPUT)
diff --git a/scripts/default_stats.R b/scripts/default_stats.R
@@ -5,6 +5,7 @@
 ### Constants and Imports ###
 
 # Importing libraries
+library(assertthat)
 library(hector)
 library(zoo)
 
@@ -121,10 +122,14 @@ CO2_mse <- get_var_mse(obs_data = obs_data,
                        var = CONCENTRATIONS_CO2(), 
                        yrs = c(1750, 1850:2014))
 OHC_mse_unc <- get_var_mse_unc(obs_data = obs_data, 
-                             hector_data = hector_data, 
-                             var = "OHC", 
-                             yrs = 1957:2014,
-                             mse_fn = mse_unc)
+                                hector_data = hector_data, 
+                                var = "OHC", 
+                                yrs = 1957:2014,
+                                mse_fn = mse_unc)
+
+assert_that(T_mse_unc == only_T_mse_unc(obs_data, hector_data))
+assert_that(CO2_mse == only_CO2_mse(obs_data, hector_data))
+assert_that(OHC_mse_unc == only_OHC_mse_unc(obs_data, hector_data))
 
 # Getting NMSEs
 T_nmse <- get_var_mse(obs_data = obs_data, 
@@ -202,10 +207,10 @@ OHC_nmae_unc <- get_var_mse_unc(obs_data = obs_data,
                                yrs = 1957:2014,
                                mse_fn = nmae_unc)
 
-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("OHC MSE with unc:", OHC_mse_unc, OUTPUT)
+write_metric("CO2 MSE:       ", CO2_mse, OUTPUT)
+write_metric("T MSE:         ", T_mse, OUTPUT)
+write_metric("T MSE w/ unc:  ", T_mse_unc, OUTPUT)
+write_metric("OHC MSE w/ unc:", OHC_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)