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

Add Angle::sinCosSnap to avoid small errors, e.g. with buffer operations #978

Merged
merged 1 commit into from Nov 7, 2023

Conversation

mwtoews
Copy link
Contributor

@mwtoews mwtoews commented Oct 30, 2023

This is a rather pedantic PR to avoid small rounding errors from std::sin(M_PI) and std::cos(M_PI_2), which are small non-zero values.

Consider a unit buffer geosop -a "POINT(0 0)" buffer 1 here with GEOS 3.12.0:
POLYGON ((1 0, 0.9807852804032304 -0.1950903220161282, 0.9238795325112867 -0.3826834323650898, 0.8314696123025452 -0.5555702330196022, 0.7071067811865476 -0.7071067811865475, 0.5555702330196023 -0.8314696123025452, 0.3826834323650898 -0.9238795325112867, 0.1950903220161283 -0.9807852804032304, 0.0000000000000001 -1, -0.1950903220161282 -0.9807852804032304, -0.3826834323650897 -0.9238795325112867, -0.555570233019602 -0.8314696123025455, -0.7071067811865475 -0.7071067811865476, -0.8314696123025453 -0.5555702330196022, -0.9238795325112867 -0.3826834323650899, -0.9807852804032304 -0.1950903220161286, -1 -0.0000000000000001, -0.9807852804032304 0.1950903220161284, -0.9238795325112868 0.3826834323650897, -0.8314696123025455 0.555570233019602, -0.7071067811865477 0.7071067811865475, -0.5555702330196022 0.8314696123025452, -0.3826834323650903 0.9238795325112865, -0.1950903220161287 0.9807852804032303, -0.0000000000000002 1, 0.1950903220161283 0.9807852804032304, 0.38268343236509 0.9238795325112866, 0.5555702330196018 0.8314696123025455, 0.7071067811865474 0.7071067811865477, 0.8314696123025452 0.5555702330196022, 0.9238795325112865 0.3826834323650904, 0.9807852804032303 0.1950903220161287, 1 0))

and since #973 for main the small values are formatted using scientific notation:
POLYGON ((1 0, 0.9807852804032304 -0.1950903220161282, 0.9238795325112867 -0.3826834323650898, 0.8314696123025452 -0.5555702330196022, 0.7071067811865476 -0.7071067811865475, 0.5555702330196023 -0.8314696123025452, 0.3826834323650898 -0.9238795325112867, 0.1950903220161283 -0.9807852804032304, 6.123233995736766e-17 -1, -0.1950903220161282 -0.9807852804032304, -0.3826834323650897 -0.9238795325112867, -0.555570233019602 -0.8314696123025455, -0.7071067811865475 -0.7071067811865476, -0.8314696123025453 -0.5555702330196022, -0.9238795325112867 -0.3826834323650899, -0.9807852804032304 -0.1950903220161286, -1 -1.2246467991473532e-16, -0.9807852804032304 0.1950903220161284, -0.9238795325112868 0.3826834323650897, -0.8314696123025455 0.555570233019602, -0.7071067811865477 0.7071067811865475, -0.5555702330196022 0.8314696123025452, -0.3826834323650903 0.9238795325112865, -0.1950903220161287 0.9807852804032303, -1.8369701987210297e-16 1, 0.1950903220161283 0.9807852804032304, 0.38268343236509 0.9238795325112866, 0.5555702330196018 0.8314696123025455, 0.7071067811865474 0.7071067811865477, 0.8314696123025452 0.5555702330196022, 0.9238795325112865 0.3826834323650904, 0.9807852804032303 0.1950903220161287, 1 0))

