def radixSortKeysOnly(self, keyBits):
     i = numpy.uint32(0)
     bitStep = self.bitStep
     
     while (keyBits > i*bitStep):
         self.radixSortStepKeysOnly(bitStep, i*bitStep)
         i+=numpy.uint32(1)

    def testIntMax(self):
        num = np.int(np.iinfo(np.int).max)
        self.assertEqual(np.int(ujson.decode(ujson.encode(num))), num)

        num = np.int8(np.iinfo(np.int8).max)
        self.assertEqual(np.int8(ujson.decode(ujson.encode(num))), num)

        num = np.int16(np.iinfo(np.int16).max)
        self.assertEqual(np.int16(ujson.decode(ujson.encode(num))), num)

        num = np.int32(np.iinfo(np.int32).max)
        self.assertEqual(np.int32(ujson.decode(ujson.encode(num))), num)

        num = np.uint8(np.iinfo(np.uint8).max)
        self.assertEqual(np.uint8(ujson.decode(ujson.encode(num))), num)

        num = np.uint16(np.iinfo(np.uint16).max)
        self.assertEqual(np.uint16(ujson.decode(ujson.encode(num))), num)

        num = np.uint32(np.iinfo(np.uint32).max)
        self.assertEqual(np.uint32(ujson.decode(ujson.encode(num))), num)

        if platform.architecture()[0] != '32bit':
            num = np.int64(np.iinfo(np.int64).max)
            self.assertEqual(np.int64(ujson.decode(ujson.encode(num))), num)

            # uint64 max will always overflow as it's encoded to signed
            num = np.uint64(np.iinfo(np.int64).max)
            self.assertEqual(np.uint64(ujson.decode(ujson.encode(num))), num)

def combine_scaled(r, g, b, a):
    """Combine components in [0, 1] to rgba uint32"""
    r2 = min(255, np.uint32(r * 255))
    g2 = min(255, np.uint32(g * 255))
    b2 = min(255, np.uint32(b * 255))
    a2 = min(255, np.uint32(a * 255))
    return np.uint32((a2 << 24) | (b2 << 16) | (g2 << 8) | r2)

def get_bucket(uhash, pieces, first_bucket_vector, second_bucket_vector):
    h_index = np.uint64(first_bucket_vector[0]) + np.uint64(second_bucket_vector[0 + 2])
    if h_index >= const.prime_default:
        h_index -= const.prime_default
    assert(h_index < const.prime_default)
    h_index = np.uint32(h_index)
    h_index = h_index % uhash.table_size

    control = np.uint64(first_bucket_vector[1]) + np.uint64(second_bucket_vector[1 + 2])
    if control >= const.prime_default:
        control -= const.prime_default
    assert(control < const.prime_default)
    control = np.uint32(control)
    
    if uhash.t == 2:
        index_hybrid = uhash.hybrid_hash_table[h_index]

        while index_hybrid:
            if index_hybrid.control_value == control:
                index_hybrid = C.pointer(index_hybrid)[1]
                return index_hybrid
            else:
                index_hybrid = C.pointer(index_hybrid)[1]
                if index_hybrid.point.is_last_bucket:
                    return None
                l = index_hybrid.point.bucket_length
                index_hybrid = C.pointer(index_hybrid)[l]
        return None

    def pad(self, max_pad):
        """
        Pad the timeseries data so that there are no missing values. We fill in
        missing power values using the previous power value in the series.
        """
        width = self.times[-1] - self.times[0] + 1
        padded_array = np.rec.array((0, 2), dtype=[('time', np.uint32),
                                                   ('power', np.float32)])
        padded_array.resize(width)

        cnt = 0
        for i in xrange(len(self.times)-1):
            padded_array[cnt] = (np.uint32(self.times[i]), self.powers[i])
            cnt += 1

            if self.times[i+1] - self.times[i] > max_pad:
                continue

            for t in xrange(self.times[i]+1, self.times[i+1]):
                padded_array[cnt] = (np.uint32(t), self.powers[i])
                cnt += 1

        padded_array[cnt] = (np.uint32(self.times[-1]), self.powers[-1])
        padded_array.resize(cnt + 1)

        self.array = padded_array

    def __call__(self, queue, tgt, src, shape):
        w, h = shape
        assert w % block_size == 0
        assert h % block_size == 0

        return self.kernel(queue, (w, h), None,
            tgt, src, numpy.uint32(w), numpy.uint32(h))

	def decode(self, server, block_header, target, job_id = None, extranonce2 = None):
		if block_header:
			job = Object()
	
			binary_data = block_header.decode('hex')
			data0 = np.zeros(64, np.uint32)
			data0 = np.insert(data0, [0] * 16, unpack('IIIIIIIIIIIIIIII', binary_data[:64]))
	
			job.target	  = np.array(unpack('IIIIIIII', target.decode('hex')), dtype=np.uint32)
			job.header	  = binary_data[:68]
			job.merkle_end  = np.uint32(unpack('I', binary_data[64:68])[0])
			job.time		= np.uint32(unpack('I', binary_data[68:72])[0])
			job.difficulty  = np.uint32(unpack('I', binary_data[72:76])[0])
			job.state	   = sha256(STATE, data0)
			job.f		   = np.zeros(8, np.uint32)
			job.state2	  = partial(job.state, job.merkle_end, job.time, job.difficulty, job.f)
			job.targetQ	 = 2**256 / int(''.join(list(chunks(target, 2))[::-1]), 16)
			job.job_id	  = job_id
			job.extranonce2 = extranonce2
			job.server	  = server
	
			calculateF(job.state, job.merkle_end, job.time, job.difficulty, job.f, job.state2)

			if job.difficulty != self.difficulty:
				self.set_difficulty(job.difficulty)
	
			return job

