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..a2a85d6 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,9 @@ 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 + 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' @@ -208,9 +208,9 @@ def get_mean_density(fdir, smin, step): def get_mean_density_fdir(fdir,init,steps): """ estimate ensemble mean """ - print '-'*60 - print 'Get 3D ensemble mean density field' - print '-'*60 + print('-'*60) + print('Get 3D ensemble mean density field') + print('-'*60) fname0,fname1=build_filelist(fdir) diff --git a/python/cosmotool/cic.py b/python/cosmotool/cic.py index ec463d8..5a3372a 100644 --- a/python/cosmotool/cic.py +++ b/python/cosmotool/cic.py @@ -4,7 +4,7 @@ import numpy as np def cicParticles(particles, L, N): - if type(N) not in [int,long]: + if type(N) not in [int,int]: raise TypeError("N must be a numeric type") def shifted(i, t): @@ -14,7 +14,7 @@ def cicParticles(particles, L, N): i =[] r = [] - for d in xrange(3): + for d in range(3): q = ne.evaluate('(p%L)*N/L', local_dict={'p':particles[d], 'L':L, 'N':N }) o = np.empty(q.size, dtype=np.int64) o[:] = np.floor(q) diff --git a/python/cosmotool/cl_cic.py b/python/cosmotool/cl_cic.py index b53c96b..bf0ab4e 100644 --- a/python/cosmotool/cl_cic.py +++ b/python/cosmotool/cl_cic.py @@ -163,8 +163,8 @@ class CIC_CL(object): # 2 dimensions translator['ndim'] = ndim translator['centered'] = '1' if centered else '0' - looperVariables = ','.join(['id%d' % d for d in xrange(ndim)]) - looperFor = '\n'.join(['for (int id{dim}=0; id{dim} < 2; id{dim}++) {{'.format(dim=d) for d in xrange(ndim)]) + looperVariables = ','.join(['id%d' % d for d in range(ndim)]) + looperFor = '\n'.join(['for (int id{dim}=0; id{dim} < 2; id{dim}++) {{'.format(dim=d) for d in range(ndim)]) looperForEnd = '}' * ndim kern = pragmas + CIC_PREKERNEL.format(**translator) + (CIC_KERNEL % {'looperVariables': looperVariables, 'looperFor': looperFor, 'looperForEnd':looperForEnd}) diff --git a/python/cosmotool/grafic.py b/python/cosmotool/grafic.py index 4ea6f9b..9310172 100644 --- a/python/cosmotool/grafic.py +++ b/python/cosmotool/grafic.py @@ -27,7 +27,7 @@ def readGrafic(filename): BoxSize = delta * Nx * h checkPoint = 4*Ny*Nz - for i in xrange(Nx): + for i in range(Nx): checkPoint = struct.unpack("I", f.read(4))[0] if checkPoint != 4*Ny*Nz: raise ValueError("Invalid unformatted access") @@ -57,7 +57,7 @@ def writeGrafic(filename, field, BoxSize, scalefac, **cosmo): cosmo['omega_M_0'], cosmo['omega_lambda_0'], 100*cosmo['h'], checkPoint)) checkPoint = 4*Ny*Nz field = field.reshape(field.shape, order='F') - for i in xrange(Nx): + for i in range(Nx): f.write(struct.pack("I", checkPoint)) f.write(field[i].astype(np.float32).tostring()) f.write(struct.pack("I", checkPoint)) @@ -73,7 +73,7 @@ def writeWhitePhase(filename, field): field = field.reshape(field.shape, order='F') checkPoint = struct.pack("I", 4*Nx*Ny) - for i in xrange(Nx): + for i in range(Nx): f.write(checkPoint) f.write(field[i].astype(np.float32).tostring()) f.write(checkPoint) @@ -87,7 +87,7 @@ def readWhitePhase(filename): checkPoint_ref = 4*Ny*Nz - for i in xrange(Nx): + for i in range(Nx): if struct.unpack("I", f.read(4))[0] != checkPoint_ref: raise ValueError("Invalid unformatted access") diff --git a/python/cosmotool/timing.py b/python/cosmotool/timing.py index c6977f1..cad8463 100644 --- a/python/cosmotool/timing.py +++ b/python/cosmotool/timing.py @@ -14,7 +14,7 @@ def time_block(name): yield te = time.time() - print ('%s %2.2f sec' % (name, te-ts)) + print('%s %2.2f sec' % (name, te-ts)) def timeit(method): """This decorator add a timing request for each call to the decorated function. @@ -28,7 +28,7 @@ def timeit(method): result = method(*args, **kw) te = time.time() - print ('%r (%r, %r) %2.2f sec' % (method.__name__, args, kw, te-ts)) + print('%r (%r, %r) %2.2f sec' % (method.__name__, args, kw, te-ts)) return result return timed @@ -46,7 +46,7 @@ def timeit_quiet(method): result = method(*args, **kw) te = time.time() - print ('%r %2.2f sec' % (method.__name__, te-ts)) + print('%r %2.2f sec' % (method.__name__, te-ts)) return result return timed diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 9737765..e7496e3 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -6,7 +6,6 @@ SET(CosmoTool_SRCS load_data.cpp loadGadget.cpp loadRamses.cpp - octTree.cpp powerSpectrum.cpp miniargs.cpp growthFactor.cpp 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; } + +};