Repository: biotite-dev/biotite

Issue:
`pdbx.compress()` may encounter integer overflows
If `BinaryCIFData` with float values given to `pdbx.compress()` contains values ranging over multiple orders of magnitude the compressed data may have false values due to integer underflows/overflows.
The reason is that the factor for the `FixedPointEncoding` is chosen so large that the large values in the array do not fit into `int32` anymore, and thus overflow.

This should be fixed in two locations: `compress()` should more carefully pick the factor and `FixedPointEncoding.encode` should raise an exception instead of only raising a warning and overflowing.

Discussion:

Solution patch:
diff --git a/src/biotite/structure/io/pdbx/compress.py b/src/biotite/structure/io/pdbx/compress.py
index 0f2e8d3d..52fefc6d 100644
--- a/src/biotite/structure/io/pdbx/compress.py
+++ b/src/biotite/structure/io/pdbx/compress.py
@@ -3,6 +3,7 @@ __name__ = "biotite.structure.io.pdbx"
 __author__ = "Patrick Kunzmann"
 
 import itertools
+import warnings
 import msgpack
 import numpy as np
 import biotite.structure.io.pdbx.bcif as bcif
@@ -17,7 +18,7 @@ from biotite.structure.io.pdbx.encoding import (
 )
 
 
-def compress(data, float_tolerance=1e-6):
+def compress(data, float_tolerance=None, rtol=1e-6, atol=1e-4):
     """
     Try to reduce the size of a *BinaryCIF* file (or block, category, etc.) by testing
     different data encodings for each data array and selecting the one, which results in
@@ -29,6 +30,12 @@ def compress(data, float_tolerance=1e-6):
         The data to compress.
     float_tolerance : float, optional
         The relative error that is accepted when compressing floating point numbers.
+        DEPRECATED: Use `rtol` instead.
+    rtol, atol : float, optional
+        The compression factor of floating point numbers is chosen such that
+        either the relative (`rtol`) or absolute (`atol`) tolerance is fulfilled
+        for each value, i.e. the difference between the compressed and uncompressed
+        value is smaller than the tolerance.
 
     Returns
     -------
@@ -58,55 +65,61 @@ def compress(data, float_tolerance=1e-6):
     >>> print(f"{len(compressed_file.read()) // 1000} KB")
     111 KB
     """
+    if float_tolerance is not None:
+        warnings.warn(
+            "The 'float_tolerance' parameter is deprecated, use 'rtol' instead",
+            DeprecationWarning,
+        )
+
     match type(data):
         case bcif.BinaryCIFFile:
-            return _compress_file(data, float_tolerance)
+            return _compress_file(data, rtol, atol)
         case bcif.BinaryCIFBlock:
-            return _compress_block(data, float_tolerance)
+            return _compress_block(data, rtol, atol)
         case bcif.BinaryCIFCategory:
-            return _compress_category(data, float_tolerance)
+            return _compress_category(data, rtol, atol)
         case bcif.BinaryCIFColumn:
-            return _compress_column(data, float_tolerance)
+            return _compress_column(data, rtol, atol)
         case bcif.BinaryCIFData:
-            return _compress_data(data, float_tolerance)
+            return _compress_data(data, rtol, atol)
         case _:
             raise TypeError(f"Unsupported type {type(data).__name__}")
 
 
-def _compress_file(bcif_file, float_tolerance):
+def _compress_file(bcif_file, rtol, atol):
     compressed_file = bcif.BinaryCIFFile()
     for block_name, bcif_block in bcif_file.items():
-        compressed_block = _compress_block(bcif_block, float_tolerance)
+        compressed_block = _compress_block(bcif_block, rtol, atol)
         compressed_file[block_name] = compressed_block
     return compressed_file
 
 
-def _compress_block(bcif_block, float_tolerance):
+def _compress_block(bcif_block, rtol, atol):
     compressed_block = bcif.BinaryCIFBlock()
     for category_name, bcif_category in bcif_block.items():
-        compressed_category = _compress_category(bcif_category, float_tolerance)
+        compressed_category = _compress_category(bcif_category, rtol, atol)
         compressed_block[category_name] = compressed_category
     return compressed_block
 
 
-def _compress_category(bcif_category, float_tolerance):
+def _compress_category(bcif_category, rtol, atol):
     compressed_category = bcif.BinaryCIFCategory()
     for column_name, bcif_column in bcif_category.items():
-        compressed_column = _compress_column(bcif_column, float_tolerance)
+        compressed_column = _compress_column(bcif_column, rtol, atol)
         compressed_category[column_name] = compressed_column
     return compressed_category
 
 
-def _compress_column(bcif_column, float_tolerance):
-    data = _compress_data(bcif_column.data, float_tolerance)
+def _compress_column(bcif_column, rtol, atol):
+    data = _compress_data(bcif_column.data, rtol, atol)
     if bcif_column.mask is not None:
-        mask = _compress_data(bcif_column.mask, float_tolerance)
+        mask = _compress_data(bcif_column.mask, rtol, atol)
     else:
         mask = None
     return bcif.BinaryCIFColumn(data, mask)
 
 
-def _compress_data(bcif_data, float_tolerance):
+def _compress_data(bcif_data, rtol, atol):
     array = bcif_data.array
     if len(array) == 1:
         # No need to compress a single value -> Use default uncompressed encoding
@@ -124,15 +137,23 @@ def _compress_data(bcif_data, float_tolerance):
 
     elif np.issubdtype(array.dtype, np.floating):
         to_integer_encoding = FixedPointEncoding(
-            10 ** _get_decimal_places(array, float_tolerance)
+            10 ** _get_decimal_places(array, rtol, atol)
         )
