Merge branch 'master' of bitbucket.org:glavaux/cosmotool
This commit is contained in:
commit
c2059ea118
16 changed files with 162 additions and 138 deletions
|
|
@ -1,5 +1,6 @@
|
|||
set(CMAKE_SHARED_MODULE_PREFIX)
|
||||
|
||||
|
||||
include_directories(${NUMPY_INCLUDE_DIRS} ${PYTHON_INCLUDE_PATH} ${CMAKE_SOURCE_DIR}/src ${CMAKE_BINARY_DIR}/src ${CMAKE_SOURCE_DIR}/python)
|
||||
|
||||
IF(CYTHON)
|
||||
|
|
@ -39,7 +40,10 @@ add_library(_project MODULE ${CMAKE_CURRENT_BINARY_DIR}/_project.cpp)
|
|||
|
||||
|
||||
|
||||
SET(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -Bsymbolic-functions")
|
||||
SET(CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -Bsymbolic-functions")
|
||||
if(APPLE)
|
||||
set(CMAKE_MODULE_LINKER_FLAGS "-undefined dynamic_lookup")
|
||||
endif()
|
||||
|
||||
target_link_libraries(_cosmotool ${CosmoTool_local} ${PYTHON_LIBRARIES} ${GSL_LIBRARIES})
|
||||
target_link_libraries(_cosmo_power ${CosmoTool_local} ${PYTHON_LIBRARIES} ${GSL_LIBRARIES})
|
||||
|
|
@ -74,6 +78,9 @@ if (NOT PYTHON_SITE_PACKAGES)
|
|||
mark_as_advanced(USER_PYTHON_SITE_PACKAGES SYSTEM_PYTHON_SITE_PACKAGES)
|
||||
endif (NOT PYTHON_SITE_PACKAGES)
|
||||
|
||||
message(STATUS "System python site: ${SYSTEM_PYTHON_SITE_PACKAGES}")
|
||||
message(STATUS "User python site: ${USER_PYTHON_SITE_PACKAGES}")
|
||||
|
||||
OPTION(INSTALL_PYTHON_LOCAL OFF)
|
||||
|
||||
IF (NOT INSTALL_PYTHON_LOCAL)
|
||||
|
|
@ -81,6 +88,7 @@ IF (NOT INSTALL_PYTHON_LOCAL)
|
|||
ELSE (NOT INSTALL_PYTHON_LOCAL)
|
||||
SET(PYTHON_SITE_PACKAGES ${USER_PYTHON_SITE_PACKAGES})
|
||||
ENDIF(NOT INSTALL_PYTHON_LOCAL)
|
||||
cmessage(STATUS "Python install location: ${PYTHON_SITE_PACKAGES}")
|
||||
|
||||
|
||||
if (WIN32 AND NOT CYGWIN)
|
||||
|
|
|
|||
|
|
@ -1,4 +1,6 @@
|
|||
#ifdef _OPENMP
|
||||
#include <omp.h>
|
||||
#endif
|
||||
#include <iostream>
|
||||
#include <boost/format.hpp>
|
||||
#include <boost/multi_array.hpp>
|
||||
|
|
@ -171,7 +173,11 @@ void CosmoTool_compute_bispectrum(
|
|||
{
|
||||
// First remap to multi_array for easy access
|
||||
size_t kNz = Nz/2+1;
|
||||
#ifdef _OPENMP
|
||||
int Ntasks = omp_get_max_threads();
|
||||
#else
|
||||
int Ntasks = 1;
|
||||
#endif
|
||||
boost::multi_array_ref<std::complex<double>, 3> a_delta(reinterpret_cast<std::complex<double>*>(delta_hat), boost::extents[Nx][Ny][kNz]);
|
||||
boost::multi_array_ref<size_t, 3> a_Nt(Ntriangles, boost::extents[Nk][Nk][Nk]);
|
||||
boost::multi_array_ref<std::complex<double>, 3> a_B(reinterpret_cast<std::complex<double>*>(B), boost::extents[Nk][Nk][Nk]);
|
||||
|
|
@ -190,6 +196,7 @@ void CosmoTool_compute_bispectrum(
|
|||
delta_mirror[n1.i1][n1.i2][n1.i3] = std::conj(delta_mirror[n2.i1][n2.i2][n2.i3]);
|
||||
}
|
||||
|
||||
#ifdef _OPENMP
|
||||
// First loop over m1
|
||||
#pragma omp parallel
|
||||
{
|
||||
|
|
@ -236,6 +243,9 @@ void CosmoTool_compute_bispectrum(
|
|||
a_B_p[q] += b_B_p[q];
|
||||
}
|
||||
}
|
||||
#else
|
||||
#warning Serial version not implemented
|
||||
#endif
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -9,6 +9,9 @@ cimport cython
|
|||
|
||||
np.import_array()
|
||||
|
||||
cdef extern from "sys/types.h":
|
||||
ctypedef np.int64_t int64_t
|
||||
|
||||
cdef extern from "loadSimu.hpp" namespace "CosmoTool":
|
||||
|
||||
cdef cppclass SimuData:
|
||||
|
|
@ -20,7 +23,7 @@ cdef extern from "loadSimu.hpp" namespace "CosmoTool":
|
|||
np.float_t Omega_Lambda
|
||||
np.int64_t TotalNumPart
|
||||
np.int64_t NumPart
|
||||
np.int64_t *Id
|
||||
int64_t *Id
|
||||
float *Pos[3]
|
||||
float *Vel[3]
|
||||
int *type
|
||||
|
|
@ -323,7 +326,7 @@ cdef object wrap_array(void *p, np.uint64_t s, int typ):
|
|||
cdef object wrap_float_array(float *p, np.uint64_t s):
|
||||
return wrap_array(<void *>p, s, np.NPY_FLOAT32)
|
||||
|
||||
cdef object wrap_int64_array(np.int64_t* p, np.uint64_t s):
|
||||
cdef object wrap_int64_array(int64_t* p, np.uint64_t s):
|
||||
return wrap_array(<void *>p, s, np.NPY_INT64)
|
||||
|
||||
cdef object wrap_int_array(int* p, np.uint64_t s):
|
||||
|
|
@ -388,6 +391,7 @@ def loadGadget(str filename, int snapshot_id, int gadgetFormat = 1, bool loadPos
|
|||
cdef int flags
|
||||
cdef SimuData *data
|
||||
cdef Simulation simu
|
||||
cdef const char *filename_bs
|
||||
|
||||
flags = 0
|
||||
if loadPosition:
|
||||
|
|
@ -401,7 +405,10 @@ def loadGadget(str filename, int snapshot_id, int gadgetFormat = 1, bool loadPos
|
|||
if loadMass:
|
||||
flags |= NEED_MASS
|
||||
|
||||
data = loadGadgetMulti(filename, snapshot_id, flags, gadgetFormat)
|
||||
filename_b = bytes(filename, 'utf-8')
|
||||
filename_bs = filename_b
|
||||
with nogil:
|
||||
data = loadGadgetMulti(filename_bs, snapshot_id, flags, gadgetFormat)
|
||||
if data == <SimuData*>0:
|
||||
return None
|
||||
|
||||
|
|
@ -497,7 +504,7 @@ def writeGadget(str filename, object simulation):
|
|||
simdata.Vel[j] = <float *>vel.data
|
||||
|
||||
ids = simulation.getIdentifiers()
|
||||
simdata.Id = <np.int64_t *>ids.data
|
||||
simdata.Id = <int64_t *>ids.data
|
||||
simdata.BoxSize = simulation.getBoxsize()
|
||||
simdata.time = simulation.getTime()
|
||||
simdata.Hubble = simulation.getHubble()
|
||||
|
|
|
|||
|
|
@ -389,7 +389,8 @@ cdef void INTERNAL_project_cic_no_mass(npx.ndarray[DTYPE_t, ndim=3] g,
|
|||
npx.ndarray[DTYPE_t, ndim=2] x, int Ngrid, double Lbox, double shifter):
|
||||
cdef double delta_Box = Ngrid/Lbox
|
||||
cdef int i
|
||||
cdef double a[3], c[3]
|
||||
cdef double a[3]
|
||||
cdef double c[3]
|
||||
cdef int b[3]
|
||||
cdef int do_not_put
|
||||
|
||||
|
|
@ -404,9 +405,6 @@ cdef void INTERNAL_project_cic_no_mass(npx.ndarray[DTYPE_t, ndim=3] g,
|
|||
if b[j] < 0 or b[j]+1 >= Ngrid:
|
||||
do_not_put = True
|
||||
|
||||
print("b = %g %g %g" % (b[0], b[1], b[2]))
|
||||
print("a = %g %g %g" % (a[0], a[1], a[2]))
|
||||
|
||||
if not do_not_put:
|
||||
g[b[0],b[1],b[2]] += c[0]*c[1]*c[2]
|
||||
g[b[0]+1,b[1],b[2]] += a[0]*c[1]*c[2]
|
||||
|
|
@ -424,8 +422,10 @@ cdef void INTERNAL_project_cic_no_mass_periodic(npx.ndarray[DTYPE_t, ndim=3] g,
|
|||
npx.ndarray[DTYPE_t, ndim=2] x, int Ngrid, double Lbox, double shifter):
|
||||
cdef double delta_Box = Ngrid/Lbox
|
||||
cdef int i
|
||||
cdef double a[3], c[3]
|
||||
cdef int b[3], b1[3]
|
||||
cdef double a[3]
|
||||
cdef double c[3]
|
||||
cdef int b[3]
|
||||
cdef int b1[3]
|
||||
cdef int do_not_put
|
||||
cdef DTYPE_t[:,:] ax
|
||||
cdef DTYPE_t[:,:,:] ag
|
||||
|
|
@ -465,7 +465,8 @@ cdef void INTERNAL_project_cic_with_mass(npx.ndarray[DTYPE_t, ndim=3] g,
|
|||
int Ngrid, double Lbox, double shifter):
|
||||
cdef double delta_Box = Ngrid/Lbox
|
||||
cdef int i
|
||||
cdef double a[3], c[3]
|
||||
cdef double a[3]
|
||||
cdef double c[3]
|
||||
cdef DTYPE_t m0
|
||||
cdef int b[3]
|
||||
|
||||
|
|
@ -502,8 +503,10 @@ cdef void INTERNAL_project_cic_with_mass_periodic(npx.ndarray[DTYPE_t, ndim=3] g
|
|||
cdef double half_Box = 0.5*Lbox, m0
|
||||
cdef double delta_Box = Ngrid/Lbox
|
||||
cdef int i
|
||||
cdef double a[3], c[3]
|
||||
cdef int b[3], b1[3]
|
||||
cdef double a[3]
|
||||
cdef double c[3]
|
||||
cdef int b[3]
|
||||
cdef int b1[3]
|
||||
|
||||
for i in range(x.shape[0]):
|
||||
|
||||
|
|
@ -550,18 +553,18 @@ def project_cic(npx.ndarray[DTYPE_t, ndim=2] x not None, npx.ndarray[DTYPE_t, nd
|
|||
if x.shape[1] != 3:
|
||||
raise ValueError("Invalid shape for x array")
|
||||
|
||||
if mass != 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")
|
||||
|
||||
g = np.zeros((Ngrid,Ngrid,Ngrid),dtype=DTYPE)
|
||||
|
||||
if not periodic:
|
||||
if mass == None:
|
||||
if mass is None:
|
||||
INTERNAL_project_cic_no_mass(g, x, Ngrid, Lbox, shifter)
|
||||
else:
|
||||
INTERNAL_project_cic_with_mass(g, x, mass, Ngrid, Lbox, shifter)
|
||||
else:
|
||||
if mass == None:
|
||||
if mass is None:
|
||||
INTERNAL_project_cic_no_mass_periodic(g, x, Ngrid, Lbox, shifter)
|
||||
else:
|
||||
INTERNAL_project_cic_with_mass_periodic(g, x, mass, Ngrid, Lbox, shifter)
|
||||
|
|
@ -629,7 +632,8 @@ cdef DTYPE_t cube_integral(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t alpha_
|
|||
@cython.cdivision(True)
|
||||
cdef DTYPE_t cube_integral_trilin(DTYPE_t u[3], DTYPE_t u0[3], int r[1], DTYPE_t vertex_value[8], DTYPE_t alpha_max) nogil:
|
||||
cdef DTYPE_t I, tmp_a
|
||||
cdef DTYPE_t v[3], term[4]
|
||||
cdef DTYPE_t v[3]
|
||||
cdef DTYPE_t term[4]
|
||||
cdef int i, j, q
|
||||
|
||||
j = 0
|
||||
|
|
@ -671,7 +675,8 @@ cdef DTYPE_t integrator1(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:
|
||||
cdef DTYPE_t vertex_value[8]
|
||||
cdef DTYPE_t d
|
||||
cdef int a[3][2], i
|
||||
cdef int a[3][2]
|
||||
cdef int i
|
||||
|
||||
for i in xrange(3):
|
||||
a[i][0] = iu0[i]
|
||||
|
|
@ -697,7 +702,10 @@ cdef DTYPE_t C_line_of_sight_projection(DTYPE_t[:,:,:] density,
|
|||
DTYPE_t min_distance,
|
||||
DTYPE_t max_distance, DTYPE_t[:] shifter, int integrator_id) nogil except? 0:
|
||||
|
||||
cdef DTYPE_t u[3], ifu0[3], u0[3], utot[3]
|
||||
cdef DTYPE_t u[3]
|
||||
cdef DTYPE_t ifu0[3]
|
||||
cdef DTYPE_t u0[3]
|
||||
cdef DTYPE_t utot[3]
|
||||
cdef int u_delta[3]
|
||||
cdef int iu0[3]
|
||||
cdef int i
|
||||
|
|
|
|||
|
|
@ -4,14 +4,23 @@ import numpy as np
|
|||
import numexpr as ne
|
||||
|
||||
class CubeFT(object):
|
||||
def __init__(self, L, N, max_cpu=-1):
|
||||
def __init__(self, L, N, max_cpu=-1, width=32):
|
||||
|
||||
if width==32:
|
||||
fourier_type='complex64'
|
||||
real_type='float32'
|
||||
elif width==64:
|
||||
fourier_type='complex128'
|
||||
real_type='float64'
|
||||
else:
|
||||
raise ValueError("Invalid bitwidth (must be 32 or 64)")
|
||||
|
||||
self.N = N
|
||||
self.align = pyfftw.simd_alignment
|
||||
self.L = float(L)
|
||||
self.max_cpu = multiprocessing.cpu_count() if max_cpu < 0 else max_cpu
|
||||
self._dhat = pyfftw.n_byte_align_empty((self.N,self.N,self.N//2+1), self.align, dtype='complex64')
|
||||
self._density = pyfftw.n_byte_align_empty((self.N,self.N,self.N), self.align, dtype='float32')
|
||||
self._dhat = pyfftw.n_byte_align_empty((self.N,self.N,self.N//2+1), self.align, dtype=fourier_type)
|
||||
self._density = pyfftw.n_byte_align_empty((self.N,self.N,self.N), self.align, dtype=real_type)
|
||||
self._irfft = pyfftw.FFTW(self._dhat, self._density, axes=(0,1,2), direction='FFTW_BACKWARD', threads=self.max_cpu)#, normalize_idft=False)
|
||||
self._rfft = pyfftw.FFTW(self._density, self._dhat, axes=(0,1,2), threads=self.max_cpu) #, normalize_idft=False)
|
||||
|
||||
|
|
|
|||
|
|
@ -14,6 +14,7 @@ class SimulationBare(PySimulationBase):
|
|||
self.positions = [q.copy() for q in s.getPositions()] if s.getPositions() is not None else None
|
||||
self.velocities = [q.copy() for q in s.getVelocities()] if s.getVelocities() is not None else None
|
||||
self.identifiers = s.getIdentifiers().copy() if s.getIdentifiers() is not None else None
|
||||
self.types = s.getTypes().copy() if s.getTypes() is not None else None
|
||||
self.boxsize = s.getBoxsize()
|
||||
self.time = s.getTime()
|
||||
self.Hubble = s.getHubble()
|
||||
|
|
@ -53,11 +54,15 @@ class SimulationBare(PySimulationBase):
|
|||
self.positions = _safe_merge(self.positions, other.getPositions())
|
||||
self.velocities = _safe_merge(self.velocities, other.getVelocities())
|
||||
self.identifiers = _safe_merge0(self.identifiers, other.getIdentifiers())
|
||||
self.types = _safe_merge0(self.types, other.getTypes())
|
||||
try:
|
||||
self.masses = _safe_merge0(self.masses, other.getMasses())
|
||||
except Exception as e:
|
||||
warnings.warn("Unexpected exception: " + repr(e));
|
||||
self.masses = None
|
||||
|
||||
def getTypes(self):
|
||||
return self.types
|
||||
|
||||
def getPositions(self):
|
||||
return self.positions
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue