# maps of Vienna (Figure 1) ----

# from 9 to 10 am weekdays
H.10 <- seq(from = 10, to = t.max, by = 24)[-which((1:365+1)%%7 > 4)] 
# from 5 to 6 pm weekdays
H.18 <- seq(from = 18, to = t.max, by = 24)[-which((1:365+1)%%7 > 4)] 

# get differences
D.10 <- D.18 <- matrix(NA, length(H.10), N.stations)
for (i in 1:N.stations) {
  H.10.i <- H.10[which(H.10 >= stations$ActiveFrom[i])]
  H.18.i <- H.18[which(H.18 >= stations$ActiveFrom[i])]
  
  ind.10 <- which(is.element(feeds[[i]]$t, H.10.i))
  ind.18 <- which(is.element(feeds[[i]]$t, H.18.i))
  
  D.10.i <- feeds[[i]]$C[ind.10] - feeds[[i]]$C[ind.10 - 1]
  D.18.i <- feeds[[i]]$C[ind.18] - feeds[[i]]$C[ind.18 - 1]
  
  D.10[which(is.element(H.10, H.10.i)), i] <- D.10.i
  D.18[which(is.element(H.18, H.18.i)), i] <- D.18.i
}

# specifications for plotting
D.10.mean <- colMeans(D.10, na.rm = TRUE)
D.18.mean <- colMeans(D.18, na.rm = TRUE)

color.10 <- color.18 <- rep("blue", N.stations)
color.10[which(D.10.mean < 0)] <- "red"
color.18[which(D.18.mean < 0)] <- "red"

shape.10 <- shape.18 <- rep(24, N.stations)
shape.10[which(D.10.mean < 0)] <- 25
shape.18[which(D.18.mean < 0)] <- 25

size.10 = abs(D.10.mean)
size.18 = abs(D.18.mean)

if(exists("API.key")){
  register_google(API.key)
} else {
  stop("Please supply a valid google API key!")
}

# Figure 1
Vienna.10 <- ggmap(get_map(location = c(16.3665, 48.21), zoom = 13, maptype = "roadmap")) + 
  geom_point(data = stations %>% slice(1:120), aes(x = Longitude, y = Latitude), color = color.10, size = 3*size.10, 
             shape = shape.10, fill = color.10) + 
  theme(axis.line = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.margin = unit(c(0, 0, -1, -1), 'lines')) 

pdf(file = "Plots/Descriptives/differences_vienna_10.pdf", width = 7/1.5, height = 7/1.5)
print(Vienna.10)
dev.off()
setEPS()
postscript("Plots/Descriptives/differences_vienna_10.eps", width = 7/1.5, height = 7/1.5)
print(Vienna.10)
dev.off()

Vienna.18 <- ggmap(get_map(location = c(16.3665, 48.21), zoom = 13, maptype = "roadmap")) + 
  geom_point(data = stations %>% slice(1:120), aes(x = Longitude, y = Latitude), color = color.18, size = 3*size.18, 
             shape = shape.18, fill = color.18) +
  theme(axis.line = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        plot.margin = unit(c(0, 0, -1, -1), 'lines'))

pdf(file = "Plots/Descriptives/differences_vienna_18.pdf", width = 7/1.5, height = 7/1.5)
print(Vienna.18)
dev.off()
setEPS()
postscript("Plots/Descriptives/differences_vienna_18.eps", width = 7/1.5, height = 7/1.5)
print(Vienna.18)
dev.off()


# counts by weekday (Figure 2) ----

