Merge pull request #67 from JuliaDecisionFocusedLearning/parametric-anticipative-solvers

Implement anticipative solvers for both stochastic benchmarks

Author: BatyLeo

src/ContextualStochasticArgmax/ContextualStochasticArgmax.jl

@@ -118,6 +118,26 @@ end
 
 include("policies.jl")
 
+"""
+$TYPEDSIGNATURES
+
+Generates the anticipative solver for the benchmark.
+"""
+function Utils.generate_anticipative_solver(::ContextualStochasticArgmaxBenchmark)
+    return AnticipativeSolver()
+end
+
+"""
+$TYPEDSIGNATURES
+
+Generates the parametric anticipative solver for the benchmark.
+"""
+function Utils.generate_parametric_anticipative_solver(
+    ::ContextualStochasticArgmaxBenchmark
+)
+    return AnticipativeSolver()
+end
+
 export ContextualStochasticArgmaxBenchmark
 
 end

src/ContextualStochasticArgmax/policies.jl

@@ -1,5 +1,3 @@
-using Statistics: mean
-
 """
 $TYPEDSIGNATURES
 
@@ -28,3 +26,35 @@ Each policy has signature `(ctx_sample, scenarios) -> Vector{DataSample}`.
 function Utils.generate_baseline_policies(::ContextualStochasticArgmaxBenchmark)
     return (; saa=Policy("SAA", "argmax of mean scenarios", csa_saa_policy))
 end
+
+"""
+$TYPEDEF
+
+A policy that acts with perfect information about the future scenario.
+"""
+struct AnticipativeSolver end
+
+function Base.show(io::IO, ::AnticipativeSolver)
+    return print(io, "Anticipative solver for ContextualStochasticArgmaxBenchmark")
+end
+
+"""
+$TYPEDSIGNATURES
+
+Evaluate the anticipative policy for a given `scenario`. 
+Returns the optimal action `one_hot_argmax(scenario)`.
+"""
+function (::AnticipativeSolver)(scenario; context...)
+    return one_hot_argmax(scenario)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Evaluate the anticipative policy with a parametric prediction `θ` and a `scenario`.
+Returns the optimal action for the combined signal `one_hot_argmax(scenario + θ)`.
+"""
+function (::AnticipativeSolver)(θ, scenario; context...)
+    ξ = scenario + θ
+    return one_hot_argmax(ξ)
+end

src/StochasticVehicleScheduling/StochasticVehicleScheduling.jl

@@ -47,6 +47,7 @@ include("solution/algorithms/mip.jl")
 include("solution/algorithms/column_generation.jl")
 include("solution/algorithms/local_search.jl")
 include("solution/algorithms/deterministic_mip.jl")
+include("solution/algorithms/anticipative_solver.jl")
 
 include("maximizer.jl")
 
@@ -113,13 +114,29 @@ end
 $TYPEDSIGNATURES
 
 Return the anticipative solver: a callable `(scenario::VSPScenario; instance, kwargs...) -> y`
