.gitignore

@@ -1,9 +1,4 @@
 *~
-.eggs/
-dist/
-wheelhouse/
-cosmotool.egg-info/
-build/
 *.o
 *.prog
 *.pyc

CMakeLists.txt

@@ -42,17 +42,6 @@ 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)
@@ -66,6 +55,9 @@ 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)
@@ -78,8 +70,8 @@ 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 "3")
+SET(CPACK_PACKAGE_VERSION_MINOR "0")
-SET(CPACK_PACKAGE_VERSION_PATCH "4${EXTRA_VERSION}")
+SET(CPACK_PACKAGE_VERSION_PATCH "0${EXTRA_VERSION}")
 SET(CPACK_PACKAGE_INSTALL_DIRECTORY "CosmoToolbox-${CPACK_PACKAGE_VERSION_MAJOR}.${CPACK_PACKAGE_VERSION_MINOR}")
 SET(CPACK_STRIP_FILES "lib/libCosmoTool.so")
 SET(CPACK_SOURCE_IGNORE_FILES 

MANIFEST.in

@@ -1,5 +1,4 @@
 include  .gitignore
-include  pyproject.toml
 include  CMakeLists.txt
 include  FindNumPy.cmake
 include  FindPyLibs.cmake
@@ -17,17 +16,14 @@ 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-8d51946.tar.gz
+include  external/libsharp-6077806.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
@@ -90,14 +86,12 @@ include  sample/testSmooth.cpp
 include  sample/test_cosmopower.cpp
 include  sample/test_fft_calls.cpp
 include  sample/test_healpix_calls.cpp
-include  sample/test_special.cpp
 include  sample/testkd.cpp
 include  sample/testkd2.cpp
 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
@@ -173,4 +167,3 @@ include  src/tf_fit.hpp
 include  src/yorick.hpp
 include  src/yorick_nc3.cpp
 include  src/yorick_nc4.cpp
-include  src/special_math.hpp

builder/build-wheels.sh

@@ -8,32 +8,28 @@ 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 libffi-devel hdf5 hdf5-devel
+yum install -y cmake3 gsl-devel zlib-devel
 
 ln -fs /usr/bin/cmake3 /usr/bin/cmake
 
 
-test -d /io/wheelhouse || mkdir /io/wheelhouse
+ALL_PYTHON="cp36-cp36m cp37-cp37m cp38-cp38" 
-test -d /io/wheelhouse/fix || mkdir /io/wheelhouse/fix
-
-ALL_PYTHON="cp39-cp39 cp310-cp310"
 
 # Compile wheels
 for pkg in $ALL_PYTHON; do
     PYBIN=/opt/python/${pkg}/bin
 #    "${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 /io/wheelhouse/
+    "${PYBIN}/pip" wheel -vvv /io/ -w wheelhouse/
 done
 
 # Bundle external shared libraries into the wheels
-for whl in /io/wheelhouse/cosmotool*linux*.whl; do
+for whl in wheelhouse/*.whl; do
-    auditwheel repair "$whl" --plat $PLAT -w /io/wheelhouse/fix
+    auditwheel repair "$whl" --plat $PLAT -w /io/wheelhouse/
 done
 
 # Install packages and test
-#for pkg in $ALL_PYTHON; do
+for pkg in $ALL_PYTHON; do
-#    PYBIN=/opt/python/${pkg}/bin
+    PYBIN=/opt/python/${pkg}/bin
-#    "${PYBIN}/pip" install cosmotool --no-index -f /io/wheelhouse
+    "${PYBIN}/pip" install cosmotool --no-index -f /io/wheelhouse
-#done
+done

builder/start.sh

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

conda/build.sh

@@ -1,4 +1 @@
-export CC=$(basename ${CC})
-export CXX=$(basename ${CXX})
-
 $PYTHON setup.py install

conda/conda_build_config.yaml

@@ -1,10 +1,7 @@
 python:
-  - 3.9
-  - 3.8
   - 3.7
+  - 3.8
 
 numpy:
+  - 1.11
   - 1.19
-
-gsl:
-  - 2.4

conda/meta.yaml

@@ -1,50 +1,38 @@
 package:
   name: cosmotool
-  version: "1.3.0"
+  version: "1.0.0a7"
 
 source:
-  git_rev: a86c9a8
+  git_rev: 7fce73e
   git_url: https://bitbucket.org/glavaux/cosmotool
 
 requirements:
   build:
-    - python                                 # [build_platform != target_platform]
+    - numpy >=1.11
-    - cross-python_{{ target_platform }}     # [build_platform != target_platform]
+    - gcc_linux-64
-    - cython                                 # [build_platform != target_platform]
+    - gxx_linux-64
-    - numpy                                  # [build_platform != target_platform]
-    - {{ compiler('c') }}
-    - {{ compiler('cxx') }}
-    - llvm-openmp                            # [osx]
-    - libgomp                                # [linux]
-
-  host:
     - python
-    - pip
+    - setuptools
-    - numpy
-    - pkgconfig
-    - numexpr
     - cython
     - healpy
+    - numexpr
     - cffi
     - pyfftw
-    - gsl {{ gsl }}
+    - gsl
     - h5py
+
   run:
+    - numpy
     - python
-    - {{ pin_compatible('numpy') }}
     - healpy
     - numexpr
     - cffi
     - pyfftw
     - h5py
-    - {{ pin_compatible('gsl') }}
-    - llvm-openmp
 
 test:
   imports:
     - cosmotool
-  requires:
-    - pip
 
 about:
   home: https://bitbucket.org/glavaux/cosmotool

external/external_build.cmake

@@ -2,19 +2,16 @@ 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.12/hdf5-1.12.2/src/hdf5-1.12.2.tar.gz" CACHE STRING "URL to download HDF5 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 "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(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.1.tar.gz" CACHE STRING "URL to download NetCDF-C++ 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.82.0/source/boost_1_82_0.tar.gz" CACHE STRING "URL to download Boost from")
+SET(BOOST_URL "http://sourceforge.net/projects/boost/files/boost/1.61.0/boost_1_61_0.tar.gz/download" 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 ")
+SET(GSL_URL "ftp://ftp.gnu.org/gnu/gsl/gsl-1.15.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) 
@@ -39,13 +36,12 @@ CHECK_CHANGE_STATE(INTERNAL_DLIB DLIB_INCLUDE_DIR DLIB_LIBRARIES)
 
 IF(ENABLE_OPENMP)
   IF (NOT OPENMP_FOUND)
-    MESSAGE(NOTICE "No known compiler option for enabling OpenMP")
+    MESSAGE(ERROR "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)
 
 
@@ -59,25 +55,17 @@ 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-8d51946.tar.gz
+    URL ${CMAKE_SOURCE_DIR}/external/libsharp-6077806.tar.gz
     PREFIX ${BUILD_PREFIX}/sharp-prefix
-    DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
     BUILD_IN_SOURCE 1 
-    CONFIGURE_COMMAND
+    CONFIGURE_COMMAND autoconf && ./configure "CC=${CMAKE_C_COMPILER}" "CXX=${CMAKE_CXX_COMPILER}" --prefix=${DEP_BUILD} ${SHARP_OPENMP}
-          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}
   )
@@ -94,23 +82,27 @@ if (INTERNAL_HDF5)
   ExternalProject_Add(hdf5
     PREFIX ${BUILD_PREFIX}/hdf5-prefix
     URL ${HDF5_URL}
-    URL_HASH MD5=30172c75e436d7f2180e274071a4ca97
+    URL_HASH MD5=29117bf488887f89888f9304c8ebea0b
-    DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
+    CMAKE_ARGS
-    CONFIGURE_COMMAND
+       -DCMAKE_INSTALL_PREFIX=${EXT_INSTALL}
-        ${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
+       -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER}
-    INSTALL_COMMAND make install
+       -DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER}
+       -DHDF5_BUILD_CPP_LIB=ON
+       -DHDF5_BUILD_TOOLS=ON
+       -DHDF5_BUILD_HL_LIB=ON
+       -DBUILD_SHARED_LIBS=OFF
   )
   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 ${RT_LIBRARY}")
+  SET(CONFIGURE_LIBS "${CONFIGURE_LIBS} -ldl")
-  set(HDF5_C_STATIC_LIBRARY ${HDF5_BIN_DIR}/lib/libhdf5.a)
+  set(HDF5_C_STATIC_LIBRARY ${HDF5_BIN_DIR}/lib/libhdf5-static.a)
-  set(HDF5_HL_STATIC_LIBRARY ${HDF5_BIN_DIR}/lib/libhdf5_hl.a)
+  set(HDF5_HL_STATIC_LIBRARY ${HDF5_BIN_DIR}/lib/libhdf5_hl-static.a)
-  set(HDF5_LIBRARIES ${HDF5_BIN_DIR}/lib/libhdf5.a CACHE STRING "HDF5 lib" FORCE)
+  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.a CACHE STRING "HDF5 HL 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.a CACHE STRING "HDF5 C++ lib" FORCE)
+  set(HDF5_CXX_LIBRARIES ${HDF5_BIN_DIR}/lib/libhdf5_cpp-static.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)
 
@@ -168,52 +160,43 @@ 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-byterange --disable-cdmremote --disable-rpc --enable-cxx-4 
+         --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 |
+  LIST_SEPARATOR |
-    CMAKE_ARGS
+	CMAKE_ARGS
 		-DCMAKE_PREFIX_PATH=${CMAKE_PREFIX_PATH_ALT_SEP}
 		-DCMAKE_C_COMPILER=${CMAKE_C_COMPILER}
 		-DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER}
-                -DNC_EXTRA_DEPS=${RT_DEP}
 		-DBUILD_SHARED_LIBS=OFF
 		-DBUILD_TESTING=OFF
 		-DCMAKE_BUILD_TYPE=Release
 		-DENABLE_NETCDF4=ON
-                -DENABLE_BYTERANGE=FALSE
+    -DCMAKE_POSITION_INDEPENDENT_CODE=ON
-		-DCMAKE_POSITION_INDEPENDENT_CODE=ON
+    -DENABLE_DAP=OFF
-		-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}
@@ -251,8 +234,6 @@ if (INTERNAL_BOOST)
   ExternalProject_Add(boost
     URL ${BOOST_URL}
     PREFIX ${BUILD_PREFIX}/boost-prefix
-    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
@@ -295,7 +276,6 @@ 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
@@ -348,7 +328,6 @@ 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}
@@ -378,7 +357,6 @@ 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}
@@ -411,7 +389,6 @@ SET(cosmotool_DEPS ${cosmotool_DEPS} omptl)
 SET(OMPTL_BUILD_DIR ${BUILD_PREFIX}/omptl-prefix/src/omptl)
 ExternalProject_Add(omptl
   PREFIX ${BUILD_PREFIX}/omptl-prefix
-  DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
   URL ${CMAKE_SOURCE_DIR}/external/omptl-20120422.tar.bz2
   CONFIGURE_COMMAND echo "No configure"
   BUILD_COMMAND echo "No build"

external/patch-omptl

@@ -1,6 +1,6 @@
-diff -ur omptl.old/omptl_algorithm omptl/omptl_algorithm
+diff -ur omptl.orig/omptl_algorithm omptl/omptl_algorithm
---- omptl.old/omptl_algorithm	2022-06-19 08:21:39.815498672 +0200
+--- omptl.orig/omptl_algorithm	2017-01-16 14:58:37.996690639 +0100
-+++ omptl/omptl_algorithm	2022-06-19 08:21:52.953544672 +0200
++++ omptl/omptl_algorithm	2017-01-16 15:00:26.678641720 +0100
 @@ -20,7 +20,7 @@
  #define OMPTL_ALGORITHM 1
  
@@ -22,14 +22,12 @@ diff -ur omptl.old/omptl_algorithm omptl/omptl_algorithm
  #endif
  
  #endif /* OMPTL_ALGORITHM */
-diff -ur omptl.old/omptl_algorithm_par.h omptl/omptl_algorithm_par.h
+diff -ur omptl.orig/omptl_algorithm_par.h omptl/omptl_algorithm_par.h
---- omptl.old/omptl_algorithm_par.h	2022-06-19 08:21:39.816498675 +0200
+--- omptl.orig/omptl_algorithm_par.h	2017-01-16 14:58:37.996690639 +0100
-+++ omptl/omptl_algorithm_par.h	2022-06-19 08:23:50.705956964 +0200
++++ omptl/omptl_algorithm_par.h	2017-01-16 14:59:57.974126410 +0100
-@@ -20,9 +20,10 @@
+@@ -21,8 +21,8 @@
- #include <utility>
  #include <cmath>
  #include <cstdlib>
-+#include <random>
  
 -#include <omptl/omptl_tools.h>
 -#include <omptl/omptl_numeric>
@@ -38,136 +36,9 @@ diff -ur omptl.old/omptl_algorithm_par.h omptl/omptl_algorithm_par.h
  
  #include <iterator>
  
-@@ -510,7 +511,7 @@
+diff -ur omptl.orig/omptl_numeric omptl/omptl_numeric
- 
+--- omptl.orig/omptl_numeric	2017-01-16 14:58:37.996690639 +0100
- 		typename std::vector<Iterator>::iterator result =
++++ omptl/omptl_numeric	2017-01-16 15:00:57.051186974 +0100
- 		    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
- 
--	const unsigned max_depth = std::floor(std::tr1::log2(P));
-+	const unsigned max_depth = unsigned(std::floor(std::log2(P)));
- 	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	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
  
@@ -192,9 +63,9 @@ diff -ur omptl.old/omptl_numeric omptl/omptl_numeric
 +#include "omptl_numeric_extensions.h"
  
  #endif /* OMPTL_NUMERIC */
-diff -ur omptl.old/omptl_numeric_extensions.h omptl/omptl_numeric_extensions.h
+diff -ur omptl.orig/omptl_numeric_extensions.h omptl/omptl_numeric_extensions.h
---- omptl.old/omptl_numeric_extensions.h	2022-06-19 08:21:39.815498672 +0200
+--- omptl.orig/omptl_numeric_extensions.h	2017-01-16 14:58:37.996690639 +0100
-+++ omptl/omptl_numeric_extensions.h	2022-06-19 08:21:52.956544683 +0200
++++ omptl/omptl_numeric_extensions.h	2017-01-16 14:59:21.549472508 +0100
 @@ -51,9 +51,9 @@
  } // namespace
  
