[Dplr-commits] r701 - in pkg/dplR: . R man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Oct 24 19:00:51 CEST 2013
Author: mvkorpel
Date: 2013-10-24 19:00:50 +0200 (Thu, 24 Oct 2013)
New Revision: 701
Added:
pkg/dplR/R/redfit.R
pkg/dplR/man/print.redfit.Rd
pkg/dplR/man/redfit.Rd
pkg/dplR/src/redfit.c
Modified:
pkg/dplR/
pkg/dplR/ChangeLog
pkg/dplR/DESCRIPTION
pkg/dplR/NAMESPACE
pkg/dplR/R/common.interval.R
pkg/dplR/R/detrend.R
pkg/dplR/R/detrend.series.R
pkg/dplR/R/pointer.R
pkg/dplR/R/powt.R
pkg/dplR/R/rwi.stats.running.R
pkg/dplR/R/wavelet.plot.R
pkg/dplR/man/bai.in.Rd
pkg/dplR/man/bai.out.Rd
pkg/dplR/man/cms.Rd
pkg/dplR/man/detrend.Rd
pkg/dplR/man/detrend.series.Rd
pkg/dplR/man/dplR-package.Rd
pkg/dplR/man/ffcsaps.Rd
pkg/dplR/man/fill.internal.NA.Rd
pkg/dplR/man/hanning.Rd
pkg/dplR/man/rwi.stats.running.Rd
pkg/dplR/man/sea.Rd
pkg/dplR/man/series.rwl.plot.Rd
pkg/dplR/man/skel.plot.Rd
pkg/dplR/man/wavelet.plot.Rd
Log:
Merge redfit branch back into trunk
Property changes on: pkg/dplR
___________________________________________________________________
Modified: svn:mergeinfo
- /branches/dplR-R-2.15:466-506
+ /branches/dplR-R-2.15:466-506
/branches/redfit:662-700
Modified: pkg/dplR/ChangeLog
===================================================================
--- pkg/dplR/ChangeLog 2013-10-24 15:59:10 UTC (rev 700)
+++ pkg/dplR/ChangeLog 2013-10-24 17:00:50 UTC (rev 701)
@@ -1,11 +1,30 @@
* CHANGES IN dplR VERSION 1.5.7
+File: DESCRIPTION
+-----------------
+
+- Import gmp (>= 0.5-2)
+
+File: NAMESPACE
+---------------
+
+- New imports from gmp and utils
+- Export redfit() and runcrit()
+
Various .R files
----------------
- Check that length of vector does not overflow integer datatype
before use of .C()
+File: common.interval.R
+-----------------------
+
+- Optimizations (for example, less subsetting of the 'rwl' data.frame)
+- Better handling of corner cases (zero dimensions etc.)
+- In the plot (make.plot = TRUE), length of lines was adjusted:
+ First year a, last year b is 'b - a + 1' years, not 'b - a' years
+
File: corr.rwl.seg.R
--------------------
@@ -30,7 +49,13 @@
every round of a loop, so no need to reinitialize
- Braces always used in if (else) constructs
+Files: redfit.R, redfit.c
+-------------------------
+- New function redfit() based on REDFIT by Schulz and Mudelsee. Also
+ another exported function runcrit().
+
+
* CHANGES IN dplR VERSION 1.5.6
File: write.tucson.R
Modified: pkg/dplR/DESCRIPTION
===================================================================
--- pkg/dplR/DESCRIPTION 2013-10-24 15:59:10 UTC (rev 700)
+++ pkg/dplR/DESCRIPTION 2013-10-24 17:00:50 UTC (rev 701)
@@ -3,19 +3,23 @@
Type: Package
Title: Dendrochronology Program Library in R
Version: 1.5.7
-Date: 2013-03-19
+Date: 2013-10-12
Authors at R: c(person(c("Andrew", "G."), "Bunn", role = c("aut", "cph",
"cre", "trl"), email = "andrew.bunn at wwu.edu"), person("Mikko",
- "Korpela", role = c("aut", "cph")), person("Franco", "Biondi",
+ "Korpela", role = c("aut", "trl")), person("Franco", "Biondi",
role = c("aut", "cph")), person("Filipe", "Campelo", role =
c("aut", "cph")), person("Pierre", "Mérian", role = c("aut",
- "cph")), person("Fares", "Qeadan", role = c("aut", "cph")),
- person("Christian", "Zang", role = c("aut", "cph")))
-Author: Andy Bunn [aut, cph, cre, trl], Mikko Korpela [aut, cph], Franco Biondi [aut, cph], Filipe Campelo [aut, cph], Pierre Mérian [aut, cph], Fares Qeadan [aut, cph], Christian Zang [aut, cph]
+ "cph")), person("Manfred", "Mudelsee", role = "aut"),
+ person("Fares", "Qeadan", role = c("aut", "cph")),
+ person("Michael", "Schulz", role = "aut"), person("Christian",
+ "Zang", role = c("aut", "cph")))
+Author: Andy Bunn [aut, cph, cre, trl], Mikko Korpela [aut, trl], Franco Biondi [aut, cph], Filipe Campelo [aut, cph], Pierre Mérian [aut, cph], Manfred Mudelsee [aut], Fares Qeadan [aut, cph], Michael Schulz [aut], Christian Zang [aut, cph]
+Copyright: Authors and Aalto University (for work of M. Korpela)
Maintainer: Andy Bunn <andy.bunn at wwu.edu>
Depends: R (>= 2.15.0)
-Imports: graphics, grDevices, grid, stats, utils, digest (>= 0.2.3),
- lattice (>= 0.13-6), stringr (>= 0.4), XML (>= 2.1-0)
+Imports: gmp (>= 0.5-2), graphics, grDevices, grid, stats, utils,
+ digest (>= 0.2.3), lattice (>= 0.13-6), stringr (>= 0.4), XML
+ (>= 2.1-0)
Suggests: foreach, iterators, RUnit (>= 0.4.25)
Description: This package contains functions for performing tree-ring
analyses, IO, and graphics.
Modified: pkg/dplR/NAMESPACE
===================================================================
--- pkg/dplR/NAMESPACE 2013-10-24 15:59:10 UTC (rev 700)
+++ pkg/dplR/NAMESPACE 2013-10-24 17:00:50 UTC (rev 701)
@@ -1,11 +1,14 @@
-useDynLib(dplR, dplR.gini=gini, dplR.mean=exactmean,
- dplR.rcompact=rcompact, dplR.sens1=sens1, dplR.sens2=sens2,
- dplR.tbrm=tbrm, rwl.readloop=readloop)
+useDynLib(dplR, dplR.gini=gini, dplR.makear1=makear1,
+ dplR.mean=exactmean, dplR.rcompact=rcompact,
+ dplR.seg50=seg50, dplR.sens1=sens1, dplR.sens2=sens2,
+ dplR.spectr=spectr, dplR.tbrm=tbrm, rwl.readloop=readloop)
import(graphics, stats)
importFrom(digest, digest)
+importFrom(gmp, as.bigq, as.bigz, chooseZ, is.bigq)
+
importFrom(grDevices, rainbow)
importFrom(grid, gpar, grid.lines, grid.newpage, grid.polygon,
@@ -17,7 +20,8 @@
importFrom(stringr, str_pad, str_trim)
-importFrom(utils, head, installed.packages, read.fwf, tail)
+importFrom(utils, head, installed.packages, read.fwf, tail,
+ packageVersion, write.table)
importFrom(XML, xmlEventParse)
@@ -25,10 +29,13 @@
combine.rwl, common.interval, corr.rwl.seg, corr.series.seg,
crn.plot, detrend, detrend.series, ffcsaps, fill.internal.NA,
gini.coef, glk, hanning, i.detrend, i.detrend.series, morlet,
- po.to.wc, pointer, powt, rcs, read.compact, read.crn, read.fh,
- read.ids, read.rwl, read.tridas, read.tucson, rwi.stats,
- rwi.stats.legacy, rwi.stats.running, rwl.stats, sea, seg.plot,
- sens1, sens2, series.rwl.plot, skel.plot, spag.plot, strip.rwl,
- tbrm, tridas.vocabulary, uuid.gen, wavelet.plot, wc.to.po,
+ po.to.wc, pointer, powt, print.redfit, rcs, read.compact,
+ read.crn, read.fh, read.ids, read.rwl, read.tridas,
+ read.tucson, redfit, runcrit, rwi.stats, rwi.stats.legacy,
+ rwi.stats.running, rwl.stats, sea, seg.plot, sens1, sens2,
+ series.rwl.plot, skel.plot, spag.plot, strip.rwl, tbrm,
+ tridas.vocabulary, uuid.gen, wavelet.plot, wc.to.po,
write.compact, write.crn, write.rwl, write.tridas,
write.tucson)
+
+S3method(print, redfit)
Modified: pkg/dplR/R/common.interval.R
===================================================================
--- pkg/dplR/R/common.interval.R 2013-10-24 15:59:10 UTC (rev 700)
+++ pkg/dplR/R/common.interval.R 2013-10-24 17:00:50 UTC (rev 701)
@@ -8,39 +8,64 @@
if (!all(vapply(rwl, is.numeric, FALSE, USE.NAMES=FALSE))) {
stop("'rwl' must have numeric columns")
}
+ rnames <- row.names(rwl)
+ if (is.null(rnames)) {
+ stop("'rwl' must have row names")
+ }
+ yrs <- as.numeric(rnames)
+ if (!is.numeric(yrs) || any(is.na(yrs)) || any(round(yrs) != yrs)) {
+ stop("row names of 'rwl' must be interpretable as years")
+ }
check.flags(make.plot)
type2 <- match.arg(type, c("series", "years", "both"))
## rm.short is a function to remove short series and keep the
- ## series with overlap
- rm.short <- function(rwl, flag=FALSE) {
+ ## series with overlaps
+ rm.short <- function(rwl, yrs, rwlNotNA, row.idx, flag=FALSE) {
n <- 0
- rwl <- rwl[!vapply(rwl, function(x) all(is.na(x)), TRUE)]
- series.range <- vapply(rwl, yr.range, numeric(2),
- yr = as.numeric(row.names(rwl)))
+ anyNotNA <- apply(rwlNotNA, 2, any)
+ which.good <- which(anyNotNA)
+ nCol.orig <- length(which.good)
+ series.range <- matrix(NA_real_, 2, nCol.orig)
+ for (k in seq_len(nCol.orig)) {
+ series.range[, k] <- yr.range(rwl[[which.good[k]]][row.idx],
+ yr.vec = yrs)
+ }
+ span.order <-
+ which.good[sort.list(series.range[2, ] - series.range[1, ])]
+ nRow.orig <- nrow(rwlNotNA)
+ keep.col <- logical(length(rwl))
+ keep.col[which.good] <- TRUE
+ keep.col.output <- keep.col
+ dontkeep.row <- rep.int(TRUE, nRow.orig)
+ keep.row.output <- rep.int(FALSE, nRow.orig)
+ nRow <- 0
+ nRow.output <- 0
+ nCol.output <- nCol.orig
+ nCol <- nCol.orig
- span.order <- order(series.range[2, ] - series.range[1, ])
- to.keep <- rep(TRUE, length(span.order))
-
- rwl.output <- rwl
-
- for (i in seq(0, max(0, length(span.order) - 2))) {
- if(i > 0) {
- to.keep[span.order[i]] <- FALSE
+ for (i in seq(0, max(0, nCol.orig - 2))) {
+ if (i > 0) {
+ keep.col[span.order[i]] <- FALSE
+ nCol <- nCol - 1
+ if (nCol * nRow.orig < n) {
+ ## to break if it is not possible to improve the
+ ## common interval
+ break
+ }
}
- rwl.short <- rwl[to.keep]
- if (ncol(rwl.short) * nrow(rwl.short) < n) {
- ## to break if it is not possible to improve the
- ## common interval
- break
- }
- rwl.short <- na.omit(rwl.short)
- n.years <- ncol(rwl.short) * nrow(rwl.short)
+ tmp <- apply(rwlNotNA[dontkeep.row, keep.col, drop = FALSE], 1, all)
+ dontkeep.row[dontkeep.row] <- !tmp
+ nRow <- nRow + sum(tmp)
+ n.years <- nCol * nRow
## to keep the rwl if has more years
if (n.years > n) {
n <- n.years
- rwl.output <- rwl.short
+ keep.col.output <- keep.col
+ keep.row.output <- !dontkeep.row
+ nCol.output <- nCol
+ nRow.output <- nRow
if (flag) {
## to give the common interval with the highest
## sample depth for the case of
@@ -49,101 +74,179 @@
}
}
}
- rwl.output
+ list(nRow.output, nCol.output, keep.row.output, keep.col.output)
}
###########
- rwl.orig <- rwl
- yrs <- as.numeric(row.names(rwl))
+ nCol.rwl <- length(rwl)
+ nRow.rwl <- nrow(rwl)
+ yrs.ordered <- all(diff(yrs) >= 0)
+ if (!yrs.ordered) {
+ order.yrs <- sort.list(yrs)
+ }
output <- 0
opt <- 0
- rwl.output <- as.data.frame(matrix(0, 0, 0))
+ keep.row.output <- numeric(0)
+ keep.col.output <- logical(nCol.rwl)
+ nCol.output <- 0
+ nRow.output <- 0
+ nCol <- 0
+ nRow <- 0
+ rwlNotNA <- !is.na(rwl)
## to get sample depth
- if (ncol(rwl) > 0) {
- tmp <- rowSums(!is.na(rwl))
+ if (nCol.rwl > 0) {
+ samp.depth <- rowSums(rwlNotNA)
} else {
- tmp <- rep(0, nrow(rwl)) # R bug number 14959
+ ## Workaround for R bug number 14959. Fixed in R >= 2.15.2.
+ samp.depth <- 0
}
- for (i in dec(max(tmp), 2)) { # dec() forces a decreasing sequence
- tmp[tmp > i] <- i
- common.range <- range(as.integer(names(tmp)[tmp %in% i]))
- rwl.common <- subset(rwl,
- yrs >= common.range[1] & yrs <= common.range[2])
- if (i * nrow(rwl.common) < output){
+ type.series <- type2 == "series"
+ type.years <- type2 == "years"
+ for (i in dec(max(samp.depth), 2)) { # dec() forces a decreasing sequence
+ if (yrs.ordered) {
+ tmp <- which(samp.depth >= i)
+ row.idx <- tmp[1]:tmp[length(tmp)]
+ } else {
+ common.range <- range(yrs[samp.depth >= i])
+ row.idx <- which(yrs >= common.range[1] & yrs <= common.range[2])
+ }
+ nRow <- length(row.idx)
+ if (i * nRow < output) {
break
}
- if (type2 == "series") {
- rwl.output <- rm.short(rwl.common, flag=TRUE)
+ if (type.series) {
+ tmp <- rm.short(rwl, yrs[row.idx],
+ rwlNotNA[row.idx, , drop = FALSE], row.idx,
+ flag = TRUE)
+ nRow.output <- tmp[[1]]
+ nCol.output <- tmp[[2]]
+ keep.row.output <- row.idx[tmp[[3]]]
+ keep.col.output <- tmp[[4]]
break
- } else if (type2 == "years") {
- rwl.common <- rm.short(rwl.common)
- opt <- ncol(rwl.common) * nrow(rwl.common)
- } else if (type2 == "both") {
- rwl.common <- rwl.common[!vapply(rwl.common,
- function(x) any(is.na(x)),
- TRUE)]
- opt <- ncol(rwl.common) * nrow(rwl.common)
+ } else if (type.years) {
+ tmp <- rm.short(rwl, yrs[row.idx],
+ rwlNotNA[row.idx, , drop = FALSE], row.idx)
+ nRow <- tmp[[1]]
+ nCol <- tmp[[2]]
+ keep.row <- tmp[[3]]
+ keep.col <- tmp[[4]]
+ } else { # type2 == "both"
+ keep.col <- apply(rwlNotNA[row.idx, , drop = FALSE], 2, all)
+ nCol <- sum(keep.col)
}
- if(opt > output) {
+ opt <- nRow * nCol
+ if (opt > output) {
output <- opt
- rwl.output <- rwl.common
+ nRow.output <- nRow
+ nCol.output <- nCol
+ if (type.years) {
+ keep.row.output <- row.idx[keep.row]
+ } else {
+ keep.row.output <- row.idx
+ }
+ keep.col.output <- keep.col
}
}
if (make.plot) {
- ## original rwl
- series.range <- vapply(rwl.orig, yr.range, numeric(2),
- yr = as.numeric(row.names(rwl)))
- ## ensure that series.range is a matrix
- dim(series.range) <- c(2, length(rwl))
- first.year <- series.range[1, ]
- yr <- as.numeric(row.names(rwl.orig))
+ op <- par(no.readonly = TRUE)
+ on.exit(par(op))
+ par(mar = c(5, 5, 2, 2) + 0.1, mgp = c(1.25, 0.25, 0), tcl = 0.25)
+ if (nRow.rwl > 0 && nCol.rwl > 0) {
+ ## original rwl
+ series.range <- vapply(rwl, yr.range, numeric(2), yr.vec = yrs)
+ ## ensure that series.range is a matrix
+ dim(series.range) <- c(2, length(rwl))
+ first.year <- series.range[1, ]
- neworder <- order(first.year, decreasing = FALSE)
- segs <- rwl.orig[neworder]
- n.col <- ncol(segs)
- seq.col <- seq_len(n.col)
- for (i in seq.col) {
- segs[[i]][!is.na(segs[[i]])] <- i
- }
-
- ## common.rwl
- yr2 <- as.numeric(row.names(rwl.output))
- segs2 <- segs
- for (j in seq_len(ncol(segs2))) {
- if (names(segs)[j] %in% colnames(rwl.output)) {
- ## get correct vector
- segs2[!(yr %in% yr2), j] <- NA
+ neworder <- sort.list(first.year, na.last = TRUE)
+ rwl.first <- first.year[neworder[1]]
+ if (is.na(rwl.first)) {
+ if (yrs.ordered) {
+ rwl.first <- yrs[1]
+ rwl.last <- yrs[nRow.rwl]
+ } else {
+ rwl.first <- min(yrs)
+ rwl.last <- max(yrs)
+ }
} else {
- segs2[, j] <- NA
+ rwl.last <- max(series.range[2, ], na.rm = TRUE)
}
+ plot(1, 1, type = "n", xlim = c(rwl.first, rwl.last + 1),
+ ylim = c(1, nCol.rwl), axes = FALSE, ylab = "",
+ xlab = gettext("Year", domain = "R-dplR"))
+ rwl.seq <- seq(from = rwl.first, to = rwl.last + 1, by = 0.5)
+ n.rwl.seq <- length(rwl.seq)
+ rwl.everyother <- seq(from = 2, by = 2, length.out = nRow.rwl)
+ } else {
+ plot(1, 1, type = "n", axes = FALSE, ylab = "", xlab = "")
}
-
sub.str1 <- gettextf("Original: %d series, %d years",
- ncol(rwl.orig), nrow(rwl.orig), domain="R-dplR")
+ nCol.rwl, nRow.rwl, domain="R-dplR")
sub.str2 <-
gettextf("Common Interval (type='%s'): %d series x %d years = %d",
- type2, ncol(rwl.output), nrow(rwl.output),
- ncol(rwl.output) * nrow(rwl.output), domain="R-dplR")
- sub.str <- paste(sub.str1, sub.str2, sep='\n')
- op <- par(no.readonly = TRUE)
- on.exit(par(op))
- par(mar = c(5, 5, 2, 2) + 0.1, mgp = c(1.25, 0.25, 0), tcl = 0.25)
- plot(yr, segs[[1]], type = "n", ylim = c(1, n.col), axes = FALSE,
- ylab = "", xlab = gettext("Year", domain = "R-dplR"))
- mtext(text=sub.str, side=1, line=3)
- apply(segs, 2, lines, x = yr, lwd = 2, col="grey")
- apply(segs2, 2, lines, x = yr, lwd = 2, col="black")
- axis(2, at = seq.col, labels = names(segs), srt = 45, tick = FALSE,
- las = 2)
- axis(1)
- range.output <- range(as.numeric(rownames(rwl.output)))
- abline(v=range.output, lty="dashed")
- axis(3, at=range.output, labels=range.output, tcl=-0.25)
+ type2, nCol.output, nRow.output,
+ nCol.output * nRow.output, domain="R-dplR")
+ sub.str <- paste(sub.str1, sub.str2, sep="\n")
+ mtext(text = sub.str, side = 1, line = 3)
+ ## common.rwl
+ yrs2 <- yrs[keep.row.output]
+ any.common <- length(yrs2) > 0
+ if (any.common) {
+ common.first <- min(yrs2)
+ common.last <- max(yrs2)
+ common.seq <- seq(from = common.first,
+ to = common.last + 1, by = 0.5)
+ n.common.seq <- length(common.seq)
+ common.everyother <- seq(from = 2, by = 2, length.out = nRow.output)
+ }
+ if (!yrs.ordered) {
+ order.yrs <- sort.list(yrs)
+ order.yrs2 <- sort.list(yrs2)
+ }
+ for (i in seq_len(nCol.rwl)) {
+ this.col <- neworder[i]
+ seg <- rwl[[this.col]]
+ seg[rwlNotNA[, this.col]] <- i
+ if (yrs.ordered) {
+ seg.ordered <- seg
+ } else {
+ seg.ordered <- seg[order.yrs]
+ }
+ seg.fill <- rep.int(i, n.rwl.seq)
+ seg.fill[rwl.everyother] <- seg.ordered
+ lines(rwl.seq, seg.fill, lwd = 2, col = "grey")
+ if (keep.col.output[this.col]) {
+ seg2 <- seg[keep.row.output]
+ if (!yrs.ordered) {
+ seg2 <- seg2[order.yrs2]
+ }
+ seg2.fill <- rep.int(i, n.common.seq)
+ seg2.fill[common.everyother] <- seg2
+ lines(common.seq, seg2.fill, lwd = 2, col = "black")
+ }
+ }
+ if (nCol.rwl > 0) {
+ axis(2, at = seq_len(nCol.rwl), labels = names(rwl)[neworder],
+ srt = 45, tick = FALSE, las = 2)
+ }
+ if (nRow.rwl > 0) {
+ axis(1)
+ }
+ if (any.common) {
+ common.at <- c(common.first, common.last + 1)
+ common.labels <- as.character(c(common.first, common.last))
+ abline(v = common.at, lty = "dashed")
+ axis(3, at = common.at, labels = common.labels, tcl = -0.25)
+ }
box()
}
- rwl.output
+ if (nRow.output < nRow.rwl || nCol.output < nCol.rwl) {
+ rwl[keep.row.output, keep.col.output, drop = FALSE]
+ } else {
+ rwl
+ }
}
Modified: pkg/dplR/R/detrend.R
===================================================================
--- pkg/dplR/R/detrend.R 2013-10-24 15:59:10 UTC (rev 700)
+++ pkg/dplR/R/detrend.R 2013-10-24 17:00:50 UTC (rev 701)
@@ -14,24 +14,29 @@
if(!make.plot &&
("Spline" %in% method2 || "ModNegExp" %in% method2) &&
!inherits(try(suppressWarnings(req.it <-
- require(iterators, quietly=TRUE)),
+ requireNamespace("iterators",
+ quietly=TRUE)),
silent = TRUE),
"try-error") && req.it &&
!inherits(try(suppressWarnings(req.fe <-
- require(foreach, quietly=TRUE)),
+ requireNamespace("foreach",
+ quietly=TRUE)),
silent = TRUE),
"try-error") && req.fe){
- it.rwl <- iter(rwl, by = "col")
+ it.rwl <- iterators::iter(rwl, by = "col")
## a way to get rid of "no visible binding" NOTE in R CMD check
rwl.i <- NULL
- out <- foreach(rwl.i=it.rwl, .packages="dplR") %dopar% {
- fits <- detrend.series(rwl.i, make.plot=FALSE,
- method=method2, nyrs=nyrs, f=f,
- pos.slope=pos.slope)
- if(is.data.frame(fits))
- row.names(fits) <- rn
- fits
- }
+ out <- foreach::"%dopar%"(foreach::foreach(rwl.i=it.rwl,
+ .packages="dplR"),
+ {
+ fits <- detrend.series(rwl.i, make.plot=FALSE,
+ method=method2,
+ nyrs=nyrs, f=f,
+ pos.slope=pos.slope)
+ if(is.data.frame(fits))
+ row.names(fits) <- rn
+ fits
+ })
} else{
out <- list()
for(i in seq_len(ncol(rwl))){
Modified: pkg/dplR/R/detrend.series.R
===================================================================
--- pkg/dplR/R/detrend.series.R 2013-10-24 15:59:10 UTC (rev 700)
+++ pkg/dplR/R/detrend.series.R 2013-10-24 17:00:50 UTC (rev 701)
@@ -3,6 +3,8 @@
method = c("Spline", "ModNegExp", "Mean"),
nyrs = NULL, f = 0.5, pos.slope = FALSE)
{
+ stopifnot(identical(make.plot, TRUE) || identical(make.plot, FALSE),
+ identical(pos.slope, FALSE) || identical(pos.slope, TRUE))
known.methods <- c("Spline", "ModNegExp", "Mean")
method2 <- match.arg(arg = method,
choices = known.methods,
@@ -37,11 +39,14 @@
ModNegExp <- try(nec.func(y2), silent=TRUE)
if(class(ModNegExp)=="try-error") {
## Straight line via linear regression
- tm <- seq_along(y2)
- lm1 <- lm(y2 ~ tm)
- ModNegExp <- predict(lm1)
- if(coef(lm1)[2] > 0 && !pos.slope)
+ tm <- cbind(1, seq_along(y2))
+ lm1 <- lm.fit(tm, y2)
+ coefs <- lm1[["coefficients"]]
+ if (all(is.finite(coefs)) && (coefs[2] <= 0 || pos.slope)) {
+ ModNegExp <- drop(tm %*% coefs)
+ } else {
ModNegExp <- rep(mean(y2), length(y2))
+ }
}
resids$ModNegExp <- y2 / ModNegExp
do.mne <- TRUE
Modified: pkg/dplR/R/pointer.R
===================================================================
--- pkg/dplR/R/pointer.R 2013-10-24 15:59:10 UTC (rev 700)
+++ pkg/dplR/R/pointer.R 2013-10-24 17:00:50 UTC (rev 701)
@@ -26,7 +26,6 @@
if (nyrs < 2) {
stop("'rwl' must have at least 2 rows")
}
- nseries <- ncol(rwl2)
gv <- rwl2[-1, , drop=FALSE] / rwl2[-nyrs, , drop=FALSE]
out <- matrix(NA_real_, nrow=nyrs - 1, ncol=7)
colnames(out) <- c("Year", "Nb.series", "Perc.pos", "Perc.neg",
Modified: pkg/dplR/R/powt.R
===================================================================
--- pkg/dplR/R/powt.R 2013-10-24 15:59:10 UTC (rev 700)
+++ pkg/dplR/R/powt.R 2013-10-24 17:00:50 UTC (rev 701)
@@ -28,8 +28,8 @@
runn.M <- (drop.1 + drop.n) / 2
runn.S <- abs(drop.1 - drop.n)
runn.S[runn.S == 0] <- prec # add minimal difference
- mod <- lm(log(runn.S) ~ log(runn.M))
- b <- coef(mod)[2]
+ mod <- lm.fit(cbind(1, log(runn.M)), log(runn.S))
+ b <- mod[["coefficients"]][2]
1 - b
}
transf <- function(x) {
Copied: pkg/dplR/R/redfit.R (from rev 700, branches/redfit/R/redfit.R)
===================================================================
--- pkg/dplR/R/redfit.R (rev 0)
+++ pkg/dplR/R/redfit.R 2013-10-24 17:00:50 UTC (rev 701)
@@ -0,0 +1,1529 @@
+### This part of dplR was (arguably non-trivially) translated and
+### adapted from public domain Fortran program REDFIT, version 3.8e
+### (Michael Schulz and Manfred Mudelsee). The possibly non-free parts
+### of REDFIT derived from Numerical Recipes were not used.
+### http://www.geo.uni-bremen.de/geomod/staff/mschulz/
+### Author of the dplR version is Mikko Korpela.
+###
+### Copyright (C) 2013 Aalto University
+###
+### This program is free software; you can redistribute it and/or modify
+### it under the terms of the GNU General Public License as published by
+### the Free Software Foundation; either version 2 of the License, or
+### (at your option) any later version.
+###
+### This program is distributed in the hope that it will be useful,
+### but WITHOUT ANY WARRANTY; without even the implied warranty of
+### MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+### GNU General Public License for more details.
+###
+### A copy of the GNU General Public License is available at
+### http://www.r-project.org/Licenses/
+
+## Comments have mostly been copied verbatim from the original version
+## (a few typos were fixed). New comments are prefixed with "dplR:".
+
+## Estimate red-noise background of an autospectrum, which is estimated from
+## an unevenly spaced time series. In addition, the program corrects for the
+## bias of Lomb-Scargle Fourier transform (correlation of Fourier components),
+## which depends on the distribution of the sampling times t(i) along the
+## time axis.
+##
+## Main Assumptions:
+## -----------------
+## - The noise background can be approximated by an AR(1) process.
+## - The distribution of data points along the time axis is not
+## too clustered.
+## - The potential effect of harmonic signal components on the
+## estimation procedure is neglected.
+## - The time series has to be weakly stationary.
+##
+## The first-order autoregressive model, AR(1) model, which is used
+## to describe the noise background in a time series x(t_i), reads
+##
+##
+## x(i) = rho(i) * x(i-1) + eps(i) (1)
+##
+##
+## with t(i) - t(i-1)
+## rho(i) = exp(- -------------)
+## tau
+##
+## and eps ~ NV(0, vareps). To ensure Var[red] = 1, we set
+##
+## 2 * (t(i) - t(i-1))
+## vareps = 1 - exp(- -------------------).
+## tau
+##
+## Stationarity of the generated AR(1) time series is assured by dropping
+## the first N generated points.
+##
+##
+## Computational Steps:
+## --------------------
+##
+## 1. Estimate autospectrum Gxx of the unevenly spaced input time series in the
+## interval [0,fNyq], using the Lomb-Scargle Fourier transform in combination
+## with the Welch-Overlapped-Segment-Averaging (WOSA) procudure, as described
+## in Schulz and Stattegger (1997).
+##
+## 2. Estimate tau from the unevenly sampled time series using the time-
+## domain algorithm of Mudelsee (200?).
+##
+## 3. Determine the area under Gxx -> estimator of data variance ==> varx.
+##
+## 4. Repeat Nsim times
+## - create AR(1) time series (red) acc. to Eq. 1, using the observation
+## times of the input data, but each time with different eps(i)
+## - estimate autospectrum of red ==> Grr
+## - scale Grr such that area under the spectrum is identical to varx
+## - sum Grr ==> GrrSum
+##
+## 5. Determine arithmetic mean of GrrSum ==> GrrAvg.
+##
+## 6. Ensure that area under GrrAvg is identical to varx (account for rounding
+## errors).
+##
+## 7. Calculate theoretical AR(1) spectrum for the estimated tau ==> GRedth.
+##
+## 8. Scale GRedth such that area under the spectrum is identical to varx (this
+## step is required since the true noise variance of the data set is
+## unknown).
+##
+## 9. Estimate the frequency-dependent correction factor (corr) for the
+## Lomb-Scargle FT from the ratio between mean of the estimated AR(1) spectra
+## (GrrAvg) and the scaled theoretical AR(1) spectrum (GRedth).
+##
+## 10. Use correction factors to eliminate the bias in the estimated spectrum
+## Gxx ==> Gxxc.
+##
+## 11. Scale theorectical AR(1) spectrum for various significance levels.
+
+## Notes:
+## ------
+## * A linear trend is subtracted from each WOSA segment.
+##
+## * tau is estimated separately for each WOSA segment and subsequently
+## averaged.
+##
+## * Default max. frequency = avg. Nyquist freq. (hifac = 1.0).
+##
+## dplR note: Authors of REDFIT
+## Authors: Michael Schulz, MARUM and Faculty of Geosciences, Univ. Bremen
+## -------- Klagenfurter Str., D-28334 Bremen
+## mschulz at marum.de
+## www.geo.uni-bremen.de/~mschulz
+##
+## Manfred Mudelsee, Inst. of Meteorology, Univ. Leipzig
+## Stephanstr. 3, D-04103 Leipzig
+## Mudelsee at rz.uni-leipzig.de
+##
+## Reference: Schulz, M. and Mudelsee, M. (2002) REDFIT: Estimating
+## ---------- red-noise spectra directly from unevenly spaced paleoclimatic
+## time series. Computers and Geosciences, 28, 421-426.
+
+
+redfit <- function(x, t, tType = c("time", "age"), nsim = 1000, mctest = TRUE,
+ ofac = 4, hifac = 1, n50 = 3, rhopre = NULL,
+ p = c(0.10, 0.05, 0.02), iwin = 2,
+ txOrdered = FALSE, verbose = FALSE, seed = NULL,
+ maxTime = 10, nLimit = 10000) {
+ cl <- match.call()
+ if (!is.null(seed)) {
+ set.seed(seed)
+ }
+ MIN_POINTS <- 2
+ WIN_NAMES <- c("rectangular", "welch i", "hanning",
+ "triangular", "blackman-harris")
+ ## dplR: 21 is the lower limit of nsim where !anyDuplicated(c(idx80,
+ ## idx90, idx95, idx99)) is TRUE. (Also, none of the indices is
+ ## 0.) For more reliable results, a much greated value is
+ ## recommended.
+ NSIM_LIMIT <- 21
+ ## dplR: Check
+ tType2 <- match.arg(tType)
+ tTime <- tType2 == "time"
+ stopifnot(is.numeric(x))
+ if (!is.null(rhopre)) {
+ stopifnot(is.numeric(rhopre), length(rhopre) == 1, rhopre <= 1)
+ }
+ stopifnot(is.numeric(ofac), length(ofac) == 1, is.finite(ofac))
+ if (ofac < 1) {
+ stop("oversampling factor 'ofac' must be >= 1")
+ }
+ stopifnot(is.numeric(hifac), length(hifac) == 1, is.finite(hifac))
+ if (hifac <= 0) {
+ stop("'hifac' must be positive")
+ }
+ stopifnot(is.numeric(n50), length(n50) == 1, is.finite(n50), n50 >= 1,
+ round(n50) == n50)
+ stopifnot(is.numeric(nsim), length(nsim) == 1, is.finite(nsim), nsim >= 1,
+ round(nsim) == nsim)
+ if (length(p) > 0) {
+ stopifnot(is.numeric(p) || is.bigq(p), p > 0, p < 1)
+ }
+ stopifnot(is.numeric(maxTime), length(maxTime) == 1, maxTime >= 0)
+ stopifnot(is.numeric(nLimit), length(nLimit) == 1, nLimit >= 0,
+ round(nLimit) == nLimit)
+ stopifnot(identical(txOrdered, TRUE) || identical(txOrdered, FALSE))
+ stopifnot(identical(verbose, TRUE) || identical(verbose, FALSE))
+ stopifnot(identical(mctest, TRUE) || identical(mctest, FALSE))
+ if (mctest && nsim < NSIM_LIMIT) {
+ stop(gettextf("if 'mctest' is TRUE, 'nsim' must be at least %.0f",
+ NSIM_LIMIT, domain = "R-dplR"),
+ domain = NA)
+ }
+ ## dplR: iwin can be a number or a string. iwin2 is a number %in% 0:4
+ if (is.numeric(iwin)) {
+ if (length(iwin) != 1 || !(iwin %in% 0:4)) {
+ stop("numeric 'iwin' must be 0, 1, 2, 3 or 4")
+ }
+ iwin2 <- iwin
+ } else if (is.character(iwin)) {
+ iwin2 <- match.arg(tolower(iwin), WIN_NAMES)
+ winvec <- 0:4
+ names(winvec) <- WIN_NAMES
+ iwin2 <- winvec[iwin2]
+ } else {
+ stop("'iwin' must be numeric or character")
+ }
+ if (is.double(x)) {
+ x2 <- x
+ } else {
+ x2 <- as.numeric(x)
+ }
+ np <- as.numeric(length(x2))
+ tGiven <- !missing(t)
+ if (tGiven) {
+ if (is.double(t)) {
+ t2 <- t
+ } else {
+ t2 <- as.numeric(t)
+ }
+ if (length(t2) != np) {
+ stop("lengths of 't' and 'x' must match")
+ }
+ } else {
+ t2 <- as.numeric(seq_len(np))
+ }
+ naidx <- is.na(x2)
+ if (tGiven) {
+ naidx <- naidx | is.na(t2)
+ }
+ if (any(naidx)) {
+ goodidx <- which(!naidx)
+ t2 <- t2[goodidx]
+ x2 <- x2[goodidx]
+ nporig <- np
+ np <- as.numeric(length(x2))
+ nna <- nporig - np
+ warning(sprintf(ngettext(nna,
+ "%.0f NA value removed",
+ "%.0f NA values removed",
+ domain = "R-dplR"), nna), domain = NA)
+ }
+ if (np < MIN_POINTS) {
+ stop(gettextf("too few points (%.0f), at least %.0f needed",
+ np, MIN_POINTS, domain = "R-dplR"), domain = NA)
+ }
+ duplT <- FALSE
+ if (tGiven && !txOrdered) {
+ idx <- order(t2)
+ t2 <- t2[idx]
+ x2 <- x2[idx]
+ dupl <- duplicated(t2)
+ if (any(dupl)) {
+ duplT <- TRUE
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/dplr -r 701
More information about the Dplr-commits
mailing list