def find_local_size(global_size, flat_local_size, threshold=0.05):
    """
    Returns a tuple of the same size as ``global_size``,
    with the product equal to ``flat_local_size``,
    and minimal difference between ``product(global_size)``
    and ``product(min_blocks(gs, ls) for gs, ls in zip(global_size, local_size))``
    (i.e. tries to minimize the amount of empty threads).
    """
    flat_global_size = product(global_size)
    if flat_local_size >= flat_global_size:
        return global_size

    threads_num = product(global_size)

    best_ratio = None
    best_local_size = None

    for local_size in get_decompositions(flat_local_size, len(global_size)):
        bounding_global_size = tuple(
            ls * min_blocks(gs, ls) for gs, ls in zip(global_size, local_size))
        empty_threads = product(bounding_global_size) - threads_num
        ratio = float(empty_threads) / threads_num

        # Stopping iteration early, because there may be a lot of elements to iterate over,
        # and we do not need the perfect solution.
        if ratio < threshold:
            return local_size

        if best_ratio is None or ratio < best_ratio:
            best_ratio = ratio
            best_local_size = local_size

    return best_local_size

def check_performance(thr_and_double, shape_and_axes):

    thr, double = thr_and_double

    dtype = numpy.complex128 if double else numpy.complex64
    dtype = dtypes.normalize_type(dtype)

    shape, axes = shape_and_axes

    data = numpy.arange(product(shape)).reshape(shape).astype(dtype)

    shift = FFTShift(data, axes=axes)
    shiftc = shift.compile(thr)

    data_dev = thr.to_device(data)
    res_dev = thr.empty_like(data)

    attempts = 10
    times = []
    for i in range(attempts):
        t1 = time.time()
        shiftc(res_dev, data_dev)
        thr.synchronize()
        times.append(time.time() - t1)

    res_ref = numpy.fft.fftshift(data, axes=axes)
    assert diff_is_negligible(res_dev.get(), res_ref)

    return min(times), product(shape) * dtype.itemsize

    def __init__(self, virtual_shape, available_shape):
        self.real_dims = {}
        self.real_strides = {}
        self.virtual_strides = {}
        self.major_vdims = {}
        self.bounding_shape = tuple()
        self.skip_thresholds = []

        v_groups, a_groups = group_dimensions(virtual_shape, available_shape)

        for v_group, a_group in zip(v_groups, a_groups):
            virtual_subshape = virtual_shape[v_group[0]:v_group[-1]+1]
            virtual_subsize = product(virtual_subshape)

            bounding_subshape = find_bounding_shape(
                virtual_subsize,
                available_shape[a_group[0]:a_group[-1]+1])

            self.bounding_shape += bounding_subshape

            if virtual_subsize < product(bounding_subshape):
                strides = [(adim, product(bounding_subshape[:i])) for i, adim in enumerate(a_group)]
                self.skip_thresholds.append((virtual_subsize, strides))

            for vdim in v_group:
                self.real_dims[vdim] = a_group
                self.real_strides[vdim] = tuple(
                    product(self.bounding_shape[a_group[0]:adim]) for adim in a_group)
                self.virtual_strides[vdim] = product(virtual_shape[v_group[0]:vdim])
                self.major_vdims[vdim] = v_group[-1]

def test_group_dimensions(virtual_shape, available_shape):
    """
    Tests that ``group_dimensions()`` obeys its contracts.
    """
    v_groups, a_groups = vsize.group_dimensions(virtual_shape, available_shape)
    v_dims = []
    a_dims = []
    for v_group, a_group in zip(v_groups, a_groups):

        # Check that axis indices in groups are actually in range
        assert any(vdim < len(virtual_shape) for vdim in v_group)
        assert any(adim < len(available_shape) for adim in a_group)

        # Check that the total number of elements (threads) in the virtual group
        # is not greater than the number of elements in the real group
        v_shape = virtual_shape[v_group[0]:v_group[-1]+1]
        a_shape = available_shape[a_group[0]:a_group[-1]+1]
        assert(product(v_shape) <= product(a_shape))

        v_dims += v_group
        a_dims += a_group

    # Check that both virtual and real groups axes add up to a successive list
    # without intersections.
    assert v_dims == list(range(len(virtual_shape)))
    assert a_dims == list(range(len(available_shape[:len(a_dims)])))

    def _add_transpose(self, plan, device_params,
            mem_out, mem_in, batch_shape, height_shape, width_shape):

        bso = self._block_width_override
        block_width = device_params.local_mem_banks if bso is None else bso

        if block_width ** 2 > device_params.max_work_group_size:
            # If it is not CPU, current solution may affect performance
            block_width = int(numpy.sqrt(device_params.max_work_group_size))

        input_height = helpers.product(height_shape)
        input_width = helpers.product(width_shape)
        batch = helpers.product(batch_shape)

        blocks_per_matrix = helpers.min_blocks(input_height, block_width)
        grid_width = helpers.min_blocks(input_width, block_width)

        render_kwds = dict(
            input_width=input_width, input_height=input_height, batch=batch,
            block_width=block_width,
            grid_width=grid_width,
            blocks_per_matrix=blocks_per_matrix,
            input_slices=[len(batch_shape), len(height_shape), len(width_shape)],
            output_slices=[len(batch_shape), len(width_shape), len(height_shape)])

        plan.kernel_call(
            TEMPLATE.get_def('transpose'), [mem_out, mem_in],
            global_size=(batch, blocks_per_matrix * block_width, grid_width * block_width),
            local_size=(1, block_width, block_width),
            render_kwds=render_kwds)

    def __init__(self, dtype, device_params, outer_shape, fft_size, curr_size,
            fft_size_real, inner_shape, pass_num, reverse_direction):

        num_passes = len(get_global_radix_info(fft_size)[0])
        real_output_shape = (pass_num == num_passes - 1 and reverse_direction)

        self.name = 'fft_global'
        self.inplace_possible = (pass_num == num_passes - 1 and num_passes % 2 == 1)
        self.output_shape = (outer_shape +
            (fft_size_real if real_output_shape else fft_size,) + inner_shape)
        if fft_size != fft_size_real and pass_num == 0 and reverse_direction:
            self.kweights = get_kweights(fft_size_real, fft_size)
        else:
            self.kweights = None

        self._fft_size = fft_size
        self._curr_size = curr_size
        self._fft_size_real = fft_size_real
        self._local_mem_size = device_params.local_mem_size
        self._itemsize = dtype.itemsize
        self._inner_batch = helpers.product(inner_shape)
        self._outer_batch = helpers.product(outer_shape)
        self._pass_num = pass_num
        self._last_pass = (pass_num == num_passes - 1)

        self._constant_kwds = get_common_kwds(dtype, device_params)
        self._constant_kwds.update(dict(
            takes_kweights=(self.kweights is not None),
            input_slices=(len(outer_shape), 1, len(inner_shape)),
            output_slices=(len(outer_shape), 1, len(inner_shape)),
            pad_in=(fft_size != fft_size_real and pass_num == 0 and not reverse_direction),
            unpad_out=(fft_size != fft_size_real and self._last_pass and reverse_direction),
            reverse_direction=reverse_direction,
            normalize=self._last_pass))

    def __call__(self, psi):
        psi_nobs = psi.get()
        psi = beam_splitter(psi_nobs)

        ns_nobs = numpy.abs(psi_nobs) ** 2 - 0.5 / product(self.ds)
        ns = numpy.abs(psi) ** 2 - 0.5 / product(self.ds)

        if len(psi.shape) == 3:
            # 1D
            n_nobs = ns_nobs.mean(1)
            n = ns.mean(1)
            Ns = (ns * self.ds[0]).sum(-1)
        elif len(psi.shape) == 4:
            # 2D
            n_nobs = ns_nobs.mean(1)
            n = ns.mean(1)
            Ns = (ns * product(self.ds)).sum(-1).sum(-1)
        elif len(psi.shape) == 5:
            # 3D
            n_nobs = (ns_nobs.mean(1) * self.ds[1] * self.ds[2]).sum(-1).sum(-1)
            n = (ns.mean(1) * self.ds[1] * self.ds[2]).sum(-1).sum(-1)
            Ns = (ns * product(self.ds)).sum(-1).sum(-1).sum(-1)

        res = dict(
            Nplus_mean=Ns[0].mean(), Nminus_mean=Ns[1].mean(),
            Nplus_std=Ns[0].std(), Nminus_std=Ns[1].std(),
            density=n,
            density_nobs=n_nobs)

        if len(psi.shape) == 5:
            res['slice_nobs'] = ns_nobs.mean(1)[:,:,:,ns.shape[-1] / 2]
            res['slice'] = ns.mean(1)[:,:,:,ns.shape[-1] / 2]

        return res

