From 2a6249527ca3af2232e3a3a451e6150d644594ae Mon Sep 17 00:00:00 2001 From: jkennel Date: Tue, 7 Jan 2020 09:07:01 -0500 Subject: [PATCH 01/15] handle multivariate input --- R/func_psdcore.R | 2 +- R/func_whiten.R | 14 +++++--------- 2 files changed, 6 insertions(+), 10 deletions(-) diff --git a/R/func_psdcore.R b/R/func_psdcore.R index d62679c..1fdce9d 100644 --- a/R/func_psdcore.R +++ b/R/func_psdcore.R @@ -60,7 +60,7 @@ psdcore <- function(X.d, ...) UseMethod("psdcore") #' @export psdcore.ts <- function(X.d, ...){ frq <- stats::frequency(X.d) - psdcore(as.vector(X.d), X.frq = frq, ...) + psdcore(as.matrix(X.d), X.frq = frq, ...) } #' @rdname psdcore diff --git a/R/func_whiten.R b/R/func_whiten.R index 5435a2a..ac80c80 100644 --- a/R/func_whiten.R +++ b/R/func_whiten.R @@ -143,19 +143,15 @@ prewhiten.ts <- function(tser, AR.max=0L, detrend=TRUE, demean=TRUE, impute=TRUE AR.max <- abs(AR.max) if (AR.max >= 0){ - # data.frame with fit params - fit.df <- data.frame(xr=base::seq_len(n.o), - xc=base::rep.int(1, n.o), - y=tser) - + X <- if (detrend){ if (verbose) message("\tdetrending (and demeaning)") - lmdfit <- stats::lm(y ~ xr, fit.df) - as.vector(stats::residuals(lmdfit)) + stats::residuals(stats::lm(tser~base::seq_len(n.o))) + } else if (demean) { if (verbose) message("\tdemeaning") - lmdfit <- stats::lm(y ~ xc, fit.df) - as.vector(stats::residuals(lmdfit)) + stats::residuals(stats::lm(tser~rep.int(1, n.o))) + } else { if (verbose) message("\tnothing was done to the timeseries object") tser From 2df9c0f107f15dd6c0d00c39519ec8fd925bca75 Mon Sep 17 00:00:00 2001 From: jkennel Date: Fri, 10 Jan 2020 14:28:40 -0500 Subject: [PATCH 02/15] addition of transfer function --- DESCRIPTION | 2 +- NAMESPACE | 8 + R/RcppExports.R | 63 +++++ R/func_coh_phase.R | 64 +++++ R/func_pilot.R | 13 +- R/func_psdcore.R | 116 ++++++--- R/func_pspectrum.R | 15 +- R/func_riedsid.R | 10 +- R/func_whiten.R | 12 +- man/coherence.Rd | 17 ++ man/det_vector.Rd | 17 ++ man/phase.Rd | 17 ++ man/pilot_spec.Rd | 3 + man/psdcore.Rd | 3 + man/pspectrum.Rd | 3 + man/resample_mvfft.Rd | 44 ++++ src/RcppExports.cpp | 56 ++++- src/resample_fft.cpp | 421 ++++++++++++++++++++++++++++---- tests/testthat/test-coh_phase.R | 9 + tests/testthat/test-riedsid.R | 26 ++ tests/testthat/test-spec-I.R | 117 ++++++++- tests/testthat/test-spec-II.R | 1 + tests/testthat/test-transfer.R | 48 ++++ 23 files changed, 984 insertions(+), 101 deletions(-) create mode 100644 R/func_coh_phase.R create mode 100644 man/coherence.Rd create mode 100644 man/det_vector.Rd create mode 100644 man/phase.Rd create mode 100644 man/resample_mvfft.Rd create mode 100644 tests/testthat/test-coh_phase.R create mode 100644 tests/testthat/test-transfer.R diff --git a/DESCRIPTION b/DESCRIPTION index e3c7b97..0780eb7 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -34,7 +34,7 @@ Suggests: RSEIS, rbenchmark, reshape2, - testthat + testthat (>= 2.1.0) Encoding: UTF-8 VignetteBuilder: knitr LinkingTo: Rcpp, RcppArmadillo diff --git a/NAMESPACE b/NAMESPACE index 68706b5..cabcd75 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -21,6 +21,7 @@ S3method(parabolic_weights,default) S3method(parabolic_weights,tapers) S3method(pgram_compare,amt) S3method(pilot_spec,default) +S3method(pilot_spec,matrix) S3method(pilot_spec,ts) S3method(plot,tapers) S3method(points,tapers) @@ -29,8 +30,10 @@ S3method(prewhiten,ts) S3method(print,summary.tapers) S3method(print,tapers) S3method(psdcore,default) +S3method(psdcore,matrix) S3method(psdcore,ts) S3method(pspectrum,default) +S3method(pspectrum,matrix) S3method(pspectrum,spec) S3method(pspectrum,ts) S3method(riedsid,default) @@ -52,6 +55,7 @@ S3method(varddiff,spec) export(.spec_confint) export(adapt_message) export(as.tapers) +export(coherence) export(colvec) export(constrain_tapers) export(create_poly) @@ -59,6 +63,7 @@ export(ctap_loess) export(ctap_simple) export(dB) export(data.frame.tapers) +export(det_vector) export(get_adapt_history) export(get_psd_env_name) export(get_psd_env_pointer) @@ -74,9 +79,11 @@ export(new_adapt_history) export(normalize) export(ones) export(parabolic_weights) +export(parabolic_weights_field) export(parabolic_weights_rcpp) export(parabolic_weights_rcpp2) export(pgram_compare) +export(phase) export(pilot_spec) export(prewhiten) export(psd_envAssign) @@ -92,6 +99,7 @@ export(pspectrum_basic) export(rcpp_ctap_simple) export(resample_fft_rcpp) export(resample_fft_rcpp2) +export(resample_mvfft) export(riedsid) export(riedsid2) export(rowvec) diff --git a/R/RcppExports.R b/R/RcppExports.R index cd0d975..9631dd8 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -129,3 +129,66 @@ riedsid_rcpp <- function(PSD, ntaper) { .Call('_psd_riedsid_rcpp', PACKAGE = 'psd', PSD, ntaper) } +#' @rdname parabolic_weights_field +#' +#' @param ntap the maximum number of tapers +#' +#' @export +parabolic_weights_field <- function(ntap) { + .Call('_psd_parabolic_weights_field', PACKAGE = 'psd', ntap) +} + +#' @title Resample an fft using varying numbers of sine tapers +#' +#' @description +#' Produce an un-normalized psd based on an fft and a vector of optimal sine +#' tapers. +#' +#' @details +#' To produce a psd estimate with our adaptive spectrum estimation method, +#' we need only make one fft calculation initially and then apply the weighting +#' factors given by \code{\link{parabolic_weights_field}}, which this function +#' does. +#' +#' @param fftz complex; a matrix representing the dual-length \code{\link{fft}}; see also the \code{dbl} argument +#' @param tapers integer; a vector of tapers +#' @param verbose logical; should messages be given? +#' @param dbl logical; should the code assume \code{fftz} is dual-length or single-length? +#' @param tapcap integer; the maximum number of tapers which can be applied; note that the length is +#' automatically limited by the length of the series. +#' +#' @return list that includes the auto and cross-spectral density, and the +#' number of tapers +#' +#' @seealso \code{\link{riedsid}} +#' +#' @examples +#' fftz <- complex(real=1:8, imaginary = 1:8) +#' taps <- 1:4 +#' try(resample_mvfft(fftz, taps)) +#' +#' @export +resample_mvfft <- function(fftz, tapers, verbose = TRUE, dbl = TRUE, tapcap = 10000L) { + .Call('_psd_resample_mvfft', PACKAGE = 'psd', fftz, tapers, verbose, dbl, tapcap) +} + +#' @title +#' det_vector +#' +#' @description +#' Determinant for an array +#' +#' @param x \code{numeric array} values to evaluate +#' +#' @return vector of determinants +#' +#' +#' @export +det_vector <- function(x) { + .Call('_psd_det_vector', PACKAGE = 'psd', x) +} + +solve_tf <- function(x) { + .Call('_psd_solve_tf', PACKAGE = 'psd', x) +} + diff --git a/R/func_coh_phase.R b/R/func_coh_phase.R new file mode 100644 index 0000000..4ae179f --- /dev/null +++ b/R/func_coh_phase.R @@ -0,0 +1,64 @@ +#' coherence +#' +#' calculate coherency from the spectra and cross-spectra +#' +#' @param pgram \code{numeric array} must be multivariate +#' +#' +#' @return list of coherency +#' +#' @export +#' +#' +coherence <- function(pgram) { + + dims = dim(pgram) + n_ser <- dims[2] + n_r <- dims[1] + + if (n_ser == 1) { + warning('Cannot calculate coherency for a single periodogram') + return(NA) + } + + coh <- matrix(NA_real_, nrow = n_r, ncol = (n_ser - 1)) + + for (i in 2L:(n_ser)) { + coh[, i-1] <- Re(Mod(pgram[, 1, i])^2 / (pgram[, 1, 1] * pgram[, i, i])) + } + + return(coh) + } + +#' phase +#' +#' calculate phase from the spectra and cross-spectra +#' +#' @param pgram \code{numeric array} must be multivariate +#' +#' +#' @return list of phase +#' +#' @export +#' +phase <- function(pgram) { + + dims = dim(pgram) + n_ser <- dims[2] + n_r <- dims[1] + + if (n_ser == 1) { + warning('Cannot calculate phase for a single periodogram') + return(NA) + } + + phase <- matrix(NA_real_, nrow = n_r, ncol = (n_ser - 1)) + + for (i in 2L:(n_ser)) { + + phase[, i-1] <- Arg(pgram[, 1, i]) + + } + + return(phase) +} \ No newline at end of file diff --git a/R/func_pilot.R b/R/func_pilot.R index 7d6ef06..f2a9065 100644 --- a/R/func_pilot.R +++ b/R/func_pilot.R @@ -62,10 +62,20 @@ pilot_spec <- function(x, ...) UseMethod("pilot_spec") #' @export pilot_spec.ts <- function(x, ...){ stopifnot(is.ts(x)) + # frq <- stats::frequency(x) + pilot_spec.default(x, ...) +} + + +#' @rdname pilot_spec +#' @aliases pilot_spec.matrix +#' @export +pilot_spec.matrix <- function(x, ...){ frq <- stats::frequency(x) - pilot_spec.default(as.vector(x), x.frequency=frq, ...) + pilot_spec(stats::ts(x, frequency=frq), ...) } + #' @rdname pilot_spec #' @export pilot_spec.default <- function(x, x.frequency=NULL, ntap=NULL, remove.AR=NULL, plot=FALSE, verbose=FALSE, fast = FALSE, ...){ @@ -84,7 +94,6 @@ pilot_spec.default <- function(x, x.frequency=NULL, ntap=NULL, remove.AR=NULL, p if (REMAR) remove.AR <- max(1, min(100, abs(remove.AR))) xprew <- prewhiten(x, x.fsamp=x.frequency, AR.max=remove.AR, detrend=TRUE, impute=TRUE, plot=FALSE, verbose=verbose) - ## Remove and AR model if (REMAR){ # AR fit diff --git a/R/func_psdcore.R b/R/func_psdcore.R index 1fdce9d..5f353be 100644 --- a/R/func_psdcore.R +++ b/R/func_psdcore.R @@ -60,7 +60,15 @@ psdcore <- function(X.d, ...) UseMethod("psdcore") #' @export psdcore.ts <- function(X.d, ...){ frq <- stats::frequency(X.d) - psdcore(as.matrix(X.d), X.frq = frq, ...) + psdcore.default(X.d, ...) +} + +#' @rdname psdcore +#' @aliases psdcore.matrix +#' @export +psdcore.matrix <- function(X.d, ...){ + frq <- stats::frequency(X.d) + psdcore(stats::ts(X.d, frequency=frq), ...) } #' @rdname psdcore @@ -88,7 +96,7 @@ psdcore.default <- function(X.d, # named series series <- deparse(substitute(X.d)) - + ## Convert to ts object if (!is.ts(X.d)){ if (is.null(X.frq)){ @@ -98,10 +106,10 @@ psdcore.default <- function(X.d, } X.d <- if (X.frq > 0){ # value represents sampling frequency - stats::ts(X.d, frequency=X.frq) + stats::ts(as.matrix(X.d), frequency=X.frq) } else if (X.frq < 0){ # value is sampling interval - stats::ts(X.d, deltat=abs(X.frq)) + stats::ts(as.matrix(X.d), deltat=abs(X.frq)) } else { stop("bad sampling information") } @@ -127,7 +135,7 @@ psdcore.default <- function(X.d, # initialize fft and other things, since this usually means a first run # # original series - n.o <- psd_envAssignGet(evars[['n.orig']], length(X.d)) + n.o <- psd_envAssignGet(evars[['n.orig']], NROW(X.d)) X <- psd_envAssignGet(evars[['series.orig']], { if (preproc){ # TODO: option for fast-detrend only, assign preproc flag in env (for plotting later) @@ -139,7 +147,7 @@ psdcore.default <- function(X.d, # Force series to be even in length (modulo division) n.e <- psd_envAssignGet(evars[['n.even']], modulo_floor(n.o)) - X.even <- ts(X[seq_len(n.e)], frequency=X.frq) + X.even <- ts(X[seq_len(n.e),], frequency=X.frq) X.even <- psd_envAssignGet(evars[['series.even']], X.even) stopifnot(is.ts(X.even)) @@ -147,8 +155,7 @@ psdcore.default <- function(X.d, nhalf <- psd_envAssignGet(evars[['n.even.half']], n.e/2) # variance of even series - varx <- psd_envAssignGet(evars[['var.even']], drop(stats::var(X.even))) - + varx <- psd_envAssignGet(evars[['var.even']], stats::var(X.even)) # create uniform tapers kseq <- psd_envAssignGet(evars[['last.taper']], { if (len_tapseq == 1){ @@ -160,16 +167,19 @@ psdcore.default <- function(X.d, } }) - ## zero pad and take double-length fft - padded <- as.numeric(c(X.even, rep.int(0, n.e))) + ## zero pad and take double-length fft + pad <- matrix(0.0, nrow = n.e, ncol = NCOL(X.even)) + padded <- rbind(as.matrix(X.even), pad) + # padded <- as.matrix(X.even) + + ## Calculate discrete Fourier tranform # Note fftw is faster for very long series but we are # using stats::fft until fftw is reliably built on CRAN - padded.fft <- psd_envAssignGet(evars[['fft.padded']], stats::fft(padded)) - + padded.fft <- psd_envAssignGet(evars[['fft.padded']], stats::mvfft(padded)) psd_envAssignGet(evars[['fft']], padded.fft) - + } else { if (verbose) warning("Working environment *not* refreshed. Results may be bogus.") @@ -194,7 +204,6 @@ psdcore.default <- function(X.d, ### Select frequencies for PSD evaluation f <- base::seq.int(0, nhalf, by=1) nfreq <- length(f) - ### Calculate the (un-normalized) PSD PSD <- psd_envAssignGet(evars[["last.psdcore"]], { @@ -206,8 +215,10 @@ psdcore.default <- function(X.d, ## resample fft with taper sequence and quadratic weighting # ( this is where the majority of the computational work is ) - kseq <- as.integer(kseq) - if(fast) { + kseq <- as.integer(kseq) + if(is.matrix(fftz)) { + reff <- resample_mvfft(fftz, tapers = kseq, verbose=verbose) + }else if(fast) { reff <- resample_fft_rcpp2(fftz, kseq, verbose=verbose) } else { reff <- resample_fft_rcpp(fftz, kseq, verbose=verbose) @@ -238,15 +249,15 @@ psdcore.default <- function(X.d, } }) - # should not be complex at this point! - stopifnot(!is.complex(PSD)) + # should not be complex at this point! - no longer true for cross-spectra -jrk + # stopifnot(!is.complex(PSD)) # there should not be any bad values here! pNAs <- is.na(PSD) if (any(pNAs)) warning("NA psd estimates?!") ## Nyquist frequencies - npsd <- length(PSD) + npsd <- NROW(PSD) frq <- as.numeric(base::seq.int(0, Nyq, length.out=npsd)) ## Update tapers for consistency @@ -258,16 +269,24 @@ psdcore.default <- function(X.d, # ( using the trapezoidal rule, the principal being that the # integrated spectrum should be equal to the variance of the signal ) # - trap.area <- base::sum(PSD, na.rm=TRUE) - mean(PSD[c(1,npsd)], na.rm=TRUE) - area.var.ratio <- varx / trap.area - PSD <- 2 * PSD * area.var.ratio * nhalf + + spec <- matrix(NA_real_, ncol = NCOL(PSD), nrow = NROW(PSD)) + for(i in 1:NCOL(PSD)) { + trap.area <- base::sum(PSD[, i, i, drop = FALSE], na.rm=TRUE) - mean(PSD[c(1,npsd), i, i, drop = FALSE], na.rm=TRUE) + area.var.ratio <- as.matrix(varx)[i,i] / trap.area + spec[, i] <- Re(2 * PSD[, i, i] * area.var.ratio * nhalf) + } ## Assemble final results mtap <- max(kseq, na.rm=TRUE) + if(NCOL(PSD) > 1) { + PSD.out <- list(freq = as.numeric(frq), - spec = as.numeric(PSD), - coh = NULL, - phase = NULL, + spec = spec, + pspec = PSD, + transfer = solve_tf(PSD), + coh = coherence(PSD), + phase = phase(PSD), kernel = NULL, # must be a scalar for plot.spec to give conf ints: df = 2 * mtap, # 2 DOF per taper, Percival and Walden eqn (370b) @@ -289,6 +308,35 @@ psdcore.default <- function(X.d, timebp = as.numeric(kseq/2), nyquist.frequency = Nyq ) + } else { + PSD.out <- list(freq = as.numeric(frq), + spec = as.numeric(spec), + pspec = NULL, + transfer = NULL, + coh = NULL, + phase = NULL, + kernel = NULL, + # must be a scalar for plot.spec to give conf ints: + df = 2 * mtap, # 2 DOF per taper, Percival and Walden eqn (370b) + numfreq = npsd, + ## bandwidth + # http://biomet.oxfordjournals.org/content/82/1/201.full.pdf + # half-width W = (K + 1)/{2(N + 1)} + # effective bandwidth ~ 2 W (accurate for many spectral windows) + bandwidth = (mtap + 1) / nhalf, + n.used = psd_envGet(evars[['n.even']]), + orig.n = psd_envGet(evars[['n.orig']]), + series = series, + snames = colnames(X), + method = sprintf("sine multitaper"), + taper = kseq, + pad = TRUE, # always! + detrend = preproc, # always true? + demean = preproc, + timebp = as.numeric(kseq/2), + nyquist.frequency = Nyq + ) + } class(PSD.out) <- c("amt","spec") if (plot) pgram_compare(PSD.out, ...) return(invisible(PSD.out)) @@ -400,13 +448,14 @@ pgram_compare.amt <- function(x, f=NULL, X=NULL, log.freq=TRUE, db.spec=TRUE, ta on.exit(par(opar)) #layout(matrix(c(1,2), ncol=1), c(1,2)) - layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow = TRUE)) + layout(matrix(c(1,1,2,2,3,4), nrow=3, ncol=2, byrow = TRUE)) par(mar=c(2, 3.5, 2.3, 1.2), oma=c(2,2,2,1)) par(las=1, tcl = -.3, mgp=c(2.2, 0.4, 0)) par(cex=0.8) ## Periodogram result first - plot(fc., sc., + + matplot(fc., sc., col="red", type="l", xaxs="i", xlim=flims, yaxs="i", ylim=slims, @@ -414,14 +463,13 @@ pgram_compare.amt <- function(x, f=NULL, X=NULL, log.freq=TRUE, db.spec=TRUE, ta ylab=slab, main="") # and adaptive result overlain - lines(f., s., type="l") + matlines(f., s., type="l", col = 'black') legend("bottomleft", c(sprintf("spec.pgram (%i%% cosine taper)", 100*tc.), sprintf("psdcore (%i - %i tapers)", min(t.), max(t.)) ), col=c("red","black"), lty=1, lwd=2, cex=0.9) mtext(flab, side=1, line=2) mtext("Spectral density estimates", adj=0, line=0.2, font=4) - ## Tapers if (is.tapers(t.)){ plot(t., f., xlim=flims, xaxs="i") @@ -434,13 +482,19 @@ pgram_compare.amt <- function(x, f=NULL, X=NULL, log.freq=TRUE, db.spec=TRUE, ta ## Original series main.pre <- ifelse(detrend | demean, "Modified e", "E") main.post <- sprintf("(dem. %s detr. %s)", detrend, demean) - plot(X, type="l", ylab="units", xlab="", xaxs="i", main="") + matplot(X, type="l", ylab="units", xlab="", xaxs="i", main="", col=c('black', 'blue', 'green', 'orange')) mtext(paste0(main.pre, "ven-length series"), line=1, font=4) mtext(main.post, cex=0.6) mtext(expression("time"), side=1, line=1.7) ## autocorrelation - acf(X, main="") + a <- acf(X, main="", plot = FALSE) + a_plot <- matrix(NA, ncol = NCOL(a$acf), nrow=NROW(a$acf)) + for(i in 1:NCOL(X)) { + a_plot[, i] <- a$acf[, i, i] + } + matplot(a$lag[,1,1], a_plot, type = 'hp', col=c('black', 'blue', 'green', 'orange'), ylab = "ACF") + mtext(expression("lag, time"), side=1, line=1.7) mtext("Auto-correlation function", line=1, font=4) diff --git a/R/func_pspectrum.R b/R/func_pspectrum.R index 2485e49..1677f64 100644 --- a/R/func_pspectrum.R +++ b/R/func_pspectrum.R @@ -43,9 +43,18 @@ pspectrum <- function(x, ...) UseMethod("pspectrum") #' @export pspectrum.ts <- function(x, ...){ frq <- stats::frequency(x) - psd::pspectrum(as.vector(x), x.frqsamp=frq, ...) + pspectrum.default(x, x.frqsamp=frq, ...) } +#' @rdname pspectrum +#' @aliases pspectrum.matrix +#' @export +pspectrum.matrix <- function(x, ...){ + frq <- stats::frequency(x) + pspectrum(stats::ts(x, frequency=frq), ...) +} + + #' @rdname pspectrum #' @aliases pspectrum.spec #' @export @@ -75,8 +84,9 @@ pspectrum.default <- function(x, verbose=TRUE, no.history=FALSE, plot=FALSE, ...){ - stopifnot(length(x)>1) + stopifnot(NROW(x)>1) + # plotting and iterations if (is.null(niter)) stopifnot(niter>=0) plotpsd_ <- FALSE @@ -124,7 +134,6 @@ pspectrum.default <- function(x, ## calculate optimal tapers kopt <- riedsid2(Pspec, verbose=rverb, fast = TRUE, ...) - # get data back for plotting, etc. if (stage==niter){ x <- psd_envGet("original_pspectrum_series") diff --git a/R/func_riedsid.R b/R/func_riedsid.R index 2c25b6d..268eab1 100644 --- a/R/func_riedsid.R +++ b/R/func_riedsid.R @@ -188,6 +188,9 @@ riedsid.default <- function(PSD, ntaper = 1L, return(kopt) } + + + #' @rdname riedsid #' @export riedsid2 <- function(PSD, ...) UseMethod("riedsid2") @@ -197,20 +200,23 @@ riedsid2 <- function(PSD, ...) UseMethod("riedsid2") riedsid2.spec <- function(PSD, ...){ pspec <- PSD[['spec']] freqs <- PSD[['freq']] + ntap <- if (is.amt(PSD)){ PSD[['taper']] } else { - rep.int(x=1L, times=length(pspec)) + rep.int(x=1L, times=NROW(pspec)) } + riedsid2(PSD=pspec, ntaper=ntap, ...) } + #' @rdname riedsid #' @export riedsid2.default <- function(PSD, ntaper=1L, constrained=TRUE, verbose=TRUE, fast=FALSE, ...){ if (fast) { - kopt <- riedsid_rcpp(PSD, ntaper) + kopt <- riedsid_rcpp(as.matrix(PSD), ntaper) } else { PSD <- as.vector(PSD) diff --git a/R/func_whiten.R b/R/func_whiten.R index ac80c80..2e91681 100644 --- a/R/func_whiten.R +++ b/R/func_whiten.R @@ -119,7 +119,7 @@ prewhiten.ts <- function(tser, AR.max=0L, detrend=TRUE, demean=TRUE, impute=TRUE stopifnot(is.ts(tser)) sps <- stats::frequency(tser) tstart <- stats::start(tser) - n.o <- length(tser) + n.o <- NROW(tser) ttime <- sps*n.o # NA action on input series @@ -146,11 +146,11 @@ prewhiten.ts <- function(tser, AR.max=0L, detrend=TRUE, demean=TRUE, impute=TRUE X <- if (detrend){ if (verbose) message("\tdetrending (and demeaning)") - stats::residuals(stats::lm(tser~base::seq_len(n.o))) + as.matrix(stats::residuals(stats::lm(tser~base::seq_len(n.o)))) } else if (demean) { if (verbose) message("\tdemeaning") - stats::residuals(stats::lm(tser~rep.int(1, n.o))) + as.matrix(stats::residuals(stats::lm(tser~rep.int(1, n.o)))) } else { if (verbose) message("\tnothing was done to the timeseries object") @@ -158,7 +158,7 @@ prewhiten.ts <- function(tser, AR.max=0L, detrend=TRUE, demean=TRUE, impute=TRUE } # TS object of residuals or equivalent - tser_prew_lm <- stats::ts(X, frequency=sps, start=tstart) + tser_prew_lm <- stats::ts(as.matrix(X), frequency=sps, start=tstart) } @@ -171,9 +171,9 @@ prewhiten.ts <- function(tser, AR.max=0L, detrend=TRUE, demean=TRUE, impute=TRUE # solve the Yule-Walker equations ardfit <- stats::ar.yw(tser_prew_lm, aic=TRUE, order.max=AR.max, demean=TRUE) - + # returns a TS object - tser_prew_ar <- ardfit[['resid']] + tser_prew_ar <- ts(as.matrix(ardfit[['resid']]), frequency=sps, start=tstart) if (impute) tser_prew_ar <- NAFUN(tser_prew_ar) } diff --git a/man/coherence.Rd b/man/coherence.Rd new file mode 100644 index 0000000..c6476b5 --- /dev/null +++ b/man/coherence.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/func_coh_phase.R +\name{coherence} +\alias{coherence} +\title{coherence} +\usage{ +coherence(pgram) +} +\arguments{ +\item{pgram}{\code{numeric array} must be multivariate} +} +\value{ +list of coherency +} +\description{ +calculate coherency from the spectra and cross-spectra +} diff --git a/man/det_vector.Rd b/man/det_vector.Rd new file mode 100644 index 0000000..f4d444e --- /dev/null +++ b/man/det_vector.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{det_vector} +\alias{det_vector} +\title{det_vector} +\usage{ +det_vector(x) +} +\arguments{ +\item{x}{\code{numeric array} values to evaluate} +} +\value{ +vector of determinants +} +\description{ +Determinant for an array +} diff --git a/man/phase.Rd b/man/phase.Rd new file mode 100644 index 0000000..1c0722f --- /dev/null +++ b/man/phase.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/func_coh_phase.R +\name{phase} +\alias{phase} +\title{phase} +\usage{ +phase(pgram) +} +\arguments{ +\item{pgram}{\code{numeric array} must be multivariate} +} +\value{ +list of phase +} +\description{ +calculate phase from the spectra and cross-spectra +} diff --git a/man/pilot_spec.Rd b/man/pilot_spec.Rd index dc49bd2..4d727ec 100644 --- a/man/pilot_spec.Rd +++ b/man/pilot_spec.Rd @@ -5,6 +5,7 @@ \alias{pilot_spectrum} \alias{spec.pilot} \alias{pilot_spec.ts} +\alias{pilot_spec.matrix} \alias{pilot_spec.default} \title{Calculate initial power spectral density estimates} \usage{ @@ -12,6 +13,8 @@ pilot_spec(x, ...) \method{pilot_spec}{ts}(x, ...) +\method{pilot_spec}{matrix}(x, ...) + \method{pilot_spec}{default}(x, x.frequency = NULL, ntap = NULL, remove.AR = NULL, plot = FALSE, verbose = FALSE, fast = FALSE, ...) diff --git a/man/psdcore.Rd b/man/psdcore.Rd index 283a861..ac83d57 100644 --- a/man/psdcore.Rd +++ b/man/psdcore.Rd @@ -3,6 +3,7 @@ \name{psdcore} \alias{psdcore} \alias{psdcore.ts} +\alias{psdcore.matrix} \alias{psdcore.default} \title{Multitaper power spectral density estimates of a series} \usage{ @@ -10,6 +11,8 @@ psdcore(X.d, ...) \method{psdcore}{ts}(X.d, ...) +\method{psdcore}{matrix}(X.d, ...) + \method{psdcore}{default}(X.d, X.frq = NULL, ntaper = as.tapers(5), preproc = TRUE, na.action = stats::na.fail, plot = FALSE, refresh = FALSE, verbose = FALSE, fast = FALSE, ndecimate, ...) diff --git a/man/pspectrum.Rd b/man/pspectrum.Rd index 06ae4ea..101e69b 100644 --- a/man/pspectrum.Rd +++ b/man/pspectrum.Rd @@ -3,6 +3,7 @@ \name{pspectrum} \alias{pspectrum} \alias{pspectrum.ts} +\alias{pspectrum.matrix} \alias{pspectrum.spec} \alias{pspectrum.default} \alias{pspectrum_basic} @@ -13,6 +14,8 @@ pspectrum(x, ...) \method{pspectrum}{ts}(x, ...) +\method{pspectrum}{matrix}(x, ...) + \method{pspectrum}{spec}(x, ...) \method{pspectrum}{default}(x, x.frqsamp = 1, ntap.init = NULL, diff --git a/man/resample_mvfft.Rd b/man/resample_mvfft.Rd new file mode 100644 index 0000000..d06755f --- /dev/null +++ b/man/resample_mvfft.Rd @@ -0,0 +1,44 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{resample_mvfft} +\alias{resample_mvfft} +\title{Resample an fft using varying numbers of sine tapers} +\usage{ +resample_mvfft(fftz, tapers, verbose = TRUE, dbl = TRUE, + tapcap = 10000L) +} +\arguments{ +\item{fftz}{complex; a matrix representing the dual-length \code{\link{fft}}; see also the \code{dbl} argument} + +\item{tapers}{integer; a vector of tapers} + +\item{verbose}{logical; should messages be given?} + +\item{dbl}{logical; should the code assume \code{fftz} is dual-length or single-length?} + +\item{tapcap}{integer; the maximum number of tapers which can be applied; note that the length is +automatically limited by the length of the series.} +} +\value{ +list that includes the auto and cross-spectral density, and the +number of tapers +} +\description{ +Produce an un-normalized psd based on an fft and a vector of optimal sine +tapers. +} +\details{ +To produce a psd estimate with our adaptive spectrum estimation method, +we need only make one fft calculation initially and then apply the weighting +factors given by \code{\link{parabolic_weights_field}}, which this function +does. +} +\examples{ +fftz <- complex(real=1:8, imaginary = 1:8) +taps <- 1:4 +try(resample_mvfft(fftz, taps)) + +} +\seealso{ +\code{\link{riedsid}} +} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 57d6295..d368ee8 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -83,17 +83,65 @@ BEGIN_RCPP END_RCPP } // riedsid_rcpp -arma::vec riedsid_rcpp(const arma::vec& PSD, const arma::ivec& ntaper); +arma::vec riedsid_rcpp(const arma::mat& PSD, const arma::ivec& ntaper); RcppExport SEXP _psd_riedsid_rcpp(SEXP PSDSEXP, SEXP ntaperSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< const arma::vec& >::type PSD(PSDSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type PSD(PSDSEXP); Rcpp::traits::input_parameter< const arma::ivec& >::type ntaper(ntaperSEXP); rcpp_result_gen = Rcpp::wrap(riedsid_rcpp(PSD, ntaper)); return rcpp_result_gen; END_RCPP } +// parabolic_weights_field +arma::field parabolic_weights_field(const int ntap); +RcppExport SEXP _psd_parabolic_weights_field(SEXP ntapSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const int >::type ntap(ntapSEXP); + rcpp_result_gen = Rcpp::wrap(parabolic_weights_field(ntap)); + return rcpp_result_gen; +END_RCPP +} +// resample_mvfft +List resample_mvfft(const arma::cx_mat& fftz, const arma::ivec& tapers, bool verbose, const bool dbl, const int tapcap); +RcppExport SEXP _psd_resample_mvfft(SEXP fftzSEXP, SEXP tapersSEXP, SEXP verboseSEXP, SEXP dblSEXP, SEXP tapcapSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::cx_mat& >::type fftz(fftzSEXP); + Rcpp::traits::input_parameter< const arma::ivec& >::type tapers(tapersSEXP); + Rcpp::traits::input_parameter< bool >::type verbose(verboseSEXP); + Rcpp::traits::input_parameter< const bool >::type dbl(dblSEXP); + Rcpp::traits::input_parameter< const int >::type tapcap(tapcapSEXP); + rcpp_result_gen = Rcpp::wrap(resample_mvfft(fftz, tapers, verbose, dbl, tapcap)); + return rcpp_result_gen; +END_RCPP +} +// det_vector +arma::cx_vec det_vector(const arma::cx_cube& x); +RcppExport SEXP _psd_det_vector(SEXP xSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::cx_cube& >::type x(xSEXP); + rcpp_result_gen = Rcpp::wrap(det_vector(x)); + return rcpp_result_gen; +END_RCPP +} +// solve_tf +arma::cx_mat solve_tf(arma::cx_cube x); +RcppExport SEXP _psd_solve_tf(SEXP xSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::cx_cube >::type x(xSEXP); + rcpp_result_gen = Rcpp::wrap(solve_tf(x)); + return rcpp_result_gen; +END_RCPP +} static const R_CallMethodDef CallEntries[] = { {"_psd_rcpp_ctap_simple", (DL_FUNC) &_psd_rcpp_ctap_simple, 2}, @@ -103,6 +151,10 @@ static const R_CallMethodDef CallEntries[] = { {"_psd_parabolic_weights_rcpp2", (DL_FUNC) &_psd_parabolic_weights_rcpp2, 1}, {"_psd_resample_fft_rcpp2", (DL_FUNC) &_psd_resample_fft_rcpp2, 5}, {"_psd_riedsid_rcpp", (DL_FUNC) &_psd_riedsid_rcpp, 2}, + {"_psd_parabolic_weights_field", (DL_FUNC) &_psd_parabolic_weights_field, 1}, + {"_psd_resample_mvfft", (DL_FUNC) &_psd_resample_mvfft, 5}, + {"_psd_det_vector", (DL_FUNC) &_psd_det_vector, 1}, + {"_psd_solve_tf", (DL_FUNC) &_psd_solve_tf, 1}, {NULL, NULL, 0} }; diff --git a/src/resample_fft.cpp b/src/resample_fft.cpp index d48fab3..647f0d6 100644 --- a/src/resample_fft.cpp +++ b/src/resample_fft.cpp @@ -118,6 +118,9 @@ List parabolic_weights_rcpp(const int ntap = 1) { return weights_out; } + + + //' @title Resample an fft using varying numbers of sine tapers //' //' @description @@ -463,16 +466,17 @@ List resample_fft_rcpp2( const arma::cx_vec& fftz, -arma::vec pad_data(const arma::vec& psd, +arma::mat pad_data(const arma::mat& psd, const int n, const int n_max, const double eps = 1e-78) { - arma::vec ii(n - 1 + (n_max) * 2); + int nc = psd.n_cols; + arma::mat ii(n - 1 + (n_max) * 2, nc); - ii.head(n_max) = arma::reverse(psd.head(n_max)); - ii.tail(n_max+1) = arma::reverse(psd.tail(n_max+1)); - ii(arma::span(n_max-1, n_max + n-2)) = psd; + ii.head_rows(n_max) = arma::reverse(psd.head_rows(n_max)); + ii.tail_rows(n_max+1) = arma::reverse(psd.tail_rows(n_max+1)); + ii(arma::span(n_max-1, n_max + n-2), arma::span::all) = psd; return(ii + eps); @@ -485,15 +489,17 @@ arma::vec pad_data(const arma::vec& psd, //' //' @return kopt vector // [[Rcpp::export]] -arma::vec riedsid_rcpp(const arma::vec& PSD, const arma::ivec& ntaper){ +arma::vec riedsid_rcpp(const arma::mat& PSD, const arma::ivec& ntaper){ double eps = 1e-78; double sc = 473.3736; double uzero, L, L2, CC; - int nf = PSD.n_elem; + int nf = PSD.n_rows; + int nc = PSD.n_cols; int nt = ntaper.n_elem; int j1, j2; + arma::vec out(nf); arma::ivec taper_vec(nf); // check for single length taper @@ -510,10 +516,10 @@ arma::vec riedsid_rcpp(const arma::vec& PSD, const arma::ivec& ntaper){ } int nadd = 1 + arma::max(nspan); - arma::vec y = arma::log(pad_data(PSD, nf, nadd, eps = 1e-78)); + arma::mat y = arma::log(pad_data(PSD, nf, nadd, eps = 1e-78)); double dy, d2y; - arma::vec yders(nf); + arma::mat yders(nf, nc); for (int j = 0; j < nf; j++) { @@ -527,78 +533,393 @@ arma::vec riedsid_rcpp(const arma::vec& PSD, const arma::ivec& ntaper){ uzero = (L2 - 1) / CC; - // first deriv - dy = as_scalar(u.t() * y(arma::span(j1, j2)) * CC / (L*(L2 - 1))); - // second deriv - d2y = as_scalar((u % u - uzero).t() * - y(arma::span(j1, j2)) * 360 / (L * (L2 - 1) * (L2 - 4))); + for (int i = 0; i < nc; i++) { + + // first deriv + dy = as_scalar(u.t() * y(arma::span(j1, j2), i) * CC / (L*(L2 - 1))); + + // second deriv + d2y = as_scalar((u % u - uzero).t() * + y(arma::span(j1, j2), i) * 360 / (L * (L2 - 1) * (L2 - 4))); + + yders(j, i) = fabs(dy*dy + d2y + eps); + } - yders(j) = fabs(dy*dy + d2y + eps); + out(j) = arma::max(yders.row(j)); } - return(arma::round(pow(sc, 0.2) / arma::pow(yders, 0.4))); + return(arma::round(pow(sc, 0.2) / arma::pow(out, 0.4))); } -/*** R -library(microbenchmark) -library(psd) -n. <- 10000 -set.seed(1234) -x <- cumsum(sample(c(-1, 1), n., TRUE)) -fftz <- fft(x) -taps <- ceiling(runif(n.,10,30)) +// This is for multivariate series -------------------------------------------- -microbenchmark( - rsz1 <- resample_fft_rcpp(fftz, taps, verbose = FALSE), - rsz2 <- resample_fft_rcpp2(fftz, taps, verbose = FALSE), - times = 10 -) -all.equal(rsz1, rsz2) +// For all but short series this should be faster +//' @rdname parabolic_weights_field +//' +//' @param ntap the maximum number of tapers +//' +//' @export +// [[Rcpp::export]] +arma::field parabolic_weights_field(const int ntap) { + + // + // return quadratic spectral weighting factors for a given number of tapers + // Barbour and Parker (2014) Equation 7 + // + arma::field f(ntap, 1); + arma::vec kseq = arma::pow(arma::regspace(0, ntap - 1), 2); + + double t3; + for (int i = 1; i < ntap+1; i++) { + t3 = log(i * (i - 0.25) * (i + 1.0)); + f(i-1,0) = arma::exp(log(1.5) + arma::log(i * i - kseq(arma::span(0, i-1))) - t3); + } + + return(f); + +} -microbenchmark( - kopt1 <- riedsid2(PSD = rsz1$psd, ntaper = taps, constrained = FALSE, verbose = FALSE), - kopt2 <- riedsid2(PSD = rsz1$psd, ntaper = taps, constrained = FALSE, verbose = FALSE, fast = TRUE), - times = 10 -) +//' @title Resample an fft using varying numbers of sine tapers +//' +//' @description +//' Produce an un-normalized psd based on an fft and a vector of optimal sine +//' tapers. +//' +//' @details +//' To produce a psd estimate with our adaptive spectrum estimation method, +//' we need only make one fft calculation initially and then apply the weighting +//' factors given by \code{\link{parabolic_weights_field}}, which this function +//' does. +//' +//' @param fftz complex; a matrix representing the dual-length \code{\link{fft}}; see also the \code{dbl} argument +//' @param tapers integer; a vector of tapers +//' @param verbose logical; should messages be given? +//' @param dbl logical; should the code assume \code{fftz} is dual-length or single-length? +//' @param tapcap integer; the maximum number of tapers which can be applied; note that the length is +//' automatically limited by the length of the series. +//' +//' @return list that includes the auto and cross-spectral density, and the +//' number of tapers +//' +//' @seealso \code{\link{riedsid}} +//' +//' @examples +//' fftz <- complex(real=1:8, imaginary = 1:8) +//' taps <- 1:4 +//' try(resample_mvfft(fftz, taps)) +//' +//' @export +// [[Rcpp::export]] +List resample_mvfft( const arma::cx_mat& fftz, + const arma::ivec& tapers, + bool verbose = true, + const bool dbl = true, + const int tapcap = 10000 ) { + + // + // resample and reweight an fft estimates for a given number of tapers + // using quadratic scaling in Barbour and Parker (2014) + // + // + // needs: + // - fftz: complex vector -- the FFT of the original signal + // - tapers: integer vector -- the number of tapers at each frequency + // + // optional args: + // - verbose: logical -- should warnings be given? + // - dbl: logical -- should the progam assume 'fftz' is for a double-length (padded) series? + // Otherwise a single-length series is assumed. + // - tapcap: integer -- the maximum number of tapers at any frequency + // + + int sc, nf, nc, nt, ne, ne2, nhalf, m2, mleft1, mleft2, Kc, ki, j1, j2; + double wi; + + + arma::cx_double fdc1, fdc2; + + + if (dbl){ + // double-length fft estimates assumed by default + sc = 2; + } else { + // but could be single-length + sc = 1; + } + + // even, double, and half lengths + nf = fftz.n_rows / sc; + nc = fftz.n_cols; + nt = tapers.n_elem; + ne = nf - (nf % 2); + + if (verbose){ + Function msg("message"); + msg(std::string("\tfft resampling")); + } + + if (ne < nf){ + warning("fft was not done on an even length series"); + } + + ne2 = 2 * ne; + nhalf = ne / 2; + + arma::ivec taper_vec(nhalf); + + if (nhalf < 1){ + stop("cannot operate on length-1 series"); + } + + + if (nt == 1){ + warning("forced taper length"); + taper_vec.fill(tapers[0]); + } else { + taper_vec = tapers; + } + + + // set the current number of tapers, limited by a few factors + arma::uvec wh = arma::find(taper_vec > nhalf); + taper_vec(wh).fill(nhalf); + wh = arma::find(taper_vec > tapcap); + taper_vec(wh).fill(tapcap); + wh = arma::find(taper_vec <= 0); + taper_vec(wh).ones(); + + + int mm = taper_vec.max(); + + arma::vec w(mm); + arma::field para = parabolic_weights_field(mm); + + + // + // Calculate the psd by averaging over tapered estimates + // + + + arma::cx_cube psd(nhalf, nc, nc); + psd.fill(arma::cx_double(0.0, 0.0)); + + // each frequency + for (int jj = 0; jj < nhalf; jj++) { + + m2 = 2 * jj; + // number of tapers applied at a given frequency (do not remove m+1 index!) + + Kc = taper_vec[jj]; // orig: m+1, but this was leading to an indexing bug + + // taper sequence and spectral weights + w = para(Kc-1, 0); + + // scan across ki to get vector of + // spectral differences based on modulo indices + for (int ik = 0; ik < Kc; ik++) { + + wi = w[ik]; + ki = ik + 1; + + mleft1 = m2 + ne2 - ki; + mleft2 = m2 + ki; + + j1 = mleft1 % ne2; + j2 = mleft2 % ne2; + + // sum as loop progresses, auto and cross-spectra + + for (int ii = 0; ii < nc; ii++) { + for (int kk = ii; kk < nc; kk++) { + + fdc1 = fftz(j1, ii) - fftz(j2, ii); + fdc2 = fftz(j1, kk) - fftz(j2, kk); + + psd(jj, ii, kk) += (fdc1 * std::conj(fdc2)) * wi; + } + } + } + } + + // Add lower triangle for use in transfer function calculation + for (int ii = 0; ii < nc; ii++) { + for (int kk = 0; kk < nc; kk++) { + if (kk < ii){ + psd(arma::span(), arma::span(ii), arma::span(kk)) = + arma::conj(psd(arma::span(), arma::span(kk), arma::span(ii))); + } + } + } + + // return list to match previous function definition - jrk + return Rcpp::List::create( + Named("freq.inds") = arma::regspace(1, nhalf), + Named("k.capped") = taper_vec(arma::span(0, nhalf-1)), + Named("psd") = psd + ); +} -all.equal(kopt1, kopt2) -microbenchmark( - psd1 <- psdcore(x, verbose = FALSE), - psd2 <- psdcore(x, verbose = FALSE, fast = TRUE), - times = 10 -) +// solve for transfer function -all.equal(psd1, psd2) +//' @title +//' det_vector +//' +//' @description +//' Determinant for an array +//' +//' @param x \code{numeric array} values to evaluate +//' +//' @return vector of determinants +//' +//' +//' @export +// [[Rcpp::export]] +arma::cx_vec det_vector(const arma::cx_cube& x) { + + std::size_t n = x.n_rows; + std::size_t nc = x.n_cols; + + arma::cx_vec output(n); + arma::cx_mat mat(nc, nc); + + for (std::size_t i = 0; i < n; i++){ + mat = x(arma::span(i), arma::span::all, arma::span::all); + output(i) = arma::det(mat); + } + + return(output); + +} +// [[Rcpp::export]] +arma::cx_mat solve_tf(arma::cx_cube x) { + + unsigned int n = x.n_rows; + unsigned int nr = x.n_cols -1; + + arma::cx_mat output(n, nr); + arma::cx_cube numer_base = x(arma::span::all, + arma::span(1, nr), + arma::span(1, nr));; + arma::cx_cube numer(n, nr, nr); + + + arma::cx_cube numer_sol = x(arma::span::all, + arma::span(1, nr), + arma::span(0)); + + arma::cx_vec denom = det_vector(numer_base); + + for (arma::uword i=0; i < nr; i++){ + + // fill new column + numer = numer_base; + numer(arma::span::all, arma::span::all, arma::span(i)) = numer_sol; + + output.col(i) = det_vector(numer) / denom; + + } + + return(output); + +} + + + + +/*** R + +library(microbenchmark) +library(psd) + +n. <- 1000000 +set.seed(1234) +nc <- 2 +x <- matrix(cumsum(sample(c(-1, 1), n., TRUE)), ncol=nc) +fftz <- mvfft(x) +psd <- Re(fftz * Conj(fftz)) +taps <- ceiling(runif(n./nc,10,300)) + microbenchmark( - pspec1 <- pspectrum(x, verbose = FALSE), - pspec2 <- pspectrum(x, verbose = FALSE, fast = TRUE), - times = 10 + a1 <- riedsid_rcpp(as.matrix(psd)[, 1, drop = FALSE], taps), + a2 <- riedsid_rcpp(psd[, 2, drop = FALSE], taps), + b <- riedsid_rcpp(as.matrix(psd), taps), + times = 2 ) -all.equal(pspec1, pspec2) +all.equal(apply(cbind(a1,a2), 1, min), b) microbenchmark( - pilot1 <- pilot_spec(x, verbose = FALSE), - pilot2 <- pilot_spec(x, verbose = FALSE, fast = TRUE), - times = 10 + rsz1 <- resample_fft_rcpp(fftz[,2], taps, verbose = FALSE), + rsz2 <- resample_fft_rcpp2(fftz[,2], taps, verbose = FALSE), + rsz3 <- resample_mvfft(fftz, taps, verbose = FALSE), + times = 2 ) -all.equal(pilot1, pilot2) +all.equal(rsz1, rsz2) +all.equal(rsz1$psd, Re(rsz3$psd[,2,2])) + + +# +# microbenchmark( +# tf <- solve_tf(rsz3$psd), +# times = 5 +# ) + + +# +# microbenchmark( +# kopt1 <- riedsid2(PSD = rsz1$psd, ntaper = taps, constrained = FALSE, verbose = FALSE), +# kopt2 <- riedsid2(PSD = rsz1$psd, ntaper = taps, constrained = FALSE, verbose = FALSE, fast = TRUE), +# times = 10 +# ) + + +# +# all.equal(kopt1, kopt2) +# +# +# microbenchmark( +# psd1 <- psdcore(x, verbose = FALSE), +# psd2 <- psdcore(x, verbose = FALSE, fast = TRUE), +# times = 10 +# ) +# +# +# all.equal(psd1, psd2) +# +# +# +# microbenchmark( +# pspec1 <- pspectrum(x, verbose = FALSE), +# pspec2 <- pspectrum(x, verbose = FALSE, fast = TRUE), +# times = 10 +# ) +# +# all.equal(pspec1, pspec2) +# +# +# microbenchmark( +# pilot1 <- pilot_spec(x, verbose = FALSE), +# pilot2 <- pilot_spec(x, verbose = FALSE, fast = TRUE), +# times = 10 +# ) +# +# all.equal(pilot1, pilot2) */ diff --git a/tests/testthat/test-coh_phase.R b/tests/testthat/test-coh_phase.R new file mode 100644 index 0000000..943c71b --- /dev/null +++ b/tests/testthat/test-coh_phase.R @@ -0,0 +1,9 @@ +test_that("coherence and phase works for single entry", { + + a <- array(rnorm(100), dim = c(100, 1, 1)) + + expect_warning(phase(a)) + expect_warning(coherence(a)) + + +}) diff --git a/tests/testthat/test-riedsid.R b/tests/testthat/test-riedsid.R index 71539d0..c0d0ad7 100644 --- a/tests/testthat/test-riedsid.R +++ b/tests/testthat/test-riedsid.R @@ -38,3 +38,29 @@ test_that("riedsid2 R-version is equal to Rcpp version",{ expect_equal(riedsid2(pa_b, fast=FALSE), riedsid2(pa_b, fast = TRUE)) }) + + +test_that("multivariate riedsid2 works",{ + + set.seed(1234) + x <- matrix(rnorm(200), ncol = 2) + taps <- ceiling(runif(200/2, 10, 300)) + + pd <- stats::spectrum(x, plot=FALSE) + + # each separately and then take the minimum number of tapers + r_s <- cbind(riedsid2(pd$spec[, 1], fast=FALSE), + riedsid2(pd$spec[, 2], fast=FALSE)) + r_s <- apply(r_s, 1, min) + + # multivariate method + r_mv <- riedsid2(pd$spec, fast=TRUE) + expect_equal(r_mv, r_s) + + # spec method works + r_mv_spec <- riedsid2(pd, fast=TRUE) + expect_equal(r_mv_spec, r_s) + + +}) + \ No newline at end of file diff --git a/tests/testthat/test-spec-I.R b/tests/testthat/test-spec-I.R index 2454e93..fd5b59a 100644 --- a/tests/testthat/test-spec-I.R +++ b/tests/testthat/test-spec-I.R @@ -79,12 +79,12 @@ test_that("pspectrum results are accurate",{ test_that("psdcore arguments are tested",{ set.seed(1234) x <- rnorm(100) - xp1 <- psdcore.default(X.d = x, X.frq = 1, plot = FALSE) - xp2 <- psdcore.default(X.d = x, X.frq = -1, plot = FALSE) + xp1 <- psdcore(X.d = x, X.frq = 1, plot = FALSE) + xp2 <- psdcore(X.d = x, X.frq = -1, plot = FALSE) expect_is(xp1, 'spec') expect_is(xp2, 'spec') expect_equal(xp1,xp2) - expect_error(psdcore.default(X.d = x, X.frq = "1", plot = FALSE)) + expect_error(psdcore(X.d = x, X.frq = "1", plot = FALSE)) }) test_that("psdcore results are accurate",{ @@ -205,4 +205,113 @@ test_that("check fast version",{ pilot_spec(xt2, verbose = FALSE, plot = FALSE, fast = TRUE)) }) -## + +test_that("check multivariate autospectra for psdcore",{ + + set.seed(1234) + x1 <- rnorm(100) + x2 <- x1 * 2 + x <- cbind(x1, x2) + ps <- psdcore(x, verbose = FALSE, plot = FALSE, fast = TRUE) + + ps1 <- psdcore(x[,1], verbose = FALSE, plot = FALSE, fast = TRUE) + ps2 <- psdcore(x[,2], verbose = FALSE, plot = FALSE, fast = TRUE) + + expect_equal(ps$spec[,1], ps1$spec) + expect_equal(ps$spec[,2], ps2$spec) + + + ps1 <- psdcore(x[,1], verbose = FALSE, plot = FALSE, fast = FALSE) + ps2 <- psdcore(x[,2], verbose = FALSE, plot = FALSE, fast = FALSE) + + expect_equal(ps$spec[,1], ps1$spec) + expect_equal(ps$spec[,2], ps2$spec) + + # ps <- pspectrum(x) +}) + + + +test_that("check multivariate autospectra for pilot_spec",{ + + set.seed(1234) + x <- matrix(rnorm(200), ncol = 2) + + ps1 <- pilot_spec(x[,1], verbose = FALSE, plot = FALSE, fast = FALSE) + ps2 <- pilot_spec(x[,2], verbose = FALSE, plot = FALSE, fast = FALSE) + + ps <- pilot_spec(x, verbose = FALSE, plot = FALSE, fast = TRUE) + + expect_equal(ps$spec[,1], ps1$spec) + expect_equal(ps$spec[,2], ps2$spec) + + + ps <- pspectrum(x, plot = FALSE) + +}) + + +test_that("check multivariate autospectra, coherence, and phase for psdcore",{ + + # compare spec.pgram to psdcore - the methods aren't equivalent but give + # similar results + + set.seed(1234) + x <- matrix(rnorm(2000), ncol = 2) + + pd <- spec.pgram(x, plot = FALSE, spans = 5) + pd <- normalize(pd, 1, "spectrum", verbose = FALSE) + pc <- psdcore(x, plot = FALSE, verbose = FALSE, ntaper = as.tapers(5)) + + # is there an indexing issue here, need to remove the first value in psdcore + # to get alignment + + pc$coh <- pc$coh[-1] + pc$phase <- pc$phase[-1] + pc$spec <- pc$spec[-1, ] + + n <- length(pd$coh) + pd$coh <- pd$coh[-n] + pd$phase <- pd$phase[-n] + pd$spec <- pd$spec[-n, ] + + + # select narrow range to avoid wrap around issues for phase + sub_range <- 1.2 + wh1 <- which(pd$phase < pi/sub_range & pd$phase > -pi/sub_range) + wh2 <- which(pc$phase < pi/sub_range & pc$phase > -pi/sub_range) + wh <- intersect(wh1, wh2) + + # do regression between methods + phase_coef <- coefficients(lm(pd$phase[wh]~pc$phase[wh])) + coh_coef <- coefficients(lm(pd$coh~pc$coh)) + spec_coef_1 <- coefficients(lm(pd$spec[,1]~pc$spec[,1])) + spec_coef_2 <- coefficients(lm(pd$spec[,2]~pc$spec[,2])) + + + # check intercept + expect_lte(abs(phase_coef[1]), 0.05) + expect_lte(abs(coh_coef[1]), 0.05) + expect_lte(abs(spec_coef_1[1]), 0.05) + expect_lte(abs(spec_coef_2[1]), 0.05) + + + # check slope + expect_lte(abs(phase_coef[2]-1), 0.05) + expect_lte(abs(coh_coef[2]-1), 0.05) + expect_lte(abs(spec_coef_1[2]-1), 0.05) + expect_lte(abs(spec_coef_2[2]-1), 0.05) + + + # plot(pd$spec[,1]) + # points(pc$spec[,1], type='l') + # plot(pd$spec[,2]) + # points(pc$spec[,2], type='l') + # plot(pd$coh) + # points(pc$coh, type='l') + # plot(pd$phase) + # points(pc$phase, type='l') + + +}) + diff --git a/tests/testthat/test-spec-II.R b/tests/testthat/test-spec-II.R index 5cfb1e4..002225d 100644 --- a/tests/testthat/test-spec-II.R +++ b/tests/testthat/test-spec-II.R @@ -18,4 +18,5 @@ test_that("'updating' spectrum is not yet implemented",{ pa_b <- pspectrum_basic(x, verbose = FALSE) expect_error(pspectrum(pa_b, plot = FALSE, verbose = FALSE)) + }) diff --git a/tests/testthat/test-transfer.R b/tests/testthat/test-transfer.R new file mode 100644 index 0000000..01810eb --- /dev/null +++ b/tests/testthat/test-transfer.R @@ -0,0 +1,48 @@ +test_that("transfer function works", { + + # can we apply and recover a constant response + set.seed(1234) + seq_length <- 10000 + x <- sin(seq(0, 8*pi, length.out = seq_length)) + rnorm(seq_length) + y <- x * 0.5 + xy <- cbind(y, x) + pc <- psdcore(xy, plot = FALSE, verbose = FALSE, ntaper = as.tapers(3)) + + expect_equal(Mod(pc$transfer), + matrix(0.5, nrow = length(pc$transfer), ncol = 1)) + + + # can we apply and recover a variable response + # generate data + set.seed(1234) + kern_length <- 101 + seq_length <- 10000 + x <- sin(seq(0, 8*pi, length.out = seq_length)) + rnorm(seq_length) + kern <- exp(seq(-10, 0, length.out = kern_length)) / kern_length + pad_kern <- c(rep(0, seq_length-kern_length), kern) + # convolution + y <- Re(fft(fft(x) * Conj(fft(pad_kern)), inverse = TRUE)) / seq_length + + # plot(x, type='l') + # points(y, col = 'red', type='l') + + xy <- cbind(y, x) + pc <- psdcore(xy, plot = FALSE, verbose = FALSE, ntaper = as.tapers(3)) + + # convert to impulse response + # this is necessary for the inverse fft as we use half length + pc$transfer <- c(pc$transfer, rev(Conj(pc$transfer))) + + # inverse fft of transfer fun + n <- floor(NROW(pc$transfer) / 2) + 1 + imp <- fft(pc$transfer, inverse = TRUE) / NROW(pc$transfer) + imp <- Mod(imp) * sign(Re(imp)) + mt_impulse <- cumsum(rev(imp)[1:n] + (imp)[1:n]) + + expect_equal(cumsum(rev(kern)), mt_impulse[1:kern_length], tolerance = 1e-3) + + # plot(cumsum(rev(kern)), type='l', col = 'red') + # points(mt_impulse[1:kern_length], cex = 0.5, pch = 20) + + +}) From e3fb1552e5b258b4df4559a44224234edc0bc904 Mon Sep 17 00:00:00 2001 From: jkennel Date: Fri, 10 Jan 2020 14:56:51 -0500 Subject: [PATCH 03/15] add coherence and phase calculation --- R/func_coherence_phase.R | 64 +++++++++++++++++++++++++++ tests/testthat/test-coherence_phase.R | 9 ++++ 2 files changed, 73 insertions(+) create mode 100644 R/func_coherence_phase.R create mode 100644 tests/testthat/test-coherence_phase.R diff --git a/R/func_coherence_phase.R b/R/func_coherence_phase.R new file mode 100644 index 0000000..4ae179f --- /dev/null +++ b/R/func_coherence_phase.R @@ -0,0 +1,64 @@ +#' coherence +#' +#' calculate coherency from the spectra and cross-spectra +#' +#' @param pgram \code{numeric array} must be multivariate +#' +#' +#' @return list of coherency +#' +#' @export +#' +#' +coherence <- function(pgram) { + + dims = dim(pgram) + n_ser <- dims[2] + n_r <- dims[1] + + if (n_ser == 1) { + warning('Cannot calculate coherency for a single periodogram') + return(NA) + } + + coh <- matrix(NA_real_, nrow = n_r, ncol = (n_ser - 1)) + + for (i in 2L:(n_ser)) { + coh[, i-1] <- Re(Mod(pgram[, 1, i])^2 / (pgram[, 1, 1] * pgram[, i, i])) + } + + return(coh) + } + +#' phase +#' +#' calculate phase from the spectra and cross-spectra +#' +#' @param pgram \code{numeric array} must be multivariate +#' +#' +#' @return list of phase +#' +#' @export +#' +phase <- function(pgram) { + + dims = dim(pgram) + n_ser <- dims[2] + n_r <- dims[1] + + if (n_ser == 1) { + warning('Cannot calculate phase for a single periodogram') + return(NA) + } + + phase <- matrix(NA_real_, nrow = n_r, ncol = (n_ser - 1)) + + for (i in 2L:(n_ser)) { + + phase[, i-1] <- Arg(pgram[, 1, i]) + + } + + return(phase) +} \ No newline at end of file diff --git a/tests/testthat/test-coherence_phase.R b/tests/testthat/test-coherence_phase.R new file mode 100644 index 0000000..943c71b --- /dev/null +++ b/tests/testthat/test-coherence_phase.R @@ -0,0 +1,9 @@ +test_that("coherence and phase works for single entry", { + + a <- array(rnorm(100), dim = c(100, 1, 1)) + + expect_warning(phase(a)) + expect_warning(coherence(a)) + + +}) From c27361ba60be10a793d9939c659a6129f055b82a Mon Sep 17 00:00:00 2001 From: jkennel Date: Fri, 10 Jan 2020 14:59:55 -0500 Subject: [PATCH 04/15] add multivariate resampling and transfer function --- src/resample_fft.cpp | 421 +++++++++++++++++++++++++++++---- tests/testthat/test-transfer.R | 48 ++++ 2 files changed, 419 insertions(+), 50 deletions(-) create mode 100644 tests/testthat/test-transfer.R diff --git a/src/resample_fft.cpp b/src/resample_fft.cpp index d48fab3..647f0d6 100644 --- a/src/resample_fft.cpp +++ b/src/resample_fft.cpp @@ -118,6 +118,9 @@ List parabolic_weights_rcpp(const int ntap = 1) { return weights_out; } + + + //' @title Resample an fft using varying numbers of sine tapers //' //' @description @@ -463,16 +466,17 @@ List resample_fft_rcpp2( const arma::cx_vec& fftz, -arma::vec pad_data(const arma::vec& psd, +arma::mat pad_data(const arma::mat& psd, const int n, const int n_max, const double eps = 1e-78) { - arma::vec ii(n - 1 + (n_max) * 2); + int nc = psd.n_cols; + arma::mat ii(n - 1 + (n_max) * 2, nc); - ii.head(n_max) = arma::reverse(psd.head(n_max)); - ii.tail(n_max+1) = arma::reverse(psd.tail(n_max+1)); - ii(arma::span(n_max-1, n_max + n-2)) = psd; + ii.head_rows(n_max) = arma::reverse(psd.head_rows(n_max)); + ii.tail_rows(n_max+1) = arma::reverse(psd.tail_rows(n_max+1)); + ii(arma::span(n_max-1, n_max + n-2), arma::span::all) = psd; return(ii + eps); @@ -485,15 +489,17 @@ arma::vec pad_data(const arma::vec& psd, //' //' @return kopt vector // [[Rcpp::export]] -arma::vec riedsid_rcpp(const arma::vec& PSD, const arma::ivec& ntaper){ +arma::vec riedsid_rcpp(const arma::mat& PSD, const arma::ivec& ntaper){ double eps = 1e-78; double sc = 473.3736; double uzero, L, L2, CC; - int nf = PSD.n_elem; + int nf = PSD.n_rows; + int nc = PSD.n_cols; int nt = ntaper.n_elem; int j1, j2; + arma::vec out(nf); arma::ivec taper_vec(nf); // check for single length taper @@ -510,10 +516,10 @@ arma::vec riedsid_rcpp(const arma::vec& PSD, const arma::ivec& ntaper){ } int nadd = 1 + arma::max(nspan); - arma::vec y = arma::log(pad_data(PSD, nf, nadd, eps = 1e-78)); + arma::mat y = arma::log(pad_data(PSD, nf, nadd, eps = 1e-78)); double dy, d2y; - arma::vec yders(nf); + arma::mat yders(nf, nc); for (int j = 0; j < nf; j++) { @@ -527,78 +533,393 @@ arma::vec riedsid_rcpp(const arma::vec& PSD, const arma::ivec& ntaper){ uzero = (L2 - 1) / CC; - // first deriv - dy = as_scalar(u.t() * y(arma::span(j1, j2)) * CC / (L*(L2 - 1))); - // second deriv - d2y = as_scalar((u % u - uzero).t() * - y(arma::span(j1, j2)) * 360 / (L * (L2 - 1) * (L2 - 4))); + for (int i = 0; i < nc; i++) { + + // first deriv + dy = as_scalar(u.t() * y(arma::span(j1, j2), i) * CC / (L*(L2 - 1))); + + // second deriv + d2y = as_scalar((u % u - uzero).t() * + y(arma::span(j1, j2), i) * 360 / (L * (L2 - 1) * (L2 - 4))); + + yders(j, i) = fabs(dy*dy + d2y + eps); + } - yders(j) = fabs(dy*dy + d2y + eps); + out(j) = arma::max(yders.row(j)); } - return(arma::round(pow(sc, 0.2) / arma::pow(yders, 0.4))); + return(arma::round(pow(sc, 0.2) / arma::pow(out, 0.4))); } -/*** R -library(microbenchmark) -library(psd) -n. <- 10000 -set.seed(1234) -x <- cumsum(sample(c(-1, 1), n., TRUE)) -fftz <- fft(x) -taps <- ceiling(runif(n.,10,30)) +// This is for multivariate series -------------------------------------------- -microbenchmark( - rsz1 <- resample_fft_rcpp(fftz, taps, verbose = FALSE), - rsz2 <- resample_fft_rcpp2(fftz, taps, verbose = FALSE), - times = 10 -) -all.equal(rsz1, rsz2) +// For all but short series this should be faster +//' @rdname parabolic_weights_field +//' +//' @param ntap the maximum number of tapers +//' +//' @export +// [[Rcpp::export]] +arma::field parabolic_weights_field(const int ntap) { + + // + // return quadratic spectral weighting factors for a given number of tapers + // Barbour and Parker (2014) Equation 7 + // + arma::field f(ntap, 1); + arma::vec kseq = arma::pow(arma::regspace(0, ntap - 1), 2); + + double t3; + for (int i = 1; i < ntap+1; i++) { + t3 = log(i * (i - 0.25) * (i + 1.0)); + f(i-1,0) = arma::exp(log(1.5) + arma::log(i * i - kseq(arma::span(0, i-1))) - t3); + } + + return(f); + +} -microbenchmark( - kopt1 <- riedsid2(PSD = rsz1$psd, ntaper = taps, constrained = FALSE, verbose = FALSE), - kopt2 <- riedsid2(PSD = rsz1$psd, ntaper = taps, constrained = FALSE, verbose = FALSE, fast = TRUE), - times = 10 -) +//' @title Resample an fft using varying numbers of sine tapers +//' +//' @description +//' Produce an un-normalized psd based on an fft and a vector of optimal sine +//' tapers. +//' +//' @details +//' To produce a psd estimate with our adaptive spectrum estimation method, +//' we need only make one fft calculation initially and then apply the weighting +//' factors given by \code{\link{parabolic_weights_field}}, which this function +//' does. +//' +//' @param fftz complex; a matrix representing the dual-length \code{\link{fft}}; see also the \code{dbl} argument +//' @param tapers integer; a vector of tapers +//' @param verbose logical; should messages be given? +//' @param dbl logical; should the code assume \code{fftz} is dual-length or single-length? +//' @param tapcap integer; the maximum number of tapers which can be applied; note that the length is +//' automatically limited by the length of the series. +//' +//' @return list that includes the auto and cross-spectral density, and the +//' number of tapers +//' +//' @seealso \code{\link{riedsid}} +//' +//' @examples +//' fftz <- complex(real=1:8, imaginary = 1:8) +//' taps <- 1:4 +//' try(resample_mvfft(fftz, taps)) +//' +//' @export +// [[Rcpp::export]] +List resample_mvfft( const arma::cx_mat& fftz, + const arma::ivec& tapers, + bool verbose = true, + const bool dbl = true, + const int tapcap = 10000 ) { + + // + // resample and reweight an fft estimates for a given number of tapers + // using quadratic scaling in Barbour and Parker (2014) + // + // + // needs: + // - fftz: complex vector -- the FFT of the original signal + // - tapers: integer vector -- the number of tapers at each frequency + // + // optional args: + // - verbose: logical -- should warnings be given? + // - dbl: logical -- should the progam assume 'fftz' is for a double-length (padded) series? + // Otherwise a single-length series is assumed. + // - tapcap: integer -- the maximum number of tapers at any frequency + // + + int sc, nf, nc, nt, ne, ne2, nhalf, m2, mleft1, mleft2, Kc, ki, j1, j2; + double wi; + + + arma::cx_double fdc1, fdc2; + + + if (dbl){ + // double-length fft estimates assumed by default + sc = 2; + } else { + // but could be single-length + sc = 1; + } + + // even, double, and half lengths + nf = fftz.n_rows / sc; + nc = fftz.n_cols; + nt = tapers.n_elem; + ne = nf - (nf % 2); + + if (verbose){ + Function msg("message"); + msg(std::string("\tfft resampling")); + } + + if (ne < nf){ + warning("fft was not done on an even length series"); + } + + ne2 = 2 * ne; + nhalf = ne / 2; + + arma::ivec taper_vec(nhalf); + + if (nhalf < 1){ + stop("cannot operate on length-1 series"); + } + + + if (nt == 1){ + warning("forced taper length"); + taper_vec.fill(tapers[0]); + } else { + taper_vec = tapers; + } + + + // set the current number of tapers, limited by a few factors + arma::uvec wh = arma::find(taper_vec > nhalf); + taper_vec(wh).fill(nhalf); + wh = arma::find(taper_vec > tapcap); + taper_vec(wh).fill(tapcap); + wh = arma::find(taper_vec <= 0); + taper_vec(wh).ones(); + + + int mm = taper_vec.max(); + + arma::vec w(mm); + arma::field para = parabolic_weights_field(mm); + + + // + // Calculate the psd by averaging over tapered estimates + // + + + arma::cx_cube psd(nhalf, nc, nc); + psd.fill(arma::cx_double(0.0, 0.0)); + + // each frequency + for (int jj = 0; jj < nhalf; jj++) { + + m2 = 2 * jj; + // number of tapers applied at a given frequency (do not remove m+1 index!) + + Kc = taper_vec[jj]; // orig: m+1, but this was leading to an indexing bug + + // taper sequence and spectral weights + w = para(Kc-1, 0); + + // scan across ki to get vector of + // spectral differences based on modulo indices + for (int ik = 0; ik < Kc; ik++) { + + wi = w[ik]; + ki = ik + 1; + + mleft1 = m2 + ne2 - ki; + mleft2 = m2 + ki; + + j1 = mleft1 % ne2; + j2 = mleft2 % ne2; + + // sum as loop progresses, auto and cross-spectra + + for (int ii = 0; ii < nc; ii++) { + for (int kk = ii; kk < nc; kk++) { + + fdc1 = fftz(j1, ii) - fftz(j2, ii); + fdc2 = fftz(j1, kk) - fftz(j2, kk); + + psd(jj, ii, kk) += (fdc1 * std::conj(fdc2)) * wi; + } + } + } + } + + // Add lower triangle for use in transfer function calculation + for (int ii = 0; ii < nc; ii++) { + for (int kk = 0; kk < nc; kk++) { + if (kk < ii){ + psd(arma::span(), arma::span(ii), arma::span(kk)) = + arma::conj(psd(arma::span(), arma::span(kk), arma::span(ii))); + } + } + } + + // return list to match previous function definition - jrk + return Rcpp::List::create( + Named("freq.inds") = arma::regspace(1, nhalf), + Named("k.capped") = taper_vec(arma::span(0, nhalf-1)), + Named("psd") = psd + ); +} -all.equal(kopt1, kopt2) -microbenchmark( - psd1 <- psdcore(x, verbose = FALSE), - psd2 <- psdcore(x, verbose = FALSE, fast = TRUE), - times = 10 -) +// solve for transfer function -all.equal(psd1, psd2) +//' @title +//' det_vector +//' +//' @description +//' Determinant for an array +//' +//' @param x \code{numeric array} values to evaluate +//' +//' @return vector of determinants +//' +//' +//' @export +// [[Rcpp::export]] +arma::cx_vec det_vector(const arma::cx_cube& x) { + + std::size_t n = x.n_rows; + std::size_t nc = x.n_cols; + + arma::cx_vec output(n); + arma::cx_mat mat(nc, nc); + + for (std::size_t i = 0; i < n; i++){ + mat = x(arma::span(i), arma::span::all, arma::span::all); + output(i) = arma::det(mat); + } + + return(output); + +} +// [[Rcpp::export]] +arma::cx_mat solve_tf(arma::cx_cube x) { + + unsigned int n = x.n_rows; + unsigned int nr = x.n_cols -1; + + arma::cx_mat output(n, nr); + arma::cx_cube numer_base = x(arma::span::all, + arma::span(1, nr), + arma::span(1, nr));; + arma::cx_cube numer(n, nr, nr); + + + arma::cx_cube numer_sol = x(arma::span::all, + arma::span(1, nr), + arma::span(0)); + + arma::cx_vec denom = det_vector(numer_base); + + for (arma::uword i=0; i < nr; i++){ + + // fill new column + numer = numer_base; + numer(arma::span::all, arma::span::all, arma::span(i)) = numer_sol; + + output.col(i) = det_vector(numer) / denom; + + } + + return(output); + +} + + + + +/*** R + +library(microbenchmark) +library(psd) + +n. <- 1000000 +set.seed(1234) +nc <- 2 +x <- matrix(cumsum(sample(c(-1, 1), n., TRUE)), ncol=nc) +fftz <- mvfft(x) +psd <- Re(fftz * Conj(fftz)) +taps <- ceiling(runif(n./nc,10,300)) + microbenchmark( - pspec1 <- pspectrum(x, verbose = FALSE), - pspec2 <- pspectrum(x, verbose = FALSE, fast = TRUE), - times = 10 + a1 <- riedsid_rcpp(as.matrix(psd)[, 1, drop = FALSE], taps), + a2 <- riedsid_rcpp(psd[, 2, drop = FALSE], taps), + b <- riedsid_rcpp(as.matrix(psd), taps), + times = 2 ) -all.equal(pspec1, pspec2) +all.equal(apply(cbind(a1,a2), 1, min), b) microbenchmark( - pilot1 <- pilot_spec(x, verbose = FALSE), - pilot2 <- pilot_spec(x, verbose = FALSE, fast = TRUE), - times = 10 + rsz1 <- resample_fft_rcpp(fftz[,2], taps, verbose = FALSE), + rsz2 <- resample_fft_rcpp2(fftz[,2], taps, verbose = FALSE), + rsz3 <- resample_mvfft(fftz, taps, verbose = FALSE), + times = 2 ) -all.equal(pilot1, pilot2) +all.equal(rsz1, rsz2) +all.equal(rsz1$psd, Re(rsz3$psd[,2,2])) + + +# +# microbenchmark( +# tf <- solve_tf(rsz3$psd), +# times = 5 +# ) + + +# +# microbenchmark( +# kopt1 <- riedsid2(PSD = rsz1$psd, ntaper = taps, constrained = FALSE, verbose = FALSE), +# kopt2 <- riedsid2(PSD = rsz1$psd, ntaper = taps, constrained = FALSE, verbose = FALSE, fast = TRUE), +# times = 10 +# ) + + +# +# all.equal(kopt1, kopt2) +# +# +# microbenchmark( +# psd1 <- psdcore(x, verbose = FALSE), +# psd2 <- psdcore(x, verbose = FALSE, fast = TRUE), +# times = 10 +# ) +# +# +# all.equal(psd1, psd2) +# +# +# +# microbenchmark( +# pspec1 <- pspectrum(x, verbose = FALSE), +# pspec2 <- pspectrum(x, verbose = FALSE, fast = TRUE), +# times = 10 +# ) +# +# all.equal(pspec1, pspec2) +# +# +# microbenchmark( +# pilot1 <- pilot_spec(x, verbose = FALSE), +# pilot2 <- pilot_spec(x, verbose = FALSE, fast = TRUE), +# times = 10 +# ) +# +# all.equal(pilot1, pilot2) */ diff --git a/tests/testthat/test-transfer.R b/tests/testthat/test-transfer.R new file mode 100644 index 0000000..01810eb --- /dev/null +++ b/tests/testthat/test-transfer.R @@ -0,0 +1,48 @@ +test_that("transfer function works", { + + # can we apply and recover a constant response + set.seed(1234) + seq_length <- 10000 + x <- sin(seq(0, 8*pi, length.out = seq_length)) + rnorm(seq_length) + y <- x * 0.5 + xy <- cbind(y, x) + pc <- psdcore(xy, plot = FALSE, verbose = FALSE, ntaper = as.tapers(3)) + + expect_equal(Mod(pc$transfer), + matrix(0.5, nrow = length(pc$transfer), ncol = 1)) + + + # can we apply and recover a variable response + # generate data + set.seed(1234) + kern_length <- 101 + seq_length <- 10000 + x <- sin(seq(0, 8*pi, length.out = seq_length)) + rnorm(seq_length) + kern <- exp(seq(-10, 0, length.out = kern_length)) / kern_length + pad_kern <- c(rep(0, seq_length-kern_length), kern) + # convolution + y <- Re(fft(fft(x) * Conj(fft(pad_kern)), inverse = TRUE)) / seq_length + + # plot(x, type='l') + # points(y, col = 'red', type='l') + + xy <- cbind(y, x) + pc <- psdcore(xy, plot = FALSE, verbose = FALSE, ntaper = as.tapers(3)) + + # convert to impulse response + # this is necessary for the inverse fft as we use half length + pc$transfer <- c(pc$transfer, rev(Conj(pc$transfer))) + + # inverse fft of transfer fun + n <- floor(NROW(pc$transfer) / 2) + 1 + imp <- fft(pc$transfer, inverse = TRUE) / NROW(pc$transfer) + imp <- Mod(imp) * sign(Re(imp)) + mt_impulse <- cumsum(rev(imp)[1:n] + (imp)[1:n]) + + expect_equal(cumsum(rev(kern)), mt_impulse[1:kern_length], tolerance = 1e-3) + + # plot(cumsum(rev(kern)), type='l', col = 'red') + # points(mt_impulse[1:kern_length], cex = 0.5, pch = 20) + + +}) From 9679a9b61a1e214c103529a78e2ba7bba125dd13 Mon Sep 17 00:00:00 2001 From: jkennel Date: Fri, 10 Jan 2020 15:02:19 -0500 Subject: [PATCH 05/15] use multivariate method --- R/func_psdcore.R | 177 ++++++++++++++++++++++------------- tests/testthat/test-spec-I.R | 119 ++++++++++++++++++++++- 2 files changed, 227 insertions(+), 69 deletions(-) diff --git a/R/func_psdcore.R b/R/func_psdcore.R index d62679c..1d63fd7 100644 --- a/R/func_psdcore.R +++ b/R/func_psdcore.R @@ -60,7 +60,15 @@ psdcore <- function(X.d, ...) UseMethod("psdcore") #' @export psdcore.ts <- function(X.d, ...){ frq <- stats::frequency(X.d) - psdcore(as.vector(X.d), X.frq = frq, ...) + psdcore.default(X.d, ...) +} + +#' @rdname psdcore +#' @aliases psdcore.matrix +#' @export +psdcore.matrix <- function(X.d, ...){ + frq <- stats::frequency(X.d) + psdcore(stats::ts(X.d, frequency=frq), ...) } #' @rdname psdcore @@ -77,7 +85,7 @@ psdcore.default <- function(X.d, fast=FALSE, ndecimate, ... - ) { +) { if (!missing(ndecimate)) warning('Support for decimation has been removed.') @@ -98,10 +106,10 @@ psdcore.default <- function(X.d, } X.d <- if (X.frq > 0){ # value represents sampling frequency - stats::ts(X.d, frequency=X.frq) + stats::ts(as.matrix(X.d), frequency=X.frq) } else if (X.frq < 0){ # value is sampling interval - stats::ts(X.d, deltat=abs(X.frq)) + stats::ts(as.matrix(X.d), deltat=abs(X.frq)) } else { stop("bad sampling information") } @@ -127,7 +135,7 @@ psdcore.default <- function(X.d, # initialize fft and other things, since this usually means a first run # # original series - n.o <- psd_envAssignGet(evars[['n.orig']], length(X.d)) + n.o <- psd_envAssignGet(evars[['n.orig']], NROW(X.d)) X <- psd_envAssignGet(evars[['series.orig']], { if (preproc){ # TODO: option for fast-detrend only, assign preproc flag in env (for plotting later) @@ -139,7 +147,7 @@ psdcore.default <- function(X.d, # Force series to be even in length (modulo division) n.e <- psd_envAssignGet(evars[['n.even']], modulo_floor(n.o)) - X.even <- ts(X[seq_len(n.e)], frequency=X.frq) + X.even <- ts(X[seq_len(n.e),], frequency=X.frq) X.even <- psd_envAssignGet(evars[['series.even']], X.even) stopifnot(is.ts(X.even)) @@ -147,8 +155,7 @@ psdcore.default <- function(X.d, nhalf <- psd_envAssignGet(evars[['n.even.half']], n.e/2) # variance of even series - varx <- psd_envAssignGet(evars[['var.even']], drop(stats::var(X.even))) - + varx <- psd_envAssignGet(evars[['var.even']], stats::var(X.even)) # create uniform tapers kseq <- psd_envAssignGet(evars[['last.taper']], { if (len_tapseq == 1){ @@ -160,16 +167,19 @@ psdcore.default <- function(X.d, } }) + ## zero pad and take double-length fft - padded <- as.numeric(c(X.even, rep.int(0, n.e))) + pad <- matrix(0.0, nrow = n.e, ncol = NCOL(X.even)) + padded <- rbind(as.matrix(X.even), pad) + # padded <- as.matrix(X.even) + ## Calculate discrete Fourier tranform # Note fftw is faster for very long series but we are # using stats::fft until fftw is reliably built on CRAN - padded.fft <- psd_envAssignGet(evars[['fft.padded']], stats::fft(padded)) - + padded.fft <- psd_envAssignGet(evars[['fft.padded']], stats::mvfft(padded)) psd_envAssignGet(evars[['fft']], padded.fft) - + } else { if (verbose) warning("Working environment *not* refreshed. Results may be bogus.") @@ -194,7 +204,6 @@ psdcore.default <- function(X.d, ### Select frequencies for PSD evaluation f <- base::seq.int(0, nhalf, by=1) nfreq <- length(f) - ### Calculate the (un-normalized) PSD PSD <- psd_envAssignGet(evars[["last.psdcore"]], { @@ -206,13 +215,10 @@ psdcore.default <- function(X.d, ## resample fft with taper sequence and quadratic weighting # ( this is where the majority of the computational work is ) - kseq <- as.integer(kseq) - if(fast) { - reff <- resample_fft_rcpp2(fftz, kseq, verbose=verbose) - } else { - reff <- resample_fft_rcpp(fftz, kseq, verbose=verbose) - } - + kseq <- as.integer(kseq) + + reff <- resample_mvfft(fftz, tapers = kseq, verbose=verbose) + # return a valid resampled fft or stop if (inherits(reff,'try-error')){ stop("Could not resample fft... inspect with psd_envGet(",evars[['fft']],"), etc.") @@ -238,15 +244,15 @@ psdcore.default <- function(X.d, } }) - # should not be complex at this point! - stopifnot(!is.complex(PSD)) + # should not be complex at this point! - no longer true for cross-spectra -jrk + # stopifnot(!is.complex(PSD)) # there should not be any bad values here! pNAs <- is.na(PSD) if (any(pNAs)) warning("NA psd estimates?!") ## Nyquist frequencies - npsd <- length(PSD) + npsd <- NROW(PSD) frq <- as.numeric(base::seq.int(0, Nyq, length.out=npsd)) ## Update tapers for consistency @@ -258,37 +264,74 @@ psdcore.default <- function(X.d, # ( using the trapezoidal rule, the principal being that the # integrated spectrum should be equal to the variance of the signal ) # - trap.area <- base::sum(PSD, na.rm=TRUE) - mean(PSD[c(1,npsd)], na.rm=TRUE) - area.var.ratio <- varx / trap.area - PSD <- 2 * PSD * area.var.ratio * nhalf + + spec <- matrix(NA_real_, ncol = NCOL(PSD), nrow = NROW(PSD)) + for(i in 1:NCOL(PSD)) { + trap.area <- base::sum(PSD[, i, i, drop = FALSE], na.rm=TRUE) - mean(PSD[c(1,npsd), i, i, drop = FALSE], na.rm=TRUE) + area.var.ratio <- as.matrix(varx)[i,i] / trap.area + spec[, i] <- Re(2 * PSD[, i, i] * area.var.ratio * nhalf) + } ## Assemble final results mtap <- max(kseq, na.rm=TRUE) - PSD.out <- list(freq = as.numeric(frq), - spec = as.numeric(PSD), - coh = NULL, - phase = NULL, - kernel = NULL, - # must be a scalar for plot.spec to give conf ints: - df = 2 * mtap, # 2 DOF per taper, Percival and Walden eqn (370b) - numfreq = npsd, - ## bandwidth - # http://biomet.oxfordjournals.org/content/82/1/201.full.pdf - # half-width W = (K + 1)/{2(N + 1)} - # effective bandwidth ~ 2 W (accurate for many spectral windows) - bandwidth = (mtap + 1) / nhalf, - n.used = psd_envGet(evars[['n.even']]), - orig.n = psd_envGet(evars[['n.orig']]), - series = series, - snames = colnames(X), - method = sprintf("sine multitaper"), - taper = kseq, - pad = TRUE, # always! - detrend = preproc, # always true? - demean = preproc, - timebp = as.numeric(kseq/2), - nyquist.frequency = Nyq - ) + if(NCOL(PSD) > 1) { + + PSD.out <- list(freq = as.numeric(frq), + spec = spec, + pspec = PSD, + transfer = solve_tf(PSD), + coh = coherence(PSD), + phase = phase(PSD), + kernel = NULL, + # must be a scalar for plot.spec to give conf ints: + df = 2 * mtap, # 2 DOF per taper, Percival and Walden eqn (370b) + numfreq = npsd, + ## bandwidth + # http://biomet.oxfordjournals.org/content/82/1/201.full.pdf + # half-width W = (K + 1)/{2(N + 1)} + # effective bandwidth ~ 2 W (accurate for many spectral windows) + bandwidth = (mtap + 1) / nhalf, + n.used = psd_envGet(evars[['n.even']]), + orig.n = psd_envGet(evars[['n.orig']]), + series = series, + snames = colnames(X), + method = sprintf("sine multitaper"), + taper = kseq, + pad = TRUE, # always! + detrend = preproc, # always true? + demean = preproc, + timebp = as.numeric(kseq/2), + nyquist.frequency = Nyq + ) + } else { + PSD.out <- list(freq = as.numeric(frq), + spec = as.numeric(spec), + pspec = NULL, + transfer = NULL, + coh = NULL, + phase = NULL, + kernel = NULL, + # must be a scalar for plot.spec to give conf ints: + df = 2 * mtap, # 2 DOF per taper, Percival and Walden eqn (370b) + numfreq = npsd, + ## bandwidth + # http://biomet.oxfordjournals.org/content/82/1/201.full.pdf + # half-width W = (K + 1)/{2(N + 1)} + # effective bandwidth ~ 2 W (accurate for many spectral windows) + bandwidth = (mtap + 1) / nhalf, + n.used = psd_envGet(evars[['n.even']]), + orig.n = psd_envGet(evars[['n.orig']]), + series = series, + snames = colnames(X), + method = sprintf("sine multitaper"), + taper = kseq, + pad = TRUE, # always! + detrend = preproc, # always true? + demean = preproc, + timebp = as.numeric(kseq/2), + nyquist.frequency = Nyq + ) + } class(PSD.out) <- c("amt","spec") if (plot) pgram_compare(PSD.out, ...) return(invisible(PSD.out)) @@ -400,28 +443,28 @@ pgram_compare.amt <- function(x, f=NULL, X=NULL, log.freq=TRUE, db.spec=TRUE, ta on.exit(par(opar)) #layout(matrix(c(1,2), ncol=1), c(1,2)) - layout(matrix(c(1,1,2,2,3,4), 3, 2, byrow = TRUE)) + layout(matrix(c(1,1,2,2,3,4), nrow=3, ncol=2, byrow = TRUE)) par(mar=c(2, 3.5, 2.3, 1.2), oma=c(2,2,2,1)) par(las=1, tcl = -.3, mgp=c(2.2, 0.4, 0)) par(cex=0.8) ## Periodogram result first - plot(fc., sc., - col="red", type="l", - xaxs="i", xlim=flims, - yaxs="i", ylim=slims, - xlab="", - ylab=slab, - main="") + + matplot(fc., sc., + col="red", type="l", + xaxs="i", xlim=flims, + yaxs="i", ylim=slims, + xlab="", + ylab=slab, + main="") # and adaptive result overlain - lines(f., s., type="l") + matlines(f., s., type="l", col = 'black') legend("bottomleft", c(sprintf("spec.pgram (%i%% cosine taper)", 100*tc.), sprintf("psdcore (%i - %i tapers)", min(t.), max(t.)) ), col=c("red","black"), lty=1, lwd=2, cex=0.9) mtext(flab, side=1, line=2) mtext("Spectral density estimates", adj=0, line=0.2, font=4) - ## Tapers if (is.tapers(t.)){ plot(t., f., xlim=flims, xaxs="i") @@ -430,17 +473,23 @@ pgram_compare.amt <- function(x, f=NULL, X=NULL, log.freq=TRUE, db.spec=TRUE, ta } mtext("Sine-tapers", adj=0, line=0.2, font=4) #mtext(flab, side=1, line=2) - + ## Original series main.pre <- ifelse(detrend | demean, "Modified e", "E") main.post <- sprintf("(dem. %s detr. %s)", detrend, demean) - plot(X, type="l", ylab="units", xlab="", xaxs="i", main="") + matplot(X, type="l", ylab="units", xlab="", xaxs="i", main="", col=c('black', 'blue', 'green', 'orange')) mtext(paste0(main.pre, "ven-length series"), line=1, font=4) mtext(main.post, cex=0.6) mtext(expression("time"), side=1, line=1.7) ## autocorrelation - acf(X, main="") + a <- acf(X, main="", plot = FALSE) + a_plot <- matrix(NA, ncol = NCOL(a$acf), nrow=NROW(a$acf)) + for(i in 1:NCOL(X)) { + a_plot[, i] <- a$acf[, i, i] + } + matplot(a$lag[,1,1], a_plot, type = 'hp', col=c('black', 'blue', 'green', 'orange'), ylab = "ACF") + mtext(expression("lag, time"), side=1, line=1.7) mtext("Auto-correlation function", line=1, font=4) diff --git a/tests/testthat/test-spec-I.R b/tests/testthat/test-spec-I.R index 2454e93..0170b0f 100644 --- a/tests/testthat/test-spec-I.R +++ b/tests/testthat/test-spec-I.R @@ -79,12 +79,12 @@ test_that("pspectrum results are accurate",{ test_that("psdcore arguments are tested",{ set.seed(1234) x <- rnorm(100) - xp1 <- psdcore.default(X.d = x, X.frq = 1, plot = FALSE) - xp2 <- psdcore.default(X.d = x, X.frq = -1, plot = FALSE) + xp1 <- psdcore(X.d = x, X.frq = 1, plot = FALSE) + xp2 <- psdcore(X.d = x, X.frq = -1, plot = FALSE) expect_is(xp1, 'spec') expect_is(xp2, 'spec') expect_equal(xp1,xp2) - expect_error(psdcore.default(X.d = x, X.frq = "1", plot = FALSE)) + expect_error(psdcore(X.d = x, X.frq = "1", plot = FALSE)) }) test_that("psdcore results are accurate",{ @@ -162,7 +162,7 @@ test_that("pilot_spec results are accurate",{ # make sure Nyquist frequencies are correct fn <- max(pc[['freq']]) fn2 <- max(pc2[['freq']]) - + expect_equal(fn, frequency(xt)/2) expect_equal(fn2, frequency(xt2)/2) @@ -205,4 +205,113 @@ test_that("check fast version",{ pilot_spec(xt2, verbose = FALSE, plot = FALSE, fast = TRUE)) }) -## + +test_that("check multivariate autospectra for psdcore",{ + + set.seed(1234) + x1 <- rnorm(100) + x2 <- x1 * 2 + x <- cbind(x1, x2) + ps <- psdcore(x, verbose = FALSE, plot = FALSE, fast = TRUE) + + ps1 <- psdcore(x[,1], verbose = FALSE, plot = FALSE, fast = TRUE) + ps2 <- psdcore(x[,2], verbose = FALSE, plot = FALSE, fast = TRUE) + + expect_equal(ps$spec[,1], ps1$spec) + expect_equal(ps$spec[,2], ps2$spec) + + + ps1 <- psdcore(x[,1], verbose = FALSE, plot = FALSE, fast = FALSE) + ps2 <- psdcore(x[,2], verbose = FALSE, plot = FALSE, fast = FALSE) + + expect_equal(ps$spec[,1], ps1$spec) + expect_equal(ps$spec[,2], ps2$spec) + + # ps <- pspectrum(x) +}) + + + +test_that("check multivariate autospectra for pilot_spec",{ + + set.seed(1234) + x <- matrix(rnorm(200), ncol = 2) + + ps1 <- pilot_spec(x[,1], verbose = FALSE, plot = FALSE, fast = FALSE) + ps2 <- pilot_spec(x[,2], verbose = FALSE, plot = FALSE, fast = FALSE) + + ps <- pilot_spec(x, verbose = FALSE, plot = FALSE, fast = TRUE) + + expect_equal(ps$spec[,1], ps1$spec) + expect_equal(ps$spec[,2], ps2$spec) + + + ps <- pspectrum(x, plot = FALSE) + +}) + + +test_that("check multivariate autospectra, coherence, and phase for psdcore",{ + + # compare spec.pgram to psdcore - the methods aren't equivalent but give + # similar results + + set.seed(1234) + x <- matrix(rnorm(2000), ncol = 2) + + pd <- spec.pgram(x, plot = FALSE, spans = 5) + pd <- normalize(pd, 1, "spectrum", verbose = FALSE) + pc <- psdcore(x, plot = FALSE, verbose = FALSE, ntaper = as.tapers(5)) + + # is there an indexing issue here, need to remove the first value in psdcore + # to get alignment + + pc$coh <- pc$coh[-1] + pc$phase <- pc$phase[-1] + pc$spec <- pc$spec[-1, ] + + n <- length(pd$coh) + pd$coh <- pd$coh[-n] + pd$phase <- pd$phase[-n] + pd$spec <- pd$spec[-n, ] + + + # select narrow range to avoid wrap around issues for phase + sub_range <- 1.2 + wh1 <- which(pd$phase < pi/sub_range & pd$phase > -pi/sub_range) + wh2 <- which(pc$phase < pi/sub_range & pc$phase > -pi/sub_range) + wh <- intersect(wh1, wh2) + + # do regression between methods + phase_coef <- coefficients(lm(pd$phase[wh]~pc$phase[wh])) + coh_coef <- coefficients(lm(pd$coh~pc$coh)) + spec_coef_1 <- coefficients(lm(pd$spec[,1]~pc$spec[,1])) + spec_coef_2 <- coefficients(lm(pd$spec[,2]~pc$spec[,2])) + + + # check intercept + expect_lte(abs(phase_coef[1]), 0.05) + expect_lte(abs(coh_coef[1]), 0.05) + expect_lte(abs(spec_coef_1[1]), 0.05) + expect_lte(abs(spec_coef_2[1]), 0.05) + + + # check slope + expect_lte(abs(phase_coef[2]-1), 0.05) + expect_lte(abs(coh_coef[2]-1), 0.05) + expect_lte(abs(spec_coef_1[2]-1), 0.05) + expect_lte(abs(spec_coef_2[2]-1), 0.05) + + + # plot(pd$spec[,1]) + # points(pc$spec[,1], type='l') + # plot(pd$spec[,2]) + # points(pc$spec[,2], type='l') + # plot(pd$coh) + # points(pc$coh, type='l') + # plot(pd$phase) + # points(pc$phase, type='l') + + +}) + From 535bc1736c1c1652152951332995507a16f747ad Mon Sep 17 00:00:00 2001 From: jkennel Date: Fri, 10 Jan 2020 15:04:12 -0500 Subject: [PATCH 06/15] multivariate riedsid method --- R/func_riedsid.R | 24 +++++++++++++++--------- tests/testthat/test-riedsid.R | 25 +++++++++++++++++++++++++ 2 files changed, 40 insertions(+), 9 deletions(-) diff --git a/R/func_riedsid.R b/R/func_riedsid.R index 2c25b6d..e823741 100644 --- a/R/func_riedsid.R +++ b/R/func_riedsid.R @@ -64,9 +64,9 @@ riedsid.spec <- function(PSD, ...){ Pspec <- PSD[['spec']] Tapseq <- PSD[['freq']] ntaper <- if (is.amt(PSD)){ - PSD[['taper']] + PSD[['taper']] } else { - rep.int(x=1L, times=length(Pspec)) + rep.int(x=1L, times=length(Pspec)) } riedsid(PSD=Pspec, ntaper=ntaper, tapseq=Tapseq, ...) } @@ -124,11 +124,11 @@ riedsid.default <- function(PSD, ntaper = 1L, rss <- if (lsmeth){ # spectral derivatives the preferred way DerivFUN <- function(j, - j1=j-nspan[j]+nadd-1, - j2=j+nspan[j]+nadd-1, - jr=j1:j2, - logY=lY[jr], - dEps=eps){ + j1=j-nspan[j]+nadd-1, + j2=j+nspan[j]+nadd-1, + jr=j1:j2, + logY=lY[jr], + dEps=eps){ u <- jr - (j1 + j2)/2 # rowvec u2 <- u*u # rowvec L <- j2-j1+1 # constant @@ -188,6 +188,9 @@ riedsid.default <- function(PSD, ntaper = 1L, return(kopt) } + + + #' @rdname riedsid #' @export riedsid2 <- function(PSD, ...) UseMethod("riedsid2") @@ -197,20 +200,23 @@ riedsid2 <- function(PSD, ...) UseMethod("riedsid2") riedsid2.spec <- function(PSD, ...){ pspec <- PSD[['spec']] freqs <- PSD[['freq']] + ntap <- if (is.amt(PSD)){ PSD[['taper']] } else { - rep.int(x=1L, times=length(pspec)) + rep.int(x=1L, times=NROW(pspec)) } + riedsid2(PSD=pspec, ntaper=ntap, ...) } + #' @rdname riedsid #' @export riedsid2.default <- function(PSD, ntaper=1L, constrained=TRUE, verbose=TRUE, fast=FALSE, ...){ if (fast) { - kopt <- riedsid_rcpp(PSD, ntaper) + kopt <- riedsid_rcpp(as.matrix(PSD), ntaper) } else { PSD <- as.vector(PSD) diff --git a/tests/testthat/test-riedsid.R b/tests/testthat/test-riedsid.R index 71539d0..64d4342 100644 --- a/tests/testthat/test-riedsid.R +++ b/tests/testthat/test-riedsid.R @@ -38,3 +38,28 @@ test_that("riedsid2 R-version is equal to Rcpp version",{ expect_equal(riedsid2(pa_b, fast=FALSE), riedsid2(pa_b, fast = TRUE)) }) + + +test_that("multivariate riedsid2 works",{ + + set.seed(1234) + x <- matrix(rnorm(200), ncol = 2) + taps <- ceiling(runif(200/2, 10, 300)) + + pd <- stats::spectrum(x, plot=FALSE) + + # each separately and then take the minimum number of tapers + r_s <- cbind(riedsid2(pd$spec[, 1], fast=FALSE), + riedsid2(pd$spec[, 2], fast=FALSE)) + r_s <- apply(r_s, 1, min) + + # multivariate method + r_mv <- riedsid2(pd$spec, fast=TRUE) + expect_equal(r_mv, r_s) + + # spec method works + r_mv_spec <- riedsid2(pd, fast=TRUE) + expect_equal(r_mv_spec, r_s) + + +}) From e3a69bb7bbebd6f95497345d882a6c6193f8521a Mon Sep 17 00:00:00 2001 From: jkennel Date: Fri, 10 Jan 2020 15:07:15 -0500 Subject: [PATCH 07/15] matrix method --- R/func_pilot.R | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/R/func_pilot.R b/R/func_pilot.R index 7d6ef06..d09fae7 100644 --- a/R/func_pilot.R +++ b/R/func_pilot.R @@ -62,10 +62,19 @@ pilot_spec <- function(x, ...) UseMethod("pilot_spec") #' @export pilot_spec.ts <- function(x, ...){ stopifnot(is.ts(x)) + pilot_spec.default(x, ...) +} + + +#' @rdname pilot_spec +#' @aliases pilot_spec.matrix +#' @export +pilot_spec.matrix <- function(x, ...){ frq <- stats::frequency(x) - pilot_spec.default(as.vector(x), x.frequency=frq, ...) + pilot_spec(stats::ts(x, frequency=frq), ...) } + #' @rdname pilot_spec #' @export pilot_spec.default <- function(x, x.frequency=NULL, ntap=NULL, remove.AR=NULL, plot=FALSE, verbose=FALSE, fast = FALSE, ...){ @@ -76,7 +85,7 @@ pilot_spec.default <- function(x, x.frequency=NULL, ntap=NULL, remove.AR=NULL, p stopifnot(length(ntap)==1) stopifnot(length(remove.AR)==1) stopifnot(length(x.frequency)==1) - + # AR spectrum or no? REMAR <- ifelse(remove.AR > 0, TRUE, FALSE) @@ -84,7 +93,6 @@ pilot_spec.default <- function(x, x.frequency=NULL, ntap=NULL, remove.AR=NULL, p if (REMAR) remove.AR <- max(1, min(100, abs(remove.AR))) xprew <- prewhiten(x, x.fsamp=x.frequency, AR.max=remove.AR, detrend=TRUE, impute=TRUE, plot=FALSE, verbose=verbose) - ## Remove and AR model if (REMAR){ # AR fit From f1a10c2e729a992b1aa0c545a170d734f120f208 Mon Sep 17 00:00:00 2001 From: jkennel Date: Fri, 10 Jan 2020 15:08:12 -0500 Subject: [PATCH 08/15] matrix methods --- R/func_whiten.R | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/R/func_whiten.R b/R/func_whiten.R index 5435a2a..1d4809a 100644 --- a/R/func_whiten.R +++ b/R/func_whiten.R @@ -119,7 +119,7 @@ prewhiten.ts <- function(tser, AR.max=0L, detrend=TRUE, demean=TRUE, impute=TRUE stopifnot(is.ts(tser)) sps <- stats::frequency(tser) tstart <- stats::start(tser) - n.o <- length(tser) + n.o <- NROW(tser) ttime <- sps*n.o # NA action on input series @@ -143,26 +143,22 @@ prewhiten.ts <- function(tser, AR.max=0L, detrend=TRUE, demean=TRUE, impute=TRUE AR.max <- abs(AR.max) if (AR.max >= 0){ - # data.frame with fit params - fit.df <- data.frame(xr=base::seq_len(n.o), - xc=base::rep.int(1, n.o), - y=tser) X <- if (detrend){ if (verbose) message("\tdetrending (and demeaning)") - lmdfit <- stats::lm(y ~ xr, fit.df) - as.vector(stats::residuals(lmdfit)) + as.matrix(stats::residuals(stats::lm(tser~base::seq_len(n.o)))) + } else if (demean) { if (verbose) message("\tdemeaning") - lmdfit <- stats::lm(y ~ xc, fit.df) - as.vector(stats::residuals(lmdfit)) + as.matrix(stats::residuals(stats::lm(tser~rep.int(1, n.o)))) + } else { if (verbose) message("\tnothing was done to the timeseries object") tser } # TS object of residuals or equivalent - tser_prew_lm <- stats::ts(X, frequency=sps, start=tstart) + tser_prew_lm <- stats::ts(as.matrix(X), frequency=sps, start=tstart) } @@ -177,7 +173,7 @@ prewhiten.ts <- function(tser, AR.max=0L, detrend=TRUE, demean=TRUE, impute=TRUE ardfit <- stats::ar.yw(tser_prew_lm, aic=TRUE, order.max=AR.max, demean=TRUE) # returns a TS object - tser_prew_ar <- ardfit[['resid']] + tser_prew_ar <- ts(as.matrix(ardfit[['resid']]), frequency=sps, start=tstart) if (impute) tser_prew_ar <- NAFUN(tser_prew_ar) } From 81538e15962d9e41902a0a1f7e975c224c867651 Mon Sep 17 00:00:00 2001 From: jkennel Date: Fri, 10 Jan 2020 15:10:02 -0500 Subject: [PATCH 09/15] matrix method --- R/func_pspectrum.R | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/R/func_pspectrum.R b/R/func_pspectrum.R index 2485e49..3157b17 100644 --- a/R/func_pspectrum.R +++ b/R/func_pspectrum.R @@ -43,9 +43,18 @@ pspectrum <- function(x, ...) UseMethod("pspectrum") #' @export pspectrum.ts <- function(x, ...){ frq <- stats::frequency(x) - psd::pspectrum(as.vector(x), x.frqsamp=frq, ...) + pspectrum.default(x, x.frqsamp=frq, ...) } +#' @rdname pspectrum +#' @aliases pspectrum.matrix +#' @export +pspectrum.matrix <- function(x, ...){ + frq <- stats::frequency(x) + pspectrum(stats::ts(x, frequency=frq), ...) +} + + #' @rdname pspectrum #' @aliases pspectrum.spec #' @export @@ -75,7 +84,8 @@ pspectrum.default <- function(x, verbose=TRUE, no.history=FALSE, plot=FALSE, ...){ - stopifnot(length(x)>1) + stopifnot(NROW(x)>1) + # plotting and iterations if (is.null(niter)) stopifnot(niter>=0) @@ -88,7 +98,7 @@ pspectrum.default <- function(x, # AR switch ordAR <- ifelse(AR, 100, 0) - + for (stage in iter_stages){ if (stage==0){ @@ -96,11 +106,11 @@ pspectrum.default <- function(x, if (niter==0 & plot) plotpsd_ <- TRUE # --- setup the environment --- psd_envRefresh(verbose=verbose) - + # --- pilot spec --- # ** normalization is here: Pspec <- psd::pilot_spec(x, x.frequency=x.frqsamp, ntap=ntap.init, - remove.AR=ordAR, verbose=verbose, plot=plotpsd_, fast=TRUE) + remove.AR=ordAR, verbose=verbose, plot=plotpsd_, fast=TRUE) kopt <- Pspec[['taper']] # ensure series is in the environment @@ -124,7 +134,6 @@ pspectrum.default <- function(x, ## calculate optimal tapers kopt <- riedsid2(Pspec, verbose=rverb, fast = TRUE, ...) - # get data back for plotting, etc. if (stage==niter){ x <- psd_envGet("original_pspectrum_series") @@ -132,7 +141,7 @@ pspectrum.default <- function(x, plotpsd_ <- TRUE } } - + # update spectrum with new tapers # TODO: here's why preproc flags are wrong... Pspec <- psdcore(X.d=x, X.frq=x.frqsamp, ntaper=kopt, From 71ce1915d17df7ff39b0229ace3f7ae435286e4b Mon Sep 17 00:00:00 2001 From: jkennel Date: Fri, 10 Jan 2020 15:18:16 -0500 Subject: [PATCH 10/15] documentation --- NAMESPACE | 8 ++++ R/RcppExports.R | 68 ++++++++++++++++++++++++++++++++++ man/coherence.Rd | 17 +++++++++ man/det_vector.Rd | 17 +++++++++ man/parabolic_weights_field.Rd | 17 +++++++++ man/phase.Rd | 17 +++++++++ man/pilot_spec.Rd | 3 ++ man/psdcore.Rd | 3 ++ man/pspectrum.Rd | 3 ++ man/resample_mvfft.Rd | 44 ++++++++++++++++++++++ src/RcppExports.cpp | 56 +++++++++++++++++++++++++++- src/resample_fft.cpp | 8 +++- 12 files changed, 257 insertions(+), 4 deletions(-) create mode 100644 man/coherence.Rd create mode 100644 man/det_vector.Rd create mode 100644 man/parabolic_weights_field.Rd create mode 100644 man/phase.Rd create mode 100644 man/resample_mvfft.Rd diff --git a/NAMESPACE b/NAMESPACE index 68706b5..cabcd75 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -21,6 +21,7 @@ S3method(parabolic_weights,default) S3method(parabolic_weights,tapers) S3method(pgram_compare,amt) S3method(pilot_spec,default) +S3method(pilot_spec,matrix) S3method(pilot_spec,ts) S3method(plot,tapers) S3method(points,tapers) @@ -29,8 +30,10 @@ S3method(prewhiten,ts) S3method(print,summary.tapers) S3method(print,tapers) S3method(psdcore,default) +S3method(psdcore,matrix) S3method(psdcore,ts) S3method(pspectrum,default) +S3method(pspectrum,matrix) S3method(pspectrum,spec) S3method(pspectrum,ts) S3method(riedsid,default) @@ -52,6 +55,7 @@ S3method(varddiff,spec) export(.spec_confint) export(adapt_message) export(as.tapers) +export(coherence) export(colvec) export(constrain_tapers) export(create_poly) @@ -59,6 +63,7 @@ export(ctap_loess) export(ctap_simple) export(dB) export(data.frame.tapers) +export(det_vector) export(get_adapt_history) export(get_psd_env_name) export(get_psd_env_pointer) @@ -74,9 +79,11 @@ export(new_adapt_history) export(normalize) export(ones) export(parabolic_weights) +export(parabolic_weights_field) export(parabolic_weights_rcpp) export(parabolic_weights_rcpp2) export(pgram_compare) +export(phase) export(pilot_spec) export(prewhiten) export(psd_envAssign) @@ -92,6 +99,7 @@ export(pspectrum_basic) export(rcpp_ctap_simple) export(resample_fft_rcpp) export(resample_fft_rcpp2) +export(resample_mvfft) export(riedsid) export(riedsid2) export(rowvec) diff --git a/R/RcppExports.R b/R/RcppExports.R index cd0d975..59e13af 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -129,3 +129,71 @@ riedsid_rcpp <- function(PSD, ntaper) { .Call('_psd_riedsid_rcpp', PACKAGE = 'psd', PSD, ntaper) } +#' @title Pre-calculation of weight vectors. +#' +#' @description +#' Generate a field of weights for subsetting. +#' +#' @param ntap the maximum number of tapers +#' +#' @return arma::field of weights +#' +#' @export +parabolic_weights_field <- function(ntap) { + .Call('_psd_parabolic_weights_field', PACKAGE = 'psd', ntap) +} + +#' @title Resample an fft using varying numbers of sine tapers +#' +#' @description +#' Produce an un-normalized psd based on an fft and a vector of optimal sine +#' tapers. +#' +#' @details +#' To produce a psd estimate with our adaptive spectrum estimation method, +#' we need only make one fft calculation initially and then apply the weighting +#' factors given by \code{\link{parabolic_weights_field}}, which this function +#' does. +#' +#' @param fftz complex; a matrix representing the dual-length \code{\link{fft}}; see also the \code{dbl} argument +#' @param tapers integer; a vector of tapers +#' @param verbose logical; should messages be given? +#' @param dbl logical; should the code assume \code{fftz} is dual-length or single-length? +#' @param tapcap integer; the maximum number of tapers which can be applied; note that the length is +#' automatically limited by the length of the series. +#' +#' @return list that includes the auto and cross-spectral density, and the +#' number of tapers +#' +#' @seealso \code{\link{riedsid}} +#' +#' @examples +#' fftz <- complex(real=1:8, imaginary = 1:8) +#' taps <- 1:4 +#' try(resample_mvfft(fftz, taps)) +#' +#' @export +resample_mvfft <- function(fftz, tapers, verbose = TRUE, dbl = TRUE, tapcap = 10000L) { + .Call('_psd_resample_mvfft', PACKAGE = 'psd', fftz, tapers, verbose, dbl, tapcap) +} + +#' @title +#' det_vector +#' +#' @description +#' Determinant for an array +#' +#' @param x \code{numeric array} values to evaluate +#' +#' @return vector of determinants +#' +#' +#' @export +det_vector <- function(x) { + .Call('_psd_det_vector', PACKAGE = 'psd', x) +} + +solve_tf <- function(x) { + .Call('_psd_solve_tf', PACKAGE = 'psd', x) +} + diff --git a/man/coherence.Rd b/man/coherence.Rd new file mode 100644 index 0000000..d539bdb --- /dev/null +++ b/man/coherence.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/func_coherence_phase.R +\name{coherence} +\alias{coherence} +\title{coherence} +\usage{ +coherence(pgram) +} +\arguments{ +\item{pgram}{\code{numeric array} must be multivariate} +} +\value{ +list of coherency +} +\description{ +calculate coherency from the spectra and cross-spectra +} diff --git a/man/det_vector.Rd b/man/det_vector.Rd new file mode 100644 index 0000000..f4d444e --- /dev/null +++ b/man/det_vector.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{det_vector} +\alias{det_vector} +\title{det_vector} +\usage{ +det_vector(x) +} +\arguments{ +\item{x}{\code{numeric array} values to evaluate} +} +\value{ +vector of determinants +} +\description{ +Determinant for an array +} diff --git a/man/parabolic_weights_field.Rd b/man/parabolic_weights_field.Rd new file mode 100644 index 0000000..ce39421 --- /dev/null +++ b/man/parabolic_weights_field.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{parabolic_weights_field} +\alias{parabolic_weights_field} +\title{Pre-calculation of weight vectors.} +\usage{ +parabolic_weights_field(ntap) +} +\arguments{ +\item{ntap}{the maximum number of tapers} +} +\value{ +arma::field of weights +} +\description{ +Generate a field of weights for subsetting. +} diff --git a/man/phase.Rd b/man/phase.Rd new file mode 100644 index 0000000..e07b5f0 --- /dev/null +++ b/man/phase.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/func_coherence_phase.R +\name{phase} +\alias{phase} +\title{phase} +\usage{ +phase(pgram) +} +\arguments{ +\item{pgram}{\code{numeric array} must be multivariate} +} +\value{ +list of phase +} +\description{ +calculate phase from the spectra and cross-spectra +} diff --git a/man/pilot_spec.Rd b/man/pilot_spec.Rd index dc49bd2..4d727ec 100644 --- a/man/pilot_spec.Rd +++ b/man/pilot_spec.Rd @@ -5,6 +5,7 @@ \alias{pilot_spectrum} \alias{spec.pilot} \alias{pilot_spec.ts} +\alias{pilot_spec.matrix} \alias{pilot_spec.default} \title{Calculate initial power spectral density estimates} \usage{ @@ -12,6 +13,8 @@ pilot_spec(x, ...) \method{pilot_spec}{ts}(x, ...) +\method{pilot_spec}{matrix}(x, ...) + \method{pilot_spec}{default}(x, x.frequency = NULL, ntap = NULL, remove.AR = NULL, plot = FALSE, verbose = FALSE, fast = FALSE, ...) diff --git a/man/psdcore.Rd b/man/psdcore.Rd index 283a861..ac83d57 100644 --- a/man/psdcore.Rd +++ b/man/psdcore.Rd @@ -3,6 +3,7 @@ \name{psdcore} \alias{psdcore} \alias{psdcore.ts} +\alias{psdcore.matrix} \alias{psdcore.default} \title{Multitaper power spectral density estimates of a series} \usage{ @@ -10,6 +11,8 @@ psdcore(X.d, ...) \method{psdcore}{ts}(X.d, ...) +\method{psdcore}{matrix}(X.d, ...) + \method{psdcore}{default}(X.d, X.frq = NULL, ntaper = as.tapers(5), preproc = TRUE, na.action = stats::na.fail, plot = FALSE, refresh = FALSE, verbose = FALSE, fast = FALSE, ndecimate, ...) diff --git a/man/pspectrum.Rd b/man/pspectrum.Rd index 06ae4ea..101e69b 100644 --- a/man/pspectrum.Rd +++ b/man/pspectrum.Rd @@ -3,6 +3,7 @@ \name{pspectrum} \alias{pspectrum} \alias{pspectrum.ts} +\alias{pspectrum.matrix} \alias{pspectrum.spec} \alias{pspectrum.default} \alias{pspectrum_basic} @@ -13,6 +14,8 @@ pspectrum(x, ...) \method{pspectrum}{ts}(x, ...) +\method{pspectrum}{matrix}(x, ...) + \method{pspectrum}{spec}(x, ...) \method{pspectrum}{default}(x, x.frqsamp = 1, ntap.init = NULL, diff --git a/man/resample_mvfft.Rd b/man/resample_mvfft.Rd new file mode 100644 index 0000000..d06755f --- /dev/null +++ b/man/resample_mvfft.Rd @@ -0,0 +1,44 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{resample_mvfft} +\alias{resample_mvfft} +\title{Resample an fft using varying numbers of sine tapers} +\usage{ +resample_mvfft(fftz, tapers, verbose = TRUE, dbl = TRUE, + tapcap = 10000L) +} +\arguments{ +\item{fftz}{complex; a matrix representing the dual-length \code{\link{fft}}; see also the \code{dbl} argument} + +\item{tapers}{integer; a vector of tapers} + +\item{verbose}{logical; should messages be given?} + +\item{dbl}{logical; should the code assume \code{fftz} is dual-length or single-length?} + +\item{tapcap}{integer; the maximum number of tapers which can be applied; note that the length is +automatically limited by the length of the series.} +} +\value{ +list that includes the auto and cross-spectral density, and the +number of tapers +} +\description{ +Produce an un-normalized psd based on an fft and a vector of optimal sine +tapers. +} +\details{ +To produce a psd estimate with our adaptive spectrum estimation method, +we need only make one fft calculation initially and then apply the weighting +factors given by \code{\link{parabolic_weights_field}}, which this function +does. +} +\examples{ +fftz <- complex(real=1:8, imaginary = 1:8) +taps <- 1:4 +try(resample_mvfft(fftz, taps)) + +} +\seealso{ +\code{\link{riedsid}} +} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 57d6295..d368ee8 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -83,17 +83,65 @@ BEGIN_RCPP END_RCPP } // riedsid_rcpp -arma::vec riedsid_rcpp(const arma::vec& PSD, const arma::ivec& ntaper); +arma::vec riedsid_rcpp(const arma::mat& PSD, const arma::ivec& ntaper); RcppExport SEXP _psd_riedsid_rcpp(SEXP PSDSEXP, SEXP ntaperSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< const arma::vec& >::type PSD(PSDSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type PSD(PSDSEXP); Rcpp::traits::input_parameter< const arma::ivec& >::type ntaper(ntaperSEXP); rcpp_result_gen = Rcpp::wrap(riedsid_rcpp(PSD, ntaper)); return rcpp_result_gen; END_RCPP } +// parabolic_weights_field +arma::field parabolic_weights_field(const int ntap); +RcppExport SEXP _psd_parabolic_weights_field(SEXP ntapSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const int >::type ntap(ntapSEXP); + rcpp_result_gen = Rcpp::wrap(parabolic_weights_field(ntap)); + return rcpp_result_gen; +END_RCPP +} +// resample_mvfft +List resample_mvfft(const arma::cx_mat& fftz, const arma::ivec& tapers, bool verbose, const bool dbl, const int tapcap); +RcppExport SEXP _psd_resample_mvfft(SEXP fftzSEXP, SEXP tapersSEXP, SEXP verboseSEXP, SEXP dblSEXP, SEXP tapcapSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::cx_mat& >::type fftz(fftzSEXP); + Rcpp::traits::input_parameter< const arma::ivec& >::type tapers(tapersSEXP); + Rcpp::traits::input_parameter< bool >::type verbose(verboseSEXP); + Rcpp::traits::input_parameter< const bool >::type dbl(dblSEXP); + Rcpp::traits::input_parameter< const int >::type tapcap(tapcapSEXP); + rcpp_result_gen = Rcpp::wrap(resample_mvfft(fftz, tapers, verbose, dbl, tapcap)); + return rcpp_result_gen; +END_RCPP +} +// det_vector +arma::cx_vec det_vector(const arma::cx_cube& x); +RcppExport SEXP _psd_det_vector(SEXP xSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::cx_cube& >::type x(xSEXP); + rcpp_result_gen = Rcpp::wrap(det_vector(x)); + return rcpp_result_gen; +END_RCPP +} +// solve_tf +arma::cx_mat solve_tf(arma::cx_cube x); +RcppExport SEXP _psd_solve_tf(SEXP xSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::cx_cube >::type x(xSEXP); + rcpp_result_gen = Rcpp::wrap(solve_tf(x)); + return rcpp_result_gen; +END_RCPP +} static const R_CallMethodDef CallEntries[] = { {"_psd_rcpp_ctap_simple", (DL_FUNC) &_psd_rcpp_ctap_simple, 2}, @@ -103,6 +151,10 @@ static const R_CallMethodDef CallEntries[] = { {"_psd_parabolic_weights_rcpp2", (DL_FUNC) &_psd_parabolic_weights_rcpp2, 1}, {"_psd_resample_fft_rcpp2", (DL_FUNC) &_psd_resample_fft_rcpp2, 5}, {"_psd_riedsid_rcpp", (DL_FUNC) &_psd_riedsid_rcpp, 2}, + {"_psd_parabolic_weights_field", (DL_FUNC) &_psd_parabolic_weights_field, 1}, + {"_psd_resample_mvfft", (DL_FUNC) &_psd_resample_mvfft, 5}, + {"_psd_det_vector", (DL_FUNC) &_psd_det_vector, 1}, + {"_psd_solve_tf", (DL_FUNC) &_psd_solve_tf, 1}, {NULL, NULL, 0} }; diff --git a/src/resample_fft.cpp b/src/resample_fft.cpp index 647f0d6..2f076d4 100644 --- a/src/resample_fft.cpp +++ b/src/resample_fft.cpp @@ -560,12 +560,16 @@ arma::vec riedsid_rcpp(const arma::mat& PSD, const arma::ivec& ntaper){ // This is for multivariate series -------------------------------------------- - // For all but short series this should be faster -//' @rdname parabolic_weights_field +//' @title Pre-calculation of weight vectors. +//' +//' @description +//' Generate a field of weights for subsetting. //' //' @param ntap the maximum number of tapers //' +//' @return arma::field of weights +//' //' @export // [[Rcpp::export]] arma::field parabolic_weights_field(const int ntap) { From 300d933b608268c06d81b90a2814ed85f179ad3e Mon Sep 17 00:00:00 2001 From: jkennel Date: Fri, 10 Jan 2020 16:16:57 -0500 Subject: [PATCH 11/15] fix travis stringi? --- .travis.yml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.travis.yml b/.travis.yml index 7378ec5..41a53d4 100644 --- a/.travis.yml +++ b/.travis.yml @@ -19,6 +19,9 @@ r_packages: - fftw - covr +r_github_packages: + - stringi + notifications: email: on_success: change From 6fa8a045b35b7b03f73826cc85b82657501ccb76 Mon Sep 17 00:00:00 2001 From: jkennel Date: Fri, 10 Jan 2020 16:25:30 -0500 Subject: [PATCH 12/15] travis fix try 2 --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 41a53d4..6681e02 100644 --- a/.travis.yml +++ b/.travis.yml @@ -20,7 +20,7 @@ r_packages: - covr r_github_packages: - - stringi + - gagolews/stringi notifications: email: From 90b556fb14071223498426cc77cc840d71206d84 Mon Sep 17 00:00:00 2001 From: jkennel Date: Mon, 13 Jan 2020 10:51:24 -0500 Subject: [PATCH 13/15] resample tests --- NAMESPACE | 1 + R/RcppExports.R | 6 ++- ...unc_coh_phase.R => func_coherence_phase.R} | 0 R/func_psdcore.R | 15 ++++--- man/coherence.Rd | 2 +- man/parabolic_weights_field.Rd | 14 ++++++ man/phase.Rd | 2 +- src/resample_fft.cpp | 6 ++- ...est-coh_phase.R => test-coherence_phase.R} | 0 tests/testthat/test-parabolic_weights.R | 10 +++++ tests/testthat/test-resample.R | 45 +++++++++++++++++++ tests/testthat/test-riedsid.R | 5 ++- 12 files changed, 94 insertions(+), 12 deletions(-) rename R/{func_coh_phase.R => func_coherence_phase.R} (100%) create mode 100644 man/parabolic_weights_field.Rd rename tests/testthat/{test-coh_phase.R => test-coherence_phase.R} (100%) create mode 100644 tests/testthat/test-parabolic_weights.R create mode 100644 tests/testthat/test-resample.R diff --git a/NAMESPACE b/NAMESPACE index cabcd75..e7fecb0 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -102,6 +102,7 @@ export(resample_fft_rcpp2) export(resample_mvfft) export(riedsid) export(riedsid2) +export(riedsid_rcpp) export(rowvec) export(spec_confint) export(spec_details) diff --git a/R/RcppExports.R b/R/RcppExports.R index 9631dd8..a35366a 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -125,15 +125,19 @@ resample_fft_rcpp2 <- function(fftz, tapers, verbose = TRUE, dbl = TRUE, tapcap #' @param ntaper scalar or vector; number of tapers to apply optimization #' #' @return kopt vector +#' @export +#' riedsid_rcpp <- function(PSD, ntaper) { .Call('_psd_riedsid_rcpp', PACKAGE = 'psd', PSD, ntaper) } -#' @rdname parabolic_weights_field +#' @title parabolic_weights_field +#' @rdname parabolic_weights #' #' @param ntap the maximum number of tapers #' #' @export +#' parabolic_weights_field <- function(ntap) { .Call('_psd_parabolic_weights_field', PACKAGE = 'psd', ntap) } diff --git a/R/func_coh_phase.R b/R/func_coherence_phase.R similarity index 100% rename from R/func_coh_phase.R rename to R/func_coherence_phase.R diff --git a/R/func_psdcore.R b/R/func_psdcore.R index 5f353be..2cd12d2 100644 --- a/R/func_psdcore.R +++ b/R/func_psdcore.R @@ -216,13 +216,14 @@ psdcore.default <- function(X.d, ## resample fft with taper sequence and quadratic weighting # ( this is where the majority of the computational work is ) kseq <- as.integer(kseq) - if(is.matrix(fftz)) { - reff <- resample_mvfft(fftz, tapers = kseq, verbose=verbose) - }else if(fast) { - reff <- resample_fft_rcpp2(fftz, kseq, verbose=verbose) - } else { - reff <- resample_fft_rcpp(fftz, kseq, verbose=verbose) - } + + # if(is.matrix(fftz)) { + reff <- resample_mvfft(fftz, tapers = kseq, verbose=verbose) + # } else if(fast) { + # reff <- resample_fft_rcpp2(fftz, kseq, verbose=verbose) + # } else { + # reff <- resample_fft_rcpp(fftz, kseq, verbose=verbose) + # } # return a valid resampled fft or stop if (inherits(reff,'try-error')){ diff --git a/man/coherence.Rd b/man/coherence.Rd index c6476b5..d539bdb 100644 --- a/man/coherence.Rd +++ b/man/coherence.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/func_coh_phase.R +% Please edit documentation in R/func_coherence_phase.R \name{coherence} \alias{coherence} \title{coherence} diff --git a/man/parabolic_weights_field.Rd b/man/parabolic_weights_field.Rd new file mode 100644 index 0000000..c5786df --- /dev/null +++ b/man/parabolic_weights_field.Rd @@ -0,0 +1,14 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{parabolic_weights_field} +\alias{parabolic_weights_field} +\title{parabolic_weights_field} +\usage{ +parabolic_weights_field(ntap) +} +\arguments{ +\item{ntap}{the maximum number of tapers} +} +\description{ +parabolic_weights_field +} diff --git a/man/phase.Rd b/man/phase.Rd index 1c0722f..e07b5f0 100644 --- a/man/phase.Rd +++ b/man/phase.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/func_coh_phase.R +% Please edit documentation in R/func_coherence_phase.R \name{phase} \alias{phase} \title{phase} diff --git a/src/resample_fft.cpp b/src/resample_fft.cpp index 647f0d6..f77bb84 100644 --- a/src/resample_fft.cpp +++ b/src/resample_fft.cpp @@ -488,6 +488,8 @@ arma::mat pad_data(const arma::mat& psd, //' @param ntaper scalar or vector; number of tapers to apply optimization //' //' @return kopt vector +//' @export +//' // [[Rcpp::export]] arma::vec riedsid_rcpp(const arma::mat& PSD, const arma::ivec& ntaper){ @@ -562,11 +564,13 @@ arma::vec riedsid_rcpp(const arma::mat& PSD, const arma::ivec& ntaper){ // For all but short series this should be faster -//' @rdname parabolic_weights_field +//' @title parabolic_weights_field +//' @rdname parabolic_weights //' //' @param ntap the maximum number of tapers //' //' @export +//' // [[Rcpp::export]] arma::field parabolic_weights_field(const int ntap) { diff --git a/tests/testthat/test-coh_phase.R b/tests/testthat/test-coherence_phase.R similarity index 100% rename from tests/testthat/test-coh_phase.R rename to tests/testthat/test-coherence_phase.R diff --git a/tests/testthat/test-parabolic_weights.R b/tests/testthat/test-parabolic_weights.R new file mode 100644 index 0000000..8117670 --- /dev/null +++ b/tests/testthat/test-parabolic_weights.R @@ -0,0 +1,10 @@ +test_that("parabolic weight methods give equivalent results", { + + w1 <- parabolic_weights_rcpp(ntap = 5)$taper_weights + w2 <- as.numeric(parabolic_weights_rcpp2(ntap = 5)) + w3 <- as.numeric(parabolic_weights_field(ntap = 5)[[5]]) + + expect_equal(w1, w2) + expect_equal(w2, w3) + +}) diff --git a/tests/testthat/test-resample.R b/tests/testthat/test-resample.R new file mode 100644 index 0000000..aabee62 --- /dev/null +++ b/tests/testthat/test-resample.R @@ -0,0 +1,45 @@ +test_that("resampling methods give same results", { + + n. <- 100000 + set.seed(1234) + nc <- 2 + x <- matrix(cumsum(sample(c(-1, 1), n., TRUE)), ncol=nc) + fftz <- mvfft(x) + psd <- Re(fftz * Conj(fftz)) + taps <- ceiling(runif(n./nc,10,300)) + + rsz1a <- resample_fft_rcpp(fftz[,1], taps, verbose = FALSE) + rsz1b <- resample_fft_rcpp2(fftz[,1], taps, verbose = FALSE) + + rsz2a <- resample_fft_rcpp(fftz[,2], taps, verbose = FALSE) + rsz2b <- resample_fft_rcpp2(fftz[,2], taps, verbose = FALSE) + + rsz3 <- resample_mvfft(fftz, taps, verbose = FALSE) + + + expect_equal(as.numeric(rsz1a$psd), as.numeric(Re(rsz3$psd[,1,1]))) + expect_equal(as.numeric(rsz2a$psd), as.numeric(Re(rsz3$psd[,2,2]))) + + expect_equal(as.numeric(rsz1b$psd), as.numeric(Re(rsz3$psd[,1,1]))) + expect_equal(as.numeric(rsz2b$psd), as.numeric(Re(rsz3$psd[,2,2]))) + +}) + +test_that("riedsid_rcpp gives minimum of matrix", { + + n. <- 100000 + set.seed(1234) + nc <- 2 + x <- matrix(cumsum(sample(c(-1, 1), n., TRUE)), ncol=nc) + fftz <- mvfft(x) + psd <- Re(fftz * Conj(fftz)) + taps <- ceiling(runif(n./nc,10,300)) + + + a1 <- riedsid_rcpp(as.matrix(psd)[, 1, drop = FALSE], taps) + a2 <- riedsid_rcpp(psd[, 2, drop = FALSE], taps) + b <- riedsid_rcpp(as.matrix(psd), taps) + + expect_equal(apply(cbind(a1,a2), 1, min), as.numeric(b)) + +}) diff --git a/tests/testthat/test-riedsid.R b/tests/testthat/test-riedsid.R index c0d0ad7..4abbba4 100644 --- a/tests/testthat/test-riedsid.R +++ b/tests/testthat/test-riedsid.R @@ -62,5 +62,8 @@ test_that("multivariate riedsid2 works",{ expect_equal(r_mv_spec, r_s) -}) +}) + + + \ No newline at end of file From 85590d068aed09fad1d75c97940c488ad6ff5cc5 Mon Sep 17 00:00:00 2001 From: jkennel Date: Mon, 13 Jan 2020 11:52:21 -0500 Subject: [PATCH 14/15] increase code coverage --- tests/testthat/test-normalize.R | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) create mode 100644 tests/testthat/test-normalize.R diff --git a/tests/testthat/test-normalize.R b/tests/testthat/test-normalize.R new file mode 100644 index 0000000..bd90834 --- /dev/null +++ b/tests/testthat/test-normalize.R @@ -0,0 +1,22 @@ +test_that("normalize works", { + + set.seed(1234) + x <- rnorm(2000) + + pd <- spec.pgram(x, plot = FALSE, spans = 5) + + # spec test + pd1a <- normalize(pd, 1, "spectrum", verbose = FALSE) + pd1b <- normalize(pd, -1, "spectrum", verbose = FALSE) + + expect_error(normalize(pd, Fsamp = 0, src = 'spectrum', verbose = FALSE)) + expect_equal(pd1a, pd1b) + + # list test + class(pd) <- NULL + pd2 <- normalize(pd, Fsamp = 1, src = "spectrum", verbose = FALSE) + class(pd2) <- 'spec' + expect_equal(pd1a, pd2) + + +}) From 45e452af7a6ca35f5d0bc5c0b49712824bf76ec9 Mon Sep 17 00:00:00 2001 From: jkennel Date: Mon, 13 Jan 2020 13:07:10 -0500 Subject: [PATCH 15/15] increase code coverage --- tests/testthat/test-resample.R | 99 +++++++++++++++++++++++++++++++++- 1 file changed, 98 insertions(+), 1 deletion(-) diff --git a/tests/testthat/test-resample.R b/tests/testthat/test-resample.R index aabee62..0ee147a 100644 --- a/tests/testthat/test-resample.R +++ b/tests/testthat/test-resample.R @@ -16,7 +16,7 @@ test_that("resampling methods give same results", { rsz3 <- resample_mvfft(fftz, taps, verbose = FALSE) - + expect_equal(as.numeric(rsz1a$psd), as.numeric(Re(rsz3$psd[,1,1]))) expect_equal(as.numeric(rsz2a$psd), as.numeric(Re(rsz3$psd[,2,2]))) @@ -43,3 +43,100 @@ test_that("riedsid_rcpp gives minimum of matrix", { expect_equal(apply(cbind(a1,a2), 1, min), as.numeric(b)) }) + + +test_that("check verbose gives message", { + + n. <- 100000 + set.seed(1234) + nc <- 2 + x <- matrix(cumsum(sample(c(-1, 1), n., TRUE)), ncol=nc) + fftz <- mvfft(x) + psd <- Re(fftz * Conj(fftz)) + taps <- ceiling(runif(n./nc,10,300)) + + expect_message(resample_fft_rcpp(fftz[,1], taps, verbose = TRUE)) + expect_message(resample_fft_rcpp2(fftz[,1], taps, verbose = TRUE)) + expect_message(resample_mvfft(fftz, taps, verbose = TRUE)) + +}) + + +test_that("check forced taper length", { + + n. <- 10000 + set.seed(1234) + nc <- 2 + x <- matrix(cumsum(sample(c(-1, 1), n., TRUE)), ncol=nc) + fftz <- mvfft(x) + psd <- Re(fftz * Conj(fftz)) + taps <- ceiling(runif(n./nc,10,300)) + + expect_warning(resample_fft_rcpp(fftz[,1], 3, verbose = FALSE)) + expect_warning(resample_fft_rcpp2(fftz[,1], 3, verbose = FALSE)) + expect_warning(resample_mvfft(fftz, 3, verbose = FALSE)) + + expect_equal(unique(resample_fft_rcpp(fftz, 3, verbose = FALSE)$k.capped), 3) + expect_equal(unique(resample_fft_rcpp2(fftz, 3, verbose = FALSE)$k.capped), 3) + expect_equal(unique(resample_mvfft(fftz, 3, verbose = FALSE)$k.capped), 3) + +}) + + +test_that("test odd length fft", { + + n. <- 204 + set.seed(1234) + nc <- 2 + x <- matrix(cumsum(sample(c(-1, 1), n., TRUE)), ncol=nc) + fftz <- mvfft(x) + psd <- Re(fftz * Conj(fftz)) + taps <- ceiling(runif(n./nc,10,300)) + + expect_warning(resample_fft_rcpp(fftz[,1], taps, verbose = FALSE)) + expect_warning(resample_fft_rcpp2(fftz[,1], taps, verbose = FALSE)) + expect_warning(resample_mvfft(fftz, taps, verbose = FALSE)) + +}) + +test_that("short series gives error", { + + n. <- 2 + set.seed(1234) + nc <- 1 + x <- matrix(cumsum(sample(c(-1, 1), n., TRUE)), ncol=nc) + fftz <- mvfft(x) + psd <- Re(fftz * Conj(fftz)) + taps <- ceiling(runif(n./nc,10,300)) + + expect_error(resample_fft_rcpp(fftz[,1], taps, verbose = FALSE)) + expect_error(resample_fft_rcpp2(fftz[,1], taps, verbose = FALSE)) + expect_error(resample_mvfft(fftz, taps, verbose = FALSE)) + +}) + + + +# single length does not appear to be correctly implemented +# test_that("test single length", { +# +# n. <- 40 +# set.seed(1234) +# nc <- 1 +# x <- matrix(cumsum(sample(c(-1, 1), n., TRUE)), ncol=nc) +# fftz <- mvfft(x) +# psd <- Re(fftz * Conj(fftz)) +# taps <- ceiling(runif(n./nc,10,50)) +# +# a1 <- resample_fft_rcpp(fftz[,1], taps, verbose = FALSE, dbl = 0) +# a2 <- resample_fft_rcpp(fftz[,1], taps, verbose = FALSE, dbl = 1) +# b1 <- resample_fft_rcpp2(fftz[,1], taps, verbose = FALSE, dbl = 0) +# b2 <- resample_fft_rcpp2(fftz[,1], taps, verbose = FALSE, dbl = 1) +# c1 <- resample_mvfft(fftz, taps, verbose = FALSE, dbl = 0) +# c2 <- resample_mvfft(fftz, taps, verbose = FALSE, dbl = 1) +# +# expect_equal(a2, b2) +# expect_equal(a1, b1) +# expect_equal(b2, c2) +# +# })