Fixup for sph

For details see https://clang.llvm.org/extra/clang-tidy/checks/modernize-replace-random-shuffle.html

CMakeLists.txt

@@ -42,6 +42,17 @@ ENDIF(BUILD_PYTHON)
 
 MESSAGE(STATUS "Using the target ${CosmoTool_local} to build python module")
 
+find_library(ZLIB_LIBRARY z)
+find_library(DL_LIBRARY dl)
+find_library(RT_LIBRARY rt)
+if (RT_LIBRARY)
+  SET(RT_DEP rt)
+  message(STATUS "RT found, linking against it")
+else()
+  SET(RT_DEP)
+endif()
+find_library(MATH_LIBRARY m)
+
 include(${CMAKE_SOURCE_DIR}/external/external_build.cmake)
 
 IF(YORICK_SUPPORT)
@@ -55,9 +66,6 @@ IF(YORICK_SUPPORT)
 ENDIF(YORICK_SUPPORT)
 
 find_program(CYTHON NAMES cython3 cython)
-find_library(ZLIB_LIBRARY z)
-find_library(DL_LIBRARY dl)
-find_library(MATH_LIBRARY m)
 
 set(NETCDF_FIND_REQUIRED ${YORICK_SUPPORT})
 set(GSL_FIND_REQUIRED TRUE)
@@ -71,7 +79,7 @@ SET(CPACK_PACKAGE_VENDOR "Guilhem Lavaux")
 SET(CPACK_RESOURCE_FILE_LICENSE "${CMAKE_CURRENT_SOURCE_DIR}/LICENCE_CeCILL_V2")
 SET(CPACK_PACKAGE_VERSION_MAJOR "1")
 SET(CPACK_PACKAGE_VERSION_MINOR "3")
-SET(CPACK_PACKAGE_VERSION_PATCH "1${EXTRA_VERSION}")
+SET(CPACK_PACKAGE_VERSION_PATCH "4${EXTRA_VERSION}")
 SET(CPACK_PACKAGE_INSTALL_DIRECTORY "CosmoToolbox-${CPACK_PACKAGE_VERSION_MAJOR}.${CPACK_PACKAGE_VERSION_MINOR}")
 SET(CPACK_STRIP_FILES "lib/libCosmoTool.so")
 SET(CPACK_SOURCE_IGNORE_FILES 

MANIFEST.in

@@ -1,4 +1,5 @@
 include  .gitignore
+include  pyproject.toml
 include  CMakeLists.txt
 include  FindNumPy.cmake
 include  FindPyLibs.cmake
@@ -16,14 +17,17 @@ include  doc/source/cpplibrary.rst
 include  doc/source/index.rst
 include  doc/source/intro.rst
 include  doc/source/pythonmodule.rst
+include  external/config.guess
+include  external/config.sub
 include  external/external_build.cmake
-include  external/libsharp-6077806.tar.gz
+include  external/libsharp-8d51946.tar.gz
 include  external/omptl-20120422.tar.bz2
 include  external/patch-omptl
 include  python/CMakeLists.txt
 include  python/_cosmo_bispectrum.cpp
 include  python/_cosmo_cic.pyx
 include  python/_cosmo_power.pyx
+include  python/_cosmomath.pyx
 include  python/_cosmotool.pyx
 include  python/_fast_interp.pyx
 include  python/_project.pyx
@@ -93,6 +97,7 @@ include  sample/testkd3.cpp
 include  setup.py
 include  src/CMakeLists.txt
 include  src/algo.hpp
+include  src/numpy_adaptors.hpp
 include  src/bqueue.hpp
 include  src/bqueue.tcc
 include  src/bsp_simple.hpp

builder/build-wheels.sh

@@ -8,12 +8,15 @@ export CC CXX
 
 # Install a system package required by our library
 #yum install -y atlas-devel
-yum install -y cmake3 gsl-devel zlib-devel fftw3-devel
+yum install -y cmake3 gsl-devel zlib-devel fftw3-devel libffi-devel hdf5 hdf5-devel
 
 ln -fs /usr/bin/cmake3 /usr/bin/cmake
 
 
-ALL_PYTHON="cp36-cp36m cp37-cp37m cp38-cp38 cp39-cp39" 
+test -d /io/wheelhouse || mkdir /io/wheelhouse
+test -d /io/wheelhouse/fix || mkdir /io/wheelhouse/fix
+
+ALL_PYTHON="cp39-cp39 cp310-cp310"
 
 # Compile wheels
 for pkg in $ALL_PYTHON; do
@@ -21,12 +24,12 @@ for pkg in $ALL_PYTHON; do
 #    "${PYBIN}/pip" install -r /io/dev-requirements.txt
     "${PYBIN}/pip" install setuptools wheel Cython 
     "${PYBIN}/pip" install -r /io/requirements.txt
-    "${PYBIN}/pip" wheel -vvv /io/ -w wheelhouse/
+    "${PYBIN}/pip" wheel -vvv /io/ -w /io/wheelhouse/
 done
 
 # Bundle external shared libraries into the wheels
-for whl in wheelhouse/cosmotool*.whl; do
-    auditwheel repair "$whl" --plat $PLAT -w /io/wheelhouse/
+for whl in /io/wheelhouse/cosmotool*linux*.whl; do
+    auditwheel repair "$whl" --plat $PLAT -w /io/wheelhouse/fix
 done
 
 # Install packages and test

builder/start.sh

@@ -9,4 +9,4 @@ if ! [ -e ${d}/setup.py ] ; then
   exit 1
 fi
 
-podman run -ti  --rm -e PLAT=manylinux2010_x86_64 -v ${d}:/io:Z quay.io/pypa/manylinux2010_x86_64 /io/builder/build-wheels.sh
+podman run -ti  --rm -e PLAT=manylinux2014_x86_64 -v ${d}:/io:Z quay.io/pypa/manylinux2014_x86_64 /io/builder/build-wheels.sh

external/external_build.cmake

@@ -2,16 +2,19 @@ include(FindOpenMP)
 
 OPTION(ENABLE_OPENMP "Set to Yes if Healpix and/or you need openMP" OFF)
 
+SET(SOURCE_PREFIX ${CMAKE_SOURCE_DIR})
 SET(FFTW_URL "http://www.fftw.org/fftw-3.3.3.tar.gz" CACHE STRING "URL to download FFTW from")
 SET(EIGEN_URL "https://gitlab.com/libeigen/eigen/-/archive/3.3.7/eigen-3.3.7.tar.bz2" CACHE STRING "URL to download Eigen from")
 SET(GENGETOPT_URL "ftp://ftp.gnu.org/gnu/gengetopt/gengetopt-2.22.5.tar.gz" CACHE STRING "URL to download gengetopt from")
-SET(HDF5_URL "https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.8/hdf5-1.8.18/src/hdf5-1.8.18.tar.bz2" CACHE STRING "URL to download HDF5 from")
-SET(NETCDF_URL "ftp://ftp.unidata.ucar.edu/pub/netcdf/netcdf-4.5.0.tar.gz" CACHE STRING "URL to download NetCDF from")
-SET(NETCDFCXX_URL "https://github.com/Unidata/netcdf-cxx4/archive/v4.3.0.tar.gz" CACHE STRING "URL to download NetCDF-C++ from")
-SET(BOOST_URL "https://boostorg.jfrog.io/artifactory/main/release/1.74.0/source/boost_1_74_0.tar.bz2" CACHE STRING "URL to download Boost from")
-SET(GSL_URL "ftp://ftp.gnu.org/gnu/gsl/gsl-1.15.tar.gz" CACHE STRING "URL to download GSL from ")
+SET(HDF5_URL "https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.12/hdf5-1.12.2/src/hdf5-1.12.2.tar.gz" CACHE STRING "URL to download HDF5 from")
+SET(NETCDF_URL "https://downloads.unidata.ucar.edu/netcdf-c/4.9.2/netcdf-c-4.9.2.tar.gz" CACHE STRING "URL to download NetCDF from")
+SET(NETCDFCXX_URL "https://github.com/Unidata/netcdf-cxx4/archive/v4.3.1.tar.gz" CACHE STRING "URL to download NetCDF-C++ from")
+SET(BOOST_URL "https://boostorg.jfrog.io/artifactory/main/release/1.82.0/source/boost_1_82_0.tar.gz" CACHE STRING "URL to download Boost from")
+SET(GSL_URL "https://ftpmirror.gnu.org/gsl/gsl-2.7.tar.gz" CACHE STRING "URL to download GSL from ")
 mark_as_advanced(FFTW_URL EIGEN_URL HDF5_URL NETCDF_URL BOOST_URL GSL_URL)
 
+file(MAKE_DIRECTORY ${SOURCE_PREFIX}/downloads)
+
 SET(all_deps)
 
 MACRO(CHECK_CHANGE_STATE VAR) 
@@ -36,12 +39,13 @@ CHECK_CHANGE_STATE(INTERNAL_DLIB DLIB_INCLUDE_DIR DLIB_LIBRARIES)
 
 IF(ENABLE_OPENMP)
   IF (NOT OPENMP_FOUND)
-    MESSAGE(ERROR "No known compiler option for enabling OpenMP")
-  ENDIF(NOT OPENMP_FOUND)
-
+    MESSAGE(NOTICE "No known compiler option for enabling OpenMP")
+  ELSE()
     SET(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
     SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
     SET(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_C_FLAGS}")
+  ENDIF(NOT OPENMP_FOUND)
+
 ENDIF(ENABLE_OPENMP)
 
 
@@ -55,17 +59,25 @@ if (ENABLE_SHARP)
   SET(DEP_BUILD ${BUILD_PREFIX}/sharp-prefix/src/sharp/auto)
   IF(NOT ENABLE_OPENMP)
     SET(SHARP_OPENMP --disable-openmp)
+  ELSE()
+    SET(SHARP_OPENMP)
   ENDIF()
   SET(CUTILS_LIBRARY ${DEP_BUILD}/lib/libc_utils.a)
   SET(FFTPACK_LIBRARY ${DEP_BUILD}/lib/libfftpack.a)
   SET(SHARP_LIBRARY ${DEP_BUILD}/lib/libsharp.a)
   SET(SHARP_LIBRARIES ${SHARP_LIBRARY} ${FFTPACK_LIBRARY} ${CUTILS_LIBRARY})
   SET(SHARP_INCLUDE_PATH ${DEP_BUILD}/include)
+  message(STATUS "Flags: ${CMAKE_C_FLAGS}")
   ExternalProject_Add(sharp
-    URL ${CMAKE_SOURCE_DIR}/external/libsharp-6077806.tar.gz
+    URL ${CMAKE_SOURCE_DIR}/external/libsharp-8d51946.tar.gz
     PREFIX ${BUILD_PREFIX}/sharp-prefix
+    DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
     BUILD_IN_SOURCE 1 
-    CONFIGURE_COMMAND autoconf && ./configure "CC=${CMAKE_C_COMPILER}" "CXX=${CMAKE_CXX_COMPILER}" --prefix=${DEP_BUILD} ${SHARP_OPENMP}
+    CONFIGURE_COMMAND
+          cp -f ${CMAKE_SOURCE_DIR}/external/config.guess . &&
+          cp -f ${CMAKE_SOURCE_DIR}/external/config.sub . &&
+          autoconf &&
+          ./configure "CC=${CMAKE_C_COMPILER}" "CXX=${CMAKE_CXX_COMPILER}" "CFLAGS=${CMAKE_C_FLAGS}" "LDFLAGS=${CMAKE_EXE_LINKER_FLAGS}" --prefix=${DEP_BUILD} ${SHARP_OPENMP}
     INSTALL_COMMAND echo "No install"
     BUILD_BYPRODUCTS ${SHARP_LIBRARIES}
   )
@@ -82,27 +94,23 @@ if (INTERNAL_HDF5)
   ExternalProject_Add(hdf5
     PREFIX ${BUILD_PREFIX}/hdf5-prefix
     URL ${HDF5_URL}
-    URL_HASH MD5=29117bf488887f89888f9304c8ebea0b
-    CMAKE_ARGS
-       -DCMAKE_INSTALL_PREFIX=${EXT_INSTALL}
-       -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER}
-       -DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER}
-       -DHDF5_BUILD_CPP_LIB=ON
-       -DHDF5_BUILD_TOOLS=ON
-       -DHDF5_BUILD_HL_LIB=ON
-       -DBUILD_SHARED_LIBS=OFF
+    URL_HASH MD5=30172c75e436d7f2180e274071a4ca97
+    DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
+    CONFIGURE_COMMAND
+        ${HDF5_SOURCE_DIR}/configure "CC=${CMAKE_C_COMPILER}" "CXX=${CMAKE_CXX_COMPILER}" "CFLAGS=${CMAKE_C_FLAGS}" "LDFLAGS=${CMAKE_EXE_LINKER_FLAGS}" --prefix=${HDF5_BIN_DIR} --disable-shared  --enable-cxx --enable-hl --enable-tools --with-pic
+    INSTALL_COMMAND make install
   )
   SET(cosmotool_DEPS ${cosmotool_DEPS} hdf5)
   SET(hdf5_built hdf5)
   SET(ENV{HDF5_ROOT} ${HDF5_BIN_DIR})
   SET(HDF5_ROOTDIR ${HDF5_BIN_DIR})
   SET(CONFIGURE_LDFLAGS "${CONFIGURE_LDFLAGS} -L${HDF5_BIN_DIR}/lib")
