diff --git a/.github/workflows/build-docker-image.yml b/.github/workflows/build-docker-image.yml
index 22a0533f..8b80e656 100644
--- a/.github/workflows/build-docker-image.yml
+++ b/.github/workflows/build-docker-image.yml
@@ -6,7 +6,7 @@ name: Build Docker image
 on:
   create:
     branches-ignore:
-      - "*"
+      - "**"
     tags:
       - 'v*.*.*'
 
diff --git a/.github/workflows/unit-tests-macos.yml b/.github/workflows/unit-tests-macos.yml
index 449db945..2db8d0f7 100644
--- a/.github/workflows/unit-tests-macos.yml
+++ b/.github/workflows/unit-tests-macos.yml
@@ -47,9 +47,18 @@ jobs:
     steps:
       - uses: actions/checkout@v3
 
-      - run: |
+      - name: Detect number available CPUs
+        run: |
+          ncpus=$(python -c 'import multiprocessing as mp; print(mp.cpu_count())')
+
+          echo "CMAKE_BUILD_PARALLEL_LEVEL=$ncpus" >> $GITHUB_ENV
+          echo "CTEST_PARALLEL_LEVEL=$ncpus" >> $GITHUB_ENV
+
+      - name: Generate requirements.txt for pip
+        run: |
           echo 'conan==${{ matrix.conan-ver }}' > requirements.txt
           echo 'cmake==${{ matrix.cmake-ver }}' >> requirements.txt
+          echo 'scipy' >> requirements.txt
 
       - uses: actions/setup-python@v4
         with:
@@ -104,19 +113,34 @@ jobs:
                 -B "${{ github.workspace }}/build"
 
       - name: Build project
-        run: cmake --build ${{ github.workspace }}/build -j $(nproc)
+        run: cmake --build ${{ github.workspace }}/build
+
+      - name: Generate R package requirement list
+        run: |
+          printf 'Package: foo\nVersion: 0.0.1\nLicense: MIT\nDescription: foo\n"' >> DESCRIPTION
+          printf 'Title: foo\nAuthor: me\nMaintainer: me\nImports: wCorr' >> DESCRIPTION
+
+      - name: Install setup-r action dependencies
+        run: brew install imagemagick@6 libgit2
+
+      - name: Setup R
+        uses: r-lib/actions/setup-r@v2
+        with:
+          r-version: '4.2.0'
+
+      - name: Install R packages
+        uses: r-lib/actions/setup-r-dependencies@v2
+        with:
+          cache: true
 
       - name: Run unit tests
         working-directory: ${{ github.workspace }}/build
         run: |
-          ctest --test-dir .             \
-                --schedule-random        \
-                --output-on-failure      \
-                --no-tests=error         \
-                --timeout 120            \
-                --repeat after-timeout:3 \
-                -E '(SciPy)|(wCorr)'     \
-                -j $(nproc)
+          ctest --test-dir .        \
+                --schedule-random   \
+                --output-on-failure \
+                --no-tests=error    \
+                --timeout 60
 
       - name: Inspect MoDLE version
         run: |
diff --git a/.github/workflows/unit-tests-ubuntu.yml b/.github/workflows/unit-tests-ubuntu.yml
index f85d78c6..9b7bfe73 100644
--- a/.github/workflows/unit-tests-ubuntu.yml
+++ b/.github/workflows/unit-tests-ubuntu.yml
@@ -234,7 +234,6 @@ jobs:
                 --output-on-failure      \
                 --no-tests=error         \
                 --timeout 120            \
-                --repeat after-timeout:3 \
                 -j $(nproc)
 
       - name: Run integration tests
@@ -359,12 +358,11 @@ jobs:
       - name: Run unit tests
         working-directory: ${{ github.workspace }}/build
         run: |
-          ctest --test-dir .             \
-                --schedule-random        \
-                --output-on-failure      \
-                --no-tests=error         \
-                --timeout 120            \
-                --repeat after-timeout:3 \
+          ctest --test-dir .        \
+                --schedule-random   \
+                --output-on-failure \
+                --no-tests=error    \
+                --timeout 60        \
                 -j $(nproc)
 
       - name: Run integration tests
diff --git a/Dockerfile b/Dockerfile
index 4db8802e..0e95d316 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -117,7 +117,6 @@ RUN cd "$src_dir/build"           \
          --output-on-failure      \
          --no-tests=error         \
          --timeout 60             \
-         --repeat after-timeout:3 \
 && rm -rf "$src_dir/test/Testing"
 
 ARG FINAL_BASE_IMAGE
