In [1]:

    

# Author: Matthew Rowley <Matt.B.Rowley@gmail.com>
# Copyright (c) 2015, Matthew Rowley
# License: MIT

from mayavi import mlab
import numpy as np
from scipy.special import sph_harm
import time
from numpy import sin
from numpy import cos


    

In [2]:

    

# 2pz functions
r_2pz = lambda r: 1/np.sqrt(24)*r*np.exp(-r/2)
tp_2pz = lambda theta, phi: sph_harm(0, 1, theta, phi).real
psi_2pz = lambda r, phi, theta: r_2pz(r) * tp_2pz(theta, phi)
# 2s functions
r_2s = lambda r: (2-r)*np.exp(-r/2)
tp_2s = lambda theta, phi: sph_harm(0, 0, theta, phi).real
psi_2s = lambda r, theta, phi: r_2s(r) * tp_2s(theta, phi)


    

In [3]:

    

def generateScatterRandomized(psi_func, r_min = 0.1, r_max = 8, r_number = 50,
                  filename = None, silent = False, threshold = 100,
                  display = True):
    global x, y, z, r, theta, phi, psi
    if not silent:
        start_time = time.time() 
    r_vals = np.linspace(r_min, r_max, r_number)
    r_step = r_vals[1]-r_vals[0]   
    r = []
    phi = []
    theta = []
    psi = []
    i = 0
    for rval in r_vals:
        diameter = rval * 2 * np.pi
        num = diameter / r_step
        phi_vals = np.linspace(0, 1 * np.pi, num / 2)
        try:
            phi_step = phi_vals[1]-phi_vals[0]
        except:
            phi_step = 0
        for phival in phi_vals:
            theta_vals = np.linspace(0, 2 * np.pi, np.sin(phival) * num)
            try:
                theta_step = theta_vals[1]-theta_vals[0]
            except:
                theta_step = 0
                theta_vals = [0]
            for thetaval in theta_vals:  
                thetaval = thetaval + (np.random.rand() - 0.5) * 0.2 * theta_step
                rval = rval + (np.random.rand() - 0.5) * 0.2 * r_step
                phival = phival + (np.random.rand() - 0.5) * 0.2 * phi_step
                psival = psi_func(rval, phival, thetaval)
                i = i + 1
                if np.random.rand() < psival * np.conj(psival) * threshold:
                    psi.append(np.real(psival))
                    r.append(rval)
                    phi.append(phival)
                    theta.append(thetaval)
    if not silent:
        generate_time = time.time()
        print("Time to Generate Points: {}s".format(generate_time-start_time))
    x = r * sin(phi) * cos(theta)
    y = r * sin(phi) * sin(theta)
    z = r * cos(phi)
    figure = mlab.figure('DensityPlot', bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(600, 600))
    mlab.clf()
    mlab.points3d(x, y, z, psi, scale_mode='none', scale_factor=0.1)
    mlab.view(0, 90, 50, (0, 0, 0)) 
    if not filename == None:
        mlab.savefig(filename)  
    if display:
        # This will allow you to interact with the figure, but stalls the program
        # until you close the window.
        # If display is False, the window will come up, but it acts like it's frozen.
        # That's fine. Just don't try to close it and it will work with each new frame
        mlab.show()
    if not silent:
        print("Total points generated: {}".format(i))
        print("Total points displayed: {}".format(len(r))) 
        print("Total time to generate image: {}s".format((time.time() - start_time)))


    

In [4]:

    

