def smoothfit(x, y, smooth=0, res=1000):
    """
    Smooth data of the form f(x) = y with a spline
    """
    z = y.copy()
    w = isnan(z)
    z[w] = 0
    spl = UnivariateSpline(x, z, w=~w)
    spl.set_smoothing_factor(smooth)
    xs = linspace(min(x), max(x), res)
    ys = spl(xs)
    ys[ys < 0] = 0
    if w[0]:
        if len(where(~w)[0]):
            first = where(~w)[0][0]
            first = x[first]
            first = where(xs >= first)[0][0] - 1
            ys[:first] = nan
    if w[-1]:
        if len(where(~w)[0]):
            last = where(~w)[0][-1]
            last = x[last]
            last = where(xs >= last)[0][0] + 1
            ys[last:] = nan
    return xs, ys

def grad_dist(curdata,ax):
    #gaussian kernel density estimation           
    x=np.linspace(0,90,91)
    kde=gaussian_kde(curdata)
    line,=ax.plot(x,kde(x), '-', label=hl[j][0:-2]+"0")
    curcolor=plt.getp(line,'color')

##    #plotting histogram
##    ax.hist(curdata,bins=int(max(curdata))+1,
##             normed=True, histtype='step',
##             color=curcolor, linewidth=0.5)

    #defining splines for kd function and corresponding 1st and seconde derivatives
    x=np.linspace(0,90,901)
    s=UnivariateSpline(x,kde(x),s=0,k=3)
    s1=UnivariateSpline(x,s(x,1),s=0,k=3)
    s2=UnivariateSpline(x,s1(x,1),s=0,k=3)

    #identifying local maxima (where s1=0, s2<0, and s>0.005)
    maxima=s1.roots()[np.where(s2(s1.roots())<0)[0]]
    maxima=maxima[np.where(s(maxima)>0.005)[0]]
    #s_max=maxima[-1]
    s_max=maxima[np.argmax(s(maxima))]
    ax.plot(s_max, s(s_max),'o', color=curcolor)

    #identifying steepest segment after maxima (where x>=maxima, s1<0)
    x2=x[np.where(x>=s_max)[0]]
    slope=s1(x2)
    s1_min=x2[np.argmin(slope)]
    ax.plot(s1_min,s(s1_min),'o', color=curcolor)

    print round(s_max,1),round(s1_min,1)

    return round(s_max,1),round(s1_min,1)

def getCurvatureForPoints(arcLengthList, fx_s, fy_s, smoothing=None):
	x, x_, x__, y, y_, y__ = getFirstAndSecondDerivForTPoints(arcLengthList, fx_s, fy_s)
	curvature = abs(x_* y__ - y_* x__) / np.power(x_** 2 + y_** 2, 3 / 2)
	fCurvature = UnivariateSpline(arcLengthList, curvature, s=smoothing)
	dxcurvature = fCurvature.derivative(1)(arcLengthList)
	dx2curvature = fCurvature.derivative(2)(arcLengthList)
	return curvature, dxcurvature, dx2curvature

def get_derivatives(xs, ys, fd=False):
    """
    return the derivatives of y(x) at the points x
    if scipy is available a spline is generated to calculate the derivatives
    if scipy is not available the left and right slopes are calculated, if both exist the average is returned
    putting fd to zero always returns the finite difference slopes
    """
    try:
        if fd:
            raise SplineInputError('no spline wanted')
        if len(xs) < 4:
            er = SplineInputError('too few data points')
            raise er
        from scipy.interpolate import UnivariateSpline
        spline = UnivariateSpline(xs, ys)
        d = spline.derivative(1)(xs)
    except (ImportError, SplineInputError):
        d = []
        m, left, right = 0, 0, 0
        for n in range(0, len(xs), 1):
            try:
                left = (ys[n] - ys[n-1]) / (xs[n] - xs[n-1])
                m += 1
            except IndexError:
                pass
            try:
                right = (ys[n+1] - ys[n]) / (xs[n+1] - xs[n])
                m += 1
            except IndexError:
                pass
            d.append(left + right / m)
    return d