@@ -207,9 +78,9 @@ diff -ur omptl.old/omptl_numeric_extensions.h omptl/omptl_numeric_extensions.h
  #endif
  
  namespace omptl
-diff -ur omptl.old/omptl_numeric_par.h omptl/omptl_numeric_par.h
+diff -ur omptl.orig/omptl_numeric_par.h omptl/omptl_numeric_par.h
---- omptl.old/omptl_numeric_par.h	2022-06-19 08:21:39.816498675 +0200
+--- omptl.orig/omptl_numeric_par.h	2017-01-16 14:58:37.996690639 +0100
-+++ omptl/omptl_numeric_par.h	2022-06-19 08:21:52.957544686 +0200
++++ omptl/omptl_numeric_par.h	2017-01-16 14:59:36.397739066 +0100
 @@ -23,8 +23,8 @@
  #include <functional>
  #include <iterator>
@@ -221,15 +92,3 @@ 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	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>
- 
--#include <tr1/cmath>
-+#include <cmath>
- 
- namespace omptl
- {

pyproject.toml

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

python/CMakeLists.txt

@@ -2,7 +2,7 @@ set(CMAKE_SHARED_MODULE_PREFIX)
 
 
 set(PYTHON_INCLUDES ${NUMPY_INCLUDE_DIRS} ${PYTHON_INCLUDE_PATH} ${CMAKE_SOURCE_DIR}/python)
-include_directories(${CMAKE_SOURCE_DIR}/python ${CMAKE_SOURCE_DIR}/src ${CMAKE_BINARY_DIR}/src)
+include_directories(${CMAKE_SOURCE_DIR}/src ${CMAKE_BINARY_DIR}/src)
 
 IF(CYTHON)
   add_custom_command(
@@ -54,10 +54,6 @@ SET(CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -Bsymbolic-functions
 if(APPLE)
   set(CMAKE_MODULE_LINKER_FLAGS "-undefined dynamic_lookup")
 endif()
-IF(NOT APPLE)
-  set(CMAKE_MODULE_LINKER_FLAGS "-Wl,--version-script=${CMAKE_CURRENT_SOURCE_DIR}/cosmotool.version")
-ENDIF()
-
 
 target_link_libraries(_cosmotool PRIVATE ${CosmoTool_local}  ${GSL_LIBRARIES})
 target_link_libraries(_cosmo_power PRIVATE ${CosmoTool_local} ${GSL_LIBRARIES})
@@ -69,10 +65,13 @@ 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")
+  message(STATUS "Building bispectrum support (path = ${Boost_INCLUDE_DIRS})")
   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()
@@ -113,17 +112,8 @@ configure_file(${CMAKE_CURRENT_SOURCE_DIR}/cosmotool/config.py.in ${CMAKE_CURREN
 
 INSTALL(TARGETS 
   ${ct_TARGETS}
-  COMPONENT python
   LIBRARY DESTINATION ${PYTHON_SITE_PACKAGES}/cosmotool
 )
 
-INSTALL(DIRECTORY cosmotool ${CMAKE_CURRENT_BINARY_DIR}/cosmotool
+INSTALL(DIRECTORY cosmotool ${CMAKE_CURRENT_BINARY_DIR}/cosmotool DESTINATION ${PYTHON_SITE_PACKAGES} 
-        COMPONENT python DESTINATION ${PYTHON_SITE_PACKAGES}
         FILES_MATCHING PATTERN "*.py")
-
-add_custom_target(
-    python-install
-    DEPENDS ${ct_TARGETS}
-    COMMAND "${CMAKE_COMMAND}" -DCMAKE_INSTALL_COMPONENT=python -P
-            "${CMAKE_BINARY_DIR}/cmake_install.cmake")
-

python/_cosmomath.pyx

@@ -1,4 +1,3 @@
-#cython: language_level=3
 import numpy as np
 cimport numpy as np
 
@@ -8,21 +7,15 @@ np.import_ufunc()
 cdef extern from "sys/types.h":
    ctypedef np.int64_t int64_t
 
-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) except + nogil
+   T log_modified_bessel_first_kind[T](T v, T z) nogil except +
-
-cdef extern from "numpy_adaptors.hpp" namespace "CosmoTool":
-   void parallel_ufunc_dd_d[T,IT](char **args, IT* dimensions, IT* steps, void *func)
 
 
 cdef np.PyUFuncGenericFunction loop_func[1]
 cdef char input_output_types[3]
 cdef void *elementwise_funcs[1]
 
-loop_func[0] = <np.PyUFuncGenericFunction>parallel_ufunc_dd_d[double,npy_intp]
+loop_func[0] = np.PyUFunc_dd_d
 
 input_output_types[0] = np.NPY_DOUBLE
 input_output_types[1] = np.NPY_DOUBLE

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) except + nogil
+  SimuData *loadGadgetMulti(const char *fname, int id, int flags, int gformat) nogil except +
   void cxx_writeGadget "CosmoTool::writeGadget" (const char * s, SimuData *data) except +
 
 cdef extern from "safe_gadget.hpp":
@@ -247,9 +247,6 @@ class PySimulationAdaptor(PySimulationBase):
   def __init__(self,sim):
     self.simu = sim
 
-  def getNumParticles(self):
-    return self.simu.numParticles
-
   def getBoxsize(self):
     return self.simu.BoxSize
     
@@ -313,13 +310,6 @@ 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
@@ -327,8 +317,7 @@ 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
@@ -421,7 +410,7 @@ def loadGadget(str filename, int snapshot_id, int gadgetFormat = 1, bool loadPos
     with nogil:
       data = loadGadgetMulti(filename_bs, snapshot_id, flags, gadgetFormat)
     if data == <SimuData*>0:
-       raise RuntimeError("File could not be read")
+        return None
 
     return PySimulationAdaptor(wrap_simudata(data, flags))
 

python/_project.pyx

@@ -533,13 +533,12 @@ cdef void INTERNAL_project_cic_with_mass_periodic(DTYPE_t[:,:,:] g,
 
        
 def project_cic(npx.ndarray[DTYPE_t, ndim=2] x not None, npx.ndarray[DTYPE_t, ndim=1] mass, int Ngrid,
-                double Lbox, bool periodic = False, centered=True, output=None):
+                double Lbox, bool periodic = False, centered=True):
     """
-    project_cic(x array (N,3), mass (may be None), Ngrid, Lbox, periodict, centered=True, output=None)
+    project_cic(x array (N,3), mass (may be None), Ngrid, Lbox, periodict, centered=True)
 
     This function does a Cloud-In-Cell projection of a 3d unstructured dataset. First argument is a Nx3 array of coordinates. 
     Second argument is an optinal mass. Ngrid is the size output grid and Lbox is the physical size of the grid.
-    if output is not None, it must be a numpy array with dimension NgridxNgridxNgrid. The result will be accumulated in that array.
     """
     cdef npx.ndarray[DTYPE_t, ndim=3] g
     cdef double shifter
@@ -559,13 +558,7 @@ def project_cic(npx.ndarray[DTYPE_t, ndim=2] x not None, npx.ndarray[DTYPE_t, nd
     if mass is not None and mass.shape[0] != x.shape[0]:
         raise ValueError("Mass array and coordinate array must have the same number of elements")
 
-    if output is None:
+    g = np.zeros((Ngrid,Ngrid,Ngrid),dtype=DTYPE)
-        g = np.zeros((Ngrid,Ngrid,Ngrid),dtype=DTYPE)
-    else:
-        if type(output) != np.ndarray:
-            raise ValueError("Invalid array type")
-        g = output
-
     cdef DTYPE_t[:,:,:] d_g = g
     cdef DTYPE_t[:,:] d_x = x
 
@@ -727,7 +720,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]

python/cosmotool.version

@@ -1,7 +0,0 @@
-CODEABI_1.0 {
-  global:
-    PyInit_*;
-    _init;
-    _fini;
-  local: *;
-};

requirements.txt

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

sample/CMakeLists.txt

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

sample/simple3DFilter.cpp

@@ -1,22 +1,22 @@
-#include "hdf5_array.hpp"
-#include "miniargs.hpp"
-#include "mykdtree.hpp"
-#include "omptl/algorithm"
 #include "openmp.hpp"
-#include "sphSmooth.hpp"
+#include "omptl/algorithm"
-#include "yorick.hpp"
-#include <H5Cpp.h>
-#include <boost/bind.hpp>
-#include <boost/format.hpp>
 #include <cassert>
+#include "yorick.hpp"
+#include "sphSmooth.hpp"
+#include "mykdtree.hpp"
+#include "miniargs.hpp"
+#include <H5Cpp.h>
+#include "hdf5_array.hpp"
 #include <iostream>
+#include <boost/format.hpp>
+#include <boost/bind.hpp>
 
 using namespace std;
 using namespace CosmoTool;
 
 #define N_SPH 32
 
-struct VCoord {
+struct VCoord{
   float v[3];
   float mass;
 };
@@ -27,87 +27,92 @@ 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>
+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,
-                              double cx, double cy, double cz,
+                              array3_type& bins, array3_type& arr, FuncT func, double rLimit2)
-                              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 collapse(3) schedule(dynamic)
+#pragma omp for schedule(dynamic) 
-    for (int rz = 0; rz < Nres; rz++) {
+    for (int rz = 0; rz < Nres; rz++)
+      {
+        double pz = (rz)*boxsize/Nres-cz;
 
-      for (int ry = 0; ry < Nres; ry++) {
+        cout << format("[%d] %d / %d") % smp_get_thread_id() % rz % Nres << endl;
-        for (int rx = 0; rx < Nres; rx++) {
+        for (int ry = 0; ry < Nres; ry++)
-          if (rz > rz_max) {
+          {
-            rz_max = rz;
+            double py = (ry)*boxsize/Nres-cy;
-            cout << format("[%d] %d / %d") % smp_get_thread_id() % rz % Nres
+            for (int rx = 0; rx < Nres; rx++)
-                 << endl;
+              {
+                double px = (rx)*boxsize/Nres-cx;
+                
+                MyTree::coords c = { float(px), float(py), float(pz) };
+
+                double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2];
+                if (r2 > rLimit2)
+                  {
+                    arr[rx][ry][rz] = 0;
+                    continue;
+                  }
+
+                uint32_t numInCell = bins[rx][ry][rz];
+                if (numInCell > N_SPH)
+                  smooth1.fetchNeighbours(c, numInCell);
+                else
+                  smooth1.fetchNeighbours(c);
+
+                arr[rx][ry][rz] = smooth1.computeSmoothedValue(c, func);
+              }
           }
-          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)};
-
-          double r2 = c[0] * c[0] + c[1] * c[1] + c[2] * c[2];
-          if (r2 > rLimit2) {
-            arr[rx][ry][rz] = 0;
-            continue;
-          }
-
-          uint32_t numInCell = bins[rx][ry][rz];
-          if (numInCell > N_SPH)
-            smooth1.fetchNeighbours(c, numInCell);
-          else
-            smooth1.fetchNeighbours(c);
-
-          arr[rx][ry][rz] = smooth1.computeSmoothedValue(c, func);
-        }
       }
-    }
   }
 }
 
-int main(int argc, char **argv) {
+int main(int argc, char **argv)
-  char *fname1, *outFile;
+{
+
+  char *fname1, *fname2;
   double rLimit, boxsize, rLimit2, cx, cy, cz;
   int Nres;
-  int periodic;
 
-  MiniArgDesc args[] = {{"INPUT DATA1", &fname1, MINIARG_STRING},
+  MiniArgDesc args[] = {
-                        {"RADIUS LIMIT", &rLimit, MINIARG_DOUBLE},
+    { "INPUT DATA1", &fname1, MINIARG_STRING },
-                        {"BOXSIZE", &boxsize, MINIARG_DOUBLE},
+    { "RADIUS LIMIT", &rLimit, MINIARG_DOUBLE },
-                        {"RESOLUTION", &Nres, MINIARG_INT},
+    { "BOXSIZE", &boxsize, MINIARG_DOUBLE },
-                        {"CX", &cx, MINIARG_DOUBLE},
+    { "RESOLUTION", &Nres, MINIARG_INT },
-                        {"CY", &cy, MINIARG_DOUBLE},
+    { "CX", &cx, MINIARG_DOUBLE },
-                        {"CZ", &cz, MINIARG_DOUBLE},
+    { "CY", &cy, MINIARG_DOUBLE },
-                        {"OUTPUT FILE", &outFile, MINIARG_STRING},
+    { "CZ", &cz, MINIARG_DOUBLE },
-                        {"PERIODIC", &periodic, MINIARG_INT},
+    { 0, 0, MINIARG_NULL }
-                        {0, 0, MINIARG_NULL}};
+  };
 
   if (!parseMiniArgs(argc, argv, args))
     return 1;
 
   H5::H5File in_f(fname1, 0);
-  H5::H5File out_f(outFile, H5F_ACC_TRUNC);
+  H5::H5File out_f("fields.h5", H5F_ACC_TRUNC);
   array_type v1_data;
-  uint64_t N1_points, N2_points;
+  uint32_t N1_points, N2_points;
   
   array3_type bins(boost::extents[Nres][Nres][Nres]);
 
-  rLimit2 = rLimit * rLimit;
+  rLimit2 = rLimit*rLimit;
 
   hdf5_read_array(in_f, "particles", v1_data);
   assert(v1_data.shape()[1] == 7);
@@ -119,93 +124,94 @@ int main(int argc, char **argv) {
   MyCell *allCells_1 = new MyCell[N1_points];
   
 #pragma omp parallel for schedule(static)
-  for (uint32_t i = 0; i < Nres * Nres * Nres; i++)
+  for (long i = 0; i < Nres*Nres*Nres; i++)
     bins.data()[i] = 0;
 
   cout << "Shuffling data in cells..." << endl;
 #pragma omp parallel for schedule(static)
-  for (uint64_t i = 0; i < N1_points; i++) {
+  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 (int j = 0; j < 3; j++)
-    for (int k = 0; k < 3; k++)
+        allCells_1[i].coord[j] = v1_data[i][j];
-      allCells_1[i].val.pValue.v[k] = v1_data[i][3 + k];
+      for (int k = 0; k < 3; k++)
-    allCells_1[i].val.pValue.mass = v1_data[i][6];
+        allCells_1[i].val.pValue.v[k] = v1_data[i][3+k];
-    allCells_1[i].active = true;
+      allCells_1[i].val.pValue.mass = v1_data[i][6];
-    allCells_1[i].val.weight = 0.0;
+      allCells_1[i].active = true;
+      allCells_1[i].val.weight = 0.0;
 
-    long rx = floor((allCells_1[i].coord[0] + cx) * Nres / boxsize + 0.5);
+      long rx = floor((allCells_1[i].coord[0]+cx)*Nres/boxsize+0.5);
-    long ry = floor((allCells_1[i].coord[1] + cy) * Nres / boxsize + 0.5);
+      long ry = floor((allCells_1[i].coord[1]+cy)*Nres/boxsize+0.5);
-    long rz = floor((allCells_1[i].coord[2] + cz) * Nres / boxsize + 0.5);
+      long rz = floor((allCells_1[i].coord[2]+cz)*Nres/boxsize+0.5);
       
-    if (rx < 0 || rx >= Nres || ry < 0 || ry >= Nres || rz < 0 || rz >= Nres)
+      if (rx < 0 || rx >= Nres || ry < 0 || ry >= Nres || rz < 0 || rz >= Nres)
-      continue;
+        continue;
 	  
-    auto &b = bins[rx][ry][rz];
+//#pragma omp atomic update
-#pragma omp atomic
+      bins[rx][ry][rz]++;
-    b++;
+    }
-  }
   v1_data.resize(boost::extents[1][1]);
   
   hdf5_write_array(out_f, "num_in_cell", bins);
 
   cout << "Building trees..." << endl;
   MyTree tree1(allCells_1, N1_points);
-  tree1.setPeriodic(periodic != 0, boxsize);
 
   cout << "Creating smoothing filter..." << endl;
 
-  //  array3_type out_rad_1(boost::extents[Nres][Nres][Nres]);
+//  array3_type out_rad_1(boost::extents[Nres][Nres][Nres]);
  
   cout << "Weighing..." << endl;
 
-  int rz_max = 0;
+#pragma omp parallel
-#pragma omp parallel shared(rz_max)
   {
     MySmooth smooth1(&tree1, N_SPH);
     
-#pragma omp for collapse(3) schedule(dynamic, 8)
+#pragma omp for schedule(dynamic) 
-    for (int rz = 0; rz < Nres; rz++) {
+    for (int rz = 0; rz < Nres; rz++)
-      for (int ry = 0; ry < Nres; ry++) {
+     {
-        for (int rx = 0; rx < Nres; rx++) {
+        double pz = (rz)*boxsize/Nres-cz;
-          if (rz > rz_max) {
+
-            rz_max = rz;
+        (cout << rz << " / " << Nres << endl).flush();
-            (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++)
+              {
+                double px = (rx)*boxsize/Nres-cx;
+                
+                MyTree::coords c = { float(px), float(py), float(pz) };
+
+                double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2];
+                if (r2 > rLimit2)
+                  {
+                    continue;
+                  }
+
+                uint32_t numInCell = bins[rx][ry][rz];
+                if (numInCell > N_SPH)
+                  smooth1.fetchNeighbours(c, numInCell);
+                else
+                  smooth1.fetchNeighbours(c);
+#pragma omp critical
+                smooth1.addGridSite(c);
+              }
           }
-          double pz = (rz)*boxsize / Nres - cz;
+        (cout << " Done " << rz << endl).flush();
-          double py = (ry)*boxsize / Nres - cy;
-          double px = (rx)*boxsize / Nres - cx;
-
-          MyTree::coords c = {float(px), float(py), float(pz)};
-
-          double r2 = c[0] * c[0] + c[1] * c[1] + c[2] * c[2];
-          if (r2 > rLimit2) {
-            continue;
-          }
-
-          uint32_t numInCell = bins[rx][ry][rz];
-          if (numInCell > N_SPH)
-            smooth1.fetchNeighbours(c, numInCell);
-          else
-            smooth1.fetchNeighbours(c);
-          smooth1.addGridSite(c);
-        }
       }  
-    }
   }
   
   cout << "Interpolating..." << endl;
 
   array3_type interpolated(boost::extents[Nres][Nres][Nres]);
   
-  computeInterpolatedField(&tree1, boxsize, Nres, cx, cy, cz, bins,
+  computeInterpolatedField(&tree1, boxsize, Nres, cx, cy, cz,
-                           interpolated, getMass, rLimit2);
+                           bins, interpolated, getMass, rLimit2);
   hdf5_write_array(out_f, "density", interpolated);
-  // out_f.flush();
+  //out_f.flush();
   for (int i = 0; i < 3; i++) {
-    computeInterpolatedField(&tree1, boxsize, Nres, cx, cy, cz, bins,
+      computeInterpolatedField(&tree1, boxsize, Nres, cx, cy, cz,
-                             interpolated, boost::bind(getVelocity, _1, i),
+                               bins, interpolated, boost::bind(getVelocity, _1, i), rLimit2);
-                             rLimit2);
+      hdf5_write_array(out_f, str(format("p%d") % i), interpolated);
-    hdf5_write_array(out_f, str(format("p%d") % i), interpolated);
   }
   
   return 0;

setup.py

@@ -2,7 +2,6 @@ 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
@@ -168,33 +167,15 @@ 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", 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,
+                    '-DCMAKE_C_COMPILER=' + os.environ["CC"], "-DCMAKE_CXX_COMPILER=" + os.environ["CXX"],
                     '-DENABLE_OPENMP=ON','-DINTERNAL_BOOST=ON','-DINTERNAL_EIGEN=ON',
                     '-DINTERNAL_HDF5=ON','-DINTERNAL_NETCDF=ON',
                     '-DBUILD_PYTHON=ON', '-DINSTALL_PYTHON_LOCAL=OFF',
                     '-DCOSMOTOOL_PYTHON_PACKAGING=ON',
                     f"-DCYTHON={cython_code}",
                     f"-DPYTHON_SITE_PACKAGES={build_dir.absolute()}/private_install",
-                    f"-DPYTHON_EXECUTABLE={sys.executable}"] + compilers)
+                    f"-DPYTHON_EXECUTABLE={sys.executable}"])
 
 
         self.announce("Building binaries", level=3)
@@ -237,10 +218,10 @@ class BuildCMakeExt(build_ext):
 CosmoTool_extension = CMakeExtension(name="cosmotool")
 
 setup(name='cosmotool',
-      version='1.3.6',
+      version='1.1.0',
       packages=["cosmotool"],
       package_dir={'cosmotool': 'python/cosmotool'},
-      install_requires=['cython','numpy','cffi','numexpr','h5py'],
+      install_requires=['numpy','cffi','numexpr','pyfftw','h5py'],
       setup_requires=['cython','cffi','numpy','numexpr'],
       ext_modules=[CosmoTool_extension],
       description='A small cosmotool box of useful functions',

src/CMakeLists.txt

@@ -38,7 +38,6 @@ if (HDF5_FOUND)
     h5_readFlash.cpp
     loadFlash.cpp
     )
-  add_dependencies(CosmoHDF5 ${cosmotool_DEPS})
   set_property(TARGET CosmoHDF5 PROPERTY POSITION_INDEPENDENT_CODE ${BUILD_SHARED_LIBS})
   target_include_directories(CosmoHDF5 BEFORE PRIVATE ${HDF5_INCLUDE_DIR})
 else(HDF5_FOUND)

src/cosmopower.cpp

@@ -212,30 +212,6 @@ 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
@@ -273,19 +249,6 @@ 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)));
@@ -379,14 +342,12 @@ void CosmoPower::normalize(double k_min, double k_max)
 
   normPower = 1;
 
