Python interface (GPU)

Quick-start examples

As mentioned in the Python GPU installation instructions, the easiest way to install the Python interface for cuFINUFFT is to run pip install cufinufft.

Note: The interface to cuFINUFFT has changed between versions 1.3 and 2.2. Please see Migration to cuFINUFFT v2.2 for details.

Assuming cuFINUFFT has been installed, we will now consider how to calculate a 1D type 1 transform. To manage the GPU and transfer to and from host and device, we will use the cupy library. Consequently, we start with a few import statements.

import cupy as cp

import cufinufft

We then proceed to setting up a few parameters.

# number of nonuniform points
M = 100000

# grid size
N = 200000

# generate positions for the nonuniform points and the coefficients
x_gpu = 2 * cp.pi * cp.random.uniform(size=M)
c_gpu = (cp.random.standard_normal(size=M)
         + 1J * cp.random.standard_normal(size=M))

Now that the data is prepared, we can call cuFINUFFT to compute the transform from the M source points onto a grid of N modes

f_gpu = cufinufft.nufft1d1(x_gpu, c_gpu, (N,))

# move results off the GPU
f = f_gpu.get()

The last line is needed here to copy the data from the device (GPU) to the host (CPU). Now f is a size-N array containing the NUDFT of c at the points x. This example, along with similar examples for other frameworks, can be found in python/cufinufft/examples/getting_started_*.py.

Other possible calculations are possible by supplying different options during plan creation. See the full API documentation below for more information.

Full documentation

cufinufft.nufft1d1(x, data, n_modes=None, out=None, eps=1e-06, isign=1, **kwargs)

1D type-1 (nonuniform to uniform) complex NUFFT

        M-1
f[k1] = SUM c[j] exp(+/-i k1 x(j))
        j=0

    for -N1/2 <= k1 <= (N1-1)/2
Parameters:
  • x (float[M]) – nonuniform points, in the interval [-pi, pi), values outside will be folded

  • c (complex[M] or complex[n_tr, M]) – source strengths.

  • n_modes (integer or integer tuple of length 1, optional) – number of uniform Fourier modes requested (N1, ). May be even or odd; in either case, modes k1 are integers satisfying -N1/2 <= k1 <= (N1-1)/2. Must be specified if out is not given.

  • out (complex[N1] or complex[n_tr, N1], optional) – output array for Fourier mode values. If n_modes is specifed, the shape must match, otherwise n_modes is inferred from out.

  • eps (float, optional) – precision requested (>1e-16).

  • isign (int, optional) – if non-negative, uses positive sign in exponential, otherwise negative sign.

  • **kwargs (optional) – other options may be specified, see the Plan constructor for details.

Note

The output is written into the out array if supplied.

Returns:

The resulting array.

Return type:

complex[N1] or complex[n_tr, N1]

Example (CuPy):

import cupy as cp
import cufinufft

# number of nonuniform points
M = 100

# the nonuniform points
x = 2 * cp.pi * cp.random.uniform(size=M)

# their complex strengths
c = (cp.random.standard_normal(size=M)
     + 1J * cp.random.standard_normal(size=M))

# desired number of Fourier modes
N1 = 50

# calculate the type-1 NUFFT
f = cufinufft.nufft1d1(x, c, (N1, ))
cufinufft.nufft1d2(x, data, out=None, eps=1e-06, isign=-1, **kwargs)

1D type-2 (uniform to nonuniform) complex NUFFT

c[j] = SUM f[k1] exp(+/-i k1 x(j))
       k1

    for j = 0, ..., M-1, where the sum is over -N1/2 <= k1 <= (N1-1)/2
Parameters:
  • x (float[M]) – nonuniform points, in the interval [-pi, pi), values outside will be folded

  • f (complex[N1] or complex[n_tr, N1]) – Fourier mode coefficients, where N1 may be even or odd. In either case the mode indices k1 satisfy -N1/2 <= k1 <= (N1-1)/2.

  • out (complex[M] or complex[n_tr, M], optional) – output array at targets.

  • eps (float, optional) – precision requested (>1e-16).

  • isign (int, optional) – if non-negative, uses positive sign in exponential, otherwise negative sign.

  • **kwargs (optional) – other options may be specified, see the Plan constructor for details.

