def calculateCentroidMeasurements(self):
        self.X[self.badFrames, :] = ma.masked
        if not self.useSmoothingFilterDerivatives:
            self.v[1:-1] = (self.X[2:, :] - self.X[0:-2])/(2.0/self.frameRate)
        else:
            # use a cubic polynomial filter to estimate the velocity
            self.v = ma.zeros(self.X.shape)
            halfWindow = int(np.round(self.filterWindow/2.*self.frameRate))
            for i in xrange(halfWindow, self.v.shape[0]-halfWindow):
                start = i-halfWindow
                mid = i
                finish = i+halfWindow+1
                if not np.any(self.X.mask[start:finish,:]):
                    px = np.polyder(np.polyfit(self.t[start:finish]-self.t[mid],
                                               self.X[start:finish, 0], 3))
                    py = np.polyder(np.polyfit(self.t[start:finish]-self.t[mid],
                                               self.X[start:finish, 1], 3))
                    self.v[i,:] = [np.polyval(px, 0), np.polyval(py, 0)]
                else:
                    self.v[i,:] = ma.masked

        self.s = ma.sqrt(ma.sum(ma.power(self.v, 2), axis=1))
        self.phi = ma.arctan2(self.v[:, 1], self.v[:, 0])
        self.t[self.badFrames] = ma.masked
        self.X[self.badFrames, :] = ma.masked
        self.v[self.badFrames, :] = ma.masked
        self.s[self.badFrames] = ma.masked
        self.phi[self.badFrames] = ma.masked

def sin_poly(P):
  poly_size = len(P)+1

  temp = np.poly1d([0])
  coeffCos = np.poly1d([0])
  coeffSin = np.poly1d([0])
  tmpCoeffCos = np.poly1d([0])
  tmpCoeffSin = np.poly1d([0])

  cos0 = math.cos(P[0])
  sin0 = math.sin(P[0])

  coeffSin[0] = 1.0
  temp[0] = sin0

  # compute the derivative of P
  dP = np.polyder(P)

  facti = 1
  for i in xrange(1,poly_size):
    facti *= i

    tmpCoeffCos = np.polyder(coeffCos) + coeffSin * dP
    tmpCoeffSin = np.polyder(coeffSin) - coeffCos * dP

    coeffCos = tmpCoeffCos
    coeffSin = tmpCoeffSin

    temp[i] = (coeffCos[0] * cos0 + coeffSin[0] * sin0)/float(facti)

  return temp

def deltafactor_polyfit(volumes, energies):
    """
    This is the routine used to compute V0, B0, B1 in the deltafactor code.

    Taken from deltafactor/eosfit.py
    """
    fitdata = np.polyfit(volumes**(-2./3.), energies, 3, full=True)
    ssr = fitdata[1]
    sst = np.sum((energies - np.average(energies))**2.)
    residuals0 = ssr/sst
    deriv0 = np.poly1d(fitdata[0])
    deriv1 = np.polyder(deriv0, 1)
    deriv2 = np.polyder(deriv1, 1)
    deriv3 = np.polyder(deriv2, 1)

    v0 = 0
    x = 0
    for x in np.roots(deriv1):
        if x > 0 and deriv2(x) > 0:
            v0 = x**(-3./2.)
            break
    else:
        raise EOSError("No minimum could be found")

    derivV2 = 4./9. * x**5. * deriv2(x)
    derivV3 = (-20./9. * x**(13./2.) * deriv2(x) - 8./27. * x**(15./2.) * deriv3(x))
    b0 = derivV2 / x**(3./2.)
    b1 = -1 - x**(-3./2.) * derivV3 / derivV2

    #print('deltafactor polyfit:')
    #print('e0, b0, b1, v0')
    #print(fitdata[0], b0, b1, v0)

    n = collections.namedtuple("DeltaFitResults", "v0 b0 b1 poly1d")
    return n(v0, b0, b1, fitdata[0])

    def fit_sjeos(self):
        """Calculate volume, energy, and bulk modulus.

        Returns the optimal volume, the minimum energy, and the bulk
        modulus.  Notice that the ASE units for the bulk modulus is
        eV/Angstrom^3 - to get the value in GPa, do this::

          v0, e0, B = eos.fit()
          print(B / kJ * 1.0e24, 'GPa')

        """

        fit0 = np.poly1d(np.polyfit(self.v**-(1 / 3), self.e, 3))
        fit1 = np.polyder(fit0, 1)
        fit2 = np.polyder(fit1, 1)

        self.v0 = None
        for t in np.roots(fit1):
            if isinstance(t, float) and t > 0 and fit2(t) > 0:
                self.v0 = t**-3
                break

        if self.v0 is None:
            raise ValueError('No minimum!')

        self.e0 = fit0(t)
        self.B = t**5 * fit2(t) / 9
        self.fit0 = fit0

        return self.v0, self.e0, self.B

