app.R

# https://stackoverflow.com/questions/45377186/increase-max-upload-file-size-in-shiny?noredirect=1&lq=1
options(shiny.maxRequestSize = 9000*1024^2)
list.of.packages <- c("tidyverse", "lubridate", "shiny", "shinythemes", "reshape2", "keras","Metrics")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

library(tidyverse)
library(lubridate)
library(shiny)
library(shinythemes)
library(reshape2)
# Sys.setenv("CUDA_VISIBLE_DEVICES" = -1)
library(keras)
library(Metrics)
# library(DT)

# https://bl.ocks.org/psychemedia/9737637
# https://stackoverflow.com/questions/29253481/data-specific-selectinput-choices-in-rmd-shiny/29255723#29255723
# https://stackoverflow.com/questions/21465411/r-shiny-passing-reactive-to-selectinput-choices?rq=1


# look up table for station ids 

look_up_table <- 
  tibble(
  names = c("CEN", "MOK", "PRE", "TSW", "WAC", "CAB", "TIH", "TIS"),
  ids = c(1,6,16,25,27,28,29,118)
  )



# Define UI for data upload app ----

ui <- fluidPage( theme = shinythemes::shinytheme("yeti"), # slate
                # shinythemes::themeSelector(),  # <--- Add this somewhere in the UI
  
  # App title ----
  titlePanel("Predictive Passenger Arrivals Models"),
  
  # Sidebar layout with input and output definitions ----
  sidebarLayout(
    # Sidebar panel for inputs ----
    sidebarPanel(
      # Input: Select a file ----
      fileInput("file1", "Choose a file",
                multiple = FALSE,
                # accept = c(""),
                # accept = c("tar.gz",
                #            " ",
                #           "text/csv",
                #            "text/comma-separated-values,text/plain",
                #            ".csv",
                #            ".tar.gz"),
                width = "80%"),
      # Horizontal line ----
      tags$hr(),
      # https://stackoverflow.com/questions/24175997/force-no-default-selection-in-selectinput
      # I() indicates it is raw JavaScript code that should be evaluated, instead
      # of a normal character string
      selectizeInput(
        'station', 'Select Station', choices = c("CEN", "MOK", "PRE", "TSW", "WAC", "CAB", "TIH", "TIS"),
        options = list(
          placeholder = 'Please select a station',
          onInitialize = I('function() { this.setValue(""); }')
        ), width="80%"
      ),
      
      # selectInput("station", "Select the station", choices = c("CEN", "MOK", "PRE", "TSW", "WAC", "CAB", "TIH", "TIS"), 
      #             selected=FALSE, width="50%", multiple=FALSE, selectize = FALSE),
      
      # tags$hr(),
      # actionButton("pri", "Print the data", width="80%"),
      tags$hr(),
      actionButton("pre", "Preprocess the data", width="80%"),
      # Horizontal line ----
      tags$hr(),
      actionButton("pred", "Make predictions", width="80%"),  
      tags$hr(),
      # Add a download button
      # downloadButton(outputId = "download_data", label = "Download predictions"),
      # downloadButton(outputId = "download_hist", label = "Download hist_avg"),
      # downloadButton(outputId = "download_y", label = "Download true demand")
      
      downloadButton(outputId = "download_outputs", label = "Download all outputs")
      
      
  ),
    
    # Main panel for displaying outputs ----
    mainPanel(
    
      plotOutput("myPlot"),
      
      verbatimTextOutput("sum"),
      tableOutput("myTable")
  
      

      
    )
    
  )
)