Note

The output is written into the out array if supplied.

Returns:

The resulting array.

Return type:

complex[M] or complex[n_tr, M]

Example (CuPy):

import cupy as cp
import cufinufft

# number of nonuniform points
M = 100

# the nonuniform points
x = 2 * cp.pi * cp.random.uniform(size=M)

# number of Fourier modes
N1 = 50

# the Fourier mode coefficients
f = (cp.random.standard_normal(size=(N1, ))
     + 1J * cp.random.standard_normal(size=(N1, )))

# calculate the type-2 NUFFT
c = cufinufft.nufft1d2(x, f)
cufinufft.nufft2d1(x, y, data, n_modes=None, out=None, eps=1e-06, isign=1, **kwargs)

2D type-1 (nonuniform to uniform) complex NUFFT

            M-1
f[k1, k2] = SUM c[j] exp(+/-i (k1 x(j) + k2 y(j)))
            j=0

    for -N1/2 <= k1 <= (N1-1)/2, -N2/2 <= k2 <= (N2-1)/2
Parameters:
  • x (float[M]) – nonuniform points, in the interval [-pi, pi), values outside will be folded

  • y (float[M]) – nonuniform points, in the interval [-pi, pi), values outside will be folded

  • c (complex[M] or complex[n_tr, M]) – source strengths.

  • n_modes (integer or integer tuple of length 2, optional) – number of uniform Fourier modes requested (N1, N2). May be even or odd; in either case, modes k1, k2 are integers satisfying -N1/2 <= k1 <= (N1-1)/2, -N2/2 <= k2 <= (N2-1)/2. Must be specified if out is not given.

  • out (complex[N1, N2] or complex[n_tr, N1, N2], optional) – output array for Fourier mode values. If n_modes is specifed, the shape must match, otherwise n_modes is inferred from out.

  • eps (float, optional) – precision requested (>1e-16).

  • isign (int, optional) – if non-negative, uses positive sign in exponential, otherwise negative sign.

  • **kwargs (optional) – other options may be specified, see the Plan constructor for details.

Note

The output is written into the out array if supplied.

Returns:

The resulting array.

Return type:

complex[N1, N2] or complex[n_tr, N1, N2]

Example (CuPy):

import cupy as cp
import cufinufft

# number of nonuniform points
M = 100

# the nonuniform points
x = 2 * cp.pi * cp.random.uniform(size=M)
y = 2 * cp.pi * cp.random.uniform(size=M)

# their complex strengths
c = (cp.random.standard_normal(size=M)
     + 1J * cp.random.standard_normal(size=M))

# desired number of Fourier modes
N1, N2 = 50, 75

# calculate the type-1 NUFFT
f = cufinufft.nufft2d1(x, y, c, (N1, N2))
cufinufft.nufft2d2(x, y, data, out=None, eps=1e-06, isign=-1, **kwargs)

2D type-2 (uniform to nonuniform) complex NUFFT

c[j] = SUM f[k1, k2] exp(+/-i (k1 x(j) + k2 y(j)))
       k1, k2

    for j = 0, ..., M-1, where the sum is over -N1/2 <= k1 <= (N1-1)/2,
    -N2/2 <= k2 <= (N2-1)/2
Parameters:
  • x (float[M]) – nonuniform points, in the interval [-pi, pi), values outside will be folded

  • y (float[M]) – nonuniform points, in the interval [-pi, pi), values outside will be folded

  • f (complex[N1, N2] or complex[n_tr, N1, N2]) – Fourier mode coefficients, where N1, N2 may be even or odd. In either case the mode indices k1, k2 satisfy -N1/2 <= k1 <= (N1-1)/2, -N2/2 <= k2 <= (N2-1)/2.

  • out (complex[M] or complex[n_tr, M], optional) – output array at targets.

  • eps (float, optional) – precision requested (>1e-16).

  • isign (int, optional) – if non-negative, uses positive sign in exponential, otherwise negative sign.

  • **kwargs (optional) – other options may be specified, see the Plan constructor for details.

Note

The output is written into the out array if supplied.

Returns:

The resulting array.

Return type:

