Add support for new MathOptVRP sets

src/Hexaly.jl

@@ -74,6 +74,9 @@ include("MOI/wrapper_constraints_singlevar.jl")
 include("MOI/wrapper_constraints_linear.jl")
 include("MOI/wrapper_constraints_cp.jl")
 include("MOI/list.jl")
+include("MOI/time_windows.jl")
+include("MOI/defined_route_values.jl")
 include("MOI/sum_distances_objective.jl")
+include("MOI/copy_to.jl")
 
 end # module Hexaly

src/MOI/copy_to.jl

@@ -0,0 +1,64 @@
+# Ordered variant of MOI.Utilities.default_copy_to. Definitions of derived
+# route values must arrive after route constructors and scalar domains, but
+# before constraints and objectives that consume their output variables.
+
+_is_definition_set(::Type{<:MathOptVRP.IsEmpty}) = true
+_is_definition_set(::Type{<:MathOptVRP.SumGetIndex}) = true
+_is_definition_set(::Type) = false
+
+function _copy_priority(::Type{MOI.ConstraintIndex{F,S}}) where {F,S}
+    if F <: MOI.VariableIndex
+        return 0
+    elseif _is_definition_set(S)
+        return 1
+    end
+    return 2
+end
+
+function _copy_priority(type_pair::Tuple{Type,Type})
+    F, S = type_pair
+    return _is_definition_set(S) ? 1 : 2
+end
+
+function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike)
+    MOI.empty!(dest)
+    index_map, vis_src, constraints_not_added =
+        MOI.Utilities._copy_variables_with_set(dest, src)
+
+    MOI.Utilities.pass_attributes(dest, src, index_map, vis_src)
+
+    # Preserve default_copy_to's constrained-variable handling. Only reorder
+    # the constraints it did not use as variable constructors.
+    sort!(constraints_not_added; by = cis -> _copy_priority(eltype(cis)))
+    for cis in constraints_not_added
+        MOI.Utilities._copy_constraints(dest, src, index_map, cis)
+    end
+
+    all_types = MOI.get(src, MOI.ListOfConstraintTypesPresent())
+    nonvariable_types = filter(all_types) do (F, _)
+        return !MOI.Utilities._is_variable_function(F)
+    end
+    sort!(nonvariable_types; by = _copy_priority)
+    for type_pair in nonvariable_types
+        MOI.Utilities.pass_nonvariable_constraints(
+            dest, src, index_map, [type_pair],
+        )
+    end
+
+    # Match default_copy_to: constraint attributes are copied only after every
+    # constructor and pending constraint has populated the index map.
+    for (F, S) in all_types
+        MOI.Utilities.pass_attributes(
+            dest,
+            src,
+            index_map,
+            MOI.get(src, MOI.ListOfConstraintIndices{F,S}()),
+        )
+    end
+
+    # Delay model attributes (particularly the objective) until definitions
+    # have been processed.
+    MOI.Utilities.pass_attributes(dest, src, index_map)
+    MOI.Utilities.final_touch(dest, index_map)
+    return index_map
+end

src/MOI/defined_route_values.jl

