Pythonic way to generate random uniformly distributed points within hollow square lamina

I've already posted a shorter and somehow unclear answer, here I've taken the time to produce what I think is a better answer.

Generalizing the OP problem, we have a "surface" composed of nsc * nsr squares disposed in nsc columns and nsr rows, and a "hole" composed of nhc * nhr squares (corresponding to the surface squares) disposed in a rectangular arrangement with nhr rows and nhc columns, with the origin placed at ohc, ohr

 +------+------+------+------+------+------+------+------+------+ nsr | | | | | | | | | | | | | | | | | | | | ...+------+------+------+------+------+------+------+------+------+ nsr-1 | | | oooo | oooo | | | | | | | | | oooo | oooo | | | | | | +------+------+------+------+------+------+------+------+------+ ... | | | oooo | oooo | | | | | | | | | oooo | oooo | | | | | | +------+------+------+------+------+------+------+------+------+ ... | | | oooo | oooo | | | | | | | | | oooo | oooo | | | | | | ohc+------+------+------+------+------+------+------+------+------+ 1 | | | | | | | | | | | | | | | | | | | | +------+------+------+------+------+------+------+------+------+ 0 0 1 2 ... ... nsc-1 nsc | | ohr=2 ohr+nhr 

Our aim is to draw n random points from the surface without the hole (the admissible surface) with a uniform distribution over the admissible surface.

We observe that the admissible area is composed of nsq = nsc*nsr -nhc*nhr equal squares: if we place a random point in an abstract unit square and then assign, with equal probability, that point to one of the squares of the admissible area, we have done our job.

In pseudocode, if random() samples from an uniformly distributed random variable over [0, 1)

(x, y) = (random(), random()) square = integer(random()*nsq) (x, y) = (x, y) + (column_of_square_origin(square), row_of_square_origin(square)) 

To speed up the process, we are using numpy and we try to avoid, as far as possible, explicit for loops.

With the names previously used, we need a listing of the origins of the admissible squares, to implement the last line of the pseudo code.

def origins_of_OK_squares(nsc, nsr, ohc, ohr, nhc, nhr): # a set of tuples, each one the origin of a square, for ALL the squares s_all = {(x, y) for x in range(nsc) for y in range(nsr)} # a set of tuples with the origin of the hole squares s_hole = {(x, y) for x in range(ohc,ohc+nhc) for y in range(ohr,ohr+nhr)} # the set of the origins of admissible squares is the difference s_adm = s_all - s_hole # return an array with all the origins --- the order is not important! # np.array doesn't like sets return np.array(list(s_adm)) 

We need to generate n random points in the unit square, organized in an array of shape (n,2)

 rand_points = np.random.random((n, 2)) 

We need an array of n admissible squares

 placements = np.random.randint(0, nsq, n) 

We translate each point in rand_points in one of the admissible squares, as specified by the elements of placements.

 rand_points += origins_of_OK_squares(nsc, nsr, ohc, ohr, nhc, nhr)[placements] 

taking advantage of the extended addressing that it is possible for numpy arrays, and its done...

In an overly compact function

import numpy as np def samples_wo_hole(n, nsc, nsr, ohc, ohr, nhc, nhr): s_all = {(x, y) for x in range(nsc) for y in range(nsr)} s_hole = {(x, y) for x in range(ohc,ohc+nhc) for y in range(ohr,ohr+nhr)} rand_points = np.random.random((n, 2)) placements = np.random.randint(0, nsc*nsr - nhc*nhr, n) return rand_points+np.array(list(s_all-s_hole))[placements] 

You have 8 equal unit squares in which it is admissible to place points, so draw as many points in an unit square as you want as in

x = np.random(n, 2) 

now it suffices to choose at random in which of the 8 admissible squares each point is placed

sq = np.random.randomint(0, 8, n) 

you need also an array of origins

delta = np.array([[0, 0], [1, 0], [2, 0], [0, 1], # no central square [2, 1], [0, 2], [1, 2] [2, 2]]) 

and finally

x = x + delta[sq] 

To generalize the solution, write a function to compute an array of the origins of the admissible squares, a possible implementation using sets being

def origins(n, hole_xor, hole_wd, hole_yor, hole_hg): all_origins = {(x,y) for x in range(n) for y in range(n)} hole_origins = {(x,y) for x in range(hole_xor, hole_xor+hole_wd) for y in range(hole_yor, hole_yor+hole_hg)} return np.array(list(all_origins-hole_origins))) 

and use it like this

delta = origins(12, 4,5, 6,2) n_squares = len(delta) # or n*n - width*height target_square = np.random.randomint(0, n_squares, size) x = x + delta[target_square]