Align Python tests to C++/Matlab variants

[janden]

@janden might want to change the tests to align to C++/matlab in a different PR.

Originally posted by @DiamonDinoia in #685 (comment)

[ahbarnett]

I wonder if chatGPT would convert my fullmathtest.m to .py ? Then the CPU+GPU Py tests could be unified, and almost non-random :)

[DiamonDinoia]

chatGPT is down, does this from gemini makes sense?

import numpy as np
import cmath
import sys
import os
import warnings
import math
import time

# Attempt to import finufft (for CPU) and cufinufft (for GPU)
_cpu_backend = None
_gpu_backend = None

try:
    import finufft
    _cpu_backend = finufft
    print("CPU backend 'finufft' found.")
except ImportError:
    warnings.warn("The 'finufft' (CPU) package not found. CPU tests will be skipped.")

try:
    import cufinufft
    # cufinufft automatically brings in cupy and manages GPU context.
    # We'll use cufinufft itself as the 'xp' (array module) when on GPU.
    # Test if a GPU is actually usable with cufinufft.
    # A simple way to check if cufinufft can use a GPU is to try a dummy operation,
    # or rely on its internal checks which might raise an error if no GPU.
    # For now, a successful import is a strong indicator.
    _gpu_backend = cufinufft
    print("GPU backend 'cufinufft' found.")
    # Check if a GPU is truly available and usable
    try:
        # A small test to see if cufinufft can allocate on GPU
        # This will raise a RuntimeError if no GPU or problem with CUDA context
        dummy_array = cufinufft.zeros(1, dtype=cufinufft.float64)
        del dummy_array
        print("cufinufft successfully detected a usable GPU.")
    except Exception as e:
        _gpu_backend = None
        warnings.warn(f"cufinufft found, but GPU not usable: {e}. GPU tests will be skipped.")

except ImportError:
    warnings.warn("The 'cufinufft' (GPU) package not found. GPU tests will be skipped.")

# --- Determine available precision and device combinations ---
test_configs = [] # List of (device_name, prec_str, finufft_module, array_module)

if _cpu_backend:
    test_configs.append(('CPU', 'single', _cpu_backend, np))
    test_configs.append(('CPU', 'double', _cpu_backend, np))

if _gpu_backend:
    test_configs.append(('GPU', 'single', _gpu_backend, _gpu_backend)) # cufinufft is its own array module
    test_configs.append(('GPU', 'double', _gpu_backend, _gpu_backend)) # cufinufft is its own array module

if not test_configs:
    warnings.warn("Found neither 'finufft' (CPU) nor a usable 'cufinufft' (GPU); testing nothing!")
    sys.exit("Exiting: No suitable FINUFFT environment found for testing.")