def color_count(image):
    '''Considering a (w,h,3) image of (dtype=uint8),
       compute the number of unique colors
    
       Encoding (i,j,k) into single index N = i+R*j+R*C*k
       Decoding N into (i,j,k) = (N-k*R*C-j*R, (N-k*R*C)/R, N/(R*C))
       using integer division\n

        Inputs:  image
        Returns: NumPy array of unique colors,
                 number of pixels of each unique color in image
    '''
    #Need to convert image to uint32 before multiplication so numbers are not truncated
    F = np.uint32(image[...,0])*256*256 + np.uint32(image[...,1])*256 + np.uint32(image[...,2])
    unique, counts = np.unique(F, return_counts=True)
    colors = np.empty(shape=(len(unique),3),dtype=np.uint32)
    numcol = np.empty(len(unique),dtype=np.uint32)
    i = 0
    for col,num in zip(unique,counts):
        R = col/(256*256)
        G = (col-R*256*256)/256
        B = (col-R*256*256-G*256)
        colors[i] = (R,G,B)
        numcol[i] = num
        i+=1
    return colors, numcol

	def calculate_pressure(self):
		# Calculate atmospheric pressure in [Pa]
		self.B6 = self.B5 - 4000
		print("B6=",self.B6)
		X1 = (self.B2 * (self.B6 * self.B6 / 2**12)) / 2**11
		print("X1=",X1)
		X2 = self.AC2 * self.B6 / 2**11
		print("X2=",X2)
		X3 = X1 + X2
		print("X3=",X3)
		self.B3 = (((self.AC1 * 4 + X3) << self.BMP183_CMD['OVERSAMPLE_3']) + 2 ) / 4
		print("B3=",self.B3)
		X1 = self.AC3 * self.B6 / 2**13
		print("X1=",X1)
		X2 = (self.B1 * (self.B6 * self.B6 / 2**12)) / 2**16
		print("X2=",X2)
		X3 = ((X1 + X2) + 2) / 2**2
		print("X3=",X3)
		self.B4 = numpy.uint32 (self.AC4 * (X3 + 32768) / 2**15)
		print("B4=",self.B4)
		self.B7 = (numpy.uint32 (self.UP) - self.B3) * (50000 >> self.BMP183_CMD['OVERSAMPLE_3'])
		print("B7=",self.B7)
		p = numpy.uint32 ((self.B7 * 2) / self.B4)
		print("p=",p)
		X1 = (p / 2**8) * ( p / 2**8)
		print("X1=",X1)
		X1 = int (X1 * 3038) / 2**16
		print("X1=",X1)
		X2 = int (-7357 * p) / 2**16
		print("X2=",X2)
		self.pressure = p + (X1 + X2 +3791) / 2**4
		print("pressure=",self.pressure)