-that solves the 1-scenario stochastic VSP via column generation.
+that solves the 1-scenario stochastic VSP.
+
+# Keyword Arguments
+- `model_builder`: a function returning an empty `JuMP.Model` with a solver attached (defaults to `scip_model`).
 """
-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
+function Utils.generate_anticipative_solver(
+    ::StochasticVehicleSchedulingBenchmark; model_builder=scip_model
+)
+    return AnticipativeSolver(; model_builder=model_builder)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Return the parametric anticipative solver: a callable `(θ, scenario::VSPScenario; instance, kwargs...) -> y`.
+
+# Keyword Arguments
+- `model_builder`: a function returning an empty `JuMP.Model` with a solver attached (defaults to `scip_model`).
+"""
+function Utils.generate_parametric_anticipative_solver(
+    ::StochasticVehicleSchedulingBenchmark; model_builder=scip_model
+)
+    return AnticipativeSolver(; model_builder=model_builder)
 end
 
 """

src/StochasticVehicleScheduling/solution/algorithms/anticipative_solver.jl

@@ -0,0 +1,22 @@
+@kwdef struct AnticipativeSolver{M,A}
+    model_builder::M = scip_model
+    single_scenario_algorithm::A = compact_mip
+end
+
+function Base.show(io::IO, ::AnticipativeSolver)
+    return print(io, "Anticipative solver for StochasticVehicleSchedulingBenchmark")
+end
+
+function (solver::AnticipativeSolver)(scenario; instance::Instance, kwargs...)
+    stochastic_inst = build_stochastic_instance(instance, [scenario])
+    return solver.single_scenario_algorithm(
+        stochastic_inst; model_builder=solver.model_builder, kwargs...
+    )
+end
+
+function (solver::AnticipativeSolver)(θ, scenario; instance::Instance, kwargs...)
+    stochastic_inst = build_stochastic_instance(instance, [scenario])
+    return solver.single_scenario_algorithm(
+        stochastic_inst, θ; model_builder=solver.model_builder, kwargs...
+    )
+end

src/StochasticVehicleScheduling/solution/algorithms/mip.jl

@@ -6,7 +6,11 @@ Quadratic constraints are linearized using Mc Cormick linearization.
 Note: If you have Gurobi, use `grb_model` as `model_builder` instead of `highs_model`.
 """
 function compact_linearized_mip(
-    instance::Instance; scenario_range=nothing, model_builder=scip_model, silent=true
+    instance::Instance,
+    θ=nothing;
+    scenario_range=nothing,
+    model_builder=scip_model,
+    silent=true,
 )
     (; graph, slacks, intrinsic_delays, vehicle_cost, delay_cost) = instance
     nb_nodes = nv(graph)
