def formal_integral_model(request, model):
    r = request.param['r']
    model.no_of_shells_i = r.shape[0] - 1
    model.inverse_time_explosion = c.c.cgs.value
    model.r_outer_i.contents = as_ctypes(r[1:])
    model.r_inner_i.contents = as_ctypes(r[:-1])
    return model

def setup_stage():
	"""docstring for setup_stage
	"""
	import ctypes
	
	e7xx = ctypes.windll.LoadLibrary('E7XX_GCS_DLL.dll')
	try:
		print "Connecting to stage"
		id = e7xx.E7XX_ConnectRS232(1, 57600)

		print "Initializing axes"
		e7xx.E7XX_INI(id, '134')

		print "initializing servos"
		err = e7xx.E7XX_SVO(id, '134', ctl.as_ctypes(ones(4, dtype=int32)))
		if err:
			print "Servos initialized OK"
		else:
			import sys
			sys.exit(e7xx.E7XX_GetError(id))
		svo = ctl.as_ctypes(ones(4, dtype=int32))
		err = e7xx.E7XX_qSVO(id, '134', svo)
		if err:
			print "Read servos OK"
		else:
			print e7xx.E7XX_GetError(id)
			time.sleep(5)
		
		while not(all(ctl.as_array(svo))):
			e7xx.E7XX_qSVO(id, '134', svo)
			print "Servo status: ", ctl.as_array(svo), ctl.as_array(svo).all()
			time.sleep(1)

	finally:
		return e7xx, id

 def dispose_gl(self):
     glDeleteTextures(self.texture_id)
     glDeleteRenderbuffers(1, as_ctypes(self.depth_buffer))
     glDeleteFramebuffers(1, as_ctypes(self.fb))
     self.fb = 0
     if self.multisample > 0:
         glDeleteTextures(self.resolve_texture_id)
         glDeleteFramebuffers(1, as_ctypes(self.resolve_fb))

def LUSolve(LU, b):
    """
    Resolve LU = b.
    """
    n = LU.shape[0]
    x = b.copy()
    lib.LUSolve(n, byref(ctypeslib.as_ctypes(LU)),
                byref(ctypeslib.as_ctypes(x)),
                byref(ctypeslib.as_ctypes(b)))
    return x

def do_some_task():
    SIZE = 1e8
    
    input_array = 1 * np.random.random(SIZE).astype(np.float32)
    output_array = np.empty_like(input_array)
    
    lib.do_some_omp_task(as_ctypes(input_array),
                         as_ctypes(output_array),
                         ct.c_size_t(input_array.size))
    return output_array

def LUCroutDecompose(A):
    """
    Implementação do método de Crout para decomposição LU.
    """
    assert A.shape[0] == A.shape[1] and type(A) is matrix, "'A' deve ser NxN."
    L = zeros(A.shape)
    n = A.shape[0]
    U = L.copy()
    lib.LUDec(n, byref(ctypeslib.as_ctypes(A)),
              byref(ctypeslib.as_ctypes(L)),
              byref(ctypeslib.as_ctypes(U)))
    return L, U

def move_and_image(mmc, e7xx, id, coords, exptime, image_queue, **kwargs):
    """
    move_and_image moves the stage to the given coordinates, takes an
    exposure for exptime seconds, then adds it to image_queue.
    """

    DEBUG = False

    for key in kwargs:
        if key == "DEBUG":
            DEBUG = True
        else:
            raise TypeError, 'Unknown argument "%s"' % key

    if DEBUG:
        print "Moving to ", coords
    err = e7xx.E7XX_MOV(id, "14", ctl.as_ctypes(array(coords, dtype=float)))

    if err:
        print "Moved OK"
    else:
        err = e7xx.E7XX_GetError(id)
        print err

    time.sleep(0.03)
    if DEBUG:
        res = ctl.as_ctypes(empty(4, dtype=float))
        e7xx.E7XX_qMOV(id, "14", res)
        print "Moved to ", ctl.as_array(res)

    noImage = True

    while noImage:
        try:
            if image_queue.qsize() < 1000:
                if DEBUG:
                    print "Snapping Image"
                mmc.snapImage()
                im1 = mmc.getImage()
                if DEBUG:
                    print "Got image"
                image_queue.put(im1)
                if DEBUG:
                    print "Queueing image"
                noImage = False
                if DEBUG:
                    print "Leaving Loop"
        except MemoryError:
            if DEBUG:
                print "Memory Error.  Going to sleep"
            time.sleep(1)
    if DEBUG:
        print "Done"

def GaussSeidel(A, b, ks=50):
    """
    Resolve Ax = b através do método iterativo de Jacobi.
    """
    assert A.shape[0] == A.shape[1] and type(A) is matrix, "'A' deve ser NxN."
    assert b.shape[0] == A.shape[0], "'b' deve ser compatível com A."
    n = A.shape[0]
    x = zeros(n)
    try:
        lib.SolveGaussSeidel(n, ks, byref(ctypeslib.as_ctypes(A)),
                             byref(ctypeslib.as_ctypes(x)),
                             byref(ctypeslib.as_ctypes(b)))
    except:
        raise
    return x

def Jacobi(A, b, ks=1000, cuda=CUDAcapable):
    """
    Resolve Ax = b através do método iterativo de Jacobi.
    """
    assert A.shape[0] == A.shape[1] and type(A) is matrix, "'A' deve ser NxN."
    assert b.shape[0] == A.shape[0], "'b' deve ser compatível com A."
    n = A.shape[0]
    x = zeros(n)
    func = kernels.CUJacobi if cuda else lib.SolveJacobi
    try:
        func(n, ks, byref(ctypeslib.as_ctypes(A)),
             byref(ctypeslib.as_ctypes(x)),
             byref(ctypeslib.as_ctypes(b)))
    except:
        raise
    return x

def GaussJordan(A, b):
    """
    Resolve Ax = b através do método de eliminação de Gauss-Jordan com
    pivotação.
    """
    assert A.shape[0] == A.shape[1] and type(A) is matrix, "'A' deve ser NxN."
    assert b.shape[0] == A.shape[0], "'b' deve ser compatível com A."
    n = A.shape[0]
    x = zeros(n)
    try:
        lib.SolveGJ(n, byref(ctypeslib.as_ctypes(A)),
                    byref(ctypeslib.as_ctypes(x)),
                    byref(ctypeslib.as_ctypes(b)))
    except:
        raise
    return x

def arrsort (arr):
    base_ptr = POINTER (as_ctypes (arr)._type_)
    functype = CFUNCTYPE (c_int, base_ptr, base_ptr)
    
    nmemb, = arr.shape
    size,  = arr.strides
    cmpfun = functype (_cmpelem)
    _lib.qsort (arr, nmemb, size, cast (cmpfun, _cmptype))

    def dedisperse_gpu(self,output_dir=".",out_bits=32,gulp=160000):
        ndm = self.get_dm_count()
        delay = self.get_max_delay()
        if gulp-delay < 2*delay:
            gulp = 2*delay

        if out_bits == 32:
            dtype = "float32"
        elif out_bits == 8:
            dtype = "ubyte"
        else:
            raise ValueError("out_bits must be 8 or 32")

        outsamps = gulp-delay

        print outsamps,out_bits,dtype
        

        dms = self.get_dm_list()
        out_files = []
        changes = {"refdm" :0, "nchans":1, "nbits" :out_bits}
        basename = self.header.basename.split("/")[-1]
        for dm in dms:
            changes["refdm"] = dm
            filename = "%s/%s_DM%08.2f.tim"%(output_dir,basename,dm)
            out_files.append(self.header.prepOutfile(filename,changes,nbits=out_bits))
            
        out_size = outsamps * ndm * out_bits/8
        output = np.empty(out_size,dtype=dtype)
        func = lib.dedisp_execute
        for nsamps,ii,data in self.readPlan(gulp=gulp,skipback=delay):
            error = func(self.plan,
                         C.c_size_t(nsamps),
                         as_ctypes(data),
                         C.c_size_t(self.header.nbits),
                         as_ctypes(output),
                         C.c_size_t(out_bits),
                         C.c_int(0));
            error_check(error)
            for ii,out_file in enumerate(out_files):
                out_file.cwrite(output[ii*outsamps:(ii+1)*outsamps])
        
        for out_file in out_files:
           out_file.close()

def LUCroutInplaceDecompose(A):
    """
    Implementação do método de Crout para decomposição LU sobrescrevendo a
    matriz original.
    """
    assert A.shape[0] == A.shape[1] and type(A) is matrix, "'A' deve ser NxN."
    LU = A.copy()
    n = A.shape[0]
    lib.LUInDec(n, byref(ctypeslib.as_ctypes(LU)))
    return LU

 def getRAM(self,ram=None):
     """This function grabs the atari RAM.
     ram MUST be a numpy array of uint8/int8. This can be initialized like so:
     ram = np.array(ram_size,dtype=uint8)
     Notice: It must be ram_size where ram_size can be retrieved via the getRAMSize function.
     If it is None, then this function will initialize it.
     """
     if(ram is None):
         ram_size = ale_lib.getRAMSize(self.obj)
         ram = np.zeros(ram_size,dtype=np.uint8)
     ale_lib.getRAM(self.obj,as_ctypes(ram))

 def getScreenGrayscale(self, screen_data=None):
     """This function fills screen_data with the data in grayscale
     screen_data MUST be a numpy array of uint8. This can be initialized like so:
     screen_data = np.empty((height,width,1), dtype=np.uint8)
     If it is None,  then this function will initialize it.
     """
     if(screen_data is None):
         width = ale_lib.getScreenWidth(self.obj)
         height = ale_lib.getScreenHeight(self.obj)
         screen_data = np.empty((height, width,1), dtype=np.uint8)
     ale_lib.getScreenGrayscale(self.obj, as_ctypes(screen_data[:]))
     return screen_data

 def getScreenRGB(self, screen_data=None):
     """This function fills screen_data with the data
     screen_data MUST be a numpy array of uint/int. This can be initialized like so:
     screen_data = np.array(w*h, dtype=np.intc)
     Notice,  it must be width*height in size also
     If it is None,  then this function will initialize it
     """
     if(screen_data is None):
         width = ale_lib.getScreenWidth(self.obj)
         height = ale_lib.getScreenHeight(self.obj)
         screen_data = np.zeros(width*height, dtype=np.intc)
     ale_lib.getScreenRGB(self.obj, as_ctypes(screen_data))
     return screen_data

 def getScreen(self, screen_data=None):
     """This function fills screen_data with the RAW Pixel data
     screen_data MUST be a numpy array of uint8/int8. This could be initialized like so:
     screen_data = np.empty(w*h, dtype=np.uint8)
     Notice,  it must be width*height in size also
     If it is None,  then this function will initialize it
     Note: This is the raw pixel values from the atari,  before any RGB palette transformation takes place
     """
     if(screen_data is None):
         width = ale_lib.getScreenWidth(self.obj)
         height = ale_lib.getScreenHeight(self.obj)
         screen_data = np.zeros(width*height, dtype=np.uint8)
     ale_lib.getScreen(self.obj, as_ctypes(screen_data))
     return screen_data

 def getScreenRGB2(self, screen_data=None):
     """This function fills screen_data with the data in RGB format.
     screen_data MUST be a numpy array of uint8. This can be initialized like so:
       screen_data = np.empty((height,width,3), dtype=np.uint8)
     If it is None,  then this function will initialize it.
     On all architectures, the channels are RGB order:
         screen_data[x, y, :] is [red, green, blue]
     There's not much error checking here: if you supply an array that's too small
     this function will produce undefined behavior.
     """
     if(screen_data is None):
         width = ale_lib.getScreenWidth(self.obj)
         height = ale_lib.getScreenHeight(self.obj)
         screen_data = np.empty((height, width, 3), dtype=np.uint8)
     assert screen_data.strides == (480, 3, 1)
     ale_lib.getScreenRGB2(self.obj, as_ctypes(screen_data[:]))
     return screen_data

 def getScreenRGB(self, screen_data=None):
     """This function fills screen_data with the data in RGB format
     screen_data MUST be a numpy array of uint8. This can be initialized like so:
     screen_data = np.empty((height,width,3), dtype=np.uint8)
     If it is None,  then this function will initialize it.
     On little-endian machines like x86, the channels are BGR order:
         screen_data[x, y, 0:3] is [blue, green, red]
     On big-endian machines (rare in 2017) the channels would be the opposite order.
     There's not much error checking here: if you supply an array that's too small
     this function will produce undefined behavior.
     """
     if(screen_data is None):
         width = ale_lib.getScreenWidth(self.obj)
         height = ale_lib.getScreenHeight(self.obj)
         screen_data = np.empty((height, width,3), dtype=np.uint8)
     ale_lib.getScreenRGB(self.obj, as_ctypes(screen_data[:]))
     return screen_data

 def from_array(self, arr, _ctype=ctype):
     "set up a descriptor from a compatible numpy array, preferably a fortran-order numpy array"
     dims = arr.shape
     rank = len(dims)
     nelems = reduce(mul, dims, 1)
     assert(len(self.dim) == rank)
     ctarr = as_ctypes(arr)
     # compatibility sanity check
     assert(arr.dtype == _ctype)
     # ctype = _base_ctype(ctarr, arr.shape)
     ct = cast(ctarr, type(self.mem))
     csize = sizeof(_ctype)
     self.mem = ct
     self.dtype_rank = rank
     self.dtype_size = csize
     self.dtype_type = GFC_DTYPE_TAB.get(_ctype, GFC_DTYPE_DERIVED)
     axes = list(self.dim)
     for width, axis, stride in zip(dims, axes, arr.strides):
         axis.stride = stride / csize
         axis.lbound = 1
         axis.ubound = width