some reductions over dimensions depend too heavily on zero()

[svs14]

Code to reproduce:

# Works
std(Float64[1,2,3], 1)
1-element Array{Float64,1}:
 1.0

# Breaks
std(FloatingPoint[1,2,3], 1)
ERROR: no method zero(Type{Any})
 in varm at statistics.jl:98
 in var at statistics.jl:109

# Julia version
versioninfo()
Julia Version 0.3.0-prerelease+2494
Commit 7d495c0 (2014-04-04 23:25 UTC)
Platform Info:
  System: Linux (x86_64-unknown-linux-gnu)
  CPU: Intel(R) Core(TM) i7-2720QM CPU @ 2.20GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY)
  LAPACK: libopenblas
  LIBM: libopenlibm

[IainNZ]

Narrowing it down...

julia> Base.sumabs2(FloatingPoint[1,2,3],1)
ERROR: no method zero(Type{Any})
 in sumabs2 at statistics.jl:76

[IainNZ]

So

julia> typeof(abs2(FloatingPoint[1,2,3]))
Array{Any,1}

which seems like a failure of abs2 because |FloatingPoint|^2 should still be FloatingPoint

[IainNZ]

Although its a pretty tough type inferenence problem I guess, because you could really have anything after you are done abs2'ing it.

[JeffBezanson]

Or simply:

julia> sum({1 2;3 4},1)
ERROR: no method zero(Type{Any})
 in sum at reducedim.jl:146

julia> mapslices(sum, {1 2;3 4}, 1)
1x2 Array{Int64,2}:
 4  6

[timholy]

For sum we could declare the type as typeof(r[1]+r[1]) (when r is not empty), but the general problem is not easy.

I'm increasingly a fan of the notion that the return type of all functions should be Void. It would make life a lot simpler 😄. There's a (very small) element of seriousness here: the following should work

A = {1 2; 3 4}
R = {0 0}
sum!(R, A, init=false)

but it doesn't (it would have, before the addition of the keyword).

[andreasnoack]

...and the solution with the use of first element doesn't really work either because you could just write

mapslices(sum, {1 2;3 4.2}, 1)
ERROR: InexactError()
 in setindex! at multidimensional.jl:68
 in anonymous at abstractarray.jl:1244
 in cartesianmap at abstractarray.jl:1147
 in mapslices at abstractarray.jl:1237
 in mapslices at abstractarray.jl:1197

It is really hard to support heterogenous arrays for all the numerical stuff. The element type of the array is not really informative when it is abstract. Even if the type inference could return Real here

julia> Base.return_types(abs2, (Vector{Real},))
1-element Array{Any,1}:
 Array{Any,1}

it wouldn't help for a reduction. Should the return type be Int, Float64, BigFloat or something else?

[timholy]

Is there really even a solution for this on the horizon? Since there's no zero(Any), I kind of think we're stuck.

Despite my Void joke above, I am semi-seriously beginning to wonder if the only real solution is a radical one: functions simply aren't allowed to allocate arrays for their return values. See #1115.

[benmoran]

I had been using Julia's matrix-vector multiplication to build up coefficient expressions for JuMP constraints. This used to work, returning an Array of JuMP.GenericAffExpr, but now it fails because there's no zero(Any).

using JuMP
m = Model()
@defVar(m, x[1:4])
vec = x[:] # vec::Array{JuMP.Variable,1}
A = eye(4)
A * vec # ERROR: no method zero(Type{Any}) in generic_matvecmul! at linalg/matmul.jl:321

I also fall foul of sum(X, dim) giving the same error when I try to work around it, though sum(X) works fine.

[JeffBezanson]

Marking as regression.

[mlubin]

@benmoran, in JuMP at least it would be more efficient to use dot as needed, we've overloaded this (and sum) explicitly to work properly. Using generic operations will be inefficient because many temporary objects are created. I agree that generic mat-vec should "just work" though if the necessary operations are defined. We should be able to fix this particular case by defining the needed promotion rules.

[andreasnoack]

@mlubin Yes, for this example, something like

promote_rule(::Type{Float64}, ::Type{Variable})=JuMP.GenericAffExpr{Float64,Variable}

is sufficient. The general case is difficult to handle with the tools we have right now.

[stevengj]

One solution would be a two-pass approach: first, look at all of the elements to figure out the correct type of the reduction, then perform the summation.

The simplest implementation might be something like:

function sum{T}(A::AbstractArray{T}, region)
    z = method_exists(zero, (T,)) ? zero(T)+zero(T) : zero(sum(A))
   _sum!(reduction_init(A, region, z), A)
end

This still doesn't work where T is an abstract type like Real, because zero(Real) is 0 (Int64). However, that could be fixed by making reduction_init accept another argument, the type of the array to allocate (i.e. you would allocate a Real array initialized to 0). In general, I think it would be reasonable for sum{T}(A::AbstractArray{T}, region) to return an Array{T}.