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 2015-01-31 15:11:51 +01:00
commit 84e841814a
9 changed files with 93 additions and 58 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
python/_project.pyx
View File

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

14
python/cosmotool/borg.py
View File

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

4
python/cosmotool/cic.py
View File

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

4
python/cosmotool/cl_cic.py
View File

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

8
python/cosmotool/grafic.py
View File

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

6
python/cosmotool/timing.py
View File

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

1
src/CMakeLists.txt
View File

@ -6,7 +6,6 @@ SET(CosmoTool_SRCS
load_data.cpp
loadGadget.cpp
loadRamses.cpp
octTree.cpp
powerSpectrum.cpp
miniargs.cpp
growthFactor.cpp

69
src/octTree.hpp
View File

@ -55,12 +55,20 @@ namespace CosmoTool
typedef octCoordType OctCoords[3];
template<class T = void>
struct OctCell
{
octPtr numberLeaves;
octPtr children[8];
T data;
};
template<typename T>
struct OctTree_defaultUpdater {
void operator()(T& d) { }
};
template<typename T_dataUpdater = OctTree_defaultUpdater<void>, class T = void>
class OctTree
{
public:
@ -103,9 +111,10 @@ namespace CosmoTool
protected:
T_dataUpdater updater;
const FCoordinates *particles;
octPtr numParticles;
OctCell *cells;
OctCell<T> *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 };
{
octPtr newNode = cells[node].children[i];
int ipos[3] = { (i&1), (i&2)>>1, (i&4)>>2 };
if (newNode == emptyOctCell || newNode == invalidOctCell)
continue;
if (newNode == emptyOctCell || newNode == invalidOctCell)
continue;
for (int j = 0; j < 3; j++)
newCoord[j] = icoord[j]+(2*ipos[j]-1)*halfNodeLength/2;
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];
}
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

19
src/octTree.cpp → src/octTree.tcc
View File

@ -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<typename Updater, typename T>
OctTree<Updater,T>::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<T>[numCells];
Lbox = (float)(octCoordTypeNorm+1);
cells[0].numberLeaves = 0;
@ -110,12 +112,14 @@ OctTree::OctTree(const FCoordinates *particles, octPtr numParticles,
//#endif
}
OctTree::~OctTree()
template<typename Updater, typename T>
OctTree<Updater,T>::~OctTree()
{
delete cells;
}
void OctTree::buildTree(uint32_t maxAbsoluteDepth)
template<typename Updater, typename T>
void OctTree<Updater,T>::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<typename Updater, typename T>
void OctTree<Updater,T>::insertParticle(octPtr node,
const OctCoords& icoord,
octCoordType halfNodeLength,
octPtr particleId,
@ -208,3 +213,5 @@ void OctTree::insertParticle(octPtr node,
cells[node].children[octPos] = newNode;
}
};