import numpy as np
from matplotlib import pyplot as plt
import numpy.random as ra

%matplotlib inline

"""
Explanation: Inverse Transform Sampling
This notebook will conceptualize how inverse transform sampling works
End of explanation
"""


spectrum = [[1, 2, 3, 4, 5, 6],[2000, 4040, 6500, 6000, 4020, 2070]]
energies = np.array(spectrum[0])
fluxes = np.array(spectrum[1])
spectrum

"""
Explanation: Below is a spectrum which follows an almost bell-curve type distribution (anyway, the specific type of distribution is not important here).
End of explanation
"""


prob = fluxes/float(sum(fluxes))
cum_prob = np.cumsum(prob)
cum_prob

"""
Explanation: Below, first we compute probabilities of flux. Afterwards, we compute the cumulative probability.
End of explanation
"""


N = 10000
R = ra.uniform(0, 1, N)
R[1:10]

"""
Explanation: We draw ten thousand numbers from uniform random distribution.
End of explanation
"""


gen_energies = [int(energies[np.argwhere(cum_prob == min(cum_prob[(cum_prob - r) > 0]))]) for r in R]
gen_energies[1:10]

"""
Explanation: We assign energies to events corresponding to the random number drawn.
Note: The command below finds bin interval using a single command. I am not sure though that it's very readble. Would
we want to split that in multiple lines and maybe use explicit loops to make it more readable? Or is it fine as it is?
Comments?
End of explanation
"""


gen_energies = ((np.array(gen_energies) - 1) / 1).astype(int)
times = np.arange(1, 6, 1)
lc = np.bincount(gen_energies, minlength=len(times))
lc

plot1, = plt.plot(lc/float(sum(lc)), 'r--', label='Assigned energies')
plot2, = plt.plot(prob,'g',label='Original Spectrum')
plt.xlabel('Energies')
plt.ylabel('Probability')
plt.legend(handles=[plot1,plot2])
plt.show()

"""
Explanation: Histogram energies to get shape approximation.
End of explanation
"""