-#if 0
   ofstream ff("PP_k.txt");
   for (int i = 0; i < 100; i++)
     {
        double k = pow(10.0, 8.0*i/100.-4);
        ff << k << " " << power(k) << endl;
     }
-#endif
 
  // gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &normVal, &abserr);
   gsl_integration_qag(&f, x_min, x_max, 0, TOLERANCE, NUM_ITERATION, GSL_INTEG_GAUSS61, w, &normVal, &abserr);
@@ -481,18 +442,12 @@ 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();
     }

src/cosmopower.hpp

@@ -89,9 +89,7 @@ namespace CosmoTool {
       POWER_SUGIYAMA,
       POWER_BDM,
       POWER_TEST,
-      HU_WIGGLES_ORIGINAL,
+      HU_WIGGLES_ORIGINAL
-      SAMPLE_WIGGLES,
-      BAO_TK
     };
 
     CosmoPower();
@@ -124,9 +122,6 @@ namespace CosmoTool {
     double powerBDM(double k);
     double powerTest(double k);
     double powerHuWigglesOriginal(double k);
-    double sample_BAO(double k);
-    double noWiggleTk(double k);
-    double BAO_Tk(double k);
   };
 
 };

src/fortran.cpp

@@ -139,22 +139,14 @@ void UnformattedRead::beginCheckpoint(bool bufferRecord)
 
 void UnformattedRead::endCheckpoint(bool autodrop)
 {
-  bool always_fail = false;
-
   if (recordBuffer != 0) {
     delete[] recordBuffer;
     recordBuffer = 0;
   }
-  if (cSize == Check_32bits) {
-    if (checkPointAccum >= 1UL<<32UL) {
-      always_fail = true;
-      checkPointAccum %= (1UL<<32UL);
-    }
-  } 
 
   if (checkPointRef != checkPointAccum)
     {
-      if (always_fail || !autodrop || checkPointAccum > checkPointRef) {
+      if (!autodrop || checkPointAccum > checkPointRef) {
 	throw InvalidUnformattedAccess();
       }
       f->seekg(checkPointRef-checkPointAccum, ios::cur);

src/fourier/fft/fftw_calls.hpp

@@ -39,112 +39,91 @@ knowledge of the CeCILL license and that you accept its terms.
 #include <fftw3.h>
 #include <complex>
 
-namespace CosmoTool {
+namespace CosmoTool
+{
 
-  static inline void init_fftw_wisdom() {
+static inline void init_fftw_wisdom()
-    fftw_import_system_wisdom();
+{
-    fftw_import_wisdom_from_filename("fft_wisdom");
+  fftw_import_system_wisdom();
-  }
+  fftw_import_wisdom_from_filename("fft_wisdom");
+}
 
-  static inline void save_fftw_wisdom() {
+static inline void save_fftw_wisdom()
-    fftw_export_wisdom_to_filename("fft_wisdom");
+{
-  }
+  fftw_export_wisdom_to_filename("fft_wisdom");
+}
 
-  template <typename T>
+template<typename T> class FFTW_Calls {};
-  class FFTW_Calls {};
 
-#define FFTW_CALLS_BASE(rtype, prefix)                                         \
-  template <>                                                                  \
-  class FFTW_Calls<rtype> {                                                    \
-  public:                                                                      \
-    typedef rtype real_type;                                                   \
-    typedef prefix##_complex complex_type;                                     \
-    typedef prefix##_plan plan_type;                                           \
-                                                                               \
-    static complex_type *alloc_complex(size_t N) {                             \
-      return prefix##_alloc_complex(N);                                        \
-    }                                                                          \
-    static real_type *alloc_real(size_t N) { return prefix##_alloc_real(N); }  \
-    static void free(void *p) { fftw_free(p); }                                \
-                                                                               \
-    static void execute(plan_type p) { prefix##_execute(p); }                  \
-    static void execute_r2c(plan_type p, real_type *in, complex_type *out) {   \
-      prefix##_execute_dft_r2c(p, in, out);                                    \
-    }                                                                          \
-    static void execute_c2r(plan_type p, complex_type *in, real_type *out) {   \
-      prefix##_execute_dft_c2r(p, in, out);                                    \
-    }                                                                          \
-    static void                                                                \
-    execute_r2c(plan_type p, real_type *in, std::complex<real_type> *out) {    \
-      prefix##_execute_dft_r2c(p, in, (complex_type *)out);                    \
-    }                                                                          \
-    static void                                                                \
-    execute_c2r(plan_type p, std::complex<real_type> *in, real_type *out) {    \
-      prefix##_execute_dft_c2r(p, (complex_type *)in, out);                    \
-    }                                                                          \
-    static void execute_c2c(                                                   \
-        plan_type p, std::complex<real_type> *in,                              \
-        std::complex<real_type> *out) {                                        \
-      prefix##_execute_dft(p, (complex_type *)in, (complex_type *)out);        \
-    }                                                                          \
-    static plan_type plan_dft_r2c_1d(                                          \
-        int Nx, real_type *in, complex_type *out, unsigned flags) {            \
-      return prefix##_plan_dft_r2c_1d(Nx, in, out, flags);                     \
-    }                                                                          \
-    static plan_type plan_dft_c2r_1d(                                          \
-        int Nx, complex_type *in, real_type *out, unsigned flags) {            \
-      return prefix##_plan_dft_c2r_1d(Nx, in, out, flags);                     \
-    }                                                                          \
-    static plan_type plan_dft_r2c_2d(                                          \
-        int Nx, int Ny, real_type *in, complex_type *out, unsigned flags) {    \
-      return prefix##_plan_dft_r2c_2d(Nx, Ny, in, out, flags);                 \
-    }                                                                          \
-    static plan_type plan_dft_c2r_2d(                                          \
-        int Nx, int Ny, complex_type *in, real_type *out, unsigned flags) {    \
-      return prefix##_plan_dft_c2r_2d(Nx, Ny, in, out, flags);                 \
-    }                                                                          \
-    static plan_type plan_dft_r2c_3d(                                          \
-        int Nx, int Ny, int Nz, real_type *in, complex_type *out,              \
-        unsigned flags) {                                                      \
-      return prefix##_plan_dft_r2c_3d(Nx, Ny, Nz, in, out, flags);             \
-    }                                                                          \
-    static plan_type plan_dft_c2r_3d(                                          \
-        int Nx, int Ny, int Nz, complex_type *in, real_type *out,              \
-        unsigned flags) {                                                      \
-      return prefix##_plan_dft_c2r_3d(Nx, Ny, Nz, in, out, flags);             \
-    }                                                                          \
-                                                                               \
-    static plan_type plan_dft_r2c(                                             \
-        int rank, const int *n, real_type *in, complex_type *out,              \
-        unsigned flags) {                                                      \
-      return prefix##_plan_dft_r2c(rank, n, in, out, flags);                   \
-    }                                                                          \
-    static plan_type plan_dft_c2r(                                             \
-        int rank, const int *n, complex_type *in, real_type *out,              \
-        unsigned flags) {                                                      \
-      return prefix##_plan_dft_c2r(rank, n, in, out, flags);                   \
-    }                                                                          \
-    static plan_type plan_dft_3d(                                              \
-        int Nx, int Ny, int Nz, complex_type *in, complex_type *out, int sign, \
-        unsigned flags) {                                                      \
-      return prefix##_plan_dft_3d(Nx, Ny, Nz, in, out, sign, flags);           \
-    }                                                                          \
-    static plan_type plan_dft_2d(                                              \
-        int Nx, int Ny, complex_type *in, complex_type *out, int sign,         \
-        unsigned flags) {                                                      \
-      return prefix##_plan_dft_2d(Nx, Ny, in, out, sign, flags);               \
-    }                                                                          \
-    static plan_type plan_dft_1d(                                              \
-        int Nx, complex_type *in, complex_type *out, int sign,                 \
-        unsigned flags) {                                                      \
-      return prefix##_plan_dft_1d(Nx, in, out, sign, flags);                   \
-    }                                                                          \
-    static void destroy_plan(plan_type plan) { prefix##_destroy_plan(plan); }  \
-  }
 
-  FFTW_CALLS_BASE(double, fftw);
+#define FFTW_CALLS_BASE(rtype, prefix) \
-  FFTW_CALLS_BASE(float, fftwf);
+  template<>				\
+class FFTW_Calls<rtype> {		\
+public: \
+  typedef rtype real_type; \
+  typedef prefix ## _complex complex_type; \
+  typedef prefix ## _plan plan_type; \
+  \
+  static complex_type *alloc_complex(size_t N) { return prefix ## _alloc_complex(N); } \
+  static real_type *alloc_real(size_t N) { return prefix ## _alloc_real(N); } \
+  static void free(void *p) { fftw_free(p); } \
+\
+  static void execute(plan_type p) { prefix ## _execute(p); } \
+  static void execute_r2c(plan_type p, real_type *in, complex_type *out) { prefix ## _execute_dft_r2c(p, in, out); } \
+  static void execute_c2r(plan_type p, complex_type *in, real_type *out) { prefix ## _execute_dft_c2r(p, in, out); } \
+  static void execute_r2c(plan_type p, real_type *in, std::complex<real_type> *out) { prefix ## _execute_dft_r2c(p, in, (complex_type*)out); } \
+  static void execute_c2r(plan_type p, std::complex<real_type> *in, real_type *out) { prefix ## _execute_dft_c2r(p, (complex_type*) in, out); } \
+  static void execute_c2c(plan_type p, std::complex<real_type> *in, std::complex<real_type> *out) { prefix ## _execute_dft(p, (complex_type *)in, (complex_type*)out); } \
+  static plan_type plan_dft_r2c_2d(int Nx, int Ny,  \
+		       real_type *in, complex_type *out, \
+		       unsigned flags) \
+  { \
+    return prefix ## _plan_dft_r2c_2d(Nx, Ny, in, out, \
+				flags); \
+  } \
+  static plan_type plan_dft_c2r_2d(int Nx, int Ny,  \
+                       complex_type *in, real_type *out, \
+                       unsigned flags) \
+  { \
+    return prefix ## _plan_dft_c2r_2d(Nx, Ny, in, out, \
+                                flags); \
+  } \
+  static plan_type plan_dft_r2c_3d(int Nx, int Ny, int Nz, \
+                                   real_type *in, complex_type *out, \
+                                   unsigned flags) \
+  { \
+    return prefix ## _plan_dft_r2c_3d(Nx, Ny, Nz, in, out, flags); \
+  } \
+  static plan_type plan_dft_c2r_3d(int Nx, int Ny, int Nz, \
+                                   complex_type *in, real_type *out, \
+                                   unsigned flags) \
+  { \
+    return prefix ## _plan_dft_c2r_3d(Nx, Ny, Nz, in, out, flags); \
+  } \
+\
+  static plan_type plan_dft_r2c(int rank, const int *n, real_type *in, \
+                                complex_type *out, unsigned flags) \
+  { \
+    return prefix ## _plan_dft_r2c(rank, n, in, out, flags); \
+  } \
+  static plan_type plan_dft_c2r(int rank, const int *n, complex_type *in, \
+                                real_type *out, unsigned flags) \
+  { \
+    return prefix ## _plan_dft_c2r(rank, n, in, out, flags); \
+  } \
+  static plan_type plan_dft_3d(int Nx, int Ny, int Nz, complex_type *in, complex_type *out, int sign, unsigned flags) { \
+    return prefix ## _plan_dft_3d(Nx, Ny, Nz, in, out, sign, flags); \
+  } \
+  static plan_type plan_dft_2d(int Nx, int Ny, complex_type *in, complex_type *out, int sign, unsigned flags) { \
+    return prefix ## _plan_dft_2d(Nx, Ny, in, out, sign, flags); \
+  } \
+  static void destroy_plan(plan_type plan) { prefix ## _destroy_plan(plan); } \
+}
+
+
+FFTW_CALLS_BASE(double, fftw);
+FFTW_CALLS_BASE(float, fftwf);
 #undef FFTW_CALLS_BASE
-}; // namespace CosmoTool
+};
 
 #endif

src/fourier/fft/fftw_calls_mpi.hpp

@@ -5,140 +5,110 @@
 #include <mpi.h>
 #include <fftw3-mpi.h>
 
-namespace CosmoTool {
+namespace CosmoTool
+{
 
-  static inline void init_fftw_mpi() { fftw_mpi_init(); }
+static inline void init_fftw_mpi()
+{
+  fftw_mpi_init();
+}
 
-  static inline void done_fftw_mpi() { fftw_mpi_cleanup(); }
+static inline void done_fftw_mpi()
+{
+  fftw_mpi_cleanup();
+}
 
-  template <typename T>
+template<typename T> class FFTW_MPI_Calls {};
-  class FFTW_MPI_Calls {};
 
-#define FFTW_MPI_CALLS_BASE(rtype, prefix)                                     \
-  template <>                                                                  \
-  class FFTW_MPI_Calls<rtype> {                                                \
-  public:                                                                      \
-    typedef rtype real_type;                                                   \
-    typedef prefix##_complex complex_type;                                     \
-    typedef prefix##_plan plan_type;                                           \
-                                                                               \
-    static complex_type *alloc_complex(size_t N) {                             \
-      return prefix##_alloc_complex(N);                                        \
-    }                                                                          \
-    static real_type *alloc_real(size_t N) { return prefix##_alloc_real(N); }  \
-    static void free(void *p) { fftw_free(p); }                                \
-                                                                               \
-    template <size_t Nd>                                                       \
-    static ptrdiff_t local_size(                                               \
-        std::array<ptrdiff_t, Nd> const &N, MPI_Comm comm,                     \
-        ptrdiff_t *local_n0, ptrdiff_t *local_0_start) {                       \
-      return prefix##_mpi_local_size(                                          \
-          Nd, N.data(), comm, local_n0, local_0_start);                        \
-    }                                                                          \
-    static ptrdiff_t local_size_2d(                                            \
-        ptrdiff_t N0, ptrdiff_t N1, MPI_Comm comm, ptrdiff_t *local_n0,        \
-        ptrdiff_t *local_0_start) {                                            \
-      return prefix##_mpi_local_size_2d(                                       \
-          N0, N1, comm, local_n0, local_0_start);                              \
-    }                                                                          \
-                                                                               \
-    static ptrdiff_t local_size_3d(                                            \
-        ptrdiff_t N0, ptrdiff_t N1, ptrdiff_t N2, MPI_Comm comm,               \
-        ptrdiff_t *local_n0, ptrdiff_t *local_0_start) {                       \
-      return prefix##_mpi_local_size_3d(                                       \
-          N0, N1, N2, comm, local_n0, local_0_start);                          \
-    }                                                                          \
-                                                                               \
-    static void execute(plan_type p) { prefix##_execute(p); }                  \
-    static void                                                                \
-    execute_c2c(plan_type p, complex_type *in, complex_type *out) {            \
-      prefix##_mpi_execute_dft(p, in, out);                                    \
-    }                                                                          \
-    static void execute_c2c(                                                   \
-        plan_type p, std::complex<real_type> *in,                              \
-        std::complex<real_type> *out) {                                        \
-      prefix##_mpi_execute_dft(p, (complex_type *)in, (complex_type *)out);    \
-    }                                                                          \
-    static void execute_r2c(plan_type p, real_type *in, complex_type *out) {   \
-      prefix##_mpi_execute_dft_r2c(p, in, out);                                \
-    }                                                                          \
-    static void                                                                \
-    execute_c2r(plan_type p, std::complex<real_type> *in, real_type *out) {    \
-      prefix##_mpi_execute_dft_c2r(p, (complex_type *)in, out);                \
-    }                                                                          \
-    static void execute_c2r(plan_type p, complex_type *in, real_type *out) {   \
-      prefix##_mpi_execute_dft_c2r(p, in, out);                                \
-    }                                                                          \
-    static void                                                                \
-    execute_r2c(plan_type p, real_type *in, std::complex<real_type> *out) {    \
-      prefix##_mpi_execute_dft_r2c(p, in, (complex_type *)out);                \
-    }                                                                          \
-                                                                               \
-    static plan_type plan_dft_r2c_1d(                                          \
-        int n, real_type *in, complex_type *out, MPI_Comm, unsigned flags) {   \
-      return prefix##_plan_dft_r2c_1d(n, in, out, flags);                      \
-    }                                                                          \
-                                                                               \
-    static plan_type plan_dft_r2c_2d(                                          \
-        int Nx, int Ny, real_type *in, complex_type *out, MPI_Comm comm,       \
-        unsigned flags) {                                                      \
-      return prefix##_mpi_plan_dft_r2c_2d(Nx, Ny, in, out, comm, flags);       \
-    }                                                                          \
-                                                                               \
-    static plan_type plan_dft_c2r_1d(                                          \
-        int n, complex_type *in, real_type *out, MPI_Comm, unsigned flags) {   \
-      return prefix##_plan_dft_c2r_1d(n, in, out, flags);                      \
-    }                                                                          \
-    static plan_type plan_dft_c2r_2d(                                          \
-        int Nx, int Ny, complex_type *in, real_type *out, MPI_Comm comm,       \
-        unsigned flags) {                                                      \
-      return prefix##_mpi_plan_dft_c2r_2d(Nx, Ny, in, out, comm, flags);       \
-    }                                                                          \
-                                                                               \
-    static plan_type plan_dft_r2c_3d(                                          \
-        int Nx, int Ny, int Nz, real_type *in, complex_type *out,              \
-        MPI_Comm comm, unsigned flags) {                                       \
-      return prefix##_mpi_plan_dft_r2c_3d(Nx, Ny, Nz, in, out, comm, flags);   \
-    }                                                                          \
-    static plan_type plan_dft_c2r_3d(                                          \
-        int Nx, int Ny, int Nz, complex_type *in, real_type *out,              \
-        MPI_Comm comm, unsigned flags) {                                       \
-      return prefix##_mpi_plan_dft_c2r_3d(Nx, Ny, Nz, in, out, comm, flags);   \
-    }                                                                          \
-                                                                               \
-    static plan_type plan_dft_r2c(                                             \
-        int rank, const ptrdiff_t *n, real_type *in, complex_type *out,        \
-        MPI_Comm comm, unsigned flags) {                                       \
-      return prefix##_mpi_plan_dft_r2c(rank, n, in, out, comm, flags);         \
-    }                                                                          \
-    static plan_type plan_dft_c2r(                                             \
-        int rank, const ptrdiff_t *n, complex_type *in, real_type *out,        \
-        MPI_Comm comm, unsigned flags) {                                       \
-      return prefix##_mpi_plan_dft_c2r(rank, n, in, out, comm, flags);         \
-    }                                                                          \
-    static plan_type plan_dft_3d(                                              \
-        int Nx, int Ny, int Nz, complex_type *in, complex_type *out,           \
-        MPI_Comm comm, int sign, unsigned flags) {                             \
-      return prefix##_mpi_plan_dft_3d(Nx, Ny, Nz, in, out, comm, sign, flags); \
-    }                                                                          \
-    static plan_type plan_dft_2d(                                              \
-        int Nx, int Ny, complex_type *in, complex_type *out, MPI_Comm comm,    \
-        int sign, unsigned flags) {                                            \
-      return prefix##_mpi_plan_dft_2d(Nx, Ny, in, out, comm, sign, flags);     \
-    }                                                                          \
-    static plan_type plan_dft_1d(                                              \
-        int Nx, complex_type *in, complex_type *out, MPI_Comm comm, int sign,  \
-        unsigned flags) {                                                      \
-      return prefix##_plan_dft_1d(Nx, in, out, sign, flags);                   \
-    }                                                                          \
-    static void destroy_plan(plan_type plan) { prefix##_destroy_plan(plan); }  \
-  }
 
-  FFTW_MPI_CALLS_BASE(double, fftw);
+#define FFTW_MPI_CALLS_BASE(rtype, prefix) \
-  FFTW_MPI_CALLS_BASE(float, fftwf);
+  template<>				\
+class FFTW_MPI_Calls<rtype> {		\
+public: \
+  typedef rtype real_type; \
+  typedef prefix ## _complex complex_type; \
+  typedef prefix ## _plan plan_type; \
+  \
+  static complex_type *alloc_complex(size_t N) { return prefix ## _alloc_complex(N); } \
+  static real_type *alloc_real(size_t N) { return prefix ## _alloc_real(N); } \
+  static void free(void *p) { fftw_free(p); } \
+\
+  template<size_t Nd> \
+  static ptrdiff_t local_size(std::array<ptrdiff_t,Nd> const& N, MPI_Comm comm, \
+                              ptrdiff_t *local_n0, ptrdiff_t *local_0_start) { \
+    return prefix ## _mpi_local_size(Nd, N.data(), comm, local_n0, local_0_start); \
+  } \
+  static ptrdiff_t local_size_2d(ptrdiff_t N0, ptrdiff_t N1, MPI_Comm comm, \
+                                ptrdiff_t *local_n0, ptrdiff_t *local_0_start) { \
+    return prefix ## _mpi_local_size_2d(N0, N1, comm, local_n0, local_0_start); \
+  } \
+\
+  static ptrdiff_t local_size_3d(ptrdiff_t N0, ptrdiff_t N1, ptrdiff_t N2, MPI_Comm comm, \
+                                ptrdiff_t *local_n0, ptrdiff_t *local_0_start) { \
+    return prefix ## _mpi_local_size_3d(N0, N1, N2, comm, local_n0, local_0_start); \
+  } \
+\
+  static void execute(plan_type p) { prefix ## _execute(p); } \
+  static void execute_c2c(plan_type p, complex_type *in, complex_type *out) { prefix ## _mpi_execute_dft(p, in, out); } \
+  static void execute_c2c(plan_type p, std::complex<real_type> *in, std::complex<real_type> *out) { prefix ## _mpi_execute_dft(p, (complex_type*)in, (complex_type*)out); } \
+  static void execute_r2c(plan_type p, real_type *in, complex_type *out) { prefix ## _mpi_execute_dft_r2c(p, in, out); } \
+  static void execute_c2r(plan_type p, std::complex<real_type> *in, real_type *out) { prefix ## _mpi_execute_dft_c2r(p, (complex_type*)in, out); } \
+  static void execute_c2r(plan_type p, complex_type *in, real_type *out) { prefix ## _mpi_execute_dft_c2r(p, in, out); } \
+  static void execute_r2c(plan_type p, real_type *in, std::complex<real_type> *out) { prefix ## _mpi_execute_dft_r2c(p, in, (complex_type*)out); } \
+\
+  static plan_type plan_dft_r2c_2d(int Nx, int Ny,  \
+		       real_type *in, complex_type *out, \
+		       MPI_Comm comm, unsigned flags) \
+  { \
+    return prefix ## _mpi_plan_dft_r2c_2d(Nx, Ny, in, out, \
+				comm, flags); \
+  } \
+  static plan_type plan_dft_c2r_2d(int Nx, int Ny,  \
+                       complex_type *in, real_type *out, \
+                       MPI_Comm comm, unsigned flags) \
+  { \
+    return prefix ## _mpi_plan_dft_c2r_2d(Nx, Ny, in, out, \
+                                comm, flags); \
+  } \
+  static plan_type plan_dft_r2c_3d(int Nx, int Ny, int Nz, \
+                                   real_type *in, complex_type *out, \
+                                   MPI_Comm comm, unsigned flags) \
+  { \
+    return prefix ## _mpi_plan_dft_r2c_3d(Nx, Ny, Nz, in, out, comm, flags); \
+  } \
+  static plan_type plan_dft_c2r_3d(int Nx, int Ny, int Nz, \
+                                   complex_type *in, real_type *out, \
+                                   MPI_Comm comm, \
+                                   unsigned flags) \
+  { \
+    return prefix ## _mpi_plan_dft_c2r_3d(Nx, Ny, Nz, in, out, comm, flags); \
+  } \
+\
+  static plan_type plan_dft_r2c(int rank, const ptrdiff_t *n, real_type *in, \
+                                complex_type *out, MPI_Comm comm, unsigned flags) \
+  { \
+    return prefix ## _mpi_plan_dft_r2c(rank, n, in, out, comm, flags); \
+  } \
+  static plan_type plan_dft_c2r(int rank, const ptrdiff_t *n, complex_type *in, \
+                                real_type *out, MPI_Comm comm, unsigned flags) \
+  { \
+    return prefix ## _mpi_plan_dft_c2r(rank, n, in, out, comm, flags); \
+  } \
+  static plan_type plan_dft_3d(int Nx, int Ny, int Nz, complex_type *in, complex_type *out, MPI_Comm comm, int sign, unsigned flags) { \
+    return prefix ## _mpi_plan_dft_3d(Nx, Ny, Nz, in, out, comm, sign, flags); \
+  } \
+  static plan_type plan_dft_2d(int Nx, int Ny, complex_type *in, complex_type *out, MPI_Comm comm, int sign, unsigned flags) { \
+    return prefix ## _mpi_plan_dft_2d(Nx, Ny, in, out, comm, sign, flags); \
+  } \
+  static void destroy_plan(plan_type plan) { prefix ## _destroy_plan(plan); } \
+}
+
+
+FFTW_MPI_CALLS_BASE(double, fftw);
+FFTW_MPI_CALLS_BASE(float, fftwf);
 
 #undef FFTW_MPI_CALLS_BASE
 
-}; // namespace CosmoTool
+};
 
 #endif

src/h5_readFlash.cpp

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

src/hdf5_array.hpp

@@ -144,8 +144,8 @@ namespace CosmoTool {
                         bool useBases = false)
   {
       std::vector<hsize_t> memdims(data.shape(), data.shape() + data.num_dimensions());
-      H5::DataSpace dataspace(int(dimensions.size()), dimensions.data());
+      H5::DataSpace dataspace(dimensions.size(), dimensions.data());
-      H5::DataSpace memspace(int(memdims.size()), memdims.data());
+      H5::DataSpace memspace(memdims.size(), memdims.data());
 
       if (useBases) {
         std::vector<hsize_t> offsets(data.index_bases(), data.index_bases() + data.num_dimensions());
@@ -398,7 +398,7 @@ namespace CosmoTool {
             hdf5_weak_check_array(data, dimensions);
 
             std::vector<hsize_t> memdims(data.shape(), data.shape() + data.num_dimensions());
-            H5::DataSpace memspace(int(memdims.size()), memdims.data());
+            H5::DataSpace memspace(memdims.size(), memdims.data());
             
             std::vector<hsize_t> offsets(data.index_bases(), data.index_bases() + data.num_dimensions());            
             dataspace.selectHyperslab(H5S_SELECT_SET, memdims.data(), offsets.data());
@@ -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 = offsetof(STRUCT, BOOST_PP_TUPLE_ELEM(2, 1, element)); \
+       position = HOFFSET(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()); \
      } 
@@ -439,6 +439,7 @@ namespace CosmoTool { \
   \
       TNAME() : type(sizeof(STRUCT)) \
       {  \
+        long position; \
         BOOST_PP_SEQ_FOR_EACH(CTOOL_HDF5_INSERT_ELEMENT, STRUCT, ATTRIBUTES) \
       }  \
       \
@@ -470,6 +471,7 @@ namespace CosmoTool { \
   \
       TNAME() : type(sizeof(STRUCT)) \
       {  \
+        long position; \
         BOOST_PP_SEQ_FOR_EACH(CTOOL_HDF5_INSERT_ENUM_ELEMENT, STRUCT, ATTRIBUTES) \
       }  \
       \

src/interpolate.cpp

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

src/loadGadget.cpp

@@ -221,7 +221,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
       cerr << "Invalid format while reading header" << endl;
       delete data;
       delete f;
-      throw;
+      return 0;
     }
 
 
@@ -275,7 +275,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
         cerr << "Invalid format while reading positions" << endl;
         delete f;
         delete data;
-        throw;
+        return 0;
       }
     
   } else {
@@ -292,7 +292,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
 	  {
 	    delete f;
 	    delete data;
-            throw;
+	    return 0;
 	  }
       }
     
@@ -317,7 +317,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
 	cerr << "Invalid format while reading velocities" << endl;
 	delete f;
 	delete data;
-        throw;
+	return 0;
       }
     
     // THE VELOCITIES ARE IN PHYSICAL COORDINATES
@@ -367,7 +367,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id,
 	cerr << "Invalid unformatted access while reading ID" << endl;
 	delete f;
 	delete data;
-        throw;
+	return 0;
       }
   } else {
     f->skip(2*4);

src/mykdtree.hpp

@@ -1,5 +1,5 @@
 /*+
-This is CosmoTool (./src/mykdtree.hpp) -- Copyright (C) Guilhem Lavaux (2007-2022)
+This is CosmoTool (./src/mykdtree.hpp) -- Copyright (C) Guilhem Lavaux (2007-2014)
 
 guilhem.lavaux@gmail.com
 
@@ -99,8 +99,8 @@ namespace CosmoTool {
     typename KDDef<N,CType>::CoordType r, r2;
     KDCell<N, ValType,CType> **cells;
     typename KDDef<N,CType>::CoordType *distances;
-    uint64_t currentRank;
+    uint32_t currentRank;
-    uint64_t numCells;
+    uint32_t numCells;
   };
   
 
@@ -114,7 +114,7 @@ namespace CosmoTool {
 
     RecursionMultipleInfo(const typename KDDef<N,CType>::KDCoordinates& rx,
 			  KDCell<N,ValType,CType> **cells,
-			  uint64_t numCells)
+			  uint32_t numCells)
       : queue(cells, numCells, INFINITY),traversed(0)
     {
       std::copy(rx, rx+N, x);
@@ -153,20 +153,20 @@ namespace CosmoTool {
        std::copy(replicate, replicate+N, this->replicate);
      }
 
-    uint64_t getIntersection(const coords& x, CoordType r,
+    uint32_t getIntersection(const coords& x, CoordType r, 
 			     Cell **cells,
-			     uint64_t numCells);
+			     uint32_t numCells);
-    uint64_t getIntersection(const coords& x, CoordType r,
+    uint32_t getIntersection(const coords& x, CoordType r, 
 			     Cell **cells,
 			     CoordType *distances,
-			     uint64_t numCells);
+			     uint32_t numCells);
-    uint64_t countCells(const coords& x, CoordType r);
+    uint32_t countCells(const coords& x, CoordType r);
 
     Cell *getNearestNeighbour(const coords& x);
 
-    void getNearestNeighbours(const coords& x, uint64_t NumCells,
+    void getNearestNeighbours(const coords& x, uint32_t NumCells,
 			      Cell **cells);
-    void getNearestNeighbours(const coords& x, uint64_t NumCells,
+    void getNearestNeighbours(const coords& x, uint32_t NumCells,
 			      Cell **cells,
 			      CoordType *distances);
 

src/mykdtree.tcc

@@ -80,9 +80,9 @@ namespace CosmoTool {
   }
 
   template<int N, typename ValType, typename CType, typename CellSplitter>
-  uint64_t KDTree<N,ValType,CType,CellSplitter>::getIntersection(const coords& x, CoordType r,
+  uint32_t KDTree<N,ValType,CType,CellSplitter>::getIntersection(const coords& x, CoordType r, 
 						    KDTree<N,ValType,CType,CellSplitter>::Cell **cells,
-						    uint64_t numCells)
+						    uint32_t numCells)
   {
     RecursionInfoCells<N,ValType,CType> info;
 
@@ -112,10 +112,10 @@ namespace CosmoTool {
   }
 
   template<int N, typename ValType, typename CType, typename CellSplitter>
-  uint64_t KDTree<N,ValType,CType,CellSplitter>::getIntersection(const coords& x, CoordType r,
+  uint32_t KDTree<N,ValType,CType,CellSplitter>::getIntersection(const coords& x, CoordType r, 
 						    Cell **cells,
 						    CoordType *distances,
-						    uint64_t numCells)
+						    uint32_t numCells)
   {
     RecursionInfoCells<N,ValType> info;
 
@@ -144,7 +144,7 @@ namespace CosmoTool {
   }
 
   template<int N, typename ValType, typename CType, typename CellSplitter>
-  uint64_t KDTree<N,ValType,CType,CellSplitter>::countCells(const coords& x, CoordType r)
+  uint32_t KDTree<N,ValType,CType,CellSplitter>::countCells(const coords& x, CoordType r)
   {
     RecursionInfoCells<N,ValType> info;
 
@@ -501,7 +501,7 @@ namespace CosmoTool {
   }
 
  template<int N, typename ValType, typename CType, typename CellSplitter>
- void KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbours(const coords& x, uint64_t N2,
+ void KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbours(const coords& x, uint32_t N2,
 						    Cell **cells)
  {
    RecursionMultipleInfo<N,ValType> info(x, cells, N2);
@@ -527,7 +527,7 @@ namespace CosmoTool {
  }
 
  template<int N, typename ValType, typename CType, typename CellSplitter>
- void KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbours(const coords& x, uint64_t N2,
+ void KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbours(const coords& x, uint32_t N2,
 						    Cell **cells,
 						    CoordType *distances)
  {
@@ -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] >= 0 && node_on_disk.children_node[0] > lastNode) || node_on_disk.children_node[0] < -1 ||
+	    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)
+	    node_on_disk.children_node[1] > lastNode || node_on_disk.children_node[1] < -1)
 	  {
 	    delete[] nodes;
             std::cerr << "Invalid cell id or children node id invalid" << std::endl;

src/numpy_adaptors.hpp

@@ -1,30 +0,0 @@
-#ifndef __COSMOTOOL_NUMPY_ADAPTOR_HPP
-#define __COSMOTOOL_NUMPY_ADAPTOR_HPP
-
-namespace CosmoTool {
-
-  template<typename T, typename IT>
-  void parallel_ufunc_dd_d(char **args, IT* dimensions, IT* steps, void *func) {
-    IT i;
-    IT n = dimensions[0];
-    char *in = args[0], *in2 = args[1], *out = args[2];
-    IT in_step = steps[0], in2_step = steps[1], out_step = steps[2];
-
-    double tmp;
-    typedef double (*F_t)(double,double);
-
-    F_t f = (F_t)func;
-
-#pragma omp parallel for schedule(static)
-    for (i = 0; i < n; i++) {
-        T *out_t = (T *)(out + i * out_step);
-        T *in_t = (T *)(in + i * in_step);
-        T *in2_t = (T *)(in2 + i * in2_step);
-        *out_t = f(*in_t, *in2_t);
-    }
-  }
-
-
-}
-
-#endif

src/powerSpectrum.cpp

@@ -58,7 +58,6 @@ using namespace std;
 #define POWER_TEST 8
 
 #define POWER_SPECTRUM POWER_EFSTATHIOU
-#define SAMPLE_WIGGLES 9
 
 namespace Cosmology {
 
@@ -207,31 +206,6 @@ 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;

src/sphSmooth.hpp

@@ -1,5 +1,5 @@
 /*+
-This is CosmoTool (./src/sphSmooth.hpp) -- Copyright (C) Guilhem Lavaux (2007-2022)
+This is CosmoTool (./src/sphSmooth.hpp) -- Copyright (C) Guilhem Lavaux (2007-2014)
 
 guilhem.lavaux@gmail.com
 
@@ -39,26 +39,30 @@ 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 KDTree<Ndims,FullType> SPHTree;
-    typedef KDTreeNode<Ndims, FullType> SPHNode;
+    typedef KDTreeNode<Ndims,FullType> SPHNode;
-    typedef KDCell<Ndims, FullType> SPHCell;
+    typedef KDCell<Ndims,FullType>  SPHCell;
-    typedef typename KDTree<Ndims, FullType>::CoordType CoordType;
+    typedef typename KDTree<Ndims,FullType>::CoordType CoordType;
 
-    typedef ComputePrecision (*computeParticleValue)(const ValType &t);
+    typedef ComputePrecision (*computeParticleValue)(const ValType& t);
-    typedef void (*runParticleValue)(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;
@@ -66,58 +70,45 @@ namespace CosmoTool {
       ComputePrecision smoothRadius;
     };
     
+  
     SPHSmooth(SPHTree *tree, uint32_t Nsph);
     virtual ~SPHSmooth();   
 
-    void
+    void fetchNeighbours(const typename SPHTree::coords& c, SPHState *state = 0);
-    fetchNeighbours(const typename SPHTree::coords &c, SPHState *state = 0);
     
-    void fetchNeighbours(
+    void fetchNeighbours(const typename SPHTree::coords& c, uint32_t newNsph);
-        const typename SPHTree::coords &c, uint32_t newNsph,
+    void fetchNeighboursOnVolume(const typename SPHTree::coords& c, ComputePrecision radius);
-        SPHState *state = 0);
+    const typename SPHTree::coords& getCurrentCenter() const
-    void fetchNeighboursOnVolume(
+    {
-        const typename SPHTree::coords &c, ComputePrecision radius);
-    const typename SPHTree::coords &getCurrentCenter() const {
       return internal.currentCenter;
     }
 
-    /** This is the pure SPH smoothing function. It does not reweight by the
+    template<typename FuncT>
-     * value computed at each grid site.
+    ComputePrecision computeSmoothedValue(const typename SPHTree::coords& c,
-     */
+					  FuncT fun, SPHState *state = 0);    
-    template <typename FuncT>
-    ComputePrecision computeSmoothedValue(
-        const typename SPHTree::coords &c, FuncT fun, SPHState *state = 0);
 	
-    /** This is the weighted SPH smoothing function. It does reweight by the
+	template<typename FuncT>
-     * value computed at each grid site. This ensures the total sum of the interpolated
+    ComputePrecision computeInterpolatedValue(const typename SPHTree::coords& c,
-     * quantity is preserved by interpolating to the target mesh.
+					      FuncT fun, SPHState *state = 0);
-     */
-    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
+    ComputePrecision getMaxDistance(const typename SPHTree::coords& c, 
-     * array. FuncT is expected to have the following prototype:
+				    SPHNode *node) const;
-     *    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
+    ComputePrecision getSmoothingLen() const
-    getMaxDistance(const typename SPHTree::coords &c, SPHNode *node) const;
+    {
-
+      return internal.smoothRadius;
-    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, SPHState *state);
+    void addGridSite(const typename SPHTree::coords& c);
-    void addGridSite(const typename SPHTree::coords &c);
 
     bool hasNeighbours() const;
 
@@ -142,21 +133,21 @@ namespace CosmoTool {
     uint32_t maxNgb;
     SPHTree *tree;
     
-    template <typename FuncT>
+    template<typename FuncT>
-    ComputePrecision computeWValue(
+    ComputePrecision computeWValue(const typename SPHTree::coords & c,
-        const typename SPHTree::coords &c, SPHCell &cell, CoordType d,
+				   SPHCell& cell,
-        FuncT fun, SPHState *state);
+				   CoordType d,
+				   FuncT fun, SPHState *state);
     
-    template <typename FuncT>
+    template<typename FuncT>
-    void runUnrollNode(SPHNode *node, FuncT fun);
+    void runUnrollNode(SPHNode *node,
+		       FuncT fun);
   };
 
-  template <class ValType1, class ValType2, int Ndims>
+  template<class ValType1, class ValType2, int Ndims>
-  bool operator<(
+  bool operator<(const SPHSmooth<ValType1, Ndims>& s1, const SPHSmooth<ValType2, Ndims>& s2);
-      const SPHSmooth<ValType1, Ndims> &s1,
-      const SPHSmooth<ValType2, Ndims> &s2);
 
-}; // namespace CosmoTool
+};
 
 #include "sphSmooth.tcc"
 

src/sphSmooth.tcc

@@ -3,248 +3,225 @@
 
 namespace CosmoTool {
 
-  template <typename ValType, int Ndims>
+template<typename ValType, int Ndims>
-  SPHSmooth<ValType, Ndims>::SPHSmooth(SPHTree *tree, uint32_t Nsph) {
+SPHSmooth<ValType,Ndims>::SPHSmooth(SPHTree *tree, uint32_t Nsph)
-    this->Nsph = Nsph;
+{
-    this->tree = tree;
+  this->Nsph = Nsph;
-    internal.currentNgb = 0;
+  this->tree = tree;
+  internal.currentNgb = 0;
 
-    this->maxNgb = Nsph;
+  this->maxNgb = Nsph;
-    internal.ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[maxNgb]);
+  internal.ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[maxNgb]);  
-    internal.distances = boost::shared_ptr<CoordType[]>(new CoordType[maxNgb]);
+  internal.distances = boost::shared_ptr<CoordType[]>(new CoordType[maxNgb]);
-  }
+}
 
-  template <typename ValType, int Ndims>
+template<typename ValType, int Ndims>
-  SPHSmooth<ValType, Ndims>::~SPHSmooth() {}
+SPHSmooth<ValType,Ndims>::~SPHSmooth()
+{
+}
 
-  template <typename ValType, int Ndims>
+template<typename ValType, int Ndims>
-  template <typename FuncT>
+template<typename FuncT>
-  ComputePrecision SPHSmooth<ValType, Ndims>::computeWValue(
+ComputePrecision SPHSmooth<ValType,Ndims>::computeWValue(const typename SPHTree::coords& c,
-      const typename SPHTree::coords &c, SPHCell &cell, CoordType d, FuncT fun,
+							 SPHCell& cell,
-      SPHState *state) {
+							 CoordType d,
-    CoordType weight;
+							 FuncT fun, SPHState *state)
+{
+  CoordType weight;
 
-    d /= state->smoothRadius;
+  d /= state->smoothRadius;
-    weight = getKernel(d);
+  weight = getKernel(d);
 
-    if (cell.val.weight != 0)
+  if (cell.val.weight != 0)
-      return weight * fun(cell.val.pValue) / cell.val.weight;
+    return weight * fun(cell.val.pValue) / cell.val.weight;
-    else
+  else
-      return 0;
+    return 0;
-  }
+}
 
-  template <typename ValType, int Ndims>
+template<typename ValType, int Ndims>
-  void SPHSmooth<ValType, Ndims>::fetchNeighbours(
+void
-      const typename SPHTree::coords &c, uint32_t newNngb, SPHState *state) {
+SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c, uint32_t newNngb)
-    ComputePrecision d2, max_dist = 0;
+{
-    uint32_t requested = newNngb;
+  ComputePrecision d2, max_dist = 0;
+  uint32_t requested = newNngb;
 
-    if (state != 0) {
+  if (requested > maxNgb)
-      state->distances = boost::shared_ptr<CoordType[]>(new CoordType[newNngb]);
+    {
-      state->ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[newNngb]);
+      maxNgb = requested;
-    } else {
+      internal.ngb = boost::shared_ptr<P_SPHCell[]>(new P_SPHCell[maxNgb]);
-      state = &internal;
+      internal.distances = boost::shared_ptr<CoordType[]>(new CoordType[maxNgb]);
-      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(state->currentCenter, c, sizeof(c));
+  memcpy(internal.currentCenter, c, sizeof(c));
-    tree->getNearestNeighbours(
+  tree->getNearestNeighbours(c, requested, (SPHCell **)internal.ngb.get(), (CoordType*)internal.distances.get());
-        c, requested, (SPHCell **)state->ngb.get(),
-        (CoordType *)state->distances.get());
 
-    state->currentNgb = 0;
+  internal.currentNgb = 0;
-    for (uint32_t i = 0; i < requested && (state->ngb)[i] != 0;
+  for (uint32_t i = 0; i < requested && (internal.ngb)[i] != 0; i++,internal.currentNgb++)
-         i++, state->currentNgb++) {
+    {
-      state->distances[i] = sqrt(state->distances[i]);
+      internal.distances[i] = sqrt(internal.distances[i]);
-      d2 = state->distances[i];
+      d2 = internal.distances[i];
       if (d2 > max_dist)
         max_dist = d2;
     }
 
-    state->smoothRadius = max_dist / 2;
+  internal.smoothRadius = max_dist / 2;  
-  }
+}
 
-  template <typename ValType, int Ndims>
+template<typename ValType, int Ndims>
-  void SPHSmooth<ValType, Ndims>::fetchNeighbours(
+void SPHSmooth<ValType,Ndims>::fetchNeighbours(const typename SPHTree::coords& c, SPHState *state)
-      const typename SPHTree::coords &c, SPHState *state) {
+{
-    ComputePrecision d2, max_dist = 0;
+  ComputePrecision d2, max_dist = 0;
-    uint32_t requested = Nsph;
+  uint32_t requested = Nsph;
 
-    if (state != 0) {
+  if (state != 0) {
-      state->distances = boost::shared_ptr<CoordType[]>(new CoordType[Nsph]);
+    state->distances = boost::shared_ptr<CoordType[]>(new CoordType[Nsph]);
-      state->ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[Nsph]);
+    state->ngb = boost::shared_ptr<SPHCell *[]>(new SPHCell *[Nsph]);
-    } else
+  } else
-      state = &internal;
+    state = &internal;
     
-    memcpy(state->currentCenter, c, sizeof(c));
+  memcpy(state->currentCenter, c, sizeof(c));
   
-    tree->getNearestNeighbours(
+  tree->getNearestNeighbours(c, requested, state->ngb.get(), state->distances.get());
-        c, requested, state->ngb.get(), state->distances.get());
 
-    state->currentNgb = 0;
+  state->currentNgb = 0;
-    for (uint32_t i = 0; i < requested && state->ngb[i] != 0;
+  for (uint32_t i = 0; i < requested && state->ngb[i] != 0; i++,state->currentNgb++)
-         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;
-  }
 
-  template <typename ValType, int Ndims>
+  state->smoothRadius = max_dist / 2;  
-  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++) {
+template<typename ValType, int Ndims>
-      outputValue +=
+void
-          computeWValue(c, *state->ngb[i], state->distances[i], fun, state);
+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);
     }
 
-    return outputValue / r3;
+  return outputValue / r3;
-  }
+}
 
-  template <typename ValType>
+template<typename ValType>
-  ComputePrecision interpolateOne(const ValType &t) {
+ComputePrecision interpolateOne(const ValType& t)
-    return 1.0;
+{
-  }
+  return 1.0;
+}
 
-  template <typename ValType, int Ndims>
+// WARNING ! Cell's weight must be 1 !!!
-  template <typename FuncT>
+template<typename ValType, int Ndims>
-  void SPHSmooth<ValType, Ndims>::computeAdjointGradientSmoothedValue(
+template<typename FuncT>
-      const typename SPHTree::coords &c, ComputePrecision ag_value, FuncT fun,
+ComputePrecision SPHSmooth<ValType,Ndims>::computeInterpolatedValue(const typename SPHTree::coords& c,
-      SPHState *state) {
+								    FuncT fun, SPHState *state)
-    if (state == 0)
+{
-      state = &internal;
+  if (state == 0)
+    state = &internal;
     
-    ComputePrecision outputValue = 0;
+  ComputePrecision outputValue = 0;
-    ComputePrecision max_dist = 0;
+  ComputePrecision max_dist = 0;
-    ComputePrecision weight = 0;
+  ComputePrecision weight = 0;
 
-    for (uint32_t i = 0; i < state->currentNgb; i++) {
+  for (uint32_t i = 0; i < state->currentNgb; i++)
-      weight +=
+    {
-          computeWValue(c, *state->ngb[i], state->distances[i], interpolateOne);
-    }
-
-    for (uint32_t i = 0; i < state->currentNgb; i++) {
-      auto &cell = *state->ngb[i];
-      double partial_ag =
-          computeWValue(
-              c, cell, state->distances[i],
-              [ag_value](ComputePrecision) { return ag_value; }) /
-          weight;
-      fun(cell.val.pValue, ag_value);
-    }
-  }
-
-  // WARNING ! Cell's weight must be 1 !!!
-  template <typename ValType, int Ndims>
-  template <typename FuncT>
-  ComputePrecision SPHSmooth<ValType, Ndims>::computeInterpolatedValue(
-      const typename SPHTree::coords &c, FuncT fun, SPHState *state) {
-    if (state == 0)
-      state = &internal;
-
-    ComputePrecision outputValue = 0;
-    ComputePrecision max_dist = 0;
-    ComputePrecision weight = 0;
-
-    for (uint32_t i = 0; i < state->currentNgb; i++) {
       outputValue += computeWValue(c, *state->ngb[i], state->distances[i], fun);
-      weight +=
+      weight += computeWValue(c, *state->ngb[i], state->distances[i], interpolateOne);
-          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 ValType, int Ndims>
-  template <typename FuncT>
+template<typename FuncT>
-  void
+void SPHSmooth<ValType,Ndims>::runForEachNeighbour(FuncT fun, SPHState *state)
-  SPHSmooth<ValType, Ndims>::runForEachNeighbour(FuncT fun, SPHState *state) {
+{
-    if (state == 0)
+  if (state == 0)
-      state = &internal;
+    state = &internal;
     
-    for (uint32_t i = 0; i < state->currentNgb; i++) {
+  for (uint32_t i = 0; i < state->currentNgb; i++)
+    {
       fun(state->ngb[i]);
     }
-  }
+}
 
-  template <typename ValType, int Ndims>
-  void SPHSmooth<ValType, Ndims>::addGridSite(
-      const typename SPHTree::coords &c, SPHState *state) {
-    ComputePrecision outputValue = 0;
-    ComputePrecision max_dist = 0;
-    ComputePrecision r3 = cube(state->smoothRadius);
 
-    for (uint32_t i = 0; i < state->currentNgb; i++) {
+template<typename ValType, int Ndims>
-      ComputePrecision d = state->distances[i];
+void SPHSmooth<ValType,Ndims>::addGridSite(const typename SPHTree::coords& c)
-      SPHCell &cell = *(state->ngb[i]);
+{
-      double kernel_value = getKernel(d / state->smoothRadius) / r3;
+  ComputePrecision outputValue = 0;
-#pragma omp atomic update
+  ComputePrecision max_dist = 0;
-      cell.val.weight += kernel_value;
+  
+  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;
     }
-  }
+}
 
-  template <typename ValType, int Ndims>
+template<typename ValType, int Ndims>
-  void
+ComputePrecision 
-  SPHSmooth<ValType, Ndims>::addGridSite(const typename SPHTree::coords &c) {
+SPHSmooth<ValType,Ndims>::getKernel(ComputePrecision x) const
-    addGridSite(c, &internal);
+{
-  }
+  // WARNING !!! This is an unnormalized version of the kernel.
-
+  if (x < 1)
-  template <typename ValType, int Ndims>
+    return 1 - 1.5 * x * x + 0.75 * x * x * x;
-  ComputePrecision
+  else if (x < 2)
-  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>
+template<typename ValType, int Ndims>
-  bool SPHSmooth<ValType, Ndims>::hasNeighbours() const {
+bool SPHSmooth<ValType,Ndims>::hasNeighbours() const
-    return (internal.currentNgb != 0);
+{
-  }
+  return (internal.currentNgb != 0);
+}
 
-  template <class ValType1, class ValType2, int Ndims>
+template<class ValType1, class ValType2, int Ndims>
-  bool operator<(
+bool operator<(const SPHSmooth<ValType1, Ndims>& s1, const SPHSmooth<ValType2, Ndims>& s2)
-      const SPHSmooth<ValType1, Ndims> &s1,
+{
-      const SPHSmooth<ValType2, Ndims> &s2) {
+  return (s1.getSmoothingLen() < s2.getSmoothingLen());
-    return (s1.getSmoothingLen() < s2.getSmoothingLen());
+}
-  }
+
-}; // namespace CosmoTool
+};

src/yorick.hpp

@@ -222,7 +222,8 @@ namespace CosmoTool
 
   template<typename T>
   void loadArray(const std::string& fname,
-		 T*& array, uint32_t *& dimList, uint32_t& rank);
+		 T*& array, uint32_t *& dimList, uint32_t& rank)
+    throw (NoSuchFileException);
 
   ProgressiveDoubleOutput saveDoubleArrayProgressive(const char *fname, uint32_t *dimList, uint32_t rank);  
 };

src/yorick_nc3.cpp

@@ -260,6 +260,7 @@ 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);
     

src/yorick_nc4.cpp

@@ -224,6 +224,7 @@ 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::read);