def group_dimensions(virtual_shape, available_shape):
    """
    Returns two lists, one of tuples with numbers of grouped virtual dimensions, the other
    one of tuples with numbers of corresponding group of available dimensions,
    such that for any group of virtual dimensions, the total number of elements they cover
    does not exceed the number of elements covered by the
    corresponding group of available dimensions.
    """
    assert product(virtual_shape) <= product(available_shape)
    return _group_dimensions(0, virtual_shape, 0, available_shape)

 def __init__(self, dtype, shape=None, strides=None):
     self.shape = tuple() if shape is None else wrap_in_tuple(shape)
     self.size = product(self.shape)
     self.dtype = dtypes.normalize_type(dtype)
     self.ctype = dtypes.ctype_module(self.dtype)
     if strides is None:
         self.strides = tuple([
             self.dtype.itemsize * product(self.shape[i+1:]) for i in range(len(self.shape))])
     else:
         self.strides = strides
     self._cast = dtypes.cast(self.dtype)

    def _build_plan(self, plan_factory, device_params, output, matrix_a, matrix_b):
        bwo = self._block_width_override

        if bwo is not None:
            block_widths = [bwo]
        else:
            nbanks = device_params.local_mem_banks
            block_widths = [2 ** n for n in range(helpers.log2(nbanks), -1, -1)]

        a_batch = helpers.product(matrix_a.shape[:-2])
        b_batch = helpers.product(matrix_b.shape[:-2])
        batch = max(a_batch, b_batch)

        for block_width in block_widths:

            plan = plan_factory()

            if block_width ** 2 > device_params.max_work_group_size:
                continue

            num_steps = helpers.min_blocks(self._convolution_size, block_width)
            a_blocks = helpers.min_blocks(self._a_outer_size, block_width)
            b_blocks = helpers.min_blocks(self._b_outer_size, block_width)

            render_kwds = dict(
                batched_a=(a_batch != 1),
                batched_b=(b_batch != 1),
                transposed_a=self._transposed_a,
                transposed_b=self._transposed_b,
                num_steps=num_steps,
                a_slices=(len(matrix_a.shape) - 2, 1, 1),
                b_slices=(len(matrix_b.shape) - 2, 1, 1),
                output_slices=(len(output.shape) - 2, 1, 1),
                block_width=block_width,
                mul=functions.mul(matrix_a.dtype, matrix_b.dtype, out_dtype=output.dtype))

            try:
                plan.kernel_call(
                    TEMPLATE.get_def('matrixmul'),
                    [output, matrix_a, matrix_b],
                    kernel_name="kernel_matrixmul",
                    global_size=(
                        batch,
                        a_blocks * block_width,
                        b_blocks * block_width),
                    local_size=(1, block_width, block_width),
                    render_kwds=render_kwds)
            except OutOfResourcesError:
                continue

            return plan

        raise ValueError("Could not find suitable call parameters for the kernel")

    def compatible_with(self, other):
        if self.dtype != other.dtype:
            return False

        common_shape_len = min(len(self.shape), len(other.shape))
        if self.shape[-common_shape_len:] != other.shape[-common_shape_len:]:
            return False
        if self.strides[-common_shape_len:] != other.strides[-common_shape_len:]:
            return False
        if helpers.product(self.shape[:-common_shape_len]) != 1:
            return False
        if helpers.product(other.shape[:-common_shape_len]) != 1:
            return False

        return True

