Reference implementation of multiharmonic least-squares periodogram (AoVMH(W) / FastChi²) using blocked NuFFT and a fast Toeplitz solver for complex coefficients.
The codebase is currently in a pre-release state (versioned as 0.0.0), so breaking changes may occur and no release builds are available yet. As such compilation from source is recommended.
To install the Python wrapper:
git clone https://github.com/kkrutkowski/Toeplitz-ls.git --depth=1cd ./Toeplitz_lsmake native # use 'make generic' if build time is a priority — this may take a few minutescd ./toeplitz-lspip3 install .The toeplitz-ls package exposes three submodules, each targeting a different precision level. The single-precision (tlsf) and double-precision (tls) submodules accept and return NumPy arrays. The double-double module (tlsdd) relies on compensated arithmetic and operates on a 106-bit mantissa; it uses Mpmath's mpf type to bridge the gap between the library's internal representation and Python's native numeric types. The get_peaks helper can turn a periodogram into a sorted, interpolated peak list for printing or saving.
Free-threaded (GIL-free) Python support
toeplitz-ls is compatible with free-threaded CPython (PEP 703, available as an optional build since Python 3.13). With the GIL disabled, ThreadPoolExecutor-based batch processing achieves true parallelism without the overhead of multiprocessing — no serialisation, no inter-process memory duplication.
Note that Mpmath is not thread-safe by default; when using tlsdd in a multithreaded context, each worker should maintain its own mp context via mpmath.workprec() rather than mutating the global mp.prec.
When the system matrix at a given frequency is numerically unstable, the solver returns NaN rather than raising an exception. Use np.nanargmax() and np.nanmax() in place of their standard counterparts to safely recover the peak frequency and power without short-circuiting on degenerate bins.
For manual filtering, call power() or autopower() with autonan=False to receive rough conditionality estimates alongside the periodogram. These estimates are solver-specific: LDLT reports a diagonal-ratio estimate, while the three Toeplitz solvers (levinson, zohar, and bareiss) return a reflection-coefficient-based upper bound. The Toeplitz estimates are generally much higher (and more reliable) for the same frequencies, so threshold values should be chosen with the solver family in mind.
importnumpyasnpimportmatplotlib.pyplotaspltfromtoeplitz_lsimporttlsrng=np.random.default_rng(seed=123)
N=1000t=np.sort(rng.uniform(0, 1000, size=N))
y=np.sin(50*t) +0.5*np.sin(150*t) +1+rng.poisson(size=N)
# Compute the generalized Scargle periodogram.freq, r2=tls.autopower(
t, y, dy=None,
fmax=12, #default valuenterms=1, # generalized Scargle periodogramoversampling=5, #frequency density, relative to a dense, uniform FFTbackend="pswf43", # or "pswf21", "lra"normalization="standard"# or "asymptotic"
)
best_idx=np.nanargmax(r2)
best_freq=freq[best_idx]
best_power=r2[best_idx]
print(f"Peak frequency: {best_freq:.6g} (R^2 = {best_power:.6g})")
plt.figure(figsize=(12, 5))
plt.plot(freq, r2, label='degree = 1', color='red')
plt.axvline(best_freq, linestyle='--', color='black', alpha=0.7)
plt.xlabel('Frequency')
plt.ylabel('Coefficient of determination')
plt.legend()
plt.tight_layout()
plt.show()frompathlibimportPathfromtempfileimportTemporaryDirectoryfromurllib.requestimporturlretrieveimportnumpyasnpfromtoeplitz_lsimporttlsf, get_peaksurl="https://www.astrouw.edu.pl/ogle/ogle4/OCVS/BLAP/phot/phot_ogle4/I/OGLE-BLAP-035.dat"withTemporaryDirectory() astmpdir:
data_path=Path(tmpdir) /"OGLE-BLAP-035.dat"urlretrieve(url, str(data_path))
t, y, _dy=np.loadtxt(data_path, unpack=True)
t=t-t.min()
# BLAPs require high-frequency search. Ignore the photometric error estimates here.freq, power, cond=tlsf.autopower(
t, y, dy=None,
fmax=100,
nterms=3,
oversampling=5,
solver="levinson", # or 'zohar', 'bareiss', 'ldlt', 'svd'normalization="asymptotic",
autonan=False,
)
peaks=get_peaks(freq, power, cond, threshold=10)
peaks.print(10)