tboudreaux/libcomposition
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #9 from tboudreaux/perf/vec

Browse Source 
Cache Locality
This commit is contained in:
Emily Boudreaux
2025-12-08 11:52:50 -05:00
committed by GitHub
parent 184df676ca a49dc32fdb
commit e3b2a19de6
19 changed files with 911 additions and 477 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
Doxyfile
View File

@@ -48,7 +48,7 @@ PROJECT_NAME = fourdst::libcomposition
# could be handy for archiving the generated documentation or if some version
# control system is used.
PROJECT_NUMBER = v2.2.2
PROJECT_NUMBER = v2.3.0
# Using the PROJECT_BRIEF tag one can provide an optional one line description
# for a project that appears at the top of each page and should give viewers a

123
benchmarks/ConstructionAndIteration/benchmark_composition_construction_and_iteration.cpp Normal file
View File

@@ -0,0 +1,123 @@
#include "benchmark_utils.h"
#include "fourdst/composition/composition.h"
#include "fourdst/atomic/species.h"
#include <chrono>
#include <random>
#include <ranges>
std::chrono::duration<double, std::nano> benchmark_construction(const size_t iterations, const size_t nSpecies) {
using namespace fourdst::composition;
using namespace fourdst::atomic;
// Setup random machine to get random double between 0 and 1 for molar abundances
std::random_device rd;
std::mt19937 gen(rd());
std::uniform_real_distribution<> dis(0.0, 1.0);
std::vector<Species> species_to_register;
std::vector<double> molarAbundances;
size_t count = 0;
for (const auto& sp : species | std::views::values) {
if (count >= nSpecies) {
break;
}
species_to_register.push_back(sp);
molarAbundances.push_back(dis(gen));
count++;
}
const auto duration = fdst_benchmark_function([&]() {
for (size_t i = 0; i < iterations; ++i) {
fourdst::composition::Composition comp(species_to_register, molarAbundances);
}
});
return duration / static_cast<double>(iterations);
}
std::chrono::duration<double, std::nano> benchmark_access(const size_t iterations, const size_t nSpecies) {
using namespace fourdst::composition;
using namespace fourdst::atomic;
// Setup random machine to get random double between 0 and 1 for molar abundances
std::random_device rd;
std::mt19937 gen(rd());
std::uniform_real_distribution<> dis(0.0, 1.0);
std::vector<Species> species_to_register;
std::vector<double> molarAbundances;
size_t count = 0;
for (const auto& sp : species | std::views::values) {
if (count >= nSpecies) {
break;
}
species_to_register.push_back(sp);
molarAbundances.push_back(dis(gen));
count++;
}
const Composition comp(species_to_register, molarAbundances);
std::uniform_int_distribution<> sIDDis(0, nSpecies - 1);
std::vector<Species> random_lookup_species;
for (size_t i = 0; i < iterations; ++i) {
random_lookup_species.push_back(species_to_register[sIDDis(gen)]);
}
const auto duration = fdst_benchmark_function([&]() {
for (size_t i = 0; i < iterations; ++i) {
volatile double y = comp.getMolarAbundance(random_lookup_species[i]);
do_not_optimize(y);
}
});
return duration / static_cast<double>(iterations);
}
int main () {
constexpr size_t nIterations = 1000;
constexpr size_t nSpecies = 100;
std::vector<double> durations;
durations.resize(nIterations);
for (size_t i = 0; i < nIterations; ++i) {
std::print("Iteration {}/{}\r", i + 1, nIterations);
auto duration = benchmark_construction(10, nSpecies);
durations[i] = duration.count();
}
std::println("");
std::println("Average time to construct composition over {} iterations: {} ns", nIterations,
std::accumulate(durations.begin(), durations.end(), 0.0) / nIterations);
std::println("Max time to construct composition over {} iterations: {} ns", nIterations,
*std::ranges::max_element(durations));
std::println("Min time to construct composition over {} iterations: {} ns", nIterations,
*std::ranges::min_element(durations));
std::println("{}", plot_ascii_histogram(durations, "Composition Construction Time Histogram"));
durations.clear();
durations.resize(nIterations);
for (size_t i = 0; i < nIterations; ++i) {
std::print("Iteration {}/{}\r", i + 1, nIterations);
auto duration = benchmark_access(1000, nSpecies);
durations[i] = duration.count();
}
std::println("");
std::println("Average time to access composition over {} iterations: {} ns", nIterations,
std::accumulate(durations.begin(), durations.end(), 0.0) / nIterations);
std::println("Max time to access composition over {} iterations: {} ns", nIterations,
*std::ranges::max_element(durations));
std::println("Min time to access composition over {} iterations: {} ns", nIterations,
*std::ranges::min_element(durations));
std::println("{}", plot_ascii_histogram(durations, "Composition Access Time Histogram"));
}

1
benchmarks/ConstructionAndIteration/meson.build Normal file
View File

@@ -0,0 +1 @@
executable('construction_and_iteration_bench', 'benchmark_composition_construction_and_iteration.cpp', dependencies: [composition_dep], include_directories: [benchmark_utils_includes])

61
benchmarks/hashing/benchmark_composition_hash.cpp
View File

