From f84e2373bce2c436efeedf79065053043ba8aa9b Mon Sep 17 00:00:00 2001 From: Magnus Hiie Date: Thu, 15 Sep 2022 15:09:35 +0300 Subject: [PATCH] Implement ShortestPathManyToMany --- bidirectional_ch_n_to_n.go | 177 ++++++++++++++++++++++++++++++++ bidirectional_ch_n_to_n_test.go | 97 +++++++++++++++++ dijkstra_bidirectional.go | 8 ++ 3 files changed, 282 insertions(+) create mode 100644 bidirectional_ch_n_to_n.go create mode 100644 bidirectional_ch_n_to_n_test.go diff --git a/bidirectional_ch_n_to_n.go b/bidirectional_ch_n_to_n.go new file mode 100644 index 0000000..93a8d60 --- /dev/null +++ b/bidirectional_ch_n_to_n.go @@ -0,0 +1,177 @@ +package ch + +import ( + "container/heap" +) + +// ShortestPath Computes and returns shortest path and it's cost (extended Dijkstra's algorithm) +// +// If there are some errors then function returns '-1.0' as cost and nil as shortest path +// +// source User's definied ID of source vertex +// target User's definied ID of target vertex +func (graph *Graph) ShortestPathManyToMany(sources, targets []int64) ([][]float64, [][][]int64) { + endpoints := [directionsCount][]int64{sources, targets} + for d, directionEndpoints := range endpoints { + for i, endpoint := range directionEndpoints { + var ok bool + if endpoints[d][i], ok = graph.mapping[endpoint]; !ok { + endpoints[d][i] = -1 + } + } + } + return graph.shortestPathManyToMany(endpoints) +} + +func (graph *Graph) initShortestPathManyToMany(endpoints [directionsCount][]int64) ( + queryDist [directionsCount][][]float64, + processed [directionsCount][][]bool, + queues [directionsCount][]*vertexDistHeap, +) { + for d := forward; d < directionsCount; d++ { + queryDist[d] = make([][]float64, len(endpoints[d])) + processed[d] = make([][]bool, len(endpoints[d])) + queues[d] = make([]*vertexDistHeap, len(endpoints[d])) + for endpointIdx := range endpoints[d] { + queryDist[d][endpointIdx] = make([]float64, len(graph.Vertices)) + + for i := range queryDist[d][endpointIdx] { + queryDist[d][endpointIdx][i] = Infinity + } + + processed[d][endpointIdx] = make([]bool, len(graph.Vertices)) + + queues[d][endpointIdx] = &vertexDistHeap{} + + heap.Init(queues[d][endpointIdx]) + } + } + + return +} + +func (graph *Graph) shortestPathManyToMany(endpoints [directionsCount][]int64) ([][]float64, [][][]int64) { + queryDist, processed, queues := graph.initShortestPathManyToMany(endpoints) + + for d := forward; d < directionsCount; d++ { + + for endpointIdx, endpoint := range endpoints[d] { + processed[d][endpointIdx][endpoint] = true + + queryDist[d][endpointIdx][endpoint] = 0 + + heapEndpoint := &vertexDist{ + id: endpoint, + dist: 0, + } + heap.Push(queues[d][endpointIdx], heapEndpoint) + } + } + + return graph.shortestPathManyToManyCore(queryDist, processed, queues) +} + +func (graph *Graph) shortestPathManyToManyCore( + queryDist [directionsCount][][]float64, + processed [directionsCount][][]bool, + queues [directionsCount][]*vertexDistHeap, +) ([][]float64, [][][]int64) { + var prev [directionsCount][]map[int64]int64 + for d := forward; d < directionsCount; d++ { + prev[d] = make([]map[int64]int64, len(queues[d])) + for endpointIdx := range queues[d] { + prev[d][endpointIdx] = make(map[int64]int64) + } + } + estimates := make([][]float64, len(queues[forward])) + middleIDs := make([][]int64, len(queues[forward])) + for sourceEndpointIdx := range queues[forward] { + sourceEstimates := make([]float64, len(queues[backward])) + sourceMiddleIDs := make([]int64, len(queues[backward])) + estimates[sourceEndpointIdx] = sourceEstimates + middleIDs[sourceEndpointIdx] = sourceMiddleIDs + for targetEndpointIdx := range queues[backward] { + sourceEstimates[targetEndpointIdx] = Infinity + sourceMiddleIDs[targetEndpointIdx] = int64(-1) + } + } + + for { + queuesProcessed := false + for d := forward; d < directionsCount; d++ { + reverseDirection := (d + 1) % directionsCount + for endpointIdx := range queues[d] { + if queues[d][endpointIdx].Len() == 0 { + continue + } + queuesProcessed = true + graph.directionalSearchManyToMany(d, endpointIdx, queues[d][endpointIdx], processed[d][endpointIdx], processed[reverseDirection], queryDist[d][endpointIdx], queryDist[reverseDirection], prev[d][endpointIdx], estimates, middleIDs) + } + } + if !queuesProcessed { + break + } + } + paths := make([][][]int64, len(estimates)) + for sourceEndpointIdx, targetEstimates := range estimates { + targetPaths := make([][]int64, len(targetEstimates)) + paths[sourceEndpointIdx] = targetPaths + for targetEndpointIdx, estimate := range targetEstimates { + if estimate == Infinity { + targetEstimates[targetEndpointIdx] = -1 + continue + } + targetPaths[targetEndpointIdx] = graph.ComputePath(middleIDs[sourceEndpointIdx][targetEndpointIdx], prev[forward][sourceEndpointIdx], prev[backward][targetEndpointIdx]) + } + } + return estimates, paths +} + +func (graph *Graph) directionalSearchManyToMany( + d direction, endpointIndex int, q *vertexDistHeap, + localProcessed []bool, reverseProcessed [][]bool, + localQueryDist []float64, reverseQueryDist [][]float64, + prev map[int64]int64, estimates [][]float64, middleIDs [][]int64) { + + vertex := heap.Pop(q).(*vertexDist) + // if vertex.dist <= *estimate { // TODO: move to another place + localProcessed[vertex.id] = true + // Edge relaxation in a forward propagation + var vertexList []*incidentEdge + if d == forward { + vertexList = graph.Vertices[vertex.id].outIncidentEdges + } else { + vertexList = graph.Vertices[vertex.id].inIncidentEdges + } + for i := range vertexList { + temp := vertexList[i].vertexID + cost := vertexList[i].weight + if graph.Vertices[vertex.id].orderPos < graph.Vertices[temp].orderPos { + alt := localQueryDist[vertex.id] + cost + if localQueryDist[temp] > alt { + localQueryDist[temp] = alt + prev[temp] = vertex.id + node := &vertexDist{ + id: temp, + dist: alt, + } + heap.Push(q, node) + } + } + } + // } + for revEndpointIdx, revEndpointProcessed := range reverseProcessed { + if revEndpointProcessed[vertex.id] { + var sourceEndpoint, targetEndpoint int + if d == forward { + sourceEndpoint, targetEndpoint = endpointIndex, revEndpointIdx + } else { + targetEndpoint, sourceEndpoint = endpointIndex, revEndpointIdx + } + if vertex.dist+reverseQueryDist[revEndpointIdx][vertex.id] < estimates[sourceEndpoint][targetEndpoint] { + middleIDs[sourceEndpoint][targetEndpoint] = vertex.id + estimates[sourceEndpoint][targetEndpoint] = vertex.dist + reverseQueryDist[revEndpointIdx][vertex.id] + } + } + } +} diff --git a/bidirectional_ch_n_to_n_test.go b/bidirectional_ch_n_to_n_test.go new file mode 100644 index 0000000..1526180 --- /dev/null +++ b/bidirectional_ch_n_to_n_test.go @@ -0,0 +1,97 @@ +package ch + +import ( + "fmt" + "math" + "testing" +) + +func TestManyToManyShortestPath(t *testing.T) { + g := Graph{} + err := graphFromCSV(&g, "./data/pgrouting_osm.csv") + if err != nil { + t.Error(err) + return + } + t.Log("Please wait until contraction hierarchy is prepared") + g.PrepareContractionHierarchies() + t.Log("TestShortestPath is starting...") + u := []int64{106600, 106600, 69618} + v := []int64{5924, 81611, 69618, 68427, 68490} + correctAns := [][]float64{ + {61089.42195558673, 94961.78959757874, 78692.8292369651, 61212.00481622628, 71101.1080090782}, + {61089.42195558673, 94961.78959757874, 78692.8292369651, 61212.00481622628, 71101.1080090782}, + {19135.6581215226, -2, -2, -2, -2}, + } + correctPath := [][]int{ + {418, 866, 591, 314, 353}, + {418, 866, 591, 314, 353}, + {160, -2, -2, -2, -2}, + } + ans, path := g.ShortestPathManyToMany(u, v) + // t.Log("ShortestPathManyToMany returned", ans, path) + for sourceIdx := range u { + for targetIdx := range v { + if correctPath[sourceIdx][targetIdx] != -2 && len(path[sourceIdx][targetIdx]) != correctPath[sourceIdx][targetIdx] { + t.Errorf("Num of vertices in path should be %d, but got %d", correctPath[sourceIdx][targetIdx], len(path[sourceIdx][targetIdx])) + return + } + if correctAns[sourceIdx][targetIdx] != -2 && Round(ans[sourceIdx][targetIdx], 0.00005) != Round(correctAns[sourceIdx][targetIdx], 0.00005) { + t.Errorf("Length of path should be %f, but got %f", correctAns[sourceIdx][targetIdx], ans[sourceIdx][targetIdx]) + return + } + } + } + + t.Log("TestShortestPath is Ok!") +} + +func BenchmarkShortestPathManyToMany(b *testing.B) { + g := Graph{} + err := graphFromCSV(&g, "./data/pgrouting_osm.csv") + if err != nil { + b.Error(err) + } + b.Log("Please wait until contraction hierarchy is prepared") + g.PrepareContractionHierarchies() + b.Log("BenchmarkShortestPathManyToMany is starting...") + b.ResetTimer() + + for k := 0.; k <= 12; k++ { + n := int(math.Pow(2, k)) + b.Run(fmt.Sprintf("%s/%d/vertices-%d", "CH shortest path (many to many)", n, len(g.Vertices)), func(b *testing.B) { + for i := 0; i < b.N; i++ { + u := []int64{106600} + v := []int64{5924, 81611, 69618, 68427, 68490} + ans, path := g.ShortestPathManyToMany(u, v) + _, _ = ans, path + } + }) + } +} + +func BenchmarkOldWayShortestPathManyToMany(b *testing.B) { + g := Graph{} + err := graphFromCSV(&g, "data/pgrouting_osm.csv") + if err != nil { + b.Error(err) + } + b.Log("Please wait until contraction hierarchy is prepared") + g.PrepareContractionHierarchies() + b.Log("BenchmarkOldWayShortestPathManyToMany is starting...") + b.ResetTimer() + + for k := 0.; k <= 12; k++ { + n := int(math.Pow(2, k)) + b.Run(fmt.Sprintf("%s/%d/vertices-%d", "CH shortest path (many to many)", n, len(g.Vertices)), func(b *testing.B) { + for i := 0; i < b.N; i++ { + u := int64(106600) + v := []int64{5924, 81611, 69618, 68427, 68490} + for vv := range v { + ans, path := g.ShortestPath(u, v[vv]) + _, _ = ans, path + } + } + }) + } +} diff --git a/dijkstra_bidirectional.go b/dijkstra_bidirectional.go index 2709798..776c31f 100644 --- a/dijkstra_bidirectional.go +++ b/dijkstra_bidirectional.go @@ -5,6 +5,14 @@ import ( "container/list" ) +type direction int + +const ( + forward direction = iota + backward + directionsCount +) + // ShortestPath Computes and returns shortest path and it's cost (extended Dijkstra's algorithm) // // If there are some errors then function returns '-1.0' as cost and nil as shortest path