problem with cumsum

[my-little-repository]

There seems to be a type problem with the cumsum function when passed an union type.

 julia> f(x) = (x == 0 ? 0 : 0.5) 
 f (generic function with 1 method)

 julia> cumsum([f(x) for x = 0:1])     # Array{Union{Float64,Int64},1}
 ERROR: InexactError()
  in convert at int.jl:209
  in cumsum_pairwise! at arraymath.jl:431
  in cumsum! at arraymath.jl:445
  in anonymous at no file

 julia> sum([f(x) for x = 0:1])        # no problem here
 0.5

The error does not appear with the function f(x) = (x == 0 ? 0 : 1.0). An Int array is returned.

[aaneto]

Well, I'm no core developer but,
I've looked on int.jl:209 and found this:

for to in (Int8, Int16, Int32, Int64)
    @eval begin
        convert(::Type{$to}, x::Float32) = box($to,checked_fptosi($to,unbox(Float32,x)))
        convert(::Type{$to}, x::Float64) = box($to,checked_fptosi($to,unbox(Float64,x))) <- line 209
    end
end

Well, julia is on the case of convert(Int64, x::Float64) setting the function to box(Int64,checked_fptosi(Int64,unbox(Float64,x))) which I believe is what is happening, since cumsum tries to sum 0 to 0.5, then the convert function is failing:

convert(Int64, 0.5)
ERROR: InexactError()
 in convert at ./int.jl:209

unbox(Float64, 0.5)
0.5
checked_fptosi(Int64, unbox(Float64, 0.5))
ERROR: InexactError()
 in anonymous at ./no file

For what it looks like, unbox(Float64, 0.5) is returning a float, and it shouldn't(since it's inside a checked_fptosi).
Interesting fact that julia tries to convert to Integer, instead of converting 0 to Float64.

I also got errors when doing this:

x = Array(Real, 2);
x[1] = 0;
x[2] = 0.5;
cumsum(x)
ERROR: InexactError()
 in convert at ./int.jl:209

Note: I'm on 'v"0.5.0-dev+289"'.

[simonster]

It's intended behavior that converting a non-integer floating point number to an integer throws an InexactError. I'm pretty sure that the problem is that cumsum decides the type of the output based on the first element of the input.

[tkelman]

The first element thing isn't right for heterogeneous arrays, we should maybe do some better specialization there to only do that when the element type is concrete?

There was an argument a while back in #9665 about whether cumsum should promote to larger types like sum does to make it less likely to overflow, I personally don't think that's right since it returns an array rather than a scalar and the storage could grow pretty substantially.

[simonster]

Yes, we should do something better with heterogeneous arrays, both here and with mapreducedim. The right solution is probably to do what map does and start with an array that is typed according to the first element and then widen its type as necessary.

[my-little-repository]

When it comes to array, cumsum sends [x_1,....x_n] to [x_1, x_1+x_2, ...., x_1+x_2+..+x_n]. Since + overflows instead of converting to big for int and floats, my guess is that the type of the array returned by cumsum should be equal to the type of the array passed as an argument. + can be overloaded though.

I must say that the following line, referenced above, is a mystery to me. Why adding v[1] to itself and taking its type?

  _cumsum_type(v) = typeof(v[1]+v[1])

[KristofferC]

I guess all types don't have to be closed during addition.

[andreasnoack]

E.g. Bool

Den onsdag den 23. september 2015 skrev Kristoffer Carlsson <
notifications@github.com>:

I guess all types doesn't have to be closed during addition.

—
Reply to this email directly or view it on GitHub
#13244 (comment).

[hayd]

Maybe cumsum (and friends) could optionally take a first argument of the type you want to reduce on.

[tkelman]

You can do cumsum! into an array with element type of your choosing

[hayd]

Interestingly I think this actually boils down to:

julia> typeof([f(x) for x = 0:1])
Array{Union{Float64,Int64},1}

julia> Base._cumsum_type([f(x) for x = 0:1])  # IMO expect Float64
Int64

julia> typeof(zero(Union{Float64,Int64}))  # stems from
Int64

which is due to this definition (since zero just converts 0):

convert{T}(::Type{T}, x::T) = x

The following definition fixes it but with ambiguity warnings galore:

Base.convert{T}(U::UnionType, x::T) = Base.convert(promote_type(U.types...), x)

gives:

julia> cumsum([f(x) for x = 0:1])
2-element Array{Float64,1}:
 0.0
 0.5

[stevengj]

I think this is a bug in _cumsum_type; should be fixed in #18364.

[stevengj]

Closed by #18364.