Merge pull request #26 from BatyLeo/better-routing

Update stochastic routing and log algorithm info

Author: BatyLeo

.gitignore

@@ -6,3 +6,4 @@
 scripts
 .vscode
 debug.jl
+*.jld2

CITATION.bib

@@ -2,7 +2,7 @@ @misc{ConstrainedShortestPaths.jl
 	author  = {Léo Baty and contributors},
 	title   = {ConstrainedShortestPaths.jl},
 	url     = {https://github.com/BatyLeo/ConstrainedShortestPaths.jl},
-	version = {v0.4.0},
+	version = {v0.5.0},
 	year    = {2024},
-	month   = {8}
+	month   = {9}
 }

Project.toml

@@ -1,7 +1,7 @@
 name = "ConstrainedShortestPaths"
 uuid = "b3798467-87dc-4d99-943d-35a1bd39e395"
 authors = ["Léo Baty and contributors"]
-version = "0.4.0"
+version = "0.5.0"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
@@ -15,7 +15,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 DataStructures = "0.18"
 DocStringExtensions = "0.9"
 Graphs = "1"
-PiecewiseLinearFunctions = "0.1"
+PiecewiseLinearFunctions = "0.3"
 SparseArrays = "<0.0.1,1"
 Statistics = "<0.0.1,1"
 julia = "1.10"
@@ -26,6 +26,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
@@ -36,4 +37,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 
 [targets]
-test = ["Aqua", "Documenter", "GLPK", "Graphs", "JET", "JuMP", "JuliaFormatter", "Plots", "Random", "SparseArrays", "StableRNGs", "Test", "UnicodePlots"]
+test = ["Aqua", "Documenter", "GLPK", "Graphs", "JET", "JLD2", "JuMP", "JuliaFormatter", "Plots", "Random", "SparseArrays", "StableRNGs", "Test", "UnicodePlots"]

src/algorithms.jl

@@ -48,6 +48,9 @@ function generalized_a_star(instance::CSPInstance, bounds; kwargs...)
     c_star = Inf
     p_star = [origin_vertex]
 
+    nb_cuts_with_bounds = 0
+    nb_cuts_with_dominance = 0
+
     while !isempty(L)
         p = dequeue!(L)
         v = p[end]
@@ -64,6 +67,7 @@ function generalized_a_star(instance::CSPInstance, bounds; kwargs...)
             end
             forward_resources[q] = rq
             c = instance.cost_function(rq, bounds[w])
+
             if c < c_star # cut using bounds
                 if w == destination_vertex # if destination is reached
                     c_star = c
@@ -72,11 +76,16 @@ function generalized_a_star(instance::CSPInstance, bounds; kwargs...)
                     remove_dominated!(M[w], rq)
                     push!(M[w], rq)
                     enqueue!(L, q => c)
+                else
+                    nb_cuts_with_dominance += 1
                 end
+            else
+                nb_cuts_with_bounds += 1
             end
         end
     end
-    return (; p_star, c_star)
+    info = (; nb_cuts_with_bounds, nb_cuts_with_dominance)
+    return (; p_star, c_star, info, bounds)
 end
 
 """
@@ -100,6 +109,7 @@ function generalized_a_star(instance::ForwardCSPInstance; kwargs...)
     push!(M[origin_vertex], forward_resources[empty_path])
     c_star = Inf
     p_star = [origin_vertex]
+    nb_cuts_with_dominance = 0
 
     while !isempty(L)
         p = dequeue!(L)
@@ -126,10 +136,13 @@ function generalized_a_star(instance::ForwardCSPInstance; kwargs...)
                 remove_dominated!(M[w], rq)
                 push!(M[w], rq)
                 enqueue!(L, q => c)
+            else
+                nb_cuts_with_dominance += 1
             end
         end
     end
-    return (; p_star, c_star)
+    info = (; nb_cuts_with_dominance)
+    return (; p_star, c_star, info)
 end
 
 """

src/examples/stochastic_routing.jl

@@ -54,27 +54,27 @@ function (f::StochasticForwardFunction)(q::StochasticForwardResource)
     return StochasticForwardResource(new_c, new_xi, new_λ), true
 end
 
-function _backward_scenario(g::PiecewiseLinearFunction, delay::Float64, slack::Float64)
+function _backward_scenario(g::PiecewiseLinearFunction, delay::Float64, slack::Float64, λᵥ)
     f = if slack == Inf
         piecewise_linear()
     else
         piecewise_linear(1.0, slack, delay)
     end
-    return f + g ∘ f
+    return f + g ∘ f - λᵥ
 end
 
 function (f::StochasticBackwardFunction)(q::StochasticBackwardResource)
     return StochasticBackwardResource(
         [
-            _backward_scenario(g, delay, slack) for
+            _backward_scenario(g, delay, slack, f.λ_value) for
             (g, delay, slack) in zip(q.g, f.delays, f.slacks)
         ],
         f.λ_value + q.λ,
     )
 end
 
 function stochastic_cost(fr::StochasticForwardResource, br::StochasticBackwardResource)
-    λ_sum = fr.λ + br.λ
+    λ_sum = fr.λ
     cp = fr.c + mean(gj(Rj) for (gj, Rj) in zip(br.g, fr.xi))
     return cp - λ_sum
 end

test/examples/stochastic_routing.jl

@@ -1,6 +1,7 @@
 using ConstrainedShortestPaths
 using Graphs
 using SparseArrays
+using JLD2
 
 @testset "Path digraph" begin
     m = 5  # number of scenarios
@@ -132,9 +133,9 @@ end
 @testset "Random graphs" begin
     n = 5
     for nb_scenarios in 1:5
-        for nb_vertices in 10:15
-            for i in 1:n
-                rng = StableRNG(i)
+        for nb_vertices in 10:1:15
+            for seed in 1:n
+                rng = StableRNG(seed)
                 graph = random_acyclic_digraph(
                     nb_vertices, rng; all_connected_to_source_and_destination=true
                 )
@@ -163,6 +164,9 @@ end
                     graph, slack_matrix, delays, initial_paths; bounding=false
                 )
                 @test obj2 ≈ _obj2
+                # if !(obj2 ≈ _obj2)
+                #     jldsave("debug.jld2"; graph, slack_matrix, delays, initial_paths)
+                # end
 
                 # Exact resolution
                 obj, sol = solve_scenarios(graph, slack_matrix, delays)

test/utils.jl

@@ -101,12 +101,19 @@ function stochastic_PLNE(g, slacks, delays, initial_paths; bounding=true)
     new_paths = Vector{Int}[]
     cons = []
 
+    # c = 0
     while true
+        # c += 1
         optimize!(model)
         λ_val = value.(λ)
         (; c_star, p_star) = stochastic_routing_shortest_path(
             g, slacks, delays, λ_val; bounding
         )
+        # if c == 2
+        #     jldsave("debug2.jld2"; graph, slacks, delays, λ_val, c_star, p_star)
+        # end
+        # @info λ_val
+        # @info (; c_star, p_star)
         full_cost =
             c_star + vehicle_cost + sum(λ_val[v] for v in job_indices if v in p_star)
         @assert path_cost(p_star, slacks, delays) + vehicle_cost ≈ full_cost