Edit models, working shiny app

R/combmodels.R

@@ -138,7 +138,7 @@ for (clus in clusters) {
     clusdf <- filter(aggdf, cluster == clus) %>% 
         dplyr::select(read_time, kwh = kwh_tot_mean) %>% 
         left_join(mtempdf, by = "read_time") %>% cbind(harm.y, harm.w, harm.d)
-    largm <- lm(nform.comp, data = clusdf, model = FALSE, qr = FALSE)
+    largm <- lm(nform.comp, data = clusdf, model = FALSE, qr = TRUE)
     allmods[[clus]] <- largm
 }
 

R/weathmod.R

@@ -80,3 +80,6 @@ ctmplot
 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)
+
+saveRDS(hmod, "../models/weatherMinharmonic.rds")
+saveRDS(maxhmod, "../models/weatherMaxharmonic.rds")

shiny/app.R

@@ -0,0 +1,135 @@
+library(reticulate)
+library(shiny)
+library(dplyr)
+library(ggplot2)
+library(shinycssloaders)
+library(TSA)
+theme_set(theme_bw())
+
+p <- import("pandas")
+aggdf <- p$read_pickle("../data/9-clusters.agg.pkl")
+aggdf$cluster <- factor(aggdf$cluster)
+sns <- import("seaborn")
+
+models <- readRDS('../models/1kmods.rds')
+clusters = names(models)
+cbp <- as.character(p$Series(sns$color_palette("colorblind", as.integer(length(clusters)))$as_hex()))
+names(cbp) <- clusters
+minhwm <- readRDS("../models/weatherMinharmonic.rds")
+maxhwm <- readRDS("../models/weatherMaxharmonic.rds")
+monthchoices <- 1:12
+names(monthchoices) <- month.name
+weekdaychoices = 1:7
+names(weekdaychoices) = c("Sunday", "Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday")
+
+
+ui <- navbarPage("Counties Power demand modelling", collapsible = TRUE, selected = "Prediction",
+      tabPanel("Data and model",
+               plotOutput("dataPlot") %>% withSpinner(type = 5),
+               fluidRow(
+                       column(4,
+               selectInput("selDR", "Cluster", choices = clusters, selected = clusters[1])
+               ),
+                       column(4,
+               dateRangeInput("dpDR", "Date range", start = min(as.Date(aggdf$read_time)),
+                              end = max(as.Date(aggdf$read_time)),
+                              min = min(as.Date(aggdf$read_time)),
+                              max = max(as.Date(aggdf$read_time)),
+                              format = "dd/mm/yyyy")
+               ),
+                       column(4, p(paste("Point represents mean value of cluster, coloured line represents model prediction of",
+                                          "mean value of cluster."))
+               )
+               )),
+      tabPanel("Prediction",
+               plotOutput("predPlot") %>% withSpinner(type = 5),
+               fluidRow(
+                        column(8,
+               sliderInput("temprange", "Temperature range", min = 0, max = 30, value = c(10, 18), width = "100%", post = " °C",
+                           step = 0.5, dragRange = TRUE)),
+               column(2,
+                      selectInput("predmon", "Month (approximate)", monthchoices, selected = 1, width = "100%")
+                      ),
+               column(2,
+                      selectInput("predday", "Day of week", weekdaychoices, selected = 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%"))
+                )
+               
+               )
+      )
+
+
+server <- function(input, output) {
+    output$dataPlot <- renderPlot({
+        clus <- input$selDR
+        fildf <- filter(aggdf, cluster %in% clus) %>% 
+            mutate(r = resid(models[[clus]]), f = fitted(models[[clus]])) %>%
+            filter(read_time >= as.POSIXct(input$dpDR[1], tz = "UTC"), read_time <= as.POSIXct(input$dpDR[2], tz = "UTC"))
+
+        yrange = range(c(range(fildf$f), range(fildf$kwh_tot_mean)))
+        dplot <- ggplot(fildf, aes(x = read_time, y = f, color = cluster))
+
+        dplot <- dplot + geom_line(size = 2)
+
+        dplot <- dplot + geom_point(aes(y = kwh_tot_mean), color = "black")
+
+        dplot <- dplot + coord_cartesian(xlim = as.POSIXct(input$dpDR, tz = "UTC"), 
+                                         ylim = c(yrange[1] - (yrange[2] - yrange[1]) * 0.05,
+                                                  yrange[2] + (yrange[2] - yrange[1]) * 0.05),
+                                         expand = FALSE)
+
+        dplot <- dplot + scale_color_manual(values = cbp)
+        dplot <- dplot + theme(legend.position = "none")
+
+        dplot <- dplot + labs(y = "kwh", x = "Date/Time", caption = sprintf("Cluster: %s", clus))
+
+        return(dplot)
+    })
+
+    prediction <- reactive({
+        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.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))
+        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)
+        }
+        return(data.frame(x = 1:48, y = predvec))
+    })
+
+    output$predPlot <- renderPlot({
+        predf <- prediction()
+        ggplot(predf, aes(x, y)) + geom_line()
+    })
+
+}
+
+
+shinyApp(ui = ui, server = server)