from random import randint
from collections import Counter
import numpy as np
import matplotlib.pyplot as plt

NUM_EQUIVALENT_VOLUMES = 1000
MAX_CIVILIZATIONS = 5000
TRIALS = 1000
CIVILIZATION_STEP_SIZE = 100

x = []
y = []

for num_of_civs in range(2, MAX_CIVILIZATIONS + 2, CIVILIZATION_STEP_SIZE):
    civs_per_volume = num_of_civs / NUM_EQUIVALENT_VOLUMES
    num_single_civs = 0
    for trial in range(TRIALS):
        locations = []
        while len(locations) < num_of_civs:
            location = randint(1, NUM_EQUIVALENT_VOLUMES)
            locations.append(location)
        overlap_count = Counter(locations)
        overlap_rollup = Counter(overlap_count.values())
        num_single_civs += overlap_rollup[1]

    probability = 1 - (num_single_civs / (num_of_civs * TRIALS))

    print("{:.4f} {:.4f}".format(civs_per_volume, probability))
    x.append(civs_per_volume)
    y.append(probability)


coefficients = np.polyfit(x, y, 4)
p = np.poly1d(coefficients)
print("\n{}".format(p))
xp = np.linspace(0, 5)
_ = plt.plot(x, y, '.', xp, p(xp), '-')
plt.ylim(-0.5, 1.5)
plt.show()