# Path: modality/util/ApproxGaussianKDE.py
# class ApproxGaussianKDE(object):
#     """
#         Approximate Gaussian kernel density estimate, using binned data. 
#         For each bin, the data in the bin is replaced by its mean. 
#         The bin width is chosen so that the maximal absolute error is 
#         bounded by tol/(h*sqrt(2*pi). This is shown in the accompanying
#         documentation "Approxmation of kernel density estimator".
# 
#         Apart from binning, the implementation mimics that of 
#         KernelDensity from sklearn.neighbors.
#     """
# 
#     def __init__(self, data, bandwidth, tol=1e-4):
#         self.h = bandwidth
#         self.kernel = lambda u: np.exp(-u**2/2.0)
#         self.data_orig = data
#         data = np.sort(data)
#         bin_width = np.sqrt(2*tol)*self.h
#         bins = np.arange(data[0]-bin_width/2.0, data[-1]+3*bin_width/2.0, bin_width)
#         ind_bins = np.searchsorted(data, bins)
#         bin_sizes = np.diff(ind_bins)
#         empty = bin_sizes == 0
#         self.weights = bin_sizes[~empty]
#         #print "len(self.weights) = {}".format(len(self.weights))
#         self.datah = np.array([np.mean(ary) for (ary, empty_) in zip(np.split(data, ind_bins[1:-1]), empty) if not empty_])
#         self.datah /= self.h
#         self._norm_factor = np.sqrt(2*np.pi)*self.h*sum(self.weights)
# 
#     def evaluate_prop(self, x):  # returns values proportional to kde.
#         xh = np.array(x).reshape(-1)/self.h
#         res = np.zeros(len(xh))
#         if len(xh) > len(self.datah):  # loop over data
#             for data_, weight in zip(self.datah, self.weights):
#                 res += weight*self.kernel(data_-xh)
#         else:  # loop over x
#             for i, x_ in enumerate(xh):
#                 res[i] = np.sum(self.weights*self.kernel(self.datah-x_))
#         return res
# 
#     def evaluate(self, x):
#         res = self.evaluate_prop(x)
#         return res / self._norm_factor
# 
#     def score_samples(self, x):
#         return np.log(self.evaluate(x))
# 
#     def distr(self, x):
#         return np.mean(norm.cdf(x - self.data_orig, scale=self.h))


import matplotlib.pyplot as plt
import numpy as np
        import time
from scipy.signal import argrelextrema
from .util import ApproxGaussianKDE as KDE
    from sklearn.neighbors import KernelDensity
        from .bandwidth_fm_test import find_reference_distr


def merge_into(z_new, z):
    if len(z) == 0:
        return z_new
    z_merged = np.zeros((2*len(z)-1,))
    z_merged[np.arange(0, len(z_merged), 2)] = z
    z_merged[np.arange(1, len(z_merged), 2)] = z_new
    return z_merged

if __name__ == '__main__':

    if 0:
        N = 1000
        data = np.hstack([np.random.randn(N/2), np.random.randn(N/4)+4])
        h = 0.1
        print("is_unimodal_kde(h, data) = {}".format(is_unimodal_kde(h, data)))
        #plt.show()
        t0 = time.time()
        h_crit = critical_bandwidth(data)
        print("critical_bandwidth(data) = {}".format(h_crit))
        t1 = time.time()
        print("pval_silverman(data) = {}".format(pval_silverman(data)))
        t2 = time.time()
        print("Critical bandwidth computation time: {}".format(t1-t0))
        print("Silverman test computation time: {}".format(t2-t1))

        fig, ax = plt.subplots()
        ax.hist(data, bins=50, normed=True)
        x_grid = np.linspace(np.min(data)-2, np.max(data)+2, 100)

# Path: modality/util/ApproxGaussianKDE.py
# class ApproxGaussianKDE(object):
#     """
#         Approximate Gaussian kernel density estimate, using binned data. 
#         For each bin, the data in the bin is replaced by its mean. 
#         The bin width is chosen so that the maximal absolute error is 
#         bounded by tol/(h*sqrt(2*pi). This is shown in the accompanying
#         documentation "Approxmation of kernel density estimator".
# 
#         Apart from binning, the implementation mimics that of 
#         KernelDensity from sklearn.neighbors.
#     """
# 
#     def __init__(self, data, bandwidth, tol=1e-4):
#         self.h = bandwidth
#         self.kernel = lambda u: np.exp(-u**2/2.0)
#         self.data_orig = data
#         data = np.sort(data)
#         bin_width = np.sqrt(2*tol)*self.h
#         bins = np.arange(data[0]-bin_width/2.0, data[-1]+3*bin_width/2.0, bin_width)
#         ind_bins = np.searchsorted(data, bins)
#         bin_sizes = np.diff(ind_bins)
#         empty = bin_sizes == 0
#         self.weights = bin_sizes[~empty]
#         #print "len(self.weights) = {}".format(len(self.weights))
#         self.datah = np.array([np.mean(ary) for (ary, empty_) in zip(np.split(data, ind_bins[1:-1]), empty) if not empty_])
#         self.datah /= self.h
#         self._norm_factor = np.sqrt(2*np.pi)*self.h*sum(self.weights)
# 
#     def evaluate_prop(self, x):  # returns values proportional to kde.
#         xh = np.array(x).reshape(-1)/self.h
#         res = np.zeros(len(xh))
#         if len(xh) > len(self.datah):  # loop over data
#             for data_, weight in zip(self.datah, self.weights):
#                 res += weight*self.kernel(data_-xh)
#         else:  # loop over x
#             for i, x_ in enumerate(xh):
#                 res[i] = np.sum(self.weights*self.kernel(self.datah-x_))
#         return res
# 
#     def evaluate(self, x):
#         res = self.evaluate_prop(x)
#         return res / self._norm_factor
# 
#     def score_samples(self, x):
#         return np.log(self.evaluate(x))
# 
#     def distr(self, x):
#         return np.mean(norm.cdf(x - self.data_orig, scale=self.h))


# Path: modality/critical_bandwidth.py
import matplotlib.pyplot as plt
import numpy as np
        import time
from scipy.signal import argrelextrema
from .util import ApproxGaussianKDE as KDE
    from sklearn.neighbors import KernelDensity
        from .bandwidth_fm_test import find_reference_distr



def merge_into(z_new, z):
    if len(z) == 0:
        return z_new
    z_merged = np.zeros((2*len(z)-1,))
    z_merged[np.arange(0, len(z_merged), 2)] = z
    z_merged[np.arange(1, len(z_merged), 2)] = z_new
    return z_merged

if __name__ == '__main__':

    if 0:
        N = 1000
        data = np.hstack([np.random.randn(N/2), np.random.randn(N/4)+4])
        h = 0.1
        print("is_unimodal_kde(h, data) = {}".format(is_unimodal_kde(h, data)))
        #plt.show()
        t0 = time.time()
        h_crit = critical_bandwidth(data)
        print("critical_bandwidth(data) = {}".format(h_crit))
        t1 = time.time()
        print("pval_silverman(data) = {}".format(pval_silverman(data)))
        t2 = time.time()
        print("Critical bandwidth computation time: {}".format(t1-t0))
        print("Silverman test computation time: {}".format(t2-t1))

        fig, ax = plt.subplots()
        ax.hist(data, bins=50, normed=True)
        x_grid = np.linspace(np.min(data)-2, np.max(data)+2, 100)

        ax.plot(x_grid, KDE(data, h_crit).evaluate(x_grid), linewidth=2, color='black')