Leon's shooting range problem 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.



Task

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)
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 &quot;open area&quot; 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 &lt;= self_range[0] &lt;= hi &lt; self_range[1]:
                new_ranges.append((hi, self_range[1]))
            elif self_range[0] &lt; lo &lt;= self_range[1] &lt;= hi:
                new_ranges.append((self_range[0], lo))
            elif self_range[0] &lt; lo &lt; hi &lt; self_range[1]:
                new_ranges.append((self_range[0], lo))
                new_ranges.append((hi, self_range[1]))
            elif lo &lt;= self_range[0] &lt; self_range[1] &lt;= hi:
                pass
            elif hi &lt;= self_range[0] or self_range[1] &lt;= 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 &quot;I&quot; + 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&quot;Rect({self.left}, {self.right}, {self.bottom}, {self.top})&quot;

    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 &lt; ends[-1].top:
                active_rect = starts.pop()
                if not active_rects and active_rect.bottom &gt; 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 &lt; 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) &lt; 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 &lt; 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(&quot;Rects:&quot;, rects)
n = 10000
# print(&quot;Points:&quot;, 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.