def imageSlice(image, xc, yc, width):
	ylen, xlen = image.shape

	xstart = max(0, xc - width / 2)
	xend = min(xlen, xc + width / 2)

	# I cannot get array slicing to work for the life of me
	slice = []
	for i in range(int(xstart), int(xend)):
		slice.append(image[yc, i])

	# https://stackoverflow.com/questions/10582795/finding-the-full-width-half-maximum-of-a-peak/10583774#10583774

	shiftedSlice = []

	halfMax = max(slice) / 2
	baseline = numpy.mean(image)

	for y in slice:
		shiftedSlice.append(y - halfMax - baseline)

	x = numpy.linspace(0, width, width)
	spline = UnivariateSpline(x, shiftedSlice, s=0)
	r1, r2 = spline.roots()
	#r1 = 0
	#r2 = 0

	return (slice, r2 - r1)

def smoothing(x,y,err=None,k=5,s=None,newx=None,derivative_order=0):
  # remove NaNs
  idx = np.isfinite(x) & np.isfinite(y)
  if idx.sum() != len(x): x=x[idx]; y=y[idx]

  # if we don't need to interpolate, use same x as input
  if newx is None: newx=x

  if err is None:
    w=None
  elif err == "auto":
    n=len(x)
    imin = int(max(0,n/2-20))
    imax = imin + 20
    idx = range(imin,imax)
    p = np.polyfit(x[idx],y[idx],2)
    e = np.std( y[idx] - np.polyval(p,x[idx] ) )
    w = np.ones_like(x)/e
  else:
    w=np.ones_like(x)/err
  from scipy.interpolate import UnivariateSpline
  if (s is not None):
    s = len(x)*s
  s = UnivariateSpline(x, y,w=w, k=k,s=s)
  if (derivative_order==0):
    return s(newx)
  else:
    try:
      len(derivative_order)
      return np.asarray([s.derivative(d)(newx) for d in derivative_order])
    except:
      return s.derivative(derivative_order)(newx)

def fwhm(x, y, k=10, ret_roots=False):
    """
    Determine full-with-half-maximum of a peaked set of points, x and y.

    Assumes that there is only one peak present in the dataset.  The function
    uses a spline interpolation with smoothing parameter k ('s' in scipy.interpolate.UnivariateSpline).
    """

    class MultiplePeaks(Exception):
        pass

    class NoPeaksFound(Exception):
        pass

    half_max = np.max(y) / 2.0
    s = UnivariateSpline(x, y - half_max, s=k)
    roots = s.roots()

    if len(roots) > 2:
        # Multiple peaks. Use the two that straddle the maximum value
        maxvel = x[np.argmax(y)]
        left_idx = np.argmin(maxvel - roots)
        right_idx = np.argmin(roots - maxvel)
        roots = np.array((roots[left_idx], roots[right_idx]))
    elif len(roots) < 2:
        raise NoPeaksFound("No proper peaks were found in the data set; likely "
                           "the dataset is flat (e.g. all zeros).")
    if ret_roots:
        return roots[0], roots[1]

    return abs(roots[1] - roots[0])

def xs_interp (inp_ene, inp_xs, inp_ene_interp, plot_cs):

    inp_ene = inp_ene # energies from Talys
    inp_xs  = inp_xs  # xs from talys
    inp_ene_interps = inp_ene_interp # energies for interpolation
    out_xs_A = []
    out_xs = np.array([]) # iterpolated xs
    plot_fig = plot_cs
    
    x_ene = np.linspace (0,660,3301)
    
    spl = UnivariateSpline(inp_ene, inp_xs, s = 0.25)
    y_xs = spl(x_ene)

    for inp_ene_interp in inp_ene_interps:
        out_xs_A.append(spl.__call__(inp_ene_interp))
    
    out_xs = np.append(out_xs, out_xs_A)

    # optional_plot

    if plot_fig:
        plt.plot (inp_ene, inp_xs, 'ro', ms = 5)
        plt.plot (x_ene, y_xs, lw = 3, c = 'g', alpha = 0.6)
        plt.plot (inp_ene_interps, out_xs, 'o', ms = 3)
        plt.show()
    
    return out_xs

def calc_conductance_curve(V_list,T,R_T,C_sigma):
    #test_voltages = arange(-v_max,v_max,v_step)
    test_currents = []
    for V in V_list:
        test_currents.append(calc_current(V,T,R_T,C_sigma))
        #print "V: %g, current %g"%(V,test_currents[-1])

    ## calc conductances manually
    #test_conductances = []
    #for idx,V in enumerate (test_currents[1:-2]):
    #    if idx==0:
    #        print idx
    #    test_conductances.append((test_currents[idx+2]-test_currents[idx])/(2.0*v_step))
    #
    #test_voltages_G = test_voltages[1:-2]

    #
    # SPLINE
    #
    spline = UnivariateSpline(V_list,test_currents,s=0)
    #print "test_conductances"
    #indices = [x for x, y in enumerate(col1) if (y >0.7 or y<-0.7)]
    test_conductances = []
    for v_iter in V_list:
        test_conductances.append(spline.derivatives(v_iter)[1])
    return test_conductances

def fwhm(x, y, bg=[0, 100, 150, 240]):
    """
    Evaluates the full width half maximum of y in units of x.
  
    Parameters
    ----------
    x : numpy.array
    y : numpy.array
    bg : list
      Background sampling limits
  
    Returns
    -------
    fwhm : number
      Full width half maximum
    """

    #  xnew = copy(x)
    #  ynew = copy(y)

    #  xc = x[(x>bg[0])&(x<bg[1]) | (x>bg[2])&(x<bg[3])]
    #  yc = y[(x>bg[0])&(x<bg[1]) | (x>bg[2])&(x<bg[3])]

    #  bgfit = polyfit(xc,yc,1)

    #  ynew = ynew-polyval(bgfit,xnew)

    xnew, ynew = rmbg(x, y, bg)

    f = UnivariateSpline(xnew, ynew / max(ynew) - 0.5, s=0)
    fwhm = f.roots()[1] - f.roots()[0]

    return fwhm

def integrated_rate_test(mx=100., annih_prod='BB'):
    # This currently doesn't work
    file_path = MAIN_PATH + "/Spectrum/"
    file_path += '{}'.format(int(mx)) + 'GeV_' + annih_prod + '_DMspectrum.dat'

    spectrum = np.loadtxt(file_path)
    imax = 0
    for i in range(len(spectrum)):
        if spectrum[i, 1] < 10 or i == (len(spectrum) - 1):
            imax = i
            break
    spectrum = spectrum[0:imax, :]
    Nevents = 10. ** 5.
    spectrum[:, 1] /= Nevents
    test = interp1d(np.log10(spectrum[:, 0] / mx), np.log10(mx * np.log(10.) * spectrum[:, 1]), kind='cubic', bounds_error=False, fill_value=0.)
    test2 = interp1d(spectrum[:, 0], spectrum[:, 0] * spectrum[:, 1], kind='cubic', bounds_error=False, fill_value=0.)
    e_gamma_tab = np.logspace(0., np.log10(spectrum[-1, 0]), 200)
    print np.column_stack((np.log10(spectrum[:, 0] / mx), np.log10(mx * np.log(10.) * spectrum[:, 1])))
    xtab = np.linspace(np.log10(1. / mx), 0., 200)
    ng2 = np.trapz(10.**test(xtab) / 10. ** xtab, xtab) / np.log(10.)
    mean_e2 = np.trapz(test2(e_gamma_tab), e_gamma_tab)
    rate_interp = UnivariateSpline(spectrum[:, 0], spectrum[:, 1])
    avg_e_interp = UnivariateSpline(spectrum[:, 0], spectrum[:, 0] * spectrum[:, 1])
    num_gamma = rate_interp.integral(1., spectrum[-1, 0])
    mean_e = avg_e_interp.integral(1., spectrum[-1, 0])


    print 'DM Mass: ', mx
    print 'Annihilation Products: ', annih_prod
    print 'Number of Gammas > 1 GeV: ', num_gamma, ng2
    print '<E> Gamma: ', mean_e, mean_e2

    return

    def response(self, disturbance_vector):
        """ Returns the response of the sensor due to a disturbance

            The acceleration imposed to the sensor is estimatad with the following equations:
                acceleration(t) = d2stress(t)/t2 * material_depth/material_prop['modulus']
            And the sensor response takes into account the frequency response, estimated by a normal curve
                freq_response(f) = norm(scale = self.bandwidth/2, loc=self.resonant_freq).pdf(f_array)
                freq_response(f) /= max(freq_response)
                response(t) = ifft(fft(acceleration) * freq_response)


            Args:
                disturbance_vector (list): list with a temporal array, in the 0 index, and a 
                stress array, in the 1 index.

            Returns:
                list: with two arrays, the temporal array and the voltage response array.
        """
        const = self.material_depth / self.material_prop['modulus']
        t_vector = disturbance_vector[0]
        # using the scipy UnivariateSpline to compute the second derivative
        data_spl = UnivariateSpline(t_vector, disturbance_vector[1], s=0, k=3)
        acceleration = data_spl.derivative(n=2)(t_vector) * const
        # we need to take the frequency response of the acceleration stimuli
        N = len(disturbance_vector[1])
        T = t_vector[1] - t_vector[0]
        f_array = np.fft.fftfreq(N, T)
        freq_acc = np.fft.fft(acceleration)
        # we need to apply a filter factor related to the frequency response of the sensor
        freq_response = self.frequency_response(N, (0, max(f_array)), mirror=True)[1]
        voltage = np.fft.ifft(freq_acc * freq_response) * self.sensitivity
        return voltage