# --- Main loop for precision and device ---
for dev_name, prec_str, finufft_lib, xp in test_configs:
    # Set precision-specific dtypes
    if prec_str == 'single':
        dtype = xp.float32
        cdtype = xp.complex64
        tol = 1e-3  # requested relative accuracy
    else: # double
        dtype = xp.float64
        cdtype = xp.complex128
        tol = 1e-9

    print(f'{dev_name}\tprec={prec_str}\ttol={tol:.3g}', end='')
    err_check = 10 * tol  # acceptable rel l2 error norm

    # Choose small problems suitable for direct NUDFT computation...
    M = 1000  # # of NU pts (in all dims, and for type 3 targs too)
    Ntot = 1000 # # of modes (approx total, used in all dims)
    ntr = 3   # # transforms
    isign = +1  # sign of imaginary unit in exponential
    # FINUFFT options (some are for internal library use)
    # These options are generally passed as kwargs
    opts = {'debug': 0, 'spread_debug': 0, 'upsampfac': 0}

    print(f', M={M}, Ntot={Ntot}, ntrans={ntr}')

    # random NU pts on whichever device, all dims
    # Data is created directly on the correct device using the chosen 'xp' module
    x = 2 * np.pi * xp.random.rand(M).astype(dtype)
    y = 2 * np.pi * xp.random.rand(M).astype(dtype)
    z = 2 * np.pi * xp.random.rand(M).astype(dtype)

    # complex strengths, possibly stacked in M*ntr array...
    c = (2 * xp.random.rand(M, ntr).astype(dtype) - 1) + \
        1j * (2 * xp.random.rand(M, ntr).astype(dtype) - 1)
    c = c.astype(cdtype)

    # ----------------------------------------------- 1D -----------------
    N = Ntot
    k = xp.arange(math.ceil(-N / 2), math.floor((N - 1) / 2) + 1).astype(dtype)
    # Using the finufft_lib (either finufft or cufinufft) directly
    f = finufft_lib.nufft1d1(x, c, isign, tol, N, isign=isign, dtype=cdtype, eps=tol, **opts)

    A_exact = xp.exp(1j * isign * xp.outer(k, x)).astype(cdtype)
    fe_exact = xp.dot(A_exact, c)
    err = xp.linalg.norm(f.flatten() - fe_exact.flatten()) / xp.linalg.norm(fe_exact.flatten())
    print(f'Rel l2 errs:\t1D type 1:\t{err:.3g}')
    if err > err_check:
        raise RuntimeError(f'{dev_name} {prec_str} 1D type 1 error fail')

    C = finufft_lib.nufft1d2(x, f, isign, tol, isign=isign, dtype=cdtype, eps=tol, **opts)
    Ce_exact = xp.dot(A_exact.T, f)
    err = xp.linalg.norm(C - Ce_exact) / xp.linalg.norm(Ce_exact)
    print(f'\t\t1D type 2:\t{err:.3g}')
    if err > err_check:
        raise RuntimeError(f'{dev_name} {prec_str} 1D type 2 error fail')

    s = N * xp.random.rand(M).astype(dtype)
    f = finufft_lib.nufft1d3(x, c, s, isign, tol, isign=isign, dtype=cdtype, eps=tol, **opts)
    fe_exact = xp.dot(xp.exp(1j * isign * xp.outer(s, x)).astype(cdtype), c)
    err = xp.linalg.norm(f - fe_exact) / xp.linalg.norm(fe_exact)
    print(f'\t\t1D type 3:\t{err:.3g}')
    if err > err_check:
        raise RuntimeError(f'{dev_name} {prec_str} 1D type 3 error fail')

    # ------------------------------------------------- 2D --------------
    N1 = round(math.sqrt(2 * Ntot))
    N2 = round(Ntot / N1)

    kx = xp.arange(math.ceil(-N1 / 2), math.floor((N1 - 1) / 2) + 1).astype(dtype)
    ky = xp.arange(math.ceil(-N2 / 2), math.floor((N2 - 1) / 2) + 1).astype(dtype)
    kx_grid, ky_grid = xp.meshgrid(kx, ky, indexing='ij')
    f = finufft_lib.nufft2d1(x, y, c, isign, tol, (N1, N2), isign=isign, dtype=cdtype, eps=tol, **opts)

    A_exact = xp.exp(1j * isign * (xp.outer(kx_grid.flatten(), x) + \
                                    xp.outer(ky_grid.flatten(), y))).astype(cdtype)
    fe_exact = xp.dot(A_exact, c)
    err = xp.linalg.norm(f.flatten() - fe_exact.flatten()) / xp.linalg.norm(fe_exact.flatten())
    print(f'\t\t2D type 1:\t{err:.3g}')
    if err > err_check:
        raise RuntimeError(f'{dev_name} {prec_str} 2D type 1 error fail')

    C = finufft_lib.nufft2d2(x, y, f, isign, tol, isign=isign, dtype=cdtype, eps=tol, **opts)
    Ce_exact = xp.dot(A_exact.T, f.reshape(N1 * N2, ntr))
    err = xp.linalg.norm(C - Ce_exact) / xp.linalg.norm(Ce_exact)
    print(f'\t\t2D type 2:\t{err:.3g}')
    if err > err_check:
        raise RuntimeError(f'{dev_name} {prec_str} 2D type 2 error fail')

    s = N1 * xp.random.rand(M).astype(dtype)
    t = N2 * xp.random.rand(M).astype(dtype)
    f = finufft_lib.nufft2d3(x, y, c, s, t, isign, tol, isign=isign, dtype=cdtype, eps=tol, **opts)
    fe_exact = xp.dot(xp.exp(1j * isign * (xp.outer(s, x) + xp.outer(t, y))).astype(cdtype), c)
    err = xp.linalg.norm(f - fe_exact) / xp.linalg.norm(fe_exact)
    print(f'\t\t2D type 3:\t{err:.3g}')
    if err > err_check:
        raise RuntimeError(f'{dev_name} {prec_str} 2D type 3 error fail')

    # ------------------------------------------------- 3D --------------
    N1 = round((2 * Ntot)**(1/3))
    N2 = round(Ntot**(1/3))
    N3 = round(Ntot / N1 / N2)

    kx = xp.arange(math.ceil(-N1 / 2), math.floor((N1 - 1) / 2) + 1).astype(dtype)
    ky = xp.arange(math.ceil(-N2 / 2), math.floor((N2 - 1) / 2) + 1).astype(dtype)
    kz = xp.arange(math.ceil(-N3 / 2), math.floor((N3 - 1) / 2) + 1).astype(dtype)
    kx_grid, ky_grid, kz_grid = xp.meshgrid(kx, ky, kz, indexing='ij')
    f = finufft_lib.nufft3d1(x, y, z, c, isign, tol, (N1, N2, N3), isign=isign, dtype=cdtype, eps=tol, **opts)

    A_exact = xp.exp(1j * isign * (xp.outer(kx_grid.flatten(), x) + \
                                    xp.outer(ky_grid.flatten(), y) + \
                                    xp.outer(kz_grid.flatten(), z))).astype(cdtype)
    fe_exact = xp.dot(A_exact, c)
    err = xp.linalg.norm(f.flatten() - fe_exact.flatten()) / xp.linalg.norm(fe_exact.flatten())
    print(f'\t\t3D type 1:\t{err:.3g}')
    if err > err_check:
        raise RuntimeError(f'{dev_name} {prec_str} 3D type 1 error fail')

    C = finufft_lib.nufft3d2(x, y, z, f, isign, tol, isign=isign, dtype=cdtype, eps=tol, **opts)
    Ce_exact = xp.dot(A_exact.T, f.reshape(N1 * N2 * N3, ntr))
    err = xp.linalg.norm(C - Ce_exact) / xp.linalg.norm(Ce_exact)
    print(f'\t\t3D type 2:\t{err:.3g}')
    if err > err_check:
        raise RuntimeError(f'{dev_name} {prec_str} 3D type 2 error fail')

    s = N1 * xp.random.rand(M).astype(dtype)
    t = N2 * xp.random.rand(M).astype(dtype)
    u = N3 * xp.random.rand(M).astype(dtype)
    f = finufft_lib.nufft3d3(x, y, z, c, s, t, u, isign, tol, isign=isign, dtype=cdtype, eps=tol, **opts)
    fe_exact = xp.dot(xp.exp(1j * isign * (xp.outer(s, x) + xp.outer(t, y) + xp.outer(u, z))).astype(cdtype), c)
    err = xp.linalg.norm(f - fe_exact) / xp.linalg.norm(fe_exact)
    print(f'\t\t3D type 3:\t{err:.3g}')
    if err > err_check:
        raise RuntimeError(f'{dev_name} {prec_str} 3D type 3 error fail')