this PR would use exact values if the angles are a multiple of 45 degrees (pi/4) on the unit circle:
POLYGON ((1 0, 0.9807852804032304 -0.1950903220161282, 0.9238795325112867 -0.3826834323650898, 0.8314696123025452 -0.5555702330196022, 0.7071067811865476 -0.7071067811865476, 0.5555702330196023 -0.8314696123025452, 0.3826834323650898 -0.9238795325112867, 0.1950903220161283 -0.9807852804032304, 0 -1, -0.1950903220161282 -0.9807852804032304, -0.3826834323650897 -0.9238795325112867, -0.555570233019602 -0.8314696123025455, -0.7071067811865476 -0.7071067811865476, -0.8314696123025453 -0.5555702330196022, -0.9238795325112867 -0.3826834323650899, -0.9807852804032304 -0.1950903220161286, -1 0, -0.9807852804032304 0.1950903220161284, -0.9238795325112868 0.3826834323650897, -0.8314696123025455 0.555570233019602, -0.7071067811865476 0.7071067811865476, -0.5555702330196022 0.8314696123025452, -0.3826834323650903 0.9238795325112865, -0.1950903220161287 0.9807852804032303, 0 1, 0.1950903220161283 0.9807852804032304, 0.38268343236509 0.9238795325112866, 0.5555702330196018 0.8314696123025455, 0.7071067811865476 0.7071067811865476, 0.8314696123025452 0.5555702330196022, 0.9238795325112865 0.3826834323650904, 0.9807852804032303 0.1950903220161287, 1 0))


Other changes in this PR use angle constants a bit more consistently, an also share Angle::sinCos in other places, where both calculations are performed.

Note that an alternative proposal for a much simpler Angle::sinCos would only zap small values from std::sin and std::cos to zero, rather than a lookup approach on 45 degrees of a unit circle currently in this PR. (Update: that's actually the one that ultimately ended up in this PR.)

@mwtoews mwtoews changed the title Add Angle::SinCos to find exact angles on unit circle, avoid small errors Add Angle::SinCos to avoid small errors, e.g. with buffer operations Nov 2, 2023
@pramsey
Copy link
Member

pramsey commented Nov 2, 2023

Are these called infrequently? Is there any advantage to inlining?

@mwtoews
Copy link
Contributor Author

mwtoews commented Nov 2, 2023

A typical point buffer will call these trig functions 33 times. I'll investigate inlining... (I'm still novice with C++)

@mwtoews
Copy link
Contributor Author

mwtoews commented Nov 5, 2023

I've been checking the performance of this PR with mpgrid.txt (a 20x20 MultiPoint grid):

./bin/geosop -a mpgrid.txt buffer 0.1 -v -r 100 > /dev/null

and on my oldish laptop the results show a minor slowdown, e.g. 810,667 usec with main vs. 822,922 usec with this PR.

I'm unsure how to restructure this PR as an inline function to check the function overhead. Tips welcome!

@dbaston
Copy link
Member

dbaston commented Nov 6, 2023

Hi @mwtoews , you can move the implementation of Angle::sinCos under its declaration in Angle.h to inline it.

@mwtoews
Copy link
Contributor Author

mwtoews commented Nov 6, 2023

Thanks @dbaston that was helpful! I've not pushed the inline version here, but can if needed.

Running the same multipoint grid example shown above, I'm seeing pretty much the same range of benchmarks for three groups: (1) from main, (2) this PR and (3) an inline version of this PR. There isn't any clear outlier. I'm inclined to stick with the current form. Further review or ideas are welcome.

@mwtoews
Copy link
Contributor Author

mwtoews commented Nov 6, 2023

Another idea is to cache 32 pairs of results. Running ctest and capturing the frequency of angles shows some prominent hits, Here is the number of counts for the top 200 hits:

