library(tidygraph)
library(igraph)
library(stplanr)
library(cancensus)
library(osmdata)
library(dplyr)
library(sf)
library(sfnetworks)
library(ggplot2)
library(stringr)

study_area <- get_census(dataset='CA16', regions=list(CSD="5915"),
                     level='CSD', use_cache = FALSE,geo_format = "sf") %>%
  filter(name=="Vancouver (CY)") %>%
  st_transform(crs = 26910) 

ggplot(study_area) +
  geom_sf() + 
  coord_sf() +
  theme_void()

bbx <- getbb("Vancouver, BC")

streets <- bbx %>%
  opq()%>%
  add_osm_feature(key = "highway", 
                  value=c("residential", "living_street",
                            "service","unclassified",
                            "pedestrian", "footway",
                            "track","path","motorway", "trunk",
                          "primary","secondary", 
                          "tertiary","motorway_link",
                          "trunk_link","primary_link",
                          "secondary_link",
                          "tertiary_link")) %>%
  osmdata_sf()
  
  
streets_sf <- st_transform(streets$osm_lines,crs=26910) %>% 
  st_intersection(.,st_geometry(study_area)) |> 
  st_cast(to = "MULTILINESTRING") |> 
  st_cast(to = "LINESTRING")

ggplot() +
  geom_sf(data = streets_sf,
          aes(color=highway),
          size = .4,
          alpha = .65)+
  theme_void()

# core base + tidy read
url  <- "https://geodash.vpd.ca/opendata/crimedata_download/crimedata_csv_all_years.zip"
zipf <- tempfile(fileext = ".zip")
csvd <- tempfile(fileext = ".csv")

download.file(url, destfile = zipf, mode = "wb")
csv_name <- unzip(zipf, list = TRUE)$Name[grepl("\\.csv$", unzip(zipf, list=TRUE)$Name)]
unzip(zipf, files = csv_name, exdir = dirname(csvd))
file.rename(file.path(dirname(csvd), csv_name), csvd)

# fast, friendly read
library(readr)
crime <- readr::read_csv(csvd, show_col_types = FALSE)
dplyr::glimpse(crime)

collisions <- crime %>% 
  filter(YEAR==2025 & str_detect(TYPE,"Vehicle Collision"))%>% 
  filter(X!=0) %>%
  st_as_sf(.,coords = c("X","Y"),crs=26910)

ggplot(collisions) +
  geom_sf() + 
  coord_sf() + 
  theme_void()

split_lines <- function(input_lines, max_length, id) {

  input_lines <- input_lines %>% ungroup()

  geom_column <- attr(input_lines, "sf_column")

  input_crs <- sf::st_crs(input_lines)

  input_lines[["geom_len"]] <- sf::st_length(input_lines[[geom_column]])

  attr(input_lines[["geom_len"]], "units") <- NULL
  input_lines[["geom_len"]] <- as.numeric(input_lines[["geom_len"]])

  too_short <- filter(select(all_of(input_lines),all_of(id), all_of(geom_column), geom_len), geom_len < max_length) %>% select(-geom_len)

  too_long <- filter(select(all_of(input_lines),all_of(id), all_of(geom_column), geom_len), geom_len >= max_length)

  rm(input_lines) # just to control memory usage in case this is big.

  too_long <- mutate(too_long,
                     pieces = ceiling(geom_len / max_length),
                     fID = 1:nrow(too_long)) %>%
    select(-geom_len)

  split_points <- sf::st_set_geometry(too_long, NULL)[rep(seq_len(nrow(too_long)), too_long[["pieces"]]),] %>%
    select(-pieces)

  split_points <- mutate(split_points, split_fID = row.names(split_points)) %>%
    group_by(fID) %>%
    mutate(piece = 1:n()) %>%
    mutate(start = (piece - 1) / n(),
           end = piece / n()) %>%
    ungroup()

  new_line <- function(i, f, t) {
    lwgeom::st_linesubstring(x = too_long[[geom_column]][i], from = f, to = t)[[1]]
  }

  split_lines <- apply(split_points[c("fID", "start", "end")], 1,
                       function(x) new_line(i = x[["fID"]], f = x[["start"]], t = x[["end"]]))

  rm(too_long)

  split_lines <- st_sf(split_points[c(id)], geometry = st_sfc(split_lines, crs = input_crs))

  lixel <- rbind(split_lines,too_short) %>% mutate(LIXID = row_number())

  return(lixel)
}