-  SET(CONFIGURE_LIBS "${CONFIGURE_LIBS} -ldl")
-  set(HDF5_C_STATIC_LIBRARY ${HDF5_BIN_DIR}/lib/libhdf5-static.a)
-  set(HDF5_HL_STATIC_LIBRARY ${HDF5_BIN_DIR}/lib/libhdf5_hl-static.a)
-  set(HDF5_LIBRARIES ${HDF5_BIN_DIR}/lib/libhdf5-static.a CACHE STRING "HDF5 lib" FORCE)
-  set(HDF5_HL_LIBRARIES ${HDF5_BIN_DIR}/lib/libhdf5_hl-static.a CACHE STRING "HDF5 HL lib" FORCE)
-  set(HDF5_CXX_LIBRARIES ${HDF5_BIN_DIR}/lib/libhdf5_cpp-static.a CACHE STRING "HDF5 C++ lib" FORCE)
+  SET(CONFIGURE_LIBS "${CONFIGURE_LIBS} -ldl ${RT_LIBRARY}")
+  set(HDF5_C_STATIC_LIBRARY ${HDF5_BIN_DIR}/lib/libhdf5.a)
+  set(HDF5_HL_STATIC_LIBRARY ${HDF5_BIN_DIR}/lib/libhdf5_hl.a)
+  set(HDF5_LIBRARIES ${HDF5_BIN_DIR}/lib/libhdf5.a CACHE STRING "HDF5 lib" FORCE)
+  set(HDF5_HL_LIBRARIES ${HDF5_BIN_DIR}/lib/libhdf5_hl.a CACHE STRING "HDF5 HL lib" FORCE)
+  set(HDF5_CXX_LIBRARIES ${HDF5_BIN_DIR}/lib/libhdf5_cpp.a CACHE STRING "HDF5 C++ lib" FORCE)
   SET(HDF5_INCLUDE_DIRS ${HDF5_BIN_DIR}/include CACHE STRING "HDF5 include path" FORCE)
   mark_as_advanced(HDF5_LIBRARIES HDF5_CXX_LIBRARIES HDF5_INCLUDE_DIRS)
 
@@ -160,24 +168,28 @@ if (INTERNAL_NETCDF)
   SET(NETCDF_CONFIG_COMMAND ${NETCDF_SOURCE_DIR}/configure
          --prefix=${NETCDF_BIN_DIR} --libdir=${NETCDF_BIN_DIR}/lib
          --enable-netcdf-4  --with-pic --disable-shared --disable-dap 
-         --disable-cdmremote --disable-rpc --enable-cxx-4 
+         --disable-byterange --disable-cdmremote --disable-rpc --enable-cxx-4 
          --disable-examples ${EXTRA_NC_FLAGS} CC=${CMAKE_C_COMPILER}
          CXX=${CMAKE_CXX_COMPILER})
   list(INSERT CMAKE_PREFIX_PATH 0 ${EXT_INSTALL})
   string(REPLACE ";" "|" CMAKE_PREFIX_PATH_ALT_SEP "${CMAKE_PREFIX_PATH}")
   ExternalProject_Add(netcdf
     DEPENDS ${hdf5_built}
+    URL_HASH MD5=f48ee01534365006934f0c63d4055ea0
     PREFIX ${BUILD_PREFIX}/netcdf-prefix
+    DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
     URL ${NETCDF_URL}
     LIST_SEPARATOR |
     CMAKE_ARGS
 		-DCMAKE_PREFIX_PATH=${CMAKE_PREFIX_PATH_ALT_SEP}
 		-DCMAKE_C_COMPILER=${CMAKE_C_COMPILER}
 		-DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER}
