How can I set a minimum distance constraint for generating points with numpy.random.rand?

Question

I am trying to generate an efficient code for generating a number of random position vectors which I then use to calculate a pair correlation function. I am wondering if there is straightforward way to set a constraint on the minimum distance allowed between any two points placed in my box.

My code currently is as follows:

def pointRun(number, dr):
"""
Compute the 3D pair correlation function
for a random distribution of 'number' particles
placed into a 1.0x1.0x1.0 box.
"""
## Create array of distances over which to calculate.   
    r = np.arange(0., 1.0+dr, dr)

## Generate list of arrays to define the positions of all points,
##    and calculate number density.
    a = np.random.rand(number, 3)
    numberDensity = len(a)/1.0**3

## Find reference points within desired region to avoid edge effects. 
    b = [s for s in a if all(s > 0.4) and all(s < 0.6) ]

## Compute pairwise correlation for each reference particle
    dist = scipy.spatial.distance.cdist(a, b, 'euclidean')
    allDists = dist[(dist < np.sqrt(3))]

## Create histogram to generate radial distribution function, (RDF) or R(r)
    Rr, bins = np.histogram(allDists, bins=r, density=False)

## Make empty containers to hold radii and pair density values.
    radii = []
    rhor = []

## Normalize RDF values by distance and shell volume to get pair density.
    for i in range(len(Rr)):
        y = (r[i] + r[i+1])/2.
        radii.append(y)
        x = np.average(Rr[i])/(4./3.*np.pi*(r[i+1]**3 - r[i]**3))
        rhor.append(x)

## Generate normalized pair density function, by total number density
    gr = np.divide(rhor, numberDensity)
    return radii, gr

I have previously tried using a loop that calculated all distances for each point as it was made and then accepted or rejected. This method was very slow if I use a lot of points.

Spellbind

You might be able to use a variant of this algorithm, but the caveat is that it has a minimum resolution (and the higher it is the more costly it is). – Carilyn

Answer 1

Here is a scalable O(n) solution using numpy. It works by initially specifying an equidistant grid of points and then perturbing the points by some amount keeping the distance between the points at most min_dist.

You'll want to tweak the number of points, box shape and perturbation sensitivity to get the min_dist you want.

Note: If you fix the size of a box and specify a minimum distance between every point, it makes sense that there will be a limit to the number of points you can draw satisfying the minimum distance.

import numpy as np
import matplotlib.pyplot as plt

# specify params
n = 500
shape = np.array([64, 64])
sensitivity = 0.8 # 0 means no movement, 1 means max distance is init_dist

# compute grid shape based on number of points
width_ratio = shape[1] / shape[0]
num_y = np.int32(np.sqrt(n / width_ratio)) + 1
num_x = np.int32(n / num_y) + 1

# create regularly spaced neurons
x = np.linspace(0., shape[1]-1, num_x, dtype=np.float32)
y = np.linspace(0., shape[0]-1, num_y, dtype=np.float32)
coords = np.stack(np.meshgrid(x, y), -1).reshape(-1,2)

# compute spacing
init_dist = np.min((x[1]-x[0], y[1]-y[0]))
min_dist = init_dist * (1 - sensitivity)

assert init_dist >= min_dist
print(min_dist)

# perturb points
max_movement = (init_dist - min_dist)/2
noise = np.random.uniform(
    low=-max_movement,
    high=max_movement,
    size=(len(coords), 2))
coords += noise

# plot
plt.figure(figsize=(10*width_ratio,10))
plt.scatter(coords[:,0], coords[:,1], s=3)
plt.show()

Alienate

Answer 2

Based on @Samir 's answer, and make it a callable function for your convenience :)

import numpy as np
import matplotlib.pyplot as plt

def generate_points_with_min_distance(n, shape, min_dist):
    # compute grid shape based on number of points
    width_ratio = shape[1] / shape[0]
    num_y = np.int32(np.sqrt(n / width_ratio)) + 1
    num_x = np.int32(n / num_y) + 1

    # create regularly spaced neurons
    x = np.linspace(0., shape[1]-1, num_x, dtype=np.float32)
    y = np.linspace(0., shape[0]-1, num_y, dtype=np.float32)
    coords = np.stack(np.meshgrid(x, y), -1).reshape(-1,2)

    # compute spacing
    init_dist = np.min((x[1]-x[0], y[1]-y[0]))

    # perturb points
    max_movement = (init_dist - min_dist)/2
    noise = np.random.uniform(low=-max_movement,
                                high=max_movement,
                                size=(len(coords), 2))
    coords += noise

    return coords

coords = generate_points_with_min_distance(n=8, shape=(2448,2448), min_dist=256)

# plot
plt.figure(figsize=(10,10))
plt.scatter(coords[:,0], coords[:,1], s=3)
plt.show()

Farmergeneral

Answer 3

As I understood, you're looking for an algorithm to create many random points in a box such that no two points are closer than some minimum distance. If this is your problem, then you can take advantage of statistical physics, and solve it using molecular dynamics software. Moreover, you do need molecular dynamics or Monte Carlo to obtain exact solution of this problem.

You place N atoms in a rectangular box, create a repulsive interaction of a fixed radius between them (such as shifted Lennard-Jones interaction), and run simulation for some time (untill you see that the points spread out uniformly throughout the box). By laws of statistical physics you can show that positions of the points would be maximally random given the constraint that points cannot be close than some distance. This would not be true if you use iterative algorithm, such as placing points one-by-one and rejecting them if they overlap

I would estimate a runtime of several seconds for 10000 points, and several minutes for 100k. I use OpenMM for all my moelcular dynamics simulations.

Nomanomad

Answer 4

#example of generating 50 points in a square of 4000x4000 and with minimum distance of 400

import numpy as np
import random as rnd

n_points=50
x,y = np.zeros(n_points),np.zeros(n_points)
x[0],y[0]=np.round(rnd.uniform(0,4000)),np.round(rnd.uniform(0,4000))
min_distances=[]
i=1
while i<n_points :
    x_temp,y_temp=np.round(rnd.uniform(0,4000)),np.round(rnd.uniform(0,4000))
    distances = []
    for j in range(0,i):
        distances.append(np.sqrt((x_temp-x[j])**2+(y_temp-y[j])**2))
    min_distance = np.min(distances)
    if min_distance>400 :
        min_distances.append(min_distance)
        x[i]=x_temp
        y[i]=y_temp
        i = i+1

print(x,y)

Pretty