Skip to content

Commit

Permalink
Introduce promotion arithmetic
Browse files Browse the repository at this point in the history
  • Loading branch information
@andreasnoack
andreasnoack committed Mar 12, 2014
1 parent bc1ced6 commit 1257643
Show file tree
Hide file tree
Showing 3 changed files with 81 additions and 7 deletions.
80 changes: 76 additions & 4 deletions base/linalg/matmul.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -266,11 +266,83 @@ end
# strides != 1 cases in libalg_blas.jl
(*){T,S}(A::StridedMatrix{T}, B::StridedVector{S}) = generic_matvecmul('N', A, B)

arithtype(T) = T
arithtype(::Type{Bool}) = Int
*(::Type{Bool},::Type{Bool}) = Bool
*(::Type{Bool},::Type{Int}) = Int
*(::Type{Bool},::Type{Float32}) = Float32
*(::Type{Bool},::Type{Float64}) = Float64
*(::Type{Bool},::Type{Complex{Float32}}) = Complex{Float32}
*(::Type{Bool},::Type{Complex{Float64}}) = Complex{Float64}
*(::Type{Int},::Type{Bool}) = Int
*(::Type{Int},::Type{Int}) = Int
*(::Type{Int},::Type{Float32}) = Float32
*(::Type{Int},::Type{Float64}) = Float64
*(::Type{Int},::Type{Complex{Float32}}) = Complex{Float32}
*(::Type{Int},::Type{Complex{Float64}}) = Complex{Float64}
*(::Type{Float32},::Type{Bool}) = Float32
*(::Type{Float32},::Type{Int}) = Float32
*(::Type{Float32},::Type{Float32}) = Float32
*(::Type{Float32},::Type{Float64}) = Float64
*(::Type{Float32},::Type{Complex{Float32}}) = Complex{Float32}
*(::Type{Float32},::Type{Complex{Float64}}) = Complex{Float64}
*(::Type{Float64},::Type{Bool}) = Float64
*(::Type{Float64},::Type{Int}) = Float64
*(::Type{Float64},::Type{Float32}) = Float64
*(::Type{Float64},::Type{Float64}) = Float64
*(::Type{Float64},::Type{Complex{Float32}}) = Complex{Float64}
*(::Type{Float64},::Type{Complex{Float64}}) = Complex{Float64}
*(::Type{Complex{Float32}},::Type{Bool}) = Complex{Float32}
*(::Type{Complex{Float32}},::Type{Int}) = Complex{Float32}
*(::Type{Complex{Float32}},::Type{Float32}) = Complex{Float32}
*(::Type{Complex{Float32}},::Type{Float64}) = Complex{Float64}
*(::Type{Complex{Float32}},::Type{Complex{Float32}}) = Complex{Float32}
*(::Type{Complex{Float32}},::Type{Complex{Float64}}) = Complex{Float64}
*(::Type{Complex{Float64}},::Type{Bool}) = Complex{Float64}
*(::Type{Complex{Float64}},::Type{Int}) = Complex{Float64}
*(::Type{Complex{Float64}},::Type{Float32}) = Complex{Float64}
*(::Type{Complex{Float64}},::Type{Float64}) = Complex{Float64}
*(::Type{Complex{Float64}},::Type{Complex{Float32}}) = Complex{Float64}
*(::Type{Complex{Float64}},::Type{Complex{Float64}}) = Complex{Float64}
*{T<:Number,S,N}(::Type{T},::Type{Array{S,N}}) = Array{T*S,N}