def BM(energies):

    fitdata = np.polyfit(energies[:, 0] ** (-2.0 / 3.0), energies[:, 1], 3, full=True)
    ssr = fitdata[1]
    sst = np.sum((energies[:, 1] - np.average(energies[:, 1])) ** 2.0)
    residuals0 = ssr / sst
    deriv0 = np.poly1d(fitdata[0])
    deriv1 = np.polyder(deriv0, 1)
    deriv2 = np.polyder(deriv1, 1)
    deriv3 = np.polyder(deriv2, 1)

    volume0 = 0
    x = 0
    for x in np.roots(deriv1):
        if x > 0 and deriv2(x) > 0:
            volume0 = x ** (-3.0 / 2.0)
            break

    if volume0 == 0:
        print("Error: No minimum could be found")
        exit()

    derivV2 = 4.0 / 9.0 * x ** 5.0 * deriv2(x)
    derivV3 = -20.0 / 9.0 * x ** (13.0 / 2.0) * deriv2(x) - 8.0 / 27.0 * x ** (15.0 / 2.0) * deriv3(x)
    bulk_modulus0 = derivV2 / x ** (3.0 / 2.0)
    bulk_deriv0 = -1 - x ** (-3.0 / 2.0) * derivV3 / derivV2

    return volume0, bulk_modulus0, bulk_deriv0, residuals0

def findStepWithPoly(x,data,kind="stepUp",excludePoints=100,order=20,fitrange=100):
  """ Look for a step in the data
      Data is 1D array
      The 'kind' keywords should be either 'stepUp' or 'stepDown'
      the 'excludePoints' keyword is used to limit the search 'excludePoints' 
      away from the extremes
      The data are fit with a polynomial of order 'order', then the maximum
      (or minimum) of derivative is located.
      After this first attempt, the position is refined by a second polyfit
      in a range [-fitrange,+fitrange] around the first guess
  """
  if (kind == "stepUp"):
    use = "max"
  else:
    use = "min"
  poly    = np.polyfit(x,data,order)
  polyder = np.polyder(poly)
  x_poly1 = findDerPeak(x,np.polyval(polyder,x),use=use,excludePoints=excludePoints)
  # find closest
  idx = np.abs(x - x_poly1).argmin()
  idx = slice(idx-fitrange,idx+fitrange)
  poly    = np.polyfit(x[idx],data[idx],order)
  polyder = np.polyder(poly)
  x_poly2 = findDerPeak(x[idx],np.polyval(polyder,x[idx]),use=use,excludePoints=10)
  return x_poly2

def calc_max_curvature_point(x_arr, y_arr):
    rot_x = []
    rot_y = []

    # Rotate all of the x and y co-ordinates around the origin by 90deg
    for current_x, current_y in zip(x_arr, y_arr):
        out_x, out_y = rotate(current_x, current_y, 90)
        rot_x.append(out_x)
        rot_y.append(out_y)

    pol = get_best_poly_fit(rot_x, rot_y, 7)

    # Differentiate the polynomial
    deriv_one = numpy.polyder(pol)
    deriv_two = numpy.polyder(deriv_one)

    # Roots of the 1st derivative - that is, X-values for where 1st deriv = 0
    roots = numpy.roots(deriv_one)

    if roots.size == 1:
        result = rotate(roots, numpy.polyval(pol, roots), -90)
    else:
        prev_val = 0
        selected_root = -1

        for root in roots:
            val = numpy.polyval(deriv_two, root)
            if val > prev_val:
                selected_root = root

        result = rotate(selected_root, numpy.polyval(pol, selected_root), -90)

    return [float(result[0]), float(result[1])]

