guilhem_lavaux/cosmotool
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Compare commits

base: guilhem_lavaux:master
Branches Tags
guilhem_lavaux:master
guilhem_lavaux:eleni/BAO
guilhem_lavaux:andrija/mpi
guilhem_lavaux:natalia/qlpt_feature
guilhem_lavaux:v1.3.6
guilhem_lavaux:v1.3.5
guilhem_lavaux:v1.3.4
guilhem_lavaux:v1.3.3
guilhem_lavaux:v1.3.2
guilhem_lavaux:v1.Z
guilhem_lavaux:v1.3.1
guilhem_lavaux:v1.3.0
guilhem_lavaux:v1.2.3
guilhem_lavaux:v1.2.2
guilhem_lavaux:v1.2.1
guilhem_lavaux:v1.2.0
guilhem_lavaux:v1.1.1
guilhem_lavaux:v1.1.0
guilhem_lavaux:v1.0.0a8
guilhem_lavaux:v1.0.0a7
guilhem_lavaux:v1.0.0a6
..
compare: guilhem_lavaux:natalia/qlpt_feature
Branches Tags
guilhem_lavaux:master
guilhem_lavaux:eleni/BAO
guilhem_lavaux:andrija/mpi
guilhem_lavaux:natalia/qlpt_feature
guilhem_lavaux:v1.3.6
guilhem_lavaux:v1.3.5
guilhem_lavaux:v1.3.4
guilhem_lavaux:v1.3.3
guilhem_lavaux:v1.3.2
guilhem_lavaux:v1.Z
guilhem_lavaux:v1.3.1
guilhem_lavaux:v1.3.0
guilhem_lavaux:v1.2.3
guilhem_lavaux:v1.2.2
guilhem_lavaux:v1.2.1
guilhem_lavaux:v1.2.0
guilhem_lavaux:v1.1.1
guilhem_lavaux:v1.1.0
guilhem_lavaux:v1.0.0a8
guilhem_lavaux:v1.0.0a7
guilhem_lavaux:v1.0.0a6

No commits in common. "master" and "natalia/qlpt_feature" have entirely different histories.
master ... natalia/ql

42 changed files with 996 additions and 5092 deletions
Show Stats Expand all files Collapse all files

5
.gitignore vendored
View File

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

18
CMakeLists.txt
View File

