From 22146f6ee0aac92b328f071e900a70599ba39222 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 27 Mar 2026 11:15:41 -0400 Subject: [PATCH 01/15] Manual cherry pick of chunk logic fix to devel --- src/Component.cc | 24 ++++++++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/src/Component.cc b/src/Component.cc index dfcf394e1..c50884ce6 100644 --- a/src/Component.cc +++ b/src/Component.cc @@ -2501,13 +2501,23 @@ void Component::write_HDF5(HighFive::Group& group, bool masses, bool IDs) auto dcplI = HighFive::DataSetCreateProps{}; auto dcplD = HighFive::DataSetCreateProps{}; - if (H5compress or H5chunk) { + // Compression and chunking. Do not set chunk size larger than + // nbodies. Turn off compression altogether if nbodies = 0 to avoid + // HDF5 errors. + // + if ((H5compress or H5chunk) and nbodies > 0) { int chunk = H5chunk; - // Sanity + // Clamp chunk to [1, nbodies]: use nbodies/8 as a downsize when + // H5chunk would exceed the dataset extent, then ensure at least 1 if (H5chunk >= nbodies) { chunk = nbodies/8; } + if (chunk < 1) { + chunk = 1; + } else if (static_cast(chunk) > nbodies) { + chunk = static_cast(nbodies); + } dcpl1.add(HighFive::Chunking(chunk)); if (H5shuffle) dcpl1.add(HighFive::Shuffle()); @@ -2636,14 +2646,20 @@ void Component::write_H5(H5::Group& group) // This could be generalized by registering a user filter, like // blosc. Right now, we're using the default (which is gzip) - if (H5compress or H5chunk) { + // + // Do not set chunk size larger than number of particles. If the + // particle number is zero, do not compress. + // + if ((H5compress or H5chunk) and h5_particles.size() > 0) { // Set chunking if (H5chunk) { - // Sanity + // Clamp chunk to [1, nbodies]: use nbodies/8 as a downsize when + // H5chunk would exceed the dataset extent, then ensure at least 1 int chunk = H5chunk; if (H5chunk >= nbodies) { chunk = nbodies/8; } + chunk = std::clamp(chunk, 1, static_cast(nbodies)); hsize_t chunk_dims[1] = {static_cast(chunk)}; dcpl.setChunk(1, chunk_dims); } From 1952df71ab4641f010c7ae5a446a8811652da5ae Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 27 Mar 2026 12:14:43 -0400 Subject: [PATCH 02/15] Version bump for bug fix --- CMakeLists.txt | 2 +- doc/exp.cfg | 2 +- doc/exp.cfg.breathe | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index e2d49dba5..9d3dd2947 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.1" + VERSION "7.10.2" HOMEPAGE_URL https://github.com/EXP-code/EXP LANGUAGES C CXX Fortran) diff --git a/doc/exp.cfg b/doc/exp.cfg index d6a169162..c930e7812 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.1 +PROJECT_NUMBER = 7.10.2 # 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 95375995d..96b4b616b 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.1 +PROJECT_NUMBER = 7.10.2 # 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 From 3d4ef6848066bf9bd52b811ab3fbd4d51339a43c Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 18 Apr 2026 12:15:25 -0400 Subject: [PATCH 03/15] Attempt to reinstate covariance support in pyEXP after, what looks like, a merge error --- expui/BasisFactory.H | 67 ++++++ expui/BiorthBasis.H | 62 ------ pyEXP/BasisWrappers.cc | 493 +++++++++++++++++++++++++++++++---------- 3 files changed, 442 insertions(+), 180 deletions(-) diff --git a/expui/BasisFactory.H b/expui/BasisFactory.H index 9ea572aff..35ac359f0 100644 --- a/expui/BasisFactory.H +++ b/expui/BasisFactory.H @@ -10,6 +10,7 @@ #include "OrthoFunction.H" #include "Coefficients.H" #include "PseudoAccel.H" +#include "Covariance.H" #include "YamlCheck.H" #include "localmpi.H" #include "exputils.H" @@ -135,6 +136,32 @@ namespace BasisClasses Eigen::MatrixXd p_accel; //@} + //@{ + //! Sample counts and masses for covariance computation + Eigen::VectorXi sampleCounts; + Eigen::VectorXd sampleMasses; + //@} + + //! Covariance storage instance + std::shared_ptr covarStore; + + //! Coefficient variance computation enabled + bool pcavar = false; + + //! Round time key to emulated fixed-point arithmetic + double roundTime(double time) + { + // Eight decimal places should be enough here... + const double multiplier = 1.0e+08; // std::pow(10.0, 8); + return std::floor(time * multiplier + 0.5) / multiplier; + } + + //! Store covariance matrix? + bool covar = true; + + //! Store summed covariance only? + bool scovr = true; + //! Number of center points in acceleration estimator int Naccel = 0; @@ -291,6 +318,46 @@ namespace BasisClasses //! Get the basis expansion center RowMatrix3d getRotation() { return coefrot; } + + //!! Sample counts, masses, coefficients, and covariance + virtual SubsampleCovariance::CovarData getCoefCovariance(double time) + { + if (!covarStore) + throw std::runtime_error("Basis::getCoefCovariance: covariance storage not initialized"); + return covarStore->getCoefCovariance(time); + } + + //! Write coefficient covariance data to an HDF5 file + virtual void writeCoefCovariance(const std::string& compname, const std::string& runtag, + double time=0.0) + { + // Must be overriden; base implementation throws error + throw std::runtime_error("Basis::writeCoefCovariance: " + "Not implemented for this basis"); + } + + //! Make covariance after accumulation + virtual void makeCoefCovariance(void) {} + + //! Enable covariance computation with optional sample time + virtual void enableCoefCovariance(bool pcavar, int sampT_in, bool ftype, bool total, bool covar_in) + { + // Must be overriden; base implementation throws error + throw std::runtime_error("Basis::enableCoefCovariance: " + "Not implemented for this basis"); + } + + //! HDF5 compression settings + void setCovarH5Compress(unsigned level, unsigned chunksize, bool shuffle, + bool szip=false) + { + if (covarStore) { + covarStore->setCovarH5Compress(level, chunksize, shuffle, szip); + } else { + throw std::runtime_error("BiorthBasis::setCovarH5Compress: covariance storage not initialized"); + } + } + }; using BasisPtr = std::shared_ptr; diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index f772ad9ba..3bf23168a 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -21,7 +21,6 @@ #include "BiorthBess.H" #include "BasisFactory.H" #include "BiorthCube.H" -#include "Covariance.H" #include "sltableMP2.H" #include "SLGridMP2.H" #include "YamlCheck.H" @@ -104,32 +103,6 @@ namespace BasisClasses RowMatrixXd varray; //@} - //! Coefficient variance computation enabled - bool pcavar = false; - - //@{ - //! Sample counts and masses for covariance computation - Eigen::VectorXi sampleCounts; - Eigen::VectorXd sampleMasses; - //@} - - //! Covariance storage instance - std::shared_ptr covarStore; - - //! Round time key to emulated fixed-point arithmetic - double roundTime(double time) - { - // Eight decimal places should be enough here... - const double multiplier = 1.0e+08; // std::pow(10.0, 8); - return std::floor(time * multiplier + 0.5) / multiplier; - } - - //! Store covariance matrix? - bool covar = true; - - //! Store summed covariance only? - bool scovr = true; - public: //! Basis identifier string @@ -273,41 +246,6 @@ namespace BasisClasses return varray; } - //!! Sample counts, masses, coefficients, and covariance - SubsampleCovariance::CovarData getCoefCovariance(double time) - { return covarStore->getCoefCovariance(time); } - - //! Write coefficient covariance data to an HDF5 file - virtual void writeCoefCovariance(const std::string& compname, const std::string& runtag, - double time=0.0) - { - // Must be overriden; base implementation throws error - throw std::runtime_error("BiorthBasis::writeCoefCovariance: " - "Not implemented for this basis"); - } - - //! Make covariance after accumulation - virtual void makeCoefCovariance(void) {} - - //! Enable covariance computation with optional sample time - virtual void enableCoefCovariance(bool pcavar, int sampT_in, bool ftype, bool covar_in) - { - // Must be overriden; base implementation throws error - throw std::runtime_error("BiorthBasis::enableCoefCovariance: " - "Not implemented for this basis"); - } - - //! HDF5 compression settings - void setCovarH5Compress(unsigned level, unsigned chunksize, bool shuffle, - bool szip=false) - { - if (covarStore) { - covarStore->setCovarH5Compress(level, chunksize, shuffle, szip); - } else { - throw std::runtime_error("BiorthBasis::setCovarH5Compress: covariance storage not initialized"); - } - } - //! Evaluate acceleration in Cartesian coordinates in centered //! coordinate system for a collections of points virtual RowMatrixXd& getAccel(Eigen::Ref pos) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 1a997a005..a31ce4987 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -311,6 +311,26 @@ void BasisFactoryClasses(py::module &m) addFromArray(Eigen::VectorXd& m, RowMatrixXd& p, bool roundrobin, bool posvelrows) override { PYBIND11_OVERRIDE_PURE(void, Basis, addFromArray, m, p, roundrobin, posvelrows); } + + virtual SubsampleCovariance::CovarData + getCoefCovariance(double time) override { + PYBIND11_OVERRIDE(SubsampleCovariance::CovarData, Basis, getCoefCovariance, time); + } + + virtual void writeCoefCovariance + (const std::string& file, const std::string& runtag, double time) override { + PYBIND11_OVERRIDE(void, Basis, writeCoefCovariance, file, runtag, time); + } + + virtual void makeCoefCovariance(void) override { + PYBIND11_OVERRIDE(void, Basis, makeCoefCovariance,); + } + + virtual void + enableCoefCovariance(bool pcavar, int nsamples, bool ftype, bool total, bool covar) { + PYBIND11_OVERRIDE(void, Basis, enableCoefCovariance, pcavar, nsamples, ftype, total, covar); + } + }; class PyFieldBasis : public FieldBasis @@ -505,8 +525,15 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE_PURE(std::vector, Spherical, orthoCheck, knots); } - void enableCoefCovariance(bool pcavar, int nsamples, bool ftype, bool covar) override { - PYBIND11_OVERRIDE(void, Spherical, enableCoefCovariance, pcavar, nsamples, ftype, covar); + virtual void writeCoefCovariance + (const std::string& file, const std::string& runtag, double time) override { + PYBIND11_OVERRIDE(void, Spherical, writeCoefCovariance, file, runtag, time); + } + + + virtual void enableCoefCovariance + (bool pcavar, int nsamples, bool ftype, bool total, bool covar) override { + PYBIND11_OVERRIDE(void, Spherical, enableCoefCovariance, pcavar, nsamples, ftype, total, covar); } }; @@ -560,8 +587,14 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE(void, Cylindrical, make_coefs,); } - void enableCoefCovariance(bool pcavar, int nsamples, bool ftype, bool covar) override { - PYBIND11_OVERRIDE(void, Cylindrical, enableCoefCovariance, pcavar, nsamples, ftype, covar); + virtual void writeCoefCovariance + (const std::string& file, const std::string& runtag, double time) override { + PYBIND11_OVERRIDE(void, Cylindrical, writeCoefCovariance, file, runtag, time); + } + + virtual void enableCoefCovariance + (bool pcavar, int nsamples, bool ftype, bool total, bool covar) override { + PYBIND11_OVERRIDE(void, Cylindrical, enableCoefCovariance, pcavar, nsamples, ftype, total, covar); } }; @@ -636,6 +669,17 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE(void, FlatDisk, make_coefs,); } + virtual void writeCoefCovariance + (const std::string& file, const std::string& runtag, double time) override { + PYBIND11_OVERRIDE(void, FlatDisk, writeCoefCovariance, file, runtag, time); + } + + virtual void enableCoefCovariance + (bool pcavar, int nsamples, bool ftype, bool total, bool covar) override { + PYBIND11_OVERRIDE(void, FlatDisk, enableCoefCovariance, pcavar, nsamples, ftype, total, covar); + } + + }; @@ -855,8 +899,14 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE(void, Cube, make_coefs,); } - void enableCoefCovariance(bool pcavar, int nsamples, bool ftype, bool covar) override { - PYBIND11_OVERRIDE(void, Cube, enableCoefCovariance, pcavar, nsamples, ftype, covar); + virtual void writeCoefCovariance + (const std::string& file, const std::string& runtag, double time) override { + PYBIND11_OVERRIDE(void, Cube, writeCoefCovariance, file, runtag, time); + } + + virtual void enableCoefCovariance + (bool pcavar, int nsamples, bool ftype, bool total, bool covar) override { + PYBIND11_OVERRIDE(void, Cube, enableCoefCovariance, pcavar, nsamples, ftype, total, covar); } }; @@ -1333,82 +1383,7 @@ void BasisFactoryClasses(py::module &m) makeFromArray : create coefficients contributions )", py::arg("mass"), py::arg("pos")) - .def("getCoefCovariance", - [](BasisClasses::BiorthBasis& A, double time) - { - auto [cnts, mass, coef, covr] = A.getCoefCovariance(time); - py::array_t> cf = make_ndarray3>(coef); - py::array_t> vr = make_ndarray4>(covr); - return std::make_tuple(cnts, mass, cf, vr); - }, - R"( - Get the covariance matrices for the basis coefficients - - Parameters - ---------- - time : float - the evaluation time - - Returns - ------- - tuple(numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray) - tuple of counts, masses, partitioned coefficients and their covariance - matrices for each subsample. The returns are complex-valued. - )", - py::arg("time")) - .def("writeCoefCovariance", &BasisClasses::BiorthBasis::writeCoefCovariance, - R"( - Write the partitioned coefficient vectors and covariance matrices - to a specified HDF5 file. The number of partitions is set by the - configuration parameter 'sampT' and defaults to 100. - - On first call, the file is created, metadata is written, and the - coefficient vectors and covariance matrices are stored. On subsequent - calls, the file is updated with new covariance datasets. - - The file will be called 'coefcovar...h5'. - - Parameters - ---------- - compname : str - the component/basis name segment used in the output HDF5 filename - runtag : str - the run identifier tag - time : float - the snapshot time - - Returns - ------- - None - - See also - -------- - getCoefCovariance : get the counts, mass, coefficient vectors and covariance matrices - )", py::arg("compname"), py::arg("runtag"), py::arg("time")=0.0) - .def("enableCoefCovariance", &BasisClasses::BiorthBasis::enableCoefCovariance, - R"( - Enable or disable the coefficient covariance computation and set the - default number of partitions to use for the covariance computation. - - Parameters - ---------- - pcavar : bool - enable (true) or disable (false) the covariance computation - nsamples : int - number of time partitions to use for covariance computation - ftype: bool - if true, use float32 for covariance storage; if false, - use float64 (default: false) - covar: bool - if true, compute and save covariance to the HDF5 file; if false, - save mean and variance vectors only (default: true) - - Returns - ------- - None - )", py::arg("pcavar"), py::arg("nsamples")=100, py::arg("ftype")=false, - py::arg("covar")=true) - .def("setCovarH5Compress", &BasisClasses::BiorthBasis::setCovarH5Compress, + .def("setCovarH5Compress", &BasisClasses::Basis::setCovarH5Compress, R"( Set the HDF5 compression level for covariance storage in HDF5. The Szip compression algorithm may also be enabled but seems to not have better @@ -1857,7 +1832,85 @@ void BasisFactoryClasses(py::module &m) dict({tag: value},...) cache parameters )", - py::arg("cachefile")); + py::arg("cachefile")) + .def("getCoefCovariance", + [](BasisClasses::Cylindrical& A, double time) + { + auto [cnts, mass, coef, covr] = A.getCoefCovariance(time); + py::array_t> cf = make_ndarray3>(coef); + py::array_t> vr = make_ndarray4>(covr); + return std::make_tuple(cnts, mass, cf, vr); + }, + R"( + Get the covariance matrices for the basis coefficients + + Parameters + ---------- + time : float + the evaluation time + + Returns + ------- + tuple(numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray) + tuple of counts, masses, partitioned coefficients and their covariance + matrices for each subsample. The returns are complex-valued. + )", + py::arg("time")) + .def("writeCoefCovariance", &BasisClasses::Cylindrical::writeCoefCovariance, + R"( + Write the partitioned coefficient vectors and covariance matrices + to a specified HDF5 file. The number of partitions is set by the + configuration parameter 'sampT' and defaults to 100. + + On first call, the file is created, metadata is written, and the + coefficient vectors and covariance matrices are stored. On subsequent + calls, the file is updated with new covariance datasets. + + The file will be called 'coefcovar...h5'. + + Parameters + ---------- + compname : str + the component/basis name segment used in the output HDF5 filename + runtag : str + the run identifier tag + time : float + the snapshot time + + Returns + ------- + None + + See also + -------- + getCoefCovariance : get the counts, mass, coefficient vectors and covariance matrices + )", py::arg("compname"), py::arg("runtag"), py::arg("time")=0.0) + .def("enableCoefCovariance", &BasisClasses::Cylindrical::enableCoefCovariance, + R"( + Enable or disable the coefficient covariance computation and set the + default number of partitions to use for the covariance computation. + + Parameters + ---------- + pcavar : bool + enable (true) or disable (false) the covariance computation + nsamples : int + number of time partitions to use for covariance computation + ftype: bool + if true, use float32 for covariance storage; if false, + use float64 (default: false) + total: bool + if true, also compute the total covariance matrix; if false, save only + the partitioned covariance matrices (default: true) + covar: bool + if true, compute and save covariance to the HDF5 file; if false, + save mean and variance vectors only (default: true) + + Returns + ------- + None + )", py::arg("pcavar"), py::arg("nsamples")=100, py::arg("ftype")=false, + py::arg("total")=true, py::arg("covar")=true); py::class_, PySpherical, BasisClasses::BiorthBasis>(m, "Spherical") @@ -2091,7 +2144,85 @@ void BasisFactoryClasses(py::module &m) list(numpy.ndarray) list of numpy.ndarrays from [0, ... , Lmax] )", - py::arg("knots")=40); + py::arg("knots")=40) + .def("getCoefCovariance", + [](BasisClasses::SphericalSL& A, double time) + { + auto [cnts, mass, coef, covr] = A.getCoefCovariance(time); + py::array_t> cf = make_ndarray3>(coef); + py::array_t> vr = make_ndarray4>(covr); + return std::make_tuple(cnts, mass, cf, vr); + }, + R"( + Get the covariance matrices for the basis coefficients + + Parameters + ---------- + time : float + the evaluation time + + Returns + ------- + tuple(numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray) + tuple of counts, masses, partitioned coefficients and their covariance + matrices for each subsample. The returns are complex-valued. + )", + py::arg("time")) + .def("writeCoefCovariance", &BasisClasses::SphericalSL::writeCoefCovariance, + R"( + Write the partitioned coefficient vectors and covariance matrices + to a specified HDF5 file. The number of partitions is set by the + configuration parameter 'sampT' and defaults to 100. + + On first call, the file is created, metadata is written, and the + coefficient vectors and covariance matrices are stored. On subsequent + calls, the file is updated with new covariance datasets. + + The file will be called 'coefcovar...h5'. + + Parameters + ---------- + compname : str + the component/basis name segment used in the output HDF5 filename + runtag : str + the run identifier tag + time : float + the snapshot time + + Returns + ------- + None + + See also + -------- + getCoefCovariance : get the counts, mass, coefficient vectors and covariance matrices + )", py::arg("compname"), py::arg("runtag"), py::arg("time")=0.0) + .def("enableCoefCovariance", &BasisClasses::SphericalSL::enableCoefCovariance, + R"( + Enable or disable the coefficient covariance computation and set the + default number of partitions to use for the covariance computation. + + Parameters + ---------- + pcavar : bool + enable (true) or disable (false) the covariance computation + nsamples : int + number of time partitions to use for covariance computation + ftype: bool + if true, use float32 for covariance storage; if false, + use float64 (default: false) + total: bool + if true, also compute the total covariance matrix; if false, save only + the partitioned covariance matrices (default: true) + covar: bool + if true, compute and save covariance to the HDF5 file; if false, + save mean and variance vectors only (default: true) + + Returns + ------- + None + )", py::arg("pcavar"), py::arg("nsamples")=100, py::arg("ftype")=false, + py::arg("total")=true, py::arg("covar")=true); py::class_, BasisClasses::Spherical>(m, "Bessel") @@ -2221,7 +2352,86 @@ void BasisFactoryClasses(py::module &m) out : dict({tag: value}) cache parameters )", - py::arg("cachefile")); + py::arg("cachefile")) + .def("getCoefCovariance", + [](BasisClasses::FlatDisk& A, double time) + { + auto [cnts, mass, coef, covr] = A.getCoefCovariance(time); + py::array_t> cf = make_ndarray3>(coef); + py::array_t> vr = make_ndarray4>(covr); + return std::make_tuple(cnts, mass, cf, vr); + }, + R"( + Get the covariance matrices for the basis coefficients + + Parameters + ---------- + time : float + the evaluation time + + Returns + ------- + tuple(numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray) + tuple of counts, masses, partitioned coefficients and their covariance + matrices for each subsample. The returns are complex-valued. + )", + py::arg("time")) + .def("writeCoefCovariance", &BasisClasses::FlatDisk::writeCoefCovariance, + R"( + Write the partitioned coefficient vectors and covariance matrices + to a specified HDF5 file. The number of partitions is set by the + configuration parameter 'sampT' and defaults to 100. + + On first call, the file is created, metadata is written, and the + coefficient vectors and covariance matrices are stored. On subsequent + calls, the file is updated with new covariance datasets. + + The file will be called 'coefcovar...h5'. + + Parameters + ---------- + compname : str + the component/basis name segment used in the output HDF5 filename + runtag : str + the run identifier tag + time : float + the snapshot time + + Returns + ------- + None + + See also + -------- + getCoefCovariance : get the counts, mass, coefficient vectors and covariance matrices + )", py::arg("compname"), py::arg("runtag"), py::arg("time")=0.0) + .def("enableCoefCovariance", &BasisClasses::FlatDisk::enableCoefCovariance, + R"( + Enable or disable the coefficient covariance computation and set the + default number of partitions to use for the covariance computation. + + Parameters + ---------- + pcavar : bool + enable (true) or disable (false) the covariance computation + nsamples : int + number of time partitions to use for covariance computation + ftype: bool + if true, use float32 for covariance storage; if false, + use float64 (default: false) + total: bool + if true, also compute the total covariance matrix; if false, save only + the partitioned covariance matrices (default: true) + covar: bool + if true, compute and save covariance to the HDF5 file; if false, + save mean and variance vectors only (default: true) + + Returns + ------- + None + )", py::arg("pcavar"), py::arg("nsamples")=100, py::arg("ftype")=false, + py::arg("total")=true, py::arg("covar")=true); + py::class_, PyCBDisk, BasisClasses::BiorthBasis>(m, "CBDisk") .def(py::init(), @@ -2369,7 +2579,7 @@ void BasisFactoryClasses(py::module &m) py::arg("zmin")=-1.0, py::arg("zmax")=1.0, py::arg("numz")=400) - .def("orthoCheck", [](BasisClasses::Cube& A) + .def("orthoCheck", [](BasisClasses::Slab& A) { return A.orthoCheck(); }, @@ -2393,7 +2603,6 @@ void BasisFactoryClasses(py::module &m) )" ); - py::class_, PyCube, BasisClasses::BiorthBasis>(m, "Cube") .def(py::init(), R"( @@ -2412,36 +2621,6 @@ void BasisFactoryClasses(py::module &m) Cube the new instance )", py::arg("YAMLstring")) - .def("enableCoefCovariance", &BasisClasses::BiorthBasis::enableCoefCovariance, - R"( - Enable or disable the coefficient covariance computation and set the - default number of partitions to use for the covariance computation. - - Parameters - ---------- - pcavar : bool - enable (true) or disable (false) the covariance computation - nsamples : int - number of time partitions to use for covariance computation - ftype: bool - if true, use float32 for covariance storage; if false, - use float64 (default: false) - covar: bool - if true, compute and save covariance to the HDF5 file; if false, - save mean and variance vectors only (default: false) - - Returns - ------- - None - - Notes - ----- - The covariance computation for the Cube can be expensive in both time and memory - because the number of basis functions can be large. To save disk space, covariance - computation is disabled by default. The user can enable it by calling this member - function with covar=True. - )", py::arg("pcavar"), py::arg("nsamples")=100, py::arg("ftype")=false, - py::arg("covar")=false) .def("index1D", &BasisClasses::Cube::index1D, R"( Returns a flattened 1-d index into the arrays and matrices returned by the @@ -2507,7 +2686,85 @@ void BasisFactoryClasses(py::module &m) numpy.ndarray) list of numpy.ndarrays from [0, ... , dx*dy*dz] )" - ); + ) + .def("getCoefCovariance", + [](BasisClasses::Cube& A, double time) + { + auto [cnts, mass, coef, covr] = A.getCoefCovariance(time); + py::array_t> cf = make_ndarray3>(coef); + py::array_t> vr = make_ndarray4>(covr); + return std::make_tuple(cnts, mass, cf, vr); + }, + R"( + Get the covariance matrices for the basis coefficients + + Parameters + ---------- + time : float + the evaluation time + + Returns + ------- + tuple(numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray) + tuple of counts, masses, partitioned coefficients and their covariance + matrices for each subsample. The returns are complex-valued. + )", + py::arg("time")) + .def("writeCoefCovariance", &BasisClasses::Cube::writeCoefCovariance, + R"( + Write the partitioned coefficient vectors and covariance matrices + to a specified HDF5 file. The number of partitions is set by the + configuration parameter 'sampT' and defaults to 100. + + On first call, the file is created, metadata is written, and the + coefficient vectors and covariance matrices are stored. On subsequent + calls, the file is updated with new covariance datasets. + + The file will be called 'coefcovar...h5'. + + Parameters + ---------- + compname : str + the component/basis name segment used in the output HDF5 filename + runtag : str + the run identifier tag + time : float + the snapshot time + + Returns + ------- + None + + See also + -------- + getCoefCovariance : get the counts, mass, coefficient vectors and covariance matrices + )", py::arg("compname"), py::arg("runtag"), py::arg("time")=0.0) + .def("enableCoefCovariance", &BasisClasses::Cube::enableCoefCovariance, + R"( + Enable or disable the coefficient covariance computation and set the + default number of partitions to use for the covariance computation. + + Parameters + ---------- + pcavar : bool + enable (true) or disable (false) the covariance computation + nsamples : int + number of time partitions to use for covariance computation + ftype: bool + if true, use float32 for covariance storage; if false, + use float64 (default: false) + total: bool + if true, also compute the total covariance matrix; if false, save only + the partitioned covariance matrices (default: true) + covar: bool + if true, compute and save covariance to the HDF5 file; if false, + save mean and variance vectors only (default: true) + + Returns + ------- + None + )", py::arg("pcavar"), py::arg("nsamples")=100, py::arg("ftype")=false, + py::arg("total")=false, py::arg("covar")=false); py::class_, PyFieldBasis, BasisClasses::Basis>(m, "FieldBasis") From 3d4bf3f4e499ed22ec2903d2e126e3779a44e00f Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 19 Apr 2026 16:09:55 +0000 Subject: [PATCH 04/15] Fix Cube enableCoefCovariance docstring defaults Agent-Logs-Url: https://github.com/EXP-code/EXP/sessions/7396ab22-dac6-467e-9e74-f2fb2f26be91 Co-authored-by: The9Cat <25960766+The9Cat@users.noreply.github.com> --- pyEXP/BasisWrappers.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index a31ce4987..9447de231 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -1901,10 +1901,10 @@ void BasisFactoryClasses(py::module &m) use float64 (default: false) total: bool if true, also compute the total covariance matrix; if false, save only - the partitioned covariance matrices (default: true) + the partitioned covariance matrices (default: false) covar: bool if true, compute and save covariance to the HDF5 file; if false, - save mean and variance vectors only (default: true) + save mean and variance vectors only (default: false) Returns ------- From 8cb217fdbf52a40e161b80261e4b8e65424a6904 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 19 Apr 2026 12:15:25 -0400 Subject: [PATCH 05/15] Minor documenation updates --- expui/BasisFactory.H | 2 +- pyEXP/BasisWrappers.cc | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/expui/BasisFactory.H b/expui/BasisFactory.H index 35ac359f0..e270fae1e 100644 --- a/expui/BasisFactory.H +++ b/expui/BasisFactory.H @@ -354,7 +354,7 @@ namespace BasisClasses if (covarStore) { covarStore->setCovarH5Compress(level, chunksize, shuffle, szip); } else { - throw std::runtime_error("BiorthBasis::setCovarH5Compress: covariance storage not initialized"); + throw std::runtime_error("Basis::setCovarH5Compress: covariance storage not initialized"); } } diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index a31ce4987..b79f15244 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -1404,7 +1404,7 @@ void BasisFactoryClasses(py::module &m) Returns ------- None - )", py::arg("compress")=5, py::arg("chunkSize")=1024*1024, py::arg("shuffle")=true, py::arg("azip")=false) + )", py::arg("compress")=5, py::arg("chunkSize")=1024*1024, py::arg("shuffle")=true, py::arg("szip")=false) .def("makeFromFunction", &BasisClasses::BiorthBasis::makeFromFunction, py::call_guard(), R"( @@ -2755,10 +2755,10 @@ void BasisFactoryClasses(py::module &m) use float64 (default: false) total: bool if true, also compute the total covariance matrix; if false, save only - the partitioned covariance matrices (default: true) + the partitioned covariance matrices (default: false) covar: bool if true, compute and save covariance to the HDF5 file; if false, - save mean and variance vectors only (default: true) + save mean and variance vectors only (default: false) Returns ------- From ac35cf78c956cebff7bab05fa12dabffa2870c3f Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Apr 2026 12:52:48 -0400 Subject: [PATCH 06/15] Update pyEXP/BasisWrappers.cc Copy tensor contents into a new numpy-owned buffer. Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- pyEXP/BasisWrappers.cc | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index a2dd0a500..9a22c7b41 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -1837,8 +1837,10 @@ void BasisFactoryClasses(py::module &m) [](BasisClasses::Cylindrical& A, double time) { auto [cnts, mass, coef, covr] = A.getCoefCovariance(time); - py::array_t> cf = make_ndarray3>(coef); - py::array_t> vr = make_ndarray4>(covr); + py::array_t> cf = + make_ndarray3>(coef).attr("copy")().cast>>(); + py::array_t> vr = + make_ndarray4>(covr).attr("copy")().cast>>(); return std::make_tuple(cnts, mass, cf, vr); }, R"( From ab46047c71af362d94725d8cbca22729d4cf9ae8 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Apr 2026 12:53:59 -0400 Subject: [PATCH 07/15] Update pyEXP/BasisWrappers.cc copying tensor contents into a new numpy-owned buffer Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- pyEXP/BasisWrappers.cc | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 9a22c7b41..faae819d2 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -2359,8 +2359,11 @@ void BasisFactoryClasses(py::module &m) [](BasisClasses::FlatDisk& A, double time) { auto [cnts, mass, coef, covr] = A.getCoefCovariance(time); - py::array_t> cf = make_ndarray3>(coef); - py::array_t> vr = make_ndarray4>(covr); + py::object numpy = py::module_::import("numpy"); + py::array cf = numpy.attr("array")(make_ndarray3>(coef), + py::arg("copy") = true); + py::array vr = numpy.attr("array")(make_ndarray4>(covr), + py::arg("copy") = true); return std::make_tuple(cnts, mass, cf, vr); }, R"( From f4a313bda036cba9dd39ed71f8bd32586b9d8c99 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Apr 2026 12:54:49 -0400 Subject: [PATCH 08/15] Update pyEXP/BasisWrappers.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- pyEXP/BasisWrappers.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index faae819d2..a6334ecb9 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -327,7 +327,7 @@ void BasisFactoryClasses(py::module &m) } virtual void - enableCoefCovariance(bool pcavar, int nsamples, bool ftype, bool total, bool covar) { + enableCoefCovariance(bool pcavar, int nsamples, bool ftype, bool total, bool covar) override { PYBIND11_OVERRIDE(void, Basis, enableCoefCovariance, pcavar, nsamples, ftype, total, covar); } From fafc6b0a993a112b1e2060588cc2c6cc9f08ad38 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 19 Apr 2026 16:55:01 +0000 Subject: [PATCH 09/15] Fix covariance ndarray lifetime in pybind getCoefCovariance Agent-Logs-Url: https://github.com/EXP-code/EXP/sessions/cd11a5d9-ffd0-441c-bb19-d13d99a616d9 Co-authored-by: The9Cat <25960766+The9Cat@users.noreply.github.com> --- pyEXP/BasisWrappers.cc | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index a6334ecb9..f58936f9c 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -2151,8 +2151,10 @@ void BasisFactoryClasses(py::module &m) [](BasisClasses::SphericalSL& A, double time) { auto [cnts, mass, coef, covr] = A.getCoefCovariance(time); - py::array_t> cf = make_ndarray3>(coef); - py::array_t> vr = make_ndarray4>(covr); + py::array_t> cf = + make_ndarray3>(coef).attr("copy")().cast>>(); + py::array_t> vr = + make_ndarray4>(covr).attr("copy")().cast>>(); return std::make_tuple(cnts, mass, cf, vr); }, R"( @@ -2696,8 +2698,10 @@ void BasisFactoryClasses(py::module &m) [](BasisClasses::Cube& A, double time) { auto [cnts, mass, coef, covr] = A.getCoefCovariance(time); - py::array_t> cf = make_ndarray3>(coef); - py::array_t> vr = make_ndarray4>(covr); + py::array_t> cf = + make_ndarray3>(coef).attr("copy")().cast>>(); + py::array_t> vr = + make_ndarray4>(covr).attr("copy")().cast>>(); return std::make_tuple(cnts, mass, cf, vr); }, R"( @@ -3189,8 +3193,10 @@ void BasisFactoryClasses(py::module &m) [](BasisClasses::SubsampleCovariance& A, double time) { auto [cnts, mass, coef, covr] = A.getCoefCovariance(time); - py::array_t> cf = make_ndarray3>(coef); - py::array_t> vr = make_ndarray4>(covr); + py::array_t> cf = + make_ndarray3>(coef).attr("copy")().cast>>(); + py::array_t> vr = + make_ndarray4>(covr).attr("copy")().cast>>(); return std::make_tuple(cnts, mass, cf, vr); }, R"( From d8ee9c9f438ee9141162a1e23ca48b8dcfce5244 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Apr 2026 12:57:05 -0400 Subject: [PATCH 10/15] Update pyEXP/BasisWrappers.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- pyEXP/BasisWrappers.cc | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index f58936f9c..7c32ebc9d 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -314,7 +314,10 @@ void BasisFactoryClasses(py::module &m) virtual SubsampleCovariance::CovarData getCoefCovariance(double time) override { - PYBIND11_OVERRIDE(SubsampleCovariance::CovarData, Basis, getCoefCovariance, time); + // Do not use PYBIND11_OVERRIDE here: SubsampleCovariance::CovarData + // contains tensor-based data and this file does not define a pybind11 + // caster for that return type. Fall back to the C++ implementation. + return Basis::getCoefCovariance(time); } virtual void writeCoefCovariance From 472eca2bc2e4a08628c76841d3e4abd24c9bc8ad Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 19 Apr 2026 16:58:38 +0000 Subject: [PATCH 11/15] Normalize covariance numpy-copy pattern across bindings Agent-Logs-Url: https://github.com/EXP-code/EXP/sessions/cd11a5d9-ffd0-441c-bb19-d13d99a616d9 Co-authored-by: The9Cat <25960766+The9Cat@users.noreply.github.com> --- pyEXP/BasisWrappers.cc | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 7c32ebc9d..902ccea71 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -2364,11 +2364,10 @@ void BasisFactoryClasses(py::module &m) [](BasisClasses::FlatDisk& A, double time) { auto [cnts, mass, coef, covr] = A.getCoefCovariance(time); - py::object numpy = py::module_::import("numpy"); - py::array cf = numpy.attr("array")(make_ndarray3>(coef), - py::arg("copy") = true); - py::array vr = numpy.attr("array")(make_ndarray4>(covr), - py::arg("copy") = true); + py::array_t> cf = + make_ndarray3>(coef).attr("copy")().cast>>(); + py::array_t> vr = + make_ndarray4>(covr).attr("copy")().cast>>(); return std::make_tuple(cnts, mass, cf, vr); }, R"( From 2791c39385e59926c5528bb6a266b6ddc0073288 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 19 Apr 2026 14:35:54 -0400 Subject: [PATCH 12/15] Clean up pyEXP interface and add a bunch of docstrings about covariance generation --- pyEXP/BasisWrappers.cc | 416 +++++++++++++++++++++-------------------- 1 file changed, 213 insertions(+), 203 deletions(-) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 902ccea71..da17f3e22 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -19,10 +19,10 @@ void BasisFactoryClasses(py::module &m) BasisFactory class bindings This module provides a factory class that will create biorthogonal - bases from input YAML configuration files. Each basis can then be - used to compute coefficients, provide field quantities such as - forces and, together with the FieldGenerator, surfaces and fields - for visualization. + bases from input YAML configuration files. Each basis can then be used + to compute coefficients, provide field quantities such as forces and, + together with the FieldGenerator, surfaces and fields for + visualization. Eight bases are currently implemented: @@ -31,16 +31,16 @@ void BasisFactoryClasses(py::module &m) 2. Bessel, the classic spherical biorthogonal constructed from the eigenfunctions of the spherical Laplacian; - 3. Cylindrical, created created by computing empirical orthogonal functions - over a densely sampled SphericalSL basis; + 3. Cylindrical, created created by computing empirical orthogonal + functions over a densely sampled SphericalSL basis; 4. FlatDisk, an EOF rotation of the finite Bessel basis; 5. CBDisk, the Clutton-Brock disk basis for testing; - 6. Slab, a biorthogonal basis for a slab geometry with a finite - finite vertical extent. The basis is constructed from direct - solution of the Sturm-Liouville equation. + 6. Slab, a biorthogonal basis for a slab geometry with a finite finite + vertical extent. The basis is constructed from direct solution of + the Sturm-Liouville equation. 7. Cube, a periodic cube basis whose functions are the Cartesian eigenfunctions of the Cartesian Laplacian: sines and cosines. @@ -48,23 +48,24 @@ void BasisFactoryClasses(py::module &m) 8. FieldBasis, for computing user-provided quantities from a phase-space snapshot. - 9. VelocityBasis, for computing the mean field velocity fields from - a phase-space snapshot. This is a specialized version of FieldBasis. - - Each of these bases take a YAML configuration file as input. These parameter - lists are as subset of and have the same structure as those used by EXP. - The factory and the individual constructors will check the parameters keys - and warn of mismatches for safety. See the EXP documentation and the pyEXP - examples for more detail. The first four bases are the most often used bi- - orthogonal basis types used for computing the potential and forces from - density distributions. Other biorthogonal bases in EXP but not in pyEXP - include those for cubic and slab geometries and other special-purpose bases - such as the Hernquist, Clutton-Brock sphere and two-dimensional disk basis. - These will be made available in a future release if there is demand. Note - that the Hernquist and Clutton-Brock spheres can be constructed using + 9. VelocityBasis, for computing the mean field velocity fields from a + phase-space snapshot. This is a specialized version of FieldBasis. + + Each of these bases take a YAML configuration file as input. These + parameter lists are as subset of and have the same structure as those + used by EXP. The factory and the individual constructors will check + the parameters keys and warn of mismatches for safety. See the EXP + documentation and the pyEXP examples for more detail. The first four + bases are the most often used bi- orthogonal basis types used for + computing the potential and forces from density distributions. Other + biorthogonal bases in EXP but not in pyEXP include those for cubic and + slab geometries and other special-purpose bases such as the Hernquist, + Clutton-Brock sphere and two-dimensional disk basis. These will be + made available in a future release if there is demand. Note that the + Hernquist and Clutton-Brock spheres can be constructed using SphericalSL with a Hernquist of modified Plummer model as input. The - FieldBasis and VelocityBasis are designed for producing summary data for - post-production analysis (using mSSA or eDMD, for example) and for + FieldBasis and VelocityBasis are designed for producing summary data + for post-production analysis (using mSSA or eDMD, for example) and for simulation cross-comparison. The primary functions of these basis classes are: @@ -82,24 +83,24 @@ void BasisFactoryClasses(py::module &m) Introspection ------------- The first two bases have a 'cacheInfo(str)' member that reports the - parameters used to create the cached basis. This may be used to - grab the parameters for creating a basis. Cache use ensures that - your analyses are computed with the same bases used in a simulation - or with the same basis used on previous pyEXP invocations. At this - point, you must create the YAML configuration for the basis even if - the basis is cached. This is a safety and consistency feature that - may be relaxed in a future version. + parameters used to create the cached basis. This may be used to grab + the parameters for creating a basis. Cache use ensures that your + analyses are computed with the same bases used in a simulation or with + the same basis used on previous pyEXP invocations. At this point, you + must create the YAML configuration for the basis even if the basis is + cached. This is a safety and consistency feature that may be relaxed + in a future version. Coefficient creation -------------------- - The Basis class creates coefficients from phase space with two - methods: 'createFromReader()' and 'createFromArray()'. The first - uses a ParticleReader, see help(pyEXP.read), and the second uses - arrays of mass and 3d position vectors. Both methods take an - optional center vector (default: 0, 0, 0). You may also register - and an optional boolean functor used to select which particles to - using the 'setSelector(functor)' member. An example functor - would be defined in Python as follows: + The Basis class creates coefficients from phase space with two methods: + 'createFromReader()' and 'createFromArray()'. The first uses a + ParticleReader, see help(pyEXP.read), and the second uses arrays of + mass and 3d position vectors. Both methods take an optional center + vector (default: 0, 0, 0). You may also register and an optional + boolean functor used to select which particles to using the + 'setSelector(functor)' member. An example functor would be defined in + Python as follows: def myFunctor(m, pos, vel, index): ret = False # Default return value @@ -107,14 +108,14 @@ void BasisFactoryClasses(py::module &m) # integer index that sets ret to True if desired . . . return ret - If you are using 'createFromArray()', you will only have access to - the mass and position vector. You may clear and turn off the - selector using the 'clrSelector()' member. + If you are using 'createFromArray()', you will only have access to the + mass and position vector. You may clear and turn off the selector + using the 'clrSelector()' member. The FieldBasis class requires a user-specified phase-space field - functor that produces an list of quantities derived from the - phase space for each particle. For example, to get a total - velocity field, we could use: + functor that produces an list of quantities derived from the phase + space for each particle. For example, to get a total velocity field, + we could use: def totalVelocity(m, pos, vel): # Some caculation with scalar mass, pos array, vel array. @@ -123,41 +124,100 @@ void BasisFactoryClasses(py::module &m) This function is registered with the FieldBasis using: - basis->addPSFunction(totalVelocity, ['total velocity']) + basis.addPSFunction(totalVelocity, ['total velocity']) The VelocityBasis is a FieldBasis that automatically sets the phase-space field functor to cylindrical or spherical velocities based on the 'dof' parameter. More on 'dof' below. + Covariance creation + ------------------- + + The SphericalSL, FlatDisk, Cylindrical and Cube bases have built-in + support for computing the coefficieint covariance from subsamples of + particles. This is implemented by the enableCoefCovariance() method + for each of these supported bases. The force configuration must + contain the parameters 'pcavar' (boolean) and 'subsamp' (integer) keys. + The 'pcavar' parameter turns on the covariance computation. The + 'subsamp' parameter sets the number of partions or subsamples for each + coefficient creation. There are two additional control parameters that + may be optionally specified with the enableCoefCovariance() call. The + 'total' parameters enables computing the total covariance matrices only + rather than the covariance for each subsample which costs more + memory. The 'covar' parameter toggles the computation of the + off-diagonal terms covariance. Set to 'False' to save internal memory + and disk space. When enabled, the covariance data is stored in an HDF5 + file with a name based on the component name and run tag. The + covariance data must be written to the file at each time step with the + writeCoefCovariance call(). + + Typical usage might be: + + basis = pyEXP.basis.Basis.factory('config.yml')) + basis.enableCoefCovariance(True, 100) + . + . + . + # Loop over snapshots + for group in batches: + reader = pyEXP.read.ParticleReader.createReader('PSPhdf5', group) + reader.SelectType('dark') + coefs = basis.createFromReader(reader) + basis.writeCoefCovariance('dark', 'myrun', reader.CurrentTime()) + + Each call to writeCoefCovariance() writes the covariance data for the + current snapshot to the HDF5 file. + + The covariance data is not available to the user until it is read from + the HDF5 file. This is done by creating a CovarianceReader instance. + Once the covariance data is read, it is available to the user through + the getCoefCovariance() method of the CovarianceReader instance. The + covariance data is returned as a tuple of four elements: the sample + counts, the sample masses, the coefficient means and the coefficient + covariance matrices. The first two elements are vectors of length + equal to the number of subsamples. The third element is a 4D array + with dimensions (subsamp, Nlm, nmax) where subsamp is the number of + subsamples, Nlm, is the number of harmonics (e.g. l, m values or m + values for the spherical and cylindrical bases), and nmax is the number + of coefficients. The fourth element is a 4D array with dimensions + (subsamp, Nlm, nmax, namx) containing the covariance matrices for each + subsample in the last two dimensions. + + Typical usage might be: + covarReader = pyEXP.covar.CovarianceReader('dark', 'myrun') + covarData = covarReader.getCoefCovariance(time) + + Scalability ------------ + createFromArray() is a convenience method allows you to transform coordinates and preprocess phase space using your own methods and readers. Inside this method are three member functions calls that separately initialize, accumulate the coefficient contributions from the provided vectors, and finally construct and return the new coeffi- - cient instance (Coefs). For scalability, we provide access to each - of these three methods so that the phase space may be partitioned into - any number of smaller pieces. These three members are: initFromArray(), + cient instance (Coefs). For scalability, we provide access to each of + these three methods so that the phase space may be partitioned into any + number of smaller pieces. These three members are: initFromArray(), addFromArray(), makeFromArray(). The initFromArray() is called once to begin the creation and the makeFromArray() method is called once to build the final set of coefficients. The addFromArray() may be called any number of times in between. For example, the addFromArray() call can be inside of a loop that iterates over any partition of phase space from your own pipeline. The underlying computation is identical to - createFromArray(). However, access to the three underlying steps allows - you to scale your phase-space processing to snapshots of any size. - For reference, the createFromReader() method uses a producer-consumer - pattern internally to provide scalability. These three methods allow - you to provide the same pattern in your own pipeline. Finally, the - makeFromFunction() creates coefficients from a user-supplied - density or potential field function. + createFromArray(). However, access to the three underlying steps + allows you to scale your phase-space processing to snapshots of any + size. For reference, the createFromReader() method uses a + producer-consumer pattern internally to provide scalability. These + three methods allow you to provide the same pattern in your own + pipeline. Finally, the makeFromFunction() creates coefficients from a + user-supplied density or potential field function. Coordinate systems ------------------- - Each basis is assigned a natural coordinate system for field evaluation - as follows: + Each basis is assigned a natural coordinate system for field + evaluation as follows: 1. SphericalSL uses spherical coordinates @@ -177,45 +237,46 @@ void BasisFactoryClasses(py::module &m) the 'dof' parameter. These use cylindrical and spherical coordinates, respectively, by default. - These default choices may be overridden by passing a string argument - to the 'setFieldType()' member. The argument is case insensitive and only - distinguishing characters are necessary. E.g. for 'Cylindrical', the - argument 'cyl' or even 'cy' is sufficient. The argument 'c' is clearly + These default choices may be overridden by passing a string argument to + the 'setFieldType()' member. The argument is case insensitive and only + distinguishing characters are necessary. E.g. for 'Cylindrical', the + argument 'cyl' or even 'cy' is sufficient. The argument 'c' is clearly not enough. Orbit integration ----------------- The IntegrateOrbits routine uses a fixed time step leap frog integrator to advance orbits from tinit to tfinal with time step h. The initial - positions and velocities are supplied in an nx6 NumPy array. Tuples - of the basis (a Basis instance) and coefficient database (a Coefs + positions and velocities are supplied in an nx6 NumPy array. Tuples of + the basis (a Basis instance) and coefficient database (a Coefs instance) for each component is supplied to IntegrateOrbtis as a list. - Finally, the type of acceleration is an instance of the AccelFunc class. - The acceleration at each time step is computed by setting a coefficient - set in Basis and evaluating and accumulating the acceleration for each - phase-space point. The coefficient are handled by implementing the - evalcoefs() method of AccelFunc. We supply two implemented derived - classes, AllTimeFunc and SingleTimeFunc. The first interpolates on the - Coefs data base and installs the interpolated coefficients for the - current time in the basis instance. The SingleTimeFunc interpolates on - the Coefs data base for a single fixed time and sets the interpolated - coefficients once at the beginning of the integration. This implements - a fixed potential model. AccelFunc can be inherited by a native Python - class and the evalcoefs() may be implemented in Python and passed to - IntegrateOrbits in the same way as a native C++ class. + Finally, the type of acceleration is an instance of the AccelFunc + class. The acceleration at each time step is computed by setting a + coefficient set in Basis and evaluating and accumulating the + acceleration for each phase-space point. The coefficient are handled + by implementing the evalcoefs() method of AccelFunc. We supply two + implemented derived classes, AllTimeFunc and SingleTimeFunc. The first + interpolates on the Coefs data base and installs the interpolated + coefficients for the current time in the basis instance. The + SingleTimeFunc interpolates on the Coefs data base for a single fixed + time and sets the interpolated coefficients once at the beginning of + the integration. This implements a fixed potential model. AccelFunc + can be inherited by a native Python class and the evalcoefs() may be + implemented in Python and passed to IntegrateOrbits in the same way as + a native C++ class. Non-inertial frames of reference -------------------------------- - Each component of a multiple component simulation may have its own expansion - center. Orbit integration in the frame of reference of the expansion is - accomplished by defining a moving frame of reference using the setNonInertial() - call with either an array of n times and center positions (as an nx3 array) - or by initializing with an EXP orient file. + Each component of a multiple component simulation may have its own + expansion center. Orbit integration in the frame of reference of the + expansion is accomplished by defining a moving frame of reference using + the setNonInertial() call with either an array of n times and center + positions (as an nx3 array) or by initializing with an EXP orient file. - We provide a member function, setNonInertialAccel(t), to estimate the frame - acceleration at a given time. This may be useful for user-defined acceleration - routines. This is automatically called default C++ evalcoefs() routine. - )"; + We provide a member function, setNonInertialAccel(t), to estimate the + frame acceleration at a given time. This may be useful for + user-defined acceleration routines. This is automatically called + default C++ evalcoefs() routine. )"; using namespace BasisClasses; @@ -312,14 +373,6 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE_PURE(void, Basis, addFromArray, m, p, roundrobin, posvelrows); } - virtual SubsampleCovariance::CovarData - getCoefCovariance(double time) override { - // Do not use PYBIND11_OVERRIDE here: SubsampleCovariance::CovarData - // contains tensor-based data and this file does not define a pybind11 - // caster for that return type. Fall back to the C++ implementation. - return Basis::getCoefCovariance(time); - } - virtual void writeCoefCovariance (const std::string& file, const std::string& runtag, double time) override { PYBIND11_OVERRIDE(void, Basis, writeCoefCovariance, file, runtag, time); @@ -1836,31 +1889,6 @@ void BasisFactoryClasses(py::module &m) cache parameters )", py::arg("cachefile")) - .def("getCoefCovariance", - [](BasisClasses::Cylindrical& A, double time) - { - auto [cnts, mass, coef, covr] = A.getCoefCovariance(time); - py::array_t> cf = - make_ndarray3>(coef).attr("copy")().cast>>(); - py::array_t> vr = - make_ndarray4>(covr).attr("copy")().cast>>(); - return std::make_tuple(cnts, mass, cf, vr); - }, - R"( - Get the covariance matrices for the basis coefficients - - Parameters - ---------- - time : float - the evaluation time - - Returns - ------- - tuple(numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray) - tuple of counts, masses, partitioned coefficients and their covariance - matrices for each subsample. The returns are complex-valued. - )", - py::arg("time")) .def("writeCoefCovariance", &BasisClasses::Cylindrical::writeCoefCovariance, R"( Write the partitioned coefficient vectors and covariance matrices @@ -1886,9 +1914,21 @@ void BasisFactoryClasses(py::module &m) ------- None + Notes + ----- + The coefficient covariance computation is enabled with the enableCoefCovariance() + member function. This member function also sets the number of partitions to use for + the covariance computation and whether to save the covariance matrices or just the + mean and variance vectors. + + Each call to writeCoefCovariance() saves the covariance data for 'time' to an HDF5 + database file. This file can be read with the CovarianceReader helper class. + See also -------- - getCoefCovariance : get the counts, mass, coefficient vectors and covariance matrices + enableCoefCovariance : enable the coefficient covariance computation and set parameters + CovarianceReader : helper class for reading and accessing the covariance data stored + in the HDF5 file )", py::arg("compname"), py::arg("runtag"), py::arg("time")=0.0) .def("enableCoefCovariance", &BasisClasses::Cylindrical::enableCoefCovariance, R"( @@ -2150,31 +2190,6 @@ void BasisFactoryClasses(py::module &m) list of numpy.ndarrays from [0, ... , Lmax] )", py::arg("knots")=40) - .def("getCoefCovariance", - [](BasisClasses::SphericalSL& A, double time) - { - auto [cnts, mass, coef, covr] = A.getCoefCovariance(time); - py::array_t> cf = - make_ndarray3>(coef).attr("copy")().cast>>(); - py::array_t> vr = - make_ndarray4>(covr).attr("copy")().cast>>(); - return std::make_tuple(cnts, mass, cf, vr); - }, - R"( - Get the covariance matrices for the basis coefficients - - Parameters - ---------- - time : float - the evaluation time - - Returns - ------- - tuple(numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray) - tuple of counts, masses, partitioned coefficients and their covariance - matrices for each subsample. The returns are complex-valued. - )", - py::arg("time")) .def("writeCoefCovariance", &BasisClasses::SphericalSL::writeCoefCovariance, R"( Write the partitioned coefficient vectors and covariance matrices @@ -2182,8 +2197,9 @@ void BasisFactoryClasses(py::module &m) configuration parameter 'sampT' and defaults to 100. On first call, the file is created, metadata is written, and the - coefficient vectors and covariance matrices are stored. On subsequent - calls, the file is updated with new covariance datasets. + coefficient vectors and covariance matrices are stored. On + subsequent calls, the file is updated with new covariance + datasets. The file will be called 'coefcovar...h5'. @@ -2200,9 +2216,27 @@ void BasisFactoryClasses(py::module &m) ------- None + Notes + ----- + The coefficient covariance computation is enabled with the + enableCoefCovariance() member function. This member function also + sets the number of partitions to use for the covariance + computation and whether to save the covariance matrices or just + the mean and variance vectors. + + Each call to writeCoefCovariance() saves the covariance data for + 'time' to an HDF5 database file. This file can be read with the + CovarianceReader helper class. + See also -------- - getCoefCovariance : get the counts, mass, coefficient vectors and covariance matrices + + enableCoefCovariance : enable the coefficient covariance + computation and set parameters + + CovarianceReader : helper class for reading and accessing the + covariance data stored in the HDF5 file + )", py::arg("compname"), py::arg("runtag"), py::arg("time")=0.0) .def("enableCoefCovariance", &BasisClasses::SphericalSL::enableCoefCovariance, R"( @@ -2360,31 +2394,6 @@ void BasisFactoryClasses(py::module &m) cache parameters )", py::arg("cachefile")) - .def("getCoefCovariance", - [](BasisClasses::FlatDisk& A, double time) - { - auto [cnts, mass, coef, covr] = A.getCoefCovariance(time); - py::array_t> cf = - make_ndarray3>(coef).attr("copy")().cast>>(); - py::array_t> vr = - make_ndarray4>(covr).attr("copy")().cast>>(); - return std::make_tuple(cnts, mass, cf, vr); - }, - R"( - Get the covariance matrices for the basis coefficients - - Parameters - ---------- - time : float - the evaluation time - - Returns - ------- - tuple(numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray) - tuple of counts, masses, partitioned coefficients and their covariance - matrices for each subsample. The returns are complex-valued. - )", - py::arg("time")) .def("writeCoefCovariance", &BasisClasses::FlatDisk::writeCoefCovariance, R"( Write the partitioned coefficient vectors and covariance matrices @@ -2410,9 +2419,22 @@ void BasisFactoryClasses(py::module &m) ------- None + Notes + ----- + The coefficient covariance computation is enabled with the enableCoefCovariance() + member function. This member function also sets the number of partitions to use for + the covariance computation and whether to save the covariance matrices or just the + mean and variance vectors. + + Each call to writeCoefCovariance() saves the covariance data for 'time' to an HDF5 + database file. This file can be read with the CovarianceReader helper class. + See also -------- - getCoefCovariance : get the counts, mass, coefficient vectors and covariance matrices + enableCoefCovariance : enable the coefficient covariance computation and set parameters + CovarianceReader : helper class for reading and accessing the covariance data stored + in the HDF5 file + )", py::arg("compname"), py::arg("runtag"), py::arg("time")=0.0) .def("enableCoefCovariance", &BasisClasses::FlatDisk::enableCoefCovariance, R"( @@ -2633,7 +2655,7 @@ void BasisFactoryClasses(py::module &m) .def("index1D", &BasisClasses::Cube::index1D, R"( Returns a flattened 1-d index into the arrays and matrices returned by the - getCoefCovariance() routines from wave number indexing. + CovarianceReader::getCoefCovariance() routines from wave number indexing. Parameters ---------- @@ -2659,7 +2681,7 @@ void BasisFactoryClasses(py::module &m) .def("index3D", &BasisClasses::Cube::index3D, R"( Returns a tuple of indices for the wave-numbers (kx, ky, kz) from the flattened - 1-d index for the arrays and matrices returned by the getCoefCovariance() routines + 1-d index for the arrays and matrices returned by the SubsampleCovariance::getCoefCovariance() routines Parameters ---------- @@ -2696,31 +2718,6 @@ void BasisFactoryClasses(py::module &m) list of numpy.ndarrays from [0, ... , dx*dy*dz] )" ) - .def("getCoefCovariance", - [](BasisClasses::Cube& A, double time) - { - auto [cnts, mass, coef, covr] = A.getCoefCovariance(time); - py::array_t> cf = - make_ndarray3>(coef).attr("copy")().cast>>(); - py::array_t> vr = - make_ndarray4>(covr).attr("copy")().cast>>(); - return std::make_tuple(cnts, mass, cf, vr); - }, - R"( - Get the covariance matrices for the basis coefficients - - Parameters - ---------- - time : float - the evaluation time - - Returns - ------- - tuple(numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray) - tuple of counts, masses, partitioned coefficients and their covariance - matrices for each subsample. The returns are complex-valued. - )", - py::arg("time")) .def("writeCoefCovariance", &BasisClasses::Cube::writeCoefCovariance, R"( Write the partitioned coefficient vectors and covariance matrices @@ -2746,9 +2743,22 @@ void BasisFactoryClasses(py::module &m) ------- None + Notes + ----- + The coefficient covariance computation is enabled with the enableCoefCovariance() + member function. This member function also sets the number of partitions to use for + the covariance computation and whether to save the covariance matrices or just the + mean and variance vectors. + + Each call to writeCoefCovariance() saves the covariance data for 'time' to an HDF5 + database file. This file can be read with the CovarianceReader helper class. + See also -------- - getCoefCovariance : get the counts, mass, coefficient vectors and covariance matrices + enableCoefCovariance : enable the coefficient covariance computation and set parameters + CovarianceReader : helper class for reading and accessing the covariance data stored + in the HDF5 file + )", py::arg("compname"), py::arg("runtag"), py::arg("time")=0.0) .def("enableCoefCovariance", &BasisClasses::Cube::enableCoefCovariance, R"( From e3e8b9e81ee188a34f66e8000ee9c939887a9985 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Mon, 20 Apr 2026 11:12:16 +0100 Subject: [PATCH 13/15] add quick test --- pyEXP/BasisWrappers.cc | 2 +- tests/CMakeLists.txt | 2 +- tests/Halo/createCoefs.py | 12 +++++++++++- 3 files changed, 13 insertions(+), 3 deletions(-) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index da17f3e22..1429c6249 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -184,7 +184,7 @@ void BasisFactoryClasses(py::module &m) subsample in the last two dimensions. Typical usage might be: - covarReader = pyEXP.covar.CovarianceReader('dark', 'myrun') + covarReader = pyEXP.basis.CovarianceReader(filename) covarData = covarReader.getCoefCovariance(time) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 30dedb2d4..45225dc3e 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -110,7 +110,7 @@ if(ENABLE_NBODY) COMMAND ${CMAKE_COMMAND} -E remove config.run0.yml current.processor.rates.run0 new.bods OUTLOG.run0 run0.levels SLGridSph.cache.run0 test.grid - outcoef.halo.run0 SLGridSph.cache.run0 + outcoef.halo.run0 SLGridSph.cache.run0 coefcovar.halo.test_covar.h5 WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) # Remove the temporary files diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py index c72b5566d..5d6bb89d0 100644 --- a/tests/Halo/createCoefs.py +++ b/tests/Halo/createCoefs.py @@ -16,6 +16,8 @@ rmapping : 0.0667 modelname: SLGridSph.model cachename: SLGridSph.cache.run0 + pcavar : true + subsamp : 10 ... """ @@ -24,8 +26,9 @@ # Construct the basis instances # basis = pyEXP.basis.Basis.factory(config) +basis.enableCoefCovariance(True, 100) # new for version 7.9.3 -print("---- created basis") +print("---- created basis and enabled covariance") # Create a coefficient structure # @@ -53,6 +56,13 @@ coefs.add(coef1) +# do a covariance calculation to test that the covariance options are working : new for version 7.9.3 +basis.writeCoefCovariance('halo','test_covar', coef1.time ) + +# read back in +testcovar = pyEXP.basis.CovarianceReader('coefcovar.halo.test_covar.h5') +testcovar.getCoefCovariance(0.0) + print("Times:", coefs.Times()) print("---- creating array data from list data") From eeb802f1905b15d44270e83ca18c51e38bc0a017 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Mon, 20 Apr 2026 11:36:58 +0100 Subject: [PATCH 14/15] this is why we have tests --- tests/Halo/createCoefs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/Halo/createCoefs.py b/tests/Halo/createCoefs.py index 5d6bb89d0..67e2a455e 100644 --- a/tests/Halo/createCoefs.py +++ b/tests/Halo/createCoefs.py @@ -61,7 +61,7 @@ # read back in testcovar = pyEXP.basis.CovarianceReader('coefcovar.halo.test_covar.h5') -testcovar.getCoefCovariance(0.0) +testcovar.getCoefCovariance(coef1.time) print("Times:", coefs.Times()) From e880c6484ed0d165d4f75455f8e404640fb12200 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 20 Apr 2026 08:25:59 -0400 Subject: [PATCH 15/15] Bug fix version bump --- CMakeLists.txt | 2 +- doc/exp.cfg | 2 +- doc/exp.cfg.breathe | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 9d3dd2947..ee224bde3 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.2" + VERSION "7.9.3" HOMEPAGE_URL https://github.com/EXP-code/EXP LANGUAGES C CXX Fortran) diff --git a/doc/exp.cfg b/doc/exp.cfg index c930e7812..f76157769 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.2 +PROJECT_NUMBER = 7.9.3 # 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 96b4b616b..740430c50 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.2 +PROJECT_NUMBER = 7.9.3 # 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