Floating point errors

[caspervdw]

In #485 a large amount of floating point errors were solved. However there are some that remained. I listed the ones that popped up in pygeos (with GEOS 3.10.1):

• binary predicates on empty geometries
• voronoi_polygons on a point
• hausdorff_distance (in our seemingly normal test)
• offset_curve on an empty linestring

All through the CAPI. See pygeos ​pygeos/pygeos#441 for the specific tests where I encountered the FP “invalid” flags.

[pramsey]

I'm not actually seeing fe errors when mimicking these tests...

template<>
template<>
void object::test<5> ()
{
    runTest(
        "LINESTRING (0 0, 100 0, 10 100, 10 100)",
        "LINESTRING (0 100, 0 10, 80 10)", 0.001, 47.89);

    // check for floating point overflow exceptions
    int raised = fetestexcept(FE_OVERFLOW);
    ensure_equals(raised & FE_OVERFLOW, 0);
}

[caspervdw]

Thanks for looking into this issue. Sorry I didn’t provide a reproducing code snippet, I am not so comfortable writing in C++.

Looking at your example, I think you might reproduce the issue if you check for FE_INVALID or FE_ALL_EXCEPT flags.

[pramsey]

Yes, it's FE_INVALID, but I cannot track down where it is being set.

[pramsey]

Even when I track it down, the placement helps not at all. Somehow this is an overflow?

* thread #1, queue = 'com.apple.main-thread', stop reason = EXC_ARITHMETIC (code=EXC_I386_SSEEXTERR, subcode=0x20)
    frame #0: 0x0000000103948014 libgeos.3.11.0.dylib`geos::geom::LineSegment::projectionFactor(this=0x00007ffeefbfeb70, p=0x00000001051b3960) const at LineSegment.cpp:79:56
   76  	    double dx = p1.x - p0.x;
   77  	    double dy = p1.y - p0.y;
   78  	    double len2 = dx * dx + dy * dy;
-> 79  	    double r = ((p.x - p0.x) * dx + (p.y - p0.y) * dy) / len2;
   80  	    return r;
   81  	}
   82  	
Target 0: (test_geos_unit) stopped.
(lldb) p len2
(double) $0 = 8100
(lldb) 

[pramsey]

This ate my afternoon, and while I learned some fun things about architectures and debuggers, I'm no closer to eliminating these soft signals. You're just going to have to suppress/ignore them on your end. Make sure you're not turning on signals at some point in your code!

[caspervdw]

Thanks for investigating this. It is certainly possible to hide these errors from the user in pygeos. I was mostly afraid something undefined was happening under the hood but now that you traced the origin it seems all right.

This FE_INVALID flag is mostly occuring with NaN values in my (little) experience. Could it be possible that one of the coordinates in (p) contains nans?

[pramsey]

Not any of the coordinates being accessed. The z ordinate would have a NaN (since the input is 2d) but it's not read.

[caspervdw]

Were you able to narrow down the floating point operation that resulted in the FE_INVALID?

[pramsey]

Yeah, but it also doesn't make any sense. It said the problem was this division. Again, the inputs were perfectly normal doubles, no nans or anything.

https://github.com/libgeos/geos/blob/main/src/algorithm/distance/DiscreteHausdorffDistance.cpp#L88

[caspervdw]

Thanks again for diving into this!

[pramsey]

If you can track down any specific locations or changes, I'm happy to apply, just not having any luck on this snipe hunt thus far.