def _get_init_guess(strsa, strsb, nroots, hdiag, orbsym, wfnsym=0):
    airreps = numpy.zeros(strsa.size, dtype=numpy.int32)
    birreps = numpy.zeros(strsb.size, dtype=numpy.int32)
    for i, ir in enumerate(orbsym):
        airreps[numpy.bitwise_and(strsa, 1<<i) > 0] ^= ir
        birreps[numpy.bitwise_and(strsb, 1<<i) > 0] ^= ir
    na = len(strsa)
    nb = len(strsb)

    ci0 = []
    iroot = 0
    for addr in numpy.argsort(hdiag):
        x = numpy.zeros((na*nb))
        addra = addr // nb
        addrb = addr % nb
        if airreps[addra] ^ birreps[addrb] == wfnsym:
            x[addr] = 1
            ci0.append(x)
            iroot += 1
            if iroot >= nroots:
                break
    try:
        # Add noise
        ci0[0][0 ] += 1e-5
        ci0[0][-1] -= 1e-5
    except IndexError:
        raise IndexError('Configuration of required symmetry (wfnsym=%d) not found' % wfnsym)
    return ci0

    def test(self, X, y, verbose=True):
        # if we don't need 3d inputs...
        if type(self.model.input_shape) is tuple:
            X = np.array(X)
            if len(self.model.input_shape) == 2:
                X = X.reshape((X.shape[0], -1))
        else:
            raise LanguageClassifierException('Mult-input models are not supported yet')

        if verbose:
            print("Getting predictions on the test set")
        predictions = self.predict(X)

        if len(predictions) != len(y):
            raise LanguageClassifierException("Non comparable arrays")

        if self.binary:
            acc = (predictions == y).mean()
            prec = np.sum(np.bitwise_and(predictions, y)) * 1.0 / np.sum(predictions)
            recall = np.sum(np.bitwise_and(predictions, y)) * 1.0 / np.sum(y)
            if verbose:
                print("Test set accuracy of {0:.3f}%".format(acc * 100.0))
                print("Test set error of {0:.3f}%".format((1 - acc) * 100.0))
                print("Precision for class=1: {0:.3f}".format(prec))
                print("Recall for class=1: {0:.3f}".format(recall))

            return (acc, prec, recall)
        else:
            # TODO: Obtain more metrics for the multiclass problem
            acc = (predictions == y).mean()
            if verbose:
                print("Test set accuracy of {0:.3f}%".format(acc * 100.0))
                print("Test set error of {0:.3f}%".format((1 - acc) * 100.0))
            return acc

    def _read_symbology_block(self, buf2):
        """ Read symbology block. """
        # Read and decode symbology header
        self.symbology_header = _unpack_from_buf(buf2, 0, SYMBOLOGY_HEADER)

        # Read radial packets
        packet_code = struct.unpack('>h', buf2[16:18])[0]
        assert packet_code in SUPPORTED_PACKET_CODES
        self.packet_header = _unpack_from_buf(buf2, 16, RADIAL_PACKET_HEADER)
        self.radial_headers = []
        nbins = self.packet_header['nbins']
        nradials = self.packet_header['nradials']
        nbytes = _unpack_from_buf(buf2, 30, RADIAL_HEADER)['nbytes']
        if packet_code == 16 and nbytes != nbins:
            nbins = nbytes  # sometimes these do not match, use nbytes
        self.raw_data = np.empty((nradials, nbins), dtype='uint8')
        pos = 30

        for radial in self.raw_data:
            radial_header = _unpack_from_buf(buf2, pos, RADIAL_HEADER)
            pos += 6
            if packet_code == 16:
                radial[:] = np.fromstring(buf2[pos:pos+nbins], '>u1')
                pos += radial_header['nbytes']
            else:
                assert packet_code == AF1F
                # decode run length encoding
                rle_size = radial_header['nbytes'] * 2
                rle = np.fromstring(buf2[pos:pos+rle_size], dtype='>u1')
                colors = np.bitwise_and(rle, 0b00001111)
                runs = np.bitwise_and(rle, 0b11110000) // 16
                radial[:] = np.repeat(colors, runs)
                pos += rle_size
            self.radial_headers.append(radial_header)

def moments(data):
    """Returns (height, x, y, width_x, width_y,offset)
    the gaussian parameters of a 2D distribution by calculating its
    moments """
    total = data.sum()
    X, Y = np.indices(data.shape)
    x = (X*data).sum()/total
    y = (Y*data).sum()/total
    height = data.max()
    firstq = np.median(data[data < np.median(data)])
    thirdq = np.median(data[data > np.median(data)])
    offset = np.median(data[np.where(np.bitwise_and(data > firstq,
                                                    data < thirdq))])
    places = np.where((data-offset) > 4*np.std(data[np.where(np.bitwise_and(
                                      data > firstq, data < thirdq))]))
    width_y = np.std(places[0])
    width_x = np.std(places[1])
    # These if statements take into account there might only be one significant
    # point above the background when that is the case it is assumend the width
    # of the gaussian must be smaller than one pixel
    if width_y == 0.0:
        width_y = 0.5
    if width_x == 0.0:
        width_x = 0.5

    height -= offset
    return height, x, y, width_x, width_y, offset

def process_cloudmask(mod09a1_file_name, cloudmask_output_name):
    fn_mod09a1 = mod09a1_file_name

    stateflags, geoTransform, proj = return_band(fn_mod09a1, 11)  # band 11 -- 500m State Flags

    goodpix_mask = numpy.where(stateflags == 65535, 1, 0)

    cloud = numpy.bitwise_and(stateflags, 3) + 1
    numpy.putmask(cloud, goodpix_mask, 0)

    # print cloudmask

    not_set = numpy.where(cloud == 4, 1, 0)
    shadow1 = numpy.where(cloud == 1, 1, 0)
    shadow2 = numpy.where(numpy.bitwise_and(stateflags, 4) == 4, 1, 0)
    shadow = numpy.logical_and(shadow1, shadow2)
    numpy.putmask(cloud, shadow, 4)
    numpy.putmask(cloud, not_set, 1)

    blue = return_band(fn_mod09a1, 3)[0]
    too_blue1 = numpy.where(cloud == 1, 1, 0)
    too_blue2 = numpy.where(blue > 2000, 1, 0)
    too_blue = numpy.logical_and(too_blue1, too_blue2)
    numpy.putmask(cloud, too_blue, 5)

    # print cloud

    output_file(cloudmask_output_name, cloud, geoTransform, proj)

    stateflags = None
    cloud = None

 def place_ship(self, position, length, orientation):
     """
     Return None if ship cannot be placed
     """
     ship = None
     if orientation == 'H':
         zeros = np.zeros(self.width * self.height, dtype='int8')
         if (position[0] + length) > self.width:
             return None
         for i in range(length):
             zeros[position[1] * self.width + position[0]+i] = 1
         if np.all(np.bitwise_and(self._layout, zeros) == 0):
             self._layout = np.bitwise_or(self._layout, zeros)
             ship = Ship(position, length, orientation)
     elif orientation == 'V':
         zeros = np.zeros(self.width * self.height, dtype='int8')
         if (position[1] + length) > self.height:
             return None
         for i in range(length):
             zeros[(position[1] + i) * self.width + position[0]] = 1
         if np.all(np.bitwise_and(self._layout, zeros) == 0):
             self._layout = np.bitwise_or(self._layout, zeros)
             ship = Ship(position, length, orientation)
     if ship:
         self._ships.append(ship)
         return ship

    def rearrange_bits(array):
        # Do bit rearrangement for the 10-bit lytro raw format
        # Normalize output to 1.0 as float64
        t0 = array[0::5]
        t1 = array[1::5]
        t2 = array[2::5]
        t3 = array[3::5]
        lsb = array[4::5]

        t0 = np.left_shift(t0, 2) + np.bitwise_and(lsb, 3)
        t1 = np.left_shift(t1, 2) + np.right_shift(np.bitwise_and(lsb, 12), 2)
        t2 = np.left_shift(t2, 2) + np.right_shift(np.bitwise_and(lsb, 48), 4)
        t3 = np.left_shift(t3, 2) + np.right_shift(np.bitwise_and(lsb, 192), 6)

        image = np.zeros(LYTRO_ILLUM_IMAGE_SIZE, dtype=np.uint16)
        image[:, 0::4] = t0.reshape(
            (LYTRO_ILLUM_IMAGE_SIZE[0], LYTRO_ILLUM_IMAGE_SIZE[1] // 4)
        )
        image[:, 1::4] = t1.reshape(
            (LYTRO_ILLUM_IMAGE_SIZE[0], LYTRO_ILLUM_IMAGE_SIZE[1] // 4)
        )
        image[:, 2::4] = t2.reshape(
            (LYTRO_ILLUM_IMAGE_SIZE[0], LYTRO_ILLUM_IMAGE_SIZE[1] // 4)
        )
        image[:, 3::4] = t3.reshape(
            (LYTRO_ILLUM_IMAGE_SIZE[0], LYTRO_ILLUM_IMAGE_SIZE[1] // 4)
        )

        # Normalize data to 1.0 as 64-bit float.
        # Division is by 1023 as the Lytro Illum saves 10-bit raw data.
        return np.divide(image, 1023.0).astype(np.float64)

def simpleAdd(exp0, exp1, badPixelMask):
    """Add two exposures, avoiding bad pixels
    """
    imArr0, maskArr0, varArr0 = exp0.getMaskedImage().getArrays()
    imArr1, maskArr1, varArr1 = exp1.getMaskedImage().getArrays()
    expRes = exp0.Factory(exp0, True)
    miRes = expRes.getMaskedImage()
    imArrRes, maskArrRes, varArrRes = miRes.getArrays()

    weightMap = afwImage.ImageF(exp0.getDimensions())
    weightArr = weightMap.getArray()

    good0 = np.bitwise_and(maskArr0, badPixelMask) == 0
    good1 = np.bitwise_and(maskArr1, badPixelMask) == 0

    imArrRes[:, :] = np.where(good0,  imArr0, 0) + np.where(good1,  imArr1, 0)
    varArrRes[:, :] = np.where(good0, varArr0, 0) + np.where(good1, varArr1, 0)
    maskArrRes[:, :] = np.bitwise_or(np.where(good0, maskArr0, 0), np.where(good1, maskArr1, 0))
    weightArr[:, :] = np.where(good0, 1, 0) + np.where(good1, 1, 0)

    miRes /= weightMap
    miRes *= 2  # want addition, not mean, where both pixels are valid

    setCoaddEdgeBits(miRes.getMask(), weightMap)

    return expRes

    def check_criterion(self, compiled_record, trial_record, **kwargs):
        trial_number = np.asarray(compiled_record['trial_number'])
        current_step = np.asarray(compiled_record['current_step'])
        correct = np.asarray(compiled_record['correct'])
        protocol_name = np.asarray(compiled_record['protocol_name'])
        protocol_ver = np.asarray(compiled_record['protocol_version_number'])

        # filter out trial_numbers for current protocol_name and protocol_ver
        current_step = current_step[np.bitwise_and(protocol_name==protocol_name[-1],protocol_ver==protocol_ver[-1])]
        trial_number = trial_number[np.bitwise_and(protocol_name==protocol_name[-1],protocol_ver==protocol_ver[-1])]
        correct = correct[np.bitwise_and(protocol_name==protocol_name[-1],protocol_ver==protocol_ver[-1])]

        if self.num_trials_mode == 'consecutive':
            jumps = np.where(np.diff(trial_number)!=1) # jumps in trial number
            if not jumps[0]:
                which_trials = trial_number
            else:
                which_trials = trial_number[jump[0][-1]:] # from the last jump
        else:
            which_trials = trial_number

        if np.size(which_trials)<self.num_trials:
            graduate = False # dont graduate if the number of trials less than num required
        else:
            which_trials = which_trials[-self.num_trials:]
            filter =  np.isin(trial_number,which_trials)
            correct = correct[filter]
            perf = np.sum(correct)/np.size(correct)
            if perf >self.pct_correct:
                graduate = True
            else:
                graduate = False

        return graduate

def computeState(isFix,md):
    ''' generic function that determines event start and end
        isFix - 1d array, time series with one element for each
            gaze data point, 1 indicates the event is on, 0 - off
        md - minimum event duration
        returns
            list with tuples with start and end for each
                event (values in frames)
            timeseries analogue to isFix but the values
                correspond to the list
    '''
    fixations=[]
    if isFix.sum()==0: return np.int32(isFix),[]
    fixon = np.bitwise_and(isFix,
        np.bitwise_not(np.roll(isFix,1))).nonzero()[0].tolist()
    fixoff=np.bitwise_and(np.roll(isFix,1),
        np.bitwise_not(isFix)).nonzero()[0].tolist()
    if len(fixon)==0 and len(fixoff)==0: fixon=[0]; fixoff=[isFix.size-1]
    if fixon[-1]>fixoff[-1]:fixoff.append(isFix.shape[0]-1)
    if fixon[0]>fixoff[0]:fixon.insert(0,0)
    if len(fixon)!=len(fixoff): print 'invalid fixonoff';raise TypeError
    for f in range(len(fixon)):
        fs=fixon[f];fe=(fixoff[f]+1);dur=fe-fs
        if  dur<md[0] or dur>md[1]:
            isFix[fs:fe]=False
        else: fixations.append([fs,fe-1])
    #fixations=np.array(fixations)
    return isFix,fixations

def interpolateBlinks(t,d,hz):
    ''' Interpolate short missing intervals
        d - 1d array, time series with gaze data, np.nan indicates blink
        hz - gaze data recording rate
    '''
    isblink= np.isnan(d)
    if isblink.sum()<2 or isblink.sum()>(isblink.size-2): return d
    blinkon = np.bitwise_and(isblink,np.bitwise_not(
        np.roll(isblink,1))).nonzero()[0].tolist()
    blinkoff=np.bitwise_and(np.roll(isblink,1),
        np.bitwise_not(isblink)).nonzero()[0].tolist()
    if len(blinkon)==0 and len(blinkoff)==0: return d
    #print 'bla',len(blinkon), len(blinkoff)
    if blinkon[-1]>blinkoff[-1]: blinkoff.append(t.size-1)
    if blinkon[0]>blinkoff[0]: blinkon.insert(0,0)
    if len(blinkon)!=len(blinkoff):
        print 'Blink Interpolation Failed'
        raise TypeError
    f=interp1d(t[~isblink],d[~isblink],bounds_error=False)
    for b in range(len(blinkon)):
        bs=blinkon[b]-1
        be=(blinkoff[b])
        if (be-bs)<INTERPMD*hz:
            d[bs:be]=f(t[bs:be])
            #for c in [7,8]: tser[bs:be,c]=np.nan
    return d

def mean_average_precision(distances, labels):
    """
    Calculate mean average precision and precision-recall breakeven.

    Returns
    -------
    mean_average_precision, mean_prb, ap_dict : float, float, dict
        The dict gives the per-type average precisions.
    """

    label_matches = generate_matches_array(labels)  # across all tokens

    ap_dict = {}
    prbs = []
    for target_type in sorted(set(labels)):
        if len(np.where(np.asarray(labels) == target_type)[0]) == 1:
            continue
        type_matches = generate_type_matches_array(labels, target_type)
        swtt_matches = np.bitwise_and(
            label_matches == True, type_matches == True
            ) # same word, target type
        dwtt_matches = np.bitwise_and(
            label_matches == False, type_matches == True
            ) # different word, target type
        ap, prb = average_precision(
            distances[swtt_matches], distances[dwtt_matches]
            )
        prbs.append(prb)
        ap_dict[target_type] = ap
    return np.mean(ap_dict.values()), np.mean(prbs), ap_dict

    def _record_data(self):
        '''
            Reads raw event data from SRAM and splits data stream into events
            ----------
            Returns:
                event_data : np.ndarray
                    Numpy array of single event numpy arrays 
        '''
        
        self.count_lost = self.dut['fadc0_rx'].get_count_lost()
#         print 'count_lost is %d' % self.count_lost
#         print 'event_count is %d' % self.event_count
        if self.count_lost > 0:
            logging.error('SRAM FIFO overflow number %d. Skip data.', self.count_lost)
            self.dut['fadc0_rx'].reset()
            self.set_adc_eventsize(self.sample_count, self.sample_delay)
            #return
        
        single_data = self.dut['DATA_FIFO'].get_data()                          # Read raw data from SRAM
        
        try:
            if single_data.shape[0] > 200:
                selection = np.where(single_data & 0x10000000 == 0x10000000)[0]         # Make mask from new-event-bit
                event_data = np.bitwise_and(single_data, 0x00003fff).astype(np.uint32)  # Remove new-event-bit from data       
                event_data = np.split(event_data, selection)                            # Split data into events by means of mask
                event_data = event_data[1:-1]                                           # Remove first and last event in case of chopping
                event_data = np.vstack(event_data)                                      # Stack events together
            else:
                event_data = np.asarray([np.bitwise_and(single_data, 0x00003fff).astype(np.uint32)])

            if event_data.shape[1] == self.sample_count:
                return event_data
        except ValueError as e:
            logging.error('_record_data() experienced a ValueError: ' + str(e))
            return

def check_fills_complement(figure1, figure2):
    figure1_bw = get_bw_image(figure1)
    figure2_bw = get_bw_image(figure2)
    #figure1_gray = figure1.convert('L')
    #figure1_bw = numpy.asarray(figure1_gray).copy()
    #figure2_gray = figure2.convert('L')
    #figure2_bw = numpy.asarray(figure2_gray).copy()
    #
    ## Dark area = 1, Light area = 0
    #figure1_bw[figure1_bw < 128] = 1
    #figure1_bw[figure1_bw >= 128] = 0
    #figure2_bw[figure2_bw < 128] = 1
    #figure2_bw[figure2_bw >= 128] = 0

    difference21 = figure2_bw - figure1_bw
    difference12 = figure1_bw - figure2_bw
    # 0 - 1 = 255 because of rollover.
    difference21[difference21 > 10] = 0
    difference12[difference12 > 10] = 0

    percent_diff = get_percent_diff(difference21, figure2_bw)
    if percent_diff < PERCENT_DIFF_THRESHOLD and \
        get_percent_diff(numpy.bitwise_and(figure1_bw, figure2_bw),
                         figure1_bw) < PERCENT_DIFF_THRESHOLD:
        return percent_diff
    percent_diff = get_percent_diff(difference12, figure1_bw)
    if percent_diff < PERCENT_DIFF_THRESHOLD and \
        get_percent_diff(numpy.bitwise_and(figure1_bw, figure2_bw),
                         figure2_bw) < PERCENT_DIFF_THRESHOLD:
        return -percent_diff

    return 0

def compute_dice_with_transfo(img_fixed, img_moving, transfo):
    
    #first transform 
    toolsPaths = ['CIP_PATH'];
    path=dict()
    for path_name in toolsPaths:
        path[path_name]=os.environ.get(path_name,False)
        if path[path_name] == False:
            print path_name + " environment variable is not set"
            exit()
    temp_out = "/Users/rolaharmouche/Documents/Data/temp_reg.nrrd"        
    resamplecall = os.path.join(path['CIP_PATH'], "ResampleCT")    
        
    sys_call = resamplecall+" -d "+img_fixed+" -r "+ temp_out+\
            " -t "+transfo+" -l "+img_moving
    os.system(sys_call) 
    
    print(" computing ssd between "+img_fixed+" and registered"+ img_moving)
                    
    img1_data, info = nrrd.read(temp_out)
    img2_data, info = nrrd.read(img_fixed)
    
    
    
    #careful reference image has labels =2 and 3
    added_images = img1_data
    np.bitwise_and(img1_data, img2_data, added_images)
    Dice_calculation = sum(added_images[:])*2.0/(sum(img1_data[:])+sum(img2_data[:]))

    return Dice_calculation

def process_t3records(t3records, time_bit=10, dtime_bit=15,
                      ch_bit=6, special_bit=True, ovcfunc=None):
    """Extract the different fields from the raw t3records array (.ht3).

    Returns:
        3 arrays representing detectors, timestamps and nanotimes.
    """
    if special_bit:
        ch_bit += 1
    assert ch_bit <= 8
    assert dtime_bit <= 16

    detectors = np.bitwise_and(
        np.right_shift(t3records, time_bit + dtime_bit), 2**ch_bit - 1).astype('uint8')
    nanotimes = np.bitwise_and(
        np.right_shift(t3records, time_bit), 2**dtime_bit - 1).astype('uint16')

    assert time_bit <= 16
    dt = np.dtype([('low16', 'uint16'), ('high16', 'uint16')])

    t3records_low16 = np.frombuffer(t3records, dt)['low16']     # View
    timestamps = t3records_low16.astype(np.int64)               # Copy
    np.bitwise_and(timestamps, 2**time_bit - 1, out=timestamps)

    overflow_ch = 2**ch_bit - 1
    overflow = 2**time_bit
    if ovcfunc is None:
        ovcfunc = _correct_overflow
    ovcfunc(timestamps, detectors, overflow_ch, overflow)
    return detectors, timestamps, nanotimes

 def convert_to_complex_samples(self, byte_arr):
     '''
     Generic parsing of complex samples.
     Handles 4-bit case
     `samples` is a byte array 
     TODO add others?
     Throws error if `self.bit_depth` is unsupported.
     '''
     # typically real is upper nibble, complex is lower nibble, but TODO make this generic
     if self.bit_depth == 4:
         lsb = (bitwise_and(byte_arr, 0x0f) << 4).astype(int8) >> 4
         msb = bitwise_and(byte_arr, 0xf0).astype(int8) >> 4
     elif self.bit_depth == 8:
         msb = byte_arr[0::2]  # TODO this might not be working
         lsb = byte_arr[1::2]
     elif self.bit_depth == 16:
         sample_arr = byte_arr.view(int16)
         msb = sample_arr[0::2]
         lsb = sample_arr[1::2]
     else:
         raise Error('Bit depth not supported for complex samples')
     if self.i_msb:
         return msb + 1j * lsb
     else:
         return lsb + 1j * msb

def center_mask_at(mask, pos, indexes, toroidal=False):
    ndim = len(mask)
    #print 'ndim:', ndim
    new = range(ndim)
    #print 'new:', new
    toKeep = np.ones(len(mask[0]), dtype=bool)
    #print 'toKeep:', toKeep
    for ic in xrange(ndim):
        #print 'ic:', ic
        #print 'pos[ic]:', pos[ic]
        comp = mask[ic] + pos[ic]
        #print 'comp:', comp
        if not toroidal:
            insiders = np.bitwise_and(comp>=0, comp<indexes.shape[ic])
            np.bitwise_and(toKeep, insiders, toKeep)
        else:
            comp[np.where(comp<0)] = indexes.shape[ic]-1
            comp[np.where(comp>=indexes.shape[ic])] = 0
            #print 'indexes.shape[ic]:', indexes.shape[ic]
            #print 'comp after modif:', comp
        new[ic] = comp
        #print 'new[ic]:', new[ic]
    #m = tuple( [n[toKeep] for n in new] )
    #np.bitwise_and(toKeep, indexes[m]!=-1, toKeep)
    #print 'toKeep:', toKeep
    #print 'new:', new
    return tuple( [n[toKeep] for n in new] )

    def checkDistribution(self,catalog):
        '''
        Separate CCD in 8 quadrants and check that are stars on all of them.

        :param catalog:
        :return:
        '''
        camera = self.getCam()
        width,heigth = camera.getPixelSize()
        ngrid = 8
        wgrid = np.linspace(0, width,ngrid)
        hgrid = np.linspace(0,heigth,ngrid)

        star_per_grid = np.zeros((ngrid-1,ngrid-1))

        for i in range(ngrid-1):
            mask_w = np.bitwise_and(catalog['X_IMAGE'] > wgrid[i],
                                    catalog['X_IMAGE'] < wgrid[i+1])
            for j in range(ngrid-1):
                mask_h = np.bitwise_and(catalog['Y_IMAGE'] > hgrid[i],
                                        catalog['Y_IMAGE'] < hgrid[i+1])
                mask = np.bitwise_and(mask_w,
                                      mask_h)

                star_per_grid[i][j] += np.sum(mask)

        nstar = len(catalog)/2/ngrid**2

        mask_starpg = star_per_grid < nstar

        if np.any(mask_starpg):
            raise StarDistributionException('Stellar distribution not suitable for optical alignment.')

def load_spc(fname):
    """Load data from Becker&Hickl SPC files.

    Returns:
        3 numpy arrays: timestamps, detector, nanotime
    """
    spc_dtype = np.dtype([("field0", "<u2"), ("b", "<u1"), ("c", "<u1"), ("a", "<u2")])
    data = np.fromfile(fname, dtype=spc_dtype)

    nanotime = 4095 - np.bitwise_and(data["field0"], 0x0FFF)
    detector = data["c"]

    # Build the macrotime (timestamps) using in-place operation for efficiency
    timestamps = data["b"].astype("int64")
    np.left_shift(timestamps, 16, out=timestamps)
    timestamps += data["a"]

    # extract the 13-th bit from data['field0']
    overflow = np.bitwise_and(np.right_shift(data["field0"], 13), 1)
    overflow = np.cumsum(overflow, dtype="int64")

    # Add the overflow bits
    timestamps += np.left_shift(overflow, 24)

    return timestamps, detector, nanotime