@@ -0,0 +1,88 @@
+# MathOptVRP definition constraints. The first MOI variable is substituted by
+# a native Hexaly expression; no Hexaly decision is created for it.
+
+for SetType in (MathOptVRP.IsEmpty, MathOptVRP.SumGetIndex)
+    @eval begin
+        function MOI.supports_constraint(
+            ::Optimizer,
+            ::Type{<:Union{MOI.VectorOfVariables,MOI.VectorAffineFunction}},
+            ::Type{<:$SetType},
+        )
+            return true
+        end
+
+        function MOI.supports_add_constrained_variables(
+            ::Optimizer,
+            ::Type{<:$SetType},
+        )
+            return false
+        end
+    end
+end
+
+function _definition_output_and_route(
+    m::Optimizer,
+    f::Union{MOI.VectorOfVariables,MOI.VectorAffineFunction},
+    dimension::Int,
+    name::AbstractString,
+)
+    items = _normalize_sum_distances_items(f)
+    length(items) == dimension || error(
+        "$name expected $dimension variables; got $(length(items)).",
+    )
+    all(it -> it isa MOI.VariableIndex, items) || error(
+        "$name expects only `MOI.VariableIndex` values.",
+    )
+    output = items[1]
+    route_vars = items[2:end]
+    route = _info(m, route_vars[1]).parent_list
+    route !== nothing || error("$name route variables have no parent Hexaly list.")
+    all(_info(m, vi).parent_list === route for vi in route_vars) || error(
+        "$name route variables must belong to the same Hexaly list.",
+    )
+    return output, route
+end
+
+function MOI.add_constraint(
+    m::Optimizer,
+    f::Union{MOI.VectorOfVariables,MOI.VectorAffineFunction},
+    s::MathOptVRP.IsEmpty,
+)
+    output, route = _definition_output_and_route(
+        m, f, MOI.dimension(s), "MathOptVRP.IsEmpty",
+    )
+    info = _info(m, output)
+    info.is_binary || error(
+        "MathOptVRP.IsEmpty output must be created in `MOI.ZeroOne`.",
+    )
+    _define_expression!(m, output, eq(m.model, count_(m.model, route), 0))
+    cindex = MOI.ConstraintIndex{typeof(f),typeof(s)}(
+        length(m.constraint_info) + 1,
+    )
+    m.constraint_info[cindex] = ConstraintInfo(cindex, nothing, f, s)
+    return cindex
+end
+
+function MOI.add_constraint(
+    m::Optimizer,
+    f::Union{MOI.VectorOfVariables,MOI.VectorAffineFunction},
+    s::MathOptVRP.SumGetIndex,
+)
+    output, route = _definition_output_and_route(
+        m, f, MOI.dimension(s), "MathOptVRP.SumGetIndex",
+    )
+    terms = HxExpression[
+        prod(m.model, s.values[i], contains_(m.model, route, i - 1))
+        for i in eachindex(s.values) if !iszero(s.values[i])
+    ]
+    expression = isempty(terms) ? create_constant(m.model, 0) :
+        sum(m.model, terms...)
+    info = _info(m, output)
+    info.is_integer = all(isinteger, s.values)
+    _define_expression!(m, output, expression)
+    cindex = MOI.ConstraintIndex{typeof(f),typeof(s)}(
+        length(m.constraint_info) + 1,
+    )
+    m.constraint_info[cindex] = ConstraintInfo(cindex, nothing, f, s)
+    return cindex
+end

src/MOI/list.jl

