How to use sigpy - 10 common examples

>>> input = np.array([0, 1, 2, 3, 4, 5])
        >>> print(array_to_blocks(input, [2], [2]))
        [[0, 1],
         [2, 3],
         [4, 5]]

    """
    ndim = input.ndim

    if ndim != len(blk_shape) or ndim != len(blk_strides):
        raise ValueError('Input must have same dimensions as blocks.')

    num_blks = [(i - b + s) // s for i, b,
                s in zip(input.shape, blk_shape, blk_strides)]
    device = backend.get_device(input)
    xp = device.xp
    with device:
        output = xp.zeros(num_blks + blk_shape, dtype=input.dtype)

        if ndim == 1:
            if device == backend.cpu_device:
                _array_to_blocks1(output, input,
                                  blk_shape[-1],
                                  blk_strides[-1],
                                  num_blks[-1])
            else:  # pragma: no cover
                _array_to_blocks1_cuda(input,
                                       blk_shape[-1],
                                       blk_strides[-1],
                                       num_blks[-1],
                                       output,

os_shape[a] = os_i
            idx = xp.arange(i, dtype=input.dtype)

            # Calculate apodization
            apod = (beta**2 - (np.pi * width * (idx - i // 2) / os_i)**2)**0.5
            apod /= xp.sinh(apod)

            # Swap axes
            output = output.swapaxes(a, -1)
            os_shape[a], os_shape[-1] = os_shape[-1], os_shape[a]

            # Apodize
            output *= apod

            # Oversampled FFT
            output = util.resize(output, os_shape)
            output = fft.fft(output, axes=[-1], norm=None)
            output /= i**0.5

            # Swap back
            output = output.swapaxes(a, -1)
            os_shape[a], os_shape[-1] = os_shape[-1], os_shape[a]

        coord = _scale_coord(backend.to_device(coord, device), input.shape, oversamp)
        table = backend.to_device(
            _kb(np.arange(n, dtype=coord.dtype) / n, width, beta, dtype=coord.dtype), device)

        output = interp.interp(output, width, table, coord)

        return output

def _apodize(input, ndim, oversamp, width, beta):
    xp = backend.get_array_module(input)
    output = input
    for a in range(-ndim, 0):
        i = output.shape[a]
        os_i = ceil(oversamp * i)
        idx = xp.arange(i, dtype=output.dtype)

        # Calculate apodization
        apod = (beta**2 - (np.pi * width * (idx - i // 2) / os_i)**2)**0.5
        apod /= xp.sinh(apod)
        output *= apod.reshape([i] + [1] * (-a - 1))

    return output

def _scale_coord(coord, shape, oversamp):
    ndim = coord.shape[-1]
    device = backend.get_device(coord)
    scale = backend.to_device([_get_ugly_number(oversamp * i) / i for i in shape[-ndim:]], device)
    shift = backend.to_device([_get_ugly_number(oversamp * i) // 2 for i in shape[-ndim:]], device)

    with device:
        coord = scale * coord + shift

    return coord

def _scale_coord(coord, shape, oversamp):
    ndim = coord.shape[-1]
    device = backend.get_device(coord)
    scale = backend.to_device([_get_ugly_number(oversamp * i) / i for i in shape[-ndim:]], device)
    shift = backend.to_device([_get_ugly_number(oversamp * i) // 2 for i in shape[-ndim:]], device)

    with device:
        coord = scale * coord + shift

    return coord

def fwt(input, wave_name='db4', axes=None, level=None):
    """Forward wavelet transform.

    Args:
        input (array): Input array.
        axes (None or tuple of int): Axes to perform wavelet transform.
        wave_name (str): Wavelet name.
        level (None or int): Number of wavelet levels.
    """
    device = backend.get_device(input)
    input = backend.to_device(input, backend.cpu_device)

    zshape = [((i + 1) // 2) * 2 for i in input.shape]
    zinput = util.resize(input, zshape)

    coeffs = pywt.wavedecn(
        zinput, wave_name, mode='zero', axes=axes, level=level)
    output, _ = pywt.coeffs_to_array(coeffs, axes=axes)

    output = backend.to_device(output, device)
    return output

def _get_params(self):
        self.device = sp.Device(self.device)
        self.dtype = self.y.dtype
        self.data_ndim = self.y.ndim - self.multi_channel - 1
        if self.checkpoint_path is not None:
            self.checkpoint_path = pathlib.Path(self.checkpoint_path)
            self.checkpoint_path.mkdir(parents=True, exist_ok=True)

        self.batch_size = min(len(self.y), self.batch_size)
        self.num_batches = len(self.y) // self.batch_size

        self.L_shape = [self.num_filters] + [self.filt_width] * self.data_ndim
        if self.multi_channel:
            self.L_shape = [self.y.shape[1]] + self.L_shape

        if self.mode == 'full':
            self.R_t_shape = [self.batch_size, self.num_filters] + [i - self.filt_width + 1
                                                                    for i in self.y.shape[-self.data_ndim:]]

output = interp.gridding(input, os_shape, width, table, coord)

        for a in range(-ndim, 0):
            i = oshape[a]
            os_i = os_shape[a]
            idx = xp.arange(i, dtype=input.dtype)
            os_shape[a] = i

            # Swap axes
            output = output.swapaxes(a, -1)
            os_shape[a], os_shape[-1] = os_shape[-1], os_shape[a]

            # Oversampled IFFT
            output = fft.ifft(output, axes=[-1], norm=None)
            output *= os_i / i**0.5
            output = util.resize(output, os_shape)

            # Calculate apodization
            apod = (beta**2 - (np.pi * width * (idx - i // 2) / os_i)**2)**0.5
            apod /= xp.sinh(apod)

            # Apodize
            output *= apod

            # Swap back
            output = output.swapaxes(a, -1)
            os_shape[a], os_shape[-1] = os_shape[-1], os_shape[a]

        return output

def test_stspa_3d_nonexplicit(self):
        nz = 3
        target, sens = self.problem_3d(3, nz)
        dim = target.shape[0]

        g, k1, t, s = rf.spiral_arch(0.24, dim, 4e-6, 200, 0.035)
        k1 = k1 / dim

        k1 = rf.stack_of(k1, nz, 0.1)
        A = sp.mri.linop.Sense(sens, k1, weights=None, tseg=None,
                               ishape=target.shape).H

        pulses = sp.mri.rf.stspa(target, sens, st=None,
                                 coord=k1,
                                 dt=4e-6, max_iter=30, alpha=10, tol=1E-3,
                                 phase_update_interval=200, explicit=False)

        npt.assert_array_almost_equal(A*pulses, target, 1E-3)

def test_MaxEig(self):
        n = 5
        mat = util.randn([n, n])
        A = linop.MatMul([n, 1], mat)
        s = np.linalg.svd(mat, compute_uv=False)

        npt.assert_allclose(app.MaxEig(
            A.H * A, max_iter=1000, show_pbar=False).run(), s[0]**2, atol=1e-3)