Python numpy.diff() Examples

def disassemble_sparse_matrix(matrix: scipy.sparse.spmatrix) -> dict:
    """
    Transform a scipy.sparse matrix into the serializable collection of \
    :class:`numpy.ndarray`-s. :func:`assemble_sparse_matrix()` does the inverse.

    :param matrix: :mod:`scipy.sparse` matrix; csr, csc and coo formats are \
                   supported.
    :return: :class:`dict` with "shape", "format" and "data" - :class:`tuple` \
             of :class:`numpy.ndarray`.
    """
    fmt = matrix.getformat()
    if fmt not in ("csr", "csc", "coo"):
        raise ValueError("Unsupported scipy.sparse matrix format: %s." % fmt)
    result = {
        "shape": matrix.shape,
        "format": fmt
    }
    if isinstance(matrix, (scipy.sparse.csr_matrix, scipy.sparse.csc_matrix)):
        lengths = numpy.concatenate(([0], numpy.diff(matrix.indptr)))
        result["data"] = [matrix.data, squeeze_bits(matrix.indices), squeeze_bits(lengths)]
    elif isinstance(matrix, scipy.sparse.coo_matrix):
        result["data"] = [matrix.data, (squeeze_bits(matrix.row), squeeze_bits(matrix.col))]
    return result 

def _ecg_findpeaks_pantompkins(signal, sampling_rate=1000):
    """From https://github.com/berndporr/py-ecg-detectors/

    - Jiapu Pan and Willis J. Tompkins. A Real-Time QRS Detection Algorithm.
    In: IEEE Transactions on Biomedical Engineering BME-32.3 (1985), pp. 230–236.

    """
    diff = np.diff(signal)

    squared = diff * diff

    N = int(0.12 * sampling_rate)
    mwa = _ecg_findpeaks_MWA(squared, N)
    mwa[: int(0.2 * sampling_rate)] = 0

    mwa_peaks = _ecg_findpeaks_peakdetect(mwa, sampling_rate)

    mwa_peaks = np.array(mwa_peaks, dtype="int")
    return mwa_peaks


# =============================================================================
# Hamilton (2002)
# ============================================================================= 

def _ecg_rsa_cycles(signals):
    """Extract respiratory cycles."""
    inspiration_onsets = np.intersect1d(
        np.where(signals["RSP_Phase"] == 1)[0], np.where(signals["RSP_Phase_Completion"] == 0)[0], assume_unique=True
    )

    expiration_onsets = np.intersect1d(
        np.where(signals["RSP_Phase"] == 0)[0], np.where(signals["RSP_Phase_Completion"] == 0)[0], assume_unique=True
    )

    cycles_length = np.diff(inspiration_onsets)

    return {
        "RSP_Inspiration_Onsets": inspiration_onsets,
        "RSP_Expiration_Onsets": expiration_onsets,
        "RSP_Cycles_Length": cycles_length,
    } 

def make_sessions(data, session_th=30 * 60, is_ordered=False, user_key='user_id', item_key='item_id', time_key='ts'):
    """Assigns session ids to the events in data without grouping keys"""
    if not is_ordered:
        # sort data by user and time
        data.sort_values(by=[user_key, time_key], ascending=True, inplace=True)
    # compute the time difference between queries
    tdiff = np.diff(data[time_key].values)
    # check which of them are bigger then session_th
    split_session = tdiff > session_th
    split_session = np.r_[True, split_session]
    # check when the user chenges is data
    new_user = data['user_id'].values[1:] != data['user_id'].values[:-1]
    new_user = np.r_[True, new_user]
    # a new sessions stars when at least one of the two conditions is verified
    new_session = np.logical_or(new_user, split_session)
    # compute the session ids
    session_ids = np.cumsum(new_session)
    data['session_id'] = session_ids
    return data 

