R/func_riedsid.R

In psd: Adaptive, Sine-Multitaper Power Spectral Density and Cross Spectrum Estimation

#' Constrained, optimal tapers using the Riedel & Sidorenko--Parker method
#' 
#' Estimates the optimal number of tapers at each frequency of
#' given PSD, using a modified Riedel-Sidorenko MSE recipe (RS-RLP).
#' 
#' @details
#' The optimization is as follows. First, weighted derivatives of the 
#' input PSD are computed.
#' Using those derivatives the optimal number of tapers is found through the 
#' RS-RLP formulation.
#' Constraints are then placed on the practicable number of tapers.
#' 
#' \code{\link{riedsid2}} is a new (faster) implementation which does not allow 
#' for multiple constraint methods; this is the preferred function to use.
#'
#' \subsection{Taper constraints}{
#' The parameter \code{c.method} provides an option to change the method
#' of taper constraints.  A description of each may be found in 
#' the documentation for \code{\link{constrain_tapers}}.
#'
#' Once can use \code{constrained=FALSE} to turn off all taper constraints; this
#' could lead to strange behavior though.
#' }
#'
#' \subsection{Spectral derivatives}{
#' The parameter \code{Deriv.method} determines which method is used
#' to estimate derivatives.
#' \itemize{
#' \item{\code{"local_qls"}}{ (\strong{default}) uses quadratic weighting and
#' local least-squares estimation; this can be slower than \code{"spg"}.}
#' \item{\code{"spg"}}{ uses \code{\link{splineGrad}}; then, additional arguments
#' may be passed to control the smoothness of the derivatives
#' (e.g \code{spar} in \code{\link{smooth.spline}}).}
#' }
#' }
#'
#' @section Warning:
#' The \code{"spg"} can become numerically unstable, and it's not clear when it will
#' be the preferred over the \code{"local_qls"} method, other than for efficiency's sake.
#'
#' @export
#' @author A.J. Barbour adapted original by R.L. Parker
#' 
#' @param PSD vector or class \code{'amt'} or \code{'spec'}; the spectral values used to optimize taper numbers
#' @param ntaper scalar or vector; number of tapers to apply optimization
#' @param tapseq vector; representing positions or frequencies (same length as \code{PSD})
#' @param Deriv.method character; choice of gradient estimation method 
#' @param constrained logical; apply constraints with \code{\link{constrain_tapers}}; \code{FALSE} turns off constraints
#' @param c.method string; constraint method to use with \code{\link{constrain_tapers}}, only if \code{constrained=TRUE}
#' @param verbose logical; should messages be printed?
#' @param fast logical; use faster method?
#' @param riedsid_column scalar integer; which column to use in multivariate optimization. If the value is 0 the maximum number of tapers for all columns is chosen. If the value is < 0 the minimum number of tapers for all columns is chosen. If the value is 1, 2, 3, etc. the number of tapers is based on the column selected.
#' @param ... optional arguments passed to \code{\link{constrain_tapers}}
#' @return Object with class \code{'tapers'}
#' 
#' @seealso \code{\link{constrain_tapers}}, \code{\link{resample_fft_rcpp}}, \code{\link{psdcore}}, \code{\link{pspectrum}}
#' @example inst/Examples/rdex_riedsid.R
riedsid <- function(PSD, ...){
  UseMethod("riedsid")
}

#' @rdname riedsid
#' @export
riedsid.spec <- function(PSD, ...){
  Pspec <- PSD[['spec']]
  Tapseq <- PSD[['freq']]
  ntaper <- if (is.amt(PSD)){
    PSD[['taper']]
  } else {
    rep.int(x=1L, times=length(Pspec))
  }
  riedsid(PSD=Pspec, ntaper=ntaper, tapseq=Tapseq, ...)
}


#' @rdname riedsid
#' @export
riedsid.default <- function(PSD, ntaper = 1L, 
                            tapseq=NULL, 
                            Deriv.method=c("local_qls","spg"),
                            constrained=TRUE,
                            c.method=NULL,
                            verbose=TRUE, ...) {
  .Deprecated('riedsid2', package='psd', old='riedsid')
  ## spectral values
  PSD <- as.vector(PSD)
  # num freqs
  nf <- psd_envAssignGet("num_freqs", length(PSD))
  # prelims
  #  A small number to protect against zeros
  eps <- 1e-78 
  # vectorize initial estimate
  Zeros <- zeros(nf)
  nt <- length(ntaper)
  ntap <- if (nt==1){
    ntaper + Zeros
  } else {
    ntaper
  }
  if (!is.tapers(ntap)) ntap <- as.tapers(ntap)
  
  # Set the number of tapers to within the range: 1/2 nf, 7/5 ntap
  # rowMins produces a rowvec of rowwise minimums; convert to colvec
  nspan <- minspan(ntap, nf)
  
  # The spectral gradients should be in log-space, so
  # create a log spec, and pad to handle begnning and end values
  nadd <- 1 + max(nspan)
  Y <- c(PSD[nadd:2], PSD, PSD[(nf-1):(nf-nadd)])
  Y[Y <= 0] <- eps
  lY <- log(Y)
  dY <- d2Y <- Zeros
  #
  kseq <- if (is.null(tapseq) | (length(tapseq) != nf)){
    seq.int(from=0, to=1/2, length.out=length(PSD))
  } else {
    tapseq
  }
  
  #
  # Smooth spectral derivatives
  #
  lsmeth <- switch(match.arg(Deriv.method), local_qls=TRUE, spg=FALSE)
  stopifnot(exists("lsmeth"))
  
  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){
      u <- jr - (j1 + j2)/2 # rowvec 
      u2 <- u*u             # rowvec
      L <- j2-j1+1          # constant
      L2 <- L*L             # constant
      LL2 <- L*L2           # constant
      LL2L <- LL2 - L       # constant
      #
      CC <- 12
      uzero <- (L2 - 1)/CC  # constant
      #
      # first deriv
      dY <- u %*% logY * CC / LL2L
      # second deriv
      d2Y <- (u2 - uzero) %*% logY * 360 / LL2L / (L2 - 4)
      return(c(fdY2=dY*dY, fd2Y=d2Y, fdEps=dEps))
    }
    DX <- seq_len(nf)
    ##
    ## Bottleneck:
    RSS <- vapply(X=DX, FUN=DerivFUN, FUN.VALUE=c(1,1,1))
    ##
    attr(RSS, which="lsderiv") <- lsmeth
    RSS <- psd_envAssignGet("spectral_derivatives.ls", RSS)
    msg <- "local quadratic regression"
    # sums:
    #[ ,1] fdY2
    #[ ,2] fd2Y
    #[ ,3] fdEps
    abs(colSums(RSS))
  } else {
    RSS <- splineGrad(dseq=log(0.5+kseq), #seq.int(0,.5,length.out=length(PSD))), 
                      dsig=log(PSD),
                      plot.derivs=FALSE, ...) #, spar=1)
    attr(RSS, which="lsderiv") <- lsmeth
    RSS <- psd_envAssignGet("spectral_derivatives", RSS) 
    #returns log
    RSS[,2:4] <- exp(RSS[,2:4])
    msg <- "weighted cubic spline"
    abs(eps + RSS[,4] + RSS[,3]**2)
  }
  if (verbose) message(sprintf("Using spectral derivatives from  %s", msg))
  #
  #(480)^0.2*abs(PSD/d2psd)^0.4
  # Original form:  kopt = 3.428*abs(PSD ./ d2psd).^0.4;
  # kopt = round( 3.428 ./ abs(eps + d2Y + dY.^2).^0.4 );
  #
  sc <- ifelse(TRUE, 473.3736, 480)
  kopt <- (sc ** 0.2)/(rss ** 0.4)
  #
  # Constrain tapers
  kopt <- if (constrained){
    constrain_tapers(tapvec = kopt, tapseq=kseq, constraint.method=c.method, verbose=verbose, ...)
  } else {
    as.tapers(kopt)
  }
  #
  return(kopt)
} 




#' @rdname riedsid
#' @export
riedsid2 <- function(PSD, ...) UseMethod("riedsid2")

#' @rdname riedsid
#' @export
riedsid2.spec <- function(PSD, ...){
  pspec <- PSD[['spec']]
  freqs <- PSD[['freq']]
  
  ntap <- if (is.amt(PSD)){
    PSD[['taper']]
  } else {
    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, 
                             riedsid_column = 0L, 
                             ...){
  
  if (fast) {
    kopt <- riedsid_rcpp(as.matrix(PSD), ntaper, riedsid_column)
  } else {
    
    PSD <- as.vector(PSD)
    ntaper <- as.vector(ntaper)
    
    #  A small number to protect against zeros
    eps <- 1e-78
    nf <- length(PSD)
    nt <- length(ntaper)
    if (nt == 1) ntaper <- rep.int(x=ntaper, times=nf)
    # some constraints
    nspan <- ceiling( pmin.int( nf/2, 7*ntaper/5 ) )
    nadd <- 1 + max(nspan)
    # Create log psd, and pad to handle beginning and end values
    ist <- nadd:2
    iend <- seq(from=(nf - 1), to=(nf - nadd))
    S <- as.numeric(c(PSD[ist], PSD, PSD[iend])) + eps
    Y <- log(S)
    DerivFUN <- function(j){
      j1 <- j - nspan[j] + nadd - 1
      j2 <- j + nspan[j] + nadd - 1
      jseq <- j1:j2
      u <- jseq - (j1 + j2)/2
      L <- j2 - j1 + 1
      CC <- 12
      #
      uzero <- (L^2 - 1)/CC
      #
      # first deriv
      dY <- u  %*%  Y[jseq] * CC / (L*(L*L - 1))
      # second deriv
      d2Y <- (u*u - uzero)  %*%  Y[jseq] * 360 / (L*(L^2 - 1)*(L^2 - 4))
      #
      return(c(eps=eps, d2Y=d2Y, dYsq=dY*dY))
    }
    # Calculate derivatives
    yders <- vapply(X=seq_len(nf), FUN=DerivFUN, FUN.VALUE=c(1,1,1))
    # and optimal tapers
    sc <- ifelse(TRUE, 473.3736, 480)
    kopt <- round( sc**0.2 / abs(colSums(yders))**0.4 )
  }
  
  kopt <- if (constrained){
    constrain_tapers(tapvec = kopt, verbose = verbose, ...)
  } else {
    as.tapers(kopt)
  }
  
  return(kopt)
}