From 327d6d55401e5b3f65a68947709462f59cd2807c Mon Sep 17 00:00:00 2001 From: aho Date: Thu, 1 Apr 2021 15:14:36 +0200 Subject: [PATCH 1/8] Add Spectrum from s2dverification --- R/Spectrum.R | 72 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 72 insertions(+) create mode 100644 R/Spectrum.R diff --git a/R/Spectrum.R b/R/Spectrum.R new file mode 100644 index 0000000..e99d323 --- /dev/null +++ b/R/Spectrum.R @@ -0,0 +1,72 @@ +#'Estimates Frequency Spectrum +#' +#'This function estimates the frequency spectrum of the data array together +#'with its 95\% and 99\% significance level. The output is provided as an +#'array with dimensions c(number of frequencies, 4). The column contains the +#'frequency values, the power, the 95\% significance level and the 99\% one.\cr +#'The spectrum estimation relies on a R built-in function and the significance +#'levels are estimated by a Monte-Carlo method. +#' +#'@param data A vector or numeric array of which the frequency spectrum is +#' required. If it's a vector, it should be a time series. If it's an array, +#' the dimensions must have at least 'time_dim'. The data is assumed to be +#' evenly spaced in time. +#' +#'@return Frequency spectrum with dimensions c(number of frequencies, 4). The +#' column contains the frequency values, the power, the 95\% significance +#' level and the 99\% one. +#' +#'@examples +#'# Load sample data as in Load() example: +#'example(Load) +#' +#'ensmod <- Mean1Dim(sampleData$mod, 2) +#'for (jstartdate in 1:3) { +#' spectrum <- Spectrum(ensmod[1, jstartdate, ]) +#' for (jlen in 1:dim(spectrum)[1]) { +#' if (spectrum[jlen, 2] > spectrum[jlen, 4]) { +#' ensmod[1, jstartdate, ] <- Filter(ensmod[1, jstartdate, ], +#' spectrum[jlen, 1]) +#' } +#' } +#'} +#' \donttest{ +#'PlotAno(InsertDim(ensmod, 2, 1), sdates = startDates, fileout = +#' 'filtered_ensemble_mean.eps') +#' } +#' +#'@importFrom stats spectrum cor rnorm sd quantile +#'@export +Spectrum <- function(data) { + + data <- data[is.na(data) == FALSE] + ndat <- length(data) + + if (ndat >= 3) { + tmp <- spectrum(data, plot = FALSE) + output <- array(dim = c(length(tmp$spec), 4)) + output[, 1] <- tmp$freq + output[, 2] <- tmp$spec + ntir <- 100 + store <- array(dim = c(ntir, length(tmp$spec))) + for (jt in 1:ntir) { + toto <- mean(data) + alpha1 <- cor(data[2:ndat], data[1:(ndat - 1)]) + for (ind in 2:ndat) { + b <- rnorm(1, mean(data) * (1 - alpha1), sd(data) * sqrt(1 - + alpha1 ^ 2)) + toto <- c(toto, toto[ind - 1] * alpha1 + b) + } + toto2 <- spectrum(toto, plot = FALSE) + store[jt, ] <- toto2$spec + } + for (jx in 1:length(tmp$spec)) { + output[jx, 3] <- quantile(store[, jx], 0.95) + output[jx, 4] <- quantile(store[, jx], 0.99) + } + } else { + output <- NA + } + + output +} -- GitLab From aebdc0a50949c8993da2d9e3e3ec2d2dad88756d Mon Sep 17 00:00:00 2001 From: aho Date: Thu, 8 Apr 2021 17:34:26 +0200 Subject: [PATCH 2/8] Add Spectrum() --- NAMESPACE | 2 + R/Spectrum.R | 93 ++++++++++++++++++++++------- man/Spectrum.Rd | 52 ++++++++++++++++ tests/testthat/test-Spectrum.R | 105 +++++++++++++++++++++++++++++++++ 4 files changed, 230 insertions(+), 22 deletions(-) create mode 100644 man/Spectrum.Rd create mode 100644 tests/testthat/test-Spectrum.R diff --git a/NAMESPACE b/NAMESPACE index b12cd82..341ae3a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -45,6 +45,7 @@ export(Reorder) export(SPOD) export(Season) export(Smoothing) +export(Spectrum) export(TPI) export(ToyModel) export(Trend) @@ -98,5 +99,6 @@ importFrom(stats,qt) importFrom(stats,quantile) importFrom(stats,rnorm) importFrom(stats,sd) +importFrom(stats,spectrum) importFrom(stats,ts) importFrom(stats,window) diff --git a/R/Spectrum.R b/R/Spectrum.R index e99d323..6890536 100644 --- a/R/Spectrum.R +++ b/R/Spectrum.R @@ -1,43 +1,92 @@ -#'Estimates Frequency Spectrum +#'Estimate frequency spectrum #' -#'This function estimates the frequency spectrum of the data array together -#'with its 95\% and 99\% significance level. The output is provided as an -#'array with dimensions c(number of frequencies, 4). The column contains the -#'frequency values, the power, the 95\% significance level and the 99\% one.\cr -#'The spectrum estimation relies on a R built-in function and the significance -#'levels are estimated by a Monte-Carlo method. +#'Estimate the frequency spectrum of the data array together with its 95\% and +#'99\% significance level. The output is provided as an array with dimensions +#"c(number of frequencies, 4). The column contains the frequency values, the +#'power, the 95\% significance level and the 99\% one.\cr +#'The spectrum estimation relies on a R built-in function \code{spectrum()} +#'and the significance levels are estimated by a Monte-Carlo method. #' #'@param data A vector or numeric array of which the frequency spectrum is #' required. If it's a vector, it should be a time series. If it's an array, #' the dimensions must have at least 'time_dim'. The data is assumed to be #' evenly spaced in time. +#'@param time_dim A character string indicating the dimension along which to +#' compute the frequency spectrum. The default value is 'ftime'. +#'@param ncores An integer indicating the number of cores to use for parallel +#' computation. The default value is NULL. #' -#'@return Frequency spectrum with dimensions c(number of frequencies, 4). The -#' column contains the frequency values, the power, the 95\% significance -#' level and the 99\% one. +#'@return A numeric array of the frequency spectrum with dimensions +#' c( = number of frequencies, stats = 4, the rest of the +#' dimensions of 'data'). The 'stats' dimension contains the frequency values, +#' the power, the 95\% significance level and the 99\% one. #' #'@examples #'# Load sample data as in Load() example: #'example(Load) #' -#'ensmod <- Mean1Dim(sampleData$mod, 2) -#'for (jstartdate in 1:3) { -#' spectrum <- Spectrum(ensmod[1, jstartdate, ]) -#' for (jlen in 1:dim(spectrum)[1]) { -#' if (spectrum[jlen, 2] > spectrum[jlen, 4]) { -#' ensmod[1, jstartdate, ] <- Filter(ensmod[1, jstartdate, ], -#' spectrum[jlen, 1]) +#'ensmod <- MeanDims(sampleData$mod, c(1, 2)) +#'spectrum <- Spectrum(ensmod) +#'for (jsdate in 1:dim(spectrum)['sdate']) { +#' for (jlen in 1:dim(spectrum)['ftime']) { +#' if (spectrum[jlen, 2, jsdate] > spectrum[jlen, 4, jsdate]) { +#' ensmod[jsdate, ] <- Filter(ensmod[jsdate, ], +#' spectrum[jlen, 1, jsdate]) #' } #' } #'} #' \donttest{ -#'PlotAno(InsertDim(ensmod, 2, 1), sdates = startDates, fileout = -#' 'filtered_ensemble_mean.eps') +#'PlotAno(InsertDim(ensmod, 2, 1), sdates = startDates) #' } #' #'@importFrom stats spectrum cor rnorm sd quantile #'@export -Spectrum <- function(data) { +Spectrum <- function(data, time_dim = 'ftime', ncores = NULL) { + + # Check inputs + ## data + if (is.null(data)) { + stop("Parameter 'data' cannot be NULL.") + } + if (!is.numeric(data)) { + stop("Parameter 'data' must be a numeric array.") + } + if (is.null(dim(data))) { #is vector + dim(data) <- c(length(data)) + names(dim(data)) <- time_dim + } + if(any(is.null(names(dim(data))))| any(nchar(names(dim(data))) == 0)) { + stop("Parameter 'data' must have dimension names.") + } + ## time_dim + if (!is.character(time_dim) | length(time_dim) > 1) { + stop("Parameter 'time_dim' must be a character string.") + } + if (!time_dim %in% names(dim(data))) { + stop("Parameter 'time_dim' is not found in 'data' dimension.") + } + ## ncores + if (!is.null(ncores)) { + if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 | + length(ncores) > 1) { + stop("Parameter 'ncores' must be a positive integer.") + } + } + + ############################### + # Calculate Spectrum + + output <- Apply(list(data), + target_dims = time_dim, + fun = .Spectrum, + output_dims = c(time_dim, 'stats'), + ncores = ncores)$output1 + + return(output) +} + +.Spectrum <- function(data) { + # data: [time] data <- data[is.na(data) == FALSE] ndat <- length(data) @@ -67,6 +116,6 @@ Spectrum <- function(data) { } else { output <- NA } - - output + + return(invisible(output)) } diff --git a/man/Spectrum.Rd b/man/Spectrum.Rd new file mode 100644 index 0000000..7f23d40 --- /dev/null +++ b/man/Spectrum.Rd @@ -0,0 +1,52 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Spectrum.R +\name{Spectrum} +\alias{Spectrum} +\title{Estimate frequency spectrum} +\usage{ +Spectrum(data, time_dim = "ftime", ncores = NULL) +} +\arguments{ +\item{data}{A vector or numeric array of which the frequency spectrum is +required. If it's a vector, it should be a time series. If it's an array, +the dimensions must have at least 'time_dim'. The data is assumed to be +evenly spaced in time.} + +\item{time_dim}{A character string indicating the dimension along which to +compute the frequency spectrum. The default value is 'ftime'.} + +\item{ncores}{An integer indicating the number of cores to use for parallel +computation. The default value is NULL.} +} +\value{ +A numeric array of the frequency spectrum with dimensions + c( = number of frequencies, stats = 4, the rest of the + dimensions of 'data'). The 'stats' dimension contains the frequency values, + the power, the 95\% significance level and the 99\% one. +} +\description{ +Estimate the frequency spectrum of the data array together with its 95\% and +99\% significance level. The output is provided as an array with dimensions +power, the 95\% significance level and the 99\% one.\cr +The spectrum estimation relies on a R built-in function \code{spectrum()} +and the significance levels are estimated by a Monte-Carlo method. +} +\examples{ +# Load sample data as in Load() example: +example(Load) + +ensmod <- MeanDims(sampleData$mod, c(1, 2)) +spectrum <- Spectrum(ensmod) +for (jsdate in 1:dim(spectrum)['sdate']) { + for (jlen in 1:dim(spectrum)['ftime']) { + if (spectrum[jlen, 2, jsdate] > spectrum[jlen, 4, jsdate]) { + ensmod[jsdate, ] <- Filter(ensmod[jsdate, ], + spectrum[jlen, 1, jsdate]) + } + } +} + \donttest{ +PlotAno(InsertDim(ensmod, 2, 1), sdates = startDates) + } + +} diff --git a/tests/testthat/test-Spectrum.R b/tests/testthat/test-Spectrum.R new file mode 100644 index 0000000..0a54101 --- /dev/null +++ b/tests/testthat/test-Spectrum.R @@ -0,0 +1,105 @@ +context("s2dv::Spectrum tests") + +############################################## + # dat1 + set.seed(1) + dat1 <- array(rnorm(10), dim = c(dat = 1, ftime = 5, sdate = 2)) + + # dat2 + set.seed(10) + dat2 <- c(1:10) + rnorm(10) + + # dat3 + dat3 <- dat2 + dat3[2] <- NA +############################################## + +test_that("1. Input checks", { + + # data + expect_error( + Spectrum(c()), + "Parameter 'data' cannot be NULL." + ) + expect_error( + Spectrum(c(NA, NA)), + "Parameter 'data' must be a numeric array." + ) + expect_error( + Spectrum(array(1:10, dim = c(2, 5))), + "Parameter 'data' must have dimension names." + ) + # time_dim + expect_error( + Spectrum(array(c(1:25), dim = c(dat = 1, date = 5, sdate = 5))), + "Parameter 'time_dim' is not found in 'data' dimension." + ) + expect_error( + Spectrum(dat1, time_dim = 2), + "Parameter 'time_dim' must be a character string." + ) + # ncores + expect_error( + Spectrum(dat1, ncore = 3.5), + "Parameter 'ncores' must be a positive integer." + ) +}) + +############################################## +test_that("2. Output checks: dat1", { + + expect_equal( + dim(Spectrum(dat1)), + c(ftime = 2, stats = 4, dat = 1, sdate = 2) + ) + expect_equal( + Spectrum(dat1)[, 1, 1, 2], + c(0.2, 0.4), + tolerance = 0.0001 + ) + expect_equal( + Spectrum(dat1)[, 2, 1, 2], + c(0.89583007, 0.05516983), + tolerance = 0.0001 + ) + +}) + +############################################## +test_that("3. Output checks: dat2", { + + expect_equal( + dim(Spectrum(dat2)), + c(ftime = 5, stats = 4) + ) + expect_equal( + Spectrum(dat2)[, 1], + c(0.1, 0.2, 0.3, 0.4, 0.5), + tolerance = 0.0001 + ) + expect_equal( + Spectrum(dat2)[, 2], + c(0.1767994, 1.0113808, 0.3341372, 0.1807377, 1.0594528), + tolerance = 0.0001 + ) + +}) + +############################################## +test_that("4. Output checks: dat3", { + expect_equal( + dim(Spectrum(dat3)), + c(ftime = 4, stats = 4) + ) + expect_equal( + Spectrum(dat3)[, 1], + c(0.1111111, 0.2222222, 0.3333333, 0.4444444), + tolerance = 0.0001 + ) + expect_equal( + Spectrum(dat3)[, 2], + c(0.7204816, 0.6529411, 0.2605188, 0.7009824), + tolerance = 0.0001 + ) + +}) -- GitLab From 8c5bb47bf6326b745e2b504c5dcb7ab3d47b64de Mon Sep 17 00:00:00 2001 From: aho Date: Thu, 8 Apr 2021 17:34:35 +0200 Subject: [PATCH 3/8] Update doc --- man/AMV.Rd | 29 ++++++++++++++--------------- man/GMST.Rd | 51 ++++++++++++++++++++++++--------------------------- man/GSAT.Rd | 30 +++++++++++++++--------------- man/SPOD.Rd | 30 +++++++++++++++--------------- man/TPI.Rd | 30 +++++++++++++++--------------- 5 files changed, 83 insertions(+), 87 deletions(-) diff --git a/man/AMV.Rd b/man/AMV.Rd index 0f81652..3cc4113 100644 --- a/man/AMV.Rd +++ b/man/AMV.Rd @@ -18,16 +18,15 @@ AMV( indices_for_clim = NULL, year_dim = "year", month_dim = "month", - member_dim = "member", + na.rm = TRUE, ncores = NULL ) } \arguments{ -\item{data}{A numerical array to be used for the index computation with the -dimensions: 1) latitude, longitude, start date, forecast month, and member -(in case of decadal predictions), 2) latitude, longitude, year, month and -member (in case of historical simulations), or 3) latitude, longitude, year -and month (in case of observations or reanalyses). This data has to be +\item{data}{A numerical array to be used for the index computation with, at least, the +dimensions: 1) latitude, longitude, start date and forecast month +(in case of decadal predictions), 2) latitude, longitude, year and month +(in case of historical simulations or observations). This data has to be provided, at least, over the whole region needed to compute the index.} \item{data_lats}{A numeric vector indicating the latitudes of the data.} @@ -70,7 +69,7 @@ climatology is calculated over the whole period. If the data are already anomalies, set it to FALSE. The default value is NULL.\cr In case of parameter 'type' is 'dcpp', 'indices_for_clim' must be relative to the first forecast year, and the climatology is automatically computed -over the actual common period for the different forecast years.} +over the common calendar period for the different forecast years.} \item{year_dim}{A character string indicating the name of the year dimension The default value is 'year'. Only used if parameter 'type' is 'hist' or @@ -80,18 +79,17 @@ The default value is 'year'. Only used if parameter 'type' is 'hist' or dimension. The default value is 'month'. Only used if parameter 'type' is 'hist' or 'obs'.} -\item{member_dim}{A character string indicating the name of the member -dimension. The default value is 'member'. Only used if parameter 'type' is -'dcpp' or 'hist'.} +\item{na.rm}{A logical value indicanting whether to remove NA values. The default +value is TRUE.} \item{ncores}{An integer indicating the number of cores to use for parallel computation. The default value is NULL.} } \value{ -A numerical array of the AMV index with the dimensions of: - 1) sdate, forecast year, and member (in case of decadal predictions); - 2) year and member (in case of historical simulations); or - 3) year (in case of observations or reanalyses). +A numerical array with the AMV index with the same dimensions as data except + the lat_dim, lon_dim and fmonth_dim (month_dim) in case of decadal predictions + (historical simulations or observations). In case of decadal predictions, a new dimension + 'fyear' is added. } \description{ The Atlantic Multidecadal Variability (AMV), also known as Atlantic @@ -100,7 +98,8 @@ surface temperatures (SST) over the North Atlantic Ocean on multi-decadal time scales. The AMV index is computed as the difference of weighted-averaged SST anomalies over the North Atlantic region (0ºN-60ºN, 280ºE-360ºE) and the weighted-averaged SST anomalies over 60ºS-60ºN, 0ºE-360ºE (Trenberth & -Dennis, 2005; Doblas-Reyes et al., 2013). +Dennis, 2005; Doblas-Reyes et al., 2013). If different members and/or datasets are provided, +the climatology (used to calculate the anomalies) is computed individually for all of them. } \examples{ ## Observations or reanalyses diff --git a/man/GMST.Rd b/man/GMST.Rd index 8021943..f23e60d 100644 --- a/man/GMST.Rd +++ b/man/GMST.Rd @@ -21,28 +21,25 @@ GMST( indices_for_clim = NULL, year_dim = "year", month_dim = "month", - member_dim = "member", + na.rm = TRUE, ncores = NULL ) } \arguments{ -\item{data_tas}{A numerical array indicating the surface air temperature data -to be used for the index computation with the dimensions: 1) latitude, -longitude, start date, forecast month, and member (in case of decadal -predictions), 2) latitude, longitude, year, month and member (in case of -historical simulations), or 3) latitude, longitude, year and month (in case -of observations or reanalyses). This data has to be provided, at least, -over the whole region needed to compute the index. The dimensions must be -identical to those of data_tos.} - -\item{data_tos}{A numerical array indicating the sea surface temperature data -to be used for the index computation with the dimensions: 1) latitude, -longitude, start date, forecast month, and member (in case of decadal -predictions), 2) latitude, longitude, year, month and member (in case of -historical simulations), or 3) latitude, longitude, year and month (in case -of observations or reanalyses). This data has to be provided, at least, -over the whole region needed to compute the index. The dimensions must be -identical to those of data_tas.} +\item{data_tas}{A numerical array with the surface air temperature data +to be used for the index computation with, at least, the +dimensions: 1) latitude, longitude, start date and forecast month +(in case of decadal predictions), 2) latitude, longitude, year and month +(in case of historical simulations or observations). This data has to be +provided, at least, over the whole region needed to compute the index. +The dimensions must be identical to thos of data_tos. +#'@param data_tos A numerical array with the sea surface temperature data +to be used for the index computation with, at least, the +dimensions: 1) latitude, longitude, start date and forecast month +(in case of decadal predictions), 2) latitude, longitude, year and month +(in case of historical simulations or observations). This data has to be +provided, at least, over the whole region needed to compute the index. +The dimensions must be identical to thos of data_tas.} \item{data_lats}{A numeric vector indicating the latitudes of the data.} @@ -90,7 +87,7 @@ climatology is calculated over the whole period. If the data are already anomalies, set it to FALSE. The default value is NULL.\cr In case of parameter 'type' is 'dcpp', 'indices_for_clim' must be relative to the first forecast year, and the climatology is automatically computed -over the actual common period for the different forecast years.} +over the common calendar period for the different forecast years.} \item{year_dim}{A character string indicating the name of the year dimension The default value is 'year'. Only used if parameter 'type' is 'hist' or @@ -100,23 +97,23 @@ The default value is 'year'. Only used if parameter 'type' is 'hist' or dimension. The default value is 'month'. Only used if parameter 'type' is 'hist' or 'obs'.} -\item{member_dim}{A character string indicating the name of the member -dimension. The default value is 'member'. Only used if parameter 'type' is -'dcpp' or 'hist'.} +\item{na.rm}{A logical value indicanting whether to remove NA values. The default +value is TRUE.} \item{ncores}{An integer indicating the number of cores to use for parallel computation. The default value is NULL.} } \value{ -A numerical array of the GMST anomalies with the dimensions of: - 1) sdate, forecast year, and member (in case of decadal predictions); - 2) year and member (in case of historical simulations); or - 3) year (in case of observations or reanalyses). +A numerical array with the GMST anomalies with the same dimensions as data_tas except + the lat_dim, lon_dim and fmonth_dim (month_dim) in case of decadal predictions + (historical simulations or observations). In case of decadal predictions, a new dimension + 'fyear' is added. } \description{ The Global Mean Surface Temperature (GMST) anomalies are computed as the weighted-averaged surface air temperature anomalies over land and sea surface -temperature anomalies over the ocean. +temperature anomalies over the ocean. If different members and/or datasets are provided, +the climatology (used to calculate the anomalies) is computed individually for all of them. } \examples{ ## Observations or reanalyses diff --git a/man/GSAT.Rd b/man/GSAT.Rd index d7fe3cf..9f03ff2 100644 --- a/man/GSAT.Rd +++ b/man/GSAT.Rd @@ -18,16 +18,15 @@ GSAT( indices_for_clim = NULL, year_dim = "year", month_dim = "month", - member_dim = "member", + na.rm = TRUE, ncores = NULL ) } \arguments{ -\item{data}{A numerical array to be used for the index computation with the -dimensions: 1) latitude, longitude, start date, forecast month, and member -(in case of decadal predictions), 2) latitude, longitude, year, month and -member (in case of historical simulations), or 3) latitude, longitude, year -and month (in case of observations or reanalyses). This data has to be +\item{data}{A numerical array to be used for the index computation with, at least, the +dimensions: 1) latitude, longitude, start date and forecast month +(in case of decadal predictions), 2) latitude, longitude, year and month +(in case of historical simulations or observations). This data has to be provided, at least, over the whole region needed to compute the index.} \item{data_lats}{A numeric vector indicating the latitudes of the data.} @@ -70,7 +69,7 @@ climatology is calculated over the whole period. If the data are already anomalies, set it to FALSE. The default value is NULL.\cr In case of parameter 'type' is 'dcpp', 'indices_for_clim' must be relative to the first forecast year, and the climatology is automatically computed -over the actual common period for the different forecast years.} +over the common calendar period for the different forecast years.} \item{year_dim}{A character string indicating the name of the year dimension The default value is 'year'. Only used if parameter 'type' is 'hist' or @@ -80,22 +79,23 @@ The default value is 'year'. Only used if parameter 'type' is 'hist' or dimension. The default value is 'month'. Only used if parameter 'type' is 'hist' or 'obs'.} -\item{member_dim}{A character string indicating the name of the member -dimension. The default value is 'member'. Only used if parameter 'type' is -'dcpp' or 'hist'.} +\item{na.rm}{A logical value indicanting whether to remove NA values. The default +value is TRUE.} \item{ncores}{An integer indicating the number of cores to use for parallel computation. The default value is NULL.} } \value{ -A numerical array of the GSAT anomalies with the dimensions of: - 1) sdate, forecast year, and member (in case of decadal predictions); - 2) year and member (in case of historical simulations); or - 3) year (in case of observations or reanalyses). +A numerical array with the GSAT anomalies with the same dimensions as data except + the lat_dim, lon_dim and fmonth_dim (month_dim) in case of decadal predictions + (historical simulations or observations). In case of decadal predictions, a new dimension + 'fyear' is added. } \description{ The Global Surface Air Temperature (GSAT) anomalies are computed as the -weighted-averaged surface air temperature anomalies over the global region. +weighted-averaged surface air temperature anomalies over the global region. +If different members and/or datasets are provided, the climatology (used to +calculate the anomalies) is computed individually for all of them. } \examples{ ## Observations or reanalyses diff --git a/man/SPOD.Rd b/man/SPOD.Rd index 0491739..4c4ed29 100644 --- a/man/SPOD.Rd +++ b/man/SPOD.Rd @@ -18,16 +18,15 @@ SPOD( indices_for_clim = NULL, year_dim = "year", month_dim = "month", - member_dim = "member", + na.rm = TRUE, ncores = NULL ) } \arguments{ -\item{data}{A numerical array to be used for the index computation with the -dimensions: 1) latitude, longitude, start date, forecast month, and member -(in case of decadal predictions), 2) latitude, longitude, year, month and -member (in case of historical simulations), or 3) latitude, longitude, year -and month (in case of observations or reanalyses). This data has to be +\item{data}{A numerical array to be used for the index computation with, at least, the +dimensions: 1) latitude, longitude, start date and forecast month +(in case of decadal predictions), 2) latitude, longitude, year and month +(in case of historical simulations or observations). This data has to be provided, at least, over the whole region needed to compute the index.} \item{data_lats}{A numeric vector indicating the latitudes of the data.} @@ -70,7 +69,7 @@ climatology is calculated over the whole period. If the data are already anomalies, set it to FALSE. The default value is NULL.\cr In case of parameter 'type' is 'dcpp', 'indices_for_clim' must be relative to the first forecast year, and the climatology is automatically computed -over the actual common period for the different forecast years.} +over the common calendar period for the different forecast years.} \item{year_dim}{A character string indicating the name of the year dimension The default value is 'year'. Only used if parameter 'type' is 'hist' or @@ -80,25 +79,26 @@ The default value is 'year'. Only used if parameter 'type' is 'hist' or dimension. The default value is 'month'. Only used if parameter 'type' is 'hist' or 'obs'.} -\item{member_dim}{A character string indicating the name of the member -dimension. The default value is 'member'. Only used if parameter 'type' is -'dcpp' or 'hist'.} +\item{na.rm}{A logical value indicanting whether to remove NA values. The default +value is TRUE.} \item{ncores}{An integer indicating the number of cores to use for parallel computation. The default value is NULL.} } \value{ -A numerical array of the SPOD index with the dimensions of: - 1) sdate, forecast year, and member (in case of decadal predictions); - 2) year and member (in case of historical simulations); or - 3) year (in case of observations or reanalyses). +A numerical array with the SPOD index with the same dimensions as data except + the lat_dim, lon_dim and fmonth_dim (month_dim) in case of decadal predictions + (historical simulations or observations). In case of decadal predictions, a new dimension + 'fyear' is added. } \description{ The South Pacific Ocean Dipole (SPOD) index is related to the El Nino-Southern Oscillation (ENSO) and the Inderdecadal Pacific Oscillation (IPO). The SPOD index is computed as the difference of weighted-averaged SST anomalies over 20ºS-48ºS, 165ºE-190ºE (NW pole) and the weighted-averaged SST -anomalies over 44ºS-65ºS, 220ºE-260ºE (SE pole) (Saurral et al., 2020). +anomalies over 44ºS-65ºS, 220ºE-260ºE (SE pole) (Saurral et al., 2020). +If different members and/or datasets are provided, the climatology (used to +calculate the anomalies) is computed individually for all of them. } \examples{ ## Observations or reanalyses diff --git a/man/TPI.Rd b/man/TPI.Rd index 3bdc17c..5e8f716 100644 --- a/man/TPI.Rd +++ b/man/TPI.Rd @@ -18,16 +18,15 @@ TPI( indices_for_clim = NULL, year_dim = "year", month_dim = "month", - member_dim = "member", + na.rm = TRUE, ncores = NULL ) } \arguments{ -\item{data}{A numerical array to be used for the index computation with the -dimensions: 1) latitude, longitude, start date, forecast month, and member -(in case of decadal predictions), 2) latitude, longitude, year, month and -member (in case of historical simulations), or 3) latitude, longitude, year -and month (in case of observations or reanalyses). This data has to be +\item{data}{A numerical array to be used for the index computation with, at least, the +dimensions: 1) latitude, longitude, start date and forecast month +(in case of decadal predictions), 2) latitude, longitude, year and month +(in case of historical simulations or observations). This data has to be provided, at least, over the whole region needed to compute the index.} \item{data_lats}{A numeric vector indicating the latitudes of the data.} @@ -70,7 +69,7 @@ climatology is calculated over the whole period. If the data are already anomalies, set it to FALSE. The default value is NULL.\cr In case of parameter 'type' is 'dcpp', 'indices_for_clim' must be relative to the first forecast year, and the climatology is automatically computed -over the actual common period for the different forecast years.} +over the common calendar period for the different forecast years.} \item{year_dim}{A character string indicating the name of the year dimension The default value is 'year'. Only used if parameter 'type' is 'hist' or @@ -80,24 +79,25 @@ The default value is 'year'. Only used if parameter 'type' is 'hist' or dimension. The default value is 'month'. Only used if parameter 'type' is 'hist' or 'obs'.} -\item{member_dim}{A character string indicating the name of the member -dimension. The default value is 'member'. Only used if parameter 'type' is -'dcpp' or 'hist'.} +\item{na.rm}{A logical value indicanting whether to remove NA values. The default +value is TRUE.} \item{ncores}{An integer indicating the number of cores to use for parallel computation. The default value is NULL.} } \value{ -A numerical array of the TPI index with the dimensions of: - 1) sdate, forecast year, and member (in case of decadal predictions); - 2) year and member (in case of historical simulations); or - 3) year (in case of observations or reanalyses). +A numerical array with the TPI index with the same dimensions as data except + the lat_dim, lon_dim and fmonth_dim (month_dim) in case of decadal predictions + (historical simulations or observations). In case of decadal predictions, a new dimension + 'fyear' is added. } \description{ The Tripole Index (TPI) for the Interdecadal Pacific Oscillation (IPO) is computed as the difference of weighted-averaged SST anomalies over 10ºS-10ºN, 170ºE-270ºE minus the mean of the weighted-averaged SST anomalies over -25ºN-45ºN, 140ºE-215ºE and 50ºS-15ºS, 150ºE-200ºE (Henley et al., 2015). +25ºN-45ºN, 140ºE-215ºE and 50ºS-15ºS, 150ºE-200ºE (Henley et al., 2015). +If different members and/or datasets are provided, the climatology (used to +calculate the anomalies) is computed individually for all of them. } \examples{ ## Observations or reanalyses -- GitLab From e705eafb1984f4473b287bd74ae5fafa7d82c494 Mon Sep 17 00:00:00 2001 From: aho Date: Fri, 9 Apr 2021 13:31:23 +0200 Subject: [PATCH 4/8] Add Filter() --- NAMESPACE | 1 + R/Filter.R | 120 +++++++++++++++++++++++++++++++++++ man/Filter.Rd | 51 +++++++++++++++ tests/testthat/test-Filter.R | 101 +++++++++++++++++++++++++++++ 4 files changed, 273 insertions(+) create mode 100644 R/Filter.R create mode 100644 man/Filter.Rd create mode 100644 tests/testthat/test-Filter.R diff --git a/NAMESPACE b/NAMESPACE index 341ae3a..ec36da5 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -21,6 +21,7 @@ export(ConfigShowSimilarEntries) export(ConfigShowTable) export(Corr) export(Eno) +export(Filter) export(GMST) export(GSAT) export(Histo2Hindcast) diff --git a/R/Filter.R b/R/Filter.R new file mode 100644 index 0000000..5293d54 --- /dev/null +++ b/R/Filter.R @@ -0,0 +1,120 @@ +#'Filter frequency peaks from an array +#' +#'Filter out the selected frequency from a time series. The filtering is +#'performed by dichotomy, seeking for a frequency around the parameter 'freq' +#'and the phase that maximizes the signal to subtract from the time series. +#'The maximization of the signal to subtract relies on a minimization of the +#'mean square differences between the time series ('data') and the cosine of +#'the specified frequency and phase. +#' +#'@param data A numeric vector or array of the data to be filtered. +#' If it's a vector, it should be a time series. If it's an array, +#' the dimensions must have at least 'time_dim'. +#'@param freq A number of the frequency to filter. +#'@param time_dim A character string indicating the dimension along which to +#' compute the filtering. The default value is 'ftime'. +#'@param ncores An integer indicating the number of cores to use for parallel +#' computation. The default value is NULL. +#' +#'@return A numeric vector or array of the filtered data with the dimensions +#' the same as 'data'. +#' +#'@examples +#'# Load sample data as in Load() example: +#'example(Load) +#'ensmod <- MeanDims(sampleData$mod, c(1, 2)) +#'spectrum <- Spectrum(ensmod) +#'for (jsdate in 1:dim(spectrum)['sdate']) { +#' for (jlen in 1:dim(spectrum)['ftime']) { +#' if (spectrum[jlen, 2, jsdate] > spectrum[jlen, 4, jsdate]) { +#' ensmod[jsdate, ] <- Filter(ensmod[jsdate, ], +#' spectrum[jlen, 1, jsdate]) +#' } +#' } +#'} +#' \donttest{ +#'PlotAno(InsertDim(ensmod, 2, 1), sdates = startDates) +#' } +#' +#'@importFrom stats lm +#'@export +Filter <- function(data, freq, time_dim = 'ftime', ncores = NULL) { + + # Check inputs + ## data + if (is.null(data)) { + stop("Parameter 'data' cannot be NULL.") + } + if (!is.numeric(data)) { + stop("Parameter 'data' must be a numeric array.") + } + if (is.null(dim(data))) { #is vector + dim(data) <- c(length(data)) + names(dim(data)) <- time_dim + } + if(any(is.null(names(dim(data))))| any(nchar(names(dim(data))) == 0)) { + stop("Parameter 'data' must have dimension names.") + } + ## freq + if (is.null(freq)) { + stop("Parameter 'freq' cannot be NULL.") + } + if (!is.numeric(freq) | length(freq) != 1) { + stop("Parameter 'freq' must be a number.") + } + ## time_dim + if (!is.character(time_dim) | length(time_dim) > 1) { + stop("Parameter 'time_dim' must be a character string.") + } + if (!time_dim %in% names(dim(data))) { + stop("Parameter 'time_dim' is not found in 'data' dimension.") + } + ## ncores + if (!is.null(ncores)) { + if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 | + length(ncores) > 1) { + stop("Parameter 'ncores' must be a positive integer.") + } + } + + ############################### + # Calculate Filter + output <- Apply(list(data), + target_dims = time_dim, + fun = .Filter, + freq = freq, + output_dims = time_dim, + ncores = ncores)$output1 + + return(output) +} + +.Filter <- function(data, freq) { + # data: [ftime] + + fac1 <- 1 + fac2 <- 1 + ndat <- length(data) + ndat2 <- length(which(!is.na(data))) + maxi <- 0 + endphase <- 0 + + for (jfreq in seq(freq - 0.5 / ndat2, freq + 0.5 / ndat2, 0.1 / (ndat2 * fac1))) { + for (phase in seq(0, pi, (pi / (10 * fac2)))) { + xtest <- cos(phase + c(1:ndat) * jfreq * 2 * pi) + test <- lm(data[is.na(data) == FALSE] ~ xtest[ + is.na(data) == FALSE])$fitted.values + if (sum(test ^ 2) > maxi) { + endphase <- phase + endfreq <- jfreq + } + maxi <- max(sum(test ^ 2), maxi) + } + } + xend <- cos(endphase + c(1:ndat) * endfreq * 2 * pi) + data[is.na(data) == FALSE] <- data[is.na(data) == FALSE] - lm( + data[is.na(data) == FALSE] ~ xend[is.na(data) == FALSE] + )$fitted.values + + return(invisible(data)) +} diff --git a/man/Filter.Rd b/man/Filter.Rd new file mode 100644 index 0000000..42f374c --- /dev/null +++ b/man/Filter.Rd @@ -0,0 +1,51 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Filter.R +\name{Filter} +\alias{Filter} +\title{Filter frequency peaks from an array} +\usage{ +Filter(data, freq, time_dim = "ftime", ncores = NULL) +} +\arguments{ +\item{data}{A numeric vector or array of the data to be filtered. +If it's a vector, it should be a time series. If it's an array, +the dimensions must have at least 'time_dim'.} + +\item{freq}{A number of the frequency to filter.} + +\item{time_dim}{A character string indicating the dimension along which to +compute the filtering. The default value is 'ftime'.} + +\item{ncores}{An integer indicating the number of cores to use for parallel +computation. The default value is NULL.} +} +\value{ +A numeric vector or array of the filtered data with the dimensions + the same as 'data'. +} +\description{ +Filter out the selected frequency from a time series. The filtering is +performed by dichotomy, seeking for a frequency around the parameter 'freq' +and the phase that maximizes the signal to subtract from the time series. +The maximization of the signal to subtract relies on a minimization of the +mean square differences between the time series ('data') and the cosine of +the specified frequency and phase. +} +\examples{ +# Load sample data as in Load() example: +example(Load) +ensmod <- MeanDims(sampleData$mod, c(1, 2)) +spectrum <- Spectrum(ensmod) +for (jsdate in 1:dim(spectrum)['sdate']) { + for (jlen in 1:dim(spectrum)['ftime']) { + if (spectrum[jlen, 2, jsdate] > spectrum[jlen, 4, jsdate]) { + ensmod[jsdate, ] <- Filter(ensmod[jsdate, ], + spectrum[jlen, 1, jsdate]) + } + } +} + \donttest{ +PlotAno(InsertDim(ensmod, 2, 1), sdates = startDates) + } + +} diff --git a/tests/testthat/test-Filter.R b/tests/testthat/test-Filter.R new file mode 100644 index 0000000..cf271e1 --- /dev/null +++ b/tests/testthat/test-Filter.R @@ -0,0 +1,101 @@ +context("s2dv::Filter tests") + +############################################## + # dat1 + set.seed(1) + dat1 <- array(rnorm(10), dim = c(dat = 1, ftime = 5, sdate = 2)) + freq1 <- 0.015 + # dat2 + set.seed(10) + dat2 <- c(1:10) + rnorm(10) + freq2 <- freq1 + + # dat3 + dat3 <- dat2 + dat3[2] <- NA + freq3 <- freq1 +############################################## + +test_that("1. Input checks", { + + # data + expect_error( + Filter(c()), + "Parameter 'data' cannot be NULL." + ) + expect_error( + Filter(c(NA, NA)), + "Parameter 'data' must be a numeric array." + ) + expect_error( + Filter(array(1:10, dim = c(2, 5))), + "Parameter 'data' must have dimension names." + ) + # freq + expect_error( + Filter(dat1, c()), + "Parameter 'freq' cannot be NULL." + ) + expect_error( + Filter(dat1, c(0.1, 0.2)), + "Parameter 'freq' must be a number." + ) + # time_dim + expect_error( + Filter(array(c(1:25), dim = c(dat = 1, date = 5, sdate = 5)), freq1), + "Parameter 'time_dim' is not found in 'data' dimension." + ) + expect_error( + Filter(dat1, freq1, time_dim = 2), + "Parameter 'time_dim' must be a character string." + ) + # ncores + expect_error( + Filter(dat1, freq1, ncore = 3.5), + "Parameter 'ncores' must be a positive integer." + ) +}) + +############################################## +test_that("2. Output checks: dat1", { + + expect_equal( + dim(Filter(dat1, freq1)), + c(ftime = 5, dat = 1, sdate = 2) + ) + expect_equal( + Filter(dat1, freq1)[, 1, 2], + c(-0.080093110, 0.141328669, -0.105230299, -0.004168101, 0.048162841), + tolerance = 0.0001 + ) + +}) + +############################################## +test_that("3. Output checks: dat2", { + + expect_equal( + dim(Filter(dat2, freq2)), + c(ftime = 10) + ) + expect_equal( + as.vector(Filter(dat2, freq2)[2:5]), + c(0.1215244, -1.0229749, -0.2053940, 0.7375181), + tolerance = 0.0001 + ) + +}) + +############################################## +test_that("4. Output checks: dat3", { + expect_equal( + dim(Filter(dat3, freq3)), + c(ftime = 10) + ) + expect_equal( + as.vector(Filter(dat3, freq3)[2:5]), + c(NA, -0.9414294, -0.1265448, 0.7910344), + tolerance = 0.0001 + ) + +}) -- GitLab From 9d86a39a5376d74c6ba36676077808cc5c8d29ac Mon Sep 17 00:00:00 2001 From: aho Date: Fri, 9 Apr 2021 13:47:26 +0200 Subject: [PATCH 5/8] Correct examples --- R/Filter.R | 8 ++++---- R/Spectrum.R | 9 ++++----- man/Filter.Rd | 8 ++++---- man/Spectrum.Rd | 9 ++++----- 4 files changed, 16 insertions(+), 18 deletions(-) diff --git a/R/Filter.R b/R/Filter.R index 5293d54..77efd11 100644 --- a/R/Filter.R +++ b/R/Filter.R @@ -22,13 +22,13 @@ #'@examples #'# Load sample data as in Load() example: #'example(Load) -#'ensmod <- MeanDims(sampleData$mod, c(1, 2)) +#'ensmod <- MeanDims(sampleData$mod, 2) #'spectrum <- Spectrum(ensmod) +#' #'for (jsdate in 1:dim(spectrum)['sdate']) { #' for (jlen in 1:dim(spectrum)['ftime']) { -#' if (spectrum[jlen, 2, jsdate] > spectrum[jlen, 4, jsdate]) { -#' ensmod[jsdate, ] <- Filter(ensmod[jsdate, ], -#' spectrum[jlen, 1, jsdate]) +#' if (spectrum[jlen, 2, 1, jsdate] > spectrum[jlen, 4, 1, jsdate]) { +#' ensmod[1, jsdate, ] <- Filter(ensmod[1, jsdate, ], spectrum[jlen, 1, 1, jsdate]) #' } #' } #'} diff --git a/R/Spectrum.R b/R/Spectrum.R index 6890536..6cb2113 100644 --- a/R/Spectrum.R +++ b/R/Spectrum.R @@ -24,14 +24,13 @@ #'@examples #'# Load sample data as in Load() example: #'example(Load) -#' -#'ensmod <- MeanDims(sampleData$mod, c(1, 2)) +#'ensmod <- MeanDims(sampleData$mod, 2) #'spectrum <- Spectrum(ensmod) +#' #'for (jsdate in 1:dim(spectrum)['sdate']) { #' for (jlen in 1:dim(spectrum)['ftime']) { -#' if (spectrum[jlen, 2, jsdate] > spectrum[jlen, 4, jsdate]) { -#' ensmod[jsdate, ] <- Filter(ensmod[jsdate, ], -#' spectrum[jlen, 1, jsdate]) +#' if (spectrum[jlen, 2, 1, jsdate] > spectrum[jlen, 4, 1, jsdate]) { +#' ensmod[1, jsdate, ] <- Filter(ensmod[1, jsdate, ], spectrum[jlen, 1, 1, jsdate]) #' } #' } #'} diff --git a/man/Filter.Rd b/man/Filter.Rd index 42f374c..814acf1 100644 --- a/man/Filter.Rd +++ b/man/Filter.Rd @@ -34,13 +34,13 @@ the specified frequency and phase. \examples{ # Load sample data as in Load() example: example(Load) -ensmod <- MeanDims(sampleData$mod, c(1, 2)) +ensmod <- MeanDims(sampleData$mod, 2) spectrum <- Spectrum(ensmod) + for (jsdate in 1:dim(spectrum)['sdate']) { for (jlen in 1:dim(spectrum)['ftime']) { - if (spectrum[jlen, 2, jsdate] > spectrum[jlen, 4, jsdate]) { - ensmod[jsdate, ] <- Filter(ensmod[jsdate, ], - spectrum[jlen, 1, jsdate]) + if (spectrum[jlen, 2, 1, jsdate] > spectrum[jlen, 4, 1, jsdate]) { + ensmod[1, jsdate, ] <- Filter(ensmod[1, jsdate, ], spectrum[jlen, 1, 1, jsdate]) } } } diff --git a/man/Spectrum.Rd b/man/Spectrum.Rd index 7f23d40..d1ebe9e 100644 --- a/man/Spectrum.Rd +++ b/man/Spectrum.Rd @@ -34,14 +34,13 @@ and the significance levels are estimated by a Monte-Carlo method. \examples{ # Load sample data as in Load() example: example(Load) - -ensmod <- MeanDims(sampleData$mod, c(1, 2)) +ensmod <- MeanDims(sampleData$mod, 2) spectrum <- Spectrum(ensmod) + for (jsdate in 1:dim(spectrum)['sdate']) { for (jlen in 1:dim(spectrum)['ftime']) { - if (spectrum[jlen, 2, jsdate] > spectrum[jlen, 4, jsdate]) { - ensmod[jsdate, ] <- Filter(ensmod[jsdate, ], - spectrum[jlen, 1, jsdate]) + if (spectrum[jlen, 2, 1, jsdate] > spectrum[jlen, 4, 1, jsdate]) { + ensmod[1, jsdate, ] <- Filter(ensmod[1, jsdate, ], spectrum[jlen, 1, 1, jsdate]) } } } -- GitLab From 2fa6de5c0bea630c1a32fae49995ee855310dc1d Mon Sep 17 00:00:00 2001 From: aho Date: Mon, 12 Apr 2021 10:19:20 +0200 Subject: [PATCH 6/8] Add @import multiApply --- R/Filter.R | 1 + R/Spectrum.R | 1 + 2 files changed, 2 insertions(+) diff --git a/R/Filter.R b/R/Filter.R index 77efd11..94472a2 100644 --- a/R/Filter.R +++ b/R/Filter.R @@ -36,6 +36,7 @@ #'PlotAno(InsertDim(ensmod, 2, 1), sdates = startDates) #' } #' +#'@import multiApply #'@importFrom stats lm #'@export Filter <- function(data, freq, time_dim = 'ftime', ncores = NULL) { diff --git a/R/Spectrum.R b/R/Spectrum.R index 6cb2113..019b421 100644 --- a/R/Spectrum.R +++ b/R/Spectrum.R @@ -38,6 +38,7 @@ #'PlotAno(InsertDim(ensmod, 2, 1), sdates = startDates) #' } #' +#'@import multiApply #'@importFrom stats spectrum cor rnorm sd quantile #'@export Spectrum <- function(data, time_dim = 'ftime', ncores = NULL) { -- GitLab From 3911f043bd4a31d0416996f75321cdd9ccee4924 Mon Sep 17 00:00:00 2001 From: aho Date: Thu, 6 May 2021 18:08:38 +0200 Subject: [PATCH 7/8] Syntax fix --- R/Spectrum.R | 2 +- man/Spectrum.Rd | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/R/Spectrum.R b/R/Spectrum.R index 019b421..7bf07be 100644 --- a/R/Spectrum.R +++ b/R/Spectrum.R @@ -2,7 +2,7 @@ #' #'Estimate the frequency spectrum of the data array together with its 95\% and #'99\% significance level. The output is provided as an array with dimensions -#"c(number of frequencies, 4). The column contains the frequency values, the +#'c(number of frequencies, 4). The column contains the frequency values, the #'power, the 95\% significance level and the 99\% one.\cr #'The spectrum estimation relies on a R built-in function \code{spectrum()} #'and the significance levels are estimated by a Monte-Carlo method. diff --git a/man/Spectrum.Rd b/man/Spectrum.Rd index d1ebe9e..d5e4ed5 100644 --- a/man/Spectrum.Rd +++ b/man/Spectrum.Rd @@ -27,6 +27,7 @@ A numeric array of the frequency spectrum with dimensions \description{ Estimate the frequency spectrum of the data array together with its 95\% and 99\% significance level. The output is provided as an array with dimensions +c(number of frequencies, 4). The column contains the frequency values, the power, the 95\% significance level and the 99\% one.\cr The spectrum estimation relies on a R built-in function \code{spectrum()} and the significance levels are estimated by a Monte-Carlo method. -- GitLab From 2e727394d592842d59496955fcee525b503986ae Mon Sep 17 00:00:00 2001 From: aho Date: Thu, 20 May 2021 21:05:18 +0200 Subject: [PATCH 8/8] Add param 'conf.lev' in Spectrum() --- R/Filter.R | 2 +- R/Spectrum.R | 36 +++++++++++++++++++++------------- man/Filter.Rd | 2 +- man/Spectrum.Rd | 25 +++++++++++++---------- tests/testthat/test-Spectrum.R | 12 +++++++++--- 5 files changed, 48 insertions(+), 29 deletions(-) diff --git a/R/Filter.R b/R/Filter.R index 94472a2..c4e76bf 100644 --- a/R/Filter.R +++ b/R/Filter.R @@ -27,7 +27,7 @@ #' #'for (jsdate in 1:dim(spectrum)['sdate']) { #' for (jlen in 1:dim(spectrum)['ftime']) { -#' if (spectrum[jlen, 2, 1, jsdate] > spectrum[jlen, 4, 1, jsdate]) { +#' if (spectrum[jlen, 2, 1, jsdate] > spectrum[jlen, 3, 1, jsdate]) { #' ensmod[1, jsdate, ] <- Filter(ensmod[1, jsdate, ], spectrum[jlen, 1, 1, jsdate]) #' } #' } diff --git a/R/Spectrum.R b/R/Spectrum.R index 7bf07be..2cbb167 100644 --- a/R/Spectrum.R +++ b/R/Spectrum.R @@ -1,11 +1,13 @@ #'Estimate frequency spectrum #' -#'Estimate the frequency spectrum of the data array together with its 95\% and -#'99\% significance level. The output is provided as an array with dimensions -#'c(number of frequencies, 4). The column contains the frequency values, the -#'power, the 95\% significance level and the 99\% one.\cr -#'The spectrum estimation relies on a R built-in function \code{spectrum()} -#'and the significance levels are estimated by a Monte-Carlo method. +#'Estimate the frequency spectrum of the data array together with a +#'user-specified confidence level. The output is provided as an array with +#'dimensions c(number of frequencies, stats = 3, other margin dimensions of +#'data). The 'stats' dimension contains the frequencies at which the spectral +#'density is estimated, the estimates of the spectral density, and the +#'significance level.\cr +#'The spectrum estimation relies on an R built-in function \code{spectrum()} +#'and the confidence interval is estimated by the Monte-Carlo method. #' #'@param data A vector or numeric array of which the frequency spectrum is #' required. If it's a vector, it should be a time series. If it's an array, @@ -13,13 +15,15 @@ #' evenly spaced in time. #'@param time_dim A character string indicating the dimension along which to #' compute the frequency spectrum. The default value is 'ftime'. +#'@param conf.lev A numeric indicating the confidence level for the Monte-Carlo +#' significance test. The default value is 0.95. #'@param ncores An integer indicating the number of cores to use for parallel #' computation. The default value is NULL. #' #'@return A numeric array of the frequency spectrum with dimensions -#' c( = number of frequencies, stats = 4, the rest of the +#' c( = number of frequencies, stats = 3, the rest of the #' dimensions of 'data'). The 'stats' dimension contains the frequency values, -#' the power, the 95\% significance level and the 99\% one. +#' the spectral density, and the confidence interval. #' #'@examples #'# Load sample data as in Load() example: @@ -29,7 +33,7 @@ #' #'for (jsdate in 1:dim(spectrum)['sdate']) { #' for (jlen in 1:dim(spectrum)['ftime']) { -#' if (spectrum[jlen, 2, 1, jsdate] > spectrum[jlen, 4, 1, jsdate]) { +#' if (spectrum[jlen, 2, 1, jsdate] > spectrum[jlen, 3, 1, jsdate]) { #' ensmod[1, jsdate, ] <- Filter(ensmod[1, jsdate, ], spectrum[jlen, 1, 1, jsdate]) #' } #' } @@ -41,7 +45,7 @@ #'@import multiApply #'@importFrom stats spectrum cor rnorm sd quantile #'@export -Spectrum <- function(data, time_dim = 'ftime', ncores = NULL) { +Spectrum <- function(data, time_dim = 'ftime', conf.lev = 0.95, ncores = NULL) { # Check inputs ## data @@ -65,6 +69,10 @@ Spectrum <- function(data, time_dim = 'ftime', ncores = NULL) { if (!time_dim %in% names(dim(data))) { stop("Parameter 'time_dim' is not found in 'data' dimension.") } + ## conf.lev + if (!is.numeric(conf.lev) | conf.lev < 0 | conf.lev > 1 | length(conf.lev) > 1) { + stop("Parameter 'conf.lev' must be a numeric number between 0 and 1.") + } ## ncores if (!is.null(ncores)) { if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 | @@ -80,12 +88,13 @@ Spectrum <- function(data, time_dim = 'ftime', ncores = NULL) { target_dims = time_dim, fun = .Spectrum, output_dims = c(time_dim, 'stats'), + conf.lev = conf.lev, ncores = ncores)$output1 return(output) } -.Spectrum <- function(data) { +.Spectrum <- function(data, conf.lev = 0.95) { # data: [time] data <- data[is.na(data) == FALSE] @@ -93,7 +102,7 @@ Spectrum <- function(data, time_dim = 'ftime', ncores = NULL) { if (ndat >= 3) { tmp <- spectrum(data, plot = FALSE) - output <- array(dim = c(length(tmp$spec), 4)) + output <- array(dim = c(length(tmp$spec), 3)) output[, 1] <- tmp$freq output[, 2] <- tmp$spec ntir <- 100 @@ -110,8 +119,7 @@ Spectrum <- function(data, time_dim = 'ftime', ncores = NULL) { store[jt, ] <- toto2$spec } for (jx in 1:length(tmp$spec)) { - output[jx, 3] <- quantile(store[, jx], 0.95) - output[jx, 4] <- quantile(store[, jx], 0.99) + output[jx, 3] <- quantile(store[, jx], conf.lev) } } else { output <- NA diff --git a/man/Filter.Rd b/man/Filter.Rd index 814acf1..f98fe0a 100644 --- a/man/Filter.Rd +++ b/man/Filter.Rd @@ -39,7 +39,7 @@ spectrum <- Spectrum(ensmod) for (jsdate in 1:dim(spectrum)['sdate']) { for (jlen in 1:dim(spectrum)['ftime']) { - if (spectrum[jlen, 2, 1, jsdate] > spectrum[jlen, 4, 1, jsdate]) { + if (spectrum[jlen, 2, 1, jsdate] > spectrum[jlen, 3, 1, jsdate]) { ensmod[1, jsdate, ] <- Filter(ensmod[1, jsdate, ], spectrum[jlen, 1, 1, jsdate]) } } diff --git a/man/Spectrum.Rd b/man/Spectrum.Rd index d5e4ed5..84b39c0 100644 --- a/man/Spectrum.Rd +++ b/man/Spectrum.Rd @@ -4,7 +4,7 @@ \alias{Spectrum} \title{Estimate frequency spectrum} \usage{ -Spectrum(data, time_dim = "ftime", ncores = NULL) +Spectrum(data, time_dim = "ftime", conf.lev = 0.95, ncores = NULL) } \arguments{ \item{data}{A vector or numeric array of which the frequency spectrum is @@ -15,22 +15,27 @@ evenly spaced in time.} \item{time_dim}{A character string indicating the dimension along which to compute the frequency spectrum. The default value is 'ftime'.} +\item{conf.lev}{A numeric indicating the confidence level for the Monte-Carlo +significance test. The default value is 0.95.} + \item{ncores}{An integer indicating the number of cores to use for parallel computation. The default value is NULL.} } \value{ A numeric array of the frequency spectrum with dimensions - c( = number of frequencies, stats = 4, the rest of the + c( = number of frequencies, stats = 3, the rest of the dimensions of 'data'). The 'stats' dimension contains the frequency values, - the power, the 95\% significance level and the 99\% one. + the spectral density, and the confidence interval. } \description{ -Estimate the frequency spectrum of the data array together with its 95\% and -99\% significance level. The output is provided as an array with dimensions -c(number of frequencies, 4). The column contains the frequency values, the -power, the 95\% significance level and the 99\% one.\cr -The spectrum estimation relies on a R built-in function \code{spectrum()} -and the significance levels are estimated by a Monte-Carlo method. +Estimate the frequency spectrum of the data array together with a +user-specified confidence level. The output is provided as an array with +dimensions c(number of frequencies, stats = 3, other margin dimensions of +data). The 'stats' dimension contains the frequencies at which the spectral +density is estimated, the estimates of the spectral density, and the +significance level.\cr +The spectrum estimation relies on an R built-in function \code{spectrum()} +and the confidence interval is estimated by the Monte-Carlo method. } \examples{ # Load sample data as in Load() example: @@ -40,7 +45,7 @@ spectrum <- Spectrum(ensmod) for (jsdate in 1:dim(spectrum)['sdate']) { for (jlen in 1:dim(spectrum)['ftime']) { - if (spectrum[jlen, 2, 1, jsdate] > spectrum[jlen, 4, 1, jsdate]) { + if (spectrum[jlen, 2, 1, jsdate] > spectrum[jlen, 3, 1, jsdate]) { ensmod[1, jsdate, ] <- Filter(ensmod[1, jsdate, ], spectrum[jlen, 1, 1, jsdate]) } } diff --git a/tests/testthat/test-Spectrum.R b/tests/testthat/test-Spectrum.R index 0a54101..caf53d3 100644 --- a/tests/testthat/test-Spectrum.R +++ b/tests/testthat/test-Spectrum.R @@ -38,6 +38,12 @@ test_that("1. Input checks", { Spectrum(dat1, time_dim = 2), "Parameter 'time_dim' must be a character string." ) + # conf.lev + expect_error( + Spectrum(dat1, conf.lev = -1), + "Parameter 'conf.lev' must be a numeric number between 0 and 1.", + fixed = T + ) # ncores expect_error( Spectrum(dat1, ncore = 3.5), @@ -50,7 +56,7 @@ test_that("2. Output checks: dat1", { expect_equal( dim(Spectrum(dat1)), - c(ftime = 2, stats = 4, dat = 1, sdate = 2) + c(ftime = 2, stats = 3, dat = 1, sdate = 2) ) expect_equal( Spectrum(dat1)[, 1, 1, 2], @@ -70,7 +76,7 @@ test_that("3. Output checks: dat2", { expect_equal( dim(Spectrum(dat2)), - c(ftime = 5, stats = 4) + c(ftime = 5, stats = 3) ) expect_equal( Spectrum(dat2)[, 1], @@ -89,7 +95,7 @@ test_that("3. Output checks: dat2", { test_that("4. Output checks: dat3", { expect_equal( dim(Spectrum(dat3)), - c(ftime = 4, stats = 4) + c(ftime = 4, stats = 3) ) expect_equal( Spectrum(dat3)[, 1], -- GitLab