如何在 Python 中最大限度地分离立方体中的 n 个点?如何使用 scipy.integrate.ode?

How do I maximally separate n points in a cube in Python? How do I use scipy.integrate.ode?

我需要在单位立方体中找到 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]

Poisson Disk Sampling 是执行此操作的线性时间算法。参考文献:http://devmag.org.za/2009/05/03/poisson-disk-sampling/, http://bost.ocks.org/mike/algorithms/

布里德森,罗伯特。 "Fast Poisson disk sampling in arbitrary dimensions." ACM SIGGRAPH。卷。 2007. 2007.