In [1]:
#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
In [3]:
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)
Out[3]:
Note the large number of dots in the poles, there is no uniform distriución.
In [5]:
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)
Out[5]:
In [6]:
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)
Out[6]:
In [7]:
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)
Out[7]:
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$.
In [15]:
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,:])
Out[15]:
In [54]:
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)
In [58]:
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?
In [7]:
num_points = 300000
t = time.time()
points = generate_random_points_in_sphere(num_points, 5)
print("time =", time.time() - t)
In [8]:
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])
Out[8]:
In [9]:
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])
Out[9]:
In [10]:
#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())
Out[10]:
El código esta licenciado bajo MIT.
La documentación bajo:
<span xmlns:dct="http://purl.org/dc/terms/" property="dct:title">TEN</span> by <span xmlns:cc="http://creativecommons.org/ns#" property="cc:attributionName">Laboratorio de Microscopia Óptica Avanzada - UNRC</span> is licensed under a Creative Commons Attribution 4.0 International License.
Based on a work at <a xmlns:dct="http://purl.org/dc/terms/" href="https://github.com/pewen/ten" rel="dct:source">https://github.com/pewen/ten</a>