lixelize_network <- function(sf_network,max_lixel_length,uid){
  print("Splitting input spatial lines by lixel length...")
  target_lixel <- split_lines(input_lines = sf_network,max_length = max_lixel_length,id = uid)
  print("Create corresponding shortest distance network...")
  shortest_distance_network <- split_lines(input_lines = target_lixel,max_length = max_lixel_length/2,id = uid)
  return(list(target_lixel=target_lixel,shortest_distance_network=shortest_distance_network))
}

split_lines <- function(input_lines, max_length, id) {

  input_lines <- input_lines %>% ungroup()

  geom_column <- attr(input_lines, "sf_column")

  input_crs <- sf::st_crs(input_lines)

  input_lines[["geom_len"]] <- sf::st_length(input_lines[[geom_column]])

  attr(input_lines[["geom_len"]], "units") <- NULL
  input_lines[["geom_len"]] <- as.numeric(input_lines[["geom_len"]])

  too_short <- filter(select(input_lines,all_of(id), all_of(geom_column), geom_len), geom_len < max_length) %>% select(-geom_len)

  too_long <- filter(select(input_lines,all_of(id), all_of(geom_column), geom_len), geom_len >= max_length)

  rm(input_lines) # just to control memory usage in case this is big.

  too_long <- mutate(too_long,
                     pieces = ceiling(geom_len / max_length),
                     fID = 1:nrow(too_long)) %>%
    select(-geom_len)

  split_points <- sf::st_set_geometry(too_long, NULL)[rep(seq_len(nrow(too_long)), too_long[["pieces"]]),] %>%
    select(-pieces)

  split_points <- mutate(split_points, split_fID = row.names(split_points)) %>%
    group_by(fID) %>%
    mutate(piece = 1:n()) %>%
    mutate(start = (piece - 1) / n(),
           end = piece / n()) %>%
    ungroup()

  new_line <- function(i, f, t) {
    lwgeom::st_linesubstring(x = too_long[[geom_column]][i], from = f, to = t)[[1]]
  }

  split_lines <- apply(split_points[c("fID", "start", "end")], 1,
                       function(x) new_line(i = x[["fID"]], f = x[["start"]], t = x[["end"]]))

  rm(too_long)

  split_lines <- st_sf(split_points[c(id)], geometry = st_sfc(split_lines, crs = input_crs))

  lixel <- rbind(split_lines,too_short) %>% mutate(LIXID = row_number())

  return(lixel)
}

lixelize_network <- function(sf_network,max_lixel_length,uid){
  print("Splitting input spatial lines by lixel length...")
  target_lixel <- split_lines(input_lines = sf_network,max_length = max_lixel_length,id = uid)
  print("Create corresponding shortest distance network...")
  shortest_distance_network <- split_lines(input_lines = target_lixel,max_length = max_lixel_length/2,id = uid)
  return(list(target_lixel=target_lixel,shortest_distance_network=shortest_distance_network))
}

lixel_list_25m <- lixelize_network(
  sf_network = streets_sf,
  max_lixel_length = 25,
  uid = "osm_id"
    
)

str(lixel_list_25m)

st_snap_points = function(x, y, max_dist = 1000) {
  
  if (inherits(x, "sf")) n = nrow(x)
  if (inherits(x, "sfc")) n = length(x)
  
  out = do.call(c,
                lapply(seq(n), function(i) {
                  nrst = st_nearest_points(st_geometry(x)[i], y)
                  nrst_len = st_length(nrst)
                  nrst_mn = which.min(nrst_len)
                  if (as.vector(nrst_len[nrst_mn]) > max_dist) return(st_geometry(x)[i])
                  return(st_cast(nrst[nrst_mn], "POINT")[2])
                })
  )
  return(out)
}

coll_snp <- st_snap_points(collisions,streets_sf,max_dist = 20) #this only returns the geometry - doesn't preserve attributes

coll_snp <- st_sf(collisions %>%
                  st_drop_geometry() %>%
                  mutate(geom=coll_snp)) #rerturn the attributes to the snapped locations

#function for calculating the centre of 
st_line_midpoints <- function(sf_lines = NULL) {
  
  g <- st_geometry(sf_lines)
  
  g_mids <- lapply(g, function(x) {
    
    coords <- as.matrix(x)
    
    # this is just a copypaste of maptools:::getMidpoints()):
    get_mids <- function (coords) {
      dist <- sqrt((diff(coords[, 1])^2 + (diff(coords[, 2]))^2))
      dist_mid <- sum(dist)/2
      dist_cum <- c(0, cumsum(dist))
      end_index <- which(dist_cum > dist_mid)[1]
      start_index <- end_index - 1
      start <- coords[start_index, ]
      end <- coords[end_index, ]
      dist_remaining <- dist_mid - dist_cum[start_index]
      mid <- start + (end - start) * (dist_remaining/dist[start_index])
      return(mid)
    }
    
    mids <- st_point(get_mids(coords))
  })
  
  out <- st_sfc(g_mids, crs = st_crs(sf_lines))
  out <- st_sf(out)
  out <- bind_cols(out,sf_lines %>% st_drop_geometry()) %>% rename(geom=out)
}