@@ -121,130 +121,122 @@ function MOI.add_constrained_variables(m::Optimizer, set::MathOptVRP.PartitionPD
     return indices, cindex
 end
 
-# ── MathOptVRP.TimeWindows ─────────────────────────────────────────────
-# Layout: `[t; depot_start; nodes...; depot_end]`. Posts the per-customer
-# time-window constraint and the makespan linkage `t >= total_time`.
-
-function MOI.supports_constraint(
-    ::Optimizer,
-    ::Type{<:Union{MOI.VectorOfVariables,MOI.VectorAffineFunction}},
-    ::Type{<:MathOptVRP.TimeWindows},
+function _mathoptvrp_parent_list(
+    m::Optimizer,
+    f::Union{MOI.VectorOfVariables,MOI.VectorAffineFunction},
+    dimension::Int,
+    set_name::AbstractString,
 )
-    return true
+    items = _normalize_sum_distances_items(f)
+    length(items) == dimension || error(
+        "$set_name expected $dimension route-node variables; got $(length(items)).",
+    )
+    _shift_to_zero_based!(items)
+    all(it -> it isa MOI.VariableIndex, items) || error(
+        "$set_name: every item must be a `MOI.VariableIndex` backed by a Hexaly list.",
+    )
+    seq = _info(m, items[1]).parent_list
+    seq !== nothing || error("$set_name: variables have no parent Hexaly list.")
+    all(_info(m, vi).parent_list === seq for vi in items) || error(
+        "$set_name: all node variables must belong to the same Hexaly list.",
+    )
+    return seq
 end
 
-# These sets constrain node variables that a `Partition` / `PartitionPD`
-# already created; they cannot create any. Without this, MOI's default
-# infers `supports_add_constrained_variables` from the `VectorOfVariables`
-# constraint above and `copy_to` builds the variables from this set instead
-# of from the partition, leaving them without a backing Hexaly list.
-function MOI.supports_add_constrained_variables(
-    ::Optimizer,
-    ::Type{<:MathOptVRP.TimeWindows},
+for SetType in (
+    MathOptVRP.RouteCompatibility,
+    MathOptVRP.RouteOrder,
+    MathOptVRP.RouteExtremities,
 )
-    return false
+    @eval begin
+        function MOI.supports_constraint(
+            ::Optimizer,
+            ::Type{<:Union{MOI.VectorOfVariables,MOI.VectorAffineFunction}},
+            ::Type{<:$SetType},
+        )
+            return true
+        end
+
+        function MOI.supports_add_constrained_variables(
+            ::Optimizer,
+            ::Type{<:$SetType},
+        )
+            return false
+        end
+    end
 end
 
 function MOI.add_constraint(
     m::Optimizer,
     f::Union{MOI.VectorOfVariables,MOI.VectorAffineFunction},
-    s::MathOptVRP.TimeWindows,
+    s::MathOptVRP.RouteCompatibility,
 )
-    items = _normalize_sum_distances_items(f)
-    length(items) == MOI.dimension(s) || error(
-        "MathOptVRP.TimeWindows expected `length([t; depot_start; nodes; depot_end]) ",
-        "== $(MOI.dimension(s))`; got $(length(items)).",
-    )
-    # Item 1 is the user's `t` variable, items 2..end are node values.
-    _shift_to_zero_based!(items, 2)
-    items[1] isa MOI.VariableIndex || error(
-        "MathOptVRP.TimeWindows: first item must be the total-time `t` variable.",
-    )
-    items[2] isa Real || error(
-        "MathOptVRP.TimeWindows: second item must be the constant `depot_start` index.",
-    )
-    items[end] isa Real || error(
-        "MathOptVRP.TimeWindows: last item must be the constant `depot_end` index.",
-    )
-    t_vi = items[1]
-    depot_start = round(Int, items[2])
-    depot_end = round(Int, items[end])
-    var_items = @view items[3:(end-1)]
-    all(it -> it isa MOI.VariableIndex, var_items) || error(
-        "MathOptVRP.TimeWindows: items 3..end-1 must be node `MOI.VariableIndex` values ",
-        "backed by a Hexaly list.",
-    )
-    first_pl = _info(m, var_items[1]).parent_list
-    first_pl !== nothing || error(
-        "MathOptVRP.TimeWindows: node variables have no parent Hexaly list.",
+    seq = _mathoptvrp_parent_list(
+        m, f, MOI.dimension(s), "MathOptVRP.RouteCompatibility",
     )
-    all(_info(m, vi).parent_list === first_pl for vi in var_items) || error(
-        "MathOptVRP.TimeWindows: all node variables must belong to the same Hexaly list.",
+    for node in findall(!, s.allowed)
+        _add_hexaly_constraint!(m, eq(m.model, contains_(m.model, seq, node - 1), 0))
+    end
+    cindex = MOI.ConstraintIndex{typeof(f),typeof(s)}(
+        length(m.constraint_info) + 1,
     )
+    m.constraint_info[cindex] = ConstraintInfo(cindex, nothing, f, s)
+    return cindex
+end
 
-    md = m.model
-    t_var = _info(m, t_vi).variable
-    seq = first_pl
-    c = count_(md, seq)
-    service = round(Int, s.service)
-    n_rows = size(s.travel, 1)
-    dist_arr = array(md,
-        [array(md, round.(Int, s.travel[i, :])) for i = 1:n_rows])
-    earliest_arr = array(md, round.(Int, s.earliest))
-    latest_arr = array(md, round.(Int, s.latest))
-
-    end_time = array(md,
-        range_(md, 0, c),
-        lambda_function(md,
-            (i, prev) -> iif(md,
-                eq(md, i, 0),
-                sum(md,
-                    max(md,
-                        at(md, earliest_arr, at(md, seq, 0)),
-                        at(md, dist_arr, depot_start, at(md, seq, 0)),
-                    ),
-                    service,
-                ),
-                sum(md,
-                    max(md,
-                        at(md, earliest_arr, at(md, seq, i)),
-                        sum(md, prev, at(md, dist_arr, at(md, seq, sub(md, i, 1)), at(md, seq, i))),
-                    ),
-                    service,
-                ),
-            ); nargs = 2,
-        ),
-        0,
+function MOI.add_constraint(
+    m::Optimizer,
+    f::Union{MOI.VectorOfVariables,MOI.VectorAffineFunction},
+    s::MathOptVRP.RouteOrder,
+)
+    seq = _mathoptvrp_parent_list(
+        m, f, MOI.dimension(s), "MathOptVRP.RouteOrder",
     )
-
-    _add_hexaly_constraint!(m,
-        and_(md,
-            range_(md, 0, c),
-            lambda_function(md,
-                i -> leq(md,
-                    sub(md, at(md, end_time, i), service),
-                    at(md, latest_arr, at(md, seq, i)),
-                ); nargs = 1,
-            ),
-        ),
+    md = m.model
+    for p in findall(s.before), q in findall(s.after)
+        both = and_(md, contains_(md, seq, p - 1), contains_(md, seq, q - 1))
+        _add_hexaly_constraint!(m, or_(md, not_(md, both),
+            lt(md, index_of(md, seq, p - 1), index_of(md, seq, q - 1))))
+    end
+    cindex = MOI.ConstraintIndex{typeof(f),typeof(s)}(
+        length(m.constraint_info) + 1,
     )
+    m.constraint_info[cindex] = ConstraintInfo(cindex, nothing, f, s)
+    return cindex
+end
 
-    total_time = iif(md,
-        gt(md, c, 0),
-        sum(md, at(md, end_time, sub(md, c, 1)),
-            at(md, dist_arr, at(md, seq, sub(md, c, 1)), depot_end)),
-        0,
+function MOI.add_constraint(
+    m::Optimizer,
+    f::Union{MOI.VectorOfVariables,MOI.VectorAffineFunction},
+    s::MathOptVRP.RouteExtremities,
+)
+    seq = _mathoptvrp_parent_list(
+        m, f, MOI.dimension(s), "MathOptVRP.RouteExtremities",
     )
-    _add_hexaly_constraint!(m, geq(md, t_var, total_time))
-
+    md = m.model
+    others = findall(!, s.members)
+    for p in findall(s.members)
+        isempty(others) && continue
+        pairwise_before = [or_(md, not_(md, contains_(md, seq, q - 1)),
+            lt(md, index_of(md, seq, p - 1), index_of(md, seq, q - 1)))
+            for q in others]
+        pairwise_after = [or_(md, not_(md, contains_(md, seq, q - 1)),
+            lt(md, index_of(md, seq, q - 1), index_of(md, seq, p - 1)))
+            for q in others]
+        before = length(pairwise_before) == 1 ? pairwise_before[1] :
+            and_(md, pairwise_before...)
+        after = length(pairwise_after) == 1 ? pairwise_after[1] :
+            and_(md, pairwise_after...)
+        _add_hexaly_constraint!(m, or_(md,
+            not_(md, contains_(md, seq, p - 1)), before, after))
+    end
     cindex = MOI.ConstraintIndex{typeof(f),typeof(s)}(
         length(m.constraint_info) + 1,
     )
     m.constraint_info[cindex] = ConstraintInfo(cindex, nothing, f, s)
     return cindex
 end
 
-# ── MathOptVRP.Capacity ────────────────────────────────────────────────
 
 function MOI.supports_constraint(
     ::Optimizer,
@@ -388,7 +380,7 @@ function MOI.add_constraint(
     )
 
     md = m.model
-    t_var = _info(m, t_vi).variable
+    t_var = _expression!(m, t_vi)
     seq = first_pl
     c = count_(md, seq)
     capacity = round(Int, s.capacity)