class AlphaInterpolator(object):
    def __init__(self, a, x, y):

        # Drop NaN values to avoid fitpack errors
        self._data = pd.DataFrame(np.array([a, x, y]).T, columns=["a", "x", "y"])
        self._data.dropna(inplace=True)

        self._create_interpolating_polynomials()
        self._find_path_length()

    def _create_interpolating_polynomials(self):
        self.x_interp = UnivariateSpline(self._data.a, self._data.x, s=0)
        self.y_interp = UnivariateSpline(self._data.a, self._data.y, s=0)

    def _find_path_length(self):
        dx_interp = self.x_interp.derivative()
        dy_interp = self.y_interp.derivative()

        ts = np.linspace(0, 1, 200)
        line_length = cumtrapz(np.sqrt(dx_interp(ts) ** 2 + dy_interp(ts) ** 2), x=ts, initial=0.0)

        line_length /= line_length.max()

        # Here we invert the line_length (ts) function, in order to evenly
        # sample the pareto front
        self.l_interp = UnivariateSpline(line_length, ts, s=0)

    def sample(self, num):
        """ Return estimates of alpha values that evenly sample the pareto
        front """

        out = self.l_interp(np.linspace(0, 1, num))
        out[0] = 0.0
        out[-1] = 1.0
        return out

def kde_minmode(data,x,max_num_mode,min_mode_pdf):
    kde=gaussian_kde(data)
    f=kde.factor
    f_list=np.linspace(f,(data.max()-data.min()),100)
    s=UnivariateSpline(x,kde(x),s=0)
    s1=UnivariateSpline(x,s(x,1),s=0)
    s2=UnivariateSpline(x,s1(x,1),s=0)
    extrema=s1.roots()
    
    maxima=extrema[np.where((s2(extrema)<0)*(s(extrema)>=min_mode_pdf))]
    
    if len(maxima)>max_num_mode:
        for q in range(1,len(f_list)):
            f=f_list[q]
            kde2=gaussian_kde(data,bw_method=f)
            s=UnivariateSpline(x,kde2(x),s=0)
            s1=UnivariateSpline(x,s(x,1),s=0)
            s2=UnivariateSpline(x,s1(x,1),s=0)
            extrema=s1.roots()
            maxima=extrema[np.where((s2(extrema)<0)*(s(extrema)>=min_mode_pdf))]
            if len(maxima)<=max_num_mode:
##                print 'modes: ',maxima
                break
        kde=gaussian_kde(data,bw_method=f)
##    else:
##        print maxima

    return kde,maxima

    def smooth(self, genome, which_x, which_y):
        interpolationPointsQty = SMOOTHING_WINDOW
        which_y_InterpolationNeighborhood = interpolationPointsQty / 2
        minimunInterpolationNeighborhoodSize = interpolationPointsQty / 4

        if which_y - interpolationPointsQty / 2 < 0:
            interpolationPointsQty -= abs(which_y - which_y_InterpolationNeighborhood) * 2
            which_y_InterpolationNeighborhood = interpolationPointsQty / 2

        elif which_y + interpolationPointsQty / 2 > genome.getHeight() - 1:
            interpolationPointsQty -= (which_y + which_y_InterpolationNeighborhood - (genome.getHeight() - 1)) * 2
            which_y_InterpolationNeighborhood = interpolationPointsQty / 2

        if which_y_InterpolationNeighborhood >= minimunInterpolationNeighborhoodSize:
            x = np.ndarray(interpolationPointsQty)
            y = np.ndarray(interpolationPointsQty)

            for k in xrange(interpolationPointsQty):
                poseToSmooth = which_y - which_y_InterpolationNeighborhood + k
                x[k] = poseToSmooth
                y[k] = genome[poseToSmooth][which_x]

            spl = UnivariateSpline(x, y)
            spl.set_smoothing_factor(SPLINE_SMOOTHING_FACTOR_SPLINE/10)

            for k in xrange(interpolationPointsQty):
                if y[k] != sysConstants.JOINT_SENTINEL:
                    newValue = spl(int(x[k]))
                    genome.setItem(int(x[k]), which_x, newValue)

