Same expression returning inconsistent data types brings considerable inconvenience

I mean if there is a function whose input arguments are allowed to degenerate from a matrix to a vector or from a vector to a scalar. I try to give an example and sincerely ask the experienced practitioners what is the most recommended solution:

using JuMP, HiGHS
function fun(a, b, c)
    model = Model(HiGHS.Optimizer)
    @variable(model, x[eachindex(a)])
    @variable(model, y[eachindex(b)])
    @variable(model, z[eachindex(c)])
    @constraint(model, a' * x + b' * y + c' * z >= 0)
    optimize!(model)
    return value.(x) <= rand(size(a)...), value.(y) <= rand(size(b)...), value.(z) <= rand(size(c)...)
end
fun([1, 2], [3, 4], [5, 6]) # (true, true, true)
fun(1, [3, 4], [5, 6]) # ERROR: Addition between an array and a JuMP scalar is not supported: instead of `x + y`, do `x .+ y` for element-wise addition.

OK, it doesn’t matter to write the constraint in a broadcast form instead of the more accurate vector form:

function fun(a, b, c)
    model = Model(HiGHS.Optimizer)
    @variable(model, x[eachindex(a)])
    @variable(model, y[eachindex(b)])
    @variable(model, z[eachindex(c)])
    @constraint(model, a' * x .+ b' * y .+ c' * z >= 0)
    optimize!(model)
    return value.(x) <= rand(size(a)...), value.(y) <= rand(size(b)...), value.(z) <= rand(size(c)...)
end
fun([1, 2], [3, 4], [5, 6]) # (true, true, true)
fun(1, [3, 4], [5, 6]) # ERROR: MethodError: no method matching isless(::Vector{Float64}, ::Float64)

But now it brings an error similar to the [0] == 0 fault.

Then, what should I do now? Should I use some means to turn variable x into a scalar when it has only one element,

function fun(a, b, c)
    model = Model(HiGHS.Optimizer)
    x = reduce(vcat, @variable(model, x[eachindex(a)]))
    y = reduce(vcat, @variable(model, y[eachindex(b)]))
    z = reduce(vcat, @variable(model, z[eachindex(c)]))
    @constraint(model, a' * x .+ b' * y .+ c' * z >= 0)
    optimize!(model)
    return value.(x) <= rand(size(a)...), value.(y) <= rand(size(b)...), value.(z) <= rand(size(b)...)
end
fun([1, 2], [3, 4], [5, 6]) # (true, true, true)
fun(1, [3, 4], [5, 6]) # (true, true, true)

or should I use multiple dispatch (where size has to be replaced by length)?

function fun(a, b, c)
    model = Model(HiGHS.Optimizer)
    @variable(model, x[eachindex(a)])
    @variable(model, y[eachindex(b)])
    @variable(model, z[eachindex(c)])
    @constraint(model, a' * x .+ b' * y .+ c' * z >= 0)
    optimize!(model)
    return value.(x) <= rand(length(a)), value.(y) <= rand(length(b)), value.(z) <= rand(length(c)) # `size` has to be replaced by `length`
end
fun(a::Number, b, c) = fun([a], b, c)
fun(a, b::Number, c) = fun(a, [b], c)
fun(a, b, c::Number) = fun(a, b, [c]) 
fun([1, 2], [3, 4], [5, 6]) # (true, true, true)
fun(1, [3, 4], [5, 6]) # (true, true, true)

How about the function has many input arguments: fun(a, b, c, d, e, f)? And how about the function’s input arguments are matrices by default (do I have to write methods for each argument in case that it is a vector and a scalar)?

If you experience this a lot, you may be interested in the only function.

2 Likes

Are you sure you want <= used like this? It’s doing a lexicographic comparison

julia> [0,0] <= [1,1]
true

julia> [2,0] <= [1,1]
false

julia> [0,2] <= [1,1]
true

This is used sometimes, but not often. You might be looking for .<= instead, maybe combined with all. I.e., value.(x) .<= rand.() to return something of the same size as x with element-wise comparisons or all(value.(x) .<= rand.()) or all(w -> value(w) <= rand(), x) to check every inequality and return a scalar.

3 Likes

Oh…Thank you for reminding me of this. I apologize for this bad example. Perhaps my mind was confused by a series of bugs at that time. :handshake: :handshake:

Perhaps the question I really want to ask is very simple:

That is, how to add methods to a function for its degenerate forms more conveniently? I wonder how everyone deals with similar issues.

1 Like

I seldom find the need to write for degenerate forms. Broadcasting handles many cases and careful use of functions covers many of the rest. Sometimes I find the need to define such helper functions myself, but they’re usually quite simple and much easier than redefining and dispatching on arbitrary combinations of inputs to the main function.

For example, in your example you write

@constraint(model, a' * x .+ b' * y .+ c' * z >= 0)

I suspect the .+ here may be a mistake. I don’t think you want matrices to produce matrices and then broadcast-add them and then do an array inequality. Although it’s a tad hard to say for sure since the macro can transform this a bit in a way I may not fully understand.

In any case, this looks like it should be a scalar inequality. The generalization of a'*x for scalars or for vectors (where it is called the scalar product) to matrices or higher-dimensional arrays is usually dot (be sure to using LinearAlgebra to make dot accessible). So I would write this as dot(a, x) + dot(b, y) + dot(c, z) >= 0 if I wanted to handle scalars and arrays of arbitrary dimension.

If you provide a specific example of some code where you find the need to tediously handle degenerate cases, people here can probably offer some suggestions.

3 Likes

You have a few choices here. You could just declare the following. It’s now up to the user to convert. Unless you have a good idea what the user wants, this may be the best idea. You should also consider if AbstractMatrix might be appropriate, but know you cannot assume 1-based indexing in that case.

"""
    fun(a::Matrix, b::Matrix, c::Matrix, d::Matrix, e::Matrix, f::Matrix)

This function requires all input arguments to be of type `Matrix`.
"""
function fun(a::Matrix, b::Matrix, c::Matrix, d::Matrix, e::Matrix, f::Matrix)
    @info "fun arguments" a b c d e f
    # do stuff with matrices
end

If you want to help the user out more and you think you have a reasonable way to convert arguments of a different type you can add a helper function.

_to_matrix(x) = convert(Matrix, x)::Matrix
_to_matrix(x::Number) = [x;;]
_to_matrix(x::Vector) = reshape(x, length(x), 1)
"""
    fun(a, b, c, d, e, f,)

Arguments will be converted to matrices via `convert`.
Scalars will also be converted into 1x1 matrices.
Vectors will be converted into column matrices.
"""
function fun(args...)
    length(args) == 6 || throw(MethodError(fun, args))
    fun(_to_matrix.(args)...) # call the matrix version above
end

Here’s a demonstration.

julia> fun(1,2,3,4,5,6)
┌ Info: fun arguments
│   a =
│    1×1 Matrix{Int64}:
│     1
│   b =
│    1×1 Matrix{Int64}:
│     2
│   c =
│    1×1 Matrix{Int64}:
│     3
│   d =
│    1×1 Matrix{Int64}:
│     4
│   e =
│    1×1 Matrix{Int64}:
│     5
│   f =
│    1×1 Matrix{Int64}:
└     6

julia> fun(1)
ERROR: MethodError: no method matching fun(::Int64)
3 Likes

Thanks. In fact, I didn’t consider matrices in that example. I just considered vector inputs and their degenerate scalar froms.

I also think .+ is not always so proper not only here but also in some other situations. Frankly, I prefer to write mathematical code in a vectorized form once possible just as what they are in the math books or papers, and at the same time I hope they can degenerate into more simple versions as the degeneration of some main parameters. In my impression MATLAB seems to be more friendly with this issue, if it is not a wrong impression.

Thank you very much, @mkitti! I think this is really one of the best solutions and I was inspired by it a lot. :handshake:

BTW, may I ask in which situation should we have the function’s name begin with a _ like _to_matrix?

" Names sometimes are given a _ prefix (or suffix) to further suggest that something is “internal” or an implementation-detail, but it is not a rule."
https://docs.julialang.org/en/v1/manual/style-guide/

1 Like

A loose convention is that names starting with _ are meant to be private internal variables or functions. In other words, other people or packages should not use them.

1 Like