Merge commit 'f77d4f5db9cbab4488a85f078dc7b5aaa08d0251' into update-h5z-sperr

Author: t20100

Diff for: src/H5Z-SPERR/.gitignore

@@ -1,3 +1,6 @@
+# vim temp files
+*.sw?
+
 # Prerequisites
 *.d
 

Diff for: src/H5Z-SPERR/CMakeLists.txt

@@ -1,15 +1,17 @@
 cmake_minimum_required(VERSION 3.14)
 
-project(H5Z-MD5 VERSION 0.1.3 LANGUAGES C DESCRIPTION "HDF5 plugin for SPERR compression")
+project(H5Z-SPERR VERSION 0.2.3 LANGUAGES C CXX DESCRIPTION "HDF5 plugin for SPERR compression")
 
 set(CMAKE_C_STANDARD 11)
+set(CMAKE_CXX_STANDARD 14)
 if(NOT CMAKE_BUILD_TYPE)
     set(CMAKE_BUILD_TYPE "Release" CACHE STRING "Choose the type of build." FORCE)
     set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "Debug" "Release" "RelWithDebInfo")
 endif()
 
 option( BUILD_SHARED_LIBS "Build shared libraries" ON )
 option( BUILD_CLI_UTILITIES "Build a set of command line utilities" ON )
