Merge pull request #28 from BatyLeo/convex-meet-option

Add an option for using convex meet in bounding

Author: BatyLeo

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.5.1},
+	version = {v0.6.0},
 	year    = {2024},
-	month   = {10}
+	month   = {12}
 }

Project.toml

@@ -1,7 +1,7 @@
 name = "ConstrainedShortestPaths"
 uuid = "b3798467-87dc-4d99-943d-35a1bd39e395"
 authors = ["Léo Baty and contributors"]
-version = "0.5.1"
+version = "0.6.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.4"
+PiecewiseLinearFunctions = "0.4.2"
 SparseArrays = "<0.0.1,1"
 Statistics = "<0.0.1,1"
 julia = "1.10"

src/ConstrainedShortestPaths.jl

@@ -13,7 +13,8 @@ using Graphs:
     edges,
     outneighbors,
     induced_subgraph
-using PiecewiseLinearFunctions: PiecewiseLinearFunction
+using PiecewiseLinearFunctions:
+    PiecewiseLinearFunction, convex_meet, remove_redundant_breakpoints
 using SparseArrays: sparse
 using Statistics: mean
 

src/algorithms.jl

@@ -199,6 +199,56 @@ end
 """
 $TYPEDSIGNATURES
 
+Compute all paths below threshold.
+"""
+function generalized_a_star_with_threshold(
+    instance::ForwardCSPInstance, threshold::Float64; kwargs...
+)
+    (; graph, origin_vertex, destination_vertex, is_useful) = instance
+
+    empty_path = [origin_vertex]
+
+    forward_resources = Dict(empty_path => instance.origin_forward_resource)
+    L = PriorityQueue{Vector{Int},Float64}(
+        empty_path => instance.cost_function(forward_resources[empty_path])
+    )
+
+    c_star = Float64[]
+    p_star = Vector{Int}[]
+
+    while !isempty(L)
+        p = dequeue!(L)
+        v = p[end]
+        for w in outneighbors(graph, v)
+            if !is_useful[w]
+                continue
+            end
+            q = copy(p)
+            push!(q, w)
+            rp = forward_resources[p]
+            rq, is_feasible = instance.forward_functions[v, w](rp; kwargs...)
+            if !is_feasible
+                continue
+            end
+            forward_resources[q] = rq
+            c = instance.cost_function(rq)
+            if w == destination_vertex # if destination is reached
+                if c < threshold
+                    push!(p_star, copy(q))
+                    push!(c_star, c)
+                else # else add path to queue
+                    enqueue!(L, q => c)
+                end
+            end
+            # else, discard path (i.e. do nothing)
+        end
+    end
+    return (; p_star, c_star)
+end
+
+"""
+$TYPEDSIGNATURES
+
 Compute the shortest path of `instance`.
 """
 function generalized_constrained_shortest_path(instance::CSPInstance; kwargs...)
@@ -216,8 +266,14 @@ $TYPEDSIGNATURES
 Compute shortest path between first and last nodes of `instance`
 """
 function generalized_constrained_shortest_path_with_threshold(
-    instance::CSPInstance{T,G}, threshold::Float64; kwargs...
-) where {T,G<:AbstractGraph}
+    instance::CSPInstance, threshold::Float64; kwargs...
+)
     bounds = compute_bounds(instance; kwargs...)
     return generalized_a_star_with_threshold(instance, bounds, threshold; kwargs...)
 end
+
+function generalized_constrained_shortest_path_with_threshold(
+    instance::ForwardCSPInstance, threshold::Float64; kwargs...
+)
+    return generalized_a_star_with_threshold(instance, threshold; kwargs...)
+end

src/examples/stochastic_routing.jl

@@ -21,7 +21,7 @@ $TYPEDEF
 # Fields
 $TYPEDFIELDS
 """
-struct StochasticBackwardResource
+struct StochasticBackwardResource{is_convex}
     "piecewise linear function for each scenario, take as input propagated delay and outputs total delay of partail path"
     g::Vector{PiecewiseLinearFunction{Float64}}
 end
@@ -70,9 +70,19 @@ function Base.:<=(r1::StochasticForwardResource, r2::StochasticForwardResource)
     return true
 end
 
-function Base.min(r1::StochasticBackwardResource, r2::StochasticBackwardResource)
+function Base.min(
+    r1::StochasticBackwardResource{true}, r2::StochasticBackwardResource{true}
+)
+    new_g = convex_meet.(r1.g, r2.g)
+    # new_g = remove_redundant_breakpoints.(convex_meet.(r1.g, r2.g); atol=1e-8)
+    return StochasticBackwardResource{true}(new_g)
+end
+
+function Base.min(
+    r1::StochasticBackwardResource{false}, r2::StochasticBackwardResource{false}
+)
     new_g = min.(r1.g, r2.g)
-    return StochasticBackwardResource(new_g)
+    return StochasticBackwardResource{false}(new_g)
 end
 
 function (f::StochasticForwardFunction)(q::StochasticForwardResource)
