Significance of pulsation searches

This page explains how HENzsearch and HENaccelsearch turn the peaks of a periodogram into probabilities, and how the numbers they use were measured. The recipes for using them are in the Pulsation searches and Targeted searches around a known ephemeris tutorials.

Every candidate table reports three probabilities:

p_1trial

The probability that pure noise gives a peak at least this high in a single, specified point of the search.

p_ntrial

The same, corrected for the naive number of trials: one per resolution element of the search.

p_ntrial_adj

The same, corrected for the calibrated number of trials described below. This is the one to use for a blind search.

The correction for \(N\) trials is \(p_N = 1 - (1 - p_1)^N\), which is about \(N p_1\) for small probabilities.

Measuring the number of trials

The points of a search grid are correlated. Two frequencies closer than the resolution \(1/T\) see almost the same noise, but not exactly the same, and the maximum of the search picks up the fluctuations between them. The number of independent trials is therefore not the number of grid points, nor the number of resolution elements: it has to be measured.

We simulate many data sets of pure noise, run the search on each, and record the smallest single-trial probability \(p_{\min}\) of the whole search. If the search behaved as \(N\) independent trials, \(p_{\min}\) would be smaller than \(p\) with probability \(1 - (1 - p)^N\). Calling \(p_\alpha\) the value that \(p_{\min}\) falls below in a fraction \(\alpha\) of the simulations, the effective number of trials at false alarm probability \(\alpha\) is

\[N_{\rm eff}(\alpha) = \frac{\ln(1 - \alpha)}{\ln(1 - p_\alpha)}.\]

The simulations use observations of \(T = 1000\) s with 5000 events of Poisson noise, and the uncertainties come from bootstrap resampling of the simulations. They are run by notebooks/trials_calibration.py, for example:

$ python notebooks/trials_calibration.py qffa --nreal 3000 --outdir calib
$ python notebooks/trials_calibration.py plot --outdir calib

The summary tables used on this page are in docs/images/trials_calibration.

For every configuration, HENDRICS uses the largest estimate at 10% and 1% false alarm probability, rounded up to 0.1, and forced not to decrease with the number of harmonics or the oversampling. Between tabulated values it uses the next larger one. All these choices overestimate the number of trials, which makes the significances conservative.

Independent trials measured on pure noise at 1% false alarm probability. Left: HENzsearch --fast, per resolution element, with the tabulated values as horizontal ticks (frequency only) and crosses (frequency and frequency derivative). Right: HENaccelsearch, per frequency bin per z row, with the tabulated values as black ticks.

The fast folding search (HENzsearch --fast)

The naive number of trials is the number of grid points divided by --oversample for each searched axis. The resolution elements are \(1/T\) in frequency and \(4/T^2\) in frequency derivative (the latter gives the same phase error at the edges of the observation as \(1/T\) in frequency, with phases measured from its middle). The calibrated number of trials is the naive one times:

Trials per resolution element, frequency only

Harmonics	--oversample 1	2	4	8
1	1.1	1.9	3.0	3.1
2	1.1	2.1	4.7	8.0
4	1.1	2.1	5.0	8.0

Trials per resolution element, frequency and frequency derivative

Harmonics	--oversample 1	2	4
1	0.8	3.7	6.1
2	0.9	3.7	11.0
4	1.0	4.3	13.9

Grids beyond the table use its largest value, with a warning: their significances may be overestimated. On the same grids, the slow folding search (HENzsearch without --fast) gives the same numbers within the uncertainties (summary_folding.ecsv).

The tables were measured on bands of 0.25 Hz (frequency only) and 0.05 Hz (with the frequency derivative), with 32 phase bins and 3000 simulations each. Wider bands (1 and 4 Hz) and 128 phase bins give the same numbers within the uncertainties, or smaller ones. A run of 20000 simulations of the \(Z^2_2\) search with --oversample 4 gives 3.4, 3.7 and 2.9 trials per resolution element at 10%, 1% and 0.1% false alarm probability, below the tabulated 4.7: the calibration holds in the tail.

Why oversample, then? Because a signal between two grid points loses power. The figure below shows the power recovered by the search relative to the power at the exact frequency, for a \(Z^2_2\) signal: 80% without oversampling, 95% with --oversample 2, and all of it from --oversample 4. The price is in trials; the default of 2 * N points per resolution element for \(Z^2_N\) is a compromise between the two.

