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 ee224bde3..e5c7cf6ad 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.25) # Needed for CUDA, MPI, and CTest features project( EXP - VERSION "7.9.3" + VERSION "7.10.0" HOMEPAGE_URL https://github.com/EXP-code/EXP LANGUAGES C CXX Fortran) @@ -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/doc/exp.cfg b/doc/exp.cfg index f76157769..c83ce4a1c 100644 --- a/doc/exp.cfg +++ b/doc/exp.cfg @@ -48,7 +48,7 @@ PROJECT_NAME = EXP # could be handy for archiving the generated documentation or if some version # control system is used. -PROJECT_NUMBER = 7.9.3 +PROJECT_NUMBER = 7.10.0 # Using the PROJECT_BRIEF tag one can provide an optional one line description # for a project that appears at the top of each page and should give viewer a diff --git a/doc/exp.cfg.breathe b/doc/exp.cfg.breathe index 740430c50..c8d8a0781 100644 --- a/doc/exp.cfg.breathe +++ b/doc/exp.cfg.breathe @@ -48,7 +48,7 @@ PROJECT_NAME = EXP # could be handy for archiving the generated documentation or if some version # control system is used. -PROJECT_NUMBER = 7.9.3 +PROJECT_NUMBER = 7.10.0 # Using the PROJECT_BRIEF tag one can provide an optional one line description # for a project that appears at the top of each page and should give viewer a diff --git a/expui/BasisFactory.H b/expui/BasisFactory.H index e270fae1e..6f3351497 100644 --- a/expui/BasisFactory.H +++ b/expui/BasisFactory.H @@ -196,15 +196,25 @@ namespace BasisClasses operator()(double x1, double x2, double x3, const Coord ctype=Coord::Spherical); - //! Evaluate fields at a point + //! Evaluate fields at a point in the expansion-origin frame virtual std::vector getFields(double x, double y, double z); + + //! Evaluate fields at a point in the coefficient-origin frame + virtual std::vector getFieldsOrigin(double x, double y, double z); - //! Evaluate fields at a point for all coefficients sets + //! Evaluate fields at a point for all coefficients sets in the + //! expansion-origin frame virtual std::tuple, Eigen::VectorXd> getFieldsCoefs (double x, double y, double z, std::shared_ptr coefs); + + //! Evaluate fields at a point for all coefficients sets in the + //! coefficient-origin frame + virtual std::tuple, + Eigen::VectorXd> getFieldsCoefsOrigin + (double x, double y, double z, std::shared_ptr coefs); - //! Evaluate fields at a point, and provide field lables + //! Evaluate fields at a point, and provide field labels virtual std::tuple, std::vector> evaluate(double x, double y, double z) { return {getFields(x, y, z), getFieldLabels(coordinates)}; } diff --git a/expui/BasisFactory.cc b/expui/BasisFactory.cc index aba544cc3..4afa92617 100644 --- a/expui/BasisFactory.cc +++ b/expui/BasisFactory.cc @@ -229,15 +229,20 @@ namespace BasisClasses } std::vector Basis::getFields(double x, double y, double z) + { return crt_eval(x, y, z); } + + std::vector Basis::getFieldsOrigin(double x, double y, double z) { - return crt_eval(x, y, z); + Eigen::Vector3d p(x, y, z); + p = coefrot * (p - coefctr); + return crt_eval(p(0), p(1), p(2)); } std::tuple, Eigen::VectorXd> Basis::getFieldsCoefs (double x, double y, double z, std::shared_ptr coefs) { - // Python dictonary for return + // Python dictionary for return std::map ret; // Times for the coefficients @@ -249,9 +254,46 @@ namespace BasisClasses // Make the return dictionary of arrays for (int i=0; igetCoefStruct(times[i])); + // The field evaluation auto v = crt_eval(x, y, z); + + // Pack the fields into the dictionary + for (int j=0; j(times.data(), times.size()); + + // Return the dictionary and the time array + return {ret, T}; + } + + std::tuple, Eigen::VectorXd> + Basis::getFieldsCoefsOrigin + (double x, double y, double z, std::shared_ptr coefs) + { + // Python dictionary for return + std::map ret; + + // Times for the coefficients + auto times = coefs->Times(); + + // Initialize the dictionary/map + auto fields = getFieldLabels(coordinates); + for (auto s : fields) ret[s].resize(times.size()); + + // Make the return dictionary of arrays + for (int i=0; igetCoefStruct(times[i])); + + // The field evaluation + auto v = getFieldsOrigin(x, y, z); + // Pack the fields into the dictionary for (int j=0; j 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..9c6d21a30 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 @@ -270,9 +275,19 @@ namespace BasisClasses if (conf["M0_ONLY"]) M0_only = conf["M0_ONLY"].as(); if (conf["pcavar"]) pcavar = conf["pcavar"].as(); if (conf["subsamp"]) sampT = conf["subsamp"].as(); + if (conf["samplesz"]) sampT = conf["samplesz"].as(); sampT = std::max(1, sampT); // Sanity - } + + // Deprecation warning + if (conf["subsamp"]) { + if (myid==0) + std::cout << "---- Spherical: parameter 'subsamp' is deprecated. " + << "It works, but will be removed in version >= 7.11. " + << "Please use 'samplesz' instead." + << std::endl; + } + } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing parameter stanza for <" << name << ">: " @@ -818,6 +833,14 @@ namespace BasisClasses void Spherical::computeAccel(double x, double y, double z, Eigen::Ref acc) { + // Shift and rotate to expansion frame + Eigen::Vector3d pos {x, y, z}; + pos = coefrot * (pos - coefctr); + + x = pos(0); + y = pos(1); + z = pos(2); + // Get polar coordinates double R2 = x*x + y*y; double r2 = R2 + z*z; @@ -920,9 +943,9 @@ namespace BasisClasses double tpoty = (potr - pott*costh/r)*y/r + potp*x/R2; double tpotz = potr*costh + pott*sinth*sinth/r; - // Return force not potential gradient + // Return force not potential gradient and rotate to caller frame // - acc << tpotx, tpoty, tpotz; + acc = coefrot.transpose() * Eigen::Vector3d(tpotx, tpoty, tpotz); } @@ -1180,6 +1203,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 +1246,9 @@ namespace BasisClasses "ashift", "expcond", "ignore", - "deproject", "logr", + "dmodel", + "ppow", "pcavar", "pcaeof", "pcavtk", @@ -1245,6 +1277,7 @@ namespace BasisClasses "coefCompute", "coefMaster", "pyname", + "pyproj", "nint", "totalCovar", "fullCovar" @@ -1391,7 +1424,6 @@ namespace BasisClasses cachename = ".eof_cache_file"; oldcache = false; Ignore = false; - deproject = false; rnum = 200; pnum = 1; @@ -1402,7 +1434,7 @@ namespace BasisClasses EVEN_M = false; cmapR = 1; cmapZ = 1; - mtype = "Exponential"; + mtype = "exponential"; dtype = "exponential"; vflag = 0; @@ -1481,8 +1513,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,13 +1524,15 @@ 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(); + if (conf["samplesz" ]) sampT = conf["samplesz" ].as(); // Sanity sampT = std::max(1, sampT); @@ -1515,14 +1548,23 @@ namespace BasisClasses // Deprecation warning if (conf["eof_file"]) { if (myid==0) - std::cout << "Cylinder: parameter 'eof_file' is deprecated. " - << "and will be removed in a future release. Please " + std::cout << "Cylindrical: parameter 'eof_file' is deprecated. " + << "and will be removed in version >= 7.11. Please " << "use 'cachename' instead." << std::endl; conf["cachename"] = conf["eof_file"]; } + // Deprecation warning + if (conf["subsamp"]) { + if (myid==0) + std::cout << "---- Cylindrical: parameter 'subsamp' is deprecated. " + << "It works, but will be removed in version >= 7.11. " + << "Please use 'samplesz' instead." + << std::endl; + } + } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing 'force' for Component <" @@ -1564,38 +1606,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 +1679,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 +1907,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 +1947,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 +1958,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"; + } + + // 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"); } - else // Default to exponential - model = std::make_shared(1.0, H); + 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 +2029,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(); @@ -1792,8 +2162,12 @@ namespace BasisClasses tdens = sl->accumulated_dens_eval(R, z, phi, tdens0); - double tpotx = tpotR*x/R - tpotp*y/R ; - double tpoty = tpotR*y/R + tpotp*x/R ; + double tpotx = 0.0; + double tpoty = 0.0; + if (R > 0.0) { + tpotx = tpotR*x/R - tpotp*y/R; + tpoty = tpotR*y/R + tpotp*x/R; + } return {tdens0, tdens - tdens0, tdens, @@ -1804,6 +2178,13 @@ namespace BasisClasses void Cylindrical::computeAccel(double x, double y, double z, Eigen::Ref acc) { + // Shift and rotate to expansion frame + Eigen::Vector3d pos {x, y, z}; + pos = coefrot * (pos - coefctr); + x = pos(0); + y = pos(1); + z = pos(2); + double R = sqrt(x*x + y*y); double phi = atan2(y, x); @@ -1813,11 +2194,15 @@ namespace BasisClasses tdens = sl->accumulated_dens_eval(R, z, phi, tdens0); - double tpotx = tpotR*x/R - tpotp*y/R ; - double tpoty = tpotR*y/R + tpotp*x/R ; + double tpotx = 0.0; + double tpoty = 0.0; + if (R > 0.0) { + tpotx = tpotR*x/R - tpotp*y/R; + tpoty = tpotR*y/R + tpotp*x/R; + } - // Apply G to forces on return - acc << tpotx*G, tpoty*G, tpotz*G; + // Apply G and rotate forces on return + acc = coefrot.transpose() * Eigen::Vector3d(tpotx*G, tpoty*G, tpotz*G); } // Evaluate in cylindrical coordinates @@ -2117,8 +2502,18 @@ namespace BasisClasses if (conf["M0_ONLY"]) M0_only = conf["M0_ONLY"].as(); if (conf["pcavar"]) pcavar = conf["pcavar"].as(); if (conf["subsamp"]) sampT = conf["subsamp"].as(); + if (conf["samplesz"]) sampT = conf["samplesz"].as(); sampT = std::max(1, sampT); // Sanity + + if (conf["subsamp"]) { + if (myid==0) + std::cout << "---- FlatDisk: parameter 'subsamp' is deprecated. " + << "It works, but will be removed in version >= 7.11. " + << "Please use 'samplesz' instead." + << std::endl; + } + } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing parameter stanza for <" @@ -2481,6 +2876,13 @@ namespace BasisClasses void FlatDisk::computeAccel(double x, double y, double z, Eigen::Ref acc) { + // Shift and rotate to expansion frame + Eigen::Vector3d pos {x, y, z}; + pos = coefrot * (pos - coefctr); + x = pos(0); + y = pos(1); + z = pos(2); + // Get thread id int tid = omp_get_thread_num(); @@ -2504,7 +2906,7 @@ namespace BasisClasses rpot = -G*totalMass*R/(r*r2 + 10.0*std::numeric_limits::min()); zpot = -G*totalMass*z/(r*r2 + 10.0*std::numeric_limits::min()); - acc << rpot, zpot, ppot; + acc = coefrot.transpose() * Eigen::Vector3d(rpot, ppot, zpot); } // Get the basis fields @@ -2565,10 +2967,14 @@ namespace BasisClasses zpot *= -G; ppot *= -G; - double potx = rpot*x/R - ppot*y/R; - double poty = rpot*y/R + ppot*x/R; + double potx = 0.0; + double poty = 0.0; + if (R > 0.0) { + potx = rpot*x/R - ppot*y/R; + poty = rpot*y/R + ppot*x/R; + } - acc << potx, poty, zpot; + acc = coefrot.transpose() * Eigen::Vector3d(potx, poty, zpot); } @@ -3270,6 +3676,13 @@ namespace BasisClasses void CBDisk::computeAccel(double x, double y, double z, Eigen::Ref acc) { + // Shift and rotate to expansion frame + Eigen::Vector3d pos {x, y, z}; + pos = coefrot * (pos - coefctr); + x = pos(0); + y = pos(1); + z = pos(2); + // Get thread id int tid = omp_get_thread_num(); @@ -3330,10 +3743,14 @@ namespace BasisClasses rpot *= -G; ppot *= -G; - double potx = rpot*x/R - ppot*y/R; - double poty = rpot*y/R + ppot*x/R; + double potx = 0.0; + double poty = 0.0; + if (R > 0.0) { + potx = rpot*x/R - ppot*y/R; + poty = rpot*y/R + ppot*x/R; + } - acc << potx, poty, zpot; + acc = coefrot.transpose() * Eigen::Vector3d(potx, poty, zpot); } @@ -3756,6 +4173,13 @@ namespace BasisClasses void Slab::computeAccel(double x, double y, double z, Eigen::Ref acc) { + // Shift and rotate to expansion frame + Eigen::Vector3d pos {x, y, z}; + pos = coefrot * (pos - coefctr); + x = pos(0); + y = pos(1); + z = pos(2); + // Loop indices // int ix, iy, iz; @@ -3827,8 +4251,9 @@ namespace BasisClasses } } - // Apply G to forces on return - acc << G*accx.real(), G*accy.real(), G*accz.real(); + // Apply G to forces on return and rotate to caller frame + acc = coefrot.transpose() * + Eigen::Vector3d(G*accx.real(), G*accy.real(), G*accz.real()); } @@ -3967,8 +4392,9 @@ namespace BasisClasses "verbose", "check", "method", - "pcavar," + "pcavar", "subsamp", + "samplesz", "nint", "totalCovar", "fullCovar" @@ -4027,20 +4453,30 @@ namespace BasisClasses // Assign values from YAML // try { - if (conf["nminx"]) nminx = conf["nminx" ].as(); - if (conf["nminy"]) nminy = conf["nminy" ].as(); - if (conf["nminz"]) nminz = conf["nminz" ].as(); + if (conf["nminx"]) nminx = conf["nminx" ].as(); + if (conf["nminy"]) nminy = conf["nminy" ].as(); + if (conf["nminz"]) nminz = conf["nminz" ].as(); - if (conf["nmaxx"]) nmaxx = conf["nmaxx" ].as(); - if (conf["nmaxy"]) nmaxy = conf["nmaxy" ].as(); - if (conf["nmaxz"]) nmaxz = conf["nmaxz" ].as(); + if (conf["nmaxx"]) nmaxx = conf["nmaxx" ].as(); + if (conf["nmaxy"]) nmaxy = conf["nmaxy" ].as(); + if (conf["nmaxz"]) nmaxz = conf["nmaxz" ].as(); - if (conf["knots"]) knots = conf["knots" ].as(); + if (conf["knots"]) knots = conf["knots" ].as(); - if (conf["check"]) check = conf["check" ].as(); + if (conf["check"]) check = conf["check" ].as(); + + if (conf["pcavar"]) pcavar = conf["pcavar" ].as(); + if (conf["subsamp"]) sampT = conf["subsamp" ].as(); + if (conf["samplesz"]) sampT = conf["samplesz"].as(); + + if (conf["subsamp"]) { + if (myid==0) + std::cout << "---- Cube: parameter 'subsamp' is deprecated. " + << "It works, but will be removed in version >= 7.11. " + << "Please use 'samplesz' instead." + << std::endl; + } - if (conf["pcavar"]) pcavar = conf["pcavar" ].as(); - if (conf["subsamp"]) sampT = conf["subsamp"].as(); } catch (YAML::Exception & error) { if (myid==0) std::cout << "Error parsing parameter stanza for <" @@ -4328,17 +4764,22 @@ namespace BasisClasses void Cube::computeAccel(double x, double y, double z, Eigen::Ref acc) { + // Shift and rotate to expansion frame + Eigen::Vector3d pos {x, y, z}; + pos = coefrot * (pos - coefctr); + x = pos(0); + y = pos(1); + z = pos(2); + // Get thread id int tid = omp_get_thread_num(); - // Position vector - Eigen::Vector3d pos {x, y, z}; - // Get the basis fields auto frc = ortho->get_force(expcoef, pos); - // Apply G to forces on return - acc << -G*frc(0).real(), -G*frc(1).real(), -G*frc(2).real(); + // Apply G to forces on return and rotate to caller frame + acc = coefrot.transpose() * + Eigen::Vector3d(-G*frc(0).real(), -G*frc(1).real(), -G*frc(2).real()); } std::vector Cube::cyl_eval(double R, double z, double phi) diff --git a/expui/CMakeLists.txt b/expui/CMakeLists.txt index f25d0d7b7..15d7e4509 100644 --- a/expui/CMakeLists.txt +++ b/expui/CMakeLists.txt @@ -97,11 +97,15 @@ target_sources(expui PUBLIC FILE_SET HEADERS ${CMAKE_CURRENT_SOURCE_DIR}/../include/DiskWithHalo.H ${CMAKE_CURRENT_SOURCE_DIR}/../include/EmpCyl2d.H ${CMAKE_CURRENT_SOURCE_DIR}/../include/EXPmath.H + # ${CMAKE_CURRENT_SOURCE_DIR}/../include/EXPversion.H ${CMAKE_CURRENT_SOURCE_DIR}/../include/libvars.H ${CMAKE_CURRENT_SOURCE_DIR}/../include/Timer.H ${CMAKE_CURRENT_SOURCE_DIR}/../include/coef.H ${CMAKE_CURRENT_SOURCE_DIR}/../include/Covariance.H ${CMAKE_CURRENT_SOURCE_DIR}/../include/ExpDeproj.H + ${CMAKE_CURRENT_SOURCE_DIR}/../include/cudaUtil.cuH + ${CMAKE_CURRENT_SOURCE_DIR}/../include/cudaParticle.cuH + ${CMAKE_CURRENT_SOURCE_DIR}/../include/cudaMappingConstants.cuH ) install(TARGETS expui FILE_SET HEADERS DESTINATION include/EXP) 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 4357d0cf7..2bc068a51 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/BiorthCube.H b/include/BiorthCube.H index a89ebb87d..808ed7f49 100644 --- a/include/BiorthCube.H +++ b/include/BiorthCube.H @@ -17,8 +17,9 @@ #include #if HAVE_LIBCUDA==1 -#include -#include +#include "cudaUtil.cuH" +#include "cudaParticle.cuH" +#include "cudaMappingConstants.cuH" #endif // For reading and writing cache file @@ -28,11 +29,6 @@ #include #include -#if HAVE_LIBCUDA==1 -#include -#include -#endif - #include "EmpCyl2d.H" //!! BiorthCube grid class diff --git a/include/BiorthCyl.H b/include/BiorthCyl.H index 687433586..226120441 100644 --- a/include/BiorthCyl.H +++ b/include/BiorthCyl.H @@ -17,8 +17,9 @@ #include #if HAVE_LIBCUDA==1 -#include -#include +#include "cudaUtil.cuH" +#include "cudaParticle.cuH" +#include "cudaMappingConstants.cuH" #endif // For reading and writing cache file @@ -28,11 +29,6 @@ #include #include -#if HAVE_LIBCUDA==1 -#include -#include -#endif - #include "EmpCyl2d.H" //!! BiorthCyl grid class 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..dff248f2b 100644 --- a/include/EmpCylSL.H +++ b/include/EmpCylSL.H @@ -24,8 +24,8 @@ #include "coef.H" #if HAVE_LIBCUDA==1 -#include -#include +#include "cudaParticle.cuH" +#include "cudaMappingConstants.cuH" #endif #include "libvars.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/SLGridMP2.H b/include/SLGridMP2.H index da209403c..dd46e28cf 100644 --- a/include/SLGridMP2.H +++ b/include/SLGridMP2.H @@ -19,8 +19,8 @@ using namespace __EXP__; #if HAVE_LIBCUDA==1 -#include -#include +#include "cudaUtil.cuH" +#include "cudaMappingConstants.cuH" #endif diff --git a/include/cudaMappingConstants.cuH b/include/cudaMappingConstants.cuH index 067d299b6..3e845436d 100644 --- a/include/cudaMappingConstants.cuH +++ b/include/cudaMappingConstants.cuH @@ -1,9 +1,7 @@ -// -*- C++ -*- - #ifndef CUDA_MAPPING_CONSTANTS_H #define CUDA_MAPPING_CONSTANTS_H -#include +#include "cudaUtil.cuH" struct cudaMappingConstants { @@ -23,3 +21,5 @@ struct cudaMappingConstants }; #endif + +// -*- C++ -*- diff --git a/include/cudaParticle.cuH b/include/cudaParticle.cuH index 5147b01dd..ccc2102b1 100644 --- a/include/cudaParticle.cuH +++ b/include/cudaParticle.cuH @@ -1,5 +1,3 @@ -// -*- C++ -*- - #ifndef PARTICLE_CUH #define PARTICLE_CUH @@ -16,11 +14,9 @@ #include #include -#include - -#include - -#include +#include "config_exp.h" +#include "cudaUtil.cuH" +#include "Particle.H" //! Simplified particle structure for use in CUDA kernel code struct cudaParticle @@ -125,3 +121,5 @@ struct cuPartToChange }; #endif + +// -*- C++ -*- 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/pyEXP/BasisWrappers.cc b/pyEXP/BasisWrappers.cc index 1429c6249..cd4f650cb 100644 --- a/pyEXP/BasisWrappers.cc +++ b/pyEXP/BasisWrappers.cc @@ -137,9 +137,9 @@ void BasisFactoryClasses(py::module &m) support for computing the coefficieint covariance from subsamples of particles. This is implemented by the enableCoefCovariance() method for each of these supported bases. The force configuration must - contain the parameters 'pcavar' (boolean) and 'subsamp' (integer) keys. + contain the parameters 'pcavar' (boolean) and 'samplesz' (integer) keys. The 'pcavar' parameter turns on the covariance computation. The - 'subsamp' parameter sets the number of partions or subsamples for each + 'samplesz' parameter sets the number of partitions or subsamples for each coefficient creation. There are two additional control parameters that may be optionally specified with the enableCoefCovariance() call. The 'total' parameters enables computing the total covariance matrices only @@ -176,11 +176,11 @@ void BasisFactoryClasses(py::module &m) counts, the sample masses, the coefficient means and the coefficient covariance matrices. The first two elements are vectors of length equal to the number of subsamples. The third element is a 4D array - with dimensions (subsamp, Nlm, nmax) where subsamp is the number of + with dimensions (samplesz, Nlm, nmax) where samplesz is the number of subsamples, Nlm, is the number of harmonics (e.g. l, m values or m values for the spherical and cylindrical bases), and nmax is the number of coefficients. The fourth element is a 4D array with dimensions - (subsamp, Nlm, nmax, namx) containing the covariance matrices for each + (samplesz, Nlm, nmax, namx) containing the covariance matrices for each subsample in the last two dimensions. Typical usage might be: @@ -562,7 +562,11 @@ void BasisFactoryClasses(py::module &m) using Spherical::Spherical; std::vector getFields(double x, double y, double z) override { - PYBIND11_OVERRIDE(std::vector, Spherical, getFields, x, y, z); + py::function py_override = py::get_override(static_cast(this), "getFields"); + if (py_override) { + return py::cast>(py_override(x, y, z)); + } + return Spherical::getFields(x, y, z); } void accumulate(double x, double y, double z, double mass, unsigned long int indx) override { @@ -1522,7 +1526,7 @@ void BasisFactoryClasses(py::module &m) potential evaluations are separated into full, axisymmetric and non-axisymmetric contributions. - You can get the fields labels by using the __call__ method of the + You can get the field labels by using the __call__ method of the basis object. This is equivalent to a tuple of the getFields() output with a list of field labels. @@ -1541,13 +1545,38 @@ void BasisFactoryClasses(py::module &m) See also -------- + getFieldsOrigin: get fields in coefficient-origin frame getFieldsCoefs : get fields for each coefficient set __call__ : same as getFields() but provides field labels in a tuple )", py::arg("x"), py::arg("y"), py::arg("z")) + .def("getFieldsOrigin", &BasisClasses::BiorthBasis::getFieldsOrigin, + R"( + Return the field evaluations for a given Cartesian position in + the frame defined by the current coefficients. + + Parameters + ---------- + x : float + x-axis position + y : float + y-axis position + z : float + z-axis position + + Returns + ------- + fields: numpy.ndarray + + See also + -------- + getFields : get fields in expansion-origin frame + getFieldsCoefs : get fields for each coefficient set + )", + py::arg("x"), py::arg("y"), py::arg("z")) .def("getAccel", py::overload_cast(&BasisClasses::BiorthBasis::getAccel), R"( - Return the acceleration for a given Cartesian position + Return the acceleration for a given Cartesian position in the frame defined by the coefficients. Parameters ---------- @@ -1571,7 +1600,7 @@ void BasisFactoryClasses(py::module &m) py::arg("x"), py::arg("y"), py::arg("z")) .def("getAccel", py::overload_cast, Eigen::Ref, Eigen::Ref>(&BasisClasses::BiorthBasis::getAccel), R"( - Return the acceleration for a given Cartesian position + Return the acceleration for a given Cartesian position in the frame defined by the coefficients. Parameters ---------- @@ -1595,7 +1624,7 @@ void BasisFactoryClasses(py::module &m) py::arg("x"), py::arg("y"), py::arg("z")) .def("getAccel", py::overload_cast>(&BasisClasses::BiorthBasis::getAccel), R"( - Return the acceleration for an array of Cartesian positions + Return the acceleration for an array of Cartesian positions in the frame defined by the coefficients. Parameters ---------- @@ -1615,7 +1644,7 @@ void BasisFactoryClasses(py::module &m) py::arg("pos")) .def("getAccelArray", py::overload_cast, Eigen::Ref, Eigen::Ref>(&BasisClasses::BiorthBasis::getAccel), R"( - Return the acceleration for a given Cartesian position + Return the acceleration for a given Cartesian position in the frame defined by the coefficients. Parameters ---------- @@ -1665,8 +1694,36 @@ void BasisFactoryClasses(py::module &m) See also -------- - getFields : get fields for the currently assigned coefficients - __call__ : same getFields() but provides field labels in a tuple + getFields : get fields for the currently assigned coefficients + getFieldsOrigin : get fields in coefficient-origin frame + __call__ : same getFields() but provides field labels in a tuple + )", + py::arg("x"), py::arg("y"), py::arg("z"), py::arg("coefs")) + .def("getFieldsCoefsOrigin", &BasisClasses::BiorthBasis::getFieldsCoefsOrigin, + R"( + Return the field evaluations for a given Cartesian position + for every frame in a coefficient set in the frame defined by + each coefficient structure. + + Parameters + ---------- + x : float + x-axis position + y : float + y-axis position + z : float + z-axis position + coefs: CoefClasses::Coefs + the coefficient set + + Returns + ------- + tuple of a dictionary of fields of array values, and an + array of evaluation times + + See also + -------- + getFieldsCoefs : get fields in expansion-origin frame )", py::arg("x"), py::arg("y"), py::arg("z"), py::arg("coefs")) .def("setFieldType", &BasisClasses::BiorthBasis::setFieldType, 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..b8bd5033a 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 + ".py"); + } 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 + ".py")}; + + 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/src/global_key_set.H b/src/global_key_set.H index cfa4a9ef0..c31745af5 100644 --- a/src/global_key_set.H +++ b/src/global_key_set.H @@ -8,6 +8,7 @@ global_valid_keys = { "dbthresh", "time", "dtime", + "maxMindt", "PFbufsz", "NICE", "VERBOSE", 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 +#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 implementation 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..e925c6c9d --- /dev/null +++ b/utils/Test/testEmp.cc @@ -0,0 +1,306 @@ +#include +#include +#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..870f8f6e9 --- /dev/null +++ b/utils/Test/testEmpDeproj.cc @@ -0,0 +1,176 @@ +#include +#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..3518f9eea --- /dev/null +++ b/utils/Test/testmd5.cc @@ -0,0 +1,56 @@ +// This test verifies that QuickDigest5 correctly computes the MD5 +// hash of a file + +#include +#include +#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::cerr << "File <" << filePath << "> not found." << std::endl + << "Usage: " << argv[0] << " [file_path]" << std::endl; + return 1; + } + + // 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; + return 1; + } + + // 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; + return 1; + } + } catch (const std::exception& e) { + std::cerr << "Error computing system md5sum: " << e.what() << std::endl; + return 1; + } + + return 0; +}