network_kde <- function(lixel_list,point_process,bandwidth = 100,n_cores=1,attribute=1,point_process_is_lixel_midpoint=FALSE){
  

  #3. Create a network of lixels by establishing the network topology between lixels as well as between lixels and lxnodes.

  #Define topology for calculating shortest distances by converting lixel with half the length to calculate distances
  #The shorter the lixel length the more accurate the calculation of network distance from nearest node of source lixel to nearest node of target lixel
  print("Defining network topology...")
  
  require(parallel)
  require(doParallel)
  require(igraph)
  require(tidygraph)
  require(sf)  
  require(dplyr)


  # Create nodes at the start and end point of each edge

  nodes <- lixel_list$shortest_distance_network %>%
    st_coordinates() %>%
    as_tibble() %>%
    rename(LIXID = L1) %>%
    group_by(LIXID) %>%
    slice(c(1, n())) %>%
    ungroup() %>%
    mutate(start_end = rep(c('start', 'end'), times = n()/2))

  # Give each node a unique index

  nodes <- nodes %>%
    mutate(xy = paste(.$X, .$Y),
           xy = factor(xy, levels = unique(xy))) %>%
    group_by(xy)%>%
    mutate(nodeID = cur_group_id()) %>%
    ungroup() %>%
    select(-xy)

  # Combine the node indices with the edges

  start_nodes <- nodes %>%
    filter(start_end == 'start') %>%
    pull(nodeID)

  end_nodes <- nodes %>%
    filter(start_end == 'end') %>%
    pull(nodeID)

  lixel_list$shortest_distance_network = lixel_list$shortest_distance_network %>%
    mutate(from = start_nodes, to = end_nodes)

  # Remove duplicate nodes
  nodes <- nodes %>%
    distinct(nodeID, .keep_all = TRUE) %>%
    select(-c(LIXID, start_end)) %>%
    st_as_sf(coords = c('X', 'Y')) %>%
    st_set_crs(st_crs(lixel_list$shortest_distance_network))

  # Convert to tbl_graph
  graph <- tbl_graph(nodes = nodes, edges = as_tibble(lixel_list$shortest_distance_network), directed = FALSE)

  graph <- graph %>%
    activate(edges) %>%
    mutate(length = st_length(geometry))

  lixel_list$shortest_distance_network <- NULL

  #4. Create the center points of all the lixels for the target lixel (lxcenters)
  print("Calculating lixel midpoints (lxcenters)...")

  lxcenters <- st_line_midpoints(lixel_list$target_lixel)

  #5. Select a point process occurring within the road network

  #input as function parameter

  #6. For each point find its nearest lixel. Count the number of points nearest to a lixel and assigned as property of lixel.

  #Points should be snapped to network within some distance threshold prior to

  if(point_process_is_lixel_midpoint==FALSE){

  print("Counting number of events on each lixel...")

  point_process <- st_join(point_process,lixel_list$target_lixel["LIXID"],
                           join = st_is_within_distance, dist = 0.001) #for each point assign the LIXID that it falls on

  source_lixels <- point_process %>% #summarize the number of points by LIXID (e.g. count the points for each LIXID)
    filter(!is.na(LIXID)) %>%
    group_by(LIXID) %>%
    summarise(n_events = n(),`.groups`="drop") %>%
    st_drop_geometry()

  source_lixels <- inner_join(lxcenters,source_lixels,by="LIXID") #define geometry for source lixels
  }

  if(point_process_is_lixel_midpoint==TRUE){

    source_lixels <- lxcenters %>% mutate(n_events = 1)
  }

  print(paste0(sum(source_lixels$n_events)," events on ",nrow(source_lixels)," source lixels"))

  #7. Define a search bandwidth, measured with the short path network distance

  #input as function parameter

  #8. Calculate the shortest-path network distance from each source lixel to lxcenter of all its neighouring lixels within the seach bandwidth


  nearest_node_to_source_lixel <- st_nearest_feature(source_lixels,nodes)#find nodes from shortest distance network associated with each source lixel
  nearest_node_to_lxcentre <- st_nearest_feature(lxcenters,nodes)#find nodes from shortest distance network associated with each lxcenters

  print("Calculating distances from each source lixel to all other lixel centers... ")
  

  cl <- makeCluster(n_cores)
  registerDoParallel(cores=cl) #parallel computing

  distances <- foreach::foreach(i = 1:length(nearest_node_to_source_lixel),.packages = c("magrittr","igraph","tidygraph","sf")) %dopar% {
    temp <- distances(
      graph = graph,
      weights = graph %>% activate(edges) %>% pull(length),
      v = nearest_node_to_source_lixel[i]
    )[,nearest_node_to_lxcentre]

    data.frame(LIXID = lxcenters[temp<=max(bandwidth),]$LIXID,
               dist = temp[temp<=max(bandwidth)])
  }

 stopCluster(cl)

  rm("graph")


  # gauss <-function(x) {
  #
  #   t1 <- 1/(sqrt(2*pi))
  #   t2 <- exp(-((x^2)/(2*bandwidth^2)))
  #   n <- (1/bandwidth)*(t1*t2)
  #   n <- ifelse(x>bandwidth,0,n)
  #
  #
  #   return(n)
  #
  # }
  #
  # plot(gauss(0:(bandwidth*2)))

  quartic <-function(x,r) {
    K <- 3/pi
    t1 <- 1-(x^2/r^2)
    q <- ifelse(x>r,0,K*t1)
    q <- (1/r)*q

    return(q)
  }

  LIXID <- unlist(lapply(lapply(distances,`[`,"LIXID"),function(x) pull(x)))
  distances_list <- lapply(lapply(distances,`[`,"dist"),function(x) pull(x))

  d_list <- list()

  for(i in 1:length(bandwidth)){

    print(paste0("Applying kernel density estimator with bandwidth of ",bandwidth[i],"m ... "))

    d <- lapply(distances_list,function(x) quartic(x,bandwidth[i]))
    d <- mapply(`*`,d,attribute) #multiply by attribute of source lixel
    d <- mapply(`*`,d,source_lixels$n_events) #sum over number of events that occured on the source lixel

    d_list[[i]] <- data.frame(unlist(d))
  }

  d_cols <- as.data.frame(do.call(cbind,d_list))

  names(d_cols) <- paste0("kde_bw_",bandwidth)

  #sum densities over all lixels
  density <- cbind(LIXID,d_cols) %>%
    group_by(LIXID) %>%
    summarise(across(everything(),sum),`.groups`="drop")

  network_kde <- left_join(lixel_list$target_lixel,density,by = "LIXID") %>%
    mutate(length = st_length(.)) %>%
    replace(., is.na(.), 0)

  return(list(network_kde=network_kde,neighbours = distances))

}