def generateScatter(psi_func, r_min = 0.1, r_max = 10, r_number = 30,
                  filename = None, silent = False, threshold = 100,
                  display = True):
    global x, y, z, r, theta, phi, psi
    if not silent:
        start_time = time.time() 
    r_vals = np.linspace(r_min, r_max, r_number)
    r_step = r_vals[1]-r_vals[0]   
    r = []
    phi = []
    theta = []
    psi = []
    i = 0
    for rval in r_vals:
        diameter = rval * 2 * np.pi
        num = diameter / r_step
        phi_vals = np.linspace(0, 1 * np.pi, num / 2)        
        for phival in phi_vals:
            theta_vals = np.linspace(0, 2 * np.pi, np.sin(phival) * num)
            if phival == 0:
                theta_vals = [0]
            theta_offset = np.random.rand()
            for thetaval in theta_vals:
                psival = psi_func(rval, phival, thetaval)
                i = i + 1
                if np.random.rand() < psival * np.conj(psival) * threshold:
                    psi.append(np.real(psival))
                    r.append(rval)
                    phi.append(phival)
                    theta.append(thetaval+theta_offset)
    if not silent:
        generate_time = time.time()
        print("Time to Generate Points: {}s".format(generate_time-start_time))
    x = r * sin(phi) * cos(theta)
    y = r * sin(phi) * sin(theta)
    z = r * cos(phi)
    figure = mlab.figure('DensityPlot', bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(600, 600))
    mlab.clf()
    mlab.points3d(x, y, z, psi, scale_mode='none', scale_factor=0.05)
    mlab.view(5, 85, 50, (0, 0, 0)) 
    if not filename == None:
        mlab.savefig(filename)  
    if display:
        # This will allow you to interact with the figure, but stalls the program
        # until you close the window.
        # If display is False, the window will come up, but it acts like it's frozen.
        # That's fine. Just don't try to close it and it will work with each new frame
        mlab.show()
    if not silent:
        print("Total points generated: {}".format(i))
        print("Total points displayed: {}".format(len(r))) 
        print("Total time to generate image: {}s".format((time.time() - start_time)))


    

In [5]:

    

def generateScatterSlice(psi_func, r_min = 0.1, r_max = 8, r_number = 30,
                  filename = None, silent = False, threshold = 100,
                  display = True):
    global x, y, z, r, theta, phi, psi
    if not silent:
        start_time = time.time() 
    r_vals = np.linspace(r_min, r_max, r_number)
    r_step = r_vals[1]-r_vals[0]   
    r = []
    phi = []
    theta = []
    psi = []
    i = 0
    for rval in r_vals:
        diameter = rval * 2 * np.pi
        num = diameter / r_step
        phi_vals = np.linspace(0, 1 * np.pi, num / 2)        
        for phival in phi_vals:
            theta_vals = np.array([0,np.pi])                   
            for thetaval in theta_vals:
                psival = psi_func(rval, phival, thetaval)
                i = i + 1
                if np.random.rand() < psival * np.conj(psival) * threshold:
                    psi.append(np.real(psival))
                    r.append(rval)
                    phi.append(phival)
                    theta.append(thetaval)
    if not silent:
        generate_time = time.time()
        print("Time to Generate Points: {}s".format(generate_time-start_time))
    x = r * sin(phi) * cos(theta)
    y = r * sin(phi) * sin(theta)
    z = r * cos(phi)
    figure = mlab.figure('DensityPlot', bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(600, 600))
    mlab.clf()
    mlab.points3d(x, y, z, psi, scale_mode='none', scale_factor=0.05)
    mlab.view(5, 85, 50, (0, 0, 0)) 
    if not filename == None:
        mlab.savefig(filename)  
    if display:
        # This will allow you to interact with the figure, but stalls the program
        # until you close the window.
        # If display is False, the window will come up, but it acts like it's frozen.
        # That's fine. Just don't try to close it and it will work with each new frame
        mlab.show()
    if not silent:
        print("Total points generated: {}".format(i))
        print("Total points displayed: {}".format(len(r))) 
        print("Total time to generate image: {}s".format((time.time() - start_time)))


    

In [6]:

    

def generateIso(psi_func, r_min = 0.2, r_max = 5, r_number = 40,
                  filename = None, silent = False, threshold = 0.0015,
                  display = True):
    global x, y, z, psi
    if not silent:
        start_time = time.time()
    r_vals = np.linspace(r_min, r_max, r_number)
    r_range = r_max - r_min
    step_size = r_range / r_number
    r = []
    phi = []
    theta = []
    psi = []
    i = 0
    for rval in r_vals:
        diameter = rval * 2 * np.pi
        num = diameter / step_size
        phi_vals = np.linspace(0, 1 * np.pi, num / 2)
        for phival in phi_vals:
            theta_vals = np.linspace(0, 2 * np.pi, np.sin(phival) * num)
            if phival == 0:
                theta_vals = [0]
            theta_offset = np.random.rand()
            for thetaval in theta_vals:  
                psival = psi_func(rval, phival, thetaval)
                i = i + 1
                if psival * psival >= 0.8 * threshold:# and psival * np.conj(psival) <= 2 * threshold :
                    psi.append(np.real(psival))
                    r.append(rval)
                    phi.append(phival)
                    theta.append(thetaval + theta_offset)
    if not silent:
        generate_time = time.time()
        print("Time to Generate Points: {}s".format(generate_time-start_time))
    x = r * sin(phi) * cos(theta)
    y = r * sin(phi) * sin(theta)
    z = r * cos(phi)
    psi = psi
    figure = mlab.figure('DensityPlot', bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(600, 600))
    mlab.clf()
    mlab.points3d(x, y, z, psi, scale_mode='none', scale_factor=1.0)
    #mlab.mesh(x, y, z, scalars = psi)
    mlab.view(0, 90, 50, (0, 0, 0)) 
    if not filename == None:
        mlab.savefig(filename)  
    if display:
        # This will allow you to interact with the figure, but stalls the program
        # until you close the window.
        # If display is False, the window will come up, but it acts like it's frozen.
        # That's fine. Just don't try to close it and it will work with each new frame
        mlab.show()
    if not silent:
        print("Total points generated: {}".format(i))
        print("Total points displayed: {}".format(len(r))) 
        print("Total time to generate image: {}s".format((time.time() - start_time)))


    

