Initial commit

CHANGES

+1.03 : added cross morlet

DESCRIPTION

+Package: gtseries
+Type: Package
+Title: Time Series analysis for galcial cycles
+Version: 1.03
+Date: 2015-02-12
+Author: Michel Crucifix
+Maintainer: Michel Crucifix <michel.crucifix@uclouvain.be>
+Description: More about what it does (maybe more than one line)
+License: What license is it under?
+LazyLoad: yes

NAMESPACE

+exportPattern("^[[:alpha:]]+")

R/approx_ts.R

+  approx_ts <- function (data,tcoord="x",dcoord="y",scale=1,n=2048,thin=FALSE,spline=FALSE,xlim=NULL)
+  {
+  local({
+   x <- data[tcoord][[1]] 
+   y <- data[dcoord][[1]] 
+   if (!is.null(xlim)) {T <- which(x>=xlim[1] & x <= xlim[2]) ; x<-x[T]; y <- y[T]}
+   deltat <- diff(range(x))/n
+   ## thin data according with thin factor deltat
+
+   if (thin) {
+   thin_f <- 0.5*deltat
+   tmp <-x
+   while (TRUE) { dd <- abs(diff(tmp))
+    if (any(drop <- dd < thin_f))
+      tmp <- tmp[-(1 + which(drop)[1])]
+    else
+      break
+      }
+
+    y <- y[match(tmp,x)]
+    x <- x[match(tmp,x)]
+    }
+
+    ## then interpolates
+
+    dummy <-  eval(call(ifelse(spline,"spline","approx"),x,y,n=n))
+    out <<- ts(dummy$y,start=dummy$x[1]/scale, deltat=diff(dummy$x[1:2])/scale)
+   })
+
+  out
+  }

R/arspec.R

+## Calculates the theoretical density frequency of an
+## autoregressive process, according to equation
+## (3.3) of Akaike (1969)
+
+arspec <- function (f,ar){
+Exp <- exp(-1i*2*pi*outer(f,seq_along(ar),"*"))
+1./abs(1-Exp%*%ar)^2
+}
+

R/cwt_morlet.R