Power recovered by the search, relative to the power at the true frequency, averaged over signal frequencies. The horizontal lines are the accelerated search, with and without interbinning.

The accelerated search (HENaccelsearch)

Stingray searches \(n_{\rm freq}\) Fourier bins for each of the \(n_z\) rows of np.arange(-zmax, zmax, delta_z), where \(z\) is the drift in Fourier bins over the observation. The calibrated number of trials is \(n_{\rm freq} \times \max(n_z, 1) \times c\), with \(c\) from:

Trials per frequency bin per z row

Search	Without --interbin	With --interbin
No z search (zmax 0)	1.0	2.1
delta_z 1	0.9	1.1
delta_z 0.5	0.7	0.9

Steps between 0.5 and 1 are interpolated linearly, finer steps use the 0.5 value, and coarser steps count the rows as independent. For the default search (zmax 100, delta_z 1) this is 180 trials per frequency bin, which raises the detection threshold by about 10 in Leahy power compared with the old count of one trial per frequency bin. The numbers were measured with 10000 simulations, with zmax 0, 10 and 100.

HENaccelsearch delegates the search to Stingray, which only returns the candidates above a power threshold computed with \(n_{\rm freq}\) trials. HENDRICS asks instead for the threshold corresponding to the requested false alarm probability with the calibrated number of trials (lowered, if needed, so that it also catches the in-between bins of interbinning, see below), and then keeps the candidates with p_ntrial_adj below the false alarm probability. With a known ephemeris it asks for the threshold of a single trial, and keeps the candidates by their corrected probability.

--pad_to_double was not simulated: its significances are only approximate.

Interbinning

Interbinning (--interbin) recovers part of the power lost between Fourier bins by adding a bin between each pair,

\[A_{k+1/2} = \frac{\pi}{4} \left(A_{k+1} - A_k\right).\]

The regular bins keep the usual noise distribution of Leahy powers, a \(\chi^2\) with two degrees of freedom: \(p_1 = e^{-P/2}\). The in-between bins do not. They are the difference of two neighbouring bins, which after the acceleration correction are correlated. Their noise is still exponential, but stretched by

\[s(z) = \frac{\pi^2}{8}\left(1 - {\rm Re}\,\rho_1(z)\right), \qquad \rho_1(z) = \frac{\sum_q w_{q+1} w_q^*}{\sum_q |w_q|^2},\]

where \(w\) is the response used by Stingray to correct for the drift \(z\), and \(\rho_1\) the correlation between neighbouring bins that it introduces. The single-trial probability of an in-between bin is \(p_1 = e^{-P/(2 s)}\). For a plain Fourier search (\(z = 0\)) the bins are independent and \(s = \pi^2/8 \approx 1.23\): without this correction, a noise power \(P\) in an in-between bin looked as significant as a power \(1.23\,P\) really is.

HENDRICS recognizes in-between bins from their frequency (Stingray reports \(r/T\), so \(r\) is a half integer), and computes \(s(z)\) from Stingray’s responses (hendrics.known_ephemeris.interbin_stretch). The tests check, on simulated white noise passed through Stingray’s own convolution and interbinning, that the resulting probabilities are uniformly distributed for \(z\) = 0, 0.25, 0.5, 1, 5, 10 and 100.

Stretch of the noise distribution of the in-between bins as a function of z. HENaccelsearch uses integer z values with the default --delta-z.

Targeted searches

With a known ephemeris (see Targeted searches around a known ephemeris), each candidate is charged a number of trials that grows with its distance from the prediction. The cells of the search are ranked, before looking at the data, by that distance. The charge of a candidate, ntrial_eff, is the calibrated number of trials of the blind search times the fraction of the grid that is at least as close to the prediction, with a floor ntrial_min covering the uncertainty region of the ephemeris (three standard deviations), and never more than the blind search. Its probability is p_value.

Picking the best candidate

p_value is honest for a candidate chosen before looking at the data. But the reported candidate is the most significant one after the correction, and that choice has a price. Think of the false alarm probability \(\alpha\) as a budget: a cell charged \(n\) trials uses up about \(\alpha / n\) of it, and the chance that any cell gives a false alarm is at most the sum over all cells. For a blind search, \(N\) cells charged \(N\) trials each use exactly \(\alpha\). Ranking the cells charges the \(r\)-th closest cell about \(r\) trials, and the sum becomes \(\alpha \sum_r 1/r \approx \alpha \left(1 + \ln(N / n_{\min})\right)\), where the first \(n_{\min}\) cells, inside the floor, contribute \(\alpha\) together.