In [15]:

    

x = [[1,0,0,0],[1,1,1,0]]
y = [[0,0,0,1],[0,0,1,1]]
z = [[0,0,1,1],[0,1,1,1]]
figure = mlab.figure('DensityPlot', bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(600, 600))
mlab.clf()
mlab.mesh(x, y, z)
mlab.axes()
mlab.view(0, 90, 10, (0, 0, 0)) 
mlab.show()


    

In [14]:

    

start_time = time.time()
for i, t in enumerate(np.linspace(0,np.pi/2,50)):    
    my_psi = lambda r, phi, theta: (np.sin(t)*psi_2pz(r,phi,theta)*np.exp(1j*21*t)+np.cos(t)*psi_2s(r,phi,theta)**np.exp(1j*1*t))
    generateScatter(my_psi, r_number = 45, filename = "image{}.png".format(i), silent = True,
                    threshold = 30 + 125 * np.sin(t) * np.sin(t), display = False)
print("Total Time to Animate: {}s".format(time.time()-start_time))


    

    




Total Time to Animate: 7294.11182594s

In [ ]:

    

start_time = time.time()
for i, t in enumerate(np.linspace(0,np.pi/2,50)):    
    my_psi = lambda r, phi, theta: (np.sin(t)*psi_2pz(r,phi,theta)*np.exp(1j*21*t)+np.cos(t)*psi_2s(r,phi,theta)**np.exp(1j*1*t))
    generateIso(my_psi, r_number = 40, filename = "Iso_image{}.png".format(i), silent = True, display = False, threshold = 0.0005)
print("Total Time to Animate: {}s".format(time.time()-start_time))


    

In [47]:

    

generateIso(psi_2pz, r_number = 20)
print(len(x))
print(np.max(psi))
print(np.max(x))


    

    




Time to Generate Points: 5.99455118179s
Total points generated: 37118
Total points displayed: 10383
Total time to generate image: 21.8142981529s
10383
0.0543487519161
3.08457609554

In [64]:

    

#print(theta)
print(np.min(phi))


    

    




-0.0151497540467

In [28]:

    

generateScatter(psi_2pz, r_number = 20, threshold = 100, display = True)
print(np.max(psi)*np.max(psi)*100)


    

    




Time to Generate Points: 7.38190317154s
Total points generated: 30584
Total points displayed: 602
Total time to generate image: 14.5747830868s
0.53419345367

In [144]:

    

pi = np.pi
cos = np.cos
sin = np.sin
dphi, dtheta = pi / 4, pi / 4
[phi, theta] = np.mgrid[0:pi + dphi * 1.5:dphi,
                           0:2 * pi + dtheta * 1.5:dtheta]
m0 = 4
m1 = 3
m2 = 2
m3 = 3
m4 = 6
m5 = 2
m6 = 6
m7 = 4
r = sin(m0 * phi) ** m1 + cos(m2 * phi) ** m3 + \
    sin(m4 * theta) ** m5 + cos(m6 * theta) ** m7
x = r * sin(phi) * cos(theta)
y = r * cos(phi)
z = r * sin(phi) * sin(theta)

mlab.mesh(x, y, z, colormap="bone")


    

    
Out[144]:





<mayavi.modules.surface.Surface at 0x7f57587a31d0>

In [148]:

    

print(z)


    

    




[[  0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00
    0.00000000e+00  -0.00000000e+00  -0.00000000e+00  -0.00000000e+00
   -0.00000000e+00   0.00000000e+00]
 [  0.00000000e+00   5.00000000e-01   7.07106781e-01   5.00000000e-01
    8.65956056e-17  -5.00000000e-01  -7.07106781e-01  -5.00000000e-01
   -1.73191211e-16   5.00000000e-01]
 [  0.00000000e+00   8.05180865e-64   0.00000000e+00   6.52196500e-62
    0.00000000e+00  -3.72922161e-59  -0.00000000e+00  -4.09224917e-62
   -0.00000000e+00   9.78304585e-59]
 [  0.00000000e+00   5.00000000e-01   7.07106781e-01   5.00000000e-01
    8.65956056e-17  -5.00000000e-01  -7.07106781e-01  -5.00000000e-01
   -1.73191211e-16   5.00000000e-01]
 [  0.00000000e+00   1.73191211e-16   2.44929360e-16   1.73191211e-16
    2.99951957e-32  -1.73191211e-16  -2.44929360e-16  -1.73191211e-16
   -5.99903913e-32   1.73191211e-16]
 [ -0.00000000e+00  -5.00000000e-01  -7.07106781e-01  -5.00000000e-01
   -8.65956056e-17   5.00000000e-01   7.07106781e-01   5.00000000e-01
    1.73191211e-16  -5.00000000e-01]]

In [41]:

    

start_time = time.time()
r_vals = np.linspace(0.1,5,40)
r = []
phi = []
theta = []
psi = []
i=0
for rval in r_vals:
    num = rval*51
    phi_vals = np.linspace(0,1*np.pi,num/2)
    for phival in phi_vals:
        theta_vals = np.linspace(0,2*np.pi,np.sin(phival)*num)
        for thetaval in theta_vals:
            psival = psi_2pz(rval, phival, thetaval)
            i = i+1
            if np.random.rand() < psival * psival * 30:#True:#
                psi.append(psival)
                r.append(rval)
                phi.append(phival)
                theta.append(thetaval + np.random.rand())
r = np.array(r)
phi = np.array(phi)
theta = np.array(theta)
psi = np.array(psi)
print("Total points evaluated: {}".format(i))
print("Total points displayed: {}".format(len(r)))
x = r * sin(phi) * cos(theta)
y = r * sin(phi) * sin(theta)
z = r * cos(phi)
#density = np.ones_like(z)
figure = mlab.figure('DensityPlot', bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(600, 600))
#mlab.clf()
mlab.points3d(x, y, z, psi, scale_mode='none', scale_factor=0.07)
#mlab.axes()
mlab.view(0, 90, 30, (0, 0, 0)) 
print("Time to generate image: {}s".format((time.time() - start_time)))
mlab.show()
# Without np.array() - 49.3 - 46.67 - 49.3685
# With np.array() - 49.95 - 49.98 - 49.55


    

    




Total points evaluated: 278843
Total points displayed: 8844
Time to generate image: 51.3174600601s

In [ ]:

    

r, phi, theta = np.mgrid[0:20:100j,0:pi:101j, 0:2 * pi:101j]

phi_vals = np.linspace(0,np.pi,100)
theta_vals = np.linspace(0,2*np.pi,100)
r_vals = np.empty(100)
rstep = lambda r: 1/(r*r)
step = 0.01
value = 0
for i in range(len(r_vals)):
    r_vals[i] = value
r_steps = np.linspace(0,20,100)
# The strategy is to first define the psi functions. Then generate
# points randomly with the proper density distribution. To do that
# we evaluate the density distribution at discrete points then
# create that number of new points randomly generated about the
# parent point. Next, calculate the value of psi at the new points.
# Finally, convert the new point coordinates into x,y,z space and
# plot them shaded by the value of psi.

# 2s functions
r_2s = lambda r: (2-r)*np.exp(-r/2)
tp_2s = lambda theta, phi: sph_harm(0, 0, theta, phi).real
psi_2s = lambda r, theta, phi: r_2s(r) * tp_2s(theta, phi)

# 2pz functions
r_2pz = lambda r: 1/np.sqrt(24)*r*np.exp(-r/2)
tp_2pz = lambda theta, phi: sph_harm(1, 0, theta, phi).real
psi_2pz = lambda r, theta, phi: r_2pz(r) * tp_2pz(theta, phi)
#psi_2pz = lambda r,theta: 1/(4*np.sqrt(np.pi))r*np.exp(-r/2)*np.cos(theta)

# Now for the density functions

dr_2s = lambda r: psi_2s(r)*psi_2s(r)*r*r*4*np.pi
def Orbital2s(r, phi, theta):
    d = 1 / (4 * np.sqrt(2*np.pi))*(2- r)*np.exp(-1*r/2)
    return d
