In [16]:
import os, os.path
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline
from scipy.special import sph_harm
import weave
from weave import converters
from numba import jit, vectorize, guvectorize
import numexpr
%load_ext Cython
In [2]:
from colloids import boo, particles
In [3]:
path = '/data/mleocmach/Documents/Tsurusawa/0_Data_retracked/3_DenseGel/163A/2_ageing/'
pos = np.loadtxt(os.path.join(path, '163A_1415_ageing_t000.dat'), skiprows=2)
bonds = np.loadtxt(os.path.join(path, '163A_1415_ageing_t000.bonds'), dtype=int)
inside = np.min((pos-pos.min(0)>14) & (pos.max()-pos>14), -1)
ngbs = particles.bonds2ngbs(bonds, len(pos))
In [13]:
%timeit boo.weave_qlm(pos, ngbs, inside)
In [4]:
%timeit boo.bonds2qlm(pos, bonds)
In [22]:
import math
def python_qlm(pos, ngbs, inside, l=6):
qlm = np.zeros([len(pos), l+1], np.complex128)
cart = np.zeros(3)
sph = np.zeros(3)
for i in range(ngbs.shape[0]):
if not inside[i]:
continue
nb = 0
for j in range(ngbs.shape[1]):
q = ngbs[i,j]
if q < 0 or q >= len(pos):
continue
nb += 1
cart[:] = pos[i] - pos[q]
sph[0] = cart[0]*cart[0] + cart[1]*cart[1] + cart[2]*cart[2]
if cart[2]**2 == sph[0] or sph[0]+1.0 == 1.0:
sph[1:] = 0
else:
sph[0] = math.sqrt(sph[0])
sph[1] = math.acos(cart[2]/sph[0])
sph[2] = math.atan2(cart[1], cart[0])
if sph[2] < 0:
sph[2] += 2.0*math.pi
qlm[i] += sph_harm(
np.arange(l+1), l,
sph[2],
sph[1]
)
if nb>0:
qlm[i] /= nb
return qlm
In [23]:
%timeit python_qlm(pos, ngbs, inside)
In [4]:
%%cython --annotate
import numpy as np
from scipy.special.cython_special cimport sph_harm
cimport numpy as np
from libc.math cimport sqrt, acos, atan2, M_PI
cimport cython
@cython.boundscheck(False) # turn off bounds-checking for entire function
@cython.wraparound(False) # turn off negative index wrapping for entire function
@cython.cdivision(True) # turn off division by zero checking for entire function
def cython_qlm(
np.ndarray[np.float64_t, ndim=2] pos,
np.ndarray[np.int64_t, ndim=2] ngbs,
np.ndarray[np.uint8_t, ndim=1, cast=True] inside,
int l=6
):
cdef int i,j, q, nb, m, d
cdef double r, r2, theta, phi, z
cdef long Npos = pos.shape[0]
cdef long Nngb = ngbs.shape[1]
cdef np.ndarray[dtype=np.complex128_t, ndim=2] qlm = np.zeros([len(pos), l+1], np.complex128)
#cdef np.ndarray[dtype=double, ndim=1] cart = np.zeros(3, float)
cdef double[3] cart
for i in range(Npos):
if not inside[i]:
continue
nb = 0
for j in range(Nngb):
q = ngbs[i,j]
if q < 0 or q >= Npos:
continue
nb += 1
for d in range(3):
cart[d] = pos[i,d] - pos[q,d]
r2 = cart[0]**2 + cart[1]**2 + cart[2]**2
if cart[2]**2 == r2 or r2+1.0 == 1.0:
theta = 0
phi = 0
else:
r = sqrt(r2)
z = cart[2]
phi = acos(z/r)
theta = atan2(cart[1], cart[0])
if theta < 0:
theta += 2.0*M_PI
for m in range(l+1):
qlm[i,m] = qlm[i,m] + sph_harm(m, l, theta, phi)
if nb>0:
for m in range(l+1):
qlm[i,m] = qlm[i,m] / nb
return qlm
Out[4]:
In [6]:
%timeit cython_qlm(pos, ngbs, inside)
In [94]:
sph_harm?
In [134]:
qlm_w = boo.weave_qlm(pos, ngbs, inside)
In [201]:
qlm_c = cython_qlm(pos, ngbs, inside)
In [202]:
np.mean(np.abs(qlm_w - qlm_c)<1e-16)
Out[202]:
In [187]:
qlm_w[0]
Out[187]:
In [188]:
qlm_c[0]
Out[188]:
In [183]:
sph_harm(0, 6, 3.685749891321086, 1.4481876153314692)
Out[183]:
In [182]:
sph_harm?
In [117]:
from colloids import periodic
def cart2sph(cartesian):
"""Convert Cartesian coordinates [[x,y,z],] to spherical coordinates [[r,phi,theta],]
phi is cologitudinal and theta azimutal"""
spherical = np.zeros_like(cartesian)
#distances
c2 = cartesian**2
r2 = c2.sum(-1)
spherical[:,0] = np.sqrt(r2)
#work only on non-zero, non purely z vectors
sel = (r2 > c2[:,0]) | (r2+1.0 > 1.0)
x, y, z = cartesian[sel].T
r = spherical[sel,0]
#colatitudinal phi [0, pi[
spherical[sel,1] = np.arccos(z/r)
#azimutal (longitudinal) theta [0, 2pi[
theta = np.arctan2(y, x)
theta[theta<0] += 2*np.pi
spherical[sel,2] = theta
return spherical
def vect2Ylm(v, l):
"""Projects vectors v on the base of spherical harmonics of degree l."""
spherical = cart2sph(v)
return sph_harm(
np.arange(l+1)[:,None], l,
spherical[:,2][None,:],
spherical[:,1][None,:]
)
def bonds2qlm(pos, bonds, l=6, periods=-1):
"""Returns the qlm for every particle"""
qlm = np.zeros((len(pos), l+1), np.complex128)
#spherical harmonic coefficients for each bond
Ylm = vect2Ylm(
periodic.periodify(
pos[bonds[:,0]],
pos[bonds[:,1]],
periods
),
l
).T
#bin bond into each particle belonging to it
np.add.at(qlm, bonds[:,0], Ylm)
np.add.at(qlm, bonds[:,1], Ylm)
#divide by the number of bonds each particle belongs to
Nngb = np.zeros(len(pos), int)
np.add.at(Nngb, bonds.ravel(), 1)
return qlm / np.maximum(1, Nngb)[:,None]
In [118]:
%timeit bonds2qlm(pos, bonds)
In [115]:
from colloids import periodic
from math import acos, atan2
@jit(['float64[:,:](float64[:,:])'], nopython=True)
def cart2sph(cartesian):
"""Convert Cartesian coordinates [[x,y,z],] to spherical coordinates [[r,phi,theta],]
phi is cologitudinal and theta azimutal"""
spherical = np.zeros_like(cartesian)
#distances
c2 = cartesian**2
r2 = c2.sum(-1)
spherical[:,0] = np.sqrt(r2)
for i in range(len(cartesian)):
#phi and theta are defined only on non-zero, non purely z vectors
if r2[i] > c2[i,2] or r2[i]+1.0 > 1.0:
#colatitudinal phi [0, pi[
spherical[i,1] = acos(cartesian[i,2]/spherical[i,0])
#azimutal (longitudinal) theta [0, 2pi[
theta = atan2(cartesian[i,1], cartesian[i,0])
if theta<0:
theta += 2*np.pi
spherical[i,2] = theta
return spherical
#@jit(['complex128[:,:](float64[:,:], int64)'], nopython=True)
def vect2Ylm(v, l):
"""Projects vectors v on the base of spherical harmonics of degree l."""
spherical = cart2sph(v)
return sph_harm(
np.arange(l+1)[:,None], l,
spherical[:,2][None,:],
spherical[:,1][None,:]
)
def bonds2qlm(pos, bonds, l=6, periods=-1):
"""Returns the qlm for every particle"""
qlm = np.zeros((len(pos), l+1), np.complex128)
#spherical harmonic coefficients for each bond
Ylm = vect2Ylm(
periodic.periodify(
pos[bonds[:,0]],
pos[bonds[:,1]],
periods
),
l
).T
#bin bond into each particle belonging to it
np.add.at(qlm, bonds[:,0], Ylm)
np.add.at(qlm, bonds[:,1], Ylm)
#divide by the number of bonds each particle belongs to
Nngb = np.zeros(len(pos), int)
np.add.at(Nngb, bonds.ravel(), 1)
return qlm / np.maximum(1, Nngb)[:,None]
In [116]:
%timeit bonds2qlm(pos, bonds)
In [5]:
qlm = boo.bonds2qlm(pos, bonds)
In [30]:
%timeit boo.wl(qlm)
In [28]:
from colloids.boo import get_qlm, get_w3j
def python_wl(qlm):
l = qlm.shape[1]-1
w = np.zeros(qlm.shape[0])
for m1 in range(-l, l+1):
for m2 in range(-l, l+1):
m3 = -m1-m2
if -l<=m3 and m3<=l:
w+= get_w3j(l, [m1, m2, m3]) * (get_qlm(qlm, m1) * get_qlm(qlm, m2) * get_qlm(qlm, m3)).real
return w.real
In [29]:
%timeit python_wl(qlm)
In [18]:
np.sum(boo.wl(qlm) - python_wl(qlm).real > 1e-12)
Out[18]:
In [23]:
@jit
def numba_wl(qlm):
l = qlm.shape[1]-1
w = np.zeros(qlm.shape[0])
for m1 in range(-l, l+1):
for m2 in range(-l, l+1):
m3 = -m1-m2
if -l<=m3 and m3<=l:
w+= boo.get_w3j(l, [m1, m2, m3]) * (boo.get_qlm(qlm, m1) * boo.get_qlm(qlm, m2) * boo.get_qlm(qlm, m3)).real
return w
In [24]:
%timeit numba_wl(qlm)
In [13]:
[len(it) for it in _w3j]
Out[13]:
In [11]:
w3j = np.zeros([6,36])
w3j[0,0] = 1
w3j[1,:4] = np.sqrt([2/35., 1/70., 2/35., 3/35.])*[-1,1,1,-1]
w3j[2,:9] = np.sqrt([
2/1001., 1/2002., 11/182., 5/1001.,
7/286., 5/143., 14/143., 35/143., 5/143.,
])*[3, -3, -1/3.0, 2, 1, -1/3.0, 1/3.0, -1/3.0, 1]
w3j[3,:16] = np.sqrt([
1/46189., 1/46189.,
11/4199., 105/46189.,
1/46189., 21/92378.,
1/46189., 35/46189., 14/46189.,
11/4199., 21/4199., 7/4199.,
11/4199., 77/8398., 70/4199., 21/4199.
])*[-20, 10, 1, -2, -43/2.0, 3, 4, 2.5, -6, 2.5, -1.5, 1, 1, -1, 1, -2]
w3j[4,:25] = np.sqrt([
10/96577., 5/193154.,
1/965770., 14/96577.,
1/965770., 66/482885.,
5/193154., 3/96577., 77/482885.,
65/14858., 5/7429., 42/37145.,
65/14858., 0.0, 3/7429., 66/37145.,
13/74290., 78/7429., 26/37145., 33/37145.,
26/37145., 13/37145., 273/37145., 429/37145., 11/7429.,
])*[
7, -7, -37, 6, 73, -3,
-5, -8, 6, -1, 3, -1,
1, 0, -3, 2, 7, -1, 3, -1,
1, -3, 1, -1, 3]
w3j[5,:36] = np.sqrt([
7/33393355., 7/33393355.,
7/33393355., 462/6678671.,
7/33393355., 1001/6678671.,
1/233753485., 77/6678671., 6006/6678671.,
1/233753485., 55/46750697., 1155/13357342.,
1/233753485., 2926/1757545., 33/46750697., 3003/6678671.,
119/1964315., 22/2750041., 1914/474145., 429/5500082.,
17/13750205., 561/2750041., 77/392863., 143/27500410., 2002/392863.,
323/723695., 1309/20677., 374/144739., 143/144739., 1001/206770.,
323/723695., 7106/723695., 561/723695., 2431/723695., 2002/103385., 1001/103385.
])*[
-126, 63, 196/3.0, -7, -259/2.0, 7/3.0,
1097/3.0, 59/6.0, -2,
4021/6.0, -113/2.0, 3,
-914, 1/3.0, 48, -3,
-7/3.0, 65/3.0, -1, 3,
214/3.0, -3, -2/3.0, 71/3.0, -1,
3, -1/3.0, 5/3.0, -2, 1/3.0,
2/3.0, -1/3.0, 2, -4/3.0, 2/3.0, -1]
In [12]:
_w3j_m1_offset = np.array([0,1,2,4,6,9,12,16,20,25,30], int)
In [45]:
@jit#(['float64(int64, int64[:])'], nopython=True)
def get_w3j(l, ms):
sm = np.sort(np.abs(ms))
return w3j[l//2, _w3j_m1_offset[sm[-1]] + sm[0]]
@jit#(nopython=True)
#@guvectorize(['(complex128[:], int64[:], complex128[:])'], '(n),()->()', nopython=True)
#@vectorize(['complex128(complex128[:], int64)'], nopython=True)
def get_qlm(qlm, m):
if m>=0:
return qlm[...,m]
elif (-m)%2 == 0:
return np.conj(qlm[...,-m])
else:
return -np.conj(qlm[...,-m])
@jit#(nopython=True)
def wl(qlm):
"""Third order rotational invariant of the bond orientational order of l-fold symmetry
$$ w_\ell = \sum_{m_1+m_2+m_3=0}
\left( \begin{array}{ccc}
\ell & \ell & \ell \\
m_1 & m_2 & m_3
\end{array} \right)
q_{\ell m_1} q_{\ell m_2} q_{\ell m_3}
$$"""
l = qlm.shape[-1]-1
w = np.zeros(qlm.shape[:-1])
for m1 in range(-l, l+1):
for m2 in range(-l, l+1):
m3 = -m1-m2
if -l<=m3 and m3<=l:
w+= get_w3j(l, np.array([m1, m2, m3])) * (get_qlm(qlm, m1) * get_qlm(qlm, m2) * get_qlm(qlm, m3)).real
return w
In [59]:
%timeit wl(qlm)
In [49]:
wl(qlm[0]),wl(qlm[:1])
Out[49]:
In [69]:
reload(boo)
Out[69]:
In [70]:
%timeit boo.wl(qlm)
In [22]:
def boo_product(qlm1, qlm2):
"""Product between two qlm"""
n = np.atleast_2d(numexpr.evaluate(
"""real(complex(real(a), -imag(a)) * b)""",
{'a':qlm1, 'b':qlm2}
))
p = numexpr.evaluate(
"""4*pi/(2*l+1)*(2*na + nb)""",
{
'na': n[:,1:].sum(-1),
'nb': n[:,0],
'l': n.shape[1]-1,
'pi': np.pi
})
return p
In [23]:
boo_product(qlm[0], qlm[1])
Out[23]:
In [26]:
%timeit boo_product(qlm[0], qlm[1])
In [27]:
%timeit boo_product(qlm[:1000], qlm[-1000:])
In [30]:
def numpy_boo_product(qlm1, qlm2):
"""Product between two qlm"""
n = np.atleast_2d((qlm1 * np.conj(qlm2)).real)
l = n.shape[1]-1
return 4*np.pi/(2*l+1)*(2*n[:,1:].sum(-1) + n[:,0])
In [31]:
numpy_boo_product(qlm[0], qlm[1])
Out[31]:
In [32]:
%timeit numpy_boo_product(qlm[0], qlm[1])
In [33]:
%timeit numpy_boo_product(qlm[:1000], qlm[-1000:])
In [119]:
@jit([
"float64[:](complex128[:,:], complex128[:,:])",
"float64[:](complex128[:], complex128[:])",
"float64[:](complex128[:], complex128[:,:])",
"float64[:](complex128[:,:], complex128[:])"
], nopython=True)
def numba_boo_product(qlm1, qlm2):
"""Product between two qlm"""
n = np.atleast_2d((qlm1 * np.conj(qlm2)).real)
l = n.shape[1]-1
return 4*np.pi/(2*l+1)*(2*n[:,1:].sum(-1) + n[:,0])
In [120]:
numba_boo_product(qlm[0], qlm[1])
Out[120]:
In [121]:
%timeit numba_boo_product(qlm[0], qlm[1])
In [122]:
%timeit numba_boo_product(qlm[:1000], qlm[-1000:])
In [123]:
%timeit numba_boo_product(qlm[0], qlm[-1000:])
In [124]:
%timeit numba_boo_product(qlm[:1000], qlm[-1])
In [125]:
from numba import guvectorize
@guvectorize(['void(complex128[:], complex128[:], float64[:])'], '(n),(n)->()', nopython=True)
def numba2_boo_product(qlm1, qlm2, prod):
"""Product between two qlm"""
l = qlm1.shape[0]-1
prod[0] = (qlm1[0] * qlm2[0].conjugate()).real
for i in range(1, len(qlm1)):
prod[0] += 2 * (qlm1[i] * qlm2[i].conjugate()).real
prod[0] *= 4*np.pi/(2*l+1)
In [126]:
numba2_boo_product(qlm[0], qlm[1])
Out[126]:
In [127]:
%timeit numba2_boo_product(qlm[0], qlm[1])
In [128]:
%timeit numba2_boo_product(qlm[:1000], qlm[-1000:])
In [129]:
%timeit numba2_boo_product(qlm[0], qlm[-1000:])
In [130]:
%timeit numba2_boo_product(qlm[:1000], qlm[-1])
In [149]:
@jit(nopython=True)
def numba3_boo_product(qlm1, qlm2):
"""Product between two qlm"""
l = qlm1.shape[-1]-1
prod = (qlm1[...,0] * np.conj(qlm2[...,0])).real
prod += 2 * (qlm1[...,1:] * np.conj(qlm2[...,1:])).real.sum(-1)
return 4 * np.pi / (2 * l + 1) * prod
In [150]:
numba3_boo_product(qlm[0], qlm[1])
Out[150]:
In [151]:
%timeit numba3_boo_product(qlm[0], qlm[1])
In [153]:
%timeit numba3_boo_product(qlm[:1000], qlm[-1000:])
In [147]:
%timeit numba3_boo_product(qlm[0], qlm[-1000:])
In [148]:
%timeit numba3_boo_product(qlm[:1000], qlm[-1])
In [81]:
@vectorize#(['float64(complex128)', 'float32(complex64)'])
def abs2(x):
return x.real**2 + x.imag**2
def ql(qlm):
"""Second order rotational invariant of the bond orientational order of l-fold symmetry
$$ q_\ell = \sqrt{\frac{4\pi}{2l+1} \sum_{m=-\ell}^{\ell} |q_{\ell m}|^2 } $$"""
q = abs2(qlm[:,0])
q += 2*abs2(qlm[:,1:]).sum(-1)
l = qlm.shape[1]-1
return np.sqrt(4*np.pi / (2*l+1) * q)
In [83]:
%timeit ql(qlm)
In [87]:
@jit(nopython=True)
def numba_ql(qlm):
"""Second order rotational invariant of the bond orientational order of l-fold symmetry
$$ q_\ell = \sqrt{\frac{4\pi}{2l+1} \sum_{m=-\ell}^{\ell} |q_{\ell m}|^2 } $$"""
q = abs2(qlm[...,0])
q += 2*abs2(qlm[...,1:]).sum(-1)
l = qlm.shape[-1]-1
return np.sqrt(4*np.pi / (2*l+1) * q)
In [89]:
%timeit numba_ql(qlm)
In [90]:
reload(boo)
Out[90]:
In [92]:
%timeit boo.ql(qlm)
In [6]:
def coarsegrain_qlm(qlm, bonds, inside):
"""Coarse grain the bond orientational order on the neighbourhood of a particle
$$Q_{\ell m}(i) = \frac{1}{N_i+1}\left( q_{\ell m}(i) + \sum_{j=0}^{N_i} q_{\ell m}(j)\right)$$
Returns Qlm and the mask of the valid particles
"""
#Valid particles must be valid themselves have only valid neighbours
inside2 = np.copy(inside)
np.bitwise_and.at(inside2, bonds[:,0], inside[bonds[:,1]])
np.bitwise_and.at(inside2, bonds[:,1], inside[bonds[:,0]])
#number of neighbours
Nngb = np.zeros(len(qlm), int)
np.add.at(Nngb, bonds.ravel(), 1)
#sum the boo coefficients of all the neighbours
Qlm = np.zeros_like(qlm)
np.add.at(Qlm, bonds[:,0], qlm[bonds[:,1]])
np.add.at(Qlm, bonds[:,1], qlm[bonds[:,0]])
Qlm[np.bitwise_not(inside2)] = 0
return Qlm / np.maximum(1, Nngb)[:,None], inside2
In [7]:
Qlm, inside2 = coarsegrain_qlm(qlm, bonds, inside)
In [8]:
%timeit coarsegrain_qlm(qlm, bonds, inside)
In [14]:
@jit
def numba_coarsegrain_qlm(qlm, bonds, inside):
"""Coarse grain the bond orientational order on the neighbourhood of a particle
$$Q_{\ell m}(i) = \frac{1}{N_i+1}\left( q_{\ell m}(i) + \sum_{j=0}^{N_i} q_{\ell m}(j)\right)$$
Returns Qlm and the mask of the valid particles
"""
#Valid particles must be valid themselves have only valid neighbours
inside2 = np.copy(inside)
np.bitwise_and.at(inside2, bonds[:,0], inside[bonds[:,1]])
np.bitwise_and.at(inside2, bonds[:,1], inside[bonds[:,0]])
#number of neighbours
Nngb = np.zeros(len(qlm), int)
np.add.at(Nngb, bonds.ravel(), 1)
#sum the boo coefficients of all the neighbours
Qlm = np.zeros_like(qlm)
np.add.at(Qlm, bonds[:,0], qlm[bonds[:,1]])
np.add.at(Qlm, bonds[:,1], qlm[bonds[:,0]])
Qlm[np.bitwise_not(inside2)] = 0
return Qlm / np.maximum(1, Nngb)[:,None], inside2
In [15]:
%timeit numba_coarsegrain_qlm(qlm, bonds, inside)
In [16]:
maxdist = 30.0
bounds = np.vstack((pos[inside2].min(0)+maxdist, pos[inside2].max(0)-maxdist))
is_center = (pos>bounds[0]).min(1) & (pos<bounds[1]).min(1)
In [59]:
hq, hQ, g = boo.gG_l(pos, qlm, Qlm, is_center, Nbins=250, maxdist=30.0)
In [60]:
plt.plot(hq, label='hq')
plt.plot(hQ, label='hQ')
plt.plot(g, label='g')
plt.legend()
Out[60]:
In [71]:
good = g>1
plt.plot(np.where(good)[0], hq[good]/g[good], label='hq/g')
plt.plot(np.where(good)[0], hQ[good]/g[good], label='hQ/g')
plt.yscale('log')
plt.legend()
Out[71]:
In [255]:
%timeit boo.gG_l(pos, qlm, Qlm, is_center, Nbins=250, maxdist=30.0)
In [66]:
maxdist = 60.0
bounds = np.vstack((pos[inside2].min(0)+maxdist, pos[inside2].max(0)-maxdist))
is_center2 = (pos>bounds[0]).min(1) & (pos<bounds[1]).min(1)
In [67]:
hq2, hQ2, g2 = boo.gG_l(pos, qlm, Qlm, is_center2, Nbins=250, maxdist=60.0)
In [68]:
plt.plot(hq2, label='hq')
plt.plot(hQ2, label='hQ')
plt.plot(g2, label='g')
plt.legend()
Out[68]:
In [70]:
good = g2>1
plt.plot(np.where(good)[0], hq2[good]/g2[good], label='hq/g')
plt.plot(np.where(good)[0], hQ2[good]/g2[good], label='hQ/g')
plt.yscale('log')
plt.legend()
Out[70]:
In [45]:
%timeit boo.gG_l(pos, qlm, Qlm, is_center2, Nbins=1000, maxdist=120.0)
In [106]:
from scipy.spatial import cKDTree as KDTree
from colloids.boo import boo_product
def gG_l(pos, qlms, Qlms, is_center, Nbins, maxdist):
"""Spatial correlation of the qlms and the Qlms (non normalized.
For each particle tagged as is_center, do the cross product between their qlm, their Qlm and count,
then bin each quantity with respect to distance.
The two first sums need to be normalised by the last one.
- pos is a Nxd array of coordinates, with d the dimension of space
- qlm is a Nx(2l+1) array of boo coordinates for l-fold symmetry
- Qlm is the coarse-grained version of qlm
- is_center is a N array of booleans. For example all particles further away than maxdist from any edge of the box.
- Nbins is the number of bins along r
- maxdist is the maximum distance considered"""
assert len(pos) == len(qlms)
assert len(qlms) == len(Qlms)
assert len(is_center) == len(pos)
maxsq = float(maxdist**2)
hQ = np.zeros(Nbins)
hq = np.zeros(Nbins)
g = np.zeros(Nbins, int)
#spatial indexing
tree = KDTree(pos, 12)
for i in np.where(is_center)[0]:
js = np.array(tree.query_ball_point(pos[i], maxdist))
js = js[js!=i]
rs = np.sqrt(np.sum((pos[js] - pos[i])**2, -1)) / maxdist * g.shape[0]
pqs = boo_product(qlms[i][None,:], qlms[js])
pQs = boo_product(Qlms[i][None,:], Qlms[js])
np.add.at(g, rs.astype(int), 1)
np.add.at(hq, rs.astype(int), pqs)
np.add.at(hQ, rs.astype(int), pQs)
return hq, hQ, g
In [102]:
hq1, hQ1, g1 = gG_l(pos, qlm, Qlm, is_center, Nbins=250, maxdist=30.0)
In [84]:
plt.plot(hq, label='hq')
plt.plot(hQ, label='hQ')
plt.plot(g, label='g')
plt.plot(hq1, ls=':', label='hq1')
plt.plot(hQ1, ls=':', label='hQ1')
plt.plot(g1, ls=':', label='g1')
plt.legend()
Out[84]:
In [86]:
np.abs(g - g1).max(), np.abs(hq - hq1).max(), np.abs(hQ - hQ1).max()
Out[86]:
In [96]:
%timeit gG_l(pos, qlm, Qlm, is_center, Nbins=250, maxdist=30.0)
In [92]:
%timeit boo.gG_l(pos, qlm, Qlm, is_center2, Nbins=1000, maxdist=120.0)
In [235]:
def gG_l_tree(pos, qlms, Qlms, is_center, Nbins, maxdist):
"""Spatial correlation of the qlms and the Qlms (non normalized.
For each particle tagged as is_center, do the cross product between their qlm, their Qlm and count,
then bin each quantity with respect to distance.
The two first sums need to be normalised by the last one.
- pos is a Nxd array of coordinates, with d the dimension of space
- qlm is a Nx(2l+1) array of boo coordinates for l-fold symmetry
- Qlm is the coarse-grained version of qlm
- is_center is a N array of booleans. For example all particles further away than maxdist from any edge of the box.
- Nbins is the number of bins along r
- maxdist is the maximum distance considered"""
assert len(pos) == len(qlms)
assert len(qlms) == len(Qlms)
assert len(is_center) == len(pos)
rsq2r = np.sqrt(np.arange(Nbins**2)).astype(int)
l2r = float((Nbins/maxdist)**2)
#maxsq = float(maxdist**2)
hQ = np.zeros(Nbins)
hq = np.zeros(Nbins)
g = np.zeros(Nbins, int)
#spatial indexing
tree = KDTree(pos, 12)
centertree = KDTree(pos[is_center], 12)
for i, js in zip(np.where(is_center)[0], centertree.query_ball_tree(tree, maxdist)):
#for i in np.where(is_center)[0]:
#js = np.array(tree.query_ball_point(pos[i], maxdist))
js = np.array(js)
js = js[js!=i]
#rs = np.sqrt(np.sum((pos[js] - pos[i])**2, -1)) / maxdist * g.shape[0]
#rsqs = numexpr.evaluate('sum((a-b)**2, axis=1)', {'a':pos[js], 'b':pos[i]}) * l2r
rsqs = np.sum((pos[js] - pos[i])**2, -1) * l2r
rs = rsq2r[rsqs.astype(int)]
pqs = boo_product(qlms[i][None,:], qlms[js])
#pQs = boo_product(Qlms[i][None,:], Qlms[js])
np.add.at(g, rs, 1)
np.add.at(hq, rs.astype(int), pqs)
#np.add.at(hQ, rs.astype(int), pQs)
return hq, hQ, g
In [236]:
%timeit gG_l_tree(pos, qlm, Qlm, is_center, Nbins=250, maxdist=30.0)
In [266]:
def gG_l_record(pos, qlms, Qlms, is_center, Nbins, maxdist):
"""Spatial correlation of the qlms and the Qlms (non normalized.
For each particle tagged as is_center, do the cross product between their qlm, their Qlm and count,
then bin each quantity with respect to distance.
The two first sums need to be normalised by the last one.
- pos is a Nxd array of coordinates, with d the dimension of space
- qlm is a Nx(2l+1) array of boo coordinates for l-fold symmetry
- Qlm is the coarse-grained version of qlm
- is_center is a N array of booleans. For example all particles further away than maxdist from any edge of the box.
- Nbins is the number of bins along r
- maxdist is the maximum distance considered"""
assert len(pos) == len(qlms)
assert len(qlms) == len(Qlms)
assert len(is_center) == len(pos)
#conversion factor between indices and bins
l2r = Nbins/maxdist
#result containers
hqQ = np.zeros((Nbins,2))
g = np.zeros(Nbins, int)
#spatial indexing
tree = KDTree(pos, 12)
centertree = KDTree(pos[is_center], 12)
#all pairs of points closer than maxdist with their distances in a record array
query = centertree.sparse_distance_matrix(tree, maxdist, output_type='ndarray')
#keep only pairs where the points are distinct
centerindex = np.where(is_center)[0]
query['i'] = centerindex[query['i']]
good = query['i'] != query['j']
query = query[good]
#binning of distances
rs = (query['v'] * l2r).astype(int)
np.add.at(g, rs, 1)
#binning of boo cross products
pqQs = np.empty((len(rs),2))
pqQs[:,0] = boo_product(qlms[query['i']], qlms[query['j']])
pqQs[:,1] = boo_product(Qlms[query['i']], Qlms[query['j']])
np.add.at(hqQ, rs, pqQs)
return hqQ[:,0], hqQ[:,1], g
In [267]:
hq3, hQ3, g3 = gG_l_record(pos, qlm, Qlm, is_center, Nbins=250, maxdist=30.0)
plt.subplot(1,3,1)
plt.plot(g)
plt.plot(g3, ls=':')
plt.subplot(1,3,2)
plt.plot(hq, label='hq')
plt.plot(hq3, ls=':', label='hq3')
plt.subplot(1,3,3)
plt.plot(hQ, label='hQ')
plt.plot(hQ3, ls=':', label='hQ3')
Out[267]:
In [268]:
%timeit gG_l_record(pos, qlm, Qlm, is_center, Nbins=250, maxdist=30.0)
In [212]:
%debug
In [51]:
%%cython --annotate
import numpy as np
cimport numpy as np
cimport cython
from libc.math cimport sqrt
from scipy.spatial import cKDTree as KDTree
#from scipy.spatial cimport cKDTree as KDTree
from colloids import particles, boo
from colloids.boo import boo_product
#_dtype = [('i',np.intp),('j',np.intp),('v',np.float64)]
#res_dtype = np.dtype(_dtype, align = True)
cdef packed struct ijdist:
np.intp_t i,j
np.float64_t v
@cython.boundscheck(False) # turn off bounds-checking for entire function
@cython.wraparound(False) # turn off negative index wrapping for entire function
@cython.cdivision(True) # turn off division by zero checking for entire function
def cython_gG_l(
np.ndarray[np.float64_t, ndim=2] pos,
qlms, Qlms,
is_center,
int Nbins, float maxdist):
"""Spatial correlation of the qlms and the Qlms (non normalized.
For each particle tagged as is_center, do the cross product between their qlm, their Qlm and count,
then bin each quantity with respect to distance.
The two first sums need to be normalised by the last one.
- pos is a Nxd array of coordinates, with d the dimension of space
- qlm is a Nx(2l+1) array of boo coordinates for l-fold symmetry
- Qlm is the coarse-grained version of qlm
- is_center is a N array of booleans. For example all particles further away than maxdist from any edge of the box.
- Nbins is the number of bins along r
- maxdist is the maximum distance considered"""
cdef int i,it,j,r
cdef float rsq
assert len(pos) == len(qlms)
assert len(qlms) == len(Qlms)
assert len(is_center) == len(pos)
#cdef np.ndarray[np.int64_t, ndim=1] rsq2r = np.sqrt(np.arange(Nbins**2)).astype(int)
cdef float l2r = (Nbins/maxdist)#**2
#maxsq = float(maxdist**2)
hQ = np.zeros(Nbins)
hq = np.zeros(Nbins)
g = np.zeros(Nbins, int)
#spatial indexing
tree = KDTree(pos, 12)
centertree = KDTree(pos[is_center], 12)
centerindex = np.where(is_center)[0]
cdef np.ndarray[ijdist, ndim=1] query = centertree.sparse_distance_matrix(tree, maxdist, output_type='ndarray')
for it in range(query.shape[0]):
i = centerindex[query[it].i]
if i==query[it].j:
continue
r = int(l2r * query[it].v)
g[r] += 1
return hq, hQ, g
Out[51]:
In [52]:
%timeit cython_gG_l(pos, qlm, Qlm, is_center, Nbins=250, maxdist=30.0)
Test the final version
In [104]:
reload(boo)
Out[104]:
In [105]:
hqQ, g = boo.gG_l(pos, [qlm, Qlm], is_center, Nbins=250, maxdist=30.0)
In [106]:
plt.subplot(1,3,1)
plt.plot(g)
plt.subplot(1,3,2)
plt.plot(hqQ[:,0], label='hq')
plt.subplot(1,3,3)
plt.plot(hqQ[:,1], label='hQ')
Out[106]:
In [107]:
%timeit boo.gG_l(pos, [qlm, Qlm], is_center, Nbins=250, maxdist=30.0)
In [7]:
In [ ]:
hq, g = boo.steinhardt_g_l(pos, bonds, is_center, 250, maxdist)
In [ ]: