Can I make my percolation in python function more efficient?

Question

We've been tasked with finding out how many randomly distributed discs with radius = 1 it takes to form a chain of discs from the left side to the other of a square with side length = sqrt(n), and then how many discs are required to cover the whole square.

We have code that does the job. The problem is that it takes way too long for large values of n. This is due to the function checking conditions for each disc before adding another disc over and over. When we get to, say 5000 discs, this becomes incredibly slow since it checks over 5000 discs every time it adds a new one.

Are there ways of making this code more efficient?

Here is the code:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle

#generates randomly spaced disks
def generate_disc(size):
    return np.random.uniform(0, size), np.random.uniform(0, size), 1.0  # Generate a disk with radius 1

#finds out of two disks touch each other
def touch_disc(disc1, disc2):
    x1, y1, r1 = disc1
    x2, y2, r2 = disc2
    distance = np.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)
    if distance <= (r1 + r2):
        return True
    return False

#find clusters of disks and returns the chains that exist
def find_chains(discs):
    clusters = []
    copy = list(discs)
    while copy != []:
        current_disc = copy[0]
        connected_discs = [current_disc]
        copy.remove(current_disc)
        i = 0
        while i < len(connected_discs):
            current_disc = connected_discs[i]
            for disc in copy[:]:
                if touch_disc(current_disc, disc):
                    connected_discs.append(disc)
                    copy.remove(disc)
            i += 1
        clusters.append(connected_discs)
    
    chains = []
    for cluster in clusters:
        if len(cluster) > 1:
            chains.append(cluster)
                
    return chains

#find out if one chain touches both the left and the right side
def complete_chain(chains):
    x = []
    for i, chain in enumerate(chains):
        x_chain = []
        if x != []:
                break
        for k in range(len(chain)):
            current_x = chain[k][0]
            x_chain.append(current_x)
        
        if (min(x_chain)-1) <= 0 and (max(x_chain)+1)>=size:
            x = x_chain
            return True
        else:
            return False

def plot_discs(discs, size):
    plt.figure(figsize=(8, 8))  # Set a square figure size
    plt.xlim(0, size)  # Set x-axis limit to [0, size]
    plt.ylim(0, size)  # Set y-axis limit to [0, size]
    plt.gca().set_aspect('equal')  # Set aspect ratio to 'equal' for square axes
    
    for disc in discs:
        x, y, r = disc
        # Ensure that the disc is plotted within the range [0, size]
        x = max(0, min(x, size))
        y = max(0, min(y, size))
        circle = plt.Circle((x, y), r, edgecolor='blue', fill=True, alpha = 0.5)
        plt.gca().add_patch(circle)
    
    plt.title('Randomly Plotted Discs')
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.show()

def place_discs_until_covered(size, max_discs):
    """Place discs with a radius of 1 until the entire square is covered.
    size is the size of the square. max_discs is the maximum number of discs we are willing to create."""
    
    # Initialize list to hold discs
    discs = []
    discs2 = []
    # Convert size to an integer
    size_int = int(size)
    tol = 0
    # Loop until the entire square is covered or maximum number of discs is reached
    for _ in range(max_discs):
        # Generate a new disc
        new_disc = generate_disc(size_int)
        # Add the new disc to the list
        discs.append(new_disc)
        discs2.append(new_disc)
        
        # Check if all points in the square are covered
        all_covered = True
        for x in range(size_int + 1):
            for y in range(size_int + 1):
                if tol <1:
                    chains = find_chains(discs2)
                    if complete_chain(chains):
                        first_chain = len(discs)
                        
                        print('First chain completed after',first_chain, 'discs.')
                        tol += 1
                        plot_discs(discs, size_int)
                # Check distance to each disc
                if not any(np.sqrt((x - disc[0]) ** 2 + (y - disc[1]) ** 2) <= 1 for disc in discs):
                    all_covered = False
                    break
            if not all_covered:
                break
        
        # If all points in the square are covered, break the loop
        if all_covered:
            break
    
    #Plot all the discs
    plot_discs(discs, size_int)
    
    # Calculate the number of discs needed
    num_discs_needed = len(discs)
    
    if num_discs_needed < max_discs:
        print(f"{num_discs_needed} discs were needed to cover the entire square.")
    else:
        print(f"The maximum number of discs ({max_discs}) was reached, but the square couldn't be fully covered.")
    
    return num_discs_needed, first_chain

n = 1000
size = np.sqrt(n)
place_discs_until_covered(size,10**3)

Any advice is appreciated.

The current state is a complete lack of ideas I'm afraid.

Answer 1

Analytical problem

This is a an analytical problem and it may be best to solve it first before coding to have an idea how best to structure your code.

If I understand you correctly:

• each disc has radius one,
• the square is of dimensions: $$ \sqrt{n} \times \sqrt{n} $$
• discs are placed sequentially at random and must cover the entire square.

