Domain error when inferring array type in 0-length degenerated case

[fverdugo]

Hi!

do you plan to support degenerated cases like the one below in the future?

Just for curiosity, how would you fix this in the library? Is it possible to fix it without calling something like Base._return_type?

Thanks for the help!

julia> using MappedArrays

julia> a = fill(1,0)
0-element Array{Int64,1}

julia> f(x) = sqrt(x-1)
f (generic function with 1 method)

julia> mappedarray(f,a)
ERROR: DomainError with -1.0:
sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
Stacktrace:
 [1] throw_complex_domainerror(::Symbol, ::Float64) at ./math.jl:31
 [2] sqrt at ./math.jl:493 [inlined]
 [3] sqrt at ./math.jl:519 [inlined]
 [4] f at ./REPL[3]:1 [inlined]
 [5] mappedarray(::typeof(f), ::Array{Int64,1}) at /home/fverdugo/.julia/packages/MappedArrays/8HESW/src/MappedArrays.jl:61
 [6] top-level scope at REPL[4]:1

[timholy]

One thing we should probably support (but don't currently) is this:

julia> m = MappedArray{Float64}(f, a)
ERROR: MethodError: no method matching MappedArray{Float64,N,A,F,Finv} where Finv where F where A<:AbstractArray where N(::typeof(f), ::Array{Int64,1})
Stacktrace:
 [1] top-level scope at REPL[4]:1

It would just be a matter of adding a couple of methods.

Without requiring the user to provide the answer, for cases like yours I am not aware of a solution that works in general that doesn't use Base._return_type. I'd be willing to switch to that, if it doesn't mess any existing usages up.

[fverdugo]

Thanks for the answer!

I also see another possible problem when inferring the eltype, now for non-empty arrays. If the input array has an eltype that is not a concrete one, then inferring the mapped eltype from the first array entry is not save.

I am right?

[timholy]

Yes, that's right. Another reason to rework this.

[fverdugo]

Just a final question:

is there any way of using an "inplace" funcion, when working with large arrays of arrays? I would like to avoid allocating memory when getting an entry from the lazy array.

something like

function f!(v::AbstractArray,a::AbstractArray)
  # mutate v from the entries in a
end

# Instead of

function f(a::AbstractArray)
# Allocate v and mutate  it from the entries in a
# return v
end

a = [rand(10,20) for i in 1:10000]
mappedarray(f,a)

Do you think, allocating memory is a problem each time an entry is accessed for large arrays?

Thanks for the help!

[timholy]

The easiest way to do that would be to put the storage and logic into the function you pass in:

julia> using MappedArrays

julia> const buf = zeros(5);

julia> f(x) = (buf .= x.^2; return buf)
f (generic function with 1 method)

julia> aa = [rand(5) for i = 1:10];

julia> m = mappedarray(f, aa);

julia> function sumalot(m)
           s = zero(first(m))
           for i = 1:10000
               s .+= m[(i % (length(m)-1)) + 1]
           end
           return s
       end
sumalot (generic function with 1 method)

julia> sumalot(m)
5-element Array{Float64,1}:
 4308.145615923303 
 4008.7673891697395
 4891.049206648661 
 2737.7823120736575
 3687.6835874596586

julia> @time sumalot(m)
  0.000576 seconds (5 allocations: 288 bytes)
5-element Array{Float64,1}:
 4308.145615923303 
 4008.7673891697395
 4891.049206648661 
 2737.7823120736575
 3687.6835874596586

You could also define a type and a corresponding call function, if you'd prefer to be more explicit about "containment" of the buffer.

[timholy]

For small arrays you can also use StaticArrays. Also note ImageAxes: StreamingContainer, if you have a really big operations you're performing.

[fverdugo]

The problem with a constant global buffer is that it's not thread-save, but certainly a solution for single-threaded applications.

The approach I like more, is to extend the AbstractArray interface with a new method getindex!(buff,index) in the same way you have done in StreamingContainer.

Thanks for the help!

This is all what I needed to know.

Best!

[timholy]

It would be great to fix some of the issues noted. If you don't want to submit a PR to fix them yourself, let's at least leave this open.

[fverdugo]

Are you thinking in a solution using Base._return_type ?

I find it quite dangerous, i.a., it hides errors:

julia> g(x) = foo(x)
g (generic function with 1 method)

julia> g(1)
ERROR: UndefVarError: foo not defined
Stacktrace:
 [1] g(::Int64) at ./REPL[1]:1
 [2] top-level scope at REPL[2]:1

julia> Base._return_type(g,(Int,))
Any

With Base._return_type, the user will never see that foo its not defined. Here not at issue, but in a large complex code it is...

I think, the best option is that the user provides test values for the function if he/she knows that zero(T) will fail, and using perhaps zero(T) by default.

For instance

foo(x) = sqrt(x-1)

# Default test value
testvalue(f::Function,::Type{T}) where T = zero(T)

# User-defined one
testvalue(::typeof(foo),,::Type{T}) where T = zero(T)+one(T)

Then using testvalue in the library to infer the return type. You can even write a meaningful error message if the default test value leads to a DomainError.

This solution is minimally invasive. What do you think?

[timholy]

Sorry for the delay. I like it!