measure_distance.py

import numpy as np
from matplotlib import pyplot as plt
import itertools

np.set_printoptions(precision=4, suppress=True)


def all_waves(N):
    for i in itertools.count():
        k = i + 1
        for func in [np.cos, np.sin]:
            quant = func(2 * np.pi * k/N * (np.arange(N)))
            quant = quant / np.linalg.norm(quant)
            yield quant

def all_basis(N):
    comps = list(itertools.islice(all_waves(N), N-1))
    for w1, w2 in itertools.product(comps, comps):
        yield np.reshape(w1, (N, 1)) @ np.reshape(w2, (1, N))


from random import gauss
def make_rand_direction(dims):
    vec = [gauss(0, 1) for i in range(dims)]
    mag = sum(x**2 for x in vec) ** .5
    return np.array([x/mag for x in vec])

        
bases = list(all_basis(60))
uniform = np.ones((60,60)) / 60**2

def distance_trial():
    accum = np.zeros((60, 60))
    d = make_rand_direction(59**2)
    for idx, i in enumerate(bases):
        accum += d[idx] * i

    to_go_over_1 = (1 - uniform) / accum
    to_go_under_0 = (0 - uniform) / accum

    upper_bound_1 = np.amin(to_go_over_1[to_go_over_1 > 0])
    lower_bound_1 = np.amax(to_go_over_1[to_go_over_1 < 0])

    upper_bound_0 = np.amin(to_go_under_0[to_go_under_0 > 0])
    lower_bound_0 = np.amax(to_go_under_0[to_go_under_0 < 0])

    upper_bound = min(upper_bound_1, upper_bound_0)
    lower_bound = max(lower_bound_1, lower_bound_0)
    
    return upper_bound - lower_bound

print(np.mean([distance_trial() for x in range(1000)]))