@@ -28,13 +32,17 @@ function compact_linearized_mip(
     @variable(model, R[v in nodes, ω in Ω] >= 0) # propagated delay of job v
     @variable(model, yR[u in nodes, v in nodes, ω in Ω; has_edge(graph, u, v)] >= 0) # yR[u, v] = y[u, v] * R[u, ω]
 
-    @objective(
-        model,
-        Min,
+    obj = (
         delay_cost * sum(sum(R[v, ω] for v in job_indices) for ω in Ω) / nb_scenarios # average total delay
-            +
-            vehicle_cost * sum(y[1, v] for v in job_indices) # nb_vehicles
+        +
+        vehicle_cost * sum(y[1, v] for v in job_indices) # nb_vehicles
     )
+    if !isnothing(θ)
+        @assert length(θ) == ne(graph)
+        obj += sum(θ[a] * y[src(edge), dst(edge)] for (a, edge) in enumerate(edges(graph)))
+    end
+
+    @objective(model, Min, obj)
 
     # Flow contraints
     @constraint(
@@ -103,7 +111,11 @@ Note: If you have Gurobi, use `grb_model` as `model_builder` instead of `highs_m
     You need to use a solver that supports quadratic constraints to use this method.
 """
 function compact_mip(
-    instance::Instance; scenario_range=nothing, model_builder=scip_model, silent=true
+    instance::Instance,
+    θ=nothing;
+    scenario_range=nothing,
+    model_builder=scip_model,
+    silent=true,
 )
     (; graph, slacks, intrinsic_delays, vehicle_cost, delay_cost) = instance
     nb_nodes = nv(graph)
@@ -124,13 +136,17 @@ function compact_mip(
     @variable(model, R[v in nodes, ω in Ω] >= 0) # propagated delay of job v
     @variable(model, yR[u in nodes, v in nodes, ω in Ω; has_edge(graph, u, v)] >= 0) # yR[u, v] = y[u, v] * R[u, ω]
 
-    @objective(
-        model,
-        Min,
+    obj = (
         delay_cost * sum(sum(R[v, ω] for v in job_indices) for ω in Ω) / nb_scenarios # average total delay
-            +
-            vehicle_cost * sum(y[1, v] for v in job_indices) # nb_vehicles
+        +
+        vehicle_cost * sum(y[1, v] for v in job_indices) # nb_vehicles
     )
+    if !isnothing(θ)
+        @assert length(θ) == ne(graph)
+        obj += sum(θ[a] * y[src(edge), dst(edge)] for (a, edge) in enumerate(edges(graph)))
+    end
+
+    @objective(model, Min, obj)
 
     # Flow contraints
     @constraint(
@@ -170,6 +186,8 @@ $TYPEDSIGNATURES
 SAA variant: build stochastic instance from `scenarios` then solve via
 [`compact_mip`](@ref).
 """
-function compact_mip(instance::Instance, scenarios::Vector{VSPScenario}; kwargs...)
-    return compact_mip(build_stochastic_instance(instance, scenarios); kwargs...)
+function compact_mip(
+    instance::Instance, scenarios::Vector{VSPScenario}, θ=nothing; kwargs...
+)
+    return compact_mip(build_stochastic_instance(instance, scenarios), θ; kwargs...)
 end

test/contextual_stochastic_argmax.jl

@@ -38,6 +38,31 @@
     @test first(d1).x ≈ first(d2).x
 end
 
+@testset "Parametric Anticipative Solver - ContextualStochasticArgmax" begin
+    using DecisionFocusedLearningBenchmarks
+
+    b = ContextualStochasticArgmaxBenchmark(; n=5, d=3, seed=0)
+    dataset = generate_dataset(b, 2; contexts_per_instance=1, nb_scenarios=1)
+    sample = first(dataset)
+    scenario = generate_scenario(b, StableRNG(0); sample.context...)
+
+    solver = generate_anticipative_solver(b)
+    parametric_solver = generate_parametric_anticipative_solver(b)
+
+    # 1. Zero perturbation equivalence
+    θ_zero = zeros(eltype(scenario), size(scenario))
+    @test parametric_solver(θ_zero, scenario; sample.context...) ==
+        solver(scenario; sample.context...)
+
+    # 2. Extreme perturbation
+    θ_extreme = zeros(eltype(scenario), size(scenario))
+    θ_extreme[1] = 1000.0  # Force dimension 1
+    y_extreme = parametric_solver(θ_extreme, scenario; sample.context...)
+
+    @test y_extreme[1] == 1.0     # Only dimension 1 should be active
+    @test sum(y_extreme) ≈ 1.0    # One-hot preserved
+end
+
 @testset "csa_saa_policy" begin
     using DecisionFocusedLearningBenchmarks
 

test/vsp.jl

@@ -95,4 +95,28 @@
     sample = unlabeled[1]
     y_anticipative = anticipative_solver(sample.scenario; sample.context...)
     @test y_anticipative isa BitVector
+
+    # Extract necessary dependencies
+    parametric_solver = generate_parametric_anticipative_solver(b)
+    nb_edges = ne(sample.instance.graph)
+
+    # 1. Zero perturbation equivalence
+    θ_zero = zeros(nb_edges)
+    y_zero = parametric_solver(θ_zero, sample.scenario; sample.context...)
+
+    @test y_zero == y_anticipative
+    @test y_zero isa BitVector
+
+    # 2. Perturbation execution
+    θ_random = randn(nb_edges)
+    y_rand = parametric_solver(θ_random, sample.scenario; sample.context...)
+
+    @test length(y_rand) == nb_edges
+    @test y_rand isa BitVector
+
+    # 3. High negative perturbation on edge forces activation (it's minimization)
+    θ_extreme = zeros(nb_edges)
+    θ_extreme[1] = -100000.0  # large negative pull for edge 1
+    y_extreme = parametric_solver(θ_extreme, sample.scenario; sample.context...)
+    @test y_extreme[1] == 1.0 # BitVector
 end