Fourier stuff

R/clusterviz.R

@@ -7,11 +7,28 @@ library(forecast)
 theme_set(theme_bw())
 use_virtualenv("../venv/")
 
+fextract <- function(x, y, keep = 1, top = TRUE) {
+    sdy <- sd(y)
+    my <- mean(y)
+    stany <- (y - my) / sdy
+    ftf <- fft(stany)
+    if (top) {
+        ftf[rank(-abs(ftf)) > keep] <- 0
+    } else {
+        ftf[(keep + 1):length(ftf)] <- 0
+    }
+    rfv <- Re(fft(ftf, inverse = TRUE))
+    dfr <- data.frame(x = x, y = y, f = (rfv - mean(rfv)) / sd(rfv) * sdy + my)
+    dfr$res <- dfr$y - dfr$f
+    return(dfr)
+}
+
+
 p <- import("pandas")
 sns <- import("seaborn")
 cbp <- as.character(p$Series(sns$color_palette("colorblind", as.integer(9))$as_hex()))
 aggdf <- p$read_pickle("../data/9-clusters.agg.pkl")
-aggdf <- as.data.frame(aggdf)
+# aggdf <- as.data.frame(aggdf)
 aggdf$cluster <- factor(aggdf$cluster)
 str(aggdf)
 clusters = levels(aggdf$cluster)
@@ -146,3 +163,33 @@ c9ts <- filter(aggdf, cluster == "9")$kwh_tot_mean
 ctsnp <- msts(c9ts, c(48, 48*7))
 ctbats <- tbats(ctsnp)
 plot(forecast(ctbats, h = 48 * 7 * 4))
+
+p <- periodogram(c1ts)
+dd <- data.frame(freq = p$freq, spec = p$spec) %>% mutate(per = 1/freq)
+dd %>% arrange(desc(spec)) %>% mutate(d = per / 48) %>% head()
+
+c9ts <- filter(aggdf, cluster == "9")
+
+ggplot(c9ts, aes(x = read_time, y = kwh_tot_mean)) + geom_line()
+
+nft <- fextract(c9ts$read_time, c9ts$kwh_tot_mean, keep = 10)
+ggplot(nft, aes(x, y)) + geom_line() + 
+    geom_line(aes(x, f), color = "blue") +
+    scale_x_datetime(date_breaks = "1 day", date_labels = "%a, %-d %B %Y") +
+    coord_cartesian(xlim = c(as.POSIXct("2017-07-16", tz = "UTC"), as.POSIXct("2017-07-23", tz = "UTC")), expand = TRUE)
+
+clus <- "9"
+kp <- 50
+cts <- filter(aggdf, cluster == clus)
+# ggplot(cts, aes(x = read_time, y = kwh_tot_mean)) + geom_line()
+nft <- fextract(cts$read_time, cts$kwh_tot_mean, keep = kp)
+ggplot(nft, aes(x, y)) + geom_line() + 
+    geom_line(aes(x, f), color = "blue") +
+    scale_x_datetime(date_breaks = "2 days", date_labels = "%a, %-d %b %Y") +
+    coord_cartesian(xlim = c(as.POSIXct("2017-09-14", tz = "UTC"), as.POSIXct("2017-09-28", tz = "UTC")), expand = TRUE)
+
+ggplot(nft, aes(x, res)) + geom_line() + scale_x_datetime(date_breaks = "1 week") +
+    coord_cartesian(xlim = c(as.POSIXct("2017-08-01", tz = "UTC"), as.POSIXct("2017-09-28", tz = "UTC")), expand = TRUE)
+
+ggplot(nft, aes(y = res)) + geom_boxplot()
+ggplot(nft, aes(x = res)) + geom_histogram()

R/fftest.R

@@ -0,0 +1,29 @@
+# Testing the reconstruction of arbitrary sequencies using fourier transforms
+library(ggplot2)
+library(scales)
+library(dplyr)
+theme_set(theme_bw())
+
+fextract <- function(x, y, keep = 1, top = TRUE) {
+    sdy <- sd(y)
+    my <- mean(y)
+    stany <- (y - my) / sdy
+    ftf <- fft(stany)
+    if (top) {
+        ftf[rank(-abs(ftf)) > keep] <- 0
+    } else {
+        ftf[(keep + 1):length(ftf)] <- 0
+    }
+    rfv <- Re(fft(ftf, inverse = TRUE))
+    return(data.frame(x = x, y = y, f = (rfv - mean(rfv)) / sd(rfv) * sdy + my))
+}
+
+n <- 100
+k <- 1
+dat1 <- data.frame(x = 1:n, yo = 10 * sin(rescale(1:n, 2 * c(-pi, pi) + rnorm(2))))
+dat1$y <- dat1$yo + rnorm(n)
+fnew <- fextract(dat1$x, dat1$y, keep = k, top = TRUE)
+# ggplot(rdat1, aes(x, y)) + geom_line()
+ggplot(fnew, aes(x, y)) + geom_line() + geom_point() + 
+    geom_line(aes(x, f), color = "blue") + geom_point(aes(x, f), color = "blue")
+