From 95becff9e48c833589ee6eeb5442f0a0dfc27124 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 4 Jul 2026 11:52:39 -0400 Subject: [PATCH 01/50] Added a development code-timing helper --- include/ScopeTimer.H | 64 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) create mode 100644 include/ScopeTimer.H diff --git a/include/ScopeTimer.H b/include/ScopeTimer.H new file mode 100644 index 000000000..fe638ad48 --- /dev/null +++ b/include/ScopeTimer.H @@ -0,0 +1,64 @@ +#pragma once + +#include +#include +#include + +class ScopeTimer { +public: + using Clock = std::chrono::steady_clock; + + // Define units for output + enum class Unit { + Seconds, + Milliseconds, + Microseconds, + Nanoseconds + }; + + // Default to Seconds if the user doesn't choose + explicit ScopeTimer(std::string task_name, Unit output_unit = Unit::Seconds) + : name(std::move(task_name)), unit(output_unit), start_time(Clock::now()) {} + + // Output on deletion of the object + ~ScopeTimer() + { + const auto end_time = Clock::now(); + const auto delta = end_time - start_time; + + // Convert natively using std::chrono based on user choice + switch (unit) { + case Unit::Seconds: { + std::chrono::duration elapsed = delta; + std::cout << "[" << name << "] completed in " << elapsed.count() << " s\n"; + break; + } + case Unit::Milliseconds: { + std::chrono::duration elapsed = delta; + std::cout << "[" << name << "] completed in " << elapsed.count() << " ms\n"; + break; + } + case Unit::Microseconds: { + std::chrono::duration elapsed = delta; + std::cout << "[" << name << "] completed in " << elapsed.count() << " µs\n"; + break; + } + case Unit::Nanoseconds: { + std::chrono::duration elapsed = delta; + std::cout << "[" << name << "] completed in " << elapsed.count() << " ns\n"; + break; + } + default: + std::cerr << "Unknown time unit for ScopeTimer\n"; + break; + } + } + + ScopeTimer(const ScopeTimer&) = delete; + ScopeTimer& operator=(const ScopeTimer&) = delete; + +private: + std::string name; + Unit unit; + std::chrono::time_point start_time; +}; From 435565de6cf8334407f047c0e8d00a19fd403f49 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 4 Jul 2026 11:52:53 -0400 Subject: [PATCH 02/50] Added bods2hdf5 --- utils/PhaseSpace/CMakeLists.txt | 3 +- utils/PhaseSpace/exp2hdf5.cc | 577 ++++++++++++++++++++++++++++++++ 2 files changed, 579 insertions(+), 1 deletion(-) create mode 100644 utils/PhaseSpace/exp2hdf5.cc diff --git a/utils/PhaseSpace/CMakeLists.txt b/utils/PhaseSpace/CMakeLists.txt index 809835df7..2f48326c2 100644 --- a/utils/PhaseSpace/CMakeLists.txt +++ b/utils/PhaseSpace/CMakeLists.txt @@ -3,7 +3,7 @@ set(bin_PROGRAMS pspinfo pspstat psp2ascii ascii2psp psp2histo psp2histo1d snap2histo1d psp2histoE psp2histoH pspreal pspstatS psp2range psp2histoT psp2histoTC pspstatT spl2psp ascii2psp2 psp2rings psp2bess psp2lagu psp2interp diffpsp psp2hdf5 pspmono - snap2ascii psporbv) + snap2ascii psporbv bods2hdf5) if(HAVE_VTK) list(APPEND bin_PROGRAMS psp2eHisto psp2vtu) @@ -70,6 +70,7 @@ add_executable(diffpsp diffpsp.cc MakeModel.cc PSP.cc KDE2d.cc) add_executable(pspmono pspmono.cc MakeModel.cc PSP.cc) add_executable(psp2hdf5 psp2hdf5.cc PSP.cc) add_executable(psporbv psporbv.cc PSP.cc) +add_executable(bods2hdf5 exp2hdf5.cc) if(HAVE_VTK) diff --git a/utils/PhaseSpace/exp2hdf5.cc b/utils/PhaseSpace/exp2hdf5.cc new file mode 100644 index 000000000..0f151bb8b --- /dev/null +++ b/utils/PhaseSpace/exp2hdf5.cc @@ -0,0 +1,577 @@ +/* +Draft EXP ascii particle phase-space data HDF5 converter + +This is a possible template for implementing HDF5 input files for +EXP. The code reads an ASCII file with particle data and writes it to +an HDF5 file using the HighFive library. It also provides a function +to read the HDF5 file back into ASCII format, allowing for round-trip +conversion. + +Compilation notes +----------------- + +To compile this code, you need to have HighFive and HDF5 development +libraries installed. It also uses the 'cxxopts.H' and 'ScopeTimer.H' +which need to be local. + +Example compilation command for Ubuntu/Debian systems: + +g++ -std=c++17 -I/usr/include/hdf5/serial exp2hdf5.cc -o exp2hdf5 -L/usr/lib/x86_64-linux-gnu/hdf5/serial -lhdf5_hl -lhdf5 -Wl,-rpath,/usr/lib/x86_64-linux-gnu/hdf5/serial + +Key features +------------ + +1. Reads ASCII header with 3 integers for particle count, aux int + fields, and aux float fields +2. Parses flexible number of auxiliary integer and floating-point fields +3. Uses HighFive for a simple C++ HDF5 interface +4. Applies GZIP compression (level 4) for reduced file size. +5. Preserves numerical precision with double-precision floats if desired. + This is less relevant for simulation initial conditions. +6. Stores metadata as HDF5 attributes for round-trip integrity +7. Comprehensive std error handling with informative messages +8. Supports both float32 and float64 precision with automatic detection +9. Round-trip conversion maintains numerical accuracy when using float64 + +Performance considerations +-------------------------- + +- Uses chunked storage with automatic chunk sizing (10th of dataset or 1024, + whichever is larger). +- Uses shuffle filter to improve compression efficiency. +- GZIP compression level 9 provides good compression ratio. Level 5 seems + good enough for most cases. +- For very large datasets, consider tuning chunk size in DataSetCreateProps +- Memory usage is O(N) where N = total number of values across all fields + +File structure in HDF5 +---------------------- +/ +├── Attributes +│ ├── num_particles (int) +│ ├── num_aux_ints (int) +│ └── num_aux_floats (int) +└── particles/ (Group) +├── m (Dataset) +├── x, y, z (Datasets) +├── u, v, w (Datasets) +├── aux_int_0, aux_int_1, ... (Datasets) +└── aux_float_0, aux_float_1, ... (Datasets) + +This schemea could (should be?) changed to use the Gadget style. On +the other hand, this design is clean and specific to exp. Also, users +can easily wrie and read this with (e.g.) h5py. I'm incliened to keep +it. +*/ + +// C++ std +#include +#include +#include +#include +#include +#include +#include +#include + +// HighFive +#include +#include +#include + +// Command-line parsing +#include "cxxopts.H" + +// For command timing +#include "ScopeTimer.H" + +using namespace HighFive; + +// Specify floating point precision +enum class FloatPrecision { FLOAT32, FLOAT64 }; + +// Type variant for flexible float handling +template +struct ParticleDataTemplate +{ + std::vector m; // mass/species index + std::vector x, y, z; // position + std::vector u, v, w; // velocity + std::vector> aux_ints; // auxiliary integer fields + std::vector> aux_floats; // auxiliary float fields + + size_t num_particles = 0; + size_t num_aux_ints = 0; + size_t num_aux_floats = 0; +}; + +// Read ASCII and write HDF5 with specified precision +template +void ascii_to_hdf5_impl(const std::string& ascii_file, + const std::string& hdf5_file, + FloatPrecision precision) +{ + std::ifstream infile(ascii_file); + if (!infile.is_open()) { + throw std::runtime_error("Could not open ASCII file: " + ascii_file); + } + + // Read header + int num_particles, num_aux_ints, num_aux_floats; + if (!(infile >> num_particles >> num_aux_ints >> num_aux_floats)) { + throw std::runtime_error("Failed to read header from ASCII file"); + } + + // Validate header + if (num_particles <= 0 || num_aux_ints < 0 || num_aux_floats < 0) { + throw std::runtime_error("Invalid header values"); + } + + // Initialize particle data structure + ParticleDataTemplate data; + + data.num_particles = num_particles; + data.num_aux_ints = num_aux_ints; + data.num_aux_floats = num_aux_floats; + + data.m.resize(num_particles); + data.x.resize(num_particles); + data.y.resize(num_particles); + data.z.resize(num_particles); + data.u.resize(num_particles); + data.v.resize(num_particles); + data.w.resize(num_particles); + + data.aux_ints.resize(num_aux_ints); + for (auto& vec : data.aux_ints) { + vec.resize(num_particles); + } + + data.aux_floats.resize(num_aux_floats); + for (auto& vec : data.aux_floats) { + vec.resize(num_particles); + } + + // Read particle data + for (int i = 0; i < num_particles; ++i) { + if (!(infile >> data.m[i] >> data.x[i] >> data.y[i] >> data.z[i] + >> data.u[i] >> data.v[i] >> data.w[i])) { + throw std::runtime_error("Failed to read particle data at row " + + std::to_string(i)); + } + + for (int j = 0; j < num_aux_ints; ++j) { + if (!(infile >> data.aux_ints[j][i])) { + throw std::runtime_error("Failed to read aux int field at row " + + std::to_string(i)); + } + } + + for (int j = 0; j < num_aux_floats; ++j) { + if (!(infile >> data.aux_floats[j][i])) { + throw std::runtime_error("Failed to read aux float field at row " + + std::to_string(i)); + } + } + } + infile.close(); + + // Create HDF5 file with compression enabled + File file(hdf5_file, File::ReadWrite | File::Create | File::Truncate); + + // Store header information and precision metadata + file.createAttribute("num_particles", DataSpace::From(num_particles)) + .write(num_particles); + file.createAttribute("num_aux_ints", DataSpace::From(num_aux_ints)) + .write(num_aux_ints); + file.createAttribute("num_aux_floats", DataSpace::From(num_aux_floats)) + .write(num_aux_floats); + + // Store precision type: 0 = float32, 1 = float64 + int precision_flag = (precision == FloatPrecision::FLOAT32) ? 0 : 1; + file.createAttribute("float_precision", DataSpace::From(precision_flag)) + .write(precision_flag); + + // Create a group for particle data + Group particles_group = file.createGroup("particles"); + + // Define compression filters + DataSetCreateProps props; + props.add(Chunking(std::vector{(hsize_t)std::max(num_particles / 10, 1024)})); + props.add(Shuffle()); + props.add(Deflate(4)); // This is a good compromise between + // speed and compression ratio + + // Write core phase-space fields + particles_group.createDataSet("m", data.m, props); + particles_group.createDataSet("x", data.x, props); + particles_group.createDataSet("y", data.y, props); + particles_group.createDataSet("z", data.z, props); + particles_group.createDataSet("u", data.u, props); + particles_group.createDataSet("v", data.v, props); + particles_group.createDataSet("w", data.w, props); + + // Write auxiliary integer fields + for (size_t j = 0; j < data.aux_ints.size(); ++j) { + std::string dset_name = "aux_int_" + std::to_string(j); + particles_group.createDataSet(dset_name, data.aux_ints[j], props); + } + + // Write auxiliary float fields + for (size_t j = 0; j < data.aux_floats.size(); ++j) { + std::string dset_name = "aux_float_" + std::to_string(j); + particles_group.createDataSet(dset_name, data.aux_floats[j], props); + } + + std::string precision_str = (precision == FloatPrecision::FLOAT32) ? "float32" : "float64"; + std::cout << "Successfully wrote " << num_particles << " particles to " + << hdf5_file << " (" << precision_str << ")" << std::endl; +} + +// Dispatcher function - user specifies precision +void ascii_to_hdf5(const std::string& ascii_file, + const std::string& hdf5_file, + FloatPrecision precision = FloatPrecision::FLOAT64) +{ + if (precision == FloatPrecision::FLOAT32) { + ascii_to_hdf5_impl(ascii_file, hdf5_file, precision); + } else { + ascii_to_hdf5_impl(ascii_file, hdf5_file, precision); + } +} + +// Variant for flexible reading +using FloatData = std::variant, std::vector>; + +struct ParticleDataVariant +{ + FloatData m, x, y, z, u, v, w; + std::vector> aux_ints; + std::vector aux_floats; + + size_t num_particles = 0; + size_t num_aux_ints = 0; + size_t num_aux_floats = 0; + FloatPrecision precision; +}; + +// Read HDF5 with automatic precision detection +ParticleDataVariant read_hdf5_data(const std::string& hdf5_file) +{ + File file(hdf5_file, File::ReadOnly); + + // Read metadata + int num_particles = file.getAttribute("num_particles").read(); + int num_aux_ints = file.getAttribute("num_aux_ints").read(); + int num_aux_floats = file.getAttribute("num_aux_floats").read(); + + // Auto-detect precision + int precision_flag = 1; // Default to float64 + try { + precision_flag = file.getAttribute("float_precision").read(); + } catch (...) { + // If attribute doesn't exist, assume float64 (backward compatibility) + } + FloatPrecision precision = (precision_flag == 0) ? FloatPrecision::FLOAT32 + : FloatPrecision::FLOAT64; + + // Open particles group + Group particles_group = file.getGroup("particles"); + + ParticleDataVariant data; + data.num_particles = num_particles; + data.num_aux_ints = num_aux_ints; + data.num_aux_floats = num_aux_floats; + data.precision = precision; + + // Lambda to read float or double based on precision + auto read_dataset = [&particles_group, precision](const std::string& name) -> FloatData { + if (precision == FloatPrecision::FLOAT32) { + return particles_group.getDataSet(name).read>(); + } else { + return particles_group.getDataSet(name).read>(); + } + }; + + // Read core PSP fields + data.m = read_dataset("m"); + data.x = read_dataset("x"); + data.y = read_dataset("y"); + data.z = read_dataset("z"); + data.u = read_dataset("u"); + data.v = read_dataset("v"); + data.w = read_dataset("w"); + + // Read auxiliary integer fields + data.aux_ints.resize(num_aux_ints); + for (int j = 0; j < num_aux_ints; ++j) { + std::string dset_name = "aux_int_" + std::to_string(j); + data.aux_ints[j] = particles_group.getDataSet(dset_name) + .read>(); + } + + // Read auxiliary float fields + data.aux_floats.resize(num_aux_floats); + for (int j = 0; j < num_aux_floats; ++j) { + std::string dset_name = "aux_float_" + std::to_string(j); + data.aux_floats[j] = read_dataset(dset_name); + } + + return data; +} + +// Write ASCII from variant data using std::visit +// +void hdf5_to_ascii(const std::string& hdf5_file, const std::string& ascii_file) +{ + ParticleDataVariant data = read_hdf5_data(hdf5_file); + + std::ofstream outfile(ascii_file); + if (!outfile.is_open()) { + throw std::runtime_error("Could not open ASCII file for writing: " + ascii_file); + } + + // Write header + outfile << data.num_particles << " " << data.num_aux_ints + << " " << data.num_aux_floats << "\n"; + + outfile.precision(16); // High precision for floats + + // Lambda visitor to handle both float and double + auto write_value = [](const FloatData& variant_data, size_t index) { + return std::visit([index](const auto& vec) { return static_cast(vec[index]); }, + variant_data); + }; + + // Write particle data + int num_particles = data.num_particles; + for (int i = 0; i < num_particles; ++i) { + // Write core phase-space fields using visitor + outfile << write_value(data.m, i) << " " + << write_value(data.x, i) << " " + << write_value(data.y, i) << " " + << write_value(data.z, i) << " " + << write_value(data.u, i) << " " + << write_value(data.v, i) << " " + << write_value(data.w, i); + + // Write auxiliary integer fields + for (int j = 0; j < data.num_aux_ints; ++j) { + outfile << " " << data.aux_ints[j][i]; + } + + // Write auxiliary float fields using visitor + for (int j = 0; j < data.num_aux_floats; ++j) { + outfile << " " << write_value(data.aux_floats[j], i); + } + outfile << "\n"; + } + outfile.close(); + + std::string precision_str = (data.precision == FloatPrecision::FLOAT32) + ? "float32" : "float64"; + std::cout << "Successfully wrote " << num_particles << " particles to " + << ascii_file << " (" << precision_str << ")" << std::endl; +} + +int main(int argc, char* argv[]) +{ + std::string prefix = "particles", output; + std::string suffix = "bods"; + std::string ascii_restored = "particles_restored.txt"; + + // Parse command-line arguments for input/output files and mode + // + cxxopts::Options options("exp2hdf5", "ASCII to HDF5 particle converter with round-trip support and precision handling"); + + options.add_options() + ("i,input", "Input prefix", cxxopts::value(prefix)->default_value("particles")) + ("o,output", "Output prefix (optional, otherwise input prefix is used)", cxxopts::value(prefix)) + ("a,suffix", "Input suffix", cxxopts::value(suffix)->default_value("bods")) + ("roundtrip", "Perform round-trip conversion (ASCII -> HDF5 -> ASCII)") + ("verify", "Verify that the restored ASCII file matches the original in 'roundtrip' mode") + ("to_hdf5", "Convert ASCII to HDF5") + ("double", "Use double precision for HDF5 output (float is default)") + ("to_ascii", "Convert HDF5 to ASCII"); + ("h,help", "Print usage"); + + cxxopts::ParseResult vm; + + try { + vm = options.parse(argc, argv); + } catch (cxxopts::OptionException& e) { + std::cout << "Option error: " << e.what() << std::endl; + exit(-1); + } + + if (vm.count("help")) { + std::cout << options.help() << std::endl; + return 0; + } + + // This is for testing and timing + // + if (vm.count("roundtrip")) { + + // Time the whole test + ScopeTimer timer("full test"); + + // Filenames + std::string ascii = prefix + "." + suffix; + std::string h5_32 = prefix + "_f32.h5"; + std::string h5_64 = prefix + "_f64.h5"; + std::string rest = prefix + "_restored." + suffix; + + try { + // Convert to HDF5 with float32 precision + std::cout << "=== Converting ASCII to HDF5 (float32) ===" << std::endl; + { + ScopeTimer timer("float 32 conversion"); + ascii_to_hdf5(ascii, h5_32, FloatPrecision::FLOAT32); + } + + // Convert to HDF5 with float64 precision + std::cout << "\n=== Converting ASCII to HDF5 (float64) ===" << std::endl; + { + ScopeTimer timer("float 64 conversion"); + ascii_to_hdf5(ascii, h5_64, FloatPrecision::FLOAT64); + } + + // Round-trip with float64 version + std::cout << "\n=== Converting HDF5 (float64) back to ASCII ===" << std::endl; + { + ScopeTimer timer("hdf5 to ascii"); + hdf5_to_ascii(h5_64, rest); + } + + std::cout << "\n=== Round-trip conversion complete! ===" << std::endl; + + if (vm.count("verify")) { + std::cout << "\n=== Verifying restored ASCII file against original ===" << std::endl; + std::ifstream original(ascii); + std::ifstream restored(rest ); + + if (!original.is_open() || !restored.is_open()) { + throw std::runtime_error("Could not open files for verification"); + } + + std::string line_orig, line_rest; + size_t line_num = 0; + bool mismatch_found = false; + + std::vector max_diff(7, 0.0); // For m, x, y, z, u, v, w + std::vector vec_orig(7), vec_rest(7); + + while (std::getline(original, line_orig) && std::getline(restored, line_rest)) { + ++line_num; + + std::istringstream ss_orig(line_orig), ss_rest(line_rest); + + if (line_num == 1) { + int total_particles_orig, total_aux_ints_orig, total_aux_floats_orig; + int total_particles_rest, total_aux_ints_rest, total_aux_floats_rest; + + ss_orig >> total_particles_orig >> total_aux_ints_orig >> total_aux_floats_orig; + ss_rest >> total_particles_rest >> total_aux_ints_rest >> total_aux_floats_rest; + + if (total_particles_orig != total_particles_rest || + total_aux_ints_orig != total_aux_ints_rest || + total_aux_floats_orig != total_aux_floats_rest) { + std::cerr << "Header mismatch at line 1:\n" + << "Original: " << line_orig << "\n" + << "Restored: " << line_rest << "\n"; + mismatch_found = true; + break; + } + + continue; + } + + for (int i = 0; i < 7; ++i) { + ss_orig >> vec_orig[i]; + ss_rest >> vec_rest[i]; + } + + for (int i = 0; i < 7; ++i) { + if (std::abs(vec_orig[i] - vec_rest[i]) > max_diff[i]) { + max_diff[i] = std::abs(vec_orig[i] - vec_rest[i]); + } + + } + } + + std::cout << "Maximum absolute differences for core fields" + << std::endl << "(m, x, y, z, u, v, w): " + << max_diff[0] << ", " + << max_diff[1] << ", " + << max_diff[2] << ", " + << max_diff[3] << ", " + << max_diff[4] << ", " + << max_diff[5] << ", " + << max_diff[6] << std::endl; + } + + } catch (const std::exception& e) { + std::cerr << "Error: " << e.what() << std::endl; + return 1; + } + + + std::cout << "\n=== Compression ratios ===" << std::endl; + + namespace fs = std::filesystem; + + try { + fs::path ascPath = ascii; + fs::path f32Path = h5_32; + fs::path f64Path = h5_64; + + // Returns size in bytes (std::uintmax_t) + std::uintmax_t ascSize = fs::file_size(ascPath); + std::uintmax_t f32Size = fs::file_size(f32Path); + std::uintmax_t f64Size = fs::file_size(f64Path); + + // Print compression ratios + std::cout << "Ascii/hdf5(float) = " + << std::setprecision(3) + << static_cast(ascSize)/f32Size + << std::endl; + + std::cout << "Ascii/hdf5(double) = " + << std::setprecision(3) + << static_cast(ascSize)/f64Size + << std::endl; + std::cout << std::endl; + } + catch (const fs::filesystem_error& e) { + std::cerr << "Error: " << e.what() << '\n'; + } + + return 0; + } + + // Production mode: convert ASCII to HDF5 or vice versa based on + // mode selection + +if (vm.count("to_hdf5")) { + ScopeTimer timer("conversion to hdf5"); + std::string ascii = prefix + "." + suffix; + std::string ofile = prefix + ".h5"; + if (vm.count("output")) ofile = output + ".h5"; + if (vm.count("double")) + ascii_to_hdf5(ascii, ofile, FloatPrecision::FLOAT64); + else + ascii_to_hdf5(ascii, ofile, FloatPrecision::FLOAT32); + } else if (vm.count("to_ascii")) { + ScopeTimer timer("conversion to ascii"); + std::string ifile = prefix + ".h5"; + std::string ofile = prefix + "." + suffix; + hdf5_to_ascii(ifile, ofile); + } else { + std::cerr << "No conversion mode specified. Use --ascii_to_hdf5 or --hdf5_to_ascii." << std::endl; + return 1; + } + + return 0; +} + From 17108898e94f1237fce7393237fdfc19af6b4dee Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 4 Jul 2026 12:11:14 -0400 Subject: [PATCH 03/50] Change the executable and source file name to hdf5bods --- utils/PhaseSpace/CMakeLists.txt | 4 ++-- utils/PhaseSpace/{exp2hdf5.cc => hdf5bods.cc} | 23 +++++++++++++++++++ 2 files changed, 25 insertions(+), 2 deletions(-) rename utils/PhaseSpace/{exp2hdf5.cc => hdf5bods.cc} (96%) diff --git a/utils/PhaseSpace/CMakeLists.txt b/utils/PhaseSpace/CMakeLists.txt index 2f48326c2..b951f1cd9 100644 --- a/utils/PhaseSpace/CMakeLists.txt +++ b/utils/PhaseSpace/CMakeLists.txt @@ -3,7 +3,7 @@ set(bin_PROGRAMS pspinfo pspstat psp2ascii ascii2psp psp2histo psp2histo1d snap2histo1d psp2histoE psp2histoH pspreal pspstatS psp2range psp2histoT psp2histoTC pspstatT spl2psp ascii2psp2 psp2rings psp2bess psp2lagu psp2interp diffpsp psp2hdf5 pspmono - snap2ascii psporbv bods2hdf5) + snap2ascii psporbv hdf5bods) if(HAVE_VTK) list(APPEND bin_PROGRAMS psp2eHisto psp2vtu) @@ -70,7 +70,7 @@ add_executable(diffpsp diffpsp.cc MakeModel.cc PSP.cc KDE2d.cc) add_executable(pspmono pspmono.cc MakeModel.cc PSP.cc) add_executable(psp2hdf5 psp2hdf5.cc PSP.cc) add_executable(psporbv psporbv.cc PSP.cc) -add_executable(bods2hdf5 exp2hdf5.cc) +add_executable(hdf5bods hdf5bods.cc) if(HAVE_VTK) diff --git a/utils/PhaseSpace/exp2hdf5.cc b/utils/PhaseSpace/hdf5bods.cc similarity index 96% rename from utils/PhaseSpace/exp2hdf5.cc rename to utils/PhaseSpace/hdf5bods.cc index 0f151bb8b..d71dae863 100644 --- a/utils/PhaseSpace/exp2hdf5.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -401,6 +401,29 @@ int main(int argc, char* argv[]) vm = options.parse(argc, argv); } catch (cxxopts::OptionException& e) { std::cout << "Option error: " << e.what() << std::endl; + + // Append custom examples + std::cout << R"( +Examples: + Convert a standard EXP input body file nameed 'mybods.bods' to hdf5 format + $ hdf5bods --to_hdf5 -i mybods + The resulting HDF5 file using float32 internally will be called 'mybods.h5' + + Do the same conversion but use full double precision + $ hdf5bods --to_hdf5 --double -i mybods + The resulting HDF5 file will use float64 internally + + Convert the HDF5 file back to the standard EXP ascii format: + $ hdf5bods --to_ascii -i input + + The suffix on input can be customized. For example, the following converts + a standard EXP input body file nameed 'mybods.asc' to hdf5 format + $ hdf5bods --to_hdf5 --suffix=asc -i mybods + The resulting HDF5 file using float32 format will be called 'mybods.h5' + +)" << std::endl; + + exit(-1); } From bea3f850d78f0c02eaf85d0ce6f600b52b4739a5 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 4 Jul 2026 13:13:39 -0400 Subject: [PATCH 04/50] Added some OpenMP parallelism for a bit more throughput --- utils/PhaseSpace/hdf5bods.cc | 300 ++++++++++++++++++++++++----------- 1 file changed, 210 insertions(+), 90 deletions(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index d71dae863..dfaeb490e 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -14,9 +14,11 @@ To compile this code, you need to have HighFive and HDF5 development libraries installed. It also uses the 'cxxopts.H' and 'ScopeTimer.H' which need to be local. -Example compilation command for Ubuntu/Debian systems: +Example standalone compilation command for Ubuntu/Debian systems: -g++ -std=c++17 -I/usr/include/hdf5/serial exp2hdf5.cc -o exp2hdf5 -L/usr/lib/x86_64-linux-gnu/hdf5/serial -lhdf5_hl -lhdf5 -Wl,-rpath,/usr/lib/x86_64-linux-gnu/hdf5/serial +g++ -std=c++17 -I/usr/include/hdf5/serial hdf5bods.cc -o hdf5bods -L/usr/lib/x86_64-linux-gnu/hdf5/serial -lhdf5_hl -lhdf5 -Wl,-rpath,/usr/lib/x86_64-linux-gnu/hdf5/serial + +It is currently build in the utils/PhaseSpace directory of the EXP source tree. Key features ------------ @@ -65,14 +67,20 @@ it. */ // C++ std +#include #include -#include -#include -#include +#include #include #include +#include +#include #include -#include +#include +#include +#include + +// For OpenMP parallelization +#include // HighFive #include @@ -105,11 +113,47 @@ struct ParticleDataTemplate size_t num_aux_floats = 0; }; -// Read ASCII and write HDF5 with specified precision +// Helper: Parse a single line of particle data template -void ascii_to_hdf5_impl(const std::string& ascii_file, - const std::string& hdf5_file, - FloatPrecision precision) +void parse_particle_line(const std::string& line, + int particle_idx, + int num_aux_ints, + int num_aux_floats, + ParticleDataTemplate& data) +{ + std::istringstream iss(line); + T val; + + // Parse core PSP fields + if (!(iss >> data.m[particle_idx] + >> data.x[particle_idx] >> data.y[particle_idx] >> data.z[particle_idx] + >> data.u[particle_idx] >> data.v[particle_idx] >> data.w[particle_idx])) { + throw std::runtime_error("Failed to parse core fields at particle " + + std::to_string(particle_idx)); + } + + // Parse auxiliary integer fields + for (int j = 0; j < num_aux_ints; ++j) { + if (!(iss >> data.aux_ints[j][particle_idx])) { + throw std::runtime_error("Failed to parse aux int at particle " + + std::to_string(particle_idx)); + } + } + + // Parse auxiliary float fields + for (int j = 0; j < num_aux_floats; ++j) { + if (!(iss >> data.aux_floats[j][particle_idx])) { + throw std::runtime_error("Failed to parse aux float at particle " + + std::to_string(particle_idx)); + } + } +} + +// Read entire ASCII file into memory +std::vector read_ascii_lines(const std::string& ascii_file, + int& num_particles, + int& num_aux_ints, + int& num_aux_floats) { std::ifstream infile(ascii_file); if (!infile.is_open()) { @@ -117,7 +161,6 @@ void ascii_to_hdf5_impl(const std::string& ascii_file, } // Read header - int num_particles, num_aux_ints, num_aux_floats; if (!(infile >> num_particles >> num_aux_ints >> num_aux_floats)) { throw std::runtime_error("Failed to read header from ASCII file"); } @@ -127,6 +170,39 @@ void ascii_to_hdf5_impl(const std::string& ascii_file, throw std::runtime_error("Invalid header values"); } + // Read all lines + std::vector lines; + lines.reserve(num_particles); + std::string line; + while (std::getline(infile, line)) { + if (!line.empty() && line[0] != '#') { // Skip empty/comment lines + lines.push_back(line); + } + } + infile.close(); + + if ((int)lines.size() != num_particles) { + throw std::runtime_error("Mismatch between header count and actual particle count"); + } + + return lines; +} + +// Read ASCII and write HDF5 with specified precision +template +void ascii_to_hdf5_impl(const std::string& ascii_file, + const std::string& hdf5_file, + FloatPrecision precision, + bool verbose) +{ + // Header values + int num_particles, num_aux_ints, num_aux_floats; + + // Read all ASCII lines into memory + if (verbose) std::cout << "Reading ASCII file..." << std::endl; + std::vector lines = read_ascii_lines(ascii_file, num_particles, + num_aux_ints, num_aux_floats); + // Initialize particle data structure ParticleDataTemplate data; @@ -152,31 +228,24 @@ void ascii_to_hdf5_impl(const std::string& ascii_file, vec.resize(num_particles); } - // Read particle data + // Parse particle data in parallel + if (verbose) std::cout << "Parsing particles (parallel)..." << std::endl; + #pragma omp parallel for schedule(dynamic, 256) \ + num_threads(omp_get_max_threads()) for (int i = 0; i < num_particles; ++i) { - if (!(infile >> data.m[i] >> data.x[i] >> data.y[i] >> data.z[i] - >> data.u[i] >> data.v[i] >> data.w[i])) { - throw std::runtime_error("Failed to read particle data at row " + - std::to_string(i)); - } - - for (int j = 0; j < num_aux_ints; ++j) { - if (!(infile >> data.aux_ints[j][i])) { - throw std::runtime_error("Failed to read aux int field at row " + - std::to_string(i)); - } - } - - for (int j = 0; j < num_aux_floats; ++j) { - if (!(infile >> data.aux_floats[j][i])) { - throw std::runtime_error("Failed to read aux float field at row " + - std::to_string(i)); + try { + parse_particle_line(lines[i], i, num_aux_ints, num_aux_floats, data); + } catch (const std::exception& e) { + // Error handling in parallel region + #pragma omp critical + { + std::cerr << "Thread error at particle " << i << ": " << e.what() << std::endl; } } } - infile.close(); - // Create HDF5 file with compression enabled + // Create HDF5 file with compression enabled (SEQUENTIAL - thread-safe) + if (verbose) std::cout << "Writing HDF5 file..." << std::endl; File file(hdf5_file, File::ReadWrite | File::Create | File::Truncate); // Store header information and precision metadata @@ -202,7 +271,7 @@ void ascii_to_hdf5_impl(const std::string& ascii_file, props.add(Deflate(4)); // This is a good compromise between // speed and compression ratio - // Write core phase-space fields + // Write datasets sequentially (HDF5 is not fully thread-safe for writing) particles_group.createDataSet("m", data.m, props); particles_group.createDataSet("x", data.x, props); particles_group.createDataSet("y", data.y, props); @@ -223,20 +292,23 @@ void ascii_to_hdf5_impl(const std::string& ascii_file, particles_group.createDataSet(dset_name, data.aux_floats[j], props); } - std::string precision_str = (precision == FloatPrecision::FLOAT32) ? "float32" : "float64"; - std::cout << "Successfully wrote " << num_particles << " particles to " - << hdf5_file << " (" << precision_str << ")" << std::endl; + if (verbose) { + std::string precision_str = (precision == FloatPrecision::FLOAT32) ? "float32" : "float64"; + std::cout << "Successfully wrote " << num_particles << " particles to " + << hdf5_file << " (" << precision_str << ")" << std::endl; + } } // Dispatcher function - user specifies precision void ascii_to_hdf5(const std::string& ascii_file, const std::string& hdf5_file, - FloatPrecision precision = FloatPrecision::FLOAT64) + FloatPrecision precision = FloatPrecision::FLOAT64, + bool verbose = true) { if (precision == FloatPrecision::FLOAT32) { - ascii_to_hdf5_impl(ascii_file, hdf5_file, precision); + ascii_to_hdf5_impl(ascii_file, hdf5_file, precision, verbose); } else { - ascii_to_hdf5_impl(ascii_file, hdf5_file, precision); + ascii_to_hdf5_impl(ascii_file, hdf5_file, precision, verbose); } } @@ -261,8 +333,8 @@ ParticleDataVariant read_hdf5_data(const std::string& hdf5_file) File file(hdf5_file, File::ReadOnly); // Read metadata - int num_particles = file.getAttribute("num_particles").read(); - int num_aux_ints = file.getAttribute("num_aux_ints").read(); + int num_particles = file.getAttribute("num_particles").read(); + int num_aux_ints = file.getAttribute("num_aux_ints").read(); int num_aux_floats = file.getAttribute("num_aux_floats").read(); // Auto-detect precision @@ -272,8 +344,9 @@ ParticleDataVariant read_hdf5_data(const std::string& hdf5_file) } catch (...) { // If attribute doesn't exist, assume float64 (backward compatibility) } + FloatPrecision precision = (precision_flag == 0) ? FloatPrecision::FLOAT32 - : FloatPrecision::FLOAT64; + : FloatPrecision::FLOAT64; // Open particles group Group particles_group = file.getGroup("particles"); @@ -322,7 +395,8 @@ ParticleDataVariant read_hdf5_data(const std::string& hdf5_file) // Write ASCII from variant data using std::visit // -void hdf5_to_ascii(const std::string& hdf5_file, const std::string& ascii_file) +void hdf5_to_ascii(const std::string& hdf5_file, const std::string& ascii_file, + bool verbose = true) { ParticleDataVariant data = read_hdf5_data(hdf5_file); @@ -332,68 +406,92 @@ void hdf5_to_ascii(const std::string& hdf5_file, const std::string& ascii_file) } // Write header + // outfile << data.num_particles << " " << data.num_aux_ints << " " << data.num_aux_floats << "\n"; + outfile.close(); // Close header write outfile.precision(16); // High precision for floats // Lambda visitor to handle both float and double - auto write_value = [](const FloatData& variant_data, size_t index) { - return std::visit([index](const auto& vec) { return static_cast(vec[index]); }, - variant_data); + // + auto write_value = [](const FloatData& variant_data, size_t index) + { + return std::visit([index](const auto& vec) + { return static_cast(vec[index]); }, variant_data); }; - // Write particle data - int num_particles = data.num_particles; - for (int i = 0; i < num_particles; ++i) { - // Write core phase-space fields using visitor - outfile << write_value(data.m, i) << " " - << write_value(data.x, i) << " " - << write_value(data.y, i) << " " - << write_value(data.z, i) << " " - << write_value(data.u, i) << " " - << write_value(data.v, i) << " " - << write_value(data.w, i); + // Strategy: Pre-build all output lines in parallel, then write sequentially + std::vector output_lines(data.num_particles); + + #pragma omp parallel for schedule(static, 256) \ + num_threads(omp_get_max_threads()) + for (int i = 0; i < (int)data.num_particles; ++i) { + std::ostringstream oss; + oss.precision(16); + + // Write core phase-space fields + oss << write_value(data.m, i) << " " + << write_value(data.x, i) << " " + << write_value(data.y, i) << " " + << write_value(data.z, i) << " " + << write_value(data.u, i) << " " + << write_value(data.v, i) << " " + << write_value(data.w, i); // Write auxiliary integer fields for (int j = 0; j < data.num_aux_ints; ++j) { - outfile << " " << data.aux_ints[j][i]; + oss << " " << data.aux_ints[j][i]; } - // Write auxiliary float fields using visitor + // Write auxiliary float fields for (int j = 0; j < data.num_aux_floats; ++j) { - outfile << " " << write_value(data.aux_floats[j], i); + oss << " " << write_value(data.aux_floats[j], i); } - outfile << "\n"; + + output_lines[i] = oss.str(); + } + + // Write all lines sequentially to avoid I/O contention + outfile.open(ascii_file, std::ios::app); + for (const auto& line : output_lines) { + outfile << line << "\n"; } outfile.close(); std::string precision_str = (data.precision == FloatPrecision::FLOAT32) ? "float32" : "float64"; - std::cout << "Successfully wrote " << num_particles << " particles to " - << ascii_file << " (" << precision_str << ")" << std::endl; + if (verbose) + std::cout << "Successfully wrote " << data.num_particles << " particles to " + << ascii_file << " (" << precision_str << ")" << std::endl; } +// Main function with command-line parsing +// int main(int argc, char* argv[]) { std::string prefix = "particles", output; std::string suffix = "bods"; std::string ascii_restored = "particles_restored.txt"; + int num_threads = omp_get_max_threads(); + bool quiet = false; // Parse command-line arguments for input/output files and mode // - cxxopts::Options options("exp2hdf5", "ASCII to HDF5 particle converter with round-trip support and precision handling"); + cxxopts::Options options(argv[0], "ASCII to HDF5 particle converter with round-trip support and precision handling"); options.add_options() - ("i,input", "Input prefix", cxxopts::value(prefix)->default_value("particles")) - ("o,output", "Output prefix (optional, otherwise input prefix is used)", cxxopts::value(prefix)) - ("a,suffix", "Input suffix", cxxopts::value(suffix)->default_value("bods")) + ("i,input", "Input prefix", cxxopts::value(prefix)->default_value("particles")) + ("o,output", "Output prefix (optional, otherwise input prefix is used)", cxxopts::value(prefix)) + ("a,suffix", "Input suffix", cxxopts::value(suffix)->default_value("bods")) + ("t,threads", "Number of OpenMP threads (default: max available)", cxxopts::value(num_threads)->default_value(std::to_string(num_threads))) ("roundtrip", "Perform round-trip conversion (ASCII -> HDF5 -> ASCII)") - ("verify", "Verify that the restored ASCII file matches the original in 'roundtrip' mode") - ("to_hdf5", "Convert ASCII to HDF5") - ("double", "Use double precision for HDF5 output (float is default)") - ("to_ascii", "Convert HDF5 to ASCII"); - ("h,help", "Print usage"); + ("verify", "Verify that the restored ASCII file matches the original in 'roundtrip' mode") + ("to_hdf5", "Convert ASCII to HDF5") + ("double", "Use double precision for HDF5 output (float is default)") + ("to_ascii", "Convert HDF5 to ASCII") + ("q,quiet", "Suppress verbose output") + ("h,help", "Print usage"); cxxopts::ParseResult vm; @@ -401,16 +499,22 @@ int main(int argc, char* argv[]) vm = options.parse(argc, argv); } catch (cxxopts::OptionException& e) { std::cout << "Option error: " << e.what() << std::endl; + exit(-1); + } + + if (vm.count("help")) { + std::cout << options.help() << std::endl; // Append custom examples std::cout << R"( Examples: + Convert a standard EXP input body file nameed 'mybods.bods' to hdf5 format $ hdf5bods --to_hdf5 -i mybods The resulting HDF5 file using float32 internally will be called 'mybods.h5' - Do the same conversion but use full double precision - $ hdf5bods --to_hdf5 --double -i mybods + Do the same conversion but use full double precision without being chatty + $ hdf5bods --quiet --to_hdf5 --double -i mybods The resulting HDF5 file will use float64 internally Convert the HDF5 file back to the standard EXP ascii format: @@ -421,23 +525,34 @@ int main(int argc, char* argv[]) $ hdf5bods --to_hdf5 --suffix=asc -i mybods The resulting HDF5 file using float32 format will be called 'mybods.h5' + Test round-trip conversion and verify that the restored ASCII file matches the + original: + $ hdf5bods --roundtrip --verify -i mybods + )" << std::endl; + return 0; + } - exit(-1); + if (vm.count("quiet")) { + quiet = true; } - if (vm.count("help")) { - std::cout << options.help() << std::endl; - return 0; + if (vm.count("threads")) { + omp_set_num_threads(num_threads); } + + if (!quiet) + std::cout << "=== Using " << num_threads << " OpenMP threads ===" + << std::endl; // This is for testing and timing // if (vm.count("roundtrip")) { // Time the whole test - ScopeTimer timer("full test"); + std::unique_ptr time_ptr; + if (!quiet) time_ptr = std::make_unique("full test"); // Filenames std::string ascii = prefix + "." + suffix; @@ -449,22 +564,25 @@ int main(int argc, char* argv[]) // Convert to HDF5 with float32 precision std::cout << "=== Converting ASCII to HDF5 (float32) ===" << std::endl; { - ScopeTimer timer("float 32 conversion"); - ascii_to_hdf5(ascii, h5_32, FloatPrecision::FLOAT32); + std::unique_ptr time_ptr; + if (!quiet) time_ptr = std::make_unique("float 32 conversion"); + ascii_to_hdf5(ascii, h5_32, FloatPrecision::FLOAT32, !quiet); } // Convert to HDF5 with float64 precision std::cout << "\n=== Converting ASCII to HDF5 (float64) ===" << std::endl; { - ScopeTimer timer("float 64 conversion"); - ascii_to_hdf5(ascii, h5_64, FloatPrecision::FLOAT64); + std::unique_ptr time_ptr; + if (!quiet) time_ptr = std::make_unique("float 64 conversion"); + ascii_to_hdf5(ascii, h5_64, FloatPrecision::FLOAT64, !quiet); } // Round-trip with float64 version std::cout << "\n=== Converting HDF5 (float64) back to ASCII ===" << std::endl; { - ScopeTimer timer("hdf5 to ascii"); - hdf5_to_ascii(h5_64, rest); + std::unique_ptr time_ptr; + if (!quiet) time_ptr = std::make_unique("hdf5 to ascii"); + hdf5_to_ascii(h5_64, rest, !quiet); } std::cout << "\n=== Round-trip conversion complete! ===" << std::endl; @@ -575,18 +693,20 @@ int main(int argc, char* argv[]) // Production mode: convert ASCII to HDF5 or vice versa based on // mode selection - -if (vm.count("to_hdf5")) { - ScopeTimer timer("conversion to hdf5"); + // + if (vm.count("to_hdf5")) { + std::unique_ptr time_ptr; + if (!quiet) time_ptr = std::make_unique("conversion to hdf5"); std::string ascii = prefix + "." + suffix; std::string ofile = prefix + ".h5"; if (vm.count("output")) ofile = output + ".h5"; if (vm.count("double")) - ascii_to_hdf5(ascii, ofile, FloatPrecision::FLOAT64); + ascii_to_hdf5(ascii, ofile, FloatPrecision::FLOAT64, !quiet); else - ascii_to_hdf5(ascii, ofile, FloatPrecision::FLOAT32); + ascii_to_hdf5(ascii, ofile, FloatPrecision::FLOAT32, !quiet); } else if (vm.count("to_ascii")) { - ScopeTimer timer("conversion to ascii"); + std::unique_ptr time_ptr; + if (!quiet) time_ptr = std::make_unique("conversion to ascii"); std::string ifile = prefix + ".h5"; std::string ofile = prefix + "." + suffix; hdf5_to_ascii(ifile, ofile); From 89dcf59077cce70b25b64c1780138fc75ff0858a Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 4 Jul 2026 16:00:43 -0400 Subject: [PATCH 05/50] Read exp-style hdf5 input files in Component --- src/Component.H | 2 + src/Component.cc | 331 ++++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 332 insertions(+), 1 deletion(-) diff --git a/src/Component.H b/src/Component.H index 923796ca1..7adbff7ae 100644 --- a/src/Component.H +++ b/src/Component.H @@ -199,6 +199,8 @@ private: //! Parallel distribute and particle io void load_balance(void); void update_indices(void); + void read_bodies_and_distribute_init(void); + void read_bodies_and_distribute_hdf5(void); void read_bodies_and_distribute_ascii(void); void read_bodies_and_distribute_binary_out(istream *); void read_bodies_and_distribute_binary_spl(istream *); diff --git a/src/Component.cc b/src/Component.cc index c50884ce6..fb0adb7e9 100644 --- a/src/Component.cc +++ b/src/Component.cc @@ -2,11 +2,20 @@ #include #include #include +#include #include #include #include #include +// HighFive +#include +#include +#include + +// HDF5 C API +#include + #include "Component.H" #include "Bessel.H" #include "Sphere.H" @@ -277,7 +286,7 @@ Component::Component(YAML::Node& CONF) initialize_cuda(); - read_bodies_and_distribute_ascii(); + read_bodies_and_distribute_init(); reset_level_lists(); @@ -1463,6 +1472,326 @@ Component::~Component(void) delete [] acc0; } +void Component::read_bodies_and_distribute_init() +{ + // Check whether input file is HDF5 + + // Note: One might need to use H5Fis_hdf5(filename.c_str()) for + // older HDF5 versions (< 1.12) + +#if H5_VERSION_GE(1, 12, 0) + bool isHDF5 = H5Fis_accessible(pfile.c_str(), H5P_DEFAULT) > 0; +#else + bool isHDF5 = H5Fis_hdf5(pfile.c_str()) > 0; +#endif + + if (isHDF5) { + if (myid==0) + std::cout << "---- Component <" << name << ">: " + << "reading bodies from HDF5 file "<< pfile << std::endl; + read_bodies_and_distribute_hdf5(); + } else { + if (myid==0) + std::cout << "---- Component <" << name << ">: " + << "reading bodies from ASCII file "<< pfile << std::endl; + read_bodies_and_distribute_ascii(); + } +} + +void Component::read_bodies_and_distribute_hdf5(void) +{ + HighFive::File file(pfile, HighFive::File::ReadOnly); + + // Read metadata (all ranks read the same metadata, but only rank 0 reads the datasets) + // + nbodies_tot = file.getAttribute("num_particles") .read(); + niattrib = file.getAttribute("num_aux_ints") .read(); + ndattrib = file.getAttribute("num_aux_floats").read(); + + // Get precision flag for float datasets (0=float32, 1=float64) + // + int precision_flag = 1; // Default to float64 + try { + precision_flag = file.getAttribute("float_precision").read(); + } catch (...) { + std::ostringstream msg; + msg << "Component ERROR: 'float_precision' attribute not found in HDF5 file <" + << pfile << ">. Please ensure the file has this attribute set to 0 (float32) or 1 (float64)."; + throw std::runtime_error(msg.str()); + } + + // Open particles group + // + HighFive::Group particles_group = file.getGroup("particles"); + + // Variant for flexible reading and access of float or double datasets + // + using FloatDataVariant = std::variant, std::vector>; + + // Wrapper struct to provide uniform access to float or double std::vectors + // + struct FloatData + { + // The underlying data can be either a vector of floats or doubles + FloatDataVariant data; + + // Default constructor + FloatData() = default; + + // Implicit constructor so you can assign vectors directly to FloatData + FloatData(FloatDataVariant v) : data(std::move(v)) {} + + // Provides a uniform interface to access elements regardless of + // underlying type. Forces a value copy, converting float to + // double if necessary. + double operator[](size_t i) const { + return std::visit([i](const auto& vec) -> double { + return vec[i]; }, data); + } + }; + + // Phase-space data + // + FloatData m, x, y, z, u, v, w; + std::vector> aux_ints; + std::vector aux_floats; + + // Lambda to read float or double based on precision flag + // + auto read_dataset = [&particles_group, precision_flag] + (const std::string& name, size_t offset, size_t batch) -> FloatData { + + HighFive::DataSet dataset = particles_group.getDataSet(name); + + // Fetch the dimensions (e.g., [rows, cols]) + size_t total = dataset.getSpace().getDimensions()[0]; + + // Sanity check to ensure we don't read beyond the dataset + size_t current = std::min(batch, total - offset); + + // Wrap returned vector in FloatData + if (precision_flag == 0) { + return FloatData(dataset.select({offset}, {current}).read>()); + } + else { + return FloatData(dataset.select({offset}, {current}).read>()); + } + }; + + + // Okay, now we can read the datasets and distribute particles + // across MPI rank, following the previously established logic. The + // following code reads the datasets in batches and distributes them + // to the appropriate ranks. + + // For radius exclusion check + double rmax1=0.0, r2; + + is_init = 1; + setup_distribution(); + is_init = 0; + + // Initialize the particle ferry instance with dynamic attribute + // sizes + if (not pf) pf = ParticleFerryPtr(new ParticleFerry(niattrib, ndattrib)); + + + if (myid==0) { + + // Offset into HDF5 datasets for each node, read in batchs for all + // datasets + size_t offset = 0; + + // Read Node 0's PSP fields + m = read_dataset("m", offset, nbodies_table[0]); + + x = read_dataset("x", offset, nbodies_table[0]); + y = read_dataset("y", offset, nbodies_table[0]); + z = read_dataset("z", offset, nbodies_table[0]); + + u = read_dataset("u", offset, nbodies_table[0]); + v = read_dataset("v", offset, nbodies_table[0]); + w = read_dataset("w", offset, nbodies_table[0]); + + // Read auxiliary integer fields + aux_ints.resize(nbodies_table[0]); + for (int j = 0; j < niattrib; ++j) { + std::string dset_name = "aux_int_" + std::to_string(j); + aux_ints[j] = particles_group.getDataSet(dset_name). + select({offset}, {nbodies_table[0]}).read>(); + } + + // Read auxiliary float fields + aux_floats.resize(nbodies_table[0]); + for (int j = 0; j < ndattrib; ++j) { + std::string dset_name = "aux_float_" + std::to_string(j); + aux_floats[j] = read_dataset(dset_name, offset, nbodies_table[0]); + } + + for (unsigned i=0; i(niattrib, ndattrib); + + // Load the particle data for root + part->mass = m[i]; + + part->pos[0] = x[i]; + part->pos[1] = y[i]; + part->pos[2] = z[i]; + + part->vel[0] = u[i]; + part->vel[1] = v[i]; + part->vel[2] = w[i]; + + for (int j=0; jiattrib[j] = aux_ints[j][i]; + } + + for (int j=0; jdattrib[j] = aux_floats[j][i]; + } + + // Make a sequential index for the particle + part->indx = i + 1; + + // Get the radius + r2 = 0.0; + for (int k=0; k<3; k++) r2 += part->pos[k]*part->pos[k]; + rmax1 = max(r2, rmax1); + + // Load the particle + particles[part->indx] = part; + + // Record top_seq + top_seq = std::max(part->indx, top_seq); + } + + nbodies = nbodies_table[0]; + + // Offset into the next node's data in the HDF5 datasets + offset += nbodies_table[0]; + + unsigned icount, ibufcount; + for (int n=1; n>(); + } + + // Read auxiliary float fields + aux_floats.resize(nbodies_table[n]); + for (int j = 0; j < ndattrib; ++j) { + std::string dset_name = "aux_float_" + std::to_string(j); + aux_floats[j] = read_dataset(dset_name, offset, nbodies_table[n]); + } + + pf->ShipParticles(n, 0, nbodies_table[n]); + + icount = 0; + ibufcount = 0; + while (icount < nbodies_table[n]) { + + PartPtr part = std::make_shared(niattrib, ndattrib); + + // Make a sequential index for the particle + part->indx = nbodies_index[n-1] + 1 + icount; + + part->mass = m[icount]; + + part->pos[0] = x[icount]; + part->pos[1] = y[icount]; + part->pos[2] = z[icount]; + + part->vel[0] = u[icount]; + part->vel[1] = v[icount]; + part->vel[2] = w[icount]; + + r2 = 0.0; + for (int k=0; k<3; k++) r2 += part->pos[k]*part->pos[k]; + rmax1 = max(r2, rmax1); + + pf->SendParticle(part); + icount++; + + // Record top_seq + top_seq = std::max(part->indx, top_seq); + } + + // Update offset for next node + offset += nbodies_table[n]; + } + + } else { + + pf->ShipParticles(myid, 0, nbodies); + +#ifdef DEBUG + int icount = 0; +#endif + + while (PartPtr part=pf->RecvParticle()) { + particles[part->indx] = part; +#ifdef DEBUG + if (icount<5) { + cout << "Process " << myid << ": received "; + cout << setw(14) << part->mass; + for (int k=0; k<3; k++) cout << setw(14) << part->pos[k]; + cout << endl; + } + icount++; +#endif + } + } + + // Default: set to max radius + // can be overriden by parameter + + rmax = sqrt(fabs(rmax1)); + MPI_Bcast(&rmax, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + + // Send top_seq to all nodes + MPI_Bcast(&top_seq, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD); + + initialize(); + +#ifdef DEBUG + if (particles.size()) { + unsigned long imin = std::numeric_limits::max(); + unsigned long imax = 0, kmin = imax, kmax = 0; + for (auto p : particles) { + imin = std::min(imin, p.first); + imax = std::max(imax, p.first); + kmin = std::min(kmin, p.second->indx); + kmax = std::max(kmax, p.second->indx); + } + cout << "read_bodies_and_distribute_ascii: process " << myid + << " name=" << name << " bodies [" << kmin << ", " + << kmax << "], [" << imin << ", " << imax << "]" + << " #=" << particles.size() << endl; + } else { + cout << "read_bodies_and_distribute_ascii: process " << myid + << " name=" << name + << " #=" << particles.size() << endl; + } +#endif +} + + void Component::read_bodies_and_distribute_ascii(void) { // Open file From 264da15eb8f1ae7098997d3ed5980b3781f466c0 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 4 Jul 2026 16:15:14 -0400 Subject: [PATCH 06/50] Give the user a hint if no conversion mode is specified --- utils/PhaseSpace/hdf5bods.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index dfaeb490e..641f4aeba 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -711,7 +711,7 @@ int main(int argc, char* argv[]) std::string ofile = prefix + "." + suffix; hdf5_to_ascii(ifile, ofile); } else { - std::cerr << "No conversion mode specified. Use --ascii_to_hdf5 or --hdf5_to_ascii." << std::endl; + std::cerr << "No conversion mode specified. Use --to_hdf5 to convert to HDF5, --to_ascii to convert from HDF5 to ascii, or the --roundtrip test." << std::endl; return 1; } From eac730224e6796b3c45d37204f92f6f1caf01e3b Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 4 Jul 2026 17:01:37 -0400 Subject: [PATCH 07/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- src/Component.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Component.cc b/src/Component.cc index fb0adb7e9..f4acae953 100644 --- a/src/Component.cc +++ b/src/Component.cc @@ -1614,15 +1614,15 @@ void Component::read_bodies_and_distribute_hdf5(void) w = read_dataset("w", offset, nbodies_table[0]); // Read auxiliary integer fields - aux_ints.resize(nbodies_table[0]); + aux_ints.resize(niattrib); for (int j = 0; j < niattrib; ++j) { std::string dset_name = "aux_int_" + std::to_string(j); - aux_ints[j] = particles_group.getDataSet(dset_name). - select({offset}, {nbodies_table[0]}).read>(); + aux_ints[j] = particles_group.getDataSet(dset_name) + .select({offset}, {nbodies_table[0]}).read>(); } // Read auxiliary float fields - aux_floats.resize(nbodies_table[0]); + aux_floats.resize(ndattrib); for (int j = 0; j < ndattrib; ++j) { std::string dset_name = "aux_float_" + std::to_string(j); aux_floats[j] = read_dataset(dset_name, offset, nbodies_table[0]); From 67ffcd784a57b643ea1f502ab33d0479b24deb09 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 4 Jul 2026 17:02:07 -0400 Subject: [PATCH 08/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- src/Component.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Component.cc b/src/Component.cc index f4acae953..4e8a4aeab 100644 --- a/src/Component.cc +++ b/src/Component.cc @@ -1686,15 +1686,15 @@ void Component::read_bodies_and_distribute_hdf5(void) w = read_dataset("w", offset, nbodies_table[n]); // Read auxiliary integer fields - aux_ints.resize(nbodies_table[n]); + aux_ints.resize(niattrib); for (int j = 0; j < niattrib; ++j) { std::string dset_name = "aux_int_" + std::to_string(j); - aux_ints[j] = particles_group.getDataSet(dset_name). - select({offset}, {nbodies_table[n]}).read>(); + aux_ints[j] = particles_group.getDataSet(dset_name) + .select({offset}, {nbodies_table[n]}).read>(); } // Read auxiliary float fields - aux_floats.resize(nbodies_table[n]); + aux_floats.resize(ndattrib); for (int j = 0; j < ndattrib; ++j) { std::string dset_name = "aux_float_" + std::to_string(j); aux_floats[j] = read_dataset(dset_name, offset, nbodies_table[n]); From 8d3dc08882e8b4b0db4bceb95084565784b937b5 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 4 Jul 2026 17:03:05 -0400 Subject: [PATCH 09/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- src/Component.cc | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/Component.cc b/src/Component.cc index 4e8a4aeab..900e4c4ea 100644 --- a/src/Component.cc +++ b/src/Component.cc @@ -1566,9 +1566,12 @@ void Component::read_bodies_and_distribute_hdf5(void) // Fetch the dimensions (e.g., [rows, cols]) size_t total = dataset.getSpace().getDimensions()[0]; + if (offset >= total) { + throw std::runtime_error("Component ERROR: dataset <" + name + "> shorter than expected"); + } + // Sanity check to ensure we don't read beyond the dataset size_t current = std::min(batch, total - offset); - // Wrap returned vector in FloatData if (precision_flag == 0) { return FloatData(dataset.select({offset}, {current}).read>()); From 8205360f4823abbfbe1f2d912248873d2c0267f5 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 4 Jul 2026 17:06:47 -0400 Subject: [PATCH 10/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- src/Component.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Component.cc b/src/Component.cc index 900e4c4ea..0e888cdd0 100644 --- a/src/Component.cc +++ b/src/Component.cc @@ -1782,12 +1782,12 @@ void Component::read_bodies_and_distribute_hdf5(void) kmin = std::min(kmin, p.second->indx); kmax = std::max(kmax, p.second->indx); } - cout << "read_bodies_and_distribute_ascii: process " << myid + cout << "read_bodies_and_distribute_hdf5: process " << myid << " name=" << name << " bodies [" << kmin << ", " << kmax << "], [" << imin << ", " << imax << "]" << " #=" << particles.size() << endl; } else { - cout << "read_bodies_and_distribute_ascii: process " << myid + cout << "read_bodies_and_distribute_hdf5: process " << myid << " name=" << name << " #=" << particles.size() << endl; } From a4442ddb19f373ee0949fb0de1dd8e8ec99d55aa Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 4 Jul 2026 17:07:58 -0400 Subject: [PATCH 11/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/PhaseSpace/hdf5bods.cc | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index 641f4aeba..e79bdd5a8 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -230,13 +230,15 @@ void ascii_to_hdf5_impl(const std::string& ascii_file, // Parse particle data in parallel if (verbose) std::cout << "Parsing particles (parallel)..." << std::endl; - #pragma omp parallel for schedule(dynamic, 256) \ + + int parse_errors = 0; + #pragma omp parallel for schedule(dynamic, 256) reduction(+:parse_errors) \ num_threads(omp_get_max_threads()) for (int i = 0; i < num_particles; ++i) { try { parse_particle_line(lines[i], i, num_aux_ints, num_aux_floats, data); } catch (const std::exception& e) { - // Error handling in parallel region + parse_errors += 1; #pragma omp critical { std::cerr << "Thread error at particle " << i << ": " << e.what() << std::endl; @@ -244,6 +246,10 @@ void ascii_to_hdf5_impl(const std::string& ascii_file, } } + if (parse_errors > 0) { + throw std::runtime_error("Encountered " + std::to_string(parse_errors) + + " particle parse errors; aborting conversion."); + } // Create HDF5 file with compression enabled (SEQUENTIAL - thread-safe) if (verbose) std::cout << "Writing HDF5 file..." << std::endl; File file(hdf5_file, File::ReadWrite | File::Create | File::Truncate); From ac77cb6765c60ebe7dcb4a2428d83e891fdd8671 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 4 Jul 2026 17:08:47 -0400 Subject: [PATCH 12/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/PhaseSpace/hdf5bods.cc | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index e79bdd5a8..a0ee877f8 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -272,7 +272,9 @@ void ascii_to_hdf5_impl(const std::string& ascii_file, // Define compression filters DataSetCreateProps props; - props.add(Chunking(std::vector{(hsize_t)std::max(num_particles / 10, 1024)})); + hsize_t chunk = static_cast(std::max(num_particles / 10, 1024)); + if (chunk > static_cast(num_particles)) chunk = static_cast(num_particles); + props.add(Chunking(std::vector{chunk})); props.add(Shuffle()); props.add(Deflate(4)); // This is a good compromise between // speed and compression ratio From d19a40a67d7851ed970d2026655102d43255e653 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 4 Jul 2026 17:09:15 -0400 Subject: [PATCH 13/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- src/Component.cc | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/Component.cc b/src/Component.cc index 0e888cdd0..da38e43d5 100644 --- a/src/Component.cc +++ b/src/Component.cc @@ -1723,6 +1723,14 @@ void Component::read_bodies_and_distribute_hdf5(void) part->vel[0] = u[icount]; part->vel[1] = v[icount]; part->vel[2] = w[icount]; + + for (int j=0; jiattrib[j] = aux_ints[j][icount]; + } + + for (int j=0; jdattrib[j] = aux_floats[j][icount]; + } r2 = 0.0; for (int k=0; k<3; k++) r2 += part->pos[k]*part->pos[k]; From c54dda69e2e6185703d40e61369f89a63cc00e8d Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 4 Jul 2026 17:09:44 -0400 Subject: [PATCH 14/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/PhaseSpace/hdf5bods.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index a0ee877f8..f73f62376 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -60,9 +60,9 @@ File structure in HDF5 ├── aux_int_0, aux_int_1, ... (Datasets) └── aux_float_0, aux_float_1, ... (Datasets) -This schemea could (should be?) changed to use the Gadget style. On -the other hand, this design is clean and specific to exp. Also, users -can easily wrie and read this with (e.g.) h5py. I'm incliened to keep +This schema could (should be?) changed to use the Gadget style. On +the other hand, this design is clean and specific to EXP. Also, users +can easily write and read this with (e.g.) h5py. I'm inclined to keep it. */ From 608b35df9ff0dbd96cb1cc11902b9b3c313f5e0e Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 4 Jul 2026 18:21:36 -0400 Subject: [PATCH 15/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/PhaseSpace/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/utils/PhaseSpace/CMakeLists.txt b/utils/PhaseSpace/CMakeLists.txt index b951f1cd9..0433b0158 100644 --- a/utils/PhaseSpace/CMakeLists.txt +++ b/utils/PhaseSpace/CMakeLists.txt @@ -71,6 +71,7 @@ add_executable(pspmono pspmono.cc MakeModel.cc PSP.cc) add_executable(psp2hdf5 psp2hdf5.cc PSP.cc) add_executable(psporbv psporbv.cc PSP.cc) add_executable(hdf5bods hdf5bods.cc) +target_include_directories(hdf5bods PUBLIC ${HighFive_SOURCE_DIR}/include) if(HAVE_VTK) From 17f9eac3cf6f887a08e3b696ef7b12f6ead42e92 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 4 Jul 2026 18:21:57 -0400 Subject: [PATCH 16/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/PhaseSpace/hdf5bods.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index f73f62376..2003b4d1b 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -490,7 +490,7 @@ int main(int argc, char* argv[]) options.add_options() ("i,input", "Input prefix", cxxopts::value(prefix)->default_value("particles")) - ("o,output", "Output prefix (optional, otherwise input prefix is used)", cxxopts::value(prefix)) + ("o,output", "Output prefix (optional, otherwise input prefix is used)", cxxopts::value(output)) ("a,suffix", "Input suffix", cxxopts::value(suffix)->default_value("bods")) ("t,threads", "Number of OpenMP threads (default: max available)", cxxopts::value(num_threads)->default_value(std::to_string(num_threads))) ("roundtrip", "Perform round-trip conversion (ASCII -> HDF5 -> ASCII)") From e5e338acc710cf03210dff8de7710d714a3a6d73 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 4 Jul 2026 18:22:15 -0400 Subject: [PATCH 17/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/PhaseSpace/hdf5bods.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index 2003b4d1b..03a3697e5 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -1,4 +1,4 @@ -/* +/* Draft EXP ascii particle phase-space data HDF5 converter This is a possible template for implementing HDF5 input files for From f5a15367b28e88cee20158db3a2c5450333fc7c7 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 4 Jul 2026 18:22:39 -0400 Subject: [PATCH 18/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/PhaseSpace/hdf5bods.cc | 3 +++ 1 file changed, 3 insertions(+) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index 03a3697e5..f2ef13224 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -462,6 +462,9 @@ void hdf5_to_ascii(const std::string& hdf5_file, const std::string& ascii_file, // Write all lines sequentially to avoid I/O contention outfile.open(ascii_file, std::ios::app); + if (!outfile.is_open()) { + throw std::runtime_error("Could not open ASCII file for appending: " + ascii_file); + } for (const auto& line : output_lines) { outfile << line << "\n"; } From 038230099405015b79e7214533a0b6d1a47a567e Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 4 Jul 2026 18:22:58 -0400 Subject: [PATCH 19/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/PhaseSpace/hdf5bods.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index f2ef13224..0b2ae760d 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -520,7 +520,7 @@ int main(int argc, char* argv[]) std::cout << R"( Examples: - Convert a standard EXP input body file nameed 'mybods.bods' to hdf5 format + Convert a standard EXP input body file named 'mybods.bods' to HDF5 format $ hdf5bods --to_hdf5 -i mybods The resulting HDF5 file using float32 internally will be called 'mybods.h5' From 6fec121e4d3bbddcf46dbb95542c867e6f07c62d Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 4 Jul 2026 18:23:11 -0400 Subject: [PATCH 20/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/PhaseSpace/hdf5bods.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index 0b2ae760d..b05ea8a35 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -532,7 +532,7 @@ int main(int argc, char* argv[]) $ hdf5bods --to_ascii -i input The suffix on input can be customized. For example, the following converts - a standard EXP input body file nameed 'mybods.asc' to hdf5 format + a standard EXP input body file named 'mybods.asc' to HDF5 format $ hdf5bods --to_hdf5 --suffix=asc -i mybods The resulting HDF5 file using float32 format will be called 'mybods.h5' From 702a03125b2049cdc7cbfd71ec727284bea3f88d Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 4 Jul 2026 19:06:40 -0400 Subject: [PATCH 21/50] Updates for automatic datatype detection --- src/Component.H | 3 ++ src/Component.cc | 49 ++++++++++++++++++---------- utils/PhaseSpace/hdf5bods.cc | 62 +++++++++++++++++++++++------------- 3 files changed, 75 insertions(+), 39 deletions(-) diff --git a/src/Component.H b/src/Component.H index 7adbff7ae..f497a0062 100644 --- a/src/Component.H +++ b/src/Component.H @@ -369,6 +369,9 @@ protected: int H5chunk = 4096; //@} + //! Detect float precision by inspecting dataset type + size_t detect_precision(const HighFive::DataSet& dataset); + public: //! Describe the phase space coordinates diff --git a/src/Component.cc b/src/Component.cc index da38e43d5..17adf2a12 100644 --- a/src/Component.cc +++ b/src/Component.cc @@ -1498,32 +1498,45 @@ void Component::read_bodies_and_distribute_init() } } +// Detect float precision by inspecting dataset type +size_t Component::detect_precision(const HighFive::DataSet& dataset) +{ + auto datatype = dataset.getDataType(); + auto class_type = datatype.getClass(); + + // Get the size in bytes + size_t type_size = datatype.getSize(); + + // Float32 is 4 bytes, Float64 is 8 bytes + if (type_size == 4) { + return type_size; + } else if (type_size == 8) { + return type_size; + } else { + throw std::runtime_error("Unexpected float type size: " + std::to_string(type_size)); + } +} + void Component::read_bodies_and_distribute_hdf5(void) { HighFive::File file(pfile, HighFive::File::ReadOnly); - // Read metadata (all ranks read the same metadata, but only rank 0 reads the datasets) + // Read metadata (all ranks read the same metadata, but only rank 0 + // reads the datasets) // nbodies_tot = file.getAttribute("num_particles") .read(); niattrib = file.getAttribute("num_aux_ints") .read(); ndattrib = file.getAttribute("num_aux_floats").read(); - // Get precision flag for float datasets (0=float32, 1=float64) - // - int precision_flag = 1; // Default to float64 - try { - precision_flag = file.getAttribute("float_precision").read(); - } catch (...) { - std::ostringstream msg; - msg << "Component ERROR: 'float_precision' attribute not found in HDF5 file <" - << pfile << ">. Please ensure the file has this attribute set to 0 (float32) or 1 (float64)."; - throw std::runtime_error(msg.str()); - } - // Open particles group // HighFive::Group particles_group = file.getGroup("particles"); + // Detect precision by inspecting first dataset "m" + // + HighFive::DataSet m_dataset = particles_group.getDataSet("m"); + size_t precision = detect_precision(m_dataset); + // Variant for flexible reading and access of float or double datasets // using FloatDataVariant = std::variant, std::vector>; @@ -1558,7 +1571,7 @@ void Component::read_bodies_and_distribute_hdf5(void) // Lambda to read float or double based on precision flag // - auto read_dataset = [&particles_group, precision_flag] + auto read_dataset = [&particles_group, precision] (const std::string& name, size_t offset, size_t batch) -> FloatData { HighFive::DataSet dataset = particles_group.getDataSet(name); @@ -1573,15 +1586,17 @@ void Component::read_bodies_and_distribute_hdf5(void) // Sanity check to ensure we don't read beyond the dataset size_t current = std::min(batch, total - offset); // Wrap returned vector in FloatData - if (precision_flag == 0) { + if (precision == 4) { return FloatData(dataset.select({offset}, {current}).read>()); } - else { + else if (precision == 8) { return FloatData(dataset.select({offset}, {current}).read>()); } + else { + throw std::runtime_error("Unsupported precision detected in dataset: " + name); + } }; - // Okay, now we can read the datasets and distribute particles // across MPI rank, following the previously established logic. The // following code reads the datasets in batches and distributes them diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index b05ea8a35..f2ca7c314 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -60,10 +60,11 @@ File structure in HDF5 ├── aux_int_0, aux_int_1, ... (Datasets) └── aux_float_0, aux_float_1, ... (Datasets) -This schema could (should be?) changed to use the Gadget style. On -the other hand, this design is clean and specific to EXP. Also, users -can easily write and read this with (e.g.) h5py. I'm inclined to keep -it. +This schema is clean, to the point, and specific to EXP while vaguely +echoing the Gadget-style. Also, users can easily write and read this with +h5py and we can easily update our EXP-specific IC generators to use this +style by porting over the ascii_to_hdf5 routine(). + */ // C++ std @@ -262,11 +263,6 @@ void ascii_to_hdf5_impl(const std::string& ascii_file, file.createAttribute("num_aux_floats", DataSpace::From(num_aux_floats)) .write(num_aux_floats); - // Store precision type: 0 = float32, 1 = float64 - int precision_flag = (precision == FloatPrecision::FLOAT32) ? 0 : 1; - file.createAttribute("float_precision", DataSpace::From(precision_flag)) - .write(precision_flag); - // Create a group for particle data Group particles_group = file.createGroup("particles"); @@ -335,6 +331,26 @@ struct ParticleDataVariant FloatPrecision precision; }; +// Detect float precision by inspecting dataset type +FloatPrecision detect_precision(const DataSet& dataset) +{ + auto datatype = dataset.getDataType(); + auto class_type = datatype.getClass(); + + // Get the size in bytes + size_t type_size = datatype.getSize(); + + // Float32 is 4 bytes, Float64 is 8 bytes + if (type_size == 4) { + return FloatPrecision::FLOAT32; + } else if (type_size == 8) { + return FloatPrecision::FLOAT64; + } else { + throw std::runtime_error("Unexpected float type size: " + std::to_string(type_size)); + } +} + + // Read HDF5 with automatic precision detection ParticleDataVariant read_hdf5_data(const std::string& hdf5_file) { @@ -345,26 +361,22 @@ ParticleDataVariant read_hdf5_data(const std::string& hdf5_file) int num_aux_ints = file.getAttribute("num_aux_ints").read(); int num_aux_floats = file.getAttribute("num_aux_floats").read(); - // Auto-detect precision - int precision_flag = 1; // Default to float64 - try { - precision_flag = file.getAttribute("float_precision").read(); - } catch (...) { - // If attribute doesn't exist, assume float64 (backward compatibility) - } - - FloatPrecision precision = (precision_flag == 0) ? FloatPrecision::FLOAT32 - : FloatPrecision::FLOAT64; - // Open particles group Group particles_group = file.getGroup("particles"); + // Detect precision by inspecting first dataset "m" + DataSet m_dataset = particles_group.getDataSet("m"); + FloatPrecision precision = detect_precision(m_dataset); + ParticleDataVariant data; data.num_particles = num_particles; data.num_aux_ints = num_aux_ints; data.num_aux_floats = num_aux_floats; data.precision = precision; + std::string precision_str = (precision == FloatPrecision::FLOAT32) ? "float32" : "float64"; + std::cout << "Detected precision: " << precision_str << std::endl; + // Lambda to read float or double based on precision auto read_dataset = [&particles_group, precision](const std::string& name) -> FloatData { if (precision == FloatPrecision::FLOAT32) { @@ -489,7 +501,8 @@ int main(int argc, char* argv[]) // Parse command-line arguments for input/output files and mode // - cxxopts::Options options(argv[0], "ASCII to HDF5 particle converter with round-trip support and precision handling"); + cxxopts::Options options(argv[0], "ASCII \u2192 HDF5 and HDF5 \u2192 ASCII particle converter for EXP body files\n" + "with built-in round-trip testing and float-size selection\n"); options.add_options() ("i,input", "Input prefix", cxxopts::value(prefix)->default_value("particles")) @@ -516,8 +529,13 @@ int main(int argc, char* argv[]) if (vm.count("help")) { std::cout << options.help() << std::endl; - // Append custom examples + // Append note and custom examples std::cout << R"( +The best performance is achieved using the default 'float' rather than 'double' +precision for the HDF5 file, and that should be sufficient for initial data in +practice. However, the round-trip conversion will be exact only if double +precision is used. + Examples: Convert a standard EXP input body file named 'mybods.bods' to HDF5 format From 9b98b997fea2ac3548ad52d6de60cf2cf7403999 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 4 Jul 2026 19:25:59 -0400 Subject: [PATCH 22/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- src/Component.cc | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/Component.cc b/src/Component.cc index 17adf2a12..5df55a68d 100644 --- a/src/Component.cc +++ b/src/Component.cc @@ -1479,11 +1479,17 @@ void Component::read_bodies_and_distribute_init() // Note: One might need to use H5Fis_hdf5(filename.c_str()) for // older HDF5 versions (< 1.12) + bool isHDF5 = false; + + if (myid == 0) { #if H5_VERSION_GE(1, 12, 0) - bool isHDF5 = H5Fis_accessible(pfile.c_str(), H5P_DEFAULT) > 0; + isHDF5 = H5Fis_accessible(pfile.c_str(), H5P_DEFAULT) > 0; #else - bool isHDF5 = H5Fis_hdf5(pfile.c_str()) > 0; + isHDF5 = H5Fis_hdf5(pfile.c_str()) > 0; #endif + } + + MPI_Bcast(&isHDF5, 1, MPI_CXX_BOOL, 0, MPI_COMM_WORLD); if (isHDF5) { if (myid==0) From 00fd85fa53349bc8eb90f379afe86cfaef7193dc Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 4 Jul 2026 19:27:06 -0400 Subject: [PATCH 23/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/PhaseSpace/hdf5bods.cc | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index f2ca7c314..9b56880ce 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -334,20 +334,22 @@ struct ParticleDataVariant // Detect float precision by inspecting dataset type FloatPrecision detect_precision(const DataSet& dataset) { - auto datatype = dataset.getDataType(); - auto class_type = datatype.getClass(); + const auto datatype = dataset.getDataType(); + if (datatype.getClass() != DataTypeClass::Float) { + throw std::runtime_error("Unexpected dataset type for precision detection (expected float)"); + } // Get the size in bytes - size_t type_size = datatype.getSize(); + const size_t type_size = datatype.getSize(); // Float32 is 4 bytes, Float64 is 8 bytes if (type_size == 4) { return FloatPrecision::FLOAT32; } else if (type_size == 8) { return FloatPrecision::FLOAT64; - } else { - throw std::runtime_error("Unexpected float type size: " + std::to_string(type_size)); } + + throw std::runtime_error("Unexpected float type size: " + std::to_string(type_size)); } From 2fd10ed415d2703dc45bf0f7a780206916ec9556 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 4 Jul 2026 19:27:40 -0400 Subject: [PATCH 24/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/PhaseSpace/hdf5bods.cc | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index 9b56880ce..14b8bf2f9 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -376,8 +376,7 @@ ParticleDataVariant read_hdf5_data(const std::string& hdf5_file) data.num_aux_floats = num_aux_floats; data.precision = precision; - std::string precision_str = (precision == FloatPrecision::FLOAT32) ? "float32" : "float64"; - std::cout << "Detected precision: " << precision_str << std::endl; + // Precision is available in data.precision; caller decides whether to print it. // Lambda to read float or double based on precision auto read_dataset = [&particles_group, precision](const std::string& name) -> FloatData { From 42e4f25d1701709336581cd41ce0e21ccbccedee Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 4 Jul 2026 23:42:28 +0000 Subject: [PATCH 25/50] Remove unused mismatch flag in hdf5bods verification --- utils/PhaseSpace/hdf5bods.cc | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index 14b8bf2f9..56716b2a2 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -628,8 +628,7 @@ precision is used. std::string line_orig, line_rest; size_t line_num = 0; - bool mismatch_found = false; - + std::vector max_diff(7, 0.0); // For m, x, y, z, u, v, w std::vector vec_orig(7), vec_rest(7); @@ -651,7 +650,6 @@ precision is used. std::cerr << "Header mismatch at line 1:\n" << "Original: " << line_orig << "\n" << "Restored: " << line_rest << "\n"; - mismatch_found = true; break; } From 0930ad628cabf6e4b1fd9354b5a6e9d5583b3cc8 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 4 Jul 2026 20:39:01 -0400 Subject: [PATCH 26/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- src/Component.cc | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/Component.cc b/src/Component.cc index 5df55a68d..e45b77a1d 100644 --- a/src/Component.cc +++ b/src/Component.cc @@ -1507,20 +1507,20 @@ void Component::read_bodies_and_distribute_init() // Detect float precision by inspecting dataset type size_t Component::detect_precision(const HighFive::DataSet& dataset) { - auto datatype = dataset.getDataType(); - auto class_type = datatype.getClass(); + const auto datatype = dataset.getDataType(); + if (datatype.getClass() != HighFive::DataTypeClass::Float) { + throw std::runtime_error("Unexpected dataset type: expected float for precision detection"); + } // Get the size in bytes - size_t type_size = datatype.getSize(); + const size_t type_size = datatype.getSize(); // Float32 is 4 bytes, Float64 is 8 bytes - if (type_size == 4) { - return type_size; - } else if (type_size == 8) { + if (type_size == 4 || type_size == 8) { return type_size; - } else { - throw std::runtime_error("Unexpected float type size: " + std::to_string(type_size)); } + + throw std::runtime_error("Unexpected float type size: " + std::to_string(type_size)); } void Component::read_bodies_and_distribute_hdf5(void) From 35a1abd8d84bcecd1745e038effd4704798a3038 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 5 Jul 2026 00:46:18 +0000 Subject: [PATCH 27/50] Verify aux fields in hdf5bods roundtrip check --- utils/PhaseSpace/hdf5bods.cc | 66 ++++++++++++++++++++++++++++++++++-- 1 file changed, 63 insertions(+), 3 deletions(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index 56716b2a2..112094b12 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -74,6 +74,7 @@ style by porting over the ascii_to_hdf5 routine(). #include #include #include +#include #include #include #include @@ -631,6 +632,11 @@ precision is used. std::vector max_diff(7, 0.0); // For m, x, y, z, u, v, w std::vector vec_orig(7), vec_rest(7); + std::vector max_aux_int_diff; + std::vector max_aux_float_diff; + std::vector aux_int_orig, aux_int_rest; + std::vector aux_float_orig, aux_float_rest; + int num_aux_ints = -1, num_aux_floats = -1; while (std::getline(original, line_orig) && std::getline(restored, line_rest)) { ++line_num; @@ -652,13 +658,23 @@ precision is used. << "Restored: " << line_rest << "\n"; break; } + + num_aux_ints = total_aux_ints_orig; + num_aux_floats = total_aux_floats_orig; + max_aux_int_diff.assign(num_aux_ints, 0.0); + max_aux_float_diff.assign(num_aux_floats, 0.0); + aux_int_orig.assign(num_aux_ints, 0); + aux_int_rest.assign(num_aux_ints, 0); + aux_float_orig.assign(num_aux_floats, 0.0); + aux_float_rest.assign(num_aux_floats, 0.0); continue; } for (int i = 0; i < 7; ++i) { - ss_orig >> vec_orig[i]; - ss_rest >> vec_rest[i]; + if (!(ss_orig >> vec_orig[i]) || !(ss_rest >> vec_rest[i])) { + throw std::runtime_error("Failed to parse core field at line " + std::to_string(line_num)); + } } for (int i = 0; i < 7; ++i) { @@ -667,6 +683,29 @@ precision is used. } } + + for (int i = 0; i < num_aux_ints; ++i) { + if (!(ss_orig >> aux_int_orig[i]) || !(ss_rest >> aux_int_rest[i])) { + throw std::runtime_error("Failed to parse auxiliary integer field at line " + std::to_string(line_num)); + } + } + + for (int i = 0; i < num_aux_ints; ++i) { + double diff = std::abs(static_cast(aux_int_orig[i]) - + static_cast(aux_int_rest[i])); + if (diff > max_aux_int_diff[i]) max_aux_int_diff[i] = diff; + } + + for (int i = 0; i < num_aux_floats; ++i) { + if (!(ss_orig >> aux_float_orig[i]) || !(ss_rest >> aux_float_rest[i])) { + throw std::runtime_error("Failed to parse auxiliary float field at line " + std::to_string(line_num)); + } + } + + for (int i = 0; i < num_aux_floats; ++i) { + double diff = std::abs(aux_float_orig[i] - aux_float_rest[i]); + if (diff > max_aux_float_diff[i]) max_aux_float_diff[i] = diff; + } } std::cout << "Maximum absolute differences for core fields" @@ -678,6 +717,28 @@ precision is used. << max_diff[4] << ", " << max_diff[5] << ", " << max_diff[6] << std::endl; + + std::cout << "Maximum absolute differences for auxiliary integer fields: "; + if (max_aux_int_diff.empty()) { + std::cout << "(none)"; + } else { + for (size_t i = 0; i < max_aux_int_diff.size(); ++i) { + if (i) std::cout << ", "; + std::cout << max_aux_int_diff[i]; + } + } + std::cout << std::endl; + + std::cout << "Maximum absolute differences for auxiliary float fields: "; + if (max_aux_float_diff.empty()) { + std::cout << "(none)"; + } else { + for (size_t i = 0; i < max_aux_float_diff.size(); ++i) { + if (i) std::cout << ", "; + std::cout << max_aux_float_diff[i]; + } + } + std::cout << std::endl; } } catch (const std::exception& e) { @@ -745,4 +806,3 @@ precision is used. return 0; } - From 2f3d782210aded136b3b874a8379ab2431a8d850 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 5 Jul 2026 00:49:43 +0000 Subject: [PATCH 28/50] Use integer-safe aux-int diff in verify mode --- utils/PhaseSpace/hdf5bods.cc | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index 112094b12..32837863d 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -71,6 +71,7 @@ style by porting over the ascii_to_hdf5 routine(). #include #include #include +#include #include #include #include @@ -691,8 +692,9 @@ precision is used. } for (int i = 0; i < num_aux_ints; ++i) { - double diff = std::abs(static_cast(aux_int_orig[i]) - - static_cast(aux_int_rest[i])); + auto diff_i = std::llabs(static_cast(aux_int_orig[i]) - + static_cast(aux_int_rest[i])); + double diff = static_cast(diff_i); if (diff > max_aux_int_diff[i]) max_aux_int_diff[i] = diff; } From a486ae0e7ff65310dd044874f083eee227cddf59 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 4 Jul 2026 20:58:22 -0400 Subject: [PATCH 29/50] White-space changes only for readability --- utils/PhaseSpace/hdf5bods.cc | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index 32837863d..faf2d5de1 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -1,5 +1,5 @@ /* -Draft EXP ascii particle phase-space data HDF5 converter +EXP ascii particle phase-space data HDF5 converter This is a possible template for implementing HDF5 input files for EXP. The code reads an ASCII file with particle data and writes it to @@ -128,9 +128,10 @@ void parse_particle_line(const std::string& line, T val; // Parse core PSP fields - if (!(iss >> data.m[particle_idx] - >> data.x[particle_idx] >> data.y[particle_idx] >> data.z[particle_idx] - >> data.u[particle_idx] >> data.v[particle_idx] >> data.w[particle_idx])) { + if (!(iss + >> data.m[particle_idx] + >> data.x[particle_idx] >> data.y[particle_idx] >> data.z[particle_idx] + >> data.u[particle_idx] >> data.v[particle_idx] >> data.w[particle_idx])) { throw std::runtime_error("Failed to parse core fields at particle " + std::to_string(particle_idx)); } @@ -271,7 +272,8 @@ void ascii_to_hdf5_impl(const std::string& ascii_file, // Define compression filters DataSetCreateProps props; hsize_t chunk = static_cast(std::max(num_particles / 10, 1024)); - if (chunk > static_cast(num_particles)) chunk = static_cast(num_particles); + if (chunk > static_cast(num_particles)) + chunk = static_cast(num_particles); props.add(Chunking(std::vector{chunk})); props.add(Shuffle()); props.add(Deflate(4)); // This is a good compromise between @@ -391,9 +393,11 @@ ParticleDataVariant read_hdf5_data(const std::string& hdf5_file) // Read core PSP fields data.m = read_dataset("m"); + data.x = read_dataset("x"); data.y = read_dataset("y"); data.z = read_dataset("z"); + data.u = read_dataset("u"); data.v = read_dataset("v"); data.w = read_dataset("w"); From b8e26d45d71b08d6a93331c9db9e98e12f3b5494 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 5 Jul 2026 11:32:27 -0400 Subject: [PATCH 30/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/PhaseSpace/hdf5bods.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index faf2d5de1..c3a1a73fd 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -804,7 +804,7 @@ precision is used. if (!quiet) time_ptr = std::make_unique("conversion to ascii"); std::string ifile = prefix + ".h5"; std::string ofile = prefix + "." + suffix; - hdf5_to_ascii(ifile, ofile); + hdf5_to_ascii(ifile, ofile, !quiet); } else { std::cerr << "No conversion mode specified. Use --to_hdf5 to convert to HDF5, --to_ascii to convert from HDF5 to ascii, or the --roundtrip test." << std::endl; return 1; From e6db8e38389f8196a53ac5a0810b72ad2a05a2c9 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 5 Jul 2026 11:32:55 -0400 Subject: [PATCH 31/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/PhaseSpace/hdf5bods.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index c3a1a73fd..08b0c1a9f 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -502,7 +502,6 @@ int main(int argc, char* argv[]) { std::string prefix = "particles", output; std::string suffix = "bods"; - std::string ascii_restored = "particles_restored.txt"; int num_threads = omp_get_max_threads(); bool quiet = false; From 2700ef8933ac397c152abd43079e6f7c2b9942a1 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 5 Jul 2026 15:38:17 +0000 Subject: [PATCH 32/50] Fix m-field comment to mass --- utils/PhaseSpace/hdf5bods.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index 08b0c1a9f..3505770b9 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -105,7 +105,7 @@ enum class FloatPrecision { FLOAT32, FLOAT64 }; template struct ParticleDataTemplate { - std::vector m; // mass/species index + std::vector m; // mass std::vector x, y, z; // position std::vector u, v, w; // velocity std::vector> aux_ints; // auxiliary integer fields From 35b37852a74f23e6fb90c083d12d37f86a35455b Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 5 Jul 2026 11:59:33 -0400 Subject: [PATCH 33/50] Added compression filter selection for testing --- utils/PhaseSpace/hdf5bods.cc | 173 +++++++++++++++++++++++++++++------ 1 file changed, 147 insertions(+), 26 deletions(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index 08b0c1a9f..31f7e188d 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -90,6 +90,9 @@ style by porting over the ascii_to_hdf5 routine(). #include #include +// Required for H5Zfilter_avail +#include + // Command-line parsing #include "cxxopts.H" @@ -98,6 +101,82 @@ style by porting over the ascii_to_hdf5 routine(). using namespace HighFive; + +// Factory designed to iterate across filter list available on Ubuntu +// 26.04 LTS and other systems with HDF5 1.14.x and HighFive 2.6.x +HighFive::DataSetCreateProps +createFilterProps(const std::vector& chunk_dims, + unsigned int filter_id, + // Scale from 1 (fast) to 9 (max compression) + int compression_level = 5) +{ + HighFive::DataSetCreateProps props; + props.add(HighFive::Chunking(chunk_dims)); + + // Best Practice: Always shuffle float/double streams BEFORE compressing + // (Note: Blosc handles shuffling internally via parameters) + if (filter_id != 3 && filter_id != 32001) { + props.add(HighFive::Shuffle()); + } + + hid_t plist_id = props.getId(); + std::vector cd_values; + + switch (filter_id) { + case 1: // Deflate / GZIP (Built-in) + // Expects 1 parameter: compression level (0-9) + cd_values = { static_cast(compression_level) }; + break; + + case 2: // SZIP (Built-in) + // SZIP is complex and highly variant; HighFive includes native wrappers. + // For raw testing, standard pixels-per-block option is typically 16. + cd_values = { 141, 16 }; + break; + + case 3: // Shuffle Filter alone (No trailing compression) + props.add(HighFive::Shuffle()); + return props; + + case 4: // Fletcher32 Checksum alone (Data validation, not compression) + // Fletcher32 uses built-in ID 4 and accepts 0 configuration arguments. + // Call the native HDF5 library directly using HighFive's internal handle. + H5Pset_filter(plist_id, 4, H5Z_FLAG_OPTIONAL, 0, NULL); + return props; + + case 307: // BZip2 + // Expects 1 parameter: Block size in 100KB steps (1-9) + cd_values = { static_cast(compression_level) }; + break; + + case 32004: // LZ4 + // Expects 1 parameter: internal chunk/block size. + // 0 falls back to the default (64KB), optimized for CPU caches. + cd_values = { 0 }; + break; + + case 32001: // Blosc (v1 Meta-Compressor) + // Expects 7 parameters. Slots 0-3 are reserved. + // [4]=level(1-9), [5]=shuffle type(1=byte, 2=bit), [6]=codec code + // Codecs: 0=blosclz, 1=lz4, 2=lz4hc, 3=snappy, 4=zlib, 5=zstd + cd_values = { 0, 0, 0, 0, + static_cast(compression_level), + 1, // 1 = Byte Shuffle (Highly recommended for float/double) + 1 }; // 1 = Redirect internal processing to LZ4 + break; + + default: + throw std::invalid_argument("Unknown or unsupported testing filter ID requested."); + } + + // Pass the custom properties down to the HDF5 backend pipeline + if (!cd_values.empty()) { + H5Pset_filter(plist_id, filter_id, H5Z_FLAG_OPTIONAL, cd_values.size(), cd_values.data()); + } + + return props; +} + // Specify floating point precision enum class FloatPrecision { FLOAT32, FLOAT64 }; @@ -197,6 +276,7 @@ template void ascii_to_hdf5_impl(const std::string& ascii_file, const std::string& hdf5_file, FloatPrecision precision, + unsigned filter_id, bool verbose) { // Header values @@ -267,17 +347,22 @@ void ascii_to_hdf5_impl(const std::string& ascii_file, .write(num_aux_floats); // Create a group for particle data + // Group particles_group = file.createGroup("particles"); - // Define compression filters - DataSetCreateProps props; - hsize_t chunk = static_cast(std::max(num_particles / 10, 1024)); - if (chunk > static_cast(num_particles)) - chunk = static_cast(num_particles); - props.add(Chunking(std::vector{chunk})); - props.add(Shuffle()); - props.add(Deflate(4)); // This is a good compromise between - // speed and compression ratio + // Calculate optimized, capped chunk size + // + hsize_t chunk_size = 262144; // Cachesafe cap (~2MB for double) + if (num_particles < chunk_size) { + chunk_size = std::max(1024, num_particles); + } + if (chunk_size > num_particles) { + chunk_size = num_particles; // Absolute safety fallback + } + + // Define compression filter + // + auto props = createFilterProps({chunk_size}, filter_id, 5); // Write datasets sequentially (HDF5 is not fully thread-safe for writing) particles_group.createDataSet("m", data.m, props); @@ -311,12 +396,13 @@ void ascii_to_hdf5_impl(const std::string& ascii_file, void ascii_to_hdf5(const std::string& ascii_file, const std::string& hdf5_file, FloatPrecision precision = FloatPrecision::FLOAT64, + unsigned filter_id = 1, bool verbose = true) { if (precision == FloatPrecision::FLOAT32) { - ascii_to_hdf5_impl(ascii_file, hdf5_file, precision, verbose); + ascii_to_hdf5_impl(ascii_file, hdf5_file, precision, filter_id, verbose); } else { - ascii_to_hdf5_impl(ascii_file, hdf5_file, precision, verbose); + ascii_to_hdf5_impl(ascii_file, hdf5_file, precision, filter_id, verbose); } } @@ -503,6 +589,7 @@ int main(int argc, char* argv[]) std::string prefix = "particles", output; std::string suffix = "bods"; int num_threads = omp_get_max_threads(); + unsigned filter_id = 1; bool quiet = false; // Parse command-line arguments for input/output files and mode @@ -511,17 +598,19 @@ int main(int argc, char* argv[]) "with built-in round-trip testing and float-size selection\n"); options.add_options() - ("i,input", "Input prefix", cxxopts::value(prefix)->default_value("particles")) - ("o,output", "Output prefix (optional, otherwise input prefix is used)", cxxopts::value(output)) - ("a,suffix", "Input suffix", cxxopts::value(suffix)->default_value("bods")) - ("t,threads", "Number of OpenMP threads (default: max available)", cxxopts::value(num_threads)->default_value(std::to_string(num_threads))) - ("roundtrip", "Perform round-trip conversion (ASCII -> HDF5 -> ASCII)") - ("verify", "Verify that the restored ASCII file matches the original in 'roundtrip' mode") - ("to_hdf5", "Convert ASCII to HDF5") - ("double", "Use double precision for HDF5 output (float is default)") - ("to_ascii", "Convert HDF5 to ASCII") - ("q,quiet", "Suppress verbose output") - ("h,help", "Print usage"); + ("i,input", "Input prefix", cxxopts::value(prefix)->default_value("particles")) + ("o,output", "Output prefix (optional, otherwise input prefix is used)", cxxopts::value(output)) + ("a,suffix", "Input suffix", cxxopts::value(suffix)->default_value("bods")) + ("t,threads", "Number of OpenMP threads (default: max available)", cxxopts::value(num_threads)->default_value(std::to_string(num_threads))) + ("f,filter", "HDF5 filter ID to use (default: 1 = GZIP)", cxxopts::value(filter_id)->default_value("1")) + ("roundtrip", "Perform round-trip conversion (ASCII -> HDF5 -> ASCII)") + ("verify", "Verify that the restored ASCII file matches the original in 'roundtrip' mode") + ("to_hdf5", "Convert ASCII to HDF5") + ("double", "Use double precision for HDF5 output (float is default)") + ("filterlist", "List available HDF5 filters and exit") + ("to_ascii", "Convert HDF5 to ASCII") + ("q,quiet", "Suppress verbose output") + ("h,help", "Print usage"); cxxopts::ParseResult vm; @@ -569,6 +658,38 @@ precision is used. return 0; } + if (vm.count("filterlist")) { + + // Map of commonly registered official HDF5 Filter IDs + std::map known_filters = { + {1, "Deflate / GZIP (Built-in)"}, + {2, "SZIP (Built-in)"}, + {3, "Shuffle (Built-in)"}, + {4, "Fletcher32 Checksum (Built-in)"}, + {32001, "Blosc"}, + {32004, "LZ4"}, + {32008, "Bitshuffle"}, + {32013, "Zfp"}, + {32015, "Zstd"}, + {307, "BZip2"} + }; + + std::cout << "Checking available HDF5 filters on your system:\n"; + std::cout << "-----------------------------------------------\n"; + + for (const auto& [id, name] : known_filters) { + // H5Zfilter_avail returns a positive number if available, 0 if not + htri_t available = H5Zfilter_avail(id); + + if (available > 0) { + std::cout << "[AVAILABLE] ID: " << id << " -> " << name << "\n"; + } else { + std::cout << "[NOT FOUND] ID: " << id << " -> " << name << "\n"; + } + } + return 0; + } + if (vm.count("quiet")) { quiet = true; } @@ -601,7 +722,7 @@ precision is used. { std::unique_ptr time_ptr; if (!quiet) time_ptr = std::make_unique("float 32 conversion"); - ascii_to_hdf5(ascii, h5_32, FloatPrecision::FLOAT32, !quiet); + ascii_to_hdf5(ascii, h5_32, FloatPrecision::FLOAT32, filter_id, !quiet); } // Convert to HDF5 with float64 precision @@ -609,7 +730,7 @@ precision is used. { std::unique_ptr time_ptr; if (!quiet) time_ptr = std::make_unique("float 64 conversion"); - ascii_to_hdf5(ascii, h5_64, FloatPrecision::FLOAT64, !quiet); + ascii_to_hdf5(ascii, h5_64, FloatPrecision::FLOAT64, filter_id, !quiet); } // Round-trip with float64 version @@ -795,9 +916,9 @@ precision is used. std::string ofile = prefix + ".h5"; if (vm.count("output")) ofile = output + ".h5"; if (vm.count("double")) - ascii_to_hdf5(ascii, ofile, FloatPrecision::FLOAT64, !quiet); + ascii_to_hdf5(ascii, ofile, FloatPrecision::FLOAT64, filter_id, !quiet); else - ascii_to_hdf5(ascii, ofile, FloatPrecision::FLOAT32, !quiet); + ascii_to_hdf5(ascii, ofile, FloatPrecision::FLOAT32, filter_id, !quiet); } else if (vm.count("to_ascii")) { std::unique_ptr time_ptr; if (!quiet) time_ptr = std::make_unique("conversion to ascii"); From 06a3a5a034b89dfe5c1c35729682e100c092238a Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 5 Jul 2026 16:01:36 +0000 Subject: [PATCH 34/50] Add optional HDF5 particle index support --- src/Component.cc | 18 +++++++-- utils/PhaseSpace/hdf5bods.cc | 75 ++++++++++++++++++++++++++---------- 2 files changed, 68 insertions(+), 25 deletions(-) diff --git a/src/Component.cc b/src/Component.cc index e45b77a1d..ab305e00b 100644 --- a/src/Component.cc +++ b/src/Component.cc @@ -1572,6 +1572,8 @@ void Component::read_bodies_and_distribute_hdf5(void) // Phase-space data // FloatData m, x, y, z, u, v, w; + std::vector index; + bool has_index = particles_group.exist("index"); std::vector> aux_ints; std::vector aux_floats; @@ -1627,6 +1629,10 @@ void Component::read_bodies_and_distribute_hdf5(void) size_t offset = 0; // Read Node 0's PSP fields + if (has_index) { + index = particles_group.getDataSet("index") + .select({offset}, {nbodies_table[0]}).read>(); + } m = read_dataset("m", offset, nbodies_table[0]); x = read_dataset("x", offset, nbodies_table[0]); @@ -1675,8 +1681,8 @@ void Component::read_bodies_and_distribute_hdf5(void) part->dattrib[j] = aux_floats[j][i]; } - // Make a sequential index for the particle - part->indx = i + 1; + // Set particle index from optional dataset, otherwise sequential + part->indx = has_index ? index[i] : i + 1; // Get the radius r2 = 0.0; @@ -1699,6 +1705,10 @@ void Component::read_bodies_and_distribute_hdf5(void) for (int n=1; n>(); + } m = read_dataset("m", offset, nbodies_table[n]); x = read_dataset("x", offset, nbodies_table[n]); @@ -1732,8 +1742,8 @@ void Component::read_bodies_and_distribute_hdf5(void) PartPtr part = std::make_shared(niattrib, ndattrib); - // Make a sequential index for the particle - part->indx = nbodies_index[n-1] + 1 + icount; + // Set particle index from optional dataset, otherwise sequential + part->indx = has_index ? index[icount] : nbodies_index[n-1] + 1 + icount; part->mass = m[icount]; diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index ebefe5a9f..27776f643 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -57,6 +57,7 @@ File structure in HDF5 ├── m (Dataset) ├── x, y, z (Datasets) ├── u, v, w (Datasets) +├── index (optional Dataset) ├── aux_int_0, aux_int_1, ... (Datasets) └── aux_float_0, aux_float_1, ... (Datasets) @@ -80,6 +81,7 @@ style by porting over the ascii_to_hdf5 routine(). #include #include #include +#include #include // For OpenMP parallelization @@ -184,6 +186,7 @@ enum class FloatPrecision { FLOAT32, FLOAT64 }; template struct ParticleDataTemplate { + std::optional> index; // optional particle index std::vector m; // mass std::vector x, y, z; // position std::vector u, v, w; // velocity @@ -199,6 +202,7 @@ struct ParticleDataTemplate template void parse_particle_line(const std::string& line, int particle_idx, + bool has_index, int num_aux_ints, int num_aux_floats, ParticleDataTemplate& data) @@ -206,6 +210,16 @@ void parse_particle_line(const std::string& line, std::istringstream iss(line); T val; + // Parse optional index field + if (has_index) { + unsigned long idx = 0; + if (!(iss >> idx)) { + throw std::runtime_error("Failed to parse index at particle " + + std::to_string(particle_idx)); + } + (*data.index)[particle_idx] = idx; + } + // Parse core PSP fields if (!(iss >> data.m[particle_idx] @@ -276,6 +290,7 @@ template void ascii_to_hdf5_impl(const std::string& ascii_file, const std::string& hdf5_file, FloatPrecision precision, + bool has_index, unsigned filter_id, bool verbose) { @@ -301,6 +316,7 @@ void ascii_to_hdf5_impl(const std::string& ascii_file, data.u.resize(num_particles); data.v.resize(num_particles); data.w.resize(num_particles); + if (has_index) data.index = std::vector(num_particles); data.aux_ints.resize(num_aux_ints); for (auto& vec : data.aux_ints) { @@ -320,7 +336,7 @@ void ascii_to_hdf5_impl(const std::string& ascii_file, num_threads(omp_get_max_threads()) for (int i = 0; i < num_particles; ++i) { try { - parse_particle_line(lines[i], i, num_aux_ints, num_aux_floats, data); + parse_particle_line(lines[i], i, has_index, num_aux_ints, num_aux_floats, data); } catch (const std::exception& e) { parse_errors += 1; #pragma omp critical @@ -372,6 +388,9 @@ void ascii_to_hdf5_impl(const std::string& ascii_file, particles_group.createDataSet("u", data.u, props); particles_group.createDataSet("v", data.v, props); particles_group.createDataSet("w", data.w, props); + if (data.index.has_value()) { + particles_group.createDataSet("index", *data.index, props); + } // Write auxiliary integer fields for (size_t j = 0; j < data.aux_ints.size(); ++j) { @@ -396,13 +415,14 @@ void ascii_to_hdf5_impl(const std::string& ascii_file, void ascii_to_hdf5(const std::string& ascii_file, const std::string& hdf5_file, FloatPrecision precision = FloatPrecision::FLOAT64, + bool has_index = false, unsigned filter_id = 1, bool verbose = true) { if (precision == FloatPrecision::FLOAT32) { - ascii_to_hdf5_impl(ascii_file, hdf5_file, precision, filter_id, verbose); + ascii_to_hdf5_impl(ascii_file, hdf5_file, precision, has_index, filter_id, verbose); } else { - ascii_to_hdf5_impl(ascii_file, hdf5_file, precision, filter_id, verbose); + ascii_to_hdf5_impl(ascii_file, hdf5_file, precision, has_index, filter_id, verbose); } } @@ -411,6 +431,7 @@ using FloatData = std::variant, std::vector>; struct ParticleDataVariant { + std::optional> index; FloatData m, x, y, z, u, v, w; std::vector> aux_ints; std::vector aux_floats; @@ -478,6 +499,9 @@ ParticleDataVariant read_hdf5_data(const std::string& hdf5_file) }; // Read core PSP fields + if (particles_group.exist("index")) { + data.index = particles_group.getDataSet("index").read>(); + } data.m = read_dataset("m"); data.x = read_dataset("x"); @@ -543,7 +567,8 @@ void hdf5_to_ascii(const std::string& hdf5_file, const std::string& ascii_file, std::ostringstream oss; oss.precision(16); - // Write core phase-space fields + // Write optional index and core phase-space fields + if (data.index.has_value()) oss << (*data.index)[i] << " "; oss << write_value(data.m, i) << " " << write_value(data.x, i) << " " << write_value(data.y, i) << " " @@ -598,19 +623,20 @@ int main(int argc, char* argv[]) "with built-in round-trip testing and float-size selection\n"); options.add_options() - ("i,input", "Input prefix", cxxopts::value(prefix)->default_value("particles")) - ("o,output", "Output prefix (optional, otherwise input prefix is used)", cxxopts::value(output)) - ("a,suffix", "Input suffix", cxxopts::value(suffix)->default_value("bods")) - ("t,threads", "Number of OpenMP threads (default: max available)", cxxopts::value(num_threads)->default_value(std::to_string(num_threads))) - ("f,filter", "HDF5 filter ID to use (default: 1 = GZIP)", cxxopts::value(filter_id)->default_value("1")) - ("roundtrip", "Perform round-trip conversion (ASCII -> HDF5 -> ASCII)") - ("verify", "Verify that the restored ASCII file matches the original in 'roundtrip' mode") - ("to_hdf5", "Convert ASCII to HDF5") - ("double", "Use double precision for HDF5 output (float is default)") - ("filterlist", "List available HDF5 filters and exit") - ("to_ascii", "Convert HDF5 to ASCII") - ("q,quiet", "Suppress verbose output") - ("h,help", "Print usage"); + ("i,input", "Input prefix", cxxopts::value(prefix)->default_value("particles")) + ("o,output", "Output prefix (optional, otherwise input prefix is used)", cxxopts::value(output)) + ("a,suffix", "Input suffix", cxxopts::value(suffix)->default_value("bods")) + ("t,threads", "Number of OpenMP threads (default: max available)", cxxopts::value(num_threads)->default_value(std::to_string(num_threads))) + ("f,filter", "HDF5 filter ID to use (default: 1 = GZIP)", cxxopts::value(filter_id)->default_value("1")) + ("roundtrip", "Perform round-trip conversion (ASCII -> HDF5 -> ASCII)") + ("verify", "Verify that the restored ASCII file matches the original in 'roundtrip' mode") + ("to_hdf5", "Convert ASCII to HDF5") + ("double", "Use double precision for HDF5 output (float is default)") + ("filterlist","List available HDF5 filters and exit") + ("aindex", "Input ASCII has leading particle index column; write optional HDF5 particles/index dataset") + ("to_ascii", "Convert HDF5 to ASCII") + ("q,quiet", "Suppress verbose output") + ("h,help", "Print usage"); cxxopts::ParseResult vm; @@ -722,7 +748,7 @@ precision is used. { std::unique_ptr time_ptr; if (!quiet) time_ptr = std::make_unique("float 32 conversion"); - ascii_to_hdf5(ascii, h5_32, FloatPrecision::FLOAT32, filter_id, !quiet); + ascii_to_hdf5(ascii, h5_32, FloatPrecision::FLOAT32, vm.count("aindex"), filter_id, !quiet); } // Convert to HDF5 with float64 precision @@ -730,7 +756,7 @@ precision is used. { std::unique_ptr time_ptr; if (!quiet) time_ptr = std::make_unique("float 64 conversion"); - ascii_to_hdf5(ascii, h5_64, FloatPrecision::FLOAT64, filter_id, !quiet); + ascii_to_hdf5(ascii, h5_64, FloatPrecision::FLOAT64, vm.count("aindex"), filter_id, !quiet); } // Round-trip with float64 version @@ -796,6 +822,13 @@ precision is used. continue; } + if (vm.count("aindex")) { + unsigned long idx_orig = 0, idx_rest = 0; + if (!(ss_orig >> idx_orig) || !(ss_rest >> idx_rest)) { + throw std::runtime_error("Failed to parse index field at line " + std::to_string(line_num)); + } + } + for (int i = 0; i < 7; ++i) { if (!(ss_orig >> vec_orig[i]) || !(ss_rest >> vec_rest[i])) { throw std::runtime_error("Failed to parse core field at line " + std::to_string(line_num)); @@ -916,9 +949,9 @@ precision is used. std::string ofile = prefix + ".h5"; if (vm.count("output")) ofile = output + ".h5"; if (vm.count("double")) - ascii_to_hdf5(ascii, ofile, FloatPrecision::FLOAT64, filter_id, !quiet); + ascii_to_hdf5(ascii, ofile, FloatPrecision::FLOAT64, vm.count("aindex"), filter_id, !quiet); else - ascii_to_hdf5(ascii, ofile, FloatPrecision::FLOAT32, filter_id, !quiet); + ascii_to_hdf5(ascii, ofile, FloatPrecision::FLOAT32, vm.count("aindex"), filter_id, !quiet); } else if (vm.count("to_ascii")) { std::unique_ptr time_ptr; if (!quiet) time_ptr = std::make_unique("conversion to ascii"); From fac2cbe440645b1dfc1d42e783dc1ee2fb7a22f6 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 5 Jul 2026 12:15:50 -0400 Subject: [PATCH 35/50] Minor change to option documentation only --- utils/PhaseSpace/hdf5bods.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index 27776f643..1b04fa6bf 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -633,7 +633,7 @@ int main(int argc, char* argv[]) ("to_hdf5", "Convert ASCII to HDF5") ("double", "Use double precision for HDF5 output (float is default)") ("filterlist","List available HDF5 filters and exit") - ("aindex", "Input ASCII has leading particle index column; write optional HDF5 particles/index dataset") + ("aindex", "Input ASCII has leading particle index column; write the particle index dataset to HDF5") ("to_ascii", "Convert HDF5 to ASCII") ("q,quiet", "Suppress verbose output") ("h,help", "Print usage"); From b59afd6f4f0d028524925406b0e62e87fd92d7c1 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 5 Jul 2026 12:27:54 -0400 Subject: [PATCH 36/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/PhaseSpace/hdf5bods.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index 1b04fa6bf..d6dc1800b 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -78,6 +78,7 @@ style by porting over the ascii_to_hdf5 routine(). #include #include #include +#include #include #include #include From 50b30769867eb32bb2cf88467a605fe213117728 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 5 Jul 2026 12:35:24 -0400 Subject: [PATCH 37/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/PhaseSpace/hdf5bods.cc | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index d6dc1800b..2f2d6f906 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -116,10 +116,7 @@ createFilterProps(const std::vector& chunk_dims, HighFive::DataSetCreateProps props; props.add(HighFive::Chunking(chunk_dims)); - // Best Practice: Always shuffle float/double streams BEFORE compressing - // (Note: Blosc handles shuffling internally via parameters) - if (filter_id != 3 && filter_id != 32001) { - props.add(HighFive::Shuffle()); + if (filter_id != 3 && filter_id != 4 && filter_id != 32001) { } hid_t plist_id = props.getId(); From aa11000d5069190488f45527f8de8158636faf87 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 5 Jul 2026 12:36:14 -0400 Subject: [PATCH 38/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- src/Component.cc | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/Component.cc b/src/Component.cc index ab305e00b..68b9c4f34 100644 --- a/src/Component.cc +++ b/src/Component.cc @@ -1591,11 +1591,9 @@ void Component::read_bodies_and_distribute_hdf5(void) throw std::runtime_error("Component ERROR: dataset <" + name + "> shorter than expected"); } - // Sanity check to ensure we don't read beyond the dataset size_t current = std::min(batch, total - offset); - // Wrap returned vector in FloatData - if (precision == 4) { - return FloatData(dataset.select({offset}, {current}).read>()); + if (current != batch) { + throw std::runtime_error("Component ERROR: dataset <" + name + "> shorter than expected"); } else if (precision == 8) { return FloatData(dataset.select({offset}, {current}).read>()); From 3823cffe52f1ae9ed0d929753056f03233fa3e13 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 5 Jul 2026 12:50:42 -0400 Subject: [PATCH 39/50] Added some Doxygen documentation to the Component header [no ci] --- src/Component.H | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/src/Component.H b/src/Component.H index f497a0062..1ae7ba5f2 100644 --- a/src/Component.H +++ b/src/Component.H @@ -184,6 +184,43 @@ struct loadb_datum specifying alternative parameters, using an old-style PSP file, or starting a new job with a previous output from another simulation) + +
Input data format Component now supports reading from + HDF5 files on input. The scheme for HDF5 input files is: + @verbatim + / + ├── Attributes + │ ├── num_particles (int) + │ ├── num_aux_ints (int) + │ └── num_aux_floats (int) + └── particles/ (Group) + ├── m (Dataset) + ├── x, y, z (Datasets) + ├── u, v, w (Datasets) + ├── index (optional Dataset) + ├── aux_int_0, aux_int_1, ... (Datasets) + └── aux_float_0, aux_float_1, ... (Datasets) + @endverbatim + + @note +
    + +
  1. The index dataset is optional. If it is not + present, the particles are assigned indices in the order they are read. + +
  2. The aux_int_* and aux_float_* + datasets are optional. If they are not present, the corresponding + vectors are empty. The number of auxiliary integer and float + datasets is determined by the num_aux_ints and + num_aux_floats attributes, respectively. + +
  3. The datatype for the m, x, + y, z, u, v, + w and aux_float_* datasets maybe float or + double (i.e. float32 or float64) but they all must are assumed to + have the same datatype as m. + +
