Add files via upload

Author: vlarmet

README.Rmd

@@ -69,25 +69,30 @@ Data has to be a 3 columns data.frame or matrix containing from, to and a cost/d
 ###Path algorithms 
 Path algorithms proposed by the package are :  
 
-  - uni-directional Dijkstra algorithm,
-  - bi-directional Dijkstra algorithm,
-  - uni-directional A* algorithm  
-  - New bi-directional A* algorithm (Piljs & Post, 2009 : see http://repub.eur.nl/pub/16100/ei2009-10.pdf)
-  - *one-to-one* bi-directional Dijkstra adapted to contraction hierarchies (Geisberger & al., 2008) 
-  - *many-to-many* bi-directional Dijkstra adapted to contraction hierarchies (Geisberger & al., 2008) 
+  - **1** uni-directional Dijkstra algorithm,
+  - **2** bi-directional Dijkstra algorithm,
+  - **3** uni-directional A* algorithm  
+  - **4** New bi-directional A* algorithm (Piljs & Post, 2009 : see http://repub.eur.nl/pub/16100/ei2009-10.pdf)
+  - **5** *one-to-one* bi-directional Dijkstra adapted to contraction hierarchies (Geisberger & al., 2008) 
+  - **6** *many-to-many* bi-directional Dijkstra adapted to contraction hierarchies (Geisberger & al., 2008) 
+  - **7** PHAST algorithm (Hardware-accelerated shortest path trees), *one-to-all* algorithm adapted to contraction hierarchies (Delling & al., 2011)
 
 
-The choice between all the algorithms is available for *one-to-one* calculation like  `get_distance_pair` and `get_path_pair` on a non-contracted graph.  
+*1*, *2*, *3* and *4* are available for **one-to-one** calculation in `get_distance_pair` and `get_path_pair` functions on a **non-contracted** graph.
 In these functions, uni-directional Dijkstra algorithm is stopped when the destination node is reached.  
 `A*` and `NBA*` are relevant if geographic coordinates of all nodes are provided. Note that coordinates should be expressed in a **projection system**.  
 To be accurate and efficient, `A*` and `NBA*` algorithms should use an admissible heuristic function (here the Euclidean distance), i.e cost and heuristic function must be expressed in the same unit.  
 In `cppRouting`, heuristic function `h` for a node (n) is defined such that :   
 **h(n,d) = ED(n,d) / k**
 with *h* the heuristic, *ED* the Euclidean distance, *d* the destination node and a constant *k*.  
-
 So in the case where coordinates are expressed in meters and cost is expressed in time, *k* is the maximum speed allowed on the road. By default, constant is 1 and is designed for graphs with cost expressed in the same unit than coordinates (for example in meters).  
-  
-If coordinates cannot be provided, bi-directional Dijkstra algorithm is the best option in terms of performance.
+If coordinates cannot be provided, bi-directional Dijkstra algorithm is the best option in terms of performance.  
+
+*5* is used for **one-to-one** calculation in `get_distance_pair` and `get_path_pair` functions on a **contracted** graph.  
+
+*1* is used for **one-to-many** calculation in `get_distance_matrix` function on a **non-contracted** graph.  
+
+*6* and *7* are available for **one-to-many** calculation in `get_distance_matrix` function on a **contracted** graph.  
 
 ##Examples
 ###Prepare data
@@ -314,7 +319,7 @@ Initially created for *one-to-one* queries, it has been extended to *many-to-man
 This technique is composed of two phases:  
 
 - preprocessing phase called *contraction* with `cpp_contract` function
-- query phase : a slightly modified version of bidirectional search for `one-to-one` query, available in `get_distance_pair` and `get_path_pair`; and a `many-to-many` algorithm using buckets available in `get_distance_matrix` function.  
+- query phase : a slightly modified version of bidirectional search for `one-to-one` query, available in `get_distance_pair` and `get_path_pair`; PHAST algorithm and a `many-to-many` algorithm using buckets available in `get_distance_matrix` function.  
 
 Contraction phase consists of iteratively removing a vertex **v** from the graph and creating a shortcut for each pair **(u,w)** of **v**'s neighborhood if the shortest path from **u** to **w** contains **v**. To be efficient and avoid creating too much shortcuts, vertices have to be ordered according to several heuristics. The two heuristics used by `cppRouting` are :  
 
@@ -362,7 +367,7 @@ kable(bench, align = "c",row.names = F) %>%
                  "nba : new bidirectional A* on the original graph"), notation="none")
 
 ```
-Here are the plots (in log-log) of query time improvement factor :  
+Here are the plots (in log-log) of query time improvement factor of *one to one CH* algorithm compared to bidirectional Dijkstra and NBA :  
 ```{r,echo=FALSE,message=FALSE,warning=FALSE}
 library(ggthemes)
 bench<-read.csv2("benchmark_contract_pair.csv")
