Enable operator-sensitive extension of element-type promotion

[timholy]

This fixes the long-standing #8027, providing much better generic fallbacks for elementwise operations involving arrays and array/scalar operations. The proximal motivation for tackling this was to fix the ambiguity warnings triggered by #12115, but this will also (given suitable type-stability improvements to SIUnits) finally allow eltype([1Meter]+1Meter) == typeof(1Meter) and eltype([1Meter]*(1Meter)) == typeof(1Meter^2). In some ways this turns out to be the (much easier) companion of #10525, defining operations that involve arrays in terms of operations defined for scalars. I think it will allow cleanup of a certain amount of torturous code.

I had forgotten why we had promote_array_type; to save any reviewers the trouble, it turns out to be because we've decided that eltype([1.0f0]+1.0) == Float32 (i.e., the eltype of the Array wins, when working with Real numbers). This rule obviously doesn't work for more generic types, but the exception does make sense so I've kept it.

If this meets with approval, I'll see about applying it to #12115.

[timholy]

I'd love a few more words of explanation. (Besides, I enjoy danger.)

[JeffBezanson]

Yes, with apologies, we cannot do this sort of thing. Program behavior cannot depend on return_types. There is no guarantee as to what that function will return.

[jakebolewski]

I know, it was a tongue-in-cheek response :-)

We already have return_types sprinkled in a couple places in Base, @carnaval and I had fun uncovering their existence.

[JeffBezanson]

Oh wow, I see there is one in subarray. Looks like that is the only one. We need to remove that too. Don't know how I missed that.

[jakebolewski]

Well, it leads to interesting debugging sessions.

[carnaval]

heh, I had a bad feeling reading the title of the PR in my mail client. Unfortunately I feel that this restriction will only seem more arbitrary as inference gets better, but we have to be firm about that...

[carnaval]

and yeah, this one in a stagedfunction is incorrect if inference returns a non leaf in that case

[timholy]

When inference breaks, in my experience it returns something wider than optimal. But at least that's still correct. In contrast, here's current master (a stripped-down version of SIUnits that hasn't gone out of its way to be type-unstable as a way of avoiding this bug):

immutable MyType{T,N} <: Number
    val::T
end
*{T}(x::MyType{T,1}, y::MyType{T,1}) = MyType{T,2}(x.val*y.val)
a = MyType{Int,1}[MyType{Int,1}(1)]

julia> a*MyType{Int,1}(1)
ERROR: MethodError: `convert` has no method matching convert(::Type{MyType{Int64,1}}, ::MyType{Int64,2})
This may have arisen from a call to the constructor MyType{Int64,1}(...),
since type constructors fall back to convert methods.
Closest candidates are:
  MyType{T,N}(::Any)
  call{T}(::Type{T}, ::Any)
  convert{T<:Number}(::Type{T<:Number}, ::Base.Dates.Period)
  ...
 in .* at arraymath.jl:110
 in * at abstractarraymath.jl:55

In contrast, with this PR:

julia> a*MyType{Int,1}(1)
1-element Array{MyType{Int64,2},1}:
 MyType{Int64,2}(1)

which is exactly what you want.

To make life more interesting, let's try to break it following @carnaval's concern. Master:

immutable MyType{T,N} <: Number
    val::T
end
# This is deliberately type-unstable
*{S,T}(x::MyType{S,1}, y::MyType{T,1}) = x.val > y.val ? MyType{T,2}(x.val*y.val) : MyType{S,2}(x.val*y.val)
a = MyType{Int,1}[MyType{Int,1}(1)]

julia> a*MyType{Float32,1}(1.0f0)
ERROR: no promotion exists for MyType{Float32,1} and MyType{Int64,1}
 in .* at arraymath.jl:108
 in * at abstractarraymath.jl:55

This PR:

julia> a*MyType{Float32,1}(1.0f0)
1-element Array{Union{MyType{Float32,2},MyType{Int64,2}},1}:
 MyType{Int64,2}(1)

There's absolutely nothing wrong with that, either. (And it's not a leaf type.)