Rank, count, angle, turns
  1. count 1495 at -90.0° or -1/4 turns
  2. count 1453 at 90.0° or 1/4 turns
  3. count 1375 at -45.0° or -1/8 turns
  4. count 1372 at -112.5° or -5/16 turns
  5. count 1371 at -101.25° or -9/32 turns
  6. count 1371 at -123.75° or -11/32 turns
  7. count 1370 at -135.0° or -3/8 turns
  8. count 1368 at -67.5° or -3/16 turns
  9. count 1366 at -157.5° or -7/16 turns
  10. count 1364 at -56.25° or -5/32 turns
  11. count 1364 at -78.75° or -7/32 turns
  12. count 1363 at -146.25° or -13/32 turns
  13. count 1363 at -168.75° or -15/32 turns
  14. count 1362 at -180.0° or -1/2 turns
  15. count 1360 at -22.5° or -1/16 turns
  16. count 1358 at -11.25° or -1/32 turns
  17. count 1358 at -33.75° or -3/32 turns
  18. count 1289 at 135.0° or 3/8 turns
  19. count 1277 at 112.5° or 5/16 turns
  20. count 1275 at 123.75° or 11/32 turns
  21. count 1275 at 101.25° or 9/32 turns
  22. count 1266 at 157.5° or 7/16 turns
  23. count 1264 at 168.75° or 15/32 turns
  24. count 1264 at 146.25° or 13/32 turns
  25. count 1247 at 0.0° or 0 turns
  26. count 1236 at 45.0° or 1/8 turns
  27. count 1224 at 67.5° or 3/16 turns
  28. count 1222 at 78.75° or 7/32 turns
  29. count 1222 at 56.25° or 5/32 turns
  30. count 1222 at 33.75° or 3/32 turns
  31. count 1222 at 22.5° or 1/16 turns
  32. count 1222 at 11.25° or 1/32 turns
  33. count 68 at 180.0° or 1/2 turns
  34. count 55 at 225.0° or 5/8 turns
  35. count 52 at 270.0° or 3/4 turns
  36. count 50 at 213.75° or 19/32 turns
  37. count 50 at 202.5° or 9/16 turns
  38. count 50 at 191.25° or 17/32 turns
  39. count 45 at 258.75° or 23/32 turns
  40. count 45 at 247.5° or 11/16 turns
  41. count 45 at 236.25° or 21/32 turns
  42. count 41 at -57.3808° or -5742669988849355/36028797018963968 turns
  43. count 40 at 3.57633° or 5726705228570285/576460752303423488 turns
  44. count 37 at 59.3493° or 1484919363137221/9007199254740992 turns
  45. count 31 at 68.1986° or 3412657661635439/18014398509481984 turns
  46. count 30 at -202.5° or -9/16 turns
  47. count 29 at -191.25° or -17/32 turns
  48. count 29 at -213.75° or -19/32 turns
  49. count 28 at 60.2551° or 3015164954524133/18014398509481984 turns
  50. count 27 at -68.2538° or -3415419869406893/18014398509481984 turns
  51. count 27 at -79.1269° or -3959509748388695/18014398509481984 turns
  52. count 27 at -120.651° or -1509343885116327/4503599627370496 turns
  53. count 26 at -119.745° or -5992039304216445/18014398509481984 turns
  54. count 25 at 123.69° or 1547361771970713/4503599627370496 turns
  55. count 24 at 183.576° or 2296535569983795/4503599627370496 turns
  56. count 22 at -114.407° or -5724925806317515/18014398509481984 turns
  57. count 22 at 77.7965° or 7785873075793973/36028797018963968 turns
  58. count 22 at 65.593° or 6564546896846955/36028797018963968 turns
  59. count 22 at 53.3894° or 1335802677475191/9007199254740992 turns
  60. count 22 at 7.12502° or 5704575540801303/288230376151711744 turns
  61. count 21 at -270.0° or -3/4 turns
  62. count 21 at -102.204° or -5114287736841935/18014398509481984 turns
  63. count 21 at -126.611° or -3167806957894477/9007199254740992 turns
  64. count 20 at 122.619° or 6135854252317143/18014398509481984 turns
  65. count 20 at 145.491° or 3640184518809783/9007199254740992 turns
  66. count 20 at -17.354° or -6947152705190007/144115188075855872 turns
  67. count 20 at -72.2553° or -7231309825678743/36028797018963968 turns
  68. count 20 at -136.848° or -6847873353404419/18014398509481984 turns
  69. count 19 at -225.0° or -5/8 turns
  70. count 19 at 112.834° or 5646212892830251/18014398509481984 turns
  71. count 18 at -315.0° or -7/8 turns
  72. count 18 at -247.5° or -11/16 turns
  73. count 18 at 43.1524° or 8637383669339671/72057594037927936 turns
  74. count 17 at -292.5° or -13/16 turns
  75. count 17 at -337.5° or -15/16 turns
  76. count 17 at -236.25° or -21/32 turns
  77. count 17 at -258.75° or -23/32 turns
  78. count 17 at -66.6148° or -6666808632385781/36028797018963968 turns
  79. count 17 at -111.801° or -5594521577107209/18014398509481984 turns
  80. count 17 at -138.814° or -6946251985264533/18014398509481984 turns
  81. count 16 at -281.25° or -25/32 turns
  82. count 16 at -303.75° or -27/32 turns
  83. count 16 at -326.25° or -29/32 turns
  84. count 16 at -348.75° or -31/32 turns
  85. count 16 at 79.4642° or 1988194119496081/9007199254740992 turns
  86. count 16 at 68.9283° or 6898343693228483/36028797018963968 turns
  87. count 16 at 101.417° or 5074906260100373/18014398509481984 turns
  88. count 16 at 58.3925° or 5843920916471815/36028797018963968 turns
  89. count 16 at 144.293° or 3610210561289839/9007199254740992 turns
  90. count 16 at -67.1663° or -6722002747818999/36028797018963968 turns
  91. count 15 at 41.1859° or 8243769061907489/72057594037927936 turns
  92. count 15 at 42.1844° or 8443628805371019/72057594037927936 turns
  93. count 15 at 5.35583° or 8576182752803457/576460752303423488 turns
  94. count 15 at -78.3074° or -3918501971781693/18014398509481984 turns
  95. count 14 at -56.3099° or -5635494325717109/36028797018963968 turns
  96. count 14 at 79.7831° or 7984691987343623/36028797018963968 turns
  97. count 14 at 69.5662° or 6962184719946253/36028797018963968 turns
  98. count 14 at 112.46° or 5627497934378733/18014398509481984 turns
  99. count 14 at 101.23° or 5065548780874615/18014398509481984 turns
  100. count 14 at -137.816° or -6896312069396581/18014398509481984 turns
  101. count 14 at -81.1277° or -4059629772104727/18014398509481984 turns
  102. count 13 at 26.5651° or 2658634988023555/36028797018963968 turns
  103. count 13 at 10.6197° or 4251277952247685/144115188075855872 turns
  104. count 13 at -100.217° or -5014858265068767/18014398509481984 turns
  105. count 13 at -110.434° or -5526116902767037/18014398509481984 turns
  106. count 12 at 50.2676° or 2515390495875657/18014398509481984 turns
  107. count 12 at -169.38° or -8475774498711273/18014398509481984 turns
  108. count 12 at -121.608° or -3042631908251507/9007199254740992 turns
  109. count 12 at -78.5832° or -3932303002639791/18014398509481984 turns
  110. count 11 at -132.322° or -3310696166071771/9007199254740992 turns
  111. count 11 at 79.0993° or 7916257289005935/36028797018963968 turns
  112. count 11 at 111.746° or 5591769377334927/18014398509481984 turns
  113. count 11 at 100.873° or 630960562794089/2251799813685248 turns
  114. count 11 at -129.732° or -6491788742866991/18014398509481984 turns
  115. count 10 at 51.3402° or 5138126790869483/36028797018963968 turns
  116. count 10 at 80.085° or 2003726534210923/9007199254740992 turns
  117. count 10 at 70.1701° or 3511311513472783/18014398509481984 turns
  118. count 10 at 132.436° or 6627096891671545/18014398509481984 turns
  119. count 10 at 162.646° or 508675322912015/1125899906842624 turns
  120. count 10 at 113.385° or 2836892465274465/9007199254740992 turns
  121. count 10 at -35.7067° or -3573526240325113/36028797018963968 turns
  122. count 10 at -67.54° or -3379701320362259/18014398509481984 turns
  123. count 10 at -78.77° or -3941650473866377/18014398509481984 turns
  124. count 10 at -100.536° or -2515410511874001/9007199254740992 turns
  125. count 10 at -111.072° or -1389510605031377/4503599627370496 turns
  126. count 10 at -127.972° or -400232396884623/1125899906842624 turns
  127. count 10 at -100.901° or -2524542811118391/9007199254740992 turns
  128. count 10 at -100.803° or -2522090851321267/9007199254740992 turns
  129. count 10 at -111.606° or -5584763777914573/18014398509481984 turns
  130. count 10 at -122.409° or -1531336463296653/4503599627370496 turns
  131. count 10 at -133.212° or -3332963964229325/9007199254740992 turns
  132. count 10 at -144.015° or -28150429702073/70368744177664 turns
  133. count 10 at -154.818° or -968386509875341/2251799813685248 turns
  134. count 10 at -165.621° or -4143837077137383/9007199254740992 turns
  135. count 9 at 315.0° or 7/8 turns
  136. count 9 at 79.197° or 7926035104196915/36028797018963968 turns
  137. count 9 at 68.3941° or 6844880961652009/36028797018963968 turns
  138. count 9 at 57.5911° or 1440929202776983/9007199254740992 turns
  139. count 9 at 46.7882° or 1170640667140757/9007199254740992 turns
  140. count 9 at 35.9852° or 7202797036037901/72057594037927936 turns
  141. count 9 at 25.1823° or 5040488750948091/72057594037927936 turns
  142. count 9 at 14.3793° or 5756320899719873/144115188075855872 turns
  143. count 9 at 190.62° or 4769312005385355/9007199254740992 turns
  144. count 9 at 143.13° or 7162224607394879/18014398509481984 turns
  145. count 9 at -112.306° or -5619791775016343/18014398509481984 turns
  146. count 9 at -99.915° or -2499873093159573/9007199254740992 turns
  147. count 9 at -109.83° or -5495892745267795/18014398509481984 turns
  148. count 9 at 187.125° or 2340933556310289/4503599627370496 turns
  149. count 9 at -34.5085° or -6907220788493989/72057594037927936 turns
  150. count 9 at -100.581° or -5033072823561687/18014398509481984 turns
  151. count 9 at -111.161° or -5562495979757019/18014398509481984 turns
  152. count 9 at -121.742° or -3045984587974105/9007199254740992 turns
  153. count 9 at -142.903° or -7150865528334733/18014398509481984 turns
  154. count 9 at -153.483° or -7680288684530065/18014398509481984 turns
  155. count 9 at -164.064° or -1026220235090157/2251799813685248 turns
  156. count 8 at 32.0054° or 100096880217947/1125899906842624 turns
  157. count 8 at 89.8461° or 4495898472007693/18014398509481984 turns
  158. count 8 at -29.7449° or -1488434672846363/18014398509481984 turns
  159. count 8 at -40.7004° or -4073295694974003/36028797018963968 turns
  160. count 8 at -51.6558° or -1292428009064027/9007199254740992 turns
  161. count 8 at -62.6113° or -6266138385537385/36028797018963968 turns
  162. count 8 at -73.5668° or -3681282367409331/18014398509481984 turns
  163. count 8 at -84.5223° or -2114747771024985/9007199254740992 turns
  164. count 8 at -95.4777° or -2388851856345511/9007199254740992 turns
  165. count 8 at -106.433° or -5325906879332489/18014398509481984 turns
  166. count 8 at -117.389° or -5874145073971057/18014398509481984 turns
  167. count 8 at -128.344° or -3211166614306883/9007199254740992 turns
  168. count 8 at -139.3° or -6970571423252335/18014398509481984 turns
  169. count 8 at 209.745° or 5247819465793471/9007199254740992 turns
  170. count 8 at 171.198° or 4283373605592079/9007199254740992 turns
  171. count 8 at 159.948° or 4001898628881423/9007199254740992 turns
  172. count 8 at 148.698° or 3720423652170767/9007199254740992 turns
  173. count 8 at 137.448° or 3438948675460111/9007199254740992 turns
  174. count 8 at 126.198° or 6314947397498909/18014398509481984 turns
  175. count 8 at 114.948° or 5751997444077597/18014398509481984 turns
  176. count 8 at 103.698° or 5189047490656285/18014398509481984 turns
  177. count 8 at 92.4481° or 4626102541234559/18014398509481984 turns
  178. count 8 at 81.1981° or 8126305175626495/36028797018963968 turns
  179. count 8 at 69.9481° or 7000405268783871/36028797018963968 turns
  180. count 8 at 58.6981° or 5874505361941247/36028797018963968 turns
  181. count 8 at 47.4481° or 4748605455098623/36028797018963968 turns
  182. count 8 at 36.1981° or 7245411096511997/72057594037927936 turns
  183. count 8 at 24.9481° or 2496805641413375/36028797018963968 turns
  184. count 8 at 13.6981° or 1370905734570751/36028797018963968 turns
  185. count 8 at 2.44812° or 7840250538494761/1152921504606846976 turns
  186. count 8 at -8.80188° or -3523568310058651/144115188075855872 turns
  187. count 8 at 77.539° or 3880051238963121/18014398509481984 turns
  188. count 8 at 65.078° or 1628251425277873/9007199254740992 turns
  189. count 8 at 52.617° or 2632954462148371/18014398509481984 turns
  190. count 8 at 77.3428° or 3870233391775453/18014398509481984 turns
  191. count 8 at 64.6856° or 6473734312360821/36028797018963968 turns
  192. count 8 at 140.906° or 881366957075371/2251799813685248 turns
  193. count 8 at 78.8469° or 7890997099095973/36028797018963968 turns
  194. count 8 at -15.7086° or -6288466231689971/144115188075855872 turns
  195. count 8 at -8.12071° or -812720589555197/36028797018963968 turns
  196. count 8 at -19.8177° or -3966710503792901/72057594037927936 turns
  197. count 8 at -31.5148° or -6308001846073587/72057594037927936 turns
  198. count 8 at -43.2118° or -1081159146544491/9007199254740992 turns
  199. count 8 at -54.9089° or -1373820564329577/9007199254740992 turns
  200. count 8 at -66.6059° or -6665917920459479/36028797018963968 turns

