import numpy as np
from scipy import stats

import matplotlib.pyplot as plt
%matplotlib inline

import evapy

"""
Explanation: Basic statistical analysis of time series
In this example, basic time series statistical analysis is demonstrated.
End of explanation
"""


t = np.arange(0., 3*3600., 0.1)

"""
Explanation: Generating a Gaussian stochastic signal
Lets start by generating a Gaussian time series signal. The duration is set to 3 hours and sampling rate is 10Hz.
End of explanation
"""


start_period = 15.
stop_period = 5.

discrete_size = 1000

amp_norm = stats.norm.rvs(size=discrete_size)
phase = stats.uniform.rvs(size=discrete_size)*2.*np.pi

signal_mean = 0.

freq = np.linspace(2.*np.pi*1./start_period, 2.*np.pi*1./stop_period, num=discrete_size)

"""
Explanation: The signal will be generated using a rectangular power spectrum between 5 and 15 seconds of period. It is important to use sufficiently many discretizations to obtain a non-repeating approximation.
End of explanation
"""


signal = np.zeros(len(t)) + signal_mean

for i, freq_i in enumerate(freq):
    signal = signal + amp_norm[i]*np.cos(freq_i*t + phase[i])

plt.figure('Snapshot of signal', figsize=(16,9))
plt.title('Snapshot of the signal')
plt.plot(t[1500:2000], signal[1500:2000])
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.grid()
plt.show()

"""
Explanation: Generating the signal by super-positioning a large number of cosine functions with stochastic amplitudes and phases.
End of explanation
"""


descr_stats = stats.describe(signal)

print('mean: {}'.format(descr_stats.mean))
print('std: {}'.format(np.sqrt(descr_stats.variance)))
print('skewness: {}'.format(descr_stats.skewness))

"""
Explanation: Statistical analysis of the signal
Lets see if the signal resambles a Guassian signal with zero mean.
End of explanation
"""


params = stats.norm.fit(signal)
signal_dist = stats.norm(*params)

x = np.linspace(-100, 100, num=200)

plt.figure('Signal histogram', figsize=(16,9))
plt.title('Signal histogram and Normal PDF')
plt.plot(x, signal_dist.pdf(x), color='r')
plt.hist(signal, bins=50, density=True)
plt.show()

"""
Explanation: Estimate the normal distribution paramters from data and plot the normal pdf with the histogram for data.
End of explanation
"""


peak_index = evapy.evstats.argrelmax(signal)
peak = signal[peak_index]

peak_dc_index = evapy.evstats.argrelmax_decluster(signal, x_up=np.mean(signal))
peak_dc = signal[peak_dc_index]

plt.figure('Finding the peaks', figsize=(16,9))
plt.title('Finding the peaks')
plt.plot(t, signal, 'k-', label='signal')
plt.plot(t[peak_dc_index], peak_dc, 'r^', markersize=12, label='peak declustered')
plt.plot(t[peak_index], peak, 'bo', label='peak')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.xlim(100,200)
plt.grid()
plt.legend()
plt.show()

"""
Explanation: It is clear that the signal is indeed Gaussian!
Extreme value analysis
Finding the peaks
evapy package provides robust and fast methods to identify upcrossings and peaks. The declustered peaks refer to the largest peak per upcrossing. hus, for a narrow-banded signal, all peaks are declustered. In this example, declustering will occur for mean-upcrossing.
End of explanation
"""


uc_0 = evapy.evstats.argupcross(signal, x_up=np.mean(signal))

print('number of upcrossings above the mean: {}'.format(len(uc_0)))

"""
Explanation: Counting upcrossings
Lets find the number of upcrossings for different levels using the evapy package.
End of explanation
"""


t[-1]/len(uc_0)

"""
Explanation: Since the signal is 10800 seconds long, this mean that the Tz ((mean or) zero crossing period) of the signal is
End of explanation
"""


uc_1 = evapy.evstats.argupcross(signal, x_up=np.mean(signal)+np.std(signal))

print('number of upcrossings above the mean+std: {}'.format(len(uc_1)))

"""
Explanation: Now, lets find the upcrossing rate at mean+std level.
End of explanation
"""


len(uc_1)/t[-1]

"""
Explanation: The upcrossing rate is then:
End of explanation
"""


eak_index = evapy.evstats.argrelmax(signal)
peak = signal[peak_index]

peak_dc_index = evapy.evstats.argrelmax_decluster(signal, x_up=np.mean(signal))
peak_dc = signal[peak_dc_index]

r = np.linspace(-10, 100, num=100)

"""
Explanation: Statistical analysis of peaks
Now that we have established the peaks, lets investigate their distribution. Using different distribution functions the evapy package offers.
End of explanation
"""


params_ray = evapy.distributions.rayleigh.fit(peak_dc)
peak_dc_ray = evapy.distributions.rayleigh(*params_ray)

print('location: {}'.format(params_ray[0]))
print('scale: {}'.format(params_ray[1]))