def allVL1leq(n, L1):
    res = []
    for i_L1 in range(L1+1):
        # Chose (n-1) the splitting points of the array [0:(n+L1)]
        s = vchoosek(np.linspace(1,n+i_L1-1,n+i_L1-1)*np.ones([1,n+i_L1-1]),n-1) # 1:n+L1-1
        
        m = s.shape[0]

        s1 = np.zeros([m,1])
        s2 = (n+i_L1)+s1

        v = np.diff(np.hstack([s1, s, s2])) 
        v = v-1

        if i_L1 == 0:         
            res = v
        else:
            res = np.vstack([res,v])
        
    return res.astype(int) 

def test_simpSample_trivial(self):
        """ Test simple rejection sampler with trivial inputs

        Test method: set up sampling with equal upper and lower bounds
        """
        
        ulim = [0,1]
        ufun = lambda x: 1.0/np.diff(ulim)

        n = 10000
        
        for mod in self.mods:
            print('Testing trivial input for sampler: %s'%mod.__name__)
            sampler = mod(ufun,0.5,0.5)
            sample = sampler(n)
            
            self.assertEqual(len(sample),n,'Sampler %s does not return all same value'%mod.__name__)
            self.assertTrue(np.all(sample == 0.5),'Sampler %s does not return all values at 0.5'%mod.__name__) 

def hfd(X, kmax):
    N = len(X)
    Nm1 = float(N - 1)
    L = np.empty((kmax,))
    L[0] = np.sum(abs(np.diff(X, n=1))) # shortcut :)
    for k in xrange(2, kmax + 1):
        Lmks = np.empty((k,))
        for m in xrange(1, k + 1):
            i_end = (N - m) / k # int
            Lmk_sum = np.sum(abs(np.diff(X[np.arange(m - 1, m + (i_end + 1) * k - 1, k)], n=1)))
            Lmk = Lmk_sum * Nm1 / (i_end * k)
            Lmks[m-1] = Lmk

        L[k - 1] = np.mean(Lmks)

    a = np.empty((kmax, 2))
    a[:, 0] = np.log(1.0 / np.arange(1.0, kmax + 1.0))
    a[:, 1] = 1.0

    b = np.log(L)

    # find x by solving for ax = b
    x, residues, rank, s = np.linalg.lstsq(a, b)
    return x[0] 

def calc_axon_contribution(self, axons):
        xyret = np.column_stack((self.grid.xret.ravel(),
                                 self.grid.yret.ravel()))
        # Only include axon segments that are < `max_d2` from the soma. These
        # axon segments will have `sensitivity` > `self.min_ax_sensitivity`:
        max_d2 = -2.0 * self.axlambda ** 2 * np.log(self.min_ax_sensitivity)
        axon_contrib = []
        for xy, bundle in zip(xyret, axons):
            idx = np.argmin((bundle[:, 0] - xy[0]) ** 2 +
                            (bundle[:, 1] - xy[1]) ** 2)
            # Cut off the part of the fiber that goes beyond the soma:
            axon = np.flipud(bundle[0: idx + 1, :])
            # Add the exact location of the soma:
            axon = np.insert(axon, 0, xy, axis=0)
            # For every axon segment, calculate distance from soma by
            # summing up the individual distances between neighboring axon
            # segments (by "walking along the axon"):
            d2 = np.cumsum(np.diff(axon[:, 0], axis=0) ** 2 +
                           np.diff(axon[:, 1], axis=0) ** 2)
            idx_d2 = d2 < max_d2
            sensitivity = np.exp(-d2[idx_d2] / (2.0 * self.axlambda ** 2))
            idx_d2 = np.insert(idx_d2, 0, False)
            contrib = np.column_stack((axon[idx_d2, :], sensitivity))
            axon_contrib.append(contrib)
        return axon_contrib 