+#This program performs the Continuous Wavelet Transform (CWT) 
+#of the input time series y. 
+#It plots the series in normalized form 
+#and displays the modulus (amplitude) of the CWT 
+#in the time-period space. 
+#The period is expressed in unit of time.
+#
+#
+#   are in the vector named "period" 
+
+#reference:
+#
+#Mallat, S. 1998: A wavelet Tour of Signal Processing. 
+#Academic Press, 577 pp.
+#
+## for ridge extraction : 
+## Tchawitchia P., Wavelets, functions and operators, ch. 3 in 
+## Wavelets : theory and applications, Erlenbacher et al. eds, Oxford University Press 1996
+
+col_wavelet <- colorRampPalette(c('darkblue','lightblue','grey','white','red'))(100)
+
+search_ridge <- function (A,first_guess,k0=5.6,tol=1e-6, itmax=20,B=seq(along=A),window=c(0,Inf))
+{
+  n <- length(A)
+  Amp <- A
+  Amp[] <- NA
+  Ome <- Amp
+  Scale <- Amp
+  deltat <- deltat(A)
+   scale <- first_guess
+   for (b in B)
+   {
+     delta = Inf
+     it <- 0
+     while ((delta > tol) & (it <= itmax)) {
+         w <- cwt_morlet(A,k0=k0,calcmask=FALSE,scale=scale,deriv=TRUE)
+         p <- attr(w,"deriv")[b]
+         expected_scale = k0/Re(p)
+         delta <- scale - expected_scale 
+         scale <- expected_scale 
+         it <- it+1
+         if (is.na(scale)) {it = itmax}
+       }
+       if (it < itmax &  scale > window[1] & scale < window[2]  ) { 
+           Amp[b] = w[b]
+           Ome[b] = p 
+           Scale[b] = scale
+           }
+       else {     Amp[b] = NA
+                  Ome[b] = NA 
+                  Scale[b] = NA 
+                  if (scale > window[2]) scale = window[2] 
+                  if (scale < window[1]) scale = window[1] 
+            }      
+   } 
+
+    
+   attr(Amp,"period") = 2*pi*deltat/Ome
+   attr(Amp,"scale") = Scale
+   attr(Amp,"k0") = k0
+   attr(Amp,"class") = "ridge"
+   Amp
+}
+
+enhance_ridge <- function(A,R,confidence=0.95,inter=0.08)
+{
+ Scale = attr(R,"scale")
+ Amp = R
+ period=matrix(NA,nrow=length(A),ncol=3)
+ k0 = attr(R,"k0")
+ deltat = deltat(A)
+
+
+ width=log(qnorm(0.5+confidence/2)*sqrt(2)/k0+1)
+
+ for (b in seq(along=R))
+  {
+    print(b)
+    if (!is.na(Scale[b]))
+    {
+      scale=exp(log(Scale[b])+seq(-width*1.03214,+width,inter))
+      ## 1.03214 is an empirical factor to take into account the ridge asymmetry 
+      w <- cwt_morlet(A,k0=k0,calcmask=FALSE,scale=scale,deriv=TRUE)
+      Amp[b]=sum(w[b,])*k0/sqrt(2*pi)*inter
+      period[b,] = 2*pi*deltat/k0*c(Scale[b],min(scale),max(scale))
+    }
+   }
+    attr(Amp,"period") <- period
+    Amp
+}
+
+
+reconstruct_morlet  <- function(W,scales=c(-Inf,Inf), periods=NULL)
+{
+ if (!(attr(W,"class")=="wavelet")) stop ("object is not a wavelet transform")
+ if (!(attr(W,"wavelet")=="morlet")) stop ("object is not a MORLET  wavelet transform")
+ a <- attr(W,"scale")
+ j <- which(a > scales[1] & a < scales[2])
+ # if period is set, then overrides scales
+ if (!is.null(periods))
+ {
+   a <- attr(W,"period")
+   j <- which(a > periods[1] & a < periods[2])
+ }
+
+ inter <- attr(W,"parameters")$inter
+ k0 <- attr(W,"parameters")$k0
+
+ T <- rowSums(W[,j])*log(2)/inter/(sqrt(2*pi))*k0
+ T <- ts(T, start=attr(W,'time')[1], deltat = diff(attr(W,'time'))[1])
+}
+
+cross_morlet <- function(A, B, ...)
+{
+  CA <- cwt_morlet(A, ...)
+  CB <- cwt_morlet(B, ...)
+  (CA * Conj(CB)) / (Mod(CA) * Mod(CB) )
+}
+
+cwt_morlet <- function (A,inter=20,k0=5.6,amin=1,amax=Inf,calcmask=TRUE,scale=NA,deriv=FALSE)
+{
+   y <- A
+   xx <- time(A)
+   deltat<-deltat(A)
+   ny<-length(y);
+
+  if (is.na(scale)) {
+    local({
+      ny2<-round(ny/2);
+      exp1<-log2(amin)+2;
+      exp2<-min(round(log2(ny2))+1,amax);
+
+      scale <- vector()
+
+      j<-0;
+      for (m in seq(exp1,exp2-1))
+      {
+        jj<-inter-1;
+        for (n in seq(0,jj)) 
+        {
+          a<-2^(m+n/inter);
+          j<-j+1;
+          scale[j]<<-a;
+        }
+    }})
+    } 
+   
+
+## now infers the corresponding periods
+ ##  omega0<-1/2*(k0/aa+sqrt(2+k0*k0)/aa);
+  omega0<-k0/scale;
+  period<-1./omega0*2*pi*deltat;
+
+  x<-y
+  n<-length(x);
+
+  k<-(0:(n-1))*2*pi/n
+
+  f<-fft(x);
+
+  J<-length(scale);
+
+  wave<-matrix(as.complex(0),nrow=n,ncol=J);
+
+  if (calcmask) mask<-matrix(TRUE,nrow=n,ncol=J);
+  if (deriv) dwave <- wave
+
+  nn <-length(k);
+  for(a1 in seq(along = scale)) {
+     expnt<- -(scale[a1]*k - k0) ^2/2; ## psi_hat(a*k)
+ ##   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;
+    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;
+    }
+  }
+    if (deriv)  {dwave <- dwave / (-1i*wave) }
+    xx <-  as.array(xx);
+    yy <- as.array(period)
+    attr(wave,"time") <- xx
+    if (calcmask) attr(wave,"mask") <- mask
+    attr(wave,"period") <- yy
+    attr(wave,"scale") <- scale
+    attr(wave,"class") <- "wavelet"
+    attr(wave,"wavelet") <- "morlet"
+    attr(wave,"parameters") <- list(k0=k0,inter=inter,deltat=deltat)
+    if (deriv) attr(wave,"deriv") <- dwave
+    wave
+}
+
+plot.wavelet <- function (wave,resx=400,resy=300,xlab="Time",ylab="Period",scaling_correction=0,col=col_wavelet,legend=FALSE,Mode=Mod,...)
+
+{
+  require(fields)
+  xx <- attr(wave,"time")
+  period <- attr(wave,"period")
+  mask <- attr(wave,"mask")
+  aa <- attr(wave,"scale")
+  thin_factor_xx <- max(ceiling(length(xx)/ resx),1)
+  thin_factor_yy <- max(ceiling(length(aa)/ resy),1)
+  subx <- seq(thin_factor_xx,length(xx),thin_factor_xx)
+  suba <- seq(thin_factor_yy,length(aa),thin_factor_yy)
+  wave_scaled <- Mode(wave[subx,suba])*mask[subx,suba]
+  for (i in 1:length(suba)) wave_scaled[,i] <- wave_scaled[,i]/(aa[i]^scaling_correction)
+  if (legend) par(oma=c(2,2,2,5))
+  image(xx[subx],period[suba],wave_scaled,log="y",ylab=ylab,xlab=xlab,axes=FALSE,col=col,...)
+  axis(1)
+  axis.log10(2,"")
+  if (legend) {par(oma=c(2,2,2,2))
+  image.plot(xx[subx],period[suba],wave_scaled,legend.only=TRUE,log="y",ylab=ylab,xlab=xlab,axes=FALSE,col=col,...)
+  par(oma=c(2,2,2,5))}
+}
+
+
+powerspectrum.wavelet <- function (wave,...)
+{
+  aa <- attr(wave,"scale")
+  yy <- attr(wave,"period")
+  xx <- attr(wave,"time")
+  f <- 1./yy
+  P <- vector("double",length(f))
+for (j in seq_along(f)) P[j] <- mean(Mod(wave[,j])^2,na.rm=TRUE)
+  data.frame("frequency"=f,"power"=P)
+}
+
+analytic.ridge <- function (wave,ridge_only = FALSE,plim=c(-Inf,+Inf))
+{
+  require("signal")
+  p <- attr(wave,"period")
+  P <- length(p)
+  pout <- which(p<plim[1] | p>plim[2])
+  s <- attr(wave,"scale")
+  t <- attr(wave,"time")
+  deltat <- diff(attr(wave,"time"))[1]
+  WR <- wave
+  WR[,pout] <- NA
+  # scale for scalogram
+  for (i in seq(1,P)) WR[,i] <- WR[,i]/sqrt(s[i])
+
+  WP <- Arg(WR)
+  for (i in seq(1,P) ) WP[,i] <- unwrap(WP[,i])
+  for (i in seq(1,P) ) WP[,i] <- c(diff(WP[,i]),NA)/(2*pi)/deltat*p[i]
+  WR[which(abs(WP-1) > 0.025)] <- NA
+  if (ridge_only) 
+  { 
+   R <- vector()
+   for (i in seq(1,length(t))) 
+    {M <- which.max(Mod(WR[i,])) ; if (any(M)) R[i] <- WR[i,M] }
+   R
+  } else 
+  WR
+
+}