If you assume, for simplicity, that the treatment of discs placed near the perimeter is to repeat/mirror rather than to truncate, then the probability that any point z=(x,y) in the square is covered by a single disc is the area of the disc divided by the area of the square:

$$ P(z \in Disc) = 1 - P(z \notin Disc) = \frac{\pi}{n} $$

Adding discs to reach a total of m discs independently and randomly means that,

$$ P(z \in \cup Discs) = 1 - P(z \notin \cup Discs) = 1 - \left ( P(z \notin Disc_1) ... P(z \notin Disc_m) \right ) = 1 - (1-\frac{\pi}{n})^m $$

This is a stochastic problem. Technically you could keep adding discs and might never cover the square but that is unlikely. If you estimated to a 99.99% likelihood that the square was covered, i.e. all points z are in the union of the discs, you can express the expected number of required discs relative to the dimension, n:

$$ m = \frac{ln(0.0001)}{ln(1-\frac{\pi}{n})} $$

Structuring code

I will ignore the disc chaining becuase that is a second task with not dissimilar concerns as here.

What is the optimal solution to ensure that the square is covered given some discs?

a) Don't exhaustively search every point every time you add a sphere. Once a point is covered it will be covered indefinitely. You do not need to check that point again.

b) To achieve the same 99.99% probability of coverage you need to create 10,000 gridpoints, i.e. 100 x 100.

c) Once you have discovered just one point that is not covered you need to add a disc and keep adding until that exact single point is covered, Once it is move on to check the next point.

d) I would optimise the stepping to the next point by doing two things. Firsly I would start in the bottom left corner and traverse up and down moving right at the boundary. Secondly when I stepped to a new point I would run a fast check that the same disc that covered the preceding point did not also cover that new point, becuase it is likely that a disc close by will cover multiple points. Once you find an uncovered point keep adding discs until it is covered and then move on.

This will most likely take a fraction of the time than it is currently taking for your code to run.

Example

The following code does exactly what I described above and it will measure upto 50,000 discs in about 1 second. It also follows the statistical analysis described above, giving me confidence it is free from bugs. I hope this skeletal structrue helps.

import numpy as np

SQRT_N = 100.0  # length of side
RV_DISC_CENTRES = np.random.rand(80000, 2) * SQRT_N
GRID_SIZE = 0.01

def get_next_gridpoint(current, direction):
    "Traverse the grid up and down moving to a new point"
    if direction == "up":
        if current[1] >= SQRT_N:
            return (current[0] + GRID_SIZE, SQRT_N), "down"
        else:
            return (current[0], current[1] + GRID_SIZE), "up"
    else:
        if current[1] <= 0.0:
            return (current[0] + GRID_SIZE, 0.0), "up"
        else:
            return (current[0], current[1] - GRID_SIZE), "down"

def is_gridpoint_in_disc(gridpoint, discpoint):
    """Ensure the euclidean distance between points is less than disc radius"""
    return ((gridpoint[0] - discpoint[0])**2 + (gridpoint[1] - discpoint[1])**2) < 1.0

def search_until_covered(gridpoint):
    """Iterate through all discs until the gridpoint is covered and return the disc num"""
    for i in range(RV_DISC_CENTRES.shape[0]):
        if is_gridpoint_in_disc(gridpoint, RV_DISC_CENTRES[i, :]):
            return True, i
    return False, None

curr_gridpoint, direction = (0.0, 0.0), "up"
disc_num, num_discs = None, 0.0
while curr_gridpoint[0] <= SQRT_N and curr_gridpoint[1] <= SQRT_N:
    if disc_num and is_gridpoint_in_disc(curr_gridpoint, RV_DISC_CENTRES[disc_num, :]):
        curr_gridpoint, direction = get_next_gridpoint(curr_gridpoint, direction)
    else:
        is_covered, disc_num = search_until_covered(curr_gridpoint)
        if not is_covered:
            raise ValueError(f"Gridpoint {curr_gridpoint} cannot be covered.")
        if disc_num > num_discs:
            num_discs = disc_num
        curr_gridpoint, direction = get_next_gridpoint(curr_gridpoint, direction)

print(f"Number of discs required for coverage: {num_discs}")

Profiling

Notice that the exhaustive search is only called 99 times for the 10000 gridpoints:

Answer 2

One idea is to split the square into a grid[][] of smaller 2x2 squares. Then, for each of these squares store all circles whose center is in that square. This means that when you insert a new circle you only need to check for intersections with any circle in at most 9 adjacent squares (in expectation far fewer than checking all of the circles).

As for the connectivity problem, you can use the union-find data structure. Each time you find two intersecting circles you unite their components. For each component you keep track whether it touches the left or the right edge of the square. Once you unite two components where one touches the left and the other touches the right, you are done. (This data structure is quite efficient and relatively simple to implement)