|
import time |
|
import numpy as np |
|
import ducc0 |
|
|
|
def bench_fft(): |
|
res = {} |
|
# 1D FFT |
|
a1 = (np.random.randn(1000000) + 1j * np.random.randn(1000000)).astype(np.complex128) |
|
for nthreads in [1, 4]: |
|
t0 = time.perf_counter() |
|
iters = 50 |
|
for _ in range(iters): |
|
ducc0.fft.c2c(a1, nthreads=nthreads) |
|
el = (time.perf_counter() - t0) / iters |
|
res[f"fft_1d_1M_threads_{nthreads}"] = el |
|
|
|
# 2D FFT |
|
a2 = (np.random.randn(1024, 1024) + 1j * np.random.randn(1024, 1024)).astype(np.complex128) |
|
for nthreads in [1, 4]: |
|
t0 = time.perf_counter() |
|
iters = 30 |
|
for _ in range(iters): |
|
ducc0.fft.c2c(a2, axes=(0, 1), nthreads=nthreads) |
|
el = (time.perf_counter() - t0) / iters |
|
res[f"fft_2d_1024x1024_threads_{nthreads}"] = el |
|
|
|
# 3D FFT |
|
a3 = (np.random.randn(128, 128, 128) + 1j * np.random.randn(128, 128, 128)).astype(np.complex128) |
|
for nthreads in [1, 4]: |
|
t0 = time.perf_counter() |
|
iters = 30 |
|
for _ in range(iters): |
|
ducc0.fft.c2c(a3, axes=(0, 1, 2), nthreads=nthreads) |
|
el = (time.perf_counter() - t0) / iters |
|
res[f"fft_3d_128x128x128_threads_{nthreads}"] = el |
|
|
|
return res |
|
|
|
def bench_healpix(): |
|
res = {} |
|
nside = 512 |
|
npix = 12 * nside * nside |
|
vlen = 100000 |
|
|
|
rng = np.random.default_rng(42) |
|
pix = rng.integers(low=0, high=npix - 1, size=vlen, dtype=np.int64) |
|
ptg = np.empty((vlen, 2), dtype=np.float64) |
|
ptg[:, 0] = np.arccos((rng.random(vlen) - 0.5) * 2) |
|
ptg[:, 1] = rng.random(vlen) * 2 * np.pi |
|
|
|
base = ducc0.healpix.Healpix_Base(nside, "RING") |
|
|
|
# pix2ang ring |
|
t0 = time.perf_counter() |
|
iters = 20 |
|
for _ in range(iters): |
|
base.pix2ang(pix) |
|
res["healpix_pix2ang_ring"] = (time.perf_counter() - t0) / iters |
|
|
|
# ang2pix ring |
|
t0 = time.perf_counter() |
|
iters = 20 |
|
for _ in range(iters): |
|
base.ang2pix(ptg) |
|
res["healpix_ang2pix_ring"] = (time.perf_counter() - t0) / iters |
|
|
|
# ring2nest |
|
t0 = time.perf_counter() |
|
iters = 20 |
|
for _ in range(iters): |
|
base.ring2nest(pix) |
|
res["healpix_ring2nest"] = (time.perf_counter() - t0) / iters |
|
|
|
# nest2ring |
|
t0 = time.perf_counter() |
|
iters = 20 |
|
for _ in range(iters): |
|
base.nest2ring(pix) |
|
res["healpix_nest2ring"] = (time.perf_counter() - t0) / iters |
|
|
|
return res |
|
|
|
def bench_sht(): |
|
res = {} |
|
lmax = 512 |
|
mmax = 512 |
|
nlat = lmax + 1 |
|
nlon = 2 * lmax + 1 |
|
|
|
nalm = ((mmax + 1) * (mmax + 2)) // 2 + (mmax + 1) * (lmax - mmax) |
|
rng = np.random.default_rng(42) |
|
alm = (rng.uniform(-1.0, 1.0, nalm) + 1j * rng.uniform(-1.0, 1.0, nalm)).astype(np.complex128) |
|
alm[0 : lmax + 1].imag = 0.0 |
|
alm = alm.reshape((1, -1)) |
|
|
|
for nthreads in [1, 4]: |
|
t0 = time.perf_counter() |
|
iters = 5 |
|
for _ in range(iters): |
|
ducc0.sht.synthesis_2d( |
|
alm=alm, |
|
ntheta=nlat, |
|
nphi=nlon, |
|
lmax=lmax, |
|
mmax=mmax, |
|
spin=0, |
|
geometry="GL", |
|
nthreads=nthreads, |
|
) |
|
el = (time.perf_counter() - t0) / iters |
|
res[f"sht_synthesis_2d_gl_lmax512_threads_{nthreads}"] = el |
|
|
|
return res |
|
|
|
def bench_nufft(): |
|
res = {} |
|
npoints = 500000 |
|
ndim = 2 |
|
shape = (512, 512) |
|
rng = np.random.default_rng(42) |
|
|
|
coord = (2 * np.pi * rng.uniform(size=(npoints, ndim)) - np.pi).astype(np.float64) |
|
points = (rng.uniform(size=npoints) - 0.5 + 1j * (rng.uniform(size=npoints) - 0.5)).astype(np.complex128) |
|
|
|
for nthreads in [1, 4]: |
|
out = np.empty(shape, dtype=np.complex128) |
|
t0 = time.perf_counter() |
|
iters = 10 |
|
for _ in range(iters): |
|
ducc0.nufft.nu2u(points=points, coord=coord, forward=True, epsilon=1e-5, nthreads=nthreads, out=out) |
|
el = (time.perf_counter() - t0) / iters |
|
res[f"nufft_type1_2d_500k_threads_{nthreads}"] = el |
|
|
|
return res |
|
|
|
def bench_wgridder(): |
|
res = {} |
|
rng = np.random.default_rng(42) |
|
nrow = 100000 |
|
nchan = 4 |
|
uvw = (rng.random((nrow, 3)) - 0.5) * 5000.0 |
|
freq = np.linspace(1e9, 1.4e9, nchan) |
|
vis = (rng.standard_normal((nrow, nchan)) + 1j * rng.standard_normal((nrow, nchan))).astype(np.complex128) |
|
wgt = rng.random((nrow, nchan)).astype(np.float64) |
|
npix_x = 512 |
|
npix_y = 512 |
|
pixsize_x = 1e-4 |
|
pixsize_y = 1e-4 |
|
|
|
for nthreads in [1, 4]: |
|
t0 = time.perf_counter() |
|
iters = 5 |
|
for _ in range(iters): |
|
ducc0.wgridder.vis2dirty( |
|
uvw=uvw, |
|
freq=freq, |
|
vis=vis, |
|
wgt=wgt, |
|
npix_x=npix_x, |
|
npix_y=npix_y, |
|
pixsize_x=pixsize_x, |
|
pixsize_y=pixsize_y, |
|
epsilon=1e-5, |
|
nthreads=nthreads, |
|
do_wgridding=True |
|
) |
|
el = (time.perf_counter() - t0) / iters |
|
res[f"wgridder_vis2dirty_threads_{nthreads}"] = el |
|
|
|
return res |
|
|
|
def bench_totalconvolve(): |
|
res = {} |
|
lmax = 256 |
|
kmax = 13 |
|
ncomp = 1 |
|
separate = True |
|
nptg = 500000 |
|
epsilon = 1e-4 |
|
|
|
rng = np.random.default_rng(42) |
|
def nalm_fn(l, m): |
|
return ((m+1)*(m+2))//2 + (m+1)*(l-m) |
|
|
|
slm = (rng.uniform(-1., 1., (ncomp, nalm_fn(lmax, lmax))) + 1j*rng.uniform(-1., 1., (ncomp, nalm_fn(lmax, lmax)))).astype(np.complex128) |
|
slm[:, 0:lmax+1].imag = 0. |
|
blm = (rng.uniform(-1., 1., (ncomp, nalm_fn(lmax, kmax))) + 1j*rng.uniform(-1., 1., (ncomp, nalm_fn(lmax, kmax)))).astype(np.complex128) |
|
blm[:, 0:lmax+1].imag = 0. |
|
|
|
ptg = rng.uniform(0., 1., (nptg, 3)).astype(np.float64) |
|
ptg[:, 0] *= np.pi |
|
ptg[:, 1] *= 2*np.pi |
|
ptg[:, 2] *= 2*np.pi |
|
|
|
for nthreads in [1, 4]: |
|
t0 = time.perf_counter() |
|
iters = 5 |
|
for _ in range(iters): |
|
interp = ducc0.totalconvolve.Interpolator(slm, blm, separate, lmax, kmax, epsilon=epsilon, npoints=nptg, nthreads=nthreads) |
|
_ = interp.interpol(ptg) |
|
el = (time.perf_counter() - t0) / iters |
|
res[f"totalconvolve_interp_lmax256_500k_threads_{nthreads}"] = el |
|
|
|
return res |
|
|
|
def run_all_benchmarks(): |
|
print(f"ducc0 module: {ducc0.__file__}") |
|
print(f"ducc0 version: {getattr(ducc0, '__version__', 'unknown')}") |
|
all_res = {} |
|
print("\n--- Running FFT Benchmark ---") |
|
r = bench_fft() |
|
all_res.update(r) |
|
for k, v in r.items(): |
|
print(f" {k}: {v*1000:.3f} ms") |
|
|
|
print("\n--- Running Healpix Benchmark ---") |
|
r = bench_healpix() |
|
all_res.update(r) |
|
for k, v in r.items(): |
|
print(f" {k}: {v*1000:.3f} ms") |
|
|
|
print("\n--- Running SHT Benchmark ---") |
|
r = bench_sht() |
|
all_res.update(r) |
|
for k, v in r.items(): |
|
print(f" {k}: {v*1000:.3f} ms") |
|
|
|
print("\n--- Running NUFFT Benchmark ---") |
|
r = bench_nufft() |
|
all_res.update(r) |
|
for k, v in r.items(): |
|
print(f" {k}: {v*1000:.3f} ms") |
|
|
|
print("\n--- Running WGridder Benchmark ---") |
|
r = bench_wgridder() |
|
all_res.update(r) |
|
for k, v in r.items(): |
|
print(f" {k}: {v*1000:.3f} ms") |
|
|
|
print("\n--- Running TotalConvolve Benchmark ---") |
|
r = bench_totalconvolve() |
|
all_res.update(r) |
|
for k, v in r.items(): |
|
print(f" {k}: {v*1000:.3f} ms") |
|
|
|
return all_res |
|
|
|
if __name__ == "__main__": |
|
import json |
|
import sys |
|
results = run_all_benchmarks() |
|
if len(sys.argv) > 1: |
|
with open(sys.argv[1], "w") as f: |
|
json.dump(results, f, indent=2) |
|
print(f"\nResults saved to {sys.argv[1]}") |