p_value_best pays for this, charging \(\min\left(k \times {\tt ntrial\_eff},\ N\right)\) trials with \(k = 1 + \ln(N / n_{\min})\). \(k\) is about 10 for a typical search and no uncertainty on the ephemeris, and 1 when the uncertainty region covers the whole search. HENzsearch uses p_value_best to decide whether the candidate is detected, and HENaccelsearch sorts its candidates by it.

The table below shows the rate of false alarms measured on 10000 simulations of pure noise, relative to the nominal false alarm probability (1 is exact, less is conservative), with the calibrated number of trials of each tool: 2 harmonics, --oversample 4 and a frequency derivative search for HENzsearch --fast (5632 trials), zmax 10 for HENaccelsearch (18000 trials).

False alarms relative to nominal, at 10%, 1% and 0.1%

Tool	Floor ntrial_min	p_value	p_value_best	Blind search
HENzsearch --fast	1	3.5, 5.8, 7.0	1.15, 1.33, 1.7	0.79, 0.85, 1.2
HENzsearch --fast	100	2.2, 3.3, 3.7	1.06, 1.24, 1.6	0.79, 0.85, 1.2
HENaccelsearch	1	4.0, 5.6, 5.2	1.16, 1.34, 1.7	0.79, 1.0, 1.1
HENaccelsearch	100	2.3, 2.8, 2.7	1.05, 1.22, 1.5	0.79, 1.0, 1.1

The figures at 0.1% rest on about ten events each, and are uncertain by \(\pm 0.4\).

p_value_best removes most of the excess, but not all of it: both tools still give false alarms about 1.1 times too often at 10%, and 1.2-1.3 times at 1%. The reason is the promise that a candidate never costs more than in a blind search. Most cells of the grid are far enough from the prediction that \(k \times {\tt ntrial\_eff}\) would exceed \(N\), so they are charged exactly \(N\): together, they use up almost the whole budget, as a blind search would. The cells near the prediction, charged less, add their share on top. For the HENaccelsearch grid above, the sum of the charges allows up to 1.8 times the nominal rate of false alarms for independent cells (1.6 with a floor of 100). Splitting the budget in advance between the cells near the prediction and a blind search would remove this, at the price of charging distant candidates several times more than a blind search. We chose to keep the promise and document the residual excess.

False alarm rate of the targeted searches relative to nominal, at 10% (left) and 1% (right) false alarm probability, as a function of the floor on the trials, with and without the price of picking the best candidate.

GPU acceleration

HENDRICS can compute some of its histograms on an NVIDIA GPU (Graphics Processing Unit) through CuPy. This is optional and off by default: HENDRICS installs, imports and passes its tests on machines without a GPU.

Installation

pip install "hendrics[gpu]" installs the cupy package, which is compiled from source and needs the CUDA toolkit. It is usually easier to install the pre-built CuPy wheel matching your CUDA version, e.g. pip install cupy-cuda12x or pip install "cupy-cuda13x[ctk]", where [ctk] also installs the CUDA libraries (see the CuPy installation guide).

What runs on the GPU

The GPU is used only when explicitly requested. From the command line, the --use-gpu option of HENzsearch affects:

• --fast: the whole search runs on the GPU. The event times are copied to the GPU once. Then, at each step of the frequency and frequency-derivative grid, the pulse phases, the sub-profiles (the 2D histogram of pulse phase vs. time) and the Z^2 statistics of all trial shifts are computed on the GPU, and only the Z^2 values are copied back. The results are identical to the CPU ones, to the last bit (see “Design” below).

• --transient: the 2D histogram (pulse phase vs. time), at each trial frequency.

• --ffa: the 1D histogram of the whole event list (computed once).

HENefsearch accepts --use-gpu with --transient and --ffa too (--fast is only available for Z searches).

The standard folding search (without --fast, --ffa or --transient) folds the events with Stingray on the CPU, so --use-gpu has no effect there and a warning is issued.

