sph_ws_natural_hypothesis.py

import numpy as np
from scipy.spatial.distance import pdist, squareform
from itertools import combinations

sph = 2000

def monte_plot(sph):
    X = 100 * np.random.random(sph * 2).reshape(sph, 2)
    radial_dists = np.sqrt(np.sum((X - np.array([50,50]))**2, axis=1))
    plot_X = X[radial_dists <= 3.999]
    best = 0
    for count in range(plot_X.shape[0]):
        for comb in combinations(range(plot_X.shape[0]), r=count):
            pairwise_dists = pdist(plot_X[comb,:])
            if pairwise_dists.size > 0 and np.min(pairwise_dists) >= 0.5:
                if count > best:
                    best = count
    return best

import matplotlib.pyplot as plt

results = []
for i in range(1000):
    results.append(monte_plot(2000))
    print(i)
plt.hist(results)
plt.show()