From d827151c4ca87464151a06c58c3aa284e23c61b6 Mon Sep 17 00:00:00 2001 From: mcrucifix <michel.crucifix@uclouvain.be> Date: Thu, 26 Oct 2023 15:44:26 +0200 Subject: [PATCH] cleanup and check --- DESCRIPTION | 4 +- NAMESPACE | 11 +++ NAMESPACE.old | 10 --- R/cwt_morlet.R | 10 ++- R/mem.R | 31 ++++++-- R/mfft.R | 7 +- R/phase_randomize.R | 3 +- R/prettylab.R | 183 +++++++++++++++++++++++--------------------- R/ssa.R | 4 + README.md | 6 -- cleanup | 3 + inst/REFERENCES.bib | 45 +++++++++++ man/axTexpr.Rd | 67 ---------------- man/mem.Rd | 22 +++++- man/pretty10exp.Rd | 32 -------- src/fmft.c | 9 +++ src/nrutil.c | 21 ++++- 17 files changed, 241 insertions(+), 227 deletions(-) delete mode 100644 NAMESPACE.old create mode 100755 cleanup delete mode 100644 man/axTexpr.Rd delete mode 100644 man/pretty10exp.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 882f36c..2e7d650 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: gtseries Type: Package -Title: Time Series analysis for galcial cycles +Title: Time Series analysis for paleoclimate and cyclostrat purposes Version: 1.04 Date: 2022-05-28 Author: Michel Crucifix [aut, cre] (<https://orcid.org/0000-0002-3437-4911>) @@ -9,5 +9,5 @@ Description: More about what it does (maybe more than one line) License: What license is it under? LazyLoad: yes RoxygenNote: 7.2.3 -Imports: Rdpack +Imports: Rdpack, sfsmisc RdMacros: Rdpack diff --git a/NAMESPACE b/NAMESPACE index 3d31d16..c738dde 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,6 +1,7 @@ # Generated by roxygen2: do not edit by hand S3method(plot,SSAObject) +S3method(plot,memObject) S3method(plot,wavelet) export(approx_ts) export(arspec) @@ -10,4 +11,14 @@ export(mfft) export(powerspectrum.wavelet) export(ssa) importFrom(Rdpack,reprompt) +importFrom(graphics,mtext) +importFrom(sfsmisc,axTexpr) +importFrom(sfsmisc,eaxis) +importFrom(stats,ar) +importFrom(stats,as.ts) +importFrom(stats,coef) +importFrom(stats,fft) +importFrom(stats,runif) +importFrom(stats,start) +importFrom(stats,time) useDynLib(gtseries) diff --git a/NAMESPACE.old b/NAMESPACE.old deleted file mode 100644 index e300e16..0000000 --- a/NAMESPACE.old +++ /dev/null @@ -1,10 +0,0 @@ -export(cwt_morlet) -export(mfft) -export(approx_ts) -export(mem) -export(phase_randomize) -export(ssa) -useDynLib(gtseries) -S3method(plot,wavelet) -S3method(plot,memObject) -importFrom(grDevices, colorRampPalette) diff --git a/R/cwt_morlet.R b/R/cwt_morlet.R index e6c871b..ea5d3de 100644 --- a/R/cwt_morlet.R +++ b/R/cwt_morlet.R @@ -120,11 +120,13 @@ cross_morlet <- function(A, B, ...) #' Continous Morlet Wavelet Transform +#' @importFrom stats fft +#' @importFrom stats time #' @export cwt_morlet cwt_morlet <- function (A,inter=20,k0=5.6,amin=1,amax=Inf,calcmask=TRUE,scale=NA,deriv=FALSE) { y <- A - xx <- time(A) + xx <- stats::time(A) deltat<-deltat(A) ny<-length(y); @@ -160,7 +162,7 @@ cwt_morlet <- function (A,inter=20,k0=5.6,amin=1,amax=Inf,calcmask=TRUE,scale=NA k<-(0:(n-1))*2*pi/n - f<-fft(x); + f<-stats::fft(x); J<-length(scale); @@ -175,8 +177,8 @@ cwt_morlet <- function (A,inter=20,k0=5.6,amin=1,amax=Inf,calcmask=TRUE,scale=NA ## norm=sqrt(scale[a1]*k[2])*(pi^(-0.25))/sqrt(nn); norm=2; daughter=norm*exp(expnt); - wave[,a1]=fft(f*daughter ,inverse=TRUE)/n - if (deriv) dwave[,a1]=fft(f*daughter*(-(1i*k)) ,inverse=TRUE)/n; + wave[,a1]=stats::fft(f*daughter ,inverse=TRUE)/n + if (deriv) dwave[,a1]=stats::fft(f*daughter*(-(1i*k)) ,inverse=TRUE)/n; if (calcmask) { mask[1:min(n,ceiling(sqrt(2)*scale[a1])),a1] = NA; mask[(n-min(n,ceiling(sqrt(2)*scale[a1]))):n,a1] = NA; diff --git a/R/mem.R b/R/mem.R index 36fb6b4..03dad80 100644 --- a/R/mem.R +++ b/R/mem.R @@ -1,15 +1,28 @@ - #' Maximum entropy method - -#' @export mem -mem <- function (A,method='burg',order.max=90,...) -{ +#' +#' Maximum entropy spectral analysis (MESA) is the statistical estimation of the power spectrum of a stationary time series using the maximum entropy (ME) method. The resulting ME spectral estimator is parametric and equivalent to an autoregressive (AR) spectral estimator (copied from \\insertCite{pardo-iguzquiza21aa}{gtseries}). The method relies heavily on the standard-R `ar` function, of which it is essentially a wrapper associated with a plot class function (authors: Martyn Plummer. C code for `ar.burg` by B. D. Ripley.) +#' +#' @param A input, a time-series object +#' @param method Estimation method, passed to the standard `ar` function, which is one of "yule-walker", "burg", "ols", "mle", "yw". Defaults to "burg". +#' @param order.max maximum order of the autoregressive model being fitted +#' @return a memObject +#' @references +#' \insertRef{Burg75aa}{gtseries} +#' \insertRef{percival98}{gtseries} +#' \insertRef{Ghil02aa}{gtseries} +#' \insertRef{pardoiguzquiza06aa}{gtseries} +#' \insertRef{pardo-iguzquiza21aa}{gtseries} +#' @importFrom graphics mtext +#' @importFrom stats ar +#' @importFrom stats coef +#' @export mem +mem <- function (A,method='burg',order.max=90,...) { deltat<-deltat(A) minf <- 2./length(A) f <- 10^(seq(log10(minf),log10(0.5),length.out=400)) - U <- ar(A,method=method,order.max=order.max,...) + U <- stats::ar(A,method=method,order.max=order.max,...) ## v.maice multiplied by deltat to have power density (I still need to really undersant ## why, but this works) @@ -24,6 +37,8 @@ mem <- function (A,method='burg',order.max=90,...) } +#' @rdname mem +#' @export plot.memObject <- function (memObject,period=FALSE,xaxp=NULL,yaxp=NULL,...) { @@ -56,7 +71,7 @@ axis.log10 <- function(side=1,title="Power",factor=1,line=3,outer.at=TRUE,small. small_ticks <- small_ticks[which(small_ticks > range[1] & small_ticks < range[2])] axis(side,at=small_ticks^factor,labels=FALSE,tcl=3/5*par("tcl"),...) } - mtext(title,side=side,line=3,cex=par("cex.lab")*par("cex")) + graphics::mtext(title,side=side,line=3,cex=par("cex.lab")*par("cex")) } add.intercept <- function (memObject, xlim = range(memObject$frequency)[2]*c(5.e-2,0.5),annote=TRUE) @@ -67,7 +82,7 @@ print(range(memObject$frequency)[2]) t <- which (frequency > xlim[1] & frequency < xlim[2]) fm1 <- lm(log10(power) ~ log10(frequency),as.data.frame(memObject[t,])) -lines(xlim,10^(coef(fm1)[1]+log10(xlim)*coef(fm1)[2]),lty=2) +lines(xlim,10^(stats::coef(fm1)[1]+log10(xlim)*stats::coef(fm1)[2]),lty=2) #print(xlim) #print(10^(coef(fm1)[1]+log10(xlim)*coef(fm1)[2])) if (annote) { diff --git a/R/mfft.R b/R/mfft.R index e0a1579..055f6ed 100644 --- a/R/mfft.R +++ b/R/mfft.R @@ -11,6 +11,8 @@ #' estimates the frequencies (f_j), amplitudes (A_j) and phases #' (psi_j) in its decomposition. #' @useDynLib gtseries +#' @importFrom stats start +#' @importFrom stats as.ts #' @param xdata The data provided either as a time series (advised), or as a vector. #' may be complex #' @param minfreq,maxfreq If provided, bracket the frequencies to be probed. Note this are @@ -25,7 +27,6 @@ #' The computed frequencies are in the range given by minfreq and maxfreq. #' @param nfreq is the number of frequencies returned, must be smaller that the length of xdata. #' @return STILL NEED TO BE DESCRIBED -#' @references #' @author Michel Crucifix for the R code, and David Nesvorny for most of the supporting C code doing the #' actual computations #' \insertRef{sidlichovsky97aa}{gtseries} @@ -34,11 +35,11 @@ mfft <- function(xdata, minfreq=NULL, maxfreq=NULL, flag=1, nfreq=30) { - xdata = as.ts(xdata) + xdata = stats::as.ts(xdata) dt = deltat(xdata) print('deltat') print(dt) - startx = start(xdata) + startx = stats::start(xdata) ydata = Im(xdata) xdata = Re(xdata) ndata = length(xdata); diff --git a/R/phase_randomize.R b/R/phase_randomize.R index 1164e74..efdc3ac 100644 --- a/R/phase_randomize.R +++ b/R/phase_randomize.R @@ -5,6 +5,7 @@ #' complex numbers, but preserving the complex conjugate relationship guarantee #' that the inverse transform will be real #' @param x the input time series +#' @importFrom stats runif #' @return x the output time series. #' phase_randomize <- function (x) @@ -15,7 +16,7 @@ X <- fft(x) k <- seq(2,floor(N/2)+1) -X[k] <- Mod(X[k])*exp(2*pi*1i*runif(length(k))) +X[k] <- Mod(X[k])*exp(2*pi*1i*stats::runif(length(k))) X[N-k+2] <- Conj(X[k]) if (N == 2*floor(N/2)) {k <- N/2+1 diff --git a/R/prettylab.R b/R/prettylab.R index 42f7bed..5e694bf 100644 --- a/R/prettylab.R +++ b/R/prettylab.R @@ -2,93 +2,98 @@ ### Help files: ../man/pretty10exp.Rd ../man/axTexpr.Rd ### -------------- ---------- +#' @importFrom sfsmisc axTexpr +#' @importFrom sfsmisc prettylab +#' @importFrom sfsmisc pretty10exp +#' @importFrom sfsmisc eaxis -pretty10exp <- function(x, drop.1 = FALSE) -{ - ## Purpose: produce "a 10^k" label expressions instead of "a e<k>" - ## ---------------------------------------------------------------------- - ## Arguments: x: numeric vector (e.g. axis tick locations) - ## ---------------------------------------------------------------------- - ## Author: Martin Maechler, Date: 7 May 2004; 24 Jan 2006 - eT <- floor(log10(abs(x))) # x == 0 case is dealt with below - mT <- x / 10^eT - ss <- vector("list", length(x)) - for(i in seq(along = x)) - ss[[i]] <- - if(x[i] == 0) quote(0) - else if(drop.1 && (abs(eT[i])< 3)) substitute(E, list(E=x[i])) - else if(drop.1 && mT[i] == 1) substitute( 10^E, list(E = eT[i])) - else if(drop.1 && mT[i] == -1) substitute(-10^E, list(E = eT[i])) - else substitute(A %*% 10^E, list(A = mT[i], E = eT[i])) - do.call("expression", ss) -} - -axTexpr <- function(side, at = axTicks(side, axp=axp, usr=usr, log=log), - axp = NULL, usr = NULL, log = NULL, drop.1 = FALSE) -{ - ## Purpose: Do "a 10^k" labeling instead of "a e<k>" - ## ------------------------------------------------- - ## Arguments: as for axTicks() - pretty10exp(at, drop.1) -} - -### TODO: -### -### Myaxis(.) function with at least two options ("engineering/not") -### Really wanted: allow xaxt = "p" (pretty) or "P" (pretty, "Engineer") - -eaxis <- function(side, at = axTicks(side, log=log), labels = NULL, log = NULL, - f.smalltcl = 3/5, at.small = NULL, small.mult = NULL, - outer.at = TRUE, drop.1 = TRUE, las = 1, ...) -{ - ## Purpose: "E"xtended, "E"ngineer-like (log-)axis - ## ---------------------------------------------------------------------- - ## Author: Martin Maechler, Date: 13 Oct 2007 - is.x <- side%%2 == 1 - if(is.null(log)) { - XY <- function(ch) paste(if (is.x) "x" else "y", ch, sep = "") - log <- par(XY("log")) - } - ## use expression (i.e. plotmath) if 'log' or exponential format: - use.expr <- log || format.info(as.numeric(at), digits=7)[3] > 0 - if(is.null(labels)) - labels <- if(use.expr) pretty10exp(at, drop.1=drop.1) else TRUE - else if(length(labels) == 1 && is.na(labels)) # no 'plotmath' - labels <- TRUE - axis(side, at = at, labels = labels, las=las, ...) - if(is.null(at.small)) { ## create smart default, using small.mult - at.small <- - if(log) { - if(is.null(small.mult)) small.mult <- 9 - at. <- at[log10(at) %% 1 < 1e-3] ## the 10^k ones: - if(length(at.)) - outer(2:small.mult, c(if(outer.at) at.[1]/10, at.)) - } else { - ## assumes that 'at' is equidistant - d <- diff(at <- sort(at)) - if(any(abs(diff(d)) > 1e-3 * (dd <- mean(d)))) - stop("'at' is not equidistant") - if(is.null(small.mult)) { - ## look at 'dd' , e.g. in {5, 50, 0.05, 0.02 ..} - d. <- dd / 10^floor(log10(dd)) - small.mult <- { - if(d. %% 5 == 0) 5 - else if(d. %% 4 == 0) 4 - else if(d. %% 2 == 0) 2 - else if(d. %% 3 == 0) 3 - else if(d. %% 0.5 == 0) 5 - else 2 } - } - outer(1:(small.mult-1)/small.mult * dd, - c(if(outer.at) at[1]-dd, at), "+") - } - ## - if(outer.at) { # make sure 'at.small' remain inside "usr" - p.u <- sort(par("usr")[if(is.x) 1:2 else 3:4]) - if(log) p.u <- 10^p.u - at.small <- at.small[p.u[1] <= at.small & at.small <= p.u[2]] - } - } - if (outer.at) {axis(side, at = at.small, label = FALSE, - tcl = f.smalltcl * par("tcl"))} -} +## pretty10exp <- function(x, drop.1 = FALSE) +## { +## ## Purpose: produce "a 10^k" label expressions instead of "a e<k>" +## ## ---------------------------------------------------------------------- +## ## Arguments: x: numeric vector (e.g. axis tick locations) +## ## ---------------------------------------------------------------------- +## ## Author: Martin Maechler, Date: 7 May 2004; 24 Jan 2006 +## eT <- floor(log10(abs(x))) # x == 0 case is dealt with below +## mT <- x / 10^eT +## ss <- vector("list", length(x)) +## for(i in seq(along = x)) +## ss[[i]] <- +## if(x[i] == 0) quote(0) +## else if(drop.1 && (abs(eT[i])< 3)) substitute(E, list(E=x[i])) +## else if(drop.1 && mT[i] == 1) substitute( 10^E, list(E = eT[i])) +## else if(drop.1 && mT[i] == -1) substitute(-10^E, list(E = eT[i])) +## else substitute(A %*% 10^E, list(A = mT[i], E = eT[i])) +## do.call("expression", ss) +## } +## +## #' export axTexpr +## axTexpr <- function(side, at = axTicks(side, axp=axp, usr=usr, log=log), +## axp = NULL, usr = NULL, log = NULL, drop.1 = FALSE) +## { +## ## Purpose: Do "a 10^k" labeling instead of "a e<k>" +## ## ------------------------------------------------- +## ## Arguments: as for axTicks() +## pretty10exp(at, drop.1) +## } +## +## ### TODO: +## ### +## ### Myaxis(.) function with at least two options ("engineering/not") +## ### Really wanted: allow xaxt = "p" (pretty) or "P" (pretty, "Engineer") +## +## eaxis <- function(side, at = axTicks(side, log=log), labels = NULL, log = NULL, +## f.smalltcl = 3/5, at.small = NULL, small.mult = NULL, +## outer.at = TRUE, drop.1 = TRUE, las = 1, ...) +## { +## ## Purpose: "E"xtended, "E"ngineer-like (log-)axis +## ## ---------------------------------------------------------------------- +## ## Author: Martin Maechler, Date: 13 Oct 2007 +## is.x <- side%%2 == 1 +## if(is.null(log)) { +## XY <- function(ch) paste(if (is.x) "x" else "y", ch, sep = "") +## log <- par(XY("log")) +## } +## ## use expression (i.e. plotmath) if 'log' or exponential format: +## use.expr <- log || format.info(as.numeric(at), digits=7)[3] > 0 +## if(is.null(labels)) +## labels <- if(use.expr) pretty10exp(at, drop.1=drop.1) else TRUE +## else if(length(labels) == 1 && is.na(labels)) # no 'plotmath' +## labels <- TRUE +## axis(side, at = at, labels = labels, las=las, ...) +## if(is.null(at.small)) { ## create smart default, using small.mult +## at.small <- +## if(log) { +## if(is.null(small.mult)) small.mult <- 9 +## at. <- at[log10(at) %% 1 < 1e-3] ## the 10^k ones: +## if(length(at.)) +## outer(2:small.mult, c(if(outer.at) at.[1]/10, at.)) +## } else { +## ## assumes that 'at' is equidistant +## d <- diff(at <- sort(at)) +## if(any(abs(diff(d)) > 1e-3 * (dd <- mean(d)))) +## stop("'at' is not equidistant") +## if(is.null(small.mult)) { +## ## look at 'dd' , e.g. in {5, 50, 0.05, 0.02 ..} +## d. <- dd / 10^floor(log10(dd)) +## small.mult <- { +## if(d. %% 5 == 0) 5 +## else if(d. %% 4 == 0) 4 +## else if(d. %% 2 == 0) 2 +## else if(d. %% 3 == 0) 3 +## else if(d. %% 0.5 == 0) 5 +## else 2 } +## } +## outer(1:(small.mult-1)/small.mult * dd, +## c(if(outer.at) at[1]-dd, at), "+") +## } +## ## +## if(outer.at) { # make sure 'at.small' remain inside "usr" +## p.u <- sort(par("usr")[if(is.x) 1:2 else 3:4]) +## if(log) p.u <- 10^p.u +## at.small <- at.small[p.u[1] <= at.small & at.small <= p.u[2]] +## } +## } +## if (outer.at) {axis(side, at = at.small, label = FALSE, +## tcl = f.smalltcl * par("tcl"))} +## } diff --git a/R/ssa.R b/R/ssa.R index d08d2d5..d1cc93b 100644 --- a/R/ssa.R +++ b/R/ssa.R @@ -19,6 +19,9 @@ #' @param N if provided, the length of the input series to beconsiderd #' @param M the length of the sliding window #' @param Nrec: number of components to be output +#' @importFrom sfsmisc axTexpr +#' @importFrom sfsmisc eaxis + #' @references #' \insertRef{Ghil02aa}{gtseries} #' \insertRef{VAUTARD89aa}{gtseries} @@ -156,6 +159,7 @@ ssa <- function(X=X,N=length(X),M=M,Nrec=10) { #' @export plot.SSAObject <- function (SSAObject,...) { + Nrec = length(SSAAbject$LA) yat <- outer(seq(1:Nrec),10^(seq(-3,5)),"*") ylab <- 10^(-2:5) diff --git a/README.md b/README.md index 03a2f25..e69de29 100644 --- a/README.md +++ b/README.md @@ -1,6 +0,0 @@ - -Work in progress. Do not cite, yet. - -_Acknowledgement:_ - -`prettylab.R ` is taken from _fsmisc_: [Utilities from 'Seminar fuer Statistik' ETH Zurich](https://rdrr.io/cran/sfsmisc/src/R/prettylab.R) diff --git a/cleanup b/cleanup new file mode 100755 index 0000000..5c8d055 --- /dev/null +++ b/cleanup @@ -0,0 +1,3 @@ +#!/bin/sh + +rm -f config.* src/Makevars src/config.h diff --git a/inst/REFERENCES.bib b/inst/REFERENCES.bib index 1fff25d..574a6ac 100644 --- a/inst/REFERENCES.bib +++ b/inst/REFERENCES.bib @@ -33,3 +33,48 @@ publisher = {Springer Netherlands}, isbn = {9789401155106}, } +@phdthesis{Burg75aa, + author = {Burg, J. P. }, + editor = {}, + school = "Department of Geophysics, Stanford University", + year = {1975}, + title = {Maximum Entropy Spectral Analysis}, + url = "https://sepwww.stanford.edu/data/media/public/oldreports/sep06/" +} + + +@article{pardo-iguzquiza21aa, + author = {Pardo-Igúzquiza, E. and RodrÃguez-Tovar, F. J. }, + editor = {}, + year = {2021}, + title = {Maximum Entropy Spectral Analysis}, + pages = {1–8}, + journal = {Encyclopedia of Earth Sciences Series}, + doi = {10.1007/978-3-030-26050-7_197-1}, + url = {http://dx.doi.org/10.1007/978-3-030-26050-7_197-1}, + publisher = {Springer International Publishing}, + issn = {1871-756X}, + isbn = {9783030260507}, +} +@article{pardoiguzquiza06aa, + author = {Pardoâ€Igúzquiza, E. and RodrÃguezâ€Tovar, F. J. }, + editor = {}, + year = {2006}, + title = {Maximum entropy spectral analysis of climatic time series revisited: Assessing the statistical significance of estimated spectral peaks}, + journal = {Journal of Geophysical Research: Atmospheres}, + volume = {111}, + doi = {10.1029/2005jd006293}, + url = {http://dx.doi.org/10.1029/2005JD006293}, + publisher = {American Geophysical Union (AGU)}, + number = {D10}, + issn = {0148-0227}, +} + + +@book{percival98, + author = {Percival, D. B. and Walden, A. }, + editor = {}, + year = {1998}, + publisher = "Cambridge University Press", + title = {Spectral analysis for physical applications: multitaper and conventional univariate techniques}, +} diff --git a/man/axTexpr.Rd b/man/axTexpr.Rd deleted file mode 100644 index 2adb4bc..0000000 --- a/man/axTexpr.Rd +++ /dev/null @@ -1,67 +0,0 @@ -\name{axTexpr} -\alias{axTexpr} -\title{Axis Ticks Expressions in Nice 10 ** k Form} -\description{ - Produce nice \eqn{a \times 10^k}{a * 10^k} expressions for - \code{\link{axis}} labeling instead of the scientific notation - \code{"a E<k>"}. - -} -\usage{ -axTexpr(side, at = axTicks(side, axp = axp, usr = usr, log = log), - axp = NULL, usr = NULL, log = NULL, - drop.1 = FALSE) -} -\arguments{ - \item{side}{integer in 1:4 specifying the axis side, as for - \code{\link{axis}}.} - \item{at}{numeric vector; with identical default as in - \code{\link{axTicks}()}.} - \item{axp, usr, log}{as for \code{\link{axTicks}()}.} - \item{drop.1}{logical indicating if \eqn{1 \times}{1 *} should be - dropped from the resulting expressions.} -} -\details{ - This is just a utility with the same arguments as - \code{\link{axTicks}}, calling \code{\link{pretty10exp}(at, *)}. -} -\value{ - an expression of the same length as \code{x}, with elements of the - form \code{a \%*\% 10 ^ k}. -} -\author{Martin Maechler} -\seealso{\code{\link{pretty10exp}}; - \code{\link{axis}}, \code{\link{axTicks}}. -} -\examples{ -x <- 1e7*(-10:50) -y <- dnorm(x, m=10e7, s=20e7) -plot(x,y)## not really nice, the following is better: - -## For horizontal y-axis labels, need more space: -op <- par(mar= .1+ c(5,5,4,1)) -plot(x,y, axes= FALSE, frame=TRUE) -aX <- axTicks(1); axis(1, at=aX, label= axTexpr(1, aX)) -## horizontal labels on y-axis: -aY <- axTicks(2); axis(2, at=aY, label= axTexpr(2, aY), las=2) -par(op) - -### -- only 'x' and using log-scale there: -plot(x,y, xaxt= "n", log = "x") -aX <- axTicks(1); axis(1, at=aX, label= axTexpr(1, aX)) - -## Now an `` engineer's version '' ( more ticks; only label "10 ^ k" ) : - -axp <- par("xaxp") #-> powers of 10 *inside* 'usr' -axp[3] <- 1 # such that only 10^. are labeled -aX <- axTicks(1, axp = axp) -xu <- 10 ^ par("usr")[1:2] -e10 <- c(-1,1) + round(log10(axp[1:2])) ## exponents of 10 *outside* 'usr' -v <- c(outer(1:9, e10[1]:e10[2], function(x,E) x * 10 ^ E)) -v <- v[xu[1] <= v & v <= xu[2]] - -plot(x,y, xaxt= "n", log = "x", main = "engineer's version of x - axis") -axis(1, at = aX, label = axTexpr(1, aX, drop.1=TRUE)) # `default' -axis(1, at = v, label = FALSE, tcl = 2/3 * par("tcl")) -} -\keyword{dplot} diff --git a/man/mem.Rd b/man/mem.Rd index 2dde547..dadb908 100644 --- a/man/mem.Rd +++ b/man/mem.Rd @@ -2,10 +2,30 @@ % Please edit documentation in R/mem.R \name{mem} \alias{mem} +\alias{plot.memObject} \title{Maximum entropy method} \usage{ mem(A, method = "burg", order.max = 90, ...) + +\method{plot}{memObject}(memObject, period = FALSE, xaxp = NULL, yaxp = NULL, ...) +} +\arguments{ +\item{A}{input, a time-series object} + +\item{method}{Estimation method, passed to the standard `ar` function, which is one of "yule-walker", "burg", "ols", "mle", "yw". Defaults to "burg".} + +\item{order.max}{maximum order of the autoregressive model being fitted} +} +\value{ +a memObject } \description{ -Maximum entropy method +Maximum entropy spectral analysis (MESA) is the statistical estimation of the power spectrum of a stationary time series using the maximum entropy (ME) method. The resulting ME spectral estimator is parametric and equivalent to an autoregressive (AR) spectral estimator (copied from \\insertCite{pardo-iguzquiza21aa}{gtseries}). The method relies heavily on the standard-R `ar` function, of which it is essentially a wrapper associated with a plot class function (authors: Martyn Plummer. C code for `ar.burg` by B. D. Ripley.) +} +\references{ +\insertRef{Burg75aa}{gtseries} +\insertRef{percival98}{gtseries} +\insertRef{Ghil02aa}{gtseries} +\insertRef{pardoiguzquiza06aa}{gtseries} +\insertRef{pardo-iguzquiza21aa}{gtseries} } diff --git a/man/pretty10exp.Rd b/man/pretty10exp.Rd deleted file mode 100644 index ceda1e0..0000000 --- a/man/pretty10exp.Rd +++ /dev/null @@ -1,32 +0,0 @@ -\name{pretty10exp} -\alias{pretty10exp} -\title{Nice 10 ** k Label Expressions} -\description{ - Produce nice \eqn{a \times 10^k}{a * 10^k} expressions to be used - instead of the scientific notation \code{"a E<k>"}. -} -\usage{ -pretty10exp(x, drop.1 = FALSE) -} -\arguments{ - \item{x}{numeric vector (e.g. axis tick locations)} - \item{drop.1}{logical indicating if \eqn{1 \times}{1 *} should be - dropped from the resulting expressions.} -} -\value{ - an expression of the same length as \code{x}, with elements of the - form \code{a \%*\% 10 ^ k}. -} -\author{Martin Maechler} -\seealso{\code{\link{axTexpr}} which builds on \code{pretty10exp()}; - further \code{\link{axis}}, \code{\link{axTicks}}. -} -\examples{ -pretty10exp(-1:3 * 1000) -pretty10exp(-1:3 * 1000, drop.1 = TRUE) -pretty10exp(c(1,2,5,10,20,50,100,200) * 1e3) -\dontshow{ -stopifnot(identical(pretty10exp(numeric(0)), expression())) -} -} -\keyword{dplot} diff --git a/src/fmft.c b/src/fmft.c index d7e9cb3..635e7ca 100644 --- a/src/fmft.c +++ b/src/fmft.c @@ -17,6 +17,15 @@ X(t) + iY(t) = Sum_j=1^N [ A_j * exp i (f_j * t + psi_j) ] */ /*#define N2FLAG */ /*#define N2FLAG2*/ +#include <R.h> +// Note: This loads R_ext/Print.h that we need + +// Define strict headers +#define STRICT_R_HEADERS +// Map printf to Rprintf +#define printf Rprintf +#define printe Reprintf + #define FMFT_TOL 1.0e-10 /* MFT NOMINAL PRECISION */ #define FMFT_NEAR 0. /* MFT OVERLAP EXCLUSION PARAMETER */ diff --git a/src/nrutil.c b/src/nrutil.c index 38d3025..8b46fb2 100644 --- a/src/nrutil.c +++ b/src/nrutil.c @@ -1,15 +1,28 @@ /* NRUTIL: from "Numerical Recipes in C", pp. 705-709 */ +/* released in "public domain" whatever it means */ +/* modified by introducing R.h headers */ #include <stdio.h> #include <stdlib.h> #include <malloc.h> +#include <R.h> +// Note: This loads R_ext/Print.h that we need + +// Define strict headers +#define STRICT_R_HEADERS +// Map printf to Rprintf +#define printf Rprintf +#define printerr REprintf + + + /* nrerror() := Numerical Recipes standard error handler */ void nrerror(char *error_text) { - fprintf(stderr, "Numerical Recipes run-time error...\n"); - fprintf(stderr, "%s\n", error_text); - fprintf(stderr, "...now exiting to system...\n"); - exit(1); + printerr("C-code run-time error...\n"); + printerr("%s\n", error_text); + printerr("...now exiting to system...\n"); +/* exit(1);*/ } /* *vector() := Allocates a float vector with range [nl..nh] */ -- GitLab