From Python, hendrics.base.histogram() and hendrics.base.histogram2d() accept use_gpu=True, and so do hendrics.efsearch.search_with_qffa(), hendrics.efsearch.transient_search() and hendrics.efsearch.search_with_ffa(). The GPU histograms are also available directly as hendrics.gpu.histogram_gpu() and hendrics.gpu.histogram2d_gpu().

The GPU histograms give the same results as the CPU ones:

• Values equal to the upper edge of the range are excluded, as in the CPU histograms of HENDRICS. numpy.histogram() and cupy.histogram include them in the last bin instead, so the GPU functions remove these values before binning. This matters in practice: in the --fast search the upper edge of the time range is the time of the last photon.

• The output types are the same: 64-bit floats for 1D histograms and for weighted 2D histograms, 64-bit unsigned integers for unweighted 2D histograms.

• The use_memmap and tmp options (memory-mapped output files) only apply to host memory, and are ignored on the GPU.

Each call to a GPU histogram copies the input to the GPU and the histogram back to host memory. The --transient search does this at each trial frequency.

Design

hendrics.gpu contains a minimal backend registry, _BACKENDS. Each entry says whether the backend can be used, which array module implements it (NumPy or CuPy), and how to copy arrays to its memory and back to host memory. Only "cpu" and "cupy" exist now; other array libraries (e.g. JAX or PyTorch) could be added without changing the functions using the registry.

The --fast search on the GPU. At every step of the search, the event times and the time slice (sub-profile) of each event are the same; only the pulse phases change. hendrics.efsearch._FastSearchOnDevice is created once per search: it copies the times and the slice indices to the GPU. At each step, its profiles method computes the phases and the sub-profiles on the GPU (with bincount on precomputed indices), and its stats method shifts and sums the sub-profiles for every trial and computes Z^2, copying back only the Z^2 values.

The second part is a small CUDA program (a “kernel”), in hendrics/gpu.py, running one GPU thread per trial shift. It is compiled at run time by CuPy (cupy.RawKernel), so it needs no dependency beyond CuPy. Writing it with numba.cuda would have needed the separate numba-cuda package.

The GPU results are identical to those of the Numba code on the CPU, not only close:

• The phases are computed with the same operations, in the same order, as the Numba functions _fast_phase*, and the time slice and phase bin of each event with the same arithmetic as the Numba 2D histogram. The events that the Numba histogram drops (the last one, on the upper edge of the time range, and those whose phase rounds to exactly 1.0) are dropped on the GPU too.

• The shifted sub-profiles are sums of whole-number counts, which are exact in any order.

• The Z^2 sums follow the same order as the Numba code, and the kernel is compiled with --fmad=false: without this, the compiler could merge a multiplication and an addition into a single operation (a “fused multiply-add”), which rounds differently in the last bit.

• The shift ramps and the cosine and sine tables are computed once by Numba (_fast_step_constants, also used by the CPU _fast_step) and copied to the GPU: Numba’s linspace does not round exactly as NumPy’s.

When the array module is NumPy, _FastSearchOnDevice falls back to the Numba _fast_step, so the same code runs on the CPU.

Without Numba. Numba is an optional dependency. Without it, the CPU histograms of HENDRICS (hendrics.base.histogram() and hendrics.base.histogram2d()) use NumPy with the same bin arithmetic, the same treatment of the upper edge and the same output types as the Numba versions, so histograms and sub-profiles are identical with or without Numba, and on the GPU. _fast_step then runs as plain Python, where NumPy adds the Z^2 terms in a different order: Z^2 values can differ from the Numba (and GPU) ones in the last digits.

The summed profile of each trial is a small array local to each GPU thread, whose size is fixed when the kernel is compiled (one compiled kernel per number of bins). With a shared buffer in GPU memory instead, the kernel was about 4 times slower, because the compiler could not assume that nothing else writes to it.

GPU memory: at most about 49 bytes per event (466 MiB for 1e7 events, with 16 or 128 phase bins): the times and slice indices (12 bytes per event) for the whole search, plus the temporary arrays of each step. It is released at the end of the search, and CuPy keeps it in its memory pool for reuse. The first time Numba compiles _fast_step_constants in a Python session, the compiler keeps a reference to the running functions until Python’s garbage collector runs, so the memory of that first search is released a little later.