One giant pink inflatable mustache to the first person who gives a test case where this PR fails but master works. (Include shipping address and allow 6-8 weeks for delivery.)

[malmaud]

Doesn't the behavior of cfunction also depend on return type inference?

[carnaval]

@timholy the one that was broken was in the subarray code (but I remember you saw it on the typeinf-free comphrension branch), where if return_types() was too wide, the subarray construction would fail because some index type that was spliced directly into the ast would not be <:Union{Int,Range{Int}} or something like that (from memory).

But, in more generality, the problem with this kind of code is that it makes the fact that inference is enabled or not user visible. It is already the case today (typegoto, maybe cfunction ? although I'm not familiar with this one but it looks like you have to specify a return type explicitely here), but we should be removing such things not adding some.

So, a concrete problem that could arise for example, is that the behavior of my program (e.g. my promoted array's type) could change at the whim of one of the various MAX_XXX constants in inference (or the logic in type_too_complex, etc, ...).

[StefanKarpinski]

Doesn't the behavior of cfunction also depend on return type inference?

Only in the sense that whether a particular cfunction invocation works or not depends on type inference. If it works in two different versions, then it does the same thing. That's a weaker form of dependency but still not awesome. This could be fixed by automatically adding convert calls and/or type assertions, which is what people end up doing manually to use cfunction as it stands.

[timholy]

The bug you refer to was unfortunate---very sorry it caused you trouble. I'm completely fine with giving up that return_type call in the subarray code, it's nowhere near as important as this case is (and was probably just me being anal about types). Here, I don't think the same problem applies, since an Array can be created with any T. Ergo the default failure scenario is that you allocate an array of wider type than strictly necessary, but that's better than allocating one that can't contain the elements.

Aside from the fact that it usually works, the other good thing about this design is that in cases of trouble, just overload promote_op for the specific failing case(s). Currently, we're doomed if a+b and/or a*b give a type that can't be converted to that of promote(a,b). As a solution, promote_type is simply broken because it lacks information about the operator; in my view, promote_type is only suitable for containment (e.g., concatenation). You can think of this PR as an implementation of @typeof, with the added bonus that it can be overridden as needed.

I'd argue that since we're currently in bad shape every place this comes up, and this makes things better, we should just fess up to the fact that we can't currently solve this problem without leveraging inference. We can aim for a more long-term fix in a future release.

[JeffBezanson]

Yes, this "weak dependence" on type inference is not as bad. If you check the equivalent of return_types for a concrete type and give an error otherwise you're sort of locally enabling static type checking. Of course, having code stop working when type inference is disabled is not ideal.

Stefan is right about cfunction; it's not a problem because in that case the code always says exactly what type is required. It's just a convenience for the inferred type to match that type.

But there is a real problem to fix here. The leading approach is typeof(f(one(T),one(S))), but that just pushes the problem to one or zero. The best solution I can think of is to remove most vectorized functions, and have the remaining ones call map.

[JeffBezanson]

Also the promote_rule for SIUnits is arguably broken. Promoting two numbers is supposed to return two values of the same type. This is only meant to assist dispatch by mapping O(n^2) argument combinations to O(n). I agree that using promotion to determine the result of array + is broken; all it can do is convert values to a common type, and nothing more.

[carnaval]

@timholy but the problem IMO is that the failure mode of a wider array is visible behavior. Imagine we implement a heuristic that makes inference try harder on more powerful machines ? Or only when you precompile ? Or only when you pass -O10 command line ? We can't have types in the program change because of that since the user might dispatch on it or observe it in any way and therefore the "result" of the program is not uniformely defined anymore.

I agree that those are not very practical concerns, but I feel this is an area we have to be more pedantic about since in the future we might have cases where we wouldn't be able to change the behavior of inference without breakage.

(and no worries about the subarray bug obviously ;-))

[timholy]

Yes, many of the promotion rules in SIUnits pay loving attention to the type of the numeric coefficient, but completely punt on the units. That's quite ironic for a package whose raison d'etre is handling units. But I think defaulting to the generic types is basically defensive behavior necessitated by the bug that this PR fixes.

[timholy]

@carnaval, I'd lobby for two inference engines, one stripped-down version specifically for handling this problem and the more powerful one for general use.

[malmaud]

What about using return_types, but raising an error/follow a different codepath it it doesn't return a single concrete type? That would make the situation in some ways more analogous to the current cfunction situation.

[timholy]

@malmaud, this PR basically does that, but allows a non-concrete type (which is sometimes needed, see example above). We could explore having it return Any instead of an error, but I worry about that in the context of #265.

[carnaval]

@timholy why not, but the stripped down version would have to have the characteristics of a static type system inference, that is, formal deduction rules for checking that only depends on syntax & an inference engine that bails out statically if the syntax is not explicit enough forcing you to add explicit type annotations

[timholy]

@carnaval, that doesn't sound too bad to me. Especially since I am skeptical that there is another practical solution.

[JeffBezanson]

The other practical solution is to do what map does.

[timholy]

Yeah, that would work too (EDIT: except for empty arrays, as I seemed to remember above but forgot here untill @mbauman reminded me below), but it does stink from a performance perspective. I imagine that operations like z+0.5 (where length(z)==3) would suffer quite a lot.

[carnaval]

there is a way to write that so that if inference figures it out then it elides the slow path just as map does right ?

[timholy]

I'm not convinced there is (currently), but I'd love to be proven wrong:

julia> function foo{T}(x::Vector{T}, y::Number)
           if isleaftype(T)
               y
           else
               error("oops")
           end
       end
foo (generic function with 1 method)

julia> foo(rand(2), 3.5)
3.5

julia> @code_llvm foo(rand(2), 3.5)

define double @julia_foo_21033(%jl_value_t*, double) {
top:
  %2 = alloca [5 x %jl_value_t*], align 8
  %.sub = getelementptr inbounds [5 x %jl_value_t*]* %2, i64 0, i64 0
  %3 = getelementptr [5 x %jl_value_t*]* %2, i64 0, i64 2
  store %jl_value_t* inttoptr (i64 6 to %jl_value_t*), %jl_value_t** %.sub, align 8
  %4 = getelementptr [5 x %jl_value_t*]* %2, i64 0, i64 1
  %5 = load %jl_value_t*** @jl_pgcstack, align 8
  %.c = bitcast %jl_value_t** %5 to %jl_value_t*
  store %jl_value_t* %.c, %jl_value_t** %4, align 8
  store %jl_value_t** %.sub, %jl_value_t*** @jl_pgcstack, align 8
  store %jl_value_t* null, %jl_value_t** %3, align 8
  %6 = getelementptr [5 x %jl_value_t*]* %2, i64 0, i64 3
  store %jl_value_t* null, %jl_value_t** %6, align 8
  %7 = getelementptr [5 x %jl_value_t*]* %2, i64 0, i64 4
  store %jl_value_t* null, %jl_value_t** %7, align 8
  %8 = load %jl_value_t** %4, align 8
  %9 = getelementptr inbounds %jl_value_t* %8, i64 0, i32 0
  store %jl_value_t** %9, %jl_value_t*** @jl_pgcstack, align 8
  ret double %1
}

I haven't looked at the code produced by map, but I'd be surprised (& pleased!) if it's any better.

[carnaval]

I don't think there is any reason for all this gcframe traffic, this should end up as ret %1

[timholy]

dest = zeros(10^6)
src = rand(10^6)

julia> @time Base.promote_to!(Base.IdFun(), 1, dest, src);
   7.229 milliseconds (4 allocations: 144 bytes)

julia> @time copy!(dest, src);
   2.877 milliseconds (4 allocations: 144 bytes)

Not as bad as I'd feared, but it's still more than two-fold.

[yuyichao]

Somehow type inference doesn't know anything about isleaftype.

[mbauman]

Oh man, I'd love to see this solved.

I think that the major downside to using map's solution is that it extends the empty map problem to arithmetic (what is Int[] + 1.0?)… which (at least currently) would make arithmetic type-unstable.

[carnaval]

@yuyichao even if we would get tighter type info if inference special cased leaftype for branches, codegen should still not generate an obviously useless gcframe here

[timholy]

It's a nice scheme. I suspect that the main limitation would be handling types that take multiple parameters (more than just degree). Unless you have a plan that I missed?

Long-term, I suspect the best way to handle this would be a separate ReturnType package that solves the problem by inference (not using any of julia's native inference code). But I don't have any immediate plans to write such a package.

[davidagold]

The parameter N in the type Polytypic{N} denotes the location of the degree parameter in a type that has multiple parameters. So, for instance, you could do

immutable OtherMeter{T, N} end
typicity{M1<:OtherMeter, M2<:OtherMeter}(::Base.MulFun, ::Type{M1}, ::Type{M2}) = Polytypic{2}()

to obtain the desired behavior:

julia> promote_op(*, OtherMeter{Int, 1}, OtherMeter{Int, -1})
Unity

An alternative approach may become necessary if one is intent on treating argument signatures in which the location of the relevant parameter (be it degree or something else) differs amongst arguments. In that case, probably one would have to define an internal degree_ind method that takes a type argument and gives the location of the degree parameter:

typicity{M1<:Meter, M2<:OtherMeter}(::Base.MulFun, ::Type{M1}, ::Type{M2}) = Polytypic{0}()
degree_ind{T<:Meter}(::Type{T}) = 1
degree_ind{T<:OtherMeter}(::Type{T}) = 2

function promote_type_radical{R<:Meter, S<:OtherMeter}(::Type{R}, ::Type{S})
    T = S.parameters[1]
    return OtherMeter{T}
end

function promote_typicity{R, S}(::Polytypic{0}, ::Type{R}, ::Type{S})
    degree = R.parameters[degree_ind(R)] + S.parameters[degree_ind(S)]
    degree == 0 ? Unity : promote_type_radical(R, S){degree}
end

julia> promote_op(*, Meter{1}, OtherMeter{Int, 1})
OtherMeter{Int64,2}

Note that I've replaced the use of promote_type with a new method promote_type_radical. The impetus for this is the problem that derives from relying on promote_type as in my previous post in cases such as promote_op(*, Meter{1}, Meter{1}). (Though that's not a problem here.)

I suspect that most cases would be able to get by fine using the first approach, i.e. using the parameter in Polytypic{N} to indicate the location of the degree parameter. But the latter approach demonstrates that it's possible to do more complex things, if one is willing to extend the interface a bit.

[timholy]

But suppose you have a type with 3 parameters and the degree parameter is the second one. In your scheme, how do the constructors work out?

My overall view is that this might be appropriate for a package but probably is not sufficiently general for Base. One might make this work for SIUnits-like types, but 20a0b84 doesn't seem to be implementable by this scheme.

[timholy]

Barring a "give me a little more time to review" request or further discussion of the core elements of this PR, I've inserted a merge reminder into by calendar for first thing Wednesday morning. That gives folks another 48 hours.

[tkelman]

Looks like Images (and maybe other packages?) is using promote_array_type, but you're the right person to deal with that.

[davidagold]

@timholy that's a good point. I do think it's possible to make it work -- one would need a complement to promote_type_radical to fill in the rest of the parameters. But that is an even further extension of the interface.

I wouldn't expect to cover all patterns with just the Homotypic and Polytypic schemes. One would require Typicitys to cover other patterns, like

promote_typicity{R, S}(::Monotypic, ::Type{R}, ::Type{S}) = R

to replace something like

Base.promote_op{P<:Period,R<:Real}(::$F, ::Type{P}, ::Type{R}) = P

with

typicity{P<:Period, R<:Real}(::$F, ::Type{P}, Type{R}) = Monotypic()

Of course, some cases wouldn't fall into a well-defined, pervasive return type promotion pattern. In these cases, it's probably easiest just to override Base.promote_op directly. In general, I think you're right that I'd need to work on the abstraction before I'd be able to demonstrate suitable generality.

[lindahua]

Being able to infer the type of some arithmetic operations is useful.

We already have functors here: https://github.com/JuliaLang/julia/blob/master/base/functors.jl

What about we just add result_type methods to those functors (instead of introducing a whole series of similar things like GenericNFunc)?

result_type{X,Y}(::AddFun, ::Type{X}, ::Type{Y}) = #...

Also, the NumFunctors package (https://github.com/lindahua/NumericFuns.jl) already provides similar things. The development there was stopped a while ago. If necessary, I can clean-up that package and make it work with latest Julia.

[timholy]

Rebased on static functors (and on the latest master). Not a @generated function in sight. I added a test.

Since there's no substantive change, I'll merge later today once CI passes.

[timholy]

Hmm, I haven't seen this travis failure before:

$ cp /tmp/julia/lib/julia/sys.ji local.ji && /tmp/julia/bin/julia -J local.ji -e 'true' && /tmp/julia/bin/julia-debug -J local.ji -e 'true' && rm local.ji
cp: cannot stat `/tmp/julia/lib/julia/sys.ji': No such file or directory
The command "cp /tmp/julia/lib/julia/sys.ji local.ji && /tmp/julia/bin/julia -J local.ji -e 'true' && /tmp/julia/bin/julia-debug -J local.ji -e 'true' && rm local.ji" exited with 1.
0.01s$ /tmp/julia/bin/julia -e 'versioninfo()'
/home/travis/build.sh: line 41: /tmp/julia/bin/julia: No such file or directory

Oh, how I dream of world where CI failures make sense.

I'll wait a few minutes before restarting it, in case anyone wants to poke at it.

[tkelman]

Scroll up. Download failure, I think. The caching server hasn't been doing all that great a job lately.

  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
  0     0    0     0    0     0      0      0 --:--:--  0:00:15 --:--:--     0
curl: (28) Operation too slow. Less than 1 bytes/sec transferred the last 15 seconds
  0     0    0     0    0     0      0      0 --:--:--  0:00:15 --:--:--     0
curl: (28) Operation too slow. Less than 1 bytes/sec transferred the last 15 seconds

[mbauman]

More specifically:

curl: (22) The requested URL returned error: 403
make[2]: *** [libgit2-159061a8ce206b694448313a84387600408f6029.tar.gz] Error 22
...
curl: (22) The requested URL returned error: 403
make[2]: *** [utf8proc-85789180158ac7fff85b9f008828d6ac44f072ea.tar.gz] Error 22

Errors from parallel make are always terrible. I just search for the word "error" and try to find the first meaningful one.

[tkelman]

Aha, good eyes Matt. We could turn off parallel make on Travis if it would help make the error messages any more sensible, but the queue has actually been worse on Travis than on AppVeyor ever since the win64 freeze bug was resolved. We could also turn OSX Travis back off, it might help things along a little faster with the queue if it's worth potentially not automatically catching platform differences.

[timholy]

Thanks for the tips, folks.

At the current moment, I'm not sure Traivs-OSX is notably slower than Travis-Linux.

[yuyichao]

@tkelman How hard it is to implement on travis the logic to skip the build if there's a new commit on the same PR like what we have on AppVeyor

[tkelman]

@yuyichao not very, let me find where I last put that code. I think I did something like that for libgit2 and it was on their master for a few days but then they reverted it.