*/ class Component { From 1f9b79789dca7c7fa20ff72bf931d9c163319f0f Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 5 Jul 2026 16:42:36 -0400 Subject: [PATCH 40/50] Looks like garbled merged; fix propsed --- src/Component.cc | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/Component.cc b/src/Component.cc index 68b9c4f34..f2dd4ea87 100644 --- a/src/Component.cc +++ b/src/Component.cc @@ -1595,6 +1595,10 @@ void Component::read_bodies_and_distribute_hdf5(void) if (current != batch) { throw std::runtime_error("Component ERROR: dataset <" + name + "> shorter than expected"); } + + if (precision == 4) { + return FloatData(dataset.select({offset}, {current}).read>()); + } else if (precision == 8) { return FloatData(dataset.select({offset}, {current}).read>()); } From 3a090f4d8935012b2067ed275a23c758f928ebad Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 6 Jul 2026 08:44:59 -0400 Subject: [PATCH 41/50] Move precision detection into the dataset reader function --- src/Component.cc | 8 ++------ utils/PhaseSpace/hdf5bods.cc | 20 +++++++------------- 2 files changed, 9 insertions(+), 19 deletions(-) diff --git a/src/Component.cc b/src/Component.cc index f2dd4ea87..0a342a9a5 100644 --- a/src/Component.cc +++ b/src/Component.cc @@ -1538,11 +1538,6 @@ void Component::read_bodies_and_distribute_hdf5(void) // HighFive::Group particles_group = file.getGroup("particles"); - // Detect precision by inspecting first dataset "m" - // - HighFive::DataSet m_dataset = particles_group.getDataSet("m"); - size_t precision = detect_precision(m_dataset); - // Variant for flexible reading and access of float or double datasets // using FloatDataVariant = std::variant, std::vector>; @@ -1579,10 +1574,11 @@ void Component::read_bodies_and_distribute_hdf5(void) // Lambda to read float or double based on precision flag // - auto read_dataset = [&particles_group, precision] + auto read_dataset = [this, &particles_group] (const std::string& name, size_t offset, size_t batch) -> FloatData { HighFive::DataSet dataset = particles_group.getDataSet(name); + size_t precision = detect_precision(dataset); // Fetch the dimensions (e.g., [rows, cols]) size_t total = dataset.getSpace().getDimensions()[0]; diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index 2f2d6f906..c61fa7a5b 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -437,7 +437,6 @@ struct ParticleDataVariant size_t num_particles = 0; size_t num_aux_ints = 0; size_t num_aux_floats = 0; - FloatPrecision precision; }; // Detect float precision by inspecting dataset type @@ -475,20 +474,17 @@ ParticleDataVariant read_hdf5_data(const std::string& hdf5_file) // Open particles group Group particles_group = file.getGroup("particles"); - // Detect precision by inspecting first dataset "m" - DataSet m_dataset = particles_group.getDataSet("m"); - FloatPrecision precision = detect_precision(m_dataset); - ParticleDataVariant data; data.num_particles = num_particles; data.num_aux_ints = num_aux_ints; data.num_aux_floats = num_aux_floats; - data.precision = precision; - - // Precision is available in data.precision; caller decides whether to print it. // Lambda to read float or double based on precision - auto read_dataset = [&particles_group, precision](const std::string& name) -> FloatData { + auto read_dataset = [&particles_group](const std::string& name) -> FloatData { + // Detect precision by inspecting dataset + DataSet dataset = particles_group.getDataSet(name); + FloatPrecision precision = detect_precision(dataset); + if (precision == FloatPrecision::FLOAT32) { return particles_group.getDataSet(name).read>(); } else { @@ -598,11 +594,9 @@ void hdf5_to_ascii(const std::string& hdf5_file, const std::string& ascii_file, } outfile.close(); - std::string precision_str = (data.precision == FloatPrecision::FLOAT32) - ? "float32" : "float64"; if (verbose) - std::cout << "Successfully wrote " << data.num_particles << " particles to " - << ascii_file << " (" << precision_str << ")" << std::endl; + std::cout << "Successfully wrote " << data.num_particles + << " particles to " << ascii_file << std::endl; } // Main function with command-line parsing From 72c66c3d53814e2872472cf7b36ee9758ef5c6c1 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Mon, 6 Jul 2026 08:54:59 -0400 Subject: [PATCH 42/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- src/Component.cc | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/Component.cc b/src/Component.cc index 0a342a9a5..3b61983f3 100644 --- a/src/Component.cc +++ b/src/Component.cc @@ -1699,16 +1699,21 @@ void Component::read_bodies_and_distribute_hdf5(void) // Offset into the next node's data in the HDF5 datasets offset += nbodies_table[0]; - unsigned icount, ibufcount; for (int n=1; nShipParticles(n, 0, 0); + continue; + } + // Read core PSP fields if (has_index) { index = particles_group.getDataSet("index") .select({offset}, {nbodies_table[n]}).read>(); } m = read_dataset("m", offset, nbodies_table[n]); - x = read_dataset("x", offset, nbodies_table[n]); y = read_dataset("y", offset, nbodies_table[n]); z = read_dataset("z", offset, nbodies_table[n]); From 4df6f296bf7e7b5d17f626abcf1e15c41eaad9a3 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Mon, 6 Jul 2026 08:57:21 -0400 Subject: [PATCH 43/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/PhaseSpace/hdf5bods.cc | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index c61fa7a5b..967783384 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -116,9 +116,11 @@ createFilterProps(const std::vector& chunk_dims, HighFive::DataSetCreateProps props; props.add(HighFive::Chunking(chunk_dims)); + // For most compressors, applying shuffle improves compression for floating-point data. + // (Blosc can do its own shuffling internally; shuffle-only/checksum-only are handled below.) if (filter_id != 3 && filter_id != 4 && filter_id != 32001) { + props.add(HighFive::Shuffle()); } - hid_t plist_id = props.getId(); std::vector cd_values; From f18ed1f283be941bacb6994530efe9c568e4df5f Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 6 Jul 2026 12:57:43 +0000 Subject: [PATCH 44/50] docs: update Component.H to reflect per-dataset independent precision detection --- src/Component.H | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/Component.H b/src/Component.H index 1ae7ba5f2..18fdc4efe 100644 --- a/src/Component.H +++ b/src/Component.H @@ -216,9 +216,10 @@ struct loadb_datum
  • The datatype for the m, x, y, z, u, v, - w and aux_float_* datasets maybe float or - double (i.e. float32 or float64) but they all must are assumed to - have the same datatype as m. + w and aux_float_* datasets may be float or + double (i.e. float32 or float64). Each dataset's precision is + detected independently, so datasets may have mixed precision within + the same file. */ From b8e2546772fe0411bee394a145f09d4c62d50817 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 6 Jul 2026 12:38:25 -0400 Subject: [PATCH 45/50] Update header reader to read the whole line rather than use the file stream to make sure that any trailing white space is ignored --- utils/PhaseSpace/hdf5bods.cc | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index c61fa7a5b..c845edb45 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -250,13 +250,21 @@ std::vector read_ascii_lines(const std::string& ascii_file, int& num_aux_ints, int& num_aux_floats) { + // Open ASCII file std::ifstream infile(ascii_file); if (!infile.is_open()) { throw std::runtime_error("Could not open ASCII file: " + ascii_file); } // Read header - if (!(infile >> num_particles >> num_aux_ints >> num_aux_floats)) { + std::string header_line; + if (std::getline(infile, header_line)) { + std::istringstream header_stream(header_line); + if (!(header_stream >> num_particles >> num_aux_ints >> num_aux_floats)) { + throw std::runtime_error("Failed to parse header line: " + header_line); + } + } + else { throw std::runtime_error("Failed to read header from ASCII file"); } From 9ec0a742d29c3ada18dd157e232c2227f5cef784 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Mon, 6 Jul 2026 12:47:25 -0400 Subject: [PATCH 46/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- src/Component.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Component.cc b/src/Component.cc index 3b61983f3..c9a2b9ce7 100644 --- a/src/Component.cc +++ b/src/Component.cc @@ -1704,7 +1704,8 @@ void Component::read_bodies_and_distribute_hdf5(void) // Some ranks may legitimately receive 0 particles (e.g., more ranks than particles). // In that case, ship an empty batch and skip all dataset reads. if (nbodies_table[n] == 0) { - pf->ShipParticles(n, 0, 0); + unsigned zero = 0; + pf->ShipParticles(n, 0, zero); continue; } From 96280c24f5aa91e6822b95d8e6b2dfb1f2f30bc1 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Mon, 6 Jul 2026 12:48:29 -0400 Subject: [PATCH 47/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/PhaseSpace/hdf5bods.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index a644ffb91..3d2b0ba1a 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -69,10 +69,10 @@ style by porting over the ascii_to_hdf5 routine(). */ // C++ std +#include #include #include #include -#include #include #include #include From c92dcadb9241ddd90a28a0cf9232cb08cea85801 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Mon, 6 Jul 2026 12:49:18 -0400 Subject: [PATCH 48/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- src/Component.H | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/Component.H b/src/Component.H index 18fdc4efe..48dc2c257 100644 --- a/src/Component.H +++ b/src/Component.H @@ -208,11 +208,10 @@ struct loadb_datum
  • The index dataset is optional. If it is not present, the particles are assigned indices in the order they are read. -
  • The aux_int_* and aux_float_* - datasets are optional. If they are not present, the corresponding - vectors are empty. The number of auxiliary integer and float - datasets is determined by the num_aux_ints and - num_aux_floats attributes, respectively. +
  • The aux_int_* and aux_float_* datasets are optional. + If omitted, the corresponding attribute counts (num_aux_ints, num_aux_floats) + should be set to 0; if the counts are non-zero, datasets aux_int_0.. and + aux_float_0.. are expected to exist.
  • The datatype for the m, x, y, z, u, v, From a956778e118ee569ea137596bb0977cd68988368 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Mon, 6 Jul 2026 12:51:16 -0400 Subject: [PATCH 49/50] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/PhaseSpace/hdf5bods.cc | 2 -- 1 file changed, 2 deletions(-) diff --git a/utils/PhaseSpace/hdf5bods.cc b/utils/PhaseSpace/hdf5bods.cc index 3d2b0ba1a..b09d0e120 100644 --- a/utils/PhaseSpace/hdf5bods.cc +++ b/utils/PhaseSpace/hdf5bods.cc @@ -208,8 +208,6 @@ void parse_particle_line(const std::string& line, ParticleDataTemplate& data) { std::istringstream iss(line); - T val; - // Parse optional index field if (has_index) { unsigned long idx = 0; From 01079aa6b2eb239dcd1245a4e218e74ffaca7667 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 6 Jul 2026 13:16:05 -0400 Subject: [PATCH 50/50] Some previous reworking during review left these variable undefined --- src/Component.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Component.cc b/src/Component.cc index c9a2b9ce7..526f65a01 100644 --- a/src/Component.cc +++ b/src/Component.cc @@ -1740,8 +1740,8 @@ void Component::read_bodies_and_distribute_hdf5(void) pf->ShipParticles(n, 0, nbodies_table[n]); - icount = 0; - ibufcount = 0; + unsigned icount = 0; + unsigned ibufcount = 0; while (icount < nbodies_table[n]) { PartPtr part = std::make_shared(niattrib, ndattrib);