AcyclePy - Cyclostratigraphy and Time Series Analysis Toolkit

AcyclePy is the Python library companion to the Acycle desktop application for cyclostratigraphy, time series analysis, and paleoclimate research. It provides a programmatic API for scripting and batch processing, and bundles the full Acycle desktop GUI tools.

Reference: Li, M., Hinnov, L., & Kump, L. (2019). Acycle: Time-series analysis software for paleoclimate research and education. Computers & Geosciences, 127, 12-22.

Installation

pip install acycle

Python >= 3.8. Dependencies (numpy, scipy, pandas, matplotlib, PySide6, scikit-image, scikit-learn, Pillow, qt-material, astropy) are installed automatically.

pip install 'acycle[dev]'      # pytest, black, flake8
pip install 'acycle[full]'     # sounddevice, h5py, netCDF4
pip install statsmodels        # faster lowess/loess detrending

Verify installation:

import acycle as ac
s = ac.Series(x=[1,2,3], y=[10,20,30])
print(s)  # Series(n=3, x=[1, 3] , y_name='y')

Quick Start: Reading and Processing Data

Read time series from whitespace, comma, or tab-delimited files. Auto-detection kicks in when delimiter is omitted:

import acycle as ac

s = ac.read_series('data.txt')                      # whitespace
s = ac.read_series('data.csv', delimiter=',')       # comma
s = ac.read_series('data.tsv', delimiter='\t')     # tab
lr04 = ac.load_example('lr04')                      # built-in dataset

The read_series function delegates to Series.from_file, which auto-detects the delimiter by scanning the first 5 lines, then loads into a Series object that carries its own metadata:

print(s.dt)           # median sampling interval
print(s.n)            # number of points
print(s.history)      # chain of operations applied

Operations return new Series, enabling a fluent chaining API:

s2 = s.clean(sort=True).interpolate(step=1.0).detrend(window=400).standardize()

Each method records itself in s2.history, so the full processing pipeline is traceable.

Core Concept: Series and Chainable Operations

A Series wraps two numpy arrays (x and y) with metadata (names, units, history). Under the hood:

class Series:
    def __init__(self, x, y, *, x_name='x', x_unit='', y_name='y', y_unit='', metadata=None):
        self.x = np.asarray(x, dtype=float).ravel()
        self.y = np.asarray(y, dtype=float).ravel()
        self.x_name, self.y_name = x_name, y_name
        self.x_unit, self.y_unit = x_unit, y_unit
        self.history = []

    def _record(self, method, **params):
        self.history.append({'method': method, 'params': params})

Every transform (interpolate, detrend, standardize, etc.) calls _record to append to the history chain, creating a fully auditable processing provenance.

Built-in Datasets

Result Objects

All analysis functions return structured result objects with a consistent interface. Each supports three methods:

result.to_dataframe()      # export to pandas DataFrame
result.save('prefix')      # write CSV to disk
result.plot()              # generate matplotlib figure
result.settings            # dict of computation parameters

The result classes are:

PSD.save() writes prefix_spectrum.csv (a custom override). All other classes use the base _ResultBase.save() which writes prefix.csv.

Example 1: Spectral Analysis on the LR04 Benthic Stack

This example loads the Lisiecki & Raymo (2005) benthic d18O record, preprocesses it, computes a multitaper power spectrum, and identifies Milankovitch periodicities.

Load the LR04 stack. It contains 2115 unevenly-spaced data points from 0 to 5320 ka:

import acycle as ac
import numpy as np

lr04 = ac.load_example('lr04')
print(f'n={lr04.n}, span={lr04.x_min:.0f}-{lr04.x_max:.0f} ka')
# n=2115, span=0-5320 ka

Interpolate to a uniform 1-kyr grid. The original data has a median spacing of ~2.5 kyr; linear interpolation gives 5321 evenly-spaced points:

s = lr04.interpolate(step=1.0, method='linear')
print(f'interpolated: {s.n} pts, dt={s.dt:.2f} ka')
# interpolated: 5321 pts, dt=1.00 ka

Detrend with a 400-kyr lowess window to isolate orbital-scale variability from long-term climate trends:

s_dt = s.detrend(window=400, window_unit='value', method='lowess')
print(f'detrended: mean={s_dt.y.mean():.4f}, std={s_dt.y.std():.4f}')

Compute the multitaper (MTM) power spectrum. nw=3 gives 2*nw-1=5 Slepian tapers, which balances frequency resolution against variance reduction. Each taper produces an independent spectral estimate; averaging across tapers suppresses noise while preserving signal:

from acycle import spectral, PSD

freq, power = spectral._mtm_spectrum(s_dt.y, dt=s_dt.dt, nw=3)
period = 1.0 / freq  # convert frequency to period (ka)

