From 494c619cc1cfc7e7d38bc48cdb8a5593f6d416ff Mon Sep 17 00:00:00 2001 From: mdw Date: Thu, 12 Mar 2026 10:00:36 -0400 Subject: [PATCH 01/78] Added a Toomre-type disk for deprojection primarily --- include/DiskModels.H | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/include/DiskModels.H b/include/DiskModels.H index 3b46da803..273638d4f 100644 --- a/include/DiskModels.H +++ b/include/DiskModels.H @@ -196,6 +196,35 @@ public: } }; +//! Toomre disk +class Toomre : public EmpCylSL::AxiDisk +{ +private: + double a, h, n; + double norm; + +public: + + Toomre(double a, double h, double n=3, double M=1) : + a(a), h(h), n(n), EmpCylSL::AxiDisk(M, "Toomre") + { + params.push_back(a); + params.push_back(h); + params.push_back(n); + if (n <= 2) throw std::runtime_error("Toomre index must be > 2"); + norm = (n - 2.0)/(2.0*M_PI*a*a); + } + + double operator()(double R, double z, double phi=0.) + { + double sigma = std::pow(1.0 * (R/a)*(R/a), -0.5*n) * norm; + double Z = std::fabs(z); + double Q = std::exp(-Z/h); + double sech = 2.0*Q/(1.0 + Q*Q); + return sigma * sech*sech * 0.25/h; + } +}; + //! Truncate a AxiDisk class Truncated : public EmpCylSL::AxiDisk { From c7bf5eeca05a244a3a78fe758660fb475aa46a4d Mon Sep 17 00:00:00 2001 From: mdw Date: Thu, 12 Mar 2026 10:01:30 -0400 Subject: [PATCH 02/78] Added an enum and reflector for deprojection disk types and included the Toomre model in the updated list of valid types --- expui/BiorthBasis.H | 19 +++++++++++++++-- expui/BiorthBasis.cc | 50 ++++++++++++++++++++++++++++++++++++++------ 2 files changed, 61 insertions(+), 8 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index f40875985..054e28a35 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -999,8 +999,9 @@ namespace BasisClasses std::string pyname; //! DiskType support - // - enum DiskType { constant, gaussian, mn, exponential, doubleexpon, diskbulge, python }; + //! + enum class DiskType + { constant, gaussian, mn, exponential, doubleexpon, diskbulge, python }; DiskType DTYPE; std::string mtype, dtype, dmodel; @@ -1009,6 +1010,20 @@ namespace BasisClasses double dcond(double R, double z, double phi, int M); std::shared_ptr pyDens; + //@{ + //! DeprojType support + + //! Disk models used for deprojection + enum class DeprojType + { mn, toomre, python, exponential}; + + //! Current model + DeprojType PTYPE; + + //! Look up by string + static const std::map dplookup; + //@} + protected: //! Evaluate basis in spherical coordinates diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index d8574afcb..fda83a3ef 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1198,6 +1198,14 @@ namespace BasisClasses {"python", DiskType::python} }; + // Dprojection model for cylindrical basis construction + const std::map Cylindrical::dplookup = + { {"mn", DeprojType::mn}, + {"exponential", DeprojType::exponential}, + {"toomre", DeprojType::toomre}, + {"python", DeprojType::python} + }; + const std::set Cylindrical::valid_keys = { "tk_type", @@ -1233,6 +1241,7 @@ namespace BasisClasses "expcond", "ignore", "deproject", + "dmodel", "logr", "pcavar", "pcaeof", @@ -1494,7 +1503,7 @@ namespace BasisClasses if (conf["cmapz" ]) cmapZ = conf["cmapz" ].as(); if (conf["ignore" ]) Ignore = conf["ignore" ].as(); if (conf["deproject" ]) deproject = conf["deproject" ].as(); - if (conf["dmodel" ]) dmodel = conf["dmodel" ].as(); + if (conf["dmodel" ]) dmodel = conf["dmodel" ].as(); if (conf["aratio" ]) aratio = conf["aratio" ].as(); if (conf["hratio" ]) hratio = conf["hratio" ].as(); @@ -1697,16 +1706,45 @@ namespace BasisClasses // EmpCylSL::AxiDiskPtr model; - if (dmodel.compare("MN")==0) // Miyamoto-Nagai + // Convert dmodel string to lower case + // + std::transform(dmodel.begin(), dmodel.end(), dmodel.begin(), + [](unsigned char c){ return std::tolower(c); }); + + + // Check for map entry + try { + PTYPE = dplookup.at(dmodel); + + // Report DeprojType + if (myid==0) { + std::cout << "---- Deprojection type is <" << dmodel + << ">" << std::endl; + } + } + catch (const std::out_of_range& err) { + if (myid==0) { + std::cout << "DeprojType error in configuraton file" << std::endl; + std::cout << "Valid options are: "; + for (auto v : dplookup) std::cout << v.first << " "; + std::cout << std::endl; + } + throw std::runtime_error("Cylindrical::initialize: invalid DiskModel"); + } + + if (PTYPE == DeprojType::mn) // Miyamoto-Nagai model = std::make_shared(1.0, H); - else if (DTYPE == DiskType::python) { + else if (PTYPE == DeprojType::toomre) { + model = std::make_shared(1.0, H); + } else if (PTYPE == DeprojType::python and + DTYPE == DiskType::python) { model = std::make_shared(pyname, acyl); std::cout << "Using AxiSymPyModel for deprojection from Python function <" << pyname << ">" << std::endl; - } - else // Default to exponential + } else { // Default to exponential model = std::make_shared(1.0, H); - + } + if (rwidth>0.0) { model = std::make_shared(rtrunc/acyl, rwidth/acyl, From 6efc02486a20147ba7d3594f96aec4e031217601 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 13 Mar 2026 12:05:29 -0400 Subject: [PATCH 03/78] Addded a test routine for Abel deprojection; updated EmpCylSL for derivative form, which is more accurate near the center --- exputil/EmpCylSL.cc | 33 ++++-- include/EmpCylSL.H | 6 + utils/Test/CMakeLists.txt | 8 +- utils/Test/CubicSpline.H | 31 +++++ utils/Test/CubicSpline.cc | 73 ++++++++++++ utils/Test/Deprojector.H | 58 +++++++++ utils/Test/Deprojector.cc | 200 ++++++++++++++++++++++++++++++++ utils/Test/testDeprojPlummer.cc | 76 ++++++++++++ utils/Test/testDeproject.cc | 63 ++++++++++ 9 files changed, 539 insertions(+), 9 deletions(-) create mode 100644 utils/Test/CubicSpline.H create mode 100644 utils/Test/CubicSpline.cc create mode 100644 utils/Test/Deprojector.H create mode 100644 utils/Test/Deprojector.cc create mode 100644 utils/Test/testDeprojPlummer.cc create mode 100644 utils/Test/testDeproject.cc diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index a89afa68e..e5bba4a83 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -487,28 +487,45 @@ void EmpCylSL::create_deprojection(double H, double Rf, int NUMR, int NINT, // Now, compute Abel inverion integral // for (int i=0; i rho(NUMR), mass(NUMR); - - Linear1d intgr(rl, rhoI); + Linear1d intgr; + if (abel_type == AbelType::IBP) intgr = Linear1d(rl, rhoI); - for (int i=0; i rho(NUMR), mass(NUMR); + for (int i=0; i + +namespace Deproject +{ + + class CubicSpline { + public: + CubicSpline() = default; + CubicSpline(const std::vector& x_in, const std::vector& y_in); + + // set data (x must be strictly increasing) + void set_data(const std::vector& x_in, const std::vector& y_in); + + // evaluate spline and its derivative (xx should lie within [xmin(), xmax()], but endpoints are clamped) + double eval(double xx) const; + double deriv(double xx) const; + + double xmin() const; + double xmax() const; + + private: + std::vector x_, y_, y2_; + int locate(double xx) const; + }; + +} + +#endif // CUBIC_SPLINE_H diff --git a/utils/Test/CubicSpline.cc b/utils/Test/CubicSpline.cc new file mode 100644 index 000000000..c88f899db --- /dev/null +++ b/utils/Test/CubicSpline.cc @@ -0,0 +1,73 @@ +#include + +#include "CubicSpline.H" + +namespace Deproject +{ + + CubicSpline::CubicSpline(const std::vector& x_in, const std::vector& y_in) { + set_data(x_in, y_in); + } + + void CubicSpline::set_data(const std::vector& x_in, const std::vector& y_in) { + x_ = x_in; + y_ = y_in; + int n = (int)x_.size(); + if (n < 2 || (int)y_.size() != n) throw std::runtime_error("CubicSpline: need at least two points and equal-sized x,y."); + + y2_.assign(n, 0.0); + std::vector u(n - 1, 0.0); + + // natural spline boundary conditions (second derivatives at endpoints = 0) + for (int i = 1; i < n - 1; ++i) { + double sig = (x_[i] - x_[i-1]) / (x_[i+1] - x_[i-1]); + double p = sig * y2_[i-1] + 2.0; + y2_[i] = (sig - 1.0) / p; + double dY1 = (y_[i+1] - y_[i]) / (x_[i+1] - x_[i]); + double dY0 = (y_[i] - y_[i-1]) / (x_[i] - x_[i-1]); + u[i] = (6.0 * (dY1 - dY0) / (x_[i+1] - x_[i-1]) - sig * u[i-1]) / p; + } + + for (int k = n - 2; k >= 0; --k) y2_[k] = y2_[k] * y2_[k+1] + u[k]; + } + + int CubicSpline::locate(double xx) const { + int n = (int)x_.size(); + if (xx <= x_.front()) return 0; + if (xx >= x_.back()) return n - 2; + int lo = 0, hi = n - 1; + while (hi - lo > 1) { + int mid = (lo + hi) >> 1; + if (x_[mid] > xx) hi = mid; else lo = mid; + } + return lo; + } + + double CubicSpline::eval(double xx) const { + int klo = locate(xx); + int khi = klo + 1; + double h = x_[khi] - x_[klo]; + if (h <= 0.0) throw std::runtime_error("CubicSpline::eval: non-increasing x."); + double A = (x_[khi] - xx) / h; + double B = (xx - x_[klo]) / h; + double val = A * y_[klo] + B * y_[khi] + + ((A*A*A - A) * y2_[klo] + (B*B*B - B) * y2_[khi]) * (h*h) / 6.0; + return val; + } + + double CubicSpline::deriv(double xx) const { + int klo = locate(xx); + int khi = klo + 1; + double h = x_[khi] - x_[klo]; + if (h <= 0.0) throw std::runtime_error("CubicSpline::deriv: non-increasing x."); + double A = (x_[khi] - xx) / h; + double B = (xx - x_[klo]) / h; + double dy = (y_[khi] - y_[klo]) / h + + ( - (3.0*A*A - 1.0) * y2_[klo] + (3.0*B*B - 1.0) * y2_[khi] ) * (h / 6.0); + return dy; + } + + double CubicSpline::xmin() const { return x_.front(); } + double CubicSpline::xmax() const { return x_.back(); } + +} diff --git a/utils/Test/Deprojector.H b/utils/Test/Deprojector.H new file mode 100644 index 000000000..ba542d34f --- /dev/null +++ b/utils/Test/Deprojector.H @@ -0,0 +1,58 @@ +#ifndef DEPROJECTOR_H +#define DEPROJECTOR_H + +#include +#include + +#include "CubicSpline.H" + +namespace Deproject +{ + + class Deprojector { + public: + // Construct from analytic Sigma functor. If dSigmaFunc is empty, numerical derivative is used. + // R_data_min and R_data_max define where SigmaFunc is trusted for tail matching. + Deprojector(std::function SigmaFunc, + std::function dSigmaFunc, + double R_data_min, + double R_data_max, + double R_max_extend = -1.0, + double tail_power = -3.0, + int Ngrid = 4000); + + // Construct from sampled data arrays (R must be positive and will be sorted) + Deprojector(const std::vector& R_in, + const std::vector& Sigma_in, + double R_max_extend = -1.0, + double tail_power = -3.0, + int Ngrid = 4000); + + // Evaluate rho at a single radius + double rho_at(double r) const; + + // Evaluate rho for a vector of r + std::vector rho(const std::vector& r_eval) const; + + private: + void build_grid(int Ngrid); + + std::function sigma_func_; + std::function dsigma_func_; // may be empty + CubicSpline spline_; // used when constructed from sampled data + + double Rdata_min_, Rdata_max_; + double Rmin_, Rmax_; + double tail_power_; + + std::vector dx_; // grid spacings + + std::vector fineR_; + std::vector Sigma_f_; + std::vector dSigma_f_; + }; + +} + +#endif // DEPROJECTOR_H + diff --git a/utils/Test/Deprojector.cc b/utils/Test/Deprojector.cc new file mode 100644 index 000000000..d778a29b3 --- /dev/null +++ b/utils/Test/Deprojector.cc @@ -0,0 +1,200 @@ +#include +#include +#include + +#include "Deprojector.H" + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +namespace Deproject +{ + + Deprojector::Deprojector(std::function SigmaFunc, + std::function dSigmaFunc, + double R_data_min, + double R_data_max, + double R_max_extend, + double tail_power, + int Ngrid) + : sigma_func_(SigmaFunc), dsigma_func_(dSigmaFunc), + Rdata_min_(R_data_min), Rdata_max_(R_data_max), + tail_power_(tail_power) + { + if (Rdata_min_ <= 0.0 || Rdata_max_ <= Rdata_min_) + throw std::runtime_error("Invalid R_data_min/R_data_max."); + Rmin_ = Rdata_min_; + Rmax_ = (R_max_extend > Rdata_max_) ? R_max_extend : Rdata_max_; + build_grid(Ngrid); + } + + Deprojector::Deprojector(const std::vector& R_in, + const std::vector& Sigma_in, + double R_max_extend, + double tail_power, + int Ngrid) + : Rdata_min_(0.0), Rdata_max_(0.0), tail_power_(tail_power) + { + if (R_in.size() != Sigma_in.size() || R_in.size() < 2) throw std::runtime_error("Input R and Sigma must be same size and >=2."); + // copy & sort + std::vector> pairs; + pairs.reserve(R_in.size()); + for (size_t i=0;i R(R_in.size()), S(R_in.size()); + for (size_t i=0;i Rdata_max_) ? R_max_extend : Rdata_max_; + + sigma_func_ = [this](double rr){ return this->spline_.eval(rr); }; + dsigma_func_ = [this](double rr){ return this->spline_.deriv(rr); }; + + build_grid(Ngrid); + } + + + // --- updated build_grid --- + void Deprojector::build_grid(int Ngrid) { + if (Ngrid < 3) Ngrid = 3; + fineR_.resize(Ngrid); + for (int i = 0; i < Ngrid; ++i) { + double t = (double)i / (Ngrid - 1); + fineR_[i] = Rmin_ + t * (Rmax_ - Rmin_); + } + + // precompute spacing + dx_.resize(Ngrid - 1); + for (int i = 0; i < Ngrid - 1; ++i) dx_[i] = fineR_[i+1] - fineR_[i]; + + Sigma_f_.assign(Ngrid, 0.0); + dSigma_f_.assign(Ngrid, 0.0); + + bool have_dsf = static_cast(dsigma_func_); + + if (have_dsf) { + for (int i = 0; i < Ngrid; ++i) { + double rr = fineR_[i]; + if (rr <= Rdata_max_) { + Sigma_f_[i] = sigma_func_(rr); + dSigma_f_[i] = dsigma_func_(rr); + } else { + double Sig_at = sigma_func_(Rdata_max_); + double factor = std::pow(rr / Rdata_max_, tail_power_); + Sigma_f_[i] = Sig_at * factor; + if (rr > 0.0) + dSigma_f_[i] = Sig_at * tail_power_ * std::pow(rr, tail_power_ - 1.0) / std::pow(Rdata_max_, tail_power_); + else + dSigma_f_[i] = 0.0; + } + } + } else { + // compute Sigma on grid, then finite-difference derivative using neighbor spacing + for (int i = 0; i < Ngrid; ++i) { + double rr = fineR_[i]; + if (rr <= Rdata_max_) Sigma_f_[i] = sigma_func_(rr); + else { + double Sig_at = sigma_func_(Rdata_max_); + double factor = std::pow(rr / Rdata_max_, tail_power_); + Sigma_f_[i] = Sig_at * factor; + } + } + + for (int i = 0; i < Ngrid; ++i) { + if (i > 0 && i < Ngrid - 1) { + // centered difference using grid neighbors (robust) + double x1 = fineR_[i-1], x2 = fineR_[i+1]; + double y1 = Sigma_f_[i-1], y2 = Sigma_f_[i+1]; + dSigma_f_[i] = (y2 - y1) / (x2 - x1); + } else if (i == 0) { + // forward diff + double x1 = fineR_[1]; + dSigma_f_[i] = (Sigma_f_[1] - Sigma_f_[0]) / (x1 - fineR_[0]); + } else { // i == Ngrid-1 + double x1 = fineR_[Ngrid-2]; + dSigma_f_[i] = (Sigma_f_[Ngrid-1] - Sigma_f_[Ngrid-2]) / (fineR_[Ngrid-1] - x1); + } + } + } + } + + // --- updated rho_at --- + double Deprojector::rho_at(double r) const { + if (r >= Rmax_) return 0.0; + + // find index near r + // choose local offset delta = 0.5 * local grid spacing to avoid singularity + auto it0 = std::lower_bound(fineR_.begin(), fineR_.end(), r); + int idx0 = (int)std::distance(fineR_.begin(), it0); + // clamp idx0 + if (idx0 <= 0) idx0 = 1; + if (idx0 >= (int)dx_.size()) idx0 = (int)dx_.size() - 1; + + double local_dx = dx_[std::max(0, idx0 - 1)]; + double delta = 0.5 * local_dx; // half a local cell + double rstart = r + delta; + // ensure rstart >= fineR_[0] + if (rstart < fineR_[0]) rstart = fineR_[0]; + + // locate starting index after rstart + auto it = std::lower_bound(fineR_.begin(), fineR_.end(), rstart); + int idx = (int)std::distance(fineR_.begin(), it); + if (idx >= (int)fineR_.size()) return 0.0; + + double integral = 0.0; + int N = (int)fineR_.size(); + + // integrate using trapezoid on intervals [R_{i-1}, R_i] for all i such that R_{i-1} >= rstart or partial first + if (idx > 0) { + // partial segment from a = max(fineR_[idx-1], rstart) to b = fineR_[idx] + int i0 = idx - 1; + double R0 = fineR_[i0], R1 = fineR_[i0+1]; + double a = std::max(R0, rstart); + double b = R1; + if (b > a) { + // linear interpolation for derivative at 'a' + double t = (a - R0) / (R1 - R0); + double dSigma_a = dSigma_f_[i0] + t * (dSigma_f_[i0+1] - dSigma_f_[i0]); + double f_a = - dSigma_a / std::sqrt(std::max(1e-300, a*a - r*r)); + double f_b = - dSigma_f_[i0+1] / std::sqrt(std::max(1e-300, b*b - r*r)); + integral += 0.5 * (f_a + f_b) * (b - a); + } + // full segments from i = idx+1 .. N-1 + for (int i = idx + 1; i < N; ++i) { + double Rlo = fineR_[i-1], Rhi = fineR_[i]; + double f_lo = - dSigma_f_[i-1] / std::sqrt(std::max(1e-300, Rlo*Rlo - r*r)); + double f_hi = - dSigma_f_[i] / std::sqrt(std::max(1e-300, Rhi*Rhi - r*r)); + integral += 0.5 * (f_lo + f_hi) * (Rhi - Rlo); + } + } else { + // idx == 0 case: integrate over full segments starting at fineR_[0] + for (int i = 1; i < N; ++i) { + double Rlo = fineR_[i-1], Rhi = fineR_[i]; + double f_lo = - dSigma_f_[i-1] / std::sqrt(std::max(1e-300, Rlo*Rlo - r*r)); + double f_hi = - dSigma_f_[i] / std::sqrt(std::max(1e-300, Rhi*Rhi - r*r)); + integral += 0.5 * (f_lo + f_hi) * (Rhi - Rlo); + } + } + + return integral / M_PI; + } + + std::vector Deprojector::rho(const std::vector& r_eval) const { + std::vector out; + out.reserve(r_eval.size()); + for (double r : r_eval) out.push_back(rho_at(r)); + return out; + } + +} diff --git a/utils/Test/testDeprojPlummer.cc b/utils/Test/testDeprojPlummer.cc new file mode 100644 index 000000000..f72ecf63d --- /dev/null +++ b/utils/Test/testDeprojPlummer.cc @@ -0,0 +1,76 @@ +#include +#include +#include +#include +#include + +#include "Deprojector.H" + +using namespace Deproject; + +int main() +{ + std::function SigmaFunc, dSigmaFunc, RhoFunc; + + enum class Type { Plummer, Gaussian, Toomre }; + Type which = Type::Toomre; + + switch (which) { + case Type::Toomre: + // Test function + SigmaFunc = [](double R)->double + { return 1.0 / std::pow(1.0 + R*R, 1.5); }; + // Analytic derivative + dSigmaFunc = [](double R)->double + { return -3.0 * R / std::pow(1.0 + R*R, 2.5); }; + // Expected result + RhoFunc = [](double r)->double + { return 2.0 / std::pow(1.0 + r*r, 2.0) / M_PI; }; + break; + case Type::Gaussian: + // Test function + SigmaFunc = [](double R)->double + { return exp(-0.5*R*R); }; + // Analytic derivative + dSigmaFunc = [](double R)->double + { return -R*exp(-0.5*R*R); }; + // Expected result + RhoFunc = [](double r)->double + { return exp(-0.5*r*r)/sqrt(2.0*M_PI); }; + break; + default: + case Type::Plummer: + // Test function + SigmaFunc = [](double R)->double + { return 4.0 / 3.0 / std::pow(1.0 + R*R, 2.0); }; + // Analytic derivative + dSigmaFunc = [](double R)->double + { return -16.0 * R / 3.0 / std::pow(1.0 + R*R, 3.0); }; + // Expected result + RhoFunc = [](double r)->double + { return 1.0 / std::pow(1.0 + r*r, 2.5); }; + break; + } + + Deprojector D(SigmaFunc, dSigmaFunc, /*R_data_min=*/0.01, /*R_data_max=*/10.0, + /*R_max_extend=*/50.0, /*tail_power=*/-4.0, /*Ngrid=*/6000); + + std::vector r_eval; + int Nr = 150; + for (int i = 0; i < Nr; ++i) { + double t = (double)i / (Nr - 1); + r_eval.push_back(0.01 + t * 8.0); + } + auto rho = D.rho(r_eval); + + std::ofstream ofs("rho_test.txt"); + for (size_t i = 0; i < r_eval.size(); ++i) + ofs << std::setw(16) << r_eval[i] + << std::setw(16) << rho[i] + << std::setw(16) << RhoFunc(r_eval[i]) + << std::endl; + ofs.close(); + std::cout << "Wrote rho_test.txt\n"; + + return 0; +} diff --git a/utils/Test/testDeproject.cc b/utils/Test/testDeproject.cc new file mode 100644 index 000000000..c376adb6d --- /dev/null +++ b/utils/Test/testDeproject.cc @@ -0,0 +1,63 @@ +#include +#include +#include +#include + +#include "Deprojector.H" + +using namespace Deproject; + +int main() +{ + // Example A: construct from sampled data + { + std::vector Rdata, Sigma; + int Ndata = 2000; + double Rmin = 0.01, Rmax = 10.0; + for (int i = 0; i < Ndata; ++i) { + double t = (double)i / (Ndata - 1); + double r = Rmin + t * (Rmax - Rmin); + Rdata.push_back(r); + Sigma.push_back(1.0 / std::pow(1.0 + r*r, -1.5)); + } + + Deprojector D(Rdata, Sigma, /*R_max_extend=*/50.0, /*tail_power=*/-4.0, /*Ngrid=*/6000); + + std::vector r_eval; + int Nr = 150; + for (int i = 0; i < Nr; ++i) { + double t = (double)i / (Nr - 1); + r_eval.push_back(0.01 + t * 8.0); + } + auto rho = D.rho(r_eval); + + std::ofstream ofs("rho_from_sampled.txt"); + for (size_t i = 0; i < r_eval.size(); ++i) ofs << r_eval[i] << " " << rho[i] << "\n"; + ofs.close(); + std::cout << "Wrote rho_from_sampled.txt\n"; + } + + // Example B: construct from analytic functor + analytic derivative + { + auto SigmaFunc = [](double R)->double { return 1.0 / std::pow(1.0 + R*R, -1.5); }; + auto dSigmaFunc = [](double R)->double { return -3.0 * R / std::pow(1.0 + R*R, -2.5); }; // analytic derivative + + Deprojector D(SigmaFunc, dSigmaFunc, /*R_data_min=*/0.01, /*R_data_max=*/10.0, + /*R_max_extend=*/50.0, /*tail_power=*/-4.0, /*Ngrid=*/6000); + + std::vector r_eval; + int Nr = 150; + for (int i = 0; i < Nr; ++i) { + double t = (double)i / (Nr - 1); + r_eval.push_back(0.01 + t * 8.0); + } + auto rho = D.rho(r_eval); + + std::ofstream ofs("rho_from_functor.txt"); + for (size_t i = 0; i < r_eval.size(); ++i) ofs << r_eval[i] << " " << rho[i] << "\n"; + ofs.close(); + std::cout << "Wrote rho_from_functor.txt\n"; + } + + return 0; +} From 22adb601bc86f192539370402fc42febb12605b6 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 13 Mar 2026 12:09:46 -0400 Subject: [PATCH 04/78] Missing source driver file --- utils/Test/testEmpDeproj.cc | 167 ++++++++++++++++++++++++++++++++++++ 1 file changed, 167 insertions(+) create mode 100644 utils/Test/testEmpDeproj.cc diff --git a/utils/Test/testEmpDeproj.cc b/utils/Test/testEmpDeproj.cc new file mode 100644 index 000000000..6c344a2cd --- /dev/null +++ b/utils/Test/testEmpDeproj.cc @@ -0,0 +1,167 @@ +#include +#include +#include +#include +#include + +#include "Deprojector.H" +#include "EmpDeproj.H" +#include "cxxopts.H" + +using namespace Deproject; + +int main(int argc, char* argv[]) +{ + + // Parameters + // + std::string type, abel, fname; + double H, Rmin, Rmax, Rcut, Rwid; + int Nr, Ngrid, NumR, Nint; + + // Parse command-line options + // + cxxopts::Options options("testEmpDeproj", + "Test the EmpDeproj class against Deprojector" + "for various surface density profiles."); + options.add_options() + ("h,help", "Print help") + ("type", "Surface density type (plummer, gaussian, toomre)", cxxopts::value()->default_value("toomre")) + ("abel", "Abel inversion method (derivative, subtraction, ibp)", cxxopts::value()->default_value("derivative")) + ("H", "Scale height for empirical deprojection", cxxopts::value(H)->default_value("0.1")) + ("Nr", "Number of radial points to evaluate", cxxopts::value(Nr)->default_value("150")) + ("o,output", "Output file name", cxxopts::value(fname)->default_value("rho_test.txt")) + ("Rmin", "Minimum radius for evaluation", cxxopts::value(Rmin)->default_value("0.01")) + ("Rmax", "Maximum radius for evaluation", cxxopts::value(Rmax)->default_value("10.0")) + ("Rcut", "Inner cutoff for donut test", cxxopts::value(Rcut)->default_value("-1.0")) + ("Rwid", "Width of transition region to inner donut", cxxopts::value(Rwid)->default_value("0.2")) + ("Ngrid", "Number of grid points for Deprojector", cxxopts::value(Ngrid)->default_value("6000")) + ("NumR", "Number of radial points for EmpDeproj", cxxopts::value(NumR)->default_value("1000")) + ("Nint", "Number of integration points for EmpDeproj", cxxopts::value(Nint)->default_value("800")); + + auto result = options.parse(argc, argv); + + if (result.count("help")) { + std::cout << options.help() << std::endl; + return 0; + } + + // Map Abel type string to enum + // + EmpDeproj::AbelType type_enum; + std::map abel_type_map = { + {"derivative", EmpDeproj::AbelType::Derivative}, + {"subtraction", EmpDeproj::AbelType::Subtraction}, + {"ibp", EmpDeproj::AbelType::IBP} + }; + + // Convert type string to lower case + // + std::transform(abel.begin(), abel.end(), abel.begin(), + [](unsigned char c){ return std::tolower(c); }); + + auto it_abel = abel_type_map.find(result["abel"].as()); + + if (it_abel != abel_type_map.end()) { + type_enum = it_abel->second; + } else { + throw std::runtime_error("Unknown Abel type: " + result["abel"].as()); + } + + std::function SigmaFunc, dSigmaFunc, RhoFunc; + enum class Type { Plummer, Gaussian, Toomre }; + Type which = Type::Toomre; + + std::map type_map = { + {"plummer", Type::Plummer}, + {"gaussian", Type::Gaussian}, + {"toomre", Type::Toomre} + }; + + // Convert type string to lower case + // + std::transform(type.begin(), type.end(), type.begin(), + [](unsigned char c){ return std::tolower(c); }); + + auto it = type_map.find(result["type"].as()); + + if (it != type_map.end()) { + which = it->second; + } else { + throw std::runtime_error("Unknown type: " + result["type"].as()); + } + + switch (which) { + case Type::Toomre: + // Test function + SigmaFunc = [](double R)->double + { return 1.0 / std::pow(1.0 + R*R, 1.5); }; + // Analytic derivative + dSigmaFunc = [](double R)->double + { return -3.0 * R / std::pow(1.0 + R*R, 2.5); }; + // Expected result + RhoFunc = [](double r)->double + { return 2.0 / std::pow(1.0 + r*r, 2.0) / M_PI; }; + break; + case Type::Gaussian: + // Test function + SigmaFunc = [](double R)->double + { return exp(-0.5*R*R); }; + // Analytic derivative + dSigmaFunc = [](double R)->double + { return -R*exp(-0.5*R*R); }; + // Expected result + RhoFunc = [](double r)->double + { return exp(-0.5*r*r)/sqrt(2.0*M_PI); }; + break; + default: + case Type::Plummer: + // Test function + SigmaFunc = [](double R)->double + { return 4.0 / 3.0 / std::pow(1.0 + R*R, 2.0); }; + // Analytic derivative + dSigmaFunc = [](double R)->double + { return -16.0 * R / 3.0 / std::pow(1.0 + R*R, 3.0); }; + // Expected result + RhoFunc = [](double r)->double + { return 1.0 / std::pow(1.0 + r*r, 2.5); }; + break; + } + + Deprojector D(SigmaFunc, dSigmaFunc, /*R_data_min=*/0.01, /*R_data_max=*/10.0, + /*R_max_extend=*/50.0, /*tail_power=*/-4.0, /*Ngrid=*/6000); + + auto SigmaZFunc = [SigmaFunc, H, Rcut, Rwid](double R, double z)->double + { double Q = exp(-std::fabs(z)/(2.0*H)); + double sech = 2.0*Q / (1.0 + Q*Q); + double hole = 1.0; + if (Rcut > 0.0) { + double x = (R - Rcut) / Rwid; + hole = 0.5 * (1.0 + std::tanh(x)); + } + return SigmaFunc(R)*sech*sech*hole/(4.0*H); }; + + EmpDeproj E(H, Rmin, Rmax, NumR, Nint, SigmaZFunc, type_enum); + + std::vector r_eval; + for (int i = 0; i < Nr; ++i) { + double t = (double)i / (Nr - 1); + r_eval.push_back(0.01 + t * 8.0); + } + auto rho = D.rho(r_eval); + + std::ofstream ofs(fname); + for (size_t i = 0; i < r_eval.size(); ++i) + ofs << std::setw(16) << r_eval[i] + << std::setw(16) << rho[i] + << std::setw(16) << E.density(r_eval[i]) + << std::setw(16) << RhoFunc(r_eval[i]) + << std::setw(16) << SigmaFunc(r_eval[i]) + << std::setw(16) << E.surfaceDensity(r_eval[i]) + << std::endl; + + ofs.close(); + std::cout << "Wrote " << fname << std::endl; + + return 0; +} From b12a1f8f38ef448933bf15abad990cec2890e684 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 13 Mar 2026 12:10:52 -0400 Subject: [PATCH 05/78] Missing test class --- utils/Test/EmpDeproj.H | 40 +++++++++++++ utils/Test/EmpDeproj.cc | 123 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 163 insertions(+) create mode 100644 utils/Test/EmpDeproj.H create mode 100644 utils/Test/EmpDeproj.cc diff --git a/utils/Test/EmpDeproj.H b/utils/Test/EmpDeproj.H new file mode 100644 index 000000000..402f547cf --- /dev/null +++ b/utils/Test/EmpDeproj.H @@ -0,0 +1,40 @@ +#pragma once + +#include +#include "interp.H" + +class EmpDeproj +{ +private: + //! Interpolators for density and mass + Linear1d densRg, massRg, surfRg; + +public: + //! Abel type + enum class AbelType { Derivative, Subtraction, IBP }; + + //! Construct from analytic Sigma functor + EmpDeproj(double H, double Rmin, double Rmax, int NUMR, int NINT, + std::function func, + AbelType type = AbelType::Derivative); + + //! Destructor + ~EmpDeproj() {} + + //! Evaluate deprojected density at a single radius + double density(double R) + { + return densRg.eval(log(R)); + } + + double surfaceDensity(double R) + { + return surfRg.eval(log(R)); + } + + //! Evaluate deprojected mass at a single radius + double mass(double R) + { + return massRg.eval(log(R)); + } +}; diff --git a/utils/Test/EmpDeproj.cc b/utils/Test/EmpDeproj.cc new file mode 100644 index 000000000..ceea6a7d9 --- /dev/null +++ b/utils/Test/EmpDeproj.cc @@ -0,0 +1,123 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "EmpDeproj.H" +#include "gaussQ.H" + +// EmpCylSL: Empirical cylindrical deprojection by numerical +// integration and finite difference + +EmpDeproj::EmpDeproj(double H, double RMIN, double RMAX, int NUMR, int NINT, + std::function func, AbelType type) +{ + LegeQuad lq(NINT); + + std::vector rr(NUMR), rl(NUMR), sigI(NUMR), rhoI(NUMR, 0.0); + + double Rmin = log(RMIN); + double Rmax = log(RMAX); + + double dr = (Rmax - Rmin)/(NUMR-1); + + // Compute surface mass density, Sigma(R) + // + for (int i=0; i rho(NUMR), mass(NUMR); + for (int i=0; i Date: Fri, 13 Mar 2026 12:42:38 -0400 Subject: [PATCH 06/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- include/DiskModels.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/DiskModels.H b/include/DiskModels.H index 273638d4f..ec18dc6b4 100644 --- a/include/DiskModels.H +++ b/include/DiskModels.H @@ -217,7 +217,7 @@ public: double operator()(double R, double z, double phi=0.) { - double sigma = std::pow(1.0 * (R/a)*(R/a), -0.5*n) * norm; + double sigma = norm * std::pow(1.0 + (R/a)*(R/a), -0.5*n); double Z = std::fabs(z); double Q = std::exp(-Z/h); double sech = 2.0*Q/(1.0 + Q*Q); From a03bbe692698f2513ac37c43d15651568a48ae88 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Fri, 13 Mar 2026 12:43:11 -0400 Subject: [PATCH 07/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- include/DiskModels.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/DiskModels.H b/include/DiskModels.H index ec18dc6b4..f7b631503 100644 --- a/include/DiskModels.H +++ b/include/DiskModels.H @@ -212,7 +212,7 @@ public: params.push_back(h); params.push_back(n); if (n <= 2) throw std::runtime_error("Toomre index must be > 2"); - norm = (n - 2.0)/(2.0*M_PI*a*a); + norm = M*(n - 2.0)/(2.0*M_PI*a*a); } double operator()(double R, double z, double phi=0.) From 40425440d7d1b348aee565965de30c2c103d1352 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Fri, 13 Mar 2026 12:43:34 -0400 Subject: [PATCH 08/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- include/DiskModels.H | 1 + 1 file changed, 1 insertion(+) diff --git a/include/DiskModels.H b/include/DiskModels.H index f7b631503..54edc17bb 100644 --- a/include/DiskModels.H +++ b/include/DiskModels.H @@ -3,6 +3,7 @@ #include "EmpCylSL.H" #include "DiskDensityFunc.H" +#include //! A Ferrers Ellipsoid + Evacuated Exponential Disc (semi-realistic //! bar+disk model) From db0e02d07fadf94006d41588492c2c3ae9425949 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Fri, 13 Mar 2026 12:45:35 -0400 Subject: [PATCH 09/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/Test/testEmpDeproj.cc | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/utils/Test/testEmpDeproj.cc b/utils/Test/testEmpDeproj.cc index 6c344a2cd..a3f75d11d 100644 --- a/utils/Test/testEmpDeproj.cc +++ b/utils/Test/testEmpDeproj.cc @@ -3,6 +3,11 @@ #include #include #include +#include +#include +#include +#include +#include #include "Deprojector.H" #include "EmpDeproj.H" From dfc99f9d2f8b7ac66dc109ed23543a4cf5fa17f0 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Fri, 13 Mar 2026 12:46:13 -0400 Subject: [PATCH 10/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/Test/testDeprojPlummer.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/Test/testDeprojPlummer.cc b/utils/Test/testDeprojPlummer.cc index f72ecf63d..e4fb296af 100644 --- a/utils/Test/testDeprojPlummer.cc +++ b/utils/Test/testDeprojPlummer.cc @@ -13,7 +13,7 @@ int main() std::function SigmaFunc, dSigmaFunc, RhoFunc; enum class Type { Plummer, Gaussian, Toomre }; - Type which = Type::Toomre; + Type which = Type::Plummer; switch (which) { case Type::Toomre: From 2e6fee4adf8b9ceb21b7a814773582cddcb7a804 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Fri, 13 Mar 2026 12:49:09 -0400 Subject: [PATCH 11/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/Test/EmpDeproj.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/utils/Test/EmpDeproj.cc b/utils/Test/EmpDeproj.cc index ceea6a7d9..cfef407f1 100644 --- a/utils/Test/EmpDeproj.cc +++ b/utils/Test/EmpDeproj.cc @@ -7,6 +7,7 @@ #include #include #include +#include #include "EmpDeproj.H" #include "gaussQ.H" From 204d43d9b6f5493f44b8c4d88668b7d6325053c3 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 13 Mar 2026 12:56:28 -0400 Subject: [PATCH 12/78] Add a test routine from EmpDeproj alone; fix parsing error in testEmpDeproj for method variables --- utils/Test/CMakeLists.txt | 3 +- utils/Test/testEmp.cc | 145 ++++++++++++++++++++++++++++++++++++ utils/Test/testEmpDeproj.cc | 8 +- 3 files changed, 151 insertions(+), 5 deletions(-) create mode 100644 utils/Test/testEmp.cc diff --git a/utils/Test/CMakeLists.txt b/utils/Test/CMakeLists.txt index 78d7e4952..1cd889a1f 100644 --- a/utils/Test/CMakeLists.txt +++ b/utils/Test/CMakeLists.txt @@ -1,6 +1,6 @@ set(bin_PROGRAMS testBarrier expyaml orthoTest testDeproj - testDeprojPlum testEmpDeproj) + testDeprojPlum testEmpDeproj testEmp) set(common_LINKLIB OpenMP::OpenMP_CXX MPI::MPI_CXX expui exputil yaml-cpp ${VTK_LIBRARIES}) @@ -38,6 +38,7 @@ add_executable(testDeprojPlum testDeprojPlummer.cc CubicSpline.cc Deprojector.cc) add_executable(testEmpDeproj testEmpDeproj.cc CubicSpline.cc Deprojector.cc EmpDeproj.cc) +add_executable(testEmp testEmp.cc EmpDeproj.cc) foreach(program ${bin_PROGRAMS}) target_link_libraries(${program} ${common_LINKLIB}) diff --git a/utils/Test/testEmp.cc b/utils/Test/testEmp.cc new file mode 100644 index 000000000..796d016d9 --- /dev/null +++ b/utils/Test/testEmp.cc @@ -0,0 +1,145 @@ +#include +#include +#include +#include +#include + +#include "Deprojector.H" +#include "EmpDeproj.H" +#include "cxxopts.H" + +using namespace Deproject; + +int main(int argc, char* argv[]) +{ + + // Parameters + // + std::string type, abel, fname; + double H, Rmin, Rmax, Rcut, Rwid; + int Nr, NumR, Nint; + + // Parse command-line options + // + cxxopts::Options options("testDonut", + "Test the EmpDeproj class for an inner donut-shaped " + "density distribution, using the Toomre profile as " + "a test case."); + + options.add_options() + ("h,help", "Print help") + ("type", "Surface density type (plummer, gaussian, toomre)", cxxopts::value()->default_value("toomre")) + ("abel", "Abel inversion method (derivative, subtraction, ibp)", cxxopts::value()->default_value("derivative")) + ("H", "Scale height for empirical deprojection", cxxopts::value(H)->default_value("0.1")) + ("Nr", "Number of radial points to evaluate", cxxopts::value(Nr)->default_value("150")) + ("o,output", "Output file name", cxxopts::value(fname)->default_value("rho_test.txt")) + ("Rmin", "Minimum radius for evaluation", cxxopts::value(Rmin)->default_value("0.01")) + ("Rmax", "Maximum radius for evaluation", cxxopts::value(Rmax)->default_value("10.0")) + ("Rcut", "Inner cutoff for donut test", cxxopts::value(Rcut)->default_value("-1.0")) + ("Rwid", "Width of transition region to inner donut", cxxopts::value(Rwid)->default_value("0.2")) + ("NumR", "Number of radial points for EmpDeproj", cxxopts::value(NumR)->default_value("1000")) + ("Nint", "Number of integration points for EmpDeproj", cxxopts::value(Nint)->default_value("800")); + + auto result = options.parse(argc, argv); + + if (result.count("help")) { + std::cout << options.help() << std::endl; + return 0; + } + + // Map Abel type string to enum + // + EmpDeproj::AbelType type_enum; + std::map abel_type_map = { + {"derivative", EmpDeproj::AbelType::Derivative}, + {"subtraction", EmpDeproj::AbelType::Subtraction}, + {"ibp", EmpDeproj::AbelType::IBP} + }; + + // Convert type string to lower case + // + std::transform(abel.begin(), abel.end(), abel.begin(), + [](unsigned char c){ return std::tolower(c); }); + + auto it_abel = abel_type_map.find(result["abel"].as()); + + if (it_abel != abel_type_map.end()) { + type_enum = it_abel->second; + } else { + throw std::runtime_error("Unknown Abel type: " + result["abel"].as()); + } + + std::function SigmaFunc, dSigmaFunc, RhoFunc; + enum class Type { Plummer, Gaussian, Toomre }; + Type which = Type::Toomre; + + std::map type_map = { + {"plummer", Type::Plummer}, + {"gaussian", Type::Gaussian}, + {"toomre", Type::Toomre} + }; + + // Convert type string to lower case + // + std::transform(type.begin(), type.end(), type.begin(), + [](unsigned char c){ return std::tolower(c); }); + + auto it = type_map.find(result["type"].as()); + + if (it != type_map.end()) { + which = it->second; + } else { + throw std::runtime_error("Unknown type: " + result["type"].as()); + } + + switch (which) { + case Type::Toomre: + // Test function + SigmaFunc = [](double R)->double + { return 1.0 / std::pow(1.0 + R*R, 1.5); }; + // Analytic derivative + break; + case Type::Gaussian: + // Test function + SigmaFunc = [](double R)->double + { return exp(-0.5*R*R); }; + // Analytic derivative + break; + default: + case Type::Plummer: + // Test function + SigmaFunc = [](double R)->double + { return 4.0 / 3.0 / std::pow(1.0 + R*R, 2.0); }; + break; + } + + auto SigmaZFunc = [SigmaFunc, H, Rcut, Rwid](double R, double z)->double + { double Q = exp(-std::fabs(z)/(2.0*H)); + double sech = 2.0*Q / (1.0 + Q*Q); + double hole = 1.0; + if (Rcut > 0.0) { + double x = (R - Rcut) / Rwid; + hole = 0.5 * (1.0 + std::tanh(x)); + } + return SigmaFunc(R)*sech*sech*hole/(4.0*H); }; + + EmpDeproj E(H, Rmin, Rmax, NumR, Nint, SigmaZFunc, type_enum); + + std::vector r_eval; + for (int i = 0; i < Nr; ++i) { + double t = (double)i / (Nr - 1); + r_eval.push_back(0.01 + t * 8.0); + } + + std::ofstream ofs(fname); + for (size_t i = 0; i < r_eval.size(); ++i) + ofs << std::setw(16) << r_eval[i] + << std::setw(16) << E.density(r_eval[i]) + << std::setw(16) << E.surfaceDensity(r_eval[i]) + << std::endl; + + ofs.close(); + std::cout << "Wrote " << fname << std::endl; + + return 0; +} diff --git a/utils/Test/testEmpDeproj.cc b/utils/Test/testEmpDeproj.cc index a3f75d11d..966313028 100644 --- a/utils/Test/testEmpDeproj.cc +++ b/utils/Test/testEmpDeproj.cc @@ -31,8 +31,8 @@ int main(int argc, char* argv[]) "for various surface density profiles."); options.add_options() ("h,help", "Print help") - ("type", "Surface density type (plummer, gaussian, toomre)", cxxopts::value()->default_value("toomre")) - ("abel", "Abel inversion method (derivative, subtraction, ibp)", cxxopts::value()->default_value("derivative")) + ("type", "Surface density type (plummer, gaussian, toomre)", cxxopts::value(type)->default_value("toomre")) + ("abel", "Abel inversion method (derivative, subtraction, ibp)", cxxopts::value(abel)->default_value("derivative")) ("H", "Scale height for empirical deprojection", cxxopts::value(H)->default_value("0.1")) ("Nr", "Number of radial points to evaluate", cxxopts::value(Nr)->default_value("150")) ("o,output", "Output file name", cxxopts::value(fname)->default_value("rho_test.txt")) @@ -65,7 +65,7 @@ int main(int argc, char* argv[]) std::transform(abel.begin(), abel.end(), abel.begin(), [](unsigned char c){ return std::tolower(c); }); - auto it_abel = abel_type_map.find(result["abel"].as()); + auto it_abel = abel_type_map.find(abel); if (it_abel != abel_type_map.end()) { type_enum = it_abel->second; @@ -88,7 +88,7 @@ int main(int argc, char* argv[]) std::transform(type.begin(), type.end(), type.begin(), [](unsigned char c){ return std::tolower(c); }); - auto it = type_map.find(result["type"].as()); + auto it = type_map.find(type); if (it != type_map.end()) { which = it->second; From 0e82cb6d409a3f11a1a7dc8e5586f35b2be41b7a Mon Sep 17 00:00:00 2001 From: mdw Date: Fri, 13 Mar 2026 12:58:44 -0400 Subject: [PATCH 13/78] Updated parsing in BiorthBasis::Cylindrical for deprojection --- expui/BiorthBasis.cc | 2 +- include/DiskModels.H | 11 +++++------ 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index fda83a3ef..386d3b136 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1735,7 +1735,7 @@ namespace BasisClasses if (PTYPE == DeprojType::mn) // Miyamoto-Nagai model = std::make_shared(1.0, H); else if (PTYPE == DeprojType::toomre) { - model = std::make_shared(1.0, H); + model = std::make_shared(1.0, H, 5.0); } else if (PTYPE == DeprojType::python and DTYPE == DiskType::python) { model = std::make_shared(pyname, acyl); diff --git a/include/DiskModels.H b/include/DiskModels.H index 54edc17bb..a951a2267 100644 --- a/include/DiskModels.H +++ b/include/DiskModels.H @@ -213,16 +213,15 @@ public: params.push_back(h); params.push_back(n); if (n <= 2) throw std::runtime_error("Toomre index must be > 2"); - norm = M*(n - 2.0)/(2.0*M_PI*a*a); + norm = M*(n - 2.0)/(2.0*M_PI*a*a*4.0*h); } double operator()(double R, double z, double phi=0.) { - double sigma = norm * std::pow(1.0 + (R/a)*(R/a), -0.5*n); - double Z = std::fabs(z); - double Q = std::exp(-Z/h); - double sech = 2.0*Q/(1.0 + Q*Q); - return sigma * sech*sech * 0.25/h; + double sigma = std::pow(1.0 + (R/a)*(R/a), -0.5*n) * norm; + double Q = std::exp(-std::fabs(z)/h); + double sech = 2.0*Q/(1.0 + Q*Q); // Prevent overflow + return sigma * sech*sech; } }; From cb5e150f5060d613b9d465db61c99d7df80777bc Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Fri, 13 Mar 2026 14:32:22 -0400 Subject: [PATCH 14/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 386d3b136..d6cba31a6 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1711,6 +1711,10 @@ namespace BasisClasses std::transform(dmodel.begin(), dmodel.end(), dmodel.begin(), [](unsigned char c){ return std::tolower(c); }); + // Map legacy/short model names to canonical keys expected by dplookup + if (dmodel == "exp") { + dmodel = "exponential"; + } // Check for map entry try { From 5f151bc0a2c385276d08ee75667301a53c14f3eb Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Fri, 13 Mar 2026 14:32:45 -0400 Subject: [PATCH 15/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index d6cba31a6..819ce37bb 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1728,7 +1728,7 @@ namespace BasisClasses } catch (const std::out_of_range& err) { if (myid==0) { - std::cout << "DeprojType error in configuraton file" << std::endl; + std::cout << "DeprojType error in configuration file" << std::endl; std::cout << "Valid options are: "; for (auto v : dplookup) std::cout << v.first << " "; std::cout << std::endl; From b7610507d79b2fbe7ef6c702802c5d844e0c950f Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Fri, 13 Mar 2026 14:39:35 -0400 Subject: [PATCH 16/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/Test/Deprojector.cc | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/utils/Test/Deprojector.cc b/utils/Test/Deprojector.cc index d778a29b3..10288b8b6 100644 --- a/utils/Test/Deprojector.cc +++ b/utils/Test/Deprojector.cc @@ -46,10 +46,14 @@ namespace Deproject for (size_t i=0;i Date: Fri, 13 Mar 2026 14:40:20 -0400 Subject: [PATCH 17/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/Test/testDeproject.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/utils/Test/testDeproject.cc b/utils/Test/testDeproject.cc index c376adb6d..309c9cf88 100644 --- a/utils/Test/testDeproject.cc +++ b/utils/Test/testDeproject.cc @@ -18,7 +18,7 @@ int main() double t = (double)i / (Ndata - 1); double r = Rmin + t * (Rmax - Rmin); Rdata.push_back(r); - Sigma.push_back(1.0 / std::pow(1.0 + r*r, -1.5)); + Sigma.push_back(std::pow(1.0 + r*r, -1.5)); } Deprojector D(Rdata, Sigma, /*R_max_extend=*/50.0, /*tail_power=*/-4.0, /*Ngrid=*/6000); @@ -39,7 +39,7 @@ int main() // Example B: construct from analytic functor + analytic derivative { - auto SigmaFunc = [](double R)->double { return 1.0 / std::pow(1.0 + R*R, -1.5); }; + auto SigmaFunc = [](double R)->double { return std::pow(1.0 + R*R, -1.5); }; auto dSigmaFunc = [](double R)->double { return -3.0 * R / std::pow(1.0 + R*R, -2.5); }; // analytic derivative Deprojector D(SigmaFunc, dSigmaFunc, /*R_data_min=*/0.01, /*R_data_max=*/10.0, From f9a6ae297f5ec0b73db46412817cc3bed728ab6a Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Fri, 13 Mar 2026 14:40:34 -0400 Subject: [PATCH 18/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/Test/testEmp.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/utils/Test/testEmp.cc b/utils/Test/testEmp.cc index 796d016d9..c2af16a16 100644 --- a/utils/Test/testEmp.cc +++ b/utils/Test/testEmp.cc @@ -3,6 +3,7 @@ #include #include #include +#include #include "Deprojector.H" #include "EmpDeproj.H" From 5a66623825241789c16e4c485769fcf103aadd7a Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Fri, 13 Mar 2026 14:41:16 -0400 Subject: [PATCH 19/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/Test/EmpDeproj.cc | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/utils/Test/EmpDeproj.cc b/utils/Test/EmpDeproj.cc index cfef407f1..b43962a5e 100644 --- a/utils/Test/EmpDeproj.cc +++ b/utils/Test/EmpDeproj.cc @@ -104,18 +104,20 @@ EmpDeproj::EmpDeproj(double H, double RMIN, double RMAX, int NUMR, int NINT, // Debug // - if (true) { +#ifdef EMPDEPROJ_DEBUG + { std::string fname("deproject_sl.txt"); std::ofstream out(fname); if (out) { for (int i=0; i Date: Fri, 13 Mar 2026 14:42:09 -0400 Subject: [PATCH 20/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/Test/testDeprojPlummer.cc | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/utils/Test/testDeprojPlummer.cc b/utils/Test/testDeprojPlummer.cc index e4fb296af..0e34cff9e 100644 --- a/utils/Test/testDeprojPlummer.cc +++ b/utils/Test/testDeprojPlummer.cc @@ -64,12 +64,13 @@ int main() auto rho = D.rho(r_eval); std::ofstream ofs("rho_test.txt"); - for (size_t i = 0; i < r_eval.size(); ++i) + for (size_t i = 0; i < r_eval.size(); ++i) { ofs << std::setw(16) << r_eval[i] << std::setw(16) << rho[i] << std::setw(16) << RhoFunc(r_eval[i]) << std::endl; - ofs.close(); + } + ofs.close(); std::cout << "Wrote rho_test.txt\n"; return 0; From 31dad722f0ce52ae8041a355b4f374020cfc9976 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Fri, 13 Mar 2026 14:42:45 -0400 Subject: [PATCH 21/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 819ce37bb..92c3343af 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1503,7 +1503,7 @@ namespace BasisClasses if (conf["cmapz" ]) cmapZ = conf["cmapz" ].as(); if (conf["ignore" ]) Ignore = conf["ignore" ].as(); if (conf["deproject" ]) deproject = conf["deproject" ].as(); - if (conf["dmodel" ]) dmodel = conf["dmodel" ].as(); + if (conf["dmodel" ]) dmodel = conf["dmodel" ].as(); if (conf["aratio" ]) aratio = conf["aratio" ].as(); if (conf["hratio" ]) hratio = conf["hratio" ].as(); From 9c2f36f96518986f7fac46d0e7f21b1b18dd2345 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Fri, 13 Mar 2026 14:43:19 -0400 Subject: [PATCH 22/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 92c3343af..c7a679a6d 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1198,7 +1198,7 @@ namespace BasisClasses {"python", DiskType::python} }; - // Dprojection model for cylindrical basis construction + // Deprojection model for cylindrical basis construction const std::map Cylindrical::dplookup = { {"mn", DeprojType::mn}, {"exponential", DeprojType::exponential}, From 22b6493b22adfe096db0576f1435a1bbf6747b8e Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 13 Mar 2026 15:08:58 -0400 Subject: [PATCH 23/78] Use new sech^2 scaling --- include/DiskModels.H | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/DiskModels.H b/include/DiskModels.H index 54edc17bb..8272e4b03 100644 --- a/include/DiskModels.H +++ b/include/DiskModels.H @@ -213,16 +213,16 @@ public: params.push_back(h); params.push_back(n); if (n <= 2) throw std::runtime_error("Toomre index must be > 2"); - norm = M*(n - 2.0)/(2.0*M_PI*a*a); + norm = M*(n - 2.0)/(2.0*M_PI*a*a*4.0*h); } double operator()(double R, double z, double phi=0.) { double sigma = norm * std::pow(1.0 + (R/a)*(R/a), -0.5*n); double Z = std::fabs(z); - double Q = std::exp(-Z/h); + double Q = std::exp(-0.5*Z/h); double sech = 2.0*Q/(1.0 + Q*Q); - return sigma * sech*sech * 0.25/h; + return sigma * sech*sech; } }; From 9439daa5ce94885e9ec4f6290310291800e7a58c Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Fri, 13 Mar 2026 19:39:09 -0400 Subject: [PATCH 24/78] Python deprojection test routine; add possibility of a separate Python functor for deprojection in EmpCylSL --- expui/BiorthBasis.H | 2 +- expui/BiorthBasis.cc | 11 +- utils/Test/CMakeLists.txt | 3 + utils/Test/testEmp.cc | 306 ++++++++++++++++++++++++++++---------- 4 files changed, 241 insertions(+), 81 deletions(-) diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index 054e28a35..d52341ab6 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -996,7 +996,7 @@ namespace BasisClasses //! Basis construction parameters double aratio, hratio, dweight, rwidth, ashift, rfactor, rtrunc, ppow, Mfac, HERNA; bool Ignore, deproject; - std::string pyname; + std::string pyname, pyproj; //! DiskType support //! diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index c7a679a6d..24608fb21 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1270,7 +1270,8 @@ namespace BasisClasses "playback", "coefCompute", "coefMaster", - "pyname" + "pyname", + "pyproj" }; Cylindrical::Cylindrical(const YAML::Node& CONF) : @@ -1519,6 +1520,7 @@ namespace BasisClasses if (conf["dtype" ]) dtype = conf["dtype" ].as(); if (conf["vflag" ]) vflag = conf["vflag" ].as(); if (conf["pyname" ]) pyname = conf["pyname" ].as(); + if (conf["pyproj" ]) pyname = conf["pyproj" ].as(); if (conf["pcavar"] ) pcavar = conf["pcavar" ].as(); if (conf["subsamp"] ) sampT = conf["subsamp" ].as(); @@ -1712,11 +1714,13 @@ namespace BasisClasses [](unsigned char c){ return std::tolower(c); }); // Map legacy/short model names to canonical keys expected by dplookup + // if (dmodel == "exp") { dmodel = "exponential"; } // Check for map entry + // try { PTYPE = dplookup.at(dmodel); @@ -1740,9 +1744,8 @@ namespace BasisClasses model = std::make_shared(1.0, H); else if (PTYPE == DeprojType::toomre) { model = std::make_shared(1.0, H, 5.0); - } else if (PTYPE == DeprojType::python and - DTYPE == DiskType::python) { - model = std::make_shared(pyname, acyl); + } else if (PTYPE == DeprojType::python) { + model = std::make_shared(pyproj, acyl); std::cout << "Using AxiSymPyModel for deprojection from Python function <" << pyname << ">" << std::endl; } else { // Default to exponential diff --git a/utils/Test/CMakeLists.txt b/utils/Test/CMakeLists.txt index 1cd889a1f..45531c968 100644 --- a/utils/Test/CMakeLists.txt +++ b/utils/Test/CMakeLists.txt @@ -47,4 +47,7 @@ foreach(program ${bin_PROGRAMS}) # install(TARGETS ${program} DESTINATION bin) endforeach() +# For Python interpreter... +target_link_libraries(testEmp ${common_LINKLIB} pybind11::embed) + install(TARGETS expyaml DESTINATION bin) diff --git a/utils/Test/testEmp.cc b/utils/Test/testEmp.cc index c2af16a16..9ee505732 100644 --- a/utils/Test/testEmp.cc +++ b/utils/Test/testEmp.cc @@ -1,36 +1,43 @@ #include #include #include +#include +#include +#include #include +#include #include -#include +#include -#include "Deprojector.H" #include "EmpDeproj.H" #include "cxxopts.H" -using namespace Deproject; +// pybind11 embedding +#include +#include +namespace py = pybind11; int main(int argc, char* argv[]) { - - // Parameters - // - std::string type, abel, fname; - double H, Rmin, Rmax, Rcut, Rwid; - int Nr, NumR, Nint; + // Default parameters + std::string type_opt; + std::string abel_opt; + std::string fname; + double H = 0.1, Rmin = 0.01, Rmax = 10.0, Rcut = -1.0, Rwid = 0.2; + int Nr = 150, NumR = 1000, Nint = 800; // Parse command-line options - // cxxopts::Options options("testDonut", - "Test the EmpDeproj class for an inner donut-shaped " - "density distribution, using the Toomre profile as " - "a test case."); + "Test the EmpDeproj class for an inner donut-shaped " + "density distribution, using the Toomre profile as " + "a test case."); options.add_options() ("h,help", "Print help") - ("type", "Surface density type (plummer, gaussian, toomre)", cxxopts::value()->default_value("toomre")) - ("abel", "Abel inversion method (derivative, subtraction, ibp)", cxxopts::value()->default_value("derivative")) + ("type", "Surface density type (plummer, gaussian, toomre) - ignored if Python function supplied", + cxxopts::value(type_opt)->default_value("toomre")) + ("abel", "Abel inversion method (derivative, subtraction, ibp)", + cxxopts::value(abel_opt)->default_value("derivative")) ("H", "Scale height for empirical deprojection", cxxopts::value(H)->default_value("0.1")) ("Nr", "Number of radial points to evaluate", cxxopts::value(Nr)->default_value("150")) ("o,output", "Output file name", cxxopts::value(fname)->default_value("rho_test.txt")) @@ -39,7 +46,10 @@ int main(int argc, char* argv[]) ("Rcut", "Inner cutoff for donut test", cxxopts::value(Rcut)->default_value("-1.0")) ("Rwid", "Width of transition region to inner donut", cxxopts::value(Rwid)->default_value("0.2")) ("NumR", "Number of radial points for EmpDeproj", cxxopts::value(NumR)->default_value("1000")) - ("Nint", "Number of integration points for EmpDeproj", cxxopts::value(Nint)->default_value("800")); + ("Nint", "Number of integration points for EmpDeproj", cxxopts::value(Nint)->default_value("800")) + // Python integration options + ("pymodule", "Python module name OR path to a .py file containing a function", cxxopts::value()) + ("pyfunc", "Function name inside module/file (default: 'Sigma')", cxxopts::value()->default_value("Sigma")); auto result = options.parse(argc, argv); @@ -49,7 +59,6 @@ int main(int argc, char* argv[]) } // Map Abel type string to enum - // EmpDeproj::AbelType type_enum; std::map abel_type_map = { {"derivative", EmpDeproj::AbelType::Derivative}, @@ -57,90 +66,235 @@ int main(int argc, char* argv[]) {"ibp", EmpDeproj::AbelType::IBP} }; - // Convert type string to lower case - // - std::transform(abel.begin(), abel.end(), abel.begin(), - [](unsigned char c){ return std::tolower(c); }); - - auto it_abel = abel_type_map.find(result["abel"].as()); + std::string abel_l = result["abel"].as(); + std::transform(abel_l.begin(), abel_l.end(), abel_l.begin(), + [](unsigned char c){ return std::tolower(c); }); + auto it_abel = abel_type_map.find(abel_l); if (it_abel != abel_type_map.end()) { type_enum = it_abel->second; } else { throw std::runtime_error("Unknown Abel type: " + result["abel"].as()); } - std::function SigmaFunc, dSigmaFunc, RhoFunc; - enum class Type { Plummer, Gaussian, Toomre }; - Type which = Type::Toomre; + // Prepare SigmaZFunc to pass into EmpDeproj + std::function SigmaZFunc; - std::map type_map = { - {"plummer", Type::Plummer}, - {"gaussian", Type::Gaussian}, - {"toomre", Type::Toomre} - }; + // For pybind11 embedding, we need to keep the interpreter alive for + // the duration of the SigmaZFunc usage. We can use a unique_ptr to + // manage the lifetime of the interpreter guard. If no Python is + // used, this will just be an empty guard that does nothing. + // + std::unique_ptr pyguard; - // Convert type string to lower case + // If user supplied a Python module/file and function, embed Python + // and load it here // - std::transform(type.begin(), type.end(), type.begin(), - [](unsigned char c){ return std::tolower(c); }); + if (result.count("pymodule")) { + std::string pymod = result["pymodule"].as(); + std::string pyfuncname = result["pyfunc"].as(); - auto it = type_map.find(result["type"].as()); + // Start the Python interpreter + pyguard = std::make_unique(); - if (it != type_map.end()) { - which = it->second; - } else { - throw std::runtime_error("Unknown type: " + result["type"].as()); - } + py::object py_module; - switch (which) { - case Type::Toomre: - // Test function - SigmaFunc = [](double R)->double - { return 1.0 / std::pow(1.0 + R*R, 1.5); }; - // Analytic derivative - break; - case Type::Gaussian: - // Test function - SigmaFunc = [](double R)->double - { return exp(-0.5*R*R); }; - // Analytic derivative - break; - default: - case Type::Plummer: - // Test function - SigmaFunc = [](double R)->double - { return 4.0 / 3.0 / std::pow(1.0 + R*R, 2.0); }; - break; - } - - auto SigmaZFunc = [SigmaFunc, H, Rcut, Rwid](double R, double z)->double - { double Q = exp(-std::fabs(z)/(2.0*H)); - double sech = 2.0*Q / (1.0 + Q*Q); - double hole = 1.0; - if (Rcut > 0.0) { - double x = (R - Rcut) / Rwid; - hole = 0.5 * (1.0 + std::tanh(x)); + // If pymod ends with .py, treat it as filepath: insert its + // directory to sys.path and import by stem + std::filesystem::path p(pymod); + try { + if (p.has_extension() && p.extension() == ".py") { + std::string dir = p.parent_path().string(); + std::string modname = p.stem().string(); + if (!dir.empty()) { + py::module_ sys = py::module_::import("sys"); + // insert at front so local dir is found first + sys.attr("path").attr("insert")(0, dir); + } + py_module = py::module_::import(modname.c_str()); + } else { + // treat as module name + py_module = py::module_::import(pymod.c_str()); + } + } catch (const py::error_already_set &e) { + throw std::runtime_error(std::string("Failed to import Python module '") + + pymod + "': " + e.what()); + } + + py::object pyfunc = py_module.attr(pyfuncname.c_str()); + + if (!py::isinstance(pyfunc) && !py::hasattr(pyfunc, "__call__")) { + throw std::runtime_error("Python object " + pyfuncname + + " is not callable."); + } + + // Inspect function argument count to decide whether it's Sigma(R) + // or SigmaZ(R,z). This is probably overkill and might not be + // robust for all callables, but it allows some flexibility for + // users. + // + int argcount = 0; + try { + // functions have __code__.co_argcount; builtins may not — in + // that case prefer calling and checking. I'm still not sure if + // this is the best way to do it, but it should cover most cases + // (functions, builtins, callables). Python experts? + // + if (py::hasattr(pyfunc, "__code__")) { + argcount = pyfunc.attr("__code__").attr("co_argcount").cast(); + } else if (py::hasattr(pyfunc, "__call__") && py::hasattr(pyfunc.attr("__call__"), "__code__")) { + argcount = pyfunc.attr("__call__").attr("__code__").attr("co_argcount").cast() - 1; // bound method + } else { + // fallback, try calling with 2 args first and catch + argcount = -1; + } + } catch (...) { + argcount = -1; } - return SigmaFunc(R)*sech*sech*hole/(4.0*H); }; + if (argcount == 1) { + // User provided Sigma(R). Wrap it and add vertical profile + + // hole logic here. Keep a copy of pyfunc alive by capturing + // py::object by value. + std::function Sigma = [pyfunc](double R)->double { + py::gil_scoped_acquire gil; + py::object out = pyfunc(R); + return out.cast(); + }; + + SigmaZFunc = [Sigma, H, Rcut, Rwid](double R, double z)->double { + // vertical profile: sech^2 with scale H (same form as original) + double Q = std::exp(-std::fabs(z) / (2.0 * H)); + double sech = 2.0 * Q / (1.0 + Q * Q); + double hole = 1.0; + if (Rcut > 0.0) { + double x = (R - Rcut) / Rwid; + hole = 0.5 * (1.0 + std::tanh(x)); + } + double s = Sigma(R); + return s * sech * sech * hole / (4.0 * H); + }; + + } else if (argcount == 2) { + // User provided SigmaZ(R,z) directly. Use it as-is (no extra + // vertical/hole logic). + py::object pyfunc2 = pyfunc; + SigmaZFunc = [pyfunc2](double R, double z)->double { + py::gil_scoped_acquire gil; + py::object out = pyfunc2(R, z); + return out.cast(); + }; + + } else { + // ambiguous: try calling with 2 args; if that fails try 1 arg + try { + // test call with dummy values + py::gil_scoped_acquire gil; + pyfunc(1.0, 0.0); + // succeeded: treat as SigmaZ(R,z) + py::object pyfunc2 = pyfunc; + SigmaZFunc = [pyfunc2](double R, double z)->double { + py::gil_scoped_acquire gil2; + py::object out = pyfunc2(R, z); + return out.cast(); + }; + } catch (const py::error_already_set &) { + // fallback: try as Sigma(R) + py::object pyfunc1 = pyfunc; + std::function Sigma = [pyfunc1](double R)->double { + py::gil_scoped_acquire gil2; + py::object out = pyfunc1(R); + return out.cast(); + }; + SigmaZFunc = [Sigma, H, Rcut, Rwid](double R, double z)->double { + double Q = std::exp(-std::fabs(z) / (2.0 * H)); + double sech = 2.0 * Q / (1.0 + Q * Q); + double hole = 1.0; + if (Rcut > 0.0) { + double x = (R - Rcut) / Rwid; + hole = 0.5 * (1.0 + std::tanh(x)); + } + double s = Sigma(R); + return s * sech * sech * hole / (4.0 * H); + }; + } + } + } // end python handling + else { + // No Python supplied: use internal choices (plummer/gaussian/toomre) + std::function SigmaFunc; + enum class Type { Plummer, Gaussian, Toomre }; + Type which = Type::Toomre; + + std::map type_map = { + {"plummer", Type::Plummer}, + {"gaussian", Type::Gaussian}, + {"toomre", Type::Toomre} + }; + + std::string type_l = result["type"].as(); + std::transform(type_l.begin(), type_l.end(), type_l.begin(), + [](unsigned char c){ return std::tolower(c); }); + + auto it = type_map.find(type_l); + if (it != type_map.end()) { + which = it->second; + } else { + throw std::runtime_error("Unknown type: " + result["type"].as()); + } + + switch (which) { + case Type::Toomre: + SigmaFunc = [](double R)->double { return 1.0 / std::pow(1.0 + R*R, 1.5); }; + break; + case Type::Gaussian: + SigmaFunc = [](double R)->double { return std::exp(-0.5*R*R); }; + break; + default: + case Type::Plummer: + SigmaFunc = [](double R)->double { return 4.0 / 3.0 / std::pow(1.0 + R*R, 2.0); }; + break; + } + + // Build SigmaZFunc from SigmaFunc, using the same vertical/hole + // logic + SigmaZFunc = [SigmaFunc, H, Rcut, Rwid](double R, double z)->double { + double Q = std::exp(-std::fabs(z) / (2.0 * H)); + double sech = 2.0 * Q / (1.0 + Q * Q); + double hole = 1.0; + if (Rcut > 0.0) { + double x = (R - Rcut) / Rwid; + hole = 0.5 * (1.0 + std::tanh(x)); + } + return SigmaFunc(R) * sech * sech * hole / (4.0 * H); + }; + } // end else internal choice + + // Construct EmpDeproj and evaluate EmpDeproj E(H, Rmin, Rmax, NumR, Nint, SigmaZFunc, type_enum); + // radial evaluation points std::vector r_eval; for (int i = 0; i < Nr; ++i) { double t = (double)i / (Nr - 1); - r_eval.push_back(0.01 + t * 8.0); + r_eval.push_back(Rmin + t * (Rmax - Rmin)); } - + std::ofstream ofs(fname); - for (size_t i = 0; i < r_eval.size(); ++i) + if (!ofs) { + std::cerr << "Failed to open output file: " << fname << std::endl; + return 1; + } + + for (size_t i = 0; i < r_eval.size(); ++i) { ofs << std::setw(16) << r_eval[i] - << std::setw(16) << E.density(r_eval[i]) - << std::setw(16) << E.surfaceDensity(r_eval[i]) - << std::endl; + << std::setw(16) << E.density(r_eval[i]) + << std::setw(16) << E.surfaceDensity(r_eval[i]) + << std::endl; + } ofs.close(); std::cout << "Wrote " << fname << std::endl; - + return 0; } From 1d50ae5c1c9f7107d518cb3e988abd1774799d83 Mon Sep 17 00:00:00 2001 From: mdw Date: Sat, 14 Mar 2026 09:05:53 -0400 Subject: [PATCH 25/78] Fix typo in parameter assignment --- expui/BiorthBasis.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 24608fb21..ca7b01472 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1520,7 +1520,7 @@ namespace BasisClasses if (conf["dtype" ]) dtype = conf["dtype" ].as(); if (conf["vflag" ]) vflag = conf["vflag" ].as(); if (conf["pyname" ]) pyname = conf["pyname" ].as(); - if (conf["pyproj" ]) pyname = conf["pyproj" ].as(); + if (conf["pyproj" ]) pyproj = conf["pyproj" ].as(); if (conf["pcavar"] ) pcavar = conf["pcavar" ].as(); if (conf["subsamp"] ) sampT = conf["subsamp" ].as(); From eb1f9eb2d612e4adb302fd27bd9f59a20a459b88 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 14 Mar 2026 10:54:02 -0400 Subject: [PATCH 26/78] Remove the initial test function; add doc strings to the help stanza on the standalone tests --- utils/Test/CMakeLists.txt | 4 +- utils/Test/testDeprojPlummer.cc | 77 ---------------- utils/Test/testDeproject.cc | 153 ++++++++++++++++++++++---------- utils/Test/testEmp.cc | 7 +- utils/Test/testEmpDeproj.cc | 4 +- 5 files changed, 113 insertions(+), 132 deletions(-) delete mode 100644 utils/Test/testDeprojPlummer.cc diff --git a/utils/Test/CMakeLists.txt b/utils/Test/CMakeLists.txt index 45531c968..123e6120d 100644 --- a/utils/Test/CMakeLists.txt +++ b/utils/Test/CMakeLists.txt @@ -1,6 +1,6 @@ set(bin_PROGRAMS testBarrier expyaml orthoTest testDeproj - testDeprojPlum testEmpDeproj testEmp) + testEmpDeproj testEmp) set(common_LINKLIB OpenMP::OpenMP_CXX MPI::MPI_CXX expui exputil yaml-cpp ${VTK_LIBRARIES}) @@ -34,8 +34,6 @@ add_executable(testBarrier test_barrier.cc) add_executable(expyaml test_config.cc) add_executable(orthoTest orthoTest.cc Biorth2Ortho.cc) add_executable(testDeproj testDeproject.cc CubicSpline.cc Deprojector.cc) -add_executable(testDeprojPlum testDeprojPlummer.cc CubicSpline.cc - Deprojector.cc) add_executable(testEmpDeproj testEmpDeproj.cc CubicSpline.cc Deprojector.cc EmpDeproj.cc) add_executable(testEmp testEmp.cc EmpDeproj.cc) diff --git a/utils/Test/testDeprojPlummer.cc b/utils/Test/testDeprojPlummer.cc deleted file mode 100644 index 0e34cff9e..000000000 --- a/utils/Test/testDeprojPlummer.cc +++ /dev/null @@ -1,77 +0,0 @@ -#include -#include -#include -#include -#include - -#include "Deprojector.H" - -using namespace Deproject; - -int main() -{ - std::function SigmaFunc, dSigmaFunc, RhoFunc; - - enum class Type { Plummer, Gaussian, Toomre }; - Type which = Type::Plummer; - - switch (which) { - case Type::Toomre: - // Test function - SigmaFunc = [](double R)->double - { return 1.0 / std::pow(1.0 + R*R, 1.5); }; - // Analytic derivative - dSigmaFunc = [](double R)->double - { return -3.0 * R / std::pow(1.0 + R*R, 2.5); }; - // Expected result - RhoFunc = [](double r)->double - { return 2.0 / std::pow(1.0 + r*r, 2.0) / M_PI; }; - break; - case Type::Gaussian: - // Test function - SigmaFunc = [](double R)->double - { return exp(-0.5*R*R); }; - // Analytic derivative - dSigmaFunc = [](double R)->double - { return -R*exp(-0.5*R*R); }; - // Expected result - RhoFunc = [](double r)->double - { return exp(-0.5*r*r)/sqrt(2.0*M_PI); }; - break; - default: - case Type::Plummer: - // Test function - SigmaFunc = [](double R)->double - { return 4.0 / 3.0 / std::pow(1.0 + R*R, 2.0); }; - // Analytic derivative - dSigmaFunc = [](double R)->double - { return -16.0 * R / 3.0 / std::pow(1.0 + R*R, 3.0); }; - // Expected result - RhoFunc = [](double r)->double - { return 1.0 / std::pow(1.0 + r*r, 2.5); }; - break; - } - - Deprojector D(SigmaFunc, dSigmaFunc, /*R_data_min=*/0.01, /*R_data_max=*/10.0, - /*R_max_extend=*/50.0, /*tail_power=*/-4.0, /*Ngrid=*/6000); - - std::vector r_eval; - int Nr = 150; - for (int i = 0; i < Nr; ++i) { - double t = (double)i / (Nr - 1); - r_eval.push_back(0.01 + t * 8.0); - } - auto rho = D.rho(r_eval); - - std::ofstream ofs("rho_test.txt"); - for (size_t i = 0; i < r_eval.size(); ++i) { - ofs << std::setw(16) << r_eval[i] - << std::setw(16) << rho[i] - << std::setw(16) << RhoFunc(r_eval[i]) - << std::endl; - } - ofs.close(); - std::cout << "Wrote rho_test.txt\n"; - - return 0; -} diff --git a/utils/Test/testDeproject.cc b/utils/Test/testDeproject.cc index 309c9cf88..046855022 100644 --- a/utils/Test/testDeproject.cc +++ b/utils/Test/testDeproject.cc @@ -1,63 +1,122 @@ #include +#include #include #include #include #include "Deprojector.H" +#include "cxxopts.H" using namespace Deproject; -int main() +int main(int argc, char* argv[]) { - // Example A: construct from sampled data - { - std::vector Rdata, Sigma; - int Ndata = 2000; - double Rmin = 0.01, Rmax = 10.0; - for (int i = 0; i < Ndata; ++i) { - double t = (double)i / (Ndata - 1); - double r = Rmin + t * (Rmax - Rmin); - Rdata.push_back(r); - Sigma.push_back(std::pow(1.0 + r*r, -1.5)); - } - - Deprojector D(Rdata, Sigma, /*R_max_extend=*/50.0, /*tail_power=*/-4.0, /*Ngrid=*/6000); - - std::vector r_eval; - int Nr = 150; - for (int i = 0; i < Nr; ++i) { - double t = (double)i / (Nr - 1); - r_eval.push_back(0.01 + t * 8.0); - } - auto rho = D.rho(r_eval); - - std::ofstream ofs("rho_from_sampled.txt"); - for (size_t i = 0; i < r_eval.size(); ++i) ofs << r_eval[i] << " " << rho[i] << "\n"; - ofs.close(); - std::cout << "Wrote rho_from_sampled.txt\n"; + // Command-line options + std::string type; + cxxopts::Options options("testDeprojector", + "Test the Deproject class for various surface density profiles.\n" + "Independent implemenation for comparison with the EmpCylSL-based\n" + "EmpDeproj class.\n"); + + options.add_options() + ("h,help", "Print help") + ("type", "Surface density type (plummer, gaussian, toomre)", cxxopts::value(type)->default_value("plummer")); + + auto result = options.parse(argc, argv); + + if (result.count("help")) { + std::cout << options.help() << std::endl; + return 0; + } + + // Density function and derivative functors + // + std::function SigmaFunc, dSigmaFunc, RhoFunc; + + // Convert type id string to lower case + // + std::transform(type.begin(), type.end(), type.begin(), + [](unsigned char c){ return std::tolower(c); }); + + // Test densities + // + enum class Type { Plummer, Gaussian, Toomre }; + + // Reflect type string to enum + // + std::map type_map = { + {"plummer", Type::Plummer}, + {"gaussian", Type::Gaussian}, + {"toomre", Type::Toomre} + }; + + Type which = Type::Plummer; + + auto it = type_map.find(type); + + if (it != type_map.end()) { + which = it->second; + } else { + throw std::runtime_error("Unknown type: " + type); + } + + switch (which) { + case Type::Toomre: + // Test function + SigmaFunc = [](double R)->double + { return 1.0 / std::pow(1.0 + R*R, 1.5); }; + // Analytic derivative + dSigmaFunc = [](double R)->double + { return -3.0 * R / std::pow(1.0 + R*R, 2.5); }; + // Expected result + RhoFunc = [](double r)->double + { return 2.0 / std::pow(1.0 + r*r, 2.0) / M_PI; }; + break; + case Type::Gaussian: + // Test function + SigmaFunc = [](double R)->double + { return exp(-0.5*R*R); }; + // Analytic derivative + dSigmaFunc = [](double R)->double + { return -R*exp(-0.5*R*R); }; + // Expected result + RhoFunc = [](double r)->double + { return exp(-0.5*r*r)/sqrt(2.0*M_PI); }; + break; + default: + case Type::Plummer: + // Test function + SigmaFunc = [](double R)->double + { return 4.0 / 3.0 / std::pow(1.0 + R*R, 2.0); }; + // Analytic derivative + dSigmaFunc = [](double R)->double + { return -16.0 * R / 3.0 / std::pow(1.0 + R*R, 3.0); }; + // Expected result + RhoFunc = [](double r)->double + { return 1.0 / std::pow(1.0 + r*r, 2.5); }; + break; + } + + Deprojector D(SigmaFunc, dSigmaFunc, /*R_data_min=*/0.01, /*R_data_max=*/10.0, + /*R_max_extend=*/50.0, /*tail_power=*/-4.0, /*Ngrid=*/6000); + + std::vector r_eval; + int Nr = 150; + for (int i = 0; i < Nr; ++i) { + double t = (double)i / (Nr - 1); + r_eval.push_back(0.01 + t * 8.0); } + auto rho = D.rho(r_eval); - // Example B: construct from analytic functor + analytic derivative - { - auto SigmaFunc = [](double R)->double { return std::pow(1.0 + R*R, -1.5); }; - auto dSigmaFunc = [](double R)->double { return -3.0 * R / std::pow(1.0 + R*R, -2.5); }; // analytic derivative - - Deprojector D(SigmaFunc, dSigmaFunc, /*R_data_min=*/0.01, /*R_data_max=*/10.0, - /*R_max_extend=*/50.0, /*tail_power=*/-4.0, /*Ngrid=*/6000); - - std::vector r_eval; - int Nr = 150; - for (int i = 0; i < Nr; ++i) { - double t = (double)i / (Nr - 1); - r_eval.push_back(0.01 + t * 8.0); - } - auto rho = D.rho(r_eval); - - std::ofstream ofs("rho_from_functor.txt"); - for (size_t i = 0; i < r_eval.size(); ++i) ofs << r_eval[i] << " " << rho[i] << "\n"; - ofs.close(); - std::cout << "Wrote rho_from_functor.txt\n"; + std::ofstream ofs("rho_test.txt"); + for (size_t i = 0; i < r_eval.size(); ++i) { + ofs << std::setw(16) << r_eval[i] + << std::setw(16) << rho[i] + << std::setw(16) << RhoFunc(r_eval[i]) + << std::endl; } + ofs.close(); + std::cout << "Wrote rho_test.txt\n"; return 0; } diff --git a/utils/Test/testEmp.cc b/utils/Test/testEmp.cc index 9ee505732..c3e8e9309 100644 --- a/utils/Test/testEmp.cc +++ b/utils/Test/testEmp.cc @@ -28,9 +28,10 @@ int main(int argc, char* argv[]) // Parse command-line options cxxopts::Options options("testDonut", - "Test the EmpDeproj class for an inner donut-shaped " - "density distribution, using the Toomre profile as " - "a test case."); + "Test the EmpDeproj class for any Python-supplied density function.\n" + "If the Python module is not supplied, one of the hard-coded\n" + "functions: Plummer, Gaussian, Toomre may be selected. The internal\n" + "Abel method: Derivative, Substraction, or IBP may be chosen as well.\n"); options.add_options() ("h,help", "Print help") diff --git a/utils/Test/testEmpDeproj.cc b/utils/Test/testEmpDeproj.cc index 966313028..54a3a8cc2 100644 --- a/utils/Test/testEmpDeproj.cc +++ b/utils/Test/testEmpDeproj.cc @@ -27,8 +27,8 @@ int main(int argc, char* argv[]) // Parse command-line options // cxxopts::Options options("testEmpDeproj", - "Test the EmpDeproj class against Deprojector" - "for various surface density profiles."); + "Test the EmpDeproj class against Deproject" + "for various surface density profiles.\n"); options.add_options() ("h,help", "Print help") ("type", "Surface density type (plummer, gaussian, toomre)", cxxopts::value(type)->default_value("toomre")) From 6402d49b314168df35a6f44e646cf847c4535c83 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Sat, 14 Mar 2026 13:14:30 -0400 Subject: [PATCH 27/78] Only print deprojection info from one node --- expui/BiorthBasis.cc | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index ca7b01472..f8a2e426d 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1745,9 +1745,10 @@ namespace BasisClasses else if (PTYPE == DeprojType::toomre) { model = std::make_shared(1.0, H, 5.0); } else if (PTYPE == DeprojType::python) { - model = std::make_shared(pyproj, acyl); - std::cout << "Using AxiSymPyModel for deprojection from Python function <" - << pyname << ">" << std::endl; + model = std::make_shared(pyproj, 1.0); + if (myid==0) + std::cout << "---- Using AxiSymPyModel for deprojection from " + << "Python module <" << pyproj << ">" << std::endl; } else { // Default to exponential model = std::make_shared(1.0, H); } From c934a7961c5a9d26fcdb94d36567262c2b5d418c Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 14 Mar 2026 18:12:59 -0400 Subject: [PATCH 28/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index dd12aaf70..b294d4da9 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1253,7 +1253,7 @@ namespace BasisClasses "coefCompute", "coefMaster", "pyname", - "pyproj" + "pyproj", "nint", "totalCovar", "fullCovar" From 5a95b504e5b6288eae2e952a7768c65305f6acf0 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 14 Mar 2026 18:14:38 -0400 Subject: [PATCH 29/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/Test/testDeproject.cc | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/utils/Test/testDeproject.cc b/utils/Test/testDeproject.cc index 046855022..c9ae00e40 100644 --- a/utils/Test/testDeproject.cc +++ b/utils/Test/testDeproject.cc @@ -7,6 +7,10 @@ #include "Deprojector.H" #include "cxxopts.H" +namespace { + constexpr double pi = std::acos(-1.0); +} + using namespace Deproject; int main(int argc, char* argv[]) @@ -70,7 +74,7 @@ int main(int argc, char* argv[]) { return -3.0 * R / std::pow(1.0 + R*R, 2.5); }; // Expected result RhoFunc = [](double r)->double - { return 2.0 / std::pow(1.0 + r*r, 2.0) / M_PI; }; + { return 2.0 / std::pow(1.0 + r*r, 2.0) / pi; }; break; case Type::Gaussian: // Test function @@ -81,7 +85,7 @@ int main(int argc, char* argv[]) { return -R*exp(-0.5*R*R); }; // Expected result RhoFunc = [](double r)->double - { return exp(-0.5*r*r)/sqrt(2.0*M_PI); }; + { return exp(-0.5*r*r)/sqrt(2.0*pi); }; break; default: case Type::Plummer: From 402dc2511b6f17a90c24a00c56f2d73a5553a930 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 14 Mar 2026 18:15:04 -0400 Subject: [PATCH 30/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/Test/testEmpDeproj.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/Test/testEmpDeproj.cc b/utils/Test/testEmpDeproj.cc index 54a3a8cc2..c8201c01e 100644 --- a/utils/Test/testEmpDeproj.cc +++ b/utils/Test/testEmpDeproj.cc @@ -27,7 +27,7 @@ int main(int argc, char* argv[]) // Parse command-line options // cxxopts::Options options("testEmpDeproj", - "Test the EmpDeproj class against Deproject" + "Test the EmpDeproj class against Deproject " "for various surface density profiles.\n"); options.add_options() ("h,help", "Print help") From 4513d2a5dffa4776a022e63cf99b2cb28d485248 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 14 Mar 2026 18:15:30 -0400 Subject: [PATCH 31/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/Test/testEmp.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/Test/testEmp.cc b/utils/Test/testEmp.cc index c3e8e9309..81c293fa1 100644 --- a/utils/Test/testEmp.cc +++ b/utils/Test/testEmp.cc @@ -31,7 +31,7 @@ int main(int argc, char* argv[]) "Test the EmpDeproj class for any Python-supplied density function.\n" "If the Python module is not supplied, one of the hard-coded\n" "functions: Plummer, Gaussian, Toomre may be selected. The internal\n" - "Abel method: Derivative, Substraction, or IBP may be chosen as well.\n"); + "Abel method: Derivative, Subtraction, or IBP may be chosen as well.\n"); options.add_options() ("h,help", "Print help") From e719f9f2503b598e9d53205586909fe5b9d86dea Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 14 Mar 2026 18:16:14 -0400 Subject: [PATCH 32/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/Test/CubicSpline.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/Test/CubicSpline.H b/utils/Test/CubicSpline.H index f9fde2f9a..38bd77b11 100644 --- a/utils/Test/CubicSpline.H +++ b/utils/Test/CubicSpline.H @@ -14,7 +14,7 @@ namespace Deproject // set data (x must be strictly increasing) void set_data(const std::vector& x_in, const std::vector& y_in); - // evaluate spline and its derivative (xx should lie within [xmin(), xmax()], but endpoints are clamped) + // evaluate spline and its derivative (xx is expected in [xmin(), xmax()]; values outside are extrapolated using the end intervals) double eval(double xx) const; double deriv(double xx) const; From 627de499b1c611b8f9d6e71f257bcd105f172c65 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 14 Mar 2026 18:16:42 -0400 Subject: [PATCH 33/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/Test/testEmpDeproj.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/Test/testEmpDeproj.cc b/utils/Test/testEmpDeproj.cc index c8201c01e..9e9ef0374 100644 --- a/utils/Test/testEmpDeproj.cc +++ b/utils/Test/testEmpDeproj.cc @@ -134,7 +134,7 @@ int main(int argc, char* argv[]) } Deprojector D(SigmaFunc, dSigmaFunc, /*R_data_min=*/0.01, /*R_data_max=*/10.0, - /*R_max_extend=*/50.0, /*tail_power=*/-4.0, /*Ngrid=*/6000); + /*R_max_extend=*/50.0, /*tail_power=*/-4.0, /*Ngrid=*/Ngrid); auto SigmaZFunc = [SigmaFunc, H, Rcut, Rwid](double R, double z)->double { double Q = exp(-std::fabs(z)/(2.0*H)); From f5a39b41f0bfc4f9b2d5256e5659fadee4ef3608 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 14 Mar 2026 18:17:28 -0400 Subject: [PATCH 34/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/Test/testEmpDeproj.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/utils/Test/testEmpDeproj.cc b/utils/Test/testEmpDeproj.cc index 9e9ef0374..e0217cb7e 100644 --- a/utils/Test/testEmpDeproj.cc +++ b/utils/Test/testEmpDeproj.cc @@ -150,8 +150,8 @@ int main(int argc, char* argv[]) std::vector r_eval; for (int i = 0; i < Nr; ++i) { - double t = (double)i / (Nr - 1); - r_eval.push_back(0.01 + t * 8.0); + double t = (Nr > 1) ? static_cast(i) / (Nr - 1) : 0.0; + r_eval.push_back(Rmin + t * (Rmax - Rmin)); } auto rho = D.rho(r_eval); From b6f5c037b794ba0f0726a9b2efa0a79d7930e934 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 14 Mar 2026 18:17:46 -0400 Subject: [PATCH 35/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/Test/testEmpDeproj.cc | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/utils/Test/testEmpDeproj.cc b/utils/Test/testEmpDeproj.cc index e0217cb7e..1ae2b0158 100644 --- a/utils/Test/testEmpDeproj.cc +++ b/utils/Test/testEmpDeproj.cc @@ -24,6 +24,9 @@ int main(int argc, char* argv[]) double H, Rmin, Rmax, Rcut, Rwid; int Nr, Ngrid, NumR, Nint; + // Define pi in a portable way instead of relying on non-standard M_PI + constexpr double pi = std::acos(-1.0); + // Parse command-line options // cxxopts::Options options("testEmpDeproj", @@ -105,8 +108,8 @@ int main(int argc, char* argv[]) dSigmaFunc = [](double R)->double { return -3.0 * R / std::pow(1.0 + R*R, 2.5); }; // Expected result - RhoFunc = [](double r)->double - { return 2.0 / std::pow(1.0 + r*r, 2.0) / M_PI; }; + RhoFunc = [pi](double r)->double + { return 2.0 / std::pow(1.0 + r*r, 2.0) / pi; }; break; case Type::Gaussian: // Test function @@ -116,8 +119,8 @@ int main(int argc, char* argv[]) dSigmaFunc = [](double R)->double { return -R*exp(-0.5*R*R); }; // Expected result - RhoFunc = [](double r)->double - { return exp(-0.5*r*r)/sqrt(2.0*M_PI); }; + RhoFunc = [pi](double r)->double + { return exp(-0.5*r*r)/sqrt(2.0*pi); }; break; default: case Type::Plummer: From 3dcafad23fa220b6c4309514ba426fbc20cc04b6 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sat, 14 Mar 2026 18:18:03 -0400 Subject: [PATCH 36/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/Test/EmpDeproj.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/Test/EmpDeproj.cc b/utils/Test/EmpDeproj.cc index b43962a5e..51244508b 100644 --- a/utils/Test/EmpDeproj.cc +++ b/utils/Test/EmpDeproj.cc @@ -53,7 +53,7 @@ EmpDeproj::EmpDeproj(double H, double RMIN, double RMAX, int NUMR, int NINT, Linear1d surf(rl, sigI); surfRg = surf; - // Now, compute Abel inverion integral + // Now, compute Abel inversion integral // for (int i=0; i Date: Sat, 14 Mar 2026 18:18:27 -0400 Subject: [PATCH 37/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/Test/EmpDeproj.cc | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/utils/Test/EmpDeproj.cc b/utils/Test/EmpDeproj.cc index 51244508b..7c6d2d5a2 100644 --- a/utils/Test/EmpDeproj.cc +++ b/utils/Test/EmpDeproj.cc @@ -9,6 +9,10 @@ #include #include +#ifndef M_PI +#define M_PI std::acos(-1.0) +#endif + #include "EmpDeproj.H" #include "gaussQ.H" From 79e0233c0b1635b9bf00f273ecce26dd900b2056 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 14 Mar 2026 22:46:59 +0000 Subject: [PATCH 38/78] Initial plan From 9c0d024e6f38aba68da09c8c9a551c7fead715b3 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 14 Mar 2026 22:49:17 +0000 Subject: [PATCH 39/78] Fix clang build: change constexpr to const for pi = std::acos(-1.0) Co-authored-by: The9Cat <25960766+The9Cat@users.noreply.github.com> --- utils/Test/testDeproject.cc | 2 +- utils/Test/testEmpDeproj.cc | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/utils/Test/testDeproject.cc b/utils/Test/testDeproject.cc index c9ae00e40..8675b9643 100644 --- a/utils/Test/testDeproject.cc +++ b/utils/Test/testDeproject.cc @@ -8,7 +8,7 @@ #include "cxxopts.H" namespace { - constexpr double pi = std::acos(-1.0); + const double pi = std::acos(-1.0); } using namespace Deproject; diff --git a/utils/Test/testEmpDeproj.cc b/utils/Test/testEmpDeproj.cc index 1ae2b0158..a502331ba 100644 --- a/utils/Test/testEmpDeproj.cc +++ b/utils/Test/testEmpDeproj.cc @@ -25,7 +25,7 @@ int main(int argc, char* argv[]) int Nr, Ngrid, NumR, Nint; // Define pi in a portable way instead of relying on non-standard M_PI - constexpr double pi = std::acos(-1.0); + const double pi = std::acos(-1.0); // Parse command-line options // From 45a578f1ac29d591463bd7043af6c2b167aa1df8 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 15 Mar 2026 12:26:17 -0400 Subject: [PATCH 40/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index b294d4da9..f3fa14c46 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1224,7 +1224,7 @@ namespace BasisClasses "ignore", "deproject", "dmodel", - "logr", + "ppow", "pcavar", "pcaeof", "pcavtk", From 36120a312e42ef9b03cb53d4e14129420504a5b0 Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 15 Mar 2026 12:26:47 -0400 Subject: [PATCH 41/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index f3fa14c46..09e285e1e 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1500,7 +1500,7 @@ namespace BasisClasses if (conf["ashift" ]) ashift = conf["ashift" ].as(); if (conf["rfactor" ]) rfactor = conf["rfactor" ].as(); if (conf["rtrunc" ]) rtrunc = conf["rtrunc" ].as(); - if (conf["pow" ]) ppow = conf["ppow" ].as(); + if (conf["ppow" ]) ppow = conf["ppow" ].as(); if (conf["mtype" ]) mtype = conf["mtype" ].as(); if (conf["dtype" ]) dtype = conf["dtype" ].as(); if (conf["vflag" ]) vflag = conf["vflag" ].as(); From 1b1f297a91bc5ad37ce455e0dd7fd2bf66e68eeb Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 15 Mar 2026 12:27:23 -0400 Subject: [PATCH 42/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- expui/BiorthBasis.cc | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 09e285e1e..93245b321 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1731,6 +1731,16 @@ namespace BasisClasses else if (PTYPE == DeprojType::toomre) { model = std::make_shared(1.0, H, 5.0); } else if (PTYPE == DeprojType::python) { + if (pyproj.empty()) { + if (myid==0) { + std::cout << "DeprojType is set to 'python' but no Python " + << "projection module name (pyname/pyproj) was provided." + << std::endl; + } + throw std::runtime_error( + "Cylindrical::initialize: DeprojType 'python' requires a " + "non-empty Python module name (pyname/pyproj)."); + } model = std::make_shared(pyproj, 1.0); if (myid==0) std::cout << "---- Using AxiSymPyModel for deprojection from " From fdb16edd212d78d59092c0cc546db8f9532cc43a Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 15 Mar 2026 12:27:55 -0400 Subject: [PATCH 43/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/Test/EmpDeproj.cc | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/utils/Test/EmpDeproj.cc b/utils/Test/EmpDeproj.cc index 7c6d2d5a2..46f28835a 100644 --- a/utils/Test/EmpDeproj.cc +++ b/utils/Test/EmpDeproj.cc @@ -8,6 +8,7 @@ #include #include #include +#include #ifndef M_PI #define M_PI std::acos(-1.0) @@ -22,6 +23,23 @@ EmpDeproj::EmpDeproj(double H, double RMIN, double RMAX, int NUMR, int NINT, std::function func, AbelType type) { + // Validate input arguments to avoid domain errors and division by zero + if (H <= 0.0) { + throw std::invalid_argument("EmpDeproj: H must be > 0"); + } + if (RMIN <= 0.0) { + throw std::invalid_argument("EmpDeproj: RMIN must be > 0"); + } + if (RMAX <= RMIN) { + throw std::invalid_argument("EmpDeproj: RMAX must be > RMIN > 0"); + } + if (NUMR < 2) { + throw std::invalid_argument("EmpDeproj: NUMR must be >= 2"); + } + if (NINT < 1) { + throw std::invalid_argument("EmpDeproj: NINT must be >= 1"); + } + LegeQuad lq(NINT); std::vector rr(NUMR), rl(NUMR), sigI(NUMR), rhoI(NUMR, 0.0); From d2406792371199d2cb7629a1230049ee9f4a0a3b Mon Sep 17 00:00:00 2001 From: Martin Weinberg Date: Sun, 15 Mar 2026 12:28:23 -0400 Subject: [PATCH 44/78] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- utils/Test/testEmp.cc | 3 +++ 1 file changed, 3 insertions(+) diff --git a/utils/Test/testEmp.cc b/utils/Test/testEmp.cc index 81c293fa1..524d72e52 100644 --- a/utils/Test/testEmp.cc +++ b/utils/Test/testEmp.cc @@ -275,6 +275,9 @@ int main(int argc, char* argv[]) EmpDeproj E(H, Rmin, Rmax, NumR, Nint, SigmaZFunc, type_enum); // radial evaluation points + if (Nr < 2) { + throw std::runtime_error("Nr must be at least 2 (received " + std::to_string(Nr) + ")"); + } std::vector r_eval; for (int i = 0; i < Nr; ++i) { double t = (double)i / (Nr - 1); From 18b14555ab8ced6b606bc998ed00772d977e77db Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 16 Mar 2026 16:41:46 -0400 Subject: [PATCH 45/78] Accidentally deleted the 'logr' key; restored --- expui/BiorthBasis.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 93245b321..528f354e5 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1223,6 +1223,7 @@ namespace BasisClasses "expcond", "ignore", "deproject", + "logr", "dmodel", "ppow", "pcavar", From ab919512abc3c0bd962fd9e9fd313ce4ad96cf12 Mon Sep 17 00:00:00 2001 From: Michael Petersen Date: Mon, 23 Mar 2026 14:06:28 +0000 Subject: [PATCH 46/78] Apply suggestions from code review This applies a couple more copilot suggestions. Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- exputil/EmpCylSL.cc | 2 +- utils/Test/EmpDeproj.H | 17 +++++++++++------ utils/Test/testDeproject.cc | 1 + 3 files changed, 13 insertions(+), 7 deletions(-) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 08f1c4f20..8fc5e5ba2 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -485,7 +485,7 @@ void EmpCylSL::create_deprojection(double H, double Rf, int NUMR, int NINT, Linear1d surf(rl, sigI); - // Now, compute Abel inverion integral + // Now, compute Abel inversion integral // for (int i=0; i +#include +#include #include "interp.H" class EmpDeproj @@ -22,19 +24,22 @@ public: ~EmpDeproj() {} //! Evaluate deprojected density at a single radius - double density(double R) + double density(double R) const { - return densRg.eval(log(R)); + assert(R > 0.0); + return densRg.eval(std::log(R)); } - double surfaceDensity(double R) + double surfaceDensity(double R) const { - return surfRg.eval(log(R)); + assert(R > 0.0); + return surfRg.eval(std::log(R)); } //! Evaluate deprojected mass at a single radius - double mass(double R) + double mass(double R) const { - return massRg.eval(log(R)); + assert(R > 0.0); + return massRg.eval(std::log(R)); } }; diff --git a/utils/Test/testDeproject.cc b/utils/Test/testDeproject.cc index 8675b9643..ed606a312 100644 --- a/utils/Test/testDeproject.cc +++ b/utils/Test/testDeproject.cc @@ -3,6 +3,7 @@ #include #include #include +#include #include "Deprojector.H" #include "cxxopts.H" From 446d52bc79139a46cca0079b8a90f796091ae9b1 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Mon, 23 Mar 2026 14:12:14 +0000 Subject: [PATCH 47/78] add version number increment proposal --- include/EmpCylSL.H | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/include/EmpCylSL.H b/include/EmpCylSL.H index a05179d1f..e6822ac01 100644 --- a/include/EmpCylSL.H +++ b/include/EmpCylSL.H @@ -239,7 +239,9 @@ protected: bool ReadH5Cache(); //! Cache versioning - inline static const std::string Version = "1.0"; + // Version 1.0 corresponds to bases generated with the first public release of EXP (7.8). + // Version 1.1 corresponds to bases generated with improved Abel inversion and deprojection methods (7.10). + inline static const std::string Version = "1.1"; //! The cache file name std::string cachefile; From cb3c977ff3d2bdfa41143acf1fb4d06fc0dfb954 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Mon, 23 Mar 2026 14:24:59 +0000 Subject: [PATCH 48/78] erroneous const breaking compile --- utils/Test/EmpDeproj.H | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/utils/Test/EmpDeproj.H b/utils/Test/EmpDeproj.H index 6c6c2c818..0120021fb 100644 --- a/utils/Test/EmpDeproj.H +++ b/utils/Test/EmpDeproj.H @@ -24,20 +24,20 @@ public: ~EmpDeproj() {} //! Evaluate deprojected density at a single radius - double density(double R) const + double density(double R) { assert(R > 0.0); return densRg.eval(std::log(R)); } - double surfaceDensity(double R) const + double surfaceDensity(double R) { assert(R > 0.0); return surfRg.eval(std::log(R)); } //! Evaluate deprojected mass at a single radius - double mass(double R) const + double mass(double R) { assert(R > 0.0); return massRg.eval(std::log(R)); From ccd4630ac91ee88c0588ab8151b685cdd90feb3d Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 23 Mar 2026 18:34:04 -0400 Subject: [PATCH 49/78] Added first implemenation for saving and reading DTYPE and Python md5sum info --- expui/BiorthBasis.cc | 135 ++++++++++++++++++++++++++++++++++++++++- exputil/CMakeLists.txt | 2 +- exputil/getmd5sum.cc | 35 +++++++++++ include/exputils.H | 4 ++ 4 files changed, 174 insertions(+), 2 deletions(-) create mode 100644 exputil/getmd5sum.cc diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 528f354e5..70c251481 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1596,7 +1596,105 @@ namespace BasisClasses // Attempt to read EOF cache // - if (sl->read_cache() == 0) { + int cache_status = sl->read_cache(); + + // Define the HighFive mapping for DTYPE + // + HighFive::EnumType disk_type({ + {"constant", DiskType::constant}, + {"gaussian", DiskType::gaussian}, + {"mn", DiskType::mn}, + {"exponential", DiskType::exponential}, + {"doubleexpon", DiskType::doubleexpon}, + {"diskbulge", DiskType::diskbulge}, + {"python", DiskType::python} + }); + + std::map disk_type_string = { + {DiskType::constant, "constant"}, + {DiskType::gaussian, "gaussian"}, + {DiskType::mn, "mn"}, + {DiskType::exponential, "exponential"}, + {DiskType::doubleexpon, "doubleexpon"}, + {DiskType::diskbulge, "diskbulge"}, + {DiskType::python, "python"} + }; + + + // Checking for cache consistency with current DiskType and Python module (if applicable) + // + if (cache_status == 1) { + + try { + // Open the cache file for reading the Python metadata + // + HighFive::File file(cachename, HighFive::File::ReadOnly); + + // Sanity: check that the DiskType attribute exists + // + if (!file.hasAttribute("DiskType")) { + if (myid==0) { + std::cout << "---- Cylindrical: DiskType attribute not found in cache file <" << cachename << ">. " + << "---- This may indicate an old cache file created before DiskType metadata was added. " + << "---- We will continue...but consider remaking the cache to avoid confusion." << std::endl; + } + } else { + // Open existing DiskType attribute + // + auto read_attr = file.getAttribute("DiskType"); + + DiskType loaded_dtype; + read_attr.read(loaded_dtype); + + if (loaded_dtype != DTYPE) { + if (myid==0) { + std::cout << "---- Cylindrical: DiskType for cache file <" << cachename << "> is <" + << disk_type_string[loaded_dtype] << ">, which does not match the requested DiskType <" + << disk_type_string[DTYPE] << ">. Forcing cache recomputation." << std::endl; + } + // Force cache recomputation + cache_status = 0; + } + else if (loaded_dtype == DiskType::python) { + // Get the pyname attribute + auto read_attr = file.getAttribute("pyname"); + std::string loaded_pyname; + read_attr.read(loaded_pyname); + + std::string loaded_md5, current_md5; + + // Get the md5sum for requested Python module + try { + current_md5 = get_md5sum(pyname); + loaded_md5 = get_md5sum(loaded_pyname); + } catch (const std::runtime_error& e) { + std::cerr << "Error: " << e.what() << std::endl; + } + + // Check that the md5sums match for the current Python + // module and the loaded Python module used to create the + // cache. If they do not match, force cache recomputation + // to ensure consistency with the current Python module. + if (current_md5 != loaded_md5) { + if (myid==0) { + std::cout << "---- Cylindrical: Python module for disk density has changed since cache creation." << std::endl + << "---- Current module: <" << pyname << ">, md5sum: " << current_md5 << std::endl + << "---- Loaded module: <" << loaded_pyname << ">, md5sum: " << loaded_md5 << std::endl + << "---- Forcing cache recomputation to ensure consistency with current Python module." << std::endl; + } + cache_status = 0; + } + } + } + } + catch (const HighFive::Exception& err) { + std::cerr << "---- BiorthBasis inconsistency in Cylindrical cache: " << err.what() << std::endl; + std::cerr << "---- Forcing cache recomputation." << std::endl; + cache_status = 0; // Fallback... + } + } + + if (cache_status == 0) { // Remake cylindrical basis // @@ -1771,6 +1869,41 @@ namespace BasisClasses }; sl->generate_eof(rnum, pnum, tnum, f); + + try { + // Open the cache file for writing the Python metadata + // + HighFive::File file(cachename, HighFive::File::ReadWrite); + + // Create a Scalar DataSpace for the single value + HighFive::DataSpace space = HighFive::DataSpace::From(DTYPE); + + // Create the attribute using the explicit (Name, DataSpace, + // DataType) signature using the base class 'DataType' to + // match the expected signature + auto attr = file.createAttribute("DiskType", space, (HighFive::DataType)disk_type); + + // Write the actual enum value + attr.write(DTYPE); + + // Get the md5sum for the Python module + if (DTYPE == DiskType::python) { + try { + std::string md5_hash = get_md5sum(pyname); + file.createAttribute + ("pyname", + HighFive::DataSpace::From(pyname)).write(pyname); + } catch (const std::runtime_error& e) { + std::cerr << "Error: " << e.what() << std::endl; + std::cerr << "Can not write the md5 hash to HDF5" << std::endl; + } + } + + } catch (const HighFive::Exception& err) { + std::cerr << err.what() << std::endl; + std::cerr << "Error writing metadata to cache file <" << cachename + << std::endl; + } } // Orthogonality sanity check diff --git a/exputil/CMakeLists.txt b/exputil/CMakeLists.txt index 0a029c2ec..3ba43a308 100644 --- a/exputil/CMakeLists.txt +++ b/exputil/CMakeLists.txt @@ -14,7 +14,7 @@ set(UTIL_SRC nrutil.cc elemfunc.cc euler.cc euler_slater.cc TransformFFT.cc QDHT.cc YamlCheck.cc # Hankel.cc parseVersionString.cc EXPmath.cc laguerre_polynomial.cpp YamlConfig.cc orthoTest.cc OrthoFunction.cc VtkGrid.cc - Sutils.cc fpetrap.cc) + Sutils.cc fpetrap.cc getmd5sum.cc) if(HAVE_VTK) list(APPEND UTIL_SRC VtkPCA.cc) diff --git a/exputil/getmd5sum.cc b/exputil/getmd5sum.cc new file mode 100644 index 000000000..b8ab2ca4a --- /dev/null +++ b/exputil/getmd5sum.cc @@ -0,0 +1,35 @@ +#include +#include +#include +#include +#include +#include + +std::string get_md5sum(const std::string& filename) +{ + // Command to execute: md5sum + std::string command = "md5sum " + filename; + std::array buffer; + std::string result = ""; + + // Use popen to execute the command and read its output + // "r" mode opens the pipe for reading + std::unique_ptr + pipe(popen(command.c_str(), "r"), pclose); + if (!pipe) { + throw std::runtime_error("popen() failed!"); + } + + // Read the output line by line + while (fgets(buffer.data(), buffer.size(), pipe.get()) != nullptr) { + result += buffer.data(); + } + + // The stanard GNU/Linux md5sum output format is: "32-char-hash + // filename". We extract the 32-character hash part... + if (result.length() >= 32) { + return result.substr(0, 32); + } else { + throw std::runtime_error("Failed to parse md5sum output."); + } +} diff --git a/include/exputils.H b/include/exputils.H index 986aed037..c530e8e68 100644 --- a/include/exputils.H +++ b/include/exputils.H @@ -12,4 +12,8 @@ orthoCompute(const std::vector& tests); void orthoTest(const std::vector& tests, const std::string& classname, const std::string& indexname); +// Signature for md5sum function computation +std::string get_md5sum(const std::string&); + + #endif From 6acbe733953460856bca78231c63869d41234e20 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 23 Mar 2026 21:15:40 -0400 Subject: [PATCH 50/78] Added ptype support to the Cylindrical cache --- expui/BiorthBasis.cc | 153 ++++++++++++++++++++++++++++++------------- 1 file changed, 109 insertions(+), 44 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 70c251481..2f996cb14 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1598,19 +1598,7 @@ namespace BasisClasses // int cache_status = sl->read_cache(); - // Define the HighFive mapping for DTYPE - // - HighFive::EnumType disk_type({ - {"constant", DiskType::constant}, - {"gaussian", DiskType::gaussian}, - {"mn", DiskType::mn}, - {"exponential", DiskType::exponential}, - {"doubleexpon", DiskType::doubleexpon}, - {"diskbulge", DiskType::diskbulge}, - {"python", DiskType::python} - }); - - std::map disk_type_string = { + std::map disk_type = { {DiskType::constant, "constant"}, {DiskType::gaussian, "gaussian"}, {DiskType::mn, "mn"}, @@ -1620,7 +1608,6 @@ namespace BasisClasses {DiskType::python, "python"} }; - // Checking for cache consistency with current DiskType and Python module (if applicable) // if (cache_status == 1) { @@ -1634,8 +1621,8 @@ namespace BasisClasses // if (!file.hasAttribute("DiskType")) { if (myid==0) { - std::cout << "---- Cylindrical: DiskType attribute not found in cache file <" << cachename << ">. " - << "---- This may indicate an old cache file created before DiskType metadata was added. " + std::cout << "---- Cylindrical: DiskType attribute not found in cache file <" << cachename << ">. " << std::endl + << "---- This may indicate an old cache file created before DiskType metadata was added. " << std::endl << "---- We will continue...but consider remaking the cache to avoid confusion." << std::endl; } } else { @@ -1643,30 +1630,31 @@ namespace BasisClasses // auto read_attr = file.getAttribute("DiskType"); - DiskType loaded_dtype; + std::string loaded_dtype; read_attr.read(loaded_dtype); - if (loaded_dtype != DTYPE) { + DiskType disktype = dtlookup.at(loaded_dtype); + + if (disktype != DTYPE) { if (myid==0) { std::cout << "---- Cylindrical: DiskType for cache file <" << cachename << "> is <" - << disk_type_string[loaded_dtype] << ">, which does not match the requested DiskType <" - << disk_type_string[DTYPE] << ">. Forcing cache recomputation." << std::endl; + << loaded_dtype << ">, which does not match the requested DiskType <" + << dtype << ">. Forcing cache recomputation." << std::endl; } // Force cache recomputation cache_status = 0; } - else if (loaded_dtype == DiskType::python) { + else if (disktype == DiskType::python) { // Get the pyname attribute - auto read_attr = file.getAttribute("pyname"); - std::string loaded_pyname; - read_attr.read(loaded_pyname); + std::vector pyinfo; + auto read_attr = file.getAttribute("pythonDiskType"); + read_attr.read(pyinfo); - std::string loaded_md5, current_md5; - + std::string current_md5; + // Get the md5sum for requested Python module try { current_md5 = get_md5sum(pyname); - loaded_md5 = get_md5sum(loaded_pyname); } catch (const std::runtime_error& e) { std::cerr << "Error: " << e.what() << std::endl; } @@ -1675,16 +1663,80 @@ namespace BasisClasses // module and the loaded Python module used to create the // cache. If they do not match, force cache recomputation // to ensure consistency with the current Python module. - if (current_md5 != loaded_md5) { + if (current_md5 != pyinfo[1]) { if (myid==0) { std::cout << "---- Cylindrical: Python module for disk density has changed since cache creation." << std::endl << "---- Current module: <" << pyname << ">, md5sum: " << current_md5 << std::endl - << "---- Loaded module: <" << loaded_pyname << ">, md5sum: " << loaded_md5 << std::endl + << "---- Loaded module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl << "---- Forcing cache recomputation to ensure consistency with current Python module." << std::endl; } cache_status = 0; } } + + // Get the deproject attribute + // + read_attr = file.getAttribute("deproject"); + bool loaded_deproject; + read_attr.read(loaded_deproject); + + if (deproject != loaded_deproject) { + if (myid==0) { + std::cout << "---- Cylindrical: deproject flag for cache file <" << cachename << "> is <" + << std::boolalpha << loaded_deproject << std::noboolalpha + << ">, which does not match the requested deproject flag <" + << std::boolalpha << deproject << std::noboolalpha + << ">. Forcing cache recomputation." << std::endl; + } + // Force cache recomputation + cache_status = 0; + } + + if (cache_status == 1 and deproject) { + // Get the dmodel attribute + // + read_attr = file.getAttribute("dmodel"); + std::string loaded_dmodel; + read_attr.read(loaded_dmodel); + + if (loaded_dmodel != dmodel) { + if (myid==0) { + std::cout << "---- Cylindrical: dmodel for cache file <" << cachename << "> is <" + << loaded_dmodel << ">, which does not match the requested dmodel <" + << dmodel << ">. Forcing cache recomputation." << std::endl; + } + // Force cache recomputation + cache_status = 0; + } + + if (cache_status == 1 and dmodel == "python") { + // Get the Python info + // + std::vector pyinfo; + read_attr = file.getAttribute("pythonProjType"); + read_attr.read(pyinfo); + + std::string current_md5; + + // Get the md5sum for requested Python projection module + try { + current_md5 = get_md5sum(pyproj); + } catch (const std::runtime_error& e) { + std::cerr << "Error: " << e.what() << std::endl; + } + // Check that the md5sums match for the current Python projection + // + if (current_md5 != pyinfo[1]) { + if (myid==0) { + std::cout << "---- Cylindrical: Python module for deprojection has changed since cache creation." << std::endl + << "---- Current module: <" << pyproj << ">, md5sum: " << current_md5 << std::endl + << "---- Loaded module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl + << "---- Forcing cache recomputation to ensure consistency with current Python projection module." << std::endl; + } + cache_status = 0; + } + } + } } } catch (const HighFive::Exception& err) { @@ -1743,7 +1795,7 @@ namespace BasisClasses // Set DiskType. This is the functional form for the disk used to // condition the basis. // - try { // Check for map entry, will through if the + try { // Check for map entry, will throw if the DTYPE = dtlookup.at(dtype); // key is not in the map. if (myid==0) { // Report DiskType @@ -1875,30 +1927,43 @@ namespace BasisClasses // HighFive::File file(cachename, HighFive::File::ReadWrite); - // Create a Scalar DataSpace for the single value - HighFive::DataSpace space = HighFive::DataSpace::From(DTYPE); - - // Create the attribute using the explicit (Name, DataSpace, - // DataType) signature using the base class 'DataType' to - // match the expected signature - auto attr = file.createAttribute("DiskType", space, (HighFive::DataType)disk_type); - - // Write the actual enum value - attr.write(DTYPE); + file.createAttribute("DiskType", + HighFive::DataSpace::From(dtype)).write(dtype); - // Get the md5sum for the Python module + // Write the md5sum for the Python module if (DTYPE == DiskType::python) { try { - std::string md5_hash = get_md5sum(pyname); + std::vector pyinfo = {pyname, get_md5sum(pyname)}; file.createAttribute - ("pyname", - HighFive::DataSpace::From(pyname)).write(pyname); + ("pythonDiskType", HighFive::DataSpace::From(pyinfo)).write(pyinfo); } catch (const std::runtime_error& e) { std::cerr << "Error: " << e.what() << std::endl; std::cerr << "Can not write the md5 hash to HDF5" << std::endl; } } + + // Save the deprojection flag + file.createAttribute + ("deproject", + HighFive::DataSpace::From(deproject)).write(deproject); + + // Reopen the DataSpace for DMODEL since we need to write it as a string + if (deproject) { + file.createAttribute + ("ProjType", HighFive::DataSpace::From(dmodel)).write(dmodel); + if (PTYPE == DeprojType::python) { + try { + std::vector pyinfo = {pyproj, get_md5sum(pyproj)}; + file.createAttribute + ("pythonProjType", + HighFive::DataSpace::From(pyinfo)).write(pyinfo); + } catch (const std::runtime_error& e) { + std::cerr << "Error: " << e.what() << std::endl; + std::cerr << "Can not writine the md5 hash to HDF5" << std::endl; + } + } + } } catch (const HighFive::Exception& err) { std::cerr << err.what() << std::endl; std::cerr << "Error writing metadata to cache file <" << cachename From 8dcc0a81b6c4e6e1b0b1ce2346d3342cb4d4466c Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Mon, 23 Mar 2026 23:15:01 -0400 Subject: [PATCH 51/78] Clean up the unique_ptr using a custom deleter --- exputil/getmd5sum.cc | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/exputil/getmd5sum.cc b/exputil/getmd5sum.cc index b8ab2ca4a..f15b8a9ff 100644 --- a/exputil/getmd5sum.cc +++ b/exputil/getmd5sum.cc @@ -5,8 +5,20 @@ #include #include +// Define the custom deleter +struct PcloseDeleter { + void operator()(FILE* f) const { + if (f) pclose(f); + } +}; + std::string get_md5sum(const std::string& filename) { + // Check if md5sum is available on the system + if (std::system("which md5sum > /dev/null 2>&1") != 0) { + throw std::runtime_error("md5sum command not found. Please ensure it is installed and in your PATH."); + } + // Command to execute: md5sum std::string command = "md5sum " + filename; std::array buffer; @@ -14,8 +26,7 @@ std::string get_md5sum(const std::string& filename) // Use popen to execute the command and read its output // "r" mode opens the pipe for reading - std::unique_ptr - pipe(popen(command.c_str(), "r"), pclose); + std::unique_ptr pipe(popen(command.c_str(), "r")); if (!pipe) { throw std::runtime_error("popen() failed!"); } From 9d11213bb1f121bcc635431391daabb46507f8a1 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Tue, 24 Mar 2026 10:40:34 +0000 Subject: [PATCH 52/78] restore testED binary --- utils/Test/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/utils/Test/CMakeLists.txt b/utils/Test/CMakeLists.txt index 6382664cb..01570787d 100644 --- a/utils/Test/CMakeLists.txt +++ b/utils/Test/CMakeLists.txt @@ -1,6 +1,6 @@ set(bin_PROGRAMS testBarrier expyaml orthoTest testDeproj - testEmpDeproj testEmp) + testEmpDeproj testEmp testED) set(common_LINKLIB OpenMP::OpenMP_CXX MPI::MPI_CXX expui exputil @@ -38,7 +38,7 @@ add_executable(testDeproj testDeproject.cc CubicSpline.cc Deprojector.cc) add_executable(testEmpDeproj testEmpDeproj.cc CubicSpline.cc Deprojector.cc EmpDeproj.cc) add_executable(testEmp testEmp.cc EmpDeproj.cc) - +add_executable(testED testED.cc EmpDeproj.cc) foreach(program ${bin_PROGRAMS}) target_link_libraries(${program} ${common_LINKLIB}) From cc7f53e1b7cb031cbc42e69b3fcc4fe99bdf8cd0 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Mar 2026 10:37:12 -0400 Subject: [PATCH 53/78] Added standalone md5 hash code to avoid problems with Mac and Windows down the line --- .gitmodules | 3 +++ CMakeLists.txt | 1 + expui/BiorthBasis.cc | 11 +++++++---- exputil/CMakeLists.txt | 14 ++++++++------ exputil/getmd5sum.cc | 3 +++ extern/QuickDigest5 | 1 + utils/Test/CMakeLists.txt | 4 +++- 7 files changed, 26 insertions(+), 11 deletions(-) create mode 160000 extern/QuickDigest5 diff --git a/.gitmodules b/.gitmodules index 3775929f0..cae935cf6 100644 --- a/.gitmodules +++ b/.gitmodules @@ -8,3 +8,6 @@ [submodule "extern/HighFive"] path = extern/HighFive url = https://github.com/BlueBrain/HighFive.git +[submodule "extern/QuickDigest5"] + path = extern/QuickDigest5 + url = https://github.com/nthnn/QuickDigest5.git diff --git a/CMakeLists.txt b/CMakeLists.txt index e2d49dba5..9ed3fa60c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -259,6 +259,7 @@ execute_process( include_directories(${PROJECT_SOURCE_DIR}/extern/yaml-cpp/include) include_directories(${PROJECT_SOURCE_DIR}/extern/pybind11/include) +include_directories(${PROJECT_SOURCE_DIR}/extern/QuickDigest5/include) # Report to the user message("Configuring build for ${GIT_BRANCH}/${GIT_COMMIT} at ${COMPILE_TIME}") diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 2f996cb14..97d94e694 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1,5 +1,6 @@ #include +#include "quickdigest5.hpp" #include "YamlCheck.H" #include "EXPException.H" #include "BiorthBasis.H" @@ -1654,7 +1655,7 @@ namespace BasisClasses // Get the md5sum for requested Python module try { - current_md5 = get_md5sum(pyname); + current_md5 = QuickDigest5::fileToHash(pyname); } catch (const std::runtime_error& e) { std::cerr << "Error: " << e.what() << std::endl; } @@ -1720,7 +1721,7 @@ namespace BasisClasses // Get the md5sum for requested Python projection module try { - current_md5 = get_md5sum(pyproj); + current_md5 = QuickDigest5::fileToHash(pyproj); } catch (const std::runtime_error& e) { std::cerr << "Error: " << e.what() << std::endl; } @@ -1933,7 +1934,8 @@ namespace BasisClasses // Write the md5sum for the Python module if (DTYPE == DiskType::python) { try { - std::vector pyinfo = {pyname, get_md5sum(pyname)}; + std::string hash = QuickDigest5::fileToHash(pyname); + std::vector pyinfo = {pyname, hash}; file.createAttribute ("pythonDiskType", HighFive::DataSpace::From(pyinfo)).write(pyinfo); } catch (const std::runtime_error& e) { @@ -1954,7 +1956,8 @@ namespace BasisClasses if (PTYPE == DeprojType::python) { try { - std::vector pyinfo = {pyproj, get_md5sum(pyproj)}; + std::string hash = QuickDigest5::fileToHash(pyproj); + std::vector pyinfo = {pyproj, hash}; file.createAttribute ("pythonProjType", HighFive::DataSpace::From(pyinfo)).write(pyinfo); diff --git a/exputil/CMakeLists.txt b/exputil/CMakeLists.txt index 3ba43a308..7d70b5c02 100644 --- a/exputil/CMakeLists.txt +++ b/exputil/CMakeLists.txt @@ -42,19 +42,21 @@ set(SLEDGE_SRC sledge.f) set(PARTICLE_SRC Particle.cc ParticleReader.cc header.cc) set(CUDA_SRC cudaParticle.cu cudaSLGridMP2.cu) set(PYWRAP_SRC DiskDensityFunc.cc) +set(QD5_SRC ${PROJECT_SOURCE_DIR}/extern/QuickDigest5/src/quickdigest5.cpp) -set(exputil_SOURCES ${ODE_SRC} ${ROOT_SRC} ${QUAD_SRC} - ${RANDOM_SRC} ${UTIL_SRC} ${SPECFUNC_SRC} - ${PHASE_SRC} ${SYMP_SRC} ${INTERP_SRC} ${MASSMODEL_SRC} - ${ORBIT_SRC} ${BIORTH_SRC} ${POLY_SRC} ${GAUSS_SRC} - ${QPDISTF_SRC} ${BESSEL_SRC} ${OPTIMIZATION_SRC} - ${SLEDGE_SRC} ${PARTICLE_SRC} ${CUDA_SRC} ${PYWRAP_SRC}) + +set(exputil_SOURCES ${ODE_SRC} ${ROOT_SRC} ${QUAD_SRC} ${RANDOM_SRC} + ${UTIL_SRC} ${SPECFUNC_SRC} ${PHASE_SRC} ${SYMP_SRC} ${INTERP_SRC} + ${MASSMODEL_SRC} ${ORBIT_SRC} ${BIORTH_SRC} ${POLY_SRC} ${GAUSS_SRC} + ${QPDISTF_SRC} ${BESSEL_SRC} ${OPTIMIZATION_SRC} ${SLEDGE_SRC} + ${PARTICLE_SRC} ${CUDA_SRC} ${PYWRAP_SRC} ${QD5_SRC}) set(common_INCLUDE_DIRS $ $ $ ${CMAKE_BINARY_DIR} $ $ + $ ${DEP_INC} ${EIGEN3_INCLUDE_DIR} ${HDF5_INCLUDE_DIRS} ${FFTW_INCLUDE_DIRS}) diff --git a/exputil/getmd5sum.cc b/exputil/getmd5sum.cc index f15b8a9ff..bea4b3673 100644 --- a/exputil/getmd5sum.cc +++ b/exputil/getmd5sum.cc @@ -1,3 +1,6 @@ +// Initial reference implementation. Better to use QuickDigest5 which +// is currently a git submodule. + #include #include #include diff --git a/extern/QuickDigest5 b/extern/QuickDigest5 new file mode 160000 index 000000000..1d61aece0 --- /dev/null +++ b/extern/QuickDigest5 @@ -0,0 +1 @@ +Subproject commit 1d61aece0e023e4c41cf24547fa1ac2f8028219d diff --git a/utils/Test/CMakeLists.txt b/utils/Test/CMakeLists.txt index 123e6120d..1d006a1a5 100644 --- a/utils/Test/CMakeLists.txt +++ b/utils/Test/CMakeLists.txt @@ -1,6 +1,6 @@ set(bin_PROGRAMS testBarrier expyaml orthoTest testDeproj - testEmpDeproj testEmp) + testEmpDeproj testEmp testmd5) set(common_LINKLIB OpenMP::OpenMP_CXX MPI::MPI_CXX expui exputil yaml-cpp ${VTK_LIBRARIES}) @@ -11,6 +11,7 @@ endif() set(common_INCLUDE $ + $ $ $ ${CMAKE_BINARY_DIR} ${DEP_INC} @@ -37,6 +38,7 @@ add_executable(testDeproj testDeproject.cc CubicSpline.cc Deprojector.cc) add_executable(testEmpDeproj testEmpDeproj.cc CubicSpline.cc Deprojector.cc EmpDeproj.cc) add_executable(testEmp testEmp.cc EmpDeproj.cc) +add_executable(testmd5 testmd5.cc) foreach(program ${bin_PROGRAMS}) target_link_libraries(${program} ${common_LINKLIB}) From 11987e503f7dfacd9d0a569865c70e66bb0eb6a2 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Mar 2026 10:45:30 -0400 Subject: [PATCH 54/78] Added missing md5 hash test code --- utils/Test/testmd5.cc | 48 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 utils/Test/testmd5.cc diff --git a/utils/Test/testmd5.cc b/utils/Test/testmd5.cc new file mode 100644 index 000000000..d8895000b --- /dev/null +++ b/utils/Test/testmd5.cc @@ -0,0 +1,48 @@ +// This test verifies that QuickDigest5 correctly computes the MD5 +// hash of a file + +#include +#include + +#include "quickdigest5.hpp" + +#include "exputils.H" // For get_md5sum which uses the system's md5sum + // command for verification + +int main(int argc, char* argv[]) +{ + // Default to example.txt if no argument + std::string filePath = argc > 1 ? argv[1] : "example.txt"; + + // Check if the file exists before trying to hash it + std::filesystem::path p(filePath); + if (!std::filesystem::exists(p)) { + std::cout << "File <" << filePath << "> not found." << std::endl + << "Usage: " << argv[0] << " [file_path]" << std::endl + << "Defaulting to example.txt if it exists." << std::endl; + } + + // One-line method to get the hex digest of a file + std::string hash = QuickDigest5::fileToHash(filePath); + + if (!hash.empty()) { + std::cout << "MD5: " << hash << std::endl; + } else { + std::cerr << "Error: Could not process file." << std::endl; + } + + // System version of md5sum for comparison + try { + std::string systemHash = get_md5sum(filePath); + std::cout << "System md5sum: " << systemHash << std::endl; + if (hash == systemHash) { + std::cout << "Success: hashes match!" << std::endl; + } else { + std::cerr << "Error: hashes do not match!" << std::endl; + } + } catch (const std::exception& e) { + std::cerr << "Error computing system md5sum: " << e.what() << std::endl; + } + + return 0; +} From fc16ed95c5ba0d1689fd745f4295ba7c48411d25 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Mar 2026 11:14:39 -0400 Subject: [PATCH 55/78] Additional cache failure info for debug --- expui/BiorthBasis.cc | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 97d94e694..69cf5ab26 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1638,9 +1638,15 @@ namespace BasisClasses if (disktype != DTYPE) { if (myid==0) { - std::cout << "---- Cylindrical: DiskType for cache file <" << cachename << "> is <" - << loaded_dtype << ">, which does not match the requested DiskType <" - << dtype << ">. Forcing cache recomputation." << std::endl; + std::cout << "---- Cylindrical: DiskType for cache file <" + << cachename << "> is <" + << loaded_dtype << ">," << std::endl + << "which does not match the requested DiskType <" + << dtype << ">, disktype=" + << (int)disktype << ", DTYPE=" + << (int)DTYPE << std::endl + << "---- Cylindrical: forcing cache recomputation" + << std::endl; } // Force cache recomputation cache_status = 0; @@ -1669,7 +1675,7 @@ namespace BasisClasses std::cout << "---- Cylindrical: Python module for disk density has changed since cache creation." << std::endl << "---- Current module: <" << pyname << ">, md5sum: " << current_md5 << std::endl << "---- Loaded module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl - << "---- Forcing cache recomputation to ensure consistency with current Python module." << std::endl; + << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; } cache_status = 0; } @@ -1687,7 +1693,8 @@ namespace BasisClasses << std::boolalpha << loaded_deproject << std::noboolalpha << ">, which does not match the requested deproject flag <" << std::boolalpha << deproject << std::noboolalpha - << ">. Forcing cache recomputation." << std::endl; + << ">" << std::endl + << "---- Cylindrical: forcing cache recomputation" << std::endl; } // Force cache recomputation cache_status = 0; @@ -1704,7 +1711,8 @@ namespace BasisClasses if (myid==0) { std::cout << "---- Cylindrical: dmodel for cache file <" << cachename << "> is <" << loaded_dmodel << ">, which does not match the requested dmodel <" - << dmodel << ">. Forcing cache recomputation." << std::endl; + << dmodel << ">" << std::endl + << "---- Cylindirical: forcing cache recomputation" << std::endl; } // Force cache recomputation cache_status = 0; @@ -1732,7 +1740,7 @@ namespace BasisClasses std::cout << "---- Cylindrical: Python module for deprojection has changed since cache creation." << std::endl << "---- Current module: <" << pyproj << ">, md5sum: " << current_md5 << std::endl << "---- Loaded module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl - << "---- Forcing cache recomputation to ensure consistency with current Python projection module." << std::endl; + << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; } cache_status = 0; } @@ -1742,7 +1750,7 @@ namespace BasisClasses } catch (const HighFive::Exception& err) { std::cerr << "---- BiorthBasis inconsistency in Cylindrical cache: " << err.what() << std::endl; - std::cerr << "---- Forcing cache recomputation." << std::endl; + std::cerr << "---- Cylindrical: forcing cache recomputation" << std::endl; cache_status = 0; // Fallback... } } From 31925de857c6a49feffc4456d90d6a880680c3f9 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Mar 2026 11:20:36 -0400 Subject: [PATCH 56/78] Need to define DTYPE before cache check --- expui/BiorthBasis.cc | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 69cf5ab26..2e182d6b2 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1609,7 +1609,17 @@ namespace BasisClasses {DiskType::python, "python"} }; - // Checking for cache consistency with current DiskType and Python module (if applicable) + // Convert dtype string to lower case + // + std::transform(dtype.begin(), dtype.end(), dtype.begin(), + [](unsigned char c){ return std::tolower(c); }); + + + // Check for map entry, will throw if the key is not in the map. + DTYPE = dtlookup.at(dtype); + + // Checking for cache consistency with current DiskType and Python + // module (if applicable) // if (cache_status == 1) { @@ -1796,17 +1806,10 @@ namespace BasisClasses throw std::runtime_error("Cylindrical:initialize: EmpCylSL bad parameter"); } - // Convert dtype string to lower case - // - std::transform(dtype.begin(), dtype.end(), dtype.begin(), - [](unsigned char c){ return std::tolower(c); }); - // Set DiskType. This is the functional form for the disk used to // condition the basis. // - try { // Check for map entry, will throw if the - DTYPE = dtlookup.at(dtype); // key is not in the map. - + try { if (myid==0) { // Report DiskType std::cout << "---- DiskType is <" << dtype << ">" << std::endl; From f40f1ed89bcfd83df20bc8ebcdbf1f54663e9ab8 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Mar 2026 13:01:39 -0400 Subject: [PATCH 57/78] Updates on the n-body side to handle conditioning types --- expui/BiorthBasis.cc | 103 +++++++++------- src/CMakeLists.txt | 1 + src/Cylinder.H | 51 +++++++- src/Cylinder.cc | 286 ++++++++++++++++++++++++++++++++++++++++--- 4 files changed, 380 insertions(+), 61 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 2e182d6b2..3e56af745 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1673,7 +1673,8 @@ namespace BasisClasses try { current_md5 = QuickDigest5::fileToHash(pyname); } catch (const std::runtime_error& e) { - std::cerr << "Error: " << e.what() << std::endl; + if (myid==0) + std::cerr << "Error: " << e.what() << std::endl; } // Check that the md5sums match for the current Python @@ -1722,7 +1723,7 @@ namespace BasisClasses std::cout << "---- Cylindrical: dmodel for cache file <" << cachename << "> is <" << loaded_dmodel << ">, which does not match the requested dmodel <" << dmodel << ">" << std::endl - << "---- Cylindirical: forcing cache recomputation" << std::endl; + << "---- Cylindrical: forcing cache recomputation" << std::endl; } // Force cache recomputation cache_status = 0; @@ -1741,7 +1742,8 @@ namespace BasisClasses try { current_md5 = QuickDigest5::fileToHash(pyproj); } catch (const std::runtime_error& e) { - std::cerr << "Error: " << e.what() << std::endl; + if (myid==0) + std::cerr << "Error: " << e.what() << std::endl; } // Check that the md5sums match for the current Python projection // @@ -1759,8 +1761,10 @@ namespace BasisClasses } } catch (const HighFive::Exception& err) { - std::cerr << "---- BiorthBasis inconsistency in Cylindrical cache: " << err.what() << std::endl; - std::cerr << "---- Cylindrical: forcing cache recomputation" << std::endl; + if (myid==0) { + std::cerr << "---- Cylindrical: " << err.what() << std::endl; + std::cerr << "---- Cylindrical: forcing cache recomputation" << std::endl; + } cache_status = 0; // Fallback... } } @@ -1934,55 +1938,66 @@ namespace BasisClasses sl->generate_eof(rnum, pnum, tnum, f); - try { - // Open the cache file for writing the Python metadata - // - HighFive::File file(cachename, HighFive::File::ReadWrite); - - file.createAttribute("DiskType", - HighFive::DataSpace::From(dtype)).write(dtype); - - // Write the md5sum for the Python module - if (DTYPE == DiskType::python) { - try { - std::string hash = QuickDigest5::fileToHash(pyname); - std::vector pyinfo = {pyname, hash}; - file.createAttribute - ("pythonDiskType", HighFive::DataSpace::From(pyinfo)).write(pyinfo); - } catch (const std::runtime_error& e) { - std::cerr << "Error: " << e.what() << std::endl; - std::cerr << "Can not write the md5 hash to HDF5" << std::endl; - } - } - - // Save the deprojection flag - file.createAttribute - ("deproject", - HighFive::DataSpace::From(deproject)).write(deproject); + if (myid==0) { + try { + // Open the cache file for writing the Python metadata + // + HighFive::File file(cachename, HighFive::File::ReadWrite); - // Reopen the DataSpace for DMODEL since we need to write it as a string - if (deproject) { - file.createAttribute - ("ProjType", HighFive::DataSpace::From(dmodel)).write(dmodel); + file.createAttribute("DiskType", + HighFive::DataSpace::From(dtype)).write(dtype); - if (PTYPE == DeprojType::python) { + // Write the md5sum for the Python module + if (DTYPE == DiskType::python) { try { - std::string hash = QuickDigest5::fileToHash(pyproj); - std::vector pyinfo = {pyproj, hash}; + std::vector pyinfo = + {pyname, QuickDigest5::fileToHash(pyname)}; file.createAttribute - ("pythonProjType", + ("pythonDiskType", HighFive::DataSpace::From(pyinfo)).write(pyinfo); } catch (const std::runtime_error& e) { - std::cerr << "Error: " << e.what() << std::endl; - std::cerr << "Can not writine the md5 hash to HDF5" << std::endl; + if (myid==0) { + std::cerr << "Error: " << e.what() << std::endl; + std::cerr << "Can not write the md5 hash to HDF5" << std::endl; + } } } + + // Save the deprojection flag + file.createAttribute + ("deproject", + HighFive::DataSpace::From(deproject)).write(deproject); + + // Reopen the DataSpace for DMODEL since we need to write it as a string + if (deproject) { + file.createAttribute + ("ProjType", HighFive::DataSpace::From(dmodel)).write(dmodel); + + if (PTYPE == DeprojType::python) { + try { + std::vector pyinfo = + {pyproj, QuickDigest5::fileToHash(pyproj)}; + file.createAttribute + ("pythonProjType", + HighFive::DataSpace::From(pyinfo)).write(pyinfo); + } catch (const std::runtime_error& e) { + if (myid==0) { + std::cerr << "Error: " << e.what() << std::endl; + std::cerr << "Can not writine the md5 hash to HDF5" << std::endl; + } + } + } + } + } catch (const HighFive::Exception& err) { + if (myid==0) { + std::cerr << err.what() << std::endl; + std::cerr << "Error writing metadata to cache file <" << cachename + << std::endl; + } } - } catch (const HighFive::Exception& err) { - std::cerr << err.what() << std::endl; - std::cerr << "Error writing metadata to cache file <" << cachename - << std::endl; + // Errors will prevent metadata from being written to the cache } + // Only the root process should be updating the cache } // Orthogonality sanity check diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index fa7ea2057..87df43c8f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -25,6 +25,7 @@ set(common_INCLUDE_DIRS $ $ $ + $ ${CMAKE_BINARY_DIR} ${DEP_INC} ${CMAKE_CURRENT_SOURCE_DIR} ${EIGEN3_INCLUDE_DIR}) diff --git a/src/Cylinder.H b/src/Cylinder.H index 5cadda2f8..49faa5de6 100644 --- a/src/Cylinder.H +++ b/src/Cylinder.H @@ -7,6 +7,7 @@ #include "Basis.H" #include "CylEXP.H" #include "Coefficients.H" +#include "DiskDensityFunc.H" #include "config_exp.h" @@ -37,6 +38,12 @@ class MixtureBasis; @param hcyl is the scale height + @param aratio is the scale length ratio for disk basis construction (default: 1.0) + + @param hratio is the scale height ratio for disk basis construction (default: 1.0) + + @param dweight is the relative weight of the second component in the double-exponential disk basis construction (default: 1.0) + @param lmaxfid is the maximum spherical harmonic index for EOF construction @param nmaxfid is the maximum radial index for EOF construction @@ -150,8 +157,11 @@ private: void initialize(void); - // Parameters + // Parameters + // double rcylmin, rcylmax, zmax, acyl; + double aratio, hratio, dweight; // Disk basis construction parameters for double-exponential disk + double Mfac, HERNA; // Disk bulge parameters: mass fraction and Hernquist scale length int nmaxfid, lmaxfid, mmax, mlim, nint; int ncylnx, ncylny, ncylr; double hcyl, hexp, snr, rem; @@ -275,7 +285,6 @@ protected: //@} - #endif /** Test change level counts for deep debugging enabled by setting @@ -336,6 +345,44 @@ protected: //! Write basis-specific parameters to HDF5 covariance file void writeCovarH5Params(HighFive::File& file); + enum class DiskType + { constant, gaussian, mn, exponential, doubleexpon, diskbulge, python }; + + //@{ + //! Disk target types + const std::map dtlookup = + { {"constant", DiskType::constant}, + {"gaussian", DiskType::gaussian}, + {"mn", DiskType::mn}, + {"exponential", DiskType::exponential}, + {"doubleexpon", DiskType::doubleexpon}, + {"diskbulge", DiskType::diskbulge}, + {"python", DiskType::python} + }; + + const std::map dtstring = + { {DiskType::constant, "constant"}, + {DiskType::gaussian, "gaussian"}, + {DiskType::mn, "mn"}, + {DiskType::exponential, "exponential"}, + {DiskType::doubleexpon, "doubleexpon"}, + {DiskType::diskbulge, "diskbulge"}, + {DiskType::python, "python"} + }; + + DiskType DTYPE = DiskType::exponential; + + std::string dtype; + + //! Dtype cache check for consistency between cache and current parameters + bool checkDtype(); + void saveDtype(); + + double DiskDens(double R, double z, double phi); + std::shared_ptr pyDens; + //@} + + public: //! Mutexes for multithreading diff --git a/src/Cylinder.cc b/src/Cylinder.cc index 9c7f84203..8fe9be50b 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -11,9 +11,10 @@ #include "Cylinder.H" #include "MixtureBasis.H" #include "DiskDensityFunc.H" -#include "Timer.H" -#include "exputils.H" -#include "NVTX.H" +#include "Timer.H" // for timing code sections +#include "exputils.H" // utility functions +#include "NVTX.H" // for NVTX profiling of CUDA code +#include "quickdigest5.hpp" // for md5 hashing of Python modules //@{ //! These are for testing exclusively (should be set false for production) @@ -31,6 +32,11 @@ Cylinder::valid_keys = { "hcyl", "sech2", "hexp", + "aratio", + "hratio", + "dweight", + "Mfac", + "HERNA", "snr", "evcut", "nmaxfid", @@ -74,6 +80,7 @@ Cylinder::valid_keys = { "playback", "coefCompute", "coefMaster", + "dtype", "pyname", "dumpbasis", "fullCovar", @@ -116,6 +123,16 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : ncylodd = 9; ncylrecomp = -1; + // For disk basis construction with doubleexpon + // + aratio = 1.0; // Radial scale length ratio + hratio = 1.0; // Vertical scale height ratio + dweight = 1.0; // Ratio of second disk relative to the first disk + + // For disk bulge + Mfac = 1.0; // mass fraction for disk + HERNA = 0.10; // Hernquist scale a disk bulge + rnum = 200; pnum = 1; tnum = 80; @@ -182,6 +199,14 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : ortho = std::make_shared (nmaxfid, lmaxfid, mmax, nmax, acyl, hcyl, ncylodd, cachename); + // Convert dtype string to lower case + // + std::transform(dtype.begin(), dtype.end(), dtype.begin(), + [](unsigned char c){ return std::tolower(c); }); + + // Check for map entry, will throw if the key is not in the map. + DTYPE = dtlookup.at(dtype); + // Set azimuthal harmonic order restriction? // if (mlim>=0) ortho->set_mlim(mlim); @@ -220,6 +245,10 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : if (std::filesystem::exists(cachename)) { cache_ok = ortho->read_cache(); + + if (cache_ok) { + cache_ok = checkDtype(); + } // If new cache is requested, backup existing basis file // @@ -269,19 +298,44 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : << "this step will take many minutes." << std::endl; - std::shared_ptr dfunc; - if (pyname.size()) dfunc = std::make_shared(pyname); - // Use either the user-supplied Python target or the default - // exponential disk model - auto DiskDens = [&](double R, double z, double phi) - { - if (dfunc) return (*dfunc)(R, z, phi); - double h = sech2 ? 0.5*hcyl : hcyl; - double f = exp(-fabs(z)/h); // Overflow prevention - double s = 2.0*f/(1.0 + f*f); // in sech computation - return exp(-R/acyl)*s*s/(4.0*M_PI*acyl*acyl*h); - }; + // Set DiskType. This is the functional form for the disk used to + // condition the basis. + // + try { + if (myid==0) { // Report DiskType + std::cout << "---- DiskType is <" << dtype << ">" << std::endl; + + if (not sech2) { + switch (DTYPE) { + case DiskType::doubleexpon: + case DiskType::exponential: + case DiskType::diskbulge: + std::cout << "---- pyEXP assumes sech^2(z/(2h)) by default in v7.9.0 and later" << std::endl + << "---- Use the 'sech2: true' in your YAML config to use sech^2(z/(2h))" << std::endl + << "---- This warning will be removed in v7.10.0." << std::endl; + break; + default: + break; + } + } + } + } + catch (const std::out_of_range& err) { + if (myid==0) { + std::cout << "DiskType error in configuraton file" << std::endl; + std::cout << "Valid options are: "; + for (auto v : dtlookup) std::cout << v.first << " "; + std::cout << std::endl; + } + throw std::runtime_error("Cylindrical::initialize: invalid DiskType"); + } + + // Check for and initialize the Python density type + // + if (DTYPE == DiskType::python) { + pyDens = std::make_shared(pyname); + } // The conditioning function for the EOF with an optional shift // for M>0 @@ -310,6 +364,7 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : }; ortho->generate_eof(rnum, pnum, tnum, dcond); + saveDtype(); } firstime = false; @@ -444,6 +499,12 @@ void Cylinder::initialize() if (conf["eof_file" ]) cachename = conf["eof_file" ].as(); if (conf["samplesz" ]) defSampT = conf["samplesz" ].as(); + if (conf["aratio" ]) aratio = conf["aratio" ].as(); + if (conf["hratio" ]) hratio = conf["hratio" ].as(); + if (conf["dweight" ]) dweight = conf["dweight" ].as(); + if (conf["Mfac" ]) Mfac = conf["Mfac" ].as(); + if (conf["HERNA" ]) HERNA = conf["HERNA" ].as(); + if (conf["rnum" ]) rnum = conf["rnum" ].as(); if (conf["pnum" ]) pnum = conf["pnum" ].as(); if (conf["tnum" ]) tnum = conf["tnum" ].as(); @@ -463,6 +524,8 @@ void Cylinder::initialize() if (conf["cmapr" ]) cmapR = conf["cmapr" ].as(); if (conf["cmapz" ]) cmapZ = conf["cmapz" ].as(); if (conf["vflag" ]) vflag = conf["vflag" ].as(); + if (conf["dtype" ]) dtype = conf["dtype" ].as(); + if (conf["pyname" ]) pyname = conf["pyname" ].as(); // Deprecation warning if (not sech2 and not conf["pyname"]) { @@ -598,6 +661,83 @@ void Cylinder::initialize() } } +double Cylinder::DiskDens(double R, double z, double phi) +{ + double ans = 0.0; + + switch (DTYPE) { + + case DiskType::constant: + if (R < acyl && fabs(z) < hcyl) + ans = 1.0/(2.0*hcyl*M_PI*acyl*acyl); + break; + + case DiskType::gaussian: + if (fabs(z) < hcyl) + ans = 1.0/(2.0*hcyl*2.0*M_PI*acyl*acyl)* + exp(-R*R/(2.0*acyl*acyl)); + break; + + case DiskType::mn: + { + double Z2 = z*z + hcyl*hcyl; + double Z = sqrt(Z2); + double Q2 = (acyl + Z)*(acyl + Z); + ans = 0.25*hcyl*hcyl/M_PI*(acyl*R*R + (acyl + 3.0*Z)*Q2)/( pow(R*R + Q2, 2.5) * Z*Z2 ); + } + break; + + case DiskType::doubleexpon: + { + double a1 = acyl; + double a2 = acyl*aratio; + double h1 = hcyl; + double h2 = hcyl*hratio; + double w1 = 1.0/(1.0+dweight); + double w2 = dweight/(1.0+dweight); + + if (sech2) { h1 *= 0.5; h2 *= 0.5; } + + double f1 = cosh(z/h1); + double f2 = cosh(z/h2); + + ans = + w1*exp(-R/a1)/(4.0*M_PI*a1*a1*h1*f1*f1) + + w2*exp(-R/a2)/(4.0*M_PI*a2*a2*h2*f2*f2) ; + } + break; + + case DiskType::diskbulge: + { + double h = sech2 ? 0.5*hcyl : hcyl; + double f = cosh(z/h); + double rr = pow(pow(R, 2) + pow(z,2), 0.5); + double w1 = Mfac; + double w2 = 1.0 - Mfac; + double as = HERNA; + + ans = w1*exp(-R/acyl)/(4.0*M_PI*acyl*acyl*h*f*f) + + w2*pow(as, 4)/(2.0*M_PI*rr)*pow(rr+as,-3.0) ; + } + break; + case DiskType::python: + ans = (*pyDens)(R, z, phi); + break; + case DiskType::exponential: + default: + { + double h = sech2 ? 0.5*hcyl : hcyl; + double f = cosh(z/h); + ans = exp(-R/acyl)/(4.0*M_PI*acyl*acyl*h*f*f); + } + break; + } + + // if (rwidth>0.0) ans *= erf((rtrunc-R)/rwidth); + + return ans; +} + void Cylinder::get_acceleration_and_potential(Component* C) { nvTracerPtr tPtr; @@ -919,6 +1059,10 @@ void Cylinder::determine_coefficients_particles(void) bool cache_ok = false; if (try_cache || restart) { cache_ok = ortho->read_cache(); + + if (cache_ok) { + cache_ok = checkDtype(); + } // For a restart, cache must be read // otherwise, abort if (restart && !cache_ok) { @@ -1878,3 +2022,115 @@ void Cylinder::writeCovarH5Params(HighFive::File& file) file.createAttribute("hcyl", HighFive::DataSpace::From(hcyl)).write(hcyl); } + +bool Cylinder::checkDtype() +{ + bool cache_status = true; + + // Open the cache file for reading the Python metadata + // + HighFive::File file(cachename, HighFive::File::ReadOnly); + + // Sanity: check that the DiskType attribute exists + // + if (!file.hasAttribute("DiskType")) { + if (myid==0) { + std::cout << "---- Cylinder: DiskType attribute not found in cache file <" << cachename << ">. " << std::endl + << "---- This may indicate an old cache file created before DiskType metadata was added. " << std::endl + << "---- We will continue...but consider remaking the cache to avoid confusion." << std::endl; + } + } else { + // Open existing DiskType attribute + // + auto read_attr = file.getAttribute("DiskType"); + + std::string loaded_dtype; + read_attr.read(loaded_dtype); + + DiskType disktype = dtlookup.at(loaded_dtype); + + if (disktype != DTYPE) { + if (myid==0) { + std::cout << "---- Cylindrical: DiskType for cache file <" + << cachename << "> is <" + << loaded_dtype << ">," << std::endl + << "which does not match the requested DiskType <" + << dtype << ">, disktype=" + << (int)disktype << ", DTYPE=" + << (int)DTYPE << std::endl + << "---- Cylindrical: forcing cache recomputation" + << std::endl; + } + // Force cache recomputation + cache_status = false; + } + else if (disktype == DiskType::python) { + // Get the pyname attribute + std::vector pyinfo; + auto read_attr = file.getAttribute("pythonDiskType"); + read_attr.read(pyinfo); + + std::string current_md5; + + // Get the md5sum for requested Python module + try { + current_md5 = QuickDigest5::fileToHash(pyname); + } catch (const std::runtime_error& e) { + if (myid==0) + std::cerr << "Error: " << e.what() << std::endl; + } + + // Check that the md5sums match for the current Python + // module and the loaded Python module used to create the + // cache. If they do not match, force cache recomputation + // to ensure consistency with the current Python module. + if (current_md5 != pyinfo[1]) { + if (myid==0) { + std::cout << "---- Cylindrical: Python module for disk density has changed since cache creation." << std::endl + << "---- Current module: <" << pyname << ">, md5sum: " << current_md5 << std::endl + << "---- Loaded module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl + << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; + } + cache_status = false; + } + } + } + + return cache_status; +} + +void Cylinder::saveDtype() +{ + // Only one process needs to write the dtype metadata + if (myid) return; + + std::string dtype = dtstring.at(DTYPE); + + try { + // Open the cache file for writing the Python metadata + // + HighFive::File file(cachename, HighFive::File::ReadWrite); + + file.createAttribute + ("DiskType", + HighFive::DataSpace::From(dtype)).write(dtype); + + // Write the md5sum for the Python module + if (DTYPE == DiskType::python) { + try { + std::vector pyinfo = + {pyname, QuickDigest5::fileToHash(pyname)}; + file.createAttribute + ("pythonDiskType", + HighFive::DataSpace::From(pyinfo)).write(pyinfo); + } catch (const std::runtime_error& e) { + std::cerr << "Error: " << e.what() << std::endl; + std::cerr << "Can not write the md5 hash to HDF5" << std::endl; + } + } + } catch (const HighFive::Exception& err) { + std::cerr << err.what() << std::endl; + std::cerr << "Error writing metadata to cache file <" << cachename + << std::endl; + } +} From 8f524056e166349556dd961e3409347779828388 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Mar 2026 13:14:30 -0400 Subject: [PATCH 58/78] A few minor checks for md5 tests; requires a GNU/Linux system --- exputil/getmd5sum.cc | 4 ++-- utils/Test/testmd5.cc | 4 ++++ 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/exputil/getmd5sum.cc b/exputil/getmd5sum.cc index bea4b3673..b7c234175 100644 --- a/exputil/getmd5sum.cc +++ b/exputil/getmd5sum.cc @@ -1,5 +1,5 @@ -// Initial reference implementation. Better to use QuickDigest5 which -// is currently a git submodule. +// Initial reference implementation for testing. Production code uses +// QuickDigest5 as a git submodule. #include #include diff --git a/utils/Test/testmd5.cc b/utils/Test/testmd5.cc index d8895000b..7490a8d2d 100644 --- a/utils/Test/testmd5.cc +++ b/utils/Test/testmd5.cc @@ -33,6 +33,10 @@ int main(int argc, char* argv[]) // System version of md5sum for comparison try { + if (std::system("which md5sum > /dev/null 2>&1") != 0) { + std::cerr << "Warning: md5sum command not found. Skipping system md5sum comparison." << std::endl; + return 0; + } std::string systemHash = get_md5sum(filePath); std::cout << "System md5sum: " << systemHash << std::endl; if (hash == systemHash) { From f608704b0e56168f0a494cfd9377f2ad7d42f140 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Mar 2026 13:50:32 -0400 Subject: [PATCH 59/78] Add more HDF5 diagnostics --- expui/BiorthBasis.cc | 189 ++++++++++++++++++++++++++----------------- 1 file changed, 116 insertions(+), 73 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 3e56af745..1ff0488c1 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1595,6 +1595,17 @@ namespace BasisClasses if (oldcache) sl->AllowOldCache(); + // Temporary muting of HDF5 error messages for EOF cache reading + H5E_auto2_t old_func; + void* old_client_data; + H5Eget_auto2(H5E_DEFAULT, &old_func, &old_client_data); + + // Mute HDF5 error messages + H5Eset_auto2(H5E_DEFAULT, NULL, NULL); + + // Unmute (Restore) + // H5Eset_auto2(H5E_DEFAULT, old_func, old_client_data); + // Attempt to read EOF cache // int cache_status = sl->read_cache(); @@ -1652,9 +1663,7 @@ namespace BasisClasses << cachename << "> is <" << loaded_dtype << ">," << std::endl << "which does not match the requested DiskType <" - << dtype << ">, disktype=" - << (int)disktype << ", DTYPE=" - << (int)DTYPE << std::endl + << dtype << ">" << std::endl << "---- Cylindrical: forcing cache recomputation" << std::endl; } @@ -1662,99 +1671,130 @@ namespace BasisClasses cache_status = 0; } else if (disktype == DiskType::python) { - // Get the pyname attribute - std::vector pyinfo; - auto read_attr = file.getAttribute("pythonDiskType"); - read_attr.read(pyinfo); - - std::string current_md5; - // Get the md5sum for requested Python module - try { - current_md5 = QuickDigest5::fileToHash(pyname); - } catch (const std::runtime_error& e) { - if (myid==0) - std::cerr << "Error: " << e.what() << std::endl; - } - - // Check that the md5sums match for the current Python - // module and the loaded Python module used to create the - // cache. If they do not match, force cache recomputation - // to ensure consistency with the current Python module. - if (current_md5 != pyinfo[1]) { + if (file.hasAttribute("pythonDiskType") == false) { if (myid==0) { - std::cout << "---- Cylindrical: Python module for disk density has changed since cache creation." << std::endl - << "---- Current module: <" << pyname << ">, md5sum: " << current_md5 << std::endl - << "---- Loaded module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl - << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; + std::cout << "---- Cylindrical: pythonDiskType attribute not found in cache file <" << cachename << ">. " << std::endl; + std::cout << "---- This may indicate an old cache file created before pythonDiskType metadata was added. " << std::endl; } cache_status = 0; + } else { + // Get the pyname attribute + auto read_attr = file.getAttribute("pythonDiskType"); + + std::vector pyinfo; + read_attr.read(pyinfo); + + std::string current_md5; + + // Get the md5sum for requested Python module + try { + current_md5 = QuickDigest5::fileToHash(pyname); + } catch (const std::runtime_error& e) { + if (myid==0) + std::cerr << "Error: " << e.what() << std::endl; + } + + // Check that the md5sums match for the current Python + // module and the loaded Python module used to create the + // cache. If they do not match, force cache recomputation + // to ensure consistency with the current Python module. + if (current_md5 != pyinfo[1]) { + if (myid==0) { + std::cout << "---- Cylindrical: Python module for disk density has changed since cache creation." << std::endl + << "---- Current module: <" << pyname << ">, md5sum: " << current_md5 << std::endl + << "---- Loaded module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl + << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; + } + cache_status = 0; + } } } // Get the deproject attribute // - read_attr = file.getAttribute("deproject"); - bool loaded_deproject; - read_attr.read(loaded_deproject); - - if (deproject != loaded_deproject) { + if (!file.hasAttribute("deproject")) { if (myid==0) { - std::cout << "---- Cylindrical: deproject flag for cache file <" << cachename << "> is <" - << std::boolalpha << loaded_deproject << std::noboolalpha - << ">, which does not match the requested deproject flag <" - << std::boolalpha << deproject << std::noboolalpha - << ">" << std::endl - << "---- Cylindrical: forcing cache recomputation" << std::endl; + std::cout << "---- Cylindrical: deproject attribute not found in cache file <" << cachename << ">. " << std::endl; + std::cout << "---- This may indicate an old cache file created before deproject metadata was added. " << std::endl; } - // Force cache recomputation - cache_status = 0; - } - - if (cache_status == 1 and deproject) { - // Get the dmodel attribute - // - read_attr = file.getAttribute("dmodel"); - std::string loaded_dmodel; - read_attr.read(loaded_dmodel); - - if (loaded_dmodel != dmodel) { + cache_status = 0; + } else { + read_attr = file.getAttribute("deproject"); + bool loaded_deproject; + read_attr.read(loaded_deproject); + + if (deproject != loaded_deproject) { if (myid==0) { - std::cout << "---- Cylindrical: dmodel for cache file <" << cachename << "> is <" - << loaded_dmodel << ">, which does not match the requested dmodel <" - << dmodel << ">" << std::endl + std::cout << "---- Cylindrical: deproject flag for cache file <" << cachename << "> is <" + << std::boolalpha << loaded_deproject << std::noboolalpha + << ">, which does not match the requested deproject flag <" + << std::boolalpha << deproject << std::noboolalpha + << ">" << std::endl << "---- Cylindrical: forcing cache recomputation" << std::endl; } // Force cache recomputation - cache_status = 0; - } + cache_status = 0; + } - if (cache_status == 1 and dmodel == "python") { - // Get the Python info + // Now check for deproject flag + // + if (cache_status == 1 and deproject) { + + // Get the dmodel attribute // - std::vector pyinfo; - read_attr = file.getAttribute("pythonProjType"); - read_attr.read(pyinfo); + read_attr = file.getAttribute("dmodel"); + std::string loaded_dmodel; + read_attr.read(loaded_dmodel); - std::string current_md5; - - // Get the md5sum for requested Python projection module - try { - current_md5 = QuickDigest5::fileToHash(pyproj); - } catch (const std::runtime_error& e) { - if (myid==0) - std::cerr << "Error: " << e.what() << std::endl; + if (loaded_dmodel != dmodel) { + if (myid==0) { + std::cout << "---- Cylindrical: dmodel for cache file <" << cachename << "> is <" + << loaded_dmodel << ">, which does not match the requested dmodel <" + << dmodel << ">" << std::endl + << "---- Cylindrical: forcing cache recomputation" << std::endl; + } + // Force cache recomputation + cache_status = 0; } - // Check that the md5sums match for the current Python projection + } + + if (cache_status == 1 and dmodel == "python") { + // Get the Python info // - if (current_md5 != pyinfo[1]) { + if (!file.hasAttribute("pythonProjType")) { if (myid==0) { - std::cout << "---- Cylindrical: Python module for deprojection has changed since cache creation." << std::endl - << "---- Current module: <" << pyproj << ">, md5sum: " << current_md5 << std::endl - << "---- Loaded module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl - << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; + std::cout << "---- Cylindrical: pythonProjType attribute not found in cache file <" << cachename << ">. " << std::endl; + std::cout << "---- This may indicate an old cache file created before pythonProjType metadata was added. " << std::endl; } + cache_status = 0; + + } else { + read_attr = file.getAttribute("pythonProjType"); + std::vector pyinfo; + read_attr.read(pyinfo); + + std::string current_md5; + + // Get the md5sum for requested Python projection module + try { + current_md5 = QuickDigest5::fileToHash(pyproj); + } catch (const std::runtime_error& e) { + if (myid==0) + std::cerr << "Error: " << e.what() << std::endl; + } + // Check that the md5sums match for the current Python projection + // + if (current_md5 != pyinfo[1]) { + if (myid==0) { + std::cout << "---- Cylindrical: Python module for deprojection has changed since cache creation." << std::endl + << "---- Current module: <" << pyproj << ">, md5sum: " << current_md5 << std::endl + << "---- Loaded module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl + << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; + } + cache_status = 0; + } } } } @@ -2000,6 +2040,9 @@ namespace BasisClasses // Only the root process should be updating the cache } + // Unmute HDF5 error messages (Restore) + H5Eset_auto2(H5E_DEFAULT, old_func, old_client_data); + // Orthogonality sanity check // if (myid==0) orthoTest(); From e02eb789d0c0a961218621db0aebb5ace3c54b21 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Mar 2026 13:56:40 -0400 Subject: [PATCH 60/78] Fix mistake in storing pyname info vector --- expui/BiorthBasis.cc | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 1ff0488c1..2696f19a5 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -12,6 +12,8 @@ #include #endif +#define HDF5_QUIET + namespace BasisClasses { const std::set @@ -1596,6 +1598,7 @@ namespace BasisClasses // Temporary muting of HDF5 error messages for EOF cache reading +#ifdef HDF5_QUIET H5E_auto2_t old_func; void* old_client_data; H5Eget_auto2(H5E_DEFAULT, &old_func, &old_client_data); @@ -1605,6 +1608,7 @@ namespace BasisClasses // Unmute (Restore) // H5Eset_auto2(H5E_DEFAULT, old_func, old_client_data); +#endif // Attempt to read EOF cache // @@ -1992,7 +1996,7 @@ namespace BasisClasses try { std::vector pyinfo = {pyname, QuickDigest5::fileToHash(pyname)}; - file.createAttribute + file.createAttribute> ("pythonDiskType", HighFive::DataSpace::From(pyinfo)).write(pyinfo); } catch (const std::runtime_error& e) { @@ -2017,7 +2021,7 @@ namespace BasisClasses try { std::vector pyinfo = {pyproj, QuickDigest5::fileToHash(pyproj)}; - file.createAttribute + file.createAttribute> ("pythonProjType", HighFive::DataSpace::From(pyinfo)).write(pyinfo); } catch (const std::runtime_error& e) { @@ -2040,9 +2044,10 @@ namespace BasisClasses // Only the root process should be updating the cache } +#ifdef HDF5_QUIET // Unmute HDF5 error messages (Restore) - H5Eset_auto2(H5E_DEFAULT, old_func, old_client_data); - + H5Eset_auto2(H5E_DEFAULT, old_func, old_client_data); +#endif // Orthogonality sanity check // if (myid==0) orthoTest(); From 8cba2a72c4cfe2cb5f51232c54842f4679bd0ac7 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Mar 2026 14:13:39 -0400 Subject: [PATCH 61/78] Additional comments [no CI] --- expui/BiorthBasis.cc | 24 +++++++++++++++++++++--- src/Cylinder.cc | 24 ++++++++++++++++++------ 2 files changed, 39 insertions(+), 9 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 2696f19a5..811607469 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1646,6 +1646,11 @@ namespace BasisClasses // Sanity: check that the DiskType attribute exists // if (!file.hasAttribute("DiskType")) { + // If the DiskType attribute is missing, this may indicate + // an old cache file created before the DiskType metadata + // was added. We will continue with a warning, but trigger + // cache recomputation to ensure consistency with the + // current DiskType. if (myid==0) { std::cout << "---- Cylindrical: DiskType attribute not found in cache file <" << cachename << ">. " << std::endl << "---- This may indicate an old cache file created before DiskType metadata was added. " << std::endl @@ -1676,10 +1681,17 @@ namespace BasisClasses } else if (disktype == DiskType::python) { + // If the DiskType is python, we also need to check that + // the Python module used to create the cache matches the + // current Python module to ensure consistency. This tag + // should be present in the cache if it was created with a + // recent version of the code, but if it is missing we + // will trigger cache recomputation to ensure consistency + // with the current Python module. if (file.hasAttribute("pythonDiskType") == false) { if (myid==0) { std::cout << "---- Cylindrical: pythonDiskType attribute not found in cache file <" << cachename << ">. " << std::endl; - std::cout << "---- This may indicate an old cache file created before pythonDiskType metadata was added. " << std::endl; + std::cout << "---- Cylindrical: this may be a logic error, trigger recomputation." << std::endl; } cache_status = 0; } else { @@ -1718,9 +1730,13 @@ namespace BasisClasses // Get the deproject attribute // if (!file.hasAttribute("deproject")) { + // We should not be able to get here since the deproject + // attribute is required for cache creation, but if it is + // missing we will trigger cache recomputation to ensure + // consistency. if (myid==0) { std::cout << "---- Cylindrical: deproject attribute not found in cache file <" << cachename << ">. " << std::endl; - std::cout << "---- This may indicate an old cache file created before deproject metadata was added. " << std::endl; + std::cout << "---- Cylindrical: this may be a logic error, trigger recomputation." << std::endl; } cache_status = 0; } else { @@ -1767,9 +1783,11 @@ namespace BasisClasses // Get the Python info // if (!file.hasAttribute("pythonProjType")) { + // We should not be able to get here since the pythonProjType + // attribute is required for cache creation with the Python if (myid==0) { std::cout << "---- Cylindrical: pythonProjType attribute not found in cache file <" << cachename << ">. " << std::endl; - std::cout << "---- This may indicate an old cache file created before pythonProjType metadata was added. " << std::endl; + std::cout << "---- Cylindrical: this may be a logic error, trigger recomputation." << std::endl; } cache_status = 0; diff --git a/src/Cylinder.cc b/src/Cylinder.cc index 8fe9be50b..4f34c64a2 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -2055,9 +2055,7 @@ bool Cylinder::checkDtype() << cachename << "> is <" << loaded_dtype << ">," << std::endl << "which does not match the requested DiskType <" - << dtype << ">, disktype=" - << (int)disktype << ", DTYPE=" - << (int)DTYPE << std::endl + << dtype << ">" << std::endl << "---- Cylindrical: forcing cache recomputation" << std::endl; } @@ -2065,9 +2063,19 @@ bool Cylinder::checkDtype() cache_status = false; } else if (disktype == DiskType::python) { + // Sanity check: if DiskType is python, then the pythonDiskType + // attribute must exist + if (!file.hasAttribute("pythonDiskType")) { + if (myid==0) { + std::cout << "---- Cylinder: pythonDiskType attribute not found in cache file <" << cachename << ">. " << std::endl; + std::cout << "---- Cylindier: this may indicate a logic error. Forcing cache recomputation." << std::endl; + } + // Force cache recomputation + cache_status = false; + } + auto read_attr = file.getAttribute("pythonDiskType"); // Get the pyname attribute std::vector pyinfo; - auto read_attr = file.getAttribute("pythonDiskType"); read_attr.read(pyinfo); std::string current_md5; @@ -2101,7 +2109,7 @@ bool Cylinder::checkDtype() void Cylinder::saveDtype() { - // Only one process needs to write the dtype metadata + // Only one process must write the dtype metadata if (myid) return; std::string dtype = dtstring.at(DTYPE); @@ -2111,11 +2119,15 @@ void Cylinder::saveDtype() // HighFive::File file(cachename, HighFive::File::ReadWrite); + // Write the DiskType attribute (as a string for human readability) file.createAttribute ("DiskType", HighFive::DataSpace::From(dtype)).write(dtype); - // Write the md5sum for the Python module + // Write the md5sum for the Python module, if Python disk type is + // used. This will allow us to check for consistency between the + // Python module used to create the cache and the current Python + // module, and force cache recomputation if they do not match. if (DTYPE == DiskType::python) { try { std::vector pyinfo = From 20eb5af46dd820a18a366ce6338f78351c7271f1 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Mar 2026 14:26:24 -0400 Subject: [PATCH 62/78] Missing '.py' suffix in md5 file call --- expui/BiorthBasis.cc | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 811607469..a270941ae 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -12,6 +12,8 @@ #include #endif +// Suppress HDF5 diagonostic messages from base layer when using +// HighFive #define HDF5_QUIET namespace BasisClasses @@ -1705,7 +1707,7 @@ namespace BasisClasses // Get the md5sum for requested Python module try { - current_md5 = QuickDigest5::fileToHash(pyname); + current_md5 = QuickDigest5::fileToHash(pyname + ".py"); } catch (const std::runtime_error& e) { if (myid==0) std::cerr << "Error: " << e.what() << std::endl; @@ -1801,7 +1803,7 @@ namespace BasisClasses // Get the md5sum for requested Python projection module try { - current_md5 = QuickDigest5::fileToHash(pyproj); + current_md5 = QuickDigest5::fileToHash(pyproj + ".py"); } catch (const std::runtime_error& e) { if (myid==0) std::cerr << "Error: " << e.what() << std::endl; @@ -2013,7 +2015,7 @@ namespace BasisClasses if (DTYPE == DiskType::python) { try { std::vector pyinfo = - {pyname, QuickDigest5::fileToHash(pyname)}; + {pyname, QuickDigest5::fileToHash(pyname + ".py")}; file.createAttribute> ("pythonDiskType", HighFive::DataSpace::From(pyinfo)).write(pyinfo); @@ -2038,7 +2040,7 @@ namespace BasisClasses if (PTYPE == DeprojType::python) { try { std::vector pyinfo = - {pyproj, QuickDigest5::fileToHash(pyproj)}; + {pyproj, QuickDigest5::fileToHash(pyproj + ".py")}; file.createAttribute> ("pythonProjType", HighFive::DataSpace::From(pyinfo)).write(pyinfo); From d60f5e34cb597f32def0aa38b9112e68ad5085d1 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Mar 2026 14:37:41 -0400 Subject: [PATCH 63/78] Identify member function errors in Cylinder [no CI] --- src/Cylinder.cc | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/Cylinder.cc b/src/Cylinder.cc index 4f34c64a2..c395d751d 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -2035,7 +2035,7 @@ bool Cylinder::checkDtype() // if (!file.hasAttribute("DiskType")) { if (myid==0) { - std::cout << "---- Cylinder: DiskType attribute not found in cache file <" << cachename << ">. " << std::endl + std::cout << "---- Cylinder:checkDtype: DiskType attribute not found in cache file <" << cachename << ">. " << std::endl << "---- This may indicate an old cache file created before DiskType metadata was added. " << std::endl << "---- We will continue...but consider remaking the cache to avoid confusion." << std::endl; } @@ -2051,7 +2051,7 @@ bool Cylinder::checkDtype() if (disktype != DTYPE) { if (myid==0) { - std::cout << "---- Cylindrical: DiskType for cache file <" + std::cout << "---- Cylinder::checkDtype: DiskType for cache file <" << cachename << "> is <" << loaded_dtype << ">," << std::endl << "which does not match the requested DiskType <" @@ -2067,8 +2067,8 @@ bool Cylinder::checkDtype() // attribute must exist if (!file.hasAttribute("pythonDiskType")) { if (myid==0) { - std::cout << "---- Cylinder: pythonDiskType attribute not found in cache file <" << cachename << ">. " << std::endl; - std::cout << "---- Cylindier: this may indicate a logic error. Forcing cache recomputation." << std::endl; + std::cout << "---- Cylinder::checkDtype: pythonDiskType attribute not found in cache file <" << cachename << ">. " << std::endl; + std::cout << "---- Cylindier:checkDtype: this may indicate a logic error. Forcing cache recomputation." << std::endl; } // Force cache recomputation cache_status = false; @@ -2094,7 +2094,7 @@ bool Cylinder::checkDtype() // to ensure consistency with the current Python module. if (current_md5 != pyinfo[1]) { if (myid==0) { - std::cout << "---- Cylindrical: Python module for disk density has changed since cache creation." << std::endl + std::cout << "---- Cylinder::checkDtype: Python module for disk density has changed since cache creation." << std::endl << "---- Current module: <" << pyname << ">, md5sum: " << current_md5 << std::endl << "---- Loaded module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; @@ -2115,11 +2115,12 @@ void Cylinder::saveDtype() std::string dtype = dtstring.at(DTYPE); try { - // Open the cache file for writing the Python metadata + // Reopen the cache file for writing the Python metadata // HighFive::File file(cachename, HighFive::File::ReadWrite); // Write the DiskType attribute (as a string for human readability) + // file.createAttribute ("DiskType", HighFive::DataSpace::From(dtype)).write(dtype); @@ -2128,6 +2129,7 @@ void Cylinder::saveDtype() // used. This will allow us to check for consistency between the // Python module used to create the cache and the current Python // module, and force cache recomputation if they do not match. + // if (DTYPE == DiskType::python) { try { std::vector pyinfo = @@ -2136,13 +2138,13 @@ void Cylinder::saveDtype() ("pythonDiskType", HighFive::DataSpace::From(pyinfo)).write(pyinfo); } catch (const std::runtime_error& e) { - std::cerr << "Error: " << e.what() << std::endl; + std::cerr << "Cylinder::saveDtype error: " << e.what() << std::endl; std::cerr << "Can not write the md5 hash to HDF5" << std::endl; } } } catch (const HighFive::Exception& err) { std::cerr << err.what() << std::endl; - std::cerr << "Error writing metadata to cache file <" << cachename - << std::endl; + std::cerr << "Cylinder::saveDtype: error writing metadata to cache file <" + << cachename << std::endl; } } From 26a896c3962d6c2f38688d29f7be12dd12dd3e18 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Mar 2026 15:07:10 -0400 Subject: [PATCH 64/78] Added more diagnostics --- expui/BiorthBasis.cc | 23 ++++++++++++++++------- src/Cylinder.cc | 6 +++--- 2 files changed, 19 insertions(+), 10 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index a270941ae..20e4bc076 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -2011,14 +2011,21 @@ namespace BasisClasses file.createAttribute("DiskType", HighFive::DataSpace::From(dtype)).write(dtype); + std::cout << "---- Cylindrical: writing DiskType <" << dtype + << "> to cache file <" << cachename << ">" << std::endl; + // Write the md5sum for the Python module if (DTYPE == DiskType::python) { try { std::vector pyinfo = {pyname, QuickDigest5::fileToHash(pyname + ".py")}; - file.createAttribute> - ("pythonDiskType", - HighFive::DataSpace::From(pyinfo)).write(pyinfo); + + file.createAttribute("pythonDiskType", pyinfo); + + std::cout << "---- Cylindrical: writing pythonDiskType <" << pyname + ".py" + << "> to cache file <" << cachename << ">" << std::endl; + + } catch (const std::runtime_error& e) { if (myid==0) { std::cerr << "Error: " << e.what() << std::endl; @@ -2041,13 +2048,15 @@ namespace BasisClasses try { std::vector pyinfo = {pyproj, QuickDigest5::fileToHash(pyproj + ".py")}; - file.createAttribute> - ("pythonProjType", - HighFive::DataSpace::From(pyinfo)).write(pyinfo); + + file.createAttribute("pythonProjType", pyinfo); + + std::cout << "---- Cylindrical: writing pythonProjType <" << pyproj + ".py" + << "> to cache file <" << cachename << ">" << std::endl; } catch (const std::runtime_error& e) { if (myid==0) { std::cerr << "Error: " << e.what() << std::endl; - std::cerr << "Can not writine the md5 hash to HDF5" << std::endl; + std::cerr << "Can not write the md5 hash to HDF5" << std::endl; } } } diff --git a/src/Cylinder.cc b/src/Cylinder.cc index c395d751d..162a75b95 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -2134,9 +2134,9 @@ void Cylinder::saveDtype() try { std::vector pyinfo = {pyname, QuickDigest5::fileToHash(pyname)}; - file.createAttribute - ("pythonDiskType", - HighFive::DataSpace::From(pyinfo)).write(pyinfo); + + file.createAttribute("pythonDiskType", pyinfo); + } catch (const std::runtime_error& e) { std::cerr << "Cylinder::saveDtype error: " << e.what() << std::endl; std::cerr << "Can not write the md5 hash to HDF5" << std::endl; From 82062e97367bc8ebce9589c497369b6f706aac41 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Mar 2026 15:28:30 -0400 Subject: [PATCH 65/78] More comments and output update with debugging info --- expui/BiorthBasis.cc | 43 ++++++++++++++++++++++++++----------------- src/Cylinder.cc | 19 ++++++++++--------- 2 files changed, 36 insertions(+), 26 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 20e4bc076..c7849a913 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1,7 +1,7 @@ #include -#include "quickdigest5.hpp" -#include "YamlCheck.H" +#include "quickdigest5.hpp" // for md5 hashing of Python modules +#include "YamlCheck.H" // for YAML configuration checking #include "EXPException.H" #include "BiorthBasis.H" #include "DiskModels.H" @@ -13,8 +13,8 @@ #endif // Suppress HDF5 diagonostic messages from base layer when using -// HighFive -#define HDF5_QUIET +// HighFive. This should be enabled unless one is debugging. +// #define HDF5_QUIET namespace BasisClasses { @@ -1608,7 +1608,7 @@ namespace BasisClasses // Mute HDF5 error messages H5Eset_auto2(H5E_DEFAULT, NULL, NULL); - // Unmute (Restore) + // For unmute, use: // H5Eset_auto2(H5E_DEFAULT, old_func, old_client_data); #endif @@ -1710,18 +1710,21 @@ namespace BasisClasses current_md5 = QuickDigest5::fileToHash(pyname + ".py"); } catch (const std::runtime_error& e) { if (myid==0) - std::cerr << "Error: " << e.what() << std::endl; + std::cerr << "BiorthBasis::Cylindrical error: " + << e.what() << ", error compuing pyname md5sum" + << std::endl; } // Check that the md5sums match for the current Python - // module and the loaded Python module used to create the - // cache. If they do not match, force cache recomputation - // to ensure consistency with the current Python module. + // module and the Python module used to create the + // cache. If they do not match, force cache + // recomputation to ensure consistency with the current + // Python module. if (current_md5 != pyinfo[1]) { if (myid==0) { std::cout << "---- Cylindrical: Python module for disk density has changed since cache creation." << std::endl << "---- Current module: <" << pyname << ">, md5sum: " << current_md5 << std::endl - << "---- Loaded module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl + << "---- Cached module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; } cache_status = 0; @@ -1806,7 +1809,9 @@ namespace BasisClasses current_md5 = QuickDigest5::fileToHash(pyproj + ".py"); } catch (const std::runtime_error& e) { if (myid==0) - std::cerr << "Error: " << e.what() << std::endl; + std::cerr << "BiorthBasis::Cylindrical error: " + << e.what() << ", error computing pyproj md5sum" + << std::endl; } // Check that the md5sums match for the current Python projection // @@ -1814,7 +1819,7 @@ namespace BasisClasses if (myid==0) { std::cout << "---- Cylindrical: Python module for deprojection has changed since cache creation." << std::endl << "---- Current module: <" << pyproj << ">, md5sum: " << current_md5 << std::endl - << "---- Loaded module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl + << "---- Cached module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; } cache_status = 0; @@ -2028,8 +2033,10 @@ namespace BasisClasses } catch (const std::runtime_error& e) { if (myid==0) { - std::cerr << "Error: " << e.what() << std::endl; - std::cerr << "Can not write the md5 hash to HDF5" << std::endl; + std::cerr << "BiorthBasis::Cylindrical error: " + << e.what() + << ", can not write the pyname and md5 hash to HDF5" + << std::endl; } } } @@ -2055,8 +2062,10 @@ namespace BasisClasses << "> to cache file <" << cachename << ">" << std::endl; } catch (const std::runtime_error& e) { if (myid==0) { - std::cerr << "Error: " << e.what() << std::endl; - std::cerr << "Can not write the md5 hash to HDF5" << std::endl; + std::cerr << "BiorthBasis::Cylindrical error: " + << e.what() + << ", can not write the pyinfo and md5 hash to HDF5" + << std::endl; } } } @@ -2074,7 +2083,7 @@ namespace BasisClasses } #ifdef HDF5_QUIET - // Unmute HDF5 error messages (Restore) + // Unmute HDF5 error messages H5Eset_auto2(H5E_DEFAULT, old_func, old_client_data); #endif // Orthogonality sanity check diff --git a/src/Cylinder.cc b/src/Cylinder.cc index 162a75b95..9bcc81fa2 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -2080,7 +2080,7 @@ bool Cylinder::checkDtype() std::string current_md5; - // Get the md5sum for requested Python module + // Get the md5sum for requested Python module source file try { current_md5 = QuickDigest5::fileToHash(pyname); } catch (const std::runtime_error& e) { @@ -2088,10 +2088,10 @@ bool Cylinder::checkDtype() std::cerr << "Error: " << e.what() << std::endl; } - // Check that the md5sums match for the current Python - // module and the loaded Python module used to create the - // cache. If they do not match, force cache recomputation - // to ensure consistency with the current Python module. + // Check that the md5sums match for the current Python module + // source files and the loaded Python module used to create the + // cache. If they do not match, force cache recomputation to + // ensure consistency with the current Python module. if (current_md5 != pyinfo[1]) { if (myid==0) { std::cout << "---- Cylinder::checkDtype: Python module for disk density has changed since cache creation." << std::endl @@ -2125,10 +2125,11 @@ void Cylinder::saveDtype() ("DiskType", HighFive::DataSpace::From(dtype)).write(dtype); - // Write the md5sum for the Python module, if Python disk type is - // used. This will allow us to check for consistency between the - // Python module used to create the cache and the current Python - // module, and force cache recomputation if they do not match. + // Write the md5sum for the Python module source file, if Python + // disk type is used. This will allow us to check for consistency + // between the Python module used to create the cache and the + // current Python module, and force cache recomputation if they do + // not match. // if (DTYPE == DiskType::python) { try { From e641cc5ebf2afe863dc97f7277f708ea267c58e7 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Mar 2026 17:13:41 -0400 Subject: [PATCH 66/78] Wrong tag type on read --- expui/BiorthBasis.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 5f1a1aef4..fd40b0fd3 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1816,7 +1816,7 @@ namespace BasisClasses // Get the dmodel attribute // - read_attr = file.getAttribute("dmodel"); + read_attr = file.getAttribute("ProjType"); std::string loaded_dmodel; read_attr.read(loaded_dmodel); From b10c2effffbf547e38b6494a3ef76c204e28cd5f Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Mar 2026 17:26:54 -0400 Subject: [PATCH 67/78] Try to autodeduce the HighFive types for simplicity --- expui/BiorthBasis.cc | 10 +++------- src/Cylinder.cc | 5 ++--- 2 files changed, 5 insertions(+), 10 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index fd40b0fd3..81bac7fae 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -2061,8 +2061,7 @@ namespace BasisClasses // HighFive::File file(cachename, HighFive::File::ReadWrite); - file.createAttribute("DiskType", - HighFive::DataSpace::From(dtype)).write(dtype); + file.createAttribute("DiskType", dtype); std::cout << "---- Cylindrical: writing DiskType <" << dtype << "> to cache file <" << cachename << ">" << std::endl; @@ -2090,14 +2089,11 @@ namespace BasisClasses } // Save the deprojection flag - file.createAttribute - ("deproject", - HighFive::DataSpace::From(deproject)).write(deproject); + file.createAttribute("deproject", deproject); // Reopen the DataSpace for DMODEL since we need to write it as a string if (deproject) { - file.createAttribute - ("ProjType", HighFive::DataSpace::From(dmodel)).write(dmodel); + file.createAttribute("ProjType", dmodel); if (PTYPE == DeprojType::python) { try { diff --git a/src/Cylinder.cc b/src/Cylinder.cc index 42cfb6095..75a1b8339 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -2118,6 +2118,7 @@ bool Cylinder::checkDtype() cache_status = false; } auto read_attr = file.getAttribute("pythonDiskType"); + // Get the pyname attribute std::vector pyinfo; read_attr.read(pyinfo); @@ -2165,9 +2166,7 @@ void Cylinder::saveDtype() // Write the DiskType attribute (as a string for human readability) // - file.createAttribute - ("DiskType", - HighFive::DataSpace::From(dtype)).write(dtype); + file.createAttribute("DiskType", dtype); // Write the md5sum for the Python module source file, if Python // disk type is used. This will allow us to check for consistency From 8e6ed66ea33ecc09d732560276b6fd06b79f2e7d Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Mar 2026 17:44:03 -0400 Subject: [PATCH 68/78] More nested checking for Cylinder --- src/Cylinder.cc | 135 +++++++++++++++++++++++++++++------------------- 1 file changed, 81 insertions(+), 54 deletions(-) diff --git a/src/Cylinder.cc b/src/Cylinder.cc index 75a1b8339..039c59594 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -2069,6 +2069,8 @@ void Cylinder::writeCovarH5Params(HighFive::File& file) bool Cylinder::checkDtype() { + // All processes must check the dtype metadata to ensure consistency + // bool cache_status = true; // Open the cache file for reading the Python metadata @@ -2079,75 +2081,96 @@ bool Cylinder::checkDtype() // if (!file.hasAttribute("DiskType")) { if (myid==0) { - std::cout << "---- Cylinder:checkDtype: DiskType attribute not found in cache file <" << cachename << ">. " << std::endl - << "---- This may indicate an old cache file created before DiskType metadata was added. " << std::endl - << "---- We will continue...but consider remaking the cache to avoid confusion." << std::endl; + std::cout << "---- Cylinder:checkDtype: DiskType attribute not found in cache file <" << cachename << ">" << std::endl + << "---- This may indicate an old cache file created before DiskType metadata was added" << std::endl + << "---- We will continue...but consider remaking the cache to avoid confusion" << std::endl; } } else { - // Open existing DiskType attribute + // Check for existence of DiskType attribute // - auto read_attr = file.getAttribute("DiskType"); + if (!file.hasAttribute("DiskType")) { + if (myid==0) { + std::cout << "---- Cylinder::checkDtype: DiskType attribute not found in cache file <" << cachename << ">" << std::endl + << "---- This may indicate an old cache file created before DiskType metadata was added" << std::endl + << "---- We will continue...but consider remaking the cache to avoid confusion" << std::endl; + } + } else { + // Open existing DiskType attribute + // + auto read_attr = file.getAttribute("DiskType"); - std::string loaded_dtype; - read_attr.read(loaded_dtype); + std::string loaded_dtype; + read_attr.read(loaded_dtype); - DiskType disktype = dtlookup.at(loaded_dtype); + // Map the loaded dtype string to a DiskType enum value + // + DiskType disktype = dtlookup.at(loaded_dtype); - if (disktype != DTYPE) { - if (myid==0) { - std::cout << "---- Cylinder::checkDtype: DiskType for cache file <" - << cachename << "> is <" - << loaded_dtype << ">," << std::endl - << "which does not match the requested DiskType <" - << dtype << ">" << std::endl - << "---- Cylindrical: forcing cache recomputation" - << std::endl; - } - // Force cache recomputation - cache_status = false; - } - else if (disktype == DiskType::python) { - // Sanity check: if DiskType is python, then the pythonDiskType - // attribute must exist - if (!file.hasAttribute("pythonDiskType")) { + if (disktype != DTYPE) { if (myid==0) { - std::cout << "---- Cylinder::checkDtype: pythonDiskType attribute not found in cache file <" << cachename << ">. " << std::endl; - std::cout << "---- Cylindier:checkDtype: this may indicate a logic error. Forcing cache recomputation." << std::endl; + std::cout << "---- Cylinder::checkDtype: DiskType for cache file <" + << cachename << "> is <" + << loaded_dtype << ">," << std::endl + << "which does not match the requested DiskType <" + << dtype << ">" << std::endl + << "---- Cylindrical: forcing cache recomputation" + << std::endl; } // Force cache recomputation - cache_status = false; + cache_status = false; } - auto read_attr = file.getAttribute("pythonDiskType"); - - // Get the pyname attribute - std::vector pyinfo; - read_attr.read(pyinfo); + else if (disktype == DiskType::python) { + // Sanity check: if DiskType is python, then the + // pythonDiskType attribute must exist + // + if (!file.hasAttribute("pythonDiskType")) { + if (myid==0) { + std::cout << "---- Cylinder::checkDtype: pythonDiskType attribute not found in cache file <" << cachename << ">. " << std::endl; + std::cout << "---- Cylindier:checkDtype: this may indicate a logic error. Forcing cache recomputation." << std::endl; + } + // Force cache recomputation + cache_status = false; + } else { + auto read_attr = file.getAttribute("pythonDiskType"); + + // Get the pyname attribute + std::vector pyinfo; + read_attr.read(pyinfo); - std::string current_md5; - - // Get the md5sum for requested Python module source file - try { - current_md5 = QuickDigest5::fileToHash(pyname); - } catch (const std::runtime_error& e) { - if (myid==0) - std::cerr << "Error: " << e.what() << std::endl; - } + std::string current_md5; - // Check that the md5sums match for the current Python module - // source files and the loaded Python module used to create the - // cache. If they do not match, force cache recomputation to - // ensure consistency with the current Python module. - if (current_md5 != pyinfo[1]) { - if (myid==0) { - std::cout << "---- Cylinder::checkDtype: Python module for disk density has changed since cache creation." << std::endl - << "---- Current module: <" << pyname << ">, md5sum: " << current_md5 << std::endl - << "---- Loaded module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl - << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; + // Get the md5sum for requested Python module source file + try { + current_md5 = QuickDigest5::fileToHash(pyname); + } catch (const std::runtime_error& e) { + if (myid==0) + std::cerr << "Cylinder::CheckDtype error: " + << e.what() << std::endl; + } + + // Check that the md5sums match for the current Python + // module source files and the loaded Python module used to + // create the cache. If they do not match, force cache + // recomputation to ensure consistency with the current + // Python module. + // + if (current_md5 != pyinfo[1]) { + if (myid==0) { + std::cout << "---- Cylinder::checkDtype: Python module for disk density has changed since cache creation." << std::endl + << "---- Current module: <" << pyname << ">, md5sum: " << current_md5 << std::endl + << "---- Loaded module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl + << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; + } + cache_status = false; + } } - cache_status = false; + // End: have Python disk type, check md5 hashes } + // End: Python disk type check } + // End: DiskType attribute exists, check value } + // End: DiskType attribute check return cache_status; } @@ -2155,6 +2178,7 @@ bool Cylinder::checkDtype() void Cylinder::saveDtype() { // Only one process must write the dtype metadata + // if (myid) return; std::string dtype = dtstring.at(DTYPE); @@ -2168,7 +2192,7 @@ void Cylinder::saveDtype() // file.createAttribute("DiskType", dtype); - // Write the md5sum for the Python module source file, if Python + // Write the md5sum for the Python module source file if Python // disk type is used. This will allow us to check for consistency // between the Python module used to create the cache and the // current Python module, and force cache recomputation if they do @@ -2186,6 +2210,9 @@ void Cylinder::saveDtype() std::cerr << "Can not write the md5 hash to HDF5" << std::endl; } } + + // Deprojection metadata could be added here + } catch (const HighFive::Exception& err) { std::cerr << err.what() << std::endl; std::cerr << "Cylinder::saveDtype: error writing metadata to cache file <" From e40d6ad73b175c963cb5d34cc075826ea37e77f8 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Mar 2026 17:48:08 -0400 Subject: [PATCH 69/78] Remove duplicateed DiskType check --- src/Cylinder.cc | 129 ++++++++++++++++++++++-------------------------- 1 file changed, 60 insertions(+), 69 deletions(-) diff --git a/src/Cylinder.cc b/src/Cylinder.cc index 039c59594..baf1efe4d 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -2086,89 +2086,80 @@ bool Cylinder::checkDtype() << "---- We will continue...but consider remaking the cache to avoid confusion" << std::endl; } } else { - // Check for existence of DiskType attribute + // Open existing DiskType attribute // - if (!file.hasAttribute("DiskType")) { - if (myid==0) { - std::cout << "---- Cylinder::checkDtype: DiskType attribute not found in cache file <" << cachename << ">" << std::endl - << "---- This may indicate an old cache file created before DiskType metadata was added" << std::endl - << "---- We will continue...but consider remaking the cache to avoid confusion" << std::endl; - } - } else { - // Open existing DiskType attribute - // - auto read_attr = file.getAttribute("DiskType"); + auto read_attr = file.getAttribute("DiskType"); - std::string loaded_dtype; - read_attr.read(loaded_dtype); + std::string loaded_dtype; + read_attr.read(loaded_dtype); - // Map the loaded dtype string to a DiskType enum value - // - DiskType disktype = dtlookup.at(loaded_dtype); + // Map the loaded dtype string to a DiskType enum value + // + DiskType disktype = dtlookup.at(loaded_dtype); - if (disktype != DTYPE) { + if (disktype != DTYPE) { + if (myid==0) { + std::cout << "---- Cylinder::checkDtype: DiskType for cache file <" + << cachename << "> is <" + << loaded_dtype << ">," << std::endl + << "which does not match the requested DiskType <" + << dtype << ">" << std::endl + << "---- Cylindrical: forcing cache recomputation" + << std::endl; + } + // Force cache recomputation + cache_status = false; + } + else if (disktype == DiskType::python) { + // Sanity check: if DiskType is python, then the + // pythonDiskType attribute must exist + // + if (!file.hasAttribute("pythonDiskType")) { if (myid==0) { - std::cout << "---- Cylinder::checkDtype: DiskType for cache file <" - << cachename << "> is <" - << loaded_dtype << ">," << std::endl - << "which does not match the requested DiskType <" - << dtype << ">" << std::endl - << "---- Cylindrical: forcing cache recomputation" - << std::endl; + std::cout << "---- Cylinder::checkDtype: pythonDiskType attribute not found in cache file <" << cachename << ">. " << std::endl; + std::cout << "---- Cylindier:checkDtype: this may indicate a logic error. Forcing cache recomputation." << std::endl; } // Force cache recomputation - cache_status = false; - } - else if (disktype == DiskType::python) { - // Sanity check: if DiskType is python, then the - // pythonDiskType attribute must exist + cache_status = false; + } else { + auto read_attr = file.getAttribute("pythonDiskType"); + + // Get the pyname attribute + std::vector pyinfo; + read_attr.read(pyinfo); + + std::string current_md5; + + // Get the md5sum for requested Python module source file + try { + current_md5 = QuickDigest5::fileToHash(pyname); + } catch (const std::runtime_error& e) { + if (myid==0) + std::cerr << "Cylinder::CheckDtype error: " + << e.what() << std::endl; + } + + // Check that the md5sums match for the current Python + // module source files and the loaded Python module used to + // create the cache. If they do not match, force cache + // recomputation to ensure consistency with the current + // Python module. // - if (!file.hasAttribute("pythonDiskType")) { + if (current_md5 != pyinfo[1]) { if (myid==0) { - std::cout << "---- Cylinder::checkDtype: pythonDiskType attribute not found in cache file <" << cachename << ">. " << std::endl; - std::cout << "---- Cylindier:checkDtype: this may indicate a logic error. Forcing cache recomputation." << std::endl; + std::cout << "---- Cylinder::checkDtype: Python module for disk density has changed since cache creation." << std::endl + << "---- Current module: <" << pyname << ">, md5sum: " << current_md5 << std::endl + << "---- Loaded module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl + << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; } - // Force cache recomputation cache_status = false; - } else { - auto read_attr = file.getAttribute("pythonDiskType"); - - // Get the pyname attribute - std::vector pyinfo; - read_attr.read(pyinfo); - - std::string current_md5; - - // Get the md5sum for requested Python module source file - try { - current_md5 = QuickDigest5::fileToHash(pyname); - } catch (const std::runtime_error& e) { - if (myid==0) - std::cerr << "Cylinder::CheckDtype error: " - << e.what() << std::endl; - } - - // Check that the md5sums match for the current Python - // module source files and the loaded Python module used to - // create the cache. If they do not match, force cache - // recomputation to ensure consistency with the current - // Python module. - // - if (current_md5 != pyinfo[1]) { - if (myid==0) { - std::cout << "---- Cylinder::checkDtype: Python module for disk density has changed since cache creation." << std::endl - << "---- Current module: <" << pyname << ">, md5sum: " << current_md5 << std::endl - << "---- Loaded module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl - << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; - } - cache_status = false; - } } - // End: have Python disk type, check md5 hashes } - // End: Python disk type check + // End: have Python disk type, check md5 hashes } - // End: DiskType attribute exists, check value + // End: Python disk type check + + // Could add deprojection checks here in the future } // End: DiskType attribute check From 2f874d764690bf853dfd015c3c1b1904b4b42776 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Mar 2026 17:51:58 -0400 Subject: [PATCH 70/78] Additional comments only [No CI] --- src/Cylinder.H | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/Cylinder.H b/src/Cylinder.H index 5247196de..b40b2425c 100644 --- a/src/Cylinder.H +++ b/src/Cylinder.H @@ -361,11 +361,11 @@ protected: //! Write basis-specific parameters to HDF5 covariance file void writeCovarH5Params(HighFive::File& file); + //! Disk density function types for basis conditioning enum class DiskType { constant, gaussian, mn, exponential, doubleexpon, diskbulge, python }; - //@{ - //! Disk target types + //! Map human strings to disk target typenums const std::map dtlookup = { {"constant", DiskType::constant}, {"gaussian", DiskType::gaussian}, @@ -376,6 +376,7 @@ protected: {"python", DiskType::python} }; + //! Disk target type strings for reporting const std::map dtstring = { {DiskType::constant, "constant"}, {DiskType::gaussian, "gaussian"}, @@ -386,18 +387,25 @@ protected: {DiskType::python, "python"} }; + //! The default disk density function type for basis conditioning DiskType DTYPE = DiskType::exponential; + //! The string name of the disk density function type for basis conditioning std::string dtype; //! Dtype cache check for consistency between cache and current parameters bool checkDtype(); + + //! Reopen and write the current Dtype into the cache metadata void saveDtype(); + //! Disk density functions for basis conditioning double DiskDens(double R, double z, double phi); - std::shared_ptr pyDens; - //@} + //! If using a Python disk density function, this is the instance of + //! the wrapper class that allows us to call the Python function + //! from C++ + std::shared_ptr pyDens; public: From 82de0308d0c61c32b1f3d3fc0488e0b54fd9f359 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Mar 2026 20:14:19 -0400 Subject: [PATCH 71/78] Attempted to fix garbled deprojection logic in BiorthBasis and EmpCylSL --- expui/BiorthBasis.cc | 247 ++++++++++++++++---------------------- exputil/EmpCylSL.cc | 54 +++++---- include/EmpCylSL.H | 8 +- src/Cylinder.H | 4 +- src/Cylinder.cc | 199 ++++++++++++++++-------------- utils/ICs/check_coefs.cc | 45 +++---- utils/ICs/check_coefs2.cc | 36 +++--- utils/ICs/initial.cc | 44 +++---- 8 files changed, 313 insertions(+), 324 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 81bac7fae..ff35f86a4 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1228,7 +1228,6 @@ namespace BasisClasses "ashift", "expcond", "ignore", - "deproject", "logr", "dmodel", "ppow", @@ -1407,7 +1406,6 @@ namespace BasisClasses cachename = ".eof_cache_file"; oldcache = false; Ignore = false; - deproject = false; rnum = 200; pnum = 1; @@ -1497,7 +1495,6 @@ namespace BasisClasses if (conf["cmapr" ]) cmapR = conf["cmapr" ].as(); if (conf["cmapz" ]) cmapZ = conf["cmapz" ].as(); if (conf["ignore" ]) Ignore = conf["ignore" ].as(); - if (conf["deproject" ]) deproject = conf["deproject" ].as(); if (conf["dmodel" ]) dmodel = conf["dmodel" ].as(); if (conf["aratio" ]) aratio = conf["aratio" ].as(); @@ -1581,39 +1578,33 @@ namespace BasisClasses EmpCylSL::logarithmic = logarithmic; EmpCylSL::VFLAG = vflag; - // Set deprojected model type - // + // Set EmpCylSL mtype. This is the spherical function used to + // generate the EOF basis. If "deproject" is set, this will be + // auto-set to the same as the deprojection model. Otherwise, + // this can be set independently to allow for different spherical + // functions for the EOF basis + // Convert mtype string to lower case + // std::transform(mtype.begin(), mtype.end(), mtype.begin(), [](unsigned char c){ return std::tolower(c); }); // Set EmpCylSL mtype. This is the spherical function used to - // generate the EOF basis. If "deproject" is set, this will be - // overriden in EmpCylSL. + // generate the EOF basis. // - EmpCylSL::mtype = EmpCylSL::Exponential; // Default - if (mtype.compare("exponential")==0) { - EmpCylSL::mtype = EmpCylSL::Exponential; + auto itm = EmpCylSL::EmpModelMap.find(mtype); + + if (itm == EmpCylSL::EmpModelMap.end()) { if (myid==0) { - std::cout << "---- Cylindrical: using original exponential deprojected disk for EOF conditioning" << std::endl; - std::cout << "---- Cylindrical: consider using the exact, spherically deprojected exponential with 'mtype: ExpSphere'" << std::endl; + std::cout << "No EmpCylSL EmpModel named <" + << mtype << ">, valid types are: " + << "Exponential, ExpSphere, Gaussian, Plummer, Power, Deproject " + << "(not case sensitive)" << std::endl; } - } else if (mtype.compare("expsphere")==0) - EmpCylSL::mtype = EmpCylSL::ExpSphere; - else if (mtype.compare("gaussian")==0) - EmpCylSL::mtype = EmpCylSL::Gaussian; - else if (mtype.compare("plummer")==0) - EmpCylSL::mtype = EmpCylSL::Plummer; - else if (mtype.compare("power")==0) { - EmpCylSL::mtype = EmpCylSL::Power; - EmpCylSL::PPOW = ppow; - } else { - if (myid==0) std::cout << "No EmpCylSL EmpModel named <" - << mtype << ">, valid types are: " - << "Exponential, ExpSphere, Gaussian, Plummer, Power " - << "(not case sensitive)" << std::endl; throw std::runtime_error("Cylindrical:initialize: EmpCylSL bad parameter"); } + + EmpCylSL::mtype = itm->second; // Check for non-null cache file name. This must be specified // to prevent recomputation and unexpected behavior. @@ -1693,18 +1684,41 @@ namespace BasisClasses // HighFive::File file(cachename, HighFive::File::ReadOnly); + if (!file.hasAttribute("EmpModel")) { + // If the EmpModel attribute is missing, we have an old cache type + if (myid==0) { + std::cout << "---- Cylindrical: EmpModel attribute not found in cache file <" << cachename << ">. " << std::endl; + std::cout << "---- Cylindrical: this may indicate an old cache file created before EmpModel metadata was added. " << std::endl; + std::cout << "---- Cylindrical: we will continue...but consider remaking the cache to avoid confusion." << std::endl; + } + } else { + auto read_attr = file.getAttribute("EmpModel"); + + std::string empmodel; + read_attr.read(empmodel); + + if (empmodel != mtype) { + if (myid==0) { + std::cout << "---- Cylindrical: EmpModel for cache file <" + << cachename << "> is <" + << empmodel << ">," << std::endl + << "which does not match the requested EmpModel <" + << mtype << ">" << std::endl + << "---- Cylindrical: forcing cache recomputation" + << std::endl; + } + // Force cache recomputation + cache_status = 0; + } + } + // Sanity: check that the DiskType attribute exists // if (!file.hasAttribute("DiskType")) { - // If the DiskType attribute is missing, this may indicate - // an old cache file created before the DiskType metadata - // was added. We will continue with a warning, but trigger - // cache recomputation to ensure consistency with the - // current DiskType. if (myid==0) { std::cout << "---- Cylindrical: DiskType attribute not found in cache file <" << cachename << ">. " << std::endl - << "---- This may indicate an old cache file created before DiskType metadata was added. " << std::endl - << "---- We will continue...but consider remaking the cache to avoid confusion." << std::endl; + << "---- This is a logic error since the DiskType attribute is required for cache creation" << std::endl + << "---- We will trigger cache recomputation to be safe." << std::endl; } } else { // Open existing DiskType attribute @@ -1782,96 +1796,66 @@ namespace BasisClasses // Get the deproject attribute // - if (!file.hasAttribute("deproject")) { - // We should not be able to get here since the deproject - // attribute is required for cache creation, but if it is - // missing we will trigger cache recomputation to ensure - // consistency. - if (myid==0) { - std::cout << "---- Cylindrical: deproject attribute not found in cache file <" << cachename << ">. " << std::endl; - std::cout << "---- Cylindrical: this may be a logic error, trigger recomputation." << std::endl; - } - cache_status = 0; - } else { - read_attr = file.getAttribute("deproject"); - bool loaded_deproject; - read_attr.read(loaded_deproject); + if (cache_status == 1 and + EmpCylSL::mtype == EmpCylSL::EmpModel::Deproject) { + + // Get the dmodel attribute + // + read_attr = file.getAttribute("ProjType"); + std::string loaded_dmodel; + read_attr.read(loaded_dmodel); - if (deproject != loaded_deproject) { + if (loaded_dmodel != dmodel) { if (myid==0) { - std::cout << "---- Cylindrical: deproject flag for cache file <" << cachename << "> is <" - << std::boolalpha << loaded_deproject << std::noboolalpha - << ">, which does not match the requested deproject flag <" - << std::boolalpha << deproject << std::noboolalpha - << ">" << std::endl + std::cout << "---- Cylindrical: dmodel for cache file <" << cachename << "> is <" + << loaded_dmodel << ">, which does not match the requested dmodel <" + << dmodel << ">" << std::endl << "---- Cylindrical: forcing cache recomputation" << std::endl; } // Force cache recomputation - cache_status = 0; - } - - // Now check for deproject flag + cache_status = 0; + } + } + + if (cache_status == 1 and dmodel == "python") { + // Get the Python info // - if (cache_status == 1 and deproject) { - - // Get the dmodel attribute - // - read_attr = file.getAttribute("ProjType"); - std::string loaded_dmodel; - read_attr.read(loaded_dmodel); - - if (loaded_dmodel != dmodel) { - if (myid==0) { - std::cout << "---- Cylindrical: dmodel for cache file <" << cachename << "> is <" - << loaded_dmodel << ">, which does not match the requested dmodel <" - << dmodel << ">" << std::endl - << "---- Cylindrical: forcing cache recomputation" << std::endl; - } - // Force cache recomputation - cache_status = 0; + if (!file.hasAttribute("pythonProjType")) { + // We should not be able to get here since the pythonProjType + // attribute is required for cache creation with the Python + if (myid==0) { + std::cout << "---- Cylindrical: pythonProjType attribute not found in cache file <" << cachename << ">. " << std::endl; + std::cout << "---- Cylindrical: this may be a logic error, trigger recomputation." << std::endl; } - } + + cache_status = 0; + + } else { + read_attr = file.getAttribute("pythonProjType"); + std::vector pyinfo; + read_attr.read(pyinfo); + + std::string current_md5; - if (cache_status == 1 and dmodel == "python") { - // Get the Python info + // Get the md5sum for requested Python projection module + try { + current_md5 = QuickDigest5::fileToHash(pyproj + ".py"); + } catch (const std::runtime_error& e) { + if (myid==0) + std::cerr << "BiorthBasis::Cylindrical error: " + << e.what() << ", error computing pyproj md5sum" + << std::endl; + } + // Check that the md5sums match for the current Python projection // - if (!file.hasAttribute("pythonProjType")) { - // We should not be able to get here since the pythonProjType - // attribute is required for cache creation with the Python + if (current_md5 != pyinfo[1]) { if (myid==0) { - std::cout << "---- Cylindrical: pythonProjType attribute not found in cache file <" << cachename << ">. " << std::endl; - std::cout << "---- Cylindrical: this may be a logic error, trigger recomputation." << std::endl; + std::cout << "---- Cylindrical: Python module for deprojection has changed since cache creation." << std::endl + << "---- Current module: <" << pyproj << ">, md5sum: " << current_md5 << std::endl + << "---- Cached module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl + << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; } - cache_status = 0; - - } else { - read_attr = file.getAttribute("pythonProjType"); - std::vector pyinfo; - read_attr.read(pyinfo); - - std::string current_md5; - - // Get the md5sum for requested Python projection module - try { - current_md5 = QuickDigest5::fileToHash(pyproj + ".py"); - } catch (const std::runtime_error& e) { - if (myid==0) - std::cerr << "BiorthBasis::Cylindrical error: " - << e.what() << ", error computing pyproj md5sum" - << std::endl; - } - // Check that the md5sums match for the current Python projection - // - if (current_md5 != pyinfo[1]) { - if (myid==0) { - std::cout << "---- Cylindrical: Python module for deprojection has changed since cache creation." << std::endl - << "---- Current module: <" << pyproj << ">, md5sum: " << current_md5 << std::endl - << "---- Cached module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl - << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; - } - cache_status = 0; - } } } } @@ -1900,32 +1884,6 @@ namespace BasisClasses << "---- remove 'ignore' from your YAML configuration." << std::endl; } - // Convert mtype string to lower case - // - std::transform(mtype.begin(), mtype.end(), mtype.begin(), - [](unsigned char c){ return std::tolower(c); }); - - // Set EmpCylSL mtype. This is the spherical function used to - // generate the EOF basis. If "deproject" is set, this will be - // overriden in EmpCylSL. - // - EmpCylSL::mtype = EmpCylSL::Exponential; - if (mtype.compare("exponential")==0) - EmpCylSL::mtype = EmpCylSL::Exponential; - else if (mtype.compare("gaussian")==0) - EmpCylSL::mtype = EmpCylSL::Gaussian; - else if (mtype.compare("plummer")==0) - EmpCylSL::mtype = EmpCylSL::Plummer; - else if (mtype.compare("power")==0) { - EmpCylSL::mtype = EmpCylSL::Power; - EmpCylSL::PPOW = ppow; - } else { - if (myid==0) std::cout << "No EmpCylSL EmpModel named <" - << mtype << ">, valid types are: " - << "Exponential, Gaussian, Plummer, Power " - << "(not case sensitive)" << std::endl; - throw std::runtime_error("Cylindrical:initialize: EmpCylSL bad parameter"); - } // Set DiskType. This is the functional form for the disk used to // condition the basis. @@ -1967,7 +1925,7 @@ namespace BasisClasses // Use these user models to deproject for the EOF spherical basis // - if (deproject) { + if (EmpCylSL::mtype == EmpCylSL::EmpModel::Deproject) { // The scale in EmpCylSL is assumed to be 1 so we compute the // height relative to the length // @@ -2061,6 +2019,12 @@ namespace BasisClasses // HighFive::File file(cachename, HighFive::File::ReadWrite); + file.createAttribute("EmpModel", mtype); + + std::cout << "---- Cylindrical: writing EmpModel <" + << EmpCylSL::EmpModelLabs.at(EmpCylSL::mtype) + << "> to cache file <" << cachename << ">" << std::endl; + file.createAttribute("DiskType", dtype); std::cout << "---- Cylindrical: writing DiskType <" << dtype @@ -2088,11 +2052,8 @@ namespace BasisClasses } } - // Save the deprojection flag - file.createAttribute("deproject", deproject); - - // Reopen the DataSpace for DMODEL since we need to write it as a string - if (deproject) { + if (EmpCylSL::mtype == EmpCylSL::EmpModel::Deproject) { + file.createAttribute("ProjType", dmodel); if (PTYPE == DeprojType::python) { diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 721d2c5a5..7b69bdc08 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -70,17 +70,25 @@ double EmpCylSL::PPOW = 4.0; bool EmpCylSL::NewCoefs = true; -EmpCylSL::EmpModel EmpCylSL::mtype = Exponential; - -std::map EmpCylSL::EmpModelLabs = - { {Exponential, "Exponential"}, - {ExpSphere, "ExpSphere" }, - {Gaussian, "Gaussian" }, - {Plummer, "Plummer" }, - {Power, "Power" }, - {Deproject, "Deproject" } - }; +EmpCylSL::EmpModel EmpCylSL::mtype = EmpCylSL::EmpModel::Exponential; + +std::map EmpCylSL::EmpModelMap = { + {"exponential", EmpModel::Exponential}, + {"expsphere", EmpModel::ExpSphere}, + {"gaussian", EmpModel::Gaussian}, + {"plummer", EmpModel::Plummer}, + {"power", EmpModel::Power}, + {"deproject", EmpModel::Deproject} +}; +std::map EmpCylSL::EmpModelLabs = { + {EmpModel::Exponential, "Exponential"}, + {EmpModel::ExpSphere, "ExpSphere"}, + {EmpModel::Gaussian, "Gaussian"}, + {EmpModel::Plummer, "Plummer"}, + {EmpModel::Power, "Power"}, + {EmpModel::Deproject, "Deproject"} +}; EmpCylSL::EmpCylSL() { @@ -555,7 +563,7 @@ void EmpCylSL::create_deprojection(double H, double Rf, int NUMR, int NINT, // densRg = Linear1d(rl, rho); massRg = Linear1d(rl, mass); - mtype = Deproject; + mtype = EmpModel::Deproject; } @@ -568,21 +576,21 @@ double EmpCylSL::massR(double R) double ans=0.0, fac, arg; switch (mtype) { - case Exponential: + case EmpModel::Exponential: ans = 1.0 - (1.0 + R)*exp(-R); break; - case ExpSphere: + case EmpModel::ExpSphere: ans = expdeproj.mass(R); break; - case Gaussian: + case EmpModel::Gaussian: arg = 0.5*R*R; ans = 1.0 - exp(-arg); break; - case Plummer: + case EmpModel::Plummer: fac = R/(1.0+R); ans = pow(fac, 3.0); break; - case Power: + case EmpModel::Power: { double z = R + 1.0; double a1 = PPOW - 1.0; @@ -593,7 +601,7 @@ double EmpCylSL::massR(double R) (1.0 - pow(z, -a1))/a1 ); } break; - case Deproject: + case EmpModel::Deproject: if (R < RMIN) ans = 0.0; else if (R>=RMAX) ans = massRg.eval(RMAX); else ans = massRg.eval(log(R)); @@ -608,21 +616,21 @@ double EmpCylSL::densR(double R) double ans=0.0, fac, arg; switch (mtype) { - case Exponential: + case EmpModel::Exponential: ans = exp(-R)/(4.0*M_PI*R); break; - case ExpSphere: + case EmpModel::ExpSphere: ans = expdeproj.density(R); break; - case Gaussian: + case EmpModel::Gaussian: arg = 0.5*R*R; ans = exp(-arg)/(4.0*M_PI*R); break; - case Plummer: + case EmpModel::Plummer: fac = 1.0/(1.0+R); ans = 3.0*pow(fac, 4.0)/(4.0*M_PI); break; - case Power: + case EmpModel::Power: { double z = R + 1.0; double a1 = PPOW - 1.0; @@ -631,7 +639,7 @@ double EmpCylSL::densR(double R) ans = 0.125*a1*a2*a3/M_PI * pow(z, -PPOW); } break; - case Deproject: + case EmpModel::Deproject: if (R < RMIN) ans = densRg.eval(RMIN); else if (R>=RMAX) ans = 0.0; else ans = densRg.eval(log(R)); diff --git a/include/EmpCylSL.H b/include/EmpCylSL.H index 445fcf713..cbbe7f308 100644 --- a/include/EmpCylSL.H +++ b/include/EmpCylSL.H @@ -374,7 +374,7 @@ public: }; //! Type of density model to use - enum EmpModel { + enum class EmpModel { Exponential, ExpSphere, Gaussian, @@ -383,6 +383,9 @@ public: Deproject, }; + static std::map EmpModelMap; + static std::map EmpModelLabs; + //! Axisymmetric disk density function for deprojection class AxiDisk { @@ -488,9 +491,6 @@ public: //! Use YAML header in coefficient file static bool NewCoefs; - //! Convert EmpModel to ascii label - static std::map EmpModelLabs; - //! Fraction of table range for basis images (for debug) static double HFAC; diff --git a/src/Cylinder.H b/src/Cylinder.H index b40b2425c..653b673e1 100644 --- a/src/Cylinder.H +++ b/src/Cylinder.H @@ -394,10 +394,10 @@ protected: std::string dtype; //! Dtype cache check for consistency between cache and current parameters - bool checkDtype(); + bool checkMetaData(); //! Reopen and write the current Dtype into the cache metadata - void saveDtype(); + void saveMetaData(); //! Disk density functions for basis conditioning double DiskDens(double R, double z, double phi); diff --git a/src/Cylinder.cc b/src/Cylinder.cc index baf1efe4d..30f00604b 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -198,39 +198,28 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : if (cachename.size()==0) throw std::runtime_error("EmpCylSL: you must specify a cachename"); - - // Set deprojected model type - // // Convert mtype string to lower case + // std::transform(mtype.begin(), mtype.end(), mtype.begin(), [](unsigned char c){ return std::tolower(c); }); // Set EmpCylSL mtype. This is the spherical function used to - // generate the EOF basis. + // generate the EOF basis. If "deproject" is set, this will be + // overriden in EmpCylSL. // - EmpCylSL::mtype = EmpCylSL::Exponential; // Default - if (mtype.compare("exponential")==0) { - EmpCylSL::mtype = EmpCylSL::Exponential; + auto itm = EmpCylSL::EmpModelMap.find(mtype); + + if (itm == EmpCylSL::EmpModelMap.end()) { if (myid==0) { - std::cout << "---- Cylindrical: using original exponential deprojected disk for EOF conditioning" << std::endl; - std::cout << "---- Cylindrical: consider using the exact, spherically deprojected exponential with 'mtype: ExpSphere'" << std::endl; + std::cout << "No EmpCylSL EmpModel named <" + << mtype << ">, valid types are: " + << "Exponential, ExpSphere, Gaussian, Plummer, Power, Deproject " + << "(not case sensitive)" << std::endl; } - } else if (mtype.compare("expsphere")==0) - EmpCylSL::mtype = EmpCylSL::ExpSphere; - else if (mtype.compare("gaussian")==0) - EmpCylSL::mtype = EmpCylSL::Gaussian; - else if (mtype.compare("plummer")==0) - EmpCylSL::mtype = EmpCylSL::Plummer; - else if (mtype.compare("power")==0) { - EmpCylSL::mtype = EmpCylSL::Power; - EmpCylSL::PPOW = ppow; - } else { - if (myid==0) std::cout << "No EmpCylSL EmpModel named <" - << mtype << ">, valid types are: " - << "Exponential, ExpSphere, Gaussian, Plummer, Power " - << "(not case sensitive)" << std::endl; - throw std::runtime_error("Cylinder:initialize: EmpCylSL bad parameter"); + throw std::runtime_error("Cylindrical:initialize: EmpCylSL bad parameter"); } + + EmpCylSL::mtype = itm->second; // Make the empirical orthogonal basis instance // @@ -285,7 +274,7 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : cache_ok = ortho->read_cache(); if (cache_ok) { - cache_ok = checkDtype(); + cache_ok = checkMetaData(); } // If new cache is requested, backup existing basis file @@ -402,7 +391,7 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : }; ortho->generate_eof(rnum, pnum, tnum, dcond); - saveDtype(); + saveMetaData(); } firstime = false; @@ -1105,7 +1094,7 @@ void Cylinder::determine_coefficients_particles(void) cache_ok = ortho->read_cache(); if (cache_ok) { - cache_ok = checkDtype(); + cache_ok = checkMetaData(); } // For a restart, cache must be read // otherwise, abort @@ -2067,7 +2056,7 @@ void Cylinder::writeCovarH5Params(HighFive::File& file) } -bool Cylinder::checkDtype() +bool Cylinder::checkMetaData() { // All processes must check the dtype metadata to ensure consistency // @@ -2077,96 +2066,130 @@ bool Cylinder::checkDtype() // HighFive::File file(cachename, HighFive::File::ReadOnly); - // Sanity: check that the DiskType attribute exists + // Sanity: check that the EmpModel attribute exists // - if (!file.hasAttribute("DiskType")) { + if (!file.hasAttribute("EmpModel")) { if (myid==0) { - std::cout << "---- Cylinder:checkDtype: DiskType attribute not found in cache file <" << cachename << ">" << std::endl - << "---- This may indicate an old cache file created before DiskType metadata was added" << std::endl + std::cout << "---- Cylinder:checkDtype: EmpModel attribute not found in cache file <" << cachename << ">" << std::endl + << "---- This may indicate an old cache file created before EmpModel metadata was added" << std::endl << "---- We will continue...but consider remaking the cache to avoid confusion" << std::endl; - } + } } else { - // Open existing DiskType attribute + // Check that the EmpModel attribute matches the current empirical model // - auto read_attr = file.getAttribute("DiskType"); - - std::string loaded_dtype; - read_attr.read(loaded_dtype); + auto read_attr = file.getAttribute("EmpModel"); + + std::string loaded_mtype; + read_attr.read(loaded_mtype); - // Map the loaded dtype string to a DiskType enum value - // - DiskType disktype = dtlookup.at(loaded_dtype); - - if (disktype != DTYPE) { + if (loaded_mtype != mtype) { if (myid==0) { - std::cout << "---- Cylinder::checkDtype: DiskType for cache file <" - << cachename << "> is <" - << loaded_dtype << ">," << std::endl - << "which does not match the requested DiskType <" - << dtype << ">" << std::endl + std::cout << "---- Cylinder::checkDtype: EmpModel for cache file <" + << cachename << "> is <" << loaded_mtype << ">," + << " which does not match the requested empirical model <" + << mtype << ">" << std::endl << "---- Cylindrical: forcing cache recomputation" << std::endl; } // Force cache recomputation cache_status = false; } - else if (disktype == DiskType::python) { - // Sanity check: if DiskType is python, then the - // pythonDiskType attribute must exist - // - if (!file.hasAttribute("pythonDiskType")) { + else { + + if (!file.hasAttribute("DiskType")) { if (myid==0) { - std::cout << "---- Cylinder::checkDtype: pythonDiskType attribute not found in cache file <" << cachename << ">. " << std::endl; - std::cout << "---- Cylindier:checkDtype: this may indicate a logic error. Forcing cache recomputation." << std::endl; + std::cout << "---- Cylinder:checkDtype: DiskType attribute not found in cache file <" << cachename << ">" << std::endl + << "---- Cylinder:checkDtype: this may indicate a logic error. Forcing cache recomputation." << std::endl; } // Force cache recomputation cache_status = false; - } else { - auto read_attr = file.getAttribute("pythonDiskType"); - - // Get the pyname attribute - std::vector pyinfo; - read_attr.read(pyinfo); + } + else { + // Open existing DiskType attribute + // + auto read_attr = file.getAttribute("DiskType"); + + std::string loaded_dtype; + read_attr.read(loaded_dtype); - std::string current_md5; + // Map the loaded dtype string to a DiskType enum value + // + DiskType disktype = dtlookup.at(loaded_dtype); - // Get the md5sum for requested Python module source file - try { - current_md5 = QuickDigest5::fileToHash(pyname); - } catch (const std::runtime_error& e) { - if (myid==0) - std::cerr << "Cylinder::CheckDtype error: " - << e.what() << std::endl; + if (disktype != DTYPE) { + if (myid==0) { + std::cout << "---- Cylinder::checkDtype: DiskType for cache file <" + << cachename << "> is <" + << loaded_dtype << ">," << std::endl + << "which does not match the requested DiskType <" + << dtype << ">" << std::endl + << "---- Cylindrical: forcing cache recomputation" + << std::endl; + } + // Force cache recomputation + cache_status = false; } + else if (disktype == DiskType::python) { + // Sanity check: if DiskType is python, then the + // pythonDiskType attribute must exist + // + if (!file.hasAttribute("pythonDiskType")) { + if (myid==0) { + std::cout << "---- Cylinder::checkDtype: pythonDiskType attribute not found in cache file <" << cachename << ">. " << std::endl; + std::cout << "---- Cylindier:checkDtype: this may indicate a logic error. Forcing cache recomputation." << std::endl; + } + // Force cache recomputation + cache_status = false; + } else { + auto read_attr = file.getAttribute("pythonDiskType"); - // Check that the md5sums match for the current Python - // module source files and the loaded Python module used to - // create the cache. If they do not match, force cache - // recomputation to ensure consistency with the current - // Python module. - // - if (current_md5 != pyinfo[1]) { - if (myid==0) { - std::cout << "---- Cylinder::checkDtype: Python module for disk density has changed since cache creation." << std::endl - << "---- Current module: <" << pyname << ">, md5sum: " << current_md5 << std::endl - << "---- Loaded module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl - << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; + // Get the pyname attribute + std::vector pyinfo; + read_attr.read(pyinfo); + + std::string current_md5; + + // Get the md5sum for requested Python module source file + try { + current_md5 = QuickDigest5::fileToHash(pyname); + } catch (const std::runtime_error& e) { + if (myid==0) + std::cerr << "Cylinder::CheckDtype error: " + << e.what() << std::endl; + } + + // Check that the md5sums match for the current Python + // module source files and the loaded Python module used to + // create the cache. If they do not match, force cache + // recomputation to ensure consistency with the current + // Python module. + // + if (current_md5 != pyinfo[1]) { + if (myid==0) { + std::cout << "---- Cylinder::checkDtype: Python module for disk density has changed since cache creation." << std::endl + << "---- Current module: <" << pyname << ">, md5sum: " << current_md5 << std::endl + << "---- Loaded module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl + << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; + } + cache_status = false; + } } - cache_status = false; + // End: have Python disk type, check md5 hashes } + // End: Python disk type check } - // End: have Python disk type, check md5 hashes + // End: DiskType attribute check + + // Could add deprojection checks here in the future } - // End: Python disk type check - - // Could add deprojection checks here in the future + // End: DiskType attribute check } - // End: DiskType attribute check + // End: EmpModel attribute check return cache_status; } -void Cylinder::saveDtype() +void Cylinder::saveMetaData() { // Only one process must write the dtype metadata // diff --git a/utils/ICs/check_coefs.cc b/utils/ICs/check_coefs.cc index b8556090b..9d1a9f5f2 100644 --- a/utils/ICs/check_coefs.cc +++ b/utils/ICs/check_coefs.cc @@ -471,35 +471,36 @@ main(int ac, char **av) } } + // Set EmpCylSL mtype. This is the spherical function used to + // generate the EOF basis. If "deproject" is set, this will be + // overriden in EmpCylSL. + // + // Convert mtype string to lower case // std::transform(mtype.begin(), mtype.end(), mtype.begin(), [](unsigned char c){ return std::tolower(c); }); - - // Set EmpCylSL mtype - // - EmpCylSL::mtype = EmpCylSL::Exponential; - if (vm.count("mtype")) { - if (mtype.compare("exponential")==0) - EmpCylSL::mtype = EmpCylSL::Exponential; - else if (mtype.compare("expsphere")==0) - EmpCylSL::mtype = EmpCylSL::ExpSphere; - else if (mtype.compare("gaussian")==0) - EmpCylSL::mtype = EmpCylSL::Gaussian; - else if (mtype.compare("plummer")==0) - EmpCylSL::mtype = EmpCylSL::Plummer; - else if (mtype.compare("power")==0) { - EmpCylSL::mtype = EmpCylSL::Power; - EmpCylSL::PPOW = ppower; - } else { - if (myid==0) std::cout << "No EmpCylSL EmpModel named <" - << mtype << ">, valid types are: " - << "Exponential, ExpSphere, Gaussian, Plummer, Power" << std::endl; - MPI_Finalize(); - return -1; + + // Set EmpCylSL mtype. This is the spherical function used to + // generate the EOF basis. If "deproject" is set, this will be + // overriden in EmpCylSL. + // + auto itm = EmpCylSL::EmpModelMap.find(mtype); + + if (itm == EmpCylSL::EmpModelMap.end()) { + if (myid==0) { + std::cout << "No EmpCylSL EmpModel named <" + << mtype << ">, valid types are: " + << "Exponential, ExpSphere, Gaussian, Plummer, Power, Deproject " + << "(not case sensitive)" << std::endl; } + + MPI_Finalize(); + return 0; } + EmpCylSL::mtype = itm->second; + // Set DiskType // std::transform(disktype.begin(), disktype.end(), disktype.begin(), diff --git a/utils/ICs/check_coefs2.cc b/utils/ICs/check_coefs2.cc index 923a0edb4..225eabe6c 100644 --- a/utils/ICs/check_coefs2.cc +++ b/utils/ICs/check_coefs2.cc @@ -479,30 +479,26 @@ main(int ac, char **av) std::transform(mtype.begin(), mtype.end(), mtype.begin(), [](unsigned char c){ return std::tolower(c); }); - // Set EmpCylSL mtype + // Set EmpCylSL mtype. This is the spherical function used to + // generate the EOF basis. If "deproject" is set, this will be + // overriden in EmpCylSL. // - EmpCylSL::mtype = EmpCylSL::Exponential; - if (vm.count("mtype")) { - if (mtype.compare("exponential")==0) - EmpCylSL::mtype = EmpCylSL::Exponential; - else if (mtype.compare("expsphere")==0) - EmpCylSL::mtype = EmpCylSL::ExpSphere; - else if (mtype.compare("gaussian")==0) - EmpCylSL::mtype = EmpCylSL::Gaussian; - else if (mtype.compare("plummer")==0) { - EmpCylSL::mtype = EmpCylSL::Plummer; - } else if (mtype.compare("power")==0) { - EmpCylSL::mtype = EmpCylSL::Power; - EmpCylSL::PPOW = ppower; - } else { - if (myid==0) std::cout << "No EmpCylSL EmpModel named <" - << mtype << ">, valid types are: " - << "Exponential, ExpSphere, Gaussian, Plummer, Power" << std::endl; - MPI_Finalize(); - return -1; + auto itm = EmpCylSL::EmpModelMap.find(mtype); + + if (itm == EmpCylSL::EmpModelMap.end()) { + if (myid==0) { + std::cout << "No EmpCylSL EmpModel named <" + << mtype << ">, valid types are: " + << "Exponential, ExpSphere, Gaussian, Plummer, Power, Deproject " + << "(not case sensitive)" << std::endl; } + + MPI_Finalize(); + return 0; } + EmpCylSL::mtype = itm->second; + // Set DiskType // std::transform(disktype.begin(), disktype.end(), disktype.begin(), diff --git a/utils/ICs/initial.cc b/utils/ICs/initial.cc index 5f2201aa6..c94929b46 100644 --- a/utils/ICs/initial.cc +++ b/utils/ICs/initial.cc @@ -677,41 +677,41 @@ main(int ac, char **av) SphericalModelTable::linear = 1; } - // Convert mtype string to lower case + // Convert mtype string to lower case + // std::transform(mtype.begin(), mtype.end(), mtype.begin(), [](unsigned char c){ return std::tolower(c); }); - // Convert gentype string to lower case + // Convert gentype string to lower case + // std::transform(gentype.begin(), gentype.end(), gentype.begin(), [](unsigned char c){ return std::tolower(c); }); + // Convert mtype string to lower case + // + std::transform(mtype.begin(), mtype.end(), mtype.begin(), + [](unsigned char c){ return std::tolower(c); }); + // Set EmpCylSL mtype. This is the spherical function used to // generate the EOF basis. If "deproject" is set, this will be // overriden in EmpCylSL. // - EmpCylSL::mtype = EmpCylSL::Exponential; - if (vm.count("mtype")) { - if (mtype.compare("exponential")==0) - EmpCylSL::mtype = EmpCylSL::Exponential; - else if (mtype.compare("expsphere")==0) - EmpCylSL::mtype = EmpCylSL::ExpSphere; - else if (mtype.compare("gaussian")==0) - EmpCylSL::mtype = EmpCylSL::Gaussian; - else if (mtype.compare("plummer")==0) - EmpCylSL::mtype = EmpCylSL::Plummer; - else if (mtype.compare("power")==0) { - EmpCylSL::mtype = EmpCylSL::Power; - EmpCylSL::PPOW = PPower; - } else { - if (myid==0) std::cout << "No EmpCylSL EmpModel named <" - << mtype << ">, valid types are: " - << "Exponential, Gaussian, Plummer, Power " - << "(not case sensitive)" << std::endl; - MPI_Finalize(); - return 0; + auto itm = EmpCylSL::EmpModelMap.find(mtype); + + if (itm == EmpCylSL::EmpModelMap.end()) { + if (myid==0) { + std::cout << "No EmpCylSL EmpModel named <" + << mtype << ">, valid types are: " + << "Exponential, ExpSphere, Gaussian, Plummer, Power, Deproject " + << "(not case sensitive)" << std::endl; } + + MPI_Finalize(); + return 0; } + EmpCylSL::mtype = itm->second; + // Set DiskType. This is the functional form for the disk used to // condition the basis. // From b23224b3dbcb31eea314a4469c1a96c757df242a Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Mar 2026 22:20:05 -0400 Subject: [PATCH 72/78] Missing case transformation for dmodel --- expui/BiorthBasis.cc | 31 +++++++++++++++---------------- 1 file changed, 15 insertions(+), 16 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index ff35f86a4..7cc70a1e0 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1578,17 +1578,27 @@ namespace BasisClasses EmpCylSL::logarithmic = logarithmic; EmpCylSL::VFLAG = vflag; + // Convert dmodel string to lower case + // + std::transform(dmodel.begin(), dmodel.end(), dmodel.begin(), + [](unsigned char c){ return std::tolower(c); }); + + // Convert mtype string to lower case + // + std::transform(mtype.begin(), mtype.end(), mtype.begin(), + [](unsigned char c){ return std::tolower(c); }); + + // Convert dtype string to lower case + // + std::transform(dtype.begin(), dtype.end(), dtype.begin(), + [](unsigned char c){ return std::tolower(c); }); + // Set EmpCylSL mtype. This is the spherical function used to // generate the EOF basis. If "deproject" is set, this will be // auto-set to the same as the deprojection model. Otherwise, // this can be set independently to allow for different spherical // functions for the EOF basis - // Convert mtype string to lower case - // - std::transform(mtype.begin(), mtype.end(), mtype.begin(), - [](unsigned char c){ return std::tolower(c); }); - // Set EmpCylSL mtype. This is the spherical function used to // generate the EOF basis. // @@ -1665,12 +1675,6 @@ namespace BasisClasses {DiskType::python, "python"} }; - // Convert dtype string to lower case - // - std::transform(dtype.begin(), dtype.end(), dtype.begin(), - [](unsigned char c){ return std::tolower(c); }); - - // Check for map entry, will throw if the key is not in the map. DTYPE = dtlookup.at(dtype); @@ -1936,11 +1940,6 @@ namespace BasisClasses // EmpCylSL::AxiDiskPtr model; - // Convert dmodel string to lower case - // - std::transform(dmodel.begin(), dmodel.end(), dmodel.begin(), - [](unsigned char c){ return std::tolower(c); }); - // Map legacy/short model names to canonical keys expected by dplookup // if (dmodel == "exp") { From fe2e667d9bf9f6e11f9163cb261b19758796c70d Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Mar 2026 22:28:36 -0400 Subject: [PATCH 73/78] Additional comments only [no ci] --- expui/BiorthBasis.cc | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 7cc70a1e0..def3d1e1d 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1578,17 +1578,20 @@ namespace BasisClasses EmpCylSL::logarithmic = logarithmic; EmpCylSL::VFLAG = vflag; - // Convert dmodel string to lower case + // Convert dmodel string to lower case (deprojection model for EOF + // basis construction) // std::transform(dmodel.begin(), dmodel.end(), dmodel.begin(), [](unsigned char c){ return std::tolower(c); }); - // Convert mtype string to lower case + // Convert mtype string to lower case (EmpCylSL spherical function + // for EOF basis construction) // std::transform(mtype.begin(), mtype.end(), mtype.begin(), [](unsigned char c){ return std::tolower(c); }); - // Convert dtype string to lower case + // Convert dtype string to lower case (disk density function for + // EOF conditioning) // std::transform(dtype.begin(), dtype.end(), dtype.begin(), [](unsigned char c){ return std::tolower(c); }); From 567f0050c5aea77b238f74fe0a83995c96f64590 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Tue, 24 Mar 2026 22:50:45 -0400 Subject: [PATCH 74/78] Minor sanity-check logic tweak --- expui/BiorthBasis.cc | 158 +++++++++++++++++++++++-------------------- 1 file changed, 85 insertions(+), 73 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index def3d1e1d..cc71539c1 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1692,18 +1692,26 @@ namespace BasisClasses HighFive::File file(cachename, HighFive::File::ReadOnly); if (!file.hasAttribute("EmpModel")) { - // If the EmpModel attribute is missing, we have an old cache type + // If the EmpModel attribute is missing, we have an old + // cache type. We will allow this for now to avoid breaking + // old caches, but we will trigger a cache build in the future + // if (myid==0) { std::cout << "---- Cylindrical: EmpModel attribute not found in cache file <" << cachename << ">. " << std::endl; std::cout << "---- Cylindrical: this may indicate an old cache file created before EmpModel metadata was added. " << std::endl; std::cout << "---- Cylindrical: we will continue...but consider remaking the cache to avoid confusion." << std::endl; } } else { + // Open existing EmpModel attribute and the the mtype string + // auto read_attr = file.getAttribute("EmpModel"); std::string empmodel; read_attr.read(empmodel); + // Check against the current mtype string. If they do not + // match, force cache recomputation to ensure consistency + // with the current mtype. if (empmodel != mtype) { if (myid==0) { std::cout << "---- Cylindrical: EmpModel for cache file <" @@ -1719,84 +1727,89 @@ namespace BasisClasses } } - // Sanity: check that the DiskType attribute exists - // - if (!file.hasAttribute("DiskType")) { - if (myid==0) { - std::cout << "---- Cylindrical: DiskType attribute not found in cache file <" << cachename << ">. " << std::endl - << "---- This is a logic error since the DiskType attribute is required for cache creation" << std::endl - << "---- We will trigger cache recomputation to be safe." << std::endl; - } - } else { - // Open existing DiskType attribute - // - auto read_attr = file.getAttribute("DiskType"); - - std::string loaded_dtype; - read_attr.read(loaded_dtype); - - DiskType disktype = dtlookup.at(loaded_dtype); + if (cache_status == 1) { - if (disktype != DTYPE) { + // Sanity: check that the DiskType attribute exists + // + if (!file.hasAttribute("DiskType")) { if (myid==0) { - std::cout << "---- Cylindrical: DiskType for cache file <" - << cachename << "> is <" - << loaded_dtype << ">," << std::endl - << "which does not match the requested DiskType <" - << dtype << ">" << std::endl - << "---- Cylindrical: forcing cache recomputation" - << std::endl; + std::cout << "---- Cylindrical: DiskType attribute not found in cache file <" << cachename << ">. " << std::endl + << "---- This is a logic error since the DiskType attribute is required for cache creation" << std::endl + << "---- We will trigger cache recomputation to be safe." << std::endl; } // Force cache recomputation - cache_status = 0; - } - else if (disktype == DiskType::python) { - - // If the DiskType is python, we also need to check that - // the Python module used to create the cache matches the - // current Python module to ensure consistency. This tag - // should be present in the cache if it was created with a - // recent version of the code, but if it is missing we - // will trigger cache recomputation to ensure consistency - // with the current Python module. - if (file.hasAttribute("pythonDiskType") == false) { - if (myid==0) { - std::cout << "---- Cylindrical: pythonDiskType attribute not found in cache file <" << cachename << ">. " << std::endl; - std::cout << "---- Cylindrical: this may be a logic error, trigger recomputation." << std::endl; - } - cache_status = 0; - } else { - // Get the pyname attribute - auto read_attr = file.getAttribute("pythonDiskType"); + cache_status = 0; + } else { + // Open existing DiskType attribute + // + auto read_attr = file.getAttribute("DiskType"); - std::vector pyinfo; - read_attr.read(pyinfo); - - std::string current_md5; + std::string loaded_dtype; + read_attr.read(loaded_dtype); - // Get the md5sum for requested Python module - try { - current_md5 = QuickDigest5::fileToHash(pyname + ".py"); - } catch (const std::runtime_error& e) { - if (myid==0) - std::cerr << "BiorthBasis::Cylindrical error: " - << e.what() << ", error compuing pyname md5sum" - << std::endl; + DiskType disktype = dtlookup.at(loaded_dtype); + + if (disktype != DTYPE) { + if (myid==0) { + std::cout << "---- Cylindrical: DiskType for cache file <" + << cachename << "> is <" + << loaded_dtype << ">," << std::endl + << "which does not match the requested DiskType <" + << dtype << ">" << std::endl + << "---- Cylindrical: forcing cache recomputation" + << std::endl; } - - // Check that the md5sums match for the current Python - // module and the Python module used to create the - // cache. If they do not match, force cache - // recomputation to ensure consistency with the current - // Python module. - if (current_md5 != pyinfo[1]) { + // Force cache recomputation + cache_status = 0; + } + else if (disktype == DiskType::python) { + + // If the DiskType is python, we also need to check that + // the Python module used to create the cache matches the + // current Python module to ensure consistency. This tag + // should be present in the cache if it was created with a + // recent version of the code, but if it is missing we + // will trigger cache recomputation to ensure consistency + // with the current Python module. + if (file.hasAttribute("pythonDiskType") == false) { if (myid==0) { - std::cout << "---- Cylindrical: Python module for disk density has changed since cache creation." << std::endl - << "---- Current module: <" << pyname << ">, md5sum: " << current_md5 << std::endl - << "---- Cached module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl - << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; + std::cout << "---- Cylindrical: pythonDiskType attribute not found in cache file <" << cachename << ">. " << std::endl; + std::cout << "---- Cylindrical: this may be a logic error, trigger recomputation." << std::endl; } cache_status = 0; + } else { + // Get the pyname attribute + auto read_attr = file.getAttribute("pythonDiskType"); + + std::vector pyinfo; + read_attr.read(pyinfo); + + std::string current_md5; + + // Get the md5sum for requested Python module + try { + current_md5 = QuickDigest5::fileToHash(pyname + ".py"); + } catch (const std::runtime_error& e) { + if (myid==0) + std::cerr << "BiorthBasis::Cylindrical error: " + << e.what() << ", error compuing pyname md5sum" + << std::endl; + } + + // Check that the md5sums match for the current Python + // module and the Python module used to create the + // cache. If they do not match, force cache + // recomputation to ensure consistency with the current + // Python module. + if (current_md5 != pyinfo[1]) { + if (myid==0) { + std::cout << "---- Cylindrical: Python module for disk density has changed since cache creation." << std::endl + << "---- Current module: <" << pyname << ">, md5sum: " << current_md5 << std::endl + << "---- Cached module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl + << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; + } + cache_status = 0; + } } } } @@ -1808,7 +1821,7 @@ namespace BasisClasses // Get the dmodel attribute // - read_attr = file.getAttribute("ProjType"); + auto read_attr = file.getAttribute("ProjType"); std::string loaded_dmodel; read_attr.read(loaded_dmodel); @@ -1838,7 +1851,8 @@ namespace BasisClasses cache_status = 0; } else { - read_attr = file.getAttribute("pythonProjType"); + // Get the pyproj attribute and md5 hash from the cache + auto read_attr = file.getAttribute("pythonProjType"); std::vector pyinfo; read_attr.read(pyinfo); @@ -1881,7 +1895,6 @@ namespace BasisClasses // Remake cylindrical basis // - if (Ignore and myid==0) { std::cout << "---- BasisFactory: We have deprecated the 'ignore' parameter for the" << std::endl << "---- Cylindrical class. If the cache file exists but does" << std::endl @@ -1891,7 +1904,6 @@ namespace BasisClasses << "---- remove 'ignore' from your YAML configuration." << std::endl; } - // Set DiskType. This is the functional form for the disk used to // condition the basis. // From 5fc11984a0147e2b7d7a348eb39d407b2d6721ba Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 25 Mar 2026 09:28:00 -0400 Subject: [PATCH 75/78] Add defaults for mtype and dtype to Cylinder --- expui/BiorthBasis.cc | 2 +- src/Cylinder.cc | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index cc71539c1..82c4d6eaa 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1416,7 +1416,7 @@ namespace BasisClasses EVEN_M = false; cmapR = 1; cmapZ = 1; - mtype = "Exponential"; + mtype = "exponential"; dtype = "exponential"; vflag = 0; diff --git a/src/Cylinder.cc b/src/Cylinder.cc index 30f00604b..f80618cfe 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -127,6 +127,9 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : ncylodd = 9; ncylrecomp = -1; + mtype = "exponential"; + dtype = "exponential"; + // For disk basis construction with doubleexpon // aratio = 1.0; // Radial scale length ratio From 0c49ed056c788e027c00903186c11da109c1b894 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Wed, 25 Mar 2026 10:43:25 -0400 Subject: [PATCH 76/78] Deduplicate cache tags, oops --- expui/BiorthBasis.cc | 78 +++++++------------- src/Cylinder.cc | 165 ++++++++++++++++++------------------------- 2 files changed, 94 insertions(+), 149 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 82c4d6eaa..7191ef5df 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1691,54 +1691,27 @@ namespace BasisClasses // HighFive::File file(cachename, HighFive::File::ReadOnly); - if (!file.hasAttribute("EmpModel")) { - // If the EmpModel attribute is missing, we have an old - // cache type. We will allow this for now to avoid breaking - // old caches, but we will trigger a cache build in the future - // - if (myid==0) { - std::cout << "---- Cylindrical: EmpModel attribute not found in cache file <" << cachename << ">. " << std::endl; - std::cout << "---- Cylindrical: this may indicate an old cache file created before EmpModel metadata was added. " << std::endl; - std::cout << "---- Cylindrical: we will continue...but consider remaking the cache to avoid confusion." << std::endl; - } - } else { - // Open existing EmpModel attribute and the the mtype string - // - auto read_attr = file.getAttribute("EmpModel"); - - std::string empmodel; - read_attr.read(empmodel); - - // Check against the current mtype string. If they do not - // match, force cache recomputation to ensure consistency - // with the current mtype. - if (empmodel != mtype) { - if (myid==0) { - std::cout << "---- Cylindrical: EmpModel for cache file <" - << cachename << "> is <" - << empmodel << ">," << std::endl - << "which does not match the requested EmpModel <" - << mtype << ">" << std::endl - << "---- Cylindrical: forcing cache recomputation" - << std::endl; - } - // Force cache recomputation - cache_status = 0; - } - } - - if (cache_status == 1) { + // mtype is already cached as "model" attribute in the HDF5 + // cache file, so we do not need to write it here. We just + // need to write the DiskType and Python metadata (if + // applicable). + // - // Sanity: check that the DiskType attribute exists - // - if (!file.hasAttribute("DiskType")) { - if (myid==0) { + // Check that the DiskType for the cache matches the requested + // DiskType If the DiskType attribute is missing from the + // cache, this indicates thet the cache was created with an + // older versions of the code that did not include the + // DiskType attribute in the cache. We should trigger a cache + // recomputation in this case to be safe, but will let this + // pass for backward compatibility with older caches. + // + if (!file.hasAttribute("DiskType")) { + if (myid==0) { std::cout << "---- Cylindrical: DiskType attribute not found in cache file <" << cachename << ">. " << std::endl - << "---- This is a logic error since the DiskType attribute is required for cache creation" << std::endl - << "---- We will trigger cache recomputation to be safe." << std::endl; - } - // Force cache recomputation - cache_status = 0; + << "---- This indicates that the cache was created with an older version of the" << std::endl + << "---- code that did not include the DiskType attribute in the cache. We" << std::endl + << "---- suggest deleting the old cache file to prevent this warning and rebuild" << std::endl + << "---- the cache, but we will continue..." << std::endl; } else { // Open existing DiskType attribute // @@ -1880,7 +1853,9 @@ namespace BasisClasses } } } + // End: DeprojType and Python module consistency checks with cache } + // End: DiskType and Python module consistency checks with cache } catch (const HighFive::Exception& err) { if (myid==0) { @@ -2033,12 +2008,11 @@ namespace BasisClasses // HighFive::File file(cachename, HighFive::File::ReadWrite); - file.createAttribute("EmpModel", mtype); - - std::cout << "---- Cylindrical: writing EmpModel <" - << EmpCylSL::EmpModelLabs.at(EmpCylSL::mtype) - << "> to cache file <" << cachename << ">" << std::endl; - + // mtype is already cached as "model" attribute in the HDF5 + // cache file, so we do not need to write it here. We just + // need to write the DiskType and Python metadata (if + // applicable). + // file.createAttribute("DiskType", dtype); std::cout << "---- Cylindrical: writing DiskType <" << dtype diff --git a/src/Cylinder.cc b/src/Cylinder.cc index f80618cfe..b50ebdf12 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -2069,125 +2069,92 @@ bool Cylinder::checkMetaData() // HighFive::File file(cachename, HighFive::File::ReadOnly); - // Sanity: check that the EmpModel attribute exists + // Sanity: check that the DiskType attribute exists // - if (!file.hasAttribute("EmpModel")) { + if (!file.hasAttribute("DiskType")) { if (myid==0) { - std::cout << "---- Cylinder:checkDtype: EmpModel attribute not found in cache file <" << cachename << ">" << std::endl + std::cout << "---- Cylinder:checkMetaData: DiskType attribute not found in cache file <" << cachename << ">" << std::endl << "---- This may indicate an old cache file created before EmpModel metadata was added" << std::endl << "---- We will continue...but consider remaking the cache to avoid confusion" << std::endl; } - } else { - // Check that the EmpModel attribute matches the current empirical model + } + else { + // Open existing DiskType attribute // - auto read_attr = file.getAttribute("EmpModel"); - - std::string loaded_mtype; - read_attr.read(loaded_mtype); - - if (loaded_mtype != mtype) { + auto read_attr = file.getAttribute("DiskType"); + + std::string loaded_dtype; + read_attr.read(loaded_dtype); + + // Map the loaded dtype string to a DiskType enum value + // + DiskType disktype = dtlookup.at(loaded_dtype); + + if (disktype != DTYPE) { if (myid==0) { - std::cout << "---- Cylinder::checkDtype: EmpModel for cache file <" - << cachename << "> is <" << loaded_mtype << ">," - << " which does not match the requested empirical model <" - << mtype << ">" << std::endl - << "---- Cylindrical: forcing cache recomputation" + std::cout << "---- Cylinder::checkMetaData: DiskType for cache file <" + << cachename << "> is <" + << loaded_dtype << ">," << std::endl + << "which does not match the requested DiskType <" + << dtype << ">" << std::endl + << "---- Cylinder::checkMetaData: forcing cache recomputation" << std::endl; } // Force cache recomputation cache_status = false; } - else { - - if (!file.hasAttribute("DiskType")) { + else if (disktype == DiskType::python) { + // Sanity check: if DiskType is python, then the + // pythonDiskType attribute must exist + // + if (!file.hasAttribute("pythonDiskType")) { if (myid==0) { - std::cout << "---- Cylinder:checkDtype: DiskType attribute not found in cache file <" << cachename << ">" << std::endl - << "---- Cylinder:checkDtype: this may indicate a logic error. Forcing cache recomputation." << std::endl; + std::cout << "---- Cylinder::checkMetaData: pythonDiskType attribute not found in cache file <" << cachename << ">. " << std::endl; + std::cout << "---- Cylinder::checkMetaData: this may indicate a logic error. Forcing cache recomputation." << std::endl; } // Force cache recomputation cache_status = false; - } - else { - // Open existing DiskType attribute - // - auto read_attr = file.getAttribute("DiskType"); - - std::string loaded_dtype; - read_attr.read(loaded_dtype); - - // Map the loaded dtype string to a DiskType enum value - // - DiskType disktype = dtlookup.at(loaded_dtype); - - if (disktype != DTYPE) { - if (myid==0) { - std::cout << "---- Cylinder::checkDtype: DiskType for cache file <" - << cachename << "> is <" - << loaded_dtype << ">," << std::endl - << "which does not match the requested DiskType <" - << dtype << ">" << std::endl - << "---- Cylindrical: forcing cache recomputation" - << std::endl; - } - // Force cache recomputation - cache_status = false; - } - else if (disktype == DiskType::python) { - // Sanity check: if DiskType is python, then the - // pythonDiskType attribute must exist - // - if (!file.hasAttribute("pythonDiskType")) { - if (myid==0) { - std::cout << "---- Cylinder::checkDtype: pythonDiskType attribute not found in cache file <" << cachename << ">. " << std::endl; - std::cout << "---- Cylindier:checkDtype: this may indicate a logic error. Forcing cache recomputation." << std::endl; - } - // Force cache recomputation - cache_status = false; - } else { - auto read_attr = file.getAttribute("pythonDiskType"); + } else { + auto read_attr = file.getAttribute("pythonDiskType"); - // Get the pyname attribute - std::vector pyinfo; - read_attr.read(pyinfo); + // Get the pyname attribute + std::vector pyinfo; + read_attr.read(pyinfo); - std::string current_md5; + std::string current_md5; - // Get the md5sum for requested Python module source file - try { - current_md5 = QuickDigest5::fileToHash(pyname); - } catch (const std::runtime_error& e) { - if (myid==0) - std::cerr << "Cylinder::CheckDtype error: " - << e.what() << std::endl; - } + // Get the md5sum for requested Python module source file + try { + current_md5 = QuickDigest5::fileToHash(pyname); + } catch (const std::runtime_error& e) { + if (myid==0) + std::cerr << "Cylinder::checkMetaData error: " + << e.what() << std::endl; + } - // Check that the md5sums match for the current Python - // module source files and the loaded Python module used to - // create the cache. If they do not match, force cache - // recomputation to ensure consistency with the current - // Python module. - // - if (current_md5 != pyinfo[1]) { - if (myid==0) { - std::cout << "---- Cylinder::checkDtype: Python module for disk density has changed since cache creation." << std::endl - << "---- Current module: <" << pyname << ">, md5sum: " << current_md5 << std::endl - << "---- Loaded module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl - << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; - } - cache_status = false; - } + // Check that the md5sums match for the current Python + // module source files and the loaded Python module used to + // create the cache. If they do not match, force cache + // recomputation to ensure consistency with the current + // Python module. + // + if (current_md5 != pyinfo[1]) { + if (myid==0) { + std::cout << "---- Cylinder::checkMetaData: Python module for disk density has changed since cache creation." << std::endl + << "---- Current module: <" << pyname << ">, md5sum: " << current_md5 << std::endl + << "---- Loaded module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl + << "---- Cylinder:checkMetaData: forcing cache recomputation to ensure consistency" << std::endl; } - // End: have Python disk type, check md5 hashes + cache_status = false; } - // End: Python disk type check } - // End: DiskType attribute check - - // Could add deprojection checks here in the future + // End: have Python disk type, check md5 hashes } - // End: DiskType attribute check + // End: Python disk type check + + // Could add deprojection checks here in the future } - // End: EmpModel attribute check + // End: DiskType attribute check return cache_status; } @@ -2199,7 +2166,11 @@ void Cylinder::saveMetaData() if (myid) return; std::string dtype = dtstring.at(DTYPE); - + + // mtype is already cached as "model" attribute in the HDF5 cache + // file, so we do not need to write it here. We just need to write + // the DiskType and Python metadata (if applicable). + // try { // Reopen the cache file for writing the Python metadata // @@ -2223,7 +2194,7 @@ void Cylinder::saveMetaData() file.createAttribute("pythonDiskType", pyinfo); } catch (const std::runtime_error& e) { - std::cerr << "Cylinder::saveDtype error: " << e.what() << std::endl; + std::cerr << "Cylinder::saveMetaData error: " << e.what() << std::endl; std::cerr << "Can not write the md5 hash to HDF5" << std::endl; } } @@ -2232,7 +2203,7 @@ void Cylinder::saveMetaData() } catch (const HighFive::Exception& err) { std::cerr << err.what() << std::endl; - std::cerr << "Cylinder::saveDtype: error writing metadata to cache file <" + std::cerr << "Cylinder::saveMetaData: error writing metadata to cache file <" << cachename << std::endl; } } From 2ef3aba295e4e4d8cbae16a713133b6f1cbe2eca Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 26 Mar 2026 11:11:34 -0400 Subject: [PATCH 77/78] Added automatic output of valid EmpModel types --- expui/BiorthBasis.cc | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/expui/BiorthBasis.cc b/expui/BiorthBasis.cc index 7191ef5df..3592a6c67 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1605,19 +1605,19 @@ namespace BasisClasses // Set EmpCylSL mtype. This is the spherical function used to // generate the EOF basis. // - auto itm = EmpCylSL::EmpModelMap.find(mtype); - - if (itm == EmpCylSL::EmpModelMap.end()) { + try { + auto itm = EmpCylSL::EmpModelMap.find(mtype); + EmpCylSL::mtype = itm->second; + } + catch (const std::out_of_range& err) { if (myid==0) { - std::cout << "No EmpCylSL EmpModel named <" - << mtype << ">, valid types are: " - << "Exponential, ExpSphere, Gaussian, Plummer, Power, Deproject " - << "(not case sensitive)" << std::endl; + std::cout << "Cylindrical::initialize error parsing 'mtype' parameter in YAML config" << std::endl; + std::cout << "Valid options are: "; + for (auto p : EmpCylSL::EmpModelLabs) std::cout << p.second << " "; + std::cout << "(not case sensitive)" << std::endl; } - throw std::runtime_error("Cylindrical:initialize: EmpCylSL bad parameter"); + throw std::runtime_error("Cylindrical::initialize: invalid 'mtype' parameter in YAML config"); } - - EmpCylSL::mtype = itm->second; // Check for non-null cache file name. This must be specified // to prevent recomputation and unexpected behavior. From 73a47f657640f3835578298b8a864c5ff62b29a1 Mon Sep 17 00:00:00 2001 From: "Martin D. Weinberg" Date: Thu, 26 Mar 2026 11:11:52 -0400 Subject: [PATCH 78/78] Added deprojection support to Cylinder --- src/Cylinder.H | 43 +++++++- src/Cylinder.cc | 253 +++++++++++++++++++++++++++++++++++++++++------- 2 files changed, 258 insertions(+), 38 deletions(-) diff --git a/src/Cylinder.H b/src/Cylinder.H index 653b673e1..6103321a5 100644 --- a/src/Cylinder.H +++ b/src/Cylinder.H @@ -46,6 +46,18 @@ class MixtureBasis; @param dweight is the relative weight of the second component in the double-exponential disk basis construction (default: 1.0) + @param dmodel is the disk model type for deprojection and EOF conditioning. See Cylinder::DiskType for options. + + @param mtype is the spherical model type for conditioning the basis. See EmpCylSL::EmpModel for options. + + @param dtype is the density model type for conditioning the basis. See Cylinder::DiskType for options. + + @param rwidth is the radial width of the error function truncation for disk basis construction (default: 0.0, which implies no truncation) + + @param rtrunc is the radial truncation radius for disk basis construction (default: 0.1) + + @param rfactor is the radial scale factor for disk basis construction (default: 1.0, which implies no scaling) + @param lmaxfid is the maximum spherical harmonic index for EOF construction @param nmaxfid is the maximum radial index for EOF construction @@ -121,7 +133,7 @@ class MixtureBasis; @param cmapz is the vertical coordinate mapping type - @param mtype is the deprojected model type for conditioning the + @param mtype is the spherical model type for conditioning the basis. See ExpCylSL::mtype for options. @param ppower is the power-law index for the Power spherical model @@ -170,6 +182,7 @@ private: // double rcylmin, rcylmax, zmax, acyl, bias; double aratio, hratio, dweight; // Disk basis construction parameters for double-exponential disk + double rwidth, rtrunc, rfactor; // Disk basis construction parameters for numerical deprojection double Mfac, HERNA; // Disk bulge parameters: mass fraction and Hernquist scale length int nmaxfid, lmaxfid, mmax, mlim, nint; @@ -355,7 +368,7 @@ protected: //! Power-law index for deprojected spherical model for basis conditioning double ppow = 2.0; - //! Deprojection model for basis conditioning + //! Spherical model for basis conditioning std::string mtype = "Exponential"; //! Write basis-specific parameters to HDF5 covariance file @@ -365,6 +378,9 @@ protected: enum class DiskType { constant, gaussian, mn, exponential, doubleexpon, diskbulge, python }; + //! The string name of the deprojection model + std::string dmodel; + //! Map human strings to disk target typenums const std::map dtlookup = { {"constant", DiskType::constant}, @@ -393,6 +409,29 @@ protected: //! The string name of the disk density function type for basis conditioning std::string dtype; + //@{ + //! DeprojType support + + //! Disk models used for deprojection + enum class DeprojType + { mn, toomre, python, exponential}; + + //! Current model + DeprojType PTYPE; + + //! Python module name for disk density function if using Python deprojection + std::string pyproj; + + //! Look up by string + const std::map dplookup = + { {"mn", DeprojType::mn}, + {"toomre", DeprojType::toomre}, + {"python", DeprojType::python}, + {"exponential", DeprojType::exponential} + }; + //@} + + //! Dtype cache check for consistency between cache and current parameters bool checkMetaData(); diff --git a/src/Cylinder.cc b/src/Cylinder.cc index b50ebdf12..50aa2dd8f 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -15,6 +15,7 @@ #include "exputils.H" // utility functions #include "NVTX.H" // for NVTX profiling of CUDA code #include "quickdigest5.hpp" // for md5 hashing of Python modules +#include "DiskModels.H" // //@{ //! These are for testing exclusively (should be set false for production) @@ -61,6 +62,9 @@ Cylinder::valid_keys = { "pnum", "tnum", "ashift", + "rwidth", + "rtrunc", + "rfactor", "expcond", "precond", "logr", @@ -129,6 +133,7 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : mtype = "exponential"; dtype = "exponential"; + dmodel = "EXP"; // For disk basis construction with doubleexpon // @@ -144,7 +149,9 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : pnum = 1; tnum = 80; ashift = 0.0; - + rwidth = 0.0; // Width of error function truncation (ignored if zero) + rfactor = 1.0; // Radial scale factor for numerical basis construction + rtrunc = 0.1; // Radial truncation for numerical basis construction vflag = 0; eof = 1; npca = std::numeric_limits::max(); @@ -214,10 +221,10 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : if (itm == EmpCylSL::EmpModelMap.end()) { if (myid==0) { - std::cout << "No EmpCylSL EmpModel named <" - << mtype << ">, valid types are: " - << "Exponential, ExpSphere, Gaussian, Plummer, Power, Deproject " - << "(not case sensitive)" << std::endl; + std::cout << "No EmpCylSL EmpModel named <" + << mtype << ">, valid types are: "; + for (auto p : EmpCylSL::EmpModelLabs) std::cout << p.second << " "; + std::cout << "(not case sensitive)" << std::endl; } throw std::runtime_error("Cylindrical:initialize: EmpCylSL bad parameter"); } @@ -234,8 +241,24 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : std::transform(dtype.begin(), dtype.end(), dtype.begin(), [](unsigned char c){ return std::tolower(c); }); + // Convert dmodel to lower case + // + std::transform(dmodel.begin(), dmodel.end(), dmodel.begin(), + [](unsigned char c){ return std::tolower(c); }); + // Check for map entry, will throw if the key is not in the map. - DTYPE = dtlookup.at(dtype); + try { + DTYPE = dtlookup.at(dtype); + } + catch (const std::out_of_range& err) { + if (myid==0) { + std::cout << "DiskType error in configuraton file" << std::endl; + std::cout << "Valid options are: "; + for (auto v : dtlookup) std::cout << v.first << " "; + std::cout << std::endl; + } + throw std::runtime_error("Cylinder: invalid DiskType"); + } // Set azimuthal harmonic order restriction? // @@ -329,37 +352,26 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : << std::endl; - // Set DiskType. This is the functional form for the disk used to - // condition the basis. + // Warning about DiskType and pyEXP assumptions for backward compatibility with + // previous versions of pyEXP. // - try { - if (myid==0) { // Report DiskType - std::cout << "---- DiskType is <" << dtype << ">" << std::endl; - - if (not sech2) { - switch (DTYPE) { - case DiskType::doubleexpon: - case DiskType::exponential: - case DiskType::diskbulge: - std::cout << "---- pyEXP assumes sech^2(z/(2h)) by default in v7.9.0 and later" << std::endl - << "---- Use the 'sech2: true' in your YAML config to use sech^2(z/(2h))" << std::endl - << "---- This warning will be removed in v7.10.0." << std::endl; - break; - default: - break; - } + if (myid==0) { + std::cout << "---- DiskType is <" << dtype << ">" << std::endl; + + if (not sech2) { + switch (DTYPE) { + case DiskType::doubleexpon: + case DiskType::exponential: + case DiskType::diskbulge: + std::cout << "---- pyEXP assumes sech^2(z/(2h)) by default in v7.9.0 and later" << std::endl + << "---- Use the 'sech2: true' in your YAML config to use sech^2(z/(2h))" << std::endl + << "---- This warning will be removed in v7.10.0." << std::endl; + break; + default: + break; } } } - catch (const std::out_of_range& err) { - if (myid==0) { - std::cout << "DiskType error in configuraton file" << std::endl; - std::cout << "Valid options are: "; - for (auto v : dtlookup) std::cout << v.first << " "; - std::cout << std::endl; - } - throw std::runtime_error("Cylindrical::initialize: invalid DiskType"); - } // Check for and initialize the Python density type // @@ -367,6 +379,82 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : pyDens = std::make_shared(pyname); } + // Use these user models to deproject for the EOF spherical basis + // + if (EmpCylSL::mtype == EmpCylSL::EmpModel::Deproject) { + // The scale in EmpCylSL is assumed to be 1 so we compute the + // height relative to the length + // + double H = sech2 ? 0.5*hcyl/acyl : hcyl/acyl; + + // The model instance (you can add others in DiskModels.H). + // It's MN or Exponential if not MN. + // + EmpCylSL::AxiDiskPtr model; + + // Map legacy/short model names to canonical keys expected by dplookup + // + if (dmodel == "exp") { + dmodel = "exponential"; + } + + // Check for map entry + // + try { + PTYPE = dplookup.at(dmodel); + + // Report DeprojType + if (myid==0) { + std::cout << "---- Deprojection type is <" << dmodel + << ">" << std::endl; + } + } + catch (const std::out_of_range& err) { + if (myid==0) { + std::cout << "DeprojType error in configuration file" << std::endl; + std::cout << "Valid options are: "; + for (auto v : dplookup) std::cout << v.first << " "; + std::cout << std::endl; + } + throw std::runtime_error("Cylinder: invalid DiskModel"); + } + + if (PTYPE == DeprojType::mn) // Miyamoto-Nagai + model = std::make_shared(1.0, H); + else if (PTYPE == DeprojType::toomre) { + model = std::make_shared(1.0, H, 5.0); + } else if (PTYPE == DeprojType::python) { + if (pyproj.empty()) { + if (myid==0) { + std::cout << "DeprojType is set to 'python' but no Python " + << "projection module name (pyname/pyproj) was provided." + << std::endl; + } + throw std::runtime_error( + "Cylindrical::initialize: DeprojType 'python' requires a " + "non-empty Python module name (pyname/pyproj)."); + } + model = std::make_shared(pyproj, 1.0); + if (myid==0) + std::cout << "---- Using AxiSymPyModel for deprojection from " + << "Python module <" << pyproj << ">" << std::endl; + } else { // Default to exponential + model = std::make_shared(1.0, H); + } + + if (rwidth>0.0) { + model = std::make_shared(rtrunc/acyl, + rwidth/acyl, + model); + if (myid==0) + std::cout << "Made truncated model with R=" << rtrunc/acyl + << " and W=" << rwidth/acyl << std::endl; + } + + ortho->create_deprojection(H, rfactor, rnum, ncylr, model); + } + + // The conditioning function for the EOF with an optional shift // for M>0 auto dcond = [&](double R, double z, double phi, int M) @@ -542,6 +630,9 @@ void Cylinder::initialize() if (conf["pnum" ]) pnum = conf["pnum" ].as(); if (conf["tnum" ]) tnum = conf["tnum" ].as(); if (conf["ashift" ]) ashift = conf["ashift" ].as(); + if (conf["rwidth" ]) rwidth = conf["rwidth" ].as(); + if (conf["rfactor" ]) rfactor = conf["rfactor" ].as(); + if (conf["rtrunc" ]) rtrunc = conf["rtrunc" ].as(); if (conf["expcond" ]) precond = conf["expcond" ].as(); if (conf["precond" ]) precond = conf["precond" ].as(); if (conf["logr" ]) logarithmic = conf["logr" ].as(); @@ -557,8 +648,11 @@ void Cylinder::initialize() if (conf["cmapr" ]) cmapR = conf["cmapr" ].as(); if (conf["cmapz" ]) cmapZ = conf["cmapz" ].as(); if (conf["vflag" ]) vflag = conf["vflag" ].as(); + if (conf["dtype" ]) dtype = conf["dtype" ].as(); + if (conf["dmodel" ]) dmodel = conf["dmodel" ].as(); if (conf["pyname" ]) pyname = conf["pyname" ].as(); + if (conf["pyproj" ]) pyproj = conf["pyproj" ].as(); if (conf["mtype" ]) mtype = conf["mtype" ].as(); if (conf["ppower" ]) ppow = conf["ppower" ].as(); @@ -2143,7 +2237,7 @@ bool Cylinder::checkMetaData() std::cout << "---- Cylinder::checkMetaData: Python module for disk density has changed since cache creation." << std::endl << "---- Current module: <" << pyname << ">, md5sum: " << current_md5 << std::endl << "---- Loaded module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl - << "---- Cylinder:checkMetaData: forcing cache recomputation to ensure consistency" << std::endl; + << "---- Cylinder::checkMetaData: forcing cache recomputation to ensure consistency" << std::endl; } cache_status = false; } @@ -2152,7 +2246,72 @@ bool Cylinder::checkMetaData() } // End: Python disk type check - // Could add deprojection checks here in the future + // Deprojection consistency checks with cache + if (EmpCylSL::mtype == EmpCylSL::EmpModel::Deproject) { + + // Get the dmodel attribute + // + auto read_attr = file.getAttribute("ProjType"); + std::string loaded_dmodel; + read_attr.read(loaded_dmodel); + + if (loaded_dmodel != dmodel) { + if (myid==0) { + std::cout << "---- Cylinder::checkMetaData: dmodel for cache file <" << cachename << "> is <" + << loaded_dmodel << ">, which does not match the requested dmodel <" + << dmodel << ">" << std::endl + << "---- Cylinder::checkMetaData: forcing cache recomputation" << std::endl; + } + // Force cache recomputation + cache_status = 0; + } + } + + if (cache_status == 1 and dmodel == "python") { + // Get the Python info + // + if (!file.hasAttribute("pythonProjType")) { + // We should not be able to get here since the pythonProjType + // attribute is required for cache creation with the Python + if (myid==0) { + std::cout << "---- Cylinder::checkMetaData: pythonProjType attribute not found in cache file <" << cachename << ">. " << std::endl; + std::cout << "---- Cylinder::checkMetaData: this may be a logic error, trigger recomputation." << std::endl; + } + + cache_status = 0; + + } else { + // Get the pyproj attribute and md5 hash from the cache + auto read_attr = file.getAttribute("pythonProjType"); + std::vector pyinfo; + read_attr.read(pyinfo); + + std::string current_md5; + + // Get the md5sum for requested Python projection module + try { + current_md5 = QuickDigest5::fileToHash(pyproj + ".py"); + } catch (const std::runtime_error& e) { + if (myid==0) + std::cerr << "BiorthBasis::Cylindrical error: " + << e.what() << ", error computing pyproj md5sum" + << std::endl; + } + // Check that the md5sums match for the current Python projection + // + if (current_md5 != pyinfo[1]) { + if (myid==0) { + std::cout << "---- Cylinder::checkMetaData: Python module for deprojection has changed since cache creation." << std::endl + << "---- Current module: <" << pyproj << ">, md5sum: " << current_md5 << std::endl + << "---- Cached module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl + << "---- Cylinder::checkMetaData: forcing cache recomputation to ensure consistency" << std::endl; + } + cache_status = 0; + } + } + // End: DeprojType and Python module consistency checks with cache + } + // End: DiskType and Python module consistency checks with cache } // End: DiskType attribute check @@ -2199,8 +2358,30 @@ void Cylinder::saveMetaData() } } - // Deprojection metadata could be added here + // Deprojection metadata + if (EmpCylSL::mtype == EmpCylSL::EmpModel::Deproject) { + file.createAttribute("ProjType", dmodel); + + if (PTYPE == DeprojType::python) { + try { + std::vector pyinfo = + {pyproj, QuickDigest5::fileToHash(pyproj + ".py")}; + + file.createAttribute("pythonProjType", pyinfo); + + std::cout << "---- Cylinder::saveMetaData: writing pythonProjType <" << pyproj + ".py" + << "> to cache file <" << cachename << ">" << std::endl; + } catch (const std::runtime_error& e) { + if (myid==0) { + std::cerr << "Cylinder::saveMetaData error: " + << e.what() + << ", can not write the pyinfo and md5 hash to HDF5" + << std::endl; + } + } + } + } } catch (const HighFive::Exception& err) { std::cerr << err.what() << std::endl; std::cerr << "Cylinder::saveMetaData: error writing metadata to cache file <"