N.weekday <- matrix(0, 24, 7)
for (s in 1:24) {
  N.weekday[s, 3] = sum(N.wo.P$dyadic[seq(0*24+s, t.max, 7*24), , ])/52 # Wed
  N.weekday[s, 4] = sum(N.wo.P$dyadic[seq(1*24+s, t.max, 7*24), , ])/52 # Thu
  N.weekday[s, 5] = sum(N.wo.P$dyadic[seq(2*24+s, t.max, 7*24), , ])/52 # Fri
  N.weekday[s, 6] = sum(N.wo.P$dyadic[seq(3*24+s, t.max, 7*24), , ])/52 # Sat
  N.weekday[s, 7] = sum(N.wo.P$dyadic[seq(4*24+s, t.max, 7*24), , ])/52 # Sun
  N.weekday[s, 1] = sum(N.wo.P$dyadic[seq(5*24+s, t.max, 7*24), , ])/51 # Mon
  N.weekday[s, 2] = sum(N.wo.P$dyadic[seq(6*24+s, t.max, 7*24), , ])/51 # Tue
}

weekday.df <- tibble(h = rep(1:24, 7), N = as.vector(N.weekday), 
                         Weekday = c(rep("Monday", 24), rep("Tuesday", 24), rep("Wednesday", 24), rep("Thursday", 24), 
                                     rep("Friday", 24), rep("Saturday", 24), rep("Sunday", 24)))
weekday.df$Weekday <- factor(weekday.df$Weekday, levels = c("Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", "Sunday"))

# Figure 2
g <- ggplot(weekday.df, aes(x = h, y = N)) + geom_line(aes(color = Weekday, linetype = Weekday)) + 
  labs(x = "h", y = "Average Trip Count per Hour", 
       title = "Average Trip Count depending on the Weekday") + 
  scale_x_continuous(breaks = 1:24) +
  theme(legend.justification = c(0.01, 0.99), legend.position = c(0.01, 0.99)) 

pdf(file = "Plots/Descriptives/N_Weekday.pdf", width = 6, height = 4)
print(g)
dev.off()
setEPS()
postscript("Plots/Descriptives/N_Weekday.eps", width = 6, height = 4)
print(g)
dev.off()


# counts by temperature (Figure 2)----

range.temp <- range(round(weather$Temperature))
bins.temp <- seq(min(range.temp), max(range.temp), 1)

# counts by temperature and different hours of the day
N.temp.8.12 <- N.temp.12.16 <- N.temp.16.20 <- rep(0, length(bins.temp))
freq.temp.8.12 <- freq.temp.12.16 <- freq.temp.16.20 <- rep(0, length(bins.temp))

ind.8.12 <- sort(c(seq(9, t.max, 24), seq(10, t.max, 24),seq(11, t.max, 24),seq(12, t.max, 24)))
ind.12.16 <- sort(c(seq(13, t.max, 24), seq(14, t.max, 24),seq(15, t.max, 24),seq(16, t.max, 24)))
ind.16.20 <- sort(c(seq(17, t.max, 24), seq(18, t.max, 24),seq(19, t.max, 24),seq(20, t.max, 24)))

# from 8 am to 12 pm
for (t in ind.8.12) {
  temp <- round(weather$Temperature[t], 1)
  N.temp.8.12[which(bins.temp == round(temp))] <- N.temp.8.12[which(bins.temp == round(temp))] + sum(N.wo.P$dyadic[t, , ])
  freq.temp.8.12[which(bins.temp == round(temp))] <- freq.temp.8.12[which(bins.temp == round(temp))] + 1
}
# from 12 pm to 4 pm
for (t in ind.12.16) {
  temp <- round(weather$Temperature[t], 1)
  N.temp.12.16[which(bins.temp == round(temp))] <- N.temp.12.16[which(bins.temp == round(temp))] + sum(N.wo.P$dyadic[t, , ])
  freq.temp.12.16[which(bins.temp == round(temp))] <- freq.temp.12.16[which(bins.temp == round(temp))] + 1
}
# from 4 pm to 8 pm
for (t in ind.16.20) {
  temp <- round(weather$Temperature[t], 1)
  N.temp.16.20[which(bins.temp == round(temp))] <- N.temp.16.20[which(bins.temp == round(temp))] + sum(N.wo.P$dyadic[t, , ])
  freq.temp.16.20[which(bins.temp == round(temp))] <- freq.temp.16.20[which(bins.temp == round(temp))] + 1
}

