# -*- coding: utf-8 -*-
"""Q3.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1nuVOMaV1Wr9RXpnbiyiNker8qGPQ86Nl
"""

import os 
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# get 2-nd norms
def getNorms(arr):
  norm = []
  for i in range(len(arr)):
    for j in range(i + 1, len(arr)):
      temp = np.linalg.norm(arr[i] - arr[j])
      norm.append(temp)
  return norm

# create and initialize covariance matrix
def buildCovMatrix(d):
  cov = [[0] * d for i in range(d)]
  for i in range(0, d):
    cov[i][i] = 1
  return cov

"""d = 3 and sample_size = 10000"""

d = 3
mean = [0] * d
cov = buildCovMatrix(d)

# generate 10000 random samples, 10000 is too much, here use 1000 as an example
sample_size = 1000
mul_normal = np.random.multivariate_normal(mean, cov, sample_size)

# calculate the 2-nd norm of two data
norm = getNorms(mul_normal)
print(f'when d = 3, the norms are {norm}')

# calculate average and standard deviation of norm
norm = np.array(norm)
average = np.mean(norm)
standard = np.std(norm)
print(f'when d = 3, the average of norms is {average}')
print(f'when d = 3, the standard deviation of norms is {standard}')

# plot the histogram of Euclidean norms
plt.hist(norm, bins=100, density=0, facecolor="blue")

"""d = 3, 50, 100, 200, 500, 1000"""

d = [3, 50, 100, 200, 500, 1000]

# d = [3, 50]
average_ans = []
standard_ans = []

for x in range(0, len(d)):
  mean = [0] * d[x]
  cov = buildCovMatrix(d[x])

  # get norms of each d
  sample_size = 1000
  arr = np.random.multivariate_normal(mean, cov, sample_size)
  norm = np.array(getNorms(arr))
  print(norm)
  
  # calculate mean and standard deviation
  average_ans.append(np.mean(norm))
  standard_ans.append(np.std(norm))

print(average_ans)
print(standard_ans)

ax1 = plt.subplot(2,1,1)
ax2 = plt.subplot(2,1,2)

# plot the average and standard deviation
plt.sca(ax1)
plt.plot(d, average_ans)
plt.xlabel('d', fontsize=10)
plt.ylabel('average', fontsize=10)

plt.sca(ax2)
plt.plot(d, standard_ans)
plt.xlabel('d', fontsize=10)
plt.ylabel('standard deviation', fontsize=10)

"""what I find is that the average is going up but the standard deviation tends to be convergence"""