Connecting to random points in a 2d numpy array based on distance

I have a 2d numpy array and select a random coordinate position (say 10×10 array and start at position 2,3). I want to randomly connect to 40% of the other points in the 2d array effectively generating a list of tuples [(x1, y1), (x2, y2) …] where the list is 40% of the other coordinates.
An additional constraint, however, is the goal is to reduce connection probability the farther the points are away from one another (so point 2,3 is far more likely to connect to 2,2 than 9, 8 but yet still be random so there is a chance albeit small of connecting to 9, 8).

I believe I need to create some sort of Guassian function centered on 2,3 and use this to select the points, but any Gaussian I create will generate non-integer values – requiring additional logic as well as presents problem of needing to handle x and y dimensions separately.

Currently, I am trying to use np.meshgrid with gauss = np.exp(-(dst2 / (2.0 * sigma2)))

Is there an easier way to do this or a different approach someone might recommend?

Answer

This problem is well suited for rejection sampling. Basically you randomly choose a point, and select if a connection should be made based on a defined probability. You have to take in account that there are many more points at further distance than at # closer distances (its number grows with radius), so maybe you have to introduce an extra weighing in the probability function. In this case I choose to use an exponential decay probability.

This code is not optimal in terms of speed, particularly for higher connectivity percents, but the ideas are better presented this way: see below for a better option.

import numpy as np
from numpy.random import default_rng

rng = default_rng()
board = np.zeros((100, 100), dtype=bool)
percent_connected = 4
N_points = round((board.size - 1) * percent_connected/100)
center = np.array((20, 30))
board[tuple(center)] = True # remove the center point from the pool
dist_char = 35  # characteristic distance where probability decays to 1/e

endpoints = []
while N_points:
    point = rng.integers(board.shape)
    if not board[tuple(point)]:
        dist = np.sqrt(np.sum((center-point)**2))
        P = np.exp(-dist / dist_char)
        if rng.random() < P:
            board[tuple(point)] = True
            endpoints.append(point)
            N_points -= 1
board[tuple(center)] = False # clear the center point

# Graphical test
import matplotlib.pyplot as plt

plt.figure()
for ep in endpoints:
    plt.plot(*zip(center, ep), c="blue")

A slightly faster approach is much faster at higher connectivity:

rng = default_rng()
board = np.zeros((100, 100), dtype=bool)
percent_connected = 4
N_points = round((board.size - 1) * percent_connected/100)
center = np.array((20, 30))
board[tuple(center)] = True # remove the center point from the pool
dist_char = 35  # characteristic distance where probability decays to 1/e
flat_board = board.ravel()
endpoints = []
while N_points:
    idx = rng.integers(flat_board.size)
    while flat_board[idx]:
        idx += 1
        if idx >= flat_board.size:
            idx = 0
    if not flat_board[idx]:
        point = np.array((idx // board.shape[0], idx % board.shape[0]))
        dist = np.sqrt(np.sum((center-point)**2))
        P = np.exp(-dist / dist_char)
        if rng.random() < P:
            flat_board[idx] = True
            endpoints.append(point)
            N_points -= 1
board[tuple(center)] = False # clear the center point


plt.figure()
for ep in endpoints:
    plt.plot(*zip(center, ep), c="blue")