def period(self, line, edge="rising"):
        """
        Returns a dictionary with avg, min, max, and st of period for a line.
        """
        bit = self._line_to_bit(line)

        if edge.lower() == "rising":
            edges = self.get_rising_edges(bit)
        elif edge.lower() == "falling":
            edges = self.get_falling_edges(bit)

        if len(edges) > 2:

            timebase_freq = self.meta_data['ni_daq']['counter_output_freq']
            avg_period = np.mean(np.ediff1d(edges[1:]))/timebase_freq
            max_period = np.max(np.ediff1d(edges[1:]))/timebase_freq
            min_period = np.min(np.ediff1d(edges[1:]))/timebase_freq
            period_sd = np.std(avg_period)

        else:
            raise IndexError("Not enough edges for period: %i" % len(edges))

        return {
            'avg': avg_period,
            'max': max_period,
            'min': min_period,
            'sd': period_sd,
        }

    def turning_angle(self,consecutive=False):
        """ Returns collection of turning angles in the trace (in radians)
            If consecutive=True, only counts angles between consecutive flights
        """
        if self.N_flights < 2:
            return [False]
        else:
            if not consecutive:
                dp = np.array([g for g in self.iter_flights('StartEnd')]) #(pairs of points)
                angs = np.ediff1d(tuple(geom.angle_vect(g[0][:-1],g[1][:-1]) for g in dp))
            else:
                angs = []
                angs_d = []
                for f in self.iter_all():
                    if not isinstance(f,Stop):
                        angs_d.append(geom.angle_vect(f.StartEnd[0][:-1],f.StartEnd[-1][:-1]))
                    else:
                        if len(angs_d)>1:
	                        angs.extend(np.ediff1d(angs_d))
                        angs_d = []
                if len(angs_d)>1:
                    angs.extend(np.ediff1d(angs_d))
            for i,a in enumerate(angs):
                if a>np.pi:
                    angs[i] = - (2*np.pi-a)
                if a<-np.pi:
                    angs[i] = 2*np.pi+a
            return np.array(angs)

    def test_x2_in_x1(self):
        """
        x2 only has one bin, and it is surrounded by x1 range
        """
        # old size
        m = 4

        # new size
        n = 1

        # bin edges
        x_old = np.linspace(0., 1., m + 1)
        x_new = np.linspace(0.3, 0.65, n + 1)

        # some arbitrary distribution
        y_old = 1. + np.sin(x_old[:-1] * np.pi) / np.ediff1d(x_old)

        # rebin
        y_new = rebin.rebin(x_old, y_old, x_new,
                            interp_kind='piecewise_constant')

        # compute answer here to check rebin
        y_old_ave  = y_old / np.ediff1d(x_old)
        y_new_here = (y_old_ave[1] * (x_old[2] - x_new[0])
                      + y_old_ave[2] * (x_new[1] - x_old[2]) )

        assert_allclose(y_new, y_new_here)

    def test_x2_surrounds_x1(self):
        """
        x2 range surrounds x1 range
        """
        # old size
        m = 2

        # new size
        n = 3

        # bin edges
        x_old = np.linspace(0., 1., m + 1)
        x_new = np.linspace(-0.1, 1.2, n + 1)

        # some arbitrary distribution
        y_old = 1. + np.sin(x_old[:-1] * np.pi) / np.ediff1d(x_old)

        # rebin
        y_new = rebin.rebin(x_old, y_old, x_new,
                            interp_kind='piecewise_constant')

        # compute answer here to check rebin
        y_old_ave  = y_old / np.ediff1d(x_old)
        y_new_here = [y_old_ave[0] * (x_new[1] - 0.),
                      y_old_ave[0] * (x_old[1]- x_new[1])
                      + y_old_ave[1] * (x_new[2] - x_old[1]),
                      y_old_ave[1] * (x_old[-1] - x_new[-2])]

        assert_allclose(y_new, y_new_here)
        assert_allclose(y_new.sum(), y_old.sum())