def test_gen_sma(self):
        r"""Test gen_sma method.

        Approach: Ensures the output is set, of the correct type, length, and units.
        Check that they are in the correct range and follow the distribution.
        """

        plan_pop = self.fixture
        n = 10000
        sma = plan_pop.gen_sma(n)

        # ensure the units are length
        self.assertEqual((sma/u.km).decompose().unit, u.dimensionless_unscaled)
        # sma > 0
        self.assertTrue(np.all(sma.value >= 0))
        # sma >= arange[0], sma <= arange[1]
        self.assertTrue(np.all(sma - plan_pop.arange[0] >= 0))
        self.assertTrue(np.all(plan_pop.arange[1] - sma >= 0))

        h = np.histogram(sma.to('AU').value,100,density=True)
        hx = np.diff(h[1])/2.+h[1][:-1]
        hp = plan_pop.dist_sma(hx)

        chi2 = scipy.stats.chisquare(h[0],hp)
        self.assertGreaterEqual(chi2[1],0.95) 

def test_gen_plan_params(self):
        r"""Test generated planet parameters:
        Expected: all 1 R_E, all p = 0.67, e = 0, and uniform a,e in arange,erange
        """

        obj = EarthTwinHabZone2(constrainOrbits=False,erange=[0.1,0.5],**self.spec)

        x = 10000

        a, e, p, Rp = obj.gen_plan_params(x)
        
        assert(np.all(p == 0.367))
        assert(np.all(Rp == 1.0*u.R_earth))

        for param,param_range in zip([a.value,e],[obj.arange.value,obj.erange]):
            h = np.histogram(param,100,density=True)
            chi2 = scipy.stats.chisquare(h[0],[1.0/np.diff(param_range)[0]]*len(h[0]))
            self.assertGreater(chi2[1], 0.95) 

def test_gen_plan_params(self):
        r"""Test generated planet parameters:
        Expected: all 1 R_E, all p = 0.67, e = 0, and uniform a in arange
        """

        obj = EarthTwinHabZone1(**self.spec)

        x = 10000

        a, e, p, Rp = obj.gen_plan_params(x)
        
        assert(np.all(e == 0))
        assert(np.all(p == 0.367))
        assert(np.all(Rp == 1.0*u.R_earth))

        h = np.histogram(a.to('AU').value,100,density=True)
        chi2 = scipy.stats.chisquare(h[0],[1.0/np.diff(obj.arange.to('AU').value)[0]]*len(h[0]))
        self.assertGreater(chi2[1], 0.95) 

def test_gen_sma(self):
        r"""Test gen_sma method.

        Approach: Ensures the output is set, of the correct type, length, and units.
        Check that they are in the correct range and follow the distribution.
        """

        plan_pop = self.fixture
        n = 10000
        sma = plan_pop.gen_sma(n)

        # ensure the units are length
        self.assertEqual((sma/u.km).decompose().unit, u.dimensionless_unscaled)
        # sma > 0
        self.assertTrue(np.all(sma.value >= 0))
        # sma >= arange[0], sma <= arange[1]
        self.assertTrue(np.all(sma - plan_pop.arange[0] >= 0))
        self.assertTrue(np.all(plan_pop.arange[1] - sma >= 0))

        h = np.histogram(sma.to('AU').value,100,density=True)
        hx = np.diff(h[1])/2.+h[1][:-1]
        hp = plan_pop.dist_sma(hx)

        chi2 = scipy.stats.chisquare(h[0],hp)
        self.assertGreaterEqual(chi2[1],0.95) 

def test_AxonMapModel_calc_axon_contribution(engine):
    model = AxonMapModel(xystep=2, engine=engine, n_axons=10,
                         xrange=(-20, 20), yrange=(-15, 15),
                         axons_range=(-30, 30))
    model.build()
    xyret = np.column_stack((model.spatial.grid.xret.ravel(),
                             model.spatial.grid.yret.ravel()))
    bundles = model.spatial.grow_axon_bundles()
    axons = model.spatial.find_closest_axon(bundles)
    contrib = model.spatial.calc_axon_contribution(axons)

    # Check lambda math:
    for ax, xy in zip(contrib, xyret):
        axon = np.insert(ax, 0, list(xy) + [0], axis=0)
        d2 = np.cumsum(np.diff(axon[:, 0], axis=0) ** 2 +
                       np.diff(axon[:, 1], axis=0) ** 2)
        sensitivity = np.exp(-d2 / (2.0 * model.spatial.axlambda ** 2))
        npt.assert_almost_equal(sensitivity, ax[:, 2]) 