@@ -90,8 +100,10 @@ function _backward_scenario(g::PiecewiseLinearFunction, delay::Float64, slack::F
     return f + g ∘ f - λᵥ
 end
 
-function (f::StochasticBackwardFunction)(q::StochasticBackwardResource)
-    return StochasticBackwardResource([
+function (f::StochasticBackwardFunction)(
+    q::StochasticBackwardResource{is_convex}
+) where {is_convex}
+    return StochasticBackwardResource{is_convex}([
         _backward_scenario(g, delay, slack, f.λ_value) for
         (g, delay, slack) in zip(q.g, f.intrinsic_delays, f.slacks)
     ])
@@ -116,6 +128,7 @@ function create_instance(
     origin_vertex::T=one(T),
     destination_vertex::T=nv(graph),
     bounding=true,
+    use_convex_resources=true,
 ) where {T}
     @assert λ_values[origin_vertex] == 0.0 && λ_values[destination_vertex] == 0.0
     @assert all(intrinsic_delays[origin_vertex] .== 0.0)
@@ -139,7 +152,7 @@ function create_instance(
             StochasticBackwardFunction(slacks[u, v], intrinsic_delays[v, :], λ_values[v])
             for (u, v) in zip(I, J)
         ]
-        destination_backward_resource = StochasticBackwardResource([
+        destination_backward_resource = StochasticBackwardResource{use_convex_resources}([
             piecewise_linear() for _ in 1:nb_scenarios
         ])
 
@@ -189,6 +202,7 @@ function stochastic_routing_shortest_path(
     origin_vertex::T=one(T),
     destination_vertex::T=nv(graph),
     bounding=true,
+    use_convex_resources=true,
 ) where {T}
     instance = create_instance(
         graph,
@@ -198,6 +212,7 @@ function stochastic_routing_shortest_path(
         origin_vertex=origin_vertex,
         destination_vertex=destination_vertex,
         bounding=bounding,
+        use_convex_resources=use_convex_resources,
     )
     return generalized_constrained_shortest_path(instance)
 end
@@ -215,6 +230,7 @@ function stochastic_routing_shortest_path_with_threshold(
     origin_vertex::T=one(T),
     destination_vertex::T=nv(graph),
     bounding=true,
+    use_convex_resources=true,
     threshold,
 ) where {T}
     instance = create_instance(
@@ -225,6 +241,7 @@ function stochastic_routing_shortest_path_with_threshold(
         bounding=bounding,
         origin_vertex=origin_vertex,
         destination_vertex=destination_vertex,
+        use_convex_resources=use_convex_resources,
     )
     return generalized_constrained_shortest_path_with_threshold(instance, threshold)
 end

test/examples/stochastic_routing.jl

@@ -124,6 +124,12 @@ end
         @test c_star == 5
         @test p_star == [1, 3, 4, 5]
         @test path_cost(p_star, slack_matrix, delays) == c_star
+        (; c_star, p_star) = stochastic_routing_shortest_path(
+            graph, slack_matrix, delays; use_convex_resources=true
+        )
+        @test c_star == 5
+        @test p_star == [1, 3, 4, 5]
+        @test path_cost(p_star, slack_matrix, delays) == c_star
     end
 
     # @testset "Detour with slack" begin
@@ -140,7 +146,7 @@ end
 @testset "Random graphs" begin
     n = 5
     for nb_scenarios in 1:5
-        for nb_vertices in 10:1:15
+        for nb_vertices in 10:1:20
             for seed in 1:n
                 rng = StableRNG(seed)
                 graph = random_acyclic_digraph(
@@ -169,12 +175,9 @@ end
                     graph, slack_matrix, delays, initial_paths
                 )
                 _value, _obj2, _paths, _dual, _dual_new = stochastic_PLNE(
-                    graph, slack_matrix, delays, initial_paths; bounding=false
+                    graph, slack_matrix, delays, initial_paths; bounding=true
                 )
                 @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)
@@ -198,6 +201,10 @@ end
                     additional_paths, costs = stochastic_routing_shortest_path_with_threshold(
                         graph, slack_matrix, delays, λ_val; threshold
                     )
+                    _additional_paths, _costs = stochastic_routing_shortest_path_with_threshold(
+                        graph, slack_matrix, delays, λ_val; threshold, bounding=true
+                    )
+                    @test costs ≈ _costs
                     full_paths = unique(vcat(initial_paths, paths, additional_paths))
                     obj4, y4 = column_generation(graph, slack_matrix, delays, full_paths)
                     @test obj4 ≈ cupp || obj4 < cupp

test/utils.jl

@@ -80,7 +80,9 @@ function path_cost(path, slacks, delays, λ=zeros(size(delays, 1)))
     return C + vehicle_cost
 end
 
-function stochastic_PLNE(g, slacks, delays, initial_paths; bounding=true)
+function stochastic_PLNE(
+    g, slacks, delays, initial_paths; bounding=true, use_convex_resources=true
+)
     nb_nodes = nv(g)
     job_indices = 2:(nb_nodes - 1)
 
@@ -107,7 +109,7 @@ function stochastic_PLNE(g, slacks, delays, initial_paths; bounding=true)
         optimize!(model)
         λ_val = value.(λ)
         (; c_star, p_star) = stochastic_routing_shortest_path(
-            g, slacks, delays, λ_val; bounding
+            g, slacks, delays, λ_val; bounding, use_convex_resources
         )
         # if c == 2
         #     jldsave("debug2.jld2"; graph, slacks, delays, λ_val, c_star, p_star)