def halbwertsbreite(x, y):
    spline = UnivariateSpline(x, y-np.max(y)/2, s=0)
    r1, r2 = spline.roots() # find the roots

    lambda1 = 2*d*np.sin(np.deg2rad(r1))
    lambda2 = 2*d*np.sin(np.deg2rad(r2))
    E1 = h*c/lambda1
    E2 = h*c/lambda2
    DE = E1 - E2
    print ('Halbwertswinkel: {0:.5e} deg, {1:.5e} deg'.format(r1, r2))
    print ('Halbwertsbreite: {0:.5e}'.format(np.abs(r1-r2)))
    print (u'Energieaufloesung: {0:.5e} eV'.format(DE))

    xnew = np.linspace(min(x), max(x))
    ynew = spline(xnew)

    plt.plot(x, y, 'rx', label='Messdaten')
    plt.plot(xnew, ynew+np.max(y)/2,'b-', label='Interpolation')
    plt.axvline(r1)
    plt.axvline(r2)

    plt.grid()
    plt.legend()
    plt.xlabel("doppelter Kristallwinkel in Grad")
    plt.ylabel(u"Zählrate")

	def __init__(self, yp, workdir, scale, sm=200):
		'''
		Constructor
		'''
		yp = np.array(yp)
		self.l = len(yp)/2
		self.xPos = (self.l-2)/2 #fPos = (self.l-2)/2 + 2
		tnsc = 2/scale
		print tnsc
		plt.rcParams['font.size'] = 24
		plt.rcParams['lines.linewidth'] = 2.4
		self.workdir = workdir
		
		avProfilePoints = yp[:self.l]
		self.avx = np.append(np.append([0], np.sort(np.tanh(tnsc*avProfilePoints[:self.xPos]))),[1])
		self.av = avProfilePoints[self.xPos:]
		
		sigmaProfilePoints = yp[self.l:]
		self.sigmax = np.append(np.append([0], np.sort(np.tanh(tnsc*sigmaProfilePoints[:self.xPos]))),[1])
		self.sigma = sigmaProfilePoints[self.xPos:]
		
		self.m = UnivariateSpline(self.avx, self.av)
		print "Created spline with " + str(len(self.m.get_knots())) + " knots"

		self.s = UnivariateSpline(self.sigmax, self.sigma)
		print "Created spline with " + str(len(self.s.get_knots())) + " knots"

    def get_t_for_vols(self, vols, t_max=1000):
        """
        Find the temperatures corresponding to a specific volume.
        The search is performed interpolating the V(T) dependence with a spline and
        finding the roots with of V(t) - v.
        It may return more than one temperature for a volume in case of non monotonic behavior.

        Args:
            vols: list of volumes
            t_max: maximum temperature considered for the fit

        Returns:
            A list of lists of temperatures. For each volume more than one temperature can
            be identified.
        """

        if not isinstance(vols, (list, tuple, np.ndarray)):
            vols = [vols]


        f = self.fit_energies(0, t_max, t_max+1)

        temps = []
        for v in vols:
            spline = UnivariateSpline(f.temp, f.min_vol - v, s=0)
            temps.append(spline.roots())

        return temps

def find_extrema_spline(x,y,k=4,s=0,**kwargs):
    """
    find local extrema of y(x) by taking derivative of spline

    Parameters
    ----------

    x,y : array-like
      find extrema of y(x)

    k : int
      order of spline interpolation [must be >3]

    s : number
      s parameter sent to scipy spline interpolation (used for smoothing)

    **kwargs : extra arguments to UnivariateSpline


    Returns
    -------
    sp : UnivariateSpline object

    x_max,y_max : array
      value of x and y at extrema(s)
    """


    sp = UVSpline(x,y,k=k,s=s,**kwargs)

    x_max = sp.derivative().roots()

    y_max = sp(x_max)

    return sp,x_max,y_max

	def getspline_Sold(self):
		"""Cubic spline interpolation of entropy and convective velocity.
		"""
		want = self.mass < max(self.mass)*self.mass_cut
		S_old = UnivariateSpline(self.mass[want], self.Sgas[want], k=self.spline_k, s=self.spline_s, ext=self.spline_ext)
		dS_old = S_old.derivative()
		vconv_Sold = UnivariateSpline(self.mass[want], self.vconv[want], k=self.spline_k, s=self.spline_s, ext=self.spline_ext)
		return [S_old, dS_old, vconv_Sold]