ampli/shiny/app.R

library(reticulate)
library(shiny)
library(dplyr)
library(ggplot2)
library(shinycssloaders)
library(TSA)
theme_set(theme_bw())
use_virtualenv("../venv/")

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"),
               fluidRow(
                        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,
                      selectInput("predmon", "Month (approximate)", monthchoices, selected = 1, width = "100%")
                      ),
               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 = 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%"))
                )
               
               )
      )


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 <- 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 * (-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))
        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]], newdata = predinp, type = 'response')
        }
        return(data.frame(x = 1:48, y = predvec))
    }

    output$predPlot <- renderPlot({
        predf <- prediction()
        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
    })

}


shinyApp(ui = ui, server = server)