From 9e1583285f8776a61f5f0083f1d4c8987857edf7 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 24 Jun 2024 23:25:00 -0700 Subject: [PATCH 01/40] Add new construction to create instance from a cachefile [no ci] --- exputil/EmpCylSL.cc | 118 +++++++++++++++++++++++++++++++++++++++++++- include/EmpCylSL.H | 3 ++ 2 files changed, 120 insertions(+), 1 deletion(-) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 23d579c0f..7ded57acd 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -216,6 +216,122 @@ EmpCylSL::EmpCylSL(int nmax, int lmax, int mmax, int nord, maxSNR = 0.0; } +template +U getH5(std::string name, HighFive::File& file) +{ + if (file.hasAttribute(name)) { + U v; + HighFive::Attribute vv = file.getAttribute(name); + vv.read(v); + return v; + } else { + std::ostringstream sout; + sout << "EmpCylSL: could not find <" << name << ">"; + throw std::runtime_error(sout.str()); + } +} + +EmpCylSL::EmpCylSL(int mlim, std::string cachename) +{ + // Use default name? + // + if (cachename.size()==0) + throw std::runtime_error("EmpCylSL: you must specify a cachename"); + + // Open and read the cache file to get the needed input parameters + // + + try { + // Silence the HDF5 error stack + // + HighFive::SilenceHDF5 quiet; + + // Open the hdf5 file + // + HighFive::File file(cachefile, HighFive::File::ReadOnly); + + // For basis ID + std::string forceID("Cylinder"), geometry("cylinder"); + std::string model = EmpModelLabs[mtype]; + + // Version check + // + if (file.hasAttribute("Version")) { + auto ver = getH5(std::string("Version"), file); + if (ver.compare(Version)) + throw std::runtime_error("EmpCylSL: version mismatch"); + } else { + throw std::runtime_error("EmpCylSL: outdated cache"); + } + + NMAX = getH5("nmaxfid", file); + MMAX = getH5("mmax", file); + LMAX = getH5("lmax", file); + NORDER = getH5("nmax", file); + MMIN = 0; + MLIM = std::min(mlim, MMAX); + NMIN = 0; + NLIM = std::numeric_limits::max(); + Neven = getH5("neven", file); + Nodd = getH5("nodd", file); + ASCALE = getH5("ascl", file); + HSCALE = getH5("hscl", file); + } + catch (YAML::Exception& error) { + std::ostringstream sout; + sout << "EmpCylSL::getHeader: invalid cache file <" << cachefile << ">. "; + sout << "YAML error in getHeader: " << error.what() << std::endl; + throw GenericError(sout.str(), __FILE__, __LINE__, 1038, false); + } + + // Set EvenOdd if values seem sane + // + EvenOdd = false; + if (Nodd>=0 and Nodd<=NORDER and Nodd+Neven==NORDER) { + EvenOdd = true; + } + + pfac = 1.0/sqrt(ASCALE); + ffac = pfac/ASCALE; + dfac = ffac/ASCALE; + + EVEN_M = false; + + // Check whether MPI is initialized + // + int flag; + MPI_Initialized(&flag); + if (flag) use_mpi = true; + else use_mpi = false; + + // Enable MPI code in SLGridSph + // + if (use_mpi and numprocs>1) SLGridSph::mpi = 1; + + // Choose table dimension + // + MPItable = 4; + + // Initialize storage and values + // + coefs_made.resize(multistep+1); + std::fill(coefs_made.begin(), coefs_made.end(), false); + + eof_made = false; + + sampT = 1; + defSampT = 1; + tk_type = None; + + cylmass = 0.0; + cylmass1 = std::vector(nthrds); + cylmass_made = false; + + hallfile = ""; + minSNR = std::numeric_limits::max(); + maxSNR = 0.0; +} + void EmpCylSL::reset(int numr, int lmax, int mmax, int nord, double ascale, double hscale, int nodd, @@ -815,7 +931,7 @@ std::string compare_out(std::string str, U one, U two) return sout.str(); } -int EmpCylSL::cache_grid(int readwrite, string cachename) +int EmpCylSL::cache_grid(int readwrite, std::string cachename) { // Option to reset cache file name diff --git a/include/EmpCylSL.H b/include/EmpCylSL.H index 30a4a6587..233039ef7 100644 --- a/include/EmpCylSL.H +++ b/include/EmpCylSL.H @@ -501,6 +501,9 @@ public: double ascale, double hscale, int Nodd, std::string cachename); + //! Construct from cache file + EmpCylSL(int mlim, const std::string cache); + //! Destructor ~EmpCylSL(void); From 9cdc6e009243ff742d396096ad801ef85fa37dc6 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 25 Jun 2024 14:56:55 -0700 Subject: [PATCH 02/40] Fix for cache-only constructor [no ci] --- exputil/EmpCylSL.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 7ded57acd..f24bb9fc7 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -248,7 +248,7 @@ EmpCylSL::EmpCylSL(int mlim, std::string cachename) // Open the hdf5 file // - HighFive::File file(cachefile, HighFive::File::ReadOnly); + HighFive::File file(cachename, HighFive::File::ReadOnly); // For basis ID std::string forceID("Cylinder"), geometry("cylinder"); @@ -266,7 +266,7 @@ EmpCylSL::EmpCylSL(int mlim, std::string cachename) NMAX = getH5("nmaxfid", file); MMAX = getH5("mmax", file); - LMAX = getH5("lmax", file); + LMAX = getH5("lmaxfid", file); NORDER = getH5("nmax", file); MMIN = 0; MLIM = std::min(mlim, MMAX); From ecd09f2307da2575a251c8e7517c51405344623b Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 25 Jun 2024 17:26:17 -0700 Subject: [PATCH 03/40] Missing initialization of accumulation structures [no ci] --- exputil/EmpCylSL.cc | 2 ++ 1 file changed, 2 insertions(+) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index f24bb9fc7..d3d414f42 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -330,6 +330,8 @@ EmpCylSL::EmpCylSL(int mlim, std::string cachename) hallfile = ""; minSNR = std::numeric_limits::max(); maxSNR = 0.0; + + setup_accumulation(); } From 532c267ff4a341ca0d13d0d2e08a1246028303c4 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 25 Jun 2024 21:22:10 -0700 Subject: [PATCH 04/40] Call read_cache() to set up grids [no ci] --- exputil/EmpCylSL.cc | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index d3d414f42..a28022364 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -238,6 +238,8 @@ EmpCylSL::EmpCylSL(int mlim, std::string cachename) if (cachename.size()==0) throw std::runtime_error("EmpCylSL: you must specify a cachename"); + cachefile = cachename; + // Open and read the cache file to get the needed input parameters // @@ -331,7 +333,7 @@ EmpCylSL::EmpCylSL(int mlim, std::string cachename) minSNR = std::numeric_limits::max(); maxSNR = 0.0; - setup_accumulation(); + read_cache(); } From c7c631101c41ed97ac1e951af9d6c4bc8d3229a1 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 29 Jun 2024 14:52:21 -0700 Subject: [PATCH 05/40] Added constructor from cache file for SLGridSph [no ci] --- expui/CMakeLists.txt | 2 +- expui/Koopman.H | 2 +- expui/Koopman.cc | 2 +- expui/expMSSA.H | 15 +++ expui/expMSSA.cc | 276 +++++++++++++++++++++++++++++++++++++++++- exputil/EmpCylSL.cc | 2 + exputil/SLGridMP2.cc | 36 ++++++ include/SLGridMP2.H | 9 ++ pyEXP/MSSAWrappers.cc | 43 ++++++- 9 files changed, 380 insertions(+), 7 deletions(-) diff --git a/expui/CMakeLists.txt b/expui/CMakeLists.txt index 8690817f7..d506bea8d 100644 --- a/expui/CMakeLists.txt +++ b/expui/CMakeLists.txt @@ -34,7 +34,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) + Koopman.cc BiorthBess.cc SGSmooth.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/Koopman.H b/expui/Koopman.H index a84b31dc7..921ac4089 100644 --- a/expui/Koopman.H +++ b/expui/Koopman.H @@ -119,7 +119,7 @@ namespace MSSA return L; } - //! Return the EDMD modes + //! Return the EDMD modes, an approximation to the Koopman eigenfunctions Eigen::MatrixXcd getModes() { if (not computed) koopman_analysis(); diff --git a/expui/Koopman.cc b/expui/Koopman.cc index b32f00b92..cc85cca1f 100644 --- a/expui/Koopman.cc +++ b/expui/Koopman.cc @@ -194,7 +194,7 @@ namespace MSSA { Eigen::VectorXcd B = Phi.inverse() * xx; Eigen::MatrixXcd LL = I.asDiagonal(); - // Propate the solution with the operator + // Propagate the solution using the approximate Koopman operator // for (int i=0; i + getKoopmanModes(const double tol, int window, bool debug); + + //! Return the reconstructed Koopman modes + std::map getReconstructedKoopman(int mode); + }; diff --git a/expui/expMSSA.cc b/expui/expMSSA.cc index 0aca81d71..3ff1c71c7 100644 --- a/expui/expMSSA.cc +++ b/expui/expMSSA.cc @@ -1536,13 +1536,283 @@ namespace MSSA { reconstructed = true; } - } catch (HighFive::Exception& err) { + } catch (HighFive::Exception& err) { // Number of channels + // + nkeys = mean.size(); + + // Make sure parameters are sane + // + if (numW<=0) numW = numT/2; + if (numW > numT/2) numW = numT/2; + + numK = numT - numW + 1; + + std::cerr << "**** Error opening or reading H5 file ****" << std::endl; throw; } } + std::tuple + expMSSA::getKoopmanModes(double tol, int D, bool debug) + { + bool use_fullKh = true; // Use the non-reduced computation of + // Koopman/eDMD + // Number of channels + // + nkeys = mean.size(); + + // Make sure parameters are sane + // + if (numW<=0) numW = numT/2; + if (numW > numT/2) numW = numT/2; + + numK = numT - numW + 1; + + Eigen::VectorXd S1; + Eigen::MatrixXd Y1; + Eigen::MatrixXd V1; + Eigen::MatrixXd VT1; + Eigen::MatrixXd VT2; + + // Make a new trajetory matrix with smoothing + // + Y1.resize(numK, numW*nkeys + D*(nkeys-1)); + + size_t n=0, offset=0; + for (auto k : mean) { + for (int i=0; i 0) { + // Back blending + for (int j=0; j(D-j)/D; + } + } + // Main series + for (int j=0; j(D-j)/D; + } + } + } + offset += numW + D; + n++; + } + + double Scale = Y1.norm(); + + // auto YY = Y1/Scale; + auto YY = Y1; + + // Use one of the built-in Eigen3 algorithms + // + /* + if (params["Jacobi"]) { + // -->Using Jacobi + Eigen::JacobiSVD + svd(YY, Eigen::ComputeThinU | Eigen::ComputeThinV); + S1 = svd.singularValues(); + V1 = svd.matrixV(); + } else if (params["BDCSVD"]) { + */ + // -->Using BDC + Eigen::BDCSVD + svd(YY, Eigen::ComputeFullU | Eigen::ComputeFullV); + // svd(YY, Eigen::ComputeThinU | Eigen::ComputeThinV); + S1 = svd.singularValues(); + V1 = svd.matrixV(); + /* + } else { + // -->Use Random approximation algorithm from Halko, Martinsson, + // and Tropp + int srank = std::min(YY.cols(), YY.rows()); + RedSVD::RedSVD svd(YY, srank); + S1 = svd.singularValues(); + V1 = svd.matrixV(); + } + */ + + std::cout << "shape V1 = " << V1.rows() << " x " + << V1.cols() << std::endl; + + std::cout << "shape Y1 = " << Y1.rows() << " x " + << Y1.cols() << std::endl; + + int lags = V1.rows(); + int rank = V1.cols(); + + std::ofstream out; + if (debug) out.open("debug.txt"); + + if (out) out << "rank=" << rank << " lags=" << lags << std::endl; + + VT1.resize(rank, lags-1); + VT2.resize(rank, lags-1); + + for (int j=0; j uu; + for (int i=0; iUsing Jacobi + Eigen::JacobiSVD + // svd(VT1, Eigen::ComputeThinU | Eigen::ComputeThinV); + svd(VT1, Eigen::ComputeFullU | Eigen::ComputeFullV); + SS = svd.singularValues(); + UU = svd.matrixU(); + VV = svd.matrixV(); + } else if (params["BDCSVD"]) { + */ + { + // -->Using BDC + Eigen::BDCSVD + // svd(VT1, Eigen::ComputeThinU | Eigen::ComputeThinV); + svd(VT1, Eigen::ComputeFullU | Eigen::ComputeFullV); + SS = svd.singularValues(); + UU = svd.matrixU(); + VV = svd.matrixV(); + } + /* + } else { + // -->Use Random approximation algorithm from Halko, Martinsson, + // and Tropp + // RedSVD::RedSVD svd(VT1, std::min(rank, numK-1)); + RedSVD::RedSVD svd(VT1, std::max(VT1.rows(), VT2.cols())); + SS = svd.singularValues(); + UU = svd.matrixU(); + VV = svd.matrixV(); + } + */ + + if (out) out << "Singular values" << std::endl << SS << std::endl; + + // Compute inverse + for (int i=0; i tol) SS(i) = 1.0/SS(i); + else SS(i) = 0.0; + } + + // Compute full Koopman operator + if (use_fullKh) { + + Eigen::MatrixXd DD(VV.cols(), UU.cols()); + DD.setZero(); + for (int i=0; i es(AT); + + L = es.eigenvalues(); + Phi = es.eigenvectors(); + + if (out) { + out << std::endl << "Eigenvalues" << std::endl << L << std::endl + << std::endl << "Eigenvectors" << std::endl << Phi << std::endl; + } + + } + // Compute the reduced Koopman operator + else { + + Eigen::MatrixXd AT = UU.transpose() * (VT2 * VV) * SS.asDiagonal(); + + // Compute spectrum + Eigen::EigenSolver es(AT, true); + + L = es.eigenvalues(); + auto W = es.eigenvectors(); + + // Compute the EDMD modes + // + Eigen::VectorXcd Linv(L); + for (int i=0; i tol) Linv(i) = 1.0/Linv(i); + else Linv(i) = 0.0; + } + + Phi = VT2 * VV * SS.asDiagonal() * W * Linv.asDiagonal(); + + if (out) { + out << std::endl << "Eigenvalues" << std::endl << L << std::endl + << std::endl << "Eigenvectors" << std::endl << Phi << std::endl; + } + } + + // Cache window size + // + window = D; + + return {L, Phi}; + } + + std::map + expMSSA::getReconstructedKoopman(int mode) + { + // Copy the original map for return + // + auto newdata = data; + + size_t n=0, offset=0; + + for (auto u : mean) { + + double disp = totVar; + if (type == TrendType::totPow) disp = totPow; + if (disp==0.0) disp = var[u.first]; + + std::complex phase = 1.0; + for (int i=0; i(); - std::cout << "Trajectory is " << std::boolalpha << trajectory - << std::endl; + // std::cout << "Trajectory is " << std::boolalpha << trajectory + // << std::endl; // Eigen OpenMP reporting // diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index a28022364..e69be02aa 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -270,6 +270,8 @@ EmpCylSL::EmpCylSL(int mlim, std::string cachename) MMAX = getH5("mmax", file); LMAX = getH5("lmaxfid", file); NORDER = getH5("nmax", file); + CMAPR = getH5("cmapr", file); + CMAPZ = getH5("cmapz", file); MMIN = 0; MLIM = std::min(mlim, MMAX); NMIN = 0; diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc index 0984fff00..7924551b7 100644 --- a/exputil/SLGridMP2.cc +++ b/exputil/SLGridMP2.cc @@ -168,6 +168,42 @@ SLGridSph::SLGridSph(std::shared_ptr mod, } +SLGridSph::SLGridSph(std::string cachename) +{ + if (cachename.size()) sph_cache_name = cachename; + else throw std::runtime_error("SLGridSph: you must specify a cachename"); + + tbdbg = false; + + int LMAX, NMAX, NUMR, CMAP, DIVERGE=0; + double RMIN, RMAX, RMAP, DFAC=1.0; + + try { + + auto node = getHeader(cachename); + + LMAX = node["lmax"].as(); + NMAX = node["nmax"].as(); + NUMR = node["numr"].as(); + CMAP = node["cmap"].as(); + RMIN = node["rmin"].as(); + RMAX = node["rmax"].as(); + RMAP = node["rmapping"].as(); + + model_file_name = node["model"].as(); + model = SphModTblPtr(new SphericalModelTable(model_file_name, diverge, dfac)); + } + catch (YAML::Exception& error) { + std::ostringstream sout; + sout << "SLGridMP2: error parsing parameters from getHeader: " + << error.what(); + throw GenericError(sout.str(), __FILE__, __LINE__, 1039, false); + } + + initialize(LMAX, NMAX, NUMR, RMIN, RMAX, false, CMAP, RMAP); +} + + std::map SLGridSph::cacheInfo(const std::string& cachefile, bool verbose) { diff --git a/include/SLGridMP2.H b/include/SLGridMP2.H index 95a9bf6b7..f83d98fe9 100644 --- a/include/SLGridMP2.H +++ b/include/SLGridMP2.H @@ -116,6 +116,10 @@ public: std::string cachename=".slgrid_sph_cache", bool Verbose=false); + //! Constructor (from cache file) + SLGridSph(std::string cachename); + + //! Destructor virtual ~SLGridSph(); @@ -216,6 +220,11 @@ public: double getRmax() { return rmax; } //@} + //@{ + //! Get expansion limits + int getLmax() { return lmax; } + int getNmax() { return nmax; } + //@} }; diff --git a/pyEXP/MSSAWrappers.cc b/pyEXP/MSSAWrappers.cc index 9fdc5658f..aa5ef5c57 100644 --- a/pyEXP/MSSAWrappers.cc +++ b/pyEXP/MSSAWrappers.cc @@ -326,7 +326,6 @@ void MSSAtoolkitClasses(py::module &m) { dict({id: Coefs},...) reconstructed time series in the original coefficient form - Notes ----- The reconstructed data will overwrite the memory of the original coefficient @@ -543,6 +542,48 @@ void MSSAtoolkitClasses(py::module &m) { )"); + f.def("getKoopmanModes", &expMSSA::getKoopmanModes, + R"( + Compute the Koopman mode estimate from the right-singular vectors + + Uses eDMD to estimate the modes + + Parameters + ---------- + tol : double + singular value truncation level + window: int + Smoothing between serialized channels (0 for no smoothing) + debug : bool + flag indicating whether to print debug information + + Notes + ----- + Use getReconstructedKoopman() to copy the reconstruction for a + particular mode back to the coefficient db + + Returns + ------- + tuple(numpy.ndarray, numpy.ndarray) + vector of eigenvalues and modes + )", py::arg("tol")=1.0e-12, py::arg("window")=0, py::arg("debug")=false); + + f.def("getReconstructedKoopman", &expMSSA::getReconstructedKoopman, + R"( + Reconstruct the coefficients for a particular Koopman mode + + Parameters + ---------- + mode: int + The index of the mode to be reconstructed + + Returns + ------- + dict({id: Coefs},...) + reconstructed time series in the original coefficient form + + )", py::arg("mode")); + f.def("getRC", &expMSSA::getRC, R"( Access the detrended reconstructed channel series by internal key From 60957b5ac0c761ee5b1144920970b97b9fe37b4d Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 29 Jun 2024 14:55:48 -0700 Subject: [PATCH 06/40] Added Savitsky-Golay smoothing class [no ci] --- expui/SGSmooth.H | 13 ++ expui/SGSmooth.cc | 549 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 562 insertions(+) create mode 100644 expui/SGSmooth.H create mode 100644 expui/SGSmooth.cc diff --git a/expui/SGSmooth.H b/expui/SGSmooth.H new file mode 100644 index 000000000..6615bffa2 --- /dev/null +++ b/expui/SGSmooth.H @@ -0,0 +1,13 @@ +#ifndef __SGSMOOTH_HPP__ +#define __SGSMOOTH_HPP__ + +#include + +//! savitzky golay smoothing. +std::vector sg_smooth(const std::vector &v, const int w, const int deg); + +//! numerical derivative based on savitzky golay smoothing. +std::vector sg_derivative(const std::vector &v, const int w, + const int deg, const double h=1.0); + +#endif // __SGSMOOTH_HPP__ diff --git a/expui/SGSmooth.cc b/expui/SGSmooth.cc new file mode 100644 index 000000000..1b5150b53 --- /dev/null +++ b/expui/SGSmooth.cc @@ -0,0 +1,549 @@ +//! +// Sliding window signal processing (and linear algebra toolkit). +// +// supported operations: +//
    +//
  • Savitzky-Golay smoothing. +//
  • computing a numerical derivative based of Savitzky-Golay smoothing. +//
  • required linear algebra support for SG smoothing using STL based +// vector/matrix classes +//
+// +// \brief Linear Algebra "Toolkit". +// +// modified by Rob Patro, 2016 +// modified by MDW, 2024 + +// system headers +#include +#include +#include // for size_t +#include // for fabs +#include + +//! default convergence +static const double TINY_FLOAT = 1.0e-300; + +//! comfortable array of doubles +using float_vect = std::vector; +//! comfortable array of ints; +using int_vect = std::vector; + +/*! matrix class. + * + * This is a matrix class derived from a vector of float_vects. Note that + * the matrix elements indexed [row][column] with indices starting at 0 (c + * style). Also note that because of its design looping through rows should + * be faster than looping through columns. + * + * \brief two dimensional floating point array + */ +class float_mat : public std::vector { +private: + //! disable the default constructor + explicit float_mat() {}; + //! disable assignment operator until it is implemented. + float_mat &operator =(const float_mat &) { return *this; }; +public: + //! constructor with sizes + float_mat(const size_t rows, const size_t cols, const double def=0.0); + //! copy constructor for matrix + float_mat(const float_mat &m); + //! copy constructor for vector + float_mat(const float_vect &v); + + //! use default destructor + // ~float_mat() {}; + + //! get size + size_t nr_rows(void) const { return size(); }; + //! get size + size_t nr_cols(void) const { return front().size(); }; +}; + + + +// constructor with sizes +float_mat::float_mat(const size_t rows,const size_t cols,const double defval) + : std::vector(rows) { + int i; + for (i = 0; i < rows; ++i) { + (*this)[i].resize(cols, defval); + } + if ((rows < 1) || (cols < 1)) { + std::ostringstream msg; + msg << "cannot build matrix with " << rows << " rows and " << cols + << "columns"; + throw std::runtime_error(msg.str()); + } +} + +// copy constructor for matrix +float_mat::float_mat(const float_mat &m) : std::vector(m.size()) { + + float_mat::iterator inew = begin(); + float_mat::const_iterator iold = m.begin(); + for (/* empty */; iold < m.end(); ++inew, ++iold) { + const size_t oldsz = iold->size(); + inew->resize(oldsz); + const float_vect oldvec(*iold); + *inew = oldvec; + } +} + +// copy constructor for vector +float_mat::float_mat(const float_vect &v) + : std::vector(1) { + + const size_t oldsz = v.size(); + front().resize(oldsz); + front() = v; +} + +////////////////////// +// Helper functions // +////////////////////// + +//! permute() orders the rows of A to match the integers in the index array. +void permute(float_mat &A, int_vect &idx) +{ + int_vect i(idx.size()); + int j,k; + + for (j = 0; j < A.nr_rows(); ++j) { + i[j] = j; + } + + // loop over permuted indices + for (j = 0; j < A.nr_rows(); ++j) { + if (i[j] != idx[j]) { + + // search only the remaining indices + for (k = j+1; k < A.nr_rows(); ++k) { + if (i[k] ==idx[j]) { + std::swap(A[j],A[k]); // swap the rows and + i[k] = i[j]; // the elements of + i[j] = idx[j]; // the ordered index. + break; // next j + } + } + } + } +} + +/*! \brief Implicit partial pivoting. + * + * The function looks for pivot element only in rows below the current + * element, A[idx[row]][column], then swaps that row with the current one in + * the index map. The algorithm is for implicit pivoting (i.e., the pivot is + * chosen as if the max coefficient in each row is set to 1) based on the + * scaling information in the vector scale. The map of swapped indices is + * recorded in swp. The return value is +1 or -1 depending on whether the + * number of row swaps was even or odd respectively. */ +static int partial_pivot(float_mat &A, const size_t row, const size_t col, + float_vect &scale, int_vect &idx, double tol) +{ + if (tol <= 0.0) + tol = TINY_FLOAT; + + int swapNum = 1; + + // default pivot is the current position, [row,col] + int pivot = row; + double piv_elem = fabs(A[idx[row]][col]) * scale[idx[row]]; + + // loop over possible pivots below current + int j; + for (j = row + 1; j < A.nr_rows(); ++j) { + + const double tmp = fabs(A[idx[j]][col]) * scale[idx[j]]; + + // if this elem is larger, then it becomes the pivot + if (tmp > piv_elem) { + pivot = j; + piv_elem = tmp; + } + } + +#if 0 + if (piv_elem < tol) { + throw std::runtime_error("partial_pivot(): Zero pivot encountered."); +#endif + + if(pivot > row) { // bring the pivot to the diagonal + j = idx[row]; // reorder swap array + idx[row] = idx[pivot]; + idx[pivot] = j; + swapNum = -swapNum; // keeping track of odd or even swap + } + return swapNum; +} + +/*! \brief Perform backward substitution. + * + * Solves the system of equations A*b=a, ASSUMING that A is upper + * triangular. If diag==1, then the diagonal elements are additionally + * assumed to be 1. Note that the lower triangular elements are never + * checked, so this function is valid to use after a LU-decomposition in + * place. A is not modified, and the solution, b, is returned in a. */ +static void lu_backsubst(float_mat &A, float_mat &a, bool diag=false) +{ + int r,c,k; + + for (r = (A.nr_rows() - 1); r >= 0; --r) { + for (c = (A.nr_cols() - 1); c > r; --c) { + for (k = 0; k < A.nr_cols(); ++k) { + a[r][k] -= A[r][c] * a[c][k]; + } + } + if(!diag) { + for (k = 0; k < A.nr_cols(); ++k) { + a[r][k] /= A[r][r]; + } + } + } +} + +/*! \brief Perform forward substitution. + * + * Solves the system of equations A*b=a, ASSUMING that A is lower + * triangular. If diag==1, then the diagonal elements are additionally + * assumed to be 1. Note that the upper triangular elements are never + * checked, so this function is valid to use after a LU-decomposition in + * place. A is not modified, and the solution, b, is returned in a. */ +static void lu_forwsubst(float_mat &A, float_mat &a, bool diag=true) +{ + int r,k,c; + for (r = 0;r < A.nr_rows(); ++r) { + for(c = 0; c < r; ++c) { + for (k = 0; k < A.nr_cols(); ++k) { + a[r][k] -= A[r][c] * a[c][k]; + } + } + if(!diag) { + for (k = 0; k < A.nr_cols(); ++k) { + a[r][k] /= A[r][r]; + } + } + } +} + +/*! \brief Performs LU factorization in place. + * + * This is Crout's algorithm (cf., Num. Rec. in C, Section 2.3). The map of + * swapped indeces is recorded in idx. The return value is +1 or -1 + * depending on whether the number of row swaps was even or odd + * respectively. idx must be preinitialized to a valid set of indices + * (e.g., {1,2, ... ,A.nr_rows()}). */ +static int lu_factorize(float_mat &A, int_vect &idx, double tol=TINY_FLOAT) +{ + if ( tol <= 0.0) + tol = TINY_FLOAT; + + if ((A.nr_rows() == 0) || (A.nr_rows() != A.nr_cols())) { + throw std::runtime_error("lu_factorize(): cannot handle empty " + "or nonsquare matrices."); + } + + float_vect scale(A.nr_rows()); // implicit pivot scaling + int i,j; + for (i = 0; i < A.nr_rows(); ++i) { + double maxval = 0.0; + for (j = 0; j < A.nr_cols(); ++j) { + if (fabs(A[i][j]) > maxval) + maxval = fabs(A[i][j]); + } + if (maxval == 0.0) { + throw std::runtime_error("lu_factorize(): zero pivot found."); + } + scale[i] = 1.0 / maxval; + } + + int swapNum = 1; + int c,r; + for (c = 0; c < A.nr_cols() ; ++c) { // loop over columns + swapNum *= partial_pivot(A, c, c, scale, idx, tol); // bring pivot to diagonal + for(r = 0; r < A.nr_rows(); ++r) { // loop over rows + int lim = (r < c) ? r : c; + for (j = 0; j < lim; ++j) { + A[idx[r]][c] -= A[idx[r]][j] * A[idx[j]][c]; + } + if (r > c) + A[idx[r]][c] /= A[idx[c]][c]; + } + } + permute(A,idx); + return swapNum; +} + +/*! \brief Solve a system of linear equations. + * Solves the inhomogeneous matrix problem with lu-decomposition. Note that + * inversion may be accomplished by setting a to the identity_matrix. */ +static float_mat lin_solve(const float_mat &A, const float_mat &a, + double tol=TINY_FLOAT) +{ + float_mat B(A); + float_mat b(a); + int_vect idx(B.nr_rows()); + int j; + + for (j = 0; j < B.nr_rows(); ++j) { + idx[j] = j; // init row swap label array + } + lu_factorize(B,idx,tol); // get the lu-decomp. + permute(b,idx); // sort the inhomogeneity to match the lu-decomp + lu_forwsubst(B,b); // solve the forward problem + lu_backsubst(B,b); // solve the backward problem + return b; +} + +/////////////////////// +// related functions // +/////////////////////// + +//! Returns the inverse of a matrix using LU-decomposition. +static float_mat invert(const float_mat &A) +{ + const int n = A.size(); + float_mat E(n, n, 0.0); + float_mat B(A); + int i; + + for (i = 0; i < n; ++i) { + E[i][i] = 1.0; + } + + return lin_solve(B, E); +} + +//! returns the transposed matrix. +static float_mat transpose(const float_mat &a) +{ + float_mat res(a.nr_cols(), a.nr_rows()); + int i,j; + + for (i = 0; i < a.nr_rows(); ++i) { + for (j = 0; j < a.nr_cols(); ++j) { + res[j][i] = a[i][j]; + } + } + return res; +} + +//! matrix multiplication. +float_mat operator *(const float_mat &a, const float_mat &b) +{ + float_mat res(a.nr_rows(), b.nr_cols()); + if (a.nr_cols() != b.nr_rows()) { + throw std::runtime_error("incompatible matrices in multiplication"); + } + + int i,j,k; + + for (i = 0; i < a.nr_rows(); ++i) { + for (j = 0; j < b.nr_cols(); ++j) { + double sum(0.0); + for (k = 0; k < a.nr_cols(); ++k) { + sum += a[i][k] * b[k][j]; + } + res[i][j] = sum; + } + } + return res; +} + + +//! calculate savitzky golay coefficients. +static float_vect sg_coeff(const float_vect &b, const size_t deg) +{ + const size_t rows(b.size()); + const size_t cols(deg + 1); + float_mat A(rows, cols); + float_vect res(rows); + + // generate input matrix for least squares fit + int i,j; + for (i = 0; i < rows; ++i) { + for (j = 0; j < cols; ++j) { + A[i][j] = pow(double(i), double(j)); + } + } + + float_mat c(invert(transpose(A) * A) * (transpose(A) * transpose(b))); + + for (i = 0; i < b.size(); ++i) { + res[i] = c[0][0]; + for (j = 1; j <= deg; ++j) { + res[i] += c[j][0] * pow(double(i), double(j)); + } + } + return res; +} + +/*! \brief savitzky golay smoothing. + * + * This method means fitting a polynome of degree 'deg' to a sliding window + * of width 2w+1 throughout the data. The needed coefficients are + * generated dynamically by doing a least squares fit on a "symmetric" unit + * vector of size 2w+1, e.g. for w=2 b=(0,0,1,0,0). evaluating the polynome + * yields the sg-coefficients. at the border non symmectric vectors b are + * used. */ +float_vect sg_smooth(const float_vect &v, const int width, const int deg) +{ + float_vect res(v.size(), 0.0); + if ((width < 1) || (deg < 0) || (v.size() < (2 * width + 2))) { + throw std::runtime_error("sgsmooth: parameter error."); + } + + const int window = 2 * width + 1; + const int endidx = v.size() - 1; + + // do a regular sliding window average + int i,j; + if (deg == 0) { + // handle border cases first because we need different coefficients +#if defined(_OPENMP) +#pragma omp parallel for private(i,j) schedule(static) +#endif + for (i = 0; i < width; ++i) { + const double scale = 1.0/double(i+1); + const float_vect c1(width, scale); + for (j = 0; j <= i; ++j) { + res[i] += c1[j] * v[j]; + res[endidx - i] += c1[j] * v[endidx - j]; + } + } + + // now loop over rest of data. reusing the "symmetric" coefficients. + const double scale = 1.0/double(window); + const float_vect c2(window, scale); +#if defined(_OPENMP) +#pragma omp parallel for private(i,j) schedule(static) +#endif + for (i = 0; i <= (v.size() - window); ++i) { + for (j = 0; j < window; ++j) { + res[i + width] += c2[j] * v[i + j]; + } + } + return res; + } + + // handle border cases first because we need different coefficients +#if defined(_OPENMP) +#pragma omp parallel for private(i,j) schedule(static) +#endif + for (i = 0; i < width; ++i) { + float_vect b1(window, 0.0); + b1[i] = 1.0; + + const float_vect c1(sg_coeff(b1, deg)); + for (j = 0; j < window; ++j) { + res[i] += c1[j] * v[j]; + res[endidx - i] += c1[j] * v[endidx - j]; + } + } + + // now loop over rest of data. reusing the "symmetric" coefficients. + float_vect b2(window, 0.0); + b2[width] = 1.0; + const float_vect c2(sg_coeff(b2, deg)); + +#if defined(_OPENMP) +#pragma omp parallel for private(i,j) schedule(static) +#endif + for (i = 0; i <= (v.size() - window); ++i) { + for (j = 0; j < window; ++j) { + res[i + width] += c2[j] * v[i + j]; + } + } + return res; +} + +/*! least squares fit a polynome of degree 'deg' to data in 'b'. + * then calculate the first derivative and return it. */ +static float_vect lsqr_fprime(const float_vect &b, const int deg) +{ + const int rows(b.size()); + const int cols(deg + 1); + float_mat A(rows, cols); + float_vect res(rows); + + // generate input matrix for least squares fit + int i,j; + for (i = 0; i < rows; ++i) { + for (j = 0; j < cols; ++j) { + A[i][j] = pow(double(i), double(j)); + } + } + + float_mat c(invert(transpose(A) * A) * (transpose(A) * transpose(b))); + + for (i = 0; i < b.size(); ++i) { + res[i] = c[1][0]; + for (j = 1; j < deg; ++j) { + res[i] += c[j + 1][0] * double(j+1) + * pow(double(i), double(j)); + } + } + return res; +} + +/*! \brief savitzky golay smoothed numerical derivative. + * + * This method means fitting a polynome of degree 'deg' to a sliding window + * of width 2w+1 throughout the data. + * + * In contrast to the sg_smooth function we do a brute force attempt by + * always fitting the data to a polynome of degree 'deg' and using the + * result. */ +float_vect sg_derivative(const float_vect &v, const int width, + const int deg, const double h) +{ + float_vect res(v.size(), 0.0); + if ((width < 1) || (deg < 1) || (v.size() < (2 * width + 2))) { + throw std::runtime_error("sgsderiv: parameter error"); + } + + const int window = 2 * width + 1; + + // handle border cases first because we do not repeat the fit + // lower part + float_vect b(window, 0.0); + int i,j; + + for (i = 0; i < window; ++i) { + b[i] = v[i] / h; + } + const float_vect c(lsqr_fprime(b, deg)); + for (j = 0; j <= width; ++j) { + res[j] = c[j]; + } + // upper part. direction of fit is reversed + for (i = 0; i < window; ++i) { + b[i] = v[v.size() - 1 - i] / h; + } + const float_vect d(lsqr_fprime(b, deg)); + for (i = 0; i <= width; ++i) { + res[v.size() - 1 - i] = -d[i]; + } + + // now loop over rest of data. wasting a lot of least squares calcs + // since we only use the middle value. +#if defined(_OPENMP) +#pragma omp parallel for private(i,j) schedule(static) +#endif + for (i = 1; i < (v.size() - window); ++i) { + for (j = 0; j < window; ++j) { + b[j] = v[i + j] / h; + } + res[i + width] = lsqr_fprime(b, deg)[width]; + } + return res; +} + +// Local Variables: +// mode: c++ +// c-basic-offset: 4 +// fill-column: 76 +// indent-tabs-mode: nil +// End: From 6825d70bed149663d0e3b8738bd23c5e28e9f9fc Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 2 Jul 2024 16:54:25 -0700 Subject: [PATCH 07/40] A draft CONTRIBUTING document for the EXP source tree ported from the readthedocs version --- CONTRIBUTING.md | 240 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 240 insertions(+) create mode 100644 CONTRIBUTING.md diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 000000000..3e5c6ac52 --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,240 @@ +# Contributing to EXP + +There are many ways to contribute to EXP. Here are some of them: + +- Blog about EXP. Cite the EXP published papers. Tell the world how + you're using EXP. This will help newcomers with more examples and + will help the EXP project to increase its visibility. +- Report bugs and request features in the [issue + tracker](https://github.com/EXP-code/EXP/issues), trying to follow + the guidelines detailed in **Reporting bugs** below. +- Submit new examples to the [EXP examples + repo](https://github.com/EXP-code/EXP-examples) or the [pyEXP + examples repo](https://github.com/EXP-code/pyEXP-examples). + +# Stable and development branches + +Our current branch policy changed as of May 1, 2023. Rather than a +`main` branch for a stable EXP/pyEXP release with development taking +place on the `devel` branch, the only official EXP branch will be +`main`. All development will take place through a [pull +request](https://help.github.com/en/github/collaborating-with-issues-and-pull-requests/creating-a-pull-request) +to `main`. All other branches are considered to be temporary. The +current stable release will be tagged in the GitHub release menu. The +HEAD of `main` will contain the latest new features and fixes. + +# Reporting bugs + +Well-written bug reports are very helpful, so keep in mind the following +guidelines when you're going to report a new bug. + +- check the [FAQ](https://exp-docs.readthedocs.io) first to see if + your issue is addressed in a well-known question +- check the [open issues](https://github.com/EXP-code/EXP/issues) + to see if the issue has already been reported. If it has, don't + dismiss the report, but check the ticket history and comments. If + you have additional useful information, please leave a comment, or + consider sending a pull request with a fix. +- write **complete, reproducible, specific bug reports**. The smaller + the test case, the better. Remember that other developers won't + have your project to reproduce the bug, so please include all + relevant files required to reproduce it. +- the most awesome way to provide a complete reproducible example is + to send a pull request which adds a failing test case to the EXP + testing suite. This is helpful even if you don't have an + intention to fix the issue yourselves. +- include the output of `exp -v` so developers working on your bug + know exactly which version and platform it occurred on, which is + often very helpful for reproducing it, or knowing if it was already + fixed. + +# Contributing to code development + +We are now using the [pull +request](https://help.github.com/en/github/collaborating-with-issues-and-pull-requests/creating-a-pull-request) +method for all EXP development, in including the code authors and +maintainers. In essence, a pull request is code patch that allows us +to easily review, test, and discuss the proposed change. The better a +patch is written, the higher the chances that it'll get accepted and +the sooner it will be merged. + +Well-written patches should: + +- contain the minimum amount of code required for the specific change. + Small patches are easier to review and merge. So, if you're doing + more than one change (or bug fix), please consider submitting one + patch per change. Do not collapse multiple changes into a single + patch. For big changes consider using a patch queue. +- pass all EXP basic example tests. See **Running tests** below. +- if you're contributing a feature, especially one that changes a + public (documented) API, please include the documentation changes + in the same patch. See **Documentation strategy** below. + +# Submitting patches + +The best way to submit a patch is to issue a [pull +request](https://help.github.com/en/github/collaborating-with-issues-and-pull-requests/creating-a-pull-request) +on GitHub. The work flow for this is: + +1. Fork the [EXP-code/EXP](https://github.com/EXP-code/EXP) repo in GitHub +2. Create a branch with the proposed change +3. Compile and test your changes +4. Once you are satisfied, create the [pull +request](https://help.github.com/en/github/collaborating-with-issues-and-pull-requests/creating-a-pull-request) +on the GitHub EXP code repo + +Remember to explain what was fixed or the new functionality (what it +is, why it's needed, etc). The more info you include, the easier will +be for core developers to understand and accept your patch. + +You can also discuss the new functionality (or bug fix) before creating +the patch, but it's always good to have a patch ready to illustrate +your arguments and show that you have put some additional thought into +the subject. A good starting point is to send a pull request on GitHub. +It can be simple enough to illustrate your idea, and leave +documentation/tests for later, after the idea has been validated and +proven useful. All functionality (including new features and bug fixes) +must include a test case to check that it works as expected, so please +include tests for your patches if you want them to get accepted sooner. + +There might be an existing pull request for the problem you'd like to +solve, so please do take a look at the existing pull request list. For +example, a pull request might be a good start, but changes are +requested by EXP maintainers, and the original pull request author +hasn't had time to address them. In this case consider picking up this +pull request: open a new pull request with all commits from the +original pull request, as well as additional changes to address the +raised issues. Doing so helps us a lot; it is not considered rude as +long as the original author is acknowledged by keeping his/her +commits. + +You can pull an existing pull request to a local branch by running +`git fetch https://github.com/EXP-code/code pull/$PR_NUMBER/head:$BRANCH_NAME_TO_CREATE` +(replace `$PR_NUMBER` with an ID of the pull request, and +`$BRANCH_NAME_TO_CREATE` with a name of the branch you want to create +locally). See also: +. + +When writing GitHub pull requests, try to keep titles short but +descriptive. Complete titles make it easy to skim through the issue +tracker. + +Finally, we all appreciate improving the code's readability and though +formatting and improved comments. But please, try to keep aesthetic +changes in separate commits from functional changes. This will make pull +requests easier to review and more likely to get merged. + +## Adding a Git Commit + +All code updates are done by pull request (PR) as described above. +The GitHub EXP-code includes Continuous Integration (CI) support which +can get very chatty and take up unnecessary compute resources. When +your changes only affect documentation (e.g. docstring additions or +software comments) or provide code snippets that you know still need +work, you may add a [no ci] in your commit message. For example: bash + +``` +> git commit -m "Updated some docstrings [no ci]" +``` + +When this commit is pushed out to your fork and branch associated with a +pull request, all CI will be skipped. + +# Documentation strategy + +EXP and pyEXP has three types of documentation: + +1. Inline Doxygen comments in C/C++ class headers. If you have not + used [Doxgyen](http://doxygen.org) in the past, we recommend that + you simply copy the style in existing headers in the `include` + directory in the EXP source tree. In short, Doxygen uses stylized + comments to extract documentation. Lines that start with `//!` or + blocks that start with `/**` and end with `*/` will end up + describing the immediately following class or member function. +2. Python wrapper code has embedded Python docstrings. For some + examples, check the C++ code in the `pyEXP`. +3. The ReadTheDocs manual that you are currently reading is designed to + provide an overview and tutorial for using EXP/pyEXP. You can fork + and issue pull requests against the `EXP-code/docs` repo just for + the EXP source code as described in + **Submitting-patches** above. + +## Contributing Documentation + +Our goal is a set of consistent documentation that provides users and +developers a shallow learning curve for using and adding to EXP. For end +users, we strive to write simple step-by-step instructions for common +tasks and give a clear description of the features and capabilities of +the EXP software. + +However, it *is* hard to know what everyone needs. As you work with this +package, we would love help with this and encourage all of your +contributions. + +This section is an attempt to provide some stylistic guidelines for +documentation writers and developers. For EXP and the documentation +overall, we hope that the existing documentation is a good starting +point. For internal Python documentation in pyEXP, we are trying to +follow the now familiar style of code documentation of both the Astropy +and Numpy projects. In particular, we have adopted the Numpy style +guidelines +[here](https://numpy.org/doc/devdocs/docs/howto_document.html). + +## Adding pyEXP docstrings + +pyEXP is a set of bindings implemented by +[pybind11](https://pybind11.readthedocs.io/) and a small amount of +native Python code. It has a full set of docstrings to provide users +with easy access to interactive tips and call signatures. If you would +like to contribute new code, please try to follow the following +guidelines: + +- Please write docstrings for all public classes, methods, and + functions +- Please consult the + [numpydoc](https://numpy.org/doc/devdocs/docs/howto_document.html) + format from the link above whenever possible +- We would like to encourage references to internal EXP project links + in docstrings when useful. This is still a work in progress. +- Examples and/or tutorials are strongly encouraged for typical + use-cases of a particular module or class. These can be included in + the [EXP examples + repository](https://github.com/EXP-code/EXP-examples) or the [pyEXP + examples repository](https://github.com/EXP-code/pyEXP-examples)) as + appropriate. + +# Writing examples + +We strongly encourage to contribute interesting examples and workflows +to one of the two example repositories: + +1. Use the [EXP examples + repo](https://github.com/EXP-code/EXP-examples) to illustrate + simulations. Check the existing ones for guidance. It's best if + your examples contain sample configuration files and phase-space + files, or possibly instructions for computing the phase-space files + if they are large. +2. Use the [pyEXP examples + repo](https://github.com/EXP-code/pyEXP-examples). Either documented + Python or Jupyter notebooks are ideal. + +# Running tests + +Please check any bug fixes and proposed new functionality works with +existing examples. Here are some suggested guidelines for EXP and pyEXP +changes, respectively: + +1. For EXP, clone the [EXP examples + repo](https://github.com/EXP-code/EXP-examples) to test changes to + the EXP simulation code. At the very least, please try the `Nbody` + example, but please try as many as you can to make sure that your + change will not break an existing use case for someone else. If your + change introduces a feature, please try to devise and contribute + example that demonstrates the new feature. That way, your new + feature will be tested by all future proposed changes and will help + others understand how to use your new feature. +2. The drill for pyEXP is similar. Clone [pyEXP examples + repo](https://github.com/EXP-code/pyEXP-examples) to test changes to + the pyEXP N-body analysis code. There are many work flows here, we + don't expect anyone to try all of them. But please use your best + judgment to try those affected by your proposed change. From c7948acb021a85743703631c2510019ada63aa2c Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 3 Jul 2024 21:55:40 -0700 Subject: [PATCH 08/40] Basis expects a node with the 'parameters' key for parsing [no ci] --- src/OutVel.cc | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/OutVel.cc b/src/OutVel.cc index 78225f9d2..a10a732ce 100644 --- a/src/OutVel.cc +++ b/src/OutVel.cc @@ -62,7 +62,10 @@ OutVel::OutVel(const YAML::Node& conf) : Output(conf) // Create the basis // - basis = std::make_shared(conf); + YAML::Node node; + node["parameters"] = conf; + + basis = std::make_shared(node); // Create the coefficient container based on the dimensionality. // Currently, these are spherical and polar grids. From b70207ff314ffe54110ef4cbbb075941ee7f1fdd Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 4 Jul 2024 13:15:58 -0700 Subject: [PATCH 09/40] Fixes for HDF5 libraries with CMake --- utils/ICs/CMakeLists.txt | 3 ++- utils/PhaseSpace/CMakeLists.txt | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/utils/ICs/CMakeLists.txt b/utils/ICs/CMakeLists.txt index b0ec76de6..06cbb8b18 100644 --- a/utils/ICs/CMakeLists.txt +++ b/utils/ICs/CMakeLists.txt @@ -5,7 +5,8 @@ set(bin_PROGRAMS shrinkics gensph gendisk gendisk2d gsphere pstmod cubeics zangics slabics) set(common_LINKLIB OpenMP::OpenMP_CXX MPI::MPI_CXX yaml-cpp exputil - ${VTK_LIBRARIES} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES}) + ${VTK_LIBRARIES} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES} + ${HDF5_CXX_LIBRARIES}) if(ENABLE_CUDA) list(APPEND common_LINKLIB CUDA::cudart CUDA::nvToolsExt) diff --git a/utils/PhaseSpace/CMakeLists.txt b/utils/PhaseSpace/CMakeLists.txt index 3076d0b12..71d1c2141 100644 --- a/utils/PhaseSpace/CMakeLists.txt +++ b/utils/PhaseSpace/CMakeLists.txt @@ -14,7 +14,7 @@ if(HAVE_XDR) endif() set(common_LINKLIB OpenMP::OpenMP_CXX MPI::MPI_CXX yaml-cpp expui - exputil ${VTK_LIBRARIES} ${HDF5_LIBRARIES}) + exputil ${VTK_LIBRARIES} ${HDF5_LIBRARIES} ${HDF5_CXX_LIBRARIES}) if(PNG_FOUND) list(APPEND common_LINKLIB PNG::PNG) From 77a1ba953f7adf2d193b521ff9a9f7f1ad841ed3 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 4 Jul 2024 14:50:28 -0700 Subject: [PATCH 10/40] Fixes for velocity field analysis --- expui/CoefContainer.H | 13 +++ expui/CoefContainer.cc | 182 ++++++++++++++++++++++++++++++++ expui/FieldBasis.cc | 8 +- utils/ICs/CMakeLists.txt | 3 +- utils/PhaseSpace/CMakeLists.txt | 2 +- 5 files changed, 202 insertions(+), 6 deletions(-) diff --git a/expui/CoefContainer.H b/expui/CoefContainer.H index 66a3c9951..0c8bf39d0 100644 --- a/expui/CoefContainer.H +++ b/expui/CoefContainer.H @@ -91,6 +91,7 @@ namespace MSSA void pack_sphere(); void unpack_sphere(); void restore_background_sphere(); + //@} //@{ //! Cylindrical coefficients @@ -99,6 +100,18 @@ namespace MSSA void restore_background_cylinder(); //@} + //@{ + //! Spherical field coefficients + void pack_sphfld(); + void unpack_sphfld(); + //@} + + //@{ + //! Cylindrical coefficients + void pack_cylfld(); + void unpack_cylfld(); + //@} + //@{ //! Slab coefficients void pack_slab(); diff --git a/expui/CoefContainer.cc b/expui/CoefContainer.cc index 03af71d92..a0ec98abe 100644 --- a/expui/CoefContainer.cc +++ b/expui/CoefContainer.cc @@ -3,6 +3,7 @@ #include #include #include +#include #include @@ -64,6 +65,12 @@ namespace MSSA else if (dynamic_cast(coefs.get())) { unpack_table(); } + else if (dynamic_cast(coefs.get())) { + unpack_sphfld(); + } + else if (dynamic_cast(coefs.get())) { + unpack_cylfld(); + } else { throw std::runtime_error("CoefDB::unpack_channels(): can not reflect coefficient type"); } @@ -79,6 +86,10 @@ namespace MSSA restore_background_cylinder(); else if (dynamic_cast(coefs.get())) { } // Do nothing + else if (dynamic_cast(coefs.get())) + { } // Do nothing + else if (dynamic_cast(coefs.get())) + { } // Do nothing else { throw std::runtime_error("CoefDB::background(): can not reflect coefficient type"); } @@ -86,6 +97,8 @@ namespace MSSA void CoefDB::pack_channels() { + std::cout << "Type: " << typeid(*coefs.get()).name() << std::endl; + if (dynamic_cast(coefs.get())) pack_sphere(); else if (dynamic_cast(coefs.get())) @@ -96,6 +109,10 @@ namespace MSSA pack_cube(); else if (dynamic_cast(coefs.get())) pack_table(); + else if (dynamic_cast(coefs.get())) + pack_sphfld(); + else if (dynamic_cast(coefs.get())) + pack_cylfld(); else { throw std::runtime_error("CoefDB::pack_channels(): can not reflect coefficient type"); } @@ -201,6 +218,85 @@ namespace MSSA // END time loop } + void CoefDB::pack_cylfld() + { + auto cur = dynamic_cast(coefs.get()); + + times = cur->Times(); + complexKey = true; + + auto cf = dynamic_cast( cur->getCoefStruct(times[0]).get() ); + + int nfld = cf->nfld; + int mmax = cf->mmax; + int nmax = cf->nmax; + int ntimes = times.size(); + + // Promote desired keys into c/s pairs + // + keys.clear(); + for (auto v : keys0) { + // Sanity check rank + // + if (v.size() != 3) { + std::ostringstream sout; + sout << "CoefDB::pack_cylfld: key vector should have rank 3; " + << "found rank " << v.size() << " instead"; + throw std::runtime_error(sout.str()); + } + // Sanity check values + // + else if (v[0] >= 0 and v[0] < nfld and + v[1] >= 0 and v[1] <= mmax and + v[2] >= 0 and v[2] <= nmax ) { + auto c = v, s = v; + c.push_back(0); + s.push_back(1); + keys.push_back(c); + if (v[0]) keys.push_back(s); + } else { + throw std::runtime_error("CoefDB::pack_cylfld: key is out of bounds"); + } + } + + bkeys.clear(); // No background fields + + // Only pack the keys in the list + // + for (auto k : keys) { + data[k].resize(ntimes); + } + + for (int t=0; t( cur->getCoefStruct(times[t]).get() ); + for (auto k : keys) { + if (k[3]==0) + data[k][t] = (*cf->coefs)(k[0], k[1], k[2]).real(); + else + data[k][t] = (*cf->coefs)(k[0], k[1], k[3]).imag(); + } + } + } + + void CoefDB::unpack_cylfld() + { + for (int i=0; i( coefs->getCoefStruct(times[i]).get() ); + + for (auto k : keys0) { + auto c = k, s = k; + c.push_back(0); s.push_back(1); + + int f = k[1], m = k[1], n = k[2]; + + if (m==0) (*cf->coefs)(f, m, n) = {data[c][i], 0.0}; + else (*cf->coefs)(f, m, n) = {data[c][i], data[s][i]}; + } + // END key loop + } + // END time loop + } + void CoefDB::pack_sphere() { auto cur = dynamic_cast(coefs.get()); @@ -309,6 +405,92 @@ namespace MSSA } // END time loop } + + void CoefDB::pack_sphfld() + { + auto cur = dynamic_cast(coefs.get()); + + times = cur->Times(); + complexKey = true; + + auto cf = dynamic_cast( cur->getCoefStruct(times[0]).get() ); + + int nfld = cf->nfld; + int lmax = cf->lmax; + int nmax = cf->nmax; + int ntimes = times.size(); + + // Make extended key list + // + keys.clear(); + for (auto k : keys0) { + // Sanity check rank + // + if (k.size() != 4) { + std::ostringstream sout; + sout << "CoefDB::pack_sphfld: key vector should have rank 4; " + << "found rank " << k.size() << " instead"; + throw std::runtime_error(sout.str()); + } + // Sanity check values + // + else if (k[0] < nfld and k[0] >= 0 and + k[1] <= lmax and k[1] >= 0 and + k[2] <= k[1] and k[2] >= 0 and + k[3] <= nmax and k[3] >= 0 ) { + + auto v = k; + v.push_back(0); + keys.push_back(v); + data[v].resize(ntimes); + + if (k[2]>0) { + v[4] = 1; + keys.push_back(v); + data[v].resize(ntimes); + } + } + else { + throw std::runtime_error("CoefDB::pack_sphfld: key is out of bounds"); + } + } + + bkeys.clear(); // No background for field quantities + + auto I = [](const Key& k) { return k[1]*(k[1]+1)/2 + k[2]; }; + + for (int t=0; t( cur->getCoefStruct(times[t]).get() ); + for (auto k : keys) { + auto c = (*cf->coefs)(k[0], I(k), k[3]); + data[k][t] = c.real(); + if (k[4]) data[k][t] = c.imag(); + } + } + } + + void CoefDB::unpack_sphfld() + { + auto I = [](const Key& k) { return k[1]*(k[1]+1)/2 + k[2]; }; + + for (int i=0; i( coefs->getCoefStruct(times[i]).get() ); + + for (auto k : keys0) { + auto c = k, s = k; + c.push_back(0); + s.push_back(1); + + int f = k[0], l = k[1], m = k[2], n = k[3]; + + if (m==0) (*cf->coefs)(f, I(k), n) = {data[c][i], 0.0 }; + else (*cf->coefs)(f, I(k), n) = {data[c][i], data[s][i]}; + } + // END key loop + } + // END time loop + } void CoefDB::pack_slab() { diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc index eecdcc520..ae0a46278 100644 --- a/expui/FieldBasis.cc +++ b/expui/FieldBasis.cc @@ -304,15 +304,15 @@ namespace BasisClasses (*coefs[tid])(0, 0, 0) += mass*p(0); - for (int l=0, k=0; l<=lmax; l++) { - for (int m=0; m<=l; m++, k++) { + for (int l=0, lm=0; l<=lmax; l++) { + for (int m=0; m<=l; m++, lm++) { std::complex P = std::exp(I*(phi*m))*Ylm01(l, m)*plgndr(l, m, cth); for (int n=0; n Date: Thu, 4 Jul 2024 16:15:39 -0700 Subject: [PATCH 11/40] Fix HighFive type mismatch that was causing an exception in writing the mSSA cache --- expui/expMSSA.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/expMSSA.cc b/expui/expMSSA.cc index 0aca81d71..5adf38f57 100644 --- a/expui/expMSSA.cc +++ b/expui/expMSSA.cc @@ -1339,7 +1339,7 @@ namespace MSSA { // Save trend state // int trend = static_cast::type>(type); - file.createAttribute("trendType", HighFive::DataSpace::From(trend)).write(trend); + file.createAttribute("trendType", HighFive::DataSpace::From(trend)).write(trend); // Save the key list // From 3f8e317cf8b75486b95c311903d2a277e5858d5a Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 4 Jul 2024 22:52:03 -0700 Subject: [PATCH 12/40] Fix H5 type for correct cache writing --- expui/expMSSA.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/expMSSA.cc b/expui/expMSSA.cc index 3ff1c71c7..7a19c6f6d 100644 --- a/expui/expMSSA.cc +++ b/expui/expMSSA.cc @@ -1339,7 +1339,7 @@ namespace MSSA { // Save trend state // int trend = static_cast::type>(type); - file.createAttribute("trendType", HighFive::DataSpace::From(trend)).write(trend); + file.createAttribute("trendType", HighFive::DataSpace::From(trend)).write(trend); // Save the key list // From aef82b5ac3539d0c7f1482d791a93bc4fae7ab34 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 4 Jul 2024 22:52:36 -0700 Subject: [PATCH 13/40] Implement proper coefficient parsing for velocity field basis --- expui/CoefContainer.cc | 85 ++++++++++++++++++++++++++++-------------- expui/FieldBasis.cc | 7 ++++ 2 files changed, 65 insertions(+), 27 deletions(-) diff --git a/expui/CoefContainer.cc b/expui/CoefContainer.cc index a0ec98abe..d5deb13da 100644 --- a/expui/CoefContainer.cc +++ b/expui/CoefContainer.cc @@ -72,7 +72,7 @@ namespace MSSA unpack_cylfld(); } else { - throw std::runtime_error("CoefDB::unpack_channels(): can not reflect coefficient type"); + throw std::runtime_error(std::string("CoefDB::unpack_channels(): can not reflect coefficient type=") + typeid(*coefs.get()).name()); } return coefs; @@ -91,14 +91,12 @@ namespace MSSA else if (dynamic_cast(coefs.get())) { } // Do nothing else { - throw std::runtime_error("CoefDB::background(): can not reflect coefficient type"); + throw std::runtime_error(std::string("CoefDB::background(): can not reflect coefficient type=") + typeid(*coefs.get()).name()); } } void CoefDB::pack_channels() { - std::cout << "Type: " << typeid(*coefs.get()).name() << std::endl; - if (dynamic_cast(coefs.get())) pack_sphere(); else if (dynamic_cast(coefs.get())) @@ -114,7 +112,7 @@ namespace MSSA else if (dynamic_cast(coefs.get())) pack_cylfld(); else { - throw std::runtime_error("CoefDB::pack_channels(): can not reflect coefficient type"); + throw std::runtime_error(std::string("CoefDB::pack_channels(): can not reflect coefficient type=") + typeid(*coefs.get()).name()); } } @@ -125,7 +123,8 @@ namespace MSSA times = cur->Times(); complexKey = true; - auto cf = dynamic_cast( cur->getCoefStruct(times[0]).get() ); + auto cf = dynamic_cast + ( cur->getCoefStruct(times[0]).get() ); int mmax = cf->mmax; int nmax = cf->nmax; @@ -182,7 +181,9 @@ namespace MSSA } for (int t=0; t( cur->getCoefStruct(times[t]).get() ); + cf = dynamic_cast + ( cur->getCoefStruct(times[t]).get() ); + for (auto k : keys) { if (k[2]==0) data[k][t] = (*cf->coefs)(k[0], k[1]).real(); @@ -202,7 +203,8 @@ namespace MSSA void CoefDB::unpack_cylinder() { for (int i=0; i( coefs->getCoefStruct(times[i]).get() ); + auto cf = dynamic_cast + ( coefs->getCoefStruct(times[i]).get() ); for (auto k : keys0) { auto c = k, s = k; @@ -225,7 +227,8 @@ namespace MSSA times = cur->Times(); complexKey = true; - auto cf = dynamic_cast( cur->getCoefStruct(times[0]).get() ); + auto cf = dynamic_cast + ( cur->getCoefStruct(times[0]).get() ); int nfld = cf->nfld; int mmax = cf->mmax; @@ -268,7 +271,9 @@ namespace MSSA } for (int t=0; t( cur->getCoefStruct(times[t]).get() ); + cf = dynamic_cast + ( cur->getCoefStruct(times[t]).get() ); + for (auto k : keys) { if (k[3]==0) data[k][t] = (*cf->coefs)(k[0], k[1], k[2]).real(); @@ -281,7 +286,8 @@ namespace MSSA void CoefDB::unpack_cylfld() { for (int i=0; i( coefs->getCoefStruct(times[i]).get() ); + auto cf = dynamic_cast + ( coefs->getCoefStruct(times[i]).get() ); for (auto k : keys0) { auto c = k, s = k; @@ -304,7 +310,8 @@ namespace MSSA times = cur->Times(); complexKey = true; - auto cf = dynamic_cast( cur->getCoefStruct(times[0]).get() ); + auto cf = dynamic_cast + ( cur->getCoefStruct(times[0]).get() ); int lmax = cf->lmax; int nmax = cf->nmax; @@ -368,7 +375,9 @@ namespace MSSA auto I = [](const Key& k) { return k[0]*(k[0]+1)/2 + k[1]; }; for (int t=0; t( cur->getCoefStruct(times[t]).get() ); + cf = dynamic_cast + ( cur->getCoefStruct(times[t]).get() ); + for (auto k : keys) { auto c = (*cf->coefs)(I(k), k[2]); data[k][t] = c.real(); @@ -389,7 +398,8 @@ namespace MSSA for (int i=0; i( coefs->getCoefStruct(times[i]).get() ); + auto cf = dynamic_cast + ( coefs->getCoefStruct(times[i]).get() ); for (auto k : keys0) { auto c = k, s = k; @@ -413,7 +423,8 @@ namespace MSSA times = cur->Times(); complexKey = true; - auto cf = dynamic_cast( cur->getCoefStruct(times[0]).get() ); + auto cf = dynamic_cast + ( cur->getCoefStruct(times[0]).get() ); int nfld = cf->nfld; int lmax = cf->lmax; @@ -437,7 +448,7 @@ namespace MSSA else if (k[0] < nfld and k[0] >= 0 and k[1] <= lmax and k[1] >= 0 and k[2] <= k[1] and k[2] >= 0 and - k[3] <= nmax and k[3] >= 0 ) { + k[3] < nmax and k[3] >= 0 ) { auto v = k; v.push_back(0); @@ -460,7 +471,10 @@ namespace MSSA auto I = [](const Key& k) { return k[1]*(k[1]+1)/2 + k[2]; }; for (int t=0; t( cur->getCoefStruct(times[t]).get() ); + + cf = dynamic_cast + ( cur->getCoefStruct(times[t]).get() ); + for (auto k : keys) { auto c = (*cf->coefs)(k[0], I(k), k[3]); data[k][t] = c.real(); @@ -475,10 +489,13 @@ namespace MSSA for (int i=0; i( coefs->getCoefStruct(times[i]).get() ); + auto cf = dynamic_cast + ( coefs->getCoefStruct(times[i]).get() ); for (auto k : keys0) { + auto c = k, s = k; + c.push_back(0); s.push_back(1); @@ -499,7 +516,8 @@ namespace MSSA times = cur->Times(); complexKey = true; - auto cf = dynamic_cast( cur->getCoefStruct(times[0]).get() ); + auto cf = dynamic_cast + ( cur->getCoefStruct(times[0]).get() ); int nmaxx = cf->nmaxx; int nmaxy = cf->nmaxy; @@ -558,7 +576,9 @@ namespace MSSA } for (int t=0; t( cur->getCoefStruct(times[t]).get() ); + cf = dynamic_cast + ( cur->getCoefStruct(times[t]).get() ); + for (auto k : keys) { auto c = (*cf->coefs)(k[0], k[1], k[2]); if (k[3]) data[k][t] = c.imag(); @@ -580,7 +600,8 @@ namespace MSSA times = cur->Times(); complexKey = true; - auto cf = dynamic_cast( cur->getCoefStruct(times[0]).get() ); + auto cf = dynamic_cast + ( cur->getCoefStruct(times[0]).get() ); int nmaxx = cf->nmaxx; int nmaxy = cf->nmaxy; @@ -639,7 +660,9 @@ namespace MSSA } for (int t=0; t( cur->getCoefStruct(times[t]).get() ); + cf = dynamic_cast + ( cur->getCoefStruct(times[t]).get() ); + for (auto k : keys) { auto c = (*cf->coefs)(k[0], k[1], k[2]); if (k[3]) data[k][t] = c.imag(); @@ -658,7 +681,8 @@ namespace MSSA { for (int i=0; i( coefs->getCoefStruct(times[i]).get() ); + auto cf = dynamic_cast + ( coefs->getCoefStruct(times[i]).get() ); for (auto k : keys0) { auto c = k, s = k; @@ -676,7 +700,8 @@ namespace MSSA { for (int i=0; i( coefs->getCoefStruct(times[i]).get() ); + auto cf = dynamic_cast + ( coefs->getCoefStruct(times[i]).get() ); for (auto k : keys0) { auto c = k, s = k; @@ -697,7 +722,8 @@ namespace MSSA times = cur->Times(); complexKey = false; - auto cf = dynamic_cast( cur->getCoefStruct(times[0]).get() ); + auto cf = dynamic_cast + ( cur->getCoefStruct(times[0]).get() ); int cols = cf->cols; int ntimes = times.size(); @@ -734,7 +760,9 @@ namespace MSSA for (unsigned c=0; c( cur->getCoefStruct(times[t]).get() ); + cf = dynamic_cast + ( cur->getCoefStruct(times[t]).get() ); + data[key][t] = (*cf->coefs)(c).real(); } } @@ -744,7 +772,8 @@ namespace MSSA { for (int i=0; i( coefs->getCoefStruct(times[i]).get() ); + auto cf = dynamic_cast + ( coefs->getCoefStruct(times[i]).get() ); int cols = cf->cols; @@ -862,6 +891,7 @@ namespace MSSA auto I = [](const Key& k) { return k[0]*(k[0]+1)/2 + k[1]; }; for (int t=0; t (cur->getCoefStruct(times[t]).get()); @@ -882,6 +912,7 @@ namespace MSSA auto cur = dynamic_cast(coefs.get()); for (int t=0; t (cur->getCoefStruct(times[t]).get()); diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc index ae0a46278..35a0089a5 100644 --- a/expui/FieldBasis.cc +++ b/expui/FieldBasis.cc @@ -290,9 +290,13 @@ namespace BasisClasses (*coefs[tid])(0, 0, 0) += mass*p(0); for (int m=0; m<=lmax; m++) { + std::complex P = std::exp(I*(phi*m)); + for (int n=0; n P = std::exp(I*(phi*m))*Ylm01(l, m)*plgndr(l, m, cth); for (int n=0; n Date: Fri, 5 Jul 2024 14:54:55 -0700 Subject: [PATCH 14/40] Parameter name changes for consistency with main N-body code --- expui/BasisFactory.cc | 2 + expui/FieldBasis.H | 8 ++-- expui/FieldBasis.cc | 103 ++++++++++++++++++++++++++++-------------- src/OutVel.H | 6 +-- src/OutVel.cc | 8 ++-- 5 files changed, 82 insertions(+), 45 deletions(-) diff --git a/expui/BasisFactory.cc b/expui/BasisFactory.cc index 5b574a8e2..ee6d53d4d 100644 --- a/expui/BasisFactory.cc +++ b/expui/BasisFactory.cc @@ -84,6 +84,8 @@ namespace BasisClasses // try { conf = node["parameters"]; + YAML::Emitter y; y << conf; + std::cout << "YAML in Basis: " << y.c_str() << std::endl; } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing Basis parameters for <" diff --git a/expui/FieldBasis.H b/expui/FieldBasis.H index d405c3179..25ef5dbe3 100644 --- a/expui/FieldBasis.H +++ b/expui/FieldBasis.H @@ -53,9 +53,9 @@ namespace BasisClasses //@{ //! Parameters - std::string model, filename; + std::string model, modelname; int lmax, nmax; - double rmin, rmax, ascl, scale, delta; + double rmin, rmax, ascl, rmapping, delta; //@} //@{ @@ -116,8 +116,8 @@ namespace BasisClasses //! Register phase-space functionoid void addPSFunction(PSFunction func, std::vector& labels); - //! Prescaling factor - void set_scale(const double scl) { scale = scl; } + //! Coordinate mapping factor + void set_scale(const double scl) { rmapping = scl; } //! Zero out coefficients to prepare for a new expansion void reset_coefs(void); diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc index 35a0089a5..75c6d0a23 100644 --- a/expui/FieldBasis.cc +++ b/expui/FieldBasis.cc @@ -11,16 +11,23 @@ #include #endif -extern double Ylm01(int ll, int mm); extern double plgndr(int l, int m, double x); +static double Ylm(int l, int m) +{ + m = abs(m); + return sqrt( (2.0*l+1)/(4.0*M_PI) ) * + exp(0.5*(lgamma(1.0+l-m) - lgamma(1.0+l+m))); +} + + namespace BasisClasses { const std::set FieldBasis::valid_keys = { - "filename", + "modelname", "dof", - "scale", + "rmapping", "rmin", "rmax", "ascl", @@ -60,18 +67,18 @@ namespace BasisClasses // void FieldBasis::configure() { - nfld = 2; // Weight and density fields by default - lmax = 4; - nmax = 10; - rmin = 1.0e-4; - rmax = 2.0; - ascl = 0.01; - delta = 0.005; - scale = 0.05; - dof = 3; - model = "file"; - name = "field"; - filename = "SLGridSph.model"; + nfld = 2; // Weight and density fields by default + lmax = 4; + nmax = 10; + rmin = 1.0e-4; + rmax = 2.0; + ascl = 0.01; + delta = 0.005; + rmapping = 0.05; + dof = 3; + model = "file"; + name = "field"; + modelname = "SLGridSph.model"; initialize(); @@ -111,8 +118,8 @@ namespace BasisClasses // if (model == "file") { std::vector r, d; - std::ifstream in(filename); - if (not in) throw std::runtime_error("Error opening file: " + filename); + std::ifstream in(modelname); + if (not in) throw std::runtime_error("Error opening file: " + modelname); std::string line; while (std::getline(in, line)) { @@ -156,13 +163,28 @@ namespace BasisClasses // Generate the orthogonal function instance // - ortho = std::make_shared(nmax, densfunc, rmin, rmax, scale, dof); + ortho = std::make_shared(nmax, densfunc, rmin, rmax, rmapping, dof); // Initialize fieldlabels // fieldLabels.clear(); fieldLabels.push_back("weight"); fieldLabels.push_back("density"); + + // Debug + // + if (true) { + auto tst = ortho->testOrtho(); + double worst = 0.0; + for (int i=0; i(worst, fabs(1.0 - tst(i, j))); + else worst = std::max(worst, fabs(tst(i, j))); + } + } + if (myid==0) + std::cout << "FieldBasis: ortho test worst <" << worst << ">" << std::endl; + } } void FieldBasis::allocateStore() @@ -187,17 +209,17 @@ namespace BasisClasses // Assign values from YAML // try { - if (conf["filename"]) filename = conf["filename"].as(); - if (conf["nfld" ]) nfld = conf["nfld" ].as(); - if (conf["lmax" ]) lmax = conf["lmax" ].as(); - if (conf["nmax" ]) nmax = conf["nmax" ].as(); - if (conf["dof" ]) dof = conf["dof" ].as(); - if (conf["rmin" ]) rmin = conf["rmin" ].as(); - if (conf["rmax" ]) rmax = conf["rmax" ].as(); - if (conf["ascl" ]) ascl = conf["ascl" ].as(); - if (conf["delta" ]) delta = conf["delta" ].as(); - if (conf["scale" ]) scale = conf["scale" ].as(); - if (conf["model" ]) model = conf["model" ].as(); + if (conf["modelname"]) modelname = conf["modelname"].as(); + if (conf["model" ]) model = conf["model" ].as(); + if (conf["nfld" ]) nfld = conf["nfld" ].as(); + if (conf["lmax" ]) lmax = conf["lmax" ].as(); + if (conf["nmax" ]) nmax = conf["nmax" ].as(); + if (conf["dof" ]) dof = conf["dof" ].as(); + if (conf["rmin" ]) rmin = conf["rmin" ].as(); + if (conf["rmax" ]) rmax = conf["rmax" ].as(); + if (conf["ascl" ]) ascl = conf["ascl" ].as(); + if (conf["delta" ]) delta = conf["delta" ].as(); + if (conf["rmapping" ]) rmapping = conf["rmapping" ].as(); } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing parameters in FieldBasis: " @@ -291,7 +313,7 @@ namespace BasisClasses for (int m=0; m<=lmax; m++) { - std::complex P = std::exp(I*(phi*m)); + std::complex P = std::exp(-I*(phi*m)); for (int n=0; n P = - std::exp(I*(phi*m))*Ylm01(l, m)*plgndr(l, m, cth); - + std::exp(-I*(phi*m))*Ylm(l, m)*plgndr(l, m, cth) * s; + + s *= -1.0; // Flip sign for next evaluation + for (int n=0; n0 ? 2.0 : 1.0; for (int n=0; n0 ? 2.0 : 1.0; for (int n=0; n::max(); Component *tcomp; CoefClasses::CoefsPtr coefs; diff --git a/src/OutVel.cc b/src/OutVel.cc index a10a732ce..cd877e4d3 100644 --- a/src/OutVel.cc +++ b/src/OutVel.cc @@ -9,12 +9,12 @@ const std::set OutVel::valid_keys = { - "filename", + "modelname", "nint", "nintsub", "name", "dof", - "scale", + "rmapping", "rmin", "rmax", "ascl", @@ -97,7 +97,7 @@ void OutVel::initialize() model = conf["model"].as(); else { std::string message = "OutVel: no model specified. Please specify " - "either 'file' with the model 'filename' or 'expon' for the\n" + "either 'file' with the model 'modelname' or 'expon' for the\n" "exponential disk model (i.e. Laguerre polynomials)"; throw std::runtime_error(message); } @@ -119,7 +119,7 @@ void OutVel::initialize() throw std::runtime_error(message); } - if (conf["filename"]) filename = conf["filename"].as(); + if (conf["modelname"]) modelname = conf["modelname"].as(); } catch (YAML::Exception & error) { From eb998ddac5ff7727f306b6e0dc56959cfeb606bd Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 5 Jul 2024 15:25:55 -0700 Subject: [PATCH 15/40] Removed some debugging cruft --- expui/BasisFactory.cc | 2 -- expui/FieldBasis.cc | 4 ---- 2 files changed, 6 deletions(-) diff --git a/expui/BasisFactory.cc b/expui/BasisFactory.cc index ee6d53d4d..5b574a8e2 100644 --- a/expui/BasisFactory.cc +++ b/expui/BasisFactory.cc @@ -84,8 +84,6 @@ namespace BasisClasses // try { conf = node["parameters"]; - YAML::Emitter y; y << conf; - std::cout << "YAML in Basis: " << y.c_str() << std::endl; } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing Basis parameters for <" diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc index 75c6d0a23..29adb26d9 100644 --- a/expui/FieldBasis.cc +++ b/expui/FieldBasis.cc @@ -755,10 +755,6 @@ namespace BasisClasses VelocityBasis::VelocityBasis(const YAML::Node& conf) : FieldBasis(conf, "VelocityBasis") { - YAML::Emitter y; y << conf; - std::cout << "YAML configuration in VelocityBasis: " - << y.c_str() << std::endl; - name = "velocity"; assignFunc(); } From 971f33d764a1be110e9fa460c51d3580978973a0 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 6 Jul 2024 11:19:33 -0700 Subject: [PATCH 16/40] Pre-release version bump --- CMakeLists.txt | 2 +- doc/exp.cfg | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 598bf4ad9..cd65e7a1a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.21) # Needed for CUDA, MPI, and CTest features project( EXP - VERSION "7.7.30" + VERSION "7.7.99" HOMEPAGE_URL https://github.com/EXP-code/EXP LANGUAGES C CXX Fortran) diff --git a/doc/exp.cfg b/doc/exp.cfg index b7ccb4716..32f7444ff 100644 --- a/doc/exp.cfg +++ b/doc/exp.cfg @@ -48,7 +48,7 @@ PROJECT_NAME = EXP # could be handy for archiving the generated documentation or if some version # control system is used. -PROJECT_NUMBER = 7.7.30 +PROJECT_NUMBER = 7.7.99 # Using the PROJECT_BRIEF tag one can provide an optional one line description # for a project that appears at the top of each page and should give viewer a From 1ccff70d7f980afd5f5a1daa3b666a9f8575a963 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 6 Jul 2024 15:33:29 -0700 Subject: [PATCH 17/40] Added a routine to dump the basis [no ci] --- exputil/OrthoFunction.cc | 27 +++++++++++++++++++++++++++ include/OrthoFunction.H | 3 +++ 2 files changed, 30 insertions(+) diff --git a/exputil/OrthoFunction.cc b/exputil/OrthoFunction.cc index 7a1badc20..e0bf4f1f9 100644 --- a/exputil/OrthoFunction.cc +++ b/exputil/OrthoFunction.cc @@ -1,3 +1,4 @@ +#include #include OrthoFunction::OrthoFunction @@ -93,6 +94,32 @@ Eigen::MatrixXd OrthoFunction::testOrtho() return ret; } +void OrthoFunction::dumpOrtho(const std::string& filename) +{ + std::ofstream fout(filename); + + if (fout) { + fout << "# OrthoFunction dump" << std::endl; + + const int number = 1000; + + for (int i=0; i Date: Sat, 6 Jul 2024 15:33:55 -0700 Subject: [PATCH 18/40] Added a factor to make the weight function normalized [no ci] --- expui/FieldBasis.cc | 48 ++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 45 insertions(+), 3 deletions(-) diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc index 29adb26d9..32533400d 100644 --- a/expui/FieldBasis.cc +++ b/expui/FieldBasis.cc @@ -182,8 +182,10 @@ namespace BasisClasses else worst = std::max(worst, fabs(tst(i, j))); } } - if (myid==0) - std::cout << "FieldBasis: ortho test worst <" << worst << ">" << std::endl; + if (myid==0) { + std::cout << "FieldBasis::orthoTest: worst=" << worst << std::endl; + ortho->dumpOrtho("fieldbasis_ortho.dump"); + } } } @@ -276,9 +278,33 @@ namespace BasisClasses if (dof==2) { auto p = dynamic_cast(c.get()); coefs[0] = p->coefs; + store[0] = p->store; + + // Sanity test dimensions + if (nfld!=p->nfld || lmax!=p->mmax || nmax!=p->nmax) { + std::ostringstream serr; + serr << "FieldBasis::set_coefs: dimension error! " + << " nfld [" << nfld << "!= " << p->nfld << "]" + << " mmax [" << lmax << "!= " << p->mmax << "]" + << " nmax [" << nmax << "!= " << p->nmax << "]"; + throw std::runtime_error(serr.str()); + } + } else { auto p = dynamic_cast(c.get()); coefs[0] = p->coefs; + store[0] = p->store; + + // Sanity test dimensions + if (nfld!=p->nfld || lmax!=p->lmax || nmax!=p->nmax) { + std::ostringstream serr; + serr << "FieldBasis::set_coefs: dimension error! " + << " nfld [" << nfld << "!= " << p->nfld << "]" + << " lmax [" << lmax << "!= " << p->lmax << "]" + << " nmax [" << nmax << "!= " << p->nmax << "]"; + throw std::runtime_error(serr.str()); + } + } } @@ -287,6 +313,8 @@ namespace BasisClasses double u, double v, double w) { constexpr std::complex I(0, 1); + constexpr double fac0 = 1.0/sqrt(4*M_PI); + int tid = omp_get_thread_num(); PS3 pos{x, y, z}, vel{u, v, w}; @@ -309,7 +337,7 @@ namespace BasisClasses auto p = (*ortho)(R); - (*coefs[tid])(0, 0, 0) += mass*p(0); + (*coefs[tid])(0, 0, 0) += mass*p(0)*fac0; for (int m=0; m<=lmax; m++) { @@ -431,6 +459,20 @@ namespace BasisClasses } return ret; } else { + static bool firstime = true; + if (firstime) { + int tid = omp_get_thread_num(); + std::ostringstream file; + file << "field.bin." << tid; + std::ofstream fout(file.str()); + const auto& d = coefs[0]->dimensions(); + fout << "Dim size: " << d.size(); + for (int i=0; i ret(nfld, 0); auto p = (*ortho)(r); for (int i=0; i Date: Sun, 7 Jul 2024 10:57:28 -0700 Subject: [PATCH 19/40] Removed experimental extensions from previous mixed bug-fix and feature merge --- expui/CMakeLists.txt | 2 +- expui/SGSmooth.H | 13 - expui/SGSmooth.cc | 549 ------------------------------------------ expui/expMSSA.cc | 259 -------------------- pyEXP/MSSAWrappers.cc | 43 ---- 5 files changed, 1 insertion(+), 865 deletions(-) delete mode 100644 expui/SGSmooth.H delete mode 100644 expui/SGSmooth.cc diff --git a/expui/CMakeLists.txt b/expui/CMakeLists.txt index d506bea8d..8690817f7 100644 --- a/expui/CMakeLists.txt +++ b/expui/CMakeLists.txt @@ -34,7 +34,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 SGSmooth.cc) + Koopman.cc BiorthBess.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/SGSmooth.H b/expui/SGSmooth.H deleted file mode 100644 index 6615bffa2..000000000 --- a/expui/SGSmooth.H +++ /dev/null @@ -1,13 +0,0 @@ -#ifndef __SGSMOOTH_HPP__ -#define __SGSMOOTH_HPP__ - -#include - -//! savitzky golay smoothing. -std::vector sg_smooth(const std::vector &v, const int w, const int deg); - -//! numerical derivative based on savitzky golay smoothing. -std::vector sg_derivative(const std::vector &v, const int w, - const int deg, const double h=1.0); - -#endif // __SGSMOOTH_HPP__ diff --git a/expui/SGSmooth.cc b/expui/SGSmooth.cc deleted file mode 100644 index 1b5150b53..000000000 --- a/expui/SGSmooth.cc +++ /dev/null @@ -1,549 +0,0 @@ -//! -// Sliding window signal processing (and linear algebra toolkit). -// -// supported operations: -//
    -//
  • Savitzky-Golay smoothing. -//
  • computing a numerical derivative based of Savitzky-Golay smoothing. -//
  • required linear algebra support for SG smoothing using STL based -// vector/matrix classes -//
-// -// \brief Linear Algebra "Toolkit". -// -// modified by Rob Patro, 2016 -// modified by MDW, 2024 - -// system headers -#include -#include -#include // for size_t -#include // for fabs -#include - -//! default convergence -static const double TINY_FLOAT = 1.0e-300; - -//! comfortable array of doubles -using float_vect = std::vector; -//! comfortable array of ints; -using int_vect = std::vector; - -/*! matrix class. - * - * This is a matrix class derived from a vector of float_vects. Note that - * the matrix elements indexed [row][column] with indices starting at 0 (c - * style). Also note that because of its design looping through rows should - * be faster than looping through columns. - * - * \brief two dimensional floating point array - */ -class float_mat : public std::vector { -private: - //! disable the default constructor - explicit float_mat() {}; - //! disable assignment operator until it is implemented. - float_mat &operator =(const float_mat &) { return *this; }; -public: - //! constructor with sizes - float_mat(const size_t rows, const size_t cols, const double def=0.0); - //! copy constructor for matrix - float_mat(const float_mat &m); - //! copy constructor for vector - float_mat(const float_vect &v); - - //! use default destructor - // ~float_mat() {}; - - //! get size - size_t nr_rows(void) const { return size(); }; - //! get size - size_t nr_cols(void) const { return front().size(); }; -}; - - - -// constructor with sizes -float_mat::float_mat(const size_t rows,const size_t cols,const double defval) - : std::vector(rows) { - int i; - for (i = 0; i < rows; ++i) { - (*this)[i].resize(cols, defval); - } - if ((rows < 1) || (cols < 1)) { - std::ostringstream msg; - msg << "cannot build matrix with " << rows << " rows and " << cols - << "columns"; - throw std::runtime_error(msg.str()); - } -} - -// copy constructor for matrix -float_mat::float_mat(const float_mat &m) : std::vector(m.size()) { - - float_mat::iterator inew = begin(); - float_mat::const_iterator iold = m.begin(); - for (/* empty */; iold < m.end(); ++inew, ++iold) { - const size_t oldsz = iold->size(); - inew->resize(oldsz); - const float_vect oldvec(*iold); - *inew = oldvec; - } -} - -// copy constructor for vector -float_mat::float_mat(const float_vect &v) - : std::vector(1) { - - const size_t oldsz = v.size(); - front().resize(oldsz); - front() = v; -} - -////////////////////// -// Helper functions // -////////////////////// - -//! permute() orders the rows of A to match the integers in the index array. -void permute(float_mat &A, int_vect &idx) -{ - int_vect i(idx.size()); - int j,k; - - for (j = 0; j < A.nr_rows(); ++j) { - i[j] = j; - } - - // loop over permuted indices - for (j = 0; j < A.nr_rows(); ++j) { - if (i[j] != idx[j]) { - - // search only the remaining indices - for (k = j+1; k < A.nr_rows(); ++k) { - if (i[k] ==idx[j]) { - std::swap(A[j],A[k]); // swap the rows and - i[k] = i[j]; // the elements of - i[j] = idx[j]; // the ordered index. - break; // next j - } - } - } - } -} - -/*! \brief Implicit partial pivoting. - * - * The function looks for pivot element only in rows below the current - * element, A[idx[row]][column], then swaps that row with the current one in - * the index map. The algorithm is for implicit pivoting (i.e., the pivot is - * chosen as if the max coefficient in each row is set to 1) based on the - * scaling information in the vector scale. The map of swapped indices is - * recorded in swp. The return value is +1 or -1 depending on whether the - * number of row swaps was even or odd respectively. */ -static int partial_pivot(float_mat &A, const size_t row, const size_t col, - float_vect &scale, int_vect &idx, double tol) -{ - if (tol <= 0.0) - tol = TINY_FLOAT; - - int swapNum = 1; - - // default pivot is the current position, [row,col] - int pivot = row; - double piv_elem = fabs(A[idx[row]][col]) * scale[idx[row]]; - - // loop over possible pivots below current - int j; - for (j = row + 1; j < A.nr_rows(); ++j) { - - const double tmp = fabs(A[idx[j]][col]) * scale[idx[j]]; - - // if this elem is larger, then it becomes the pivot - if (tmp > piv_elem) { - pivot = j; - piv_elem = tmp; - } - } - -#if 0 - if (piv_elem < tol) { - throw std::runtime_error("partial_pivot(): Zero pivot encountered."); -#endif - - if(pivot > row) { // bring the pivot to the diagonal - j = idx[row]; // reorder swap array - idx[row] = idx[pivot]; - idx[pivot] = j; - swapNum = -swapNum; // keeping track of odd or even swap - } - return swapNum; -} - -/*! \brief Perform backward substitution. - * - * Solves the system of equations A*b=a, ASSUMING that A is upper - * triangular. If diag==1, then the diagonal elements are additionally - * assumed to be 1. Note that the lower triangular elements are never - * checked, so this function is valid to use after a LU-decomposition in - * place. A is not modified, and the solution, b, is returned in a. */ -static void lu_backsubst(float_mat &A, float_mat &a, bool diag=false) -{ - int r,c,k; - - for (r = (A.nr_rows() - 1); r >= 0; --r) { - for (c = (A.nr_cols() - 1); c > r; --c) { - for (k = 0; k < A.nr_cols(); ++k) { - a[r][k] -= A[r][c] * a[c][k]; - } - } - if(!diag) { - for (k = 0; k < A.nr_cols(); ++k) { - a[r][k] /= A[r][r]; - } - } - } -} - -/*! \brief Perform forward substitution. - * - * Solves the system of equations A*b=a, ASSUMING that A is lower - * triangular. If diag==1, then the diagonal elements are additionally - * assumed to be 1. Note that the upper triangular elements are never - * checked, so this function is valid to use after a LU-decomposition in - * place. A is not modified, and the solution, b, is returned in a. */ -static void lu_forwsubst(float_mat &A, float_mat &a, bool diag=true) -{ - int r,k,c; - for (r = 0;r < A.nr_rows(); ++r) { - for(c = 0; c < r; ++c) { - for (k = 0; k < A.nr_cols(); ++k) { - a[r][k] -= A[r][c] * a[c][k]; - } - } - if(!diag) { - for (k = 0; k < A.nr_cols(); ++k) { - a[r][k] /= A[r][r]; - } - } - } -} - -/*! \brief Performs LU factorization in place. - * - * This is Crout's algorithm (cf., Num. Rec. in C, Section 2.3). The map of - * swapped indeces is recorded in idx. The return value is +1 or -1 - * depending on whether the number of row swaps was even or odd - * respectively. idx must be preinitialized to a valid set of indices - * (e.g., {1,2, ... ,A.nr_rows()}). */ -static int lu_factorize(float_mat &A, int_vect &idx, double tol=TINY_FLOAT) -{ - if ( tol <= 0.0) - tol = TINY_FLOAT; - - if ((A.nr_rows() == 0) || (A.nr_rows() != A.nr_cols())) { - throw std::runtime_error("lu_factorize(): cannot handle empty " - "or nonsquare matrices."); - } - - float_vect scale(A.nr_rows()); // implicit pivot scaling - int i,j; - for (i = 0; i < A.nr_rows(); ++i) { - double maxval = 0.0; - for (j = 0; j < A.nr_cols(); ++j) { - if (fabs(A[i][j]) > maxval) - maxval = fabs(A[i][j]); - } - if (maxval == 0.0) { - throw std::runtime_error("lu_factorize(): zero pivot found."); - } - scale[i] = 1.0 / maxval; - } - - int swapNum = 1; - int c,r; - for (c = 0; c < A.nr_cols() ; ++c) { // loop over columns - swapNum *= partial_pivot(A, c, c, scale, idx, tol); // bring pivot to diagonal - for(r = 0; r < A.nr_rows(); ++r) { // loop over rows - int lim = (r < c) ? r : c; - for (j = 0; j < lim; ++j) { - A[idx[r]][c] -= A[idx[r]][j] * A[idx[j]][c]; - } - if (r > c) - A[idx[r]][c] /= A[idx[c]][c]; - } - } - permute(A,idx); - return swapNum; -} - -/*! \brief Solve a system of linear equations. - * Solves the inhomogeneous matrix problem with lu-decomposition. Note that - * inversion may be accomplished by setting a to the identity_matrix. */ -static float_mat lin_solve(const float_mat &A, const float_mat &a, - double tol=TINY_FLOAT) -{ - float_mat B(A); - float_mat b(a); - int_vect idx(B.nr_rows()); - int j; - - for (j = 0; j < B.nr_rows(); ++j) { - idx[j] = j; // init row swap label array - } - lu_factorize(B,idx,tol); // get the lu-decomp. - permute(b,idx); // sort the inhomogeneity to match the lu-decomp - lu_forwsubst(B,b); // solve the forward problem - lu_backsubst(B,b); // solve the backward problem - return b; -} - -/////////////////////// -// related functions // -/////////////////////// - -//! Returns the inverse of a matrix using LU-decomposition. -static float_mat invert(const float_mat &A) -{ - const int n = A.size(); - float_mat E(n, n, 0.0); - float_mat B(A); - int i; - - for (i = 0; i < n; ++i) { - E[i][i] = 1.0; - } - - return lin_solve(B, E); -} - -//! returns the transposed matrix. -static float_mat transpose(const float_mat &a) -{ - float_mat res(a.nr_cols(), a.nr_rows()); - int i,j; - - for (i = 0; i < a.nr_rows(); ++i) { - for (j = 0; j < a.nr_cols(); ++j) { - res[j][i] = a[i][j]; - } - } - return res; -} - -//! matrix multiplication. -float_mat operator *(const float_mat &a, const float_mat &b) -{ - float_mat res(a.nr_rows(), b.nr_cols()); - if (a.nr_cols() != b.nr_rows()) { - throw std::runtime_error("incompatible matrices in multiplication"); - } - - int i,j,k; - - for (i = 0; i < a.nr_rows(); ++i) { - for (j = 0; j < b.nr_cols(); ++j) { - double sum(0.0); - for (k = 0; k < a.nr_cols(); ++k) { - sum += a[i][k] * b[k][j]; - } - res[i][j] = sum; - } - } - return res; -} - - -//! calculate savitzky golay coefficients. -static float_vect sg_coeff(const float_vect &b, const size_t deg) -{ - const size_t rows(b.size()); - const size_t cols(deg + 1); - float_mat A(rows, cols); - float_vect res(rows); - - // generate input matrix for least squares fit - int i,j; - for (i = 0; i < rows; ++i) { - for (j = 0; j < cols; ++j) { - A[i][j] = pow(double(i), double(j)); - } - } - - float_mat c(invert(transpose(A) * A) * (transpose(A) * transpose(b))); - - for (i = 0; i < b.size(); ++i) { - res[i] = c[0][0]; - for (j = 1; j <= deg; ++j) { - res[i] += c[j][0] * pow(double(i), double(j)); - } - } - return res; -} - -/*! \brief savitzky golay smoothing. - * - * This method means fitting a polynome of degree 'deg' to a sliding window - * of width 2w+1 throughout the data. The needed coefficients are - * generated dynamically by doing a least squares fit on a "symmetric" unit - * vector of size 2w+1, e.g. for w=2 b=(0,0,1,0,0). evaluating the polynome - * yields the sg-coefficients. at the border non symmectric vectors b are - * used. */ -float_vect sg_smooth(const float_vect &v, const int width, const int deg) -{ - float_vect res(v.size(), 0.0); - if ((width < 1) || (deg < 0) || (v.size() < (2 * width + 2))) { - throw std::runtime_error("sgsmooth: parameter error."); - } - - const int window = 2 * width + 1; - const int endidx = v.size() - 1; - - // do a regular sliding window average - int i,j; - if (deg == 0) { - // handle border cases first because we need different coefficients -#if defined(_OPENMP) -#pragma omp parallel for private(i,j) schedule(static) -#endif - for (i = 0; i < width; ++i) { - const double scale = 1.0/double(i+1); - const float_vect c1(width, scale); - for (j = 0; j <= i; ++j) { - res[i] += c1[j] * v[j]; - res[endidx - i] += c1[j] * v[endidx - j]; - } - } - - // now loop over rest of data. reusing the "symmetric" coefficients. - const double scale = 1.0/double(window); - const float_vect c2(window, scale); -#if defined(_OPENMP) -#pragma omp parallel for private(i,j) schedule(static) -#endif - for (i = 0; i <= (v.size() - window); ++i) { - for (j = 0; j < window; ++j) { - res[i + width] += c2[j] * v[i + j]; - } - } - return res; - } - - // handle border cases first because we need different coefficients -#if defined(_OPENMP) -#pragma omp parallel for private(i,j) schedule(static) -#endif - for (i = 0; i < width; ++i) { - float_vect b1(window, 0.0); - b1[i] = 1.0; - - const float_vect c1(sg_coeff(b1, deg)); - for (j = 0; j < window; ++j) { - res[i] += c1[j] * v[j]; - res[endidx - i] += c1[j] * v[endidx - j]; - } - } - - // now loop over rest of data. reusing the "symmetric" coefficients. - float_vect b2(window, 0.0); - b2[width] = 1.0; - const float_vect c2(sg_coeff(b2, deg)); - -#if defined(_OPENMP) -#pragma omp parallel for private(i,j) schedule(static) -#endif - for (i = 0; i <= (v.size() - window); ++i) { - for (j = 0; j < window; ++j) { - res[i + width] += c2[j] * v[i + j]; - } - } - return res; -} - -/*! least squares fit a polynome of degree 'deg' to data in 'b'. - * then calculate the first derivative and return it. */ -static float_vect lsqr_fprime(const float_vect &b, const int deg) -{ - const int rows(b.size()); - const int cols(deg + 1); - float_mat A(rows, cols); - float_vect res(rows); - - // generate input matrix for least squares fit - int i,j; - for (i = 0; i < rows; ++i) { - for (j = 0; j < cols; ++j) { - A[i][j] = pow(double(i), double(j)); - } - } - - float_mat c(invert(transpose(A) * A) * (transpose(A) * transpose(b))); - - for (i = 0; i < b.size(); ++i) { - res[i] = c[1][0]; - for (j = 1; j < deg; ++j) { - res[i] += c[j + 1][0] * double(j+1) - * pow(double(i), double(j)); - } - } - return res; -} - -/*! \brief savitzky golay smoothed numerical derivative. - * - * This method means fitting a polynome of degree 'deg' to a sliding window - * of width 2w+1 throughout the data. - * - * In contrast to the sg_smooth function we do a brute force attempt by - * always fitting the data to a polynome of degree 'deg' and using the - * result. */ -float_vect sg_derivative(const float_vect &v, const int width, - const int deg, const double h) -{ - float_vect res(v.size(), 0.0); - if ((width < 1) || (deg < 1) || (v.size() < (2 * width + 2))) { - throw std::runtime_error("sgsderiv: parameter error"); - } - - const int window = 2 * width + 1; - - // handle border cases first because we do not repeat the fit - // lower part - float_vect b(window, 0.0); - int i,j; - - for (i = 0; i < window; ++i) { - b[i] = v[i] / h; - } - const float_vect c(lsqr_fprime(b, deg)); - for (j = 0; j <= width; ++j) { - res[j] = c[j]; - } - // upper part. direction of fit is reversed - for (i = 0; i < window; ++i) { - b[i] = v[v.size() - 1 - i] / h; - } - const float_vect d(lsqr_fprime(b, deg)); - for (i = 0; i <= width; ++i) { - res[v.size() - 1 - i] = -d[i]; - } - - // now loop over rest of data. wasting a lot of least squares calcs - // since we only use the middle value. -#if defined(_OPENMP) -#pragma omp parallel for private(i,j) schedule(static) -#endif - for (i = 1; i < (v.size() - window); ++i) { - for (j = 0; j < window; ++j) { - b[j] = v[i + j] / h; - } - res[i + width] = lsqr_fprime(b, deg)[width]; - } - return res; -} - -// Local Variables: -// mode: c++ -// c-basic-offset: 4 -// fill-column: 76 -// indent-tabs-mode: nil -// End: diff --git a/expui/expMSSA.cc b/expui/expMSSA.cc index 7a19c6f6d..b136d2e7d 100644 --- a/expui/expMSSA.cc +++ b/expui/expMSSA.cc @@ -1554,265 +1554,6 @@ namespace MSSA { } - std::tuple - expMSSA::getKoopmanModes(double tol, int D, bool debug) - { - bool use_fullKh = true; // Use the non-reduced computation of - // Koopman/eDMD - // Number of channels - // - nkeys = mean.size(); - - // Make sure parameters are sane - // - if (numW<=0) numW = numT/2; - if (numW > numT/2) numW = numT/2; - - numK = numT - numW + 1; - - Eigen::VectorXd S1; - Eigen::MatrixXd Y1; - Eigen::MatrixXd V1; - Eigen::MatrixXd VT1; - Eigen::MatrixXd VT2; - - // Make a new trajetory matrix with smoothing - // - Y1.resize(numK, numW*nkeys + D*(nkeys-1)); - - size_t n=0, offset=0; - for (auto k : mean) { - for (int i=0; i 0) { - // Back blending - for (int j=0; j(D-j)/D; - } - } - // Main series - for (int j=0; j(D-j)/D; - } - } - } - offset += numW + D; - n++; - } - - double Scale = Y1.norm(); - - // auto YY = Y1/Scale; - auto YY = Y1; - - // Use one of the built-in Eigen3 algorithms - // - /* - if (params["Jacobi"]) { - // -->Using Jacobi - Eigen::JacobiSVD - svd(YY, Eigen::ComputeThinU | Eigen::ComputeThinV); - S1 = svd.singularValues(); - V1 = svd.matrixV(); - } else if (params["BDCSVD"]) { - */ - // -->Using BDC - Eigen::BDCSVD - svd(YY, Eigen::ComputeFullU | Eigen::ComputeFullV); - // svd(YY, Eigen::ComputeThinU | Eigen::ComputeThinV); - S1 = svd.singularValues(); - V1 = svd.matrixV(); - /* - } else { - // -->Use Random approximation algorithm from Halko, Martinsson, - // and Tropp - int srank = std::min(YY.cols(), YY.rows()); - RedSVD::RedSVD svd(YY, srank); - S1 = svd.singularValues(); - V1 = svd.matrixV(); - } - */ - - std::cout << "shape V1 = " << V1.rows() << " x " - << V1.cols() << std::endl; - - std::cout << "shape Y1 = " << Y1.rows() << " x " - << Y1.cols() << std::endl; - - int lags = V1.rows(); - int rank = V1.cols(); - - std::ofstream out; - if (debug) out.open("debug.txt"); - - if (out) out << "rank=" << rank << " lags=" << lags << std::endl; - - VT1.resize(rank, lags-1); - VT2.resize(rank, lags-1); - - for (int j=0; j uu; - for (int i=0; iUsing Jacobi - Eigen::JacobiSVD - // svd(VT1, Eigen::ComputeThinU | Eigen::ComputeThinV); - svd(VT1, Eigen::ComputeFullU | Eigen::ComputeFullV); - SS = svd.singularValues(); - UU = svd.matrixU(); - VV = svd.matrixV(); - } else if (params["BDCSVD"]) { - */ - { - // -->Using BDC - Eigen::BDCSVD - // svd(VT1, Eigen::ComputeThinU | Eigen::ComputeThinV); - svd(VT1, Eigen::ComputeFullU | Eigen::ComputeFullV); - SS = svd.singularValues(); - UU = svd.matrixU(); - VV = svd.matrixV(); - } - /* - } else { - // -->Use Random approximation algorithm from Halko, Martinsson, - // and Tropp - // RedSVD::RedSVD svd(VT1, std::min(rank, numK-1)); - RedSVD::RedSVD svd(VT1, std::max(VT1.rows(), VT2.cols())); - SS = svd.singularValues(); - UU = svd.matrixU(); - VV = svd.matrixV(); - } - */ - - if (out) out << "Singular values" << std::endl << SS << std::endl; - - // Compute inverse - for (int i=0; i tol) SS(i) = 1.0/SS(i); - else SS(i) = 0.0; - } - - // Compute full Koopman operator - if (use_fullKh) { - - Eigen::MatrixXd DD(VV.cols(), UU.cols()); - DD.setZero(); - for (int i=0; i es(AT); - - L = es.eigenvalues(); - Phi = es.eigenvectors(); - - if (out) { - out << std::endl << "Eigenvalues" << std::endl << L << std::endl - << std::endl << "Eigenvectors" << std::endl << Phi << std::endl; - } - - } - // Compute the reduced Koopman operator - else { - - Eigen::MatrixXd AT = UU.transpose() * (VT2 * VV) * SS.asDiagonal(); - - // Compute spectrum - Eigen::EigenSolver es(AT, true); - - L = es.eigenvalues(); - auto W = es.eigenvectors(); - - // Compute the EDMD modes - // - Eigen::VectorXcd Linv(L); - for (int i=0; i tol) Linv(i) = 1.0/Linv(i); - else Linv(i) = 0.0; - } - - Phi = VT2 * VV * SS.asDiagonal() * W * Linv.asDiagonal(); - - if (out) { - out << std::endl << "Eigenvalues" << std::endl << L << std::endl - << std::endl << "Eigenvectors" << std::endl << Phi << std::endl; - } - } - - // Cache window size - // - window = D; - - return {L, Phi}; - } - - std::map - expMSSA::getReconstructedKoopman(int mode) - { - // Copy the original map for return - // - auto newdata = data; - - size_t n=0, offset=0; - - for (auto u : mean) { - - double disp = totVar; - if (type == TrendType::totPow) disp = totPow; - if (disp==0.0) disp = var[u.first]; - - std::complex phase = 1.0; - for (int i=0; i Date: Sun, 7 Jul 2024 11:20:51 -0700 Subject: [PATCH 20/40] Missing Koopman removal from the header file [no ci] --- expui/expMSSA.H | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/expui/expMSSA.H b/expui/expMSSA.H index a5cdb1606..1200f0d97 100644 --- a/expui/expMSSA.H +++ b/expui/expMSSA.H @@ -50,13 +50,6 @@ namespace MSSA //! Right singular vectors Eigen::MatrixXd U; - //@{ - //! Koopman modes - Eigen::VectorXcd L; - Eigen::MatrixXcd Phi; - int window; - //@} - //! Parameters //@{ bool flip, verbose, powerf; @@ -287,13 +280,6 @@ namespace MSSA return ret; } - //! Estimate Koopman modes from the trajectory eigenvectors - std::tuple - getKoopmanModes(const double tol, int window, bool debug); - - //! Return the reconstructed Koopman modes - std::map getReconstructedKoopman(int mode); - }; From 1dfb00c1c8b0ba0f62294cb990e27b6556581bb0 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Mon, 8 Jul 2024 13:37:51 -0700 Subject: [PATCH 21/40] sqrt not a constexpr in clang --- expui/FieldBasis.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc index 32533400d..158180dfa 100644 --- a/expui/FieldBasis.cc +++ b/expui/FieldBasis.cc @@ -313,7 +313,7 @@ namespace BasisClasses double u, double v, double w) { constexpr std::complex I(0, 1); - constexpr double fac0 = 1.0/sqrt(4*M_PI); + double fac0 = 1.0/sqrt(4*M_PI); int tid = omp_get_thread_num(); PS3 pos{x, y, z}, vel{u, v, w}; From 76c68cfc28a38eb5a62116f7fe5c1e6fb3e7349f Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Mon, 8 Jul 2024 17:53:55 -0700 Subject: [PATCH 22/40] Add C++17 workaround for sqrt constexpr --- expui/FieldBasis.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/FieldBasis.cc b/expui/FieldBasis.cc index 158180dfa..7c85e0aee 100644 --- a/expui/FieldBasis.cc +++ b/expui/FieldBasis.cc @@ -313,7 +313,7 @@ namespace BasisClasses double u, double v, double w) { constexpr std::complex I(0, 1); - double fac0 = 1.0/sqrt(4*M_PI); + constexpr double fac0 = 0.25*M_2_SQRTPI; int tid = omp_get_thread_num(); PS3 pos{x, y, z}, vel{u, v, w}; From 23fbb21334233f203f5616acd21eb26d846a5916 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 8 Jul 2024 23:25:21 -0700 Subject: [PATCH 23/40] Updates to 'make tests' scripts and config [no ci] --- tests/CMakeLists.txt | 6 +++--- tests/Halo/check.py | 7 +++++-- tests/Halo/config.yml | 5 +++-- 3 files changed, 11 insertions(+), 7 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 803800330..5928b4395 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -30,14 +30,14 @@ if(ENABLE_PYEXP) WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) add_test(NAME removeCylCache - COMMAND ${CMAKE_COMMAND} -E remove .eof.cache.run0t test_adddisk_sl.0 + COMMAND ${CMAKE_COMMAND} -E remove .eof.cache.run0t WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Disk) set_tests_properties(removeSphCache PROPERTIES DEPENDS pyexpSphBasisTest REQUIRED_FILES ".slgrid_sph_cache") set_tests_properties(removeCylCache PROPERTIES DEPENDS pyexpCylBasisTest - REQUIRED_FILES ".eof.cache.run0t;test_adddisk_sl.0") + REQUIRED_FILES ".eof.cache.run0t") # Other tests for pyEXP go here ... @@ -56,7 +56,7 @@ if(ENABLE_NBODY) # Makes some spherical ICs using utils/ICs/gensph add_test(NAME makeICTest - COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/utils/ICs/gensph -N 1000 -i SLGridSph.model + COMMAND ${EXP_MPI_LAUNCH} ${CMAKE_BINARY_DIR}/utils/ICs/gensph -N 10000 -i SLGridSph.model WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) # Runs those ICs using exp diff --git a/tests/Halo/check.py b/tests/Halo/check.py index 4299f403b..639b02b27 100644 --- a/tests/Halo/check.py +++ b/tests/Halo/check.py @@ -3,12 +3,15 @@ # Read the output log file data = np.loadtxt("OUTLOG.run0", skiprows=6, delimiter="|") -# Column 16 is -2T/VC. The mean should be 1 with a std dev < 0.03 +# Column 16 is -2T/VC. The mean should be 1 with some small number of std dev mean = np.mean(data[:,16]) stdv = np.std (data[:,16]) +# print("Mean = " + str(mean) + ", std dev = " + str(stdv)) +# print("Check = {}".format(np.abs(mean-1.0) - 5.0*stdv)) + # If the values are within 3 sigma, assume that the simulation worked -if np.abs(mean - 1.0) > 3.0*stdv: +if np.abs(mean - 1.0) > 5.0*stdv: exit(1) else: exit(0) diff --git a/tests/Halo/config.yml b/tests/Halo/config.yml index 3743a65e6..bc83d3dc0 100644 --- a/tests/Halo/config.yml +++ b/tests/Halo/config.yml @@ -8,9 +8,9 @@ # ------------------------------------------------------------------------ Global: nthrds : 1 - dtime : 0.005 + dtime : 0.002 runtag : run0 - nsteps : 200 + nsteps : 500 multistep : 4 dynfracV : 0.01 dynfracA : 0.03 @@ -40,6 +40,7 @@ Components: rmapping : 0.0667 self_consistent: true modelname: SLGridSph.model + cachename: SLGridSph.cache.run0 # ------------------------------------------------------------------------ # This is a sequence of outputs From 2e89c5dd856bcbf0eb321c32bcb1657772850637 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Tue, 9 Jul 2024 07:38:50 -0700 Subject: [PATCH 24/40] Testing testing 1 2 3 --- .github/workflows/build.yml | 50 +++++-------------------------------- 1 file changed, 6 insertions(+), 44 deletions(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index a37377ef8..acdf4f4d2 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -8,7 +8,7 @@ on: branches: - main jobs: - pyexp: + exp: strategy: matrix: os: [ubuntu-latest] @@ -31,14 +31,14 @@ jobs: git submodule update --init --recursive mkdir -p build/install - - name: Compile pyEXP + - name: Compile EXP if: runner.os == 'Linux' env: CC: ${{ matrix.cc }} working-directory: ./build run: >- cmake - -DENABLE_NBODY=NO + -DENABLE_NBODY=YES -DENABLE_PYEXP=YES -DCMAKE_BUILD_TYPE=Release -DEigen3_DIR=/usr/include/eigen3/share/eigen3/cmake @@ -50,46 +50,8 @@ jobs: working-directory: ./build run: make -j 2 - # ----------------------------------------------------------------------------------- - - exp: - strategy: - matrix: - os: [ubuntu-latest] - cc: [gcc, mpicc] - - name: "Test Full EXP Build" - runs-on: ${{ matrix.os }} - steps: - - name: Checkout - uses: actions/checkout@v3 - - - name: Install core dependencies - ubuntu - if: runner.os == 'Linux' - run: | - sudo apt-get update - sudo apt-get install -y build-essential libeigen3-dev libfftw3-dev libhdf5-dev libopenmpi-dev - - - name: Setup submodule and build - run: | - git submodule update --init --recursive - mkdir -p build/install - - - name: Compile Full EXP - Linux - if: runner.os == 'Linux' - env: - CC: ${{ matrix.cc }} + - name: CTest working-directory: ./build - run: >- - cmake - -DENABLE_NBODY=YES - -DENABLE_PYEXP=NO - -DCMAKE_BUILD_TYPE=Release - -DEigen3_DIR=/usr/include/eigen3/share/eigen3/cmake - -DCMAKE_INSTALL_PREFIX=./install - -Wno-dev - .. + run: ctest -j 2 -L quick + - - name: Make - working-directory: ./build - run: make -j 2 From 93c48f285934b4fd5b633adb6b17a28318e4d877 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Tue, 9 Jul 2024 07:47:29 -0700 Subject: [PATCH 25/40] rm old build action [no ci] --- .github/workflows/buildfull.yml | 146 -------------------------------- 1 file changed, 146 deletions(-) delete mode 100644 .github/workflows/buildfull.yml diff --git a/.github/workflows/buildfull.yml b/.github/workflows/buildfull.yml deleted file mode 100644 index eb9c86552..000000000 --- a/.github/workflows/buildfull.yml +++ /dev/null @@ -1,146 +0,0 @@ -name: "Test Builds" - -on: - -jobs: - pyexp: - strategy: - matrix: - os: [macos-latest, ubuntu-latest] - cc: [gcc, mpicc] - - name: "Test pyEXP Build" - runs-on: ${{ matrix.os }} - steps: - - name: Checkout - uses: actions/checkout@v3 - - - name: Install core dependencies - ubuntu - if: runner.os == 'Linux' - run: | - sudo apt-get update - sudo apt-get install -y build-essential libeigen3-dev libfftw3-dev libhdf5-dev libopenmpi-dev - - - name: Install core dependencies - mac - if: startsWith(matrix.os, 'mac') - run: | - brew update - brew reinstall gcc - brew install eigen fftw hdf5 open-mpi libomp - - - name: Setup submodule and build - run: | - git submodule update --init --recursive - mkdir -p build/install - - - name: Compile pyEXP - Linux - if: runner.os == 'Linux' - env: - CC: ${{ matrix.cc }} - working-directory: ./build - run: >- - cmake - -DENABLE_NBODY=NO - -DENABLE_PYEXP=YES - -DCMAKE_BUILD_TYPE=Release - -DEigen3_DIR=/usr/include/eigen3/share/eigen3/cmake - -DCMAKE_INSTALL_PREFIX=./install - -Wno-dev - .. - - # Note for future: The homebrew paths are for intel only. Once ARM macs are - # supported in here, we'll need to update to /opt/homebrew/... instead - - name: Compile pyEXP - Mac - if: startsWith(matrix.os, 'mac') - env: - CC: ${{ matrix.cc }} - LDFLAGS: -L/usr/local/opt/libomp/lib - CPPFLAGS: -I/usr/local/opt/libomp/include - working-directory: ./build - run: >- - cmake - -DENABLE_NBODY=NO - -DENABLE_PYEXP=YES - -DCMAKE_BUILD_TYPE=Release - -DEigen3_DIR=/usr/local/share/eigen3/cmake - -DCMAKE_INSTALL_PREFIX=./install - -DOpenMP_CXX_INCLUDE_DIR=/usr/local/opt/libomp/include - -DOpenMP_C_INCLUDE_DIR=/usr/local/opt/libomp/include - -Wno-dev - .. - - - name: Make - working-directory: ./build - run: make -j 2 - - # ----------------------------------------------------------------------------------- - - exp: - strategy: - matrix: - os: [macos-latest, ubuntu-latest] - cc: [gcc, mpicc] - - name: "Test Full EXP Build" - runs-on: ${{ matrix.os }} - steps: - - name: Checkout - uses: actions/checkout@v3 - - - name: Install core dependencies - ubuntu - if: runner.os == 'Linux' - run: | - sudo apt-get update - sudo apt-get install -y build-essential libeigen3-dev libfftw3-dev libhdf5-dev libopenmpi-dev - - - name: Install core dependencies - mac - if: startsWith(matrix.os, 'mac') - run: | - brew update - brew reinstall gcc - brew install eigen fftw hdf5 open-mpi libomp - - - name: Setup submodule and build - run: | - git submodule update --init --recursive - mkdir -p build/install - - - name: Compile Full EXP - Linux - if: runner.os == 'Linux' - env: - CC: ${{ matrix.cc }} - working-directory: ./build - run: >- - cmake - -DENABLE_NBODY=YES - -DENABLE_PYEXP=NO - -DCMAKE_BUILD_TYPE=Release - -DEigen3_DIR=/usr/include/eigen3/share/eigen3/cmake - -DCMAKE_INSTALL_PREFIX=./install - -Wno-dev - .. - - # Note for future: The homebrew paths are for intel only. Once ARM macs are - # supported in here, we'll need to update to /opt/homebrew/... instead - - name: Compile Full EXP - Mac - if: startsWith(matrix.os, 'mac') - env: - CC: ${{ matrix.cc }} - LDFLAGS: -L/usr/local/opt/libomp/lib - CPPFLAGS: -I/usr/local/opt/libomp/include - working-directory: ./build - run: >- - cmake - -DENABLE_NBODY=YES - -DENABLE_PYEXP=NO - -DCMAKE_BUILD_TYPE=Release - -DEigen3_DIR=/usr/local/share/eigen3/cmake - -DCMAKE_INSTALL_PREFIX=./install - -DOpenMP_CXX_INCLUDE_DIR=/usr/local/opt/libomp/include - -DOpenMP_C_INCLUDE_DIR=/usr/local/opt/libomp/include - -Wno-dev - .. - - - name: Make - working-directory: ./build - run: make -j 2 From b02f556c0c0f636e608a10c5ae1d257dd457529e Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 9 Jul 2024 08:52:50 -0700 Subject: [PATCH 26/40] Remove temporary files automatically --- utils/ICs/gensph.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/utils/ICs/gensph.cc b/utils/ICs/gensph.cc index f685133c1..d921c883d 100644 --- a/utils/ICs/gensph.cc +++ b/utils/ICs/gensph.cc @@ -809,6 +809,7 @@ main(int argc, char **argv) if (myid==0) { std::ostringstream sout; sout << "cat " << OUTPS << ".* > " << OUTPS; + sout << "; rm " << OUTPS << ".*"; int ret = system(sout.str().c_str()); } From b2041ad709e3ef4f694b31ad64d8e556a38bc7f6 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 9 Jul 2024 08:53:19 -0700 Subject: [PATCH 27/40] A few minor updates and clean ups [no ci] --- tests/CMakeLists.txt | 6 +++--- tests/Halo/check.py | 9 +++------ 2 files changed, 6 insertions(+), 9 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 5928b4395..f528ca5b0 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -92,14 +92,14 @@ if(ENABLE_NBODY) # perhaps there is a better way? add_test(NAME removeTempFiles COMMAND ${CMAKE_COMMAND} -E remove - config.run0.yml current.processor.rates.run0 new.bods new.bods.0 + config.run0.yml current.processor.rates.run0 new.bods OUTLOG.run0 run0.levels SLGridSph.cache.run0 test.grid 'outcoef.dark halo.run0' WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) # Remove the temporary files set_tests_properties(removeTempFiles PROPERTIES DEPENDS expNbodyCheck2TW - REQUIRED_FILES "config.run0.yml;current.processor.rates.run0;new.bods;new.bods.0;OUTLOG.run0;run0.levels;SLGridSph.cache.run0;test.grid;" + REQUIRED_FILES "config.run0.yml;current.processor.rates.run0;new.bods;run0.levels;SLGridSph.cache.run0;test.grid;" ) # Makes some cube ICs using utils/ICs/cubeics @@ -137,7 +137,7 @@ if(ENABLE_NBODY) # Set labels for pyEXP tests set_tests_properties(expExecuteTest PROPERTIES LABELS "quick") set_tests_properties(makeICTest expNbodyTest expNbodyCheck2TW - removeTempFiles expCubeTest removeCubeFiles + removeTempFiles makeCubeICTest expCubeTest removeCubeFiles PROPERTIES LABELS "long") endif() diff --git a/tests/Halo/check.py b/tests/Halo/check.py index 639b02b27..54a499349 100644 --- a/tests/Halo/check.py +++ b/tests/Halo/check.py @@ -3,15 +3,12 @@ # Read the output log file data = np.loadtxt("OUTLOG.run0", skiprows=6, delimiter="|") -# Column 16 is -2T/VC. The mean should be 1 with some small number of std dev +# Column 16 is -2T/VC. The mean should be 1 mean = np.mean(data[:,16]) stdv = np.std (data[:,16]) -# print("Mean = " + str(mean) + ", std dev = " + str(stdv)) -# print("Check = {}".format(np.abs(mean-1.0) - 5.0*stdv)) - -# If the values are within 3 sigma, assume that the simulation worked -if np.abs(mean - 1.0) > 5.0*stdv: +# If the values are within 6 sigma of 1, assume that the simulation worked +if np.abs(mean - 1.0) > 6.0*stdv: exit(1) else: exit(0) From 78b41f0e362049d3e3c1f569a5eae38a9fe221c0 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 9 Jul 2024 09:02:11 -0700 Subject: [PATCH 28/40] Change the component name to prevent a coefficient file name that contains a space [no ci] --- tests/CMakeLists.txt | 2 +- tests/Halo/config.yml | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index f528ca5b0..42d84db61 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -94,7 +94,7 @@ if(ENABLE_NBODY) COMMAND ${CMAKE_COMMAND} -E remove config.run0.yml current.processor.rates.run0 new.bods OUTLOG.run0 run0.levels SLGridSph.cache.run0 test.grid - 'outcoef.dark halo.run0' + outcoef.halo.run0 WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) # Remove the temporary files diff --git a/tests/Halo/config.yml b/tests/Halo/config.yml index bc83d3dc0..596d51f36 100644 --- a/tests/Halo/config.yml +++ b/tests/Halo/config.yml @@ -26,7 +26,7 @@ Global: # Each indented stanza beginning with '-' is a component # ------------------------------------------------------------------------ Components: - - name : dark halo + - name : halo parameters : {nlevel: 1, indexing: true} bodyfile : new.bods force : @@ -49,7 +49,7 @@ Output: - id : outlog parameters : {nint: 10} - id : outcoef - parameters : {nint: 1, name: dark halo} + parameters : {nint: 1, name: halo} # ------------------------------------------------------------------------ # This is a sequence of external forces From 60d061e9e8741ad436cc398c68bd67b1ed959707 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 9 Jul 2024 09:32:50 -0700 Subject: [PATCH 29/40] Change component name to remove spaces from the file name; tested all label categories successfully [no ci] --- tests/CMakeLists.txt | 4 ++-- tests/Halo/readCoefs.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 42d84db61..81efdca36 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -94,12 +94,12 @@ if(ENABLE_NBODY) COMMAND ${CMAKE_COMMAND} -E remove config.run0.yml current.processor.rates.run0 new.bods OUTLOG.run0 run0.levels SLGridSph.cache.run0 test.grid - outcoef.halo.run0 + outcoef.halo.run0 halo.cache WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) # Remove the temporary files set_tests_properties(removeTempFiles PROPERTIES DEPENDS expNbodyCheck2TW - REQUIRED_FILES "config.run0.yml;current.processor.rates.run0;new.bods;run0.levels;SLGridSph.cache.run0;test.grid;" + REQUIRED_FILES "config.run0.yml;current.processor.rates.run0;new.bods;run0.levels;halo.cache;test.grid;" ) # Makes some cube ICs using utils/ICs/cubeics diff --git a/tests/Halo/readCoefs.py b/tests/Halo/readCoefs.py index ad7ab60e2..0004001f7 100644 --- a/tests/Halo/readCoefs.py +++ b/tests/Halo/readCoefs.py @@ -3,7 +3,7 @@ import pyEXP -coefs = pyEXP.coefs.Coefs.factory('outcoef.dark halo.run0') +coefs = pyEXP.coefs.Coefs.factory('outcoef.halo.run0') data = coefs.getAllCoefs() print(data.shape) print(coefs.getName()) From a78ded51f7c41bd8f2a8a07cb0668dd95d392ead Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 9 Jul 2024 10:13:17 -0700 Subject: [PATCH 30/40] Typo in temp file removal list. CTest runs should be good to go at this point. [no ci] --- tests/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 81efdca36..9aee83328 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -94,12 +94,12 @@ if(ENABLE_NBODY) COMMAND ${CMAKE_COMMAND} -E remove config.run0.yml current.processor.rates.run0 new.bods OUTLOG.run0 run0.levels SLGridSph.cache.run0 test.grid - outcoef.halo.run0 halo.cache + outcoef.halo.run0 SLGridSph.cache.run0 WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/Halo) # Remove the temporary files set_tests_properties(removeTempFiles PROPERTIES DEPENDS expNbodyCheck2TW - REQUIRED_FILES "config.run0.yml;current.processor.rates.run0;new.bods;run0.levels;halo.cache;test.grid;" + REQUIRED_FILES "config.run0.yml;current.processor.rates.run0;new.bods;run0.levels;SLGridSph.cache.run0;test.grid;" ) # Makes some cube ICs using utils/ICs/cubeics From 824a13e1f15be8d37257d60b453ccd904f9c8eca Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Tue, 9 Jul 2024 10:35:02 -0700 Subject: [PATCH 31/40] updates to ctest suite --- .github/workflows/build.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index acdf4f4d2..a6514ca58 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -12,7 +12,7 @@ jobs: strategy: matrix: os: [ubuntu-latest] - cc: [gcc, mpicc] + cc: [gcc] name: "Test pyEXP Build" runs-on: ${{ matrix.os }} From 68ed00e81a5a183ac9b48090b17f7ea145b2754d Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Tue, 9 Jul 2024 16:49:30 -0700 Subject: [PATCH 32/40] numpy dependence tweaks --- .github/workflows/build.yml | 1 + tests/Disk/cyl_basis.py | 1 - tests/Halo/sph_basis.py | 1 - 3 files changed, 1 insertion(+), 2 deletions(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index a6514ca58..b966d0d4a 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -25,6 +25,7 @@ jobs: run: | sudo apt-get update sudo apt-get install -y build-essential libeigen3-dev libfftw3-dev libhdf5-dev libopenmpi-dev + sudo pip install numpy - name: Setup submodule and build run: | diff --git a/tests/Disk/cyl_basis.py b/tests/Disk/cyl_basis.py index 227b22f4a..c3c55b4e4 100644 --- a/tests/Disk/cyl_basis.py +++ b/tests/Disk/cyl_basis.py @@ -1,7 +1,6 @@ #!/usr/bin/env python # coding: utf-8 -import numpy as np import pyEXP # Make the disk basis config diff --git a/tests/Halo/sph_basis.py b/tests/Halo/sph_basis.py index 1cee5d8ac..f7b6a18b9 100644 --- a/tests/Halo/sph_basis.py +++ b/tests/Halo/sph_basis.py @@ -1,7 +1,6 @@ #!/usr/bin/env python # coding: utf-8 -import numpy as np import pyEXP # Make the halo basis config From 7ffd4b8554a0d8ef204bfaa208cdf2423c49469a Mon Sep 17 00:00:00 2001 From: mdw Date: Tue, 9 Jul 2024 20:34:01 -0400 Subject: [PATCH 33/40] Fixes to remove numpy --- tests/Disk/cyl_basis.py | 2 +- tests/Halo/check.py | 22 ++++++++++++++-------- tests/Halo/sph_basis.py | 2 +- 3 files changed, 16 insertions(+), 10 deletions(-) diff --git a/tests/Disk/cyl_basis.py b/tests/Disk/cyl_basis.py index 227b22f4a..2d2b6e0d3 100644 --- a/tests/Disk/cyl_basis.py +++ b/tests/Disk/cyl_basis.py @@ -1,7 +1,7 @@ #!/usr/bin/env python # coding: utf-8 -import numpy as np +# import numpy as np import pyEXP # Make the disk basis config diff --git a/tests/Halo/check.py b/tests/Halo/check.py index 54a499349..f181081d9 100644 --- a/tests/Halo/check.py +++ b/tests/Halo/check.py @@ -1,14 +1,20 @@ -import numpy as np +# Open the output log file +file = open("OUTLOG.run0") -# Read the output log file -data = np.loadtxt("OUTLOG.run0", skiprows=6, delimiter="|") +n = 0 # Count lines +mean = 0.0 # Accumulate 2T/VC valkues -# Column 16 is -2T/VC. The mean should be 1 -mean = np.mean(data[:,16]) -stdv = np.std (data[:,16]) +# Open the output log file +while (line := file.readline()) != "": + if n >= 6: # Skip the header stuff + v = [float(x) for x in line.split('|')] + mean += v[16] # This is the 2T/VC column + n = n + 1 # Count lines -# If the values are within 6 sigma of 1, assume that the simulation worked -if np.abs(mean - 1.0) > 6.0*stdv: +if n>6: mean /= n-6 # Sanity check + +# Check closeness to 1.0 +if (mean-1.0)*(mean-1.0) > 0.003: exit(1) else: exit(0) diff --git a/tests/Halo/sph_basis.py b/tests/Halo/sph_basis.py index 1cee5d8ac..e72ca54e5 100644 --- a/tests/Halo/sph_basis.py +++ b/tests/Halo/sph_basis.py @@ -1,7 +1,7 @@ #!/usr/bin/env python # coding: utf-8 -import numpy as np +# import numpy as np import pyEXP # Make the halo basis config From 8c100f02ce3ea6906deeacadabe173f757cc53a1 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 9 Jul 2024 17:59:51 -0700 Subject: [PATCH 34/40] Toggle off HIGH_UNIT_TESTS manually --- CMakeLists.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index cd65e7a1a..dbede3fe3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -228,7 +228,8 @@ add_subdirectory(extern/pybind11) # Set options for the HighFive git submodule in extern set(HIGHFIVE_EXAMPLES OFF CACHE BOOL "Do not build the examples") set(HIGHFIVE_BUILD_DOCS OFF CACHE BOOL "Do not build the documentation") -set(HIGHFIVE_USE_BOOST OFF CACHE BOOL "Do not use Boost in HighFIve") +set(HIGHFIVE_USE_BOOST OFF CACHE BOOL "Do not use Boost in HighFive") +set(HIGHFIVE_UNIT_TESTS OFF CACHE BOOL "Turn off internal testing for HighFIve") set(H5_USE_EIGEN TRUE CACHE BOOL "Eigen3 support in HighFive") add_subdirectory(extern/HighFive EXCLUDE_FROM_ALL) From 0de5e0d3a27031cfdb48697eaeb32d5c33a3e67b Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 9 Jul 2024 19:10:19 -0700 Subject: [PATCH 35/40] Fix comment spelling typo [no ci] --- tests/Halo/check.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tests/Halo/check.py b/tests/Halo/check.py index f181081d9..d05381c3d 100644 --- a/tests/Halo/check.py +++ b/tests/Halo/check.py @@ -2,9 +2,10 @@ file = open("OUTLOG.run0") n = 0 # Count lines -mean = 0.0 # Accumulate 2T/VC valkues +mean = 0.0 # Accumulate 2T/VC values # Open the output log file +# while (line := file.readline()) != "": if n >= 6: # Skip the header stuff v = [float(x) for x in line.split('|')] @@ -14,6 +15,7 @@ if n>6: mean /= n-6 # Sanity check # Check closeness to 1.0 +# if (mean-1.0)*(mean-1.0) > 0.003: exit(1) else: From af674567bd8705123bbc01a1620a46610b368de8 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Tue, 9 Jul 2024 20:22:27 -0700 Subject: [PATCH 36/40] adjustments to build action --- .github/workflows/build.yml | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index b966d0d4a..82c76f90e 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -7,6 +7,7 @@ on: pull_request: branches: - main + - devel jobs: exp: strategy: @@ -49,10 +50,12 @@ jobs: - name: Make working-directory: ./build - run: make -j 2 + run: make -j 4 - - name: CTest + - name: CTest Quick working-directory: ./build - run: ctest -j 2 -L quick + run: ctest -j 4 -L quick - + - name: CTest Long + working-directory: ./build + run: ctest -j 4 -L long From 65dd03aabe4e38a1b2f6cc152cf9fa9f488c185f Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 9 Jul 2024 21:05:25 -0700 Subject: [PATCH 37/40] Remove numpy dependence from Cube OUTLOG check [no ci] --- tests/Cube/check.py | 29 ++++++++++++++++++++--------- 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/tests/Cube/check.py b/tests/Cube/check.py index 1956e0ee6..b6b0a3eef 100644 --- a/tests/Cube/check.py +++ b/tests/Cube/check.py @@ -1,15 +1,26 @@ -import numpy as np +# open the output log file +file = open("OUTLOG.runS") -# Read the output log file -data = np.loadtxt("OUTLOG.runS", skiprows=6, delimiter="|") +n = 0 # Count lines +mean = [0.0, 0.0, 0.0] # Mean positions -# Columns 4, 5, 6 is mean position -x = np.mean(data[:,3]) -y = np.mean(data[:,4]) -z = np.mean(data[:,5]) +# Open the output log file +# +while (line := file.readline()) != "": + if n >= 6: # Skip the header stuff + v = [float(x) for x in line.split('|')] + mean[0] += v[3] # x pos + mean[1] += v[4] # y pos + mean[2] += v[5] # z pos + n = n + 1 # Count lines -# If the values are close to 0.5, assume it worked -if np.abs(x - 0.5) > 0.15 or np.abs(y - 0.5) > 0.15 or np.abs(z - 0.5) > 0.15: +if n>6: # Sanity check + for i in range(3): + mean[i] = mean[i]/(n-6) - 0.5 + +# If the squared values are close to 0.0, assume it worked +# +if mean[0]*mean[0] > 0.03 or mean[1]*mean[1] > 0.03 or mean[2]*mean[2] > 0.03: exit(1) else: exit(0) From 2a42b53dace10c0532fbf9ca9d6453f0d63e91e5 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Wed, 10 Jul 2024 07:13:40 -0700 Subject: [PATCH 38/40] Disable long tests for actions --- .github/workflows/build.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 82c76f90e..71dfb58ba 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -56,6 +56,6 @@ jobs: working-directory: ./build run: ctest -j 4 -L quick - - name: CTest Long - working-directory: ./build - run: ctest -j 4 -L long + #- name: CTest Long + #working-directory: ./build + #run: ctest -j 4 -L long From 4572df803227c43af29b46460a7d2ae3a5b80166 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Wed, 10 Jul 2024 11:36:06 -0700 Subject: [PATCH 39/40] Rename runner --- .github/workflows/build.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 71dfb58ba..bfad99340 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -1,4 +1,4 @@ -name: "Test Builds" +name: "Build and Test" on: push: From 7758e867268bcd5fbad5ea5cc0db9ad60c29d6c0 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Wed, 10 Jul 2024 11:55:53 -0700 Subject: [PATCH 40/40] run tests in serial [no ci] --- .github/workflows/build.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index bfad99340..b9ddde223 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -54,8 +54,8 @@ jobs: - name: CTest Quick working-directory: ./build - run: ctest -j 4 -L quick + run: ctest -L quick #- name: CTest Long #working-directory: ./build - #run: ctest -j 4 -L long + #run: ctest -L long