(+)(::Type{Bool},::Type{Bool}) = Int
(+)(::Type{Bool},::Type{Int}) = Int
(+)(::Type{Bool},::Type{Float32}) = Float32
(+)(::Type{Bool},::Type{Float64}) = Float64
(+)(::Type{Bool},::Type{Complex{Float32}}) = Complex{Float32}
(+)(::Type{Bool},::Type{Complex{Float64}}) = Complex{Float64}
(+)(::Type{Int},::Type{Bool}) = Int
(+)(::Type{Int},::Type{Int}) = Int
(+)(::Type{Int},::Type{Float32}) = Float32
(+)(::Type{Int},::Type{Float64}) = Float64
(+)(::Type{Int},::Type{Complex{Float32}}) = Complex{Float32}
(+)(::Type{Int},::Type{Complex{Float64}}) = Complex{Float64}
(+)(::Type{Float32},::Type{Bool}) = Float32
(+)(::Type{Float32},::Type{Int}) = Float32
(+)(::Type{Float32},::Type{Float32}) = Float32
(+)(::Type{Float32},::Type{Float64}) = Float64
(+)(::Type{Float32},::Type{Complex{Float32}}) = Complex{Float32}
(+)(::Type{Float32},::Type{Complex{Float64}}) = Complex{Float64}
(+)(::Type{Float64},::Type{Bool}) = Float64
(+)(::Type{Float64},::Type{Int}) = Float64
(+)(::Type{Float64},::Type{Float32}) = Float64
(+)(::Type{Float64},::Type{Float64}) = Float64
(+)(::Type{Float64},::Type{Complex{Float32}}) = Complex{Float64}
(+)(::Type{Float64},::Type{Complex{Float64}}) = Complex{Float64}
(+)(::Type{Complex{Float32}},::Type{Bool}) = Complex{Float32}
(+)(::Type{Complex{Float32}},::Type{Int}) = Complex{Float32}
(+)(::Type{Complex{Float32}},::Type{Float32}) = Complex{Float32}
(+)(::Type{Complex{Float32}},::Type{Float64}) = Complex{Float64}
(+)(::Type{Complex{Float32}},::Type{Complex{Float32}}) = Complex{Float32}
(+)(::Type{Complex{Float32}},::Type{Complex{Float64}}) = Complex{Float64}
(+)(::Type{Complex{Float64}},::Type{Bool}) = Complex{Float64}
(+)(::Type{Complex{Float64}},::Type{Int}) = Complex{Float64}
(+)(::Type{Complex{Float64}},::Type{Float32}) = Complex{Float64}
(+)(::Type{Complex{Float64}},::Type{Float64}) = Complex{Float64}
(+)(::Type{Complex{Float64}},::Type{Complex{Float32}}) = Complex{Float64}
(+)(::Type{Complex{Float64}},::Type{Complex{Float64}}) = Complex{Float64}

function generic_matvecmul{T,S}(tA::Char, A::StridedMatrix{T}, B::StridedVector{S})
C = Array(promote_type(arithtype(T),arithtype(S)), size(A, tA=='N' ? 1 : 2))
C = Array(T*S + T*S, size(A, tA=='N' ? 1 : 2))
generic_matvecmul(C, tA, A, B)
end

Expand Down Expand Up @@ -323,7 +395,7 @@ end
function generic_matmatmul{T,S}(tA, tB, A::StridedVecOrMat{T}, B::StridedMatrix{S})
mA, nA = lapack_size(tA, A)
mB, nB = lapack_size(tB, B)
C = Array(promote_type(arithtype(T),arithtype(S)), mA, nB)
C = Array(T*S + T*S, mA, nB)
generic_matmatmul(C, tA, tB, A, B)
end

Expand Down
6 changes: 4 additions & 2 deletions base/operators.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,8 +143,10 @@ At_ldiv_Bt(a,b) = transpose(a)\transpose(b)
oftype{T}(::Type{T},c) = convert(T,c)
oftype{T}(x::T,c) = convert(T,c)

zero(x) = oftype(x,0)
one(x) = oftype(x,1)
zero(x::Number) = oftype(x,0)
zero{T<:Number}(::Type{T}) = oftype(T,0)
one(x::Number) = oftype(x,1)
one{T<:Number}(::Type{T}) = oftype(T,1)

sizeof(T::Type) = error(string("size of type ",T," unknown"))
sizeof(T::DataType) = if isleaftype(T) T.size else error("type does not have a native size") end
Expand Down
2 changes: 1 addition & 1 deletion test/numbers.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1575,7 +1575,7 @@ approx_eq(a, b) = approx_eq(a, b, 1e-6)
for T in [Int,BigInt], n = [1:1000,1000000]
n = convert(T,n)
f = factor(n)
@test n == prod([p^k for (p,k)=f])
@test n == prod(T[p^k for (p,k)=f])
prime = n!=1 && length(f)==1 && get(f,n,0)==1
@test isprime(n) == prime

Expand Down

21 comments on commit 1257643

@JeffBezanson
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure this will fly.

@StefanKarpinski
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed, this doesn't seem nearly generic enough. We definitely need to do something but I don't think this is it.

@andreasnoack
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wasn't aware that others were watching. I am just playing with it. The motivation was the work on the generic linear algebra, but also the promotion rules for multiplication that have been discussed today. Solving linear systems with LU is stable under division, with QR it is stable under normalisation with 2-norm and for multiplication it is T*S+T*S stability and therefore I wanted to try to define promotion with respect to some arithmetic operation. Does it make sense?

@jiahao
Copy link
Member

@jiahao jiahao commented on 1257643 Mar 12, 2014

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's a very cogent summary of the type stability mess. Ugh.

@JeffBezanson
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think there are 2 big issues:

  1. Some types that you think of as groups are not really groups, since T+T does not give T
  2. The type of a matrix is not sufficient to specify what ring you're in, since it doesn't include the size.

@ivarne
Copy link
Member

@ivarne ivarne commented on 1257643 Mar 12, 2014

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Jeff and Stefan gets a email any time someone pushes to julia.