+                -DNC_EXTRA_DEPS=${RT_DEP}
 		-DBUILD_SHARED_LIBS=OFF
 		-DBUILD_TESTING=OFF
 		-DCMAKE_BUILD_TYPE=Release
 		-DENABLE_NETCDF4=ON
+                -DENABLE_BYTERANGE=FALSE
 		-DCMAKE_POSITION_INDEPENDENT_CODE=ON
 		-DENABLE_DAP=OFF
 		-DCMAKE_INSTALL_PREFIX=${NETCDF_BIN_DIR}
@@ -191,12 +203,17 @@ if (INTERNAL_NETCDF)
   ExternalProject_Add(netcdf-c++
     DEPENDS ${hdf5_built} netcdf
     PREFIX ${BUILD_PREFIX}/netcdf-c++-prefix
+    DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
     URL ${NETCDFCXX_URL}
 	CMAKE_ARGS
+		-DCMAKE_PREFIX_PATH=${CMAKE_PREFIX_PATH_ALT_SEP}
 		-DCMAKE_C_COMPILER=${CMAKE_C_COMPILER}
 		-DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER}
 		-DBUILD_SHARED_LIBS=OFF
                 -DCMAKE_POSITION_INDEPENDENT_CODE=ON
+                -DHDF5_C_LIBRARY=${HDF5_C_STATIC_LIBRARY}
+                -DHDF5_HL_LIBRARY=${HDF5_HL_STATIC_LIBRARY}
+                -DHDF5_INCLUDE_DIR=${HDF5_INCLUDE_DIRS}
 		-DBUILD_TESTING=OFF
 		-DCMAKE_BUILD_TYPE=Release
 		-DCMAKE_INSTALL_PREFIX=${NETCDF_BIN_DIR}
