Improve 2D compression performance by taking advantage of the int_8 msb info. Also bump version

CMakeLists.txt

@@ -2,7 +2,7 @@
 
 cmake_minimum_required(VERSION 3.14)
 
-project(SPERR VERSION 0.8.4 LANGUAGES CXX DESCRIPTION "Lossy Scientific Compression with SPERR")
+project(SPERR VERSION 0.8.5 LANGUAGES CXX DESCRIPTION "Lossy Scientific Compression with SPERR")
 
 if(NOT CMAKE_CXX_STANDARD)
     set(CMAKE_CXX_STANDARD "20" CACHE STRING "Choose the C++ Standard to use." FORCE)
@@ -109,7 +109,7 @@ if( BUILD_CLI_UTILITIES )
   set( CLI11_SINGLE_FILE OFF CACHE INTERNAL "Don't use single file CLI11")
   FetchContent_Declare( cli11
     GIT_REPOSITORY https://github.com/CLIUtils/CLI11
-    GIT_TAG        6c7b07a878ad834957b98d0f9ce1dbe0cb204fc9 # v2.4.2
+    GIT_TAG        37bb6edc5317e99af72ef48405e65d9ca5218861 # v2.6.2
   )
   FetchContent_MakeAvailable(cli11)
 

include/SPECK2D_INT.h

@@ -44,6 +44,7 @@ class SPECK2D_INT : public SPECK_INT<T> {
   virtual void m_process_S(size_t idx1, size_t idx2, size_t& counter, bool need_decide) = 0;
   virtual void m_process_P(size_t idx, size_t& counter, bool need_decide) = 0;
   virtual void m_process_I(bool need_decide) = 0;
+  virtual void m_additional_initialization() {};
 
   auto m_partition_S(Set2D) const -> std::array<Set2D, 4>;
   auto m_partition_I() -> std::array<Set2D, 3>;

include/SPECK2D_INT_ENC.h

@@ -29,9 +29,17 @@ class SPECK2D_INT_ENC final : public SPECK2D_INT<T> {
   void m_process_S(size_t idx1, size_t idx2, size_t& counter, bool need_decide) final;
   void m_process_P(size_t idx, size_t& counter, bool need_decide) final;
   void m_process_I(bool need_decide) final;
+  void m_additional_initialization() final;
+  void m_bitplane_init() final;
+  void m_refinement_extra() final;
 
   auto m_decide_S_significance(const Set2D&) const -> bool;
   auto m_decide_I_significance() const -> bool;
+
+  // `m_msb_buf` stores the MSB bit position of each coefficient, in the same order as
+  // m_coeff_buf. Significance tests compare entries against `m_msb_threshold`.
+  std::vector<int8_t> m_msb_buf;
+  int8_t m_msb_threshold = -1;
 };
 
 };  // namespace sperr

include/SPECK3D_INT_ENC.h

@@ -47,8 +47,6 @@ class SPECK3D_INT_ENC final : public SPECK3D_INT<T> {
   // m_bitplane_init(). Significance tests compare m_morton_buf entries against this value.
   int8_t m_morton_threshold = -1;
   void m_deposit_set(Set3D);
-  // Returns the bit position of the most significant bit (0-based), or -1 for zero.
-  auto m_msb_position(uint_type v) const -> int8_t;
 };
 
 };  // namespace sperr

include/sperr_helper.h

@@ -166,6 +166,10 @@ auto calc_stats(const T* arr1, const T* arr2, size_t arr_len, size_t omp_nthread
 template <typename T>
 auto kahan_summation(const T*, size_t) -> T;
 
+// Returns the bit position of the most significant bit (0-based), or -1 for zero.
+template <typename T>
+auto msb_position(T v) -> int8_t;
+
 // Given a whole volume size and a desired chunk size, this helper function
 // returns a list of chunks specified by 6 integers:
 // chunk[0], [2], [4]: starting index of this chunk in X, Y, and Z;

src/SPECK2D_INT.cpp

@@ -212,6 +212,9 @@ void sperr::SPECK2D_INT<T>::m_initialize_lists()
   m_I.length_x = m_dims[0];
   m_I.length_y = m_dims[1];
   m_I.part_level = num_of_xforms;
+
+  // Encoder and decoder might have different additional tasks.
+  m_additional_initialization();
 }
 
 template class sperr::SPECK2D_INT<uint8_t>;

