New array views based on stagedfunctions

[timholy]

Here's my updated proposal for an Array view type. With many different implementations and names floating around (ArrayViews, View, ArrayViewsAPL), to distinguish this one I decided on a bold and innovative name: SubArray.

However, these aren't the SubArrays your grandma played with out behind the old tire swing:

• They are based on stagedfunctions. Both construction and cartesian indexing should be fast. (Haven't had time to benchmark yet.) I expect the performance to still be slightly behind arrays, and it will likely stay that way until tuples get inlined in types.
• In addition to supporting "strided" views, these also support non-strided views (a Vector{Int} along any axis)
• Unlike our current SubArrays, they support all kinds of whacky stuff, like changing dimensionality (A = reshape(1:120,5,4,6); B = sub(A, 1:5, 3:18)).

I realized only recently that these last two properties are required if we're going to consistently return views from indexing. Getting this working properly was the trickiest part.

This container also shares a number of features with our current SubArrays:

• They support arbitrary AbstractArrays
• They support slices (dropping scalar-indexed dimensions) as well as sub/getindex style indexing.

This doesn't make the transition to returning views from indexing. I think we should merge the container and #8432 first, give them a little shaking out, and see where we stand. I already know of a couple likely performance bottlenecks, which I've annotated with comments.

This gets much of the way through all of our tests, but seems to uncover a bug with stagedfunctions in the broadcast tests. So this shouldn't be merged until that gets fixed. (I'll document that later after I've had a little time to investigate.)

@simonster, merge_indexes_div would be fertile territory for your contemplated div speedups.

Also, over at https://github.com/timholy/ArrayViewsAPL.jl (where all the development actually happened), I have a much-expanded set of tests. If desired, I can merge these into base.

[timholy]

Oh, also, I really don't think it's practical to switch to returning views until we have #8227.

[timholy]

I should add that many features of the design are similar to #8235, which you can read for an overview. But basing it on stagedfunctions allows one to avoid the difficult choices between providing enough types and increasing the load time.

[andreasnoack]

@timholy I think it is time to synchronize our efforts towards array views in base. As far as I can see, your and @jakebolewski's and my efforts have been roughly complementary, but we are approaching overlapping work. Our work has been based on @lindahua's package, but much of the work doesn't depend on the choice of view type.

It has proven to be two different problems to introduce one-dimensional and multidimensional range indexing in base. Much of the string related code and the code in inference.jl critically depend on the implicit copy when indexing. Try changing the getindex(A::Array, I::UnitRange{Int}) method in array.jl to return a view. At this point @jakebolewski and I have added enough copys in our branches anj/arrayviews2 and jcb/arrayviews to get Julia to build with one dimensional UnitRange views.

We have also started work on multidimensional views and here we are probably overlapping with your efforts. An important step for getting this to work for @lindahua's implementation is to disable the colon rewrite in julia-syntax.scm because his code dispatches on Colon. This change requires some changes in your code for multidimensional indexing. I tried, but lost to cartesian.jl and right now I'm using a hack to bypass @nloops and just build the expressions for getindex in a loop.

We should try to figure out if it is easier to merge this pull request and add @jakebolewski's and mine on top on that or the other way around. I think that it is better to merge your work and then add @jakebolewski's and mine on top, i.e. return views from getindex and use @lindahua's views when the underlying array allows. However, it would be extremely useful if you could include the changes to julia-syntax.scm and the required adjustments of multidimensional.jl to handle Colon. Right now, I'm using

decolon(A::AbstractArray, i::Colon, dim::Real) = 1:size(A, dim)
decolon(A::AbstractArray, i, dim::Real) = i
decolon1(A::AbstractArray, i::Colon, dim::Real) = 1:length(A)
decolon1(A::AbstractArray, i, dim::Real) = i

but there is probably a better solution.

[timholy]

If we have both types, are you planning on writing an implementation of slice for a StridedView input? What are the obstacles you see with having a single container for all cases?

Getting the colon translation out of the parser might still be a good idea, but with this approach it isn't actually necessary unless you want to extend it to handle Contiguous views (which you could do by adding a 5th boolean type parameter). My hope, though, is that this may soon not be needed anymore (see #8432).

However, the work you folks did to make strings.jl & inference.jl work will indeed be general. Exciting to hear!

[andreasnoack]

It appears to me that the contiguous property is pretty useful to be able to dispatch on, although I'm not sure that I'll use it that much. I'm mainly interested in avoiding copying when doing BLAS operations on a matrix.

I must admit that I don't fully understand cartesian.jl, and maybe I'm not the only one. I can see the usefulness in having a way to avoid the splatting argument, but I just don't get this linear versus cartesian indexing even after reading through the issues. Googling "cartesian indexing" only gives links to your issues and something I believe is irrelevant.

First of all, is what you call the general case a situation where you want to compute for something like?

foo(x)
for i = 1:length(x)
bar(x[i])
end

If so, it appears to me that this is a special kind of indexing where you want to go through all elements one by one.

Next, could you explain (or point to an explanation of) why your approach is "Cartesian"? I don't understand the reference.

Regarding the slice question, I haven't really thought that much about it to be honest. I think I prefer slices, but I have mainly been concerned about getting getindex to return some kind of view as soon as possible such that we can start the transition and people can get used to it.

[timholy]

Sorry, I should have explained better; I didn't realize that was the issue. By cartesian indexing, I don't mean "uses Base.Cartesian," I mean accessing the elements of a 2d array like this:

for j = 1:size(A,2), i = 1:size(A,1)
    s += A[i,j]
end

rather than as

for i = 1:length(A)
    s += A[i]
end

The reason is simple: in the general case (meaning, where A may not be contiguous), the latter gets translated into

for i = 1:length(A)
    i_2, i_1 = divrem(i-1, size(A,1))
    s += A[i_1+1, i_2]
end

and divrem is appallingly slow.

EDIT:
This is why I think any "general purpose" container that can handle AbstractArrays, many of which do not have efficient linear indexing, needs to be based on cartesian indexing.

This is also why #8432 is important: it gives users "easy" cartesian indexing without needing to know a darned thing about Base.Cartesian.

[timholy]

Regarding slice, I simply mean that we currently have both sub and slice; independently of what gets decided for how getindex should behave, both sub and slice can be (and are) directly called by users. So we need to support both operations for all of our view containers, no matter what. This is why I think it's better to use a single container, if possible, to avoid the annoyance of having to implement both sub and slice for 3 different container types.

I use this all the time in my image processing code, because sometimes you want to preserve dimensions (sub) and sometimes you want to discard them (slice). If both are available, then even the people who are unhappy with how the getindex debate works out can still have something that works for them, just with ever-so-slightly-longer notation.

[simonster]

But if we have only one type of view, that means we give up performance for linear indexing of contiguous views of arrays, right? I think we should either do everything we can to make linear indexing as fast as possible or get rid of it entirely.

[timholy]

We should definitely make linear indexing as fast as possible. But independently, you can integrate ArrayViews's ContiguousView into this SubArray variant: just add a 5th boolean parameter. The stagedfunction can do the type analysis to determine whether to set it to true.

Doing this automatically for A[:, i] would require the parser changes that @jakebolewski and @andreasnoack have been working on, since you'd need to pass Colon to the constructor.

[timholy]

But my recommendation would be to fix #8504, then merge this, then implement contiguous views, then work on the parser.

[simonster]

That sounds like a good plan to me.

[timholy]

We finally have green. So, I've removed the WIP tag. Now is a good time to complain about anything you don't like. Along those lines, I'll point out that this does not support the sub(A, ::Tuple) syntax. I always assumed that was an internal interface (it's not documented), but I can restore it if necessary.

