implement ranking

Author: BatyLeo

src/Ranking/Ranking.jl

@@ -0,0 +1,84 @@
+module Ranking
+
+using ..Utils
+using DocStringExtensions: TYPEDEF, TYPEDFIELDS, TYPEDSIGNATURES
+using Flux: Chain, Dense
+using Random
+
+"""
+$TYPEDEF
+
+Benchmark problem with an argmax as the CO algorithm.
+
+# Fields
+$TYPEDFIELDS
+"""
+struct RankingBenchmark <: AbstractBenchmark
+    "iinstances dimension, total number of classes"
+    instance_dim::Int
+    "number of features"
+    nb_features::Int
+end
+
+function Base.show(io::IO, bench::RankingBenchmark)
+    (; instance_dim, nb_features) = bench
+    return print(
+        io, "RankingBenchmark(instance_dim=$instance_dim, nb_features=$nb_features)"
+    )
+end
+
+function RankingBenchmark(; instance_dim::Int=10, nb_features::Int=5)
+    return RankingBenchmark(instance_dim, nb_features)
+end
+
+"""
+$TYPEDSIGNATURES
+
+Compute the vector `r` such that `rᵢ` is the rank of `θᵢ` in `θ`.
+"""
+function ranking(θ::AbstractVector; rev::Bool=false, kwargs...)
+    return invperm(sortperm(θ; rev=rev))
+end
+
+"""
+$TYPEDSIGNATURES
+
+Return a top k maximizer.
+"""
+function Utils.generate_maximizer(bench::RankingBenchmark)
+    return ranking
+end
+
+"""
+$TYPEDSIGNATURES
+
+Generate a dataset of labeled instances for the subset selection problem.
+The mapping between features and cost is identity.
+"""
+function Utils.generate_dataset(bench::RankingBenchmark, dataset_size::Int=10; seed::Int=0)
+    (; instance_dim, nb_features) = bench
+    rng = MersenneTwister(seed)
+    features = [randn(rng, Float32, nb_features, instance_dim) for _ in 1:dataset_size]
+    mapping = Chain(Dense(nb_features => 1; bias=false), vec)
+    costs = mapping.(features)
+    solutions = ranking.(costs)
+    return [
+        DataSample(; x, θ_true, y_true) for
+        (x, θ_true, y_true) in zip(features, costs, solutions)
+    ]
+end
+
+"""
+$TYPEDSIGNATURES
+
+Initialize a linear model for `bench` using `Flux`.
+"""
+function Utils.generate_statistical_model(bench::RankingBenchmark; seed=0)
+    Random.seed!(seed)
+    (; nb_features) = bench
+    return Chain(Dense(nb_features => 1; bias=false), vec)
+end
+
+export RankingBenchmark
+
+end

src/StochasticVehicleScheduling/solution/heuristic_algorithms/deterministic_vsp.jl

@@ -0,0 +1,52 @@
+"""
+    solve_deterministic_VSP(instance::Instance; include_delays=true)
+
+Return the optimal solution of the deterministic VSP problem associated to `instance`.
+The objective function is `vehicle_cost * nb_vehicles + include_delays * delay_cost * sum_of_travel_times`
+Note: If you have Gurobi, use `grb_model` as `model_builder` instead od `cbc_model`.
+"""
+function solve_deterministic_VSP(
+    instance::Instance; include_delays=true, model_builder=cbc_model
+)
+    (; city, graph) = instance
+
+    travel_times = [
+        distance(task1.end_point, task2.start_point) for task1 in city.tasks,
+        task2 in city.tasks
+    ]
+
+    model = model_builder()
+
+    nb_nodes = nv(graph)
+    job_indices = 2:(nb_nodes - 1)
+
+    @variable(model, x[i=1:nb_nodes, j=1:nb_nodes; has_edge(graph, i, j)], Bin)
+
+    @objective(
+        model,
+        Min,
+        instance.city.vehicle_cost * sum(x[1, j] for j in job_indices) +
+            include_delays *
+        instance.city.delay_cost *
+        sum(
+            travel_times[i, j] * x[i, j] for i in 1:nb_nodes for
+            j in 1:nb_nodes if has_edge(graph, i, j)
+        )
+    )
+
+    @constraint(
+        model,
+        flow[i in job_indices],
+        sum(x[j, i] for j in inneighbors(graph, i)) ==
+            sum(x[i, j] for j in outneighbors(graph, i))
+    )
+    @constraint(
+        model, demand[i in job_indices], sum(x[j, i] for j in inneighbors(graph, i)) == 1
+    )
+
+    optimize!(model)
+
+    solution = solution_from_JuMP_array(value.(x), graph)
+
+    return objective_value(model), solution
+end