def test_find_local_size(global_size, flat_local_size, expected_local_size):
    """
    Checking that ``find_local_size`` finds the sizes we expect from it.
    """
    local_size = vsize.find_local_size(global_size, flat_local_size)
    assert product(local_size) == flat_local_size
    assert local_size == expected_local_size

    def _build_plan(self, plan_factory, device_params, output, input_, inverse):

        if helpers.product([input_.shape[i] for i in self._axes]) == 1:
            return self._build_trivial_plan(plan_factory, output, input_)

        # While resource consumption of GlobalFFTKernel can be made lower by passing
        # lower value to prepare_for(), LocalFFTKernel may have to be split into several kernels.
        # Therefore, if GlobalFFTKernel.prepare_for() raises OutOfResourcesError,
        # we just call prepare_for() with lower limit, but if LocalFFTKernel.prepare_for()
        # does that, we have to recreate the whole chain.
        local_kernel_limit = device_params.max_work_group_size

        while local_kernel_limit >= 1:
            try:
                plan = self._build_limited_plan(
                    plan_factory, device_params, local_kernel_limit, output, input_, inverse)
            except LocalKernelFail:
            # One of LocalFFTKernels was out of resources.
            # Reduce the limit and try to create operations from scratch again.
                local_kernel_limit //= 2
                continue
            except GlobalKernelFail:
                raise ValueError(
                    "Could not find suitable call parameters for one of the global kernels")

            return plan

        raise ValueError("Could not find suitable call parameters for one of the local kernels")

def find_bounding_shape(virtual_size, available_shape):
    """
    Finds a tuple of the same length as ``available_shape``, with every element
    not greater than the corresponding element of ``available_shape``,
    and product not lower than ``virtual_size``.
    """
    assert virtual_size <= product(available_shape)

    free_size = virtual_size
    free_dims = set(range(len(available_shape)))
    bounding_shape = [None] * len(available_shape)

    while len(free_dims) > 0:
        guess = ceiling_root(free_size, len(free_dims))
        for fdim in free_dims:
            bounding_shape[fdim] = guess

        for fdim in free_dims:
            if bounding_shape[fdim] > available_shape[fdim]:
                bounding_shape[fdim] = available_shape[fdim]
                free_dims.remove(fdim)
                free_size = min_blocks(free_size, bounding_shape[fdim])
                break
        else:
            return tuple(bounding_shape)

    return tuple(available_shape)

    def _build_plan(self, plan_factory, device_params, output, input_, inverse):

        if helpers.product([input_.shape[i] for i in self._axes]) == 1:
            return self._build_trivial_plan(plan_factory, output, input_)

        plan = plan_factory()

        axes = tuple(sorted(self._axes))
        shape = list(input_.shape)

        if all(shape[axis] % 2 == 0 for axis in axes):
        # If all shift axes have even length, it is possible to perform the shift inplace
        # (by swapping pairs of elements).
        # Note that the inplace fftshift is its own inverse.
            shape[axes[0]] //= 2
            plan.kernel_call(
                TEMPLATE.get_def('fftshift_inplace'), [output, input_],
                kernel_name="kernel_fftshift_inplace",
                global_size=shape,
                render_kwds=dict(axes=axes))
        else:
        # Resort to an out-of-place shift to a temporary array and then copy.
            temp = plan.temp_array_like(output)
            plan.kernel_call(
                TEMPLATE.get_def('fftshift_outplace'), [temp, input_, inverse],
                kernel_name="kernel_fftshift_outplace",
                global_size=shape,
                render_kwds=dict(axes=axes))

            copy_trf = copy(input_, out_arr_t=output)
            copy_comp = PureParallel.from_trf(copy_trf, copy_trf.input)
            plan.computation_call(copy_comp, output, temp)

        return plan