def roc(x_vals, y_vals):
    speeds = []
    # x_diff = np.gradient(x_vals)
    x_diff = np.ediff1d(x_vals)
    # y_diff = np.gradient(y_vals)
    y_diff = np.ediff1d(y_vals)
    dy = y_diff / x_diff
    # dy_diff = np.gradient(dy)
    dy_diff = np.ediff1d(dy, to_begin=dy[1]-dy[0])
    d2y = dy_diff / x_diff
    R = np.power(1 + np.square(dy), 1.5) / d2y
    R[np.abs(R) > 32000] = 0
    print R
    for i in range(1,len(R)):
        if x_vals[i] - x_vals[i-1] < 0:
            if i == 1:
                R[i-1] = -R[i-1]
            R[i] = -R[i]
    for i in range(0, len(R)-1):
        dist = np.sqrt((x_vals[i+1] - x_vals[i])**2 + (y_vals[i+1] - y_vals[i])**2)
        theta = np.arccos(1 - dist**2 / (2*R[i]**2))
        if np.isnan(theta) or theta == 0:
            speeds.append(dist/time_interval)
        else:
            speeds.append(R[i]*theta / time_interval)
    R = R[:-1]
    return R, speeds

def speeds_by_arc_length_to_speeds_by_time(speeds_by_arc_length,
                                    arc_lengths, time_step_size):
    times_by_arc_length = speeds_by_arc_length_to_times_by_arc_length(
                                        speeds_by_arc_length, arc_lengths)
    trip_time = times_by_arc_length[-1]
    time_differences = np.ediff1d(times_by_arc_length)
    speed_differences = np.ediff1d(speeds_by_arc_length)
    slopes = np.divide(speed_differences, time_differences)
    num_time_steps = int(trip_time / time_step_size)
    time_steps = np.empty(num_time_steps)
    time_steps.fill(time_step_size)
    cumulative_time_steps = np.cumsum(time_steps)
    cumulative_time_steps = np.append(cumulative_time_steps, trip_time)
    selected_indices = np.searchsorted(times_by_arc_length,
                                       cumulative_time_steps)
    selected_times_by_arc_length = times_by_arc_length[selected_indices]
    selected_speeds_by_arc_length = speeds_by_arc_length[
                                                    selected_indices]
    selected_slopes = slopes[selected_indices - 1]
    excess_times = np.subtract(cumulative_time_steps,
                               selected_times_by_arc_length)
    excess_speeds = np.multiply(excess_times, selected_slopes)
    speeds_by_time = np.add(excess_speeds, 
                                selected_speeds_by_arc_length)
    return [speeds_by_time, cumulative_time_steps]

def write_griddata_fits(x0, x1, y, name_or_obj):
    """Write a 2-d grid of data to a FITS file

    The grid coordinates are encoded in the FITS-WCS header keywords.

    Parameters
    ----------
    x0 : numpy.ndarray
        1-d array.
    x1 : numpy.ndarray
        1-d array.
    y : numpy.ndarray
        2-d array of shape (len(x0), len(x1)).
    name_or_obj : str or file-like object
        Filename to write to or open file.
    """

    d0, d1 = np.ediff1d(x0), np.ediff1d(x1)
    if not (np.allclose(d0, d0[0]) and np.allclose(d1, d1[0])):
        raise ValueError('grid must be regularly spaced in both x0 and x1')
    if not (len(x0), len(x1)) == y.shape:
        raise ValueError('length of x0 and x1 do not match shape of y')

    w = wcs.WCS(naxis=2)
    w.wcs.crpix = [1, 1]
    w.wcs.crval = [x1[0], x0[0]]
    w.wcs.cdelt = [d1[0], d0[0]]
    hdu = fits.PrimaryHDU(y, header=w.to_header())
    hdu.writeto(name_or_obj)

def _jackknife_2d_random(rbins, box_size, jackknife_nside):
    def corner_area(x, y):
        a = math.sqrt(1.0-x*x)-y
        b = math.sqrt(1.0-y*y)-x
        theta = math.asin(math.sqrt(a*a+b*b)*0.5)*2.0
        return (a*b + theta - math.sin(theta))*0.5

    def segment_without_corner_area(x, r):
        half_chord = math.sqrt(1.0-x*x)
        return math.acos(x) - x*half_chord \
                - quad(corner_area, 0, min(half_chord, 1.0/r), (x,))[0]*r

    def overlapping_circular_areas(r):
        if r*r >= 2: return 1.0
        return (math.pi - quad(segment_without_corner_area, 0, min(1, 1.0/r), \
                (r,))[0]*4.0*r)*r*r

    overlapping_circular_areas_vec = np.vectorize(overlapping_circular_areas, \
            [float])

    side_length = box_size/float(jackknife_nside)
    square_area = 1.0/float(jackknife_nside*jackknife_nside)
    rbins_norm = rbins/side_length
    annulus_areas = np.ediff1d(overlapping_circular_areas_vec(rbins_norm))
    annulus_areas /= np.ediff1d(rbins_norm*rbins_norm)*math.pi
    return 1.0 - square_area * (2.0 - annulus_areas)