src/StochasticVehicleScheduling/solution/heuristic_algorithms/local_search.jl

@@ -0,0 +1,101 @@
+"""
+    move_one_random_task!(path_value::BitMatrix, graph::AbstractGraph)
+
+Select one random (uniform) task and move it to another random (uniform) feasible vehicle
+"""
+function move_one_random_task!(path_value::BitMatrix, graph::AbstractGraph)
+    nb_tasks = size(path_value, 2)
+    selected_task = rand(DiscreteUniform(1, nb_tasks))
+    selected_vehicle = find_first_one(@view path_value[:, selected_task])
+
+    can_be_inserted = Int[]
+    # do not empty if already empty
+    empty_encountered = false  #sum(@view path_value[selected_vehicle, :]) == 1 ? true : false
+    for i in 1:nb_tasks
+        if i == selected_vehicle
+            continue
+        end
+        # else
+        is_empty = false
+        if selected_task > 1
+            before = @view path_value[i, 1:(selected_task - 1)]
+            if any(before)
+                precedent_task = selected_task - find_first_one(reverse(before))
+                if !has_edge(graph, precedent_task + 1, selected_task + 1)
+                    continue
+                end
+            elseif empty_encountered
+                continue
+            else # if !empty_encountered
+                is_empty = true
+            end
+        end
+
+        if selected_task < nb_tasks
+            after = @view path_value[i, (selected_task + 1):end]
+            if any(after)
+                next_task =
+                    selected_task +
+                    find_first_one(@view path_value[i, (selected_task + 1):end])
+                if !has_edge(graph, selected_task + 1, next_task + 1)
+                    continue
+                end
+            elseif empty_encountered
+                continue
+            elseif !empty_encountered && is_empty
+                empty_encountered = true
+            end
+        end
+
+        push!(can_be_inserted, i)
+    end
+    if length(can_be_inserted) == 0
+        @warn "No space to be inserted" selected_task path_value
+        return nothing
+    end
+    new_vehicle = rand(can_be_inserted)
+    path_value[selected_vehicle, selected_task] = false
+    path_value[new_vehicle, selected_task] = true
+    return nothing
+end
+
+"""
+    local_search(solution::Solution, instance::AbstractInstance; nb_it::Integer=100)
+
+Very simple local search heuristic, using the neighborhood defined by `move_one_random_task`
+"""
+function local_search(solution::Solution, instance::AbstractInstance; nb_it::Integer=100)
+    best_solution = copy(solution.path_value)
+    best_value = evaluate_solution(solution, instance)
+    history_x = [0]
+    history_y = [best_value]
+
+    candidate_solution = copy(solution.path_value)
+    for it in 1:nb_it
+        move_one_random_task!(candidate_solution, instance.graph)
+
+        value = evaluate_solution(candidate_solution, instance)
+        if value <= best_value # keep changes
+            best_solution = copy(candidate_solution)
+            best_value = value
+            push!(history_x, it)
+            push!(history_y, best_value)
+        else # revert changes
+            candidate_solution = copy(best_solution)
+        end
+    end
+
+    return Solution(best_solution, instance), best_value, history_x, history_y
+end
+
+"""
+    heuristic_solution(instance::AbstractInstance; nb_it=100)
+
+Very simple heuristic, using [`local_search`](@ref)
+    initialised with the solution of the deterministic Linear program
+"""
+function heuristic_solution(instance::AbstractInstance; nb_it=100)
+    _, initial_solution = solve_deterministic_VSP(instance)
+    sol, _, _, _ = local_search(initial_solution, instance; nb_it=nb_it)
+    return sol
+end

test/ranking.jl

@@ -0,0 +1,31 @@
+@testitem "Ranking" begin
+    using DecisionFocusedLearningBenchmarks
+
+    instance_dim = 10
+    nb_features = 5
+
+    b = RankingBenchmark(; instance_dim=instance_dim, nb_features=nb_features)
+
+    io = IOBuffer()
+    show(io, b)
+    @test String(take!(io)) == "RankingBenchmark(instance_dim=10, nb_features=5)"
+
+    dataset = generate_dataset(b, 50)
+    model = generate_statistical_model(b)
+    maximizer = generate_maximizer(b)
+
+    for (i, sample) in enumerate(dataset)
+        (; x, θ_true, y_true) = sample
+        @test size(x) == (nb_features, instance_dim)
+        @test length(θ_true) == instance_dim
+        @test length(y_true) == instance_dim
+        @test isnothing(sample.instance)
+        @test all(y_true .== maximizer(θ_true))
+
+        θ = model(x)
+        @test length(θ) == instance_dim
+
+        y = maximizer(θ)
+        @test length(y) == instance_dim
+    end
+end