Metadata-Version: 2.4
Name: acycle
Version: 0.7.1
Summary: AcyclePy: Advanced cyclostratigraphy, time series analysis, and paleoclimate toolkit
Author-email: Meng Wang <grampus90@gmail.com>
License-Expression: MIT
Project-URL: Homepage, https://github.com/mingsongli/AcyclePy
Project-URL: Repository, https://github.com/mingsongli/AcyclePy.git
Project-URL: Issue Tracker, https://github.com/mingsongli/AcyclePy.git/issues
Project-URL: Documentation, https://github.com/mingsongli/AcyclePy/tree/main/doc
Keywords: cyclostratigraphy,time-series,spectral,wavelet,geology
Classifier: Programming Language :: Python :: 3
Requires-Python: >=3.8
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.20.0
Requires-Dist: pandas>=1.3.0
Requires-Dist: PySide6>=6.0.0
Requires-Dist: matplotlib>=3.3.0
Requires-Dist: scipy>=1.6.0
Requires-Dist: scikit-image>=0.18.0
Requires-Dist: scikit-learn>=0.24.0
Requires-Dist: Pillow>=8.0.0
Requires-Dist: qt-material>=2.0.0
Requires-Dist: astropy>=5.0
Provides-Extra: dev
Requires-Dist: pytest>=7.0.0; extra == "dev"
Requires-Dist: black>=22.0.0; extra == "dev"
Requires-Dist: flake8>=4.0.0; extra == "dev"
Provides-Extra: full
Requires-Dist: sounddevice>=0.4.0; extra == "full"
Requires-Dist: h5py>=3.0.0; extra == "full"
Requires-Dist: netCDF4>=1.5.0; extra == "full"
Dynamic: license-file

# AcyclePy - Cyclostratigraphy and Time Series Analysis Toolkit