def findlevel(inwave, threshold, direction='both'):
    temp = inwave - threshold
    if (direction.find("up")+1):
        crossings = np.nonzero(np.ediff1d(np.sign(temp), to_begin=0)>0)
    elif (direction.find("down")+1):
        crossings = np.nonzero(np.ediff1d(np.sign(temp), to_begin=0)<0)
    else:
        crossings = np.nonzero(np.ediff1d(np.sign(temp), to_begin=0))
    return crossings[0][0]

def velocity(x, y):
	"""Returns array of velocity at each time step"""
	if isinstance(x, pd.Series):
		x = x.values
		y = y.values
	vx = np.ediff1d(x)
	vy = np.ediff1d(y)
	vel = np.sqrt( vx**2 + vy **2 ) # Pythagoras
	return vel

def ediff1d(ary, to_end=None, to_begin=None):
    if not isinstance(ary, Quantity):
        return np.ediff1d(ary, to_end, to_begin)

    return Quantity(
        np.ediff1d(ary.magnitude, to_end, to_begin),
        ary.dimensionality,
        copy=False
    )

def rr_1_update(rr_1, NFound, Found):
    if np.logical_and(NFound <= 7, NFound > 0):
        rr_1[0:NFound - 1] = np.ediff1d(Found[0:NFound, 0])
    elif NFound > 7:
        rr_1 = np.ediff1d((Found[NFound - 7:NFound - 1, 0]))

    rr_average_1 = np.mean(rr_1)

    return rr_1, rr_average_1

def check_sync(file_obj, wps, word_index, pattern_name):
    '''
    Check sync words in the array.

    Performs analysis of sync words in the file to check if there are no
    sync problems.
    
    :param file_obj: File object containing byte-aligned data.
    :param wps: Expected words per second
    :type wps: int
    :param word_index: First index of a sync word within the data.
    :type word_index: int
    :param pattern_name: Sync word pattern name, either 'Standard or Reverse'.
    :type pattern_name: str
    '''
    array = np.fromfile(file_obj, dtype=np.short)

    array &= 0xFFF

    s1 = np.array(
        np.nonzero(array == SYNC_PATTERNS[pattern_name][0])).flatten()
    s2 = np.array(
        np.nonzero(array == SYNC_PATTERNS[pattern_name][1])).flatten()
    s3 = np.array(
        np.nonzero(array == SYNC_PATTERNS[pattern_name][2])).flatten()
    s4 = np.array(
        np.nonzero(array == SYNC_PATTERNS[pattern_name][3])).flatten()
    syncs = np.concatenate((s1, s2, s3, s4))
    syncs = np.sort(syncs)
    syncs_i = iter(syncs)
    # print 'Sync words\n', syncs
    prev_sync = next(syncs_i)
    ix = prev_sync
    while ix < array.size:
        next_sync = ix + wps
        found = syncs == next_sync
        if np.any(found):
            ix = next_sync
        else:
            try:
                last_ix = ix
                ix = next(syncs_i)
                while ix < next_sync:
                    ix = next(syncs_i)
                logger.warning(
                    'Sync lost at word %d, next sync not found at %d, '
                    'found at %d instead', last_ix, next_sync, ix)
            except StopIteration:
                break

    syncs = np.ediff1d(syncs, to_begin=0, to_end=0)
    syncs[0] = syncs[1]
    syncs[-1] = syncs[-2]
    # print 'Distances\n', syncs
    syncs = np.ediff1d(syncs, to_begin=0, to_end=0)

def distancesFromArrays(xData, yData):
    """Returns distances between each points

    :param numpy.ndarray xData: X coordinate of points
    :param numpy.ndarray yData: Y coordinate of points
    :rtype: numpy.ndarray
    """
    deltas = numpy.dstack((
        numpy.ediff1d(xData, to_begin=numpy.float32(0.)),
        numpy.ediff1d(yData, to_begin=numpy.float32(0.))))[0]
    return numpy.cumsum(numpy.sqrt(numpy.sum(deltas ** 2, axis=1)))