def communion_encode(msg):
    assert msg["mode"] in ("request", "response")
    m = 'SEAMLESS'.encode()
    tip = b'\x00' if msg["mode"] == "request" else b'\x01'
    m += tip

    m += np.uint32(msg["id"]).tobytes()
    remainder = msg.copy()
    remainder.pop("mode")
    remainder.pop("id")
    remainder.pop("content")    
    if len(remainder.keys()):
        rem = json.dumps(remainder).encode()
        nrem = np.uint32(len(rem))
        m += nrem
        m += rem
    else:
        m += b'\x00\x00\x00\x00'
    content = msg["content"]
    if content is None:
        m += b'\x00'
    else:
        assert isinstance(content, (str, bytes, bool)), content
        if isinstance(content, bool):
            is_str = b'\x01'
        else:
            is_str = b'\x03' if isinstance(content, str) else b'\x02'
        m += is_str
        if isinstance(content, str):
            content = content.encode()
        elif isinstance(content, bool):
            content = b'\x01' if content else b'\x00'
        m += content
    #assert communion_decode(m) == msg, (communion_decode(m), msg)
    return m

 def set_demod(self, data):
     @command()
     def set_demod_buffer(self, data): 
         pass
     data1 = np.uint32(np.mod(np.floor(8192 * data[0, :]) + 8192,16384) + 8192)
     data2 = np.uint32(np.mod(np.floor(8192 * data[1, :]) + 8192,16384) + 8192)
     set_demod_buffer(self, data1 + data2 * 2**16)

 def set_dac(self):
     @command()
     def set_dac_data(self, data):
         pass
     dac_data_1 = np.uint32(np.mod(np.floor(8192 * self.dac[0, :]) + 8192, 16384) + 8192)
     dac_data_2 = np.uint32(np.mod(np.floor(8192 * self.dac[1, :]) + 8192, 16384) + 8192)
     set_dac_data(self, dac_data_1 + 65536 * dac_data_2)

 def createSplitPerm(self, size, subset_ratio=0.8, seed=None):
     rng = np.random.RandomState(np.uint32((time.time())))
     if seed is not None:
         rng.seed(np.uint32(seed))
     perm = rng.permutation(size)
     k = int(size * subset_ratio)
     return perm[:k]

 def index_to_array(self, index, dim, result = None):
     '''
     index_to_array(float* index, float* result, uint r_col, uint r_row)
     '''
     if result is None:
         result = gpuarray.GPUArray([index.size, dim], dtype = self.dtype)
           
         self.index_to_array_kernel(index.gpudata, result.gpudata, \
                                    np.uint32(r_col), np.uint32(r_row), \
                                    block = self._2d_block, \
                                    grid = self._2d_grid(r_col, r_row) \
                                    )     
         return result   
     
     else:
         r_col = result.shape[0]
         r_row = result.shape[1]
         
         if r_col != index.size or r_row != dim:
             raise ValueError('index_to_array: the dim of the result dont match the input')
           
         self.index_to_array_kernel(index.gpudata, result.gpudata, \
                                    np.uint32(r_col), np.uint32(r_row), \
                                    block = self._2d_block, \
                                    grid = self._2d_grid(r_col, r_row) \
                                    )      

 def getTimeBaseIndices(self, name, tBegin, tEnd):
     """ Return time indices of name corresponding to tBegin and tEnd """
     if not self.status:
         raise Exception('Shotfile not open!')
     try:
         sigName = ctypes.c_char_p(name)
     except TypeError:
         sigName = ctypes.c_char_p(name.encode())
     error = ctypes.c_int32(0)
     info = self.getTimeBaseInfo(name)
     if tEnd < tBegin:
         temp = tEnd
         tEnd = tBegin
         tBegin = temp
     if tBegin < info.tBegin:
         tBegin = info.tBegin
     if tEnd > info.tEnd:
         tEnd = info.tEnd
     try:
         time1 = ctypes.c_float(tBegin)
     except TypeError:
         time1 = ctypes.c_float(tBegin.value)
     try:
         time2 = ctypes.c_float(tEnd)
     except TypeError:
         time2 = ctypes.c_float(tEnd.value)
     k1 = ctypes.c_uint32(0)
     k2 = ctypes.c_uint32(0)
     lname = ctypes.c_uint64(len(name))
     __libddww__.ddtindex_(ctypes.byref(error), ctypes.byref(self.diaref), sigName, ctypes.byref(time1), 
                           ctypes.byref(time2), ctypes.byref(k1), ctypes.byref(k2), lname)
     getError(error.value)
     return numpy.uint32(k1.value), numpy.uint32(k2.value)

    def set_data(self, time, lcids=None, pbids=None, nsamples=None, exptimes=None):
        mf = cl.mem_flags

        if self._b_time is not None:
            self._b_time.release()
            self._b_lcids.release()
            self._b_pbids.release()
            self._b_nsamples.release()
            self._b_etimes.release()

        self.nlc = uint32(1 if lcids is None else unique(lcids).size)
        self.npb = uint32(1 if pbids is None else unique(pbids).size)
        self.nptb = time.size

        self.time = asarray(time, dtype='float32')
        self.lcids = zeros(time.size, 'uint32') if lcids is None else asarray(lcids, dtype='uint32')
        self.pbids = zeros(self.nlc, 'uint32') if pbids is None else asarray(pbids, dtype='uint32')
        self.nsamples = ones(self.nlc, 'uint32') if nsamples is None else asarray(nsamples, dtype='uint32')
        self.exptimes = ones(self.nlc, 'float32') if exptimes is None else asarray(exptimes, dtype='float32')

        self._b_time = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.time)
        self._b_lcids = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.lcids)
        self._b_pbids = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.pbids)
        self._b_nsamples = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.nsamples)
        self._b_etimes = cl.Buffer(self.ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=self.exptimes)

 def __extractRegionIndices(self):
     ''' Function generates, and returns, a dictionary of indices for regions in
     self.regionFile.'''
     # Create variable to store region indices and loop through file
     regionIndices = {}
     with open(self.regionFile, 'r') as infile:
         for line in infile:
             # Extract region data and find indices of matching bins
             chrom, start, end, region = line.strip().split('\t')
             acceptableBins = (
                 (self.binChr == chrom)
                 & (self.binStart >= np.uint32(start))
                 & (self.binEnd <= np.uint32(end)))
             indices = np.where(acceptableBins)[0]
             # Add region indices to dictionary
             if region in regionIndices:
                 regionIndices[region] = np.concatenate(
                     (regionIndices[region], indices))
             else:
                 regionIndices[region] = indices
     # Check region dictionary for absent or erronous regions
     for region in regionIndices:
         # Extract indices, sort, and update dictionary
         indices = regionIndices[region]
         indices.sort()
         regionIndices[region] = indices
         # Check for absent or duplicate indices
         if len(indices) == 0:
             raise IOError('{} has no bins'.format(region))
         if len(set(indices)) != len(indices):
             raise IOError('{} has overlapping segments'.format(region))
     # Return region index data
     return(regionIndices)

    def test_valid(self):
        prop = bcpp.Int()

        assert prop.is_valid(None)

        assert prop.is_valid(0)
        assert prop.is_valid(1)

        assert prop.is_valid(np.int8(0))
        assert prop.is_valid(np.int8(1))
        assert prop.is_valid(np.int16(0))
        assert prop.is_valid(np.int16(1))
        assert prop.is_valid(np.int32(0))
        assert prop.is_valid(np.int32(1))
        assert prop.is_valid(np.int64(0))
        assert prop.is_valid(np.int64(1))
        assert prop.is_valid(np.uint8(0))
        assert prop.is_valid(np.uint8(1))
        assert prop.is_valid(np.uint16(0))
        assert prop.is_valid(np.uint16(1))
        assert prop.is_valid(np.uint32(0))
        assert prop.is_valid(np.uint32(1))
        assert prop.is_valid(np.uint64(0))
        assert prop.is_valid(np.uint64(1))

        # TODO (bev) should fail
        assert prop.is_valid(False)
        assert prop.is_valid(True)