def mask_to_rle(img, mask_value=255, transpose=True):
    img = np.int32(img)
    if transpose:
      img = img.T
    img = img.flatten()
    img[img == mask_value] = 1
    pimg = np.pad(img, 1, mode='constant')
    diff = np.diff(pimg)
    starts = np.where(diff == 1)[0]
    ends = np.where(diff == -1)[0]
    rle = []
    previous_end = 0
    for start, end in zip(starts, ends):
      relative_start = start - previous_end
      length = end - start
      previous_end = end
      rle.append(str(relative_start))
      rle.append(str(length))
    if len(rle) == 0:
      return "-1"
    return " ".join(rle) 

def _rsp_findpeaks_biosppy(rsp_cleaned, sampling_rate):

    extrema = _rsp_findpeaks_extrema(rsp_cleaned)
    extrema, amplitudes = _rsp_findpeaks_outliers(rsp_cleaned, extrema, amplitude_min=0)

    peaks, troughs = _rsp_findpeaks_sanitize(extrema, amplitudes)

    # Apply minimum period outlier-criterion (exclude inter-breath-intervals
    # that produce breathing rate larger than 35 breaths per minute.
    outlier_idcs = np.where((np.diff(peaks) / sampling_rate) < 1.7)[0]

    peaks = np.delete(peaks, outlier_idcs)
    troughs = np.delete(troughs, outlier_idcs)

    info = {"RSP_Peaks": peaks, "RSP_Troughs": troughs}
    return info 

def assemble_sparse_matrix(subtree: dict) -> scipy.sparse.spmatrix:
    """
    Transform a dictionary with "shape", "format" and "data" into the \
    :mod:`scipy.sparse` matrix. \
    Opposite to :func:`disassemble_sparse_matrix()`.

    :param subtree: :class:`dict` which describes the :mod:`scipy.sparse` \
                    matrix.
    :return: :mod:`scipy.sparse` matrix of the specified format.
    """
    matrix_class = getattr(scipy.sparse, "%s_matrix" % subtree["format"])
    if subtree["format"] in ("csr", "csc"):
        indptr = subtree["data"][2]
        if indptr[-1] != subtree["data"][0].shape[0]:
            # indptr is diff-ed
            subtree["data"][2] = indptr.cumsum()
    matrix = matrix_class(tuple(subtree["data"]), shape=subtree["shape"])
    return matrix 

def test_cmag(self):
        '''
        test_cmag() ensures that the neuropythy.vision cortical magnification function is working.
        '''
        import neuropythy.vision as vis
        logging.info('neuropythy: Testing areal cortical magnification...')
        dset = ny.data['benson_winawer_2018']
        sub = dset.subjects['S1202']
        hem = [sub.lh, sub.rh][np.random.randint(2)]
        cm = vis.areal_cmag(hem.midgray_surface, 'prf_',
                            mask=('inf-prf_visual_area', 1),
                            weight='prf_variance_explained')
        # cmag should get smaller in general
        ths = np.arange(0, 2*np.pi, np.pi/3)
        es = [0.5, 1, 2, 4]
        x = np.diff([np.mean(cm(e*np.cos(ths), e*np.sin(ths))) for e in es])
        self.assertTrue((x < 0).all()) 

def subcurve(self, t0, t1):
        '''
        curve.subcurve(t0, t1) yields a curve-spline object that is equivalent to the given
          curve but that extends from curve(t0) to curve(t1) only.
        '''
        # if t1 is less than t0, then we want to actually do this in reverse...
        if t1 == t0: raise ValueError('Cannot take subcurve of a point')
        if t1 < t0:
            tt = self.curve_length()
            return self.reverse().subcurve(tt - t0, tt - t1)
        idx = [ii for (ii,t) in enumerate(self.t) if t0 < t and t < t1]
        pt0 = self(t0)
        pt1 = self(t1)
        coords = np.vstack([[pt0], self.coordinates.T[idx], [pt1]])
        ts = np.concatenate([[t0], self.t[idx], [t1]])
        dists  = None if self.distances is None else np.diff(ts)
        return CurveSpline(
            coords.T,
            order=self.order,
            smoothing=self.smoothing,
            periodic=False,
            distances=dists,
            meta_data=self.meta_data) 