def get_mids_and_lengths(x0, y0, x1, y1, gx, gy):
    """Return the midpoints and intersection lengths of a line and a grid.

    Parameters
    ----------
    x0,y0,x1,y1 : float
        Two points which define the line. Points must be outside the grid
    gx,gy : :py:class:`np.array`
        Defines positions for the gridlines

    Return
    ------
    xm,ym : :py:class:`np.array`
        Coordinates along the line within each intersected grid pixel.
    dist : :py:class:`np.array`
        Lengths of the line segments crossing each pixel

    """
    # avoid upper-right boundary errors
    if (x1 - x0) == 0:
        x0 += 1e-6
    if (y1 - y0) == 0:
        y0 += 1e-6

    # vector lengths (ax, ay)
    ax = (gx - x0) / (x1 - x0)
    ay = (gy - y0) / (y1 - y0)

    # edges of alpha (a0, a1)
    ax0 = min(ax[0], ax[-1])
    ax1 = max(ax[0], ax[-1])
    ay0 = min(ay[0], ay[-1])
    ay1 = max(ay[0], ay[-1])
    a0 = max(max(ax0, ay0), 0)
    a1 = min(min(ax1, ay1), 1)

    # sorted alpha vector
    cx = (ax >= a0) & (ax <= a1)
    cy = (ay >= a0) & (ay <= a1)
    alpha = np.sort(np.r_[ax[cx], ay[cy]])

    # lengths
    xv = x0 + alpha * (x1 - x0)
    yv = y0 + alpha * (y1 - y0)
    lx = np.ediff1d(xv)
    ly = np.ediff1d(yv)
    dist = np.sqrt(lx**2 + ly**2)

    # indexing
    mid = alpha[:-1] + np.ediff1d(alpha) / 2.
    xm = x0 + mid * (x1 - x0)
    ym = y0 + mid * (y1 - y0)

    return xm, ym, dist

	def check_move(self):
		'Checks to see if move available'
		if np.count_nonzero(self.data == 0):
			return True
		for row in self.data:
			if np.count_nonzero(np.ediff1d(row) == 0):
				return True
		for row in self.data.T:
			if np.count_nonzero(np.ediff1d(row) == 0):
				return True
		return False

def euc_dist(x, y):
    """
    Calculate Euclidean distances travelled by the trajectory at each time step.
    :param x: (array) x-coordinates
    :param y: (array) y-coordinates
    :return: Euclidean distance
    """
    xdist = np.ediff1d(x)**2
    ydist = np.ediff1d(y)**2
    dist = np.sqrt(xdist+ydist)
    return dist

def find_angles(core, var1, var2):
    """
    Calculates the "bends"/angles of variables graphed against each other
    (e.g. depth v age to look for sharp elbows)
    """
    x, y = graphlist(core, var1, var2)
    x1 = np.ediff1d(x)
    y1 = np.ediff1d(y)
    a = x1[:-1] ** 2 + y1[:-1] ** 2
    b = x1[1:] ** 2 + y1[1:] ** 2
    c = (x[2:] - x[:-2]) ** 2 + (y[2:] - y[:-2]) ** 2
    return np.degrees(np.arccos((a + b - c) / np.sqrt(4 * a * b)))

 def compute_derivative(self, x_vals, y_vals):
     y_diffs = np.ediff1d(y_vals)
     x_diffs = np.ediff1d(x_vals)
     quotients = np.divide(y_diffs, x_diffs)
     quotients_a = quotients[:-1]
     quotients_b = quotients[1:]
     mean_quotients = (quotients_a + quotients_b) / 2.0
     derivative = np.empty(x_vals.shape[0])
     derivative[1:-1] = mean_quotients
     derivative[0] = quotients[0]
     derivative[-1] = quotients[-1]
     return derivative

def jitter(vals):
    totVals = len(vals)
    B = numpy.unique(vals)
    N = WhereIs(vals, B)
    (n, B) = numpy.histogram(vals, B, new=False)
    fwdSpace = numpy.ediff1d(B, to_end=0.0)
    prevSpace = numpy.flipud(numpy.ediff1d(numpy.flipud(B), to_end=0.0))

    baseJit = ((fwdSpace[N] - prevSpace[N]) * numpy.random.rand(totVals)) - prevSpace[N]
    baseJit[n[N] <= 1] = 0.0

    return(vals + baseJit)