However, I'm not sure how to implement this caching/lookup. Or if it would make single point buffers slower. It might be more worthwhile adding more prepared trig values beyond 0, 45, 90 etc. already implemented, e.g. exact trigonometric values.

Looking at the results a bit more, it could make sense to try adding a static const lookup table to Angle.h (similar to (src/deps/ryu/d2s_full_table.h) with 32 values that are each 1/32 of a turn.

@pramsey
Copy link
Member

pramsey commented Nov 6, 2023

Don't polish the 💩 too much, you've already shown that it's only a small performance change, hard to measure, I think know that it's a super critical path, so perhaps just go ahead and commit.

@mwtoews
Copy link
Contributor Author

mwtoews commented Nov 7, 2023

Ok, I think I'm done the pedanticness. I tried the look-up-of-32-angles approach and didn't see any performance gains, so I'm going with a much simpler inline approach to just clip trig values within 5e-16 to zero. The performance is the same.

@mwtoews mwtoews merged commit e4a6e1a into libgeos:main Nov 7, 2023
28 checks passed
@mwtoews mwtoews deleted the clipsincos branch November 7, 2023 21:05
@mwtoews mwtoews changed the title Add Angle::SinCos to avoid small errors, e.g. with buffer operations Add Angle::sinCos to avoid small errors, e.g. with buffer operations Nov 8, 2023
@mwtoews mwtoews changed the title Add Angle::sinCos to avoid small errors, e.g. with buffer operations Add Angle::sinCosSnap to avoid small errors, e.g. with buffer operations Nov 16, 2023
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.

None yet

3 participants