{"action":"created","issue":{"url":"https://api.github.com/repos/JuliaLang/julia/issues/15648","repository_url":"https://api.github.com/repos/JuliaLang/julia","labels_url":"https://api.github.com/repos/JuliaLang/julia/issues/15648/labels{/name}","comments_url":"https://api.github.com/repos/JuliaLang/julia/issues/15648/comments","events_url":"https://api.github.com/repos/JuliaLang/julia/issues/15648/events","html_url":"https://github.com/JuliaLang/julia/issues/15648","id":143882891,"number":15648,"title":"Julep: solving tricky iteration problems","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":20,"created_at":"2016-03-28T03:21:02Z","updated_at":"2016-03-30T13:42:43Z","closed_at":null,"body":"# The problem (a brief summary)\r\n\r\nIn a [PR for a blog post](https://github.com/JuliaLang/julialang.github.com/pull/366), I described a number of challenges for several potential new `AbstractArray` types.  Very briefly, the types and the challenges they exhibit are:\r\n\r\n- `ReshapedArrays`: the only efficient way of indexing some types is by indexing the parent array directly (i.e., with an iterator optimized for the pre-reshaping array).  This means that indexing a 2D array as `B[i,j]`, for integer `i` and `j`, may not always be efficient.\r\n\r\n- `TransposedMatrix` (for #4774): iteration is most efficient when it produces a cache-friendly sequence of array accesses. Unfortunately, algorithms that mix column-major and row-major arrays imply that we can't naively allow `eachindex(A)` to always return the most efficient iterator.\r\n\r\n- `OffsetArrays` (sometimes called \"Fortran arrays\"): these violate the assumptions that array indexes run from `1:size(A,d)`, and more generally question how one should match elements of two different\r\n  arrays.\r\n\r\n# Possible APIs\r\n\r\nIn the blog post, only two concrete algorithms were considered: matrix-vector multiplication and `copy!`.  At the risk of being insufficiently general, let's explore possible APIs in terms of these two functions.\r\n\r\n## The \"raw\" API\r\n\r\nFor the moment, let's worry about correctness and efficiency more than elegance.  Correctness trumps all other considerations.  Efficiency means that we have to be able to iterate over the matrix in\r\ncache-friendly order.  That means that the array `B` has to be \"in charge\" of the sequence of values and pick its most efficient iterator.\r\n\r\nI should emphasize that I haven't taken even the first step towards attempting to implement this API.  So this is crystal-ball development, trying to imagine what might work while hoping to shake\r\nout likely problems in advance via community feedback.  That caveat understood, here's a possible implementation of matrix-vector multiplication:\r\n\r\n```jl\r\nRB = eachindex(B)\r\nRv = eachindex(RB[*,:], v)\r\nRdest = eachindex(RB[:,*], dest)\r\nfor (idest, iB, iv) in zip(Rdest, RB, Rv)\r\n    dest[idest] += B[iB] * v[iv]\r\nend\r\n```\r\n\r\n`eachindex(iter, obj)` should perform bounds-checking, to ensure that the indexes of `iter` align with `obj` (even handling things like `OffsetArrays`).  The notation `RB[*,:]` means \"corresponding to the\r\nsecond dimension of `RB`\"; interestingly, this expression parses without error, attempting to dispatch to `getindex` with argument types `Tuple{typeof(RB), Function, Colon}`. Hopefully there isn't another pending meaning for indexing with a `Function`.\r\n\r\nBecause either `idest` or `iv` is constant over the \"inner\" loop (depending on storage order of `B`), ideally one would like to hoist the corresponding access `dest[idest]` or `v[iv]`.  I'd argue that's a\r\njob for the compiler, and not something we should have to bake into the API.\r\n\r\nBecause `zip(X, Y)` does not allow `Y` to see the state of `X`, it seems quite possible that this could not be made efficient in general (consider the case described [for `ReshapedArray`s](https://github.com/JuliaLang/julialang.github.com/blob/teh/arrays/blog/_posts/2016-03-27-arrays-iteration.md#reshapedarrays-15449), where the\r\nefficiently-parametrized iterator for the first dimension depends on the state of the index for the second dimension).  If this proves to be the case, we have to abandon `zip` and write this more like\r\n\r\n```jl\r\nRB = eachindex(B)\r\nfor (iB, idest, iv) in couple(RB, (RB[:,*], dest), (RB[*,:], v))\r\n    dest[idest] += B[iB] * v[iv]\r\nend\r\n```\r\n\r\nwhere `couple` is a new exported function.\r\n\r\nOne interesting thing about this API is that sparse multiplication might be implemented efficiently just by replacing the first line with\r\n\r\n```jl\r\nRB = eachstored(B)\r\n```\r\n\r\nso that only \"stored\" elements of `B` are visited.  (@tkelman may have been the first to propose an `eachstored` iterator.)  It should be pointed out that several recent discussions have addressed some of the complications that arise from `Inf` and `NaN`, and it remains to be determined whether a version that [doesn't ignore these complications](https://github.com/JuliaLang/julia/pull/15579#issuecomment-200604174) can be written with this (or similar) API.\r\n\r\nLikewise, `copy!` might be implemented\r\n```jl\r\nRdest = eachindex(dest)\r\nfor (idest, isrc) in couple(Rdest, (Rdest, src))\r\n    dest[idest] = src[isrc]\r\nend\r\n```\r\n\r\nAt least for the array types described here, this API seems plausibly capable of generating both correct and (fairly) efficient code.  One might ideally want to write generic \"tiled\" implementations, but that seems beyond the scope of this initial effort.\r\n\r\n## Making it prettier with \"sufficiently smart iteration\"\r\n\r\nThis API is still more complicated than naive iteration.  One might possibly be able to hide a lot of this with a magical `@iterate` macro,\r\n\r\n```jl\r\n@iterate B dest[i] += B[i,j]*v[j]\r\n```\r\n\r\nfor matrix-vector multiplication and\r\n\r\n```jl\r\n@iterate dest dest[i] = src[i]\r\n```\r\n\r\nfor `copy!`. The idea is that the `@iterate` macro would expand these expressions into their longer forms above. Note that the first argument to `@iterate` specifies which object is \"in charge\" of the\r\nsequence of operations.  `@iterate` might be too limited for many things, but may be worth exploring. The [KernelTools](https://github.com/timholy/KernelTools.jl) repository may have some useful tidbits.\r\n"},"comment":{"url":"https://api.github.com/repos/JuliaLang/julia/issues/comments/203439381","html_url":"https://github.com/JuliaLang/julia/issues/15648#issuecomment-203439381","issue_url":"https://api.github.com/repos/JuliaLang/julia/issues/15648","id":203439381,"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},"created_at":"2016-03-30T13:42:43Z","updated_at":"2016-03-30T13:42:43Z","body":"Thanks to #13412, I'm noticing that this is also possible: `couple(any, x, y)` or `couple(all, x, y)`. I'm not sure you'd really use the `any` and `all` functions directly, but you can dispatch on them:\r\n```jl\r\n@eval foo(::$(typeof(any))) = 1\r\n@eval foo(::$(typeof(all))) = 2\r\n```\r\n"}}