diff --git a/.github/workflows/test_examples.yml b/.github/workflows/test_examples.yml index 62ae5a09a..1fe212a99 100644 --- a/.github/workflows/test_examples.yml +++ b/.github/workflows/test_examples.yml @@ -32,7 +32,15 @@ jobs: run: | mkdir build cd build - cmake .. -DCMAKE_INSTALL_PREFIX=$HOME/.local -DENABLE_PYTHON=True -DPython_EXECUTABLE=${{ matrix.config.py }} -DENABLE_TESTING=Off -DENABLE_SWIG_BUILTIN=${{ matrix.config.swig_builtin }} -DSIMD_EXTENSIONS=native -DPython_INSTALL_PACKAGE_DIR=/home/runner/.local/ + cmake .. -DCMAKE_INSTALL_PREFIX=$HOME/.local \ + -DENABLE_PYTHON=True \ + -DPython_EXECUTABLE=${{ matrix.config.py }} \ + -DENABLE_TESTING=Off \ + -DENABLE_SWIG_BUILTIN=${{ matrix.config.swig_builtin }} \ + -DSIMD_EXTENSIONS=native \ + -DPython_INSTALL_PACKAGE_DIR=/home/runner/.local/ \ + -DCMAKE_CXX_STANDARD=17 \ + -DCMAKE_CXX_FLAGS="-std=c++17 -I${{ github.workspace }}/libs/nanoflann/include -Wall -Wextra" - name: Build CRPropa run: | cd build diff --git a/.github/workflows/testing_OSX.yml b/.github/workflows/testing_OSX.yml index c86e2398e..c6f9fc834 100644 --- a/.github/workflows/testing_OSX.yml +++ b/.github/workflows/testing_OSX.yml @@ -54,7 +54,7 @@ jobs: CRPROPA_DATA_PATH: "/Users/runner/work/CRPropa3/CRPropa3/build/data" run: | cd build - make test + make test ARGS="--output-on-failure" - name: Archive test results if: always() uses: actions/upload-artifact@v4 diff --git a/.github/workflows/testing_ubuntu2 b/.github/workflows/testing_ubuntu2 new file mode 100644 index 000000000..e69de29bb diff --git a/.github/workflows/testing_ubuntu22.yml b/.github/workflows/testing_ubuntu22.yml index 49b60a65b..e79bcddb8 100644 --- a/.github/workflows/testing_ubuntu22.yml +++ b/.github/workflows/testing_ubuntu22.yml @@ -33,7 +33,13 @@ jobs: run: | mkdir build cd build - cmake .. -DENABLE_PYTHON=True -DENABLE_TESTING=On -DPython_EXECUTABLE=${{ matrix.config.py }} -DENABLE_SWIG_BUILTIN=${{ matrix.config.swig_builtin }} -DSIMD_EXTENSIONS=native + cmake .. -DENABLE_PYTHON=True \ + -DENABLE_TESTING=On \ + -DPython_EXECUTABLE=${{ matrix.config.py }} \ + -DENABLE_SWIG_BUILTIN=${{ matrix.config.swig_builtin }} \ + -DSIMD_EXTENSIONS=native \ + -DCMAKE_CXX_STANDARD=17 \ + -DCMAKE_CXX_FLAGS="-std=c++17 -Wall -Wextra -I${{ github.workspace }}/libs/nanoflann/include" - name: Build CRPropa run: | cd build @@ -41,7 +47,7 @@ jobs: - name: Run tests run: | cd build - make test + make test ARGS="--output-on-failure" - name: Archive test results if: always() uses: actions/upload-artifact@v4 diff --git a/.gitignore b/.gitignore index 75599931a..d38055167 100644 --- a/.gitignore +++ b/.gitignore @@ -26,3 +26,4 @@ cmake/CMakeFiles/ *.dat *.ipynb_checkpoints *.idea +.DS_Store diff --git a/CHANGELOG.md b/CHANGELOG.md index 344caca60..e75a68683 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,15 @@ +## PR #535 position dependent photon field for EM interactions (optimised) + +### fixes after JanNiklasB comments: +* The amount of includes has been optimised +* The splitFilename has been added to Common.h +* Comments have been added on functions +* The flags for the PhotonBackground are more explicit and efficient (the redundant if sentences have been removed) +* The CMakeFiles.txt has been implemented following the comments + +### in the github PR the reasons why in the new design the InteractionRates class have been implemented for both homogeneous and inhomogeneous case have been elaborated. + + ## CRPropa vNext ### Bug fixes: diff --git a/CMakeLists.txt b/CMakeLists.txt index d5424f4db..0bd664743 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,7 +3,16 @@ cmake_minimum_required(VERSION 3.14) project(CRPropa Fortran C CXX) set(CRPROPA_RELEASE_VERSION 3.2.1+) # Update for new release -set(CMAKE_CXX_STANDARD 11) +# Detect platform and set C++ standard +if(APPLE) + message(STATUS "Building on macOS, using C++11") + set(CMAKE_CXX_STANDARD 11) + set(CMAKE_CXX_STANDARD_REQUIRED ON) +else() + message(STATUS "Building on Linux/Ubuntu, using C++17") + set(CMAKE_CXX_STANDARD 17) + set(CMAKE_CXX_STANDARD_REQUIRED ON) +endif() set(CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH}) set(CRPROPA_EXTRA_SOURCES) @@ -57,7 +66,6 @@ if(NOT CMAKE_BUILD_TYPE) endif(NOT CMAKE_BUILD_TYPE) message(STATUS "Build Type: ${CMAKE_BUILD_TYPE}") - # ---------------------------------------------------------------------------- # Version info from Git # ---------------------------------------------------------------------------- @@ -146,6 +154,9 @@ add_subdirectory(libs/kiss) list(APPEND CRPROPA_EXTRA_LIBRARIES kiss) list(APPEND CRPROPA_EXTRA_INCLUDES libs/kiss/include) +# nanoflann (provided) +list(APPEND CRPROPA_EXTRA_INCLUDES libs/nanoflann/include) + # HepID (provided) add_subdirectory(libs/HepPID) list(APPEND CRPROPA_EXTRA_LIBRARIES HepPID) @@ -335,7 +346,7 @@ CRPropa should compile, but will likely not work properly! Please install data f endif() # ---------------------------------------------------------------------------- -# Library and Binary +# Library and Binary # ---------------------------------------------------------------------------- file(GLOB_RECURSE CRPROPA_INCLUDES RELATIVE ${CMAKE_SOURCE_DIR} include/*.h) include_directories(include ${CRPROPA_EXTRA_INCLUDES}) @@ -355,6 +366,7 @@ add_library(crpropa SHARED src/ParticleMass.cpp src/ParticleState.cpp src/PhotonBackground.cpp + src/InteractionRates.cpp src/ProgressBar.cpp src/Random.cpp src/Source.cpp @@ -441,7 +453,10 @@ else() MESSAGE(STATUS "Build of documentation disabeled. Enable with BUILD_DOC=On") endif(BUILD_DOC) - +# Only add stdc++fs for older GCC (needed for Ubuntu 20) +if(CMAKE_CXX_COMPILER_ID STREQUAL "GNU" AND CMAKE_CXX_COMPILER_VERSION VERSION_LESS 10) + target_link_libraries(crpropa stdc++fs) +endif() # ---------------------------------------------------------------------------- # Python @@ -546,7 +561,6 @@ if(ENABLE_PYTHON AND Python_FOUND) endif(ENABLE_PYTHON AND Python_FOUND) - # ---------------------------------------------------------------------------- # Install # ---------------------------------------------------------------------------- diff --git a/doc/pages/example_notebooks/photon_propagation/galactic_propagation_gammarays_ISRF.ipynb b/doc/pages/example_notebooks/photon_propagation/galactic_propagation_gammarays_ISRF.ipynb new file mode 100644 index 000000000..91b303035 --- /dev/null +++ b/doc/pages/example_notebooks/photon_propagation/galactic_propagation_gammarays_ISRF.ipynb @@ -0,0 +1,268 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Galactic propagation of highly-energetic gamma rays\n", + "\n", + "Through this example, one can test the propagation of galactic gamma rays from point-like sources, also including the electromagnetic processes with the 3D spatial models of the interstellar radiation field (ISRF) and the cosmic microwave background (CMB). More details on this new code implementation at https://arxiv.org/abs/2507.11475." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } + }, + "outputs": [], + "source": [ + "from crpropa import *\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "# gamma-ray source cartesian position (kpc).\n", + "x = -7.\n", + "y = 6.8\n", + "z = 0.\n", + "\n", + "# source distance from the Earth, placed at (-8.5, 0, 0) kpc in the JF12 galactic magnetic field model. \n", + "D = np.sqrt((x - (-8.5))**2. + y * y + z * z)\n", + "\n", + "# aperture angle of the source (rad).\n", + "Acon = 0.1\n", + "\n", + "# radius of the spherical observer placed at the Earth position (kpc).\n", + "Orad = np.tan(Scon) * D \n", + "\n", + "# number of gamma rays injected by the source.\n", + "nEvents = 100000\n", + "\n", + "# gamma-ray source injection spectrum (eV).\n", + "Emin = 1e12 \n", + "Emax = 1e15 \n", + "specIndex = -1\n", + "\n", + "# particles track\n", + "photons = True\n", + "electrons = True\n", + "thinning = 1.\n", + "minEnergy = Emin #eV\n", + "\n", + "# simulation steps. \n", + "minStep = 1e-3 #kpc\n", + "maxStep = 1e-2 #kpc\n", + "tol = 1e-4\n", + "\n", + "# magnetic field setup.\n", + "B = JF12Field()\n", + "randomSeed = 691342\n", + "B.randomStriated(randomSeed)\n", + "B.randomTurbulent(randomSeed)\n", + "\n", + "# output coordinates\n", + "path = \"...\"\n", + "outputName = 'galactic_propagation_3Disrf.txt'\n", + "outputName_inj = 'galactic_propagation_3Disrf_inj.txt'\n", + "\n", + "# photon backgrounds:\n", + "# homogeneuos, \n", + "cmb = CMB()\n", + "ebl = IRB_Gilmore12()\n", + "\n", + "# position dependent. (Sphere to load only the ISRF nodes close to observer and the source.)\n", + "isrf = ISRF_Freudenreich98(Sphere(Vector3d(x, y, z) * kpc, D * kpc + 1 * kpc))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot the background densities\n", + "\n", + "A plot to check the ISRF densities at different positions in the Galaxy (cartesian coordinates). The code takes the ISRF density as the one of the closest node to each position. These are compared to the CMB and extragalactic background light (EBL) homogeneous fields. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# positions where to compute the ISRF densities. \n", + "pos = Vector3d(-9, -1.5, 0) * kpc\n", + "pos1 = Vector3d(-9, 1.5, 0) * kpc\n", + "pos2 = Vector3d(-9, 2, -2) * kpc\n", + "pos3 = Vector3d(6, 1., 1.) * kpc\n", + "\n", + "eps = np.logspace(-4, 5, 1000)\n", + "\n", + "y = np.zeros(len(eps))\n", + "y1 = np.zeros(len(eps))\n", + "y2 = np.zeros(len(eps))\n", + "y3 = np.zeros(len(eps))\n", + "\n", + "ycmb = np.zeros(len(eps))\n", + "yebl = np.zeros(len(eps))\n", + "\n", + "\n", + "for i in range(len(eps)):\n", + " y[i] = isrf.getPhotonDensity(eps[i] * eV, 0., pos)\n", + " y1[i] = isrf.getPhotonDensity(eps[i] * eV, 0., pos1)\n", + " y2[i] = isrf.getPhotonDensity(eps[i] * eV, 0., pos2)\n", + " y3[i] = isrf.getPhotonDensity(eps[i] * eV, 0., pos3)\n", + " ycmb[i] = cmb.getPhotonDensity(eps[i] * eV, 0., pos)\n", + " yebl[i] = ebl.getPhotonDensity(eps[i] * eV, 0., pos)\n", + " \n", + "plt.plot(eps, y, label = \"ISRF (F98) at (-9, -1.5, 0) kpc\", c='darkgrey')\n", + "plt.plot(eps, y1, label = \"at (-9, 1.5, 0) kpc\", c='y')\n", + "plt.plot(eps, y2, label = \"at (-9, 2, -2) kpc\", c='orange')\n", + "plt.plot(eps, y3, label = \"at (6, 1, 1) kpc\", c='darksalmon')\n", + "plt.plot(eps, ycmb, label = \"CMB\", c='b', linestyle='-.')\n", + "plt.plot(eps, yebl, label = \"EBL\", c='g', linestyle=':')\n", + "plt.loglog()\n", + "plt.ylim(1e0,1e10)\n", + "plt.xlim(2e-4,1e2)\n", + "plt.legend(ncol = 3, loc= \"upper center\", bbox_to_anchor=(0.5, 1.12))\n", + "plt.xlabel(\"$\\epsilon$ [eV]\")\n", + "plt.ylabel(r\"n [$m^{-3}$]\")\n", + "plt.savefig(path + \"bkgDensities.png\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Defining the galactic gamma-ray source and the observer\n", + "\n", + "The source points towards the spherical osberver at the Earth. An second tiny spherical observer is placed at the source position to get the injection spectrum. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# versor pointing to the observer.\n", + "vx = -(x + 8.5)\n", + "vy = -y\n", + "vz = -z\n", + "v = Vector3d(vx, vy, vz)\n", + "v /= v.getR() \n", + "\n", + "source = Source()\n", + "source.add(SourcePosition(Vector3d(x, y, z) * kpc))\n", + "source.add(SourceEmissionCone(v, Acon))\n", + "source.add(SourcePowerLawSpectrum(Emin * eV, Emax * eV, specIndex))\n", + "source.add(SourceParticleType(22))\n", + "\n", + "# observer surrounding the source.\n", + "obs1 = Observer()\n", + "obs1.add(ObserverSurface(Sphere(Vector3d(x, y, z) * kpc, 1e-6 * kpc)))\n", + "output1 = TextOutput(path + outputName_inj, Output.Everything) \n", + "output1.enable(Output.WeightColumn)\n", + "output1.enable(Output.CandidateTagColumn)\n", + "output1.setEnergyScale(EeV)\n", + "obs1.onDetection(output1)\n", + "obs1.setDeactivateOnDetection(False)\n", + "\n", + "\n", + "# observer at the Earth. \n", + "obs = Observer()\n", + "obs.add(ObserverSurface(Sphere(Vector3d(-8.5, 0., 0.) * kpc, Orad * kpc)))\n", + "output = TextOutput(path + outputName, Output.Everything)\n", + "output.enable(Output.WeightColumn)\n", + "output.enable(Output.CandidateTagColumn)\n", + "output.setEnergyScale(EeV)\n", + "obs.onDetection(output)\n", + "obs.setDeactivateOnDetection(True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Running the simulation\n", + "\n", + "The full development of the electromagnetic cascade is simulated, considering also the deflection of the electrons due to the galactic magnetic field. The processes at play are pair production and inverse Compton scatterings, both on the CMB and the ISRF.\n", + "\n", + "The ISRF interaction tables are restricted to the ones in the region of interest, i.e. the space between the observer and the source. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# module setup\n", + "sim = ModuleList()\n", + "sim.add(PropagationBP(B, tol, minStep * kpc, maxStep * kpc))\n", + "\n", + "sim.add(EMPairProduction(cmb, electrons, thinning))\n", + "sim.add(EMPairProduction(isrf, electrons, thinning, Sphere(Vector3d(x, y, z) * kpc, D * kpc + 1 * kpc)))\n", + "\n", + "sim.add(EMInverseComptonScattering(cmb, electrons, thinning))\n", + "sim.add(EMInverseComptonScattering(cmb, electrons, thinning, Sphere(Vector3d(x, y, z) * kpc, D * kpc + 1 * kpc)))\n", + "\n", + "sim.add(MinimumEnergy(minEnergy * eV))\n", + "sim.add(obs1)\n", + "sim.add(obs)\n", + "sim.add(MaximumTrajectoryLength((D + 5) * kpc))\n", + "\n", + "# run simulation\n", + "sim.setShowProgress(True)\n", + "sim.run(source, nEvents)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Resulting energy spectra (LHAASO one?)\n", + "\n", + "Once performed the simulations, the energy spectra can be plotted (after being re-weighted to a certain power-law index). \n", + " " + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "crp_docu", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.2" + }, + "vscode": { + "interpreter": { + "hash": "c416687c884a42c367c2f4b19e8bea2627679ca3202fbf20d972b7cd00ee0b77" + } + }, + "widgets": { + "application/vnd.jupyter.widget-state+json": { + "state": {}, + "version_major": 2, + "version_minor": 0 + } + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/include/crpropa/Common.h b/include/crpropa/Common.h index ac7bc3cd3..ef8636f00 100644 --- a/include/crpropa/Common.h +++ b/include/crpropa/Common.h @@ -49,8 +49,10 @@ double interpolateEquidistant(double x, double lo, double hi, // Find index of value in a sorted vector X that is closest to x size_t closestIndex(double x, const std::vector &X); -/** @}*/ +// Takes the filename from a full data path. Used in EM* modules. +std::string splitFilename(const std::string str); +/** @}*/ // pow implementation as template for integer exponents pow_integer<2>(x) // evaluates to x*x diff --git a/include/crpropa/Geometry.h b/include/crpropa/Geometry.h index 47d800610..aaae70d2b 100644 --- a/include/crpropa/Geometry.h +++ b/include/crpropa/Geometry.h @@ -32,6 +32,7 @@ class Surface : public Referenced { @param point vector corresponding to the point to which compute the normal vector */ virtual Vector3d normal(const Vector3d& point) const = 0; + virtual bool isInside(const Vector3d& point) const = 0; virtual std::string getDescription() const {return "Surface without description.";}; }; @@ -48,6 +49,7 @@ class Plane: public Surface { Plane(const Vector3d& x0, const Vector3d& n); virtual double distance(const Vector3d &x) const; virtual Vector3d normal(const Vector3d& point) const; + virtual bool isInside(const Vector3d& point) const; virtual std::string getDescription() const; }; @@ -64,6 +66,7 @@ class Sphere: public Surface { Sphere(const Vector3d& center, double radius); virtual double distance(const Vector3d &point) const; virtual Vector3d normal(const Vector3d& point) const; + virtual bool isInside(const Vector3d& point) const; virtual std::string getDescription() const; }; @@ -79,6 +82,7 @@ class ParaxialBox: public Surface { ParaxialBox(const Vector3d& corner, const Vector3d& size); virtual double distance(const Vector3d &point) const; virtual Vector3d normal(const Vector3d& point) const; + virtual bool isInside(const Vector3d& point) const; virtual std::string getDescription() const; }; diff --git a/include/crpropa/InteractionRates.h b/include/crpropa/InteractionRates.h new file mode 100644 index 000000000..4118dd569 --- /dev/null +++ b/include/crpropa/InteractionRates.h @@ -0,0 +1,203 @@ +#ifndef CRPROPA_INTERACTIONRATES_H +#define CRPROPA_INTERACTIONRATES_H + +#include "crpropa/Common.h" +#include "crpropa/Referenced.h" +#include "crpropa/Vector3.h" +#include "crpropa/Geometry.h" + +#include + +#include +#include +#include + +namespace crpropa { + +struct PointCloud { + std::vector points; + std::vector ids; + + inline size_t kdtree_get_point_count() const { return points.size(); } + + inline double kdtree_get_pt(const size_t idx, const size_t dim) const { + if (dim == 0) return points[idx].x; + if (dim == 1) return points[idx].y; + return points[idx].z; + } + + // optional bounding-box computation (required by nanoflann) + template + bool kdtree_get_bbox(BBOX& /*bb*/) const { + return false; // no bounding box optimization + } + +}; + +using KDTree = nanoflann::KDTreeSingleIndexAdaptor< + nanoflann::L2_Simple_Adaptor, + PointCloud, + 3 +>; + +/** + * \addtogroup InteractionRates + * @{ + */ + +/** + @class Interaction Rates + @brief Abstract base class for photon fields interaction rates. + */ +class InteractionRates: public Referenced { +public: + virtual double getProcessRate(const double E, const Vector3d &position) const = 0; + virtual void loadPerformInteractionTabs(const Vector3d &position, std::vector &tabE, std::vector &tabs, std::vector> &tabCDF) const = 0; + + std::string getRatesName() const { + return this->ratesName; + } + + bool hasPositionDependence() const { + return this->isPositionDependent; + } + + void setRatesName(std::string ratesName) { + this->ratesName = ratesName; + } + + virtual void initRate(std::string path) = 0; + virtual void initCumulativeRate(std::string path) = 0; + +protected: + + std::string ratesName = "AbstractInteractionRates"; + bool isPositionDependent = false; + +}; + +/** + @class InteractionRateHomogeneous + @brief Interaction rates decorator for tabulated homogeneous photon fields. + */ +class InteractionRatesHomogeneous: public InteractionRates { +public: + /** Constructor of InteractionRatesHomogeneous + * @param RateFile Path to the file containing the interaction rate data + * @param CumulativeRateFile Path to the file containing the cumulative interaction rate data + */ + InteractionRatesHomogeneous(std::string RateFile = "", std::string CumulativeRateFile = ""); + + std::vector getTabulatedEnergy() const; + std::vector getTabulatedRate() const; + std::vector getTabulatedE() const; + std::vector getTabulateds() const; + std::vector> getTabulatedCDF() const; + + double getProcessRate(const double E, const Vector3d &position) const; + void loadPerformInteractionTabs(const Vector3d &position, std::vector &tabE, std::vector &tabs, std::vector> &tabCDF) const; + + void setTabulatedEnergy (std::vector& tabEnergy); + void setTabulatedRate (std::vector& tabRate); + void setTabulatedE (std::vector& tabE); + void setTabulateds (std::vector& tabs); + void setTabulatedCDF (std::vector>& tabCDF); + + /** Loads the interaction rate + * This function loads the interaction rate + * @param filename The name of the file containing the interaction rates + */ + void initRate(std::string filename); + /** Loads the cumulative interaction rate + * This function loads the interaction rate + * @param filename The name of the file containing the interaction rates + */ + void initCumulativeRate(std::string filename); + +protected: + + // tabulated interaction rates 1/lambda(E) + std::vector tabEnergy; //!< electron energy in [J] + std::vector tabRate; //!< interaction rate in [1/m] + + // tabulated CDF(s_kin, E) = cumulative differential interaction rate + std::vector tabE; //!< electron energy in [J] + std::vector tabs; //!< s_kin = s - m^2 in [J**2] + std::vector> tabCDF; //!< cumulative interaction rate + +}; + +/** + @class InteractionRatePositionDependent + @brief Interaction rates decorator for tabulated position dependent photon fields. + */ +class InteractionRatesPositionDependent: public InteractionRates { + +public: + /** Constructor of InteractioNRatesPositionDependent + * @param RateFilePath Path containing the interaction rates files (* /Rate) + * @param CumulativeRateFilePath Path containing the cumulative interaction rates files (* /CumulativeRate) + * @param surface Closed surface to confine the grid nodes to be uploaded (optional) + */ + InteractionRatesPositionDependent(std::string RateFilePath = "", std::string CumulativeRateFilePath = "", ref_ptr surface = NULL); + + int findClosestGridPoint(const Vector3d &position) const; + + std::vector getTabulatedEnergy() const; + std::vector> getTabulatedRate() const; + std::vector getTabulatedE() const; + std::vector> getTabulateds() const; + std::vector>> getTabulatedCDF() const; + std::unordered_map getPhotonDict() const; + std::vector getClosestRate(const Vector3d &position) const; + std::vector getClosests(const Vector3d &position) const; + std::vector> getClosestCDF(const Vector3d &position) const; + + double getProcessRate(const double E, const Vector3d &position) const; + void loadPerformInteractionTabs(const Vector3d &position, std::vector &tabE, std::vector &tabs, std::vector> &tabCDF) const; + + void setTabulatedEnergy (std::vector& tabEnergy); + void setTabulatedRate (std::vector>& tabRate); + void setTabulatedE (std::vector& tabE); + void setTabulateds (std::vector>& tabs); + void setTabulatedCDF (std::vector>>& tabCDF); + void setPhotonDict (std::unordered_map& photonDict); + + /** Apply a surface that confine the position dependent photon field + * @param surface Closed surface to confine the grid nodes to be uploaded + */ + void setSurface(ref_ptr surface); + ref_ptr getSurface() const; + + /** Loads the interaction rate + * This function loads the position dependent interaction rate + * @param filepath The name of the folder containing the interaction rates (* /Rate) + */ + void initRate(std::string filepath); + /** Loads the interaction rate + * This function loads the cumulative position dependent interaction rate + * @param filepath The name of the folder containing the cumulative interaction rates (* /CumulativeRate) + */ + void initCumulativeRate(std::string filepath); + +protected: + + // tabulated interaction rates 1/lambda(E) + std::vector tabEnergy; //!< electron energy in [J], assuming the same energy binning in each node + std::vector> tabRate; //!< interaction rate in [1/m] + + // tabulated CDF(s_kin, E) = cumulative differential interaction rate + std::vector tabE; //!< electron energy in [J], assuming the same energy binning in each node + std::vector> tabs; //!< s_kin = s - m^2 in [J**2] + std::vector>> tabCDF; //!< cumulative interaction rate + std::unordered_map photonDict; //!< dictionary to link tables to spatial coordinates + + PointCloud cloud; //!< point cloud for nanoflann KD-tree + KDTree* tree = nullptr; //!< pointer to the KD Tree + ref_ptr surface; + +}; + +} // namespace crpropa + +#endif // CRPROPA_INTERACTIONRATES_H diff --git a/include/crpropa/PhotonBackground.h b/include/crpropa/PhotonBackground.h index 1048cda06..c9a8fcef6 100644 --- a/include/crpropa/PhotonBackground.h +++ b/include/crpropa/PhotonBackground.h @@ -3,9 +3,12 @@ #include "crpropa/Common.h" #include "crpropa/Referenced.h" +#include "crpropa/Vector3.h" +#include "crpropa/Geometry.h" #include #include +#include namespace crpropa { /** @@ -18,25 +21,29 @@ namespace crpropa { @brief Abstract base class for photon fields. */ class PhotonField: public Referenced { + public: + PhotonField() { this->fieldName = "AbstractPhotonField"; this->isRedshiftDependent = false; + this->isPositionDependent = false; + this->surface = nullptr; } - + /** returns comoving photon density [1/m^3]. multiply with (1+z^3) for physical number density. @param ePhoton photon energy [J] @param z redshift (if redshift dependent, default = 0.) */ - virtual double getPhotonDensity(double ePhoton, double z = 0.) const = 0; - virtual double getMinimumPhotonEnergy(double z) const = 0; - virtual double getMaximumPhotonEnergy(double z) const = 0; + virtual double getPhotonDensity(double ePhoton, double z = 0., const Vector3d &pos = Vector3d(0.,0.,0.)) const = 0; + virtual double getMinimumPhotonEnergy(double z, const Vector3d &pos = Vector3d(0.,0.,0.)) const = 0; + virtual double getMaximumPhotonEnergy(double z, const Vector3d &pos = Vector3d(0.,0.,0.)) const = 0; virtual std::string getFieldName() const { return this->fieldName; } - + /** returns overall comoving scaling factor (cf. CRPropa3-data/calc_scaling.py) @@ -45,18 +52,29 @@ class PhotonField: public Referenced { virtual double getRedshiftScaling(double z) const { return 1.; }; - + bool hasRedshiftDependence() const { return this->isRedshiftDependent; } - + + bool hasPositionDependence() const { + return this->isPositionDependent; + } + + ref_ptr getSurface() const { + return this->surface; + } + void setFieldName(std::string fieldName) { this->fieldName = fieldName; } - + protected: std::string fieldName; bool isRedshiftDependent; + bool isPositionDependent; + ref_ptr surface; + }; /** @@ -66,16 +84,16 @@ class PhotonField: public Referenced { This class reads photon field data from files; The first file must be a list of photon energies [J], named fieldName_photonEnergy.txt The second file must be a list of comoving photon field densities [1/m^3], named fieldName_photonDensity.txt - Optionally, a third file contains redshifts, named fieldName_redshift.txt + Optionally, a third file contains redshifts, named fieldName_redshift.txt. */ class TabularPhotonField: public PhotonField { public: TabularPhotonField(const std::string fieldName, const bool isRedshiftDependent = true); - double getPhotonDensity(double ePhoton, double z = 0.) const; + double getPhotonDensity(double ePhoton, double z = 0., const Vector3d &pos = Vector3d(0.,0.,0.)) const; double getRedshiftScaling(double z) const; - double getMinimumPhotonEnergy(double z) const; - double getMaximumPhotonEnergy(double z) const; + double getMinimumPhotonEnergy(double z, const Vector3d &pos = Vector3d(0.,0.,0.)) const; + double getMaximumPhotonEnergy(double z, const Vector3d &pos = Vector3d(0.,0.,0.)) const; protected: void readPhotonEnergy(std::string filePath); @@ -90,6 +108,39 @@ class TabularPhotonField: public PhotonField { std::vector redshiftScalings; }; +/** + @class TabularSpatialPhotonField + @brief Position dependent photon field decorator for tabulated photon fields. + + This class reads photon field data from files in the appropriate directory; + The first files must be lists of photon energies [J], named fieldName_photonEnergy.txt and contained in the subdirectory /photonEnegy/; + The second files must be lists of comoving photon field densities [1/m^3], named fieldName_photonDensity.txt and contained in the subdirectory /photonDensity/; + The generated files through the CRPropa procedure (https://crpropa.github.io/CRPropa3/pages/example_notebooks/custom_photonfield/custom-photon-field.html) have a different ordering: the energy bins from the larger to the lower. + No redshift dependence is available. + The surface is defined to include the nodes of the grid contained within. + */ +class TabularSpatialPhotonField: public PhotonField { +public: + TabularSpatialPhotonField(const std::string fieldName, ref_ptr surface = nullptr); + + double getPhotonDensity(double ePhoton = 0., double z = 0., const Vector3d &pos = Vector3d(0.,0.,0.)) const; + double getMinimumPhotonEnergy(double z, const Vector3d &pos = Vector3d(0.,0.,0.)) const; + double getMaximumPhotonEnergy(double z, const Vector3d &pos = Vector3d(0.,0.,0.)) const; + +protected: + std::vector readPhotonEnergy(std::string filePath); + std::vector readPhotonDensity(std::string filePath); + void checkInputData() const; + + /** Apply a surface that confine the position dependent photon field + * @param surface closed surface to confine the nodes of grid to be uploaded */ + void setSurface(ref_ptr surface); + + std::vector photonEnergies; // assuming all the nodes in the grid have the same energy binning + std::vector> photonDensity; + std::unordered_map photonDict; +}; + /** @class IRB_Kneiske04 @brief Extragalactic background light model from Kneiske et al. 2004 @@ -287,6 +338,32 @@ class URB_Nitu21: public TabularPhotonField { URB_Nitu21() : TabularPhotonField("URB_Nitu21", false) {} }; +/** + @class ISRF + @brief Interstellar radiation field model by Freudenreich et al. (1998) implemented in Porter et al. (2017) + + Source info: + DOI: + https://iopscience.iop.org/article/10.3847/1538-4357/aa844d + */ +class ISRF_Freudenreich98: public TabularSpatialPhotonField { +public: + ISRF_Freudenreich98(ref_ptr surface) : TabularSpatialPhotonField("ISRF_Freudenreich98", surface) {} +}; + +/** + @class ISRF + @brief Interstellar radiation field model by Robitaille et al. (2012) implemented in Porter et al. (2017) + + Source info: + DOI: + https://iopscience.iop.org/article/10.3847/1538-4357/aa844d + */ +class ISRF_Robitaille12: public TabularSpatialPhotonField { +public: + ISRF_Robitaille12(ref_ptr surface) : TabularSpatialPhotonField("ISRF_Robitaille12", surface) {} +}; + /** @class BlackbodyPhotonField @brief Photon field decorator for black body photon fields. @@ -294,9 +371,9 @@ class URB_Nitu21: public TabularPhotonField { class BlackbodyPhotonField: public PhotonField { public: BlackbodyPhotonField(const std::string fieldName, const double blackbodyTemperature); - double getPhotonDensity(double ePhoton, double z = 0.) const; - double getMinimumPhotonEnergy(double z) const; - double getMaximumPhotonEnergy(double z) const; + double getPhotonDensity(double ePhoton, double z = 0., const Vector3d &pos = Vector3d(0.,0.,0.)) const; + double getMinimumPhotonEnergy(double z, const Vector3d &pos = Vector3d(0.,0.,0.)) const; + double getMaximumPhotonEnergy(double z, const Vector3d &pos = Vector3d(0.,0.,0.)) const; void setQuantile(double q); protected: diff --git a/include/crpropa/module/EMDoublePairProduction.h b/include/crpropa/module/EMDoublePairProduction.h index 42b136e51..770bbbc3a 100644 --- a/include/crpropa/module/EMDoublePairProduction.h +++ b/include/crpropa/module/EMDoublePairProduction.h @@ -6,6 +6,8 @@ #include "crpropa/Module.h" #include "crpropa/PhotonBackground.h" +#include "crpropa/InteractionRates.h" +#include "crpropa/Geometry.h" namespace crpropa { /** @@ -23,6 +25,7 @@ namespace crpropa { Thinning is available. A thinning of 0 means that all particles are tracked. For the maximum thinning of 1, only a few representative particles are added to the list of secondaries. Note that for thinning>0 the output must contain the column "weights", which should be included in the post-processing. + The surface is defined to include the nodes of the grid contained within. */ class EMDoublePairProduction: public Module { private: @@ -30,24 +33,23 @@ class EMDoublePairProduction: public Module { bool haveElectrons; double limit; double thinning; + ref_ptr surface; std::string interactionTag = "EMDP"; - - // tabulated interaction rate 1/lambda(E) - std::vector tabEnergy; //!< electron energy in [J] - std::vector tabRate; //!< interaction rate in [1/m] - + ref_ptr interactionRates; public: /** Constructor + The object used to load, store and access to the interaction rates of the process is the interactionRates pointer. @param photonField target photon field @param haveElectrons if true, add secondary electrons as candidates @param thinning weighted sampling of secondaries (0: all particles are tracked; 1: maximum thinning) @param limit step size limit as fraction of mean free path + @param surface suface to enclose the grid nodes to be loaded */ - EMDoublePairProduction(ref_ptr photonField, bool haveElectrons = false, double thinning = 0, double limit = 0.1); - + EMDoublePairProduction(ref_ptr photonField, bool haveElectrons = false, double thinning = 0, double limit = 0.1, ref_ptr surface = nullptr); + // set the target photon field void setPhotonField(ref_ptr photonField); - + // decide if secondary electrons are added to the simulation void setHaveElectrons(bool haveElectrons); @@ -55,21 +57,42 @@ class EMDoublePairProduction: public Module { * @param limit fraction of the mean free path */ void setLimit(double limit); - + /** Apply thinning with a given thinning factor * @param thinning factor of thinning (0: no thinning, 1: maximum thinning) */ void setThinning(double thinning); + /** Apply a surface that confine the position dependent photon field region. + * @param surface closed surface to confine the grid to be uploaded + */ + void setSurface(ref_ptr surface); + ref_ptr getSurface() const; + + /** set a custom interaction rate + * With this function you can change the type of interaction rate, + * if you would for example like to change from a homogeneous to a position + * dependent interaction rate. + * @param intRates ref_ptr to a InteractionRates class + */ + void setInteractionRates(ref_ptr intRates); + ref_ptr getInteractionRates() const; + /** set a custom interaction tag to trace back this interaction * @param tag string that will be added to the candidate and output */ void setInteractionTag(std::string tag); std::string getInteractionTag() const; - void initRate(std::string filename); + /** Loads the interaction rate + * This function loads the interaction rate from a given file/folder. + * @param path The name of the file/folder containing the interaction rates + */ + void initRate(std::string path); + void process(Candidate *candidate) const; void performInteraction(Candidate *candidate) const; + }; /** @}*/ diff --git a/include/crpropa/module/EMInverseComptonScattering.h b/include/crpropa/module/EMInverseComptonScattering.h index f2a2c888e..baa4b7b4a 100644 --- a/include/crpropa/module/EMInverseComptonScattering.h +++ b/include/crpropa/module/EMInverseComptonScattering.h @@ -6,6 +6,8 @@ #include "crpropa/Module.h" #include "crpropa/PhotonBackground.h" +#include "crpropa/InteractionRates.h" +#include "crpropa/Geometry.h" namespace crpropa { /** @@ -23,6 +25,7 @@ namespace crpropa { Thinning is available. A thinning of 0 means that all particles are tracked. For the maximum thinning of 1, only a few representative particles are added to the list of secondaries. Note that for thinning>0 the output must contain the column "weights", which should be included in the post-processing. + The surface is defined to include the nodes of the grid contained within. */ class EMInverseComptonScattering: public Module { private: @@ -30,53 +33,73 @@ class EMInverseComptonScattering: public Module { bool havePhotons; double limit; double thinning; + ref_ptr surface; std::string interactionTag = "EMIC"; - - // tabulated interaction rate 1/lambda(E) - std::vector tabEnergy; //!< electron energy in [J] - std::vector tabRate; //!< interaction rate in [1/m] - - // tabulated CDF(s_kin, E) = cumulative differential interaction rate - std::vector tabE; //!< electron energy in [J] - std::vector tabs; //!< s_kin = s - m^2 in [J**2] - std::vector< std::vector > tabCDF; //!< cumulative interaction rate - + ref_ptr interactionRates; public: /** Constructor + The object used to load, store and access to the interaction rates of the process is the interactionRates pointer. @param photonField target photon field @param havePhotons if true, add secondary photons as candidates @param thinning weighted sampling of secondaries (0: all particles are tracked; 1: maximum thinning) @param limit step size limit as fraction of mean free path + @param surface suface to enclose the grid nodes to be loaded */ - EMInverseComptonScattering(ref_ptr photonField, bool havePhotons = false, double thinning = 0, double limit = 0.1); - - // set the target photon field + EMInverseComptonScattering(ref_ptr photonField, bool havePhotons = false, double thinning = 0, double limit = 0.1, ref_ptr surface = nullptr); + + // set the target photon field void setPhotonField(ref_ptr photonField); - + // decide if secondary photons are added to the simulation void setHavePhotons(bool havePhotons); - + /** limit the step to a fraction of the mean free path @param limit fraction of the mean free path, should be between 0 and 1 - */ + */ void setLimit(double limit); - + /** Apply thinning with a given thinning factor * @param thinning factor of thinning (0: no thinning, 1: maximum thinning) */ void setThinning(double thinning); - + + /** Apply a surface that confine the position dependent photon field + * @param surface closed surface to confine the grid nodes to be uploaded + */ + void setSurface(ref_ptr surface); + ref_ptr getSurface() const; + /** set a custom interaction tag to trace back this interaction * @param tag string that will be added to the candidate and output */ void setInteractionTag(std::string tag); std::string getInteractionTag() const; - void initRate(std::string filename); - void initCumulativeRate(std::string filename); + /** set a custom interaction rate + * With this function you can change the type of interaction rate, + * if you would for example like to change from a homogeneous to a position + * dependent interaction rate. + * @param intRates ref_ptr to a InteractionRates class + */ + void setInteractionRates(ref_ptr intRates); + ref_ptr getInteractionRates() const; + + /** Loads the interaction rate + * This function loads the interaction rate from a given file/folder. + * @param path The name of the file/folder containing the interaction rates + */ + void initRate(std::string path); + + /** Loads the cumulative interaction rate + * This function loads the interaction rate from a given file/folder. + * @param path The name of the file/folder containing the interaction rates + */ + void initCumulativeRate(std::string path); + void process(Candidate *candidate) const; void performInteraction(Candidate *candidate) const; + }; /** @}*/ diff --git a/include/crpropa/module/EMPairProduction.h b/include/crpropa/module/EMPairProduction.h index f2ab8a1f7..bf060868f 100644 --- a/include/crpropa/module/EMPairProduction.h +++ b/include/crpropa/module/EMPairProduction.h @@ -6,6 +6,8 @@ #include "crpropa/Module.h" #include "crpropa/PhotonBackground.h" +#include "crpropa/Geometry.h" +#include "crpropa/InteractionRates.h" namespace crpropa { @@ -25,6 +27,7 @@ namespace crpropa { Thinning is available. A thinning of 0 means that all particles are tracked. For the maximum thinning of 1, only a few representative particles are added to the list of secondaries. Note that for thinning>0 the output must contain the column "weights", which should be included in the post-processing. + The surface is defined to include the nodes of the grid contained within. */ class EMPairProduction: public Module { private: @@ -32,25 +35,19 @@ class EMPairProduction: public Module { bool haveElectrons; // add secondary electrons to simulation double limit; // limit the step to a fraction of the mean free path double thinning; // factor of the thinning (0: no thinning, 1: maximum thinning) + ref_ptr surface; // surface that includes the nodes in the photonField grid to be included std::string interactionTag = "EMPP"; - - // tabulated interaction rate 1/lambda(E) - std::vector tabEnergy; //!< electron energy in [J] - std::vector tabRate; //!< interaction rate in [1/m] - - // tabulated CDF(s_kin, E) = cumulative differential interaction rate - std::vector tabE; //!< electron energy in [J] - std::vector tabs; //!< s_kin = s - m^2 in [J**2] - std::vector< std::vector > tabCDF; //!< cumulative interaction rate - + ref_ptr interactionRates; public: /** Constructor + The object used to load, store and access to the interaction rates of the process is the interactionRates pointer. @param photonField target photon field @param haveElectrons if true, add secondary electrons as candidates @param thinning weighted sampling of secondaries (0: all particles are tracked; 1: maximum thinning) @param limit step size limit as fraction of mean free path + @param surface suface to enclose the grid nodes to be loaded */ - EMPairProduction(ref_ptr photonField, bool haveElectrons = false, double thinning = 0,double limit = 0.1); + EMPairProduction(ref_ptr photonField, bool haveElectrons = false, double thinning = 0, double limit = 0.1, ref_ptr surface = nullptr); // set the target photon field void setPhotonField(ref_ptr photonField); @@ -68,20 +65,42 @@ class EMPairProduction: public Module { */ void setThinning(double thinning); + /** Apply a surface that confine the position dependent photon field + * @param surface closed surface to confine the grid to be uploaded */ + void setSurface(ref_ptr surface); + ref_ptr getSurface() const; + /** set a custom interaction tag to trace back this interaction * @param tag string that will be added to the candidate and output - */ + */ void setInteractionTag(std::string tag); std::string getInteractionTag() const; - void initRate(std::string filename); - void initCumulativeRate(std::string filename); - + /** set a custom interaction rate + * With this function you can change the type of interaction rate, + * if you would for example like to change from a homogeneous to a position + * dependent interaction rate. + * @param intRates ref_ptr to a InteractionRates class + */ + void setInteractionRates(ref_ptr intRates); + ref_ptr getInteractionRates() const; + + /** Loads the interaction rate + * This function loads the interaction rate from a given file/folder. + * @param path The name of the file/folder containing the interaction rates + */ + void initRate(std::string path); + + /** Loads the cumulative interaction rate + * This function loads the interaction rate from a given file/folder. + * @param path The name of the file/folder containing the interaction rates + */ + void initCumulativeRate(std::string path); + void performInteraction(Candidate *candidate) const; void process(Candidate *candidate) const; + }; -/** @}*/ - } // namespace crpropa diff --git a/include/crpropa/module/EMTripletPairProduction.h b/include/crpropa/module/EMTripletPairProduction.h index c689fc2ba..959a52b65 100644 --- a/include/crpropa/module/EMTripletPairProduction.h +++ b/include/crpropa/module/EMTripletPairProduction.h @@ -6,6 +6,8 @@ #include "crpropa/Module.h" #include "crpropa/PhotonBackground.h" +#include "crpropa/Geometry.h" +#include "crpropa/InteractionRates.h" namespace crpropa { /** @@ -23,6 +25,7 @@ namespace crpropa { Thinning is available. A thinning of 0 means that all particles are tracked. For the maximum thinning of 1, only a few representative particles are added to the list of secondaries. Note that for thinning>0 the output must contain the column "weights", which should be included in the post-processing. + The surface is defined to include the nodes of the grid contained within. */ class EMTripletPairProduction: public Module { private: @@ -30,25 +33,19 @@ class EMTripletPairProduction: public Module { bool haveElectrons; double limit; double thinning; + ref_ptr surface; std::string interactionTag = "EMTP"; - - // tabulated interaction rate 1/lambda(E) - std::vector tabEnergy; //!< electron energy in [J] - std::vector tabRate; //!< interaction rate in [1/m] - - // tabulated CDF(s_kin, E) = cumulative differential interaction rate - std::vector tabE; //!< electron energy in [J] - std::vector tabs; //!< s_kin = s - m^2 in [J**2] - std::vector< std::vector > tabCDF; //!< cumulative interaction rate - + ref_ptr interactionRates; public: /** Constructor + The object used to load, store and access to the interaction rates of the process is the interactionRates pointer. @param photonField target photon field @param haveElectrons if true, add secondary electrons as candidates @param thinning weighted sampling of secondaries (0: all particles are tracked; 1: maximum thinning) @param limit step size limit as fraction of mean free path + @param surface suface to enclose the grid nodes to be loaded */ - EMTripletPairProduction(ref_ptr photonField, bool haveElectrons = false, double thinning = 0, double limit = 0.1); + EMTripletPairProduction(ref_ptr photonField, bool haveElectrons = false, double thinning = 0, double limit = 0.1, ref_ptr surface = nullptr); // set the target photon field void setPhotonField(ref_ptr photonField); @@ -65,16 +62,40 @@ class EMTripletPairProduction: public Module { * @param thinning factor of thinning (0: no thinning, 1: maximum thinning) */ void setThinning(double thinning); - + + /** Apply a surface that confine the position dependent photon field + * @param surface closed surface to confine the grid to be uploaded + */ + void setSurface(ref_ptr surface); + ref_ptr getSurface() const; + /** set a custom interaction tag to trace back this interaction * @param tag string that will be added to the candidate and output */ void setInteractionTag(std::string tag); std::string getInteractionTag() const; - void initRate(std::string filename); - void initCumulativeRate(std::string filename); - + /** set a custom interaction rate + * With this function you can change the type of interaction rate, + * if you would for example like to change from a homogeneous to a position + * dependent interaction rate. + * @param intRates ref_ptr to a InteractionRates class + */ + void setInteractionRates(ref_ptr intRates); + ref_ptr getInteractionRates() const; + + /** Loads the interaction rate + * This function loads the interaction rate from a given file/folder. + * @param path The name of the file/folder containing the interaction rates + */ + void initRate(std::string path); + + /** Loads the cumulative interaction rate + * This function loads the interaction rate from a given file/folder. + * @param path The name of the file/folder containing the interaction rates + */ + void initCumulativeRate(std::string path); + void process(Candidate *candidate) const; void performInteraction(Candidate *candidate) const; diff --git a/libs/nanoflann/include/nanoflann.hpp b/libs/nanoflann/include/nanoflann.hpp new file mode 100644 index 000000000..503f9fae1 --- /dev/null +++ b/libs/nanoflann/include/nanoflann.hpp @@ -0,0 +1,2720 @@ +/*********************************************************************** + * Software License Agreement (BSD License) + * + * Copyright 2008-2009 Marius Muja (mariusm@cs.ubc.ca). All rights reserved. + * Copyright 2008-2009 David G. Lowe (lowe@cs.ubc.ca). All rights reserved. + * Copyright 2011-2025 Jose Luis Blanco (joseluisblancoc@gmail.com). + * All rights reserved. + * + * THE BSD LICENSE + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR + * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. + * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT + * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF + * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + *************************************************************************/ + +/** \mainpage nanoflann C++ API documentation + * nanoflann is a C++ header-only library for building KD-Trees, mostly + * optimized for 2D or 3D point clouds. + * + * nanoflann does not require compiling or installing, just an + * #include in your code. + * + * See: + * - [Online README](https://github.com/jlblancoc/nanoflann) + * - [C++ API documentation](https://jlblancoc.github.io/nanoflann/) + */ + +#pragma once + +#include +#include +#include +#include +#include // for abs() +#include +#include // for abs() +#include // std::reference_wrapper +#include +#include +#include // std::numeric_limits +#include +#include +#include +#include + +/** Library version: 0xMmP (M=Major,m=minor,P=patch) */ +#define NANOFLANN_VERSION 0x170 + +// Avoid conflicting declaration of min/max macros in Windows headers +#if !defined(NOMINMAX) && \ + (defined(_WIN32) || defined(_WIN32_) || defined(WIN32) || defined(_WIN64)) +#define NOMINMAX +#ifdef max +#undef max +#undef min +#endif +#endif +// Avoid conflicts with X11 headers +#ifdef None +#undef None +#endif + +namespace nanoflann +{ +/** @addtogroup nanoflann_grp nanoflann C++ library for KD-trees + * @{ */ + +/** the PI constant (required to avoid MSVC missing symbols) */ +template +T pi_const() +{ + return static_cast(3.14159265358979323846); +} + +/** + * Traits if object is resizable and assignable (typically has a resize | assign + * method) + */ +template +struct has_resize : std::false_type +{ +}; + +template +struct has_resize().resize(1), 0)> + : std::true_type +{ +}; + +template +struct has_assign : std::false_type +{ +}; + +template +struct has_assign().assign(1, 0), 0)> + : std::true_type +{ +}; + +/** + * Free function to resize a resizable object + */ +template +inline typename std::enable_if::value, void>::type resize( + Container& c, const size_t nElements) +{ + c.resize(nElements); +} + +/** + * Free function that has no effects on non resizable containers (e.g. + * std::array) It raises an exception if the expected size does not match + */ +template +inline typename std::enable_if::value, void>::type + resize(Container& c, const size_t nElements) +{ + if (nElements != c.size()) + throw std::logic_error("Try to change the size of a std::array."); +} + +/** + * Free function to assign to a container + */ +template +inline typename std::enable_if::value, void>::type assign( + Container& c, const size_t nElements, const T& value) +{ + c.assign(nElements, value); +} + +/** + * Free function to assign to a std::array + */ +template +inline typename std::enable_if::value, void>::type + assign(Container& c, const size_t nElements, const T& value) +{ + for (size_t i = 0; i < nElements; i++) c[i] = value; +} + +/** operator "<" for std::sort() */ +struct IndexDist_Sorter +{ + /** PairType will be typically: ResultItem */ + template + bool operator()(const PairType& p1, const PairType& p2) const + { + return p1.second < p2.second; + } +}; + +/** + * Each result element in RadiusResultSet. Note that distances and indices + * are named `first` and `second` to keep backward-compatibility with the + * `std::pair<>` type used in the past. In contrast, this structure is ensured + * to be `std::is_standard_layout` so it can be used in wrappers to other + * languages. + * See: https://github.com/jlblancoc/nanoflann/issues/166 + */ +template +struct ResultItem +{ + ResultItem() = default; + ResultItem(const IndexType index, const DistanceType distance) + : first(index), second(distance) + { + } + + IndexType first; //!< Index of the sample in the dataset + DistanceType second; //!< Distance from sample to query point +}; + +/** @addtogroup result_sets_grp Result set classes + * @{ */ + +/** Result set for KNN searches (N-closest neighbors) */ +template < + typename _DistanceType, typename _IndexType = size_t, + typename _CountType = size_t> +class KNNResultSet +{ + public: + using DistanceType = _DistanceType; + using IndexType = _IndexType; + using CountType = _CountType; + + private: + IndexType* indices; + DistanceType* dists; + CountType capacity; + CountType count; + + public: + explicit KNNResultSet(CountType capacity_) + : indices(nullptr), dists(nullptr), capacity(capacity_), count(0) + { + } + + void init(IndexType* indices_, DistanceType* dists_) + { + indices = indices_; + dists = dists_; + count = 0; + } + + CountType size() const { return count; } + bool empty() const { return count == 0; } + bool full() const { return count == capacity; } + + /** + * Called during search to add an element matching the criteria. + * @return true if the search should be continued, false if the results are + * sufficient + */ + bool addPoint(DistanceType dist, IndexType index) + { + CountType i; + for (i = count; i > 0; --i) + { + /** If defined and two points have the same distance, the one with + * the lowest-index will be returned first. */ +#ifdef NANOFLANN_FIRST_MATCH + if ((dists[i - 1] > dist) || + ((dist == dists[i - 1]) && (indices[i - 1] > index))) + { +#else + if (dists[i - 1] > dist) + { +#endif + if (i < capacity) + { + dists[i] = dists[i - 1]; + indices[i] = indices[i - 1]; + } + } + else + break; + } + if (i < capacity) + { + dists[i] = dist; + indices[i] = index; + } + if (count < capacity) count++; + + // tell caller that the search shall continue + return true; + } + + //! Returns the worst distance among found solutions if the search result is + //! full, or the maximum possible distance, if not full yet. + DistanceType worstDist() const + { + return count < capacity ? std::numeric_limits::max() + : dists[count - 1]; + } + + void sort() + { + // already sorted + } +}; + +/** Result set for RKNN searches (N-closest neighbors with a maximum radius) */ +template < + typename _DistanceType, typename _IndexType = size_t, + typename _CountType = size_t> +class RKNNResultSet +{ + public: + using DistanceType = _DistanceType; + using IndexType = _IndexType; + using CountType = _CountType; + + private: + IndexType* indices; + DistanceType* dists; + CountType capacity; + CountType count; + DistanceType maximumSearchDistanceSquared; + + public: + explicit RKNNResultSet( + CountType capacity_, DistanceType maximumSearchDistanceSquared_) + : indices(nullptr), + dists(nullptr), + capacity(capacity_), + count(0), + maximumSearchDistanceSquared(maximumSearchDistanceSquared_) + { + } + + void init(IndexType* indices_, DistanceType* dists_) + { + indices = indices_; + dists = dists_; + count = 0; + if (capacity) dists[capacity - 1] = maximumSearchDistanceSquared; + } + + CountType size() const { return count; } + bool empty() const { return count == 0; } + bool full() const { return count == capacity; } + + /** + * Called during search to add an element matching the criteria. + * @return true if the search should be continued, false if the results are + * sufficient + */ + bool addPoint(DistanceType dist, IndexType index) + { + CountType i; + for (i = count; i > 0; --i) + { + /** If defined and two points have the same distance, the one with + * the lowest-index will be returned first. */ +#ifdef NANOFLANN_FIRST_MATCH + if ((dists[i - 1] > dist) || + ((dist == dists[i - 1]) && (indices[i - 1] > index))) + { +#else + if (dists[i - 1] > dist) + { +#endif + if (i < capacity) + { + dists[i] = dists[i - 1]; + indices[i] = indices[i - 1]; + } + } + else + break; + } + if (i < capacity) + { + dists[i] = dist; + indices[i] = index; + } + if (count < capacity) count++; + + // tell caller that the search shall continue + return true; + } + + //! Returns the worst distance among found solutions if the search result is + //! full, or the maximum possible distance, if not full yet. + DistanceType worstDist() const + { + return count < capacity ? maximumSearchDistanceSquared + : dists[count - 1]; + } + + void sort() + { + // already sorted + } +}; + +/** + * A result-set class used when performing a radius based search. + */ +template +class RadiusResultSet +{ + public: + using DistanceType = _DistanceType; + using IndexType = _IndexType; + + public: + const DistanceType radius; + + std::vector>& m_indices_dists; + + explicit RadiusResultSet( + DistanceType radius_, + std::vector>& indices_dists) + : radius(radius_), m_indices_dists(indices_dists) + { + init(); + } + + void init() { clear(); } + void clear() { m_indices_dists.clear(); } + + size_t size() const { return m_indices_dists.size(); } + size_t empty() const { return m_indices_dists.empty(); } + + bool full() const { return true; } + + /** + * Called during search to add an element matching the criteria. + * @return true if the search should be continued, false if the results are + * sufficient + */ + bool addPoint(DistanceType dist, IndexType index) + { + if (dist < radius) m_indices_dists.emplace_back(index, dist); + return true; + } + + DistanceType worstDist() const { return radius; } + + /** + * Find the worst result (farthest neighbor) without copying or sorting + * Pre-conditions: size() > 0 + */ + ResultItem worst_item() const + { + if (m_indices_dists.empty()) + throw std::runtime_error( + "Cannot invoke RadiusResultSet::worst_item() on " + "an empty list of results."); + auto it = std::max_element( + m_indices_dists.begin(), m_indices_dists.end(), IndexDist_Sorter()); + return *it; + } + + void sort() + { + std::sort( + m_indices_dists.begin(), m_indices_dists.end(), IndexDist_Sorter()); + } +}; + +/** @} */ + +/** @addtogroup loadsave_grp Load/save auxiliary functions + * @{ */ +template +void save_value(std::ostream& stream, const T& value) +{ + stream.write(reinterpret_cast(&value), sizeof(T)); +} + +template +void save_value(std::ostream& stream, const std::vector& value) +{ + size_t size = value.size(); + stream.write(reinterpret_cast(&size), sizeof(size_t)); + stream.write(reinterpret_cast(value.data()), sizeof(T) * size); +} + +template +void load_value(std::istream& stream, T& value) +{ + stream.read(reinterpret_cast(&value), sizeof(T)); +} + +template +void load_value(std::istream& stream, std::vector& value) +{ + size_t size; + stream.read(reinterpret_cast(&size), sizeof(size_t)); + value.resize(size); + stream.read(reinterpret_cast(value.data()), sizeof(T) * size); +} +/** @} */ + +/** @addtogroup metric_grp Metric (distance) classes + * @{ */ + +struct Metric +{ +}; + +/** Manhattan distance functor (generic version, optimized for + * high-dimensionality data sets). Corresponding distance traits: + * nanoflann::metric_L1 + * + * \tparam T Type of the elements (e.g. double, float, uint8_t) + * \tparam DataSource Source of the data, i.e. where the vectors are stored + * \tparam _DistanceType Type of distance variables (must be signed) + * \tparam IndexType Type of the arguments with which the data can be + * accessed (e.g. float, double, int64_t, T*) + */ +template < + class T, class DataSource, typename _DistanceType = T, + typename IndexType = uint32_t> +struct L1_Adaptor +{ + using ElementType = T; + using DistanceType = _DistanceType; + + const DataSource& data_source; + + L1_Adaptor(const DataSource& _data_source) : data_source(_data_source) {} + + DistanceType evalMetric( + const T* a, const IndexType b_idx, size_t size, + DistanceType worst_dist = -1) const + { + DistanceType result = DistanceType(); + const T* last = a + size; + const T* lastgroup = last - 3; + size_t d = 0; + + /* Process 4 items with each loop for efficiency. */ + while (a < lastgroup) + { + const DistanceType diff0 = + std::abs(a[0] - data_source.kdtree_get_pt(b_idx, d++)); + const DistanceType diff1 = + std::abs(a[1] - data_source.kdtree_get_pt(b_idx, d++)); + const DistanceType diff2 = + std::abs(a[2] - data_source.kdtree_get_pt(b_idx, d++)); + const DistanceType diff3 = + std::abs(a[3] - data_source.kdtree_get_pt(b_idx, d++)); + result += diff0 + diff1 + diff2 + diff3; + a += 4; + if ((worst_dist > 0) && (result > worst_dist)) { return result; } + } + /* Process last 0-3 components. Not needed for standard vector lengths. + */ + while (a < last) + { + result += std::abs(*a++ - data_source.kdtree_get_pt(b_idx, d++)); + } + return result; + } + + template + DistanceType accum_dist(const U a, const V b, const size_t) const + { + return std::abs(a - b); + } +}; + +/** **Squared** Euclidean distance functor (generic version, optimized for + * high-dimensionality data sets). Corresponding distance traits: + * nanoflann::metric_L2 + * + * \tparam T Type of the elements (e.g. double, float, uint8_t) + * \tparam DataSource Source of the data, i.e. where the vectors are stored + * \tparam _DistanceType Type of distance variables (must be signed) + * \tparam IndexType Type of the arguments with which the data can be + * accessed (e.g. float, double, int64_t, T*) + */ +template < + class T, class DataSource, typename _DistanceType = T, + typename IndexType = uint32_t> +struct L2_Adaptor +{ + using ElementType = T; + using DistanceType = _DistanceType; + + const DataSource& data_source; + + L2_Adaptor(const DataSource& _data_source) : data_source(_data_source) {} + + DistanceType evalMetric( + const T* a, const IndexType b_idx, size_t size, + DistanceType worst_dist = -1) const + { + DistanceType result = DistanceType(); + const T* last = a + size; + const T* lastgroup = last - 3; + size_t d = 0; + + /* Process 4 items with each loop for efficiency. */ + while (a < lastgroup) + { + const DistanceType diff0 = + a[0] - data_source.kdtree_get_pt(b_idx, d++); + const DistanceType diff1 = + a[1] - data_source.kdtree_get_pt(b_idx, d++); + const DistanceType diff2 = + a[2] - data_source.kdtree_get_pt(b_idx, d++); + const DistanceType diff3 = + a[3] - data_source.kdtree_get_pt(b_idx, d++); + result += + diff0 * diff0 + diff1 * diff1 + diff2 * diff2 + diff3 * diff3; + a += 4; + if ((worst_dist > 0) && (result > worst_dist)) { return result; } + } + /* Process last 0-3 components. Not needed for standard vector lengths. + */ + while (a < last) + { + const DistanceType diff0 = + *a++ - data_source.kdtree_get_pt(b_idx, d++); + result += diff0 * diff0; + } + return result; + } + + template + DistanceType accum_dist(const U a, const V b, const size_t) const + { + return (a - b) * (a - b); + } +}; + +/** **Squared** Euclidean (L2) distance functor (suitable for low-dimensionality + * datasets, like 2D or 3D point clouds) Corresponding distance traits: + * nanoflann::metric_L2_Simple + * + * \tparam T Type of the elements (e.g. double, float, uint8_t) + * \tparam DataSource Source of the data, i.e. where the vectors are stored + * \tparam _DistanceType Type of distance variables (must be signed) + * \tparam IndexType Type of the arguments with which the data can be + * accessed (e.g. float, double, int64_t, T*) + */ +template < + class T, class DataSource, typename _DistanceType = T, + typename IndexType = uint32_t> +struct L2_Simple_Adaptor +{ + using ElementType = T; + using DistanceType = _DistanceType; + + const DataSource& data_source; + + L2_Simple_Adaptor(const DataSource& _data_source) + : data_source(_data_source) + { + } + + DistanceType evalMetric( + const T* a, const IndexType b_idx, size_t size) const + { + DistanceType result = DistanceType(); + for (size_t i = 0; i < size; ++i) + { + const DistanceType diff = + a[i] - data_source.kdtree_get_pt(b_idx, i); + result += diff * diff; + } + return result; + } + + template + DistanceType accum_dist(const U a, const V b, const size_t) const + { + return (a - b) * (a - b); + } +}; + +/** SO2 distance functor + * Corresponding distance traits: nanoflann::metric_SO2 + * + * \tparam T Type of the elements (e.g. double, float, uint8_t) + * \tparam DataSource Source of the data, i.e. where the vectors are stored + * \tparam _DistanceType Type of distance variables (must be signed) (e.g. + * float, double) orientation is constrained to be in [-pi, pi] + * \tparam IndexType Type of the arguments with which the data can be + * accessed (e.g. float, double, int64_t, T*) + */ +template < + class T, class DataSource, typename _DistanceType = T, + typename IndexType = uint32_t> +struct SO2_Adaptor +{ + using ElementType = T; + using DistanceType = _DistanceType; + + const DataSource& data_source; + + SO2_Adaptor(const DataSource& _data_source) : data_source(_data_source) {} + + DistanceType evalMetric( + const T* a, const IndexType b_idx, size_t size) const + { + return accum_dist( + a[size - 1], data_source.kdtree_get_pt(b_idx, size - 1), size - 1); + } + + /** Note: this assumes that input angles are already in the range [-pi,pi] + */ + template + DistanceType accum_dist(const U a, const V b, const size_t) const + { + DistanceType result = DistanceType(); + DistanceType PI = pi_const(); + result = b - a; + if (result > PI) + result -= 2 * PI; + else if (result < -PI) + result += 2 * PI; + return result; + } +}; + +/** SO3 distance functor (Uses L2_Simple) + * Corresponding distance traits: nanoflann::metric_SO3 + * + * \tparam T Type of the elements (e.g. double, float, uint8_t) + * \tparam DataSource Source of the data, i.e. where the vectors are stored + * \tparam _DistanceType Type of distance variables (must be signed) (e.g. + * float, double) + * \tparam IndexType Type of the arguments with which the data can be + * accessed (e.g. float, double, int64_t, T*) + */ +template < + class T, class DataSource, typename _DistanceType = T, + typename IndexType = uint32_t> +struct SO3_Adaptor +{ + using ElementType = T; + using DistanceType = _DistanceType; + + L2_Simple_Adaptor + distance_L2_Simple; + + SO3_Adaptor(const DataSource& _data_source) + : distance_L2_Simple(_data_source) + { + } + + DistanceType evalMetric( + const T* a, const IndexType b_idx, size_t size) const + { + return distance_L2_Simple.evalMetric(a, b_idx, size); + } + + template + DistanceType accum_dist(const U a, const V b, const size_t idx) const + { + return distance_L2_Simple.accum_dist(a, b, idx); + } +}; + +/** Metaprogramming helper traits class for the L1 (Manhattan) metric */ +struct metric_L1 : public Metric +{ + template + struct traits + { + using distance_t = L1_Adaptor; + }; +}; +/** Metaprogramming helper traits class for the L2 (Euclidean) **squared** + * distance metric */ +struct metric_L2 : public Metric +{ + template + struct traits + { + using distance_t = L2_Adaptor; + }; +}; +/** Metaprogramming helper traits class for the L2_simple (Euclidean) + * **squared** distance metric */ +struct metric_L2_Simple : public Metric +{ + template + struct traits + { + using distance_t = L2_Simple_Adaptor; + }; +}; +/** Metaprogramming helper traits class for the SO3_InnerProdQuat metric */ +struct metric_SO2 : public Metric +{ + template + struct traits + { + using distance_t = SO2_Adaptor; + }; +}; +/** Metaprogramming helper traits class for the SO3_InnerProdQuat metric */ +struct metric_SO3 : public Metric +{ + template + struct traits + { + using distance_t = SO3_Adaptor; + }; +}; + +/** @} */ + +/** @addtogroup param_grp Parameter structs + * @{ */ + +enum class KDTreeSingleIndexAdaptorFlags +{ + None = 0, + SkipInitialBuildIndex = 1 +}; + +inline std::underlying_type::type operator&( + KDTreeSingleIndexAdaptorFlags lhs, KDTreeSingleIndexAdaptorFlags rhs) +{ + using underlying = + typename std::underlying_type::type; + return static_cast(lhs) & static_cast(rhs); +} + +/** Parameters (see README.md) */ +struct KDTreeSingleIndexAdaptorParams +{ + KDTreeSingleIndexAdaptorParams( + size_t _leaf_max_size = 10, + KDTreeSingleIndexAdaptorFlags _flags = + KDTreeSingleIndexAdaptorFlags::None, + unsigned int _n_thread_build = 1) + : leaf_max_size(_leaf_max_size), + flags(_flags), + n_thread_build(_n_thread_build) + { + } + + size_t leaf_max_size; + KDTreeSingleIndexAdaptorFlags flags; + unsigned int n_thread_build; +}; + +/** Search options for KDTreeSingleIndexAdaptor::findNeighbors() */ +struct SearchParameters +{ + SearchParameters(float eps_ = 0, bool sorted_ = true) + : eps(eps_), sorted(sorted_) + { + } + + float eps; //!< search for eps-approximate neighbours (default: 0) + bool sorted; //!< only for radius search, require neighbours sorted by + //!< distance (default: true) +}; +/** @} */ + +/** @addtogroup memalloc_grp Memory allocation + * @{ */ + +/** + * Pooled storage allocator + * + * The following routines allow for the efficient allocation of storage in + * small chunks from a specified pool. Rather than allowing each structure + * to be freed individually, an entire pool of storage is freed at once. + * This method has two advantages over just using malloc() and free(). First, + * it is far more efficient for allocating small objects, as there is + * no overhead for remembering all the information needed to free each + * object or consolidating fragmented memory. Second, the decision about + * how long to keep an object is made at the time of allocation, and there + * is no need to track down all the objects to free them. + * + */ +class PooledAllocator +{ + static constexpr size_t WORDSIZE = 16; // WORDSIZE must >= 8 + static constexpr size_t BLOCKSIZE = 8192; + + /* We maintain memory alignment to word boundaries by requiring that all + allocations be in multiples of the machine wordsize. */ + /* Size of machine word in bytes. Must be power of 2. */ + /* Minimum number of bytes requested at a time from the system. Must be + * multiple of WORDSIZE. */ + + using Size = size_t; + + Size remaining_ = 0; //!< Number of bytes left in current block of storage + void* base_ = nullptr; //!< Pointer to base of current block of storage + void* loc_ = nullptr; //!< Current location in block to next allocate + + void internal_init() + { + remaining_ = 0; + base_ = nullptr; + usedMemory = 0; + wastedMemory = 0; + } + + public: + Size usedMemory = 0; + Size wastedMemory = 0; + + /** + Default constructor. Initializes a new pool. + */ + PooledAllocator() { internal_init(); } + + /** + * Destructor. Frees all the memory allocated in this pool. + */ + ~PooledAllocator() { free_all(); } + + /** Frees all allocated memory chunks */ + void free_all() + { + while (base_ != nullptr) + { + // Get pointer to prev block + void* prev = *(static_cast(base_)); + ::free(base_); + base_ = prev; + } + internal_init(); + } + + /** + * Returns a pointer to a piece of new memory of the given size in bytes + * allocated from the pool. + */ + void* malloc(const size_t req_size) + { + /* Round size up to a multiple of wordsize. The following expression + only works for WORDSIZE that is a power of 2, by masking last bits + of incremented size to zero. + */ + const Size size = (req_size + (WORDSIZE - 1)) & ~(WORDSIZE - 1); + + /* Check whether a new block must be allocated. Note that the first + word of a block is reserved for a pointer to the previous block. + */ + if (size > remaining_) + { + wastedMemory += remaining_; + + /* Allocate new storage. */ + const Size blocksize = + size > BLOCKSIZE ? size + WORDSIZE : BLOCKSIZE + WORDSIZE; + + // use the standard C malloc to allocate memory + void* m = ::malloc(blocksize); + if (!m) + { + fprintf(stderr, "Failed to allocate memory.\n"); + throw std::bad_alloc(); + } + + /* Fill first word of new block with pointer to previous block. */ + static_cast(m)[0] = base_; + base_ = m; + + remaining_ = blocksize - WORDSIZE; + loc_ = static_cast(m) + WORDSIZE; + } + void* rloc = loc_; + loc_ = static_cast(loc_) + size; + remaining_ -= size; + + usedMemory += size; + + return rloc; + } + + /** + * Allocates (using this pool) a generic type T. + * + * Params: + * count = number of instances to allocate. + * Returns: pointer (of type T*) to memory buffer + */ + template + T* allocate(const size_t count = 1) + { + T* mem = static_cast(this->malloc(sizeof(T) * count)); + return mem; + } +}; +/** @} */ + +/** @addtogroup nanoflann_metaprog_grp Auxiliary metaprogramming stuff + * @{ */ + +/** Used to declare fixed-size arrays when DIM>0, dynamically-allocated vectors + * when DIM=-1. Fixed size version for a generic DIM: + */ +template +struct array_or_vector +{ + using type = std::array; +}; +/** Dynamic size version */ +template +struct array_or_vector<-1, T> +{ + using type = std::vector; +}; + +/** @} */ + +/** kd-tree base-class + * + * Contains the member functions common to the classes KDTreeSingleIndexAdaptor + * and KDTreeSingleIndexDynamicAdaptor_. + * + * \tparam Derived The name of the class which inherits this class. + * \tparam DatasetAdaptor The user-provided adaptor, which must be ensured to + * have a lifetime equal or longer than the instance of this class. + * \tparam Distance The distance metric to use, these are all classes derived + * from nanoflann::Metric + * \tparam DIM Dimensionality of data points (e.g. 3 for 3D points) + * \tparam IndexType Type of the arguments with which the data can be + * accessed (e.g. float, double, int64_t, T*) + */ +template < + class Derived, typename Distance, class DatasetAdaptor, int32_t DIM = -1, + typename index_t = uint32_t> +class KDTreeBaseClass +{ + public: + /** Frees the previously-built index. Automatically called within + * buildIndex(). */ + void freeIndex(Derived& obj) + { + obj.pool_.free_all(); + obj.root_node_ = nullptr; + obj.size_at_index_build_ = 0; + } + + using ElementType = typename Distance::ElementType; + using DistanceType = typename Distance::DistanceType; + using IndexType = index_t; + + /** + * Array of indices to vectors in the dataset_. + */ + std::vector vAcc_; + + using Offset = typename decltype(vAcc_)::size_type; + using Size = typename decltype(vAcc_)::size_type; + using Dimension = int32_t; + + /*--------------------------- + * Internal Data Structures + * --------------------------*/ + struct Node + { + /** Union used because a node can be either a LEAF node or a non-leaf + * node, so both data fields are never used simultaneously */ + union + { + struct leaf + { + Offset left, right; //!< Indices of points in leaf node + } lr; + struct nonleaf + { + Dimension divfeat; //!< Dimension used for subdivision. + /// The values used for subdivision. + DistanceType divlow, divhigh; + } sub; + } node_type; + + /** Child nodes (both=nullptr mean its a leaf node) */ + Node *child1 = nullptr, *child2 = nullptr; + }; + + using NodePtr = Node*; + using NodeConstPtr = const Node*; + + struct Interval + { + ElementType low, high; + }; + + NodePtr root_node_ = nullptr; + + Size leaf_max_size_ = 0; + + /// Number of thread for concurrent tree build + Size n_thread_build_ = 1; + /// Number of current points in the dataset + Size size_ = 0; + /// Number of points in the dataset when the index was built + Size size_at_index_build_ = 0; + Dimension dim_ = 0; //!< Dimensionality of each data point + + /** Define "BoundingBox" as a fixed-size or variable-size container + * depending on "DIM" */ + using BoundingBox = typename array_or_vector::type; + + /** Define "distance_vector_t" as a fixed-size or variable-size container + * depending on "DIM" */ + using distance_vector_t = typename array_or_vector::type; + + /** The KD-tree used to find neighbours */ + BoundingBox root_bbox_; + + /** + * Pooled memory allocator. + * + * Using a pooled memory allocator is more efficient + * than allocating memory directly when there is a large + * number small of memory allocations. + */ + PooledAllocator pool_; + + /** Returns number of points in dataset */ + Size size(const Derived& obj) const { return obj.size_; } + + /** Returns the length of each point in the dataset */ + Size veclen(const Derived& obj) { return DIM > 0 ? DIM : obj.dim; } + + /// Helper accessor to the dataset points: + ElementType dataset_get( + const Derived& obj, IndexType element, Dimension component) const + { + return obj.dataset_.kdtree_get_pt(element, component); + } + + /** + * Computes the inde memory usage + * Returns: memory used by the index + */ + Size usedMemory(Derived& obj) + { + return obj.pool_.usedMemory + obj.pool_.wastedMemory + + obj.dataset_.kdtree_get_point_count() * + sizeof(IndexType); // pool memory and vind array memory + } + + void computeMinMax( + const Derived& obj, Offset ind, Size count, Dimension element, + ElementType& min_elem, ElementType& max_elem) + { + min_elem = dataset_get(obj, vAcc_[ind], element); + max_elem = min_elem; + for (Offset i = 1; i < count; ++i) + { + ElementType val = dataset_get(obj, vAcc_[ind + i], element); + if (val < min_elem) min_elem = val; + if (val > max_elem) max_elem = val; + } + } + + /** + * Create a tree node that subdivides the list of vecs from vind[first] + * to vind[last]. The routine is called recursively on each sublist. + * + * @param left index of the first vector + * @param right index of the last vector + */ + NodePtr divideTree( + Derived& obj, const Offset left, const Offset right, BoundingBox& bbox) + { + assert(left < obj.dataset_.kdtree_get_point_count()); + + NodePtr node = obj.pool_.template allocate(); // allocate memory + const auto dims = (DIM > 0 ? DIM : obj.dim_); + + /* If too few exemplars remain, then make this a leaf node. */ + if ((right - left) <= static_cast(obj.leaf_max_size_)) + { + node->child1 = node->child2 = nullptr; /* Mark as leaf node. */ + node->node_type.lr.left = left; + node->node_type.lr.right = right; + + // compute bounding-box of leaf points + for (Dimension i = 0; i < dims; ++i) + { + bbox[i].low = dataset_get(obj, obj.vAcc_[left], i); + bbox[i].high = dataset_get(obj, obj.vAcc_[left], i); + } + for (Offset k = left + 1; k < right; ++k) + { + for (Dimension i = 0; i < dims; ++i) + { + const auto val = dataset_get(obj, obj.vAcc_[k], i); + if (bbox[i].low > val) bbox[i].low = val; + if (bbox[i].high < val) bbox[i].high = val; + } + } + } + else + { + Offset idx; + Dimension cutfeat; + DistanceType cutval; + middleSplit_(obj, left, right - left, idx, cutfeat, cutval, bbox); + + node->node_type.sub.divfeat = cutfeat; + + BoundingBox left_bbox(bbox); + left_bbox[cutfeat].high = cutval; + node->child1 = this->divideTree(obj, left, left + idx, left_bbox); + + BoundingBox right_bbox(bbox); + right_bbox[cutfeat].low = cutval; + node->child2 = this->divideTree(obj, left + idx, right, right_bbox); + + node->node_type.sub.divlow = left_bbox[cutfeat].high; + node->node_type.sub.divhigh = right_bbox[cutfeat].low; + + for (Dimension i = 0; i < dims; ++i) + { + bbox[i].low = std::min(left_bbox[i].low, right_bbox[i].low); + bbox[i].high = std::max(left_bbox[i].high, right_bbox[i].high); + } + } + + return node; + } + + /** + * Create a tree node that subdivides the list of vecs from vind[first] to + * vind[last] concurrently. The routine is called recursively on each + * sublist. + * + * @param left index of the first vector + * @param right index of the last vector + * @param thread_count count of std::async threads + * @param mutex mutex for mempool allocation + */ + NodePtr divideTreeConcurrent( + Derived& obj, const Offset left, const Offset right, BoundingBox& bbox, + std::atomic& thread_count, std::mutex& mutex) + { + std::unique_lock lock(mutex); + NodePtr node = obj.pool_.template allocate(); // allocate memory + lock.unlock(); + + const auto dims = (DIM > 0 ? DIM : obj.dim_); + + /* If too few exemplars remain, then make this a leaf node. */ + if ((right - left) <= static_cast(obj.leaf_max_size_)) + { + node->child1 = node->child2 = nullptr; /* Mark as leaf node. */ + node->node_type.lr.left = left; + node->node_type.lr.right = right; + + // compute bounding-box of leaf points + for (Dimension i = 0; i < dims; ++i) + { + bbox[i].low = dataset_get(obj, obj.vAcc_[left], i); + bbox[i].high = dataset_get(obj, obj.vAcc_[left], i); + } + for (Offset k = left + 1; k < right; ++k) + { + for (Dimension i = 0; i < dims; ++i) + { + const auto val = dataset_get(obj, obj.vAcc_[k], i); + if (bbox[i].low > val) bbox[i].low = val; + if (bbox[i].high < val) bbox[i].high = val; + } + } + } + else + { + Offset idx; + Dimension cutfeat; + DistanceType cutval; + middleSplit_(obj, left, right - left, idx, cutfeat, cutval, bbox); + + node->node_type.sub.divfeat = cutfeat; + + std::future right_future; + + BoundingBox right_bbox(bbox); + right_bbox[cutfeat].low = cutval; + if (++thread_count < n_thread_build_) + { + // Concurrent right sub-tree + right_future = std::async( + std::launch::async, &KDTreeBaseClass::divideTreeConcurrent, + this, std::ref(obj), left + idx, right, + std::ref(right_bbox), std::ref(thread_count), + std::ref(mutex)); + } + else { --thread_count; } + + BoundingBox left_bbox(bbox); + left_bbox[cutfeat].high = cutval; + node->child1 = this->divideTreeConcurrent( + obj, left, left + idx, left_bbox, thread_count, mutex); + + if (right_future.valid()) + { + // Block and wait for concurrent right sub-tree + node->child2 = right_future.get(); + --thread_count; + } + else + { + node->child2 = this->divideTreeConcurrent( + obj, left + idx, right, right_bbox, thread_count, mutex); + } + + node->node_type.sub.divlow = left_bbox[cutfeat].high; + node->node_type.sub.divhigh = right_bbox[cutfeat].low; + + for (Dimension i = 0; i < dims; ++i) + { + bbox[i].low = std::min(left_bbox[i].low, right_bbox[i].low); + bbox[i].high = std::max(left_bbox[i].high, right_bbox[i].high); + } + } + + return node; + } + + void middleSplit_( + const Derived& obj, const Offset ind, const Size count, Offset& index, + Dimension& cutfeat, DistanceType& cutval, const BoundingBox& bbox) + { + const auto dims = (DIM > 0 ? DIM : obj.dim_); + const auto EPS = static_cast(0.00001); + ElementType max_span = bbox[0].high - bbox[0].low; + for (Dimension i = 1; i < dims; ++i) + { + ElementType span = bbox[i].high - bbox[i].low; + if (span > max_span) { max_span = span; } + } + ElementType max_spread = -1; + cutfeat = 0; + ElementType min_elem = 0, max_elem = 0; + for (Dimension i = 0; i < dims; ++i) + { + ElementType span = bbox[i].high - bbox[i].low; + if (span >= (1 - EPS) * max_span) + { + ElementType min_elem_, max_elem_; + computeMinMax(obj, ind, count, i, min_elem_, max_elem_); + ElementType spread = max_elem_ - min_elem_; + if (spread > max_spread) + { + cutfeat = i; + max_spread = spread; + min_elem = min_elem_; + max_elem = max_elem_; + } + } + } + // split in the middle + DistanceType split_val = (bbox[cutfeat].low + bbox[cutfeat].high) / 2; + + if (split_val < min_elem) + cutval = min_elem; + else if (split_val > max_elem) + cutval = max_elem; + else + cutval = split_val; + + Offset lim1, lim2; + planeSplit(obj, ind, count, cutfeat, cutval, lim1, lim2); + + if (lim1 > count / 2) + index = lim1; + else if (lim2 < count / 2) + index = lim2; + else + index = count / 2; + } + + /** + * Subdivide the list of points by a plane perpendicular on the axis + * corresponding to the 'cutfeat' dimension at 'cutval' position. + * + * On return: + * dataset[ind[0..lim1-1]][cutfeat]cutval + */ + void planeSplit( + const Derived& obj, const Offset ind, const Size count, + const Dimension cutfeat, const DistanceType& cutval, Offset& lim1, + Offset& lim2) + { + /* Move vector indices for left subtree to front of list. */ + Offset left = 0; + Offset right = count - 1; + for (;;) + { + while (left <= right && + dataset_get(obj, vAcc_[ind + left], cutfeat) < cutval) + ++left; + while (right && left <= right && + dataset_get(obj, vAcc_[ind + right], cutfeat) >= cutval) + --right; + if (left > right || !right) + break; // "!right" was added to support unsigned Index types + std::swap(vAcc_[ind + left], vAcc_[ind + right]); + ++left; + --right; + } + /* If either list is empty, it means that all remaining features + * are identical. Split in the middle to maintain a balanced tree. + */ + lim1 = left; + right = count - 1; + for (;;) + { + while (left <= right && + dataset_get(obj, vAcc_[ind + left], cutfeat) <= cutval) + ++left; + while (right && left <= right && + dataset_get(obj, vAcc_[ind + right], cutfeat) > cutval) + --right; + if (left > right || !right) + break; // "!right" was added to support unsigned Index types + std::swap(vAcc_[ind + left], vAcc_[ind + right]); + ++left; + --right; + } + lim2 = left; + } + + DistanceType computeInitialDistances( + const Derived& obj, const ElementType* vec, + distance_vector_t& dists) const + { + assert(vec); + DistanceType dist = DistanceType(); + + for (Dimension i = 0; i < (DIM > 0 ? DIM : obj.dim_); ++i) + { + if (vec[i] < obj.root_bbox_[i].low) + { + dists[i] = + obj.distance_.accum_dist(vec[i], obj.root_bbox_[i].low, i); + dist += dists[i]; + } + if (vec[i] > obj.root_bbox_[i].high) + { + dists[i] = + obj.distance_.accum_dist(vec[i], obj.root_bbox_[i].high, i); + dist += dists[i]; + } + } + return dist; + } + + static void save_tree( + const Derived& obj, std::ostream& stream, const NodeConstPtr tree) + { + save_value(stream, *tree); + if (tree->child1 != nullptr) { save_tree(obj, stream, tree->child1); } + if (tree->child2 != nullptr) { save_tree(obj, stream, tree->child2); } + } + + static void load_tree(Derived& obj, std::istream& stream, NodePtr& tree) + { + tree = obj.pool_.template allocate(); + load_value(stream, *tree); + if (tree->child1 != nullptr) { load_tree(obj, stream, tree->child1); } + if (tree->child2 != nullptr) { load_tree(obj, stream, tree->child2); } + } + + /** Stores the index in a binary file. + * IMPORTANT NOTE: The set of data points is NOT stored in the file, so + * when loading the index object it must be constructed associated to the + * same source of data points used while building it. See the example: + * examples/saveload_example.cpp \sa loadIndex */ + void saveIndex(const Derived& obj, std::ostream& stream) const + { + save_value(stream, obj.size_); + save_value(stream, obj.dim_); + save_value(stream, obj.root_bbox_); + save_value(stream, obj.leaf_max_size_); + save_value(stream, obj.vAcc_); + if (obj.root_node_) save_tree(obj, stream, obj.root_node_); + } + + /** Loads a previous index from a binary file. + * IMPORTANT NOTE: The set of data points is NOT stored in the file, so + * the index object must be constructed associated to the same source of + * data points used while building the index. See the example: + * examples/saveload_example.cpp \sa loadIndex */ + void loadIndex(Derived& obj, std::istream& stream) + { + load_value(stream, obj.size_); + load_value(stream, obj.dim_); + load_value(stream, obj.root_bbox_); + load_value(stream, obj.leaf_max_size_); + load_value(stream, obj.vAcc_); + load_tree(obj, stream, obj.root_node_); + } +}; + +/** @addtogroup kdtrees_grp KD-tree classes and adaptors + * @{ */ + +/** kd-tree static index + * + * Contains the k-d trees and other information for indexing a set of points + * for nearest-neighbor matching. + * + * The class "DatasetAdaptor" must provide the following interface (can be + * non-virtual, inlined methods): + * + * \code + * // Must return the number of data poins + * size_t kdtree_get_point_count() const { ... } + * + * + * // Must return the dim'th component of the idx'th point in the class: + * T kdtree_get_pt(const size_t idx, const size_t dim) const { ... } + * + * // Optional bounding-box computation: return false to default to a standard + * bbox computation loop. + * // Return true if the BBOX was already computed by the class and returned + * in "bb" so it can be avoided to redo it again. + * // Look at bb.size() to find out the expected dimensionality (e.g. 2 or 3 + * for point clouds) template bool kdtree_get_bbox(BBOX &bb) const + * { + * bb[0].low = ...; bb[0].high = ...; // 0th dimension limits + * bb[1].low = ...; bb[1].high = ...; // 1st dimension limits + * ... + * return true; + * } + * + * \endcode + * + * \tparam DatasetAdaptor The user-provided adaptor, which must be ensured to + * have a lifetime equal or longer than the instance of this class. + * \tparam Distance The distance metric to use: nanoflann::metric_L1, + * nanoflann::metric_L2, nanoflann::metric_L2_Simple, etc. \tparam DIM + * Dimensionality of data points (e.g. 3 for 3D points) \tparam IndexType Will + * be typically size_t or int + */ +template < + typename Distance, class DatasetAdaptor, int32_t DIM = -1, + typename index_t = uint32_t> +class KDTreeSingleIndexAdaptor + : public KDTreeBaseClass< + KDTreeSingleIndexAdaptor, + Distance, DatasetAdaptor, DIM, index_t> +{ + public: + /** Deleted copy constructor*/ + explicit KDTreeSingleIndexAdaptor( + const KDTreeSingleIndexAdaptor< + Distance, DatasetAdaptor, DIM, index_t>&) = delete; + + /** The data source used by this index */ + const DatasetAdaptor& dataset_; + + const KDTreeSingleIndexAdaptorParams indexParams; + + Distance distance_; + + using Base = typename nanoflann::KDTreeBaseClass< + nanoflann::KDTreeSingleIndexAdaptor< + Distance, DatasetAdaptor, DIM, index_t>, + Distance, DatasetAdaptor, DIM, index_t>; + + using Offset = typename Base::Offset; + using Size = typename Base::Size; + using Dimension = typename Base::Dimension; + + using ElementType = typename Base::ElementType; + using DistanceType = typename Base::DistanceType; + using IndexType = typename Base::IndexType; + + using Node = typename Base::Node; + using NodePtr = Node*; + + using Interval = typename Base::Interval; + + /** Define "BoundingBox" as a fixed-size or variable-size container + * depending on "DIM" */ + using BoundingBox = typename Base::BoundingBox; + + /** Define "distance_vector_t" as a fixed-size or variable-size container + * depending on "DIM" */ + using distance_vector_t = typename Base::distance_vector_t; + + /** + * KDTree constructor + * + * Refer to docs in README.md or online in + * https://github.com/jlblancoc/nanoflann + * + * The KD-Tree point dimension (the length of each point in the datase, e.g. + * 3 for 3D points) is determined by means of: + * - The \a DIM template parameter if >0 (highest priority) + * - Otherwise, the \a dimensionality parameter of this constructor. + * + * @param inputData Dataset with the input features. Its lifetime must be + * equal or longer than that of the instance of this class. + * @param params Basically, the maximum leaf node size + * + * Note that there is a variable number of optional additional parameters + * which will be forwarded to the metric class constructor. Refer to example + * `examples/pointcloud_custom_metric.cpp` for a use case. + * + */ + template + explicit KDTreeSingleIndexAdaptor( + const Dimension dimensionality, const DatasetAdaptor& inputData, + const KDTreeSingleIndexAdaptorParams& params, Args&&... args) + : dataset_(inputData), + indexParams(params), + distance_(inputData, std::forward(args)...) + { + init(dimensionality, params); + } + + explicit KDTreeSingleIndexAdaptor( + const Dimension dimensionality, const DatasetAdaptor& inputData, + const KDTreeSingleIndexAdaptorParams& params = {}) + : dataset_(inputData), indexParams(params), distance_(inputData) + { + init(dimensionality, params); + } + + private: + void init( + const Dimension dimensionality, + const KDTreeSingleIndexAdaptorParams& params) + { + Base::size_ = dataset_.kdtree_get_point_count(); + Base::size_at_index_build_ = Base::size_; + Base::dim_ = dimensionality; + if (DIM > 0) Base::dim_ = DIM; + Base::leaf_max_size_ = params.leaf_max_size; + if (params.n_thread_build > 0) + { + Base::n_thread_build_ = params.n_thread_build; + } + else + { + Base::n_thread_build_ = + std::max(std::thread::hardware_concurrency(), 1u); + } + + if (!(params.flags & + KDTreeSingleIndexAdaptorFlags::SkipInitialBuildIndex)) + { + // Build KD-tree: + buildIndex(); + } + } + + public: + /** + * Builds the index + */ + void buildIndex() + { + Base::size_ = dataset_.kdtree_get_point_count(); + Base::size_at_index_build_ = Base::size_; + init_vind(); + this->freeIndex(*this); + Base::size_at_index_build_ = Base::size_; + if (Base::size_ == 0) return; + computeBoundingBox(Base::root_bbox_); + // construct the tree + if (Base::n_thread_build_ == 1) + { + Base::root_node_ = + this->divideTree(*this, 0, Base::size_, Base::root_bbox_); + } + else + { +#ifndef NANOFLANN_NO_THREADS + std::atomic thread_count(0u); + std::mutex mutex; + Base::root_node_ = this->divideTreeConcurrent( + *this, 0, Base::size_, Base::root_bbox_, thread_count, mutex); +#else /* NANOFLANN_NO_THREADS */ + throw std::runtime_error("Multithreading is disabled"); +#endif /* NANOFLANN_NO_THREADS */ + } + } + + /** \name Query methods + * @{ */ + + /** + * Find set of nearest neighbors to vec[0:dim-1]. Their indices are stored + * inside the result object. + * + * Params: + * result = the result object in which the indices of the + * nearest-neighbors are stored vec = the vector for which to search the + * nearest neighbors + * + * \tparam RESULTSET Should be any ResultSet + * \return True if the requested neighbors could be found. + * \sa knnSearch, radiusSearch + * + * \note If L2 norms are used, all returned distances are actually squared + * distances. + */ + template + bool findNeighbors( + RESULTSET& result, const ElementType* vec, + const SearchParameters& searchParams = {}) const + { + assert(vec); + if (this->size(*this) == 0) return false; + if (!Base::root_node_) + throw std::runtime_error( + "[nanoflann] findNeighbors() called before building the " + "index."); + float epsError = 1 + searchParams.eps; + + // fixed or variable-sized container (depending on DIM) + distance_vector_t dists; + // Fill it with zeros. + auto zero = static_cast(0); + assign(dists, (DIM > 0 ? DIM : Base::dim_), zero); + DistanceType dist = this->computeInitialDistances(*this, vec, dists); + searchLevel(result, vec, Base::root_node_, dist, dists, epsError); + + if (searchParams.sorted) result.sort(); + + return result.full(); + } + + /** + * Find the "num_closest" nearest neighbors to the \a query_point[0:dim-1]. + * Their indices and distances are stored in the provided pointers to + * array/vector. + * + * \sa radiusSearch, findNeighbors + * \return Number `N` of valid points in the result set. + * + * \note If L2 norms are used, all returned distances are actually squared + * distances. + * + * \note Only the first `N` entries in `out_indices` and `out_distances` + * will be valid. Return is less than `num_closest` only if the + * number of elements in the tree is less than `num_closest`. + */ + Size knnSearch( + const ElementType* query_point, const Size num_closest, + IndexType* out_indices, DistanceType* out_distances) const + { + nanoflann::KNNResultSet resultSet(num_closest); + resultSet.init(out_indices, out_distances); + findNeighbors(resultSet, query_point); + return resultSet.size(); + } + + /** + * Find all the neighbors to \a query_point[0:dim-1] within a maximum + * radius. The output is given as a vector of pairs, of which the first + * element is a point index and the second the corresponding distance. + * Previous contents of \a IndicesDists are cleared. + * + * If searchParams.sorted==true, the output list is sorted by ascending + * distances. + * + * For a better performance, it is advisable to do a .reserve() on the + * vector if you have any wild guess about the number of expected matches. + * + * \sa knnSearch, findNeighbors, radiusSearchCustomCallback + * \return The number of points within the given radius (i.e. indices.size() + * or dists.size() ) + * + * \note If L2 norms are used, search radius and all returned distances + * are actually squared distances. + */ + Size radiusSearch( + const ElementType* query_point, const DistanceType& radius, + std::vector>& IndicesDists, + const SearchParameters& searchParams = {}) const + { + RadiusResultSet resultSet( + radius, IndicesDists); + const Size nFound = + radiusSearchCustomCallback(query_point, resultSet, searchParams); + return nFound; + } + + /** + * Just like radiusSearch() but with a custom callback class for each point + * found in the radius of the query. See the source of RadiusResultSet<> as + * a start point for your own classes. \sa radiusSearch + */ + template + Size radiusSearchCustomCallback( + const ElementType* query_point, SEARCH_CALLBACK& resultSet, + const SearchParameters& searchParams = {}) const + { + findNeighbors(resultSet, query_point, searchParams); + return resultSet.size(); + } + + /** + * Find the first N neighbors to \a query_point[0:dim-1] within a maximum + * radius. The output is given as a vector of pairs, of which the first + * element is a point index and the second the corresponding distance. + * Previous contents of \a IndicesDists are cleared. + * + * \sa radiusSearch, findNeighbors + * \return Number `N` of valid points in the result set. + * + * \note If L2 norms are used, all returned distances are actually squared + * distances. + * + * \note Only the first `N` entries in `out_indices` and `out_distances` + * will be valid. Return is less than `num_closest` only if the + * number of elements in the tree is less than `num_closest`. + */ + Size rknnSearch( + const ElementType* query_point, const Size num_closest, + IndexType* out_indices, DistanceType* out_distances, + const DistanceType& radius) const + { + nanoflann::RKNNResultSet resultSet( + num_closest, radius); + resultSet.init(out_indices, out_distances); + findNeighbors(resultSet, query_point); + return resultSet.size(); + } + + /** @} */ + + public: + /** Make sure the auxiliary list \a vind has the same size than the current + * dataset, and re-generate if size has changed. */ + void init_vind() + { + // Create a permutable array of indices to the input vectors. + Base::size_ = dataset_.kdtree_get_point_count(); + if (Base::vAcc_.size() != Base::size_) Base::vAcc_.resize(Base::size_); + for (IndexType i = 0; i < static_cast(Base::size_); i++) + Base::vAcc_[i] = i; + } + + void computeBoundingBox(BoundingBox& bbox) + { + const auto dims = (DIM > 0 ? DIM : Base::dim_); + resize(bbox, dims); + if (dataset_.kdtree_get_bbox(bbox)) + { + // Done! It was implemented in derived class + } + else + { + const Size N = dataset_.kdtree_get_point_count(); + if (!N) + throw std::runtime_error( + "[nanoflann] computeBoundingBox() called but " + "no data points found."); + for (Dimension i = 0; i < dims; ++i) + { + bbox[i].low = bbox[i].high = + this->dataset_get(*this, Base::vAcc_[0], i); + } + for (Offset k = 1; k < N; ++k) + { + for (Dimension i = 0; i < dims; ++i) + { + const auto val = + this->dataset_get(*this, Base::vAcc_[k], i); + if (val < bbox[i].low) bbox[i].low = val; + if (val > bbox[i].high) bbox[i].high = val; + } + } + } + } + + /** + * Performs an exact search in the tree starting from a node. + * \tparam RESULTSET Should be any ResultSet + * \return true if the search should be continued, false if the results are + * sufficient + */ + template + bool searchLevel( + RESULTSET& result_set, const ElementType* vec, const NodePtr node, + DistanceType mindist, distance_vector_t& dists, + const float epsError) const + { + /* If this is a leaf node, then do check and return. */ + if ((node->child1 == nullptr) && (node->child2 == nullptr)) + { + DistanceType worst_dist = result_set.worstDist(); + for (Offset i = node->node_type.lr.left; + i < node->node_type.lr.right; ++i) + { + const IndexType accessor = Base::vAcc_[i]; // reorder... : i; + DistanceType dist = distance_.evalMetric( + vec, accessor, (DIM > 0 ? DIM : Base::dim_)); + if (dist < worst_dist) + { + if (!result_set.addPoint(dist, Base::vAcc_[i])) + { + // the resultset doesn't want to receive any more + // points, we're done searching! + return false; + } + } + } + return true; + } + + /* Which child branch should be taken first? */ + Dimension idx = node->node_type.sub.divfeat; + ElementType val = vec[idx]; + DistanceType diff1 = val - node->node_type.sub.divlow; + DistanceType diff2 = val - node->node_type.sub.divhigh; + + NodePtr bestChild; + NodePtr otherChild; + DistanceType cut_dist; + if ((diff1 + diff2) < 0) + { + bestChild = node->child1; + otherChild = node->child2; + cut_dist = + distance_.accum_dist(val, node->node_type.sub.divhigh, idx); + } + else + { + bestChild = node->child2; + otherChild = node->child1; + cut_dist = + distance_.accum_dist(val, node->node_type.sub.divlow, idx); + } + + /* Call recursively to search next level down. */ + if (!searchLevel(result_set, vec, bestChild, mindist, dists, epsError)) + { + // the resultset doesn't want to receive any more points, we're done + // searching! + return false; + } + + DistanceType dst = dists[idx]; + mindist = mindist + cut_dist - dst; + dists[idx] = cut_dist; + if (mindist * epsError <= result_set.worstDist()) + { + if (!searchLevel( + result_set, vec, otherChild, mindist, dists, epsError)) + { + // the resultset doesn't want to receive any more points, we're + // done searching! + return false; + } + } + dists[idx] = dst; + return true; + } + + public: + /** Stores the index in a binary file. + * IMPORTANT NOTE: The set of data points is NOT stored in the file, so + * when loading the index object it must be constructed associated to the + * same source of data points used while building it. See the example: + * examples/saveload_example.cpp \sa loadIndex */ + void saveIndex(std::ostream& stream) const + { + Base::saveIndex(*this, stream); + } + + /** Loads a previous index from a binary file. + * IMPORTANT NOTE: The set of data points is NOT stored in the file, so + * the index object must be constructed associated to the same source of + * data points used while building the index. See the example: + * examples/saveload_example.cpp \sa loadIndex */ + void loadIndex(std::istream& stream) { Base::loadIndex(*this, stream); } + +}; // class KDTree + +/** kd-tree dynamic index + * + * Contains the k-d trees and other information for indexing a set of points + * for nearest-neighbor matching. + * + * The class "DatasetAdaptor" must provide the following interface (can be + * non-virtual, inlined methods): + * + * \code + * // Must return the number of data poins + * size_t kdtree_get_point_count() const { ... } + * + * // Must return the dim'th component of the idx'th point in the class: + * T kdtree_get_pt(const size_t idx, const size_t dim) const { ... } + * + * // Optional bounding-box computation: return false to default to a standard + * bbox computation loop. + * // Return true if the BBOX was already computed by the class and returned + * in "bb" so it can be avoided to redo it again. + * // Look at bb.size() to find out the expected dimensionality (e.g. 2 or 3 + * for point clouds) template bool kdtree_get_bbox(BBOX &bb) const + * { + * bb[0].low = ...; bb[0].high = ...; // 0th dimension limits + * bb[1].low = ...; bb[1].high = ...; // 1st dimension limits + * ... + * return true; + * } + * + * \endcode + * + * \tparam DatasetAdaptor The user-provided adaptor (see comments above). + * \tparam Distance The distance metric to use: nanoflann::metric_L1, + * nanoflann::metric_L2, nanoflann::metric_L2_Simple, etc. + * \tparam DIM Dimensionality of data points (e.g. 3 for 3D points) + * \tparam IndexType Type of the arguments with which the data can be + * accessed (e.g. float, double, int64_t, T*) + */ +template < + typename Distance, class DatasetAdaptor, int32_t DIM = -1, + typename IndexType = uint32_t> +class KDTreeSingleIndexDynamicAdaptor_ + : public KDTreeBaseClass< + KDTreeSingleIndexDynamicAdaptor_< + Distance, DatasetAdaptor, DIM, IndexType>, + Distance, DatasetAdaptor, DIM, IndexType> +{ + public: + /** + * The dataset used by this index + */ + const DatasetAdaptor& dataset_; //!< The source of our data + + KDTreeSingleIndexAdaptorParams index_params_; + + std::vector& treeIndex_; + + Distance distance_; + + using Base = typename nanoflann::KDTreeBaseClass< + nanoflann::KDTreeSingleIndexDynamicAdaptor_< + Distance, DatasetAdaptor, DIM, IndexType>, + Distance, DatasetAdaptor, DIM, IndexType>; + + using ElementType = typename Base::ElementType; + using DistanceType = typename Base::DistanceType; + + using Offset = typename Base::Offset; + using Size = typename Base::Size; + using Dimension = typename Base::Dimension; + + using Node = typename Base::Node; + using NodePtr = Node*; + + using Interval = typename Base::Interval; + /** Define "BoundingBox" as a fixed-size or variable-size container + * depending on "DIM" */ + using BoundingBox = typename Base::BoundingBox; + + /** Define "distance_vector_t" as a fixed-size or variable-size container + * depending on "DIM" */ + using distance_vector_t = typename Base::distance_vector_t; + + /** + * KDTree constructor + * + * Refer to docs in README.md or online in + * https://github.com/jlblancoc/nanoflann + * + * The KD-Tree point dimension (the length of each point in the datase, e.g. + * 3 for 3D points) is determined by means of: + * - The \a DIM template parameter if >0 (highest priority) + * - Otherwise, the \a dimensionality parameter of this constructor. + * + * @param inputData Dataset with the input features. Its lifetime must be + * equal or longer than that of the instance of this class. + * @param params Basically, the maximum leaf node size + */ + KDTreeSingleIndexDynamicAdaptor_( + const Dimension dimensionality, const DatasetAdaptor& inputData, + std::vector& treeIndex, + const KDTreeSingleIndexAdaptorParams& params = + KDTreeSingleIndexAdaptorParams()) + : dataset_(inputData), + index_params_(params), + treeIndex_(treeIndex), + distance_(inputData) + { + Base::size_ = 0; + Base::size_at_index_build_ = 0; + for (auto& v : Base::root_bbox_) v = {}; + Base::dim_ = dimensionality; + if (DIM > 0) Base::dim_ = DIM; + Base::leaf_max_size_ = params.leaf_max_size; + if (params.n_thread_build > 0) + { + Base::n_thread_build_ = params.n_thread_build; + } + else + { + Base::n_thread_build_ = + std::max(std::thread::hardware_concurrency(), 1u); + } + } + + /** Explicitly default the copy constructor */ + KDTreeSingleIndexDynamicAdaptor_( + const KDTreeSingleIndexDynamicAdaptor_& rhs) = default; + + /** Assignment operator definiton */ + KDTreeSingleIndexDynamicAdaptor_ operator=( + const KDTreeSingleIndexDynamicAdaptor_& rhs) + { + KDTreeSingleIndexDynamicAdaptor_ tmp(rhs); + std::swap(Base::vAcc_, tmp.Base::vAcc_); + std::swap(Base::leaf_max_size_, tmp.Base::leaf_max_size_); + std::swap(index_params_, tmp.index_params_); + std::swap(treeIndex_, tmp.treeIndex_); + std::swap(Base::size_, tmp.Base::size_); + std::swap(Base::size_at_index_build_, tmp.Base::size_at_index_build_); + std::swap(Base::root_node_, tmp.Base::root_node_); + std::swap(Base::root_bbox_, tmp.Base::root_bbox_); + std::swap(Base::pool_, tmp.Base::pool_); + return *this; + } + + /** + * Builds the index + */ + void buildIndex() + { + Base::size_ = Base::vAcc_.size(); + this->freeIndex(*this); + Base::size_at_index_build_ = Base::size_; + if (Base::size_ == 0) return; + computeBoundingBox(Base::root_bbox_); + // construct the tree + if (Base::n_thread_build_ == 1) + { + Base::root_node_ = + this->divideTree(*this, 0, Base::size_, Base::root_bbox_); + } + else + { +#ifndef NANOFLANN_NO_THREADS + std::atomic thread_count(0u); + std::mutex mutex; + Base::root_node_ = this->divideTreeConcurrent( + *this, 0, Base::size_, Base::root_bbox_, thread_count, mutex); +#else /* NANOFLANN_NO_THREADS */ + throw std::runtime_error("Multithreading is disabled"); +#endif /* NANOFLANN_NO_THREADS */ + } + } + + /** \name Query methods + * @{ */ + + /** + * Find set of nearest neighbors to vec[0:dim-1]. Their indices are stored + * inside the result object. + * This is the core search function, all others are wrappers around this + * one. + * + * \param result The result object in which the indices of the + * nearest-neighbors are stored. + * \param vec The vector of the query point for which to search the + * nearest neighbors. + * \param searchParams Optional parameters for the search. + * + * \tparam RESULTSET Should be any ResultSet + * \return True if the requested neighbors could be found. + * + * \sa knnSearch(), radiusSearch(), radiusSearchCustomCallback() + * + * \note If L2 norms are used, all returned distances are actually squared + * distances. + */ + template + bool findNeighbors( + RESULTSET& result, const ElementType* vec, + const SearchParameters& searchParams = {}) const + { + assert(vec); + if (this->size(*this) == 0) return false; + if (!Base::root_node_) return false; + float epsError = 1 + searchParams.eps; + + // fixed or variable-sized container (depending on DIM) + distance_vector_t dists; + // Fill it with zeros. + assign( + dists, (DIM > 0 ? DIM : Base::dim_), + static_cast(0)); + DistanceType dist = this->computeInitialDistances(*this, vec, dists); + searchLevel(result, vec, Base::root_node_, dist, dists, epsError); + return result.full(); + } + + /** + * Find the "num_closest" nearest neighbors to the \a query_point[0:dim-1]. + * Their indices are stored inside the result object. \sa radiusSearch, + * findNeighbors + * \return Number `N` of valid points in + * the result set. + * + * \note If L2 norms are used, all returned distances are actually squared + * distances. + * + * \note Only the first `N` entries in `out_indices` and `out_distances` + * will be valid. Return may be less than `num_closest` only if the + * number of elements in the tree is less than `num_closest`. + */ + Size knnSearch( + const ElementType* query_point, const Size num_closest, + IndexType* out_indices, DistanceType* out_distances, + const SearchParameters& searchParams = {}) const + { + nanoflann::KNNResultSet resultSet(num_closest); + resultSet.init(out_indices, out_distances); + findNeighbors(resultSet, query_point, searchParams); + return resultSet.size(); + } + + /** + * Find all the neighbors to \a query_point[0:dim-1] within a maximum + * radius. The output is given as a vector of pairs, of which the first + * element is a point index and the second the corresponding distance. + * Previous contents of \a IndicesDists are cleared. + * + * If searchParams.sorted==true, the output list is sorted by ascending + * distances. + * + * For a better performance, it is advisable to do a .reserve() on the + * vector if you have any wild guess about the number of expected matches. + * + * \sa knnSearch, findNeighbors, radiusSearchCustomCallback + * \return The number of points within the given radius (i.e. indices.size() + * or dists.size() ) + * + * \note If L2 norms are used, search radius and all returned distances + * are actually squared distances. + */ + Size radiusSearch( + const ElementType* query_point, const DistanceType& radius, + std::vector>& IndicesDists, + const SearchParameters& searchParams = {}) const + { + RadiusResultSet resultSet( + radius, IndicesDists); + const size_t nFound = + radiusSearchCustomCallback(query_point, resultSet, searchParams); + return nFound; + } + + /** + * Just like radiusSearch() but with a custom callback class for each point + * found in the radius of the query. See the source of RadiusResultSet<> as + * a start point for your own classes. \sa radiusSearch + */ + template + Size radiusSearchCustomCallback( + const ElementType* query_point, SEARCH_CALLBACK& resultSet, + const SearchParameters& searchParams = {}) const + { + findNeighbors(resultSet, query_point, searchParams); + return resultSet.size(); + } + + /** @} */ + + public: + void computeBoundingBox(BoundingBox& bbox) + { + const auto dims = (DIM > 0 ? DIM : Base::dim_); + resize(bbox, dims); + + if (dataset_.kdtree_get_bbox(bbox)) + { + // Done! It was implemented in derived class + } + else + { + const Size N = Base::size_; + if (!N) + throw std::runtime_error( + "[nanoflann] computeBoundingBox() called but " + "no data points found."); + for (Dimension i = 0; i < dims; ++i) + { + bbox[i].low = bbox[i].high = + this->dataset_get(*this, Base::vAcc_[0], i); + } + for (Offset k = 1; k < N; ++k) + { + for (Dimension i = 0; i < dims; ++i) + { + const auto val = + this->dataset_get(*this, Base::vAcc_[k], i); + if (val < bbox[i].low) bbox[i].low = val; + if (val > bbox[i].high) bbox[i].high = val; + } + } + } + } + + /** + * Performs an exact search in the tree starting from a node. + * \tparam RESULTSET Should be any ResultSet + */ + template + void searchLevel( + RESULTSET& result_set, const ElementType* vec, const NodePtr node, + DistanceType mindist, distance_vector_t& dists, + const float epsError) const + { + /* If this is a leaf node, then do check and return. */ + if ((node->child1 == nullptr) && (node->child2 == nullptr)) + { + DistanceType worst_dist = result_set.worstDist(); + for (Offset i = node->node_type.lr.left; + i < node->node_type.lr.right; ++i) + { + const IndexType index = Base::vAcc_[i]; // reorder... : i; + if (treeIndex_[index] == -1) continue; + DistanceType dist = distance_.evalMetric( + vec, index, (DIM > 0 ? DIM : Base::dim_)); + if (dist < worst_dist) + { + if (!result_set.addPoint( + static_cast(dist), + static_cast( + Base::vAcc_[i]))) + { + // the resultset doesn't want to receive any more + // points, we're done searching! + return; // false; + } + } + } + return; + } + + /* Which child branch should be taken first? */ + Dimension idx = node->node_type.sub.divfeat; + ElementType val = vec[idx]; + DistanceType diff1 = val - node->node_type.sub.divlow; + DistanceType diff2 = val - node->node_type.sub.divhigh; + + NodePtr bestChild; + NodePtr otherChild; + DistanceType cut_dist; + if ((diff1 + diff2) < 0) + { + bestChild = node->child1; + otherChild = node->child2; + cut_dist = + distance_.accum_dist(val, node->node_type.sub.divhigh, idx); + } + else + { + bestChild = node->child2; + otherChild = node->child1; + cut_dist = + distance_.accum_dist(val, node->node_type.sub.divlow, idx); + } + + /* Call recursively to search next level down. */ + searchLevel(result_set, vec, bestChild, mindist, dists, epsError); + + DistanceType dst = dists[idx]; + mindist = mindist + cut_dist - dst; + dists[idx] = cut_dist; + if (mindist * epsError <= result_set.worstDist()) + { + searchLevel(result_set, vec, otherChild, mindist, dists, epsError); + } + dists[idx] = dst; + } + + public: + /** Stores the index in a binary file. + * IMPORTANT NOTE: The set of data points is NOT stored in the file, so + * when loading the index object it must be constructed associated to the + * same source of data points used while building it. See the example: + * examples/saveload_example.cpp \sa loadIndex */ + void saveIndex(std::ostream& stream) { saveIndex(*this, stream); } + + /** Loads a previous index from a binary file. + * IMPORTANT NOTE: The set of data points is NOT stored in the file, so + * the index object must be constructed associated to the same source of + * data points used while building the index. See the example: + * examples/saveload_example.cpp \sa loadIndex */ + void loadIndex(std::istream& stream) { loadIndex(*this, stream); } +}; + +/** kd-tree dynaimic index + * + * class to create multiple static index and merge their results to behave as + * single dynamic index as proposed in Logarithmic Approach. + * + * Example of usage: + * examples/dynamic_pointcloud_example.cpp + * + * \tparam DatasetAdaptor The user-provided adaptor (see comments above). + * \tparam Distance The distance metric to use: nanoflann::metric_L1, + * nanoflann::metric_L2, nanoflann::metric_L2_Simple, etc. \tparam DIM + * Dimensionality of data points (e.g. 3 for 3D points) \tparam IndexType + * Will be typically size_t or int + */ +template < + typename Distance, class DatasetAdaptor, int32_t DIM = -1, + typename IndexType = uint32_t> +class KDTreeSingleIndexDynamicAdaptor +{ + public: + using ElementType = typename Distance::ElementType; + using DistanceType = typename Distance::DistanceType; + + using Offset = typename KDTreeSingleIndexDynamicAdaptor_< + Distance, DatasetAdaptor, DIM>::Offset; + using Size = typename KDTreeSingleIndexDynamicAdaptor_< + Distance, DatasetAdaptor, DIM>::Size; + using Dimension = typename KDTreeSingleIndexDynamicAdaptor_< + Distance, DatasetAdaptor, DIM>::Dimension; + + protected: + Size leaf_max_size_; + Size treeCount_; + Size pointCount_; + + /** + * The dataset used by this index + */ + const DatasetAdaptor& dataset_; //!< The source of our data + + /** treeIndex[idx] is the index of tree in which point at idx is stored. + * treeIndex[idx]=-1 means that point has been removed. */ + std::vector treeIndex_; + std::unordered_set removedPoints_; + + KDTreeSingleIndexAdaptorParams index_params_; + + Dimension dim_; //!< Dimensionality of each data point + + using index_container_t = KDTreeSingleIndexDynamicAdaptor_< + Distance, DatasetAdaptor, DIM, IndexType>; + std::vector index_; + + public: + /** Get a const ref to the internal list of indices; the number of indices + * is adapted dynamically as the dataset grows in size. */ + const std::vector& getAllIndices() const + { + return index_; + } + + private: + /** finds position of least significant unset bit */ + int First0Bit(IndexType num) + { + int pos = 0; + while (num & 1) + { + num = num >> 1; + pos++; + } + return pos; + } + + /** Creates multiple empty trees to handle dynamic support */ + void init() + { + using my_kd_tree_t = KDTreeSingleIndexDynamicAdaptor_< + Distance, DatasetAdaptor, DIM, IndexType>; + std::vector index( + treeCount_, + my_kd_tree_t(dim_ /*dim*/, dataset_, treeIndex_, index_params_)); + index_ = index; + } + + public: + Distance distance_; + + /** + * KDTree constructor + * + * Refer to docs in README.md or online in + * https://github.com/jlblancoc/nanoflann + * + * The KD-Tree point dimension (the length of each point in the datase, e.g. + * 3 for 3D points) is determined by means of: + * - The \a DIM template parameter if >0 (highest priority) + * - Otherwise, the \a dimensionality parameter of this constructor. + * + * @param inputData Dataset with the input features. Its lifetime must be + * equal or longer than that of the instance of this class. + * @param params Basically, the maximum leaf node size + */ + explicit KDTreeSingleIndexDynamicAdaptor( + const int dimensionality, const DatasetAdaptor& inputData, + const KDTreeSingleIndexAdaptorParams& params = + KDTreeSingleIndexAdaptorParams(), + const size_t maximumPointCount = 1000000000U) + : dataset_(inputData), index_params_(params), distance_(inputData) + { + treeCount_ = static_cast(std::log2(maximumPointCount)) + 1; + pointCount_ = 0U; + dim_ = dimensionality; + treeIndex_.clear(); + if (DIM > 0) dim_ = DIM; + leaf_max_size_ = params.leaf_max_size; + init(); + const size_t num_initial_points = dataset_.kdtree_get_point_count(); + if (num_initial_points > 0) { addPoints(0, num_initial_points - 1); } + } + + /** Deleted copy constructor*/ + explicit KDTreeSingleIndexDynamicAdaptor( + const KDTreeSingleIndexDynamicAdaptor< + Distance, DatasetAdaptor, DIM, IndexType>&) = delete; + + /** Add points to the set, Inserts all points from [start, end] */ + void addPoints(IndexType start, IndexType end) + { + const Size count = end - start + 1; + int maxIndex = 0; + treeIndex_.resize(treeIndex_.size() + count); + for (IndexType idx = start; idx <= end; idx++) + { + const int pos = First0Bit(pointCount_); + maxIndex = std::max(pos, maxIndex); + treeIndex_[pointCount_] = pos; + + const auto it = removedPoints_.find(idx); + if (it != removedPoints_.end()) + { + removedPoints_.erase(it); + treeIndex_[idx] = pos; + } + + for (int i = 0; i < pos; i++) + { + for (int j = 0; j < static_cast(index_[i].vAcc_.size()); + j++) + { + index_[pos].vAcc_.push_back(index_[i].vAcc_[j]); + if (treeIndex_[index_[i].vAcc_[j]] != -1) + treeIndex_[index_[i].vAcc_[j]] = pos; + } + index_[i].vAcc_.clear(); + } + index_[pos].vAcc_.push_back(idx); + pointCount_++; + } + + for (int i = 0; i <= maxIndex; ++i) + { + index_[i].freeIndex(index_[i]); + if (!index_[i].vAcc_.empty()) index_[i].buildIndex(); + } + } + + /** Remove a point from the set (Lazy Deletion) */ + void removePoint(size_t idx) + { + if (idx >= pointCount_) return; + removedPoints_.insert(idx); + treeIndex_[idx] = -1; + } + + /** + * Find set of nearest neighbors to vec[0:dim-1]. Their indices are stored + * inside the result object. + * + * Params: + * result = the result object in which the indices of the + * nearest-neighbors are stored vec = the vector for which to search the + * nearest neighbors + * + * \tparam RESULTSET Should be any ResultSet + * \return True if the requested neighbors could be found. + * \sa knnSearch, radiusSearch + * + * \note If L2 norms are used, all returned distances are actually squared + * distances. + */ + template + bool findNeighbors( + RESULTSET& result, const ElementType* vec, + const SearchParameters& searchParams = {}) const + { + for (size_t i = 0; i < treeCount_; i++) + { + index_[i].findNeighbors(result, &vec[0], searchParams); + } + return result.full(); + } +}; + +/** An L2-metric KD-tree adaptor for working with data directly stored in an + * Eigen Matrix, without duplicating the data storage. You can select whether a + * row or column in the matrix represents a point in the state space. + * + * Example of usage: + * \code + * Eigen::Matrix mat; + * + * // Fill out "mat"... + * using my_kd_tree_t = nanoflann::KDTreeEigenMatrixAdaptor< + * Eigen::Matrix>; + * + * const int max_leaf = 10; + * my_kd_tree_t mat_index(mat, max_leaf); + * mat_index.index->... + * \endcode + * + * \tparam DIM If set to >0, it specifies a compile-time fixed dimensionality + * for the points in the data set, allowing more compiler optimizations. + * \tparam Distance The distance metric to use: nanoflann::metric_L1, + * nanoflann::metric_L2, nanoflann::metric_L2_Simple, etc. + * \tparam row_major If set to true the rows of the matrix are used as the + * points, if set to false the columns of the matrix are used as the + * points. + */ +template < + class MatrixType, int32_t DIM = -1, class Distance = nanoflann::metric_L2, + bool row_major = true> +struct KDTreeEigenMatrixAdaptor +{ + using self_t = + KDTreeEigenMatrixAdaptor; + using num_t = typename MatrixType::Scalar; + using IndexType = typename MatrixType::Index; + using metric_t = typename Distance::template traits< + num_t, self_t, IndexType>::distance_t; + + using index_t = KDTreeSingleIndexAdaptor< + metric_t, self_t, + row_major ? MatrixType::ColsAtCompileTime + : MatrixType::RowsAtCompileTime, + IndexType>; + + index_t* index_; //! The kd-tree index for the user to call its methods as + //! usual with any other FLANN index. + + using Offset = typename index_t::Offset; + using Size = typename index_t::Size; + using Dimension = typename index_t::Dimension; + + /// Constructor: takes a const ref to the matrix object with the data points + explicit KDTreeEigenMatrixAdaptor( + const Dimension dimensionality, + const std::reference_wrapper& mat, + const int leaf_max_size = 10, const unsigned int n_thread_build = 1) + : m_data_matrix(mat) + { + const auto dims = row_major ? mat.get().cols() : mat.get().rows(); + if (static_cast(dims) != dimensionality) + throw std::runtime_error( + "Error: 'dimensionality' must match column count in data " + "matrix"); + if (DIM > 0 && static_cast(dims) != DIM) + throw std::runtime_error( + "Data set dimensionality does not match the 'DIM' template " + "argument"); + index_ = new index_t( + dims, *this /* adaptor */, + nanoflann::KDTreeSingleIndexAdaptorParams( + leaf_max_size, nanoflann::KDTreeSingleIndexAdaptorFlags::None, + n_thread_build)); + } + + public: + /** Deleted copy constructor */ + KDTreeEigenMatrixAdaptor(const self_t&) = delete; + + ~KDTreeEigenMatrixAdaptor() { delete index_; } + + const std::reference_wrapper m_data_matrix; + + /** Query for the \a num_closest closest points to a given point (entered as + * query_point[0:dim-1]). Note that this is a short-cut method for + * index->findNeighbors(). The user can also call index->... methods as + * desired. + * + * \note If L2 norms are used, all returned distances are actually squared + * distances. + */ + void query( + const num_t* query_point, const Size num_closest, + IndexType* out_indices, num_t* out_distances) const + { + nanoflann::KNNResultSet resultSet(num_closest); + resultSet.init(out_indices, out_distances); + index_->findNeighbors(resultSet, query_point); + } + + /** @name Interface expected by KDTreeSingleIndexAdaptor + * @{ */ + + const self_t& derived() const { return *this; } + self_t& derived() { return *this; } + + // Must return the number of data points + Size kdtree_get_point_count() const + { + if (row_major) + return m_data_matrix.get().rows(); + else + return m_data_matrix.get().cols(); + } + + // Returns the dim'th component of the idx'th point in the class: + num_t kdtree_get_pt(const IndexType idx, size_t dim) const + { + if (row_major) + return m_data_matrix.get().coeff(idx, IndexType(dim)); + else + return m_data_matrix.get().coeff(IndexType(dim), idx); + } + + // Optional bounding-box computation: return false to default to a standard + // bbox computation loop. + // Return true if the BBOX was already computed by the class and returned + // in "bb" so it can be avoided to redo it again. Look at bb.size() to + // find out the expected dimensionality (e.g. 2 or 3 for point clouds) + template + bool kdtree_get_bbox(BBOX& /*bb*/) const + { + return false; + } + + /** @} */ + +}; // end of KDTreeEigenMatrixAdaptor +/** @} */ + +/** @} */ // end of grouping +} // namespace nanoflann diff --git a/python/2_headers.i b/python/2_headers.i index bc06d7d85..0b07700a4 100644 --- a/python/2_headers.i +++ b/python/2_headers.i @@ -41,6 +41,7 @@ %ignore operator crpropa::ObserverFeature*; %ignore operator crpropa::MagneticField*; %ignore operator crpropa::PhotonField*; +%ignore operator crpropa::InteractionRates*; %ignore operator crpropa::AdvectionField*; %ignore operator crpropa::ParticleCollector*; %ignore operator crpropa::Density*; @@ -270,6 +271,11 @@ %feature("director") crpropa::PhotonField; %include "crpropa/PhotonBackground.h" +%implicitconv crpropa::ref_ptr; +%template(InteractionRatesRefPtr) crpropa::ref_ptr; +%feature("director") crpropa::InteractionRates; +%include "crpropa/InteractionRates.h" + %implicitconv crpropa::ref_ptr; %template(AdvectionFieldRefPtr) crpropa::ref_ptr; %feature("director") crpropa::AdvectionField; diff --git a/src/Common.cpp b/src/Common.cpp index 591ededc8..03a77ae55 100644 --- a/src/Common.cpp +++ b/src/Common.cpp @@ -68,7 +68,7 @@ std::string getInstallPrefix() { std::string _path = ""; #ifdef CRPROPA_INSTALL_PREFIX - _path += CRPROPA_INSTALL_PREFIX; + _path += CRPROPA_INSTALL_PREFIX; #endif return _path; }; @@ -140,5 +140,11 @@ size_t closestIndex(double x, const std::vector &X) { return i1; } +std::string splitFilename(const std::string str) { + std::size_t found = str.find_last_of("/\\"); + std::string s = str.substr(found + 1); + return s; +} + } // namespace crpropa diff --git a/src/Geometry.cpp b/src/Geometry.cpp index 8e3ef1b57..6f1d25e9a 100644 --- a/src/Geometry.cpp +++ b/src/Geometry.cpp @@ -29,11 +29,14 @@ std::string Plane::getDescription() const { return ss.str(); }; +bool Plane::isInside(const Vector3d& point) const { + return distance(point) >= 0; // if the point is on the "positive" or "negative" side of the infinite plane +} + Vector3d Plane::normal(const Vector3d& point) const { return n; } - // Sphere ------------------------------------------------------------------ Sphere::Sphere(const Vector3d& _center, double _radius) : center(_center), radius(_radius) { }; @@ -48,6 +51,18 @@ Vector3d Sphere::normal(const Vector3d& point) const { return d.getUnitVector(); } +bool Sphere::isInside(const Vector3d& point) const { + Vector3d dR = point - center; + double dist2 = dR.getR2(); + + return dist2 <= radius * radius; + /**if (distanceSquared <= radius * radius) { + return true; + } else { + return false; + }*/ +} + std::string Sphere::getDescription() const { std::stringstream ss; ss << "Sphere: " << std::endl @@ -112,6 +127,10 @@ Vector3d ParaxialBox::normal(const Vector3d& point) const { return n; } +bool ParaxialBox::isInside(const Vector3d& point) const { + return distance(point) >= 0; +} + std::string ParaxialBox::getDescription() const { std::stringstream ss; ss << "ParaxialBox: " << std::endl diff --git a/src/InteractionRates.cpp b/src/InteractionRates.cpp new file mode 100644 index 000000000..acf71825d --- /dev/null +++ b/src/InteractionRates.cpp @@ -0,0 +1,507 @@ +#include "crpropa/InteractionRates.h" +#include "crpropa/Units.h" +#include "crpropa/Random.h" + +#include + +#include +#include +#include +#include +#include + +#if defined(__APPLE__) && defined(_LIBCPP_VERSION) + namespace fs = std::__fs::filesystem; +#else + namespace fs = std::filesystem; +#endif + +namespace crpropa { + +InteractionRatesHomogeneous::InteractionRatesHomogeneous(std::string RateFile, std::string CumulativeRateFile) { + + this->ratesName = "interactionRatesHomogeneous"; + this->isPositionDependent = false; + + if (RateFile!="") + initRate(RateFile); + if (CumulativeRateFile!="") + initCumulativeRate(CumulativeRateFile); +} + +std::vector InteractionRatesHomogeneous::getTabulatedEnergy() const { + return tabEnergy; +} + +std::vector InteractionRatesHomogeneous::getTabulatedRate() const { + return tabRate; +} + +std::vector InteractionRatesHomogeneous::getTabulatedE() const { + return tabE; +} + +std::vector InteractionRatesHomogeneous::getTabulateds() const { + return tabs; +} + +std::vector> InteractionRatesHomogeneous::getTabulatedCDF() const { + return tabCDF; +} + +double InteractionRatesHomogeneous::getProcessRate(const double E, const Vector3d &position) const { + if (!this->isPositionDependent) { + + // compute the interaction rate for the given candidate energy, E + double rate = interpolate(E, this->tabEnergy, this->tabRate); + return rate; + + } else { + + throw std::runtime_error("Error in boolean isPositionDependent!"); + + } +} + +void InteractionRatesHomogeneous::loadPerformInteractionTabs(const Vector3d &position, std::vector &tabE, std::vector &tabs, std::vector> &tabCDF) const { + tabE = this->getTabulatedE(); + tabs = this->getTabulateds(); + tabCDF = this->getTabulatedCDF(); +} + +void InteractionRatesHomogeneous::setTabulatedEnergy(std::vector& tabEnergy) { + this->tabEnergy = tabEnergy; +} + +void InteractionRatesHomogeneous::setTabulatedRate(std::vector& tabRate) { + this->tabRate = tabRate; +} + +void InteractionRatesHomogeneous::setTabulatedE(std::vector& tabE) { + this->tabE = tabE; +} + +void InteractionRatesHomogeneous::setTabulateds(std::vector& tabs) { + this->tabs = tabs; +} + +void InteractionRatesHomogeneous::setTabulatedCDF(std::vector>& tabCDF) { + this->tabCDF = tabCDF; +} + +void InteractionRatesHomogeneous::initRate(std::string filename){ + if(!fs::is_regular_file(filename)) + throw std::runtime_error("InteractionRatesHomogeneous: The given filename " + filename + " is a directory!\ + If you wanted to use position dependent photon fields instead, set the correct InteractionRates class first!"); + + std::ifstream infile(filename.c_str()); + + std::vector tabEnergy; + std::vector tabRate; + + if (!infile.good()) + throw std::runtime_error("InteractionRatesHomogeneous: could not open file " + filename); + + while (infile.good()) { + + if (infile.peek() != '#') { + double a, b; + infile >> a >> b; + if (infile) { + tabEnergy.push_back(pow(10, a) * eV); + tabRate.push_back(b / Mpc); + } + } + infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n'); + + } + + infile.close(); + + setTabulatedEnergy(tabEnergy); + setTabulatedRate(tabRate); +} + +void InteractionRatesHomogeneous::initCumulativeRate(std::string filename){ + if(!fs::is_regular_file(filename)) + throw std::runtime_error("InteractionRatesHomogeneous: The given filename " + filename + " is a directory!\ + If you wanted to use position dependent photon fields instead, set the correct InteractionRates class first!"); + + std::ifstream infile(filename.c_str()); + + std::vector tabE; + std::vector tabs; + std::vector> tabCDF; + + if (!infile.good()) + throw std::runtime_error("InteractionRatesHomogeneous: could not open file " + filename); + + // skip header + while (infile.peek() == '#') + infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n'); + + // read s values in first line + double a; + infile >> a; // skip first value + while (infile.good() and (infile.peek() != '\n')) { + infile >> a; + tabs.push_back(pow(10, a) * eV * eV); + } + + // read all following lines: E, cdf values + while (infile.good()) { + infile >> a; + if (!infile) + break; // end of file + tabE.push_back(pow(10, a) * eV); + std::vector cdf; + for (int i = 0; i < tabs.size(); i++) { + infile >> a; + cdf.push_back(a / Mpc); + } + tabCDF.push_back(cdf); + } + infile.close(); + + setTabulatedE(tabE); + setTabulateds(tabs); + setTabulatedCDF(tabCDF); +} + + +InteractionRatesPositionDependent::InteractionRatesPositionDependent( + std::string RateFilePath, std::string CumulativeRateFilePath, ref_ptr surface) { + + this->ratesName = "interactionRatesPositionDependent"; + this->isPositionDependent = true; + this->surface = surface; + + if (RateFilePath!="") + initRate(RateFilePath); + if (CumulativeRateFilePath!="") + initCumulativeRate(CumulativeRateFilePath); +} + +int InteractionRatesPositionDependent::findClosestGridPoint(const Vector3d &position) const { + + if (!tree) { + throw std::runtime_error("KD-Tree not initialized!"); + } + + unsigned int closestIndex; + double closestDistSquared; + double queryPoint[3] = { position.x, position.y, position.z }; + + this->tree->knnSearch(queryPoint, 1, &closestIndex, &closestDistSquared); + return this->cloud.ids[closestIndex]; + +} + +std::vector InteractionRatesPositionDependent::getTabulatedEnergy() const { + return tabEnergy; +} + +std::vector> InteractionRatesPositionDependent::getTabulatedRate() const { + return tabRate; +} + +std::vector InteractionRatesPositionDependent::getTabulatedE() const { + return tabE; +} + +std::vector> InteractionRatesPositionDependent::getTabulateds() const { + return tabs; +} + +std::vector>> InteractionRatesPositionDependent::getTabulatedCDF() const { + return tabCDF; +} + +std::unordered_map InteractionRatesPositionDependent::getPhotonDict() const { + return photonDict; +} + +std::vector InteractionRatesPositionDependent::getClosestRate(const Vector3d &position) const { + int iMin = findClosestGridPoint(position); + return tabRate[iMin]; +} + +std::vector InteractionRatesPositionDependent::getClosests(const Vector3d &position) const { + int iMin = findClosestGridPoint(position); + return tabs[iMin]; +} + +std::vector> InteractionRatesPositionDependent::getClosestCDF(const Vector3d &position) const { + int iMin = findClosestGridPoint(position); + return tabCDF[iMin]; +} + +double InteractionRatesPositionDependent::getProcessRate(const double E, const Vector3d &position) const { + if (!this->isPositionDependent) { + + throw std::runtime_error("Error in boolean isPositionDependent!"); + + } else { + + std::vector tabRate = this->getClosestRate(position); + + // compute the interaction rate for the given candidate energy, E + double rate = interpolate(E, this->tabEnergy, tabRate); + return rate; + + } +} + +void InteractionRatesPositionDependent::loadPerformInteractionTabs(const Vector3d &position, std::vector &tabE, std::vector &tabs, std::vector> &tabCDF) const { + + std::vector E = this->getTabulatedE(); + std::vector s = this->getClosests(position); + std::vector> CDF = this->getClosestCDF(position); + + tabE = E; + tabs = s; + tabCDF = CDF; +} + +void InteractionRatesPositionDependent::setTabulatedEnergy(std::vector& tabEnergy) { + this->tabEnergy = tabEnergy; +} + +void InteractionRatesPositionDependent::setTabulatedRate(std::vector>& tabRate) { + this->tabRate = tabRate; +} + +void InteractionRatesPositionDependent::setTabulatedE(std::vector& tabE) { + this->tabE = tabE; +} + +void InteractionRatesPositionDependent::setTabulateds(std::vector>& tabs) { + this->tabs = tabs; +} + +void InteractionRatesPositionDependent::setTabulatedCDF(std::vector>>& tabCDF) { + this->tabCDF = tabCDF; +} + +void InteractionRatesPositionDependent::setPhotonDict(std::unordered_map& photonDict) { + this->photonDict = photonDict; + + // delete old clouds + this->cloud.points.clear(); + this->cloud.ids.clear(); + + for (const auto& el : this->photonDict) { + this->cloud.ids.push_back(el.first); + this->cloud.points.push_back(el.second); + } + + // delete old tree + if (this->tree) { + delete this->tree; + } + + int maxLeafTree = 20; + int nThreads = 4; + nanoflann::KDTreeSingleIndexAdaptorFlags flag; + + this->tree = new KDTree(3, cloud, nanoflann::KDTreeSingleIndexAdaptorParams(maxLeafTree, flag, nThreads)); + this->tree->buildIndex(); + +} + +void InteractionRatesPositionDependent::setSurface(ref_ptr surface) { + this->surface = surface; +} + +ref_ptr InteractionRatesPositionDependent::getSurface() const { + return this->surface; +} + +void InteractionRatesPositionDependent::initRate(std::string filepath){ + if(!fs::is_directory(filepath)) + throw std::runtime_error("InteractionRatesPositionDependent: The given path " + filepath + " is not a directory!\ + If you wanted to use homogeneous photon fields instead, set the correct InteractionRates class first!"); + + std::vector> tabRate; + + fs::path dir = filepath; + std::unordered_map photonDict; + int iFile = 0; + + if (!fs::exists(dir)) { + std::cout << "Photon tables not found in " << dir << std::endl; + return; + } + + for (auto const& dir_entry : fs::directory_iterator{dir}) { + + // the input filename here should be a string + //check if it is correct, i.e. a proper filename string + std::string filename = dir_entry.path().string(); + std::ifstream infile(filename.c_str()); + + std::vector vecEnergy; + std::vector vecRate; + + if (!infile.good()) + throw + std::runtime_error("InteractionRatesPositionDependent: could not open file " + filename); + + double x, y, z; + std::string str; + std::stringstream ss; + + std::string filename_split = splitFilename(dir_entry.path().string()); + ss << filename_split; + + int iLine = 0; + + std::locale::global(std::locale("C")); + + while (getline(ss, str, '_')) { + if (iLine == 3) { + x = -std::stod(str) * kpc; + } + if (iLine == 4) { + y = std::stod(str) * kpc; + } + if (iLine == 5) { + z = std::stod(str) * kpc; + } + iLine = iLine + 1; + } + + Vector3d vPos(x, y, z); + + if (getSurface() and !getSurface()->isInside(vPos)) + continue; + + photonDict[iFile] = vPos; + + while (infile.good()) { + if (infile.peek() != '#') { + double a, b; + infile >> a >> b; + if (infile) { + if (iFile == 0) { + vecEnergy.push_back(pow(10, a) * eV); + setTabulatedEnergy(vecEnergy); + } + vecRate.push_back(b / Mpc); + } + } + infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n'); + } + + tabRate.push_back(vecRate); + + iFile = iFile + 1; + infile.close(); + + } + + if (tabRate.empty()) + throw std::runtime_error("Rate's table empty! Check if the surface is properly set."); + + setTabulatedRate(tabRate); + setPhotonDict(photonDict); +} + +void InteractionRatesPositionDependent::initCumulativeRate(std::string filepath){ + if(!fs::is_directory(filepath)) + throw std::runtime_error("InteractionRatesPositionDependent: The given path " + filepath + " is not a directory!\ + If you wanted to use homogeneous photon fields instead, set the correct InteractionRates class first!"); + + std::vector> tabs; + std::vector>> tabCDF; + + fs::path dir = filepath; + int iFile = 0; + + if (!fs::exists(dir)) { + std::cout << "Photon tables not found in " << dir << std::endl; + return; + } + + for (auto const& dir_entry : fs::directory_iterator{dir}) { + + std::vector vecE; + std::vector vecs; + std::vector> vecCDF; + + std::string filename = dir_entry.path().string(); + std::ifstream infile(filename.c_str()); + + if (!infile.good()) + throw std::runtime_error("InteractionRatesPositionDependent: could not open file " + filename); + + double x, y, z; + std::string str; + std::stringstream ss; + + std::string filename_split = splitFilename(dir_entry.path().string()); + ss << filename_split; + + int iLine = 0; + + std::locale::global(std::locale("C")); + + while (getline(ss, str, '_')) { + if (iLine == 3) { + x = -std::stod(str) * kpc; + } + if (iLine == 4) { + y = std::stod(str) * kpc; + } + if (iLine == 5) { + z = std::stod(str) * kpc; + } + iLine = iLine + 1; + } + + Vector3d vPos(x, y, z); + + if (getSurface() and !getSurface()->isInside(vPos)) + continue; + + // skip header + while (infile.peek() == '#') + infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n'); + + // read s values in first line + double a; + infile >> a; // skip first value + while (infile.good() and (infile.peek() != '\n')) { + infile >> a; + vecs.push_back(pow(10, a) * eV * eV); + } + + // read all following lines: E, cdf values + while (infile.good()) { + infile >> a; + if (!infile) + break; // end of file + if (iFile == 0) { + vecE.push_back(pow(10, a) * eV); + setTabulatedE(vecE); + } + std::vector cdf; + for (int i = 0; i < tabs.size(); i++) { + infile >> a; + cdf.push_back(a / Mpc); + } + vecCDF.push_back(cdf); + } + + iFile = iFile + 1; + + tabs.push_back(vecs); + tabCDF.push_back(vecCDF); + infile.close(); + } + + setTabulateds(tabs); + setTabulatedCDF(tabCDF); +} + +} //namespace crpropa diff --git a/src/PhotonBackground.cpp b/src/PhotonBackground.cpp index b29ffd7da..08a29492f 100644 --- a/src/PhotonBackground.cpp +++ b/src/PhotonBackground.cpp @@ -1,52 +1,79 @@ #include "crpropa/PhotonBackground.h" #include "crpropa/Units.h" #include "crpropa/Random.h" +#include "crpropa/Common.h" #include "kiss/logger.h" #include +#include #include #include #include +#include +#include +#include + +#if defined(__APPLE__) && defined(_LIBCPP_VERSION) + #include + namespace fs = std::__fs::filesystem; +#elif defined(__GNUC__) && (__GNUC__ < 10) + #include + namespace fs = std::experimental::filesystem; +#else + #include + namespace fs = std::filesystem; +#endif namespace crpropa { TabularPhotonField::TabularPhotonField(std::string fieldName, bool isRedshiftDependent) { this->fieldName = fieldName; this->isRedshiftDependent = isRedshiftDependent; - + this->isPositionDependent = false; + this->surface = NULL; + readPhotonEnergy(getDataPath("") + "Scaling/" + this->fieldName + "_photonEnergy.txt"); readPhotonDensity(getDataPath("") + "Scaling/" + this->fieldName + "_photonDensity.txt"); + if (this->isRedshiftDependent) readRedshift(getDataPath("") + "Scaling/" + this->fieldName + "_redshift.txt"); - + checkInputData(); - + if (this->isRedshiftDependent) initRedshiftScaling(); + } -double TabularPhotonField::getPhotonDensity(double Ephoton, double z) const { - if ((this->isRedshiftDependent)) { +double TabularPhotonField::getPhotonDensity(double Ephoton, double z, const Vector3d &pos) const { + if (this->isRedshiftDependent) { // fix behaviour for future redshift. See issue #414 - // with redshift < 0 the photon density is set to 0 in interpolate2d. + // with redshift < 0 the photon density is set to 0 in interpolate2d. // Therefore it is assumed that the photon density does not change from values at z = 0. This is only valid for small changes in redshift. double zMin = this->redshifts[0]; + if(z < zMin){ + if(z < -1) { KISS_LOG_WARNING << "Photon Field " << fieldName << " uses FutureRedshift with z < -1. The photon density is set to n(Ephoton, z=0). \n"; + } return getPhotonDensity(Ephoton, zMin); + } else { + return interpolate2d(Ephoton, z, this->photonEnergies, this->redshifts, this->photonDensity); + } + } else { return interpolate(Ephoton, this->photonEnergies, this->photonDensity); + } } - double TabularPhotonField::getRedshiftScaling(double z) const { if (!this->isRedshiftDependent) return 1.; @@ -60,29 +87,34 @@ double TabularPhotonField::getRedshiftScaling(double z) const { return interpolate(z, this->redshifts, this->redshiftScalings); } -double TabularPhotonField::getMinimumPhotonEnergy(double z) const{ +double TabularPhotonField::getMinimumPhotonEnergy(double z, const Vector3d &pos) const{ return photonEnergies[0]; } -double TabularPhotonField::getMaximumPhotonEnergy(double z) const{ +double TabularPhotonField::getMaximumPhotonEnergy(double z, const Vector3d &pos) const{ return photonEnergies[photonEnergies.size() -1]; } void TabularPhotonField::readPhotonEnergy(std::string filePath) { + std::ifstream infile(filePath.c_str()); + if (!infile.good()) throw std::runtime_error("TabularPhotonField::readPhotonEnergy: could not open " + filePath); - + std::string line; while (std::getline(infile, line)) { if ((line.size() > 0) & (line[0] != '#') ) this->photonEnergies.push_back(std::stod(line)); } + infile.close(); + } void TabularPhotonField::readPhotonDensity(std::string filePath) { std::ifstream infile(filePath.c_str()); + if (!infile.good()) throw std::runtime_error("TabularPhotonField::readPhotonDensity: could not open " + filePath); @@ -91,27 +123,34 @@ void TabularPhotonField::readPhotonDensity(std::string filePath) { if ((line.size() > 0) & (line[0] != '#') ) this->photonDensity.push_back(std::stod(line)); } + infile.close(); } void TabularPhotonField::readRedshift(std::string filePath) { std::ifstream infile(filePath.c_str()); + if (!infile.good()) throw std::runtime_error("TabularPhotonField::initRedshift: could not open " + filePath); - + std::string line; while (std::getline(infile, line)) { if ((line.size() > 0) & (line[0] != '#') ) this->redshifts.push_back(std::stod(line)); } + infile.close(); + } void TabularPhotonField::initRedshiftScaling() { double n0 = 0.; + for (int i = 0; i < this->redshifts.size(); ++i) { + double z = this->redshifts[i]; double n = 0.; + for (int j = 0; j < this->photonEnergies.size()-1; ++j) { double e_j = this->photonEnergies[j]; double e_j1 = this->photonEnergies[j+1]; @@ -120,22 +159,30 @@ void TabularPhotonField::initRedshiftScaling() { n0 += (getPhotonDensity(e_j, 0) + getPhotonDensity(e_j1, 0)) / 2. * deltaLogE; n += (getPhotonDensity(e_j, z) + getPhotonDensity(e_j1, z)) / 2. * deltaLogE; } + this->redshiftScalings.push_back(n / n0); + } } void TabularPhotonField::checkInputData() const { if (this->isRedshiftDependent) { - if (this->photonDensity.size() != this->photonEnergies.size() * this-> redshifts.size()) + + if (this->photonDensity.size() != this->photonEnergies.size() * this-> redshifts.size()) throw std::runtime_error("TabularPhotonField::checkInputData: length of photon density input is unequal to length of photon energy input times length of redshift input"); + } else { + if (this->photonEnergies.size() != this->photonDensity.size()) throw std::runtime_error("TabularPhotonField::checkInputData: length of photon energy input is unequal to length of photon density input"); + } for (int i = 0; i < this->photonEnergies.size(); ++i) { + double ePrevious = 0.; double e = this->photonEnergies[i]; + if (e <= 0.) throw std::runtime_error("TabularPhotonField::checkInputData: a value in the photon energy input is not positive"); if (e <= ePrevious) @@ -144,17 +191,21 @@ void TabularPhotonField::checkInputData() const { } for (int i = 0; i < this->photonDensity.size(); ++i) { + if (this->photonDensity[i] < 0.) throw std::runtime_error("TabularPhotonField::checkInputData: a value in the photon density input is negative"); + } if (this->isRedshiftDependent) { + if (this->redshifts[0] != 0.) throw std::runtime_error("TabularPhotonField::checkInputData: redshift input must start with zero"); for (int i = 0; i < this->redshifts.size(); ++i) { double zPrevious = -1.; double z = this->redshifts[i]; + if (z < 0.) throw std::runtime_error("TabularPhotonField::checkInputData: a value in the redshift input is negative"); if (z <= zPrevious) @@ -163,24 +214,234 @@ void TabularPhotonField::checkInputData() const { } for (int i = 0; i < this->redshiftScalings.size(); ++i) { + double scalingFactor = this->redshiftScalings[i]; if (scalingFactor <= 0.) throw std::runtime_error("TabularPhotonField::checkInputData: initRedshiftScaling has created a non-positive scaling factor"); + } } } -BlackbodyPhotonField::BlackbodyPhotonField(std::string fieldName, double blackbodyTemperature) { +TabularSpatialPhotonField::TabularSpatialPhotonField(std::string fieldName, ref_ptr surface) { + this->fieldName = fieldName; - this->blackbodyTemperature = blackbodyTemperature; - this->quantile = 0.0001; // tested to be sufficient, only used for extreme values of primary energy or temperature + this->isRedshiftDependent = isRedshiftDependent; + this->isPositionDependent = true; + this->surface = surface; + + fs::path dirE = getDataPath("") + "Scaling/" + this->fieldName + "/photonEnergy/"; + if (!fs::exists(dirE)) { + std::cout << "Photon tables not found in " << dirE << std::endl; + return; + } + + std::unordered_map photonDict; + int iFile = 0; + + for (auto const& dir_entry : fs::directory_iterator{dirE}) { + std::vector vE = readPhotonEnergy(dir_entry.path().string()); + + this->photonEnergies = vE; + break; + + } + + fs::path dirD = getDataPath("") + "Scaling/" + this->fieldName + "/photonDensity/"; + + if (!fs::exists(dirD)) { + std::cout << "Photon tables not found in " << dirD << std::endl; + return; + } + + for (auto const& dir_entry : fs::directory_iterator{dirD}) { + + double x, y, z; + std::string str; + std::stringstream ss; + + + std::string filename = splitFilename(dir_entry.path().string()); + ss << filename; + + //Getline function to take and store the x, y, z coordinates of each node + int iLine = 0; + // it ensures the double numbers are of the type 1.00329, with the . for the decimal part + std::locale::global(std::locale("C")); + + while (getline(ss, str, '_')) { + if (iLine == 2) { + x = -stod(str) * kpc; + } + if (iLine == 3) { + y = stod(str) * kpc; + } + if (iLine == 4) { + z = stod(str) * kpc; + } + iLine = iLine + 1; + } + + Vector3d vPos(x, y, z); + + if (getSurface() and !getSurface()->isInside(vPos)) + continue; + + photonDict[iFile] = vPos; + + iFile = iFile + 1; + + std::vector vD = readPhotonDensity(dir_entry.path().string()); + this->photonDensity.push_back(vD); + + } + + if (this->photonDensity.empty()) + throw std::runtime_error("Tabular spatial photon field for " + fieldName + " empty! Check if the surface is properly set."); + + this->photonDict = photonDict; + checkInputData(); + +} + +void TabularSpatialPhotonField::setSurface(ref_ptr surface) { + this->surface = surface; +} + +double TabularSpatialPhotonField::getPhotonDensity(const double ePhoton, double z, const Vector3d &pos) const { + + double dMin = 1000.; + int iMin = -1; + + for (const auto& el : this->photonDict) { + + Vector3d posNode = el.second; + double d; + d = sqrt((- posNode.x / kpc - pos.x / kpc) * (- posNode.x / kpc - pos.x / kpc) + (posNode.y / kpc - pos.y / kpc) * (posNode.y / kpc - pos.y / kpc ) + (posNode.z / kpc - pos.z / kpc) * (posNode.z / kpc - pos.z / kpc)); + + if (d photonEnergies[photonEnergies.size() - 1])) { + return 0; + + } else { + + std::vector rowE = this->photonEnergies; // assuming all the nodes have the same energy binning + std::vector rowD = this->photonDensity[iMin]; + return interpolate(ePhoton, rowE, rowD); + + } + } +} + +double TabularSpatialPhotonField::getMinimumPhotonEnergy(double z, const Vector3d &pos) const { + return photonEnergies[0]; // assuming all the nodes have the same energy bins +} + +double TabularSpatialPhotonField::getMaximumPhotonEnergy(double z, const Vector3d &pos) const { + return photonEnergies[photonEnergies.size() - 1]; // assuming all the nodes have the same energy bins +} + +std::vector TabularSpatialPhotonField::readPhotonEnergy(std::string filePath) { + std::ifstream infile(filePath.c_str()); + + if (!infile.good()) + throw std::runtime_error("TabularPhotonField::readPhotonEnergy: could not open " + filePath); + + + std::string line; + std::vector vE; + + while (std::getline(infile, line)) { + + if ((line.size() > 0) & (line[0] != '#') ) { + vE.insert(vE.begin(),std::stod(line)); + + } + } + + infile.close(); + return vE; + +} + +std::vector TabularSpatialPhotonField::readPhotonDensity(std::string filePath) { + + std::ifstream infile(filePath.c_str()); + + if (!infile.good()) + throw std::runtime_error("TabularPhotonField::readPhotonDensity: could not open " + filePath); + + + std::string line; + std::vector vD; + + + while (std::getline(infile, line)) { + + if ((line.size() > 0) & (line[0] != '#') ) + vD.insert(vD.begin(),std::stod(line)); + + } + + infile.close(); + return vD; + +} + +void TabularSpatialPhotonField::checkInputData() const { + + std::size_t numRowsDens = this->photonEnergies.size(); + std::size_t numRowsEn = this->photonDensity.size(); + + for (int j = 0; j < this->photonDensity.size(); ++j) { //take the proper row size! + + if (this->photonEnergies.size() != this->photonDensity[j].size()) + throw std::runtime_error("TabularPhotonField::checkInputData: length of photon energy input is unequal to length of photon density input"); + + for (int i = 0; i < this->photonEnergies.size(); ++i) { + + double ePrevious = 0.; + double e = this->photonEnergies[i]; + + if (e <= 0.) + throw std::runtime_error("TabularSpatialPhotonField::checkInputData: a value in the photon energy input is not positive"); + + if (e <= ePrevious) + throw std::runtime_error("TabularSpatialPhotonField::checkInputData: photon energy values are not strictly increasing"); + + ePrevious = e; + + } + + for (int i = 0; i < this->photonDensity[j].size(); ++i) { + + if (this->photonDensity[j][i] < 0.) + throw std::runtime_error("TabularSpatialPhotonField::checkInputData: a value in the photon density input is negative"); + + } + } +} + +BlackbodyPhotonField::BlackbodyPhotonField(std::string fieldName, double blackbodyTemperature) { + this->fieldName = fieldName; + this->blackbodyTemperature = blackbodyTemperature; + this->quantile = 0.0001; // tested to be sufficient, only used for extreme values of primary energy or temperature } -double BlackbodyPhotonField::getPhotonDensity(double Ephoton, double z) const { +double BlackbodyPhotonField::getPhotonDensity(double Ephoton, double z, const Vector3d &pos) const { return 8 * M_PI * pow_integer<3>(Ephoton / (h_planck * c_light)) / std::expm1(Ephoton / (k_boltzmann * this->blackbodyTemperature)); } -double BlackbodyPhotonField::getMinimumPhotonEnergy(double z) const { +double BlackbodyPhotonField::getMinimumPhotonEnergy(double z, const Vector3d &pos) const { double A; int quantile_int = 10000 * quantile; switch (quantile_int) @@ -201,7 +462,7 @@ double BlackbodyPhotonField::getMinimumPhotonEnergy(double z) const { return A * this -> blackbodyTemperature; } -double BlackbodyPhotonField::getMaximumPhotonEnergy(double z) const { +double BlackbodyPhotonField::getMaximumPhotonEnergy(double z, const Vector3d &pos) const { double factor = std::max(1., blackbodyTemperature / 2.73); return 0.1 * factor * eV; // T dependent scaling, starting at 0.1 eV as suitable for CMB } diff --git a/src/module/EMDoublePairProduction.cpp b/src/module/EMDoublePairProduction.cpp index 7477cf81a..8c028c9d6 100644 --- a/src/module/EMDoublePairProduction.cpp +++ b/src/module/EMDoublePairProduction.cpp @@ -1,25 +1,57 @@ #include "crpropa/module/EMDoublePairProduction.h" #include "crpropa/Units.h" #include "crpropa/Random.h" +#include "crpropa/Common.h" #include +#include #include #include +#include +#include +#include +#include +#include + +#if defined(__APPLE__) && defined(_LIBCPP_VERSION) + namespace fs = std::__fs::filesystem; +#else + namespace fs = std::filesystem; +#endif namespace crpropa { -EMDoublePairProduction::EMDoublePairProduction(ref_ptr photonField, bool haveElectrons, double thinning, double limit) { +EMDoublePairProduction::EMDoublePairProduction(ref_ptr photonField, bool haveElectrons, double thinning, double limit, ref_ptr surface) { + + setSurface(surface); setPhotonField(photonField); setHaveElectrons(haveElectrons); setLimit(limit); setThinning(thinning); + } void EMDoublePairProduction::setPhotonField(ref_ptr photonField) { + this->photonField = photonField; std::string fname = photonField->getFieldName(); setDescription("EMDoublePairProduction: " + fname); - initRate(getDataPath("EMDoublePairProduction/rate_" + fname + ".txt")); + + // choose the right interaction rates for the used photon field + if (!this->photonField->hasPositionDependence()) { + this->interactionRates = new InteractionRatesHomogeneous( + getDataPath("EMDoublePairProduction/rate_" + fname + ".txt") + ); + + } else { + this->interactionRates = new InteractionRatesPositionDependent( + getDataPath("EMDoublePairProduction/"+fname+"/Rate/"), + "", + this->surface + ); + + } + } void EMDoublePairProduction::setHaveElectrons(bool haveElectrons) { @@ -34,37 +66,27 @@ void EMDoublePairProduction::setThinning(double thinning) { this->thinning = thinning; } -void EMDoublePairProduction::initRate(std::string filename) { - std::ifstream infile(filename.c_str()); +void EMDoublePairProduction::setSurface(ref_ptr surface) { + this->surface = surface; +} - if (!infile.good()) - throw std::runtime_error("EMDoublePairProduction: could not open file " + filename); +ref_ptr EMDoublePairProduction::getSurface() const { + return this->surface; +} - // clear previously loaded interaction rates - tabEnergy.clear(); - tabRate.clear(); +void EMDoublePairProduction::setInteractionRates(ref_ptr intRates) { + this->interactionRates = intRates; +} - while (infile.good()) { - if (infile.peek() != '#') { - double a, b; - infile >> a >> b; - if (infile) { - tabEnergy.push_back(pow(10, a) * eV); - tabRate.push_back(b / Mpc); - } - } - infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n'); - } - infile.close(); +ref_ptr EMDoublePairProduction::getInteractionRates() const { + return this->interactionRates; } +void EMDoublePairProduction::initRate(std::string path) { + this->interactionRates->initRate(path); +} void EMDoublePairProduction::performInteraction(Candidate *candidate) const { - // the photon is lost after the interaction - candidate->setActive(false); - - if (not haveElectrons) - return; // Use assumption of Lee 96 arXiv:9604098 // Energy is equally shared between one e+e- pair, but take mass of second e+e- pair into account. @@ -73,6 +95,12 @@ void EMDoublePairProduction::performInteraction(Candidate *candidate) const { double E = candidate->current.getEnergy() * (1 + z); double Ee = (E - 2 * mass_electron * c_squared) / 2; + // the photon is lost after the interaction + candidate->setActive(false); + + if (not haveElectrons) + return; + Random &random = Random::instance(); Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(), candidate->current.getPosition()); @@ -86,9 +114,11 @@ void EMDoublePairProduction::performInteraction(Candidate *candidate) const { double w = 1. / pow(f, thinning); candidate->addSecondary(-11, Ee / (1 + z), pos, w, interactionTag); } + } void EMDoublePairProduction::process(Candidate *candidate) const { + // check if photon if (candidate->current.getId() != 22) return; @@ -96,13 +126,14 @@ void EMDoublePairProduction::process(Candidate *candidate) const { // scale the electron energy instead of background photons double z = candidate->getRedshift(); double E = (1 + z) * candidate->current.getEnergy(); - - // check if in tabulated energy range - if (E < tabEnergy.front() or (E > tabEnergy.back())) - return; + Vector3d position = candidate->current.getPosition(); // interaction rate - double rate = interpolate(E, tabEnergy, tabRate); + double rate = this->interactionRates->getProcessRate(E, position); + + if (rate < 0) + return; + rate *= pow_integer<2>(1 + z) * photonField->getRedshiftScaling(z); // check for interaction diff --git a/src/module/EMInverseComptonScattering.cpp b/src/module/EMInverseComptonScattering.cpp index c15d689d6..e85515f8a 100644 --- a/src/module/EMInverseComptonScattering.cpp +++ b/src/module/EMInverseComptonScattering.cpp @@ -3,27 +3,58 @@ #include "crpropa/Random.h" #include "crpropa/Common.h" +#include "crpropa/InteractionRates.h" + #include +#include +#include // Required for std::fixed and std::setprecision #include #include +#include +#include +#include +#include +#include + +#if defined(__APPLE__) && defined(_LIBCPP_VERSION) + namespace fs = std::__fs::filesystem; +#else + namespace fs = std::filesystem; +#endif namespace crpropa { static const double mec2 = mass_electron * c_squared; -EMInverseComptonScattering::EMInverseComptonScattering(ref_ptr photonField, bool havePhotons, double thinning, double limit) { +EMInverseComptonScattering::EMInverseComptonScattering(ref_ptr photonField, bool havePhotons, double thinning, double limit, ref_ptr surface) { + + setSurface(surface); setPhotonField(photonField); setHavePhotons(havePhotons); setLimit(limit); setThinning(thinning); + } void EMInverseComptonScattering::setPhotonField(ref_ptr photonField) { + this->photonField = photonField; std::string fname = photonField->getFieldName(); setDescription("EMInverseComptonScattering: " + fname); - initRate(getDataPath("EMInverseComptonScattering/rate_" + fname + ".txt")); - initCumulativeRate(getDataPath("EMInverseComptonScattering/cdf_" + fname + ".txt")); + + if (!this->photonField->hasPositionDependence()) { + this->interactionRates = new InteractionRatesHomogeneous( + getDataPath("EMInverseComptonScattering/rate_" + fname + ".txt"), + getDataPath("EMInverseComptonScattering/cdf_" + fname + ".txt") + ); + + } else { + this->interactionRates = new InteractionRatesPositionDependent( + getDataPath("EMInverseComptonScattering/" + fname + "/Rate/"), + getDataPath("EMInverseComptonScattering/" + fname + "/CumulativeRate/"), + this->surface + ); + } } void EMInverseComptonScattering::setHavePhotons(bool havePhotons) { @@ -38,67 +69,28 @@ void EMInverseComptonScattering::setThinning(double thinning) { this->thinning = thinning; } -void EMInverseComptonScattering::initRate(std::string filename) { - std::ifstream infile(filename.c_str()); - - if (!infile.good()) - throw std::runtime_error("EMInverseComptonScattering: could not open file " + filename); - - // clear previously loaded tables - tabEnergy.clear(); - tabRate.clear(); - - while (infile.good()) { - if (infile.peek() != '#') { - double a, b; - infile >> a >> b; - if (infile) { - tabEnergy.push_back(pow(10, a) * eV); - tabRate.push_back(b / Mpc); - } - } - infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n'); - } - infile.close(); +void EMInverseComptonScattering::setSurface(ref_ptr surface) { + this->surface = surface; } -void EMInverseComptonScattering::initCumulativeRate(std::string filename) { - std::ifstream infile(filename.c_str()); +ref_ptr EMInverseComptonScattering::getSurface() const { + return this->surface; +} - if (!infile.good()) - throw std::runtime_error("EMInverseComptonScattering: could not open file " + filename); +void EMInverseComptonScattering::setInteractionRates(ref_ptr intRates) { + this->interactionRates = intRates; +} - // clear previously loaded tables - tabE.clear(); - tabs.clear(); - tabCDF.clear(); - - // skip header - while (infile.peek() == '#') - infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n'); +ref_ptr EMInverseComptonScattering::getInteractionRates() const { + return this->interactionRates; +} - // read s values in first line - double a; - infile >> a; // skip first value - while (infile.good() and (infile.peek() != '\n')) { - infile >> a; - tabs.push_back(pow(10, a) * eV * eV); - } +void EMInverseComptonScattering::initRate(std::string path) { + this->interactionRates->initRate(path); +} - // read all following lines: E, cdf values - while (infile.good()) { - infile >> a; - if (!infile) - break; // end of file - tabE.push_back(pow(10, a) * eV); - std::vector cdf; - for (int i = 0; i < tabs.size(); i++) { - infile >> a; - cdf.push_back(a / Mpc); - } - tabCDF.push_back(cdf); - } - infile.close(); +void EMInverseComptonScattering::initCumulativeRate(std::string path) { + this->interactionRates->initCumulativeRate(path); } // Class to calculate the energy distribution of the ICS photon and to sample from it @@ -170,37 +162,47 @@ class ICSSecondariesEnergyDistribution { }; void EMInverseComptonScattering::performInteraction(Candidate *candidate) const { + // scale the particle energy instead of background photons double z = candidate->getRedshift(); double E = candidate->current.getEnergy() * (1 + z); - + + Vector3d position = candidate->current.getPosition(); + + std::vector tabE; + std::vector tabs; + std::vector> tabCDF; + + this->interactionRates->loadPerformInteractionTabs(position, tabE, tabs, tabCDF); + if (E < tabE.front() or E > tabE.back()) return; - + // sample the value of s Random &random = Random::instance(); size_t i = closestIndex(E, tabE); size_t j = random.randBin(tabCDF[i]); double s_kin = pow(10, log10(tabs[j]) + (random.rand() - 0.5) * 0.1); double s = s_kin + mec2 * mec2; - + // sample electron energy after scattering static ICSSecondariesEnergyDistribution distribution; double Enew = distribution.sample(E, s); - + // add up-scattered photon + double Esecondary = E - Enew; + double f = Enew / E; if (havePhotons) { - double Esecondary = E - Enew; - double f = Enew / E; if (random.rand() < pow(1 - f, thinning)) { double w = 1. / pow(1 - f, thinning); Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(), candidate->current.getPosition()); candidate->addSecondary(22, Esecondary / (1 + z), pos, w, interactionTag); } } - + // update the primary particle energy; do this after adding the secondary to correctly set the secondary's parent candidate->current.setEnergy(Enew / (1 + z)); + } void EMInverseComptonScattering::process(Candidate *candidate) const { @@ -208,34 +210,37 @@ void EMInverseComptonScattering::process(Candidate *candidate) const { int id = candidate->current.getId(); if (abs(id) != 11) return; - + // scale the particle energy instead of background photons double z = candidate->getRedshift(); double E = candidate->current.getEnergy() * (1 + z); - - if (E < tabEnergy.front() or (E > tabEnergy.back())) - return; - + Vector3d position = candidate->current.getPosition(); + // interaction rate - double rate = interpolate(E, tabEnergy, tabRate); + double rate = this->interactionRates->getProcessRate(E, position); + + if (rate < 0) + return; + rate *= pow_integer<2>(1 + z) * photonField->getRedshiftScaling(z); - + // run this loop at least once to limit the step size double step = candidate->getCurrentStep(); Random &random = Random::instance(); do { double randDistance = -log(random.rand()) / rate; - + // check for interaction; if it doesn't ocurr, limit next step if (step < randDistance) { candidate->limitNextStep(limit / rate); return; } performInteraction(candidate); - + // repeat with remaining step step -= randDistance; } while (step > 0); + } void EMInverseComptonScattering::setInteractionTag(std::string tag) { diff --git a/src/module/EMPairProduction.cpp b/src/module/EMPairProduction.cpp index 717c62cc4..454532eac 100644 --- a/src/module/EMPairProduction.cpp +++ b/src/module/EMPairProduction.cpp @@ -4,27 +4,55 @@ #include "kiss/logger.h" #include +#include +#include #include #include - +#include +#include +#include +#include +#include + +#if defined(__APPLE__) && defined(_LIBCPP_VERSION) + namespace fs = std::__fs::filesystem; +#else + namespace fs = std::filesystem; +#endif namespace crpropa { static const double mec2 = mass_electron * c_squared; -EMPairProduction::EMPairProduction(ref_ptr photonField, bool haveElectrons, double thinning, double limit) { +EMPairProduction::EMPairProduction(ref_ptr photonField, bool haveElectrons, double thinning, double limit, ref_ptr surface) { + + setSurface(surface); setPhotonField(photonField); setThinning(thinning); setLimit(limit); setHaveElectrons(haveElectrons); + } void EMPairProduction::setPhotonField(ref_ptr photonField) { + this->photonField = photonField; std::string fname = photonField->getFieldName(); setDescription("EMPairProduction: " + fname); - initRate(getDataPath("EMPairProduction/rate_" + fname + ".txt")); - initCumulativeRate(getDataPath("EMPairProduction/cdf_" + fname + ".txt")); + + if (!this->photonField->hasPositionDependence()){ + this->interactionRates = new InteractionRatesHomogeneous( + getDataPath("EMPairProduction/rate_" + fname + ".txt"), + getDataPath("EMPairProduction/cdf_" + fname + ".txt") + ); + + } else { + this->interactionRates = new InteractionRatesPositionDependent( + getDataPath("EMPairProduction/"+fname+"/Rate/"), + getDataPath("EMPairProduction/"+fname+"/CumulativeRate/"), + this->surface + ); + } } void EMPairProduction::setHaveElectrons(bool haveElectrons) { @@ -39,153 +67,119 @@ void EMPairProduction::setThinning(double thinning) { this->thinning = thinning; } -void EMPairProduction::initRate(std::string filename) { - std::ifstream infile(filename.c_str()); - - if (!infile.good()) - throw std::runtime_error("EMPairProduction: could not open file " + filename); - - // clear previously loaded interaction rates - tabEnergy.clear(); - tabRate.clear(); - - while (infile.good()) { - if (infile.peek() != '#') { - double a, b; - infile >> a >> b; - if (infile) { - tabEnergy.push_back(pow(10, a) * eV); - tabRate.push_back(b / Mpc); - } - } - infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n'); - } - infile.close(); +void EMPairProduction::setSurface(ref_ptr surface) { + this->surface = surface; } -void EMPairProduction::initCumulativeRate(std::string filename) { - std::ifstream infile(filename.c_str()); +ref_ptr EMPairProduction::getSurface() const { + return this->surface; +} - if (!infile.good()) - throw std::runtime_error("EMPairProduction: could not open file " + filename); +void EMPairProduction::setInteractionRates(ref_ptr intRates) { + this->interactionRates = intRates; +} - // clear previously loaded tables - tabE.clear(); - tabs.clear(); - tabCDF.clear(); - - // skip header - while (infile.peek() == '#') - infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n'); +ref_ptr EMPairProduction::getInteractionRates() const { + return this->interactionRates; +} - // read s values in first line - double a; - infile >> a; // skip first value - while (infile.good() and (infile.peek() != '\n')) { - infile >> a; - tabs.push_back(pow(10, a) * eV * eV); - } +void EMPairProduction::initRate(std::string path) { + this->interactionRates->initRate(path); +} - // read all following lines: E, cdf values - while (infile.good()) { - infile >> a; - if (!infile) - break; // end of file - tabE.push_back(pow(10, a) * eV); - std::vector cdf; - for (int i = 0; i < tabs.size(); i++) { - infile >> a; - cdf.push_back(a / Mpc); - } - tabCDF.push_back(cdf); - } - infile.close(); +void EMPairProduction::initCumulativeRate(std::string path) { + this->interactionRates->initCumulativeRate(path); } // Hold an data array to interpolate the energy distribution on class PPSecondariesEnergyDistribution { - private: - std::vector tab_s; - std::vector< std::vector > data; - size_t N; - - public: - // differential cross section for pair production for x = Epositron/Egamma, compare Lee 96 arXiv:9604098 - double dSigmadE_PPx(double x, double beta) { - double A = (x / (1. - x) + (1. - x) / x ); - double B = (1. / x + 1. / (1. - x) ); - double y = (1 - beta * beta); - return A + y * B - y * y / 4 * B * B; - } - - PPSecondariesEnergyDistribution() { - N = 1000; - size_t Ns = 1000; - double s_min = 4 * mec2 * mec2; - double s_max = 1e23 * eV * eV; - double dls = log(s_max / s_min) / Ns; - data = std::vector< std::vector >(Ns, std::vector(N)); - tab_s = std::vector(Ns + 1); - - for (size_t i = 0; i < Ns + 1; ++i) - tab_s[i] = s_min * exp(i*dls); // tabulate s bin borders - - for (size_t i = 0; i < Ns; i++) { - double s = s_min * exp(i*dls + 0.5*dls); - double beta = sqrt(1 - s_min/s); - double x0 = (1 - beta) / 2; - double dx = log((1 + beta) / (1 - beta)) / N; - - // cumulative midpoint integration - std::vector data_i(1000); - data_i[0] = dSigmadE_PPx(x0, beta) * expm1(dx); - for (size_t j = 1; j < N; j++) { - double x = x0 * exp(j*dx + 0.5*dx); - double binWidth = exp((j+1)*dx)-exp(j*dx); - data_i[j] = dSigmadE_PPx(x, beta) * binWidth + data_i[j-1]; - } - data[i] = data_i; - } - } - - // sample positron energy from cdf(E, s_kin) - double sample(double E0, double s) { - // get distribution for given s - size_t idx = std::lower_bound(tab_s.begin(), tab_s.end(), s) - tab_s.begin(); - if (idx > data.size()) - return NAN; - - std::vector s0 = data[idx]; - - // draw random bin - Random &random = Random::instance(); - size_t j = random.randBin(s0) + 1; - - double s_min = 4. * mec2 * mec2; - double beta = sqrtl(1. - s_min / s); - double x0 = (1. - beta) / 2.; +private: + std::vector tab_s; + std::vector< std::vector > data; + size_t N; + +public: + // differential cross section for pair production for x = Epositron/Egamma, compare Lee 96 arXiv:9604098 + double dSigmadE_PPx(double x, double beta) { + double A = (x / (1. - x) + (1. - x) / x ); + double B = (1. / x + 1. / (1. - x) ); + double y = (1 - beta * beta); + return A + y * B - y * y / 4 * B * B; + } + + PPSecondariesEnergyDistribution() { + N = 1000; + size_t Ns = 1000; + double s_min = 4 * mec2 * mec2; + double s_max = 1e23 * eV * eV; + double dls = log(s_max / s_min) / Ns; + data = std::vector< std::vector >(Ns, std::vector(N)); + tab_s = std::vector(Ns + 1); + + for (size_t i = 0; i < Ns + 1; ++i) + tab_s[i] = s_min * exp(i*dls); // tabulate s bin borders + + for (size_t i = 0; i < Ns; i++) { + double s = s_min * exp(i*dls + 0.5*dls); + double beta = sqrt(1 - s_min/s); + double x0 = (1 - beta) / 2; double dx = log((1 + beta) / (1 - beta)) / N; - double binWidth = x0 * (exp(j*dx) - exp((j-1)*dx)); - if (random.rand() < 0.5) - return E0 * (x0 * exp((j-1) * dx) + binWidth); - else - return E0 * (1 - (x0 * exp((j-1) * dx) + binWidth)); + + // cumulative midpoint integration + std::vector data_i(1000); + data_i[0] = dSigmadE_PPx(x0, beta) * expm1(dx); + for (size_t j = 1; j < N; j++) { + double x = x0 * exp(j*dx + 0.5*dx); + double binWidth = exp((j+1)*dx)-exp(j*dx); + data_i[j] = dSigmadE_PPx(x, beta) * binWidth + data_i[j-1]; + } + data[i] = data_i; } + } + + // sample positron energy from cdf(E, s_kin) + double sample(double E0, double s) { + // get distribution for given s + size_t idx = std::lower_bound(tab_s.begin(), tab_s.end(), s) - tab_s.begin(); + if (idx > data.size()) + return NAN; + + std::vector s0 = data[idx]; + + // draw random bin + Random &random = Random::instance(); + size_t j = random.randBin(s0) + 1; + + double s_min = 4. * mec2 * mec2; + double beta = sqrtl(1. - s_min / s); + double x0 = (1. - beta) / 2.; + double dx = log((1 + beta) / (1 - beta)) / N; + double binWidth = x0 * (exp(j*dx) - exp((j-1)*dx)); + if (random.rand() < 0.5) + return E0 * (x0 * exp((j-1) * dx) + binWidth); + else + return E0 * (1 - (x0 * exp((j-1) * dx) + binWidth)); + } + }; void EMPairProduction::performInteraction(Candidate *candidate) const { - // photon is lost after interacting - candidate->setActive(false); - - // check if secondary electron pair needs to be produced - if (not haveElectrons) - return; - // scale particle energy instead of background photon energy double z = candidate->getRedshift(); double E = candidate->current.getEnergy() * (1 + z); - + Vector3d position = candidate->current.getPosition(); + + // check if secondary electron pair needs to be produced + if (not haveElectrons) + return; + + std::vector tabE; + std::vector tabs; + std::vector> tabCDF; + + this->interactionRates->loadPerformInteractionTabs(position, tabE, tabs, tabCDF); + // check if in tabulated energy range if (E < tabE.front() or (E > tabE.back())) { KISS_LOG_WARNING @@ -207,13 +201,13 @@ void EMPairProduction::performInteraction(Candidate *candidate) const { double lo = std::max(4 * mec2 * mec2, tabs[j-1]); // first s-tabulation point below min(s_kin) = (2 me c^2)^2; ensure physical value double hi = tabs[j]; double s = lo + random.rand() * (hi - lo); - + // sample electron / positron energy static PPSecondariesEnergyDistribution interpolation; double Ee = interpolation.sample(E, s); double Ep = E - Ee; double f = Ep / E; - + // for some backgrounds Ee=nan due to precision limitations. if (not std::isfinite(Ee) || not std::isfinite(Ep)) { KISS_LOG_WARNING @@ -230,44 +224,49 @@ void EMPairProduction::performInteraction(Candidate *candidate) const { double w = 1. / pow(f, thinning); candidate->addSecondary(11, Ep / (1 + z), pos, w, interactionTag); } - if (random.rand() < pow(1 - f, thinning)){ + if (random.rand() < pow(1 - f, thinning)) { double w = 1. / pow(1 - f, thinning); - candidate->addSecondary(-11, Ee / (1 + z), pos, w, interactionTag); + candidate->addSecondary(-11, Ee / (1 + z), pos, w, interactionTag); } + + // cosmic ray photon is lost after interacting + candidate->setActive(false); + } void EMPairProduction::process(Candidate *candidate) const { + // check if photon if (candidate->current.getId() != 22) return; - + // scale particle energy instead of background photon energy double z = candidate->getRedshift(); double E = candidate->current.getEnergy() * (1 + z); - - // check if in tabulated energy range - if ((E < tabEnergy.front()) or (E > tabEnergy.back())) + Vector3d position = candidate->current.getPosition(); + + double rate = this->interactionRates->getProcessRate(E, position); + + if (rate < 0) return; - - // interaction rate - double rate = interpolate(E, tabEnergy, tabRate); + rate *= pow_integer<2>(1 + z) * photonField->getRedshiftScaling(z); - - // run this loop at least once to limit the step size + + // run this loop at least once to limit the step size double step = candidate->getCurrentStep(); Random &random = Random::instance(); do { double randDistance = -log(random.rand()) / rate; // check for interaction; if it doesn't ocurr, limit next step - if (step < randDistance) { + if (step < randDistance) { candidate->limitNextStep(limit / rate); } else { performInteraction(candidate); return; } - step -= randDistance; + step -= randDistance; } while (step > 0.); - + } void EMPairProduction::setInteractionTag(std::string tag) { diff --git a/src/module/EMTripletPairProduction.cpp b/src/module/EMTripletPairProduction.cpp index f66764374..bcaaafa05 100644 --- a/src/module/EMTripletPairProduction.cpp +++ b/src/module/EMTripletPairProduction.cpp @@ -1,28 +1,51 @@ #include "crpropa/module/EMTripletPairProduction.h" #include "crpropa/Units.h" #include "crpropa/Random.h" +#include "crpropa/Common.h" #include +#include #include #include +#include + +#if defined(__APPLE__) && defined(_LIBCPP_VERSION) + namespace fs = std::__fs::filesystem; +#else + namespace fs = std::filesystem; +#endif namespace crpropa { static const double mec2 = mass_electron * c_squared; -EMTripletPairProduction::EMTripletPairProduction(ref_ptr photonField, bool haveElectrons, double thinning, double limit) { +EMTripletPairProduction::EMTripletPairProduction(ref_ptr photonField, bool haveElectrons, double thinning, double limit, ref_ptr surface) { + + setSurface(surface); setPhotonField(photonField); setHaveElectrons(haveElectrons); setLimit(limit); setThinning(thinning); + } void EMTripletPairProduction::setPhotonField(ref_ptr photonField) { this->photonField = photonField; std::string fname = photonField->getFieldName(); setDescription("EMTripletPairProduction: " + fname); - initRate(getDataPath("EMTripletPairProduction/rate_" + fname + ".txt")); - initCumulativeRate(getDataPath("EMTripletPairProduction/cdf_" + fname + ".txt")); + if (!this->photonField->hasPositionDependence()){ + this->interactionRates = new InteractionRatesHomogeneous( + getDataPath("EMTripletPairProduction/rate_" + fname + ".txt"), + getDataPath("EMTripletPairProduction/cdf_" + fname + ".txt") + ); + + } else { + this->interactionRates = new InteractionRatesPositionDependent( + getDataPath("EMTripletPairProduction/"+fname+"/Rate/"), + getDataPath("EMTripletPairProduction/"+fname+"/CumulativeRate/"), + this->surface + ); + } } void EMTripletPairProduction::setHaveElectrons(bool haveElectrons) { @@ -37,96 +60,60 @@ void EMTripletPairProduction::setThinning(double thinning) { this->thinning = thinning; } -void EMTripletPairProduction::initRate(std::string filename) { - std::ifstream infile(filename.c_str()); - - if (!infile.good()) - throw std::runtime_error("EMTripletPairProduction: could not open file " + filename); - - // clear previously loaded interaction rates - tabEnergy.clear(); - tabRate.clear(); - - while (infile.good()) { - if (infile.peek() != '#') { - double a, b; - infile >> a >> b; - if (infile) { - tabEnergy.push_back(pow(10, a) * eV); - tabRate.push_back(b / Mpc); - } - } - infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n'); - } - infile.close(); +void EMTripletPairProduction::setSurface(ref_ptr surface) { + this->surface = surface; } -void EMTripletPairProduction::initCumulativeRate(std::string filename) { - std::ifstream infile(filename.c_str()); +ref_ptr EMTripletPairProduction::getSurface() const { + return this->surface; +} - if (!infile.good()) - throw std::runtime_error( - "EMTripletPairProduction: could not open file " + filename); +void EMTripletPairProduction::setInteractionRates(ref_ptr intRates) { + this->interactionRates = intRates; +} - // clear previously loaded tables - tabE.clear(); - tabs.clear(); - tabCDF.clear(); - - // skip header - while (infile.peek() == '#') - infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n'); +ref_ptr EMTripletPairProduction::getInteractionRates() const { + return this->interactionRates; +} - // read s values in first line - double a; - infile >> a; // skip first value - while (infile.good() and (infile.peek() != '\n')) { - infile >> a; - tabs.push_back(pow(10, a) * eV * eV); - } +void EMTripletPairProduction::initRate(std::string path) { + this->interactionRates->initRate(path); +} - // read all following lines: E, cdf values - while (infile.good()) { - infile >> a; - if (!infile) - break; // end of file - tabE.push_back(pow(10, a) * eV); - std::vector cdf; - for (int i = 0; i < tabs.size(); i++) { - infile >> a; - cdf.push_back(a / Mpc); - } - tabCDF.push_back(cdf); - } - infile.close(); +void EMTripletPairProduction::initCumulativeRate(std::string path) { + this->interactionRates->initCumulativeRate(path); } void EMTripletPairProduction::performInteraction(Candidate *candidate) const { - int id = candidate->current.getId(); - if (abs(id) != 11) - return; - + // scale the particle energy instead of background photons double z = candidate->getRedshift(); double E = candidate->current.getEnergy() * (1 + z); - + Vector3d position = candidate->current.getPosition(); + + std::vector tabE; + std::vector tabs; + std::vector> tabCDF; + + this->interactionRates->loadPerformInteractionTabs(position, tabE, tabs, tabCDF); + if (E < tabE.front() or E > tabE.back()) return; - + // sample the value of eps Random &random = Random::instance(); size_t i = closestIndex(E, tabE); size_t j = random.randBin(tabCDF[i]); double s_kin = pow(10, log10(tabs[j]) + (random.rand() - 0.5) * 0.1); double eps = s_kin / 4. / E; // random background photon energy - + // Use approximation from A. Mastichiadis et al., Astroph. Journ. 300:178-189 (1986), eq. 30. // This approx is valid only for alpha >=100 where alpha = p0*eps*costheta - E0*eps // For our purposes, me << E0 --> p0~E0 --> alpha = E0*eps*(costheta - 1) >= 100 double Epp = 5.7e-1 * pow(eps / mec2, -0.56) * pow(E / mec2, 0.44) * mec2; - + double f = Epp / E; - + if (haveElectrons) { Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(), candidate->current.getPosition()); if (random.rand() < pow(1 - f, thinning)) { @@ -138,9 +125,11 @@ void EMTripletPairProduction::performInteraction(Candidate *candidate) const { candidate->addSecondary(-11, Epp / (1 + z), pos, w, interactionTag); } } + // Update the primary particle energy. // This is done after adding the secondaries to correctly set the secondaries parent candidate->current.setEnergy((E - 2 * Epp) / (1. + z)); + } void EMTripletPairProduction::process(Candidate *candidate) const { @@ -148,32 +137,34 @@ void EMTripletPairProduction::process(Candidate *candidate) const { int id = candidate->current.getId(); if (abs(id) != 11) return; - + // scale the particle energy instead of background photons double z = candidate->getRedshift(); double E = (1 + z) * candidate->current.getEnergy(); - - // check if in tabulated energy range - if ((E < tabEnergy.front()) or (E > tabEnergy.back())) - return; - - // cosmological scaling of interaction distance (comoving) + Vector3d position = candidate->current.getPosition(); + double scaling = pow_integer<2>(1 + z) * photonField->getRedshiftScaling(z); - double rate = scaling * interpolate(E, tabEnergy, tabRate); - + double rate = this->interactionRates->getProcessRate(E, position); + + if (rate < 0) + return; + + rate *= scaling; + // run this loop at least once to limit the step size double step = candidate->getCurrentStep(); Random &random = Random::instance(); do { double randDistance = -log(random.rand()) / rate; // check for interaction; if it doesn't occur, limit next step - if (step < randDistance) { + if (step < randDistance) { candidate->limitNextStep(limit / rate); return; } performInteraction(candidate); - step -= randDistance; + step -= randDistance; } while (step > 0.); + } void EMTripletPairProduction::setInteractionTag(std::string tag) { diff --git a/test/testInteraction.cpp b/test/testInteraction.cpp index 86758bdba..83864ba15 100644 --- a/test/testInteraction.cpp +++ b/test/testInteraction.cpp @@ -17,6 +17,14 @@ #include +#if defined(__APPLE__) && defined(_LIBCPP_VERSION) + #include + namespace fs = std::__fs::filesystem; +#else + #include + namespace fs = std::filesystem; +#endif + namespace crpropa { // ElectronPairProduction ----------------------------------------------------- @@ -40,7 +48,7 @@ TEST(ElectronPairProduction, allBackgrounds) { epp.setPhotonField(irb); irb = new IRB_Stecker16_lower(); epp.setPhotonField(irb); - irb = new IRB_Finke22(); + irb = new IRB_Finke22(); epp.setPhotonField(irb); } @@ -562,7 +570,7 @@ TEST(PhotoPionProduction, allBackgrounds) { ppp.setPhotonField(irb); irb = new IRB_Stecker16_lower(); ppp.setPhotonField(irb); - irb = new IRB_Finke22(); + irb = new IRB_Finke22(); ppp.setPhotonField(irb); ref_ptr urb = new URB_Protheroe96(); ppp.setPhotonField(urb); @@ -710,7 +718,16 @@ TEST(EMPairProduction, allBackgrounds) { em.setPhotonField(ebl); ref_ptr urb = new URB_Protheroe96(); em.setPhotonField(urb); - ebl = new IRB_Stecker05(); + + fs::path dir = getDataPath("") + "Scaling/ISRF_Freudenreich98"; + if (!fs::exists(dir)) { + std::cout << "Photon background tables not available" << std::endl; + } else { + ref_ptr isrf = new ISRF_Freudenreich98(nullptr); + em.setPhotonField(isrf); + } + + ebl = new IRB_Stecker05(); em.setPhotonField(ebl); ebl = new IRB_Franceschini08(); em.setPhotonField(ebl); @@ -783,8 +800,8 @@ TEST(EMPairProduction, secondaries) { Etot += s.current.getEnergy(); } - // test energy conservation - EXPECT_DOUBLE_EQ(Ep, Etot); + // test energy conservation + EXPECT_NEAR(Ep, Etot, 1E-9); } } } @@ -815,6 +832,15 @@ TEST(EMDoublePairProduction, allBackgrounds) { em.setPhotonField(ebl); ref_ptr urb = new URB_Protheroe96(); em.setPhotonField(urb); + + fs::path dir = getDataPath("") + "Scaling/ISRF_Freudenreich98"; + if (!fs::exists(dir)) { + std::cout << "Photon background tables not available" << std::endl; + } else { + ref_ptr isrf = new ISRF_Freudenreich98(nullptr); + em.setPhotonField(isrf); + } + ebl = new IRB_Stecker05(); em.setPhotonField(ebl); ebl = new IRB_Franceschini08(); @@ -921,6 +947,15 @@ TEST(EMTripletPairProduction, allBackgrounds) { em.setPhotonField(ebl); ref_ptr urb = new URB_Protheroe96(); em.setPhotonField(urb); + + fs::path dir = getDataPath("") + "Scaling/ISRF_Freudenreich98"; + if (!fs::exists(dir)) { + std::cout << "Photon background tables not available" << std::endl; + } else { + ref_ptr isrf = new ISRF_Freudenreich98(nullptr); + em.setPhotonField(isrf); + } + ebl = new IRB_Stecker05(); em.setPhotonField(ebl); ebl = new IRB_Franceschini08(); @@ -1028,6 +1063,15 @@ TEST(EMInverseComptonScattering, allBackgrounds) { em.setPhotonField(ebl); ref_ptr urb = new URB_Protheroe96(); em.setPhotonField(urb); + + fs::path dir = getDataPath("") + "Scaling/ISRF_Freudenreich98"; + if (!fs::exists(dir)) { + std::cout << "Photon background tables not available" << std::endl; + } else { + ref_ptr isrf = new ISRF_Freudenreich98(nullptr); + em.setPhotonField(isrf); + } + ebl = new IRB_Stecker05(); em.setPhotonField(ebl); ebl = new IRB_Franceschini08(); diff --git a/test/testPythonExtension.py b/test/testPythonExtension.py index 21c2f8d55..913721a92 100644 --- a/test/testPythonExtension.py +++ b/test/testPythonExtension.py @@ -155,8 +155,9 @@ def __init__(self, density): self.density = density self.fieldName = 'testCustomPhotonField' self.isRedshiftDependent = True + self.isPositionDependent = False - def getPhotonDensity(self, energy, z): + def getPhotonDensity(self, energy, z, pos): return self.density def getFieldName(self): @@ -165,15 +166,18 @@ def getFieldName(self): def hasRedshiftDependence(self): return self.isRedshiftDependent + def hasPositionDependence(self): + return self.isPositionDependent + photonDensity = 10 photonField = CustomPhotonField(photonDensity) energy = 10*crp.GeV z = 0 self.assertEqual(photonDensity, photonField.getPhotonDensity(energy, z)) self.assertEqual('testCustomPhotonField', photonField.getFieldName()) - self.assertEqual(True, photonField.hasRedshiftDependence()) - - + self.assertEqual(True, photonField.hasRedshiftDependence(), photonField.hasPositionDependence()) + self.assertEqual(True, photonField.hasPositionDependence(), photonField.hasPositionDependence()) + def testCustomPhotonField(self): class CustomPhotonField(crp.PhotonField): def __init__(self, val): @@ -183,7 +187,7 @@ def __init__(self, val): def getFieldName(self): return 'CMB' - def getPhotonDensity(self, ePhoton, z): + def getPhotonDensity(self, ePhoton, z, pos): return self.val constDensity = 1