def _fold_exp_and_coh(t_array, w, tz, tau_arr):
    if tz != 0.:
        t_array -= tz

    shape = t_array.shape
    t_array = t_array.astype(np.float32)

    t_arr_gpu = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=t_array)
    tau_buf = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR,
                        hostbuf=(1/tau_arr).astype(np.float32))
    shape = (shape[0], shape[1], tau_arr.size)
    shape_coh = (shape[0], shape[1], 3)
    out = cl_array.empty(queue, shape=shape, dtype=np.float32)
    out_coh = cl_array.empty(queue, shape=shape_coh, dtype=np.float32)

    global_work_size = t_array.size + (work_size[0] - t_array.size % work_size[0])

    prg.fold_exp(queue, (global_work_size, tau_arr.size), work_size, t_arr_gpu, np.float32(w),
                 tau_buf, out.data, np.uint32(t_array.size))

    coh_no_div.coh_gauss(queue, (global_work_size, 3), work_size, t_arr_gpu,
                   np.float32(w/1.4142), out_coh.data, np.uint32(t_array.size))

    queue.finish()
    a = out.get(async_=True)
    b = out_coh.get(async_=True)
    b /= np.abs(b).max(0)
    queue.finish()
    return a, b

    def setup_arrays(self, nrays, nsamples, cutoff):

        prog_params = (nrays, nsamples, cutoff)

        if prog_params in self.array_cache:
            return self.array_cache[prog_params]

        else:
            arrays = ArraySet()
            arrays.scratch = cla.empty(self.queue,
                                 (nsamples, nrays),
                                 dtype=np.float32,
                                 allocator=self.memory_pool)

            arrays.result = cla.empty(self.queue,
                                (nrays,),
                                dtype=np.int32,
                                allocator=self.memory_pool)

            arrays.pre_cutoff = cla.empty(self.queue,
                                    (nrays, cutoff),
                                    dtype=np.float32,
                                    allocator=self.memory_pool)

            arrays.pre_cutoff_squared = cla.empty_like(arrays.pre_cutoff)

            arrays.idx = cla.arange(self.queue, 0, cutoff * nrays, 1,
                                    dtype=np.int32,
                                    allocator=self.memory_pool)

            self.array_cache[prog_params] = arrays
            return arrays

def test_numpy_integer_shape(ctx_factory):
    try:
        list(np.int32(17))
    except:
        pass
    else:
        from pytest import skip
        skip("numpy implementation does not handle scalar correctly.")
    context = ctx_factory()
    queue = cl.CommandQueue(context)

    cl_array.empty(queue, np.int32(17), np.float32)
    cl_array.empty(queue, (np.int32(17), np.int32(17)), np.float32)

    def _prep_gpu():
        """ Set up GPU calculation dependencies """

        # try to import the necessary libraries
        fallback = False
        try:
            import gpu
            import string
            import pyopencl as cl
            import pyopencl.array as cla
            from pyfft.cl import Plan
        except ImportError:
            fallback = True
            
        # check gpu_info
        try:
            assert gpu.valid(gpu_info),\
            "gpu_info in propagate_distances improperly specified"
            
            context, device, queue, platform = gpu_info
        except AssertionError:
            fallback = True
            
        if fallback:
            propagate_distances(data, distances, energy_or_wavelength,
                                pixel_pitch, subregion=subregion,
                                silent=silent, band_limit=band_limit,
                                gpu_info=None, im_convert=im_convert)
    
        # if everything is OK, allocate memory and build kernels
        kp = string.join(gpu.__file__.split('/')[:-1], '/')+'/kernels/'
        build = _build_helper(context, device, kp)
        phase_multiply = build('propagate_phase_multiply.cl')
        copy_to_buffer = build('propagate_copy_to_save_buffer.cl')
        fftplan = Plan((N, N), queue=queue)

        # put the signals onto the gpu along with buffers for the
        # various operations
        rarray = cla.to_device(queue, r.astype(np.float32))
        fourier = cla.to_device(queue, data.astype(np.complex64))
        phase = cla.empty(queue, (N, N), np.complex64)
        back = cla.empty(queue, (N, N), np.complex64)
        store = cla.empty(queue, (nf, rows, cols), np.complex64)
        
        # precompute the fourier transform of data. 
        fftplan.execute(fourier.data, wait_for_finish=True)

        return phase_multiply, copy_to_buffer, fftplan, rarray, fourier,\
               phase, back, store, build

    def build_scratch(self, imshape):

        self.scratch = []
        self.index_scratch = []

        l = np.prod(imshape)
        self.array_indices = cla.arange(self.queue, 0, l, 1, dtype=np.int32)
        if l % self.runlen != 0:
            l += l % self.runlen
        while l > 1:
            l /= self.runlen
            self.scratch.append(cla.empty(self.queue, (l,), np.float32))
            self.index_scratch.append(cla.empty(self.queue, (l,), np.int32))

        self.imshape = imshape

def axis_convolve(X, h, axis=0, queue=None, output=None):
    """Filter along an of *X* using filter vector *h*.  If *h* has odd length, each
    output sample is aligned with each input sample and *Y* is the same size as
    *X*.  If *h* has even length, each output sample is aligned with the mid point
    of each pair of input samples, and the output matrix's shape is increased
    by one along the convolution axis.

    After convolution, the :py:class:`pyopencl.array.Array` instance holding the
    device-side output is returned. This may be accessed on the host via
    :py:func:`to_array`.

    The axis of convolution is specified by *axis*. The default direction of
    convolution is column-wise.

    If *queue* is non-``None``, it should be a :py:class:`pyopencl.CommandQueue`
    instance which is used to perform the computation. If ``None``, a default
    global queue is used.

    If *output* is non-``None``, it should be a :py:class:`pyopencl.array.Array`
    instance which the result is written into. If ``None``, an output array is
    created.
    """

    _check_cl()
    queue = to_queue(queue)
    kern = _convolve_kernel_for_queue(queue.context)

    # Create output if not specified
    if output is None:
        output_shape = list(X.shape)
        if h.shape[0] % 2 == 0:
            output_shape[axis] += 1
        output = cl_array.empty(queue, output_shape, np.float32)

    return _apply_kernel(X, h, kern, output, axis=axis)

def test_index_preservation(ctx_factory):
    from pytest import importorskip
    importorskip("mako")

    context = ctx_factory()
    queue = cl.CommandQueue(context)

    from pyopencl.scan import GenericScanKernel, GenericDebugScanKernel
    classes = [GenericScanKernel]

    dev = context.devices[0]
    if dev.type & cl.device_type.CPU:
        classes.append(GenericDebugScanKernel)

    for cls in classes:
        for n in scan_test_counts:
            knl = cls(
                    context, np.int32,
                    arguments="__global int *out",
                    input_expr="i",
                    scan_expr="b", neutral="0",
                    output_statement="""
                        out[i] = item;
                        """)

            out = cl_array.empty(queue, n, dtype=np.int32)
            knl(out)

            assert (out.get() == np.arange(n)).all()
            from gc import collect
            collect()

 def _dev_array(self):
     if not hasattr(self, '__dev_array'):
         setattr(self, '__dev_array',
                 array.empty(_queue,
                             self.sparsity.nz,
                             self.dtype))
     return getattr(self, '__dev_array')

def get_array(data, queue=None):
    """Get pyopencl.array.Array from *data* which can be a numpy array, a pyopencl.array.Array or a
    pyopencl.Image. *queue* is an OpenCL command queue.
    """
    if not queue:
        queue = cfg.OPENCL.queue

    if isinstance(data, cl_array.Array):
        result = data
    elif isinstance(data, np.ndarray):
        if data.dtype.kind == 'c':
            if data.dtype.itemsize != cfg.PRECISION.cl_cplx:
                data = data.astype(cfg.PRECISION.np_cplx)
            result = cl_array.to_device(queue, data.astype(cfg.PRECISION.np_cplx))
        else:
            if data.dtype.kind != 'f' or data.dtype.itemsize != cfg.PRECISION.cl_float:
                data = data.astype(cfg.PRECISION.np_float)
            result = cl_array.to_device(queue, data.astype(cfg.PRECISION.np_float))
    elif isinstance(data, cl.Image):
        result = cl_array.empty(queue, data.shape[::-1], np.float32)
        cl.enqueue_copy(queue, result.data, data, offset=0, origin=(0, 0),
                        region=result.shape[::-1])
        if result.dtype.itemsize != cfg.PRECISION.cl_float:
            result = result.astype(cfg.PRECISION.np_float)
    else:
        raise TypeError('Unsupported data type {}'.format(type(data)))

    return result

    def uniform(self, *args, **kwargs):
        a = kwargs.pop("a", 0)
        b = kwargs.pop("b", 1)

        result = cl_array.empty(*args, **kwargs)

        self.fill_uniform(result, queue=result.queue, a=a, b=b)
        return result

    def normal(self, *args, **kwargs):
        mu = kwargs.pop("mu", 0)
        sigma = kwargs.pop("sigma", 1)

        result = cl_array.empty(*args, **kwargs)

        self.fill_normal(result, queue=result.queue, mu=mu, sigma=sigma)
        return result

 def _allocate_device(self):
     if self.state is DeviceDataMixin.DEVICE_UNALLOCATED:
         if self.soa:
             shape = self._data.T.shape
         else:
             shape = self._data.shape
         self._device_data = array.empty(_queue, shape=shape, dtype=self.dtype)
         self.state = DeviceDataMixin.HOST

    def __init__(self, queue, num_work_items,
            luxury=None, seed=None, no_warmup=False,
            use_legacy_init=False, max_work_items=None):
        if luxury is None:
            luxury = 4

        if seed is None:
            from time import time
            seed = int(time()*1e6) % 2<<30

        self.context = queue.context
        self.luxury = luxury
        self.num_work_items = num_work_items

        from pyopencl.characterize import has_double_support
        self.support_double = has_double_support(queue.device)

        self.no_warmup = no_warmup
        self.use_legacy_init = use_legacy_init
        self.max_work_items = max_work_items

        src = """
            %(defines)s

            #include <pyopencl-ranluxcl.cl>

            kernel void init_ranlux(unsigned seeds, global ranluxcl_state_t *ranluxcltab)
            {
              if (get_global_id(0) < %(num_work_items)d)
                ranluxcl_initialization(seeds, ranluxcltab);
            }
            """ % {
                    "defines": self.generate_settings_defines(),
                    "num_work_items": num_work_items
                }
        prg = cl.Program(queue.context, src).build()

        # {{{ compute work group size

        wg_size = None

        import sys
        import platform
        if ("darwin" in sys.platform
                and "Apple" in queue.device.platform.vendor
                and platform.mac_ver()[0].startswith("10.7")
                and queue.device.type == cl.device_type.CPU):
            wg_size = (1,)

        self.wg_size = wg_size

        # }}}

        self.state = cl_array.empty(queue, (num_work_items, 112), dtype=np.uint8)
        self.state.fill(17)

        prg.init_ranlux(queue, (num_work_items,), self.wg_size, np.uint32(seed),
                self.state.data)

 def _allocate(self, size, dtype, name=None):
     """ Wrapper to define new arrays whether gpu or cpu path"""
     if self.use_gpu:
         import pyopencl.array as cla
         x = cla.empty(self.queue, size, dtype)
         y = arrayWrapper(x, name)
         return y
     else:
         return np.zeros(size, dtype)

def multi_dot(a_gpu, c):
    #a_gpu = cl_array.to_device(queue, a.astype(np.float32))
    c_gpu = cl_array.to_device(queue, c.astype(np.float32))
    out = cl_array.empty(queue, shape=(a_gpu.shape[0], a_gpu.shape[1]), dtype=np.float32)

    prg3.multi_dot(queue, out.shape, (128,1), a_gpu.data, c_gpu.data,
               out.data, np.uint32(a_gpu.shape[-1]), np.uint32(30)).wait()
    ax = out.get()
    return ax

def q2c(X1, X2, X3, queue=None, output=None):
    _check_cl()
    queue = to_queue(queue)
    kern = _q2c_kernel_for_queue(queue.context)

    if X1.shape != X2.shape or X2.shape != X3.shape:
        raise ValueError('All three X matrices must have the same shape.')

    # Create output if not specified
    if output is None:
        output_shape = [1,1,1]
        output_shape[:len(X1.shape[:2])] = X1.shape[:2]
        output_shape[0] >>= 1
        output_shape[1] >>= 1
        output_shape[2] = 6
        output = cl_array.empty(queue, tuple(output_shape), np.complex64)

    # If necessary, convert X
    X1_device = to_device(X1, queue)
    X2_device = to_device(X2, queue)
    X3_device = to_device(X3, queue)

    # Work out size of work group taking into account element step
    work_shape = np.array(output.shape[:3])

    # Work out optimum group size
    if work_shape.shape[0] >= 2 and np.all(work_shape[:2] > 1):
        local_shape = (int(np.floor(np.sqrt(queue.device.max_work_group_size))),) * 2 + (1,1,)
    else:
        local_shape = (queue.device.max_work_group_size, 1, 1)
    local_shape = local_shape[:len(work_shape)]

    global_shape = list(int(np.ceil(x/float(y))*y) for x, y in zip(work_shape, local_shape))

    X_shape = struct.pack('iiii', *(tuple(X1_device.shape) + (1,1,1,1))[:4])

    X1_strides = struct.pack('iiii', *(tuple(s//X1_device.dtype.itemsize for s in X1_device.strides) + (0,0,0,0))[:4])
    X1_offset = np.int32(X1_device.offset)
    X2_strides = struct.pack('iiii', *(tuple(s//X2_device.dtype.itemsize for s in X2_device.strides) + (0,0,0,0))[:4])
    X2_offset = np.int32(X2_device.offset)
    X3_strides = struct.pack('iiii', *(tuple(s//X3_device.dtype.itemsize for s in X3_device.strides) + (0,0,0,0))[:4])
    X3_offset = np.int32(X3_device.offset)

    Y_strides = struct.pack('iiii', *(tuple(s//output.dtype.itemsize for s in output.strides) + (0,0,0,0))[:4])
    Y_shape = struct.pack('iiii', *(tuple(output.shape) + (1,1,1,1))[:4])
    Y_offset = np.int32(output.offset)

    # Perform actual convolution
    kern(queue, global_shape, local_shape,
            X_shape,
            X1_device.base_data, X1_strides, X1_offset,
            X2_device.base_data, X2_strides, X2_offset,
            X3_device.base_data, X3_strides, X3_offset,
            output.base_data, Y_strides, Y_shape, Y_offset)

    return output

    def normal(self, *args, **kwargs):
        """Make a new empty array, apply :meth:`fill_normal` to it.
        """
        mu = kwargs.pop("mu", 0)
        sigma = kwargs.pop("sigma", 1)

        result = cl_array.empty(*args, **kwargs)

        self.fill_normal(result, queue=result.queue, mu=mu, sigma=sigma)
        return result

    def uniform(self, *args, **kwargs):
        """Make a new empty array, apply :meth:`fill_uniform` to it.
        """
        a = kwargs.pop("a", 0)
        b = kwargs.pop("b", 1)

        result = cl_array.empty(*args, **kwargs)

        self.fill_uniform(result, queue=result.queue, a=a, b=b)
        return result

def axis_convolve_ifilter(X, h, axis=0, queue=None, output=None):
    _check_cl()
    queue = to_queue(queue)
    kern = _ifilter_kernel_for_queue(queue.context)

    # Create output if not specified
    if output is None:
        output_shape = list(X.shape)
        output_shape[axis] <<= 1
        output = cl_array.empty(queue, output_shape, np.float32)

    return _apply_kernel(X, h, kern, output, axis=axis, elementstep=0.5)

    def set_ggr(self,ggr):
        
        assert self.can_has_domains, "must set domains before ggr"
        assert isinstance(ggr,tuple) and len(ggr) == 2, "ggr must be a 2-tuple"
        
        growth_rate,ncrossings = ggr
        
        window_length      = 10  # can be changed but not exposed for simplicity      
        rate               = (1+growth_rate)**(1./window_length)-1
        self.plan          = self._ggr_make_plan(self.m0,rate,0.02,50)
        self.target        = 0
        self.optimized_spa = 0.05

        if not self.can_has_ggr:
            self.next_crossing = 0.0
            self.crossed       = False
            self.ggr_tracker   = np.zeros((len(self.plan),3),float)
            self.spa_buffer    = cla.empty(self.queue,(self.N,self.N),np.float32)
            self.whenflipped   = cla.empty(self.queue,(self.N,self.N),np.int32)

            # build the lookup table for the recency enforcement
            # these parameters can be changed but are not exposed to the user to keep things simple
            rmin, rmax, rrate = 0.05, 2., 0.5
            x = np.arange(len(self.plan)).astype('float')
            recency_need = rmin*rmax*np.exp(rrate*x)/(rmax+rmin*np.exp(rrate*x))
            self.recency_need = cla.to_device(self.queue,recency_need.astype(np.float32))
    
            self.set_zero(self.whenflipped)
        
            # self.crossings are the values of m_out which, when crossed over, generate a signal
            # to save the output to make a movie out of or whatever
            if isinstance(ncrossings,(int,float)): self.crossings = np.arange(0,1,1./ncrossings)[1:]
            if isinstance(ncrossings,(list,tuple,np.ndarray)): self.crossings = ncrossings
            if ncrossings != None: self.next_crossing = self.crossings[-1]
        
        self.direction = np.sign(self.m0-self.plan[-1])
        
        self.can_has_ggr = True