[Analogue-commits] r135 - in pkg: R inst src

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

Author: gsimpson
Date: 2009-06-09 00:11:35 +0200 (Tue, 09 Jun 2009)
New Revision: 135

Modified:
   pkg/R/predict.wa.R
   pkg/R/wa.R
   pkg/inst/ChangeLog
   pkg/src/wastats.c
Log:
C code in WATpred plus changes to predict.wa and wa to make use of the new C code

Modified: pkg/R/predict.wa.R
===================================================================
--- pkg/R/predict.wa.R	2009-06-08 19:23:50 UTC (rev 134)
+++ pkg/R/predict.wa.R	2009-06-08 22:11:35 UTC (rev 135)
@@ -32,7 +32,9 @@
         want <- names(object$wa.optima)[want]
         if(object$tol.dw) {
             pred <- WATpred(newdata[,want], object$wa.optima[want],
-                            object$model.tol[want])
+                            object$model.tol[want],
+                            NROW(newdata[,want]),
+                            NCOL(newdata[,want]))
         } else {
             pred <- WApred(newdata[,want], object$wa.optima[want])
         }
@@ -46,6 +48,10 @@
             useN2 <- object$options.tol$useN2
             want <- names(object$wa.optima) %in% colnames(newdata)
             want <- names(object$wa.optima)[want]
+            nr <- NROW(X) - 1
+            nr.new <- NROW(newdata)
+            nc <- NCOL(X)
+            nc.want <- length(want)
             for(i in seq_len(n.train)) {
                 if(verbose && ((i %% 10) == 0)) {
                     cat(paste("Leave one out sample", i, "\n"))
@@ -59,9 +65,11 @@
                                 O$min.tol, O$f, ENV[-i])
                 ## CV for the training set
                 if(object$tol.dw) {
-                    wa.env <- WATpred(X[-i,], wa.optima, tol)
+                    wa.env <- WATpred(X[-i,], wa.optima, tol,
+                                      nr, nc)
                     mod.pred[i] <- WATpred(X[i,,drop=FALSE],
-                                           wa.optima, tol)
+                                           wa.optima, tol,
+                                           1, nc)
                 } else {
                     wa.env <- WApred(X[-i,], wa.optima)
                     mod.pred[i] <- WApred(X[i,,drop=FALSE],
@@ -75,7 +83,7 @@
                 ## newdata predictions
                 pred <- if(object$tol.dw) {
                     WATpred(newdata[,want], wa.optima[want],
-                            tol[want])
+                            tol[want], nr.new, nc.want)
                 } else {
                     WApred(newdata[,want], wa.optima[want])
                 }
@@ -237,10 +245,10 @@
     ((X %*% optima) / rsum)
 }
 
-WATpred <- function(X, optima, tol) {
-    ones <- rep.int(1, length(optima))
+WATpred <- function(X, optima, tol, nr, nc) {
+    #ones <- rep.int(1, length(optima))
     miss <- is.na(optima)
-    ones[miss] <- 0
+    #ones[miss] <- 0
     optima[miss] <- 0
     tol[miss] <- 1
     tol2 <- tol^2
@@ -252,5 +260,11 @@
     #                    check.margin = FALSE)[,!miss, drop = FALSE])
     #tmp / (sweep(X, 2, tol2, "/",check.margin = FALSE) %*% ones)
     ##tmp / (t(t(X) / tol2) %*% ones)
-    RowSums(t(t(X) * optima / tol2)) / (t(t(X) / tol2) %*% ones)
+    #RowSums(t(t(X) * optima / tol2)) / (t(t(X) / tol2) %*% ones)
+    res <- .C("WATpred", spp = as.double(X), opt = as.double(optima),
+              tol2 = as.double(tol2), nr = as.integer(nr),
+              nc = as.integer(nc), want = as.integer(miss),
+              stat = double(nr),
+              NAOK = TRUE, PACKAGE="analogue")$stat
+    res
 }

Modified: pkg/R/wa.R
===================================================================
--- pkg/R/wa.R	2009-06-08 19:23:50 UTC (rev 134)
+++ pkg/R/wa.R	2009-06-08 22:11:35 UTC (rev 135)
@@ -41,7 +41,7 @@
                     min.tol = min.tol, f = f, env = env)
     ## calculate WA estimate of env for each site
     if(tol.dw) {
-        wa.env <- WATpred(x, wa.optima, tol)
+        wa.env <- WATpred(x, wa.optima, tol, n.samp, n.spp)
     } else {
         wa.env <- WApred(x, wa.optima)
     }

Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog	2009-06-08 19:23:50 UTC (rev 134)
+++ pkg/inst/ChangeLog	2009-06-08 22:11:35 UTC (rev 135)
@@ -4,6 +4,9 @@
 
 	* tran: now converts input data to matrix using 'data.matrix' which
 	deals with factor variables appropriately.
+
+	* predict.wa, wa, WATpred: Now use faster C code for computing
+	predictions from WA models with tolerance down-weighting.
 	
 Version 0.6-10
 

Modified: pkg/src/wastats.c
===================================================================
--- pkg/src/wastats.c	2009-06-08 19:23:50 UTC (rev 134)
+++ pkg/src/wastats.c	2009-06-08 22:11:35 UTC (rev 135)
@@ -43,3 +43,28 @@
 	stat[i] = sqrt(sumWX / sumW);
     }
 }
+
+void WATpred(double *spp, double *opt, double *tol2, 
+	     int *nr, int *nc, int *want, double *stat)
+{
+    int i, j, ij;
+    /* spp == species data
+       opt == species optima
+       tol == species tolerances^2
+       nr, nc == n rows, n cols
+    */
+    
+    double nomin, denom;
+
+    for (i=0; i<*nr; i++) {
+	nomin = 0.0;
+	denom= 0.0;
+	for(j=0; j<*nc; j++) {
+	    ij = i + (j * *nr);
+	    nomin += (spp[ij] * opt[j]) / tol2[j];
+	    if(*want == 0)
+		denom += spp[ij] / tol2[j];
+	}
+	stat[i] = nomin / denom;
+    }
+}