temp.df <- tibble(temp = rep(bins.temp, 3), 
                      N = c(N.temp.8.12, N.temp.12.16, N.temp.16.20),
                      freq = c(freq.temp.8.12, freq.temp.12.16, freq.temp.16.20),
                      time = c(rep("8am - 12pm", length(bins.temp)), rep("12pm - 4pm", length(bins.temp)), rep("4pm - 8pm", length(bins.temp))))
temp.df$time <- factor(temp.df$time, levels = c("8am - 12pm", "12pm - 4pm", "4pm - 8pm"))

# Figure 2
g <- ggplot(temp.df, aes(x = temp, y = N/freq)) + geom_line(aes(color = time, linetype = time)) + 
  labs(x = "Outside Temperature in Degrees Celsius", y = "Average Trip Count per Hour", 
       title = "Average Trip Count depending on the Outside Temperature", color = "Time Window", linetype = "Time Window") + 
  theme(legend.justification = c(0.01, 0.99), legend.position = c(0.01, 0.99)) + 
  scale_x_continuous(breaks = seq(-10, 35, 5))

pdf(file = "Plots/Descriptives/N_Temp.pdf", width = 6, height = 4)
print(g)
dev.off()
setEPS()
postscript("Plots/Descriptives/N_Temp.eps", width = 6, height = 4)
print(g)
dev.off()

# trips by distance (Figure 3) ----

bins.distance <- seq(0, 12, 0.25)
N.distance <- rep(0, 49)
freq.distance <- rep(0, 49)
for (k in 1:N.stations) {
  for (l in 1:N.stations) {
    dist <- ceiling(distances[k, l]*4)
    N.distance[dist + 1] <- N.distance[dist + 1] + sum(N.wo.P$dyadic[, k, l])
    freq.distance[dist + 1] <- freq.distance[dist + 1] + 1
  }
}

dist.df <- tibble(dist = bins.distance[2:49], N = N.distance[2:49]/freq.distance[2:49])

# Figure 3
g <- ggplot(dist.df, aes(x = dist, y = N)) + geom_bar(stat = "identity") + 
  labs(x = "Distance between two Stations in Kilometers", y = "Average Trip Count per Hour", 
       title = "Average Trip Count depending on the Distance between two Stations") + 
  scale_x_continuous(breaks = 0:12)

pdf(file = "Plots/Descriptives/N_Dist.pdf", width = 6, height = 4)
print(g)
dev.off()
setEPS()
postscript("Plots/Descriptives/N_Dist.eps", width = 6, height = 4)
print(g)
dev.off()

# function f(alpha) (Figure 3) ----

X <- seq(0, 12, 0.01)
Y <- matrix(0, length(X), 6)
for (alpha in seq(0.5, 3, 0.5)) {
  Y[, 2*alpha] <- X^alpha*exp(-X)/(alpha^alpha*exp(-alpha))
}
f.alpha.df <- tibble(x = rep(X, 6), 
                     y = as.vector(Y),
                     alpha = factor(c(rep(0.5, length(X)), rep(1, length(X)), rep(1.5, length(X)),
                                      rep(2, length(X)), rep(2.5, length(X)), rep(3, length(X)))))

# Figure 3
g <- ggplot(f.alpha.df, aes(x = x, y = y)) + geom_line(aes(color = alpha, linetype = alpha)) +
  labs(x = "x", title = expression(paste("Plot of ", tilde(f)[alpha](x), " for different ", alpha, sep = "")), 
       y = expression(tilde(f)[alpha](x)), color = expression(alpha), linetype = expression(alpha)) + 
  theme(legend.justification = c(0.99, 0.99), legend.position = c(0.99, 0.99)) + 
  scale_x_continuous(breaks = 0:12) 

pdf(file = "Plots/Descriptives/f_alpha.pdf", width = 6, height = 4)
print(g) 
dev.off() 
setEPS()
postscript("Plots/Descriptives/f_alpha.eps", width = 6, height = 4)
print(g)
dev.off()