def Orbital2pz(r, phi, theta):
    d = 1 / (4 * np.sqrt(2*np.pi))*r*np.exp(-1*r/2)*np.cos(theta)
    return d
orbital2s = Orbital2s(r,phi,theta)
x = r * sin(phi) * cos(theta)
y = r * sin(phi) * sin(theta)
z = r * cos(phi)
#mlab.figure(1, bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(400, 300))
#mlab.clf()
#pts = mlab.points3d(x, y, z, orbital2s, scale_mode='none', scale_factor=0.07)
#mlab.axes()
#mlab.show()


    

In [137]:

    

mlab.figure(1, bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(400, 300))
mlab.clf()
# Represent spherical harmonics on the surface of the sphere
for n in range(0, 2):
    for m in range(n):
        s = sph_harm(m, n, theta, phi).real

        mlab.mesh(x - m, y - n, z, scalars=s, colormap='jet')

        s[s < 0] *= 0.97

        s /= s.max()
        mlab.mesh(s * x - m, s * y - n, s * z + 1.3,
                  scalars=s, colormap='Spectral')

mlab.view(90, 70, 6.2, (-1.3, -2.9, 0.25))
mlab.show()


    

    




---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-137-ac9f11d3eafd> in <module>()
      6         s = sph_harm(m, n, theta, phi).real
      7 
----> 8         mlab.mesh(x - m, y - n, z, scalars=s, colormap='jet')
      9 
     10         s[s < 0] *= 0.97

/usr/lib/python2.7/dist-packages/mayavi/tools/helper_functions.pyc in the_function(*args, **kwargs)
     32 def document_pipeline(pipeline):
     33     def the_function(*args, **kwargs):
---> 34         return pipeline(*args, **kwargs)
     35 
     36     if hasattr(pipeline, 'doc'):

/usr/lib/python2.7/dist-packages/mayavi/tools/helper_functions.pyc in __call__(self, *args, **kwargs)
     77             scene.disable_render = True
     78         # Then call the real logic
---> 79         output = self.__call_internal__(*args, **kwargs)
     80         # And re-enable the rendering, if needed.
     81         if scene is not None:

