matbesancon / MathOptSetDistances.jl

Distances to sets for MathOptInterface
MIT License
24 stars 4 forks source link

Checking if a point is a feasible solution to a JuMP model #40

Open oxinabox opened 3 years ago

oxinabox commented 3 years ago

I don't know that this actually belongs here. The is_feasible(constraint_set::MOI.AbstractSet, point) probably does. but all the stuff that wrangles through a JuMP.Model probably doesn't.

But here is the code that can detect if a set of values meets the constraints.

using JuMP
using MathOptSetDistances
using MathOptSetDistances: DefaultDistance
const MOI = JuMP.MOI
using JuMP.MOI.Utilities

"""
    is_feasible(model, variable_values::AbstractDict)

Detects is `variable_values` reperesents a feasible solution to the constraints of the
`model`.
Any variables present in the model that values are not provided for are assume to take
allowable values.

 - model: a `JuMP.Model`
 - `variable_values` a dictionary mapping `variable_ref => value` where each
   `variable_ref` is a variable that is part of the `model`, and `value` is the value it is
   set to in this candidate solution point.
"""
function is_feasible(model::Model, variable_values::AbstractDict{VariableRef})
    # Check inputs are sensible
    for (variable_ref, variable_val) in variable_values
        check_belongs_to_model(variable_ref, model)
    end

    var_map = construct_varmap(variable_values)

    # now we check all the constraints
    for (f, s) in list_of_constraint_types(model)
        for constraint in constraint_object.(all_constraints(model, f, s))
            constraint_moi_func = moi_function(constraint)
            point = substitute(var_map, constraint_moi_func)
            # We have not provided a value for one of the terms, thus this is not binding
            point === nothing && continue
            constraint_set = moi_set(constraint)
            is_feasible(constraint_set, point) || return false
        end
    end
    return true
end

"""
    is_feasible(constraint_set::MOI.AbstractSet, point)

Determines if the `point` is within the `constraint_set`.
"""
function is_feasible(constraint_set::MOI.AbstractSet, point)
    return distance_to_set(DefaultDistance(), point, constraint_set) <= 0
end

"""
    substitute(varmap, moi_func::MOI.AbstractFunction)

Given a funtion `varmap` that maps from some `MOI.VariableIndex` to a value,
and a MOI Function `moi_func`:
evalutes the `moi_func`, substituting in variable values according to `varmap`/
Returns `nothing` if a `KeyError` is thrown by `varmap` indicating that no value was
provided for one ot the variables that was used.
"""
function substitute(varmap, moi_func::MOI.AbstractFunction)
    try
        return MOI.Utilities.eval_variables(varmap, moi_func)
    catch err
        err isa KeyError || rethrow()
        return nothing
    end
end

"""
    construct_varmap(variable_values::AbstractDict{VariableRef})

Given a `variable_values` a dictionary mapping `variable_ref => value` where each
`variable_ref` is a JuMP variable that is part of the `model`, and `value` being it's value,
this constructs `varmap` a function mapping from `MOI.VariableIndex` to the value.
"""
function construct_varmap(variable_values::AbstractDict{VariableRef})
    moi_variable_values = Dict(JuMP.index(k) => v for (k, v) in variable_values)
    return k->getindex(moi_variable_values, k)
end

Demo

using Test

model = Model()
@variable(model, 0 <= y <= 30)
@variable(model, x <= 6)
@constraint(model, x <= y)

@variable(model, z)
m = @constraint(model, z <= x^2+y^2)

@test is_feasible(model, Dict(x=>1, y=>2))  # Good
@test !is_feasible(model, Dict(x=>2, y=>1)) # x>y
@test !is_feasible(model, Dict(x=>2, y=>100)) # y>30
@test !is_feasible(model, Dict(x=>7, y=>1)) # x > 6
@test !is_feasible(model, Dict(x=>-2, y=>-1)) # y < 0

@test is_feasible(model, Dict(x=>-2)) # y not given.
@test !is_feasible(model, Dict(z=>10, y=>2, x=>1))
@test is_feasible(model, Dict(z=>10, y=>5, x=>1))
matbesancon commented 3 years ago

So a complete feasibility checker is out of the scope for this package, see https://github.com/joaquimg/FeasibilityOptInterface.jl for that.

I think the function we can implement here is

function set_member(::AbstractDistance, v::V, set::S) -> bool
# can default to verifying that distance is exactly zero
end

# can be useful to have
Base.in(v::V, s::S) = set_member(DefaultDistance(), v, set)
oxinabox commented 3 years ago

Checking a point is feasible is not feasibility checking though, it is exactly what you decribed. Possibly the name is_feasible is not good. Checking if the problem is feasible is checking if there exists any point that is feasible right? which is much harder. The code above just checks points.

matbesancon commented 3 years ago

As a first step, the set membership that you describe definitely does. For the model-level feasibility of a solution, we can implement it with MOI, using eval_variables.

Fun-fact: this set membership part was the initial motivation for the PR that was replaced by this package

oxinabox commented 3 years ago

As a first step, the set membership that you describe definitely does.

Does what?

For the model-level feasibility of a solution, we can implement it with MOI, using eval_variables.

I can't find MOI.eval_variables. Would that replace my substitute function?

matbesancon commented 3 years ago

Does fit in the scope of this package, sorry.

The function is in MOI.Utilities: https://jump.dev/MathOptInterface.jl/dev/apireference/#MathOptInterface.Utilities.eval_variables

Yes that's the substitute you have

oxinabox commented 3 years ago

Nice, I have updated the code in the top post to use that. Just cos I need somewhere to keep this code til FeasibilityOptInterface is ready.