diff --git a/.gitignore b/.gitignore index c2c3e22..48b1d01 100644 --- a/.gitignore +++ b/.gitignore @@ -73,10 +73,12 @@ subprojects/libconfig/ subprojects/libconstants/ subprojects/liblogging/ subprojects/.wraplock +subprojects/tomlplusplus-*/ qhull.wrap quill.wrap yaml-cpp.wrap +tomlplusplus.wrap .vscode/ diff --git a/benchmarks/hashing/benchmark_composition_hash.cpp b/benchmarks/hashing/benchmark_composition_hash.cpp new file mode 100644 index 0000000..0eeb108 --- /dev/null +++ b/benchmarks/hashing/benchmark_composition_hash.cpp @@ -0,0 +1,139 @@ +#include "fourdst/composition/composition.h" +#include "fourdst/composition/utils/composition_hash.h" +#include "fourdst/composition/utils.h" +#include "fourdst/atomic/atomicSpecies.h" +#include "fourdst/atomic/species.h" + +#include +#include +#include +#include +#include +#include +#include + +template +void do_not_optimize(T&& datum) { + asm volatile("" : "+r" (datum)); +} + +uint32_t calc_num_bins(const std::vector& data) { + // Use Sturges' formula + const size_t n = data.size(); + return static_cast(std::ceil(std::log2(n) + 1)); +} + +std::string plot_ascii_histogram(std::vector data, std::string title) { + // Use std::format + const uint32_t nBins = calc_num_bins(data); + const double minVal = *std::ranges::min_element(data); + const double maxVal = *std::ranges::max_element(data); + + std::string histogram; + histogram += std::format("{:^60}\n", title); + histogram += std::string(60, '=') + "\n"; + std::vector bins(nBins, 0); + const double binWidth = (maxVal - minVal) / nBins; + for (const auto& value : data) { + const uint32_t binIndex = static_cast((value - minVal) / binWidth); + if (binIndex < nBins) { + bins[binIndex]++; + } else { + bins[nBins - 1]++; + } + } + const uint32_t maxBinCount = *std::ranges::max_element(bins); + for (uint32_t i = 0; i < nBins; ++i) { + const double binStart = minVal + i * binWidth; + const double binEnd = binStart + binWidth; + const uint32_t barLength = static_cast(std::round((static_cast(bins[i]) / maxBinCount) * 50.0)); + histogram += std::format("[{:.2e}, {:.2e}): {:>15} | {:}\n", + binStart, binEnd, bins[i], std::string(barLength, '*')); + } + return histogram; + +} + +std::chrono::duration build_and_hash_compositions(const size_t iter, const size_t nSpecies = 8) { + using namespace fourdst::composition; + using namespace fourdst::atomic; + + Composition comp; + size_t count = 0; + for (const auto& sp : species | std::views::values) { + if (count >= nSpecies) { + break; + } + comp.registerSpecies(sp); + comp.setMolarAbundance(sp, 0.1); + count++; + } + + const auto start = std::chrono::high_resolution_clock::now(); + for (size_t i = 0; i < iter; ++i) { + uint64_t hashValue = utils::CompositionHash::hash_exact(comp); + do_not_optimize(hashValue); + } + const auto end = std::chrono::high_resolution_clock::now(); + + const std::chrono::duration duration = (end - start)/iter; + return duration; +} + +int main() { + using namespace fourdst::composition; + using namespace fourdst::atomic; + + + const size_t nIterations = 1000; + std::vector durations; + durations.resize(nIterations); + for (size_t i = 0; i < nIterations; ++i) { + std::print("Iteration {}/{}\r", i + 1, nIterations); + auto duration = build_and_hash_compositions(1000, 100); + durations[i] = duration.count(); + } + std::println(""); + + std::println("Average time to build and hash composition over {} iterations: {} ns", nIterations, + std::accumulate(durations.begin(), durations.end(), 0.0) / nIterations); + std::println("Max time to build and hash composition over {} iterations: {} ns", nIterations, + *std::ranges::max_element(durations)); + std::println("Min time to build and hash composition over {} iterations: {} ns", nIterations, + *std::ranges::min_element(durations)); + std::println("Standard deviation of time to build and hash composition over {} iterations: {} ns", nIterations, + [] (const std::vector& data, const double mean) { + double sum = 0.0; + for (const auto& d : data) { + sum += (d - mean) * (d - mean); + } + return std::sqrt(sum / data.size()); + } (durations, std::accumulate(durations.begin(), durations.end(), 0.0) / nIterations) + ); + std::println("Index of max time: {}", std::distance(durations.begin(), + std::ranges::max_element(durations))); + std::println("Index of min time: {}", std::distance(durations.begin(), + std::ranges::min_element(durations))); + + std::vector log_duration = durations; + std::ranges::transform(log_duration, log_duration.begin(), [](const double d) { + return std::log10(d); + }); + + std::vector filtered_durations; + const double mean = std::accumulate(durations.begin(), durations.end(), 0.0) / durations.size(); + const double stddev = [] (const std::vector& data, const double mean) { + double sum = 0.0; + for (const auto& d : data) { + sum += (d - mean) * (d - mean); + } + return std::sqrt(sum / data.size()); + } (durations, mean); + + for (const auto& d : durations) { + if (std::abs(d - mean) <= 3 * stddev) { + filtered_durations.push_back(d); + } + } + std::println("{}", plot_ascii_histogram(filtered_durations, "Build and Hash Composition Times (ns)")); +} \ No newline at end of file diff --git a/benchmarks/hashing/meson.build b/benchmarks/hashing/meson.build new file mode 100644 index 0000000..d92b6d5 --- /dev/null +++ b/benchmarks/hashing/meson.build @@ -0,0 +1 @@ +executable('hashing_bench', 'benchmark_composition_hash.cpp', dependencies: [composition_dep]) \ No newline at end of file diff --git a/benchmarks/meson.build b/benchmarks/meson.build new file mode 100644 index 0000000..4b32e51 --- /dev/null +++ b/benchmarks/meson.build @@ -0,0 +1 @@ +subdir('hashing') \ No newline at end of file diff --git a/meson.build b/meson.build index f8cab4d..279f9b0 100644 --- a/meson.build +++ b/meson.build @@ -18,7 +18,7 @@ # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA # # *********************************************************************** # -project('libcomposition', 'cpp', version: 'v2.2.5', default_options: ['cpp_std=c++23'], meson_version: '>=1.5.0') +project('libcomposition', 'cpp', version: 'v2.2.6', default_options: ['cpp_std=c++23'], meson_version: '>=1.5.0') # Add default visibility for all C++ targets add_project_arguments('-fvisibility=default', language: 'cpp') @@ -39,6 +39,10 @@ if get_option('build_examples') subdir('examples') endif +if get_option('build_benchmarks') + subdir('benchmarks') +endif + if get_option('pkg_config') message('Generating pkg-config file for libcomposition...') pkg = import('pkgconfig') diff --git a/meson_options.txt b/meson_options.txt index cdd7739..07eb2ed 100644 --- a/meson_options.txt +++ b/meson_options.txt @@ -1,3 +1,4 @@ option('pkg_config', type: 'boolean', value: true, description: 'generate pkg-config file for libcomposition (fourdst_composition.pc)') option('build_tests', type: 'boolean', value: true, description: 'build unit tests (uses gtest)') option('build_examples', type: 'boolean', value: true, description: 'build example programs') +option('build_benchmarks', type: 'boolean', value: false, description: 'build benchmark programs') diff --git a/src/composition/include/fourdst/composition/utils/composition_hash.h b/src/composition/include/fourdst/composition/utils/composition_hash.h index f2d7e0a..a9b3542 100644 --- a/src/composition/include/fourdst/composition/utils/composition_hash.h +++ b/src/composition/include/fourdst/composition/utils/composition_hash.h @@ -6,141 +6,95 @@ #include #include "xxhash64.h" +#include "fourdst/composition/composition.h" +#include "fourdst/composition/composition_abstract.h" namespace fourdst::composition::utils { struct CompositionHash { - static constexpr uint64_t kSeed = 0xC04D5EEDBEEFull; - static constexpr char kTag[] = "4DSTAR:Composition"; - template static uint64_t hash_exact(const CompositionT& comp) { - std::vector buf; - reserve_bytes(comp, buf); - write_header(comp, buf); + uint64_t h0 = kSeed; + uint64_t h1 = kSeed ^ kPrime1; + uint64_t h2 = kSeed ^ kPrime2; + uint64_t h3 = kSeed ^ kPrime3; - for (auto it = comp.begin(); it != comp.end(); ++it) { - const auto& species = it->first; - const double abundance = it->second; + auto it = comp.begin(); + size_t remaining = comp.size(); - const std::uint32_t spWord = pack_species(species); - push_le32(buf, spWord); + while (remaining >= 4) { + const auto& p0 = *it; + ++it; + h0 ^= pack_species_id(p0.first); + h0 = mum(h0, kPrime1); + h0 ^= normalize_double_bits(p0.second); + h0 = mum(h0, kPrime2); - const std::uint64_t bits = normalize_double_bits(abundance); - push_le64(buf, bits); + const auto& p1 = *it; + ++it; + h1 ^= pack_species_id(p1.first); + h1 = mum(h1, kPrime1); + h1 ^= normalize_double_bits(p1.second); + h1 = mum(h1, kPrime2); + + const auto& p2 = *it; + ++it; + h2 ^= pack_species_id(p2.first); + h2 = mum(h2, kPrime1); + h2 ^= normalize_double_bits(p2.second); + h2 = mum(h2, kPrime2); + + const auto& p3 = *it; + ++it; + h3 ^= pack_species_id(p3.first); + h3 = mum(h3, kPrime1); + h3 ^= normalize_double_bits(p3.second); + h3 = mum(h3, kPrime2); + + remaining -= 4; } - return XXHash64::hash(buf.data(), buf.size(), kSeed); - } - - static inline bool is_finite(double v) noexcept { - return std::isfinite(v); - } - - static inline std::int64_t quantize_index(double v, double eps) noexcept { - const auto ld_v = static_cast(v); - const auto ld_eps = static_cast(eps); - - const long double scaled = ld_v / ld_eps; - const long long idx = std::llroundl(scaled); - return static_cast(idx); - } - - template - static uint64_t hash_quantized(const CompositionT& comp, double eps) noexcept { - std::vector buf; - reserve_bytes(comp, buf); - write_header(comp, buf); - push_bytes(buf, reinterpret_cast("quantized"), 9); - push_le64(buf, encode_fp64(eps)); - - for (auto it = comp.begin(); it != comp.end(); ++it) { - const auto& species = it->first; - const double abundance = it->second; - - const std::uint32_t spWord = pack_species(species); - push_le32(buf, spWord); - - if (!is_finite(abundance) || eps <= 0.0) { - const std::uint64_t bits = normalize_double_bits(abundance); - push_le64(buf, bits); - } else { - const std::int64_t idx = quantize_index(abundance, eps); - push_le64(buf, static_cast(idx)); - } + while (remaining > 0) { + const auto& p = *it; + ++it; + h0 ^= pack_species_id(p.first); + h0 = mum(h0, kPrime1); + h0 ^= normalize_double_bits(p.second); + h0 = mum(h0, kPrime2); + --remaining; } - return XXHash64::hash(buf.data(), buf.size(), kSeed ^ 0x7319'BEEF'1234ull); + return mum(h0 ^ h1 ^ h2 ^ h3, kPrime3); } - private: - template - static std::uint32_t pack_species(const SpeciesT& s) noexcept { - // Adjust accessors if your Species API differs. - const auto z = static_cast(s.z()); - const auto a = static_cast(s.a()); - return (static_cast(z) << 16) | static_cast(a); + static constexpr uint64_t kSeed = 0xC04D5EEDBEEFull; + static constexpr uint64_t kPrime1 = 0xa0761d6478bd642fULL; + static constexpr uint64_t kPrime2 = 0xe7037ed1a0b428dbULL; + static constexpr uint64_t kPrime3 = 0x8ebc6af09c88c6e3ULL; + + // --- Helper: Fast integer mixing --- + static inline uint64_t mum(const uint64_t a, const uint64_t b) noexcept { + const unsigned __int128 r = static_cast(a) * static_cast(b); + return static_cast(r) ^ static_cast(r >> 64); } - static inline std::uint64_t normalize_double_bits(double v) noexcept { + static inline uint64_t mix(const uint64_t h) noexcept { + return mum(h, kPrime1); + } + + // --- Normalization Logic --- + static inline uint64_t normalize_double_bits(double v) noexcept { if (v == 0.0) v = 0.0; // fold -0.0 -> +0.0 if (std::isnan(v)) { return 0x7ff8000000000000ULL; // canonical quiet NaN } - return std::bit_cast(v); + return std::bit_cast(v); } - static inline double quantize(double v, double eps) noexcept { - if (!std::isfinite(v) || eps <= 0.0) return v; - const double q = std::nearbyint(v / eps) * eps; - return (q == 0.0) ? 0.0 : q; - } - - static inline std::uint64_t encode_fp64(double v) noexcept { - return std::bit_cast(v); - } - - // ---------- byte helpers (explicit little-endian) ---------- - static inline void push_le32(std::vector& b, std::uint32_t x) { - b.push_back(static_cast( x & 0xFF)); - b.push_back(static_cast((x >> 8 ) & 0xFF)); - b.push_back(static_cast((x >> 16) & 0xFF)); - b.push_back(static_cast((x >> 24) & 0xFF)); - } - - static inline void push_le64(std::vector& b, std::uint64_t x) noexcept { - b.push_back(static_cast( x & 0xFF)); - b.push_back(static_cast((x >> 8 ) & 0xFF)); - b.push_back(static_cast((x >> 16) & 0xFF)); - b.push_back(static_cast((x >> 24) & 0xFF)); - b.push_back(static_cast((x >> 32) & 0xFF)); - b.push_back(static_cast((x >> 40) & 0xFF)); - b.push_back(static_cast((x >> 48) & 0xFF)); - b.push_back(static_cast((x >> 56) & 0xFF)); - } - - static inline void push_bytes(std::vector& b, const std::uint8_t* p, std::size_t n) noexcept{ - b.insert(b.end(), p, p + n); - } - - template - static void write_header(const CompositionT& comp, std::vector& buf) noexcept { - push_bytes(buf, reinterpret_cast(kTag), sizeof(kTag) - 1); - - const std::size_t nRegistered = comp.getRegisteredSpecies().size(); - std::size_t nMolar = 0; - for (auto it = comp.begin(); it != comp.end(); ++it) { ++nMolar; } - - push_le64(buf, static_cast(nRegistered)); - push_le64(buf, static_cast(nMolar)); - } - - template - static void reserve_bytes(const CompositionT& comp, std::vector& buf) noexcept { - std::size_t nMolar = 0; - for (auto it = comp.begin(); it != comp.end(); ++it) { ++nMolar; } - const std::size_t approx = (sizeof(kTag) - 1) + 16 + nMolar * (4 + 8 + 0 /*quantized flag optional*/); - buf.reserve(approx); + static inline uint32_t pack_species_id(const auto& s) noexcept { + const auto z = static_cast(s.z()); + const auto a = static_cast(s.a()); + return (static_cast(z) << 16) | static_cast(a); } }; } diff --git a/subprojects/libconfig.wrap b/subprojects/libconfig.wrap index 21f05b2..ba46330 100644 --- a/subprojects/libconfig.wrap +++ b/subprojects/libconfig.wrap @@ -1,6 +1,6 @@ [wrap-git] url = https://github.com/4D-STAR/libconfig.git -revision = v1.1.4 +revision = v2.0.2 depth = 1 [provide] diff --git a/tests/composition/compositionTest.cpp b/tests/composition/compositionTest.cpp index 70048ad..6e962b5 100644 --- a/tests/composition/compositionTest.cpp +++ b/tests/composition/compositionTest.cpp @@ -367,43 +367,6 @@ TEST_F(compositionTest, negativeZeroEqualsPositiveZero) { fourdst::composition::utils::CompositionHash::hash_exact(b)); } -TEST_F(compositionTest, quantizedIgnoresTinyJitter) { - fourdst::composition::Composition a, b; - a.registerSymbol("H-1"); b.registerSymbol("H-1"); - a.setMolarAbundance("H-1", 1.0); - b.setMolarAbundance("H-1", 1.0 + 5e-14); // smaller than eps - - const double eps = 1e-12; - const auto hqa = fourdst::composition::utils::CompositionHash::hash_quantized(a, eps); - const auto hqb = fourdst::composition::utils::CompositionHash::hash_quantized(b, eps); - ASSERT_EQ(hqa, hqb); - - // But exact should differ if the bit pattern changes - const auto hea = fourdst::composition::utils::CompositionHash::hash_exact(a); - const auto heb = fourdst::composition::utils::CompositionHash::hash_exact(b); - ASSERT_NE(hea, heb); -} - -TEST_F(compositionTest, quantizedDetectsAboveEps) { - fourdst::composition::Composition a, b; - a.registerSymbol("H-1"); b.registerSymbol("H-1"); - a.setMolarAbundance("H-1", 1.0); - b.setMolarAbundance("H-1", 1.0 + 2e-12); // larger than eps - - const double eps = 1e-12; - ASSERT_NE(fourdst::composition::utils::CompositionHash::hash_quantized(a, eps), - fourdst::composition::utils::CompositionHash::hash_quantized(b, eps)); -} - -TEST_F(compositionTest, exactVsQuantizedDifferentSeeds) { - fourdst::composition::Composition a; - a.registerSymbol("H-1"); a.setMolarAbundance("H-1", 0.5); - - const auto he = fourdst::composition::utils::CompositionHash::hash_exact(a); - const auto hq = fourdst::composition::utils::CompositionHash::hash_quantized(a, 1e-12); - ASSERT_NE(he, hq); -} - TEST_F(compositionTest, cloneAndCopyStable) { std::vector symbols = {"H-1", "He-4"}; std::vector abundances = {0.6, 0.4}; @@ -442,15 +405,12 @@ TEST_F(compositionTest, canonicalizeNaNIfAllowed) { } TEST_F(compositionTest, hash) { + using namespace fourdst::atomic; fourdst::composition::Composition a, b; a.registerSymbol("H-1"); b.registerSymbol("H-1"); a.setMolarAbundance("H-1", 0.6); b.setMolarAbundance("H-1", 0.6); EXPECT_NO_THROW((void)a.hash()); EXPECT_EQ(a.hash(), b.hash()); -} - -TEST_F(compositionTest, hashCaching) { - using namespace fourdst::atomic; const std::unordered_map abundances ={ {H_1, 0.702}, {He_4, 0.06}, @@ -458,19 +418,8 @@ TEST_F(compositionTest, hashCaching) { {N_14, 0.0005}, {O_16, 0.22}, }; - const fourdst::composition::Composition a(abundances); - - const auto first_hash_clock_start = std::chrono::high_resolution_clock::now(); - (void)a.hash(); - const auto first_hash_clock_end = std::chrono::high_resolution_clock::now(); - const auto first_hash_duration = std::chrono::duration_cast(first_hash_clock_end - first_hash_clock_start).count(); - - const auto second_hash_clock_start = std::chrono::high_resolution_clock::now(); - (void)a.hash(); - const auto second_hash_clock_end = std::chrono::high_resolution_clock::now(); - const auto second_hash_duration = std::chrono::duration_cast(second_hash_clock_end - second_hash_clock_start).count(); - - EXPECT_LT(second_hash_duration, first_hash_duration); + fourdst::composition::Composition c(abundances), d(abundances); + EXPECT_EQ(c.hash(), d.hash()); } TEST_F(compositionTest, hashStaleing) {