src/SPECK2D_INT_ENC.cpp

@@ -2,6 +2,7 @@
 
 #include <algorithm>
 #include <cassert>
+#include <numeric>
 
 template <typename T>
 void sperr::SPECK2D_INT_ENC<T>::m_process_S(size_t idx1,
@@ -31,13 +32,12 @@ void sperr::SPECK2D_INT_ENC<T>::m_process_P(size_t idx, size_t& counter, bool ne
   bool is_sig = true;
 
   if (need_decide) {
-    is_sig = (m_coeff_buf[idx] >= m_threshold);
+    is_sig = (m_msb_buf[idx] >= m_msb_threshold);
     m_bit_buffer.wbit(is_sig);
   }
 
   if (is_sig) {
     counter++;
-    m_coeff_buf[idx] -= m_threshold;
     m_bit_buffer.wbit(m_sign_array.rbit(idx));
     m_LSP_new.push_back(idx);
     m_LIP_mask.wfalse(idx);
@@ -65,8 +65,9 @@ auto sperr::SPECK2D_INT_ENC<T>::m_decide_S_significance(const Set2D& set) const
   assert(!set.is_empty());
 
   for (auto y = set.start_y; y < (set.start_y + set.length_y); y++) {
-    auto first = m_coeff_buf.data() + y * m_dims[0] + set.start_x;
-    if (std::any_of(first, first + set.length_x, [th = m_threshold](auto v) { return v >= th; }))
+    auto first = m_msb_buf.cbegin() + y * m_dims[0] + set.start_x;
+    if (std::any_of(first, first + set.length_x,
+                    [thld = m_msb_threshold](auto v) { return v >= thld; }))
       return true;
   }
   return false;
@@ -79,22 +80,46 @@ auto sperr::SPECK2D_INT_ENC<T>::m_decide_I_significance() const -> bool
   // It's stored in a contiguous chunk of memory till the buffer end.
   //
   assert(m_I.length_x == m_dims[0]);
-  auto first = m_coeff_buf.data() + size_t{m_I.start_y} * size_t{m_I.length_x};
-  auto len = m_coeff_buf.size() - size_t{m_I.start_y} * size_t{m_I.length_x};
-  if (std::any_of(first, first + len, [thld = m_threshold](auto v) { return v >= thld; }))
+  auto first = m_msb_buf.cbegin() + size_t{m_I.start_y} * size_t{m_I.length_x};
+  auto len = m_msb_buf.size() - size_t{m_I.start_y} * size_t{m_I.length_x};
+  if (std::any_of(first, first + len, [thld = m_msb_threshold](auto v) { return v >= thld; }))
     return true;
 
   // Second, test the rectangle that's directly to the right of the missing top-left corner.
   //
   len = m_dims[0] - m_I.start_x;
   for (auto y = 0u; y < m_I.start_y; y++) {
-    first = m_coeff_buf.data() + y * m_dims[0] + m_I.start_x;
-    if (std::any_of(first, first + len, [thld = m_threshold](auto v) { return v >= thld; }))
+    first = m_msb_buf.cbegin() + y * m_dims[0] + m_I.start_x;
+    if (std::any_of(first, first + len, [thld = m_msb_threshold](auto v) { return v >= thld; }))
       return true;
   }
   return false;
 }
 
+template <typename T>
+void sperr::SPECK2D_INT_ENC<T>::m_additional_initialization()
+{
+  const auto len = m_dims[0] * m_dims[1] * m_dims[2];
+  m_msb_buf.resize(len);
+  std::transform(m_coeff_buf.cbegin(), m_coeff_buf.cend(), m_msb_buf.begin(),
+                 [](auto v) { return sperr::msb_position(v); });
+}
+
+template <typename T>
+void sperr::SPECK2D_INT_ENC<T>::m_bitplane_init()
+{
+  m_msb_threshold = sperr::msb_position(m_threshold);
+}
+
+template <typename T>
+void sperr::SPECK2D_INT_ENC<T>::m_refinement_extra()
+{
+  for (auto idx : m_LSP_new) {
+    assert(m_coeff_buf[idx] >= m_threshold);
+    m_coeff_buf[idx] -= m_threshold;
+  }
+}
+
 template class sperr::SPECK2D_INT_ENC<uint8_t>;
 template class sperr::SPECK2D_INT_ENC<uint16_t>;
 template class sperr::SPECK2D_INT_ENC<uint32_t>;

src/SPECK3D_INT_ENC.cpp

@@ -5,10 +5,6 @@
 #include <cstring>  // std::memcpy()
 #include <numeric>
 
-#if __cplusplus >= 202002L
-#include <bit>
-#endif
-
 template <typename T>
 void sperr::SPECK3D_INT_ENC<T>::m_deposit_set(Set3D set)
 {
@@ -17,7 +13,7 @@ void sperr::SPECK3D_INT_ENC<T>::m_deposit_set(Set3D set)
       return;
     case 1: {
       auto id = set.start_z * m_dims[0] * m_dims[1] + set.start_y * m_dims[0] + set.start_x;
-      m_morton_buf[set.morton_idx] = m_msb_position(m_coeff_buf[id]);
+      m_morton_buf[set.morton_idx] = sperr::msb_position(m_coeff_buf[id]);
       return;
     }
     case 2: {
@@ -26,7 +22,7 @@ void sperr::SPECK3D_INT_ENC<T>::m_deposit_set(Set3D set)
       // Deposit the 1st element.
       auto id = set.start_z * m_dims[0] * m_dims[1] + set.start_y * m_dims[0] + set.start_x;
       auto morton_id = set.morton_idx;
-      m_morton_buf[morton_id] = m_msb_position(m_coeff_buf[id]);
+      m_morton_buf[morton_id] = sperr::msb_position(m_coeff_buf[id]);
 
       // Deposit the 2nd element.
       if (set.length_x == 2)
@@ -35,7 +31,7 @@ void sperr::SPECK3D_INT_ENC<T>::m_deposit_set(Set3D set)
         id += m_dims[0];
       else
         id += m_dims[0] * m_dims[1];
-      m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id]);
+      m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id]);
 
       return;
     }