-        integer_array = to_integer_encoding.encode(array)
-        best_encoding, size_compressed = _find_best_integer_compression(integer_array)
-        if size_compressed < _data_size_in_file(bcif.BinaryCIFData(array)):
-            return bcif.BinaryCIFData(array, [to_integer_encoding] + best_encoding)
-        else:
-            # The float array is smaller -> encode it directly as bytes
+        try:
+            integer_array = to_integer_encoding.encode(array)
+        except ValueError:
+            # With the given tolerances integer underflow/overflow would occur
+            # -> do not use integer encoding
             return bcif.BinaryCIFData(array, [ByteArrayEncoding()])
+        else:
+            best_encoding, size_compressed = _find_best_integer_compression(
+                integer_array
+            )
+            if size_compressed < _data_size_in_file(bcif.BinaryCIFData(array)):
+                return bcif.BinaryCIFData(array, [to_integer_encoding] + best_encoding)
+            else:
+                # The float array is smaller -> encode it directly as bytes
+                return bcif.BinaryCIFData(array, [ByteArrayEncoding()])
 
     elif np.issubdtype(array.dtype, np.integer):
         array = _to_smallest_integer_type(array)
@@ -273,7 +294,7 @@ def _data_size_in_file(data):
     return len(bytes_in_file)
 
 
-def _get_decimal_places(array, tol):
+def _get_decimal_places(array, rtol, atol):
     """
     Get the number of decimal places in a floating point array.
 
@@ -281,21 +302,25 @@ def _get_decimal_places(array, tol):
     ----------
     array : numpy.ndarray
         The array to analyze.
-    tol : float, optional
-        The relative tolerance allowed when the values are cut off after the returned
-        number of decimal places.
+    rtol, atol : float, optional
+        The relative and absolute tolerance allowed when the values are cut off after
+        the returned number of decimal places.
 
     Returns
     -------
     decimals : int
         The number of decimal places.
     """
+    if rtol <= 0 and atol <= 0:
+        raise ValueError("At least one of 'rtol' and 'atol' must be greater than 0")
     # Decimals of NaN or infinite values do not make sense
     # and 0 would give NaN when rounding on decimals
     array = array[np.isfinite(array) & (array != 0)]
-    for decimals in itertools.count(start=-_order_magnitude(array)):
+    for decimals in itertools.count(start=min(0, -_order_magnitude(array))):
         error = np.abs(np.round(array, decimals) - array)
-        if np.all(error < tol * np.abs(array)):
+        if decimals == 100:
+            raise
+        if np.all((error < rtol * np.abs(array)) | (error < atol)):
             return decimals
 
 
diff --git a/src/biotite/structure/io/pdbx/encoding.pyx b/src/biotite/structure/io/pdbx/encoding.pyx
index 55d5d51a..af7d3a96 100644
--- a/src/biotite/structure/io/pdbx/encoding.pyx
+++ b/src/biotite/structure/io/pdbx/encoding.pyx
@@ -325,7 +325,16 @@ class FixedPointEncoding(Encoding):
                 )
 
         # Round to avoid wrong values due to floating point inaccuracies
-        return np.round(data * self.factor).astype(np.int32)
+        scaled_data = np.round(data * self.factor)
+        # Check if an integer underflow/overflow would occur during conversion
+        if (
+            np.max(scaled_data) > np.iinfo(np.int32).max or
+            np.min(scaled_data) < np.iinfo(np.int32).min
+        ):  # fmt: skip
+            raise ValueError(
+                "With the given factor some values get too large for integer encoding"
+            )
+        return scaled_data.astype(np.int32)
 
     def decode(self, data):
         return (data / self.factor).astype(


Test patch:
diff --git a/tests/structure/io/test_pdbx.py b/tests/structure/io/test_pdbx.py
index 262a03a6..97e3686b 100644
--- a/tests/structure/io/test_pdbx.py
+++ b/tests/structure/io/test_pdbx.py
@@ -1000,15 +1000,39 @@ def test_compress_file(path):
     assert _file_size(compressed_file) <= _file_size(orig_file)
 
 
+@pytest.mark.parametrize("value", [1e10, 1e-10])
+def test_extreme_float_compression(value):
+    """
+    Check if :func:`compress()` correctly falls back to direct byte encoding of floats
+    in extreme cases where fixed point encoding would lead to integer
+    underflow/overflow.
+    """
+    # Not only very small/large values, but a large difference between the values are
+    # required, to make fixed point encoding fail
+    ref_array = np.array([value, 1.0])
+
+    compressed_data = pdbx.compress(pdbx.BinaryCIFData(ref_array), atol=0)
+    serialized_compressed_data = compressed_data.serialize()
+    data = pdbx.BinaryCIFData.deserialize(serialized_compressed_data)
+
+    # Check that no fixed point encoding was used
+    assert len(data.encoding) == 1
+    assert type(data.encoding[0]) is pdbx.ByteArrayEncoding
+    assert np.all(data.array == ref_array)
+
+
 @pytest.mark.parametrize(
     "number, ref_decimals",
     [
         (1.0, 0),
-        (1.2345, 4),
-        (0.00012345, 8),
+        (1.23, 2),
+        (0.001, 3),
+        (0.0012345, 4),
         (12300, -2),
-        (123.456, 3),
-        (123.0000000001, 0),
+        (123.45, 2),
+        (123.45678, 4),
+        (123.00001, 0),
+        (0.00001, 0),
         (0.0, 0),
     ],
 )
@@ -1017,7 +1041,7 @@ def test_decimal_places(number, ref_decimals):
     Check if :func`:_get_decimal_places()` returns the correct number of decimal places
     for known examples.
     """
-    test_decimals = get_decimal_places(np.array([number]), 1e-6)
+    test_decimals = get_decimal_places(np.array([number]), 1e-6, 1e-4)
     assert test_decimals == ref_decimals