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

2024-06-02 12:32:11 发布

您现在位置:Python中文网/ 问答频道 /正文

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

Tags: sizeclipreturndefnpdimensionsdeltaprint