Area-based algorithm, \$O(k^2+n\log(k))\$
Our general strategy is to instead select points (outside the sticker rectangles) from the square (-1,-1) to (1,1), then repeat if we don't get a point within the circle. On average, we will have to select \$\frac{4}{\pi}\$ points from the square for each point in the circle.
The runtime is not depend on rectangle coverage, but it does depend on measures of rectangle intersection (denoted L in the complexity analysis in the code). For the listed complexity, I use \$L=O(k)\$, which is probably not tight. If \$L=O(\log(k))\$ can be proven, then the complexity is \$O((k+n)\log(k))\$ (see get_heights for an algorithmic alternative).
Define "open area" to be the area that points could be generated in. Recall that we are using a square instead of a circle, so the total area is piece-wise linear in x.
Precompute the total open height function \$h(x)\$ (this is piece-wise constant, so store the change at each transition: either the start or end of a new rectangle). Similarly, pre-compute the cumulative areas \$\int_{-1}^{x} h(x)\$, which is piece-wise linear, so once again store the value at each transition. This can be done in \$O(k \log(k))\$ time
Precompute the intervals on which y coordinates can be distributed for a given x (again, only storing this interval on each transition). This might be able to be done in \$O(k \log(k))\$, but this algorithm only succeeds at \$O(k^2)\$ time
Generate points. For each point: (average \$\frac{4}{\pi} n\$ times for a total runtime of \$O(n\log(k))\$
a. Choose a random value \$r\$ for the area (\$O(1)\$)
b. Find \$x\$ such that \$\int_{-1}^{x} h(t) dt = r\$. This is accomplished through binary search on the pre-computed heights \$h\$ and cumulative areas \$\int h\$, so this is done in \$O(\log(k))\$ time
c. Choose a random value \$y\$ for the given x such that (x,y) is in an open area. Through similar binary search as 3b, this is done in \$O(\log(k))\$ time
d. Reject the point if it is outside the circle (constant factor 4/pi)
The total runtime is \$O(k\log(k)) + O(k^2) + O(n\log(k)) = O(k^2+n\log(k))\$
In the following example images, I used the same n=10000, so more rectangle coverage meant a higher density of points.
Full python code below, including matplotlib display:
import math
import random
# for binary search
import bisect
# matplotlib for plotting the output points
import matplotlib.pyplot as plt
import matplotlib
def find_le(ls, x):
# find_le is O(log(|ls|))
# Find rightmost value less than or equal to x in ls
i = bisect.bisect_right(ls, x)
return i-1, ls[i-1]
class Interval:
def __init__(self, lo, hi):
# each range is tuple (lo, hi)
# ranges shall be sorted in ascending order
# each range does not overlap with other ranges
self.ranges = [(lo, hi)]
def remove(self, lo, hi):
new_ranges = []
for self_range in self.ranges:
if lo <= self_range[0] <= hi < self_range[1]:
new_ranges.append((hi, self_range[1]))
elif self_range[0] < lo <= self_range[1] <= hi:
new_ranges.append((self_range[0], lo))
elif self_range[0] < lo < hi < self_range[1]:
new_ranges.append((self_range[0], lo))
new_ranges.append((hi, self_range[1]))
elif lo <= self_range[0] < self_range[1] <= hi:
pass
elif hi <= self_range[0] or self_range[1] <= lo:
# no overlap
new_ranges.append(self_range)
self.ranges = new_ranges
def cache_lengths(self):
# If L is the number of subintervals of this Interval
# then cache_lengths is O(L) ⊂ O(k)
# in cum_lengths, (0.6, 0.7) means that the total length to the
# left of x=0.7 in this interval is 0.6
self.cum_lengths = [(0, -1)]
last_length = 0
for self_range in self.ranges:
length = self_range[1] - self_range[0]
self.cum_lengths.append((last_length, self_range[0]))
last_length += length
self.total_length = last_length
def random_point(self):
# If L is the number of subintervals of this Interval
# then random_point is O(log(L)) ⊂ O(log(k))
r = random.random()*self.total_length
# This find_le is O(log(L))
i, (length_left, x_left) = find_le(self.cum_lengths, (r,))
return (r - length_left) + x_left
def __repr__(self):
return "I" + str(self.ranges)
def interval_complement(rects):
# O(L) where L is the length of rects list
interval = Interval(-1, 1)
for rect in rects:
interval.remove(rect.bottom, rect.top)
return interval
class Rect:
def __init__(self, left, right, bottom, top):
self.left = min(left, right)
self.right = max(left, right)
self.bottom = min(bottom, top)
self.top = max(bottom, top)
self.width = self.right - self.left
self.height = self.top - self.bottom
def __repr__(self):
return f"Rect({self.left}, {self.right}, {self.bottom}, {self.top})"
def overlap(self, others):
# returns vertical (y) overlap between self and other rects
total_overlap = 0
# https://pyinterval.readthedocs.io/en/latest/guide.html#operations is
# looking so tempting right now, but unions/merges would
# require Principle Inclusion-Exclusion which can be k^2 worst case
starts = sorted(others, key=lambda rect: rect.bottom, reverse=True)
ends = sorted(others, key=lambda rect: rect.top, reverse=True)
active_rects = set()
overlap_start = self.bottom
while ends:
if starts and starts[-1].bottom < ends[-1].top:
active_rect = starts.pop()
if not active_rects and active_rect.bottom > self.bottom:
# This is the start of an overlap interval
overlap_start = active_rect.bottom
active_rects.add(active_rect)
else:
active_rect = ends.pop()
active_rects.remove(active_rect)
if not active_rects:
# This is the end of an overlap interval
# Rare cases of negative = started and ended below start of self
total_overlap += max(0, min(active_rect.top, self.top) - overlap_start)
return total_overlap
def as_matplotlib_rect(self):
return matplotlib.patches.Rectangle((self.left, self.bottom), self.width, self.height)
def get_heights(rects):
# get_heights is O(k log(k)) total
# Sorting: O(k log(k))
# Use these lists as a queue
starts = sorted(rects, key=lambda rect: rect.left, reverse=True)
ends = sorted(rects, key=lambda rect: rect.right, reverse=True)
# list of tuples (open height starting at x, x, open intervals starting at x)
heights = [(2, -1, Interval(-1,1))]
active_rects = set()
while ends:
# inner loop performed once for k left edges and k right edges,
# so this is O(k)
if starts and starts[-1].left < ends[-1].right:
# next x-pos is the start of a new rect
rect = starts.pop()
height_diff = rect.height - rect.overlap(active_rects)
height = heights[-1][0] - height_diff
active_rects.add(rect)
# interval_complement is O(L) so this is O(k^2), not sure if it is O(k log(k))
heights.append((height, rect.left, interval_complement(active_rects)))
else:
# next x-pos is the end of an old rect
rect = ends.pop()
active_rects.remove(rect)
# rect.overlap is O(L) as well
height_diff = rect.height - rect.overlap(active_rects)
# Two choices to get O(k log(k)):
# 1. prove that |active_rects| is O(log(k)) under the given approaches
# 2. avoid using interval_complement on each transition using some increments approach, and avoid doing rect.overlap on each transition
height = heights[-1][0] + height_diff
heights.append((height, rect.right, interval_complement(active_rects)))
heights.append((2, 1, Interval(-1,1)))
# cum_areas is the cumulative areas from -1 to x (of not rects)
# Not subtracting curve because we will be sampling from square
cum_areas = [(0, -1, Interval(-1, 1))]
last_height = 2
for height, x, interval in heights[1:]:
last_area, last_x, _ = cum_areas[-1]
new_area = (x-last_x)*last_height
cum_areas.append((last_area + new_area, x, interval))
last_height = height
return (heights, cum_areas)
def random_pick_points(rects, num_points):
# random_pick_points is O((n+k) log(k))
# Step 1. Precompute heights: O(k log (k))
heights, cum_areas = get_heights(rects)
total_area = cum_areas[-1][0]
points = []
# Step 2. Precompute vertical intervals: O(k)
for _, _, y_interval in cum_areas:
y_interval.cache_lengths()
# Step 3. Sample points outside rects and within circle
# Total time complexity of this loop is O(n log(k))
while len(points) < num_points:
# This loop is reached average (4/pi)*n = O(n) times
# Each occurence has complexity O(log(k))
# Step 3a. Choose a random area value
area = random.random()*total_area
# Step 3b. Choose the x value such that area is the total area
# outside rects and within 2x2 square to the left of x
# There are 2k+1 total cum_areas so this find_le is O(log(k))
i, (area_left, x_left, y_interval) = find_le(cum_areas, (area,))
height, _, _ = heights[i]
x = x_left + (area - area_left) / height
# Step 3c. Choose a y value uniformly in the open spaces of that x value
# Interval.random_point() is O(log(k))
y = y_interval.random_point()
# Step 3d. Rejection sample until points are inside the circle
# On average, takes constant 4/pi random points to land within the circle
if x*x + y*y < 1:
points.append((x, y))
# Done!
return points
def display_points(rects, points):
fig, ax = plt.subplots(1)
ax.plot(*zip(*points), 'ko', markersize=1, alpha=0.5)
ax.axis([-1, 1, -1, 1])
ax.add_collection(
matplotlib.collections.PatchCollection(
[rect.as_matplotlib_rect() for rect in rects],
facecolor = 'r', edgecolor='None', alpha=0.5
)
)
plt.show()
def display_random_points(rects, num_points):
points = random_pick_points(rects, num_points)
display_points(rects, points)
def random_rects(k):
return [Rect(*[random.random()*2-1 for _ in range(4)]) for i in range(k)]
k = 20
# rects = random_rects(k)
# print("Rects:", rects)
n = 10000
# print("Points:", random_pick_points(rects, n))
display_random_points(random_rects(k), n)
Try it online I guess, if you can read a list of points.
k,nand rectangles' placement or whatever variable is significant for your approach \$\endgroup\$nandk? It seems to me that whether you have many points or many rectangles affects how effective it is to try to preprocess the rectangles and their area before generating points. \$\endgroup\$