diff --git a/docs/src/apireference.md b/docs/src/apireference.md index dfecf81c80..e09ba87928 100644 --- a/docs/src/apireference.md +++ b/docs/src/apireference.md @@ -660,6 +660,7 @@ Utilities.mode The following utilities are available for functions: ```@docs Utilities.eval_variables +Utilities.substitute_variables Utilities.remove_variable Utilities.all_coefficients Utilities.unsafe_add @@ -680,6 +681,7 @@ Utilities.scalar_type Utilities.promote_operation Utilities.operate Utilities.operate! +Utilities.operate_output_index! Utilities.vectorize ``` diff --git a/src/Utilities/functions.jl b/src/Utilities/functions.jl index c010cb7ece..466c77cdf0 100644 --- a/src/Utilities/functions.jl +++ b/src/Utilities/functions.jl @@ -3,7 +3,10 @@ using Test """ eval_variables(varval::Function, f::AbstractFunction) -Returns the value of function `f` if each variable index `vi` is evaluated as `varval(vi)`. +Returns the value of function `f` if each variable index `vi` is evaluated as +`varval(vi)`. Note that `varval` should return a number, see +[`substitute_variables`](@ref) for a similar function where `varval` returns a +function. """ function eval_variables end eval_variables(varval::Function, f::SVF) = varval(f.variable) @@ -83,6 +86,97 @@ function mapvariables(varmap, f::MOI.AbstractFunctionModification) return mapvariables(vi -> varmap[vi], f) end +# For performance reason, we assume that the type of the function does not +# change in `substitute_variables`. +""" + substitute_variables(variable_map::Function, x) + +Substitute any [`MOI.VariableIndex`](@ref) in `x` by `variable_map(x)`. The +`variable_map` function returns either [`MOI.SingleVariable`](@ref) or +[`MOI.ScalarAffineFunction`](@ref), see [`eval_variables`](@ref) for a similar +function where `variable_map` returns a number. + +This function is used by bridge optimizers on constraint functions, attribute +values and submittable values when at least one variable bridge is used hence it +needs to be implemented for custom types that are meant to be used as attribute +or submittable value. +""" +function substitute_variables end + +function substitute_variables( + ::Function, x::Union{Number, Enum, AbstractArray{<:Union{Number, Enum}}}) + return x +end + +substitute_variables(::Function, set::MOI.AbstractSet) = set + +function substitute_variables(variable_map::Function, + term::MOI.ScalarQuadraticTerm{T}) where T + f1 = variable_map(term.variable_index_1) + f2 = variable_map(term.variable_index_2) + f12 = operate(*, T, f1, f2)::MOI.ScalarQuadraticFunction{T} + coef = term.coefficient + # The quadratic terms are evaluated as x'Qx/2 so a diagonal term should + # be divided by 2 while an off-diagonal term appears twice in the matrix + # and is divided by 2 so it stays the same. + if term.variable_index_1 == term.variable_index_2 + coef /= 2 + end + return operate!(*, T, f12, coef) +end +function substitute_variables(variable_map::Function, + term::MOI.ScalarAffineTerm{T}) where T + return operate(*, T, term.coefficient, + variable_map(term.variable_index))::MOI.ScalarAffineFunction{T} +end +function substitute_variables( + variable_map::Function, + func::MOI.ScalarAffineFunction{T}) where T + g = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{T}[], MOI.constant(func)) + for term in func.terms + operate!(+, T, g, substitute_variables(variable_map, term))::typeof(func) + end + return g +end +function substitute_variables( + variable_map::Function, + func::MOI.VectorAffineFunction{T}) where T + g = MOI.VectorAffineFunction(MOI.VectorAffineTerm{T}[], copy(MOI.constant(func))) + for term in func.terms + sub = substitute_variables(variable_map, term.scalar_term) + operate_output_index!(+, T, term.output_index, g, sub)::typeof(func) + end + return g +end +function substitute_variables( + variable_map::Function, + func::MOI.ScalarQuadraticFunction{T}) where T + g = MOI.ScalarQuadraticFunction( + MOI.ScalarAffineTerm{T}[], MOI.ScalarQuadraticTerm{T}[], MOI.constant(func)) + for term in func.affine_terms + operate!(+, T, g, substitute_variables(variable_map, term))::typeof(func) + end + for term in func.quadratic_terms + operate!(+, T, g, substitute_variables(variable_map, term))::typeof(func) + end + return g +end +function substitute_variables( + variable_map::Function, + func::MOI.VectorQuadraticFunction{T}) where T + g = MOI.VectorQuadraticFunction( + MOI.VectorAffineTerm{T}[], MOI.VectorQuadraticTerm{T}[], copy(MOI.constant(func))) + for term in func.affine_terms + sub = substitute_variables(variable_map, term.scalar_term) + operate_output_index!(+, T, term.output_index, g, sub)::typeof(func) + end + for term in func.quadratic_terms + sub = substitute_variables(variable_map, term.scalar_term) + operate_output_index!(+, T, term.output_index, g, sub)::typeof(func) + end + return g +end + # Vector of constants constant_vector(f::Union{SAF, SQF}) = [f.constant] constant_vector(f::Union{VAF, VQF}) = f.constants @@ -597,6 +691,20 @@ can be modified. The return type is the same than the method """ function operate! end +""" + operate_output_index!( + op::Function, ::Type{T}, output_index::Integer, + func::MOI.AbstractVectorFunction + args::Union{T, MOI.AbstractScalarFunction}...)::MOI.AbstractFunction where T + +Returns an `MOI.AbstractVectorFunction` where the function at `output_index` +is the result of the operation `op` applied to the function at `output_index` +of `func` and `args`. The functions at output index different to `output_index` +are the same as the functions at the same output index in `func`. The first +argument can be modified. +""" +function operate_output_index! end + """ promote_operation(op::Function, ::Type{T}, ArgsTypes::Type{<:Union{T, MOI.AbstractFunction}}...) where T @@ -751,8 +859,7 @@ end ## operate! # + with at least 3 arguments function operate!(op::typeof(+), ::Type{T}, f, g, h, args...) where T - operate!(op, T, f, g) - return operate!(+, T, f, h, args...) + return operate!(+, T, operate!(op, T, f, g), h, args...) end # Unary - @@ -967,6 +1074,13 @@ function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T}, return operate(op, T, f, g) end # Vector Affine +/-! ... +function operate_output_index!( + op::Union{typeof(+), typeof(-)}, ::Type{T}, + output_index::Integer, + f::MOI.VectorAffineFunction{T}, α::T) where T + f.constants[output_index] = op(f.constants[output_index], α) + return f +end function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T}, f::MOI.VectorAffineFunction{T}, g::Vector{T}) where T @@ -984,6 +1098,15 @@ function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T}, MOI.ScalarAffineTerm.(op(one(T)), g.variables))) return f end +function operate_output_index!( + op::Union{typeof(+), typeof(-)}, ::Type{T}, + output_index::Integer, + f::MOI.VectorAffineFunction{T}, + g::MOI.ScalarAffineFunction{T}) where T + append!(f.terms, MOI.VectorAffineTerm.( + output_index, operate_terms(op, g.terms))) + return operate_output_index!(op, T, output_index, f, MOI.constant(g)) +end function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T}, f::MOI.VectorAffineFunction{T}, g::MOI.VectorAffineFunction{T}) where T @@ -997,6 +1120,13 @@ function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T}, return operate(op, T, f, g) end # Vector Quadratic +/-! ... +function operate_output_index!( + op::Union{typeof(+), typeof(-)}, ::Type{T}, + output_index::Integer, + f::MOI.VectorQuadraticFunction{T}, α::T) where T + f.constants[output_index] = op(f.constants[output_index], α) + return f +end function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T}, f::MOI.VectorQuadraticFunction{T}, g::Vector{T}) where T @@ -1012,6 +1142,15 @@ function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T}, append!(f.affine_terms, MOI.VectorAffineTerm.(collect(1:d), MOI.ScalarAffineTerm.(op(one(T)), g.variables))) return f end +function operate_output_index!( + op::Union{typeof(+), typeof(-)}, ::Type{T}, + output_index::Integer, + f::MOI.VectorQuadraticFunction{T}, + g::MOI.ScalarAffineFunction{T}) where T + append!(f.affine_terms, MOI.VectorAffineTerm.( + output_index, operate_terms(op, g.terms))) + return operate_output_index!(op, T, output_index, f, MOI.constant(g)) +end function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T}, f::MOI.VectorQuadraticFunction{T}, g::MOI.VectorAffineFunction{T}) where T @@ -1019,6 +1158,17 @@ function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T}, f.constants .= op.(f.constants, g.constants) return f end +function operate_output_index!( + op::Union{typeof(+), typeof(-)}, ::Type{T}, + output_index::Integer, + f::MOI.VectorQuadraticFunction{T}, + g::MOI.ScalarQuadraticFunction{T}) where T + append!(f.affine_terms, MOI.VectorAffineTerm.( + output_index, operate_terms(op, g.affine_terms))) + append!(f.quadratic_terms, MOI.VectorQuadraticTerm.( + output_index, operate_terms(op, g.quadratic_terms))) + return operate_output_index!(op, T, output_index, f, MOI.constant(g)) +end function operate!(op::Union{typeof(+), typeof(-)}, ::Type{T}, f::MOI.VectorQuadraticFunction{T}, g::MOI.VectorQuadraticFunction{T}) where T @@ -1192,22 +1342,29 @@ end function operate(::typeof(*), ::Type{T}, f::MOI.SingleVariable, g::MOI.SingleVariable) where T - return MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{T}[], - [MOI.ScalarQuadraticTerm(one(T), - f.variable, - g.variable)], - zero(T)) + return MOI.ScalarQuadraticFunction( + MOI.ScalarAffineTerm{T}[], + [MOI.ScalarQuadraticTerm(f.variable == g.variable ? 2one(T) : one(T), + f.variable, g.variable)], + zero(T)) end function operate(::typeof(*), ::Type{T}, f::MOI.ScalarAffineFunction{T}, g::MOI.SingleVariable) where T - aff_terms = [MOI.ScalarAffineTerm(f.constant, g.variable)] - quad_terms = map(t -> MOI.ScalarQuadraticTerm(t.coefficient, - t.variable_index, - g.variable), - f.terms) + if iszero(f.constant) + aff_terms = MOI.ScalarAffineTerm{T}[] + else + aff_terms = [MOI.ScalarAffineTerm(f.constant, g.variable)] + end + quad_terms = map(t -> MOI.ScalarQuadraticTerm( + t.variable_index == g.variable ? 2t.coefficient : t.coefficient, + t.variable_index, g.variable), f.terms) return MOI.ScalarQuadraticFunction(aff_terms, quad_terms, zero(T)) end +function operate(::typeof(*), ::Type{T}, f::MOI.SingleVariable, + g::MOI.ScalarAffineFunction{T}) where T + return operate(*, T, g, f) +end function operate(::typeof(*), ::Type{T}, f::MOI.ScalarAffineFunction{T}, g::MOI.ScalarAffineFunction{T}) where T diff --git a/src/sets.jl b/src/sets.jl index 9ff4b316ea..b396b8531a 100644 --- a/src/sets.jl +++ b/src/sets.jl @@ -5,7 +5,9 @@ """ AbstractSet -Abstract supertype for set objects used to encode constraints. +Abstract supertype for set objects used to encode constraints. A set object +should not contain any [`VariableIndex`](@ref) or [`ConstraintIndex`](@ref) +as the set is passed unmodifed during [`copy_to`](@ref). """ abstract type AbstractSet end diff --git a/test/Utilities/functions.jl b/test/Utilities/functions.jl index ce12ad3035..9b4f101855 100644 --- a/test/Utilities/functions.jl +++ b/test/Utilities/functions.jl @@ -4,9 +4,13 @@ const MOI = MathOptInterface const MOIU = MOI.Utilities w = MOI.VariableIndex(0) +fw = MOI.SingleVariable(w) x = MOI.VariableIndex(1) +fx = MOI.SingleVariable(x) y = MOI.VariableIndex(2) +fy = MOI.SingleVariable(y) z = MOI.VariableIndex(3) +fz = MOI.SingleVariable(z) @testset "Vectorization" begin g = MOI.VectorAffineFunction(MOI.VectorAffineTerm.([3, 1], MOI.ScalarAffineTerm.([5, 2], @@ -105,6 +109,37 @@ end @test MOIU.eval_variables(vi -> vals[vi], fvq) ≈ [13, 1] @test MOIU.eval_variables(vi -> vals[vi], fvq) ≈ [13, 1] end +@testset "substitute_variables" begin + # We do tests twice to make sure the function is not modified + subs = Dict(w => 1.0fy + 1.0fz, x => 2.0fy + 1.0, y => 1.0fy, z => -1.0fw) + vals = Dict(w => 0.0, x => 3.0, y => 1.0, z => 5.0) + subs_vals = Dict(w => 6.0, x => 3.0, y => 1.0, z => 0.0) + fsa = fx + 3.0fz + 2.0fy + 2.0 + subs_sa = -3.0fw + 4.0fy + 3.0 + @test MOIU.eval_variables(vi -> subs_vals[vi], fsa) == MOIU.eval_variables(vi -> vals[vi], subs_sa) + @test MOIU.substitute_variables(vi -> subs[vi], fsa) ≈ subs_sa + @test MOIU.substitute_variables(vi -> subs[vi], fsa) ≈ subs_sa + fva = MOIU.operate(vcat, Float64, 3.0fz - 3.0, fx + 2.0fy + 2.0) + subs_va = MOIU.operate(vcat, Float64, -3.0fw - 3.0, 4.0fy + 3.0) + @test MOIU.eval_variables(vi -> subs_vals[vi], fva) == MOIU.eval_variables(vi -> vals[vi], subs_va) + @test MOIU.substitute_variables(vi -> subs[vi], fva) ≈ subs_va + @test MOIU.substitute_variables(vi -> subs[vi], fva) ≈ subs_va + fsq = 1.0fx + 1.0fy + 1.0fx * fz + 1.0fw * fz + 1.0fw * fy + 2.0fw * fw - 3.0 + # 2y + 1 + y - 2yw - w + (y + z) * (-w) + (y + z) * y + 2 (y + z) * (y + z) - 3 + # 2y + y - 2yw - w - yw - zw + y^2 + zy + 2y^2 + 2z^2 + 4yz - 2 + # 3y - w + 3y^2 + 2z^2 - 3yw - zw + 5yz - 2 + subs_sq = 3.0fy - 1.0fw + 3.0fy * fy + 2.0fz * fz - 3.0fy * fw - 1.0fz * fw + 5.0fy * fz - 2.0 + @test MOIU.eval_variables(vi -> subs_vals[vi], fsq) == MOIU.eval_variables(vi -> vals[vi], subs_sq) + @test MOIU.substitute_variables(vi -> subs[vi], fsq) ≈ subs_sq + @test MOIU.substitute_variables(vi -> subs[vi], fsq) ≈ subs_sq + fvq = MOIU.operate(vcat, Float64, 1.0fy + 1.0fx*fz - 3.0, 1.0fx + 1.0fw*fz + 1.0fw*fy - 2.0) + subs_vq = MOIU.operate(vcat, Float64, + 1.0fy - fw - 2.0fy*fw - 3.0, + 2.0fy + 1.0fy*fy - 1.0fw*fy - 1.0fw*fz + 1.0fy*fz - 1.0) + @test MOIU.eval_variables(vi -> subs_vals[vi], fvq) == MOIU.eval_variables(vi -> vals[vi], subs_vq) + @test MOIU.substitute_variables(vi -> subs[vi], fvq) ≈ subs_vq + @test MOIU.substitute_variables(vi -> subs[vi], fvq) ≈ subs_vq +end @testset "mapvariables" begin fsq = MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm.(1.0, [x, y]), MOI.ScalarQuadraticTerm.(1.0, [x, w, w], [z, z, y]), -3.0) @@ -336,8 +371,6 @@ end MOI.ScalarQuadraticFunction{Float64}, Float64) == MOI.ScalarQuadraticFunction{Float64} end - fx = MOI.SingleVariable(x) - fy = MOI.SingleVariable(y) f = 7 + 3fx + 1fx * fx + 2fy * fy + 3fx * fy MOIU.canonicalize!(f) @test MOI.output_dimension(f) == 1 @@ -365,6 +398,29 @@ end @test convert(MOI.SingleVariable, g) == fx end @testset "operate" begin + @testset "No zero affine term" begin + # Test that no affine 0y term is created when multiplying 1fx by fy + for fxfy in [1fx * fy, fx * 1fy] + @test isempty(fxfy.affine_terms) + @test length(fxfy.quadratic_terms) == 1 + @test fxfy.quadratic_terms[1] == MOI.ScalarQuadraticTerm(1, x, y) || + fxfy.quadratic_terms[1] == MOI.ScalarQuadraticTerm(1, y, x) + end + for fxfx in [1fx * fx, fx * 1fx] + @test isempty(fxfx.affine_terms) + @test length(fxfx.quadratic_terms) == 1 + @test fxfx.quadratic_terms[1] == MOI.ScalarQuadraticTerm(2, x, x) + end + end + @testset "operate!" begin + q = 1.0fx + 1.0fy + (1.0fx) * fz + (1.0fw) * fz + @test q ≈ 1.0fx + 1.0fy + (1.0fw) * fz + (1.0fx) * fz + # This calls + aff = 1.0fx + 1.0fy + # which tries to mutate `aff`, gets a quadratic expression + # and mutate it with the remaining term + @test MOIU.operate!(+, Float64, aff, (1.0fx) * fz, (1.0fw) * fz) ≈ q + end @test f ≈ 7 + (fx + 2fy) * (1fx + fy) + 3fx @test f ≈ -(-7 - 3fx) + (fx + 2fy) * (1fx + fy) @test f ≈ -((fx + 2fy) * (MOIU.operate(-, Int, fx) - fy)) + 3fx + 7 @@ -372,46 +428,46 @@ end @test f ≈ (fx + 2) * (fx + 1) + (fy + 1) * (2fy + 3fx) + (5 - 3fx - 2fy) @test f ≈ begin MOI.ScalarQuadraticFunction([MOI.ScalarAffineTerm(3, x)], - MOI.ScalarQuadraticTerm.([1], [x], [x]), 4) + + MOI.ScalarQuadraticTerm.([2], [x], [x]), 4) + MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Int}[], - MOI.ScalarQuadraticTerm.([2, 3], [y, x], [y, y]), 3) + MOI.ScalarQuadraticTerm.([4, 3], [y, x], [y, y]), 3) end @test f ≈ begin MOI.ScalarQuadraticFunction([MOI.ScalarAffineTerm(3, x)], - MOI.ScalarQuadraticTerm.([1], [x], [x]), 10) - + MOI.ScalarQuadraticTerm.([2], [x], [x]), 10) - MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Int}[], - MOI.ScalarQuadraticTerm.([-2, -3], [y, x], [y, y]), 3) + MOI.ScalarQuadraticTerm.([-4, -3], [y, x], [y, y]), 3) end @test f ≈ begin MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(3, x)], 5) + MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Int}[], - MOI.ScalarQuadraticTerm.([1, 2, 3], [x, y, x], [x, y, y]), 2) + MOI.ScalarQuadraticTerm.([2, 4, 3], [x, y, x], [x, y, y]), 2) end @test f ≈ begin MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(3, x)], 5) - MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Int}[], - MOI.ScalarQuadraticTerm.([-1, -2, -3], [x, y, x], [x, y, y]), -2) + MOI.ScalarQuadraticTerm.([-2, -4, -3], [x, y, x], [x, y, y]), -2) end @test f ≈ begin MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Int}[], - MOI.ScalarQuadraticTerm.([1, 2, 3], [x, y, x], [x, y, y]), 2) + + MOI.ScalarQuadraticTerm.([2, 4, 3], [x, y, x], [x, y, y]), 2) + MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(3, x)], 5) end @test f ≈ begin MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Int}[], - MOI.ScalarQuadraticTerm.([1, 2, 3], [x, y, x], [x, y, y]), 12) - + MOI.ScalarQuadraticTerm.([2, 4, 3], [x, y, x], [x, y, y]), 12) - MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(-3, x)], 5) end @test f ≈ begin MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm.([2], [x]), - MOI.ScalarQuadraticTerm.([1, 2, 3], [x, y, x], [x, y, y]), 7) + + MOI.ScalarQuadraticTerm.([2, 4, 3], [x, y, x], [x, y, y]), 7) + MOI.SingleVariable(x) end @test f ≈ MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm.([3], [x]), - MOI.ScalarQuadraticTerm.([1, 2, 3], [x, y, x], [x, y, y]), 10) - 3 - @test f ≈ - 2.0 * MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm.([3.0], [x]), - MOI.ScalarQuadraticTerm.([1.0, 2.0, 3.0], [x, y, x], [x, y, y]), 7.0) / 2.0 + MOI.ScalarQuadraticTerm.([2, 4, 3], [x, y, x], [x, y, y]), 10) - 3 + @test f ≈ 2.0 * MOI.ScalarQuadraticFunction( + MOI.ScalarAffineTerm.([3.0], [x]), + MOI.ScalarQuadraticTerm.([2.0, 4.0, 3.0], [x, y, x], [x, y, y]), 7.0) / 2.0 end @testset "modification" begin f = MOIU.modify_function(f, MOI.ScalarConstantChange(9))