diff --git a/CMakeLists.txt b/CMakeLists.txt index e5c7cf6ad..9ed3fa60c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.25) # Needed for CUDA, MPI, and CTest features project( EXP - VERSION "7.10.0" + VERSION "7.10.1" HOMEPAGE_URL https://github.com/EXP-code/EXP LANGUAGES C CXX Fortran) diff --git a/doc/exp.cfg b/doc/exp.cfg index c83ce4a1c..d6a169162 100644 --- a/doc/exp.cfg +++ b/doc/exp.cfg @@ -48,7 +48,7 @@ PROJECT_NAME = EXP # could be handy for archiving the generated documentation or if some version # control system is used. -PROJECT_NUMBER = 7.10.0 +PROJECT_NUMBER = 7.10.1 # Using the PROJECT_BRIEF tag one can provide an optional one line description # for a project that appears at the top of each page and should give viewer a diff --git a/doc/exp.cfg.breathe b/doc/exp.cfg.breathe index c8d8a0781..95375995d 100644 --- a/doc/exp.cfg.breathe +++ b/doc/exp.cfg.breathe @@ -48,7 +48,7 @@ PROJECT_NAME = EXP # could be handy for archiving the generated documentation or if some version # control system is used. -PROJECT_NUMBER = 7.10.0 +PROJECT_NUMBER = 7.10.1 # Using the PROJECT_BRIEF tag one can provide an optional one line description # for a project that appears at the top of each page and should give viewer a diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 9c6d21a30..0367512c5 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -513,16 +513,15 @@ namespace BasisClasses cf->store.resize(ldim*nmax); // Make the coefficient map - cf->coefs = std::make_shared - (cf->store.data(), ldim, nmax); + cf->initCoefMap(ldim, nmax); for (int l=0, L0=0, L1=0; l<=lmax; l++) { for (int m=0; m<=l; m++) { for (int n=0; ncoefs)(L0, n) = {expcoef(L1, n), 0.0}; + cf->getCoefs()(L0, n) = {expcoef(L1, n), 0.0}; else - (*cf->coefs)(L0, n) = {expcoef(L1, n), expcoef(L1+1, n)}; + cf->getCoefs()(L0, n) = {expcoef(L1, n), expcoef(L1+1, n)}; } L0 += 1; if (m==0) L1 += 1; @@ -543,8 +542,8 @@ namespace BasisClasses { auto & cf = *dynamic_cast(coef.get()); - int rows = (*cf.coefs).rows(); - int cols = (*cf.coefs).cols(); + int rows = cf.getCoefs().rows(); + int cols = cf.getCoefs().cols(); int rexp = (lmax+1)*(lmax+2)/2; if (rows != rexp or cols != nmax) { std::ostringstream sout; @@ -575,10 +574,10 @@ namespace BasisClasses for (int m=0; m<=l; m++) { for (int n=0; ncoefs)(L0, n).real(); + expcoef(L1, n) = cf->getCoefs()(L0, n).real(); else { - expcoef(L1, n) = (*cf->coefs)(L0, n).real(); - expcoef(L1+1, n) = (*cf->coefs)(L0, n).imag(); + expcoef(L1, n) = cf->getCoefs()(L0, n).real(); + expcoef(L1+1, n) = cf->getCoefs()(L0, n).imag(); } } L0 += 1; @@ -2266,14 +2265,13 @@ namespace BasisClasses cf->store.resize((mmax+1)*nmax); // Create a new instance - cf->coefs = std::make_shared - (cf->store.data(), mmax+1, nmax); + cf->initCoefMap(mmax+1, nmax); for (int m=0; m<=mmax; m++) { sl->get_coefs(m, cos1, sin1); for (int n=0; ncoefs)(m, n) = {cos1(n), sin1(n)}; + cf->getCoefs()(m, n) = {cos1(n), sin1(n)}; } } } @@ -2294,7 +2292,7 @@ namespace BasisClasses coefret = coef; for (int m=0; m<=mmax; m++) { // Set to zero on m=0 call only--------+ - sl->set_coefs(m, (*cf->coefs).row(m).real(), (*cf->coefs).row(m).imag(), m==0); + sl->set_coefs(m, cf->getCoefs().row(m).real(), cf->getCoefs().row(m).imag(), m==0); } // Assign center if need be @@ -2612,15 +2610,14 @@ namespace BasisClasses cf->store.resize((mmax+1)*nmax); // Make the coefficient map - cf->coefs = std::make_shared - (cf->store.data(), mmax+1, nmax); + cf->initCoefMap(mmax+1, nmax); for (int m=0, m0=0; m<=mmax; m++) { for (int n=0; ncoefs)(m, n) = {expcoef(m0, n), 0.0}; + cf->getCoefs()(m, n) = {expcoef(m0, n), 0.0}; else - (*cf->coefs)(m, n) = {expcoef(m0, n), expcoef(m0+1, n)}; + cf->getCoefs()(m, n) = {expcoef(m0, n), expcoef(m0+1, n)}; } if (m==0) m0 += 1; else m0 += 2; @@ -2638,9 +2635,9 @@ namespace BasisClasses // { auto cc = dynamic_cast(coef.get()); - auto cf = cc->coefs; - int rows = cf->rows(); - int cols = cf->cols(); + auto cf = cc->getCoefs(); + int rows = cf.rows(); + int cols = cf.cols(); if (rows != mmax+1 or cols != nmax) { std::ostringstream sout; sout << "FlatDisk::set_coefs: the basis has (mmax+1, nmax)=(" @@ -2653,7 +2650,7 @@ namespace BasisClasses } CoefClasses::CylStruct* cf = dynamic_cast(coef.get()); - auto & cc = *cf->coefs; + auto & cc = cf->getCoefs(); // Set gravitational constant // @@ -3458,15 +3455,14 @@ namespace BasisClasses cf->store.resize((mmax+1)*nmax); // Make the coefficient map - cf->coefs = std::make_shared - (cf->store.data(), mmax+1, nmax); + cf->initCoefMap(mmax+1, nmax); for (int m=0, m0=0; m<=mmax; m++) { for (int n=0; ncoefs)(m, n) = {expcoef(m0, n), 0.0}; + cf->getCoefs()(m, n) = {expcoef(m0, n), 0.0}; else - (*cf->coefs)(m, n) = {expcoef(m0, n), expcoef(m0+1, n)}; + cf->getCoefs()(m, n) = {expcoef(m0, n), expcoef(m0+1, n)}; } if (m==0) m0 += 1; else m0 += 2; @@ -3484,9 +3480,9 @@ namespace BasisClasses // { auto cc = dynamic_cast(coef.get()); - auto cf = cc->coefs; - int rows = cf->rows(); - int cols = cf->cols(); + auto cf = cc->getCoefs(); + int rows = cf.rows(); + int cols = cf.cols(); if (rows != mmax+1 or cols != nmax) { std::ostringstream sout; sout << "CBDisk::set_coefs: the basis has (mmax+1, nmax)=(" @@ -3499,7 +3495,7 @@ namespace BasisClasses } CoefClasses::CylStruct* cf = dynamic_cast(coef.get()); - auto & cc = *cf->coefs; + auto & cc = cf->getCoefs(); // Set gravitational constant // @@ -3966,8 +3962,7 @@ namespace BasisClasses G = cf->getGravConstant(); cf->allocate(); - - *cf->coefs = expcoef; + cf->setCoefs(expcoef); } void Slab::set_coefs(CoefClasses::CoefStrPtr coef) @@ -3981,7 +3976,7 @@ namespace BasisClasses // { auto cc = dynamic_cast(coef.get()); - auto d = cc->coefs->dimensions(); + auto d = cc->getCoefs().dimensions(); if (d[0] != 2*nmaxx+1 or d[1] != 2*nmaxy+1 or d[2] != nmaxz) { std::ostringstream sout; sout << "Slab::set_coefs: the basis has (2*nmaxx+1, 2*nmaxy+1, nmaxz)=(" @@ -3996,7 +3991,7 @@ namespace BasisClasses } auto cf = dynamic_cast(coef.get()); - expcoef = *cf->coefs; + expcoef = cf->getCoefs(); // Set gravitational constant // @@ -4540,8 +4535,7 @@ namespace BasisClasses G = cf->getGravConstant(); cf->allocate(); - - *cf->coefs = expcoef; + cf->setCoefs(expcoef); } void Cube::set_coefs(CoefClasses::CoefStrPtr coef) @@ -4555,7 +4549,7 @@ namespace BasisClasses // { auto cc = dynamic_cast(coef.get()); - auto d = cc->coefs->dimensions(); + auto d = cc->getCoefs().dimensions(); if (d[0] != 2*nmaxx+1 or d[1] != 2*nmaxy+1 or d[2] != 2*nmaxz+1) { std::ostringstream sout; sout << "Cube::set_coefs: the basis has (2*nmaxx+1, 2*nmaxy+1, 2*nmaxz+1)=(" @@ -4570,7 +4564,7 @@ namespace BasisClasses } auto cf = dynamic_cast(coef.get()); - expcoef = *cf->coefs; + expcoef = cf->getCoefs(); // Set gravitational constant // diff --git a/expui/CoefContainer.cc b/expui/CoefContainer.cc index 9dda72b9d..55c49e5d6 100644 --- a/expui/CoefContainer.cc +++ b/expui/CoefContainer.cc @@ -188,16 +188,16 @@ namespace MSSA for (auto k : keys) { if (k[2]==0) - data[k][t] = (*cf->coefs)(k[0], k[1]).real(); + data[k][t] = (cf->getCoefs())(k[0], k[1]).real(); else - data[k][t] = (*cf->coefs)(k[0], k[1]).imag(); + data[k][t] = (cf->getCoefs())(k[0], k[1]).imag(); } for (auto k : bkeys) { if (k[2]==0) - data[k][t] = (*cf->coefs)(k[0], k[1]).real(); + data[k][t] = (cf->getCoefs())(k[0], k[1]).real(); else - data[k][t] = (*cf->coefs)(k[0], k[1]).imag(); + data[k][t] = (cf->getCoefs())(k[0], k[1]).imag(); } } } @@ -214,8 +214,8 @@ namespace MSSA int m = k[0], n = k[1]; - if (m==0) (*cf->coefs)(m, n) = {data[c][i], 0.0}; - else (*cf->coefs)(m, n) = {data[c][i], data[s][i]}; + if (m==0) (cf->getCoefs())(m, n) = {data[c][i], 0.0}; + else (cf->getCoefs())(m, n) = {data[c][i], data[s][i]}; } // END key loop } @@ -278,9 +278,9 @@ namespace MSSA for (auto k : keys) { if (k[3]==0) - data[k][t] = (*cf->coefs)(k[0], k[1], k[2]).real(); + data[k][t] = (cf->getCoefs())(k[0], k[1], k[2]).real(); else - data[k][t] = (*cf->coefs)(k[0], k[1], k[3]).imag(); + data[k][t] = (cf->getCoefs())(k[0], k[1], k[2]).imag(); } } } @@ -295,10 +295,10 @@ namespace MSSA auto c = k, s = k; c.push_back(0); s.push_back(1); - int f = k[1], m = k[1], n = k[2]; + int f = k[0], m = k[1], n = k[2]; - if (m==0) (*cf->coefs)(f, m, n) = {data[c][i], 0.0}; - else (*cf->coefs)(f, m, n) = {data[c][i], data[s][i]}; + if (m==0) (cf->getCoefs())(f, m, n) = {data[c][i], 0.0}; + else (cf->getCoefs())(f, m, n) = {data[c][i], data[s][i]}; } // END key loop } @@ -381,13 +381,13 @@ namespace MSSA ( cur->getCoefStruct(times[t]).get() ); for (auto k : keys) { - auto c = (*cf->coefs)(I(k), k[2]); + auto c = (cf->getCoefs())(I(k), k[2]); data[k][t] = c.real(); if (k[3]) data[k][t] = c.imag(); } for (auto k : bkeys) { - auto c = (*cf->coefs)(I(k), k[2]); + auto c = (cf->getCoefs())(I(k), k[2]); data[k][t] = c.real(); if (k[3]) data[k][t] = c.imag(); } @@ -410,8 +410,8 @@ namespace MSSA int m = k[1], n = k[2]; - if (m==0) (*cf->coefs)(I(k), n) = {data[c][i], 0.0 }; - else (*cf->coefs)(I(k), n) = {data[c][i], data[s][i]}; + if (m==0) (cf->getCoefs())(I(k), n) = {data[c][i], 0.0 }; + else (cf->getCoefs())(I(k), n) = {data[c][i], data[s][i]}; } // END key loop } @@ -478,7 +478,7 @@ namespace MSSA ( cur->getCoefStruct(times[t]).get() ); for (auto k : keys) { - auto c = (*cf->coefs)(k[0], I(k), k[3]); + auto c = (cf->getCoefs())(k[0], I(k), k[3]); data[k][t] = c.real(); if (k[4]) data[k][t] = c.imag(); } @@ -503,8 +503,8 @@ namespace MSSA int f = k[0], l = k[1], m = k[2], n = k[3]; - if (m==0) (*cf->coefs)(f, I(k), n) = {data[c][i], 0.0 }; - else (*cf->coefs)(f, I(k), n) = {data[c][i], data[s][i]}; + if (m==0) (cf->getCoefs())(f, I(k), n) = {data[c][i], 0.0 }; + else (cf->getCoefs())(f, I(k), n) = {data[c][i], data[s][i]}; } // END key loop } @@ -582,13 +582,13 @@ namespace MSSA ( cur->getCoefStruct(times[t]).get() ); for (auto k : keys) { - auto c = (*cf->coefs)(k[0], k[1], k[2]); + auto c = (cf->getCoefs())(k[0], k[1], k[2]); if (k[3]) data[k][t] = c.imag(); else data[k][t] = c.real(); } for (auto k : bkeys) { - auto c = (*cf->coefs)(k[0], k[1], k[2]); + auto c = (cf->getCoefs())(k[0], k[1], k[2]); if (k[3]) data[k][t] = c.imag(); else data[k][t] = c.real(); } @@ -666,13 +666,13 @@ namespace MSSA ( cur->getCoefStruct(times[t]).get() ); for (auto k : keys) { - auto c = (*cf->coefs)(k[0], k[1], k[2]); + auto c = (cf->getCoefs())(k[0], k[1], k[2]); if (k[3]) data[k][t] = c.imag(); else data[k][t] = c.real(); } for (auto k : bkeys) { - auto c = (*cf->coefs)(k[0], k[1], k[2]); + auto c = (cf->getCoefs())(k[0], k[1], k[2]); if (k[3]) data[k][t] = c.imag(); else data[k][t] = c.real(); } @@ -691,7 +691,7 @@ namespace MSSA c.push_back(0); s.push_back(1); - (*cf->coefs)(k[0], k[1], k[2]) = {data[c][i], data[s][i]}; + (cf->getCoefs())(k[0], k[1], k[2]) = {data[c][i], data[s][i]}; } // END key loop } @@ -710,7 +710,7 @@ namespace MSSA c.push_back(0); s.push_back(1); - (*cf->coefs)(k[0], k[1], k[2]) = {data[c][i], data[s][i]}; + (cf->getCoefs())(k[0], k[1], k[2]) = {data[c][i], data[s][i]}; } // END key loop } @@ -765,7 +765,7 @@ namespace MSSA cf = dynamic_cast ( cur->getCoefStruct(times[t]).get() ); - data[key][t] = (*cf->coefs)(c).real(); + data[key][t] = (cf->getCoefs())(c).real(); } } } @@ -781,7 +781,7 @@ namespace MSSA for (unsigned c=0; ccoefs)(c) = data[key][i]; + (cf->getCoefs())(c) = data[key][i]; } // End field loop } @@ -840,7 +840,7 @@ namespace MSSA ( cur->getCoefStruct(times[t]).get() ); for (int n=0; ncoefs)(m, n).real(); + data[key][t*rank+n] = (cf->getCoefs())(m, n).real(); } } @@ -864,7 +864,7 @@ namespace MSSA for (unsigned m=0; mcoefs)(m, n) = data[key][t*rank+n]; + (cf->getCoefs())(m, n) = data[key][t*rank+n]; } // End field loop } @@ -982,7 +982,7 @@ namespace MSSA (cur->getCoefStruct(times[t]).get()); // Get the coefficient map - auto & ar = *(cf->coefs); + auto & ar = cf->getCoefs(); for (auto k : bkeys) { auto c = ar(I(k), k[2]); @@ -1003,7 +1003,7 @@ namespace MSSA (cur->getCoefStruct(times[t]).get()); // Get the coefficient map - auto & ar = *(cf->coefs); + auto & ar = cf->getCoefs(); for (auto k : bkeys) { if (k[2]==0) diff --git a/expui/CoefStruct.H b/expui/CoefStruct.H index 9f9468790..8a65d975c 100644 --- a/expui/CoefStruct.H +++ b/expui/CoefStruct.H @@ -119,9 +119,6 @@ namespace CoefClasses //! Coefficient map type using coefType = Eigen::Map; - //! Coefficient map - std::shared_ptr coefs; - //! Angular and radial dimension int lmax, nmax; @@ -130,7 +127,12 @@ namespace CoefClasses //! Is basis unit normed? bool normed; - + + private: + //! Cached coefficient map (mutable for lazy initialization) + mutable std::shared_ptr coefs; + + public: //! Constructor SphStruct(Coefs* C=nullptr) : lmax(0), nmax(0), scale(1.0), normed(true), CoefStruct(C) @@ -145,13 +147,36 @@ namespace CoefClasses //! Copy all of data contents to a new instance std::shared_ptr deepcopy(); + //! Lazy-create and cache the coefficient map + const coefType& getCoefMap() const + { + if (!coefs) { + int rows = (lmax+1)*(lmax+2)/2; + int cols = nmax; + coefs = std::make_shared + (const_cast*>(store.data()), rows, cols); + } + return *coefs; + } + + //! Non-const version for assignment + coefType& getCoefMapMutable() + { + if (!coefs) { + int rows = (lmax+1)*(lmax+2)/2; + int cols = nmax; + coefs = std::make_shared(store.data(), rows, cols); + } + return *coefs; + } + //! Allocate storage arrays void allocate() { int rows = (lmax+1)*(lmax+2)/2; int cols = nmax; + coefs.reset(); // Clear map BEFORE resize store.resize(rows*cols); - coefs = std::make_shared(store.data(), rows, cols); } //! Assign matrix @@ -160,11 +185,35 @@ namespace CoefClasses lmax = Lmax; nmax = Nmax; allocate(); - *coefs = mat; + getCoefMapMutable() = mat; } + //! Initialize coefficient map with dimensions + void initCoefMap(int rows, int cols) + { + coefs.reset(); + coefs = std::make_shared(store.data(), rows, cols); + } + + //! Set the underlying coefficient data + void setCoefs(const Eigen::MatrixXcd& mat) + { + int rows = (lmax+1)*(lmax+2)/2; + int cols = nmax; + if (mat.rows() != rows || mat.cols() != cols) + throw std::invalid_argument("setCoefs: size mismatch"); + + getCoefMapMutable() = mat; + } + + //@{ + //! Direct access to underlying Eigen Map (for client code) + const coefType& getCoefs() const { return getCoefMap(); } + coefType& getCoefs() { return getCoefMapMutable(); } + //@} + }; - + //! Specialization of CoefStruct for cylinders class CylStruct : public CoefStruct { @@ -172,12 +221,14 @@ namespace CoefClasses //! Coefficient map type using coefType = Eigen::Map; - //! Coefficient map - std::shared_ptr coefs; - //! Angular and radial dimension int mmax, nmax; - + + private: + //! Cached coefficient map (mutable for lazy initialization) + mutable std::shared_ptr coefs; + + public: //! Constructor CylStruct(Coefs* C=nullptr) : mmax(0), nmax(0), CoefStruct(C) { geom = "cylinder"; } @@ -191,13 +242,35 @@ namespace CoefClasses //! Copy all of data contents to a new instance std::shared_ptr deepcopy(); + //! Lazy-create and cache the coefficient map + const coefType& getCoefMap() const + { + if (!coefs) { + int rows = mmax+1; + int cols = nmax; + coefs->resize(rows, cols); + } + return *coefs; + } + + //! Non-const version for assignment + coefType& getCoefMapMutable() + { + int rows = mmax+1; + int cols = nmax; + coefs.reset(); // Clear old map before creating new one + coefs = std::make_shared(store.data(), rows, cols); + return *coefs; + } + //! Allocate storage arrays void allocate() { int rows = mmax+1; int cols = nmax; + + coefs.reset(); // Clear map BEFORE resize store.resize(rows*cols); - coefs = std::make_shared(store.data(), rows, cols); } //! Assign matrix @@ -206,8 +279,33 @@ namespace CoefClasses mmax = Mmax; nmax = Nmax; allocate(); - *coefs = mat; + getCoefMapMutable() = mat; } + + //! Initialize coefficient map with dimensions + void initCoefMap(int rows, int cols) + { + coefs.reset(); + coefs = std::make_shared(store.data(), rows, cols); + } + + //! Set the underlying coefficient data + void setCoefs(const Eigen::MatrixXcd& mat) + { + int rows = mmax+1; + int cols = nmax; + + if (mat.rows() != rows || mat.cols() != cols) + throw std::invalid_argument("setCoefs: size mismatch"); + getCoefMapMutable() = mat; + } + + //@{ + //! Direct access to underlying Eigen Map (for client code) + const coefType& getCoefs() const { return getCoefMap(); } + coefType& getCoefs() { return getCoefMapMutable(); } + //@} + }; //! Specialization of CoefStruct for a slab @@ -217,15 +315,17 @@ namespace CoefClasses //! Coefficient map type using coefType = Eigen::TensorMap, 3>>; - //! Coefficient map - std::shared_ptr coefs; - //! Angular and radial dimension int nmaxx, nmaxy, nmaxz; //! Basis dimensions int nx, ny, nz, dim; + private: + //! Cached coefficient map (mutable for lazy initialization) + mutable std::shared_ptr coefs; + + public: //! Constructor SlabStruct(Coefs* C=nullptr) : nmaxx(0), nmaxy(0), nmaxz(0), CoefStruct(C) { geom = "slab"; } @@ -239,6 +339,30 @@ namespace CoefClasses //! Copy all of data contents to a new instance std::shared_ptr deepcopy(); + +// For Map-based classes: +/**/ + + + //! Lazy-create and cache the coefficient tensor map + const coefType& getCoefTensorMap() const + { + if (!coefs) { + coefs = std::make_shared + (const_cast*>(store.data()), nx, ny, nz); + } + return *coefs; + } + + + //! Non-const version for assignment + coefType& getCoefTensorMapMutable() + { + coefs.reset(); // Clear old map before creating new one + coefs = std::make_shared(store.data(), nx, ny, nz); + return *coefs; + } + //! Allocate storage arrays void allocate() { @@ -246,11 +370,9 @@ namespace CoefClasses ny = 2*nmaxy + 1; nz = nmaxz; dim = nx * ny * nz; + coefs.reset(); // Clear map BEFORE resize store.resize(dim); - coefs = std::make_shared(store.data(), nx, ny, nz); } - - //! Assign matrix void assign(const Eigen::Tensor, 3>& dat) { @@ -258,15 +380,32 @@ namespace CoefClasses nx = d[0]; ny = d[1]; nz = d[2]; - nmaxx = (nx - 1)/2; nmaxy = (ny - 1)/2; nmaxz = nz; - allocate(); - *coefs = dat; + getCoefTensorMapMutable() = dat; + } + + //! Initialize coefficient tensor map with dimensions + void initCoefMap(int nx, int ny, int nz) + { + coefs.reset(); + coefs = std::make_shared(store.data(), nx, ny, nz); + } + + //! Set the underlying coefficient data + void setCoefs(const Eigen::Tensor, 3>& tensor) + { + const auto& d = tensor.dimensions(); + if (d[0] != nx || d[1] != ny || d[2] != nz) + throw std::invalid_argument("setCoefs: dimension mismatch"); + getCoefTensorMapMutable() = tensor; } + //! Direct access to underlying Eigen TensorMap (for client code) + const coefType& getCoefs() const { return getCoefTensorMap(); } + coefType& getCoefs() { return getCoefTensorMapMutable(); } }; //! Specialization of CoefStruct for a cube @@ -276,16 +415,16 @@ namespace CoefClasses //! Array dimensions int nx, ny, nz, dim; + //! Cached coefficient map (mutable for lazy initialization) + mutable std::shared_ptr, 3>>> coefs; + public: //! Coefficient map type using coefType = Eigen::TensorMap, 3>>; - //! Coefficient map - std::shared_ptr coefs; - //! Angular and radial dimension int nmaxx, nmaxy, nmaxz; - + //! Constructor CubeStruct(Coefs* C=nullptr) : nmaxx(0), nmaxy(0), nmaxz(0), nx(0), ny(0), nz(0), dim(0), CoefStruct(C) @@ -300,6 +439,24 @@ namespace CoefClasses //! Copy all of data contents to a new instance std::shared_ptr deepcopy(); + //! Lazy-create and cache the coefficient tensor map + const coefType& getCoefTensorMap() const + { + if (!coefs) { + coefs = std::make_shared + (const_cast*>(store.data()), nx, ny, nz); + } + return *coefs; + } + + //! Non-const version for assignment + coefType& getCoefTensorMapMutable() + { + coefs.reset(); // Clear old map before creating new one + coefs = std::make_shared(store.data(), nx, ny, nz); + return *coefs; + } + //! Allocate storage arrays void allocate() { @@ -310,15 +467,14 @@ namespace CoefClasses dim = nx * ny * nz; // Allocate + coefs.reset(); // Clear map BEFORE resize store.resize(dim); - coefs = std::make_shared(store.data(), nx, ny, nz); } //! Assign matrix void assign(const Eigen::Tensor, 3>& dat) { const auto& d = dat.dimensions(); - // Assign grid dimensions nx = d[0]; ny = d[1]; @@ -329,13 +485,31 @@ namespace CoefClasses nmaxx = (d[0]-1)/2; nmaxy = (d[1]-1)/2; nmaxz = (d[2]-1)/2; - // Allocate storage and assign allocate(); - *coefs = dat; + getCoefTensorMapMutable() = dat; + } + + //! Initialize coefficient tensor map with dimensions + void initCoefMap(int nx, int ny, int nz) + { + coefs.reset(); + coefs = std::make_shared(store.data(), nx, ny, nz); } + + //! Set the underlying coefficient data + void setCoefs(const Eigen::Tensor, 3>& tensor) + { + const auto& d = tensor.dimensions(); + if (d[0] != nx || d[1] != ny || d[2] != nz) + throw std::invalid_argument("setCoefs: dimension mismatch"); + getCoefTensorMapMutable() = tensor; + } + + //! Direct access to underlying Eigen TensorMap (for client code) + const coefType& getCoefs() const { return getCoefTensorMap(); } + coefType& getCoefs() { return getCoefTensorMapMutable(); } }; - //! Specialization of CoefStruct for a data table class TblStruct : public CoefStruct @@ -344,12 +518,14 @@ namespace CoefClasses //! Coefficient map type using coefType = Eigen::Map; - //! Coefficient map - std::shared_ptr coefs; - //! Number of columns int cols; - + + private: + //! Cached coefficient map (mutable for lazy initialization) + mutable std::shared_ptr coefs; + + public: //! Constructor TblStruct(Coefs* C=nullptr) : cols(0), CoefStruct(C) { geom = "table"; } @@ -362,23 +538,62 @@ namespace CoefClasses //! Copy all of data contents to a new instance std::shared_ptr deepcopy(); + //! Lazy-create and cache the coefficient map + const coefType& getCoefMap() const + { + if (!coefs) { + coefs = std::make_shared + (const_cast*>(store.data()), store.size()); + } + return *coefs; + } + + //! Non-const version for assignment + coefType& getCoefMapMutable() + { + coefs.reset(); // Clear old map before creating new one + coefs = std::make_shared(store.data(), store.size()); + return *coefs; + } + //! Assign array void assign(const Eigen::VectorXcd& arr) { - cols = arr.size(); + cols = arr.size(); + coefs.reset(); // Clear map BEFORE resize store = arr; - coefs = std::make_shared(store.data(), store.size()); } //! Allocate storage arrays void allocate() { + coefs.reset(); // Clear map BEFORE resize store.resize(cols); - coefs = std::make_shared(store.data(), cols); } + //! Initialize coefficient map + void initCoefMap() + { + coefs.reset(); + coefs = std::make_shared(store.data(), store.size()); + } + + //! Set the underlying coefficient data + void setCoefs(const Eigen::VectorXcd& vec) + { + if (vec.size() != cols) + throw std::invalid_argument("setCoefs: size mismatch"); + getCoefMapMutable() = vec; + } + + //@{ + //! Direct access to underlying Eigen Map (for client code) + const coefType& getCoefs() const { return getCoefMap(); } + coefType& getCoefs() { return getCoefMapMutable(); } + //@} + }; - + //! Specialization of CoefStruct for collection of trajectories class TrajStruct : public CoefStruct { @@ -386,15 +601,17 @@ namespace CoefClasses //! Coefficient map type using coefType = Eigen::Map; - //! Coefficient map - std::shared_ptr coefs; - //! Number of trajectories int traj; //! Phase-space rank int rank; - + +private: + //! Cached coefficient map (mutable for lazy initialization) + mutable std::shared_ptr coefs; + + public: //! Constructor TrajStruct(Coefs* C=nullptr) : traj(0), rank(0), CoefStruct(C) { geom = "trajectory"; } @@ -408,23 +625,65 @@ namespace CoefClasses //! Copy all of data contents to a new instance std::shared_ptr deepcopy(); + //! Lazy-create and cache the coefficient map + const coefType& getCoefMap() const + { + if (!coefs) { + coefs = std::make_shared + (const_cast*>(store.data()), traj, rank); + } + return *coefs; + } + + //! Non-const version for assignment + coefType& getCoefMapMutable() + { + coefs.reset(); // Clear old map before creating new one + coefs = std::make_shared(store.data(), traj, rank); + + return *coefs; + } + //! Assign array void assign(const Eigen::MatrixXd& arr) { traj = arr.rows(); rank = arr.cols(); - + allocate(); - *coefs = arr; + getCoefMapMutable() = arr; } //! Allocate storage arrays void allocate() { + coefs.reset(); // Clear map BEFORE resize store.resize(traj*rank); - coefs = std::make_shared(store.data(), traj, rank); } + //! Initialize coefficient map with dimensions + void initCoefMap(int rows, int cols) + { + coefs.reset(); + coefs = std::make_shared(store.data(), rows, cols); + } + + //! Set the underlying coefficient data + void setCoefs(const Eigen::MatrixXcd& mat) + { + traj = mat.rows(); + rank = mat.cols(); + + allocate(); + getCoefMapMutable() = mat; + } + + //@{ + //! Direct access to underlying Eigen Map (for client code) + const coefType& getCoefs() const { return getCoefMap(); } + coefType& getCoefs() { return getCoefMapMutable(); } + //@} + }; //! Specialization of CoefStruct for spheres for general fields @@ -437,18 +696,19 @@ namespace CoefClasses //! Coefficient map type using coefType = Eigen::TensorMap; - //! Coefficient map - std::shared_ptr coefs; - //! Number of data fields int nfld; - //! Angular and radial dimension int lmax, nmax; //! Radial scale factor double scale; - + + private: + //! Cached coefficient map (mutable for lazy initialization) + mutable std::shared_ptr coefs; + + public: //! Constructor SphFldStruct(Coefs* C=nullptr) : nfld(0), lmax(0), nmax(0), scale(1.0), CoefStruct(C) @@ -463,13 +723,35 @@ namespace CoefClasses //! Copy all of data contents to a new instance std::shared_ptr deepcopy(); + //! Lazy-create and cache the coefficient tensor map + const coefType& getCoefTensorMap() const + { + if (!coefs) { + int rows = (lmax+1)*(lmax+2)/2; + int cols = nmax; + coefs = std::make_shared + (const_cast*>(store.data()), nfld, rows, cols); + } + return *coefs; + } + + //! Non-const version for assignment + coefType& getCoefTensorMapMutable() + { + int rows = (lmax+1)*(lmax+2)/2; + int cols = nmax; + coefs.reset(); // Clear old map before creating new one + coefs = std::make_shared(store.data(), nfld, rows, cols); + return *coefs; + } + //! Allocate storage arrays void allocate() { int rows = (lmax+1)*(lmax+2)/2; int cols = nmax; + coefs.reset(); // Clear map BEFORE resize store.resize(nfld*rows*cols); - coefs = std::make_shared(store.data(), nfld, rows, cols); } //! Assign data @@ -482,8 +764,30 @@ namespace CoefClasses store = data; } + //! Initialize coefficient tensor map with dimensions + void initCoefMap(int nfld, int rows, int cols) + { + coefs.reset(); + coefs = std::make_shared(store.data(), nfld, rows, cols); + } + + //! Set the underlying coefficient data + void setCoefs(const Eigen::Tensor, 3>& tensor) + { + int rows = (lmax+1)*(lmax+2)/2; + int cols = nmax; + + const auto& d = tensor.dimensions(); + if (d[0] != nfld || d[1] != rows || d[2] != cols) + throw std::invalid_argument("setCoefs: dimension mismatch"); + + getCoefTensorMapMutable() = tensor; + } + + //! Direct access to underlying Eigen TensorMap (for client code) + const coefType& getCoefs() const { return getCoefTensorMap(); } + coefType& getCoefs() { return getCoefTensorMapMutable(); } }; - //! Specialization of CoefStruct for disks for general fields class CylFldStruct : public CoefStruct @@ -495,9 +799,6 @@ namespace CoefClasses //! Coefficient map type using coefType = Eigen::TensorMap; - //! Coefficient map - std::shared_ptr coefs; - //! Number of data fields int nfld; @@ -506,7 +807,12 @@ namespace CoefClasses //! Radial scale factor double scale; - + + private: + //! Cached coefficient map (mutable for lazy initialization) + mutable std::shared_ptr coefs; + + public: //! Constructor CylFldStruct(Coefs* C=nullptr) : nfld(0), mmax(0), nmax(0), scale(1.0), CoefStruct(C) @@ -521,15 +827,35 @@ namespace CoefClasses //! Copy all of data contents to a new instance std::shared_ptr deepcopy(); + //! Lazy-create and cache the coefficient tensor map + const coefType& getCoefTensorMap() const + { + if (!coefs) { + coefs = std::make_shared + (const_cast*>(store.data()), nfld, mmax+1, nmax); + } + return *coefs; + } + + //! Non-const version for assignment + coefType& getCoefTensorMapMutable() + { + int rows = mmax+1; + int cols = nmax; + coefs.reset(); // Clear old map before creating new one + coefs = std::make_shared(store.data(), nfld, rows, cols); + return *coefs; + } + //! Allocate storage arrays void allocate() { int rows = mmax+1; int cols = nmax; + coefs.reset(); // Clear map BEFORE resize store.resize(nfld*rows*cols); - coefs = std::make_shared(store.data(), nfld, rows, cols); } - + //! Assign matrix void assign(const Eigen::VectorXcd& data, int Nfld, int Mmax, int Nmax) { @@ -540,6 +866,29 @@ namespace CoefClasses store = data; } + //! Initialize coefficient tensor map with dimensions + void initCoefMap(int nfld, int rows, int cols) + { + coefs.reset(); + coefs = std::make_shared(store.data(), nfld, rows, cols); + } + + //! Set the underlying coefficient data + void setCoefs(const Eigen::Tensor, 3>& tensor) + { + int rows = mmax+1; + int cols = nmax; + + const auto& d = tensor.dimensions(); + if (d[0] != nfld || d[1] != rows || d[2] != cols) + throw std::invalid_argument("setCoefs: dimension mismatch"); + + getCoefTensorMapMutable() = tensor; + } + + //! Direct access to underlying Eigen TensorMap (for client code) + const coefType& getCoefs() const { return getCoefTensorMap(); } + coefType& getCoefs() { return getCoefTensorMapMutable(); } }; @@ -558,4 +907,3 @@ namespace CoefClasses #endif - diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc index da76c6f28..4613a80c9 100644 --- a/expui/Coefficients.cc +++ b/expui/Coefficients.cc @@ -325,7 +325,7 @@ namespace CoefClasses coef->setGravConstant(G); coef->allocate(); - *coef->coefs = in; + coef->setCoefs(in); coefs[roundTime(Time)] = coef; } @@ -705,8 +705,7 @@ namespace CoefClasses throw std::runtime_error(str.str()); } else { it->second->store = dat; - it->second->coefs = std::make_shared - (it->second->store.data(), (Lmax+1)*(Lmax+2)/2, Nmax); + it->second->initCoefMap((Lmax+1)*(Lmax+2)/2, Nmax); } } @@ -720,7 +719,7 @@ namespace CoefClasses throw std::runtime_error(str.str()); } else { it->second->allocate(); - *it->second->coefs = dat; + it->second->setCoefs(dat); } } @@ -735,7 +734,7 @@ namespace CoefClasses ret.resize((Lmax+1)*(Lmax+2)/2, Nmax, ntim); for (int t=0; tcoefs); + auto & cof = coefs[roundTime(times[t])]->getCoefs(); for (int l=0; l<(Lmax+2)*(Lmax+1)/2; l++) { for (int n=0; nlmax, nmax= C->nmax; - Eigen::MatrixXcd out(*C->coefs); + Eigen::MatrixXcd out = C->getCoefs(); HighFive::DataSet dataset = stanza.createDataSet("coefficients", out); } @@ -959,7 +958,7 @@ namespace CoefClasses unsigned I = 0; if (lmin>0) I += lmin*lmin; - auto & cof = *(c.second->coefs); + auto & cof = c.second->getCoefs(); for (int ll=lmin; ll<=std::min(lmax, Lmax); ll++) { for (int mm=0; mm<=ll; mm++) { @@ -1012,8 +1011,8 @@ namespace CoefClasses << std::endl; for (auto v : coefs) { auto it = other->coefs.find(v.first); - auto & cv = *(v.second->coefs); - auto & ci = *(it->second->coefs); + auto & cv = v.second->getCoefs(); + auto & ci = it->second->getCoefs(); for (int i=0; icoefs->row(L); + auto rad = v.second->getCoefs().row(L); for (int n=std::max(0, min); n(nmax, max); n++) { double val = std::abs(rad(n)); power(T, l) += val * val; @@ -1212,8 +1211,7 @@ namespace CoefClasses throw std::runtime_error(str.str()); } else { it->second->store = dat; - it->second->coefs = std::make_shared - (it->second->store.data(), Mmax+1, Nmax); + it->second->initCoefMap(Mmax+1, Nmax); } } @@ -1227,7 +1225,7 @@ namespace CoefClasses throw std::runtime_error(str.str()); } else { it->second->allocate(); - *it->second->coefs = dat; + it->second->setCoefs(dat); } } @@ -1243,7 +1241,7 @@ namespace CoefClasses ret.resize(Mmax+1, Nmax, ntim); for (int t=0; tcoefs; + auto & cof = coefs[roundTime(times[t])]->getCoefs(); for (int m=0; mcoefs); + Eigen::MatrixXcd out = C->getCoefs(); HighFive::DataSet dataset = stanza.createDataSet("coefficients", out); } @@ -1419,7 +1417,7 @@ namespace CoefClasses for (auto c : coefs) { // The coefficient matrix - auto & cof = *c.second->coefs; + auto & cof = c.second->getCoefs(); // Index loop for (int mm=mmin; mm<=std::min(mmax, c.second->mmax); mm++) { @@ -1453,7 +1451,7 @@ namespace CoefClasses int T=0; for (auto v : coefs) { for (int m=0; m<=mmax; m++) { - auto rad = (*v.second->coefs).row(m); + auto rad = v.second->getCoefs().row(m); for (int n=std::max(0, min); n(nmax, max); n++) { double val = std::abs(rad(n)); power(T, m) += val * val; @@ -1515,7 +1513,7 @@ namespace CoefClasses int T=0; for (auto v : coefs) { for (int m=0; m<=mmax; m++) { - auto rad = (*v.second->coefs).row(m); + auto rad = v.second->getCoefs().row(m); // Even for (int n=std::max(0, min); n(nmax-nodd, max); n++) { double val = std::abs(rad(n)); @@ -1621,8 +1619,7 @@ namespace CoefClasses throw std::runtime_error(str.str()); } else { it->second->store = dat; - it->second->coefs = std::make_shared - (it->second->store.data(), 2*NmaxX+1, 2*NmaxY+1, NmaxZ); + it->second->initCoefMap(2*NmaxX+1, 2*NmaxY+1, NmaxZ); } } @@ -1635,8 +1632,8 @@ namespace CoefClasses str << "SlabCoefs::setTensor: requested time=" << time << " not found"; throw std::runtime_error(str.str()); } else { - it->second->allocate(); // Assign storage for the flattened tensor - *it->second->coefs = dat; // Populate using the tensor map + it->second->allocate(); // Assign storage for the flattened tensor + it->second->setCoefs(dat); // Populate using the tensor map } } @@ -1670,7 +1667,7 @@ namespace CoefClasses for (int ix=0; ix<=2*NmaxX; ix++) { for (int iy=0; iy<=2*NmaxY; iy++) { for (int iz=0; izcoefs)(ix, iy, iz); + ret(ix, iy, iz, t) = (cof->getCoefs())(ix, iy, iz); } } } @@ -1798,8 +1795,8 @@ namespace CoefClasses << std::setw(5) << ix << std::setw(5) << iy << std::setw(5) << iz - << std::setw(18) << (*c.second->coefs)(ix, iy, iz).real() - << std::setw(18) << (*c.second->coefs)(ix, iy, iz).imag() + << std::setw(18) << c.second->getCoefs()(ix, iy, iz).real() + << std::setw(18) << c.second->getCoefs()(ix, iy, iz).imag() << std::endl; } // Z loop @@ -1844,7 +1841,7 @@ namespace CoefClasses if (abs(iy - nmaxY) < min) continue; for (int iz=0; izcoefs)(ix, iy, iz)); + double val = std::abs(v.second->getCoefs()(ix, iy, iz)); power(T, ix) += val * val; } } @@ -1856,7 +1853,7 @@ namespace CoefClasses if (abs(ix - nmaxX) < min) continue; for (int iz=0; izcoefs)(ix, iy, iz)); + double val = std::abs(v.second->getCoefs()(ix, iy, iz)); power(T, iy) += val * val; } } @@ -1868,7 +1865,7 @@ namespace CoefClasses if (abs(ix - nmaxX) < min) continue; for (int iy=0; iy<=2*NmaxY; iy++) { if (abs(iy - nmaxY) < min) continue; - double val = std::abs((*v.second->coefs)(ix, iy, iz)); + double val = std::abs(v.second->getCoefs()(ix, iy, iz)); power(T, iz) += val * val; } } @@ -1926,8 +1923,8 @@ namespace CoefClasses << std::endl; for (auto v : coefs) { auto it = other->coefs.find(v.first); - auto & cv = *(v.second->coefs); - auto & ci = *(it->second->coefs); + auto & cv = v.second->getCoefs(); + auto & ci = it->second->getCoefs(); auto dim = cv.dimensions(); // This is an Eigen::Tensor map for (int i=0; isecond->store = dat; - it->second->coefs = std::make_shared - (it->second->store.data(), 2*NmaxX+1, 2*NmaxY+1, 2*NmaxZ+1); + it->second->initCoefMap(2*NmaxX+1, 2*NmaxY+1, 2*NmaxZ+1); } } @@ -2045,8 +2041,8 @@ namespace CoefClasses str << "CubeCoefs::setTensor: requested time=" << time << " not found"; throw std::runtime_error(str.str()); } else { - it->second->allocate(); // Assign storage for the flattened tensor - *it->second->coefs = dat; // Populate using the tensor map + it->second->allocate(); // Assign storage for the flattened tensor + it->second->setCoefs(dat); // Populate using the tensor map } } @@ -2080,10 +2076,10 @@ namespace CoefClasses for (int ix=0; ix<=2*NmaxX; ix++) { for (int iy=0; iy<=2*NmaxY; iy++) { for (int iz=0; iz<=2*NmaxZ; iz++) { - if (std::isnan(std::abs((*cof->coefs)(ix, iy, iz)))) + if (std::isnan(std::abs(cof->getCoefs()(ix, iy, iz)))) ret(ix, iy, iz, t) = 0.0; else - ret(ix, iy, iz, t) = (*cof->coefs)(ix, iy, iz); + ret(ix, iy, iz, t) = cof->getCoefs()(ix, iy, iz); } } } @@ -2211,8 +2207,8 @@ namespace CoefClasses << std::setw(5) << ix << std::setw(5) << iy << std::setw(5) << iz - << std::setw(18) << (*c.second->coefs)(ix, iy, iz).real() - << std::setw(18) << (*c.second->coefs)(ix, iy, iz).imag() + << std::setw(18) << c.second->getCoefs()(ix, iy, iz).real() + << std::setw(18) << c.second->getCoefs()(ix, iy, iz).imag() << std::endl; } // Z loop @@ -2252,7 +2248,7 @@ namespace CoefClasses if (abs(iy - nmaxY) < min) continue; for (int iz=0; iz<=2*NmaxZ; iz++) { if (abs(iz - nmaxZ) < min) continue; - double val = std::abs((*v.second->coefs)(ix, iy, iz)); + double val = std::abs(v.second->getCoefs()(ix, iy, iz)); power(T, ix) += val * val; } } @@ -2264,7 +2260,7 @@ namespace CoefClasses if (abs(ix - nmaxX) < min) continue; for (int iz=0; iz<=2*NmaxZ; iz++) { if (abs(iz - nmaxZ) < min) continue; - double val = std::abs((*v.second->coefs)(ix, iy, iz)); + double val = std::abs(v.second->getCoefs()(ix, iy, iz)); power(T, iy) += val * val; } } @@ -2276,7 +2272,7 @@ namespace CoefClasses if (abs(ix - nmaxX) < min) continue; for (int iy=0; iy<=2*NmaxY; iy++) { if (abs(iy - nmaxY) < min) continue; - double val = std::abs((*v.second->coefs)(ix, iy, iz)); + double val = std::abs(v.second->getCoefs()(ix, iy, iz)); power(T, iz) += val * val; } } @@ -2334,8 +2330,8 @@ namespace CoefClasses << std::endl; for (auto v : coefs) { auto it = other->coefs.find(v.first); - auto & cv = *(v.second->coefs); - auto & ci = *(it->second->coefs); + auto & cv = v.second->getCoefs(); + auto & ci = it->second->getCoefs(); auto dim = cv.dimensions(); // This is an Eigen::Tensor map for (int i=0; i& Times, const std::vector& data, bool verbose) : @@ -2379,16 +2374,13 @@ namespace CoefClasses c->traj = traj; c->rank = rank; c->store.resize(c->traj*c->rank); - c->coefs = std::make_shared(c->store.data(), c->traj, c->rank); - for (int m=0; mcoefs)(m, n) = data[m](i, n); - } - } - coefs[roundTime(c->time)] = c; + c->initCoefMap(traj, rank); + c->setCoefs(data[i]); + + coefs[roundTime(times[i])] = c; } } - + TrajectoryData::TrajectoryData(std::string& file, bool verbose) : Coefs("trajectory", verbose) { @@ -2440,14 +2432,13 @@ namespace CoefClasses coef->rank = rank; coef->time = times[n]; coef->store.resize(traj*rank); - coef->coefs = std::make_shared - (coef->store.data(), traj, rank); - *coef->coefs = data[n]; + coef->initCoefMap(traj, rank); + coef->setCoefs(data[n]); coefs[roundTime(times[n])] = coef; } } - + Eigen::VectorXcd& TrajectoryData::getData(double time) { auto it = coefs.find(roundTime(time)); @@ -2475,8 +2466,7 @@ namespace CoefClasses throw std::runtime_error(str.str()); } else { it->second->store = dat; - it->second->coefs = std::make_shared - (it->second->store.data(), it->second->traj, it->second->rank); + it->second->initCoefMap(it->second->traj, it->second->rank); } } @@ -2610,7 +2600,7 @@ namespace CoefClasses auto it = other->coefs.find(v.first); for (int m=0; mtraj; m++) { for (int n=0; nrank; n++) { - if ((*v.second->coefs)(m, n) != (*it->second->coefs)(m, n)) { + if (v.second->getCoefs()(m, n) != it->second->getCoefs()(m, n)) { ret = false; } } @@ -2640,7 +2630,7 @@ namespace CoefClasses auto cof = coefs[roundTime(times[t])]; for (int m=0; mcoefs)(m, n).real(); + ret(n, m, t) = cof->getCoefs()(m, n).real(); } } } @@ -2670,8 +2660,8 @@ namespace CoefClasses c->time = times[i]; c->cols = data[i].size(); c->store.resize(c->cols); - c->coefs = std::make_shared(c->store.data(), c->cols); - for (int j=0; jcols; j++) (*c->coefs)(j) = data[i][j]; + c->initCoefMap(); + for (int j=0; jcols; j++) c->getCoefs()(j) = data[i][j]; coefs[roundTime(c->time)] = c; } } @@ -2724,8 +2714,8 @@ namespace CoefClasses coef->cols = cols; coef->time = times[n]; coef->store.resize(cols); - coef->coefs = std::make_shared(coef->store.data(), cols); - for (int i=0; icoefs)(0, i) = data[n][i]; + coef->initCoefMap(); + for (int i=0; igetCoefs()(i) = data[n][i]; coefs[roundTime(times[n])] = coef; } @@ -2758,12 +2748,10 @@ namespace CoefClasses throw std::runtime_error(str.str()); } else { it->second->store = dat; - it->second->coefs = std::make_shared - (it->second->store.data(), dat.size()); - + it->second->initCoefMap(); } } - + void TableData::readNativeCoefs(const std::string& file, int stride, double tmin, double tmax) { @@ -2875,10 +2863,11 @@ namespace CoefClasses for (auto v : coefs) { auto it = other->coefs.find(v.first); for (int m=0; mcols; m++) { - if ((*v.second->coefs)(m) != (*it->second->coefs)(m)) { + if (v.second->getCoefs()(m) != it->second->getCoefs()(m)) { ret = false; } } + it->second->initCoefMap(); } } @@ -2900,7 +2889,7 @@ namespace CoefClasses for (int t=0; tcoefs)(c).real(); + ret(c, t) = cof->getCoefs()(c).real(); } } @@ -3086,7 +3075,7 @@ namespace CoefClasses auto it = other->coefs.find(v.first); for (int m=0; mmmax; m++) { for (int n=0; nnmax; n++) { - if ((*v.second->coefs)(m, n) != (*it->second->coefs)(m, n)) { + if (v.second->getCoefs()(m, n) != it->second->getCoefs()(m, n)) { ret = false; } } @@ -3320,8 +3309,8 @@ namespace CoefClasses } return *mat; - } - + } + void SphFldCoefs::setData(double time, Eigen::VectorXcd& dat) { auto it = coefs.find(roundTime(time)); @@ -3332,8 +3321,7 @@ namespace CoefClasses throw std::runtime_error(str.str()); } else { it->second->store = dat; - it->second->coefs = std::make_shared - (it->second->store.data(), Nfld, (Lmax+1)*(Lmax+2)/2, Nmax); + it->second->initCoefMap(Nfld, (Lmax+1)*(Lmax+2)/2, Nmax); } } @@ -3347,7 +3335,7 @@ namespace CoefClasses throw std::runtime_error(str.str()); } else { it->second->allocate(); - *it->second->coefs = dat; + it->second->setCoefs(dat); } } @@ -3362,7 +3350,7 @@ namespace CoefClasses ret.resize(Nfld, (Lmax+1)*(Lmax+2)/2, Nmax, ntim); for (int t=0; tcoefs); + auto & cof = coefs[roundTime(times[t])]->getCoefs(); for (int i=0; icoefs->dimensions(); + const auto& d = C->getCoefs().dimensions(); std::array shape {d[0], d[1], d[2]}; stanza.createAttribute("shape", HighFive::DataSpace::From(shape)).write(shape); @@ -3597,8 +3585,8 @@ namespace CoefClasses << std::endl; for (auto v : coefs) { auto it = other->coefs.find(v.first); - auto & cv = *(v.second->coefs); - auto & ci = *(it->second->coefs); + auto & cv = v.second->getCoefs(); + auto & ci = it->second->getCoefs(); const auto & d = cv.dimensions(); for (int n=0; ncoefs->dimensions(); + const auto& d = C->getCoefs().dimensions(); std::array shape {d[0], d[1], d[2]}; stanza.createAttribute("shape", HighFive::DataSpace::From(shape)).write(shape); @@ -3780,8 +3768,8 @@ namespace CoefClasses << std::endl; for (auto v : coefs) { auto it = other->coefs.find(v.first); - auto & cv = *(v.second->coefs); - auto & ci = *(it->second->coefs); + auto & cv = v.second->getCoefs(); + auto & ci = it->second->getCoefs(); const auto & d = cv.dimensions(); for (int n=0; n(0, min); n(nmax, max); n++) { - power(T, l) += std::norm((*v.second->coefs)(i, l, n)); + power(T, l) += std::norm(v.second->getCoefs()(i, l, n)); } } } @@ -3847,7 +3835,7 @@ namespace CoefClasses for (int i=1; i(0, min); n(nmax, max); n++) { - power(T, m) += std::norm((*v.second->coefs)(i, m, n)); + power(T, m) += std::norm(v.second->getCoefs()(i, m, n)); } } } @@ -3899,8 +3887,7 @@ namespace CoefClasses throw std::runtime_error(str.str()); } else { it->second->store = dat; - it->second->coefs = std::make_shared - (it->second->store.data(), Nfld, Mmax+1, Nmax); + it->second->initCoefMap(Nfld, Mmax+1, Nmax); } } @@ -3914,7 +3901,7 @@ namespace CoefClasses throw std::runtime_error(str.str()); } else { it->second->allocate(); - *it->second->coefs = dat; + it->second->setCoefs(dat); } } @@ -3929,7 +3916,7 @@ namespace CoefClasses ret.resize(Nfld, Mmax+1, Nmax, ntim); for (int t=0; tcoefs); + auto & cof = coefs[roundTime(times[t])]->getCoefs(); for (int i=0; i store; using dataType = Eigen::Tensor, 3>; using coefType = Eigen::TensorMap; - std::vector> coefs; + std::vector> coefs; void allocateStore(); //@} diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc index 7f193fc03..dec17b60b 100644 --- a/expui/FieldBasis.cc +++ b/expui/FieldBasis.cc @@ -199,10 +199,10 @@ namespace BasisClasses for (int t=0; t(store[t].data(), nfld, lmax+1, nmax); + coefs[t] = std::make_unique(store[t].data(), nfld, lmax+1, nmax); } else { store[t].resize(nfld*(lmax+1)*(lmax+2)/2*nmax); - coefs[t] = std::make_shared(store[t].data(), nfld, (lmax+1)*(lmax+2)/2, nmax); + coefs[t] = std::make_unique(store[t].data(), nfld, (lmax+1)*(lmax+2)/2, nmax); } } } @@ -285,8 +285,16 @@ namespace BasisClasses { if (dof==2) { auto p = dynamic_cast(c.get()); - coefs[0] = p->coefs; - store[0] = p->store; + + auto& src = p->getCoefs(); + + // Check that the size of the coefficient data matches the expected size + if (src.size() != store[0].size()) + throw std::runtime_error( + "FieldBasis::set_coefs: CylFldStruct data size mismatch with FieldBasis store"); + + // Copy the data from the input structure to the internal store + std::copy(src.data(), src.data() + src.size(), store[0].data()); // Sanity test dimensions if (nfld!=p->nfld || lmax!=p->mmax || nmax!=p->nmax) { @@ -300,8 +308,15 @@ namespace BasisClasses } else { auto p = dynamic_cast(c.get()); - coefs[0] = p->coefs; - store[0] = p->store; + + // Check that the size of the coefficient data matches the expected size + if (p->getCoefs().size() == store[0].size()) + throw std::runtime_error + ("FieldBasis::set_coefs: SphFldStruct data size mismatch with FieldBasis store"); + + // Copy the data from the input structure to the internal store + std::copy(store[0].data(), store[0].data() + store[0].size(), + p->getCoefs().data()); // Sanity test dimensions if (nfld!=p->nfld || lmax!=p->lmax || nmax!=p->nmax) { diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 2bc068a51..5030c42df 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -2412,6 +2412,11 @@ void EmpCylSL::generate_eof(int numr, int nump, int numt, orthoTest(ortho->orthoCheck(std::max(NMAX*50, 200)), "EmpCylSL[SLGridSph]", "l"); } + // Save quadrature parameters for provenance + // + eof_vars.numr = numr; + eof_vars.nump = nump; + eof_vars.numt = numt; // Initialize fixed variables and storage // @@ -7413,29 +7418,51 @@ void EmpCylSL::WriteH5Cache() std::string forceID("Cylinder"), geometry("cylinder"); std::string model = EmpModelLabs[mtype]; - file.createAttribute("geometry", HighFive::DataSpace::From(geometry)).write(geometry); - file.createAttribute("forceID", HighFive::DataSpace::From(forceID)).write(forceID); - file.createAttribute("Version", HighFive::DataSpace::From(Version)).write(Version); + file.createAttribute("geometry", HighFive::DataSpace::From(geometry)).write(geometry); + file.createAttribute("forceID", HighFive::DataSpace::From(forceID)). write(forceID); + file.createAttribute("Version", HighFive::DataSpace::From(Version)). write(Version); - // Write the specific parameters + // Write the grid specific parameters // - file.createAttribute("model", HighFive::DataSpace::From(model)). write(model); - file.createAttribute ("mmax", HighFive::DataSpace::From(MMAX)). write(MMAX); - file.createAttribute ("numx", HighFive::DataSpace::From(NUMX)). write(NUMX); - file.createAttribute ("numy", HighFive::DataSpace::From(NUMY)). write(NUMY); - file.createAttribute ("nmax", HighFive::DataSpace::From(NORDER)). write(NORDER); - file.createAttribute ("lmaxfid", HighFive::DataSpace::From(LMAX)). write(LMAX); - file.createAttribute ("nmaxfid", HighFive::DataSpace::From(NMAX)). write(NMAX); - file.createAttribute ("neven", HighFive::DataSpace::From(Neven)). write(Neven); - file.createAttribute ("nodd", HighFive::DataSpace::From(Nodd)). write(Nodd); - file.createAttribute ("cmapr", HighFive::DataSpace::From(CMAPR)). write(CMAPR); - file.createAttribute ("cmapz", HighFive::DataSpace::From(CMAPZ)). write(CMAPZ); - file.createAttribute ("rmin", HighFive::DataSpace::From(RMIN)). write(RMIN); - file.createAttribute ("rmax", HighFive::DataSpace::From(RMAX)). write(RMAX); - file.createAttribute ("ascl", HighFive::DataSpace::From(ASCALE)). write(ASCALE); - file.createAttribute ("hscl", HighFive::DataSpace::From(HSCALE)). write(HSCALE); - file.createAttribute ("cmass", HighFive::DataSpace::From(cylmass)). write(cylmass); - + file.createAttribute("model", HighFive::DataSpace::From(model)). write(model); + file.createAttribute ("mmax", HighFive::DataSpace::From(MMAX)). write(MMAX); + file.createAttribute ("numx", HighFive::DataSpace::From(NUMX)). write(NUMX); + file.createAttribute ("numy", HighFive::DataSpace::From(NUMY)). write(NUMY); + file.createAttribute ("nmax", HighFive::DataSpace::From(NORDER)). write(NORDER); + file.createAttribute ("lmaxfid", HighFive::DataSpace::From(LMAX)). write(LMAX); + file.createAttribute ("nmaxfid", HighFive::DataSpace::From(NMAX)). write(NMAX); + file.createAttribute ("neven", HighFive::DataSpace::From(Neven)). write(Neven); + file.createAttribute ("nodd", HighFive::DataSpace::From(Nodd)). write(Nodd); + file.createAttribute ("cmapr", HighFive::DataSpace::From(CMAPR)). write(CMAPR); + file.createAttribute ("cmapz", HighFive::DataSpace::From(CMAPZ)). write(CMAPZ); + file.createAttribute ("rmin", HighFive::DataSpace::From(RMIN)). write(RMIN); + file.createAttribute ("rmax", HighFive::DataSpace::From(RMAX)). write(RMAX); + file.createAttribute ("ascl", HighFive::DataSpace::From(ASCALE)). write(ASCALE); + file.createAttribute ("hscl", HighFive::DataSpace::From(HSCALE)). write(HSCALE); + file.createAttribute ("cmass", HighFive::DataSpace::From(cylmass)). write(cylmass); + + // Additional unchecked variables + // + file.createAttribute ("PCAVAR", HighFive::DataSpace::From(PCAVAR)). write(PCAVAR); + file.createAttribute ("PCAVTK", HighFive::DataSpace::From(PCAVTK)). write(PCAVTK); + file.createAttribute ("PCAEOF", HighFive::DataSpace::From(PCAEOF)). write(PCAEOF); + file.createAttribute ("numr", HighFive::DataSpace::From(NUMR)). write(NUMR); + file.createAttribute ("hexp", HighFive::DataSpace::From(HEXP)). write(HEXP); + file.createAttribute ("hfac", HighFive::DataSpace::From(HFAC)). write(HFAC); + file.createAttribute ("ppow", HighFive::DataSpace::From(PPOW)). write(PPOW); + + // Options + // + file.createAttribute ("logarithmic", HighFive::DataSpace::From(logarithmic)). write(logarithmic); + file.createAttribute ("enforce_limits", HighFive::DataSpace::From(enforce_limits)).write(enforce_limits); + + // EOF quadrature parameters + // + file.createAttribute ("eof_numr", HighFive::DataSpace::From(eof_vars.numr)).write(eof_vars.numr); + file.createAttribute ("eof_nump", HighFive::DataSpace::From(eof_vars.nump)).write(eof_vars.nump); + file.createAttribute ("eof_numt", HighFive::DataSpace::From(eof_vars.numt)).write(eof_vars.numt); + + // Cosine functions auto cosine = file.createGroup("Cosine"); diff --git a/include/EmpCylSL.H b/include/EmpCylSL.H index dff248f2b..3b92175fc 100644 --- a/include/EmpCylSL.H +++ b/include/EmpCylSL.H @@ -78,6 +78,11 @@ protected: double HSCALE; double pfac, dfac, ffac; + // For provenance +struct eof_variables { + int numr{-1}, nump{-1}, numt{-1}; +} eof_vars; + std::vector facC, facS; int rank2, rank3; diff --git a/include/ScopeTimer.H b/include/ScopeTimer.H new file mode 100644 index 000000000..fe638ad48 --- /dev/null +++ b/include/ScopeTimer.H @@ -0,0 +1,64 @@ +#pragma once + +#include +#include +#include + +class ScopeTimer { +public: + using Clock = std::chrono::steady_clock; + + // Define units for output + enum class Unit { + Seconds, + Milliseconds, + Microseconds, + Nanoseconds + }; + + // Default to Seconds if the user doesn't choose + explicit ScopeTimer(std::string task_name, Unit output_unit = Unit::Seconds) + : name(std::move(task_name)), unit(output_unit), start_time(Clock::now()) {} + + // Output on deletion of the object + ~ScopeTimer() + { + const auto end_time = Clock::now(); + const auto delta = end_time - start_time; + + // Convert natively using std::chrono based on user choice + switch (unit) { + case Unit::Seconds: { + std::chrono::duration elapsed = delta; + std::cout << "[" << name << "] completed in " << elapsed.count() << " s\n"; + break; + } + case Unit::Milliseconds: { + std::chrono::duration elapsed = delta; + std::cout << "[" << name << "] completed in " << elapsed.count() << " ms\n"; + break; + } + case Unit::Microseconds: { + std::chrono::duration elapsed = delta; + std::cout << "[" << name << "] completed in " << elapsed.count() << " µs\n"; + break; + } + case Unit::Nanoseconds: { + std::chrono::duration elapsed = delta; + std::cout << "[" << name << "] completed in " << elapsed.count() << " ns\n"; + break; + } + default: + std::cerr << "Unknown time unit for ScopeTimer\n"; + break; + } + } + + ScopeTimer(const ScopeTimer&) = delete; + ScopeTimer& operator=(const ScopeTimer&) = delete; + +private: + std::string name; + Unit unit; + std::chrono::time_point start_time; +}; diff --git a/src/Component.H b/src/Component.H index 923796ca1..48dc2c257 100644 --- a/src/Component.H +++ b/src/Component.H @@ -184,6 +184,43 @@ struct loadb_datum specifying alternative parameters, using an old-style PSP file, or starting a new job with a previous output from another simulation) + +
Input data format Component now supports reading from + HDF5 files on input. The scheme for HDF5 input files is: + @verbatim + / + ├── Attributes + │ ├── num_particles (int) + │ ├── num_aux_ints (int) + │ └── num_aux_floats (int) + └── particles/ (Group) + ├── m (Dataset) + ├── x, y, z (Datasets) + ├── u, v, w (Datasets) + ├── index (optional Dataset) + ├── aux_int_0, aux_int_1, ... (Datasets) + └── aux_float_0, aux_float_1, ... (Datasets) + @endverbatim + + @note +
    + +
  1. The index dataset is optional. If it is not + present, the particles are assigned indices in the order they are read. + +
  2. The aux_int_* and aux_float_* datasets are optional. + If omitted, the corresponding attribute counts (num_aux_ints, num_aux_floats) + should be set to 0; if the counts are non-zero, datasets aux_int_0.. and + aux_float_0.. are expected to exist. + +
  3. The datatype for the m, x, + y, z, u, v, + w and aux_float_* datasets may be float or + double (i.e. float32 or float64). Each dataset's precision is + detected independently, so datasets may have mixed precision within + the same file. + +
