Conditional indexing issue with float precision

General Usage
1

Hi!
I came accross a float-precision issue that I want to share and ask if there is a solution for my problem.

I define an array of coordinates, an equally-sized array for the result, and I want to fill the second array with numbers in segmentwise fashion through conditional indexing. This works fine for integer coordinates, but gets totally confusing (for me) with float coordinates.
I realize, isapprox points into the right direction, but its not fixing the issue entirely.

Here’s my code

#integer version
x = 0:1:20
seg_size = 2
seg_period = 4

I = zeros(Int64,size(x))
I[((x .% seg_period).>=0) .& ((x .% seg_period).<seg_size)] .= 2
I[((x .% seg_period).>=seg_size) .& ((x .% seg_period).<seg_period)] .= 1
println(I)

#float version
x = 0:.025:.5
seg_size = .05
seg_period = .1

I = zeros(Int64,size(x))
I[((x .% seg_period).>=0) .& ((x .% seg_period).<seg_size)] .= 2
I[((x .% seg_period).>=seg_size) .& ((x .% seg_period).<seg_period)] .= 1
println(I)

#float with isapprox
I = zeros(Int64,size(x))
I[(((x .% seg_period).>0) .| ((x .% seg_period).≈0)) .& ((x .% seg_period).<seg_size)] .= 2
I[(((x .% seg_period).>seg_size) .| ((x .% seg_period).≈seg_size)) .& ((x .% seg_period).<seg_period)] .= 1
println(I)

output:
[2, 2, 1, 1, 2, 2, 1, 1, 2, 2, 1, 1, 2, 2, 1, 1, 2, 2, 1, 1, 2]
[2, 2, 1, 1, 2, 2, 2, 1, 2, 2, 2, 1, 1, 2, 2, 1, 2, 2, 2, 1, 1]
[2, 2, 1, 1, 2, 2, 1, 1, 2, 2, 1, 1, 1, 2, 1, 1, 2, 2, 1, 1, 1]

The integer version (first output line) is nicely alternating between 1 and 2, which is the expected result.
The float version (second output line) has several mis-assignments due to float precision.
The float version using isapprox (third output line) can solve some issues, where the >= condition failed due to float precision, but cannot solve the issue where the < condition becomes true although the arguments are approximately equal.

So my question:

  • Is there any better way to do things like that while keeping the code quite readable?
  • Is there or should there be a approximate-version of >, <, >=, <= so that such an example (with appearantly very simple numbers) would work as naively expected?

Looking forward to your suggestions.

2

This is to be expected because of floating point error:

If there is an outcome for a calculation that follows from theory, you should generate that directly.

1 Like
3

Hi Tamas_Papp,
thank you for your quick reply.

In fact, I’m fully aware that this is a limitation of the floating point representation and not a Julia bug!

However, the question is, how do I better address segments on a floating-point coordinate system without gaps and without deviations in segment width?

I cannot fully understand your last sentence suggestion.

4

I don’t understand what you are trying do, since

(x .% seg_period) .< seg_period

should always be true in theory by construction, and when it isn’t, that’s just floating point error.

In this particular case, having non-exactly representable numbers on the boundary is also problematic, so could use rational numbers.

