[Dplr-commits] r744 - in pkg/dplR: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Mar 27 12:35:13 CET 2014
Author: mvkorpel
Date: 2014-03-27 12:35:13 +0100 (Thu, 27 Mar 2014)
New Revision: 744
Modified:
pkg/dplR/ChangeLog
pkg/dplR/DESCRIPTION
pkg/dplR/R/rwi.stats.running.R
Log:
rwi.stats.running.R:
* 'n.trees' and 'n.cores' should now be OK
* 'c.eff' is now 0 if no correlations were computed
* When using period = "common", 'n' is now 0 instead of the full
number in case there are no complete cases in the running window.
* Technical optimizations
Modified: pkg/dplR/ChangeLog
===================================================================
--- pkg/dplR/ChangeLog 2014-03-27 00:42:17 UTC (rev 743)
+++ pkg/dplR/ChangeLog 2014-03-27 11:35:13 UTC (rev 744)
@@ -5,8 +5,17 @@
- Added signal-to-noise ratio as an output. Followed pg 109 in
Cook's chapter in Hughes et al. 2011.
+- New outputs 'n.cores' and 'n.trees' show the total number of cores and
+ trees in the window, respectively. At least one non-missing value is
+ required for a core and tree to be counted.
+- When using period = "common" in rwi.stats() or
+ rwi.stats.running(), the number of trees 'n' in the return value
+ is now 0 instead of the full number in case there are no complete
+ cases in the running window.
+- 'c.eff' in the return value is now 0 if no correlations were
+ computed
+- Optimizations
-
* CHANGES IN dplR VERSION 1.5.9
Files: dplR.h, rcompact.c, redfit.c
Modified: pkg/dplR/DESCRIPTION
===================================================================
--- pkg/dplR/DESCRIPTION 2014-03-27 00:42:17 UTC (rev 743)
+++ pkg/dplR/DESCRIPTION 2014-03-27 11:35:13 UTC (rev 744)
@@ -3,7 +3,7 @@
Type: Package
Title: Dendrochronology Program Library in R
Version: 1.6.0
-Date: 2014-03-26
+Date: 2014-03-27
Authors at R: c(person("Andy", "Bunn", role = c("aut", "cph",
"cre", "trl"), email = "andy.bunn at wwu.edu"), person("Mikko",
"Korpela", role = c("aut", "trl")), person("Franco", "Biondi",
Modified: pkg/dplR/R/rwi.stats.running.R
===================================================================
--- pkg/dplR/R/rwi.stats.running.R 2014-03-27 00:42:17 UTC (rev 743)
+++ pkg/dplR/R/rwi.stats.running.R 2014-03-27 11:35:13 UTC (rev 744)
@@ -137,13 +137,15 @@
rwi3 <- rwi2
}
}
+ rwiNotNA <- !is.na(rwi3)
n.years <- nrow(rwi3)
if (running.window && window.length > n.years) {
stop("'window.length' is larger than the number of years in 'rwi'")
}
- unique.trees <- unique(ids3$tree)
+ treeIds <- ids3$tree
+ unique.trees <- unique(treeIds)
n.trees <- length(unique.trees)
if (n.trees < 2) {
stop("at least 2 trees are needed")
@@ -151,7 +153,7 @@
cores.of.tree <- list()
seq.tree <- seq_len(n.trees)
for (i in seq.tree) {
- cores.of.tree[[i]] <- which(ids3$tree==unique.trees[i])
+ cores.of.tree[[i]] <- which(treeIds==unique.trees[i])
}
## n.trees.by.year is recorded before setting rows with missing
@@ -159,21 +161,22 @@
tree.any <- matrix(FALSE, n.years, n.trees)
for (i in seq.tree) {
tree.any[, i] <-
- apply(!is.na(rwi3[, ids3$tree == unique.trees[i], drop=FALSE]),
- 1, any)
+ apply(rwiNotNA[, treeIds == unique.trees[i], drop=FALSE], 1, any)
}
n.trees.by.year <- rowSums(tree.any)
- good.rows <- which(n.trees.by.year > 1)
## Easy way to force complete overlap of data
if (period2 == "common") {
- bad.rows <- which(apply(is.na(rwi3), 1, any))
+ bad.rows <- !apply(rwiNotNA, 1, all)
rwi3[bad.rows, ] <- NA
- good.rows <- setdiff(good.rows, bad.rows)
+ rwiNotNA[bad.rows, ] <- FALSE
+ good.rows.flag <- !bad.rows
period.common <- TRUE
} else {
+ good.rows.flag <- n.trees.by.year > 1
period.common <- FALSE
}
+ good.rows <- which(good.rows.flag)
if (length(good.rows) < min.corr.overlap) {
stop("too few years with enough trees for correlation calculations")
@@ -205,8 +208,8 @@
(n.years - offset - window.length) %/% window.advance
max.idx <-
offset + window.length + n.windows.minusone * window.advance
- n.data[i] <- sum(!is.na(rwi3[intersect(good.rows,
- (1 + offset):max.idx), ]))
+ rowIdx <- seq(1 + offset, max.idx)
+ n.data[i] <- sum(rwiNotNA[rowIdx[good.rows.flag[rowIdx]], ])
}
## In case of a tie, choose large offset.
## In practice, this prefers recent years.
@@ -294,27 +297,45 @@
rbar.bt <- rsum.bt / n.bt
}
+ coresPresent <-
+ which(apply(rwiNotNA[year.idx, , drop = FALSE], 2, any))
+ treesPresent <- unique(treeIds[coresPresent])
+ nCores <- length(coresPresent)
+ nTrees <- length(treesPresent)
if (period.common) {
## If period is "common", we are only looking at the rows
- ## with no missing values.
- n <- n.trees
+ ## with no missing values (if any, so all or nothing).
+ n <- nTrees
} else {
## Number of trees averaged over the years in the window.
## We keep this number separate of the correlation
## estimates, i.e. the data from some tree / year may
## contribute to n without taking part in the correlation
## estimates.
- n <- mean(n.trees.by.year[year.idx], na.rm=TRUE)
+ n <- mean(n.trees.by.year[year.idx])
}
## Expressed population signal
if (n.wt == 0) {
- c.eff <- 1
+ if (n.bt > 0) {
+ c.eff <- 1
+ } else {
+ c.eff <- 0
+ }
rbar.eff <- rbar.bt
} else {
- c.eff.rproc <- mean(1 / n.cores.tree, na.rm=TRUE)
- c.eff <- 1 / c.eff.rproc # bookkeeping only
- rbar.eff <- rbar.bt / (rbar.wt + (1 - rbar.wt) * c.eff.rproc)
+ nCoresTree <- na.omit(n.cores.tree)
+ uniqueNC <- unique(nCoresTree)
+ ## The branches are equivalent but optimized for numerical
+ ## precision in each situation
+ if (length(uniqueNC) == 1) {
+ c.eff <- uniqueNC
+ rbar.eff <- rbar.bt / (rbar.wt + (1 - rbar.wt) / c.eff)
+ } else {
+ c.eff.rproc <- mean(1 / nCoresTree)
+ c.eff <- 1 / c.eff.rproc # bookkeeping only
+ rbar.eff <- rbar.bt / (rbar.wt + (1 - rbar.wt) * c.eff.rproc)
+ }
}
## EPS is on page 146 of C&K.
## In our interpretation of EPS, we use the average number of trees.
@@ -324,17 +345,16 @@
snr <- n * rbar.eff / (1-rbar.eff)
if (running.window) {
- c(start.year = start.year, mid.year = mid.year, end.year = end.year,
- n.cores = n.cores, n.trees = n.trees, n = n,
- n.tot = n.tot, n.wt = n.wt, n.bt = n.bt, rbar.tot = rbar.tot,
- rbar.wt = rbar.wt, rbar.bt = rbar.bt, c.eff = c.eff,
- rbar.eff = rbar.eff, eps = eps, snr = snr)
+ out <- c(start.year = start.year,
+ mid.year = mid.year, end.year = end.year)
} else {
- c(n.cores = n.cores, n.trees = n.trees, n = n,
- n.tot = n.tot, n.wt = n.wt, n.bt = n.bt, rbar.tot = rbar.tot,
- rbar.wt = rbar.wt, rbar.bt = rbar.bt, c.eff = c.eff,
- rbar.eff = rbar.eff, eps = eps, snr = snr)
+ out <- numeric(0)
}
+ c(out,
+ n.cores = nCores, n.trees = nTrees, n = n,
+ n.tot = n.tot, n.wt = n.wt, n.bt = n.bt, rbar.tot = rbar.tot,
+ rbar.wt = rbar.wt, rbar.bt = rbar.bt, c.eff = c.eff,
+ rbar.eff = rbar.eff, eps = eps, snr = snr)
}
## Iterate over all windows
More information about the Dplr-commits
mailing list