# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#### available kernels ####
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#' @title triangle kernel
#'
#' @description Function implementing the triangle kernel.
#'
#' @param d The distance from the event
#' @param bw The bandwidth used for the kernel
#' @return The estimated density
#' @export
#' @examples
#' #This is an internal function, no example provided
triangle_kernel <- function(d, bw){
u <- d/bw
k <- 1 - abs(u)
k <- k/bw
k <- ifelse(abs(d)>bw,0,k)
return(k)
}
#' @title Uniform kernel
#'
#' @description Function implementing the uniform kernel.
#'
#' @param d The distance from the event
#' @param bw The bandwidth used for the kernel
#' @return The estimated density
#' @export
#' @examples
#' #This is an internal function, no example provided
uniform_kernel <- function(d, bw){
k <- 1/(2*bw)
k <- ifelse(abs(d)>bw,0,k)
return(k)
}
#' @title Epanechnikov kernel
#'
#' @description Function implementing the epanechnikov kernel.
#'
#' @param d The distance from the event
#' @param bw The bandwidth used for the kernel
#' @return The estimated density
#' @export
#' @examples
#' #This is an internal function, no example provided
epanechnikov_kernel <- function(d, bw){
u <- d/bw
k <- (3/4) * (1-u**2)
k <- k/bw
k <- ifelse(abs(d)>bw,0,k)
return(k)
}
#' @title Quartic kernel
#'
#' @description Function implementing the quartic kernel.
#'
#' @param d The distance from the event
#' @param bw The bandwidth used for the kernel
#' @return The estimated density
#' @export
#' @examples
#' #This is an internal function, no example provided
quartic_kernel <- function(d, bw){
u <- d/bw
k <- (15/16)*(1-u**2)**2
k <- k / bw
k <- ifelse(abs(d)>bw,0,k)
return(k)
}
#' @title Triweight kernel
#'
#' @description Function implementing the triweight kernel.
#'
#' @param d The distance from the event
#' @param bw The bandwidth used for the kernel
#' @return The estimated density
#' @export
#' @examples
#' #This is an internal function, no example provided
triweight_kernel <- function(d, bw){
u <- d/bw
k <- (35/32)*(1-u**2)**3
k <- k / bw
k <- ifelse(abs(d)>bw,0,k)
return(k)
}
#' @title Tricube kernel
#'
#' @description Function implementing the tricube kernel.
#'
#' @param d The distance from the event
#' @param bw The bandwidth used for the kernel
#' @return The estimated density
#' @export
#' @examples
#' #This is an internal function, no example provided
tricube_kernel <- function(d, bw){
u <- d/bw
k <- (70/81)*(1-abs(u)**3)**3
k <- k / bw
k <- ifelse(abs(d)>bw,0,k)
return(k)
}
#' @title Cosine kernel
#'
#' @description Function implementing the cosine kernel.
#'
#' @param d The distance from the event
#' @param bw The bandwidth used for the kernel
#' @return The estimated density
#' @export
#' @examples
#' #This is an internal function, no example provided
cosine_kernel <- function(d, bw){
u <- abs(d/bw)
k <- (pi/4) * cos((pi/2)*u)
k <- k / bw
k <- ifelse(abs(d)>bw,0,k)
return(k)
}
#' @title Gaussian kernel
#'
#' @description Function implementing the gaussian kernel.
#'
#' @param d The distance from the event
#' @param bw The bandwidth used for the kernel
#' @return The estimated density
#' @export
#' @examples
#' #This is an internal function, no example provided
gaussian_kernel <- function(d, bw){
u <- d/bw
t1 <- 1/(sqrt(2*pi))
t2 <- exp(-1 * (1/2) * u**2 )
k <- t1*t2
k <- k/bw
k <- ifelse(abs(d)>bw,0,k)
return(k)
}
#' @title Scaled gaussian kernel
#'
#' @description Function implementing the scaled gaussian kernel.
#'
#' @param d The distance from the event
#' @param bw The bandwidth used for the kernel
#' @return The estimated density
#' @export
#' @examples
#' #This is an internal function, no example provided
gaussian_kernel_scaled <- function(d, bw){
bw2 <- bw/3
t1 <- 1/(bw2 * sqrt(2*pi))
t2 <- exp(-1 * (d**2 / (2*bw2**2)))
k <- t1*t2
k <- ifelse(abs(d)>bw,0,k)
return(k)
}
#' @title Select kernel function
#'
#' @description select the kernel function with its name.
#'
#' @param name The name of the kernel to use
#' @return A kernel function
#' @keywords internal
#' @examples
#' #This is an internal function, no example provided
select_kernel <- function(name){
list_kernels <- c("triangle", "gaussian", "scaled gaussian", "tricube",
"cosine" ,"triweight", "quartic", 'epanechnikov',
'uniform')
if((name %in% list_kernels)==FALSE){
allkernels <- paste(list_kernels,collapse = ", ")
stop(paste("the name of the kernel function must be one of",
allkernels, sep=" "))
}
if(name=="gaussian"){
return(gaussian_kernel)
}
if(name == "scaled gaussian"){
return(gaussian_kernel_scaled)
}
if(name=="triangle"){
return(triangle_kernel)
}
if(name=="tricube"){
return(tricube_kernel)
}
if(name=="triweight"){
return(triweight_kernel)
}
if(name=="quartic"){
return(quartic_kernel)
}
if(name=="cosine"){
return(cosine_kernel)
}
if(name=="epanechnikov"){
return(epanechnikov_kernel)
}
if(name=="uniform"){
return(uniform_kernel)
}
}
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#### Functions for adaptative bandwidth calculation ####
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#' @title Geometric mean
#'
#' @description Function to calculate the geometric mean.
#'
#' @param x A vector of numeric values
#' @param na.rm A boolean indicating if we filter the NA values
#' @return The geometric mean of x
#' @keywords internal
#' @examples
#' #This is an internal function, no example provided
gm_mean <- function(x, na.rm=TRUE){
exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
}
#' @title Gamma parameter for Abramson’s adaptive bandwidth
#'
#' @description Function to calculate the gamma parameter in Abramson’s
#' smoothing regimen.
#'
#' @param k a vector of numeric values (the estimated kernel densities)
#' @return the gamma parameter in Abramson’s smoothing regimen
#' @keywords internal
#' @examples
#' #This is an internal function, no example provided
calc_gamma <- function(k){
return(gm_mean(k**(-1/2)))
}
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#### Functions to perform the simple NKDE ####
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#' @title Simple NKDE algorithm
#'
#' @description Function to perform the simple nkde.
#'
#' @param graph a graph object from igraph representing the network
#' @param events a feature collection of points representing the events. It must be
#' snapped on the network, and be nodes of the network. A column vertex_id
#' must indicate for each event its corresponding node
#' @param samples a a feature collection of points representing the sampling points.
#' The samples must be snapped on the network. A column edge_id must indicate
#' for each sample on which edge it is snapped.
#' @param bws a vector indicating the kernel bandwidth (in meters) for each
#' event
#' @param kernel_func a function obtained with the function select_kernel
#' @param nodes a a feature collection of points representing the nodes of the network
#' @param edges a a feature collection of linestrings representing the edges of the network
#' @param div The divisor to use for the kernel. Must be "n" (the number of events within the radius around each sampling point), "bw" (the bandwidth) "none" (the simple sum).
#' @return a dataframe with two columns. sum_k is the sum for each sample point
#' of the kernel values. n is the number of events influencing each sample
#' point
#' @keywords internal
#' @examples
#' #This is an internal function, no example provided
simple_nkde <- function(graph, events, samples, bws, kernel_func, nodes, edges, div = "bw"){
## step 1 set all values to 0
base_k <- rep(0,nrow(samples))
base_count <- rep(0,nrow(samples))
#if the number of event is 0, then return only 0 values
if(nrow(events)==0){
return(data.frame("sum_k"=base_k,
"n"=base_count))
}
sample_tree <- build_quadtree(samples)
edges_tree <- build_quadtree(edges)
## step2 iterate over each event
pb <- txtProgressBar(min = 0, max = nrow(events), style = 3)
for(i in 1:nrow(events)){
setTxtProgressBar(pb, i)
bw <- bws[[i]]
#extracting the starting values
e <- events[i,]
y <- e$vertex_id
w <- e$weight
#calculating the kernel values
samples_k <- ess_kernel(graph,y,bw, kernel_func, samples,
nodes, edges, sample_tree, edges_tree)
samples_count <- ifelse(samples_k>0,1,0)
if(div == "bw"){
base_k <- (samples_k * w / bw) + base_k
}else{
base_k <- samples_k * w + base_k
}
base_count <- base_count + (samples_count * w)
}
return(data.frame("sum_k"=base_k,
"n"=base_count))
}
#' @title Worker for simple NKDE algorithm
#'
#' @description The worker function to perform the simple nkde.
#'
#' @param graph a graph object from igraph representing the network
#' @param y the index of the actual event
#' @param bw a float indicating the kernel bandwidth (in meters)
#' @param kernel_func a function obtained with the function select_kernel
#' @param samples a a feature collection of points representing the sampling points. The
#' samples must be snapped on the network. A column edge_id must indicate for
#' each sample on which edge it is snapped.
#' @param nodes a a feature collection of points representing the nodes of the network
#' @param edges a a feature collection of linestrings representing the edges of the network
#' @param sample_tree a quadtree object, the spatial index of the samples
#' @param edges_tree a quadtree object, the spatial index of the edges
#' @importFrom igraph get.edge.attribute
#' @importFrom igraph ends distances
#' @keywords internal
#' @examples
#' #This is an internal function, no example provided
ess_kernel <- function(graph, y, bw, kernel_func, samples, nodes, edges, sample_tree, edges_tree){
samples_k <- rep(0,nrow(samples))
event_node <- nodes[y,]
buff <- st_buffer(event_node, dist = bw)
## step1 find all the samples in the radius
ok_samples <- spatial_request(buff,sample_tree,samples)
if(nrow(ok_samples)==0){
return(samples_k)
}
## step2 find all the edges in the radius
buff2 <- st_buffer(event_node, dist = bw+0.1*bw)
ok_edges <- spatial_request(buff2,edges_tree,edges)$edge_id
## Step3 for each edge, find its two vertices
vertices <- ends(graph,ok_edges,names = FALSE)
## step4 calculate the distance between the start node and each edge vertex
un_vertices <- unique(c(vertices[,1],vertices[,2]))
dist1 <- as.numeric(distances(graph,y,to=un_vertices,mode="out"))
dist_table <- data.frame("vertex"=un_vertices,
"distance" = dist1)
## step5 aggregate all the data
df_edges <- data.frame("edge_id" = ok_edges,
"node1" = vertices[,1],
"node2" = vertices[,2]
)
A1 <- data.table(df_edges)
A2 <- data.table(df_edges)
B <- data.table(dist_table)
df_edges$d1 <- A1[B, on = c("node1" = "vertex"),
names(B) := mget(paste0("i.", names(B)))]$distance
df_edges$d2 <- A2[B, on = c("node2" = "vertex"),
names(B) := mget(paste0("i.", names(B)))]$distance
#df_edges$d1 <- left_join(df_edges,dist_table,by=c("node1"="vertex"))$distance
#df_edges$d2 <- left_join(df_edges,dist_table,by=c("node2"="vertex"))$distance
## now, we will calculate for each sample, the minimum distance for the two distances
df1 <- df_edges[c("edge_id","node1","d1")]
df2 <- df_edges[c("edge_id","node2","d2")]
## getting the first part of the distances
df_samples1 <- data.table(st_drop_geometry(ok_samples))
df_samples2 <- data.table(st_drop_geometry(ok_samples))
B <- data.table(df1)
C <- data.table(df2)
df_samples1[B, on = "edge_id", names(B) := mget(paste0("i.", names(B)))]
df_samples2[C, on = "edge_id", names(C) := mget(paste0("i.", names(C)))]
#df_samples1 <- left_join(ok_samples@data,df1,by="edge_id")
#df_samples2 <- left_join(ok_samples@data,df2,by="edge_id")
## getting the second part of the distance
nodes_ng <- st_drop_geometry(nodes)
start_nodes1 <- nodes_ng[df_samples1$node1,]
start_nodes2 <- nodes_ng[df_samples2$node2,]
dist1_b <- sqrt((df_samples1$X_coords - start_nodes1$X_coords)**2 +
(df_samples1$Y_coords - start_nodes1$Y_coords)**2)
dist2_b <- sqrt((df_samples2$X_coords - start_nodes2$X_coords)**2 +
(df_samples2$Y_coords - start_nodes2$Y_coords)**2)
## getting the full distances
full_dist1 <- dist1_b + df_samples1$d1
full_dist2 <- dist2_b + df_samples2$d2
ok_distances <-ifelse(full_dist1 < full_dist2, full_dist1, full_dist2)
##and then, the final distance
k <- kernel_func(ok_distances,bw)
samples_k[df_samples1$oid] <- k
return(samples_k)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.