@@ -391,40 +396,44 @@ p1
 As we can see on the plot, the larger is the graph, the higher is the benefit of using contraction hierarchies. For OSM Europe, query time can be faster by a factor of 1000 compared to bidirectional Dijkstra and 600 to NBA.  
 
 #####Distance matrix
-Here are the measurements of contraction time and query time (in second) of contraction hierarchies on different graphs.  
-Matrix are square (i.e the sets of source and target nodes are of equal length)
+Here are the measurements of query time (in second) of contraction hierarchies on different graphs.  
+We compare *PHAST* and *many to many CH* to Dijkstra algorithm on square matrix (i.e the sets of source and target nodes are of equal length).  
+
 ```{r,echo=FALSE,message=FALSE,warning=FALSE}
 bench<-read.csv2("benchmark_contract_matrix.csv",stringsAsFactors = F)
-bench<-bench[bench$algorithm %in% c("ch","normal"),]
+bench<-bench[bench$algorithm %in% c("mch","phast","normal"),]
 colnames(bench)<-gsub("\\."," ",colnames(bench))
 colnames(bench)<-gsub("X","",colnames(bench))
 bench$algorithm[bench$algorithm=="normal"]<-"Dijkstra"
 
 kable(bench, align = "c",row.names = F) %>%
   kable_styling(full_width = F) %>%
   column_spec(1, bold = T) %>%
-  row_spec(c(1,3,5), bold = TRUE) %>%
+  row_spec(c(1,4,7), bold = TRUE) %>%
   collapse_rows(columns = 1:2, valign = "middle") %>%
   add_header_above(c(" " = 2, "|S|=|T|" = 4))%>%
-  add_footnote(c("ch : many-to-many bidirectional search on the contracted graph",
+  add_footnote(c("mch : many-to-many bidirectional search on the contracted graph",
+                 "phast : phast algorithm on the contracted graph",
                  "Dijkstra : Dijkstra search on the original graph",
                  "|S| : number of source nodes",
                  "|T| : number of target nodes"), notation="none")
 
 ```
-Here are the plots (in log-log) of query time improvement factor for matrix :  
+Here are the plots (in log-log) of query time improvement factor of *PHAST* and *many to many CH* compared to Dijkstra algorithm :  
+
 ```{r,echo=FALSE,message=FALSE,warning=FALSE}
 library(ggthemes)
 bench<-read.csv2("benchmark_contract_matrix.csv")
-bench<-bench[bench$algorithm %in% c("ratio"),c(1,3:ncol(bench))]
+bench<-bench[bench$algorithm %in% c("ratio mch","ratio phast"),]
 
-test2<-melt(bench,id.vars = c("Graph"))
+test2<-melt(bench,id.vars = c("Graph","algorithm"))
 test2$variable<-as.character(test2$variable)
 test2$variable<-as.numeric(gsub("X","",test2$variable))
+test2$algorithm<-gsub("ratio ","",test2$algorithm)
 
 options(scipen=999)
 p2<-ggplot(data=test2,
-          aes(x=variable,y=value,colour=Graph))+
+          aes(x=variable,y=value,colour=Graph,linetype=algorithm))+
   geom_line(size=1)+
   geom_rangeframe(colour="black")+
   theme_tufte()+
@@ -438,7 +447,45 @@ p2=p2+scale_y_log10(labels=function(x) format(x, big.mark = ",", scientific = FA
 p2
 ```
 