diff --git a/conanfile.py b/conanfile.py
index c768488f..65b4ae1b 100644
--- a/conanfile.py
+++ b/conanfile.py
@@ -17,7 +17,7 @@ class MoDLE(ConanFile):
   requires = ["abseil/20211102.0",
               "boost/1.79.0",
               "bzip2/1.0.8",
-              "catch2/2.13.9",
+              "catch2/3.0.1",
               "cli11/2.2.0",
               "concurrentqueue/1.0.3",
               "cpp-sort/1.13.0",
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index b2c70fbe..6b479e55 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -16,8 +16,7 @@ add_executable(test_main)
 
 target_sources(
   test_main
-  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/test_main.cpp
-          ${CMAKE_CURRENT_SOURCE_DIR}/units/common/cli_utils_test.cpp
+  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/units/common/cli_utils_test.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/units/common/const_map_test.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/units/common/dna_test.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/units/contact_matrix/contact_matrix_test.cpp
@@ -32,7 +31,6 @@ target_sources(
           ${CMAKE_CURRENT_SOURCE_DIR}/units/simulation_cpu/simulation_complex_unit_test.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/units/simulation_cpu/simulation_simple_unit_test.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/units/simulation_internal/extrusion_barriers_test.cpp
-          ${CMAKE_CURRENT_SOURCE_DIR}/units/stats/common.hpp
           ${CMAKE_CURRENT_SOURCE_DIR}/units/stats/correlation_test.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/units/stats/correlation_utils_test.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/units/stats/descriptive_test.cpp
@@ -67,7 +65,7 @@ target_link_system_libraries(
   absl::strings
   Boost::filesystem
   Boost::headers
-  Catch2::Catch2
+  Catch2::Catch2WithMain
   fmt::fmt
   HDF5::HDF5
   spdlog::spdlog
@@ -85,12 +83,8 @@ catch_discover_tests(
   " - SHORT"
   WORKING_DIRECTORY
   ${CMAKE_CURRENT_SOURCE_DIR}/../
-  REPORTER
-  xml
   OUTPUT_DIR
   ${CMAKE_CURRENT_SOURCE_DIR}/Testing/
-  OUTPUT_SUFFIX
-  .xml
   EXTRA_ARGS
   --success)
 
@@ -102,12 +96,8 @@ catch_discover_tests(
   " - MEDIUM"
   WORKING_DIRECTORY
   ${CMAKE_CURRENT_SOURCE_DIR}/../
-  REPORTER
-  xml
   OUTPUT_DIR
   ${CMAKE_CURRENT_SOURCE_DIR}/Testing/
-  OUTPUT_SUFFIX
-  .xml
   EXTRA_ARGS
   --success)
 
@@ -119,7 +109,7 @@ catch_discover_tests(
   " - LONG"
   WORKING_DIRECTORY
   ${CMAKE_CURRENT_SOURCE_DIR}/../
-  REPORTER
-  compact # This dramatically speeds up the test suit
+  OUTPUT_DIR
+  ${CMAKE_CURRENT_SOURCE_DIR}/Testing/
   EXTRA_ARGS
   --success)
diff --git a/test/test_main.cpp b/test/test_main.cpp
deleted file mode 100644
index 72a5985b..00000000
--- a/test/test_main.cpp
+++ /dev/null
@@ -1,19 +0,0 @@
-// Copyright (C) 2022 Roberto Rossini <roberros@uio.no>
-//
-// SPDX-License-Identifier: MIT
-
-#define CATCH_CONFIG_MAIN
-
-// clang-format off
-#include "modle/common/suppress_compiler_warnings.hpp"
-DISABLE_WARNING_PUSH
-DISABLE_WARNING_MAYBE_UNINITIALIZED
-#include <catch2/catch.hpp>
-DISABLE_WARNING_POP
-// clang-format on
-
-#include "modle/test/self_deleting_folder.hpp"  // for SelfDeletingFolder
-
-namespace modle::test {
-[[maybe_unused]] inline const SelfDeletingFolder testdir{true};  // NOLINT(cert-err58-cpp)
-}  // namespace modle::test
diff --git a/test/units/common/cli_utils_test.cpp b/test/units/common/cli_utils_test.cpp
index 9097df0c..d03c0c29 100644
--- a/test/units/common/cli_utils_test.cpp
+++ b/test/units/common/cli_utils_test.cpp
@@ -2,7 +2,7 @@
 //
 // SPDX-License-Identifier: MIT
 
-#include <catch2/catch.hpp>  // for operator""_catch_sr, AssertionHandler
+#include <catch2/catch_test_macros.hpp>
 // clang-format off
 #include "modle/common/suppress_compiler_warnings.hpp"
 DISABLE_WARNING_PUSH
diff --git a/test/units/common/const_map_test.cpp b/test/units/common/const_map_test.cpp
index 78281085..aa8657e3 100644
--- a/test/units/common/const_map_test.cpp
+++ b/test/units/common/const_map_test.cpp
@@ -4,7 +4,7 @@
 
 #include "modle/common/const_map.hpp"
 
-#include <catch2/catch.hpp>
+#include <catch2/catch_test_macros.hpp>
 #include <string_view>
 
 namespace modle::test::utils {
diff --git a/test/units/common/dna_test.cpp b/test/units/common/dna_test.cpp
index 36fc6bcb..dbee204c 100644
--- a/test/units/common/dna_test.cpp
+++ b/test/units/common/dna_test.cpp
@@ -4,7 +4,7 @@
 
 #include "modle/common/dna.hpp"
 
-#include <catch2/catch.hpp>
+#include <catch2/catch_test_macros.hpp>
 
 #include "modle/common/common.hpp"
 
diff --git a/test/units/contact_matrix/contact_matrix_test.cpp b/test/units/contact_matrix/contact_matrix_test.cpp
index 6dc3dfa8..5c0db577 100644
--- a/test/units/contact_matrix/contact_matrix_test.cpp
+++ b/test/units/contact_matrix/contact_matrix_test.cpp
@@ -8,7 +8,9 @@
 #include <algorithm>                                // for generate, max
 #include <boost/dynamic_bitset/dynamic_bitset.hpp>  // for dynamic_bitset, dynamic_bitset<>::ref...
 #include <boost/process.hpp>
-#include <catch2/catch.hpp>             // for operator""_catch_sr, AssertionHandler
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_string.hpp>
 #include <filesystem>                   // for path
 #include <stdexcept>                    // for runtime_error
 #include <string>                       // for string
@@ -171,12 +173,14 @@ TEST_CASE("CMatrix in/decrement", "[cmatrix][short]") {
 
   if constexpr (utils::ndebug_not_defined()) {
     CHECK_THROWS_WITH(m.increment(25, 25),
-                      Catch::Contains("Detected an out-of-bound read: attempt to access item at"));
+                      Catch::Matchers::ContainsSubstring(
+                          "Detected an out-of-bound read: attempt to access item at"));
     CHECK(m.get_n_of_missed_updates() == 1);
     CHECK(m.get_tot_contacts() == 1);
 
     CHECK_THROWS_WITH(m.decrement(25, 25),
-                      Catch::Contains("Detected an out-of-bound read: attempt to access item at"));
+                      Catch::Matchers::ContainsSubstring(
+                          "Detected an out-of-bound read: attempt to access item at"));
     CHECK(m.get_n_of_missed_updates() == 1);
     CHECK(m.get_tot_contacts() == 1);
   }
@@ -368,7 +372,7 @@ TEST_CASE("CMatrix blur (SciPy)", "[cmatrix][long]") {
 
     for (usize j = 4; j < input_matrix.nrows(); ++j) {
       for (auto k = j; k < input_matrix.ncols() - 4; ++k) {
-        CHECK(Approx(reference_matrix.get(j, k)) == blurred_matrix.get(j, k));
+        CHECK(Catch::Approx(reference_matrix.get(j, k)) == blurred_matrix.get(j, k));
       }
     }
   }
@@ -408,7 +412,7 @@ TEST_CASE("CMatrix blur parallel (SciPy)", "[cmatrix][long]") {
 
     for (usize j = 4; j < input_matrix.nrows(); ++j) {
       for (auto k = j; k < input_matrix.ncols() - 4; ++k) {
-        CHECK(Approx(reference_matrix.get(j, k)) == blurred_matrix.get(j, k));
+        CHECK(Catch::Approx(reference_matrix.get(j, k)) == blurred_matrix.get(j, k));
       }
     }
   }
@@ -450,8 +454,8 @@ TEST_CASE("CMatrix difference of gaussians (SciPy)", "[cmatrix][long]") {
 
   for (usize j = 4; j < input_matrix.nrows(); ++j) {
     for (auto k = j; k < input_matrix.ncols() - 4; ++k) {
-      CHECK(Approx(reference_matrix.get(j, k)) == gauss_diff_matrix.get(j, k));
-      CHECK(Approx(m1.get(j, k) - m2.get(j, k)) == gauss_diff_matrix.get(j, k));
+      CHECK(Catch::Approx(reference_matrix.get(j, k)) == gauss_diff_matrix.get(j, k));
+      CHECK(Catch::Approx(m1.get(j, k) - m2.get(j, k)) == gauss_diff_matrix.get(j, k));
     }
   }
 }
@@ -496,8 +500,8 @@ TEST_CASE("CMatrix difference of gaussians - parallel (SciPy)", "[cmatrix][long]
 
   for (usize j = 4; j < input_matrix.nrows(); ++j) {
     for (auto k = j; k < input_matrix.ncols() - 4; ++k) {
-      CHECK(Approx(reference_matrix.get(j, k)) == gauss_diff_matrix.get(j, k));
-      CHECK(Approx(m1.get(j, k) - m2.get(j, k)) == gauss_diff_matrix.get(j, k));
+      CHECK(Catch::Approx(reference_matrix.get(j, k)) == gauss_diff_matrix.get(j, k));
+      CHECK(Catch::Approx(m1.get(j, k) - m2.get(j, k)) == gauss_diff_matrix.get(j, k));
     }
   }
 }
diff --git a/test/units/include/modle/test/self_deleting_folder.hpp b/test/units/include/modle/test/self_deleting_folder.hpp
index 7871eb48..2944b3a5 100644
--- a/test/units/include/modle/test/self_deleting_folder.hpp
+++ b/test/units/include/modle/test/self_deleting_folder.hpp
@@ -27,8 +27,6 @@ class SelfDeletingFolder {
       tmpdir = std::filesystem::temp_directory_path();
     } catch (const std::filesystem::filesystem_error& e) {
       // Workaround spurious CI failures due to missing /tmp folder exception
-      assert(absl::StartsWith(e.what(), "std::filesystem::temp_directory_path: Not a directory"));
-      assert(absl::EndsWith(e.what(), "\"/tmp\""));
       tmpdir = "test/data/unit_tests/scratch";
     }
 
diff --git a/test/units/interval_tree/interval_tree_test.cpp b/test/units/interval_tree/interval_tree_test.cpp
index 9d202a9e..238d6483 100644
--- a/test/units/interval_tree/interval_tree_test.cpp
+++ b/test/units/interval_tree/interval_tree_test.cpp
@@ -7,12 +7,12 @@
 #include <absl/strings/str_split.h>  // for StrSplit, Splitter, ByAnyChar
 #include <fmt/format.h>              // for format
 
-#include <algorithm>         // for max
-#include <catch2/catch.hpp>  // for AssertionHandler, operator""_catch_sr, SourceLineInfo
-#include <filesystem>        // for path
-#include <string>            // for string, basic_string, operator==
-#include <string_view>       // for string_view
-#include <vector>            // for vector, allocator
+#include <algorithm>  // for max
+#include <catch2/catch_test_macros.hpp>
+#include <filesystem>   // for path
+#include <string>       // for string, basic_string, operator==
+#include <string_view>  // for string_view
+#include <vector>       // for vector, allocator
 
 #include "modle/common/numeric_utils.hpp"         // for parse_numeric_or_throw
 #include "modle/compressed_io/compressed_io.hpp"  // for Reader
diff --git a/test/units/libmodle_io/bed_parser_test.cpp b/test/units/libmodle_io/bed_parser_test.cpp
index eca96df7..fe978481 100644
--- a/test/units/libmodle_io/bed_parser_test.cpp
+++ b/test/units/libmodle_io/bed_parser_test.cpp
@@ -5,12 +5,12 @@
 #include <absl/strings/str_split.h>  // for StrSplit, Splitter
 #include <fmt/format.h>              // for to_string
 
-#include <algorithm>         // for sort, max
-#include <catch2/catch.hpp>  // for AssertionHandler, operator""_catch_sr, SourceLineInfo
-#include <filesystem>        // for path
-#include <string>            // for string, basic_string, operator==, char_traits, stoull
-#include <string_view>       // for operator!=, basic_string_view, string_view, operator<
-#include <vector>            // for vector
+#include <algorithm>  // for sort, max
+#include <catch2/catch_test_macros.hpp>
+#include <filesystem>   // for path
+#include <string>       // for string, basic_string, operator==, char_traits, stoull
+#include <string_view>  // for operator!=, basic_string_view, string_view, operator<
+#include <vector>       // for vector
 
 #include "absl/strings/match.h"                   // for StrContains
 #include "modle/bed/bed.hpp"                      // for BED, Parser, formatter<>::format, BED::BED3
diff --git a/test/units/libmodle_io/bed_tree_test.cpp b/test/units/libmodle_io/bed_tree_test.cpp
index d117b0f7..1a242e3f 100644
--- a/test/units/libmodle_io/bed_tree_test.cpp
+++ b/test/units/libmodle_io/bed_tree_test.cpp
@@ -6,12 +6,12 @@
 #include <absl/types/span.h>         // for Span
 #include <fmt/format.h>              // for format
 
-#include <algorithm>         // for max, min
-#include <catch2/catch.hpp>  // for AssertionHandler, operator""_catch_sr, SourceLineInfo
-#include <filesystem>        // for path
-#include <string>            // for basic_string, operator==, string
-#include <string_view>       // for string_view
-#include <vector>            // for vector
+#include <algorithm>  // for max, min
+#include <catch2/catch_test_macros.hpp>
+#include <filesystem>   // for path
+#include <string>       // for basic_string, operator==, string
+#include <string_view>  // for string_view
+#include <vector>       // for vector
 
 #include "modle/bed/bed.hpp"        // for BED, BED_tree, BED_tree::contains, BED_tree::count_...
 #include "modle/common/common.hpp"  // for usize, u8
diff --git a/test/units/libmodle_io/compressed_io_test.cpp b/test/units/libmodle_io/compressed_io_test.cpp
index 8ec6e939..295feb34 100644
--- a/test/units/libmodle_io/compressed_io_test.cpp
+++ b/test/units/libmodle_io/compressed_io_test.cpp
@@ -4,11 +4,11 @@
 
 #include "modle/compressed_io/compressed_io.hpp"  // for Reader
 
-#include <catch2/catch.hpp>  // for AssertionHandler, operator""_catch_sr, Source...
-#include <filesystem>        // for path, operator/
-#include <fstream>           // for ifstream, basic_ios, basic_istream, operator<<
-#include <iostream>          // for cerr
-#include <string>            // for operator==, string, basic_string, getline
+#include <catch2/catch_test_macros.hpp>
+#include <filesystem>  // for path, operator/
+#include <fstream>     // for ifstream, basic_ios, basic_istream, operator<<
+#include <iostream>    // for cerr
+#include <string>      // for operator==, string, basic_string, getline
 
 #include "modle/test/self_deleting_folder.hpp"  // for SelfDeletingFolder
 
diff --git a/test/units/libmodle_io/cooler_test.cpp b/test/units/libmodle_io/cooler_test.cpp
index 1c91e556..6e51bf27 100644
--- a/test/units/libmodle_io/cooler_test.cpp
+++ b/test/units/libmodle_io/cooler_test.cpp
@@ -11,13 +11,14 @@
 #include <spdlog/sinks/null_sink.h>  // for null_sink_mt
 #include <spdlog/spdlog.h>           // for set_default_logger
 
-#include <algorithm>         // for max, max_element, transform
-#include <catch2/catch.hpp>  // for operator""_catch_sr, AssertionHandler, Source...
-#include <exception>         // for exception
-#include <filesystem>        // for operator/, path
-#include <memory>            // for make_shared, allocator_traits<>::value_type
-#include <stdexcept>         // for runtime_error
-#include <string_view>       // for string_view
+#include <algorithm>  // for max, max_element, transform
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_string.hpp>
+#include <exception>    // for exception
+#include <filesystem>   // for operator/, path
+#include <memory>       // for make_shared, allocator_traits<>::value_type
+#include <stdexcept>    // for runtime_error
+#include <string_view>  // for string_view
 
 #include "modle/common/common.hpp"                // for u64, u32, usize, i64, u8
 #include "modle/common/utils.hpp"                 // for parse_numeric_or_throw
@@ -51,7 +52,7 @@ TEST_CASE("cooler ctor", "[io][cooler][short]") {
   CHECK(Cooler::validate_file_format(f, Cooler::FLAVOR::COOL));
   CHECK_THROWS_WITH(!Cooler::validate_file_format(f, Cooler::FLAVOR::MCOOL,
                                                   Cooler::IO_MODE::READ_ONLY, 1'000'000),
-                    Catch::Matchers::Contains("Expected format flavor MCOOL, found COOL"));
+                    Catch::Matchers::ContainsSubstring("Expected format flavor MCOOL, found COOL"));
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
diff --git a/test/units/libmodle_io/hdf5_test.cpp b/test/units/libmodle_io/hdf5_test.cpp
index e4d4b778..89309ddb 100644
--- a/test/units/libmodle_io/hdf5_test.cpp
+++ b/test/units/libmodle_io/hdf5_test.cpp
@@ -11,14 +11,16 @@
 #include <H5Spublic.h>     // for H5S_UNLIMITED
 #include <fmt/format.h>    // for format
 
-#include <algorithm>         // for max
-#include <catch2/catch.hpp>  // for operator""_catch_sr, AssertionHandler, Source...
-#include <filesystem>        // for operator/, path
-#include <limits>            // for numeric_limits
-#include <string>            // for string, basic_string, allocator, operator==
-#include <string_view>       // for string_view
-#include <tuple>             // for ignore
-#include <vector>            // for vector
+#include <algorithm>  // for max
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_contains.hpp>
+#include <catch2/matchers/catch_matchers_string.hpp>
+#include <filesystem>   // for operator/, path
+#include <limits>       // for numeric_limits
+#include <string>       // for string, basic_string, allocator, operator==
+#include <string_view>  // for string_view
+#include <tuple>        // for ignore
+#include <vector>       // for vector
 
 #include "H5Ppublic.h"              // for H5F_ACC_TRUNC, H5T_CSET_ASCII, H5T_STR_NULLPAD
 #include "modle/common/common.hpp"  // for i64, usize
@@ -169,11 +171,11 @@ TEST_CASE("check_dataset_type HDF5", "[io][hdf5][short]") {
   CHECK(!check_dataset_type(str_dataset, H5::PredType::NATIVE_INT64, false));
 
   CHECK_THROWS_WITH(check_dataset_type(int_dataset, H5::PredType::NATIVE_INT),
-                    Catch::Matchers::Contains("incorrect datasize"));
+                    Catch::Matchers::ContainsSubstring("incorrect datasize"));
   CHECK_THROWS_WITH(check_dataset_type(int_dataset, H5::PredType::NATIVE_FLOAT),
                     Catch::Matchers::EndsWith("incorrect datatype"));
   CHECK_THROWS_WITH(check_dataset_type(int_dataset, H5::PredType::NATIVE_UINT64),
-                    Catch::Matchers::Contains("incorrect signedness"));
+                    Catch::Matchers::ContainsSubstring("incorrect signedness"));
 
   std::filesystem::remove_all(testdir());
   if (const auto& p = testdir().parent_path(); std::filesystem::is_empty(p)) {
@@ -192,8 +194,9 @@ TEST_CASE("read_write_string HDF5 - long string", "[io][hdf5][short]") {
   auto dataset = init_test_str_dataset(f);
 
 #ifndef NDEBUG
-  CHECK_THROWS_WITH(write_str(s, dataset, init_str_type(), 0),
-                    Catch::Contains("string does not fit in the receiving dataset"));
+  CHECK_THROWS_WITH(
+      write_str(s, dataset, init_str_type(), 0),
+      Catch::Matchers::ContainsSubstring("string does not fit in the receiving dataset"));
 #else
   CHECK(write_str(s, dataset, init_str_type(), 0) == 1);
 #endif
diff --git a/test/units/simulation_cpu/collision_encoding_test.cpp b/test/units/simulation_cpu/collision_encoding_test.cpp
index 9e778b83..c88bfb93 100644
--- a/test/units/simulation_cpu/collision_encoding_test.cpp
+++ b/test/units/simulation_cpu/collision_encoding_test.cpp
@@ -6,7 +6,9 @@
 
 #include <fmt/format.h>
 
-#include <catch2/catch.hpp>  // for SourceLineInfo, StringRef, TEST_CASE
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_string.hpp>
+
 namespace modle::test::libmodle {
 
 [[maybe_unused]] static std::string event_to_str(const CollisionEvent<> event) {
diff --git a/test/units/simulation_cpu/common.hpp b/test/units/simulation_cpu/common.hpp
index 686a3aee..ad2e19ec 100644
--- a/test/units/simulation_cpu/common.hpp
+++ b/test/units/simulation_cpu/common.hpp
@@ -8,7 +8,7 @@
 #include <fmt/format.h>
 
 #include <algorithm>
-#include <catch2/catch.hpp>
+#include <catch2/catch_test_macros.hpp>
 #include <cstddef>
 #include <vector>
 
diff --git a/test/units/simulation_cpu/simulation_complex_unit_test.cpp b/test/units/simulation_cpu/simulation_complex_unit_test.cpp
index 117490eb..a703d165 100644
--- a/test/units/simulation_cpu/simulation_complex_unit_test.cpp
+++ b/test/units/simulation_cpu/simulation_complex_unit_test.cpp
@@ -6,10 +6,10 @@
 
 #include <boost/dynamic_bitset/dynamic_bitset.hpp>  // for dynamic_bitset
 #include <cassert>                                  // for assert
-#include <catch2/catch.hpp>                         // for SourceLineInfo, StringRef, TEST_CASE
-#include <memory>                                   // for allocator, allocator_traits<>::value_...
-#include <string_view>                              // for string_view
-#include <vector>                                   // for vector
+#include <catch2/catch_test_macros.hpp>
+#include <memory>       // for allocator, allocator_traits<>::value_...
+#include <string_view>  // for string_view
+#include <vector>       // for vector
 
 #include "./common.hpp"                        // for construct_lef, NO_COLLISION, check_si...
 #include "modle/common/common.hpp"             // for bp_t, usize
diff --git a/test/units/simulation_cpu/simulation_simple_unit_test.cpp b/test/units/simulation_cpu/simulation_simple_unit_test.cpp
index 204f893e..ba4cd74d 100644
--- a/test/units/simulation_cpu/simulation_simple_unit_test.cpp
+++ b/test/units/simulation_cpu/simulation_simple_unit_test.cpp
@@ -4,14 +4,14 @@
 
 #include <absl/types/span.h>  // for MakeSpan, Span, MakeConstSpan
 
-#include <algorithm>         // for all_of, equal, copy, for_each, gene...
-#include <cassert>           // for assert
-#include <catch2/catch.hpp>  // for AssertionHandler, operator""_catch_sr
-#include <iterator>          // for back_insert_iterator, back_inserter
-#include <memory>            // for allocator, allocator_traits<>::valu...
-#include <numeric>           // for iota
-#include <string_view>       // for string_view
-#include <vector>            // for vector
+#include <algorithm>  // for all_of, equal, copy, for_each, gene...
+#include <cassert>    // for assert
+#include <catch2/catch_test_macros.hpp>
+#include <iterator>     // for back_insert_iterator, back_inserter
+#include <memory>       // for allocator, allocator_traits<>::valu...
+#include <numeric>      // for iota
+#include <string_view>  // for string_view
+#include <vector>       // for vector
 
 #include "./common.hpp"                        // for construct_lef, NO_COLLISION, requir...
 #include "modle/common/common.hpp"             // for usize, bp_t
diff --git a/test/units/simulation_internal/extrusion_barriers_test.cpp b/test/units/simulation_internal/extrusion_barriers_test.cpp
index d99b7a4d..3f0e0db4 100644
--- a/test/units/simulation_internal/extrusion_barriers_test.cpp
+++ b/test/units/simulation_internal/extrusion_barriers_test.cpp
@@ -4,7 +4,8 @@
 
 #include "modle/extrusion_barriers.hpp"  // for ExtrusionBarrier
 
-#include <catch2/catch.hpp>  // for AssertionHandler, operator""_catch_sr
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
 
 #include "modle/common/common.hpp"  // for dna::Direction
 
@@ -46,8 +47,8 @@ TEST_CASE("Extrusion barriers - occupancy", "[barriers][simulation][short]") {
                                           {stp_inactive1, stp_inactive2},
                                           {State::ACTIVE, State::ACTIVE}};
 
-  CHECK(barriers.occupancy(0) == Approx(occupancy1));
-  CHECK(barriers.occupancy(1) == Approx(occupancy2));
+  CHECK(barriers.occupancy(0) == Catch::Approx(occupancy1));
+  CHECK(barriers.occupancy(1) == Catch::Approx(occupancy2));
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
@@ -79,7 +80,8 @@ TEST_CASE("Extrusion barriers - compute_occupancy", "[barriers][simulation][shor
 
   stp_active = 0.7;
   stp_inactive = 0.5;
-  CHECK(ExtrusionBarrier::compute_occupancy_from_stp(stp_active, stp_inactive) == Approx(0.625));
+  CHECK(ExtrusionBarrier::compute_occupancy_from_stp(stp_active, stp_inactive) ==
+        Catch::Approx(0.625));
 
   stp_active = 1.0;
   stp_inactive = 0.5;
diff --git a/test/units/stats/common.hpp b/test/units/stats/common.hpp
deleted file mode 100644
index c7dad974..00000000
--- a/test/units/stats/common.hpp
+++ /dev/null
@@ -1,237 +0,0 @@
-// Copyright (C) 2022 Roberto Rossini <roberros@uio.no>
-//
-// SPDX-License-Identifier: MIT
-
-#pragma once
-#include <absl/container/flat_hash_set.h>  // for flat_hash_set, operator!=
-#include <absl/hash/hash.h>                // for Hash
-#include <fmt/compile.h>
-#include <fmt/format.h>
-
-#include <algorithm>  // for max, generate
-#include <array>      // for array
-#include <boost/container_hash/hash.hpp>
-#include <boost/process.hpp>
-#include <cstdio>      // for pclose, fgets, popen, FILE
-#include <filesystem>  // for create_directories, exists, path, remove
-#include <fstream>     // for basic_ofstream, operator<<, basic_ostream
-#include <memory>      // for allocator, unique_ptr
-#include <numeric>     // for iota
-#include <stdexcept>   // for runtime_error
-#include <string>      // for string, operator+, char_traits, stod
-#include <string_view>
-#include <type_traits>  // for enable_if, is_arithmetic
-#include <utility>      // for pair
-#include <vector>       // for vector
-
-#include "modle/common/common.hpp"  // for u32, i32, i64
-#include "modle/common/const_map.hpp"
-#include "modle/common/numeric_utils.hpp"
-#include "modle/common/random.hpp"
-
-namespace modle::test::stats {
-using namespace std::string_view_literals;
-
-// clang-format off
-#define MAKE_CORR_TEST_CASE(METHOD)                           \
-  std::make_pair(std::string_view{METHOD},                    \
-     "#!/usr/bin/env python3\n"                               \
-     "from scipy.stats import " METHOD "r\n"                  \
-     "from numpy import fromstring\n"                         \
-     "import fileinput\n"                                     \
-                                                              \
-     "for line in fileinput.input():\n"                       \
-     "  v1, _, v2 = line.strip().partition(\"\\t\")\n"        \
-     "  v1 = fromstring(v1, sep=\",\")\n"                     \
-     "  v2 = fromstring(v2, sep=\",\")\n"                     \
-                                                              \
-     "  corr, pv = " METHOD "r(v1, v2)\n"                     \
-     "  print(f\"{corr:.16e}\\t{pv:.16e}\", flush=True)\n"sv)
-
-#define MAKE_WEIGHTED_CORR_TEST_CASE(METHOD)                       \
-  std::make_pair(std::string_view{METHOD},                         \
-     "#!/usr/bin/env Rscript\n"                                    \
-     "library(\"wCorr\")\n"                                        \
-     "f <- file(\"stdin\")\n"                                      \
-     "open(f)\n"                                                   \
-     "method=sub(\"weighted_\", \"\", \"" METHOD "\")\n"           \
-     "while(length(line <- readLines(f, n=1)) > 0) {\n"            \
-     "   toks <- strsplit(line, '\\t')[[1]]\n"                     \
-     "   v1 <- as.numeric(strsplit(toks[[1]], ',')[[1]])\n"        \
-     "   v2 <- as.numeric(strsplit(toks[[2]], ',')[[1]])\n"        \
-     "   w <- as.numeric(strsplit(toks[[3]], ',')[[1]])\n"         \
-                                                                   \
-     "   corr <- weightedCorr(v1, v2, weights=w, method=method)\n" \
-     "   write(sprintf(\"%.22f\\t-1\", corr), stdout())\n"         \
-     "   flush.console()\n"                                        \
-     "}"sv)
-
-static constexpr utils::ConstMap<std::string_view, std::string_view, 4> external_cmds{
-    MAKE_CORR_TEST_CASE("pearson"),
-    MAKE_CORR_TEST_CASE("spearman"),
-    MAKE_WEIGHTED_CORR_TEST_CASE("weighted_pearson"),
-    MAKE_WEIGHTED_CORR_TEST_CASE("weighted_spearman")
-};  // clang-format on
-
-template <class N, class = std::enable_if_t<std::is_arithmetic_v<N>>>
-inline void generate_random_vect(random::PRNG_t& rand_eng, std::vector<N>& buff, N min, N max,
-                                 bool allow_duplicates = true) {
-  using dist_t =
-      typename std::conditional<std::is_floating_point_v<N>, random::uniform_real_distribution<N>,
-                                random::uniform_int_distribution<N>>::type;
-  dist_t dist(min, max);
-  if (allow_duplicates) {
-    std::generate(buff.begin(), buff.end(), [&]() { return dist(rand_eng); });
-  } else {
-    absl::flat_hash_set<u32> s;
-    while (s.size() < buff.size()) {
-      s.insert(dist(rand_eng));
-    }
-    std::copy(s.begin(), s.end(), buff.begin());
-  }
-}
-
-template <class N = u32, class = std::enable_if_t<std::is_arithmetic_v<N>>>
-[[nodiscard]] inline std::vector<N> generate_random_vect(random::PRNG_t& rand_eng, usize size,
-                                                         N min, N max,
-                                                         bool allow_duplicates = true) {
-  std::vector<N> v(size);
-  generate_random_vect(rand_eng, v, min, max, allow_duplicates);
-  return v;
-}
-
-inline std::pair<std::vector<u32>, std::vector<u32>> generate_correlated_vects(
-    random::PRNG_t& rand_eng, u32 size) {
-  random::uniform_int_distribution<i32> dist(static_cast<i32>(size) / -50,
-                                             static_cast<i32>(size / 50));
-  std::vector<u32> v1(size);
-  std::vector<u32> v2(size);
-  std::iota(v1.begin(), v1.end(), 0);
-  std::iota(v2.begin(), v2.end(), 0);
-  for (usize i = 0; i < size; ++i) {
-    auto n = static_cast<i64>(v1[i]) + dist(rand_eng);
-    v1[i] = static_cast<u32>(std::max(i64(0), n));
-    n = static_cast<i64>(v2[i]) + dist(rand_eng);
-    v2[i] = static_cast<u32>(std::max(i64(0), n));
-  }
-  return {v1, v2};
-}
-
-template <typename N, typename = typename std::enable_if<std::is_arithmetic<N>::value, N>::type>
-inline void run_external_corr(std::string_view method, std::vector<N>& v1, std::vector<N>& v2,
-                              std::vector<double>& weights, std::atomic<double>& rho,
-                              std::atomic<double>& pv, std::mutex& data_mutex,
-                              std::condition_variable& input_data_cv,
-                              std::condition_variable& output_data_cv,
-                              std::atomic<bool>& input_data_ready,
-                              std::atomic<bool>& output_data_ready) {
-  boost::process::ipstream stdout_stream;
-  boost::process::ipstream stderr_stream;
-  boost::process::opstream stdin_stream;
-
-  try {
-    auto c = [&]() {
-      if (absl::StartsWith(method, "weighted_"sv)) {
-        return boost::process::child(
-            boost::process::search_path("Rscript").string(), "--quiet", "-e",
-            std::string{external_cmds.at(method)},
-            boost::process::std_in<stdin_stream, boost::process::std_out> stdout_stream,
-            boost::process::std_err > stderr_stream);
-      }
-
-      return boost::process::child(
-          boost::process::search_path("python3").string(), "-c",
-          std::string{external_cmds.at(method)},
-          boost::process::std_in<stdin_stream, boost::process::std_out> stdout_stream,
-          boost::process::std_err > stderr_stream);
-    }();
-    assert(c.running());
-
-    std::string sbuff;
-    while (true) {
-      {
-        std::unique_lock<std::mutex> l(data_mutex);
-        input_data_cv.wait(l, [&]() { return input_data_ready.load(); });
-        input_data_ready = false;
-      }
-
-      if (v1.empty()) {  // EOQ signal
-        assert(v2.empty());
-        stdin_stream.pipe().close();
-        c.wait();
-        return;
-      }
-
-      sbuff = [&]() {
-        if (absl::StartsWith(method, "weighted_")) {
-          return fmt::format(FMT_COMPILE("{}\t{}\t{}\n"), fmt::join(v1, ","), fmt::join(v2, ","),
-                             fmt::join(weights, ","));
-        }
-        return fmt::format(FMT_COMPILE("{}\t{}\n"), fmt::join(v1, ","), fmt::join(v2, ","));
-      }();
-      stdin_stream.write(sbuff.data(), static_cast<std::streamsize>(sbuff.size()));
-      stdin_stream.flush();
-
-      std::getline(stdout_stream, sbuff);
-      const auto sep_idx = sbuff.find('\t');
-      assert(sep_idx < sbuff.size());
-      rho = utils::parse_numeric_or_throw<double>(std::string_view(sbuff.data(), sep_idx));
-      pv = utils::parse_numeric_or_throw<double>(
-          std::string_view(sbuff.data() + sep_idx + 1, sbuff.size() - sep_idx));
-      output_data_ready = true;
-      output_data_cv.notify_one();
-    }
-  } catch (const std::exception& e) {
-    std::string buff1;
-    std::string buff2;
-    while (std::getline(stderr_stream, buff1)) {
-      buff2.append(buff1);
-      buff2.append("\n");
-    }
-
-    if (!buff2.empty()) {
-      throw std::runtime_error(fmt::format(FMT_STRING("{}:\n{}"), e.what(), buff2));
-    }
-    throw;
-  }
-}
-
-template <typename N, typename = typename std::enable_if<std::is_arithmetic<N>::value, N>::type>
-[[nodiscard]] inline std::pair<double, double> external_corr(const std::vector<N>& v1,
-                                                             const std::vector<N>& v2,
-                                                             const std::vector<double>& w,
-                                                             std::string_view method) {
-  boost::process::ipstream stdout_stream;
-  boost::process::opstream stdin_stream;
-  if (absl::StartsWith(method, "weighted_"sv)) {
-    boost::process::spawn(
-        boost::process::search_path("Rscript").string(), "--quiet", "-e",
-        std::string{external_cmds.at(method)},
-        boost::process::std_in<stdin_stream, boost::process::std_out> stdout_stream);
-  } else {
-    boost::process::spawn(
-        boost::process::search_path("python3").string(), "-c",
-        std::string{external_cmds.at(method)},
-        boost::process::std_in<stdin_stream, boost::process::std_out> stdout_stream);
-  }
-
-  std::string buff;
-  if (absl::StartsWith(method, "weighted_"sv)) {
-    buff = fmt::format(FMT_COMPILE("{}\t{}\t{}\n"), fmt::join(v1, ","), fmt::join(v2, ","),
-                       fmt::join(w, ","));
-  } else {
-    buff = fmt::format(FMT_COMPILE("{}\t{}\n"), fmt::join(v1, ","), fmt::join(v2, ","));
-  }
-  stdin_stream.write(buff.data(), static_cast<std::streamsize>(buff.size()));
-  stdin_stream.close();
-
-  std::getline(stdout_stream, buff);
-  const auto sep_idx = buff.find('\t');
-  assert(sep_idx < buff.size());
-  const auto rho = utils::parse_numeric_or_throw<double>(std::string_view(buff.data(), sep_idx));
-  const auto pv = utils::parse_numeric_or_throw<double>(
-      std::string_view(buff.data() + sep_idx + 1, buff.size() - sep_idx));
-
-  return {rho, pv};
-}
-}  // namespace modle::test::stats
diff --git a/test/units/stats/correlation_test.cpp b/test/units/stats/correlation_test.cpp
index 4f01b516..31e59576 100644
--- a/test/units/stats/correlation_test.cpp
+++ b/test/units/stats/correlation_test.cpp
@@ -6,102 +6,298 @@
 
 #include <absl/container/flat_hash_set.h>  // for BitMask, operator!=
 #include <absl/strings/match.h>            // for StartsWith
+#include <fmt/compile.h>                   // for FMT_COMPILE
 #include <fmt/format.h>                    // for format
 
-#include <atomic>                         // for atomic
-#include <boost/exception/exception.hpp>  // for clone_base
-#include <cassert>                        // for assert
-#include <catch2/catch.hpp>               // for Approx, operator==, AssertionHandler, operator...
-#include <condition_variable>             // for condition_variable
-#include <exception>                      // for exception
-#include <mutex>                          // for mutex, unique_lock
-#include <random>                         // for random_device
-#include <stdexcept>                      // for overflow_error
-#include <string>                         // for char_traits
-#include <string_view>                    // for operator==, operator""sv, basic_string_view
-#include <thread>                         // for thread
-#include <utility>                        // for make_pair
-#include <vector>                         // for vector
-
-#include "./common.hpp"             // for generate_random_vect, external_corr, run_exter...
+#include <boost/exception/exception.hpp>
+#include <boost/process.hpp>
+#include <cassert>  // for assert
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <exception>    // for exception
+#include <ostream>      // for ostream
+#include <random>       // for random_device
+#include <string>       // for char_traits
+#include <string_view>  // for operator==, operator""sv, basic_string_view
+#include <utility>      // for make_pair
+#include <vector>       // for vector
+
 #include "modle/common/common.hpp"  // for u32, usize
+#include "modle/common/const_map.hpp"
+#include "modle/common/numeric_utils.hpp"
 #include "modle/common/random.hpp"  // for PRNG_t
-#include "modle/test/self_deleting_folder.hpp"  // for SelfDeletingFolder
+#include "modle/test/self_deleting_folder.hpp"
+
+namespace modle::test {
+inline const SelfDeletingFolder testdir{true};  // NOLINT(cert-err58-cpp)
+}  // namespace modle::test
 
 namespace modle::test::stats {
 using namespace modle::stats;
+using namespace std::string_view_literals;
+
+// clang-format off
+#define MAKE_CORR_TEST_CASE(METHOD)                           \
+  std::make_pair(std::string_view{METHOD},                    \
+     "#!/usr/bin/env python3\n"                               \
+     "from scipy.stats import " METHOD "r\n"                  \
+     "from numpy import fromstring\n"                         \
+     "import fileinput\n"                                     \
+                                                              \
+     "for line in fileinput.input():\n"                       \
+     "  v1, _, v2 = line.strip().partition(\"\\t\")\n"        \
+     "  v1 = fromstring(v1, sep=\",\")\n"                     \
+     "  v2 = fromstring(v2, sep=\",\")\n"                     \
+                                                              \
+     "  corr, pv = " METHOD "r(v1, v2)\n"                     \
+     "  print(f\"{corr:.16e}\\t{pv:.16e}\", flush=True)\n"sv)
+
+#define MAKE_WEIGHTED_CORR_TEST_CASE(METHOD)                       \
+  std::make_pair(std::string_view{METHOD},                         \
+     "#!/usr/bin/env Rscript\n"                                    \
+     "library(\"wCorr\")\n"                                        \
+     "f <- file(\"stdin\")\n"                                      \
+     "open(f)\n"                                                   \
+     "method <- sub(\"weighted_\", \"\", \"" METHOD "\")\n"        \
+     "while(length(line <- readLines(f, n=1)) > 0) {\n"            \
+     "   toks <- strsplit(line, '\\t')[[1]]\n"                     \
+     "   v1 <- as.numeric(strsplit(toks[[1]], ',')[[1]])\n"        \
+     "   v2 <- as.numeric(strsplit(toks[[2]], ',')[[1]])\n"        \
+     "   w <- as.numeric(strsplit(toks[[3]], ',')[[1]])\n"         \
+                                                                   \
+     "   corr <- weightedCorr(v1, v2, weights=w, method=method)\n" \
+     "   write(sprintf(\"%.22f\\t-1\", corr), stdout())\n"         \
+     "   flush.console()\n"                                        \
+     "}"sv)
+
+static constexpr utils::ConstMap<std::string_view, std::string_view, 4> external_cmds{
+    MAKE_CORR_TEST_CASE("pearson"),
+    MAKE_CORR_TEST_CASE("spearman"),
+    MAKE_WEIGHTED_CORR_TEST_CASE("weighted_pearson"),
+    MAKE_WEIGHTED_CORR_TEST_CASE("weighted_spearman")
+};  // clang-format on
+
+template <class N>
+void generate_random_vector(random::PRNG_t& rand_eng, std::vector<N>& buff, N min, N max,
+                            bool allow_duplicates = true) {
+  using DistrT =
+      typename std::conditional<std::is_floating_point_v<N>, random::uniform_real_distribution<N>,
+                                random::uniform_int_distribution<N>>::type;
+  DistrT dist(min, max);
+  if (allow_duplicates) {
+    std::generate(buff.begin(), buff.end(), [&]() { return dist(rand_eng); });
+    return;
+  }
+
+  absl::flat_hash_set<u32> s(buff.size());
+  while (s.size() < buff.size()) {
+    s.insert(dist(rand_eng));
+  }
+  std::copy(s.begin(), s.end(), buff.begin());
+}
 
-template <typename N>
-static void test_correlation_w_random_vector(std::string_view method, usize vector_size,
-                                             usize iterations, N min, N max) {
-  static_assert(std::is_arithmetic_v<N>, "N should be an arithmetic type.");
-  random::PRNG_t rand_eng{std::random_device{}()};
-
-  std::mutex data_mutex;
-  std::condition_variable input_data_cv;
-  std::condition_variable output_data_cv;
-  std::atomic<bool> input_data_ready{false};
-  std::atomic<bool> output_data_ready{false};
-
-  std::vector<N> v1(vector_size);
-  std::vector<N> v2(vector_size);
-  std::vector<double> weights(vector_size);
-
-  std::atomic<double> cfx_reference{};
-  std::atomic<double> pv_reference{};
-
-  std::thread t([&]() {
-    run_external_corr(method, v1, v2, weights, cfx_reference, pv_reference, data_mutex,
-                      input_data_cv, output_data_cv, input_data_ready, output_data_ready);
-  });
-
-  auto pearson = Pearson<>{};
-  auto spearman = Spearman<>{};
-  for (usize i = 0; i < iterations; ++i) {
-    assert(!input_data_ready);
-    generate_random_vect(rand_eng, v1, min, max);
-    generate_random_vect(rand_eng, v2, min, max);
-    if (absl::StartsWith(method, "weighted_")) {
-      generate_random_vect(rand_eng, weights, 0.0, 1.0);
-      // fmt::print(stderr, FMT_STRING("{}\n"), fmt::join(weights, ","));
+template <class N>
+[[nodiscard]] std::vector<N> generate_random_vector(random::PRNG_t& rand_eng, usize size, N min,
+                                                    N max, bool allow_duplicates = true) {
+  std::vector<N> v(size);
+  generate_random_vector(rand_eng, v, min, max, allow_duplicates);
+  return v;
+}
+
+template <class N1, class N2>
+class CorrelationBuff {
+  static_assert(std::is_arithmetic_v<N1>);
+  static_assert(std::is_arithmetic_v<N2>);
+  std::vector<N1> _v1;
+  std::vector<N1> _v2;
+  std::vector<N2> _weights;
+
+ public:
+  CorrelationBuff() = delete;
+  CorrelationBuff(usize size, bool with_weights)
+      : _v1(size, N1(0)), _v2(size, N1(0)), _weights(with_weights ? size : 0, N2(0)) {}
+  [[nodiscard]] const std::vector<N1>& v1() const noexcept { return this->_v1; }
+  [[nodiscard]] const std::vector<N1>& v2() const noexcept { return this->_v2; }
+  [[nodiscard]] const std::vector<N2>& weights() const noexcept { return this->_weights; }
+
+  void generate(random::PRNG_t& rand_eng, N1 min, N1 max) {
+    assert(this->_v1.size() == this->_v2.size());
+    generate_random_vector(rand_eng, this->_v1, min, max);
+    generate_random_vector(rand_eng, this->_v2, min, max);
+    assert(this->_v1.size() == this->_v2.size());
+
+    if (!this->_weights.empty()) {
+      assert(this->_v1.size() == this->_weights.size());
+      generate_random_vector(rand_eng, this->_weights, N2(0), N2(1));
+      assert(this->_v1.size() == this->_weights.size());
     }
-    input_data_ready = true;
-    input_data_cv.notify_one();
-    {
-      std::unique_lock<std::mutex> l(data_mutex);
-      output_data_cv.wait(l, [&]() { return output_data_ready.load(); });
-      output_data_ready = false;
+  }
+
+  void serialize(std::string& buff) const {
+    if (!this->_weights.empty()) {
+      buff = fmt::format(FMT_COMPILE("{}\t{}\t{}\n"), fmt::join(this->_v1, ","),
+                         fmt::join(this->_v2, ","), fmt::join(this->_weights, ","));
+      return;
     }
 
-    const auto [cfx, pv] = [&]() {
-      if (method == "pearson"sv) {
-        const auto res = pearson(v1, v2);
-        return std::make_pair(res.pcc, res.pvalue);
-      }
+    buff =
+        fmt::format(FMT_COMPILE("{}\t{}\n"), fmt::join(this->_v1, ","), fmt::join(this->_v2, ","));
+  }
+  [[nodiscard]] std::string serialize() const {
+    std::string buff;
+    this->serialize(buff);
+    return buff;
+  }
+};
+
+template <class N1, class N2>
+[[nodiscard]] static std::pair<double, double> compute_correlation(
+    std::string_view method, const CorrelationBuff<N1, N2>& data) {
+  const auto is_weighted = absl::StartsWith(method, "weighted_");
+  const auto is_pearson = absl::EndsWith(method, "pearson");
+  const auto is_spearman = absl::EndsWith(method, "spearman");
+
+  if (is_pearson && is_weighted) {
+    const auto res = Pearson<>{}(data.v1(), data.v2(), data.weights());
+    return std::make_pair(res.pcc, res.pvalue);
+  }
 
-      if (method == "weighted_pearson"sv) {
-        const auto res = pearson(v1, v2, weights);
-        return std::make_pair(res.pcc, res.pvalue);
-      }
+  if (is_spearman && is_weighted) {
+    const auto res = Spearman<>{}(data.v1(), data.v2(), data.weights());
+    return std::make_pair(res.rho, res.pvalue);
+  }
+
+  if (is_pearson) {
+    const auto res = Pearson<>{}(data.v1(), data.v2());
+    return std::make_pair(res.pcc, res.pvalue);
+  }
+
+  assert(is_spearman);
+
+  const auto res = Spearman<>{}(data.v1(), data.v2());
+  return std::make_pair(res.rho, res.pvalue);
+}
+
+class ExternalCorrelationRunner {
+  boost::process::ipstream _stdout{};
+  boost::process::ipstream _stderr{};
+  boost::process::opstream _stdin{};
+
+  boost::process::child _c{};
+  std::filesystem::path _script{};
+
+ public:
+  ExternalCorrelationRunner() = delete;
+  explicit ExternalCorrelationRunner(std::string_view method) {
+    _script = testdir() / boost::filesystem::unique_path().string();
+    std::ofstream f(_script);
+    f << external_cmds.at(method);
+    f.close();
+    std::filesystem::permissions(_script, std::filesystem::perms::owner_exec,
+                                 std::filesystem::perm_options::add);
+
+    const auto exec =
+        boost::process::search_path(absl::StartsWith(method, "weighted_") ? "Rscript" : "python3");
+
+    this->_c = boost::process::child(
+        exec.string(), _script.string(),
+        boost::process::std_in<this->_stdin, boost::process::std_out> this->_stdout,
+        boost::process::std_err > this->_stderr);
+
+    assert(this->_c.running());
+  }
 
-      if (method == "spearman"sv) {
-        const auto res = spearman(v1, v2);
-        return std::make_pair(res.rho, res.pvalue);
+  ~ExternalCorrelationRunner() noexcept {
+    std::filesystem::remove(_script);
+    try {
+      if (!this->_c) {
+        std::ostringstream buff;
+        buff << this->_stderr.rdbuf();
+        if (const auto s = buff.str(); !s.empty()) {
+          throw std::runtime_error(fmt::format(FMT_STRING("{}\n"), s));
+        }
       }
+    } catch (const std::exception& e) {
+      fmt::print(
+          stderr,
+          FMT_STRING("An exception was raised while calling ~ExternalCorrelationRunner():  {}"),
+          e.what());
+    } catch (...) {
+      fmt::print(
+          stderr,
+          FMT_STRING("An unknown error occurred while calling ~ExternalCorrelationRunner()\n"));
+    }
+  }
+
+  void wait(std::chrono::milliseconds duration = std::chrono::seconds(15)) {
+    this->_c.wait_for(duration);
+  }
+
+  void write_to_stdin(const std::string& msg) {
+    assert(!msg.empty() && msg.back() == '\n');
+    this->_stdin.write(msg.data(), static_cast<std::streamsize>(msg.size()));
+    this->_stdin.flush();
+  }
 
-      assert(method == "weighted_spearman"sv);
-      const auto res = spearman(v1, v2, weights);
-      return std::make_pair(res.rho, res.pvalue);
-    }();
+  void read_from_stdout(std::string& buff) { std::getline(this->_stdout, buff); }
 
-    CHECK(Approx(cfx) == cfx_reference);
-    CHECK(Approx(pv) == pv_reference);
+  [[nodiscard]] std::string read_from_stdout() {
+    std::string buff;
+    this->read_from_stdout(buff);
+    return buff;
+  }
+
+  void read_from_stderr(std::string& buff) { std::getline(this->_stderr, buff); }
+
+  [[nodiscard]] std::string read_from_stderr() {
+    std::string buff;
+    this->read_from_stderr(buff);
+    return buff;
+  }
+
+  void signal_end_of_input() { this->_stdin.pipe().close(); }
+};
+
+template <class N1, class N2 = double, usize iterations = 1>
+static void run_correlation_test(std::string_view method, usize vector_size, N1 min, N1 max,
+                                 u64 seed = std::random_device{}()) {
+  static_assert(std::is_arithmetic_v<N1>);
+  static_assert(std::is_arithmetic_v<N2>);
+
+  const auto weighted_corr = absl::StartsWith(method, "weighted_");
+  CorrelationBuff<N1, N2> buff(vector_size, weighted_corr);
+  std::string sbuff;
+  std::vector<std::string_view> tok_buff;
+
+  auto rand_eng = random::PRNG(seed);
+
+  ExternalCorrelationRunner runner(method);
+  try {
+    for (usize i = 0; i < iterations; ++i) {
+      buff.generate(rand_eng, min, max);
+      buff.serialize(sbuff);
+
+      runner.write_to_stdin(sbuff);
+      const auto [corr, pv] = compute_correlation(method, buff);
+      runner.read_from_stdout(sbuff);
+      tok_buff = absl::StrSplit(sbuff, '\t');
+      REQUIRE(tok_buff.size() == 2);
+
+      const auto expected_corr = utils::parse_numeric_or_throw<double>(tok_buff.front());
+      const auto expected_pv = utils::parse_numeric_or_throw<double>(tok_buff.back());
+
+      CHECK(Catch::Approx(expected_corr) == corr);
+      CHECK(Catch::Approx(expected_pv) == pv);
+    }
+    runner.signal_end_of_input();
+    runner.wait();
+  } catch (const std::exception& e) {
+    const auto stderr_ = runner.read_from_stderr();
+    if (!stderr_.empty()) {
+      throw std::runtime_error(fmt::format(FMT_STRING("{}:\n{}"), e.what(), stderr_));
+    }
+    throw;
   }
-  v1.clear();
-  v2.clear();
-  input_data_ready = true;
-  input_data_cv.notify_one();
-  t.join();
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
@@ -109,8 +305,8 @@ TEST_CASE("Corr. test: Pearson wo ties", "[correlation][pearson][short]") {
   const std::vector<u32> v1{17, 86, 60, 77, 47, 3, 70, 87, 88, 92};
   const std::vector<u32> v2{70, 29, 85, 61, 80, 34, 60, 31, 73, 66};
   const auto [pcc, pv] = Pearson<>{}(v1, v2);
-  CHECK(Approx(pcc) == -0.033621194725622014);
-  CHECK(Approx(pv) == 0.926536715854247);
+  CHECK(Catch::Approx(pcc) == -0.033621194725622014);
+  CHECK(Catch::Approx(pv) == 0.926536715854247);
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
@@ -119,8 +315,8 @@ TEST_CASE("Corr. test: Weighted Pearson wo ties", "[correlation][pearson][short]
   const std::vector<u32> v2{70, 29, 85, 61, 80, 34, 60, 31, 73, 66};
   const std::vector<double> w{0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0};
   const auto [pcc, pv] = Pearson<>{}(v1, v2, w);
-  CHECK(Approx(pcc) == 0.1892337717235999250409);
-  CHECK(Approx(pv) == -1.0);
+  CHECK(Catch::Approx(pcc) == 0.1892337717235999250409);
+  CHECK(Catch::Approx(pv) == -1.0);
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
@@ -128,8 +324,8 @@ TEST_CASE("Corr. test: Pearson w ties", "[correlation][pearson][short]") {
   std::vector<u32> v1{17, 86, 60, 77, 47, 3, 70, 47, 88, 92};
   std::vector<u32> v2{70, 29, 85, 61, 80, 34, 60, 31, 73, 66};
   const auto [pcc, pv] = Pearson<>{}(v1, v2);
-  CHECK(Approx(pcc) == 0.16426413174421572);
-  CHECK(Approx(pv) == 0.6502118872600098);
+  CHECK(Catch::Approx(pcc) == 0.16426413174421572);
+  CHECK(Catch::Approx(pv) == 0.6502118872600098);
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
@@ -138,60 +334,46 @@ TEST_CASE("Corr. test: Weighted Pearson w ties", "[correlation][pearson][short]"
   const std::vector<u32> v2{70, 29, 85, 61, 80, 34, 60, 31, 73, 66};
   const std::vector<double> w{0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0};
   const auto [pcc, pv] = Pearson<>{}(v1, v2, w);
-  CHECK(Approx(pcc) == 0.5009581087644285890548);
-  CHECK(Approx(pv) == -1.0);
+  CHECK(Catch::Approx(pcc) == 0.5009581087644285890548);
+  CHECK(Catch::Approx(pv) == -1.0);
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
 TEST_CASE("Corr. test: Pearson (SciPy)", "[correlation][pearson][short]") {
-  random::PRNG_t rand_eng{2427588200550938527ULL};
-  const auto v1 = generate_random_vect(rand_eng, 1'000, 0, 15'000);
-  const auto v2 = generate_random_vect(rand_eng, 1'000, 0, 15'000);
-  std::vector<double> dummy_weights{};
-  const auto [pcc, pv] = Pearson<>{}(v1, v2);
-  const auto [pcc_reference, pv_reference] = external_corr(v1, v2, dummy_weights, "pearson");
-  CHECK(Approx(pcc) == pcc_reference);
-  CHECK(Approx(pv) == pv_reference);
+  run_correlation_test("pearson", 1'000, 0, 15'000, 2427588200550938527ULL);
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
 TEST_CASE("Corr. test: Weighted Pearson (wCorr)", "[correlation][pearson][short]") {
-  random::PRNG_t rand_eng{17383284879759537016ULL};
-  const auto v1 = generate_random_vect(rand_eng, 1'000, 0, 15'000);
-  const auto v2 = generate_random_vect(rand_eng, 1'000, 0, 15'000);
-  const auto w = generate_random_vect(rand_eng, 1'000, 0.0, 1.0);
-  const auto [pcc, pv] = Pearson<>{}(v1, v2, w);
-  const auto [pcc_reference, pv_reference] = external_corr(v1, v2, w, "weighted_pearson");
-  CHECK(Approx(pcc) == pcc_reference);
-  CHECK(Approx(pv) == pv_reference);
+  run_correlation_test("weighted_pearson", 1'000, 0, 15'000, 17383284879759537016ULL);
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
 TEST_CASE("Corr. test: Pearson long (SciPy)", "[correlation][pearson][medium]") {
-  test_correlation_w_random_vector("pearson", 1'000, 250, 0U, 15'000U);
-  test_correlation_w_random_vector("pearson", 1'000, 250, -7'250, 7'250);
-  test_correlation_w_random_vector("pearson", 1'000, 250, -7'250.0, 7'250.0);
+  run_correlation_test<u32, double, 250>("pearson", 1'000, 0, 15'000);
+  run_correlation_test<i32, double, 250>("pearson", 1'000, -7'250, 7'250);
+  run_correlation_test<double, double, 250>("pearson", 1'000, -7'250.0, 7'250.0);
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
 TEST_CASE("Corr. test: Pearson long vect. (SciPy)", "[correlation][pearson][long]") {
-  test_correlation_w_random_vector("pearson", 500'000, 2, 0U, 15'000U);
-  test_correlation_w_random_vector("pearson", 500'000, 2, -7'250, 7'250);
-  test_correlation_w_random_vector("pearson", 500'000, 2, -7'250.0, 7'250.0);
+  run_correlation_test<u32, double, 5>("pearson", 100'000, 0U, 15'000U);
+  run_correlation_test<i32, double, 5>("pearson", 100'000, -7'250, 7'250);
+  run_correlation_test<double, double, 5>("pearson", 100'000, -7'250.0, 7'250.0);
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
 TEST_CASE("Corr. test: Weighted Pearson long (wCorr)", "[correlation][pearson][medium]") {
-  test_correlation_w_random_vector("weighted_pearson", 1'000, 250, 0U, 15'000U);
-  test_correlation_w_random_vector("weighted_pearson", 1'000, 250, -7'250, 7'250);
-  test_correlation_w_random_vector("weighted_pearson", 1'000, 250, -7'250.0, 7'250.0);
+  run_correlation_test<u32, double, 250>("weighted_pearson", 1'000, 0U, 15'000U);
+  run_correlation_test<i32, double, 250>("weighted_pearson", 1'000, -7'250, 7'250);
+  run_correlation_test<double, double, 250>("weighted_pearson", 1'000, -7'250.0, 7'250.0);
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
 TEST_CASE("Corr. test: Weighted Pearson long vect. (wCorr)", "[correlation][pearson][long]") {
-  test_correlation_w_random_vector("weighted_pearson", 500'000, 2, 0U, 15'000U);
-  test_correlation_w_random_vector("weighted_pearson", 500'000, 2, -7'250, 7'250);
-  test_correlation_w_random_vector("weighted_pearson", 500'000, 2, -7'250.0, 7'250.0);
+  run_correlation_test<u32, double, 5>("weighted_pearson", 100'000, 0U, 15'000U);
+  run_correlation_test<i32, double, 5>("weighted_pearson", 100'000, -7'250, 7'250);
+  run_correlation_test<double, double, 5>("weighted_pearson", 100'000, -7'250.0, 7'250.0);
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
@@ -199,8 +381,8 @@ TEST_CASE("Corr. test: Spearman wo ties", "[correlation][spearman][short]") {
   const std::vector<u32> v1{17, 86, 60, 77, 47, 3, 70, 87, 88, 92};
   const std::vector<u32> v2{70, 29, 85, 61, 80, 34, 60, 31, 73, 66};
   const auto [rho, pv] = Spearman<>{}(v1, v2);
-  CHECK(Approx(rho) == -0.16363636363636364);
-  CHECK(Approx(pv) == 0.6514773427962428);
+  CHECK(Catch::Approx(rho) == -0.16363636363636364);
+  CHECK(Catch::Approx(pv) == 0.6514773427962428);
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
@@ -208,60 +390,46 @@ TEST_CASE("Corr. test: Spearman w ties", "[correlation][spearman][short]") {
   const std::vector<u32> v1{17, 86, 60, 77, 47, 3, 70, 47, 88, 92};
   const std::vector<u32> v2{70, 29, 85, 61, 80, 34, 60, 31, 73, 66};
   const auto [rho, pv] = Spearman<>{}(v1, v2);
-  CHECK(Approx(rho) == 0.024316221747202587);
-  CHECK(Approx(pv) == 0.9468397049085097);
+  CHECK(Catch::Approx(rho) == 0.024316221747202587);
+  CHECK(Catch::Approx(pv) == 0.9468397049085097);
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
 TEST_CASE("Corr. test: Spearman (SciPy)", "[correlation][spearman][short]") {
-  random::PRNG_t rand_eng{3860329809333667103ULL};
-  const auto v1 = generate_random_vect(rand_eng, 1'000, 0, 15'000);
-  const auto v2 = generate_random_vect(rand_eng, 1'000, 0, 15'000);
-  std::vector<double> dummy_weights{};
-  const auto [rho, pv] = Spearman<>{}(v1, v2);
-  const auto [rho_reference, pv_reference] = external_corr(v1, v2, dummy_weights, "spearman");
-  CHECK(Approx(rho) == rho_reference);
-  CHECK(Approx(pv) == pv_reference);
+  run_correlation_test("spearman", 1'000, 0, 15'000, 3860329809333667103ULL);
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
 TEST_CASE("Corr. test: Weighted Spearman (wCorr)", "[correlation][spearman][short]") {
-  random::PRNG_t rand_eng{8469130800688738654ULL};
-  const auto v1 = generate_random_vect(rand_eng, 1'000, 0, 15'000);
-  const auto v2 = generate_random_vect(rand_eng, 1'000, 0, 15'000);
-  const auto w = generate_random_vect(rand_eng, 1'000, 0.0, 1.0);
-  const auto [rho, pv] = Spearman<>{}(v1, v2, w);
-  const auto [rho_reference, pv_reference] = external_corr(v1, v2, w, "weighted_spearman");
-  CHECK(Approx(rho) == rho_reference);
-  CHECK(Approx(pv) == pv_reference);
+  run_correlation_test("weighted_spearman", 1'000, 0, 15'000, 8469130800688738654ULL);
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
 TEST_CASE("Corr. test: Spearman long (SciPy)", "[correlation][spearman][long]") {
-  test_correlation_w_random_vector("spearman", 1'000, 250, 0U, 15'000U);
-  test_correlation_w_random_vector("spearman", 1'000, 250, -7'250, 7'250);
-  test_correlation_w_random_vector("spearman", 1'000, 250, -7'250.0, 7'250.0);
+  run_correlation_test<u32, double, 250>("spearman", 1'000, 0, 15'000);
+  run_correlation_test<i32, double, 250>("spearman", 1'000, -7'250, 7'250);
+  run_correlation_test<double, double, 250>("spearman", 1'000, -7'250.0, 7'250.0);
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
 TEST_CASE("Corr. test: Spearman long vect. (SciPy)", "[correlation][spearman][long]") {
-  test_correlation_w_random_vector("spearman", 500'000, 2, 0U, 15'000U);
-  test_correlation_w_random_vector("spearman", 500'000, 2, -7'250, 7'250);
-  test_correlation_w_random_vector("spearman", 500'000, 2, -7'250.0, 7'250.0);
+  run_correlation_test<u32, double, 5>("spearman", 100'000, 0U, 15'000U);
+  run_correlation_test<i32, double, 5>("spearman", 100'000, -7'250, 7'250);
+  run_correlation_test<double, double, 5>("spearman", 100'000, -7'250.0, 7'250.0);
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
 TEST_CASE("Corr. test: Weighted Spearman long (wCorr)", "[correlation][spearman][medium]") {
-  test_correlation_w_random_vector("weighted_spearman", 1'000, 250, 0U, 15'000U);
-  test_correlation_w_random_vector("weighted_spearman", 1'000, 250, -7'250, 7'250);
-  test_correlation_w_random_vector("weighted_spearman", 1'000, 250, -7'250.0, 7'250.0);
+  run_correlation_test<u32, double, 250>("weighted_spearman", 1'000, 0U, 15'000U);
+  run_correlation_test<i32, double, 250>("weighted_spearman", 1'000, -7'250, 7'250);
+  run_correlation_test<double, double, 250>("weighted_spearman", 1'000, -7'250.0, 7'250.0);
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
 TEST_CASE("Corr. test: Weighted Spearman long vect. (wCorr)", "[correlation][spearman][long]") {
-  test_correlation_w_random_vector("weighted_spearman", 500'000, 2, 0U, 15'000U);
-  test_correlation_w_random_vector("weighted_spearman", 500'000, 2, -7'250, 7'250);
-  test_correlation_w_random_vector("weighted_spearman", 500'000, 2, -7'250.0, 7'250.0);
+  run_correlation_test<u32, double, 5>("weighted_spearman", 100'000, 0U, 15'000U);
+  run_correlation_test<i32, double, 5>("weighted_spearman", 100'000, -7'250, 7'250);
+  run_correlation_test<double, double, 5>("weighted_spearman", 100'000, -7'250.0, 7'250.0);
 }
 
 }  // namespace modle::test::stats
diff --git a/test/units/stats/correlation_utils_test.cpp b/test/units/stats/correlation_utils_test.cpp
index e546bcb8..ddeaaaec 100644
--- a/test/units/stats/correlation_utils_test.cpp
+++ b/test/units/stats/correlation_utils_test.cpp
@@ -1,8 +1,8 @@
 // Copyright (C) 2022 Roberto Rossini <roberros@uio.no>
 //
 // SPDX-License-Identifier: MIT
-#include <catch2/catch.hpp>  // for AssertionHandler, operator""_catch_sr, SourceLineInfo
-#include <vector>            // for vector
+#include <catch2/catch_test_macros.hpp>
+#include <vector>  // for vector
 
 #include "modle/common/common.hpp"  // for u32
 #include "modle/stats/correlation.hpp"
diff --git a/test/units/stats/descriptive_test.cpp b/test/units/stats/descriptive_test.cpp
index b63ad9ad..9f77b941 100644
--- a/test/units/stats/descriptive_test.cpp
+++ b/test/units/stats/descriptive_test.cpp
@@ -4,10 +4,11 @@
 
 #include "modle/stats/descriptive.hpp"
 
-#include <algorithm>         // for transform
-#include <catch2/catch.hpp>  // for Approx, operator==, AssertionHandler, operator""_catc...
-#include <string_view>       // for string_view_literals
-#include <vector>            // for vector
+#include <algorithm>  // for transform
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <string_view>  // for string_view_literals
+#include <vector>       // for vector
 
 #include "modle/common/common.hpp"  // for u8, usize
 #include "modle/common/utils.hpp"   // for identity::operator(), identity
@@ -27,7 +28,7 @@ TEST_CASE("Mean", "[stats][short]") {
   std::transform(v1.begin(), v1.end(), v2.begin(),
                  [](const auto n) { return FP{static_cast<float>(n)}; });
 
-  const auto result = Approx(5.0);
+  const auto result = Catch::Approx(5.0);
 
   CHECK(stats::mean(v1.begin(), v1.end()) == result);
   CHECK(stats::mean(v2.begin(), v2.end(), [](const auto& fp) { return fp.n; }) == result);
@@ -48,7 +49,7 @@ TEST_CASE("Moving average", "[stats][short]") {
           v1.size() - window_size);
 
   for (usize i = 0; i < results.size(); ++i) {
-    CHECK(results[i] == Approx(output[i]));
+    CHECK(results[i] == Catch::Approx(output[i]));
   }
 
   output.clear();
@@ -57,13 +58,13 @@ TEST_CASE("Moving average", "[stats][short]") {
                                 [](const auto& fp) { return fp.n; }) == v1.size() - window_size);
 
   for (usize i = 0; i < results.size(); ++i) {
-    CHECK(results[i] == Approx(output[i]));
+    CHECK(results[i] == Catch::Approx(output[i]));
   }
 
   output.clear();
   output.resize(1);
   REQUIRE(stats::moving_average(v1.begin(), v1.end(), output.begin(), v1.size() + 1) == 1);
-  CHECK(stats::mean(v1.begin(), v1.end()) == Approx(output.front()));
+  CHECK(stats::mean(v1.begin(), v1.end()) == Catch::Approx(output.front()));
 }
 
 // NOLINTNEXTLINE(readability-function-cognitive-complexity)
@@ -73,7 +74,7 @@ TEST_CASE("Sum of squared deviations", "[stats][short]") {
   std::transform(v1.begin(), v1.end(), v2.begin(),
                  [](const auto n) { return FP{static_cast<float>(n)}; });
 
-  const auto result = Approx(110.0);
+  const auto result = Catch::Approx(110.0);
 
   CHECK(stats::sum_of_squared_deviations(v1.begin(), v1.end()) == result);
   CHECK(stats::sum_of_squared_deviations(v2.begin(), v2.end(),
@@ -87,7 +88,7 @@ TEST_CASE("Variance", "[stats][short]") {
   std::transform(v1.begin(), v1.end(), v2.begin(),
                  [](const auto n) { return FP{static_cast<float>(n)}; });
 
-  const auto result = Approx(10.0);
+  const auto result = Catch::Approx(10.0);
 
   CHECK(stats::variance(v1.begin(), v1.end()) == result);
   CHECK(stats::variance(v2.begin(), v2.end(), [](const auto& fp) { return fp.n; }) == result);
@@ -100,7 +101,7 @@ TEST_CASE("Standard Deviation", "[stats][short]") {
   std::transform(v1.begin(), v1.end(), v2.begin(),
                  [](const auto n) { return FP{static_cast<float>(n)}; });
 
-  const auto result = Approx(3.1622776601683795);
+  const auto result = Catch::Approx(3.1622776601683795);
 
   CHECK(stats::standard_dev(v1.begin(), v1.end()) == result);
   CHECK(stats::standard_dev(v2.begin(), v2.end(), [](const auto& fp) { return fp.n; }) == result);
@@ -116,8 +117,8 @@ TEST_CASE("SED", "[stats][short]") {
   const std::vector<double> weights{0.97302005, 0.05226173, 0.15995629, 0.31495018, 0.95241483,
                                     0.87420081, 0.21360278, 0.0822476,  0.26402032, 0.49666325};
   // computed with scipy.spatial.distance.euclidean
-  const auto result = Approx(154.65977964758991);
-  const auto result_weighted = Approx(99.94033647108283);
+  const auto result = Catch::Approx(154.65977964758991);
+  const auto result_weighted = Catch::Approx(99.94033647108283);
 
   CHECK(stats::sed(v1.begin(), v1.end(), v2.begin()) == result);
   CHECK(stats::weighted_sed(v1.begin(), v1.end(), v2.begin(), weights.begin()) == result_weighted);
@@ -134,8 +135,8 @@ TEST_CASE("RMSE", "[stats][short]") {
                                     0.87420081, 0.21360278, 0.0822476,  0.26402032, 0.49666325};
 
   // computed with mean_squared_error from sklearn.metrics
-  const auto result = Approx(48.90771661061378);
-  const auto result_weighted = Approx(47.73515476405454);
+  const auto result = Catch::Approx(48.90771661061378);
+  const auto result_weighted = Catch::Approx(47.73515476405454);
 
   CHECK(stats::rmse(v1.begin(), v1.end(), v2.begin()) == result);
   CHECK(stats::weighted_rmse(v1.begin(), v1.end(), v2.begin(), weights.begin()) == result_weighted);
diff --git a/test/units/stats/misc_test.cpp b/test/units/stats/misc_test.cpp
index fbfa20fe..0ac93f14 100644
--- a/test/units/stats/misc_test.cpp
+++ b/test/units/stats/misc_test.cpp
@@ -8,7 +8,8 @@
 #include <fmt/format.h>
 
 #include <boost/process.hpp>
-#include <catch2/catch.hpp>
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
 #include <vector>
 
 #include "modle/common/numeric_utils.hpp"
@@ -60,7 +61,7 @@ TEST_CASE("Gaussian kernel (SciPy)", "[stats][short]") {
       const auto kernel = compute_gauss_kernel(size, sigma);
       REQUIRE(ref_kernel.size() == kernel.size());
       for (usize i = 0; i < kernel.size(); ++i) {
-        CHECK(Approx(ref_kernel[i]) == kernel[i]);
+        CHECK(Catch::Approx(ref_kernel[i]) == kernel[i]);
       }
     }
   }
diff --git a/test/units/stats/tests_test.cpp b/test/units/stats/tests_test.cpp
index 558511d6..c618076c 100644
--- a/test/units/stats/tests_test.cpp
+++ b/test/units/stats/tests_test.cpp
@@ -19,20 +19,21 @@
 #include <boost/process/pipe.hpp>                                 // for opstream, ipstream
 #include <boost/process/search_path.hpp>                          // for search_path
 #include <cassert>                                                // for assert
-#include <catch2/catch.hpp>                                       // for Approx, operator==, Ass...
-#include <cmath>                                                  // for trunc
-#include <condition_variable>                                     // for condition_variable
-#include <filesystem>                                             // for path
-#include <memory>                                                 // for make_shared
-#include <mutex>                                                  // for mutex, unique_lock
-#include <ostream>                                                // for basic_ostream::flush
-#include <random>                                                 // for random_device
-#include <stdexcept>                                              // for overflow_error
-#include <string>                                                 // for string, getline
-#include <string_view>                                            // for string_view, basic_stri...
-#include <thread>                                                 // for thread
-#include <utility>                                                // for make_pair
-#include <vector>                                                 // for vector
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <cmath>               // for trunc
+#include <condition_variable>  // for condition_variable
+#include <filesystem>          // for path
+#include <memory>              // for make_shared
+#include <mutex>               // for mutex, unique_lock
+#include <ostream>             // for basic_ostream::flush
+#include <random>              // for random_device
+#include <stdexcept>           // for overflow_error
+#include <string>              // for string, getline
+#include <string_view>         // for string_view, basic_stri...
+#include <thread>              // for thread
+#include <utility>             // for make_pair
+#include <vector>              // for vector
 
 #include "modle/common/common.hpp"  // for i64, usize
 #include "modle/common/const_map.hpp"
@@ -116,8 +117,8 @@ TEST_CASE("Binom test - two-sided", "[stats][short]") {
   const auto n1 = 143;
   const auto n2 = 220;
 
-  const auto result_two_sided_n1 = Approx(2.1004568301653535e-06);
-  const auto result_two_sided_n2 = Approx(0.20009222902827126);
+  const auto result_two_sided_n1 = Catch::Approx(2.1004568301653535e-06);
+  const auto result_two_sided_n2 = Catch::Approx(0.20009222902827126);
 
   CHECK(stats::binomial_test<TWO_SIDED>(k, n1) == result_two_sided_n1);
   CHECK(stats::binomial_test<TWO_SIDED>(k, n2) == result_two_sided_n2);
@@ -128,7 +129,7 @@ TEST_CASE("Binom test - less", "[stats][short]") {
   const auto k = 100;
   const auto n = 143;
 
-  const auto result_less = Approx(0.9999995628468261);
+  const auto result_less = Catch::Approx(0.9999995628468261);
 
   CHECK(stats::binomial_test<LESS>(k, n) == result_less);
 }
@@ -138,7 +139,7 @@ TEST_CASE("Binom test - greater", "[stats][short]") {
   const auto k = 100;
   const auto n = 143;
 
-  const auto result_greater = Approx(1.0502284150826767e-06);
+  const auto result_greater = Catch::Approx(1.0502284150826767e-06);
 
   CHECK(stats::binomial_test<GREATER>(k, n) == result_greater);
 }
@@ -180,7 +181,7 @@ TEST_CASE("Binom test - two-sided randomized (SciPy)", "[stats][long]") {
       output_data_ready = false;
     }
     const auto pv = stats::binomial_test<TWO_SIDED>(k, n);
-    CHECK(Approx(pv).margin(1.0e-250) == pv_py);
+    CHECK(Catch::Approx(pv).margin(1.0e-250) == pv_py);
   }
   k = 0;
   n = 0;