diff --git a/.gitmodules b/.gitmodules index 3775929f0..cae935cf6 100644 --- a/.gitmodules +++ b/.gitmodules @@ -8,3 +8,6 @@ [submodule "extern/HighFive"] path = extern/HighFive url = https://github.com/BlueBrain/HighFive.git +[submodule "extern/QuickDigest5"] + path = extern/QuickDigest5 + url = https://github.com/nthnn/QuickDigest5.git diff --git a/CMakeLists.txt b/CMakeLists.txt index e2d49dba5..9ed3fa60c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -259,6 +259,7 @@ execute_process( include_directories(${PROJECT_SOURCE_DIR}/extern/yaml-cpp/include) include_directories(${PROJECT_SOURCE_DIR}/extern/pybind11/include) +include_directories(${PROJECT_SOURCE_DIR}/extern/QuickDigest5/include) # Report to the user message("Configuring build for ${GIT_BRANCH}/${GIT_COMMIT} at ${COMPILE_TIME}") diff --git a/expui/BiorthBasis.H b/expui/BiorthBasis.H index f772ad9ba..7681fc61b 100644 --- a/expui/BiorthBasis.H +++ b/expui/BiorthBasis.H @@ -1071,11 +1071,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; @@ -1084,6 +1085,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 b7723992f..3592a6c67 100644 --- a/expui/BiorthBasis.cc +++ b/expui/BiorthBasis.cc @@ -1,6 +1,7 @@ #include -#include "YamlCheck.H" +#include "quickdigest5.hpp" // for md5 hashing of Python modules +#include "YamlCheck.H" // for YAML configuration checking #include "EXPException.H" #include "BiorthBasis.H" #include "DiskModels.H" @@ -11,6 +12,10 @@ #include #endif +// Suppress HDF5 diagonostic messages from base layer when using +// HighFive. This should be enabled unless one is debugging. +// #define HDF5_QUIET + namespace BasisClasses { const std::set @@ -1180,6 +1185,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", @@ -1215,8 +1228,9 @@ namespace BasisClasses "ashift", "expcond", "ignore", - "deproject", "logr", + "dmodel", + "ppow", "pcavar", "pcaeof", "pcavtk", @@ -1245,6 +1259,7 @@ namespace BasisClasses "coefCompute", "coefMaster", "pyname", + "pyproj", "nint", "totalCovar", "fullCovar" @@ -1391,7 +1406,6 @@ namespace BasisClasses cachename = ".eof_cache_file"; oldcache = false; Ignore = false; - deproject = false; rnum = 200; pnum = 1; @@ -1402,7 +1416,7 @@ namespace BasisClasses EVEN_M = false; cmapR = 1; cmapZ = 1; - mtype = "Exponential"; + mtype = "exponential"; dtype = "exponential"; vflag = 0; @@ -1481,8 +1495,7 @@ namespace BasisClasses if (conf["cmapr" ]) cmapR = conf["cmapr" ].as(); if (conf["cmapz" ]) cmapZ = conf["cmapz" ].as(); if (conf["ignore" ]) Ignore = conf["ignore" ].as(); - if (conf["deproject" ]) deproject = conf["deproject" ].as(); - if (conf["dmodel" ]) dmodel = conf["dmodel" ].as(); + if (conf["dmodel" ]) dmodel = conf["dmodel" ].as(); if (conf["aratio" ]) aratio = conf["aratio" ].as(); if (conf["hratio" ]) hratio = conf["hratio" ].as(); @@ -1493,11 +1506,12 @@ 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(); 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(); @@ -1564,38 +1578,45 @@ namespace BasisClasses EmpCylSL::logarithmic = logarithmic; EmpCylSL::VFLAG = vflag; - // Set deprojected model type + // Convert dmodel string to lower case (deprojection model for EOF + // basis construction) // - // Convert mtype string to lower case - std::transform(mtype.begin(), mtype.end(), mtype.begin(), + std::transform(dmodel.begin(), dmodel.end(), dmodel.begin(), [](unsigned char c){ return std::tolower(c); }); + // Convert mtype string to lower case (EmpCylSL spherical function + // for EOF basis construction) + // + std::transform(mtype.begin(), mtype.end(), mtype.begin(), + [](unsigned char c){ return std::tolower(c); }); + + // Convert dtype string to lower case (disk density function for + // EOF conditioning) + // + std::transform(dtype.begin(), dtype.end(), dtype.begin(), + [](unsigned char c){ return std::tolower(c); }); + // Set EmpCylSL mtype. This is the spherical function used to // generate the EOF basis. If "deproject" is set, this will be - // overriden in EmpCylSL. + // auto-set to the same as the deprojection model. Otherwise, + // this can be set independently to allow for different spherical + // functions for the EOF basis + + // Set EmpCylSL mtype. This is the spherical function used to + // generate the EOF basis. // - EmpCylSL::mtype = EmpCylSL::Exponential; // Default - if (mtype.compare("exponential")==0) { - EmpCylSL::mtype = EmpCylSL::Exponential; + try { + auto itm = EmpCylSL::EmpModelMap.find(mtype); + EmpCylSL::mtype = itm->second; + } + catch (const std::out_of_range& err) { if (myid==0) { - std::cout << "---- Cylindrical: using original exponential deprojected disk for EOF conditioning" << std::endl; - std::cout << "---- Cylindrical: consider using the exact, spherically deprojected exponential with 'mtype: ExpSphere'" << std::endl; + std::cout << "Cylindrical::initialize error parsing 'mtype' parameter in YAML config" << std::endl; + std::cout << "Valid options are: "; + for (auto p : EmpCylSL::EmpModelLabs) std::cout << p.second << " "; + std::cout << "(not case sensitive)" << std::endl; } - } else if (mtype.compare("expsphere")==0) - EmpCylSL::mtype = EmpCylSL::ExpSphere; - else if (mtype.compare("gaussian")==0) - EmpCylSL::mtype = EmpCylSL::Gaussian; - else if (mtype.compare("plummer")==0) - EmpCylSL::mtype = EmpCylSL::Plummer; - else if (mtype.compare("power")==0) { - EmpCylSL::mtype = EmpCylSL::Power; - EmpCylSL::PPOW = ppow; - } else { - if (myid==0) std::cout << "No EmpCylSL EmpModel named <" - << mtype << ">, valid types are: " - << "Exponential, ExpSphere, Gaussian, Plummer, Power " - << "(not case sensitive)" << std::endl; - throw std::runtime_error("Cylindrical:initialize: EmpCylSL bad parameter"); + throw std::runtime_error("Cylindrical::initialize: invalid 'mtype' parameter in YAML config"); } // Check for non-null cache file name. This must be specified @@ -1630,13 +1651,225 @@ namespace BasisClasses if (oldcache) sl->AllowOldCache(); + // Temporary muting of HDF5 error messages for EOF cache reading +#ifdef HDF5_QUIET + H5E_auto2_t old_func; + void* old_client_data; + H5Eget_auto2(H5E_DEFAULT, &old_func, &old_client_data); + + // Mute HDF5 error messages + H5Eset_auto2(H5E_DEFAULT, NULL, NULL); + + // For unmute, use: + // H5Eset_auto2(H5E_DEFAULT, old_func, old_client_data); +#endif + // Attempt to read EOF cache // - if (sl->read_cache() == 0) { + int cache_status = sl->read_cache(); + + std::map disk_type = { + {DiskType::constant, "constant"}, + {DiskType::gaussian, "gaussian"}, + {DiskType::mn, "mn"}, + {DiskType::exponential, "exponential"}, + {DiskType::doubleexpon, "doubleexpon"}, + {DiskType::diskbulge, "diskbulge"}, + {DiskType::python, "python"} + }; + + // Check for map entry, will throw if the key is not in the map. + DTYPE = dtlookup.at(dtype); + + // Checking for cache consistency with current DiskType and Python + // module (if applicable) + // + if (cache_status == 1) { + + try { + // Open the cache file for reading the Python metadata + // + HighFive::File file(cachename, HighFive::File::ReadOnly); + + // mtype is already cached as "model" attribute in the HDF5 + // cache file, so we do not need to write it here. We just + // need to write the DiskType and Python metadata (if + // applicable). + // + + // Check that the DiskType for the cache matches the requested + // DiskType If the DiskType attribute is missing from the + // cache, this indicates thet the cache was created with an + // older versions of the code that did not include the + // DiskType attribute in the cache. We should trigger a cache + // recomputation in this case to be safe, but will let this + // pass for backward compatibility with older caches. + // + if (!file.hasAttribute("DiskType")) { + if (myid==0) { + std::cout << "---- Cylindrical: DiskType attribute not found in cache file <" << cachename << ">. " << std::endl + << "---- This indicates that the cache was created with an older version of the" << std::endl + << "---- code that did not include the DiskType attribute in the cache. We" << std::endl + << "---- suggest deleting the old cache file to prevent this warning and rebuild" << std::endl + << "---- the cache, but we will continue..." << std::endl; + } else { + // Open existing DiskType attribute + // + auto read_attr = file.getAttribute("DiskType"); + + std::string loaded_dtype; + read_attr.read(loaded_dtype); + + DiskType disktype = dtlookup.at(loaded_dtype); + + if (disktype != DTYPE) { + if (myid==0) { + std::cout << "---- Cylindrical: DiskType for cache file <" + << cachename << "> is <" + << loaded_dtype << ">," << std::endl + << "which does not match the requested DiskType <" + << dtype << ">" << std::endl + << "---- Cylindrical: forcing cache recomputation" + << std::endl; + } + // Force cache recomputation + cache_status = 0; + } + else if (disktype == DiskType::python) { + + // If the DiskType is python, we also need to check that + // the Python module used to create the cache matches the + // current Python module to ensure consistency. This tag + // should be present in the cache if it was created with a + // recent version of the code, but if it is missing we + // will trigger cache recomputation to ensure consistency + // with the current Python module. + if (file.hasAttribute("pythonDiskType") == false) { + if (myid==0) { + std::cout << "---- Cylindrical: pythonDiskType attribute not found in cache file <" << cachename << ">. " << std::endl; + std::cout << "---- Cylindrical: this may be a logic error, trigger recomputation." << std::endl; + } + cache_status = 0; + } else { + // Get the pyname attribute + auto read_attr = file.getAttribute("pythonDiskType"); + + std::vector pyinfo; + read_attr.read(pyinfo); + + std::string current_md5; + + // Get the md5sum for requested Python module + try { + current_md5 = QuickDigest5::fileToHash(pyname + ".py"); + } catch (const std::runtime_error& e) { + if (myid==0) + std::cerr << "BiorthBasis::Cylindrical error: " + << e.what() << ", error compuing pyname md5sum" + << std::endl; + } + + // Check that the md5sums match for the current Python + // module and the Python module used to create the + // cache. If they do not match, force cache + // recomputation to ensure consistency with the current + // Python module. + if (current_md5 != pyinfo[1]) { + if (myid==0) { + std::cout << "---- Cylindrical: Python module for disk density has changed since cache creation." << std::endl + << "---- Current module: <" << pyname << ">, md5sum: " << current_md5 << std::endl + << "---- Cached module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl + << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; + } + cache_status = 0; + } + } + } + } + + // Get the deproject attribute + // + if (cache_status == 1 and + EmpCylSL::mtype == EmpCylSL::EmpModel::Deproject) { + + // Get the dmodel attribute + // + auto read_attr = file.getAttribute("ProjType"); + std::string loaded_dmodel; + read_attr.read(loaded_dmodel); + + if (loaded_dmodel != dmodel) { + if (myid==0) { + std::cout << "---- Cylindrical: dmodel for cache file <" << cachename << "> is <" + << loaded_dmodel << ">, which does not match the requested dmodel <" + << dmodel << ">" << std::endl + << "---- Cylindrical: forcing cache recomputation" << std::endl; + } + // Force cache recomputation + cache_status = 0; + } + } + + if (cache_status == 1 and dmodel == "python") { + // Get the Python info + // + if (!file.hasAttribute("pythonProjType")) { + // We should not be able to get here since the pythonProjType + // attribute is required for cache creation with the Python + if (myid==0) { + std::cout << "---- Cylindrical: pythonProjType attribute not found in cache file <" << cachename << ">. " << std::endl; + std::cout << "---- Cylindrical: this may be a logic error, trigger recomputation." << std::endl; + } + + cache_status = 0; + + } else { + // Get the pyproj attribute and md5 hash from the cache + auto read_attr = file.getAttribute("pythonProjType"); + std::vector pyinfo; + read_attr.read(pyinfo); + + std::string current_md5; + + // Get the md5sum for requested Python projection module + try { + current_md5 = QuickDigest5::fileToHash(pyproj + ".py"); + } catch (const std::runtime_error& e) { + if (myid==0) + std::cerr << "BiorthBasis::Cylindrical error: " + << e.what() << ", error computing pyproj md5sum" + << std::endl; + } + // Check that the md5sums match for the current Python projection + // + if (current_md5 != pyinfo[1]) { + if (myid==0) { + std::cout << "---- Cylindrical: Python module for deprojection has changed since cache creation." << std::endl + << "---- Current module: <" << pyproj << ">, md5sum: " << current_md5 << std::endl + << "---- Cached module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl + << "---- Cylindrical: forcing cache recomputation to ensure consistency" << std::endl; + } + cache_status = 0; + } + } + } + // End: DeprojType and Python module consistency checks with cache + } + // End: DiskType and Python module consistency checks with cache + } + catch (const HighFive::Exception& err) { + if (myid==0) { + std::cerr << "---- Cylindrical: " << err.what() << std::endl; + std::cerr << "---- Cylindrical: forcing cache recomputation" << std::endl; + } + cache_status = 0; // Fallback... + } + } + + if (cache_status == 0) { // Remake cylindrical basis // - if (Ignore and myid==0) { std::cout << "---- BasisFactory: We have deprecated the 'ignore' parameter for the" << std::endl << "---- Cylindrical class. If the cache file exists but does" << std::endl @@ -1646,17 +1879,10 @@ namespace BasisClasses << "---- remove 'ignore' from your YAML configuration." << std::endl; } - // Convert dtype string to lower case - // - std::transform(dtype.begin(), dtype.end(), dtype.begin(), - [](unsigned char c){ return std::tolower(c); }); - // Set DiskType. This is the functional form for the disk used to // condition the basis. // - try { // Check for map entry, will through if the - DTYPE = dtlookup.at(dtype); // key is not in the map. - + try { if (myid==0) { // Report DiskType std::cout << "---- DiskType is <" << dtype << ">" << std::endl; @@ -1693,7 +1919,7 @@ namespace BasisClasses // Use these user models to deproject for the EOF spherical basis // - if (deproject) { + if (EmpCylSL::mtype == EmpCylSL::EmpModel::Deproject) { // The scale in EmpCylSL is assumed to be 1 so we compute the // height relative to the length // @@ -1704,16 +1930,56 @@ namespace BasisClasses // EmpCylSL::AxiDiskPtr model; - if (dmodel.compare("MN")==0) // Miyamoto-Nagai - model = std::make_shared(1.0, H); - else if (DTYPE == DiskType::python) { - model = std::make_shared(pyname, acyl); - std::cout << "Using AxiSymPyModel for deprojection from Python function <" - << pyname << ">" << std::endl; + // Map legacy/short model names to canonical keys expected by dplookup + // + if (dmodel == "exp") { + dmodel = "exponential"; } - else // Default to exponential - model = std::make_shared(1.0, H); + // 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 (PTYPE == DeprojType::toomre) { + model = std::make_shared(1.0, H, 5.0); + } else if (PTYPE == DeprojType::python) { + if (pyproj.empty()) { + if (myid==0) { + std::cout << "DeprojType is set to 'python' but no Python " + << "projection module name (pyname/pyproj) was provided." + << std::endl; + } + throw std::runtime_error( + "Cylindrical::initialize: DeprojType 'python' requires a " + "non-empty Python module name (pyname/pyproj)."); + } + model = std::make_shared(pyproj, 1.0); + if (myid==0) + std::cout << "---- Using AxiSymPyModel for deprojection from " + << "Python module <" << pyproj << ">" << std::endl; + } else { // Default to exponential + model = std::make_shared(1.0, H); + } + if (rwidth>0.0) { model = std::make_shared(rtrunc/acyl, rwidth/acyl, @@ -1735,8 +2001,84 @@ namespace BasisClasses }; sl->generate_eof(rnum, pnum, tnum, f); + + if (myid==0) { + try { + // Open the cache file for writing the Python metadata + // + HighFive::File file(cachename, HighFive::File::ReadWrite); + + // mtype is already cached as "model" attribute in the HDF5 + // cache file, so we do not need to write it here. We just + // need to write the DiskType and Python metadata (if + // applicable). + // + file.createAttribute("DiskType", dtype); + + std::cout << "---- Cylindrical: writing DiskType <" << dtype + << "> to cache file <" << cachename << ">" << std::endl; + + // Write the md5sum for the Python module + if (DTYPE == DiskType::python) { + try { + std::vector pyinfo = + {pyname, QuickDigest5::fileToHash(pyname + ".py")}; + + file.createAttribute("pythonDiskType", pyinfo); + + std::cout << "---- Cylindrical: writing pythonDiskType <" << pyname + ".py" + << "> to cache file <" << cachename << ">" << std::endl; + + + } catch (const std::runtime_error& e) { + if (myid==0) { + std::cerr << "BiorthBasis::Cylindrical error: " + << e.what() + << ", can not write the pyname and md5 hash to HDF5" + << std::endl; + } + } + } + + if (EmpCylSL::mtype == EmpCylSL::EmpModel::Deproject) { + + file.createAttribute("ProjType", dmodel); + + if (PTYPE == DeprojType::python) { + try { + std::vector pyinfo = + {pyproj, QuickDigest5::fileToHash(pyproj + ".py")}; + + file.createAttribute("pythonProjType", pyinfo); + + std::cout << "---- Cylindrical: writing pythonProjType <" << pyproj + ".py" + << "> to cache file <" << cachename << ">" << std::endl; + } catch (const std::runtime_error& e) { + if (myid==0) { + std::cerr << "BiorthBasis::Cylindrical error: " + << e.what() + << ", can not write the pyinfo and md5 hash to HDF5" + << std::endl; + } + } + } + } + } catch (const HighFive::Exception& err) { + if (myid==0) { + std::cerr << err.what() << std::endl; + std::cerr << "Error writing metadata to cache file <" << cachename + << std::endl; + } + } + // Errors will prevent metadata from being written to the cache + } + // Only the root process should be updating the cache } +#ifdef HDF5_QUIET + // Unmute HDF5 error messages + H5Eset_auto2(H5E_DEFAULT, old_func, old_client_data); +#endif // Orthogonality sanity check // if (myid==0) orthoTest(); diff --git a/exputil/CMakeLists.txt b/exputil/CMakeLists.txt index 1124b101e..3f8862a7b 100644 --- a/exputil/CMakeLists.txt +++ b/exputil/CMakeLists.txt @@ -14,7 +14,7 @@ set(UTIL_SRC nrutil.cc elemfunc.cc euler.cc euler_slater.cc TransformFFT.cc QDHT.cc YamlCheck.cc ExpDeproj.cc # Hankel.cc parseVersionString.cc EXPmath.cc laguerre_polynomial.cpp YamlConfig.cc orthoTest.cc OrthoFunction.cc VtkGrid.cc - Sutils.cc fpetrap.cc) + Sutils.cc fpetrap.cc getmd5sum.cc) if(HAVE_VTK) list(APPEND UTIL_SRC VtkPCA.cc) @@ -42,19 +42,21 @@ set(SLEDGE_SRC sledge.f) set(PARTICLE_SRC Particle.cc ParticleReader.cc header.cc) set(CUDA_SRC cudaParticle.cu cudaSLGridMP2.cu) set(PYWRAP_SRC DiskDensityFunc.cc) +set(QD5_SRC ${PROJECT_SOURCE_DIR}/extern/QuickDigest5/src/quickdigest5.cpp) -set(exputil_SOURCES ${ODE_SRC} ${ROOT_SRC} ${QUAD_SRC} - ${RANDOM_SRC} ${UTIL_SRC} ${SPECFUNC_SRC} - ${PHASE_SRC} ${SYMP_SRC} ${INTERP_SRC} ${MASSMODEL_SRC} - ${ORBIT_SRC} ${BIORTH_SRC} ${POLY_SRC} ${GAUSS_SRC} - ${QPDISTF_SRC} ${BESSEL_SRC} ${OPTIMIZATION_SRC} - ${SLEDGE_SRC} ${PARTICLE_SRC} ${CUDA_SRC} ${PYWRAP_SRC}) + +set(exputil_SOURCES ${ODE_SRC} ${ROOT_SRC} ${QUAD_SRC} ${RANDOM_SRC} + ${UTIL_SRC} ${SPECFUNC_SRC} ${PHASE_SRC} ${SYMP_SRC} ${INTERP_SRC} + ${MASSMODEL_SRC} ${ORBIT_SRC} ${BIORTH_SRC} ${POLY_SRC} ${GAUSS_SRC} + ${QPDISTF_SRC} ${BESSEL_SRC} ${OPTIMIZATION_SRC} ${SLEDGE_SRC} + ${PARTICLE_SRC} ${CUDA_SRC} ${PYWRAP_SRC} ${QD5_SRC}) set(common_INCLUDE_DIRS $ $ $ ${CMAKE_BINARY_DIR} $ $ + $ ${DEP_INC} ${EIGEN3_INCLUDE_DIR} ${HDF5_INCLUDE_DIRS} ${FFTW_INCLUDE_DIRS}) diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc index 2b84bee38..7b69bdc08 100644 --- a/exputil/EmpCylSL.cc +++ b/exputil/EmpCylSL.cc @@ -70,17 +70,25 @@ double EmpCylSL::PPOW = 4.0; bool EmpCylSL::NewCoefs = true; -EmpCylSL::EmpModel EmpCylSL::mtype = Exponential; - -std::map EmpCylSL::EmpModelLabs = - { {Exponential, "Exponential"}, - {ExpSphere, "ExpSphere" }, - {Gaussian, "Gaussian" }, - {Plummer, "Plummer" }, - {Power, "Power" }, - {Deproject, "Deproject" } - }; +EmpCylSL::EmpModel EmpCylSL::mtype = EmpCylSL::EmpModel::Exponential; + +std::map EmpCylSL::EmpModelMap = { + {"exponential", EmpModel::Exponential}, + {"expsphere", EmpModel::ExpSphere}, + {"gaussian", EmpModel::Gaussian}, + {"plummer", EmpModel::Plummer}, + {"power", EmpModel::Power}, + {"deproject", EmpModel::Deproject} +}; +std::map EmpCylSL::EmpModelLabs = { + {EmpModel::Exponential, "Exponential"}, + {EmpModel::ExpSphere, "ExpSphere"}, + {EmpModel::Gaussian, "Gaussian"}, + {EmpModel::Plummer, "Plummer"}, + {EmpModel::Power, "Power"}, + {EmpModel::Deproject, "Deproject"} +}; EmpCylSL::EmpCylSL() { @@ -487,31 +495,48 @@ void EmpCylSL::create_deprojection(double H, double Rf, int NUMR, int NINT, Linear1d surf(rl, sigI); - // Now, compute Abel inverion integral + // Now, compute Abel inversion integral // for (int i=0; i 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=RMAX) ans = massRg.eval(RMAX); else ans = massRg.eval(log(R)); @@ -591,21 +616,21 @@ double EmpCylSL::densR(double R) double ans=0.0, fac, arg; switch (mtype) { - case Exponential: + case EmpModel::Exponential: ans = exp(-R)/(4.0*M_PI*R); break; - case ExpSphere: + case EmpModel::ExpSphere: ans = expdeproj.density(R); break; - case Gaussian: + case EmpModel::Gaussian: arg = 0.5*R*R; ans = exp(-arg)/(4.0*M_PI*R); break; - case Plummer: + case EmpModel::Plummer: fac = 1.0/(1.0+R); ans = 3.0*pow(fac, 4.0)/(4.0*M_PI); break; - case Power: + case EmpModel::Power: { double z = R + 1.0; double a1 = PPOW - 1.0; @@ -614,7 +639,7 @@ double EmpCylSL::densR(double R) ans = 0.125*a1*a2*a3/M_PI * pow(z, -PPOW); } break; - case Deproject: + case EmpModel::Deproject: if (R < RMIN) ans = densRg.eval(RMIN); else if (R>=RMAX) ans = 0.0; else ans = densRg.eval(log(R)); diff --git a/exputil/getmd5sum.cc b/exputil/getmd5sum.cc new file mode 100644 index 000000000..b7c234175 --- /dev/null +++ b/exputil/getmd5sum.cc @@ -0,0 +1,49 @@ +// Initial reference implementation for testing. Production code uses +// QuickDigest5 as a git submodule. + +#include +#include +#include +#include +#include +#include + +// Define the custom deleter +struct PcloseDeleter { + void operator()(FILE* f) const { + if (f) pclose(f); + } +}; + +std::string get_md5sum(const std::string& filename) +{ + // Check if md5sum is available on the system + if (std::system("which md5sum > /dev/null 2>&1") != 0) { + throw std::runtime_error("md5sum command not found. Please ensure it is installed and in your PATH."); + } + + // Command to execute: md5sum + std::string command = "md5sum " + filename; + std::array buffer; + std::string result = ""; + + // Use popen to execute the command and read its output + // "r" mode opens the pipe for reading + std::unique_ptr pipe(popen(command.c_str(), "r")); + if (!pipe) { + throw std::runtime_error("popen() failed!"); + } + + // Read the output line by line + while (fgets(buffer.data(), buffer.size(), pipe.get()) != nullptr) { + result += buffer.data(); + } + + // The stanard GNU/Linux md5sum output format is: "32-char-hash + // filename". We extract the 32-character hash part... + if (result.length() >= 32) { + return result.substr(0, 32); + } else { + throw std::runtime_error("Failed to parse md5sum output."); + } +} diff --git a/extern/QuickDigest5 b/extern/QuickDigest5 new file mode 160000 index 000000000..1d61aece0 --- /dev/null +++ b/extern/QuickDigest5 @@ -0,0 +1 @@ +Subproject commit 1d61aece0e023e4c41cf24547fa1ac2f8028219d diff --git a/include/DiskModels.H b/include/DiskModels.H index 3b46da803..4b409c392 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) @@ -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 9811a0059..cbbe7f308 100644 --- a/include/EmpCylSL.H +++ b/include/EmpCylSL.H @@ -243,7 +243,9 @@ protected: bool ReadH5Cache(); //! Cache versioning - inline static const std::string Version = "1.0"; + // Version 1.0 corresponds to bases generated with the first public release of EXP (7.8). + // Version 1.1 corresponds to bases generated with improved Abel inversion and deprojection methods (7.10). + inline static const std::string Version = "1.1"; //! The cache file name std::string cachefile; @@ -355,6 +357,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 @@ -366,7 +374,7 @@ public: }; //! Type of density model to use - enum EmpModel { + enum class EmpModel { Exponential, ExpSphere, Gaussian, @@ -375,6 +383,9 @@ public: Deproject, }; + static std::map EmpModelMap; + static std::map EmpModelLabs; + //! Axisymmetric disk density function for deprojection class AxiDisk { @@ -480,9 +491,6 @@ public: //! Use YAML header in coefficient file static bool NewCoefs; - //! Convert EmpModel to ascii label - static std::map EmpModelLabs; - //! Fraction of table range for basis images (for debug) static double HFAC; diff --git a/include/exputils.H b/include/exputils.H index 986aed037..c530e8e68 100644 --- a/include/exputils.H +++ b/include/exputils.H @@ -12,4 +12,8 @@ orthoCompute(const std::vector& tests); void orthoTest(const std::vector& tests, const std::string& classname, const std::string& indexname); +// Signature for md5sum function computation +std::string get_md5sum(const std::string&); + + #endif diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index fa7ea2057..87df43c8f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -25,6 +25,7 @@ set(common_INCLUDE_DIRS $ $ $ + $ ${CMAKE_BINARY_DIR} ${DEP_INC} ${CMAKE_CURRENT_SOURCE_DIR} ${EIGEN3_INCLUDE_DIR}) diff --git a/src/Cylinder.H b/src/Cylinder.H index 897427d01..6103321a5 100644 --- a/src/Cylinder.H +++ b/src/Cylinder.H @@ -7,6 +7,7 @@ #include "Basis.H" #include "CylEXP.H" #include "Coefficients.H" +#include "DiskDensityFunc.H" #include "config_exp.h" @@ -39,6 +40,24 @@ class MixtureBasis; @param hcyl is the scale height + @param aratio is the scale length ratio for disk basis construction (default: 1.0) + + @param hratio is the scale height ratio for disk basis construction (default: 1.0) + + @param dweight is the relative weight of the second component in the double-exponential disk basis construction (default: 1.0) + + @param dmodel is the disk model type for deprojection and EOF conditioning. See Cylinder::DiskType for options. + + @param mtype is the spherical model type for conditioning the basis. See EmpCylSL::EmpModel for options. + + @param dtype is the density model type for conditioning the basis. See Cylinder::DiskType for options. + + @param rwidth is the radial width of the error function truncation for disk basis construction (default: 0.0, which implies no truncation) + + @param rtrunc is the radial truncation radius for disk basis construction (default: 0.1) + + @param rfactor is the radial scale factor for disk basis construction (default: 1.0, which implies no scaling) + @param lmaxfid is the maximum spherical harmonic index for EOF construction @param nmaxfid is the maximum radial index for EOF construction @@ -114,7 +133,7 @@ class MixtureBasis; @param cmapz is the vertical coordinate mapping type - @param mtype is the deprojected model type for conditioning the + @param mtype is the spherical model type for conditioning the basis. See ExpCylSL::mtype for options. @param ppower is the power-law index for the Power spherical model @@ -159,8 +178,13 @@ private: void initialize(void); - // Parameters + // Parameters + // double rcylmin, rcylmax, zmax, acyl, bias; + double aratio, hratio, dweight; // Disk basis construction parameters for double-exponential disk + double rwidth, rtrunc, rfactor; // Disk basis construction parameters for numerical deprojection + double Mfac, HERNA; // Disk bulge parameters: mass fraction and Hernquist scale length + int nmaxfid, lmaxfid, mmax, mlim, nint; int ncylnx, ncylny, ncylr; double hcyl, hexp, snr, rem; @@ -284,7 +308,6 @@ protected: //@} - #endif /** Test change level counts for deep debugging enabled by setting @@ -345,12 +368,84 @@ protected: //! Power-law index for deprojected spherical model for basis conditioning double ppow = 2.0; - //! Deprojection model for basis conditioning + //! Spherical model for basis conditioning std::string mtype = "Exponential"; //! Write basis-specific parameters to HDF5 covariance file void writeCovarH5Params(HighFive::File& file); + //! Disk density function types for basis conditioning + enum class DiskType + { constant, gaussian, mn, exponential, doubleexpon, diskbulge, python }; + + //! The string name of the deprojection model + std::string dmodel; + + //! Map human strings to disk target typenums + const std::map dtlookup = + { {"constant", DiskType::constant}, + {"gaussian", DiskType::gaussian}, + {"mn", DiskType::mn}, + {"exponential", DiskType::exponential}, + {"doubleexpon", DiskType::doubleexpon}, + {"diskbulge", DiskType::diskbulge}, + {"python", DiskType::python} + }; + + //! Disk target type strings for reporting + const std::map dtstring = + { {DiskType::constant, "constant"}, + {DiskType::gaussian, "gaussian"}, + {DiskType::mn, "mn"}, + {DiskType::exponential, "exponential"}, + {DiskType::doubleexpon, "doubleexpon"}, + {DiskType::diskbulge, "diskbulge"}, + {DiskType::python, "python"} + }; + + //! The default disk density function type for basis conditioning + DiskType DTYPE = DiskType::exponential; + + //! The string name of the disk density function type for basis conditioning + std::string dtype; + + //@{ + //! DeprojType support + + //! Disk models used for deprojection + enum class DeprojType + { mn, toomre, python, exponential}; + + //! Current model + DeprojType PTYPE; + + //! Python module name for disk density function if using Python deprojection + std::string pyproj; + + //! Look up by string + const std::map dplookup = + { {"mn", DeprojType::mn}, + {"toomre", DeprojType::toomre}, + {"python", DeprojType::python}, + {"exponential", DeprojType::exponential} + }; + //@} + + + //! Dtype cache check for consistency between cache and current parameters + bool checkMetaData(); + + //! Reopen and write the current Dtype into the cache metadata + void saveMetaData(); + + //! Disk density functions for basis conditioning + double DiskDens(double R, double z, double phi); + + //! If using a Python disk density function, this is the instance of + //! the wrapper class that allows us to call the Python function + //! from C++ + std::shared_ptr pyDens; + public: //! Mutexes for multithreading diff --git a/src/Cylinder.cc b/src/Cylinder.cc index 63a5ca1d0..50aa2dd8f 100644 --- a/src/Cylinder.cc +++ b/src/Cylinder.cc @@ -11,9 +11,11 @@ #include "Cylinder.H" #include "MixtureBasis.H" #include "DiskDensityFunc.H" -#include "Timer.H" -#include "exputils.H" -#include "NVTX.H" +#include "Timer.H" // for timing code sections +#include "exputils.H" // utility functions +#include "NVTX.H" // for NVTX profiling of CUDA code +#include "quickdigest5.hpp" // for md5 hashing of Python modules +#include "DiskModels.H" // //@{ //! These are for testing exclusively (should be set false for production) @@ -32,6 +34,11 @@ Cylinder::valid_keys = { "hcyl", "sech2", "hexp", + "aratio", + "hratio", + "dweight", + "Mfac", + "HERNA", "snr", "evcut", "nmaxfid", @@ -55,6 +62,9 @@ Cylinder::valid_keys = { "pnum", "tnum", "ashift", + "rwidth", + "rtrunc", + "rfactor", "expcond", "precond", "logr", @@ -77,6 +87,7 @@ Cylinder::valid_keys = { "playback", "coefCompute", "coefMaster", + "dtype", "pyname", "dumpbasis", "fullCovar", @@ -120,11 +131,27 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : ncylodd = 9; ncylrecomp = -1; + mtype = "exponential"; + dtype = "exponential"; + dmodel = "EXP"; + + // For disk basis construction with doubleexpon + // + aratio = 1.0; // Radial scale length ratio + hratio = 1.0; // Vertical scale height ratio + dweight = 1.0; // Ratio of second disk relative to the first disk + + // For disk bulge + Mfac = 1.0; // mass fraction for disk + HERNA = 0.10; // Hernquist scale a disk bulge + rnum = 200; pnum = 1; tnum = 80; ashift = 0.0; - + rwidth = 0.0; // Width of error function truncation (ignored if zero) + rfactor = 1.0; // Radial scale factor for numerical basis construction + rtrunc = 0.1; // Radial truncation for numerical basis construction vflag = 0; eof = 1; npca = std::numeric_limits::max(); @@ -181,45 +208,58 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : if (cachename.size()==0) throw std::runtime_error("EmpCylSL: you must specify a cachename"); - - // Set deprojected model type - // // Convert mtype string to lower case + // std::transform(mtype.begin(), mtype.end(), mtype.begin(), [](unsigned char c){ return std::tolower(c); }); // Set EmpCylSL mtype. This is the spherical function used to - // generate the EOF basis. + // generate the EOF basis. If "deproject" is set, this will be + // overriden in EmpCylSL. // - EmpCylSL::mtype = EmpCylSL::Exponential; // Default - if (mtype.compare("exponential")==0) { - EmpCylSL::mtype = EmpCylSL::Exponential; + auto itm = EmpCylSL::EmpModelMap.find(mtype); + + if (itm == EmpCylSL::EmpModelMap.end()) { if (myid==0) { - std::cout << "---- Cylindrical: using original exponential deprojected disk for EOF conditioning" << std::endl; - std::cout << "---- Cylindrical: consider using the exact, spherically deprojected exponential with 'mtype: ExpSphere'" << std::endl; + std::cout << "No EmpCylSL EmpModel named <" + << mtype << ">, valid types are: "; + for (auto p : EmpCylSL::EmpModelLabs) std::cout << p.second << " "; + std::cout << "(not case sensitive)" << std::endl; } - } else if (mtype.compare("expsphere")==0) - EmpCylSL::mtype = EmpCylSL::ExpSphere; - else if (mtype.compare("gaussian")==0) - EmpCylSL::mtype = EmpCylSL::Gaussian; - else if (mtype.compare("plummer")==0) - EmpCylSL::mtype = EmpCylSL::Plummer; - else if (mtype.compare("power")==0) { - EmpCylSL::mtype = EmpCylSL::Power; - EmpCylSL::PPOW = ppow; - } else { - if (myid==0) std::cout << "No EmpCylSL EmpModel named <" - << mtype << ">, valid types are: " - << "Exponential, ExpSphere, Gaussian, Plummer, Power " - << "(not case sensitive)" << std::endl; - throw std::runtime_error("Cylinder:initialize: EmpCylSL bad parameter"); + throw std::runtime_error("Cylindrical:initialize: EmpCylSL bad parameter"); } + + EmpCylSL::mtype = itm->second; // Make the empirical orthogonal basis instance // ortho = std::make_shared (nmaxfid, lmaxfid, mmax, nmax, bias*acyl, hcyl, ncylodd, cachename); + // Convert dtype string to lower case + // + std::transform(dtype.begin(), dtype.end(), dtype.begin(), + [](unsigned char c){ return std::tolower(c); }); + + // Convert dmodel to lower case + // + std::transform(dmodel.begin(), dmodel.end(), dmodel.begin(), + [](unsigned char c){ return std::tolower(c); }); + + // Check for map entry, will throw if the key is not in the map. + try { + DTYPE = dtlookup.at(dtype); + } + catch (const std::out_of_range& err) { + if (myid==0) { + std::cout << "DiskType error in configuraton file" << std::endl; + std::cout << "Valid options are: "; + for (auto v : dtlookup) std::cout << v.first << " "; + std::cout << std::endl; + } + throw std::runtime_error("Cylinder: invalid DiskType"); + } + // Set azimuthal harmonic order restriction? // if (mlim>=0) ortho->set_mlim(mlim); @@ -258,6 +298,10 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : if (std::filesystem::exists(cachename)) { cache_ok = ortho->read_cache(); + + if (cache_ok) { + cache_ok = checkMetaData(); + } // If new cache is requested, backup existing basis file // @@ -307,19 +351,109 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : << "this step will take many minutes." << std::endl; - std::shared_ptr dfunc; - if (pyname.size()) dfunc = std::make_shared(pyname); - // Use either the user-supplied Python target or the default - // exponential disk model - auto DiskDens = [&](double R, double z, double phi) - { - if (dfunc) return (*dfunc)(R, z, phi); - double h = sech2 ? 0.5*hcyl : hcyl; - double f = exp(-fabs(z)/h); // Overflow prevention - double s = 2.0*f/(1.0 + f*f); // in sech computation - return exp(-R/acyl)*s*s/(4.0*M_PI*acyl*acyl*h); - }; + // Warning about DiskType and pyEXP assumptions for backward compatibility with + // previous versions of pyEXP. + // + if (myid==0) { + std::cout << "---- DiskType is <" << dtype << ">" << std::endl; + + if (not sech2) { + switch (DTYPE) { + case DiskType::doubleexpon: + case DiskType::exponential: + case DiskType::diskbulge: + std::cout << "---- pyEXP assumes sech^2(z/(2h)) by default in v7.9.0 and later" << std::endl + << "---- Use the 'sech2: true' in your YAML config to use sech^2(z/(2h))" << std::endl + << "---- This warning will be removed in v7.10.0." << std::endl; + break; + default: + break; + } + } + } + + // Check for and initialize the Python density type + // + if (DTYPE == DiskType::python) { + pyDens = std::make_shared(pyname); + } + + // Use these user models to deproject for the EOF spherical basis + // + if (EmpCylSL::mtype == EmpCylSL::EmpModel::Deproject) { + // The scale in EmpCylSL is assumed to be 1 so we compute the + // height relative to the length + // + double H = sech2 ? 0.5*hcyl/acyl : hcyl/acyl; + + // The model instance (you can add others in DiskModels.H). + // It's MN or Exponential if not MN. + // + EmpCylSL::AxiDiskPtr model; + + // Map legacy/short model names to canonical keys expected by dplookup + // + if (dmodel == "exp") { + dmodel = "exponential"; + } + + // Check for map entry + // + try { + PTYPE = dplookup.at(dmodel); + + // Report DeprojType + if (myid==0) { + std::cout << "---- Deprojection type is <" << dmodel + << ">" << std::endl; + } + } + catch (const std::out_of_range& err) { + if (myid==0) { + std::cout << "DeprojType error in configuration file" << std::endl; + std::cout << "Valid options are: "; + for (auto v : dplookup) std::cout << v.first << " "; + std::cout << std::endl; + } + throw std::runtime_error("Cylinder: invalid DiskModel"); + } + + if (PTYPE == DeprojType::mn) // Miyamoto-Nagai + model = std::make_shared(1.0, H); + else if (PTYPE == DeprojType::toomre) { + model = std::make_shared(1.0, H, 5.0); + } else if (PTYPE == DeprojType::python) { + if (pyproj.empty()) { + if (myid==0) { + std::cout << "DeprojType is set to 'python' but no Python " + << "projection module name (pyname/pyproj) was provided." + << std::endl; + } + throw std::runtime_error( + "Cylindrical::initialize: DeprojType 'python' requires a " + "non-empty Python module name (pyname/pyproj)."); + } + model = std::make_shared(pyproj, 1.0); + if (myid==0) + std::cout << "---- Using AxiSymPyModel for deprojection from " + << "Python module <" << pyproj << ">" << std::endl; + } else { // Default to exponential + model = std::make_shared(1.0, H); + } + + if (rwidth>0.0) { + model = std::make_shared(rtrunc/acyl, + rwidth/acyl, + model); + if (myid==0) + std::cout << "Made truncated model with R=" << rtrunc/acyl + << " and W=" << rwidth/acyl << std::endl; + } + + ortho->create_deprojection(H, rfactor, rnum, ncylr, model); + } + // The conditioning function for the EOF with an optional shift // for M>0 @@ -348,6 +482,7 @@ Cylinder::Cylinder(Component* c0, const YAML::Node& conf, MixtureBasis *m) : }; ortho->generate_eof(rnum, pnum, tnum, dcond); + saveMetaData(); } firstime = false; @@ -485,10 +620,19 @@ void Cylinder::initialize() if (conf["eof_file" ]) cachename = conf["eof_file" ].as(); if (conf["samplesz" ]) defSampT = conf["samplesz" ].as(); + if (conf["aratio" ]) aratio = conf["aratio" ].as(); + if (conf["hratio" ]) hratio = conf["hratio" ].as(); + if (conf["dweight" ]) dweight = conf["dweight" ].as(); + if (conf["Mfac" ]) Mfac = conf["Mfac" ].as(); + if (conf["HERNA" ]) HERNA = conf["HERNA" ].as(); + if (conf["rnum" ]) rnum = conf["rnum" ].as(); if (conf["pnum" ]) pnum = conf["pnum" ].as(); if (conf["tnum" ]) tnum = conf["tnum" ].as(); if (conf["ashift" ]) ashift = conf["ashift" ].as(); + if (conf["rwidth" ]) rwidth = conf["rwidth" ].as(); + if (conf["rfactor" ]) rfactor = conf["rfactor" ].as(); + if (conf["rtrunc" ]) rtrunc = conf["rtrunc" ].as(); if (conf["expcond" ]) precond = conf["expcond" ].as(); if (conf["precond" ]) precond = conf["precond" ].as(); if (conf["logr" ]) logarithmic = conf["logr" ].as(); @@ -504,6 +648,11 @@ void Cylinder::initialize() if (conf["cmapr" ]) cmapR = conf["cmapr" ].as(); if (conf["cmapz" ]) cmapZ = conf["cmapz" ].as(); if (conf["vflag" ]) vflag = conf["vflag" ].as(); + + if (conf["dtype" ]) dtype = conf["dtype" ].as(); + if (conf["dmodel" ]) dmodel = conf["dmodel" ].as(); + if (conf["pyname" ]) pyname = conf["pyname" ].as(); + if (conf["pyproj" ]) pyproj = conf["pyproj" ].as(); if (conf["mtype" ]) mtype = conf["mtype" ].as(); if (conf["ppower" ]) ppow = conf["ppower" ].as(); @@ -642,6 +791,83 @@ void Cylinder::initialize() } } +double Cylinder::DiskDens(double R, double z, double phi) +{ + double ans = 0.0; + + switch (DTYPE) { + + case DiskType::constant: + if (R < acyl && fabs(z) < hcyl) + ans = 1.0/(2.0*hcyl*M_PI*acyl*acyl); + break; + + case DiskType::gaussian: + if (fabs(z) < hcyl) + ans = 1.0/(2.0*hcyl*2.0*M_PI*acyl*acyl)* + exp(-R*R/(2.0*acyl*acyl)); + break; + + case DiskType::mn: + { + double Z2 = z*z + hcyl*hcyl; + double Z = sqrt(Z2); + double Q2 = (acyl + Z)*(acyl + Z); + ans = 0.25*hcyl*hcyl/M_PI*(acyl*R*R + (acyl + 3.0*Z)*Q2)/( pow(R*R + Q2, 2.5) * Z*Z2 ); + } + break; + + case DiskType::doubleexpon: + { + double a1 = acyl; + double a2 = acyl*aratio; + double h1 = hcyl; + double h2 = hcyl*hratio; + double w1 = 1.0/(1.0+dweight); + double w2 = dweight/(1.0+dweight); + + if (sech2) { h1 *= 0.5; h2 *= 0.5; } + + double f1 = cosh(z/h1); + double f2 = cosh(z/h2); + + ans = + w1*exp(-R/a1)/(4.0*M_PI*a1*a1*h1*f1*f1) + + w2*exp(-R/a2)/(4.0*M_PI*a2*a2*h2*f2*f2) ; + } + break; + + case DiskType::diskbulge: + { + double h = sech2 ? 0.5*hcyl : hcyl; + double f = cosh(z/h); + double rr = pow(pow(R, 2) + pow(z,2), 0.5); + double w1 = Mfac; + double w2 = 1.0 - Mfac; + double as = HERNA; + + ans = w1*exp(-R/acyl)/(4.0*M_PI*acyl*acyl*h*f*f) + + w2*pow(as, 4)/(2.0*M_PI*rr)*pow(rr+as,-3.0) ; + } + break; + case DiskType::python: + ans = (*pyDens)(R, z, phi); + break; + case DiskType::exponential: + default: + { + double h = sech2 ? 0.5*hcyl : hcyl; + double f = cosh(z/h); + ans = exp(-R/acyl)/(4.0*M_PI*acyl*acyl*h*f*f); + } + break; + } + + // if (rwidth>0.0) ans *= erf((rtrunc-R)/rwidth); + + return ans; +} + void Cylinder::get_acceleration_and_potential(Component* C) { nvTracerPtr tPtr; @@ -963,6 +1189,10 @@ void Cylinder::determine_coefficients_particles(void) bool cache_ok = false; if (try_cache || restart) { cache_ok = ortho->read_cache(); + + if (cache_ok) { + cache_ok = checkMetaData(); + } // For a restart, cache must be read // otherwise, abort if (restart && !cache_ok) { @@ -1922,3 +2152,239 @@ void Cylinder::writeCovarH5Params(HighFive::File& file) file.createAttribute("hcyl", HighFive::DataSpace::From(hcyl)).write(hcyl); } + +bool Cylinder::checkMetaData() +{ + // All processes must check the dtype metadata to ensure consistency + // + bool cache_status = true; + + // Open the cache file for reading the Python metadata + // + HighFive::File file(cachename, HighFive::File::ReadOnly); + + // Sanity: check that the DiskType attribute exists + // + if (!file.hasAttribute("DiskType")) { + if (myid==0) { + std::cout << "---- Cylinder:checkMetaData: DiskType attribute not found in cache file <" << cachename << ">" << std::endl + << "---- This may indicate an old cache file created before EmpModel metadata was added" << std::endl + << "---- We will continue...but consider remaking the cache to avoid confusion" << std::endl; + } + } + else { + // Open existing DiskType attribute + // + auto read_attr = file.getAttribute("DiskType"); + + std::string loaded_dtype; + read_attr.read(loaded_dtype); + + // Map the loaded dtype string to a DiskType enum value + // + DiskType disktype = dtlookup.at(loaded_dtype); + + if (disktype != DTYPE) { + if (myid==0) { + std::cout << "---- Cylinder::checkMetaData: DiskType for cache file <" + << cachename << "> is <" + << loaded_dtype << ">," << std::endl + << "which does not match the requested DiskType <" + << dtype << ">" << std::endl + << "---- Cylinder::checkMetaData: forcing cache recomputation" + << std::endl; + } + // Force cache recomputation + cache_status = false; + } + else if (disktype == DiskType::python) { + // Sanity check: if DiskType is python, then the + // pythonDiskType attribute must exist + // + if (!file.hasAttribute("pythonDiskType")) { + if (myid==0) { + std::cout << "---- Cylinder::checkMetaData: pythonDiskType attribute not found in cache file <" << cachename << ">. " << std::endl; + std::cout << "---- Cylinder::checkMetaData: this may indicate a logic error. Forcing cache recomputation." << std::endl; + } + // Force cache recomputation + cache_status = false; + } else { + auto read_attr = file.getAttribute("pythonDiskType"); + + // Get the pyname attribute + std::vector pyinfo; + read_attr.read(pyinfo); + + std::string current_md5; + + // Get the md5sum for requested Python module source file + try { + current_md5 = QuickDigest5::fileToHash(pyname); + } catch (const std::runtime_error& e) { + if (myid==0) + std::cerr << "Cylinder::checkMetaData error: " + << e.what() << std::endl; + } + + // Check that the md5sums match for the current Python + // module source files and the loaded Python module used to + // create the cache. If they do not match, force cache + // recomputation to ensure consistency with the current + // Python module. + // + if (current_md5 != pyinfo[1]) { + if (myid==0) { + std::cout << "---- Cylinder::checkMetaData: Python module for disk density has changed since cache creation." << std::endl + << "---- Current module: <" << pyname << ">, md5sum: " << current_md5 << std::endl + << "---- Loaded module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl + << "---- Cylinder::checkMetaData: forcing cache recomputation to ensure consistency" << std::endl; + } + cache_status = false; + } + } + // End: have Python disk type, check md5 hashes + } + // End: Python disk type check + + // Deprojection consistency checks with cache + if (EmpCylSL::mtype == EmpCylSL::EmpModel::Deproject) { + + // Get the dmodel attribute + // + auto read_attr = file.getAttribute("ProjType"); + std::string loaded_dmodel; + read_attr.read(loaded_dmodel); + + if (loaded_dmodel != dmodel) { + if (myid==0) { + std::cout << "---- Cylinder::checkMetaData: dmodel for cache file <" << cachename << "> is <" + << loaded_dmodel << ">, which does not match the requested dmodel <" + << dmodel << ">" << std::endl + << "---- Cylinder::checkMetaData: forcing cache recomputation" << std::endl; + } + // Force cache recomputation + cache_status = 0; + } + } + + if (cache_status == 1 and dmodel == "python") { + // Get the Python info + // + if (!file.hasAttribute("pythonProjType")) { + // We should not be able to get here since the pythonProjType + // attribute is required for cache creation with the Python + if (myid==0) { + std::cout << "---- Cylinder::checkMetaData: pythonProjType attribute not found in cache file <" << cachename << ">. " << std::endl; + std::cout << "---- Cylinder::checkMetaData: this may be a logic error, trigger recomputation." << std::endl; + } + + cache_status = 0; + + } else { + // Get the pyproj attribute and md5 hash from the cache + auto read_attr = file.getAttribute("pythonProjType"); + std::vector pyinfo; + read_attr.read(pyinfo); + + std::string current_md5; + + // Get the md5sum for requested Python projection module + try { + current_md5 = QuickDigest5::fileToHash(pyproj + ".py"); + } catch (const std::runtime_error& e) { + if (myid==0) + std::cerr << "BiorthBasis::Cylindrical error: " + << e.what() << ", error computing pyproj md5sum" + << std::endl; + } + // Check that the md5sums match for the current Python projection + // + if (current_md5 != pyinfo[1]) { + if (myid==0) { + std::cout << "---- Cylinder::checkMetaData: Python module for deprojection has changed since cache creation." << std::endl + << "---- Current module: <" << pyproj << ">, md5sum: " << current_md5 << std::endl + << "---- Cached module: <" << pyinfo[0] << ">, md5sum: " << pyinfo[1] << std::endl + << "---- Cylinder::checkMetaData: forcing cache recomputation to ensure consistency" << std::endl; + } + cache_status = 0; + } + } + // End: DeprojType and Python module consistency checks with cache + } + // End: DiskType and Python module consistency checks with cache + } + // End: DiskType attribute check + + return cache_status; +} + +void Cylinder::saveMetaData() +{ + // Only one process must write the dtype metadata + // + if (myid) return; + + std::string dtype = dtstring.at(DTYPE); + + // mtype is already cached as "model" attribute in the HDF5 cache + // file, so we do not need to write it here. We just need to write + // the DiskType and Python metadata (if applicable). + // + try { + // Reopen the cache file for writing the Python metadata + // + HighFive::File file(cachename, HighFive::File::ReadWrite); + + // Write the DiskType attribute (as a string for human readability) + // + file.createAttribute("DiskType", dtype); + + // Write the md5sum for the Python module source file if Python + // disk type is used. This will allow us to check for consistency + // between the Python module used to create the cache and the + // current Python module, and force cache recomputation if they do + // not match. + // + if (DTYPE == DiskType::python) { + try { + std::vector pyinfo = + {pyname, QuickDigest5::fileToHash(pyname)}; + + file.createAttribute("pythonDiskType", pyinfo); + + } catch (const std::runtime_error& e) { + std::cerr << "Cylinder::saveMetaData error: " << e.what() << std::endl; + std::cerr << "Can not write the md5 hash to HDF5" << std::endl; + } + } + + // Deprojection metadata + if (EmpCylSL::mtype == EmpCylSL::EmpModel::Deproject) { + + file.createAttribute("ProjType", dmodel); + + if (PTYPE == DeprojType::python) { + try { + std::vector pyinfo = + {pyproj, QuickDigest5::fileToHash(pyproj + ".py")}; + + file.createAttribute("pythonProjType", pyinfo); + + std::cout << "---- Cylinder::saveMetaData: writing pythonProjType <" << pyproj + ".py" + << "> to cache file <" << cachename << ">" << std::endl; + } catch (const std::runtime_error& e) { + if (myid==0) { + std::cerr << "Cylinder::saveMetaData error: " + << e.what() + << ", can not write the pyinfo and md5 hash to HDF5" + << std::endl; + } + } + } + } + } catch (const HighFive::Exception& err) { + std::cerr << err.what() << std::endl; + std::cerr << "Cylinder::saveMetaData: error writing metadata to cache file <" + << cachename << std::endl; + } +} diff --git a/utils/ICs/check_coefs.cc b/utils/ICs/check_coefs.cc index b8556090b..9d1a9f5f2 100644 --- a/utils/ICs/check_coefs.cc +++ b/utils/ICs/check_coefs.cc @@ -471,35 +471,36 @@ main(int ac, char **av) } } + // Set EmpCylSL mtype. This is the spherical function used to + // generate the EOF basis. If "deproject" is set, this will be + // overriden in EmpCylSL. + // + // Convert mtype string to lower case // std::transform(mtype.begin(), mtype.end(), mtype.begin(), [](unsigned char c){ return std::tolower(c); }); - - // Set EmpCylSL mtype - // - EmpCylSL::mtype = EmpCylSL::Exponential; - if (vm.count("mtype")) { - if (mtype.compare("exponential")==0) - EmpCylSL::mtype = EmpCylSL::Exponential; - else if (mtype.compare("expsphere")==0) - EmpCylSL::mtype = EmpCylSL::ExpSphere; - else if (mtype.compare("gaussian")==0) - EmpCylSL::mtype = EmpCylSL::Gaussian; - else if (mtype.compare("plummer")==0) - EmpCylSL::mtype = EmpCylSL::Plummer; - else if (mtype.compare("power")==0) { - EmpCylSL::mtype = EmpCylSL::Power; - EmpCylSL::PPOW = ppower; - } else { - if (myid==0) std::cout << "No EmpCylSL EmpModel named <" - << mtype << ">, valid types are: " - << "Exponential, ExpSphere, Gaussian, Plummer, Power" << std::endl; - MPI_Finalize(); - return -1; + + // Set EmpCylSL mtype. This is the spherical function used to + // generate the EOF basis. If "deproject" is set, this will be + // overriden in EmpCylSL. + // + auto itm = EmpCylSL::EmpModelMap.find(mtype); + + if (itm == EmpCylSL::EmpModelMap.end()) { + if (myid==0) { + std::cout << "No EmpCylSL EmpModel named <" + << mtype << ">, valid types are: " + << "Exponential, ExpSphere, Gaussian, Plummer, Power, Deproject " + << "(not case sensitive)" << std::endl; } + + MPI_Finalize(); + return 0; } + EmpCylSL::mtype = itm->second; + // Set DiskType // std::transform(disktype.begin(), disktype.end(), disktype.begin(), diff --git a/utils/ICs/check_coefs2.cc b/utils/ICs/check_coefs2.cc index 923a0edb4..225eabe6c 100644 --- a/utils/ICs/check_coefs2.cc +++ b/utils/ICs/check_coefs2.cc @@ -479,30 +479,26 @@ main(int ac, char **av) std::transform(mtype.begin(), mtype.end(), mtype.begin(), [](unsigned char c){ return std::tolower(c); }); - // Set EmpCylSL mtype + // Set EmpCylSL mtype. This is the spherical function used to + // generate the EOF basis. If "deproject" is set, this will be + // overriden in EmpCylSL. // - EmpCylSL::mtype = EmpCylSL::Exponential; - if (vm.count("mtype")) { - if (mtype.compare("exponential")==0) - EmpCylSL::mtype = EmpCylSL::Exponential; - else if (mtype.compare("expsphere")==0) - EmpCylSL::mtype = EmpCylSL::ExpSphere; - else if (mtype.compare("gaussian")==0) - EmpCylSL::mtype = EmpCylSL::Gaussian; - else if (mtype.compare("plummer")==0) { - EmpCylSL::mtype = EmpCylSL::Plummer; - } else if (mtype.compare("power")==0) { - EmpCylSL::mtype = EmpCylSL::Power; - EmpCylSL::PPOW = ppower; - } else { - if (myid==0) std::cout << "No EmpCylSL EmpModel named <" - << mtype << ">, valid types are: " - << "Exponential, ExpSphere, Gaussian, Plummer, Power" << std::endl; - MPI_Finalize(); - return -1; + auto itm = EmpCylSL::EmpModelMap.find(mtype); + + if (itm == EmpCylSL::EmpModelMap.end()) { + if (myid==0) { + std::cout << "No EmpCylSL EmpModel named <" + << mtype << ">, valid types are: " + << "Exponential, ExpSphere, Gaussian, Plummer, Power, Deproject " + << "(not case sensitive)" << std::endl; } + + MPI_Finalize(); + return 0; } + EmpCylSL::mtype = itm->second; + // Set DiskType // std::transform(disktype.begin(), disktype.end(), disktype.begin(), diff --git a/utils/ICs/initial.cc b/utils/ICs/initial.cc index 5f2201aa6..c94929b46 100644 --- a/utils/ICs/initial.cc +++ b/utils/ICs/initial.cc @@ -677,41 +677,41 @@ main(int ac, char **av) SphericalModelTable::linear = 1; } - // Convert mtype string to lower case + // Convert mtype string to lower case + // std::transform(mtype.begin(), mtype.end(), mtype.begin(), [](unsigned char c){ return std::tolower(c); }); - // Convert gentype string to lower case + // Convert gentype string to lower case + // std::transform(gentype.begin(), gentype.end(), gentype.begin(), [](unsigned char c){ return std::tolower(c); }); + // Convert mtype string to lower case + // + std::transform(mtype.begin(), mtype.end(), mtype.begin(), + [](unsigned char c){ return std::tolower(c); }); + // Set EmpCylSL mtype. This is the spherical function used to // generate the EOF basis. If "deproject" is set, this will be // overriden in EmpCylSL. // - EmpCylSL::mtype = EmpCylSL::Exponential; - if (vm.count("mtype")) { - if (mtype.compare("exponential")==0) - EmpCylSL::mtype = EmpCylSL::Exponential; - else if (mtype.compare("expsphere")==0) - EmpCylSL::mtype = EmpCylSL::ExpSphere; - else if (mtype.compare("gaussian")==0) - EmpCylSL::mtype = EmpCylSL::Gaussian; - else if (mtype.compare("plummer")==0) - EmpCylSL::mtype = EmpCylSL::Plummer; - else if (mtype.compare("power")==0) { - EmpCylSL::mtype = EmpCylSL::Power; - EmpCylSL::PPOW = PPower; - } else { - if (myid==0) std::cout << "No EmpCylSL EmpModel named <" - << mtype << ">, valid types are: " - << "Exponential, Gaussian, Plummer, Power " - << "(not case sensitive)" << std::endl; - MPI_Finalize(); - return 0; + auto itm = EmpCylSL::EmpModelMap.find(mtype); + + if (itm == EmpCylSL::EmpModelMap.end()) { + if (myid==0) { + std::cout << "No EmpCylSL EmpModel named <" + << mtype << ">, valid types are: " + << "Exponential, ExpSphere, Gaussian, Plummer, Power, Deproject " + << "(not case sensitive)" << std::endl; } + + MPI_Finalize(); + return 0; } + EmpCylSL::mtype = itm->second; + // Set DiskType. This is the functional form for the disk used to // condition the basis. // diff --git a/utils/Test/CMakeLists.txt b/utils/Test/CMakeLists.txt index 5a8a6fc4b..cc735fe08 100644 --- a/utils/Test/CMakeLists.txt +++ b/utils/Test/CMakeLists.txt @@ -1,5 +1,6 @@ -set(bin_PROGRAMS testBarrier expyaml orthoTest testED) +set(bin_PROGRAMS testBarrier expyaml orthoTest testDeproj + testEmpDeproj testEmp testED testmd5) set(common_LINKLIB OpenMP::OpenMP_CXX MPI::MPI_CXX expui exputil yaml-cpp ${VTK_LIBRARIES}) @@ -10,6 +11,7 @@ endif() set(common_INCLUDE $ + $ $ $ ${CMAKE_BINARY_DIR} ${DEP_INC} @@ -32,7 +34,12 @@ endif() add_executable(testBarrier test_barrier.cc) add_executable(expyaml test_config.cc) add_executable(orthoTest orthoTest.cc Biorth2Ortho.cc) -add_executable(testED testED.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) +add_executable(testmd5 testmd5.cc) +add_executable(testED testED.cc EmpDeproj.cc) foreach(program ${bin_PROGRAMS}) target_link_libraries(${program} ${common_LINKLIB}) @@ -41,4 +48,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..38bd77b11 --- /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 is expected in [xmin(), xmax()]; values outside are extrapolated using the end intervals) + 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..0120021fb --- /dev/null +++ b/utils/Test/EmpDeproj.H @@ -0,0 +1,45 @@ +#pragma once + +#include +#include +#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) + { + assert(R > 0.0); + return densRg.eval(std::log(R)); + } + + double surfaceDensity(double R) + { + assert(R > 0.0); + return surfRg.eval(std::log(R)); + } + + //! Evaluate deprojected mass at a single radius + double mass(double R) + { + assert(R > 0.0); + return massRg.eval(std::log(R)); + } +}; diff --git a/utils/Test/EmpDeproj.cc b/utils/Test/EmpDeproj.cc new file mode 100644 index 000000000..46f28835a --- /dev/null +++ b/utils/Test/EmpDeproj.cc @@ -0,0 +1,148 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifndef M_PI +#define M_PI std::acos(-1.0) +#endif + +#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) +{ + // 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); + + 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 + +#include "Deprojector.H" +#include "cxxopts.H" + +namespace { + const double pi = std::acos(-1.0); +} + +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) / 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*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..524d72e52 --- /dev/null +++ b/utils/Test/testEmp.cc @@ -0,0 +1,304 @@ +#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, Subtraction, 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 + 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); + 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..a502331ba --- /dev/null +++ b/utils/Test/testEmpDeproj.cc @@ -0,0 +1,175 @@ +#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; + + // Define pi in a portable way instead of relying on non-standard M_PI + const double pi = std::acos(-1.0); + + // 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 = [pi](double r)->double + { return 2.0 / std::pow(1.0 + r*r, 2.0) / 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 = [pi](double r)->double + { return exp(-0.5*r*r)/sqrt(2.0*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=*/Ngrid); + + 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 = (Nr > 1) ? static_cast(i) / (Nr - 1) : 0.0; + r_eval.push_back(Rmin + t * (Rmax - Rmin)); + } + 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; +} diff --git a/utils/Test/testmd5.cc b/utils/Test/testmd5.cc new file mode 100644 index 000000000..7490a8d2d --- /dev/null +++ b/utils/Test/testmd5.cc @@ -0,0 +1,52 @@ +// This test verifies that QuickDigest5 correctly computes the MD5 +// hash of a file + +#include +#include + +#include "quickdigest5.hpp" + +#include "exputils.H" // For get_md5sum which uses the system's md5sum + // command for verification + +int main(int argc, char* argv[]) +{ + // Default to example.txt if no argument + std::string filePath = argc > 1 ? argv[1] : "example.txt"; + + // Check if the file exists before trying to hash it + std::filesystem::path p(filePath); + if (!std::filesystem::exists(p)) { + std::cout << "File <" << filePath << "> not found." << std::endl + << "Usage: " << argv[0] << " [file_path]" << std::endl + << "Defaulting to example.txt if it exists." << std::endl; + } + + // One-line method to get the hex digest of a file + std::string hash = QuickDigest5::fileToHash(filePath); + + if (!hash.empty()) { + std::cout << "MD5: " << hash << std::endl; + } else { + std::cerr << "Error: Could not process file." << std::endl; + } + + // System version of md5sum for comparison + try { + if (std::system("which md5sum > /dev/null 2>&1") != 0) { + std::cerr << "Warning: md5sum command not found. Skipping system md5sum comparison." << std::endl; + return 0; + } + std::string systemHash = get_md5sum(filePath); + std::cout << "System md5sum: " << systemHash << std::endl; + if (hash == systemHash) { + std::cout << "Success: hashes match!" << std::endl; + } else { + std::cerr << "Error: hashes do not match!" << std::endl; + } + } catch (const std::exception& e) { + std::cerr << "Error computing system md5sum: " << e.what() << std::endl; + } + + return 0; +}