# Leon's story

Leon is a professional sling shooter and he comes to a shooting range everyday to practice. A casual target is not a challenge for him anymore so before shooting he first covers the target of radius 1.0 with k rectangle stickers. He then fires n shots that cannot hit the stickers.

What's special about Leon is that he always shoots randomly but with a uniform distribution.

Leon everyday takes on a bigger challenge and puts more stickers. However he started noticing that shooting takes too much time now so he decided to ask you for the best strategy.

Implement a program which will randomly pick n points with the lowest time complexity possible. Each point should be inside a circle target, and none of them lands in any of k rectangle stickers.

Input:

List of k rectangles in whatever form fits you.

Eg. where each rectangle is a list of:

[x_min, x_max, y_min, y_max]


The list could be:

[[0.1, 0.3, -0.2, 0.1], [-0.73, -0.12, 0.8, 1.4]]


Output:

List of coordinates of n points in whatever form fits you.

Rules:

• Points have to be picked randomly with uniform distribution
• Every point lays inside a circle
• None of the points lays inside a rectangle sticker
• Stickers can overlap each other
• Circle's radius is equal to 1.0 and its center is at coordinates (0.0, 0.0)
• Every coordinate is in form of floating point number
• For edge cases - choose whatever fits you

EDIT: Thanks to @Luis Mendo and @xnor it turned out that one more rule might be useful.

• With each new sticker 1/4 of the remaining board will be covered. This means that having k stickers the uncovered part of the target will be (3/4)^k of the whole.

## Rating

The best time complexity wins!

EDIT:

For the sake of rating we will assume that k and n are linearly proportional, meaning that O(k^2 + n) will be rated equally to O(k + n^2)

• Do you mean asymptotic time complexity? Also, how do you define time complexity for algorithms with random running time (like rejection samping)? May 12, 2020 at 15:06
• @LuisMendo Yes, asymptotic time complexity. For random running time take average number of operations for different k, n and rectangles' placement or whatever variable is significant for your approach May 12, 2020 at 15:24
• What precision is required on the chosen floats? May 12, 2020 at 16:50
• One possible strategy is rejection sampling, that is picking a random point in the circle, and if it overlaps a rectangle, then repeating. How should one treat the expected number of trials for the purposes of time complexity. For a "typical case", this is a constant, but if the rectangles are made to cover all but a tiny part of the circle, it can be arbitrarily high.
– xnor
May 13, 2020 at 5:59
• How do we compare run-times of solutions that make different tradeoffs in n and k? 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.
– xnor
May 13, 2020 at 6:05

# 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.

1. 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

2. 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

3. 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
# 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
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)

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
# 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])
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.