Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Use error bound filter in isInCircle predicate instead of long double. #1162

Open
wants to merge 1 commit into
base: main
Choose a base branch
from

Conversation

tinko92
Copy link
Contributor

@tinko92 tinko92 commented Sep 14, 2024

This PR is edited to separate out the change to the orientation filter, now in #1184 .

For the incircle predicate, this PR proposes to not generally use long double (which increases precision a little but is not robust as the code notes) but to use a similar error bound filter instead (bound computation based on https://doi.org/10.1007/s10543-023-00975-x , I don't know the source for alternative expression of the incircle determinant, I have first seen it in CGAL). The rational for that is higher performance (roughly 4-10% speedup on my machine for perf_voronoi) and that a more robust computation here is probably motivated by preventing infinite loops in Delaunay Triangulation due to false cases of INTERIOR, which should be impossible with this version (except possibly for inputs near the lower end of the representable magnitude of double that produce denormalized results/underflow).

PS: I also have an implementation of exact predicates that I could port to geos to be run after failures of the fast filters if it would be considered beneficial to the library by its maintainers. I don't know enough about its use cases to make a guess about that. Exact predicates guarantee consistency wrt to predicate calls. They would not guarantee consistency between predicates and computed geometries (e.g. predicate might say some polygons overlap but the computed intersection may then be empty after rounding to double coordinates in the output geometry, even if all intermediate computations were exact).

@pramsey
Copy link
Member

pramsey commented Sep 19, 2024

I don't have a voronoi test yet, but my performance suite shows no differences, with one exception: the prepared point-in-poly test seems to be consistently about 10% faster.

@tinko92
Copy link
Contributor Author

tinko92 commented Sep 20, 2024

@pramsey Thanks for the feedback. Interesting, point-in-poly may benefit from a faster orientation but I would have expected that to be negligible if its implemented with some ray counting then I'd expect runtime to be dominated by walking over segments that don't intersect the ray and need no orientation check (at least for polygons with many segments). I will try to look into that when I find time.

I would expect the change in the orientation test to be measurable in e.g. convex hull or triangulation without Delaunay requirement and synthetic benchmarks (I ran perf_orientation from the benchmark directory and compared the number of Iterations since the ns timing seems too coarse for such a micro benchmark, I think running it in a larger algorithm is a more meaningful benchmark for this) and the change in the incircle test to be measurable in benchmarks like Delaunay Triangulation or Voronoi diagram (I ran perf_voronoi from the benchmark directory in the geos repository with release config for the numbers in the opening post). The incircle test maybe very significantly on non-x86-architectures where no machine hardware accelerated 80-bit floats are available and long double would compile to function calls to some soft float implementation, e.g. https://godbolt.org/z/EosY5ehjq .

@pramsey
Copy link
Member

pramsey commented Sep 26, 2024

I hooked a Delaunay build into my test harness but couldn't measure any different between this PR and 3.13.

@tinko92
Copy link
Contributor Author

tinko92 commented Sep 29, 2024

I hooked a Delaunay build into my test harness but couldn't measure any different between this PR and 3.13.

I tested as follows using the Delaunay/Voronoi benchmarks in the benchmark folder from the main repo:
I checked out the main branch, went into the directory and:

cd geos
mkdir build_gcc && build_gcc
cmake -DCMAKE_BUILD_TYPE=Release -DBUILD_BENCHMARKS=ON -DCMAKE_CXX_FLAGS="-march=native" .. 
make #will complain about missing override in another benchmark but doesn't affect the build of perf_voronoi
echo 1 | sudo tee /sys/devices/system/cpu/intel_pstate/no_turbo #only on intel
cpupower frequency-set -g performance #only on intel 
for i in {1..10}; do taskset -c 2 ./bin/perf_voronoi >> baseline; done;
grep "BM_DelaunayFromSeq/1000000" baseline
BM_DelaunayFromSeq/1000000        5050644908 ns   5045715175 ns            1
BM_DelaunayFromSeq/1000000        5049845175 ns   5045245900 ns            1
BM_DelaunayFromSeq/1000000        6646025308 ns   6638513791 ns            1
BM_DelaunayFromSeq/1000000        6649783153 ns   6642201348 ns            1
BM_DelaunayFromSeq/1000000        6648571528 ns   6641021393 ns            1
BM_DelaunayFromSeq/1000000        5053053815 ns   5048554170 ns            1
BM_DelaunayFromSeq/1000000        5049133353 ns   5044510485 ns            1
BM_DelaunayFromSeq/1000000        5052486875 ns   5047828776 ns            1
BM_DelaunayFromSeq/1000000        5050604032 ns   5046051328 ns            1
BM_DelaunayFromSeq/1000000        5051929265 ns   5047293562 ns            1

#checkout the branch with this patch, same steps, then
for i in {1..10}; do taskset -c 2 ./bin/perf_voronoi >> patched; done;
grep "BM_DelaunayFromSeq/1000000" patched
BM_DelaunayFromSeq/1000000        4659588868 ns   4654149493 ns            1
BM_DelaunayFromSeq/1000000        4654536370 ns   4649391288 ns            1
BM_DelaunayFromSeq/1000000        4653135324 ns   4648135157 ns            1
BM_DelaunayFromSeq/1000000        6127312010 ns   6119523080 ns            1
BM_DelaunayFromSeq/1000000        6125711504 ns   6118105844 ns            1
BM_DelaunayFromSeq/1000000        4655963070 ns   4651342713 ns            1
BM_DelaunayFromSeq/1000000        4652838792 ns   4648184104 ns            1
BM_DelaunayFromSeq/1000000        4654345257 ns   4649771922 ns            1
BM_DelaunayFromSeq/1000000        4666023790 ns   4661365412 ns            1
BM_DelaunayFromSeq/1000000        4653681302 ns   4649277328 ns            1

I honestly wish I knew, where the outliers come from with turbo disabled, maybe the P-core/E-core (it's a mobile i7-12800H) thing can affect it despite fixing the CPU core with taskset. Same steps with CXX=/usr/bin/clang++ before the cmake command yield for me:

grep "BM_DelaunayFromSeq/1000000" baseline
BM_DelaunayFromSeq/1000000        4329365684 ns   4323526473 ns            1
BM_DelaunayFromSeq/1000000        5433784294 ns   5427474646 ns            1
BM_DelaunayFromSeq/1000000        4148863139 ns   4144892923 ns            1
BM_DelaunayFromSeq/1000000        4139543832 ns   4135714625 ns            1
BM_DelaunayFromSeq/1000000        5428346793 ns   5422042258 ns            1
BM_DelaunayFromSeq/1000000        5439176317 ns   5432822740 ns            1
BM_DelaunayFromSeq/1000000        5432680958 ns   5426345764 ns            1
BM_DelaunayFromSeq/1000000        4140026534 ns   4136116714 ns            1
BM_DelaunayFromSeq/1000000        4149153847 ns   4145284595 ns            1
BM_DelaunayFromSeq/1000000        5431052965 ns   5424639386 ns            1

grep "BM_DelaunayFromSeq/1000000" patched 
BM_DelaunayFromSeq/1000000        3767896479 ns   3764021258 ns            1
BM_DelaunayFromSeq/1000000        4923848243 ns   4917679386 ns            1
BM_DelaunayFromSeq/1000000        3764235827 ns   3760602930 ns            1
BM_DelaunayFromSeq/1000000        3762541532 ns   3758951324 ns            1
BM_DelaunayFromSeq/1000000        3765034916 ns   3761276158 ns            1
BM_DelaunayFromSeq/1000000        3764714564 ns   3761178183 ns            1
BM_DelaunayFromSeq/1000000        4928664355 ns   4922710265 ns            1
BM_DelaunayFromSeq/1000000        4922105205 ns   4916260250 ns            1
BM_DelaunayFromSeq/1000000        3767456227 ns   3763747766 ns            1
BM_DelaunayFromSeq/1000000        3760288885 ns   3756814473 ns            1

Similar results for the other workloads.

And on an EPYC-Milan I don't have these outlier timings, so I post just the first sample for each workload:

grep 100000 baseline | head -n5
BM_DelaunayFromSeq/1000000        4916505592 ns   4916348118 ns            1
BM_DelaunayFromGeom/1000000       4671339867 ns   4671204309 ns            1
BM_VoronoiFromSeq/1000000         5506559255 ns   5506240288 ns            1
BM_VoronoiFromGeom/1000000        5435366604 ns   5435120888 ns            1
BM_OrderedVoronoiFromGeom/1000000 7680563331 ns   7680220398 ns            1

grep 100000 patched | head -n5 
BM_DelaunayFromSeq/1000000        4427675159 ns   4427496208 ns            1
BM_DelaunayFromGeom/1000000       4462223273 ns   4461665059 ns            1
BM_VoronoiFromSeq/1000000         5045797628 ns   5045404499 ns            1
BM_VoronoiFromGeom/1000000        4988454644 ns   4987171839 ns            1
BM_OrderedVoronoiFromGeom/1000000 7392607964 ns   7392126698 ns            1

@pramsey
Copy link
Member

pramsey commented Oct 25, 2024

Thanks for the nice writeup of your perf_* results. I did the same tests on my MacBook Pro and I'm afraid I see no consistent advantage for this particular optimization.

./bin/perf_voronoi --benchmark_filter=1000000 --benchmark_repetitions=10 > baseline
grep mean baseline

BM_DelaunayFromSeq/1000000_mean          1885685129 ns   1885561400 ns           10
BM_DelaunayFromGeom/1000000_mean         1859082321 ns   1858690300 ns           10
BM_VoronoiFromSeq/1000000_mean           2238724242 ns   2238228000 ns           10
BM_VoronoiFromGeom/1000000_mean          2254024629 ns   2253934700 ns           10
BM_OrderedVoronoiFromGeom/1000000_mean   2975758075 ns   2975618100 ns           10

./bin/perf_voronoi --benchmark_filter=1000000 --benchmark_repetitions=10 > patched
grep mean patched

BM_DelaunayFromSeq/1000000_mean          1906208617 ns   1906192700 ns           10
BM_DelaunayFromGeom/1000000_mean         1914398154 ns   1913137200 ns           10
BM_VoronoiFromSeq/1000000_mean           2288991796 ns   2288465000 ns           10
BM_VoronoiFromGeom/1000000_mean          2295409392 ns   2294846200 ns           10
BM_OrderedVoronoiFromGeom/1000000_mean   2993584858 ns   2993345900 ns           10

@pramsey
Copy link
Member

pramsey commented Oct 25, 2024

Would it make sense to gate this behind something like

#if defined(__x86_64__) || defined(__x86_64) || defined(__amd64__) || defined(__amd64) || defined(_M_X64) || defined(_M_AMD64)
#   define ARCHITECTURE_X86_64

@tinko92
Copy link
Contributor Author

tinko92 commented Oct 26, 2024

@pramsey Thanks for getting back to me on this.

The test results are interesting. I do not have a mac to test. But using clangs target option, I find that it seems to default to compile long double to 8 byte, so same as double:

echo "int s = sizeof(long double);" | clang++ -x c++ -S - -O3 -target aarch64-apple-darwin -o - | grep '.long'
	.long	8                               ; 0x8

So my patch would not benefit the circle predicate performance there from switching down to double. Instead it will cost a little performance due to the slightly more expensive error bound computation (visible in your benchmark, seems to be ~1-2%).

On the other hand, for architectures that are already compiling the current code to 64-bit precision/double, this also means that the current code will be a lot less robust. If I generate 100k degenerate cases, the current code will produce an opposing result (interior when the truth is exterior or the other way around) in 8 cases with 80-bit precision but in ~14k cases with 64-bit precision. All code for the test is in this gist: https://gist.github.com/tinko92/86126244ed04b65aeb09c9223ec4c37c#file-test-sh (note that this test code will exclusively produce degenerate cases where the test point is near the circle boundary up to roughly machine precision).

The patched code will never confuse interior and exterior, regardless of availability of 80-bit precision doubles but always produce the correct sign or BOUNDARY. If you run that gist code on your Macbook, I would expect that you see many results in which the 2nd and 3rd column show opposite signs.

So overall my opinion would be to not gate it. The architectures that lack 80-bit long doubles suffer a 2% performance hit but will benefit from it a lot robustness-wise (with an understanding of robustness where "boundary" is a safe result while confusing exterior and interior is bad) and the architectures that have 80-bit long doubles should benefit from it somewhat performance-wise and still a little robustness-wise (guarantee of not confusing exterior and interior).

This understanding of robustness is based on the idea that confusing interior and exterior could lead to infinite loops in triangle flipping for Delaunay Triangulation and should be avoided for this reason. But I don't know the specifics of geos' Delaunay Triangulation algorithm or the general motivation behind having an incircle-predicate with increased robustness in geos in the first place.

Edit: Discussion of orientation, now in #1184, removed.

@dr-jts
Copy link
Contributor

dr-jts commented Oct 27, 2024

It would be better if this was split into 2 PRs, one for each algorithm. This would separate the different discussions, tests, and changes.

@tinko92
Copy link
Contributor Author

tinko92 commented Oct 27, 2024

That makes sense, since most of the discussion is around the incircle predicate, I'll split out the orientation predicate into a separate PR.

Edit: If there is interest in making this method robust as the name suggests, I can add a commit with a slightly transformed version of Shewchuk's public domain implementation, where all defines are changed to constexpr functions. https://gist.github.com/tinko92/8229b4d90ed0dddc410ec22ce5d75b2f . This would guarantee correct predicate results for every input that does not cause over or underflow (so anything with coordinates whose absolute values are roughly in the [1e-75,1e75] plus {-0., 0.} range. It would only be executed in case of filter failure and therefore not contribute to runtime for non-degenerate inputs.

Edit2: After reading more from the docs, I understand that I should probably prose that to JTS first, this may take some time, though to port to Java.

@tinko92 tinko92 changed the title Use tighter bounds and faster filters in orientation and isInCircle p… Use error bound filter in isInCircle predicate instead of long double. Oct 27, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants