def test_init(self):
        import numpy as np
        import math
        import sys

        assert np.intp() == np.intp(0)
        assert np.intp("123") == np.intp(123)
        raises(TypeError, np.intp, None)
        assert np.float64() == np.float64(0)
        assert math.isnan(np.float64(None))
        assert np.bool_() == np.bool_(False)
        assert np.bool_("abc") == np.bool_(True)
        assert np.bool_(None) == np.bool_(False)
        assert np.complex_() == np.complex_(0)
        # raises(TypeError, np.complex_, '1+2j')
        assert math.isnan(np.complex_(None))
        for c in ["i", "I", "l", "L", "q", "Q"]:
            assert np.dtype(c).type().dtype.char == c
        for c in ["l", "q"]:
            assert np.dtype(c).type(sys.maxint) == sys.maxint
        for c in ["L", "Q"]:
            assert np.dtype(c).type(sys.maxint + 42) == sys.maxint + 42
        assert np.float32(np.array([True, False])).dtype == np.float32
        assert type(np.float32(np.array([True]))) is np.ndarray
        assert type(np.float32(1.0)) is np.float32
        a = np.array([True, False])
        assert np.bool_(a) is a

 def test_intp(self):
     # Ticket #99
     i_width = np.int_(0).nbytes*2 - 1
     np.intp('0x' + 'f'*i_width, 16)
     assert_raises(OverflowError, np.intp, '0x' + 'f'*(i_width+1), 16)
     assert_raises(ValueError, np.intp, '0x1', 32)
     assert_equal(255, np.intp('0xFF', 16))

def Connect3D(EToV):
    """Build global connectivity arrays for grid based on
    standard EToV input array from grid generator.
    """

    EToV = EToV.astype(np.intp)
    Nfaces = 4
    # Find number of elements and vertices
    K = EToV.shape[0]
    #Nv = EToV.max()+1

    # Create face to node connectivity matrix
    #TotalFaces = Nfaces*K

    # List of local face to local vertex connections
    vn = np.int32([[0,1,2],[0,1,3],[1,2,3],[0,2,3]])

    # Build global face to node connectivity
    g_face_no = 0
    vert_indices_to_face_numbers = {}
    face_numbers = xrange(Nfaces)
    for k in xrange(K):
        for face in face_numbers:
            vert_indices_to_face_numbers.setdefault(
                    frozenset(EToV[k,vn[face]]), []).append(g_face_no)
            g_face_no += 1

    faces1 = []
    faces2 = []

    # check this
    for i in vert_indices_to_face_numbers.itervalues():
        if len(i) == 2:
            faces1.append(i[0])
            faces2.append(i[1])
            faces2.append(i[0])
            faces1.append(i[1])

    faces1 = np.intp(faces1)
    faces2 = np.intp(faces2)

    # Convert faceglobal number to element and face numbers
    element1, face1 = divmod(faces1, Nfaces)
    element2, face2 = divmod(faces2, Nfaces)

    # Rearrange into Nelements x Nfaces sized arrays
    ind = element1*Nfaces + face1

    EToE = np.outer(np.arange(K), np.ones((1, Nfaces)))
    EToF = np.outer(np.ones((K,1)), np.arange(Nfaces))
    EToE = EToE.reshape(K*Nfaces)
    EToF = EToF.reshape(K*Nfaces)

    EToE[np.int32(ind)] = element2
    EToF[np.int32(ind)] = face2

    EToE = EToE.reshape(K, Nfaces)
    EToF = EToF.reshape(K, Nfaces)

    return  EToE, EToF

    def __init__(self, ldis, Nv, VX, VY, K, EToV):
        l = self.ldis = ldis

        self.dimensions = ldis.dimensions

        self.Nv = Nv
        self.VX   = VX
        self.K  = K

        va = np.intp(EToV[:, 0].T)
        vb = np.intp(EToV[:, 1].T)
        vc = np.intp(EToV[:, 2].T)

        x = self.x = 0.5*(
                -np.outer(VX[va], l.r+l.s, )
                +np.outer(VX[vb], 1+l.r)
                +np.outer(VX[vc], 1+l.s))
        y = self.y = 0.5*(
                -np.outer(VY[va], l.r+l.s)
                +np.outer(VY[vb], 1+l.r)
                +np.outer(VY[vc], 1+l.s))

        self.rx, self.sx, self.ry, self.sy, self.J = GeometricFactors2D(x, y, l.Dr, l.Ds)
        self.nx, self.ny, self.sJ = Normals2D(l, x, y, K)
        self.Fscale = self.sJ/self.J[:, l.FmaskF]

        # element-to-element, element-to-face connectivity
        self.EToE, self.EToF = Connect2D(EToV)

        self.mapM, self.mapP, self.vmapM, self.vmapP, self.vmapB, self.mapB = \
                BuildMaps2D(l, l.Fmask, VX, VY, EToV, self.EToE, self.EToF, K, l.N, x, y)

    def from_range(cat_comp, lo, hi):
        """
        Utility function to help construct the Roi from a range.

        :param cat_comp: Anything understood by ._categorical_helper ... array, list or component
        :param lo: lower bound of the range
        :param hi: upper bound of the range
        :return: CategoricalROI object
        """

        # Convert lo and hi to integers. Note that if lo or hi are negative,
        # which can happen if the user zoomed out, we need to reset the to zero
        # otherwise they will have strange effects when slicing the categories.

        # Note that we used ceil for lo, because if lo is 0.9 then we should
        # only select 1 and above.

        lo = np.intp(np.ceil(lo) if lo > 0 else 0)
        hi = np.intp(np.ceil(hi) if hi > 0 else 0)

        roi = CategoricalROI()
        cat_data = cat_comp.categories
        roi.update_categories(cat_data[lo:hi])

        return roi

 def test_init(self):
     import numpy as np
     import math
     import sys
     assert np.intp() == np.intp(0)
     assert np.intp('123') == np.intp(123)
     raises(TypeError, np.intp, None)
     assert np.float64() == np.float64(0)
     assert math.isnan(np.float64(None))
     assert np.bool_() == np.bool_(False)
     assert np.bool_('abc') == np.bool_(True)
     assert np.bool_(None) == np.bool_(False)
     assert np.complex_() == np.complex_(0)
     #raises(TypeError, np.complex_, '1+2j')
     assert math.isnan(np.complex_(None))
     for c in ['i', 'I', 'l', 'L', 'q', 'Q']:
         assert np.dtype(c).type().dtype.char == c
     for c in ['l', 'q']:
         assert np.dtype(c).type(sys.maxint) == sys.maxint
     for c in ['L', 'Q']:
         assert np.dtype(c).type(sys.maxint + 42) == sys.maxint + 42
     assert np.float32(np.array([True, False])).dtype == np.float32
     assert type(np.float32(np.array([True]))) is np.ndarray
     assert type(np.float32(1.0)) is np.float32
     a = np.array([True, False])
     assert np.bool_(a) is a

    def _accumDiffStims(self, d_resp_tmp, diffV1GausBuf, sizes, orderX,
                        orderY, orderT):
        """ Gets the responses of the filters specified in d_v1popDirs by
            interpolation.
            This is basically what shSwts.m did in the original S&H code."""

        # a useful list of factorials for computing the scaling factors for
        # the derivatives
        factorials = (1, 1, 2, 6)

        # the scaling factor for this directional derivative
        # similar to the binomial coefficients
        scale = 6/factorials[orderX]/factorials[orderY]/factorials[orderT]

        gdim = (int(iDivUp(sizes[0] * sizes[1], 256)), 1)
        bdim = (256, 1, 1)
        self.dev_accumDiffStims(
            np.intp(d_resp_tmp),
            np.intp(diffV1GausBuf),
            np.int32(sizes[0] * sizes[1]),
            np.int32(scale),
            np.int32(orderX),
            np.int32(orderY),
            np.int32(orderT),
            block=bdim, grid=gdim)

 def getrows(self, privateKey, senderID, mType, params, extra):
     """ 
     Receive a single row or list of rows and store in a class property variable 
     
     Use **client.crow** or **client.rowlist** to access the indices of 
     **previously broadcasted** rows 
     """
     filen = params['url']
     if mType == 'table.highlight.row':
         idx   = np.intp(params['row'])
         print '[SAMP] Selected row %s from %s' % (idx,filen)
         print '[SAMP] Row index stored in property -> crow'
         self.crow = idx
     elif mType == 'table.select.rowList':
         idx   = np.intp(params['row-list'])
         print '[SAMP] Selected %s rows from %s' % (len(idx),filen)
         print '[SAMP] List stored in property -> rowlist'
         self.rowlist = idx
     
     self.lastMessage = {'label':'Selected Rows',
                         'privateKey':privateKey, 
                         'senderID': senderID, 
                         'mType': mType, 
                         'params': params, 
                         'extra': extra }        

 def __init__(self, code, point, struct_ptr):
     self.code = cuda.to_device(code)
     self.point = cuda.to_device(point)
     self.code_shape, self.code_dtype = code.shape, code.dtype
     self.point_shape, self.point_dtype = point.shape, point.dtype
     cuda.memcpy_htod(int(struct_ptr), np.int32(code.size))
     cuda.memcpy_htod(int(struct_ptr) + 8, np.intp(int(self.code)))
     cuda.memcpy_htod(int(struct_ptr) + 8 + np.intp(0).nbytes, np.intp(int(self.point)))

 def test_intp(self,level=rlevel):
     """Ticket #99"""
     i_width = np.int_(0).nbytes*2 - 1
     np.intp('0x' + 'f'*i_width,16)
     self.assertRaises(OverflowError,np.intp,'0x' + 'f'*(i_width+1),16)
     self.assertRaises(ValueError,np.intp,'0x1',32)
     assert_equal(255,np.intp('0xFF',16))
     assert_equal(1024,np.intp(1024))

 def check_intp(self,level=rlevel):
     """Ticket #99"""
     i_width = N.int_(0).nbytes*2 - 1
     N.intp('0x' + 'f'*i_width,16)
     self.failUnlessRaises(OverflowError,N.intp,'0x' + 'f'*(i_width+1),16)
     self.failUnlessRaises(ValueError,N.intp,'0x1',32)
     assert_equal(255,N.intp('0xFF',16))
     assert_equal(1024,N.intp(1024))

    def getSRT(self, gmap,  store_srt=False):
        """
        Computes the sample rank templates for the expression matrix(on instantiation) and 
        gmap

        gmap is a 1d numpy array where gmap[2*i] and gmap[2*i +1] are 
            gene indices for comparison i

        b_size is the block size to use in gpu computation

        store_srt - determines what is returned,
            False(default) - returns the srt numpy array (npairs,nsamp)
            True - returns the srt_gpu object and the object's padded shape (npairs,nsamp)
        """
        #the x coords in the gpu map to sample_ids
        #the y coords to gmap
        #sample blocks
        b_size = self.b_size
        exp = self.exp 
        g_y_sz = self.getGrid( exp.shape[1] )
        #pair blocks
        g_x_sz = self.getGrid( gmap.shape[0]/2 )

        #put gene map on gpu
        gmap_buffer = self.gmap_buffer= self.getBuff(gmap, 2*(g_x_sz*b_size), 1,np.int32)
 
        gmap_gpu = np.intp(gmap_buffer.base.get_device_pointer()) #cuda.mem_alloc(gmap_buffer.nbytes)
        #cuda.memcpy_htod(gmap_gpu,gmap_buffer)
        #make room for srt
        srt_shape = (g_x_sz*b_size , g_y_sz*b_size)
        srt_buffer = self.srt_buffer = self.getBuff(np.zeros(srt_shape, dtype=np.int32),srt_shape[0],srt_shape[1], np.int32)
        srt_gpu = np.intp(srt_buffer.base.get_device_pointer()) #cuda.mem_alloc(srt_shape[0]*srt_shape[1]*np.int32(1).nbytes)
        srtKern = self.getsrtKern()
        

        exp_gpu = self.exp_gpu
        nsamp = np.uint32( g_y_sz * b_size )
        ngenes = np.uint32( self.exp.shape[0] )
        npairs = np.uint32( g_x_sz * b_size )

        block = (b_size,b_size,1)
        grid = (g_x_sz, g_y_sz)

        srtKern(exp_gpu, nsamp, ngenes, gmap_gpu, npairs, srt_gpu, block=block, grid=grid)

        #gmap_gpu.free()
        
       
        if store_srt:
            #this is in case we want to run further stuff without 
            #transferring back and forth
            return (srt_gpu, npairs , nsamp)
        else:
            #srt_buffer = np.zeros(srt_shape, dtype=np.int32)
            #cuda.memcpy_dtoh(srt_buffer, srt_gpu)
            #srt_gpu.free()

            return srt_buffer[:gmap.shape[0]/2,:self.exp.shape[1]]

 def _default_norm(self, layer):
     vals = np.sort(layer.ravel())
     vals = vals[np.isfinite(vals)]
     result = DS9Normalize()
     result.stretch = 'arcsinh'
     result.clip = True
     if vals.size > 0:
         result.vmin = vals[np.intp(.01 * vals.size)]
         result.vmax = vals[np.intp(.99 * vals.size)]
     return result

def readout(mesh, pos, mode="raise", period=None, transform=None, out=None):
    """ CIC approximation, reading out mesh values at pos,
        see document of paint. 
    """
    pos = numpy.array(pos)
    if out is None:
        out = numpy.zeros(len(pos), dtype='f8')
    else:
        out[:] = 0
    chunksize = 1024 * 16 * 4
    Ndim = pos.shape[-1]
    Np = pos.shape[0]
    if transform is None:
        transform = lambda x: x

    neighbours = ((numpy.arange(2 ** Ndim)[:, None] >> \
            numpy.arange(Ndim)[None, :]) & 1)
    for start in range(0, Np, chunksize):
        chunk = slice(start, start+chunksize)

        if mode == 'raise':
            gridpos = transform(pos[chunk])
            rmi_mode = 'raise'
            intpos = numpy.intp(numpy.floor(gridpos))
        elif mode == 'ignore':
            gridpos = transform(pos[chunk])
            rmi_mode = 'raise'
            intpos = numpy.intp(numpy.floor(gridpos))

        for i, neighbour in enumerate(neighbours):
            neighbour = neighbour[None, :]

            targetpos = intpos + neighbour

            kernel = (1.0 - numpy.abs(gridpos - targetpos)).prod(axis=-1)

            if period is not None:
                period = numpy.int32(period)
                numpy.remainder(targetpos, period, targetpos)

            if mode == 'ignore':
                # filter out those outside of the mesh
                mask = (targetpos >= 0).all(axis=-1)
                for d in range(Ndim):
                    mask &= (targetpos[..., d] < mesh.shape[d])
                targetpos = targetpos[mask]
                kernel = kernel[mask]
            else:
                mask = Ellipsis

            if len(targetpos) > 0:
                targetindex = numpy.ravel_multi_index(
                        targetpos.T, mesh.shape, mode=rmi_mode)
                out[chunk][mask] += kernel * mesh.flat[targetindex]
    return out

    def _loadInput(self, stim):
        logging.debug('loadInput')

        # shortcuts
        nrXY = self.nrX * self.nrY
        nrXYD = self.nrX * self.nrY * self.nrDirs

        # parse input
        assert type(stim).__module__ == "numpy", "stim must be numpy array"
        assert type(stim).__name__ == "ndarray", "stim must be numpy.ndarray"
        assert stim.size > 0, "stim cannot be []"
        stim = stim.astype(np.ubyte)

        rows, cols = stim.shape
        logging.debug("- stim shape={0}x{1}".format(rows, cols))

        # shift d_stimBuf in time by 1 frame, from frame i to frame i-1
        # write our own memcpy kernel... :-(
        gdim = (int(iDivUp(nrXY, 128)), 1)
        bdim = (128, 1, 1)
        for i in xrange(1, self.nrT):
            stimBufPt_dst = np.intp(self.d_stimBuf) + self.szXY * (i - 1)
            stimBufPt_src = np.intp(self.d_stimBuf) + self.szXY * i
            self.dev_memcpy_dtod(
                stimBufPt_dst,
                stimBufPt_src,
                np.int32(nrXY),
                block=bdim, grid=gdim)

        # index into d_stimBuf array to place the new stim at the end
        # (newest frame at pos: nrT-1)
        d_stimBufPt = np.intp(self.d_stimBuf) + self.szXY * (self.nrT-1)

        # \TODO implement RGB support
        self.dev_split_gray(
            d_stimBufPt,
            cuda.In(stim),
            np.int32(stim.size),
            block=bdim, grid=gdim)

        # create working copy of d_stimBuf
        cuda.memcpy_dtod(self.d_scalingStimBuf, self.d_stimBuf,
                         self.szXY*self.nrT)

        # reset V1complex responses to 0
        # \FIXME not sure how to use memset...doesn't seem to give expected
        # result
        tmp = np.zeros(nrXYD).astype(np.float32)
        cuda.memcpy_htod(self.d_respV1c, tmp)

        # allocate d_resp, which will contain the response to all 28
        # (nrFilters) space-time orientations at 3 (nrScales) scales for
        # every pixel location (nrX*nrY)
        tmp = np.zeros(nrXY*self.nrFilters*self.nrScales).astype(np.float32)
        cuda.memcpy_htod(self.d_resp, tmp)

    def getRT(self, s_map, srt_gpu, srt_nsamp, srt_npairs, npairs, store_rt=False):
        """
        Computes the rank template

        s_map(Sample Map) -  an list of 1s and 0s of length nsamples where 1 means use this sample
            to compute rank template
        srt_gpu - cuda memory object containing srt(sample rank template) array on gpu
        srt_nsamp, srt_npairs - shape(buffered) of srt_gpu object
        npairs - true number of gene pairs being compared
        b_size - size of the blocks for computation
        store_rt - determines the RETURN value
            False(default) = returns an numpy array shape(npairs) of the rank template
            True = returns the rt_gpu object and the padded size of the rt_gpu objet (rt_obj, npairs_padded)
        """

        b_size = self.b_size
        s_map_buff = self.s_map_buff = cuda.pagelocked_zeros((int(srt_nsamp),), np.int32,  mem_flags=cuda.host_alloc_flags.DEVICEMAP)

        s_map_buff[:len(s_map)] =  np.array(s_map,dtype=np.int32)

        s_map_gpu = np.intp(s_map_buff.base.get_device_pointer())
        #cuda.memcpy_htod(s_map_gpu, s_map_buff)
        
        #sample blocks
        g_y_sz = self.getGrid( srt_nsamp)
        #pair blocks
        g_x_sz = self.getGrid( srt_npairs )
        
        block_rt_gpu =  cuda.mem_alloc(int(g_y_sz*srt_npairs*(np.uint32(1).nbytes)) ) 


        grid = (g_x_sz, g_y_sz)

        func1,func2 = self.getrtKern(g_y_sz)

        shared_size = b_size*b_size*np.uint32(1).nbytes

        func1( srt_gpu, np.uint32(srt_nsamp), np.uint32(srt_npairs), s_map_gpu, block_rt_gpu, np.uint32(g_y_sz), block=(b_size,b_size,1), grid=grid, shared=shared_size)

        rt_buffer =self.rt_buffer = cuda.pagelocked_zeros((int(srt_npairs),), np.int32, mem_flags=cuda.host_alloc_flags.DEVICEMAP)
        rt_gpu = np.intp(rt_buffer.base.get_device_pointer())

        func2( block_rt_gpu, rt_gpu, np.int32(s_map_buff.sum()), block=(b_size,1,1), grid=(g_x_sz,))

        
        if store_rt:
            #this is in case we want to run further stuff without 
            #transferring back and forth
            return (rt_gpu, srt_npairs)
        else:
            #rt_buffer = np.zeros((srt_npairs ,), dtype=np.int32)
            #cuda.memcpy_dtoh(rt_buffer, rt_gpu)
            #rt_gpu.free()
            return rt_buffer[:npairs]

 def coarsegrain(self,flowy,deltafy,Ny):
     if not hasattr(self,'Moments'):
         self.getMultipleMoments(msign='pn')
     Nx = self.Length
     flowx = self.FreqOffset
     deltafx = self.Cadence
     if ( (deltafx <= 0) | (deltafy <= 0) | (Ny <= 0) ):
         raise ValueError, 'bad input argument'
     if ( deltafy < deltafx ):
         raise ValueError, 'deltaf coarse-grain < deltaf fine-grain'
     if ( (flowy - 0.5*deltafy) < (flowx - 0.5*deltafx) ):
         raise ValueError, 'desired coarse-grained start frequency is too low'
     fhighx = flowx + (Nx-1)*deltafx
     fhighy = flowy + (Ny-1)*deltafy
     if ( (fhighy + 0.5*deltafy) > (fhighx + 0.5*deltafx) ):
         raise ValueError, 'desired coarse-rained stop frequency is too high'
     i = numpy.arange(Ny)
     jlow = numpy.intp(
         1 + numpy.floor((flowy + (i-0.5)*deltafy - flowx - 0.5*deltafx)/deltafx))
     jhigh = numpy.intp(
         1 + numpy.floor((flowy + (i+0.5)*deltafy - flowx - 0.5*deltafx)/deltafx))
     index1 = jlow[0]
     index2 = jhigh[-1]
     fraclow = (flowx + (jlow+0.5)*deltafx - flowy - (i-0.5)*deltafy)/deltafx
     frachigh = (flowy + (i+0.5)*deltafy - flowx - (jhigh-0.5)*deltafx)/deltafx
     frac1 = fraclow[0]
     frac2 = frachigh[-1]
     jtemp = jlow + 1
     coarseMoments = numpy.zeros( (numpy.shape(self.Moments)[0],Ny) , complex)
     for lm in range(numpy.shape(self.Moments)[0]):
         midsum = sumTerms(self.Moments[lm,:], jtemp, jhigh)
         ya = (deltafx/deltafy)*(self.Moments[lm,:][jlow[:-1]]*fraclow[:-1] +
                                 self.Moments[lm,:][jhigh[:-1]]*frachigh[:-1] +
                                 midsum[:-1])
         if (jhigh[-1] > Nx-1):
             yb = (deltafx/deltafy)*(self.Moments[lm,:][jlow[-1]]*fraclow[-1] + midsum[-1])
         else:
             yb = (deltafx/deltafy)*(self.Moments[lm,:][jlow[-1]]*fraclow[-1] +
                                     self.Moments[lm,:][jhigh[-1]]*frachigh[-1] +
                                     midsum[-1])
         coarseMoments[lm,:] = numpy.array( list(ya) + [yb] )
     self.coarseMoments = coarseMoments
     self.coarseFreqOffset = flowy
     self.coarseCadence = deltafy
     self.coarseLength = numpy.shape(self.coarseMoments)[1]
     self.index1 = index1
     self.index2 = index2
     self.frac1 = frac1
     self.frac2 = frac2
     return coarseMoments

    def _build_arg_buf(args):
        handlers = []

        arg_data = []
        format = ""
        for i, arg in enumerate(args):
            if isinstance(arg, np.number):
                arg_data.append(arg)
                format += arg.dtype.char
            elif isinstance(arg, (DeviceAllocation, PooledDeviceAllocation)):
                arg_data.append(int(arg))
                format += "P"
            elif isinstance(arg, ArgumentHandler):
                handlers.append(arg)
                arg_data.append(int(arg.get_device_alloc()))
                format += "P"
            elif isinstance(arg, np.ndarray):
                arg_data.append(arg)
                format += "%ds" % arg.nbytes
            else:
                try:
                    gpudata = np.intp(arg.gpudata)
                except AttributeError:
                    raise TypeError("invalid type on parameter #%d (0-based)" % i)
                else:
                    # for gpuarrays
                    arg_data.append(int(gpudata))
                    format += "P"

        from pycuda._pvt_struct import pack
        return handlers, pack(format, *arg_data)

 def __new__(cls, x=0):
     if isinstance(x, afnumpy.ndarray):
         return x.astype(cls)
     elif isinstance(x, numbers.Number):
         return numpy.intp(x)
     else:
         return afnumpy.array(x).astype(cls)

 def test_numpy(self):
     """NumPy objects get serialized to readable JSON."""
     l = [
         np.float32(12.5),
         np.float64(2.0),
         np.float16(0.5),
         np.bool(True),
         np.bool(False),
         np.bool_(True),
         np.unicode_("hello"),
         np.byte(12),
         np.short(12),
         np.intc(-13),
         np.int_(0),
         np.longlong(100),
         np.intp(7),
         np.ubyte(12),
         np.ushort(12),
         np.uintc(13),
         np.ulonglong(100),
         np.uintp(7),
         np.int8(1),
         np.int16(3),
         np.int32(4),
         np.int64(5),
         np.uint8(1),
         np.uint16(3),
         np.uint32(4),
         np.uint64(5),
     ]
     l2 = [l, np.array([1, 2, 3])]
     roundtripped = loads(dumps(l2, cls=EliotJSONEncoder))
     self.assertEqual([l, [1, 2, 3]], roundtripped)