diff --git a/Project.toml b/Project.toml index b19c4aa..09e60e4 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.2.4" [deps] Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" +Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" ConstrainedShortestPaths = "b3798467-87dc-4d99-943d-35a1bd39e395" DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" @@ -13,7 +14,10 @@ Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0" +InferOpt = "4846b161-c94e-4150-8dac-c7ae193c601f" Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" +IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e" +JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -31,6 +35,7 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" [compat] Colors = "0.13.1" +Combinatorics = "1.0.3" ConstrainedShortestPaths = "0.6.0" DataDeps = "0.7" Distributions = "0.25" @@ -39,7 +44,10 @@ Flux = "0.14, 0.15, 0.16" Graphs = "1.11" HiGHS = "1.9" Images = "0.26.1" +InferOpt = "0.7.0" Ipopt = "1.6" +IterTools = "1.10.0" +JSON = "0.21.4" JuMP = "1.22" LaTeXStrings = "1.4.0" LinearAlgebra = "1" diff --git a/docs/make.jl b/docs/make.jl index 21f5480..1946e7f 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -12,8 +12,8 @@ api_dir = joinpath(@__DIR__, "src", "api") api_files = map(x -> joinpath("api", x), readdir(api_dir)) tutorial_files = readdir(tutorial_dir) md_tutorial_files = [split(file, ".")[1] * ".md" for file in tutorial_files] -benchmark_files = readdir(benchmarks_dir) -md_benchmark_files = [split(file, ".")[1] * ".md" for file in benchmark_files] +benchmark_files = [joinpath("benchmarks", e) for e in readdir(benchmarks_dir)] +# md_benchmark_files = [split(file, ".")[1] * ".md" for file in benchmark_files] include_tutorial = true @@ -25,20 +25,14 @@ if include_tutorial end makedocs(; - modules=[DecisionFocusedLearningBenchmarks, DecisionFocusedLearningBenchmarks.Warcraft], + modules=[DecisionFocusedLearningBenchmarks], authors="Members of JuliaDecisionFocusedLearning", sitename="DecisionFocusedLearningBenchmarks.jl", - format=Documenter.HTML(), + format=Documenter.HTML(; size_threshold=typemax(Int)), pages=[ "Home" => "index.md", "Tutorials" => include_tutorial ? md_tutorial_files : [], - "Benchmark problems list" => [ - "benchmarks/subset_selection.md", - "benchmarks/fixed_size_shortest_path.md", - "benchmarks/warcraft.md", - "benchmarks/portfolio_optimization.md", - "benchmarks/vsp.md", - ], + "Benchmark problems list" => benchmark_files, "API reference" => api_files, ], ) diff --git a/docs/src/api/0_interface.md b/docs/src/api/0_interface.md index 1b0a22c..6363833 100644 --- a/docs/src/api/0_interface.md +++ b/docs/src/api/0_interface.md @@ -1,3 +1,7 @@ +```@meta +CollapsedDocStrings = true +``` + # Interface ## Public diff --git a/docs/src/api/argmax.md b/docs/src/api/argmax.md index 6ea12e4..d3b8d29 100644 --- a/docs/src/api/argmax.md +++ b/docs/src/api/argmax.md @@ -1,3 +1,7 @@ +```@meta +CollapsedDocStrings = true +``` + # Argmax ## Public diff --git a/docs/src/api/argmax_2d.md b/docs/src/api/argmax_2d.md index 1b6b44e..ce28b54 100644 --- a/docs/src/api/argmax_2d.md +++ b/docs/src/api/argmax_2d.md @@ -1,3 +1,7 @@ +```@meta +CollapsedDocStrings = true +``` + # Argmax2D ## Public diff --git a/docs/src/api/dvsp.md b/docs/src/api/dvsp.md new file mode 100644 index 0000000..2922696 --- /dev/null +++ b/docs/src/api/dvsp.md @@ -0,0 +1,19 @@ +```@meta +CollapsedDocStrings = true +``` + +# Dynamic Vehicle Scheduling + +## Public + +```@autodocs +Modules = [DecisionFocusedLearningBenchmarks.DynamicVehicleScheduling] +Private = false +``` + +## Private + +```@autodocs +Modules = [DecisionFocusedLearningBenchmarks.DynamicVehicleScheduling] +Public = false +``` diff --git a/docs/src/api/dynamic_assorment.md b/docs/src/api/dynamic_assorment.md new file mode 100644 index 0000000..847d184 --- /dev/null +++ b/docs/src/api/dynamic_assorment.md @@ -0,0 +1,19 @@ +```@meta +CollapsedDocStrings = true +``` + +# Dynamic Assortment + +## Public + +```@autodocs +Modules = [DecisionFocusedLearningBenchmarks.DynamicAssortment] +Private = false +``` + +## Private + +```@autodocs +Modules = [DecisionFocusedLearningBenchmarks.DynamicAssortment] +Public = false +``` diff --git a/docs/src/api/fixed_shortest_path.md b/docs/src/api/fixed_shortest_path.md index df50a9f..36a03b2 100644 --- a/docs/src/api/fixed_shortest_path.md +++ b/docs/src/api/fixed_shortest_path.md @@ -1,3 +1,7 @@ +```@meta +CollapsedDocStrings = true +``` + # Fixed-size shortest path ## Public diff --git a/docs/src/api/portfolio_optimization.md b/docs/src/api/portfolio_optimization.md index 5b0102b..6d198ac 100644 --- a/docs/src/api/portfolio_optimization.md +++ b/docs/src/api/portfolio_optimization.md @@ -1,3 +1,7 @@ +```@meta +CollapsedDocStrings = true +``` + # Subset selection ## Public diff --git a/docs/src/api/ranking.md b/docs/src/api/ranking.md index f249a48..82d0719 100644 --- a/docs/src/api/ranking.md +++ b/docs/src/api/ranking.md @@ -1,3 +1,7 @@ +```@meta +CollapsedDocStrings = true +``` + # Ranking ## Public diff --git a/docs/src/api/subset_selection.md b/docs/src/api/subset_selection.md index 76b686d..946eb3c 100644 --- a/docs/src/api/subset_selection.md +++ b/docs/src/api/subset_selection.md @@ -1,3 +1,7 @@ +```@meta +CollapsedDocStrings = true +``` + # Subset selection ## Public diff --git a/docs/src/api/vsp.md b/docs/src/api/vsp.md index 96e4cdb..119c9ba 100644 --- a/docs/src/api/vsp.md +++ b/docs/src/api/vsp.md @@ -1,3 +1,7 @@ +```@meta +CollapsedDocStrings = true +``` + # Stochastic Vehicle Scheduling ## Public diff --git a/docs/src/api/warcraft.md b/docs/src/api/warcraft.md index 3ff6824..c3bd480 100644 --- a/docs/src/api/warcraft.md +++ b/docs/src/api/warcraft.md @@ -1,3 +1,7 @@ +```@meta +CollapsedDocStrings = true +``` + # Warcraft ## Public diff --git a/docs/src/benchmarks/argmax.md b/docs/src/benchmarks/argmax.md index e69de29..1ab74f9 100644 --- a/docs/src/benchmarks/argmax.md +++ b/docs/src/benchmarks/argmax.md @@ -0,0 +1 @@ +# Argmax diff --git a/docs/src/benchmarks/dvsp.md b/docs/src/benchmarks/dvsp.md new file mode 100644 index 0000000..2b96c67 --- /dev/null +++ b/docs/src/benchmarks/dvsp.md @@ -0,0 +1,3 @@ +# Dynamic Vehicle Scheduling + +[`DynamicVehicleSchedulingBenchmark`](@ref). diff --git a/docs/src/benchmarks/dynamic_assorment.md b/docs/src/benchmarks/dynamic_assorment.md new file mode 100644 index 0000000..dcf3243 --- /dev/null +++ b/docs/src/benchmarks/dynamic_assorment.md @@ -0,0 +1,3 @@ +# Dynamic Assortment + +[`DynamicAssortmentBenchmark`](@ref). diff --git a/docs/src/benchmarks/ranking.md b/docs/src/benchmarks/ranking.md index e69de29..5bfcaeb 100644 --- a/docs/src/benchmarks/ranking.md +++ b/docs/src/benchmarks/ranking.md @@ -0,0 +1 @@ +# Ranking diff --git a/docs/src/tutorials/warcraft.jl b/docs/src/tutorials/warcraft.jl index 13f21ba..2d41563 100644 --- a/docs/src/tutorials/warcraft.jl +++ b/docs/src/tutorials/warcraft.jl @@ -86,3 +86,6 @@ final_gap = compute_gap(b, test_dataset, model, maximizer) θ = model(x) y = maximizer(θ) plot_data(b, DataSample(; x, θ_true=θ, y_true=y)) + +using Test #src +@test final_gap < starting_gap #src diff --git a/docs/src/warcraft.md b/docs/src/warcraft.md new file mode 100644 index 0000000..c3400e7 --- /dev/null +++ b/docs/src/warcraft.md @@ -0,0 +1,155 @@ +```@meta +EditURL = "tutorials/warcraft.jl" +``` + +# Path-finding on image maps + +In this tutorial, we showcase DecisionFocusedLearningBenchmarks.jl capabilities on one of its main benchmarks: the Warcraft benchmark. +This benchmark problem is a simple path-finding problem where the goal is to find the shortest path between the top left and bottom right corners of a given image map. +The map is represented as a 2D image representing a 12x12 grid, each cell having an unknown travel cost depending on the terrain type. + +First, let's load the package and create a benchmark object as follows: + +````@example warcraft +using DecisionFocusedLearningBenchmarks +b = WarcraftBenchmark() +```` + +## Dataset generation + +These benchmark objects behave as generators that can generate various needed elements in order to build an algorithm to tackle the problem. +First of all, all benchmarks are capable of generating datasets as needed, using the [`generate_dataset`](@ref) method. +This method takes as input the benchmark object for which the dataset is to be generated, and a second argument specifying the number of samples to generate: + +````@example warcraft +dataset = generate_dataset(b, 50); +nothing #hide +```` + +We obtain a vector of [`DataSample`](@ref) objects, containing all needed data for the problem. +Subdatasets can be created through regular slicing: + +````@example warcraft +train_dataset, test_dataset = dataset[1:45], dataset[46:50] +```` + +And getting an individual sample will return a [`DataSample`](@ref) with four fields: `x`, `instance`, `θ`, and `y`: + +````@example warcraft +sample = test_dataset[1] +```` + +`x` correspond to the input features, i.e. the input image (3D array) in the Warcraft benchmark case: + +````@example warcraft +x = sample.x +```` + +`θ_true` correspond to the true unknown terrain weights. We use the opposite of the true weights in order to formulate the optimization problem as a maximization problem: + +````@example warcraft +θ_true = sample.θ_true +```` + +`y_true` correspond to the optimal shortest path, encoded as a binary matrix: + +````@example warcraft +y_true = sample.y_true +```` + +`instance` is not used in this benchmark, therefore set to nothing: + +````@example warcraft +isnothing(sample.instance) +```` + +For some benchmarks, we provide the following plotting method [`plot_data`](@ref) to visualize the data: + +````@example warcraft +plot_data(b, sample) +```` + +We can see here the terrain image, the true terrain weights, and the true shortest path avoiding the high cost cells. + +## Building a pipeline + +DecisionFocusedLearningBenchmarks also provides methods to build an hybrid machine learning and combinatorial optimization pipeline for the benchmark. +First, the [`generate_statistical_model`](@ref) method generates a machine learning predictor to predict cell weights from the input image: + +````@example warcraft +model = generate_statistical_model(b) +```` + +In the case of the Warcraft benchmark, the model is a convolutional neural network built using the Flux.jl package. + +````@example warcraft +θ = model(x) +```` + +Note that the model is not trained yet, and its parameters are randomly initialized. + +Finally, the [`generate_maximizer`](@ref) method can be used to generate a combinatorial optimization algorithm that takes the predicted cell weights as input and returns the corresponding shortest path: + +````@example warcraft +maximizer = generate_maximizer(b; dijkstra=true) +```` + +In the case o fthe Warcraft benchmark, the method has an additional keyword argument to chose the algorithm to use: Dijkstra's algorithm or Bellman-Ford algorithm. + +````@example warcraft +y = maximizer(θ) +```` + +As we can see, currently the pipeline predicts random noise as cell weights, and therefore the maximizer returns a straight line path. + +````@example warcraft +plot_data(b, DataSample(; x, θ_true=θ, y_true=y)) +```` + +We can evaluate the current pipeline performance using the optimality gap metric: + +````@example warcraft +starting_gap = compute_gap(b, test_dataset, model, maximizer) +```` + +## Using a learning algorithm + +We can now train the model using the InferOpt.jl package: + +````@example warcraft +using InferOpt +using Flux +using Plots + +perturbed_maximizer = PerturbedMultiplicative(maximizer; ε=0.2, nb_samples=100) +loss = FenchelYoungLoss(perturbed_maximizer) + +starting_gap = compute_gap(b, test_dataset, model, maximizer) + +opt_state = Flux.setup(Adam(1e-3), model) +loss_history = Float64[] +for epoch in 1:50 + val, grads = Flux.withgradient(model) do m + sum(loss(m(x), y_true) for (; x, y_true) in train_dataset) / length(train_dataset) + end + Flux.update!(opt_state, model, grads[1]) + push!(loss_history, val) +end + +plot(loss_history; xlabel="Epoch", ylabel="Loss", title="Training loss") +```` + +````@example warcraft +final_gap = compute_gap(b, test_dataset, model, maximizer) +```` + +````@example warcraft +θ = model(x) +y = maximizer(θ) +plot_data(b, DataSample(; x, θ_true=θ, y_true=y)) +```` + +--- + +*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* + diff --git a/src/Argmax/Argmax.jl b/src/Argmax/Argmax.jl index c4d07bf..fa6ddda 100644 --- a/src/Argmax/Argmax.jl +++ b/src/Argmax/Argmax.jl @@ -67,22 +67,18 @@ end """ $TYPEDSIGNATURES -Generate a dataset of labeled instances for the argmax problem. +Generate a data sample for the argmax benchmark. +This function generates a random feature matrix, computes the costs using the encoder, +and adds noise to the costs before computing a target solution. """ -function Utils.generate_dataset( - bench::ArgmaxBenchmark, dataset_size::Int=10; seed::Int=0, noise_std=0.0 +function Utils.generate_sample( + bench::ArgmaxBenchmark, rng::AbstractRNG; noise_std::Float32=0.0f0 ) (; instance_dim, nb_features, encoder) = bench - rng = MersenneTwister(seed) - features = [randn(rng, Float32, nb_features, instance_dim) for _ in 1:dataset_size] - costs = encoder.(features) - noisy_solutions = [ - one_hot_argmax(θ + noise_std * randn(rng, Float32, instance_dim)) for θ in costs - ] - return [ - DataSample(; x, θ_true, y_true) for - (x, θ_true, y_true) in zip(features, costs, noisy_solutions) - ] + features = randn(rng, Float32, nb_features, instance_dim) + costs = encoder(features) + noisy_solution = one_hot_argmax(costs + noise_std * randn(rng, Float32, instance_dim)) + return DataSample(; x=features, θ_true=costs, y_true=noisy_solution) end """ diff --git a/src/Argmax2D/Argmax2D.jl b/src/Argmax2D/Argmax2D.jl index a1c76c4..169c403 100644 --- a/src/Argmax2D/Argmax2D.jl +++ b/src/Argmax2D/Argmax2D.jl @@ -7,7 +7,7 @@ using Flux: Chain, Dense using LaTeXStrings: @L_str using LinearAlgebra: dot, norm using Plots: Plots -using Random: Random, MersenneTwister +using Random: Random, MersenneTwister, AbstractRNG include("polytope.jl") @@ -53,20 +53,16 @@ maximizer(θ; instance, kwargs...) = instance[argmax(dot(θ, v) for v in instanc """ $TYPEDSIGNATURES -Generate a dataset for the [`Argmax2DBenchmark`](@ref). +Generate a sample for the [`Argmax2DBenchmark`](@ref). """ -function Utils.generate_dataset( - bench::Argmax2DBenchmark, dataset_size=10; seed=nothing, rng=MersenneTwister(seed) -) +function Utils.generate_sample(bench::Argmax2DBenchmark, rng::AbstractRNG) (; nb_features, encoder, polytope_vertex_range) = bench - return map(1:dataset_size) do _ - x = randn(rng, Float32, nb_features) - θ_true = encoder(x) - θ_true ./= 2 * norm(θ_true) - instance = build_polytope(rand(rng, polytope_vertex_range); shift=rand(rng)) - y_true = maximizer(θ_true; instance) - return DataSample(; x=x, θ_true=θ_true, y_true=y_true, instance=instance) - end + x = randn(rng, Float32, nb_features) + θ_true = encoder(x) + θ_true ./= 2 * norm(θ_true) + instance = build_polytope(rand(rng, polytope_vertex_range); shift=rand(rng)) + y_true = maximizer(θ_true; instance) + return DataSample(; x=x, θ_true=θ_true, y_true=y_true, instance=instance) end """ diff --git a/src/DecisionFocusedLearningBenchmarks.jl b/src/DecisionFocusedLearningBenchmarks.jl index 0a6b168..be2c500 100644 --- a/src/DecisionFocusedLearningBenchmarks.jl +++ b/src/DecisionFocusedLearningBenchmarks.jl @@ -4,17 +4,41 @@ using DataDeps using Requires: @require function __init__() + function _euro_neurips_unpack(local_filepath) + directory = dirname(local_filepath) + unpack(local_filepath) + # Move instances and delete the rest + for filepath in readdir( + joinpath(directory, "euro-neurips-vrp-2022-quickstart-main", "instances"); + join=true, + ) + if endswith(filepath, ".txt") + mv(filepath, joinpath(directory, basename(filepath))) + end + end + rm(joinpath(directory, "euro-neurips-vrp-2022-quickstart-main"); recursive=true) + return nothing + end # Register the Warcraft dataset ENV["DATADEPS_ALWAYS_ACCEPT"] = "true" register( DataDep( "warcraft", - "This is the warcraft dataset", + "Warcraft shortest path dataset", "http://cermics.enpc.fr/~bouvierl/warcraft_TP/data.zip"; post_fetch_method=unpack, ), ) + register( + DataDep( + "dvrptw", + "EURO-NeurIPS challenge 2022 dataset for the dynamic vehicle routing problem with time windows", + "https://github.com/ortec/euro-neurips-vrp-2022-quickstart/archive/refs/heads/main.zip"; + post_fetch_method=_euro_neurips_unpack, + ), + ) + # Gurobi setup @info "If you have Gurobi installed and want to use it, make sure to `using Gurobi` in order to enable it." @require Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b" include("gurobi_setup.jl") @@ -31,34 +55,50 @@ include("Warcraft/Warcraft.jl") include("FixedSizeShortestPath/FixedSizeShortestPath.jl") include("PortfolioOptimization/PortfolioOptimization.jl") include("StochasticVehicleScheduling/StochasticVehicleScheduling.jl") +include("DynamicVehicleScheduling/DynamicVehicleScheduling.jl") +include("DynamicAssortment/DynamicAssortment.jl") using .Utils -using .Argmax -using .Argmax2D -using .Ranking -using .SubsetSelection -using .Warcraft -using .FixedSizeShortestPath -using .PortfolioOptimization -using .StochasticVehicleScheduling # Interface -export AbstractBenchmark, DataSample -export generate_dataset +export AbstractBenchmark, AbstractStochasticBenchmark, AbstractDynamicBenchmark, DataSample +export AbstractEnvironment, get_seed, is_terminated, observe, reset!, step! + +export Policy, evaluate_policy! + +export generate_sample, generate_dataset, generate_environments, generate_environment +export generate_scenario +export generate_policies export generate_statistical_model export generate_maximizer, maximizer_kwargs +export generate_anticipative_solution +export is_exogenous, is_endogenous + export objective_value export plot_data, plot_instance, plot_solution export compute_gap # Export all benchmarks -export ArgmaxBenchmark +using .Argmax +using .Argmax2D +using .Ranking +using .SubsetSelection +using .Warcraft +using .FixedSizeShortestPath +using .PortfolioOptimization +using .StochasticVehicleScheduling +using .DynamicVehicleScheduling +using .DynamicAssortment + export Argmax2DBenchmark -export RankingBenchmark -export SubsetSelectionBenchmark -export WarcraftBenchmark +export ArgmaxBenchmark +export DynamicAssortmentBenchmark +export DynamicVehicleSchedulingBenchmark export FixedSizeShortestPathBenchmark export PortfolioOptimizationBenchmark +export RankingBenchmark export StochasticVehicleSchedulingBenchmark +export SubsetSelectionBenchmark +export WarcraftBenchmark end # module DecisionFocusedLearningBenchmarks diff --git a/src/DynamicAssortment/DynamicAssortment.jl b/src/DynamicAssortment/DynamicAssortment.jl new file mode 100644 index 0000000..2c61c5f --- /dev/null +++ b/src/DynamicAssortment/DynamicAssortment.jl @@ -0,0 +1,97 @@ +module DynamicAssortment + +using ..Utils + +using DocStringExtensions: TYPEDEF, TYPEDFIELDS, TYPEDSIGNATURES +using Distributions: Uniform, Categorical +using LinearAlgebra: dot +using Random: Random, AbstractRNG, MersenneTwister +using Statistics: mean + +using Flux: Chain, Dense +using Combinatorics: combinations + +""" +$TYPEDEF + +Benchmark for the dynamic assortment problem. + +# Fields +$TYPEDFIELDS +""" +struct DynamicAssortmentBenchmark{exogenous,M} <: AbstractDynamicBenchmark{exogenous} + "customer choice model (price, hype, saturation, and features)" + customer_choice_model::M + "number of items" + N::Int + "dimension of feature vectors (in addition to hype, satisfaction, and price)" + d::Int + "assortment size constraint" + K::Int + "number of steps per episode" + max_steps::Int +end + +function DynamicAssortmentBenchmark(; + N=20, + d=2, + K=4, + max_steps=80, + customer_choice_model=Chain(Dense([-0.8 0.6 -0.4 0.3 0.5]), vec), + exogenous=false, +) + return DynamicAssortmentBenchmark{exogenous,typeof(customer_choice_model)}( + customer_choice_model, N, d, K, max_steps + ) +end + +include("instance.jl") +include("environment.jl") +include("policies.jl") + +customer_choice_model(b::DynamicAssortmentBenchmark) = b.customer_choice_model +item_count(b::DynamicAssortmentBenchmark) = b.N +feature_count(b::DynamicAssortmentBenchmark) = b.d +assortment_size(b::DynamicAssortmentBenchmark) = b.K +max_steps(b::DynamicAssortmentBenchmark) = b.max_steps + +function Utils.generate_sample( + b::DynamicAssortmentBenchmark, rng::AbstractRNG=MersenneTwister(0) +) + return DataSample(; instance=Instance(b, rng)) +end + +function Utils.generate_statistical_model(b::DynamicAssortmentBenchmark; seed=nothing) + Random.seed!(seed) + d = feature_count(b) + return Chain(Dense(d + 8 => 5), Dense(5 => 1), vec) +end + +function Utils.generate_maximizer(b::DynamicAssortmentBenchmark) + return TopKMaximizer(assortment_size(b)) +end + +function Utils.generate_environment( + ::DynamicAssortmentBenchmark, instance::Instance, rng::AbstractRNG; kwargs... +) + seed = rand(rng, 1:typemax(Int)) + return Environment(instance; seed) +end + +function Utils.generate_policies(b::DynamicAssortmentBenchmark) + greedy = Policy( + "Greedy", + "policy that selects the assortment with items with the highest prices", + greedy_policy, + ) + expert = Policy( + "Expert", + "policy that selects the assortment with the highest expected revenue", + expert_policy, + ) + return (expert, greedy) +end + +export DynamicAssortmentBenchmark + +end diff --git a/src/DynamicAssortment/environment.jl b/src/DynamicAssortment/environment.jl new file mode 100644 index 0000000..8389a0c --- /dev/null +++ b/src/DynamicAssortment/environment.jl @@ -0,0 +1,169 @@ +""" +$TYPEDEF + +Environment for the dynamic assortment problem. + +# Fields +$TYPEDFIELDS +""" +@kwdef mutable struct Environment{I<:Instance,R<:AbstractRNG,S<:Union{Nothing,Int}} <: + Utils.AbstractEnvironment + "associated instance" + instance::I + "current step" + step::Int + "purchase history (used to update hype feature)" + purchase_hist::Vector{Int} + "rng" + rng::R + "seed for RNG" + seed::S + "customer utility for each item" + utility::Vector{Float64} + "current full features" + features::Matrix{Float64} + "satisfaction + hype feature change from the last step" + d_features::Matrix{Float64} +end + +function Environment(instance::Instance; seed=0, rng::AbstractRNG=MersenneTwister(seed)) + N = item_count(instance) + (; prices, features, starting_hype_and_saturation) = instance + full_features = vcat( + reshape(prices[1:(end - 1)], 1, :), starting_hype_and_saturation, features + ) + model = customer_choice_model(instance) + env = Environment(; + instance, + step=1, + purchase_hist=Int[], + rng=rng, + seed=seed, + utility=model(full_features), + features=full_features, + d_features=zeros(2, N), + ) + Utils.reset!(env; reset_rng=true) + return env +end + +Utils.get_seed(env::Environment) = env.seed +customer_choice_model(b::Environment) = customer_choice_model(b.instance) +item_count(b::Environment) = item_count(b.instance) +feature_count(b::Environment) = feature_count(b.instance) +assortment_size(b::Environment) = assortment_size(b.instance) +max_steps(b::Environment) = max_steps(b.instance) +prices(b::Environment) = b.instance.prices + +## Basic operations of environment + +# Reset the environment +function Utils.reset!(env::Environment; reset_rng=false, seed=env.seed) + reset_rng && Random.seed!(env.rng, seed) + + env.step = 1 + + (; prices, starting_hype_and_saturation, features) = env.instance + features = vcat( + reshape(prices[1:(end - 1)], 1, :), starting_hype_and_saturation, features + ) + env.features .= features + + env.d_features .= 0.0 + + model = customer_choice_model(env) + env.utility .= model(features) + + empty!(env.purchase_hist) + return nothing +end + +function Utils.is_terminated(env::Environment) + return env.step > max_steps(env) +end + +function Utils.observe(env::Environment) + delta_features = env.features[2:3, :] .- env.instance.starting_hype_and_saturation + return vcat( + env.features, + env.d_features, + delta_features, + ones(1, item_count(env)) .* (env.step / max_steps(env) * 10), + ) ./ 10, + nothing +end + +# Compute the hype vector +function hype_update(env::Environment) + N = item_count(env) + hype_vector = ones(N) + hist = env.purchase_hist + + # Define decay factors for each time step + factors = [0.02, -0.005, -0.005, -0.005, -0.005] + + # Apply updates for the last 5 purchases + for (i, factor) in enumerate(factors) + if length(hist) >= i + item = hist[end - i + 1] + if item <= N + hype_vector[item] += factor + end + end + end + + return hype_vector +end + +# Step function +function buy_item!(env::Environment, item::Int) + push!(env.purchase_hist, item) + env.step += 1 + + if is_endogenous(env.instance.config) + old_features = copy(env.features[2:3, :]) + # update hype feature + hype_vector = hype_update(env) + env.features[2, :] .*= hype_vector + + # update saturation feature + if item <= item_count(env) + env.features[3, item] *= 1.01 + end + + env.utility .= customer_choice_model(env)(env.features) + env.d_features = env.features[2:3, :] - old_features + end + return nothing +end + +# Choice probabilities +function choice_probabilities(env::Environment, S) + N = item_count(env) + θ = env.utility + exp_values = [exp(θ[i]) * S[i] for i in 1:N] + push!(exp_values, 1.0) # No purchase action + denominator = sum(exp_values) + probs = exp_values ./ denominator + return probs +end + +# Purchase decision +function Utils.step!(env::Environment, assortment) + @assert !Utils.is_terminated(env) "Environment is terminated, cannot act!" + r = prices(env) + probs = choice_probabilities(env, assortment) + item = rand(env.rng, Categorical(probs)) + reward = r[item] + buy_item!(env, item) + return reward +end + +## Solution functions +# enumerate all possible assortments of size K and return the best one +function compute_expected_revenue(env::Environment, S) + r = prices(env) + probs = choice_probabilities(env, S) + expected_revenue = dot(probs, r) + return expected_revenue +end diff --git a/src/DynamicAssortment/instance.jl b/src/DynamicAssortment/instance.jl new file mode 100644 index 0000000..3250cdd --- /dev/null +++ b/src/DynamicAssortment/instance.jl @@ -0,0 +1,33 @@ +""" +$TYPEDEF + +Instance of the dynamic assortment problem. + +# Fields +$TYPEDFIELDS +""" +@kwdef struct Instance{B<:DynamicAssortmentBenchmark} + "associated benchmark" + config::B + "item prices (including no purchase action)" + prices::Vector{Float64} + "static features, size (d, N)" + features::Matrix{Float64} + "starting hype and saturation features, size (2, N)" + starting_hype_and_saturation::Matrix{Float64} +end + +function Instance(b::DynamicAssortmentBenchmark, rng::AbstractRNG) + N = item_count(b) + d = feature_count(b) + prices = vcat(rand(rng, Uniform(1.0, 10.0), N), 0.0) # last price is for no purchase action + features = rand(rng, Uniform(1.0, 10.0), (d, N)) + starting_hype_and_saturation = rand(rng, Uniform(1.0, 10.0), (2, N)) + return Instance(; config=b, prices, features, starting_hype_and_saturation) +end + +customer_choice_model(b::Instance) = customer_choice_model(b.config) +item_count(b::Instance) = item_count(b.config) +feature_count(b::Instance) = feature_count(b.config) +assortment_size(b::Instance) = assortment_size(b.config) +max_steps(b::Instance) = max_steps(b.config) diff --git a/src/DynamicAssortment/policies.jl b/src/DynamicAssortment/policies.jl new file mode 100644 index 0000000..320c501 --- /dev/null +++ b/src/DynamicAssortment/policies.jl @@ -0,0 +1,21 @@ +function expert_policy(env::Environment) + N = item_count(env) + K = assortment_size(env) + best_S = falses(N) + best_revenue = -1.0 + S_vec = falses(N) + for S in combinations(1:N, K) + S_vec .= false + S_vec[S] .= true + expected_revenue = compute_expected_revenue(env, S_vec) + if expected_revenue > best_revenue + best_S, best_revenue = copy(S_vec), expected_revenue + end + end + return best_S +end + +function greedy_policy(env::Environment) + maximizer = generate_maximizer(env.instance.config) + return maximizer(prices(env)) +end diff --git a/src/DynamicVehicleScheduling/DynamicVehicleScheduling.jl b/src/DynamicVehicleScheduling/DynamicVehicleScheduling.jl new file mode 100644 index 0000000..7421032 --- /dev/null +++ b/src/DynamicVehicleScheduling/DynamicVehicleScheduling.jl @@ -0,0 +1,120 @@ +module DynamicVehicleScheduling + +using ..Utils + +using Base: @kwdef +using DataDeps: @datadep_str +using DocStringExtensions: TYPEDEF, TYPEDFIELDS, TYPEDSIGNATURES +using Flux: Chain, Dense +using Graphs +using HiGHS +using InferOpt: LinearMaximizer +using IterTools: partition +using JSON +using JuMP +using Plots: plot, plot!, scatter! +using Printf: @printf +using Random: Random, AbstractRNG, MersenneTwister, seed!, randperm +using Requires: @require +using Statistics: mean, quantile + +include("utils.jl") + +# static vsp stuff +include("static_vsp/instance.jl") +include("static_vsp/parsing.jl") +include("static_vsp/solution.jl") +include("static_vsp/plot.jl") + +include("instance.jl") +include("state.jl") +include("scenario.jl") +include("environment.jl") +include("plot.jl") + +include("maximizer.jl") +include("anticipative_solver.jl") + +include("features.jl") +include("policy.jl") + +""" +$TYPEDEF + +Abstract type for dynamic vehicle scheduling benchmarks. + +# Fields +$TYPEDFIELDS +""" +@kwdef struct DynamicVehicleSchedulingBenchmark <: AbstractDynamicBenchmark{true} + "maximum number of customers entering the system per epoch" + max_requests_per_epoch::Int = 10 + "time between decision and dispatch of a vehicle" + Δ_dispatch::Float64 = 1.0 + "duration of an epoch" + epoch_duration::Float64 = 1.0 + "whether to use two-dimensional features" + two_dimensional_features::Bool = false +end + +function Utils.generate_dataset(b::DynamicVehicleSchedulingBenchmark, dataset_size::Int=1) + (; max_requests_per_epoch, Δ_dispatch, epoch_duration, two_dimensional_features) = b + files = readdir(datadep"dvrptw"; join=true) + dataset_size = min(dataset_size, length(files)) + return [ + DataSample(; + instance=Instance( + read_vsp_instance(files[i]); + max_requests_per_epoch, + Δ_dispatch, + epoch_duration, + two_dimensional_features, + ), + ) for i in 1:dataset_size + ] +end + +function Utils.generate_environment( + ::DynamicVehicleSchedulingBenchmark, instance::Instance, rng::AbstractRNG; kwargs... +) + seed = rand(rng, 1:typemax(Int)) + return DVSPEnv(instance; seed) +end + +function Utils.generate_maximizer(::DynamicVehicleSchedulingBenchmark) + return LinearMaximizer(oracle; g, h) +end + +function Utils.generate_scenario(b::DynamicVehicleSchedulingBenchmark, args...; kwargs...) + return Utils.generate_scenario(args...; kwargs...) +end + +function Utils.generate_anticipative_solution( + b::DynamicVehicleSchedulingBenchmark, args...; kwargs... +) + return anticipative_solver( + args...; kwargs..., two_dimensional_features=b.two_dimensional_features + ) +end + +function Utils.generate_policies(b::DynamicVehicleSchedulingBenchmark) + lazy = Policy( + "Lazy", + "Lazy policy that dispatches vehicles only when they are ready.", + lazy_policy, + ) + greedy = Policy( + "Greedy", + "Greedy policy that dispatches vehicles to the nearest customer.", + greedy_policy, + ) + return (lazy, greedy) +end + +function Utils.generate_statistical_model(b::DynamicVehicleSchedulingBenchmark) + return Chain(Dense((b.two_dimensional_features ? 2 : 14) => 1), vec) +end + +export DynamicVehicleSchedulingBenchmark + +end diff --git a/src/DynamicVehicleScheduling/anticipative_solver.jl b/src/DynamicVehicleScheduling/anticipative_solver.jl new file mode 100644 index 0000000..5847808 --- /dev/null +++ b/src/DynamicVehicleScheduling/anticipative_solver.jl @@ -0,0 +1,213 @@ +""" +$TYPEDSIGNATURES + +Retrieve anticipative routes solution from the given MIP solution `y`. +Outputs a set of routes per epoch. +""" +function retrieve_routes_anticipative( + y::AbstractArray, dvspenv::DVSPEnv, customer_index, epoch_indices +) + nb_tasks = length(customer_index) + # first_epoch = 1 + # (; last_epoch) = dvspenv.instance + job_indices = 2:(nb_tasks) + # epoch_indices = first_epoch:last_epoch + + routes = [Vector{Int}[] for _ in epoch_indices] + for (i, t) in enumerate(epoch_indices) + start = [i for i in job_indices if y[1, i, t] ≈ 1] + for task in start + route = Int[] + current_task = task + while current_task != 1 # < nb_tasks + push!(route, current_task) + local next_task + for i in 1:nb_tasks + if isapprox(y[current_task, i, t], 1; atol=0.1) + next_task = i + break + end + end + current_task = next_task + end + push!(routes[i], route) + end + end + return routes +end + +""" +$TYPEDSIGNATURES + +Solve the anticipative VSP problem for environment `env`. +For this, it uses the current environment history, so make sure that the environment is terminated before calling this method. +""" +function anticipative_solver( + env::DVSPEnv, + scenario=env.scenario; + model_builder=highs_model, + two_dimensional_features=env.instance.two_dimensional_features, + reset_env=true, + nb_epochs=typemax(Int), + seed=get_seed(env), +) + reset_env && reset!(env; reset_rng=true, seed) + + start_epoch = current_epoch(env) + end_epoch = min(last_epoch(env), start_epoch + nb_epochs - 1) + T = start_epoch:end_epoch + + request_epoch = [0] + for t in T + request_epoch = vcat(request_epoch, fill(t, length(scenario.indices[t]))) + end + customer_index = vcat(1, scenario.indices[T]...) + service_time = vcat(0.0, scenario.service_time[T]...) + start_time = vcat(0.0, scenario.start_time[T]...) + + duration = env.instance.static_instance.duration[customer_index, customer_index] + (; epoch_duration, Δ_dispatch) = env.instance + + model = model_builder() + set_silent(model) + + nb_nodes = length(customer_index) + job_indices = 2:nb_nodes + epoch_indices = T#first_epoch:last_epoch + + @variable(model, y[i=1:nb_nodes, j=1:nb_nodes, t=epoch_indices]; binary=true) + + @objective( + model, + Max, + sum( + -duration[i, j] * y[i, j, t] for i in 1:nb_nodes, j in 1:nb_nodes, + t in epoch_indices + ) + ) + + # flow constraint per epoch + for t in epoch_indices, i in 1:nb_nodes + @constraint( + model, + sum(y[j, i, t] for j in 1:nb_nodes) == sum(y[i, j, t] for j in 1:nb_nodes) + ) + end + + # each task must be done once along the horizon + @constraint( + model, + demand[i in job_indices], + sum(y[j, i, t] for j in 1:nb_nodes, t in epoch_indices) == 1 + ) + + # a trip from i can be planned only after request appeared (release times) + for i in job_indices, t in epoch_indices, j in 1:nb_nodes + if t < request_epoch[i] + @constraint(model, y[i, j, t] <= 0) + end + end + + # a trip from i can be done only before limit date + for i in job_indices, t in epoch_indices, j in 1:nb_nodes + if (t - 1) * epoch_duration + duration[1, i] + Δ_dispatch > start_time[i] + @constraint(model, y[i, j, t] <= 0) + end + end + + # trips can be planned if start, service and transport times enable it + for i in job_indices, t in epoch_indices, j in job_indices + if start_time[i] <= start_time[j] + if start_time[i] + service_time[i] + duration[i, j] > start_time[j] + @constraint(model, y[i, j, t] <= 0) + end + else + @constraint(model, y[i, j, t] <= 0) + end + end + + optimize!(model) + + obj = JuMP.objective_value(model) + epoch_routes = retrieve_routes_anticipative( + value.(y), env, customer_index, epoch_indices + ) + + epoch_indices = Vector{Int}[] + N = 1 + indices = [1] + index = 1 + for epoch in 1:last_epoch(env) + M = length(scenario.indices[epoch]) + indices = vcat(indices, (N + 1):(N + M)) + push!(epoch_indices, copy(indices)) + N = N + M + if epoch in T + dispatched = vcat(epoch_routes[index]...) + index += 1 + indices = setdiff(indices, dispatched) + end + end + + indices = vcat(1, scenario.indices...) + start_time = vcat(0.0, scenario.start_time...) + service_time = vcat(0.0, scenario.service_time...) + + dataset = map(enumerate(T)) do (i, epoch) + routes = epoch_routes[i] + epoch_customers = epoch_indices[epoch] + + y_true = + VSPSolution( + Vector{Int}[ + map(idx -> findfirst(==(idx), epoch_customers), route) for + route in routes + ]; + max_index=length(epoch_customers), + ).edge_matrix + + location_indices = indices[epoch_customers] + new_coordinates = env.instance.static_instance.coordinate[location_indices] + new_start_time = start_time[epoch_customers] + new_service_time = service_time[epoch_customers] + new_duration = env.instance.static_instance.duration[ + location_indices, location_indices + ] + static_instance = StaticInstance( + new_coordinates, new_service_time, new_start_time, new_duration + ) + + is_must_dispatch = falses(length(location_indices)) + is_postponable = falses(length(location_indices)) + + epoch_duration = env.instance.epoch_duration + Δ_dispatch = env.instance.Δ_dispatch + planning_start_time = (epoch - 1) * epoch_duration + Δ_dispatch + if epoch == last_epoch + # If we are in the last epoch, all requests must be dispatched + is_must_dispatch[2:end] .= true + else + is_must_dispatch[2:end] .= + planning_start_time .+ epoch_duration .+ @view(new_duration[1, 2:end]) .> new_start_time[2:end] + end + is_postponable[2:end] .= .!is_must_dispatch[2:end] + + state = DVSPState(; + state_instance=static_instance, + is_must_dispatch, + is_postponable, + location_indices, + current_epoch=epoch, + ) + + x = if two_dimensional_features + compute_2D_features(state, env.instance) + else + compute_features(state, env.instance) + end + + return DataSample(; instance=state, y_true, x) + end + + return obj, dataset +end diff --git a/src/DynamicVehicleScheduling/environment.jl b/src/DynamicVehicleScheduling/environment.jl new file mode 100644 index 0000000..339bd78 --- /dev/null +++ b/src/DynamicVehicleScheduling/environment.jl @@ -0,0 +1,100 @@ +mutable struct DVSPEnv{S<:DVSPState,R<:AbstractRNG,SS} <: Utils.AbstractEnvironment + "associated instance" + instance::Instance + "current state" + state::S + "scenario the environment will use when not given a specific one" + scenario::Scenario + "random number generator" + rng::R + "seed for the environment" + seed::SS +end + +""" +$TYPEDSIGNATURES + +Constructor for [`DVSPEnv`](@ref). +""" +function DVSPEnv(instance::Instance; seed=nothing) + rng = MersenneTwister(seed) + scenario = Utils.generate_scenario(instance; rng) + initial_state = DVSPState(instance; scenario[1]...) + return DVSPEnv(instance, initial_state, scenario, rng, seed) +end + +currrent_epoch(env::DVSPEnv) = current_epoch(env.state) +epoch_duration(env::DVSPEnv) = epoch_duration(env.instance) +last_epoch(env::DVSPEnv) = last_epoch(env.instance) +Δ_dispatch(env::DVSPEnv) = Δ_dispatch(env.instance) + +Utils.get_seed(env::DVSPEnv) = env.seed + +""" +$TYPEDSIGNATURES + +Get the current state of the environment. +""" +function Utils.observe(env::DVSPEnv) + if env.instance.two_dimensional_features + return compute_2D_features(env.state, env.instance), env.state + end + # else + return compute_features(env.state, env.instance), env.state +end + +current_epoch(env::DVSPEnv) = current_epoch(env.state) + +""" +$TYPEDSIGNATURES + +Get the current time of the environment, i.e. the start time of the current_epoch. +""" +time(env::DVSPEnv) = (current_epoch(env) - 1) * epoch_duration(env) + +""" +$TYPEDSIGNATURES + +Get the planning start time of the environment, i.e. the time at which vehicles routes dispatched in current epoch can depart. +""" +planning_start_time(env::DVSPEnv) = time(env) + Δ_dispatch(env) + +""" +$TYPEDSIGNATURES + +Check if the episode is terminated, i.e. if the current epoch is the last one. +""" +Utils.is_terminated(env::DVSPEnv) = current_epoch(env) > last_epoch(env) + +""" +$TYPEDSIGNATURES + +Reset the environment to its initial state. +Also reset the rng to `seed` if `reset_rng` is set to true. +""" +function Utils.reset!(env::DVSPEnv; seed=get_seed(env), reset_rng=false) + if reset_rng + Random.seed!(env.rng, seed) + end + env.scenario = Utils.generate_scenario(env; rng=env.rng) + reset_state!(env.state, env.instance; env.scenario[1]...) + return nothing +end + +""" +$TYPEDSIGNATURES + +Remove dispatched customers, advance time, and add new requests to the environment. +""" +function Utils.step!(env::DVSPEnv, routes, scenario=env.scenario) + reward = -apply_routes!(env.state, routes) + env.state.current_epoch += 1 + if !Utils.is_terminated(env) + add_new_customers!(env.state, env.instance; scenario[current_epoch(env)]...) + end + return reward +end + +function Utils.generate_scenario(env::DVSPEnv; kwargs...) + return generate_scenario(env.instance; kwargs...) +end diff --git a/src/DynamicVehicleScheduling/features.jl b/src/DynamicVehicleScheduling/features.jl new file mode 100644 index 0000000..10e0ab8 --- /dev/null +++ b/src/DynamicVehicleScheduling/features.jl @@ -0,0 +1,59 @@ +function get_features_quantileTimeToRequests(state::DVSPState, instance::Instance) + quantiles = [i * 0.1 for i in 1:9] + a = instance.static_instance.duration[state.location_indices, 2:end] + quantileTimeToRequests = mapslices(x -> quantile(x, quantiles), a; dims=2) + return quantileTimeToRequests +end + +function compute_model_free_features(state::DVSPState, instance::Instance) + (; state_instance, is_postponable) = state + + startTimes = state_instance.start_time + endTimes = startTimes .+ state_instance.service_time + timeDepotRequest = state_instance.duration[:, 1] + timeRequestDepot = state_instance.duration[1, :] + + slack_next_epoch = startTimes .- instance.epoch_duration + + model_free_features = hcat( + startTimes[is_postponable], # 1 + endTimes[is_postponable], # 2 + timeDepotRequest[is_postponable], # 3 + timeRequestDepot[is_postponable], # 4 + slack_next_epoch[is_postponable], # 5-14 + ) + return model_free_features +end + +function compute_model_aware_features(state::DVSPState, instance::Instance) + quantileTimeToRequests = get_features_quantileTimeToRequests(state, instance) + model_aware_features = quantileTimeToRequests + return model_aware_features[state.is_postponable, :] +end + +function compute_features(state::DVSPState, instance::Instance) + model_free_features = compute_model_free_features(state, instance) + model_aware_features = compute_model_aware_features(state, instance) + return hcat(model_free_features, model_aware_features)' +end + +function compute_features(env::DVSPEnv) + return compute_features(env.state, env.instance) +end + +function get_features_meanTimeToRequests(state::DVSPState, instance::Instance) + quantiles = [0.5] + a = instance.static_instance.duration[state.location_indices, 2:end] + quantileTimeToRequests = mapslices(x -> quantile(x, quantiles), a; dims=2) + return quantileTimeToRequests +end + +function compute_2D_features(state::DVSPState, instance::Instance) + timeDepotRequest = state.state_instance.duration[:, 1][state.is_postponable] + quantileTimeToRequests = get_features_meanTimeToRequests(state, instance)[state.is_postponable] + return hcat(timeDepotRequest, quantileTimeToRequests)' +end + +function compute_2D_features(env::DVSPEnv) + return compute_2D_features(env.state, env.instance) +end diff --git a/src/DynamicVehicleScheduling/instance.jl b/src/DynamicVehicleScheduling/instance.jl new file mode 100644 index 0000000..d65010c --- /dev/null +++ b/src/DynamicVehicleScheduling/instance.jl @@ -0,0 +1,48 @@ +""" +$TYPEDEF + +Instance data structure for the dynamic vehicle scheduling problem. +""" +@kwdef struct Instance{I<:StaticInstance,T} + "static instance to sample arriving requests from" + static_instance::I + "max number of new requests per epoch (rejection sampling)" + max_requests_per_epoch::Int = 10 + "time distance between epoch start and routes start" + Δ_dispatch::T = 1.0 + "duration of each epoch" + epoch_duration::T = 1.0 + "last epoch index" + last_epoch::Int + "whether to use two-dimensional features" + two_dimensional_features::Bool = false +end + +function Instance( + static_instance::StaticInstance; + max_requests_per_epoch::Int=10, + Δ_dispatch::Float64=1.0, + epoch_duration::Float64=1.0, + two_dimensional_features::Bool=false, +) + last_epoch = trunc( + Int, + ( + maximum(static_instance.start_time) - minimum(static_instance.duration[1, :]) - + Δ_dispatch + ) / epoch_duration, + ) + return Instance(; + static_instance=static_instance, + max_requests_per_epoch=max_requests_per_epoch, + Δ_dispatch=Δ_dispatch, + epoch_duration=epoch_duration, + last_epoch=last_epoch, + two_dimensional_features=two_dimensional_features, + ) +end + +Δ_dispatch(instance::Instance) = instance.Δ_dispatch +epoch_duration(instance::Instance) = instance.epoch_duration +last_epoch(instance::Instance) = instance.last_epoch +max_requests_per_epoch(instance::Instance) = instance.max_requests_per_epoch diff --git a/src/DynamicVehicleScheduling/maximizer.jl b/src/DynamicVehicleScheduling/maximizer.jl new file mode 100644 index 0000000..450ab8a --- /dev/null +++ b/src/DynamicVehicleScheduling/maximizer.jl @@ -0,0 +1,154 @@ +""" +$TYPEDSIGNATURES + +Create the acyclic digraph associated with the given VSP `instance`. +""" +function create_graph(instance::StaticInstance) + (; duration, start_time, service_time) = instance + # Initialize directed graph + nb_vertices = location_count(instance) + graph = SimpleDiGraph(nb_vertices) + + depot = 1 # depot is always index 1 + customers = 2:nb_vertices # other vertices are customers + + # Create existing edges + for i₁ in customers + # link every task to depot + add_edge!(graph, depot, i₁) + add_edge!(graph, i₁, depot) + + t₁ = start_time[i₁] + for i₂ in (i₁ + 1):nb_vertices + t₂ = start_time[i₂] + + if t₁ <= t₂ + if t₁ + service_time[i₁] + duration[i₁, i₂] <= t₂ + add_edge!(graph, i₁, i₂) + end + else + if t₂ + service_time[i₂] + duration[i₂, i₁] <= t₁ + add_edge!(graph, i₂, i₁) + end + end + end + end + + return graph +end + +""" +$TYPEDSIGNATURES + +Create the acyclic digraph associated with the given VSP `state`. +""" +function create_graph(state::DVSPState) + return create_graph(state.state_instance) +end + +""" +$TYPEDSIGNATURES + +Retrieve routes solution from the given MIP solution `y` matrix and `graph`. +""" +function retrieve_routes(y::AbstractArray, graph::AbstractGraph) + nb_tasks = nv(graph) + job_indices = 2:(nb_tasks) + routes = Vector{Int}[] + + start = [i for i in job_indices if y[1, i] ≈ 1] + for task in start + route = Int[] + current_task = task + while current_task != 1 # < nb_tasks + push!(route, current_task) + local next_task + for i in outneighbors(graph, current_task) + if isapprox(y[current_task, i], 1; atol=0.1) + next_task = i + break + end + end + current_task = next_task + end + push!(routes, route) + end + return routes +end + +""" +$TYPEDSIGNATURES + +Solve the Prize Collecting Vehicle Scheduling Problem defined by `instance` and prize vector `θ`. +""" +function prize_collecting_vsp( + θ::AbstractVector; instance::DVSPState, model_builder=highs_model, kwargs... +) + (; duration) = instance.state_instance + graph = create_graph(instance) + + model = model_builder() + set_silent(model) + + nb_nodes = nv(graph) + job_indices = 2:(nb_nodes) + + @variable(model, y[i=1:nb_nodes, j=1:nb_nodes; has_edge(graph, i, j)] >= 0) + + θ_ext = fill(0.0, location_count(instance)) # no prize for must dispatch requests, only hard constraints + θ_ext[instance.is_postponable] .= θ + + @objective( + model, + Max, + sum( + (θ_ext[dst(edge)] - duration[src(edge), dst(edge)]) * y[src(edge), dst(edge)] + for edge in edges(graph) + ) + ) + @constraint( + model, + flow[i in 2:nb_nodes], + sum(y[j, i] for j in inneighbors(graph, i)) == + sum(y[i, j] for j in outneighbors(graph, i)) + ) + @constraint( + model, demand[i in job_indices], sum(y[j, i] for j in inneighbors(graph, i)) <= 1 + ) + # must dispatch constraints + @constraint( + model, + demand_must_dispatch[i in job_indices; instance.is_must_dispatch[i]], + sum(y[j, i] for j in inneighbors(graph, i)) == 1 + ) + + optimize!(model) + + return retrieve_routes(value.(y), graph) +end + +function oracle(θ; instance::DVSPState, kwargs...) + routes = prize_collecting_vsp(θ; instance=instance, kwargs...) + return VSPSolution( + routes; max_index=location_count(instance.state_instance) + ).edge_matrix +end + +function g(y; instance, kwargs...) + return vec(sum(y[:, instance.is_postponable]; dims=1)) +end + +function h(y, duration) + value = 0.0 + N = size(duration, 1) + for i in 1:N + for j in 1:N + value -= y[i, j] * duration[i, j] + end + end + return value +end + +function h(y; instance, kwargs...) + return h(y, instance.state_instance.duration) +end diff --git a/src/DynamicVehicleScheduling/plot.jl b/src/DynamicVehicleScheduling/plot.jl new file mode 100644 index 0000000..adb0fa6 --- /dev/null +++ b/src/DynamicVehicleScheduling/plot.jl @@ -0,0 +1,138 @@ +function plot_instance(env::DVSPEnv; kwargs...) + return plot_instance(env.instance.static_instance; kwargs...) +end + +# """ +# $TYPEDSIGNATURES + +# Plot the environment of a DVSPEnv, restricted to the given `epoch_indices` (all epoch if not given). +# """ +# function plot_environment( +# env::DVSPEnv; +# customer_markersize=4, +# depot_markersize=7, +# alpha_depot=0.8, +# depot_color=:lightgreen, +# epoch_indices=nothing, +# kwargs..., +# ) +# draw_all_epochs!(env) + +# epoch_appearance = env.request_epoch +# coordinates = coordinate(get_state(env)) + +# epoch_indices = isnothing(epoch_indices) ? get_epoch_indices(env) : epoch_indices + +# xlims = (minimum(c.x for c in coordinates), maximum(c.x for c in coordinates)) +# ylims = (minimum(c.y for c in coordinates), maximum(c.y for c in coordinates)) + +# fig = plot(; +# legend=:topleft, +# xlabel="x coordinate", +# ylabel="y coordinate", +# xlims, +# ylims, +# kwargs..., +# ) + +# for epoch in epoch_indices +# requests = findall(epoch_appearance .== epoch) +# x = [coordinates[request].x for request in requests] +# y = [coordinates[request].y for request in requests] +# scatter!( +# fig, x, y; label="Epoch $epoch", marker=:circle, markersize=customer_markersize +# ) +# end +# scatter!( +# fig, +# [coordinates[1].x], +# [coordinates[1].y]; +# label="Depot", +# markercolor=depot_color, +# marker=:rect, +# markersize=depot_markersize, +# alpha=alpha_depot, +# ) + +# return fig +# end + +# """ +# $TYPEDSIGNATURES + +# Plot the given `routes`` for a VSP `state`. +# """ +# function plot_epoch(state::DVSPState, routes; kwargs...) +# (; coordinate, start_time) = state.instance +# x_depot = coordinate[1].x +# y_depot = coordinate[1].y +# X = [p.x for p in coordinate] +# Y = [p.y for p in coordinate] +# markersize = 5 +# fig = plot(; +# legend=:topleft, xlabel="x", ylabel="y", clim=(0.0, maximum(start_time)), kwargs... +# ) +# for route in routes +# x_points = vcat(x_depot, X[route], x_depot) +# y_points = vcat(y_depot, Y[route], y_depot) +# plot!(fig, x_points, y_points; label=nothing) +# end +# scatter!( +# fig, +# [x_depot], +# [y_depot]; +# label="depot", +# markercolor=:lightgreen, +# markersize, +# marker=:rect, +# ) +# if sum(state.is_postponable) > 0 +# scatter!( +# fig, +# X[state.is_postponable], +# Y[state.is_postponable]; +# label="Postponable customers", +# marker_z=start_time[state.is_postponable], +# markersize, +# colormap=:turbo, +# marker=:utriangle, +# ) +# end +# if sum(state.is_must_dispatch) > 0 +# scatter!( +# fig, +# X[state.is_must_dispatch], +# Y[state.is_must_dispatch]; +# label="Must-dispatch customers", +# marker_z=start_time[state.is_must_dispatch], +# markersize, +# colormap=:turbo, +# marker=:star5, +# ) +# end +# return fig +# end + +# """ +# $TYPEDSIGNATURES + +# Create a plot of routes for each epoch. +# """ +# function plot_routes(env::DVSPEnv, routes; epoch_indices=nothing, kwargs...) +# reset!(env) +# epoch_indices = isnothing(epoch_indices) ? get_epoch_indices(env) : epoch_indices + +# coordinates = env.config.static_instance.coordinate +# xlims = (minimum(c.x for c in coordinates), maximum(c.x for c in coordinates)) +# ylims = (minimum(c.y for c in coordinates), maximum(c.y for c in coordinates)) + +# figs = map(epoch_indices) do epoch +# s = next_epoch!(env) +# fig = plot_epoch( +# s, state_route_from_env_routes(env, routes[epoch]); xlims, ylims, kwargs... +# ) +# apply_decision!(env, routes[epoch]) +# return fig +# end +# return figs +# end diff --git a/src/DynamicVehicleScheduling/policy.jl b/src/DynamicVehicleScheduling/policy.jl new file mode 100644 index 0000000..244dc66 --- /dev/null +++ b/src/DynamicVehicleScheduling/policy.jl @@ -0,0 +1,36 @@ +function greedy_policy(env::DVSPEnv; model_builder=highs_model) + _, state = observe(env) + (; is_postponable) = state + nb_postponable_requests = sum(is_postponable) + θ = ones(nb_postponable_requests) * 1e9 + routes = prize_collecting_vsp(θ; instance=state, model_builder) + @assert is_feasible(state, routes) + return routes +end + +function lazy_policy(env::DVSPEnv; model_builder=highs_model) + _, state = observe(env) + nb_postponable_requests = sum(state.is_postponable) + θ = ones(nb_postponable_requests) * -1e9 + routes = prize_collecting_vsp(θ; instance=state, model_builder) + @assert is_feasible(state, routes) + return routes +end + +""" +$TYPEDEF + +Kleopatra policy for the Dynamic Vehicle Scheduling Problem. +""" +struct KleopatraVSPPolicy{P} + prize_predictor::P +end + +function (π::KleopatraVSPPolicy)(env::DVSPEnv; model_builder=highs_model) + x, state = observe(env) + (; prize_predictor) = π + # x = has_2D_features ? compute_2D_features(env) : compute_features(env) + θ = prize_predictor(x) + routes = prize_collecting_vsp(θ; instance=state, model_builder) + return routes +end diff --git a/src/DynamicVehicleScheduling/scenario.jl b/src/DynamicVehicleScheduling/scenario.jl new file mode 100644 index 0000000..9059477 --- /dev/null +++ b/src/DynamicVehicleScheduling/scenario.jl @@ -0,0 +1,51 @@ + +struct Scenario + "indices of the new requests in each epoch" + indices::Vector{Vector{Int}} + "service times of the new requests in each epoch" + service_time::Vector{Vector{Float64}} + "start times of the new requests in each epoch" + start_time::Vector{Vector{Float64}} +end + +function Base.getindex(scenario::Scenario, idx::Integer) + return (; + indices=scenario.indices[idx], + service_time=scenario.service_time[idx], + start_time=scenario.start_time[idx], + ) +end + +function Utils.generate_scenario( + instance::Instance; seed=nothing, rng::AbstractRNG=MersenneTwister(seed) +) + (; Δ_dispatch, static_instance, last_epoch, epoch_duration, max_requests_per_epoch) = + instance + (; duration, start_time, service_time) = static_instance + N = customer_count(static_instance) + depot = 1 + + new_indices = Vector{Int}[] + new_service_time = Vector{Float64}[] + new_start_time = Vector{Float64}[] + + for epoch in 1:last_epoch + time = epoch_duration * (epoch - 1) + Δ_dispatch + + coordinate_indices = sample_indices(rng, max_requests_per_epoch, N) + start_time_indices = sample_indices(rng, max_requests_per_epoch, N) + service_time_indices = sample_indices(rng, max_requests_per_epoch, N) + + is_feasible = + time .+ duration[depot, coordinate_indices] .<= start_time[start_time_indices] + + push!(new_indices, coordinate_indices[is_feasible]) + push!(new_service_time, service_time[service_time_indices[is_feasible]]) + push!(new_start_time, start_time[start_time_indices[is_feasible]]) + end + return Scenario(new_indices, new_service_time, new_start_time) +end + +function Utils.generate_scenario(sample::DataSample; kwargs...) + return Utils.generate_scenario(sample.instance; kwargs...) +end diff --git a/src/DynamicVehicleScheduling/state.jl b/src/DynamicVehicleScheduling/state.jl new file mode 100644 index 0000000..0d0a177 --- /dev/null +++ b/src/DynamicVehicleScheduling/state.jl @@ -0,0 +1,225 @@ +""" +$TYPEDSIGNATURES + +State data structure for the Dynamic Vehicle Scheduling Problem. +""" +@kwdef mutable struct DVSPState{I} + "current epoch number" + current_epoch::Int = 1 + "list of location indices from the upper instance (useful for adding new customers)" + location_indices::Vector{Int} = Int[] + "associated (static) vehicle scheduling instance" + state_instance::I = StaticInstance() + "for each location, 1 if the request must be dispatched, 0 otherwise. The depot is always 0." + is_must_dispatch::BitVector = falses(0) + "for each location, 1 if the request can be postponed, 0 otherwise. The depot is always 0." + is_postponable::BitVector = falses(0) +end + +function Base.show(io::IO, state::DVSPState) + return print( + io, + "DVSPState(", + "current_epoch=", + state.current_epoch, + ", ", + "location_indices=", + state.location_indices, + ", ", + "is_must_dispatch=", + state.is_must_dispatch, + ", ", + "is_postponable=", + state.is_postponable, + ")", + ) +end + +function reset_state!( + state::DVSPState, instance::Instance; indices, service_time, start_time +) + (; epoch_duration, Δ_dispatch, static_instance) = instance + indices_with_depot = vcat(1, indices) + service_time_with_depot = vcat(0.0, service_time) + start_time_with_depot = vcat(0.0, start_time) + + coordinates = coordinate(static_instance)[indices_with_depot] + duration_matrix = duration(static_instance)[indices_with_depot, indices_with_depot] + + is_must_dispatch = falses(length(indices_with_depot)) + is_must_dispatch[2:end] .= + Δ_dispatch .+ epoch_duration .+ @view(duration_matrix[1, 2:end]) .> start_time + + is_postponable = falses(length(is_must_dispatch)) + is_postponable[2:end] .= .!is_must_dispatch[2:end] + + state.current_epoch = 1 + state.state_instance = StaticInstance(; + service_time=service_time_with_depot, + start_time=start_time_with_depot, + coordinate=coordinates, + duration=duration_matrix, + ) + state.is_must_dispatch = is_must_dispatch + state.is_postponable = is_postponable + state.location_indices = indices_with_depot + return nothing +end + +function DVSPState(instance::Instance; indices, service_time, start_time) + state = DVSPState() + reset_state!(state, instance; indices=indices, service_time=service_time, start_time) + return state +end + +current_epoch(state::DVSPState) = state.current_epoch + +""" +$TYPEDSIGNATURES + +Return the number of locations in `state` (customers + depot). +""" +location_count(state::DVSPState) = location_count(state.state_instance) + +""" +$TYPEDSIGNATURES + +Return the number of customers in `state`. +""" +customer_count(state::DVSPState) = customer_count(state.state_instance) + +""" +$TYPEDSIGNATURES + +Get the service time vector +""" +service_time(state::DVSPState) = service_time(state.state_instance) + +""" +$TYPEDSIGNATURES + +Get the coordinates vector. +""" +coordinate(state::DVSPState) = coordinate(state.state_instance) + +""" +$TYPEDSIGNATURES + +Get the duration matrix. +""" +duration(state::DVSPState) = duration(state.state_instance) + +""" +$TYPEDSIGNATURES + +Get the start time vector. +""" +start_time(state::DVSPState) = start_time(state.state_instance) + +""" +$TYPEDSIGNATURES + +Check if the given routes are feasible. +Routes should be given with global indexation. +Use `env_routes_from_state_routes` if needed to convert the indices beforehand. +""" +function is_feasible(state::DVSPState, routes::Vector{Vector{Int}}; verbose::Bool=false) + (; is_must_dispatch, state_instance) = state + (; duration, start_time, service_time) = state_instance + is_dispatched = falses(length(is_must_dispatch)) + + # Check that routes follow time constraints + for route in routes + is_dispatched[route] .= true + current = 1 # start at the depot + current_time = start_time[current] + for next in route + current_time += duration[current, next] + if current_time > start_time[next] + verbose && + @warn "Route $route is infeasible: time constraint violated at location $next" + return false + end + current_time += service_time[next] + current = next + end + end + + # Check that all must dispatch requests are dispatched + if all(is_dispatched[is_must_dispatch]) + return true + else + verbose && @warn "Not all must-dispatch requests are dispatched" + return false + end +end + +""" +remove dispatched customers, and update must-dispatch and postponable flags. +""" +function apply_routes!( + state::DVSPState, routes::Vector{Vector{Int}}; check_feasibility::Bool=true +) + check_feasibility && @assert is_feasible(state, routes; verbose=true) + (; is_must_dispatch, is_postponable, state_instance, location_indices) = state + c = cost(state, routes) + + # Remove dispatched customers + N = location_count(state_instance) + undispatched_indices = trues(N) + undispatched_indices[vcat(routes...)] .= false + state.state_instance = StaticInstance(; + coordinate=state_instance.coordinate[undispatched_indices], + service_time=state_instance.service_time[undispatched_indices], + start_time=state_instance.start_time[undispatched_indices], + duration=state_instance.duration[undispatched_indices, undispatched_indices], + ) + state.is_must_dispatch = is_must_dispatch[undispatched_indices] + state.is_postponable = is_postponable[undispatched_indices] + state.location_indices = location_indices[undispatched_indices] + return c +end + +function cost(state::DVSPState, routes::Vector{Vector{Int}}) + return cost(routes, duration(state.state_instance)) +end + +function add_new_customers!( + state::DVSPState, instance::Instance; indices, service_time, start_time +) + (; state_instance, is_must_dispatch, is_postponable, location_indices) = state + + updated_indices = vcat(location_indices, indices) + updated_service_time = vcat(state_instance.service_time, service_time) + updated_start_time = vcat(state_instance.start_time, start_time) + updated_coordinates = instance.static_instance.coordinate[updated_indices] + updated_duration = instance.static_instance.duration[updated_indices, updated_indices] + is_must_dispatch = falses(length(updated_indices)) + is_postponable = falses(length(updated_indices)) + + state.state_instance = StaticInstance(; + coordinate=updated_coordinates, + service_time=updated_service_time, + start_time=updated_start_time, + duration=updated_duration, + ) + + # Compute must-dispatch flags + epoch_duration = instance.epoch_duration + Δ_dispatch = instance.Δ_dispatch + planning_start_time = (state.current_epoch - 1) * epoch_duration + Δ_dispatch + if state.current_epoch == last_epoch(instance) + # If we are in the last epoch, all requests must be dispatched + is_must_dispatch[2:end] .= true + else + is_must_dispatch[2:end] .= + planning_start_time .+ epoch_duration .+ @view(updated_duration[1, 2:end]) .> + updated_start_time[2:end] + end + is_postponable[2:end] .= .!is_must_dispatch[2:end] + + state.is_must_dispatch = is_must_dispatch + state.is_postponable = is_postponable + state.location_indices = updated_indices + return nothing +end diff --git a/src/DynamicVehicleScheduling/static_vsp/instance.jl b/src/DynamicVehicleScheduling/static_vsp/instance.jl new file mode 100644 index 0000000..97091a0 --- /dev/null +++ b/src/DynamicVehicleScheduling/static_vsp/instance.jl @@ -0,0 +1,65 @@ +""" +$TYPEDEF + +Instance data structure for the (deterministic and static) Vehicle Scheduling Problem. + +# Fields +$TYPEDFIELDS +""" +@kwdef struct StaticInstance{T} + "coordinates of the locations. The first one is always the depot." + coordinate::Vector{Point{T}} = Point{Float64}[] + "service time at each location" + service_time::Vector{T} = Float64[] + "start time at each location" + start_time::Vector{T} = Float64[] + "duration matrix between locations" + duration::Matrix{T} = zeros(Float64, 0, 0) +end + +function Base.show(io::IO, instance::StaticInstance) + N = customer_count(instance) + return print(io, "VSPInstance with $N customers") +end + +""" +$TYPEDSIGNATURES + +Return the number of locations in `instance` (customers + depot). +""" +location_count(instance::StaticInstance) = length(instance.coordinate) + +""" +$TYPEDSIGNATURES + +Return the number of customers in `instance` (excluding the depot). +""" +customer_count(instance::StaticInstance) = location_count(instance) - 1 + +""" +$TYPEDSIGNATURES + +Get the service time vector. +""" +service_time(instance::StaticInstance) = instance.service_time + +""" +$TYPEDSIGNATURES + +Get the coordinates vector. +""" +coordinate(instance::StaticInstance) = instance.coordinate + +""" +$TYPEDSIGNATURES + +Get the duration matrix. +""" +duration(instance::StaticInstance) = instance.duration + +""" +$TYPEDSIGNATURES + +Get the start time vector. +""" +start_time(instance::StaticInstance) = instance.start_time diff --git a/src/DynamicVehicleScheduling/static_vsp/parsing.jl b/src/DynamicVehicleScheduling/static_vsp/parsing.jl new file mode 100644 index 0000000..7bd7f92 --- /dev/null +++ b/src/DynamicVehicleScheduling/static_vsp/parsing.jl @@ -0,0 +1,95 @@ +""" +$TYPEDSIGNATURES + +Create a `VSPInstance` from file `filepath` containing a VRPTW instance. +It uses time window values to compute task times as the middle of the interval. + +Round all values to `Int` if `rounded=true`. +Normalize all time values by the `normalization` parameter. +""" +function read_vsp_instance(filepath::String; rounded::Bool=false, normalization=3600.0) + type = rounded ? Int : Float64 + mode = "" + local edge_weight_type + local edge_weight_format + duration_matrix = Vector{type}[] + nb_locations = 0 + local demand + local service_time + local coordinates + local start_time + + file = open(filepath, "r") + for line in eachline(file) + line = strip(line, [' ', '\n', '\t']) + if line == "" + continue + elseif startswith(line, "DIMENSION") + nb_locations = parse(Int, split(line, " : ")[2]) + demand = zeros(type, nb_locations) + service_time = zeros(type, nb_locations) + coordinates = zeros(type, (nb_locations, 2)) + start_time = zeros(type, nb_locations) + elseif startswith(line, "EDGE_WEIGHT_TYPE") + edge_weight_type = split(line, " : ")[2] + elseif startswith(line, "EDGE_WEIGHT_FORMAT") + edge_weight_format = split(line, " : ")[2] + elseif startswith(line, "NODE_COORD_SECTION") + mode = "coord" + elseif line == "DEMAND_SECTION" + mode = "demand" + elseif line == "DEPOT_SECTION" + mode = "depot" + elseif line == "EDGE_WEIGHT_SECTION" + mode = "edge_weights" + @assert edge_weight_type == "EXPLICIT" + @assert edge_weight_format == "FULL_MATRIX" + elseif line == "TIME_WINDOW_SECTION" + mode = "time_windows" + elseif line == "SERVICE_TIME_SECTION" + mode = "service_t" + elseif line == "EOF" + break + elseif mode == "coord" + node, x, y = split(line) # Split by whitespace or \t, skip duplicate whitespace + node = parse(Int, node) + x, y = (parse(type, x), parse(type, y)) + coordinates[node, :] = [x, y] + elseif mode == "demand" + node, d = split(line) + node, d = parse(Int, node), parse(type, d) + if node == 1 # depot + @assert d == 0 + end + demand[node] = d + elseif mode == "edge_weights" + push!(duration_matrix, [parse(type, e) for e in split(line)]) + elseif mode == "service_t" + node, t = split(line) + node = parse(Int, node) + t = parse(type, t) + if node == 1 # depot + @assert t == 0 + end + service_time[node] = t + elseif mode == "time_windows" + node, l, u = split(line) + node = parse(Int, node) + l, u = parse(type, l), parse(type, u) + start_time[node] = (u + l) / 2 + end + end + close(file) + + duration = mapreduce(permutedims, vcat, duration_matrix) + + coordinate = [ + Point(x / normalization, y / normalization) for + (x, y) in zip(coordinates[:, 1], coordinates[:, 2]) + ] + service_time ./= normalization + start_time ./= normalization + duration ./= normalization + + return StaticInstance(; coordinate, service_time, start_time, duration) +end diff --git a/src/DynamicVehicleScheduling/static_vsp/plot.jl b/src/DynamicVehicleScheduling/static_vsp/plot.jl new file mode 100644 index 0000000..515ab3d --- /dev/null +++ b/src/DynamicVehicleScheduling/static_vsp/plot.jl @@ -0,0 +1,39 @@ +""" +$TYPEDSIGNATURES + +Plot the given static VSP `instance`. +""" +function plot_instance( + instance::StaticInstance; + customer_markersize=4, + depot_markersize=7, + alpha_depot=0.8, + customer_color=:lightblue, + depot_color=:lightgreen, + kwargs..., +) + x = [p.x for p in instance.coordinate] + y = [p.y for p in instance.coordinate] + + fig = plot(; legend=:topleft, xlabel="x coordinate", ylabel="y coordinate", kwargs...) + scatter!( + fig, + x[2:end], + y[2:end]; + label="Customers", + markercolor=customer_color, + marker=:circle, + markersize=customer_markersize, + ) + scatter!( + fig, + [x[1]], + [y[1]]; + label="Depot", + markercolor=depot_color, + marker=:rect, + markersize=depot_markersize, + alpha=alpha_depot, + ) + return fig +end diff --git a/src/DynamicVehicleScheduling/static_vsp/solution.jl b/src/DynamicVehicleScheduling/static_vsp/solution.jl new file mode 100644 index 0000000..d6fb25e --- /dev/null +++ b/src/DynamicVehicleScheduling/static_vsp/solution.jl @@ -0,0 +1,50 @@ +""" +$TYPEDEF + +Solution for the static Vehicle Scheduling Problem. + +# Fields +$TYPEDFIELDS +""" +struct VSPSolution + "list of routes, each route being a list of request indices in corresponding instance (excluding the depot)." + routes::Vector{Vector{Int}} + "size (nb_locations, nb_locations). `edge_matrix[i, j]` is equal to 1 if a route takes edge `(i, j)`." + edge_matrix::BitMatrix +end + +""" +$TYPEDSIGNATURES + +Get routes from `solution`. +""" +routes(solution::VSPSolution) = solution.routes + +""" +$TYPEDSIGNATURES + +Get edge matrix from `solution`. +""" +edge_matrix(solution::VSPSolution) = solution.edge_matrix + +""" +$TYPEDSIGNATURES + +Build a `VSPSolution` from routes. Set `max_index` to manually define the size of the `edge_index` matrix. +""" +function VSPSolution(routes::Vector{Vector{Int}}; max_index=nothing) + if length(routes) == 0 && isnothing(max_index) + return VSPSolution(routes, falses(0, 0)) + end + N = isnothing(max_index) ? maximum(maximum(route) for route in routes) : max_index + edge_matrix = falses(N, N) + for route in routes + old = 1 + for r in route + edge_matrix[old, r] = true + old = r + end + edge_matrix[old, 1] = true + end + return VSPSolution(routes, edge_matrix) +end diff --git a/src/DynamicVehicleScheduling/utils.jl b/src/DynamicVehicleScheduling/utils.jl new file mode 100644 index 0000000..bd1dfe8 --- /dev/null +++ b/src/DynamicVehicleScheduling/utils.jl @@ -0,0 +1,37 @@ +""" +$TYPEDSIGNATURES + +Sample k random different indices from 2 to N+1. +""" +sample_indices(rng::AbstractRNG, k, N) = randperm(rng, N)[1:k] .+ 1 + +""" +$TYPEDSIGNATURES + +Compute the total cost of a set of routes given a distance matrix, i.e. the sum of the distances between each location in the route. +Note that the first location is implicitly assumed to be the depot, and should not appear in the route. +""" +function cost(routes::Vector{Vector{Int}}, duration::AbstractMatrix) + total = zero(eltype(duration)) + for route in routes + current_location = 1 + for r in route + total += duration[current_location, r] + current_location = r + end + total += duration[current_location, 1] + end + return total +end + +""" +$TYPEDEF + +Basic point structure. +""" +struct Point{T} + x::T + y::T +end + +Base.show(io::IO, p::Point) = print(io, "($(p.x), $(p.y))") diff --git a/src/FixedSizeShortestPath/FixedSizeShortestPath.jl b/src/FixedSizeShortestPath/FixedSizeShortestPath.jl index fd60de2..46a22fe 100644 --- a/src/FixedSizeShortestPath/FixedSizeShortestPath.jl +++ b/src/FixedSizeShortestPath/FixedSizeShortestPath.jl @@ -106,42 +106,26 @@ end """ $TYPEDSIGNATURES -Generate dataset for the shortest path problem. +Generate a labeled sample for the fixed size shortest path benchmark. """ -function Utils.generate_dataset( - bench::FixedSizeShortestPathBenchmark, - dataset_size::Int=10; - seed::Int=0, - type::Type=Float32, +function Utils.generate_sample( + bench::FixedSizeShortestPathBenchmark, rng::AbstractRNG; type::Type=Float32 ) - # Set seed - rng = MersenneTwister(seed) (; graph, p, deg, ν) = bench - + features = randn(rng, Float32, bench.p) E = Graphs.ne(graph) - - # Features - features = [randn(rng, type, p) for _ in 1:dataset_size] - # True weights B = rand(rng, Bernoulli(0.5), E, p) ξ = if ν == 0.0 - [ones(type, E) for _ in 1:dataset_size] + ones(type, E) else - [rand(rng, Uniform{type}(1 - ν, 1 + ν), E) for _ in 1:dataset_size] + rand(rng, Uniform{type}(1 - ν, 1 + ν), E) end - costs = [ - -(1 .+ (3 .+ B * zᵢ ./ type(sqrt(p))) .^ deg) .* ξᵢ for (ξᵢ, zᵢ) in zip(ξ, features) - ] - - shortest_path_maximizer = Utils.generate_maximizer(bench) - - # Label solutions - solutions = shortest_path_maximizer.(costs) - return [ - DataSample(; x, θ_true, y_true) for - (x, θ_true, y_true) in zip(features, costs, solutions) - ] + costs = -(1 .+ (3 .+ B * features ./ type(sqrt(p))) .^ deg) .* ξ + + maximizer = Utils.generate_maximizer(bench) + solution = maximizer(costs) + return DataSample(; x=features, θ_true=costs, y_true=solution) end """ diff --git a/src/PortfolioOptimization/PortfolioOptimization.jl b/src/PortfolioOptimization/PortfolioOptimization.jl index 308770a..5fc7e8f 100644 --- a/src/PortfolioOptimization/PortfolioOptimization.jl +++ b/src/PortfolioOptimization/PortfolioOptimization.jl @@ -7,7 +7,7 @@ using Flux: Chain, Dense using Ipopt: Ipopt using JuMP: @variable, @objective, @constraint, optimize!, value, Model, set_silent using LinearAlgebra: I -using Random: MersenneTwister +using Random: AbstractRNG, MersenneTwister """ $TYPEDEF @@ -85,35 +85,21 @@ end """ $TYPEDSIGNATURES -Generate a dataset of labeled instances for the portfolio optimization problem. +Generate a labeled sample for the portfolio optimization problem. """ -function Utils.generate_dataset( - bench::PortfolioOptimizationBenchmark, - dataset_size::Int=10; - seed::Int=0, - type::Type=Float32, +function Utils.generate_sample( + bench::PortfolioOptimizationBenchmark, rng::AbstractRNG; type::Type=Float32 ) (; d, p, deg, ν, L, f) = bench - rng = MersenneTwister(seed) - - # Features - features = [randn(rng, type, p) for _ in 1:dataset_size] - - # True weights + features = randn(rng, type, p) B = rand(rng, Bernoulli(0.5), d, p) - c̄ = [ - (0.05 / type(sqrt(p)) .* B * features[i] .+ 0.1^(1 / deg)) .^ deg for - i in 1:dataset_size - ] - costs = [c̄ᵢ .+ L * f .+ 0.01 .* ν .* randn(rng, type, d) for c̄ᵢ in c̄] + c̄ = (0.05 / type(sqrt(p)) .* B * features .+ 0.1^(1 / deg)) .^ deg + costs = c̄ .+ L * f .+ 0.01 * ν * randn(rng, type, d) maximizer = Utils.generate_maximizer(bench) - solutions = maximizer.(costs) + solution = maximizer(costs) - return [ - DataSample(; x, θ_true, y_true) for - (x, θ_true, y_true) in zip(features, costs, solutions) - ] + return DataSample(; x=features, θ_true=costs, y_true=solution) end """ diff --git a/src/Ranking/Ranking.jl b/src/Ranking/Ranking.jl index 8b93b8a..c6ec398 100644 --- a/src/Ranking/Ranking.jl +++ b/src/Ranking/Ranking.jl @@ -61,22 +61,16 @@ end """ $TYPEDSIGNATURES -Generate a dataset of labeled instances for the ranking problem. +Generate a labeled sample for the ranking problem. """ -function Utils.generate_dataset( - bench::RankingBenchmark, dataset_size::Int=10; seed::Int=0, noise_std=0.0 +function Utils.generate_sample( + bench::RankingBenchmark, rng::AbstractRNG; noise_std::Float32=0.0f0 ) (; instance_dim, nb_features, encoder) = bench - rng = MersenneTwister(seed) - features = [randn(rng, Float32, nb_features, instance_dim) for _ in 1:dataset_size] - costs = encoder.(features) - noisy_solutions = [ - ranking(θ .+ noise_std * randn(rng, Float32, instance_dim)) for θ in costs - ] - return [ - DataSample(; x, θ_true, y_true) for - (x, θ_true, y_true) in zip(features, costs, noisy_solutions) - ] + features = randn(rng, Float32, nb_features, instance_dim) + costs = encoder(features) + noisy_solution = ranking(costs .+ noise_std * randn(rng, Float32, instance_dim)) + return DataSample(; x=features, θ_true=costs, y_true=noisy_solution) end """ diff --git a/src/StochasticVehicleScheduling/StochasticVehicleScheduling.jl b/src/StochasticVehicleScheduling/StochasticVehicleScheduling.jl index 150f147..41801c5 100644 --- a/src/StochasticVehicleScheduling/StochasticVehicleScheduling.jl +++ b/src/StochasticVehicleScheduling/StochasticVehicleScheduling.jl @@ -73,43 +73,32 @@ end """ $TYPEDSIGNATURES -Create a dataset of `dataset_size` instances for the given `StochasticVehicleSchedulingBenchmark`. -If you want to not add label solutions in the dataset, set `compute_solutions=false`. +Generate a sample for the given `StochasticVehicleSchedulingBenchmark`. +If you want to not add label solutions in the sample, set `compute_solutions=false`. By default, they will be computed using column generation. Note that computing solutions can be time-consuming, especially for large instances. You can also use instead `compact_mip` or `compact_linearized_mip` as the algorithm to compute solutions. If you want to provide a custom algorithm to compute solutions, you can pass it as the `algorithm` keyword argument. If `algorithm` takes keyword arguments, you can pass them as well directly in `kwargs...`. -If `store_city=false`, the coordinates and unnecessary information about instances will not be stored in the dataset. +If `store_city=false`, the coordinates and unnecessary information about instances will not be stored in the sample. """ -function Utils.generate_dataset( +function Utils.generate_sample( benchmark::StochasticVehicleSchedulingBenchmark, - dataset_size::Int; + rng::AbstractRNG; + store_city=true, compute_solutions=true, - seed=nothing, - rng=MersenneTwister(0), algorithm=column_generation_algorithm, - store_city=true, kwargs..., ) (; nb_tasks, nb_scenarios) = benchmark - Random.seed!(rng, seed) - instances = [ - Instance(; nb_tasks, nb_scenarios, rng, store_city) for _ in 1:dataset_size - ] - features = get_features.(instances) - if compute_solutions - solutions = [algorithm(instance; kwargs...).value for instance in instances] - return [ - DataSample(; x=feature, instance, y_true=solution) for - (instance, feature, solution) in zip(instances, features, solutions) - ] + instance = Instance(; nb_tasks, nb_scenarios, rng, store_city) + x = get_features(instance) + y_true = if compute_solutions + algorithm(instance; kwargs...) + else + nothing end - # else - return [ - DataSample(; x=feature, instance) for - (instance, feature) in zip(instances, features) - ] + return DataSample(; x, instance, y_true) end """ @@ -126,7 +115,7 @@ end $TYPEDSIGNATURES """ function Utils.generate_maximizer( - bench::StochasticVehicleSchedulingBenchmark; model_builder=highs_model + ::StochasticVehicleSchedulingBenchmark; model_builder=highs_model ) return StochasticVechicleSchedulingMaximizer(model_builder) end @@ -156,10 +145,7 @@ end $TYPEDSIGNATURES """ function plot_instance( - ::StochasticVehicleSchedulingBenchmark, - sample::DataSample{<:Instance{City}}; - color_scheme=:lightrainbow, - kwargs..., + ::StochasticVehicleSchedulingBenchmark, sample::DataSample{<:Instance{City}}; kwargs... ) (; tasks, district_width, width) = sample.instance.city ticks = 0:district_width:width @@ -208,7 +194,6 @@ function plot_instance( marker_z=task.end_time, colormap=:turbo, label=nothing, - # color=palette[max(floor(Int, task.end_time), 1)], ) annotate!(fig, (points[1]..., text("$(i_task)", 10))) end diff --git a/src/StochasticVehicleScheduling/solution/algorithms/column_generation.jl b/src/StochasticVehicleScheduling/solution/algorithms/column_generation.jl index dbd2fd6..1bfbe1f 100644 --- a/src/StochasticVehicleScheduling/solution/algorithms/column_generation.jl +++ b/src/StochasticVehicleScheduling/solution/algorithms/column_generation.jl @@ -189,5 +189,5 @@ function column_generation_algorithm( end col_solution = solution_from_paths(sol, instance) - return col_solution + return col_solution.value end diff --git a/src/StochasticVehicleScheduling/solution/algorithms/deterministic_mip.jl b/src/StochasticVehicleScheduling/solution/algorithms/deterministic_mip.jl index 5f68190..9e14861 100644 --- a/src/StochasticVehicleScheduling/solution/algorithms/deterministic_mip.jl +++ b/src/StochasticVehicleScheduling/solution/algorithms/deterministic_mip.jl @@ -41,5 +41,5 @@ function deterministic_mip(instance::Instance; model_builder=highs_model, silent solution = value.(y) sol = solution_from_JuMP_array(solution, graph) - return sol + return sol.value end diff --git a/src/StochasticVehicleScheduling/solution/algorithms/local_search.jl b/src/StochasticVehicleScheduling/solution/algorithms/local_search.jl index b4f0f0f..49ae00c 100644 --- a/src/StochasticVehicleScheduling/solution/algorithms/local_search.jl +++ b/src/StochasticVehicleScheduling/solution/algorithms/local_search.jl @@ -153,5 +153,5 @@ Very simple heuristic, using [`local_search`](@ref) function local_search(instance::Instance; num_iterations=1000) _, initial_solution = solve_deterministic_VSP(instance) sol, _, _, _ = _local_search(initial_solution, instance; nb_it=num_iterations) - return sol + return sol.value end diff --git a/src/StochasticVehicleScheduling/solution/algorithms/mip.jl b/src/StochasticVehicleScheduling/solution/algorithms/mip.jl index e202569..10b0b40 100644 --- a/src/StochasticVehicleScheduling/solution/algorithms/mip.jl +++ b/src/StochasticVehicleScheduling/solution/algorithms/mip.jl @@ -78,7 +78,7 @@ function compact_linearized_mip( solution = value.(y) sol = solution_from_JuMP_array(solution, graph) - return sol + return sol.value end """ @@ -149,5 +149,5 @@ function compact_mip( solution = value.(y) sol = solution_from_JuMP_array(solution, graph) - return sol + return sol.value end diff --git a/src/SubsetSelection/SubsetSelection.jl b/src/SubsetSelection/SubsetSelection.jl index 0e738a5..085324d 100644 --- a/src/SubsetSelection/SubsetSelection.jl +++ b/src/SubsetSelection/SubsetSelection.jl @@ -17,11 +17,13 @@ without knowing their values, but only observing some features. # Fields $TYPEDFIELDS """ -struct SubsetSelectionBenchmark <: AbstractBenchmark +struct SubsetSelectionBenchmark{M} <: AbstractBenchmark "total number of items" n::Int "number of items to select" k::Int + "hidden unknown mapping from features to costs" + mapping::M end function Base.show(io::IO, bench::SubsetSelectionBenchmark) @@ -29,9 +31,14 @@ function Base.show(io::IO, bench::SubsetSelectionBenchmark) return print(io, "SubsetSelectionBenchmark(n=$n, k=$k)") end -function SubsetSelectionBenchmark(; n::Int=25, k::Int=5) +function SubsetSelectionBenchmark(; n::Int=25, k::Int=5, identity_mapping::Bool=true) @assert n >= k "number of items n must be greater than k" - return SubsetSelectionBenchmark(n, k) + mapping = if identity_mapping + copy + else + Dense(n => n; bias=false) + end + return SubsetSelectionBenchmark(n, k, mapping) end function top_k(v::AbstractVector, k::Int) @@ -54,29 +61,14 @@ end """ $TYPEDSIGNATURES -Generate a dataset of labeled instances for the subset selection problem. -The mapping between features and cost is identity. +Generate a labeled instance for the subset selection problem. """ -function Utils.generate_dataset( - bench::SubsetSelectionBenchmark, - dataset_size::Int=10; - seed::Int=0, - identity_mapping=true, -) - (; n, k) = bench - rng = MersenneTwister(seed) - features = [randn(rng, Float32, n) for _ in 1:dataset_size] - costs = if identity_mapping - copy(features) # we assume that the cost is the same as the feature - else - mapping = Dense(n => n; bias=false) - mapping.(features) - end - solutions = top_k.(costs, k) - return [ - DataSample(; x, θ_true, y_true) for - (x, θ_true, y_true) in zip(features, costs, solutions) - ] +function Utils.generate_sample(bench::SubsetSelectionBenchmark, rng::AbstractRNG) + (; n, k, mapping) = bench + features = randn(rng, Float32, n) + costs = mapping(features) + solution = top_k(costs, k) + return DataSample(; x=features, θ_true=costs, y_true=solution) end """ diff --git a/src/Utils/Utils.jl b/src/Utils/Utils.jl index 60b5b92..d738e31 100644 --- a/src/Utils/Utils.jl +++ b/src/Utils/Utils.jl @@ -5,26 +5,41 @@ using Flux: softplus using HiGHS: HiGHS using JuMP: Model using LinearAlgebra: dot +using Random: Random, MersenneTwister using SCIP: SCIP using SimpleWeightedGraphs: SimpleWeightedDiGraph using StatsBase: StatsBase using Statistics: mean include("data_sample.jl") +include("maximizers.jl") +include("environment.jl") +include("policy.jl") include("interface.jl") include("grid_graph.jl") include("misc.jl") include("model_builders.jl") -export DataSample +export DataSample, Policy +export evaluate_policy! +export TopKMaximizer -export AbstractBenchmark -export generate_dataset, - generate_statistical_model, generate_maximizer, plot_data, compute_gap +export AbstractEnvironment, get_seed, is_terminated, observe, reset!, step! + +export AbstractBenchmark, AbstractStochasticBenchmark, AbstractDynamicBenchmark +export generate_sample, generate_dataset +export generate_statistical_model, generate_maximizer +export generate_scenario +export generate_environment, generate_environments +export generate_policies +export generate_anticipative_solution + +export plot_data, compute_gap export maximizer_kwargs export grid_graph, get_path, path_to_matrix export neg_tensor, squeeze_last_dims, average_tensor export scip_model, highs_model export objective_value +export is_exogenous, is_endogenous end diff --git a/src/Utils/data_sample.jl b/src/Utils/data_sample.jl index e9a8a3c..d0cccc6 100644 --- a/src/Utils/data_sample.jl +++ b/src/Utils/data_sample.jl @@ -7,10 +7,13 @@ Data sample data structure. $TYPEDFIELDS """ @kwdef struct DataSample{ - I,F<:AbstractArray,S<:Union{AbstractArray,Nothing},C<:Union{AbstractArray,Nothing} + I, + F<:Union{AbstractArray,Nothing}, + S<:Union{AbstractArray,Nothing}, + C<:Union{AbstractArray,Nothing}, } "features" - x::F + x::F = nothing "target cost parameters (optional)" θ_true::C = nothing "target solution (optional)" @@ -19,6 +22,23 @@ $TYPEDFIELDS instance::I = nothing end +function Base.show(io::IO, d::DataSample) + fields = String[] + if !isnothing(d.x) + push!(fields, "x=$(d.x)") + end + if !isnothing(d.θ_true) + push!(fields, "θ_true=$(d.θ_true)") + end + if !isnothing(d.y_true) + push!(fields, "y_true=$(d.y_true)") + end + if !isnothing(d.instance) + push!(fields, "instance=$(d.instance)") + end + return print(io, "DataSample(", join(fields, ", "), ")") +end + """ $TYPEDSIGNATURES diff --git a/src/Utils/environment.jl b/src/Utils/environment.jl new file mode 100644 index 0000000..88eff08 --- /dev/null +++ b/src/Utils/environment.jl @@ -0,0 +1,52 @@ +""" +$TYPEDEF + +Abstract type for environments in decision-focused learning benchmarks. +""" +abstract type AbstractEnvironment end + +""" +$TYPEDSIGNATURES + +Seed accessor for environments. +By default, environments have no seed. +Override this method to provide a seed for the environment. +""" +function get_seed(::AbstractEnvironment) + return nothing +end + +""" + is_terminated(env::AbstractEnvironment) --> Bool + +Check if the environment has reached a terminal state. +""" +function is_terminated end + +""" + observe(env::AbstractEnvironment) --> Tuple + +Get the current observation from the environment. +This function should return a tuple of two elements: + 1. An array of features representing the current state of the environment. + 2. An internal state of the environment, which can be used for further processing (return `nothing` if not needed). +""" +function observe end + +""" + reset!(env::AbstractEnvironment; reset_rng::Bool, seed=get_seed(env)) --> Nothing + +Reset the environment to its initial state. +If `reset_rng` is true, the random number generator is reset to the given `seed`. +""" +function reset! end + +""" + step!(env::AbstractEnvironment, action) --> Float64 + +Perform a step in the environment with the given action. +Returns the reward received after taking the action. +This function may also update the internal state of the environment. +If the environment is terminated, it should raise an error. +""" +function step! end diff --git a/src/Utils/interface.jl b/src/Utils/interface.jl index 96a2a39..1fa1f65 100644 --- a/src/Utils/interface.jl +++ b/src/Utils/interface.jl @@ -1,10 +1,10 @@ """ $TYPEDEF -Abstract type interface for a benchmark problem. +Abstract type interface for benchmark problems. The following methods are mandatory for benchmarks: -- [`generate_dataset`](@ref) +- [`generate_dataset`](@ref) or [`generate_sample`](@ref) - [`generate_statistical_model`](@ref) - [`generate_maximizer`](@ref) @@ -15,13 +15,36 @@ The following methods are optional: """ abstract type AbstractBenchmark end +""" + generate_sample(::AbstractBenchmark, rng::AbstractRNG; kwargs...) -> DataSample + +Generate a single [`DataSample`](@ref) for given benchmark. +This is a low-level function that is used by [`generate_dataset`](@ref) to create +a dataset of samples. It is not mandatory to implement this method, but it is +recommended for benchmarks that have a well-defined way to generate individual samples. +An alternative is to directly implement [`generate_dataset`](@ref) to create a dataset +without generating individual samples. +""" +function generate_sample end + """ generate_dataset(::AbstractBenchmark, dataset_size::Int; kwargs...) -> Vector{<:DataSample} -Generate a `Vector` of [`DataSample`](@ref) of length `dataset_size` for given benchmark. +Generate a `Vector` of [`DataSample`](@ref) of length `dataset_size` for given benchmark. Content of the dataset can be visualized using [`plot_data`](@ref), when it applies. + +By default, it uses [`generate_sample`](@ref) to create each sample in the dataset, and passes any keyword arguments to it. """ -function generate_dataset end +function generate_dataset( + bench::AbstractBenchmark, + dataset_size::Int; + seed=nothing, + rng=MersenneTwister(seed), + kwargs..., +) + Random.seed!(rng, seed) + return [generate_sample(bench, rng; kwargs...) for _ in 1:dataset_size] +end """ generate_maximizer(::AbstractBenchmark; kwargs...) @@ -39,6 +62,11 @@ It's usually a Flux model, that takes a feature matrix x as input, and returns a """ function generate_statistical_model end +""" + generate_policies(::AbstractBenchmark) -> Vector{Policy} +""" +function generate_policies end + """ plot_data(::AbstractBenchmark, ::DataSample; kwargs...) @@ -154,3 +182,63 @@ function compute_gap( end, ) end + +""" +$TYPEDEF + +Abstract type interface for stochastic benchmark problems. +This type should be used for benchmarks that involve single stage stochastic optimization problems. + +It follows the same interface as [`AbstractBenchmark`](@ref), with the addition of the following methods: +- TODO +""" +abstract type AbstractStochasticBenchmark{exogenous} <: AbstractBenchmark end + +is_exogenous(::AbstractStochasticBenchmark{exogenous}) where {exogenous} = exogenous +is_endogenous(::AbstractStochasticBenchmark{exogenous}) where {exogenous} = !exogenous + +""" + generate_scenario(::AbstractStochasticBenchmark{true}, instance; kwargs...) +""" +function generate_scenario end + +""" + generate_anticipative_solution(::AbstractStochasticBenchmark{true}, instance, scenario; kwargs...) +""" +function generate_anticipative_solution end + +""" +$TYPEDEF + +Abstract type interface for dynamic benchmark problems. +This type should be used for benchmarks that involve multi-stage stochastic optimization problems. + +It follows the same interface as [`AbstractStochasticBenchmark`](@ref), with the addition of the following methods: +TODO +""" +abstract type AbstractDynamicBenchmark{exogenous} <: AbstractStochasticBenchmark{exogenous} end + +""" + generate_environment(::AbstractDynamicBenchmark, instance, rng::AbstractRNG; kwargs...) + +Initialize an environment for the given dynamic benchmark instance. +""" +function generate_environment end + +""" +$TYPEDSIGNATURES + +Generate a vector of environments for the given dynamic benchmark and dataset. +""" +function generate_environments( + bench::AbstractDynamicBenchmark, + dataset::Vector{<:DataSample}; + seed=nothing, + rng=MersenneTwister(seed), + kwargs..., +) + Random.seed!(rng, seed) + return map(dataset) do sample + generate_environment(bench, sample.instance, rng; kwargs...) + end +end diff --git a/src/Utils/maximizers.jl b/src/Utils/maximizers.jl new file mode 100644 index 0000000..ee5ceea --- /dev/null +++ b/src/Utils/maximizers.jl @@ -0,0 +1,22 @@ +""" +$TYPEDEF + +Top k maximizer. +""" +struct TopKMaximizer + k::Int +end + +""" +$TYPEDSIGNATURES + +Return the top k indices of `θ`. +""" +function (m::TopKMaximizer)(θ; kwargs...) + N = length(θ) + @assert N >= m.k "The length of θ must be at least k" + solution = partialsortperm(θ, 1:(m.k); rev=true) + res = falses(N) + res[solution] .= 1 + return res +end diff --git a/src/Utils/model_builders.jl b/src/Utils/model_builders.jl index 95df58b..4f0c838 100644 --- a/src/Utils/model_builders.jl +++ b/src/Utils/model_builders.jl @@ -5,7 +5,6 @@ Initialize a HiGHS model (with disabled logging). """ function highs_model() model = Model(HiGHS.Optimizer) - # set_attribute(model, "log_to_console", false) return model end diff --git a/src/Utils/policy.jl b/src/Utils/policy.jl new file mode 100644 index 0000000..2b3c8e5 --- /dev/null +++ b/src/Utils/policy.jl @@ -0,0 +1,102 @@ +""" +$TYPEDEF + +Policy type for decision-focused learning benchmarks. +""" +struct Policy{P} + "policy name" + name::String + "policy description" + description::String + "policy run function" + policy::P +end + +function Base.show(io::IO, p::Policy) + println(io, "$(p.name): $(p.description)") + return nothing +end +""" +$TYPEDSIGNATURES + +Run the policy and get the next decision on the given environment/instance. +""" +function (p::Policy)(args...; kwargs...) + return p.policy(args...; kwargs...) +end + +""" +$TYPEDSIGNATURES + +Run the policy on the environment and return the total reward and a dataset of observations. +By default, the environment is reset before running the policy. +""" +function evaluate_policy!(policy, env::AbstractEnvironment; kwargs...) + total_reward = 0.0 + reset!(env; reset_rng=false) + local labeled_dataset + while !is_terminated(env) + y = policy(env; kwargs...) + features, state = observe(env) + if @isdefined labeled_dataset + push!( + labeled_dataset, + DataSample(; x=features, y_true=y, instance=deepcopy(state)), + ) + else + labeled_dataset = [DataSample(; x=features, y_true=y, instance=deepcopy(state))] + end + reward = step!(env, y) + total_reward += reward + end + return total_reward, labeled_dataset +end + +# function evaluate_policy!(policy, envs::Vector{<:AbstractEnvironment}; kwargs...) +# E = length(envs) +# rewards = zeros(Float64, E) +# datasets = map(1:E) do e +# reward, dataset = evaluate_policy!(policy, envs[e]; kwargs...) +# rewards[e] = reward +# return dataset +# end +# return rewards, vcat(datasets...) +# end + +""" +$TYPEDSIGNATURES + +Evaluate the policy on the environment and return the total reward and a dataset of observations. +By default, the environment is reset before running the policy. +""" +function evaluate_policy!( + policy, env::AbstractEnvironment, episodes::Int; seed=get_seed(env), kwargs... +) + reset!(env; reset_rng=true, seed) + total_reward = 0.0 + datasets = map(1:episodes) do _i + reward, dataset = evaluate_policy!(policy, env; kwargs...) + total_reward += reward + return dataset + end + return total_reward / episodes, vcat(datasets...) +end + +""" +$TYPEDSIGNATURES + +Run the policy on the environments and return the total rewards and a dataset of observations. +By default, the environments are reset before running the policy. +""" +function evaluate_policy!( + policy, envs::Vector{<:AbstractEnvironment}, episodes::Int=1; kwargs... +) + E = length(envs) + rewards = zeros(Float64, E) + datasets = map(1:E) do e + reward, dataset = evaluate_policy!(policy, envs[e], episodes; kwargs...) + rewards[e] = reward + return dataset + end + return rewards, vcat(datasets...) +end diff --git a/src/Warcraft/Warcraft.jl b/src/Warcraft/Warcraft.jl index 669a828..c4dcbae 100644 --- a/src/Warcraft/Warcraft.jl +++ b/src/Warcraft/Warcraft.jl @@ -2,7 +2,7 @@ module Warcraft using ..Utils -using DataDeps +using DataDeps: @datadep_str using DocStringExtensions: TYPEDEF, TYPEDFIELDS, TYPEDSIGNATURES using Flux using Graphs diff --git a/test/dynamic_assortment.jl b/test/dynamic_assortment.jl new file mode 100644 index 0000000..54030fe --- /dev/null +++ b/test/dynamic_assortment.jl @@ -0,0 +1,32 @@ +@testitem "dynamic Assortment" begin + using DecisionFocusedLearningBenchmarks + using Statistics: mean + + b = DynamicAssortmentBenchmark() + + @test is_endogenous(b) + @test !is_exogenous(b) + + dataset = generate_dataset(b, 10; seed=0) + environments = generate_environments(b, dataset) + + env = environments[1] + get_seed(env) + env.seed + + policies = generate_policies(b) + expert = policies[1] + greedy = policies[2] + + r_expert, d = evaluate_policy!(expert, environments) + r_greedy, _ = evaluate_policy!(greedy, environments) + + @test mean(r_expert) >= mean(r_greedy) + + model = generate_statistical_model(b) + maximizer = generate_maximizer(b) + sample = d[1] + x = sample.x + θ = model(x) + y = maximizer(θ) +end diff --git a/test/dynamic_vsp.jl b/test/dynamic_vsp.jl new file mode 100644 index 0000000..0f890c0 --- /dev/null +++ b/test/dynamic_vsp.jl @@ -0,0 +1,49 @@ +@testitem "DVSP" begin + using DecisionFocusedLearningBenchmarks.DynamicVehicleScheduling + using Statistics: mean + + b = DynamicVehicleSchedulingBenchmark(; two_dimensional_features=true) + b2 = DynamicVehicleSchedulingBenchmark(; two_dimensional_features=false) + + @test is_exogenous(b) + @test !is_endogenous(b) + + dataset = generate_dataset(b, 10) + environments = generate_environments(b, dataset; seed=0) + + env = environments[1] + get_seed(env) + + policies = generate_policies(b) + lazy = policies[1] + greedy = policies[2] + + d = evaluate_policy!(lazy, env, 1; seed=0)[2] + + r_lazy, d = evaluate_policy!(lazy, environments, 10) + r_greedy, d = evaluate_policy!(greedy, environments, 10) + + @test mean(r_lazy) <= mean(r_greedy) + + env = environments[1] + instance = dataset[1].instance + scenario = generate_scenario(b, instance) + v, y = generate_anticipative_solution(b, env, scenario; nb_epochs=2, reset_env=true) + + maximizer = generate_maximizer(b) + + x, instance = observe(env) + model = generate_statistical_model(b) + θ = model(x) + y = maximizer(θ; instance) + + dataset2 = generate_dataset(b2, 10) + environments2 = generate_environments(b2, dataset2; seed=0) + env2 = environments2[1] + x2, instance2 = observe(env2) + model2 = generate_statistical_model(b2) + θ2 = model2(x2) + y2 = maximizer(θ2; instance=instance2) + @test size(x, 1) == 2 + @test size(x2, 1) == 14 +end diff --git a/test/subset_selection.jl b/test/subset_selection.jl index 694f7f4..d59ae54 100644 --- a/test/subset_selection.jl +++ b/test/subset_selection.jl @@ -4,14 +4,15 @@ n = 25 k = 5 - b = SubsetSelectionBenchmark(; n=n, k=k) + b_identity = SubsetSelectionBenchmark(; n=n, k=k) + b = SubsetSelectionBenchmark(; n=n, k=k, identity_mapping=false) io = IOBuffer() show(io, b) @test String(take!(io)) == "SubsetSelectionBenchmark(n=25, k=5)" - dataset = generate_dataset(b, 50) - dataset2 = generate_dataset(b, 50; identity_mapping=false) + dataset = generate_dataset(b_identity, 50) + dataset2 = generate_dataset(b, 50) model = generate_statistical_model(b) maximizer = generate_maximizer(b) diff --git a/test/utils.jl b/test/utils.jl index b071bd8..4fd4b4f 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -23,3 +23,103 @@ @test max(h, w) <= length(path) <= h + w end end + +@testitem "DataSample" begin + using DecisionFocusedLearningBenchmarks + using StableRNGs + + rng = StableRNG(1234) + + function random_sample() + return DataSample(; + x=randn(rng, 10, 5), + θ_true=rand(rng, 5), + y_true=rand(rng, 10), + instance="this is an instance", + ) + end + + sample = random_sample() + @test sample isa DataSample + + io = IOBuffer() + show(io, sample) + @test String(take!(io)) == + "DataSample(x=$(sample.x), θ_true=$(sample.θ_true), y_true=$(sample.y_true), instance=$(sample.instance))" + + # Test StatsBase methods + using StatsBase: + ZScoreTransform, + UnitRangeTransform, + fit, + transform, + transform!, + reconstruct, + reconstruct! + + # Create a dataset for testing + N = 5 + dataset = [random_sample() for _ in 1:N] + + # Test fit with ZScoreTransform + zt = fit(ZScoreTransform, dataset) + @test zt isa ZScoreTransform + + # Test fit with UnitRangeTransform + ut = fit(UnitRangeTransform, dataset) + @test ut isa UnitRangeTransform + + # Test transform (non-mutating) + dataset_zt = transform(zt, dataset) + @test length(dataset_zt) == length(dataset) + @test all(d -> d isa DataSample, dataset_zt) + + # Check that other fields are preserved + for i in 1:N + @test dataset_zt[i].θ_true == dataset[i].θ_true + @test dataset_zt[i].y_true == dataset[i].y_true + @test dataset_zt[i].instance == dataset[i].instance + end + + # Check that features are actually transformed + @test dataset_zt[1].x != dataset[1].x + + # Test transform! (mutating) + dataset_copy = deepcopy(dataset) + original_x = copy(dataset_copy[1].x) + transform!(ut, dataset_copy) + @test dataset_copy[1].x != original_x + + # Check that other fields remain unchanged after transform! + for i in 1:N + @test dataset_copy[i].θ_true == dataset[i].θ_true + @test dataset_copy[i].y_true == dataset[i].y_true + @test dataset_copy[i].instance == dataset[i].instance + end + + # Test reconstruct (non-mutating) + dataset_reconstructed = reconstruct(zt, dataset_zt) + @test length(dataset_reconstructed) == length(dataset) + + # Test round-trip consistency (should be close to original) + for i in 1:N + @test dataset_reconstructed[i].x ≈ dataset[i].x atol = 1e-10 + @test dataset_reconstructed[i].θ_true == dataset[i].θ_true + @test dataset_reconstructed[i].y_true == dataset[i].y_true + @test dataset_reconstructed[i].instance == dataset[i].instance + end + + # Test reconstruct! (mutating) + reconstruct!(zt, dataset_zt) + for i in 1:N + @test dataset_zt[i].x ≈ dataset[i].x atol = 1e-10 + end +end + +@testitem "Maximizers" begin + using DecisionFocusedLearningBenchmarks.Utils: TopKMaximizer + top_k = TopKMaximizer(3) + @test top_k([1, 2, 3, 4, 5]) == [0, 0, 1, 1, 1] + @test top_k([5, 4, 3, 2, 1]) == [1, 1, 1, 0, 0] + @test_throws(AssertionError, top_k([1, 2])) +end