{"action":"created","issue":{"url":"https://api.github.com/repos/JuliaLang/julia/issues/18777","repository_url":"https://api.github.com/repos/JuliaLang/julia","labels_url":"https://api.github.com/repos/JuliaLang/julia/issues/18777/labels{/name}","comments_url":"https://api.github.com/repos/JuliaLang/julia/issues/18777/comments","events_url":"https://api.github.com/repos/JuliaLang/julia/issues/18777/events","html_url":"https://github.com/JuliaLang/julia/pull/18777","id":180684412,"number":18777,"title":"Make LinSpace generic","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,"assignees":[],"milestone":null,"comments":28,"created_at":"2016-10-03T16:25:43Z","updated_at":"2016-10-06T11:36:40Z","closed_at":null,"pull_request":{"url":"https://api.github.com/repos/JuliaLang/julia/pulls/18777","html_url":"https://github.com/JuliaLang/julia/pull/18777","diff_url":"https://github.com/JuliaLang/julia/pull/18777.diff","patch_url":"https://github.com/JuliaLang/julia/pull/18777.patch"},"body":"This is my proposed alternative to #18492 for fixing the issues with LinSpace. The main attractions of this PR are:\r\n- closes #14420 \r\n- you can now use `linspace` with any type `T` for which `(1-t)*a::T + t*b::T` makes sense.\r\n\r\nGiven the large amount of work & discussion that has gone into `linspace` and friends with respect to numerical precision (at least #2333, #9637, #18492), it would be remiss not to comment (extensively :wink:) on how this does and does not address previous concerns or implementations. Let me address the issues one-by-one.\r\n\r\n## Mind-reading\r\n\r\nFirst and foremost, this implementation does not try to read your mind, whereas the implementation currently in master does. The easiest way to see what I mean is to compare master's behavior:\r\n```jl\r\njulia> r = linspace(0.0, 0.3, 4)\r\n4-element LinSpace{Float64}:\r\n 0.0,0.1,0.2,0.3\r\n\r\njulia> r[2]\r\n0.1\r\n```\r\nwith the behavior of this PR:\r\n```jl\r\njulia> r = linspace(0, 0.3, 4)\r\n4-element LinSpace{Float64}:\r\n 0.0,0.1,0.2,0.3\r\n\r\njulia> r[2]\r\n0.09999999999999999\r\n```\r\nNote that, **given the specified endpoints**, that's the correct answer:\r\n```jl\r\njulia> 0.3/3\r\n0.09999999999999999\r\n\r\njulia> Float64(BigFloat(0.3)/3)\r\n0.09999999999999999\r\n```\r\nThe version in master returns 0.1 essentially because it interprets the right endpoint 0.3 as `3//10`. This is what I refer to as mind-reading, and while it has some undeniable attractions, it's (1) only applicable to special types (e.g., Float32/64), (2) fairly fragile (it's the origin of #14420), and (3) arguably wrong even in principle (the user said 0.3, not `3//10`).\r\n\r\nNote that, now, people who want mind reading can be told to use `Rational`s directly:\r\n```jl\r\njulia> r = linspace(0, 3//10, 4)\r\n4-element LinSpace{Rational{Int64}}:\r\n 0//1,1//10,1//5,3//10\r\n\r\njulia> r[2]\r\n1//10\r\n```\r\n\r\n## Precision\r\n\r\nAll that said, there's still a valid concern about precision: linear interpolation is famously subject to roundoff error. **This implementation returns results with full numeric precision** by using two calls to `fma`. The key function (which other types can specialize) is:\r\n```jl\r\n# High-precision interpolation. Accurate for t ∈ [0.5,1], so that 1-t is exact.\r\nfunction lerpi{T}(j::Integer, d::Integer, a::T, b::T)\r\n    t = j/d\r\n    # computes (1-t)*a + t*b\r\n    T(fma(t, b, fma(-t, a, a)))\r\nend\r\n```\r\nAs long as `j >= d/2`, it turns out you don't need to compute `t` to higher precision for this to return a numerically-precise result (as long as you use `fma`). \r\n\r\nAll the (new) tests compare the result against what happens when you take the two endpoints, convert them to BigFloats, perform all interpolation at BigFloat precision, and then cast back to the element type. That's what I mean by \"full precision\" results---there is no roundoff error to within the precision of the numerical type, as determined by the inputs upon construction.\r\n\r\n## Test changes and deletions\r\n\r\nFirst, I have to say that whoever wrote the old tests (mostly Stefan?) did an amazing job. While my final implementation is dirt-simple, I confess it was not the first (or second, or third...) thing I tried; it was very demanding to find an implementation that could pass the \"good\" ones (spelled out below). The extensive tests saved me from many implementations that seemed equivalent but ran into trouble in corner cases. :100: \r\n\r\nThat said, anyone who looks at this PR will see many test changes and a few deletions, because the new implementation does *not* pass many of our old tests. For the vast majority of these I am unconcerned and changed them guilt-free, because most of these turned out to rely on mind-reading: in particular, they don't pass the \"convert endpoints to BigFloat before performing operations\" test. Even the so-called \"additive identity\" tests (e.g., `(r - 1) + 1`) fall in that category when combined with `collect`, due to the non-associativity of floating-point addition. (The old tests happened to pass for the provided inputs, but it doesn't take long at all to discover that tweaking the endpoints in those tests, in just about any way, breaks them. So we were testing for an \"identity\" that doesn't actually hold in general.)\r\n\r\nHowever, a few of the test deletions make me sad or are worthy of public confession:\r\n- Behavior at `u = eps(0.0)`: most tests pass, but ones that specifically test for generation of `-0.0` don't. I'm basically unbothered by this---we have several concerns about `-0.0` anyway (e.g., #18341). If you think of it as a bug, it's probably best thought of as an `fma` bug, not a problem with this implementation of `linspace`.\r\n- `linspace(-3u,3u,7)` doesn't pass. I'd rather this one passed, because it's seemingly unrelated to nonsense like `-0.0`. I suspect this is an `fma` bug, so I kept the test and marked it `@test_skip`. Note this test passes if `u = 2*eps(0.0)`, so the failure is specific to the smallest-magnitude floating point number.\r\n\r\n## Performance\r\n\r\nI'm hoping that on machines with CPU-support for `fma`, this will be competitive with our current implementation. `fma` support was introduced around 2012-2013, so there many be a good number of machines that don't have such support. Without it, this ends up `ccalling` C code for the `fma` and is therefore glacially slow.\r\n\r\nIf we are concerned about performance on older machines, one potential way to fix it would be to [Dekker-split](https://github.com/simonbyrne/DoubleDouble.jl) the two endpoints upon construction; then you'd only have to split `t` at runtime, and perhaps along with some tweaks that take advantage of the symmetry of the problem you might be able to largely close the gap.\r\n\r\nWe might also consider adding a `@fastmath` specialization.\r\n\r\nBut before getting too speculative about performance, let's see what `@nanosoldier runbenchmarks(ALL, vs=:master)` thinks.\r\n\r\n## Future changes to Ranges\r\n\r\nOne source of concern is that `LinSpace` and `FloatRange` now have different behavior. I think the right way to handle this is to cause `start:step:stop` to return a `LinSpace` when any of those entries are `AbstractFloat`s, and move `FloatRange` to a package. That way there won't be any inconsistencies. `StepRange` can still exist, only now (with this more generic `LinSpace`) I think it can safely be confined to `Integer` (thus fixing all kinds of problems like those I'm attempting to address in  #18744).\r\n"},"comment":{"url":"https://api.github.com/repos/JuliaLang/julia/issues/comments/251937467","html_url":"https://github.com/JuliaLang/julia/pull/18777#issuecomment-251937467","issue_url":"https://api.github.com/repos/JuliaLang/julia/issues/18777","id":251937467,"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-10-06T11:36:40Z","updated_at":"2016-10-06T11:36:40Z","body":"@c42f, regarding integer linspaces, I think the only thing that's confusing currently is `Base.print_ranges`, other than that I think it already does what you want:\r\n```jl\r\njulia> using Ranges\r\n\r\njulia> r = Ranges.LinSpace{Int}(1,3,3)\r\n3-element Ranges.LinSpace{Int64}:\r\n 1.0,2.0,3.0                                             # this is confusing, should fix\r\n\r\njulia> r[1]\r\n1\r\n\r\njulia> typeof(r[1])\r\nInt64\r\n\r\njulia> r[2]\r\n2\r\n\r\njulia> r[3]\r\n3\r\n\r\njulia> r = Ranges.LinSpace{Int}(1,7,3)\r\n3-element Ranges.LinSpace{Int64}:\r\n 1.0,4.0,7.0\r\n\r\njulia> r[2]\r\n4\r\n\r\njulia> r = Ranges.LinSpace(1,10,3)\r\n3-element Ranges.LinSpace{Float64}:\r\n 1.0,5.5,10.0\r\n\r\njulia> r[2]\r\n5.5\r\n\r\njulia> r = Ranges.LinSpace{Int}(1,10,3)\r\n3-element Ranges.LinSpace{Int64}:\r\n 1.0,5.5,10.0\r\n\r\njulia> r[1]\r\n1\r\n\r\njulia> r[2]\r\nERROR: InexactError()\r\n in unsafe_getindex(::Ranges.LinSpace{Int64}, ::Int64) at /home/tim/.julia/v0.5/Ranges/src/linspace.jl:76\r\n in getindex(::Ranges.LinSpace{Int64}, ::Int64) at /home/tim/.julia/v0.5/Ranges/src/linspace.jl:69\r\n```\r\n"}}
{"actor":{"display_login":"timholy"}}