From 494c619cc1cfc7e7d38bc48cdb8a5593f6d416ff Mon Sep 17 00:00:00 2001 From: mdw Date: Thu, 12 Mar 2026 10:00:36 -0400 Subject: [PATCH 01/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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/44] 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);