import cupy succeeds even on machines without a GPU, so the presence of a CUDA device is checked only when the GPU is requested. If CuPy or a device is missing, a RuntimeError explains what to install. The choice between CPU and GPU is never automatic.

Testing without a GPU

• hendrics/tests/test_gpu.py replaces the CuPy entry of the backend registry with a stand-in that uses NumPy and counts the copies to the GPU and back to host memory. The GPU code is mostly a thin layer on top of the array module, so this checks argument handling, the edge convention, the output types, the command line options and the number of copies: in the --fast search, the events are copied to the GPU once per search, and one array comes back per step. The CUDA kernel of _fast_step cannot run on NumPy: with the stand-in, the Numba _fast_step is used instead.

• The tests using the real CuPy are skipped without a GPU. With a GPU, they check that sub-profiles, Z^2 values and whole --fast searches are identical to the CPU ones (np.array_equal, not a tolerance), for all phase formulas, several numbers of bins and harmonics, and with and without the frequency derivative search. The exact comparisons of Z^2 values with the real CuPy also need Numba.

• hendrics/tests/test_base.py checks that the NumPy fallback histograms follow the conventions of the Numba ones, and are identical to them when Numba is installed. The CI environments without the all extra run without Numba.

• hendrics/tests/test_gpu_kernel_sim.py runs the benchmark code (below) with numba’s CUDA simulator, which executes numba.cuda kernels on the CPU. The simulator is enabled with NUMBA_ENABLE_CUDASIM=1 before numba.cuda is imported, so the checks run in a separate Python process. They compare the custom kernel with the CPU histogram (including the case that made the original prototype read past the end of the input array), and the averaged power spectrum loops with Stingray, using NumPy in place of CuPy.

Benchmarks

benchmarks/gpu_histogram_benchmark.py is run by hand on a machine with a GPU:

python benchmarks/gpu_histogram_benchmark.py --sizes 100000 10000000

It reports the best time of each method, its speedup and whether its result agrees with the CPU version, for:

• 1D histograms: HENDRICS Numba (CPU), a custom numba.cuda kernel, and cupy.histogram;

• a 1D histogram followed by an FFT, all on the GPU;

• 2D histograms with the shapes used by the --fast search: HENDRICS Numba (CPU) and cupy.histogram2d;

• averaged power spectra from events: Stingray (CPU); an “unfused” GPU loop, that copies each light curve and each FFT back to host memory, as would happen by replacing the histogram and the FFT used by Stingray with GPU versions; and a “fused” GPU loop, where binning, FFT and averaging all stay on the GPU with a single copy at the end;

• one step of the --fast search (--benchmarks qffa), split into its parts: phases, 2D histogram and _fast_step on the CPU; copying the events to the GPU (once per search), sub-profiles and _fast_step on the GPU.

