Merge pull request #1578 from phargogh/bugfix/1562-tests-failing-on-m…

…ain-viewshed

Correcting viewshed behavior on M1 with a compiler flag

setup.py

@@ -2,8 +2,8 @@
 import platform
 import subprocess
 
-from Cython.Build import cythonize
 import numpy
+from Cython.Build import cythonize
 from setuptools import setup
 from setuptools.command.build_py import build_py as _build_py
 from setuptools.extension import Extension
@@ -50,17 +50,26 @@ def run(self):
         Extension(
             name=f'natcap.invest.{package}.{module}',
             sources=[f'src/natcap/invest/{package}/{module}.pyx'],
-            extra_compile_args=compiler_and_linker_args,
+            extra_compile_args=compiler_args + compiler_and_linker_args,
             extra_link_args=compiler_and_linker_args,
             language='c++',
             define_macros=[("NPY_NO_DEPRECATED_API", "NPY_1_7_API_VERSION")]
-        ) for package, module in [
-            ('delineateit', 'delineateit_core'),
-            ('recreation', 'out_of_core_quadtree'),
-            ('scenic_quality', 'viewshed'),
-            ('ndr', 'ndr_core'),
-            ('sdr', 'sdr_core'),
-            ('seasonal_water_yield', 'seasonal_water_yield_core')
+        ) for package, module, compiler_args in [
+            ('delineateit', 'delineateit_core', []),
+            ('recreation', 'out_of_core_quadtree', []),
+            # clang-14 defaults to -ffp-contract=on, which causes the
+            # arithmetic of A*B+C to be implemented using a contraction, which
+            # causes an unexpected change in the precision in some viewshed
+            # tests on ARM64 (mac M1).  See these issues for more details:
+            #  * https://github.com/llvm/llvm-project/issues/91824
+            #  * https://github.com/natcap/invest/issues/1562
+            #  * https://github.com/natcap/invest/pull/1564/files
+            # Using this flag on gcc and on all versions of clang should work
+            # as expected, with consistent results.
+            ('scenic_quality', 'viewshed', ['-ffp-contract=off']),
+            ('ndr', 'ndr_core', []),
+            ('sdr', 'sdr_core', []),
+            ('seasonal_water_yield', 'seasonal_water_yield_core', [])
         ]
     ], compiler_directives={'language_level': '3'}),
     include_dirs=[numpy.get_include()],

src/natcap/invest/scenic_quality/viewshed.pyx

@@ -813,23 +813,7 @@ def viewshed(dem_raster_path_band,
             if target_distance > max_visible_radius:
                 break
 
-            # This is a weird platform-specific workaround addressing
-            # https://github.com/natcap/invest/issues/1562
-            # On M1 macs, the all-in-one-line addition of _product and r_v
-            # would create small but noticeable numerical error.  Breaking the
-            # calculation onto two lines eliminates the numerical error.  This
-            # behavior is reproducible in C, outside of Cython on an M1 mac.
-            # So, this calculation would introduce error:
-            #   z = (((previous_height-r_v)/slope_distance) * target_distance) + r_v
-            # while the formlation below does not.
-            # For the script used for testing, see
-            #  https://gist.github.com/phargogh/c4264b37e7f0beed31661eacce53d14a
-            #
-            # Some of this may be related to the fact that x86 chips have
-            # extended precision for FPU-based calculations while M1 ARM chips
-            # do not.  Still, that doesn't explain why the error is introduced.
-            _product = (((previous_height-r_v)/slope_distance) * target_distance)
-            z = _product + r_v
+            z = (((previous_height-r_v)/slope_distance) * target_distance) + r_v
 
             # add on refractivity/curvature-of-earth calculations.
             adjustment = 0.0  # increase in required height due to curvature