*/ class Component { @@ -199,6 +236,8 @@ private: //! Parallel distribute and particle io void load_balance(void); void update_indices(void); + void read_bodies_and_distribute_init(void); + void read_bodies_and_distribute_hdf5(void); void read_bodies_and_distribute_ascii(void); void read_bodies_and_distribute_binary_out(istream *); void read_bodies_and_distribute_binary_spl(istream *); @@ -367,6 +406,9 @@ protected: int H5chunk = 4096; //@} + //! Detect float precision by inspecting dataset type + size_t detect_precision(const HighFive::DataSet& dataset); + public: //! Describe the phase space coordinates diff --git a/src/Component.cc b/src/Component.cc index c50884ce6..526f65a01 100644 --- a/src/Component.cc +++ b/src/Component.cc @@ -2,11 +2,20 @@ #include #include #include +#include #include #include #include #include +// HighFive +#include +#include +#include + +// HDF5 C API +#include + #include "Component.H" #include "Bessel.H" #include "Sphere.H" @@ -277,7 +286,7 @@ Component::Component(YAML::Node& CONF) initialize_cuda(); - read_bodies_and_distribute_ascii(); + read_bodies_and_distribute_init(); reset_level_lists(); @@ -1463,6 +1472,372 @@ Component::~Component(void) delete [] acc0; } +void Component::read_bodies_and_distribute_init() +{ + // Check whether input file is HDF5 + + // Note: One might need to use H5Fis_hdf5(filename.c_str()) for + // older HDF5 versions (< 1.12) + + bool isHDF5 = false; + + if (myid == 0) { +#if H5_VERSION_GE(1, 12, 0) + isHDF5 = H5Fis_accessible(pfile.c_str(), H5P_DEFAULT) > 0; +#else + isHDF5 = H5Fis_hdf5(pfile.c_str()) > 0; +#endif + } + + MPI_Bcast(&isHDF5, 1, MPI_CXX_BOOL, 0, MPI_COMM_WORLD); + + if (isHDF5) { + if (myid==0) + std::cout << "---- Component <" << name << ">: " + << "reading bodies from HDF5 file "<< pfile << std::endl; + read_bodies_and_distribute_hdf5(); + } else { + if (myid==0) + std::cout << "---- Component <" << name << ">: " + << "reading bodies from ASCII file "<< pfile << std::endl; + read_bodies_and_distribute_ascii(); + } +} + +// Detect float precision by inspecting dataset type +size_t Component::detect_precision(const HighFive::DataSet& dataset) +{ + const auto datatype = dataset.getDataType(); + if (datatype.getClass() != HighFive::DataTypeClass::Float) { + throw std::runtime_error("Unexpected dataset type: expected float for precision detection"); + } + + // Get the size in bytes + const size_t type_size = datatype.getSize(); + + // Float32 is 4 bytes, Float64 is 8 bytes + if (type_size == 4 || type_size == 8) { + return type_size; + } + + throw std::runtime_error("Unexpected float type size: " + std::to_string(type_size)); +} + +void Component::read_bodies_and_distribute_hdf5(void) +{ + HighFive::File file(pfile, HighFive::File::ReadOnly); + + // Read metadata (all ranks read the same metadata, but only rank 0 + // reads the datasets) + // + nbodies_tot = file.getAttribute("num_particles") .read(); + niattrib = file.getAttribute("num_aux_ints") .read(); + ndattrib = file.getAttribute("num_aux_floats").read(); + + // Open particles group + // + HighFive::Group particles_group = file.getGroup("particles"); + + // Variant for flexible reading and access of float or double datasets + // + using FloatDataVariant = std::variant, std::vector>; + + // Wrapper struct to provide uniform access to float or double std::vectors + // + struct FloatData + { + // The underlying data can be either a vector of floats or doubles + FloatDataVariant data; + + // Default constructor + FloatData() = default; + + // Implicit constructor so you can assign vectors directly to FloatData + FloatData(FloatDataVariant v) : data(std::move(v)) {} + + // Provides a uniform interface to access elements regardless of + // underlying type. Forces a value copy, converting float to + // double if necessary. + double operator[](size_t i) const { + return std::visit([i](const auto& vec) -> double { + return vec[i]; }, data); + } + }; + + // Phase-space data + // + FloatData m, x, y, z, u, v, w; + std::vector index; + bool has_index = particles_group.exist("index"); + std::vector> aux_ints; + std::vector aux_floats; + + // Lambda to read float or double based on precision flag + // + auto read_dataset = [this, &particles_group] + (const std::string& name, size_t offset, size_t batch) -> FloatData { + + HighFive::DataSet dataset = particles_group.getDataSet(name); + size_t precision = detect_precision(dataset); + + // Fetch the dimensions (e.g., [rows, cols]) + size_t total = dataset.getSpace().getDimensions()[0]; + + if (offset >= total) { + throw std::runtime_error("Component ERROR: dataset <" + name + "> shorter than expected"); + } + + size_t current = std::min(batch, total - offset); + if (current != batch) { + throw std::runtime_error("Component ERROR: dataset <" + name + "> shorter than expected"); + } + + if (precision == 4) { + return FloatData(dataset.select({offset}, {current}).read>()); + } + else if (precision == 8) { + return FloatData(dataset.select({offset}, {current}).read>()); + } + else { + throw std::runtime_error("Unsupported precision detected in dataset: " + name); + } + }; + + // Okay, now we can read the datasets and distribute particles + // across MPI rank, following the previously established logic. The + // following code reads the datasets in batches and distributes them + // to the appropriate ranks. + + // For radius exclusion check + double rmax1=0.0, r2; + + is_init = 1; + setup_distribution(); + is_init = 0; + + // Initialize the particle ferry instance with dynamic attribute + // sizes + if (not pf) pf = ParticleFerryPtr(new ParticleFerry(niattrib, ndattrib)); + + + if (myid==0) { + + // Offset into HDF5 datasets for each node, read in batchs for all + // datasets + size_t offset = 0; + + // Read Node 0's PSP fields + if (has_index) { + index = particles_group.getDataSet("index") + .select({offset}, {nbodies_table[0]}).read>(); + } + m = read_dataset("m", offset, nbodies_table[0]); + + x = read_dataset("x", offset, nbodies_table[0]); + y = read_dataset("y", offset, nbodies_table[0]); + z = read_dataset("z", offset, nbodies_table[0]); + + u = read_dataset("u", offset, nbodies_table[0]); + v = read_dataset("v", offset, nbodies_table[0]); + w = read_dataset("w", offset, nbodies_table[0]); + + // Read auxiliary integer fields + aux_ints.resize(niattrib); + for (int j = 0; j < niattrib; ++j) { + std::string dset_name = "aux_int_" + std::to_string(j); + aux_ints[j] = particles_group.getDataSet(dset_name) + .select({offset}, {nbodies_table[0]}).read>(); + } + + // Read auxiliary float fields + aux_floats.resize(ndattrib); + for (int j = 0; j < ndattrib; ++j) { + std::string dset_name = "aux_float_" + std::to_string(j); + aux_floats[j] = read_dataset(dset_name, offset, nbodies_table[0]); + } + + for (unsigned i=0; i(niattrib, ndattrib); + + // Load the particle data for root + part->mass = m[i]; + + part->pos[0] = x[i]; + part->pos[1] = y[i]; + part->pos[2] = z[i]; + + part->vel[0] = u[i]; + part->vel[1] = v[i]; + part->vel[2] = w[i]; + + for (int j=0; jiattrib[j] = aux_ints[j][i]; + } + + for (int j=0; jdattrib[j] = aux_floats[j][i]; + } + + // Set particle index from optional dataset, otherwise sequential + part->indx = has_index ? index[i] : i + 1; + + // Get the radius + r2 = 0.0; + for (int k=0; k<3; k++) r2 += part->pos[k]*part->pos[k]; + rmax1 = max(r2, rmax1); + + // Load the particle + particles[part->indx] = part; + + // Record top_seq + top_seq = std::max(part->indx, top_seq); + } + + nbodies = nbodies_table[0]; + + // Offset into the next node's data in the HDF5 datasets + offset += nbodies_table[0]; + + for (int n=1; nShipParticles(n, 0, zero); + continue; + } + + // Read core PSP fields + if (has_index) { + index = particles_group.getDataSet("index") + .select({offset}, {nbodies_table[n]}).read>(); + } + m = read_dataset("m", offset, nbodies_table[n]); + x = read_dataset("x", offset, nbodies_table[n]); + y = read_dataset("y", offset, nbodies_table[n]); + z = read_dataset("z", offset, nbodies_table[n]); + + u = read_dataset("u", offset, nbodies_table[n]); + v = read_dataset("v", offset, nbodies_table[n]); + w = read_dataset("w", offset, nbodies_table[n]); + + // Read auxiliary integer fields + aux_ints.resize(niattrib); + for (int j = 0; j < niattrib; ++j) { + std::string dset_name = "aux_int_" + std::to_string(j); + aux_ints[j] = particles_group.getDataSet(dset_name) + .select({offset}, {nbodies_table[n]}).read>(); + } + + // Read auxiliary float fields + aux_floats.resize(ndattrib); + for (int j = 0; j < ndattrib; ++j) { + std::string dset_name = "aux_float_" + std::to_string(j); + aux_floats[j] = read_dataset(dset_name, offset, nbodies_table[n]); + } + + pf->ShipParticles(n, 0, nbodies_table[n]); + + unsigned icount = 0; + unsigned ibufcount = 0; + while (icount < nbodies_table[n]) { + + PartPtr part = std::make_shared(niattrib, ndattrib); + + // Set particle index from optional dataset, otherwise sequential + part->indx = has_index ? index[icount] : nbodies_index[n-1] + 1 + icount; + + part->mass = m[icount]; + + part->pos[0] = x[icount]; + part->pos[1] = y[icount]; + part->pos[2] = z[icount]; + + part->vel[0] = u[icount]; + part->vel[1] = v[icount]; + part->vel[2] = w[icount]; + + for (int j=0; jiattrib[j] = aux_ints[j][icount]; + } + + for (int j=0; jdattrib[j] = aux_floats[j][icount]; + } + + r2 = 0.0; + for (int k=0; k<3; k++) r2 += part->pos[k]*part->pos[k]; + rmax1 = max(r2, rmax1); + + pf->SendParticle(part); + icount++; + + // Record top_seq + top_seq = std::max(part->indx, top_seq); + } + + // Update offset for next node + offset += nbodies_table[n]; + } + + } else { + + pf->ShipParticles(myid, 0, nbodies); + +#ifdef DEBUG + int icount = 0; +#endif + + while (PartPtr part=pf->RecvParticle()) { + particles[part->indx] = part; +#ifdef DEBUG + if (icount<5) { + cout << "Process " << myid << ": received "; + cout << setw(14) << part->mass; + for (int k=0; k<3; k++) cout << setw(14) << part->pos[k]; + cout << endl; + } + icount++; +#endif + } + } + + // Default: set to max radius + // can be overriden by parameter + + rmax = sqrt(fabs(rmax1)); + MPI_Bcast(&rmax, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + + // Send top_seq to all nodes + MPI_Bcast(&top_seq, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD); + + initialize(); + +#ifdef DEBUG + if (particles.size()) { + unsigned long imin = std::numeric_limits::max(); + unsigned long imax = 0, kmin = imax, kmax = 0; + for (auto p : particles) { + imin = std::min(imin, p.first); + imax = std::max(imax, p.first); + kmin = std::min(kmin, p.second->indx); + kmax = std::max(kmax, p.second->indx); + } + cout << "read_bodies_and_distribute_hdf5: process " << myid + << " name=" << name << " bodies [" << kmin << ", " + << kmax << "], [" << imin << ", " << imax << "]" + << " #=" << particles.size() << endl; + } else { + cout << "read_bodies_and_distribute_hdf5: process " << myid + << " name=" << name + << " #=" << particles.size() << endl; + } +#endif +} + + void Component::read_bodies_and_distribute_ascii(void) { // Open file diff --git a/src/Cube.cc b/src/Cube.cc index f4e184dca..eee48ca51 100644 --- a/src/Cube.cc +++ b/src/Cube.cc @@ -820,7 +820,7 @@ void Cube::dump_coefs_h5(const std::string& file) cur->allocate(); // Set the storage and copy the // coefficients through the map - *cur->coefs = expcoef; + cur->setCoefs(expcoef); // Check if file exists // diff --git a/src/Cylinder.cc b/src/Cylinder.cc index b8bd5033a..c53997ff3 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -1865,7 +1865,7 @@ void Cylinder::dump_coefs_h5(const std::string& file) cur->allocate(); - auto & cof = *cur->coefs; // Reference to the coefficient map + auto & cof = cur->getCoefs(); // Reference to the coefficient map Eigen::VectorXd cos1(nmax), sin1(nmax); diff --git a/src/PolarBasis.cc b/src/PolarBasis.cc index 1a34864a2..ba2ef6d91 100644 --- a/src/PolarBasis.cc +++ b/src/PolarBasis.cc @@ -1840,7 +1840,7 @@ void PolarBasis::dump_coefs_h5(const std::string& file) cur->allocate(); // Assign storage and make the map - auto & cof = *cur->coefs; // Reference for convenience + auto & cof = cur->getCoefs(); // Reference for convenience for (int ir=0; irallocate(); // Set the storage and copy the // coefficients through the map - *cur->coefs = expccof[0]; + cur->setCoefs(expccof[0]); // Check if file exists // diff --git a/src/SphericalBasis.cc b/src/SphericalBasis.cc index e69cba4af..b20cfab6e 100644 --- a/src/SphericalBasis.cc +++ b/src/SphericalBasis.cc @@ -1922,7 +1922,7 @@ void SphericalBasis::dump_coefs_h5(const std::string& file) cur->allocate(); - auto & cof = *cur->coefs; // Reference for convenience + auto & cof = cur->getCoefs(); // Reference for convenience for (int ir=0; ir +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// For OpenMP parallelization +#include + +// HighFive +#include +#include +#include + +// Required for H5Zfilter_avail +#include + +// Command-line parsing +#include "cxxopts.H" + +// For command timing +#include "ScopeTimer.H" + +using namespace HighFive; + + +// Factory designed to iterate across filter list available on Ubuntu +// 26.04 LTS and other systems with HDF5 1.14.x and HighFive 2.6.x +HighFive::DataSetCreateProps +createFilterProps(const std::vector& chunk_dims, + unsigned int filter_id, + // Scale from 1 (fast) to 9 (max compression) + int compression_level = 5) +{ + HighFive::DataSetCreateProps props; + props.add(HighFive::Chunking(chunk_dims)); + + // For most compressors, applying shuffle improves compression for floating-point data. + // (Blosc can do its own shuffling internally; shuffle-only/checksum-only are handled below.) + if (filter_id != 3 && filter_id != 4 && filter_id != 32001) { + props.add(HighFive::Shuffle()); + } + hid_t plist_id = props.getId(); + std::vector cd_values; + + switch (filter_id) { + case 1: // Deflate / GZIP (Built-in) + // Expects 1 parameter: compression level (0-9) + cd_values = { static_cast(compression_level) }; + break; + + case 2: // SZIP (Built-in) + // SZIP is complex and highly variant; HighFive includes native wrappers. + // For raw testing, standard pixels-per-block option is typically 16. + cd_values = { 141, 16 }; + break; + + case 3: // Shuffle Filter alone (No trailing compression) + props.add(HighFive::Shuffle()); + return props; + + case 4: // Fletcher32 Checksum alone (Data validation, not compression) + // Fletcher32 uses built-in ID 4 and accepts 0 configuration arguments. + // Call the native HDF5 library directly using HighFive's internal handle. + H5Pset_filter(plist_id, 4, H5Z_FLAG_OPTIONAL, 0, NULL); + return props; + + case 307: // BZip2 + // Expects 1 parameter: Block size in 100KB steps (1-9) + cd_values = { static_cast(compression_level) }; + break; + + case 32004: // LZ4 + // Expects 1 parameter: internal chunk/block size. + // 0 falls back to the default (64KB), optimized for CPU caches. + cd_values = { 0 }; + break; + + case 32001: // Blosc (v1 Meta-Compressor) + // Expects 7 parameters. Slots 0-3 are reserved. + // [4]=level(1-9), [5]=shuffle type(1=byte, 2=bit), [6]=codec code + // Codecs: 0=blosclz, 1=lz4, 2=lz4hc, 3=snappy, 4=zlib, 5=zstd + cd_values = { 0, 0, 0, 0, + static_cast(compression_level), + 1, // 1 = Byte Shuffle (Highly recommended for float/double) + 1 }; // 1 = Redirect internal processing to LZ4 + break; + + default: + throw std::invalid_argument("Unknown or unsupported testing filter ID requested."); + } + + // Pass the custom properties down to the HDF5 backend pipeline + if (!cd_values.empty()) { + H5Pset_filter(plist_id, filter_id, H5Z_FLAG_OPTIONAL, cd_values.size(), cd_values.data()); + } + + return props; +} + +// Specify floating point precision +enum class FloatPrecision { FLOAT32, FLOAT64 }; + +// Type variant for flexible float handling +template +struct ParticleDataTemplate +{ + std::optional> index; // optional particle index + std::vector m; // mass + std::vector x, y, z; // position + std::vector u, v, w; // velocity + std::vector> aux_ints; // auxiliary integer fields + std::vector> aux_floats; // auxiliary float fields + + size_t num_particles = 0; + size_t num_aux_ints = 0; + size_t num_aux_floats = 0; +}; + +// Helper: Parse a single line of particle data +template +void parse_particle_line(const std::string& line, + int particle_idx, + bool has_index, + int num_aux_ints, + int num_aux_floats, + ParticleDataTemplate& data) +{ + std::istringstream iss(line); + // Parse optional index field + if (has_index) { + unsigned long idx = 0; + if (!(iss >> idx)) { + throw std::runtime_error("Failed to parse index at particle " + + std::to_string(particle_idx)); + } + (*data.index)[particle_idx] = idx; + } + + // Parse core PSP fields + if (!(iss + >> data.m[particle_idx] + >> data.x[particle_idx] >> data.y[particle_idx] >> data.z[particle_idx] + >> data.u[particle_idx] >> data.v[particle_idx] >> data.w[particle_idx])) { + throw std::runtime_error("Failed to parse core fields at particle " + + std::to_string(particle_idx)); + } + + // Parse auxiliary integer fields + for (int j = 0; j < num_aux_ints; ++j) { + if (!(iss >> data.aux_ints[j][particle_idx])) { + throw std::runtime_error("Failed to parse aux int at particle " + + std::to_string(particle_idx)); + } + } + + // Parse auxiliary float fields + for (int j = 0; j < num_aux_floats; ++j) { + if (!(iss >> data.aux_floats[j][particle_idx])) { + throw std::runtime_error("Failed to parse aux float at particle " + + std::to_string(particle_idx)); + } + } +} + +// Read entire ASCII file into memory +std::vector read_ascii_lines(const std::string& ascii_file, + int& num_particles, + int& num_aux_ints, + int& num_aux_floats) +{ + // Open ASCII file + std::ifstream infile(ascii_file); + if (!infile.is_open()) { + throw std::runtime_error("Could not open ASCII file: " + ascii_file); + } + + // Read header + std::string header_line; + if (std::getline(infile, header_line)) { + std::istringstream header_stream(header_line); + if (!(header_stream >> num_particles >> num_aux_ints >> num_aux_floats)) { + throw std::runtime_error("Failed to parse header line: " + header_line); + } + } + else { + throw std::runtime_error("Failed to read header from ASCII file"); + } + + // Validate header + if (num_particles <= 0 || num_aux_ints < 0 || num_aux_floats < 0) { + throw std::runtime_error("Invalid header values"); + } + + // Read all lines + std::vector lines; + lines.reserve(num_particles); + std::string line; + while (std::getline(infile, line)) { + if (!line.empty() && line[0] != '#') { // Skip empty/comment lines + lines.push_back(line); + } + } + infile.close(); + + if ((int)lines.size() != num_particles) { + throw std::runtime_error("Mismatch between header count and actual particle count"); + } + + return lines; +} + +// Read ASCII and write HDF5 with specified precision +template +void ascii_to_hdf5_impl(const std::string& ascii_file, + const std::string& hdf5_file, + FloatPrecision precision, + bool has_index, + unsigned filter_id, + bool verbose) +{ + // Header values + int num_particles, num_aux_ints, num_aux_floats; + + // Read all ASCII lines into memory + if (verbose) std::cout << "Reading ASCII file..." << std::endl; + std::vector lines = read_ascii_lines(ascii_file, num_particles, + num_aux_ints, num_aux_floats); + + // Initialize particle data structure + ParticleDataTemplate data; + + data.num_particles = num_particles; + data.num_aux_ints = num_aux_ints; + data.num_aux_floats = num_aux_floats; + + data.m.resize(num_particles); + data.x.resize(num_particles); + data.y.resize(num_particles); + data.z.resize(num_particles); + data.u.resize(num_particles); + data.v.resize(num_particles); + data.w.resize(num_particles); + if (has_index) data.index = std::vector(num_particles); + + data.aux_ints.resize(num_aux_ints); + for (auto& vec : data.aux_ints) { + vec.resize(num_particles); + } + + data.aux_floats.resize(num_aux_floats); + for (auto& vec : data.aux_floats) { + vec.resize(num_particles); + } + + // Parse particle data in parallel + if (verbose) std::cout << "Parsing particles (parallel)..." << std::endl; + + int parse_errors = 0; + #pragma omp parallel for schedule(dynamic, 256) reduction(+:parse_errors) \ + num_threads(omp_get_max_threads()) + for (int i = 0; i < num_particles; ++i) { + try { + parse_particle_line(lines[i], i, has_index, num_aux_ints, num_aux_floats, data); + } catch (const std::exception& e) { + parse_errors += 1; + #pragma omp critical + { + std::cerr << "Thread error at particle " << i << ": " << e.what() << std::endl; + } + } + } + + if (parse_errors > 0) { + throw std::runtime_error("Encountered " + std::to_string(parse_errors) + + " particle parse errors; aborting conversion."); + } + // Create HDF5 file with compression enabled (SEQUENTIAL - thread-safe) + if (verbose) std::cout << "Writing HDF5 file..." << std::endl; + File file(hdf5_file, File::ReadWrite | File::Create | File::Truncate); + + // Store header information and precision metadata + file.createAttribute("num_particles", DataSpace::From(num_particles)) + .write(num_particles); + file.createAttribute("num_aux_ints", DataSpace::From(num_aux_ints)) + .write(num_aux_ints); + file.createAttribute("num_aux_floats", DataSpace::From(num_aux_floats)) + .write(num_aux_floats); + + // Create a group for particle data + // + Group particles_group = file.createGroup("particles"); + + // Calculate optimized, capped chunk size + // + hsize_t chunk_size = 262144; // Cachesafe cap (~2MB for double) + if (num_particles < chunk_size) { + chunk_size = std::max(1024, num_particles); + } + if (chunk_size > num_particles) { + chunk_size = num_particles; // Absolute safety fallback + } + + // Define compression filter + // + auto props = createFilterProps({chunk_size}, filter_id, 5); + + // Write datasets sequentially (HDF5 is not fully thread-safe for writing) + particles_group.createDataSet("m", data.m, props); + particles_group.createDataSet("x", data.x, props); + particles_group.createDataSet("y", data.y, props); + particles_group.createDataSet("z", data.z, props); + particles_group.createDataSet("u", data.u, props); + particles_group.createDataSet("v", data.v, props); + particles_group.createDataSet("w", data.w, props); + if (data.index.has_value()) { + particles_group.createDataSet("index", *data.index, props); + } + + // Write auxiliary integer fields + for (size_t j = 0; j < data.aux_ints.size(); ++j) { + std::string dset_name = "aux_int_" + std::to_string(j); + particles_group.createDataSet(dset_name, data.aux_ints[j], props); + } + + // Write auxiliary float fields + for (size_t j = 0; j < data.aux_floats.size(); ++j) { + std::string dset_name = "aux_float_" + std::to_string(j); + particles_group.createDataSet(dset_name, data.aux_floats[j], props); + } + + if (verbose) { + std::string precision_str = (precision == FloatPrecision::FLOAT32) ? "float32" : "float64"; + std::cout << "Successfully wrote " << num_particles << " particles to " + << hdf5_file << " (" << precision_str << ")" << std::endl; + } +} + +// Dispatcher function - user specifies precision +void ascii_to_hdf5(const std::string& ascii_file, + const std::string& hdf5_file, + FloatPrecision precision = FloatPrecision::FLOAT64, + bool has_index = false, + unsigned filter_id = 1, + bool verbose = true) +{ + if (precision == FloatPrecision::FLOAT32) { + ascii_to_hdf5_impl(ascii_file, hdf5_file, precision, has_index, filter_id, verbose); + } else { + ascii_to_hdf5_impl(ascii_file, hdf5_file, precision, has_index, filter_id, verbose); + } +} + +// Variant for flexible reading +using FloatData = std::variant, std::vector>; + +struct ParticleDataVariant +{ + std::optional> index; + FloatData m, x, y, z, u, v, w; + std::vector> aux_ints; + std::vector aux_floats; + + size_t num_particles = 0; + size_t num_aux_ints = 0; + size_t num_aux_floats = 0; +}; + +// Detect float precision by inspecting dataset type +FloatPrecision detect_precision(const DataSet& dataset) +{ + const auto datatype = dataset.getDataType(); + if (datatype.getClass() != DataTypeClass::Float) { + throw std::runtime_error("Unexpected dataset type for precision detection (expected float)"); + } + + // Get the size in bytes + const size_t type_size = datatype.getSize(); + + // Float32 is 4 bytes, Float64 is 8 bytes + if (type_size == 4) { + return FloatPrecision::FLOAT32; + } else if (type_size == 8) { + return FloatPrecision::FLOAT64; + } + + throw std::runtime_error("Unexpected float type size: " + std::to_string(type_size)); +} + + +// Read HDF5 with automatic precision detection +ParticleDataVariant read_hdf5_data(const std::string& hdf5_file) +{ + File file(hdf5_file, File::ReadOnly); + + // Read metadata + int num_particles = file.getAttribute("num_particles").read(); + int num_aux_ints = file.getAttribute("num_aux_ints").read(); + int num_aux_floats = file.getAttribute("num_aux_floats").read(); + + // Open particles group + Group particles_group = file.getGroup("particles"); + + ParticleDataVariant data; + data.num_particles = num_particles; + data.num_aux_ints = num_aux_ints; + data.num_aux_floats = num_aux_floats; + + // Lambda to read float or double based on precision + auto read_dataset = [&particles_group](const std::string& name) -> FloatData { + // Detect precision by inspecting dataset + DataSet dataset = particles_group.getDataSet(name); + FloatPrecision precision = detect_precision(dataset); + + if (precision == FloatPrecision::FLOAT32) { + return particles_group.getDataSet(name).read>(); + } else { + return particles_group.getDataSet(name).read>(); + } + }; + + // Read core PSP fields + if (particles_group.exist("index")) { + data.index = particles_group.getDataSet("index").read>(); + } + data.m = read_dataset("m"); + + data.x = read_dataset("x"); + data.y = read_dataset("y"); + data.z = read_dataset("z"); + + data.u = read_dataset("u"); + data.v = read_dataset("v"); + data.w = read_dataset("w"); + + // Read auxiliary integer fields + data.aux_ints.resize(num_aux_ints); + for (int j = 0; j < num_aux_ints; ++j) { + std::string dset_name = "aux_int_" + std::to_string(j); + data.aux_ints[j] = particles_group.getDataSet(dset_name) + .read>(); + } + + // Read auxiliary float fields + data.aux_floats.resize(num_aux_floats); + for (int j = 0; j < num_aux_floats; ++j) { + std::string dset_name = "aux_float_" + std::to_string(j); + data.aux_floats[j] = read_dataset(dset_name); + } + + return data; +} + +// Write ASCII from variant data using std::visit +// +void hdf5_to_ascii(const std::string& hdf5_file, const std::string& ascii_file, + bool verbose = true) +{ + ParticleDataVariant data = read_hdf5_data(hdf5_file); + + std::ofstream outfile(ascii_file); + if (!outfile.is_open()) { + throw std::runtime_error("Could not open ASCII file for writing: " + ascii_file); + } + + // Write header + // + outfile << data.num_particles << " " << data.num_aux_ints + << " " << data.num_aux_floats << "\n"; + outfile.close(); // Close header write + + outfile.precision(16); // High precision for floats + + // Lambda visitor to handle both float and double + // + auto write_value = [](const FloatData& variant_data, size_t index) + { + return std::visit([index](const auto& vec) + { return static_cast(vec[index]); }, variant_data); + }; + + // Strategy: Pre-build all output lines in parallel, then write sequentially + std::vector output_lines(data.num_particles); + + #pragma omp parallel for schedule(static, 256) \ + num_threads(omp_get_max_threads()) + for (int i = 0; i < (int)data.num_particles; ++i) { + std::ostringstream oss; + oss.precision(16); + + // Write optional index and core phase-space fields + if (data.index.has_value()) oss << (*data.index)[i] << " "; + oss << write_value(data.m, i) << " " + << write_value(data.x, i) << " " + << write_value(data.y, i) << " " + << write_value(data.z, i) << " " + << write_value(data.u, i) << " " + << write_value(data.v, i) << " " + << write_value(data.w, i); + + // Write auxiliary integer fields + for (int j = 0; j < data.num_aux_ints; ++j) { + oss << " " << data.aux_ints[j][i]; + } + + // Write auxiliary float fields + for (int j = 0; j < data.num_aux_floats; ++j) { + oss << " " << write_value(data.aux_floats[j], i); + } + + output_lines[i] = oss.str(); + } + + // Write all lines sequentially to avoid I/O contention + outfile.open(ascii_file, std::ios::app); + if (!outfile.is_open()) { + throw std::runtime_error("Could not open ASCII file for appending: " + ascii_file); + } + for (const auto& line : output_lines) { + outfile << line << "\n"; + } + outfile.close(); + + if (verbose) + std::cout << "Successfully wrote " << data.num_particles + << " particles to " << ascii_file << std::endl; +} + +// Main function with command-line parsing +// +int main(int argc, char* argv[]) +{ + std::string prefix = "particles", output; + std::string suffix = "bods"; + int num_threads = omp_get_max_threads(); + unsigned filter_id = 1; + bool quiet = false; + + // Parse command-line arguments for input/output files and mode + // + cxxopts::Options options(argv[0], "ASCII \u2192 HDF5 and HDF5 \u2192 ASCII particle converter for EXP body files\n" + "with built-in round-trip testing and float-size selection\n"); + + options.add_options() + ("i,input", "Input prefix", cxxopts::value(prefix)->default_value("particles")) + ("o,output", "Output prefix (optional, otherwise input prefix is used)", cxxopts::value(output)) + ("a,suffix", "Input suffix", cxxopts::value(suffix)->default_value("bods")) + ("t,threads", "Number of OpenMP threads (default: max available)", cxxopts::value(num_threads)->default_value(std::to_string(num_threads))) + ("f,filter", "HDF5 filter ID to use (default: 1 = GZIP)", cxxopts::value(filter_id)->default_value("1")) + ("roundtrip", "Perform round-trip conversion (ASCII -> HDF5 -> ASCII)") + ("verify", "Verify that the restored ASCII file matches the original in 'roundtrip' mode") + ("to_hdf5", "Convert ASCII to HDF5") + ("double", "Use double precision for HDF5 output (float is default)") + ("filterlist","List available HDF5 filters and exit") + ("aindex", "Input ASCII has leading particle index column; write the particle index dataset to HDF5") + ("to_ascii", "Convert HDF5 to ASCII") + ("q,quiet", "Suppress verbose output") + ("h,help", "Print usage"); + + cxxopts::ParseResult vm; + + try { + vm = options.parse(argc, argv); + } catch (cxxopts::OptionException& e) { + std::cout << "Option error: " << e.what() << std::endl; + exit(-1); + } + + if (vm.count("help")) { + std::cout << options.help() << std::endl; + + // Append note and custom examples + std::cout << R"( +The best performance is achieved using the default 'float' rather than 'double' +precision for the HDF5 file, and that should be sufficient for initial data in +practice. However, the round-trip conversion will be exact only if double +precision is used. + +Examples: + + Convert a standard EXP input body file named 'mybods.bods' to HDF5 format + $ hdf5bods --to_hdf5 -i mybods + The resulting HDF5 file using float32 internally will be called 'mybods.h5' + + Do the same conversion but use full double precision without being chatty + $ hdf5bods --quiet --to_hdf5 --double -i mybods + The resulting HDF5 file will use float64 internally + + Convert the HDF5 file back to the standard EXP ascii format: + $ hdf5bods --to_ascii -i input + + The suffix on input can be customized. For example, the following converts + a standard EXP input body file named 'mybods.asc' to HDF5 format + $ hdf5bods --to_hdf5 --suffix=asc -i mybods + The resulting HDF5 file using float32 format will be called 'mybods.h5' + + Test round-trip conversion and verify that the restored ASCII file matches the + original: + $ hdf5bods --roundtrip --verify -i mybods + +)" << std::endl; + + return 0; + } + + if (vm.count("filterlist")) { + + // Map of commonly registered official HDF5 Filter IDs + std::map known_filters = { + {1, "Deflate / GZIP (Built-in)"}, + {2, "SZIP (Built-in)"}, + {3, "Shuffle (Built-in)"}, + {4, "Fletcher32 Checksum (Built-in)"}, + {32001, "Blosc"}, + {32004, "LZ4"}, + {32008, "Bitshuffle"}, + {32013, "Zfp"}, + {32015, "Zstd"}, + {307, "BZip2"} + }; + + std::cout << "Checking available HDF5 filters on your system:\n"; + std::cout << "-----------------------------------------------\n"; + + for (const auto& [id, name] : known_filters) { + // H5Zfilter_avail returns a positive number if available, 0 if not + htri_t available = H5Zfilter_avail(id); + + if (available > 0) { + std::cout << "[AVAILABLE] ID: " << id << " -> " << name << "\n"; + } else { + std::cout << "[NOT FOUND] ID: " << id << " -> " << name << "\n"; + } + } + return 0; + } + + if (vm.count("quiet")) { + quiet = true; + } + + if (vm.count("threads")) { + omp_set_num_threads(num_threads); + } + + if (!quiet) + std::cout << "=== Using " << num_threads << " OpenMP threads ===" + << std::endl; + + // This is for testing and timing + // + if (vm.count("roundtrip")) { + + // Time the whole test + std::unique_ptr time_ptr; + if (!quiet) time_ptr = std::make_unique("full test"); + + // Filenames + std::string ascii = prefix + "." + suffix; + std::string h5_32 = prefix + "_f32.h5"; + std::string h5_64 = prefix + "_f64.h5"; + std::string rest = prefix + "_restored." + suffix; + + try { + // Convert to HDF5 with float32 precision + std::cout << "=== Converting ASCII to HDF5 (float32) ===" << std::endl; + { + std::unique_ptr time_ptr; + if (!quiet) time_ptr = std::make_unique("float 32 conversion"); + ascii_to_hdf5(ascii, h5_32, FloatPrecision::FLOAT32, vm.count("aindex"), filter_id, !quiet); + } + + // Convert to HDF5 with float64 precision + std::cout << "\n=== Converting ASCII to HDF5 (float64) ===" << std::endl; + { + std::unique_ptr time_ptr; + if (!quiet) time_ptr = std::make_unique("float 64 conversion"); + ascii_to_hdf5(ascii, h5_64, FloatPrecision::FLOAT64, vm.count("aindex"), filter_id, !quiet); + } + + // Round-trip with float64 version + std::cout << "\n=== Converting HDF5 (float64) back to ASCII ===" << std::endl; + { + std::unique_ptr time_ptr; + if (!quiet) time_ptr = std::make_unique("hdf5 to ascii"); + hdf5_to_ascii(h5_64, rest, !quiet); + } + + std::cout << "\n=== Round-trip conversion complete! ===" << std::endl; + + if (vm.count("verify")) { + std::cout << "\n=== Verifying restored ASCII file against original ===" << std::endl; + std::ifstream original(ascii); + std::ifstream restored(rest ); + + if (!original.is_open() || !restored.is_open()) { + throw std::runtime_error("Could not open files for verification"); + } + + std::string line_orig, line_rest; + size_t line_num = 0; + + std::vector max_diff(7, 0.0); // For m, x, y, z, u, v, w + std::vector vec_orig(7), vec_rest(7); + std::vector max_aux_int_diff; + std::vector max_aux_float_diff; + std::vector aux_int_orig, aux_int_rest; + std::vector aux_float_orig, aux_float_rest; + int num_aux_ints = -1, num_aux_floats = -1; + + while (std::getline(original, line_orig) && std::getline(restored, line_rest)) { + ++line_num; + + std::istringstream ss_orig(line_orig), ss_rest(line_rest); + + if (line_num == 1) { + int total_particles_orig, total_aux_ints_orig, total_aux_floats_orig; + int total_particles_rest, total_aux_ints_rest, total_aux_floats_rest; + + ss_orig >> total_particles_orig >> total_aux_ints_orig >> total_aux_floats_orig; + ss_rest >> total_particles_rest >> total_aux_ints_rest >> total_aux_floats_rest; + + if (total_particles_orig != total_particles_rest || + total_aux_ints_orig != total_aux_ints_rest || + total_aux_floats_orig != total_aux_floats_rest) { + std::cerr << "Header mismatch at line 1:\n" + << "Original: " << line_orig << "\n" + << "Restored: " << line_rest << "\n"; + break; + } + + num_aux_ints = total_aux_ints_orig; + num_aux_floats = total_aux_floats_orig; + max_aux_int_diff.assign(num_aux_ints, 0.0); + max_aux_float_diff.assign(num_aux_floats, 0.0); + aux_int_orig.assign(num_aux_ints, 0); + aux_int_rest.assign(num_aux_ints, 0); + aux_float_orig.assign(num_aux_floats, 0.0); + aux_float_rest.assign(num_aux_floats, 0.0); + + continue; + } + + if (vm.count("aindex")) { + unsigned long idx_orig = 0, idx_rest = 0; + if (!(ss_orig >> idx_orig) || !(ss_rest >> idx_rest)) { + throw std::runtime_error("Failed to parse index field at line " + std::to_string(line_num)); + } + } + + for (int i = 0; i < 7; ++i) { + if (!(ss_orig >> vec_orig[i]) || !(ss_rest >> vec_rest[i])) { + throw std::runtime_error("Failed to parse core field at line " + std::to_string(line_num)); + } + } + + for (int i = 0; i < 7; ++i) { + if (std::abs(vec_orig[i] - vec_rest[i]) > max_diff[i]) { + max_diff[i] = std::abs(vec_orig[i] - vec_rest[i]); + } + + } + + for (int i = 0; i < num_aux_ints; ++i) { + if (!(ss_orig >> aux_int_orig[i]) || !(ss_rest >> aux_int_rest[i])) { + throw std::runtime_error("Failed to parse auxiliary integer field at line " + std::to_string(line_num)); + } + } + + for (int i = 0; i < num_aux_ints; ++i) { + auto diff_i = std::llabs(static_cast(aux_int_orig[i]) - + static_cast(aux_int_rest[i])); + double diff = static_cast(diff_i); + if (diff > max_aux_int_diff[i]) max_aux_int_diff[i] = diff; + } + + for (int i = 0; i < num_aux_floats; ++i) { + if (!(ss_orig >> aux_float_orig[i]) || !(ss_rest >> aux_float_rest[i])) { + throw std::runtime_error("Failed to parse auxiliary float field at line " + std::to_string(line_num)); + } + } + + for (int i = 0; i < num_aux_floats; ++i) { + double diff = std::abs(aux_float_orig[i] - aux_float_rest[i]); + if (diff > max_aux_float_diff[i]) max_aux_float_diff[i] = diff; + } + } + + std::cout << "Maximum absolute differences for core fields" + << std::endl << "(m, x, y, z, u, v, w): " + << max_diff[0] << ", " + << max_diff[1] << ", " + << max_diff[2] << ", " + << max_diff[3] << ", " + << max_diff[4] << ", " + << max_diff[5] << ", " + << max_diff[6] << std::endl; + + std::cout << "Maximum absolute differences for auxiliary integer fields: "; + if (max_aux_int_diff.empty()) { + std::cout << "(none)"; + } else { + for (size_t i = 0; i < max_aux_int_diff.size(); ++i) { + if (i) std::cout << ", "; + std::cout << max_aux_int_diff[i]; + } + } + std::cout << std::endl; + + std::cout << "Maximum absolute differences for auxiliary float fields: "; + if (max_aux_float_diff.empty()) { + std::cout << "(none)"; + } else { + for (size_t i = 0; i < max_aux_float_diff.size(); ++i) { + if (i) std::cout << ", "; + std::cout << max_aux_float_diff[i]; + } + } + std::cout << std::endl; + } + + } catch (const std::exception& e) { + std::cerr << "Error: " << e.what() << std::endl; + return 1; + } + + + std::cout << "\n=== Compression ratios ===" << std::endl; + + namespace fs = std::filesystem; + + try { + fs::path ascPath = ascii; + fs::path f32Path = h5_32; + fs::path f64Path = h5_64; + + // Returns size in bytes (std::uintmax_t) + std::uintmax_t ascSize = fs::file_size(ascPath); + std::uintmax_t f32Size = fs::file_size(f32Path); + std::uintmax_t f64Size = fs::file_size(f64Path); + + // Print compression ratios + std::cout << "Ascii/hdf5(float) = " + << std::setprecision(3) + << static_cast(ascSize)/f32Size + << std::endl; + + std::cout << "Ascii/hdf5(double) = " + << std::setprecision(3) + << static_cast(ascSize)/f64Size + << std::endl; + std::cout << std::endl; + } + catch (const fs::filesystem_error& e) { + std::cerr << "Error: " << e.what() << '\n'; + } + + return 0; + } + + // Production mode: convert ASCII to HDF5 or vice versa based on + // mode selection + // + if (vm.count("to_hdf5")) { + std::unique_ptr time_ptr; + if (!quiet) time_ptr = std::make_unique("conversion to hdf5"); + std::string ascii = prefix + "." + suffix; + std::string ofile = prefix + ".h5"; + if (vm.count("output")) ofile = output + ".h5"; + if (vm.count("double")) + ascii_to_hdf5(ascii, ofile, FloatPrecision::FLOAT64, vm.count("aindex"), filter_id, !quiet); + else + ascii_to_hdf5(ascii, ofile, FloatPrecision::FLOAT32, vm.count("aindex"), filter_id, !quiet); + } else if (vm.count("to_ascii")) { + std::unique_ptr time_ptr; + if (!quiet) time_ptr = std::make_unique("conversion to ascii"); + std::string ifile = prefix + ".h5"; + std::string ofile = prefix + "." + suffix; + hdf5_to_ascii(ifile, ofile, !quiet); + } else { + std::cerr << "No conversion mode specified. Use --to_hdf5 to convert to HDF5, --to_ascii to convert from HDF5 to ascii, or the --roundtrip test." << std::endl; + return 1; + } + + return 0; +}