R/mem.R

+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,...)
+
+  ##  v.maice multiplied by deltat to have power density (I still need to really undersant
+  ##    why, but this works)
+  arc <- sapply (U$ar,function (x) {x})
+  mem <- data.frame(frequency=f/deltat,power= arspec(f,arc)*U$var.pred*(deltat),row.names=NULL)
+
+  attr(mem,"var") <- U$v.maice*(deltat)
+  attr(mem,"L") <- U$order
+  attr(mem,"class") <- c("memObject","data.frame")
+
+  mem
+}
+
+
+plot.memObject <- function (memObject,period=FALSE,xaxp=NULL,yaxp=NULL,...)
+{
+
+  plot(memObject$frequency,memObject$power,
+       frame=FALSE,axes=FALSE,log="xy",type="l",xlab="",ylab="Power density",...)
+  local({
+  if (is.null(yaxp)) {yaxp <- par("yaxp"); yaxp[3] <- 1; par(yaxp=yaxp)}
+  if (is.null(xaxp)) {xaxp <- par("xaxp"); xaxp[3] <- 1; par(xaxp=xaxp)}
+   eaxis(2,outer.at=FALSE)
+  
+  if (period) {axis.period()} else {axis.frequency()}
+  })
+}
+
+axis.period    <- function(side=1,...) { axis.log10(side,"Period",factor=-1,...) }
+axis.frequency <- function(side=1,...) { axis.log10(side,"Frequency",factor=1,...) }
+
+
+axis.log10 <- function(side=1,title="Power",factor=1,line=3,outer.at=TRUE,small.tcl=3/5,las=1,...)
+{
+   range10 <- sort(factor*par("usr")[c(3,4)-2*(side %% 2)])
+   range <- 10^range10
+
+   big_ticks <- 10^seq(floor(range10[1]),ceiling(range10[2]))
+   big_ticks <- big_ticks[which(big_ticks > range[1] & big_ticks < range[2])]
+   axis(side,at=big_ticks^factor,labels=axTexpr(side=side,at=big_ticks,drop.1=TRUE),las=las,...)
+
+   if (outer.at ) {
+         small_ticks <- outer(seq(1,9),10^seq(floor(range10[1]),ceiling(range10[2])),"*")
+         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"))
+}
+
+add.intercept <- function (memObject, xlim = range(memObject$frequency)[2]*c(5.e-2,0.5),annote=TRUE)
+{
+frequency <- memObject$frequency
+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)
+#print(xlim)
+#print(10^(coef(fm1)[1]+log10(xlim)*coef(fm1)[2]))
+if (annote) {
+text(mean(frequency)+0.1*diff(range(frequency)),mean(memObject$power)+0.0*diff(range(memObject$power)),
+    paste("slope = ",format(coef(fm1)[[2]],digit=2)))}
+coef(fm1)[[2]]
+}
+
+
+
+
+

