Floating point comparisons: inform or protect the user?

[trexfeathers]

I recently promised to look at changing Constraints to use approximate comparisons, to avoid floating point comparison problems (stemming from #4700). I now think there are compelling arguments for not doing this, but to instead lean on / enhance the existing warning in the user guide.

What does everyone else think? There was a lot of support for changing to approximate equality in a recent team discussion.

---

1. Principle of least surprise - Iris should be consistent with the rest of the scientific Python community.
   • This is a well known programming concept.
   • AFAICT no other scientific Python hides this problem from end users, rather offering specific solutions and/or guidance (np.isclose(), pd.testing.assert_frame_equal(check_exact=False)).
   • A particular problem since we know many of our users start their Python journey with Iris.
   • This point is similar to our import debate - Populate namespace: lazy_loader #4710.
2. The behaviour change won't be limited to Constraints - much of Constraint equality uses Cell equality.
   • Can we be sure this is desirable?
   • If it is desirable, do we also need to be consistent by implementing approximate equality throughout Iris? This becomes a much bigger question!
3. Further opportunity for confusion/controversy: should we also make similar changes to > / >= / <= / < operations on Cells?
   • These are used when users pass functions to Constraints, consistent behaviour might be expected.
   • IMO this would be consistent if we were already modifying Cell.__eq__() - see point 2 - but conversation with @bjlittle suggests that some of us might see this as a line to not cross.
4. Code structure increases change of future divergence.
   • Constraints approach equality in several different ways depending on the specific circumstance, so == and in are used on many different lines, some calling down to Cell.__eq__() while others are doing raw comparisons of values.
   • Very easy for future misunderstanding/skim-reading the code and not being consistent with use of exact or approximate equality.
5. Performance reduction
   • The current code makes it clear that performance is a significant concern for Constraints.
   • Approximate equality involves more operations, so will likely be somewhat slower than exact equality.

Further reading:

• BUG: Tolerant crs comparison #2062
• _CoordConstraint.extract()
• Cell.__eq__()

[rcomer]

I know we’ve said we don’t want to proliferate warnings, and in general I am in agreement with that. But if a cell doesn’t have bounds and other is a float, perhaps it is reasonable for cell equality to throw a warning that this comparison is ill-defined?

[trexfeathers]

Ugh I didn't even think about the significance of bounds here.

But in some cases isn't one of the bounds identical to the point? What if both those numbers are infinitesimally higher/lower than the Constraint value?

[rcomer]

😬

[trexfeathers]

Team have discussed an alternative API something like ApproxConstraint, which explicitly handles the floating point cases.

• Slower, easier version: convert the Constraint into a range function Constraint
• Faster, trickier version: work out what the actual floating point representation of the provided number should be, accounting for the type of the coordinate in question, then use that value in the Constraint instead.

[rcomer]

Constraints on DimCoords already find the nearest matching cell if you are matching to anything other than a Cell

iris/lib/iris/_constraints.py

Lines 308 to 321 in caae450

	try_quick = isinstance(coord, DimCoord) and not isinstance(
	self._coord_thing, Cell
	)
	# Simple, yet dramatic, optimisation for the monotonic case.
	if try_quick:
	try:
	i = coord.nearest_neighbour_index(self._coord_thing)
	except TypeError:
	try_quick = False
	if try_quick:
	r = np.zeros(coord.shape, dtype=np.bool_)
	if coord.cell(i) == self._coord_thing:
	r[i] = True

I've not thought through the details, but it seems like we should be able to use this.

[trexfeathers]

When I read this code I assumed this meant that nearest_neighbour_index() needs monotonicity to work / work efficiently.

Alternatively, not using this for AuxCoords could just be an oversight.

[rcomer]

Yep, there's nothing in the docstring for nearest_neighbour_index to suggest it has to be monotonic. It has to be 1D, but extract only works on 1D coordinates anyway...