Make things work for general AbstractArrays

[tkoolen]

A lot of functions currently require Array types, while they would work just fine for any other subtype of AbstractArray. This PR simply replaces most ::Arrays in function signatures by ::AbstractArray (and likewise for Vector and Matrix).

I defined norm only for AbstractVector and AbstractMatrix, not for general AbstractArrays, because Base only defines it for AbstractVector and AbstractMatrix, and defining it for AbstractArray results in an ambiguity error.

Please be extra careful reviewing my changes to _multiply! and _multiplyt! in operators.jl.

This came up while I was playing around with storing JuMP Variables in AxisArrays for a trajectory optimization style problem, because I had trouble keeping track of the order of the dimensions of the variable arrays (this is proving to be quite a nice combination of tools by the way).

[tkoolen]

cc: @rdeits

[mlubin]

Not needed, use getvalue.(arr) instead.

[tkoolen]

So delete the method?

[mlubin]

It's been previously discussed. The conclusion for now is to allow getvalue on objects that JuMP creates, and JuMP can create Array{Variable}, e.g., PSD matrices.

[mlubin]

Also not needed, use setvalue.(set,val)

[mlubin]

I don't think this is needed

[mlubin]

This is an internal function, it's unclear how it could be called with an AbstractArray

[mlubin]

Internal, won't be called with AbstractVector

[mlubin]

Internal, won't be called with AbstractVector

[mlubin]

Internal, won't be called with AbstractVector

[mlubin]

Internal, won't be called with AbstractVector

[mlubin]

Internal, won't be called with AbstractVector

[mlubin]

Internal, won't be called with AbstractVector

[mlubin]

Internal, won't be called with AbstractVector

[mlubin]

None of the changes in this file are needed

[mlubin]

These methods are designed for dense matrices and shouldn't be widened

[mlubin]

These changes don't seem relevant

[mlubin]

I'm willing to consider this use case, but the changes here are way too broad and touch a number of internal functions that would never be called with AbstractArray. Could you try to make the changes minimal with respect to what you need, and also add unit tests covering them?

[tkoolen]

Sure. Thanks for the quick review. I'd like to argue that AbstractArray should always be the default (unless things are specialized for dense arrays, like in _multiply!, of course), and restricting to Array is kind of arbitrary, even in internal code that is currently not being called with other types. But I'll try to make the changes minimal and add tests.

[mlubin]

Also getvalue.(arr) here

[rdeits]

To motivate why this would be nice to have, here are some cool things that @tkoolen and I have been doing w.r.t. AxisArrays of JuMP variables: https://gist.github.com/rdeits/3371af92c9c4aa61353f38c9a91c02e9

[mlubin]

AxisArrays look like a potential replacement for JuMPArrays.

[mlubin]

restricting to Array is kind of arbitrary, even in internal code that is currently not being called with other types.

Personally I don't see much value in putting AbstractArray over no type restriction at all if we're already going to loosen the type beyond what's actually used and has been tested.

[rdeits]

AxisArrays do seem like a great replacement for the categorical indexing in JuMPArray. I think they don't cover the ragged and non-1-based indexing that JuMPArray supports, though.

[tkoolen]

I've addressed your comments. If there are more places where you want things to be more specific, let me know.

I limited the test to cases that resulted in MethodErrors before this PR. There's quite a few cases that 'work' before this PR, but don't use the (I guess optimized) JuMP code, even though it would be applicable. I'm not sure how to test for these. For example,

m = Model()
x = view(@variable(m, [1: 3]), :)
A = rand(3, 3)

@show @which sum(x)
@show @which diagm(x)
@show @which norm(x)
@show @which A * x
@show @which x + 1
@show @which x * 1

prints

@which(sum(x)) = sum(a) at reduce.jl:229
@which(diagm(x)) = diagm{T}(v::AbstractArray{T,1}) at linalg/dense.jl:121
@which(norm(x)) = norm(x::AbstractArray{T<:Any,1}) at linalg/generic.jl:201
@which(A * x) = *{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},S}(A::Union{Base.ReshapedArray{T,2,A<:DenseArray,MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N<:Any}}},DenseArray{T,2},SubArray{T,2,A<:Union{Base.ReshapedArray{T<:Any,N<:Any,A<:DenseArray,MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N<:Any}}},DenseArray},I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex,Colon,Int64,Range{Int64}},N<:Any}},L<:Any}}, x::Union{Base.ReshapedArray{S,1,A<:DenseArray,MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N<:Any}}},DenseArray{S,1},SubArray{S,1,A<:Union{Base.ReshapedArray{T<:Any,N<:Any,A<:DenseArray,MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N<:Any}}},DenseArray},I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex,Colon,Int64,Range{Int64}},N<:Any}},L<:Any}}) at linalg/matmul.jl:79
@which(x + 1) = +(A::AbstractArray, x::Number) at arraymath.jl:94
@which(x * 1) = *(A::AbstractArray, B::Number) at abstractarraymath.jl:97

