import numpy as np

# parameters of Gaussian distribution
mu = 1
sigma = 1
# dimensions
d = 100
k = 10
# Gaussian matrix
G = np.random.normal(mu, sigma, (k, d))

# isotropic vector on the unit sphere
x = np.random.normal(mu, sigma, (d, 1))
# L2 Norm of X
l2NormX = np.sqrt(np.dot(np.transpose(x), x))
# re-sample x if norm too small to avoid floating point errors due to finite precision
while l2NormX < 0.001:
    x = np.random.normal(mu, sigma, d)
    l2NormX = np.sqrt(np.dot(np.transpose(x), x))

# normalize x
x /= l2NormX
# mapping function phi
phi_x = 1 / np.sqrt(k) * np.matmul(G, x)  # matrix multiplication
groundTruthValue = np.dot(np.transpose(x), x)[0][0]
estimate = np.sqrt(np.dot(np.transpose(phi_x), phi_x))[0][0]
print(groundTruthValue)
print(estimate)