WIP: Disentangle arrays in iteration #15434

timholy · 2016-03-10T03:44:43Z

This PR is designed to allow us to expand our collection of array iterators (most proximally, to solve serious performance problems with ReshapedArrays).
As a motivating example, suppose copy! is implemented like this:

function copy!(dest, src)
    for I in eachindex(dest, src)
        dest[I] = src[I]
    end
    dest
end

Let's suppose that src and dest don't have an efficient iterator in common; in that case, the loop will be slow. This PR replaces such implementations with

function copy!(dest, src)
    for (Idest, Isrc) in zip(eachindex(dest), eachindex(src))
        dest[Idest] = src[Isrc]
    end
    dest
end

leveraging recent improvements (and planned additional improvements) in the performance of zip.

This first commit addresses every use of eachindex in Base. I still have the classic for i = 1:n cases to go. (Where's that :janitor: emoticon when you need it?) I'm putting this out now in case there are reservations. It's also an almost guarantee that some workloads will see performance regressions. I'd rather find out now how bad they are, so: runbenchmarks(ALL, vs = "JuliaLang/julia:master").

tkelman · 2016-03-10T09:43:42Z

base/io.jl

-    for i in eachindex(a)
-        nb += write(s, a[i])
+    for b in a
+        nb += write(s, b)


Would it be a little more readable to rename the input A::AbstractArray then do for a in A ? Having nb and b nearby is a little offputting (nb could be a few characters longer of course). Similar above with the for y in x comprehensions.

I was tempted to do that, but decided to try to blend in with the surrounding code. I'll reverse that decision and change it.

This also replaces several instances of `for i in eachindex(A); a = A[i]` with `for a in A`. The latter is presumably easier for automatic bounds-checking removal.

timholy · 2016-03-10T12:27:27Z

For those who arrived here looking to help, here are some guidelines I'm using:

When you don't actually need the index, for a in A iteration is often nicer (see below for further discussion). In theory, it's also easier to handle in automatic bounds-checking removal, so it might lead to performance improvements without the need for @inbounds annotations.
When possible, avoid assuming that arrays are efficiently indexed by integers (use eachindex where it makes sense).
When possible, avoid assuming that arrays are indexed from 1:length(n) (yes, I'm hoping to implement Fortran-style OffsetArrays someday). This will probably never change for julia's base arrays, and in some cases it's impractical to be general about this, but in many situations you don't actually care what the values of the indexes are, as long as they are appropriate for that particular array.
When the algorithm explicitly requires that the index i in for i = 1:n is an integer with values in that range, don't change it. (In primes.jl I switched some eachindex usages back to for i = 1:n).
When a function both creates the array and iterates over it, and you can guarantee that the arrays are plain-old Arrays, no changes are needed (there are no concerns about generality). inference.jl is a good example (I don't plan on changing any cases of iteration there). By contrast, when arrays are created via similar, then you can't be certain that you know their most efficient style of iteration, so in such cases you should generalize the indexing.

There's one important sticking point: comprehensions. See #10523. Currently [f(i,j) for i = 1:m, j = 1:n] creates a 2d array but [f(I) for I in CartesianRange((m,n))] creates a 1d array. That might change, and there's some reason to suspect that this cleanup would be better done once that decision is made.

Finally, it's a very minor issue, but I'm torn about the best choice in some situations. Should it be like this:

for (idest, iA, iB) in zip(eachindex(dest), eachindex(A), eachindex(B))
    dest[idest] = f(A[iA], B[iB])
end

or

for (idest, a, b) in zip(eachindex(dest), A, B)
    dest[idest] = f(a, b)
end

The second is more readable, but I'm tending to do the first because of the issues discussed in #15356.

timholy · 2016-03-10T12:41:04Z

Looks like force-pushing nixed the still-pending benchmark report. Hopefully it will still appear at https://github.com/JuliaCI/BaseBenchmarkReports?

nalimilan · 2016-03-10T13:16:46Z

I suggest you make something like 5-10 groups of sources files (e.g. in lexical order), so that each person willing to help can reserve a group without stepping on others' toes. I'm OK with doing some of this work.

KristofferC · 2016-03-10T13:48:55Z

Should people open PRs against this branch or just make a new PR against master?

eschnett · 2016-03-10T16:08:33Z

@timholy Fortran-style array indexing: See https://github.com/eschnett/FlexibleArrays.jl. From the examples:

# A 3d array with lower index bounds 0
typealias Arr3d_lb0 FlexArray(0, 0, 0)
a3 = Arr3d_lb0{Float64}(9, 9, 9)

Here the lower index is part of the type, i.e. has no run-time overhead. You can also write:

# A generic array, all bounds determined at creation time
typealias Arr4d_generic FlexArray(:, :, :, :)
a4 = Arr4d_generic{Float64}(1:10, 0:10, -1:10, 15:15)

Fixed and open indices can also be mixed. Suggestions for improved syntax are always welcome.

nanosoldier · 2016-03-10T17:46:46Z

Your benchmark job has completed, but something went wrong when trying to upload the result data. cc @jrevels

timholy · 2016-03-10T18:13:23Z

Remarkable. No cases that seem like real regressions.

timholy · 2016-03-10T18:30:25Z

I would make separate PRs---they don't need to intersect at all. This can be the "organizer" and just hold the eachindex changes.

For comprehensions, maybe best for now would be to mark any that need fixing with # fixmecomp.

Blocks of tasks (first check box is "claimed", second checkbox is "done"):

musm · 2016-03-10T19:58:12Z

base/primes.jl

@@ -39,7 +39,7 @@ function _primesmask(lo::Int, hi::Int)
    sieve = ones(Bool, whi - wlo)
    hi < 49 && return sieve
    small_sieve = _primesmask(isqrt(hi))
-    @inbounds for i in eachindex(small_sieve)
+    @inbounds for i = 1:length(small_sieve)  # don't use eachindex here


why not here?

It calls a function, wheel_prime(i), on i. This means that you need to insist on a particular meaning for i, in which case it's better to be explicit about what i needs to be.

musm · 2016-03-10T19:59:37Z

@timholy recommend an easy one I could try? (noob contributor, but I think some of this seems easy enough)

timholy · 2016-03-10T20:34:09Z

I can't really predict without looking, but as a random guess how about iobuffer.jl, io.jl, iostream.jl?

Thanks for volunteering!

dcarrera · 2016-03-11T13:43:02Z

Maybe we should announce what file we are going to work on. I call dibs io.jl

Let me see if I understand. In io.jl I see the following on line 145:

function write(s::IO, a::AbstractArray)
    nb = 0
    for i in eachindex(a)
        nb += write(s, a[i])
    end
    return nb
end

The for loop can be trivially rewritten as:

    for i in a
        nb += write(s, i)
    end

Is this the sort of thing we are trying to do? Assume I don't know anything. What do I do now?

How do I pull this branch? A simple git pull will just give me master.
To compile Julia do I just run make?
How do I run the benchmarks? The function runbenchmarks does not seem to be defined.

The idea then is that I'd do the benchmarks before and after my change and report the result here? How do I submit a patch?

KristofferC · 2016-03-11T14:05:12Z

@dcarrera This has some good general information: https://github.com/JuliaLang/julia/blob/master/CONTRIBUTING.md

You pull the branch by first running git fetch and then checking out this branch, for example git checkout -b branch_name origin/teh/separate_indexes
Yes. You might need a make clean depending on what version you built before.
Regarding benchmarks, the benchmarks are actually defined in a separate repo (https://github.com/JuliaCI/BaseBenchmarks.jl) and are usually run when a PR has been opened.

meggart · 2016-03-11T15:06:28Z

I volunteer, too. I start working on the first block, abstractarray.jl ... arraymath.jl.

timholy · 2016-03-11T15:15:59Z

I'll put the "claimed" checkmarks next to those two. Thanks!

@dcarrera, any loop using eachindex is already addressed in this PR. It's just loops written as for i = 1:n that need attention. To reduce confusion, I'll merge this once I re-run that one weird travis result.

WIP: Disentangle arrays in iteration

timholy · 2016-03-11T16:28:27Z

OK, this is merged to reduce the likelihood of conflicts.

I encourage folks to check out #15459: new syntax eachindex(A, d) is coming, and it might be useful. (Don't worry about it if you've already gotten started, just let me know which files need a second look and I'll do it.) I deliberately excluded linalg/ and sparse/ from the list above, because I could tell I would need to do a bit more thinking to figure out how to handle those. I'm pretty excited about the solution, though we'll see how others like it.

dfdx · 2016-03-11T18:30:24Z

I've been always interested in REPL-related stuff. I'll take REPLCompletions.jl ,..., sharedarray.jl.

timholy · 2016-03-11T19:39:09Z

Great!

meggart · 2016-03-11T20:47:39Z

A question how far we should push this. For example in this function:

function ipermutedims(A::AbstractArray,perm)
    iperm = Array(Int,length(perm))
    for i = 1:length(perm)
        iperm[perm[i]] = i
    end
    return permutedims(A,iperm)
end

Here, we obviously assume that perm is indexed from 1:length(perm), which might be wrong if someone uses a kind of OffSetArray for perm. What is the better solution here:

just call perm=collect(perm) before the loop.
Rewrite the loop to something like:

function ipermutedims(A::AbstractArray,perm)
    iperm = Array(Int,length(perm))
    for (i,j) = zip(1:length(perm),eachindex(perm))
        iperm[perm[j]] = i
    end
    return permutedims(A,iperm)
end

meggart · 2016-03-11T21:35:35Z

Sorry for asking again, here is another quite common pattern: What should we do about Iterations with a given offset. For example:

function findnext(A, start::Integer)
    for i = start:length(A)
        if A[i] != 0
            return i
        end
    end
    return 0
end

Is there a version of eachindex that takes an offset as an optional argument? Any other solutions?

pkofod · 2016-03-11T21:50:47Z

irrationals.jl, iterator.jl, libdl.jl
done

I think this is already "done". Very few loops, but

irrationals.jl only has for a in A style loops
iterator.jl only has one loop which is

for i in 1:it.n
    if done(it.xs, xs_state)
        break
    end

    _, xs_state = next(it.xs, xs_state)
end

where i is not used explicitly

libdl.jl has for a in A style loops and one loop where i is used as input to a function

for i in 1:numImages-1
    name = bytestring(ccall(:_dyld_get_image_name, Cstring, (UInt32,), i))
    push!(dynamic_libraries, name)
end

Didn't even need a PR, lovely!

dcarrera · 2016-03-12T19:56:25Z

@dfdx Thanks for the tips! I got master to compile, as well as origin/teh/separate_indexes. It looks like the compile error was in something I did. I reverted all my changes and I'm now going to test them one by one.

dcarrera · 2016-03-12T20:19:10Z

Ok. I made a fairly trivial set of changes to test.jl that compile correctly. The changes are too trivial to improve performance, but it's a step. I pushed the changes up to my repository and created my first ever pull request: #15471

Does this look ok?

dfdx · 2016-03-13T20:35:43Z

I think I can take everything from show.jl and to the end (except for test.jl which is already done by @dcarrera ).

dfdx · 2016-03-13T22:06:44Z

I've got an interesting one here. Here's an implementation of cov2cor! from statistics.jl:

function cov2cor!(C::AbstractMatrix, xsd::AbstractArray, ysd::AbstractArray)
    nx, ny = size(C)
    ... 
    for j = 1:ny
        for i = 1:nx
            C[i,j] /= (xsd[i] * ysd[j])
        end
    end
    return C
end

So we iterate first over columns and then over rows of C, assuming it is efficiently indexed this way. It would be great to be able to write something like:

function cov2cor!(C::AbstractMatrix, xsd::AbstractArray, ysd::AbstractArray)
    ...
    for (i, j) in eachfullindex(C)
        C[i,j] /= (xsd[i] * ysd[j])
    end
    return C
end

Which seems to be both - more readable and efficient. Is there anything in Base we can use instead of eachfullindex now? Would eachindex(A, dims) support this use case?

meggart · 2016-03-14T14:09:29Z

@dfdx I also cam across on of those, so far I just marked it with #Fixme iter As far as I understand it this is exactly what #15459 wants to solve.

Another case I came across was reverse iteration, so what to do about:

for i=length(A):-1:1

Maybe this is something else to consider for the new eachindex API? @timholy

timholy · 2016-03-14T23:05:53Z

You can do this:

for I in CartesianRange(size(C))

However, the more serious issue is that xsd and ysd may not be efficiently indexed by integers. A better option is

for (i, x) in enumerate(xsd)
    for (j, y) in enumerate(ysd)
        C[i,j] /= x*y
    end
end

The only problem with that is what if C has unusual indexing? So I'd suggesting making this change but also marking with a #fixme.

timholy · 2016-03-14T23:06:46Z

@meggart, yeah, mark reversed ones too. This is proving to be a good exercise in finding out what we need to consider in generalizing our iteration!

dfdx · 2016-03-19T10:27:53Z

#15564 covering show.jl...uv_constants.jl is ready for review.

Meanwhile, I notices that make test takes terribly long for me. Is it expected behavior or there's a swithc / option to make things faster?

kswietli · 2016-03-19T11:56:33Z

I intend to do RNG work on Julia in the summer, if someone from GSoC won't work on it instead. Probably a good idea would be to get familiar with the code for those, and I'd like the two sections that include random.jl and dSFMT.jl, if they have not been claimed yet. I may take some others later, but for now I'd like to reserve those two tickboxes to get done over the next few days.

timholy · 2016-03-21T04:35:11Z

Thanks, @dfdx and @kswietli!

meggart · 2016-03-21T13:28:16Z

I would also start another block, so I'd reserve ascii.jl ... collections.jl.

dfdx · 2016-03-26T00:13:30Z

I have some more spare time and can take points from markdown/ to pkg/.

timholy · 2016-03-26T09:41:25Z

Sounds great, thanks!

dfdx · 2016-03-26T15:16:41Z

I see that eachindex isn't defined for SimpleVector, but I can't find a definition of this type and see only one field length in it. Should we define eachindex for SimpleVector or it's strictly internal and better not bother with it?

timholy · 2016-03-27T11:41:24Z

It's strictly internal, so just use an integer.

dfdx · 2016-04-27T21:33:31Z

Wow, this one isn't finished yet. Let's fix it.

I see that there are several items that have been processed already, but aren't marked as done. I also scanned through a couple of other items, so now we can mark everything below the line libgit2/strarray.jl, libgit2/oid.jl, libgit2.jl.

I think I will be able to cover the rest within a week or so and submit one final PR.

timholy · 2016-04-27T23:44:03Z

Done, and thanks in advance! Please do mention this issue in your PR, as back-references will be useful to me if I find there are future changes needed (e.g., https://github.com/timholy/ArrayIteration.jl).

meggart · 2016-04-28T08:12:35Z

Good you revived this thread, I will also submit my next PR for the second block soon.

dfdx · 2016-04-29T22:06:57Z

EDIT: my bad, I forgot about #10523, so I will just leave a comment in the code. Below is my old comment.

While working on floatfuncs.jl I've got an error that I can't understand. The error happens when I change the following lines:

for f in (:trunc,:floor,:ceil,:round)
    @eval begin
        ...
        function ($f){T,R}(::Type{T}, x::AbstractArray{R,2})
            # [ ($f)(T, x[i,j])::T for i = 1:size(x,1), j = 1:size(x,2) ]      # old one
            [ ($f)(T, x[I])::T for I in CartesianRange(size(x)) ]           # new one
        end
        ...
    end
end

With this line changed tests in test/linalg/tridiag.jl fail with such a stack trace:

LoadError: There was an error during testing
 in record at ./test.jl:322
 in do_test at ./test.jl:220
 [inlined code] from ./test.jl:137
 in anonymous at ./<no file>:4294967295
 [inlined code] from ./int.jl:33
 in include_string at ./loading.jl:380
 in include_from_node1 at ./loading.jl:429
 [inlined code] from ./util.jl:179
 in runtests at /home/dfdx/work/repos/julia/test/testdefs.jl:7
 in #16 at /home/dfdx/work/repos/julia/test/runtests.jl:36
 [inlined code] from ./promotion.jl:229
 in #308 at ./multi.jl:1017
 in run_work_thunk at ./multi.jl:747
 [inlined code] from ./multi.jl:1017
 in #307 at ./event.jl:46
while loading /home/dfdx/work/repos/julia/test/linalg/tridiag.jl, in expression starting on line 224

I know that it is either Tridiagonal or SymTridiagonal matrix that causes the error, but I couldn't identify exact case. Does anybody know why CartesianRange may not work for these types of matrices and these functions?

dfdx · 2016-05-04T23:12:49Z

I've just posted a PR (see above) that covers all previously unclaimed lines. Surprisingly, there were only 4 files that needed changes, so I'd say Base is already in a good state anyway.

@timholy Also note that lines from libgit2/ to pkg/ have been processed already (see previous PRs/comments), so feel free to mark "done" too.

timholy · 2016-05-06T13:28:51Z

It's great that we're looking so good.

tkelman reviewed Mar 10, 2016
View reviewed changes

eachindex: use separate indexes for each array

3820cf3

This also replaces several instances of `for i in eachindex(A); a = A[i]` with `for a in A`. The latter is presumably easier for automatic bounds-checking removal.

timholy mentioned this pull request Mar 10, 2016

RFC (C=collaborator): multidimensional comprehensions that use CartesianRange iterators #10523

Closed

timholy force-pushed the teh/separate_indexes branch from fb359e1 to 3820cf3 Compare March 10, 2016 12:14

musm reviewed Mar 10, 2016
View reviewed changes

timholy added a commit that referenced this pull request Mar 11, 2016

Merge pull request #15434 from JuliaLang/teh/separate_indexes

beac4d3

WIP: Disentangle arrays in iteration

timholy merged commit beac4d3 into master Mar 11, 2016

timholy deleted the teh/separate_indexes branch March 11, 2016 16:25

dfdx mentioned this pull request Mar 11, 2016

Disentangled arrays in iteration: repl, serialization and shared arrays #15463

Merged

meggart mentioned this pull request Mar 14, 2016

fix for loop for array indexing #15497

Merged

dfdx mentioned this pull request Mar 19, 2016

More general indexes [2] #15564

Merged

timholy mentioned this pull request Mar 21, 2016

ReshapedArrays, take2 #15449

Merged

timholy mentioned this pull request Mar 26, 2016

RFC: Facilitate not paying double lookup cost when mutating a value in a LinearSlow AbstractArray #15630

Open

dfdx mentioned this pull request Mar 29, 2016

[WIP] More generalindexing [3] #15679

Merged

timholy mentioned this pull request Mar 30, 2016

Julep: solving tricky iteration problems #15648

Closed

vtjnash mentioned this pull request Apr 1, 2016

invalid character index in findfirst #15723

Closed

dfdx mentioned this pull request May 4, 2016

More general indexing [4] #16203

Merged

pabloferz mentioned this pull request Jun 10, 2016

misc. benchmark regressions since 0.4 #16128

Closed

WIP: Disentangle arrays in iteration #15434

WIP: Disentangle arrays in iteration #15434

Conversation

timholy commented Mar 10, 2016

tkelman Mar 10, 2016

Choose a reason for hiding this comment

timholy Mar 10, 2016

Choose a reason for hiding this comment

timholy commented Mar 10, 2016

timholy commented Mar 10, 2016

nalimilan commented Mar 10, 2016

KristofferC commented Mar 10, 2016

eschnett commented Mar 10, 2016

nanosoldier commented Mar 10, 2016

timholy commented Mar 10, 2016

timholy commented Mar 10, 2016 • edited Loading

musm Mar 10, 2016

Choose a reason for hiding this comment

timholy Mar 10, 2016

Choose a reason for hiding this comment

musm commented Mar 10, 2016

timholy commented Mar 10, 2016

dcarrera commented Mar 11, 2016

KristofferC commented Mar 11, 2016

meggart commented Mar 11, 2016

timholy commented Mar 11, 2016

timholy commented Mar 11, 2016

dfdx commented Mar 11, 2016

timholy commented Mar 11, 2016

meggart commented Mar 11, 2016

meggart commented Mar 11, 2016

pkofod commented Mar 11, 2016

dcarrera commented Mar 12, 2016

dcarrera commented Mar 12, 2016

dfdx commented Mar 13, 2016

dfdx commented Mar 13, 2016

meggart commented Mar 14, 2016

timholy commented Mar 14, 2016

timholy commented Mar 14, 2016

dfdx commented Mar 19, 2016

kswietli commented Mar 19, 2016

timholy commented Mar 21, 2016

meggart commented Mar 21, 2016

dfdx commented Mar 26, 2016

timholy commented Mar 26, 2016

dfdx commented Mar 26, 2016

timholy commented Mar 27, 2016

dfdx commented Apr 27, 2016

timholy commented Apr 27, 2016

meggart commented Apr 28, 2016

dfdx commented Apr 29, 2016 • edited Loading

dfdx commented May 4, 2016

timholy commented May 6, 2016

timholy commented Mar 10, 2016 •

edited

Loading

dfdx commented Apr 29, 2016 •

edited

Loading