@@ -45,51 +41,51 @@ void sperr::SPECK3D_INT_ENC<T>::m_deposit_set(Set3D set)
 
       if (set.length_x == 2 && set.length_y == 2) {
         // Element (0, 0, 0)
-        m_morton_buf[morton_id] = m_msb_position(m_coeff_buf[id]);
+        m_morton_buf[morton_id] = sperr::msb_position(m_coeff_buf[id]);
 
         // Element (1, 0, 0)
-        m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id + 1]);
+        m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id + 1]);
 
         // Element (0, 1, 0)
         auto id2 = id + m_dims[0];
-        m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id2]);
+        m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id2]);
 
         // Element (1, 1, 0)
-        m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[++id2]);
+        m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[++id2]);
 
         return;
       }
       else if (set.length_x == 2 && set.length_z == 2) {
         // Element (0, 0, 0)
-        m_morton_buf[morton_id] = m_msb_position(m_coeff_buf[id]);
+        m_morton_buf[morton_id] = sperr::msb_position(m_coeff_buf[id]);
 
         // Element (1, 0, 0)
-        m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id + 1]);
+        m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id + 1]);
 
         // Element (0, 0, 1)
         auto id2 = id + m_dims[0] * m_dims[1];
-        m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id2]);
+        m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id2]);
 
         // Element (1, 0, 1)
-        m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[++id2]);
+        m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[++id2]);
 
         return;
       }
       else if (set.length_y == 2 && set.length_z == 2) {
         // Element (0, 0, 0)
-        m_morton_buf[morton_id] = m_msb_position(m_coeff_buf[id]);
+        m_morton_buf[morton_id] = sperr::msb_position(m_coeff_buf[id]);
 
         // Element (0, 1, 0)
         auto id2 = id + m_dims[0];
-        m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id2]);
+        m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id2]);
 
         // Element (0, 0, 1)
         id2 = id + m_dims[0] * m_dims[1];
-        m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id2]);
+        m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id2]);
 
         // Element (0, 1, 1)
         id2 += m_dims[0];
-        m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id2]);
+        m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id2]);
 
         return;
       }
