From 06eb7f2439912eb8767eda8b55b8edc6fa9a4f07 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 22 Jul 2026 23:23:52 -0400 Subject: [PATCH 1/7] Fix heap-buffer-overflow in Eigen Map/TensorMap cached pointers --- expui/BiorthBasis.cc | 68 +++--- expui/CoefContainer.cc | 58 ++--- expui/CoefStruct.H | 466 +++++++++++++++++++++++++++++++++++------ expui/Coefficients.cc | 165 +++++++-------- expui/FieldBasis.H | 2 +- expui/FieldBasis.cc | 26 ++- src/Cube.cc | 2 +- src/Cylinder.cc | 2 +- src/PolarBasis.cc | 2 +- src/SlabSL.cc | 2 +- src/SphericalBasis.cc | 2 +- 11 files changed, 568 insertions(+), 227 deletions(-) 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..047654349 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[3]).imag(); } } } @@ -297,8 +297,8 @@ namespace MSSA int f = k[1], 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..4c2f692d4 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,34 @@ 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->resize(rows, cols); + } + return *coefs; + } + + //! Non-const version for assignment + coefType& getCoefMapMutable() + { + 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(), 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 +183,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 +219,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 +240,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 +277,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 +313,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 +337,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 +368,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 +378,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 +413,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 +437,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 +465,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 +483,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 +516,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 +536,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::MatrixXcd& mat) + { + if (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 collection of trajectories class TrajStruct : public CoefStruct { @@ -386,15 +599,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 +623,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 +694,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 +721,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 +762,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 +797,6 @@ namespace CoefClasses //! Coefficient map type using coefType = Eigen::TensorMap; - //! Coefficient map - std::shared_ptr coefs; - //! Number of data fields int nfld; @@ -506,7 +805,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 +825,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 +864,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 +905,3 @@ namespace CoefClasses #endif - diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc index da76c6f28..0a4ad2463 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()(0, 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..e142be086 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,15 @@ namespace BasisClasses { if (dof==2) { 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: CylFldStruct 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->mmax || nmax!=p->nmax) { @@ -300,8 +307,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/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 Date: Thu, 23 Jul 2026 09:12:48 -0400 Subject: [PATCH 2/7] Typo fix in TableData update for new interface --- expui/Coefficients.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc index 0a4ad2463..4613a80c9 100644 --- a/expui/Coefficients.cc +++ b/expui/Coefficients.cc @@ -2715,7 +2715,7 @@ namespace CoefClasses coef->time = times[n]; coef->store.resize(cols); coef->initCoefMap(); - for (int i=0; igetCoefs()(0, i) = data[n][i]; + for (int i=0; igetCoefs()(i) = data[n][i]; coefs[roundTime(times[n])] = coef; } From f87530a2349cf1a32d8d76b5c1182f9240cc4e7f Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Thu, 23 Jul 2026 14:42:32 -0400 Subject: [PATCH 3/7] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- expui/CoefStruct.H | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/expui/CoefStruct.H b/expui/CoefStruct.H index 4c2f692d4..903c3d706 100644 --- a/expui/CoefStruct.H +++ b/expui/CoefStruct.H @@ -151,9 +151,10 @@ namespace CoefClasses const coefType& getCoefMap() const { if (!coefs) { - int rows = (lmax+1)*(lmax+2)/2; - int cols = nmax; - coefs->resize(rows, cols); + int rows = (lmax+1)*(lmax+2)/2; + int cols = nmax; + coefs = std::make_shared + (const_cast*>(store.data()), rows, cols); } return *coefs; } @@ -161,10 +162,11 @@ namespace CoefClasses //! Non-const version for assignment coefType& getCoefMapMutable() { - 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(), rows, cols); + if (!coefs) { + int rows = (lmax+1)*(lmax+2)/2; + int cols = nmax; + coefs = std::make_shared(store.data(), rows, cols); + } return *coefs; } From 1c6f7b72fe585c6096a932f48cbd9627519920f0 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Thu, 23 Jul 2026 14:43:15 -0400 Subject: [PATCH 4/7] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- expui/FieldBasis.cc | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc index e142be086..dec17b60b 100644 --- a/expui/FieldBasis.cc +++ b/expui/FieldBasis.cc @@ -286,14 +286,15 @@ namespace BasisClasses if (dof==2) { auto p = dynamic_cast(c.get()); + auto& src = p->getCoefs(); + // 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: CylFldStruct data size mismatch with FieldBasis store"); + 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(store[0].data(), store[0].data() + store[0].size(), - p->getCoefs().data()); + std::copy(src.data(), src.data() + src.size(), store[0].data()); // Sanity test dimensions if (nfld!=p->nfld || lmax!=p->mmax || nmax!=p->nmax) { From c01f44b310f46ccdb7b52f49a82e7b4cb1033016 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Thu, 23 Jul 2026 15:09:37 -0400 Subject: [PATCH 5/7] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- expui/CoefContainer.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/CoefContainer.cc b/expui/CoefContainer.cc index 047654349..7b99230ff 100644 --- a/expui/CoefContainer.cc +++ b/expui/CoefContainer.cc @@ -295,7 +295,7 @@ 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->getCoefs())(f, m, n) = {data[c][i], 0.0}; else (cf->getCoefs())(f, m, n) = {data[c][i], data[s][i]}; From 0cf62062dd1aab0211b163b2b858e70d2b69d81f Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 23 Jul 2026 19:09:20 +0000 Subject: [PATCH 6/7] Fix TblStruct setCoefs vector type mismatch --- expui/CoefStruct.H | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/expui/CoefStruct.H b/expui/CoefStruct.H index 903c3d706..8a65d975c 100644 --- a/expui/CoefStruct.H +++ b/expui/CoefStruct.H @@ -579,11 +579,11 @@ namespace CoefClasses } //! Set the underlying coefficient data - void setCoefs(const Eigen::MatrixXcd& mat) + void setCoefs(const Eigen::VectorXcd& vec) { - if (mat.cols() != cols) + if (vec.size() != cols) throw std::invalid_argument("setCoefs: size mismatch"); - getCoefMapMutable() = mat; + getCoefMapMutable() = vec; } //@{ From 1bc38a91110657fbaef75b585602526a05ceb3bc Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 23 Jul 2026 19:11:34 +0000 Subject: [PATCH 7/7] Fix cylfld imag index typo in pack_cylfld --- expui/CoefContainer.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/CoefContainer.cc b/expui/CoefContainer.cc index 7b99230ff..55c49e5d6 100644 --- a/expui/CoefContainer.cc +++ b/expui/CoefContainer.cc @@ -280,7 +280,7 @@ namespace MSSA if (k[3]==0) data[k][t] = (cf->getCoefs())(k[0], k[1], k[2]).real(); else - data[k][t] = (cf->getCoefs())(k[0], k[1], k[3]).imag(); + data[k][t] = (cf->getCoefs())(k[0], k[1], k[2]).imag(); } } }