Fix some stuff with the shiny app so that it actually works

R/weathmod.R

     geom_line(aes(y = runmin), color = "blue") +
     geom_line(aes(y = runmax), color = "red")
 
-tplot
+# tplot
 
-tplot + coord_cartesian(xlim = c(as.POSIXct("2016-12-01"), as.POSIXct("2017-03-01")))
+# tplot + coord_cartesian(xlim = c(as.POSIXct("2016-12-01"), as.POSIXct("2017-03-01")))
 
 # Create a harmonic (sine wave) model for minimum temperature
-yharm <- harmonic(ts(1:nrow(fulltemp), frequency = floor(365.25 * 48)), 2)
+yharm <- harmonic(ts(1:nrow(fulltemp), frequency = floor(365.25 * 48)), 2) %>% as.data.frame()
 #dharm <- harmonic(ts(1:nrow(fulltemp), frequency = floor(48)), 1)
-hmod <- lm(fulltemp$runmin ~ yharm)
+hmod <- lm(fulltemp$runmin ~ ., data = yharm)
 summary(hmod)
 
 hmdf <- data.frame(x = fulltemp$temp_timestamp, y = fulltemp$runmin, f = fitted(hmod), r = resid(hmod))
 tmplot <- ggplot(hmdf, aes(x = x, y = y)) + geom_line(aes(y = f), color = "blue", size = 2) + geom_point() +
     geom_point(aes(y = r), color = "darkgreen")
 
-tmplot
+# tmplot
 
-tmplot + coord_cartesian(xlim = c(as.POSIXct("2017-05-01", tz = "UTC"), as.POSIXct("2017-06-01", tz = "UTC")))
+# tmplot + coord_cartesian(xlim = c(as.POSIXct("2017-05-01", tz = "UTC"), as.POSIXct("2017-06-01", tz = "UTC")))
 
 
 
-maxhmod <- lm(fulltemp$runmax ~ yharm)
+maxhmod <- lm(fulltemp$runmax ~ ., data = yharm)
 summary(maxhmod)
 mhmdf <- data.frame(x = fulltemp$temp_timestamp, y = fulltemp$runmax, f = fitted(maxhmod), r = resid(maxhmod))
 
 mtmplot <- ggplot(mhmdf, aes(x = x, y = y)) + geom_line(aes(y = f), color = "blue", size = 2) + geom_point() +
     geom_point(aes(y = r), color = "darkgreen")
 
-mtmplot
+# mtmplot
 
-mtmplot + coord_cartesian(xlim = c(as.POSIXct("2017-05-01", tz = "UTC"), as.POSIXct("2017-06-01", tz = "UTC")))
+# mtmplot + coord_cartesian(xlim = c(as.POSIXct("2017-05-01", tz = "UTC"), as.POSIXct("2017-06-01", tz = "UTC")))
 
 hmdf.comb <- left_join(hmdf, mhmdf, by = "x", suffix = c(".min", ".max"))
 
     geom_point(aes(y = y.max), color = "magenta") +
     geom_point(aes(y = y.min), color = "lightblue")
 
-ctmplot
+# ctmplot
 
-ctmplot + coord_cartesian(xlim = c(as.POSIXct("2017-05-01", tz = "UTC"), as.POSIXct("2017-06-01", tz = "UTC")))
+# ctmplot + coord_cartesian(xlim = c(as.POSIXct("2017-05-01", tz = "UTC"), as.POSIXct("2017-06-01", tz = "UTC")))
 
 write.csv(hmdf.comb, "../data/weatherharm.csv", row.names = FALSE)
 

shiny/app.R

 library(shinycssloaders)
 library(TSA)
 theme_set(theme_bw())
