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

Merge branch 'master' of bitbucket.org:glavaux/cosmotool

Browse Source
This commit is contained in:
Guilhem Lavaux 2022-01-27 13:36:25 +01:00
commit 0b77b010e4
11 changed files with 243 additions and 226 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
CMakeLists.txt
View File

@ -70,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 "2") SET(CPACK_PACKAGE_VERSION_MINOR "3")
SET(CPACK_PACKAGE_VERSION_PATCH "3${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

2
conda/meta.yaml
View File

@ -1,6 +1,6 @@
package: package:
name: cosmotool name: cosmotool
version: "1.2.3" version: "1.3.0"
source: source:
git_rev: a86c9a8 git_rev: a86c9a8

35
python/CMakeLists.txt
View File

@ -8,32 +8,32 @@ IF(CYTHON)
add_custom_command( add_custom_command(
OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/_cosmotool.cpp OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/_cosmotool.cpp
COMMAND ${CYTHON} --cplus -o ${CMAKE_CURRENT_BINARY_DIR}/_cosmotool.cpp ${CMAKE_CURRENT_SOURCE_DIR}/_cosmotool.pyx COMMAND ${CYTHON} --cplus -o ${CMAKE_CURRENT_BINARY_DIR}/_cosmotool.cpp ${CMAKE_CURRENT_SOURCE_DIR}/_cosmotool.pyx
DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/_cosmotool.pyx) DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/_cosmotool.pyx)
add_custom_command( add_custom_command(
OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/_cosmo_power.cpp OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/_cosmo_power.cpp
COMMAND ${CYTHON} --cplus -o ${CMAKE_CURRENT_BINARY_DIR}/_cosmo_power.cpp ${CMAKE_CURRENT_SOURCE_DIR}/_cosmo_power.pyx COMMAND ${CYTHON} --cplus -o ${CMAKE_CURRENT_BINARY_DIR}/_cosmo_power.cpp ${CMAKE_CURRENT_SOURCE_DIR}/_cosmo_power.pyx
DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/_cosmo_power.pyx) DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/_cosmo_power.pyx)
add_custom_command( add_custom_command(
OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/_fast_interp.cpp OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/_fast_interp.cpp
COMMAND ${CYTHON} --cplus -o ${CMAKE_CURRENT_BINARY_DIR}/_fast_interp.cpp ${CMAKE_CURRENT_SOURCE_DIR}/_fast_interp.pyx COMMAND ${CYTHON} --cplus -o ${CMAKE_CURRENT_BINARY_DIR}/_fast_interp.cpp ${CMAKE_CURRENT_SOURCE_DIR}/_fast_interp.pyx
DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/_fast_interp.pyx) DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/_fast_interp.pyx)
add_custom_command( add_custom_command(
OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/_cosmo_cic.cpp OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/_cosmo_cic.cpp
COMMAND ${CYTHON} --cplus -o ${CMAKE_CURRENT_BINARY_DIR}/_cosmo_cic.cpp ${CMAKE_CURRENT_SOURCE_DIR}/_cosmo_cic.pyx COMMAND ${CYTHON} --cplus -o ${CMAKE_CURRENT_BINARY_DIR}/_cosmo_cic.cpp ${CMAKE_CURRENT_SOURCE_DIR}/_cosmo_cic.pyx
DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/_cosmo_cic.pyx) DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/_cosmo_cic.pyx)
add_custom_command( add_custom_command(
OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/_project.cpp OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/_project.cpp
COMMAND ${CYTHON} --cplus -o ${CMAKE_CURRENT_BINARY_DIR}/_project.cpp ${CMAKE_CURRENT_SOURCE_DIR}/_project.pyx COMMAND ${CYTHON} --cplus -o ${CMAKE_CURRENT_BINARY_DIR}/_project.cpp ${CMAKE_CURRENT_SOURCE_DIR}/_project.pyx
DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/_project.pyx ${CMAKE_CURRENT_SOURCE_DIR}/project_tool.hpp ) DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/_project.pyx ${CMAKE_CURRENT_SOURCE_DIR}/project_tool.hpp )
add_custom_command( add_custom_command(
OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/_cosmomath.cpp OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/_cosmomath.cpp
COMMAND ${CYTHON} --cplus -o ${CMAKE_CURRENT_BINARY_DIR}/_cosmomath.cpp ${CMAKE_CURRENT_SOURCE_DIR}/_cosmomath.pyx COMMAND ${CYTHON} --cplus -o ${CMAKE_CURRENT_BINARY_DIR}/_cosmomath.cpp ${CMAKE_CURRENT_SOURCE_DIR}/_cosmomath.pyx
DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/_cosmomath.pyx ) DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/_cosmomath.pyx )
ENDIF(CYTHON) ENDIF(CYTHON)
@ -43,11 +43,11 @@ add_library(_cosmo_cic MODULE ${CMAKE_CURRENT_BINARY_DIR}/_cosmo_cic.cpp)
add_library(_fast_interp MODULE ${CMAKE_CURRENT_BINARY_DIR}/_fast_interp.cpp) add_library(_fast_interp MODULE ${CMAKE_CURRENT_BINARY_DIR}/_fast_interp.cpp)
add_library(_project MODULE ${CMAKE_CURRENT_BINARY_DIR}/_project.cpp) add_library(_project MODULE ${CMAKE_CURRENT_BINARY_DIR}/_project.cpp)
add_library(_cosmomath MODULE ${CMAKE_CURRENT_BINARY_DIR}/_cosmomath.cpp) add_library(_cosmomath MODULE ${CMAKE_CURRENT_BINARY_DIR}/_cosmomath.cpp)
target_include_directories(_cosmotool PRIVATE ${PYTHON_INCLUDES}) target_include_directories(_cosmotool PRIVATE ${PYTHON_INCLUDES})
target_include_directories(_cosmo_power PRIVATE ${PYTHON_INCLUDES}) target_include_directories(_cosmo_power PRIVATE ${PYTHON_INCLUDES})
target_include_directories(_cosmo_cic PRIVATE ${PYTHON_INCLUDES}) target_include_directories(_cosmo_cic PRIVATE ${PYTHON_INCLUDES})
target_include_directories(_fast_interp PRIVATE ${PYTHON_INCLUDES}) target_include_directories(_fast_interp PRIVATE ${PYTHON_INCLUDES})
target_include_directories(_project PRIVATE ${PYTHON_INCLUDES}) target_include_directories(_project PRIVATE ${PYTHON_INCLUDES})
target_include_directories(_cosmomath PRIVATE ${PYTHON_INCLUDES}) target_include_directories(_cosmomath PRIVATE ${PYTHON_INCLUDES})
SET(CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -Bsymbolic-functions") SET(CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -Bsymbolic-functions")
@ -114,10 +114,19 @@ endif (WIN32 AND NOT CYGWIN)
configure_file(${CMAKE_CURRENT_SOURCE_DIR}/cosmotool/config.py.in ${CMAKE_CURRENT_BINARY_DIR}/cosmotool/config.py @ONLY) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/cosmotool/config.py.in ${CMAKE_CURRENT_BINARY_DIR}/cosmotool/config.py @ONLY)
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 DESTINATION ${PYTHON_SITE_PACKAGES} INSTALL(DIRECTORY cosmotool ${CMAKE_CURRENT_BINARY_DIR}/cosmotool
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")

137
python/_cosmotool.pyx
View File

@ -31,7 +31,7 @@ cdef extern from "loadSimu.hpp" namespace "CosmoTool":
bool noAuto bool noAuto
cdef const int NEED_GADGET_ID cdef const int NEED_GADGET_ID
cdef const int NEED_POSITION cdef const int NEED_POSITION
cdef const int NEED_VELOCITY cdef const int NEED_VELOCITY
cdef const int NEED_TYPE cdef const int NEED_TYPE
@ -57,36 +57,36 @@ class PySimulationBase(object):
""" """
This is the base class to representation Simulation in CosmoTool/python. This is the base class to representation Simulation in CosmoTool/python.
""" """
def getPositions(self): def getPositions(self):
""" """
getPositions(self) getPositions(self)
Returns: Returns:
A list of three arrays holding the positions of the particles. A list of three arrays holding the positions of the particles.
The i-th element is the i-th coordinate of each particle. The i-th element is the i-th coordinate of each particle.
It may be None if the positions were not requested. It may be None if the positions were not requested.
""" """
raise NotImplementedError("getPositions is not implemented") raise NotImplementedError("getPositions is not implemented")
def getVelocities(self): def getVelocities(self):
""" """
getVelocities(self) getVelocities(self)
Returns: Returns:
A list of three arrays holding the velocities of the particles. A list of three arrays holding the velocities of the particles.
The i-th element is the i-th coordinate of each particle. The i-th element is the i-th coordinate of each particle.
It may be None if the velocities were not requested. It may be None if the velocities were not requested.
""" """
raise NotImplementedError("getVelocities is not implemented") raise NotImplementedError("getVelocities is not implemented")
def getIdentifiers(self): def getIdentifiers(self):
""" """
getIdentifiers(self) getIdentifiers(self)
Returns: Returns:
Returns an integer array that hold the unique identifiers of Returns an integer array that hold the unique identifiers of
each particle. each particle.
It may be None if the identifiers were not requested. It may be None if the identifiers were not requested.
""" """
raise NotImplementedError("getIdentifiers is not implemented") raise NotImplementedError("getIdentifiers is not implemented")
@ -94,10 +94,10 @@ class PySimulationBase(object):
def getTypes(self): def getTypes(self):
""" """
getTypes(self) getTypes(self)
Returns: Returns:
Returns an integer array that hold the type of Returns an integer array that hold the type of
each particle. each particle.
It may be None if the types were not requested. It may be None if the types were not requested.
""" """
raise NotImplementedError("getTypes is not implemented") raise NotImplementedError("getTypes is not implemented")
@ -105,27 +105,27 @@ class PySimulationBase(object):
def getOmega_M(self): def getOmega_M(self):
""" """
getOmega_M(self) getOmega_M(self)
Returns: Returns:
the mean matter density in the simulation, with respect the mean matter density in the simulation, with respect
to the critical density. to the critical density.
""" """
raise NotImplementedError("getOmega_M is not implemented") raise NotImplementedError("getOmega_M is not implemented")
def getOmega_Lambda(self): def getOmega_Lambda(self):
""" """
getOmega_Lambda(self) getOmega_Lambda(self)
Returns: Returns:
the mean dark energy density in the simulation, with respect the mean dark energy density in the simulation, with respect
to the critical density. to the critical density.
""" """
raise NotImplementedError("getOmega_Lambda is not implemented") raise NotImplementedError("getOmega_Lambda is not implemented")
def getTime(self): def getTime(self):
""" """
getTime(self) getTime(self)
Returns: Returns:
the time the snapshot was taken in the simulation. It can the time the snapshot was taken in the simulation. It can
have various units depending on the file format. have various units depending on the file format.
@ -135,7 +135,7 @@ class PySimulationBase(object):
def getHubble(self): def getHubble(self):
""" """
getHubble(self) getHubble(self)
Returns: Returns:
the hubble constant in unit of 100 km/s/Mpc the hubble constant in unit of 100 km/s/Mpc
""" """
@ -144,7 +144,7 @@ class PySimulationBase(object):
def getBoxsize(self): def getBoxsize(self):
""" """
getBoxsize(self) getBoxsize(self)
Returns: Returns:
the size of the simulation box. The length unit is not fixed, the size of the simulation box. The length unit is not fixed,
though it is customary to have it in Mpc/h if the loader has though it is customary to have it in Mpc/h if the loader has
@ -155,7 +155,7 @@ class PySimulationBase(object):
def getMasses(self): def getMasses(self):
""" """
getMasses(self) getMasses(self)
Returns: Returns:
an array with the masses of each particles, in unspecified unit that an array with the masses of each particles, in unspecified unit that
depend on the loader. depend on the loader.
@ -165,7 +165,7 @@ class PySimulationBase(object):
cdef class Simulation: cdef class Simulation:
""" """
Simulation() Simulation()
Class that directly manages internal loaded data obtained from a loader Class that directly manages internal loaded data obtained from a loader
""" """
@ -180,7 +180,7 @@ cdef class Simulation:
property BoxSize: property BoxSize:
def __get__(Simulation self): def __get__(Simulation self):
return self.data.BoxSize return self.data.BoxSize
property time: property time:
def __get__(Simulation self): def __get__(Simulation self):
return self.data.time return self.data.time
@ -192,15 +192,15 @@ cdef class Simulation:
property Omega_M: property Omega_M:
def __get__(Simulation self): def __get__(Simulation self):
return self.data.Omega_M return self.data.Omega_M
property Omega_Lambda: property Omega_Lambda:
def __get__(Simulation self): def __get__(Simulation self):
return self.data.Omega_Lambda return self.data.Omega_Lambda
property positions: property positions:
def __get__(Simulation self): def __get__(Simulation self):
return self.positions return self.positions
property velocities: property velocities:
def __get__(Simulation self): def __get__(Simulation self):
return self.velocities return self.velocities
@ -216,7 +216,7 @@ cdef class Simulation:
property masses: property masses:
def __get__(Simulation self): def __get__(Simulation self):
return self.masses return self.masses
property numParticles: property numParticles:
def __get__(Simulation self): def __get__(Simulation self):
return self.data.NumPart return self.data.NumPart
@ -228,7 +228,7 @@ cdef class Simulation:
def __cinit__(Simulation self): def __cinit__(Simulation self):
self.data = <SimuData *>0 self.data = <SimuData *>0
def __dealloc__(Simulation self): def __dealloc__(Simulation self):
if self.data != <SimuData *>0: if self.data != <SimuData *>0:
del self.data del self.data
@ -237,40 +237,43 @@ cdef class Simulation:
class PySimulationAdaptor(PySimulationBase): class PySimulationAdaptor(PySimulationBase):
""" """
PySimulationAdaptor(PySimulationBase_) PySimulationAdaptor(PySimulationBase_)
This class is an adaptor for an internal type to the loader. It defines This class is an adaptor for an internal type to the loader. It defines
all the methods of PySimulationBase. all the methods of PySimulationBase.
Attributes: Attributes:
simu: a Simulation_ object simu: a Simulation_ object
""" """
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
def getPositions(self): def getPositions(self):
return self.simu.positions return self.simu.positions
def getTypes(self): def getTypes(self):
return self.simu.types return self.simu.types
def getVelocities(self): def getVelocities(self):
return self.simu.velocities return self.simu.velocities
def getIdentifiers(self): def getIdentifiers(self):
return self.simu.identifiers return self.simu.identifiers
def getTime(self): def getTime(self):
return self.simu.time return self.simu.time
def getHubble(self): def getHubble(self):
return self.simu.Hubble return self.simu.Hubble
def getOmega_M(self): def getOmega_M(self):
return self.simu.Omega_M return self.simu.Omega_M
def getOmega_Lambda(self): def getOmega_Lambda(self):
return self.simu.Omega_Lambda return self.simu.Omega_Lambda
@ -281,7 +284,7 @@ cdef class ArrayWrapper:
cdef void* data_ptr cdef void* data_ptr
cdef np.uint64_t size cdef np.uint64_t size
cdef int type_array cdef int type_array
cdef set_data(self, np.uint64_t size, int type_array, void* data_ptr): cdef set_data(self, np.uint64_t size, int type_array, void* data_ptr):
""" Set the data of the array """ Set the data of the array
This cannot be done in the constructor as it must recieve C-level This cannot be done in the constructor as it must recieve C-level
@ -294,22 +297,22 @@ cdef class ArrayWrapper:
self.data_ptr = data_ptr self.data_ptr = data_ptr
self.size = size self.size = size
self.type_array = type_array self.type_array = type_array
def __array__(self): def __array__(self):
""" Here we use the __array__ method, that is called when numpy """ Here we use the __array__ method, that is called when numpy
tries to get an array from the object.""" tries to get an array from the object."""
cdef np.npy_intp shape[1] cdef np.npy_intp shape[1]
shape[0] = <np.npy_intp> self.size shape[0] = <np.npy_intp> self.size
# Create a 1D array, of length 'size' # Create a 1D array, of length 'size'
ndarray = np.PyArray_SimpleNewFromData(1, shape, self.type_array, self.data_ptr) ndarray = np.PyArray_SimpleNewFromData(1, shape, self.type_array, self.data_ptr)
return ndarray return ndarray
def __dealloc__(self): def __dealloc__(self):
""" Frees the array. This is called by Python when all the """ Frees the array. This is called by Python when all the
references to the object are gone. """ references to the object are gone. """
pass pass
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
@ -319,7 +322,7 @@ cdef object wrap_array(void *p, np.uint64_t s, int typ):
ndarray = np.array(wrapper, copy=False) ndarray = np.array(wrapper, copy=False)
ndarray.base = <PyObject*> wrapper ndarray.base = <PyObject*> wrapper
Py_INCREF(wrapper) Py_INCREF(wrapper)
return ndarray return ndarray
@ -368,7 +371,7 @@ def loadGadget(str filename, int snapshot_id, int gadgetFormat = 1, bool loadPos
"""loadGadget(filename, snapshot_id, gadgetFormat = 1, loadPosition=True, loadVelocity=False, loadId=False, loadType=False) """loadGadget(filename, snapshot_id, gadgetFormat = 1, loadPosition=True, loadVelocity=False, loadId=False, loadType=False)
This function loads Gadget-1 snapshot format. This function loads Gadget-1 snapshot format.
If snapshot_id is negative then the snapshot is considered not to be part of If snapshot_id is negative then the snapshot is considered not to be part of
a set of snapshots written by different cpu. Otherwise the filename is modified a set of snapshots written by different cpu. Otherwise the filename is modified
to reflect the indicated snapshot_id. to reflect the indicated snapshot_id.
@ -376,16 +379,16 @@ def loadGadget(str filename, int snapshot_id, int gadgetFormat = 1, bool loadPos
Arguments: Arguments:
filename (str): input filename filename (str): input filename
snapshot_id (int): identifier of the gadget file if it is a multi-file snapshot snapshot_id (int): identifier of the gadget file if it is a multi-file snapshot
Keyword arguments: Keyword arguments:
loadPosition (bool): whether to load positions loadPosition (bool): whether to load positions
loadVelocity (bool): whether to load velocities loadVelocity (bool): whether to load velocities
loadId (bool): whether to load unique identifiers loadId (bool): whether to load unique identifiers
loadType (bool): whether to set types to particles loadType (bool): whether to set types to particles
loadMass (bool): whether to set the mass of particles loadMass (bool): whether to set the mass of particles
Returns: Returns:
an PySimulationAdaptor instance. an PySimulationAdaptor instance.
""" """
cdef int flags cdef int flags
@ -419,13 +422,13 @@ def loadParallelGadget(object filename_list, int gadgetFormat = 1, bool loadPosi
Arguments: Arguments:
filename (list): a list or tuple of filenames to load in parallel filename (list): a list or tuple of filenames to load in parallel
Keyword arguments: Keyword arguments:
loadPosition (bool): indicate to load positions loadPosition (bool): indicate to load positions
loadVelocity (bool): indicate to load velocities loadVelocity (bool): indicate to load velocities
loadId (bool): indicate to load id loadId (bool): indicate to load id
loadType (bool): indicate to load particle types loadType (bool): indicate to load particle types
Returns: Returns:
It loads a gadget-1 snapshot and return a cosmotool.PySimulationBase_ object. It loads a gadget-1 snapshot and return a cosmotool.PySimulationBase_ object.
""" """
@ -453,16 +456,16 @@ def loadParallelGadget(object filename_list, int gadgetFormat = 1, bool loadPosi
data = alloc_simudata(num_files) data = alloc_simudata(num_files)
for i,l in enumerate(filename_list): for i,l in enumerate(filename_list):
filenames[i] = l.encode('utf-8') filenames[i] = l.encode('utf-8')
with nogil: with nogil:
for i in prange(num_files): for i in prange(num_files):
local_data = loadGadgetMulti_safe(filenames[i], flags, gadgetFormat) local_data = loadGadgetMulti_safe(filenames[i], flags, gadgetFormat)
data[i] = local_data data[i] = local_data
# data[i] = loadGadgetMulti(filenames[i].c_str(), -1, flags) # data[i] = loadGadgetMulti(filenames[i].c_str(), -1, flags)
out_arrays = [] out_arrays = []
for i in xrange(num_files): for i in xrange(num_files):
if data[i] == <SimuData*>0: if data[i] == <SimuData*>0:
out_arrays.append(None) out_arrays.append(None)
else: else:
out_arrays.append(PySimulationAdaptor(wrap_simudata(data[i], flags))) out_arrays.append(PySimulationAdaptor(wrap_simudata(data[i], flags)))
@ -473,10 +476,10 @@ def loadParallelGadget(object filename_list, int gadgetFormat = 1, bool loadPosi
def writeGadget(str filename, object simulation): def writeGadget(str filename, object simulation):
"""writeGadget(filename, simulation) """writeGadget(filename, simulation)
This function attempts to write the content of the simulation object into This function attempts to write the content of the simulation object into
a file named `filename` using a Gadget-1 file format. a file named `filename` using a Gadget-1 file format.
Arguments: Arguments:
filename (str): output filename filename (str): output filename
simulation (PySimulationBase): a simulation object simulation (PySimulationBase): a simulation object
@ -486,23 +489,23 @@ def writeGadget(str filename, object simulation):
cdef np.ndarray[np.int64_t, ndim=1] ids cdef np.ndarray[np.int64_t, ndim=1] ids
cdef np.int64_t NumPart cdef np.int64_t NumPart
cdef int j cdef int j
if not isinstance(simulation,PySimulationBase): if not isinstance(simulation,PySimulationBase):
raise TypeError("Second argument must be of type SimulationBase") raise TypeError("Second argument must be of type SimulationBase")
NumPart = simulation.positions[0].size NumPart = simulation.positions[0].size
simdata.noAuto = True simdata.noAuto = True
for j in xrange(3): for j in xrange(3):
pos = simulation.getPositions()[j] pos = simulation.getPositions()[j]
vel = simulation.getVelocities()[j] vel = simulation.getVelocities()[j]
if pos.size != NumPart or vel.size != NumPart: if pos.size != NumPart or vel.size != NumPart:
raise ValueError("Invalid number of particles") raise ValueError("Invalid number of particles")
simdata.Pos[j] = <float *>pos.data simdata.Pos[j] = <float *>pos.data
simdata.Vel[j] = <float *>vel.data simdata.Vel[j] = <float *>vel.data
ids = simulation.getIdentifiers() ids = simulation.getIdentifiers()
simdata.Id = <int64_t *>ids.data simdata.Id = <int64_t *>ids.data
simdata.BoxSize = simulation.getBoxsize() simdata.BoxSize = simulation.getBoxsize()
@ -519,21 +522,21 @@ def writeGadget(str filename, object simulation):
def loadRamses(str basepath, int snapshot_id, int cpu_id, bool doublePrecision = False, bool loadPosition = True, bool loadVelocity = False, bool loadId = False, bool loadMass = False): def loadRamses(str basepath, int snapshot_id, int cpu_id, bool doublePrecision = False, bool loadPosition = True, bool loadVelocity = False, bool loadId = False, bool loadMass = False):
""" loadRamses(basepath, snapshot_id, cpu_id, doublePrecision = False, loadPosition = True, loadVelocity = False) """ loadRamses(basepath, snapshot_id, cpu_id, doublePrecision = False, loadPosition = True, loadVelocity = False)
Loads the indicated snapshot based on the cpu id, snapshot id and basepath. It is important to specify the correct precision in doublePrecision otherwise the loading will fail. There is no way of auto-detecting properly the precision of the snapshot file. Loads the indicated snapshot based on the cpu id, snapshot id and basepath. It is important to specify the correct precision in doublePrecision otherwise the loading will fail. There is no way of auto-detecting properly the precision of the snapshot file.
Args: Args:
basepath (str): the base directory of the snapshot basepath (str): the base directory of the snapshot
snapshot_id (int): the snapshot id snapshot_id (int): the snapshot id
cpu_id (int): the cpu id of the file to load cpu_id (int): the cpu id of the file to load
Keyword args: Keyword args:
doublePrecision (bool): By default it is False, thus singlePrecision doublePrecision (bool): By default it is False, thus singlePrecision
loadPosition (bool): Whether to load positions loadPosition (bool): Whether to load positions
loadVelocity (bool): Whether to load velocities loadVelocity (bool): Whether to load velocities
loadId (bool): Whether to load identifiers loadId (bool): Whether to load identifiers
loadMass (bool): Whether to load mass value loadMass (bool): Whether to load mass value
Returns: Returns:
An object derived from PySimulationBase_. An object derived from PySimulationBase_.
""" """
cdef int flags cdef int flags
cdef SimuData *data cdef SimuData *data
@ -549,7 +552,7 @@ def loadRamses(str basepath, int snapshot_id, int cpu_id, bool doublePrecision =
if loadMass: if loadMass:
flags |= NEED_MASS flags |= NEED_MASS
encpath = basepath.encode('utf-8') encpath = basepath.encode('utf-8')
try: try:
data = loadRamsesSimu(encpath, snapshot_id, cpu_id, doublePrecision, flags) data = loadRamsesSimu(encpath, snapshot_id, cpu_id, doublePrecision, flags)
if data == <SimuData*>0: if data == <SimuData*>0:
@ -558,4 +561,4 @@ def loadRamses(str basepath, int snapshot_id, int cpu_id, bool doublePrecision =
raise RuntimeError(str(e) + ' (check the float precision in snapshot)') raise RuntimeError(str(e) + ' (check the float precision in snapshot)')
return PySimulationAdaptor(wrap_simudata(data, flags)) return PySimulationAdaptor(wrap_simudata(data, flags))

143
python/_project.pyx
View File

@ -25,8 +25,8 @@ cdef extern from "openmp.hpp" namespace "CosmoTool":
@cython.boundscheck(False) @cython.boundscheck(False)
@cython.cdivision(True) @cython.cdivision(True)
@cython.wraparound(False) @cython.wraparound(False)
cdef void interp3d_INTERNAL_periodic(DTYPE_t x, DTYPE_t y, cdef void interp3d_INTERNAL_periodic(DTYPE_t x, DTYPE_t y,
DTYPE_t z, DTYPE_t z,
DTYPE_t[:,:,:] d, DTYPE_t Lbox, DTYPE_t *retval) nogil: DTYPE_t[:,:,:] d, DTYPE_t Lbox, DTYPE_t *retval) nogil:
cdef int Ngrid = d.shape[0] cdef int Ngrid = d.shape[0]
@ -84,8 +84,8 @@ cdef void interp3d_INTERNAL_periodic(DTYPE_t x, DTYPE_t y,
@cython.boundscheck(False) @cython.boundscheck(False)
@cython.cdivision(True) @cython.cdivision(True)
@cython.wraparound(False) @cython.wraparound(False)
cdef void ngp3d_INTERNAL_periodic(DTYPE_t x, DTYPE_t y, cdef void ngp3d_INTERNAL_periodic(DTYPE_t x, DTYPE_t y,
DTYPE_t z, DTYPE_t z,
DTYPE_t[:,:,:] d, DTYPE_t Lbox, DTYPE_t *retval) nogil: DTYPE_t[:,:,:] d, DTYPE_t Lbox, DTYPE_t *retval) nogil:
cdef int Ngrid = d.shape[0] cdef int Ngrid = d.shape[0]
@ -108,14 +108,14 @@ cdef void ngp3d_INTERNAL_periodic(DTYPE_t x, DTYPE_t y,
iy = iy%Ngrid iy = iy%Ngrid
iz = iz%Ngrid iz = iz%Ngrid
retval[0] = d[ix ,iy ,iz ] retval[0] = d[ix ,iy ,iz ]
@cython.boundscheck(False) @cython.boundscheck(False)
@cython.cdivision(True) @cython.cdivision(True)
@cython.wraparound(False) @cython.wraparound(False)
cdef void ngp3d_INTERNAL(DTYPE_t x, DTYPE_t y, cdef void ngp3d_INTERNAL(DTYPE_t x, DTYPE_t y,
DTYPE_t z, DTYPE_t z,
DTYPE_t[:,:,:] d, DTYPE_t Lbox, DTYPE_t *retval, DTYPE_t inval) nogil: DTYPE_t[:,:,:] d, DTYPE_t Lbox, DTYPE_t *retval, DTYPE_t inval) nogil:
cdef int Ngrid = d.shape[0] cdef int Ngrid = d.shape[0]
@ -137,16 +137,16 @@ cdef void ngp3d_INTERNAL(DTYPE_t x, DTYPE_t y,
retval[0] = inval retval[0] = inval
return return
retval[0] = d[ix ,iy ,iz ] retval[0] = d[ix ,iy ,iz ]
@cython.boundscheck(False) @cython.boundscheck(False)
@cython.cdivision(True) @cython.cdivision(True)
@cython.wraparound(False) @cython.wraparound(False)
cdef void interp3d_INTERNAL(DTYPE_t x, DTYPE_t y, cdef void interp3d_INTERNAL(DTYPE_t x, DTYPE_t y,
DTYPE_t z, DTYPE_t z,
DTYPE_t[:,:,:] d, DTYPE_t Lbox, DTYPE_t *retval, DTYPE_t inval) nogil: DTYPE_t[:,:,:] d, DTYPE_t Lbox, DTYPE_t *retval, DTYPE_t inval) nogil:
cdef int Ngrid = d.shape[0] cdef int Ngrid = d.shape[0]
cdef DTYPE_t inv_delta = Ngrid/Lbox cdef DTYPE_t inv_delta = Ngrid/Lbox
cdef int ix, iy, iz cdef int ix, iy, iz
@ -193,13 +193,13 @@ cdef void interp3d_INTERNAL(DTYPE_t x, DTYPE_t y,
d[ix+1,iy+1,iz+1] * f[1][1][1] d[ix+1,iy+1,iz+1] * f[1][1][1]
@cython.boundscheck(False) @cython.boundscheck(False)
def interp3d(x not None, y not None, def interp3d(x not None, y not None,
z not None, z not None,
npx.ndarray[DTYPE_t, ndim=3] d not None, DTYPE_t Lbox, npx.ndarray[DTYPE_t, ndim=3] d not None, DTYPE_t Lbox,
bool periodic=False, bool centered=True, bool ngp=False, DTYPE_t inval = 0): bool periodic=False, bool centered=True, bool ngp=False, DTYPE_t inval = 0):
""" interp3d(x,y,z,d,Lbox,periodic=False,centered=True,ngp=False) -> interpolated values """ interp3d(x,y,z,d,Lbox,periodic=False,centered=True,ngp=False) -> interpolated values
Compute the tri-linear interpolation of the given field (d) at the given position (x,y,z). It assumes that they are box-centered coordinates. So (x,y,z) == (0,0,0) is equivalent to the pixel at (Nx/2,Ny/2,Nz/2) with Nx,Ny,Nz = d.shape. If periodic is set, it assumes the box is periodic Compute the tri-linear interpolation of the given field (d) at the given position (x,y,z). It assumes that they are box-centered coordinates. So (x,y,z) == (0,0,0) is equivalent to the pixel at (Nx/2,Ny/2,Nz/2) with Nx,Ny,Nz = d.shape. If periodic is set, it assumes the box is periodic
""" """
cdef npx.ndarray[DTYPE_t] out cdef npx.ndarray[DTYPE_t] out
cdef DTYPE_t[:] out_slice cdef DTYPE_t[:] out_slice
@ -227,12 +227,12 @@ def interp3d(x not None, y not None,
if type(x) == np.ndarray or type(y) == np.ndarray or type(z) == np.ndarray: if type(x) == np.ndarray or type(y) == np.ndarray or type(z) == np.ndarray:
if type(x) != np.ndarray or type(y) != np.ndarray or type(z) != np.ndarray: if type(x) != np.ndarray or type(y) != np.ndarray or type(z) != np.ndarray:
raise ValueError("All or no array. No partial arguments") raise ValueError("All or no array. No partial arguments")
ax = x ax = x
ay = y ay = y
az = z az = z
assert ax.size == ay.size and ax.size == az.size assert ax.size == ay.size and ax.size == az.size
out = np.empty(x.shape, dtype=DTYPE) out = np.empty(x.shape, dtype=DTYPE)
out_slice = out out_slice = out
in_slice = d in_slice = d
@ -280,10 +280,10 @@ cdef DTYPE_t interp2d_INTERNAL_periodic(DTYPE_t x, DTYPE_t y,
rx = (inv_delta*x + Ngrid/2) rx = (inv_delta*x + Ngrid/2)
ry = (inv_delta*y + Ngrid/2) ry = (inv_delta*y + Ngrid/2)
ix = int(floor(rx)) ix = int(floor(rx))
iy = int(floor(ry)) iy = int(floor(ry))
rx -= ix rx -= ix
ry -= iy ry -= iy
@ -291,13 +291,13 @@ cdef DTYPE_t interp2d_INTERNAL_periodic(DTYPE_t x, DTYPE_t y,
ix += Ngrid ix += Ngrid
while iy < 0: while iy < 0:
iy += Ngrid iy += Ngrid
jx = (ix+1)%Ngrid jx = (ix+1)%Ngrid
jy = (iy+1)%Ngrid jy = (iy+1)%Ngrid
assert ((ix >= 0) and ((jx) < Ngrid)) assert ((ix >= 0) and ((jx) < Ngrid))
assert ((iy >= 0) and ((jy) < Ngrid)) assert ((iy >= 0) and ((jy) < Ngrid))
f[0][0] = (1-rx)*(1-ry) f[0][0] = (1-rx)*(1-ry)
f[1][0] = ( rx)*(1-ry) f[1][0] = ( rx)*(1-ry)
f[0][1] = (1-rx)*( ry) f[0][1] = (1-rx)*( ry)
@ -314,7 +314,7 @@ cdef DTYPE_t interp2d_INTERNAL_periodic(DTYPE_t x, DTYPE_t y,
@cython.cdivision(True) @cython.cdivision(True)
cdef DTYPE_t interp2d_INTERNAL(DTYPE_t x, DTYPE_t y, cdef DTYPE_t interp2d_INTERNAL(DTYPE_t x, DTYPE_t y,
npx.ndarray[DTYPE_t, ndim=2] d, DTYPE_t Lbox) except? 0: npx.ndarray[DTYPE_t, ndim=2] d, DTYPE_t Lbox) except? 0:
cdef int Ngrid = d.shape[0] cdef int Ngrid = d.shape[0]
cdef DTYPE_t inv_delta = Ngrid/Lbox cdef DTYPE_t inv_delta = Ngrid/Lbox
cdef int ix, iy cdef int ix, iy
@ -348,7 +348,7 @@ cdef DTYPE_t interp2d_INTERNAL(DTYPE_t x, DTYPE_t y,
d[ix+1,iy ] * f[1][0] + \ d[ix+1,iy ] * f[1][0] + \
d[ix ,iy+1] * f[0][1] + \ d[ix ,iy+1] * f[0][1] + \
d[ix+1,iy+1] * f[1][1] d[ix+1,iy+1] * f[1][1]
def interp2d(x not None, y not None, def interp2d(x not None, y not None,
npx.ndarray[DTYPE_t, ndim=2] d not None, DTYPE_t Lbox, npx.ndarray[DTYPE_t, ndim=2] d not None, DTYPE_t Lbox,
bool periodic=False): bool periodic=False):
@ -362,11 +362,11 @@ def interp2d(x not None, y not None,
if type(x) == np.ndarray or type(y) == np.ndarray: if type(x) == np.ndarray or type(y) == np.ndarray:
if type(x) != np.ndarray or type(y) != np.ndarray: if type(x) != np.ndarray or type(y) != np.ndarray:
raise ValueError("All or no array. No partial arguments") raise ValueError("All or no array. No partial arguments")
ax = x ax = x
ay = y ay = y
assert ax.size == ay.size assert ax.size == ay.size
out = np.empty(x.shape, dtype=DTYPE) out = np.empty(x.shape, dtype=DTYPE)
if periodic: if periodic:
for i in range(ax.size): for i in range(ax.size):
@ -381,8 +381,8 @@ def interp2d(x not None, y not None,
return interp2d_INTERNAL_periodic(x, y, d, Lbox) return interp2d_INTERNAL_periodic(x, y, d, Lbox)
else: else:
return interp2d_INTERNAL(x, y, d, Lbox) return interp2d_INTERNAL(x, y, d, Lbox)
@cython.boundscheck(False) @cython.boundscheck(False)
@cython.cdivision(True) @cython.cdivision(True)
cdef void INTERNAL_project_cic_no_mass(DTYPE_t[:,:,:] g, cdef void INTERNAL_project_cic_no_mass(DTYPE_t[:,:,:] g,
@ -450,7 +450,7 @@ cdef void INTERNAL_project_cic_no_mass_periodic(DTYPE_t[:,:,:] g,
ag[b1[0],b[1],b[2]] += a[0]*c[1]*c[2] ag[b1[0],b[1],b[2]] += a[0]*c[1]*c[2]
ag[b[0],b1[1],b[2]] += c[0]*a[1]*c[2] ag[b[0],b1[1],b[2]] += c[0]*a[1]*c[2]
ag[b1[0],b1[1],b[2]] += a[0]*a[1]*c[2] ag[b1[0],b1[1],b[2]] += a[0]*a[1]*c[2]
ag[b[0],b[1],b1[2]] += c[0]*c[1]*a[2] ag[b[0],b[1],b1[2]] += c[0]*c[1]*a[2]
ag[b1[0],b[1],b1[2]] += a[0]*c[1]*a[2] ag[b1[0],b[1],b1[2]] += a[0]*c[1]*a[2]
ag[b[0],b1[1],b1[2]] += c[0]*a[1]*a[2] ag[b[0],b1[1],b1[2]] += c[0]*a[1]*a[2]
@ -525,20 +525,21 @@ cdef void INTERNAL_project_cic_with_mass_periodic(DTYPE_t[:,:,:] g,
g[b1[0],b[1],b[2]] += a[0]*c[1]*c[2]*m0 g[b1[0],b[1],b[2]] += a[0]*c[1]*c[2]*m0
g[b[0],b1[1],b[2]] += c[0]*a[1]*c[2]*m0 g[b[0],b1[1],b[2]] += c[0]*a[1]*c[2]*m0
g[b1[0],b1[1],b[2]] += a[0]*a[1]*c[2]*m0 g[b1[0],b1[1],b[2]] += a[0]*a[1]*c[2]*m0
g[b[0],b[1],b1[2]] += c[0]*c[1]*a[2]*m0 g[b[0],b[1],b1[2]] += c[0]*c[1]*a[2]*m0
g[b1[0],b[1],b1[2]] += a[0]*c[1]*a[2]*m0 g[b1[0],b[1],b1[2]] += a[0]*c[1]*a[2]*m0
g[b[0],b1[1],b1[2]] += c[0]*a[1]*a[2]*m0 g[b[0],b1[1],b1[2]] += c[0]*a[1]*a[2]*m0
g[b1[0],b1[1],b1[2]] += a[0]*a[1]*a[2]*m0 g[b1[0],b1[1],b1[2]] += a[0]*a[1]*a[2]*m0
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):
"""
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. 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):
"""
project_cic(x array (N,3), mass (may be None), Ngrid, Lbox, periodict, centered=True, output=None)
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
@ -558,7 +559,13 @@ 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")
g = np.zeros((Ngrid,Ngrid,Ngrid),dtype=DTYPE) if output is None:
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
@ -569,7 +576,7 @@ def project_cic(npx.ndarray[DTYPE_t, ndim=2] x not None, npx.ndarray[DTYPE_t, nd
else: else:
d_mass = mass d_mass = mass
with nogil: with nogil:
INTERNAL_project_cic_with_mass(d_g, d_x, d_mass, Ngrid, Lbox, shifter) INTERNAL_project_cic_with_mass(d_g, d_x, d_mass, Ngrid, Lbox, shifter)
else: else:
if mass is None: if mass is None:
with nogil: with nogil:
@ -577,13 +584,13 @@ def project_cic(npx.ndarray[DTYPE_t, ndim=2] x not None, npx.ndarray[DTYPE_t, nd
else: else:
d_mass = mass d_mass = mass
with nogil: with nogil:
INTERNAL_project_cic_with_mass_periodic(d_g, d_x, d_mass, Ngrid, Lbox, shifter) INTERNAL_project_cic_with_mass_periodic(d_g, d_x, d_mass, Ngrid, Lbox, shifter)
return g return g
def tophat_fourier_internal(npx.ndarray[DTYPE_t, ndim=1] x not None): def tophat_fourier_internal(npx.ndarray[DTYPE_t, ndim=1] x not None):
cdef int i cdef int i
cdef npx.ndarray[DTYPE_t] y cdef npx.ndarray[DTYPE_t] y
cdef DTYPE_t x0 cdef DTYPE_t x0
y = np.empty(x.size, dtype=DTYPE) y = np.empty(x.size, dtype=DTYPE)
@ -609,7 +616,7 @@ def tophat_fourier(x not None):
return b.reshape(x.shape) return b.reshape(x.shape)
@cython.boundscheck(False) @cython.boundscheck(False)
@cython.cdivision(True) @cython.cdivision(True)
@ -659,25 +666,25 @@ cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t
if tmp_a < alpha_max: if tmp_a < alpha_max:
alpha_max = tmp_a alpha_max = tmp_a
j = i j = i
I = compute_projection(vertex_value, u, u0, alpha_max) I = compute_projection(vertex_value, u, u0, alpha_max)
for i in xrange(3): for i in xrange(3):
u0[i] += u[i]*alpha_max u0[i] += u[i]*alpha_max
# alpha_max is the integration length # alpha_max is the integration length
# we integrate between 0 and alpha_max (curvilinear coordinates) # we integrate between 0 and alpha_max (curvilinear coordinates)
r[0] = j r[0] = j
return I return I
@cython.boundscheck(False) @cython.boundscheck(False)
cdef DTYPE_t integrator0(DTYPE_t[:,:,:] density, cdef DTYPE_t integrator0(DTYPE_t[:,:,:] density,
DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1], DTYPE_t alpha_max) nogil: DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1], DTYPE_t alpha_max) nogil:
cdef DTYPE_t d cdef DTYPE_t d
d = density[iu0[0], iu0[1], iu0[2]] d = density[iu0[0], iu0[1], iu0[2]]
return cube_integral(u, u0, jumper, alpha_max)*d return cube_integral(u, u0, jumper, alpha_max)*d
@cython.boundscheck(False) @cython.boundscheck(False)
@ -687,7 +694,7 @@ cdef DTYPE_t integrator1(DTYPE_t[:,:,:] density,
cdef DTYPE_t d cdef DTYPE_t d
cdef int a[3][2] cdef int a[3][2]
cdef int i cdef int i
for i in xrange(3): for i in xrange(3):
a[i][0] = iu0[i] a[i][0] = iu0[i]
a[i][1] = iu0[i]+1 a[i][1] = iu0[i]+1
@ -705,14 +712,14 @@ cdef DTYPE_t integrator1(DTYPE_t[:,:,:] density,
return cube_integral_trilin(u, u0, jumper, vertex_value, alpha_max) return cube_integral_trilin(u, u0, jumper, vertex_value, alpha_max)
@cython.boundscheck(False) @cython.boundscheck(False)
cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density, cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density,
DTYPE_t a_u[3], DTYPE_t a_u[3],
DTYPE_t min_distance, DTYPE_t min_distance,
DTYPE_t max_distance, DTYPE_t[:] shifter, int integrator_id) nogil except? 0: DTYPE_t max_distance, DTYPE_t[:] shifter, int integrator_id) nogil except? 0:
cdef DTYPE_t u[3] cdef DTYPE_t u[3]
cdef DTYPE_t ifu0[3] cdef DTYPE_t ifu0[3]
cdef DTYPE_t u0[3] cdef DTYPE_t u0[3]
cdef DTYPE_t utot[3] cdef DTYPE_t utot[3]
@ -724,7 +731,7 @@ cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density,
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]
cdef DTYPE_t (*integrator)(DTYPE_t[:,:,:], cdef DTYPE_t (*integrator)(DTYPE_t[:,:,:],
DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1], DTYPE_t alpha_max) nogil DTYPE_t u[3], DTYPE_t u0[3], int u_delta[3], int iu0[3], int jumper[1], DTYPE_t alpha_max) nogil
@ -747,7 +754,7 @@ cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density,
if (not ((iu0[i]>= 0) and (iu0[i] < N))): if (not ((iu0[i]>= 0) and (iu0[i] < N))):
with gil: with gil:
raise RuntimeError("iu0[%d] = %d !!" % (i,iu0[i])) raise RuntimeError("iu0[%d] = %d !!" % (i,iu0[i]))
if (not (u0[i]>=0 and u0[i]<=1)): if (not (u0[i]>=0 and u0[i]<=1)):
with gil: with gil:
raise RuntimeError("u0[%d] = %g !" % (i,u0[i])) raise RuntimeError("u0[%d] = %g !" % (i,u0[i]))
@ -756,7 +763,7 @@ cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density,
if ((iu0[0] >= N-1) or (iu0[0] <= 0) or if ((iu0[0] >= N-1) or (iu0[0] <= 0) or
(iu0[1] >= N-1) or (iu0[1] <= 0) or (iu0[1] >= N-1) or (iu0[1] <= 0) or
(iu0[2] >= N-1) or (iu0[2] <= 0)): (iu0[2] >= N-1) or (iu0[2] <= 0)):
completed = 1 completed = 1
I0 = 0 I0 = 0
jumper[0] = 0 jumper[0] = 0
@ -771,8 +778,8 @@ cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density,
iu0[jumper[0]] += 1 iu0[jumper[0]] += 1
u0[jumper[0]] = 0 u0[jumper[0]] = 0
if ((iu0[0] >= N-1) or (iu0[0] <= 0) or if ((iu0[0] >= N-1) or (iu0[0] <= 0) or
(iu0[1] >= N-1) or (iu0[1] <= 0) or (iu0[1] >= N-1) or (iu0[1] <= 0) or
(iu0[2] >= N-1) or (iu0[2] <= 0)): (iu0[2] >= N-1) or (iu0[2] <= 0)):
completed = 1 completed = 1
@ -787,7 +794,7 @@ cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density,
#delta = sqrt(dist2) - max_distance #delta = sqrt(dist2) - max_distance
#I0 -= d*delta #I0 -= d*delta
completed = 1 completed = 1
return I0 return I0
def line_of_sight_projection(DTYPE_t[:,:,:] density not None, def line_of_sight_projection(DTYPE_t[:,:,:] density not None,
@ -795,18 +802,18 @@ def line_of_sight_projection(DTYPE_t[:,:,:] density not None,
DTYPE_t min_distance, DTYPE_t min_distance,
DTYPE_t max_distance, DTYPE_t[:] shifter not None, int integrator_id=0): DTYPE_t max_distance, DTYPE_t[:] shifter not None, int integrator_id=0):
cdef DTYPE_t u[3] cdef DTYPE_t u[3]
u[0] = a_u[0] u[0] = a_u[0]
u[1] = a_u[1] u[1] = a_u[1]
u[2] = a_u[2] u[2] = a_u[2]
return C_line_of_sight_projection(density, return C_line_of_sight_projection(density,
u, u,
min_distance, min_distance,
max_distance, shifter, integrator_id) max_distance, shifter, integrator_id)
cdef double _spherical_projloop(double theta, double phi, DTYPE_t[:,:,:] density, cdef double _spherical_projloop(double theta, double phi, DTYPE_t[:,:,:] density,
double min_distance, double max_distance, double min_distance, double max_distance,
DTYPE_t[:] shifter, int integrator_id) nogil: DTYPE_t[:] shifter, int integrator_id) nogil:
cdef DTYPE_t u0[3] cdef DTYPE_t u0[3]
@ -814,7 +821,7 @@ cdef double _spherical_projloop(double theta, double phi, DTYPE_t[:,:,:] density
u0[0] = cos(phi)*stheta u0[0] = cos(phi)*stheta
u0[1] = sin(phi)*stheta u0[1] = sin(phi)*stheta
u0[2] = cos(theta) u0[2] = cos(theta)
return C_line_of_sight_projection(density, u0, min_distance, max_distance, shifter, integrator_id) return C_line_of_sight_projection(density, u0, min_distance, max_distance, shifter, integrator_id)
@ -825,20 +832,20 @@ def spherical_projection(int Nside,
DTYPE_t max_distance, int progress=1, int integrator_id=0, DTYPE_t[:] shifter = None, int booster=-1): 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) spherical_projection(Nside, density, min_distance, max_distance, progress=1, integrator_id=0, shifter=None, booster=-1)
Keyword arguments: Keyword arguments:
progress (int): show progress if it is equal to 1 progress (int): show progress if it is equal to 1
integrator_id (int): specify the order of integration along the line of shift 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 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. booster (int): what is the frequency of refreshment of the progress bar. Small number decreases performance by locking the GIL.
Arguments: Arguments:
Nside (int): Nside of the returned map Nside (int): Nside of the returned map
density (NxNxN array): this is the density field, expressed as a cubic array density (NxNxN array): this is the density field, expressed as a cubic array
min_distance (float): lower bound of the integration min_distance (float): lower bound of the integration
max_distance (float): upper bound of the integration max_distance (float): upper bound of the integration
Returns: Returns:
an healpix map, as a 1-dimensional array. an healpix map, as a 1-dimensional array.
""" """
@ -853,11 +860,11 @@ def spherical_projection(int Nside,
cdef long N, N0 cdef long N, N0
cdef double stheta cdef double stheta
cdef int tid cdef int tid
if shifter is None: if shifter is None:
shifter = view.array(shape=(3,), format=FORMAT_DTYPE, itemsize=sizeof(DTYPE_t)) shifter = view.array(shape=(3,), format=FORMAT_DTYPE, itemsize=sizeof(DTYPE_t))
shifter[:] = 0 shifter[:] = 0
print("allocating map") print("allocating map")
outm_array = np.empty(hp.nside2npix(Nside),dtype=DTYPE) outm_array = np.empty(hp.nside2npix(Nside),dtype=DTYPE)
print("initializing views") print("initializing views")
@ -870,10 +877,10 @@ def spherical_projection(int Nside,
N = smp_get_max_threads() N = smp_get_max_threads()
N0 = outm.size N0 = outm.size
if booster < 0: if booster < 0:
booster = 1#000 booster = 1#000
job_done = view.array(shape=(N,), format="i", itemsize=sizeof(int)) job_done = view.array(shape=(N,), format="i", itemsize=sizeof(int))
job_done[:] = 0 job_done[:] = 0
theta,phi = hp.pix2ang(Nside, np.arange(N0)) theta,phi = hp.pix2ang(Nside, np.arange(N0))

6
sample/simple3DFilter.cpp
View File

@ -108,7 +108,7 @@ int main(int argc, char **argv)
H5::H5File in_f(fname1, 0); H5::H5File in_f(fname1, 0);
H5::H5File out_f(outFile, H5F_ACC_TRUNC); H5::H5File out_f(outFile, H5F_ACC_TRUNC);
array_type v1_data; array_type v1_data;
uint32_t N1_points, N2_points; uint64_t N1_points, N2_points;
array3_type bins(boost::extents[Nres][Nres][Nres]); array3_type bins(boost::extents[Nres][Nres][Nres]);
@ -124,12 +124,12 @@ 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 (long i = 0; i < Nres*Nres*Nres; i++) for (uint32_t 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 (int i = 0 ; i < N1_points; i++) for (uint64_t 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];

2
setup.py
View File

@ -229,7 +229,7 @@ class BuildCMakeExt(build_ext):
CosmoTool_extension = CMakeExtension(name="cosmotool") CosmoTool_extension = CMakeExtension(name="cosmotool")
setup(name='cosmotool', setup(name='cosmotool',
version='1.2.3', version='1.3.0',
packages=["cosmotool"], packages=["cosmotool"],
package_dir={'cosmotool': 'python/cosmotool'}, package_dir={'cosmotool': 'python/cosmotool'},
install_requires=['numpy','cffi','numexpr','pyfftw','h5py'], install_requires=['numpy','cffi','numexpr','pyfftw','h5py'],

20
src/mykdtree.hpp
View File

@ -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;
uint32_t currentRank; uint64_t currentRank;
uint32_t numCells; uint64_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,
uint32_t numCells) uint64_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);
} }
uint32_t getIntersection(const coords& x, CoordType r, uint64_t getIntersection(const coords& x, CoordType r,
Cell **cells, Cell **cells,
uint32_t numCells); uint64_t numCells);
uint32_t getIntersection(const coords& x, CoordType r, uint64_t getIntersection(const coords& x, CoordType r,
Cell **cells, Cell **cells,
CoordType *distances, CoordType *distances,
uint32_t numCells); uint64_t numCells);
uint32_t countCells(const coords& x, CoordType r); uint64_t countCells(const coords& x, CoordType r);
Cell *getNearestNeighbour(const coords& x); Cell *getNearestNeighbour(const coords& x);
void getNearestNeighbours(const coords& x, uint32_t NumCells, void getNearestNeighbours(const coords& x, uint64_t NumCells,
Cell **cells); Cell **cells);
void getNearestNeighbours(const coords& x, uint32_t NumCells, void getNearestNeighbours(const coords& x, uint64_t NumCells,
Cell **cells, Cell **cells,
CoordType *distances); CoordType *distances);

14
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>
uint32_t KDTree<N,ValType,CType,CellSplitter>::getIntersection(const coords& x, CoordType r, uint64_t KDTree<N,ValType,CType,CellSplitter>::getIntersection(const coords& x, CoordType r,
KDTree<N,ValType,CType,CellSplitter>::Cell **cells, KDTree<N,ValType,CType,CellSplitter>::Cell **cells,
uint32_t numCells) uint64_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>
uint32_t KDTree<N,ValType,CType,CellSplitter>::getIntersection(const coords& x, CoordType r, uint64_t KDTree<N,ValType,CType,CellSplitter>::getIntersection(const coords& x, CoordType r,
Cell **cells, Cell **cells,
CoordType *distances, CoordType *distances,
uint32_t numCells) uint64_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>
uint32_t KDTree<N,ValType,CType,CellSplitter>::countCells(const coords& x, CoordType r) uint64_t KDTree<N,ValType,CType,CellSplitter>::countCells(const coords& x, CoordType r)
{ {
RecursionInfoCells<N,ValType> info; 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, uint32_t N2, void KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbours(const coords& x, uint64_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, uint32_t N2, void KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbours(const coords& x, uint64_t N2,
Cell **cells, Cell **cells,
CoordType *distances) CoordType *distances)
{ {

55
src/yorick.hpp
View File

@ -7,16 +7,16 @@ This software is a computer program whose purpose is to provide a toolbox for co
data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
This software is governed by the CeCILL license under French law and This software is governed by the CeCILL license under French law and
abiding by the rules of distribution of free software. You can use, abiding by the rules of distribution of free software. You can use,
modify and/ or redistribute the software under the terms of the CeCILL modify and/ or redistribute the software under the terms of the CeCILL
license as circulated by CEA, CNRS and INRIA at the following URL license as circulated by CEA, CNRS and INRIA at the following URL
"http://www.cecill.info". "http://www.cecill.info".
As a counterpart to the access to the source code and rights to copy, As a counterpart to the access to the source code and rights to copy,
modify and redistribute granted by the license, users are provided only modify and redistribute granted by the license, users are provided only
with a limited warranty and the software's author, the holder of the with a limited warranty and the software's author, the holder of the
economic rights, and the successive licensors have only limited economic rights, and the successive licensors have only limited
liability. liability.
In this respect, the user's attention is drawn to the risks associated In this respect, the user's attention is drawn to the risks associated
with loading, using, modifying and/or developing or reproducing the with loading, using, modifying and/or developing or reproducing the
@ -25,9 +25,9 @@ that may mean that it is complicated to manipulate, and that also
therefore means that it is reserved for developers and experienced therefore means that it is reserved for developers and experienced
professionals having in-depth computer knowledge. Users are therefore professionals having in-depth computer knowledge. Users are therefore
encouraged to load and test the software's suitability as regards their encouraged to load and test the software's suitability as regards their
requirements in conditions enabling the security of their systems and/or requirements in conditions enabling the security of their systems and/or
data to be ensured and, more generally, to use and operate it in the data to be ensured and, more generally, to use and operate it in the
same conditions as regards security. same conditions as regards security.
The fact that you are presently reading this means that you have had The fact that you are presently reading this means that you have had
knowledge of the CeCILL license and that you accept its terms. knowledge of the CeCILL license and that you accept its terms.
@ -42,7 +42,7 @@ knowledge of the CeCILL license and that you accept its terms.
#include <string> #include <string>
namespace CosmoTool namespace CosmoTool
{ {
class ProgressiveDoubleOutputImpl class ProgressiveDoubleOutputImpl
@ -75,7 +75,7 @@ namespace CosmoTool
bool initialized; bool initialized;
int *ref; int *ref;
ProgressiveDoubleOutputImpl *impl; ProgressiveDoubleOutputImpl *impl;
friend ProgressiveDoubleOutput saveDoubleArrayProgressive(const char *fname, uint32_t *dimList, uint32_t rank); friend ProgressiveDoubleOutput saveDoubleArrayProgressive(const char *fname, uint32_t *dimList, uint32_t rank);
void decRef(); void decRef();
public: public:
@ -86,7 +86,7 @@ namespace CosmoTool
virtual void addDouble(double a); virtual void addDouble(double a);
const ProgressiveDoubleOutput& operator=(const ProgressiveDoubleOutput& b); const ProgressiveDoubleOutput& operator=(const ProgressiveDoubleOutput& b);
}; };
template<class T> template<class T>
@ -95,12 +95,12 @@ namespace CosmoTool
private: private:
int *ref; int *ref;
ProgressiveInputImpl<T> *impl; ProgressiveInputImpl<T> *impl;
void decRef() void decRef()
{ {
if (ref == 0) if (ref == 0)
return; return;
(*ref)--; (*ref)--;
if (*ref == 0) if (*ref == 0)
{ {
@ -112,17 +112,17 @@ namespace CosmoTool
} }
public: public:
static ProgressiveInput<T> static ProgressiveInput<T>
loadArrayProgressive(const std::string& fname, uint32_t *&dimList, loadArrayProgressive(const std::string& fname, uint32_t *&dimList,
uint32_t& rank); uint32_t& rank);
ProgressiveInput() { ProgressiveInput() {
impl = 0; impl = 0;
ref = 0; ref = 0;
} }
ProgressiveInput(ProgressiveInputImpl<T> *i) { ProgressiveInput(ProgressiveInputImpl<T> *i) {
impl = i; impl = i;
ref = new int; ref = new int;
*ref = 1; *ref = 1;
} }
ProgressiveInput(const ProgressiveInput<T>& o) { ProgressiveInput(const ProgressiveInput<T>& o) {
@ -161,12 +161,12 @@ namespace CosmoTool
private: private:
int *ref; int *ref;
ProgressiveOutputImpl<T> *impl; ProgressiveOutputImpl<T> *impl;
void decRef() void decRef()
{ {
if (ref == 0) if (ref == 0)
return; return;
(*ref)--; (*ref)--;
if (*ref == 0) if (*ref == 0)
{ {
@ -178,17 +178,17 @@ namespace CosmoTool
} }
public: public:
static ProgressiveOutput<T> static ProgressiveOutput<T>
saveArrayProgressive(const std::string& fname, uint32_t *dimList, saveArrayProgressive(const std::string& fname, uint32_t *dimList,
uint32_t rank); uint32_t rank);
ProgressiveOutput() { ProgressiveOutput() {
impl = 0; impl = 0;
ref = 0; ref = 0;
} }
ProgressiveOutput(ProgressiveOutputImpl<T> *i) { ProgressiveOutput(ProgressiveOutputImpl<T> *i) {
impl = i; impl = i;
ref = new int; ref = new int;
*ref = 1; *ref = 1;
} }
ProgressiveOutput(const ProgressiveOutput<T>& o) { ProgressiveOutput(const ProgressiveOutput<T>& o) {
@ -222,10 +222,9 @@ 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);
}; };
#endif #endif

51
src/yorick_nc4.cpp
View File

@ -7,16 +7,16 @@ This software is a computer program whose purpose is to provide a toolbox for co
data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...)
This software is governed by the CeCILL license under French law and This software is governed by the CeCILL license under French law and
abiding by the rules of distribution of free software. You can use, abiding by the rules of distribution of free software. You can use,
modify and/ or redistribute the software under the terms of the CeCILL modify and/ or redistribute the software under the terms of the CeCILL
license as circulated by CEA, CNRS and INRIA at the following URL license as circulated by CEA, CNRS and INRIA at the following URL
"http://www.cecill.info". "http://www.cecill.info".
As a counterpart to the access to the source code and rights to copy, As a counterpart to the access to the source code and rights to copy,
modify and redistribute granted by the license, users are provided only modify and redistribute granted by the license, users are provided only
with a limited warranty and the software's author, the holder of the with a limited warranty and the software's author, the holder of the
economic rights, and the successive licensors have only limited economic rights, and the successive licensors have only limited
liability. liability.
In this respect, the user's attention is drawn to the risks associated In this respect, the user's attention is drawn to the risks associated
with loading, using, modifying and/or developing or reproducing the with loading, using, modifying and/or developing or reproducing the
@ -25,9 +25,9 @@ that may mean that it is complicated to manipulate, and that also
therefore means that it is reserved for developers and experienced therefore means that it is reserved for developers and experienced
professionals having in-depth computer knowledge. Users are therefore professionals having in-depth computer knowledge. Users are therefore
encouraged to load and test the software's suitability as regards their encouraged to load and test the software's suitability as regards their
requirements in conditions enabling the security of their systems and/or requirements in conditions enabling the security of their systems and/or
data to be ensured and, more generally, to use and operate it in the data to be ensured and, more generally, to use and operate it in the
same conditions as regards security. same conditions as regards security.
The fact that you are presently reading this means that you have had The fact that you are presently reading this means that you have had
knowledge of the CeCILL license and that you accept its terms. knowledge of the CeCILL license and that you accept its terms.
@ -53,7 +53,7 @@ public:
vector<size_t> counts; vector<size_t> counts;
vector<NcDim> dimList; vector<NcDim> dimList;
uint32_t rank; uint32_t rank;
NetCDF_handle(NcFile *f, NcVar v, vector<NcDim>& dimList, uint32_t rank); NetCDF_handle(NcFile *f, NcVar v, vector<NcDim>& dimList, uint32_t rank);
virtual ~NetCDF_handle(); virtual ~NetCDF_handle();
}; };
@ -86,14 +86,14 @@ public:
InputGenCDF(NcFile *f, NcVar v, vector<NcDim>& dimList, uint32_t rank) InputGenCDF(NcFile *f, NcVar v, vector<NcDim>& dimList, uint32_t rank)
: NetCDF_handle(f,v,dimList,rank) : NetCDF_handle(f,v,dimList,rank)
{} {}
virtual ~InputGenCDF() {} virtual ~InputGenCDF() {}
virtual T read() virtual T read()
{ {
T a; T a;
curVar.getVar(curPos, counts, &a); curVar.getVar(curPos, counts, &a);
curPos[rank-1]++; curPos[rank-1]++;
for (long i = rank-1; i >= 1; i--) for (long i = rank-1; i >= 1; i--)
{ {
@ -102,7 +102,7 @@ public:
curPos[i-1]++; curPos[i-1]++;
curPos[i] = 0; curPos[i] = 0;
} }
} }
return a; return a;
} }
@ -120,12 +120,12 @@ public:
OutputGenCDF(NcFile *f, NcVar v, vector<NcDim>& dimList, uint32_t rank) OutputGenCDF(NcFile *f, NcVar v, vector<NcDim>& dimList, uint32_t rank)
: NetCDF_handle(f,v,dimList,rank) : NetCDF_handle(f,v,dimList,rank)
{} {}
virtual ~OutputGenCDF() {} virtual ~OutputGenCDF() {}
virtual void put(T a) virtual void put(T a)
{ {
curVar.putVar(curPos, counts, &a); curVar.putVar(curPos, counts, &a);
curPos[rank-1]++; curPos[rank-1]++;
for (long i = rank-1; i >= 1; i--) for (long i = rank-1; i >= 1; i--)
{ {
@ -134,7 +134,7 @@ public:
curPos[i-1]++; curPos[i-1]++;
curPos[i] = 0; curPos[i] = 0;
} }
} }
} }
}; };
@ -159,17 +159,17 @@ namespace CosmoTool {
uint32_t rank) uint32_t rank)
{ {
NcFile *f = new NcFile(fname, NcFile::replace); NcFile *f = new NcFile(fname, NcFile::replace);
vector<NcDim> dimArray; vector<NcDim> dimArray;
for (uint32_t i = 0; i < rank; i++) for (uint32_t i = 0; i < rank; i++)
{ {
char dimName[255]; char dimName[255];
sprintf(dimName, "dim%d", i); sprintf(dimName, "dim%d", i);
dimArray.push_back(f->addDim(dimName, dimList[rank-1-i])); dimArray.push_back(f->addDim(dimName, dimList[rank-1-i]));
} }
NcVar v = f->addVar("array", get_NetCDF_type<T>(), dimArray); NcVar v = f->addVar("array", get_NetCDF_type<T>(), dimArray);
vector<NcDim> ldimList; vector<NcDim> ldimList;
@ -177,7 +177,7 @@ namespace CosmoTool {
ldimList.push_back(dimArray[i]); ldimList.push_back(dimArray[i]);
OutputGenCDF<T> *impl = new OutputGenCDF<T>(f, v, ldimList, rank); OutputGenCDF<T> *impl = new OutputGenCDF<T>(f, v, ldimList, rank);
return ProgressiveOutput<T>(impl); return ProgressiveOutput<T>(impl);
} }
template<typename T> template<typename T>
@ -206,32 +206,31 @@ namespace CosmoTool {
const T *array, uint32_t *dimList, uint32_t rank) const T *array, uint32_t *dimList, uint32_t rank)
{ {
NcFile f(fname.c_str(), NcFile::replace); NcFile f(fname.c_str(), NcFile::replace);
vector<NcDim> dimArray; vector<NcDim> dimArray;
for (uint32_t i = 0; i < rank; i++) for (uint32_t i = 0; i < rank; i++)
{ {
char dimName[255]; char dimName[255];
sprintf(dimName, "dim%d", i); sprintf(dimName, "dim%d", i);
dimArray.push_back(f.addDim(dimName, dimList[i])); dimArray.push_back(f.addDim(dimName, dimList[i]));
} }
NcVar v = f.addVar("array", get_NetCDF_type<T>(), dimArray); NcVar v = f.addVar("array", get_NetCDF_type<T>(), dimArray);
v.putVar(array); v.putVar(array);
} }
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);
//if (!f.is_valid()) //if (!f.is_valid())
// throw NoSuchFileException(fname); // throw NoSuchFileException(fname);
NcVar v = f.getVar("array"); NcVar v = f.getVar("array");
vector<NcDim> dims = v.getDims(); vector<NcDim> dims = v.getDims();
rank = v.getDimCount(); rank = v.getDimCount();
uint32_t fullSize = 1; uint32_t fullSize = 1;
@ -268,5 +267,5 @@ namespace CosmoTool {
const float *array, uint32_t *dimList, uint32_t rank); const float *array, uint32_t *dimList, uint32_t rank);
template void saveArray<double>(const std::string& fname, template void saveArray<double>(const std::string& fname,
const double *array, uint32_t *dimList, uint32_t rank); const double *array, uint32_t *dimList, uint32_t rank);
} }