def splrep(coordinates, t, order, weights, smoothing, periodic):
        from scipy import interpolate
        (x,y) = coordinates
        # we need to skip anything where t[i] and t[i+1] are too close
        wh = np.where(np.isclose(np.diff(t), 0))[0]
        if len(wh) > 0:
            (t,x,y) = [np.array(u) for u in (t,x,y)]
            ii = np.arange(len(t))
            for i in reversed(wh):
                ii[i+1:-1] = ii[i+2:]
                for u in (t,x,y):
                    u[i] = np.mean(u[i:i+2])
            ii = ii[:-len(wh)]
            (t,x,y) = [u[ii] for u in (t,x,y)]
        xtck = interpolate.splrep(t, x, k=order, s=smoothing, w=weights, per=periodic)
        ytck = interpolate.splrep(t, y, k=order, s=smoothing, w=weights, per=periodic)
        return tuple([tuple([pimms.imm_array(u) for u in tck])
                      for tck in (xtck,ytck)]) 

def get_stages_from_blocks(self, widths):
        """Gets widths/stage_blocks of network at each stage.

        Args:
            widths (list[int]): Width in each stage.

        Returns:
            tuple(list): width and depth of each stage
        """
        width_diff = [
            width != width_prev
            for width, width_prev in zip(widths + [0], [0] + widths)
        ]
        stage_widths = [
            width for width, diff in zip(widths, width_diff[:-1]) if diff
        ]
        stage_blocks = np.diff([
            depth for depth, diff in zip(range(len(width_diff)), width_diff)
            if diff
        ]).tolist()
        return stage_widths, stage_blocks 

def test_disassemble_sparse_matrix(self):
        arr = numpy.zeros((10, 10), dtype=numpy.float32)
        numpy.random.seed(0)
        arr[numpy.random.randint(0, 10, (50, 2))] = 1
        mat = csr_matrix(arr)
        dis = disassemble_sparse_matrix(mat)
        self.assertIsInstance(dis, dict)
        self.assertIn("shape", dis)
        self.assertIn("format", dis)
        self.assertIn("data", dis)
        self.assertEqual(dis["shape"], arr.shape)
        self.assertEqual(dis["format"], "csr")
        self.assertIsInstance(dis["data"], list)
        self.assertEqual(len(dis["data"]), 3)
        self.assertTrue((dis["data"][0] == mat.data).all())
        self.assertTrue((dis["data"][1] == mat.indices).all())
        self.assertTrue((dis["data"][2] == [0] + list(numpy.diff(mat.indptr))).all())
        self.assertEqual(dis["data"][2].dtype, numpy.uint8) 

def timeDomain(NN):
    
    L = len(NN)    
    ANN = np.mean(NN)
    SDNN = np.std(NN)
    SDSD = np.std(np.diff(NN))    
    NN50 = len(np.where(np.diff(NN) > 0.05)[0])    
    pNN50 = NN50/L    
    NN20 = len(np.where(np.diff(NN) > 0.02)[0])
    pNN20 = NN20/L
    rMSSD = np.sqrt((1/L) * sum(np.diff(NN) ** 2))        
    MedianNN = np.median(NN)
    
    timeDomainFeats = {'ANN': ANN, 'SDNN': SDNN,
                       'SDSD': SDSD, 'NN50': NN50,
                       'pNN50': pNN50, 'NN20': NN20,
                       'pNN20': pNN20, 'rMSSD': rMSSD,
                       'MedianNN':MedianNN}
                       
    return timeDomainFeats 