@@ -101,31 +97,31 @@ void sperr::SPECK3D_INT_ENC<T>::m_deposit_set(Set3D set)
         // Element (0, 0, 0)
         const auto id = set.start_z * m_dims[0] * m_dims[1] + set.start_y * m_dims[0] + set.start_x;
         auto morton_id = set.morton_idx;
-        m_morton_buf[morton_id] = m_msb_position(m_coeff_buf[id]);
+        m_morton_buf[morton_id] = sperr::msb_position(m_coeff_buf[id]);
 
         // Element (1, 0, 0)
-        m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id + 1]);
+        m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id + 1]);
 
         // Element (0, 1, 0)
         auto id2 = id + m_dims[0];
-        m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id2]);
+        m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id2]);
 
         // Element (1, 1, 0)
-        m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[++id2]);
+        m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[++id2]);
 
         // Element (0, 0, 1)
         id2 = id + m_dims[0] * m_dims[1];
-        m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id2]);
+        m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id2]);
 
         // Element (1, 0, 1)
-        m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[++id2]);
+        m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[++id2]);
 
         // Element (0, 1, 1)
         id2 = id + m_dims[0] * (m_dims[1] + 1);
-        m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[id2]);
+        m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[id2]);
 
         // Element (1, 1, 1)
-        m_morton_buf[++morton_id] = m_msb_position(m_coeff_buf[++id2]);
+        m_morton_buf[++morton_id] = sperr::msb_position(m_coeff_buf[++id2]);
 
         return;
       }
@@ -189,7 +185,7 @@ void sperr::SPECK3D_INT_ENC<T>::m_process_P(size_t idx, size_t morton, size_t& c
   bool is_sig = true;
 
   if (output) {
-    assert(m_msb_position(m_coeff_buf[idx]) == m_morton_buf[morton]);
+    assert(sperr::msb_position(m_coeff_buf[idx]) == m_morton_buf[morton]);
     is_sig = (m_morton_buf[morton] >= m_morton_threshold);
     m_bit_buffer.wbit(is_sig);
   }
@@ -215,25 +211,10 @@ void sperr::SPECK3D_INT_ENC<T>::m_process_P_lite(size_t idx)
   }
 }
 
-template <typename T>
-auto sperr::SPECK3D_INT_ENC<T>::m_msb_position(uint_type v) const -> int8_t
-{
-#if __cplusplus >= 202002L
-  return static_cast<int8_t>(sizeof(uint_type) * 8 - 1 - std::countl_zero(v));
-#else
-  int8_t pos = -1;
-  while (v) {
-    v >>= 1;
-    pos++;
-  }
-  return pos;
-#endif
-}
-
 template <typename T>
 void sperr::SPECK3D_INT_ENC<T>::m_bitplane_init()
 {
-  m_morton_threshold = m_msb_position(m_threshold);
+  m_morton_threshold = sperr::msb_position(m_threshold);
 }
 
 template <typename T>

src/sperr_helper.cpp

@@ -7,6 +7,10 @@
 #include <cstring>
 #include <numeric>
 
+#if __cplusplus >= 202002L
+#include <bit>
+#endif
+
 #ifdef __AVX2__
 #include <immintrin.h>
 #endif
@@ -641,3 +645,22 @@ auto sperr::calc_mean_var(const T* arr, size_t len, size_t omp_nthreads) -> std:
 }
 template auto sperr::calc_mean_var(const float*, size_t, size_t) -> std::array<float, 2>;
 template auto sperr::calc_mean_var(const double*, size_t, size_t) -> std::array<double, 2>;
+
+template <typename T>
+auto sperr::msb_position(T v) -> int8_t
+{
+#if __cplusplus >= 202002L
+  return static_cast<int8_t>(sizeof(T) * 8 - 1 - std::countl_zero(v));
+#else
+  int8_t pos = -1;
+  while (v) {
+    v >>= 1;
+    pos++;
+  }
+  return pos;
+#endif
+}
+template auto sperr::msb_position(uint8_t) -> int8_t;
+template auto sperr::msb_position(uint16_t) -> int8_t;
+template auto sperr::msb_position(uint32_t) -> int8_t;
+template auto sperr::msb_position(uint64_t) -> int8_t;