@@ -234,7 +251,8 @@ if (INTERNAL_BOOST)
   ExternalProject_Add(boost
     URL ${BOOST_URL}
     PREFIX ${BUILD_PREFIX}/boost-prefix
-    URL_HASH MD5=da07ca30dd1c0d1fdedbd487efee01bd
+    DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
+    URL_HASH MD5=f7050f554a65f6a42ece221eaeec1660
     CONFIGURE_COMMAND 
            ${BOOST_SOURCE_DIR}/bootstrap.sh --prefix=${CMAKE_BINARY_DIR}/ext_build/boost
     BUILD_IN_SOURCE 1
@@ -277,6 +295,7 @@ IF(INTERNAL_GSL)
   ExternalProject_Add(gsl
     URL ${GSL_URL}
     PREFIX ${BUILD_PREFIX}/gsl-prefix
+    DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
     CONFIGURE_COMMAND ${GSL_SOURCE_DIR}/configure
            --prefix=${EXT_INSTALL} --disable-shared
            --with-pic
@@ -329,6 +348,7 @@ IF(INTERNAL_FFTW)
   ExternalProject_Add(fftw
      URL ${FFTW_URL}
      PREFIX ${BUILD_PREFIX}/fftw-prefix
+    DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
      CONFIGURE_COMMAND 
            ${FFTW_SOURCE}/configure
                  --prefix=${EXT_INSTALL}
@@ -358,6 +378,7 @@ IF (INTERNAL_EIGEN)
   ExternalProject_Add(eigen
      URL ${EIGEN_URL}
      URL_HASH MD5=b9e98a200d2455f06db9c661c5610496
+    DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
      PREFIX ${BUILD_PREFIX}/eigen-prefix
      CMAKE_ARGS -DCMAKE_INSTALL_PREFIX=${EXT_INSTALL}
        -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER}
@@ -390,6 +411,7 @@ SET(cosmotool_DEPS ${cosmotool_DEPS} omptl)
 SET(OMPTL_BUILD_DIR ${BUILD_PREFIX}/omptl-prefix/src/omptl)
 ExternalProject_Add(omptl
   PREFIX ${BUILD_PREFIX}/omptl-prefix
+  DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
   URL ${CMAKE_SOURCE_DIR}/external/omptl-20120422.tar.bz2
   CONFIGURE_COMMAND echo "No configure"
   BUILD_COMMAND echo "No build"

external/patch-omptl

@@ -1,6 +1,6 @@
 diff -ur omptl.old/omptl_algorithm omptl/omptl_algorithm
---- omptl.old/omptl_algorithm	2012-04-22 16:29:41.000000000 +0200
-+++ omptl/omptl_algorithm	2021-06-20 15:40:29.000000000 +0200
+--- omptl.old/omptl_algorithm	2022-06-19 08:21:39.815498672 +0200
++++ omptl/omptl_algorithm	2022-06-19 08:21:52.953544672 +0200
 @@ -20,7 +20,7 @@
  #define OMPTL_ALGORITHM 1
  
@@ -23,11 +23,13 @@ diff -ur omptl.old/omptl_algorithm omptl/omptl_algorithm
  
  #endif /* OMPTL_ALGORITHM */
 diff -ur omptl.old/omptl_algorithm_par.h omptl/omptl_algorithm_par.h
---- omptl.old/omptl_algorithm_par.h	2012-04-22 16:29:41.000000000 +0200
-+++ omptl/omptl_algorithm_par.h	2021-06-20 15:40:50.000000000 +0200
-@@ -21,8 +21,8 @@
+--- omptl.old/omptl_algorithm_par.h	2022-06-19 08:21:39.816498675 +0200
++++ omptl/omptl_algorithm_par.h	2022-06-19 08:23:50.705956964 +0200
+@@ -20,9 +20,10 @@
+ #include <utility>
  #include <cmath>
  #include <cstdlib>
++#include <random>
  
 -#include <omptl/omptl_tools.h>
 -#include <omptl/omptl_numeric>
@@ -36,7 +38,78 @@ diff -ur omptl.old/omptl_algorithm_par.h omptl/omptl_algorithm_par.h
  
  #include <iterator>
  
-@@ -1700,7 +1700,7 @@
+@@ -510,7 +511,7 @@
+ 
+ 		typename std::vector<Iterator>::iterator result =
+ 		    std::find_if(results.begin(),results.end(),
+-			std::bind2nd(std::not_equal_to<Iterator>(), last) );
++			std::bind(std::not_equal_to<Iterator>(), std::placeholders::_1, last) );
+ 
+ 		if ( result != results.end() )
+ 			return *result;
+@@ -569,7 +570,7 @@
+ 
+ 		const typename std::vector<Iterator>::iterator result
+ 			 = std::find_if(results.begin(), results.end(),
+-			std::bind2nd(std::not_equal_to<Iterator>(), last) );
++			std::bind(std::not_equal_to<Iterator>(), std::placeholders::_1, last) );
+ 
+ 		if ( result != results.end() )
+ 			return *result;
+@@ -654,7 +655,7 @@
+ 
+ 		const typename std::vector<Iterator>::iterator
+ 		result = std::find_if(results.begin(), results.end(),
+-					std::bind2nd(std::not_equal_to<Iterator>(), last1));
++					std::bind(std::not_equal_to<Iterator>(), std::placeholders::_1, last1));
+ 
+ 		if ( result != results.end() )
+ 			return *result;
+@@ -953,7 +954,7 @@
+ 		results[t] = std::lower_bound(partitions[t].first, partitions[t].second, value, comp);
+ 
+ 	const typename std::vector<ForwardIterator>::iterator result =
+-		std::find_if(results.begin(), results.end(), std::bind2nd(std::not_equal_to<ForwardIterator>(), last) );
++		std::find_if(results.begin(), results.end(), std::bind(std::not_equal_to<ForwardIterator>(), std::placeholders::_1, last) );
+ 
+ 	if (result != results.end())
+ 		return *result;
+@@ -1179,7 +1180,7 @@
+ 
+ namespace detail
+ {
+-	
++
+ template<typename Iterator, class StrictWeakOrdering>
+ Iterator _pivot_range(Iterator first, Iterator last,
+ 	const typename std::iterator_traits<Iterator>::value_type pivot,
+@@ -1309,14 +1310,14 @@
+ void random_shuffle(RandomAccessIterator first, RandomAccessIterator last,
+                     const unsigned P)
+ {
+-	std::random_shuffle(first, last);
++	std::shuffle(first, last, std::mt19937(std::random_device()()));
+ }
+ 
+ template <class RandomAccessIterator, class RandomNumberGenerator>
+ void random_shuffle(RandomAccessIterator first, RandomAccessIterator last,
+                     RandomNumberGenerator& rgen, const unsigned P)
+ {
+-	std::random_shuffle(first, last, rgen);
++	std::shuffle(first, last, std::mt19937(std::random_device()()));
+ }
+ 
+ // Not (yet) parallelized, not straightforward due to possible dependencies
+@@ -1472,7 +1473,7 @@
+                             const T& new_value, const unsigned P)
+ {
+ 	return ::omptl::replace_copy_if(first, last, result,
+-			std::bind2nd(std::equal_to<T>(), old_value), new_value, P);
++			std::bind(std::equal_to<T>(), std::placeholders::_1, old_value), new_value, P);
+ }
+ 
+ template <class ForwardIterator, class Predicate, class T>
+@@ -1700,7 +1701,7 @@
  
  	std::vector<char> pivot_used(pivots.size(), false); // can't be bool due to parallel write
  
@@ -45,9 +118,56 @@ diff -ur omptl.old/omptl_algorithm_par.h omptl/omptl_algorithm_par.h
  	assert(1u << max_depth <= P);
  	for (unsigned i = 0; i < max_depth; ++i)
  	{
+@@ -1781,7 +1782,7 @@
+ std::cout  << std::endl;
+ 
+ std::cout << borders.size() << " " << partitions.size() << " " << P << std::endl;
+-*/	
++*/
+ 	// Round one: sort final partitions, split remaining
+ 	#pragma omp parallel for
+ 	for (int i = 0; i < int(partitions.size()); ++i)
+@@ -1814,7 +1815,7 @@
+ 
+ 			const RandomAccessIterator begin  = partitions[i].first;
+ 			const RandomAccessIterator end    = partitions[i].second;
+-			
++
+ 			const RandomAccessIterator middle =
+ 				detail::_pivot_range(begin, end, pivots[pivot_index], comp);
+ 			partitions[i  ] = std::make_pair(begin, middle);
+diff -ur omptl.old/omptl_algorithm_ser.h omptl/omptl_algorithm_ser.h
+--- omptl.old/omptl_algorithm_ser.h	2022-06-19 08:21:39.815498672 +0200
++++ omptl/omptl_algorithm_ser.h	2022-06-19 08:21:52.960544697 +0200
+@@ -14,7 +14,7 @@
+ // You should have received a copy of the GNU Lesser General Public
+ // License along with this library; if not, write to the Free Software
+ // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+-
++#include <random>
+ 
+ namespace omptl
+ {
+@@ -463,14 +463,14 @@
+ void random_shuffle(RandomAccessIterator first, RandomAccessIterator last,
+       const unsigned P)
+ {
+-	return ::std::random_shuffle(first, last);
++	return ::std::shuffle(first, last, std::mt19937(std::random_device()()));
+ }
+ 
+ template <class RandomAccessIterator, class RandomNumberGenerator>
+ void random_shuffle(RandomAccessIterator first, RandomAccessIterator last,
+                     RandomNumberGenerator &rgen, const unsigned P)
+ {
+-	return ::std::random_shuffle(first, last, rgen);
++	return ::std::shuffle(first, last, std::mt19937(std::random_device()()));
+ }
+ 
+ template <class ForwardIterator, class T>
 diff -ur omptl.old/omptl_numeric omptl/omptl_numeric
---- omptl.old/omptl_numeric	2012-04-22 16:29:41.000000000 +0200
-+++ omptl/omptl_numeric	2021-06-20 15:40:29.000000000 +0200
+--- omptl.old/omptl_numeric	2022-06-19 08:21:39.816498675 +0200
++++ omptl/omptl_numeric	2022-06-19 08:21:52.955544679 +0200
 @@ -19,7 +19,7 @@
  #define OMPTL_NUMERIC 1
  
@@ -73,8 +193,8 @@ diff -ur omptl.old/omptl_numeric omptl/omptl_numeric
  
  #endif /* OMPTL_NUMERIC */
 diff -ur omptl.old/omptl_numeric_extensions.h omptl/omptl_numeric_extensions.h
---- omptl.old/omptl_numeric_extensions.h	2012-04-22 16:29:41.000000000 +0200
-+++ omptl/omptl_numeric_extensions.h	2021-06-20 15:40:29.000000000 +0200
+--- omptl.old/omptl_numeric_extensions.h	2022-06-19 08:21:39.815498672 +0200
++++ omptl/omptl_numeric_extensions.h	2022-06-19 08:21:52.956544683 +0200
 @@ -51,9 +51,9 @@
  } // namespace
  
@@ -88,8 +208,8 @@ diff -ur omptl.old/omptl_numeric_extensions.h omptl/omptl_numeric_extensions.h
  
  namespace omptl
 diff -ur omptl.old/omptl_numeric_par.h omptl/omptl_numeric_par.h
---- omptl.old/omptl_numeric_par.h	2012-04-22 16:29:41.000000000 +0200
-+++ omptl/omptl_numeric_par.h	2021-06-20 15:40:29.000000000 +0200
+--- omptl.old/omptl_numeric_par.h	2022-06-19 08:21:39.816498675 +0200
++++ omptl/omptl_numeric_par.h	2022-06-19 08:21:52.957544686 +0200
 @@ -23,8 +23,8 @@
  #include <functional>
  #include <iterator>
@@ -102,8 +222,8 @@ diff -ur omptl.old/omptl_numeric_par.h omptl/omptl_numeric_par.h
  namespace omptl
  {
 diff -ur omptl.old/omptl_tools.h omptl/omptl_tools.h
---- omptl.old/omptl_tools.h	2012-04-22 16:29:41.000000000 +0200
-+++ omptl/omptl_tools.h	2021-06-20 15:40:42.000000000 +0200
+--- omptl.old/omptl_tools.h	2022-06-19 08:21:39.816498675 +0200
++++ omptl/omptl_tools.h	2022-06-19 08:21:52.957544686 +0200
 @@ -25,7 +25,7 @@
  #include <climits>
  #include <iterator>

pyproject.toml

@@ -0,0 +1,4 @@
+[build-system]
+requires = ["setuptools","wheel","cython"]
+build-backend = "setuptools.build_meta"
+

python/CMakeLists.txt

@@ -2,7 +2,7 @@ set(CMAKE_SHARED_MODULE_PREFIX)
 
 
 set(PYTHON_INCLUDES ${NUMPY_INCLUDE_DIRS} ${PYTHON_INCLUDE_PATH} ${CMAKE_SOURCE_DIR}/python)
-include_directories(${CMAKE_SOURCE_DIR}/src ${CMAKE_BINARY_DIR}/src)
+include_directories(${CMAKE_SOURCE_DIR}/python ${CMAKE_SOURCE_DIR}/src ${CMAKE_BINARY_DIR}/src)
 
 IF(CYTHON)
   add_custom_command(
@@ -69,13 +69,10 @@ target_link_libraries(_fast_interp PRIVATE ${CosmoTool_local} )
 SET(ct_TARGETS   _cosmotool _project _cosmo_power _cosmo_cic _fast_interp _cosmomath)
 
 if (Boost_FOUND)
-  message(STATUS "Building bispectrum support (path = ${Boost_INCLUDE_DIRS})")
+  message(STATUS "Building bispectrum support")
   include_directories(${Boost_INCLUDE_DIRS})
   add_library(_cosmo_bispectrum MODULE _cosmo_bispectrum.cpp)
   target_link_libraries(_cosmo_bispectrum ${MATH_LIBRARY})
-  if(ENABLE_OPENMP)
-    set_target_properties(_cosmo_bispectrum PROPERTIES COMPILE_FLAGS "${OpenMP_CXX_FLAGS}" LINK_FLAGS "${OpenMP_CXX_FLAGS}")
-  endif()
   if (Boost_DEP)
     add_dependencies(_cosmo_bispectrum ${Boost_DEP})
   endif()

python/_cosmo_power.pyx

@@ -39,14 +39,10 @@ cdef extern from "cosmopower.hpp" namespace "CosmoTool":
     double OmegaEff
     double Gamma0
     double normPower
-    double A_BAO
-    double r_BAO
-    double k_D_BAO
     
         
     CosmoPower()
     void setFunction(CosmoFunction)
-    void setFunction_BAO(CosmoFunction,double,double,double)
     void updateCosmology()
     void updatePhysicalCosmology()
     void normalize(double,double)
@@ -106,6 +102,7 @@ cdef class CosmologyPower:
     self.power.SIGMA8 = s8
     self.power.normalize(k_min, k_max)
     
+    
   def setFunction(self,funcname):
     """setFunction(self, funcname)
     

python/_cosmomath.pyx

@@ -12,7 +12,7 @@ cdef extern from "numpy/npy_common.h":
    ctypedef npy_intp
 
 cdef extern from "special_math.hpp" namespace "CosmoTool":
-   T log_modified_bessel_first_kind[T](T v, T z) nogil except +
+   T log_modified_bessel_first_kind[T](T v, T z) except + nogil
 
 cdef extern from "numpy_adaptors.hpp" namespace "CosmoTool":
    void parallel_ufunc_dd_d[T,IT](char **args, IT* dimensions, IT* steps, void *func)

python/_cosmotool.pyx

@@ -39,7 +39,7 @@ cdef extern from "loadSimu.hpp" namespace "CosmoTool":
 
 cdef extern from "loadGadget.hpp" namespace "CosmoTool":
 
-  SimuData *loadGadgetMulti(const char *fname, int id, int flags, int gformat) nogil except +
+  SimuData *loadGadgetMulti(const char *fname, int id, int flags, int gformat) except + nogil
   void cxx_writeGadget "CosmoTool::writeGadget" (const char * s, SimuData *data) except +
 
 cdef extern from "safe_gadget.hpp":
@@ -313,6 +313,13 @@ tries to get an array from the object."""
 references to the object are gone. """
         pass
 
+
+cdef extern from "numpy/arrayobject.h":
+    # a little bit awkward: the reference to obj will be stolen
+    # using PyObject*  to signal that Cython cannot handle it automatically
+    int PyArray_SetBaseObject(np.ndarray arr, PyObject *obj) except -1 # -1 means there was an error
+
+
 cdef object wrap_array(void *p, np.uint64_t s, int typ):
     cdef np.ndarray ndarray
     cdef ArrayWrapper wrapper
@@ -320,7 +327,8 @@ cdef object wrap_array(void *p, np.uint64_t s, int typ):
     wrapper = ArrayWrapper()
     wrapper.set_data(s, typ, p)
     ndarray = np.array(wrapper, copy=False)
-    ndarray.base = <PyObject*> wrapper
+    #ndarray.base = <PyObject*> wrapper
+    PyArray_SetBaseObject(ndarray, <PyObject*> wrapper)
     Py_INCREF(wrapper)
 
     return ndarray

python/_project.pyx

@@ -727,7 +727,7 @@ cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density,
     cdef int iu0[3]
     cdef int i
     cdef int N = density.shape[0]
-    cdef int half_N = density.shape[0]/2
+    cdef int half_N = density.shape[0]//2
     cdef int completed
     cdef DTYPE_t I0, d, dist2, delta, s, max_distance2
     cdef int jumper[1]

requirements.txt

@@ -1,5 +1,4 @@
-numpy<1.19
+numpy<1.22
 cffi
 numexpr
-pyfftw
-cython
+cython<3

sample/CMakeLists.txt

@@ -1,4 +1,8 @@
 SET(tolink  ${CosmoTool_local} ${CosmoTool_LIBS} ${GSL_LIBRARIES} ${DL_LIBRARY})
+if (RT_LIBRARY)
+  SET(tolink ${tolink} ${RT_LIBRARY})
+endif()
+
 include_directories(${CMAKE_SOURCE_DIR}/src)
 include_directories(${FFTW3_INCLUDE_DIRS} ${EIGEN3_INCLUDE_DIRS} ${GSL_INCLUDE_PATH})
 if(YORICK_SUPPORT)

sample/simple3DFilter.cpp

@@ -84,6 +84,7 @@ int main(int argc, char **argv) {
   char *fname1, *outFile;
   double rLimit, boxsize, rLimit2, cx, cy, cz;
   int Nres;
+  int periodic;
 
   MiniArgDesc args[] = {{"INPUT DATA1", &fname1, MINIARG_STRING},
                         {"RADIUS LIMIT", &rLimit, MINIARG_DOUBLE},
@@ -93,6 +94,7 @@ int main(int argc, char **argv) {
                         {"CY", &cy, MINIARG_DOUBLE},
                         {"CZ", &cz, MINIARG_DOUBLE},
                         {"OUTPUT FILE", &outFile, MINIARG_STRING},
+                        {"PERIODIC", &periodic, MINIARG_INT},
                         {0, 0, MINIARG_NULL}};
 
   if (!parseMiniArgs(argc, argv, args))
@@ -148,6 +150,7 @@ int main(int argc, char **argv) {
 
   cout << "Building trees..." << endl;
   MyTree tree1(allCells_1, N1_points);
+  tree1.setPeriodic(periodic != 0, boxsize);
 
   cout << "Creating smoothing filter..." << endl;
 

setup.py

@@ -2,6 +2,7 @@ import stat
 import os
 import sys
 import shutil
+from sysconfig import get_config_var
 from distutils.command.install_data import install_data
 import pathlib
 from setuptools import find_packages, setup, Extension
@@ -167,17 +168,24 @@ class BuildCMakeExt(build_ext):
         # Change your cmake arguments below as necessary
         # Below is just an example set of arguments for building Blender as a Python module
 
+        c_compiler=os.environ.get('CC', get_config_var("CC"))
+        cxx_compiler=os.environ.get('CXX', get_config_var("CXX"))
+
         compilers=[]
         fill_up_settings=[
-              ("CMAKE_C_COMPILER", "CC"),
-              ("CMAKE_CXX_COMPILER", "CXX"),
-              ("CMAKE_C_FLAGS", "CFLAGS"),
-              ("CMAKE_EXE_LINKER_FLAGS_INIT", "LDFLAGS"),
-              ("CMAKE_SHARED_LINKER_FLAGS_INIT", "LDFLAGS"),
-              ("CMAKE_MODULE_LINKER_FLAGS_INIT", "LDFLAGS")]
-        for cmake_flag, os_flag in fill_up_settings:
-          if os_flag in os.environ:
-             compilers.append(f"-D{cmake_flag}={os.environ[os_flag]}")
+              ("CMAKE_C_COMPILER", "CC", get_config_var("CC")),
+              ("CMAKE_CXX_COMPILER", "CXX", get_config_var("CXX")),
+              ("CMAKE_C_FLAGS", "CFLAGS", None),
+              ("CMAKE_CXX_FLAGS", "CXXFLAGS", None),
+              ("CMAKE_EXE_LINKER_FLAGS_INIT", "LDFLAGS", None),
+              ("CMAKE_SHARED_LINKER_FLAGS_INIT", "LDFLAGS", None),
+              ("CMAKE_MODULE_LINKER_FLAGS_INIT", "LDFLAGS", None)]
+        for cmake_flag, os_flag, default_flag in fill_up_settings:
+          if os_flag in os.environ or default_flag is not None:
+             c = os.environ.get(os_flag, default_flag)
+             if not c is None:
+               compilers.append(f"-D{cmake_flag}={c}")
+        print(compilers)
 
         self.spawn(['cmake', '-H'+SOURCE_DIR, '-B'+self.build_temp,
                     '-DENABLE_OPENMP=ON','-DINTERNAL_BOOST=ON','-DINTERNAL_EIGEN=ON',
@@ -229,10 +237,10 @@ class BuildCMakeExt(build_ext):
 CosmoTool_extension = CMakeExtension(name="cosmotool")
 
 setup(name='cosmotool',
-      version='1.3.1',
+      version='1.3.6',
       packages=["cosmotool"],
       package_dir={'cosmotool': 'python/cosmotool'},
-      install_requires=['numpy','cffi','numexpr','pyfftw','h5py'],
+      install_requires=['cython','numpy','cffi','numexpr','h5py'],
       setup_requires=['cython','cffi','numpy','numexpr'],
       ext_modules=[CosmoTool_extension],
       description='A small cosmotool box of useful functions',

src/cosmopower.cpp

@@ -212,6 +212,30 @@ double CosmoPower::powerHuWiggles(double k)
   return normPower * pow(k,n) * T_k * T_k;
 }
 
+double CosmoPower::BAO_Tk(double k){
+
+  double no_wiggle_tk = noWiggleTk(k);
+
+  double A = 0;
+  double r_s = 10;
+  double k_D = 2 * M_PI / 100;
+
+  //sqrt as we want to make the parameterization part of the transfer function
+  double param = sqrt(1 + A * sin(k * r_s) * exp(- k / k_D)); 
+  return no_wiggle_tk * param;
+
+}
+
+double CosmoPower::sample_BAO(double k)
+{
+  // BAO wiggle parameterization for reconstruction
+  // Babic et al. 2022, https://arxiv.org/abs/2203.06177
+
+  double T_k = BAO_Tk(k); 
+
+  return normPower * pow(k,n) * T_k * T_k;
+  }
+
 double CosmoPower::primordialPowerSpectrum(double k)
 {
   //Primordial power spectrum, needed for PNG
@@ -457,14 +481,18 @@ void CosmoPower::setFunction(CosmoFunction f)
     case POWER_SUGIYAMA:
       eval = &CosmoPower::powerSugiyama;
       break;
+    case SAMPLE_WIGGLES:
+      eval = &CosmoPower::sample_BAO;
+      break;
     case POWER_BDM:
       eval = &CosmoPower::powerBDM;
       break;
-    case NOWIGGLE_TK:
-      eval = &CosmoPower::noWiggleTk;
     case POWER_TEST:
       eval = &CosmoPower::powerTest;
       break;
+    case BAO_TK:
+      eval = &CosmoPower::BAO_Tk;
+      break;
     default:
       abort();
     }

src/cosmopower.hpp

@@ -66,10 +66,6 @@ namespace CosmoTool {
     double Gamma0;
     double normPower;
 
-    double A_BAO;
-    double r_BAO;
-    double k_D_BAO;
-
     struct EHuParams {
       double k_silk;
       double s;
@@ -94,13 +90,15 @@ namespace CosmoTool {
       POWER_BDM,
       POWER_TEST,
       HU_WIGGLES_ORIGINAL,
-      NOWIGGLE_TK
+      SAMPLE_WIGGLES,
+      BAO_TK
     };
 
     CosmoPower();
     ~CosmoPower();
 
     void setFunction(CosmoFunction f);
+
     void updateCosmology();
     void updatePhysicalCosmology();
     void normalize(double k_min = -1, double k_max = -1);
@@ -112,7 +110,6 @@ namespace CosmoTool {
     double power(double k);
 
     double integrandNormalize(double k);
-
   private:
     double (CosmoPower::*eval)(double);
 
@@ -127,8 +124,9 @@ namespace CosmoTool {
     double powerBDM(double k);
     double powerTest(double k);
     double powerHuWigglesOriginal(double k);
+    double sample_BAO(double k);
     double noWiggleTk(double k);
-
+    double BAO_Tk(double k);
   };
 
 };

src/h5_readFlash.cpp

@@ -123,10 +123,7 @@ void h5_read_runtime_parameters
   DataSpace memspace(rank, &dimens_1d);
   //memspace = H5Screate_simple(rank, &dimens_1d, NULL);
 
-  dataset.read(int_rt_parms, int_rt_type, memspace, dataspace,
-		    H5P_DEFAULT);
-  //status = H5Dread(dataset, int_rt_type, memspace, dataspace,
-	//	    H5P_DEFAULT, int_rt_parms);
+  dataset.read(int_rt_parms, int_rt_type, memspace, dataspace);
 
 
   for (i = 0; i < nint_runtime_parameters; i++) {
@@ -167,10 +164,7 @@ void h5_read_runtime_parameters
    memspace = DataSpace(rank, &dimens_1d);
   //memspace = H5Screate_simple(rank, &dimens_1d, NULL);
 
-  dataset.read(real_rt_parms, real_rt_type, memspace, dataspace,
-		    H5P_DEFAULT);
-  //status = H5Dread(dataset, real_rt_type, memspace, dataspace,
-//		    H5P_DEFAULT, real_rt_parms);
+  dataset.read(real_rt_parms, real_rt_type, memspace, dataspace);
 
 
   for (i = 0; i < nreal_runtime_parameters; i++) {
@@ -289,9 +283,7 @@ void h5_read_flash3_particles (H5File* file,
   DataSpace memspace(rank, dimens_2d);
   //memspace = H5Screate_simple(rank, dimens_2d, NULL);
 
-  dataset.read(part_names, string_type, H5S_ALL, H5S_ALL, H5P_DEFAULT);
-  //status = H5Dread(dataset, string_type, H5S_ALL, H5S_ALL,
-  //              H5P_DEFAULT, part_names);
+  dataset.read(part_names, string_type);
 
 
   string_type.close();
@@ -374,8 +366,7 @@ void h5_read_flash3_particles (H5File* file,
     //memspace   = H5Screate_simple(rank, dimens_2d, maxdimens_2d);
 
     /* read data from the dataset */
-   dataset.read(partBuffer, datatype, memspace, dataspace, H5P_DEFAULT);
-   //status = H5Dread(dataset, datatype, memspace, dataspace, H5P_DEFAULT, partBuffer);
+   dataset.read(partBuffer, datatype, memspace, dataspace);
 
     /* convert buffer into particle struct */
 
@@ -485,8 +476,7 @@ void h5_read_flash3_header_info(H5File* file,
           PredType::NATIVE_DOUBLE);
 
   // read the data into 'real_list'
-  dataset.read( real_list, real_list_type, memspace, dataspace,
-                H5P_DEFAULT);
+  dataset.read( real_list, real_list_type, memspace, dataspace);
 
 
   if (status < 0) {

src/hdf5_array.hpp

@@ -426,7 +426,7 @@ namespace CosmoTool {
 #define CTOOL_HDF5_INSERT_ELEMENT(r, STRUCT, element) \
      { \
        ::CosmoTool::get_hdf5_data_type<BOOST_PP_TUPLE_ELEM(2, 0, element)> t; \
-       long position = HOFFSET(STRUCT, BOOST_PP_TUPLE_ELEM(2, 1, element)); \
+       long position = offsetof(STRUCT, BOOST_PP_TUPLE_ELEM(2, 1, element)); \
        const char *field_name = BOOST_PP_STRINGIZE(BOOST_PP_TUPLE_ELEM(2, 1, element)); \
        type.insertMember(field_name,  position, t.type()); \
      } 

src/interpolate.cpp

@@ -157,6 +157,7 @@ const Interpolate& Interpolate::operator=(const Interpolate& a)
   gsl_spline_init(spline, a.spline->x, a.spline->y, a.spline->size);
   logx = a.logx;
   logy = a.logy;
+  return *this;
 }
 
 double Interpolate::getMaxX() const

src/mykdtree.hpp

@@ -1,5 +1,5 @@
 /*+
-This is CosmoTool (./src/mykdtree.hpp) -- Copyright (C) Guilhem Lavaux (2007-2014)
+This is CosmoTool (./src/mykdtree.hpp) -- Copyright (C) Guilhem Lavaux (2007-2022)
 
 guilhem.lavaux@gmail.com
 

src/mykdtree.tcc

@@ -643,8 +643,8 @@ namespace CosmoTool {
            throw InvalidOnDiskKDTree();
         }
         if (node_on_disk.cell_id > numNodes || node_on_disk.cell_id < 0 ||
-	    node_on_disk.children_node[0] > lastNode || node_on_disk.children_node[0] < -1 ||
-	    node_on_disk.children_node[1] > lastNode || node_on_disk.children_node[1] < -1)
+	    (node_on_disk.children_node[0] >= 0 && node_on_disk.children_node[0] > lastNode) || node_on_disk.children_node[0] < -1 ||
+	    (node_on_disk.children_node[1] >= 0 && node_on_disk.children_node[1] > lastNode) || node_on_disk.children_node[1] < -1)
 	  {
 	    delete[] nodes;
             std::cerr << "Invalid cell id or children node id invalid" << std::endl;

src/sphSmooth.hpp

@@ -1,5 +1,5 @@
 /*+
-This is CosmoTool (./src/sphSmooth.hpp) -- Copyright (C) Guilhem Lavaux (2007-2014)
+This is CosmoTool (./src/sphSmooth.hpp) -- Copyright (C) Guilhem Lavaux (2007-2022)
 
 guilhem.lavaux@gmail.com
 
@@ -39,14 +39,11 @@ knowledge of the CeCILL license and that you accept its terms.
 #include "config.hpp"
 #include "mykdtree.hpp"
 
-namespace CosmoTool
-{
+namespace CosmoTool {
   template <typename ValType, int Ndims = NUMDIMS>
-  class SPHSmooth
-  {
+  class SPHSmooth {
   public:
-    typedef struct
-    {
+    typedef struct {
       ComputePrecision weight;
       ValType pValue;
     } FullType;
@@ -61,8 +58,7 @@ namespace CosmoTool
 
   public:
     typedef SPHCell *P_SPHCell;
-    struct SPHState
-    {
+    struct SPHState {
       boost::shared_ptr<P_SPHCell[]> ngb;
       boost::shared_ptr<CoordType[]> distances;
       typename SPHTree::coords currentCenter;
@@ -70,44 +66,57 @@ namespace CosmoTool
       ComputePrecision smoothRadius;
     };
 
-  
     SPHSmooth(SPHTree *tree, uint32_t Nsph);
     virtual ~SPHSmooth();
 
-    void fetchNeighbours(const typename SPHTree::coords& c, SPHState *state = 0);
+    void
+    fetchNeighbours(const typename SPHTree::coords &c, SPHState *state = 0);
 
-    void fetchNeighbours(const typename SPHTree::coords& c, uint32_t newNsph);
-    void fetchNeighboursOnVolume(const typename SPHTree::coords& c, ComputePrecision radius);
-    const typename SPHTree::coords& getCurrentCenter() const
-    {
+    void fetchNeighbours(
+        const typename SPHTree::coords &c, uint32_t newNsph,
+        SPHState *state = 0);
+    void fetchNeighboursOnVolume(
+        const typename SPHTree::coords &c, ComputePrecision radius);
+    const typename SPHTree::coords &getCurrentCenter() const {
       return internal.currentCenter;
     }
 
+    /** This is the pure SPH smoothing function. It does not reweight by the
+     * value computed at each grid site.
+     */
     template <typename FuncT>
-    ComputePrecision computeSmoothedValue(const typename SPHTree::coords& c,
-					  FuncT fun, SPHState *state = 0);    
+    ComputePrecision computeSmoothedValue(
+        const typename SPHTree::coords &c, FuncT fun, SPHState *state = 0);
 
+    /** This is the weighted SPH smoothing function. It does reweight by the
+     * value computed at each grid site. This ensures the total sum of the interpolated
+     * quantity is preserved by interpolating to the target mesh.
+     */
     template <typename FuncT>
-    ComputePrecision computeInterpolatedValue(const typename SPHTree::coords& c,
-					      FuncT fun, SPHState *state = 0);
+    ComputePrecision computeInterpolatedValue(
+        const typename SPHTree::coords &c, FuncT fun, SPHState *state = 0);
 
-    ComputePrecision getMaxDistance(const typename SPHTree::coords& c, 
-				    SPHNode *node) const;
+    /** This is the adjoint gradient of computeInterpolatedValue w.r.t. to the value
+     * array. FuncT is expected to have the following prototype:
+     *    void((CellValue defined by the user), ComputePrecision weighted_ag_value)
+     */
+    template <typename FuncT>
+    void computeAdjointGradientSmoothedValue(
+        const typename SPHTree::coords &c, ComputePrecision ag_value, FuncT fun,
+        SPHState *state = 0);
 
-    ComputePrecision getSmoothingLen() const
-    {
-      return internal.smoothRadius;
-    }
+    ComputePrecision
+    getMaxDistance(const typename SPHTree::coords &c, SPHNode *node) const;
+
+    ComputePrecision getSmoothingLen() const { return internal.smoothRadius; }
 
     // TO USE WITH EXTREME CARE !
-    void setSmoothingLen(ComputePrecision len)
-    {
-      internal.smoothRadius = len;
-    }
+    void setSmoothingLen(ComputePrecision len) { internal.smoothRadius = len; }
     // END
 
     template <typename FuncT>
     void runForEachNeighbour(FuncT fun, SPHState *state = 0);
+    void addGridSite(const typename SPHTree::coords &c, SPHState *state);
     void addGridSite(const typename SPHTree::coords &c);
 
     bool hasNeighbours() const;
@@ -134,20 +143,20 @@ namespace CosmoTool
     SPHTree *tree;
 
     template <typename FuncT>
-    ComputePrecision computeWValue(const typename SPHTree::coords & c,
-				   SPHCell& cell,
-				   CoordType d,
+    ComputePrecision computeWValue(
+        const typename SPHTree::coords &c, SPHCell &cell, CoordType d,
         FuncT fun, SPHState *state);
 
     template <typename FuncT>
-    void runUnrollNode(SPHNode *node,
-		       FuncT fun);
+    void runUnrollNode(SPHNode *node, FuncT fun);
   };
 
   template <class ValType1, class ValType2, int Ndims>
-  bool operator<(const SPHSmooth<ValType1, Ndims>& s1, const SPHSmooth<ValType2, Ndims>& s2);
+  bool operator<(
+      const SPHSmooth<ValType1, Ndims> &s1,
+      const SPHSmooth<ValType2, Ndims> &s2);
 
-};
+}; // namespace CosmoTool
 
 #include "sphSmooth.tcc"
 

src/sphSmooth.tcc

@@ -4,8 +4,7 @@
 namespace CosmoTool {
 
   template <typename ValType, int Ndims>
-SPHSmooth<ValType,Ndims>::SPHSmooth(SPHTree *tree, uint32_t Nsph)
-{
+  SPHSmooth<ValType, Ndims>::SPHSmooth(SPHTree *tree, uint32_t Nsph) {
     this->Nsph = Nsph;
     this->tree = tree;
     internal.currentNgb = 0;
@@ -16,17 +15,13 @@ SPHSmooth<ValType,Ndims>::SPHSmooth(SPHTree *tree, uint32_t Nsph)
   }
 
   template <typename ValType, int Ndims>
-SPHSmooth<ValType,Ndims>::~SPHSmooth()
-{
-}
+  SPHSmooth<ValType, Ndims>::~SPHSmooth() {}
 
   template <typename ValType, int Ndims>
   template <typename FuncT>
-ComputePrecision SPHSmooth<ValType,Ndims>::computeWValue(const typename SPHTree::coords& c,
-							 SPHCell& cell,
-							 CoordType d,
-							 FuncT fun, SPHState *state)
-{
+  ComputePrecision SPHSmooth<ValType, Ndims>::computeWValue(
+      const typename SPHTree::coords &c, SPHCell &cell, CoordType d, FuncT fun,
+      SPHState *state) {
     CoordType weight;
 
     d /= state->smoothRadius;
@@ -39,37 +34,44 @@ ComputePrecision SPHSmooth<ValType,Ndims>::computeWValue(const typename SPHTree:
   }
 
   template <typename ValType, int Ndims>
-void
-SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c, uint32_t newNngb)
-{
+  void SPHSmooth<ValType, Ndims>::fetchNeighbours(
+      const typename SPHTree::coords &c, uint32_t newNngb, SPHState *state) {
     ComputePrecision d2, max_dist = 0;
     uint32_t requested = newNngb;
 
-  if (requested > maxNgb)
-    {
+    if (state != 0) {
+      state->distances = boost::shared_ptr<CoordType[]>(new CoordType[newNngb]);
+      state->ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[newNngb]);
+    } else {
+      state = &internal;
+      if (requested > maxNgb) {
         maxNgb = requested;
         internal.ngb = boost::shared_ptr<P_SPHCell[]>(new P_SPHCell[maxNgb]);
-      internal.distances = boost::shared_ptr<CoordType[]>(new CoordType[maxNgb]);
+        internal.distances =
+            boost::shared_ptr<CoordType[]>(new CoordType[maxNgb]);
+      }
     }
 
-  memcpy(internal.currentCenter, c, sizeof(c));
-  tree->getNearestNeighbours(c, requested, (SPHCell **)internal.ngb.get(), (CoordType*)internal.distances.get());
+    memcpy(state->currentCenter, c, sizeof(c));
+    tree->getNearestNeighbours(
+        c, requested, (SPHCell **)state->ngb.get(),
+        (CoordType *)state->distances.get());
 
-  internal.currentNgb = 0;
-  for (uint32_t i = 0; i < requested && (internal.ngb)[i] != 0; i++,internal.currentNgb++)
-    {
-      internal.distances[i] = sqrt(internal.distances[i]);
-      d2 = internal.distances[i];
+    state->currentNgb = 0;
+    for (uint32_t i = 0; i < requested && (state->ngb)[i] != 0;
+         i++, state->currentNgb++) {
+      state->distances[i] = sqrt(state->distances[i]);
+      d2 = state->distances[i];
       if (d2 > max_dist)
         max_dist = d2;
     }
 
-  internal.smoothRadius = max_dist / 2;  
+    state->smoothRadius = max_dist / 2;
   }
 
   template <typename ValType, int Ndims>
-void SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c, SPHState *state)
-{
+  void SPHSmooth<ValType, Ndims>::fetchNeighbours(
+      const typename SPHTree::coords &c, SPHState *state) {
     ComputePrecision d2, max_dist = 0;
     uint32_t requested = Nsph;
 
@@ -81,12 +83,13 @@ void SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c
 
     memcpy(state->currentCenter, c, sizeof(c));
 
-  tree->getNearestNeighbours(c, requested, state->ngb.get(), state->distances.get());
+    tree->getNearestNeighbours(
+        c, requested, state->ngb.get(), state->distances.get());
 
     state->currentNgb = 0;
-  for (uint32_t i = 0; i < requested && state->ngb[i] != 0; i++,state->currentNgb++)
-    {
-      d2 = internal.distances[i] = sqrt(internal.distances[i]);
+    for (uint32_t i = 0; i < requested && state->ngb[i] != 0;
+         i++, state->currentNgb++) {
+      d2 = state->distances[i] = sqrt(state->distances[i]);
       if (d2 > max_dist)
         max_dist = d2;
     }
@@ -94,22 +97,18 @@ void SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c
     state->smoothRadius = max_dist / 2;
   }
 
-
   template <typename ValType, int Ndims>
-void
-SPHSmooth<ValType,Ndims>::fetchNeighboursOnVolume(const typename SPHTree::coords& c,
-                                                  ComputePrecision radius)
-{
+  void SPHSmooth<ValType, Ndims>::fetchNeighboursOnVolume(
+      const typename SPHTree::coords &c, ComputePrecision radius) {
     uint32_t numPart;
     ComputePrecision d2, max_dist = 0;
 
     memcpy(internal.currentCenter, c, sizeof(c));
 
-  internal.currentNgb = tree->getIntersection(c, radius, internal.ngb, internal.distances,
-				  maxNgb);
+    internal.currentNgb = tree->getIntersection(
+        c, radius, internal.ngb, internal.distances, maxNgb);
 
-  for (uint32_t i = 0; i < internal.currentNgb; i++)
-    {
+    for (uint32_t i = 0; i < internal.currentNgb; i++) {
       d2 = internal.distances[i] = sqrt(internal.distances[i]);
       if (d2 > max_dist)
         max_dist = d2;
@@ -119,10 +118,8 @@ SPHSmooth<ValType,Ndims>::fetchNeighboursOnVolume(const typename SPHTree::coords
 
   template <typename ValType, int Ndims>
   template <typename FuncT>
-ComputePrecision 
-SPHSmooth<ValType,Ndims>::computeSmoothedValue(const typename SPHTree::coords& c,
-                                               FuncT fun, SPHState *state)
-{
+  ComputePrecision SPHSmooth<ValType, Ndims>::computeSmoothedValue(
+      const typename SPHTree::coords &c, FuncT fun, SPHState *state) {
     if (state == 0)
       state = &internal;
 
@@ -130,26 +127,24 @@ SPHSmooth<ValType,Ndims>::computeSmoothedValue(const typename SPHTree::coords& c
     ComputePrecision max_dist = 0;
     ComputePrecision r3 = cube(state->smoothRadius);
 
-  for (uint32_t i = 0; i < state->currentNgb; i++)
-    {
-      outputValue += computeWValue(c, *state->ngb[i], state->distances[i], fun, state);
+    for (uint32_t i = 0; i < state->currentNgb; i++) {
+      outputValue +=
+          computeWValue(c, *state->ngb[i], state->distances[i], fun, state);
     }
 
     return outputValue / r3;
   }
 
   template <typename ValType>
-ComputePrecision interpolateOne(const ValType& t)
-{
+  ComputePrecision interpolateOne(const ValType &t) {
     return 1.0;
   }
 
-// WARNING ! Cell's weight must be 1 !!!
   template <typename ValType, int Ndims>
   template <typename FuncT>
-ComputePrecision SPHSmooth<ValType,Ndims>::computeInterpolatedValue(const typename SPHTree::coords& c,
-								    FuncT fun, SPHState *state)
-{
+  void SPHSmooth<ValType, Ndims>::computeAdjointGradientSmoothedValue(
+      const typename SPHTree::coords &c, ComputePrecision ag_value, FuncT fun,
+      SPHState *state) {
     if (state == 0)
       state = &internal;
 
@@ -157,10 +152,38 @@ ComputePrecision SPHSmooth<ValType,Ndims>::computeInterpolatedValue(const typena
     ComputePrecision max_dist = 0;
     ComputePrecision weight = 0;
 
-  for (uint32_t i = 0; i < state->currentNgb; i++)
-    {
+    for (uint32_t i = 0; i < state->currentNgb; i++) {
+      weight +=
+          computeWValue(c, *state->ngb[i], state->distances[i], interpolateOne);
+    }
+
+    for (uint32_t i = 0; i < state->currentNgb; i++) {
+      auto &cell = *state->ngb[i];
+      double partial_ag =
+          computeWValue(
+              c, cell, state->distances[i],
+              [ag_value](ComputePrecision) { return ag_value; }) /
+          weight;
+      fun(cell.val.pValue, ag_value);
+    }
+  }
+
+  // WARNING ! Cell's weight must be 1 !!!
+  template <typename ValType, int Ndims>
+  template <typename FuncT>
+  ComputePrecision SPHSmooth<ValType, Ndims>::computeInterpolatedValue(
+      const typename SPHTree::coords &c, FuncT fun, SPHState *state) {
+    if (state == 0)
+      state = &internal;
+
+    ComputePrecision outputValue = 0;
+    ComputePrecision max_dist = 0;
+    ComputePrecision weight = 0;
+
+    for (uint32_t i = 0; i < state->currentNgb; i++) {
       outputValue += computeWValue(c, *state->ngb[i], state->distances[i], fun);
-      weight += computeWValue(c, *state->ngb[i], state->distances[i], interpolateOne);
+      weight +=
+          computeWValue(c, *state->ngb[i], state->distances[i], interpolateOne);
     }
 
     return (outputValue == 0) ? 0 : (outputValue / weight);
@@ -168,62 +191,60 @@ ComputePrecision SPHSmooth<ValType,Ndims>::computeInterpolatedValue(const typena
 
   template <typename ValType, int Ndims>
   template <typename FuncT>
-void SPHSmooth<ValType,Ndims>::runForEachNeighbour(FuncT fun, SPHState *state)
-{
+  void
+  SPHSmooth<ValType, Ndims>::runForEachNeighbour(FuncT fun, SPHState *state) {
     if (state == 0)
       state = &internal;
 
-  for (uint32_t i = 0; i < state->currentNgb; i++)
-    {
+    for (uint32_t i = 0; i < state->currentNgb; i++) {
       fun(state->ngb[i]);
     }
   }
 
-
   template <typename ValType, int Ndims>
-void SPHSmooth<ValType,Ndims>::addGridSite(const typename SPHTree::coords& c)
-{
+  void SPHSmooth<ValType, Ndims>::addGridSite(
+      const typename SPHTree::coords &c, SPHState *state) {
     ComputePrecision outputValue = 0;
     ComputePrecision max_dist = 0;
+    ComputePrecision r3 = cube(state->smoothRadius);
 
-  ComputePrecision r3 = cube(internal.smoothRadius);
-
-  for (uint32_t i = 0; i < internal.currentNgb; i++)
-    {
-      ComputePrecision d = internal.distances[i];
-      SPHCell& cell = *(internal.ngb[i]);
-      double kernel_value = getKernel(d/internal.smoothRadius) / r3;
-#pragma omp atomic
+    for (uint32_t i = 0; i < state->currentNgb; i++) {
+      ComputePrecision d = state->distances[i];
+      SPHCell &cell = *(state->ngb[i]);
+      double kernel_value = getKernel(d / state->smoothRadius) / r3;
+#pragma omp atomic update
       cell.val.weight += kernel_value;
     }
   }
 
+  template <typename ValType, int Ndims>
+  void
+  SPHSmooth<ValType, Ndims>::addGridSite(const typename SPHTree::coords &c) {
+    addGridSite(c, &internal);
+  }
+
   template <typename ValType, int Ndims>
   ComputePrecision
-SPHSmooth<ValType,Ndims>::getKernel(ComputePrecision x) const
-{
+  SPHSmooth<ValType, Ndims>::getKernel(ComputePrecision x) const {
     // WARNING !!! This is an unnormalized version of the kernel.
     if (x < 1)
       return 1 - 1.5 * x * x + 0.75 * x * x * x;
-  else if (x < 2)
-    {
+    else if (x < 2) {
       ComputePrecision d = 2 - x;
       return 0.25 * d * d * d;
-    }
-  else
+    } else
       return 0;
   }
 
   template <typename ValType, int Ndims>
-bool SPHSmooth<ValType,Ndims>::hasNeighbours() const
-{
+  bool SPHSmooth<ValType, Ndims>::hasNeighbours() const {
     return (internal.currentNgb != 0);
   }
 
   template <class ValType1, class ValType2, int Ndims>
-bool operator<(const SPHSmooth<ValType1, Ndims>& s1, const SPHSmooth<ValType2, Ndims>& s2)
-{
+  bool operator<(
+      const SPHSmooth<ValType1, Ndims> &s1,
+      const SPHSmooth<ValType2, Ndims> &s2) {
     return (s1.getSmoothingLen() < s2.getSmoothingLen());
   }
-
-};
+}; // namespace CosmoTool