def solve_for_nearest( px,py,rx,ry ):
  dpx = polyder(px)
  dpy = polyder(py)
  cp = polymul( dpx, px ) + polymul( dpy, py )
  cp = polyadd( cp, -rx*dpx )      
  cp = polyadd( cp, -ry*dpy )
  t = roots(cp)
  t = real(t[isreal(t)])
  t = t[ (t>=0) * (t<=1) ]

  ##tt = linspace(0,1,100)
  ##from pylab import plot
  ##plot( polyval(px,tt), polyval(py,tt), 'k', hold = 0 )
  ##plot( [rx],[ry], 'r.' )
  ##plot( polyval(px,t[isreal(t)*(real(t)>=0)*(real(t)<=1)]),
  ##      polyval(py,t[isreal(t)*(real(t)>=0)*(real(t)<=1)]), 'o' )
  ##pdb.set_trace()

  if len(t):
    if len(t) == 1:
      return t[0]
    else:
      ux = polyval( px, t )
      uy = polyval( py, t )
      d = hypot( ux - rx, uy - ry )
      return t[ d==d.min() ][0]
  else:
    t = array([0.0,1.0])
    ux = polyval( px, t )
    uy = polyval( py, t )
    d = hypot( ux - rx, uy - ry )
    if d[0] < d[1]:
      return 0.0
    else:
      return 1.0

def solve_pathlength_func_bendrad(p, edge_y, edge_dydx, D, L0, BendRad, n_grid=100):
    """A helper function for solve_pathlength where a fixed minimum bend radius 
    is desired.
    
    Parameters
    ----------
    BendRad: float
        Enforced minimum bend radius.
    """
    pathlength_poly = np.poly1d(p)
    p_d = np.polyder(pathlength_poly)
    p_d_d = np.polyder(p_d)
    L = integrate.quad(polynomial_pathlength, 0, D, args=(p_d))

    # Find the minimum radius of curvature along the curve.
    # Given the the curve is non-differentiable, try brute-force with n_grid points
    # along the curve.
    x_vect = np.meshgrid(0, D, n_grid)

    a = bend_radius(x_vect, p_d, p_d_d)

    # Find the minimum radius of curvature.
    SmallCurve = np.min(a)

    retpar = [
        pathlength_poly(0) - edge_y[0],
        pathlength_poly(D) - edge_y[1],
        p_d(0) - edge_dydx[0],
        p_d(D) - edge_dydx[1],
        L[0] - L0,
        SmallCurve - BendRad,
    ]
    # print(retpar); pdb.set_trace()
    return retpar

    def splitBimodal(self, x, y, largepoly=30):
        p = np.polyfit(x, y, largepoly) # polynomial coefficients for fit

        extrema = np.roots(np.polyder(p))
        extrema = extrema[np.isreal(extrema)]
        extrema = extrema[(extrema - x[1]) * (x[-2] - extrema) > 0] # exclude the endpoints due false maxima during fitting
        try:
            root_vals = [sum([p[::-1][i]*(root**i) for i in range(len(p))]) for root in extrema]
            peaks = extrema[np.argpartition(root_vals, -2)][-2:] # find two peaks of bimodal distribution

            mid, = np.where((x - peaks[0])* (peaks[1] - x) > 0)
             # want data points between the peaks
        except:
            warnings.warn("Peak finding failed!")
            return None

        try:
            p_mid = np.polyfit(x[mid], y[mid], 2) # fit middle section to a parabola
            midpoint = np.roots(np.polyder(p_mid))[0]
        except:
            warnings.warn("Polynomial fit between peaks of distribution poorly conditioned. Falling back on using the minimum! May result in inaccurate split determination.")
            if len(mid) == 0:
                return None

            midx = np.argmin(y[mid])
            midpoint = x[mid][midx]

        return midpoint