n_cores <- bigstatsr::nb_cores()

system.time({
  collisions_network_kde <- network_kde(lixel_list = lixel_list_25m,
                                 point_process = coll_snp,
                                 bandwidth = c(50,100,150,250),
                                 n_cores = n_cores,
                                 point_process_is_lixel_midpoint = FALSE)
})

glimpse(collisions_network_kde$network_kde)

ggplot() +
  geom_sf(data=study_area,color=NA,fill="gray5")+
  geom_sf(data = collisions_network_kde$network_kde,aes(color=kde_bw_150))+
  scale_color_viridis_c(option = "C",direction = 1,name="Density") + 
  # geom_sf(data=mask,color=NA,fill="white")+
  theme_void()+
  theme(panel.background= element_rect(fill = NA,color=NA),
        legend.position = "top",
        legend.title.align = 0,
        legend.key.size = unit(1, 'cm'), #change legend key size
        legend.key.height = unit(0.5, 'cm'), #change legend key height
        legend.key.width = unit(1.5, 'cm'))+
  coord_sf(xlim = c(483690.9,499329.2),
           ylim = c(5450534.8,5463381.0))

ggplot() +
  geom_sf(data=study_area,color=NA,fill="gray5")+
  geom_sf(data = collisions_network_kde$network_kde,aes(color=kde_bw_150))+
  scale_color_viridis_c(option = "C",direction = 1,name="Density") + 
  theme_void()+
  theme(panel.background= element_rect(fill = NA,color=NA),
        legend.position = "top",
        legend.title.align = 0,
        legend.key.size = unit(1, 'cm'), #change legend key size
        legend.key.height = unit(0.5, 'cm'), #change legend key height
        legend.key.width = unit(1.5, 'cm'))+
  coord_sf(xlim = c(491321.39853389-1000,491321.39853389+1100),
           ylim = c(5459015.689874-1000,5459015.689874+1000))