float ranges: if a+n*s == b, length(a:s:b) should be n+1

[StefanKarpinski]

See http://stackoverflow.com/questions/41985718/why-the-values-are-not-correctly-iterated-using-range-arrays/41986875#41986875.

Example:

julia> a, s, b = 3*0.05, 0.05, 4*0.05
(0.15000000000000002,0.05,0.2)

julia> a + 1*s == b
true

julia> length(a:s:b)
1

Rational lifting doesn't address this case since 0.15000000000000002 != 3/20 == 0.15.

[timholy]

Interestingly, this does work with #18777 (comment). That might introduce problems of its own, of course. A more conservative version is

function mycolon(a, s, b)
    len_est = (b-a)/s + 1
    len_test = floor(Int, len_est+10*eps(len_est))
    len = a+(len_test-1)*s <= b ? len_test : len_test-1  # needs modification for s < 0
    StepRangeLen(a, s, len)
end

[StefanKarpinski]

I can't find the comment that links to – there's way too many. Which proposal fixes it?

[StefanKarpinski]

Why not this:

n = ceil(Int, (b-a)/s)
a + n*s <= b || (n -= 1)
StepRangeLen(a, s, n+1)

In case you haven't figured it out, I'm vehemently against arbitrary ulp fuzzing. Why is 10 the correct number of ulps?

[JeffBezanson]

For this case would it be enough to check whether adding 1 to the length exactly hits b, no eps needed?

[StefanKarpinski]

Agree – that's all that's needed, I think.

[timholy]

The problem comes if you've already committed to rational approximation of the inputs---at that point you may need some ulp fuzzing.

[timholy]

I don't disagree with your vehement distaste for arbitrary constants, however.

[StefanKarpinski]

If you've already committed to rational approximation, then you do that – that guarantees hitting the endpoint exactly or you would bail out. This is a slightly better way of picking the length in the case where rational lifting has already failed.

[timholy]

Yeah, you're probably right. Seems like it should work, and it would be a very simple tweak. Are you trying it, or should I?

[StefanKarpinski]

I'm trying it already. In other words, just this:

diff --git a/base/twiceprecision.jl b/base/twiceprecision.jl
index 80d1917ea..05c7d1153 100644
--- a/base/twiceprecision.jl
+++ b/base/twiceprecision.jl
@@ -149,6 +149,7 @@ function colon{T<:Union{Float16,Float32,Float64}}(start::T, step::T, stop::T)
         end
     end
     # Fallback, taking start and step literally
+    start + len*step <= stop && (len += 1)
     StepRangeLen(TwicePrecision(start, zero(T)), twiceprecision(step, nbitslen(T, len, 1)), len)
 end

Not sure where else I'd need to apply this, though.

[timholy]

That's probably it. I think the only other place with similar logic is for range and linspace, but there len is provided directly.

[StefanKarpinski]

I also just noticed that the definition of len at the top is dead code...

[timholy]

? If you're looking at the same function, here's a resetting of len later but you only get there if early steps of rational approximation succeed. It would be clearer if moved to the end, however, and might indeed avoid corner cases where the re-computation of len overflows...so that would be a strict improvement.

[StefanKarpinski]

So there's good news and bad news. The good news is that I fixed this problem. The bad news is that in writing tests for it, I found that there's a fair number of start, step, stop with n = round(Int, (stop-start)/step) such that one or more of the following fail:

@test first(start:step:stop) == start
@test last(start:step:stop) == stop
@test first(linspace(start,stop,n)) == start
@test last(linspace(start,stop,n)) == stop

The test cases are randomly generated, so not terribly typical and it does take a lot of randomly generated cases to hit a bad one, but this probably should not happen. The tests that sometimes fail are included as @test_skip.

[JeffBezanson]

last(start:step:stop) == stop

I think this being false is expected behavior, unavoidable given the interface.

With linspace though, my understanding is that part of the point of it is to hit start and stop exactly.

