From 080aa4437861560fa1983313706b9b0c6f3f746b Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 18 Oct 2025 14:50:17 -0400 Subject: [PATCH 01/82] A preliminary implementation of partition covariance --- expui/BiorthBasis.H | 53 ++++++++++++++++++++++++++ expui/BiorthBasis.cc | 88 ++++++++++++++++++++++++++++++++++++++++++-- exputil/EmpCylSL.cc | 38 +++++++++++++++++++ include/EmpCylSL.H | 9 +++++ 4 files changed, 185 insertions(+), 3 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index fac810bf3..2bc952e24 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -96,6 +96,15 @@ namespace BasisClasses RowMatrixXd varray; //@} + //! Coefficient variance computation enabled + bool pcavar = false; + + //@{ + //! Sample counts and masses for covariance computation + Eigen::VectorXi sampleCounts; + Eigen::VectorXd sampleMasses; + //@} + public: //! Constructor from YAML node @@ -213,6 +222,23 @@ namespace BasisClasses return varray; } + + //! Return a vector of tuples of basis functions and the + //! covariance matrix for subsamples of particles + using CoefCovarType = std::tuple; + virtual std::vector> + getCoefCovariance() + { + throw std::runtime_error("BiorthBasis::getCoefCovariance: Not implemented for this basis"); + } + + //! Get sample counts for the covariance computation + virtual std::tuple + getCovarSamples() + { + return std::make_tuple(sampleCounts, sampleMasses); + } + }; /** @@ -308,6 +334,16 @@ namespace BasisClasses virtual void computeAccel(double x, double y, double z, Eigen::Ref vstore); + //@{ + //! Covariance structures. First index is T, second is harmonic + //! (l,m). + std::vector> meanV; + std::vector> covrV; + void init_covariance(); + void zero_covariance(); + int sampT = 100; + //@} + public: //! Constructor from YAML node @@ -363,6 +399,10 @@ namespace BasisClasses return ret; } + //! Return a vector of tuples of basis functions and the + //! covariance matrix for subsamples of particles + std::vector> + getCoefCovariance(); }; /** @@ -920,6 +960,19 @@ namespace BasisClasses return ret; } + //! Return a vector of tuples of basis functions and the + //! covariance matrix for subsamples of particles + std::vector> + getCoefCovariance() + { + return sl->getCoefCovariance(); + } + + std::tuple + getCovarSamples() + { + return sl->getCovarSamples(); + } }; /** diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 9963284cf..a373af843 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -26,6 +26,7 @@ namespace BasisClasses "pcavtk", "pcaeof", "subsamp", + "sampT", "snr", "samplesz", "vtkfreq", @@ -262,6 +263,8 @@ namespace BasisClasses if (conf["EVEN_L"]) EVEN_L = conf["EVEN_L"].as(); if (conf["EVEN_M"]) EVEN_M = conf["EVEN_M"].as(); if (conf["M0_ONLY"]) M0_only = conf["M0_ONLY"].as(); + if (conf["pcavar"]) pcavar = conf["pcavar"].as(); + if (conf["sampT"]) sampT = conf["sampT"].as(); } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing parameter stanza for <" @@ -312,11 +315,39 @@ namespace BasisClasses used = 0; + // Initialize covariance + // + if (pcavar) init_covariance(); + // Set spherical coordindates // coordinates = Coord::Spherical; } + void Spherical::init_covariance() + { + if (pcavar) { + meanV.resize(sampT); + for (auto& v : meanV) v.resize(nmax); + covrV.resize(sampT); + for (auto& v : covrV) v.resize(nmax); + sampleCounts.resize(sampT); + sampleMasses.resize(sampT); + zero_covariance(); + } + } + + void Spherical::zero_covariance() + { + for (int T=0; T0 && expcoef.cols()>0) expcoef.setZero(); totalMass = 0.0; used = 0; + if (pcavar) zero_covariance(); } @@ -535,6 +567,14 @@ namespace BasisClasses legendre_R(lmax, costh, legs[tid]); + // Sample index for pcavar + int T = 0; + if (pcavar) { + T = used % sampT; + sampleCounts(T) += 1; + sampleMasses(T) += mass; + } + // L loop for (int l=0, loffset=0; l<=lmax; loffset+=(2*l+1), l++) { @@ -549,6 +589,16 @@ namespace BasisClasses for (int n=0; n> + Spherical::getCoefCovariance() + { + std::vector> ret; + if (pcavar) { + ret.resize(sampT); + for (int T=0; T(ret[T][l]) = meanV[T][l]; + std::get<1>(ret[T][l]) = covrV[T][l]; + } + } + } + + return ret; + } + #define MINEPS 1.0e-10 @@ -1350,9 +1431,9 @@ namespace BasisClasses if (conf["aratio" ]) aratio = conf["aratio" ].as(); if (conf["hratio" ]) hratio = conf["hratio" ].as(); - if (conf["dweight" ]) dweight = conf["dweight" ].as(); - if (conf["Mfac" ]) Mfac = conf["Mfac" ].as(); - if (conf["HERNA" ]) HERNA = conf["HERNA" ].as(); + if (conf["dweight" ]) dweight = conf["dweight" ].as(); + if (conf["Mfac" ]) Mfac = conf["Mfac" ].as(); + if (conf["HERNA" ]) HERNA = conf["HERNA" ].as(); if (conf["rwidth" ]) rwidth = conf["rwidth" ].as(); if (conf["ashift" ]) ashift = conf["ashift" ].as(); if (conf["rfactor" ]) rfactor = conf["rfactor" ].as(); @@ -1362,6 +1443,7 @@ namespace BasisClasses if (conf["dtype" ]) dtype = conf["dtype" ].as(); if (conf["vflag" ]) vflag = conf["vflag" ].as(); if (conf["pyname" ]) pyname = conf["pyname" ].as(); + if (conf["pcavar"] ) pcavar = conf["pcavar" ].as(); // Deprecation warning if (conf["density" ]) { diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 9024c0e47..383cfbae5 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -4812,6 +4812,44 @@ void EmpCylSL::pca_hall(bool compute, bool subsamp) } +/** Return a vector of tuples of basis functions and the + covariance matrix for subsamples of particles */ +std::vector> +EmpCylSL::getCoefCovariance() +{ + std::vector>> ret; + + if (PCAVAR) { + ret.resize(sampT); + for (unsigned T=0; T(ret[T][M]) = covV[0][T][M]; + std::get<1>(ret[T][M]) = covM[0][T][M]; + } + } + } + + return ret; +} + +std::tuple +EmpCylSL::getCovarSamples() +{ + std::tuple ret; + + if (PCAVAR) { + std::get<0>(ret).resize(sampT); + std::get<1>(ret).resize(sampT); + for (int T=0; T(ret)(T) = numbT[T]; + std::get<1>(ret)(T) = massT[T]; + } + } + + return ret; +} + void EmpCylSL::get_trimmed (double snr, std::vector& ac_cos, std::vector& ac_sin, diff --git a/include/EmpCylSL.H b/include/EmpCylSL.H index f745297f7..39d04a108 100644 --- a/include/EmpCylSL.H +++ b/include/EmpCylSL.H @@ -930,6 +930,15 @@ public: //! Check orthogonality for basis (pyEXP style) std::vector orthoCheck(); + //! Return a vector of tuples of basis functions and the + //! covariance matrix for subsamples of particles + using CoefCovarType = std::tuple; + std::vector> getCoefCovariance(); + + //! Covariance statistics + std::tuple getCovarSamples(); + + #if HAVE_LIBCUDA==1 cudaMappingConstants getCudaMappingConstants(); From d8d04c595bf65038e02c6f7e16ceb95fd6ff8f20 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 18 Oct 2025 22:33:30 -0400 Subject: [PATCH 02/82] Cleanup up bindings for new variance members --- pyEXP/BasisWrappers.cc | 57 ++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 55 insertions(+), 2 deletions(-) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 2a3378389..5db2e8b0c 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -435,6 +435,19 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE_PURE(void, BiorthBasis, set_coefs, coefs); } + using CCelement = std::tuple; + using CCreturn = std::vector>; + + CCreturn getCoefCovariance(void) override { + PYBIND11_OVERRIDE(CCreturn, BiorthBasis, getCoefCovariance,); + } + + using Selement = std::tuple; + + Selement getCovarSamples() override { + PYBIND11_OVERRIDE(Selement, BiorthBasis, getCovarSamples,); + } + }; class PySpherical : public Spherical @@ -1316,6 +1329,44 @@ void BasisFactoryClasses(py::module &m) makeFromArray : create coefficients contributions )", py::arg("mass"), py::arg("pos")) + .def("getCoefCovariance", &BasisClasses::BiorthBasis::getCoefCovariance, + R"( + Return the coefficient vectors and covariance matrices for a partition + of the accumulated particles. The number of partitions is set by the + configuration parameters 'sampT' and defaults to 100. + + Parameters + ---------- + None + + Returns + ------- + arrays: a 2-dimensional array of tuples of numpy.ndarray representing the + coefficient vectors and coefficient variance matrices + + See also + -------- + getCovarSamples : get counts and mass in each partition + )") + .def("getCovarSamples", &BasisClasses::BiorthBasis::getCovarSamples, + R"( + Return a vector of counts and mass used to compute the partitioned + vectors and covariance. The rank of the vectors is set by 'sampT' + + Parameters + ---------- + None + + Returns + ------- + arrays: a tuple of arrays containing the counts and mass in each + partitioned sample + + See also + -------- + getCoefCovariance : get the coefficient vectors and covariance matrices + for the partitioned phase space. + )") .def("getFields", &BasisClasses::BiorthBasis::getFields, R"( Return the field evaluations for a given cartesian position. The @@ -1616,6 +1667,7 @@ void BasisFactoryClasses(py::module &m) )", py::arg("cachefile")); + py::class_, PySpherical, BasisClasses::BiorthBasis>(m, "Spherical") .def(py::init(), R"( @@ -1702,6 +1754,7 @@ void BasisFactoryClasses(py::module &m) cache parameters )", py::arg("cachefile")) + .def_static("I", [](int l, int m) { if (l<0) throw std::runtime_error("l must be greater than 0"); @@ -1725,6 +1778,7 @@ void BasisFactoryClasses(py::module &m) index array packing index )", py::arg("l"), py::arg("m")) + .def_static("invI", [](int I) { if (I<0) std::runtime_error("I must be an interger greater than or equal to 0"); @@ -1744,8 +1798,7 @@ void BasisFactoryClasses(py::module &m) ------- (l, m) : tuple the harmonic indices (l, m). - )", py::arg("I")); - + )", py::arg("I")); py::class_, BasisClasses::Spherical>(m, "SphericalSL") .def(py::init(), From 5de90613973cc2d54dd6520e1b19283393f567c4 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 19 Oct 2025 12:31:53 -0400 Subject: [PATCH 03/82] Implementation of HDF5 writing --- expui/BiorthBasis.H | 24 ++++ expui/BiorthBasis.cc | 254 ++++++++++++++++++++++++++++++++++++++++- pyEXP/BasisWrappers.cc | 39 ++++++- 3 files changed, 314 insertions(+), 3 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 2bc952e24..0938ddcf2 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -105,8 +105,23 @@ namespace BasisClasses Eigen::VectorXd sampleMasses; //@} + //! Write H5 covariance; returns updated group count + virtual unsigned writeCovarH5(HighFive::Group& group, unsigned count, double time); + + //! Write H5 parameters for coefficient covariance + virtual void writeCovarH5Params(HighFive::File& file) {} + + //! Add coefficient covariance data to an HDF5 file + virtual void extendCoefCovariance(const std::string& filename, double time); + + //! Variance file versioning + inline static const std::string CovarianceFileVersion = "1.0"; + public: + //! Basis identifier string + std::string BasisID = "BiorthBasis"; + //! Constructor from YAML node BiorthBasis(const YAML::Node& conf, const std::string& name="BiorthBasis") : Basis(conf, name) {} @@ -239,6 +254,11 @@ namespace BasisClasses return std::make_tuple(sampleCounts, sampleMasses); } + //! Write coefficient covariance data to an HDF5 file + void writeCoefCovariance(const std::string& compname, + const std::string& runtag, + double time=0.0); + }; /** @@ -341,6 +361,7 @@ namespace BasisClasses std::vector> covrV; void init_covariance(); void zero_covariance(); + void writeCovarH5Params(HighFive::File& file); int sampT = 100; //@} @@ -901,6 +922,9 @@ namespace BasisClasses void computeAccel(double x, double y, double z, Eigen::Ref vstore); + //! Needed for covariance writing + void writeCovarH5Params(HighFive::File& file); + public: //! Constructor from YAML node diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index a373af843..2352951b4 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -11,6 +11,9 @@ #include #endif +#include +#include + namespace BasisClasses { const std::set @@ -350,6 +353,9 @@ namespace BasisClasses void SphericalSL::initialize() { + // Identifier + // + BasisID = "sphereSL"; // Assign some defaults // @@ -411,6 +417,10 @@ namespace BasisClasses void Bessel::initialize() { + // Identifier + // + BasisID = "bessel"; + try { if (conf["rnum"]) rnum = conf["rnum"].as(); @@ -598,7 +608,6 @@ namespace BasisClasses fac * potd[tid](l, o) * norm * mass; } } - } moffset++; @@ -631,6 +640,7 @@ namespace BasisClasses void Spherical::make_coefs() { + // MPI reduction of coefficients if (use_mpi) { MPI_Allreduce(MPI_IN_PLACE, &used, 1, MPI_INT, @@ -642,6 +652,27 @@ namespace BasisClasses MPI_SUM, MPI_COMM_WORLD); expcoef.row(l) = work; } + + if (pcavar) { + + MPI_Allreduce(MPI_IN_PLACE, &sampleCounts, 1, MPI_INT, + MPI_SUM, MPI_COMM_WORLD); + + MPI_Allreduce(MPI_IN_PLACE, &sampleMasses, 1, MPI_DOUBLE, + MPI_SUM, MPI_COMM_WORLD); + + for (int T=0; T::max(); nminy = std::numeric_limits::max(); nminz = std::numeric_limits::max(); @@ -4730,5 +4779,208 @@ namespace BasisClasses return {times, ret}; } + void Spherical::writeCovarH5Params(HighFive::File& file) + { + file.createAttribute("lmax", HighFive::DataSpace::From(lmax)).write(lmax); + file.createAttribute("nmax", HighFive::DataSpace::From(nmax)).write(nmax); + file.createAttribute("scale", HighFive::DataSpace::From(scale)).write(scale); + file.createAttribute("rmin", HighFive::DataSpace::From(rmin)).write(rmin); + file.createAttribute("rmax", HighFive::DataSpace::From(rmax)).write(rmax); + } + + void Cylindrical::writeCovarH5Params(HighFive::File& file) + { + file.createAttribute("mmax", HighFive::DataSpace::From(mmax)).write(mmax); + file.createAttribute("nmax", HighFive::DataSpace::From(nmax)).write(nmax); + file.createAttribute("scale", HighFive::DataSpace::From(rscl)).write(rscl); + file.createAttribute("rmin", HighFive::DataSpace::From(rmin)).write(rmin); + file.createAttribute("rmax", HighFive::DataSpace::From(rmax)).write(rmax); + file.createAttribute("acyl", HighFive::DataSpace::From(rmin)).write(acyl); + file.createAttribute("hcyl", HighFive::DataSpace::From(rmax)).write(hcyl); + } + + unsigned BiorthBasis::writeCovarH5(HighFive::Group& snaps, unsigned count, double time) + { + std::ostringstream stim; + stim << std::setw(8) << std::setfill('0') << std::right << count++; + + // Make a new group for this time + // + HighFive::Group stanza = snaps.createGroup(stim.str()); + + // Add a time attribute + // + stanza.createAttribute("Time", HighFive::DataSpace::From(time)).write(time); + + // Add the sample statistics + // + HighFive::DataSet s1data = stanza.createDataSet("sampleCounts", sampleCounts); + HighFive::DataSet s2data = stanza.createDataSet("sampleMasses", sampleMasses); + + auto covar = getCoefCovariance(); + + // Add the samples + // + for (size_t T=0; T(covar[T][0]).rows(); + + Eigen::VectorXd data(lmax*nmax); + for (size_t l=0, c=0; l(covar[T][l])(n); + } + } + + HighFive::DataSet coefdata = sample.createDataSet("coefficients", data); + + // Pack the covariance data in a upper triangular format + // + size_t diagonalSize = nmax*(nmax + 1)/2; + data.resize(lmax*diagonalSize); + for (size_t l=0, c=0; l(covar[T][l])(n1, n2); + } + } + } + + HighFive::DataSet covdata = sample.createDataSet("covariance", data); + } + // END: sample loop + + return count; + } + + void BiorthBasis::writeCoefCovariance(const std::string& compname, const std::string& runtag, double time) + { + // Check check that variance computation is on + // + if (not pcavar) { + std::cout << "BiorthBasis::writeCoefCovariance: covariance computation is disabled. " + << "Set 'pcavar: true' to enable." << std::endl; + return; + } + + // Only root process writes + // + if (myid) return; + + // Check that there is something to write + // + int totalCount = 0; + for (auto c : sampleCounts) totalCount += c; + if (totalCount==0) { + std::cout << "BiorthBasis::writeCoefCovariance: no data" << std::endl; + return; + } + + // The H5 filename + std::string fname = "coefcovar." + compname + "." + runtag + ".h5"; + + // Check if file exists? + // + try { + // Open the HDF5 file in read-only mode + HighFive::File file(fname, HighFive::File::ReadOnly); + + // Check for version string + std::string path = "CovarianceFileVersion"; + + // If exists, use extension method + if (file.exist(path)) { + extendCoefCovariance(fname, time); + return; + } + + } catch (const HighFive::Exception& err) { + // Handle HighFive specific errors (e.g., file not found) + std::cerr << "HighFive Error: " << err.what() << std::endl; + } catch (const std::exception& err) { + // Handle other general exceptions + std::cerr << "Error: " << err.what() << std::endl; + } + + // Write coefficients + // + try { + // Create a new hdf5 file + // + HighFive::File file(fname, + HighFive::File::ReadWrite | + HighFive::File::Create); + + // Write the Version string + // + file.createAttribute("CovarianceFileVersion", HighFive::DataSpace::From(CovarianceFileVersion)).write(CovarianceFileVersion); + + // We write the basis identifier string + // + file.createAttribute("BasisID", HighFive::DataSpace::From(BasisID)).write(BasisID); + + // Write the specific parameters + // + writeCovarH5Params(file); + + // Group count variable + // + unsigned count = 0; + HighFive::DataSet dataset = file.createDataSet("count", count); + + // Create a new group for coefficient snapshots + // + HighFive::Group group = file.createGroup("snapshots"); + + // Write the coefficients + // + count = writeCovarH5(group, count, time); + + // Update the count + // + dataset.write(count); + + } catch (HighFive::Exception& err) { + std::cerr << err.what() << std::endl; + } + } + + void BiorthBasis::extendCoefCovariance(const std::string& fname, double time) + { + try { + // Open an hdf5 file + // + HighFive::File file(fname, HighFive::File::ReadWrite); + + // Get the dataset + HighFive::DataSet dataset = file.getDataSet("count"); + + unsigned count; + dataset.read(count); + + HighFive::Group group = file.getGroup("snapshots"); + + // Write the coefficients + // + count = writeCovarH5(group, count, time); + + // Update the count + // + dataset.write(count); + + } catch (HighFive::Exception& err) { + std::cerr << err.what() << std::endl; + } + } + } // END namespace BasisClasses diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 5db2e8b0c..714429b62 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -1346,7 +1346,9 @@ void BasisFactoryClasses(py::module &m) See also -------- - getCovarSamples : get counts and mass in each partition + getCovarSamples : get counts and mass in each partition + writeCoefCovariance : write the coefficient vectors, covariance + matrices, and metadata to an HDF5 file. )") .def("getCovarSamples", &BasisClasses::BiorthBasis::getCovarSamples, R"( @@ -1362,11 +1364,44 @@ void BasisFactoryClasses(py::module &m) arrays: a tuple of arrays containing the counts and mass in each partitioned sample + See also + -------- + getCoefCovariance : get the coefficient vectors and covariance matrices + for the partitioned phase space. + writeCoefCovariance : write the coefficient vectors, covariance + matrices, and metadata to an HDF5 file. + )") + .def("writeCoefCovariance", &BasisClasses::BiorthBasis::writeCoefCovariance, + R"( + Write the partitioned coefficient vectors and covariance matrices + 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 output HDF5 file name + runtag : str + the run identifier tag + time : float + the snapshot time + + Returns + ------- + None + See also -------- getCoefCovariance : get the coefficient vectors and covariance matrices for the partitioned phase space. - )") + getCovarSamples : get counts and mass in each partition + )", py::arg("compname"), py::arg("runtag"), py::arg("time")=0.0) .def("getFields", &BasisClasses::BiorthBasis::getFields, R"( Return the field evaluations for a given cartesian position. The From 8db952bc6310fd1c4b8a8ee8f745ea001ef36735 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 19 Oct 2025 13:00:39 -0400 Subject: [PATCH 04/82] Fix storage assignment mistake in covariance arrays --- expui/BiorthBasis.cc | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 2352951b4..806da61be 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -330,12 +330,23 @@ namespace BasisClasses void Spherical::init_covariance() { if (pcavar) { + int Ltot = (lmax+1)*(lmax+2)/2; + meanV.resize(sampT); - for (auto& v : meanV) v.resize(nmax); + for (auto& v : meanV) { + v.resize(Ltot); + for (auto& vec : v) vec.resize(nmax); + } + covrV.resize(sampT); - for (auto& v : covrV) v.resize(nmax); + for (auto& v : covrV) { + v.resize(Ltot); + for (auto& mat : v) mat.resize(nmax, nmax); + } + sampleCounts.resize(sampT); sampleMasses.resize(sampT); + zero_covariance(); } } From 14da31abec956cee8f87123a058dfecc04c6d1aa Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 13:16:13 -0400 Subject: [PATCH 05/82] Update expui/BiorthBasis.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 806da61be..f0b7650da 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -4806,8 +4806,8 @@ namespace BasisClasses file.createAttribute("scale", HighFive::DataSpace::From(rscl)).write(rscl); file.createAttribute("rmin", HighFive::DataSpace::From(rmin)).write(rmin); file.createAttribute("rmax", HighFive::DataSpace::From(rmax)).write(rmax); - file.createAttribute("acyl", HighFive::DataSpace::From(rmin)).write(acyl); - file.createAttribute("hcyl", HighFive::DataSpace::From(rmax)).write(hcyl); + file.createAttribute("acyl", HighFive::DataSpace::From(acyl)).write(acyl); + file.createAttribute("hcyl", HighFive::DataSpace::From(hcyl)).write(hcyl); } unsigned BiorthBasis::writeCovarH5(HighFive::Group& snaps, unsigned count, double time) From 8433db5647f4dabee0d2a6e6534ecff5b74b3df0 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 13:16:49 -0400 Subject: [PATCH 06/82] 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 714429b62..c1adfdcd4 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -1386,7 +1386,7 @@ void BasisFactoryClasses(py::module &m) Parameters ---------- compname : str - the output HDF5 file name + the component/basis name segment used in the output HDF5 filename runtag : str the run identifier tag time : float From 9f5b52c4f4f0325fe24310e47fe02d9072ae8730 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 13:17:22 -0400 Subject: [PATCH 07/82] 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 c1adfdcd4..a716a90d1 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -1353,7 +1353,7 @@ void BasisFactoryClasses(py::module &m) .def("getCovarSamples", &BasisClasses::BiorthBasis::getCovarSamples, R"( Return a vector of counts and mass used to compute the partitioned - vectors and covariance. The rank of the vectors is set by 'sampT' + vectors and covariance. The length of both returned vectors equals 'sampT' Parameters ---------- From 6223a588b578122d54aaff95a912401cc9176da1 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 19 Oct 2025 13:18:13 -0400 Subject: [PATCH 08/82] Minor updates following Copilot review --- expui/BiorthBasis.H | 3 +++ expui/BiorthBasis.cc | 41 +++++++++++++++++------------------------ 2 files changed, 20 insertions(+), 24 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 0938ddcf2..1bb1ed7c2 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -8,6 +8,9 @@ #include // For 3d rectangular grids #include +#include // HDF5 interface +#include + #include #include #include diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 806da61be..d466acc2b 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -11,9 +11,6 @@ #include #endif -#include -#include - namespace BasisClasses { const std::set @@ -612,12 +609,10 @@ namespace BasisClasses expcoef(loffset+moffset, n) += fac4 * norm * mass; if (pcavar) { - for (int n=0; n Date: Sun, 19 Oct 2025 13:46:40 -0400 Subject: [PATCH 09/82] Update expui/BiorthBasis.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 5d30c2f93..2ad2fa9d4 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -4847,7 +4847,7 @@ namespace BasisClasses HighFive::DataSet coefdata = sample.createDataSet("coefficients", data); - // Pack the covariance data in a upper triangular format + // Pack the covariance data in an upper triangular format // size_t diagonalSize = nmax*(nmax + 1)/2; data.resize(lmax*diagonalSize); From 239aa5e5eaff096a9d3b87a4443a37b8196d48af Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 19 Oct 2025 13:47:25 -0400 Subject: [PATCH 10/82] Change square to triangular allocation for storage efficiency --- expui/BiorthBasis.cc | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 5d30c2f93..0b902c7db 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -298,12 +298,12 @@ namespace BasisClasses for (auto & v : dlegs ) v.resize(lmax+1, lmax+1); for (auto & v : d2legs) v.resize(lmax+1, lmax+1); - expcoef.resize((lmax+1)*(lmax+1), nmax); + expcoef.resize((lmax+1)*(lmax+2)/2, nmax); expcoef.setZero(); work.resize(nmax); - factorial.resize(lmax+1, lmax+1); + factorial.resize(lmax+1, lmax+1); // Wasteful but simple for (int l=0; l<=lmax; l++) { for (int m=0; m<=l; m++) { @@ -476,7 +476,7 @@ namespace BasisClasses int ldim = (lmax+1)*(lmax+2)/2; // Allocate the coefficient storage - cf->store.resize((lmax+1)*(lmax+2)/2*nmax); + cf->store.resize(ldim*nmax); // Make the coefficient map cf->coefs = std::make_shared @@ -609,9 +609,9 @@ namespace BasisClasses expcoef(loffset+moffset, n) += fac4 * norm * mass; if (pcavar) { - meanV[T][l](n) += fac4 * norm * mass; + meanV[T][loffset+moffset](n) += fac4 * norm * mass; for (int o=0; o Date: Sun, 19 Oct 2025 14:42:36 -0400 Subject: [PATCH 11/82] More dimensioning consistency fixes --- expui/BiorthBasis.cc | 37 +++++++++++++++++++++++++------------ 1 file changed, 25 insertions(+), 12 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 6fe080c7f..bebf52619 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -294,16 +294,24 @@ namespace BasisClasses for (auto & v : dpt2) v.resize(lmax+1, nmax); for (auto & v : dend) v.resize(lmax+1, nmax); + // Wasteful but simple factorial table. Could be done with a + // triangular indexing scheme... + // for (auto & v : legs ) v.resize(lmax+1, lmax+1); for (auto & v : dlegs ) v.resize(lmax+1, lmax+1); for (auto & v : d2legs) v.resize(lmax+1, lmax+1); - expcoef.resize((lmax+1)*(lmax+2)/2, nmax); + // Allocate coefficient storage; stores real and imaginary parts as reals + // + expcoef.resize((lmax+1)*(lmax+1), nmax); expcoef.setZero(); work.resize(nmax); - factorial.resize(lmax+1, lmax+1); // Wasteful but simple + // Wasteful but simple factorial table. Could be done with a + // triangular indexing scheme... + // + factorial.resize(lmax+1, lmax+1); for (int l=0; l<=lmax; l++) { for (int m=0; m<=l; m++) { @@ -327,6 +335,7 @@ namespace BasisClasses void Spherical::init_covariance() { if (pcavar) { + // Triangular for l,m>=0 int Ltot = (lmax+1)*(lmax+2)/2; meanV.resize(sampT); @@ -472,7 +481,7 @@ namespace BasisClasses cf->time = time; cf->normed = true; - // Angular storage dimension + // Angular storage dimension; triangular number for complex coefs int ldim = (lmax+1)*(lmax+2)/2; // Allocate the coefficient storage @@ -594,13 +603,13 @@ namespace BasisClasses } // L loop - for (int l=0, loffset=0; l<=lmax; loffset+=(2*l+1), l++) { + for (int l=0, loffset=0, L=0; l<=lmax; loffset+=(2*l+1), l++) { Eigen::VectorXd workE; int esize = (l+1)*nmax; // M loop - for (int m=0, moffset=0, moffE=0; m<=l; m++) { + for (int m=0, moffset=0, moffE=0; m<=l; m++, L++) { if (m==0) { fac = factorial(l, m) * legs[tid](l, m); @@ -609,9 +618,9 @@ namespace BasisClasses expcoef(loffset+moffset, n) += fac4 * norm * mass; if (pcavar) { - meanV[T][loffset+moffset](n) += fac4 * norm * mass; + meanV[T][L](n) += fac4 * norm * mass; for (int o=0; o Date: Sun, 19 Oct 2025 16:57:48 -0400 Subject: [PATCH 12/82] Implementation of the CovarianceReader --- expui/BiorthBasis.H | 63 ++++++++++++++++++++++ expui/BiorthBasis.cc | 119 +++++++++++++++++++++++++++++++++++++++++ pyEXP/BasisWrappers.cc | 109 +++++++++++++++++++++++++++++++++++++ 3 files changed, 291 insertions(+) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 1bb1ed7c2..494fdba6b 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -1285,6 +1285,69 @@ namespace BasisClasses }; + class CovarianceReader + { + public: + + using CoefCovarType = std::tuple; + + protected: + + std::vector times; + std::vector sampleCounts; + std::vector sampleMasses; + std::vector>> covarData; + + const double fixedPointPrecision = 1.0e04; + std::map timeMap; + unsigned timeIndex(double time) + { + int itime = static_cast(time * fixedPointPrecision + 0.5); + auto it = timeMap.find(itime); + if (it == timeMap.end()) + throw std::runtime_error("CovarianceReader: time not found"); + return it->second; + } + + std::string basisID; + + public: + + //! Constructor from HDF5 file + CovarianceReader(const std::string& filename, int stride=1); + + //! Get time list + std::vector Times() { return times; } + + //@{ + //! Return a vector of tuples of basis functions and the + //! covariance matrix for subsamples of particles + std::vector> + getCoefCovariance(unsigned index) { return covarData[index]; } + + std::vector> + getCoefCovariance(double time) { return covarData[timeIndex(time)]; } + //@} + + + //@{ + //! Get sample counts for the covariance computation + std::tuple + getCovarSamples(unsigned index) + { + return {sampleCounts[index], sampleMasses[index]}; + } + + std::tuple + getCovarSamples(double time) + { + unsigned index = timeIndex(time); + return {sampleCounts[index], sampleMasses[index]}; + } + //@} + + std::string basisIDname() { return basisID; } + }; //! Time-dependent potential-density model using BasisCoef = std::tuple, std::shared_ptr>; diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index bebf52619..1d0dcada6 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -4,6 +4,7 @@ #include #include #include +#include #include #include @@ -5001,5 +5002,123 @@ namespace BasisClasses } } + // Read covariance data + CovarianceReader::CovarianceReader(const std::string& filename, int stride) + { + try { + // Create a new hdf5 file + // + HighFive::File file(filename, HighFive::File::ReadOnly); + + // Write the Version string + // + std::string version; + file.getAttribute("CovarianceFileVersion").read(version); + if (version != std::string("1.0")) { + throw std::runtime_error("CovarianceReader: unsupported file version, " + + version); + } + + // Read the basis identifier string + // + file.getAttribute("BasisID").read(basisID); + + int lmax, nmax, ltot; + + if (basisID == "Spherical") { + file.getAttribute("lmax").read(lmax); + file.getAttribute("nmax").read(nmax); + ltot = (lmax+1)*(lmax+2)/2; + } else if (basisID == "Cylindrical") { + file.getAttribute("mmax").read(lmax); + file.getAttribute("nmax").read(nmax); + ltot = lmax + 1; + } else { + throw std::runtime_error("CovarianceReader: unknown basis type, " + basisID); + } + + // Group count variable + // + unsigned count = 0; + file.getDataSet("count").read(count); + + // Open the snapshot group + // + auto snaps = file.getGroup("snapshots"); + + for (unsigned n=0; n(Time * fixedPointPrecision + 0.5); + timeMap[itime] = times.size(); + times.push_back(Time); + + // Get sample properties + // + sampleCounts.push_back(Eigen::VectorXi()); + file.getDataSet("sampleCounts").read(sampleCounts.back()); + + sampleMasses.push_back(Eigen::VectorXd()); + file.getDataSet("sampleMasses").read(sampleMasses.back()); + + size_t nT = sampleCounts.back().size(); + for (int T=0; T elem(ltot); + for (auto & e : elem) { + std::get<0>(e).resize(nmax); + std::get<1>(e).resize(nmax, nmax); + } + + // Get the flattened coefficient array + data = sample.getDataSet("coefficients").read(); + + // Pack the coefficient data + for (size_t l=0, c=0; l(elem[l])(n) = data(c); + } + } + + // Get the flattened covariance array + data = sample.getDataSet("covariance").read(); + + // Pack the coefficient data + for (size_t l=0, c=0; l(elem[l])(n1, n2) = data(c); + if (n1 != n2) + std::get<1>(elem[l])(n1, n2) = std::get<1>(elem[l])(n2, n1); + } + } + } + } + // END: sample loop + + } + // END: snapshot loop + } catch (HighFive::Exception& err) { + std::cerr << err.what() << std::endl; + } + } + } // END namespace BasisClasses diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index a716a90d1..971e46477 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -2628,4 +2628,113 @@ void BasisFactoryClasses(py::module &m) py::arg("tinit"), py::arg("tfinal"), py::arg("h"), py::arg("ps"), py::arg("basiscoef"), py::arg("func"), py::arg("nout")=0); + + py::class_> + (m, "CovarianceReader") + .def(py::init(), + R"( + Read a covariance database from an HDF5 file + + Parameters + ---------- + filename : str + The HDF5 filename + stride : int, default=1 + The stride for reading datasets + + Returns + ------- + CovarianceRearder + the new instance + + )", py::arg("filename"), py::arg("stride")=1) + .def("Times", &BasisClasses::CovarianceReader::Times, + R"( + Get the list of evaluation times + + Parameters + ---------- + None + + Returns + ------- + times : list(float) + a list of evaluation times + )") + .def("getCoefCovariance", static_cast>> + (CovarianceReader::*)(unsigned)>(&BasisClasses::CovarianceReader::getCoefCovariance), + R"( + Get the covariance matrices for the basis coefficients + + Parameters + ---------- + index : int + the time index + Returns + ------- + list(list(tuple(numpy.ndarray, numpy.ndarray))) + list of covariance matrices for each subsample + )", + py::arg("index")=0) + .def("getCoefCovariance", static_cast>> + (CovarianceReader::*)(double)>(&BasisClasses::CovarianceReader::getCoefCovariance), + R"( + Get the covariance matrices for the basis coefficients + + Parameters + ---------- + time : float + the evaluation time + + Returns + ------- + list(list(tuple(numpy.ndarray, numpy.ndarray))) + list of covariance matrices for each subsample + )", + py::arg("time")=0.0) + .def("getCovarSamples", static_cast + (CovarianceReader::*)(unsigned)>(&BasisClasses::CovarianceReader::getCovarSamples), + R"( + Get sample counts for the covariance computation + + Parameters + ---------- + index : unsigned int + the time index + + Returns + ------- + tuple(numpy.ndarray, numpy.ndarray) + sample counts and masses for the covariance computation + )", + py::arg("index")=0) + .def("getCovarSamples", static_cast + (CovarianceReader::*)(double)>(&BasisClasses::CovarianceReader::getCovarSamples), + R"( + Get sample counts for the covariance computation + + Parameters + ---------- + time : float + the evaluation time + + Returns + ------- + tuple(numpy.ndarray, numpy.ndarray) + sample counts and masses for the covariance computation + )", + py::arg("time")=0) + .def("basisIDname", &BasisClasses::CovarianceReader::basisIDname, + R"( + Get the basis ID name + + Parameters + ---------- + None + + Returns + ------- + BasisID : str + the basis ID name + )"); } From fef0ebde9ad9745caa00176cde12729d9e0158e9 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 16:59:36 -0400 Subject: [PATCH 13/82] Update pyEXP/BasisWrappers.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- pyEXP/BasisWrappers.cc | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index a716a90d1..36dda4094 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -1341,8 +1341,9 @@ void BasisFactoryClasses(py::module &m) Returns ------- - arrays: a 2-dimensional array of tuples of numpy.ndarray representing the - coefficient vectors and coefficient variance matrices + arrays: list[list[tuple[np.ndarray, np.ndarray]]] + Each element is a tuple (coef_vector, coef_covariance_matrix), + where coef_vector and coef_covariance_matrix are numpy.ndarray. See also -------- From b82d305e3150a61f37e9118194daa61ab6088848 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 17:00:12 -0400 Subject: [PATCH 14/82] Update expui/BiorthBasis.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index bebf52619..a8763f60d 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -671,7 +671,7 @@ namespace BasisClasses if (pcavar) { - // Triangular number of angular coefficients in complex from + // Triangular number of angular coefficients in complex form int ltot = (lmax+1)*(lmax+2)/2; MPI_Allreduce(MPI_IN_PLACE, sampleCounts.data(), sampleCounts.size(), From d16eddf95481daf3460f1a0a92283ed61b5597a5 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 19 Oct 2025 17:01:54 -0400 Subject: [PATCH 15/82] Some docstring corrections --- pyEXP/BasisWrappers.cc | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 971e46477..4e5974128 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -1361,8 +1361,9 @@ void BasisFactoryClasses(py::module &m) Returns ------- - arrays: a tuple of arrays containing the counts and mass in each - partitioned sample + arrays: tuple(ndarray, ndarray) + a tuple of arrays containing the counts and mass in each + partitioned sample See also -------- From e4833ff79247f934bd6ff176b698c31edcf3aa40 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 17:13:37 -0400 Subject: [PATCH 16/82] 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 0833d6cbe..dceefce02 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -2646,7 +2646,7 @@ void BasisFactoryClasses(py::module &m) Returns ------- - CovarianceRearder + CovarianceReader the new instance )", py::arg("filename"), py::arg("stride")=1) From 8b59de36dbff55efbd83070144aaac7da18acae9 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 17:16:10 -0400 Subject: [PATCH 17/82] Update expui/BiorthBasis.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index a0afe7bc8..6b697b6b4 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -5103,10 +5103,10 @@ namespace BasisClasses // Pack the coefficient data for (size_t l=0, c=0; l(elem[l])(n1, n2) = data(c); if (n1 != n2) - std::get<1>(elem[l])(n1, n2) = std::get<1>(elem[l])(n2, n1); + std::get<1>(elem[l])(n2, n1) = std::get<1>(elem[l])(n1, n2); } } } From 0cc2ba5aa793c69069f7ea892d4dfe7182187645 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 19 Oct 2025 17:16:58 -0400 Subject: [PATCH 18/82] More minor fixes --- expui/BiorthBasis.cc | 4 ++-- pyEXP/BasisWrappers.cc | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index a0afe7bc8..85b34db71 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -5063,10 +5063,10 @@ namespace BasisClasses // Get sample properties // sampleCounts.push_back(Eigen::VectorXi()); - file.getDataSet("sampleCounts").read(sampleCounts.back()); + stanza.getDataSet("sampleCounts").read(sampleCounts.back()); sampleMasses.push_back(Eigen::VectorXd()); - file.getDataSet("sampleMasses").read(sampleMasses.back()); + stanza.getDataSet("sampleMasses").read(sampleMasses.back()); size_t nT = sampleCounts.back().size(); for (int T=0; T>> - (CovarianceReader::*)(unsigned)>(&BasisClasses::CovarianceReader::getCoefCovariance), + (BasisClasses::CovarianceReader::*)(unsigned)>(&BasisClasses::CovarianceReader::getCoefCovariance), R"( Get the covariance matrices for the basis coefficients @@ -2679,7 +2679,7 @@ void BasisFactoryClasses(py::module &m) )", py::arg("index")=0) .def("getCoefCovariance", static_cast>> - (CovarianceReader::*)(double)>(&BasisClasses::CovarianceReader::getCoefCovariance), + (BasisClasses::CovarianceReader::*)(double)>(&BasisClasses::CovarianceReader::getCoefCovariance), R"( Get the covariance matrices for the basis coefficients @@ -2695,7 +2695,7 @@ void BasisFactoryClasses(py::module &m) )", py::arg("time")=0.0) .def("getCovarSamples", static_cast - (CovarianceReader::*)(unsigned)>(&BasisClasses::CovarianceReader::getCovarSamples), + (BasisClasses::CovarianceReader::*)(unsigned)>(&BasisClasses::CovarianceReader::getCovarSamples), R"( Get sample counts for the covariance computation @@ -2711,7 +2711,7 @@ void BasisFactoryClasses(py::module &m) )", py::arg("index")=0) .def("getCovarSamples", static_cast - (CovarianceReader::*)(double)>(&BasisClasses::CovarianceReader::getCovarSamples), + (BasisClasses::CovarianceReader::*)(double)>(&BasisClasses::CovarianceReader::getCovarSamples), R"( Get sample counts for the covariance computation From f5d839d7053318d32339b4be988e0836bf0a20a3 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 19 Oct 2025 17:29:13 -0400 Subject: [PATCH 19/82] Data was not copied to CovarianceReader internal array --- expui/BiorthBasis.cc | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 8719e7453..b1bad42de 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -5068,7 +5068,11 @@ namespace BasisClasses sampleMasses.push_back(Eigen::VectorXd()); stanza.getDataSet("sampleMasses").read(sampleMasses.back()); - size_t nT = sampleCounts.back().size(); + int nT = sampleCounts.back().size(); + + // Allocate vector for current time + covarData.push_back(std::vector>(nT)); + for (int T=0; T Date: Sun, 19 Oct 2025 17:33:57 -0400 Subject: [PATCH 20/82] Update expui/BiorthBasis.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index b1bad42de..9a1ac70f4 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -5006,7 +5006,7 @@ namespace BasisClasses CovarianceReader::CovarianceReader(const std::string& filename, int stride) { try { - // Create a new hdf5 file + // Open an existing hdf5 file for reading // HighFive::File file(filename, HighFive::File::ReadOnly); From 028149c22d6c20da08c07ba14d234ae9e9e03c7a Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 17:34:33 -0400 Subject: [PATCH 21/82] Update expui/BiorthBasis.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 9a1ac70f4..2ff18bffa 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -4917,8 +4917,8 @@ namespace BasisClasses // Check for version string std::string path = "CovarianceFileVersion"; - // Check for valid HDF file - if (file.exist(path)) { + // Check for valid HDF file by attribute + if (file.hasAttribute(path)) { extendCoefCovariance(fname, time); return; } From 3d6644ee596c64d1b45ae5651cc92a65206456b6 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 19 Oct 2025 17:40:23 -0400 Subject: [PATCH 22/82] Generalize type detection for covariance geometry --- expui/BiorthBasis.cc | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index b1bad42de..60ea8a25d 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1,4 +1,5 @@ #include +#include #include #include @@ -5025,16 +5026,22 @@ namespace BasisClasses int lmax, nmax, ltot; - if (basisID == "Spherical") { + // Current implemented spherical types + const std::set sphereType = {"Spherical", "SphereSL", "bessel"}; + + // Currently implemented cylindrical types + const std::set cylinderType = {"Cylindrical"}; + + if (sphereType.find(basisID) != sphereType.end()) { file.getAttribute("lmax").read(lmax); file.getAttribute("nmax").read(nmax); ltot = (lmax+1)*(lmax+2)/2; - } else if (basisID == "Cylindrical") { + } else if (cylinderType.find(basisID) != cylinderType.end()) { file.getAttribute("mmax").read(lmax); file.getAttribute("nmax").read(nmax); ltot = lmax + 1; } else { - throw std::runtime_error("CovarianceReader: unknown basis type, " + basisID); + throw std::runtime_error("CovarianceReader: unknown or unimplemented covariance for basis type, " + basisID); } // Group count variable From c6e4921375391deb87524b118f8d08554f7d8379 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 19 Oct 2025 17:57:19 -0400 Subject: [PATCH 23/82] More minor fixes based on a Copilot review --- expui/BiorthBasis.H | 2 ++ expui/BiorthBasis.cc | 11 +++++------ pyEXP/BasisWrappers.cc | 6 ++++-- 3 files changed, 11 insertions(+), 8 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 494fdba6b..a64dea90b 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -3,6 +3,8 @@ #include #include +#include +#include #include #include // For 3d rectangular grids diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 7e272be25..7722850e0 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1,5 +1,4 @@ #include -#include #include #include @@ -374,7 +373,7 @@ namespace BasisClasses { // Identifier // - BasisID = "sphereSL"; + BasisID = "SphereSL"; // Assign some defaults // @@ -438,7 +437,7 @@ namespace BasisClasses { // Identifier // - BasisID = "bessel"; + BasisID = "Bessel"; try { if (conf["rnum"]) @@ -1026,8 +1025,8 @@ namespace BasisClasses int ltot = (lmax+1)*(lmax+2)/2; ret[T].resize(ltot); for (int l=0; l(ret[T][l]) = meanV[T][l]; - std::get<1>(ret[T][l]) = covrV[T][l]; + std::get<0>(ret[T][l]) = std::move(meanV[T][l]); + std::get<1>(ret[T][l]) = std::move(covrV[T][l]); } } } @@ -5027,7 +5026,7 @@ namespace BasisClasses int lmax, nmax, ltot; // Current implemented spherical types - const std::set sphereType = {"Spherical", "SphereSL", "bessel"}; + const std::set sphereType = {"Spherical", "SphereSL", "Bessel"}; // Currently implemented cylindrical types const std::set cylinderType = {"Cylindrical"}; diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index e6bcf3138..8f0e4cbc7 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -2675,7 +2675,8 @@ void BasisFactoryClasses(py::module &m) Returns ------- list(list(tuple(numpy.ndarray, numpy.ndarray))) - list of covariance matrices for each subsample + list of partitioned coefficients and their covariance matrices for + each subsample )", py::arg("index")=0) .def("getCoefCovariance", static_cast>> @@ -2691,7 +2692,8 @@ void BasisFactoryClasses(py::module &m) Returns ------- list(list(tuple(numpy.ndarray, numpy.ndarray))) - list of covariance matrices for each subsample + list of partitioned coefficients and their covariance matrices for + each subsample )", py::arg("time")=0.0) .def("getCovarSamples", static_cast From 1d0b5d976cb2a294214d5ab648823e8b3d7fb41d Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 18:06:13 -0400 Subject: [PATCH 24/82] Update expui/BiorthBasis.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 7722850e0..739a49544 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1025,8 +1025,8 @@ namespace BasisClasses int ltot = (lmax+1)*(lmax+2)/2; ret[T].resize(ltot); for (int l=0; l(ret[T][l]) = std::move(meanV[T][l]); - std::get<1>(ret[T][l]) = std::move(covrV[T][l]); + std::get<0>(ret[T][l]) = meanV[T][l]; + std::get<1>(ret[T][l]) = covrV[T][l]; } } } From 3a178a7f2b5a2c1251847288600cebb68be85d28 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 18:06:51 -0400 Subject: [PATCH 25/82] Update expui/BiorthBasis.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 739a49544..0aee6fa37 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -4899,7 +4899,7 @@ namespace BasisClasses // Check that there is something to write // int totalCount = 0; - for (auto c : sampleCounts) totalCount += c; + for (int i = 0; i < sampleCounts.size(); ++i) totalCount += sampleCounts(i); if (totalCount==0) { std::cout << "BiorthBasis::writeCoefCovariance: no data" << std::endl; return; From 57656edbbeeadd51688a9b855e12c2cfa38a8cff Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 18:07:27 -0400 Subject: [PATCH 26/82] Update expui/BiorthBasis.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 0aee6fa37..26fd3f172 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -331,6 +331,7 @@ namespace BasisClasses // Set spherical coordindates // coordinates = Coord::Spherical; + BasisID = "Spherical"; } void Spherical::init_covariance() From b9d4122fa3f768bc8f0080346806e39002c30737 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 18:08:01 -0400 Subject: [PATCH 27/82] Update expui/BiorthBasis.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 26fd3f172..5575a86c5 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -5016,8 +5016,7 @@ namespace BasisClasses std::string version; file.getAttribute("CovarianceFileVersion").read(version); if (version != std::string("1.0")) { - throw std::runtime_error("CovarianceReader: unsupported file version, " - + version); + throw std::runtime_error(std::string("CovarianceReader: unsupported file version, ") + version); } // Read the basis identifier string From 28464543fcd3e1e23b637d215862743adc25f140 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 18:08:25 -0400 Subject: [PATCH 28/82] Update expui/BiorthBasis.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 5575a86c5..35fd1972b 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -5040,7 +5040,7 @@ namespace BasisClasses file.getAttribute("nmax").read(nmax); ltot = lmax + 1; } else { - throw std::runtime_error("CovarianceReader: unknown or unimplemented covariance for basis type, " + basisID); + throw std::runtime_error(std::string("CovarianceReader: unknown or unimplemented covariance for basis type, ") + basisID); } // Group count variable From 50dc8613e56c14796026ef1c9636a108b5bef8e0 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 18:08:48 -0400 Subject: [PATCH 29/82] 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 8f0e4cbc7..c9daa0ff4 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -2727,7 +2727,7 @@ void BasisFactoryClasses(py::module &m) tuple(numpy.ndarray, numpy.ndarray) sample counts and masses for the covariance computation )", - py::arg("time")=0) + py::arg("time")=0.0) .def("basisIDname", &BasisClasses::CovarianceReader::basisIDname, R"( Get the basis ID name From cca94d89895aaa7548d50d40c3b92752b3829301 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 18:09:19 -0400 Subject: [PATCH 30/82] Update exputil/EmpCylSL.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- exputil/EmpCylSL.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 383cfbae5..9f69ab90c 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -4817,7 +4817,7 @@ void EmpCylSL::pca_hall(bool compute, bool subsamp) std::vector> EmpCylSL::getCoefCovariance() { - std::vector>> ret; + std::vector> ret; if (PCAVAR) { ret.resize(sampT); From a336cdaf4ccffeeb6b5ff32c126323880bd9b34a Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 18:09:37 -0400 Subject: [PATCH 31/82] Update expui/BiorthBasis.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 35fd1972b..3ea23f1c3 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -4885,7 +4885,7 @@ namespace BasisClasses void BiorthBasis::writeCoefCovariance(const std::string& compname, const std::string& runtag, double time) { - // Check check that variance computation is on + // Check that variance computation is on // if (not pcavar) { std::cout << "BiorthBasis::writeCoefCovariance: covariance computation is disabled. " From 4ce5ee3ffef5068a6379cd202e79ce649dc594fb Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 19:26:20 -0400 Subject: [PATCH 32/82] Update pyEXP/BasisWrappers.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- pyEXP/BasisWrappers.cc | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index c9daa0ff4..f5c6f1597 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -1331,19 +1331,21 @@ void BasisFactoryClasses(py::module &m) py::arg("mass"), py::arg("pos")) .def("getCoefCovariance", &BasisClasses::BiorthBasis::getCoefCovariance, R"( - Return the coefficient vectors and covariance matrices for a partition - of the accumulated particles. The number of partitions is set by the - configuration parameters 'sampT' and defaults to 100. + Return the coefficient vectors and covariance matrices for each partition + of the accumulated particles. The number of partitions is set by the + configuration parameter 'sampT' (default: 100). - Parameters - ---------- - None + The coefficients are returned in the [specify: real basis (separate cos/sin for m>0; total angular dimension (lmax+1)^2) or complex triangular packing (l,m≥0; angular dimension (lmax+1)(lmax+2)/2)]. + Each coefficient vector has shape (angular_dimension,), and each covariance matrix has shape (angular_dimension, angular_dimension). + The ordering of harmonic indices within the coefficient vector and covariance matrix is [specify: e.g., for real basis, (l, m) ordered with m from -l to l for each l from 0 to lmax; for complex packing, (l, m) with m from 0 to l for each l from 0 to lmax]. + This allows consumers to unambiguously map (l, m) to positions in the arrays. Returns ------- arrays: list[list[tuple[np.ndarray, np.ndarray]]] - Each element is a tuple (coef_vector, coef_covariance_matrix), - where coef_vector and coef_covariance_matrix are numpy.ndarray. + Each element is a tuple (coef_vector, coef_covariance_matrix), + where coef_vector and coef_covariance_matrix are numpy.ndarray. + Each coef_covariance_matrix is of shape (angular_dimension, angular_dimension). See also -------- From 15b004fa2b2c1675aed205dce8db139dddc6f1b9 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 19:27:04 -0400 Subject: [PATCH 33/82] Update pyEXP/BasisWrappers.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- pyEXP/BasisWrappers.cc | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index f5c6f1597..42a959258 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -1347,6 +1347,18 @@ void BasisFactoryClasses(py::module &m) where coef_vector and coef_covariance_matrix are numpy.ndarray. Each coef_covariance_matrix is of shape (angular_dimension, angular_dimension). + Shape and Indexing + ------------------ + - Coefficient vectors are returned in the real basis, with separate cosine and sine + components for m>0. The total angular dimension is (lmax+1)^2, corresponding to + all (l, m) pairs with -l ≤ m ≤ l for each l = 0, ..., lmax. + - Each covariance matrix is of shape (nmax, nmax), where nmax is the number of basis + functions (typically (lmax+1)^2). + - The ordering of harmonic indices within the outer list is such that for each l from + 0 to lmax, m runs from -l to l, with cosine and sine components for m>0 stored + separately (e.g., [Y_{l,0}, Y_{l,1}^{cos}, Y_{l,1}^{sin}, ..., Y_{l,l}^{cos}, Y_{l,l}^{sin}]). + This allows consumers to unambiguously map (l, m) to positions in the coefficient vector + and covariance matrix. See also -------- getCovarSamples : get counts and mass in each partition From 2406860dfb72fafe3059e39228706f27f53b6d21 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 19 Oct 2025 19:49:23 -0400 Subject: [PATCH 34/82] More minor fixes based on a Copilot review --- expui/BiorthBasis.H | 2 ++ pyEXP/BasisWrappers.cc | 31 ++++++++++++++----------------- 2 files changed, 16 insertions(+), 17 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index a64dea90b..e00179b35 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -1000,6 +1000,8 @@ namespace BasisClasses std::tuple getCovarSamples() { + // Copy to base class members + std::tie(sampleCounts, sampleMasses) = sl->getCovarSamples(); return sl->getCovarSamples(); } }; diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 42a959258..dfe82ce6a 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -1331,34 +1331,31 @@ void BasisFactoryClasses(py::module &m) py::arg("mass"), py::arg("pos")) .def("getCoefCovariance", &BasisClasses::BiorthBasis::getCoefCovariance, R"( - Return the coefficient vectors and covariance matrices for each partition - of the accumulated particles. The number of partitions is set by the - configuration parameter 'sampT' (default: 100). + Return the coefficient vectors and covariance matrices for each partition of the + accumulated particles. The number of partitions is set by the configuration + parameter 'sampT' (default: 100). - The coefficients are returned in the [specify: real basis (separate cos/sin for m>0; total angular dimension (lmax+1)^2) or complex triangular packing (l,m≥0; angular dimension (lmax+1)(lmax+2)/2)]. - Each coefficient vector has shape (angular_dimension,), and each covariance matrix has shape (angular_dimension, angular_dimension). - The ordering of harmonic indices within the coefficient vector and covariance matrix is [specify: e.g., for real basis, (l, m) ordered with m from -l to l for each l from 0 to lmax; for complex packing, (l, m) with m from 0 to l for each l from 0 to lmax]. - This allows consumers to unambiguously map (l, m) to positions in the arrays. + The first dimension are the time samples. The second dimension is the angular + index. The coefficients are the moduli of the original complex coefficient vectors. Returns ------- arrays: list[list[tuple[np.ndarray, np.ndarray]]] Each element is a tuple (coef_vector, coef_covariance_matrix), where coef_vector and coef_covariance_matrix are numpy.ndarray. - Each coef_covariance_matrix is of shape (angular_dimension, angular_dimension). + Each coef_covariance_matrix is of shape (nmax, nmax) Shape and Indexing ------------------ - - Coefficient vectors are returned in the real basis, with separate cosine and sine - components for m>0. The total angular dimension is (lmax+1)^2, corresponding to - all (l, m) pairs with -l ≤ m ≤ l for each l = 0, ..., lmax. + - Coefficient vectors are moduli of the originall complex coefficient vectors + - The first list index is the number of time samples. + - The second list index is the angular elements. For spherical bases, all (l, m) pairs + are in triangular index order with l in [0,...,lmax] and m in [0,...,l] for a total of + (lmax+1)*(lmax+2)/2 entries. For cylindrical bases, there are (mmax+1) harmonic + entries entries for each value m in [0,...,mmax]. - Each covariance matrix is of shape (nmax, nmax), where nmax is the number of basis - functions (typically (lmax+1)^2). - - The ordering of harmonic indices within the outer list is such that for each l from - 0 to lmax, m runs from -l to l, with cosine and sine components for m>0 stored - separately (e.g., [Y_{l,0}, Y_{l,1}^{cos}, Y_{l,1}^{sin}, ..., Y_{l,l}^{cos}, Y_{l,l}^{sin}]). - This allows consumers to unambiguously map (l, m) to positions in the coefficient vector - and covariance matrix. + functions + See also -------- getCovarSamples : get counts and mass in each partition From b44d8bab65b9e1159024f425296af511a2fd2e98 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 19:59:29 -0400 Subject: [PATCH 35/82] 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 dfe82ce6a..ae7394d32 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -2706,7 +2706,7 @@ void BasisFactoryClasses(py::module &m) list of partitioned coefficients and their covariance matrices for each subsample )", - py::arg("time")=0.0) + py::arg("time")) .def("getCovarSamples", static_cast (BasisClasses::CovarianceReader::*)(unsigned)>(&BasisClasses::CovarianceReader::getCovarSamples), R"( From 044e9e4ccd225fb94045b2e834cb008a469b61e3 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 19:59:50 -0400 Subject: [PATCH 36/82] 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 ae7394d32..8efeae9b8 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -1347,7 +1347,7 @@ void BasisFactoryClasses(py::module &m) Shape and Indexing ------------------ - - Coefficient vectors are moduli of the originall complex coefficient vectors + - Coefficient vectors are moduli of the original complex coefficient vectors - The first list index is the number of time samples. - The second list index is the angular elements. For spherical bases, all (l, m) pairs are in triangular index order with l in [0,...,lmax] and m in [0,...,l] for a total of From c18ea29ff391072dd4997c31d0f6e24e3ab47590 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 20:00:12 -0400 Subject: [PATCH 37/82] 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 8efeae9b8..686ed4f3e 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -1352,7 +1352,7 @@ void BasisFactoryClasses(py::module &m) - The second list index is the angular elements. For spherical bases, all (l, m) pairs are in triangular index order with l in [0,...,lmax] and m in [0,...,l] for a total of (lmax+1)*(lmax+2)/2 entries. For cylindrical bases, there are (mmax+1) harmonic - entries entries for each value m in [0,...,mmax]. + entries for each value m in [0,...,mmax]. - Each covariance matrix is of shape (nmax, nmax), where nmax is the number of basis functions From ad8a5d467404c2bf27fb96dbbdfaa4c0fe426d74 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 20:00:55 -0400 Subject: [PATCH 38/82] Update include/EmpCylSL.H Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- include/EmpCylSL.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/EmpCylSL.H b/include/EmpCylSL.H index 39d04a108..22704a507 100644 --- a/include/EmpCylSL.H +++ b/include/EmpCylSL.H @@ -930,7 +930,7 @@ public: //! Check orthogonality for basis (pyEXP style) std::vector orthoCheck(); - //! Return a vector of tuples of basis functions and the + //! Return a vector of tuples of basis coefficients and the //! covariance matrix for subsamples of particles using CoefCovarType = std::tuple; std::vector> getCoefCovariance(); From 32b9a12b05850092232e0ef184fe50ac541ed81a Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 20:57:32 -0400 Subject: [PATCH 39/82] Update expui/BiorthBasis.H Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index e00179b35..2ba1f451f 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -1002,7 +1002,7 @@ namespace BasisClasses { // Copy to base class members std::tie(sampleCounts, sampleMasses) = sl->getCovarSamples(); - return sl->getCovarSamples(); + return {sampleCounts, sampleMasses}; } }; From 4accfdfdc4159992b76e53d7b93be9f81528b544 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 19 Oct 2025 22:23:52 -0400 Subject: [PATCH 40/82] Simplify pcavar computation by separating the computation from expcoef --- expui/BiorthBasis.cc | 30 +++++++++++++----------------- 1 file changed, 13 insertions(+), 17 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 3ea23f1c3..f26d7973d 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -613,44 +613,39 @@ namespace BasisClasses // M loop for (int m=0, moffset=0, moffE=0; m<=l; m++, L++) { + fac = factorial(l, m) * legs[tid](l, m); + if (m==0) { - fac = factorial(l, m) * legs[tid](l, m); for (int n=0; n Date: Sun, 19 Oct 2025 22:27:01 -0400 Subject: [PATCH 41/82] Use Eigen reduction for totalCounts accumulation --- expui/BiorthBasis.cc | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index f26d7973d..06102a1b1 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -4895,8 +4895,9 @@ namespace BasisClasses // Check that there is something to write // int totalCount = 0; - std::tie(sampleCounts, sampleMasses) = getCovarSamples() - for (int i = 0; i < sampleCounts.size(); ++i) totalCount += sampleCounts(i); + std::tie(sampleCounts, sampleMasses) = getCovarSamples(); + totalCount += sampleCounts.sum(); + if (totalCount==0) { std::cout << "BiorthBasis::writeCoefCovariance: no data" << std::endl; return; From b0a0cbb212ef687f449b92a78ec9ea37152d9337 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 19 Oct 2025 22:30:16 -0400 Subject: [PATCH 42/82] Remove default value for pybind11 overloads for covariance access --- pyEXP/BasisWrappers.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 686ed4f3e..c00feb849 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -2689,7 +2689,7 @@ void BasisFactoryClasses(py::module &m) list of partitioned coefficients and their covariance matrices for each subsample )", - py::arg("index")=0) + py::arg("index")) .def("getCoefCovariance", static_cast>> (BasisClasses::CovarianceReader::*)(double)>(&BasisClasses::CovarianceReader::getCoefCovariance), R"( @@ -2722,7 +2722,7 @@ void BasisFactoryClasses(py::module &m) tuple(numpy.ndarray, numpy.ndarray) sample counts and masses for the covariance computation )", - py::arg("index")=0) + py::arg("index")) .def("getCovarSamples", static_cast (BasisClasses::CovarianceReader::*)(double)>(&BasisClasses::CovarianceReader::getCovarSamples), R"( @@ -2738,7 +2738,7 @@ void BasisFactoryClasses(py::module &m) tuple(numpy.ndarray, numpy.ndarray) sample counts and masses for the covariance computation )", - py::arg("time")=0.0) + py::arg("time")) .def("basisIDname", &BasisClasses::CovarianceReader::basisIDname, R"( Get the basis ID name From ae7941f51785248923a4d031e89ca09e10eafc7e Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 19 Oct 2025 22:39:59 -0400 Subject: [PATCH 43/82] Update expui/BiorthBasis.cc Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 06102a1b1..2de5fe2e3 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -328,7 +328,7 @@ namespace BasisClasses // if (pcavar) init_covariance(); - // Set spherical coordindates + // Set spherical coordinates // coordinates = Coord::Spherical; BasisID = "Spherical"; From e4af5bfc6b89c628116808f95f7dfd0d35402ffb Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 20 Oct 2025 17:41:25 -0400 Subject: [PATCH 44/82] Added particle index to accumulate() members. Fixed Cylindrical calls with pcavar --- expui/BasisFactory.H | 8 ++-- expui/BiorthBasis.H | 43 ++++++++++++++++------ expui/BiorthBasis.cc | 83 ++++++++++++++++++++++++------------------ expui/FieldBasis.H | 9 +++-- expui/FieldBasis.cc | 23 ++++++------ exputil/EmpCylSL.cc | 10 ++--- include/EmpCylSL.H | 2 +- pyEXP/BasisWrappers.cc | 45 ++++++++++++----------- src/CylEXP.cc | 2 +- 9 files changed, 130 insertions(+), 95 deletions(-) diff --git a/expui/BasisFactory.H b/expui/BasisFactory.H index b0480da39..f5a7f80f6 100644 --- a/expui/BasisFactory.H +++ b/expui/BasisFactory.H @@ -184,13 +184,13 @@ namespace BasisClasses { return coefret; } //! Accumulate new coefficients from full phase space - virtual void accumulate(double mass, - double x, double y, double z, - double u, double v, double w) = 0; + virtual void accumulate(double x, double y, double z, + double u, double v, double w, + double mass, unsigned long int indx) = 0; //! Accumulate new coefficients from coordinates only virtual void accumulate(double x, double y, double z, - double mass) = 0; + double mass, unsigned long int indx) = 0; //! Get mass on grid double getMass(void) { return totalMass; } diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 2ba1f451f..37e53ccd9 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -110,6 +110,14 @@ namespace BasisClasses Eigen::VectorXd sampleMasses; //@} + //! 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; + } + //! Write H5 covariance; returns updated group count virtual unsigned writeCovarH5(HighFive::Group& group, unsigned count, double time); @@ -174,13 +182,14 @@ namespace BasisClasses //@} //! Accumulate new coefficients - virtual void accumulate(double x, double y, double z, double mass) = 0; + virtual void accumulate(double x, double y, double z, double mass, + unsigned long indx = 0) = 0; //! Accumulate new coefficients - virtual void accumulate(double mass, - double x, double y, double z, - double u, double v, double w) - { accumulate(x, y, z, mass); } + virtual void accumulate(double x, double y, double z, + double u, double v, double w, + double mass, unsigned long indx = 0) + { accumulate(x, y, z, mass, indx); } //! Get mass on grid double getMass(void) { return totalMass; } @@ -264,6 +273,9 @@ namespace BasisClasses const std::string& runtag, double time=0.0); + //! Make covariance after accumulation + virtual void makeCoefCovariance(void) {} + }; /** @@ -398,7 +410,8 @@ namespace BasisClasses void make_coefs(void); //! Accumulate new coefficients - virtual void accumulate(double x, double y, double z, double mass); + virtual void accumulate(double x, double y, double z, double mass, + unsigned long indx = 0); //! Return current maximum harmonic order in expansion int getLmax() { return lmax; } @@ -643,7 +656,8 @@ namespace BasisClasses void make_coefs(void); //! Accumulate new coefficients - virtual void accumulate(double x, double y, double z, double mass); + virtual void accumulate(double x, double y, double z, double mass, + unsigned long indx = 0); //! Return current maximum harmonic order in expansion int getMmax() { return mmax; } @@ -800,7 +814,8 @@ namespace BasisClasses void make_coefs(void); //! Accumulate new coefficients - virtual void accumulate(double x, double y, double z, double mass); + virtual void accumulate(double x, double y, double z, double mass, + unsigned long indx = 0); //! Return current maximum harmonic order in expansion int getMmax() { return mmax; } @@ -958,7 +973,8 @@ namespace BasisClasses void make_coefs(void); //! Accumulate new coefficients - virtual void accumulate(double x, double y, double z, double mass); + virtual void accumulate(double x, double y, double z, double mass, + unsigned long indx = 0); //! Return current maximum harmonic order in expansion int getMmax() { return mmax; } @@ -1004,6 +1020,7 @@ namespace BasisClasses std::tie(sampleCounts, sampleMasses) = sl->getCovarSamples(); return {sampleCounts, sampleMasses}; } + }; /** @@ -1120,7 +1137,8 @@ namespace BasisClasses void make_coefs(void); //! Accumulate new coefficients - virtual void accumulate(double x, double y, double z, double mass); + virtual void accumulate(double x, double y, double z, double mass, + unsigned long indx = 0); //! Return current maximum harmonic order in expansion Eigen::Vector3i getNmax() { return {nmaxx, nmaxy, nmaxz}; } @@ -1258,7 +1276,8 @@ namespace BasisClasses void make_coefs(void); //! Accumulate new coefficients - virtual void accumulate(double x, double y, double z, double mass); + virtual void accumulate(double x, double y, double z, double mass, + unsigned long indx = 0); //! Return current maximum harmonic order in expansion Eigen::Vector3i getNmax() { return {nmaxx, nmaxy, nmaxz}; } @@ -1302,7 +1321,7 @@ namespace BasisClasses std::vector sampleMasses; std::vector>> covarData; - const double fixedPointPrecision = 1.0e04; + const double fixedPointPrecision = 1.0e08; std::map timeMap; unsigned timeIndex(double time) { diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 2de5fe2e3..12ef002fc 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -265,7 +265,7 @@ namespace BasisClasses if (conf["EVEN_M"]) EVEN_M = conf["EVEN_M"].as(); if (conf["M0_ONLY"]) M0_only = conf["M0_ONLY"].as(); if (conf["pcavar"]) pcavar = conf["pcavar"].as(); - if (conf["sampT"]) sampT = conf["sampT"].as(); + if (conf["subsamp"]) sampT = conf["subsamp"].as(); } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing parameter stanza for <" @@ -568,7 +568,8 @@ namespace BasisClasses coefctr = {0.0, 0.0, 0.0}; } - void Spherical::accumulate(double x, double y, double z, double mass) + void Spherical::accumulate(double x, double y, double z, double mass, + unsigned long indx) { double fac, fac1, fac2, fac4; double norm = -4.0*M_PI; @@ -1444,6 +1445,8 @@ namespace BasisClasses if (unmatched.size()) throw YamlConfigError("Basis::Basis::Cylindrical", "parameter", unmatched, __FILE__, __LINE__); + int sampT = 0; + // Assign values from YAML // try { @@ -1495,6 +1498,7 @@ namespace BasisClasses if (conf["vflag" ]) vflag = conf["vflag" ].as(); if (conf["pyname" ]) pyname = conf["pyname" ].as(); if (conf["pcavar"] ) pcavar = conf["pcavar" ].as(); + if (conf["subsamp"] ) sampT = conf["subsamp" ].as(); // Deprecation warning if (conf["density" ]) { @@ -1544,6 +1548,7 @@ namespace BasisClasses EmpCylSL::CMAPZ = cmapZ; EmpCylSL::logarithmic = logarithmic; EmpCylSL::VFLAG = vflag; + EmpCylSL::PCAVAR = pcavar; // Check for non-null cache file name. This must be specified // to prevent recomputation and unexpected behavior. @@ -1566,7 +1571,11 @@ namespace BasisClasses // if (mlim>=0) sl->set_mlim(mlim); if (EVEN_M) sl->setEven(EVEN_M); - + if (sampT>0) { + sl->setSampT(sampT); + sl->init_pca(); + } + // Cache override for old Eigen cache // if (oldcache) sl->AllowOldCache(); @@ -1790,11 +1799,12 @@ namespace BasisClasses } } - void Cylindrical::accumulate(double x, double y, double z, double mass) + void Cylindrical::accumulate(double x, double y, double z, double mass, + unsigned long indx) { double R = sqrt(x*x + y*y); double phi = atan2(y, x); - sl->accumulate(R, z, phi, mass, 0, 0); + sl->accumulate(R, z, phi, mass, indx, 0, 0, pcavar); } void Cylindrical::reset_coefs(void) @@ -1857,7 +1867,7 @@ namespace BasisClasses void Cylindrical::make_coefs(void) { - sl->make_coefficients(); + sl->make_coefficients(pcavar); } @@ -2170,7 +2180,8 @@ namespace BasisClasses coefctr = {0.0, 0.0, 0.0}; } - void FlatDisk::accumulate(double x, double y, double z, double mass) + void FlatDisk::accumulate(double x, double y, double z, double mass, + unsigned long indx) { // Normalization factors // @@ -2974,7 +2985,8 @@ namespace BasisClasses coefctr = {0.0, 0.0, 0.0}; } - void CBDisk::accumulate(double x, double y, double z, double mass) + void CBDisk::accumulate(double x, double y, double z, double mass, + unsigned long int idx) { // Normalization factors // @@ -3432,7 +3444,8 @@ namespace BasisClasses coefctr = {0.0, 0.0, 0.0}; } - void Slab::accumulate(double x, double y, double z, double mass) + void Slab::accumulate(double x, double y, double z, double mass, + unsigned long int idx) { // Truncate to slab with sides in [0,1] if (x<0.0) @@ -3949,7 +3962,8 @@ namespace BasisClasses coefctr = {0.0, 0.0, 0.0}; } - void Cube::accumulate(double x, double y, double z, double mass) + void Cube::accumulate(double x, double y, double z, double mass, + unsigned long int indx) { // Truncate to cube with sides in [0,1] if (x<0.0) @@ -4179,7 +4193,7 @@ namespace BasisClasses } if (use) { - accumulate(pp(0), pp(1), pp(2), p->mass); + accumulate(pp(0), pp(1), pp(2), p->mass, p->indx); } } make_coefs(); @@ -4298,7 +4312,7 @@ namespace BasisClasses coefindx++; if (use) { - accumulate(pp(0), pp(1), pp(2), m(n)); + accumulate(pp(0), pp(1), pp(2), m(n), n); } } } @@ -4336,7 +4350,7 @@ namespace BasisClasses coefindx++; if (use) { - accumulate(pp(0), pp(1), pp(2), m(n)); + accumulate(pp(0), pp(1), pp(2), m(n), n); } } } @@ -4826,6 +4840,7 @@ namespace BasisClasses // Add a time attribute // + time = roundTime(time); stanza.createAttribute("Time", HighFive::DataSpace::From(time)).write(time); // Add the sample statistics @@ -4903,14 +4918,22 @@ namespace BasisClasses return; } + // Round time + // + time = roundTime(time); + // The H5 filename + // std::string fname = "coefcovar." + compname + "." + runtag + ".h5"; // Check if file exists? // try { - // Open the HDF5 file in read-only mode - HighFive::File file(fname, HighFive::File::ReadOnly); + // Open the HDF5 file in read-write mode, creating if it doesn't + // exist + HighFive::File file(fname, + HighFive::File::ReadWrite | + HighFive::File::Create); // Check for version string std::string path = "CovarianceFileVersion"; @@ -4921,23 +4944,6 @@ namespace BasisClasses return; } - } catch (const HighFive::Exception& err) { - // Handle HighFive specific errors (e.g., file not found) - std::cerr << "HighFive Error: " << err.what() << std::endl; - } catch (const std::exception& err) { - // Handle other general exceptions - std::cerr << "Error: " << err.what() << std::endl; - } - - // Write coefficients - // - try { - // Create a new hdf5 file - // - HighFive::File file(fname, - HighFive::File::ReadWrite | - HighFive::File::Create); - // Write the Version string // file.createAttribute("CovarianceFileVersion", HighFive::DataSpace::From(CovarianceFileVersion)).write(CovarianceFileVersion); @@ -4967,8 +4973,14 @@ namespace BasisClasses // dataset.write(count); - } catch (HighFive::Exception& err) { - std::cerr << err.what() << std::endl; + } catch (const HighFive::Exception& err) { + // Handle HighFive specific errors (e.g., file not found) + throw std::runtime_error + (std::string("BiorthBasis::writeCoefCovariance HighFive Error: ") + err.what()); + } catch (const std::exception& err) { + // Handle other general exceptions + throw std::runtime_error + (std::string("BiorthBasis::writeCoefCovariance Error: ") + err.what()); } } @@ -4996,7 +5008,8 @@ namespace BasisClasses dataset.write(count); } catch (HighFive::Exception& err) { - std::cerr << err.what() << std::endl; + throw std::runtime_error + (std::string("BiorthBasis::extendCoefCovariance: HighFive error: ") + err.what()); } } diff --git a/expui/FieldBasis.H b/expui/FieldBasis.H index 37fc20044..1e2a2ba50 100644 --- a/expui/FieldBasis.H +++ b/expui/FieldBasis.H @@ -156,12 +156,13 @@ namespace BasisClasses (Eigen::VectorXd& m, RowMatrixXd& p, bool roundrobin=true, bool posvelrows=false); //! Accumulate new coefficients - virtual void accumulate(double mass, - double x, double y, double z, - double u, double v, double w); + virtual void accumulate(double x, double y, double z, + double u, double v, double w, + double mass, unsigned long int indx = 0); [[deprecated("not relevant for this class")]] - virtual void accumulate(double x, double y, double z, double mass) + virtual void accumulate(double x, double y, double z, double mass, + unsigned long indx = 0) { } diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc index 66a928f34..168d44470 100644 --- a/expui/FieldBasis.cc +++ b/expui/FieldBasis.cc @@ -316,9 +316,9 @@ namespace BasisClasses } } - void FieldBasis::accumulate(double mass, - double x, double y, double z, - double u, double v, double w) + void FieldBasis::accumulate(double x, double y, double z, + double u, double v, double w, + double mass, unsigned long int indx) { constexpr double fac2 = 0.5*M_2_SQRTPI/M_SQRT2; // 1/sqrt(2*pi)=0.3989422804 @@ -580,13 +580,13 @@ namespace BasisClasses use = true; } - if (use) accumulate(p->mass, - p->pos[0]-ctr[0], + if (use) accumulate(p->pos[0]-ctr[0], p->pos[1]-ctr[1], p->pos[2]-ctr[2], p->vel[0], p->vel[1], - p->vel[2]); + p->vel[2], + p->mass, p->indx); } make_coefs(); @@ -672,13 +672,13 @@ namespace BasisClasses } coefindx++; - if (use) accumulate(m(n), - p(0, n)-coefctr[0], + if (use) accumulate(p(0, n)-coefctr[0], p(1, n)-coefctr[1], p(2, n)-coefctr[2], p(3, n), p(4, n), - p(5, n)); + p(5, n), + m(n), n); } } @@ -707,13 +707,12 @@ namespace BasisClasses } coefindx++; - if (use) accumulate(m(n), - p(n, 0)-coefctr[0], + if (use) accumulate(p(n, 0)-coefctr[0], p(n, 1)-coefctr[1], p(n, 2)-coefctr[2], p(n, 3), p(n, 4), - p(n, 5)); + p(n, 5), m(n), n); } } } diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 9f69ab90c..5c6aec27a 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -1923,7 +1923,7 @@ void EmpCylSL::setup_accumulation(int mlevel) } } - if (PCAVAR and sampT>0) { + if (PCAVAR and sampT>1) { for (int nth=0; nth0) accum_sin[m].setZero(); } - if ( (PCAVAR or PCAEOF) and mlevel==0 and sampT>0) { + if ( (PCAVAR or PCAEOF) and mlevel==0 and sampT>1) { for (int nth=0; nth coefs_made; + std::vector coefs_made; bool eof_made; SphModTblPtr make_sl(); diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index c00feb849..fb87c00f0 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -279,8 +279,8 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE(FCReturn, Basis, getFieldsCoefs, x, y, z, coefs); } - void accumulate(double x, double y, double z, double mass) override { - PYBIND11_OVERRIDE_PURE(void, Basis, accumulate, mass, x, y, z, mass); + void accumulate(double x, double y, double z, double mass, unsigned long int indx) override { + PYBIND11_OVERRIDE_PURE(void, Basis, accumulate, mass, x, y, z, mass, indx); } void reset_coefs(void) override { @@ -360,9 +360,10 @@ void BasisFactoryClasses(py::module &m) } void accumulate(double m, double x, double y, double z, - double u, double v, double w) override + double u, double v, double w, + unsigned long int indx) override { - PYBIND11_OVERRIDE(void, FieldBasis, accumulate, m, x, y, z, u, v, w); + PYBIND11_OVERRIDE(void, FieldBasis, accumulate, m, x, y, z, u, v, w, indx); } void reset_coefs(void) override { @@ -419,8 +420,8 @@ void BasisFactoryClasses(py::module &m) // Inherit the constructors using BasisClasses::BiorthBasis::BiorthBasis; - void accumulate(double x, double y, double z, double mass) override { - PYBIND11_OVERRIDE_PURE(void, BiorthBasis, accumulate, x, y, z, mass); + void accumulate(double x, double y, double z, double mass, unsigned long int indx) override { + PYBIND11_OVERRIDE_PURE(void, BiorthBasis, accumulate, x, y, z, mass, indx); } void reset_coefs(void) override { @@ -497,8 +498,8 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE(std::vector, Spherical, getFields, x, y, z); } - void accumulate(double x, double y, double z, double mass) override { - PYBIND11_OVERRIDE(void, Spherical, accumulate, x, y, z, mass); + void accumulate(double x, double y, double z, double mass, unsigned long int indx) override { + PYBIND11_OVERRIDE(void, Spherical, accumulate, x, y, z, mass, indx); } void reset_coefs(void) override { @@ -552,8 +553,8 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE(std::vector, Cylindrical, getFields, x, y, z); } - void accumulate(double x, double y, double z, double mass) override { - PYBIND11_OVERRIDE(void, Cylindrical, accumulate, x, y, z, mass); + void accumulate(double x, double y, double z, double mass, unsigned long int indx) override { + PYBIND11_OVERRIDE(void, Cylindrical, accumulate, x, y, z, mass, indx); } void reset_coefs(void) override { @@ -621,9 +622,9 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE(std::vector, FlatDisk, getFields, x, y, z); } - void accumulate(double x, double y, double z, double mass) override + void accumulate(double x, double y, double z, double mass, unsigned long int indx) override { - PYBIND11_OVERRIDE(void, FlatDisk, accumulate, x, y, z, mass); + PYBIND11_OVERRIDE(void, FlatDisk, accumulate, x, y, z, mass, indx); } void reset_coefs(void) override @@ -692,9 +693,9 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE(std::vector, CBDisk, getFields, x, y, z); } - void accumulate(double x, double y, double z, double mass) override + void accumulate(double x, double y, double z, double mass, unsigned long int indx) override { - PYBIND11_OVERRIDE(void, CBDisk, accumulate, x, y, z, mass); + PYBIND11_OVERRIDE(void, CBDisk, accumulate, x, y, z, mass, indx); } void reset_coefs(void) override @@ -766,9 +767,9 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE(std::vector, Slab, getFields, x, y, z); } - void accumulate(double x, double y, double z, double mass) override + void accumulate(double x, double y, double z, double mass, unsigned long int indx) override { - PYBIND11_OVERRIDE(void, Slab, accumulate, x, y, z, mass); + PYBIND11_OVERRIDE(void, Slab, accumulate, x, y, z, mass, indx); } void reset_coefs(void) override @@ -840,9 +841,9 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE(std::vector, Cube, getFields, x, y, z); } - void accumulate(double x, double y, double z, double mass) override + void accumulate(double x, double y, double z, double mass, unsigned long int indx) override { - PYBIND11_OVERRIDE(void, Cube, accumulate, x, y, z, mass); + PYBIND11_OVERRIDE(void, Cube, accumulate, x, y, z, mass, indx); } void reset_coefs(void) override @@ -1557,9 +1558,9 @@ void BasisFactoryClasses(py::module &m) ------- None )") - .def("accumulate", [](BasisClasses::BiorthBasis& A, double x, double y, double z, double mass) + .def("accumulate", [](BasisClasses::BiorthBasis& A, double x, double y, double z, double mass, unsigned long int indx) { - return A.accumulate(x, y, z, mass); + return A.accumulate(x, y, z, mass, indx); }, R"( Add the contribution of a single particle to the coefficients @@ -1574,12 +1575,14 @@ void BasisFactoryClasses(py::module &m) z-axis position mass : float particle mass + indx : unsigned long int + particle index for selection functor Returns ------- None )", - py::arg("x"), py::arg("y"), py::arg("z"), py::arg("mass")) + py::arg("x"), py::arg("y"), py::arg("z"), py::arg("mass"), py::arg("indx")=0) .def("getMass", &BasisClasses::BiorthBasis::getMass, R"( Return the total mass of particles contributing the current coefficient set diff --git a/src/CylEXP.cc b/src/CylEXP.cc index d76abd1e0..5256f2cf5 100644 --- a/src/CylEXP.cc +++ b/src/CylEXP.cc @@ -278,7 +278,7 @@ void CylEXP::compute_multistep_coefficients(unsigned mlevel) } } - coefs_made = vector(multistep+1, true); + coefs_made = std::vector(multistep+1, true); } From 18c4c262ae6a0ff88c0559626d05b9db5974274b Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Mon, 20 Oct 2025 18:20:19 -0400 Subject: [PATCH 45/82] 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 fb87c00f0..cdd9fee35 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -280,7 +280,7 @@ void BasisFactoryClasses(py::module &m) } void accumulate(double x, double y, double z, double mass, unsigned long int indx) override { - PYBIND11_OVERRIDE_PURE(void, Basis, accumulate, mass, x, y, z, mass, indx); + PYBIND11_OVERRIDE_PURE(void, Basis, accumulate, x, y, z, mass, indx); } void reset_coefs(void) override { From d9c68c9860b60b9422f6d193ab943da21c8dadab Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 20 Oct 2025 19:48:32 -0400 Subject: [PATCH 46/82] Minor formatting changes --- exputil/EmpCylSL.cc | 2 +- pyEXP/BasisWrappers.cc | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 5c6aec27a..c2f6a992f 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -4841,7 +4841,7 @@ EmpCylSL::getCovarSamples() if (PCAVAR) { std::get<0>(ret).resize(sampT); std::get<1>(ret).resize(sampT); - for (int T=0; T(ret)(T) = numbT[T]; std::get<1>(ret)(T) = massT[T]; } diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index fb87c00f0..06189d282 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -1575,7 +1575,7 @@ void BasisFactoryClasses(py::module &m) z-axis position mass : float particle mass - indx : unsigned long int + indx : unsigned long int particle index for selection functor Returns From cfa363b04a64f0fffa4b409e20c618d5eed57315 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 24 Oct 2025 17:47:14 -0400 Subject: [PATCH 47/82] Some minor logic fixes for EmpCylSL --- expui/BiorthBasis.H | 8 +- expui/BiorthBasis.cc | 108 +++++++++++++++----- exputil/EmpCylSL.cc | 218 ++++++++++++++++++++++++++++++++--------- include/EmpCylSL.H | 18 +++- pyEXP/BasisWrappers.cc | 40 +++++--- 5 files changed, 306 insertions(+), 86 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index c71551bfe..ce46ca72f 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -101,6 +101,9 @@ namespace BasisClasses //! Coefficient variance computation enabled bool pcavar = false; + //! Data type for covariance + bool floatType = false; + //@{ //! Sample counts and masses for covariance computation Eigen::VectorXi sampleCounts; @@ -270,6 +273,9 @@ namespace BasisClasses const std::string& runtag, double time=0.0); + //! Choose between float and double storage for covariance + void setCovarFloatType(bool flt) { floatType = flt; } + //! Make covariance after accumulation virtual void makeCoefCovariance(void) {} @@ -1328,7 +1334,7 @@ namespace BasisClasses { public: - using CoefCovarType = std::tuple; + using CoefCovarType = std::tuple; protected: diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 3f12e9746..356240cec 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -637,9 +637,11 @@ namespace BasisClasses } if (pcavar) { - Eigen::VectorXd g = potd[tid].row(l).transpose() * fac * norm; + Eigen::VectorXcd g = std::exp(std::complex(0.0, m*phi)) * + potd[tid].row(l).transpose() * fac * norm; + meanV[T][L].noalias() += g * mass; - covrV[T][L].noalias() += g * g.transpose() * mass; + covrV[T][L].noalias() += g * g.adjoint() * mass; } } @@ -1547,8 +1549,7 @@ namespace BasisClasses EmpCylSL::CMAPZ = cmapZ; EmpCylSL::logarithmic = logarithmic; EmpCylSL::VFLAG = vflag; - EmpCylSL::PCAVAR = pcavar; - + // Check for non-null cache file name. This must be specified // to prevent recomputation and unexpected behavior. // @@ -1570,9 +1571,9 @@ namespace BasisClasses // if (mlim>=0) sl->set_mlim(mlim); if (EVEN_M) sl->setEven(EVEN_M); - if (sampT>0) { + if (pcavar) { sl->setSampT(sampT); - sl->init_pca(); + sl->set_covar(true); } // Cache override for old Eigen cache @@ -4869,28 +4870,56 @@ namespace BasisClasses size_t lmax = covar[T].size(); size_t nmax = std::get<0>(covar[T][0]).rows(); - Eigen::VectorXd data(lmax*nmax); - for (size_t l=0, c=0; l(covar[T][l])(n); + if (floatType) { + Eigen::VectorXcf data(lmax*nmax); + + for (size_t l=0, c=0; l> + (std::get<0>(covar[T][l])(n)); + } } - } - HighFive::DataSet coefdata = sample.createDataSet("coefficients", data); + HighFive::DataSet coefdata = sample.createDataSet("coefficients", data); - // Pack the covariance data in an upper triangular format - // - size_t diagonalSize = nmax*(nmax + 1)/2; - data.resize(lmax*diagonalSize); - for (size_t l=0, c=0; l(covar[T][l])(n1, n2); + // Pack the covariance data in an upper triangular format + // + size_t diagonalSize = nmax*(nmax + 1)/2; + data.resize(lmax*diagonalSize); + for (size_t l=0, c=0; l> + (std::get<1>(covar[T][l])(n1, n2)); + } + } + } + + HighFive::DataSet covdata = sample.createDataSet("covariance", data); + } else { + Eigen::VectorXcd data(lmax*nmax); + for (size_t l=0, c=0; l(covar[T][l])(n); } } - } - HighFive::DataSet covdata = sample.createDataSet("covariance", data); + HighFive::DataSet coefdata = sample.createDataSet("coefficients", data); + + // Pack the covariance data in an upper triangular format + // + size_t diagonalSize = nmax*(nmax + 1)/2; + data.resize(lmax*diagonalSize); + for (size_t l=0, c=0; l(covar[T][l])(n1, n2); + } + } + } + + HighFive::DataSet covdata = sample.createDataSet("covariance", data); + } } // END: sample loop @@ -4952,10 +4981,15 @@ namespace BasisClasses // file.createAttribute("CovarianceFileVersion", HighFive::DataSpace::From(CovarianceFileVersion)).write(CovarianceFileVersion); - // We write the basis identifier string + // Write the basis identifier string // file.createAttribute("BasisID", HighFive::DataSpace::From(BasisID)).write(BasisID); + // Write the data type size + // + int sz = 8; if (floatType) sz = 4; + file.createAttribute("FloatSize", HighFive::DataSpace::From(sz)).write(sz); + // Write the specific parameters // writeCovarH5Params(file); @@ -5037,6 +5071,16 @@ namespace BasisClasses // file.getAttribute("BasisID").read(basisID); + // Get the float size + int sz = 8; + file.getAttribute("FloatSize").read(sz); + if (sz != 4 and sz != 8) { + std::ostringstream sout; + sout << "CovarianceReader: unsupported float size, " << sz; + throw std::runtime_error(sout.str()); + } + std::cout << "Float size is " << sz << std::endl; + int lmax, nmax, ltot; // Current implemented spherical types @@ -5102,7 +5146,7 @@ namespace BasisClasses HighFive::Group sample = stanza.getGroup(sT.str()); // Storage - Eigen::VectorXd data; + Eigen::VectorXcd data; // Repack the data std::vector elem(ltot); @@ -5112,7 +5156,13 @@ namespace BasisClasses } // Get the flattened coefficient array - data = sample.getDataSet("coefficients").read(); + if (sz==4) { + Eigen::VectorXcf dataF; + dataF = sample.getDataSet("coefficients").read(); + data = dataF.cast>(); + } else { + data = sample.getDataSet("coefficients").read(); + } // Pack the coefficient data for (size_t l=0, c=0; l(); + if (sz==4) { + Eigen::VectorXcf dataF; + dataF = sample.getDataSet("covariance").read(); + data = dataF.cast>(); + } else { + data = sample.getDataSet("covariance").read(); + } // Pack the coefficient data for (size_t l=0, c=0; l::max(); maxSNR = 0.0; for (unsigned T=0; T(mcos, msin) * vc[id].row(mm).transpose() * norm; + VC[id][whch][mm] += mass * vec; + MV[id][whch][mm] += mass * (vec * vec.adjoint()); + } + if (compute and PCAVAR) { double hc1 = vc[id](mm, nn)*mcos, hs1 = 0.0; if (mm) hs1 = vs[id](mm, nn)*msin; @@ -4137,38 +4199,69 @@ void EmpCylSL::make_coefficients(unsigned M0, bool compute) tvar[0][mm](nn, oo) += tvar[nth][mm](nn, oo); } - for (unsigned T=0; T> ret; - if (PCAVAR) { + if (covar) { ret.resize(sampT); for (unsigned T=0; T(ret[T][M]) = covV[0][T][M]; - std::get<1>(ret[T][M]) = covM[0][T][M]; + std::get<0>(ret[T][M]) = VC[0][T][M]; + std::get<1>(ret[T][M]) = MV[0][T][M]; } } } @@ -4838,7 +4964,7 @@ EmpCylSL::getCovarSamples() { std::tuple ret; - if (PCAVAR) { + if (covar) { std::get<0>(ret).resize(sampT); std::get<1>(ret).resize(sampT); for (unsigned T=0; T > > >; VarMat SC, SS, SCe, SCo, SSe, SSo; + + //! Complex covariance implemenation + using CVarMat = std::vector< std::vector< std::vector< Eigen::MatrixXcd> > >; + using CVarVec = std::vector< std::vector< std::vector< Eigen::VectorXcd> > >; + CVarMat MV; + CVarVec VC; + bool covar = false; //@} std::vector var, varE, varO; @@ -630,6 +637,15 @@ public: //! Initialize PCA work space void init_pca(); + //! Initialize covariance work space + void init_covar(); + + //! Turn on/off covariance computation + void set_covar(bool tf) { + covar = tf; + if (tf) init_covar(); + } + //! Necessary member function currently unused (change design?) void determine_coefficients() {}; @@ -932,7 +948,7 @@ public: //! Return a vector of tuples of basis coefficients and the //! covariance matrix for subsamples of particles - using CoefCovarType = std::tuple; + using CoefCovarType = std::tuple; std::vector> getCoefCovariance(); //! Covariance statistics diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 133601e12..fc5f76ab5 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -436,7 +436,7 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE_PURE(void, BiorthBasis, set_coefs, coefs); } - using CCelement = std::tuple; + using CCelement = std::tuple; using CCreturn = std::vector>; CCreturn getCoefCovariance(void) override { @@ -1332,23 +1332,24 @@ void BasisFactoryClasses(py::module &m) py::arg("mass"), py::arg("pos")) .def("getCoefCovariance", &BasisClasses::BiorthBasis::getCoefCovariance, R"( - Return the coefficient vectors and covariance matrices for each partition of the + Return the coefficient vectors and covariance matrices for each partition of the accumulated particles. The number of partitions is set by the configuration parameter 'sampT' (default: 100). The first dimension are the time samples. The second dimension is the angular - index. The coefficients are the moduli of the original complex coefficient vectors. + index. Each element is a tuple of the coefficient vector and covariance. + The values are complex128 and represet the full amplitude and phase information. Returns ------- arrays: list[list[tuple[np.ndarray, np.ndarray]]] Each element is a tuple (coef_vector, coef_covariance_matrix), where coef_vector and coef_covariance_matrix are numpy.ndarray. - Each coef_covariance_matrix is of shape (nmax, nmax) + Each coef_covariance_matrix is of shape (nmax, nmax). All values are + complex128. Shape and Indexing ------------------ - - Coefficient vectors are moduli of the original complex coefficient vectors - The first list index is the number of time samples. - The second list index is the angular elements. For spherical bases, all (l, m) pairs are in triangular index order with l in [0,...,lmax] and m in [0,...,l] for a total of @@ -1365,7 +1366,7 @@ void BasisFactoryClasses(py::module &m) )") .def("getCovarSamples", &BasisClasses::BiorthBasis::getCovarSamples, R"( - Return a vector of counts and mass used to compute the partitioned + Return a vector of counts and mass used to compute the partitioned vectors and covariance. The length of both returned vectors equals 'sampT' Parameters @@ -1387,7 +1388,7 @@ void BasisFactoryClasses(py::module &m) )") .def("writeCoefCovariance", &BasisClasses::BiorthBasis::writeCoefCovariance, R"( - Write the partitioned coefficient vectors and covariance matrices + 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. @@ -1415,7 +1416,22 @@ void BasisFactoryClasses(py::module &m) getCoefCovariance : get the coefficient vectors and covariance matrices for the partitioned phase space. getCovarSamples : get counts and mass in each partition - )", py::arg("compname"), py::arg("runtag"), py::arg("time")=0.0) + )", py::arg("compname"), py::arg("runtag"), py::arg("time")=0.0) + .def("setCovarFloatType", &BasisClasses::BiorthBasis::setCovarFloatType, + R"( + Set the floating point type used for covariance storage in HDF5. Set to + true for float4. Type float8 is the default (false). This does not affect + the return type in getCoefCovariance(). + + Parameters + ---------- + floatType : bool + Use 'float32' rather than 'float64' if true + + Returns + ------- + None + )") .def("getFields", &BasisClasses::BiorthBasis::getFields, R"( Return the field evaluations for a given cartesian position. The @@ -2721,7 +2737,7 @@ void BasisFactoryClasses(py::module &m) times : list(float) a list of evaluation times )") - .def("getCoefCovariance", static_cast>> + .def("getCoefCovariance", static_cast>> (BasisClasses::CovarianceReader::*)(unsigned)>(&BasisClasses::CovarianceReader::getCoefCovariance), R"( Get the covariance matrices for the basis coefficients @@ -2737,7 +2753,7 @@ void BasisFactoryClasses(py::module &m) each subsample )", py::arg("index")) - .def("getCoefCovariance", static_cast>> + .def("getCoefCovariance", static_cast>> (BasisClasses::CovarianceReader::*)(double)>(&BasisClasses::CovarianceReader::getCoefCovariance), R"( Get the covariance matrices for the basis coefficients @@ -2750,8 +2766,8 @@ void BasisFactoryClasses(py::module &m) Returns ------- list(list(tuple(numpy.ndarray, numpy.ndarray))) - list of partitioned coefficients and their covariance matrices for - each subsample + List of partitioned coefficients and their covariance matrices for + each subsample. The returns are complex-valued. )", py::arg("time")) .def("getCovarSamples", static_cast From feb7a16044f7437fc0400a2117acdda94d827ca2 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 24 Oct 2025 22:52:47 -0400 Subject: [PATCH 48/82] Remove an overlooked debug line --- expui/BiorthBasis.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 356240cec..4a7cdddd1 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -5079,7 +5079,6 @@ namespace BasisClasses sout << "CovarianceReader: unsupported float size, " << sz; throw std::runtime_error(sout.str()); } - std::cout << "Float size is " << sz << std::endl; int lmax, nmax, ltot; From 1da2921a9d2bf0920b20e565cb89bbe5848028c8 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 25 Oct 2025 22:40:20 -0400 Subject: [PATCH 49/82] Implemente HDF5 compression --- expui/BiorthBasis.H | 22 ++++++- expui/BiorthBasis.cc | 136 ++++++++++++++++++++++++++++++++++++----- pyEXP/BasisWrappers.cc | 17 ++++++ 3 files changed, 159 insertions(+), 16 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index ce46ca72f..0bab9c5e7 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -10,6 +10,11 @@ #include // For 3d rectangular grids #include +#include +#include +#include +#include + #include "DiskDensityFunc.H" #include "ParticleReader.H" #include "Coefficients.H" @@ -110,6 +115,13 @@ namespace BasisClasses Eigen::VectorXd sampleMasses; //@} + //@{ + //! HDF5 compression settings + unsigned H5compress = 5; + unsigned H5chunk = 4096; // 1024*1024; + bool H5shuffle = true; + //@} + //! Round time key to emulated fixed-point arithmetic double roundTime(double time) { @@ -128,7 +140,7 @@ namespace BasisClasses virtual void extendCoefCovariance(const std::string& filename, double time); //! Variance file versioning - inline static const std::string CovarianceFileVersion = "1.0"; + inline static const std::string CovarianceFileVersion = "1.01"; public: @@ -276,6 +288,14 @@ namespace BasisClasses //! Choose between float and double storage for covariance void setCovarFloatType(bool flt) { floatType = flt; } + //! HDF5 compression settings + void setCovarH5Compress(unsigned level, unsigned chunksize, bool shuffle) + { + H5compress = level; + H5chunk = chunksize; + H5shuffle = shuffle; + } + //! Make covariance after accumulation virtual void makeCoefCovariance(void) {} diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 4a7cdddd1..d7d0fd21d 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -4848,13 +4848,45 @@ namespace BasisClasses time = roundTime(time); stanza.createAttribute("Time", HighFive::DataSpace::From(time)).write(time); + // Enable compression + // + auto dcpl1 = HighFive::DataSetCreateProps{}; + auto dcpl2 = HighFive::DataSetCreateProps{}; + auto dcpl3 = HighFive::DataSetCreateProps{}; + + // Properties for sample data + if (H5compress) { + unsigned int csz = sampleCounts.size(); + dcpl1.add(HighFive::Chunking({csz, 1})); + if (H5shuffle) dcpl1.add(HighFive::Shuffle()); + dcpl1.add(HighFive::Deflate(H5compress)); + } + // Add the sample statistics // - HighFive::DataSet s1data = stanza.createDataSet("sampleCounts", sampleCounts); - HighFive::DataSet s2data = stanza.createDataSet("sampleMasses", sampleMasses); + HighFive::DataSet s1data = stanza.createDataSet("sampleCounts", sampleCounts, dcpl1); + HighFive::DataSet s2data = stanza.createDataSet("sampleMasses", sampleMasses, dcpl1); auto covar = getCoefCovariance(); + if (H5compress) { + // Properties for coefficients + unsigned int csz2 = std::get<0>(covar[0][0]).size(); + HighFive::Chunking data_dims2{std::min(csz2, H5chunk), 1}; + + dcpl2.add(data_dims2); + if (H5shuffle) dcpl2.add(HighFive::Shuffle()); + dcpl2.add(HighFive::Deflate(H5compress)); + + // Properties for covariance + unsigned int csz3 = std::get<1>(covar[0][0]).size(); + HighFive::Chunking data_dims3{std::min(csz3, H5chunk), 1}; + + dcpl3.add(data_dims3); + if (H5shuffle) dcpl3.add(HighFive::Shuffle()); + dcpl3.add(HighFive::Deflate(H5compress)); + } + // Add the samples // for (size_t T=0; T(); - data = dataF.cast>(); + // Get the real and imaginary parts + Eigen::VectorXf data_real = + sample.getDataSet("coefficients_real").read(); + + Eigen::VectorXf data_imag = + sample.getDataSet("coefficients_imag").read(); + + // Resize the complex array and assign + data.resize(data_real.size()); + data.real() = data_real.cast(); + data.imag() = data_imag.cast(); } else { - data = sample.getDataSet("coefficients").read(); + // Get the real and imaginary parts + Eigen::VectorXd data_real = + sample.getDataSet("coefficients_real").read(); + + Eigen::VectorXd data_imag = + sample.getDataSet("coefficients_imag").read(); + + // Resize the complex array and assign + data.resize(data_real.size()); + data.real() = data_real; + data.imag() = data_imag; } // Pack the coefficient data @@ -5172,11 +5258,31 @@ namespace BasisClasses // Get the flattened covariance array if (sz==4) { - Eigen::VectorXcf dataF; - dataF = sample.getDataSet("covariance").read(); - data = dataF.cast>(); + // Real parts + Eigen::VectorXf data_real = + sample.getDataSet("covariance_real").read(); + + // Imaginary parts + Eigen::VectorXf data_imag = + sample.getDataSet("covariance_imag").read(); + + // Resize the final storage + data.resize(data_real.size()); + data.real() = data_real.cast(); + data.imag() = data_real.cast(); } else { - data = sample.getDataSet("covariance").read(); + // Real parts + Eigen::VectorXd data_real = + sample.getDataSet("covariance_real").read(); + + // Imaginary parts + Eigen::VectorXd data_imag = + sample.getDataSet("covariance_imag").read(); + + // Resize the final storage + data.resize(data_real.size()); + data.real() = data_real; + data.imag() = data_real; } // Pack the coefficient data diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index fc5f76ab5..58fc80536 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -1432,6 +1432,23 @@ void BasisFactoryClasses(py::module &m) ------- None )") + .def("setCovarH5Compress", &BasisClasses::BiorthBasis::setCovarH5Compress, + R"( + Set the HDF5 compression level for covariance storage in HDF5. + + Parameters + ---------- + compress : int + HDF5 compression level 0-9 (default: 5) + chunkSize : int + HDF5 chunk size for dataset storage (default: 1024*1024) + shuffle : bool + Use shuffle filter if true (default: true) + + Returns + ------- + None + )", py::arg("compress")=5, py::arg("chunkSize")=1024*1024, py::arg("shuffle")=true) .def("getFields", &BasisClasses::BiorthBasis::getFields, R"( Return the field evaluations for a given cartesian position. The From c37ceaa27665c54ef49904280fb873e22802eed1 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 25 Oct 2025 23:22:29 -0400 Subject: [PATCH 50/82] Remove debug lines --- expui/BiorthBasis.cc | 6 ------ 1 file changed, 6 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index d7d0fd21d..1f92ced91 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -4953,9 +4953,6 @@ namespace BasisClasses Eigen::VectorXd real_part = data.real(); Eigen::VectorXd imag_part = data.imag(); - std::cout << "Real part size = " << real_part.size() << std::endl; - std::cout << "Imag part size = " << real_part.size() << std::endl; - // Create two separate, compressed datasets sample.createDataSet("coefficients_real", real_part, dcpl2); sample.createDataSet("coefficients_imag", imag_part, dcpl2); @@ -4977,10 +4974,7 @@ namespace BasisClasses imag_part = data.imag(); // Create two separate, compressed datasets - std::cout << "Making real coefs [double]" << std::endl; sample.createDataSet("covariance_real", real_part, dcpl3); - - std::cout << "Making imag coefs [double]" << std::endl; sample.createDataSet("covariance_imag", imag_part, dcpl3); } } From ed50705d19a5d9c69848ee03a57eed4b8ffff294 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 26 Oct 2025 11:39:59 -0400 Subject: [PATCH 51/82] Decreased granularity of HDF5 stanzas to allow the compressor to be efficient --- expui/BiorthBasis.cc | 298 +++++++++++++++++++++++-------------------- 1 file changed, 158 insertions(+), 140 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 1f92ced91..853760673 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -4850,11 +4850,11 @@ namespace BasisClasses // Enable compression // - auto dcpl1 = HighFive::DataSetCreateProps{}; - auto dcpl2 = HighFive::DataSetCreateProps{}; - auto dcpl3 = HighFive::DataSetCreateProps{}; + auto dcpl1 = HighFive::DataSetCreateProps{}; // sample stats + auto dcpl2 = HighFive::DataSetCreateProps{}; // coefficients + auto dcpl3 = HighFive::DataSetCreateProps{}; // covariance - // Properties for sample data + // Properties for sample stats if (H5compress) { unsigned int csz = sampleCounts.size(); dcpl1.add(HighFive::Chunking({csz, 1})); @@ -4867,11 +4867,32 @@ namespace BasisClasses HighFive::DataSet s1data = stanza.createDataSet("sampleCounts", sampleCounts, dcpl1); HighFive::DataSet s2data = stanza.createDataSet("sampleMasses", sampleMasses, dcpl1); + // Covariance data + // auto covar = getCoefCovariance(); + // Number of samples + // + unsigned sampleSize = covar.size(); + unsigned ltot = covar[0].size(); + unsigned nmax = std::get<0>(covar[0][0]).rows(); + unsigned diagonalSize = nmax*(nmax + 1)/2; + + // Add data dimensions + // + stanza.createAttribute + ("sampleSize", HighFive::DataSpace::From(sampleSize)).write(sampleSize); + + stanza.createAttribute + ("angularSize", HighFive::DataSpace::From(ltot)).write(ltot); + + stanza.createAttribute + ("rankSize", HighFive::DataSpace::From(nmax)).write(nmax); + if (H5compress) { // Properties for coefficients - unsigned int csz2 = std::get<0>(covar[0][0]).size(); + // + unsigned int csz2 = nmax * ltot * sampleSize; HighFive::Chunking data_dims2{std::min(csz2, H5chunk), 1}; dcpl2.add(data_dims2); @@ -4879,7 +4900,8 @@ namespace BasisClasses dcpl2.add(HighFive::Deflate(H5compress)); // Properties for covariance - unsigned int csz3 = std::get<1>(covar[0][0]).size(); + // + unsigned int csz3 = ltot * diagonalSize * sampleSize; HighFive::Chunking data_dims3{std::min(csz3, H5chunk), 1}; dcpl3.add(data_dims3); @@ -4887,44 +4909,33 @@ namespace BasisClasses dcpl3.add(HighFive::Deflate(H5compress)); } - // Add the samples + // Pack the coefficient data // - for (size_t T=0; T(covar[T][0]).rows(); + if (floatType) { + Eigen::VectorXcf data(nmax*ltot*sampleSize); - if (floatType) { - Eigen::VectorXcf data(lmax*nmax); - - for (size_t l=0, c=0; l> (std::get<0>(covar[T][l])(n)); } } + } - // Create a vector of doubles for the real and imaginary parts - Eigen::VectorXf real_part = data.real(); - Eigen::VectorXf imag_part = data.imag(); - - // Create two separate, compressed datasets - sample.createDataSet("coefficients_real", real_part, dcpl2); - sample.createDataSet("coefficients_imag", imag_part, dcpl2); + // Create a vector of doubles for the real and imaginary parts + Eigen::VectorXf real_part = data.real(); + Eigen::VectorXf imag_part = data.imag(); - // Pack the covariance data in an upper triangular format - // - size_t diagonalSize = nmax*(nmax + 1)/2; - data.resize(lmax*diagonalSize); - for (size_t l=0, c=0; l> @@ -4932,51 +4943,59 @@ namespace BasisClasses } } } + } - // Create a vector of doubles for the real and imaginary parts - real_part = data.real(); - imag_part = data.imag(); + // Create a vector of doubles for the real and imaginary parts + real_part = data.real(); + imag_part = data.imag(); - // Create two separate, compressed datasets - sample.createDataSet("covariance_real", real_part, dcpl3); - sample.createDataSet("covariance_imag", imag_part, dcpl3); + // Create two separate, compressed datasets + stanza.createDataSet("covariance_real", real_part, dcpl3); + stanza.createDataSet("covariance_imag", imag_part, dcpl3); + + } else { + Eigen::VectorXcd data(ltot*nmax*sampleSize); - } else { - Eigen::VectorXcd data(lmax*nmax); - for (size_t l=0, c=0; l(covar[T][l])(n); } } + } - // Create a vector of doubles for the real and imaginary parts - Eigen::VectorXd real_part = data.real(); - Eigen::VectorXd imag_part = data.imag(); - - // Create two separate, compressed datasets - sample.createDataSet("coefficients_real", real_part, dcpl2); - sample.createDataSet("coefficients_imag", imag_part, dcpl2); + // Create a vector of doubles for the real and imaginary parts + // + Eigen::VectorXd real_part = data.real(); + Eigen::VectorXd imag_part = data.imag(); - // Pack the covariance data in an upper triangular format - // - size_t diagonalSize = nmax*(nmax + 1)/2; - data.resize(lmax*diagonalSize); - for (size_t l=0, c=0; l(covar[T][l])(n1, n2); } } } + } - // Create a vector of doubles for the real and imaginary parts - real_part = data.real(); - imag_part = data.imag(); + // Create a vector of doubles for the real and imaginary parts + // + real_part = data.real(); + imag_part = data.imag(); - // Create two separate, compressed datasets - sample.createDataSet("covariance_real", real_part, dcpl3); - sample.createDataSet("covariance_imag", imag_part, dcpl3); - } + // Create two separate, compressed datasets + // + stanza.createDataSet("covariance_real", real_part, dcpl3); + stanza.createDataSet("covariance_imag", imag_part, dcpl3); } // END: sample loop @@ -5193,108 +5212,107 @@ namespace BasisClasses sampleMasses.push_back(Eigen::VectorXd()); stanza.getDataSet("sampleMasses").read(sampleMasses.back()); - int nT = sampleCounts.back().size(); + // Get data attributes + // + int nT, lSize, rank; + stanza.getAttribute("sampleSize") .read(nT); + stanza.getAttribute("angularSize").read(lSize); + stanza.getAttribute("rankSize") .read(rank); - // Allocate vector for current time + // Allocate sample vector for current time covarData.push_back(std::vector>(nT)); - for (int T=0; T(); - // Repack the data - std::vector elem(ltot); - for (auto & e : elem) { - std::get<0>(e).resize(nmax); - std::get<1>(e).resize(nmax, nmax); - } - - // Get the flattened coefficient array - if (sz==4) { - // Get the real and imaginary parts - Eigen::VectorXf data_real = - sample.getDataSet("coefficients_real").read(); + Eigen::VectorXf data_imag = + stanza.getDataSet("coefficients_imag").read(); + + // Resize the complex array and assign + data0.resize(data_real.size()); + data0.real() = data_real.cast(); + data0.imag() = data_imag.cast(); - Eigen::VectorXf data_imag = - sample.getDataSet("coefficients_imag").read(); + data_real = + stanza.getDataSet("covariance_real").read(); - // Resize the complex array and assign - data.resize(data_real.size()); - data.real() = data_real.cast(); - data.imag() = data_imag.cast(); - } else { - // Get the real and imaginary parts - Eigen::VectorXd data_real = - sample.getDataSet("coefficients_real").read(); + data_imag = + stanza.getDataSet("covariance_imag").read(); + + // Resize the complex array and assign + data1.resize(data_real.size()); + data1.real() = data_real.cast(); + data1.imag() = data_imag.cast(); + } else { + // Get the real and imaginary parts + Eigen::VectorXd data_real = + stanza.getDataSet("coefficients_real").read(); + + Eigen::VectorXd data_imag = + stanza.getDataSet("coefficients_imag").read(); + + // Resize the complex array and assign + data0.resize(data_real.size()); + data0.real() = data_real; + data0.imag() = data_imag; - Eigen::VectorXd data_imag = - sample.getDataSet("coefficients_imag").read(); + // Get the real and imaginary parts + data_real = + stanza.getDataSet("covariance_real").read(); - // Resize the complex array and assign - data.resize(data_real.size()); - data.real() = data_real; - data.imag() = data_imag; - } + data_imag = + stanza.getDataSet("covariance_imag").read(); - // Pack the coefficient data - for (size_t l=0, c=0; l(elem[l])(n) = data(c); - } - } + // Resize the complex array and assign + data1.resize(data_real.size()); + data1.real() = data_real; + data1.imag() = data_imag; + } - // Get the flattened covariance array - if (sz==4) { - // Real parts - Eigen::VectorXf data_real = - sample.getDataSet("covariance_real").read(); + // Positions in data stanzas + int sCof = 0, sCov = 0; - // Imaginary parts - Eigen::VectorXf data_imag = - sample.getDataSet("covariance_imag").read(); + for (int T=0; T(); - data.imag() = data_real.cast(); - } else { - // Real parts - Eigen::VectorXd data_real = - sample.getDataSet("covariance_real").read(); - - // Imaginary parts - Eigen::VectorXd data_imag = - sample.getDataSet("covariance_imag").read(); - - // Resize the final storage - data.resize(data_real.size()); - data.real() = data_real; - data.imag() = data_real; + // Repack the data + std::vector elem(lSize); + for (auto & e : elem) { + std::get<0>(e).resize(rank); + std::get<1>(e).resize(rank, rank); } - + // Pack the coefficient data - for (size_t l=0, c=0; l(elem[l])(n1, n2) = data(c); + int c = 0; + for (size_t l=0; l(elem[l])(n) = data0(c+sCof); + } + } + sCof += c; + + // Pack the covariance data + c = 0; + for (size_t l=0; l(elem[l])(n1, n2) = data1(c+sCov); if (n1 != n2) std::get<1>(elem[l])(n2, n1) = std::get<1>(elem[l])(n1, n2); } } } + sCov += c; // Add the data covarData.back()[T] = std::move(elem); } // END: sample loop - } // END: snapshot loop } catch (HighFive::Exception& err) { From e8ac21decb5dbad0d3fc32228dcb88ec832635c0 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 26 Oct 2025 12:31:56 -0400 Subject: [PATCH 52/82] Removed some debugging output; checked against original version --- expui/BiorthBasis.H | 5 ++-- expui/BiorthBasis.cc | 70 ++++++++++++++++++++------------------------ 2 files changed, 35 insertions(+), 40 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 0bab9c5e7..e0e8a246f 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -139,8 +139,9 @@ namespace BasisClasses //! Add coefficient covariance data to an HDF5 file virtual void extendCoefCovariance(const std::string& filename, double time); - //! Variance file versioning - inline static const std::string CovarianceFileVersion = "1.01"; + //! Variance file versioning (Version 1.0 has more granularity and + //! poor compression) + inline static const std::string CovarianceFileVersion = "1.1"; public: diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 853760673..e819ddb0b 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -4912,63 +4912,56 @@ namespace BasisClasses // Pack the coefficient data // if (floatType) { - Eigen::VectorXcf data(nmax*ltot*sampleSize); + // Create a vector of doubles for the real and imaginary parts + Eigen::VectorXf real_part(nmax*ltot*sampleSize); + Eigen::VectorXf imag_part(nmax*ltot*sampleSize); for (size_t T=0, c=0; T> - (std::get<0>(covar[T][l])(n)); + real_part(c) = std::real(std::get<0>(covar[T][l])(n)); + imag_part(c) = std::imag(std::get<0>(covar[T][l])(n)); } } } - // Create a vector of doubles for the real and imaginary parts - Eigen::VectorXf real_part = data.real(); - Eigen::VectorXf imag_part = data.imag(); - // Create two separate, compressed datasets stanza.createDataSet("coefficients_real", real_part, dcpl2); stanza.createDataSet("coefficients_imag", imag_part, dcpl2); // Pack the covariance data in an upper triangular format // - data.resize(ltot*diagonalSize*sampleSize); + real_part.resize(ltot*diagonalSize*sampleSize); + imag_part.resize(ltot*diagonalSize*sampleSize); + for (size_t T=0, c=0; T> - (std::get<1>(covar[T][l])(n1, n2)); + real_part(c) = std::real(std::get<1>(covar[T][l])(n1, n2)); + imag_part(c) = std::imag(std::get<1>(covar[T][l])(n1, n2)); } } } } - // Create a vector of doubles for the real and imaginary parts - real_part = data.real(); - imag_part = data.imag(); - // Create two separate, compressed datasets stanza.createDataSet("covariance_real", real_part, dcpl3); stanza.createDataSet("covariance_imag", imag_part, dcpl3); } else { - Eigen::VectorXcd data(ltot*nmax*sampleSize); + Eigen::VectorXd real_part(ltot*nmax*sampleSize); + Eigen::VectorXd imag_part(ltot*nmax*sampleSize); - for (size_t T=0; T(covar[T][l])(n); + real_part(c) = std::real(std::get<0>(covar[T][l])(n)); + imag_part(c) = std::imag(std::get<0>(covar[T][l])(n)); } } } - // Create a vector of doubles for the real and imaginary parts - // - Eigen::VectorXd real_part = data.real(); - Eigen::VectorXd imag_part = data.imag(); - // Create two separate, compressed datasets // stanza.createDataSet("coefficients_real", real_part, dcpl2); @@ -4976,22 +4969,20 @@ namespace BasisClasses // Pack the covariance data in an upper triangular format // - data.resize(ltot*diagonalSize*sampleSize); - for (size_t T=0; T(covar[T][l])(n1, n2); + real_part(c) = std::real(std::get<1>(covar[T][l])(n1, n2)); + imag_part(c) = std::imag(std::get<1>(covar[T][l])(n1, n2)); } } } } - // Create a vector of doubles for the real and imaginary parts - // - real_part = data.real(); - imag_part = data.imag(); - // Create two separate, compressed datasets // stanza.createDataSet("covariance_real", real_part, dcpl3); @@ -5144,7 +5135,7 @@ namespace BasisClasses throw std::runtime_error("CovarianceReader: this is an early alpha test version. Please remake your files"); } // Test for current version - if (version != std::string("1.01")) { + if (version != std::string("1.1")) { throw std::runtime_error(std::string("CovarianceReader: unsupported file version, ") + version); } @@ -5278,20 +5269,23 @@ namespace BasisClasses // Positions in data stanzas int sCof = 0, sCov = 0; + // Loop through all indices and repack for (int T=0; T elem(lSize); for (auto & e : elem) { + // Coefficients std::get<0>(e).resize(rank); + // Covariance matrix std::get<1>(e).resize(rank, rank); } // Pack the coefficient data int c = 0; for (size_t l=0; l(elem[l])(n) = data0(c+sCof); + for (size_t n=0; n(elem[l])(n) = data0(sCof + c++); } } sCof += c; @@ -5300,8 +5294,8 @@ namespace BasisClasses c = 0; for (size_t l=0; l(elem[l])(n1, n2) = data1(c+sCov); + for (size_t n2=n1; n2(elem[l])(n1, n2) = data1(sCov + c++); if (n1 != n2) std::get<1>(elem[l])(n2, n1) = std::get<1>(elem[l])(n1, n2); } From 7bdd491e7d91fd759d62707afae72088de9cdd33 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 26 Oct 2025 12:34:46 -0400 Subject: [PATCH 53/82] Restore 1MB as the default chunk size --- expui/BiorthBasis.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index e0e8a246f..aa2ce30b9 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -118,7 +118,7 @@ namespace BasisClasses //@{ //! HDF5 compression settings unsigned H5compress = 5; - unsigned H5chunk = 4096; // 1024*1024; + unsigned H5chunk = 1024*1024; // (1 MB) bool H5shuffle = true; //@} From d2b1eccefb271da021b1d3b535d287d12744818d Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 26 Oct 2025 14:06:57 -0400 Subject: [PATCH 54/82] Added Python doc strings --- expui/BiorthBasis.H | 4 +++- expui/BiorthBasis.cc | 16 ++++++++++++++-- pyEXP/BasisWrappers.cc | 11 ++++++++--- 3 files changed, 25 insertions(+), 6 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index aa2ce30b9..8cea78905 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -120,6 +120,7 @@ namespace BasisClasses unsigned H5compress = 5; unsigned H5chunk = 1024*1024; // (1 MB) bool H5shuffle = true; + bool H5szip = false; //@} //! Round time key to emulated fixed-point arithmetic @@ -290,11 +291,12 @@ namespace BasisClasses void setCovarFloatType(bool flt) { floatType = flt; } //! HDF5 compression settings - void setCovarH5Compress(unsigned level, unsigned chunksize, bool shuffle) + void setCovarH5Compress(unsigned level, unsigned chunksize, bool shuffle, bool szip=false) { H5compress = level; H5chunk = chunksize; H5shuffle = shuffle; + H5szip = szip; } //! Make covariance after accumulation diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index e819ddb0b..e002a1b57 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -4890,6 +4890,10 @@ namespace BasisClasses ("rankSize", HighFive::DataSpace::From(nmax)).write(nmax); if (H5compress) { + // Szip parameters + const int options_mask = H5_SZIP_NN_OPTION_MASK; + const int pixels_per_block = 8; + // Properties for coefficients // unsigned int csz2 = nmax * ltot * sampleSize; @@ -4897,7 +4901,11 @@ namespace BasisClasses dcpl2.add(data_dims2); if (H5shuffle) dcpl2.add(HighFive::Shuffle()); - dcpl2.add(HighFive::Deflate(H5compress)); + if (H5szip) { + dcpl2.add(HighFive::Szip(options_mask, pixels_per_block)); + } else { + dcpl2.add(HighFive::Deflate(H5compress)); + } // Properties for covariance // @@ -4906,7 +4914,11 @@ namespace BasisClasses dcpl3.add(data_dims3); if (H5shuffle) dcpl3.add(HighFive::Shuffle()); - dcpl3.add(HighFive::Deflate(H5compress)); + if (H5szip) { + dcpl3.add(HighFive::Szip(options_mask, pixels_per_block)); + } else { + dcpl3.add(HighFive::Deflate(H5compress)); + } } // Pack the coefficient data diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 58fc80536..3549aeae9 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -1434,21 +1434,26 @@ void BasisFactoryClasses(py::module &m) )") .def("setCovarH5Compress", &BasisClasses::BiorthBasis::setCovarH5Compress, R"( - Set the HDF5 compression level for covariance storage in HDF5. + Set the HDF5 compression level for covariance storage in HDF5. The Szip + compression algorithm may also be enabled but seems to not have better + performance for these data than the standard Gzip compression. Gzip with + level 5 and shuffling is enabled by default. Parameters ---------- compress : int - HDF5 compression level 0-9 (default: 5) + HDF5 Gzip compression level 0-9 (default: 5) chunkSize : int HDF5 chunk size for dataset storage (default: 1024*1024) shuffle : bool Use shuffle filter if true (default: true) + szip : bool + Use Szip compression algorithm (default: false) Returns ------- None - )", py::arg("compress")=5, py::arg("chunkSize")=1024*1024, py::arg("shuffle")=true) + )", py::arg("compress")=5, py::arg("chunkSize")=1024*1024, py::arg("shuffle")=true, py::arg("azip")=false) .def("getFields", &BasisClasses::BiorthBasis::getFields, R"( Return the field evaluations for a given cartesian position. The From 2c919807553ff3d2c1a7475da9927754209f2163 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 27 Oct 2025 18:31:53 -0400 Subject: [PATCH 55/82] Fixed loop misorder in complex covariance accumulation --- exputil/EmpCylSL.cc | 42 +++++++++++++++++++++++++----------------- 1 file changed, 25 insertions(+), 17 deletions(-) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 079351659..fb315af7c 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -1946,7 +1946,7 @@ void EmpCylSL::setup_accumulation(int mlevel) if (m>0) accum_sin[m].setZero(); } - if ( (PCAVAR or PCAEOF) and mlevel==0 and sampT>1) { + if ( (covar or PCAVAR or PCAEOF) and mlevel==0 and sampT>0) { for (int nth=0; nth(mcos, msin) * vc[id].row(mm).transpose() * norm; - VC[id][whch][mm] += mass * vec; - MV[id][whch][mm] += mass * (vec * vec.adjoint()); - } - if (compute and PCAVAR) { double hc1 = vc[id](mm, nn)*mcos, hs1 = 0.0; if (mm) hs1 = vs[id](mm, nn)*msin; @@ -4107,9 +4102,17 @@ void EmpCylSL::accumulate(double r, double z, double phi, double mass, } } + if (compute and covar) { + Eigen::VectorXcd vec = std::complex(mcos, msin) * + vc[id].row(mm).transpose() * norm; + + VC[id][whch][mm] += mass * vec; + MV[id][whch][mm] += mass * (vec * vec.adjoint()); + } + cylmass1[id] += mass; } - + // End: M loop } @@ -4352,6 +4355,10 @@ void EmpCylSL::make_coefficients(bool compute) MPIin2 .resize(rank3*rank3*(MMAX+1)); MPIout2.resize(rank3*rank3*(MMAX+1)); } + + std::cout << "Process " << myid << ": completed coefficient accumulation, " + << "multistep=" << multistep << ", sampT=" << sampT << std::endl; + // Sum up over threads // for (unsigned M=0; M<=multistep; M++) { @@ -4363,6 +4370,7 @@ void EmpCylSL::make_coefficients(bool compute) howmany1[M][0] += howmany1[M][nth]; if (compute and (covar or PCAVAR)) { + for (unsigned T=0; T(ret).resize(sampT); std::get<1>(ret).resize(sampT); + std::get<0>(ret).setZero(); + std::get<1>(ret).setZero(); for (unsigned T=0; T(ret)(T) = numbT[T]; std::get<1>(ret)(T) = massT[T]; From 4dc27cc076352dbd36dfaa3bfc2876e4df37241c Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 28 Oct 2025 15:14:49 -0400 Subject: [PATCH 56/82] Remove a noisy DEBUG line --- exputil/EmpCylSL.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index fb315af7c..628a82224 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -1962,7 +1962,6 @@ void EmpCylSL::setup_accumulation(int mlevel) } if (covar) { - std::cout << "SANITY: cleaning up" << std::endl; for (unsigned T=0; T Date: Wed, 29 Oct 2025 15:44:05 -0400 Subject: [PATCH 57/82] More debug output removal --- exputil/EmpCylSL.cc | 3 --- 1 file changed, 3 deletions(-) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 628a82224..af33cf7b4 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -4355,9 +4355,6 @@ void EmpCylSL::make_coefficients(bool compute) MPIout2.resize(rank3*rank3*(MMAX+1)); } - std::cout << "Process " << myid << ": completed coefficient accumulation, " - << "multistep=" << multistep << ", sampT=" << sampT << std::endl; - // Sum up over threads // for (unsigned M=0; M<=multistep; M++) { From e4e032fba904aca904ec4d82d0760e3242614b7e Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 30 Oct 2025 11:04:27 -0400 Subject: [PATCH 58/82] Remove some unused Cylindrical parameters and update the attribute output --- expui/BiorthBasis.H | 7 +------ expui/BiorthBasis.cc | 5 ++--- 2 files changed, 3 insertions(+), 9 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 8cea78905..3b8c2c30b 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -910,8 +910,7 @@ namespace BasisClasses int lmaxfid, nmaxfid, mmax, mlim, nmax; int ncylodd, ncylnx, ncylny, ncylr, cmap, cmapR, cmapZ, vflag; int rnum, pnum, tnum; - double rmin, rmax, rcylmin, rcylmax; - double acyl, hcyl; + double rcylmin, rcylmax, acyl, hcyl; bool expcond, logarithmic, density, EVEN_M, sech2 = false; std::vector potd, dpot, dpt2, dend; @@ -920,7 +919,6 @@ namespace BasisClasses bool coefs_defined; Eigen::MatrixXd factorial; Eigen::MatrixXd expcoef; - double rscl; int used; std::string cachename; @@ -1008,9 +1006,6 @@ namespace BasisClasses return EmpCylSL::cacheInfo(cachefile); } - //! Prescaling factor - void set_scale(const double scl) { rscl = scl; } - //! Zero out coefficients to prepare for a new expansion void reset_coefs(void); diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index e002a1b57..2f0bfb919 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -4827,9 +4827,8 @@ namespace BasisClasses { file.createAttribute("mmax", HighFive::DataSpace::From(mmax)).write(mmax); file.createAttribute("nmax", HighFive::DataSpace::From(nmax)).write(nmax); - file.createAttribute("scale", HighFive::DataSpace::From(rscl)).write(rscl); - file.createAttribute("rmin", HighFive::DataSpace::From(rmin)).write(rmin); - file.createAttribute("rmax", HighFive::DataSpace::From(rmax)).write(rmax); + file.createAttribute("rcylmin", HighFive::DataSpace::From(rcylmin)).write(rcylmin); + file.createAttribute("rcylmax", HighFive::DataSpace::From(rcylmax)).write(rcylmax); file.createAttribute("acyl", HighFive::DataSpace::From(acyl)).write(acyl); file.createAttribute("hcyl", HighFive::DataSpace::From(hcyl)).write(hcyl); } From 968090a4e02c344b7e4b750ed5fc6c6d3e5ec1d8 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 2 Nov 2025 12:14:49 -0500 Subject: [PATCH 59/82] Add KD density estimation utility --- expui/CMakeLists.txt | 2 +- expui/KDdensity.H | 55 ++++++++++++++++++ expui/KDdensity.cc | 126 ++++++++++++++++++++++++++++++++++++++++++ pyEXP/UtilWrappers.cc | 55 ++++++++++++++++++ 4 files changed, 237 insertions(+), 1 deletion(-) create mode 100644 expui/KDdensity.H create mode 100644 expui/KDdensity.cc diff --git a/expui/CMakeLists.txt b/expui/CMakeLists.txt index 30102d192..171003115 100644 --- a/expui/CMakeLists.txt +++ b/expui/CMakeLists.txt @@ -38,7 +38,7 @@ endif() set(expui_SOURCES BasisFactory.cc BiorthBasis.cc FieldBasis.cc CoefContainer.cc CoefStruct.cc FieldGenerator.cc expMSSA.cc Coefficients.cc KMeans.cc Centering.cc ParticleIterator.cc - Koopman.cc BiorthBess.cc SvdSignChoice.cc) + Koopman.cc BiorthBess.cc SvdSignChoice.cc KDdensity.cc) add_library(expui ${expui_SOURCES}) set_target_properties(expui PROPERTIES OUTPUT_NAME expui) target_include_directories(expui PUBLIC ${common_INCLUDE}) diff --git a/expui/KDdensity.H b/expui/KDdensity.H new file mode 100644 index 000000000..46c812178 --- /dev/null +++ b/expui/KDdensity.H @@ -0,0 +1,55 @@ +#ifndef _KDdensity_H_ +#define _KDdensity_H_ + +#include +#include + +#include + +#include "ParticleReader.H" +#include "permutation.H" +#include "localmpi.H" +#include "KDtree.H" + +namespace Utility +{ + class KDdensity + { + protected: + using point3 = KDtree::point ; + using tree3 = KDtree::kdtree; + + //! Point list + std::vector points; + + //! Density array + std::vector KDdens; + + //! Accumulated mass of phase space + double KDmass = 0.0; + + //! Number of points per ball + int Ndens = 32; + + //! The KD tree instance + std::shared_ptr kdtree_; + + //! The point index to particle ID map + std::map indx; + + public: + using RowMatrixXd = Eigen::Matrix; + + //! Create KD tree from a particle reader + KDdensity(PR::PRptr reader, int Ndens=32); + + //! Create KD tree from a arrays of masses and positions + KDdensity(const Eigen::VectorXd& m, const RowMatrixXd& d, int Ndens=32); + + //! Get density for particle index + double getDensity(unsigned long index); + + }; +} + +#endif diff --git a/expui/KDdensity.cc b/expui/KDdensity.cc new file mode 100644 index 000000000..7b995ecee --- /dev/null +++ b/expui/KDdensity.cc @@ -0,0 +1,126 @@ +#include "KDdensity.H" + +namespace Utility +{ + KDdensity::KDdensity(PR::PRptr reader, int Ndens) : Ndens(Ndens) + { + // Check for MPI and initialize if needed + // + int flag; + MPI_Initialized(&flag); + bool mpi_enabled = (flag != 0); + + // Get number of particles + // + int nbod = reader->CurrentNumber(); + + // Every node needs to make the tree (a parallel share could be + // implemented in KDtree.H). Build the ID map. + // + KDmass = 0.0; + for (auto part=reader->firstParticle(); part!=0; part=reader->nextParticle()) { + KDmass += part->mass; + points.push_back(point3({part->pos[0], part->pos[1], part->pos[2]}, part->mass)); + indx[part->indx] = points.size()-1; + } + + // Build the k-d tree + // + kdtree_ = std::make_shared(points.begin(), points.end()); + + + // Share the density computation among the nodes + // + KDdens.resize(nbod, 0.0); + unsigned long badVol = 0; + for (int k=0; knearestN(points[k], Ndens); + double volume = 4.0*M_PI/3.0*std::pow(std::get<2>(ret), 3.0); + if (volume>0.0 and KDmass>0.0) + KDdens[k] = std::get<1>(ret)/volume/KDmass; + else badVol++; + } + } + + if (mpi_enabled) + MPI_Allreduce(MPI_IN_PLACE, KDdens.data(), nbod, + MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + + if (myid==0) { + std::cout << "---- KDdensity: finished density estimate with " << badVol + << " undetermined densities, mass=" << KDmass << std::endl; + std::cout << "---- KDdensity: a few densities are " << KDdens.front() + << ", " << KDdens[nbod/2] << ", " << KDdens.back() << std::endl; + std::cout << "---- KDdensity: bad density count is " << badVol + << " of " << nbod << " particles." << std::endl; + } + } + + KDdensity::KDdensity(const Eigen::VectorXd& m, + const KDdensity::RowMatrixXd& pos, + int Ndens) : Ndens(Ndens) + { + // Check for MPI and initialize if needed + // + int flag; + MPI_Initialized(&flag); + bool mpi_enabled = (flag != 0); + + // Get number of particles + // + int nbod = m.size(); + + // Sanity check + // + if (pos.rows() != nbod or pos.cols() != 3) + throw std::runtime_error("KDdensity: position array has wrong dimensions"); + + // Every node needs to make the tree (a parallel share could be + // implemented in KDtree.H). Build the ID map. + // + KDmass = 0.0; + for (int i=0; i(points.begin(), points.end()); + + + // Share the density computation among the nodes + // + KDdens.resize(nbod, 0.0); + unsigned long badVol = 0; + for (int k=0; knearestN(points[k], Ndens); + double volume = 4.0*M_PI/3.0*std::pow(std::get<2>(ret), 3.0); + if (volume>0.0 and KDmass>0.0) + KDdens[k] = std::get<1>(ret)/volume/KDmass; + else badVol++; + } + } + + if (mpi_enabled) + MPI_Allreduce(MPI_IN_PLACE, KDdens.data(), nbod, + MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + + if (myid==0) { + std::cout << "---- KDdensity: finished density estimate with " << badVol + << " undetermined densities, mass=" << KDmass << std::endl; + std::cout << "---- KDdensity: a few densities are " << KDdens.front() + << ", " << KDdens[nbod/2] << ", " << KDdens.back() << std::endl; + std::cout << "---- KDdensity: bad density count is " << badVol + << " of " << nbod << " particles." << std::endl; + } + } + + double KDdensity::getDensity(unsigned long index) + { + return KDdens[indx[index]]; + } +} diff --git a/pyEXP/UtilWrappers.cc b/pyEXP/UtilWrappers.cc index 2dc43a195..682809e7c 100644 --- a/pyEXP/UtilWrappers.cc +++ b/pyEXP/UtilWrappers.cc @@ -6,6 +6,7 @@ namespace py = pybind11; #include "Centering.H" +#include "KDdensity.H" #include "ParticleIterator.H" void UtilityClasses(py::module &m) { @@ -243,4 +244,58 @@ void UtilityClasses(py::module &m) { None )", py::arg("verbose") = false); + + py::class_(m, "KDdensity") + .def(py::init(), + R"( + Create a k-d tree density estimator for a particle set + + Parameters + ---------- + reader : ParticleReader + the particle-reader class instance + Ndens : int, default=32 + number of particles per sample ball (32 is a good + compromise between accuracy and runtime; 16 is okay if + you are trying to shave off runtime. + + Returns + ------- + KDdensity + the KDdensity instance + )", py::arg("reader"), py::arg("Ndens")=32) + .def(py::init(), + R"( + Create a k-d tree density estimator for a particle set + + Parameters + ---------- + mass : ndarray(N) + mass array + pos : ndarray(N,3) + x, y, z array + Ndens : int, default=32 + number of particles per sample ball (32 is a good + compromise between accuracy and runtime; 16 is okay if + you are trying to shave off runtime. + + Returns + ------- + KDdensity + the KDdensity instance + )", py::arg("mass"), py::arg("pos"), py::arg("Ndens")=32) + .def("getDensity", &Utility::KDdensity::getDensity, + R"( + Get the density estimate at a given particle index + + Parameters + ---------- + index : int + the particle index + + Returns + ------- + float + the density estimate at the particle position + )", py::arg("index")); } From cfad2a5bf633e7f61c790f45d4aa82f7c6db44c7 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 2 Nov 2025 12:15:06 -0500 Subject: [PATCH 60/82] Fixed a few errors in the PSPhdf5 version of the ParticleReader --- exputil/ParticleReader.cc | 2 +- include/ParticleReader.H | 168 ++++++++++++++++---------------- pyEXP/ParticleReaderWrappers.cc | 123 ++++++++++++++++++++--- 3 files changed, 195 insertions(+), 98 deletions(-) diff --git a/exputil/ParticleReader.cc b/exputil/ParticleReader.cc index 52a33dd0c..857ba367b 100644 --- a/exputil/ParticleReader.cc +++ b/exputil/ParticleReader.cc @@ -705,7 +705,7 @@ namespace PR { } - PSPhdf5::PSPhdf5(const std::vector& files, bool verbose) : ParticleReader() + PSPhdf5::PSPhdf5(const std::vector& files, bool verbose) : PSP(verbose) { _files = files; _verbose = verbose; diff --git a/include/ParticleReader.H b/include/ParticleReader.H index b96986446..b57322dfa 100644 --- a/include/ParticleReader.H +++ b/include/ParticleReader.H @@ -234,90 +234,6 @@ namespace PR }; - class PSPhdf5 : public ParticleReader - { - protected: - - //@{ - //! Parameters - double time; - std::vector particles; - bool _verbose; - std::vector _files; - - std::vector mass; - std::vector Niattrib, Ndattrib; - std::vector npart, nptot; - unsigned long totalCount; - bool real4, gadget4; - int ntypes, nfiles; - //@} - - //! Current file - std::vector::iterator curfile; - - //! Components - std::vector comps; - - //@{ - //! Current component and index - std::string curcomp; - int curindx; - //@} - - unsigned pcount; - template void read_and_load(); - template void read_and_load_gadget4(); - template void read_and_load_psp(); - - void getInfo(); - void packParticle(); - bool nextFile(); - - public: - - //! Constructor - PSPhdf5(const std::vector& file, bool verbose=false); - - //! Select a particular particle type and reset the iterator - virtual void SelectType(const std::string& type) { - auto it = std::find(comps.begin(), comps.end(), type); - if (it != comps.end()) { - curcomp = type; - curindx = it - comps.begin(); - } - else { - std::cerr << "PSPhdf5 error: could not find particle component <" << type << ">" << std::endl; - std::cerr << "Available particle components are:"; - for (auto s : comps) std::cerr << " " << s; - std::cerr << std::endl; - throw std::runtime_error("PSPhdf5: non-existent component"); - } - - curfile = _files.begin(); // Set to first file and open - nextFile(); - } - - //! Number of particles in the chosen type - virtual unsigned long CurrentNumber() { return totalCount; } - - //! Return list of particle types - virtual std::vector GetTypes() { return comps; } - - //! Get current time - virtual double CurrentTime() { return time; } - - //! Reset to beginning of particles for this component - virtual const Particle* firstParticle(); - - //! Get the next particle - virtual const Particle* nextParticle(); - - //! Get file cound per snapshot - virtual int NumFiles() { return nfiles; } - - }; - class PSPstanza { public: @@ -646,6 +562,90 @@ namespace PR //@} }; + class PSPhdf5 : public PSP + { + protected: + + //@{ + //! Parameters + double time; + std::vector particles; + bool _verbose; + std::vector _files; + + std::vector mass; + std::vector Niattrib, Ndattrib; + std::vector npart, nptot; + unsigned long totalCount; + bool real4, gadget4; + int ntypes, nfiles; + //@} + + //! Current file + std::vector::iterator curfile; + + //! Components + std::vector comps; + + //@{ + //! Current component and index + std::string curcomp; + int curindx; + //@} + + unsigned pcount; + template void read_and_load(); + template void read_and_load_gadget4(); + template void read_and_load_psp(); + + void getInfo(); + void packParticle(); + bool nextFile(); + + public: + + //! Constructor + PSPhdf5(const std::vector& file, bool verbose=false); + + //! Select a particular particle type and reset the iterator + virtual void SelectType(const std::string& type) { + auto it = std::find(comps.begin(), comps.end(), type); + if (it != comps.end()) { + curcomp = type; + curindx = it - comps.begin(); + } + else { + std::cerr << "PSPhdf5 error: could not find particle component <" << type << ">" << std::endl; + std::cerr << "Available particle components are:"; + for (auto s : comps) std::cerr << " " << s; + std::cerr << std::endl; + throw std::runtime_error("PSPhdf5: non-existent component"); + } + + curfile = _files.begin(); // Set to first file and open + nextFile(); + } + + //! Number of particles in the chosen type + virtual unsigned long CurrentNumber() { return nptot[curindx]; } + + //! Return list of particle types + virtual std::vector GetTypes() { return comps; } + + //! Get current time + virtual double CurrentTime() { return time; } + + //! Reset to beginning of particles for this component + virtual const Particle* firstParticle(); + + //! Get the next particle + virtual const Particle* nextParticle(); + + //! Get file cound per snapshot + virtual int NumFiles() { return nfiles; } + + }; + typedef std::shared_ptr PRptr; } diff --git a/pyEXP/ParticleReaderWrappers.cc b/pyEXP/ParticleReaderWrappers.cc index 2e01b46d5..770d01aeb 100644 --- a/pyEXP/ParticleReaderWrappers.cc +++ b/pyEXP/ParticleReaderWrappers.cc @@ -15,11 +15,12 @@ void ParticleReaderClasses(py::module &m) { "The available particle readers are:\n" " 1. PSPout The monolithic EXP phase-space snapshot format\n" " 2. PSPspl Like PSPout, but split into multiple file chunks\n" - " 3. GadgetNative The original Gadget native format\n" - " 4 GadgetHDF5 The newer HDF5 Gadget format\n" - " 5. TipsyNative The original Tipsy format\n" - " 6. TipsyXDR The original XDR Tipsy format\n" - " 7. Bonsai This is the Bonsai varient of Tipsy files\n\n" + " 3. PSPhdf5 The HDF5 version of the EXP phase-space format\n" + " 4. GadgetNative The original Gadget native format\n" + " 5 GadgetHDF5 The newer HDF5 Gadget format\n" + " 6. TipsyNative The original Tipsy format\n" + " 7. TipsyXDR The original XDR Tipsy format\n" + " 8. Bonsai This is the Bonsai varient of Tipsy files\n\n" "We have a helper function, getReaders, to get a list to help you\n" "remember. Try: pyEXP.read.ParticleReader.getReaders()\n\n" "Each reader can manage snapshots split into many files by parallel,\n" @@ -216,6 +217,38 @@ void ParticleReaderClasses(py::module &m) { } }; + class PyPSPhdf5 : public PSPhdf5 + { + public: + + // Inherit the constructors + using PSPhdf5::PSPhdf5; + + void SelectType(const std::string& type) override { + PYBIND11_OVERRIDE(void, PSPhdf5, SelectType, type); + } + + std::vector GetTypes() override { + PYBIND11_OVERRIDE(std::vector, PSPhdf5, GetTypes,); + } + + unsigned long CurrentNumber() override { + PYBIND11_OVERRIDE(unsigned long, PSPhdf5, CurrentNumber,); + } + + double CurrentTime() override { + PYBIND11_OVERRIDE(double, PSPhdf5, CurrentTime,); + } + + const Particle* firstParticle() override { + PYBIND11_OVERRIDE(const Particle*, PSPhdf5, firstParticle,); + } + + const Particle* nextParticle() override { + PYBIND11_OVERRIDE(const Particle*, PSPhdf5, nextParticle,); + } + }; + class PyTipsy : public Tipsy { public: @@ -401,8 +434,8 @@ void ParticleReaderClasses(py::module &m) { pr.def_static("getReaders", []() { const std::vector formats = { - "PSPout", "PSPspl", "GadgetNative", "GadgetHDF5", "TipsyNative", - "TipsyXDR", "Bonsai"}; + "PSPout", "PSPspl", "PSPhdf5", "GadgetNative", "GadgetHDF5", + "TipsyNative", "TipsyXDR", "Bonsai"}; return formats; }, @@ -431,7 +464,7 @@ void ParticleReaderClasses(py::module &m) { Returns ------- ParticleReader - )"); + )", py::arg("files"), py::arg("verbose")=false); py::class_, PyGadgetNative, ParticleReader>(m, "GadgetNative") .def(py::init&, bool>(), @@ -448,7 +481,7 @@ void ParticleReaderClasses(py::module &m) { Returns ------- ParticleReader - )"); + )", py::arg("files"), py::arg("verbose")=false); py::class_, PyPSP, ParticleReader>(m, "PSP") .def(py::init(), "Base class for PSP reader") @@ -464,7 +497,20 @@ void ParticleReaderClasses(py::module &m) { py::class_, PyPSPout, PSP>(m, "PSPout") .def(py::init&, bool>(), - "Reader for monolitic PSP format (single file written by root process)") + R"( + Read PSP ascii monolithic format snapshot files + + Parameters + ---------- + files : list(str) + List of files with phase-space segments comprising a single snapshot + verbose : bool, default=False + Verbose, diagnostic output + + Returns + ------- + ParticleReader + )", py::arg("files"), py::arg("verbose")=false) .def("CurrentNumber", &PSPout::CurrentNumber) .def("GetTypes", &PSPout::GetTypes) .def("CurrentTime", &PSPout::CurrentTime); @@ -472,11 +518,46 @@ void ParticleReaderClasses(py::module &m) { py::class_, PyPSPspl, PSP>(m, "PSPspl") .def(py::init&, bool>(), - "Reader for split PSP files (multiple files written by each process)") + R"( + Reader for split PSP files (multiple files written by each process) + + Parameters + ---------- + files : list(str) + List of files with phase-space segments comprising a single snapshot + verbose : bool, default=False + Verbose, diagnostic output + + Returns + ------- + ParticleReader + )", py::arg("files"), py::arg("verbose")=false) .def("CurrentNumber", &PSPspl::CurrentNumber) .def("GetTypes", &PSPspl::GetTypes) .def("CurrentTime", &PSPspl::CurrentTime); + py::class_, PyPSPhdf5, PSP>(m, "PSPhdf5") + .def(py::init&, bool>(), + R"( + Reader for HDF5 PSP files (both single and multiple files) + + Parameters + ---------- + files : list(str) + List of files with phase-space segments comprising a single snapshot + verbose : bool, default=False + Verbose, diagnostic output + + Returns + ------- + ParticleReader + )", py::arg("files"), py::arg("verbose")=false) + .def("SelectType", &PSPhdf5::SelectType) + .def("NumFiles", &PSPhdf5::NumFiles) + .def("CurrentNumber", &PSPhdf5::CurrentNumber) + .def("GetTypes", &PSPhdf5::GetTypes) + .def("CurrentTime", &PSPhdf5::CurrentTime); + py::class_, PyTipsy, ParticleReader> tipsy(m, "Tipsy"); py::enum_(tipsy, "TipsyType") @@ -485,7 +566,24 @@ void ParticleReaderClasses(py::module &m) { .value("bonsai", Tipsy::TipsyType::bonsai) .export_values(); - tipsy.def(py::init()) + tipsy.def(py::init(), + R"( + Read Tipsy format snapshots + + Parameters + ---------- + file : str + The Tipsy snapshot file + type : TipsyType + The Tipsy file type (native, xdr, bonsai) + verbose : bool, default=False + Verbose, diagnostic output + + Returns + ------- + ParticleReader + + )", py::arg("file"), py::arg("type"), py::arg("verbose")=false) .def(py::init&, Tipsy::TipsyType, bool>()) .def("SelectType", &Tipsy::SelectType) .def("CurrentNumber", &Tipsy::CurrentNumber) @@ -493,4 +591,3 @@ void ParticleReaderClasses(py::module &m) { .def("CurrentTime", &Tipsy::CurrentTime); } - From ffd89df46b2d02529c918deb0abcfa568ceaed8e Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 2 Nov 2025 22:21:25 -0500 Subject: [PATCH 61/82] Added defaults for 'verbose' --- pyEXP/CoefWrappers.cc | 40 ++++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/pyEXP/CoefWrappers.cc b/pyEXP/CoefWrappers.cc index d4f273d2c..dad2d54c4 100644 --- a/pyEXP/CoefWrappers.cc +++ b/pyEXP/CoefWrappers.cc @@ -1053,7 +1053,7 @@ void CoefficientClasses(py::module &m) { type : str type of coefficient container verbose : bool - display verbose information. + display verbose information (default=False) Returns ------- @@ -1066,7 +1066,7 @@ void CoefficientClasses(py::module &m) { simulation )", py::arg("type"), - py::arg("verbose")) + py::arg("verbose")=false) .def("__call__", &CoefClasses::Coefs::getData, R"( @@ -1358,12 +1358,12 @@ void CoefficientClasses(py::module &m) { Parameters ---------- verbose : bool - display verbose information. + display verbose information (default=False) Returns ------- SphCoefs instance - )") + )", py::arg("verbose")=false) .def("__call__", &CoefClasses::SphCoefs::getMatrix, R"( @@ -1433,12 +1433,12 @@ void CoefficientClasses(py::module &m) { Parameters ---------- verbose : bool - display verbose information. + display verbose information (default=False) Returns ------- CylCoefs instance - )") + )", py::arg("verbose")=false) .def("__call__", &CoefClasses::CylCoefs::getMatrix, R"( @@ -1535,12 +1535,12 @@ void CoefficientClasses(py::module &m) { Parameters ---------- verbose : bool - display verbose information. + display verbose information (default=False) Returns ------- SphFldCoefs instance - )") + )", py::arg("verbose")=false) .def("__call__", [](CoefClasses::SphFldCoefs& A, double time) { @@ -1620,12 +1620,12 @@ void CoefficientClasses(py::module &m) { Parameters ---------- verbose : bool - display verbose information. + display verbose information (default=False) Returns ------- CylFldCoefs instance - )") + )", py::arg("verbose")=false) .def("__call__", [](CoefClasses::CylFldCoefs& A, double time) { @@ -1699,12 +1699,12 @@ void CoefficientClasses(py::module &m) { Parameters ---------- verbose : bool - display verbose information. + display verbose information (default=False) Returns ------- SlabCoefs instance - )") + )", py::arg("verbose")=false) .def("__call__", &CoefClasses::SlabCoefs::getTensor, R"( @@ -1793,12 +1793,12 @@ void CoefficientClasses(py::module &m) { Parameters ---------- verbose : bool - display verbose information. + display verbose information (default=False) Returns ------- CubeCoefs instance - )") + )", py::arg("verbose")=false) .def("__call__", &CoefClasses::CubeCoefs::getTensor, R"( @@ -1886,7 +1886,7 @@ void CoefficientClasses(py::module &m) { Parameters ---------- verbose : bool - display verbose information. + display verbose information (default=True) Returns ------- @@ -1914,7 +1914,7 @@ void CoefficientClasses(py::module &m) { type : str ascii table data file verbose : bool - display verbose information. + display verbose information (default=True) Returns ------- @@ -1931,7 +1931,7 @@ void CoefficientClasses(py::module &m) { array : ndarray data columns verbose : bool - display verbose information. + display verbose information (default=True) Returns ------- @@ -1956,7 +1956,7 @@ void CoefficientClasses(py::module &m) { Parameters ---------- verbose : bool - display verbose information. + display verbose information (default=True) Returns ------- @@ -1984,7 +1984,7 @@ void CoefficientClasses(py::module &m) { type : str ascii table data file verbose : bool - display verbose information. + display verbose information (default=True) Returns ------- @@ -2001,7 +2001,7 @@ void CoefficientClasses(py::module &m) { array : ndarray data columns verbose : bool - display verbose information. + display verbose information (default=True) Returns ------- From 10ad90076c3b07b6836d4e955c466139c0d7e707 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 3 Nov 2025 14:14:21 -0500 Subject: [PATCH 62/82] Fix override binding for KDdensity::getDensityAtPoint() --- expui/KDdensity.H | 5 ++++- expui/KDdensity.cc | 14 +++++++++++++- pyEXP/UtilWrappers.cc | 37 ++++++++++++++++++++++++++++++++++++- 3 files changed, 53 insertions(+), 3 deletions(-) diff --git a/expui/KDdensity.H b/expui/KDdensity.H index 46c812178..7afcd9c40 100644 --- a/expui/KDdensity.H +++ b/expui/KDdensity.H @@ -46,8 +46,11 @@ namespace Utility //! Create KD tree from a arrays of masses and positions KDdensity(const Eigen::VectorXd& m, const RowMatrixXd& d, int Ndens=32); + //! Get density at a given point + double getDensityAtPoint(double x, double y, double z); + //! Get density for particle index - double getDensity(unsigned long index); + double getDensityByIndex(unsigned long index); }; } diff --git a/expui/KDdensity.cc b/expui/KDdensity.cc index 7b995ecee..8a5252a1d 100644 --- a/expui/KDdensity.cc +++ b/expui/KDdensity.cc @@ -119,7 +119,19 @@ namespace Utility } } - double KDdensity::getDensity(unsigned long index) + double KDdensity::getDensityAtPoint(double x, double y, double z) + { + point3 query_point({x, y, z}, 0.0); + auto ret = kdtree_->nearestN(query_point, Ndens); + double volume = 4.0*M_PI/3.0*std::pow(std::get<2>(ret), 3.0); + if (volume>0.0 and KDmass>0.0) + return std::get<1>(ret)/volume/KDmass; + else + return 0.0; + } + + + double KDdensity::getDensityByIndex(unsigned long index) { return KDdens[indx[index]]; } diff --git a/pyEXP/UtilWrappers.cc b/pyEXP/UtilWrappers.cc index 682809e7c..7970fa95f 100644 --- a/pyEXP/UtilWrappers.cc +++ b/pyEXP/UtilWrappers.cc @@ -284,7 +284,42 @@ void UtilityClasses(py::module &m) { KDdensity the KDdensity instance )", py::arg("mass"), py::arg("pos"), py::arg("Ndens")=32) - .def("getDensity", &Utility::KDdensity::getDensity, + .def("getDensityAtPoint", &Utility::KDdensity::getDensityAtPoint, + R"( + Get the density estimate at a given point + + Parameters + ---------- + x : float + the x coordinate of the evaluation point + y : float + the y coordinate of the evaluation point + z : float + the z coordinate of the evaluation point + + Returns + ------- + float + the density estimate at the particle position + )", py::arg("x"), py::arg("y"), py::arg("z")) + .def("getDensityAtPoint", + [](Utility::KDdensity& A, std::vector& pos) { + return A.getDensityAtPoint(pos[0], pos[1], pos[2]); + }, + R"( + Get the density estimate at a given point + + Parameters + ---------- + pos : list(float) + x, y, z coordinates of the evaluation point + + Returns + ------- + float + the density estimate at the particle position + )", py::arg("pos")) + .def("getDensityByIndex", &Utility::KDdensity::getDensityByIndex, R"( Get the density estimate at a given particle index From fcf4fc211d7968bb9c0c3e731fb74e5d2f85628b Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 3 Nov 2025 15:14:44 -0500 Subject: [PATCH 63/82] Fix up restarts of ParticleReader to restore initial file --- exputil/ParticleReader.cc | 37 +++++++++++++++++++++++++++++++++++-- include/ParticleReader.H | 1 + 2 files changed, 36 insertions(+), 2 deletions(-) diff --git a/exputil/ParticleReader.cc b/exputil/ParticleReader.cc index 857ba367b..8ce042841 100644 --- a/exputil/ParticleReader.cc +++ b/exputil/ParticleReader.cc @@ -294,6 +294,17 @@ namespace PR { const Particle* GadgetNative::firstParticle() { + if (done) { + std::cout << "---- GadgetNative: restarting at the beginning" + << std::endl; + done = false; + curfile = _files.begin(); + if (not nextFile()) { // Try opening + std::cerr << "GadgetNative: no files found; logic error?" + << std::endl; + } + } + pcount = 0; return &particles[pcount++]; @@ -305,6 +316,7 @@ namespace PR { return &particles[pcount++]; } else { if (nextFile()) return firstParticle(); + done = true; return 0; } } @@ -650,6 +662,15 @@ namespace PR { const Particle* GadgetHDF5::firstParticle() { + if (done) { + std::cout << "---- GadgetHDF5: restarting at the beginning" << std::endl; + done = false; + curfile = _files.begin(); + if (not nextFile()) { + std::cerr << "GadgetHDF5: no files found; logic error?" << std::endl; + } + } + pcount = 0; return & particles[pcount++]; @@ -661,7 +682,8 @@ namespace PR { return & particles[pcount++]; } else { if (nextFile()) return firstParticle(); - else return 0; + done = true; + return 0; } } @@ -1221,6 +1243,16 @@ namespace PR { const Particle* PSPhdf5::firstParticle() { + if (done) { + std::cout << "---- PSPhdf5: restarting at the beginning" + << std::endl; + done = false; + curfile = _files.begin(); + if (not nextFile()) { + std::cerr << "PSPhdf5: no files found; logic error?" << std::endl; + } + } + pcount = 0; return & particles[pcount++]; @@ -1232,7 +1264,8 @@ namespace PR { return & particles[pcount++]; } else { if (nextFile()) return firstParticle(); - else return 0; + done = true; + return 0; } } diff --git a/include/ParticleReader.H b/include/ParticleReader.H index b57322dfa..b9b02fad6 100644 --- a/include/ParticleReader.H +++ b/include/ParticleReader.H @@ -35,6 +35,7 @@ namespace PR static std::vector readerTypes; int numprocs, myid; bool use_mpi; + bool done = false; //! Check for a list of directories or files and not both static bool parseDirectoryList(const std::vector& files); From f4db58dd31048e79fa91297be14454137a13bfe4 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 3 Nov 2025 16:22:28 -0500 Subject: [PATCH 64/82] Compute KDE density in real units; allow verbose flag to toggle restart message --- expui/KDdensity.cc | 6 +++--- exputil/ParticleReader.cc | 14 +++++++++----- 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/expui/KDdensity.cc b/expui/KDdensity.cc index 8a5252a1d..caf02bf4a 100644 --- a/expui/KDdensity.cc +++ b/expui/KDdensity.cc @@ -38,7 +38,7 @@ namespace Utility auto ret = kdtree_->nearestN(points[k], Ndens); double volume = 4.0*M_PI/3.0*std::pow(std::get<2>(ret), 3.0); if (volume>0.0 and KDmass>0.0) - KDdens[k] = std::get<1>(ret)/volume/KDmass; + KDdens[k] = std::get<1>(ret)/volume; else badVol++; } } @@ -100,7 +100,7 @@ namespace Utility auto ret = kdtree_->nearestN(points[k], Ndens); double volume = 4.0*M_PI/3.0*std::pow(std::get<2>(ret), 3.0); if (volume>0.0 and KDmass>0.0) - KDdens[k] = std::get<1>(ret)/volume/KDmass; + KDdens[k] = std::get<1>(ret)/volume; else badVol++; } } @@ -125,7 +125,7 @@ namespace Utility auto ret = kdtree_->nearestN(query_point, Ndens); double volume = 4.0*M_PI/3.0*std::pow(std::get<2>(ret), 3.0); if (volume>0.0 and KDmass>0.0) - return std::get<1>(ret)/volume/KDmass; + return std::get<1>(ret)/volume; else return 0.0; } diff --git a/exputil/ParticleReader.cc b/exputil/ParticleReader.cc index 8ce042841..9a177bdd2 100644 --- a/exputil/ParticleReader.cc +++ b/exputil/ParticleReader.cc @@ -295,8 +295,9 @@ namespace PR { const Particle* GadgetNative::firstParticle() { if (done) { - std::cout << "---- GadgetNative: restarting at the beginning" - << std::endl; + if (myid==0 and _verbose) + std::cout << "---- ParticleReader::GadgetNative: " + << "restarting reader" << std::endl; done = false; curfile = _files.begin(); if (not nextFile()) { // Try opening @@ -663,7 +664,9 @@ namespace PR { const Particle* GadgetHDF5::firstParticle() { if (done) { - std::cout << "---- GadgetHDF5: restarting at the beginning" << std::endl; + if (myid==0 and _verbose) + std::cout << "---- ParticleReader::GadgetHDF5: " + << "restarting reader" << std::endl; done = false; curfile = _files.begin(); if (not nextFile()) { @@ -1244,8 +1247,9 @@ namespace PR { const Particle* PSPhdf5::firstParticle() { if (done) { - std::cout << "---- PSPhdf5: restarting at the beginning" - << std::endl; + if (myid==0 and _verbose) + std::cout << "---- ParticleReader::PSPhdf5: " + << "restarting reader" << std::endl; done = false; curfile = _files.begin(); if (not nextFile()) { From 5ca913515f9ccfd081cb0b517125c88bbb704356 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 4 Nov 2025 11:09:38 -0500 Subject: [PATCH 65/82] Add sech2 toggle for gendisk --- utils/ICs/initial.cc | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/utils/ICs/initial.cc b/utils/ICs/initial.cc index 1366cfda2..e8a3b927e 100644 --- a/utils/ICs/initial.cc +++ b/utils/ICs/initial.cc @@ -264,6 +264,7 @@ double HRATIO = 1.0; double DWEIGHT = 1.0; double Mfac = 1.0; double HERNA = 0.10; +bool sech2 = true; #include "Particle.H" @@ -298,8 +299,8 @@ double DiskDens(double R, double z, double phi) { double a1 = ASCALE; double a2 = ASCALE*ARATIO; - double h1 = HSCALE; - double h2 = HSCALE*HRATIO; + double h1 = sech2 ? 0.5*HSCALE : HSCALE; + double h2 = sech2 ? 0.5*HSCALE*HRATIO : HSCALE*HRATIO; double w1 = 1.0/(1.0+DWEIGHT); double w2 = DWEIGHT/(1.0+DWEIGHT); @@ -315,7 +316,7 @@ double DiskDens(double R, double z, double phi) case DiskType::diskbulge: { double acyl = ASCALE; - double hcyl = HSCALE; + double hcyl = sech2 ? 0.5*HSCALE : HSCALE; double f = cosh(z/HSCALE); double rr = pow(pow(R, 2) + pow(z,2), 0.5); double w1 = Mfac; @@ -335,8 +336,9 @@ double DiskDens(double R, double z, double phi) case DiskType::exponential: default: { - double f = cosh(z/HSCALE); - ans = exp(-R/ASCALE)/(4.0*M_PI*ASCALE*ASCALE*HSCALE*f*f); + double h = sech2 ? 0.5*HSCALE : HSCALE; + double f = cosh(z/h); + ans = exp(-R/ASCALE)/(4.0*M_PI*ASCALE*ASCALE*h*f*f); } break; } @@ -409,6 +411,12 @@ main(int ac, char **av) options.add_options() ("h,help", "Print this help message") ("T,template", "Write template options file with current and all default values") + ("sech2", "Use sech^2 vertical profile for disk density", + cxxopts::value()->default_value("true")) + ("mtype", "EmpCylSL spherical model type (one of: Exponential, Gaussian, Plummer, Power)", + cxxopts::value(mtype)->default_value("Exponential")) + ("ppow", "Power-law index for EmpCylSL Power spherical model", + cxxopts::value(PPower)->default_value("2.0")) ("c,config", "Parameter configuration file", cxxopts::value(config)) ("deproject", "The EmpCylSL deprojection from specified disk model (EXP or MN)", From 001a701cf3e7f87acbe5d544d1ed3f859c7fd3d0 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 4 Nov 2025 17:06:07 -0500 Subject: [PATCH 66/82] Change to keep gcc 15.x happy --- include/tipsy.H | 1 + 1 file changed, 1 insertion(+) diff --git a/include/tipsy.H b/include/tipsy.H index 01be7e3df..b598fcc7e 100644 --- a/include/tipsy.H +++ b/include/tipsy.H @@ -1,6 +1,7 @@ #ifndef _tipsy_H #define _tipsy_H +#include #include // For MPI awareness From 776f7feaf08a174fff7b6008f0b6f8131f71d826 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 5 Nov 2025 09:23:48 -0500 Subject: [PATCH 67/82] Fix sampT default for Cylindrical --- expui/BiorthBasis.cc | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 2f0bfb919..06eb08875 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -26,7 +26,6 @@ namespace BasisClasses "pcavtk", "pcaeof", "subsamp", - "sampT", "snr", "samplesz", "vtkfreq", @@ -265,6 +264,8 @@ namespace BasisClasses if (conf["M0_ONLY"]) M0_only = conf["M0_ONLY"].as(); if (conf["pcavar"]) pcavar = conf["pcavar"].as(); if (conf["subsamp"]) sampT = conf["subsamp"].as(); + + sampT = std::max(1, sampT); // Sanity } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing parameter stanza for <" @@ -1446,7 +1447,7 @@ namespace BasisClasses if (unmatched.size()) throw YamlConfigError("Basis::Basis::Cylindrical", "parameter", unmatched, __FILE__, __LINE__); - int sampT = 0; + int sampT = 1; // Assign values from YAML // @@ -1501,6 +1502,9 @@ namespace BasisClasses if (conf["pcavar"] ) pcavar = conf["pcavar" ].as(); if (conf["subsamp"] ) sampT = conf["subsamp" ].as(); + // Sanity + sampT = std::max(1, sampT); + // Deprecation warning if (conf["density" ]) { if (myid==0) From f94c32ab97cfaf23897407e3fe530c65375f7993 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 5 Nov 2025 17:34:29 -0500 Subject: [PATCH 68/82] Added enableCoefCovariance() member function to override the 'pcavar' and 'subsamp' config tags. Also, double checked 'subsamp' logic to enforce nonzero-sized samples. --- expui/BiorthBasis.H | 27 ++++++++++++++++++++++++++- expui/BiorthBasis.cc | 2 -- exputil/EmpCylSL.cc | 2 +- pyEXP/BasisWrappers.cc | 20 ++++++++++++++++++++ 4 files changed, 47 insertions(+), 4 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 3b8c2c30b..00e8b7828 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -272,6 +272,7 @@ namespace BasisClasses virtual std::vector> getCoefCovariance() { + // Must be overriden; base implementation throws error throw std::runtime_error("BiorthBasis::getCoefCovariance: Not implemented for this basis"); } @@ -302,6 +303,14 @@ namespace BasisClasses //! Make covariance after accumulation virtual void makeCoefCovariance(void) {} + //! Enable covariance computation with optional sample time + virtual void enableCoefCovariance(bool pcavar, int sampT_in) + { + // Must be overriden; base implementation throws error + throw std::runtime_error("BiorthBasis::enableCoefCovariance: " + "Not implemented for this basis"); + } + //! Evaluate acceleration in Cartesian coordinates in centered //! coordinate system for a collections of points virtual RowMatrixXd& getAccel(RowMatrixXd& pos) @@ -487,6 +496,13 @@ namespace BasisClasses //! covariance matrix for subsamples of particles std::vector> getCoefCovariance(); + + virtual void enableCoefCovariance(bool pcavar, int sampT_in=100) + { + pcavar = true; + sampT = sampT_in; + if (pcavar) init_covariance(); + } }; /** @@ -919,7 +935,7 @@ namespace BasisClasses bool coefs_defined; Eigen::MatrixXd factorial; Eigen::MatrixXd expcoef; - int used; + int used, sampT = 1; std::string cachename; bool oldcache = false; @@ -1061,6 +1077,15 @@ namespace BasisClasses return {sampleCounts, sampleMasses}; } + virtual void enableCoefCovariance(bool pcavar_in, int sampT_in=100) + { + pcavar = pcavar_in; + sampT = sampT_in; + if (pcavar) { + sl->setSampT(std::max(1, sampT)); + sl->set_covar(true); + } + } }; /** diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 06eb08875..c3f35939c 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1447,8 +1447,6 @@ namespace BasisClasses if (unmatched.size()) throw YamlConfigError("Basis::Basis::Cylindrical", "parameter", unmatched, __FILE__, __LINE__); - int sampT = 1; - // Assign values from YAML // try { diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index af33cf7b4..b99ddbb3e 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -332,7 +332,7 @@ EmpCylSL::EmpCylSL(int mlim, std::string cachename) eof_made = false; sampT = 1; - defSampT = 1; + defSampT = 0; tk_type = None; cylmass = 0.0; diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 3549aeae9..a775ae6e0 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -449,6 +449,10 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE(Selement, BiorthBasis, getCovarSamples,); } + void enableCoefCovariance(bool pcavar, int nsamples) override { + PYBIND11_OVERRIDE(void, BiorthBasis, enableCoefCovariance, pcavar, nsamples); + } + }; class PySpherical : public Spherical @@ -1432,6 +1436,22 @@ void BasisFactoryClasses(py::module &m) ------- None )") + .def("enableCoefCovariance", &BasisClasses::BiorthBasis::enableCoefCovariance, + R"( + Enable or disable the coefficient covariance computation and set the + default number of paritions 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 + + Returns + ------- + None + )", py::arg("pcavar"), py::arg("nsamples")=100) .def("setCovarH5Compress", &BasisClasses::BiorthBasis::setCovarH5Compress, R"( Set the HDF5 compression level for covariance storage in HDF5. The Szip From 240d69c8f4748d96e3b0b65c169617f3cfa9b940 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 6 Nov 2025 12:09:29 -0500 Subject: [PATCH 69/82] Fix typo in enableCoefCovariance for Spherical --- expui/BiorthBasis.H | 4 ++-- pyEXP/BasisWrappers.cc | 8 ++++++++ 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 00e8b7828..2ca9095e2 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -497,9 +497,9 @@ namespace BasisClasses std::vector> getCoefCovariance(); - virtual void enableCoefCovariance(bool pcavar, int sampT_in=100) + virtual void enableCoefCovariance(bool pcavar_in, int sampT_in=100) { - pcavar = true; + pcavar = pcavar_in; sampT = sampT_in; if (pcavar) init_covariance(); } diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index a775ae6e0..662bd91f8 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -518,6 +518,10 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE_PURE(std::vector, Spherical, orthoCheck, knots); } + void enableCoefCovariance(bool pcavar, int nsamples) override { + PYBIND11_OVERRIDE(void, Spherical, enableCoefCovariance, pcavar, nsamples); + } + }; @@ -569,6 +573,10 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE(void, Cylindrical, make_coefs,); } + void enableCoefCovariance(bool pcavar, int nsamples) override { + PYBIND11_OVERRIDE(void, Cylindrical, enableCoefCovariance, pcavar, nsamples); + } + }; From 227be01e6718a730c5be86f9e24b862c971cda37 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 6 Nov 2025 13:55:52 -0500 Subject: [PATCH 70/82] Restore the original default of TRUE for ENABLE_UTILS to allow the EXP unit tests to run more fully --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index a14956616..5d7e5cb8c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -49,7 +49,7 @@ endif() option(ENABLE_NBODY "Enable EXP n-body" ON) option(ENABLE_PYEXP "Enable the Python bindings" ON) option(ENABLE_PYEXP_ONLY "Python bindings and support libraries only" OFF) -option(ENABLE_UTILS "Enable build of the EXP standalone utilities" OFF) +option(ENABLE_UTILS "Enable build of the EXP standalone utilities" ON) option(ENABLE_PNG "Enable PNG graphics support" FALSE) option(ENABLE_CUDA "Enable CUDA" FALSE) option(ENABLE_SLURM "Enable SLURM checkpointing support" FALSE) From 8cb1426168b05a8cfdc88d8ffa91eeeae56a91fa Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 10 Nov 2025 11:21:19 -0500 Subject: [PATCH 71/82] Remove some cruft --- expui/BiorthBasis.cc | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index c3f35939c..21ff0807c 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -609,10 +609,9 @@ namespace BasisClasses for (int l=0, loffset=0, L=0; l<=lmax; loffset+=(2*l+1), l++) { Eigen::VectorXd workE; - int esize = (l+1)*nmax; // M loop - for (int m=0, moffset=0, moffE=0; m<=l; m++, L++) { + for (int m=0, moffset=0; m<=l; m++, L++) { fac = factorial(l, m) * legs[tid](l, m); From 720d9aec942e5f1c280f35a6619b1c77940a69ff Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 10 Nov 2025 23:05:26 -0500 Subject: [PATCH 72/82] Fixed halo covariance MPI type mistake --- expui/BiorthBasis.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 21ff0807c..95474cd34 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -611,7 +611,7 @@ namespace BasisClasses Eigen::VectorXd workE; // M loop - for (int m=0, moffset=0; m<=l; m++, L++) { + for (int m=0, moffset=0, moffE=0; m<=l; m++, L++) { fac = factorial(l, m) * legs[tid](l, m); @@ -684,10 +684,10 @@ namespace BasisClasses for (int l=0; l Date: Thu, 13 Nov 2025 16:16:32 -0500 Subject: [PATCH 73/82] Added covariance computation for the Cube --- expui/BiorthBasis.H | 31 +++++++ expui/BiorthBasis.cc | 183 +++++++++++++++++++++++++++++++++++++++-- pyEXP/BasisWrappers.cc | 18 ++++ utils/ICs/initial.cc | 34 ++++---- 4 files changed, 240 insertions(+), 26 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 2ca9095e2..ef877f49a 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -1278,6 +1278,7 @@ namespace BasisClasses unsigned used; using coefType = Eigen::Tensor, 3>; + using covrType = Eigen::Tensor, 4>; coefType expcoef; @@ -1322,6 +1323,18 @@ namespace BasisClasses void computeAccel(double x, double y, double z, Eigen::Ref vstore); + //@{ + //! Covariance structures. First index is T, second is the + //! flattened 3-d k vector + std::tuple index3D(unsigned k); + unsigned index1D(int kx, int ky, int kz); + std::vector meanV; + std::vector covrV; + void init_covariance(); + void zero_covariance(); + void writeCovarH5Params(HighFive::File& file); + int sampT = 100; + //@} public: @@ -1371,6 +1384,24 @@ namespace BasisClasses return true; } + //! Return a vector of tuples of basis functions and the + //! covariance matrix for subsamples of particles + std::vector> getCoefCovariance(); + + //! Get sample count and mass arrays + std::tuple + getCovarSamples() + { + return {sampleCounts, sampleMasses}; + } + + //! Enable coefficient covariance computation + virtual void enableCoefCovariance(bool pcavar_in, int sampT_in=100) + { + pcavar = pcavar_in; + sampT = sampT_in; + } + }; class CovarianceReader diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 95474cd34..8422f7dd0 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3819,17 +3819,28 @@ namespace BasisClasses "knots", "verbose", "check", - "method" + "method", + "pcavar," + "subsamp" }; Cube::Cube(const YAML::Node& CONF) : BiorthBasis(CONF, "cube") { initialize(); + + // Initialize covariance + // + if (pcavar) init_covariance(); } Cube::Cube(const std::string& confstr) : BiorthBasis(confstr, "cube") { initialize(); + + + // Initialize covariance + // + if (pcavar) init_covariance(); } void Cube::initialize() @@ -3866,17 +3877,20 @@ namespace BasisClasses // Assign values from YAML // try { - if (conf["nminx"]) nminx = conf["nminx"].as(); - if (conf["nminy"]) nminy = conf["nminy"].as(); - if (conf["nminz"]) nminz = conf["nminz"].as(); + if (conf["nminx"]) nminx = conf["nminx"].as(); + if (conf["nminy"]) nminy = conf["nminy"].as(); + if (conf["nminz"]) nminz = conf["nminz"].as(); - if (conf["nmaxx"]) nmaxx = conf["nmaxx"].as(); - if (conf["nmaxy"]) nmaxy = conf["nmaxy"].as(); - if (conf["nmaxz"]) nmaxz = conf["nmaxz"].as(); + if (conf["nmaxx"]) nmaxx = conf["nmaxx"].as(); + if (conf["nmaxy"]) nmaxy = conf["nmaxy"].as(); + if (conf["nmaxz"]) nmaxz = conf["nmaxz"].as(); - if (conf["knots"]) knots = conf["knots"].as(); + if (conf["knots"]) knots = conf["knots"].as(); - if (conf["check"]) check = conf["check"].as(); + if (conf["check"]) check = conf["check"].as(); + + if (conf["pcavar"]) pcavar = conf["pcavar"].as(); + if (conf["subsamp"]) sampT = conf["subsamp"].as(); } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing parameter stanza for <" @@ -3916,6 +3930,7 @@ namespace BasisClasses expcoef.setZero(); totalMass = 0.0; used = 0; + if (pcavar) zero_covariance(); } @@ -3998,6 +4013,9 @@ namespace BasisClasses std::exp(-kfac*(y*nmaxy)), std::exp(-kfac*(z*nmaxz))}; + unsigned Itot = (2*nmaxx+1)*(2*nmaxy+1)*(2*nmaxz+1); + Eigen::VectorXcd g(Itot); + Eigen::Vector3cd curr(init); for (int ix=0; ix<=2*nmaxx; ix++, curr(0)*=step(0)) { curr(1) = init(1); @@ -4016,9 +4034,25 @@ namespace BasisClasses double norm = 1.0/sqrt(M_PI*(ii*ii + jj*jj + kk*kk));; expcoef(ix, iy, iz) += - mass * curr(0)*curr(1)*curr(2) * norm; + + if (pcavar) + g[index1D(ix, iy, iz)] = curr(0)*curr(1)*curr(2) * norm; } } } + + if (pcavar) { + // Sample index for pcavar + int T = 0; + T = used % sampT; + sampleCounts(T) += 1; + sampleMasses(T) += mass; + + meanV[T].noalias() += g * mass; + covrV[T].noalias() += g * g.adjoint() * mass; + } + + used++; } void Cube::make_coefs() @@ -4030,9 +4064,50 @@ namespace BasisClasses MPI_Allreduce(MPI_IN_PLACE, expcoef.data(), expcoef.size(), MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD); + + + if (pcavar) { + + MPI_Allreduce(MPI_IN_PLACE, sampleCounts.data(), sampleCounts.size(), + MPI_INT, MPI_SUM, MPI_COMM_WORLD); + + MPI_Allreduce(MPI_IN_PLACE, sampleMasses.data(), sampleMasses.size(), + MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + + for (int T=0; T> + Cube::getCoefCovariance() + { + std::vector> ret; + if (pcavar) { + ret.resize(sampT); + for (int T=0; T(ret[T][0]) = meanV[T]; + std::get<1>(ret[T][0]) = covrV[T]; + } + } + + return ret; + } + + std::vector Cube::crt_eval(double x, double y, double z) { // Get thread id @@ -4139,6 +4214,86 @@ namespace BasisClasses return ret; } + void Cube::init_covariance() + { + if (pcavar) { + int Itot = (2*nmaxx+1)*(2*nmaxy+1)*(2*nmaxz+1); + + meanV.resize(sampT); + for (auto& v : meanV) { + v.resize(Itot); + } + + covrV.resize(sampT); + for (auto& v : covrV) { + v.resize(Itot, Itot); + } + + sampleCounts.resize(sampT); + sampleMasses.resize(sampT); + + zero_covariance(); + } + } + + + void Cube::zero_covariance() + { + for (int T=0; T nmaxx) { + std::ostringstream sout; + sout << "Cube::index1d: x index [" << kx << "] must be in [" << -nmaxx << ", " << nmaxx << "]"; + throw std::runtime_error(sout.str()); + } + + if (ky <-nmaxy or ky > nmaxy) { + std::ostringstream sout; + sout << "Cube::index1d: y index [" << ky << "] must be in [" << -nmaxy << ", " << nmaxy << "]"; + throw std::runtime_error(sout.str()); + } + + if (kz <-nmaxz or kx > nmaxz) { + std::ostringstream sout; + sout << "Cube::index1d: z index [" << kz << "] must be in [" << -nmaxz << ", " << nmaxz << "]"; + throw std::runtime_error(sout.str()); + } + + return (kx + nmaxx)*(2*nmaxy+1)*(2*nmaxz+1) + (ky + nmaxy)*(2*nmaxz+1) + kz + nmaxz; + } + + std::tuple Cube::index3D(unsigned indx) + { + int Itot = (2*nmaxx+1)*(2*nmaxy+1)*(2*nmaxz+1); + + // Sanity check + // + if (indx >= Itot) { + std::ostringstream sout; + sout << "Cube::index3d: index [" << indx << "] must be in 0 <= indx < " << Itot; + throw std::runtime_error(sout.str()); + } + + // Compute the 3d index + // + int ix = indx/((2*nmaxy+1)*(2*nmaxz+1)); + int iy = (indx - ix*(2*nmaxy+1)*(2*nmaxz+1))/(2*nmaxz+1); + int iz = indx - ix*(2*nmaxy+1)*(2*nmaxz+1)/(2*nmaxz+1) - iy*(2*nmaxz+1); + + return {ix - nmaxx, iy - nmaxy, iz - nmaxz}; + } + + // Generate coeffients from a particle reader CoefClasses::CoefStrPtr BiorthBasis::createFromReader (PR::PRptr reader, Eigen::Vector3d ctr, RowMatrix3d rot) @@ -4834,6 +4989,16 @@ namespace BasisClasses file.createAttribute("hcyl", HighFive::DataSpace::From(hcyl)).write(hcyl); } + void Cube::writeCovarH5Params(HighFive::File& file) + { + file.createAttribute("nminx", HighFive::DataSpace::From(nminx)).write(nminx); + file.createAttribute("nminy", HighFive::DataSpace::From(nminy)).write(nminy); + file.createAttribute("nminz", HighFive::DataSpace::From(nminz)).write(nminz); + file.createAttribute("nmaxx", HighFive::DataSpace::From(nmaxx)).write(nmaxx); + file.createAttribute("nmaxy", HighFive::DataSpace::From(nmaxy)).write(nmaxy); + file.createAttribute("nmaxz", HighFive::DataSpace::From(nmaxz)).write(nmaxz); + } + unsigned BiorthBasis::writeCovarH5(HighFive::Group& snaps, unsigned count, double time) { std::ostringstream stim; diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 662bd91f8..842ca2002 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -868,6 +868,24 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE(void, Cube, make_coefs,); } + using CCelement = std::tuple; + using CCreturn = std::vector>; + + CCreturn getCoefCovariance(void) override { + PYBIND11_OVERRIDE(CCreturn, Cube, getCoefCovariance,); + } + + using Selement = std::tuple; + + Selement getCovarSamples() override { + PYBIND11_OVERRIDE(Selement, Cube, getCovarSamples,); + } + + void enableCoefCovariance(bool pcavar, int nsamples) override { + PYBIND11_OVERRIDE(void, Cube, enableCoefCovariance, pcavar, nsamples); + } + + }; diff --git a/utils/ICs/initial.cc b/utils/ICs/initial.cc index e8a3b927e..fc4f43c40 100644 --- a/utils/ICs/initial.cc +++ b/utils/ICs/initial.cc @@ -401,7 +401,7 @@ main(int ac, char **av) bool const_height, images, multi, basis, zeropos, zerovel; bool report, ignore, evolved, diskmodel; int nhalo, ndisk, ngas, ngparam; - std::string hbods, dbods, gbods, suffix, centerfile, halofile1, halofile2; + std::string hbods, dbods, gbods, outtag, runtag, centerfile, halofile1, halofile2; std::string cachefile, config, gentype, dtype, dmodel, mtype, ctype; const std::string mesg("Generates a Monte Carlo realization of a halo with an\n embedded disk using Jeans' equations\n"); @@ -607,10 +607,10 @@ main(int ac, char **av) cxxopts::value(Tmin)->default_value("500.0")) ("centerfile", "File containing phase-space center", cxxopts::value(centerfile)) - ("runtag", "Prefix for output files", - cxxopts::value(runtag)->default_value("gendisk")) - ("suffix", "Suffix for output files (none by default)", - cxxopts::value(suffix)) + ("suffix", "Prefix for output files", + cxxopts::value(outtag)->default_value("gendisk")) + ("runtag", "Suffix for output files (none by default)", + cxxopts::value(runtag)) ("threads", "Number of threads to run", cxxopts::value(nthrds)->default_value("1")) ("allow", "Allow multimass algorithm to generature negative masses for testing") @@ -807,10 +807,10 @@ main(int ac, char **av) int n_particlesH, n_particlesD, n_particlesG; - if (suffix.size()>0) { - hbods = hbods + "." + suffix; - dbods = dbods + "." + suffix; - gbods = gbods + "." + suffix; + if (runtag.size()>0) { + hbods = hbods + "." + runtag; + dbods = dbods + "." + runtag; + gbods = gbods + "." + runtag; } // Divvy up the particles by core @@ -874,7 +874,7 @@ main(int ac, char **av) if (vm.count("itmax")) DiskHalo::ITMAX = itmax; if (vm.count("allow")) DiskHalo::ALLOW = true; if (vm.count("nomono")) DiskHalo::use_mono = false; - if (suffix.size()) DiskHalo::RUNTAG = suffix; + if (runtag.size()) DiskHalo::RUNTAG = runtag; AddDisk::use_mpi = true; AddDisk::Rmin = RMIN; @@ -986,7 +986,7 @@ main(int ac, char **av) // Basis orthgonality check // if (vm.count("ortho")) { - std::ofstream out(runtag + ".ortho_check"); + std::ofstream out(outtag + ".ortho_check"); expandd->ortho_check(out); } } @@ -1187,7 +1187,7 @@ main(int ac, char **av) std::cout << "Dumping coefficients halo . . . " << std::flush; ostringstream sout; sout << "halo_coefs."; - if (suffix.size()>0) sout << suffix; + if (outtag.size()>0) sout << outtag; else sout << "dump"; ofstream out(sout.str()); if (out) expandh->dump_coefs(out, false); @@ -1197,7 +1197,7 @@ main(int ac, char **av) std::cout << "Dumping a probe through the halo . . . " << std::flush; ostringstream sout; sout << "halo_probe."; - if (suffix.size()>0) sout << suffix; + if (outtag.size()>0) sout << outtag; else sout << "dump"; ofstream out(sout.str()); if (out) { @@ -1287,7 +1287,7 @@ main(int ac, char **av) std::cout << "Dumping coefficients . . . " << std::flush; ostringstream sout; sout << "disk_coefs."; - if (suffix.size()>0) sout << suffix; + if (outtag.size()>0) sout << outtag; else sout << "dump"; ofstream out(sout.str()); if (out) expandd->dump_coefs(out); @@ -1333,11 +1333,11 @@ main(int ac, char **av) if (ndisk) { int nout = 200; - string dumpstr = runtag + ".dump"; + string dumpstr = outtag + ".dump"; expandd->dump_basis(dumpstr, 0); - expandd->dump_images(runtag, 5.0*scale_length, 5.0*scale_height, + expandd->dump_images(outtag, 5.0*scale_length, 5.0*scale_height, nout, nout, false); - expandd->dump_images_basis(runtag, 5.0*scale_length, 5.0*scale_height, + expandd->dump_images_basis(outtag, 5.0*scale_length, 5.0*scale_height, nout, nout, false, 0, MMAX, 0, NMAXD-1); } From e823e4d6c917d361fb7c6f525da9439144b1d6dd Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 13 Nov 2025 16:54:54 -0500 Subject: [PATCH 74/82] Missing stanza for deducing Cube coefficients --- expui/BiorthBasis.cc | 1 - expui/Coefficients.cc | 2 ++ 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 8422f7dd0..2a53be3f6 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3837,7 +3837,6 @@ namespace BasisClasses { initialize(); - // Initialize covariance // if (pcavar) init_covariance(); diff --git a/expui/Coefficients.cc b/expui/Coefficients.cc index 6f91a1129..f79ad70dd 100644 --- a/expui/Coefficients.cc +++ b/expui/Coefficients.cc @@ -2543,6 +2543,8 @@ namespace CoefClasses ret = std::make_shared(); } else if (dynamic_cast(coef.get())) { ret = std::make_shared(); + } else if (dynamic_cast(coef.get())) { + ret = std::make_shared(); } else if (dynamic_cast(coef.get())) { ret = std::make_shared(); } else if (dynamic_cast(coef.get())) { From e8da520aa69d3faa3a107dba579ca61964e68a3b Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 13 Nov 2025 17:37:27 -0500 Subject: [PATCH 75/82] More Cube updates for BiorthBasis --- expui/BiorthBasis.H | 9 ++++++++- expui/BiorthBasis.cc | 43 +++++++++++++++++++++++++++++-------------- 2 files changed, 37 insertions(+), 15 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index ef877f49a..46a2dfe15 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -1271,6 +1271,9 @@ namespace BasisClasses //! Maximum expansion order for construction int nmaxx, nmaxy, nmaxz; + //! Total number of modes + unsigned Itot; + //! Number of integration knots for orthogonal check int knots; @@ -1282,7 +1285,7 @@ namespace BasisClasses coefType expcoef; - //! Notal mass on grid + //! Total mass on grid double totalMass; //! Number of particles @@ -1400,6 +1403,10 @@ namespace BasisClasses { pcavar = pcavar_in; sampT = sampT_in; + if (pcavar) { + init_covariance(); + zero_covariance(); + } } }; diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 2a53be3f6..626d63958 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -590,6 +590,7 @@ namespace BasisClasses if (r < rmin or r > rmax) return; + // Update counters used++; totalMass += mass; @@ -2206,6 +2207,7 @@ namespace BasisClasses if (R < ortho->getRtable() and fabs(z) < ortho->getRtable()) { + // Update counters used++; totalMass += mass; @@ -3009,6 +3011,7 @@ namespace BasisClasses // Get thread id int tid = omp_get_thread_num(); + // Update counters used++; totalMass += mass; @@ -3464,6 +3467,7 @@ namespace BasisClasses else y -= std::floor( y); + // Update counters used++; // Storage for basis evaluation @@ -3914,10 +3918,17 @@ namespace BasisClasses // int nthrds = omp_get_max_threads(); + // Total number of wavenumbers + // + Itot = (2*nmaxx + 1) * (2*nmaxy + 1) * (2*nmaxz + 1); + expcoef.resize(2*nmaxx+1, 2*nmaxy+1, 2*nmaxz+1); expcoef.setZero(); + // Counters + // used = 0; + totalMass = 0.0; // Set cartesian coordindates // @@ -4001,6 +4012,9 @@ namespace BasisClasses else z -= std::floor( z); + // Update counters + used++; + totalMass += mass; // Recursion multipliers Eigen::Vector3cd step @@ -4012,8 +4026,8 @@ namespace BasisClasses std::exp(-kfac*(y*nmaxy)), std::exp(-kfac*(z*nmaxz))}; - unsigned Itot = (2*nmaxx+1)*(2*nmaxy+1)*(2*nmaxz+1); - Eigen::VectorXcd g(Itot); + Eigen::VectorXcd g; + if (pcavar) g.resize(Itot); Eigen::Vector3cd curr(init); for (int ix=0; ix<=2*nmaxx; ix++, curr(0)*=step(0)) { @@ -4035,7 +4049,7 @@ namespace BasisClasses expcoef(ix, iy, iz) += - mass * curr(0)*curr(1)*curr(2) * norm; if (pcavar) - g[index1D(ix, iy, iz)] = curr(0)*curr(1)*curr(2) * norm; + g[index1D(ii, jj, kk)] = curr(0)*curr(1)*curr(2) * norm; } } } @@ -4050,8 +4064,6 @@ namespace BasisClasses meanV[T].noalias() += g * mass; covrV[T].noalias() += g * g.adjoint() * mass; } - - used++; } void Cube::make_coefs() @@ -4061,6 +4073,9 @@ namespace BasisClasses MPI_Allreduce(MPI_IN_PLACE, &used, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &totalMass, 1, MPI_DOUBLE, + MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, expcoef.data(), expcoef.size(), MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD); @@ -4216,8 +4231,7 @@ namespace BasisClasses void Cube::init_covariance() { if (pcavar) { - int Itot = (2*nmaxx+1)*(2*nmaxy+1)*(2*nmaxz+1); - + meanV.resize(sampT); for (auto& v : meanV) { v.resize(Itot); @@ -4252,34 +4266,35 @@ namespace BasisClasses { if (kx <-nmaxx or kx > nmaxx) { std::ostringstream sout; - sout << "Cube::index1d: x index [" << kx << "] must be in [" << -nmaxx << ", " << nmaxx << "]"; + sout << "Cube::index1D: x index [" << kx << "] must be in [" << -nmaxx << ", " << nmaxx << "]"; throw std::runtime_error(sout.str()); } if (ky <-nmaxy or ky > nmaxy) { std::ostringstream sout; - sout << "Cube::index1d: y index [" << ky << "] must be in [" << -nmaxy << ", " << nmaxy << "]"; + sout << "Cube::index1D: y index [" << ky << "] must be in [" << -nmaxy << ", " << nmaxy << "]"; throw std::runtime_error(sout.str()); } if (kz <-nmaxz or kx > nmaxz) { std::ostringstream sout; - sout << "Cube::index1d: z index [" << kz << "] must be in [" << -nmaxz << ", " << nmaxz << "]"; + sout << "Cube::index1D: z index [" << kz << "] must be in [" << -nmaxz << ", " << nmaxz << "]"; throw std::runtime_error(sout.str()); } - return (kx + nmaxx)*(2*nmaxy+1)*(2*nmaxz+1) + (ky + nmaxy)*(2*nmaxz+1) + kz + nmaxz; + return + (kx + nmaxx)*(2*nmaxy+1)*(2*nmaxz+1) + + (ky + nmaxy)*(2*nmaxz+1) + + (kz + nmaxz); } std::tuple Cube::index3D(unsigned indx) { - int Itot = (2*nmaxx+1)*(2*nmaxy+1)*(2*nmaxz+1); - // Sanity check // if (indx >= Itot) { std::ostringstream sout; - sout << "Cube::index3d: index [" << indx << "] must be in 0 <= indx < " << Itot; + sout << "Cube::index3D: index [" << indx << "] must be in 0 <= indx < " << Itot; throw std::runtime_error(sout.str()); } From a7cda15c51f58c84e53deef465ed5d914764ca92 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 13 Nov 2025 17:44:18 -0500 Subject: [PATCH 76/82] Comment in typo only [no CI] --- expui/BiorthBasis.H | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 46a2dfe15..a4b8b83f9 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -1139,7 +1139,7 @@ namespace BasisClasses coefType expcoef; - //! Notal mass on grid + //! Total mass on grid double totalMass; //! Number of particles @@ -1281,7 +1281,6 @@ namespace BasisClasses unsigned used; using coefType = Eigen::Tensor, 3>; - using covrType = Eigen::Tensor, 4>; coefType expcoef; From b0ac36c77665af23d42829580ad34dee78781d95 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 14 Nov 2025 09:19:40 -0500 Subject: [PATCH 77/82] Added pyEXP bindings for index helpers for Cube covariance --- expui/BiorthBasis.H | 7 +++++-- expui/BiorthBasis.cc | 21 ++++++++++++--------- pyEXP/BasisWrappers.cc | 40 +++++++++++++++++++++++++++++++++++++++- 3 files changed, 56 insertions(+), 12 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index a4b8b83f9..91c251d47 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -1328,8 +1328,6 @@ namespace BasisClasses //@{ //! Covariance structures. First index is T, second is the //! flattened 3-d k vector - std::tuple index3D(unsigned k); - unsigned index1D(int kx, int ky, int kz); std::vector meanV; std::vector covrV; void init_covariance(); @@ -1408,6 +1406,11 @@ namespace BasisClasses } } + //! Return wave-number indices from flattened index + std::tuple index3D(unsigned k); + + //! Get flattened index from wave-number indicies + unsigned index1D(int kx, int ky, int kz); }; class CovarianceReader diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 626d63958..17e7611d9 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -4049,7 +4049,7 @@ namespace BasisClasses expcoef(ix, iy, iz) += - mass * curr(0)*curr(1)*curr(2) * norm; if (pcavar) - g[index1D(ii, jj, kk)] = curr(0)*curr(1)*curr(2) * norm; + g[index1D(ix, iy, iz)] = curr(0)*curr(1)*curr(2) * norm; } } } @@ -4264,28 +4264,31 @@ namespace BasisClasses unsigned Cube::index1D(int kx, int ky, int kz) { - if (kx <-nmaxx or kx > nmaxx) { + if (kx <0 or kx > 2*nmaxx) { std::ostringstream sout; - sout << "Cube::index1D: x index [" << kx << "] must be in [" << -nmaxx << ", " << nmaxx << "]"; + sout << "Cube::index1D: x index [" << kx << "] must be in [0, " + << 2*nmaxx << "]"; throw std::runtime_error(sout.str()); } if (ky <-nmaxy or ky > nmaxy) { std::ostringstream sout; - sout << "Cube::index1D: y index [" << ky << "] must be in [" << -nmaxy << ", " << nmaxy << "]"; + sout << "Cube::index1D: y index [" << ky << "] must be in [0, " + << 2*nmaxy << "]"; throw std::runtime_error(sout.str()); } if (kz <-nmaxz or kx > nmaxz) { std::ostringstream sout; - sout << "Cube::index1D: z index [" << kz << "] must be in [" << -nmaxz << ", " << nmaxz << "]"; + sout << "Cube::index1D: z index [" << kz << "] must be in [0, " + << 2*nmaxz << "]"; throw std::runtime_error(sout.str()); } return - (kx + nmaxx)*(2*nmaxy+1)*(2*nmaxz+1) + - (ky + nmaxy)*(2*nmaxz+1) + - (kz + nmaxz); + kx*(2*nmaxy+1)*(2*nmaxz+1) + + ky*(2*nmaxz+1) + + kz; } std::tuple Cube::index3D(unsigned indx) @@ -4304,7 +4307,7 @@ namespace BasisClasses int iy = (indx - ix*(2*nmaxy+1)*(2*nmaxz+1))/(2*nmaxz+1); int iz = indx - ix*(2*nmaxy+1)*(2*nmaxz+1)/(2*nmaxz+1) - iy*(2*nmaxz+1); - return {ix - nmaxx, iy - nmaxy, iz - nmaxz}; + return {ix, iy, iz} } diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 842ca2002..d5e813665 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -2281,7 +2281,7 @@ void BasisFactoryClasses(py::module &m) Returns ------- - Cube + Slab the new instance )", py::arg("YAMLstring")) .def("getBasis", &BasisClasses::Slab::getBasis, @@ -2366,6 +2366,44 @@ void BasisFactoryClasses(py::module &m) Cube the new instance )", py::arg("YAMLstring")) + .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. + + Parameters + ---------- + kx : int + wave-number index in [0, 2*nmax_x] + kx : int + wave-number index in [0, 2*nmax_y] + kx : int + wave-number index in [0, 2*nmax_z] + + Returns + ------- + indx : int + flattened index for coefficient array and covariance matrix + )", + py::arg("kx"), py::arg("ky"), py::arg("kz")) + .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. + + Parameters + ---------- + index : int + flattened index for coefficient array and covariance matrix + + Returns + ------- + indices : tuple(int, int, int) + wave-number indices for each dimension in the ranges [0, 2*nmax-x], + [0, 2*nmax_y], [0, 2*nmax_z] + )", + py::arg("index")); + .def("orthoCheck", [](BasisClasses::Cube& A) { return A.orthoCheck(); From 5b2ec358001c0388e064cb18030bfd36f761466e Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 14 Nov 2025 11:04:28 -0500 Subject: [PATCH 78/82] Expose index manipulators from Cube basis in pyEXP --- expui/BiorthBasis.cc | 22 +++++++++++++++++----- pyEXP/BasisWrappers.cc | 24 ++++++++++++++---------- 2 files changed, 31 insertions(+), 15 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 17e7611d9..e22de62e6 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3852,9 +3852,9 @@ namespace BasisClasses // BasisID = "Cube"; - nminx = std::numeric_limits::max(); - nminy = std::numeric_limits::max(); - nminz = std::numeric_limits::max(); + nminx = 0; + nminy = 0; + nminz = 0; nmaxx = 6; nmaxy = 6; @@ -4027,7 +4027,10 @@ namespace BasisClasses std::exp(-kfac*(z*nmaxz))}; Eigen::VectorXcd g; - if (pcavar) g.resize(Itot); + if (pcavar) { + g.resize(Itot); + g.setZero(); + } Eigen::Vector3cd curr(init); for (int ix=0; ix<=2*nmaxx; ix++, curr(0)*=step(0)) { @@ -4043,6 +4046,9 @@ namespace BasisClasses int jj = iy-nmaxy; int kk = iz-nmaxz; + // Limit to minimum wave number + if (abs(ii) Date: Fri, 14 Nov 2025 11:27:09 -0500 Subject: [PATCH 79/82] Fix some minor doc-string typos --- pyEXP/BasisWrappers.cc | 58 +++++++++++++++++++++--------------------- 1 file changed, 29 insertions(+), 29 deletions(-) diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index d41002533..be739b112 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -51,12 +51,12 @@ void BasisFactoryClasses(py::module &m) 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 thosed used by EXP. + 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 biorthgonal bases in EXP but not in pyEXP + 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 @@ -102,7 +102,7 @@ void BasisFactoryClasses(py::module &m) def myFunctor(m, pos, vel, index): ret = False # Default return value - # some caculation with scalar mass, pos array, vel array and + # some calculation with scalar mass, pos array, vel array and # integer index that sets ret to True if desired . . . return ret @@ -128,7 +128,7 @@ void BasisFactoryClasses(py::module &m) phase-space field functor to cylindrical or spherical velocities based on the 'dof' parameter. More on 'dof' below. - Scalablility + Scalability ------------ createFromArray() is a convenience method allows you to transform coordinates and preprocess phase space using your own methods and @@ -195,7 +195,7 @@ void BasisFactoryClasses(py::module &m) 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 implementes + 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. @@ -1132,7 +1132,7 @@ void BasisFactoryClasses(py::module &m) See also -------- setNonInertial : set non-inertial data from an Orient file - setNonInertialAccel : set the non-inertial acceration + setNonInertialAccel : set the non-inertial acceleration )", py::arg("N"), py::arg("times"), py::arg("pos") ) @@ -1159,7 +1159,7 @@ void BasisFactoryClasses(py::module &m) See also -------- setNonInertial : set non-inertial data from a time series of values - setNonInertialAccel : set the non-inertial acceration + setNonInertialAccel : set the non-inertial acceleration )", py::arg("N"), py::arg("orient") ) @@ -1368,7 +1368,7 @@ void BasisFactoryClasses(py::module &m) The first dimension are the time samples. The second dimension is the angular index. Each element is a tuple of the coefficient vector and covariance. - The values are complex128 and represet the full amplitude and phase information. + The values are complex128 and represents the full amplitude and phase information. Returns ------- @@ -1465,7 +1465,7 @@ void BasisFactoryClasses(py::module &m) .def("enableCoefCovariance", &BasisClasses::BiorthBasis::enableCoefCovariance, R"( Enable or disable the coefficient covariance computation and set the - default number of paritions to use for the covariance computation. + default number of partitions to use for the covariance computation. Parameters ---------- @@ -1502,13 +1502,13 @@ void BasisFactoryClasses(py::module &m) )", py::arg("compress")=5, py::arg("chunkSize")=1024*1024, py::arg("shuffle")=true, py::arg("azip")=false) .def("getFields", &BasisClasses::BiorthBasis::getFields, R"( - Return the field evaluations for a given cartesian position. The + Return the field evaluations for a given Cartesian position. The fields include density, potential, and force. The density and potential evaluations are separated into full, axisymmetric and non-axisymmetric contributions. You can get the fields labels by using the __call__ method of the - basis object. This is equilevalent to a tuple of the getFields() + basis object. This is equivalent to a tuple of the getFields() output with a list of field labels. Parameters @@ -1532,7 +1532,7 @@ void BasisFactoryClasses(py::module &m) py::arg("x"), py::arg("y"), py::arg("z")) .def("getAccel", py::overload_cast(&BasisClasses::BiorthBasis::getAccel), R"( - Return the acceleration for a given cartesian position + Return the acceleration for a given Cartesian position Parameters ---------- @@ -1556,7 +1556,7 @@ void BasisFactoryClasses(py::module &m) py::arg("x"), py::arg("y"), py::arg("z")) .def("getAccel", py::overload_cast(&BasisClasses::BiorthBasis::getAccel), R"( - Return the acceleration for a given cartesian position + Return the acceleration for a given Cartesian position Parameters ---------- @@ -1580,7 +1580,7 @@ void BasisFactoryClasses(py::module &m) py::arg("x"), py::arg("y"), py::arg("z")) .def("getAccel", py::overload_cast(&BasisClasses::BiorthBasis::getAccel), R"( - Return the acceleration for an array of cartesian positions + Return the acceleration for an array of Cartesian positions Parameters ---------- @@ -1600,7 +1600,7 @@ void BasisFactoryClasses(py::module &m) py::arg("pos")) .def("getAccelArray", py::overload_cast(&BasisClasses::BiorthBasis::getAccel), R"( - Return the acceleration for a given cartesian position + Return the acceleration for a given Cartesian position Parameters ---------- @@ -1624,7 +1624,7 @@ void BasisFactoryClasses(py::module &m) py::arg("x"), py::arg("y"), py::arg("z")) .def("getFieldsCoefs", &BasisClasses::BiorthBasis::getFieldsCoefs, R"( - Return the field evaluations for a given cartesian position + Return the field evaluations for a given Cartesian position for every frame in a coefficient set. The field evaluations are produced by a call to getFields(). @@ -1661,7 +1661,7 @@ void BasisFactoryClasses(py::module &m) coordinates for SphericalSL and Bessel, cylindrical coordinates for Cylindrical, FlatDisk, and CBDisk, and Cartesian coordinates for the Slab and Cube. This member function can be used to override the - default. The available coorindates are: 'spherical', 'cylindrical', + default. The available coordinates are: 'spherical', 'cylindrical', 'cartesian'. Parameters @@ -1730,7 +1730,7 @@ void BasisFactoryClasses(py::module &m) )") .def("make_coefs", &BasisClasses::BiorthBasis::make_coefs, R"( - Create the coefficients after particle accumuluation is complete + Create the coefficients after particle accumulation is complete Returns ------- @@ -1812,7 +1812,7 @@ void BasisFactoryClasses(py::module &m) return A.orthoCheck(knots); }, R"( - Check orthgonality of basis functions by quadrature + Check orthogonality of basis functions by quadrature Inner-product matrix of Sturm-Liouville solutions indexed by harmonic order used to assess fidelity. @@ -1899,7 +1899,7 @@ void BasisFactoryClasses(py::module &m) return A.orthoCheck(knots); }, R"( - Check orthgonality of basis functions by quadrature + Check orthogonality of basis functions by quadrature Inner-product matrix of Sturm-Liouville solutions indexed by harmonic order used to assess fidelity. @@ -1960,7 +1960,7 @@ void BasisFactoryClasses(py::module &m) .def_static("invI", [](int I) { - if (I<0) std::runtime_error("I must be an interger greater than or equal to 0"); + if (I<0) std::runtime_error("I must be an integer greater than or equal to 0"); int l = std::floor(0.5*(-1.0 + std::sqrt(1.0 + 8.0 * I))); int m = I - int(l * (l + 1) / 2); return std::tuple(l, m); @@ -2030,7 +2030,7 @@ void BasisFactoryClasses(py::module &m) return A.orthoCheck(knots); }, R"( - Check orthgonality of basis functions by quadrature + Check orthogonality of basis functions by quadrature Inner-product matrix of Sturm-Liouville solutions indexed by harmonic order used to assess fidelity. @@ -2072,7 +2072,7 @@ void BasisFactoryClasses(py::module &m) return A.orthoCheck(knots); }, R"( - Check orthgonality of basis functions by quadrature + Check orthogonality of basis functions by quadrature Inner-product matrix of Sturm-Liouville solutions indexed by harmonic order used to assess fidelity. @@ -2142,7 +2142,7 @@ void BasisFactoryClasses(py::module &m) return A.orthoCheck(); }, R"( - Check orthgonality of basis functions by quadrature + Check orthogonality of basis functions by quadrature Inner-product matrix of Sturm-Liouville solutions indexed by harmonic order used to assess fidelity. @@ -2230,7 +2230,7 @@ void BasisFactoryClasses(py::module &m) return A.orthoCheck(); }, R"( - Check orthgonality of basis functions by quadrature + Check orthogonality of basis functions by quadrature Inner-product matrix of Sturm-Liouville solutions indexed by harmonic order used to assess fidelity. @@ -2328,7 +2328,7 @@ void BasisFactoryClasses(py::module &m) return A.orthoCheck(); }, R"( - Check orthgonality of basis functions by quadrature + Check orthogonality of basis functions by quadrature Inner-product matrix of indexed by flattened wave number (nx, ny, nz) where each of nx is in [-nmaxx, nmaxx], ny is in [-nmaxy, nmaxy] and nz is in @@ -2413,13 +2413,13 @@ void BasisFactoryClasses(py::module &m) return A.orthoCheck(); }, R"( - Check orthgonality of basis functions by quadrature + Check orthogonality of basis functions by quadrature Inner-product matrix of indexed by flattened wave number (nx, ny, nz) where each of nx is in [-nmaxx, nmaxx], and so on for ny and nz. Each dimension has dx=2*nmaxx+1 wave numbers and similarly for dy and dz. The index into the array is index=(nx+nmaxx)*dx*dy + (ny+nmaxy)*dy + (nz+nmaxz). This is an - analyic basis so the orthogonality matrix is not a check of any numerical + analytic basis so the orthogonality matrix is not a check of any numerical computation other than the quadrature itself. It is included for completeness. Parameters @@ -2652,7 +2652,7 @@ void BasisFactoryClasses(py::module &m) return A.orthoCheck(); }, R"( - Check orthgonality of basis functions by quadrature + Check orthogonality of basis functions by quadrature Inner-product matrix of orthogonal functions From b049baffef32868941b343888401de9119d0693b Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 15 Nov 2025 19:36:30 -0500 Subject: [PATCH 80/82] Typo in index bounds fixed --- expui/BiorthBasis.cc | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index e22de62e6..d04e91f10 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -3880,19 +3880,19 @@ namespace BasisClasses // Assign values from YAML // try { - if (conf["nminx"]) nminx = conf["nminx"].as(); - if (conf["nminy"]) nminy = conf["nminy"].as(); - if (conf["nminz"]) nminz = conf["nminz"].as(); + if (conf["nminx"]) nminx = conf["nminx" ].as(); + if (conf["nminy"]) nminy = conf["nminy" ].as(); + if (conf["nminz"]) nminz = conf["nminz" ].as(); - if (conf["nmaxx"]) nmaxx = conf["nmaxx"].as(); - if (conf["nmaxy"]) nmaxy = conf["nmaxy"].as(); - if (conf["nmaxz"]) nmaxz = conf["nmaxz"].as(); + if (conf["nmaxx"]) nmaxx = conf["nmaxx" ].as(); + if (conf["nmaxy"]) nmaxy = conf["nmaxy" ].as(); + if (conf["nmaxz"]) nmaxz = conf["nmaxz" ].as(); - if (conf["knots"]) knots = conf["knots"].as(); + if (conf["knots"]) knots = conf["knots" ].as(); - if (conf["check"]) check = conf["check"].as(); + if (conf["check"]) check = conf["check" ].as(); - if (conf["pcavar"]) pcavar = conf["pcavar"].as(); + if (conf["pcavar"]) pcavar = conf["pcavar" ].as(); if (conf["subsamp"]) sampT = conf["subsamp"].as(); } catch (YAML::Exception & error) { @@ -4055,7 +4055,7 @@ namespace BasisClasses expcoef(ix, iy, iz) += - mass * curr(0)*curr(1)*curr(2) * norm; if (pcavar) - g[index1D(ix, iy, iz)] = curr(0)*curr(1)*curr(2) * norm; + g[index1D(ix, iy, iz)] = - mass * curr(0)*curr(1)*curr(2) * norm; } } } @@ -4270,21 +4270,21 @@ namespace BasisClasses unsigned Cube::index1D(int kx, int ky, int kz) { - if (kx <0 or kx > 2*nmaxx) { + if (kx < 0 or kx > 2*nmaxx) { std::ostringstream sout; sout << "Cube::index1D: x index [" << kx << "] must be in [0, " << 2*nmaxx << "]"; throw std::runtime_error(sout.str()); } - if (ky <-nmaxy or ky > nmaxy) { + if (ky < 0 or ky > 2*nmaxy) { std::ostringstream sout; sout << "Cube::index1D: y index [" << ky << "] must be in [0, " << 2*nmaxy << "]"; throw std::runtime_error(sout.str()); } - if (kz <-nmaxz or kx > nmaxz) { + if (kz < 0 or kx > 2*nmaxz) { std::ostringstream sout; sout << "Cube::index1D: z index [" << kz << "] must be in [0, " << 2*nmaxz << "]"; From af32999cb7a6979353321527dda2c8f6760637b9 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 16 Nov 2025 15:36:25 -0500 Subject: [PATCH 81/82] Remove unused loop variable --- expui/BiorthBasis.cc | 148 +++++++++++++++++++++++++------------------ 1 file changed, 87 insertions(+), 61 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index d04e91f10..3607adf24 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -612,7 +612,7 @@ namespace BasisClasses Eigen::VectorXd workE; // M loop - for (int m=0, moffset=0, moffE=0; m<=l; m++, L++) { + for (int m=0, moffset=0; m<=l; m++, L++) { fac = factorial(l, m) * legs[tid](l, m); @@ -4055,7 +4055,7 @@ namespace BasisClasses expcoef(ix, iy, iz) += - mass * curr(0)*curr(1)*curr(2) * norm; if (pcavar) - g[index1D(ix, iy, iz)] = - mass * curr(0)*curr(1)*curr(2) * norm; + g[index1D(ix, iy, iz)] = - curr(0)*curr(1)*curr(2) * norm; } } } @@ -4068,7 +4068,7 @@ namespace BasisClasses sampleMasses(T) += mass; meanV[T].noalias() += g * mass; - covrV[T].noalias() += g * g.adjoint() * mass; + if (covar) covrV[T].noalias() += g * g.adjoint() * mass; } } @@ -4099,8 +4099,9 @@ namespace BasisClasses MPI_Allreduce(MPI_IN_PLACE, meanV[T].data(), meanV[T].size(), MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(MPI_IN_PLACE, covrV[T].data(), covrV[T].size(), - MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD); + if (covar) + MPI_Allreduce(MPI_IN_PLACE, covrV[T].data(), covrV[T].size(), + MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD); } // END: sample loop } @@ -4120,7 +4121,7 @@ namespace BasisClasses for (int T=0; T(ret[T][0]) = meanV[T]; - std::get<1>(ret[T][0]) = covrV[T]; + if (covar) std::get<1>(ret[T][0]) = covrV[T]; } } @@ -4243,10 +4244,12 @@ namespace BasisClasses v.resize(Itot); } - covrV.resize(sampT); - for (auto& v : covrV) { - v.resize(Itot, Itot); - } + if (covar) { + covrV.resize(sampT); + for (auto& v : covrV) { + v.resize(Itot, Itot); + } + } else covrV.clear(); sampleCounts.resize(sampT); sampleMasses.resize(sampT); @@ -4260,7 +4263,7 @@ namespace BasisClasses { for (int T=0; T(covar[0][0]).size(); + // Pack the coefficient data // if (floatType) { @@ -5131,23 +5139,26 @@ namespace BasisClasses // Pack the covariance data in an upper triangular format // - real_part.resize(ltot*diagonalSize*sampleSize); - imag_part.resize(ltot*diagonalSize*sampleSize); + if (varsz) { - for (size_t T=0, c=0; T(covar[T][l])(n1, n2)); - imag_part(c) = std::imag(std::get<1>(covar[T][l])(n1, n2)); + real_part.resize(ltot*diagonalSize*sampleSize); + imag_part.resize(ltot*diagonalSize*sampleSize); + + for (size_t T=0, c=0; T(covar[T][l])(n1, n2)); + imag_part(c) = std::imag(std::get<1>(covar[T][l])(n1, n2)); + } } } } - } - // Create two separate, compressed datasets - stanza.createDataSet("covariance_real", real_part, dcpl3); - stanza.createDataSet("covariance_imag", imag_part, dcpl3); + // Create two separate, compressed datasets + stanza.createDataSet("covariance_real", real_part, dcpl3); + stanza.createDataSet("covariance_imag", imag_part, dcpl3); + } } else { Eigen::VectorXd real_part(ltot*nmax*sampleSize); @@ -5169,24 +5180,26 @@ namespace BasisClasses // Pack the covariance data in an upper triangular format // - real_part.resize(ltot*diagonalSize*sampleSize); - imag_part.resize(ltot*diagonalSize*sampleSize); - - for (size_t T=0, c=0; T(covar[T][l])(n1, n2)); - imag_part(c) = std::imag(std::get<1>(covar[T][l])(n1, n2)); + if (varsz) { + real_part.resize(ltot*diagonalSize*sampleSize); + imag_part.resize(ltot*diagonalSize*sampleSize); + + for (size_t T=0, c=0; T(covar[T][l])(n1, n2)); + imag_part(c) = std::imag(std::get<1>(covar[T][l])(n1, n2)); + } } } } - } - // Create two separate, compressed datasets - // - stanza.createDataSet("covariance_real", real_part, dcpl3); - stanza.createDataSet("covariance_imag", imag_part, dcpl3); + // Create two separate, compressed datasets + // + stanza.createDataSet("covariance_real", real_part, dcpl3); + stanza.createDataSet("covariance_imag", imag_part, dcpl3); + } } // END: sample loop @@ -5436,16 +5449,21 @@ namespace BasisClasses data0.real() = data_real.cast(); data0.imag() = data_imag.cast(); - data_real = - stanza.getDataSet("covariance_real").read(); + // Check for existence of covariance + // + if (stanza.exist("covariance_real")) { + + data_real = + stanza.getDataSet("covariance_real").read(); - data_imag = - stanza.getDataSet("covariance_imag").read(); + data_imag = + stanza.getDataSet("covariance_imag").read(); - // Resize the complex array and assign - data1.resize(data_real.size()); - data1.real() = data_real.cast(); - data1.imag() = data_imag.cast(); + // Resize the complex array and assign + data1.resize(data_real.size()); + data1.real() = data_real.cast(); + data1.imag() = data_imag.cast(); + } } else { // Get the real and imaginary parts Eigen::VectorXd data_real = @@ -5459,17 +5477,22 @@ namespace BasisClasses data0.real() = data_real; data0.imag() = data_imag; - // Get the real and imaginary parts - data_real = - stanza.getDataSet("covariance_real").read(); + // Check for existence of covariance + // + if (stanza.exist("covariance_real")) { - data_imag = - stanza.getDataSet("covariance_imag").read(); - - // Resize the complex array and assign - data1.resize(data_real.size()); - data1.real() = data_real; - data1.imag() = data_imag; + // Get the real and imaginary parts + data_real = + stanza.getDataSet("covariance_real").read(); + + data_imag = + stanza.getDataSet("covariance_imag").read(); + + // Resize the complex array and assign + data1.resize(data_real.size()); + data1.real() = data_real; + data1.imag() = data_imag; + } } // Positions in data stanzas @@ -5484,7 +5507,7 @@ namespace BasisClasses // Coefficients std::get<0>(e).resize(rank); // Covariance matrix - std::get<1>(e).resize(rank, rank); + if (data1.size()) std::get<1>(e).resize(rank, rank); } // Pack the coefficient data @@ -5499,11 +5522,14 @@ namespace BasisClasses // Pack the covariance data c = 0; for (size_t l=0; l(elem[l])(n1, n2) = data1(sCov + c++); - if (n1 != n2) - std::get<1>(elem[l])(n2, n1) = std::get<1>(elem[l])(n1, n2); + + if (data1.size()) { + for (size_t n1=0; n1(elem[l])(n1, n2) = data1(sCov + c++); + if (n1 != n2) + std::get<1>(elem[l])(n2, n1) = std::get<1>(elem[l])(n1, n2); + } } } } From f3e41ae7e8e4a69045aaba38791e8f3188d640ef Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sun, 16 Nov 2025 15:45:20 -0500 Subject: [PATCH 82/82] Added covariance toggle to Cube --- expui/BiorthBasis.H | 6 +++++- pyEXP/BasisWrappers.cc | 27 ++++++++++++++++++++++----- 2 files changed, 27 insertions(+), 6 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 91c251d47..b3d21eb27 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -1333,7 +1333,8 @@ namespace BasisClasses void init_covariance(); void zero_covariance(); void writeCovarH5Params(HighFive::File& file); - int sampT = 100; + bool covar = false; + int sampT = 100; //@} public: @@ -1406,6 +1407,9 @@ namespace BasisClasses } } + //! Turn on/off covariance computation + void setCovar(bool flag=true) { covar = flag; } + //! Return wave-number indices from flattened index std::tuple index3D(unsigned k); diff --git a/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index be739b112..8bb11f8a1 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -885,7 +885,6 @@ void BasisFactoryClasses(py::module &m) PYBIND11_OVERRIDE(void, Cube, enableCoefCovariance, pcavar, nsamples); } - }; @@ -2366,6 +2365,25 @@ void BasisFactoryClasses(py::module &m) Cube the new instance )", py::arg("YAMLstring")) + .def("setCovar", &BasisClasses::Cube::setCovar, + R"( + Turn on coefficient covariance matrix computation. Set to OFF by default to + save disk space, computatation time, and memory. + + Parameters + ---------- + flag : bool + True to enable covariance matrix computation, False to disable it. + + Returns + ------- + None + + Notes + ----- + The covariance matrix computation can use significant memory and disk space + for the Cube basis especially and is disabled by default. + )", py::arg("flag")=true) .def("index1D", &BasisClasses::Cube::index1D, R"( Returns a flattened 1-d index into the arrays and matrices returned by the @@ -2880,8 +2898,8 @@ void BasisFactoryClasses(py::module &m) each subsample. The returns are complex-valued. )", py::arg("time")) - .def("getCovarSamples", static_cast - (BasisClasses::CovarianceReader::*)(unsigned)>(&BasisClasses::CovarianceReader::getCovarSamples), + .def("getCovarSamples", + py::overload_cast(&BasisClasses::CovarianceReader::getCovarSamples), R"( Get sample counts for the covariance computation @@ -2896,8 +2914,7 @@ void BasisFactoryClasses(py::module &m) sample counts and masses for the covariance computation )", py::arg("index")) - .def("getCovarSamples", static_cast - (BasisClasses::CovarianceReader::*)(double)>(&BasisClasses::CovarianceReader::getCovarSamples), + .def("getCovarSamples", py::overload_cast(&BasisClasses::CovarianceReader::getCovarSamples), R"( Get sample counts for the covariance computation