Working around floating-point inaccuracies
Hello everyone! For a library I write at the moment, I made a class for approximate geometric equality comparisons (and dubbed it `ApproxEq`). It's not a very lawful class, but there is one invariant that should hold absolutely always - if two entities coincide but have different representations, they should compare "equal". All instances for my geometric entities seem to uphold this invariant (btw, QuickCheck is awesome). Except for the `Plane` type, which gives me a false negative once in about 50000 runs. I've chosen to represent planes by two vectors: an arbitrary point on a plane, and a normal. To test this invariant I translate the plane by a vector parallel to the plane, and then compare the translated plane with the original. The general pattern to these failures is that one of the coordinates of the origin point of the untranslated plane is huge compared to the other two - by 4-6 orders of magnitude, while the translation vector is not particularly big and has about the same magnitude as the two smaller coordinates. I suspect that this difference throws off the function that tests if a point belongs to a plane (which is in turn used in the equality comparison) by filling lower digits of participating numbers with floating-point garbage. I can live with this false negative, but I'd prefer not to. There probably is a clever way to work around it, I just can't find it. Any ideas? Thanks in advance. The relevant code is below. ``` -- | A plane is represented by an arbitrary point on it and a normal vector. data Plane a = Plane (V3 a) (V3 a) deriving (Show) instance (Epsilon a, Floating a) => ApproxEq (Plane a) where p1@(Plane o1 n1) ~== p2@(Plane o2 n2) = res where n1' = normalize n1 n2' = normalize n2 parallel = n1' ~== n2' antiparallel = (-n1') ~== n2' originsAreTheSame = o1 ~== o2 originsBelong = pointOnPlane p1 o2 || pointOnPlane p2 o1 res = (parallel || antiparallel) && (originsAreTheSame || originsBelong) -- | Return True if a point lies on a plane. pointOnPlane :: (Epsilon a, Floating a) => Plane a -> V3 a -> Bool pointOnPlane (Plane o n) v = nearZero $ delta `dot` n' where delta = normalize $ o - v n' = normalize n ``` -- Michail.
Am 05.08.2018 um 14:30 schrieb Michail Pevnev:
The general pattern to these failures is that one of the coordinates of the origin point of the untranslated plane is huge compared to the other two - by 4-6 orders of magnitude, while the translation vector is not particularly big and has about the same magnitude as the two smaller coordinates.
I haven't checked whether this is really the case, but such failures tend to happen when the math is "badly conditioned" - i.e. if you're doing the equivalent of a/(b-c) somewhere, where b and c are so close to each other that most or all significant digits cancel out. Geometrically, a typical task would be trying to find the intersection of almost-parallel linear objects. Either your test code has the same instabilities as the tested code. In this case, your tested code isn't suitable for certain classes of situations; you'll want to find out which these are, and document them. Or try and find a different algorithm. In this case: Congratulations, testing revealed a real problem in your code :-) Or your test code has different instabilities. In that case, rewrite your test code. A different way to phrase this is: Your test code makes different assumptions (about numeric stability, that is) than the code under test. The *clean* way to deal with this is to calculate the error bounds, document them, and write test code that checks whether the code under test is behaving within the calculated error bounds. However, this is a pretty hairy thing to do; expect to spend several weeks to get the math worked out - or a few days with a numerics expert who already has seen it all. So maybe you'll want to think about alternatives; which of these are relevant depends a lot on what you're using your code for, and what kinds of qualities you want to achieve. Yet another approach might be using Rational instead of float. Rational is a net win only if you don't have iterations in your code, otherwise your numerators and denominators will grow on each iteration, eating up your memory. Regards, Jo
I've used an alternative algorithm as you suggested, and it works like a charm. Thanks for helping me out of my tunnel vision. -- Michail.
Michail Pevnev
Hello everyone!
For a library I write at the moment, I made a class for approximate geometric equality comparisons (and dubbed it `ApproxEq`). It's not a very lawful class, but there is one invariant that should hold absolutely always - if two entities coincide but have different representations, they should compare "equal".
The approximate-equality package (https://hackage.haskell.org/package/approximate-equality) provides a nicer way of doing that IMO. I've used that a few times when I needed an aprpimxate Floating type. Generally I think it is just better to use exact number types, like Rational, when implementing geometric algorithms though. That is what I tend to use in hgeometry as well ( that quite often a Fractional constraint is sufficient anyway). -- - Frank
participants (3)
-
Frank Staals -
Joachim Durchholz -
Michail Pevnev