/usr/lib/python2.7/dist-packages/mayavi/tools/helper_functions.pyc in __call_internal__(self, *args, **kwargs)
    830         filters.
    831         """
--> 832         self.source = self._source_function(*args, **kwargs)
    833         kwargs.pop('name', None)
    834         self.store_kwargs(kwargs)

/usr/lib/python2.7/dist-packages/mayavi/tools/sources.pyc in grid_source(x, y, z, **kwargs)
   1259     x, y, z, scalars = convert_to_arrays((x, y, z, scalars))
   1260     data_source = MGridSource()
-> 1261     data_source.reset(x=x, y=y, z=z, scalars=scalars)
   1262 
   1263     name = kwargs.pop('name', 'GridSource')

/usr/lib/python2.7/dist-packages/mayavi/tools/sources.pyc in reset(self, **traits)
    687         x, y, z = self.x, self.y, self.z
    688 
--> 689         assert len(x.shape) == 2, "Array x must be 2 dimensional."
    690         assert len(y.shape) == 2, "Array y must be 2 dimensional."
    691         assert len(z.shape) == 2, "Array z must be 2 dimensional."

AssertionError: Array x must be 2 dimensional.

In [2]:

    

s = sph_harm(0, 0, theta, phi).real
abss = np.abs(s)
s_orbital = [abss*x,abss*y,abss*z, s]


    

In [3]:

    

p = sph_harm(0, 1, theta, phi).real
absp = np.abs(p)
p_orbital = [absp*x, absp*y, absp*z, p]


    

In [4]:

    

d = sph_harm(1,2,theta,phi).real
absd = np.abs(d)
d_orbital = [absd*x, absd*y, absd*z, d]


    

In [13]:

    

def superposition(time):
    s_coeff = np.sin(time)*np.sin(time)
    p_coeff = 1-s_coeff
    super_orbital = [1,1,1,1]
    for i in range(len(super_orbital)):
        value = s_coeff*p_orbital[i]+p_coeff*d_orbital[i]
        super_orbital[i] = 16/np.sqrt(2)* value * value
    return super_orbital


    

In [8]:

    

mlab.figure(1, bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(400, 300))
mlab.clf()
mlab.mesh(super_orbital[0], super_orbital[1], super_orbital[2], scalars = super_orbital[3])
mlab.view(90, 70, 1, (0, 0, 0))
mlab.show()


    

    




---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-8-4fa6478b2938> in <module>()
      1 mlab.figure(1, bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(400, 300))
      2 mlab.clf()
----> 3 mlab.mesh(super_orbital[0], super_orbital[1], super_orbital[2], scalars = super_orbital[3])
      4 mlab.view(90, 70, 1, (0, 0, 0))
      5 mlab.show()

NameError: name 'super_orbital' is not defined

In [40]:

    

print(super_orbital)


    

    




[]

In [14]:

    

my_orbital = superposition(0)
mlab.figure(1, bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(400, 300))
mlab.clf()
fig = mlab.mesh(my_orbital[0], my_orbital[1], my_orbital[2], scalars = my_orbital[3])
mlab.view(90, 70, 1, (0, 0, 0))
fig_source = fig.mlab_source
for time in np.arange(0,6,0.1):
    my_orbital = superposition(time)
    fig_source.set(x=my_orbital[0], y=my_orbital[1], z=my_orbital[2], scalars = my_orbital[3])
    mlab.savefig("Time{:1f}.png".format(time))
#mlab.clf()


    

In [23]:

    

import numpy as np
from scipy import stats
from mayavi import mlab

mu, sigma = 0, 0.1 
x = 10*np.random.normal(mu, sigma, 5000)
y = 10*np.random.normal(mu, sigma, 5000)
z = 10*np.random.normal(mu, sigma, 5000)

xyz = np.vstack([x,y,z])
kde = stats.gaussian_kde(xyz)
density = kde(xyz)

# Plot scatter with mayavi
figure = mlab.figure('DensityPlot')
pts = mlab.points3d(x, y, z, density, scale_mode='none', scale_factor=0.07)
mlab.axes()
mlab.show()


    

In [24]:

    

# Plot scatter with mayavi
figure = mlab.figure('DensityPlot', bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(400, 300))
figure.scene.disable_render = True

pts = mlab.points3d(x, y, z, density, scale_mode='none', scale_factor=0.07) 
mask = pts.glyph.mask_points
mask.maximum_number_of_points = x.size
mask.on_ratio = 1
pts.glyph.mask_input_points = True

figure.scene.disable_render = False 
mlab.axes()
mlab.show()


    

In [22]:

    

mu, sigma = 0, 0.1 
x = 10*np.random.normal(mu, sigma, 5000)
y = 10*np.random.normal(mu, sigma, 5000)    
z = 10*np.random.normal(mu, sigma, 5000)

xyz = np.vstack([x,y,z])
kde = stats.gaussian_kde(xyz)

# Evaluate kde on a grid
xmin, ymin, zmin = x.min(), y.min(), z.min()
xmax, ymax, zmax = x.max(), y.max(), z.max()
xi, yi, zi = np.mgrid[xmin:xmax:30j, ymin:ymax:30j, zmin:zmax:30j]
coords = np.vstack([item.ravel() for item in [xi, yi, zi]]) 
density = kde(coords).reshape(xi.shape)

# Plot scatter with mayavi
figure = mlab.figure('DensityPlot', bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(400, 300))

grid = mlab.pipeline.scalar_field(xi, yi, zi, density)
min = density.min()
max=density.max()
mlab.pipeline.volume(grid, vmin=min, vmax=min + .5*(max-min))

mlab.axes()
mlab.show()


    

In [ ]:

    

import multiprocessing

def calc_kde(data):
    return kde(data.T)

xyz = np.vstack([x,y,z])
kde = stats.gaussian_kde(xyz)

# Evaluate kde on a grid
xmin, ymin, zmin = x.min(), y.min(), z.min()
xmax, ymax, zmax = x.max(), y.max(), z.max()
xi, yi, zi = np.mgrid[xmin:xmax:30j, ymin:ymax:30j, zmin:zmax:30j]
coords = np.vstack([item.ravel() for item in [xi, yi, zi]]) 

# Multiprocessing
cores = multiprocessing.cpu_count()
pool = multiprocessing.Pool(processes=cores)
results = pool.map(calc_kde, np.array_split(coords.T, 2))
density = np.concatenate(results).reshape(xi.shape)

# Plot scatter with mayavi
figure = mlab.figure('DensityPlot')

grid = mlab.pipeline.scalar_field(xi, yi, zi, density)
min = density.min()
max=density.max()
mlab.pipeline.volume(grid, vmin=min, vmax=min + .5*(max-min))

mlab.axes()
mlab.show()


    

In [ ]: