AcyclePy -- Cyclostratigraphy and Time Series Analysis Toolkit

AcyclePy is the Python library companion to the Acycle desktop application, providing a programmatic API for cyclostratigraphy, time series analysis, and paleoclimate research. It also 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

ash pip install acycle

Requirements: Python >= 3.8. Core dependencies are installed automatically.

Optional dependencies:

ash pip install "acycle[dev]" # pytest, black, flake8 pip install "acycle[full]" # sounddevice, h5py, netCDF4 pip install statsmodels # better lowess performance

Verify:

python 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')

Result Objects

All analysis functions return structured result objects with o_dataframe(), save(), and plot() methods:

python psd = PSD(frequency=freq, power=power) psd.to_dataframe() # pandas DataFrame psd.save("my_spectrum") # writes my_spectrum_spectrum.csv psd.plot() # matplotlib figure psd.settings # parameters used

Example 1: Spectral Analysis of the LR04 Benthic Stack

This example walks through loading a classic paleoclimate dataset, preprocessing it, computing its power spectrum, and identifying Milankovitch cycles.

`python import acycle as ac import numpy as np

Step 1: Load the LR04 benthic d18O stack (Lisiecki & Raymo, 2005).

This is a 0-5320 ka record with 2115 unevenly-spaced data points.

lr04 = ac.load_example("lr04") print(f"Loaded: {lr04.n} points, span {lr04.x_min:.0f}-{lr04.x_max:.0f} ka")

-> Loaded: 2115 points, span 0-5320 ka

Step 2: Interpolate to a uniform 1-kyr grid.

The original data has variable spacing (median ~2.5 kyr).

Linear interpolation onto a 1-kyr grid gives 5321 evenly-spaced points.

lr04_uniform = lr04.interpolate(step=1.0, method="linear") print(f"Interpolated: {lr04_uniform.n} points, dt={lr04_uniform.dt:.2f} ka")

-> Interpolated: 5321 points, dt=1.00 ka

Step 3: Detrend with a 400-kyr lowess window.

This removes long-term trends (>400 kyr) while preserving orbital cycles.

lr04_dt = lr04_uniform.detrend(window=400, window_unit="value", method="lowess") print(f"Detrended: mean={lr04_dt.y.mean():.4f}, std={lr04_dt.y.std():.4f}")

Step 4: Compute the multitaper power spectrum.

nw=3 gives 5 tapers (2*nw-1), balancing frequency resolution and variance.

from acycle import spectral, PSD freq, power = spectral._mtm_spectrum(lr04_dt.y, dt=lr04_dt.dt, nw=3) period = 1.0 / freq # convert frequency to period (kyr)

Step 5: Identify the top 10 spectral peaks.

peak_idx = np.argsort(power)[-10:][::-1] print("\nTop 10 spectral peaks:") for i in peak_idx: if period[i] > 1: # skip the zero-frequency component print(f" Period: {period[i]:7.1f} kyr Power: {power[i]:.4f}")

Typical output shows peaks near 405, 100, 41, 23, 19 kyr.

Step 6: Save the spectrum.

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

Example 2: Continuous Wavelet Transform

This example shows how to perform wavelet analysis to see how spectral power evolves through time, useful for detecting non-stationary behavior in climate records.

`python import acycle as ac import numpy as np

Step 1: Load and prepare the LR04 0-2000 ka section.

lr04 = ac.load_example("lr04") s = lr04.select(0, 2000).interpolate(step=1.0)

Step 2: Detrend to remove the long-term d18O trend.

s_dt = s.detrend(window=200, window_unit="value", method="lowess")

Step 3: Compute the continuous wavelet transform (CWT).

Morlet wavelet (default), scanning periods from 2-500 kyr.

dj=0.05 gives fine period resolution (~100 scale steps per octave).

pad=True zero-pads to the next power of 2 for faster FFT.

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 )

Step 4: Wrap in a WaveletResult for easy access.

wr = WaveletResult( x=s_dt.x, period=period, power=power, coi=coi, significance=significance, mother="MORLET" )

Step 5: Export to DataFrame (periods as rows, time as columns).

df = wr.to_dataframe() print(f"Wavelet power shape: {df.shape}")

Step 6: Extract the time-averaged (global) wavelet spectrum.

global_power = power.mean(axis=1)

Step 7: Find dominant periods.

peak_idx = np.argsort(global_power)[-5:][::-1] print("\nDominant periods (global wavelet):") for i in peak_idx: print(f" {period[i]:.1f} kyr: power={global_power[i]:.4f}")

Step 8: Plot and save.

wr.plot(cmap="viridis") # time-period heatmap with cone of influence wr.save("lr04_wavelet") # saves to lr04_wavelet.csv

Step 9: Cross-wavelet coherence between two series segments.