# Identify the top 10 peaks (skip zero-frequency)
peaks = np.argsort(power)[-10:][::-1]
for i in peaks:
    if period[i] > 1:
        print(f'  {period[i]:7.1f} kyr  power={power[i]:.4f}')

Wrap the result in a PSD object for saving and plotting:

psd = PSD(frequency=freq, power=power, period=period, method='MTM(nw=3)')
psd.save('lr04_spectrum')   # writes lr04_spectrum_spectrum.csv
psd.plot()                  # displays the power spectrum

For unevenly-spaced data without interpolation, use the Lomb-Scargle periodogram instead of MTM:

freq, power = spectral._lomb_scargle_spectrum(lr04.x, lr04.y, fmin=0.0002, fmax=0.5, pad=4000)

This calls astropy.timeseries.LombScargle if available (faster, more features), falling back to scipy.signal.lombscargle.

Example 2: Continuous Wavelet Transform

Wavelet analysis reveals how spectral power evolves through time, useful for detecting non-stationary behavior such as changes in dominant orbital forcing or sedimentation rate variations.

Load and prepare the 0-2000 ka section of LR04:

lr04 = ac.load_example('lr04')
s = lr04.select(0, 2000).interpolate(step=1.0)
s_dt = s.detrend(window=200, window_unit='value', method='lowess')

Run the continuous wavelet transform using the Morlet wavelet (a complex sinusoid modulated by a Gaussian envelope). The CWT convolves the signal with scaled and translated versions of the mother wavelet. Parameters:

• period_min=2, period_max=500: scan periods from 2 to 500 kyr
• dj=0.05: scale resolution (~100 steps per octave)
• pad=True: zero-pad the signal to the next power of two for faster FFT computation

from acycle import wavelet, WaveletResult

period, power, coi, significance = wavelet.cwt(
    s_dt.y, dt=s_dt.dt,
    period_min=2, period_max=500,
    dj=0.05, mother='MORLET', pad=True
)

The cone of influence (COI) marks the region where edge effects become significant. Data outside the COI should be interpreted with caution.

Wrap the results and export:

wr = WaveletResult(x=s_dt.x, period=period, power=power,
                   coi=coi, significance=significance, mother='MORLET')
wr.to_dataframe()               # periods as rows, time as columns
wr.plot(cmap='viridis')         # heatmap with COI overlay
wr.save('lr04_wavelet')         # writes lr04_wavelet.csv

Extract the time-averaged (global) wavelet spectrum:

global_power = power.mean(axis=1)
peaks = np.argsort(global_power)[-5:][::-1]
for i in peaks:
    print(f'  {period[i]:.1f} kyr: power={global_power[i]:.4f}')

Cross-wavelet coherence compares two series in time-frequency space:

s1, s2 = s.select(0, 1000), s.select(500, 1500)
coh = wavelet.wavelet_coherence(s1.y[:500], s2.y[-500:], dt=s1.dt)
# coh has keys: 'period', 'coherence', 'phase', 'coi', 'dt'

Example 3: Age Modeling and Depth-to-Time Tuning

This example builds an age model by counting 405-kyr eccentricity cycles, then tunes a depth-domain proxy to the time domain. The key concept: orbital cycles in sedimentary records provide a natural chronometer. By identifying cycle boundaries, we pin absolute ages to depth positions.

Generate a synthetic depth-domain series with a 405-kyr eccentricity signal. At an average sedimentation rate of ~2 cm/kyr, one 405-kyr cycle occupies ~8.1 meters of stratigraphy:

np.random.seed(42)
depth = np.arange(0, 200, 0.1)  # 200 meters at 10 cm spacing
signal = np.sin(2 * np.pi * depth / 8.1)  # 405-kyr cycle at 2 cm/kyr
y = signal + 0.3 * np.random.randn(len(depth))
s = ac.Series(x=depth, y=y, x_name='Depth', x_unit='m')
s = s.interpolate(step=0.1)

build_age_model identifies cycle peaks (or troughs) and assigns ages. Each detected peak becomes a tie point, separated by exactly 405 kyr:

from acycle import age, AgeModel

model = age.build_age_model(s.x, s.y, cycle_period=405,
    anchor='max', start_age=0, age_direction='increasing')
print(f"tie points: {len(model['tie_points'])}, sed rate: {model['sed_rate'].mean():.3f} m/kyr")
# tie points: ~24, sed rate: ~0.020 m/kyr (= 2 cm/kyr)

The sedimentation rate between consecutive tie points reveals changes in accumulation rate through time.

Tune the depth-domain series to the time domain using scipy.interpolate.interp1d:

tuned_x, tuned_y = age.tune(s.x, s.y, model)
tuned = ac.Series(x=tuned_x, y=tuned_y, x_name='Age', x_unit='ka')