[StefanKarpinski]

In this case, I've constructed the values so that start + (len-1)*step == stop, so it's definitely not expected behavior, it's a bug.

[timholy]

It can't possibly miss the endpoints if you use LinSpace(start, stop, len), but to implement mind-reading the lowercase version tries to be more clever (and obviously misses sometimes). The tricky part of the logic is here.

[StefanKarpinski]

There are failures for the linspace case as well.

[StefanKarpinski]

Oh, I see, I missed "lowercase version" in there.

[timholy]

Some numbers (obtained by making the obvious modifications to Stefan's new test script):
Master (with #20377):

Fraction that fail:
colon:
  len 0.0, start 6.103515625e-6, stop 0.0545379638671875
linspace:
  len 0.0, start 2.13623046875e-5, stop 1.52587890625e-5

julia-0.5:

Fraction that fail:
colon:
  len 0.085711669921875, start 0.0284698486328125, stop 0.085711669921875
linspace:
  len 0.0, start 0.02374267578125, stop 0.028363037109375

So we're better off in every respect, typically by several orders of magnitude, but still not perfect.

[StefanKarpinski]

Definitely not a show-stopper, but dayum, this is a hard problem!

[timholy]

I'm pretty sure the framework is capable of hitting all of these cases exactly. There's an exception for really long ranges, because of this which deliberately breaks @simonbyrne's criterion

step_hi contains sufficient trailing zeros such that step_hi * (i - offset) is exact for all 1 <= i <= length

It turns out to be counterproductive to truncate more than half the bits of the mantissa of step_hi.

But that doesn't apply here. It's presumably "just" a matter of tweaking those TwicePrecision constants. The mathematics is very clear (it's just 2 linear equations with 2 unknowns), but the handling of twice-precision arithmetic and nonassociativity makes it a little more "exciting." Fortunately, you're good at that kind of thing, so I'm sure you'll pull it off 😉.

[StefanKarpinski]

See #20396 – with that patch the 0.5 float ranges always hit the start, length and endpoint exactly.

[timholy]

In this case, I've constructed the values so that start + (len-1)*step == stop, so it's definitely not expected behavior, it's a bug.

I don't think that's true. For example, my impression is that you think it's desirable that currently we have last(0.0:0.1:0.3) == 0.3 despite the fact that stp = stp = 3*0.1 == 0.30000000000000004. If we went with your claim here, then 0.0:0.1:0.3 would return a different range from 0.0:0.1:stp. What would 0.0:0.1:0.31 return?

[StefanKarpinski]

IIRC (it's been a while), I was excluding such cases and only testing cases where the end point is hit exactly by the start + (len-1)*step. Can you give me a little more context on this necrocomment?

[timholy]

Sure, but what should collect(0.1:0.1:0.30000000000000004) return? The answer has to be [0.0, 0.1, 0.2, 0.3] or we're being inconsistent, and hence last(r) != stop even though stop = start + (len-1)*step.

[StefanKarpinski]

I don't get it – the last value of 0.1:0.1:0.30000000000000004 should be 0.30000000000000004 which it currently is. This is exactly what the title says:

julia> a, s, b, n = 0.1, 0.1, 0.30000000000000004, 2
(0.1, 0.1, 0.30000000000000004, 2)

julia> a + n*s == b
true

julia> length(a:s:b) == n+1
true

Perhaps I'm being dense here and not seeing what the problem is. Can you give me a problematic example?

[timholy]

We have

julia> 0.1:0.1:0.30000000000000004
0.1:0.1:0.30000000000000004                                                                                                                                   
                                                                                                                                                                            
julia> 0.1:0.1:0.3
0.1:0.1:0.3                                                                                                                                                                 
                                                                                                                                                                            
julia> 0.1:0.1:0.31
0.1:0.1:0.3                                                                                                                                                                                  

but this is fairly crazy and definitely not defensible. In #23194 you sometimes get this and sometimes don't hit the endpoint exactly, depending on whether either variant hits the rational approximation criterion.

The point being that the test you devised is too stringent to be sensible with non-associative arithmetic. We have to drop the stop criterion with the colon constructor.

[StefanKarpinski]

These seems fine to me... What's the principle that it violates?

[StefanKarpinski]

Ah, I do see what you're getting at here. The issue is that the stop value 0.31 is larger than 0.30000000000000004 yet ends up with a smaller actual stop value in the resulting range. That is rather unintuitive. It seems to me that the result of 0.1:0.1:0.31 is not really right – if the algorithm is going to consider the stop point to be inexact, then it needs to consider the start and step values to be exact since the overspecification that justifies rational "lifting" is that there are three values (start, step, stop) that can't all be simultaneously exact. But as soon as the stop point is inexact, that justification goes away and the start and step values must be taken literally.

[StefanKarpinski]

There used to be a check at the very top of the lifting code that made sure that start + round((stop-start)/step))*step ≈ stop or something to that effect, but it got removed during the great refactoring. That should be put back or replaced by a test somewhere else to the same effect. The bottom line is that if the stop point cannot be hit exactly, then the start and stop points must be taken literally.

[StefanKarpinski]

I keep looking at this code trying to figure out where that check should go, but I don't really understand the DoublePrecision stuff well enough. I suspect that this inner check should be way more stringent:

# Integer ops could overflow, so check that this makes sense
if isbetween(start, start + (len-1)*step, stop + step/2) &&
  !isbetween(start, start + len*step, stop)

[timholy]

From the standpoint of exact arithmetic (i.e., math, not floating-point approximations of it), start, step, and length specify each value of the range of exactly; moreover, many different stop values will generate the same range. Consequently, colon should become a convenience constructor for range and not something that has its own unique set of behaviors. That's what's behind this comment.

Consequently, all 3 of these should produce 0.1:0.1:0.3, and you just have to live with not hitting the endpoint exactly sometimes even though it is equal (in inexact floating-point arithmetic) to start + (len-1)*step.

[StefanKarpinski]

I don't disagree with that comment, but r = start:Δ:stop need not have last(r) == stop. However, I definitely do not think that they should all produce 0.1:0.1:0.3 – I still maintain that if start + n*step == stop then we should always give a range that ends at stop exactly. Otherwise we are taking unacceptable liberties with the meaning of floating-point numbers and are no longer "lifting" in any sense, but just making things up. On the other hand, I do think the current behavior is incorrect in that 0.1:0.1:0.31 should end at 0.30000000000000004 == 0.1 + 2*0.1 exactly.

[timholy]

So then

julia> collect(0.0:0.1:1.0)
11-element Array{Float64,1}:
 0.0
 0.1
 0.2
 0.3
 0.4
 0.5
 0.6
 0.7
 0.8
 0.9
 1.0

needs to switch to

 0.0
 0.1
 0.2
 0.30000000000000004
...

? Because of course 0.0:0.1:0.31 should produce a subset of 0.0:0.1:1.0.

[StefanKarpinski]

Because of course 0.0:0.1:0.31 should naturally produce a subset of 0.0:0.1:1.0.

Ah, here's the disagreement. I don't think that's true at all. We can think of the a:s:b syntax as meaning two very different things balled into one: the exact range construct and the inexact range construct where the endpoint is inexact but gives an upper bound for the range values. The rational lifting is only justified in the former meaning, not in the latter. This means that a:s:b and a:s:c might disagree quite badly when the two different meanings are compared – as they do in the case of 0.0:0.1:0.31 and 0.0:0.1:1.0. These two meanings should arguably have totally different input syntaxes, but for historical reasons, they don't. A very different question is whether those should always agree on prefixes when both are used in the exact endpoint sense, which does seem like a nice property if we can guarantee it. For the inexact endpoint meaning, that should happen naturally, since the values produced should only depend on the literal values of a and s and b and c should only affect the length of the ranges.

In the inexact endpoint case, when the stop point cannot be hit exactly, there's no justification for non-literal interpretation of the start and step. So while 0.0:0.1:1.0 should not include 0.30000000000000004, I think that 0.0:0.1:1.01 should include it. In the 0.0:0.1:1.0 case the stop point is hit exactly, so one needs to interpret at least one of 0.0, 0.1 and 1.0 non-literally. One interpretation is to take the start and step literally, and that works because 10*0.1 rounds to 1.0 but it is not exactly equal to 1, so that's only one possible interpretation. Another viable interpretation is to take 0.0 and 1.0 literally and infer that 0.1 is intended to mean 1/10 exactly. In the 0.0:0.1:1.01 case, no such justification can be made since no interpretation makes start + n*step round to stop. So we don't have enough conflicting constraints to justify any guessing – we only have a literal start and step value and a number of steps, implicitly given by the stop point, which we know that we must stop short of. The only justifiable interpretation is taking start and step literally.

Part of the trouble here is that there's no set of principles written down for what these float ranges should do – without that it's impossible to verify that those principles are consistent and can actually all be satisfied at the same time. That, of course, is significantly my fault for not writing them down sufficiently clearly when I first implemented the lifting behavior.

[timholy]

I worry that the behavior you want is much more difficult to formalize/document than what I was aiming for. But I want to give you space to try, so I'll focus my efforts solely on fixing linspace and leave the rest to you. Note this is a 1.0 issue since the values might change.

[StefanKarpinski]

I guess there are several interpretations we’ve “floated” (get it?) for the meaning of a:s:b, in terms of real numbers a, s and b and an integer n which are the “true start, step and stop values and number of steps” of the range intended.

The literal interpretation. Here a is taken to be the real value of a, s is the real value of s, and b is the real value of b, and n is the maximum integer such that a + n*s ≤ b. The values of the range are precisely [a + k*s for k = 0:n].

This interpretaiton has the benefit that it is very simple and easy to understand. Its most significant drawback is that 0:0.1:0.3 would not end with 0.3 but would end with 0.2 instead, and have length 2 instead of 3, which would be a very frequent irritation. It has several less severe problems where intermediate points are not what one niavely expects, such as 0.30000000000000004 in 0:0.1:1. This is the interpretation that @andreasnoack favors, if I understand correctly.

The lifted exact + literal inexact interpretation. If there exist real numbers a, s and b and an integer n such that a rounds to a, s rounds to s, b rounds to b and a + n*s == b then in this interpretation, a:s:b may be interpreted as “lifted” in the sense of computing the values of the range based on the real values a, s, b and n and then rounding exact intermediate values back to floating-point. If no such real numbers can be found, however, the range must be interpreted literally in the sense above. In particular, if no such reals exist, as in the case of 0:0.1:0.31, in this interpretation the range must be interpreted literally, which implies that 0:0.1:0.31 == [k*0.1 for k = 0:3] and in particular last(0:0.1:0.31) == 0.30000000000000004.

This is what I implemented with my FloatRange type. My implementation used rationals with "reasonably sized numerator and denominator" as a lifting space. It lifted the floats to rational values which rounded exactly the inputs a, s and b values and then used the straightforward a + k*s formula on integer floating-point values (which have exact arithmetic, even in floating-point), projecting back down at the end by dividing by the least common denominator. This implementation cannot represent LinSpace objects, which was annoying. There was also an outstanding bug where one could havefirst(a:s:b) != a because of floating-point values where d*a/d != a.

The always lifted interpretation? This is what we currently have. I’m not sure what principle this follows to be honest because I don’t fully understand it. Can you write up what the rules for this are? I know that it uses double precision numbers to allow linear interpolation to get behavior similar to my rational lifting technique. I thought it was the same as my approach but with a different lifting method, but that seems not to be the case since it no longer respects the rule that inexact ranges must be interpreted literally. E.g. last(0:0.1:0.31) == 0.3 whereas the inexact literal rule would dictate that the last point be 0.0 + 3*0.1 == 0.30000000000000004.

[StefanKarpinski]

My current best understanding of what you implemented, Tim, is that you still use rational lifting but then approximate the lifted rationals with double-precision values – which must still round to the original inputs (i.e the high part must equal the input). Intermediate values are then computed using linear interpolation on double-precision values, which allows the type to represent both float ranges and linear spaces accurately, and avoids doing an expensive division operation (although I don't recall if the end result ended up being faster or not).

My only objection to any of this is doing rational lifting in the case where the stop value cannot be hit exactly – because that's the case where deviating from the user's literal input value cannot be justified. I suspect that others (e.g. @andreasnoack) may object even in some cases where the stop value is hit exactly, because there do exist ranges where a literal interpretation hits the same stop value but gives different intermediate values. Of course, those are precisely the cases where people find the literal interpretation unintuitive and they file issues, but I can understand the objection.

[timholy]

On my phone (no wifi), but briefly I think we have an awkward hybrid now and something has to be changed. I would probably change so we decide to lift to rational based on a and s alone. Either way, compute the length in 2x precision and call range (thus discarding b from further participation).

[StefanKarpinski]

I suspect I would prefer to go the other way and do the lifting and interpolation the way you're doing it now but reinstate the exact interpretation rule in cases where the stop point can't be hit. We should figure out what the pros and cons of each direction are.

[GlenHertz]

Shouldn't a Julian user use something like d"0":d"0.1":d"1" to get numbers on a grid or rational numbers 0//1:1//10:1//1? I'm not sure if floating point can be made to avoid all surprises but I do appreciate the efforts.

[StefanKarpinski]

Decimal only solves the problem if the numbers are representable exactly in decimal, which most aren't. The rational numbers already (and have always) worked. Arguably, that's the correct way to do this, except that it's quite slow. It's a considerably better API, however, since you can at least express what you'd like to produce accurately.

[oscardssmith]

is there a way that we could spell fancy float ranges as float(rational range) or something similar and optimize it correctly?

[andreasnoack]

This is the interpretation that @andreasnoack favors, if I understand correctly.

Just to clarify my view here: yes, I do think that we should base floating point operations on the exact binary floating point values of the inputs and not on low-denominator rational numbers in small neighborhoods around the inputs. My view is not that length(0.0:0.1:0.3)==3 would be a good solution. I don't think we can get sensible behavior for float ranges as long as we use decimal literals as input for binary float values (i.e. probably forever). Hence my view remains that we should just have LinSpace.

[timholy]

The good news is that now at least one has control over what happens using the exported interface:

julia> StepRangeLen(0.0, 0.1, 4)    # literal
0.0:0.1:0.30000000000000004

julia> range(0.0, 0.1, 4)           # allows lifting 
0.0:0.1:0.3

julia> typeof(range(0.0, 0.1, 4))
StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}

[StefanKarpinski]

So it seems that the consensus is that the traditional way of inputting floating-point ranges is simply inherently flawed, and that a better way of writing them would be the only way to really fix the "float range problem". However, we don't know what that would look like and the place to experiment with it is in packages, not the base language. To that end, I think that specifying how we want floating-point ranges to behave and making sure our implementation matches that specification is the best we can do for the base language, short of taking float ranges out entirely (which I don't think is viable).

[GlenHertz]

Like `float(0:1//10:3//10)`? That could give speed and clarity.

…

On Aug 12, 2017 9:15 PM, "oscardssmith" ***@***.***> wrote: is there a way that we could spell fancy float ranges as float(rational range) or something similar and optimize it correctly? — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#20373 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABopE1WYL0-24UZ60e4qBYGlOC0KclRDks5sXk4YgaJpZM4L0Nps> .