Fixup for sph

Add possibility to have an external state
Remove support for old python
2024-04-26 08:53:38 +02:00 · 2024-04-25 18:11:06 +02:00 · 2024-04-25 17:06:15 +02:00 · 2024-02-13 10:47:28 +01:00 · 2024-02-13 10:37:52 +01:00 · 2024-02-13 10:28:58 +01:00
30 changed files with 4431 additions and 534 deletions
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@ -10,8 +10,6 @@ include(FindPkgConfig)
 include(FindPackageHandleStandardArgs)
 include(color_msg)

-find_package(MPI)
-
 option(BUILD_SHARED_LIBS "Build shared libraries." OFF)
 option(BUILD_STATIC_LIBS "Build static libraries." ON)
 option(ENABLE_SHARP "Enable SHARP support." ON)
@ -44,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)
@ -57,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)
@ -72,14 +78,12 @@ SET(CPACK_PACKAGE_DESCRIPTION_SUMMARY "A toolbox for impatient cosmologists")
 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 "2")
-SET(CPACK_PACKAGE_VERSION_PATCH "3${EXTRA_VERSION}")
+SET(CPACK_PACKAGE_VERSION_MINOR "3")
+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 "/CVS/;/\\\\.git/;/\\\\.svn/;\\\\.swp$;\\\\.#;/#;.*~;cscope.*;/CMakeFiles/;.*\\\\.cmake;Makefile")
-
-include_directories( ${MPI_C_INCLUDE_PATH})
-
+SET(CPACK_SOURCE_IGNORE_FILES 
+"/CVS/;/\\\\.git/;/\\\\.svn/;\\\\.swp$;\\\\.#;/#;.*~;cscope.*;/CMakeFiles/;.*\\\\.cmake;Makefile")

 add_subdirectory(src)
 add_subdirectory(sample)
--- a/MANIFEST.in
+++ b/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
--- a/builder/build-wheels.sh
+++ b/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
--- a/builder/start.sh
+++ b/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
--- a/conda/meta.yaml
+++ b/conda/meta.yaml
@ -1,6 +1,6 @@
 package:
  name: cosmotool
-  version: "1.2.3"
+  version: "1.3.0"

 source:
  git_rev: a86c9a8
--- a/external/config.guess
+++ b/external/config.guess
--- a/external/config.sub
+++ b/external/config.sub
--- a/external/external_build.cmake
+++ b/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")
+    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)

-  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(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,43 +168,52 @@ 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
+    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
-    -DCMAKE_POSITION_INDEPENDENT_CODE=ON
-    -DENABLE_DAP=OFF
+                -DENABLE_BYTERANGE=FALSE
+		-DCMAKE_POSITION_INDEPENDENT_CODE=ON
+		-DENABLE_DAP=OFF
 		-DCMAKE_INSTALL_PREFIX=${NETCDF_BIN_DIR}
                -DHDF5_C_LIBRARY=${HDF5_C_STATIC_LIBRARY}
                -DHDF5_HL_LIBRARY=${HDF5_HL_STATIC_LIBRARY}
                -DHDF5_INCLUDE_DIR=${HDF5_INCLUDE_DIRS}
-    -DCMAKE_INSTALL_LIBDIR=lib
+		-DCMAKE_INSTALL_LIBDIR=lib
  )
  
  SET(NETCDFCXX_SOURCE_DIR ${BUILD_PREFIX}/netcdf-c++-prefix/src/netcdf-c++)
  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
+                -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"
--- a/external/libsharp-8d51946.tar.gz
+++ b/external/libsharp-8d51946.tar.gz
--- a/external/patch-omptl
+++ b/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>
--- a/pyproject.toml
+++ b/pyproject.toml
@ -0,0 +1,4 @@
+[build-system]
+requires = ["setuptools","wheel","cython"]
+build-backend = "setuptools.build_meta"
+
--- a/python/CMakeLists.txt
+++ b/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()
--- a/python/_cosmomath.pyx
+++ b/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)
--- a/python/_cosmotool.pyx
+++ b/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
--- a/python/_project.pyx
+++ b/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]
--- a/requirements.txt
+++ b/requirements.txt
@ -1,5 +1,4 @@
-numpy<1.19
+numpy<1.22
 cffi
 numexpr
-pyfftw
-cython
+cython<3
--- a/sample/CMakeLists.txt
+++ b/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)
@ -87,7 +91,7 @@ if (Boost_FOUND)
  ENDIF (YORICK_SUPPORT)
  if (HDF5_FOUND)
    add_executable(simple3DFilter simple3DFilter.cpp)
-    target_link_libraries(simple3DFilter ${tolink}  ${MPI_C_LIBRARIES})
+    target_link_libraries(simple3DFilter ${tolink})

    add_executable(simpleDistanceFilter simpleDistanceFilter.cpp)
    target_link_libraries(simpleDistanceFilter ${tolink})
@ -101,7 +105,7 @@ if (Boost_FOUND)

    add_executable(gadgetToArray gadgetToArray.cpp)
    target_link_libraries(gadgetToArray ${tolink})
-
+    
    add_executable(testHDF5 testHDF5.cpp)
    target_link_libraries(testHDF5 ${tolink})

--- a/sample/simple3DFilter.cpp
+++ b/sample/simple3DFilter.cpp
@ -1,11 +1,3 @@
-#include <H5Cpp.h>
-#include <mpi.h>
-
-#include <boost/bind.hpp>
-#include <boost/format.hpp>
-#include <cassert>
-#include <iostream>
-
 #include "hdf5_array.hpp"
 #include "miniargs.hpp"
 #include "mykdtree.hpp"
@ -13,6 +5,11 @@
 #include "openmp.hpp"
 #include "sphSmooth.hpp"
 #include "yorick.hpp"
+#include <H5Cpp.h>
+#include <boost/bind.hpp>
+#include <boost/format.hpp>
+#include <cassert>
+#include <iostream>

 using namespace std;
 using namespace CosmoTool;
@ -30,32 +27,37 @@ typedef boost::multi_array<float, 2> array_type;
 typedef boost::multi_array<float, 3> array3_type;
 typedef boost::multi_array<float, 4> array4_type;

-ComputePrecision getVelocity(const VCoord& v, int i) { return v.mass * v.v[i]; }
+ComputePrecision getVelocity(const VCoord &v, int i) { return v.mass * v.v[i]; }

-ComputePrecision getMass(const VCoord& v) { return v.mass; }
+ComputePrecision getMass(const VCoord &v) { return v.mass; }

 typedef SPHSmooth<VCoord> MySmooth;
 typedef MySmooth::SPHTree MyTree;
 typedef MyTree::Cell MyCell;

 template <typename FuncT>
