From fb7760eeb1171a1e305a7080389ced7311cc6230 Mon Sep 17 00:00:00 2001 From: nperez Date: Fri, 28 Apr 2023 15:14:23 +0200 Subject: [PATCH 1/7] RMSS and clim scores --- modules/Skill/R/CRPS_clim.R | 13 +++++++++++++ modules/Skill/R/RPS_clim.R | 20 ++++++++++++++++++++ modules/Skill/Skill.R | 37 +++++++++++++++++++++++++++++++++++-- 3 files changed, 68 insertions(+), 2 deletions(-) create mode 100644 modules/Skill/R/CRPS_clim.R create mode 100644 modules/Skill/R/RPS_clim.R diff --git a/modules/Skill/R/CRPS_clim.R b/modules/Skill/R/CRPS_clim.R new file mode 100644 index 00000000..36db4e94 --- /dev/null +++ b/modules/Skill/R/CRPS_clim.R @@ -0,0 +1,13 @@ +CRPS_clim <- function(obs, memb_dim ='ensemble'){ + time_dim <- names(dim(obs)) + obs_time_len <- dim(obs)[time_dim] + + ref <- array(data = rep(obs, each = obs_time_len), dim = c(obs_time_len, obs_time_len)) + names(dim(ref)) <- c(time_dim, memb_dim) + # ref: [sdate, memb] + # obs: [sdate] + crps_ref <- s2dv:::.CRPS(exp = ref, obs = obs, time_dim = time_dim, memb_dim = memb_dim, + dat_dim = NULL, Fair = FALSE) + # crps_ref should be [sdate] + return(mean(crps_ref)) +} diff --git a/modules/Skill/R/RPS_clim.R b/modules/Skill/R/RPS_clim.R new file mode 100644 index 00000000..58f9da4e --- /dev/null +++ b/modules/Skill/R/RPS_clim.R @@ -0,0 +1,20 @@ +# RPS version for climatology +RPS_clim <- function(obs, indices_for_clim = NULL, prob_thresholds = c(1/3, 2/3)) { + + if (is.null(indices_for_clim)){ + indices_for_clim <- 1:length(obs) + } + + obs_probs <- .get_probs(data = obs, indices_for_quantiles = indices_for_clim, ## temporarily removed s2dv::: + prob_thresholds = prob_thresholds, weights = NULL, cross.val = T) ## here! + # clim_probs: [bin, sdate] + clim_probs <- c(prob_thresholds[1], diff(prob_thresholds), 1 - prob_thresholds[length(prob_thresholds)]) + clim_probs <- array(clim_probs, dim = dim(obs_probs)) + + # Calculate RPS for each time step + probs_clim_cumsum <- apply(clim_probs, 2, cumsum) + probs_obs_cumsum <- apply(obs_probs, 2, cumsum) + rps_ref <- apply((probs_clim_cumsum - probs_obs_cumsum)^2, 2, sum) + + return(mean(rps_ref)) +} diff --git a/modules/Skill/Skill.R b/modules/Skill/Skill.R index d3f74b13..42a277b2 100644 --- a/modules/Skill/Skill.R +++ b/modules/Skill/Skill.R @@ -111,18 +111,27 @@ compute_skill_metrics <- function(recipe, data) { time_dim = time_dim, memb_dim = memb_dim, Fair = Fair, + cross.val = recipe$Analysis$Workflow$Skill$cross.val, ncores = ncores) skill <- .drop_dims(skill) skill_metrics[[ metric ]] <- skill + skill <- lapply(skill, function(x) { + .drop_dims(x)}) + skill_metrics[[ metric ]] <- skill$rpss + skill_metrics[[ paste0(metric, "_significance") ]] <- skill$sign + rps_clim <- Apply(list(data$obs$data), + target_dims = c(time_dim, memb_dim), + RPS_clim)$output1 + rps_clim <- .drop_dims(rps_clim) + # Ranked Probability Skill Score and Fair version } else if (metric %in% c('rpss', 'frpss')) { skill <- RPSS(data$hcst$data, data$obs$data, time_dim = time_dim, memb_dim = memb_dim, Fair = Fair, + cross.val = recipe$Analysis$Workflow$Skill$cross.val, ncores = ncores) - skill <- lapply(skill, function(x) { - .drop_dims(x)}) skill_metrics[[ metric ]] <- skill$rpss skill_metrics[[ paste0(metric, "_significance") ]] <- skill$sign # Brier Skill Score - 10th percentile @@ -131,6 +140,7 @@ compute_skill_metrics <- function(recipe, data) { time_dim = time_dim, memb_dim = memb_dim, prob_thresholds = 0.1, + cross.val = recipe$Analysis$Workflow$Skill$cross.val, Fair = Fair, ncores = ncores) skill <- lapply(skill, function(x) { @@ -143,6 +153,7 @@ compute_skill_metrics <- function(recipe, data) { time_dim = time_dim, memb_dim = memb_dim, prob_thresholds = 0.9, + cross.val = recipe$Analysis$Workflow$Skill$cross.val, Fair = Fair, ncores = ncores) skill <- lapply(skill, function(x) { @@ -158,6 +169,10 @@ compute_skill_metrics <- function(recipe, data) { ncores = ncores) skill <- .drop_dims(skill) skill_metrics[[ metric ]] <- skill + crps_clim <- Apply(list(data$obs$data), target_dims = time_dim, + fun = CRPS_clim, memb_dim = memb_dim)$output1 + crps_clim <- .drop_dims(crps_clim) + skill_metrics[['crps_clim']] <- crps_clim # CRPSS and FCRPSS } else if (metric %in% c('crpss', 'fcrpss')) { skill <- CRPSS(data$hcst$data, data$obs$data, @@ -169,6 +184,24 @@ compute_skill_metrics <- function(recipe, data) { .drop_dims(x)}) skill_metrics[[ metric ]] <- skill$crpss skill_metrics[[ paste0(metric, "_significance") ]] <- skill$sign + } else if (metric == 'RMS') { + skill <- RMS(exp = data$hcst$data, + obs = data$obs$data, + time_dim = time_dim, dat_dim = NULL, comp_dim = NULL, + limits = NULL, conf = TRUE, alpha = 0.05, ncores = ncores) + skill <- lapply(skill, function(x) { + .drop_dims(x)}) + skill_metrics[[ metric ]] <- skill$rms + #skill_metrics[[ paste0(metric, "_significance") ]] <- NULL #skill$sign + } else if (metric == 'RMSS') { + skill <- RMSS(exp = data$hcst$data, + obs = data$obs$data, ref = NULL, + time_dim = time_dim, dat_dim = NULL, + memb_dim = memb_dim, pval = TRUE, sign = FALSE, alpha = 0.05, + sig_method = 'one-sided Fisher', ncores = NULL) + skill <- .drop_dims(skill) + skill_metrics[[ metric ]] <- skill + skill_metrics[[ paste0(metric, "_significance") ]] <- skill$sign # Mean bias (climatology) } else if (metric == 'mean_bias') { ## TODO: Eliminate option to compute from anomalies -- GitLab From 25f747305c2964fb8cc8ee7804f48d982d98ba7c Mon Sep 17 00:00:00 2001 From: nperez Date: Fri, 28 Apr 2023 15:17:19 +0200 Subject: [PATCH 2/7] drop_dims added back --- modules/Skill/Skill.R | 2 ++ 1 file changed, 2 insertions(+) diff --git a/modules/Skill/Skill.R b/modules/Skill/Skill.R index 42a277b2..9df1498f 100644 --- a/modules/Skill/Skill.R +++ b/modules/Skill/Skill.R @@ -132,6 +132,8 @@ compute_skill_metrics <- function(recipe, data) { Fair = Fair, cross.val = recipe$Analysis$Workflow$Skill$cross.val, ncores = ncores) + skill <- lapply(skill, function(x) { + .drop_dims(x)}) skill_metrics[[ metric ]] <- skill$rpss skill_metrics[[ paste0(metric, "_significance") ]] <- skill$sign # Brier Skill Score - 10th percentile -- GitLab From dd93e0ec4433485b9608a02b745460af9085680a Mon Sep 17 00:00:00 2001 From: nperez Date: Fri, 28 Apr 2023 16:24:05 +0200 Subject: [PATCH 3/7] fix drop_dims and source --- modules/Skill/R/RPS_clim.R | 2 +- modules/Skill/Skill.R | 19 +++++-------------- 2 files changed, 6 insertions(+), 15 deletions(-) diff --git a/modules/Skill/R/RPS_clim.R b/modules/Skill/R/RPS_clim.R index 58f9da4e..e8b6452d 100644 --- a/modules/Skill/R/RPS_clim.R +++ b/modules/Skill/R/RPS_clim.R @@ -5,7 +5,7 @@ RPS_clim <- function(obs, indices_for_clim = NULL, prob_thresholds = c(1/3, 2/3) indices_for_clim <- 1:length(obs) } - obs_probs <- .get_probs(data = obs, indices_for_quantiles = indices_for_clim, ## temporarily removed s2dv::: + obs_probs <- .GetProbs(data = obs, indices_for_quantiles = indices_for_clim, ## temporarily removed s2dv::: prob_thresholds = prob_thresholds, weights = NULL, cross.val = T) ## here! # clim_probs: [bin, sdate] clim_probs <- c(prob_thresholds[1], diff(prob_thresholds), 1 - prob_thresholds[length(prob_thresholds)]) diff --git a/modules/Skill/Skill.R b/modules/Skill/Skill.R index 9df1498f..ea7b0220 100644 --- a/modules/Skill/Skill.R +++ b/modules/Skill/Skill.R @@ -15,6 +15,8 @@ source("modules/Skill/R/tmp/RPS.R") source("modules/Skill/R/tmp/RPSS.R") source("modules/Skill/R/tmp/GetProbs.R") source("modules/Skill/R/tmp/RandomWalkTest.R") +source("modules/SKill/R/RPS_clim.R") +source("modules/Skill/R/CRPS_clim.R") ## TODO: Implement this in the future ## Which parameter are required? # if (!("obs" %in% ls()) || is.null(obs)) { @@ -113,8 +115,6 @@ compute_skill_metrics <- function(recipe, data) { Fair = Fair, cross.val = recipe$Analysis$Workflow$Skill$cross.val, ncores = ncores) - skill <- .drop_dims(skill) - skill_metrics[[ metric ]] <- skill skill <- lapply(skill, function(x) { .drop_dims(x)}) skill_metrics[[ metric ]] <- skill$rpss @@ -123,7 +123,7 @@ compute_skill_metrics <- function(recipe, data) { target_dims = c(time_dim, memb_dim), RPS_clim)$output1 rps_clim <- .drop_dims(rps_clim) - + skill_metrics[['rps_clim']] <- rps_clim # Ranked Probability Skill Score and Fair version } else if (metric %in% c('rpss', 'frpss')) { skill <- RPSS(data$hcst$data, data$obs$data, @@ -187,23 +187,14 @@ compute_skill_metrics <- function(recipe, data) { skill_metrics[[ metric ]] <- skill$crpss skill_metrics[[ paste0(metric, "_significance") ]] <- skill$sign } else if (metric == 'RMS') { + source("https://earth.bsc.es/gitlab/es/s2dv/-/raw/master/R/RMS.R") skill <- RMS(exp = data$hcst$data, obs = data$obs$data, time_dim = time_dim, dat_dim = NULL, comp_dim = NULL, - limits = NULL, conf = TRUE, alpha = 0.05, ncores = ncores) + limits = NULL, conf = FALSE, alpha = 0.05, ncores = ncores) skill <- lapply(skill, function(x) { .drop_dims(x)}) skill_metrics[[ metric ]] <- skill$rms - #skill_metrics[[ paste0(metric, "_significance") ]] <- NULL #skill$sign - } else if (metric == 'RMSS') { - skill <- RMSS(exp = data$hcst$data, - obs = data$obs$data, ref = NULL, - time_dim = time_dim, dat_dim = NULL, - memb_dim = memb_dim, pval = TRUE, sign = FALSE, alpha = 0.05, - sig_method = 'one-sided Fisher', ncores = NULL) - skill <- .drop_dims(skill) - skill_metrics[[ metric ]] <- skill - skill_metrics[[ paste0(metric, "_significance") ]] <- skill$sign # Mean bias (climatology) } else if (metric == 'mean_bias') { ## TODO: Eliminate option to compute from anomalies -- GitLab From e61b64c5bcd494c0dfdb132af81e2238033f87d6 Mon Sep 17 00:00:00 2001 From: nperez Date: Fri, 28 Apr 2023 16:31:59 +0200 Subject: [PATCH 4/7] typo --- modules/Skill/Skill.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modules/Skill/Skill.R b/modules/Skill/Skill.R index ea7b0220..122aa722 100644 --- a/modules/Skill/Skill.R +++ b/modules/Skill/Skill.R @@ -15,7 +15,7 @@ source("modules/Skill/R/tmp/RPS.R") source("modules/Skill/R/tmp/RPSS.R") source("modules/Skill/R/tmp/GetProbs.R") source("modules/Skill/R/tmp/RandomWalkTest.R") -source("modules/SKill/R/RPS_clim.R") +source("modules/Skill/R/RPS_clim.R") source("modules/Skill/R/CRPS_clim.R") ## TODO: Implement this in the future ## Which parameter are required? -- GitLab From cbae842f0a91df5e68a7a386eb6feace699b6a2e Mon Sep 17 00:00:00 2001 From: nperez Date: Fri, 28 Apr 2023 16:52:45 +0200 Subject: [PATCH 5/7] cross.val check --- modules/Skill/Skill.R | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/modules/Skill/Skill.R b/modules/Skill/Skill.R index 122aa722..7b2ef7ed 100644 --- a/modules/Skill/Skill.R +++ b/modules/Skill/Skill.R @@ -92,6 +92,11 @@ compute_skill_metrics <- function(recipe, data) { } else { na.rm = recipe$Analysis$remove_NAs } + if (is.null(recipe$Analysis$Workflow$Skill$cross.val)) { + cross.val <- FALSE + } else { + cross.val <- recipe$Analysis$Workflow$Skill$cross.val + } skill_metrics <- list() for (metric in strsplit(metrics, ", | |,")[[1]]) { # Whether the fair version of the metric is to be computed @@ -113,12 +118,10 @@ compute_skill_metrics <- function(recipe, data) { time_dim = time_dim, memb_dim = memb_dim, Fair = Fair, - cross.val = recipe$Analysis$Workflow$Skill$cross.val, + cross.val = cross.val, ncores = ncores) - skill <- lapply(skill, function(x) { - .drop_dims(x)}) - skill_metrics[[ metric ]] <- skill$rpss - skill_metrics[[ paste0(metric, "_significance") ]] <- skill$sign + skill <- .drop_dims(skill) + skill_metrics[[ metric ]] <- skill rps_clim <- Apply(list(data$obs$data), target_dims = c(time_dim, memb_dim), RPS_clim)$output1 @@ -130,7 +133,7 @@ compute_skill_metrics <- function(recipe, data) { time_dim = time_dim, memb_dim = memb_dim, Fair = Fair, - cross.val = recipe$Analysis$Workflow$Skill$cross.val, + cross.val = cross.val, ncores = ncores) skill <- lapply(skill, function(x) { .drop_dims(x)}) @@ -142,7 +145,7 @@ compute_skill_metrics <- function(recipe, data) { time_dim = time_dim, memb_dim = memb_dim, prob_thresholds = 0.1, - cross.val = recipe$Analysis$Workflow$Skill$cross.val, + cross.val = cross.val, Fair = Fair, ncores = ncores) skill <- lapply(skill, function(x) { @@ -155,7 +158,7 @@ compute_skill_metrics <- function(recipe, data) { time_dim = time_dim, memb_dim = memb_dim, prob_thresholds = 0.9, - cross.val = recipe$Analysis$Workflow$Skill$cross.val, + cross.val = cross.val, Fair = Fair, ncores = ncores) skill <- lapply(skill, function(x) { @@ -186,9 +189,11 @@ compute_skill_metrics <- function(recipe, data) { .drop_dims(x)}) skill_metrics[[ metric ]] <- skill$crpss skill_metrics[[ paste0(metric, "_significance") ]] <- skill$sign - } else if (metric == 'RMS') { + } else if (metric == 'rms') { source("https://earth.bsc.es/gitlab/es/s2dv/-/raw/master/R/RMS.R") - skill <- RMS(exp = data$hcst$data, + hcst_mean <- Apply(list(data$hcst$data), target_dims = memb_dim, + fun = mean, na.rm = na.rm, ncores = ncores)$output1 + skill <- RMS(exp = hcst_mean, obs = data$obs$data, time_dim = time_dim, dat_dim = NULL, comp_dim = NULL, limits = NULL, conf = FALSE, alpha = 0.05, ncores = ncores) -- GitLab From 3a7a39d4bbce672ffe591b271efc62fed9929e57 Mon Sep 17 00:00:00 2001 From: nperez Date: Fri, 28 Apr 2023 16:58:31 +0200 Subject: [PATCH 6/7] ensemble mean --- modules/Skill/Skill.R | 1 + 1 file changed, 1 insertion(+) diff --git a/modules/Skill/Skill.R b/modules/Skill/Skill.R index 7b2ef7ed..20eec5c3 100644 --- a/modules/Skill/Skill.R +++ b/modules/Skill/Skill.R @@ -193,6 +193,7 @@ compute_skill_metrics <- function(recipe, data) { source("https://earth.bsc.es/gitlab/es/s2dv/-/raw/master/R/RMS.R") hcst_mean <- Apply(list(data$hcst$data), target_dims = memb_dim, fun = mean, na.rm = na.rm, ncores = ncores)$output1 + hcst_mean <- InsertDim(hcst_mean, pos = 1, len = 1, name = memb_dim) skill <- RMS(exp = hcst_mean, obs = data$obs$data, time_dim = time_dim, dat_dim = NULL, comp_dim = NULL, -- GitLab From 591234d3ec2882686a3a9844db51eba7fba8652a Mon Sep 17 00:00:00 2001 From: Victoria Agudetse Roures Date: Tue, 2 May 2023 13:20:41 +0200 Subject: [PATCH 7/7] Change name of cross.val parameter, fix _clim metric name, fix pipeline --- modules/Skill/Skill.R | 8 +++++--- tests/testthat/test-decadal_monthly_2.R | 2 +- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/modules/Skill/Skill.R b/modules/Skill/Skill.R index 20eec5c3..7eb44321 100644 --- a/modules/Skill/Skill.R +++ b/modules/Skill/Skill.R @@ -92,10 +92,12 @@ compute_skill_metrics <- function(recipe, data) { } else { na.rm = recipe$Analysis$remove_NAs } - if (is.null(recipe$Analysis$Workflow$Skill$cross.val)) { + if (is.null(recipe$Analysis$Workflow$Skill$cross_validation)) { + warn(recipe$Run$logger, + "cross_validation parameter not defined, setting it to FALSE.") cross.val <- FALSE } else { - cross.val <- recipe$Analysis$Workflow$Skill$cross.val + cross.val <- recipe$Analysis$Workflow$Skill$cross_validation } skill_metrics <- list() for (metric in strsplit(metrics, ", | |,")[[1]]) { @@ -126,7 +128,7 @@ compute_skill_metrics <- function(recipe, data) { target_dims = c(time_dim, memb_dim), RPS_clim)$output1 rps_clim <- .drop_dims(rps_clim) - skill_metrics[['rps_clim']] <- rps_clim + skill_metrics[[paste0(metric, "_clim")]] <- rps_clim # Ranked Probability Skill Score and Fair version } else if (metric %in% c('rpss', 'frpss')) { skill <- RPSS(data$hcst$data, data$obs$data, diff --git a/tests/testthat/test-decadal_monthly_2.R b/tests/testthat/test-decadal_monthly_2.R index a6ca9254..01fe0440 100644 --- a/tests/testthat/test-decadal_monthly_2.R +++ b/tests/testthat/test-decadal_monthly_2.R @@ -164,7 +164,7 @@ TRUE ) expect_equal( names(skill_metrics), -c("rpss_specs", "enscorr_specs", "frps_specs", "frpss_specs", "bss10_specs", "frps") +c("rpss_specs", "enscorr_specs", "frps_specs", "frpss_specs", "bss10_specs", "frps", "frps_clim") ) expect_equal( class(skill_metrics$rpss_specs), -- GitLab