def _minmax_impl(a_gpu, axis, min_or_max, stream=None):
    ''' Returns both max and argmax (min/argmin) along an axis.'''
    assert len(a_gpu.shape) < 3
    if iscomplextype(a_gpu.dtype):
        raise ValueError("Cannot compute min/max of complex values")

    if axis is None:  ## Note: PyCUDA doesn't have an overall argmax/argmin!
        if min_or_max == 'max':
            return gpuarray.max(a_gpu).get()
        else:
            return gpuarray.min(a_gpu).get()
    else:
        if axis < 0:
            axis += 2
    assert axis in (0, 1)

    global _global_cublas_allocator
    alloc = _global_cublas_allocator

    n, m = a_gpu.shape if a_gpu.flags.c_contiguous else (a_gpu.shape[1], a_gpu.shape[0])
    col_kernel, row_kernel = _get_minmax_kernel(a_gpu.dtype, min_or_max)
    if (axis == 0 and a_gpu.flags.c_contiguous) or (axis == 1 and a_gpu.flags.f_contiguous):
        target = gpuarray.empty(m, dtype=a_gpu.dtype, allocator=alloc)
        idx = gpuarray.empty(m, dtype=np.uint32, allocator=alloc)
        col_kernel(a_gpu, target, idx, np.uint32(m), np.uint32(n),
                   block=(32, 1, 1), grid=(m, 1, 1), stream=stream)
    else:
        target = gpuarray.empty(n, dtype=a_gpu, allocator=alloc)
        idx = gpuarray.empty(n, dtype=np.uint32, allocator=alloc)
        row_kernel(a_gpu, target, idx, np.uint32(m), np.uint32(n),
                block=(32, 1, 1), grid=(n, 1, 1), stream=stream)
    return target, idx

    def arm_timed_latch(self, latch_name, time=None, force=False):
        ''' Arm a timed latch. Use force=True to force even if already armed 
        @param fpga host
        @param latch_name: base name for latch
        @param time: time in adc samples since epoch
        @param force: force arm of latch that is already armed
        '''
        status = self.host.device_by_name('%s_status' %(latch_name)).read()['data']

        # get armed, arm and load counts
        armed_before = status['armed']
        arm_count_before = status['arm_count']
        load_count_before = status['load_count']

        # if not forcing it, check for already armed first
        if armed_before == True:
            if force == False:
                LOGGER.info('forcing arm of already armed timed latch %s' %latch_name)
            else:
                LOGGER.error('timed latch %s already armed, use force=True' %latch_name)
                return

        # we load immediate if not given time
        if time == None:
            self.host.device_by_name('%s_control0' %(latch_name)).write(arm=0, load_immediate=1)
            self.host.device_by_name('%s_control0' %(latch_name)).write(arm=1, load_immediate=1)
                
            LOGGER.info('Timed latch %s arm-for-immediate-loading attempt' %latch_name)
            
        else: 
            # TODO check time values
            time_samples = numpy.uint64(time)
            time_msw = numpy.uint32((time_samples & 0x0000FFFFFFFFFFFF) >> 32)
            time_lsw = numpy.uint32((time_samples & 0x00000000FFFFFFFF))
            self.host.device_by_name('%s_control' %(latch_name)).write(load_time_lsw=time_lsw)
            
            self.host.device_by_name('%s_control0' %(latch_name)).write(arm=0, load_immediate=0, load_time_msw=time_msw)
            self.host.device_by_name('%s_control0' %(latch_name)).write(arm=1, load_immediate=0, load_time_msw=time_msw)

            LOGGER.info('Timed latch %s arm-for-loading-at-time attempt' %latch_name)

        # TODO check that arm count increased as expected
        status = self.host.device_by_name('%s_status' %(latch_name)).read()['data']
        
        # get armed, arm and load counts
        armed_after = status['armed']
        arm_count_after = status['arm_count']
        load_count_after = status['load_count']
        
        # armed count did not succeed
        if arm_count_after != (arm_count_before+1):
            # TODO check time
            LOGGER.error('Timed latch %s arm count at %i instead of %i' %(latch_name, arm_count_after, (arm_count_before+1)))
        else:
            # check load count increased as expected
            if time == None: 
                if load_count_after != (load_count_before+1):
                    LOGGER.error('Timed latch %s load count at %i instead of %i' %(latch_name, load_count_after, (load_count_before+1)))
            else:
                LOGGER.info('Timed latch %s successfully armed' %(latch_name))