def _eog_clean_kong1998(eog_signal, sampling_rate=1000):
    """Kong, X., & Wilson, G.

    F. (1998). A new EOG-based eyeblink detection algorithm. Behavior Research Methods, Instruments, & Computers,
    30(4), 713-719.

    """
    #  The order E should be less than half of the expected eyeblink duration. For example, if
    # the expected blink duration is 200 msec (10 samples with a sampling rate of 50 Hz), the
    # order E should be less than five samples.
    eroded = scipy.ndimage.grey_erosion(eog_signal, size=int((0.2 / 2) * sampling_rate))

    # a "low-noise" Lanczos differentiation filter introduced in Hamming (1989) is employed.
    # Frequently, a first order differentiation filter is sufficient and has the familiar
    # form of symmetric difference:
    # w[k] = 0.5 * (y[k + 1] - y[k - 1])
    diff = eroded - np.concatenate([[0], 0.5 * np.diff(eroded)])

    # To reduce the effects of noise, characterized by small fluctuations around zero, a
    # median filter is also used with the order of the median filter denoted as M.
    # The median filter acts like a mean filter except that it preserves the sharp edges ofthe
    # input. The order M should be less than a quarter ofthe expected eyeblink duration.
    clean = scipy.ndimage.median_filter(diff, size=int((0.2 / 4) * sampling_rate))

    return clean 

def _hrv_get_rri(peaks=None, sampling_rate=1000, interpolate=False, **kwargs):

    rri = np.diff(peaks) / sampling_rate * 1000

    if interpolate is False:
        return rri

    else:

        # Minimum sampling rate for interpolation
        if sampling_rate < 10:
            sampling_rate = 10

        # Compute length of interpolated heart period signal at requested sampling rate.
        desired_length = int(np.rint(peaks[-1] / sampling_rate * sampling_rate))

        rri = signal_interpolate(
            peaks[1:],  # Skip first peak since it has no corresponding element in heart_period
            rri,
            x_new=np.arange(desired_length),
            **kwargs
        )
        return rri, sampling_rate 

def test_basic(self):
        x = [1, 4, 6, 7, 12]
        out = np.array([3, 2, 1, 5])
        out2 = np.array([-1, -1, 4])
        out3 = np.array([0, 5])
        assert_array_equal(diff(x), out)
        assert_array_equal(diff(x, n=2), out2)
        assert_array_equal(diff(x, n=3), out3)

        x = [1.1, 2.2, 3.0, -0.2, -0.1]
        out = np.array([1.1, 0.8, -3.2, 0.1])
        assert_almost_equal(diff(x), out)

        x = [True, True, False, False]
        out = np.array([False, True, False])
        out2 = np.array([True, True])
        assert_array_equal(diff(x), out)
        assert_array_equal(diff(x, n=2), out2) 

def pgradient(self):
        d = {}
        for parm in self.parameters:
            s = int("beta" in parm)
            # Get AOs for our spin channel only
            i0, i1 = s * self._nelec[0], self._nelec[0] + s * self._nelec[1]
            ao = self._aovals[:, :, i0:i1, :]  # (kpt, config, electron, ao)
            pgrad_shape = (ao.shape[-3],) + self.parameters[parm].shape
            pgrad = np.zeros(pgrad_shape, dtype=self.dtype)  # (nconf, coeff)
            # Compute derivatives w.r.t. MO coefficients
            if ao.shape[0] > 1:  # multiple kpts
                split_sizes = np.diff([0] + list(self.param_split[parm]))
                k = np.repeat(np.arange(self.nk), split_sizes)
                for i in range(self._nelec[s]):  # MO loop
                    pgrad[:, :, i] = self._testcol(i, s, ao[k[i]])
            else:
                ao = ao[0]
                for i in range(self._nelec[s]):  # MO loop
                    pgrad[:, :, i] = self._testcol(i, s, ao)
            d[parm] = np.asarray(pgrad)
        return d 