complex[M] or complex[n_tr, M]

Example (CuPy):

import cupy as cp
import cufinufft

# number of nonuniform points
M = 100

# the nonuniform points
x = 2 * cp.pi * cp.random.uniform(size=M)
y = 2 * cp.pi * cp.random.uniform(size=M)

# number of Fourier modes
N1, N2 = 50, 75

# the Fourier mode coefficients
f = (cp.random.standard_normal(size=(N1, N2))
     + 1J * cp.random.standard_normal(size=(N1, N2)))

# calculate the type-2 NUFFT
c = cufinufft.nufft2d2(x, y, f)
cufinufft.nufft3d1(x, y, z, data, n_modes=None, out=None, eps=1e-06, isign=1, **kwargs)

3D type-1 (nonuniform to uniform) complex NUFFT

                M-1
f[k1, k2, k3] = SUM c[j] exp(+/-i (k1 x(j) + k2 y(j) + k3 z(j)))
                j=0

    for -N1/2 <= k1 <= (N1-1)/2, -N2/2 <= k2 <= (N2-1)/2, -N3/2 <= k3 <=
    (N3-1)/2
Parameters:
  • x (float[M]) – nonuniform points, in the interval [-pi, pi), values outside will be folded

  • y (float[M]) – nonuniform points, in the interval [-pi, pi), values outside will be folded

  • z (float[M]) – nonuniform points, in the interval [-pi, pi), values outside will be folded

  • c (complex[M] or complex[n_tr, M]) – source strengths.

  • n_modes (integer or integer tuple of length 3, optional) – number of uniform Fourier modes requested (N1, N2, N3). May be even or odd; in either case, modes k1, k2, k3 are integers satisfying -N1/2 <= k1 <= (N1-1)/2, -N2/2 <= k2 <= (N2-1)/2, -N3/2 <= k3 <= (N3-1)/2. Must be specified if out is not given.

  • out (complex[N1, N2, N3] or complex[n_tr, N1, N2, N3], optional) – output array for Fourier mode values. If n_modes is specifed, the shape must match, otherwise n_modes is inferred from out.

  • eps (float, optional) – precision requested (>1e-16).

  • isign (int, optional) – if non-negative, uses positive sign in exponential, otherwise negative sign.

  • **kwargs (optional) – other options may be specified, see the Plan constructor for details.

Note

The output is written into the out array if supplied.

Returns:

The resulting array.

Return type:

complex[N1, N2, N3] or complex[n_tr, N1, N2, N3]

Example (CuPy):

import cupy as cp
import cufinufft

# number of nonuniform points
M = 100

# the nonuniform points
x = 2 * cp.pi * cp.random.uniform(size=M)
y = 2 * cp.pi * cp.random.uniform(size=M)
z = 2 * cp.pi * cp.random.uniform(size=M)

# their complex strengths
c = (cp.random.standard_normal(size=M)
     + 1J * cp.random.standard_normal(size=M))

# desired number of Fourier modes
N1, N2, N3 = 50, 75, 100

# calculate the type-1 NUFFT
f = cufinufft.nufft3d1(x, y, z, c, (N1, N2, N3))
cufinufft.nufft3d2(x, y, z, data, out=None, eps=1e-06, isign=-1, **kwargs)

3D type-2 (uniform to nonuniform) complex NUFFT

c[j] = SUM f[k1, k2, k3] exp(+/-i (k1 x(j) + k2 y(j) + k3 z(j)))
       k1, k2, k3

    for j = 0, ..., M-1, where the sum is over -N1/2 <= k1 <= (N1-1)/2,
    -N2/2 <= k2 <= (N2-1)/2, -N3/2 <= k3 <= (N3-1)/2