classify(x, seg_period, seg_size) = (x % seg_period) < seg_size ? 2 : 1
classify.(0:1:20, 4, 2)
classify.(0:(1//40):(1//2), 1//10, 1//20)
5

Thanks.

You are right that the quoted condition doesnt make sense… The example is of coarse simplified and some sense got lost in this process. However, if you remove the quoted condition, the outcome of all 3 versions is still different.

You’re example pretty much answers my question. I’ll give it a try.
Thank you!

6

Hi,

I figured out another workaround: just defining comparison operators taking into account approximate equality:

 (≲)(x,y) = <(x,y) | ≈(x,y)
 (⋦)(x,y) = <(x,y) & !≈(x,y)
 (≳)(x,y) = >(x,y) | ≈(x,y)
 (⋧)(x,y) = >(x,y) & !≈(x,y)

So, replacing >= by ≳ and < by ⋦ seems to make the range examples work with floats work as expected for rationals/integers.

7

This may fix operations which should “seemingly work” (eg 0.3 \cdot 3 = 9, while with Float64 you have 0.3 * 3 !== 9), but depending on your input and what is “expected”, you could find examples where it does not work.

If you want exact arithmetic for rationals, you should just use rationals.

8

In the example you give, values follow the following properties:

  • the range step (0.025) is not exactly represent able as a float
  • neither is the value assigned to seg_size but it is a multiple of the range step
  • the value seg_period is both representable and a multiple of the range step

Are those properties true in your real use case, or is it only true of your simplified MWE here?
If in particular it is always true that seg_size and seg_period are multiples of the range size, then you could (and should) only work with integers.

If, on the other hand, these properties do not hold in your real application (eg x values are not equally spaced), then your MWE might not be representative enough.

A key question here is: how often does it happen that one element of x falls right on the boundary between two segments? And can these situations be known in advance?

2 Likes
9

I agree that there are likely examples that still dont work “as expected”.

With regard to my problem: I try to model a physical experiment.

  • So usually range steps are dictated by the experiment but could of course (inconveniently) be mapped to integers.
  • seg_size is not necessarily multiple of the range step.
  • seg_period is also given by the real-world experiment and cannot be freely chosen.

Trying to answer your key question: as long as I dont think about the numerics, to me it makes intuitive sense to have any kind of segments start and end at an coordinate. In a floating-point numerics context, this is obviously a very bad idea.
Maybe I have to rethink the way I approach these kinds of problems to have segment borders usually right between two coordinates.

Thanks for all your input.

10

IIUC even in your real use case x values are equally spaced?
In that case, my intuitive sense would be to place segment boundaries right in the middle between two neighboring coordinates. But I guess it all depends on what you’re trying to achieve and you probably know better.

I would tend to think that if all your numbers come from physical measurements (or a model thereof), it should almost never happen that two numbers are equal (to machine precision) because measurement instruments (or models) tend to be affected by far higher uncertainties than 1e-15.

This is probably a very naive question, but have you checked that you actually have significant issues with floats in your real use case?

11

Of course, measured values are never equal. However, especially spatial coordinates or time-information is often not measured but defined/assumed and other quantities are than measured vs. these nominal values.
E.g., if you measure a wavefrom with an oscilloscope, the time information is float and equally spaced. It’s determined by the instruments operating frequency and settings. Time information is not measured but assumed to be correctly known. The measured voltage-values corresponding to these time points are definitively prone to measurement accuracy much worse than numerical accuracy, and comparing these for exact equality is of no sense.
However, modelling the experiment (with identical float time-base) or accessing certain segments of the measured data by indexing the waveform with conditions involving the time-array, is a quite usual use case, I would say.

12

Thanks for the explanation, that makes perfect sense.

Still, I the question remains: if neither the segment size nor the segment period are multiples of the range step, what are the odds of an element in x coinciding with a segment boundary?

13

seg_period usually is a multiple of step size…

14

OK. In that case, let me reformulate your initial MWE in this way:

# x = 0:step:x_max
# x[i] = (i-1) * step
# x_max = N * step
step = 0.02;
N = 20;

# seg_size is not a multiple of step
seg_size = .05;

# ... but seg_period is
P = 5;
seg_period = P*step;

Then you initial implementation looks like this:

I = zeros(Int64,N+1);
for i in 1:N+1
    xi = (i-1)*step

    if xi % seg_period < seg_size
        I[i] = 2
    else
        I[i] = 1
    end
end
println(I)
# -> [2, 2, 2, 1, 1, 2, 2, 2, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 1, 1, 2]
#                                     problem here ^

However, taking advantage of your knowledge that seg_period is a multiple of step, you know when problems will happen and avoid them:

I = zeros(Int64,N+1);
for i in 1:N+1
    xi = (i-1)*step

    if (i-1)%P == 0  ||  xi % seg_period < seg_size
        # if (i-1)%P==0, then (i-1)=k*P
        # =>  x[i] = (i-1)*step = k*P*step = k*seg_period
        # =>  x[i] % seg_period = 0
        I[i] = 2
    else
        I[i] = 1
    end
end
println(I)
# -> [2, 2, 2, 1, 1, 2, 2, 2, 1, 1, 2, 2, 2, 1, 1, 2, 2, 2, 1, 1, 2]
#                              no problem any more ^

Would something like this be applicable to your real use case?

1 Like
15

This would surely solve the problem. However, for the sake of readability, I tried to avoid having float-coordinates and integer representations of the same (either loop indices or another coordinate array). So I’ll stick to using rationals or using approximate comparisons.