while

m = Model()
x = @variable(m, [1: 3])
A = rand(3, 3)

@show @which sum(x)
@show @which diagm(x)
@show @which norm(x)
@show @which A * x
@show @which x + 1
@show @which x * 1

prints

@which(sum(x)) = sum(j::Array{JuMP.Variable,N<:Any}) at /Users/twan/code/julia/MixedIntegerExperiments/v0.5/JuMP/src/operators.jl:279
@which(diagm(x)) = diagm(x::Array{JuMP.Variable,1}) at /Users/twan/code/julia/MixedIntegerExperiments/v0.5/JuMP/src/operators.jl:351
@which(norm(x)) = norm{V<:JuMP.AbstractJuMPScalar}(x::Array{V,N<:Any}) at /Users/twan/code/julia/MixedIntegerExperiments/v0.5/JuMP/src/norms.jl:44
@which(A * x) = *{T<:Union{JuMP.AbstractJuMPScalar,JuMP.GenericNormExpr{2,Float64,JuMP.Variable},JuMP.GenericNorm{P,Float64,JuMP.Variable},JuMP.NonlinearExpression}}(A::Union{Array{T<:Any,2},SparseMatrixCSC}, x::Union{Array{T,1},Array{T,2},SparseMatrixCSC{T,Ti<:Integer}}) at /Users/twan/code/julia/MixedIntegerExperiments/v0.5/JuMP/src/operators.jl:447
@which(x + 1) = +{T<:Union{JuMP.AbstractJuMPScalar,JuMP.GenericNormExpr{2,Float64,JuMP.Variable},JuMP.GenericNorm{P,Float64,JuMP.Variable},JuMP.NonlinearExpression}}(lhs::Union{Array{T,N<:Any},SparseMatrixCSC{T,Ti<:Integer}}, rhs::Number) at /Users/twan/code/julia/MixedIntegerExperiments/v0.5/JuMP/src/operators.jl:507
@which(x * 1) = *{T<:Union{JuMP.AbstractJuMPScalar,JuMP.GenericNormExpr{2,Float64,JuMP.Variable},JuMP.GenericNorm{P,Float64,JuMP.Variable},JuMP.NonlinearExpression}}(lhs::Array{T,N<:Any}, rhs::Number) at /Users/twan/code/julia/MixedIntegerExperiments/v0.5/JuMP/src/operators.jl:538