Parameters:
  • x (float[M]) – nonuniform points, in the interval [-pi, pi), values outside will be folded

  • y (float[M]) – nonuniform points, in the interval [-pi, pi), values outside will be folded

  • z (float[M]) – nonuniform points, in the interval [-pi, pi), values outside will be folded

  • f (complex[N1, N2, N3] or complex[n_tr, N1, N2, N3]) – Fourier mode coefficients, where N1, N2, N3 may be even or odd. In either case the mode indices k1, k2, k3 satisfy -N1/2 <= k1 <= (N1-1)/2, -N2/2 <= k2 <= (N2-1)/2, -N3/2 <= k3 <= (N3-1)/2.

  • out (complex[M] or complex[n_tr, M], optional) – output array at targets.

  • eps (float, optional) – precision requested (>1e-16).

  • isign (int, optional) – if non-negative, uses positive sign in exponential, otherwise negative sign.

  • **kwargs (optional) – other options may be specified, see the Plan constructor for details.

Note

The output is written into the out array if supplied.

Returns:

The resulting array.

Return type:

complex[M] or complex[n_tr, M]

Example (CuPy):

import cupy as cp
import cufinufft

# number of nonuniform points
M = 100

# the nonuniform points
x = 2 * cp.pi * cp.random.uniform(size=M)
y = 2 * cp.pi * cp.random.uniform(size=M)
z = 2 * cp.pi * cp.random.uniform(size=M)

# number of Fourier modes
N1, N2, N3 = 50, 75, 100

# the Fourier mode coefficients
f = (cp.random.standard_normal(size=(N1, N2, N3))
     + 1J * cp.random.standard_normal(size=(N1, N2, N3)))

# calculate the type-2 NUFFT
c = cufinufft.nufft3d2(x, y, z, f)
class cufinufft.Plan(nufft_type, n_modes, n_trans=1, eps=1e-06, isign=None, dtype='complex64', **kwargs)

A non-uniform fast Fourier transform (NUFFT) plan

The Plan class lets the user exercise more fine-grained control over the execution of an NUFFT. First, the plan is created with a certain set of parameters (type, mode configuration, precision, sign, number of simultaneous transforms, and so on). Then the nonuniform points are set (source or target depending on the type). Finally, the plan is executed on some data, yielding the desired output.

Parameters:
  • nufft_type (int) – type of NUFFT (1 or 2).

  • n_modes (tuple of ints) – the number of modes in each dimension (for example (50, 100)).

  • n_trans (int, optional) – number of transforms to compute.

  • eps (float, optional) – precision requested (>1e-16).

  • isign (int, optional) – if +1, uses the positive sign exponential, otherwise the negative sign; defaults to +1 for type 1 and to -1 for type 2.

  • dtype (string, optional) – the precision of the transform, ‘complex64’ or ‘complex128’.

  • **kwargs (optional) – additional options corresponding may be specified as keyword arguments. These include upsampling ratio upsampfac (default 2.0), gpu_method (1: nonuniform points-driven, 2: shared memory), gpu_sort (for gpu_method == 1, 0: no sort, 1: sort), gpu_kerevalmeth (0: direct exp(sqrt), Horner evaluation), gpu_device_id (GPU ID), and gpu_stream (CUDA stream pointer).

setpts(x, y=None, z=None, s=None, t=None, u=None)

Set the nonuniform points

For type 1, this sets the coordinates of the M nonuniform source points and for type 2, it sets the coordinates of the M target points.

The dimension of the plan determines the number of arguments supplied. For example, if dim == 2, we provide x and y.

Parameters:
  • x (float[M]) – first coordinate of the nonuniform points (source for type 1, target for type 2).

  • y (float[M], optional) – second coordinate of the nonuniform points (source for type 1, target for type 2).

  • z (float[M], optional) – third coordinate of the nonuniform points (source for type 1, target for type 2).

execute(data, out=None)

Execute the plan

Performs the NUFFT specified at plan instantiation with the points set by setpts. For type-1 and type-3 transforms, the input is a set of source strengths, while for a type-2 transform, it consists of an array of size n_modes. If n_trans is greater than one, n_trans inputs are expected, stacked along the first axis.

Parameters:
  • data (complex[M], complex[n_tr, M], complex[n_modes], or complex[n_tr, n_modes]) – The input source strengths (type 1) or source modes (type 2).

  • out (complex[n_modes], complex[n_tr, n_modes], complex[M], or complex[n_tr, M], optional) – The array where the output is stored. Must be of the correct size.

Returns:

The output array of the transform(s).

Return type:

complex[n_modes], complex[n_tr, n_modes], complex[M], or complex[n_tr, M]