def poly_sigma(time, y):
    """
    Calculates velocity and acceleration by fitting a polynomial over the prior N elements
    defined by poly_window. The sigma value defines an averaging window inside of poly_window
    in which all points inside of sigma window are averaged and treated as a single point for the
    polynomial fitting process.
    """
    y_d = np.zeros(time.shape)
    y_dd = np.zeros(time.shape)
    window = poly_window

    for i in range(window * sigma, time.shape[0]):
        y_history = y[i - window * sigma + 1:i + 1]
        y_hist_avg = np.mean(y_history.reshape(-1, sigma), axis=1)
        params = np.polyfit(
            x=time[i - window * sigma + 1:i + 1:sigma], y=y_hist_avg, deg=degree)

        p = np.poly1d(params)
        p_d = np.polyder(p)
        p_dd = np.polyder(p_d)

        y_d[i] = p_d(time[i])
        y_dd[i] = p_dd(time[i])

    y_d = y_d / encoder_resolution
    y_dd = y_dd / encoder_resolution

    return y_d, y_dd

    def _set_params(self):
        """
        Overriden to account for the fact the fit with volume**(2/3) instead
        of volume.
        """
        deriv0 = np.poly1d(self.eos_params)
        deriv1 = np.polyder(deriv0, 1)
        deriv2 = np.polyder(deriv1, 1)
        deriv3 = np.polyder(deriv2, 1)

        for x in np.roots(deriv1):
            if x > 0 and deriv2(x) > 0:
                v0 = x**(-3./2.)
                break
        else:
            raise EOSError("No minimum could be found")

        derivV2 = 4./9. * x**5. * deriv2(x)
        derivV3 = (-20./9. * x**(13./2.) * deriv2(x) - 8./27. *
                   x**(15./2.) * deriv3(x))
        b0 = derivV2 / x**(3./2.)
        b1 = -1 - x**(-3./2.) * derivV3 / derivV2

        # e0, b0, b1, v0
        self._params = [deriv0(v0**(-2./3.)), b0, b1, v0]

def fit_edge_hist(bins, counts, fwhm_guess=10.0):
    if len(bins) == len(counts)+1: bins = bins[:-1]+0.5*(bins[1]-bins[0]) # convert bin edge to bin centers if neccesary
    pfit = np.polyfit(bins, counts, 3)
    edgeGuess = np.roots(np.polyder(pfit,2))
    try:
        preGuessX, postGuessX = np.sort(np.roots(np.polyder(pfit,1)))
    except:
        raise ValueError("failed to generate guesses")
    use = bins>(edgeGuess+2*fwhm_guess)
    if np.sum(use)>4:
        pfit2 = np.polyfit(bins[use], counts[use],1)
        slope_guess = pfit2[0]
    else:
        slope_guess=1
    pGuess = np.array([edgeGuess, np.polyval(pfit,preGuessX), np.polyval(pfit,postGuessX),fwhm_guess,slope_guess],dtype='float64')

    try:
        pOut = curve_fit(edge_model, bins, counts, pGuess)
    except:
        return (0,0,0,0,0,0)
    (edgeCenter, preHeight, postHeight, fwhm, bgSlope) = pOut[0]
    model_counts = edge_model(bins, edgeCenter, preHeight, postHeight, fwhm, bgSlope)
    num_degree_of_freedom = float(len(bins)-1-5) # num points - 1 - number of fitted parameters
    chi2 = np.sum(((counts - model_counts)**2)/model_counts)/num_degree_of_freedom
    return (edgeCenter, preHeight, postHeight, fwhm, bgSlope, chi2)

def second_deriv_from_profile(xinterp,F):
    """Calculate second derivative at extrema"""
    dFdx = np.polyder(F,m=1)
    d2Fdx2 = np.polyder(F,m=2)
    minidx, maxidx = extrema_from_profile(xinterp,F)
    omegamin = d2Fdx2(xinterp[minidx])
    omegamax = d2Fdx2(xinterp[maxidx])
    return omegamin, omegamax

def polykruemm(p,x):
    y = np.polyval(p,x)
    dy = np.polyval(np.polyder(p),x)
    ddy = np.polyval(np.polyder(p,2),x)
    kappa = np.abs(ddy/(1+dy**2)**(3/2))
    mx = x - dy*(1+dy**2)/ddy
    my = y + (1+dy**2)/ddy
    return kappa,mx,my

