# Sampling from k circles with different radii

Description: I would like to sample n_points points from k circles whose radii and centers are given as lists.

Note: This is my first question :) And k is not a very large number, usually k<7.

My code: I guess the parts where I compute the number of points for each circle, and the for loop can be improved. Perhaps with a simple trick, the for loop can be removed, but I am not really sure. Any help to beautify the code is appreciated.

# given parameters
n_points_total = 250
centers = [[3,0],[7,0],[13,0]]

# my code starts here
import numpy as np

r_sum = sum(radii) # proportional to total arc length
probabilities = [r/r_sum for r in radii] # probability of a point being in the i'th circle.
positions = np.random.choice(len(radii), n_points_total, p=probabilities) # each point is randomly assigned to a circle
n_points = np.unique(positions, return_counts = True)[1] # counts how many points each circle will contain

data = []

for n, center, r in zip(n_points, centers, radii):
t = np.random.random(n)*2*np.pi
x = r*np.cos(t) + center[0]
y = r*np.sin(t) + center[1]
data.append(np.vstack([x,y]).T)

data = np.concatenate(data)

• Are you sure that this works the way you want? zip is a "zip-shortest", and your arguments do not have equal length. centers has an outer dimension of 3; radii has 3; n_points has fewer (depending on the n_points parameter) May 14 at 12:57
• Rather the n_points_total parameter. May 14 at 13:06

Let's assume that your zip-shortest behaviour is intended.

radii and centers should immediately be made np.ndarrays. If they're fed from some external source like a file or stdin, there are easy ways to get them into arrays.

import at the very top, not after your givens.

Everything here can be vectorised. Avoid like the plague for loops, lists and iterative vstack calls.

First, recognise that your centres and radii are going to be repeated by the respective counts represented by the elements of n_points. This is exactly what np.repeat does.

Rather than indexing into the result of np.unique(), prefer tuple-unpacking.

Do not use np.random.random; that's deprecated. Prefer the new Generator stuff and call into uniform instead of post-multiplying.

Add unit tests with a predictable random seed.

## Suggested

import numpy as np
from numpy.random import default_rng
from numpy.random._generator import Generator

def new(
centers: np.array,
rand: Generator,
n_points_total: int,
) -> np.ndarray:
probabilities = radii / radii.sum()  # probability of a point being in the i'th circle.
positions = rand.choice(a=len(radii), size=n_points_total, p=probabilities)  # each point is randomly assigned to a circle
_, n_points = np.unique(positions, return_counts=True)  # counts how many points each circle will contain

t = rand.uniform(low=0, high=2*np.pi, size=n_points_total)
np.cos(t, out=cossin[:, 0])
np.sin(t, out=cossin[:, 1])

return data

def old(
centers: np.array,
rand: Generator,
n_points_total: int,
) -> np.ndarray:
r_sum = sum(radii)  # proportional to total arc length
probabilities = [r / r_sum for r in radii]  # probability of a point being in the i'th circle.
p=probabilities)  # each point is randomly assigned to a circle
n_points = np.unique(positions, return_counts=True)[1]  # counts how many points each circle will contain

data = []

for n, center, r in zip(n_points, centers, radii):
t = rand.random(n) * 2 * np.pi
x = r * np.cos(t) + center[0]
y = r * np.sin(t) + center[1]
data.append(np.vstack([x, y]).T)

data = np.concatenate(data)
return data

def test() -> None:
centers = np.array((
(3, 0), (7, 0), (13, 0),
), dtype=np.float64)
radii = np.array((3, 1, 5), dtype=np.float64)

for method in (old, new):
rand = default_rng(seed=0)
assert np.allclose(
np.array((
(4.20619233, -2.74683455),
(5.99955592,  0.05161697),
(4.87431947, -2.34241895),
(7.97781686,  0.20946168),
(6.87251924, -0.99184104),
(7.45031626,  0.89286912),
(7.65268587, -0.75762864),
(6.03374088, -0.25757196),
(6.69270516,  0.95161435),
(6.11568453,  0.46688988),
)),
)

data = method(centers, radii, rand, n_points_total=250)
assert data.shape == (250, 2)
assert np.isclose(-4.999847197562797, data.min())
assert np.isclose(17.99763106087105, data.max())

if __name__ == '__main__':
test()


I don't have much to add to the answer of @Reinderien, maybe a slightly different solution.


import numpy as np

"""
usage:

generate_points_on_circles(
centers = np.array([[3, 7, 13],
[0, 0, 0]]),
num_pts = 10
)
"""
assert(radii.shape[0] == centers.shape[1] and centers.shape[0] == 2 and num_pts >= 0)


This does assume row vectors, but you can always transpose with a .T. To keep the points grouped by the circle they are from, add a np.sort call to the line defining the circle indices.