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.

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 pycuda library. Consequently, we start with a few import statements.

import numpy as np

import pycuda.autoinit
from pycuda.gpuarray import to_gpu

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 = 2 * np.pi * np.random.uniform(size=M)
c = (np.random.standard_normal(size=M)
     + 1J * np.random.standard_normal(size=M))

Now that the data is prepared, we need to set up a cuFINUFFT plan that can be executed on that data.

# create plan
plan = cufinufft.Plan(1, (N,), dtype=np.float64)

# set the nonuniform points
plan.setpts(to_gpu(x))

With everything set up, we are now ready to execute the plan.

# execute the plan
f_gpu = plan.execute(to_gpu(c))

# 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.

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, valid only in [-3pi, 3pi].
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, valid only in [-3pi, 3pi].
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, valid only in [-3pi, 3pi].
y (float[M]) – nonuniform points, valid only in [-3pi, 3pi].
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, valid only in [-3pi, 3pi].
y (float[M]) – nonuniform points, valid only in [-3pi, 3pi].
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, valid only in [-3pi, 3pi].
y (float[M]) – nonuniform points, valid only in [-3pi, 3pi].
z (float[M]) – nonuniform points, valid only in [-3pi, 3pi].
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, valid only in [-3pi, 3pi].
y (float[M]) – nonuniform points, valid only in [-3pi, 3pi].
z (float[M]) – nonuniform points, valid only in [-3pi, 3pi].
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]

Python interface (GPU)¶

Quick-start examples¶

Full documentation¶

Table of Contents

Previous topic

Next topic

This Page