R/phase_randomize.R

+## randomize the phases of a real signal, for non-linear time series analysis
+
+phase_randomize <- function (x)
+{
+N <- length(x)
+
+X <- fft(x)
+
+k <- seq(2,floor(N/2)+1)
+
+X[k] <- Mod(X[k])*exp(2*pi*1i*runif(length(k)))
+X[N-k+2] <- Conj(X[k])
+
+if (N == 2*floor(N/2)) {k <- N/2+1
+                        X[k] <- Mod(X[k])}
+
+y <- fft(X,inverse=TRUE)
+Re(y)/N
+}
+

R/prettylab.R

+### --> these are from ~/R/MM/GRAPHICS/axis-prettylab.R
+
+### Help files: ../man/pretty10exp.Rd  ../man/axTexpr.Rd
+###                    --------------         ----------
+
+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"))}
+}

R/ssa.R

+# SSA as in Ghil et al., 
+# Advances spectral methods for climatic time series
+# Rev Geoph., 40 (1), art. 3 (2002)
+# doi 10.1029/2000RG000092
+
+# Has been validated and show to produce EXACTLY the same results
+# as the Spectra toolkin on 25.9.08
+#----------------------------------------------------------------------   
+# using astronomical forcing as sample data
+
+ssa <- function(X=X,N=length(X),M=M)
+{
+
+  # construct  Toeplitz Matrix for reference signal (eq. 6)
+  # according to Vautard and Ghil, 1989
+
+  X <- X - mean(X)
+  cij <- vector("double",M)
+
+  for (i in 1:M) {for (j in i:M)
+     ij <- j-i
+     T <- N-ij
+     cij[ij+1]=crossprod(X[1:T],X[(1+ij):(T+ij)])/T
+     }
+
+  c <- toeplitz(cij)
+
+  # calculate eigenvalues of reference signal (eq. 7)
+
+  Eig <- eigen(c)
+  lambda <- Eig$values  ## they are immediately sorted. Thanks R!
+  rhok   <- Eig$vectors
+
+  # sort Lambda of reference signal 
+
+  #  [lambda,index]=sort(Lambdak);
+
+
+  # evaluation of the principal components of the original signal (eq. 10)
+   
+  A <- matrix(0,N-M+1,10)
+
+
+
+   for (k in 1:10) { for (t in 1:(N-M+1))
+     A[t,k] <- crossprod(X[t:(t+M-1)],rhok[,k])
+   }
+
+  # reconstruction from the eigenvectors (eq. 11)
+
+  Mt <- matrix(0,N,1)
+  Lt <- matrix(0,N,1)
+  Lt1 <- matrix(0,N,1)
+  Ut <- matrix(0,N,1)
+  Ut1 <- matrix(0,N,1)
+
+  for (t in 1:(M-1)) {
+    Mt[t]<-t;
+    Lt[t]<-1;
+    Ut[t]<-t}
+
+  Mt[M:(N-M)]<-M
+  Lt[M:(N-M)]<-1
+  Ut[M:(N-M)]<-M
+
+  for (t in (N-M+1):N) {
+    Mt[t]<-N-t+1
+    Lt[t]<-t-N+M
+    Ut[t]<-M}
+
+  PCA <- matrix(0,10,N);
+
+
+  for (t in 1:N) {for (k in 1:10)
+   {
+   PCA[k,t] <- PCA[k,t] + crossprod(A[t-Lt[t]:Ut[t]+1,k],rhok[Lt[t]:Ut[t],k]);
+  }}
+
+  LA <- vector("double",10)
+  for (k in 1:10) {
+   PCA[k,] <- PCA[k,]/Mt
+   LA[k]=crossprod(PCA[k,],X[1:N])
+  }
+
+  SSAObject<- list(lambda=lambda,LA=LA,A=A,rhok=rhok,cij=cij,PCA=PCA)
+  attr(SSAObject,"class") <- "SSAObject"
+  SSAObject
+  }
+
+
+
+plot.SSAObject <- function (SSAObject,...)
+{
+  yat <- outer(seq(1:10),10^(seq(-3,5)),"*")
+  ylab <- 10^(-2:5)
+
+   ss <- lapply(seq(along = ylab),
+                    function(i) if(ylab[i] == 0) quote(0) else
+                       substitute(10^E, list(E=log10(ylab[i]))))
+
+   S <- do.call("expression",ss)
+  plot(SSAObject$lambda,frame=T,axes=F,log="y",xlab="rank",ylab="Eigenvalue",...)
+  axis(1,xlab="rank")
+  axis(2,at=yat,ylab="Eigenvalue",labels=F)
+  axis(2,at=ylab,labels=axTexpr(1,at=ylab),tick=F,las=1)
+}
+
+

man/axTexpr.Rd

+\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}

man/pretty10exp.Rd

+\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}