@@ -1,58 +1,16 @@
#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 <chrono>
#include <numeric>
#include <print>
#include <string>
#include <vector>
#include <cstdint>
#include <ranges>
#include <chrono>
template <class T>
void do_not_optimize(T&& datum) {
asm volatile("" : "+r" (datum));
}
#include "benchmark_utils.h"
uint32_t calc_num_bins(const std::vector<double>& data) {
// Use Sturges' formula
const size_t n = data.size();
return static_cast<uint32_t>(std::ceil(std::log2(n) + 1));
}
std::string plot_ascii_histogram(std::vector<double> 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<uint32_t> bins(nBins, 0);
const double binWidth = (maxVal - minVal) / nBins;
for (const auto& value : data) {
const uint32_t binIndex = static_cast<uint32_t>((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<uint32_t>(std::round((static_cast<double>(bins[i]) / maxBinCount) * 50.0));
histogram += std::format("[{:.2e}, {:.2e}): {:>15} | {:}\n",
binStart, binEnd, bins[i], std::string(barLength, '*'));
}
return histogram;
}
std::chrono::duration<double, std::nano> build_and_hash_compositions(const size_t iter, const size_t nSpecies = 8) {
using namespace fourdst::composition;
@@ -69,15 +27,14 @@ std::chrono::duration<double, std::nano> build_and_hash_compositions(const size_
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 auto duration = fdst_benchmark_function([&]() {
for (size_t i = 0; i < iter; ++i) {
uint64_t hashValue = utils::CompositionHash::hash_exact(comp);
do_not_optimize(hashValue);
}
});
const std::chrono::duration<double, std::nano> duration = (end - start)/iter;
return duration;
return duration / static_cast<double>(iter);
}
int main() {

2
benchmarks/hashing/meson.build
View File

@@ -1 +1 @@
executable('hashing_bench', 'benchmark_composition_hash.cpp', dependencies: [composition_dep])
executable('hashing_bench', 'benchmark_composition_hash.cpp', dependencies: [composition_dep], include_directories: [benchmark_utils_includes])

5
benchmarks/meson.build
View File

@@ -1 +1,4 @@
subdir('hashing')
benchmark_utils_includes = include_directories('utils')
subdir('hashing')
subdir('ConstructionAndIteration')

65
benchmarks/utils/benchmark_utils.h Normal file
View File

@@ -0,0 +1,65 @@
#pragma once
#include <algorithm>
#include <cstdint>
#include <vector>
#include <string>
#include <cmath>
#include <format>
#include <chrono>
template <class T>
void do_not_optimize(T&& datum) {
asm volatile("" : "+r" (datum));
}
inline uint32_t calc_num_bins(const std::vector<double>& data) {
const size_t n = data.size();
return static_cast<uint32_t>(std::ceil(std::log2(n) + 1));
}
inline std::string plot_ascii_histogram(std::vector<double> data, std::string title) {
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<uint32_t> bins(nBins, 0);
const double binWidth = (maxVal - minVal) / nBins;
for (const auto& value : data) {
const uint32_t binIndex = static_cast<uint32_t>((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<uint32_t>(std::round((static_cast<double>(bins[i]) / maxBinCount) * 50.0));
histogram += std::format("[{:.2e}, {:.2e}): {:>15} | {:}\n",
binStart, binEnd, bins[i], std::string(barLength, '*'));
}
return histogram;
}
template <typename Func>
auto fdst_benchmark_function(Func&& func_call) {
auto start = std::chrono::high_resolution_clock::now();
// Forward the callable
std::forward<Func>(func_call)();
auto end = std::chrono::high_resolution_clock::now();
auto duration = std::chrono::duration_cast<std::chrono::nanoseconds>(end - start);
do_not_optimize(duration.count());
return duration;
}

2
meson.build
View File

@@ -18,7 +18,7 @@
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
#
# *********************************************************************** #
project('libcomposition', 'cpp', version: 'v2.2.6', default_options: ['cpp_std=c++23'], meson_version: '>=1.5.0')
project('libcomposition', 'cpp', version: 'v2.3.0', 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')

42
src/composition/include/fourdst/atomic/atomicSpecies.h
View File

@@ -1,6 +1,7 @@
#pragma once
#include <format>
#include <string_view>
#include <string>
#include <optional>
@@ -275,8 +276,7 @@ namespace fourdst::atomic {
friend bool operator==(const Species& lhs, const Species& rhs);
friend bool operator!=(const Species& lhs, const Species& rhs);
friend bool operator<(const Species& lhs, const Species& rhs);
friend bool operator>(const Species& lhs, const Species& rhs);
friend std::partial_ordering operator<=>(const Species &lhs, const Species &rhs);
};
/**
* @brief Equality operator for Species. Compares based on name.
@@ -296,25 +296,18 @@ namespace fourdst::atomic {
inline bool operator!=(const Species& lhs, const Species& rhs) {
return (lhs.m_name != rhs.m_name);
}
/**
* @brief Less-than operator for Species. Compares based on atomic mass.
* @param lhs The left-hand side Species.
* @param rhs The right-hand side Species.
* @return `true` if lhs atomic mass is less than rhs atomic mass, `false` otherwise.
*/
inline bool operator<(const Species& lhs, const Species& rhs) {
return (lhs.m_atomicMass < rhs.m_atomicMass);
}
/**
* @brief Greater-than operator for Species. Compares based on atomic mass.
* @param lhs The left-hand side Species.
* @param rhs The right-hand side Species.
* @return `true` if lhs atomic mass is greater than rhs atomic mass, `false` otherwise.
*/
inline bool operator>(const Species& lhs, const Species& rhs) {
return (lhs.m_atomicMass > rhs.m_atomicMass);
inline std::partial_ordering operator<=>(const Species &lhs, const Species &rhs) {
if (const auto cmp = lhs.m_atomicMass <=> rhs.m_atomicMass; cmp != 0) {
return cmp;
}
return lhs.m_name <=> rhs.m_name;
}
/**
* @brief Converts a spin-parity string (JPI string) to a double-precision floating-point number representing the spin.
*
@@ -424,14 +417,15 @@ namespace fourdst::atomic {
* }
* @endcode
*/
template<>
struct std::hash<fourdst::atomic::Species> {
/**
* @brief Computes the hash for a Species object.
* @param s The Species object to hash.
* @return The hash value of the species' name.
*/
* @brief Computes the hash for a Species object.
* @param s The Species object to hash.
* @return The hash value of the species' name.
*/
size_t operator()(const fourdst::atomic::Species& s) const noexcept {
return std::hash<std::string>()(s.m_name);
}
}; // namespace std
};

50
src/composition/include/fourdst/composition/composition.h
View File

@@ -26,6 +26,7 @@
#include <optional>
#include <unordered_set>
#include <expected>
#include "fourdst/config/config.h"
#include "fourdst/logging/logging.h"
@@ -96,6 +97,9 @@ namespace fourdst::composition {
*/
// ReSharper disable once CppClassCanBeFinal
class Composition final : public CompositionAbstract {
public:
using iterator = detail::CompositionIterator<false>;
using const_iterator = detail::CompositionIterator<true>;
private:
/**
* @struct CompositionCache
@@ -138,6 +142,11 @@ namespace fourdst::composition {
!Ye.has_value() && !sortedSpecies.has_value() && !hash.has_value();
}
};
enum class SpeciesIndexLookupError : uint8_t {
NO_REGISTERED_SPECIES,
SPECIES_NOT_FOUND
};
private:
/**
* @brief Gets the logger instance for the Composition class. This is static to ensure that all composition
@@ -149,11 +158,18 @@ namespace fourdst::composition {
return logger;
}
std::set<atomic::Species> m_registeredSpecies; ///< Set of registered species in the composition.
std::map<atomic::Species, double> m_molarAbundances; ///< Map of species to their molar abundances.
// std::set<atomic::Species> m_registeredSpecies; ///< Set of registered species in the composition.
// std::map<atomic::Species, double> m_molarAbundances; ///< Map of species to their molar abundances.
std::vector<atomic::Species> m_species;
std::vector<double> m_molarAbundances;
mutable CompositionCache m_cache; ///< Cache for computed properties to avoid redundant calculations.
private:
[[nodiscard]] std::expected<std::ptrdiff_t, SpeciesIndexLookupError> findSpeciesIndex(const atomic::Species &species) const noexcept;
[[nodiscard]] static std::vector<atomic::Species> symbolVectorToSpeciesVector(const std::vector<std::string>& symbols);
public:
/**
* @brief Default constructor.
@@ -295,6 +311,7 @@ namespace fourdst::composition {
* @return A reference to this Composition.
*/
Composition& operator=(Composition const& other);
Composition& operator=(const CompositionAbstract& other);
/**
* @brief Registers a new symbol for inclusion in the composition.
@@ -546,7 +563,7 @@ namespace fourdst::composition {
* value of this method will only be valid as long as the Composition object is valid (i.e. it cannot
* outlive the Composition object it was called on).
*/
[[nodiscard]] const std::set<atomic::Species> &getRegisteredSpecies() const noexcept override;
[[nodiscard]] const std::vector<atomic::Species> &getRegisteredSpecies() const noexcept override;
/**
* @brief Gets the mass fractions of all species in the composition.
@@ -767,8 +784,8 @@ namespace fourdst::composition {
* for species are defined based on their atomic mass. When iterating over the molar abundance map, species will be
* seen in order from lightest to heaviest.
*/
[[nodiscard]] std::map<atomic::Species, double>::iterator begin() override {
return m_molarAbundances.begin();
[[nodiscard]] iterator begin() override {
return {m_species.begin(), m_molarAbundances.begin()};
}
/**
@@ -788,8 +805,8 @@ namespace fourdst::composition {
* for species are defined based on their atomic mass. When iterating over the molar abundance map, species will be
* seen in order from lightest to heaviest.
*/
[[nodiscard]] std::map<atomic::Species, double>::const_iterator begin() const override {
return m_molarAbundances.cbegin();
[[nodiscard]] const_iterator begin() const override {
return {m_species.cbegin(), m_molarAbundances.cbegin()};
}
/**
@@ -809,8 +826,8 @@ namespace fourdst::composition {
* for species are defined based on their atomic mass. When iterating over the molar abundance map, species will be
* seen in order from lightest to heaviest.
*/
[[nodiscard]] std::map<atomic::Species, double>::iterator end() override {
return m_molarAbundances.end();
[[nodiscard]] detail::CompositionIterator<false> end() override {
return {m_species.end(), m_molarAbundances.end()};
}
/**
@@ -830,8 +847,8 @@ namespace fourdst::composition {
* for species are defined based on their atomic mass. When iterating over the molar abundance map, species will be
* seen in order from lightest to heaviest.
*/
[[nodiscard]] std::map<atomic::Species, double>::const_iterator end() const override {
return m_molarAbundances.cend();
[[nodiscard]] detail::CompositionIterator<true> end() const override {
return {m_species.cend(), m_molarAbundances.cend()};
}
[[nodiscard]] std::size_t hash() const override;
@@ -841,18 +858,9 @@ namespace fourdst::composition {
inline bool operator==(const Composition& a, const Composition& b) noexcept {
if (a.size() != b.size()) return false;
// Compare species sets quickly
if (a.getRegisteredSpecies() != b.getRegisteredSpecies())
return false;
// Compare all abundances
for (auto itA = a.begin(), itB = b.begin();
itA != a.end() && itB != b.end(); ++itA, ++itB) {
if (itA->first != itB->first)
return false;
if (itA->second != itB->second)
return false;
}
return true;
return a.hash() == b.hash();
}
}; // namespace fourdst::composition

14
src/composition/include/fourdst/composition/composition_abstract.h
View File

@@ -1,10 +1,10 @@
#pragma once
#include "fourdst/atomic/atomicSpecies.h"
#include "fourdst/composition/iterators/composition_abstract_iterator.h"
#include <string>
#include <unordered_map>
#include <map>
#include <set>
#include <vector>
#include <memory>
@@ -35,6 +35,8 @@ namespace fourdst::composition {
*/
class CompositionAbstract {
public:
using iterator = detail::CompositionIterator<false>;
using const_iterator = detail::CompositionIterator<true>;
/**
* @brief Virtual destructor.
*/
@@ -66,7 +68,7 @@ namespace fourdst::composition {
* @brief Get all registered atomic species in the composition.
* @return A set of registered atomic species.
*/
[[nodiscard]] virtual const std::set<fourdst::atomic::Species> &getRegisteredSpecies() const noexcept = 0;
[[nodiscard]] virtual const std::vector<atomic::Species> &getRegisteredSpecies() const noexcept = 0;
/**
* @brief Get the mass fraction for all registered symbols.
@@ -175,10 +177,10 @@ namespace fourdst::composition {
[[nodiscard]] virtual std::unique_ptr<CompositionAbstract> clone() const = 0;
[[nodiscard]] virtual std::map<atomic::Species, double>::iterator begin() = 0;
[[nodiscard]] virtual std::map<atomic::Species, double>::iterator end() = 0;
[[nodiscard]] virtual std::map<atomic::Species, double>::const_iterator begin() const = 0;
[[nodiscard]] virtual std::map<atomic::Species, double>::const_iterator end() const = 0;
[[nodiscard]] virtual iterator begin() = 0;
[[nodiscard]] virtual iterator end() = 0;
[[nodiscard]] virtual const_iterator begin() const = 0;
[[nodiscard]] virtual const_iterator end() const = 0;
[[nodiscard]] virtual std::size_t hash() const = 0;
};

12
src/composition/include/fourdst/composition/decorators/composition_decorator_abstract.h
View File

@@ -1,6 +1,7 @@
#pragma once
#include "fourdst/atomic/atomicSpecies.h"
#include "fourdst/composition/iterators/composition_abstract_iterator.h"
#include "fourdst/composition/composition_abstract.h"
@@ -19,7 +20,7 @@ namespace fourdst::composition {
[[nodiscard]] bool contains(const std::string& symbol) const override { return m_base_composition->contains(symbol); };
[[nodiscard]] size_t size() const noexcept override { return m_base_composition->size(); };
[[nodiscard]] std::set<std::string> getRegisteredSymbols() const noexcept override { return m_base_composition->getRegisteredSymbols(); };
[[nodiscard]] const std::set<atomic::Species> &getRegisteredSpecies() const noexcept override { return m_base_composition->getRegisteredSpecies(); };
[[nodiscard]] const std::vector<atomic::Species> &getRegisteredSpecies() const noexcept override { return m_base_composition->getRegisteredSpecies(); };
[[nodiscard]] std::unordered_map<atomic::Species, double> getMassFraction() const noexcept override { return m_base_composition->getMassFraction(); };
[[nodiscard]] std::unordered_map<atomic::Species, double> getNumberFraction() const noexcept override { return m_base_composition->getNumberFraction(); };
[[nodiscard]] double getMassFraction(const std::string& symbol) const override { return m_base_composition->getMassFraction(symbol); };
@@ -36,13 +37,12 @@ namespace fourdst::composition {
[[nodiscard]] size_t getSpeciesIndex(const std::string& symbol) const override { return m_base_composition->getSpeciesIndex(symbol); };
[[nodiscard]] size_t getSpeciesIndex(const atomic::Species& species) const override { return m_base_composition->getSpeciesIndex(species); };
[[nodiscard]] atomic::Species getSpeciesAtIndex(const size_t index) const override { return m_base_composition->getSpeciesAtIndex(index); }
[[nodiscard]] size_t hash() const override { return m_base_composition->hash(); };
[[nodiscard]] std::map<atomic::Species, double>::iterator begin() override { return m_base_composition->begin(); };
[[nodiscard]] std::map<atomic::Species, double>::iterator end() override { return m_base_composition->end(); };
[[nodiscard]] detail::CompositionIterator<false> begin() override { return m_base_composition->begin(); };
[[nodiscard]] detail::CompositionIterator<false> end() override { return m_base_composition->end(); };
[[nodiscard]] std::map<atomic::Species, double>::const_iterator begin() const override { return std::as_const(*m_base_composition).begin(); };
[[nodiscard]] std::map<atomic::Species, double>::const_iterator end() const override { return std::as_const(*m_base_composition).end(); };
[[nodiscard]] detail::CompositionIterator<true> begin() const override { return std::as_const(*m_base_composition).begin(); };
[[nodiscard]] detail::CompositionIterator<true> end() const override { return std::as_const(*m_base_composition).end(); };
protected:
std::unique_ptr<CompositionAbstract> m_base_composition;
};

36
src/composition/include/fourdst/composition/decorators/composition_masked.h
View File

@@ -1,29 +1,36 @@
#pragma once
#include <vector>
#include <set>
#include <unordered_map>
#include <memory>
#include <string>
#include "fourdst/composition/composition_abstract.h"
#include "fourdst/composition/decorators/composition_decorator_abstract.h"
#include "fourdst/composition/exceptions/exceptions_composition.h"
#include "fourdst/composition/iterators/composition_abstract_iterator.h"
#include "fourdst/atomic/atomicSpecies.h"
namespace fourdst::composition {
class MaskedComposition final : public CompositionDecorator {
public:
using iterator = detail::CompositionIterator<false>;
using const_iterator = detail::CompositionIterator<true>;
MaskedComposition(
const CompositionAbstract& baseComposition,
const std::set<atomic::Species>& activeSpecies
const std::vector<atomic::Species>& activeSpecies
);
[[nodiscard]] bool contains(const atomic::Species &species) const noexcept override;
[[nodiscard]] bool contains(const std::string &symbol) const override;
[[nodiscard]] const std::set<atomic::Species>& getRegisteredSpecies() const noexcept override;
[[nodiscard]] const std::vector<atomic::Species>& getRegisteredSpecies() const noexcept override;
[[nodiscard]] std::set<std::string> getRegisteredSymbols() const noexcept override;
[[nodiscard]] size_t size() const noexcept override;
[[nodiscard]] std::unordered_map<atomic::Species, double> getMassFraction() const noexcept override;
[[nodiscard]] std::unordered_map<atomic::Species, double> getNumberFraction() const noexcept override;
[[nodiscard]] double getMassFraction(const std::string &symbol) const override;
@@ -37,29 +44,26 @@ namespace fourdst::composition {
[[nodiscard]] double getElectronAbundance() const noexcept override;
[[nodiscard]] std::vector<double> getMassFractionVector() const noexcept override;
[[nodiscard]] std::vector<double> getNumberFractionVector() const noexcept override;
[[nodiscard]] std::vector<double> getMolarAbundanceVector() const noexcept override;
[[nodiscard]] size_t getSpeciesIndex(const std::string &symbol) const override;
[[nodiscard]] size_t getSpeciesIndex(const atomic::Species &species) const override;
[[nodiscard]] atomic::Species getSpeciesAtIndex(size_t index) const override;
[[nodiscard]] std::unique_ptr<CompositionAbstract> clone() const override;
[[nodiscard]] std::map<atomic::Species, double>::iterator begin() override;
[[nodiscard]] iterator begin() override;
[[nodiscard]] iterator end() override;
[[nodiscard]] std::map<atomic::Species, double>::iterator end() override;
[[nodiscard]] const_iterator begin() const override;
[[nodiscard]] const_iterator end() const override;
[[nodiscard]] std::map<atomic::Species, double>::const_iterator begin() const override;
[[nodiscard]] size_t hash() const override;
[[nodiscard]] std::map<atomic::Species, double>::const_iterator end() const override;
private:
std::set<atomic::Species> m_activeSpecies;
std::map<atomic::Species, double> m_masked_composition;
std::vector<atomic::Species> m_activeSpecies;
std::vector<double> m_molarAbundances;
};
}

9
src/composition/include/fourdst/composition/exceptions/exceptions_composition.h
View File

@@ -2,6 +2,7 @@
#include <exception>
#include <string>
#include <utility>
namespace fourdst::composition::exceptions {
/**
@@ -22,14 +23,14 @@ namespace fourdst::composition::exceptions {
* @brief Constructs a CompositionError with an error message.
* @param message The error message.
*/
explicit CompositionError(const std::string& message)
explicit CompositionError(std::string message)
: m_message(std::move(message)) {}
/**
* @brief Returns the error message.
* @return A C-style string containing the error message.
*/
const char* what() const noexcept override{
[[nodiscard]] const char* what() const noexcept override{
return m_message.c_str();
}
};
@@ -60,10 +61,10 @@ namespace fourdst::composition::exceptions {
protected:
std::string m_message;
public:
explicit SpeciesError(const std::string& message)
explicit SpeciesError(std::string message)
: m_message(std::move(message)) {}
const char* what() const noexcept override {
[[nodiscard]] const char* what() const noexcept override {
return m_message.c_str();
}
};

106
src/composition/include/fourdst/composition/iterators/composition_abstract_iterator.h Normal file
View File

@@ -0,0 +1,106 @@
#pragma once
#include <vector>
#include <iterator>
#include <utility>
#include <compare>
#include "fourdst/atomic/atomicSpecies.h"
namespace fourdst::composition::detail {
template <bool IsConst>
class CompositionIterator {
public:
using iterator_category = std::random_access_iterator_tag;
using difference_type = std::ptrdiff_t;
// Returns a pair of references. Natively supports structured binding [sp, y]
using value_type = std::pair<const atomic::Species, double>;
// Define reference types based on const-ness
using SpeciesRef = const atomic::Species&;
using AbundRef = std::conditional_t<IsConst, const double&, double&>;
using reference = std::pair<SpeciesRef, AbundRef>;
struct ArrowProxy {
reference m_payload;
const reference* operator->() const { return &m_payload; }
};
using pointer = ArrowProxy;
private:
using SpecIt = std::vector<atomic::Species>::const_iterator;
using AbunIt = std::conditional_t<IsConst,
std::vector<double>::const_iterator,
std::vector<double>::iterator>;
SpecIt m_sIt;
AbunIt m_aIt;
public:
CompositionIterator() = default;
CompositionIterator(SpecIt sIt, AbunIt aIt) : m_sIt(sIt), m_aIt(aIt) {}
template <bool WasConst, typename = std::enable_if_t<IsConst && !WasConst>>
CompositionIterator(const CompositionIterator<WasConst>& other)
: m_sIt(other.getSpeciesIt()), m_aIt(other.getAbundanceIt()) {}
[[nodiscard]] SpecIt getSpeciesIt() const { return m_sIt; }
[[nodiscard]] AbunIt getAbundanceIt() const { return m_aIt; }
reference operator*() const {
return { *m_sIt, *m_aIt };
}
ArrowProxy operator->() const {
return ArrowProxy{ **this };
}
reference operator[](difference_type n) const {
return { *(m_sIt + n), *(m_aIt + n) };
}
// --- Movement ---
CompositionIterator& operator++() { ++m_sIt; ++m_aIt; return *this; }
CompositionIterator operator++(int) { auto tmp = *this; ++(*this); return tmp; }
CompositionIterator& operator--() { --m_sIt; --m_aIt; return *this; } // FIXED
CompositionIterator operator--(int) { auto tmp = *this; --(*this); return tmp; }
CompositionIterator& operator+=(difference_type n) { m_sIt += n; m_aIt += n; return *this; }
CompositionIterator& operator-=(difference_type n) { m_sIt -= n; m_aIt -= n; return *this; }
// --- Arithmetic ---
friend CompositionIterator operator+(CompositionIterator it, difference_type n) { return it += n; }
// Commutative addition (n + it)
friend CompositionIterator operator+(difference_type n, CompositionIterator it) { return it += n; }
friend CompositionIterator operator-(CompositionIterator it, difference_type n) { return it -= n; }
// Difference between iterators
friend difference_type operator-(const CompositionIterator& lhs, const CompositionIterator& rhs) {
return lhs.m_sIt - rhs.m_sIt;
}
template <bool R>
bool operator==(const CompositionIterator<R>& other) const { return m_sIt == other.getSpeciesIt(); }
template <bool R>
bool operator!=(const CompositionIterator<R>& other) const { return m_sIt != other.getSpeciesIt(); }
template <bool R>
bool operator<(const CompositionIterator<R>& other) const { return m_sIt < other.getSpeciesIt(); }
template <bool R>
bool operator>(const CompositionIterator<R>& other) const { return m_sIt > other.getSpeciesIt(); }
template <bool R>
bool operator<=(const CompositionIterator<R>& other) const { return m_sIt <= other.getSpeciesIt(); }
template <bool R>
bool operator>=(const CompositionIterator<R>& other) const { return m_sIt >= other.getSpeciesIt(); }
};
}

33
src/composition/include/fourdst/composition/utils/composition_hash.h
View File

@@ -1,7 +1,6 @@
#pragma once
#include <cstring>
#include <cmath>
#include <vector>
#include <bit>
@@ -10,8 +9,12 @@
#include "fourdst/composition/composition_abstract.h"
namespace fourdst::composition::utils {
// Make a concept which checks if a type inherits from CompositionAbstract
template <typename CompT>
concept CompositionType = std::is_base_of_v<CompositionAbstract, CompT>;
struct CompositionHash {
template <typename CompositionT>
template <CompositionType CompositionT>
static uint64_t hash_exact(const CompositionT& comp) {
uint64_t h0 = kSeed;
uint64_t h1 = kSeed ^ kPrime1;
@@ -99,13 +102,19 @@ namespace fourdst::composition::utils {
};
}
namespace std {
template<>
struct hash<fourdst::composition::Composition> {
std::size_t operator()(const fourdst::composition::Composition& c) const noexcept {
return static_cast<std::size_t>(
fourdst::composition::utils::CompositionHash::hash_exact(c)
);
}
};
}
template<>
struct std::hash<fourdst::composition::CompositionAbstract> {
std::size_t operator()(const fourdst::composition::CompositionAbstract& c) const noexcept {
return static_cast<std::size_t>(
fourdst::composition::utils::CompositionHash::hash_exact(c)
);
}
};
template<>
struct std::hash<fourdst::composition::Composition> {
std::size_t operator()(const fourdst::composition::Composition& c) const noexcept {
return static_cast<std::size_t>(
fourdst::composition::utils::CompositionHash::hash_exact(c)
);
}
};

694
src/composition/lib/composition.cpp
View File

@@ -34,36 +34,15 @@
#include "fourdst/atomic/atomicSpecies.h"
#include "fourdst/atomic/species.h"
#include "fourdst/composition/composition.h"
#include <numeric>
#include "fourdst/composition/utils/composition_hash.h"
#include "fourdst/composition/utils.h"
#include "fourdst/composition/exceptions/exceptions_composition.h"
namespace {
template<typename A, typename B>
std::vector<A> sortVectorBy(
std::vector<A> toSort,
const std::vector<B>& by
) {
std::vector<std::size_t> indices(by.size());
for (size_t i = 0; i < indices.size(); i++) {
indices[i] = i;
}
std::ranges::sort(indices, [&](size_t a, size_t b) {
return by[a] < by[b];
});
std::vector<A> sorted;
sorted.reserve(indices.size());
for (const auto idx: indices) {
sorted.push_back(toSort[idx]);
}
return sorted;
}
void throw_unknown_symbol(quill::Logger* logger, const std::string& symbol) {
LOG_ERROR(logger, "Symbol {} is not a valid species symbol (not in the species database)", symbol);
throw fourdst::composition::exceptions::UnknownSymbolError("Symbol " + symbol + " is not a valid species symbol (not in the species database)");
@@ -76,127 +55,142 @@ namespace {
}
namespace fourdst::composition {
Composition::Composition(
const std::vector<std::string>& symbols
) {
for (const auto& symbol : symbols) {
registerSymbol(symbol);
}
}
/////////////////////////////////////////////
/// Constructors without molar abundances ///
/// These all delegate to the ctor ///
/// vector<Species> ///
/////////////////////////////////////////////
Composition::Composition(
const std::set<std::string>& symbols
) {
for (const auto& symbol : symbols) {
registerSymbol(symbol);
}
}
) : Composition(symbols | std::ranges::to<std::vector>()) {}
Composition::Composition(
const std::set<atomic::Species> &species
) : Composition(species | std::ranges::to<std::vector>()) {}
Composition::Composition(
const std::unordered_set<std::string> &symbols
) : Composition(symbols | std::ranges::to<std::vector>()) {}
Composition::Composition(
const std::unordered_set<atomic::Species> &species
) : Composition(species | std::ranges::to<std::vector>()) {}
Composition::Composition(
const std::vector<std::string>& symbols
) : Composition(symbolVectorToSpeciesVector(symbols)) {}
Composition::Composition(
const std::vector<atomic::Species> &species
) {
for (const auto& s : species) {
registerSpecies(s);
}
m_species = species;
std::ranges::sort(m_species, [&](const atomic::Species& a, const atomic::Species& b) {
return a < b;
});
const auto last = std::ranges::unique(m_species).begin();
m_species.erase(last, m_species.end());
m_molarAbundances.resize(m_species.size(), 0.0);
}
Composition::Composition(
const std::set<atomic::Species> &species
) {
for (const auto& s : species) {
registerSpecies(s);
}
}
Composition::Composition(const std::unordered_set<std::string> &symbols) {
for (const auto& symbol : symbols) {
registerSymbol(symbol);
}
}
Composition::Composition(const std::unordered_set<atomic::Species> &species) {
for (const auto& s : species) {
registerSpecies(s);
}
}
//////////////////////////////////////////
/// Constructors with molar abundances ///
/// These all delegate to the ctor ///
/// vector<Species, vector<double>> ///
//////////////////////////////////////////
Composition::Composition(
const std::vector<std::string>& symbols,
const std::vector<double>& molarAbundances
) {
if (symbols.size() != molarAbundances.size()) {
LOG_CRITICAL(getLogger(), "The number of symbols and molarAbundances must be equal (got {} symbols and {} molarAbundances).", symbols.size(), molarAbundances.size());
throw exceptions::InvalidCompositionError("The number of symbols and fractions must be equal. Got " + std::to_string(symbols.size()) + " symbols and " + std::to_string(molarAbundances.size()) + " fractions.");
}
) : Composition(symbolVectorToSpeciesVector(symbols), molarAbundances) {}
for (const auto &[symbol, y] : std::views::zip(symbols, molarAbundances)) {
registerSymbol(symbol);
setMolarAbundance(symbol, y);
}
}
Composition::Composition(
const std::set<std::string> &symbols,
const std::vector<double> &molarAbundances
) : Composition(symbolVectorToSpeciesVector(symbols | std::ranges::to<std::vector>()), molarAbundances) {}
Composition::Composition(
const std::unordered_map<std::string, double> &symbolMolarAbundances
) : Composition(
symbolMolarAbundances | std::views::keys | std::ranges::to<std::vector>(),
symbolMolarAbundances | std::views::values | std::ranges::to<std::vector>()
) {}
Composition::Composition(
const std::map<std::string, double> &symbolMolarAbundances
) : Composition(
symbolMolarAbundances | std::views::keys | std::ranges::to<std::vector>(),
symbolMolarAbundances | std::views::values | std::ranges::to<std::vector>()
) {}
Composition::Composition(
const std::unordered_map<atomic::Species, double> &speciesMolarAbundances
) : Composition(
speciesMolarAbundances | std::views::keys | std::ranges::to<std::vector>(),
speciesMolarAbundances | std::views::values | std::ranges::to<std::vector>()
) {}
Composition::Composition(
const std::map<atomic::Species, double> &speciesMolarAbundances
) : Composition(
speciesMolarAbundances | std::views::keys | std::ranges::to<std::vector>(),
speciesMolarAbundances | std::views::values | std::ranges::to<std::vector>()
) {}
Composition::Composition(
const std::vector<atomic::Species> &species,
const std::vector<double> &molarAbundances
) {
if (species.size() != molarAbundances.size()) {
if (__builtin_expect(species.size() != molarAbundances.size(), 0)) {
LOG_CRITICAL(getLogger(), "The number of species and molarAbundances must be equal (got {} species and {} molarAbundances).", species.size(), molarAbundances.size());
throw exceptions::InvalidCompositionError("The number of species and fractions must be equal. Got " + std::to_string(species.size()) + " species and " + std::to_string(molarAbundances.size()) + " fractions.");
}
for (const auto& [s, y] : std::views::zip(species, molarAbundances)) {
registerSpecies(s);
setMolarAbundance(s, y);
const size_t numSpecies = species.size();
m_species.reserve(numSpecies);
m_molarAbundances.reserve(numSpecies);
for (size_t i = 0; i < numSpecies; ++i) {
m_species.push_back(species[i]);
if (__builtin_expect(molarAbundances[i] < 0.0, 0)) {
LOG_CRITICAL(getLogger(), "Molar abundance for species {} is negative (y = {}). Molar abundances must be non-negative.", species[i].name(), molarAbundances[i]);
throw exceptions::InvalidCompositionError("Molar abundance for species " + std::string(species[i].name()) + " is negative (y = " + std::to_string(molarAbundances[i]) + "). Molar abundances must be non-negative.");
}
m_molarAbundances.push_back(molarAbundances[i]);
}
auto combined = std::views::zip(m_species, m_molarAbundances);
std::ranges::sort(combined, [](const auto& a, const auto& b) -> bool {
const auto& spA = std::get<0>(a);
const auto& spB = std::get<0>(b);
if (spA != spB) {
return spA < spB;
}
return std::get<1>(a) > std::get<1>(b);
});
auto [first, last] = std::ranges::unique(combined, [](const auto& a, const auto& b) {
return std::get<0>(a) == std::get<0>(b);
});
const auto newEndIndex = std::distance(combined.begin(), first);
m_species.erase(m_species.begin() + newEndIndex, m_species.end());
m_molarAbundances.erase(m_molarAbundances.begin() + newEndIndex, m_molarAbundances.end());
}
Composition::Composition(
const std::set<std::string> &symbols,
const std::vector<double> &molarAbundances
) {
if (symbols.size() != molarAbundances.size()) {
LOG_CRITICAL(getLogger(), "The number of symbols and molarAbundances must be equal (got {} symbols and {} molarAbundances).", symbols.size(), molarAbundances.size());
throw exceptions::InvalidCompositionError("The number of symbols and fractions must be equal. Got " + std::to_string(symbols.size()) + " symbols and " + std::to_string(molarAbundances.size()) + " fractions.");
}
////////////////////////////////////////////
/// Copy and conversion constructors ///
////////////////////////////////////////////
for (const auto& [symbol, y] : std::views::zip(sortVectorBy<std::string>(std::vector<std::string>(symbols.begin(), symbols.end()), molarAbundances), molarAbundances)) {
registerSymbol(symbol);
setMolarAbundance(symbol, y);
}
}
Composition::Composition(const std::unordered_map<std::string, double> &symbolMolarAbundances) {
for (const auto& [symbol, y] : symbolMolarAbundances) {
registerSymbol(symbol);
setMolarAbundance(symbol, y);
}
}
Composition::Composition(const std::map<std::string, double> &symbolMolarAbundances) {
for (const auto& [symbol, y] : symbolMolarAbundances) {
registerSymbol(symbol);
setMolarAbundance(symbol, y);
}
}
Composition::Composition(const std::unordered_map<atomic::Species, double> &speciesMolarAbundances) {
for (const auto& [species, y] : speciesMolarAbundances) {
registerSpecies(species);
setMolarAbundance(species, y);
}
}
Composition::Composition(const std::map<atomic::Species, double> &speciesMolarAbundances) {
for (const auto& [species, y] : speciesMolarAbundances) {
registerSpecies(species);
setMolarAbundance(species, y);
}
}
Composition::Composition(
const Composition &composition
) {
m_registeredSpecies = composition.m_registeredSpecies;
Composition::Composition(const Composition &composition) {
m_species = composition.m_species;
m_molarAbundances = composition.m_molarAbundances;
}
@@ -211,12 +205,32 @@ namespace fourdst::composition {
const Composition &other
) {
if (this != &other) {
m_registeredSpecies = other.m_registeredSpecies;
m_species = other.m_species;
m_molarAbundances = other.m_molarAbundances;
}
m_cache.clear();
return *this;
}
Composition & Composition::operator=(const CompositionAbstract &other) {
m_species.clear();
m_molarAbundances.clear();
m_cache.clear();
for (const auto& species : other.getRegisteredSpecies()) {
registerSpecies(species);
setMolarAbundance(species, other.getMolarAbundance(species));
}
return *this;
}
std::unique_ptr<CompositionAbstract> Composition::clone() const {
return std::make_unique<Composition>(*this);
}
//------------------------------------------
// Registration methods
//------------------------------------------
void Composition::registerSymbol(
const std::string& symbol
) {
@@ -231,41 +245,190 @@ namespace fourdst::composition {
void Composition::registerSymbol(
const std::vector<std::string>& symbols
) {
for (const auto& symbol : symbols) {
registerSymbol(symbol);
}
registerSpecies(symbolVectorToSpeciesVector(symbols));
}
void Composition::registerSpecies(
const atomic::Species &species
) noexcept {
m_registeredSpecies.insert(species);
if (!m_molarAbundances.contains(species)) {
m_molarAbundances.emplace(species, 0.0);
if (const auto it = std::ranges::lower_bound(m_species, species); it == m_species.end() || *it != species) {
const auto index = std::distance(m_species.begin(), it);
m_species.insert(it, species);
m_molarAbundances.insert(m_molarAbundances.begin() + index, 0.0);
m_cache.clear();
}
}
void Composition::registerSpecies(
const std::vector<atomic::Species> &species
) noexcept {
for (const auto& s : species) {
registerSpecies(s);
// We do not simply call registerSpecies(species) here as that would have a complexity of O(n^2) due to constantly
// reinserting into the vector. Rather we build the vector once and then sort it
if (species.empty()) return;
const size_t total_size = m_species.size() + species.size();
m_species.reserve(total_size);
m_molarAbundances.reserve(total_size);
for (const auto& sp : species) {
m_species.push_back(sp);
m_molarAbundances.push_back(0.0);
}
auto combined = std::views::zip(m_species, m_molarAbundances);
std::ranges::sort(combined, [](const auto& a, const auto& b) {
const auto& speciesA = std::get<0>(a);
const auto& speciesB = std::get<0>(b);
if (speciesA != speciesB) {
return speciesA < speciesB;
}
return std::get<1>(a) > std::get<1>(b);
});
auto [first, last] = std::ranges::unique(combined, [](const auto& a, const auto& b) {
return std::get<0>(a) == std::get<0>(b);
});
const auto newEndIndex = std::distance(combined.begin(), first);
m_species.erase(m_species.begin() + newEndIndex, m_species.end());
m_molarAbundances.erase(m_molarAbundances.begin() + newEndIndex, m_molarAbundances.end());
m_cache.clear();
}
std::set<std::string> Composition::getRegisteredSymbols() const noexcept {
std::set<std::string> symbols;
for (const auto& species : m_registeredSpecies) {
for (const auto& species : m_species) {
symbols.insert(std::string(species.name()));
}
return symbols;
}
const std::set<atomic::Species> &Composition::getRegisteredSpecies() const noexcept {
return m_registeredSpecies;
const std::vector<atomic::Species> &Composition::getRegisteredSpecies() const noexcept {
return m_species;
}
//------------------------------------------
// Molar abundance setters
//------------------------------------------
void Composition::setMolarAbundance(
const std::string &symbol,
const double &molar_abundance
) {
const auto species = getSpecies(symbol);
if (__builtin_expect(!species, 0)) {
throw_unknown_symbol(getLogger(), symbol);
}
setMolarAbundance(species.value(), molar_abundance);
}
void Composition::setMolarAbundance(
const atomic::Species &species,
const double &molar_abundance
) {
if (__builtin_expect(molar_abundance < 0.0, 0)) {
LOG_ERROR(getLogger(), "Molar abundance must be non-negative for symbol {}. Currently it is {}.", species.name(), molar_abundance);
throw exceptions::InvalidCompositionError("Molar abundance must be non-negative, got " + std::to_string(molar_abundance) + " for symbol " + std::string(species.name()) + ".");
}
const std::expected<std::ptrdiff_t, SpeciesIndexLookupError> speciesIndexResult = findSpeciesIndex(species);
if (__builtin_expect(!speciesIndexResult, 0)) {
throw_unregistered_symbol(getLogger(), std::string(species.name()));
}
assert(speciesIndexResult.value() < m_molarAbundances.size());
m_molarAbundances[speciesIndexResult.value()] = molar_abundance;
m_cache.clear();
}
////----------------------------------------------
/// Methods which set multiple molar abundances
/// delegate to vector<Species>, vector<double>
///-----------------------------------------------
void Composition::setMolarAbundance(
const std::vector<std::string> &symbols,
const std::vector<double> &molar_abundances
) {
setMolarAbundance(symbolVectorToSpeciesVector(symbols), molar_abundances);
}
void Composition::setMolarAbundance(
const std::set<std::string> &symbols,
const std::vector<double> &molar_abundances
) {
setMolarAbundance(symbolVectorToSpeciesVector(symbols | std::ranges::to<std::vector>()), molar_abundances);
}
void Composition::setMolarAbundance(
const std::set<atomic::Species> &species,
const std::vector<double> &molar_abundances
) {
setMolarAbundance(species | std::ranges::to<std::vector>(), molar_abundances);
}
void Composition::setMolarAbundance(
const std::vector<atomic::Species> &species,
const std::vector<double> &molar_abundances
) {
if (__builtin_expect(species.size() != molar_abundances.size(), 0)) {
LOG_CRITICAL(getLogger(), "The number of species and molar_abundances must be equal (got {} species and {} molar_abundances).", species.size(), molar_abundances.size());
throw exceptions::InvalidCompositionError("The number of species and fractions must be equal. Got " + std::to_string(species.size()) + " species and " + std::to_string(molar_abundances.size()) + " fractions.");
}
if (species.empty()) return;
if (species.size() == m_species.size()) {
if (species == m_species) {
for (const auto& [sp, y] : std::views::zip(species, molar_abundances)) {
if (__builtin_expect(y < 0.0, 0)) {
LOG_ERROR(getLogger(), "Molar abundance must be non-negative. Instead got {} for species {}.", y, sp.name());
throw exceptions::InvalidCompositionError("Molar abundance must be non-negative. Instead got " + std::to_string(y) + " for species " + std::string(sp.name()) + ".");
}
}
m_molarAbundances = molar_abundances;
m_cache.clear();
return;
}
}
for (size_t i = 0; i < species.size(); ++i) {
const double y = molar_abundances[i];
const auto& sp = species[i];
if (__builtin_expect(y < 0.0, 0)) {
LOG_CRITICAL(getLogger(), "Molar abundance must be non-negative. Instead got {} for species {}.", y, sp.name());
throw exceptions::InvalidCompositionError("Molar abundance must be non-negative. Instead got " + std::to_string(y) + " for species " + std::string(sp.name()) + ".");
}
const std::expected<std::ptrdiff_t, SpeciesIndexLookupError> speciesIndexResult = findSpeciesIndex(sp);
if (__builtin_expect(!speciesIndexResult, 0)) {
throw_unregistered_symbol(getLogger(), std::string(sp.name()));
}
const std::ptrdiff_t speciesIndex = speciesIndexResult.value();
m_molarAbundances[speciesIndex] = y;
}
m_cache.clear();
}
//------------------------------------------
// Fraction and abundance getters
//------------------------------------------
double Composition::getMassFraction(const std::string& symbol) const {
const auto species = getSpecies(symbol);
if (!species) {
@@ -277,22 +440,26 @@ namespace fourdst::composition {
double Composition::getMassFraction(
const atomic::Species &species
) const {
if (!m_molarAbundances.contains(species)) {
const std::expected<std::ptrdiff_t, SpeciesIndexLookupError> speciesIndexResult = findSpeciesIndex(species);
if (!speciesIndexResult) {
throw_unregistered_symbol(getLogger(), std::string(species.name()));
}
std::map<atomic::Species, double> raw_mass;
double totalMass = 0;
for (const auto& [sp, y] : m_molarAbundances) {
double speciesMass = 0;
for (const auto& [sp, y] : *this) {
const double contrib = y * sp.mass();
totalMass += contrib;
raw_mass.emplace(sp, contrib);
if (sp == species) {
speciesMass = contrib;
}
}
return raw_mass.at(species) / totalMass;
return speciesMass / totalMass;
}
std::unordered_map<atomic::Species, double> Composition::getMassFraction() const noexcept {
std::unordered_map<atomic::Species, double> mass_fractions;
for (const auto &species: m_molarAbundances | std::views::keys) {
for (const auto &species: *this | std::views::keys) {
mass_fractions.emplace(species, getMassFraction(species));
}
return mass_fractions;
@@ -312,19 +479,23 @@ namespace fourdst::composition {
double Composition::getNumberFraction(
const atomic::Species &species
) const {
if (!m_molarAbundances.contains(species)) {
const std::expected<std::ptrdiff_t, SpeciesIndexLookupError> speciesIndexResult = findSpeciesIndex(species);
if (!speciesIndexResult) {
throw_unregistered_symbol(getLogger(), std::string(species.name()));
}
double total_moles_per_gram = 0.0;
for (const auto &y: m_molarAbundances | std::views::values) {
total_moles_per_gram += y;
}
return m_molarAbundances.at(species) / total_moles_per_gram;
const std::ptrdiff_t speciesIndex = speciesIndexResult.value();
const double total_moles_per_gram = std::accumulate(
m_molarAbundances.begin(),
m_molarAbundances.end(),
0.0
);
return m_molarAbundances[speciesIndex] / total_moles_per_gram;
}
std::unordered_map<atomic::Species, double> Composition::getNumberFraction() const noexcept {
std::unordered_map<atomic::Species, double> number_fractions;
for (const auto &species: m_molarAbundances | std::views::keys) {
for (const auto &species: m_species) {
number_fractions.emplace(species, getNumberFraction(species));
}
return number_fractions;
@@ -344,25 +515,34 @@ namespace fourdst::composition {
double Composition::getMolarAbundance(
const atomic::Species &species
) const {
if (!m_molarAbundances.contains(species)) {
const std::expected<std::ptrdiff_t, SpeciesIndexLookupError> speciesIndexResult = findSpeciesIndex(species);
if (!speciesIndexResult) {
throw_unregistered_symbol(getLogger(), std::string(species.name()));
}
return m_molarAbundances.at(species);
const std::ptrdiff_t speciesIndex = speciesIndexResult.value();
return m_molarAbundances[speciesIndex];
}
//------------------------------------------
// Derived property getters
//------------------------------------------
double Composition::getMeanParticleMass() const noexcept {
std::vector<double> X = getMassFractionVector();
double sum = 0.0;
for (const auto& [species, x] : std::views::zip(m_registeredSpecies, X)) {
sum += x/species.mass();
double totalMass = 0.0;
double totalMoles = 0.0;
for (size_t i = 0; i < m_species.size(); ++i) {
totalMoles += m_molarAbundances[i];
totalMass += m_molarAbundances[i] * m_species[i].mass();
}
return 1.0 / sum;
return totalMass / totalMoles;
}
double Composition::getElectronAbundance() const noexcept {
double Ye = 0.0;
for (const auto& [species, y] : m_molarAbundances) {
for (const auto& [species, y] : *this) {
Ye += species.z() * y;
}
return Ye;
@@ -377,8 +557,8 @@ namespace fourdst::composition {
return m_cache.canonicalComp.value(); // Short circuit if we have cached the canonical composition
}
CanonicalComposition canonicalComposition;
const std::set<Species> canonicalH = {H_1, H_2, H_3, H_4, H_5, H_6, H_7};
const std::set<Species> canonicalHe = {He_3, He_4, He_5, He_6, He_7, He_8, He_9, He_10};
static const std::unordered_set<Species> canonicalH = {H_1, H_2, H_3, H_4, H_5, H_6, H_7};
static const std::unordered_set<Species> canonicalHe = {He_3, He_4, He_5, He_6, He_7, He_8, He_9, He_10};
for (const auto& symbol : canonicalH) {
if (contains(symbol)) {
@@ -391,11 +571,8 @@ namespace fourdst::composition {
}
}
for (const auto& species : m_molarAbundances | std::views::keys) {
const bool isHIsotope = canonicalH.contains(species);
const bool isHeIsotope = canonicalHe.contains(species);
if (isHIsotope || isHeIsotope) {
for (const auto& species : m_species) {
if (canonicalH.contains(species) || canonicalHe.contains(species)) {
continue; // Skip canonical H and He symbols
}
@@ -412,25 +589,25 @@ namespace fourdst::composition {
return canonicalComposition;
}
//------------------------------------------
// Vector getters
//------------------------------------------
std::vector<double> Composition::getMassFractionVector() const noexcept {
if (m_cache.massFractions.has_value()) {
return m_cache.massFractions.value(); // Short circuit if we have cached the mass fractions
}
std::vector<double> massFractionVector;
std::vector<double> speciesMass;
massFractionVector.reserve(m_molarAbundances.size());
speciesMass.reserve(m_molarAbundances.size());
for (const auto &species: m_molarAbundances | std::views::keys) {
for (const auto &species: m_species) {
massFractionVector.push_back(getMassFraction(species));
speciesMass.push_back(species.mass());
}
std::vector<double> massFractions = sortVectorBy(massFractionVector, speciesMass);
m_cache.massFractions = massFractions; // Cache the result
return massFractions;
m_cache.massFractions = massFractionVector; // Cache the result
return massFractionVector;
}
@@ -440,43 +617,25 @@ namespace fourdst::composition {
}
std::vector<double> numberFractionVector;
std::vector<double> speciesMass;
numberFractionVector.reserve(m_molarAbundances.size());
speciesMass.reserve(m_molarAbundances.size());
for (const auto &species: m_molarAbundances | std::views::keys) {
for (const auto &species: m_species) {
numberFractionVector.push_back(getNumberFraction(species));
speciesMass.push_back(species.mass());
}
std::vector<double> numberFractions = sortVectorBy(numberFractionVector, speciesMass);
m_cache.numberFractions = numberFractions; // Cache the result
return numberFractions;
m_cache.numberFractions = numberFractionVector; // Cache the result
return numberFractionVector;
}
std::vector<double> Composition::getMolarAbundanceVector() const noexcept {
if (m_cache.molarAbundances.has_value()) {
return m_cache.molarAbundances.value(); // Short circuit if we have cached the molar abundances
}
std::vector<double> molarAbundanceVector;
std::vector<double> speciesMass;
molarAbundanceVector.reserve(m_molarAbundances.size());
speciesMass.reserve(m_molarAbundances.size());
for (const auto &[species, y]: m_molarAbundances) {
molarAbundanceVector.push_back(y);
speciesMass.push_back(species.mass());
}
std::vector<double> molarAbundances = sortVectorBy(molarAbundanceVector, speciesMass);
m_cache.molarAbundances = molarAbundances; // Cache the result
return molarAbundances;
return m_molarAbundances;
}
//------------------------------------------
// Species index getters and lookups
//------------------------------------------
size_t Composition::getSpeciesIndex(
const std::string &symbol
) const {
@@ -491,66 +650,36 @@ namespace fourdst::composition {
size_t Composition::getSpeciesIndex(
const atomic::Species &species
) const {
if (!m_registeredSpecies.contains(species)) {
LOG_ERROR(getLogger(), "Species {} is not in the composition.", species.name());
throw exceptions::UnregisteredSymbolError("Species " + std::string(species.name()) + " is not in the composition.");
}
if (m_cache.sortedSpecies.has_value()) {
return std::distance(
m_cache.sortedSpecies->begin(),
std::ranges::find(
m_cache.sortedSpecies.value().begin(),
m_cache.sortedSpecies.value().end(),
species
)
);
std::expected<std::ptrdiff_t, SpeciesIndexLookupError> speciesIndexResult = findSpeciesIndex(species);
if (!speciesIndexResult) {
switch (speciesIndexResult.error()) {
case SpeciesIndexLookupError::NO_REGISTERED_SPECIES:
[[fallthrough]];
case SpeciesIndexLookupError::SPECIES_NOT_FOUND:
throw_unregistered_symbol(getLogger(), std::string(species.name()));
default:
throw std::logic_error("Unhandled SpeciesIndexLookupError in Composition::getSpeciesIndex");
}
}
std::vector<atomic::Species> speciesVector;
std::vector<double> speciesMass;
speciesVector.reserve(m_molarAbundances.size());
speciesMass.reserve(m_molarAbundances.size());
for (const auto &s: m_registeredSpecies) {
speciesVector.emplace_back(s);
speciesMass.push_back(s.mass());
}
std::vector<atomic::Species> sortedSpecies = sortVectorBy(speciesVector, speciesMass);
m_cache.sortedSpecies = sortedSpecies;
return std::distance(sortedSpecies.begin(), std::ranges::find(sortedSpecies, species));
return static_cast<size_t>(speciesIndexResult.value());
}
atomic::Species Composition::getSpeciesAtIndex(
const size_t index
) const {
if (m_cache.sortedSpecies.has_value()) {
return m_cache.sortedSpecies.value().at(index);
if (index >= m_species.size()) {
LOG_ERROR(getLogger(), "Index {} is out of bounds for registered species (size {}).", index, m_species.size());
throw std::out_of_range("Index " + std::to_string(index) + " is out of bounds for registered species (size " + std::to_string(m_species.size()) + ").");
}
std::vector<atomic::Species> speciesVector;
std::vector<double> speciesMass;
speciesVector.reserve(m_molarAbundances.size());
speciesMass.reserve(m_molarAbundances.size());
for (const auto &species: m_registeredSpecies) {
speciesVector.emplace_back(species);
speciesMass.push_back(species.mass());
}
std::vector<atomic::Species> sortedSymbols = sortVectorBy(speciesVector, speciesMass);
if (index >= sortedSymbols.size()) {
LOG_ERROR(getLogger(), "Index {} is out of range for composition of size {}.", index, sortedSymbols.size());
throw std::out_of_range("Index " + std::to_string(index) + " is out of range for composition of size " + std::to_string(sortedSymbols.size()) + ".");
}
return sortedSymbols.at(index);
return m_species[index];
}
std::unique_ptr<CompositionAbstract> Composition::clone() const {
return std::make_unique<Composition>(*this);
}
//------------------------------------------
// Utility methods
//------------------------------------------
std::size_t Composition::hash() const {
if (m_cache.hash.has_value()) {
@@ -564,7 +693,7 @@ namespace fourdst::composition {
bool Composition::contains(
const atomic::Species &species
) const noexcept {
return m_registeredSpecies.contains(species);
return std::ranges::binary_search(m_species, species);
}
bool Composition::contains(
@@ -578,74 +707,41 @@ namespace fourdst::composition {
}
size_t Composition::size() const noexcept {
return m_registeredSpecies.size();
return m_species.size();
}
void Composition::setMolarAbundance(
const std::string &symbol,
const double &molar_abundance
) {
const auto species = getSpecies(symbol);
if (!species) {
throw_unknown_symbol(getLogger(), symbol);
std::expected<std::ptrdiff_t, Composition::SpeciesIndexLookupError> Composition::findSpeciesIndex(const atomic::Species &species) const noexcept {
if (m_species.empty()) return std::unexpected(SpeciesIndexLookupError::NO_REGISTERED_SPECIES);
const auto it = std::ranges::lower_bound(m_species, species);
if (it == m_species.end() || *it != species) {
return std::unexpected(SpeciesIndexLookupError::SPECIES_NOT_FOUND);
}
setMolarAbundance(species.value(), molar_abundance);
return std::distance(m_species.begin(), it);
}
void Composition::setMolarAbundance(
const atomic::Species &species,
const double &molar_abundance
) {
if (!m_registeredSpecies.contains(species)) {
throw_unregistered_symbol(getLogger(), std::string(species.name()));
}
if (molar_abundance < 0.0) {
LOG_ERROR(getLogger(), "Molar abundance must be non-negative for symbol {}. Currently it is {}.", species.name(), molar_abundance);
throw exceptions::InvalidCompositionError("Molar abundance must be non-negative, got " + std::to_string(molar_abundance) + " for symbol " + std::string(species.name()) + ".");
std::vector<atomic::Species> Composition::symbolVectorToSpeciesVector(const std::vector<std::string> &symbols) {
std::vector<atomic::Species> species;
species.reserve(symbols.size());
for (const auto& symbol : symbols) {
const auto speciesResult = getSpecies(symbol);
if (!speciesResult) {
throw_unknown_symbol(getLogger(), symbol);
}
species.push_back(speciesResult.value());
}
m_cache.clear();
m_molarAbundances.at(species) = molar_abundance;
return species;
}
void Composition::setMolarAbundance(
const std::vector<std::string> &symbols,
const std::vector<double> &molar_abundances
) {
for (const auto& [symbol, y] : std::views::zip(symbols, molar_abundances)) {
setMolarAbundance(symbol, y);
}
}
void Composition::setMolarAbundance(
const std::vector<atomic::Species> &species,
const std::vector<double> &molar_abundances
) {
for (const auto& [s, y] : std::views::zip(species, molar_abundances)) {
setMolarAbundance(s, y);
}
}
void Composition::setMolarAbundance(
const std::set<std::string> &symbols,
const std::vector<double> &molar_abundances
) {
for (const auto& [symbol, y] : std::views::zip(symbols, molar_abundances)) {
setMolarAbundance(symbol, y);
}
}
void Composition::setMolarAbundance(
const std::set<atomic::Species> &species,
const std::vector<double> &molar_abundances
) {
for (const auto& [s, y] : std::views::zip(species, molar_abundances)) {
setMolarAbundance(s, y);
}
}
/// OVERLOADS
//------------------------------------------
// Stream operator
//------------------------------------------
std::ostream& operator<<(
std::ostream& os,
@@ -653,7 +749,7 @@ namespace fourdst::composition {
) {
os << "Composition(Mass Fractions => [";
size_t count = 0;
for (const auto &species : composition.m_registeredSpecies) {
for (const auto &species : composition.m_species) {
os << species << ": " << composition.getMassFraction(species);
if (count < composition.size() - 1) {
os << ", ";

72
src/composition/lib/decorators/composition_masked.cpp
View File

@@ -1,29 +1,41 @@
#include "fourdst/composition/decorators/composition_masked.h"
#include "fourdst/composition/exceptions/exceptions_composition.h"
#include "fourdst/atomic/species.h"
#include <algorithm>
#include <memory>
#include <string>
#include <vector>
#include <set>
#include <unordered_map>
#include "fourdst/composition/utils/composition_hash.h"
namespace fourdst::composition {
MaskedComposition::MaskedComposition(
MaskedComposition::MaskedComposition(
const CompositionAbstract& baseComposition,
const std::set<atomic::Species>& activeSpecies
const std::vector<atomic::Species>& activeSpecies
) :
CompositionDecorator(baseComposition.clone()),
m_activeSpecies(activeSpecies) {
std::ranges::sort(m_activeSpecies, [](const auto &a, const auto &b) {
return a < b;
});
m_molarAbundances.reserve(m_activeSpecies.size());
for (const auto& species : m_activeSpecies) {
if (CompositionDecorator::contains(species)) {
m_masked_composition.emplace(species, CompositionDecorator::getMolarAbundance(species));
if (!CompositionDecorator::contains(species)) {
m_molarAbundances.push_back(0.0);
} else {
m_masked_composition.emplace(species, 0.0);
m_molarAbundances.push_back(CompositionDecorator::getMolarAbundance(species));
}
}
}
bool MaskedComposition::contains(const atomic::Species &species) const noexcept{
if (m_activeSpecies.contains(species)) {
return true;
}
return false;
return std::ranges::contains(m_activeSpecies, species);
}
bool MaskedComposition::contains(const std::string &symbol) const {
@@ -31,14 +43,10 @@ MaskedComposition::MaskedComposition(
throw exceptions::UnknownSymbolError("Cannot find species '" + symbol + "' in base composition");
}
const atomic::Species& species = atomic::species.at(symbol);
if (m_activeSpecies.contains(species)) {
return true;
}
return false;
return contains(species);
}
const std::set<atomic::Species>& MaskedComposition::getRegisteredSpecies() const noexcept {
const std::vector<atomic::Species>& MaskedComposition::getRegisteredSpecies() const noexcept {
return m_activeSpecies;
}
@@ -179,11 +187,19 @@ MaskedComposition::MaskedComposition(
if (!contains(symbol)) {
throw exceptions::UnregisteredSymbolError("Species '" + symbol + "' is not part of the active species in the MaskedComposition.");
}
return std::distance(m_activeSpecies.begin(), m_activeSpecies.find(atomic::species.at(symbol)));
return std::distance(
m_activeSpecies.begin(),
std::ranges::find_if(m_activeSpecies,
[&symbol](const atomic::Species& sp) {
return std::string(sp.name()) == symbol;
}));
}
size_t MaskedComposition::getSpeciesIndex(const atomic::Species &species) const {
return std::distance(m_activeSpecies.begin(), m_activeSpecies.find(species));
return std::distance(
m_activeSpecies.begin(),
std::ranges::find(m_activeSpecies, species)
);
}
atomic::Species MaskedComposition::getSpeciesAtIndex(const size_t index) const {
@@ -199,19 +215,23 @@ MaskedComposition::MaskedComposition(
return std::make_unique<MaskedComposition>(*m_base_composition, m_activeSpecies);
}
std::map<atomic::Species, double>::iterator MaskedComposition::begin() {
return m_masked_composition.begin();
MaskedComposition::iterator MaskedComposition::begin() {
return {m_activeSpecies.begin(), m_molarAbundances.begin()};
}
std::map<atomic::Species, double>::iterator MaskedComposition::end() {
return m_masked_composition.end();
MaskedComposition::iterator MaskedComposition::end() {
return {m_activeSpecies.end(), m_molarAbundances.end()};
}
std::map<atomic::Species, double>::const_iterator MaskedComposition::begin() const {
return m_masked_composition.cbegin();
MaskedComposition::const_iterator MaskedComposition::begin() const {
return {m_activeSpecies.cbegin(), m_molarAbundances.cbegin()};
}
std::map<atomic::Species, double>::const_iterator MaskedComposition::end() const {
return m_masked_composition.cend();
MaskedComposition::const_iterator MaskedComposition::end() const {
return {m_activeSpecies.cend(), m_molarAbundances.cend()};
}
size_t MaskedComposition::hash() const {
return utils::CompositionHash::hash_exact<MaskedComposition>(*this);
}
};

59
tests/composition/compositionTest.cpp
View File

@@ -3,6 +3,7 @@
#include <string>
#include <algorithm>
#include <chrono>
#include <ranges>
#include "fourdst/atomic/atomicSpecies.h"
#include "fourdst/atomic/species.h"
@@ -176,12 +177,9 @@ TEST_F(compositionTest, setGetComposition) {
TEST_F(compositionTest, getRegisteredSpecies) {
fourdst::composition::Composition comp;
comp.registerSpecies({fourdst::atomic::Be_7, fourdst::atomic::H_1, fourdst::atomic::He_4});
auto registeredSpecies = comp.getRegisteredSpecies();
EXPECT_TRUE(registeredSpecies.contains(fourdst::atomic::H_1));
EXPECT_TRUE(registeredSpecies.contains(fourdst::atomic::He_4));
EXPECT_FALSE(registeredSpecies.contains(fourdst::atomic::Li_6));
const auto it1 = registeredSpecies.begin();
EXPECT_EQ(*it1, fourdst::atomic::H_1);
EXPECT_TRUE(comp.contains(fourdst::atomic::H_1));
EXPECT_TRUE(comp.contains(fourdst::atomic::He_4));
EXPECT_FALSE(comp.contains(fourdst::atomic::Li_6));
}
/**
@@ -397,7 +395,7 @@ TEST_F(compositionTest, canonicalizeNaNIfAllowed) {
fourdst::composition::Composition a, b;
a.registerSymbol("H-1"); b.registerSymbol("H-1");
double qnan1 = std::numeric_limits<double>::quiet_NaN();
double qnan2 = std::bit_cast<double>(std::uint64_t{0x7ff80000'00000042ULL});
auto qnan2 = std::bit_cast<double>(std::uint64_t{0x7ff80000'00000042ULL});
a.setMolarAbundance("H-1", qnan1);
b.setMolarAbundance("H-1", qnan2);
ASSERT_EQ(fourdst::composition::utils::CompositionHash::hash_exact(a),
@@ -437,4 +435,51 @@ TEST_F(compositionTest, hashStaleing) {
a.setMolarAbundance("C-12", 0.002);
const std::size_t hash2 = a.hash();
EXPECT_NE(hash1, hash2);
}
TEST_F(compositionTest, speciesOrdering) {
using namespace fourdst::atomic;
const std::unordered_map<Species, double> abundances ={
{C_12, 0.001},
{H_1, 0.702},
{O_16, 0.22},
{N_14, 0.0005},
{He_4, 0.06},
{Ar_30, 0.0001},
{Li_6, 0.0002}
};
const fourdst::composition::Composition a(abundances);
static const std::vector<Species> speciesInOrder = {
H_1, He_4, Li_6, C_12, N_14, O_16, Ar_30
};
const std::vector<Species>& specieFromComposition = a.getRegisteredSpecies();
for (const auto& [spe, spo] : std::views::zip(speciesInOrder, specieFromComposition)) {
EXPECT_EQ(spe, spo);
}
}
TEST_F(compositionTest, iterationOrdering) {
using namespace fourdst::atomic;
const std::unordered_map<Species, double> abundances ={
{C_12, 0.001},
{H_1, 0.702},
{O_16, 0.22},
{N_14, 0.0005},
{He_4, 0.06},
{Ar_30, 0.0001},
{Li_6, 0.0002}
};
const fourdst::composition::Composition a(abundances);
static const std::vector<Species> speciesInOrder = {
H_1, He_4, Li_6, C_12, N_14, O_16, Ar_30
};
size_t count = 0;
for (const auto &sp: a | std::views::keys) {
EXPECT_EQ(sp, speciesInOrder[count]);
++count;
}
}