diff --git a/CMakeLists.txt b/CMakeLists.txt index ecf5f9f..1693040 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -22,11 +22,13 @@ option(INTERNAL_HDF5 "Build internal version of HDF5" OFF) option(INTERNAL_NETCDF "Build internal version of NETCDF" OFF) option(INTERNAL_BOOST "Build internal version of BOOST" OFF) -IF(BUILD_SHARED_LIBS AND BUILD_STATIC_LIBS) +IF(NOT BUILD_SHARED_LIBS AND BUILD_STATIC_LIBS) SET(CosmoTool_local CosmoTool_static) -ELSE(BUILD_SHARED_LIBS AND BUILD_STATIC_LIBS) +ELSE(NOT BUILD_SHARED_LIBS AND BUILD_STATIC_LIBS) SET(CosmoTool_local CosmoTool) -ENDIF(BUILD_SHARED_LIBS AND BUILD_STATIC_LIBS) +ENDIF(NOT BUILD_SHARED_LIBS AND BUILD_STATIC_LIBS) + +MESSAGE(STATUS "Using the target ${CosmoTool_local} to build python module") include(${CMAKE_SOURCE_DIR}/external/external_build.cmake) @@ -40,11 +42,10 @@ ENDIF(EXISTS ${NETCDFCPP_INCLUDE_PATH}/netcdf AND ${NETCDFCPP_LIBRARY} MATCHES " find_program(CYTHON cython) -set(HDF5_FIND_COMPONENTS HL CXX) if(HDF5_ROOTDIR) SET(ENV{HDF5_ROOT} ${HDF5_ROOTDIR}) endif(HDF5_ROOTDIR) -include(FindHDF5) +find_package(HDF5 COMPONENTS HL CXX) set(NETCDF_FIND_REQUIRED TRUE) diff --git a/doc/source/pythonmodule.rst b/doc/source/pythonmodule.rst index bc5e42a..11fa942 100644 --- a/doc/source/pythonmodule.rst +++ b/doc/source/pythonmodule.rst @@ -41,3 +41,13 @@ Timing .. autofunction:: time_block .. autofunction:: timeit .. autofunction:: timeit_quiet + + +Cosmology +^^^^^^^^^ + +Power spectrum +-------------- + +.. autoclass:: CosmologyPower + :members: diff --git a/external/external_build.cmake b/external/external_build.cmake index edbf6ae..6a4ec87 100644 --- a/external/external_build.cmake +++ b/external/external_build.cmake @@ -100,7 +100,7 @@ if (INTERNAL_NETCDF) CONFIGURE_COMMAND ${NETCDF_SOURCE_DIR}/configure --prefix=${NETCDF_BIN_DIR} --libdir=${NETCDF_BIN_DIR}/lib --enable-netcdf-4 --with-pic --disable-shared --disable-dap - --disable-cdmremote --disable-rpc + --disable-cdmremote --disable-rpc --enable-cxx-4 --disable-examples ${EXTRA_NC_FLAGS} CC=${CMAKE_C_COMPILER} CXX=${CMAKE_CXX_COMPILER} BUILD_IN_SOURCE 1 @@ -110,7 +110,7 @@ if (INTERNAL_NETCDF) SET(EXTRA_NC_FLAGS CPPFLAGS=${CONFIGURE_CPP_FLAGS} LDFLAGS=${CONFIGURE_CPP_LDFLAGS}) SET(cosmotool_DEPS ${cosmotool_DEPS} netcdf) SET(NETCDF_LIBRARY ${NETCDF_BIN_DIR}/lib/libnetcdf.a CACHE STRING "NetCDF lib" FORCE) - SET(NETCDFCPP_LIBRARY ${NETCDF_BIN_DIR}/lib/libnetcdf_c++.a CACHE STRING "NetCDF-C++ lib" FORCE) + SET(NETCDFCPP_LIBRARY ${NETCDF_BIN_DIR}/lib/libnetcdf_c++4.a CACHE STRING "NetCDF-C++ lib" FORCE) SET(NETCDF_INCLUDE_PATH ${NETCDF_BIN_DIR}/include CACHE STRING "NetCDF include" FORCE) SET(NETCDFCPP_INCLUDE_PATH ${NETCDF_INCLUDE_PATH} CACHE STRING "NetCDF C++ include path" FORCE) diff --git a/python/_cosmo_power.pyx b/python/_cosmo_power.pyx index 9b46c82..4e9e871 100644 --- a/python/_cosmo_power.pyx +++ b/python/_cosmo_power.pyx @@ -49,11 +49,21 @@ cdef extern from "cosmopower.hpp" namespace "CosmoTool": double power(double) cdef class CosmologyPower: + """CosmologyPower(**cosmo) + + CosmologyPower manages and compute power spectra computation according to different + approximation given in the litterature. + + Keyword arguments: + omega_B_0 (float): relative baryon density + omega_M_0 (float): relative matter density + h (float): Hubble constant relative to 100 km/s/Mpc + ns (float): power law of the large scale inflation spectrum + """ cdef CosmoPower power - def __init__(self,**cosmo): - + def __init__(self,**cosmo): self.power = CosmoPower() self.power.OMEGA_B = cosmo['omega_B_0'] self.power.OMEGA_C = cosmo['omega_M_0']-cosmo['omega_B_0'] @@ -66,11 +76,26 @@ cdef class CosmologyPower: self.power.updateCosmology() def normalize(self,s8): + """normalize(self, sigma8) + + Compute the normalization of the power spectrum using sigma8. + + Arguments: + sigma8 (float): standard deviation of density field smoothed at 8 Mpc/h + """ self.power.SIGMA8 = s8 self.power.normalize() def setFunction(self,funcname): + """setFunction(self, funcname) + + Choose an approximation to use for the computation of the power spectrum + + Arguments: + funcname (str): the name of the approximation. It can be either + EFSTATHIOU, HU_WIGGLES, HU_BARYON, BARDEEN or SUGIYAMA. + """ cdef CosmoFunction f f = POWER_EFSTATHIOU @@ -93,6 +118,19 @@ cdef class CosmologyPower: return self.power.power(k) * self.power.h**3 def compute(self, k): + """compute(self, k) + + Compute the power spectrum for mode which length k. + + Arguments: + k (float): Mode for which to evaluate the power spectrum. + It can be a scalar or a numpy array. + The units must be in 'h Mpc^{-1}'. + + Returns: + a scalar or a numpy array depending on the type of the k argument + """ + cdef np.ndarray out cdef double kval cdef tuple i diff --git a/python/_project.pyx b/python/_project.pyx index 3ccf7ad..ae5b4b9 100644 --- a/python/_project.pyx +++ b/python/_project.pyx @@ -808,7 +808,25 @@ def spherical_projection(int Nside, npx.ndarray[DTYPE_t, ndim=3] density not None, DTYPE_t min_distance, DTYPE_t max_distance, int progress=1, int integrator_id=0, DTYPE_t[:] shifter = None, int booster=-1): - + """ + spherical_projection(Nside, density, min_distance, max_distance, progress=1, integrator_id=0, shifter=None, booster=-1) + + Keyword arguments: + progress (int): show progress if it is equal to 1 + integrator_id (int): specify the order of integration along the line of shift + shifter (DTYPE_t array): this is an array of size 3. It specifies the amount of shift to apply to the center, in unit of voxel + booster (int): what is the frequency of refreshment of the progress bar. Small number decreases performance by locking the GIL. + + + Arguments: + Nside (int): Nside of the returned map + density (NxNxN array): this is the density field, expressed as a cubic array + min_distance (float): lower bound of the integration + max_distance (float): upper bound of the integration + + Returns: + an healpix map, as a 1-dimensional array. + """ import healpy as hp import progressbar as pb cdef int i diff --git a/python/cosmotool/borg.py b/python/cosmotool/borg.py index df10fae..2deac20 100644 --- a/python/cosmotool/borg.py +++ b/python/cosmotool/borg.py @@ -1,7 +1,7 @@ ### ### BORG code is from J. Jasche ### -import StringIO +import io import numpy as np from numpy import * import os.path @@ -154,9 +154,11 @@ def get_grid_values(xx,data, ranges): def get_mean_density(fdir, smin, step): """ estimate ensemble mean """ - print '-'*60 - print 'Get 3D ensemble mean density field' - print '-'*60 + import progressbar as pb + + print('-'*60) + print('Get 3D ensemble mean density field') + print('-'*60) fname0 = fdir + 'initial_density_'+str(0)+'.dat' fname1 = fdir + 'final_density_'+str(0)+'.dat' @@ -174,23 +176,23 @@ def get_mean_density(fdir, smin, step): fname1 = fdir + 'final_density_'+str(idat)+'.dat' #and (idat #include "yorick.hpp" @@ -53,7 +54,7 @@ void computeInterpolatedField(MyTree *tree1, double boxsize, int Nres, double cx { double pz = (rz)*boxsize/Nres-cz; - cout << format("[%d] %d / %d") % omp_get_thread_num() % rz % Nres << endl; + cout << format("[%d] %d / %d") % smp_get_thread_id() % rz % Nres << endl; for (int ry = 0; ry < Nres; ry++) { double py = (ry)*boxsize/Nres-cy; diff --git a/sample/testHDF5.cpp b/sample/testHDF5.cpp index d38eb01..4011889 100644 --- a/sample/testHDF5.cpp +++ b/sample/testHDF5.cpp @@ -38,10 +38,37 @@ knowledge of the CeCILL license and that you accept its terms. using namespace std; +struct MyStruct +{ + int a; + double b; + char c; +}; + +struct MyStruct2 +{ + MyStruct base; + int d; +}; + + +CTOOL_STRUCT_TYPE(MyStruct, hdf5t_MyStruct, + ((int, a)) + ((double, b)) + ((char, c)) +) + +CTOOL_STRUCT_TYPE(MyStruct2, hdf5t_MyStruct2, + ((MyStruct, base)) + ((int, d)) +) + int main() { typedef boost::multi_array array_type; typedef boost::multi_array array3_type; + typedef boost::multi_array array_mys_type; + typedef boost::multi_array array_mys2_type; typedef boost::multi_array, 2> arrayc_type; typedef array_type::index index; @@ -53,13 +80,27 @@ int main() array_type B; array3_type C(boost::extents[2][3][4]); arrayc_type D, E; + array_mys_type F(boost::extents[10]), G; + array_mys2_type H(boost::extents[10]); int values = 0; for (index i = 0; i != 2; i++) for (index j = 0; j != 3; j++) A[i][j] = values++; + for (index i = 0; i != 10; i++) + { + F[i].a = i; + F[i].b = double(i)/4.; + F[i].c = 'r'+i; + H[i].base = F[i]; + H[i].d = 2*i; + } + std::cout << " c = " << ((char *)&F[1])[offsetof(MyStruct, c)] << endl; + CosmoTool::hdf5_write_array(g, "test_data", A); + CosmoTool::hdf5_write_array(g, "test_struct", F); + CosmoTool::hdf5_write_array(g, "test_struct2", H); CosmoTool::hdf5_read_array(g, "test_data", B); @@ -95,5 +136,12 @@ int main() abort(); } + CosmoTool::hdf5_read_array(g, "test_struct", G); + for (index i = 0; i != 10; i++) + if (G[i].a != F[i].a || G[i].b != F[i].b || G[i].c != F[i].c) { + std::cout << "Invalid struct content" << endl; + abort(); + } + return 0; } diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 9737765..5d42ce4 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,4 +1,3 @@ -add_definitions(-fPIC) SET(CosmoTool_SRCS fortran.cpp @@ -6,7 +5,6 @@ SET(CosmoTool_SRCS load_data.cpp loadGadget.cpp loadRamses.cpp - octTree.cpp powerSpectrum.cpp miniargs.cpp growthFactor.cpp @@ -76,9 +74,11 @@ if (BUILD_SHARED_LIBS) target_link_libraries(CosmoTool ${CosmoTool_LIBS}) if (BUILD_STATIC_LIBS) add_library(CosmoTool_static STATIC ${CosmoTool_SRCS}) + set_target_properties(CosmoTool_static PROPERTIES COMPILE_FLAGS "${CMAKE_C_COMPILE_OPTIONS_PIC}") endif(BUILD_STATIC_LIBS) else (BUILD_SHARED_LIBS) add_library(CosmoTool STATIC ${CosmoTool_SRCS}) + set_target_properties(CosmoTool PROPERTIES COMPILE_FLAGS "${CMAKE_C_COMPILE_OPTIONS_PIC}") endif (BUILD_SHARED_LIBS) install(TARGETS CosmoTool diff --git a/src/hdf5_array.hpp b/src/hdf5_array.hpp index 455d4b7..abcdc79 100644 --- a/src/hdf5_array.hpp +++ b/src/hdf5_array.hpp @@ -44,6 +44,9 @@ knowledge of the CeCILL license and that you accept its terms. #include #include #include +#include +#include +#include #include #include @@ -59,7 +62,7 @@ namespace CosmoTool { template struct get_hdf5_data_type { - static H5::PredType type() + static H5::DataType type() { BOOST_MPL_ASSERT_MSG(0, Unknown_HDF5_data_type, ()); return H5::PredType::NATIVE_DOUBLE; @@ -68,7 +71,7 @@ namespace CosmoTool { #define HDF5_TYPE(tl, thdf5) \ template struct get_hdf5_data_type >::type> \ - { static H5::PredType type() { return H5::PredType::thdf5; }; } + { static H5::DataType type() { return H5::PredType::thdf5; }; } HDF5_TYPE(char , NATIVE_CHAR); HDF5_TYPE(long long , NATIVE_LLONG); @@ -202,6 +205,39 @@ namespace CosmoTool { hdf5_read_array(fg, data_set_name, data, hdf5_ComplexType::ctype()->type); } + +#define CTOOL_HDF5_NAME(STRUCT) BOOST_PP_CAT(hdf5_,STRUCT) +#define CTOOL_HDF5_INSERT_ELEMENT(r, STRUCT, element) \ + { \ + ::CosmoTool::get_hdf5_data_type t; \ + position = HOFFSET(STRUCT, BOOST_PP_TUPLE_ELEM(2, 1, element)); \ + const char *field_name = BOOST_PP_STRINGIZE(BOOST_PP_TUPLE_ELEM(2, 1, element)); \ + type.insertMember(field_name, position, t.type()); \ + } + +#define CTOOL_STRUCT_TYPE(STRUCT, TNAME, ATTRIBUTES) \ +namespace CosmoTool { \ + class TNAME { \ + public: \ + H5::CompType type; \ + \ + TNAME() : type(sizeof(STRUCT)) \ + { \ + long position; \ + BOOST_PP_SEQ_FOR_EACH(CTOOL_HDF5_INSERT_ELEMENT, STRUCT, ATTRIBUTES) \ + } \ + \ + static const TNAME *ctype() \ + { \ + static TNAME singleton; \ + return &singleton; \ + } \ + }; \ + template<> struct get_hdf5_data_type { \ + static H5::DataType type() { return TNAME::ctype()->type; }; \ + }; \ +}; + }; #endif diff --git a/src/octTree.hpp b/src/octTree.hpp index 877017a..87d7dfa 100644 --- a/src/octTree.hpp +++ b/src/octTree.hpp @@ -55,12 +55,20 @@ namespace CosmoTool typedef octCoordType OctCoords[3]; + template struct OctCell { octPtr numberLeaves; octPtr children[8]; + T data; }; + template + struct OctTree_defaultUpdater { + void operator()(T& d) { } + }; + + template, class T = void> class OctTree { public: @@ -103,9 +111,10 @@ namespace CosmoTool protected: + T_dataUpdater updater; const FCoordinates *particles; octPtr numParticles; - OctCell *cells; + OctCell *cells; float Lbox; octPtr lastNode; octPtr numCells; @@ -128,47 +137,47 @@ namespace CosmoTool FCoordinates center, realCenter; for (int j = 0; j < 3; j++) - { - center[j] = icoord[j]/(2.*octCoordCenter); - realCenter[j] = xMin[j] + center[j]*lenNorm; - } + { + center[j] = icoord[j]/(2.*octCoordCenter); + realCenter[j] = xMin[j] + center[j]*lenNorm; + } f(realCenter, cells[node].numberLeaves, lenNorm*halfNodeLength/(float)octCoordCenter, - cells[node].children[0] == invalidOctCell, // True if this is a meta-node - false); + cells[node].children[0] == invalidOctCell, // True if this is a meta-node + false); if (!condition(realCenter, cells[node].numberLeaves, lenNorm*halfNodeLength/(float)octCoordCenter, cells[node].children[0] == invalidOctCell)) - return; + return; for (int i = 0; i < 8; i++) - { - octPtr newNode = cells[node].children[i]; - int ipos[3] = { (i&1), (i&2)>>1, (i&4)>>2 }; - - if (newNode == emptyOctCell || newNode == invalidOctCell) - continue; - - for (int j = 0; j < 3; j++) - newCoord[j] = icoord[j]+(2*ipos[j]-1)*halfNodeLength/2; - - if (newNode & octParticleMarker) - { - for (int j = 0; j < 3; j++) - { - center[j] = newCoord[j]/(2.*octCoordCenter); - realCenter[j] = xMin[j] + lenNorm*center[j]; - } + { + octPtr newNode = cells[node].children[i]; + int ipos[3] = { (i&1), (i&2)>>1, (i&4)>>2 }; + + if (newNode == emptyOctCell || newNode == invalidOctCell) + continue; + + for (int j = 0; j < 3; j++) + newCoord[j] = icoord[j]+(2*ipos[j]-1)*halfNodeLength/2; + + if (newNode & octParticleMarker) + { + for (int j = 0; j < 3; j++) + { + center[j] = newCoord[j]/(2.*octCoordCenter); + realCenter[j] = xMin[j] + lenNorm*center[j]; + } - f(realCenter, - 1, lenNorm*halfNodeLength/(2.*octCoordCenter), - false, true); - continue; - } + f(realCenter, + 1, lenNorm*halfNodeLength/(2.*octCoordCenter), + false, true); + continue; + } - walkTreeElements(f, condition, cells[node].children[i], newCoord, halfNodeLength/2); - } + walkTreeElements(f, condition, cells[node].children[i], newCoord, halfNodeLength/2); + } } @@ -177,4 +186,6 @@ namespace CosmoTool }; +#include "octTree.tcc" + #endif diff --git a/src/octTree.cpp b/src/octTree.tcc similarity index 92% rename from src/octTree.cpp rename to src/octTree.tcc index 9f1d5f7..e685ea2 100644 --- a/src/octTree.cpp +++ b/src/octTree.tcc @@ -39,8 +39,9 @@ knowledge of the CeCILL license and that you accept its terms. #include "config.hpp" #include "octTree.hpp" +namespace CosmoTool { + using namespace std; -using namespace CosmoTool; //#define VERBOSE @@ -59,7 +60,8 @@ static uint32_t mypow(uint32_t i, uint32_t p) return j*j*i; } -OctTree::OctTree(const FCoordinates *particles, octPtr numParticles, +template +OctTree::OctTree(const FCoordinates *particles, octPtr numParticles, uint32_t maxMeanTreeDepth, uint32_t maxAbsoluteDepth, uint32_t threshold) { @@ -94,7 +96,7 @@ OctTree::OctTree(const FCoordinates *particles, octPtr numParticles, } cout << xMin[0] << " " << xMin[1] << " " << xMin[2] << " lNorm=" << lenNorm << endl; - cells = new OctCell[numCells]; + cells = new OctCell[numCells]; Lbox = (float)(octCoordTypeNorm+1); cells[0].numberLeaves = 0; @@ -110,12 +112,14 @@ OctTree::OctTree(const FCoordinates *particles, octPtr numParticles, //#endif } -OctTree::~OctTree() +template +OctTree::~OctTree() { delete cells; } -void OctTree::buildTree(uint32_t maxAbsoluteDepth) +template +void OctTree::buildTree(uint32_t maxAbsoluteDepth) { for (octPtr i = 0; i < numParticles; i++) { @@ -129,7 +133,8 @@ void OctTree::buildTree(uint32_t maxAbsoluteDepth) } -void OctTree::insertParticle(octPtr node, +template +void OctTree::insertParticle(octPtr node, const OctCoords& icoord, octCoordType halfNodeLength, octPtr particleId, @@ -208,3 +213,5 @@ void OctTree::insertParticle(octPtr node, cells[node].children[octPos] = newNode; } + +}; diff --git a/src/openmp.hpp b/src/openmp.hpp new file mode 100644 index 0000000..7ef2c76 --- /dev/null +++ b/src/openmp.hpp @@ -0,0 +1,44 @@ +#ifndef __CTOOL_OPENMP_HPP +#define __CTOOL_OPENMP_HPP + +#ifdef _OPENMP +#include +#endif + +namespace CosmoTool { + + static int smp_get_max_threads() { +#ifdef _OPENMP + return omp_get_max_threads(); +#else + return 1; +#endif + } + + static int smp_get_thread_id() { +#ifdef _OPENMP + return omp_get_thread_num(); +#else + return 0; +#endif + } + + static int smp_get_num_threads() { +#ifdef _OPENMP + return omp_get_num_threads(); +#else + return 1; +#endif + + } + + static void smp_set_nested(bool n) { +#ifdef _OPENMP + omp_set_nested(n ? 1 : 0); +#endif + } + + +}; + +#endif diff --git a/src/yorick_nc3.cpp b/src/yorick_nc3.cpp index 30a8173..b1f09e9 100644 --- a/src/yorick_nc3.cpp +++ b/src/yorick_nc3.cpp @@ -37,7 +37,7 @@ knowledge of the CeCILL license and that you accept its terms. #include "config.hpp" #ifdef NETCDFCPP4 #include -using namespace netCDF +using namespace netCDF; #else #include #endif