(note: this PR in its current form doesn't address all of these).

Personally I don't see much value in putting AbstractArray over no type restriction at all if we're already going to loosen the type beyond what's actually used and has been tested.

The main value I see is that it conveys the intention of the code. Moreover, you still get to put constraints on the dimension and element type, and you need it to avoid ambiguities. However, widening signatures to ::AbstractArray also makes it so that if things don't work for a type that satisfies the requirements of being a (reasonable) AbstractArray subtype, then users get to complain, so I understand your hesitancy.

A (possibly overly) cheap response to your argument is that for some methods, you're not currently testing all possible values of N or T for Array{T, N} either. You couldn't possibly test all possible AbstractArray subtypes, but in my opinion, that doesn't mean that using ::Array instead of ::AbstractArray is necessarily the right thing to do.

[mlubin]

@joehuchette should comment but I don't think this method was designed for AbstractArrays. The built-in fallback should work

[joehuchette]

Yes, I don't think this part of the code will work with generic abstract arrays as-is.

[mlubin]

I'm concerned about the validity of using eachindex from ret to index into lhs and rhs when these could be different types.

[joehuchette]

I think eachindex(ret,lhs,rhs) should work?

[tkoolen]

Oops, forgot this one.

[tkoolen]

Fixed (here and in a number of other locations).

[mlubin]

Where is this used and why does it need to change?

[mlubin]

Are you redefining printing for AbstractArray{Variable}? If not, this doesn't need to be touched

[mlubin]

I could change the tests to use OffsetArrays, but I do foresee problems on 0.5, where support for such arrays is only partial at the moment.

It's the opposite. Any code in JuMP that's calling size on (now) an AbstractArray with custom indices is probably incorrect. This will throw an error on 0.5 which is fine. On 0.6 it may silently do the wrong thing which is not okay.

[mlubin]

One option is to follow the suggestion in the docs:

Note that if you don’t want to be bothered supporting arrays with non-1 indexing, you can add the following line:

@assert all(x->isa(x, Base.OneTo), indices(A))

at the top of any function.

[tkoolen]

I'll push some more changes tonight.

[tkoolen]

New round of changes, please take a look.

[mlubin]

Unless things have changed recently, splatting large arrays can cause big performance issues when the dimension is large, e.g., into the thousands. collect(j) might be better here

[tkoolen]

collect(::OffsetArray) returns a copy of the OffsetArray (with the same indices). I found it surprisingly hard to convert from an OffsetArray to a regular Array. Maybe I'm missing something.

[rdeits]

How about [v[i] for i in eachindex(v)] ?

[tkoolen]

D'oh. Yeah, that's better.

[mlubin]

Are these methods now covered by built-ins in Julia?

[tkoolen]

Yes. https://github.com/JuliaLang/julia/blob/bdef3b57d34ec90781d30f62adf73d048e4a3090/base/arraymath.jl#L36 on current master, https://github.com/JuliaLang/julia/blob/v0.5.0/base/arraymath.jl#L49 on 0.5.

[mlubin]

Ok, thanks for checking

[mlubin]

Getting closer. Please correct the travis failures as well. (Appveyor failures are expected.)

[tkoolen]

More fixes.

Two random comments on getting things to work for OffsetArrays, for future reference:

• the multi-argument eachindex is pretty useless when dealing with OffsetArrays. See Multi-arg eachindex results in error JuliaArrays/OffsetArrays.jl#24.
• even in 0.6, diagm doesn't work for OffsetArrays. See diagm doesn't work JuliaArrays/OffsetArrays.jl#25.

[joehuchette]

I'm somewhat uncertain if this method should be changed. Is there a compelling reason to pass a sparse matrix of linear constraints, for example?

[tkoolen]

Without changing this method, the following fails:

m = Model()
v = @variable(m, [1:3])
x = OffsetArray(v, -3)
@constraint(m, x .== 0)

because there is no Convert method that constructs an Array from an OffsetArray (a conscious design decision in OffsetArrays).

[mlubin]

This doesn't seem like the right fix. I think constructconstraint!(x::Array, sense::Symbol) = map(c->constructconstraint!(c,sense), x) should be changed to return a flat Vector of constraints. Then there's no need to touch addVectorizedConstraint.

[joehuchette]

Since this also doesn't work:

t = OffsetArray(rand(3), -3)
t .== 0

I'm not sure this is desirable behavior.

[tkoolen]

OK, will change. @joehuchette, that seems like a bug in Base.

[tkoolen]

Done.

[tkoolen]

Travis failures with no output seem to be because Amazon S3 is having issues: https://www.traviscistatus.com/incidents/hmwq9yy5dh9d.

[mlubin]

I'll restart the builds when travis stabilizes a bit.

[joehuchette]

As long as vec(::AbstractArray{Variable}) will reliably work, I think the definition above should work fine.

[tkoolen]

I tried that first, but AffExpr is a type alias of GenericAffExpr{Float64,Variable}, so that directly calls the inner constructor of AffExpr, which expects precisely a Vector{Variable}. I tried adding an inner constructor to GenericAffExpr that takes AbstractArrays (+ additional outer constructors), but ran into some more issues that I don't remember anymore. Basically, it was going to be a pretty big change, whereas sum seemed to be the only current case where handling AbstractVectors was needed when constructing AffExprs.

[joehuchette]

Cool, thanks for the explanation.

[joehuchette]

Modulo my iffy-ness about changing addVectorizedConstraint, this PR looks good to me.

[joehuchette]

the second ::Symbol shouldn't be necessary

[joehuchette]

Looks good to merge to me. @mlubin?

[mlubin]

I'm having trouble viewing the status of travis.

[mlubin]

@tkoolen, thanks for addressing all of our concerns. Looking forward to seeing the uses of AxisArrays.

[tkoolen]

Awesome, thanks guys!

[joehuchette]

Thanks for this, @tkoolen!