The custom numba.cuda kernel comes from the original GPU prototype (PR #181). It is kept only in the benchmark: HENDRICS uses cupy.histogram and cupy.histogram2d, which are tested upstream and also support weights and 2D histograms. On the hardware below, the custom kernel was 1.5 to 2.3 times as fast as cupy.histogram, but a 1D histogram alone beat the CPU only with 1e7 events.

On the same hardware, a GPU averaged power spectrum that keeps all data on the GPU was 3.3 to 4.7 times faster than one copying each light curve and each FFT back to host memory. This loop is not yet used by HENfspec.

Results on real hardware

Measured on 2026-09-15 on an NVIDIA GeForce RTX 2070 SUPER (8 GB, driver 595.91.07, CUDA 13.2) and an Intel Core i7-10700 (8 cores, 16 threads), with Python 3.14.4, NumPy 2.4.3, Numba 0.66.0 with numba-cuda 0.30.4, CuPy 14.2.0 (cupy-cuda13x) and Stingray 2.3.2. CuPy needed the CUDA libraries from pip install "cupy-cuda13x[ctk]", and the custom kernel needed pip install "numba-cuda[cu13]". All tests passed, including those using the real CuPy, and every benchmark result agreed with the CPU.

Speedup with respect to the CPU (CPU time divided by GPU time, so values below 1 mean that the GPU is slower). Times include the copies to and from the GPU.

Benchmark (number of events)	1e5	1e6	1e7
1D histogram, 1e6 bins: numba.cuda kernel	0.15	0.74	1.49
1D histogram, 1e6 bins: cupy.histogram	0.10	0.32	0.65
1D histogram + FFT: numba.cuda kernel	8.9	6.5	2.7
1D histogram + FFT: cupy.histogram	5.4	3.0	1.4
2D histogram, 16x256 bins: cupy.histogram2d	0.18	0.88	1.42
2D histogram, 16x4096 bins: cupy.histogram2d	0.21	1.06	1.72
Averaged PDS, 64 x 524288 bins: unfused	2.9	2.7	2.5
Averaged PDS, 64 x 524288 bins: fused	13.0	12.7	8.2
Averaged PDS, 256 x 524288 bins: unfused	2.9	2.7	2.8
Averaged PDS, 256 x 524288 bins: fused	13.7	12.3	12.3

The Stingray averaged power spectrum took 1.2 s with 64 segments and 4.5-4.9 s with 256 segments, almost independently of the number of events.

The --fast search on the GPU

Measured on 2026-09-16, same hardware and software. The Z^2 values and the candidates were identical on CPU and GPU in all cases.

End to end, HENzsearch --fast -f 1.22 -F 1.25 -N 2 --find-candidates -n NBIN on a simulated event file with 1e7 events over 1e5 s and a 2% pulsation at 1.235 Hz (376 steps of the search, 32 x 32 trials each). 16 bins is the minimum for N = 2 (8 bins per harmonic); 128 is the default. Total times include reading the file and analysing the candidates; the search loop times come from the progress bar. These times were measured before the Numba compilation cache was added, so they include compiling the Numba functions, which happens during the first step of the search: about 5 s on the CPU and 1.4 s on the GPU (see “Numba compilation cache” below).

Phase bins	CPU: total (search loop)	GPU: total (search loop)	Speedup
16	27.0 s (23 s)	7.2 s (2 s)	3.8 (11)
128	64.0 s (60 s)	16.2 s (11 s)	4.0 (5.5)

One step of the search (--benchmarks qffa), 1e7 events, 1024 trials, best of 5 calls, in milliseconds:

Part of the step	16 bins	128 bins
CPU: phases	10.8	11.0
CPU: 2D histogram	28.6	40.9
CPU: _fast_step (shift, sum, Z^2)	1.9	75.7
CPU: whole step	43.8	131.6
GPU: copy the events to the GPU (once per search)	65.3	66.8
GPU: phases and sub-profiles	3.2	4.2
GPU: _fast_step and copy back	0.8	27.6
GPU: whole step	4.0	27.6

With few bins the histogram dominates the CPU step, and with many bins _fast_step does; the GPU speeds up both.

For comparison, before these changes (only the 2D histogram on the GPU, copying the event arrays at each step) the same search with 128 bins took 79 s on the CPU and 69 s with --use-gpu in the search loop. An earlier measurement with a 16 x 16 grid of trials per step (before PR #203 made --oversample count grid points per 1/T) gave 45 s and 35 s.

CPU improvements

Profiling the --fast search while preparing the GPU version showed two costs that did not need a GPU, now removed (1e7 events, 128 bins, 1024 trials):

• The Numba 2D histograms stacked their two input arrays into a new array at each call (160 MB with 1e7 events). Passing them separately made the 2D histogram go from 62.3 to 42.6 ms.

• _fast_step read the sub-profiles through a transposed view of the histogram, with the values of each sub-profile scattered in memory. Copying them so that each sub-profile is contiguous in memory made the whole step go from 161.8 to 130.0 ms, with identical results. With 16 bins, _fast_step is too short for this to matter.

Overall, one CPU step with 128 bins went from 180 ms (measured on 2026-09-15) to 131 ms.

Numba compilation cache

Numba compiles a function the first time it is called in a Python session. Each run of a command-line script like HENzsearch is a new Python session, so without a cache the same functions are compiled again at every run. The time this takes does not depend on the number of events: compiling _fast_step takes about 5 s, the largest part of the fixed cost of a --fast search.

All the njit functions of HENDRICS are compiled with cache=True: the first time, Numba saves the compiled code to disk, and later Python sessions load it instead of compiling again. The code loaded from the cache is the same, so the results do not change.

• Where. The cache files (.nbi and .nbc) are written to the __pycache__ folders next to the HENDRICS source files. If these folders cannot be written (for example, in a system-wide installation), Numba uses a folder in the user’s cache directory (~/.cache/numba on Linux). The NUMBA_CACHE_DIR environment variable selects another folder.

• When Numba compiles again. The first run after installing or upgrading HENDRICS still compiles, and so does the first call with new types of arguments. Numba also compiles again when the source file of a function changes, or with a different version of Numba or Python.

• Caveat for developers. Numba only checks the source file of the function it loads from the cache, not the files of the functions it calls. For example, _fast_step in efsearch.py calls _z_n_fast_cached in ffa.py: after editing only ffa.py, the old compiled _fast_step would still be loaded. After such edits, delete the cache files, e.g. with find hendrics -name "*.nb[ic]" -delete.

• @vectorize functions with a list of signatures are compiled when their module is imported, and caching does not help them: for a test copy of sum_arrays, caching changed that time from 0.36 s to 0.30 s, because most of it is not spent in the part that Numba can cache.

  • hendrics.ffa.sum_arrays, used by the FFA search, was such a function, and made import hendrics.efsearch (and so HENzsearch and HENefsearch) compile it even without an FFA search. It is now a cached njit function, which Numba only compiles when it is first called. The import went from 3.04 s to 2.83 s (median of 10 runs), ffa_search on 1e6 bins from 1.87 s to 1.79 s, and the FFA results did not change (the same values and types, for 64- and 32-bit floats and integers). hendrics/tests/test_numba_cache.py checks that importing hendrics.efsearch compiles no Numba code.

  • phases_from_zero_to_one and phases_around_zero in hendrics.ml_timing are still compiled at import (0.25 s), but only HENphaseogram imports that module. Compiling them only when first called would change their output types in some cases (e.g. 32-bit input giving 64-bit output), so they were left as they are.

hendrics/tests/test_numba_cache.py checks that every njit function of HENDRICS has caching turned on, and runs a small --fast search in two separate Python processes sharing an empty cache folder: in the second process, the functions of the search are loaded from the cache, and the Z^2 values are identical to those of the first process.

Measured on 2026-09-16 on the hardware and software described in “Results on real hardware” above, with HENzsearch --fast -f 1.22 -F 1.25 -N 2 --find-candidates -n NBIN on the simulated event files of that section (1e5 or 1e7 events over 1e5 s). Each command was run three times, starting from an empty cache folder: the first run compiles and fills the cache, the second and third load from it. Times are in seconds, for the whole command. All periodograms were identical.

Events	Bins	Device	Before (mean)	After: run 1	After: 2 and 3
1e5	16	CPU	10.6	11.0	4.8, 5.0
1e5	16	GPU	5.8	5.8	4.5, 4.5
1e5	128	CPU	44.8	45.8	35.2, 35.4
1e5	128	GPU	14.6	14.6	13.3, 13.3
1e7	16	CPU	27.2	28.1	22.0, 22.0
1e7	16	GPU	7.4	7.2	6.0, 6.0
1e7	128	CPU	67.4	67.8	58.3, 57.5
1e7	128	GPU	16.3	16.3	15.2, 15.0

From the second run on, each run is about 1.3 s faster on the GPU, where only _fast_step_constants is compiled, and 5 to 6 s faster on the CPU with 16 bins. With 128 bins the CPU saves about 10 s, more than the 5.8 s of compilation measured with Numba’s compilation events: in these runs, the search itself was also faster with the code loaded from the cache (32.8 s, against 34.8 to 37.5 s with the code compiled in the same process, 1e5 events). We have not found the reason.

The savings do not depend on the data. The cache depends on the source code and on the types of the arguments (e.g. arrays of 64-bit floats), not on their values. After filling the cache with the 1e5-event file above, we searched three new simulated event files: 1e6 to 5e6 events, lengths from 5e4 to 2e5 s, pulsations at 0.7 to 2.1 Hz, 1 to 3 harmonics and 16 to 64 bins. On every file, the new code ran 5.5 to 6.0 s faster on the CPU and 1.1 to 1.3 s faster on the GPU. It spent 0.08 s loading the cache, while the old code spent 5.7 to 5.9 s (CPU) and 1.4 s (GPU) compiling. The periodograms were identical. The first run with an option that needs a function not yet in the cache compiles that function once: with --mean-fdot, _fast_phase_fdot took 0.4 s the first time and was loaded from the cache afterwards.