+use_virtualenv("../venv/")
 
 p <- import("pandas")
 aggdf <- p$read_pickle("../data/9-clusters.agg.pkl")
                )
                )),
       tabPanel("Prediction",
-               plotOutput("predPlot") %>% withSpinner(type = 5),
+               plotOutput("predPlot"),
                fluidRow(
-                        column(8,
+                        column(6,
                sliderInput("temprange", "Temperature range", min = 0, max = 30, value = c(10, 18), width = "100%", post = " °C",
                            step = 0.5, dragRange = TRUE)),
                column(2,
                       ),
                column(2,
                       selectInput("predday", "Day of week", weekdaychoices, selected = 2, width = "100%")
+                      ),
+               column(2,
+                      numericInput("prl", "Reference line (kwh)", value = 2, width = "100%")
                       )
                ),
 
                h5("Number of ICPs of each cluster serviced by node:"),
                fluidRow(
                     column(1, numericInput("c1v", "1:", value = 1, min = 0, step = 1, width = "100%"), offset = 0),
-                    column(1, numericInput("c2v", "2:", value = 1, min = 0, step = 1, width = "100%")),
-                    column(1, numericInput("c3v", "3:", value = 1, min = 0, step = 1, width = "100%")),
-                    column(1, numericInput("c4v", "4:", value = 1, min = 0, step = 1, width = "100%")),
-                    column(1, numericInput("c5v", "5:", value = 1, min = 0, step = 1, width = "100%")),
-                    column(1, numericInput("c6v", "6:", value = 1, min = 0, step = 1, width = "100%")),
-                    column(1, numericInput("c7v", "7:", value = 1, min = 0, step = 1, width = "100%")),
-                    column(1, numericInput("c8v", "8:", value = 1, min = 0, step = 1, width = "100%")),
-                    column(1, numericInput("c9v", "9:", value = 1, min = 0, step = 1, width = "100%"))
+                    column(1, numericInput("c2v", "2:", value = 0, min = 0, step = 1, width = "100%")),
+                    column(1, numericInput("c3v", "3:", value = 0, min = 0, step = 1, width = "100%")),
+                    column(1, numericInput("c4v", "4:", value = 0, min = 0, step = 1, width = "100%")),
+                    column(1, numericInput("c5v", "5:", value = 0, min = 0, step = 1, width = "100%")),
+                    column(1, numericInput("c6v", "6:", value = 0, min = 0, step = 1, width = "100%")),
+                    column(1, numericInput("c7v", "7:", value = 0, min = 0, step = 1, width = "100%")),
+                    column(1, numericInput("c8v", "8:", value = 0, min = 0, step = 1, width = "100%")),
+                    column(1, numericInput("c9v", "9:", value = 0, min = 0, step = 1, width = "100%"))
                 )
                
                )
         return(dplot)
     })
 
-    prediction <- reactive({
+    prediction <- function() {
         numeclus <- c(input$c1v,input$c2v,input$c3v,input$c4v,input$c5v,input$c6v,input$c7v,input$c8v,input$c9v)
         numeclus[is.na(numeclus)] <- 0
         ystart <- floor((as.numeric(input$predmon) - 0.5) / 12 * 365.25 * 48)
         harm.y <- ts(1:48, frequency = 365.25 * 48, start = c(1, ystart)) %>% harmonic(2)
-        harm.w <- ts(1:48, frequency = 7 * 48, start = c(1, 48 * as.numeric(input$predday))) %>% harmonic(3)
+        harm.w <- ts(1:48, frequency = 7 * 48, start = c(1, 48 * (-1 + as.numeric(input$predday)))) %>% harmonic(3)
         harm.d <- ts(1:48, frequency = 48, start = c(1, 1)) %>% harmonic(3)
         colnames(harm.y) <- sprintf("%s.%s.%s", "year", rep(c("cos", "sin"), each = ncol(harm.y)/2), rep(1:(ncol(harm.y)/2), times = 2))
         colnames(harm.w) <- sprintf("%s.%s.%s", "week", rep(c("cos", "sin"), each = ncol(harm.w)/2), rep(1:(ncol(harm.w)/2), times = 2))
         colnames(harm.d) <- sprintf("%s.%s.%s", "day",  rep(c("cos", "sin"), each = ncol(harm.d)/2), rep(1:(ncol(harm.d)/2), times = 2))
-        yharm <- harmonic(ts(1:48, frequency = floor(365.25 * 48), start = c(1, ystart)), 2)
-        harnames <- paste0("yharm", colnames(yharm))
-        # yharm <- as.data.frame(yharm)
-        colnames(yharm) <- harnames
-        wmmin <- input$temprange[1] - predict(minhwm, as.data.frame(yharm))
-        wmmax <- input$temprange[2] - predict(maxhwm, as.data.frame(yharm))
+        qharm <- harmonic(ts(1:48, frequency = floor(365.25 * 48), start = c(1, ystart)), 2)
+        # harnames <- paste0("yharm", colnames(qharm))
+        qharm <- as.data.frame(qharm)
+        # names(qharm) <- harnames
+        wmmin <- input$temprange[1] - predict(minhwm, newdata = qharm, type = 'response')
+        wmmax <- input$temprange[2] - predict(maxhwm, newdata = qharm, type = 'response')
         predinp <- data.frame(resmin = wmmin, resmax = wmmax)
         predinp <- cbind(predinp, harm.y, harm.w, harm.d)
         predvec = rep(0, 48)
         for (c in 1:length(clusters)) {
-            predvec <- predvec + numeclus[c] * predict(models[[c]], predinp)
+            predvec <- predvec + numeclus[c] * predict(models[[c]], newdata = predinp, type = 'response')
         }
         return(data.frame(x = 1:48, y = predvec))
-    })
+    }
 
     output$predPlot <- renderPlot({
         predf <- prediction()
-        ggplot(predf, aes(x, y)) + geom_line()
+        predp <- ggplot(predf, aes(x, y)) + geom_line()
+        predp <- predp + scale_y_continuous("Predicted kwh", limits = c(0, NA))
+        predp <- predp + scale_x_continuous("Time", labels = function(x) {
+                                                sprintf("%02d:%02d", (x %/% 2) %% 24, x %% 2 * 30)
+                                         }, breaks = 0:12 * 4)
+        if (!is.na(input$prl)) {
+            predp <- predp + geom_hline(yintercept = input$prl, linetype = "dashed")
+        }
+        predp
     })
 
 }