-Benefits are less important than *one-to-one* queries but still interesting. For OSM Europe, query time can be faster by a factor of 90. 
+Benefits are less important than *one-to-one* queries but still interesting. For OSM Europe, query time can be faster by a factor of 90.  
+PHAST's improvement is constant since it iteratively perform an *one-to-all* search, just like original Dijkstra.  
+*many to many CH* is well adapted for **square matrix**.  
+
+Here are the plots of query time of *PHAST* and *many to many CH* on assymetric matrix (i.e. number of source and number of target are unequal) with *|S| / |T|* the number of sources divided by the number of targets :  
+
+```{r,echo=FALSE,message=FALSE,warning=FALSE}
+test<-read.csv2("assymetry_8000.csv")
+test$lab<-paste0(2^test$p2/2^test$p1)
+test<-test[,c("lab","seconds","dijkstra","phast","graph")]
+test<-melt(test,id.vars=c("lab","graph"))
+
+labb<-test$lab[1:13]
+test<-test[test$lab %in% labb[4:13],]
+labb<-labb[4:13]
+test$id<-rep(1:10,9)
+test$variable<-as.character(test$variable)
+test$variable<-ifelse(test$variable=="seconds","Many-to-many CH",test$variable)
+test$graph<-as.character(test$graph)
+test$graph<-ifelse(test$graph=="europe","OSM Europe",test$graph)
+test$graph<-ifelse(test$graph=="France OSM","OSM France",test$graph)
+test$graph<-ifelse(test$graph=="France","ROUTE 500",test$graph)
+
+options(scipen=999)
+p<-ggplot(data=test[test$variable!="dijkstra",],aes(x=id,y=value,colour=graph,linetype=variable))+
+  geom_line(size=1)+
+  scale_x_continuous(labels=labb,breaks=1:10)+
+  scale_y_log10()+
+  theme_minimal()+
+  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
+  labs(colour="Graph",linetype="Algorithm")+
+  ylab("seconds")+
+  xlab("|S| / |T|")
+  
+p 
+```
+
+Note that this plot is the same for *|T| / |S|*.  
+*PHAST* algorithm is much faster for rectangular matrix. The rate *|S| / |T|* where *many to many CH* is better varies according the graph size. For example, if we have to calculate a distance matrix between 10000 sources and 10 targets (or 10 sources and 10000 targets) on OSM France, we must use *PHAST*. On the other hand, if we want a matrix of 10000 sources and 8000 targets, we use *many to many CH* algorithm.  
 
 ###Compute isochrones
 Let's compute isochrones around Dijon city  
@@ -475,6 +522,7 @@ p
 
 ```
 
+
 #Applications
 
 **Except application 4, all indicators are calculated at the country scale but for the limited `R`'s ability to plot large shapefile, only one region is mapped.**  
@@ -538,11 +586,11 @@ The shortest travel time is computed with the `cppRouting` function `get_distanc
 In order to compute multiple distances from one source, original uni-directional Dijkstra algorithm is ran without early stopping.  
 We compute travel time from all commune nodes to all maternity ward nodes (i.e ~36000*400 distances).
 ```{r,echo=TRUE,message=FALSE,warning=FALSE,include=TRUE}
-#Distance matrix (few seconds to compute)
-dists<-get_distance_matrix(graph,
+#Distance matrix on contracted graph
+dists<-get_distance_matrix(graph3,
                            from=ndcom$id_noeud,
                            to=ndcom$id_noeud[ndcom$com %in% maternity$CODGEO],
-                           allcores=TRUE)
+                           algorithm = "phast")#because of the rectangular shape of the matrix
 #We extract each minimum travel time for all the communes
 dists2<-data.frame(node=ndcom$id_noeud,mindist=apply(dists,1,min,na.rm=T))
 #Joining commune IDs to nodes
@@ -675,7 +723,8 @@ roads2$to_id<-as.character(roads2$to_id)
 microbenchmark(igraph=test_igraph<-distances(graph_igraph,origin,to=destination,weights = E(graph_igraph)$weight,mode="out"),
                dodgr=test_dodgr<-dodgr_dists(graph=data.frame(roads2),from=origin,to=destination,parallel=FALSE),
                cpprouting=test_cpp<-get_distance_matrix(graph,origin,destination,allcores = FALSE),
-               contracted=test_cpp_cont<-get_distance_matrix(graph3,origin,destination,allcores = FALSE),
+               contr_mch=test_mch<-get_distance_matrix(graph3,origin,destination,algorithm = "mch",allcores = FALSE),
+               contr_phast=test_phast<-get_distance_matrix(graph3,origin,destination,algorithm = "phast",allcores = FALSE),
                times=1)
 
 ```
@@ -684,14 +733,15 @@ Even if we add the preprocessing time (i.e. 20s) to the query time, the whole pr
 
 **Ouput**  
 ```{r,echo=TRUE}
-head(cbind(test_igraph[,1],test_dodgr[,1],test_cpp[,1],test_cpp_cont[,1]))
+head(cbind(test_igraph[,1],test_dodgr[,1],test_cpp[,1],test_mch[,1],test_phast[,1]))
 ```
 
 **Distance matrix : parallel**  
 ```{r,echo=TRUE,warning=FALSE}
 microbenchmark(dodgr=test_dodgr<-dodgr_dists(graph=data.frame(roads2),from=origin,to=destination,parallel=TRUE),
                cpprouting=test_cpp<-get_distance_matrix(graph,origin,destination,allcores = TRUE),
-               contracted=test_cpp_cont<-get_distance_matrix(graph3,origin,destination,allcores = TRUE),
+               contr_mch=test_mch<-get_distance_matrix(graph3,origin,destination,algorithm = "mch",allcores = TRUE),
+               contr_phast=test_phast<-get_distance_matrix(graph3,origin,destination,algorithm = "phast",allcores = TRUE),
                times=1)
 ```