0.5 Adding matrices of special type incorrectly promoted to Matrix{Any}

[dlfivefifty]

abstract AbstractFoo
immutable Foo <: AbstractFoo end
import Base.+
+(::Foo,::Foo) = Foo()

A=fill(Foo(),10,10)
A+A

returns a Matrix{Any}, not Matrix{Foo} as expected.

Julia Version 0.5.0-dev+5122
Commit 2e1bcc6 (2016-07-03 00:01 UTC)
Platform Info:
  System: Darwin (x86_64-apple-darwin15.5.0)
  CPU: Intel(R) Core(TM) i7-3820QM CPU @ 2.70GHz
  WORD_SIZE: 64
  BLAS: libgfortblas
  LAPACK: liblapack
  LIBM: libopenlibm
  LLVM: libLLVM-3.7.1 (ORCJIT, ivybridge)

[pkofod]

I had a version from 1st of July, and that didn't have that error, so it must have come in the last two days,

[martinholters]

As part of #16995, the output type is not longer guessed to be equal to the input type, so at the moment, you need to provide your own promote_op method to make this work.

[pkofod]

Should maybe have had a mention of this in NEWS? @nalimilan

[nalimilan]

As part of #16995, the output type is not longer guessed to be equal to the input type, so at the moment, you need to provide your own promote_op method to make this work.

I don't think that's the plan. The idea is that we should use the same algorithm as in broadcast to automatically compute the return type. broadcast(+, A, A) correctly returns an Array{Foo} already.

@JeffBezanson Is that correct?

[nalimilan]

In the meantime, note that if you declare Foo <: Number and provide a one(::Type{Foo}) method, rules are already present to choose the correct return type for all operators.

[wbhart]

Can you provide a link to the piece of code where these rules are defined for numbers? We can't make our types belong to Number, since the top of our abstract type hierarchy is Set, and obviously sets are not always sets of numbers.

By the way, one(::Type{Foo}) is really not the right thing to ask for. You should ask for one(::Foo), since the former will not work for types like the integers modulo n, since it is not possible to create 1 mod n from its type alone, since you need to know what n is. That's only possible if you have another value modulo n to play with, since you can look inside it to see what the n is supposed to be. A lot of the Julia matrix routines fail for us because of this.

[martinholters]

Can you provide a link to the piece of code where these rules are defined for numbers?

https://github.com/JuliaLang/julia/blob/master/base/number.jl#L72

By the way, one(::Type{Foo}) is really not the right thing to ask for. You should ask for one(::Foo)

But promote_op is only given the type, so this would require a major overhaul.

[nalimilan]

It's here:

julia/base/number.jl

Line 73 in 6ab1b19

T = typeof(op(one(R), one(S)))

