Permalink
Browse files

Move hmatr() to its own file. Add plot routine.

  • Loading branch information...
1 parent c6860cb commit e17e51584e5788796b90451967a4b3fa55927798 @asl committed Apr 23, 2012
Showing with 57 additions and 26 deletions.
  1. +1 −0 NAMESPACE
  2. +56 −0 R/hmatr.R
  3. +0 −26 R/ssa.R
View
@@ -84,3 +84,4 @@ S3method("rforecast", "1d-ssa")
S3method("vforecast", "1d-ssa")
S3method("roots", "lrf")
S3method("plot", "lrf")
+S3method("plot", "hmatr")
View
@@ -0,0 +1,56 @@
+# R package for Singular Spectrum Analysis
+# Copyright (c) 2012 Anton Korobeynikov <asl@math.spbu.ru>
+#
+# 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.
+#
+# You should have received a copy of the GNU General Public
+# License along with this program; if not, write to the
+# Free Software Foundation, Inc., 675 Mass Ave, Cambridge,
+# MA 02139, USA.
+
+hmatr <- function(F, ...,
+ B = N %/% 4, T = N %/% 4, L = B %/% 2,
+ neig = 10) {
+ N <- length(F)
+
+ # Pre-calculate embedding vectors and their squared norms
+ th <- t(hankel(F, L = L))
+ cth2 <- c(0, cumsum(rowSums(th^2)))
+ cth2 <- (cth2[1:(N-T)+(T-L+1)] - cth2[1:(N-T)])
+
+ hc <- function(idx) {
+ Fb <- F[idx:(idx+B)] # Form a basis subspace
+ s <- new.ssa(Fb, L = L, ..., neig = min(2*neig, 50))
+
+ # Calculate the distance
+ U <- s$U[, 1:neig, drop = FALSE]
+ # FIXME: Can we use FFT stuff here somehow?
+ cXU2 <- c(0, cumsum(rowSums((th %*% U)^2)))
+
+ 1 - (cXU2[1:(N-T)+(T-L+1)] - cXU2[1:(N-T)]) / cth2
+ }
+
+ h <- sapply(1:(N-B), hc)
+ class(h) <- "hmatr"
+
+ invisible(h)
+}
+
+plot.hmatr <- function(x,
+ col = rev(heat.colors(256)),
+ main = "Heterogeneity Matrix",
+ xlab = "", ylab = "",
+ ...) {
+ image(1:nrow(x), 1:ncol(x), x,
+ col = col, xlab = xlab, ylab = ylab, main = main, ...)
+}
+
View
@@ -225,32 +225,6 @@ clusterify.ssa <- function(this, groups, nclust = length(groups) / 2,
out;
}
-hmatr <- function(F, ...,
- B = N %/% 4, T = N %/% 4, L = B %/% 2,
- neig = 10) {
- N <- length(F)
-
- # Pre-calculate embedding vectors and their squared norms
- th <- t(hankel(F, L = L))
- cth2 <- c(0, cumsum(rowSums(th^2)))
- cth2 <- (cth2[1:(N-T)+(T-L+1)] - cth2[1:(N-T)])
-
- hc <- function(idx) {
- Fb <- F[idx:(idx+B)] # Form a basis subspace
- s <- new.ssa(Fb, L = L, ..., neig = min(2*neig, 50))
-
- # Calculate the distance
- U <- s$U[, 1:neig, drop = FALSE]
- # FIXME: Can we use FFT stuff here somehow?
- cXU2 <- c(0, cumsum(rowSums((th %*% U)^2)))
-
- 1 - (cXU2[1:(N-T)+(T-L+1)] - cXU2[1:(N-T)]) / cth2
- }
-
- invisible(sapply(1:(N-B), hc))
-}
-
-
'$.ssa' <- function(this, name) {
if (ind <- charmatch(name, names(this), nomatch = 0))
return (this[[ind]]);

0 comments on commit e17e515

Please sign in to comment.