#All necesary import
import sys
import time
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
sys.path.append('../')
from ten.utils import generate_random_points_in_sphere
n = 3000
R = 5
theta = np.random.uniform(low = 0, high = 2*np.pi, size = n)
phi = np.random.uniform(low = 0, high = np.pi, size = n)
x = np.sin(phi)*np.cos(theta)
y = np.sin(phi)*np.sin(theta)
z = np.cos(phi)
fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z)
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fa775223ba8>
Note the large number of dots in the poles, there is no uniform distriución.
n = 3000
theta = 2*np.pi*np.random.uniform(size=n)
phi = np.arccos(2*np.random.uniform(size=n) - 1)
x = np.sin(phi)*np.cos(theta)
y = np.sin(phi)*np.sin(theta)
z = np.cos(phi)
fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z)
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fa7750cf4a8>
n = 3000
theta = 2*np.pi*np.random.uniform(size=n)
phi = np.arccos(2*np.random.uniform(size=n) - 1)
#theta = np.random.uniform(low=0, high=2*np.pi, size=n)
#phi = np.arcsin(2 * np.random.uniform(low=0, high=2*np.pi, size=n) - 1)
u = np.random.uniform(low=0, high=10, size=n)
x = np.sin(phi)*np.cos(theta)*u
y = np.sin(phi)*np.sin(theta)*u
z = np.cos(phi)*u
fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z)
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fa77510c4a8>
Obviously there is not a uniform distribution
n = 3000
R = 5
a = 0
x_l = []
y_l = []
z_l = []
t = time.time()
while a <= n:
x = np.random.uniform(low=-R, high=R, size=1)
y = np.random.uniform(low=-R, high=R, size=1)
z = np.random.uniform(low=-R, high=R, size=1)
if x*x + y*y + z*z <= R*R:
a += 1
x_l.append(x)
y_l.append(y)
z_l.append(z)
print("time =",time.time() - t)
fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x_l, y_l, z_l)
time = 0.30646824836730957
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fa774f50780>
This way it works but uses an if, which makes it extremely slow.
Thaks to Nate Eldredge to this solution in math.stackexchange
Let's say your sphere is centered at the origin $(0,0,0)$. For the distance $D$ from the origin of your random pointpoint, note that you want $P(D≤r)=(\frac{r}{R_s})^3$. Thus if $U$ is uniformly distributed between 0 and 1, taking $D=R_sU^{1/3}$ will do the trick. For the direction, a useful fact is that if $X_1$,$X_2$,$X_3$ are independent normal random variables with mean 0 and variance 1, then $$\frac{1}{\sqrt{X_1^2+X_2^2+X_3^2}}(X_1,X_2,X_3)$$ Is uniformly distributed on (the surface of) the unit sphere. You can generate normal random variables from uniform ones in various ways; the Box-Muller algorithm is a nice simple approach.
So then $$\frac{R_s U^{1/3}}{\sqrt{X_1^2+X_2^2+X_3^2}}(X_1,X_2,X_3)$$
would produce a uniformly distributed point inside the ball of radius $Rs$.
n_points = 3000
R = 5
r = 0
t = time.time()
U = np.random.random(n_points)
uniform_between_R_r = (R - r) * U**(1/3) + r
X = np.random.randn(n_points, 3)
randoms_versors = (np.sqrt(X[:, 0]**2 + X[:, 1]**2 + X[:, 2]**2))**(-1) * X.T
points = randoms_versors * uniform_between_R_r
print("time =", time.time() - t)
fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(points[0,:], points[1,:], points[2,:])
time = 0.003650188446044922
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fa774df46a0>
aceptors = list(range(1000, 10001, 1000)) + list(range(10000, 100000, 10000)) + list(range(100000, 300000, 100000))
slower = []
faster = []
R = 5
r = 0
for aceptor in aceptors:
t = time.time()
a = 0
while a <= aceptor:
x = np.random.uniform(low=-R, high=R, size=1)
y = np.random.uniform(low=-R, high=R, size=1)
z = np.random.uniform(low=-R, high=R, size=1)
if x*x + y*y + z*z <= R*R:
a += 1
x_l.append(x)
y_l.append(y)
z_l.append(z)
slower.append(time.time() - t)
t = time.time()
n_points = aceptor
U = np.random.random(n_points)
uniform_between_R_r = (R - r) * U**(1/3) + r
X = np.random.randn(n_points, 3)
randoms_versors = (np.sqrt(X[:, 0]**2 + X[:, 1]**2 + X[:, 2]**2))**(-1) * X.T
points = randoms_versors * uniform_between_R_r
faster.append(time.time() - t)
plt.figure(figsize=(10, 6))
plt.plot(aceptors, slower, aceptors, faster, linewidth=2)
plt.legend(['slower', 'faster'], loc=0)
plt.xlabel('Num. of acceptors')
plt.ylabel('Time')
plt.grid()
The question is: this code is quicker?
num_points = 300000
t = time.time()
points = generate_random_points_in_sphere(num_points, 5)
print("time =", time.time() - t)
time = 0.19737458229064941
The version with the if take 0.13 second with 3000 points.This vercion take 0.19 with 300000!! 2 orders of magnitude difference.
num_points = 3000
t = time.time()
points = generate_random_points_in_sphere(num_points, 5, 4)
print("time =", time.time() - t)
fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(points[:,0], points[:,1], points[:,2])
time = 0.0024046897888183594
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7f12940e4438>
num_points = 3000
points = generate_random_points_in_sphere(num_points, 5, 4)
mid_points = np.zeros_like(points)
for i in range(3000):
if points[i, 1] < 3 and points[i,1] > -3:
mid_points[i] = points[i]
fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(mid_points[:,0], mid_points[:,1], mid_points[:,2])
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7f12940e1160>
#Este css es trabajo de @LorenaABarba y su grupo
from IPython.core.display import HTML
css_file = 'css/personal.css'
HTML(open(css_file, "r").read())
El código esta licenciado bajo MIT.
La documentación bajo:
TEN by Laboratorio de Microscopia Óptica Avanzada - UNRC is licensed under a Creative Commons Attribution 4.0 International License.
Based on a work at https://github.com/pewen/ten