Add filter to force preservation of precision in floating point values (

#32)

* add filter to force double precision

* add tests for double precision

* add note on docs page about precision

docs/index.md

@@ -43,7 +43,7 @@ $ pytest
 
 ## Bridging the IDL-Python Gap the Bad Way...
 
-hissw is a hack, albeit a clever one. In general, the methodology employed here (round-tripping everything in memory to disk) is a *terrible idea*. There's no shared memory between Python and IDL or anything fancy like that. If you're interested in something more complicated (and harder to install), you may want to check out the more official [Python-to-IDL bridge](https://www.harrisgeospatial.com/docs/Python.html).
+`hissw` is a hack, albeit a clever one. In general, the methodology employed here (round-tripping everything in memory to disk) is a *terrible idea*. There's no shared memory between Python and IDL or anything fancy like that. If you're interested in something more complicated (and harder to install), you may want to check out the more official [Python-to-IDL bridge](https://www.harrisgeospatial.com/docs/Python.html).
 
 ## Word(s) of Caution
 
@@ -54,3 +54,18 @@ hissw is a hack, albeit a clever one. In general, the methodology employed here
 * Widgets and plotting will (likely) **not** work
 * This has **not** been tested extensively against all SSW/IDL functionality. There are likely many cases where hissw will not work. [Bug reports](https://github.com/wtbarnes/hissw/issues) and [pull requests](https://github.com/wtbarnes/hissw/pulls) welcome!
 
+## A Note on Preserving Precision between Python and IDL
+
+`hissw` relies on string representations of Python objects when inserting values into the resulting IDL script.
+However, [the default string representation of Python floats is truncated well below the actual floating point precision.](https://docs.python.org/3/tutorial/floatingpoint.html#floating-point-arithmetic-issues-and-limitations)
+This can result in a loss of precision when injecting floating point values into IDL from Python as documented in [this issue](https://github.com/wtbarnes/hissw/issues/31).
+Thus, to ensure that floating point values are accurately represented, hissw provides the `force_double_precision` filter.
+This filter can be used as follows, where `var` is a (scalar or array) variable containing a floating point value,
+
+```IDL
+var = {{ var | force_double_precision }}
+```
+
+This filter uses the builtin [`float.as_integer_ratio`](https://docs.python.org/3/library/stdtypes.html#float.as_integer_ratio) method to represent floating point values as division operations between integer values
+to ensure that precision is not lost in the string representation Python values in the resulting IDL script.
+If used in combination with other filters, this filter should be used last.

hissw/environment.py

@@ -59,6 +59,7 @@ def __init__(self, ssw_packages=None, ssw_paths=None, extra_paths=None,
         self.env.filters['to_unit'] = units_filter
         self.env.filters['log10'] = log10_filter
         self.env.filters['string_list'] = string_list_filter
+        self.env.filters['force_double_precision'] = force_double_precision_filter
         if filters is not None:
             for k, v in filters.items():
                 self.env.filters[k] = v

hissw/filters.py

@@ -3,7 +3,10 @@
 """
 import numpy as np
 
-__all__ = ['units_filter', 'log10_filter']
+__all__ = ['units_filter',
+           'log10_filter',
+           'string_list_filter',
+           'force_double_precision_filter',]
 
 
 def units_filter(quantity, unit):
@@ -31,3 +34,29 @@ def string_list_filter(string_list):
     in the list will not be quoted when passed into the template.
     """
     return [f"'{s}'" for s in string_list]
+
+
+def force_double_precision_filter(value):
+    """
+    Force a number (or array of numbers) to have double precision in IDL.
+
+    hissw relies on string representations of Python objects when inserting
+    values into the resulting IDL script. However, the default string representation
+    of Python floats is truncated well below the actual floating point precision.
+    See https://docs.python.org/3/tutorial/floatingpoint.html#floating-point-arithmetic-issues-and-limitations.
+    Thus, to preserve precision, this filter uses the `as_integer_ratio` method
+    to represent floating point values as division operations between integer values
+    to ensure that precision is preserved when passing floating point values from
+    IDL into Python. When used in combination with other filters, this filter should
+    be used last as it forces a conversion to a list of strings.
+    """
+    if isinstance(value, (np.ndarray, list)):
+        str_list = [force_double_precision_filter(x) for x in value]
+        # NOTE: this has to be done manually because each entry is formatted as a
+        # string such that the division is not evaluated in Python. However, we
+        # want this to be inserted as an array of integer divison operations.
+        return f"[{','.join(str_list)}]"
+    else:
+        # If it is neither an array or list, 
+        a, b = value.as_integer_ratio()
+        return f'{a}d / {b}d'

hissw/tests/test_hissw.py

@@ -176,3 +176,42 @@ def test_custom_header_footer(idl_home):
     result = env_custom.run(script, args=args)
     assert result['foo'] == args['a']
     assert result['bar'] == args['a'] + args['b']
+
+
+@pytest.mark.parametrize('var', [
+    1 / 3,
+    float(10),
+    [1/3, 1/7, 0.1],
+    np.random.rand(100),
+    list(np.random.rand(10).astype(np.longdouble))
+])
+def test_force_double_precision_filter(var, idl_env):
+    # This test ensures that floating point precision is conserved when passing
+    # things from Python to IDL. See https://github.com/wtbarnes/hissw/issues/31
+    result = idl_env.run(
+            '''
+            var = {{ var | force_double_precision }}
+            var_size = size(var)
+            ''',
+            args={'var': var})
+    assert u.allclose(var, result['var'], atol=0.0, rtol=0.0)
+    # The result of IDL size has a variable number of entries depending on the
+    # dimensionality of the input array
+    # NOTE: 5 corresponds to double-precision floating point. See
+    # https://www.l3harrisgeospatial.com/docs/make_array.html#TYPE
+    assert result['var_size'][len(result['var'].shape) + 1] == 5
+
+
+@pytest.mark.parametrize('var', [
+    1 / 3 * u.s,
+    float(10) * u.s,
+    [1/3, 1/7, 0.1] * u.s,
+    np.random.rand(100) * u.minute,
+])
+def test_force_double_precision_filter_with_quantity(var, idl_env):
+    # This test ensures that floating point precision is conserved when passing
+    # things from Python to IDL. See https://github.com/wtbarnes/hissw/issues/31
+    result = idl_env.run('var = {{ var | to_unit("h") | force_double_precision }}',
+                         args={'var': var})
+    assert u.allclose(var.to_value('h'), result['var'], atol=0.0, rtol=0.0)
+