# 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/util/GaussianMixture1d.py
# class GaussianMixture1d(object):
# 
#     def __init__(self, means, covs, weights):
#         self.mus = means
#         self.sigmas = np.sqrt(covs)
#         self.pis = weights
#         self.kernel = lambda u: np.exp(-u**2/2.0)
#         self._norm_factor = np.sqrt(2*np.pi)
# 
#     def evaluate_prop(self, x):  # returns values proportional to kde.
#         xh = np.array(x).reshape(-1)
#         res = np.zeros(len(xh))
#         if len(xh) > len(self.mus):  # loop over data
#             for mu, sigma, pi in zip(self.mus, self.sigmas, self.pis):
#                 res += pi/sigma*self.kernel((mu-xh)/sigma)
#         else:  # loop over x
#             for i, x_ in enumerate(xh):
#                 res[i] = np.sum(self.pis/self.sigmas*self.kernel((self.mus-x_)/self.sigmas))
#         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.sum(self.pis*norm.cdf((x - self.mus)/self.sigmas))
# 
#     def sample(self, n):
#         ns = np.random.multinomial(n, self.pis)
#         return np.hstack([np.random.randn(n_)*sigma+mu for n_, sigma, mu in zip(ns, self.sigmas, self.mus)])


import matplotlib.pyplot as plt
import numpy as np
        import time
        import time
from scipy.signal import argrelextrema
from scipy.stats import norm
from scipy.optimize import minimize, leastsq
from sklearn.neighbors import KernelDensity
from .util.ApproxGaussianKDE import ApproxGaussianKDE as KDE
from .util.GaussianMixture1d import GaussianMixture1d as GM
        from shoulder_distributions import bump_distribution
            return 0
        return mode_sizes[1]
    minval = 1
    i = 0
    max_init = 200
    while minval > mtol+m_err_tol:
        x0 = min+(max-min)*np.random.rand(1)
        res_minimizer = minimize(second_mode_size, x0=x0, bounds=[(min, max)], tol=1e-10)
        minval = res_minimizer.fun
        #print "(minval, res_minimizer.x, x0) = {}".format((minval, res_minimizer.x, x0))
        i += 1
        if i > max_init:
            raise ValueError('Unimodal function not found')
    min = res_minimizer.x
    return bisection_search_reference_distr(min, max, mtol, shoulder_ratio, shoulder_variance)


def bisection_search_reference_distr(amin, amax, mtol, shoulder_ratio, shoulder_variance, atol=1e-6):
    anew = (amin + amax)/2.0
    if amax - amin < atol:
        return amin  # ensure that resulting distribution is unimodal

    if shoulder_variance == 1:
        data = np.repeat([anew, 0], shoulder_ratio)
        passed = is_unimodal_kde(1, data, 0, mtol)
    else:
        m_err_tol = 1e-7
        covs = np.array([1, shoulder_variance])
        weights = np.array(shoulder_ratio, dtype=np.float)
        weights /= np.sum(weights)

# 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/util/GaussianMixture1d.py
# class GaussianMixture1d(object):
# 
#     def __init__(self, means, covs, weights):
#         self.mus = means
#         self.sigmas = np.sqrt(covs)
#         self.pis = weights
#         self.kernel = lambda u: np.exp(-u**2/2.0)
#         self._norm_factor = np.sqrt(2*np.pi)
# 
#     def evaluate_prop(self, x):  # returns values proportional to kde.
#         xh = np.array(x).reshape(-1)
#         res = np.zeros(len(xh))
#         if len(xh) > len(self.mus):  # loop over data
#             for mu, sigma, pi in zip(self.mus, self.sigmas, self.pis):
#                 res += pi/sigma*self.kernel((mu-xh)/sigma)
#         else:  # loop over x
#             for i, x_ in enumerate(xh):
#                 res[i] = np.sum(self.pis/self.sigmas*self.kernel((self.mus-x_)/self.sigmas))
#         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.sum(self.pis*norm.cdf((x - self.mus)/self.sigmas))
# 
#     def sample(self, n):
#         ns = np.random.multinomial(n, self.pis)
#         return np.hstack([np.random.randn(n_)*sigma+mu for n_, sigma, mu in zip(ns, self.sigmas, self.mus)])


# Path: modality/critical_bandwidth_fm.py
import matplotlib.pyplot as plt
import numpy as np
        import time
        import time
from scipy.signal import argrelextrema
from scipy.stats import norm
from scipy.optimize import minimize, leastsq
from sklearn.neighbors import KernelDensity
from .util.ApproxGaussianKDE import ApproxGaussianKDE as KDE
from .util.GaussianMixture1d import GaussianMixture1d as GM
        from shoulder_distributions import bump_distribution

            return 0
        return mode_sizes[1]
    minval = 1
    i = 0
    max_init = 200
    while minval > mtol+m_err_tol:
        x0 = min+(max-min)*np.random.rand(1)
        res_minimizer = minimize(second_mode_size, x0=x0, bounds=[(min, max)], tol=1e-10)
        minval = res_minimizer.fun
        #print "(minval, res_minimizer.x, x0) = {}".format((minval, res_minimizer.x, x0))
        i += 1
        if i > max_init:
            raise ValueError('Unimodal function not found')
    min = res_minimizer.x
    return bisection_search_reference_distr(min, max, mtol, shoulder_ratio, shoulder_variance)


def bisection_search_reference_distr(amin, amax, mtol, shoulder_ratio, shoulder_variance, atol=1e-6):
    anew = (amin + amax)/2.0
    if amax - amin < atol:
        return amin  # ensure that resulting distribution is unimodal

    if shoulder_variance == 1:
        data = np.repeat([anew, 0], shoulder_ratio)
        passed = is_unimodal_kde(1, data, 0, mtol)
    else:
        m_err_tol = 1e-7
        covs = np.array([1, shoulder_variance])
        weights = np.array(shoulder_ratio, dtype=np.float)
        weights /= np.sum(weights)

        passed = is_unimodal_from_kde(GM(np.array([0, anew]), covs, weights),