Add methods for summary statistics for Leaderboard

This commit adds functionality to find the best single model, worst
single model, and median though `Leaderboard.find_best_single_model`,
`Leaderboard.find_worst_single_model`, `Leaderboard.median`. For
handling NaNs, all NaNs are converted to positive Inf or negative Inf as
appopriate for finding the best single model or worst single model. For
finding the median, the NaNs are filtered out.

NEWS.md

@@ -363,6 +363,21 @@ ClimaAnalysis.add_unit!(rmse_var, "CliMA", "K")
 ClimaAnalysis.add_unit!(rmse_var, Dict("CliMA" => "K")) # for adding multiple units
 ```
 
+#### Summary statistics
+Comparsion between models can be done using `find_best_single_model`,
+`find_worst_single_model`, and `median`. The functions `find_best_single_model` and
+`find_worst_single_model` default to the category "ANN" (corresponding to the annual mean),
+but any category be considered using the parameter `category_name`. Furthermore, the model's
+root mean squared errors (RMSEs) and name is returned. The function `median` only return the
+model's RMSEs. Any `NaN` that appear in the data is ignored when computing the summary
+statistics. See the example below on how to use this functionality.
+
+```julia rmse_var
+ClimaAnalysis.find_best_single_model(rmse_var, category_name = "DJF")
+ClimaAnalysis.find_worst_single_model(rmse_var, category_name = "DJF")
+ClimaAnalysis.median(rmse_var)
+```
+
 ## Bug fixes
 
 - Increased the default value for `warp_string` to 72.

docs/src/api.md

@@ -88,6 +88,9 @@ Base.setindex!(rmse_var::RMSEVariable, rmse, model_name::String)
 Leaderboard.add_category
 Leaderboard.add_model
 Leaderboard.add_unit!
+Leaderboard.find_best_single_model
+Leaderboard.find_worst_single_model
+Leaderboard.median
 ```
 
 ## Utilities

docs/src/rmse_var.md

@@ -123,9 +123,23 @@ ClimaAnalysis.add_unit!(rmse_var, "CliMA", "K")
 ClimaAnalysis.add_unit!(rmse_var, Dict("CliMA" => "K")) # for adding multiple units
 ```
 
+## Summary statistics
+
+`ClimaAnalysis` provides several functions to compute summary statistics. As of now,
+`ClimaAnalysis` provides methods for find the best single model, the worst single model,
+and the median model.
+
+The functions `find_best_single_model` and `find_worst_single_model` default to the category
+"ANN" (corresponding to the annual mean), but any category can be considered using the
+parameter `category_name`. Furthermore, the model's root mean squared errors (RMSEs) and the
+model's name are returned. The function `median` only returns the median model's RMSEs.
+
+Any `NaN` that appears in the data is ignored when computing the summary statistics.
+
+See the example below using this functionality.
+
 ```@repl rmse_var
