{"action":"created","issue":{"url":"https://api.github.com/repos/JuliaLang/julia/issues/16878","repository_url":"https://api.github.com/repos/JuliaLang/julia","labels_url":"https://api.github.com/repos/JuliaLang/julia/issues/16878/labels{/name}","comments_url":"https://api.github.com/repos/JuliaLang/julia/issues/16878/comments","events_url":"https://api.github.com/repos/JuliaLang/julia/issues/16878/events","html_url":"https://github.com/JuliaLang/julia/issues/16878","id":159769638,"number":16878,"title":"julep: extended iteration API, with proof-of-principle for fixing performance of Cartesian iteration","user":{"login":"timholy","id":1525481,"avatar_url":"https://avatars.githubusercontent.com/u/1525481?v=3","gravatar_id":"","url":"https://api.github.com/users/timholy","html_url":"https://github.com/timholy","followers_url":"https://api.github.com/users/timholy/followers","following_url":"https://api.github.com/users/timholy/following{/other_user}","gists_url":"https://api.github.com/users/timholy/gists{/gist_id}","starred_url":"https://api.github.com/users/timholy/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/timholy/subscriptions","organizations_url":"https://api.github.com/users/timholy/orgs","repos_url":"https://api.github.com/users/timholy/repos","events_url":"https://api.github.com/users/timholy/events{/privacy}","received_events_url":"https://api.github.com/users/timholy/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"milestone":null,"comments":3,"created_at":"2016-06-11T11:38:20Z","updated_at":"2016-06-11T18:39:23Z","closed_at":null,"body":"This proposes a change in our iteration API (refs #9182, #9178, #6125) that narrows the gap on #9080 (which seems to have gotten worse), following the convoluted process of discovery in #16035. Let's suppose this function:\r\n```jl\r\nfunction mysum(A)\r\n    s = zero(eltype(A))\r\n    for I in eachindex(A)\r\n        @inbounds s += A[I]\r\n    end\r\n    s\r\nend\r\n```\r\ngot expanded to something like this (note: I'm passing in the iterator as an argument so I can compare performance of existing and new implementations):\r\n```jl\r\nfunction itergeneral(A, iter)\r\n    # We'll sum over the elements of A that lie within bounds of iter\r\n    s = zero(eltype(A))\r\n    state = start(iter)\r\n    isdone = maybe_done(iter, state)\r\n    if isdone\r\n        isdone, state = done(iter, state, true)  # use dummy 3rd arg to indicate 2 returns\r\n    end\r\n    if !isdone\r\n        while true\r\n            item, state = next(iter, state)\r\n            @inbounds s += A[item]\r\n            if maybe_done(iter, state)\r\n                isdone, state = done(iter, state, true)\r\n                isdone && break\r\n            end\r\n        end\r\n    end\r\n    s\r\nend\r\n```\r\n\r\nTo avoid breakage of current code, we need the following fallback definitions:\r\n```jl\r\nimport Base: start, next, done, getindex\r\n\r\n# fallback definitions\r\nmaybe_done(iter, state) = true\r\n# For 3-argument code when only the 2-argument version is defined\r\ndone(iter, state, ::Bool) = done(iter, state), state\r\n# For explicit 2-argument code when only the 3-argument version is defined\r\ndone(iter, state) = done(iter, state, true)[1]\r\n```\r\nThe only downside of this I see is that if no valid definition of `done` exists, now you get a `StackOverflow` rather than `MethodError`. Aside from that, this appears to be non-breaking. If the compiler knows that `maybe_done` always returns `true` for `iter`, then note that the code above simplifies to\r\n```jl\r\nfunction itergeneral(A, iter)\r\n    # We'll sum over the elements of A that lie within bounds of iter\r\n    s = zero(eltype(A))\r\n    state = start(iter)\r\n    isdone, state = done(iter, state, true)  # use dummy 3rd arg to indicate 2 returns\r\n    if !isdone\r\n        while true\r\n            item, state = next(iter, state)\r\n            @inbounds s += A[item]\r\n            isdone, state = done(iter, state, true)\r\n            isdone && break\r\n        end\r\n    end\r\n    s\r\nend\r\n```\r\n\r\nAn advantage of this more symmetric interface is that either `next` or `done` can increment the `state`; therefore, this subsumes the `nextval`/`nextstate` split first proposed in #6125. But unlike the implementation in #9182, thanks to `maybe_done` this also introduces a key advantage for #9080: it ensures that typically (when dimension 1 has size larger than 1) there's only one branch per iteration, using two or more branches only when the \"carry\" operation needs to be performed.\r\n\r\nProof that this helps #9080 (with the previous definitions all loaded):\r\n```jl\r\nusing Base: tail, @propagate_inbounds\r\n\r\nimmutable CartIdx{N}\r\n    I::NTuple{N,Int}\r\n    CartIdx(index::NTuple{N,Integer}) = new(index)\r\nend\r\nCartIdx{N}(index::NTuple{N,Integer}) = CartIdx{N}(index)\r\ngetindex(ci::CartIdx, d::Integer) = ci.I[d]\r\n@propagate_inbounds getindex(A::AbstractArray, ci::CartIdx{2}) = A[ci.I[1], ci.I[2]]\r\n\r\nimmutable CartRange{I}\r\n    start::I\r\n    stop::I\r\nend\r\n\r\nstart(cr::CartRange) = cr.start\r\nnext(cr::CartRange, state::CartIdx) = state, inc1(state)\r\n@inline inc1(state) = CartIdx((state[1]+1, tail(state.I)...))  # a \"dumb\" incrementer\r\nmaybe_done(cr::CartRange, state::CartIdx) = state[1] > cr.stop[1]\r\nfunction done{N}(cr::CartRange, state::CartIdx{N}, ::Bool)\r\n    # handle the carry operation (clean up after the \"dumb\" incrementer)\r\n    state = CartIdx(carry((), state.I, cr.start.I, cr.stop.I))\r\n    state[N] > cr.stop[N], state\r\nend\r\ncarry(out, ::Tuple{}, ::Tuple{}, ::Tuple{}) = out\r\n@inline carry(out, s::Tuple{Int}, start::Tuple{Int}, stop::Tuple{Int}) = (out..., s[1])\r\n@inline function carry(out, s, start, stop)\r\n    if s[1] > stop[1]\r\n        return carry((out..., start[1]), (s[2]+1, tail(tail(s))...), tail(start), tail(stop))\r\n    end\r\n    return s\r\nend\r\n```\r\n\r\nwith the following test script (`sumcart_manual` and `sumcart_iter` are from #9080):\r\n```jl\r\nA = rand(10^4, 10^4)\r\nsumcart_manual(A)\r\nsumcart_iter(A)\r\niter = CartRange(CartIdx((1,1)), CartIdx(size(A)))\r\nR = CartesianRange(size(A))\r\nitergeneral(A, R)\r\nitergeneral(A, iter)\r\n@time 1\r\n@time sumcart_manual(A)\r\n@time sumcart_iter(A)\r\n@time itergeneral(A, R)\r\n@time itergeneral(A, iter)\r\n```\r\nwith results\r\n```jl\r\n  0.000004 seconds (150 allocations: 9.106 KB)\r\n  0.108923 seconds (5 allocations: 176 bytes)\r\n  0.589222 seconds (5 allocations: 176 bytes)\r\n  0.950276 seconds (5 allocations: 176 bytes)\r\n  0.220657 seconds (5 allocations: 176 bytes)\r\n```\r\nIt's awesome that this is almost 3x faster than our current scheme. It stinks that it's still not as good as the manual version. I profiled it, and almost all the time is spend on `@inbounds s += A[item]`. Other than the possibility that somehow the CPU isn't as good at cache prefetch with this code as with traditional for-loop code, I'm at a loss to explain the gap. \r\n\r\nNote that currently we have one other way of circumventing #9080: add `@simd`. Even when this doesn't actually turn on vectorization, the `@simd` macro splits out the cartesian iterator into \"inner index\" and \"remaining index,\" and thus achieves parity with the manual version (i.e., is better than this julep) even without vectorization. Now, `@simd` is limited, so this julep is still attractive, but I would feel stronger about pushing for it if we could achieve parity.\r\n"},"comment":{"url":"https://api.github.com/repos/JuliaLang/julia/issues/comments/225382476","html_url":"https://github.com/JuliaLang/julia/issues/16878#issuecomment-225382476","issue_url":"https://api.github.com/repos/JuliaLang/julia/issues/16878","id":225382476,"user":{"login":"mbauman","id":154641,"avatar_url":"https://avatars.githubusercontent.com/u/154641?v=3","gravatar_id":"","url":"https://api.github.com/users/mbauman","html_url":"https://github.com/mbauman","followers_url":"https://api.github.com/users/mbauman/followers","following_url":"https://api.github.com/users/mbauman/following{/other_user}","gists_url":"https://api.github.com/users/mbauman/gists{/gist_id}","starred_url":"https://api.github.com/users/mbauman/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/mbauman/subscriptions","organizations_url":"https://api.github.com/users/mbauman/orgs","repos_url":"https://api.github.com/users/mbauman/repos","events_url":"https://api.github.com/users/mbauman/events{/privacy}","received_events_url":"https://api.github.com/users/mbauman/received_events","type":"User","site_admin":false},"created_at":"2016-06-11T18:39:22Z","updated_at":"2016-06-11T18:39:22Z","body":">  I ended up writing recursive code, where an n-dimensional array is traversed via n nested function calls; after inlining, you recover the usual n nested loops.\r\n\r\nI tested this a long time ago with non-scalar indexing.  It was quite a bit slower back then, but I've not tried it since.  https://gist.github.com/mbauman/6434fff5f793cf0414559fa554822578"}}
{"actor":{"display_login":"mbauman"}}