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

Use numexpr for accumulating

Browse Source
This commit is contained in:
Guilhem Lavaux 2014-09-26 18:20:18 +02:00
parent 4e79e2eec6
commit adfc76a96d
1 changed files with 16 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
python/cosmotool/cic.py
View File

@ -11,25 +11,32 @@ def cicParticles(particles, L, N):
a = np.empty(i[0].size, dtype=np.int64)
return ne.evaluate('(i2+t2)%N + N*((i1+t1)%N + N*((i0+t0)%N) )', local_dict={'i2':i[2], 't2':t[2], 'i1':i[1], 't1':t[1], 'i0':i[0], 't0':t[0], 'N':N}, out=a)
i =[]
r = []
for d in xrange(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)
ne.evaluate('floor(q)', out=o, casting='unsafe')
o[:] = np.floor(q)
i.append(o)
r.append(ne.evaluate('q-o'))
D = {'a':r[0],'b':r[1],'c':r[2]}
N3 = N*N*N
density = np.bincount(shifted(i, (1,1,1)), weights=ne.evaluate('a*b*c', local_dict=D), minlength=N3)
density += np.bincount(shifted(i, (1,1,0)), weights=ne.evaluate('a*b*(1-c)', local_dict=D), minlength=N3)
density += np.bincount(shifted(i, (1,0,1)), weights=ne.evaluate('a*(1-b)*c', local_dict=D), minlength=N3)
density += np.bincount(shifted(i, (1,0,0)), weights=ne.evaluate('a*(1-b)*(1-c)', local_dict=D), minlength=N3)
density += np.bincount(shifted(i, (0,1,1)), weights=ne.evaluate('(1-a)*b*c', local_dict=D), minlength=N3)
density += np.bincount(shifted(i, (0,1,0)), weights=ne.evaluate('(1-a)*b*(1-c)', local_dict=D), minlength=N3)
density += np.bincount(shifted(i, (0,0,1)), weights=ne.evaluate('(1-a)*(1-b)*c', local_dict=D), minlength=N3)
density += np.bincount(shifted(i, (0,0,0)), weights=ne.evaluate('(1-a)*(1-b)*(1-c)', local_dict=D), minlength=N3)
def accum(density, ss, op):
d0 = np.bincount(shifted(i, ss), weights=ne.evaluate(op, local_dict=D), minlength=N3)
ne.evaluate('d + d0', local_dict={'d':density, 'd0':d0}, out=density)
density = np.empty(N3, dtype=np.float64)
accum(density, (1,1,1), 'a * b * c ')
accum(density, (1,1,0), 'a * b * (1-c)')
accum(density, (1,0,1), 'a * (1-b) * c ')
accum(density, (1,0,0), 'a * (1-b) * (1-c)')
accum(density, (0,1,1), '(1-a) * b * c ')
accum(density, (0,1,0), '(1-a) * b * (1-c)')
accum(density, (0,0,1), '(1-a) * (1-b) * c ')
accum(density, (0,0,0), '(1-a) * (1-b) * (1-c)')
return density.reshape((N,N,N))