Meanwhile, I'll look into what would be needed to add support for efficient linear indexing à la ContiguousViews. (We may also be able to support things like row slices.)

[ViralBShah]

Hooray!

[lindahua]

@timholy what about having a benchmark to compare these with ArrayViews and ArrayViewsAPL?

I think it is good to go if the benchmark shows that it is fast enough.

[timholy]

I've benchmarked similar things before, but yeah, I should do that, because now this is serious 😄. #8432 has some semi-relevant benchmarks, but it's with the old SubArrays.

Until we include efficient linear indexing for those types that can support it, this will not be as efficient as linear indexing on a ContiguousView. (It won't even be close.) My main hope is that with #8432, linear indexing may not be all that important anymore. But I suspect people will sometimes want to use linear indexing, and we should explore whether we can implement the strengths of your approach in a more general container type (or add both sets of types and figure out how to switch among them). Since a lot of type analysis can be done via stagedfunctions, I'm optimistic it should be possible (and perhaps even extendable to row slices, etc.)

[timholy]

OK, time for an update here. With the merger of #8432, things are looking pretty good. Here's a benchmark script. If one just tallies up the differences, I find 4 cases in which this is at least 4x faster than ArrayViews, and 1 case in which ArrayViews is 4x faster than this. There are also several factor-of-two differences that mostly favor this PR. Otherwise the two are pretty similar.

Construction
buildcols_view: elapsed time: 0.019044581 seconds (19200096 bytes allocated)
buildcols_sub: elapsed time: 0.018561561 seconds (26400096 bytes allocated)
buildcols_slice: elapsed time: 0.013784228 seconds (26400096 bytes allocated)
buildrows_view: elapsed time: 0.023215199 seconds (33600096 bytes allocated)
buildrows_sub: elapsed time: 0.022335015 seconds (36000096 bytes allocated)
buildrows_slice: elapsed time: 0.014829145 seconds (26400096 bytes allocated)
buildsubset_view: elapsed time: 0.019534433 seconds (30240096 bytes allocated)
buildsubset_sub: elapsed time: 0.031503118 seconds (51840096 bytes allocated)
buildsubset_slice: elapsed time: 0.020893927 seconds (38880096 bytes allocated)
Access
30-element ContiguousView{Float64,1,Array{Float64,2}}
    mdsum: elapsed time: 0.210575476 seconds (16 bytes allocated)
sumlinear: elapsed time: 0.15027899 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.150343726 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.150261564 seconds (16 bytes allocated)
30-element SubArray{Float64,1,Array{Float64,2},(Colon,Int64)}
    mdsum: elapsed time: 0.148767953 seconds (16 bytes allocated)
sumlinear: elapsed time: 0.109161175 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.109284359 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.110314285 seconds (16 bytes allocated)
30-element SubArray{Float64,1,Array{Float64,2},(Colon,Int64)}
    mdsum: elapsed time: 0.110716632 seconds (16 bytes allocated)
sumlinear: elapsed time: 0.10902906 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.108984881 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.109754716 seconds (16 bytes allocated)
1x30 StridedView{Float64,2,0,Array{Float64,2}}
    mdsum: elapsed time: 1.087570857 seconds (16 bytes allocated)
sumlinear: elapsed time: 2.55995489 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.686259948 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.225551411 seconds (16 bytes allocated)
1x30 SubArray{Float64,2,Array{Float64,2},(UnitRange{Int64},Colon)}
    mdsum: elapsed time: 0.251582036 seconds (16 bytes allocated)
sumlinear: elapsed time: 3.198257375 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.276703199 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.318401684 seconds (16 bytes allocated)
30-element SubArray{Float64,1,Array{Float64,2},(Int64,Colon)}
    mdsum: elapsed time: 0.154634719 seconds (16 bytes allocated)
sumlinear: elapsed time: 0.152036507 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.152012612 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.147994735 seconds (16 bytes allocated)
30x2 ContiguousView{Float64,2,Array{Float64,2}}
    mdsum: elapsed time: 0.132879825 seconds (16 bytes allocated)
sumlinear: elapsed time: 0.107304793 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.114413543 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.221376907 seconds (16 bytes allocated)
30x2 SubArray{Float64,2,Array{Float64,2},(Colon,UnitRange{Int64})}
    mdsum: elapsed time: 0.168667217 seconds (16 bytes allocated)
sumlinear: elapsed time: 3.129825562 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.107680293 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.204832289 seconds (16 bytes allocated)
2x30 StridedView{Float64,2,1,Array{Float64,2}}
    mdsum: elapsed time: 1.390568064 seconds (16 bytes allocated)
sumlinear: elapsed time: 2.45289969 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.49883066 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.193457333 seconds (16 bytes allocated)
2x30 SubArray{Float64,2,Array{Float64,2},(UnitRange{Int64},Colon)}
    mdsum: elapsed time: 0.150158097 seconds (16 bytes allocated)
sumlinear: elapsed time: 3.162968933 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.13820917 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.197927709 seconds (16 bytes allocated)
28x28 StridedView{Float64,2,1,Array{Float64,2}}
    mdsum: elapsed time: 1.186577744 seconds (16 bytes allocated)
sumlinear: elapsed time: 2.576119829 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.14964565 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.274510347 seconds (16 bytes allocated)
28x28 SubArray{Float64,2,Array{Float64,2},(UnitRange{Int64},UnitRange{Int64})}
    mdsum: elapsed time: 0.243553184 seconds (16 bytes allocated)
sumlinear: elapsed time: 3.124831148 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.15031546 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.307520535 seconds (16 bytes allocated)

By comparison with ArrayViews, the only serious weakness is linear indexing with the 30x2 ContiguousView. I'm going to get started on implementing fast linear indexing for cases that can support it. But since we now have very good alternatives to linear indexing, there is part of me that wonders if it is worth it. The main cost will be growing the container by 2 Ints (I plan to implement fast "row slice" indexing too). On balance, I suspect we'll want it; linear indexing is occasionally very useful.

[timholy]

OK, I incorporated efficient linear indexing. Seems to work, and it fixes the last cases in which performance is notably poor compared to ArrayViews:

Construction
buildcols_view: elapsed time: 0.016437897 seconds (19215872 bytes allocated)
buildcols_sub: elapsed time: 0.112973297 seconds (44640096 bytes allocated)
buildcols_slice: elapsed time: 0.112691719 seconds (44640096 bytes allocated)
buildrows_view: elapsed time: 0.020025157 seconds (33600096 bytes allocated)
buildrows_sub: elapsed time: 0.167410432 seconds (50400096 bytes allocated)
buildrows_slice: elapsed time: 0.1549283 seconds (40800096 bytes allocated)
buildsubset_view: elapsed time: 0.026403215 seconds (30240096 bytes allocated)
buildsubset_sub: elapsed time: 0.156244869 seconds (68160096 bytes allocated)
buildsubset_slice: elapsed time: 0.16171357 seconds (68160096 bytes allocated)
Access
30-element ContiguousView{Float64,1,Array{Float64,2}}
    mdsum: elapsed time: 0.224426815 seconds (16 bytes allocated)
sumlinear: elapsed time: 0.150395494 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.150353279 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.15227818 seconds (16 bytes allocated)
30-element SubArray{Float64,1,Array{Float64,2},(Colon,Int64),2}
    mdsum: elapsed time: 0.154745698 seconds (16 bytes allocated)
sumlinear: elapsed time: 0.151981558 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.150315828 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.150466393 seconds (16 bytes allocated)
30-element SubArray{Float64,1,Array{Float64,2},(Colon,Int64),2}
    mdsum: elapsed time: 0.154488798 seconds (16 bytes allocated)
sumlinear: elapsed time: 0.15030374 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.150345076 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.150308991 seconds (16 bytes allocated)
1x30 StridedView{Float64,2,0,Array{Float64,2}}
    mdsum: elapsed time: 1.522944058 seconds (16 bytes allocated)
sumlinear: elapsed time: 3.036434464 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.903976489 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.312141298 seconds (16 bytes allocated)
1x30 SubArray{Float64,2,Array{Float64,2},(UnitRange{Int64},Colon),2}
    mdsum: elapsed time: 0.269215442 seconds (16 bytes allocated)
sumlinear: elapsed time: 0.151958861 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.277800317 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.320623793 seconds (16 bytes allocated)
30-element SubArray{Float64,1,Array{Float64,2},(Int64,Colon),2}
    mdsum: elapsed time: 0.155687126 seconds (16 bytes allocated)
sumlinear: elapsed time: 0.153157405 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.151998673 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.153909207 seconds (16 bytes allocated)
30x2 ContiguousView{Float64,2,Array{Float64,2}}
    mdsum: elapsed time: 0.190875125 seconds (16 bytes allocated)
sumlinear: elapsed time: 0.15030275 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.160453244 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.310126423 seconds (16 bytes allocated)
30x2 SubArray{Float64,2,Array{Float64,2},(Colon,UnitRange{Int64}),2}
    mdsum: elapsed time: 0.237873813 seconds (16 bytes allocated)
sumlinear: elapsed time: 0.150376657 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.153586218 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.307795716 seconds (16 bytes allocated)
2x30 StridedView{Float64,2,1,Array{Float64,2}}
    mdsum: elapsed time: 1.488393786 seconds (16 bytes allocated)
sumlinear: elapsed time: 3.351965134 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.532021058 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.273943452 seconds (16 bytes allocated)
2x30 SubArray{Float64,2,Array{Float64,2},(UnitRange{Int64},Colon),1}
    mdsum: elapsed time: 0.218928902 seconds (16 bytes allocated)
sumlinear: elapsed time: 3.102060723 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.193469576 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.27718949 seconds (16 bytes allocated)
28x28 StridedView{Float64,2,1,Array{Float64,2}}
    mdsum: elapsed time: 1.49909814 seconds (16 bytes allocated)
sumlinear: elapsed time: 3.336110625 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.162210062 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.304791965 seconds (16 bytes allocated)
28x28 SubArray{Float64,2,Array{Float64,2},(UnitRange{Int64},UnitRange{Int64}),1}
    mdsum: elapsed time: 0.276495672 seconds (16 bytes allocated)
sumlinear: elapsed time: 4.145109959 seconds (16 bytes allocated)
  sumcart: elapsed time: 0.153321363 seconds (16 bytes allocated)
  sumiter: elapsed time: 0.275343035 seconds (16 bytes allocated)

Note in particular the excellent linear indexing performance for the 1x30 SubArray, which is not something that is currently matched in ArrayViews. All the magic, of course, is thanks to stagedfunctions, which make this stuff easy; at the time ArrayViews.jl was written, I have no doubt it would have been really, really hard to pull this off just in terms of dispatch. New technology is nice.

Performance aside, I suspect I should write a few more tests to make sure this is really working properly, but it does pass the test suite on my machine (and hopefully on Travis, too).

Because I don't want to become the sole maintainer of Julia's Array frameworks for perpetuity 😄, I'm also writing up a devdoc on how all this works (it's pretty darn straightforward, albeit somewhat tedious). I'll try to push those over the next couple of days, but meanwhile feel free to kick the tires and comment.