@ -42,17 +42,6 @@ ENDIF(BUILD_PYTHON)
MESSAGE(STATUS "Using the target ${CosmoTool_local} to build python module") 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) include(${CMAKE_SOURCE_DIR}/external/external_build.cmake)
IF(YORICK_SUPPORT) IF(YORICK_SUPPORT)
@ -66,6 +55,9 @@ IF(YORICK_SUPPORT)
ENDIF(YORICK_SUPPORT) ENDIF(YORICK_SUPPORT)
find_program(CYTHON NAMES cython3 cython) 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(NETCDF_FIND_REQUIRED ${YORICK_SUPPORT})
set(GSL_FIND_REQUIRED TRUE) 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_PACKAGE_VENDOR "Guilhem Lavaux")
SET(CPACK_RESOURCE_FILE_LICENSE "${CMAKE_CURRENT_SOURCE_DIR}/LICENCE_CeCILL_V2") SET(CPACK_RESOURCE_FILE_LICENSE "${CMAKE_CURRENT_SOURCE_DIR}/LICENCE_CeCILL_V2")
SET(CPACK_PACKAGE_VERSION_MAJOR "1") 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_PACKAGE_INSTALL_DIRECTORY "CosmoToolbox-${CPACK_PACKAGE_VERSION_MAJOR}.${CPACK_PACKAGE_VERSION_MINOR}")
SET(CPACK_STRIP_FILES "lib/libCosmoTool.so") SET(CPACK_STRIP_FILES "lib/libCosmoTool.so")
SET(CPACK_SOURCE_IGNORE_FILES SET(CPACK_SOURCE_IGNORE_FILES

9
MANIFEST.in
View File

@ -1,5 +1,4 @@
include .gitignore include .gitignore
include pyproject.toml
include CMakeLists.txt include CMakeLists.txt
include FindNumPy.cmake include FindNumPy.cmake
include FindPyLibs.cmake include FindPyLibs.cmake
@ -17,17 +16,14 @@ include doc/source/cpplibrary.rst
include doc/source/index.rst include doc/source/index.rst
include doc/source/intro.rst include doc/source/intro.rst
include doc/source/pythonmodule.rst include doc/source/pythonmodule.rst
include external/config.guess
include external/config.sub
include external/external_build.cmake 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/omptl-20120422.tar.bz2
include external/patch-omptl include external/patch-omptl
include python/CMakeLists.txt include python/CMakeLists.txt
include python/_cosmo_bispectrum.cpp include python/_cosmo_bispectrum.cpp
include python/_cosmo_cic.pyx include python/_cosmo_cic.pyx
include python/_cosmo_power.pyx include python/_cosmo_power.pyx
include python/_cosmomath.pyx
include python/_cosmotool.pyx include python/_cosmotool.pyx
include python/_fast_interp.pyx include python/_fast_interp.pyx
include python/_project.pyx include python/_project.pyx
@ -90,14 +86,12 @@ include sample/testSmooth.cpp
include sample/test_cosmopower.cpp include sample/test_cosmopower.cpp
include sample/test_fft_calls.cpp include sample/test_fft_calls.cpp
include sample/test_healpix_calls.cpp include sample/test_healpix_calls.cpp
include sample/test_special.cpp
include sample/testkd.cpp include sample/testkd.cpp
include sample/testkd2.cpp include sample/testkd2.cpp
include sample/testkd3.cpp include sample/testkd3.cpp
include setup.py include setup.py
include src/CMakeLists.txt include src/CMakeLists.txt
include src/algo.hpp include src/algo.hpp
include src/numpy_adaptors.hpp
include src/bqueue.hpp include src/bqueue.hpp
include src/bqueue.tcc include src/bqueue.tcc
include src/bsp_simple.hpp include src/bsp_simple.hpp
@ -173,4 +167,3 @@ include src/tf_fit.hpp
include src/yorick.hpp include src/yorick.hpp
include src/yorick_nc3.cpp include src/yorick_nc3.cpp
include src/yorick_nc4.cpp include src/yorick_nc4.cpp
include src/special_math.hpp

22
builder/build-wheels.sh
View File

@ -8,32 +8,28 @@ export CC CXX
# Install a system package required by our library # Install a system package required by our library
#yum install -y atlas-devel #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 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 # Compile wheels
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 -r /io/dev-requirements.txt # "${PYBIN}/pip" install -r /io/dev-requirements.txt
"${PYBIN}/pip" install setuptools wheel Cython
"${PYBIN}/pip" install -r /io/requirements.txt "${PYBIN}/pip" install -r /io/requirements.txt
"${PYBIN}/pip" wheel -vvv /io/ -w /io/wheelhouse/ "${PYBIN}/pip" wheel -vvv /io/ -w wheelhouse/
done done
# Bundle external shared libraries into the wheels # 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 done
# Install packages and test # 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

2
builder/start.sh
View File

@ -9,4 +9,4 @@ if ! [ -e ${d}/setup.py ] ; then
exit 1 exit 1
fi 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

3
conda/build.sh
View File

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

7
conda/conda_build_config.yaml
View File

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

32
conda/meta.yaml
View File

@ -1,50 +1,38 @@
package: package:
name: cosmotool name: cosmotool
version: "1.3.0" version: "1.0.0a7"
source: source:
git_rev: a86c9a8 git_rev: 7fce73e
git_url: https://bitbucket.org/glavaux/cosmotool git_url: https://bitbucket.org/glavaux/cosmotool
requirements: requirements:
build: 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 - python
- pip - setuptools
- numpy
- pkgconfig
- numexpr
- cython - cython
- healpy - healpy
- numexpr
- cffi - cffi
- pyfftw - pyfftw
- gsl {{ gsl }} - gsl
- h5py - h5py
run: run:
- numpy
- python - python
- {{ pin_compatible('numpy') }}
- healpy - healpy
- numexpr - numexpr
- cffi - cffi
- pyfftw - pyfftw
- h5py - h5py
- {{ pin_compatible('gsl') }}
- llvm-openmp
test: test:
imports: imports:
- cosmotool - cosmotool
requires:
- pip
about: about:
home: https://bitbucket.org/glavaux/cosmotool home: https://bitbucket.org/glavaux/cosmotool

1754
external/config.guess vendored
View File

File diff suppressed because it is too large Load Diff

1890
external/config.sub vendored
View File

File diff suppressed because it is too large Load Diff

75
external/external_build.cmake vendored
View File

@ -2,19 +2,16 @@ include(FindOpenMP)
OPTION(ENABLE_OPENMP "Set to Yes if Healpix and/or you need openMP" OFF) 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(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(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(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) mark_as_advanced(FFTW_URL EIGEN_URL HDF5_URL NETCDF_URL BOOST_URL GSL_URL)
file(MAKE_DIRECTORY ${SOURCE_PREFIX}/downloads)
SET(all_deps) SET(all_deps)
MACRO(CHECK_CHANGE_STATE VAR) MACRO(CHECK_CHANGE_STATE VAR)
@ -39,13 +36,12 @@ CHECK_CHANGE_STATE(INTERNAL_DLIB DLIB_INCLUDE_DIR DLIB_LIBRARIES)
IF(ENABLE_OPENMP) IF(ENABLE_OPENMP)
IF (NOT OPENMP_FOUND) IF (NOT OPENMP_FOUND)
MESSAGE(NOTICE "No known compiler option for enabling OpenMP") MESSAGE(ERROR "No known compiler option for enabling OpenMP")
ELSE() ENDIF(NOT OPENMP_FOUND)
SET(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") SET(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
SET(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_C_FLAGS}") SET(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_C_FLAGS}")
ENDIF(NOT OPENMP_FOUND)
ENDIF(ENABLE_OPENMP) ENDIF(ENABLE_OPENMP)
@ -59,25 +55,17 @@ if (ENABLE_SHARP)
SET(DEP_BUILD ${BUILD_PREFIX}/sharp-prefix/src/sharp/auto) SET(DEP_BUILD ${BUILD_PREFIX}/sharp-prefix/src/sharp/auto)
IF(NOT ENABLE_OPENMP) IF(NOT ENABLE_OPENMP)
SET(SHARP_OPENMP --disable-openmp) SET(SHARP_OPENMP --disable-openmp)
ELSE()
SET(SHARP_OPENMP)
ENDIF() ENDIF()
SET(CUTILS_LIBRARY ${DEP_BUILD}/lib/libc_utils.a) SET(CUTILS_LIBRARY ${DEP_BUILD}/lib/libc_utils.a)
SET(FFTPACK_LIBRARY ${DEP_BUILD}/lib/libfftpack.a) SET(FFTPACK_LIBRARY ${DEP_BUILD}/lib/libfftpack.a)
SET(SHARP_LIBRARY ${DEP_BUILD}/lib/libsharp.a) SET(SHARP_LIBRARY ${DEP_BUILD}/lib/libsharp.a)
SET(SHARP_LIBRARIES ${SHARP_LIBRARY} ${FFTPACK_LIBRARY} ${CUTILS_LIBRARY}) SET(SHARP_LIBRARIES ${SHARP_LIBRARY} ${FFTPACK_LIBRARY} ${CUTILS_LIBRARY})
SET(SHARP_INCLUDE_PATH ${DEP_BUILD}/include) SET(SHARP_INCLUDE_PATH ${DEP_BUILD}/include)
message(STATUS "Flags: ${CMAKE_C_FLAGS}")
ExternalProject_Add(sharp 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 PREFIX ${BUILD_PREFIX}/sharp-prefix
DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
BUILD_IN_SOURCE 1 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" INSTALL_COMMAND echo "No install"
BUILD_BYPRODUCTS ${SHARP_LIBRARIES} BUILD_BYPRODUCTS ${SHARP_LIBRARIES}
) )
@ -94,23 +82,27 @@ if (INTERNAL_HDF5)
ExternalProject_Add(hdf5 ExternalProject_Add(hdf5
PREFIX ${BUILD_PREFIX}/hdf5-prefix PREFIX ${BUILD_PREFIX}/hdf5-prefix
URL ${HDF5_URL} 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(cosmotool_DEPS ${cosmotool_DEPS} hdf5)
SET(hdf5_built hdf5) SET(hdf5_built hdf5)
SET(ENV{HDF5_ROOT} ${HDF5_BIN_DIR}) SET(ENV{HDF5_ROOT} ${HDF5_BIN_DIR})
SET(HDF5_ROOTDIR ${HDF5_BIN_DIR}) SET(HDF5_ROOTDIR ${HDF5_BIN_DIR})
SET(CONFIGURE_LDFLAGS "${CONFIGURE_LDFLAGS} -L${HDF5_BIN_DIR}/lib") 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) 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) mark_as_advanced(HDF5_LIBRARIES HDF5_CXX_LIBRARIES HDF5_INCLUDE_DIRS)
@ -168,28 +160,24 @@ if (INTERNAL_NETCDF)
SET(NETCDF_CONFIG_COMMAND ${NETCDF_SOURCE_DIR}/configure SET(NETCDF_CONFIG_COMMAND ${NETCDF_SOURCE_DIR}/configure
--prefix=${NETCDF_BIN_DIR} --libdir=${NETCDF_BIN_DIR}/lib --prefix=${NETCDF_BIN_DIR} --libdir=${NETCDF_BIN_DIR}/lib
--enable-netcdf-4 --with-pic --disable-shared --disable-dap --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} --disable-examples ${EXTRA_NC_FLAGS} CC=${CMAKE_C_COMPILER}
CXX=${CMAKE_CXX_COMPILER}) CXX=${CMAKE_CXX_COMPILER})
list(INSERT CMAKE_PREFIX_PATH 0 ${EXT_INSTALL}) list(INSERT CMAKE_PREFIX_PATH 0 ${EXT_INSTALL})
string(REPLACE ";" "|" CMAKE_PREFIX_PATH_ALT_SEP "${CMAKE_PREFIX_PATH}") string(REPLACE ";" "|" CMAKE_PREFIX_PATH_ALT_SEP "${CMAKE_PREFIX_PATH}")
ExternalProject_Add(netcdf ExternalProject_Add(netcdf
DEPENDS ${hdf5_built} DEPENDS ${hdf5_built}
URL_HASH MD5=f48ee01534365006934f0c63d4055ea0
PREFIX ${BUILD_PREFIX}/netcdf-prefix PREFIX ${BUILD_PREFIX}/netcdf-prefix
DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
URL ${NETCDF_URL} URL ${NETCDF_URL}
LIST_SEPARATOR | LIST_SEPARATOR |
CMAKE_ARGS CMAKE_ARGS
-DCMAKE_PREFIX_PATH=${CMAKE_PREFIX_PATH_ALT_SEP} -DCMAKE_PREFIX_PATH=${CMAKE_PREFIX_PATH_ALT_SEP}
-DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER}
-DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER} -DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER}
-DNC_EXTRA_DEPS=${RT_DEP}
-DBUILD_SHARED_LIBS=OFF -DBUILD_SHARED_LIBS=OFF
-DBUILD_TESTING=OFF -DBUILD_TESTING=OFF
-DCMAKE_BUILD_TYPE=Release -DCMAKE_BUILD_TYPE=Release
-DENABLE_NETCDF4=ON -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} -DCMAKE_INSTALL_PREFIX=${NETCDF_BIN_DIR}
@ -203,17 +191,12 @@ if (INTERNAL_NETCDF)
ExternalProject_Add(netcdf-c++ ExternalProject_Add(netcdf-c++
DEPENDS ${hdf5_built} netcdf DEPENDS ${hdf5_built} netcdf
PREFIX ${BUILD_PREFIX}/netcdf-c++-prefix PREFIX ${BUILD_PREFIX}/netcdf-c++-prefix
DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
URL ${NETCDFCXX_URL} URL ${NETCDFCXX_URL}
CMAKE_ARGS CMAKE_ARGS
-DCMAKE_PREFIX_PATH=${CMAKE_PREFIX_PATH_ALT_SEP}
-DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER}
-DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER} -DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER}
-DBUILD_SHARED_LIBS=OFF -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 -DBUILD_TESTING=OFF
-DCMAKE_BUILD_TYPE=Release -DCMAKE_BUILD_TYPE=Release
-DCMAKE_INSTALL_PREFIX=${NETCDF_BIN_DIR} -DCMAKE_INSTALL_PREFIX=${NETCDF_BIN_DIR}
@ -251,8 +234,6 @@ if (INTERNAL_BOOST)
ExternalProject_Add(boost ExternalProject_Add(boost
URL ${BOOST_URL} URL ${BOOST_URL}
PREFIX ${BUILD_PREFIX}/boost-prefix PREFIX ${BUILD_PREFIX}/boost-prefix
DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
URL_HASH MD5=f7050f554a65f6a42ece221eaeec1660
CONFIGURE_COMMAND CONFIGURE_COMMAND
${BOOST_SOURCE_DIR}/bootstrap.sh --prefix=${CMAKE_BINARY_DIR}/ext_build/boost ${BOOST_SOURCE_DIR}/bootstrap.sh --prefix=${CMAKE_BINARY_DIR}/ext_build/boost
BUILD_IN_SOURCE 1 BUILD_IN_SOURCE 1
@ -295,7 +276,6 @@ IF(INTERNAL_GSL)
ExternalProject_Add(gsl ExternalProject_Add(gsl
URL ${GSL_URL} URL ${GSL_URL}
PREFIX ${BUILD_PREFIX}/gsl-prefix PREFIX ${BUILD_PREFIX}/gsl-prefix
DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
CONFIGURE_COMMAND ${GSL_SOURCE_DIR}/configure CONFIGURE_COMMAND ${GSL_SOURCE_DIR}/configure
--prefix=${EXT_INSTALL} --disable-shared --prefix=${EXT_INSTALL} --disable-shared
--with-pic --with-pic
@ -348,7 +328,6 @@ IF(INTERNAL_FFTW)
ExternalProject_Add(fftw ExternalProject_Add(fftw
URL ${FFTW_URL} URL ${FFTW_URL}
PREFIX ${BUILD_PREFIX}/fftw-prefix PREFIX ${BUILD_PREFIX}/fftw-prefix
DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
CONFIGURE_COMMAND CONFIGURE_COMMAND
${FFTW_SOURCE}/configure ${FFTW_SOURCE}/configure
--prefix=${EXT_INSTALL} --prefix=${EXT_INSTALL}
@ -378,7 +357,6 @@ IF (INTERNAL_EIGEN)
ExternalProject_Add(eigen ExternalProject_Add(eigen
URL ${EIGEN_URL} URL ${EIGEN_URL}
URL_HASH MD5=b9e98a200d2455f06db9c661c5610496 URL_HASH MD5=b9e98a200d2455f06db9c661c5610496
DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
PREFIX ${BUILD_PREFIX}/eigen-prefix PREFIX ${BUILD_PREFIX}/eigen-prefix
CMAKE_ARGS -DCMAKE_INSTALL_PREFIX=${EXT_INSTALL} CMAKE_ARGS -DCMAKE_INSTALL_PREFIX=${EXT_INSTALL}
-DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} -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) SET(OMPTL_BUILD_DIR ${BUILD_PREFIX}/omptl-prefix/src/omptl)
ExternalProject_Add(omptl ExternalProject_Add(omptl
PREFIX ${BUILD_PREFIX}/omptl-prefix PREFIX ${BUILD_PREFIX}/omptl-prefix
DOWNLOAD_DIR ${SOURCE_PREFIX}/downloads
URL ${CMAKE_SOURCE_DIR}/external/omptl-20120422.tar.bz2 URL ${CMAKE_SOURCE_DIR}/external/omptl-20120422.tar.bz2
CONFIGURE_COMMAND echo "No configure" CONFIGURE_COMMAND echo "No configure"
BUILD_COMMAND echo "No build" BUILD_COMMAND echo "No build"

BIN
external/libsharp-8d51946.tar.gz vendored
View File

Binary file not shown.

173
external/patch-omptl vendored
View File

@ -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 @@ @@ -20,7 +20,7 @@
#define OMPTL_ALGORITHM 1 #define OMPTL_ALGORITHM 1
@ -22,14 +22,12 @@ diff -ur omptl.old/omptl_algorithm omptl/omptl_algorithm
#endif #endif
#endif /* OMPTL_ALGORITHM */ #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 <cmath>
#include <cstdlib> #include <cstdlib>
+#include <random>
-#include <omptl/omptl_tools.h> -#include <omptl/omptl_tools.h>
-#include <omptl/omptl_numeric> -#include <omptl/omptl_numeric>
@ -38,136 +36,9 @@ diff -ur omptl.old/omptl_algorithm_par.h omptl/omptl_algorithm_par.h
#include <iterator> #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 @@ @@ -19,7 +19,7 @@
#define OMPTL_NUMERIC 1 #define OMPTL_NUMERIC 1
@ -192,9 +63,9 @@ diff -ur omptl.old/omptl_numeric omptl/omptl_numeric
+#include "omptl_numeric_extensions.h" +#include "omptl_numeric_extensions.h"
#endif /* OMPTL_NUMERIC */ #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 @@ @@ -51,9 +51,9 @@
} // namespace } // namespace
@ -207,9 +78,9 @@ diff -ur omptl.old/omptl_numeric_extensions.h omptl/omptl_numeric_extensions.h
#endif #endif
namespace omptl 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 @@ @@ -23,8 +23,8 @@
#include <functional> #include <functional>
#include <iterator> #include <iterator>
@ -221,15 +92,3 @@ diff -ur omptl.old/omptl_numeric_par.h omptl/omptl_numeric_par.h
namespace omptl 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
{

4
pyproject.toml
View File

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

22
python/CMakeLists.txt
View File

@ -2,7 +2,7 @@ set(CMAKE_SHARED_MODULE_PREFIX)
set(PYTHON_INCLUDES ${NUMPY_INCLUDE_DIRS} ${PYTHON_INCLUDE_PATH} ${CMAKE_SOURCE_DIR}/python) 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) IF(CYTHON)
add_custom_command( add_custom_command(
@ -54,10 +54,6 @@ SET(CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -Bsymbolic-functions
if(APPLE) if(APPLE)
set(CMAKE_MODULE_LINKER_FLAGS "-undefined dynamic_lookup") set(CMAKE_MODULE_LINKER_FLAGS "-undefined dynamic_lookup")
endif() 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(_cosmotool PRIVATE ${CosmoTool_local} ${GSL_LIBRARIES})
target_link_libraries(_cosmo_power 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) SET(ct_TARGETS _cosmotool _project _cosmo_power _cosmo_cic _fast_interp _cosmomath)
if (Boost_FOUND) if (Boost_FOUND)
message(STATUS "Building bispectrum support") message(STATUS "Building bispectrum support (path = ${Boost_INCLUDE_DIRS})")
include_directories(${Boost_INCLUDE_DIRS}) include_directories(${Boost_INCLUDE_DIRS})
add_library(_cosmo_bispectrum MODULE _cosmo_bispectrum.cpp) add_library(_cosmo_bispectrum MODULE _cosmo_bispectrum.cpp)
target_link_libraries(_cosmo_bispectrum ${MATH_LIBRARY}) 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) if (Boost_DEP)
add_dependencies(_cosmo_bispectrum ${Boost_DEP}) add_dependencies(_cosmo_bispectrum ${Boost_DEP})
endif() endif()
@ -113,17 +112,8 @@ configure_file(${CMAKE_CURRENT_SOURCE_DIR}/cosmotool/config.py.in ${CMAKE_CURREN
INSTALL(TARGETS INSTALL(TARGETS
${ct_TARGETS} ${ct_TARGETS}
COMPONENT python
LIBRARY DESTINATION ${PYTHON_SITE_PACKAGES}/cosmotool 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") 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")

11
python/_cosmomath.pyx
View File

@ -1,4 +1,3 @@
#cython: language_level=3
import numpy as np import numpy as np
cimport numpy as np cimport numpy as np
@ -8,21 +7,15 @@ np.import_ufunc()
cdef extern from "sys/types.h": cdef extern from "sys/types.h":
ctypedef np.int64_t int64_t ctypedef np.int64_t int64_t
cdef extern from "numpy/npy_common.h":
ctypedef npy_intp
cdef extern from "special_math.hpp" namespace "CosmoTool": 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 np.PyUFuncGenericFunction loop_func[1]
cdef char input_output_types[3] cdef char input_output_types[3]
cdef void *elementwise_funcs[1] 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[0] = np.NPY_DOUBLE
input_output_types[1] = np.NPY_DOUBLE input_output_types[1] = np.NPY_DOUBLE

17
python/_cosmotool.pyx
View File

@ -39,7 +39,7 @@ cdef extern from "loadSimu.hpp" namespace "CosmoTool":
cdef extern from "loadGadget.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 + void cxx_writeGadget "CosmoTool::writeGadget" (const char * s, SimuData *data) except +
cdef extern from "safe_gadget.hpp": cdef extern from "safe_gadget.hpp":
@ -247,9 +247,6 @@ class PySimulationAdaptor(PySimulationBase):
def __init__(self,sim): def __init__(self,sim):
self.simu = sim self.simu = sim
def getNumParticles(self):
return self.simu.numParticles
def getBoxsize(self): def getBoxsize(self):
return self.simu.BoxSize return self.simu.BoxSize
@ -313,13 +310,6 @@ tries to get an array from the object."""
references to the object are gone. """ references to the object are gone. """
pass 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 object wrap_array(void *p, np.uint64_t s, int typ):
cdef np.ndarray ndarray cdef np.ndarray ndarray
cdef ArrayWrapper wrapper cdef ArrayWrapper wrapper
@ -327,8 +317,7 @@ cdef object wrap_array(void *p, np.uint64_t s, int typ):
wrapper = ArrayWrapper() wrapper = ArrayWrapper()
wrapper.set_data(s, typ, p) wrapper.set_data(s, typ, p)
ndarray = np.array(wrapper, copy=False) ndarray = np.array(wrapper, copy=False)
#ndarray.base = <PyObject*> wrapper ndarray.base = <PyObject*> wrapper
PyArray_SetBaseObject(ndarray, <PyObject*> wrapper)
Py_INCREF(wrapper) Py_INCREF(wrapper)
return ndarray return ndarray
@ -421,7 +410,7 @@ def loadGadget(str filename, int snapshot_id, int gadgetFormat = 1, bool loadPos
with nogil: with nogil:
data = loadGadgetMulti(filename_bs, snapshot_id, flags, gadgetFormat) data = loadGadgetMulti(filename_bs, snapshot_id, flags, gadgetFormat)
if data == <SimuData*>0: if data == <SimuData*>0:
raise RuntimeError("File could not be read") return None
return PySimulationAdaptor(wrap_simudata(data, flags)) return PySimulationAdaptor(wrap_simudata(data, flags))

13
python/_project.pyx
View File

@ -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, 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. 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. 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 npx.ndarray[DTYPE_t, ndim=3] g
cdef double shifter 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]: 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") 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_g = g
cdef DTYPE_t[:,:] d_x = x 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 iu0[3]
cdef int i cdef int i
cdef int N = density.shape[0] 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 int completed
cdef DTYPE_t I0, d, dist2, delta, s, max_distance2 cdef DTYPE_t I0, d, dist2, delta, s, max_distance2
cdef int jumper[1] cdef int jumper[1]

7
python/cosmotool.version
View File

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

5
requirements.txt
View File

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

4
sample/CMakeLists.txt
View File

@ -1,8 +1,4 @@
SET(tolink ${CosmoTool_local} ${CosmoTool_LIBS} ${GSL_LIBRARIES} ${DL_LIBRARY}) 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(${CMAKE_SOURCE_DIR}/src)
include_directories(${FFTW3_INCLUDE_DIRS} ${EIGEN3_INCLUDE_DIRS} ${GSL_INCLUDE_PATH}) include_directories(${FFTW3_INCLUDE_DIRS} ${EIGEN3_INCLUDE_DIRS} ${GSL_INCLUDE_PATH})
if(YORICK_SUPPORT) if(YORICK_SUPPORT)

170
sample/simple3DFilter.cpp
View File

@ -1,22 +1,22 @@
#include "hdf5_array.hpp"
#include "miniargs.hpp"
#include "mykdtree.hpp"
#include "omptl/algorithm"
#include "openmp.hpp" #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 <cassert>
#include "yorick.hpp"
#include "sphSmooth.hpp"
#include "mykdtree.hpp"
#include "miniargs.hpp"
#include <H5Cpp.h>
#include "hdf5_array.hpp"
#include <iostream> #include <iostream>
#include <boost/format.hpp>
#include <boost/bind.hpp>
using namespace std; using namespace std;
using namespace CosmoTool; using namespace CosmoTool;
#define N_SPH 32 #define N_SPH 32
struct VCoord { struct VCoord{
float v[3]; float v[3];
float mass; float mass;
}; };
@ -27,42 +27,46 @@ typedef boost::multi_array<float, 2> array_type;
typedef boost::multi_array<float, 3> array3_type; typedef boost::multi_array<float, 3> array3_type;
typedef boost::multi_array<float, 4> array4_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 SPHSmooth<VCoord> MySmooth;
typedef MySmooth::SPHTree MyTree; typedef MySmooth::SPHTree MyTree;
typedef MyTree::Cell MyCell; 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); 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;
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)}; MyTree::coords c = { float(px), float(py), float(pz) };
double r2 = c[0] * c[0] + c[1] * c[1] + c[2] * c[2]; double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2];
if (r2 > rLimit2) { if (r2 > rLimit2)
{
arr[rx][ry][rz] = 0; arr[rx][ry][rz] = 0;
continue; continue;
} }
@ -80,34 +84,35 @@ void computeInterpolatedField(MyTree *tree1, double boxsize, int Nres,
} }
} }
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; double rLimit, boxsize, rLimit2, cx, cy, cz;
int Nres; 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)) if (!parseMiniArgs(argc, argv, args))
return 1; return 1;
H5::H5File in_f(fname1, 0); 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; array_type v1_data;
uint64_t N1_points, N2_points; uint32_t N1_points, N2_points;
array3_type bins(boost::extents[Nres][Nres][Nres]); array3_type bins(boost::extents[Nres][Nres][Nres]);
rLimit2 = rLimit * rLimit; rLimit2 = rLimit*rLimit;
hdf5_read_array(in_f, "particles", v1_data); hdf5_read_array(in_f, "particles", v1_data);
assert(v1_data.shape()[1] == 7); assert(v1_data.shape()[1] == 7);
@ -119,30 +124,30 @@ int main(int argc, char **argv) {
MyCell *allCells_1 = new MyCell[N1_points]; MyCell *allCells_1 = new MyCell[N1_points];
#pragma omp parallel for schedule(static) #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; bins.data()[i] = 0;
cout << "Shuffling data in cells..." << endl; cout << "Shuffling data in cells..." << endl;
#pragma omp parallel for schedule(static) #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++) for (int j = 0; j < 3; j++)
allCells_1[i].coord[j] = v1_data[i][j]; allCells_1[i].coord[j] = v1_data[i][j];
for (int k = 0; k < 3; k++) for (int k = 0; k < 3; k++)
allCells_1[i].val.pValue.v[k] = v1_data[i][3 + k]; allCells_1[i].val.pValue.v[k] = v1_data[i][3+k];
allCells_1[i].val.pValue.mass = v1_data[i][6]; allCells_1[i].val.pValue.mass = v1_data[i][6];
allCells_1[i].active = true; allCells_1[i].active = true;
allCells_1[i].val.weight = 0.0; 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]); v1_data.resize(boost::extents[1][1]);
@ -150,35 +155,35 @@ int main(int argc, char **argv) {
cout << "Building trees..." << endl; cout << "Building trees..." << endl;
MyTree tree1(allCells_1, N1_points); MyTree tree1(allCells_1, N1_points);
tree1.setPeriodic(periodic != 0, boxsize);
cout << "Creating smoothing filter..." << endl; 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; cout << "Weighing..." << endl;
int rz_max = 0; #pragma omp parallel
#pragma omp parallel shared(rz_max)
{ {
MySmooth smooth1(&tree1, N_SPH); 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 pz = (rz)*boxsize / Nres - cz; {
double py = (ry)*boxsize / Nres - cy; double py = (ry)*boxsize/Nres-cy;
double px = (rx)*boxsize / Nres - cx; for (int rx = 0; rx < Nres; rx++)
{
double px = (rx)*boxsize/Nres-cx;
MyTree::coords c = {float(px), float(py), float(pz)}; MyTree::coords c = { float(px), float(py), float(pz) };
double r2 = c[0] * c[0] + c[1] * c[1] + c[2] * c[2]; double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2];
if (r2 > rLimit2) { if (r2 > rLimit2)
{
continue; continue;
} }
@ -187,9 +192,11 @@ int main(int argc, char **argv) {
smooth1.fetchNeighbours(c, numInCell); smooth1.fetchNeighbours(c, numInCell);
else else
smooth1.fetchNeighbours(c); smooth1.fetchNeighbours(c);
#pragma omp critical
smooth1.addGridSite(c); smooth1.addGridSite(c);
} }
} }
(cout << " Done " << rz << endl).flush();
} }
} }
@ -197,14 +204,13 @@ int main(int argc, char **argv) {
array3_type interpolated(boost::extents[Nres][Nres][Nres]); 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); hdf5_write_array(out_f, "density", interpolated);
// out_f.flush(); //out_f.flush();
for (int i = 0; i < 3; i++) { 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);
} }

27
setup.py
View File

@ -2,7 +2,6 @@ import stat
import os import os
import sys import sys
import shutil import shutil
from sysconfig import get_config_var
from distutils.command.install_data import install_data from distutils.command.install_data import install_data
import pathlib import pathlib
from setuptools import find_packages, setup, Extension from setuptools import find_packages, setup, Extension
@ -168,33 +167,15 @@ class BuildCMakeExt(build_ext):
# Change your cmake arguments below as necessary # Change your cmake arguments below as necessary
# Below is just an example set of arguments for building Blender as a Python module # 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, 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', '-DENABLE_OPENMP=ON','-DINTERNAL_BOOST=ON','-DINTERNAL_EIGEN=ON',
'-DINTERNAL_HDF5=ON','-DINTERNAL_NETCDF=ON', '-DINTERNAL_HDF5=ON','-DINTERNAL_NETCDF=ON',
'-DBUILD_PYTHON=ON', '-DINSTALL_PYTHON_LOCAL=OFF', '-DBUILD_PYTHON=ON', '-DINSTALL_PYTHON_LOCAL=OFF',
'-DCOSMOTOOL_PYTHON_PACKAGING=ON', '-DCOSMOTOOL_PYTHON_PACKAGING=ON',
f"-DCYTHON={cython_code}", f"-DCYTHON={cython_code}",
f"-DPYTHON_SITE_PACKAGES={build_dir.absolute()}/private_install", 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) self.announce("Building binaries", level=3)
@ -237,10 +218,10 @@ class BuildCMakeExt(build_ext):
CosmoTool_extension = CMakeExtension(name="cosmotool") CosmoTool_extension = CMakeExtension(name="cosmotool")
setup(name='cosmotool', setup(name='cosmotool',
version='1.3.6', version='1.1.0',
packages=["cosmotool"], packages=["cosmotool"],
package_dir={'cosmotool': 'python/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'], setup_requires=['cython','cffi','numpy','numexpr'],
ext_modules=[CosmoTool_extension], ext_modules=[CosmoTool_extension],
description='A small cosmotool box of useful functions', description='A small cosmotool box of useful functions',

1
src/CMakeLists.txt
View File

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

45
src/cosmopower.cpp
View File

@ -212,30 +212,6 @@ double CosmoPower::powerHuWiggles(double k)
return normPower * pow(k,n) * T_k * T_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) double CosmoPower::primordialPowerSpectrum(double k)
{ {
//Primordial power spectrum, needed for PNG //Primordial power spectrum, needed for PNG
@ -273,19 +249,6 @@ double CosmoPower::matterTransferFunctionHu(double k)
return T_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 CosmoPower::powerHuBaryons(double k)
{ {
double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75))); 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; normPower = 1;
#if 0
ofstream ff("PP_k.txt"); ofstream ff("PP_k.txt");
for (int i = 0; i < 100; i++) for (int i = 0; i < 100; i++)
{ {
double k = pow(10.0, 8.0*i/100.-4); double k = pow(10.0, 8.0*i/100.-4);
ff << k << " " << power(k) << endl; ff << k << " " << power(k) << endl;
} }
#endif
// gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &normVal, &abserr); // 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); 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: case POWER_SUGIYAMA:
eval = &CosmoPower::powerSugiyama; eval = &CosmoPower::powerSugiyama;
break; break;
case SAMPLE_WIGGLES:
eval = &CosmoPower::sample_BAO;
break;
case POWER_BDM: case POWER_BDM:
eval = &CosmoPower::powerBDM; eval = &CosmoPower::powerBDM;
break; break;
case POWER_TEST: case POWER_TEST:
eval = &CosmoPower::powerTest; eval = &CosmoPower::powerTest;
break; break;
case BAO_TK:
eval = &CosmoPower::BAO_Tk;
break;
default: default:
abort(); abort();
} }

7
src/cosmopower.hpp
View File

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

10
src/fortran.cpp
View File

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

151
src/fourier/fft/fftw_calls.hpp
View File

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

190
src/fourier/fft/fftw_calls_mpi.hpp
View File

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

20
src/h5_readFlash.cpp
View File

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

10
src/hdf5_array.hpp
View File

@ -144,8 +144,8 @@ namespace CosmoTool {
bool useBases = false) bool useBases = false)
{ {
std::vector<hsize_t> memdims(data.shape(), data.shape() + data.num_dimensions()); 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) { if (useBases) {
std::vector<hsize_t> offsets(data.index_bases(), data.index_bases() + data.num_dimensions()); 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); hdf5_weak_check_array(data, dimensions);
std::vector<hsize_t> memdims(data.shape(), data.shape() + data.num_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()); std::vector<hsize_t> offsets(data.index_bases(), data.index_bases() + data.num_dimensions());
dataspace.selectHyperslab(H5S_SELECT_SET, memdims.data(), offsets.data()); dataspace.selectHyperslab(H5S_SELECT_SET, memdims.data(), offsets.data());
@ -426,7 +426,7 @@ namespace CosmoTool {
#define CTOOL_HDF5_INSERT_ELEMENT(r, STRUCT, element) \ #define CTOOL_HDF5_INSERT_ELEMENT(r, STRUCT, element) \
{ \ { \
::CosmoTool::get_hdf5_data_type<BOOST_PP_TUPLE_ELEM(2, 0, element)> t; \ ::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)); \ const char *field_name = BOOST_PP_STRINGIZE(BOOST_PP_TUPLE_ELEM(2, 1, element)); \
type.insertMember(field_name, position, t.type()); \ type.insertMember(field_name, position, t.type()); \
} }
@ -439,6 +439,7 @@ namespace CosmoTool { \
\ \
TNAME() : type(sizeof(STRUCT)) \ TNAME() : type(sizeof(STRUCT)) \
{ \ { \
long position; \
BOOST_PP_SEQ_FOR_EACH(CTOOL_HDF5_INSERT_ELEMENT, STRUCT, ATTRIBUTES) \ BOOST_PP_SEQ_FOR_EACH(CTOOL_HDF5_INSERT_ELEMENT, STRUCT, ATTRIBUTES) \
} \ } \
\ \
@ -470,6 +471,7 @@ namespace CosmoTool { \
\ \
TNAME() : type(sizeof(STRUCT)) \ TNAME() : type(sizeof(STRUCT)) \
{ \ { \
long position; \
BOOST_PP_SEQ_FOR_EACH(CTOOL_HDF5_INSERT_ENUM_ELEMENT, STRUCT, ATTRIBUTES) \ BOOST_PP_SEQ_FOR_EACH(CTOOL_HDF5_INSERT_ENUM_ELEMENT, STRUCT, ATTRIBUTES) \
} \ } \
\ \

1
src/interpolate.cpp
View File

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

10
src/loadGadget.cpp
View File

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

22
src/mykdtree.hpp
View File

@ -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 guilhem.lavaux@gmail.com
@ -99,8 +99,8 @@ namespace CosmoTool {
typename KDDef<N,CType>::CoordType r, r2; typename KDDef<N,CType>::CoordType r, r2;
KDCell<N, ValType,CType> **cells; KDCell<N, ValType,CType> **cells;
typename KDDef<N,CType>::CoordType *distances; 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, RecursionMultipleInfo(const typename KDDef<N,CType>::KDCoordinates& rx,
KDCell<N,ValType,CType> **cells, KDCell<N,ValType,CType> **cells,
uint64_t numCells) uint32_t numCells)
: queue(cells, numCells, INFINITY),traversed(0) : queue(cells, numCells, INFINITY),traversed(0)
{ {
std::copy(rx, rx+N, x); std::copy(rx, rx+N, x);
@ -153,20 +153,20 @@ namespace CosmoTool {
std::copy(replicate, replicate+N, this->replicate); 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, 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, Cell **cells,
CoordType *distances, 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); Cell *getNearestNeighbour(const coords& x);
void getNearestNeighbours(const coords& x, uint64_t NumCells, void getNearestNeighbours(const coords& x, uint32_t NumCells,
Cell **cells); Cell **cells);
void getNearestNeighbours(const coords& x, uint64_t NumCells, void getNearestNeighbours(const coords& x, uint32_t NumCells,
Cell **cells, Cell **cells,
CoordType *distances); CoordType *distances);

18
src/mykdtree.tcc
View File

@ -80,9 +80,9 @@ namespace CosmoTool {
} }
template<int N, typename ValType, typename CType, typename CellSplitter> 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, KDTree<N,ValType,CType,CellSplitter>::Cell **cells,
uint64_t numCells) uint32_t numCells)
{ {
RecursionInfoCells<N,ValType,CType> info; RecursionInfoCells<N,ValType,CType> info;
@ -112,10 +112,10 @@ namespace CosmoTool {
} }
template<int N, typename ValType, typename CType, typename CellSplitter> 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, Cell **cells,
CoordType *distances, CoordType *distances,
uint64_t numCells) uint32_t numCells)
{ {
RecursionInfoCells<N,ValType> info; RecursionInfoCells<N,ValType> info;
@ -144,7 +144,7 @@ namespace CosmoTool {
} }
template<int N, typename ValType, typename CType, typename CellSplitter> 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; RecursionInfoCells<N,ValType> info;
@ -501,7 +501,7 @@ namespace CosmoTool {
} }
template<int N, typename ValType, typename CType, typename CellSplitter> 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) Cell **cells)
{ {
RecursionMultipleInfo<N,ValType> info(x, cells, N2); RecursionMultipleInfo<N,ValType> info(x, cells, N2);
@ -527,7 +527,7 @@ namespace CosmoTool {
} }
template<int N, typename ValType, typename CType, typename CellSplitter> 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, Cell **cells,
CoordType *distances) CoordType *distances)
{ {
@ -643,8 +643,8 @@ namespace CosmoTool {
throw InvalidOnDiskKDTree(); throw InvalidOnDiskKDTree();
} }
if (node_on_disk.cell_id > numNodes || node_on_disk.cell_id < 0 || 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; delete[] nodes;
std::cerr << "Invalid cell id or children node id invalid" << std::endl; std::cerr << "Invalid cell id or children node id invalid" << std::endl;

30
src/numpy_adaptors.hpp
View File

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

26
src/powerSpectrum.cpp
View File

@ -58,7 +58,6 @@ using namespace std;
#define POWER_TEST 8 #define POWER_TEST 8
#define POWER_SPECTRUM POWER_EFSTATHIOU #define POWER_SPECTRUM POWER_EFSTATHIOU
#define SAMPLE_WIGGLES 9
namespace Cosmology { namespace Cosmology {
@ -207,31 +206,6 @@ double powG(double y)
*/ */
double powerSpectrum(double k, double normPower) 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 #if POWER_SPECTRUM == POWER_EFSTATHIOU
double a = 6.4/Gamma0; double a = 6.4/Gamma0;
double b = 3/Gamma0; double b = 3/Gamma0;

107
src/sphSmooth.hpp
View File

@ -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 guilhem.lavaux@gmail.com
@ -39,26 +39,30 @@ knowledge of the CeCILL license and that you accept its terms.
#include "config.hpp" #include "config.hpp"
#include "mykdtree.hpp" #include "mykdtree.hpp"
namespace CosmoTool { namespace CosmoTool
{
template <typename ValType, int Ndims = NUMDIMS> template <typename ValType, int Ndims = NUMDIMS>
class SPHSmooth { class SPHSmooth
{
public: public:
typedef struct { typedef struct
{
ComputePrecision weight; ComputePrecision weight;
ValType pValue; ValType pValue;
} FullType; } 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: public:
typedef SPHCell *P_SPHCell; typedef SPHCell *P_SPHCell;
struct SPHState { struct SPHState
{
boost::shared_ptr<P_SPHCell[]> ngb; boost::shared_ptr<P_SPHCell[]> ngb;
boost::shared_ptr<CoordType[]> distances; boost::shared_ptr<CoordType[]> distances;
typename SPHTree::coords currentCenter; typename SPHTree::coords currentCenter;
@ -66,58 +70,45 @@ namespace CosmoTool {
ComputePrecision smoothRadius; ComputePrecision smoothRadius;
}; };
SPHSmooth(SPHTree *tree, uint32_t Nsph); SPHSmooth(SPHTree *tree, uint32_t Nsph);
virtual ~SPHSmooth(); 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; 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 ! // TO USE WITH EXTREME CARE !
void setSmoothingLen(ComputePrecision len) { internal.smoothRadius = len; } void setSmoothingLen(ComputePrecision len)
{
internal.smoothRadius = len;
}
// END // END
template <typename FuncT> template<typename FuncT>
void runForEachNeighbour(FuncT fun, SPHState *state = 0); 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; bool hasNeighbours() const;
@ -142,21 +133,21 @@ namespace CosmoTool {
uint32_t maxNgb; uint32_t maxNgb;
SPHTree *tree; 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,
CoordType d,
FuncT fun, SPHState *state); 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" #include "sphSmooth.tcc"

247
src/sphSmooth.tcc
View File

@ -3,8 +3,9 @@
namespace CosmoTool { 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->Nsph = Nsph;
this->tree = tree; this->tree = tree;
internal.currentNgb = 0; internal.currentNgb = 0;
@ -12,16 +13,20 @@ namespace CosmoTool {
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,
FuncT fun, SPHState *state)
{
CoordType weight; CoordType weight;
d /= state->smoothRadius; d /= state->smoothRadius;
@ -31,47 +36,40 @@ namespace CosmoTool {
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; ComputePrecision d2, max_dist = 0;
uint32_t requested = newNngb; 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]);
} else {
state = &internal;
if (requested > maxNgb) {
maxNgb = requested; maxNgb = requested;
internal.ngb = boost::shared_ptr<P_SPHCell[]>(new P_SPHCell[maxNgb]); internal.ngb = boost::shared_ptr<P_SPHCell[]>(new P_SPHCell[maxNgb]);
internal.distances = internal.distances = boost::shared_ptr<CoordType[]>(new CoordType[maxNgb]);
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) if (d2 > max_dist)
max_dist = d2; 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;
@ -83,43 +81,48 @@ namespace CosmoTool {
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]); d2 = internal.distances[i] = sqrt(internal.distances[i]);
if (d2 > max_dist) if (d2 > max_dist)
max_dist = d2; max_dist = d2;
} }
state->smoothRadius = max_dist / 2; state->smoothRadius = max_dist / 2;
} }
template <typename ValType, int Ndims>
void SPHSmooth<ValType, Ndims>::fetchNeighboursOnVolume( template<typename ValType, int Ndims>
const typename SPHTree::coords &c, ComputePrecision radius) { void
SPHSmooth<ValType,Ndims>::fetchNeighboursOnVolume(const typename SPHTree::coords& c,
ComputePrecision radius)
{
uint32_t numPart; uint32_t numPart;
ComputePrecision d2, max_dist = 0; ComputePrecision d2, max_dist = 0;
memcpy(internal.currentCenter, c, sizeof(c)); memcpy(internal.currentCenter, c, sizeof(c));
internal.currentNgb = tree->getIntersection( internal.currentNgb = tree->getIntersection(c, radius, internal.ngb, internal.distances,
c, radius, internal.ngb, internal.distances, maxNgb); maxNgb);
for (uint32_t i = 0; i < internal.currentNgb; i++) { for (uint32_t i = 0; i < internal.currentNgb; i++)
{
d2 = internal.distances[i] = sqrt(internal.distances[i]); d2 = internal.distances[i] = sqrt(internal.distances[i]);
if (d2 > max_dist) if (d2 > max_dist)
max_dist = d2; max_dist = d2;
} }
internal.smoothRadius = max_dist / 2; internal.smoothRadius = max_dist / 2;
} }
template <typename ValType, int Ndims> template<typename ValType, int Ndims>
template <typename FuncT> template<typename FuncT>
ComputePrecision SPHSmooth<ValType, Ndims>::computeSmoothedValue( ComputePrecision
const typename SPHTree::coords &c, FuncT fun, SPHState *state) { SPHSmooth<ValType,Ndims>::computeSmoothedValue(const typename SPHTree::coords& c,
FuncT fun, SPHState *state)
{
if (state == 0) if (state == 0)
state = &internal; state = &internal;
@ -127,24 +130,26 @@ namespace CosmoTool {
ComputePrecision max_dist = 0; ComputePrecision max_dist = 0;
ComputePrecision r3 = cube(state->smoothRadius); ComputePrecision r3 = cube(state->smoothRadius);
for (uint32_t i = 0; i < state->currentNgb; i++) { for (uint32_t i = 0; i < state->currentNgb; i++)
outputValue += {
computeWValue(c, *state->ngb[i], state->distances[i], fun, state); 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) if (state == 0)
state = &internal; state = &internal;
@ -152,99 +157,71 @@ namespace CosmoTool {
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); 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]); fun(state->ngb[i]);
} }
} }
template <typename ValType, int Ndims>
void SPHSmooth<ValType, Ndims>::addGridSite( template<typename ValType, int Ndims>
const typename SPHTree::coords &c, SPHState *state) { void SPHSmooth<ValType,Ndims>::addGridSite(const typename SPHTree::coords& c)
{
ComputePrecision outputValue = 0; ComputePrecision outputValue = 0;
ComputePrecision max_dist = 0; ComputePrecision max_dist = 0;
ComputePrecision r3 = cube(state->smoothRadius);
for (uint32_t i = 0; i < state->currentNgb; i++) { ComputePrecision r3 = cube(internal.smoothRadius);
ComputePrecision d = state->distances[i];
SPHCell &cell = *(state->ngb[i]);
double kernel_value = getKernel(d / state->smoothRadius) / r3;
#pragma omp atomic update
cell.val.weight += kernel_value;
}
}
template <typename ValType, int Ndims> for (uint32_t i = 0; i < internal.currentNgb; i++)
void {
SPHSmooth<ValType, Ndims>::addGridSite(const typename SPHTree::coords &c) { ComputePrecision d = internal.distances[i];
addGridSite(c, &internal); SPHCell& cell = *(internal.ngb[i]);
cell.val.weight += getKernel(d/internal.smoothRadius) / r3;
} }
}
template <typename ValType, int Ndims> template<typename ValType, int Ndims>
ComputePrecision ComputePrecision
SPHSmooth<ValType, Ndims>::getKernel(ComputePrecision x) const { SPHSmooth<ValType,Ndims>::getKernel(ComputePrecision x) const
{
// WARNING !!! This is an unnormalized version of the kernel. // WARNING !!! This is an unnormalized version of the kernel.
if (x < 1) if (x < 1)
return 1 - 1.5 * x * x + 0.75 * x * x * x; return 1 - 1.5 * x * x + 0.75 * x * x * x;
else if (x < 2) { else if (x < 2)
{
ComputePrecision d = 2 - x; ComputePrecision d = 2 - x;
return 0.25 * d * d * d; return 0.25 * d * d * d;
} else }
else
return 0; 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
};

3
src/yorick.hpp
View File

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

1
src/yorick_nc3.cpp
View File

@ -260,6 +260,7 @@ namespace CosmoTool {
template<typename T> template<typename T>
void loadArray(const std::string& fname, void loadArray(const std::string& fname,
T*&array, uint32_t *&dimList, uint32_t& rank) T*&array, uint32_t *&dimList, uint32_t& rank)
throw (NoSuchFileException)
{ {
NcFile f(fname.c_str(), NcFile::ReadOnly); NcFile f(fname.c_str(), NcFile::ReadOnly);

1
src/yorick_nc4.cpp
View File

@ -224,6 +224,7 @@ namespace CosmoTool {
template<typename T> template<typename T>
void loadArray(const std::string& fname, void loadArray(const std::string& fname,
T*&array, uint32_t *&dimList, uint32_t& rank) T*&array, uint32_t *&dimList, uint32_t& rank)
throw (NoSuchFileException)
{ {
NcFile f(fname.c_str(), NcFile::read); NcFile f(fname.c_str(), NcFile::read);