print('\nfullmathtest passed.')

[ahbarnett]

It completely hallucinated this method: dummy_array = cufinufft.zeros(1, dtype=cufinufft.float64)

…

On Tue, Jun 10, 2025 at 5:22 PM Marco Barbone ***@***.***> wrote: *DiamonDinoia* left a comment (flatironinstitute/finufft#695) <#695 (comment)> chatGPT is down, does this from gemini makes sense? import numpy as npimport cmathimport sysimport osimport warningsimport mathimport time # Attempt to import finufft (for CPU) and cufinufft (for GPU)_cpu_backend = None_gpu_backend = None try: import finufft _cpu_backend = finufft print("CPU backend 'finufft' found.")except ImportError: warnings.warn("The 'finufft' (CPU) package not found. CPU tests will be skipped.") try: import cufinufft # cufinufft automatically brings in cupy and manages GPU context. # We'll use cufinufft itself as the 'xp' (array module) when on GPU. # Test if a GPU is actually usable with cufinufft. # A simple way to check if cufinufft can use a GPU is to try a dummy operation, # or rely on its internal checks which might raise an error if no GPU. # For now, a successful import is a strong indicator. _gpu_backend = cufinufft print("GPU backend 'cufinufft' found.") # Check if a GPU is truly available and usable try: # A small test to see if cufinufft can allocate on GPU # This will raise a RuntimeError if no GPU or problem with CUDA context dummy_array = cufinufft.zeros(1, dtype=cufinufft.float64) del dummy_array print("cufinufft successfully detected a usable GPU.") except Exception as e: _gpu_backend = None warnings.warn(f"cufinufft found, but GPU not usable: {e}. GPU tests will be skipped.") except ImportError: warnings.warn("The 'cufinufft' (GPU) package not found. GPU tests will be skipped.") # --- Determine available precision and device combinations ---test_configs = [] # List of (device_name, prec_str, finufft_module, array_module) if _cpu_backend: test_configs.append(('CPU', 'single', _cpu_backend, np)) test_configs.append(('CPU', 'double', _cpu_backend, np)) if _gpu_backend: test_configs.append(('GPU', 'single', _gpu_backend, _gpu_backend)) # cufinufft is its own array module test_configs.append(('GPU', 'double', _gpu_backend, _gpu_backend)) # cufinufft is its own array module if not test_configs: warnings.warn("Found neither 'finufft' (CPU) nor a usable 'cufinufft' (GPU); testing nothing!") sys.exit("Exiting: No suitable FINUFFT environment found for testing.") # --- Main loop for precision and device ---for dev_name, prec_str, finufft_lib, xp in test_configs: # Set precision-specific dtypes if prec_str == 'single': dtype = xp.float32 cdtype = xp.complex64 tol = 1e-3 # requested relative accuracy else: # double dtype = xp.float64 cdtype = xp.complex128 tol = 1e-9 print(f'{dev_name}\tprec={prec_str}\ttol={tol:.3g}', end='') err_check = 10 * tol # acceptable rel l2 error norm # Choose small problems suitable for direct NUDFT computation... M = 1000 # # of NU pts (in all dims, and for type 3 targs too) Ntot = 1000 # # of modes (approx total, used in all dims) ntr = 3 # # transforms isign = +1 # sign of imaginary unit in exponential # FINUFFT options (some are for internal library use) # These options are generally passed as kwargs opts = {'debug': 0, 'spread_debug': 0, 'upsampfac': 0} print(f', M={M}, Ntot={Ntot}, ntrans={ntr}') # random NU pts on whichever device, all dims # Data is created directly on the correct device using the chosen 'xp' module x = 2 * np.pi * xp.random.rand(M).astype(dtype) y = 2 * np.pi * xp.random.rand(M).astype(dtype) z = 2 * np.pi * xp.random.rand(M).astype(dtype) # complex strengths, possibly stacked in M*ntr array... c = (2 * xp.random.rand(M, ntr).astype(dtype) - 1) + \ 1j * (2 * xp.random.rand(M, ntr).astype(dtype) - 1) c = c.astype(cdtype) # ----------------------------------------------- 1D ----------------- N = Ntot k = xp.arange(math.ceil(-N / 2), math.floor((N - 1) / 2) + 1).astype(dtype) # Using the finufft_lib (either finufft or cufinufft) directly f = finufft_lib.nufft1d1(x, c, isign, tol, N, isign=isign, dtype=cdtype, eps=tol, **opts) A_exact = xp.exp(1j * isign * xp.outer(k, x)).astype(cdtype) fe_exact = xp.dot(A_exact, c) err = xp.linalg.norm(f.flatten() - fe_exact.flatten()) / xp.linalg.norm(fe_exact.flatten()) print(f'Rel l2 errs:\t1D type 1:\t{err:.3g}') if err > err_check: raise RuntimeError(f'{dev_name} {prec_str} 1D type 1 error fail') C = finufft_lib.nufft1d2(x, f, isign, tol, isign=isign, dtype=cdtype, eps=tol, **opts) Ce_exact = xp.dot(A_exact.T, f) err = xp.linalg.norm(C - Ce_exact) / xp.linalg.norm(Ce_exact) print(f'\t\t1D type 2:\t{err:.3g}') if err > err_check: raise RuntimeError(f'{dev_name} {prec_str} 1D type 2 error fail') s = N * xp.random.rand(M).astype(dtype) f = finufft_lib.nufft1d3(x, c, s, isign, tol, isign=isign, dtype=cdtype, eps=tol, **opts) fe_exact = xp.dot(xp.exp(1j * isign * xp.outer(s, x)).astype(cdtype), c) err = xp.linalg.norm(f - fe_exact) / xp.linalg.norm(fe_exact) print(f'\t\t1D type 3:\t{err:.3g}') if err > err_check: raise RuntimeError(f'{dev_name} {prec_str} 1D type 3 error fail') # ------------------------------------------------- 2D -------------- N1 = round(math.sqrt(2 * Ntot)) N2 = round(Ntot / N1) kx = xp.arange(math.ceil(-N1 / 2), math.floor((N1 - 1) / 2) + 1).astype(dtype) ky = xp.arange(math.ceil(-N2 / 2), math.floor((N2 - 1) / 2) + 1).astype(dtype) kx_grid, ky_grid = xp.meshgrid(kx, ky, indexing='ij') f = finufft_lib.nufft2d1(x, y, c, isign, tol, (N1, N2), isign=isign, dtype=cdtype, eps=tol, **opts) A_exact = xp.exp(1j * isign * (xp.outer(kx_grid.flatten(), x) + \ xp.outer(ky_grid.flatten(), y))).astype(cdtype) fe_exact = xp.dot(A_exact, c) err = xp.linalg.norm(f.flatten() - fe_exact.flatten()) / xp.linalg.norm(fe_exact.flatten()) print(f'\t\t2D type 1:\t{err:.3g}') if err > err_check: raise RuntimeError(f'{dev_name} {prec_str} 2D type 1 error fail') C = finufft_lib.nufft2d2(x, y, f, isign, tol, isign=isign, dtype=cdtype, eps=tol, **opts) Ce_exact = xp.dot(A_exact.T, f.reshape(N1 * N2, ntr)) err = xp.linalg.norm(C - Ce_exact) / xp.linalg.norm(Ce_exact) print(f'\t\t2D type 2:\t{err:.3g}') if err > err_check: raise RuntimeError(f'{dev_name} {prec_str} 2D type 2 error fail') s = N1 * xp.random.rand(M).astype(dtype) t = N2 * xp.random.rand(M).astype(dtype) f = finufft_lib.nufft2d3(x, y, c, s, t, isign, tol, isign=isign, dtype=cdtype, eps=tol, **opts) fe_exact = xp.dot(xp.exp(1j * isign * (xp.outer(s, x) + xp.outer(t, y))).astype(cdtype), c) err = xp.linalg.norm(f - fe_exact) / xp.linalg.norm(fe_exact) print(f'\t\t2D type 3:\t{err:.3g}') if err > err_check: raise RuntimeError(f'{dev_name} {prec_str} 2D type 3 error fail') # ------------------------------------------------- 3D -------------- N1 = round((2 * Ntot)**(1/3)) N2 = round(Ntot**(1/3)) N3 = round(Ntot / N1 / N2) kx = xp.arange(math.ceil(-N1 / 2), math.floor((N1 - 1) / 2) + 1).astype(dtype) ky = xp.arange(math.ceil(-N2 / 2), math.floor((N2 - 1) / 2) + 1).astype(dtype) kz = xp.arange(math.ceil(-N3 / 2), math.floor((N3 - 1) / 2) + 1).astype(dtype) kx_grid, ky_grid, kz_grid = xp.meshgrid(kx, ky, kz, indexing='ij') f = finufft_lib.nufft3d1(x, y, z, c, isign, tol, (N1, N2, N3), isign=isign, dtype=cdtype, eps=tol, **opts) A_exact = xp.exp(1j * isign * (xp.outer(kx_grid.flatten(), x) + \ xp.outer(ky_grid.flatten(), y) + \ xp.outer(kz_grid.flatten(), z))).astype(cdtype) fe_exact = xp.dot(A_exact, c) err = xp.linalg.norm(f.flatten() - fe_exact.flatten()) / xp.linalg.norm(fe_exact.flatten()) print(f'\t\t3D type 1:\t{err:.3g}') if err > err_check: raise RuntimeError(f'{dev_name} {prec_str} 3D type 1 error fail') C = finufft_lib.nufft3d2(x, y, z, f, isign, tol, isign=isign, dtype=cdtype, eps=tol, **opts) Ce_exact = xp.dot(A_exact.T, f.reshape(N1 * N2 * N3, ntr)) err = xp.linalg.norm(C - Ce_exact) / xp.linalg.norm(Ce_exact) print(f'\t\t3D type 2:\t{err:.3g}') if err > err_check: raise RuntimeError(f'{dev_name} {prec_str} 3D type 2 error fail') s = N1 * xp.random.rand(M).astype(dtype) t = N2 * xp.random.rand(M).astype(dtype) u = N3 * xp.random.rand(M).astype(dtype) f = finufft_lib.nufft3d3(x, y, z, c, s, t, u, isign, tol, isign=isign, dtype=cdtype, eps=tol, **opts) fe_exact = xp.dot(xp.exp(1j * isign * (xp.outer(s, x) + xp.outer(t, y) + xp.outer(u, z))).astype(cdtype), c) err = xp.linalg.norm(f - fe_exact) / xp.linalg.norm(fe_exact) print(f'\t\t3D type 3:\t{err:.3g}') if err > err_check: raise RuntimeError(f'{dev_name} {prec_str} 3D type 3 error fail') print('\nfullmathtest passed.') — Reply to this email directly, view it on GitHub <#695 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACNZRSXNI4N6GE2RPIOUN4D3C5ECVAVCNFSM6AAAAAB7AQGK7SVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDSNRQGYYTGNZXHE> . You are receiving this because you commented.Message ID: ***@***.***>

-- *-------------------------------------------------------------------~^`^~._.~' |\ Alex Barnett Center for Computational Mathematics, Flatiron Institute | \ http://users.flatironinstitute.org/~ahb 646-876-5942

[DiamonDinoia]

I guess we can't replace @janden with gemini yet :)