Transform Sto vsp benchmark into an AbstractStochasticBenchmark{true}

Author: BatyLeo

src/DynamicVehicleScheduling/scenario.jl

@@ -51,7 +51,7 @@ function Utils.generate_scenario(
 end
 
 function Utils.generate_scenario(
-    ::DynamicVehicleSchedulingBenchmark, rng::AbstractRNG; instance, kwargs...
+    ::DynamicVehicleSchedulingBenchmark, rng::AbstractRNG; instance::Instance, kwargs...
 )
     return generate_scenario(instance; rng)
 end

src/StochasticVehicleScheduling/StochasticVehicleScheduling.jl

@@ -6,6 +6,7 @@ export plot_instance, plot_solution
 export compact_linearized_mip,
     compact_mip, column_generation_algorithm, local_search, deterministic_mip
 export evaluate_solution, is_feasible
+export VSPScenario, build_stochastic_instance
 
 using ..Utils
 using DocStringExtensions: TYPEDEF, TYPEDFIELDS, TYPEDSIGNATURES
@@ -30,7 +31,7 @@ using JuMP:
 using Plots: Plots, plot, plot!, scatter!, annotate!, text
 using Printf: @printf
 using Random: Random, AbstractRNG, MersenneTwister
-using SparseArrays: sparse
+using SparseArrays: sparse, SparseMatrixCSC, nnz
 using Statistics: quantile, mean
 
 include("utils.jl")
@@ -41,6 +42,8 @@ include("instance/city.jl")
 include("instance/features.jl")
 include("instance/instance.jl")
 
+include("scenario.jl")
+
 include("solution/solution.jl")
 include("solution/algorithms/mip.jl")
 include("solution/algorithms/column_generation.jl")
@@ -57,17 +60,53 @@ Data structure for a stochastic vehicle scheduling benchmark.
 # Fields
 $TYPEDFIELDS
 """
-@kwdef struct StochasticVehicleSchedulingBenchmark <: AbstractBenchmark
+@kwdef struct StochasticVehicleSchedulingBenchmark <: AbstractStochasticBenchmark{true}
     "number of tasks in each instance"
     nb_tasks::Int = 25
-    "number of scenarios in each instance"
+    "number of scenarios in each instance (only used to compute features and objective evaluation)"
     nb_scenarios::Int = 10
 end
 
+include("policies.jl")
+
 function Utils.objective_value(
     ::StochasticVehicleSchedulingBenchmark, sample::DataSample, y::BitVector
 )
-    return evaluate_solution(y, sample.instance)
+    stoch = build_stochastic_instance(sample.instance, sample.extra.scenarios)
+    return evaluate_solution(y, stoch)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Draw a single fresh [`VSPScenario`](@ref) for the given instance.
+Requires `store_city=true` (the default) when generating instances.
+"""
+function Utils.generate_scenario(
+    ::StochasticVehicleSchedulingBenchmark, rng::AbstractRNG; instance::Instance, kwargs...
+)
+    @assert !isnothing(instance.city) "`generate_scenario` requires `store_city=true`"
+    return draw_scenario(instance.city, instance.graph, rng)
+end
+
+"""
+$TYPEDSIGNATURES
+"""
+function Utils.generate_baseline_policies(bench::StochasticVehicleSchedulingBenchmark)
+    return svs_generate_baseline_policies(bench)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Return the anticipative solver: a callable `(scenario::VSPScenario; instance, kwargs...) -> y`
+that solves the 1-scenario stochastic VSP via column generation.
+"""
+function Utils.generate_anticipative_solver(::StochasticVehicleSchedulingBenchmark)
+    return (scenario::VSPScenario; instance::Instance, kwargs...) -> begin
+        stochastic_inst = build_stochastic_instance(instance, [scenario])
+        return column_generation_algorithm(stochastic_inst)
+    end
 end
 
 """
@@ -84,7 +123,9 @@ policy = sample -> DataSample(; sample.context..., x=sample.x,
 dataset = generate_dataset(benchmark, N; target_policy=policy)
 ```
 
-If `store_city=false`, coordinates and city information are not stored in the instance.
+If `store_city=false`, coordinates and city information are not stored in the instance,
+and `generate_scenario` will not work. This can be used to save memory if you only need to evaluate
+solutions on a fixed set of scenarios.
 """
 function Utils.generate_instance(
     benchmark::StochasticVehicleSchedulingBenchmark,

src/StochasticVehicleScheduling/instance/instance.jl

@@ -26,6 +26,18 @@ end
 """
 $TYPEDSIGNATURES
 
+Display a compact summary of an [`Instance`](@ref): number of tasks, scenarios, and edges.
+"""
+function Base.show(io::IO, instance::Instance)
+    return print(
+        io,
+        "VSP Instance($(get_nb_tasks(instance)) tasks, $(get_nb_scenarios(instance)) scenarios, $(ne(instance.graph)) edges)",
+    )
+end
+
+"""
+$TYPEDSIGNATURES
+
 Return the acyclic directed graph corresponding to `city`.
 Each vertex represents a task. Vertices are ordered by start time of corresponding task.
 There is an edge from task u to task v the (end time of u + tie distance between u and v <= start time of v).

src/StochasticVehicleScheduling/policies.jl

@@ -0,0 +1,54 @@
+"""
+$TYPEDSIGNATURES
+
+SAA baseline policy: builds a stochastic instance from all K scenarios and solves
+via column generation.
+Returns a single labeled [`DataSample`](@ref) with `extra=(; scenarios)`.
+"""
+function svs_saa_policy(sample, scenarios)
+    stochastic_inst = build_stochastic_instance(sample.instance, scenarios)
+    y = column_generation_algorithm(stochastic_inst)
+    return [DataSample(; sample.context..., x=sample.x, y, extra=(; scenarios))]
+end
+
+"""
+$TYPEDSIGNATURES
+
+Deterministic baseline policy: solves the deterministic MIP (ignores scenario delays).
+Returns a single labeled [`DataSample`](@ref) with `extra=(; scenarios)`.
+"""
+function svs_deterministic_policy(sample, scenarios; model_builder=highs_model)
+    y = deterministic_mip(sample.instance; model_builder)
+    return [DataSample(; sample.context..., x=sample.x, y, extra=(; scenarios))]
+end
+
+"""
+$TYPEDSIGNATURES
+
+Local search baseline policy: builds a stochastic instance from all K scenarios and
+solves via local search heuristic.
+Returns a single labeled [`DataSample`](@ref) with `extra=(; scenarios)`.
+"""
+function svs_local_search_policy(sample, scenarios)
+    stochastic_inst = build_stochastic_instance(sample.instance, scenarios)
+    y = local_search(stochastic_inst)
+    return [DataSample(; sample.context..., x=sample.x, y, extra=(; scenarios))]
+end
+
+"""
+$TYPEDSIGNATURES
+
+Return the named baseline policies for [`StochasticVehicleSchedulingBenchmark`](@ref).
+Each policy has signature `(sample, scenarios) -> Vector{DataSample}`.
+"""
+function svs_generate_baseline_policies(::StochasticVehicleSchedulingBenchmark)
+    return (;
+        deterministic=Policy(
+            "Deterministic MIP", "Ignores delays", svs_deterministic_policy
+        ),
+        saa=Policy("SAA (col gen)", "Stochastic MIP over K scenarios", svs_saa_policy),
+        local_search=Policy(
+            "Local search", "Heuristic with K scenarios", svs_local_search_policy
+        ),
+    )
+end

src/StochasticVehicleScheduling/scenario.jl

@@ -0,0 +1,130 @@
+"""
+$TYPEDEF
+
+Represents a single scenario for the stochastic vehicle scheduling problem.
+
+# Fields
+$TYPEDFIELDS
+"""
+struct VSPScenario
+    "delays per task (length = nb_tasks + 2): scenario_end_time - nominal_end_time"
+    delays::Vector{Float64}
+    "scalar slack per edge for this scenario"
+    slacks::SparseMatrixCSC{Float64,Int}
+end
+
+"""
+$TYPEDSIGNATURES
+
+Display a compact summary of a [`VSPScenario`](@ref): number of tasks and edges.
+"""
+function Base.show(io::IO, s::VSPScenario)
+    return print(io, "VSPScenario($(length(s.delays) - 2) tasks)")
+end
+
+"""
+$TYPEDSIGNATURES
+
+Draw a single fresh scenario from the city's random distributions,
+independently of the stored scenario draws in the `City` struct.
+"""
+function draw_scenario(city::City, graph::AbstractGraph, rng::AbstractRNG)
+    tasks = city.tasks
+    N = length(tasks)
+
+    # 1. Draw inter-area factors for 24 hours (single scenario)
+    inter_area = zeros(24)
+    previous = 0.0
+    for h in 1:24
+        previous = (previous + 0.1) * rand(rng, city.random_inter_area_factor)
+        inter_area[h] = previous
+    end
+
+    # 2. Draw district delays for each district and 24 hours (single scenario)
+    nb_per_edge = size(city.districts, 1)
+    district_delays = [zeros(24) for _ in 1:nb_per_edge, _ in 1:nb_per_edge]
+    for x in 1:nb_per_edge
+        for y in 1:nb_per_edge
+            prev = 0.0
+            for h in 1:24
+                prev = scenario_next_delay(prev, city.districts[x, y].random_delay, rng)
+                district_delays[x, y][h] = prev
+            end
+        end
+    end
+
+    # 3. Draw task start times (single scenario per task)
+    scenario_start_time = [t.start_time + rand(rng, t.random_delay) for t in tasks]
+
+    # 4. Compute task end times for job tasks (indices 2:(N-1))
+    scenario_end_time = [t.end_time for t in tasks]
+    for i in 2:(N - 1)
+        task = tasks[i]
+        origin_x, origin_y = get_district(task.start_point, city)
+        dest_x, dest_y = get_district(task.end_point, city)
+
+        ξ₁ = scenario_start_time[i]
+        ξ₂ = ξ₁ + district_delays[origin_x, origin_y][hour_of(ξ₁)]
+        ξ₃ = ξ₂ + (task.end_time - task.start_time) + inter_area[hour_of(ξ₂)]
+        scenario_end_time[i] = ξ₃ + district_delays[dest_x, dest_y][hour_of(ξ₃)]
+    end
+
+    # 5. Compute delays: scenario_end_time - nominal_end_time
+    delays = scenario_end_time .- [t.end_time for t in tasks]
+
+    # 6. Compute scalar slack for each edge in this scenario
+    I_idx = [src(e) for e in edges(graph)]
+    J_idx = [dst(e) for e in edges(graph)]
+    slack_vals = map(edges(graph)) do e
+        u = src(e)
+        v = dst(e)
+        origin_x, origin_y = get_district(tasks[u].end_point, city)
+        dest_x, dest_y = get_district(tasks[v].start_point, city)
+        ξ₁ = scenario_end_time[u]
+        ξ₂ = ξ₁ + district_delays[origin_x, origin_y][hour_of(ξ₁)]
+        ξ₃ =
+            ξ₂ +
+            distance(tasks[u].end_point, tasks[v].start_point) +
+            inter_area[hour_of(ξ₂)]
+        perturbed_arrival = ξ₃ + district_delays[dest_x, dest_y][hour_of(ξ₃)]
+        perturbed_travel_time = perturbed_arrival - ξ₁
+        return (v < N ? scenario_start_time[v] : Inf) -
+               (tasks[u].end_time + perturbed_travel_time)
+    end
+    slacks = sparse(I_idx, J_idx, slack_vals, N, N)
+
+    return VSPScenario(delays, slacks)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Build a stochastic [`Instance`](@ref) from a base instance and a vector of fresh
+[`VSPScenario`](@ref)s. Each scenario contributes one column to the `intrinsic_delays`
+matrix and one entry per edge to the `slacks` sparse matrix.
+"""
+function build_stochastic_instance(instance::Instance, scenarios::Vector{VSPScenario})
+    K = length(scenarios)
+    nb_nodes = length(first(scenarios).delays)
+    intrinsic_delays = Matrix{Float64}(undef, nb_nodes, K)
+    for (k, s) in enumerate(scenarios)
+        intrinsic_delays[:, k] = s.delays
+    end
+
+    graph = instance.graph
+    N = nv(graph)
+    I_idx = [src(e) for e in edges(graph)]
+    J_idx = [dst(e) for e in edges(graph)]
+    slack_vecs = [[scenarios[k].slacks[src(e), dst(e)] for k in 1:K] for e in edges(graph)]
+    new_slacks = sparse(I_idx, J_idx, slack_vecs, N, N)
+
+    return Instance(
+        graph,
+        instance.features,
+        new_slacks,
+        intrinsic_delays,
+        instance.vehicle_cost,
+        instance.delay_cost,
+        instance.city,
+    )
+end