From d3027cf628f5590c7a9c918fb4b02b1ab60cc95b Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 29 Nov 2025 16:50:53 -0500 Subject: [PATCH 01/38] Preliminary implementation of subsample features in PotAccel and Cube to test a parallel output class to OutCoef called OutSample --- include/EmpCylSL.H | 2 +- pyEXP/CoefWrappers.cc | 3 +- src/AxisymmetricBasis.cc | 1 + src/Basis.H | 2 +- src/Basis.cc | 2 - src/CMakeLists.txt | 4 +- src/Cube.H | 19 +++ src/Cube.cc | 17 +- src/OutSample.H | 44 ++++++ src/OutSample.cc | 94 +++++++++++ src/PotAccel.H | 100 ++++++++++++ src/PotAccel.cc | 329 +++++++++++++++++++++++++++++++++++++++ src/SphericalBasis.cc | 7 +- 13 files changed, 615 insertions(+), 9 deletions(-) create mode 100644 src/OutSample.H create mode 100644 src/OutSample.cc diff --git a/include/EmpCylSL.H b/include/EmpCylSL.H index ebb2af947..386e0b975 100644 --- a/include/EmpCylSL.H +++ b/include/EmpCylSL.H @@ -205,7 +205,7 @@ protected: std::vector numbT; std::vector massT; unsigned sampT, defSampT; - + bool make_covar = false; std::vector vc, vs; diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc index 8aa420b32..59e2dbbdd 100644 --- a/pyEXP/CoefWrappers.cc +++ b/pyEXP/CoefWrappers.cc @@ -37,7 +37,8 @@ void CoefficientClasses(py::module &m) { coefficient structure using Python. To do this, use the constructor to make a blank instance, assign the dimensions and use assign() to create a data matrix with the supplied matrix or - array. The dimen- sions are: + array. The dimensions are: + 1. (lmax, nmax) for SphStruct 2. (mmax, nmax) for a CylStruct 3. (nmaxx, nmaxy, nmaxz) for a SlabStruct diff --git a/src/AxisymmetricBasis.cc b/src/AxisymmetricBasis.cc index 6409501b0..7617eca0f 100644 --- a/src/AxisymmetricBasis.cc +++ b/src/AxisymmetricBasis.cc @@ -21,6 +21,7 @@ AxisymmetricBasis::valid_keys = "pcaeof", "pcadiag", "pcavtk", + "covar", "subsamp", "hexp", "snr", diff --git a/src/Basis.H b/src/Basis.H index 119467da6..51229aa38 100644 --- a/src/Basis.H +++ b/src/Basis.H @@ -3,6 +3,7 @@ #include "PotAccel.H" #include +#include //! Defines a basis-based potential and acceleration class class Basis : public PotAccel @@ -39,7 +40,6 @@ public: double *tdens0, double *dpotl0, double *tdens, double *tpotl, double *tpotr, double *tpotz, double *tpotp) = 0; - /** @name Utility functions */ // @{ diff --git a/src/Basis.cc b/src/Basis.cc index 8eb381a1f..5d94de5a9 100644 --- a/src/Basis.cc +++ b/src/Basis.cc @@ -110,5 +110,3 @@ void Basis::sinecosine_R(int mmax, double phi, } } } - - diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 42682157f..fa7ea2057 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -8,8 +8,8 @@ set(exp_SOURCES Basis.cc Bessel.cc Component.cc Cube.cc Cylinder.cc generateRelaxation.cc HaloBulge.cc incpos.cc incvel.cc ComponentContainer.cc OutAscii.cc OutHDF5.cc OutMulti.cc OutRelaxation.cc OrbTrace.cc OutDiag.cc OutLog.cc OutVel.cc - OutCoef.cc multistep.cc parse.cc SlabSL.cc step.cc tidalField.cc - ultra.cc ultrasphere.cc MPL.cc OutFrac.cc OutCalbr.cc + OutCoef.cc OutSample.cc multistep.cc parse.cc SlabSL.cc step.cc + tidalField.cc ultra.cc ultrasphere.cc MPL.cc OutFrac.cc OutCalbr.cc ParticleFerry.cc chkSlurm.c chkTimer.cc GravKernel.cc CenterFile.cc PolarBasis.cc FlatDisk.cc signals.cc) diff --git a/src/Cube.H b/src/Cube.H index 38c27a4ec..b42e8eb59 100644 --- a/src/Cube.H +++ b/src/Cube.H @@ -250,6 +250,25 @@ private: } + /** Return a subsample of the basis coefficients and covariance + elements for matrices for PCA analysis. The default + implementation returns empty vectors. This parallels the + implementation in expui. */ + virtual CovarElement getSubsample() + { + std::vector mass; + std::vector count; + std::vector>> covar; + + return std::make_tuple(mass, count, covar); + } + + //! Write basis-specific parameters to HDF5 covariance file + void writeCovarH5Params(HighFive::File& file); + + //! Write basis-specific parameters to HDF5 covariance file + unsigned writeCovarH5(HighFive::Group& group, + unsigned count, double time); public: //! Id string diff --git a/src/Cube.cc b/src/Cube.cc index deee46575..496b1e31a 100644 --- a/src/Cube.cc +++ b/src/Cube.cc @@ -15,7 +15,8 @@ Cube::valid_keys = { "nmaxy", "nmaxz", "method", - "wrap" + "wrap", + "subsampleFloat" }; //@{ @@ -119,6 +120,11 @@ void Cube::initialize(void) if (conf["nmaxz" ]) nmaxz = conf["nmaxz" ].as(); if (conf["method"]) cuMethod = conf["method"].as(); if (conf["wrap" ]) wrap = conf["wrap" ].as(); + + if (conf["subsampleFloat"]) { + floatType = conf["subsampleFloat"].as(); + } + } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing parameters in Cube: " @@ -837,3 +843,12 @@ void Cube::compute_multistep_coefficients() } } +void Cube::writeCovarH5Params(HighFive::File& file) +{ + file.createAttribute("nminx", HighFive::DataSpace::From(nminx)).write(nminx); + file.createAttribute("nminy", HighFive::DataSpace::From(nminy)).write(nminy); + file.createAttribute("nminz", HighFive::DataSpace::From(nminz)).write(nminz); + file.createAttribute("nmaxx", HighFive::DataSpace::From(nmaxx)).write(nmaxx); + file.createAttribute("nmaxy", HighFive::DataSpace::From(nmaxy)).write(nmaxy); + file.createAttribute("nmaxz", HighFive::DataSpace::From(nmaxz)).write(nmaxz); +} diff --git a/src/OutSample.H b/src/OutSample.H new file mode 100644 index 000000000..fac9595ef --- /dev/null +++ b/src/OutSample.H @@ -0,0 +1,44 @@ +#ifndef _OutSample_H +#define _OutSample_H + +#include "OutSample.H" + +/** Dump coefficients at each interval + + @param filename of the coefficient file + + @param name of the component to dump + + @param nint is the frequency between file updates +*/ +class OutSample : public Output +{ + +private: + + std::string filename; + double prev = -std::numeric_limits::max(); + Component *tcomp; + + void initialize(void); + + //! Valid keys for YAML configurations + static const std::set valid_keys; + +public: + + //! Constructor + OutSample(const YAML::Node& conf); + + //! Generate the output + /*! + \param nstep is the current time step used to decide whether or not + to dump + \param last should be true on final step to force phase space dump + indepentently of whether or not the frequency criterion is met + */ + void Run(int nstep, int mstep, bool last); + +}; + +#endif diff --git a/src/OutSample.cc b/src/OutSample.cc new file mode 100644 index 000000000..e1c5b6d99 --- /dev/null +++ b/src/OutSample.cc @@ -0,0 +1,94 @@ +#include +#include +#include + +#include "expand.H" +#include "Basis.H" +#include "OutSample.H" + +const std::set +OutSample::valid_keys = { + "filename", + "nint", + "name" +}; + +OutSample::OutSample(const YAML::Node& conf) : Output(conf) +{ + nint = 10; + nintsub = std::numeric_limits::max(); + tcomp = NULL; + + initialize(); + + if (!(tcomp->force->hasSubsample())) { + throw std::runtime_error("OutSample: no subsampling for this force"); + } + +} + +void OutSample::initialize() +{ + // Remove matched keys + // + for (auto v : valid_keys) current_keys.erase(v); + + // Assign values from YAML + // + try { + if (conf["nint"]) nint = conf["nint"].as(); + if (conf["name"]) + { // Search for desired component + std::string tmp = conf["name"].as(); + for (auto c : comp->components) { + if (!(c->name.compare(tmp))) tcomp = c; + } + } + + if (!tcomp) { + std::string message = "OutSample: no component to trace. Please specify " + "the component name using the 'name' parameter."; + throw std::runtime_error(message); + } + + if (conf["filename"]) + { + filename = outdir + conf["filename"].as(); + } + else + { + filename = outdir + "outcoef." + tcomp->name + "." + runtag; + } + + } + catch (YAML::Exception & error) { + if (myid==0) std::cout << "Error parsing parameters in OutSample: " + << error.what() << std::endl + << std::string(60, '-') << std::endl + << "Config node" << std::endl + << std::string(60, '-') << std::endl + << conf << std::endl + << std::string(60, '-') << std::endl; + throw std::runtime_error("OutSample::initialize: error parsing YAML"); + } +} + +void OutSample::Run(int n, int mstep, bool last) +{ + // Subsampling is not available for output + // + if (not tcomp->force->subsampleReady()) return; + + // Skip this master step + // + if (n % nint != 0 && !last) return; + + // Check for repeat time + // + if (tnow <= prev) return; + + prev = tnow; + + // Write the subsample data + tcomp->force->writeSubsample(); +} diff --git a/src/PotAccel.H b/src/PotAccel.H index 39f10fd3e..67dbfa9d9 100644 --- a/src/PotAccel.H +++ b/src/PotAccel.H @@ -9,10 +9,18 @@ #include #include +#include + +#include +#include +#include +#include +#include #include "Particle.H" #include "StringTok.H" #include "YamlCheck.H" +#include "localmpi.H" #include "config_exp.h" @@ -93,6 +101,74 @@ protected: //! Current YAML keys to check configuration std::set current_keys; + //! @name Subsampling + //@{ + /** Define the type returned by getSubsample + + The outermost vectors contain subsamples of mass per sample, + counts per sample, and arrays of coefficients and covariance + elements per sample indexed by harmonic + */ + using CovarElement = + std::tuple< std::vector, std::vector, + std::vector>> >; + + /** Return a subsample of the basis coefficients and covariance + elements for matrices for PCA analysis. The default + implementation returns empty vectors. */ + virtual CovarElement getSubsample() + { + std::vector mass; + std::vector count; + std::vector>> covar; + + return std::make_tuple(mass, count, covar); + } + + //! Reset subsample data + virtual void subsampleReset() { subsampleComputed = false; } + + //! Extend the subsample data + virtual void extendSubsample(const std::string& fname, double time); + + //! Write basis-specific parameters to HDF5 covariance file + virtual void writeCovarH5Params(HighFive::File& file) + { + if (myid==0) + std::cout << "Basis::writeCovarH5Params: " + << "not implemented for basis type <" << id << ">" + << std::endl; + } + + //! Write basis-specific parameters to HDF5 covariance file + virtual unsigned writeCovarH5(HighFive::Group& group, + unsigned count, double time); + + //! Variance file versioning (Version 1.0 has more granularity and + //! poor compression) + inline static const std::string CovarianceFileVersion = "1.1"; + + //! Use 32-bit floating point type in covariance file + bool floatType = false; + + //! Save covariance subsamples to file + bool fullCovar = false; + + //! Is subsample computed and ready to write? + bool subsampleComputed = false; + + //@} + + //@{ + //! HDF5 compression settings + unsigned H5compress = 5; + unsigned H5chunk = 1024*1024; // (1 MB) + bool H5shuffle = true; + bool H5szip = false; + //@} + + public: //! For timing data @@ -181,6 +257,27 @@ public: virtual void dump_coefs(ostream &out) {}; virtual void dump_coefs_h5(const std::string &file) {}; + //@name Subsampling methods + //@{ + //! Write the subsample to file + virtual void writeSubsample(); + + //! Indicate whether or not subsamples are available + virtual bool subsampleReady() { return subsampleComputed; } + + //! Indicate whether or not subsampling is supported + virtual bool hasSubsample() { return false; } + //@} + + //! HDF5 compression settings + void setCovarH5Compress(unsigned level, unsigned chunksize, bool shuffle, bool szip=false) + { + H5compress = level; + H5chunk = chunksize; + H5shuffle = shuffle; + H5szip = szip; + } + /** Update the multi time step force algorithm when moving particle i from level cur to level next @@ -235,6 +332,9 @@ public: //! Get unaccounted keys std::set unmatched() { return current_keys; } + //! Round time key to emulated fixed-point arithmetic + double roundTime(double time); + }; //! Helper class to pass info to threaded member diff --git a/src/PotAccel.cc b/src/PotAccel.cc index 8449b2ef3..ec6cc606e 100644 --- a/src/PotAccel.cc +++ b/src/PotAccel.cc @@ -2,6 +2,10 @@ #include #include +#include + +#include + #include "expand.H" #include "PotAccel.H" @@ -363,3 +367,328 @@ bool ltdub(const pair& A, return A.first < B.first; } + +void PotAccel::writeSubsample() +{ + // Check that subsample data is available + // + if (not hasSubsample()) { + std::cout << "PotAccel::writeSubsample: " + << "subsample data is not available. " << std::endl; + return; + } + + // Only root process writes + // + if (myid) return; + + // Round time + // + double time = roundTime(tnow); + + // The H5 filename + // + std::string fname = "subsample." + cC->name + "." + runtag + ".h5"; + + // Check if file exists? + // + try { + // Open the HDF5 file in read-write mode, creating if it doesn't + // exist + HighFive::File file(fname, + HighFive::File::ReadWrite | + HighFive::File::Create); + + // Check for version string + std::string path = "CovarianceFileVersion"; + + // Check for valid HDF file by attribute + if (file.hasAttribute(path)) { + extendSubsample(fname, time); + return; + } + + // Write the Version string + // + file.createAttribute("CovarianceFileVersion", HighFive::DataSpace::From(CovarianceFileVersion)).write(CovarianceFileVersion); + + // Write the basis identifier string + // + file.createAttribute("BasisID", HighFive::DataSpace::From(id)).write(id); + + // Write the data type size + // + int sz = 8; if (floatType) sz = 4; + file.createAttribute("FloatSize", HighFive::DataSpace::From(sz)).write(sz); + + // Write the specific parameters + // + writeCovarH5Params(file); + + // Group count variable + // + unsigned count = 0; + HighFive::DataSet dataset = file.createDataSet("count", count); + + // Create a new group for coefficient snapshots + // + HighFive::Group group = file.createGroup("snapshots"); + + // Write the coefficients + // + count = writeCovarH5(group, count, time); + + // Update the count + // + dataset.write(count); + + } catch (const HighFive::Exception& err) { + // Handle HighFive specific errors (e.g., file not found) + throw std::runtime_error + (std::string("PotAccel::writeCoefCovariance HighFive Error: ") + err.what()); + } catch (const std::exception& err) { + // Handle other general exceptions + throw std::runtime_error + (std::string("PotAccel::writeCoefCovariance Error: ") + err.what()); + } + + // Reset subsample data + subsampleReset(); +} + +void PotAccel::extendSubsample(const std::string& fname, double time) +{ + try { + // Open an hdf5 file + // + HighFive::File file(fname, HighFive::File::ReadWrite); + + // Get the dataset + HighFive::DataSet dataset = file.getDataSet("count"); + + unsigned count; + dataset.read(count); + + HighFive::Group group = file.getGroup("snapshots"); + + // Write the coefficients + // + count = writeCovarH5(group, count, time); + + // Update the count + // + dataset.write(count); + + } catch (HighFive::Exception& err) { + throw std::runtime_error + (std::string("PotAccel::extendSubsample: HighFive error: ") + err.what()); + } +} + +double PotAccel::roundTime(double time) +{ + // Eight decimal places should be enough here... + const double multiplier = 1.0e+08; // std::pow(10.0, 8); + return std::floor(time * multiplier + 0.5) / multiplier; +} + +unsigned PotAccel::writeCovarH5(HighFive::Group& snaps, unsigned count, double time) +{ + std::ostringstream stim; + stim << std::setw(8) << std::setfill('0') << std::right << count++; + + // Get the subsample data + // + auto [sampleMasses, sampleCounts, covdata] = getSubsample(); + + // Make a new group for this time + // + HighFive::Group stanza = snaps.createGroup(stim.str()); + + // Add a time attribute + // + time = roundTime(time); + stanza.createAttribute("Time", HighFive::DataSpace::From(time)).write(time); + + // Enable compression + // + auto dcpl1 = HighFive::DataSetCreateProps{}; // sample stats + auto dcpl2 = HighFive::DataSetCreateProps{}; // coefficients + auto dcpl3 = HighFive::DataSetCreateProps{}; // covariance + + // Properties for sample stats + if (H5compress) { + unsigned int csz = sampleCounts.size(); + dcpl1.add(HighFive::Chunking({csz, 1})); + if (H5shuffle) dcpl1.add(HighFive::Shuffle()); + dcpl1.add(HighFive::Deflate(H5compress)); + } + + // Add the sample statistics + // + HighFive::DataSet s1data = stanza.createDataSet("sampleCounts", sampleCounts, dcpl1); + HighFive::DataSet s2data = stanza.createDataSet("sampleMasses", sampleMasses, dcpl1); + + + // Number of samples + // + unsigned sampleSize = covdata.size(); + unsigned ltot = covdata[0].size(); + unsigned nmax = std::get<0>(covdata[0][0]).rows(); + unsigned diagonalSize = nmax; + + // Save variance or full covariance + if (fullCovar) diagonalSize = nmax*(nmax + 1)/2; + + // Add data dimensions + // + stanza.createAttribute + ("sampleSize", HighFive::DataSpace::From(sampleSize)).write(sampleSize); + + stanza.createAttribute + ("angularSize", HighFive::DataSpace::From(ltot)).write(ltot); + + stanza.createAttribute + ("rankSize", HighFive::DataSpace::From(nmax)).write(nmax); + + int icov = fullCovar ? 1 : 0; + stanza.createAttribute + ("fullCovar", HighFive::DataSpace::From(icov)).write(icov); + + if (H5compress) { + // Szip parameters + const int options_mask = H5_SZIP_NN_OPTION_MASK; + const int pixels_per_block = 8; + + // Properties for coefficients + // + unsigned int csz2 = nmax * ltot * sampleSize; + HighFive::Chunking data_dims2{std::min(csz2, H5chunk), 1}; + + dcpl2.add(data_dims2); + if (H5shuffle) dcpl2.add(HighFive::Shuffle()); + if (H5szip) { + dcpl2.add(HighFive::Szip(options_mask, pixels_per_block)); + } else { + dcpl2.add(HighFive::Deflate(H5compress)); + } + + // Properties for covariance + // + unsigned int csz3 = ltot * diagonalSize * sampleSize; + HighFive::Chunking data_dims3{std::min(csz3, H5chunk), 1}; + + dcpl3.add(data_dims3); + if (H5shuffle) dcpl3.add(HighFive::Shuffle()); + if (H5szip) { + dcpl3.add(HighFive::Szip(options_mask, pixels_per_block)); + } else { + dcpl3.add(HighFive::Deflate(H5compress)); + } + } + + // Pack the coefficient data + // + if (floatType) { + // Create a vector of doubles for the real and imaginary parts + Eigen::VectorXf real_part(nmax*ltot*sampleSize); + Eigen::VectorXf imag_part(nmax*ltot*sampleSize); + + for (size_t T=0, c=0; T(covdata[T][l])(n)); + imag_part(c) = std::imag(std::get<0>(covdata[T][l])(n)); + } + } + } + + // Create two separate, compressed datasets + stanza.createDataSet("coefficients_real", real_part, dcpl2); + stanza.createDataSet("coefficients_imag", imag_part, dcpl2); + + if (std::get<1>(covdata[0][0]).size()) { + + real_part.resize(ltot*diagonalSize*sampleSize); + imag_part.resize(ltot*diagonalSize*sampleSize); + + for (size_t T=0, c=0; T(covdata[T][l])(n1, n2)); + imag_part(c) = std::imag(std::get<1>(covdata[T][l])(n1, n2)); + } + } + // Pack the diagonal only + // + else { + real_part(c ) = std::real(std::get<1>(covdata[T][l])(n1, n1)); + imag_part(c++) = std::imag(std::get<1>(covdata[T][l])(n1, n1)); + } + } + } + } + // Create two separate, compressed datasets + stanza.createDataSet("covariance_real", real_part, dcpl3); + stanza.createDataSet("covariance_imag", imag_part, dcpl3); + } + + } else { + Eigen::VectorXd real_part(ltot*nmax*sampleSize); + Eigen::VectorXd imag_part(ltot*nmax*sampleSize); + + for (size_t T=0, c=0; T(covdata[T][l])(n)); + imag_part(c) = std::imag(std::get<0>(covdata[T][l])(n)); + } + } + } + + // Create two separate, compressed datasets + // + stanza.createDataSet("coefficients_real", real_part, dcpl2); + stanza.createDataSet("coefficients_imag", imag_part, dcpl2); + + if (std::get<1>(covdata[0][0]).size()) { + + real_part.resize(ltot*diagonalSize*sampleSize); + imag_part.resize(ltot*diagonalSize*sampleSize); + + for (size_t T=0, c=0; T(covdata[T][l])(n1, n2)); + imag_part(c) = std::imag(std::get<1>(covdata[T][l])(n1, n2)); + } + } + // Pack the diagonal only + // + else { + real_part(c ) = std::real(std::get<1>(covdata[T][l])(n1, n1)); + imag_part(c++) = std::imag(std::get<1>(covdata[T][l])(n1, n1)); + } + } + } + } + + // Create two separate, compressed datasets + // + stanza.createDataSet("covariance_real", real_part, dcpl3); + stanza.createDataSet("covariance_imag", imag_part, dcpl3); + } + } + // END: sample loop + + return count; +} + diff --git a/src/SphericalBasis.cc b/src/SphericalBasis.cc index a520ff6d0..ffa632971 100644 --- a/src/SphericalBasis.cc +++ b/src/SphericalBasis.cc @@ -45,7 +45,8 @@ SphericalBasis::valid_keys = { "playback", "coefCompute", "coefMaster", - "orthocheck" + "orthocheck", + "subsampleFloat" }; SphericalBasis::SphericalBasis(Component* c0, const YAML::Node& conf, MixtureBasis *m) : @@ -132,6 +133,10 @@ SphericalBasis::SphericalBasis(Component* c0, const YAML::Node& conf, MixtureBas if (conf["noise_model_file"]) noise_model_file = conf["noise_model_file"].as(); + if (conf["subsampleFloat"]) { + floatType = conf["subsampleFloat"].as(); + } + if (conf["seedN"]) seedN = conf["seedN"].as(); if (conf["ssfrac"]) { From 5d6ea214c1b65ebff113e8f674ddb792b3c40c4f Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 30 Nov 2025 19:54:57 -0500 Subject: [PATCH 02/38] Additional implementation errors fixed and added some documentation --- src/Cube.H | 49 ++++++---- src/Cube.cc | 207 ++++++++++++++++++++++++++++++++++++++++- src/OutSample.H | 13 ++- src/OutSample.cc | 37 +++----- src/OutputContainer.cc | 5 + src/PotAccel.H | 38 +++++--- src/PotAccel.cc | 1 + 7 files changed, 291 insertions(+), 59 deletions(-) diff --git a/src/Cube.H b/src/Cube.H index b42e8eb59..657a192e8 100644 --- a/src/Cube.H +++ b/src/Cube.H @@ -65,15 +65,35 @@ private: //@{ //! Variables int imx, imy, imz, osize; - int use1, use0; + int use1, use0, sampT=0, nint=0; std::complex kfac; - double dfac; + double dfac, last=-std::numeric_limits::max(); bool wrap = true; // Enforce periodic wrapping //@} //! Valid keys for YAML configurations static const std::set valid_keys; + //@{ + //! Covariance structures. First index is T, second is the + //! flattened 3-d k vector + std::vector meanV; + std::vector covrV; + + //@{ + //! Per thread covariance structures + std::vector> meanV1; + std::vector> covrV1; + std::vector workV1; + std::vector countV1; + std::vector massV1; + //@} + + void init_covariance(); + void zero_covariance(); + //@} + + #if HAVE_LIBCUDA==1 virtual void determine_coefficients_cuda(); virtual void determine_acceleration_cuda(); @@ -250,25 +270,9 @@ private: } - /** Return a subsample of the basis coefficients and covariance - elements for matrices for PCA analysis. The default - implementation returns empty vectors. This parallels the - implementation in expui. */ - virtual CovarElement getSubsample() - { - std::vector mass; - std::vector count; - std::vector>> covar; - - return std::make_tuple(mass, count, covar); - } - //! Write basis-specific parameters to HDF5 covariance file void writeCovarH5Params(HighFive::File& file); - //! Write basis-specific parameters to HDF5 covariance file - unsigned writeCovarH5(HighFive::Group& group, - unsigned count, double time); public: //! Id string @@ -288,6 +292,15 @@ public: //! Coefficient output void dump_coefs_h5(const std::string& file); + + /** Return a subsample of the basis coefficients and covariance + elements for matrices for analysis. The default implementation + returns empty vectors. This parallels the implementation in + expui. */ + CovarElement getSubsample(); + + //! Indicate whether or not subsampling is supported + virtual bool hasSubsample() { return true; } }; diff --git a/src/Cube.cc b/src/Cube.cc index 496b1e31a..e16650372 100644 --- a/src/Cube.cc +++ b/src/Cube.cc @@ -16,6 +16,8 @@ Cube::valid_keys = { "nmaxz", "method", "wrap", + "nint", + "samplesz", "subsampleFloat" }; @@ -36,6 +38,7 @@ Cube::Cube(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf) coef_dump = true; byPlanes = true; cuMethod = "planes"; + sampT = 100; // Default parameter values // @@ -94,6 +97,15 @@ Cube::Cube(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf) // dfac = 2.0*M_PI; kfac = std::complex(0.0, dfac); + + // Initialize covariance + // + init_covariance(); + +#if HAVE_LIBCUDA==1 + cuda_initalize(); +#endif + } Cube::~Cube(void) @@ -121,10 +133,19 @@ void Cube::initialize(void) if (conf["method"]) cuMethod = conf["method"].as(); if (conf["wrap" ]) wrap = conf["wrap" ].as(); + if (conf["nint"]) { + nint = conf["nint" ].as(); + if (nint>0) computeSubsample = true; + } + if (conf["subsampleFloat"]) { floatType = conf["subsampleFloat"].as(); } + if (conf["samplesz"]) { + sampT = conf["samplesz"].as(); + } + } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing parameters in Cube: " @@ -140,12 +161,81 @@ void Cube::initialize(void) if (myid==0) std::cout << "---- Cube::initialize: wrap=" << std::boolalpha << wrap << std::endl; +} -#if HAVE_LIBCUDA==1 - cuda_initialize(); -#endif +void Cube::init_covariance() +{ + if (computeSubsample) { + + meanV.resize(sampT); + for (auto& v : meanV) { + v.resize(osize); + } + + workV1.resize(nthrds); + for (auto& v : workV1) v.resize(osize); + + if (fullCovar) { + covrV.resize(sampT); + for (auto& v : covrV) { + v.resize(osize, osize); + } + } else { + covrV.clear(); + } + + sampleCounts.resize(sampT); + sampleMasses.resize(sampT); + + meanV1.resize(nthrds); + covrV1.resize(nthrds); + countV1.resize(nthrds); + massV1.resize(nthrds); + + for (int n=0; n facx, facy, facz; int ix, iy, iz; + if (requestSubsample) workV1[id].setZero(); + for (facx=startx, ix=0; ixPart(i); if (part->indx < 10) { @@ -250,9 +347,31 @@ void * Cube::determine_coefficients_thread(void * arg) << std::endl; } } + // END: deepDebug } + // END: iz loop } + // END: iy loop } + // END: ix loop + + if (requestSubsample) { + // Which subsample bin? + // + int T = q % sampT; + + // Accumulate counts and masses + // + countV1[id](T) += 1; + massV1[id](T) += mass; + + // Accumulate subsample contributions + // + meanV1[id][T] += workV1[id] * mass; + if (fullCovar) + covrV1[id][T] += workV1[id] * workV1[id].adjoint() * mass; + } + } // END: particle loop } @@ -293,6 +412,19 @@ void Cube::determine_coefficients(void) for (int n=0; n::max()) { + if (nint>0 and this_step % nint == 0) { + if (tnow > last) { + requestSubsample = true; + last = tnow; + init_covariance(); + } + } + } else { + subsampleComputed = false; + } + #if HAVE_LIBCUDA==1 (*barrier)("Cube::entering cuda coefficients", __FILE__, __LINE__); if (component->cudaDevice>=0 and use_cuda) { @@ -334,9 +466,51 @@ void Cube::determine_coefficients(void) // Last level? // if (multistep and mlevel==multistep) { - compute_multistep_coefficients(); + compute_multistep_coefficients(); } + // Accumulate mean and covariance subsample contributions + // + if (requestSubsample) { + + // Only finalize at the last multistep level + // + if ( (multistep and mlevel==multistep) or multistep==0 ) { + + // Sum over threads + // + for (int n=1; n("nmaxy", HighFive::DataSpace::From(nmaxy)).write(nmaxy); file.createAttribute("nmaxz", HighFive::DataSpace::From(nmaxz)).write(nmaxz); } + + +PotAccel::CovarElement Cube::getSubsample() +{ + using covTuple = std::tuple; + // Prepare the covariance structure + std::vector< std::vector > covar(sampT); + + // Resize the inner vector to 1 (one component for Cube) + for (auto & v : covar) v.resize(1); + + // Fill the covariance structure with subsamples + for (int T=0; T OutSample::valid_keys = { "filename", - "nint", "name" }; OutSample::OutSample(const YAML::Node& conf) : Output(conf) { - nint = 10; nintsub = std::numeric_limits::max(); tcomp = NULL; @@ -36,14 +34,13 @@ void OutSample::initialize() // Assign values from YAML // try { - if (conf["nint"]) nint = conf["nint"].as(); - if (conf["name"]) - { // Search for desired component - std::string tmp = conf["name"].as(); - for (auto c : comp->components) { - if (!(c->name.compare(tmp))) tcomp = c; - } + // Search for desired component by name + if (conf["name"]) { + std::string tmp = conf["name"].as(); + for (auto c : comp->components) { + if (!(c->name.compare(tmp))) tcomp = c; } + } if (!tcomp) { std::string message = "OutSample: no component to trace. Please specify " @@ -51,14 +48,11 @@ void OutSample::initialize() throw std::runtime_error(message); } - if (conf["filename"]) - { - filename = outdir + conf["filename"].as(); - } - else - { - filename = outdir + "outcoef." + tcomp->name + "." + runtag; - } + if (conf["filename"]) { + filename = outdir + conf["filename"].as(); + } else { + filename = outdir + "outcoef." + tcomp->name + "." + runtag; + } } catch (YAML::Exception & error) { @@ -78,17 +72,14 @@ void OutSample::Run(int n, int mstep, bool last) // Subsampling is not available for output // if (not tcomp->force->subsampleReady()) return; - - // Skip this master step - // - if (n % nint != 0 && !last) return; - + // Check for repeat time // if (tnow <= prev) return; + // Cache the current simulation time prev = tnow; - + // Write the subsample data tcomp->force->writeSubsample(); } diff --git a/src/OutputContainer.cc b/src/OutputContainer.cc index c26235e94..64673ea3c 100644 --- a/src/OutputContainer.cc +++ b/src/OutputContainer.cc @@ -20,6 +20,7 @@ #include "OutFrac.H" #include "OutCalbr.H" #include "OutMulti.H" +#include "OutSample.H" OutputContainer::OutputContainer() { @@ -112,6 +113,10 @@ void OutputContainer::initialize(void) out.push_back(new OutCalbr(node)); } + else if ( !name.compare("outsamp") ) { + out.push_back(new OutSample(node)); + } + else { string msg("I don't know about the output type: "); msg += name; diff --git a/src/PotAccel.H b/src/PotAccel.H index 67dbfa9d9..efdc9c0ea 100644 --- a/src/PotAccel.H +++ b/src/PotAccel.H @@ -42,6 +42,16 @@ public: //! For timing typedef std::vector TList; + /** Define the type returned by getSubsample + + The outermost vectors contain subsamples of mass per sample, + counts per sample, and arrays of coefficients and covariance + elements per sample indexed by harmonic + */ + using CovarElement = + std::tuple< Eigen::VectorXd, Eigen::VectorXi, + std::vector>> >; + private: // Threading stuff @@ -103,23 +113,14 @@ protected: //! @name Subsampling //@{ - /** Define the type returned by getSubsample - - The outermost vectors contain subsamples of mass per sample, - counts per sample, and arrays of coefficients and covariance - elements per sample indexed by harmonic - */ - using CovarElement = - std::tuple< std::vector, std::vector, - std::vector>> >; /** Return a subsample of the basis coefficients and covariance elements for matrices for PCA analysis. The default implementation returns empty vectors. */ virtual CovarElement getSubsample() { - std::vector mass; - std::vector count; + Eigen::VectorXd mass; + Eigen::VectorXi count; std::vector>> covar; @@ -155,6 +156,12 @@ protected: //! Save covariance subsamples to file bool fullCovar = false; + //! Request to compute subsample data + bool requestSubsample = false; + + //! Enable computation of subsample data + bool computeSubsample = false; + //! Is subsample computed and ready to write? bool subsampleComputed = false; @@ -262,6 +269,9 @@ public: //! Write the subsample to file virtual void writeSubsample(); + //! Set flag to compute subsample data + virtual void subsampleRequest() { requestSubsample = true; } + //! Indicate whether or not subsamples are available virtual bool subsampleReady() { return subsampleComputed; } @@ -269,6 +279,12 @@ public: virtual bool hasSubsample() { return false; } //@} + //@{ + //! Sample counts and masses for covariance computation + Eigen::VectorXi sampleCounts; + Eigen::VectorXd sampleMasses; + //@} + //! HDF5 compression settings void setCovarH5Compress(unsigned level, unsigned chunksize, bool shuffle, bool szip=false) { diff --git a/src/PotAccel.cc b/src/PotAccel.cc index ec6cc606e..85acb92fb 100644 --- a/src/PotAccel.cc +++ b/src/PotAccel.cc @@ -5,6 +5,7 @@ #include #include +#include #include "expand.H" #include "PotAccel.H" From 56fb55c1c326e2ae17fba6b05fddbde848738a5e Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 1 Dec 2025 09:47:32 -0500 Subject: [PATCH 03/38] Typo in call name --- src/Cube.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Cube.cc b/src/Cube.cc index e16650372..62b1685a8 100644 --- a/src/Cube.cc +++ b/src/Cube.cc @@ -103,7 +103,7 @@ Cube::Cube(Component* c0, const YAML::Node& conf) : PotAccel(c0, conf) init_covariance(); #if HAVE_LIBCUDA==1 - cuda_initalize(); + cuda_initialize(); #endif } From 322ab22dcade118cf610092875ec0821e10a5f61 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 1 Dec 2025 17:44:52 -0500 Subject: [PATCH 04/38] Refactoring into a single covariance class for simplicity --- expui/BiorthBasis.H | 153 +++++------ expui/BiorthBasis.cc | 586 +---------------------------------------- expui/CMakeLists.txt | 8 +- expui/Covariance.cc | 559 +++++++++++++++++++++++++++++++++++++++ include/Covariance.H | 177 +++++++++++++ pyEXP/BasisWrappers.cc | 167 ++---------- src/Cube.H | 5 +- src/Cube.cc | 10 +- src/OutSample.H | 10 +- src/OutSample.cc | 41 ++- src/PotAccel.H | 87 ++---- src/PotAccel.cc | 325 ----------------------- src/SphericalBasis.cc | 4 - 13 files changed, 910 insertions(+), 1222 deletions(-) create mode 100644 expui/Covariance.cc create mode 100644 include/Covariance.H diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index dba7b1671..e01f275b6 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -21,6 +21,7 @@ #include "BiorthBess.H" #include "BasisFactory.H" #include "BiorthCube.H" +#include "Covariance.H" #include "sltableMP2.H" #include "SLGridMP2.H" #include "YamlCheck.H" @@ -106,22 +107,14 @@ namespace BasisClasses //! Coefficient variance computation enabled bool pcavar = false; - //! Data type for covariance - bool floatType = false; - //@{ //! Sample counts and masses for covariance computation Eigen::VectorXi sampleCounts; Eigen::VectorXd sampleMasses; //@} - //@{ - //! HDF5 compression settings - unsigned H5compress = 5; - unsigned H5chunk = 1024*1024; // (1 MB) - bool H5shuffle = true; - bool H5szip = false; - //@} + //! Covariance storage instance + std::shared_ptr covarStore; //! Round time key to emulated fixed-point arithmetic double roundTime(double time) @@ -131,19 +124,6 @@ namespace BasisClasses return std::floor(time * multiplier + 0.5) / multiplier; } - //! Write H5 covariance; returns updated group count - virtual unsigned writeCovarH5(HighFive::Group& group, unsigned count, double time); - - //! Write H5 parameters for coefficient covariance - virtual void writeCovarH5Params(HighFive::File& file) {} - - //! Add coefficient covariance data to an HDF5 file - virtual void extendCoefCovariance(const std::string& filename, double time); - - //! Variance file versioning (Version 1.0 has more granularity and - //! poor compression) - inline static const std::string CovarianceFileVersion = "1.1"; - //! Store full coavariance matrix? bool covar = true; @@ -290,51 +270,36 @@ namespace BasisClasses return varray; } - //! Return a vector of tuples of basis functions and the - //! covariance matrix for subsamples of particles - using CoefCovarType = std::tuple; - virtual std::vector> - getCoefCovariance() - { - // Must be overriden; base implementation throws error - throw std::runtime_error("BiorthBasis::getCoefCovariance: Not implemented for this basis"); - } - - //! Get sample counts for the covariance computation - virtual std::tuple - getCovarSamples() - { - return std::make_tuple(sampleCounts, sampleMasses); - } - //! Write coefficient covariance data to an HDF5 file - void writeCoefCovariance(const std::string& compname, - const std::string& runtag, - double time=0.0); - - //! Choose between float and double storage for covariance - void setCovarFloatType(bool flt) { floatType = flt; } - - //! HDF5 compression settings - void setCovarH5Compress(unsigned level, unsigned chunksize, bool shuffle, bool szip=false) + virtual void writeCoefCovariance(const std::string& compname, const std::string& runtag, + double time=0.0) { - H5compress = level; - H5chunk = chunksize; - H5shuffle = shuffle; - H5szip = szip; + // Must be overriden; base implementation throws error + throw std::runtime_error("BiorthBasis::writeCoefCovariance: " + "Not implemented for this basis"); } //! Make covariance after accumulation virtual void makeCoefCovariance(void) {} //! Enable covariance computation with optional sample time - virtual void enableCoefCovariance(bool pcavar, int sampT_in, bool covar_in) + virtual void enableCoefCovariance(bool pcavar, int sampT_in, bool ftype, bool covar_in) { // Must be overriden; base implementation throws error throw std::runtime_error("BiorthBasis::enableCoefCovariance: " "Not implemented for this basis"); } + //! HDF5 compression settings + void setCovarH5Compress(unsigned level, unsigned chunksize, bool shuffle, bool szip=false) + { + if (covarStore) { + covarStore->setCovarH5Compress(level, chunksize, shuffle, szip); + } else { + throw std::runtime_error("BiorthBasis::setCovarH5Compress: covariance storage not initialized"); + } + } + //! Evaluate acceleration in Cartesian coordinates in centered //! coordinate system for a collections of points virtual RowMatrixXd& getAccel(Eigen::Ref pos) @@ -518,18 +483,35 @@ namespace BasisClasses return ret; } - //! Return a vector of tuples of basis functions and the - //! covariance matrix for subsamples of particles - std::vector> - getCoefCovariance(); + //! Write coefficient covariance data to an HDF5 file + virtual void writeCoefCovariance(const std::string& compname, const std::string& runtag, double time=0.0) + { + if (covarStore) { + std::vector> covarData(meanV.size()); + for (int T=0; TwriteCoefCovariance(compname, runtag, elements, time); + } else { + throw std::runtime_error("Spherical::writeCoefCovariance: covariance storage not initialized"); + } + } virtual void enableCoefCovariance(bool pcavar_in, int sampT_in=100, - bool covar_in=true) + bool ftype=false, bool covar_in=true) { pcavar = pcavar_in; sampT = sampT_in; covar = covar_in; - if (pcavar) init_covariance(); + if (pcavar) { + init_covariance(); + covarStore = std::make_shared + ([this](HighFive::File& file){this->writeCovarH5Params(file);}, BasisID, ftype, covar); + } } //! Create and the coefficients from a function callback with the //! provided time value. Set potential=true if function is a @@ -990,7 +972,6 @@ namespace BasisClasses //! For coefficient writing typedef Eigen::Matrix EigenColMajor; - //@{ //! Basis construction parameters @@ -1103,24 +1084,8 @@ namespace BasisClasses return ret; } - //! Return a vector of tuples of basis functions and the - //! covariance matrix for subsamples of particles - std::vector> - getCoefCovariance() - { - return sl->getCoefCovariance(); - } - - std::tuple - getCovarSamples() - { - // Copy to base class members - std::tie(sampleCounts, sampleMasses) = sl->getCovarSamples(); - return {sampleCounts, sampleMasses}; - } - virtual void enableCoefCovariance(bool pcavar_in, int sampT_in=100, - bool covar_in=true) + bool ftype=false, bool covar_in=true) { pcavar = pcavar_in; sampT = sampT_in; @@ -1128,6 +1093,8 @@ namespace BasisClasses if (pcavar) { sl->setSampT(std::max(1, sampT)); sl->set_covar(true); + covarStore = std::make_shared + ([this](HighFive::File& file){this->writeCovarH5Params(file);}, BasisID, ftype, covar); } } //! Create and the coefficients from a function callback with the @@ -1443,25 +1410,35 @@ namespace BasisClasses return true; } - //! Return a vector of tuples of basis functions and the - //! covariance matrix for subsamples of particles - std::vector> getCoefCovariance(); - - //! Get sample count and mass arrays - std::tuple - getCovarSamples() + //! Write coefficient covariance data to an HDF5 file + virtual void writeCoefCovariance(const std::string& compname, const std::string& runtag, double time=0.0) { - return {sampleCounts, sampleMasses}; + if (covarStore) { + std::vector> covarData(meanV.size()); + for (int T=0; TwriteCoefCovariance(compname, runtag, elements, time); + } else { + throw std::runtime_error("Spherical::writeCoefCovariance: covariance storage not initialized"); + } } + //! Enable coefficient covariance computation virtual void enableCoefCovariance(bool pcavar_in, int sampT_in=100, - bool covar_in=false) + bool ftype=false, bool covar_in=false) { pcavar = pcavar_in; sampT = sampT_in; covar = covar_in; - if (pcavar) { init_covariance(); } + if (pcavar) { + init_covariance(); + covarStore = std::make_shared + ([this](HighFive::File& file){this->writeCovarH5Params(file);}, BasisID, ftype, covar); + } } //! Return wave-number indices from flattened index diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index e555bed01..d4e6806e0 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1017,29 +1017,6 @@ namespace BasisClasses return ret; } - - /** Return a vector of tuples of basis functions and the covariance - matrix for subsamples of particles */ - std::vector> - Spherical::getCoefCovariance() - { - std::vector> ret; - if (pcavar) { - ret.resize(sampT); - for (int T=0; T(ret[T][l]) = meanV[T][l]; - std::get<1>(ret[T][l]) = covrV[T][l]; - } - } - } - - return ret; - } - - #define MINEPS 1.0e-10 void BiorthBasis::legendre_R(int lmax, double x, Eigen::MatrixXd& p) @@ -3831,7 +3808,8 @@ namespace BasisClasses "check", "method", "pcavar," - "subsamp" + "subsamp", + "nint" }; Cube::Cube(const YAML::Node& CONF) : BiorthBasis(CONF, "cube") @@ -4131,6 +4109,7 @@ namespace BasisClasses /** Return a vector of tuples of basis functions and the covariance matrix for subsamples of particles for Cube type */ + /* std::vector> Cube::getCoefCovariance() { @@ -4151,6 +4130,7 @@ namespace BasisClasses return ret; } + */ std::vector Cube::crt_eval(double x, double y, double z) @@ -5061,564 +5041,6 @@ namespace BasisClasses file.createAttribute("nmaxy", HighFive::DataSpace::From(nmaxy)).write(nmaxy); file.createAttribute("nmaxz", HighFive::DataSpace::From(nmaxz)).write(nmaxz); } - - unsigned BiorthBasis::writeCovarH5(HighFive::Group& snaps, unsigned count, double time) - { - std::ostringstream stim; - stim << std::setw(8) << std::setfill('0') << std::right << count++; - - // Make a new group for this time - // - HighFive::Group stanza = snaps.createGroup(stim.str()); - - // Add a time attribute - // - time = roundTime(time); - stanza.createAttribute("Time", HighFive::DataSpace::From(time)).write(time); - - // Enable compression - // - auto dcpl1 = HighFive::DataSetCreateProps{}; // sample stats - auto dcpl2 = HighFive::DataSetCreateProps{}; // coefficients - auto dcpl3 = HighFive::DataSetCreateProps{}; // covariance - - // Properties for sample stats - if (H5compress) { - unsigned int csz = sampleCounts.size(); - dcpl1.add(HighFive::Chunking({csz, 1})); - if (H5shuffle) dcpl1.add(HighFive::Shuffle()); - dcpl1.add(HighFive::Deflate(H5compress)); - } - - // Add the sample statistics - // - HighFive::DataSet s1data = stanza.createDataSet("sampleCounts", sampleCounts, dcpl1); - HighFive::DataSet s2data = stanza.createDataSet("sampleMasses", sampleMasses, dcpl1); - - // Covariance data - // - auto covdata = getCoefCovariance(); - - // Number of samples - // - unsigned sampleSize = covdata.size(); - unsigned ltot = covdata[0].size(); - unsigned nmax = std::get<0>(covdata[0][0]).rows(); - unsigned diagonalSize = nmax; - - // Save variance or full covariance - if (covar) diagonalSize = nmax*(nmax + 1)/2; - - // Add data dimensions - // - stanza.createAttribute - ("sampleSize", HighFive::DataSpace::From(sampleSize)).write(sampleSize); - - stanza.createAttribute - ("angularSize", HighFive::DataSpace::From(ltot)).write(ltot); - - stanza.createAttribute - ("rankSize", HighFive::DataSpace::From(nmax)).write(nmax); - - int icov = covar ? 1 : 0; - stanza.createAttribute - ("fullCovar", HighFive::DataSpace::From(icov)).write(icov); - - if (H5compress) { - // Szip parameters - const int options_mask = H5_SZIP_NN_OPTION_MASK; - const int pixels_per_block = 8; - - // Properties for coefficients - // - unsigned int csz2 = nmax * ltot * sampleSize; - HighFive::Chunking data_dims2{std::min(csz2, H5chunk), 1}; - - dcpl2.add(data_dims2); - if (H5shuffle) dcpl2.add(HighFive::Shuffle()); - if (H5szip) { - dcpl2.add(HighFive::Szip(options_mask, pixels_per_block)); - } else { - dcpl2.add(HighFive::Deflate(H5compress)); - } - - // Properties for covariance - // - unsigned int csz3 = ltot * diagonalSize * sampleSize; - HighFive::Chunking data_dims3{std::min(csz3, H5chunk), 1}; - - dcpl3.add(data_dims3); - if (H5shuffle) dcpl3.add(HighFive::Shuffle()); - if (H5szip) { - dcpl3.add(HighFive::Szip(options_mask, pixels_per_block)); - } else { - dcpl3.add(HighFive::Deflate(H5compress)); - } - } - - // Pack the coefficient data - // - if (floatType) { - // Create a vector of doubles for the real and imaginary parts - Eigen::VectorXf real_part(nmax*ltot*sampleSize); - Eigen::VectorXf imag_part(nmax*ltot*sampleSize); - - for (size_t T=0, c=0; T(covdata[T][l])(n)); - imag_part(c) = std::imag(std::get<0>(covdata[T][l])(n)); - } - } - } - - // Create two separate, compressed datasets - stanza.createDataSet("coefficients_real", real_part, dcpl2); - stanza.createDataSet("coefficients_imag", imag_part, dcpl2); - - if (std::get<1>(covdata[0][0]).size()) { - - real_part.resize(ltot*diagonalSize*sampleSize); - imag_part.resize(ltot*diagonalSize*sampleSize); - - for (size_t T=0, c=0; T(covdata[T][l])(n1, n2)); - imag_part(c) = std::imag(std::get<1>(covdata[T][l])(n1, n2)); - } - } - // Pack the diagonal only - // - else { - real_part(c ) = std::real(std::get<1>(covdata[T][l])(n1, n1)); - imag_part(c++) = std::imag(std::get<1>(covdata[T][l])(n1, n1)); - } - } - } - } - // Create two separate, compressed datasets - stanza.createDataSet("covariance_real", real_part, dcpl3); - stanza.createDataSet("covariance_imag", imag_part, dcpl3); - } - - } else { - Eigen::VectorXd real_part(ltot*nmax*sampleSize); - Eigen::VectorXd imag_part(ltot*nmax*sampleSize); - - for (size_t T=0, c=0; T(covdata[T][l])(n)); - imag_part(c) = std::imag(std::get<0>(covdata[T][l])(n)); - } - } - } - - // Create two separate, compressed datasets - // - stanza.createDataSet("coefficients_real", real_part, dcpl2); - stanza.createDataSet("coefficients_imag", imag_part, dcpl2); - - if (std::get<1>(covdata[0][0]).size()) { - - real_part.resize(ltot*diagonalSize*sampleSize); - imag_part.resize(ltot*diagonalSize*sampleSize); - - for (size_t T=0, c=0; T(covdata[T][l])(n1, n2)); - imag_part(c) = std::imag(std::get<1>(covdata[T][l])(n1, n2)); - } - } - // Pack the diagonal only - // - else { - real_part(c ) = std::real(std::get<1>(covdata[T][l])(n1, n1)); - imag_part(c++) = std::imag(std::get<1>(covdata[T][l])(n1, n1)); - } - } - } - } - - // Create two separate, compressed datasets - // - stanza.createDataSet("covariance_real", real_part, dcpl3); - stanza.createDataSet("covariance_imag", imag_part, dcpl3); - } - } - // END: sample loop - - return count; - } - - void BiorthBasis::writeCoefCovariance(const std::string& compname, const std::string& runtag, double time) - { - // Check that variance computation is on - // - if (not pcavar) { - std::cout << "BiorthBasis::writeCoefCovariance: covariance computation is disabled. " - << "Set 'pcavar: true' to enable." << std::endl; - return; - } - - // Only root process writes - // - if (myid) return; - - // Check that there is something to write - // - int totalCount = 0; - std::tie(sampleCounts, sampleMasses) = getCovarSamples(); - totalCount += sampleCounts.sum(); - - if (totalCount==0) { - std::cout << "BiorthBasis::writeCoefCovariance: no data" << std::endl; - return; - } - - // Round time - // - time = roundTime(time); - - // The H5 filename - // - std::string fname = "coefcovar." + compname + "." + runtag + ".h5"; - - // Check if file exists? - // - try { - // Open the HDF5 file in read-write mode, creating if it doesn't - // exist - HighFive::File file(fname, - HighFive::File::ReadWrite | - HighFive::File::Create); - - // Check for version string - std::string path = "CovarianceFileVersion"; - - // Check for valid HDF file by attribute - if (file.hasAttribute(path)) { - extendCoefCovariance(fname, time); - return; - } - - // Write the Version string - // - file.createAttribute("CovarianceFileVersion", HighFive::DataSpace::From(CovarianceFileVersion)).write(CovarianceFileVersion); - - // Write the basis identifier string - // - file.createAttribute("BasisID", HighFive::DataSpace::From(BasisID)).write(BasisID); - - // Write the data type size - // - int sz = 8; if (floatType) sz = 4; - file.createAttribute("FloatSize", HighFive::DataSpace::From(sz)).write(sz); - - // Write the specific parameters - // - writeCovarH5Params(file); - - // Group count variable - // - unsigned count = 0; - HighFive::DataSet dataset = file.createDataSet("count", count); - - // Create a new group for coefficient snapshots - // - HighFive::Group group = file.createGroup("snapshots"); - - // Write the coefficients - // - count = writeCovarH5(group, count, time); - - // Update the count - // - dataset.write(count); - - } catch (const HighFive::Exception& err) { - // Handle HighFive specific errors (e.g., file not found) - throw std::runtime_error - (std::string("BiorthBasis::writeCoefCovariance HighFive Error: ") + err.what()); - } catch (const std::exception& err) { - // Handle other general exceptions - throw std::runtime_error - (std::string("BiorthBasis::writeCoefCovariance Error: ") + err.what()); - } - } - - void BiorthBasis::extendCoefCovariance(const std::string& fname, double time) - { - try { - // Open an hdf5 file - // - HighFive::File file(fname, HighFive::File::ReadWrite); - - // Get the dataset - HighFive::DataSet dataset = file.getDataSet("count"); - - unsigned count; - dataset.read(count); - - HighFive::Group group = file.getGroup("snapshots"); - - // Write the coefficients - // - count = writeCovarH5(group, count, time); - - // Update the count - // - dataset.write(count); - - } catch (HighFive::Exception& err) { - throw std::runtime_error - (std::string("BiorthBasis::extendCoefCovariance: HighFive error: ") + err.what()); - } - } - - // Read covariance data - CovarianceReader::CovarianceReader(const std::string& filename, int stride) - { - try { - // Open an existing hdf5 file for reading - // - HighFive::File file(filename, HighFive::File::ReadOnly); - - // Write the Version string - // - std::string version; - file.getAttribute("CovarianceFileVersion").read(version); - // Check for alpha version - if (version == std::string("1.0")) { - throw std::runtime_error("CovarianceReader: this is an early alpha test version. Please remake your files"); - } - // Test for current version - if (version != std::string("1.1")) { - throw std::runtime_error(std::string("CovarianceReader: unsupported file version, ") + version); - } - - // Read the basis identifier string - // - file.getAttribute("BasisID").read(basisID); - - // Get the float size - int sz = 8; - file.getAttribute("FloatSize").read(sz); - if (sz != 4 and sz != 8) { - std::ostringstream sout; - sout << "CovarianceReader: unsupported float size, " << sz; - throw std::runtime_error(sout.str()); - } - - int lmax, nmax, ltot; - - // Current implemented spherical types - const std::set sphereType = {"Spherical", "SphereSL", "Bessel"}; - - // Currently implemented cylindrical types - const std::set cylinderType = {"Cylindrical"}; - - if (sphereType.find(basisID) != sphereType.end()) { - file.getAttribute("lmax").read(lmax); - file.getAttribute("nmax").read(nmax); - ltot = (lmax+1)*(lmax+2)/2; - } else if (cylinderType.find(basisID) != cylinderType.end()) { - file.getAttribute("mmax").read(lmax); - file.getAttribute("nmax").read(nmax); - ltot = lmax + 1; - } else if (basisID == "Cube") { - int nmaxx, nmaxy, nmaxz; - file.getAttribute("nmaxx").read(nmaxx); - file.getAttribute("nmaxy").read(nmaxy); - file.getAttribute("nmaxz").read(nmaxz); - ltot = (2*nmaxx + 1) * (2*nmaxy + 1) * (2*nmaxz + 1); - } else { - throw std::runtime_error(std::string("CovarianceReader: unknown or unimplemented covariance for basis type, ") + basisID); - } - - // Group count variable - // - unsigned count = 0; - file.getDataSet("count").read(count); - - // Open the snapshot group - // - auto snaps = file.getGroup("snapshots"); - - for (unsigned n=0; n(Time * fixedPointPrecision + 0.5); - timeMap[itime] = times.size(); - times.push_back(Time); - - // Get sample properties - // - sampleCounts.push_back(Eigen::VectorXi()); - stanza.getDataSet("sampleCounts").read(sampleCounts.back()); - - sampleMasses.push_back(Eigen::VectorXd()); - stanza.getDataSet("sampleMasses").read(sampleMasses.back()); - - // Get data attributes - // - int nT, lSize, rank, icov=1; - stanza.getAttribute("sampleSize") .read(nT); - stanza.getAttribute("angularSize").read(lSize); - stanza.getAttribute("rankSize") .read(rank); - - if (stanza.hasAttribute("fullCovar")) - stanza.getAttribute("fullCovar").read(icov); - - // Full covariance or variance? - // - bool covar = icov==1 ? true : false; - - // Allocate sample vector for current time - // - covarData.push_back(std::vector>(nT)); - - // Storage - // - Eigen::VectorXcd data0, data1; - - // Get the flattened coefficient array - // - if (sz==4) { - // Get the real and imaginary parts - // - Eigen::VectorXf data_real = - stanza.getDataSet("coefficients_real").read(); - - Eigen::VectorXf data_imag = - stanza.getDataSet("coefficients_imag").read(); - - // Resize the complex array and assign - // - data0.resize(data_real.size()); - data0.real() = data_real.cast(); - data0.imag() = data_imag.cast(); - - // Check for existence of covariance - // - if (stanza.exist("covariance_real")) { - - data_real = - stanza.getDataSet("covariance_real").read(); - - data_imag = - stanza.getDataSet("covariance_imag").read(); - - // Resize the complex array and assign - data1.resize(data_real.size()); - data1.real() = data_real.cast(); - data1.imag() = data_imag.cast(); - } - } else { - // Get the real and imaginary parts - Eigen::VectorXd data_real = - stanza.getDataSet("coefficients_real").read(); - - Eigen::VectorXd data_imag = - stanza.getDataSet("coefficients_imag").read(); - - // Resize the complex array and assign - data0.resize(data_real.size()); - data0.real() = data_real; - data0.imag() = data_imag; - - // Check for existence of covariance - // - if (stanza.exist("covariance_real")) { - - // Get the real and imaginary parts - data_real = - stanza.getDataSet("covariance_real").read(); - - data_imag = - stanza.getDataSet("covariance_imag").read(); - - // Resize the complex array and assign - data1.resize(data_real.size()); - data1.real() = data_real; - data1.imag() = data_imag; - } - } - - // Positions in data stanzas - int sCof = 0, sCov = 0; - - // Loop through all indices and repack - for (int T=0; T elem(lSize); - for (auto & e : elem) { - // Coefficients - std::get<0>(e).resize(rank); - // Covariance matrix - if (data1.size()) std::get<1>(e).resize(rank, rank); - } - - // Pack the coefficient data - int c = 0; - for (size_t l=0; l(elem[l])(n) = data0(sCof + c++); - } - } - sCof += c; - - // Pack the covariance data - c = 0; - for (size_t l=0; l(elem[l])(n1, n2) = data1(sCov + c++); - if (n1 != n2) - std::get<1>(elem[l])(n2, n1) = std::get<1>(elem[l])(n1, n2); - } - else { - if (n1==n2) - std::get<1>(elem[l])(n1, n2) = data1(sCov + c++); - else - std::get<1>(elem[l])(n1, n2) = 0.0; - } - } - } - } - } - sCov += c; - - // Add the data - covarData.back()[T] = std::move(elem); - } - // END: sample loop - } - // END: snapshot loop - } catch (HighFive::Exception& err) { - std::cerr << err.what() << std::endl; - } - } CoefClasses::CoefStrPtr Spherical::makeFromFunction (std::function func, diff --git a/expui/CMakeLists.txt b/expui/CMakeLists.txt index fa993f75e..2cad9de0d 100644 --- a/expui/CMakeLists.txt +++ b/expui/CMakeLists.txt @@ -35,10 +35,12 @@ endif() # Make the expui shared library # + set(expui_SOURCES BasisFactory.cc BiorthBasis.cc FieldBasis.cc - CoefContainer.cc CoefStruct.cc FieldGenerator.cc expMSSA.cc - Coefficients.cc KMeans.cc Centering.cc ParticleIterator.cc - Koopman.cc BiorthBess.cc SvdSignChoice.cc KDdensity.cc) + Covariance.cc CoefContainer.cc CoefStruct.cc FieldGenerator.cc + expMSSA.cc Coefficients.cc KMeans.cc Centering.cc + ParticleIterator.cc Koopman.cc BiorthBess.cc SvdSignChoice.cc + KDdensity.cc) add_library(expui ${expui_SOURCES}) set_target_properties(expui PROPERTIES OUTPUT_NAME expui) target_include_directories(expui PUBLIC ${common_INCLUDE}) diff --git a/expui/Covariance.cc b/expui/Covariance.cc new file mode 100644 index 000000000..b1f5af840 --- /dev/null +++ b/expui/Covariance.cc @@ -0,0 +1,559 @@ +#include + +#include "Covariance.H" + +namespace BasisClasses +{ + + unsigned SubsampleCovariance::writeCovarH5 + (HighFive::Group& snaps, CovarData& elem, unsigned count, double time) + { + std::ostringstream stim; + stim << std::setw(8) << std::setfill('0') << std::right << count++; + + // Make a new group for this time + // + HighFive::Group stanza = snaps.createGroup(stim.str()); + + // Add a time attribute + // + time = roundTime(time); + stanza.createAttribute("Time", HighFive::DataSpace::From(time)).write(time); + + // Enable compression + // + auto dcpl1 = HighFive::DataSetCreateProps{}; // sample stats + auto dcpl2 = HighFive::DataSetCreateProps{}; // coefficients + auto dcpl3 = HighFive::DataSetCreateProps{}; // covariance + + // Number of samples + // + unsigned sampleSize = std::get<2>(elem).size(); + unsigned ltot = std::get<2>(elem)[0].size(); + unsigned nmax = std::get<0>(std::get<2>(elem)[0][0]).size(); + unsigned diagonalSize = nmax; + + // Properties for sample stats + if (H5compress) { + unsigned int csz = sampleSize; + dcpl1.add(HighFive::Chunking({csz, 1})); + if (H5shuffle) dcpl1.add(HighFive::Shuffle()); + dcpl1.add(HighFive::Deflate(H5compress)); + } + + // Add the sample statistics + // + HighFive::DataSet s1data = stanza.createDataSet("sampleCounts", std::get<0>(elem), dcpl1); + HighFive::DataSet s2data = stanza.createDataSet("sampleMasses", std::get<1>(elem), dcpl1); + + // Save variance or full covariance + if (covar) diagonalSize = nmax*(nmax + 1)/2; + + // Add data dimensions + // + stanza.createAttribute + ("sampleSize", HighFive::DataSpace::From(sampleSize)).write(sampleSize); + + stanza.createAttribute + ("angularSize", HighFive::DataSpace::From(ltot)).write(ltot); + + stanza.createAttribute + ("rankSize", HighFive::DataSpace::From(nmax)).write(nmax); + + int icov = covar ? 1 : 0; + stanza.createAttribute + ("fullCovar", HighFive::DataSpace::From(icov)).write(icov); + + if (H5compress) { + // Szip parameters + const int options_mask = H5_SZIP_NN_OPTION_MASK; + const int pixels_per_block = 8; + + // Properties for coefficients + // + unsigned int csz2 = nmax * ltot * sampleSize; + HighFive::Chunking data_dims2{std::min(csz2, H5chunk), 1}; + + dcpl2.add(data_dims2); + if (H5shuffle) dcpl2.add(HighFive::Shuffle()); + if (H5szip) { + dcpl2.add(HighFive::Szip(options_mask, pixels_per_block)); + } else { + dcpl2.add(HighFive::Deflate(H5compress)); + } + + // Properties for covariance + // + unsigned int csz3 = ltot * diagonalSize * sampleSize; + HighFive::Chunking data_dims3{std::min(csz3, H5chunk), 1}; + + dcpl3.add(data_dims3); + if (H5shuffle) dcpl3.add(HighFive::Shuffle()); + if (H5szip) { + dcpl3.add(HighFive::Szip(options_mask, pixels_per_block)); + } else { + dcpl3.add(HighFive::Deflate(H5compress)); + } + } + + // Pack the coefficient data + // + if (floatType) { + // Create a vector of doubles for the real and imaginary parts + Eigen::VectorXf real_part(nmax*ltot*sampleSize); + Eigen::VectorXf imag_part(nmax*ltot*sampleSize); + + for (size_t T=0, c=0; T(std::get<2>(elem)[T][l])(n)); + imag_part(c) = std::imag(std::get<0>(std::get<2>(elem)[T][l])(n)); + } + } + } + + // Create two separate, compressed datasets + stanza.createDataSet("coefficients_real", real_part, dcpl2); + stanza.createDataSet("coefficients_imag", imag_part, dcpl2); + + if (std::get<1>(std::get<2>(elem)[0][0]).size()) { + + real_part.resize(ltot*diagonalSize*sampleSize); + imag_part.resize(ltot*diagonalSize*sampleSize); + + for (size_t T=0, c=0; T(std::get<2>(elem)[T][l])(n1, n2)); + imag_part(c) = std::imag(std::get<1>(std::get<2>(elem)[T][l])(n1, n2)); + } + } + // Pack the diagonal only + // + else { + real_part(c ) = std::real(std::get<1>(std::get<2>(elem)[T][l])(n1, n1)); + imag_part(c++) = std::imag(std::get<1>(std::get<2>(elem)[T][l])(n1, n1)); + } + } + } + } + // Create two separate, compressed datasets + stanza.createDataSet("covariance_real", real_part, dcpl3); + stanza.createDataSet("covariance_imag", imag_part, dcpl3); + } + + } else { + Eigen::VectorXd real_part(ltot*nmax*sampleSize); + Eigen::VectorXd imag_part(ltot*nmax*sampleSize); + + for (size_t T=0, c=0; T(std::get<2>(elem)[T][l])(n)); + imag_part(c) = std::imag(std::get<0>(std::get<2>(elem)[T][l])(n)); + } + } + } + + // Create two separate, compressed datasets + // + stanza.createDataSet("coefficients_real", real_part, dcpl2); + stanza.createDataSet("coefficients_imag", imag_part, dcpl2); + + if (std::get<1>(std::get<2>(elem)[0][0]).size()) { + + real_part.resize(ltot*diagonalSize*sampleSize); + imag_part.resize(ltot*diagonalSize*sampleSize); + + for (size_t T=0, c=0; T(std::get<2>(elem)[T][l])(n1, n2)); + imag_part(c) = std::imag(std::get<1>(std::get<2>(elem)[T][l])(n1, n2)); + } + } + // Pack the diagonal only + // + else { + real_part(c ) = std::real(std::get<1>(std::get<2>(elem)[T][l])(n1, n1)); + imag_part(c++) = std::imag(std::get<1>(std::get<2>(elem)[T][l])(n1, n1)); + } + } + } + } + + // Create two separate, compressed datasets + // + stanza.createDataSet("covariance_real", real_part, dcpl3); + stanza.createDataSet("covariance_imag", imag_part, dcpl3); + } + } + // END: sample loop + + return count; + } + + void SubsampleCovariance::writeCoefCovariance + (const std::string& compname, const std::string& runtag, CovarData& elem, double time) + { + // Only root process writes + // + if (myid) return; + + // Check that there is something to write + // + int totalCount = 0; + totalCount += std::get<0>(elem).sum(); + + if (totalCount==0) { + std::cout << "SubsampleCovariance::writeCoefCovariance: no data" << std::endl; + return; + } + + // Round time + // + time = roundTime(time); + + // The H5 filename + // + std::string fname = "coefcovar." + compname + "." + runtag + ".h5"; + + // Check if file exists? + // + try { + // Open the HDF5 file in read-write mode, creating if it doesn't + // exist + HighFive::File file(fname, + HighFive::File::ReadWrite | + HighFive::File::Create); + + // Check for version string + std::string path = "CovarianceFileVersion"; + + // Check for valid HDF file by attribute + if (file.hasAttribute(path)) { + extendCoefCovariance(fname, elem, time); + return; + } + + // Write the Version string + // + file.createAttribute("CovarianceFileVersion", HighFive::DataSpace::From(CovarianceFileVersion)).write(CovarianceFileVersion); + + // Write the basis identifier string + // + file.createAttribute("BasisID", HighFive::DataSpace::From(BasisID)).write(BasisID); + + // Write the data type size + // + int sz = 8; if (floatType) sz = 4; + file.createAttribute("FloatSize", HighFive::DataSpace::From(sz)).write(sz); + + // Write the specific parameters + // + writeCovarH5Params(file); + + // Group count variable + // + unsigned count = 0; + HighFive::DataSet dataset = file.createDataSet("count", count); + + // Create a new group for coefficient snapshots + // + HighFive::Group group = file.createGroup("snapshots"); + + // Write the coefficients + // + count = writeCovarH5(group, elem, count, time); + + // Update the count + // + dataset.write(count); + + } catch (const HighFive::Exception& err) { + // Handle HighFive specific errors (e.g., file not found) + throw std::runtime_error + (std::string("SubsampleCovariance::writeCoefCovariance HighFive Error: ") + err.what()); + } catch (const std::exception& err) { + // Handle other general exceptions + throw std::runtime_error + (std::string("SubsampleCovariance::writeCoefCovariance Error: ") + err.what()); + } + } + + void SubsampleCovariance::extendCoefCovariance + (const std::string& fname, CovarData& elem, double time) + { + try { + // Open an hdf5 file + // + HighFive::File file(fname, HighFive::File::ReadWrite); + + // Get the dataset + HighFive::DataSet dataset = file.getDataSet("count"); + + unsigned count; + dataset.read(count); + + HighFive::Group group = file.getGroup("snapshots"); + + // Write the coefficients + // + count = writeCovarH5(group, elem, count, time); + + // Update the count + // + dataset.write(count); + + } catch (HighFive::Exception& err) { + throw std::runtime_error + (std::string("SubsampleCovariance::extendCoefCovariance: HighFive error: ") + err.what()); + } + } + + // Read covariance data + SubsampleCovariance::SubsampleCovariance + (const std::string& filename, int stride) + { + try { + // Open an existing hdf5 file for reading + // + HighFive::File file(filename, HighFive::File::ReadOnly); + + // Write the Version string + // + std::string version; + file.getAttribute("CovarianceFileVersion").read(version); + // Check for alpha version + if (version == std::string("1.0")) { + throw std::runtime_error("CovarianceReader: this is an early alpha test version. Please remake your files"); + } + // Test for current version + if (version != std::string("1.1")) { + throw std::runtime_error(std::string("CovarianceReader: unsupported file version, ") + version); + } + + // Read the basis identifier string + // + file.getAttribute("BasisID").read(BasisID); + + // Get the float size + int sz = 8; + file.getAttribute("FloatSize").read(sz); + if (sz != 4 and sz != 8) { + std::ostringstream sout; + sout << "CovarianceReader: unsupported float size, " << sz; + throw std::runtime_error(sout.str()); + } + + int lmax, nmax, ltot; + + // Current implemented spherical types + const std::set sphereType = {"Spherical", "SphereSL", "Bessel"}; + + // Currently implemented cylindrical types + const std::set cylinderType = {"Cylindrical"}; + + std::cout << "Covariance: reading basis type " << BasisID << std::endl; + + if (sphereType.find(BasisID) != sphereType.end()) { + file.getAttribute("lmax").read(lmax); + file.getAttribute("nmax").read(nmax); + ltot = (lmax+1)*(lmax+2)/2; + } else if (cylinderType.find(BasisID) != cylinderType.end()) { + file.getAttribute("mmax").read(lmax); + file.getAttribute("nmax").read(nmax); + ltot = lmax + 1; + } else if (BasisID == "Cube") { + int nmaxx, nmaxy, nmaxz; + file.getAttribute("nmaxx").read(nmaxx); + file.getAttribute("nmaxy").read(nmaxy); + file.getAttribute("nmaxz").read(nmaxz); + ltot = (2*nmaxx + 1) * (2*nmaxy + 1) * (2*nmaxz + 1); + } else { + throw std::runtime_error(std::string("Covariance: unknown or unimplemented covariance for basis type, ") + BasisID); + } + + // Group count variable + // + unsigned count = 0; + file.getDataSet("count").read(count); + + // Open the snapshot group + // + auto snaps = file.getGroup("snapshots"); + + for (unsigned n=0; n>(nT) }; + + // Storage + // + Eigen::VectorXcd data0, data1; + + // Get the flattened coefficient array + // + if (sz==4) { + // Get the real and imaginary parts + // + Eigen::VectorXf data_real = + stanza.getDataSet("coefficients_real").read(); + + Eigen::VectorXf data_imag = + stanza.getDataSet("coefficients_imag").read(); + + // Resize the complex array and assign + // + data0.resize(data_real.size()); + data0.real() = data_real.cast(); + data0.imag() = data_imag.cast(); + + // Check for existence of covariance + // + if (stanza.exist("covariance_real")) { + + data_real = + stanza.getDataSet("covariance_real").read(); + + data_imag = + stanza.getDataSet("covariance_imag").read(); + + // Resize the complex array and assign + data1.resize(data_real.size()); + data1.real() = data_real.cast(); + data1.imag() = data_imag.cast(); + } + } else { + // Get the real and imaginary parts + Eigen::VectorXd data_real = + stanza.getDataSet("coefficients_real").read(); + + Eigen::VectorXd data_imag = + stanza.getDataSet("coefficients_imag").read(); + + // Resize the complex array and assign + data0.resize(data_real.size()); + data0.real() = data_real; + data0.imag() = data_imag; + + // Check for existence of covariance + // + if (stanza.exist("covariance_real")) { + + // Get the real and imaginary parts + data_real = + stanza.getDataSet("covariance_real").read(); + + data_imag = + stanza.getDataSet("covariance_imag").read(); + + // Resize the complex array and assign + data1.resize(data_real.size()); + data1.real() = data_real; + data1.imag() = data_imag; + } + } + + // Positions in data stanzas + int sCof = 0, sCov = 0; + + // Loop through all indices and repack + for (int T=0; T elem(lSize); + for (auto & e : elem) { + // Coefficients + std::get<0>(e).resize(rank); + // Covariance matrix + if (data1.size()) std::get<1>(e).resize(rank, rank); + } + + // Pack the coefficient data + int c = 0; + for (size_t l=0; l(elem[l])(n) = data0(sCof + c++); + } + } + sCof += c; + + // Pack the covariance data + c = 0; + for (size_t l=0; l(elem[l])(n1, n2) = data1(sCov + c++); + if (n1 != n2) + std::get<1>(elem[l])(n2, n1) = std::get<1>(elem[l])(n1, n2); + } + else { + if (n1==n2) + std::get<1>(elem[l])(n1, n2) = data1(sCov + c++); + else + std::get<1>(elem[l])(n1, n2) = 0.0; + } + } + } + } + } + sCov += c; + + // Add the data + std::get<2>(covarData[Time])[T] = std::move(elem); + } + // END: sample loop + } + // END: snapshot loop + } catch (HighFive::Exception& err) { + std::cerr << err.what() << std::endl; + } + } + +} +// END: namespace BasisClasses diff --git a/include/Covariance.H b/include/Covariance.H new file mode 100644 index 000000000..1fdb012dd --- /dev/null +++ b/include/Covariance.H @@ -0,0 +1,177 @@ +#ifndef _Covariance_H +#define _Covariance_H + +#include +#include +#include +#include + +#include +#include // For 3d rectangular grids +#include + +#include +#include +#include +#include +#include + +/* +#include "DiskDensityFunc.H" +#include "ParticleReader.H" +#include "Coefficients.H" +#include "BiorthBess.H" +#include "BasisFactory.H" +#include "BiorthCube.H" +#include "sltableMP2.H" +#include "SLGridMP2.H" +#include "YamlCheck.H" +#include "BiorthCyl.H" +#include "EmpCylSL.H" +*/ + +#include "localmpi.H" +#include "exputils.H" + +namespace BasisClasses +{ + /** + Class to maintain subsample covariance of basis functions + */ + class SubsampleCovariance + { + public: + + using ParamCallback = std::function; + + using CoefCovarType = std::tuple; + + using CovarData = std::tuple>>; + + protected: + + ParamCallback paramCallback; + + //! Variance file versioning (Version 1.0 has more granularity and + //! poor compression) + inline static const std::string CovarianceFileVersion = "1.1"; + + //! Store full covariance matrix? + bool covar = true; + + //! Use 32-bit float storage for covariance? + bool floatType = false; + + + //! Fixed-point multiplier for time mapping. Eight decimal places + //! should be enough here... + const double multiplier = 1.0e+08; + + //! Round time key to emulated fixed-point arithmetic + double roundTime(double time) + { + return std::floor(time * multiplier + 0.5) / multiplier; + } + + //! Round time key to emulated fixed-point arithmetic + int roundTimeInt(double time) + { + return std::floor(time * multiplier + 0.5); + } + + //@{ + //! HDF5 compression settings + unsigned H5compress = 5; + unsigned H5chunk = 1024*1024; // (1 MB) + bool H5shuffle = true; + bool H5szip = false; + //@} + + //! Basis identifier string + std::string BasisID; + + //! Time list + std::vector times; + + //! Imported data members + std::map covarData; + + //! Write H5 covariance; returns updated group count + virtual unsigned writeCovarH5 + (HighFive::Group& group, CovarData& elem, unsigned count, double time); + + //! Write H5 parameters for coefficient covariance + virtual void writeCovarH5Params(HighFive::File& file) + { + if (paramCallback) { + paramCallback(file); + } else { + throw std::runtime_error("SubsampleCovariance::writeCovarH5Params: no paramCallback defined"); + } + } + + //! Add coefficient covariance data to an HDF5 file + virtual void extendCoefCovariance + (const std::string& filename, CovarData& elem, double time); + + public: + + //! Constructor from scratch + SubsampleCovariance(ParamCallback func, const std::string& BasisID, + bool flt, bool cov=false) + : paramCallback(func), BasisID(BasisID), floatType(flt), covar(cov) {} + + //! Constructor from HDF5 file + SubsampleCovariance(const std::string& filename, int stride=1); + + //! Get time list + std::vector Times() { return times; } + + //! Return a vector of tuples of basis functions and the + //! covariance matrix for subsamples of particles + virtual std::vector> + getCoefCovariance(double time) + { + auto it = covarData.find(roundTime(time)); + if (it == covarData.end()) + throw std::runtime_error("SubsampleCovariance::getCoefCovariance: time not found"); + return std::get<2>(it->second); + } + + //! Get sample counts for the covariance computation + virtual std::tuple + getCovarSamples(double time) + { + auto it = covarData.find(roundTime(time)); + if (it == covarData.end()) + throw std::runtime_error("SubsampleCovariance::getCovarSamples: time not found"); + return std::make_tuple(std::get<0>(it->second), std::get<1>(it->second)); + } + + //! Write coefficient covariance data to an HDF5 file + void writeCoefCovariance(const std::string& compname, + const std::string& runtag, + CovarData& covarElem, + double time=0.0); + + //! Choose between float and double storage for covariance + void setCovarFloatType(bool flt) { floatType = flt; } + + //! HDF5 compression settings + void setCovarH5Compress(unsigned level, unsigned chunksize, bool shuffle, bool szip=false) + { + H5compress = level; + H5chunk = chunksize; + H5shuffle = shuffle; + H5szip = szip; + } + + //! Return the basis ID string that generated this covariance sample + const std::string& basisIDname() const { return BasisID; } + }; +} +// END: namespace BasisClasses + +#endif // _Covariance_H + diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 79a4e2f0d..a8a741feb 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -439,23 +439,6 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE_PURE(void, BiorthBasis, set_coefs, coefs); } - using CCelement = std::tuple; - using CCreturn = std::vector>; - - CCreturn getCoefCovariance(void) override { - PYBIND11_OVERRIDE(CCreturn, BiorthBasis, getCoefCovariance,); - } - - using Selement = std::tuple; - - Selement getCovarSamples() override { - PYBIND11_OVERRIDE(Selement, BiorthBasis, getCovarSamples,); - } - - void enableCoefCovariance(bool pcavar, int nsamples, bool covar) override { - PYBIND11_OVERRIDE(void, BiorthBasis, enableCoefCovariance, pcavar, nsamples, covar); - } - }; class PySpherical : public Spherical @@ -521,8 +504,8 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE_PURE(std::vector, Spherical, orthoCheck, knots); } - void enableCoefCovariance(bool pcavar, int nsamples, bool covar) override { - PYBIND11_OVERRIDE(void, Spherical, enableCoefCovariance, pcavar, nsamples, covar); + void enableCoefCovariance(bool pcavar, int nsamples, bool ftype, bool covar) override { + PYBIND11_OVERRIDE(void, Spherical, enableCoefCovariance, pcavar, nsamples, ftype, covar); } }; @@ -576,8 +559,8 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE(void, Cylindrical, make_coefs,); } - void enableCoefCovariance(bool pcavar, int nsamples, bool covar) override { - PYBIND11_OVERRIDE(void, Cylindrical, enableCoefCovariance, pcavar, nsamples, covar); + void enableCoefCovariance(bool pcavar, int nsamples, bool ftype, bool covar) override { + PYBIND11_OVERRIDE(void, Cylindrical, enableCoefCovariance, pcavar, nsamples, ftype, covar); } }; @@ -871,21 +854,8 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE(void, Cube, make_coefs,); } - using CCelement = std::tuple; - using CCreturn = std::vector>; - - CCreturn getCoefCovariance(void) override { - PYBIND11_OVERRIDE(CCreturn, Cube, getCoefCovariance,); - } - - using Selement = std::tuple; - - Selement getCovarSamples() override { - PYBIND11_OVERRIDE(Selement, Cube, getCovarSamples,); - } - - void enableCoefCovariance(bool pcavar, int nsamples, bool covar) override { - PYBIND11_OVERRIDE(void, Cube, enableCoefCovariance, pcavar, nsamples, covar); + void enableCoefCovariance(bool pcavar, int nsamples, bool ftype, bool covar) override { + PYBIND11_OVERRIDE(void, Cube, enableCoefCovariance, pcavar, nsamples, ftype, covar); } }; @@ -1362,62 +1332,6 @@ void BasisFactoryClasses(py::module &m) makeFromArray : create coefficients contributions )", py::arg("mass"), py::arg("pos")) - .def("getCoefCovariance", &BasisClasses::BiorthBasis::getCoefCovariance, - R"( - Return the coefficient vectors and covariance matrices for each partition of the - accumulated particles. The number of partitions is set by the configuration - parameter 'sampT' (default: 100). - - The first dimension are the time samples. The second dimension is the angular - index. Each element is a tuple of the coefficient vector and covariance. - The values are complex128 and represents the full amplitude and phase information. - - Returns - ------- - arrays: list[list[tuple[np.ndarray, np.ndarray]]] - Each element is a tuple (coef_vector, coef_covariance_matrix), - where coef_vector and coef_covariance_matrix are numpy.ndarray. - Each coef_covariance_matrix is of shape (nmax, nmax). All values are - complex128. - - Shape and Indexing - ------------------ - - The first list index is the number of time samples. - - The second list index is the angular elements. For spherical bases, all (l, m) pairs - are in triangular index order with l in [0,...,lmax] and m in [0,...,l] for a total of - (lmax+1)*(lmax+2)/2 entries. For cylindrical bases, there are (mmax+1) harmonic - entries for each value m in [0,...,mmax]. - - Each covariance matrix is of shape (nmax, nmax), where nmax is the number of basis - functions - - See also - -------- - getCovarSamples : get counts and mass in each partition - writeCoefCovariance : write the coefficient vectors, covariance - matrices, and metadata to an HDF5 file. - )") - .def("getCovarSamples", &BasisClasses::BiorthBasis::getCovarSamples, - R"( - Return a vector of counts and mass used to compute the partitioned - vectors and covariance. The length of both returned vectors equals 'sampT' - - Parameters - ---------- - None - - Returns - ------- - arrays: tuple(ndarray, ndarray) - a tuple of arrays containing the counts and mass in each - partitioned sample - - See also - -------- - getCoefCovariance : get the coefficient vectors and covariance matrices - for the partitioned phase space. - writeCoefCovariance : write the coefficient vectors, covariance - matrices, and metadata to an HDF5 file. - )") .def("writeCoefCovariance", &BasisClasses::BiorthBasis::writeCoefCovariance, R"( Write the partitioned coefficient vectors and covariance matrices @@ -1449,21 +1363,6 @@ void BasisFactoryClasses(py::module &m) for the partitioned phase space. getCovarSamples : get counts and mass in each partition )", py::arg("compname"), py::arg("runtag"), py::arg("time")=0.0) - .def("setCovarFloatType", &BasisClasses::BiorthBasis::setCovarFloatType, - R"( - Set the floating point type used for covariance storage in HDF5. Set to - true for float4. Type float8 is the default (false). This does not affect - the return type in getCoefCovariance(). - - Parameters - ---------- - floatType : bool - Use 'float32' rather than 'float64' if true - - Returns - ------- - None - )") .def("enableCoefCovariance", &BasisClasses::BiorthBasis::enableCoefCovariance, R"( Enable or disable the coefficient covariance computation and set the @@ -1475,6 +1374,9 @@ void BasisFactoryClasses(py::module &m) enable (true) or disable (false) the covariance computation nsamples : int number of time partitions to use for covariance computation + ftype: bool + if true, use float32 for covariance storage; if false, + use float64 (default: false) covar: bool if true, compute and save covariance to the HDF5 file; if false, save mean and variance vectors only (default: true) @@ -1482,7 +1384,8 @@ void BasisFactoryClasses(py::module &m) Returns ------- None - )", py::arg("pcavar"), py::arg("nsamples")=100, py::arg("covar")=true) + )", py::arg("pcavar"), py::arg("nsamples")=100, py::arg("ftype")=false, + py::arg("covar")=true) .def("setCovarH5Compress", &BasisClasses::BiorthBasis::setCovarH5Compress, R"( Set the HDF5 compression level for covariance storage in HDF5. The Szip @@ -2498,6 +2401,9 @@ void BasisFactoryClasses(py::module &m) enable (true) or disable (false) the covariance computation nsamples : int number of time partitions to use for covariance computation + ftype: bool + if true, use float32 for covariance storage; if false, + use float64 (default: false) covar: bool if true, compute and save covariance to the HDF5 file; if false, save mean and variance vectors only (default: false) @@ -2512,7 +2418,8 @@ void BasisFactoryClasses(py::module &m) because the number of basis functions can be large. To save disk space, covariance computation is disabled by default. The user can enable it by calling this member function with covar=True. - )", py::arg("pcavar"), py::arg("nsamples")=100, py::arg("covar")=false) + )", py::arg("pcavar"), py::arg("nsamples")=100, py::arg("ftype")=false, + py::arg("covar")=false) .def("index1D", &BasisClasses::Cube::index1D, R"( Returns a flattened 1-d index into the arrays and matrices returned by the @@ -2962,7 +2869,7 @@ void BasisFactoryClasses(py::module &m) py::arg("ps"), py::arg("basiscoef"), py::arg("func"), py::arg("nout")=0); - py::class_> + py::class_> (m, "CovarianceReader") .def(py::init(), R"( @@ -2981,7 +2888,7 @@ void BasisFactoryClasses(py::module &m) the new instance )", py::arg("filename"), py::arg("stride")=1) - .def("Times", &BasisClasses::CovarianceReader::Times, + .def("Times", &BasisClasses::SubsampleCovariance::Times, R"( Get the list of evaluation times @@ -2995,23 +2902,7 @@ void BasisFactoryClasses(py::module &m) a list of evaluation times )") .def("getCoefCovariance", static_cast>> - (BasisClasses::CovarianceReader::*)(unsigned)>(&BasisClasses::CovarianceReader::getCoefCovariance), - R"( - Get the covariance matrices for the basis coefficients - - Parameters - ---------- - index : int - the time index - Returns - ------- - list(list(tuple(numpy.ndarray, numpy.ndarray))) - list of partitioned coefficients and their covariance matrices for - each subsample - )", - py::arg("index")) - .def("getCoefCovariance", static_cast>> - (BasisClasses::CovarianceReader::*)(double)>(&BasisClasses::CovarianceReader::getCoefCovariance), + (BasisClasses::SubsampleCovariance::*)(double)>(&BasisClasses::SubsampleCovariance::getCoefCovariance), R"( Get the covariance matrices for the basis coefficients @@ -3027,23 +2918,7 @@ void BasisFactoryClasses(py::module &m) each subsample. The returns are complex-valued. )", py::arg("time")) - .def("getCovarSamples", - py::overload_cast(&BasisClasses::CovarianceReader::getCovarSamples), - R"( - Get sample counts for the covariance computation - - Parameters - ---------- - index : unsigned int - the time index - - Returns - ------- - tuple(numpy.ndarray, numpy.ndarray) - sample counts and masses for the covariance computation - )", - py::arg("index")) - .def("getCovarSamples", py::overload_cast(&BasisClasses::CovarianceReader::getCovarSamples), + .def("getCovarSamples", &BasisClasses::SubsampleCovariance::getCovarSamples, R"( Get sample counts for the covariance computation @@ -3058,7 +2933,7 @@ void BasisFactoryClasses(py::module &m) sample counts and masses for the covariance computation )", py::arg("time")) - .def("basisIDname", &BasisClasses::CovarianceReader::basisIDname, + .def("basisIDname", &BasisClasses::SubsampleCovariance::basisIDname, R"( Get the basis ID name diff --git a/src/Cube.H b/src/Cube.H index 657a192e8..5a8b2a3ed 100644 --- a/src/Cube.H +++ b/src/Cube.H @@ -275,9 +275,6 @@ private: public: - //! Id string - std::string id; - //! Constructor Cube(Component* c0, const YAML::Node& conf); @@ -297,7 +294,7 @@ public: elements for matrices for analysis. The default implementation returns empty vectors. This parallels the implementation in expui. */ - CovarElement getSubsample(); + CovarData getSubsample(); //! Indicate whether or not subsampling is supported virtual bool hasSubsample() { return true; } diff --git a/src/Cube.cc b/src/Cube.cc index 62b1685a8..4b1e71779 100644 --- a/src/Cube.cc +++ b/src/Cube.cc @@ -134,14 +134,10 @@ void Cube::initialize(void) if (conf["wrap" ]) wrap = conf["wrap" ].as(); if (conf["nint"]) { - nint = conf["nint" ].as(); + nint = conf["nint"].as(); if (nint>0) computeSubsample = true; } - if (conf["subsampleFloat"]) { - floatType = conf["subsampleFloat"].as(); - } - if (conf["samplesz"]) { sampT = conf["samplesz"].as(); } @@ -1028,7 +1024,7 @@ void Cube::writeCovarH5Params(HighFive::File& file) } -PotAccel::CovarElement Cube::getSubsample() +PotAccel::CovarData Cube::getSubsample() { using covTuple = std::tuple; @@ -1046,6 +1042,6 @@ PotAccel::CovarElement Cube::getSubsample() covar[T][0] = std::make_tuple(meanV[T], Eigen::MatrixXcd::Zero(0,0)); } - return {sampleMasses, sampleCounts, covar}; + return {sampleCounts, sampleMasses, covar}; } diff --git a/src/OutSample.H b/src/OutSample.H index 7b3ca6dee..eb7407641 100644 --- a/src/OutSample.H +++ b/src/OutSample.H @@ -1,7 +1,7 @@ #ifndef _OutSample_H #define _OutSample_H -#include "OutSample.H" +#include "Covariance.H" /** @file OutSample.H @@ -28,12 +28,20 @@ private: std::string filename; double prev = -std::numeric_limits::max(); Component *tcomp; + unsigned level=5; + unsigned chunksize=1024*1024; // 1 MB + bool shuffle=true; + bool szip=false; + bool floatType = false; void initialize(void); //! Valid keys for YAML configurations static const std::set valid_keys; + //! Coefficient covariance storage instance + std::shared_ptr covarStore; + public: //! Constructor diff --git a/src/OutSample.cc b/src/OutSample.cc index 02c7035ac..52c4ebd2e 100644 --- a/src/OutSample.cc +++ b/src/OutSample.cc @@ -8,8 +8,13 @@ const std::set OutSample::valid_keys = { + "floatType", "filename", - "name" + "name", + "level", + "chunksize", + "compress", + "szip" }; OutSample::OutSample(const YAML::Node& conf) : Output(conf) @@ -23,6 +28,16 @@ OutSample::OutSample(const YAML::Node& conf) : Output(conf) throw std::runtime_error("OutSample: no subsampling for this force"); } + covarStore = std::make_shared + ([this](HighFive::File& file){this->tcomp->force->writeCovarH5Params(file);}, + this->tcomp->force->id, floatType, this->tcomp->force->FullCovar()); + + covarStore->setCovarH5Compress(level, chunksize, shuffle, szip); + + if (myid==0) + std::cout << "---- OutSample: writing subsample coefficients for component '" + << tcomp->name << "' to file " << filename << " for force '" + << tcomp->force->id << "'" << std::endl; } void OutSample::initialize() @@ -54,6 +69,26 @@ void OutSample::initialize() filename = outdir + "outcoef." + tcomp->name + "." + runtag; } + if (conf["floatType"]) { + floatType = conf["floatType"].as(); + } + + if (conf["level"]) { + level = conf["level"].as(); + } + + if (conf["chunksize"]) { + chunksize = conf["chunksize"].as(); + } + + if (conf["shuffle"]) { + shuffle = conf["shuffle"].as(); + } + + if (conf["szip"]) { + szip = conf["szip"].as(); + } + } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing parameters in OutSample: " @@ -81,5 +116,7 @@ void OutSample::Run(int n, int mstep, bool last) prev = tnow; // Write the subsample data - tcomp->force->writeSubsample(); + auto elem = tcomp->force->getSubsample(); + covarStore->writeCoefCovariance(tcomp->name, runtag, elem, tnow); + } diff --git a/src/PotAccel.H b/src/PotAccel.H index efdc9c0ea..6f1f56460 100644 --- a/src/PotAccel.H +++ b/src/PotAccel.H @@ -114,45 +114,9 @@ protected: //! @name Subsampling //@{ - /** Return a subsample of the basis coefficients and covariance - elements for matrices for PCA analysis. The default - implementation returns empty vectors. */ - virtual CovarElement getSubsample() - { - Eigen::VectorXd mass; - Eigen::VectorXi count; - std::vector>> covar; - - return std::make_tuple(mass, count, covar); - } - //! Reset subsample data virtual void subsampleReset() { subsampleComputed = false; } - - //! Extend the subsample data - virtual void extendSubsample(const std::string& fname, double time); - - //! Write basis-specific parameters to HDF5 covariance file - virtual void writeCovarH5Params(HighFive::File& file) - { - if (myid==0) - std::cout << "Basis::writeCovarH5Params: " - << "not implemented for basis type <" << id << ">" - << std::endl; - } - - //! Write basis-specific parameters to HDF5 covariance file - virtual unsigned writeCovarH5(HighFive::Group& group, - unsigned count, double time); - - //! Variance file versioning (Version 1.0 has more granularity and - //! poor compression) - inline static const std::string CovarianceFileVersion = "1.1"; - - //! Use 32-bit floating point type in covariance file - bool floatType = false; - + //! Save covariance subsamples to file bool fullCovar = false; @@ -167,15 +131,6 @@ protected: //@} - //@{ - //! HDF5 compression settings - unsigned H5compress = 5; - unsigned H5chunk = 1024*1024; // (1 MB) - bool H5shuffle = true; - bool H5szip = false; - //@} - - public: //! For timing data @@ -186,7 +141,7 @@ public: }; //! Id string - string id; + std::string id; //! Possible geometries enum Geometry {sphere, cylinder, cube, slab, table, other}; @@ -266,9 +221,6 @@ public: //@name Subsampling methods //@{ - //! Write the subsample to file - virtual void writeSubsample(); - //! Set flag to compute subsample data virtual void subsampleRequest() { requestSubsample = true; } @@ -277,6 +229,15 @@ public: //! Indicate whether or not subsampling is supported virtual bool hasSubsample() { return false; } + + //! Write basis-specific parameters to HDF5 covariance file + virtual void writeCovarH5Params(HighFive::File& file) + { + if (myid==0) + std::cout << "Basis::writeCovarH5Params: " + << "not implemented for basis type <" << id << ">" + << std::endl; + } //@} //@{ @@ -285,14 +246,21 @@ public: Eigen::VectorXd sampleMasses; //@} - //! HDF5 compression settings - void setCovarH5Compress(unsigned level, unsigned chunksize, bool shuffle, bool szip=false) + //@{ + //! Type for coefficient and covariance matrix tuple + + using CoefCovarType = std::tuple; + + using CovarData = std::tuple>>; + + virtual CovarData getSubsample() { - H5compress = level; - H5chunk = chunksize; - H5shuffle = shuffle; - H5szip = szip; + // Must be overriden; base implementation throws error + throw std::runtime_error("PotAccel::getSubsample: Not implemented for this basis"); } + //@} + /** Update the multi time step force algorithm when moving particle i from level cur to level @@ -331,6 +299,9 @@ public: //! Check whether coeffcients are not available bool NoCoefs() { return play_back and !play_cnew; } + //! Get fullCovar flag + bool FullCovar() { return fullCovar; } + //! For timing the fraction spent in threaded force methods //@{ void print_timings(const string&, TList& tlist); @@ -347,10 +318,6 @@ public: //! Get unaccounted keys std::set unmatched() { return current_keys; } - - //! Round time key to emulated fixed-point arithmetic - double roundTime(double time); - }; //! Helper class to pass info to threaded member diff --git a/src/PotAccel.cc b/src/PotAccel.cc index 85acb92fb..c1d6f9d9f 100644 --- a/src/PotAccel.cc +++ b/src/PotAccel.cc @@ -368,328 +368,3 @@ bool ltdub(const pair& A, return A.first < B.first; } - -void PotAccel::writeSubsample() -{ - // Check that subsample data is available - // - if (not hasSubsample()) { - std::cout << "PotAccel::writeSubsample: " - << "subsample data is not available. " << std::endl; - return; - } - - // Only root process writes - // - if (myid) return; - - // Round time - // - double time = roundTime(tnow); - - // The H5 filename - // - std::string fname = "subsample." + cC->name + "." + runtag + ".h5"; - - // Check if file exists? - // - try { - // Open the HDF5 file in read-write mode, creating if it doesn't - // exist - HighFive::File file(fname, - HighFive::File::ReadWrite | - HighFive::File::Create); - - // Check for version string - std::string path = "CovarianceFileVersion"; - - // Check for valid HDF file by attribute - if (file.hasAttribute(path)) { - extendSubsample(fname, time); - return; - } - - // Write the Version string - // - file.createAttribute("CovarianceFileVersion", HighFive::DataSpace::From(CovarianceFileVersion)).write(CovarianceFileVersion); - - // Write the basis identifier string - // - file.createAttribute("BasisID", HighFive::DataSpace::From(id)).write(id); - - // Write the data type size - // - int sz = 8; if (floatType) sz = 4; - file.createAttribute("FloatSize", HighFive::DataSpace::From(sz)).write(sz); - - // Write the specific parameters - // - writeCovarH5Params(file); - - // Group count variable - // - unsigned count = 0; - HighFive::DataSet dataset = file.createDataSet("count", count); - - // Create a new group for coefficient snapshots - // - HighFive::Group group = file.createGroup("snapshots"); - - // Write the coefficients - // - count = writeCovarH5(group, count, time); - - // Update the count - // - dataset.write(count); - - } catch (const HighFive::Exception& err) { - // Handle HighFive specific errors (e.g., file not found) - throw std::runtime_error - (std::string("PotAccel::writeCoefCovariance HighFive Error: ") + err.what()); - } catch (const std::exception& err) { - // Handle other general exceptions - throw std::runtime_error - (std::string("PotAccel::writeCoefCovariance Error: ") + err.what()); - } - - // Reset subsample data - subsampleReset(); -} - -void PotAccel::extendSubsample(const std::string& fname, double time) -{ - try { - // Open an hdf5 file - // - HighFive::File file(fname, HighFive::File::ReadWrite); - - // Get the dataset - HighFive::DataSet dataset = file.getDataSet("count"); - - unsigned count; - dataset.read(count); - - HighFive::Group group = file.getGroup("snapshots"); - - // Write the coefficients - // - count = writeCovarH5(group, count, time); - - // Update the count - // - dataset.write(count); - - } catch (HighFive::Exception& err) { - throw std::runtime_error - (std::string("PotAccel::extendSubsample: HighFive error: ") + err.what()); - } -} - -double PotAccel::roundTime(double time) -{ - // Eight decimal places should be enough here... - const double multiplier = 1.0e+08; // std::pow(10.0, 8); - return std::floor(time * multiplier + 0.5) / multiplier; -} - -unsigned PotAccel::writeCovarH5(HighFive::Group& snaps, unsigned count, double time) -{ - std::ostringstream stim; - stim << std::setw(8) << std::setfill('0') << std::right << count++; - - // Get the subsample data - // - auto [sampleMasses, sampleCounts, covdata] = getSubsample(); - - // Make a new group for this time - // - HighFive::Group stanza = snaps.createGroup(stim.str()); - - // Add a time attribute - // - time = roundTime(time); - stanza.createAttribute("Time", HighFive::DataSpace::From(time)).write(time); - - // Enable compression - // - auto dcpl1 = HighFive::DataSetCreateProps{}; // sample stats - auto dcpl2 = HighFive::DataSetCreateProps{}; // coefficients - auto dcpl3 = HighFive::DataSetCreateProps{}; // covariance - - // Properties for sample stats - if (H5compress) { - unsigned int csz = sampleCounts.size(); - dcpl1.add(HighFive::Chunking({csz, 1})); - if (H5shuffle) dcpl1.add(HighFive::Shuffle()); - dcpl1.add(HighFive::Deflate(H5compress)); - } - - // Add the sample statistics - // - HighFive::DataSet s1data = stanza.createDataSet("sampleCounts", sampleCounts, dcpl1); - HighFive::DataSet s2data = stanza.createDataSet("sampleMasses", sampleMasses, dcpl1); - - - // Number of samples - // - unsigned sampleSize = covdata.size(); - unsigned ltot = covdata[0].size(); - unsigned nmax = std::get<0>(covdata[0][0]).rows(); - unsigned diagonalSize = nmax; - - // Save variance or full covariance - if (fullCovar) diagonalSize = nmax*(nmax + 1)/2; - - // Add data dimensions - // - stanza.createAttribute - ("sampleSize", HighFive::DataSpace::From(sampleSize)).write(sampleSize); - - stanza.createAttribute - ("angularSize", HighFive::DataSpace::From(ltot)).write(ltot); - - stanza.createAttribute - ("rankSize", HighFive::DataSpace::From(nmax)).write(nmax); - - int icov = fullCovar ? 1 : 0; - stanza.createAttribute - ("fullCovar", HighFive::DataSpace::From(icov)).write(icov); - - if (H5compress) { - // Szip parameters - const int options_mask = H5_SZIP_NN_OPTION_MASK; - const int pixels_per_block = 8; - - // Properties for coefficients - // - unsigned int csz2 = nmax * ltot * sampleSize; - HighFive::Chunking data_dims2{std::min(csz2, H5chunk), 1}; - - dcpl2.add(data_dims2); - if (H5shuffle) dcpl2.add(HighFive::Shuffle()); - if (H5szip) { - dcpl2.add(HighFive::Szip(options_mask, pixels_per_block)); - } else { - dcpl2.add(HighFive::Deflate(H5compress)); - } - - // Properties for covariance - // - unsigned int csz3 = ltot * diagonalSize * sampleSize; - HighFive::Chunking data_dims3{std::min(csz3, H5chunk), 1}; - - dcpl3.add(data_dims3); - if (H5shuffle) dcpl3.add(HighFive::Shuffle()); - if (H5szip) { - dcpl3.add(HighFive::Szip(options_mask, pixels_per_block)); - } else { - dcpl3.add(HighFive::Deflate(H5compress)); - } - } - - // Pack the coefficient data - // - if (floatType) { - // Create a vector of doubles for the real and imaginary parts - Eigen::VectorXf real_part(nmax*ltot*sampleSize); - Eigen::VectorXf imag_part(nmax*ltot*sampleSize); - - for (size_t T=0, c=0; T(covdata[T][l])(n)); - imag_part(c) = std::imag(std::get<0>(covdata[T][l])(n)); - } - } - } - - // Create two separate, compressed datasets - stanza.createDataSet("coefficients_real", real_part, dcpl2); - stanza.createDataSet("coefficients_imag", imag_part, dcpl2); - - if (std::get<1>(covdata[0][0]).size()) { - - real_part.resize(ltot*diagonalSize*sampleSize); - imag_part.resize(ltot*diagonalSize*sampleSize); - - for (size_t T=0, c=0; T(covdata[T][l])(n1, n2)); - imag_part(c) = std::imag(std::get<1>(covdata[T][l])(n1, n2)); - } - } - // Pack the diagonal only - // - else { - real_part(c ) = std::real(std::get<1>(covdata[T][l])(n1, n1)); - imag_part(c++) = std::imag(std::get<1>(covdata[T][l])(n1, n1)); - } - } - } - } - // Create two separate, compressed datasets - stanza.createDataSet("covariance_real", real_part, dcpl3); - stanza.createDataSet("covariance_imag", imag_part, dcpl3); - } - - } else { - Eigen::VectorXd real_part(ltot*nmax*sampleSize); - Eigen::VectorXd imag_part(ltot*nmax*sampleSize); - - for (size_t T=0, c=0; T(covdata[T][l])(n)); - imag_part(c) = std::imag(std::get<0>(covdata[T][l])(n)); - } - } - } - - // Create two separate, compressed datasets - // - stanza.createDataSet("coefficients_real", real_part, dcpl2); - stanza.createDataSet("coefficients_imag", imag_part, dcpl2); - - if (std::get<1>(covdata[0][0]).size()) { - - real_part.resize(ltot*diagonalSize*sampleSize); - imag_part.resize(ltot*diagonalSize*sampleSize); - - for (size_t T=0, c=0; T(covdata[T][l])(n1, n2)); - imag_part(c) = std::imag(std::get<1>(covdata[T][l])(n1, n2)); - } - } - // Pack the diagonal only - // - else { - real_part(c ) = std::real(std::get<1>(covdata[T][l])(n1, n1)); - imag_part(c++) = std::imag(std::get<1>(covdata[T][l])(n1, n1)); - } - } - } - } - - // Create two separate, compressed datasets - // - stanza.createDataSet("covariance_real", real_part, dcpl3); - stanza.createDataSet("covariance_imag", imag_part, dcpl3); - } - } - // END: sample loop - - return count; -} - diff --git a/src/SphericalBasis.cc b/src/SphericalBasis.cc index ffa632971..51e0fcaf2 100644 --- a/src/SphericalBasis.cc +++ b/src/SphericalBasis.cc @@ -133,10 +133,6 @@ SphericalBasis::SphericalBasis(Component* c0, const YAML::Node& conf, MixtureBas if (conf["noise_model_file"]) noise_model_file = conf["noise_model_file"].as(); - if (conf["subsampleFloat"]) { - floatType = conf["subsampleFloat"].as(); - } - if (conf["seedN"]) seedN = conf["seedN"].as(); if (conf["ssfrac"]) { From 3799ed277161fe0a9961d79b38dab2d1dd00d099 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 2 Dec 2025 13:53:21 -0500 Subject: [PATCH 05/38] Added writeCoefCovariance for Cylindrical --- expui/BiorthBasis.H | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index e01f275b6..77dc61018 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -501,6 +501,7 @@ namespace BasisClasses } } + //! Request and initialize covariance computation virtual void enableCoefCovariance(bool pcavar_in, int sampT_in=100, bool ftype=false, bool covar_in=true) { @@ -513,6 +514,7 @@ namespace BasisClasses ([this](HighFive::File& file){this->writeCovarH5Params(file);}, BasisID, ftype, covar); } } + //! Create and the coefficients from a function callback with the //! provided time value. Set potential=true if function is a //! potential field. Density is assumed if false (default). @@ -1084,6 +1086,18 @@ namespace BasisClasses return ret; } + //! Write coefficient covariance data to an HDF5 file + virtual void writeCoefCovariance(const std::string& compname, const std::string& runtag, double time=0.0) + { + if (covarStore) { + auto covarData = sl->getCoefCovariance(); + SubsampleCovariance::CovarData elements = std::make_tuple(sampleCounts, sampleMasses, covarData); + covarStore->writeCoefCovariance(compname, runtag, elements, time); + } else { + throw std::runtime_error("Spherical::writeCoefCovariance: covariance storage not initialized"); + } + } + virtual void enableCoefCovariance(bool pcavar_in, int sampT_in=100, bool ftype=false, bool covar_in=true) { From 825ece21c8bbbe67a382b1277914445afe8883e9 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 2 Dec 2025 13:53:48 -0500 Subject: [PATCH 06/38] Comment changes only --- include/Covariance.H | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/include/Covariance.H b/include/Covariance.H index 1fdb012dd..d9636b66f 100644 --- a/include/Covariance.H +++ b/include/Covariance.H @@ -16,20 +16,6 @@ #include #include -/* -#include "DiskDensityFunc.H" -#include "ParticleReader.H" -#include "Coefficients.H" -#include "BiorthBess.H" -#include "BasisFactory.H" -#include "BiorthCube.H" -#include "sltableMP2.H" -#include "SLGridMP2.H" -#include "YamlCheck.H" -#include "BiorthCyl.H" -#include "EmpCylSL.H" -*/ - #include "localmpi.H" #include "exputils.H" @@ -63,7 +49,6 @@ namespace BasisClasses //! Use 32-bit float storage for covariance? bool floatType = false; - //! Fixed-point multiplier for time mapping. Eight decimal places //! should be enough here... const double multiplier = 1.0e+08; From 8bac1183eaa161d3d2a048d179059413e18e5b53 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 2 Dec 2025 13:54:07 -0500 Subject: [PATCH 07/38] Comment changes only --- src/Cube.H | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/Cube.H b/src/Cube.H index 5a8b2a3ed..bf030426d 100644 --- a/src/Cube.H +++ b/src/Cube.H @@ -47,6 +47,17 @@ dimensions at once and planes, axes, and 1d which compute wave vectors by plane. This may be extended to specific algorithms in the future. This does nothing for the CPU-only computation. + + @param wrap boolean to enforce periodic wrapping (default: true) + + @param nint integer number of steps between subsample computations + (default: 0, no subsampling) + + @param samplesz integer number of particles in subsample (default: + 100) + + @param subsampleFloat boolean to use single precision in HDF5 + files (default: false) */ class Cube : public PotAccel { From 1e5dcc15fc10256fb3b9ba7004d9f843fdd54442 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Tue, 2 Dec 2025 14:04:13 -0500 Subject: [PATCH 08/38] Update src/Cube.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- src/Cube.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Cube.cc b/src/Cube.cc index 34ad0cd15..4a82b8acb 100644 --- a/src/Cube.cc +++ b/src/Cube.cc @@ -410,7 +410,7 @@ void Cube::determine_coefficients(void) // Determine whether or not to compute a subsample if (mstep==0 or mstep==std::numeric_limits::max()) { - if (nint>0 and this_step % nint == 0) { + if (nint>0 && this_step % nint == 0) { if (tnow > last) { requestSubsample = true; last = tnow; From 311c70dbd075080aca569a2799aca621f7644700 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 2 Dec 2025 19:04:57 +0000 Subject: [PATCH 09/38] Initial plan From 1ab52bc055d994daef88fed68932d22cddfa924a Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 2 Dec 2025 19:05:13 +0000 Subject: [PATCH 10/38] Initial plan From 43f4d728fff7225704271152782d1ad46c08789d Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Tue, 2 Dec 2025 14:05:25 -0500 Subject: [PATCH 11/38] Update expui/Covariance.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- expui/Covariance.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/Covariance.cc b/expui/Covariance.cc index b1f5af840..1079d7a40 100644 --- a/expui/Covariance.cc +++ b/expui/Covariance.cc @@ -379,7 +379,7 @@ namespace BasisClasses file.getAttribute("nmaxz").read(nmaxz); ltot = (2*nmaxx + 1) * (2*nmaxy + 1) * (2*nmaxz + 1); } else { - throw std::runtime_error(std::string("Covariance: unknown or unimplemented covariance for basis type, ") + BasisID); + throw std::runtime_error(std::string("SubsampleCovariance: unknown or unimplemented covariance for basis type, ") + BasisID); } // Group count variable From e0f66bd4190ad4d7cf814d5d7eb38cf9075664c1 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 2 Dec 2025 19:05:31 +0000 Subject: [PATCH 12/38] Initial plan From ed77d7b5e378cacc6331de45a86874a601f8812d Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 2 Dec 2025 19:07:59 +0000 Subject: [PATCH 13/38] Fix error message class name from CovarianceReader to SubsampleCovariance Co-authored-by: The9Cat <25960766+The9Cat@users.noreply.github.com> --- expui/Covariance.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/expui/Covariance.cc b/expui/Covariance.cc index b1f5af840..f6a0eb9b1 100644 --- a/expui/Covariance.cc +++ b/expui/Covariance.cc @@ -334,11 +334,11 @@ namespace BasisClasses file.getAttribute("CovarianceFileVersion").read(version); // Check for alpha version if (version == std::string("1.0")) { - throw std::runtime_error("CovarianceReader: this is an early alpha test version. Please remake your files"); + throw std::runtime_error("SubsampleCovariance: this is an early alpha test version. Please remake your files"); } // Test for current version if (version != std::string("1.1")) { - throw std::runtime_error(std::string("CovarianceReader: unsupported file version, ") + version); + throw std::runtime_error(std::string("SubsampleCovariance: unsupported file version, ") + version); } // Read the basis identifier string @@ -350,7 +350,7 @@ namespace BasisClasses file.getAttribute("FloatSize").read(sz); if (sz != 4 and sz != 8) { std::ostringstream sout; - sout << "CovarianceReader: unsupported float size, " << sz; + sout << "SubsampleCovariance: unsupported float size, " << sz; throw std::runtime_error(sout.str()); } From 539a19ae875a516988ee8f9a2d4c0c1e0d9bfb17 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 2 Dec 2025 19:09:15 +0000 Subject: [PATCH 14/38] Fix class name in error messages: CovarianceReader -> SubsampleCovariance Co-authored-by: The9Cat <25960766+The9Cat@users.noreply.github.com> --- expui/Covariance.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/expui/Covariance.cc b/expui/Covariance.cc index b1f5af840..f6a0eb9b1 100644 --- a/expui/Covariance.cc +++ b/expui/Covariance.cc @@ -334,11 +334,11 @@ namespace BasisClasses file.getAttribute("CovarianceFileVersion").read(version); // Check for alpha version if (version == std::string("1.0")) { - throw std::runtime_error("CovarianceReader: this is an early alpha test version. Please remake your files"); + throw std::runtime_error("SubsampleCovariance: this is an early alpha test version. Please remake your files"); } // Test for current version if (version != std::string("1.1")) { - throw std::runtime_error(std::string("CovarianceReader: unsupported file version, ") + version); + throw std::runtime_error(std::string("SubsampleCovariance: unsupported file version, ") + version); } // Read the basis identifier string @@ -350,7 +350,7 @@ namespace BasisClasses file.getAttribute("FloatSize").read(sz); if (sz != 4 and sz != 8) { std::ostringstream sout; - sout << "CovarianceReader: unsupported float size, " << sz; + sout << "SubsampleCovariance: unsupported float size, " << sz; throw std::runtime_error(sout.str()); } From 1c17c8766f9dd0ebab74902ca21c3a79b5db0197 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 2 Dec 2025 19:09:44 +0000 Subject: [PATCH 15/38] Fix additional class name references for consistency Co-authored-by: The9Cat <25960766+The9Cat@users.noreply.github.com> --- expui/Covariance.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/expui/Covariance.cc b/expui/Covariance.cc index f6a0eb9b1..c10449081 100644 --- a/expui/Covariance.cc +++ b/expui/Covariance.cc @@ -362,7 +362,7 @@ namespace BasisClasses // Currently implemented cylindrical types const std::set cylinderType = {"Cylindrical"}; - std::cout << "Covariance: reading basis type " << BasisID << std::endl; + std::cout << "SubsampleCovariance: reading basis type " << BasisID << std::endl; if (sphereType.find(BasisID) != sphereType.end()) { file.getAttribute("lmax").read(lmax); @@ -379,7 +379,7 @@ namespace BasisClasses file.getAttribute("nmaxz").read(nmaxz); ltot = (2*nmaxx + 1) * (2*nmaxy + 1) * (2*nmaxz + 1); } else { - throw std::runtime_error(std::string("Covariance: unknown or unimplemented covariance for basis type, ") + BasisID); + throw std::runtime_error(std::string("SubsampleCovariance: unknown or unimplemented covariance for basis type, ") + BasisID); } // Group count variable From ac06b304bf09fb64a06f5117e10125590295543b Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 2 Dec 2025 19:10:31 +0000 Subject: [PATCH 16/38] Complete review checks Co-authored-by: The9Cat <25960766+The9Cat@users.noreply.github.com> --- _codeql_detected_source_root | 1 + 1 file changed, 1 insertion(+) create mode 120000 _codeql_detected_source_root diff --git a/_codeql_detected_source_root b/_codeql_detected_source_root new file mode 120000 index 000000000..945c9b46d --- /dev/null +++ b/_codeql_detected_source_root @@ -0,0 +1 @@ +. \ No newline at end of file From 0d46fd6bcd38aed6a5a3d604517f2cc643007bd5 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 2 Dec 2025 19:11:05 +0000 Subject: [PATCH 17/38] Complete class name consistency fixes Co-authored-by: The9Cat <25960766+The9Cat@users.noreply.github.com> --- _codeql_detected_source_root | 1 + 1 file changed, 1 insertion(+) create mode 120000 _codeql_detected_source_root diff --git a/_codeql_detected_source_root b/_codeql_detected_source_root new file mode 120000 index 000000000..945c9b46d --- /dev/null +++ b/_codeql_detected_source_root @@ -0,0 +1 @@ +. \ No newline at end of file From 3cf02b9820d51bdab7eeaffd4819b6afea6abed9 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 2 Dec 2025 19:13:18 +0000 Subject: [PATCH 18/38] Fix class name in error messages to SubsampleCovariance Co-authored-by: The9Cat <25960766+The9Cat@users.noreply.github.com> --- expui/Covariance.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/expui/Covariance.cc b/expui/Covariance.cc index 1079d7a40..c10449081 100644 --- a/expui/Covariance.cc +++ b/expui/Covariance.cc @@ -334,11 +334,11 @@ namespace BasisClasses file.getAttribute("CovarianceFileVersion").read(version); // Check for alpha version if (version == std::string("1.0")) { - throw std::runtime_error("CovarianceReader: this is an early alpha test version. Please remake your files"); + throw std::runtime_error("SubsampleCovariance: this is an early alpha test version. Please remake your files"); } // Test for current version if (version != std::string("1.1")) { - throw std::runtime_error(std::string("CovarianceReader: unsupported file version, ") + version); + throw std::runtime_error(std::string("SubsampleCovariance: unsupported file version, ") + version); } // Read the basis identifier string @@ -350,7 +350,7 @@ namespace BasisClasses file.getAttribute("FloatSize").read(sz); if (sz != 4 and sz != 8) { std::ostringstream sout; - sout << "CovarianceReader: unsupported float size, " << sz; + sout << "SubsampleCovariance: unsupported float size, " << sz; throw std::runtime_error(sout.str()); } @@ -362,7 +362,7 @@ namespace BasisClasses // Currently implemented cylindrical types const std::set cylinderType = {"Cylindrical"}; - std::cout << "Covariance: reading basis type " << BasisID << std::endl; + std::cout << "SubsampleCovariance: reading basis type " << BasisID << std::endl; if (sphereType.find(BasisID) != sphereType.end()) { file.getAttribute("lmax").read(lmax); From 8ea07a7b258011fb62f7c70b6c0136e2c9b67971 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 2 Dec 2025 19:14:01 +0000 Subject: [PATCH 19/38] Complete feedback implementation Co-authored-by: The9Cat <25960766+The9Cat@users.noreply.github.com> --- _codeql_detected_source_root | 1 + 1 file changed, 1 insertion(+) create mode 120000 _codeql_detected_source_root diff --git a/_codeql_detected_source_root b/_codeql_detected_source_root new file mode 120000 index 000000000..945c9b46d --- /dev/null +++ b/_codeql_detected_source_root @@ -0,0 +1 @@ +. \ No newline at end of file From 0e1f720a46225870f7157a9231e6c11c181cc4ef Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 2 Dec 2025 15:14:07 -0500 Subject: [PATCH 20/38] Name stopping with OutCoef --- src/OutSample.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OutSample.cc b/src/OutSample.cc index 52c4ebd2e..b648b0220 100644 --- a/src/OutSample.cc +++ b/src/OutSample.cc @@ -66,7 +66,7 @@ void OutSample::initialize() if (conf["filename"]) { filename = outdir + conf["filename"].as(); } else { - filename = outdir + "outcoef." + tcomp->name + "." + runtag; + filename = outdir + "coefcovar." + tcomp->name + "." + runtag; } if (conf["floatType"]) { From b30b9918e76f4ff5d13995e1e0c0706695764c2d Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 2 Dec 2025 19:43:47 -0500 Subject: [PATCH 21/38] Implemented CUDA version of subsampling but not tested --- src/Cube.H | 11 ++- src/cudaCube.cu | 205 +++++++++++++++++++++++++++++++++++++++++++++--- 2 files changed, 204 insertions(+), 12 deletions(-) diff --git a/src/Cube.H b/src/Cube.H index bf030426d..68cae68f9 100644 --- a/src/Cube.H +++ b/src/Cube.H @@ -135,7 +135,16 @@ private: thrust::device_vector> dw_coef; thrust::device_vector> df_coef; - void resize_coefs(int N, int osize, int gridSize, int stride); + thrust::device_vector> dN_tvar; + thrust::device_vector> dc_tvar; + thrust::device_vector> dw_tvar; + + thrust::device_vector u_d; + + std::vector>> T_covr; + + void resize_coefs(int N, int osize, int gridSize, + int stride, int sampT); }; //! A storage instance diff --git a/src/cudaCube.cu b/src/cudaCube.cu index e97aadba5..f3a387f74 100644 --- a/src/cudaCube.cu +++ b/src/cudaCube.cu @@ -94,6 +94,10 @@ void Cube::cuda_initialize() // byPlanes = true; + // Turn off full covariance + // + fullCovar = false; + // Make method string lower case // std::transform(cuMethod.begin(), cuMethod.end(), cuMethod.begin(), @@ -163,7 +167,7 @@ void Cube::initialize_constants() __global__ void coefKernelCube (dArray P, dArray I, dArray coef, - int stride, PII lohi) + dArray used, int stride, PII lohi) { // Thread ID // @@ -189,6 +193,8 @@ __global__ void coefKernelCube cuFP_t pos[3] = {p.pos[0], p.pos[1], p.pos[2]}; cuFP_t mm = p.mass; + used._v[i] = mm; // Mark particle as used + // Enforce periodicity in the unit cube // if (cubeWrap) { @@ -269,7 +275,7 @@ __global__ void coefKernelCube __global__ void coefKernelCubeX (dArray P, dArray I, dArray coef, - int jj, int kk, int stride, PII lohi) + dArray used, int jj, int kk, int stride, PII lohi) { // Thread ID // @@ -295,6 +301,8 @@ __global__ void coefKernelCubeX cuFP_t pos[3] = {p.pos[0], p.pos[1], p.pos[2]}; cuFP_t mm = p.mass; + used._v[i] = mm; // Mark particle as used + // Enforce periodicity in the unit cube // if (cubeWrap) { @@ -519,7 +527,8 @@ public: } }; -void Cube::cudaStorage::resize_coefs(int N, int osize, int gridSize, int stride) +void Cube::cudaStorage::resize_coefs(int N, int osize, int gridSize, int stride, + int sampT) { // Reserve space for coefficient reduction // @@ -534,6 +543,15 @@ void Cube::cudaStorage::resize_coefs(int N, int osize, int gridSize, int stride) dN_coef.resize(osize*N); dc_coef.resize(osize*gridSize); dw_coef.resize(osize); // This will stay fixed + + if (sampT) { + T_covr.resize(sampT); + for (int T=0; TCudaGetLevelRange(mlevel, mlevel), cur; @@ -660,11 +682,18 @@ void Cube::determine_coefficients_cuda() // int sMemSize = BLOCK_SIZE * sizeof(CmplxT); + std::vector::iterator> bm; + if (requestSubsample) { + for (int T=0; Tstream>>> (toKernel(cs->cuda_particles), toKernel(cs->indx1), - toKernel(cuS.dN_coef), jj, kk, stride, cur); + toKernel(cuS.dN_coef), toKernel(cuS.u_d), jj, kk, stride, cur); // Begin the reduction by blocks [perhaps this should use a // stride?] @@ -708,6 +737,77 @@ void Cube::determine_coefficients_cuda() beg, beg, thrust::plus()); thrust::advance(beg, imx); + + if (requestSubsample) { + + int sN = N/sampT; + int nT = sampT; + + if (sN==0) { // Fail-safe underrun + sN = 1; + nT = N; + } + + for (int T=0; T gridSize1*BLOCK_SIZE*stride1) gridSize1++; + */ + + unsigned int gridSize1 = s/BLOCK_SIZE; + if (s > gridSize1*BLOCK_SIZE) gridSize1++; + + // Sum reduce into gridsize1 blocks taking advantage of + // GPU warp structure + // + + // Mean computation only + // + reduceSumS + <<stream>>> + (toKernel(cuS.dc_tvar), toKernel(cuS.dN_tvar), imx, N, k, k+s); + + // Finish the reduction for this order in parallel + // + thrust::counting_iterator index_begin(0); + thrust::counting_iterator index_end(gridSize1*imx); + + thrust::reduce_by_key + ( + thrust::cuda::par.on(cs->stream), + thrust::make_transform_iterator(index_begin, key_functor(gridSize1)), + thrust::make_transform_iterator(index_end, key_functor(gridSize1)), + cuS.dc_tvar.begin(), thrust::make_discard_iterator(), cuS.dw_tvar.begin() + ); + + + thrust::transform(thrust::cuda::par.on(cs->stream), + cuS.dw_tvar.begin(), cuS.dw_tvar.end(), + bm[T], bm[T], thrust::plus()); + + thrust::advance(bm[T], imx); + } + } } // END: y wave number loop } @@ -717,7 +817,7 @@ void Cube::determine_coefficients_cuda() // Adjust cached storage, if necessary // - cuS.resize_coefs(N, osize, gridSize, stride); + cuS.resize_coefs(N, osize, gridSize, stride, sampT); // Compute the coefficient contribution for each order // @@ -725,7 +825,7 @@ void Cube::determine_coefficients_cuda() coefKernelCube<<stream>>> (toKernel(cs->cuda_particles), toKernel(cs->indx1), - toKernel(cuS.dN_coef), stride, cur); + toKernel(cuS.dN_coef), toKernel(cuS.u_d), stride, cur); // Begin the reduction by blocks [perhaps this should use a // stride?] @@ -757,6 +857,77 @@ void Cube::determine_coefficients_cuda() beg, beg, thrust::plus()); thrust::advance(beg, osize); + + if (requestSubsample) { + + int sN = N/sampT; + int nT = sampT; + + if (sN==0) { // Fail-safe underrun + sN = 1; + nT = N; + } + + for (int T=0; T gridSize1*BLOCK_SIZE*stride1) gridSize1++; + */ + + unsigned int gridSize1 = s/BLOCK_SIZE; + if (s > gridSize1*BLOCK_SIZE) gridSize1++; + + // Sum reduce into gridsize1 blocks taking advantage of + // GPU warp structure + // + + // Mean computation only + // + reduceSumS + <<stream>>> + (toKernel(cuS.dc_tvar), toKernel(cuS.dN_tvar), osize, N, k, k+s); + + // Finish the reduction for this order in parallel + // + thrust::counting_iterator index_begin(0); + thrust::counting_iterator index_end(gridSize1*osize); + + thrust::reduce_by_key + ( + thrust::cuda::par.on(cs->stream), + thrust::make_transform_iterator(index_begin, key_functor(gridSize1)), + thrust::make_transform_iterator(index_end, key_functor(gridSize1)), + cuS.dc_tvar.begin(), thrust::make_discard_iterator(), cuS.dw_tvar.begin() + ); + + + thrust::transform(thrust::cuda::par.on(cs->stream), + cuS.dw_tvar.begin(), cuS.dw_tvar.end(), + bm[T], bm[T], thrust::plus()); + + thrust::advance(bm[T], osize); + } + } } use1 += N; // Increment particle count @@ -766,6 +937,18 @@ void Cube::determine_coefficients_cuda() // host_coefs = cuS.df_coef; + if (requestSubsample) { + // T loop + // + for (int T=0; T retV = cuS.T_covr[T]; + + for (int n=0; n(retV[n]); + } + } + } + // Create a wavenumber tuple from a flattened index // auto indices = [&](int indx) @@ -1162,7 +1345,7 @@ void Cube::multistep_update_cuda() if (byPlanes) { // Adjust cached storage, if necessary // - cuS.resize_coefs(N, imx, gridSize, stride); + cuS.resize_coefs(N, imx, gridSize, stride, sampT); // Compute the coefficient contribution for each order // @@ -1176,7 +1359,7 @@ void Cube::multistep_update_cuda() // coefKernelCubeX<<stream>>> (toKernel(cs->cuda_particles), toKernel(cs->indx1), - toKernel(cuS.dN_coef), jj, kk, stride, cur); + toKernel(cuS.dN_coef), toKernel(cuS.u_d), jj, kk, stride, cur); unsigned int gridSize1 = N/BLOCK_SIZE; if (N > gridSize1*BLOCK_SIZE) gridSize1++; @@ -1214,7 +1397,7 @@ void Cube::multistep_update_cuda() else { // Adjust cached storage, if necessary // - cuS.resize_coefs(N, osize, gridSize, stride); + cuS.resize_coefs(N, osize, gridSize, stride, sampT); // Shared memory size for the reduction // @@ -1228,7 +1411,7 @@ void Cube::multistep_update_cuda() // coefKernelCube<<stream>>> (toKernel(cs->cuda_particles), toKernel(cs->indx1), - toKernel(cuS.dN_coef), stride, cur); + toKernel(cuS.dN_coef), toKernel(cuS.u_d), stride, cur); unsigned int gridSize1 = N/BLOCK_SIZE; if (N > gridSize1*BLOCK_SIZE) gridSize1++; From e580b21a3b4dfe616f849cf6b13b752ff8785170 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 2 Dec 2025 22:42:09 -0500 Subject: [PATCH 22/38] Reuse intermediate reduction variables --- src/Cube.H | 6 +----- src/cudaCube.cu | 36 ++++++++++++++++++++---------------- 2 files changed, 21 insertions(+), 21 deletions(-) diff --git a/src/Cube.H b/src/Cube.H index 68cae68f9..e6c3c2f4c 100644 --- a/src/Cube.H +++ b/src/Cube.H @@ -135,13 +135,9 @@ private: thrust::device_vector> dw_coef; thrust::device_vector> df_coef; - thrust::device_vector> dN_tvar; - thrust::device_vector> dc_tvar; - thrust::device_vector> dw_tvar; - thrust::device_vector u_d; - std::vector>> T_covr; + std::vector>> T_samp; void resize_coefs(int N, int osize, int gridSize, int stride, int sampT); diff --git a/src/cudaCube.cu b/src/cudaCube.cu index f3a387f74..cc1fa6c50 100644 --- a/src/cudaCube.cu +++ b/src/cudaCube.cu @@ -421,8 +421,6 @@ forceKernelCube(dArray P, dArray I, CmplxT acc[3] = {0.0, 0.0, 0.0}, pot = 0.0; cuFP_t pos[3] = {p.pos[0], p.pos[1], p.pos[2]}; - cuFP_t mm = p.mass; - int ind[3]; #ifdef NORECURSION // Index loop @@ -544,13 +542,12 @@ void Cube::cudaStorage::resize_coefs(int N, int osize, int gridSize, int stride, dc_coef.resize(osize*gridSize); dw_coef.resize(osize); // This will stay fixed + u_d.resize(N); + + // Resize covariance storage, if sampT>0 if (sampT) { - T_covr.resize(sampT); - for (int T=0; Tstream), cuS.df_coef.begin(), cuS.df_coef.end(), 0.0); + + if (requestSubsample) { + for (int T=0; Tstream), + cuS.T_samp[T].begin(), cuS.T_samp[T].end(), 0.0); + } + } } void Cube::determine_coefficients_cuda() @@ -685,7 +689,7 @@ void Cube::determine_coefficients_cuda() std::vector::iterator> bm; if (requestSubsample) { for (int T=0; T <<stream>>> - (toKernel(cuS.dc_tvar), toKernel(cuS.dN_tvar), imx, N, k, k+s); + (toKernel(cuS.dc_coef), toKernel(cuS.dN_coef), imx, N, k, k+s); // Finish the reduction for this order in parallel // @@ -797,12 +801,12 @@ void Cube::determine_coefficients_cuda() thrust::cuda::par.on(cs->stream), thrust::make_transform_iterator(index_begin, key_functor(gridSize1)), thrust::make_transform_iterator(index_end, key_functor(gridSize1)), - cuS.dc_tvar.begin(), thrust::make_discard_iterator(), cuS.dw_tvar.begin() + cuS.dc_coef.begin(), thrust::make_discard_iterator(), cuS.dw_coef.begin() ); thrust::transform(thrust::cuda::par.on(cs->stream), - cuS.dw_tvar.begin(), cuS.dw_tvar.end(), + cuS.dw_coef.begin(), cuS.dw_coef.end(), bm[T], bm[T], thrust::plus()); thrust::advance(bm[T], imx); @@ -905,7 +909,7 @@ void Cube::determine_coefficients_cuda() // reduceSumS <<stream>>> - (toKernel(cuS.dc_tvar), toKernel(cuS.dN_tvar), osize, N, k, k+s); + (toKernel(cuS.dc_coef), toKernel(cuS.dN_coef), osize, N, k, k+s); // Finish the reduction for this order in parallel // @@ -917,12 +921,12 @@ void Cube::determine_coefficients_cuda() thrust::cuda::par.on(cs->stream), thrust::make_transform_iterator(index_begin, key_functor(gridSize1)), thrust::make_transform_iterator(index_end, key_functor(gridSize1)), - cuS.dc_tvar.begin(), thrust::make_discard_iterator(), cuS.dw_tvar.begin() + cuS.dc_coef.begin(), thrust::make_discard_iterator(), cuS.dw_coef.begin() ); thrust::transform(thrust::cuda::par.on(cs->stream), - cuS.dw_tvar.begin(), cuS.dw_tvar.end(), + cuS.dw_coef.begin(), cuS.dw_coef.end(), bm[T], bm[T], thrust::plus()); thrust::advance(bm[T], osize); @@ -941,7 +945,7 @@ void Cube::determine_coefficients_cuda() // T loop // for (int T=0; T retV = cuS.T_covr[T]; + thrust::host_vector retV = cuS.T_samp[T]; for (int n=0; n(retV[n]); From cb1948624bc6457f4d294da5933456b95e6b2318 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 3 Dec 2025 13:32:55 -0500 Subject: [PATCH 23/38] Add base-class overrides --- include/plummer.H | 4 ++++ src/Cylinder.H | 4 ++++ src/Direct.H | 3 +++ src/ExternalForce.H | 3 +++ src/NoForce.H | 3 +++ src/Shells.H | 3 +++ src/SlabSL.H | 25 ++++++++++++++----------- src/Sphere.H | 3 +++ 8 files changed, 37 insertions(+), 11 deletions(-) diff --git a/include/plummer.H b/include/plummer.H index e1520ef92..4711b2500 100644 --- a/include/plummer.H +++ b/include/plummer.H @@ -6,6 +6,10 @@ class PlummerSphere : public AxiSymModel { + using AxiSymModel::get_mass; + using AxiSymModel::get_density; + using AxiSymModel::get_pot; + private: double rmin; diff --git a/src/Cylinder.H b/src/Cylinder.H index 0b39901e6..5143af617 100644 --- a/src/Cylinder.H +++ b/src/Cylinder.H @@ -125,6 +125,10 @@ class MixtureBasis; */ class Cylinder : public Basis { + + //! Pull in base-class overrides + using PotAccel::determine_coefficients; + private: bool precond, EVEN_M, subsamp, sech2 = false; diff --git a/src/Direct.H b/src/Direct.H index c4d705147..ec15423f1 100644 --- a/src/Direct.H +++ b/src/Direct.H @@ -48,6 +48,9 @@ class Direct : public PotAccel { + //! Pull in base-class overrides + using PotAccel::determine_coefficients; + private: int soft_indx; diff --git a/src/ExternalForce.H b/src/ExternalForce.H index c76c3d4a8..2c8744853 100644 --- a/src/ExternalForce.H +++ b/src/ExternalForce.H @@ -11,6 +11,9 @@ */ class ExternalForce : public PotAccel { + //! Bring in base-class overloads + using PotAccel::determine_coefficients; + private: void determine_coefficients(void); diff --git a/src/NoForce.H b/src/NoForce.H index 97f4f6682..5226e1003 100644 --- a/src/NoForce.H +++ b/src/NoForce.H @@ -14,6 +14,9 @@ */ class NoForce : public PotAccel { + //! Pull in base-class overrides + using PotAccel::determine_coefficients; + private: void initialize() {} diff --git a/src/Shells.H b/src/Shells.H index af15e0ad7..f97cab9bb 100644 --- a/src/Shells.H +++ b/src/Shells.H @@ -27,6 +27,9 @@ class Shells : public PotAccel { + //! Pull in base-class overrides + using PotAccel::determine_coefficients; + private: bool self_consistent, firstime_accel; diff --git a/src/SlabSL.H b/src/SlabSL.H index 573f04859..693ded977 100644 --- a/src/SlabSL.H +++ b/src/SlabSL.H @@ -23,17 +23,20 @@ class SlabSL : public PotAccel { -//! Header structure for Sturm-Liouville slab expansion -//! Used for deprecated native coefficient files -struct SlabSLCoefHeader { - double time; - double zmax; - double h; - int type; - int nmaxx, nmaxy, nmaxz; - int jmax; -}; - + //! Pull in base-class overrides + using PotAccel::determine_coefficients; + + //! Header structure for Sturm-Liouville slab expansion + //! Used for deprecated native coefficient files + struct SlabSLCoefHeader { + double time; + double zmax; + double h; + int type; + int nmaxx, nmaxy, nmaxz; + int jmax; + }; + private: std::shared_ptr grid; diff --git a/src/Sphere.H b/src/Sphere.H index b53aaa575..4c369bc78 100644 --- a/src/Sphere.H +++ b/src/Sphere.H @@ -69,6 +69,9 @@ typedef std::shared_ptr SLGridSphPtr; class Sphere : public SphericalBasis { + //! Pull in base-class overrides + using PotAccel::determine_coefficients; + private: SLGridSphPtr ortho; From 095aed8002be6336675ac5c0042c288a264d3795 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 3 Dec 2025 13:34:54 -0500 Subject: [PATCH 24/38] Fix memory allocation for subsampling; allow user to set Cuda reduction method; swap z,y,x packing for the correct x,y,z packing --- src/cudaCube.cu | 124 ++++++++++++++++++++++++------------------------ 1 file changed, 63 insertions(+), 61 deletions(-) diff --git a/src/cudaCube.cu b/src/cudaCube.cu index cc1fa6c50..949017c57 100644 --- a/src/cudaCube.cu +++ b/src/cudaCube.cu @@ -90,10 +90,6 @@ void testConstantsCube() // void Cube::cuda_initialize() { - // Default - // - byPlanes = true; - // Turn off full covariance // fullCovar = false; @@ -118,7 +114,8 @@ void Cube::cuda_initialize() if (cuMethod.find("1d") != std::string::npos) byPlanes = true; if (myid==0) - std::cout << "---- Cube::cuda_initialize: byPlanes=" + std::cout << "---- Cube::cuda_initialize: cuMethod=" << cuMethod << std::endl + << "---- Cube::cuda_initialize: byPlanes=" << std::boolalpha << byPlanes << std::endl; } @@ -275,7 +272,7 @@ __global__ void coefKernelCube __global__ void coefKernelCubeX (dArray P, dArray I, dArray coef, - dArray used, int jj, int kk, int stride, PII lohi) + dArray used, int ii, int jj, int stride, PII lohi) { // Thread ID // @@ -326,10 +323,10 @@ __global__ void coefKernelCubeX #ifdef NORECURSION // Index loop // - for (int s=0; sstream), cuS.T_samp[T].begin(), cuS.T_samp[T].end(), 0.0); } @@ -697,19 +700,19 @@ void Cube::determine_coefficients_cuda() // Adjust cached storage, if necessary // - cuS.resize_coefs(N, imx, gridSize, stride, sampT); + cuS.resize_coefs(N, osize, imz, gridSize, stride, sampT); // Compute the coefficient contribution for each order // auto beg = cuS.df_coef.begin(); - for (int kk=-nmaxz; kk<=nmaxz; kk++) { + for (int ii=-nmaxx; ii<=nmaxx; ii++) { for (int jj=-nmaxy; jj<=nmaxy; jj++) { coefKernelCubeX<<stream>>> (toKernel(cs->cuda_particles), toKernel(cs->indx1), - toKernel(cuS.dN_coef), toKernel(cuS.u_d), jj, kk, stride, cur); + toKernel(cuS.dN_coef), toKernel(cuS.u_d), ii, jj, stride, cur); // Begin the reduction by blocks [perhaps this should use a // stride?] @@ -719,12 +722,12 @@ void Cube::determine_coefficients_cuda() reduceSum <<stream>>> - (toKernel(cuS.dc_coef), toKernel(cuS.dN_coef), imx, N); + (toKernel(cuS.dc_coef), toKernel(cuS.dN_coef), imz, N); // Finish the reduction for this order in parallel // thrust::counting_iterator index_begin(0); - thrust::counting_iterator index_end(gridSize1*imx); + thrust::counting_iterator index_end(gridSize1*imz); // The key_functor indexes the sum reduced series by array index // @@ -740,8 +743,8 @@ void Cube::determine_coefficients_cuda() cuS.dw_coef.begin(), cuS.dw_coef.end(), beg, beg, thrust::plus()); - thrust::advance(beg, imx); - + thrust::advance(beg, imz); + if (requestSubsample) { int sN = N/sampT; @@ -757,16 +760,19 @@ void Cube::determine_coefficients_cuda() int s = sN; if (T==sampT-1) s = N - k; - // Mass accumulation - // - auto mbeg = cuS.u_d.begin(); - auto mend = mbeg; - thrust::advance(mbeg, sN*T); - if (T(mend - mbeg); + massV1 [0][T] += static_cast(thrust::reduce(mbeg, mend)); + } // Begin the reduction per grid block // @@ -789,12 +795,12 @@ void Cube::determine_coefficients_cuda() // reduceSumS <<stream>>> - (toKernel(cuS.dc_coef), toKernel(cuS.dN_coef), imx, N, k, k+s); + (toKernel(cuS.dc_coef), toKernel(cuS.dN_coef), imz, N, k, k+s); // Finish the reduction for this order in parallel // thrust::counting_iterator index_begin(0); - thrust::counting_iterator index_end(gridSize1*imx); + thrust::counting_iterator index_end(gridSize1*imz); thrust::reduce_by_key ( @@ -809,7 +815,7 @@ void Cube::determine_coefficients_cuda() cuS.dw_coef.begin(), cuS.dw_coef.end(), bm[T], bm[T], thrust::plus()); - thrust::advance(bm[T], imx); + thrust::advance(bm[T], imz); } } } @@ -818,10 +824,10 @@ void Cube::determine_coefficients_cuda() // END: z wave number loop } else { - + // Adjust cached storage, if necessary // - cuS.resize_coefs(N, osize, gridSize, stride, sampT); + cuS.resize_coefs(N, osize, osize, gridSize, stride, sampT); // Compute the coefficient contribution for each order // @@ -860,8 +866,6 @@ void Cube::determine_coefficients_cuda() cuS.dw_coef.begin(), cuS.dw_coef.end(), beg, beg, thrust::plus()); - thrust::advance(beg, osize); - if (requestSubsample) { int sN = N/sampT; @@ -885,8 +889,8 @@ void Cube::determine_coefficients_cuda() if (T(mend - mbeg); + massV1 [0][T] += static_cast(thrust::reduce(mbeg, mend)); // Begin the reduction per grid block // @@ -928,8 +932,6 @@ void Cube::determine_coefficients_cuda() thrust::transform(thrust::cuda::par.on(cs->stream), cuS.dw_coef.begin(), cuS.dw_coef.end(), bm[T], bm[T], thrust::plus()); - - thrust::advance(bm[T], osize); } } } @@ -1349,13 +1351,13 @@ void Cube::multistep_update_cuda() if (byPlanes) { // Adjust cached storage, if necessary // - cuS.resize_coefs(N, imx, gridSize, stride, sampT); + cuS.resize_coefs(N, osize, imz, gridSize, stride, sampT); // Compute the coefficient contribution for each order // auto beg = cuS.df_coef.begin(); - for (int kk=-nmaxz; kk<=nmaxz; kk++) { + for (int ii=-nmaxx; ii<=nmaxx; ii++) { for (int jj=-nmaxy; jj<=nmaxy; jj++) { @@ -1363,19 +1365,19 @@ void Cube::multistep_update_cuda() // coefKernelCubeX<<stream>>> (toKernel(cs->cuda_particles), toKernel(cs->indx1), - toKernel(cuS.dN_coef), toKernel(cuS.u_d), jj, kk, stride, cur); + toKernel(cuS.dN_coef), toKernel(cuS.u_d), ii, jj, stride, cur); unsigned int gridSize1 = N/BLOCK_SIZE; if (N > gridSize1*BLOCK_SIZE) gridSize1++; reduceSum <<stream>>> - (toKernel(cuS.dc_coef), toKernel(cuS.dN_coef), imx, N); + (toKernel(cuS.dc_coef), toKernel(cuS.dN_coef), imz, N); // Finish the reduction for this order in parallel // thrust::counting_iterator index_begin(0); - thrust::counting_iterator index_end(gridSize1*imx); + thrust::counting_iterator index_end(gridSize1*imz); // The key_functor indexes the sum reduced series by array index // @@ -1391,17 +1393,17 @@ void Cube::multistep_update_cuda() cuS.dw_coef.begin(), cuS.dw_coef.end(), beg, beg, thrust::plus()); - thrust::advance(beg, imx); + thrust::advance(beg, imz); } - // END: z wave numbers + // END: y wave numbers } - // END: y wave numbers + // END: x wave numbers } // END: bunch by planes else { // Adjust cached storage, if necessary // - cuS.resize_coefs(N, osize, gridSize, stride, sampT); + cuS.resize_coefs(N, osize, osize, gridSize, stride, sampT); // Shared memory size for the reduction // From 8840d329722acdfa2020c9ce7add39d5f1b3786e Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 3 Dec 2025 13:35:54 -0500 Subject: [PATCH 25/38] Add signature for updated coefficient storage; add base-class overrides --- src/Cube.H | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/Cube.H b/src/Cube.H index e6c3c2f4c..743153677 100644 --- a/src/Cube.H +++ b/src/Cube.H @@ -139,7 +139,7 @@ private: std::vector>> T_samp; - void resize_coefs(int N, int osize, int gridSize, + void resize_coefs(int N, int osize, int psize, int gridSize, int stride, int sampT); }; @@ -297,6 +297,9 @@ public: //! Destructor virtual ~Cube(); + //! Bring in base-class overloads + using PotAccel::determine_coefficients; + //! Compute the coefficients void determine_coefficients(void); From 7b3e9f06a4be4a23f426e8194340dcb98cd8ddab Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 3 Dec 2025 13:36:59 -0500 Subject: [PATCH 26/38] Save a small bit of time to zero rather than resize previously sized covariance elements --- src/Cube.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Cube.cc b/src/Cube.cc index 4a82b8acb..994c23112 100644 --- a/src/Cube.cc +++ b/src/Cube.cc @@ -414,7 +414,7 @@ void Cube::determine_coefficients(void) if (tnow > last) { requestSubsample = true; last = tnow; - init_covariance(); + zero_covariance(); } } } else { From c09d8fe842bf3ec7368d4b8c21a0510544c14cad Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 3 Dec 2025 17:38:30 -0500 Subject: [PATCH 27/38] Make GPU indexing scheme consistent with CPU --- src/cudaCube.cu | 21 ++++----------------- 1 file changed, 4 insertions(+), 17 deletions(-) diff --git a/src/cudaCube.cu b/src/cudaCube.cu index 949017c57..9f3098066 100644 --- a/src/cudaCube.cu +++ b/src/cudaCube.cu @@ -38,29 +38,16 @@ int cubeIndex(int i, int j, int k) i += cubeNumX; j += cubeNumY; k += cubeNumZ; - return k*cubeNX*cubeNY + j*cubeNX + i; -} - -/* -// Index function for modulus coefficients -// -__device__ -thrust::tuple TensorIndices(int indx) -{ - int k = indx/(cubeNX*cubeNY); - int j = (indx - k*cubeNX*cubeNY)/cubeNX; - int i = indx - (j + k*cubeNY)*cubeNX; - return {i, j, k}; + return (i*cubeNY + j)*cubeNZ + k; } -*/ __device__ thrust::tuple cubeWaveNumbers(int indx) { - int k = indx/(cubeNX*cubeNY); - int j = (indx - k*cubeNX*cubeNY)/cubeNX; - int i = indx - (j + k*cubeNY)*cubeNX; + int i = indx/(cubeNY*cubeNZ); + int j = (indx - i*cubeNY*cubeNZ)/cubeNY; + int k = indx - (j + i*cubeNY)*cubeNX; return {i-cubeNumX, j-cubeNumY, k-cubeNumZ}; } From e3030f5c7d5f43db16541cd3320fedd0e7eea667 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 3 Dec 2025 17:53:09 -0500 Subject: [PATCH 28/38] Comment addition only [no CI] --- src/Cube.H | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/Cube.H b/src/Cube.H index 743153677..8117b21af 100644 --- a/src/Cube.H +++ b/src/Cube.H @@ -58,6 +58,13 @@ @param subsampleFloat boolean to use single precision in HDF5 files (default: false) + + @note The subsampling for the cpu implemtation is round-robin + while the subsampling for the gpu implementation is a partition. + This implies that the coefficient values will agree but the + subsamples contain different particles. Since the goal of + subsampling is statistical assessment, this should not be a + limitation. */ class Cube : public PotAccel { From 53c51ae6067a98b97cc3c96e7064f1be1b5e5a8c Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 4 Dec 2025 13:18:45 -0500 Subject: [PATCH 29/38] Refactoring fix --- expui/BiorthBasis.H | 65 ++------------------------------------------- 1 file changed, 2 insertions(+), 63 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 77dc61018..85fe5b87f 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -494,6 +494,7 @@ namespace BasisClasses covarData[T][lm] = std::make_tuple(meanV[T][lm], covrV[T][lm]); } } + std::cout << "Spherical::writeCoefCovariance: #, mass=" << sampleCounts.sum() << ", " << sampleMasses.sum() << std::endl; SubsampleCovariance::CovarData elements = std::make_tuple(sampleCounts, sampleMasses, covarData); covarStore->writeCoefCovariance(compname, runtag, elements, time); } else { @@ -1091,6 +1092,7 @@ namespace BasisClasses { if (covarStore) { auto covarData = sl->getCoefCovariance(); + std::tie(sampleCounts, sampleMasses) = sl->getCovarSamples(); SubsampleCovariance::CovarData elements = std::make_tuple(sampleCounts, sampleMasses, covarData); covarStore->writeCoefCovariance(compname, runtag, elements, time); } else { @@ -1462,69 +1464,6 @@ namespace BasisClasses unsigned index1D(int kx, int ky, int kz); }; - class CovarianceReader - { - public: - - using CoefCovarType = std::tuple; - - protected: - - std::vector times; - std::vector sampleCounts; - std::vector sampleMasses; - std::vector>> covarData; - - const double fixedPointPrecision = 1.0e08; - std::map timeMap; - unsigned timeIndex(double time) - { - int itime = static_cast(time * fixedPointPrecision + 0.5); - auto it = timeMap.find(itime); - if (it == timeMap.end()) - throw std::runtime_error("CovarianceReader: time not found"); - return it->second; - } - - std::string basisID; - - public: - - //! Constructor from HDF5 file - CovarianceReader(const std::string& filename, int stride=1); - - //! Get time list - std::vector Times() { return times; } - - //@{ - //! Return a vector of tuples of basis functions and the - //! covariance matrix for subsamples of particles - std::vector> - getCoefCovariance(unsigned index) { return covarData[index]; } - - std::vector> - getCoefCovariance(double time) { return covarData[timeIndex(time)]; } - //@} - - - //@{ - //! Get sample counts for the covariance computation - std::tuple - getCovarSamples(unsigned index) - { - return {sampleCounts[index], sampleMasses[index]}; - } - - std::tuple - getCovarSamples(double time) - { - unsigned index = timeIndex(time); - return {sampleCounts[index], sampleMasses[index]}; - } - //@} - - std::string basisIDname() { return basisID; } - }; //! Time-dependent potential-density model using BasisCoef = std::tuple, std::shared_ptr>; From df26050616dd36e5acc540da977f85628d660346 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 4 Dec 2025 13:19:06 -0500 Subject: [PATCH 30/38] Preliminary implementation of subsample support for Spherical and Cylinder --- src/AxisymmetricBasis.H | 7 ++++ src/AxisymmetricBasis.cc | 17 ++++++++-- src/Cylinder.H | 2 +- src/Cylinder.cc | 15 +++++++++ src/SphericalBasis.H | 9 +++++ src/SphericalBasis.cc | 71 ++++++++++++++++++++++++++++++++++++++-- 6 files changed, 116 insertions(+), 5 deletions(-) diff --git a/src/AxisymmetricBasis.H b/src/AxisymmetricBasis.H index fdf2cbb3d..c4b473453 100644 --- a/src/AxisymmetricBasis.H +++ b/src/AxisymmetricBasis.H @@ -43,6 +43,9 @@ @param tk_type is the smoothing type, one of: Hall, VarianceCut, CumulativeCut, VarianceWeighted @param subsamp true sets partition variance computation (default: false) + + @param nint is the step cadence for subsample variance computation + (default: 0, no subsampling) */ class AxisymmetricBasis : public Basis { @@ -77,6 +80,9 @@ protected: //! First step for PCA computation int npca0; + // Frequency of subsample computation + int nint = 0; + //! Hall smoothing exponent (default: 1.0) double hexp; @@ -125,6 +131,7 @@ protected: //@{ //! Mass and counts for subsample covariance + std::vector countT, countT1; std::vector massT, massT1; unsigned sampT, defSampT; //@} diff --git a/src/AxisymmetricBasis.cc b/src/AxisymmetricBasis.cc index 7617eca0f..18fb43000 100644 --- a/src/AxisymmetricBasis.cc +++ b/src/AxisymmetricBasis.cc @@ -17,6 +17,7 @@ AxisymmetricBasis::valid_keys = "dof", "npca", "npca0", + "nint", "pcavar", "pcaeof", "pcadiag", @@ -45,6 +46,7 @@ AxisymmetricBasis:: AxisymmetricBasis(Component* c0, const YAML::Node& conf) : dof = 3; npca = 500; npca0 = 0; + nint = 0; pcavar = false; pcaeof = false; pcadiag = false; @@ -83,6 +85,15 @@ AxisymmetricBasis:: AxisymmetricBasis(Component* c0, const YAML::Node& conf) : if (conf["tk_type"]) tk_type = setTK(conf["tk_type"].as()); if (conf["Mmax"] and not conf["Lmax"]) Lmax = Mmax; + + if (conf["nint"]) { // For OutSample support + nint = conf["nint"].as(); + if (nint>0) { + pcavar = true; + subsamp = true; + } + } + } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing parameters in AxisymmetricBasis: " @@ -1207,10 +1218,12 @@ void AxisymmetricBasis::parallel_gather_coef2(void) // Report mass used [with storage sanity checks] // - if (massT1.size() == massT.size() and massT.size()==sampT) + if (massT1.size() == massT.size() and massT.size()==sampT) { + MPI_Allreduce(&countT1[0], &countT[0], sampT, + MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD); MPI_Allreduce(&massT1[0], &massT[0], sampT, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - else { + } else { std::cout << "[" << myid << "] AxisymmetricBasis: " << "coef2 out of bounds in mass" << std::endl; } diff --git a/src/Cylinder.H b/src/Cylinder.H index 5143af617..442ded197 100644 --- a/src/Cylinder.H +++ b/src/Cylinder.H @@ -149,7 +149,7 @@ private: // Parameters double rcylmin, rcylmax, zmax, acyl; - int nmaxfid, lmaxfid, mmax, mlim; + int nmaxfid, lmaxfid, mmax, mlim, nint; int ncylnx, ncylny, ncylr; double hcyl, hexp, snr, rem; int nmax, ncylodd, ncylrecomp, npca, npca0, nvtk, cmapR, cmapZ; diff --git a/src/Cylinder.cc b/src/Cylinder.cc index b33e5132e..5958bd951 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -62,6 +62,7 @@ Cylinder::valid_keys = { "pcavtk", "pcadiag", "subsamp", + "nint", "try_cache", "density", "EVEN_M", @@ -137,6 +138,7 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : pcadiag = false; pcaeof = false; subsamp = false; + nint = 0; nvtk = 1; pcainit = true; coef_dump = true; @@ -451,6 +453,7 @@ void Cylinder::initialize() if (conf["pcavtk" ]) pcavtk = conf["pcavtk" ].as(); if (conf["pcadiag" ]) pcadiag = conf["pcadiag" ].as(); if (conf["subsamp" ]) subsamp = conf["subsamp" ].as(); + if (conf["nint" ]) nint = conf["nint" ].as(); if (conf["try_cache" ]) try_cache = conf["try_cache" ].as(); if (conf["EVEN_M" ]) EVEN_M = conf["EVEN_M" ].as(); if (conf["cmap" ]) cmapR = conf["cmap" ].as(); @@ -943,6 +946,18 @@ void Cylinder::determine_coefficients_particles(void) pcainit = false; } + if (nint) { + compute = (mstep == 0) && (!(this_step % nint)); + if (compute) { + ortho->setSampT(sampT); + ortho->set_covar(true); + requestSubsample = true; + } else { + ortho->set_covar(false); + subsampleComputed = false; + } + } + if (pcavar or pcaeof) { if (this_step >= npca0) compute = (mstep == 0) && !( (this_step-npca0) % npca); diff --git a/src/SphericalBasis.H b/src/SphericalBasis.H index c4dc2730f..faaa7cb6a 100644 --- a/src/SphericalBasis.H +++ b/src/SphericalBasis.H @@ -544,6 +544,15 @@ public: //! Dump current coefficients into named HDF5 file void dump_coefs_h5(const std::string& file); + + /** Return a subsample of the basis coefficients and covariance + elements for matrices for analysis. The default implementation + returns empty vectors. This parallels the implementation in + expui. */ + CovarData getSubsample(); + + //! Indicate whether or not subsampling is supported + virtual bool hasSubsample() { return true; } }; #endif diff --git a/src/SphericalBasis.cc b/src/SphericalBasis.cc index be24ecc03..cba710164 100644 --- a/src/SphericalBasis.cc +++ b/src/SphericalBasis.cc @@ -492,7 +492,8 @@ void * SphericalBasis::determine_coefficients_thread(void * arg) if (pcavar) { whch = indx % sampT; pthread_mutex_lock(&cc_lock); - massT1[whch] += mass; + countT1[whch] += 1; + massT1 [whch] += mass; pthread_mutex_unlock(&cc_lock); } } @@ -683,8 +684,12 @@ void SphericalBasis::determine_coefficients_particles(void) if (pcavar or pcaeof) { if (this_step >= npca0) compute = (mstep == 0) && !( (this_step-npca0) % npca); - else + else if (nint) + compute = (mstep == 0) && !(this_step % nint); + else { compute = false; + subsampleComputed = false; + } } @@ -692,11 +697,15 @@ void SphericalBasis::determine_coefficients_particles(void) if (compute) { + requestSubsample = true; + if (massT.size() == 0) { // Allocate storage for subsampling if (defSampT) sampT = defSampT; else sampT = floor(sqrt(component->CurTotal())); massT .resize(sampT, 0); massT1 .resize(sampT, 0); + countT .resize(sampT, 0); + countT1 .resize(sampT, 0); expcoefT .resize(sampT); for (auto & t : expcoefT ) { @@ -735,6 +744,7 @@ void SphericalBasis::determine_coefficients_particles(void) if (pcavar) { for (auto & t : expcoefT1) { for (auto & v : t) v->setZero(); } for (auto & t : expcoefM1) { for (auto & v : t) v->setZero(); } + for (auto & v : countT1) v = 0; for (auto & v : massT1) v = 0; } @@ -893,6 +903,9 @@ void SphericalBasis::determine_coefficients_particles(void) } pca_hall(compute); + + requestSubsample = false; + subsampleComputed = true; } print_timings("SphericalBasis: coefficient timings"); @@ -2349,3 +2362,57 @@ void SphericalBasis::biorthogonality_check() MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); } } + +PotAccel::CovarData SphericalBasis::getSubsample() +{ + using covTuple = std::tuple; + // Prepare the covariance structure + std::vector< std::vector > covar(sampT); + + // Resize the inner vector to 1 (one component for Cube) + int totL = (Lmax+1)*(Lmax+2)/2; + for (auto & v : covar) v.resize(totL); + + // Fill the covariance structure with subsamples + for (int T=0; T((*expcoefT[T][L ])(n)); + else + meanV(n) = std::complex((*expcoefT[T][L ])(n), + (*expcoefT[T][L+1])(n)); + if (fullCovar) { + for (int nn=0; nn((*expcoefM[T][L ])(n, nn)); + else + covrV(n, nn) = std::complex((*expcoefM[T][L ])(n, nn), + (*expcoefM[T][L+1])(n, nn)); + } + } + } + + covar[T][L] = std::make_tuple(meanV, covrV); + + if (m==0) L++; + else L+=2; + } + // END: m loop + } + // END: l loop + + sampleCounts[T] = countT[T]; + sampleMasses[T] = massT [T]; + } + // END: T loop + + return {sampleCounts, sampleMasses, covar}; +} From 8521af78bdcb2ecba47462caf7478979c5fc867d Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 4 Dec 2025 13:22:11 -0500 Subject: [PATCH 31/38] Remove debug output line --- expui/BiorthBasis.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 85fe5b87f..26def95a6 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -494,7 +494,7 @@ namespace BasisClasses covarData[T][lm] = std::make_tuple(meanV[T][lm], covrV[T][lm]); } } - std::cout << "Spherical::writeCoefCovariance: #, mass=" << sampleCounts.sum() << ", " << sampleMasses.sum() << std::endl; + SubsampleCovariance::CovarData elements = std::make_tuple(sampleCounts, sampleMasses, covarData); covarStore->writeCoefCovariance(compname, runtag, elements, time); } else { From a662ac975a5a98d069e1f3664ded59860fa44ada Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 4 Dec 2025 17:11:41 -0500 Subject: [PATCH 32/38] Fix debug routine for EL3's use of Eigen rather than Vector --- src/Cylinder.H | 8 ++++++++ src/Cylinder.cc | 16 ++++++++++++++-- src/Orient.cc | 35 +++++++++++++++++------------------ 3 files changed, 39 insertions(+), 20 deletions(-) diff --git a/src/Cylinder.H b/src/Cylinder.H index 442ded197..27a83a562 100644 --- a/src/Cylinder.H +++ b/src/Cylinder.H @@ -484,6 +484,14 @@ public: //! Sanity check on grid: dumps SM-style images of initial field void dump_mzero(const string& name, int step); + /** Return a subsample of the basis coefficients and covariance + elements for matrices for analysis. The default implementation + returns empty vectors. This parallels the implementation in + expui. */ + CovarData getSubsample(); + + //! Indicate whether or not subsampling is supported + virtual bool hasSubsample() { return true; } }; diff --git a/src/Cylinder.cc b/src/Cylinder.cc index 5958bd951..29c323def 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -185,7 +185,11 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : if (mlim>=0) ortho->set_mlim(mlim); if (EVEN_M) ortho->setEven(EVEN_M); ortho->setSampT(defSampT); - + if (nint>0) { + EmpCylSL::PCAVAR = true; + ortho->init_pca(); + } + try { if (conf["tk_type"]) ortho->setTK(conf["tk_type"].as()); } @@ -947,9 +951,9 @@ void Cylinder::determine_coefficients_particles(void) } if (nint) { + // Determine if we need to compute covariance this step compute = (mstep == 0) && (!(this_step % nint)); if (compute) { - ortho->setSampT(sampT); ortho->set_covar(true); requestSubsample = true; } else { @@ -1831,3 +1835,11 @@ void Cylinder::occt_output() } } } + + +PotAccel::CovarData Cylinder::getSubsample() +{ + std::tie(sampleCounts, sampleMasses) = ortho->getCovarSamples(); + auto covarData = ortho->getCoefCovariance(); + return {sampleCounts, sampleMasses, covarData}; +} diff --git a/src/Orient.cc b/src/Orient.cc index a2f4b48a0..9f0eb37c6 100644 --- a/src/Orient.cc +++ b/src/Orient.cc @@ -18,21 +18,20 @@ void EL3::debug() const { if (myid==0) { - cerr << left << setfill('-'); - ostringstream ostr; - ostr << "--- EL3 [" << myid << "] "; - cerr << setw(60) << ostr.str() << endl << setfill(' '); - cerr.precision(4); - cerr << setw(12) << T - << setw(12) << M - << setw(12) << E; - cerr << " ["; - for (int i=1; i<=3; i++) cerr << setw(12) << L[i]; - cerr << "] ["; - for (int i=1; i<=3; i++) cerr << setw(12) << R[i]; - cerr << "]" << endl; - cerr << left << setfill('-') - << setw(60) << '-' << endl << setfill(' '); + std::cerr << left << setfill('-'); + std::ostringstream ostr; ostr << "--- EL3 [" << myid << "] "; + std::cerr << setw(60) << ostr.str() << std::endl << std::setfill(' '); + std::cerr.precision(4); + std::cerr << std::setw(12) << T + << std::setw(12) << M + << std::setw(12) << E; + std::cerr << " ["; + for (int i=0; i<3; i++) std::cerr << std::setw(12) << L[i]; + std::cerr << "] ["; + for (int i=0; i<3; i++) std::cerr << std::setw(12) << R[i]; + std::cerr << "]" << std::endl; + std::cerr << std::left << std::setfill('-') + << std::setw(60) << '-' << std::endl << std::setfill(' '); } } @@ -50,9 +49,9 @@ Orient::Orient(int n, int nwant, int naccel, unsigned Oflg, unsigned Cflg, damp = damping; linear = false; - pos = vector(3); - psa = vector(3); - vel = vector(3); + pos = std::vector(3); + psa = std::vector(3); + vel = std::vector(3); center .setZero(); center0.setZero(); From a0f3f584e154b12f96c575aeef0e6db45a54725e Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 4 Dec 2025 18:17:30 -0500 Subject: [PATCH 33/38] Add missing methods needed for Covariance --- src/Cube.H | 1 - src/Cylinder.H | 3 +++ src/Cylinder.cc | 11 ++++++++ src/SphericalBasis.H | 3 +++ src/SphericalBasis.cc | 61 ++++++++++++++++++++++++++++--------------- 5 files changed, 57 insertions(+), 22 deletions(-) diff --git a/src/Cube.H b/src/Cube.H index 8117b21af..80b3ad5a8 100644 --- a/src/Cube.H +++ b/src/Cube.H @@ -292,7 +292,6 @@ private: to = tmp; } - //! Write basis-specific parameters to HDF5 covariance file void writeCovarH5Params(HighFive::File& file); diff --git a/src/Cylinder.H b/src/Cylinder.H index 27a83a562..3797d6d9a 100644 --- a/src/Cylinder.H +++ b/src/Cylinder.H @@ -324,6 +324,9 @@ protected: //! Valid keys for YAML configurations static const std::set valid_keys; + //! Write basis-specific parameters to HDF5 covariance file + void writeCovarH5Params(HighFive::File& file); + public: //! Mutexes for multithreading diff --git a/src/Cylinder.cc b/src/Cylinder.cc index 29c323def..331682700 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -1843,3 +1843,14 @@ PotAccel::CovarData Cylinder::getSubsample() auto covarData = ortho->getCoefCovariance(); return {sampleCounts, sampleMasses, covarData}; } + +void Cylinder::writeCovarH5Params(HighFive::File& file) +{ + file.createAttribute("mmax", HighFive::DataSpace::From(mmax)).write(mmax); + file.createAttribute("nmax", HighFive::DataSpace::From(nmax)).write(nmax); + file.createAttribute("rcylmin", HighFive::DataSpace::From(rcylmin)).write(rcylmin); + file.createAttribute("rcylmax", HighFive::DataSpace::From(rcylmax)).write(rcylmax); + file.createAttribute("acyl", HighFive::DataSpace::From(acyl)).write(acyl); + file.createAttribute("hcyl", HighFive::DataSpace::From(hcyl)).write(hcyl); +} + diff --git a/src/SphericalBasis.H b/src/SphericalBasis.H index faaa7cb6a..24bcdb63f 100644 --- a/src/SphericalBasis.H +++ b/src/SphericalBasis.H @@ -407,6 +407,9 @@ protected: void biorthogonality_check(); //@} + //! Write basis-specific parameters to HDF5 covariance file + void writeCovarH5Params(HighFive::File& file); + public: //! Use YAML header in coefficient file diff --git a/src/SphericalBasis.cc b/src/SphericalBasis.cc index cba710164..36339f949 100644 --- a/src/SphericalBasis.cc +++ b/src/SphericalBasis.cc @@ -2370,49 +2370,68 @@ PotAccel::CovarData SphericalBasis::getSubsample() // Prepare the covariance structure std::vector< std::vector > covar(sampT); - // Resize the inner vector to 1 (one component for Cube) + // Number of real angular terms l in [0, Lmax], m in [0, l] int totL = (Lmax+1)*(Lmax+2)/2; + + // Sanity check + if (expcoefT[0].size() != totL) { + std::ostringstream ss; + ss << "SphericalBasis::getSubsample() " + << "internal error: unexpected number of absolute angular terms, " + << expcoefT[0].size() << " != " << totL; + throw std::runtime_error(ss.str()); + } + + // Resize the covariance structure for (auto & v : covar) v.resize(totL); + // Resize the sample counts and masses + sampleCounts.resize(sampT); + sampleMasses.resize(sampT); + // Fill the covariance structure with subsamples for (int T=0; T((*expcoefT[T][L ])(n)); - else - meanV(n) = std::complex((*expcoefT[T][L ])(n), - (*expcoefT[T][L+1])(n)); + meanV(n) = std::complex((*expcoefT[T][k])(n)); + if (fullCovar) { + // inner n loop for (int nn=0; nn((*expcoefM[T][L ])(n, nn)); - else - covrV(n, nn) = std::complex((*expcoefM[T][L ])(n, nn), - (*expcoefM[T][L+1])(n, nn)); + covrV(n, nn) = std::complex((*expcoefM[T][k])(n, nn)); } } } - - covar[T][L] = std::make_tuple(meanV, covrV); - - if (m==0) L++; - else L+=2; + // Assign data to return structure + covar[T][k++] = std::make_tuple(meanV, covrV); } // END: m loop } // END: l loop sampleCounts[T] = countT[T]; - sampleMasses[T] = massT [T]; + sampleMasses[T] = massT[T]; } // END: T loop return {sampleCounts, sampleMasses, covar}; } + + +void SphericalBasis::writeCovarH5Params(HighFive::File& file) +{ + file.createAttribute("lmax", HighFive::DataSpace::From(Lmax)).write(Lmax); + file.createAttribute("nmax", HighFive::DataSpace::From(nmax)).write(nmax); + file.createAttribute("scale", HighFive::DataSpace::From(scale)).write(scale); + file.createAttribute("rmin", HighFive::DataSpace::From(rmin)).write(rmin); + file.createAttribute("rmax", HighFive::DataSpace::From(rmax)).write(rmax); +} + From 92b1e8f94dcc93dadbde848c5ffb47a377887307 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 4 Dec 2025 23:02:57 -0500 Subject: [PATCH 34/38] Add a constructor from an explicit path+filename --- expui/Covariance.cc | 18 ++++++++++++++---- include/Covariance.H | 15 ++++++++++++--- 2 files changed, 26 insertions(+), 7 deletions(-) diff --git a/expui/Covariance.cc b/expui/Covariance.cc index c10449081..b7910fe24 100644 --- a/expui/Covariance.cc +++ b/expui/Covariance.cc @@ -203,6 +203,20 @@ namespace BasisClasses void SubsampleCovariance::writeCoefCovariance (const std::string& compname, const std::string& runtag, CovarData& elem, double time) + { + // Only root process writes + // + if (myid) return; + + // The H5 filename + // + std::string fname = "coefcovar." + compname + "." + runtag + ".h5"; + + writeCoefCovariance(fname, elem, time); + } + + void SubsampleCovariance::writeCoefCovariance + (const std::string& fname, CovarData& elem, double time) { // Only root process writes // @@ -222,10 +236,6 @@ namespace BasisClasses // time = roundTime(time); - // The H5 filename - // - std::string fname = "coefcovar." + compname + "." + runtag + ".h5"; - // Check if file exists? // try { diff --git a/include/Covariance.H b/include/Covariance.H index d9636b66f..d894288ae 100644 --- a/include/Covariance.H +++ b/include/Covariance.H @@ -120,7 +120,8 @@ namespace BasisClasses { auto it = covarData.find(roundTime(time)); if (it == covarData.end()) - throw std::runtime_error("SubsampleCovariance::getCoefCovariance: time not found"); + throw std::runtime_error("SubsampleCovariance::getCoefCovariance: " + "time not found"); return std::get<2>(it->second); } @@ -130,16 +131,24 @@ namespace BasisClasses { auto it = covarData.find(roundTime(time)); if (it == covarData.end()) - throw std::runtime_error("SubsampleCovariance::getCovarSamples: time not found"); + throw std::runtime_error("SubsampleCovariance::getCovarSamples: " + "time not found"); return std::make_tuple(std::get<0>(it->second), std::get<1>(it->second)); } - //! Write coefficient covariance data to an HDF5 file + //! Write coefficient covariance data to an HDF5 file, creating a + //! filename based on component name and run tag void writeCoefCovariance(const std::string& compname, const std::string& runtag, CovarData& covarElem, double time=0.0); + //! Write coefficient covariance data to an HDF5 file with a + //! supplied path + filename + void writeCoefCovariance(const std::string& filename, + CovarData& covarElem, + double time=0.0); + //! Choose between float and double storage for covariance void setCovarFloatType(bool flt) { floatType = flt; } From d913617b8e8773f0d7250f5ee23ce89ddbbe1906 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 4 Dec 2025 23:03:27 -0500 Subject: [PATCH 35/38] Use the explicit path+filename constructor for Covariance --- src/OutSample.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OutSample.cc b/src/OutSample.cc index b648b0220..b0b0ff857 100644 --- a/src/OutSample.cc +++ b/src/OutSample.cc @@ -117,6 +117,6 @@ void OutSample::Run(int n, int mstep, bool last) // Write the subsample data auto elem = tcomp->force->getSubsample(); - covarStore->writeCoefCovariance(tcomp->name, runtag, elem, tnow); + covarStore->writeCoefCovariance(filename, elem, tnow); } From f29c04d5864eb92a12e63cdf3f7fdcc5c92f949d Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 4 Dec 2025 23:03:47 -0500 Subject: [PATCH 36/38] Fix OutSample logic --- src/Cylinder.cc | 18 +++++++++++++++--- src/SphericalBasis.cc | 42 +++++++++++++++++++++++++++++++++--------- 2 files changed, 48 insertions(+), 12 deletions(-) diff --git a/src/Cylinder.cc b/src/Cylinder.cc index 331682700..4fe91d50b 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -950,10 +950,16 @@ void Cylinder::determine_coefficients_particles(void) pcainit = false; } + // No covarance by default + // + compute = false; + + // Subsample computation flag + // if (nint) { - // Determine if we need to compute covariance this step - compute = (mstep == 0) && (!(this_step % nint)); - if (compute) { + // Computed only for mstep==0 and every nint steps + if (mstep==0 and this_step % nint == 0) { + compute = true; ortho->set_covar(true); requestSubsample = true; } else { @@ -1060,6 +1066,12 @@ void Cylinder::determine_coefficients_particles(void) if ((pcavar or pcaeof) and mlevel==0) ortho->pca_hall(compute, subsamp); + // If subsample requested and computed, turn off for next time + if (nint and compute and mlevel==multistep) { + requestSubsample = false; + subsampleComputed = true; + } + //========================= // Apply Hall smoothing //========================= diff --git a/src/SphericalBasis.cc b/src/SphericalBasis.cc index 36339f949..5f2ce4319 100644 --- a/src/SphericalBasis.cc +++ b/src/SphericalBasis.cc @@ -681,17 +681,29 @@ void SphericalBasis::determine_coefficients_particles(void) // if (!self_consistent && !firstime_coef && !initializing) return; - if (pcavar or pcaeof) { - if (this_step >= npca0) - compute = (mstep == 0) && !( (this_step-npca0) % npca); - else if (nint) - compute = (mstep == 0) && !(this_step % nint); + // No covarance by default + // + compute = false; + + // Subsample computation flag + // + if (nint) { + // Computed only for mstep==0 and every nint steps + if (mstep==0 and this_step % nint == 0) { + compute = true; + requestSubsample = true; + } else { - compute = false; + requestSubsample = false; subsampleComputed = false; + } } + if (pcavar or pcaeof) { + if (this_step >= npca0) + compute = (mstep == 0) && !( (this_step-npca0) % npca); + } int loffset, moffset, use1; @@ -900,12 +912,12 @@ void SphericalBasis::determine_coefficients_particles(void) for (int i=0; i Date: Wed, 10 Dec 2025 16:12:53 -0500 Subject: [PATCH 37/38] Should be a project not system include here --- expui/Coefficients.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/Coefficients.H b/expui/Coefficients.H index 061f42ac9..9effd7c5b 100644 --- a/expui/Coefficients.H +++ b/expui/Coefficients.H @@ -18,7 +18,7 @@ #include "CoefStruct.H" // Unit validation -#include +#include "UnitValidator.H" // All coefficient classes are in this namespace namespace CoefClasses From d8f9c75185aa1482f77357631fadc69eac43484f Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 13 Dec 2025 12:33:00 -0500 Subject: [PATCH 38/38] Implemented coefficient parameter checking on reading and writing --- expui/Coefficients.H | 27 ++++ expui/Coefficients.cc | 354 ++++++++++++++++++++++++++++++++++++++++++ pyEXP/CoefWrappers.cc | 28 ++++ 3 files changed, 409 insertions(+) diff --git a/expui/Coefficients.H b/expui/Coefficients.H index 9effd7c5b..b6f21c6e1 100644 --- a/expui/Coefficients.H +++ b/expui/Coefficients.H @@ -97,6 +97,9 @@ namespace CoefClasses //! Write parameter attributes (needed for derived classes) virtual void WriteH5Params(HighFive::File& file) = 0; + //! Checking file for parameter consistency + virtual bool CheckH5Params(HighFive::File& file) = 0; + //! Write coefficient data in H5 virtual unsigned WriteH5Times(HighFive::Group& group, unsigned count) = 0; @@ -311,6 +314,9 @@ namespace CoefClasses //! Write parameter attributes virtual void WriteH5Params(HighFive::File& file); + //! Check parameter attributes + virtual bool CheckH5Params(HighFive::File& file); + //! Write coefficient data in H5 virtual unsigned WriteH5Times(HighFive::Group& group, unsigned count); @@ -447,6 +453,9 @@ namespace CoefClasses //! Write parameter attributes virtual void WriteH5Params(HighFive::File& file); + //! Check parameter attributes; returns true for success + virtual bool CheckH5Params(HighFive::File& file); + //! Write coefficient data in H5 virtual unsigned WriteH5Times(HighFive::Group& group, unsigned count); @@ -593,6 +602,9 @@ namespace CoefClasses //! Write parameter attributes virtual void WriteH5Params(HighFive::File& file); + //! Check parameter attributes + virtual bool CheckH5Params(HighFive::File& file); + //! Write coefficient data in H5 virtual unsigned WriteH5Times(HighFive::Group& group, unsigned count); @@ -724,6 +736,9 @@ namespace CoefClasses //! Write parameter attributes virtual void WriteH5Params(HighFive::File& file); + //! Check parameter attributes + virtual bool CheckH5Params(HighFive::File& file); + //! Write coefficient data in H5 virtual unsigned WriteH5Times(HighFive::Group& group, unsigned count); @@ -849,6 +864,9 @@ namespace CoefClasses //! Write parameter attributes virtual void WriteH5Params(HighFive::File& file); + //! Check parameter attributes + virtual bool CheckH5Params(HighFive::File& file); + //! Write coefficient data in H5 virtual unsigned WriteH5Times(HighFive::Group& group, unsigned count); @@ -958,6 +976,9 @@ namespace CoefClasses //! Write parameter attributes virtual void WriteH5Params(HighFive::File& file); + //! Check parameter attributes + virtual bool CheckH5Params(HighFive::File& file); + //! Write coefficient data in H5 virtual unsigned WriteH5Times(HighFive::Group& group, unsigned count); @@ -1072,6 +1093,9 @@ namespace CoefClasses //! Write parameter attributes virtual void WriteH5Params(HighFive::File& file); + //! Check parameter attributes + virtual bool CheckH5Params(HighFive::File& file); + //! Write coefficient data in H5 virtual unsigned WriteH5Times(HighFive::Group& group, unsigned count); @@ -1202,6 +1226,9 @@ namespace CoefClasses //! Write parameter attributes virtual void WriteH5Params(HighFive::File& file); + //! Check parameter attributes + virtual bool CheckH5Params(HighFive::File& file); + //! Write coefficient data in H5 virtual unsigned WriteH5Times(HighFive::Group& group, unsigned count); diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc index dc7444dca..45c7d57a1 100644 --- a/expui/Coefficients.cc +++ b/expui/Coefficients.cc @@ -283,6 +283,18 @@ namespace CoefClasses auto in = stanza.getDataSet("coefficients").read(); + // Sanity check on shape + // + if (in.rows() != (Lmax+1)*(Lmax+2)/2 or in.cols() != Nmax) { + if (myid==0) { + std::cout << "---- SphCoefs WARNING: skipping time=" << Time + << " with shape (" << in.rows() << ", " << in.cols() << ")" + << " expected (" << (Lmax+1)*(Lmax+2)/2 + << ", " << Nmax << ")" << std::endl; + } + continue; + } + // If we have a legacy set of coefficients, re-order the // coefficients to match the new HighFive/Eigen ordering // @@ -840,6 +852,50 @@ namespace CoefClasses file.createAttribute("forceID", HighFive::DataSpace::From(forceID)).write(forceID); } + bool SphCoefs::CheckH5Params(HighFive::File& file) + { + double scale0 = coefs.begin()->second->scale, scale1; + std::string forceID0(coefs.begin()->second->id), forceID1; + int Lmax1, Nmax1; + + file.getAttribute("lmax" ).read(Lmax1 ); + file.getAttribute("nmax" ).read(Nmax1 ); + file.getAttribute("scale" ).read(scale1 ); + file.getAttribute("forceID").read(forceID1); + + bool ret = true; + + if (Lmax != Lmax1) { + if (myid==0) + std::cout << "---- SphCoefs::CheckH5Params: Lmax mismatch " << Lmax + << " != " << Lmax1 << std::endl; + ret = false; + } + + if (Nmax != Nmax1) { + if (myid==0) + std::cout << "---- SphCoefs::CheckH5Params: Nmax mismatch " << Nmax + << " != " << Nmax1 << std::endl; + ret = false; + } + + if (std::abs(scale0 - scale1)/scale0 > 1.e-8) { + if (myid==0) + std::cout << "---- SphCoefs::CheckH5Params: scale mismatch " << scale0 + << " != " << scale1 << std::endl; + ret = false; + } + + if (forceID0 != forceID1) { + if (myid==0) + std::cout << "---- SphCoefs::CheckH5Params: forceID mismatch <" << forceID0 + << "> != <" << forceID1 << ">" << std::endl; + ret = false; + } + + return ret; + } + unsigned SphCoefs::WriteH5Times(HighFive::Group& snaps, unsigned count) { for (auto c : coefs) { @@ -1064,6 +1120,17 @@ namespace CoefClasses auto in = stanza.getDataSet("coefficients").read(); + // Sanity check on shape + // + if (in.rows() != Mmax+1 or in.cols() != Nmax) { + if (myid==0) + std::cout << "---- CylCoefs WARNING: skipping time=" << Time + << " with shape (" << in.rows() << ", " << in.cols() << ")" + << " expected (" << Mmax+1 << ", " << Nmax << ")" + << std::endl; + continue; + } + // If we have a legacy set of coefficients, re-order the // coefficients to match the new HighFive/Eigen ordering // @@ -1256,6 +1323,41 @@ namespace CoefClasses file.createAttribute("forceID", HighFive::DataSpace::From(forceID)).write(forceID); } + bool CylCoefs::CheckH5Params(HighFive::File& file) + { + std::string forceID0(coefs.begin()->second->id), forceID1; + int Mmax1, Nmax1; + + file.getAttribute("mmax" ).read(Mmax1 ); + file.getAttribute("nmax" ).read(Nmax1 ); + file.getAttribute("forceID").read(forceID1); + + bool ret = true; + + if (Mmax != Mmax1) { + if (myid==0) + std::cout << "---- CylCoefs::CheckH5Params: Mmax mismatch " << Mmax + << " != " << Mmax1 << std::endl; + ret = false; + } + + if (Nmax != Nmax1) { + if (myid==0) + std::cout << "---- CylCoefs::CheckH5Params: Nmax mismatch " << Nmax + << " != " << Nmax1 << std::endl; + ret = false; + } + + if (forceID0 != forceID1) { + if (myid==0) + std::cout << "---- CylCoefs::CheckH5Params: forceID mismatch <" << forceID0 + << "> != <" << forceID1 << ">" << std::endl; + ret = false; + } + + return ret; + } + unsigned CylCoefs::WriteH5Times(HighFive::Group& snaps, unsigned count) { for (auto c : coefs) { @@ -1591,6 +1693,48 @@ namespace CoefClasses file.createAttribute("forceID", HighFive::DataSpace::From(forceID)).write(forceID); } + bool SlabCoefs::CheckH5Params(HighFive::File& file) + { + std::string forceID0(coefs.begin()->second->id), forceID1; + int NmaxX1, NmaxY1, NmaxZ1; + + file.getAttribute("nmaxx" ).read(NmaxX1 ); + file.getAttribute("nmaxy" ).read(NmaxY1 ); + file.getAttribute("nmaxz" ).read(NmaxZ1 ); + file.getAttribute("forceID").read(forceID1); + + bool ret = true; + + if (NmaxX != NmaxX1) { + if (myid==0) + std::cout << "---- SlabCoefs::CheckH5Params: NmaxX mismatch " << NmaxX + << " != " << NmaxX1 << std::endl; + ret = false; + } + + if (NmaxY != NmaxY1) { + if (myid==0) + std::cout << "---- SlabCoefs::CheckH5Params: NmaxY mismatch " << NmaxY + << " != " << NmaxY1 << std::endl; + ret = false; + } + + if (NmaxZ != NmaxZ1) { + if (myid==0) + std::cout << "---- SlabCoefs::CheckH5Params: NmaxZ mismatch " << NmaxZ + << " != " << NmaxZ1 << std::endl; + ret = false; + } + + if (forceID0 != forceID1) { + if (myid==0) + std::cout << "---- SlabCoefs::CheckH5Params: forceID mismatch <" << forceID0 + << "> != <" << forceID1 << ">" << std::endl; + ret = false; + } + return ret; + } + unsigned SlabCoefs::WriteH5Times(HighFive::Group& snaps, unsigned count) { for (auto c : coefs) { @@ -1954,6 +2098,48 @@ namespace CoefClasses file.createAttribute("forceID", HighFive::DataSpace::From(forceID)).write(forceID); } + bool CubeCoefs::CheckH5Params(HighFive::File& file) + { + std::string forceID0(coefs.begin()->second->id), forceID1; + int NmaxX1, NmaxY1, NmaxZ1; + + file.getAttribute("nmaxx" ).read(NmaxX1 ); + file.getAttribute("nmaxy" ).read(NmaxY1 ); + file.getAttribute("nmaxz" ).read(NmaxZ1 ); + file.getAttribute("forceID").read(forceID1); + + bool ret = true; + + if (NmaxX != NmaxX1) { + if (myid==0) + std::cout << "---- CubeCoefs::CheckH5Params: NmaxX mismatch " << NmaxX + << " != " << NmaxX1 << std::endl; + ret = false; + } + + if (NmaxY != NmaxY1) { + if (myid==0) + std::cout << "---- CubeCoefs::CheckH5Params: NmaxY mismatch " << NmaxY + << " != " << NmaxY1 << std::endl; + ret = false; + } + + if (NmaxZ != NmaxZ1) { + if (myid==0) + std::cout << "---- CubeCoefs::CheckH5Params: NmaxZ mismatch " << NmaxZ + << " != " << NmaxZ1 << std::endl; + ret = false; + } + if (forceID0 != forceID1) { + if (myid==0) + std::cout << "---- CubeCoefs::CheckH5Params: forceID mismatch <" << forceID0 + << "> != <" << forceID1 << ">" << std::endl; + ret = false; + } + + return ret; + } + unsigned CubeCoefs::WriteH5Times(HighFive::Group& snaps, unsigned count) { for (auto c : coefs) { @@ -2308,6 +2494,35 @@ namespace CoefClasses file.createAttribute("rank", HighFive::DataSpace::From(rank)).write(rank); } + bool TrajectoryData::CheckH5Params(HighFive::File& file) + { + int traj1, rank1; + + file.getAttribute("traj").read(traj1); + file.getAttribute("rank").read(rank1); + + int traj0 = coefs.begin()->second->traj; + int rank0 = coefs.begin()->second->rank; + + bool ret = true; + + if (traj0 != traj1) { + if (myid==0) + std::cout << "---- TrajectoryData::CheckH5Params: traj mismatch " << traj0 + << " != " << traj1 << std::endl; + ret = false; + } + + if (rank0 != rank1) { + if (myid==0) + std::cout << "---- TrajectoryData::CheckH5Params: rank mismatch " << rank0 + << " != " << rank1 << std::endl; + ret = false; + } + + return ret; + } + unsigned TrajectoryData::WriteH5Times(HighFive::Group& snaps, unsigned count) { snaps.createDataSet("times", times); @@ -2552,6 +2767,26 @@ namespace CoefClasses file.createAttribute("cols", HighFive::DataSpace::From(cols)).write(cols); } + bool TableData::CheckH5Params(HighFive::File& file) + { + int cols1; + + file.getAttribute("cols").read(cols1 ); + + int cols0 = coefs.begin()->second->cols; + + bool ret = true; + + if (cols0 != cols1) { + if (myid==0) + std::cout << "---- TableData::CheckH5Params: cols mismatch " << cols0 + << " != " << cols1 << std::endl; + ret = false; + } + + return ret; + } + unsigned TableData::WriteH5Times(HighFive::Group& snaps, unsigned count) { snaps.createDataSet("times", times); @@ -2899,7 +3134,15 @@ namespace CoefClasses // ReadH5Units(file); + // Check the specific parameters + // + if (not CheckH5Params(file)) { + throw std::runtime_error("Coefs::ExtendH5Coefs: " + "H5 parameter check failed, aborting extension"); + } + // Get the dataset + // HighFive::DataSet dataset = file.getDataSet("count"); unsigned count; @@ -3181,6 +3424,61 @@ namespace CoefClasses file.createAttribute ("dof", HighFive::DataSpace::From(dof) ).write(dof); } + bool SphFldCoefs::CheckH5Params(HighFive::File& file) + { + double scale0 = coefs.begin()->second->scale, scale1; + int nfld1, lmax1, nmax1, dof1; + + file.getAttribute("nfld" ).read(nfld1 ); + file.getAttribute("lmax" ).read(lmax1 ); + file.getAttribute("nmax" ).read(nmax1 ); + file.getAttribute("scale").read(scale1); + file.getAttribute("dof" ).read(dof1 ); + + int nfld0 = Nfld; + int lmax0 = Lmax; + int nmax0 = Nmax; + int dof0 = dof; + + bool ret = true; + + if (nfld0 != nfld1) { + if (myid==0) + std::cout << "---- SphFldCoefs::CheckH5Params: nfld mismatch " << nfld0 + << " != " << nfld1 << std::endl; + ret = false; + } + + if (lmax0 != lmax1) { + if (myid==0) + std::cout << "---- SphFldCoefs::CheckH5Params: lmax mismatch " << lmax0 + << " != " << lmax1 << std::endl; + ret = false; + } + if (nmax0 != nmax1) { + if (myid==0) + std::cout << "---- SphFldCoefs::CheckH5Params: nmax mismatch " << nmax0 + << " != " << nmax1 << std::endl; + ret = false; + } + + if (scale0 != scale1) { + if (myid==0) + std::cout << "---- SphFldCoefs::CheckH5Params: scale mismatch " << scale0 + << " != " << scale1 << std::endl; + ret = false; + } + + if (dof0 != dof1) { + if (myid==0) + std::cout << "---- SphFldCoefs::CheckH5Params: dof mismatch " << dof0 + << " != " << dof1 << std::endl; + ret = false; + } + + return ret; + } + unsigned SphFldCoefs::WriteH5Times(HighFive::Group& snaps, unsigned count) { for (auto c : coefs) { @@ -3306,6 +3604,62 @@ namespace CoefClasses file.createAttribute ("dof", HighFive::DataSpace::From(dof) ).write(dof); } + bool CylFldCoefs::CheckH5Params(HighFive::File& file) + { + double scale0 = coefs.begin()->second->scale, scale1; + int nfld1, mmax1, nmax1, dof1; + + file.getAttribute("nfld" ).read(nfld1 ); + file.getAttribute("mmax" ).read(mmax1 ); + file.getAttribute("nmax" ).read(nmax1 ); + file.getAttribute("scale").read(scale1); + file.getAttribute("dof" ).read(dof1 ); + + int nfld0 = Nfld; + int mmax0 = Mmax; + int nmax0 = Nmax; + int dof0 = dof; + + bool ret = true; + + if (nfld0 != nfld1) { + if (myid==0) + std::cout << "---- CylFldCoefs::CheckH5Params: nfld mismatch " << nfld0 + << " != " << nfld1 << std::endl; + ret = false; + } + + if (mmax0 != mmax1) { + if (myid==0) + std::cout << "---- CylFldCoefs::CheckH5Params: mmax mismatch " << mmax0 + << " != " << mmax1 << std::endl; + ret = false; + } + + if (nmax0 != nmax1) { + if (myid==0) + std::cout << "---- CylFldCoefs::CheckH5Params: nmax mismatch " << nmax0 + << " != " << nmax1 << std::endl; + ret = false; + } + + if (scale0 != scale1) { + if (myid==0) + std::cout << "---- CylFldCoefs::CheckH5Params: scale mismatch " << scale0 + << " != " << scale1 << std::endl; + ret = false; + } + + if (dof0 != dof1) { + if (myid==0) + std::cout << "---- CylFldCoefs::CheckH5Params: dof mismatch " << dof0 + << " != " << dof1 << std::endl; + ret = false; + } + + return ret; + } + unsigned CylFldCoefs::WriteH5Times(HighFive::Group& snaps, unsigned count) { for (auto c : coefs) { diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc index f51500470..1164ebf6b 100644 --- a/pyEXP/CoefWrappers.cc +++ b/pyEXP/CoefWrappers.cc @@ -144,6 +144,10 @@ void CoefficientClasses(py::module &m) { PYBIND11_OVERRIDE_PURE(void, Coefs, WriteH5Params, file); } + bool CheckH5Params(HighFive::File& file) override { + PYBIND11_OVERRIDE_PURE(bool, Coefs, CheckH5Params, file); + } + unsigned WriteH5Times(HighFive::Group& group, unsigned count) override { PYBIND11_OVERRIDE_PURE(unsigned, Coefs, WriteH5Times, group, count); } @@ -233,6 +237,10 @@ void CoefficientClasses(py::module &m) { PYBIND11_OVERRIDE(void, SphCoefs, WriteH5Params, file); } + bool CheckH5Params(HighFive::File& file) override { + PYBIND11_OVERRIDE(bool, SphCoefs, CheckH5Params, file); + } + unsigned WriteH5Times(HighFive::Group& group, unsigned count) override { PYBIND11_OVERRIDE(unsigned, SphCoefs, WriteH5Times, group, count); } @@ -312,6 +320,10 @@ void CoefficientClasses(py::module &m) { PYBIND11_OVERRIDE(void, CylCoefs, WriteH5Params, file); } + bool CheckH5Params(HighFive::File& file) override { + PYBIND11_OVERRIDE(bool, CylCoefs, CheckH5Params, file); + } + unsigned WriteH5Times(HighFive::Group& group, unsigned count) override { PYBIND11_OVERRIDE(unsigned, CylCoefs, WriteH5Times, group, count); } @@ -390,6 +402,10 @@ void CoefficientClasses(py::module &m) { PYBIND11_OVERRIDE(void, SlabCoefs, WriteH5Params, file); } + bool CheckH5Params(HighFive::File& file) override { + PYBIND11_OVERRIDE(bool, SlabCoefs, CheckH5Params, file); + } + unsigned WriteH5Times(HighFive::Group& group, unsigned count) override { PYBIND11_OVERRIDE(unsigned, SlabCoefs, WriteH5Times, group, count); } @@ -469,6 +485,10 @@ void CoefficientClasses(py::module &m) { PYBIND11_OVERRIDE(void, CubeCoefs, WriteH5Params, file); } + bool CheckH5Params(HighFive::File& file) override { + PYBIND11_OVERRIDE(bool, CubeCoefs, CheckH5Params, file); + } + unsigned WriteH5Times(HighFive::Group& group, unsigned count) override { PYBIND11_OVERRIDE(unsigned, CubeCoefs, WriteH5Times, group, count); } @@ -548,6 +568,10 @@ void CoefficientClasses(py::module &m) { PYBIND11_OVERRIDE(void, TableData, WriteH5Params, file); } + bool CheckH5Params(HighFive::File& file) override { + PYBIND11_OVERRIDE(bool, TableData, CheckH5Params, file); + } + unsigned WriteH5Times(HighFive::Group& group, unsigned count) override { PYBIND11_OVERRIDE(unsigned, TableData, WriteH5Times, group, count); } @@ -626,6 +650,10 @@ void CoefficientClasses(py::module &m) { PYBIND11_OVERRIDE(void, TrajectoryData, WriteH5Params, file); } + bool CheckH5Params(HighFive::File& file) override { + PYBIND11_OVERRIDE(bool, TrajectoryData, CheckH5Params, file); + } + unsigned WriteH5Times(HighFive::Group& group, unsigned count) override { PYBIND11_OVERRIDE(unsigned, TrajectoryData, WriteH5Times, group, count); }