From d9c6f8eb98ae2f4825e2caaf1ce06c83822b24ab Mon Sep 17 00:00:00 2001 From: aho Date: Fri, 2 Oct 2020 11:19:46 +0200 Subject: [PATCH 1/3] Bugfix for NA existence exclusion. --- R/Clim.R | 29 ++++++++++++++++++----------- tests/testthat/test-Clim.R | 20 ++++++-------------- 2 files changed, 24 insertions(+), 25 deletions(-) diff --git a/R/Clim.R b/R/Clim.R index d879fc4..46423ab 100644 --- a/R/Clim.R +++ b/R/Clim.R @@ -166,23 +166,30 @@ Clim <- function(exp, obs, time_dim = 'sdate', dat_dim = c('dataset', 'member'), #---------------------------------- # Remove all sdate if not complete along dat_dim - - pos <- rep(0, length(dat_dim)) - for (i in 1:length(dat_dim)) { #[dat, sdate] - ## dat_dim: [dataset, member] - pos[i] <- which(names(dim(obs)) == dat_dim[i]) + if (ftime_dim %in% names(dim(exp))) { + pos <- rep(0, length(time_dim) + length(ftime_dim)) + ## time_dim: [sdate]; ftime_dim: [ftime] + pos[1] <- which(names(dim(obs)) == time_dim) + pos[2] <- which(names(dim(obs)) == ftime_dim) + } else { + pos <- rep(0, length(time_dim)) + ## time_dim: [sdate] + pos[1] <- which(names(dim(obs)) == time_dim) } - outrows_exp <- MeanDims(exp, pos, na.rm = FALSE) + - MeanDims(obs, pos, na.rm = FALSE) + + avg_dims <- c(1:length(dim(obs)))[-pos] + outrows_exp <- MeanDims(exp, avg_dims, na.rm = FALSE) + + MeanDims(obs, avg_dims, na.rm = FALSE) outrows_obs <- outrows_exp + # outrows: [sdate] - for (i in 1:length(pos)) { - outrows_exp <- InsertDim(outrows_exp, pos[i], dim(exp)[pos[i]]) - outrows_obs <- InsertDim(outrows_obs, pos[i], dim(obs)[pos[i]]) + for (i in 1:length(avg_dims)) { + outrows_exp <- InsertDim(outrows_exp, avg_dims[i], dim(exp)[avg_dims[i]]) + outrows_obs <- InsertDim(outrows_obs, avg_dims[i], dim(obs)[avg_dims[i]]) } + exp[which(is.na(outrows_exp))] <- NA obs[which(is.na(outrows_obs))] <- NA - #----------------------------------- if (method == 'clim') { diff --git a/tests/testthat/test-Clim.R b/tests/testthat/test-Clim.R index 1d2ff3e..d5bfdd4 100644 --- a/tests/testthat/test-Clim.R +++ b/tests/testthat/test-Clim.R @@ -10,13 +10,9 @@ context("s2dv::Clim tests") ftime = 3, lon = 2, lat = 4, sdate = 5)) # dat2 exp2 <- exp1 - set.seed(1) - na <- floor(runif(5, min = 1, max = 360)) - exp2[na] <- NA + exp2[1, 2, 1, 1, 1, 1] <- NA obs2 <- obs1 - set.seed(2) - na <- floor(runif(30, min = 1, max = 120)) - obs2[na] <- NA + obs2[1, 1, 1, 1, 1, 2] <- NA # dat3 set.seed(1) @@ -175,21 +171,21 @@ test_that("3. Output checks: dat2", { expect_equal( (Clim(exp2, obs2)$clim_obs)[1:5], - c(0.23142987, -0.60462627, -0.03669491, -0.14193572, 0.61163024), + c(0.54451161, -0.60462627, 0.06609153, 0.18601029, 0.50553522), tolerance = 0.001 ) expect_equal( (Clim(exp2, obs2)$clim_exp)[1:5], - c(0.01054951, -0.04744191, -0.03533071, 0.14149945, 0.02359945), + c(-0.14639997, 0.01180200, 0.69685184, 0.14149945, 0.02359945), tolerance = 0.001 ) expect_equal( length(which(is.na(Clim(exp2, obs2, na.rm = FALSE)$clim_exp))), - 57 + 24 ) expect_equal( length(which(is.na(Clim(exp2, obs2, na.rm = FALSE)$clim_obs))), - 19 + 8 ) expect_equal( length(which(is.na(Clim(exp2, obs2)$clim_obs))), @@ -212,9 +208,5 @@ test_that("4. Output checks: dat3", { tolerance = 0.00001 ) - - - - }) -- GitLab From 95f357d41cc467bd665c92e966e988f0201b8d20 Mon Sep 17 00:00:00 2001 From: aho Date: Fri, 2 Oct 2020 11:21:21 +0200 Subject: [PATCH 2/3] Update NEWS.md for Clim() bugifx --- NEWS.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 9979a52..1557816 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,7 +8,8 @@ - Add p-value by ANOVA in Trend(). - Change MeanDims() na.rm default to FALSE to be in line with mean() - Remove unecessary parameter checks in Clim() - +- Bugfix for Clim() per-paired NA exclusion. + # s2dv 0.0.1 (Release date: 2020-02-07) - The package is the advanced version of package 's2dverification', adopting the regime of package 'multiApply' for all the analytic functions. Most of the other functions for plotting and data retrieval in 's2dverification' are also preserved in this package. - Because of the adoption of 'multiApply' regime, the functions work well with package 'startR'. -- GitLab From 90403e37b9555d8e2c705a6a3cb3917e00bd03c1 Mon Sep 17 00:00:00 2001 From: aho Date: Fri, 23 Oct 2020 16:29:06 +0200 Subject: [PATCH 3/3] Try to solve NA per-paired problem, but it seems wrong --- R/Clim.R | 71 +++++++++++++++++++++++++++++++++++++- tests/testthat/test-Clim.R | 43 +++++++++++------------ 2 files changed, 91 insertions(+), 23 deletions(-) diff --git a/R/Clim.R b/R/Clim.R index 46423ab..7d4ab88 100644 --- a/R/Clim.R +++ b/R/Clim.R @@ -160,12 +160,81 @@ Clim <- function(exp, obs, time_dim = 'sdate', dat_dim = c('dataset', 'member'), order_obs <- match(name_exp, name_obs) obs <- Reorder(obs, order_obs) - +#print('dim(exp) 1') +#print(dim(exp)) +#print('dim(obs) 1') +#print(dim(obs)) +#print('exp[1,,,1,1,1] 1') +#print(exp[1,,,1,1,1]) +#print('obs[1,1,,,1,1] 1') +#print(obs[1,1,,,1,1]) ############################### # Calculate Clim #---------------------------------- # Remove all sdate if not complete along dat_dim + + # First, for loop through dat_dim to create per-paired NA array + dat_dim_ind <- match(dat_dim, name_exp) + na_array <- array(0, dim = dim(exp)[-dat_dim_ind]) + + na_array <- na_array + MeanDims(exp, dat_dim_ind, na.rm = FALSE) + na_array <- na_array + MeanDims(obs, dat_dim_ind, na.rm = FALSE) +#print('na_array') +#print(dim(na_array)) +#print(summary(na_array)) + +# na_array <- as.logical(na_array) # TRUE is non-NA, NA is NA + + exp_dat_dim_length <- dim(exp)[dat_dim_ind] #[dat = 2, memb = 3, dat = 2] + exp_dat_dim_length <- lapply(lapply(exp_dat_dim_length, seq, 1), sort) # $dat: 1 2; $memb: 1 2 3 + exp_ind_list_for_Subset <- expand.grid(exp_dat_dim_length) +#print('exp_ind_list_for_Subset') +#print(exp_ind_list_for_Subset) + obs_dat_dim_length <- dim(obs)[dat_dim_ind] + obs_dat_dim_length <- lapply(lapply(obs_dat_dim_length, seq, 1), sort) + obs_ind_list_for_Subset <- expand.grid(obs_dat_dim_length) +#print('obs_ind_list_for_Subset') +#print(obs_ind_list_for_Subset) + + list_exp <- vector('list', length = length(exp_ind_list_for_Subset[, 1])) + list_obs <- vector('list', length = length(obs_ind_list_for_Subset[, 1])) + + # exp + for (kk in 1:length(list_exp)) { + tmp_exp <- ClimProjDiags::Subset(exp, dat_dim, as.list(exp_ind_list_for_Subset[kk, ]), drop = 'selected') + tmp_exp[which(is.na(na_array))] <- NA + list_exp[[kk]] <- tmp_exp + } + # obs + for (kk in 1:length(list_obs)) { + tmp_obs <- ClimProjDiags::Subset(obs, dat_dim, as.list(obs_ind_list_for_Subset[kk, ]), drop = 'selected') + tmp_obs[which(is.na(na_array))] <- NA + list_obs[[kk]] <- tmp_obs + } + +#print('str(list_exp)') +#print(str(list_exp)) +#print('str(list_obs)') +#print(str(list_obs)) + + exp <- array(unlist(list_exp), dim = c(dim(exp)[-dat_dim_ind], dim(exp)[dat_dim_ind])) + obs <- array(unlist(list_obs), dim = c(dim(obs)[-dat_dim_ind], dim(obs)[dat_dim_ind])) + + # Reorder exp and obs back to the original one + ori_order <- match(name_exp, names(dim(exp))) + exp <- Reorder(exp, ori_order) + obs <- Reorder(obs, ori_order) +#print('dim(exp) 2') +#print(dim(exp)) +#print('dim(obs) 2') +#print(dim(obs)) +#print('exp[1,,,1,1,1] 2') +#print(exp[1,,,1,1,1]) +#print('obs[1,1,,,1,1] 2') +#print(obs[1,1,,,1,1]) +#-------------------------------------------- + if (ftime_dim %in% names(dim(exp))) { pos <- rep(0, length(time_dim) + length(ftime_dim)) ## time_dim: [sdate]; ftime_dim: [ftime] diff --git a/tests/testthat/test-Clim.R b/tests/testthat/test-Clim.R index d5bfdd4..d96a42a 100644 --- a/tests/testthat/test-Clim.R +++ b/tests/testthat/test-Clim.R @@ -3,14 +3,14 @@ context("s2dv::Clim tests") ############################################## # dat1 set.seed(1) - exp1 <- array(rnorm(360), dim = c(dataset = 1, member = 3, sdate = 5, - ftime = 3, lon = 2, lat = 4)) + exp1 <- array(rnorm(360), dim = c(dataset = 2, member = 3, sdate = 5, + ftime = 3, lon = 2, lat = 2)) set.seed(2) - obs1 <- array(rnorm(120), dim = c(dataset = 1, member = 1, - ftime = 3, lon = 2, lat = 4, sdate = 5)) + obs1 <- array(rnorm(60), dim = c(dataset = 1, member = 1, + ftime = 3, lon = 2, lat = 2, sdate = 5)) # dat2 exp2 <- exp1 - exp2[1, 2, 1, 1, 1, 1] <- NA + exp2[1, 2, , 1, 1, 1] <- NA obs2 <- obs1 obs2[1, 1, 1, 1, 1, 2] <- NA @@ -108,59 +108,59 @@ test_that("2. Output checks: dat1", { expect_equal( dim(Clim(exp1, obs1)$clim_exp), - c(dataset = 1, member = 3, ftime = 3, lon = 2, lat = 4) + c(dataset = 2, member = 3, ftime = 3, lon = 2, lat = 2) ) expect_equal( dim(Clim(exp1, obs1, memb = FALSE)$clim_exp), - c(dataset = 1, ftime = 3, lon = 2, lat = 4) + c(dataset = 2, ftime = 3, lon = 2, lat = 2) ) expect_equal( dim(Clim(exp1, obs1, time_dim = 'lon')$clim_exp), - c(dataset = 1, member = 3, sdate = 5, ftime = 3, lat = 4) + c(dataset = 2, member = 3, sdate = 5, ftime = 3, lat = 2) ) expect_equal( dim(Clim(exp1, obs1, method = 'kharin')$clim_exp), - c(sdate = 5, dataset = 1, member = 3, ftime = 3, lon = 2, lat = 4) + c(sdate = 5, dataset = 2, member = 3, ftime = 3, lon = 2, lat = 2) ) expect_equal( dim(Clim(exp1, obs1, method = 'NDV')$clim_exp), - c(sdate = 5, dataset = 1, member = 3, ftime = 3, lon = 2, lat = 4) + c(sdate = 5, dataset = 2, member = 3, ftime = 3, lon = 2, lat = 2) ) expect_equal( dim(Clim(exp1, obs1)$clim_obs), - c(dataset = 1, member = 1, ftime = 3, lon = 2, lat = 4) + c(dataset = 1, member = 1, ftime = 3, lon = 2, lat = 2) ) expect_equal( dim(Clim(exp1, obs1, method = 'kharin')$clim_obs), - c(dataset = 1, member = 1, ftime = 3, lon = 2, lat = 4) + c(dataset = 1, member = 1, ftime = 3, lon = 2, lat = 2) ) expect_equal( dim(Clim(exp1, obs1, method = 'NDV')$clim_obs), - c(dataset = 1, member = 1, ftime = 3, lon = 2, lat = 4) + c(dataset = 1, member = 1, ftime = 3, lon = 2, lat = 2) ) expect_equal( (Clim(exp1, obs1)$clim_obs)[1:5], - c(0.14831161, -0.60462627, 0.06609153, -0.23378059, 0.50553522), + c(-0.6389627, -1.0815757, 0.4899925, -0.4436428, 0.1289446), tolerance = 0.001 ) expect_equal( (Clim(exp1, obs1, memb = FALSE)$clim_exp)[1:5], - c(0.10084284, 0.06407350, 0.09028584, 0.17526332, 0.18004387), + c(0.24870405, -0.08378771, 0.23637962, 0.02916955, 0.18643424), tolerance = 0.001 ) expect_equal( max(Clim(exp1, obs1)$clim_exp, na.rm = T), - 1.026186, + 1.065416, tolerance = 0.001 ) expect_equal( max(Clim(exp1, obs1, method = 'kharin')$clim_exp, na.rm = T), - 2.282634, + 2.634209, tolerance = 0.001 ) expect_equal( min(Clim(exp1, obs1, method = 'NDV')$clim_exp, na.rm = T), - -4.025745, + -3.156395, tolerance = 0.001 ) @@ -171,13 +171,12 @@ test_that("3. Output checks: dat2", { expect_equal( (Clim(exp2, obs2)$clim_obs)[1:5], - c(0.54451161, -0.60462627, 0.06609153, 0.18601029, 0.50553522), + c(NaN, -1.0815757, 0.4899925, NaN, 0.1289446), tolerance = 0.001 ) expect_equal( - (Clim(exp2, obs2)$clim_exp)[1:5], - c(-0.14639997, 0.01180200, 0.69685184, 0.14149945, 0.02359945), - tolerance = 0.001 + all(is.na(Clim(exp2, obs2)$clim_exp[, , 1, 2, 2])), + TRUE ) expect_equal( length(which(is.na(Clim(exp2, obs2, na.rm = FALSE)$clim_exp))), -- GitLab