s1 = s.select(0, 1000) s2 = s.select(500, 1500) coh = wavelet.wavelet_coherence(s1.y[:500], s2.y[-500:], dt=s1.dt) print(f"Coherence keys: {list(coh.keys())}")

-> ['period', 'coherence', 'phase', 'coi', 'dt']

`

Example 3: Full Age Modeling and Tuning Workflow

This example demonstrates building an age model from cycle counting, tuning a depth-domain series to the time domain, and computing sedimentation rates.

`python import acycle as ac import numpy as np

Step 1: Generate a synthetic depth-domain series with 405-kyr cycles.

np.random.seed(42) depth = np.arange(0, 200, 0.1) # 0-200 meters at 10-cm spacing

405 kyr * ~2 cm/kyr = ~8.1 m per cycle -> ~24.7 cycles over 200 m

signal = np.sin(2 * np.pi * depth / 8.1) noise = 0.3 * np.random.randn(len(depth)) y = signal + noise s = ac.Series(x=depth, y=y, x_name="Depth", x_unit="m", y_name="Proxy")

Step 2: Interpolate to uniform spacing.

s = s.interpolate(step=0.1, method="linear")

Step 3: Build an age model by counting 405-kyr cycles.

anchor="max" counts from eccentricity maxima.

start_age=0 sets the first counted cycle to 0 ka.

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" )

Step 4: Examine the age model.

print(f"Tie points: {len(model['tie_points'])}") print(f"Depth range: {model['depth'][0]:.1f}-{model['depth'][-1]:.1f} m") print(f"Age range: {model['age'][0]:.1f}-{model['age'][-1]:.1f} ka")

Expected: ~24 tie points covering ~0-9720 ka

Step 5: Compute sedimentation rate (m/kyr).

sed_rates = model["sed_rate"] print(f"Mean sedimentation rate: {sed_rates.mean():.3f} m/kyr") print(f" = {sed_rates.mean()*100:.1f} cm/kyr")

Step 6: Tune the depth-domain series to time domain.

tuned_x, tuned_y = age.tune(s.x, s.y, model) tuned_series = ac.Series(x=tuned_x, y=tuned_y, x_name="Age", x_unit="ka") print(f"Tuned series: {tuned_series.n} pts, dt={tuned_series.dt:.2f} ka")

Step 7: Verify tuning by checking the 405 kyr peak.

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

Step 8: Create and save an AgeModel result object.

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") `

Quick Reference

Reading data: python s = ac.read_series("file.txt") # whitespace s = ac.read_series("file.csv", delimiter=",") # comma lr04 = ac.load_example("lr04") # built-in dataset

Series processing: python s.clean(sort=True).interpolate(step=1.0).detrend(window=400).standardize()

Spectral analysis: python freq, power = ac.spectral._mtm_spectrum(y, dt=1.0, nw=3) freq, power = ac.spectral._lomb_scargle_spectrum(x, y)

Wavelet analysis: python period, power, coi, sig = ac.wavelet.cwt(y, dt=1.0, period_min=2, period_max=500) coh = ac.wavelet.wavelet_coherence(y1, y2, dt=1.0)

Filtering: python result = ac.filter.apply_filter(y, dt=1.0, kind="bandpass", flow=0.01, fhigh=0.05) result = ac.filter.amplitude_modulation(x, y, flow=0.01, fhigh=0.05)

Age modeling: python model = ac.age.build_age_model(depth, y, cycle_period=405) tuned_x, tuned_y = ac.age.tune(depth, y, model)

Preprocessing: python yd, trend = ac.preprocess.detrend(x, y, window=0.35, method="lowess") x, y = ac.preprocess.multiply_series(x, y1, y2) df = ac.preprocess.merge_series([(x1,y1,"A"), (x2,y2,"B")])

GUI Tools

ash acycle-imageprocessor acycle-interpolation acycle-gap-adder acycle-plot acycle-data-extractor acycle-data-clipper acycle-image-analyzer acycle-section-remover

Notes

pandas >= 2.0: Tested with pandas 2.x. The removed delim_whitespace parameter has been replaced with sep=r"\s+".

Headless servers: The programmatic API works without any display. GUI tools require a graphical environment.

statsmodels: Install statsmodels for faster lowess/loess detrending. Falls back to pure-SciPy if unavailable.

File encoding on Windows: Chinese Windows uses cp936/GBK by default. For UTF-8 files, pass encoding="utf-8" to ead_series().

Large files: For files > 100 MB, load with pandas and construct Series: python 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)

Built-in Datasets

| ewark_depth_rank | Newark Depth Rank | -- | Olsen & Kent (1999) | | ednoise_0.7_2000 | Synthetic red noise | 2,000 | -- | | launa_loa_co2 | Mauna Loa CO2 | 722 | NOAA | | csa_extinction | CSA extinction | 10 | -- |

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