Can anyone else feel the nirvana soaking in?

[johnmyleswhite]

🍰

[timholy]

We have new SubArrays. While the timings above emphasized the comparison to ArrayViews, it's worth mentioning that compared to the SubArrays we have in julia 0.3, these are ~2-10x faster to use for indexing, and 20-50x faster to construct. Those factors should be large enough to make a difference in many applications.

@andreasnoack and @jakebolewski, as discussed I turn phase II (making A[:, 3] return a view) over to you. Sorry it has taken so long. Note that I documented many aspects of the internal design, which hopefully will help with any questions you might have (and happy to answer questions/help as appropriate).

Meanwhile, I'm off to begin the process of integrating the new iteration framework more deeply into base, as it should enable improvements in several functions. I'll probably be experimenting with moving some things forward in the bootstrap sequence (mostly, cartesian.jl and the multidimensional iterators). Hopefully that won't conflict with your own work.

[JeffBezanson]

Sorry I haven't been paying attention here, but this looks awesome.

As part of this, I'd like to request that cartesian.jl not use regexes to rename symbols. That was where I encountered #9011, and it's not a good dependency for core computational stuff like this to have.

[ViralBShah]

👏

[timholy]

I'd like to request that cartesian.jl not use regexes to rename symbols... it's not a good dependency for core computational stuff like this to have.

That sounds reasonable.

[JeffBezanson]

There seem to be lots of crashes during jl_instantiate_staged after this, e.g. https://travis-ci.org/JuliaLang/julia/jobs/41586398 . Don't know if this caused it, but it's certainly triggering it.

[lindahua]

@timholy I wholeheartedly think that this is a major step forward for Julia. Thanks for all the efforts!

[timholy]

Thanks, everyone!

There seem to be lots of crashes during jl_instantiate_staged

That's weird; I pushed and forced-pushed this branch something like 20 times, and (once I got the bugs worked out) got a green check basically every time. So I felt pretty confident about it, despite the fact that bugs in stagedfunctions (#8504, #8505, #8944, plus my grant 😄) were the main thing that held this PR up.

But yes, the SubArray work exercises stagedfunctions with a vengence, so it's not unlikely that it's uncovering problems. I also have another PR coming that will amplify their usage dramatically, and much more widely in typical usage (not that much code yet uses SubArrays, since they have been so problematic). Meanwhile, if anyone can come up with a reproducible way to make the crash happen, it would be great.