Normalise angle fix

[jhaiduce]

This modifies `eckit::geometry::normalise_angle`` so that it uses a modulo operator to find the correct output value, rather than repeatedly adding or subtracting 360 from the input. This makes the function orders of magnitude faster for very large input values.

Also added a few additional checks to test_coordinate_helpers to verify that normalise_angle works correctly on negative inputs and very large (positive or negative) inputs.

Fixes #102.

[FussyDuck]

All committers have signed the CLA.

[fmahebert]

Thanks for this iteration on the code (that IIRC I put in)... the diff itself looks fine to me but I have a few questions about the intention of this code change:

• did you check whether there's any performance implication for "small" inputs? This method would typically be called to remove "a few" 2pi of phase, not 10^many 2pi of phase, so we should make sure we didn't move the performance problem to the typical case. The lines of code as written don't look slow, so I'm not overly concerned, but would like to know if this was tested.
• does it really make sense to mod away 10^many 2pi of phase? Depending on how large the input is, there's a significant risk of modding away some/many/all of the significant digits of the input, so that the output of normalise_angle is just garbage bits. I think from an algorithmic point of view, it would make more sense to abort/throw if the angle needs more than a few unwindings. This way you can be sure the function is actually preserving ~ double precision.

The above are my concerns, but aren't a change request, as this is not my repo :)

[jhaiduce]

did you check whether there's any performance implication for "small" inputs? This method would typically be called to remove "a few" 2pi of phase, not 10^many 2pi of phase, so we should make sure we didn't move the performance problem to the typical case. The lines of code as written don't look slow, so I'm not overly concerned, but would like to know if this was tested.

I didn't actually test that. Might be worth trying though; it would be simple enough to write something that calls normalise_angle in a loop and time it.

does it really make sense to mod away 10^many 2pi of phase? Depending on how large the input is, there's a significant risk of modding away some/many/all of the significant digits of the input, so that the output of normalise_angle is just garbage bits. I think from an algorithmic point of view, it would make more sense to abort/throw if the angle needs more than a few unwindings. This way you can be sure the function is actually preserving ~ double precision.

The magnitude of the input angle is a limiting factor for having useful precision here, but I expect this kind of precision loss probably would have happened with the previous implementation too, as there would be a little precision loss with each addition/subtraction of 360. In the example I gave (input of 1e36), there is too much precision loss for the output to be usefully precise. However, if one gives an input of say 1e14 there is enough precision remaining to be potentially useful and the time to run the loop is already very large at that point.

[jhaiduce]

did you check whether there's any performance implication for "small" inputs?

I just wrote a small test that calls normalise_angle 10,000,000 times with random inputs between -720 and 720. Both implementations average about 1.68-1.73 ms per call. I'd have to do some more careful testing to know for sure if there is a performance change for "small" inputs, but it looks like the impact is on the order of 3% or less.

[jhaiduce]

I think from an algorithmic point of view, it would make more sense to abort/throw if the angle needs more than a few unwindings.

I considered adding something like this, but it would have required making an arbitrary decision on how large an angle should be allowed before the exception is triggered, and that arbitrary threshold would then be imposed on client code. Not checking this effectively leaves the decision up to the calling code how much precision loss is acceptable, which IMO is a more appropriate place since the calling code has visibility into how the output will be used.

[fmahebert]

I just wrote a small test that calls normalise_angle 10,000,000 times with random inputs between -720 and 720. Both implementations average about 1.68-1.73 ms per call. I'd have to do some more careful testing to know for sure if there is a performance change for "small" inputs, but it looks like the impact is on the order of 3% or less.

Thank you for confirming this!

I expect this kind of precision loss probably would have happened with the previous implementation too, as there would be a little precision loss with each addition/subtraction of 360

I considered adding something like this, but it would have required making an arbitrary decision on how large an angle should be allowed before the exception is triggered, and that arbitrary threshold would then be imposed on client code. Not checking this effectively leaves the decision up to the calling code how much precision loss is acceptable, which IMO is a more appropriate place since the calling code has visibility into how the output will be used.

I agree with everything you've written. Producing a low-precision output (as in this PR) for "large" inputs is definitely an improvement over hanging (as in develop), so this PR sounds like a good improvement to me. We'll see what the ECMWF team prefers regarding precision-loss vs erroring...

[wdeconinck]

Not saying you should, but a hybrid solution could also be created?

if abs(lon) < some_treshold then "old method", else "use modulo"

[jhaiduce]

@wdeconinck, implementing a switch between the two methods as you describe would be straightforward to do, but I'm not sure there's a benefit to keeping both since the modulo method produces identical output in comparable or faster time for every case I've tested.

[wdeconinck]

@jhaiduce if the output is bit identical, then I am OK with this change.
I saw you were mentioning a precision loss though. Could you clarify?

[jhaiduce]

@wdeconinck in the case of very large numbers being input to the function the output will have fewer significant figures than the input value. This is true for both methods. If the input is around 3,600 degrees for example, the output will have one less significant digit than the input for both methods but the new method will return faster. If the input is 36,000 degrees the output will have two fewer significant digit than the input did, and this will be true for both methods but the new method will give a result in about 100x less time.

This is limited by the floating-point precision of the double data type. Since a double has about 16 significant digits, you can't expect useful output for inputs larger than about 1e18 for this function because 0-360 range would then be smaller than the least significant digit of the input (the function can't provide a precise output because the precision of the input didn't include the 0-360 range in the first place). The new method doesn't improve on this precision loss, but it produces an output with the same precision loss in much less time.

[wdeconinck]

Thank you for clarifying @jhaiduce . I have also no objection. @pmaciel ?

[pmaciel]

Hi @jhaiduce, @wdeconinck, maybe there is, in the future, room to pursue a solution removing integer multiples of 360 so as not to have any loss of precision whatsoever (as argued) but the solution found is quite elegant already.

I'm of the opinion the performance motivation isn't the best though, as values wildly different from the typical ones are probably symptomatic of other (likely bigger) issues -- but it doesn't detract from that this is an improvement.

The only final thing I would add is to run clang-format, my eyes are now quite atuned to that :-)

[jhaiduce]

maybe there is, in the future, room to pursue a solution removing integer multiples of 360 so as not to have any loss of precision whatsoever.

I suspect in the extreme cases (where the last digit significant digit is already larger than 360), the only way to improve the precision is to change the type of the input arguments.

I'm of the opinion the performance motivation isn't the best though, as values wildly different from the typical ones are probably symptomatic of other (likely bigger) issues -- but it doesn't detract from that this is an improvement.

Agreed, if one passes in values much larger than 360 it probably reflects the function being used for something other than its intended use, perhaps passing in other coordinate types where an angle was expected (incidentally the way I noticed that an improvement could be made here was by accidentally passing in uninitialized data).

The only final thing I would add is to run clang-format, my eyes are now quite atuned to that :-)

Just ran clang-format and committed the change.

[wdeconinck]

Thanks @jhaiduce for your contribution! Looks perfect now.

[wdeconinck]

@jhaiduce it appears the results are not bit-identical after all. Follow-up PR #111 should address this, with a threshold where the previous formula stays applied. It's not nice but I am not sure if other solutions could be found.

[jhaiduce]

My statement that they were identical was based on passing the test built in to eckit, but those only covered a few cases. Apparently the downstream tests are more thorough.