-ClimaAnalysis.category_names(rmse_var2)
-ClimaAnalysis.model_names(rmse_var)
-ClimaAnalysis.rmse_units(rmse_var)
-rmse_var[:,:]
+ClimaAnalysis.find_best_single_model(rmse_var, category_name = "DJF")
+ClimaAnalysis.find_worst_single_model(rmse_var, category_name = "DJF")
+ClimaAnalysis.median(rmse_var)
 ```

src/Leaderboard.jl

@@ -12,7 +12,10 @@ export RMSEVariable,
     setindex!,
     add_category,
     add_model,
-    add_unit!
+    add_unit!,
+    find_best_single_model,
+    find_worst_single_model,
+    median
 
 """
     Holding root mean squared errors over multiple categories and models for a single
@@ -506,4 +509,75 @@ function add_unit!(rmse_var::RMSEVariable, model_name2unit::Dict)
     return nothing
 end
 
+"""
+    _unit_check(rmse_var::RMSEVariable)
+
+Return nothing if units are not missing and units are the same across all models. Otherwise,
+return an error.
+"""
+function _unit_check(rmse_var::RMSEVariable)
+    units = values(rmse_var.units) |> collect
+    unit_equal = all(unit -> unit == first(units), units)
+    (!unit_equal || first(units) == "") &&
+        error("Units are not the same across all models or units are missing")
+    return nothing
+end
+
+"""
+    find_best_single_model(rmse_var::RMSEVariable; category_name = "ANN")
+
+Return a tuple of the best single model and the name of the model. Find the best single
+model using the root mean squared errors of the category `category_name`.
+"""
+function find_best_single_model(rmse_var::RMSEVariable; category_name = "ANN")
+    _unit_check(rmse_var)
+    categ_names = category_names(rmse_var)
+    ann_idx = categ_names |> (x -> findfirst(y -> (y == category_name), x))
+    isnothing(ann_idx) &&
+        error("The category $category_name does not exist in $categ_names")
+    rmse_vec = rmse_var[:, ann_idx] |> copy
+    # Replace all NaN with Inf so that we do not get NaN as a result
+    # We do this instead of filtering because if we filter, then we need to keep track of
+    # original indices
+    replace!(rmse_vec, NaN => Inf)
+    _, model_idx = findmin(rmse_vec)
+    mdl_names = model_names(rmse_var)
+    return rmse_var[model_idx, :], mdl_names[model_idx]
+end
+
+"""
+    find_worst_single_model(rmse_var::RMSEVariable; category_name = "ANN")
+
+Return a tuple of the worst single model and the name of the model. Find the worst single
+model using the root mean squared errors of the category `category_name`.
+"""
+function find_worst_single_model(rmse_var::RMSEVariable; category_name = "ANN")
+    _unit_check(rmse_var)
+    categ_names = category_names(rmse_var)
+    ann_idx = categ_names |> (x -> findfirst(y -> (y == category_name), x))
+    isnothing(ann_idx) && error("Annual does not exist in $categ_names")
+    rmse_vec = rmse_var[:, ann_idx] |> copy
+    # Replace all NaN with Inf so that we do not get NaN as a result
+    # We do this instead of filtering because if we filter, then we need to keep track of
+    # original indices
+    replace!(rmse_vec, NaN => -Inf)
+    _, model_idx = findmax(rmse_vec)
+    mdl_names = model_names(rmse_var)
+    return rmse_var[model_idx, :], mdl_names[model_idx]
+end
+
+"""
+    median(rmse_var::RMSEVariable)
+
+Find the median using the root mean squared errors across all categories.
+
+Any `NaN` is ignored in computing the median.
+"""
+function median(rmse_var::RMSEVariable)
+    _unit_check(rmse_var)
+    # Drop dimension so that size is (n,) instead of (1,n) so that it is consistent with the
+    # size of the arrays returned from find_worst_single_model and find_best_single_model
+    return dropdims(nanmedian(rmse_var.RMSEs, dims = 1), dims = 1)
+end
+
 end

test/test_Leaderboard.jl

@@ -258,3 +258,45 @@ end
     @test ClimaAnalysis.rmse_units(rmse_var)["hello1"] == "units1"
     @test ClimaAnalysis.rmse_units(rmse_var)["hello2"] == "units2"
 end
+
+@testset "Finding best, worst, and median model" begin
+    csv_file_path = joinpath(@__DIR__, "sample_data/test_csv.csv")
+    rmse_var = ClimaAnalysis.read_rmses(csv_file_path, "ta")
+    rmse_var[:, :] = [[1.0 2.0 3.0 4.0 5.0]; [6.0 7.0 8.0 9.0 10.0]]
+    ClimaAnalysis.add_unit!(rmse_var, "ACCESS-CM2", "units")
+    ClimaAnalysis.add_unit!(rmse_var, "ACCESS-ESM1-5", "units")
+    val, model_name =
+        ClimaAnalysis.find_best_single_model(rmse_var, category_name = "ANN")
+    @test model_name == "ACCESS-CM2"
+    @test val == [1.0, 2.0, 3.0, 4.0, 5.0]
+    @test val |> size == (5,)
+
+    val, model_name =
+        ClimaAnalysis.find_worst_single_model(rmse_var, category_name = "ANN")
+    @test model_name == "ACCESS-ESM1-5"
+    @test val == [6.0, 7.0, 8.0, 9.0, 10.0]
+    @test val |> size == (5,)
+
+    val = ClimaAnalysis.median(rmse_var)
+    @test val == [7.0, 9.0, 11.0, 13.0, 15.0] ./ 2.0
+    @test val |> size == (5,)
+
+    # Test with NaN in RMSE array
+    rmse_var = ClimaAnalysis.add_model(rmse_var, "for adding NaN")
+    ClimaAnalysis.add_unit!(rmse_var, "for adding NaN", "units")
+    val, model_name =
+        ClimaAnalysis.find_best_single_model(rmse_var, category_name = "ANN")
+    @test model_name == "ACCESS-CM2"
+    @test val == [1.0, 2.0, 3.0, 4.0, 5.0]
+    @test val |> size == (5,)
+
+    val, model_name =
+        ClimaAnalysis.find_worst_single_model(rmse_var, category_name = "ANN")
+    @test model_name == "ACCESS-ESM1-5"
+    @test val == [6.0, 7.0, 8.0, 9.0, 10.0]
+    @test val |> size == (5,)
+
+    val = ClimaAnalysis.median(rmse_var)
+    @test val == [7.0, 9.0, 11.0, 13.0, 15.0] ./ 2.0
+    @test val |> size == (5,)
+end