def elastic_constant(name, atoms, e0, u0, mode, strain=0.007, debug=False):
    """Calculate an elastic constant.
    
    Parameters::
    
    name
        Name of the constant (a string). 
        
    atoms
        The atoms.
        
    e0
        Energy in the equilibrium configuration.
        
    u0
        Unit cell in the equilibrium configuration.
        
    mode
        The deformation mode as six numbers, giving weights
        to the six strain components.
    """
    strainlist = (-1.0, 0.5, 0, 0.5, 1.0)
    energies = []
    order = np.array([[0, 5, 4],
                      [5, 1, 3],
                      [4, 3, 2]])
    straintensor = np.array(mode)[order]
    if debug:
        print "%s unit strain tensor:" % (name,)
        print straintensor
    for s in strainlist:
        if s == 0:
            energies.append(e0)
        else:
            # Apply strain
            s = s * strain * straintensor + np.identity(3)
            atoms.set_cell(np.dot(u0, s), scale_atoms=True)
            energies.append(atoms.get_potential_energy())
    atoms.set_cell(u0, scale_atoms=True)  # Reset atoms
    fit0 = np.poly1d(np.polyfit(strain * np.array(strainlist), energies, 3))
    fit1 = np.polyder(fit0, 1)
    fit2 = np.polyder(fit1, 1)
    x0 = None
    for x in np.roots(fit1):
        if fit2(x) > 0:
            if x0 is not None:
                raise RuntimeError("More than two roots found.")
            assert x0 is None
            x0 = x
    if x0 is None:
        raise RuntimeError("No roots found.")
    if np.abs(x0) > 0.5 * strain:
        raise RuntimeError("Unreasonable root (%f): " % (x0,) +
                           "Maybe the system was not at the equilibrium configuration")
    if debug:
        print "Root:", x0
    value = fit2(x0)
    return value / atoms.get_volume()

def compute_join_length( px, py, tlow = 0.0, thigh = 1.0 ):
  from scipy.integrate import quad
  xp  = polyder( px, 1 )
  yp  = polyder( py, 1 )
  xp2 = polymul( xp, xp )
  yp2 = polymul( yp, yp )
  p   = polyadd( xp2, yp2 )
  integrand = lambda t: sqrt( polyval( p, t ) )
  return quad(integrand, tlow, thigh) [0]

    def analyse(self):
        OptimizeTask.analyse(self)

        for name, data in self.data.items():
            if 'distances' in data:
                distances = data['distances']
                energies = data['energies']
                fit0 = np.poly1d(np.polyfit(1 / distances, energies, 3))
                fit1 = np.polyder(fit0, 1)
                fit2 = np.polyder(fit1, 1)

                dmin = None
                for t in np.roots(fit1):
                    if t > 0 and fit2(t) > 0:
                        dmin = 1 / t
                        break

                if dmin is None:
                    raise ValueError('No minimum!')

                if abs(dmin) < min(distances) or abs(dmin) > max(distances):
                    raise ValueError('Fit outside of range! ' + \
                                      str(abs(dmin)) + ' not in ' + \
                                      str(distances))

                emin = fit0(t)
                k = fit2(t) * t**4
                m1, m2 = self.create_system(name).get_masses()
                m = m1 * m2 / (m1 + m2)
                hnu = units._hbar * 1e10 * sqrt(k / units._e / units._amu / m)

                data['minimum energy'] = emin
                self.results[name][1:] = [energies[2] - emin, dmin, 1000 * hnu]
            else:
                self.results[name].extend([None, None])

        for name, data in self.data.items():
            atoms = self.create_system(name)
            if len(atoms) == 1:
                self.results[name].extend([None, None])
                continue

            eatoms = 0.0
            for symbol in atoms.get_chemical_symbols():
                if symbol in self.data and symbol != name:
                    eatoms += self.data[symbol]['energy']
                else:
                    eatoms = None
                    break
            ea = None
            ea0 = None
            if eatoms is not None:
                ea = eatoms - data['energy']
                if 'minimum energy' in data:
                    ea0 = eatoms - data['minimum energy']
            self.results[name].extend([ea, ea0])

def solveEqnsSubOld(spline1,r1bar,r3bar,alphabar,offset):
    dspline1 = np.polyder(spline1)
    ddspline1 = np.polyder(dspline1)
    energy1 = morse(r1bar,alphabar,offset)
    denergy1 = morseD(r1bar,alphabar)
    f1 = np.polyval(spline1,r1bar) - energy1
    f2 = np.polyval(dspline1,r1bar) - denergy1
    f3 = np.polyval(spline1,r3bar)
    f4 = np.polyval(dspline1,r3bar)
    return np.array([f1,f2,f3,f4])

def root_angle_rad(w, side, dx, n=16):
    n = min(n, len(w.x) / 4)
    L = cumulative_path_length(w)
    tt = L / L.max()
    teval = tt[n] if side == 0 else tt[-n]
    px = np.polyfit(tt[n:-n], w.x[n:-n], 2)
    py = np.polyfit(tt[n:-n], w.y[n:-n], 2)
    xp = np.polyder(px, 1)
    yp = np.polyder(py, 1)
    return np.arctan2(dx * np.polyval(yp, teval), dx * np.polyval(xp, teval))