That rule was already present in 0.4, what's new is that we no longer assume that the return type of an operation is the same as it's input type (which isn't correct in general).

I don't understand what you mean by "the top of our abstract type hierarchy is Set, and obviously sets are not always sets of numbers." Is Foo an AbstractSet? Or do you want to store Foo inside Set? Anyway, you can simply override Base.promote_op for now if you need to support the 0.5 development version immediately.

Regarding integers modulo n, why don't you make n a type parameter? Or can't you return a dummy value for n (like 0 or 1) which would be dropped in favor of the modulo used by the other argument in operations?

[wbhart]

Thanks for the link.

Just to clarify, we hit this issue in our Nemo package, unrelated to the package maintained by the OP. It currently prevents our documentation from building, since documentation is built on the nightly build. Which is a shame, since I am currently rewriting our documentation.

What I mean is that we have an abstract type called Set, which sits at the top of our type hierarchy, i.e. FlintIntegerRing <: Nemo.Ring <: Nemo.Group <: Nemo.Set. (Our type FlintIntegerRing is the equivalent of Foo in the OP's example). So we can't just make our types belong to Number, since that doesn't make any sense.

Regarding integers modulo n, you can't make n a type parameter. Among the reasons:

1. n is a bignum, not an Int
2. multimodular algorithms can use thousands or millions of different primes n in a calculation. If n was a type parameter, the Julia compiler would recompile every function used for every different value of n, which would take months for a computation that is supposed to take seconds.
3. Many Julia developers have told us not to do this.

You can't put a dummy value for n, obviously. What happens if you get two of them? What should the operation do then? And when will the correct information get filled in? How would we ensure the integrity of our data? Or what happens if the algorithm needs to compute -a mod n. How will it compute the correct residue -a if it doesn't know n? E.g. if I give you s = 3 mod something and ask for -s, what is the value that should be stored for the residue?

It's clear upon reflection that there are only two options for implementing integers mod n. 1) put n in the type, 2) put n in the actual values. Since the former doesn't work, for the reasons outlined above, there is only one workable way of implementing integers mod n. And that has implications if you are implementing say a matrix determinant algorithm (which should work just fine over Z/nZ, so long as it doesn't do any divisions, e.g. using clow sequences or the charpoly method, etc).

There really is only one correct way of doing this sort of thing. Someone already noticed this, and Julia defines things like zero(1) so that this is possible. But the matrix code in Julia doesn't use it systematically, it seems. We are pretty much unable to use any of the sophisticated matrix algorithm implementations in Julia for this reason (determinant, multiplication, you name it). Anyway, it's off topic here. I am just mentioning it to raise awareness every time I see it.

[nalimilan]

What I mean is that we have an abstract type called Set, which sits at the top of our type hierarchy, i.e. FlintIntegerRing <: Nemo.Ring <: Nemo.Group <: Nemo.Set. (Our type FlintIntegerRing is the equivalent of Foo in the OP's example). So we can't just make our types belong to Number, since that doesn't make any sense.

OK. Since multiple inheritance isn't supported yet, the best solution is to declare promote_op methods for each supported operation. But AFAICT this should only be temporary, so you may want to wait for the final 0.5 release to avoid wasting your time.

I can't really comment on the best implementation of integers modulo n, but it could be problematic for a Number to be unable to zero(::Type) and one(::Type). I guess you should try to get rid of their uses in Base and see whether it's workable (promote_op isn't a totally legitimate use of these IMHO, so that's not the main problem).

[wbhart]

[OT] I don't see any problem with Number requiring zero(::Type) and one(::Type). It just means integers mod n are not Numbers. The problem of course is that you want to put things that are not numbers in matrices. Numbers are only one very tiny part of algebra. You can define matrices over any ring (at least). The standard example I use to illustrate why the current Julia model breaks is the integers mod n (which is a ring, even if it is not a ring of numbers). It's not a perfect example, since you could use n = 0 as you suggest and then write a pile of logic to handle all the cases where n = 0. But there are other more complex examples in algebra where there is no canonical value you can use to indicate missing information, for example quotient rings, i.e. a ring R modulo an ideal n. Here every ideal is a valid value for n, so there is no equivalent of n = 0 that you could use to indicate missing information. The problem becomes intractable if you are trying to write a generic system for abstract algebra (as we are).

[andreasnoack]

@wbhart Wouldn't the best solution be to define promote_op in for the Nemo types?

[wbhart]

@andreasnoack We'll wait to see how the issue gets resolved. It looks to me like people are suggesting this won't be required for the final 0.5 release.

[wbhart]

@martinholters Yes, this is another example of the same problem I mention in my "off topic" bit. We already define our own methods for promoting the "type" of abstract algebraic values when doing arithmetic operations in Julia, since promote_rule and friends are only given the type, not the value, so that doesn't work for us. The type based system clearly can't work in general. It's fine for numerical mathematics, but no use for any other kind of maths.

It shouldn't cause any performance issues to pass the values instead of the types to functions like this. For all the numerical types present in Julia it would then just default to calling the current type based functions, and the compiler should just optimise this out to exactly the same assembly instructions as at present. But at least people not doing numerical work could overload the value based methods directly.

[wbhart]

@martinholters Sorry my previous post was confused. Let me be more specific.

If I'm interpreting the documentation correctly, promote_op is a replacement for promote_rule on a per operation basis. I assume this is correct?

In fact, the problem for us is probably not with promote_rule and friends. The problem is actually just that we can't rely on promote_rules/promote_ops and convert functions to do our promotions for us when doing arithmetic operations.

The reason is that convert does not have enough information to actually do a conversion. For example, in adding an integer mod n to an integer, convert would only have the type of the first argument and the value of the second, whereas it really needs the value of both in order to do the conversion (since the modulus n is in the value, not the type of the first argument).

The only way around it is for us to define promote(xs..) directly for every possible combination of types that we define, as the documentation recommends we should not do [1].

So the problem actually isn't with promote_rule and promote_op as such, as far as I can see, but with convert. And in any case, as the type of the result should not depend on which operation has been performed, in our case, we would be using promote_rule, not promote_op, I think (though there are doubtlessly cases where we will need the latter, eventually).

What I am still confused about is this. We don't use promote_op, and we have promote_rules for all our promotions. Why when adding two matrices together where all entries in both matrices are of the same type, is it not using the promote_rule we defined to determine the type? Why would it even check the promote_op, since we never defined one. Something seems definitely broken if I am now understanding how this is supposed to work.

[1] http://docs.julialang.org/en/latest/manual/conversion-and-promotion/