@lindahua
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I found that some machinery that can tell you the result type of a math operation is very useful. I actually ended up writing some similar stuff in developing a package: https://github.com/lindahua/Accelereval.jl/blob/master/src/maptype.jl

This file answers the questions like: given two arrays with element type T1 and T2, what is the element type of the result array. Also note that the type behaviors are different for scalar computation and vectorized computation.

I have always been looking forward to Julia providing something with which we may get the result type of an operation. However, I hope that it is something more generic than this.

@andreasnoack
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JeffBezanson Regarding 1, I'd expect that it was covered since e.g. (+)(Bool,Bool)=Int, but maybe I'm not getting the point. My imagination and judgement is limited by my lack of knowledge about deeper levels of the type inference. My thought was that basically my "type arithmetic" (or should I call it tag arithmetic) does the same as promote_rule but without implicitly assuming the operation. Regarding 2, I was thinking about the handling of number element types only. Just curious, have you ever considered to parametrise arrays by their size?

@lindahua I can see that you have had pretty much the same idea. To me it appears pretty generic that you can compute with the types, but I'm aware that I don't really understand the scope of this.

@ivarne Good to know. I expected this for the master branch and pull requests but not the other branches. I only pushed to JuliaLang to avoid compiling and running tests on my old laptop on which I was watching Barca-City.

@ivarne
Copy link
Member

@ivarne ivarne commented on 1257643 Mar 13, 2014

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You can set up travis testing in your forked repo also. It is just a few clicks away, if you login with github and go to https://travis-ci.org/profile.

@andreasnoack
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JeffBezanson I have been thinking a little more about this. I can use #6118 as case to follow up on your second item. The main reasons why it is a problem that arrays don't have their size as parameters is the use of zero and one to determine the correct return eltype. That was the reason for #4796 because the promotion rule used to be typeof(one(T)*one(S)+one(T)*one(S)). You fixed the issue by avoiding one and instead defining a new promotion function arithtype to deal with bools. My idea is to have +(Type,Type) and *(Type,Type) instead of having promote_type(Type,Type) and arithtype(Type). If the return eltype in generic_matvecmul was determined by T*S+T*S it would solve #6118 as well #4796. There might be a smarter solution, but at least I don't think that the array size is problem for type arithmetic.

@timholy
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My idea is to have +(Type,Type) and *(Type,Type) instead of having promote_type(Type,Type) and arithtype(Type)

I agree that this is surely the direction we need to move in. SIUnits is a great example: 1Meter*1Meter is a completely different type from 1Meter/1Meter and from 1Meter+1Meter, so there's no way we'll be able to use a single type for all arithmetic operations.

@lindahua
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Until Julia provides/exposes the machinery to obtain the result type of a function, I agree that fun(Type) and fun(Type, Type) etc are the best way I can see.

@JeffBezanson
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think that's a good way to do it, because there is duplication. There's no guarantee that Type+Type will tell you the correct type of x+y for actual values. It would be better to have a non-evaluating typeof, and use @typeof(a[1]+b[1]), which would not actually access the arrays.

@andreasnoack
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you provide an example of Type+Type giving the wrong type? I don't question your statement, but I cannot come up with one.

@JeffBezanson
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's just that in the above scheme, you have to define both +(::Type{A}, ::Type{B}) and +(::A, ::B), which makes it possible for them to be different. It should be designed so that this is not even possible.

@lindahua
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be great if a macro such as @typeof is provided. However, it is unclear to me how such a macro can be implemented, as in the macro expansion stage, the type information of the involved variables has not been available.

@ivarne
Copy link
Member

@ivarne ivarne commented on 1257643 Mar 16, 2014

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@lindahua see Jeff's comment in: #6114 (comment).

@lindahua
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, then I am looking forward to this @typeof macro.

@andreasnoack
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok. I see the problem. I still like the idea of using the types instead of the instances so maybe the definition of +(x::Float64,y::Float64)=box(...) could also make the definition +(::Type{Float64},Type{Float64})=@typeof(box(...)). An example of the problem with the instance approach could be

A=[1 0;0.0 1.0]
b=[1,0.5]

here A*b would construct an Array{Int,1} for the result if the type was determined from A[1] and b[1].

@JeffBezanson
Copy link
Member

@JeffBezanson JeffBezanson commented on 1257643 Mar 16, 2014 via email

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@andreasnoack
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It was meant as a device to access the return type information from the normal definition without instances. It also seemed to me as a natural continuation of definitions like promote_type(arithtype(T),arithtype(S)).

What I would like is some way to determine the return type of array computations. The solutions right now are not sufficient as far as I can see. promote_type generally doesn't give the right type, constructing instances from one(Type) doesn't work for arrays and even for numbers they might not be correct because of the ambiguity of one(Real) and using a[1] is also problematic as the example in the last comments showed, and it also doesn't work for empty arrays.

Please sign in to comment.