options(shiny.maxRequestSize = 9000*1024^2)
server <- function(input, output, session) {
  options(shiny.maxRequestSize = 9000*1024^2)
  values <- reactiveValues(df_data = NULL, station_id= NULL, station_name= NULL, station_data=NULL, processed_data=NULL,df=NULL)
  
  


  observeEvent(input$file1, {
    # https://community.rstudio.com/t/413-request-error-while-using-shinyapps-io-how-to-modify-nginx-config/31904/3
    # a = untar(input$file1$datapath, list=TRUE)
    # untar(input$file1$datapath)
    # fname = a[2]
    # values$df_data <- read.csv(fname);
    # 
    # # values$df_data <- read.csv(input$file1$datapath);
    # output$sum <- renderPrint({
    #   print(head(values$df_data, 10))
    # }) 
    # Start the clock!
    start.time <- Sys.time()
    
    observeEvent(input$file1, {
      values$df_data <- read.csv(input$file1$datapath);
      
      # Stop the clock
      end.time <- Sys.time()
      
      output$sum <- renderPrint({
        print ("Upload completed in ")
        print(end.time - start.time)
        # print(head(values$df_data, 10))
      })
    })
  }) 



  observeEvent(input$station, {
    values$station_name <- input$station;
    values$station_id <- look_up_table %>% filter(names == input$station ) %>% select(ids) %>% as.numeric();
    output$sum <- renderPrint({
      # summary(values$df_data)
      print("Station name and ID")
      print(isolate(values$station_name));
      print(isolate(values$station_id));
      # print(isolate(head(values$df_data)));
    })
  }, ignoreNULL=FALSE, ignoreInit=TRUE
    )

  
  # https://shiny.rstudio.com/articles/action-buttons.html
  # it keeps running this function whenever the selectInput is updated. That's why isolate is needed
  # observeEvent(input$pri, {
  #   values$station_data <- values$df_data %>% filter(TRAIN_ENTRY_STN ==  values$station_id & TXN_TYPE_CO == "ENT") # & TXN_TYPE_CO == "ENT"
  #   output$sum <- renderPrint({
  #     # summary(values$df_data)
  #     # print(isolate(values$station_name));
  #     # print(isolate(values$station_id));
  #     # print(isolate(head(values$station_data)));
  #     print("Finished preprocessing the data ")
  #     
  #   })
  # })
  # 

  observeEvent(input$pre, {
    values$station_data <- values$df_data %>% filter(TRAIN_ENTRY_STN ==  values$station_id & TXN_TYPE_CO == "ENT") 
    # to fill the missing values
    all_time_bins_template <- 
      tibble(
        bins = 1:96
      )
    
    # only work with 1 day 
    dd <- 
      values$station_data %>% 
      mutate(timestamp = dmy_hms(TXN_DT, tz= "Asia/Hong_Kong") ) %>% 
      mutate(dow = wday(timestamp, label=FALSE), DayOfMonth = day(timestamp)) %>% # monday is 2
      filter(!(dow %in% c(1,7))) %>% # filter to only include weekday
      group_by(DayOfMonth) %>% 
      tally(sort = T) %>% 
      slice(1) %>% 
      select(DayOfMonth) 
    
    # extract day of the week 
    dow <-
      values$station_data  %>% 
      mutate(timestamp = dmy_hms(TXN_DT, tz= "Asia/Hong_Kong") ) %>% 
      mutate(dow = wday(timestamp, label=FALSE), DayOfMonth = day(timestamp)) %>% # monday is 2
      slice(1) %>% 
      select(dow) %>% 
      as.numeric()
    
    # bin the data and counts 
    values$processed_data <- 
      values$station_data %>% 
      mutate(timestamp = dmy_hms(TXN_DT) ) %>% 
      mutate(dow = wday(timestamp, label=FALSE), DayOfMonth = day(timestamp)) %>% # monday is 2
      # only extract the time 
      mutate(time =  timestamp %>% format("%H:%M") %>% hm()) %>% 
      # bin it to be a number btw 0-96
      mutate(bin = ceiling((time %>% period_to_seconds())/900)) %>% 
      # Having intervals in a column breaks many dplyr verbs, https://github.com/tidyverse/lubridate/issues/635
      select(-time) %>% 
      # filter to only include weekday
      filter(!(dow %in% c(1,7))) %>% 
      filter(DayOfMonth == dd$DayOfMonth) %>%  
      # mutate(bin = as.numeric(cut(timestamp, breaks = "15 mins"))) %>% 
      group_by(bin, 
               add=FALSE)  %>% 
      summarise(coalesce = n())  %>% 
      right_join(all_time_bins_template, by=c("bin"= "bins")) %>%
      replace_na(list(coalesce= 1))
    
    # prepare for keras 
    
    df <- values$processed_data %>% 
      
      mutate(t1=lag(coalesce),
             t2=lag(t1),
             t3=lag(t2),
             t4=lag(t3),
             t5=lag(t4),
             t6=lag(t5),
             t7=lag(t6),
             t8=lag(t7)
      ) %>% 
      slice(-c(1:8)) 
    

    # dummy code day of week
    target_col <- (dow - 1) %>% unique()
    dow_oneHot <- data.frame(
      matrix(0, nrow(df), 5)
    )
    # dow_oneHot <- keras::to_categorical(df$dow) %>% as_tibble()
    names(dow_oneHot) <- c("M", "T", "W", "Th", "F")
    # fill the column for the test day
    dow_oneHot[,target_col] <- 1
    # keras expects categorical to start from 0
    df$dow <- dow -2
    values$df <- bind_cols(df[,-which(names(df) == "dow")], dow_oneHot)
    
    output$sum <- renderPrint({
      # print(isolate(head(values$df)));
      print("Finished preprocessing the data ")
     
 
      
    })
    
    
  })

  
  observeEvent(input$pred, {
    ##############################
    # read the model and mean/std
    ##############################
    model_path <- paste0("./models/my_model_station_", as.character(values$station_id), ".h5")
    model <- keras::load_model_hdf5(model_path)
    
    hist_avg <- readRDS(paste0("./models/station_", as.character(values$station_id),"_hist_avg.rds"))
    values$hist <-  hist_avg$demand[-c(1:5)]
    mean_std <- readRDS(paste0("./models/mean_std_", as.character(values$station_id),".rds"))
   
    training_mean <- mean_std$m
    training_std <- mean_std$std
    
    #############################
    # make predictions 
    #############################

    
 
    y.test = values$df$coalesce
    y.test = y.test[-c(1:17)]
    

    x.test <- values$df %>% select (-c( "bin", 'coalesce'))
    x.test <- as.matrix(x.test)
    # data preprocessing 
  
 
    test_data <- scale(x.test, center = training_mean, scale = training_std)
    test_targets <- y.test
    values$y <- y.test
    test_input_3D <- array_reshape(x = test_data, dim = c(dim(test_data), 1))
    
    values$results <- model %>% predict(test_input_3D) %>% as.integer()
    
    values$results <- values$results[-c(1:17)]
    
    output$sum <- renderPrint({
      # print(isolate(head(values$df)));
      # print(isolate((values$results)));
      # print(model %>% summary())
      print ("Prediction Accuracy")

    })
    
    
    output$myPlot <- renderPlot({

  
    # plot.ts(values$results,ylim=c(0,5000),col='red')
    # par(new=T)
    # plot.ts(hist_avg$demand[-c(1:8)],ylim=c(0,5000), col='springgreen4')
    # par(new=T)
    # plot.ts(y.test,ylim=c(0,5000),col='blue')
    
    
    start_day <- as.POSIXct("2010-01-01")
    start <- as.POSIXct("2010-01-01 06:00:00")
    interval <- 60
    end <- start_day + as.difftime(1, units="days")
    
    a <- seq(from=start, by=interval*120, to=end)
    x <- seq_along(values$results)
    
    plot_data <- as.data.frame(cbind(x = x,  results = values$results, expected=test_targets,
                                      historical_average=values$hist))
    
    d2 <- melt(plot_data, id="x")
    
    p <- ggplot(d2, aes(x, value, color=variable, linetype=variable))+
      geom_line() +
      # scale_linetype_manual( values = c("solid","solid","dashed")) +
      
      labs(title="Demand Prediction") +
      scale_linetype_manual(labels=c("1 Step-ahead Predictions", "Observed", "Historical Average"),
                            values = c("solid","solid","dashed"))+
      scale_color_manual(labels=c("1 Step-ahead Predictions", "Observed", "Historical Average"),
                         values = c("#d7191c","#2b83ba","#4dac26")) + 
      # scale_size_manual(labels=c("1 Step-ahead Predictions", "Observed", "Historical Average"),
      #                   values=c(1, 1, .5, .75,.75))+
      
      ylab('Passengers per 15 minutes') +
      xlab('Time of day (one unit = 15 minutes)') +
      
      theme_bw() +
      theme(plot.title=element_text(size=rel(1.5), lineheight=.9, face="bold", colour="black", hjust = .5),
            plot.subtitle=element_text(  colour="black", hjust = .5),
            axis.text=element_text(size=12), axis.title=element_text(size=14,face="bold"),
            legend.title=element_blank(), legend.text = element_text(size=14)) 
    
    # if(length(values$results)==72){
         p <- p + scale_x_continuous(breaks=seq(1,72, by = 8), labels= lapply(a[1:length(a)-1], format,'%H:%M')) + xlab('Time of day')
      p
    # }
    
    
    })
    
    output$myTable <- renderTable({
      tibble(
        RMSE_Model = Metrics::rmse(values$results, y.test),
        MAD_Model =  Metrics::mae(values$results, y.test),
        RMSE_hist = Metrics::rmse(hist_avg$demand, y.test),
        MAD_hist =  Metrics::mae(hist_avg$demand, y.test)
      )
    })
    
  })
  
  

  # Create a download handler
  output$download_outputs  <- downloadHandler(
    filename <- "all_outputs.csv",
    content = function(file) {
      data <- 
        cbind(values$results, values$y, values$hist)
      # Write the filtered data into a CSV file
      write.csv(data, file, row.names = FALSE)
    }
  )
  # output$download_data <- downloadHandler(
  #   filename <- "predictions.csv",
  #   content = function(file) {
  #     data <- 
  #       values$results
  #       # Write the filtered data into a CSV file
  #       write.csv(data, file, row.names = FALSE)
  #   }
  # )
  # 
  # # Create a download handler
  # output$download_hist <- downloadHandler(
  #   filename <- "hist_avg.csv",
  #   content = function(file) {
  #     data <- values$hist
  #       
  #     # Write the filtered data into a CSV file
  #     write.csv(data, file, row.names = FALSE)
  #   }
  # )
  # output$download_y <- downloadHandler(
  #   filename <- "observed.csv",
  #   content = function(file) {
  #     data <- values$y
  #     
  #     # Write the filtered data into a CSV file
  #     write.csv(data, file, row.names = FALSE)
  #   }
  # )
  
}

shinyApp(ui = ui, server = server)