Codebase audit: fix wrong-result, correctness, cleanup and performance findings

.gitignore

@@ -28,3 +28,9 @@ docs/_build
 /subprojects/*
 !/subprojects/*.wrap
 .junie
+
+# NetCDF files written to the working directory by the test suite
+/angstrom_unit.nc
+/no_physical_sizes.nc
+/no_unit.nc
+/parallel_save_test.nc

SurfaceTopography/ChangeLog.md

@@ -1,6 +1,41 @@
 Change log for SurfaceTopography
 ================================
 
+v1.23.0 (not yet released)
+--------------------------
+
+- BUG: Transposed and downsampled topographies report correct `pixel_size`; affected derivatives and PSDs for anisotropic pixels
+- BUG: Masked data is normalized by the number of defined points; changes `rms_height` and moments of files with undefined data
+- BUG: Power spectrum no longer halves the `q=0` entry when folding the +q and -q branches (#282)
+- BUG: Fixed reshape order for non-square scans in the DI, EZD, MI and PS readers
+- BUG: MI reader converts lateral sizes from meters to the requested unit
+- BUG: BCR reader treats `voidpixels` as a count and masks against format markers
+- BUG: Fixed channel selection in the NMM and OIR readers
+- BUG: Fixed EZD magic check (bytes vs. str, inverted logic)
+- BUG: ZAG container reader re-opens the ZIP file for each read
+- BUG: Rewrote the C++ moment kernels in exact polynomial form (sign errors, missing prefactors, division by zero)
+- BUG: `Bicubic.__call__` handles non-contiguous input arrays
+- BUG: Fixed bounds of the periodic 2D bearing area
+- BUG: `to_uniform()` interpolates at the correct positions for scans not starting at zero
+- BUG: Slope detrending uses the length-weighted slope
+- BUG: Hann window in the nonuniform power spectrum is origin invariant
+- BUG: `fourier_synthesis` zeroes the DC mode instead of setting it to `C(q=1)`
+- BUG: `make_sphere(..., periodic=True)` passes the periodic flag on
+- BUG: Fixed `direction` semantics of `scan_line_align`
+- BUG: Nonuniform `polyfit` solves in scan-centered coordinates; detrending no longer depends on the x origin
+- BUG: Fixed Nyquist handling in Fourier interpolation and for anisotropic pixels in `power_spectrum_from_area`
+- BUG: Registering an analysis function on a subclass no longer modifies the base class
+- BUG: Fixed scale-dependent curvature on log-spaced and reliability-trimmed grids
+- ENH: Bicubic interpolation computes spline coefficients on demand, instead of 128 bytes per pixel up front
+- ENH: Container PSD integration precomputes bandwidth intervals in one pass, instead of re-reading every file per topography
+- ENH: The IBW, MI and MNT readers no longer keep file data in memory after construction
+- ENH: Nonuniform height-height autocorrelation moved to a C++ kernel
+- ENH: Rigid-sphere scan vectorized over blocks of bounded memory
+- ENH: Gaussian process regression reuses a single Cholesky factorization
+- ENH: Flood-fill stack in the patch finder is allocated once per call, not once per patch
+- TST: Added property-based invariant tests for transposition, translation, scaling and origin invariance
+- MAINT: Removed broken `pytest-flake8` from the test dependencies; renamed `tiffile` to `tifffile`
+
 v1.22.0 (08May26)
 -----------------
 

SurfaceTopography/Container/Averaging.py

@@ -34,7 +34,7 @@
 from ..Support.Regression import resample
 from .SurfaceContainer import SurfaceContainer
 
-_log = logging.Logger(__name__)
+_log = logging.getLogger(__name__)
 
 
 def log_average(self, function_name, unit, nb_points_per_decade=10, reliable=True, progress_callback=None, **kwargs):

SurfaceTopography/Container/IO/CE.py

@@ -77,7 +77,18 @@ def __init__(self, zipname, filename):
         self._filename = filename
 
     def __call__(self):
-        return ZipFile(self._zipname, mode="r").open(self._filename, mode="r")
+        zipfile = ZipFile(self._zipname, mode="r")
+        stream = zipfile.open(self._filename, mode="r")
+        # Close the containing ZipFile together with the member stream;
+        # otherwise its file descriptor leaks until garbage collection
+        original_close = stream.close
+
+        def close():
+            original_close()
+            zipfile.close()
+
+        stream.close = close
+        return stream
 
 
 class CEReader(ContainerReaderBase):

SurfaceTopography/Container/IO/ZAG.py

@@ -57,7 +57,46 @@
 from ..SurfaceContainer import LazySurfaceContainer
 from .Reader import ContainerReaderBase
 
-_log = logging.Logger(__name__)
+_log = logging.getLogger(__name__)
+
+
+class ZAGFileOpener(object):
+    """
+    Callable that (re)opens a member of the ZIP archive embedded in a ZAG
+    file. Readers constructed from ZAG containers hold on to this callable
+    instead of an open stream, so that reading topographies still works
+    after `ZAGReader.close` has been called (e.g. when the reader was used
+    as a context manager by `read_container`).
+    """
+
+    def __init__(self, fobj, filename):
+        self._fobj = fobj  # file name or stream
+        self._filename = filename
+
+    def __call__(self):
+        # `ZipFile` automatically skips the ZAG header (magic + BMP
+        # thumbnail) preceding the ZIP archive, since ZIP archives are
+        # located through their central directory at the end of the file.
+        if hasattr(self._fobj, "read"):
+            zipfile = ZipFile(self._fobj, "r")
+            close_files = [zipfile]
+        else:
+            rawfile = open(self._fobj, "rb")
+            zipfile = ZipFile(rawfile, "r")
+            close_files = [zipfile, rawfile]
+        stream = zipfile.open(self._filename, "r")
+        # Close the containing ZipFile (and raw file) together with the
+        # member stream; otherwise file descriptors leak until garbage
+        # collection
+        original_close = stream.close
+
+        def close():
+            original_close()
+            for f in close_files:
+                f.close()
+
+        stream.close = close
+        return stream
 
 
 class ZAGReader(ContainerReaderBase):
@@ -113,6 +152,7 @@ def __init__(self, fobj):
         from ...IO.common import OpenFromAny
 
         # Open if a file name is given
+        self._fobj = fobj
         if not hasattr(fobj, "read"):
             # This is a string
             self._fstream = open(fobj, "rb")
@@ -155,10 +195,16 @@ def __init__(self, fobj):
                             # Lazy import to avoid circular dependency during package initialization
                             from ...IO import ZONReader
                             data_path = data.find(self._PATH_TAG).text
+                            # We pass a callable that reopens the ZIP member
+                            # rather than an open stream: the stream `f`
+                            # is closed when this reader is closed, but the
+                            # lazy container must be able to read
+                            # topographies after that (see `read_container`).
                             readers += [
                                 ZONReader(
-                                    z.open(
-                                        f"{data_uuid}/{data_path}/{self._ZON_UUID}", "r"
+                                    ZAGFileOpener(
+                                        self._fobj,
+                                        f"{data_uuid}/{data_path}/{self._ZON_UUID}",
                                     )
                                 ).topography
                             ]

SurfaceTopography/Container/IO/__init__.py

@@ -68,7 +68,7 @@ def detect_format(fobj):
     raise CannotDetectFileFormat(msg)
 
 
-def open_container(fobj, format=None):
+def open_container(fobj, format=None, **kwargs):
     r"""
     Returns a container reader object for the file `fobj`.
 
@@ -79,6 +79,9 @@ def open_container(fobj, format=None):
     format : str, optional
         Specify in which format the file should be interpreted.
         (Default: None, which means autodetect file format)
+    **kwargs : dict
+        Additional keyword arguments passed to the reader's constructor
+        (e.g. `datafile_keys` or `ignore_filters` of `CEReader`).
 
     Returns
     -------
@@ -87,13 +90,13 @@ def open_container(fobj, format=None):
     """  # noqa: E501
     if not hasattr(fobj, 'read'):  # fobj is a path
         if not os.path.isfile(fobj):
-            raise FileExistsError("file {} not found".format(fobj))
+            raise FileNotFoundError("file {} not found".format(fobj))
 
     if format is None:
         msg = ""
         for reader in readers:
             try:
-                return reader(fobj)
+                return reader(fobj, **kwargs)
             except Exception as err:
                 msg += "tried {}: \n {}\n\n".format(reader.__name__, err)
             finally:
@@ -107,7 +110,7 @@ def open_container(fobj, format=None):
             raise UnknownFileFormat(
                 "{} not in registered container file formats {}".format(
                     fobj, lookup_reader_by_format.keys()))
-        return lookup_reader_by_format[format](fobj)
+        return lookup_reader_by_format[format](fobj, **kwargs)
 
 
 def read_container(fn, format=None, **kwargs):
@@ -116,18 +119,24 @@ def read_container(fn, format=None, **kwargs):
 
     Parameters
     ----------
-    fobj : str or filelike object
+    fn : str or filelike object
         Path of the file or file-like object.
+    format : str, optional
+        Specify in which format the file should be interpreted.
+        (Default: None, which means autodetect file format)
+    **kwargs : dict
+        Additional keyword arguments passed to the reader's constructor
+        (e.g. `datafile_keys` or `ignore_filters` of `CEReader`).
 
     Returns
     -------
     container : list of subclasses of :obj:`SurfaceContainer`
         A list of container objects.
     """
     containers = []
-    with open_container(fn, format=format) as reader:
+    with open_container(fn, format=format, **kwargs) as reader:
         for i in range(reader.nb_containers):
-            containers += [reader.container(index=i, **kwargs)]
+            containers += [reader.container(index=i)]
     return containers
 
 
@@ -147,12 +156,21 @@ def read_published_container(publication_url, **request_args):
     container : SurfaceContainer
         Surface container object read from URL.
     """
+    # Default to a finite timeout so a stale server cannot hang the caller
+    # indefinitely; can be overridden through `request_args`
+    request_args.setdefault('timeout', 60)
+
     # If we send json as a request header, then contact.engineering will response with a JSON dictionary
-    response = requests.get(publication_url, headers={'Accept': 'application/json'})
+    custom_headers = request_args.pop('headers', {})
+    response = requests.get(publication_url,
+                            headers={'Accept': 'application/json', **custom_headers},
+                            **request_args)
+    response.raise_for_status()
     data = response.json()
     download_url = data['download_url']
 
     # Then download and read container
-    container_response = requests.get(download_url, **request_args)
+    container_response = requests.get(download_url, headers=custom_headers, **request_args)
+    container_response.raise_for_status()
     container_file = io.BytesIO(container_response.content)
     return read_container(container_file)

SurfaceTopography/Container/Integration.py

@@ -28,28 +28,32 @@
 from ..Generic.Moments import compute_1d_moment, compute_iso_moment
 
 
-def _bandwidth_count_from_profile(self, qx, unit, reliable=True):
+def _bandwidth_intervals_from_profile(self, unit, reliable=True):
     r"""
-    Return number of topographies that include qx in their bandwidth.
+    Return the wavevector intervals covered by the topographies of this
+    container as a pair of sorted arrays (lower bounds, upper bounds).
+    This requires a single pass over the container; counting queries can
+    then be answered without touching the topographies again (lazy
+    containers read the data file on every element access).
 
     Parameters:
     -----------
     self : SurfaceContainer
         Collection of Height containers
-    qx: np.ndarrau of floats
-        wavevector
     unit : str
         Unit of lengths in which the wavevector is defined.
     reliable : bool, optional
         Only incorporate data deemed reliable. (Default: True)
 
     Returns
     -------
-    number: np.ndarray
-        number of topographies having qx in their bandwidth
+    lower : np.ndarray
+        Sorted lower bandwidth bounds of the individual topographies.
+    upper : np.ndarray
+        Sorted upper bandwidth bounds of the individual topographies.
     """
-    qx = np.abs(qx)
-    factor = np.zeros_like(qx)
+    lower = []
+    upper = []
     for t in self:
         t = t.to_unit(unit)
 
@@ -58,14 +62,55 @@ def _bandwidth_count_from_profile(self, qx, unit, reliable=True):
         if short_cutoff is not None:
             qxmax = min(2 * np.pi / short_cutoff, qxmax)
 
-        factor += np.logical_and(
-            qx >= 2 * np.pi / t.physical_sizes[0],
-            qx <= qxmax,
-        )
+        qxmin = 2 * np.pi / t.physical_sizes[0]
+        if qxmax >= qxmin:
+            lower += [qxmin]
+            upper += [qxmax]
+        # else: empty bandwidth interval (reliability cutoff longer than the
+        # scan); such a topography contains no wavevector and must not enter
+        # the counting arrays, where the subtraction in `_count_in_intervals`
+        # would tally it as -1 between its inverted bounds
+    return np.sort(lower), np.sort(upper)
+
+
+def _count_in_intervals(intervals, nb_topographies, qx):
+    r"""
+    Return the number of bandwidth intervals (from
+    `_bandwidth_intervals_from_profile`) that contain each wavevector in
+    `qx`.
+    """
+    lower, upper = intervals
+    qx = np.abs(qx)
+    # Number of intervals with lower <= qx minus number of intervals with
+    # upper < qx
+    count = np.searchsorted(lower, qx, side='right') - np.searchsorted(upper, qx, side='left')
     # All topographies have qx == 0 wavevector
-    factor = np.where(qx == 0, len(self), factor)
-    # TODO: This will be very slow on Topobank
-    return factor
+    return np.where(qx == 0, nb_topographies, count)
+
+
+def _bandwidth_count_from_profile(self, qx, unit, reliable=True):
+    r"""
+    Return number of topographies that include qx in their bandwidth.
+
+    Parameters:
+    -----------
+    self : SurfaceContainer
+        Collection of Height containers
+    qx: np.ndarrau of floats
+        wavevector
+    unit : str
+        Unit of lengths in which the wavevector is defined.
+    reliable : bool, optional
+        Only incorporate data deemed reliable. (Default: True)
+
+    Returns
+    -------
+    number: np.ndarray
+        number of topographies having qx in their bandwidth
+    """
+    return _count_in_intervals(
+        _bandwidth_intervals_from_profile(self, unit, reliable), len(self), qx
+    )
 
 
 def integrate_psd_from_profile(self, factor, unit, window=None, reliable=True):
@@ -84,7 +129,7 @@ def integrate_psd_from_profile(self, factor, unit, window=None, reliable=True):
 
     .. math::
 
-        \frac{1}{2 \pi} \int_0^\infty dq_x factor(q_x) C^{1D}(q_x)
+        \frac{1}{2 \pi} \int_{-\infty}^\infty dq_x factor(q_x) C^{1D}(q_x)
 
     Discrete
 
@@ -118,13 +163,16 @@ def integrate_psd_from_profile(self, factor, unit, window=None, reliable=True):
     """
     integ = 0
 
-    # TODO: faster: precompute _bandwidth_count_from_profile (defining a piecewise constant function, )
-    # a table with bins and number of topographies inside the bins
-    # this will make us only loop 2N times through the topographies instead of N^2
-    # with N the total number of topographies
+    # Precompute the bandwidth intervals in a single pass over the
+    # container; the `average` callback below is invoked once per
+    # topography, and looping over the container inside it would read
+    # every data file N times (lazy containers construct topographies
+    # from the file on each element access).
+    intervals = _bandwidth_intervals_from_profile(self, unit, reliable)
+    nb_topographies = len(self)
 
     def average(qx):
-        count = _bandwidth_count_from_profile(self, qx, unit, reliable)
+        count = _count_in_intervals(intervals, nb_topographies, qx)
         return np.where(count > 0, factor(qx) / count, 0)
 
     for t in self: