diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index f40875985..d52341ab6 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -996,11 +996,12 @@ 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 - // - 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..ca7b01472 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1198,6 +1198,14 @@ namespace BasisClasses {"python", DiskType::python} }; + // Deprojection 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", @@ -1261,7 +1270,8 @@ namespace BasisClasses "playback", "coefCompute", "coefMaster", - "pyname" + "pyname", + "pyproj" }; Cylindrical::Cylindrical(const YAML::Node& CONF) : @@ -1494,7 +1504,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(); @@ -1510,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" ]) pyproj = conf["pyproj" ].as(); if (conf["pcavar"] ) pcavar = conf["pcavar" ].as(); if (conf["subsamp"] ) sampT = conf["subsamp" ].as(); @@ -1697,16 +1708,50 @@ 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); }); + + // 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("Cylindrical::initialize: invalid DiskModel"); + } + + if (PTYPE == DeprojType::mn) // Miyamoto-Nagai model = std::make_shared(1.0, H); - else if (DTYPE == DiskType::python) { - model = std::make_shared(pyname, acyl); + 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; - } - 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, 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 //! A Ferrers Ellipsoid + Evacuated Exponential Disc (semi-realistic //! bar+disk model) @@ -196,6 +197,34 @@ 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 = 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 Q = std::exp(-0.5*std::fabs(z)/h); + double sech = 2.0*Q/(1.0 + Q*Q); // Prevent overflow + return sigma * sech*sech; + } +}; + //! Truncate a AxiDisk class Truncated : public EmpCylSL::AxiDisk { diff --git a/include/EmpCylSL.H b/include/EmpCylSL.H index ebb2af947..28fda1197 100644 --- a/include/EmpCylSL.H +++ b/include/EmpCylSL.H @@ -351,6 +351,12 @@ protected: //! with the new Eigen3 API bool allow_old_cache = false; + //! Abel type for deprojection + enum class AbelType { Derivative, Subtraction, IBP }; + + //! Deprojection method (default: derivative) + AbelType abel_type = AbelType::Derivative; + public: /*! Enum listing the possible selection algorithms for coefficient diff --git a/utils/Test/CMakeLists.txt b/utils/Test/CMakeLists.txt index 6d7ed3a45..123e6120d 100644 --- a/utils/Test/CMakeLists.txt +++ b/utils/Test/CMakeLists.txt @@ -1,5 +1,6 @@ -set(bin_PROGRAMS testBarrier expyaml orthoTest) +set(bin_PROGRAMS testBarrier expyaml orthoTest testDeproj + testEmpDeproj testEmp) set(common_LINKLIB OpenMP::OpenMP_CXX MPI::MPI_CXX expui exputil yaml-cpp ${VTK_LIBRARIES}) @@ -32,6 +33,10 @@ endif() 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(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}) @@ -40,4 +45,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/CubicSpline.H b/utils/Test/CubicSpline.H new file mode 100644 index 000000000..f9fde2f9a --- /dev/null +++ b/utils/Test/CubicSpline.H @@ -0,0 +1,31 @@ +#ifndef CUBIC_SPLINE_H +#define CUBIC_SPLINE_H + +#include + +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..10288b8b6 --- /dev/null +++ b/utils/Test/Deprojector.cc @@ -0,0 +1,204 @@ +#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/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..b43962a5e --- /dev/null +++ b/utils/Test/EmpDeproj.cc @@ -0,0 +1,126 @@ +#include +#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 +#include +#include +#include +#include + +#include "Deprojector.H" +#include "cxxopts.H" + +using namespace Deproject; + +int main(int argc, char* argv[]) +{ + // 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); + + 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 new file mode 100644 index 000000000..c3e8e9309 --- /dev/null +++ b/utils/Test/testEmp.cc @@ -0,0 +1,301 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "EmpDeproj.H" +#include "cxxopts.H" + +// pybind11 embedding +#include +#include +namespace py = pybind11; + +int main(int argc, char* argv[]) +{ + // 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 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") + ("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")) + ("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")) + // 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); + + 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} + }; + + 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()); + } + + // Prepare SigmaZFunc to pass into EmpDeproj + std::function SigmaZFunc; + + // 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; + + // If user supplied a Python module/file and function, embed Python + // and load it here + // + if (result.count("pymodule")) { + std::string pymod = result["pymodule"].as(); + std::string pyfuncname = result["pyfunc"].as(); + + // Start the Python interpreter + pyguard = std::make_unique(); + + py::object py_module; + + // 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; + } + + 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(Rmin + t * (Rmax - Rmin)); + } + + std::ofstream ofs(fname); + 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; + } + + ofs.close(); + std::cout << "Wrote " << fname << std::endl; + + return 0; +} diff --git a/utils/Test/testEmpDeproj.cc b/utils/Test/testEmpDeproj.cc new file mode 100644 index 000000000..54a3a8cc2 --- /dev/null +++ b/utils/Test/testEmpDeproj.cc @@ -0,0 +1,172 @@ +#include +#include +#include +#include +#include +#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 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")) + ("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")) + ("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(abel); + + 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(type); + + 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; +}