import numpy as np
import scipy.stats as st
import matplotlib.pyplot as plt
from tools import cdf
from scipy import io

np.random.seed(42)

# creates the un-normalized density
f = lambda u: (np.sin(6*u)**2 + 3 * np.cos(u)**2 * np.sin(4*u)**2 + 1.) * np.exp(- u**2 / 2.)
x = np.linspace(-4, 4, 1001)

#plots
plt.plot(x, f(x)/st.norm.pdf(x), linewidth = 1.)
plt.show()

#defines some hyperparameters
N = 10000
X = np.zeros(N)
i = 0
j = 0
C = 5 * np.sqrt(2 * np.pi)

# Does the AR part
while i < N:
	x_prop = st.norm.rvs()
	U = st.uniform.rvs()
	if U <= f(x_prop) / (C* st.norm.pdf(x_prop)):
		X[i] = x_prop.copy()
		i = i + 1
	j = j + 1

#computes the acceptance rate
acc = float(i) / j
k = 1. / (acc * C)
print('Acceptance rate: ' + str(acc))
print('Normalization constant: ' + str(k))

#plots results
plt.plot(x, C*st.norm.pdf(x), linewidth = 1., label = r'$Cg(x)$')
plt.plot(x, f(x), linewidth = 1., label = r'$\tilde{f}(x)$')
plt.plot(x, k * f(x), linewidth = 1., label = r'$f(x)$')
plt.hist(X, bins = 100, density = True)
plt.legend()
plt.show()

truecdf = io.loadmat('truecdf.mat')['truecdf'][0,:]
truex = np.linspace(-4, 4, 1601)
xx, Nxx = cdf(X, [-4, 4])
plt.plot(xx, Nxx, linewidth = 1., label = 'Empirical CDF')
plt.plot(truex, truecdf, linewidth = 1., label = 'True CDF')
plt.legend()
plt.show()