Verify the tuning worked by checking that the 405-kyr peak dominates the tuned spectrum:

tuned_dt = tuned.detrend(window=400, method='lowess')
freq, pwr = ac.spectral._mtm_spectrum(tuned_dt.y, dt=tuned_dt.dt, nw=2)
peak = (1.0 / freq)[np.argmax(pwr[1:]) + 1]
print(f'dominant period: {peak:.1f} kyr (expected ~405)')

Create and save the AgeModel result:

am = AgeModel(depth=model['depth'], age=model['age'],
              sed_rate=model['sed_rate'])
am.to_dataframe()
am.plot()    # Harker diagram: depth vs age
am.save('my_age_model')

Example 4: Preprocessing Pipeline

A realistic preprocessing workflow: read a noisy record, clip extreme values, remove a known bad interval, detrend, and standardize.

import acycle as ac
import numpy as np
from acycle import preprocess

x = np.arange(0, 100, 0.2)
y = np.sin(2 * np.pi * x / 20) + 0.1 * x + 0.3 * np.random.randn(len(x))
s = ac.Series(x=x, y=y, x_name='Depth', x_unit='m', y_name='GR')

clip_by_threshold removes rows where y exceeds the threshold. The mode parameter controls whether data is deleted, capped, or zeroed:

s_clipped = s.clip_by_threshold(threshold=2.5, side='above', mode='delete')

remove_sections deletes specified depth ranges from the record:

x_clean, y_clean = preprocess.remove_sections(s_clipped.x, s_clipped.y, [(20, 25)])

Detrend with a first-order polynomial to remove the linear trend from the synthetic data:

y_dt, trend = preprocess.detrend(x_clean, y_clean, method='polynomial', poly_order=1)

Standardize to zero mean and unit variance:

s_final = ac.Series(x=x_clean, y=y_dt)
s_final = s_final.standardize()
print(f'n={s_final.n}, mean={s_final.y.mean():.3f}, std={s_final.y.std():.3f}')
# n=473, mean=0.000, std=1.000

Example 5: Merging and Multiplying Series

When comparing records from different sources, resample to a common time grid and merge into a single DataFrame:

lr04 = ac.load_example('lr04').interpolate(step=1.0)
ceno = ac.load_example('cenogrid_d18o').select(0, 5320).interpolate(step=1.0)
df = preprocess.merge_series([
    (lr04.x,  lr04.y,  'LR04'),
    (ceno.x,  ceno.y,  'CENOGRID'),
], require_same_x=True)
print(f'merged: {df.shape}')  # (5321, 3)

multiply_series computes element-wise products. When require_same_x=True, both series must share the same x grid:

x, y = preprocess.multiply_series(lr04.x, lr04.y, ceno.y, require_same_x=True)

When require_same_x=False, the second series is linearly interpolated onto the first series' grid before multiplication.

Example 6: Changepoint Detection and PCA

Detect structural breaks in a time series using the PELT algorithm (via the ruptures package) or Bayesian changepoint analysis:

y = np.random.randn(500)
y[200:350] += 2.0  # insert a mean shift
result = preprocess.changepoint(y, method='mean_shift', min_size=10)
print(f"changepoints: {result['changepoints']}")

Principal component analysis on multi-column data:

import pandas as pd
data = pd.DataFrame({'x': x, 'y1': y1, 'y2': y2, 'y3': y3})
pca_result = preprocess.pca(data, x_col=0, n_components=2)
# pca_result['pcs']         -- principal component scores
# pca_result['variance_explained'] -- explained variance ratios
# pca_result['loadings']    -- component loadings

API Reference

Series

spectral

wavelet

filter

age

preprocess

io

GUI Tools

Notes

pandas >= 2.0: The delim_whitespace parameter (removed in pandas 2.0) has been replaced with sep=r"\s+" throughout the codebase.

Headless servers: The programmatic API works without any display server. GUI tools require a graphical environment (X11/Wayland/Windows GUI).

statsmodels: Installing statsmodels improves detrend(method='lowess') performance. A pure-SciPy fallback is included.

Encoding on Windows: Chinese Windows uses cp936/GBK as the default encoding. For UTF-8 files, pass encoding='utf-8' to read_series().

Large files: For datasets >100 MB, load with pandas directly:

import pandas as pd
df = pd.read_csv('large.csv', usecols=[0, 1])
s = ac.Series(x=df.iloc[:, 0].values, y=df.iloc[:, 1].values)

Citation

Li, M., Hinnov, L., & Kump, L. (2019). Acycle: Time-series analysis software for paleoclimate research and education. Computers & Geosciences, 127, 12-22.

License

MIT

Links

• PyPI
• Repository
• Documentation