def _bounds_from_array(arr, dim_name, bounds_name):
    """Get the bounds of an array given its center values.

    E.g. if lat-lon grid center lat/lon values are known, but not the
    bounds of each grid box.  The algorithm assumes that the bounds
    are simply halfway between each pair of center values.
    """
    # TODO: don't assume needed dimension is in axis=0
    # TODO: refactor to get rid of repetitive code
    spacing = arr.diff(dim_name).values
    lower = xr.DataArray(np.empty_like(arr), dims=arr.dims,
                         coords=arr.coords)
    lower.values[:-1] = arr.values[:-1] - 0.5*spacing
    lower.values[-1] = arr.values[-1] - 0.5*spacing[-1]
    upper = xr.DataArray(np.empty_like(arr), dims=arr.dims,
                         coords=arr.coords)
    upper.values[:-1] = arr.values[:-1] + 0.5*spacing
    upper.values[-1] = arr.values[-1] + 0.5*spacing[-1]
    bounds = xr.concat([lower, upper], dim='bounds')
    return bounds.T 

def test_floor_exact_64():
    # float64
    for e in range(53, 63):
        start = np.float64(2**e)
        across = start + np.arange(2048, dtype=np.float64)
        gaps = set(np.diff(across)).difference([0])
        assert_equal(len(gaps), 1)
        gap = gaps.pop()
        assert_equal(gap, int(gap))
        test_val = 2**(e+1)-1
        assert_equal(floor_exact(test_val, np.float64), 2**(e+1)-int(gap)) 

def unwrap(angle, period=2*math.pi):
    """Unwrap a phase angle to give a continuous curve

    Parameters
    ----------
    angle : array_like
        Array of angles to be unwrapped
    period : float, optional
        Period (defaults to `2*pi`)

    Returns
    -------
    angle_out : array_like
        Output array, with jumps of period/2 eliminated

    Examples
    --------
    >>> import numpy as np
    >>> theta = [5.74, 5.97, 6.19, 0.13, 0.35, 0.57]
    >>> unwrap(theta, period=2 * np.pi)
    [5.74, 5.97, 6.19, 6.413185307179586, 6.633185307179586, 6.8531853071795865]

    """
    dangle = np.diff(angle)
    dangle_desired = (dangle + period/2.) % period - period/2.
    correction = np.cumsum(dangle_desired - dangle)
    angle[1:] += correction
    return angle 

def _pandas_to_bucket(df, symbol, initial_image):
        rtn = {SYMBOL: symbol, VERSION: CHUNK_VERSION_NUMBER, COLUMNS: {}, COUNT: len(df)}
        end = to_dt(df.index[-1].to_pydatetime())
        if initial_image:
            if 'index' in initial_image:
                start = min(to_dt(df.index[0].to_pydatetime()), initial_image['index'])
            else:
                start = to_dt(df.index[0].to_pydatetime())
            image_start = initial_image.get('index', start)
            rtn[IMAGE_DOC] = {IMAGE_TIME: image_start, IMAGE: initial_image}
            final_image = TickStore._pandas_compute_final_image(df, initial_image, end)
        else:
            start = to_dt(df.index[0].to_pydatetime())
            final_image = {}
        rtn[END] = end
        rtn[START] = start

        logger.warning("NB treating all values as 'exists' - no longer sparse")
        rowmask = Binary(lz4_compressHC(np.packbits(np.ones(len(df), dtype='uint8')).tostring()))

        index_name = df.index.names[0] or "index"
        if PD_VER < '0.23.0':
            recs = df.to_records(convert_datetime64=False)
        else:
            recs = df.to_records()

        for col in df:
            array = TickStore._ensure_supported_dtypes(recs[col])
            col_data = {
                DATA: Binary(lz4_compressHC(array.tostring())),
                ROWMASK: rowmask,
                DTYPE: TickStore._str_dtype(array.dtype),
            }
            rtn[COLUMNS][col] = col_data
        rtn[INDEX] = Binary(
            lz4_compressHC(np.concatenate(
                ([recs[index_name][0].astype('datetime64[ms]').view('uint64')],
                 np.diff(
                     recs[index_name].astype('datetime64[ms]').view('uint64')))).tostring()))
        return rtn, final_image