-void computeInterpolatedField(MyTree* tree1, double boxsize, int Nres,
+void computeInterpolatedField(MyTree *tree1, double boxsize, int Nres,
                              double cx, double cy, double cz,
-                              array3_type& bins, array3_type& arr, FuncT func,
+                              array3_type &bins, array3_type &arr, FuncT func,
                              double rLimit2) {
-#pragma omp parallel
+  int rz_max = 0;
+#pragma omp parallel shared(rz_max)
  {
    MySmooth smooth1(tree1, N_SPH);

-#pragma omp for schedule(dynamic)
+#pragma omp for collapse(3) schedule(dynamic)
    for (int rz = 0; rz < Nres; rz++) {
-      double pz = (rz)*boxsize / Nres - cz;

-      cout << format("[%d] %d / %d") % smp_get_thread_id() % rz % Nres << endl;
      for (int ry = 0; ry < Nres; ry++) {
-        double py = (ry)*boxsize / Nres - cy;
        for (int rx = 0; rx < Nres; rx++) {
+          if (rz > rz_max) {
+            rz_max = rz;
+            cout << format("[%d] %d / %d") % smp_get_thread_id() % rz % Nres
+                 << endl;
+          }
          double px = (rx)*boxsize / Nres - cx;
+          double py = (ry)*boxsize / Nres - cy;
+          double pz = (rz)*boxsize / Nres - cz;

          MyTree::coords c = {float(px), float(py), float(pz)};

@ -78,18 +80,11 @@ void computeInterpolatedField(MyTree* tree1, double boxsize, int Nres,
  }
 }

-int main(int argc, char** argv) {
-  int provided;
-  MPI_Init_thread(&argc, &argv, MPI_THREAD_FUNNELED, &provided);
-  if (provided < MPI_THREAD_FUNNELED) {
-    std::cerr << "Cannot mix MPI and Threads here. Please recompile with "
-                 "OpenMP or MPI switched off."
-              << std::endl;
-    MPI_Abort(MPI_COMM_WORLD, 99);
-  }
-  char *fname1, *fname2;
+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},
@ -98,17 +93,17 @@ int main(int argc, char** argv) {
                        {"CX", &cx, MINIARG_DOUBLE},
                        {"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)) return 1;
+  if (!parseMiniArgs(argc, argv, args))
+    return 1;

-  int rank, size;
-  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
-  MPI_Comm_size(MPI_COMM_WORLD, &size);
  H5::H5File in_f(fname1, 0);
-  H5::H5File out_f(boost::str(boost::format("fields_%d.h5", rank), H5F_ACC_TRUNC);
+  H5::H5File out_f(outFile, H5F_ACC_TRUNC);
  array_type v1_data;
-  uint32_t N1_points, N2_points;
+  uint64_t N1_points, N2_points;

  array3_type bins(boost::extents[Nres][Nres][Nres]);

@ -121,15 +116,17 @@ int main(int argc, char** argv) {

  cout << "Got " << N1_points << " in the first file." << endl;

-  MyCell* allCells_1 = new MyCell[N1_points];
+  MyCell *allCells_1 = new MyCell[N1_points];

 #pragma omp parallel for schedule(static)
-  for (long i = 0; i < Nres * Nres * Nres; i++) bins.data()[i] = 0;
+  for (uint32_t i = 0; i < Nres * Nres * Nres; i++)
+    bins.data()[i] = 0;

  cout << "Shuffling data in cells..." << endl;
 #pragma omp parallel for schedule(static)
-  for (int i = 0; i < N1_points; i++) {
-    for (int j = 0; j < 3; j++) allCells_1[i].coord[j] = v1_data[i][j];
+  for (uint64_t i = 0; i < N1_points; i++) {
+    for (int j = 0; j < 3; j++)
+      allCells_1[i].coord[j] = v1_data[i][j];
    for (int k = 0; k < 3; k++)
      allCells_1[i].val.pValue.v[k] = v1_data[i][3 + k];
    allCells_1[i].val.pValue.mass = v1_data[i][6];
@ -143,8 +140,9 @@ int main(int argc, char** argv) {
    if (rx < 0 || rx >= Nres || ry < 0 || ry >= Nres || rz < 0 || rz >= Nres)
      continue;

-    //#pragma omp atomic update
-    bins[rx][ry][rz]++;
+    auto &b = bins[rx][ry][rz];
+#pragma omp atomic
+    b++;
  }
  v1_data.resize(boost::extents[1][1]);

@ -152,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;

@ -159,18 +158,21 @@ int main(int argc, char** argv) {

  cout << "Weighing..." << endl;

-#pragma omp parallel
+  int rz_max = 0;
+#pragma omp parallel shared(rz_max)
  {
    MySmooth smooth1(&tree1, N_SPH);

-#pragma omp for schedule(dynamic)
+#pragma omp for collapse(3) schedule(dynamic, 8)
    for (int rz = 0; rz < Nres; rz++) {
-      double pz = (rz)*boxsize / Nres - cz;
-
-      (cout << rz << " / " << Nres << endl).flush();
      for (int ry = 0; ry < Nres; ry++) {
-        double py = (ry)*boxsize / Nres - cy;
        for (int rx = 0; rx < Nres; rx++) {
+          if (rz > rz_max) {
+            rz_max = rz;
+            (cout << rz << " / " << Nres << endl).flush();
+          }
+          double pz = (rz)*boxsize / Nres - cz;
+          double py = (ry)*boxsize / Nres - cy;
          double px = (rx)*boxsize / Nres - cx;

          MyTree::coords c = {float(px), float(py), float(pz)};
@ -185,30 +187,11 @@ int main(int argc, char** argv) {
            smooth1.fetchNeighbours(c, numInCell);
          else
            smooth1.fetchNeighbours(c);
-#pragma omp critical
          smooth1.addGridSite(c);
        }
      }
-      (cout << " Done " << rz << endl).flush();
    }
  }
-  //
-  // Reduction on the cell.weight in the tree.
-  // MPI_Allreduce to act on contiguous arrays.
-  auto tree = smooth1.getTree();
-  auto nodes = tree->getAllNodes();
-  double *weight_array = new double[N1_points];
-  for (size_t c = 0; c < N1_points; c++) {
-    weight_array[c] = allCells[c].val.weight;
-  }
-  MPI_Allreduce(MPI_IN_PLACE, weight_array, N1_points, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
-  for (size_t c = 0; c < N1_points; c++) {
-    allCells[c].val.weight = weight_array[c];
-  }
-  delete[] weight_array;
-
-  // cell.weights  -> build a 1d array of the particles weight -> MPI_Allreduce -> resend the new weights to the particles
-  //

  cout << "Interpolating..." << endl;

--- a/setup.py
+++ b/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
@ -76,9 +77,9 @@ class InstallCMakeLibs(install_lib):
 #        # your files are moved to the appropriate location when the installation
 #        # is run
 #
-#        libs = [os.path.join(bin_dir, _lib) for _lib in 
-#                os.listdir(bin_dir) if 
-#                os.path.isfile(os.path.join(bin_dir, _lib)) and 
+#        libs = [os.path.join(bin_dir, _lib) for _lib in
+#                os.listdir(bin_dir) if
+#                os.path.isfile(os.path.join(bin_dir, _lib)) and
 #                os.path.splitext(_lib)[1] in [".dll", ".so"]
 #                and not (_lib.startswith("python") or _lib.startswith(PACKAGE_NAME))]
 #
@ -87,16 +88,16 @@ class InstallCMakeLibs(install_lib):
 #            shutil.move(lib, os.path.join(self.build_dir,
 #                                          os.path.basename(lib)))
 #
-#        # Mark the libs for installation, adding them to 
-#        # distribution.data_files seems to ensure that setuptools' record 
+#        # Mark the libs for installation, adding them to
+#        # distribution.data_files seems to ensure that setuptools' record
 #        # writer appends them to installed-files.txt in the package's egg-info
 #        #
-#        # Also tried adding the libraries to the distribution.libraries list, 
-#        # but that never seemed to add them to the installed-files.txt in the 
-#        # egg-info, and the online recommendation seems to be adding libraries 
-#        # into eager_resources in the call to setup(), which I think puts them 
-#        # in data_files anyways. 
-#        # 
+#        # Also tried adding the libraries to the distribution.libraries list,
+#        # but that never seemed to add them to the installed-files.txt in the
+#        # egg-info, and the online recommendation seems to be adding libraries
+#        # into eager_resources in the call to setup(), which I think puts them
+#        # in data_files anyways.
+#        #
 #        # What is the best way?
 #
 #        # These are the additional installation files that should be
@ -104,7 +105,7 @@ class InstallCMakeLibs(install_lib):
 #        # step; depending on the files that are generated from your cmake
 #        # build chain, you may need to modify the below code
 #
-#        self.distribution.data_files = [os.path.join(self.install_dir, 
+#        self.distribution.data_files = [os.path.join(self.install_dir,
 #                                                     os.path.basename(lib))
 #                                        for lib in libs]
 #        print(self.distribution.data_files)
@ -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',
@ -213,7 +221,7 @@ class BuildCMakeExt(build_ext):
               if  _pyd_top[0].startswith(PACKAGE_NAME):
                   if os.path.splitext(_pyd)[1] in [".pyd", ".so"] or _pyd_top[-1] == 'config.py':
                     pyd_path.append((_pyd_top,_pyd))
-               
+

        for top,p in pyd_path:
           _,n = os.path.split(p)
@ -229,10 +237,10 @@ class BuildCMakeExt(build_ext):
 CosmoTool_extension = CMakeExtension(name="cosmotool")

 setup(name='cosmotool',
-      version='1.2.3',
+      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',
--- a/src/cosmopower.cpp
+++ b/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
@ -249,6 +273,19 @@ double CosmoPower::matterTransferFunctionHu(double k)
  return T_k;
 }

+double CosmoPower::noWiggleTk(double k)
+{
+  double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75)));
+  double alpha_Gamma = 1 - 0.328 * log(431 * OmegaEff) * OMEGA_B / OMEGA_0 + 0.38 * log(22.3 * OmegaEff) * pow(OMEGA_B / OMEGA_0, 2);
+  double GammaEff = OMEGA_0 * h * (alpha_Gamma + (1 - alpha_Gamma)/(1 + pow(0.43 * k * s, 4)));
+  double q = k/(h*GammaEff) * pow(Theta_27, 2);
+  double L_0 = log(2 * M_E + 1.8 * q);
+  double C_0 = 14.2 + 731 / (1 + 62.5 * q);
+  double T0 = L_0 / (L_0 + C_0 * q * q);
+
+  return T0;
+}
+
 double CosmoPower::powerHuBaryons(double k)
 {
  double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75)));
@ -444,12 +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 POWER_TEST:
      eval = &CosmoPower::powerTest;
      break;
+    case BAO_TK:
+      eval = &CosmoPower::BAO_Tk;
+      break;
    default:
      abort();
    }
--- a/src/cosmopower.hpp
+++ b/src/cosmopower.hpp
@ -89,7 +89,9 @@ namespace CosmoTool {
      POWER_SUGIYAMA,
      POWER_BDM,
      POWER_TEST,
-      HU_WIGGLES_ORIGINAL
+      HU_WIGGLES_ORIGINAL,
+      SAMPLE_WIGGLES,
+      BAO_TK
    };

    CosmoPower();
@ -122,6 +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);
  };

 };
--- a/src/h5_readFlash.cpp
+++ b/src/h5_readFlash.cpp
@ -7,16 +7,16 @@ This software is a computer program whose purpose is to provide a toolbox for co
 data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)

 This software is governed by the CeCILL license under French law and
-abiding by the rules of distribution of free software.  You can  use, 
+abiding by the rules of distribution of free software.  You can  use,
 modify and/ or redistribute the software under the terms of the CeCILL
 license as circulated by CEA, CNRS and INRIA at the following URL
-"http://www.cecill.info". 
+"http://www.cecill.info".

 As a counterpart to the access to the source code and  rights to copy,
 modify and redistribute granted by the license, users are provided only
 with a limited warranty  and the software's author,  the holder of the
 economic rights,  and the successive licensors  have only  limited
-liability. 
+liability.

 In this respect, the user's attention is drawn to the risks associated
 with loading,  using,  modifying and/or developing or reproducing the
@ -25,9 +25,9 @@ that may mean  that it is complicated to manipulate,  and  that  also
 therefore means  that it is reserved for developers  and  experienced
 professionals having in-depth computer knowledge. Users are therefore
 encouraged to load and test the software's suitability as regards their
-requirements in conditions enabling the security of their systems and/or 
-data to be ensured and,  more generally, to use and operate it in the 
-same conditions as regards security. 
+requirements in conditions enabling the security of their systems and/or
+data to be ensured and,  more generally, to use and operate it in the
+same conditions as regards security.

 The fact that you are presently reading this means that you have had
 knowledge of the CeCILL license and that you accept its terms.
@ -59,8 +59,8 @@ int ipvz_out   = 5;

 /* xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx */

-/*  n*_runtime_parameters should be set by the caller to 
-    the maximum number of runtime parameters to read. 
+/*  n*_runtime_parameters should be set by the caller to
+    the maximum number of runtime parameters to read.
 */

 void h5_read_runtime_parameters
@ -100,12 +100,12 @@ void h5_read_runtime_parameters
  nreal_runtime_parameters = MAX_PARM;

  /* integer runtime parameters */
-  int_rt_parms = (int_runtime_params_t *) malloc(nint_runtime_parameters * sizeof(int_runtime_params_t)); 
+  int_rt_parms = (int_runtime_params_t *) malloc(nint_runtime_parameters * sizeof(int_runtime_params_t));

  rank = 1;
  DataSet dataset = file->openDataSet("integer runtime parameters");
-  
-  IntType int_rt_type = dataset.getIntType();  
+
+  IntType int_rt_type = dataset.getIntType();
  //int_rt_type = H5Dget_type(dataset);

  DataSpace dataspace = dataset.getSpace();
@ -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++) {
@ -145,14 +142,14 @@ void h5_read_runtime_parameters

  /* reals */

-  real_rt_parms = (real_runtime_params_t *) malloc(nreal_runtime_parameters * sizeof(real_runtime_params_t)); 
-  
+  real_rt_parms = (real_runtime_params_t *) malloc(nreal_runtime_parameters * sizeof(real_runtime_params_t));
+
  rank = 1;
  dataset = file->openDataSet("real runtime parameters");
-  //dataset = H5Dopen(*file_identifier, "real runtime parameters"); 
-  
+  //dataset = H5Dopen(*file_identifier, "real runtime parameters");
+
  dataspace = dataset.getSpace();
-  FloatType real_rt_type = dataset.getFloatType();  
+  FloatType real_rt_type = dataset.getFloatType();
  ndims = dataspace.getSimpleExtentDims(&dimens_1d, NULL);
  //dataspace = H5Dget_space(dataset);
  //real_rt_type = H5Dget_type(dataset);
@ -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++) {
@ -193,7 +187,7 @@ void h5_read_runtime_parameters
      *LBox = real_runtime_parameter_values[i];
    }
    if (strncmp(real_runtime_parameter_names[i],"hubbleconstant", 14) == 0 ) {
-      *hubble = real_runtime_parameter_values[i];      
+      *hubble = real_runtime_parameter_values[i];
    }
    if (strncmp(real_runtime_parameter_names[i],"omegamatter", 11) == 0 ) {
      *omegam = real_runtime_parameter_values[i];
@ -205,7 +199,7 @@ void h5_read_runtime_parameters
      *omegalambda = real_runtime_parameter_values[i];
    }
  }
-  
+
  for (i = 0; i < nint_runtime_parameters; i++) {
    if (strncmp(int_runtime_parameter_names[i],"pt_numx",7) == 0 ) {
      *numPart = int_runtime_parameter_values[i];
@ -214,7 +208,7 @@ void h5_read_runtime_parameters
  }
 }

-     
+
 /* xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
 void h5_read_flash3_particles (H5File* file,
                        int* totalparticles,
@ -241,8 +235,8 @@ void h5_read_flash3_particles (H5File* file,
  char     *propName;
  double   *partBuffer;
  char    part_names[50][OUTPUT_PROP_LENGTH];
-  int     string_size;  
-         
+  int     string_size;
+
 //  char part_names[NPART_PROPS][OUTPUT_PROP_LENGTH];
  hsize_t dimens_2d[2], maxdimens_2d[2];
  hsize_t start_2d[2], count_2d[2], stride_2d[2];
@ -265,12 +259,12 @@ void h5_read_flash3_particles (H5File* file,

  //total number of particle properties
  numProps = dimens_2d[0];
-  
+
  string_size = OUTPUT_PROP_LENGTH;
  StrType string_type = H5Tcopy(H5T_C_S1);
  string_type.setSize(string_size);
  //status = H5Tset_size(string_type, string_size);
-  
+
  rank = 2;

  start_2d[0] = 0;
@ -282,22 +276,20 @@ void h5_read_flash3_particles (H5File* file,
  count_2d[0] = dimens_2d[0];
  count_2d[1] = dimens_2d[1];

-  dataspace.selectHyperslab(H5S_SELECT_SET, count_2d, start_2d); 
+  dataspace.selectHyperslab(H5S_SELECT_SET, count_2d, start_2d);
  //status = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, start_2d,
  //                      stride_2d, count_2d, NULL);

  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();
  memspace.close();
  dataspace.close();
-  dataset.close(); 
+  dataset.close();
  //H5Tclose(string_type);
  //H5Sclose(memspace);
  //H5Sclose(dataspace);
@ -313,7 +305,7 @@ void h5_read_flash3_particles (H5File* file,
    if (strncmp(part_names[i], "velz", 4) == 0) { ipvz = i+1; }
  }

-  if ((iptag < 0) || (ipx < 0) || (ipy < 0) || (ipz < 0) || (ipvx < 0) || 
+  if ((iptag < 0) || (ipx < 0) || (ipy < 0) || (ipz < 0) || (ipvx < 0) ||
       (ipvy < 0) || (ipvz < 0) ) {
    printf("One or more required particle attributes not found in file!\n");
    return;
@ -325,7 +317,7 @@ void h5_read_flash3_particles (H5File* file,
  //read particles
  dataset = file->openDataSet("tracer particles");
  //dataset = H5Dopen(*file_identifier, "tracer particles");
-  
+
  FloatType datatype = dataset.getFloatType();
  //datatype = H5Dget_type(dataset);

@ -338,7 +330,7 @@ void h5_read_flash3_particles (H5File* file,
  ndims = dataspace.getSimpleExtentDims(dimens_2d, NULL);
  //dataspace = H5Dget_space(dataset);
  //H5Sget_simple_extent_dims(dataspace, dimens_2d, maxdimens_2d);
- 
+
  /*insert particle properties (numPartBuffer) particles at a time*/
  pstack = (*localnp);
  poffset = 0;
@ -366,7 +358,7 @@ void h5_read_flash3_particles (H5File* file,
    dimens_2d[0] = (pcount);
    dimens_2d[1] = (numProps);

-    dataspace.selectHyperslab(H5S_SELECT_SET, count_2d, start_2d); 
+    dataspace.selectHyperslab(H5S_SELECT_SET, count_2d, start_2d);
    //status     = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, start_2d,
    //                                 stride_2d, count_2d, NULL);

@ -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 */

@ -391,8 +382,8 @@ void h5_read_flash3_particles (H5File* file,
      pos3[p+poffset]  = (float) *(partBuffer+ipz-1+p*numProps);
    }
    }
-   
-    
+
+
    if (vel1 && vel2 && vel3) {
    for(p=0; p < (pcount); p++) {
      vel1[p+poffset]  = (float) *(partBuffer+ipvx-1+p*numProps);
@ -423,7 +414,7 @@ void h5_read_flash3_particles (H5File* file,
  //status = H5Sclose(dataspace);
  //status = H5Dclose(dataset);
 }
- 
+

 /*xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/

@ -451,7 +442,7 @@ void h5_read_flash3_header_info(H5File* file,

  H5std_string DATASET_NAME;

-  string_type = H5Tcopy(H5T_C_S1);  
+  string_type = H5Tcopy(H5T_C_S1);
  H5Tset_size(string_type, MAX_STRING_LENGTH);

  DataSet dataset = file->openDataSet("real scalars");
@ -467,13 +458,13 @@ void h5_read_flash3_header_info(H5File* file,
  /* malloc a pointer to a list of real_list_t's */
  real_list = (real_list_t *) malloc(dimens_1d * sizeof(real_list_t));

-  // create a new simple dataspace of 1 dimension and size of 'dimens_1d' 
+  // create a new simple dataspace of 1 dimension and size of 'dimens_1d'
  DataSpace memspace(1, &dimens_1d);

-  // create an empty vessel sized to hold one real_list_t's worth of data 
+  // create an empty vessel sized to hold one real_list_t's worth of data
  CompType real_list_type( sizeof(real_list_t) );

-  // subdivide the empty vessel into its component sections (name and value) 
+  // subdivide the empty vessel into its component sections (name and value)
  real_list_type.insertMember(
          "name",
          HOFFSET(real_list_t, name),
@ -484,9 +475,8 @@ void h5_read_flash3_header_info(H5File* file,
          HOFFSET(real_list_t, value),
          PredType::NATIVE_DOUBLE);

-  // read the data into 'real_list' 
-  dataset.read( real_list, real_list_type, memspace, dataspace,
-                H5P_DEFAULT);
+  // read the data into 'real_list'
+  dataset.read( real_list, real_list_type, memspace, dataspace);


  if (status < 0) {
--- a/src/hdf5_array.hpp
+++ b/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()); \
     } 
--- a/src/interpolate.cpp
+++ b/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
--- a/src/mykdtree.hpp
+++ b/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

@ -7,16 +7,16 @@ This software is a computer program whose purpose is to provide a toolbox for co
 data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)

 This software is governed by the CeCILL license under French law and
-abiding by the rules of distribution of free software.  You can  use, 
+abiding by the rules of distribution of free software.  You can  use,
 modify and/ or redistribute the software under the terms of the CeCILL
 license as circulated by CEA, CNRS and INRIA at the following URL
-"http://www.cecill.info". 
+"http://www.cecill.info".

 As a counterpart to the access to the source code and  rights to copy,
 modify and redistribute granted by the license, users are provided only
 with a limited warranty  and the software's author,  the holder of the
 economic rights,  and the successive licensors  have only  limited
-liability. 
+liability.

 In this respect, the user's attention is drawn to the risks associated
 with loading,  using,  modifying and/or developing or reproducing the
@ -25,9 +25,9 @@ that may mean  that it is complicated to manipulate,  and  that  also
 therefore means  that it is reserved for developers  and  experienced
 professionals having in-depth computer knowledge. Users are therefore
 encouraged to load and test the software's suitability as regards their
-requirements in conditions enabling the security of their systems and/or 
-data to be ensured and,  more generally, to use and operate it in the 
-same conditions as regards security. 
+requirements in conditions enabling the security of their systems and/or
+data to be ensured and,  more generally, to use and operate it in the
+same conditions as regards security.

 The fact that you are presently reading this means that you have had
 knowledge of the CeCILL license and that you accept its terms.
@ -48,13 +48,13 @@ namespace CosmoTool {

  typedef uint64_t NodeIntType;

-  template<int N, typename CType = ComputePrecision> 
+  template<int N, typename CType = ComputePrecision>
  struct KDDef
  {
    typedef CType CoordType;
    typedef float KDCoordinates[N];
  };
-  
+
  template<int N, typename ValType, typename CType = ComputePrecision>
  struct KDCell
  {
@ -99,10 +99,10 @@ namespace CosmoTool {
    typename KDDef<N,CType>::CoordType r, r2;
    KDCell<N, ValType,CType> **cells;
    typename KDDef<N,CType>::CoordType *distances;
-    uint32_t currentRank;
-    uint32_t numCells;
+    uint64_t currentRank;
+    uint64_t numCells;
  };
-  
+

  template<int N, typename ValType, typename CType = ComputePrecision>
  class RecursionMultipleInfo
@ -114,14 +114,14 @@ namespace CosmoTool {

    RecursionMultipleInfo(const typename KDDef<N,CType>::KDCoordinates& rx,
 			  KDCell<N,ValType,CType> **cells,
-			  uint32_t numCells)
+			  uint64_t numCells)
      : queue(cells, numCells, INFINITY),traversed(0)
    {
      std::copy(rx, rx+N, x);
    }
  };

-  template<int N, typename ValType, typename CType = ComputePrecision> 
+  template<int N, typename ValType, typename CType = ComputePrecision>
  struct KD_default_cell_splitter
  {
    void operator()(KDCell<N,ValType,CType> **cells, NodeIntType Ncells, NodeIntType& split_index, int axis, typename KDDef<N,CType>::KDCoordinates minBound, typename KDDef<N,CType>::KDCoordinates maxBound);
@ -135,7 +135,7 @@ namespace CosmoTool {
    typedef typename KDDef<N>::KDCoordinates coords;
    typedef KDCell<N,ValType,CType> Cell;
    typedef KDTreeNode<N,ValType,CType> Node;
-    
+
    CellSplitter splitter;

    KDTree(Cell *cells, NodeIntType Ncells);
@ -153,20 +153,20 @@ namespace CosmoTool {
       std::copy(replicate, replicate+N, this->replicate);
     }

-    uint32_t getIntersection(const coords& x, CoordType r, 
+    uint64_t getIntersection(const coords& x, CoordType r,
 			     Cell **cells,
-			     uint32_t numCells);
-    uint32_t getIntersection(const coords& x, CoordType r, 
+			     uint64_t numCells);
+    uint64_t getIntersection(const coords& x, CoordType r,
 			     Cell **cells,
 			     CoordType *distances,
-			     uint32_t numCells);
-    uint32_t countCells(const coords& x, CoordType r);
+			     uint64_t numCells);
+    uint64_t countCells(const coords& x, CoordType r);

    Cell *getNearestNeighbour(const coords& x);

-    void getNearestNeighbours(const coords& x, uint32_t NumCells,
+    void getNearestNeighbours(const coords& x, uint64_t NumCells,
 			      Cell **cells);
-    void getNearestNeighbours(const coords& x, uint32_t NumCells,
+    void getNearestNeighbours(const coords& x, uint64_t NumCells,
 			      Cell **cells,
 			      CoordType *distances);

@ -183,7 +183,7 @@ namespace CosmoTool {
    NodeIntType getNumberInNode(const Node *n) const { return n->numNodes; }
 #else
    NodeIntType getNumberInNode(const Node *n) const {
-      if (n == 0) 
+      if (n == 0)
        return 0;
      return 1+getNumberInNode(n->children[0])+getNumberInNode(n->children[1]);
    }
@ -211,7 +211,7 @@ namespace CosmoTool {
 		    uint32_t depth,
 		    coords minBound,
 		    coords maxBound);
-    
+
    template<bool justCount>
    void recursiveIntersectionCells(RecursionInfoCells<N,ValType, CType>& info,
 				    Node *node,
@ -224,7 +224,7 @@ namespace CosmoTool {
 			  CoordType& R2,
 			  Cell*& cell);
    void recursiveMultipleNearest(RecursionMultipleInfo<N,ValType,CType>& info, Node *node,
-				  int level);    
+				  int level);

  };

--- a/src/mykdtree.tcc
+++ b/src/mykdtree.tcc
@ -33,7 +33,7 @@ namespace CosmoTool {
  template<int N, typename ValType, typename CType, typename CellSplitter>
  KDTree<N,ValType,CType,CellSplitter>::KDTree(Cell *cells, NodeIntType Ncells)
  {
-    periodic = false; 
+    periodic = false;

    base_cell = cells;
    numNodes = Ncells;
@ -41,7 +41,7 @@ namespace CosmoTool {

    sortingHelper = new Cell *[Ncells];
    for (NodeIntType i = 0; i < Ncells; i++)
-	sortingHelper[i] = &cells[i];    
+	sortingHelper[i] = &cells[i];

    optimize();
  }
@ -73,16 +73,16 @@ namespace CosmoTool {
            absoluteMax[k] = cell->coord[k];
        }
      }
-    
+
    std::cout << "   rebuilding the tree..." << std::endl;
-    root = buildTree(sortingHelper, activeCells, 0, absoluteMin, absoluteMax);    
+    root = buildTree(sortingHelper, activeCells, 0, absoluteMin, absoluteMax);
    std::cout << "   done." << std::endl;
  }

  template<int N, typename ValType, typename CType, typename CellSplitter>
-  uint32_t KDTree<N,ValType,CType,CellSplitter>::getIntersection(const coords& x, CoordType r, 
+  uint64_t KDTree<N,ValType,CType,CellSplitter>::getIntersection(const coords& x, CoordType r,
 						    KDTree<N,ValType,CType,CellSplitter>::Cell **cells,
-						    uint32_t numCells)
+						    uint64_t numCells)
  {
    RecursionInfoCells<N,ValType,CType> info;

@ -112,10 +112,10 @@ namespace CosmoTool {
  }

  template<int N, typename ValType, typename CType, typename CellSplitter>
-  uint32_t KDTree<N,ValType,CType,CellSplitter>::getIntersection(const coords& x, CoordType r, 
+  uint64_t KDTree<N,ValType,CType,CellSplitter>::getIntersection(const coords& x, CoordType r,
 						    Cell **cells,
 						    CoordType *distances,
-						    uint32_t numCells)
+						    uint64_t numCells)
  {
    RecursionInfoCells<N,ValType> info;

@ -144,7 +144,7 @@ namespace CosmoTool {
  }

  template<int N, typename ValType, typename CType, typename CellSplitter>
-  uint32_t KDTree<N,ValType,CType,CellSplitter>::countCells(const coords& x, CoordType r)
+  uint64_t KDTree<N,ValType,CType,CellSplitter>::countCells(const coords& x, CoordType r)
  {
    RecursionInfoCells<N,ValType> info;

@ -175,7 +175,7 @@ namespace CosmoTool {

  template<int N, typename ValType, typename CType, typename CellSplitter>
  template<bool justCount>
-  void KDTree<N,ValType,CType,CellSplitter>::recursiveIntersectionCells(RecursionInfoCells<N,ValType,CType>& info, 
+  void KDTree<N,ValType,CType,CellSplitter>::recursiveIntersectionCells(RecursionInfoCells<N,ValType,CType>& info,
 							   Node *node,
 							   int level)
  {
@ -183,7 +183,7 @@ namespace CosmoTool {
    CoordType d2 = 0;

 #if __KD_TREE_ACTIVE_CELLS == 1
-    if (node->value->active) 
+    if (node->value->active)
 #endif
      {
 	for (int j = 0; j < 3; j++)
@ -250,9 +250,9 @@ namespace CosmoTool {
  }

  template<int N, typename ValType, typename CType>
-  void KD_default_cell_splitter<N,ValType,CType>::operator()(KDCell<N,ValType,CType> **cells, NodeIntType Ncells, 
+  void KD_default_cell_splitter<N,ValType,CType>::operator()(KDCell<N,ValType,CType> **cells, NodeIntType Ncells,
                                                             NodeIntType& split_index, int axis,
-                                                             typename KDDef<N,CType>::KDCoordinates minBound, 
+                                                             typename KDDef<N,CType>::KDCoordinates minBound,
                                                             typename KDDef<N,CType>::KDCoordinates maxBound)
  {
    CellCompare<N,ValType,CType> compare(axis);
@ -279,9 +279,9 @@ namespace CosmoTool {

 //#pragma omp atomic capture
    nodeId = (this->lastNode)++;
-    
+
    node = &nodes[nodeId];
-    
+
    // Isolate the environment
    splitter(cell0, Ncells, mid, axis, minBound, maxBound);

@ -297,12 +297,12 @@ namespace CosmoTool {
    {
      node->children[0] = buildTree(cell0, mid, depth, minBound, tmpBound);
    }
-    
+
    memcpy(tmpBound, minBound, sizeof(coords));
    tmpBound[axis] = node->value->coord[axis];
 #pragma omp task private(tmpBound)
    {
-      node->children[1] = buildTree(cell0+mid+1, Ncells-mid-1, depth, 
+      node->children[1] = buildTree(cell0+mid+1, Ncells-mid-1, depth,
 				  tmpBound, maxBound);
    }

@ -391,17 +391,17 @@ namespace CosmoTool {
 	  }
 	return;
      }
-    
+
    // Check if current node is not the nearest
-    CoordType thisR2 = 
+    CoordType thisR2 =
      computeDistance(node->value, x);
-    
+
    if (thisR2 < R2)
      {
 	R2 = thisR2;
 	best = node->value;
-      }    
-    
+      }
+
    // Now we found the best. We check whether the hypersphere
    // intersect the hyperplane of the other branch

@ -435,11 +435,11 @@ namespace CosmoTool {
         {
           coords x_new;
           r.getPosition(x_new);
-           recursiveNearest(root, 0, x_new, R2, best); 
+           recursiveNearest(root, 0, x_new, R2, best);
         }
       while (r.next());
     }
-   
+
   return best;
 }

@ -474,15 +474,15 @@ namespace CosmoTool {
      {
 	recursiveMultipleNearest(info, go, level+1);
      }
-    
+
    // Check if current node is not the nearest
-    CoordType thisR2 = 
+    CoordType thisR2 =
      computeDistance(node->value, info.x);
    info.queue.push(node->value, thisR2);
    info.traversed++;
    //    if (go == 0)
    //      return;
-    
+
    // Now we found the best. We check whether the hypersphere
    // intersect the hyperplane of the other branch

@ -497,15 +497,15 @@ namespace CosmoTool {
 	  {
 	    recursiveMultipleNearest(info, other, level+1);
 	  }
-      }    
+      }
  }

 template<int N, typename ValType, typename CType, typename CellSplitter>
- void KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbours(const coords& x, uint32_t N2,
+ void KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbours(const coords& x, uint64_t N2,
 						    Cell **cells)
 {
   RecursionMultipleInfo<N,ValType> info(x, cells, N2);
-   
+
   for (int i = 0; i < N2; i++)
     cells[i] = 0;

@ -527,12 +527,12 @@ namespace CosmoTool {
 }

 template<int N, typename ValType, typename CType, typename CellSplitter>
- void KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbours(const coords& x, uint32_t N2,
+ void KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbours(const coords& x, uint64_t N2,
 						    Cell **cells,
 						    CoordType *distances)
 {
   RecursionMultipleInfo<N,ValType> info(x, cells, N2);
-   
+
   for (int i = 0; i < N2; i++)
     cells[i] = 0;

@ -555,14 +555,14 @@ namespace CosmoTool {
 #ifdef __KD_TREE_SAVE_ON_DISK
 #define KDTREE_DISK_SIGNATURE "KDTREE"
 #define KDTREE_DISK_SIGNATURE_LEN 7
-  
+
  template<int N, typename CType>
  struct KDTreeOnDisk
  {
    long cell_id;
    long children_node[2];
    typename KDDef<N, CType>::KDCoordinates minBound, maxBound;
-  };  
+  };

  struct KDTreeHeader
  {
@ -619,7 +619,7 @@ namespace CosmoTool {
      {
         std::cerr << "KDTree Signature invalid" << std::endl;
         throw InvalidOnDiskKDTree();
-      } 
+      }

    if (h.numCells != Ncells || h.nodesUsed < 0) {
      std::cerr << "The number of cells has changed (" << h.numCells << " != " << Ncells << ") or nodesUsed=" << h.nodesUsed << std::endl;
@ -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;
@ -683,10 +683,10 @@ namespace CosmoTool {
      }

    root = &nodes[h.rootId];
-    
+
    sortingHelper = new Cell *[Ncells];
    for (NodeIntType i = 0; i < Ncells; i++)
-      sortingHelper[i] = &cells[i];    
+      sortingHelper[i] = &cells[i];
  }
 #endif

--- a/src/powerSpectrum.cpp
+++ b/src/powerSpectrum.cpp
@ -58,6 +58,7 @@ using namespace std;
 #define POWER_TEST 8

 #define POWER_SPECTRUM POWER_EFSTATHIOU
+#define SAMPLE_WIGGLES 9

 namespace Cosmology {

@ -206,6 +207,31 @@ double powG(double y)
 */
 double powerSpectrum(double k, double normPower)
 {
+#if POWER_SPECTRUM == SAMPLE_WIGGLES
+  // BAO wiggle parameterization for reconstruction
+  // Babic et al. 2022, https://arxiv.org/abs/2203.06177
+
+  // No-wiggle transfer function
+
+  double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75)));
+  double alpha_Gamma = 1 - 0.328 * log(431 * OmegaEff) * OMEGA_B / OMEGA_0 + 0.38 * log(22.3 * OmegaEff) * pow(OMEGA_B / OMEGA_0, 2);
+  double GammaEff = OMEGA_0 * h * (alpha_Gamma + (1 - alpha_Gamma)/(1 + pow(0.43 * k * s, 4)));
+  double q = k/(h*GammaEff) * pow(Theta_27, 2);
+  double L_0 = log(2 * M_E + 1.8 * q);
+  double C_0 = 14.2 + 731 / (1 + 62.5 * q);
+  double T0 = L_0 / (L_0 + C_0 * q * q);
+
+  // Wiggle parameterization
+
+  double A = 0;
+  double r_s = 10;
+  double k_D = 2 * M_PI / 100;
+
+  double param = 1 + A * sin(k * r_s) * exp(- k / k_D);
+
+  return normPower * pow(k, n) * T0 * T0 * param;
+#endif
+
 #if POWER_SPECTRUM == POWER_EFSTATHIOU
  double a = 6.4/Gamma0;
  double b = 3/Gamma0;
--- a/src/sphSmooth.hpp
+++ b/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

@ -7,16 +7,16 @@ This software is a computer program whose purpose is to provide a toolbox for co
 data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)

 This software is governed by the CeCILL license under French law and
-abiding by the rules of distribution of free software.  You can  use, 
+abiding by the rules of distribution of free software.  You can  use,
 modify and/ or redistribute the software under the terms of the CeCILL
 license as circulated by CEA, CNRS and INRIA at the following URL
-"http://www.cecill.info". 
+"http://www.cecill.info".

 As a counterpart to the access to the source code and  rights to copy,
 modify and redistribute granted by the license, users are provided only
 with a limited warranty  and the software's author,  the holder of the
 economic rights,  and the successive licensors  have only  limited
-liability. 
+liability.

 In this respect, the user's attention is drawn to the risks associated
 with loading,  using,  modifying and/or developing or reproducing the
@ -25,9 +25,9 @@ that may mean  that it is complicated to manipulate,  and  that  also
 therefore means  that it is reserved for developers  and  experienced
 professionals having in-depth computer knowledge. Users are therefore
 encouraged to load and test the software's suitability as regards their
-requirements in conditions enabling the security of their systems and/or 
-data to be ensured and,  more generally, to use and operate it in the 
-same conditions as regards security. 
+requirements in conditions enabling the security of their systems and/or
+data to be ensured and,  more generally, to use and operate it in the
+same conditions as regards security.

 The fact that you are presently reading this means that you have had
 knowledge of the CeCILL license and that you accept its terms.
@ -39,80 +39,89 @@ 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;

-    typedef KDTree<Ndims,FullType> SPHTree;
-    typedef KDTreeNode<Ndims,FullType> SPHNode;
-    typedef KDCell<Ndims,FullType>  SPHCell;
-    typedef typename KDTree<Ndims,FullType>::CoordType CoordType;
+    typedef KDTree<Ndims, FullType> SPHTree;
+    typedef KDTreeNode<Ndims, FullType> SPHNode;
+    typedef KDCell<Ndims, FullType> SPHCell;
+    typedef typename KDTree<Ndims, FullType>::CoordType CoordType;

-    typedef ComputePrecision (*computeParticleValue)(const ValType& t);
-    typedef void (*runParticleValue)(ValType& t);
+    typedef ComputePrecision (*computeParticleValue)(const ValType &t);
+    typedef void (*runParticleValue)(ValType &t);

  public:
    typedef SPHCell *P_SPHCell;
-    struct SPHState
-    {
+    struct SPHState {
      boost::shared_ptr<P_SPHCell[]> ngb;
      boost::shared_ptr<CoordType[]> distances;
      typename SPHTree::coords currentCenter;
      int currentNgb;
      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, uint32_t newNsph);
-    void fetchNeighboursOnVolume(const typename SPHTree::coords& c, ComputePrecision radius);
-    const typename SPHTree::coords& getCurrentCenter() const
-    {
+    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, uint32_t newNsph,
+        SPHState *state = 0);
+    void fetchNeighboursOnVolume(
+        const typename SPHTree::coords &c, ComputePrecision radius);
+    const typename SPHTree::coords &getCurrentCenter() const {
      return internal.currentCenter;
    }

-    template<typename FuncT>
-    ComputePrecision computeSmoothedValue(const typename SPHTree::coords& c,
-					  FuncT fun, SPHState *state = 0);    
-	
-	template<typename FuncT>
-    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 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 getSmoothingLen() const
-    {
-      return internal.smoothRadius;
-    }
+    /** 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);
+
+    /** 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
+    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>
+    template <typename FuncT>
    void runForEachNeighbour(FuncT fun, SPHState *state = 0);
-    void addGridSite(const typename SPHTree::coords& c);
+    void addGridSite(const typename SPHTree::coords &c, SPHState *state);
+    void addGridSite(const typename SPHTree::coords &c);

    bool hasNeighbours() const;

-    virtual ComputePrecision getKernel(ComputePrecision d) const;    
+    virtual ComputePrecision getKernel(ComputePrecision d) const;

    SPHTree *getTree() { return tree; }

@ -125,29 +134,29 @@ namespace CosmoTool
    uint32_t getCurrent() const { return internal.currentNgb; }

    uint32_t getNgb() const { return maxNgb; }
- 
+
  protected:
    SPHState internal;
    uint32_t Nsph;
    uint32_t deltaNsph;
    uint32_t maxNgb;
    SPHTree *tree;
-    
-    template<typename FuncT>
-    ComputePrecision computeWValue(const typename SPHTree::coords & c,
-				   SPHCell& cell,
-				   CoordType d,
-				   FuncT fun, SPHState *state);
-    
-    template<typename FuncT>
-    void runUnrollNode(SPHNode *node,
-		       FuncT fun);
+
+    template <typename FuncT>
+    ComputePrecision computeWValue(
+        const typename SPHTree::coords &c, SPHCell &cell, CoordType d,
+        FuncT fun, SPHState *state);
+
+    template <typename FuncT>
+    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);
+  template <class ValType1, class ValType2, int Ndims>
+  bool operator<(
+      const SPHSmooth<ValType1, Ndims> &s1,
+      const SPHSmooth<ValType2, Ndims> &s2);

-};
+}; // namespace CosmoTool

 #include "sphSmooth.tcc"

--- a/src/sphSmooth.tcc
+++ b/src/sphSmooth.tcc
@ -3,225 +3,248 @@

 namespace CosmoTool {

-template<typename ValType, int Ndims>
-SPHSmooth<ValType,Ndims>::SPHSmooth(SPHTree *tree, uint32_t Nsph)
-{
-  this->Nsph = Nsph;
-  this->tree = tree;
-  internal.currentNgb = 0;
+  template <typename ValType, int Ndims>
+  SPHSmooth<ValType, Ndims>::SPHSmooth(SPHTree *tree, uint32_t Nsph) {
+    this->Nsph = Nsph;
+    this->tree = tree;
+    internal.currentNgb = 0;

-  this->maxNgb = Nsph;
-  internal.ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[maxNgb]);  
-  internal.distances = boost::shared_ptr<CoordType[]>(new CoordType[maxNgb]);
-}
+    this->maxNgb = Nsph;
+    internal.ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[maxNgb]);
+    internal.distances = boost::shared_ptr<CoordType[]>(new CoordType[maxNgb]);
+  }

-template<typename ValType, int Ndims>
-SPHSmooth<ValType,Ndims>::~SPHSmooth()
-{
-}
+  template <typename ValType, int Ndims>
+  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)
-{
-  CoordType weight;
+  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) {
+    CoordType weight;

-  d /= state->smoothRadius;
-  weight = getKernel(d);
+    d /= state->smoothRadius;
+    weight = getKernel(d);

-  if (cell.val.weight != 0)
-    return weight * fun(cell.val.pValue) / cell.val.weight;
-  else
-    return 0;
-}
+    if (cell.val.weight != 0)
+      return weight * fun(cell.val.pValue) / cell.val.weight;
+    else
+      return 0;
+  }

-template<typename ValType, int Ndims>
-void
-SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c, uint32_t newNngb)
-{
-  ComputePrecision d2, max_dist = 0;
-  uint32_t requested = newNngb;
+  template <typename ValType, int Ndims>
+  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)
-    {
-      maxNgb = requested;
-      internal.ngb = boost::shared_ptr<P_SPHCell[]>(new P_SPHCell[maxNgb]);
-      internal.distances = boost::shared_ptr<CoordType[]>(new CoordType[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]);
+      }
    }

-  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)
-{
-  ComputePrecision d2, max_dist = 0;
-  uint32_t requested = Nsph;
+  template <typename ValType, int Ndims>
+  void SPHSmooth<ValType, Ndims>::fetchNeighbours(
+      const typename SPHTree::coords &c, SPHState *state) {
+    ComputePrecision d2, max_dist = 0;
+    uint32_t requested = Nsph;

-  if (state != 0) {
-    state->distances = boost::shared_ptr<CoordType[]>(new CoordType[Nsph]);
-    state->ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[Nsph]);
-  } else
-    state = &internal;
-    
-  memcpy(state->currentCenter, c, sizeof(c));
-  
-  tree->getNearestNeighbours(c, requested, state->ngb.get(), state->distances.get());
+    if (state != 0) {
+      state->distances = boost::shared_ptr<CoordType[]>(new CoordType[Nsph]);
+      state->ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[Nsph]);
+    } else
+      state = &internal;

-  state->currentNgb = 0;
-  for (uint32_t i = 0; i < requested && state->ngb[i] != 0; i++,state->currentNgb++)
-    {
+    memcpy(state->currentCenter, c, sizeof(c));
+
+    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 = state->distances[i] = sqrt(state->distances[i]);
+      if (d2 > max_dist)
+        max_dist = d2;
+    }
+
+    state->smoothRadius = max_dist / 2;
+  }
+
+  template <typename ValType, int Ndims>
+  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);
+
+    for (uint32_t i = 0; i < internal.currentNgb; i++) {
      d2 = internal.distances[i] = sqrt(internal.distances[i]);
      if (d2 > max_dist)
        max_dist = d2;
    }
+    internal.smoothRadius = max_dist / 2;
+  }

-  state->smoothRadius = max_dist / 2;  
-}
+  template <typename ValType, int Ndims>
+  template <typename FuncT>
+  ComputePrecision SPHSmooth<ValType, Ndims>::computeSmoothedValue(
+      const typename SPHTree::coords &c, FuncT fun, SPHState *state) {
+    if (state == 0)
+      state = &internal;

+    ComputePrecision outputValue = 0;
+    ComputePrecision max_dist = 0;
+    ComputePrecision r3 = cube(state->smoothRadius);

-template<typename ValType, int Ndims>
-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);
-
-  for (uint32_t i = 0; i < internal.currentNgb; i++)
-    {
-      d2 = internal.distances[i] = sqrt(internal.distances[i]);
-      if (d2 > max_dist)
-        max_dist = d2;
-    }
-  internal.smoothRadius = max_dist / 2;  
-}				
-
-template<typename ValType, int Ndims>
-template<typename FuncT>
-ComputePrecision 
-SPHSmooth<ValType,Ndims>::computeSmoothedValue(const typename SPHTree::coords& c,
-                                               FuncT fun, SPHState *state)
-{
-  if (state == 0)
-    state = &internal;
-    
-  ComputePrecision outputValue = 0;
-  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;
-}
+    return outputValue / r3;
+  }

-template<typename ValType>
-ComputePrecision interpolateOne(const ValType& t)
-{
-  return 1.0;
-}
+  template <typename ValType>
+  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)
-{
-  if (state == 0)
-    state = &internal;
-    
-  ComputePrecision outputValue = 0;
-  ComputePrecision max_dist = 0;
-  ComputePrecision weight = 0;
+  template <typename ValType, int Ndims>
+  template <typename FuncT>
+  void SPHSmooth<ValType, Ndims>::computeAdjointGradientSmoothedValue(
+      const typename SPHTree::coords &c, ComputePrecision ag_value, FuncT fun,
+      SPHState *state) {
+    if (state == 0)
+      state = &internal;

-  for (uint32_t i = 0; i < state->currentNgb; i++)
-    {
+    ComputePrecision outputValue = 0;
+    ComputePrecision max_dist = 0;
+    ComputePrecision weight = 0;
+
+    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);  
-}
+    return (outputValue == 0) ? 0 : (outputValue / weight);
+  }

-template<typename ValType, int Ndims>
-template<typename FuncT>
-void SPHSmooth<ValType,Ndims>::runForEachNeighbour(FuncT fun, SPHState *state)
-{
-  if (state == 0)
-    state = &internal;
-    
-  for (uint32_t i = 0; i < state->currentNgb; i++)
-    {
+  template <typename ValType, int Ndims>
+  template <typename FuncT>
+  void
+  SPHSmooth<ValType, Ndims>::runForEachNeighbour(FuncT fun, SPHState *state) {
+    if (state == 0)
+      state = &internal;
+
+    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, SPHState *state) {
+    ComputePrecision outputValue = 0;
+    ComputePrecision max_dist = 0;
+    ComputePrecision r3 = cube(state->smoothRadius);

-template<typename ValType, int Ndims>
-void SPHSmooth<ValType,Ndims>::addGridSite(const typename SPHTree::coords& c)
-{
-  ComputePrecision outputValue = 0;
-  ComputePrecision max_dist = 0;
-  
-  ComputePrecision r3 = cube(internal.smoothRadius);
-
-  for (uint32_t i = 0; i < internal.currentNgb; i++)
-    {
-      ComputePrecision d = internal.distances[i];
-      SPHCell& cell = *(internal.ngb[i]);
-      cell.val.weight += getKernel(d/internal.smoothRadius) / r3;
+    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>
-ComputePrecision 
-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)
-    {
+  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 {
+    // 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) {
      ComputePrecision d = 2 - x;
      return 0.25 * d * d * d;
-    }
-  else
-    return 0;
-}
+    } else
+      return 0;
+  }

-template<typename ValType, int Ndims>
-bool SPHSmooth<ValType,Ndims>::hasNeighbours() const
-{
-  return (internal.currentNgb != 0);
-}
+  template <typename ValType, int Ndims>
+  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)
-{
-  return (s1.getSmoothingLen() < s2.getSmoothingLen());
-}
-
-};
+  template <class ValType1, class ValType2, int Ndims>
+  bool operator<(
+      const SPHSmooth<ValType1, Ndims> &s1,
+      const SPHSmooth<ValType2, Ndims> &s2) {
+    return (s1.getSmoothingLen() < s2.getSmoothingLen());
+  }
+}; // namespace CosmoTool
--- a/src/yorick_nc3.cpp
+++ b/src/yorick_nc3.cpp
@ -260,7 +260,6 @@ namespace CosmoTool {
  template<typename T>
  void loadArray(const std::string& fname,
 		 T*&array, uint32_t *&dimList, uint32_t& rank)
-	throw (NoSuchFileException)
  {
    NcFile f(fname.c_str(), NcFile::ReadOnly);
Author	SHA1	Message	Date
Guilhem Lavaux	fb51ffd35e	Fixup for sph	2024-04-26 08:53:38 +02:00
Guilhem Lavaux	924047de22	Add possibility to have an external state	2024-04-25 18:11:06 +02:00
Guilhem Lavaux	be64c7fd7a	Remove support for old python	2024-04-25 17:06:15 +02:00
guilhem.lavaux@iap.fr	47d63a25ce	Bump version	2024-02-13 10:47:28 +01:00
guilhem.lavaux@iap.fr	4bcc5f3270	Add missing numpy_adaptors	2024-02-13 10:37:52 +01:00
guilhem.lavaux@iap.fr	1d59533b17	Add python in search path	2024-02-13 10:28:58 +01:00
guilhem.lavaux@iap.fr	1548fd8450	Add missing pyx in MANIFEST	2024-02-13 09:15:36 +01:00
guilhem.lavaux@iap.fr	1151d0c8b6	Add cython back to dependency	2024-02-13 08:25:14 +01:00
guilhem.lavaux@iap.fr	1db266b4ea	Add pyproject.toml	2024-02-13 08:19:54 +01:00
guilhem.lavaux@iap.fr	f2a5092cf1	Fix MANIFEST	2024-02-13 07:53:12 +01:00
guilhem.lavaux@iap.fr	54a59b5246	Fix manifest	2024-02-13 07:53:12 +01:00
Guilhem Lavaux	883c338c08	Fixed building	2024-01-17 14:00:09 +01:00
guilhem.lavaux@iap.fr	b8e60a7d33	Updating to 1.3.5	2024-01-16 07:50:50 +01:00
guilhem.lavaux@iap.fr	751d8a19a0	Remove pyfftw from requirements	2024-01-16 07:49:08 +01:00
guilhem.lavaux@iap.fr	97a1e34132	Strip mandatory pyfftw support	2024-01-16 07:33:45 +01:00
guilhem.lavaux@iap.fr	26e095fc71	Fix for default MacOS build with python	2024-01-16 07:21:56 +01:00
Guilhem Lavaux	6adf02b287	Attempt to fix linkage errors	2023-12-07 08:56:24 +01:00
Guilhem Lavaux	b878efb8b1	Fix requirements	2023-12-06 22:35:55 +01:00
Guilhem Lavaux	7c7ccd6f87	Changed to use configure and not CMake for hdf5.	2023-12-06 21:21:11 +00:00
Guilhem Lavaux	7a81120977	Fix setPeriodic	2023-12-06 18:40:06 +01:00
Guilhem Lavaux	c88f91ba80	Bump to 1.3.4	2023-12-06 09:29:58 +01:00
Guilhem Lavaux	0093a6aa0f	Fix cython syntax and port to more recent numpy.pxd	2023-12-06 09:28:43 +01:00
Guilhem Lavaux	4afc982dfc	Cache downloads	2023-12-06 09:28:43 +01:00
Guilhem Lavaux	d4b1742cbf	Add periodic argument	2023-12-05 16:35:54 +01:00
guilhem.lavaux@iap.fr	d2cd6bb650	Merge remote-tracking branch 'origin/master'	2023-12-03 20:30:25 +01:00
guilhem.lavaux@iap.fr	2aa7c96e48	Update for macmini	2023-12-03 20:29:56 +01:00
Guilhem Lavaux	d8b58bd8f1	Update version	2023-04-23 12:44:12 +02:00
Guilhem Lavaux	05ac03fb7a	Update bundled external libs	2023-04-22 17:21:14 +02:00
Guilhem Lavaux	522588fc1f	Fixup script for more modern pypa image	2023-02-27 14:16:36 +01:00
Guilhem Lavaux	39fe922143	Merge branch 'master' of bitbucket.org:glavaux/cosmotool	2023-02-04 09:54:41 +01:00
Guilhem Lavaux	fd16d7dc69	Add config.* updated tool for mac support	2023-02-04 09:54:32 +01:00
Guilhem Lavaux	5eb0fe210b	Fix version	2023-02-04 08:26:42 +01:00
Guilhem Lavaux	957d5187e9	Bump requirement of numpy	2023-02-03 21:44:28 +01:00
Guilhem Lavaux	f326962bc8	Add dependencies	2023-01-31 14:38:04 +01:00
Eleni	4509174d81	moved param to transfer func	2022-12-06 20:59:13 +01:00
Eleni	2cc11b3633	wiggle parameterization	2022-12-06 20:59:07 +01:00
Eleni	f79d7f6830	added BAO wiggles	2022-12-06 20:59:02 +01:00
Guilhem Lavaux	046e9a1447	Reformat and add some adjoint gradient	2022-11-17 21:50:06 +01:00
Guilhem Lavaux	fe06434619	Reformat	2022-11-16 16:48:40 +01:00
Guilhem Lavaux	4633f6edc9	Reformat	2022-11-16 16:48:31 +01:00
Guilhem Lavaux	b538d4974d	Add and fix support for parallel SPH state	2022-11-15 09:02:11 +01:00
Guilhem Lavaux	f03751907b	Push compatibility to HDF5 1.12	2022-08-17 08:57:52 +03:00
Guilhem Lavaux	31d0f19408	Fix patch	2022-06-19 10:53:45 +02:00
Guilhem Lavaux	cc94866bc3	Remove bind2nd from omptl	2022-06-19 08:24:26 +02:00
Ewoud	a309aa732b	Fix for C++17: switch out random_shuffle for shuffle. For details see https://clang.llvm.org/extra/clang-tidy/checks/modernize-replace-random-shuffle.html	2022-06-13 11:36:04 +02:00
Guilhem Lavaux	59bb99e7ee	Fix check of kdtree loading from disk	2022-05-27 16:26:59 +03:00
Guilhem Lavaux	e2a2c7287c	Remove obsolete exception specification	2022-02-10 07:20:21 +01:00
Guilhem Lavaux	d875416200	Tweaks for speed	2022-01-29 08:22:08 +01:00
Guilhem Lavaux	8ab094ad3d	Bump version	2022-01-27 16:28:16 +01:00
Guilhem Lavaux	0b77b010e4	Merge branch 'master' of bitbucket.org:glavaux/cosmotool	2022-01-27 13:36:25 +01:00
Guilhem Lavaux	c264f2c69d	Add output file argument	2022-01-27 13:36:19 +01:00
Guilhem Lavaux	533d8d0630	Merge branch 'master' of bitbucket.org:glavaux/cosmotool	2021-10-15 18:39:01 +02:00
Guilhem Lavaux	b01543a02a	Use 64 bits instead of 32 bits	2021-10-15 18:38:55 +02:00
Guilhem Lavaux	009fae2418	Bump version	2021-08-05 08:00:19 +02:00