diff --git a/.github/workflows/c-cpp.yml b/.github/workflows/c-cpp.yml index 4f2b5ea..e29f065 100644 --- a/.github/workflows/c-cpp.yml +++ b/.github/workflows/c-cpp.yml @@ -10,10 +10,10 @@ jobs: steps: - name: Clone ParticleZoo repository - uses: actions/checkout@v4 + uses: actions/checkout@v5 - name: Setup Python - uses: actions/setup-python@v5 + uses: actions/setup-python@v6 with: python-version: '3.x' @@ -68,10 +68,10 @@ jobs: steps: - name: Clone ParticleZoo repository - uses: actions/checkout@v4 + uses: actions/checkout@v5 - name: Setup Python - uses: actions/setup-python@v5 + uses: actions/setup-python@v6 with: python-version: '3.x' @@ -126,19 +126,19 @@ jobs: steps: - name: Clone ParticleZoo repository - uses: actions/checkout@v4 + uses: actions/checkout@v5 - name: Setup MSVC uses: ilammy/msvc-dev-cmd@v1 - name: Setup Python - uses: actions/setup-python@v5 + uses: actions/setup-python@v6 with: python-version: '3.x' - name: Cache ROOT framework id: cache-root - uses: actions/cache@v4 + uses: actions/cache@v5 with: path: C:\root key: root-v6.36.04-win64-python311-vc17 diff --git a/.gitignore b/.gitignore index db9123b..6ba7965 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,8 @@ config.status build tests +*.md +!README.md .vscode .venv docs/latex @@ -23,3 +25,4 @@ python/.tox/ python/.venv/ python/venv/ python/ENV/ +.DS_Store \ No newline at end of file diff --git a/Doxyfile b/Doxyfile index a9b6016..beab071 100644 --- a/Doxyfile +++ b/Doxyfile @@ -48,7 +48,7 @@ PROJECT_NAME = "ParticleZoo" # could be handy for archiving the generated documentation or if some version # control system is used. -PROJECT_NUMBER = v1.1 +PROJECT_NUMBER = v1.1.1 # 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 viewer a diff --git a/ParticleZoo Test Summary Report v1.1.1.pdf b/ParticleZoo Test Summary Report v1.1.1.pdf new file mode 100644 index 0000000..3d479a8 Binary files /dev/null and b/ParticleZoo Test Summary Report v1.1.1.pdf differ diff --git a/docs/ParticleZoo Reference Manual v1.1.pdf b/docs/ParticleZoo Reference Manual v1.1.1.pdf similarity index 58% rename from docs/ParticleZoo Reference Manual v1.1.pdf rename to docs/ParticleZoo Reference Manual v1.1.1.pdf index e3da1ab..cbd2f6e 100644 Binary files a/docs/ParticleZoo Reference Manual v1.1.pdf and b/docs/ParticleZoo Reference Manual v1.1.1.pdf differ diff --git a/docs/ParticleZoo Test Summary Report v1.1.0.pdf b/docs/ParticleZoo Test Summary Report v1.1.0.pdf deleted file mode 100644 index 9e1a7ac..0000000 Binary files a/docs/ParticleZoo Test Summary Report v1.1.0.pdf and /dev/null differ diff --git a/docs/scripts/custom_doxygen.sty b/docs/scripts/custom_doxygen.sty index afc40e9..010384c 100644 --- a/docs/scripts/custom_doxygen.sty +++ b/docs/scripts/custom_doxygen.sty @@ -78,14 +78,14 @@ % Set text at the outside edge (LE = even-left, RO = odd-right) % Version number will be replaced by gendocs.sh - \fancyfoot[RE,LO]{\bfseries ParticleZoo Reference Manual v1.1}% + \fancyfoot[RE,LO]{\bfseries ParticleZoo Reference Manual v1.1.1}% \fancyfoot[LE,RO]{\bfseries \thepage}% % Also for 'plain' pages (chapter openings) % Version number will be replaced by gendocs.sh \fancypagestyle{plain}{% \fancyhf{}% - \fancyfoot[RE,LO]{\bfseries ParticleZoo Reference Manual v1.1}% + \fancyfoot[RE,LO]{\bfseries ParticleZoo Reference Manual v1.1.1}% \fancyfoot[LE,RO]{\bfseries \thepage}% }% } \ No newline at end of file diff --git a/docs/scripts/gendocs.sh b/docs/scripts/gendocs.sh index fa40f58..1aef653 100755 --- a/docs/scripts/gendocs.sh +++ b/docs/scripts/gendocs.sh @@ -38,7 +38,7 @@ VERSION=$(gawk ' } } END { - version = "v" major "." minor + version = "v" major "." minor "." patch if (caveat != "" && caveat != " ") version = version " " caveat print version }' include/particlezoo/utilities/version.h) diff --git a/ext b/ext index 89b9f88..4178402 160000 --- a/ext +++ b/ext @@ -1 +1 @@ -Subproject commit 89b9f8867887b8b7c2e7e6c4345520167cd3c319 +Subproject commit 417840280c79c78ad63ebfe77e8bd5498bf0a23c diff --git a/include/particlezoo/IAEA/IAEAHeader.h b/include/particlezoo/IAEA/IAEAHeader.h index 73f5fb8..7d0d2bb 100644 --- a/include/particlezoo/IAEA/IAEAHeader.h +++ b/include/particlezoo/IAEA/IAEAHeader.h @@ -731,13 +731,13 @@ namespace ParticleZoo::IAEAphspFile bool wIsStored_; bool weightIsStored_; - float constantX_; - float constantY_; - float constantZ_; - float constantU_; - float constantV_; - float constantW_; - float constantWeight_; + float constantX_{}; + float constantY_{}; + float constantZ_{}; + float constantU_{}; + float constantV_{}; + float constantW_{}; + float constantWeight_{}; std::vector extraFloatData_; std::vector extraLongData_; diff --git a/include/particlezoo/IAEA/IAEAphspFile.h b/include/particlezoo/IAEA/IAEAphspFile.h index c0409fd..58536f8 100644 --- a/include/particlezoo/IAEA/IAEAphspFile.h +++ b/include/particlezoo/IAEA/IAEAphspFile.h @@ -59,6 +59,12 @@ namespace ParticleZoo::IAEAphspFile */ std::uint64_t getNumberOfOriginalHistories() const override; + /** + * @brief Check if this format stores incremental history numbers per-particle. + * @return true if the header has an extra long of type INCREMENTAL_HISTORY_NUMBER + */ + bool hasNativeIncrementalHistoryCounters() const override; + /** * @brief Get the number of particles of a specific type * @param particleType Type of particle to count @@ -241,6 +247,7 @@ namespace ParticleZoo::IAEAphspFile inline const IAEAHeader & Reader::getHeader() const { return header_; } inline std::uint64_t Reader::getNumberOfParticles() const { return header_.getNumberOfParticles(); } inline std::uint64_t Reader::getNumberOfOriginalHistories() const { return header_.getOriginalHistories(); } + inline bool Reader::hasNativeIncrementalHistoryCounters() const { return header_.hasExtraLong(IAEAHeader::EXTRA_LONG_TYPE::INCREMENTAL_HISTORY_NUMBER); } inline std::uint64_t Reader::getNumberOfParticles(ParticleType particleType) const { return header_.getNumberOfParticles(particleType); } inline std::size_t Reader::getParticleRecordStartOffset() const { return 0; } inline std::size_t Reader::getParticleRecordLength() const { return header_.getRecordLength(); } diff --git a/include/particlezoo/PhaseSpaceFileReader.h b/include/particlezoo/PhaseSpaceFileReader.h index 282ac9e..8d3fffc 100644 --- a/include/particlezoo/PhaseSpaceFileReader.h +++ b/include/particlezoo/PhaseSpaceFileReader.h @@ -112,6 +112,25 @@ namespace ParticleZoo */ virtual std::uint64_t getNumberOfRepresentedHistories() const; + /** + * @brief Check if this format can provide the number of represented histories + * without scanning the file. + * + * @return true if getNumberOfRepresentedHistories() is cheap (e.g. stored in header) + */ + virtual bool hasNativeRepresentedHistoryCount() const; + + /** + * @brief Check if this format directly stores incremental history numbers per-particle. + * + * When true, particles returned by getNextParticle() carry the + * INCREMENTAL_HISTORY_NUMBER property with file-sourced values. Otherwise + * the incremental history numbers may be indirectly derived or otherwise determined. + * + * @return true if incremental history data is directly stored per-particle in this format + */ + virtual bool hasNativeIncrementalHistoryCounters() const; + /** * @brief Get the number of Monte Carlo histories that have been read so far. * If the end of the file has been reached, this will return the total number of original histories @@ -623,6 +642,9 @@ namespace ParticleZoo throw std::runtime_error("getNumberOfRepresentedHistories() is not supported for this file format."); } + inline bool PhaseSpaceFileReader::hasNativeRepresentedHistoryCount() const { return false; } + inline bool PhaseSpaceFileReader::hasNativeIncrementalHistoryCounters() const { return false; } + inline std::size_t PhaseSpaceFileReader::getParticleRecordLength() const { throw std::runtime_error("getParticleRecordLength() must be implemented for binary formatted file writers."); } diff --git a/include/particlezoo/ROOT/ROOTphsp.h b/include/particlezoo/ROOT/ROOTphsp.h index f9dbad7..1814e53 100644 --- a/include/particlezoo/ROOT/ROOTphsp.h +++ b/include/particlezoo/ROOT/ROOTphsp.h @@ -138,6 +138,12 @@ namespace ParticleZoo::ROOT { */ std::uint64_t getNumberOfOriginalHistories() const override; + /** + * @brief Check if this format stores incremental history numbers per-particle. + * @return true if the TTree has a history number branch + */ + bool hasNativeIncrementalHistoryCounters() const override; + /** * @brief Get format-specific CLI commands for ROOT configuration. * @return Vector of ROOT-specific CLI commands @@ -209,6 +215,7 @@ namespace ParticleZoo::ROOT { inline std::uint64_t Reader::getNumberOfParticles() const { return numberOfParticles_; } inline std::uint64_t Reader::getNumberOfOriginalHistories() const { return numberOfOriginalHistories_; } + inline bool Reader::hasNativeIncrementalHistoryCounters() const { return treeHasHistoryNumber_; } /** diff --git a/include/particlezoo/TOPAS/TOPASphspFile.h b/include/particlezoo/TOPAS/TOPASphspFile.h index 263e79b..5c121c9 100644 --- a/include/particlezoo/TOPAS/TOPASphspFile.h +++ b/include/particlezoo/TOPAS/TOPASphspFile.h @@ -49,6 +49,23 @@ namespace ParticleZoo::TOPASphspFile */ std::uint64_t getNumberOfOriginalHistories() const override; + /** + * @brief Get the number of represented histories in the phase space. + * + * Only available for ASCII and BINARY sub-formats (not LIMITED). + * For LIMITED, delegates to the base class which throws. + * + * @return Number of represented histories from the header + * @throws std::runtime_error if format is LIMITED + */ + std::uint64_t getNumberOfRepresentedHistories() const override; + + /** + * @brief Check if this format can provide represented history count cheaply. + * @return true for ASCII and BINARY, false for LIMITED + */ + bool hasNativeRepresentedHistoryCount() const override; + /** * @brief Get the TOPAS format type of this file * @return TOPASFormat enum indicating ASCII, BINARY, or LIMITED @@ -165,6 +182,15 @@ namespace ParticleZoo::TOPASphspFile // Inline implementations for the Reader class inline std::uint64_t Reader::getNumberOfParticles() const { return header_.getNumberOfParticles(); } inline std::uint64_t Reader::getNumberOfOriginalHistories() const { return header_.getNumberOfOriginalHistories(); } + inline std::uint64_t Reader::getNumberOfRepresentedHistories() const { + if (formatType_ == TOPASFormat::LIMITED) { + return PhaseSpaceFileReader::getNumberOfRepresentedHistories(); // throws + } + return header_.getNumberOfRepresentedHistories(); + } + inline bool Reader::hasNativeRepresentedHistoryCount() const { + return formatType_ != TOPASFormat::LIMITED; + } inline TOPASFormat Reader::getTOPASFormat() const { return formatType_; } inline const Header & Reader::getHeader() const { return header_; } inline void Reader::setDetailedReading(bool enable) { readFullDetails_ = enable; } diff --git a/include/particlezoo/parallel/HistoryBalancedParallelReader.h b/include/particlezoo/parallel/HistoryBalancedParallelReader.h index 35754d4..f3258c2 100644 --- a/include/particlezoo/parallel/HistoryBalancedParallelReader.h +++ b/include/particlezoo/parallel/HistoryBalancedParallelReader.h @@ -8,6 +8,8 @@ #include #include #include +#include +#include namespace ParticleZoo { @@ -131,6 +133,32 @@ namespace ParticleZoo { */ std::uint64_t getHistoriesRead(size_t threadIndex) const; + /** + * @brief Gets the total number of particles processed by a specific thread. + * + * Returns the cumulative count of particles that have been read by this thread. + * + * @param threadIndex The index of the thread (0 to numThreads-1) + * @return Total number of particles processed + * + * @throws std::out_of_range If threadIndex is invalid + */ + std::uint64_t getParticlesRead(size_t threadIndex) const; + + /** + * @brief Gets the total number of particles read across all threads. + * + * @return Total number of particles read + */ + std::uint64_t getTotalParticlesRead() const; + + /** + * @brief Gets the total number of original histories read across all threads. + * + * @return Total number of original histories read + */ + std::uint64_t getTotalHistoriesRead() const; + /** * @brief Gets the total number of particles in the phase space file. * @@ -159,6 +187,27 @@ namespace ParticleZoo { */ std::uint64_t getNumberOfRepresentedHistories() const { return numberOfRepresentedHistories_; } + /** + * @brief Gets the number of threads used by this reader. + * + * @return Number of parallel threads + */ + std::size_t getNumberOfThreads() const { return readers_.size(); } + + /** + * @brief Checks if the underlying phase space format provides native represented history count. + * + * @return True if native represented history count is available, false otherwise + */ + bool hasNativeRepresentedHistoryCount() const; + + /** + * @brief Checks if the underlying phase space format provides native incremental history counters. + * + * @return True if native incremental history counters are available, false otherwise + */ + bool hasNativeIncrementalHistoryCounters() const; + /** * @brief Closes all underlying phase space file readers. * @@ -168,20 +217,28 @@ namespace ParticleZoo { void close(); private: + struct ThreadStatistics { + mutable std::mutex mutex; // For thread-safe updates + std::atomic particlesRead = 0; + std::atomic totalHistoriesRead = 0; + // Per-thread-only fields (no cross-thread access needed) + std::uint64_t historiesRead = 0; + std::uint64_t emptyHistoryError = 0; + HasMoreParticlesResult hasMoreParticlesCache = NEEDS_CHECKING; + }; + std::vector> readers_; + std::vector> threadStats_; + std::vector startingHistorys_; + bool hasNativeRepresentedHistoryCount_; + bool hasNativeIncrementalHistoryCounters_; bool hasGapsBetweenHistories_; - std::uint64_t totalHistoriesToRead_; std::uint64_t numberOfOriginalHistories_; std::uint64_t numberOfParticlesInPhsp_; std::uint64_t numberOfRepresentedHistories_; std::uint64_t emptyHistoriesBetweenEachHistory_; std::uint64_t perHistoryErrorContribution_; - std::vector emptyHistoryErrors_; - std::vector startingHistorys_; - std::vector historiesRead_; - std::vector totalHistoriesRead_; - std::vector hasMoreParticlesCache_; }; } \ No newline at end of file diff --git a/include/particlezoo/parallel/ParticleBalancedParallelReader.h b/include/particlezoo/parallel/ParticleBalancedParallelReader.h index 339c567..999c9a5 100644 --- a/include/particlezoo/parallel/ParticleBalancedParallelReader.h +++ b/include/particlezoo/parallel/ParticleBalancedParallelReader.h @@ -7,6 +7,8 @@ #include #include #include +#include +#include namespace ParticleZoo { @@ -119,6 +121,42 @@ namespace ParticleZoo { */ std::uint64_t getParticlesRead(size_t threadIndex) const; + /** + * @brief Gets the number of original histories processed by a specific thread. + * + * Returns the cumulative count of original histories (including empty ones) + * covered by particles read via getNextParticle() for this thread. + * + * Uses one of two computation methods depending on format capabilities: + * - RATIO: The represented history count is available from the header. + * Maps the per-thread represented count to original histories via the + * known O/R ratio. + * - INCREMENTAL: Uses the direct sum of per-particle incremental history + * values. For formats with native per-particle data this gives exact + * original history counts. For formats without native data the default + * of 1 per new-history particle is used, so the result equals the + * number of represented (non-empty) histories read by this thread. + * + * @param threadIndex The index of the thread (0 to numThreads-1) + * @return Number of original histories covered by this thread so far + * @throws std::out_of_range If threadIndex is invalid + */ + std::uint64_t getHistoriesRead(size_t threadIndex) const; + + /** + * @brief Gets the total number of particles read across all threads. + * + * @return Total number of particles read + */ + std::uint64_t getTotalParticlesRead() const; + + /** + * @brief Gets the total number of original histories read across all threads. + * + * @return Total number of original histories read + */ + std::uint64_t getTotalHistoriesRead() const; + /** * @brief Gets the total number of particles in the phase space file. * @@ -156,6 +194,20 @@ namespace ParticleZoo { */ std::size_t getNumberOfThreads() const { return numThreads_; } + /** + * @brief Checks if the underlying phase space format provides native represented history count. + * + * @return True if native represented history count is available, false otherwise + */ + bool hasNativeRepresentedHistoryCount() const; + + /** + * @brief Checks if the underlying phase space format provides native incremental history counters. + * + * @return True if native incremental history counters are available, false otherwise + */ + bool hasNativeIncrementalHistoryCounters() const; + /** * @brief Closes all underlying phase space file readers. * @@ -165,17 +217,30 @@ namespace ParticleZoo { void close(); private: + /// @brief Computation mode for per-thread original history counting + enum class HistoryCountMode { RATIO, INCREMENTAL }; + std::string filename_; UserOptions options_; std::size_t numThreads_; + struct ThreadStatistics { + mutable std::mutex mutex; // For thread-safe updates + std::atomic particlesRead = 0; + std::atomic representedHistoriesRead = 0; + std::atomic incrementalHistorySum = 0; + }; + std::vector> readers_; std::vector startingParticleIndex_; - std::vector particlesRead_; + std::vector> threadStats_; + bool hasNativeRepresentedHistoryCount_; + bool hasNativeIncrementalHistoryCounters_; std::uint64_t numberOfOriginalHistories_; std::uint64_t numberOfParticlesInPhsp_; std::uint64_t numberOfRepresentedHistories_; + HistoryCountMode historyCountMode_; }; } \ No newline at end of file diff --git a/include/particlezoo/peneasy/penEasyphspFile.h b/include/particlezoo/peneasy/penEasyphspFile.h index 1273e3e..b4bf9b5 100644 --- a/include/particlezoo/peneasy/penEasyphspFile.h +++ b/include/particlezoo/peneasy/penEasyphspFile.h @@ -133,6 +133,12 @@ namespace ParticleZoo * @return Sum of all delta-N values from particle records */ std::uint64_t getNumberOfOriginalHistories() const override; + + /** + * @brief Check if this format stores incremental history numbers per-particle. + * @return true — penEasy always stores delta-N per particle + */ + bool hasNativeIncrementalHistoryCounters() const override; protected: /** @@ -167,6 +173,7 @@ namespace ParticleZoo // Inline implementations for the Reader class inline std::uint64_t Reader::getNumberOfParticles() const { return numberOfParticles_; } inline std::uint64_t Reader::getNumberOfOriginalHistories() const { return numberOfOriginalHistories_; } + inline bool Reader::hasNativeIncrementalHistoryCounters() const { return true; } inline size_t Reader::getMaximumASCIILineLength() const { return MAX_ASCII_LINE_LENGTH; } } // namespace penEasyphspFile diff --git a/include/particlezoo/utilities/version.h b/include/particlezoo/utilities/version.h index 10010fb..cc689c8 100644 --- a/include/particlezoo/utilities/version.h +++ b/include/particlezoo/utilities/version.h @@ -67,7 +67,7 @@ namespace ParticleZoo { * that do not add new features. Applications should always be able to * upgrade to higher patch versions safely. */ - static constexpr const int PATCH_VERSION = 0; + static constexpr const int PATCH_VERSION = 1; /** * @brief Development status indicator. diff --git a/python/pyproject.toml b/python/pyproject.toml index 97a49ae..8ef0ebd 100644 --- a/python/pyproject.toml +++ b/python/pyproject.toml @@ -8,7 +8,7 @@ build-backend = "setuptools.build_meta" [project] name = "particlezoo" -version = "1.1.0" +version = "1.1.1" description = "Python bindings for the ParticleZoo C++ library" readme = "README.md" requires-python = ">=3.8" diff --git a/python/setup.cfg b/python/setup.cfg index a0908e6..99de51d 100644 --- a/python/setup.cfg +++ b/python/setup.cfg @@ -1,6 +1,6 @@ [metadata] name = particlezoo -version = 1.1.0 +version = 1.1.1 [options] packages = find: diff --git a/python/setup.py b/python/setup.py index 7401d7a..15b1a85 100644 --- a/python/setup.py +++ b/python/setup.py @@ -24,6 +24,9 @@ str(Path("..") / "src" / "topas" / "TOPASphspFile.cc"), str(Path("..") / "src" / "IAEA" / "IAEAHeader.cc"), str(Path("..") / "src" / "IAEA" / "IAEAphspFile.cc"), + # Parallel readers + str(Path("..") / "src" / "parallel" / "HistoryBalancedParallelReader.cc"), + str(Path("..") / "src" / "parallel" / "ParticleBalancedParallelReader.cc"), ] define_macros = [("PYBIND11_DETAILED_ERROR_MESSAGES", "1")] diff --git a/python/src/particlezoo/__init__.py b/python/src/particlezoo/__init__.py index 44133db..6edd09a 100644 --- a/python/src/particlezoo/__init__.py +++ b/python/src/particlezoo/__init__.py @@ -1,6 +1,10 @@ from ._pz import ( + __version__, + get_version_string, Particle, ParticleType, + HistoryBalancedParallelReader, + ParticleBalancedParallelReader, IntPropertyType, FloatPropertyType, BoolPropertyType, @@ -71,8 +75,12 @@ ) __all__ = [ + "__version__", + "get_version_string", "Particle", "ParticleType", + "HistoryBalancedParallelReader", + "ParticleBalancedParallelReader", "IntPropertyType", "FloatPropertyType", "BoolPropertyType", diff --git a/python/src/pybind/module.cpp b/python/src/pybind/module.cpp index aba8558..b9ffb05 100644 --- a/python/src/pybind/module.cpp +++ b/python/src/pybind/module.cpp @@ -8,6 +8,9 @@ #include "particlezoo/utilities/formats.h" #include "particlezoo/utilities/argParse.h" #include "particlezoo/utilities/units.h" +#include "particlezoo/utilities/version.h" +#include "particlezoo/parallel/HistoryBalancedParallelReader.h" +#include "particlezoo/parallel/ParticleBalancedParallelReader.h" #include "particlezoo/egs/EGSLATCH.h" #include "particlezoo/penelope/ILBArray.h" @@ -17,6 +20,19 @@ using namespace ParticleZoo; PYBIND11_MODULE(_pz, m) { m.doc() = "Python bindings for ParticleZoo core and IAEA reader"; + // ===== Version ===== + { + std::string ver = std::to_string(Version::MAJOR_VERSION) + "." + + std::to_string(Version::MINOR_VERSION) + "." + + std::to_string(Version::PATCH_VERSION); + if (Version::CAVEAT[0] != '\0') { + ver += "." + std::string(Version::CAVEAT); + } + m.attr("__version__") = ver; + } + m.def("get_version_string", &Version::GetVersionString, + "Get the full ParticleZoo version string (e.g. 'ParticleZoo v1.1.1')."); + // Ensure built-in formats are registered before any factory calls try { FormatRegistry::RegisterStandardFormats(); } catch (...) { /* idempotent; ignore */ } @@ -479,6 +495,82 @@ PYBIND11_MODULE(_pz, m) { "Returns PhaseSpaceFileWriter instance configured for the specified format. " "Use fixed_values to specify constant properties for all particles."); + // ===== Parallel Readers ===== + + py::class_(m, "HistoryBalancedParallelReader", + "Multi-threaded phase space file reader that distributes histories evenly across threads. " + "Each thread maintains its own reader positioned at a specific section of the file. " + "Histories are the fundamental unit of work; particle counts may vary between threads.") + .def(py::init(), + py::arg("filename"), py::arg("options") = UserOptions{}, py::arg("num_threads") = 1, + "Create a history-balanced parallel reader. " + "Partitions represented histories evenly across the specified number of threads.") + .def("peek_next_particle", &HistoryBalancedParallelReader::peekNextParticle, + py::arg("thread_index"), + "Peek at the next particle for a thread without consuming it.") + .def("get_next_particle", &HistoryBalancedParallelReader::getNextParticle, + py::arg("thread_index"), + "Read and return the next particle for a thread.") + .def("has_more_particles", &HistoryBalancedParallelReader::hasMoreParticles, + py::arg("thread_index"), + "Check if more particles are available for a thread.") + .def("get_histories_read", &HistoryBalancedParallelReader::getHistoriesRead, + py::arg("thread_index"), + "Get the number of original histories processed by a thread (including empty ones).") + .def("get_particles_read", &HistoryBalancedParallelReader::getParticlesRead, + py::arg("thread_index"), + "Get the number of particles processed by a thread.") + .def("get_total_histories_read", &HistoryBalancedParallelReader::getTotalHistoriesRead, + "Get the total number of original histories read across all threads.") + .def("get_number_of_particles", &HistoryBalancedParallelReader::getNumberOfParticles, + "Get the total number of particles in the phase space file.") + .def("get_number_of_original_histories", &HistoryBalancedParallelReader::getNumberOfOriginalHistories, + "Get the number of original histories in the file (including empty ones).") + .def("get_number_of_represented_histories", &HistoryBalancedParallelReader::getNumberOfRepresentedHistories, + "Get the number of histories that produced at least one particle.") + .def("get_number_of_threads", &HistoryBalancedParallelReader::getNumberOfThreads, + "Get the number of threads used by this reader.") + .def("close", &HistoryBalancedParallelReader::close, + "Close all underlying readers and release resources.") + ; + + py::class_(m, "ParticleBalancedParallelReader", + "Multi-threaded phase space file reader that distributes particles evenly across threads. " + "Each thread maintains its own reader positioned at a specific section of the file. " + "Particles are the fundamental unit of work; history counts may vary between threads.") + .def(py::init(), + py::arg("filename"), py::arg("options") = UserOptions{}, py::arg("num_threads") = 1, + "Create a particle-balanced parallel reader. " + "Partitions particles evenly across the specified number of threads.") + .def("peek_next_particle", &ParticleBalancedParallelReader::peekNextParticle, + py::arg("thread_index"), + "Peek at the next particle for a thread without consuming it.") + .def("get_next_particle", &ParticleBalancedParallelReader::getNextParticle, + py::arg("thread_index"), + "Read and return the next particle for a thread.") + .def("has_more_particles", &ParticleBalancedParallelReader::hasMoreParticles, + py::arg("thread_index"), + "Check if more particles are available for a thread.") + .def("get_particles_read", &ParticleBalancedParallelReader::getParticlesRead, + py::arg("thread_index"), + "Get the number of particles processed by a thread.") + .def("get_histories_read", &ParticleBalancedParallelReader::getHistoriesRead, + py::arg("thread_index"), + "Get the number of original histories covered by a thread.") + .def("get_total_histories_read", &ParticleBalancedParallelReader::getTotalHistoriesRead, + "Get the total number of original histories read across all threads.") + .def("get_number_of_particles", &ParticleBalancedParallelReader::getNumberOfParticles, + "Get the total number of particles in the phase space file.") + .def("get_number_of_original_histories", &ParticleBalancedParallelReader::getNumberOfOriginalHistories, + "Get the number of original histories in the file (including empty ones).") + .def("get_number_of_represented_histories", &ParticleBalancedParallelReader::getNumberOfRepresentedHistories, + "Get the number of histories that produced at least one particle.") + .def("get_number_of_threads", &ParticleBalancedParallelReader::getNumberOfThreads, + "Get the number of threads used by this reader.") + .def("close", &ParticleBalancedParallelReader::close, + "Close all underlying readers and release resources.") + ; + // Expose a simple alias to create a Particle by PDG code m.def("particle_from_pdg", [](int pdg, float kineticEnergy, float x, float y, float z, float px, float py, float pz, bool isNewHistory, float weight){ ParticleType t = getParticleTypeFromPDGID(static_cast(pdg)); diff --git a/src/PhaseSpaceFileWriter.cc b/src/PhaseSpaceFileWriter.cc index 8b578d6..c345217 100644 --- a/src/PhaseSpaceFileWriter.cc +++ b/src/PhaseSpaceFileWriter.cc @@ -44,7 +44,10 @@ namespace ParticleZoo historiesToAccountFor_(0), buffer_(BUFFER_SIZE), writeParticleDepth_(0), - fixedValues_(fixedValues) + fixedValues_(fixedValues), + flipXDirection_(false), + flipYDirection_(false), + flipZDirection_(false) { if (formatType != FormatType::NONE && !file_.is_open()) { diff --git a/src/egs/egsphspFile.cc b/src/egs/egsphspFile.cc index e3f5acd..9378c56 100644 --- a/src/egs/egsphspFile.cc +++ b/src/egs/egsphspFile.cc @@ -145,7 +145,7 @@ namespace ParticleZoo::EGSphspFile // Writer class implementation Writer::Writer(const std::string & fileName, const UserOptions & options) - : PhaseSpaceFileWriter("EGS", fileName, options) + : PhaseSpaceFileWriter("EGS", fileName, options), latchOption_(EGSLATCHOPTION::LATCH_OPTION_2) { mode_ = EGSMODE::MODE0; // Default mode, MODE2 requires source particles to include ZLAST information @@ -160,6 +160,23 @@ namespace ParticleZoo::EGSphspFile throw std::runtime_error("Unsupported EGS phase-space file mode: " + modeStr); } } + + if (options.contains(EGSLATCHOptionCommand)) { + int latchOptionInt = options.extractIntOption(EGSLATCHOptionCommand); + switch (latchOptionInt) { + case 1: + latchOption_ = EGSLATCHOPTION::LATCH_OPTION_1; + break; + case 2: + latchOption_ = EGSLATCHOPTION::LATCH_OPTION_2; + break; + case 3: + latchOption_ = EGSLATCHOPTION::LATCH_OPTION_3; + break; + default: + throw std::runtime_error("Unsupported EGS LATCH option: " + std::to_string(latchOptionInt)); + } + } } std::vector Writer::getFormatSpecificCLICommands() { diff --git a/src/parallel/HistoryBalancedParallelReader.cc b/src/parallel/HistoryBalancedParallelReader.cc index 7209591..2e2b1ee 100644 --- a/src/parallel/HistoryBalancedParallelReader.cc +++ b/src/parallel/HistoryBalancedParallelReader.cc @@ -3,8 +3,6 @@ #include "particlezoo/utilities/formats.h" -#include - namespace ParticleZoo { HistoryBalancedParallelReader::HistoryBalancedParallelReader(const std::string& filename, const UserOptions& options, size_t numThreads) @@ -24,12 +22,15 @@ namespace ParticleZoo { readers_.emplace_back(std::move(reader)); } + // Detect if the format provides native represented history count or incremental history counters + hasNativeRepresentedHistoryCount_ = readers_[0]->hasNativeRepresentedHistoryCount(); + hasNativeIncrementalHistoryCounters_ = readers_[0]->hasNativeIncrementalHistoryCounters(); + // Determine the number of represented histories - try { - // Attempt to get the number of represented histories directly + if (hasNativeRepresentedHistoryCount_) { numberOfRepresentedHistories_ = readers_[0]->getNumberOfRepresentedHistories(); - } catch (const std::runtime_error&) { - // Manually count the number of represented histories if not supported + } else { + // Manually count the number of represented histories by scanning the file numberOfRepresentedHistories_ = 0; auto reader = FormatRegistry::CreateReader(filename, options); while (reader->hasMoreParticles()) { @@ -47,9 +48,6 @@ namespace ParticleZoo { // Get the number of original histories and adjust totalHistoriesToRead_ if necessary numberOfOriginalHistories_ = readers_[0]->getNumberOfOriginalHistories(); - if (numberOfOriginalHistories_ < totalHistoriesToRead_) { - totalHistoriesToRead_ = numberOfOriginalHistories_; - } // Get the number of particles in the phase space numberOfParticlesInPhsp_ = readers_[0]->getNumberOfParticles(); @@ -64,25 +62,25 @@ namespace ParticleZoo { if (hasGapsBetweenHistories_) { emptyHistoriesBetweenEachHistory_ = numberOfEmptyHistories / numberOfRepresentedHistories_; perHistoryErrorContribution_ = numberOfEmptyHistories % numberOfRepresentedHistories_; - emptyHistoryErrors_.resize(numThreads); } // Calculate starting history for each thread startingHistorys_.reserve(numThreads); - historiesRead_.reserve(numThreads); - totalHistoriesRead_.resize(numThreads, 0); - hasMoreParticlesCache_.resize(numThreads, NEEDS_CHECKING); + threadStats_.reserve(numThreads); + for (size_t i = 0; i < numThreads; ++i) { + threadStats_.emplace_back(std::make_unique()); + } std::uint64_t historiesPerThread = numberOfRepresentedHistories_ / numThreads; std::uint64_t remainderHistories = numberOfRepresentedHistories_ % numThreads; std::uint64_t currentStartingHistory = 0; for (size_t i = 0; i < numThreads; ++i) { startingHistorys_.push_back(currentStartingHistory); - historiesRead_.push_back(currentStartingHistory); + threadStats_[i]->historiesRead = currentStartingHistory; // Calculate the correct initial error for this thread based on how many // represented histories have been "processed" before this thread's starting point if (hasGapsBetweenHistories_) { std::uint64_t initialError = numberOfRepresentedHistories_ / 2; - emptyHistoryErrors_[i] = (initialError + currentStartingHistory * perHistoryErrorContribution_) % numberOfRepresentedHistories_; + threadStats_[i]->emptyHistoryError = (initialError + currentStartingHistory * perHistoryErrorContribution_) % numberOfRepresentedHistories_; } currentStartingHistory += historiesPerThread; if (i < remainderHistories) { @@ -90,35 +88,33 @@ namespace ParticleZoo { } } - // In separate threads, move each reader to its starting history - std::vector threads; - for (size_t i = 0; i < numThreads; ++i) { - threads.emplace_back([this, i]() { - std::uint64_t targetHistory = startingHistorys_[i]; - // Thread 0 starts at the beginning, no positioning needed - if (targetHistory == 0) { - return; - } - std::uint64_t currentHistory = 0; - // Skip to the first particle of history #targetHistory (0-indexed) - // We need to consume all particles from histories 0 through targetHistory-1 - while (readers_[i]->hasMoreParticles()) { - const Particle particle = readers_[i]->peekNextParticle(); + // Single-pass scan to find the particle index at each thread's starting history. + if (numThreads > 1) { + std::vector startingParticleIndices(numThreads, 0); + { + auto scanner = FormatRegistry::CreateReader(filename, options); + std::uint64_t historyCount = 0; + std::uint64_t particleIndex = 0; + size_t nextThreadToFind = 1; // Thread 0 always starts at particle 0 + + while (nextThreadToFind < numThreads && scanner->hasMoreParticles()) { + const Particle particle = scanner->getNextParticle(); if (particle.isNewHistory()) { - if (currentHistory == targetHistory) { - // We've reached the start of the target history, stop here - break; + if (historyCount == startingHistorys_[nextThreadToFind]) { + startingParticleIndices[nextThreadToFind] = particleIndex; + nextThreadToFind++; } - currentHistory++; + historyCount++; } - readers_[i]->getNextParticle(); // Consume this particle + particleIndex++; } - }); - } + scanner->close(); + } - // Join threads - for (auto& thread : threads) { - thread.join(); + // Position each reader at its starting particle using direct seek + for (size_t i = 1; i < numThreads; ++i) { + readers_[i]->moveToParticle(startingParticleIndices[i]); + } } } @@ -129,8 +125,8 @@ namespace ParticleZoo { } // check the cache first - if (hasMoreParticlesCache_[threadIndex] != NEEDS_CHECKING) { - return hasMoreParticlesCache_[threadIndex] == HAS_MORE_PARTICLES; + if (threadStats_[threadIndex]->hasMoreParticlesCache != NEEDS_CHECKING) { + return threadStats_[threadIndex]->hasMoreParticlesCache == HAS_MORE_PARTICLES; } // Check if the reader has more particles @@ -144,7 +140,7 @@ namespace ParticleZoo { : numberOfRepresentedHistories_; // Check if we have completed all histories for this thread - if (historiesRead_[threadIndex] >= targetHistories) { + if (threadStats_[threadIndex]->historiesRead >= targetHistories) { // Check if the next particle would start a new history const Particle nextParticle = readers_[threadIndex]->peekNextParticle(); // If the next particle starts a new history, we have no more particles for this thread @@ -153,7 +149,7 @@ namespace ParticleZoo { } // Update the cache - hasMoreParticlesCache_[threadIndex] = hasMoreParticles ? HAS_MORE_PARTICLES : NO_MORE_PARTICLES; + threadStats_[threadIndex]->hasMoreParticlesCache = hasMoreParticles ? HAS_MORE_PARTICLES : NO_MORE_PARTICLES; return hasMoreParticles; } @@ -168,34 +164,37 @@ namespace ParticleZoo { } // reset the cache since we are consuming a particle - hasMoreParticlesCache_[threadIndex] = NEEDS_CHECKING; + threadStats_[threadIndex]->hasMoreParticlesCache = NEEDS_CHECKING; // Get the next particle from the appropriate reader Particle particle = readers_[threadIndex]->getNextParticle(); // Update history count if this particle starts a new history int32_t incrementalHistoryNumber = 0; - if (particle.isNewHistory()) { - if (hasGapsBetweenHistories_) { - emptyHistoryErrors_[threadIndex] += perHistoryErrorContribution_; - if (emptyHistoryErrors_[threadIndex] >= numberOfRepresentedHistories_) { - incrementalHistoryNumber = 2 + emptyHistoriesBetweenEachHistory_; - emptyHistoryErrors_[threadIndex] -= numberOfRepresentedHistories_; + { + std::lock_guard lock(threadStats_[threadIndex]->mutex); + threadStats_[threadIndex]->particlesRead++; + + if (particle.isNewHistory()) { + if (hasGapsBetweenHistories_) { + threadStats_[threadIndex]->emptyHistoryError += perHistoryErrorContribution_; + if (threadStats_[threadIndex]->emptyHistoryError >= numberOfRepresentedHistories_) { + incrementalHistoryNumber = 2 + emptyHistoriesBetweenEachHistory_; + threadStats_[threadIndex]->emptyHistoryError -= numberOfRepresentedHistories_; + } else { + incrementalHistoryNumber = 1 + emptyHistoriesBetweenEachHistory_; + } } else { - incrementalHistoryNumber = 1 + emptyHistoriesBetweenEachHistory_; + incrementalHistoryNumber = 1; } - } else { - incrementalHistoryNumber = 1; + threadStats_[threadIndex]->historiesRead++; } - historiesRead_[threadIndex]++; - } else { - incrementalHistoryNumber = 0; + + threadStats_[threadIndex]->totalHistoriesRead += incrementalHistoryNumber; } particle.setIntProperty(IntPropertyType::INCREMENTAL_HISTORY_NUMBER, incrementalHistoryNumber); - totalHistoriesRead_[threadIndex] += incrementalHistoryNumber; - return particle; } @@ -217,7 +216,28 @@ namespace ParticleZoo { if (threadIndex >= readers_.size()) { throw std::out_of_range("Thread index out of range in getHistoriesRead()"); } - return totalHistoriesRead_[threadIndex]; + return threadStats_[threadIndex]->totalHistoriesRead; + } + + std::uint64_t HistoryBalancedParallelReader::getParticlesRead(size_t threadIndex) const { + if (threadIndex >= readers_.size()) { + throw std::out_of_range("Thread index out of range in getParticlesRead()"); + } + return threadStats_[threadIndex]->particlesRead; + } + + std::uint64_t HistoryBalancedParallelReader::getTotalParticlesRead() const { + std::uint64_t total = 0; + for (size_t t = 0; t < readers_.size(); t++) + total += threadStats_[t]->particlesRead; + return total; + } + + std::uint64_t HistoryBalancedParallelReader::getTotalHistoriesRead() const { + std::uint64_t total = 0; + for (size_t t = 0; t < readers_.size(); t++) + total += threadStats_[t]->totalHistoriesRead; + return total; } void HistoryBalancedParallelReader::close() { diff --git a/src/parallel/ParticleBalancedParallelReader.cc b/src/parallel/ParticleBalancedParallelReader.cc index 6963736..98a4bc0 100644 --- a/src/parallel/ParticleBalancedParallelReader.cc +++ b/src/parallel/ParticleBalancedParallelReader.cc @@ -33,12 +33,32 @@ namespace ParticleZoo { // Leave numberOfRepresentedHistories_ undetermined for now, lazily compute it later if needed numberOfRepresentedHistories_ = 0; + // Detect if the format provides native represented history count or incremental history counters + hasNativeRepresentedHistoryCount_ = readers_[0]->hasNativeRepresentedHistoryCount(); + hasNativeIncrementalHistoryCounters_ = readers_[0]->hasNativeIncrementalHistoryCounters(); + + // Detect the appropriate history counting mode based on format capabilities + if (hasNativeRepresentedHistoryCount_) { + // R is available cheaply from the header — use ratio approach + numberOfRepresentedHistories_ = readers_[0]->getNumberOfRepresentedHistories(); + historyCountMode_ = HistoryCountMode::RATIO; + } else { + // Use incremental sums. For formats with native per-particle incremental + // history data this gives exact original history counts. For other formats + // getIncrementalHistories() defaults to 1 per new-history particle, so the + // sum equals the number of represented (non-empty) histories read. + historyCountMode_ = HistoryCountMode::INCREMENTAL; + } + std::uint64_t particlesPerThread = numberOfParticlesInPhsp_ / numThreads; std::uint64_t remainderParticles = numberOfParticlesInPhsp_ % numThreads; // Position each reader at its starting particle startingParticleIndex_.reserve(numThreads); - particlesRead_.resize(numThreads, 0); + threadStats_.reserve(numThreads); + for (size_t i = 0; i < numThreads; ++i) { + threadStats_.emplace_back(std::make_unique()); + } std::uint64_t currentParticleIndex = 0; for (size_t i = 0; i < numThreads; ++i) { std::uint64_t particlesToRead = particlesPerThread + (i < remainderParticles ? 1 : 0); @@ -71,7 +91,7 @@ namespace ParticleZoo { ? startingParticleIndex_[threadIndex + 1] : numberOfParticlesInPhsp_; // Check if we have completed all particles for this thread - if (particlesRead_[threadIndex] >= (targetParticleIndex - startingParticleIndex_[threadIndex])) { + if (threadStats_[threadIndex]->particlesRead >= (targetParticleIndex - startingParticleIndex_[threadIndex])) { hasMore = false; } } @@ -90,7 +110,17 @@ namespace ParticleZoo { // Get the next particle from the appropriate reader Particle particle = readers_[threadIndex]->getNextParticle(); - particlesRead_[threadIndex]++; + + { + std::lock_guard lock(threadStats_[threadIndex]->mutex); + threadStats_[threadIndex]->particlesRead++; + + if (particle.isNewHistory()) { + threadStats_[threadIndex]->representedHistoriesRead++; + threadStats_[threadIndex]->incrementalHistorySum += particle.getIncrementalHistories(); + } + } + return particle; } @@ -110,10 +140,10 @@ namespace ParticleZoo { std::uint64_t ParticleBalancedParallelReader::getNumberOfRepresentedHistories() { if (numberOfRepresentedHistories_ == 0) { // Lazily compute the number of represented histories - try { - // Attempt to get the number of represented histories directly + if (readers_[0]->hasNativeRepresentedHistoryCount()) { + // Get the number of represented histories directly numberOfRepresentedHistories_ = readers_[0]->getNumberOfRepresentedHistories(); - } catch (const std::runtime_error&) { + } else { // Manually count the number of represented histories if not supported numberOfRepresentedHistories_ = 0; auto reader = FormatRegistry::CreateReader(filename_, options_); @@ -134,7 +164,67 @@ namespace ParticleZoo { if (threadIndex >= readers_.size()) { throw std::out_of_range("Thread index out of range in getHistoriesRead()"); } - return particlesRead_[threadIndex]; + return threadStats_[threadIndex]->particlesRead; + } + + std::uint64_t ParticleBalancedParallelReader::getHistoriesRead(size_t threadIndex) const { + if (threadIndex >= numThreads_) { + throw std::out_of_range("Thread index out of range in getHistoriesRead()"); + } + + switch (historyCountMode_) { + case HistoryCountMode::RATIO: + { + const std::uint64_t representedHistories = threadStats_[threadIndex]->representedHistoriesRead; + return (representedHistories * numberOfOriginalHistories_) + / numberOfRepresentedHistories_; + } + case HistoryCountMode::INCREMENTAL: + return threadStats_[threadIndex]->incrementalHistorySum; + } + + return threadStats_[threadIndex]->representedHistoriesRead; // unreachable fallback + } + + std::uint64_t ParticleBalancedParallelReader::getTotalParticlesRead() const { + std::uint64_t total = 0; + for (size_t t = 0; t < numThreads_; t++) { + total += threadStats_[t]->particlesRead; + } + return total; + } + + std::uint64_t ParticleBalancedParallelReader::getTotalHistoriesRead() const { + switch (historyCountMode_) { + case HistoryCountMode::RATIO: { + std::uint64_t totalRep = 0; + for (size_t t = 0; t < numThreads_; t++) { + totalRep += threadStats_[t]->representedHistoriesRead; + } + return (totalRep * numberOfOriginalHistories_) / numberOfRepresentedHistories_; + } + case HistoryCountMode::INCREMENTAL: { + std::uint64_t total = 0; + std::uint64_t totalParticles = 0; + for (size_t t = 0; t < numThreads_; t++) { + std::lock_guard lock(threadStats_[t]->mutex); + total += threadStats_[t]->incrementalHistorySum; + totalParticles += threadStats_[t]->particlesRead; + } + if (totalParticles >= numberOfParticlesInPhsp_) + return numberOfOriginalHistories_; + return total; + } + } + return 0; + } + + bool ParticleBalancedParallelReader::hasNativeRepresentedHistoryCount() const { + return hasNativeRepresentedHistoryCount_; + } + + bool ParticleBalancedParallelReader::hasNativeIncrementalHistoryCounters() const { + return hasNativeIncrementalHistoryCounters_; } void ParticleBalancedParallelReader::close() {