plt.figure('Peaks histogram and Rayleigh distribution', figsize=(16,9))
plt.plot(r, peak_dc_ray.pdf(r), color='r')
plt.hist(peak_dc, bins=25, density=True)
plt.show()

"""
Explanation: Using Rayleigh distribution:
End of explanation
"""


params_wb = evapy.distributions.weibull.fit(peak_dc, scale=20.)
peak_dc_wb = evapy.distributions.weibull(*params_wb)

print('shape: {}'.format(params_wb[0]))
print('location: {}'.format(params_wb[1]))
print('scale: {}'.format(params_wb[2]))

"""
Explanation: Using the Weibull distribution:
Note: Sometimes, the estimation alogorithm fails if the starting values for scale parameter is too far off. This is a known issue and we are working on it. In the mean time, you may provice a reasonable starting value for the scale parameter.
End of explanation
"""


plt.figure('Peaks histogram and Weibull distribution', figsize=(16,9))
plt.plot(r, peak_dc_wb.pdf(r), color='r')
plt.hist(peak_dc, bins=25, density=True)
plt.show()

"""
Explanation: The shape parameter should be close to 2. Sometimes the estimation alogorithm fails and divereges  for 3 parameter Weibull using the MLE.
Lets see how the histogram looks.
End of explanation
"""


params_wb = evapy.distributions.weibull.fit(peak_dc, floc=np.mean(signal))
peak_dc_wb = evapy.distributions.weibull(*params_wb)

print('shape: {}'.format(params_wb[0]))
print('location: {}'.format(params_wb[1]))
print('scale: {}'.format(params_wb[2]))

plt.figure('Peaks histogram and Weibull distribution', figsize=(16,9))
plt.plot(r, peak_dc_wb.pdf(r), color='r')
plt.hist(peak_dc, bins=25, density=True)
plt.show()

"""
Explanation: Often, the to provide a good starting value for the scale parameter. Or to lock the location parameter based on prior knowledge. In our case, the signal mean can be a reasonable choice.
End of explanation
"""


peak_dc_trunc = peak_dc[peak_dc>=35.]

"""
Explanation: Using the generalized exponential tail distribution:
Note that the generalized exponential tail distribution is particularly suitable for truncated data. Lets truncate our peak data at 35.
End of explanation
"""


print('original data size: {}'.format(len(peak_dc)))
print('truncated data size: {}'.format(len(peak_dc_trunc)))

"""
Explanation: We have now discareded most of our data:
End of explanation
"""


params_get = evapy.distributions.genexptail.fit(peak_dc_trunc, scale=20., floc=np.mean(signal))
peak_dc_get = evapy.distributions.genexptail(*params_get)

print('shape: {}'.format(params_get[0]))
print('q: {}'.format(params_get[1]))
print('location: {}'.format(params_get[2]))
print('scale: {}'.format(params_get[3]))

plt.figure('Peaks histogram and generalized exponential tail distribution', figsize=(16,9))
plt.plot(r, peak_dc_get.pdf(r), color='r')
plt.hist(peak_dc_trunc, bins=15, density=True)
plt.show()

"""
Explanation: WARNING: At this point, the estimation alogrithm for the generalized exponential tail distribution is very unstable. It is recommended to lock the location parameter and provide a good starting value for the scale parameter. We are working on it!
End of explanation
"""


max_obs = np.max(signal)
print('largest maxima: {}'.format(max_obs))

"""
Explanation: Estimating the largest maxima
We can use the the relationship between the peak distributions and the maxima distributions to establish an estimate for the largest maxima. First, lets see whats the largest observed value in our signal.
End of explanation
"""


rloc, rscale = evapy.distributions.rayleigh.fit(peak_dc, floc=0.)
rN = len(peak_dc)

maxima_ray_dist = evapy.distributions.acer_o1(2., rN, loc=rloc, scale=np.sqrt(2.)*rscale)

median_maxima = maxima_ray_dist.ppf(0.5)
print('median maxima: {}'.format(median_maxima))

"""
Explanation: Using Rayleig:
End of explanation
"""


wshape, wloc, wscale = evapy.distributions.weibull.fit(peak_dc, floc=0.)
wN = len(peak_dc)

maxima_wb_dist = evapy.distributions.acer_o1(wshape, wN, loc=wloc, scale=wscale)

median_maxima = maxima_wb_dist.ppf(0.5)
print('median maxima: {}'.format(median_maxima))

"""
Explanation: Using Weibull:
End of explanation
"""


gshape, qn, gloc, gscale = evapy.distributions.genexptail.fit(peak_dc_trunc, floc=0., scale=20.)
gN = len(peak_dc_trunc)

maxima_get_dist = evapy.distributions.acer_o1(gshape, gN*qn, loc=gloc, scale=gscale)

median_maxima = maxima_get_dist.ppf(0.5)
print('median maxima: {}'.format(median_maxima))

"""
Explanation: Using generalized exponential tail:
End of explanation
"""