+option( BUILD_UNIT_TESTS "Build unit tests using GoogleTest" OFF )
 option( H5ZPLUGIN_PREFER_RPATH "Set RPATH; this can fight with package managers 
                                 so turn off when building for them" ON )
 mark_as_advanced(FORCE H5ZPLUGIN_PREFER_RPATH)
@@ -42,11 +44,39 @@ if( BUILD_CLI_UTILITIES )
   add_subdirectory( utilities ${CMAKE_BINARY_DIR}/bin )
 endif()
 
+#
+# Build unit tests
+#
+if( BUILD_UNIT_TESTS )
+  # Control internal options of GoogleTest
+  #
+  set( INSTALL_GTEST OFF CACHE INTERNAL "Not install GoogleTest")
+  set( BUILD_GMOCK ON CACHE INTERNAL "Build gmock")
+
+  # Let's use the new mechanism to incorporate GoogleTest
+  #
+  include(FetchContent)
+  if(${CMAKE_VERSION} VERSION_GREATER_EQUAL "3.24")
+    FetchContent_Declare( googletest
+      URL https://github.com/google/googletest/archive/refs/heads/main.zip
+      DOWNLOAD_EXTRACT_TIMESTAMP NEW )
+  endif()
+
+  # Prevent overriding the parent project's compiler/linker settings on Windows
+  #
+  set(gtest_force_shared_crt ON CACHE BOOL "" FORCE)
+  FetchContent_MakeAvailable(googletest)
+
+  enable_testing() # calling this function before adding subdirectory to enable 
+                   # invoking ctest from the top-level build directory.
+  add_subdirectory( test_scripts )
+endif()
+
 #
 # Start installation using GNU installation rules
 #
 include( GNUInstallDirs )
-install( TARGETS h5z-sperr LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} )
+install( TARGETS h5z-sperr h5z-clamp LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} )
 
 if( BUILD_CLI_UTILITIES )
   install( TARGETS generate_cd_values decode_cd_values

Diff for: src/H5Z-SPERR/README.md

@@ -24,6 +24,66 @@ export HDF5_PLUGIN_PATH=/path/to/install/this/plugin
 ```
 The user program does not need to link to this plugin or SPERR; it only needs to specify the plugin ID of `32028`.
 
+<!--
+## Use in NetCDF-4 APIs
+`H5Z-SPERR` also facilitates the application of SPERR compression on 
+[NetCDF-4 files](https://docs.unidata.ucar.edu/netcdf/NUG/md_filters.html#filters_enable);
+one simply needs to define the filter on a variable:
+```C
+nc_def_var_filter(ncid, varid, 32028, 1, &cd_values);
+```
+See a complete example [here](https://github.com/NCAR/H5Z-SPERR/blob/main/example/simple_xy_nc4_wr.c).
+-->
+
+## Use in Python
+`H5Z-SPERR` version `0.1.3` is supported by the Python package [hdf5plugin](https://github.com/silx-kit/hdf5plugin)
+since version `5.0.0`.
+One can install the package by issuing 
+```bash
+pip install hdf5plugin [--user]           # using pip
+conda install -c conda-forge hdf5plugin   # using conda
+```
+and use it by importing these two packages:
+```python
+import h5py         # provide general HDF5 support
+import hdf5plugin   # provide HDF5 plugin support
+```
+
+## Handling of Missing Values
+Simulation models sometimes use a special value (i.e., missing value) to indicate that there's no meaningful value at that specific location.
+For example, in an ocean simulation, all the land area is marked by missing values.
+The most often seen missing values are either `NaN`, or an extremely large value, such as `1e35` or `-9.9e35`.
+When these missing values participate in compression, they easily introduce numeric error and result in data corruption.
+
+`H5Z-SPERR` can handle common missing values with a little help from the user.
+Specifically, a user can indicate that there's potentially missing values, 
+and `H5Z-SPERR` will use a compact bitmask to keep track of where exactly those missing values are.
+`H5Z-SPERR` then replaces them with a value that is friendly to SPERR compression before passing the field to the SPERR compressor.
+
+During decompression, `H5Z-SPERR` fills the original missing value at locations indicated by the compact bitmask.
+Specifically, if the original missing values are `NaN`, then `NaN` will be filled. If the orignal missing values
+are values with a magnitude larger than `1e35`, then the *first occurance* of such values will be used to fill 
+in all missing value locations.
+
+Users use an integer to indicate the potential existance of missing values:
+- Mode `0`: no missing values;
+- Mode `1`: there are potential `NaN`s;
+- Mode `2`: there are potential values with a magnitude larger than `1e35`.
+
+`H5Z-SPERR` behaves accordingly: (`1e35` denotes the *first occurance* of such values)
+| Mode      | Actual Input Data    |  Filter Behavior |
+|-----------|----------------------|------------------|
+| 0         | No `NaN`, no `1e35`  | :heavy_check_mark: Normal SPERR compression |
+| 0         | Has `NaN` or `1e35`  | :x: Likely numeric error  |
+| 1         | No `NaN`, no `1e35`  | :heavy_check_mark: Normal SPERR compression  |
+| 1         | Has `NaN`, no `1e35` | :heavy_check_mark: Normal SPERR compression; `NaN` is restored at its exact locations  |
+| 1         | Regardless of `NaN`, has `1e35` |  :x: Likely numeric error  |
+| 2         | No `NaN`, no `1e35`  | :heavy_check_mark: Normal SPERR compression  |
+| 2         | No `NaN`, has `1e35` | :heavy_check_mark: Normal SPERR compression; `1e35` is restored at its exact locations  |
+| 2         | Has `NaN`, regardless of `1e35` | :x: Likely numeric error |
+
+**Final note:** if a variable is indicated to have missing values, but it actually does not, then there's no bitmasks involved thus no storage overhead! 
+
 ##  Find `cd_values[]`
 To apply SPERR compression using the HDF5 plugin, one needs to specify 1) what compression mode and 2)
 what compression quality to use. Supported compression modes and qualities are summarized below:
@@ -42,7 +102,8 @@ and the `Z` rank to be varying the slowest, before the data is passed to the com
 
 The HDF5 libraries takes in these compression parameters as one or more 32-bit `unsigned int` values,
 which are named `cd_values[]` in most HDF5 routines.
-In the case of `H5Z-SPERR`, there is exactly one `unsigned int` used to carry this information.
+In the case of `H5Z-SPERR`, there is exactly one `unsigned int` used to carry compression-related information, and 
+possibly one more `unsigned int` to indicate the potential existance of missing values.
 
 ### Find  `cd_values[]` Using the Programming Interface
 Using the HDF5 programming interface, `cd_values[]` carrying the compression parameters are passed
@@ -75,11 +136,17 @@ Please use this value as a single 32-bit unsigned integer in your applications.
 Note: an integer produced by `generate_cd_values` can be decoded by another command line tool, `decode_cd_values`,
 to show the coded compression parameters.
 
-## Use in NetCDF-4 APIs
-`H5Z-SPERR` also facilitates the application of SPERR compression on 
-[NetCDF-4 files](https://docs.unidata.ucar.edu/netcdf/NUG/md_filters.html#filters_enable);
-one simply needs to define the filter on a variable:
-```C
-nc_def_var_filter(ncid, varid, 32028, 1, &cd_values);
+### Examples
+Assume using the `nccopy` tool:
+```Bash
+# Compress variable VAR0, using fixed-rate compression, bitrate = 3.3, no special handling of missing values.
+nccopy -F "VAR0, 268651725u" <input_file> <output_file>
+nccopy -F "VAR0, 268651725u, 0" <input_file> <output_file>
+
+# Compress variable VAR1, using fixed-rate compression, bitrate = 3.3. VAR1 might have NaNs!
+nccopy -F "VAR1, 268651725u, 1" <input_file> <output_file>
+
+# Compress variable VAR2, using fixed-rate compression, bitrate = 3.3. VAR2 might have values such as 1e35!
+nccopy -F "VAR2, 268651725u, 2" <input_file> <output_file> 
 ```
-See a complete example [here](https://github.com/NCAR/H5Z-SPERR/blob/main/example/simple_xy_nc4_wr.c).
+

Diff for: src/H5Z-SPERR/include/compactor.h

@@ -0,0 +1,74 @@
+/* This is a set of functions that compact a bitmask.
+ * In the intended use case, a bitmask is produced by masking all
+ * "missing values" or "fill values" in a model output with zero's,
+ * whereas the locations with valid data points are marked with one's.
+ * However, this compactor is likely to be effective with any bit patterns
+ * that have lots of consecutive 0's or 1's.
+ *
+ * The bitmask compactor works in the following way:
+ * 1. Assume that we use 32-bit ints; the compactor encodes 32 bits at a time.
+ * 2. Every incoming int is encoded in one of three ways:
+ *    2.1. For an int with all 0's, use a single 0 bit.
+ *    2.2. For an int with all 1's, use two bits: 10.
+ *    2.3. For all other ints, use 34 bits: two bits 11 then followed by the
+ *         verbose presentation of the 32-bit int.
+ * 3. Obviously, it's the most economical to use a single 0 bit to present
+ *    the most frequent pattern (all 0's or all 1's). The encoder thus does
+ *    a test at the beginning and records the test result.
+ */
+
+#ifndef COMPACTOR_H
+#define COMPACTOR_H
+
+#include <stdlib.h>
+#include <stdint.h>
+#include <assert.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/* Return the compaction strategy to use:
+ * 0: compact with all 0's being the most frequent
+ * 1: compact with all 1's being the most frequent
+ * Note: `bytes` has to be a multiple of 8.
+ */
+int compactor_strategy(const void* buf, size_t bytes);
+
+/* Return the size in bytes of the resulting compacted bitstream, given an input buf.
+ * Note: `buf_bytes` has to be a multiple of 8.
+ */
+size_t compactor_comp_size(const void* buf, size_t buf_bytes);
+
+/* Return the number of useful bytes in a compacted bitstream.
+ * This value is the same as the output of `compactor_comp_size()` during encoding.
+ */
+size_t compactor_useful_bytes(const void* comp_buf);
+
+/* Return the useful size of the output bitstream,
+ * which is the same as the output of `compactor_comp_size()`.
+ * Note 1: the input bitmask length (in bytes) has to be a multiple of 8.
+ *         This requirement is inheritated from the bitstream implementation.
+ * Note 2: the output buffer length should be 1) a multiple of 8, and
+ *         2) no less than the size returned by `compactor_comp_size()`.
+ */
+size_t compactor_encode(const void* bitmask,
+                        size_t bitmask_bytes,
+                        void* compact_bitstream,
+                        size_t compact_bitstream_bytes);
+
+/* Return the number of useful bytes in the decoded bitmask.
+ * Note: The number of useful bytes might be bigger than the number of bytes being
+ *       encoded, because of the word size that the compactor operates on.
+ * Note: `compact_bitstream_bytes` should be a multiple of 8 that is no less than
+ *       the size returned by `compactor_encode()`.
+ */
+size_t compactor_decode(const void* compact_bitstream,
+                        size_t compact_bitstream_bytes,
+                        void* decoded_bitmask);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif

Diff for: src/H5Z-SPERR/include/h5zsperr_helper.h

@@ -0,0 +1,61 @@
+/*
+ * This file contains a few helper functions for the H5Z-SPERR filter.
+ */
+
+#ifndef H5ZSPERR_HELPER_H
+#define H5ZSPERR_HELPER_H
+
+#include <stdlib.h>
+
+#define LARGE_MAGNITUDE_F 1e35f
+#define LARGE_MAGNITUDE_D 1e35
+#define H5ZSPERR_COMPATIBILITY 1
+
+#ifdef __cplusplus
+namespace C_API {
+extern "C" {
+#endif
+
+/*
+ * Pack and unpack additional information about the input data into an integer.
+ * It returns the encoded unsigned int, which shouldn't be zero.
+ * The packing function is called by `set_local()` to prepare information
+ * for `H5Z_filter_sperr()`, which calls the unpack function to extract such info.
+ */
+unsigned int h5zsperr_pack_extra_info(int rank, int is_float, int missing_val_mode, int magic_num);
+void h5zsperr_unpack_extra_info(unsigned int meta,
+                                int* rank,
+                                int* is_float,
+                                int* missing_val_mode,
+                                int* magic_num);
+
+/*
+ * Check if an input array really has missing values.
+ */
+int h5zsperr_has_nan(const void* buf, size_t nelem, int is_float);
+int h5zsperr_has_large_mag(const void* buf, size_t nelem, int is_float);
+
+/*
+ * Produce a compact bitmask.
+ * `mask_buf` is already allocated with length `mask_bytes`.
+ * Returns 0 upon success.
+ */
+int h5zsperr_make_mask_nan(const void* data_buf, size_t nelem, int is_float,
+                           void* mask_buf, size_t mask_bytes, size_t* useful_bytes);
+int h5zsperr_make_mask_large_mag(const void* data_buf, size_t nelem, int is_float,
+                                 void* mask_buf, size_t mask_bytes, size_t* useful_bytes);
+
+/*
+ * Replace every missing value in the `data_buf` with the mean of the field.
+ */
+float h5zsperr_treat_nan_f32(float* data_buf, size_t nelem);
+double h5zsperr_treat_nan_f64(double* data_buf, size_t nelem);
+float h5zsperr_treat_large_mag_f32(float* data_buf, size_t nelem);
+double h5zsperr_treat_large_mag_f64(double* data_buf, size_t nelem);
+
+#ifdef __cplusplus
+} /* end of extern "C" */
+} /* end of namespace C_API */
+#endif
+
+#endif

Diff for: src/H5Z-SPERR/include/icecream.h

@@ -0,0 +1,64 @@
+/*
+ * This is a mimic of the Bitstream class in SPERR:
+ * https://github.com/NCAR/SPERR/blob/main/include/Bitstream.h
+ *
+ * The most significant difference is that bitstream here doesn't manage
+ * any memory; it reads a bit sequence from a user-provided memory buffer,
+ * or writes a bit sequence to a user-provided memory buffer.
+ * 
+ * The "object" is named `icecream` and all functions operating on it
+ * are named with a prefix `icecream`.
+ */
+
+#ifndef ICECREAM_H
+#define ICECREAM_H
+
+#include <stdlib.h>
+#include <stdint.h>
+#include <assert.h>
+
+#ifndef NDEBUG
+#include <stdio.h>
+#endif
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef struct {
+  uint64_t* begin;  /* begin of the stream */
+  uint64_t* ptr;    /* pointer to the next word to be read/written */
+  uint64_t  buffer; /* incoming/outgoing bits */
+  int       bits;   /* number of buffered bits (0 <= bits < 64) */
+} icecream;
+
+/*
+ * Specify a bitstream to use memory provided by users.
+ * NOTE: the memory length (in bytes) have to be a multiplier of 8,
+ * because the icecream class writes/reads in 64-bit integers.
+ */
+void icecream_use_mem(icecream* s, void* mem, size_t bytes);
+
+/* Position the bitstream for reading or writing at the beginning. */
+void icecream_rewind(icecream* s);
+
+/* Read a bit. Please don't read beyond the end of the stream. */
+int icecream_rbit(icecream* s);
+
+/* Write a bit (0 or 1). Please don't write beyond the end of the stream. */
+void icecream_wbit(icecream* s, int bit);
+
+/* Return the bit offset to the next bit to be read. */
+size_t icecream_rtell(icecream* s);
+
+/* Return the bit offset to the next bit to be written. */
+size_t icecream_wtell(icecream* s);
+
+/* Write any remaining buffered bits and align stream on next word boundary. */
+void icecream_flush(icecream* s);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif

Diff for: src/H5Z-SPERR/src/CMakeLists.txt

@@ -1,5 +1,18 @@
-add_library( h5z-sperr h5z-sperr.c )
+#
+# The main H5Z-SPERR filter
+#
+add_library( h5z-sperr h5z-sperr.c
+                       h5zsperr_helper.cpp
+                       icecream.c
+                       compactor.c)
 target_include_directories( h5z-sperr PUBLIC ${HDF5_INCLUDE_DIR} 
                                       PUBLIC ${SPERR_INCLUDE_DIRS}
                                       PUBLIC ${CMAKE_SOURCE_DIR}/include )
 target_link_libraries( h5z-sperr PUBLIC ${HDF5_LIBRARIES} PUBLIC PkgConfig::SPERR )
+
+#
+# Experimental clamping filter
+#
+add_library( h5z-clamp h5z-clamp.c )
+target_include_directories( h5z-clamp PUBLIC ${HDF5_INCLUDE_DIR} )
+target_link_libraries( h5z-clamp PUBLIC ${HDF5_LIBRARIES} )