def check_performance(thr_and_double, shape_and_axes, fast_math):
    thr, double = thr_and_double

    shape, axes = shape_and_axes
    dtype = numpy.complex128 if double else numpy.complex64

    data = get_test_array(shape, dtype)
    data_dev = thr.to_device(data)
    res_dev = thr.empty_like(data_dev)

    fft = FFT(data_dev, axes=axes)
    fftc = fft.compile(thr, fast_math=fast_math)

    attempts = 10
    t1 = time.time()
    for i in range(attempts):
        fftc(res_dev, data_dev)
    thr.synchronize()
    t2 = time.time()
    dev_time = (t2 - t1) / attempts

    fwd_ref = numpy.fft.fftn(data, axes=axes).astype(dtype)
    assert diff_is_negligible(res_dev.get(), fwd_ref)

    return dev_time, product(shape) * sum([numpy.log2(shape[a]) for a in axes]) * 5

def test_find_bounding_shape(virtual_size, available_shape):
    """
    Tests that ``find_bounding_shape()`` obeys its contracts.
    """
    shape = vsize.find_bounding_shape(virtual_size, available_shape)
    assert all(isinstance(d, int) for d in shape)
    assert product(shape) >= virtual_size
    assert all(d <= ad for d, ad in zip(shape, available_shape))

    def prepare_for(self, max_local_size):
        kwds = dict(self._constant_kwds)

        radix_arr, radix1_arr, radix2_arr = get_global_radix_info(self._fft_size)

        radix = radix_arr[self._pass_num]
        radix1 = radix1_arr[self._pass_num]
        radix2 = radix2_arr[self._pass_num]

        stride_out = self._inner_batch * helpers.product(radix_arr[:self._pass_num])
        stride = stride_out * radix
        stride_in = stride_out * helpers.product(radix_arr[self._pass_num+1:])

        threads_per_xform = radix2

        coalesce_width = kwds['min_mem_coalesce_width']
        local_batch = max_local_size if radix2 == 1 else coalesce_width
        local_batch = min(local_batch, stride_in)
        local_size = min(local_batch * threads_per_xform, max_local_size)
        local_batch = local_size // threads_per_xform

        workgroups_num = helpers.min_blocks(stride_in, local_batch) * self._outer_batch

        if radix2 == 1:
            lmem_size = 0
        else:
            if stride_out == 1:
                lmem_size = (radix + 1) * local_batch
            else:
                lmem_size = local_size * radix1

        if lmem_size * self._itemsize // 2 > self._local_mem_size:
            raise OutOfResourcesError

        kwds.update(self._constant_kwds)
        kwds.update(dict(
            fft_size=self._fft_size, curr_size=self._curr_size, fft_size_real=self._fft_size_real,
            pass_num=self._pass_num,
            lmem_size=lmem_size, local_batch=local_batch, local_size=local_size,
            inner_batch=self._inner_batch,
            radix_arr=radix_arr, radix1_arr=radix1_arr, radix2_arr=radix2_arr,
            radix1=radix1, radix2=radix2, radix=radix,
            stride_in=stride_in, stride_out=stride_out, stride=stride,
            last_pass=self._last_pass))

        return workgroups_num * local_size, local_size, kwds

def test_group_dimensions(virtual_shape, available_shape):
    """
    Tests that ``group_dimensions()`` obeys its contracts.
    """
    v_groups, a_groups = vsize.group_dimensions(virtual_shape, available_shape)
    v_dims = []
    a_dims = []
    for v_group, a_group in zip(v_groups, a_groups):
        v_shape = virtual_shape[v_group[0]:v_group[-1]+1]
        a_shape = available_shape[a_group[0]:a_group[-1]+1]
        assert(product(v_shape) <= product(a_shape))

        v_dims += v_group
        a_dims += a_group

    assert v_dims == list(range(len(virtual_shape)))
    assert a_dims == list(range(len(available_shape[:len(a_dims)])))

def test_large_scan_performance(thr, large_perf_shape, exclusive):
    """
    Large problem sizes.
    """
    dtype = dtypes.normalize_type(numpy.int64)
    min_time = check_scan(
        thr, large_perf_shape, dtype=dtype, axes=None, exclusive=exclusive, measure_time=True)
    return min_time, helpers.product(large_perf_shape) * dtype.itemsize