[![PyPI](https://img.shields.io/pypi/v/acycle)](https://pypi.org/project/acycle/)
[![Python](https://img.shields.io/pypi/pyversions/acycle)](https://pypi.org/project/acycle/)
[![License](https://img.shields.io/pypi/l/acycle)](LICENSE)

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

```bash
pip install acycle
```

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

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

Verify installation:

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

---

## Quick Start: Reading and Processing Data

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

```python
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:

```python
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:

```python
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:

```python
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

| Name | Description | Points | Reference |
|------|-------------|--------|-----------|
| `lr04` | LR04 benthic d18O | 2,115 | Lisiecki & Raymo (2005) |
| `cenogrid_d18o` | CENOGRID d18O | 23,659 | Westerhold et al. (2020) |
| `cenogrid_d13c` | CENOGRID d13C | 23,666 | Westerhold et al. (2020) |
| `la2004_etp` | La2004 ETP solution | 2,001 | Laskar et al. (2004) |
| `petm_logfe` | Svalbard PETM logFe | 392 | Charles et al. (2011) |
| `wayao_gr` | Wayao Carnian GR | 499 | Li et al. (2018) |
| `guandao_gr` | Guandao Anisian GR | 1,071 | Li et al. (2018) |
| `rednoise_0.7_2000` | Synthetic red noise | 2,000 | -- |
| `launa_loa_co2` | Mauna Loa CO2 | 722 | NOAA |
| `newark_depth_rank` | Newark Depth Rank | -- | Olsen & Kent (1999) |
| `csa_extinction` | CSA extinction | 10 | -- |

---

## Result Objects

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

```python
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:

| Class | Description | Key Attributes |
|-------|-------------|----------------|
| `PSD` | Power spectral density | `frequency`, `power`, `period` |
| `EvolutiveSpectrum` | Evolutionary (time-varying) spectrum | `x`, `frequency`, `power` |
| `WaveletResult` | Continuous wavelet transform | `x`, `period`, `power`, `coi` |
| `FilterResult` | Filtered signal + amplitude envelope | `filtered`, `amplitude`, `phase` |
| `AgeModel` | Age-depth model + sedimentation rate | `depth`, `age`, `sed_rate` |
| `CocoResult` | COCO/eCOCO correlation | `sed_rate`, `rho`, `p_value` |
| `SedNoiseResult` | DYNOT sedimentation noise | `age`, `median`, `quantiles` |

`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:

```python
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:

```python
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:

```python
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:

```python
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:

```python
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:

```python
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:

```python
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

```python
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:

```python
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:

```python
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:

```python
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:

```python
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:

```python
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`:

```python
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:

```python
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:

```python
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.

```python
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:

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

`remove_sections` deletes specified depth ranges from the record:

```python
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:

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

Standardize to zero mean and unit variance:

```python
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:

```python
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:

```python
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:

```python
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:

```python
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

| Method | Description |
|--------|-------------|
| `Series(x, y)` | Construct from arrays |
| `.from_file(path)` | Read from delimited text |
| `.clean(sort, duplicate, dropna)` | Sort, deduplicate, drop NaN |
| `.interpolate(step, method)` | Uniform grid interpolation |
| `.interpolate_pro(step, method)` | Advanced interpolation with gap filling |
| `.interpolate_to(ref)` | Interpolate onto reference grid |
| `.detrend(window, method)` | Remove trend (lowess/loess/polynomial/savgol) |
| `.standardize()` | Z-score standardization |
| `.log10()` | Base-10 logarithm |
| `.derivative(order)` | Numerical derivative |
| `.prewhiten(method)` | Prewhitening |
| `.select(start, stop)` | Sub-range selection |
| `.moving_average(n)` | Moving average |
| `.gaussian_smooth(n, sigma)` | Gaussian smoothing |
| `.moving_median(n)` | Moving median |
| `.multiply(other)` | Element-wise multiply with another Series |
| `.clip_by_threshold(t, side, mode)` | Clip data by threshold |
| `.to_dataframe()` | Export to pandas |
| `.copy()` | Deep copy |
| `.dt` | Median sampling interval (read-only) |
| `.n` | Number of points (read-only) |

### spectral

| Function | Description |
|----------|-------------|
| `_periodogram(y, dt)` | Classical periodogram |
| `_mtm_spectrum(y, dt, nw)` | Multitaper method (Slepian tapers) |
| `_lomb_scargle_spectrum(x, y)` | Lomb-Scargle for uneven spacing |
| `_estimate_ar1_rho(y)` | AR(1) lag-1 autocorrelation |
| `_estimate_ar1_noise(y, dt)` | AR(1) noise background |
| `_ftest_mtm(y, dt)` | F-test for peak significance |

### wavelet

| Function | Description |
|----------|-------------|
| `cwt(y, dt, period_min, period_max, dj, mother)` | Continuous wavelet transform (Torrence & Compo 1998) |
| `wavelet_coherence(y1, y2, dt)` | Wavelet coherence and cross phase |

### filter

| Function | Description |
|----------|-------------|
| `apply_filter(y, dt, kind, method)` | Bandpass/lowpass/highpass (Gaussian, Butterworth, etc.) |
| `dynamic_filter(x, y, window)` | Sliding-window dynamic filter |
| `amplitude_modulation(x, y, flow, fhigh)` | Amplitude envelope via Hilbert transform |

### age

| Function | Description |
|----------|-------------|
| `build_age_model(x, y, cycle_period)` | Age model from cycle counting |
| `sedrate_to_age_model(depth, sedrate)` | Sedimentation rate to age model |
| `tune(depth, y, age_model)` | Depth-to-time conversion |
| `stratigraphic_correlation(ref, target, ties)` | Correlate two stratigraphic sections |

### preprocess

| Function | Description |
|----------|-------------|
| `detrend(x, y, window, method)` | Remove trend from (x, y) arrays |
| `interpolate_pro(x, y, step)` | Advanced interpolation |
| `clip_by_threshold(x, y, t)` | Clip by threshold value |
| `remove_sections(x, y, sections)` | Remove data ranges |
| `add_gaps(x, y, gaps)` | Insert NaN-filled gaps |
| `remove_peaks(y, ymin, ymax)` | Remove or cap peak values |
| `multiply_series(x, y1, y2)` | Element-wise product of two series |
| `merge_series(series_list)` | Merge multiple series by common x |
| `pca(data, n_components)` | Principal component analysis |
| `changepoint(y, method)` | Changepoint detection (PELT/Bayesian) |
| `transform_xy(x, y, a,b,c,d)` | Affine coordinate transform |
| `find_extreme(x, y)` | Find max or min in range |

### io

| Function | Description |
|----------|-------------|
| `read_series(path)` | Read Series from file (auto-detect delimiter) |
| `write_series(series, path)` | Write Series to file |
| `load_example(name)` | Load built-in dataset |
| `load_lr04(start, stop)` | Load LR04 benthic stack |
| `load_cenogrid(variable)` | Load CENOGRID isotope data |
| `extract_columns(path, x_col, y_col)` | Extract columns from multi-column file |

---

## GUI Tools

| Command | Tool |
|---------|------|
| `acycle-imageprocessor` | Image digitizing: load image, calibrate coordinates, extract data points by color matching |
| `acycle-plot` | PlotPro: publication-quality plotting with multi-panel and subplot support |
| `acycle-interpolation` | Advanced interpolation with multiple methods and gap filling |
| `acycle-data-extractor` | Extract data segments by specifying start/end ranges |
| `acycle-section-remover` | Remove data sections with optional time adjustment |
| `acycle-gap-adder` | Insert NaN-filled gaps at specified positions |
| `acycle-data-clipper` | Clip data above/below thresholds |
| `acycle-image-analyzer` | Advanced image processing and analysis |

---

## 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:

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

---

## 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](https://pypi.org/project/acycle/)
- [Repository](https://github.com/mingsongli/AcyclePy)
- [Documentation](https://github.com/mingsongli/AcyclePy/tree/main/doc)
