[Dplr-commits] r760 - pkg/dplR/R

[noreply at r-forge.r-project.org]

Author: mvkorpel
Date: 2014-04-02 15:45:06 +0200 (Wed, 02 Apr 2014)
New Revision: 760

Modified:
   pkg/dplR/R/normalize.xdate.R
   pkg/dplR/R/series.rho.R
Log:
Major speedup to series.rho() by use of a modified normalize.xdate()
(leave.one.out = TRUE).


Modified: pkg/dplR/R/normalize.xdate.R
===================================================================
--- pkg/dplR/R/normalize.xdate.R	2014-04-02 12:36:30 UTC (rev 759)
+++ pkg/dplR/R/normalize.xdate.R	2014-04-02 13:45:06 UTC (rev 760)
@@ -1,24 +1,59 @@
-normalize.xdate <- function(rwl, series, n, prewhiten, biweight){
+normalize.xdate <- function(rwl, series, n, prewhiten, biweight,
+                            leave.one.out = FALSE) {
+    loo <- isTRUE(leave.one.out)
     ## Run hanning filter over the data if n isn't NULL
     ## divide by mean if n is null
     if(is.null(n)){
         master.stats <- colMeans(rwl, na.rm=TRUE)
         master.df <- sweep(rwl, 2, master.stats, "/")
-        series.out <- series / mean(series, na.rm=TRUE)
+        if (!loo) {
+            series.out <- series / mean(series, na.rm=TRUE)
+        }
     } else {
         master.stats <- apply(rwl, 2, hanning, n)
         master.df <-  rwl / master.stats
-        series.out <- series / hanning(series, n)
+        if (!loo) {
+            series.out <- series / hanning(series, n)
+        }
     }
-    ## Apply ar if prewhiten
-    if(prewhiten){
-        ## drop any columns without at least four observations
-        master.df <- master.df[, colSums(!is.na(master.df)) > 3, drop=FALSE]
-        master.df <-  apply(master.df, 2, ar.func)
-        series.out <- ar.func(series.out)
+    if (loo) {
+        nseries <- ncol(rwl)
+        ## Apply ar if prewhiten
+        if(prewhiten){
+            ## mark any columns without at least four observations
+            goodCol <- colSums(!is.na(master.df)) > 3
+            series.out <-  apply(master.df, 2, ar.func)
+        } else {
+            goodCol <- rep.int(TRUE, nseries)
+            series.out <- master.df
+        }
+        master <- series.out
+        if (!biweight) {
+            for (i in seq_len(nseries)) {
+                goodCol2 <- goodCol
+                goodCol2[i] <- FALSE
+                master[, i] <-
+                    rowMeans(series.out[, goodCol2, drop = FALSE], na.rm=TRUE)
+            }
+        } else {
+            for (i in seq_len(nseries)) {
+                goodCol2 <- goodCol
+                goodCol2[i] <- FALSE
+                master[, i] <-
+                    apply(series.out[, goodCol2, drop = FALSE], 1, tbrm, C = 9)
+            }
+        }
+    } else {
+        ## Apply ar if prewhiten
+        if(prewhiten){
+            ## drop any columns without at least four observations
+            master.df <- master.df[, colSums(!is.na(master.df)) > 3, drop=FALSE]
+            master.df <-  apply(master.df, 2, ar.func)
+            series.out <- ar.func(series.out)
+        }
+
+        if (!biweight) master <- rowMeans(master.df, na.rm=TRUE)
+        else master <- apply(master.df, 1, tbrm, C = 9)
     }
-
-    if (!biweight) master <- rowMeans(master.df, na.rm=TRUE)
-    else master <- apply(master.df, 1, tbrm, C = 9)
     list(master=master, series=series.out)
 }

Modified: pkg/dplR/R/series.rho.R
===================================================================
--- pkg/dplR/R/series.rho.R	2014-04-02 12:36:30 UTC (rev 759)
+++ pkg/dplR/R/series.rho.R	2014-04-02 13:45:06 UTC (rev 760)
@@ -3,11 +3,13 @@
     rho <- numeric(nseries)
     p.val <- numeric(nseries)
     rwl.mat <- as.matrix(rwl)
+    tmp <- normalize.xdate(rwl=rwl.mat, n=n,
+                           prewhiten=prewhiten, biweight=biweight,
+                           leave.one.out = TRUE)
+    series <- tmp[["series"]]
+    master <- tmp[["master"]]
     for (i in seq_len(nseries)) {
-        tmp <- normalize.xdate(rwl=rwl.mat[, -i, drop=FALSE],
-                               series=rwl.mat[, i], n=n,
-                               prewhiten=prewhiten, biweight=biweight)
-        tmp2 <- cor.test(tmp[["series"]], tmp[["master"]],
+        tmp2 <- cor.test(series[, i], master[, i],
                          method = "spearman", alternative = "greater")
         rho[i] <- tmp2[["estimate"]]
         p.val[i] <- tmp2[["p.value"]]