[Vegan-commits] r1509 - in branches/1.17: R inst

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

Author: jarioksa
Date: 2011-02-25 15:10:45 +0100 (Fri, 25 Feb 2011)
New Revision: 1509

Modified:
   branches/1.17/R/decorana.R
   branches/1.17/inst/ChangeLog
Log:
merge r1500,4: handling zero eigenvalues in decorana

Modified: branches/1.17/R/decorana.R
===================================================================
--- branches/1.17/R/decorana.R	2011-02-24 14:27:01 UTC (rev 1508)
+++ branches/1.17/R/decorana.R	2011-02-25 14:10:45 UTC (rev 1509)
@@ -5,6 +5,7 @@
     Const1 <- 1e-10
     Const2 <- 5
     Const3 <- 1e-11
+    ZEROEIG <- 1e-7 # consider as zero eigenvalue
     veg <- as.matrix(veg)
     if (any(rowSums(veg) <= 0)) 
         stop("All row sums must be >0 in the community matrix: remove empty sites.")
@@ -63,6 +64,11 @@
                    ix2 = as.integer(ix2), ix3 = as.integer(ix3), aidot = as.double(aidot), 
                    adotj = as.double(adotj), PACKAGE = "vegan")[c("x", "y", 
                                              "eig")]
+    ## Eigenvalues can be computed zeros: zero the results exactly to
+    ## avoid problems with rescaling and estimating final eigenvalues
+    ## in arbitrary corner cases.
+    if (s1$eig < ZEROEIG)
+        s1$x[] <- s1$y[] <- s1$eig <- 0
     if (!ira) 
         ix1 <- .Fortran("cutup", x = as.double(s1$x), ix = as.integer(ix1), 
                         mi = as.integer(cep$mi), mk = as.integer(mk), PACKAGE = "vegan")$ix
@@ -78,6 +84,8 @@
                    ix2 = as.integer(ix2), ix3 = as.integer(ix3), aidot = as.double(aidot), 
                    adotj = as.double(adotj), PACKAGE = "vegan")[c("x", "y", 
                                              "eig")]
+    if (s2$eig < ZEROEIG)
+        s2$x[] <- s2$y[] <- s2$eig <- 0
     if (!ira) 
         ix2 <- .Fortran("cutup", x = as.double(s2$x), ix = as.integer(ix2), 
                         mi = as.integer(cep$mi), mk = as.integer(mk), PACKAGE = "vegan")$ix
@@ -93,6 +101,8 @@
                    ix2 = as.integer(ix2), ix3 = as.integer(ix3), aidot = as.double(aidot), 
                    adotj = as.double(adotj), PACKAGE = "vegan")[c("x", "y", 
                                              "eig")]
+    if (s3$eig < ZEROEIG)
+        s3$x[] <- s3$y[] <- s3$eig <- 0
     if (!ira) 
         ix3 <- .Fortran("cutup", x = as.double(s3$x), ix = as.integer(ix3), 
                         mi = as.integer(cep$mi), mk = as.integer(mk), PACKAGE = "vegan")$ix
@@ -108,19 +118,25 @@
                    ix2 = as.integer(ix2), ix3 = as.integer(ix3), aidot = as.double(aidot), 
                    adotj = as.double(adotj), PACKAGE = "vegan")[c("x", "y", 
                                              "eig")]
-    s1$x <- .Fortran("yxmult", y = as.double(s1$y), x = as.double(s1$x), 
+    if (s4$eig < ZEROEIG)
+        s4$x[] <- s4$y[] <- s4$eig <- 0
+    if (s1$eig > ZEROEIG)
+        s1$x <- .Fortran("yxmult", y = as.double(s1$y), x = as.double(s1$x), 
                      as.integer(cep$mi), as.integer(cep$n), as.integer(cep$nid), 
                      as.integer(cep$ibegin), as.integer(cep$iend), as.integer(cep$idat), 
                      as.double(cep$qidat), PACKAGE = "vegan")$x/aidot
-    s2$x <- .Fortran("yxmult", y = as.double(s2$y), x = as.double(s2$x), 
+    if (s2$eig > ZEROEIG)
+        s2$x <- .Fortran("yxmult", y = as.double(s2$y), x = as.double(s2$x), 
                      as.integer(cep$mi), as.integer(cep$n), as.integer(cep$nid), 
                      as.integer(cep$ibegin), as.integer(cep$iend), as.integer(cep$idat), 
                      as.double(cep$qidat), PACKAGE = "vegan")$x/aidot
-    s3$x <- .Fortran("yxmult", y = as.double(s3$y), x = as.double(s3$x), 
+    if (s3$eig > ZEROEIG)
+        s3$x <- .Fortran("yxmult", y = as.double(s3$y), x = as.double(s3$x), 
                      as.integer(cep$mi), as.integer(cep$n), as.integer(cep$nid), 
                      as.integer(cep$ibegin), as.integer(cep$iend), as.integer(cep$idat), 
                      as.double(cep$qidat), PACKAGE = "vegan")$x/aidot
-    s4$x <- .Fortran("yxmult", y = as.double(s4$y), x = as.double(s1$x), 
+    if (s4$eig > ZEROEIG)
+        s4$x <- .Fortran("yxmult", y = as.double(s4$y), x = as.double(s4$x), 
                      as.integer(cep$mi), as.integer(cep$n), as.integer(cep$nid), 
                      as.integer(cep$ibegin), as.integer(cep$iend), as.integer(cep$idat), 
                      as.double(cep$qidat), PACKAGE = "vegan")$x/aidot
@@ -146,6 +162,8 @@
         var.c <- cov.wt(cproj, adotj)
         var.c <- diag(var.c$cov) * (1 - sum(var.c$wt^2))
         evals <- var.r/var.c
+        if (any(ze <- evals.decorana < ZEROEIG))
+            evals[ze] <- 0
     }
     CA <- list(rproj = rproj, cproj = cproj, evals = evals, evals.decorana = evals.decorana, 
                origin = origin, v = v, fraction = v.fraction, adotj = adotj, 

Modified: branches/1.17/inst/ChangeLog
===================================================================
--- branches/1.17/inst/ChangeLog	2011-02-24 14:27:01 UTC (rev 1508)
+++ branches/1.17/inst/ChangeLog	2011-02-25 14:10:45 UTC (rev 1509)
@@ -12,6 +12,8 @@
 
 	* merged r1501: as.preston 'tiesplit = TRUE' bug.
 
+	* merged r1500,4: decorana zero-eigenvalue fix.
+
 Version 1.17-7 (released February 16, 2011)
 
 	* merged r1489: nestednodf with weighted = FALSE.