我需要在单位立方体中找到n个点,它们的间距大致最大。我目前正在考虑将这个问题写成一个ODE,并使用scipy.integrate.ode
来完成这项工作。也就是说,做一个模拟,每个粒子沿着立方体的壁相互排斥。我不太在乎效率。你知道吗
不幸的是,我不能让它在大于n=5的问题上工作:
import numpy as np
def spread_points_in_cube(n, dimensions=3, rng=None):
from scipy.integrate import ode
if rng is None:
rng = np.random
size = n * dimensions
y0 = np.zeros((2 * size))
y0[:size] = rng.uniform(size=size)
t0 = 0.0
def clip_to_wall(positions, forces):
forces = np.where(positions == 0.0,
np.clip(forces, 0, np.inf),
forces)
forces = np.where(positions == 1.0,
np.clip(forces, -np.inf, 0.0),
forces)
return forces
def decode(y):
positions = np.clip(y[:size].reshape((n, dimensions)), 0, 1)
velocities = clip_to_wall(positions, y[size:].reshape((n, dimensions)))
return positions, velocities
def f(t, y):
retval = np.zeros((2 * size))
positions, velocities = decode(y)
#print("pos", positions)
delta_positions = positions[:, np.newaxis, :] - positions[np.newaxis, :, :]
# print("delta_positions", delta_positions)
distances = np.linalg.norm(delta_positions, axis=2, ord=2)
distances += 1e-5
# print("distances", distances)
pairwise_forces = delta_positions * (distances ** -3)[:, :, np.newaxis]
# print("pairwise f", pairwise_forces[0, 1])
forces = np.sum(pairwise_forces, axis=1)
forces -= 0.9 * velocities
forces = clip_to_wall(positions, forces)
#print("forces", forces)
retval[:size] = velocities.reshape(size)
retval[size:] = forces.reshape(size)
return retval
r = ode(f).set_integrator('vode', method='adams')
r.set_initial_value(y0, t0)
t_max = 40000
dt = 1
while r.successful() and r.t < t_max:
r.integrate(r.t + dt)
return decode(r.y)[0]
泊松圆盘采样是一种线性时间算法。参考文献:http://devmag.org.za/2009/05/03/poisson-disk-sampling/,http://bost.ocks.org/mike/algorithms/
布里森,罗伯特。”任意维的快速泊松圆盘取样。第2007卷。2007你知道吗
相关问题 更多 >
编程相关推荐